Updated Smoothed Particle Hydrodynamics for Simulating Bending and Compression Failure Progress of Ice

: In this paper, an updated Smoothed Particle Hydrodynamics (SPH) method based on the Simpliﬁed Finite Difference Interpolation scheme (SPH_SFDI) is presented to simulate the failure process of ice. The Drucker–Prager model is embedded into the SPH code to simulate the four point bending and uniaxial compression failure of ice. The cohesion softening elastic–plastic model is also used in the SPH_SFDI framework. To validate the proposed modeling approach, the numerical results of SPH_SFDI are compared with the standard SPH and the experimental data. The good agreement demonstrated that the proposed SPH_SFDI method including the elastic–plastic cohesion softening Drucker–Prager failure model can provide a useful numerical tool for simulating failure progress of the ice in practical ﬁeld. It is also shown that the SPH_SFDI can signiﬁcantly improve the capability and accuracy for simulating ice bending and compression failures as compared with the original SPH scheme.


Introduction
With the increasing activities in Arctic regions, the method for the accurate calculation of the corresponding ice loads on structures are crucial to the design of marine structures operating in ice [1]. To simulate ice-water or ice-ship interactions effectively, it is necessary to have a reasonable study and understanding of the ice failure progress. The bending failure is the common failure behavior of the ice and is important for ships in ice-ship interactions because of the inclined contact interfaces with the ice [2,3]. In addition, under the compression of a ship or structure, compressive failure (crushing) is also easy to occur during the failure process of the ice. Thus, it is of high importance to study the bending and compression failure progress of the ice.
In the past few years, many full-scale tests and model tests of the ice failure have been investigated [4][5][6][7][8]. However, the experimental data are highly dispersed because of different experimental equipment, different test methods, and different measurements of the ice specimens. In addition, some simplified empirical models are also used to study the ice failure and the interaction of ice and structure [9][10][11]. However, these simplified empirical models only focus on some main aspects of ice failure and are lack of the study of dynamics and some changing details during the ice failure progress. Thus, it is very important to develop a reliable numerical ice model to simulate the sea ice failure in the bending and compression process, especially the current studies on the behavior of sea ice failure are not adequate; although some obtained numerical approaches were proposed to where α and β are the Cartesian components in x, y and z directions; v is the particle velocity; ρ is the ice density; σ αβ is the tress tensor of ice particles; g is the gravitational acceleration; and D/Dt is the particle derivative following motion. The ice constitutive relation need to be applied into the system to solve the governing equations (Equations (1) and (2)). In this paper, the stress tensor can be divided into two parts, which is same with [32] and includes the hydrostatic pressure and deviatoric shear stress: in which δ αβ = Kronecker delta and satisfies the following conditions: δ αβ = 1 if α = β or δ αβ = 0 when α = β.

Ice Elasto-Plastic Constitutive Model
To simulate the ice failure process, an elasto-plastic constitutive model [33] is applied into SPH in this paper. The components of the strain rate . ε αβ are given by: For an elasto-plastic material, the strain rate λ ∂Q ∂σ αβ (6) where . λ is the plastic multiplier rate, and Q is the plastic potential function which determines the development direction of plastic strain. The plastic multiplier . λ is computed through the consistency condition, which is given by: According to Equations (5) and (6), the total strain rate tensor can be expressed as: .
The plastic multiplier rate . λ of an elasto-plastic material can be calculated by substituting Equation (9) into Equation (7)

Drucker-Prager Model
The Drucker-Prager yield criterion ( Figure 1) has been widely used in soil and rock mechanics but is almost unknown in the ice field. In this paper, the Drucker-Prager yield criterion with flow rules has been prospectively used to determine the plastic regime of the ice. The validation of numerical results in Section 5 shows that the Drucker-Prager yield criterion can be used to identify the occurrence of the plastic deformation of ice particles in SPH.
Water 2017, 9,882 4 of 24 in which α and β are free indexes, m and n are dummy indexes which denote the Cartesian components x , y with the Einstein convention applied to repeated indices; is the deviatoric shear strain rate tensor; is the elastic bulk modulus; and is the shear modulus.
The plastic multiplier rate λ  of an elasto-plastic material can be calculated by substituting Equation (9) into Equation (7) as follows: The Drucker-Prager yield criterion ( Figure 1) has been widely used in soil and rock mechanics but is almost unknown in the ice field. In this paper, the Drucker-Prager yield criterion with flow rules has been prospectively used to determine the plastic regime of the ice. The validation of numerical results in Section 5 shows that the Drucker-Prager yield criterion can be used to identify the occurrence of the plastic deformation of ice particles in SPH.

Drucker-Prager Model
In this study, the Drucker-Prager yield criterion can be expressed as following: in which c is the ice cohesion, 2 ( ) J s αβ is the second invariant of the stress tensor, and 1 ( ) I αβ σ is one third of the first invariant of the stress tensor. The parameters φ α and ξ are defined as: where φ is the friction angle. In addition, the plastic potential function is also used to completely In this study, the Drucker-Prager yield criterion can be expressed as following: in which c is the ice cohesion, J 2 (s αβ ) is the second invariant of the stress tensor, and I 1 (σ αβ ) is one third of the first invariant of the stress tensor. The parameters α φ and ξ are defined as: Water 2017, 9, 882 5 of 24 where φ is the friction angle. In addition, the plastic potential function is also used to completely define the relationship between the stress and strain. The flow rules are usually applied into SPH to simulate solid fracture, including the associated flow rule and non-associated flow rule. In the associated flow rules, the plastic potential function has the same form with the yield criterion, namely as: The non-associated plastic potential function is taken to be: where parameter η is related to dilatancy angle ϕ, which can be expressed as: Substituting Equation (13) into Equations (9) and (10), the stress-strain relationship with the associated plastic flow rule is given by: When F(σ αβ , c) < 0, it is in pure elasticity condition: where the plastic multiplier rate . λ is calculated for the ice model by: The stress-strain equation of the ice model with the non-associated flow rule is obtained by taking Equation (14) into Equations (9) and (10) as follows: where the plastic multiplier rate . λ can be written as: It can be seen from the above description that the main difference between the associative and non-associative models is reflected in the dilatancy angle. In the associated flow rule ice model, the Dilatancy angle is always equal to the friction angle, whereas dilatancy angle is optional in the non-associated flow rule. It should be noted that according to the comparative analysis in the following Section 5, non-associative flow rule yield more stable and precise numerical results than associative flow rule in this paper.

Numerical Errors in Computational Plasticity
In computation for plastic deformation of elastic-plastic material using the Drucker-Prager yield criterion, the numerical errors are easy to occur, which corresponds to the following condition: In this study, a stress-rescaling procedure based on Bui et al. [33] is adopted to modify the stress. The stress components are modified according to the following relation: This scaling factor r n at the time step n is defined by: In addition, if the condition −α φ I n 1 + ζc < 0 is satisfied at the time step n, the normal stress components need to be adjusted to the new correct values σ αβ : When −α φ I n 1 + k c < J n 2 is satisfied, the stress tensor needs to take the plastic correction, which can be expressed as:

Cohesion Softening
In this paper, the cohesion softening law [36] needs to be used in the Drucker-Prager model to simulate the reduction of the ice strength under external loading numerically. In addition, the cohesion softening model can imply the time dependency of ice failure, which is validated in Section 5.2. The model of cohesion softening is realized by making cohesion c a purely linear function of the accumulated plastic strain (Figure 2), which is similar to: The specific cohesion softening law in this paper is shown as: k is the specific softening coefficient and c R is the minimum cohesion. The accumulated plastic strain ε p can be obtained by the associative softening law as: λξ (28) Water Because the relationship between cohesion and accumulated plastic strain is a purely mathematical construct, it is difficult to obtain an exact characterization of this relationship. According to Figure 3, different cohesion softening laws can get different results of the cohesion softening and the stress-strain relationships, which can make the simulation of widespread material failure behaviors possible, and can include both brittle and ductile failure. The higher order mathematical equation for cohesion softening law may need to simulate more complex and precise material failure behaviors. More details about the cohesion softening law can be seen in Whyatt and Board [36].

The Particle Approximation and Spatial Derivatives of SPH
In SPH method, the computational domain is discretized into a set of particles which carry some variables such as pressure, stress, velocity, density, etc. The smoothing kernels are used to approximate a continuous flow field. The basic principle of SPH expression is that, for any quantity of particle i , whether a scalar or a vector, it can be approximated by the direct summation of the relevant quantities of its neighbor particles j , which is shown as: Because the relationship between cohesion and accumulated plastic strain is a purely mathematical construct, it is difficult to obtain an exact characterization of this relationship. According to Figure 3, different cohesion softening laws can get different results of the cohesion softening and the stress-strain relationships, which can make the simulation of widespread material failure behaviors possible, and can include both brittle and ductile failure. The higher order mathematical equation for cohesion softening law may need to simulate more complex and precise material failure behaviors. More details about the cohesion softening law can be seen in Whyatt and Board [36]. Because the relationship between cohesion and accumulated plastic strain is a purely mathematical construct, it is difficult to obtain an exact characterization of this relationship. According to Figure 3, different cohesion softening laws can get different results of the cohesion softening and the stress-strain relationships, which can make the simulation of widespread material failure behaviors possible, and can include both brittle and ductile failure. The higher order mathematical equation for cohesion softening law may need to simulate more complex and precise material failure behaviors. More details about the cohesion softening law can be seen in Whyatt and Board [36].

The Particle Approximation and Spatial Derivatives of SPH
In SPH method, the computational domain is discretized into a set of particles which carry some variables such as pressure, stress, velocity, density, etc. The smoothing kernels are used to approximate a continuous flow field. The basic principle of SPH expression is that, for any quantity of particle i , whether a scalar or a vector, it can be approximated by the direct summation of the relevant quantities of its neighbor particles j , which is shown as:

The Particle Approximation and Spatial Derivatives of SPH
In SPH method, the computational domain is discretized into a set of particles which carry some variables such as pressure, stress, velocity, density, etc. The smoothing kernels are used to approximate a continuous flow field. The basic principle of SPH expression is that, for any quantity of particle i, whether a scalar or a vector, it can be approximated by the direct summation of the relevant quantities of its neighbor particles j, which is shown as: and its gradient can be shown as: where i and j are the referred particle and its neighbor, respectively, and W r ij is a kernel function and has different forms. In this paper, the cubic B-spline kernel proposed by Monaghan and Lattanzio [38] was used: where q = r/h, α d = 15/ 7πh 2 for 2D cases, and h is equal to 1.2-1.4dx (dx is the initial particle spacing).
In SPH, the mass conservation Equation (1) can be approximated as follows: where ρ i is the density of particle i with velocity component v i ; and m j is the mass of particle j which has velocity component v j . The most widely used SPH approximation of the momentum equation (Equation (2)) is: where Π ij is the artificial viscosity, which was proposed by Monaghan [39]. Finally, the position of particle i in SPH is calculated based on the following equation: In addition, the XSPH method [39] is used to solve problems involving the tension. In this method, particle i is defined based on an average velocity, which is shown as:

Artificial Stress Method
An artificial stress method presented by Monaghan [40] and Gray et al. [41] was used in many papers to remove numerical instability [42] caused by the clumping of SPH particles when SPH is Water 2017, 9, 882 9 of 24 applied to solid mechanics. This method adopts an artificial repulsive force. The artificial repulsive force proposed in Gray et al. [41] is used in this paper and takes the form: where n is the variable exponent based on the smoothing kernel. f ij is defined as where ∆d is the initial distance between neighbor particles. h is set to be 1.2∆d for the cubic B-spline kernel in this paper. The R αβ i and R αβ j in Equation (36) is the artificial stress tensor of particles i and j, respectively, with the correction parameter ε (Gray et al. [41]) : The same rule applies for R ββ i with αα replaced by ββ. Where σ αα i and σ ββ i are the new components of the stress tensor in the rotated frame: where c = cos θ i and s = sin θ i . θ i is the angle of roiration for particle i, which statisfies More details about the artificial stress can be found in Gray et al. [41]. For the tests discussed in this study, the parameter ε and n are equal to 0.3 and 4, respectively, to solve the tensile instability problems in SPH.

Boundary Conditions
In this paper, we deal with boundary conditions by two types of particles: solid boundary particles and mirror particles.
The solid boundary is fixed by the particles, which may prevent the real ice particles from penetrating the solid wall ( Figure 4). The boundary particles contribute to the velocity and stress gradients for the real ice particles near the boundary. These boundary particles have the same velocity density as the solid wall and their density is set equal to reference density. The stresses of the boundary particles on the solid boundary are calculated by using: where σ αβ w is the stress of the particle w on a boundary solid boundary; i is its neighboring particle and i can only be the real ice particle; and N is the number of particles in the support domain of wall boundary particle w. . i θ is the angle of roiration for particle i ,which statisfies More details about the artificial stress can be found in Gray et al. [41]. For the tests discussed in this study, the parameter ε and n are equal to 0.3 and 4, respectively, to solve the tensile instability problems in SPH. In this paper, we deal with boundary conditions by two types of particles: solid boundary particles and mirror particles.

Boundary Conditions
The solid boundary is fixed by the particles, which may prevent the real ice particles from penetrating the solid wall ( Figure 4). The boundary particles contribute to the velocity and stress gradients for the real ice particles near the boundary. These boundary particles have the same velocity density as the solid wall and their density is set equal to reference density. The stresses of the boundary particles on the solid boundary are calculated by using: In addition, the mirror particle ( Figure 4) method following Libersky and Petschek [43] is also used to simulate the solid boundary with the free-slip condition. For each real particle i that is close to the wall, a mirror particle i mir is set by a direct reflection of particle i across the boundary. The mirror particle i mir has the same tangential velocity (v i mir ,t ) with that of real particle: v i mir ,t = v i,t to simulate the free-slip boundary condition. The normal velocity (v i mir ,n ) of i mir is set opposite to that of real particle v i mir ,n = v i,n to prevent the real particles from penetrating the boundary as shown in Figure 4. The density and stress tensors of mirror particles are set to be equal to those of real ice particles.

Corrective SPH Method
The strain rate of the tensor Equation (4) needs to be converted into the discrete form to get the stress rate based on the generalized Hooke's law. In standard SPH, the strain rate is obtained by: The standard SPH algorithm is lack of high accuracy due to kernel approximation of its first-order derivative, such as Equation (46). To overcome the shortcomings in first order derivative accuracy of the original SPH, this paper adopts the Simplified Finite Difference Interpolation (SPH_SFDI) method to calculate the strain rate of the ice particles, more details about SFDI method can be found in Ma [37]. According to the results in Zheng et al. [44], SFDI can be a very good option as a high order accuracy. For the purpose of the completion of theory, the formulas of strain rate of the tensor in 2D case can be shown as: where in which α = x, y, β = x, y and m = x, k = y or m = y, k = x, and N is the neighbor particle number of particle i, r m j is the component of the position vector in x or y direction. Similarly, the derivative of other variables can also be calculated by this corrective method.
To justify that the SFDI method is more effective than the standard SPH method for the strain rate calculation. Figure 5 shows the comparison of the bending stress in the middle of the ice beam for four-point bending of the ice beam which will be discussed in Section 5.2. In Figure 5, the standard formula is referred to Equation (46) and the SFDI is referred to Equation (47).  (46)) and the SFDI scheme (Equation (47)).
According to the comparison of Figure 5, the results from SFDI scheme can get better agreement with the theoretical value than the ones of traditional equation. Especially when the fracture failure start to occur in the ice beam at about t = 0.45 s, the stress values by standard formula deviates from theoretical results obviously. The source of discrepancy is expected to be that the accuracy of strain rate in the standard formula is less than the ones in the SFDI scheme.
In SPH_SFDI, the main procedures of numerical implementation of failure model of ice are shown as follows: (1) Calculate the values of αβ ε and αβ σ from Equations (47), (17) or (19).
(2) Calculate the stress components αβ σ based on the obtained stress rate αβ σ .
(3) Check the stress state and judge whether the corresponding stress need to be corrected: if , the stress need to be modified by Equation (25).

Numerical Simulations
In this section, we firstly use the elastic vibration of a cantilever beam to verify the feasibility of SPH_SFDI method in solid mechanics. To test the effectiveness of the SPH_SFDI for simulating the failure progress of ice, two typical tests are included: the ice four-point bending and uniaxial ice compressive test. The enhanced performance of the SPH_SFDI algorithm will be demonstrated through the quantitative comparisons with the standard SPH and experimental data.

Elastic Vibration of a Cantilever Beam
The elastic vibration of a cantilever beam is used as a benchmark test to verify the reliability of the SPH_SFDI model for the calculation of solid mechanics. The cantilever beam is shown in Figure  6, the dynamic load P is acting at the free end of the cantilever beam. The length L = 48 m, the height is D = 12 m, the elastic modulus is E = 3.0 × 10 7 N/m 2 , the Poisson's ratio is , and the mass density is P = 1 kg/m 3 . External excitation force P = 1000 ( ) g t and ( ) g t is a function related to time.  (46)) and the SFDI scheme (Equation (47)).
According to the comparison of Figure 5, the results from SFDI scheme can get better agreement with the theoretical value than the ones of traditional equation. Especially when the fracture failure start to occur in the ice beam at about t = 0.45 s, the stress values by standard formula deviates from theoretical results obviously. The source of discrepancy is expected to be that the accuracy of strain rate in the standard formula is less than the ones in the SFDI scheme.
In SPH_SFDI, the main procedures of numerical implementation of failure model of ice are shown as follows: (1) Calculate the values of . ε αβ and . σ αβ from Equations (47), (17) or (19).
(2) Calculate the stress components σ αβ based on the obtained stress rate . σ αβ .
(3) Check the stress state and judge whether the corresponding stress need to be corrected: if −α φ I n 1 + k c < J n 2 , the stress need to be modified by Equation (25). (4) Implement Cohesion softening model based on Equation (27).

Numerical Simulations
In this section, we firstly use the elastic vibration of a cantilever beam to verify the feasibility of SPH_SFDI method in solid mechanics. To test the effectiveness of the SPH_SFDI for simulating the failure progress of ice, two typical tests are included: the ice four-point bending and uniaxial ice compressive test. The enhanced performance of the SPH_SFDI algorithm will be demonstrated through the quantitative comparisons with the standard SPH and experimental data.

Elastic Vibration of a Cantilever Beam
The elastic vibration of a cantilever beam is used as a benchmark test to verify the reliability of the SPH_SFDI model for the calculation of solid mechanics. The cantilever beam is shown in Figure 6, the dynamic load P is acting at the free end of the cantilever beam. The length L = 48 m, the height is D = 12 m, the elastic modulus is E = 3.0 × 10 7 N/m 2 , the Poisson's ratio is υ = 0.3, and the mass density is P = 1 kg/m 3 . External excitation force P = 1000g(t) and g(t) is a function related to time.

Elastic Vibration of a Cantilever Beam
The elastic vibration of a cantilever beam is used as a benchmark test to verify the reliability of the SPH_SFDI model for the calculation of solid mechanics. The cantilever beam is shown in Figure  6, the dynamic load P is acting at the free end of the cantilever beam. The length L = 48 m, the height is D = 12 m, the elastic modulus is E = 3.0 × 10 7 N/m 2 , the Poisson's ratio is , and the mass density is P = 1 kg/m 3 . External excitation force P = 1000 ( ) g t and ( ) g t is a function related to time. A simple harmonic load g(t) = sin ωt is considered. ω is the frequency of harmonic load and in this case ω = 27 s −1 . Figure 7 shows the comparison of the displacement in y direction of the free end of the cantilever beam (y) between the SPH and SPH_SFDI results with 10,000 particles and the finite element method (FEM) solution from Long [45]. This shows that the displacement time histories computed by the SPH_SFDI method shares a better agreement with the FEM data than the SPH result.  Figure 7 shows the comparison of the displacement in y direction of the free end of the cantilever beam ( y ) between the SPH and SPH_SFDI results with 10,000 particles and the finite element method (FEM) solution from Long [45]. This shows that the displacement time histories computed by the SPH_SFDI method shares a better agreement with the FEM data than the SPH result. To evaluate the enhanced performance of SPH_SFDI method further, the convergence properties of the SPH and SPH_SFDI models are now examined in terms of the displacement y . For this purpose, the time histories of displacement computed by SPH and SPH_SFDI are presented in Figure  8 with the different particle numbers. Figure 9 gives the convergence tests on displacement, in which N is the total particle number and different values using 1600, 3600, 6400 and 10,000 are analyzed here. The relative error Err is defined as the errors between FEM result and SPH, SPH_SFDI results, which are calculated by  Figure 8 that the error of force decreases as the particle number increases unanimously for both the SPH and SPH_ SFDI approaches. This indicates the convergence of all numerical models. However, the error magnitude of SPH_SFDI is much smaller than that of SPH. Besides, we could also conclude from Figure 9 that the convergence of SPH_ SFDI method is much better than that of the SPH, in that the errors of the former reduce more rapidly following the refinement of spatial resolutions. To evaluate the enhanced performance of SPH_SFDI method further, the convergence properties of the SPH and SPH_SFDI models are now examined in terms of the displacement y. For this purpose, the time histories of displacement computed by SPH and SPH_SFDI are presented in Figure 8 with the different particle numbers. Figure 9 gives the convergence tests on displacement, in which N is the total particle number and different values using 1600, 3600, 6400 and 10,000 are analyzed here. The relative error Err is defined as the errors between FEM result and SPH, SPH_SFDI results, which are calculated by Err = |y 0 − y|/y 0 , where y is the computed displacement by SPH and SPH_SFDI from t = 0.0 s to t = 2.0 s, y 0 is the displacement of FEM from t = 0.0 s to t = 2.0 s. It is shown in Figure 8 that the error of force decreases as the particle number increases unanimously for both the SPH and SPH_ SFDI approaches. This indicates the convergence of all numerical models. However, the error magnitude of SPH_SFDI is much smaller than that of SPH. Besides, we could also conclude from Figure 9 that the convergence of SPH_ SFDI method is much better than that of the SPH, in that the errors of the former reduce more rapidly following the refinement of spatial resolutions.

Four-Point Pending of Ice Beam
The ice four-point bending experiment was conducted by Kujala et al. [4]. In their work, a loading rig was used to bend the ice beam upward during the experiments. In addition, they used a hydraulic cylinder to push two moving supports to produce a force, which were located 1 m apart in

Four-Point Pending of Ice Beam
The ice four-point bending experiment was conducted by Kujala et al. [4]. In their work, a loading rig was used to bend the ice beam upward during the experiments. In addition, they used a hydraulic cylinder to push two moving supports to produce a force, which were located 1 m apart in

Four-Point Pending of Ice Beam
The ice four-point bending experiment was conducted by Kujala et al. [4]. In their work, a loading rig was used to bend the ice beam upward during the experiments. In addition, they used a hydraulic cylinder to push two moving supports to produce a force, which were located 1 m apart in the middle of the beam, so it can bend the ice beam upward. At the same time, two fixed supports, which were 4 m apart, were placed at both ends of the beam to against the ice beam. The detailed resulting measurements can be found in Ehlers and Kujala [46].
In this section, the ice beam, the upper and lower supports are modeled with the particles, which is shown in Figure 10. the middle of the beam, so it can bend the ice beam upward. At the same time, two fixed supports, which were 4 m apart, were placed at both ends of the beam to against the ice beam. The detailed resulting measurements can be found in Ehlers and Kujala [46]. In this section, the ice beam, the upper and lower supports are modeled with the particles, which is shown in Figure 10  To show the fracture patterns clearly, Figure 11 gives the results obtained by SPH_SFDI using the non-associative flow rule at different time. As shown in Figure 11a, two fracture cracks obviously occur at the upper area by two moving supports and the ice beam breaks into three sections. Then the cracks in the ice beam widen and two sections of ice beams on either side of the two bottom supports sink downward as shown in Figure 12b. As the two supports move up slowly, the cracks in the ice beam widen and the ice beam eventually breaks into three sections as shown in Figure 12c. It need to highlight the point that in our numerical results, due to the complete symmetry of the characteristics of ice beam and the external loading and supporting condition, the fracture location of the ice beam is almost completely symmetric, which is not completely consistent with that in the experimental results. In addition, there are slight crushing failures at the place contacted with the upper two fixed supports, as the two upper supports are fixed and the ice beam has a tendency to move upward.  To show the fracture patterns clearly, Figure 11 gives the results obtained by SPH_SFDI using the non-associative flow rule at different time. As shown in Figure 11a, two fracture cracks obviously occur at the upper area by two moving supports and the ice beam breaks into three sections. Then the cracks in the ice beam widen and two sections of ice beams on either side of the two bottom supports sink downward as shown in Figure 12b. As the two supports move up slowly, the cracks in the ice beam widen and the ice beam eventually breaks into three sections as shown in Figure 12c. It need to highlight the point that in our numerical results, due to the complete symmetry of the characteristics of ice beam and the external loading and supporting condition, the fracture location of the ice beam is almost completely symmetric, which is not completely consistent with that in the experimental results. In addition, there are slight crushing failures at the place contacted with the upper two fixed supports, as the two upper supports are fixed and the ice beam has a tendency to move upward. the middle of the beam, so it can bend the ice beam upward. At the same time, two fixed supports, which were 4 m apart, were placed at both ends of the beam to against the ice beam. The detailed resulting measurements can be found in Ehlers and Kujala [46]. In this section, the ice beam, the upper and lower supports are modeled with the particles, which is shown in Figure 10  To show the fracture patterns clearly, Figure 11 gives the results obtained by SPH_SFDI using the non-associative flow rule at different time. As shown in Figure 11a, two fracture cracks obviously occur at the upper area by two moving supports and the ice beam breaks into three sections. Then the cracks in the ice beam widen and two sections of ice beams on either side of the two bottom supports sink downward as shown in Figure 12b. As the two supports move up slowly, the cracks in the ice beam widen and the ice beam eventually breaks into three sections as shown in Figure 12c. It need to highlight the point that in our numerical results, due to the complete symmetry of the characteristics of ice beam and the external loading and supporting condition, the fracture location of the ice beam is almost completely symmetric, which is not completely consistent with that in the experimental results. In addition, there are slight crushing failures at the place contacted with the upper two fixed supports, as the two upper supports are fixed and the ice beam has a tendency to move upward.  The snapshots of the failure path predicted by SPH and SPH_SFDI using non-associative flow rule at t = 2.75 s are shown in Figure 7. According to the results of Figure 7, the ice beam breaks into three segments which can be obtained both by the standard SPH and SPH_SFDI. The fracture points of the horizontal coordinate on lower two moving supports by standard SPH show the apparent inward deviation compared with the ones of SPH_SFDI. Furthermore, it can be easily observed that the particle distributions for the results of standard SPH are in chaotic, whereas the results of SPH_SFDI are more stable and reliable. To show the accuracy of numerical solutions of standard SPH and SPH_SFDI, Figure 13 gives the comparison of the force time histories among stand SPH, SPH_SFDI and the experimental data [46]. The relationship between external force and flexural stress according to ITTC [47] is defined as: The snapshots of the failure path predicted by SPH and SPH_SFDI using non-associative flow rule at t = 2.75 s are shown in Figure 7. According to the results of Figure 7, the ice beam breaks into three segments which can be obtained both by the standard SPH and SPH_SFDI. The fracture points of the horizontal coordinate on lower two moving supports by standard SPH show the apparent inward deviation compared with the ones of SPH_SFDI. Furthermore, it can be easily observed that the particle distributions for the results of standard SPH are in chaotic, whereas the results of SPH_SFDI are more stable and reliable. To show the accuracy of numerical solutions of standard SPH and SPH_SFDI, Figure 13 gives the comparison of the force time histories among stand SPH, SPH_SFDI and the experimental data [46]. The relationship between external force and flexural stress according to ITTC [47] is defined as: The snapshots of the failure path predicted by SPH and SPH_SFDI using non-associative flow rule at t = 2.75 s are shown in Figure 7. According to the results of Figure 7, the ice beam breaks into three segments which can be obtained both by the standard SPH and SPH_SFDI. The fracture points of the horizontal coordinate on lower two moving supports by standard SPH show the apparent inward deviation compared with the ones of SPH_SFDI. Furthermore, it can be easily observed that the particle distributions for the results of standard SPH are in chaotic, whereas the results of SPH_SFDI are more stable and reliable.
To show the accuracy of numerical solutions of standard SPH and SPH_SFDI, Figure 13 gives the comparison of the force time histories among stand SPH, SPH_SFDI and the experimental data [46]. The relationship between external force and flexural stress according to ITTC [47] is defined as: where B is the width of 3D ice beam. We use the same value of B as that in the experiment of Kujala et al. [4]. L 0 is the distance between the fixed support and the bottom moving support on the same side. σ is flexural stress generated by bending of ice beam. According to the results in Figure 13, the numerical results of SPH_SFDI have obviously better agreement with the experimental data than the ones of standard SPH results. With the accuracy improvement of the gradient approximation, the force time histories and fractured crack of the ice beam bending can get more accurate and reliable results than the ones of standard SPH.  (48) where B is the width of 3D ice beam. We use the same value of B as that in the experiment of Kujala et al. [4]. 0 L is the distance between the fixed support and the bottom moving support on the same side. σ is flexural stress generated by bending of ice beam.
According to the results in Figure 13, the numerical results of SPH_SFDI have obviously better agreement with the experimental data than the ones of standard SPH results. With the accuracy improvement of the gradient approximation, the force time histories and fractured crack of the ice beam bending can get more accurate and reliable results than the ones of standard SPH. To validate that the numerical model can simulate the failure of ice beam at different times effectively, the failure of ice beam with the same material parameters above under two extra moving velocities of upward supports, such 1 V = 0.001842 m/s and 2 V = 0.003225 m/s are also considered. Figure 14 gives the comparison of the force versus time curve under different moving velocities of two upward supports among SPH_SFDI and the experimental data [46]. In addition, in Table 1, the results from SPH_SFDI has been compared with experiment tests in terms of the failure force F′ , failure time t′ and the corresponding deflection δ . The good agreement between the numerical results and the experiment data can be obtained clearly in Figure 14 and Table 1 To validate that the numerical model can simulate the failure of ice beam at different times effectively, the failure of ice beam with the same material parameters above under two extra moving velocities of upward supports, such V 1 = 0.001842 m/s and V 2 = 0.003225 m/s are also considered. Figure 14 gives the comparison of the force versus time curve under different moving velocities of two upward supports among SPH_SFDI and the experimental data [46]. In addition, in Table 1, the results from SPH_SFDI has been compared with experiment tests in terms of the failure force F , failure time t and the corresponding deflection δ. The good agreement between the numerical results and the experiment data can be obtained clearly in Figure 14 and Table 1, although there exists some little difference. Thus, the presented SPH_SFDI model including the cohesion softening model can imply the time dependent of ice failure and get good simulated results for ice failure with different loading rates. where B is the width of 3D ice beam. We use the same value of B as that in the experiment of Kujala et al. [4]. 0 L is the distance between the fixed support and the bottom moving support on the same side. σ is flexural stress generated by bending of ice beam.
According to the results in Figure 13, the numerical results of SPH_SFDI have obviously better agreement with the experimental data than the ones of standard SPH results. With the accuracy improvement of the gradient approximation, the force time histories and fractured crack of the ice beam bending can get more accurate and reliable results than the ones of standard SPH. To validate that the numerical model can simulate the failure of ice beam at different times effectively, the failure of ice beam with the same material parameters above under two extra moving velocities of upward supports, such 1 V = 0.001842 m/s and 2 V = 0.003225 m/s are also considered. Figure 14 gives the comparison of the force versus time curve under different moving velocities of two upward supports among SPH_SFDI and the experimental data [46]. In addition, in Table 1, the results from SPH_SFDI has been compared with experiment tests in terms of the failure force F ′ , failure time t′ and the corresponding deflection δ . The good agreement between the numerical results and the experiment data can be obtained clearly in Figure 14 and Table 1, although there exists some little difference. Thus, the presented SPH_SFDI model including the cohesion softening model can imply the time dependent of ice failure and get good simulated results for ice failure with different loading rates.  To show the effects of different plastic flow rules, Figure 15 gives the comparisons of the force time histories by SPH_SFDI with associative plastic flow (Equation (17)) and non-associative plastic flow (Equation (19)). To show the difference between associative and non-associative plastic flow clearly, Figure 16 gives the snapshots of the cracks in the brittle failure process obtained by SPH_SFDI with associative and non-associative flow rule at t = 2.75 s. The force obtained by the associative rule is basically consistent with the experimental data and the failure paths are also consistent with the non-associative flow rule. According to the results of Figure 16, the particles on the cracks, especially near the left bottom support, are slightly disordered by associative flow rule. In comparison, these particles near the same domain are more regular and reliable by non-associative flow rule.   To show the effects of different plastic flow rules, Figure 15 gives the comparisons of the force time histories by SPH_SFDI with associative plastic flow (Equation (17)) and non-associative plastic flow (Equation (19)). To show the difference between associative and non-associative plastic flow clearly, Figure 16 gives the snapshots of the cracks in the brittle failure process obtained by SPH_SFDI with associative and non-associative flow rule at t = 2.75 s. The force obtained by the associative rule is basically consistent with the experimental data and the failure paths are also consistent with the non-associative flow rule. According to the results of Figure 16, the particles on the cracks, especially near the left bottom support, are slightly disordered by associative flow rule. In comparison, these particles near the same domain are more regular and reliable by non-associative flow rule.   To show the effects of different plastic flow rules, Figure 15 gives the comparisons of the force time histories by SPH_SFDI with associative plastic flow (Equation (17)) and non-associative plastic flow (Equation (19)). To show the difference between associative and non-associative plastic flow clearly, Figure 16 gives the snapshots of the cracks in the brittle failure process obtained by SPH_SFDI with associative and non-associative flow rule at t = 2.75 s. The force obtained by the associative rule is basically consistent with the experimental data and the failure paths are also consistent with the non-associative flow rule. According to the results of Figure 16, the particles on the cracks, especially near the left bottom support, are slightly disordered by associative flow rule. In comparison, these particles near the same domain are more regular and reliable by non-associative flow rule.

Uniaxial Compressive Test of Ice Specimen
In this section, it will justify the efficiency of SPH_SFDI scheme for the ice compressed behavior simulation. Uniaxial compression of ice specimen is one of the most introduced benchmarks in this field. A two-dimensional rectangle ice specimen will be considered. The width ( D ) and height ( H ) of the ice specimen are 7 cm and 17.5 cm, respectively. The schematic geometry of this model can be shown in Figure 17. An axial velocity with vertical downward is loaded on the upper platen, which is of the value 0.0034675 m/s. The experiment of the same scale model was conducted by Li et al. [7] and Zhang [48].   Figure 18a illustrates the comparisons of the axial stress-strain curves among the experimental data [7], standard SPH and SPH_SFDI with non-associative flow rule. According to the results of Figure 18a, the stress-strain relation obtained by the SPH_SFDI method is in more agreement with experimental data than the ones of standard SPH, despite some unavoidable discrepancies due to the complication of the physical problem. Figure 18b gives the comparisons of the stress-strain curve in the experimental data with the results obtained by SPH_SFDI with associative and non-associative plastic flow rules. According to the results of Figure 18, there exists a certain difference between the numerical results and the experimental results for real sea ice. The elastic-plastic model of this paper can get a reasonable agreement with the ones of experimental data. However, the nonlinear behavior

Uniaxial Compressive Test of Ice Specimen
In this section, it will justify the efficiency of SPH_SFDI scheme for the ice compressed behavior simulation. Uniaxial compression of ice specimen is one of the most introduced benchmarks in this field. A two-dimensional rectangle ice specimen will be considered. The width (D) and height (H) of the ice specimen are 7 cm and 17.5 cm, respectively. The schematic geometry of this model can be shown in Figure 17. An axial velocity with vertical downward is loaded on the upper platen, which is of the value 0.0034675 m/s. The experiment of the same scale model was conducted by Li et al. [7] and Zhang [48]. Two rigid plates on the top and bottom deal with solid boundary, which can support the cuboid. The top plate could be moved freely in the vertical direction with a certain velocity, which focuses on the compressed ice behavior. The ice specimen has the cohesion c = 0.45 Mpa, the friction angle is 22.5 • . The dilatancy angle ϕ in the non-associative plastic rule is set to be one-third of the friction angle (ϕ = φ/3).

Uniaxial Compressive Test of Ice Specimen
In this section, it will justify the efficiency of SPH_SFDI scheme for the ice compressed behavior simulation. Uniaxial compression of ice specimen is one of the most introduced benchmarks in this field. A two-dimensional rectangle ice specimen will be considered. The width ( D ) and height ( H ) of the ice specimen are 7 cm and 17.5 cm, respectively. The schematic geometry of this model can be shown in Figure 17. An axial velocity with vertical downward is loaded on the upper platen, which is of the value 0.0034675 m/s. The experiment of the same scale model was conducted by Li et al. [7] and Zhang [48]. Two rigid plates on the top and bottom deal with solid boundary, which can support the cuboid. The top plate could be moved freely in the vertical direction with a certain velocity, which focuses on the compressed ice behavior. The ice specimen has the cohesion c = 0.45 Mpa, the friction angle is 22.5°. The dilatancy angle ϕ in the non-associative plastic rule is set to be onethird of the friction angle (   Figure 18a illustrates the comparisons of the axial stress-strain curves among the experimental data [7], standard SPH and SPH_SFDI with non-associative flow rule. According to the results of Figure 18a, the stress-strain relation obtained by the SPH_SFDI method is in more agreement with experimental data than the ones of standard SPH, despite some unavoidable discrepancies due to the complication of the physical problem. Figure 18b gives the comparisons of the stress-strain curve in the experimental data with the results obtained by SPH_SFDI with associative and non-associative plastic flow rules. According to the results of Figure 18, there exists a certain difference between the numerical results and the experimental results for real sea ice. The elastic-plastic model of this paper can get a reasonable agreement with the ones of experimental data. However, the nonlinear behavior  Figure 18a illustrates the comparisons of the axial stress-strain curves among the experimental data [7], standard SPH and SPH_SFDI with non-associative flow rule. According to the results of Figure 18a, the stress-strain relation obtained by the SPH_SFDI method is in more agreement with experimental data than the ones of standard SPH, despite some unavoidable discrepancies due to the complication of the physical problem. Figure 18b gives the comparisons of the stress-strain curve in the experimental data with the results obtained by SPH_SFDI with associative and non-associative plastic flow rules. According to the results of Figure 18, there exists a certain difference between the numerical results and the experimental results for real sea ice. The elastic-plastic model of this paper can get a reasonable agreement with the ones of experimental data. However, the nonlinear behavior of stress-strain time histories cannot be captured exactly; there exists many different factors, such the ice viscosity, the anisotropy and the temperature, which should be further investigated to make it more reliable for the numerical simulation of real sea ice. Water 2017, 9,882 19 of 24 of stress-strain time histories cannot be captured exactly; there exists many different factors, such the ice viscosity, the anisotropy and the temperature, which should be further investigated to make it more reliable for the numerical simulation of real sea ice.  Figure 19 shows the comparisons of the typical fracture pattern among experimental results of Zhang [48] (Figure 19a), the standard SPH and SPH_SFDI with non-associative flow rule. According to the results of Figure 19, the shear failures in the ice sample are predicted by the standard SPH and SPH_SFDI at t = 0.69 s and t = 1.38 s. The ice specimen exists the brittle failure and there is a main crack in the fracture pattern. The upper part of the body has the trend of sliding along the main crack and falling out of the specimen. Although the standard SPH method can predict the shear failure, the position of fracture crack differs greatly from the experimental result. The results of the SPH_SFDI are in better agreement with the experimental test than the ones of the standard SPH. In addition, some irrational damage occurred where the particle distribution is obviously ill conditioned in the SPH result, which can be seen in Figure 19b1,c1. By comparison, the results of SPH_SFDI are more stable and more regular. In summary, the present simulations also provide a strong indication that the results of SPH_SFDI method could be superior to the standard SPH in predicting the compressive failure process accurately. It should be noted that with the development of shear failure, the lower part of the ice specimen tilts under the downward sliding extrusion of the upper part and deformation occurs at the lower left corner, as shown in Figure 19c.  Figure 19 shows the comparisons of the typical fracture pattern among experimental results of Zhang [48] (Figure 19a), the standard SPH and SPH_SFDI with non-associative flow rule. According to the results of Figure 19, the shear failures in the ice sample are predicted by the standard SPH and SPH_SFDI at t = 0.69 s and t = 1.38 s. The ice specimen exists the brittle failure and there is a main crack in the fracture pattern. The upper part of the body has the trend of sliding along the main crack and falling out of the specimen. Although the standard SPH method can predict the shear failure, the position of fracture crack differs greatly from the experimental result. The results of the SPH_SFDI are in better agreement with the experimental test than the ones of the standard SPH. In addition, some irrational damage occurred where the particle distribution is obviously ill conditioned in the SPH result, which can be seen in Figure 19b1,c1. By comparison, the results of SPH_SFDI are more stable and more regular. In summary, the present simulations also provide a strong indication that the results of SPH_SFDI method could be superior to the standard SPH in predicting the compressive failure process accurately. It should be noted that with the development of shear failure, the lower part of the ice specimen tilts under the downward sliding extrusion of the upper part and deformation occurs at the lower left corner, as shown in Figure 19c.  Figure 19 shows the comparisons of the typical fracture pattern among experimental results of Zhang [48] (Figure 19a), the standard SPH and SPH_SFDI with non-associative flow rule. According to the results of Figure 19, the shear failures in the ice sample are predicted by the standard SPH and SPH_SFDI at t = 0.69 s and t = 1.38 s. The ice specimen exists the brittle failure and there is a main crack in the fracture pattern. The upper part of the body has the trend of sliding along the main crack and falling out of the specimen. Although the standard SPH method can predict the shear failure, the position of fracture crack differs greatly from the experimental result. The results of the SPH_SFDI are in better agreement with the experimental test than the ones of the standard SPH. In addition, some irrational damage occurred where the particle distribution is obviously ill conditioned in the SPH result, which can be seen in Figure 19b1,c1. By comparison, the results of SPH_SFDI are more stable and more regular. In summary, the present simulations also provide a strong indication that the results of SPH_SFDI method could be superior to the standard SPH in predicting the compressive failure process accurately. It should be noted that with the development of shear failure, the lower part of the ice specimen tilts under the downward sliding extrusion of the upper part and deformation occurs at the lower left corner, as shown in Figure 19c.   Figure 20 shows a direct comparison for the brittle shear failure simulation by SPH_SFDI with the associative flow rule and non-associative flow rule respectively. Although the stress-strain curve obtained by the associative rule is basically consistent with the experimental data (seen in Figure 18b) and the failure patterns predicted by the associative rule are also consistent with that of the nonassociative flow rule, it can be seen in Figure 20a that there is a slight particle strip distribution in the SPH_SFDI results with associative flow rule. In contrast, the distribution of particles for the SPH_SFDI result using non-associative flow rule is more stable and reliable, more details can be found in Figure 20b. Therefore, with the combination of a comparative analysis of the different flow  Figure 20 shows a direct comparison for the brittle shear failure simulation by SPH_SFDI with the associative flow rule and non-associative flow rule respectively. Although the stress-strain curve obtained by the associative rule is basically consistent with the experimental data (seen in Figure 18b) and the failure patterns predicted by the associative rule are also consistent with that of the non-associative flow rule, it can be seen in Figure 20a that there is a slight particle strip distribution in the SPH_SFDI results with associative flow rule. In contrast, the distribution of particles for the SPH_SFDI result using non-associative flow rule is more stable and reliable, more details can be found in Figure 20b. Therefore, with the combination of a comparative analysis of the different flow laws in the four points bending, it can be found that the non-associative flow rule can yield the better results for simulating the fractures.
In addition, Figure 21 shows the comparison of the bulge fracture patterns among the experiment results by Zhang [48], the standard SPH and SPH_SFDI. According to the results of Figure 21b, the ice sample exhibits the ductile failure feature and the bulge failure occurs in the bottom part of the ice sample. During the compression process, the failure progress of the ice sample is slow and there is no obvious main crack at the stage of the specimen failure. The bottom part of the ice sample distends to the two outer sides and eventually fractures. It is also shown in Figure 21b that the predicted cracks by SPH_SFDI can make a better agreement with the experimental test than the ones of the standard SPH. Although the standard SPH method can predict the bulge failure, the position of bulge fracture differs greatly from the experimental results. In addition, the particles below the damage position are obviously disordered in the results of standard SPH, which can be shown in Figure 21b1. By comparison, the results of SPH_SFDI are more reliable and the particle distributions of SPH_SFDI are more regular.

Conclusions
In this paper, the SPH_SFDI model including the elastic-plastic cohesion softening Drucker-Prager failure model is proposed to simulate the bending and compression failure processes of ice. The predicted force in a four-point bending and the axial stress of a uniaxial compressive test are in a good agreement with the experimental data. The simulated fracture patterns are also reasonably close to the reality. The conducted studies disclosed that the elasto-plastic cohesion softening Drucker-Prager failure model, which originated from the soil and rock mechanics, can also be Figure 21. Comparisons of the bulge fracture pattern among (a) the bulge fracture patterns among the experiment results by Zhang [48], (b) the standard SPH and SPH_SFDI results (The color legend is the accumulated plastic strain).

Conclusions
In this paper, the SPH_SFDI model including the elastic-plastic cohesion softening Drucker-Prager failure model is proposed to simulate the bending and compression failure processes of ice. The predicted force in a four-point bending and the axial stress of a uniaxial compressive test are in a good agreement with the experimental data. The simulated fracture patterns are also reasonably close to the reality. The conducted studies disclosed that the elasto-plastic cohesion softening Drucker-Prager failure model, which originated from the soil and rock mechanics, can also be effectively used to simulate the physical destruction phenomena during the failure process of the ice. According to the comparisons between the numerical results conducted by the standard SPH and improved SPH_SFDI, the performance of the latter is found to be much better in view of the numerical accuracy and stability in the study of the bending and compression failure processes of ice.