Coupling of Finite Element Method and Peridynamics to Simulate Ship-Ice Interaction

Abstract

In this work, the finite element method (PD-FEM) coupling strategy is used to simulate ship-ice interaction. Two numerical benchmark tests are selected to validate the coupling approach and its program. During the ice-breaking process simulation, the generation and propagation of radial and circular cracks in level ice are modeled and phenomena such as the shedding of wedge ice, flipping of brash ice, and cleaning of the channel are observed to be broadly consistent with experimental observation. The influence of ship speed and ice thickness on the ice load are investigated and analyzed. The ice load obtained from the numerical simulations is in general agreement with that given by Lindqvist’s empirical formula. The boundary effect on the crack path can also be avoid with the current coupling method.

1. Introduction

In recent years, global climate change and the melting of ice in Arctic regions has raised the possibility of exploiting Arctic resources and opening an Arctic channel [1,2]. The exploitation of resources and scientific research in Arctic regions rely on icebreakers to open the necessary routes [3,4,5]. Therefore, it is great significant to simulate the icebreaking scenarios and calculate the ice load of ship–ice interaction, and it helps in improving the design and safe navigation of icebreakers. The ship-ice interaction scenarios are studied with full-scale tests, model tests, theoretical analyses, and numerical simulations. For full-scale testing, the results are reliable, but the associated cost is high. Model test is a promising candidate to study the ship–ice interaction [6,7,8,9]. However, compared with full-scale tests, models have many uncertainties, and can be expensive and time-consuming [10,11]. Theoretical analysis is still challenging in some cases, such as dealing with complicated structures [12]. Fortunately, numerical methods to study ship-ice interactions do not need to consider the structure complexity, and are not restricted by factors such as geography, cost, and time, and have been shown to be both efficient and accurate, both in theoretical research and engineering application [13,14,15,16]. Finite element method (FEM) was successfully applied to estimate the strength of ship structure problems [10,17,18]. The discrete element method (DEM) to calculate ice loads for offshore structures and ships [19,20,21,22]. Smoothed-particle hydrodynamics method (SPH) was adopted in the ice field to simulate the ice-structure interaction dynamics [23,24], and other methods [25,26].

In recent years, a mesh-free method of peridynamics was proposed [27]. This reformulation of the classical continuum mechanics is a non-local theory that does not assume the spatial differentiability of displacement fields. Based on integrodifferential equations, peridynamics can deal with discontinuous displacement fields. Therefore, it can simulate spontaneous crack nucleation and propagation, and can be used to simulate the ice-breaking and calculate ice loads [28,29,30,31,32,33,34,35,36,37]. However, as a non-local theory computational efficiency of peridynamics is far lower than that of FEM, especially for engineering applications like ship-ice interaction. To improve its computational efficiency, researchers have coupled peridynamics with FEM. Macek and Silling [38] proposed the PD-FEM coupling approach and implemented peridynamics in a commercial finite element analysis code, Liu et al. [39] introduced interface elements to calculate the coupling force in a PD-FEM approach, and Lee et al. [40] proposed a coupled PD-FEM approach to analyze impact fractures. To date, the advantages of combining PD with FEM have been demonstrated in applications to concrete and composite materials, but PD-FEM coupling has not been used to deal with the ship–ice interaction. In this work, the coupling strategy proposed by Liu et al. is employed for its easy to implement and robust theory foundation.

The following work is organized as, peridynamics theory and PD-FEM coupling scheme is introduced in Section 2 and Section 3, respectively. The proposed coupling approach is verified with both dynamic and static cases in Section 4. The ship-ice interaction is simulated in Section 5. Conclusion is drawn in Section 6.

2. Peridynamics Framework

Peridynamics assumes that the continuum body is composed of small particles. Each particle interacts with other particles within a finite distance $\delta$ called the horizon. The pairwise interaction between two particles exists despite they are not in contact. This physical interaction is referred to as a bond, which in some way has a close analogy to a mechanical spring. In bond-based peridynamics, the kinetic equation of particle x in the reference configuration at time t is

$$\rho \ddot{\mathbf{u}}\left(\mathbf{x},t\right)={\displaystyle {\int }_{H_{\mathbf{x}}}\mathbf{f}\left(\mathbf{u}\left({\mathbf{x}}^{\prime },t\right)-\mathbf{u}\left(\mathbf{x},t\right),{\mathbf{x}}^{\prime }-\mathbf{x}\right)}d{V}_{{\mathbf{x}}^{\prime }}+\mathbf{b}\left(\mathbf{x},t\right)‖{\mathbf{x}}^{\prime }-\mathbf{x}‖\le \delta$$

(1)

where $H_{\mathbf{x}}$ is the domain of integration within the horizon of particle $\mathbf{x}$, u is the displacement vector field, and b is the body force density. $\rho$ is the mass density, and $\mathbf{f}$ is a pairwise force density function defined as the force per unit volume that particle ${\mathbf{x}}^{\prime }$ exerts on particle x, which contains all the constitutive information of the materials.

To simplify the notation, the relative position in the initial configuration and its relative displacement are denoted as $\mathsf{\xi }={\mathbf{x}}^{\prime }-\mathbf{x}$ and $\mathsf{\eta }=\mathbf{u}({\mathbf{x}}^{\prime },t)-\mathbf{u}(\mathbf{x},t)$, respectively. Therefore, the relative position of the two interacting particles at t in the current configuration is $\mathsf{\xi }+\mathsf{\eta }$ and the pairwise force density function can be described as $\mathbf{f}\left(\mathsf{\eta },\mathsf{\xi }\right)$.

For the prototype micro-elastic brittle (PMB) material defined by Silling and Askari [41], the pairwise force density function can be expressed as

$$\mathbf{f}\left(\mathsf{\eta },\mathsf{\xi }\right)=\frac{\mathsf{\xi }+\mathsf{\eta }}{\left|\mathsf{\xi }+\mathsf{\eta }\right|}cs\mu \left(t,\mathsf{\eta },\mathsf{\xi }\right)\forall \mathsf{\eta },\mathsf{\xi }$$

(2)

where $c=12E/\pi {\delta }^4$ is the micro modulus, E is Young’s modulus, and s(η,ξ) is denoted as the stretch of the bond, which can be defined as

$$s=\frac{\left|\mathsf{\xi }+\mathsf{\eta }\right|-\left|\mathsf{\xi }\right|}{\left|\mathsf{\xi }\right|}$$

(3)

When the deformation stretch s exceeds a limit s₀ (described as the critical stretch for failure), the bond between the two particles breaks and no pairwise force remains. The term $\mu \left(t,\mathsf{\eta },\mathsf{\xi }\right)$ is a history-dependent scalar-valued function, which is introduced to represent the bond failure of two particles. This can be defined as

$$\mu \left(t,\mathsf{\eta },\mathsf{\xi }\right)=\left\{\begin{array}{l}1, s<s_0\\ 0, s\ge s_0\end{array}\right.$$

(4)

Accordingly, the level of damage is illustrated by the local damage at one particle, defined as

$$\phi \left(\mathbf{x},t\right)=1-\frac{{\displaystyle {\int }_{{{H}}_{\mathbf{x}}}\mu \left(t,\mathsf{\eta },\mathsf{\xi }\right)d{V}_{\mathsf{\xi }}}}{{\displaystyle {\int }_{{{H}}_{\mathbf{x}}}d{V}_{\mathsf{\xi }}}}$$

(5)

When solving the elastic problem in which the damage is not considered, the critical stretch can be set to infinity. Dealing with the damage problem, the value of s₀ can be obtained from the energy release rate.

3. Coupling of PD-FEM

3.1. Coupling Scheme

The PD-FEM coupling approach proposed by Liu et al. [39] is adopted and presented. The coupling scheme is as follows: the domain to be solved is partitioned into FEM subregions, which are modeled as a non-failure area and a PD subregion containing the area expected to be damaged. An interface element is introduced to bridge from the FEM subregion to the PD subregion. The interface element contains several peridynamic nodes for calculating the coupling forces, which are the interaction forces between embedded peridynamics nodes and peridynamics nodes outside the interface element. The coupling scheme is illustrated in Figure 1.

To implement the coupling scheme, interfaces between the peridynamics subregion and the FEM subregion should be defined prior to analysis. The coupling forces on embedded nodes are divided between the FEM nodes on the interface segment, as shown in Figure 2, by

$$f_{i,cp}=N_i\left({\xi }^{\prime },{\eta }^{\prime },{\psi }^{\prime }\right)f_p$$

(6)

where $f_{i,cp}$ is the force of the FEM nodes on the interface segment, $N_i$ is the shape function on the interface segment, $f_p$ is the coupling force on the embedded nodes, $\left({\xi }^{\prime },{\eta }^{\prime },{\psi }^{\prime }\right)$ are the natural coordinates of the projection of an embedded node onto the interface segment, i is the number of FEM nodes on the interface segment, and m is the total number of embedded nodes. Note that for FEM nodes that are not on the interface segment, $f_{i,cp}=0$.

For the FEM subregion, the equation of motion for the FEM nodes is written as

$$M_i{\ddot{u}}_i^n=F_{i,ext}+F_{i,int}$$

(7)

where $F_{i,e\mathbf{x}t}$ is the external force and the internal force is given by FEM nodes on segment FEM nodes not on segment,

$$F_{i,int}=f_{i,fem}+f_{i,cp}=\left\{\begin{array}{l}{\left[{\displaystyle \sum _e{K}^{\left(e\right)}{u}^{\left(e\right)}}\right]}_i+f_{i,cp},\hspace{0.17em}\mathrm{FEM}\hspace{0.17em}\mathrm{nodes}\hspace{0.17em}\mathrm{on}\hspace{0.17em}\mathrm{segment}\\ {\left[{\displaystyle \sum _e{K}^{\left(e\right)}{u}^{\left(e\right)}}\right]}_i ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{FEM}\hspace{0.17em}\mathrm{nodes}\hspace{0.17em}\mathrm{not}\hspace{0.17em}\mathrm{on}\hspace{0.17em}\mathrm{segment}\end{array}\right.$$

(8)

The displacements of the embedded peridynamics nodes are determined by

$$u_{ep}={\displaystyle \sum _{i=1}^8{N}_i\left(\xi ,\eta ,\psi \right)}{u}_i i=1,\dots ,8$$

(9)

where $\left(\xi ,\eta ,\psi \right)$ are the natural coordinates of an embedded peridynamics node in the interface element and $u_i$ is the nodal displacement of an interface element.

3.2. Numerical Implementation

The peridynamics equation of motion after discretization is written as

$$\rho {\ddot{\mathbf{u}}}_i^n={\displaystyle \sum _{j=1}^m\mathbf{f}\left({\mathbf{u}}_j^n-{\mathbf{u}}_i^n,{\mathbf{x}}_j-{\mathbf{x}}_i\right)V_j}+{\mathbf{b}}_i^n$$

(10)

For the FEM subregion, the equation of motion for the FEM nodes is written as

$$M_i{\ddot{u}}_i^n=F_{i,ext}^n+F_{i,\mathrm{int}}^n$$

(11)

where n denotes the number of time steps. The displacement of node i can be obtained by approximating the acceleration in Equations (12) and (13) using an explicit central difference formula

$${\ddot{\mathbf{u}}}_i^n=\frac{{\mathbf{u}}_i^{n+1}-2{\mathbf{u}}_i^n+{\mathbf{u}}_i^{n-1}}{\mathrm{\Delta }{t}^2}$$

(12)

$${\mathbf{u}}_i^{n+1}=\left\{\begin{array}{l}\frac{\mathrm{\Delta }{t}^2}{\rho }\left[{\displaystyle \sum _{j=1}^m\mathbf{f}\left({\mathbf{u}}_j^n-{\mathbf{u}}_i^n,{\mathbf{x}}_j-{\mathbf{x}}_i\right)V_j}+{\mathbf{b}}_i^n\right]+2{\mathbf{u}}_i^n-{\mathbf{u}}_i^{n-1}, \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{P}\mathrm{D} \mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{s}\\ \frac{\mathrm{\Delta }{t}^2}{\rho }\left[F_{i,ext}^n-F_{i,\mathrm{int}}^n\right]+2{\mathbf{u}}_i^n-{\mathbf{u}}_i^{n-1}, \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{F}\mathrm{E}\mathrm{M} \mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{s}\end{array}\right.$$

(13)

where $\mathrm{\Delta }t$ is the size of the time step. A stability condition derived by Silling and Askri [41] can be used to determine the time step size, $\mathrm{\Delta }t$ as

$$\mathrm{\Delta }t<\sqrt{\frac{2\rho }{{\displaystyle {\sum }_{j=1}^m\frac{c}{\left|\left(x_p-x_i\right)\right|}{V}_j}}}$$

(14)

Moreover, for the PD particles, the horizon size $\delta$ has a significant influence on the accuracy of the numerical simulations. The horizon size can be selected using the scale characteristics of the simulated object. In practice, $\delta =3\mathrm{\Delta }x$ usually works well [42]. Therefore, the horizon size is set to $\delta =3\mathrm{\Delta }x$.

The numerical code for the proposed PD-FEM coupling approach is compiled using Fortran 90. A flowchart of the PD-FEM coupling approach is shown in Figure 3.

4. Validation of PD-FEM Coupling Approach

4.1. Bending Deformation of Cantilever Beam

A three-dimensional cantilever beam subjected to a transverse loading of F = 0.64 N at the free end is examined, and the solutions given by the proposed coupling approach are compared with the FEM solutions. Because the bending deformation of a cantilever beam is a quasi-static problem, and to achieve a quantitative quasi-static calculation, the dynamic relaxation method is introduced to peridynamics [43].

The cantilever beam is 8 mm long, 2 mm wide, and 2 mm thick with Young’s modulus of 1.0 GPa, Poisson’s ratio of 0.25, and a density of 900 kg/m³. Figure 4 shows the PD-FEM coupling model of this cantilever beam, which is partitioned into one FEM subregion and one PD subregion. The FEM subregion consists of 16 hexahedral elements of size 1 mm × 1 mm × 1 mm, whereas the PD subregion is discretized uniformly into 16 × 8 × 8 = 1024 particles with the grid spacing $\mathrm{\Delta }x=0.25$ mm and horizon size $\delta =3\mathrm{\Delta }x$. Three layers of peridynamics nodes are embedded in four interface elements for the coupling force calculations.

Figure 5 shows the change in deflection along the central line of the beam. The displacement curve obtained by the numerical simulation is in good agreement with that of FEM, and its smoothness verifies the displacement coordination of the coupling approach. From Figure 5, it can be found that the ratio of the characteristic scale to the horizon size should be no less than 1 to obtain an acceptable result, however, more cases may be need to achieve this conclusion. The above results prove that the coupling algorithm achieves good accuracy and displacement coordination in the calculation of bending deformation, verifying the correctness of the proposed PD-FEM coupling approach in static problem.

Figure 4. PD-FEM coupling model of the cantilever beamThe simulated deflection at the free end of the beam and the coupling force are presented in

Table 1. Simulation results of bending beam.

	PD-FEM	FEM	Errors
Deflection	8.0492 × 10−5 m	8.19 × 10−5 m	1.74%
Force	−0.6327 N	−0.64 N	1.14%

. The deflection and force obtained from the proposed coupling approach are very close to the FEM results (error less than 2%), which indicates that the coupling approach transfers the force accurately.

Figure 4. PD-FEM coupling model of the cantilever beamThe simulated deflection at the free end of the beam and the coupling force are presented in Table 1. The deflection and force obtained from the proposed coupling approach are very close to the FEM results (error less than 2%), which indicates that the coupling approach transfers the force accurately.

4.2. Failure of 2D Plate with Central Crack

Mode-I crack is selected to simulate crack initiation and propagation along the plate. A 50 mm × 50 mm square plate with a 10 mm central pre-crack is stretched from both ends at a velocity of 50 mm/s, as shown in Figure 6.

The PMB material properties used in this example are as follows: Young’s modulus is E = 192 GPa, Poisson’s ratio is v = 0.33, mass density is 8000 kg/m³, and the critical stretch s₀ is 0.04472. These material parameters, geometry, and loading conditions are the same as the simulation example reported by Madenci and Oterkus [42]. For the coupling model, the plate is partitioned into one PD subregion and two FEM subregions (see Figure 6). The PD zone contains 100 × 200 = 20,000 particles, which are discretized regularly with a grid spacing of $\mathrm{\Delta }x=0.25$ mm and a horizon size of $\delta =3\mathrm{\Delta }x$. The FEM parts are composed of 32 four-node rectangular elements of size 6.25 mm × 6.25 mm. The interface region has three additional layers of peridynamics nodes (total of 2 × 3 × 200 = 1200 nodes). The time step is $\mathrm{\Delta }t=3.34\times {10}^{-8}$ s, which satisfies the stable time step condition.

For comparison, the PD solution is considered, as shown in Figure 6. The model size and load conditions are the same as for the coupling model, and the region is discretized regularly into 200 × 200 = 20,000 nodes with a grid spacing of $\mathrm{\Delta }x=0.25$ mm and horizon size of $\delta =3\mathrm{\Delta }x$. In the velocity boundary area, three virtual boundary layers are added, each with 3 × 200 = 600 nodes.

Figure 7 shows numerical simulation results of crack tracks using the PD-FEM coupling model and PD model (see Figure 7). Crack paths obtained from the proposed coupling method resemble the mode-I failure in brittle material, and are in good agreement with those obtained from the PD method. These results are similar to those reported by Madenci and Oterkus [42]. The displacements along the y-axis obtained from the PD-FEM and PD methods are plotted in Figure 8. We can find that they are in close agreement. Therefore, the PD-FEM coupling method can simulate crack propagation well, which verifies the correctness of the coupling approach in dynamic conditions.

5. PD-FEM Simulation of Icebreaker Navigation in Ice Level

5.1. Numerical Simulation

5.1.1. Ice Constitutive Model and Failure Criterion

Ice is a complex material that is affected by factors such as temperature, porosity, and grain size. Several laboratory tests have analyzed the brittle strength and failure patterns of ice, as well as the characteristics of ice–structure interactions [44,45]. The compressive strength and tensile strength of ice varies, with the compressive strength being around 3–4 times the tensile strength. The ductile–brittle transition of ice is another challenging topic, and ice shows different behavior at different strain rates (i.e., loading rate), as discussed by Schulson [45].

At low-speed loading rates, ice behaves as a ductile material, whereas at high-speed strain rates, the ice presents the characteristics of a linear elastic material, with a brittle mode when damage occurs. Therefore, ice is regarded as an elastic brittle material (PMB material in bond-based peridynamics) when interacting with an icebreaker, as the relatively high speed of icebreaker vessels corresponds to a high strain rate in the ice. Hence, a reasonable linear elastic constitutive model of ice is established for the bond-based peridynamics. Mechanical ice tests conducted by Schulson can be used to demonstrate the rationality of this linear elastic constitutive model.

In bond-based peridynamics, the bond force represents stress and the bond stretch represents strain. Let $s_t=s_0$ be the critical bond stretch and $\left|s_c\right|=4s_0$ be the critical bond compression. If the bond stretch exceeds $s_t$ in the tensile case or $s_c$ in the compressive case, the bond will break and the bond force becomes zero.

Based on the bond-based peridynamics theory [42], the critical bond stretch in 3D cases is defined as

$$s_0=\sqrt{\frac{5\pi G_0}{18E\delta }}$$

(15)

where $G_0$ is the energy release rate, which reflects the resistance of a material to crack propagation and can be derived from fracture mechanics. Linear elastic fracture mechanics is based on linear elastic theory and is applicable to brittle fractures. From the perspective of energy conservation, the condition for crack propagation is

$$G_0\ge G_C$$

(16)

where $G_C$ is the energy absorption rate. The primary fracture mode is tensile failures, because its compression strength is 3–4 times of tensile strength. $G_0$ can be expressed as

$$G_0=\frac{K_I^2}{E}$$

(17)

where $K_I$ is the fracture toughness, which reflects the resistance of a material to brittle fracture and can be measured experimentally. Therefore, the critical stretch can be calculated as

$$s_0=\sqrt{\frac{5\pi K_I^2}{18E^2\delta }}$$

(18)

5.1.2. The Gravity and Buoyancy Model of Ice

The gravity and buoyancy of the ice are in balance in still water. when interacting with a ship, the ice will deviate from its equilibrium position as the gravity and buoyancy become unbalanced. To simplify the influence of gravity and buoyancy, a body force density $b_z$ is introduced.

If a particle is completely under the waterline, we have

$$b_z=-g{\rho }_i+g{\rho }_w$$

(19)

If a particle is completely above the waterline, we have

$$b_z=-g{\rho }_i$$

(20)

For other particles, we have

$$b_z=-g{\rho }_i+g{\rho }_w{l}_w/d$$

(21)

where ${\rho }_i$ is the density of ice, ${\rho }_w$ is the density of water, $l_w$ is the length of particles immersed in water, and d is the size of particles.

5.1.3. Ship-Ice Contact Model

The hull is modeled with FEM. The ice sheet contacting with the hull is modeled with peridynamics. Therefore, the ship–ice contact can be transformed into the interaction between peridynamics and the FEM models. In this case, the contact model developed by Liu [30] is introduced to calculate the contact force between the triangular elements of FEM and the particles of PD. To determine contact between a particle and a triangular element, some specific points must be defined (see Figure 9), namely the vertexes of the triangular element A, B, C and the particle P.

The determination of contact is divided into two stages. In the first stage, we estimate whether the distance between the particle and the plane of triangle ABC is less than the particle radius $r=\mathrm{\Delta }x/2$. The distance is calculated with Equation (22), and the critical distance is defined by Equation (23)

$$\left|PQ\right|=PA\cdot \frac{BA\times CB}{\left|BA\times CB\right|}$$

(22)

$$\left|PQ\right|<r$$

(23)

where Q is the projection of P onto the plane of triangle ABC.

If the criterion in the first stage is not satisfied, the particle P does not contact the triangle element. When the first stage criterion is satisfied, it is necessary to further decide whether the projection point Q of P is inside the triangular region ABC, and this is the second stage. The centroid of the particle is checked to determine whether it is inside the triangle.

For any point Q in the plane ABC, the vector d = AQ can be expressed by two non-parallel vectors a = AC and b = AB in plane ABC as

$$\mathit{d}=u\mathit{a}+v\mathit{b}$$

(24)

where the coefficients u and v are defined as

$$u=\frac{\left(\mathit{b}\cdot \mathit{b}\right)\left(\mathit{d}\cdot \mathit{a}\right)-\left(\mathit{b}\cdot \mathit{a}\right)\left(\mathit{d}\cdot \mathit{b}\right)}{\left(\mathit{b}\cdot \mathit{b}\right)\left(\mathit{a}\cdot \mathit{a}\right)-\left(\mathit{b}\cdot \mathit{a}\right)\left(\mathit{a}\cdot \mathit{b}\right)}$$

(25)

$$v=\frac{\left(\mathit{a}\cdot \mathit{a}\right)\left(\mathit{d}\cdot \mathit{b}\right)-\left(\mathit{a}\cdot \mathit{b}\right)\left(\mathit{d}\cdot \mathit{a}\right)}{\left(\mathit{a}\cdot \mathit{a}\right)\left(\mathit{b}\cdot \mathit{b}\right)-\left(\mathit{a}\cdot \mathit{b}\right)\left(\mathit{b}\cdot \mathit{a}\right)}$$

(26)

If the projection point Q is inside the triangular element ABC, the two coefficients must satisfy conditions, as

$$\begin{array}{l}u\ge 0\\ v\ge 0\\ u+v\le 1\end{array}$$

(27)

If the criteria in these two stages are satisfied, the particle P is in contact with the triangular element. The contact force is defined by the repelling short-range force. The force between two particles is given as

$$f=-\frac{PQ}{\left|PQ\right|}\min \left\{0,c_{sh}\left(\frac{\left|PQ\right|}{r}-1\right)\right\}$$

(28)

where $c_{sh}$ is the short-range force coefficient, can be choose as $c_{sh}=5c$.

5.1.4. Numerical Model

This section describes numerical simulations of icebreaker navigation in level ice. Numerical models for the icebreaker and the level ice are illustrated in Figure 10.

Xuelong icebreaker is selected and its bow is modeled. Main parameters are listed in

Table 2. Main parameters of icebreaker.

Parameter	Variable	Value
Ship length	L	166.0 m
Ship breadth	B	22.6 m
Ship depth	D	13.5 m
Bow length	l	29.6 m
Bow breadth	b	22.6 m
Draft	T	8.0 m
Stem angle	α	20°
Flooding angle	β	24°
Ship-ice friction coefficient	μ	0.15

. The FEM model of the ship’s bow consists of 304 triangular elements. The ship is treated as a rigid body sailing in a straight line at a speed of 3 kn. The ice constitutive model is elastic-brittle, as described in Section 5.1.1. The ice sheet is a rectangle whose edges are fixed, except the one interacting with the icebreaker. The parameters of the ice sheet are presented in

Table 3. Parameters of level ice.

Parameter	Variable	Value
Young’s modulus	E	6.83 GPa
Poisson’s ratio	ν	0.25
Bending strength	σf	2.96 MPa
Fracture toughness	KI	115 kNm−3/2
Density	ρ	894 kg/m3
Area	A	100 × 100 m2
Thickness	h	1.0 m

, and the critical bond stretch is calculated using Equation (18). The level ice is modeled using the proposed PD-FEM coupling approach, with PD and FEM subregions. In the PD subregion (width = 40 m, length = 70 m), the grid spacing is $\mathrm{\Delta }x=0.25$ m and the horizon radius is $\delta =3\mathrm{\Delta }x$. In the interface between the PD and FEM subregions, three layers of peridynamics nodes are embedded in each edge. Therefore, there are approximately 187,840 PD nodes. The FEM subregion is discretized into 1800 hexahedral elements of size 2 × 2 × 1 m and 3936 nodes. The size of time step is $\mathrm{\Delta }t=1.0\times {10}^{-7}$ s, which is satisfies with the stability time step condition. The total time steps are 300,000,000 steps.

5.1.5. Numerical Result and Discussion

The numerical simulation of the ice-breaking process of an icebreaker navigating through level ice is illustrated in Figure 11. Xuelong model test is performed in the ice mechanics and ice engineering laboratory of Tianjin University to observe the failure model of level ice and the motion of broken ice. After the test is finished, Xuelong model is dragged to slowly retreat by the main trailer, and then a more complete ice breaking area can be shown on the water surface.

When the icebreaker first contacts the ice, the ice is subjected to a compressive force from the bow tip, mainly along the x-axis. Once the force is greater than the ice critical strength, ice particles fall from the level ice, and a notch like the tip of the bow is formed (see Figure 11a). As the ship moves, the notch expands. The three-dimensional curved hull makes the ice bear the force in three directions. The force along the y-axis causes the ice to tear, and the notch along the profile forms a radial crack, as shown in Figure 11b.

As the contact area increases, the force along the z-axis exceeds the ice critical strength, resulting in ice sheet bending deformation and failure. Accordingly, circular cracks form, as shown in Figure 11c; these are like the circumferential cracks in model test, as shown in Figure 12.

With the ship continuous movement, the contact force increases in all three directions and radial cracks and circular cracks propagate, which is in good agreement with the ice sheet failure model observed in Xuelong model test as shown in Figure 12. The front segments of the level ice form a wedge shape, as shown in Figure 11d. The further movement of the ship causes wedge ice to form on the two sides of the bow, like the real phenomenon observed in Xuelong model test. Subsequently, ice blocks fall off, flip, push away, and pile up, as shown in Figure 11g,h.

Numerical simulation results in the generation and propagation of radial and circular cracks, as well as phenomena such as the shedding of wedge ice, flipping of brash ice, and cleaning of the channel, which are in broad agreement with experimental and real phenomena. In addition, the coupling method of finite element and perdynamics can effectively suppress the boundary effect when the level ice is failure, compared with bond-based peridynamics [28,29,30,31,32,33,34,35,36,37].

To show the computational efficiency of PD-FEM coupling approach, the PD solution is considered. Except the level ice is completely modeled by PD solution, loading conditions are the same as the PD-FEM coupling approach.

Table 4. Comparison of PD and PD-FEM coupling approach in computational efficiency.

Item	PD-FEM	PD
Particle number	187,840	654,400
Element number	1800	0
Total time steps	300,000	300,000
Total CPU time	52.32 h	135.47 h

shows the comparison results of PD and PD-FEM coupling method in computational time. Both two methods are compiled using Fortran 90, and use CPU_TIME function in Fortran language to calculate the program running time. In the case of the same equipment condition and 300,000 steps of the calculation time steps, the PD method needs calculate 135.47 h, while the PD-FEM coupling approach only needs 52.32 h. Compared with PD model, PD-FEM coupling model is 2.59 times more efficient and reduces the calculation time.

5.2. Influence of Ship Speed on Ice Load

To further validate sailing speed influence on the ice load, six simulation sets are selected as speeds with 2, 3, 4, 5, 6 and 7 kn with an ice thickness of 1 m. The ice load results are shown in Figure 13. The average ice force during the period of ship-ice interaction is calculated and compared with Lindqvist’s empirical formula as shown in Figure 14.

Lindqvist’s empirical formula divides the icebreaking resistance into crushing at the stem, breaking by bending and submersion resistance. When ice is broken, the crushed ice can be cleared from both sides of the ship, and the opened channel will be larger than the ship’s width. Therefore, the ice crushing at the stem and breaking by bending parts of ice resistance are mainly caused by the bow. Therefore, the ice load obtained from numerical simulation is compared with the crushing at the stem and breaking by bending resistance calculated by Lindqvist’s empirical formula.

Figure 13 shows that the ice load curves (blue solid lines) are periodic. When the icebreaker interacts with the level ice, the ice load increases, but when wedge ice falls from the level ice and the icebreaker is not in contact with the ice, the ice load decreases. This is a reason for impulsive ice loads. Figure 13 also indicates that, as the ship speed increases, the period of the ice load becomes shorter and the ship velocity influences the peak of ice load. The ice load shows a rise trend as the velocity increases. From Figure 14, although there are some differences between the mean ice load results and those obtained from Lindqvist’s empirical formula, they are generally in a good agreement (see

Table 5. Ice resistance calculated by PD-FEM coupling model and Lindqvist empirical formula.

Method	2 kn	3 kn	4 kn	5 kn	6 kn	7 kn
Lindqvist	2.370 × 106 N	2.743 × 106 N	3.120 × 106 N	3.490 × 106 N	3.863 × 106 N	4.237 × 106 N
PD-FEM	2.202 × 106 N	3.921 × 106 N	3.008 × 106 N	4.612 × 106 N	4.307 × 106 N	4.560 × 106 N
Errors	7.1%	42.9%	3.6%	32.1%	11.5%	7.6%

and Figure 14).

5.3. Influence of Ice Thickness on Ice Load

Numerical simulations of an icebreaker moving at a fixed speed are performed with different ice thicknesses. The ice thickness is selected as 0.5 m, 0.75 m, and 1 m, and the ship’s speed is set to 3 kn.

The results of the ice load for different ice thicknesses are shown in Figure 15. It can be found that ice load presents an increasing trend in general. The mean ice load is 0.978 MN, 1.569 MN and 3.921 MN at ice thicknesses of 0.5 m, 0.75 m and 1.0 m, respectively. From Figure 16, we can find that the numerical results are in well agreement with those from Lindqvist’s empirical formula.

6. Conclusions

The coupling model of peridynamics with the finite element method is employed to simulate ship–ice interaction. The characteristics of the ice-breaking scenarios and the ice load are captured successfully. From the simulation results, the following conclusions can be drawn.

(1)

The PD-FEM coupling model can successfully simulate the generation and propagation of radial and circular cracks in level ice, as well as the phenomena of wedge ice shedding, broken ice flipping, and ice cleaning of the channel during the ice-breaking process.

(2)

Compared with bond-based peridynamics, the PD-FEM coupling model has better computational efficiency, and can effectively suppress the boundary effect when the level ice is failure.

(3)

The ice load obtained from the PD-FEM coupling model is in good agreement with that obtained from Lindqvist’s empirical formula.