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

: In this work, the ﬁnite 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, ﬂipping of brash ice, and cleaning of the channel are observed to be broadly consistent with experimental observation. The inﬂuence 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.


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 fullscale 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 where H x is the domain of integration within the horizon of particle x, u is the displacement vector field, and b is the body force density.ρ is the mass density, and f is a pairwise force density function defined as the force per unit volume that particle x ' 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 ξ = x − x and η = u(x , t) − u(x, t), respectively.Therefore, the relative position of the two interacting particles at t in the current configuration is ξ + η and the pairwise force density function can be described as f(η, ξ).
For the prototype micro-elastic brittle (PMB) material defined by Silling and Askari [41], the pairwise force density function can be expressed as where c = 12E/πδ 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 When the deformation stretch s exceeds a limit s 0 (described as the critical stretch for failure), the bond between the two particles breaks and no pairwise force remains.The term µ(t, η, ξ) is a history-dependent scalar-valued function, which is introduced to represent the bond failure of two particles.This can be defined as µ(t, η, ξ) = 1, s < s 0 0, s ≥ s 0 (4) Accordingly, the level of damage is illustrated by the local damage at one particle, defined as 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 0 can be obtained from the energy release rate.

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.
Accordingly, the level of damage is illustrated by the local damage at one particle, defined as 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 s0 can be obtained from the energy release rate.

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 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 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, (ξ , η , ψ ) 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 where F i,ext is the external force and the internal force is given by FEM nodes on segment FEM nodes not on segment, The displacements of the embedded peridynamics nodes are determined by where (ξ, η, ψ) are the natural coordinates of an embedded peridynamics node in the interface element and u i is the nodal displacement of an interface element.

Numerical Implementation
The peridynamics equation of motion after discretization is written as For the FEM subregion, the equation of motion for the FEM nodes is written as 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 ..
, for FEM nodes (13) where ∆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, ∆t as Moreover, for the PD particles, the horizon size δ 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, δ = 3∆x usually works well [42].Therefore, the horizon size is set to δ = 3∆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

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

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 3 .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 ∆x = 0.25 mm and horizon size δ = 3∆x.Three layers of peridynamics nodes are embedded in four interface elements for the coupling force calculations.1.The deflection and force obta the proposed coupling approach are very close to the FEM results (error less than 2% indicates that the coupling approach transfers the force accurately.Figure 5 shows the change in deflection along the central line of the be displacement curve obtained by the numerical simulation is in good agreement w of FEM, and its smoothness verifies the displacement coordination of the approach.From Figure 5, it can be found that the ratio of the characteristic sca horizon size should be no less than 1 to obtain an acceptable result, however, m may be need to achieve this conclusion.The above results prove that the algorithm achieves good accuracy and displacement coordination in the calcu bending deformation, verifying the correctness of the proposed PD-FEM approach in static problem.  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.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.

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.
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.

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 3 , and the critical stretch s0 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 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.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 3 , and the critical stretch s 0 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 ∆x = 0.25 mm and a horizon size of δ = 3∆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 ∆t = 3.34 × 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 ∆x = 0.25 mm and horizon size of δ = 3∆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.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 |s c | = 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 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 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 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

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 If a particle is completely above the waterline, we have For other particles, we have where ρ i is the density of ice, ρ w is the density of water, l w is the length of particles immersed in water, and d is the size of particles.

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.

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 2 rx = .The distance is calculated with Equation ( 22), and the critical distance is defined by Equation ( 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.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 = ∆x/2.The distance is calculated with Equation ( 22), and the critical distance is defined by Equation ( 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 nonparallel vectors a = AC and b = AB in plane ABC as where the coefficients u and v are defined as If the projection point Q is inside the triangular element ABC, the two coefficients must satisfy conditions, as 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 where c sh is the short-range force coefficient, can be choose as c sh = 5c.

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.
If the criteria in these two stages are satisfied, the particle P is in contact triangular element.The contact force is defined by the repelling short-range fo force between two particles is given as min 0,  Xuelong icebreaker is selected and its bow is modeled.Main parameters are Table 2.The FEM model of the ship's bow consists of 304 triangular elements.Th treated as a rigid body sailing in a straight line at a speed of 3 kn.The ice con model is elastic-brittle, as described in Section 5.1.1.The ice sheet is a rectangl edges are fixed, except the one interacting with the icebreaker.The parameters o Xuelong icebreaker is selected and its bow is modeled.Main parameters are listed in Table 2.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, 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 ∆x = 0.25 m and the horizon radius is δ = 3∆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 ∆t = 1.0 × 10 −7 s, which is satisfies with the stability time step condition.The total time steps are 300,000,000 steps.

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].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  When the icebreaker first contacts the ice, the ice is subjected to a compressive forc from the bow tip, mainly along the x-axis.Once the force is greater than the ice critic strength, ice particles fall from the level ice, and a notch like the tip of the bow is forme (see Figure 11a).As the ship moves, the notch expands.The three-dimensional curve hull makes the ice bear the force in three directions.The force along the y-axis causes th 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 critic strength, resulting in ice sheet bending deformation and failure.Accordingly, circula cracks form, as shown in Figure 11c; these are like the circumferential cracks in model tes as shown in Figure 12.With the ship continuous movement, the contact force increases in all three direction and radial cracks and circular cracks propagate, which is in good agreement with the ic sheet failure model observed in Xuelong model test as shown in Figure 12.The fron segments of the level ice form a wedge shape, as shown in Figure 11d.The furthe movement of the ship causes wedge ice to form on the two sides of the bow, like the re phenomenon observed in Xuelong model test.Subsequently, ice blocks fall off, flip, pus away, and pile up, as shown in Figure 11g,h.
Numerical simulation results in the generation and propagation of radial and circula cracks, as well as phenomena such as the shedding of wedge ice, flipping of brash ice, an cleaning of the channel, which are in broad agreement with experimental and re 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 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.

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 and Figure 14).

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.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 and Figure 14).

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.

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

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

22 Figure 2 .Figure 2 .
Figure 2. coupling scheme that divides a coupling force fp to FEM nodes on the interface segment.

Figure 3 .
Figure 3. Flow chart of the FEM-PD coupling scheme.

Figure 3 .
Figure 3. Flow chart of the FEM-PD coupling scheme.

J
Figure 4. PD-FEM coupling model of the cantilever beamThe simulated deflection at the f the beam and the coupling force are presented in Table1.The deflection and force obta the proposed coupling approach are very close to the FEM results (error less than 2% indicates that the coupling approach transfers the force accurately.

Figure 4 .
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 Table1.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 5 .
Figure 5. Vertical displacement on central axis of the beam.

Figure 5 .
Figure 5. Vertical displacement on central axis of the beam.J. Mar.Sci.Eng.2023, 11, x FOR PEER REVIEW 8 of 22

Figure 6 .
Figure 6.Geometry and loading condition of a plate with a central crack.
).The PD zone contains 100 × 200 = 20,000 particles, which are discretized regularly with a grid spacing of 0.25 x = mm and a horizon size of 3 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 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 0.25 x = mm and horizon size of 3 x  =.In the velocity boundary area, three virtual boundary layers are added, each with 3 × 200 = 600 nodes.

Figure 6 .
Figure 6.Geometry and loading condition of a plate with a central crack.


is the density of ice, w  is the density of water, w l is the length of particles immersed in water, and d is the size of particles.

Figure 9 .
Figure 9. Diagram of ship-ice search model.

Figure 9 .
Figure 9. Diagram of ship-ice search model.

c
is the short-range force coefficient, can be choose as 5This section describes numerical simulations of icebreaker navigation in l Numerical models for the icebreaker and the level ice are illustrated in Figure10

Figure 11 .
Figure 11.Ice-breaking process of icebreaker navigating in level ice simulated by PD-FEM coupling model.

Figure 12 .
Figure 12.Picture of level ice failure in the model test of Xuelong model at Tianjin University.

Figure 11 .
Figure 11.Ice-breaking process of icebreaker navigating in level ice simulated by PD-FEM coupling model.

Figure 11 .
Figure 11.Ice-breaking process of icebreaker navigating in level ice simulated by PD-FEM couplin model.

Figure 12 .
Figure 12.Picture of level ice failure in the model test of Xuelong model at Tianjin University.

Figure 12 .
Figure 12.Picture of level ice failure in the model test of Xuelong model at Tianjin University.

J 22 (Figure 13 .
Figure 13.Time history of ice load with different ship velocity.

Figure 13 .
Figure 13.Time history of ice load with different ship velocity.

Figure 14 .
Figure 14.Ice resistance obtained by Lindqvist formula and simulated with FEM-PD coupling model for different speeds.

Figure 14 .
Figure 14.Ice resistance obtained by Lindqvist formula and simulated with FEM-PD coupling model for different speeds.

Figure 15 .
Figure 15.Time history of ice load with different ice thicknesses.

Figure 16 .
Figure 16.Ice resistance obtain by Lindqvist formula and FEM-PD coupling model for different ice thicknesses.

Figure 15 . 22 (Figure 15 .
Figure 15.Time history of ice load with different ice thicknesses.

Figure 16 .
Figure 16.Ice resistance obtain by Lindqvist formula and FEM-PD coupling model for different ice thicknesses.

Figure 16 .
Figure 16.Ice resistance obtain by Lindqvist formula and FEM-PD coupling model for different ice thicknesses. .

Table 1 .
Simulation results of bending beam.

Table 1 .
Simulation results of bending beam.

Table 2 .
Main parameters of icebreaker.

Table 3 .
Parameters of level ice.

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

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

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