A MultiYieldSurface Plasticity StateBased Peridynamics Model and its Applications to Simulations of IceStructure Interactions
https://doi.org/10.1007/s11804023003448

Abstract
Due to complex mesoscopic and the distinct macroscopic evolution characteristics of ice, especially for its brittletoductile transition in dynamic response, it is still a challenging task to build an accurate ice constitutive model to predict ice loads during shipice collision. To address this, we incorporate the conventional multiyieldsurface plasticity model with the statebased peridynamics to simulate the stress and crack formation of ice under impact. Additionally, we take into account of the effects of inhomogeneous temperature distribution, strain rate, and pressure sensitivity. By doing so, we can successfully predict material failure of isotropic freshwater ice, iceberg ice, and columnar saline ice. Particularly, the proposed ice constitutive model is validated through several benchmark tests, and proved its applicability to model ice fragmentation under impacts, including drop tower tests and ballistic problems. Our results show that the proposed approach provides good computational performance to simulate shipice collision.
Keywords:
 Ice loads ·
 Multiyieldsurface ·
 Constitutive model ·
 Peridynamics ·
 Icestructure interaction
Article Highlights• Compared with the traditional ice constitutive model, we can successfully predict material failure of isotropic freshwater ice, iceberg ice, and columnar saline ice.• The proposed ice constitutive model can capture the ice failure process by taking into account of the effects of inhomogeneous temperature distribution, strain rate, and pressure sensitivity.• We successfully applied the proposed peridynamics model simulating structural failure of ice in engineering settings. 
1 Introduction
Ice is formed under extreme polar weather conditions is the major threat to the arctic marine transportation as well as polar scientific research operations. Shipice collision affects the normal operation of navigation (Liu et al., 2017), and even causes the ship to sink or capsize completely, seriously endangering personnel safety. The prediction of ice load is critically important in marine structure design in cold regions (Palmer et al., 2009).
However, using the simple uniaxial strength of ice to calculate ice loads is inaccurate due to the complex deformation before ice causes marine structures to fail and rupture (Zhang et al., 2019). To study the ice load on naval architectural structures (Zhang et al., 2017; Jia et al., 2019), most of the empirical design formulas or semiempirical formulas that sum up from a large number of tension or compression experimental test (Cui et al., 2018) data are used to construct ice constitutive model equations (Snyder et al., 2016). Since they are based on many empirical assumptions, this approach has flaws in calculating the strength of ice (Mellor and Cole, 1982). Current constitutive models (Sain and Narasimhan, 2011) are difficult to provide a unified multiple scale modeling, and it is difficult to simulate the mechanical behavior of transition regions at high strain rates. Many numerical methods have been developed to predict ice failure processes (shear, melting, splitting, and fracture), showing higher efficiency than the experimental tests (Schulson, 2001). Since ice is a multiscale material with complicated physical properties, it is a challenging task to build an appropriate single scale constitutive model to facilitate accurate prediction on the ice fragmentation.
PernasSanchez et al. (2015) used the Lagrangian finite element method, arbitrary Lagrangian Eulerian (ALE) method (Schulson and Buck, 1995) and Smoothed Particle Hydrodynamics (SPH) method (Sun et al., 2021; 2019; 2018) to model and simulate mechanical behavior of ice under high impacts (Kolsky, 1949), and they compared the advantages and disadvantages of the three modeling methods. Their results show that the SPH method is the most efficient method among the three methods, enabling the ice microcracks to reach a fluidlike state before creeping (Schulson and Buck, 1995; Carney et al., 2006). In addition to SPH method (Wang et al., 2006), the Discrete Element Method (DEM) also shows good performance in simulating ice loading and ice breaking properties from a macroscopic perspective, by taking into account cohesion and fragmentation effects (Gold et al., 1988; Jones, 1982; Jones, 1997). The above work, however, are within the framework of continuum mechanics, encounter great difficulties in dealing with mesoscale discontinuities, i.e. micro cracks, which are extensively occurring in the natural ice. In nature, ice is a composite material, and it is mixture of salt, air, and water (Liu et al., 2021). There are irregular bubbles (Zhang et al., 2023) in the ice interior. The existence of bubbles in the ice interior is a universal phenomenon, which affects the physical and mechanical properties of ice. Moreover, the properties of its microstructure affect its macroscopic mechanical properties (Xie and Li, 2021a; 2021b; Xie et al., 2022a; 2022b) which are also affected by temperature, salinity, porosity, and density, making it extra challenge to study. Although many efforts have been made to simulate the mechanical behavior of ice (Vazic et al., 2020), a suitable constitutive relation of ice has not been established for engineering (Han et al., 2022) applications in wide scenarios. Existing numerical methods (Li et al., 2019) simplify ice, as an isotropic uniform material, e. g. Wang et al. (2018), which ignore many key factors affecting the failure process of ice (Lu et al., 2020). More efforts are still needed to build a realistic macroscale ice constitutive model.
Recently, peridynamics method has been used to model ice destruction due to its advantages stemming from its intrinsic nonlocal formulation in dealing with discontinuities. The integral formulation of the nonlocal peridynamic theory (Silling, 2000) avoids the mathematical difficulties in evaluation of spatial derivation of the discontinuous displacement field (Madenci and Oterkus, 2014), so that it can more accurately simulate crack growth of the ice sheet in its damage state. Moreover, the statebased peridynamics can incorporate material constitutive models into mesoscale peridynamics formulation through the "force state" in the equations of motion (Silling, 2007), so that it provides a useful platform for building multiscale material constitutive models (Fan and Li, 2017a; 2017b; Fan et al., 2015; 2016).
In this work, we adopt the constitutive modeling theory of the statebased peridynamics (Hu et al., 2020), and incorporated the material constitutive model (Drucker and Prager, 1952; Johnson and Holmquist, 1994) into the integral equation form of the "force state", so as to establish a suitable and accurate mechanical model to capture ice mechanical properties over a range of multiscale. Different from the previous work, suitable ice constitutive model has been implemented into the peridynamic state equation, which can capture the multiscale mechanical responses of ice. For the microscale, peridynamics has advantages in capturing microscopic cracks in brittle materials, and it has been successfully applied in the simulation of many brittle materials such as glass (Lai et al., 2018). For the macroscale, the constitutive model used is the modified multisurface criterion model proposed by DerradjiAouat, the brittle failure of ice at high strain rates can be well reproduced and has many successful applications, such as shipiceberg collisions, shipice interaction, and ice impact problems. The DerradjiAouat ice constitutive model have been used to simulate the shipiceberg collisions successfully, the proposed plastic material model yields good simulation results, both of the mesh convergence and the computational time is stable. By comparing the different inelastic constitutive models, elasticplastic model, elasticbrittle model and the DerradjiAouat multisurface cohesive model, the advantages and disadvantages between the three different ice constitutive models was analyzed by using the MSC.DYTRAN, and the results indicated that the multisurface failure criterion is the most suitable one in simulation of shipice interaction owning to its consideration of the effects of strain rate, temperature, and hydrostatic pressure. As one of the highlights in the present peridynamic model, the influence of the added mass by water is ignored when simulating the shipice interaction, and the ship motion is not considered. Such issues required to be taken into account in the future work.
2 Ice constitutive model
In this section, we first briefly introduce how to incorporate ice constitutive model which is called multisurface theory developed by DerradjiAouat into the statebased peridynamics formulation. Then, we discuss the cohesive law used in corresponding peridynamics to determine the failure criterion for ice. Finally, a peridynamic numerical algorithm considering strain rate effect and the pressure sensitivity for simulating fresh water isotropic ice, iceberg ice and columnar saline ice is presented.
2.1 Multiyield surface constitutive model for ice
Ice undergoes nonlinear elastoplastic behavior at low strain rate (< 10^{−3} s^{−1}), and linear elastobrittle failure at high strain rate (>10^{−3} s^{−1}). There is a transition zone at the strain rate (≈10^{−3} s^{−1}), the failure mode of ice is from ductile to brittle. Thus, it is an essential task to establish a suitable ice constitutive model to show the complex physical evolution during the interactions of ice with structures. To address this issue, DerradjiAouat (2003) modified the traditional yield surface criteria by using the multisurface theory, considerating the effects of the strain rate and the temperature of ice. It should be mentioned that this "multisurface constitutive model" has the ability to deal with not only the isotropic fresh water ice and iceberg, but also the anisotropic sea ice, which is columnar structure, as is shown in Figure 1.
To well address this, DerradjiAouat developed a failure criterion based on the previous data for three types of ice as is shown in Table 1, developed an elliptical yield envelop is given as the following formation:
$$ \left(\frac{q\eta}{q_{\max }}\right)^2+\left(\frac{P\lambda}{P_{\max }}\right)^2=1 $$ (1) Types η λ P_{max} Freshwater ice 0 45.0 55.0 Iceberg ice 0 45.0 55.0 Grown ice 0 45.0 45.0 where η, λ are the coordinates of the center point of the ellipse, P is the hydrostatic pressure, and P_{max} is the major axes of the ellipse, respectively. q represents the octahedral shear stress, q_{max} is the apex of the ellipse, DerradjiAouat (2003) provided the q_{max} by the following equations:
$$ q_{\max }=\left[\frac{\dot{\varepsilon}}{\xi}\right]^{1 / n} $$ (2) where $ n=4, \xi=5 \times 10^{6} \exp \left[10.5 \times 10^3\left(\frac{1}{T}\frac{1}{273}\right)\right]$.
In which shows the relationship between the absolute maximum octahedral shear stress q_{max}, the strain rate ε, and the temperature T. Eq. (2) combines the effects of strain rate and the temperature on the q_{max} of ice, which shows that the value of q_{max} increases with the increasing of the strain rate ε and the decreasing of the temperature; and decreases with the decreasing of the strain rate ε and the increasing of the temperature. In above equation,
$$ q=\sqrt{\frac{1}{2}\left[\left(\sigma_{11}\sigma_{22}\right)^2+\left(\sigma_{22}\sigma_{33}\right)^2+\left(\sigma_{33}\sigma_{11}\right)^2\right]} $$ (3) $$ P=\frac{1}{3}\left(\sigma_{11}+\sigma_{22}+\sigma_{33}\right) $$ (4) Figure 3 shows the shape of the failure envelope proposed by DerradjiAouat, and has been applied to model the ice materials for many years. The "TsaiWu" yield surface on the condition of η =0 for the isotropic material is usually written as the following form:
$$ f\left(p, J_2\right)=J_2\left(a_0+a_1 P+a_2 P^2\right)=0 $$ (5) In which, J_{2} is the second invariant for the deviatoric stress tensor, P is the Hydrostatic pressure. a_{0}, a_{1}, a_{2} are constants derived from the uniaxial tensile test data. The invariant of the secondorder deviatoric stress tensor can be expressed in the following form:
$$ J_2=\frac{1}{2} \boldsymbol{s}: \boldsymbol{s}, P=\sqrt{\frac{1}{3} \boldsymbol{s}: \boldsymbol{s}} $$ (6) To describe the temperature dependent mechanical response, the corresponding change of the yield function for the ice with temperaturedependence expressed can be expressed in the following expression,
$$ f(p, q)=J_2\left(a_0(T)+a_1(T) P+a_2(T) P^2\right)=0 $$ (7) in which, T is the temperature, and linear interpolation is used to describe the temperature change from the surface of the ship and core area.
Remark 2.1 In this work, we have made several improvements over the previouslyused ice constitutive model. For example, the relationship between the yield criterion and the hydrostatic pressure P is established in Eq. (5), and the relationship between the yield criterion and the material temperature T is established in Eq. (7). By doing so, the pressure sensitivity and temperature characteristics of ice are considered, which will affect the yield limit of the material.
2.2 Constitutive model in nonordinary state based peridynamics
The multiyield surface criterion for ice was successfully implemented in LSDYNA software but not to its full advantage because it is difficult to deal with the discontinue problem in traditional continuum mechanics. Whereas, Peridynamics is a nonlocal mesh free method, especially for the statebased peridynamics providing a platform to integrate material constitutive model.
In the statebased peridynamics, a continuum domain Ω_{0} is discretized by a set of material particle x_{i} with associated mass density ρ_{xi} and volume V_{i}, where i=1, 2, 3, …, ∞ is the particle index. Assume that a material particle is only influenced by forces from its neighboring particles x_{j}, j=1, 2, 3, …, ∞ within a local region, called as a horizon H_{x(i)}. The horizon is typically chosen to be a circle (in 2D) or a sphere (in 3D), with the radius δ called as the horizon size of particle x_{i}. The relative position vector pointing from particle x_{i} to x_{j} in the reference configuration is called "bond", which can be considered as a spring in the case of elastic interaction (as shown in Figure 3), and can be denoted as:
$$ \xi_{i j}=x_jx_i $$ (8) Under certain motion or deformation χ, the continuum body deforms and the relative displacement vector in current configuration is denoted by
$$ \eta_{i j}=\boldsymbol{u}\left(x_j, t\right)\boldsymbol{u}\left(x_i, t\right) $$ (9) Then, the current relative position between these two particles can be represented by the deformation state function $ \boldsymbol{Y}\left(\xi_{i j}\right)$, which maps the undeformed bond to a deformed bond:
$$ \boldsymbol{Y}\left(\xi_{i j}\right)=y_jy_i=\left(x_jx_i\right)+\left(\boldsymbol{u}_j\boldsymbol{u}_i\right)=\xi_{i j}+\eta_{i j} $$ (10) The peridynamic governing equation of motion, for particle x_{i} at time t at a reference configuration are then can be expressed as:
$$ \rho_0 \ddot{\boldsymbol{u}}\left(x_i, t\right)=\int_V \boldsymbol{f}\left(\eta_{i j}, \xi_{i j}\right) \mathrm{d} V_{x_j}+\boldsymbol{b}\left(x_i, t\right) $$ (11) where ρ_{0} is the mass density of the solid in the reference configuration, and $\boldsymbol{f}\left(\eta_{i j}, \xi_{i j}\right) $ is a pairwise function defined as the nonlocal integration of force vector that particle x_{j} exerts on the particle x_{i} as shown in Figure 3, and $\boldsymbol{b}\left(x_i, t\right) $ is the external body force density vector, as shown in Eq. (11), such particle based nonlocal formulation gives the peridynamics to be effective to represent discontinuity everywhere.
In the Peridynamics theory, besides the deformation state $\underline{\boldsymbol{Y}} $, another concern is the force state $\underline{\boldsymbol{T}} $, with which the balance of the momentum is expressed as:
$$ \begin{aligned} \rho_0 \ddot{\boldsymbol{u}}\left(x_i, t\right)= & \int_{H_{x_i}}\left[\underline{\boldsymbol{T}}\left(x_i, t\right)\left\langle x_jx_i\right\rangle\underline{\boldsymbol{T}}\left(x_j, t\right)\left\langle x_ix_j\right\rangle\right] \\ & \mathrm{d} V_{x_j}+\boldsymbol{b}\left(x_i, t\right) \end{aligned} $$ (12) In this work, the statebased Peridynamics theory takes strainrate effects into account. The governing equation can be expressed as follows:
$$ \begin{aligned} \rho_0 \ddot{\boldsymbol{u}}\left(x_i, t\right)= & \int_{H_{x_i}}\left[\underline{\boldsymbol{T}}\left(\xi_{i j}, \underline{\boldsymbol{Y}}\left(\xi_{i j}\right)\right)\underline{\boldsymbol{T}}\left(\xi_{j i}, \underline{\boldsymbol{Y}}\left(\xi_{j i}\right)\right)\right] \\ & \mathrm{d} V_{x_j}+\boldsymbol{b}\left(x_i, t\right) \end{aligned} $$ (13) We can obtain the relation between the forcevector state T and the first PiolaKirchhoff stress tensor P_{xi} of classical continuum mechanics as:
$$ \underline{\boldsymbol{T}}\left(\xi_{i j}, \underline{\boldsymbol{Y}}\left(\xi_{i j}\right)\right)=\omega\left(\xi_{i j}\right) \boldsymbol{P}_{x_i} \cdot \boldsymbol{K}_{x_i}^{1} \cdot\left\langle\xi_{i j}\right\rangle $$ (14) where ω(ξ_{ij}) is a positive scalar influence function defined on the bond length of each material bond ξ_{ij}, which is the distance of the two particles, and K_{xi} is the reference symmetric shape tensor of particle x_{i} which can be defined as:
$$ \begin{aligned} \boldsymbol{K}_{x_i} & =X(\xi) \times X(\xi) \int_{H_{x_i}} \omega(\xi) \xi \otimes \xi \mathrm{d} V_{x_j} \\ & \approx \sum\limits_{j=1}^N \omega(\xi) \xi \otimes \xi \Delta V_{x_j} \end{aligned} $$ (15) If the horizon size of the particle x_{i} is the same as that of the particle x_{j}, assume that the deformations of all the bonds ξ_{ij} within a horizon H_{xi} is uniform, we have $ \boldsymbol{K}_{x_i}=\boldsymbol{K}_{x_j}$. Based on continuum mechanics, the relation between P_{xi} and the Cauchy stress σ_{xi} is written by
$$ \boldsymbol{P}_{x_{\mathrm{i}}}=\operatorname{det}\left[\boldsymbol{F}_{x_i}\right] \sigma_{x_i} \boldsymbol{F}_{x_i}^{\mathrm{T}} $$ (16) where F_{xi} is the nonlocal deformation gradient at the particle x_{i}.
In order to obtain the first PiolaKirchhoff stress P_{xi}, we should first obtain the deformation gradient tensor F_{xi} and its time derivative $\dot{\boldsymbol{F}}_{x_i} $, the velocity tensor, $ \boldsymbol{L}=\dot{\boldsymbol{F}} \boldsymbol{F}^{1}$, and the rate of deformation tensor $\boldsymbol{D}=\frac{1}{2}\left(\boldsymbol{L}+\boldsymbol{L}^{\mathrm{T}}\right) $. Then based the hypoelasticplastic formulation that will be discussed in the next Section, we can find the Cauchy stress σ.
Here, we first discuss the construction of the nonlocal deformation gradient. The deformation gradient in a material point from the initial configuration to the current configuration is described by the deformation gradient. For quasiuniform particle distributions, the corresponding deformed bonds can be expressed as:
$$ \xi_{i j}=\underline{\boldsymbol{Y}}\left(\xi_{i j}\right)=\boldsymbol{F}_{x_i} \cdot \xi_{i j} $$ (17) where F_{xi} is a secondorder tensor and it can be viewed as the approximated deformation gradient at the particle x_{i}. Then we can obtain,
$$ \boldsymbol{N}_{x_i}=\boldsymbol{F}_{x_i}\left[\int_{H_{x_i}} \omega\left(\xi_{i j}\right) \xi_{i j} \otimes \xi_{i j} \mathrm{~d} V_{x_j}\right]=\boldsymbol{F}_{x_i} \cdot \boldsymbol{K}_{x_i} $$ (18) Then the nonlocal deformation gradient F_{xi} defined at the particle x_{i} in terms of shape tensors can be obtained as:
$$ \boldsymbol{F}_{x_i}=\boldsymbol{N}_{x_i} \cdot \boldsymbol{K}_{x_i}^{1} $$ (19) The spatial velocity gradient satisfies the following relationship:
$$ \boldsymbol{L}=\dot{\boldsymbol{F}} \boldsymbol{F}^{1} $$ (20) in which, $ \dot{\boldsymbol{F}}$ is the time derivative of the deformation gradient, defined as:
$$ \begin{aligned} \dot{\boldsymbol{F}}_{x_i} & =\int_{H_{x_i}} \omega\left(\xi_{i j}\right) \underline{\boldsymbol{Y}}\left(\xi_{i j}\right) \otimes \xi_{i j} \mathrm{~d} V_{x_j} \boldsymbol{K}_i^{1} \\ & \approx \sum\limits_{j=1}^N\left[\omega\left(\xi_{i j}\right) \underline{\boldsymbol{Y}}\left(\xi_{i j}\right) \otimes \xi_{i j} \Delta V_{x_j}\right] \boldsymbol{K}_i^{1} \end{aligned} $$ (21) The spatial velocity gradient L can be decomposed into a symmetric deformation rate tensor D and an asymmetric rotation tensor W:
$$ \boldsymbol{D}=\frac{1}{2}\left[\boldsymbol{L}+\boldsymbol{L}^{1}\right] \quad \text { and } \quad \boldsymbol{W}=\frac{1}{2}\left[\boldsymbol{L}\boldsymbol{L}^{1}\right] $$ (22) The material frameindifference criterion imposes the objectivity condition, and the constitutive model must be invariant within a certain range. Therefore, only objective quantities can be used in the constitutive model.
In this work, we adopt the HughesWinget formula (Hughes and Winget, 1980; Rubinstein and Atluri, 1983) in the stress update under large deformation condition. First, the rate of deformation tensor d is calculated as,
$$ \boldsymbol{d}=\boldsymbol{R}_t^{\mathrm{T}} \boldsymbol{D} \boldsymbol{R}_t $$ (23) where R_{t} is an orthogonal rotation tensor, which is used to describe the rotation of the rigid body at time t, and it can be expressed by the following incremental formula:
$$ \boldsymbol{R}_t=\left[\boldsymbol{I}+\frac{\sin (\Delta t \varOmega)}{\varOmega}\frac{\cos (\Delta t \varOmega)}{\varOmega^2} \varOmega^2\right] \boldsymbol{R}_{t\Delta t} $$ (24) where $ \varOmega^2=\omega_i \omega_j \text { and } \varOmega_{i j}=\mathrm{e}_{i k j} \omega_k$ is the tensor sequence number, e_{ikj} is the axial vector can be expressed as:
$$ \omega=\boldsymbol{w}+\left[\boldsymbol{I} \operatorname{tr}\left(\boldsymbol{V}_t\right)\boldsymbol{V}_t\right]^{1} \boldsymbol{z} $$ (25) in which, the vorticity $\boldsymbol{w}=\frac{1}{2} e_{i j k} \boldsymbol{W}_{j k} \boldsymbol{e}_i, \boldsymbol{z}=e_{i k j} D_{j m} V_{m k} \boldsymbol{e}_i $, the left stretch tensor can be expressed as:
$$ \boldsymbol{V}_t=\boldsymbol{V}_{t\Delta t}+\Delta \dot{t} \dot{\boldsymbol{V}}_{\Delta t} $$ (26) where $ \dot{\boldsymbol{V}}_{\Delta t}=\boldsymbol{L} \dot{\boldsymbol{V}}_t\dot{\boldsymbol{V}}_t \varOmega$ is the left stretch rate tensor at each time step, which can be used to V and R.
After the rate of deformation tensor d is being calculated, the Cauchy stress, vonMises stress, and plastic strain can be calculated subsequently. Therefore, the Cauchy stress of the material can be obtained from the following formula:
$$ \sigma=\boldsymbol{R}_t \tau \boldsymbol{R}_t^{\mathrm{T}} $$ (27) Then we can calculate the first PiolaKirchhoff stress tensor as follows,
$$ \boldsymbol{P}_{x_i}=(J \sigma) \boldsymbol{F}_{x_i}^{\mathrm{T}} $$ (28) where $ J=\operatorname{det} \boldsymbol{F}_{x_i}$.
Remark 2.2 In this work, we incorporated Eq. (28) into the statebased peridynamics to find the bond force states, after we find the Cauchy stress that is determined by the ice constitutive model. Thus, the ice constitutive model can be incorporated through Eq. (28) into the framework of the statebased peridynamics.
2.3 Plastic response for ice in the framework of peridynamics
According to the plastic flow law, under ideal conditions, the surface of the plastic region and the yield surface coincide. Therefore, the plastic potential function based on the DerradjiAouat (DA) constitutive model can be expressed as follows,
$$ \varPhi=\\boldsymbol{s}\\left(a_0+a_1 p+a_2 p^2\right)=0 $$ (29) where $ \\boldsymbol{s}\$ denotes the scalar of deviator stress tensor; a_{0} a_{1}, a_{2} are the experimental constant.
The hydrostatic pressure can be expressed as:
$$ p=\frac{1}{3} \sigma: \boldsymbol{I} $$ (30) where I is the unit second order tensor. The deviatoric stress tensor can be obtained by the following formula,
$$ \boldsymbol{s}=\sigma+p \boldsymbol{I} $$ (31) in which, σ is the stress tensor.
In order to describe the material deformation. The deformation rate tensor is decomposed into the elastic part, the plastic part and the temperature part, respectively, and can be denoted as:
$$ \boldsymbol{d}=\boldsymbol{d}^e+\boldsymbol{d}^p+\boldsymbol{d}^T $$ (32) The relationship between the elastic deformation rate tensor and elastic stress in the above formula satisfies:
$$ \sigma^{\nabla}=\boldsymbol{C}^e: \boldsymbol{d}^e $$ (33) in which, $ \sigma^{\nabla}$ is the true objective stress rate tensor, C^{e} is the elasticity coefficient tensor. In order to facilitate the subsequent derivation, the elastic coefficient tensor can be expressed as:
$$ \boldsymbol{C}^e=2 G \delta_{i k} \delta_{j l}+\left(K\frac{2}{3} G\right) \delta_{i j} \delta_{k l} $$ (34) where K is the yield modulus of the material, G is the shear modulus of the material. Generally, the value of the true equivalent stress can be obtained from the following formula:
$$ \bar{\sigma}=\sqrt{\frac{3}{2} \boldsymbol{s}: \boldsymbol{s}} $$ (35) The true plastic strain rate can be obtained from the following equation:
$$ \dot{\bar{\varepsilon}}^p=\sqrt{\frac{2}{3} \boldsymbol{d}^p: \boldsymbol{d}^p} $$ (36) Substituting Eq. (35) into Eq. (31):
$$ \sigma=\frac{2}{3} \bar{\sigma} \cdot \boldsymbol{n}p \cdot \boldsymbol{I} $$ (37) where n is the direction vector of deviator stress tensor, and I is the secondorder unit tensor:
$$ \boldsymbol{n}=\frac{3}{2} \cdot \frac{\boldsymbol{s}}{\bar{\sigma}} $$ (38) $$ \boldsymbol{I}=\delta_{i j} \boldsymbol{e}_i \otimes \boldsymbol{e}_j $$ (39) Therefore, according to the plastic flow law, the plastic part of the deformation rate d^{p} in Eq. (32) can be expressed as the product of the plastic consistency parameter term γ and the derivative of the plastic potential function $ \frac{\partial \varPhi}{\partial \sigma}$ as follows:
$$ \begin{aligned} \boldsymbol{d}^p & =\lambda \frac{\partial \Phi}{\partial \sigma}=\gamma\left(\frac{\partial \Phi}{\partial s} \boldsymbol{n}\frac{1}{3} \frac{\partial \Phi}{\partial p} \boldsymbol{I}\right) \\ & =\gamma\left(\frac{\partial s}{\partial \sigma}a_1 \frac{\partial p}{\partial \sigma}2 a_2 p \frac{\partial p}{\partial \sigma}\right) \\ & =\gamma\left(\boldsymbol{n}\frac{1}{3}\left(a_1+2 a_2 p\right) \boldsymbol{I}\right) \end{aligned} $$ (40) where $ s=\\boldsymbol{s}\$.
Under the condition that the material is deformed by rotation, the regression mapping format is used for processing. Figure 2 shows a schematic diagram of the regression mapping method. In the process of spatial integration, the integral increment of the constitutive model can be regarded as a kind of "straindriven" process, and all the strain tensor increments at the orthogonal point will update the stress at time t. The NewtonRaphson algorithm is used to solve the nonlinear equations. Because a fully implicit solver consumes computational time, we used a semiimplicit integration format to solve the integral equations. In computations, at time step n+1, the envelope of inelastic flow can be estimated by the yield function and the plastic potential function, and the corresponding updated stress is given as:
$$ \sigma_{n+1}=\sigma_{n+1}^{\mathrm{trial}}+\Delta \sigma^{\mathrm{re}} $$ (41) In the formula, the test stress can be calculated by the following formula:
$$ \sigma_{n+1}^{\text {trial }}=\sigma_n+\boldsymbol{C}^e: \Delta \boldsymbol{d} $$ (42) where σ_{n} is the stress at time step n, Δd is the total deformation increment.
The return increment of stress Δσ^{re} can be expressed as:
$$ \Delta \sigma^{\mathrm{re}}=\boldsymbol{C}^e: \Delta \boldsymbol{d}^p=\boldsymbol{C}^e: \Delta \gamma\left[\boldsymbol{n}\frac{1}{3}\left(a_1+2 a_2 p\right) \boldsymbol{I}\right] $$ (43) $$ =\Delta \gamma\left[\boldsymbol{K}\left(a_1+2 a_2 p\right) \boldsymbol{I}2 G \boldsymbol{n}\right] $$ (44) where Δγ is the plastic consistency parameter increment. Therefore, the test stress is updated to the following expression:
$$ \sigma_{n+1}^{\text {trial }}=\sigma_n+\Delta t\left(\boldsymbol{C}^e: \boldsymbol{d}^e\right)=\sigma_n+\Delta t\left[\boldsymbol{C}^e:\left(\boldsymbol{d}\boldsymbol{d}^p\right)\right] $$ (45) $$ =\left[\sigma_n+\Delta t\left(\boldsymbol{C}^e: \boldsymbol{d}\right)\right]\Delta t\left(\boldsymbol{C}^e: \boldsymbol{d}^p\right) $$ (46) $$ =\sigma^e+\Delta t\left[K \gamma_{n+1}\left(a_1+2 a_2 p\right) \boldsymbol{I}2 G \gamma_{n+1} \boldsymbol{n}\right] $$ (47) The stress is updated to the following expression:
$$ \sigma_{n+1}=\sigma_{n+1}^{\text {trial }}+\Delta \sigma^{\mathrm{re}} $$ (48) $$ =\sigma_{n+1}^{\text {trial }}+\Delta \gamma\left[\boldsymbol{K}\left(a_1+2 a_2 p\right) \boldsymbol{I}2 G \boldsymbol{n}\right] $$ (49) In the process of the iterative solution, the deformation of the material first undergoes an elastic phase, and then plastic deformation occurs after the elastic phase. The elastic and plastic phases are judged by the yield equation. The material reaches the yield condition at the Gauss integration point. The yield equation at the n+1 time step is written as follows:
$$ f_{n+1}^{\text {trial }}=\left\\boldsymbol{s}^{\text {tr }}\right\\left[a_0+a_1\left(p_{n+1}^{\text {trial }}\right)+a_2\left(p_{n+1}^{\text {trial }}\right)^2\right] $$ (50) If $ f_{n+1}^{\text {trial }}<0$, the material enters the elastic phase; if $f_{n+1}^{\text {trial }} \geqslant 0 $, the material enters the plastic phase. As a result, the expression of the plastic consistency parameter increment Δγ used to solve the return value increment of the stress in Eq. (43) can be obtained:
$$ \Delta \gamma=\frac{f_{n+1}^{\text {trial }}}{2 G+K\left(a_1+2 a_2 p_{n+1}^{\text {trial }}\right)} $$ (51) The return increment value of the true equivalent stress can be solved by the increment of the plastic consistency parameter as:
$$ \Delta \bar{\sigma}^{\mathrm{re}}=\Delta \gamma G $$ (52) The return increment value of the pressure is solved by the plastic consistency parameter increment, which can be expressed as follows:
$$ \Delta p^{\mathrm{re}}=\frac{\Delta \sigma^{\mathrm{re}}: \boldsymbol{I}}{3}=\Delta \gamma K\left(a_1+2 a_2 p_{n+1}^{\text {trial }}\right) $$ (53) From Eq. (48) to Eq. (49), the expression of the envelope$\{p, \bar{\sigma}\} $of the twodimensional space can be obtained:
$$ K^{\mathrm{re}}=\frac{\Delta \bar{\sigma}^{\mathrm{re}}}{\Delta p^{\mathrm{re}}}=\frac{G}{K\left(a_1+2 a_2 p_{n+1}^{\text {trial }}\right)} $$ (54) Finally, the updated equivalent plastic strain is obtained by using the plastic consistency parameter increment to solve the expression:
$$ \begin{aligned} \dot{\bar{\varepsilon}}_{n+1}^p & =\sqrt{\frac{2}{3} \boldsymbol{d}_{n+1}^p \cdot \boldsymbol{d}_{n+1}^p}=\gamma_{n+1} \sqrt{\boldsymbol{I}+\left(a_1+2 a_2 p_{n+1}^{\text {trial }}\right)^2} \\ & =\frac{\Delta \gamma}{\Delta t} \sqrt{\boldsymbol{I}+\left(a_1+2 a_2 p_{n+1}^{\text {trial }}\right)^2} \end{aligned} $$ (55) The hydrostatic pressure can be expressed as follows:
$$ p_t=\frac{1}{3} \operatorname{tr}\left(\tau_t^{\mathrm{tr}}\right)+\Delta p_t $$ (56) Therefore, the unrotated Cauchy stress including deviatoric stress and hydrostatic pressure is expressed by the formula in Eq. (28). The unrotated Cauchy stress needs to be the rotating Cauchy stress. Finally, the unrotated Cauchy stress tensor of the material can be calculated by the following equation:
$$ \tau_t=\boldsymbol{S}+\boldsymbol{I} p_t $$ (57) Then the real VonMises stress at time t is updated as:
$$ \sigma_t=\boldsymbol{R}_t \tau_t \boldsymbol{R}_t^{\mathrm{T}} $$ (58) Remark 2.3 One of the main features of the proposed peridynamics model is that the updated multiyield surface criterion shown in Eq. (50) is used to calculate inelastic stress of the ice. Ice will be under plastic deformation when $f_{n+1}^{\text {trial }} \geqslant 0 $, as is shown in Figure 3, and the plastic response can be calculated by Eq. (58).
2.4 Damage model for ice
In peridynamics, the key component used to describe the break bond is the "bond stretch" as follows:
$$ s=\frac{\left\xi_{i j}+\eta_{i j}\right\left\xi_{i j}\right}{\left\xi_{i j}\right} $$ (59) The "bond stretch" is the bridge to relate peridynamics to the fracture toughness G_{c}, which is a key classical fracture parameters can be measured for ice. We define the work requred to break the bond across the fracture surface is G_{c}:
$$ G_c=\int_0^\delta \int_0^{2 \pi} \int_z^\delta \int_0^{\cos ^{1} z / \xi}\left(\frac{c s_0^2 \xi}{2}\right) \xi^2 \sin \phi \mathrm{d} \phi \mathrm{d} \xi \mathrm{d} z $$ (60) The critical "bond stretch" is given as,
$$ s_0=\sqrt{\frac{5 \pi G_0}{18 E \delta}} $$ (61) G_{0} is used to determine by the fracture toughness K_{I}:
$$ G_0=\frac{K_I^2}{E} $$ (62) From Figure 4, we can get s and ξ as follows:
$$ \overline{\mathrm{AC}}=\sqrt{\xi_n^2+\varDelta_n^2+2 \xi_n \varDelta \cos \phi} $$ (63) $$ \xi_n=\frac{\varDelta_n\left((\cos \phi)+\sqrt{\cos ^2 \phi+2 s_n+s_n^2}\right)}{s_n\left(2+s_n\right)} $$ (64) $$ s_n=\sqrt{1+2\left(\frac{\varDelta_n}{\xi_n}\right) \cos \phi+\left(\frac{\varDelta_n}{\xi_n}\right)^2}1 $$ (65) From Figure 5, we can get s and ξ as follows:
$$ \overline{\mathrm{AC}}=\sqrt{\xi_t^{\prime 2}+\varDelta_t^2+2 \xi_t \varDelta_t \sin \phi \sin \theta} $$ (66) $$ \xi_t=\frac{\varDelta_t\left((\sin \phi \sin \theta)+\sqrt{\sin ^2 \phi \sin ^2 \theta+2 s_t+s_t^2}\right)}{s_t\left(2+s_t\right)} $$ (67) $$ S_t=\sqrt{1+2\left(\frac{\varDelta_t}{\xi_t}\right) \sin \phi \sin \theta+\left(\frac{\varDelta_t}{\xi_t}\right)^2}1 $$ (68) The damage will occur when the length exceed ξ, which is depending on s and Δ.
Remark 2.4 We note that in the proposed peridynamics model, both normal and tangential deformation as well as fracture are considered to predict the quasibrittle failure of ice by using the above established equations.
2.5 Constitutive update
The computational algorithm of the constitutive update or the stress calculation for the ice constitutive model is given as follows:
1. Start to search all material points and initialization parameters;
2. Calculate shape tensor K and inverse;
3. Calculate the approximate deformation gradient $\dot{\boldsymbol{F}} $ and its time derivative;
4. Calculate the spatial velocity gradient L, the deformation rate tensor D, and the asymmetric rotation tensor W;
5. Update the left stretch tensor V and rotate tensor R;
6. Calculate the nonrotational deformation rate tensor d;
7. Select the material constitutive model, and calculate and update the nonrotational Cauchy stress, the vonMises stress and plastic strain;
(a) Calculate the nonrotational Cauchy stress $\sigma_{n+1}^{\text {trial }} $ of the elastic test;
(b) Calculate the yield function, $ f_{n+1}^{\text {trial }}$, and execute the judgment statement: If $f_{n+1}^{\text {trial }} < 0 $, enter the elastic deformation calculation; If $f_{n+1}^{\text {trial }} \geqslant 0 $, enter the plastic deformation calculation;
(c) Calculate the plastic consistency variable increment Δγ;
(d) Calculate the Cauchy stress;
(e) Update the true equivalent plastic strain;
(f) Output displacement, velocity, acceleration;
8. Calculate the rotational Cauchy stress (true stress), σ;
9. Calculate the firstorder RungeKutta stress tensor, P;
10. Cycle all material points x_{i} within the particle neighborhood δ;
11. Calculate the "force state", $\underline{\boldsymbol{T}} $;
12. End loop.
The flowchart of the entire computation of the peridynamics model and the constitutive update for the proposed ice constitutive model are shown in Figure 6.
3 Results and discussions
3.1 Drop tower test
In order to verify the validity and applicability of the proposed ice constitutive model, an ice cylinder is simulated to hit a steel plate at a given impact velocity. The geometrical model size is as shown in Figure 7. The material 1 is chosen to be ice material, and the material 2 is defined as rigid body. There are 74 225 particles and 66 436 cells discrete in total. The mesh size is Δx=1.5 mm, and the horizon size is δ =3.15Δx. Assuming that the cylindrical ice moves at a constant speed, the rigid body does not deform during the collision with the ice body, and the given cylindrical ice provides the conversion of kinetic energy at a certain speed. The calculated material parameters are shown in Table 2.
Methods T(℃) ε(s^{1}) a_{0} a_{1} a_{2} DerradjiAouat (2000) 1 4×10^{3} 22.93 2.06 0.002 3 Kierkegaard ― ― 2.588 8.63 0.163 Riska and Frederking (1987) data 1 2 2×10^{3} 1.60 4.26 0.62 Riska and Frederking (1987) data 2 10 2×10^{3} 3.1 9.20 0.83 The damage counter of the destruction process of the cylindrical ice with a collision velocity of 152 m/s is shown in Figure 8. The damage factor in the figure is in the range of 0.05−0.95. It can be seen from Figure 8 that the initial damage occurs at the position where the cylindrical ice contacts the cylindrical rigid body, and the ice damage gradually expands upward. Pernas et al. (2012) used three numerical methods: SPH method, ALE method and Largrange method. The predictions of these methods are focused on the deformation of the ice body, and there is no good reflection of the damage and cracks of the ice body. From the calculation results in this paper, it can be seen that the prediction effect of the cylindrical ice damage is that the cracks grow upward from the contact with the rigid body. From 8, the crack propagation and damage process of the cylindrical ice can be clearly shown.
Figure 9 illustrates the curve of the collision force between cylindrical ice and rigid body over time. The simulated results are compared with the experimental results of Carney et al. (2006) and the numerical results of (PernasSánchez et al., 2012). The initial speed conditions of 152 m/s and 213 m/s were selected respectively. It can be seen from Figure 9 that the current results is consistent with the experimental results within the time range of 0 ms to 0.1 ms at the initial time (the collision velocity is 152 m/s). The predicted collision force is larger than the experimental one. Before reaching the peak of the collision force, there is a certain forward effect. No elastic deformation occurs when contact occurs, so the calculation result is too large. The collision force reaches its peak at about 0.15 ms, and then gradually decreases as the energy is dissipated. Under the condition of the collision velocity of 213 m/s adopted in Figure 9(b), the calculation results are basically consistent with the experimental results, and the error is smaller than the calculation result of the collision velocity of 152 m/s adopted in Figure 9(a). To a certain extent, the calculation method in this paper is more suitable for highspeed collision conditions, thus verifying the feasibility of the nonordinary statebased peridynamics constitutive equation used in the calculation of torsional deformation of brittle ice structures. It can be seen from the calculation results of the numerical simulation methods SPH method, ALE method and Largrange method adopted by Pernas that the Peridynamics calculation method used in this paper has certain advantages over other numerical simulation methods.
3.2 Ballistic problem
To investigate the propagation characteristics of ice cracks from the mesoscale perspective, the size of the plate model established here is selected as 0.02 m×0.02 m× 0.001 m thin ice plate, the diameter of the steel ball is 0.001 5 m. A total of 62 203 particles and 50 864 elements are discretized. The rigid ball is assumed to be a rigid body moving at a constant speed, which will not deform during the collision with the ice body. The given rigid ball provides the conversion of kinetic energy at a certain speed. The calculation model is shown in Figure 10. The damage process of ice at a velocity of 150 m/s is shown in Figure 11. The comparison result of the calculated damage contour is shown in Figure 12.
The proposed constitutive has the same effectiveness in processing torsional deformation of ice material. The yield curve equation of the multiyield surface constitutive model is determined by the parameter constants obtained on the basis of a large number of experiments. It can simulate the plastic deformation of anisotropic materials, and it can simulate the ice failure in the brittle regime as well. Due to the complexity of the structure and composition of the ice material itself, and the different growth and formation processes of the ice, the crystal arrangement is also diverse, so the crack shape of the ice material under high strain rate also presents different manifestations, and has a larger the randomness. This paper selects several experimentally photographed crack shapes during ice breakage from the literature, and compares them with the numerical results in this paper, as shown in Figure 12. In Figure 12, the red area indicates the broken ice body, and the yellow area indicates the ice material in the elastoplastic stage. Compare with the experimental result, the numerical model can better simulate the crack shape. The crack appears to be broken from the impact center of the ice plate, and the debris area around the center area is approximately circular and runs along the center point, striped cracks radiated all around. The result shows the multiyield ice constitutive model can better simulate the characteristics of the diversity of ice cracks.
3.3 Shipice interaction
In order to verify the applicability of the proposed ice constitutive model in simulating shipice interaction, the icebreaker is used. Assume that the ship is a rigid body and sails through ice at a certain speed. The ship dimensions are chosen in Table 3. The numerical model is built by Abaqus software, and there are 184 225 particles and 126 436 cells discrete in total, as is shown in Figure 13. The mesh size is Δx=1 m, and the horizon size is δ=3.15Δx. Assuming that the icebreaker moves at a constant speed at 5 m/s, the rigid body does not deform during the collision with the ice sheet, and the given ice sheet provides the conversion of kinetic energy at a certain speed. The calculated ice material parameters are the same as Table 2.
Length L Width B Depth D Design draft T L_{wl} 84 22 11 6 78 It can be seen from Figures 14‒18 that the proposed method can well capture the dynamic response of the ice, and shows good performance to simulate shipice interaction.
4 Conclusions
In this paper, we developed a statebased cohesive peridynamics model that incorporates a coupled multiyield surface elastoplasticity model to simulate the stress and crack formation of ice under impact. Compared with the traditional ice constitutive model, the main advantages of the proposed method are:
1) We incorporate the multi yieldsurface elastoplasticity model for ice (Derradji, 2002) into the statebased peridynamics formulation, which is then successfully being used to simulate and predict failure of fresh water isotropic ice, iceberg ice, and columnar saline ice.
2) The simulated drop tower test demonstrates that the proposed peridynamics model can better simulate the ice fragmentation under impact loads accurately, and it captures the elasticplastic quasibrittle failure behavior in ice, and
3) The proposed method can well capture the ice compressive and tension behaviors in the high strain rate regime, and it shows good performance in simulating shipice collision in ocean engineering applications.
Despite of these advantages, there are still some limitations in the present approach. For example, the influence of the added mass by water is ignored in the present peridynamics model when simulating the shipice interaction, and moreover, the ship motion and the subsequent fluidstructure interaction are not considered. Such issues will be addressed in future work.
Nomenclature ξ_{ij} Bond vector in reference configuration η_{ij} The relative displacement vector in current configuration $ \underline{\boldsymbol{Y}}<\cdot>$ State function H_{xi} The horizon of particle x_{i} $ \underline{\boldsymbol{T}}\left(\xi_{i j}, \underline{\boldsymbol{Y}}\left(\xi_{i j}\right)\right)$ The forcevector state K_{xi} The reference symmetric shape tensor of particle x_{i} P_{xi} The first PiolaKirchhoff stress F_{xi} The nonlocal deformation gradient L The spatial velocity gradient D Symmetric deformation rate tensor W Asymmetric rotation tensor d The nonrotational deformation rate tensor R_{t} Orthogonal rotation tensor f The yield function Competing interest The authors have no competing interests to declare that are relevant to the content of this article.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author (s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 
Figure 3 Illustration of stress return in the TsaiWu yield surface (DerradjiAouat, 2003)
Figure 7 Discretized Peridynamic model (Carney et al., 2006)
Figure 18 Comparison between the peridynamics simulation results and experimental data (Xu et al., 2020)
Table 1 Elliptical failure envelope parameters for three types of ice (DerradjiAouat, 2003)
Types η λ P_{max} Freshwater ice 0 45.0 55.0 Iceberg ice 0 45.0 55.0 Grown ice 0 45.0 45.0 Table 2 Input parameter
Methods T(℃) ε(s^{1}) a_{0} a_{1} a_{2} DerradjiAouat (2000) 1 4×10^{3} 22.93 2.06 0.002 3 Kierkegaard ― ― 2.588 8.63 0.163 Riska and Frederking (1987) data 1 2 2×10^{3} 1.60 4.26 0.62 Riska and Frederking (1987) data 2 10 2×10^{3} 3.1 9.20 0.83 Table 3 Dimension of the icebreaker (m)
Length L Width B Depth D Design draft T L_{wl} 84 22 11 6 78 Nomenclature ξ_{ij} Bond vector in reference configuration η_{ij} The relative displacement vector in current configuration $ \underline{\boldsymbol{Y}}<\cdot>$ State function H_{xi} The horizon of particle x_{i} $ \underline{\boldsymbol{T}}\left(\xi_{i j}, \underline{\boldsymbol{Y}}\left(\xi_{i j}\right)\right)$ The forcevector state K_{xi} The reference symmetric shape tensor of particle x_{i} P_{xi} The first PiolaKirchhoff stress F_{xi} The nonlocal deformation gradient L The spatial velocity gradient D Symmetric deformation rate tensor W Asymmetric rotation tensor d The nonrotational deformation rate tensor R_{t} Orthogonal rotation tensor f The yield function 
Carney KS, Benson DJ, DuBois P, Lee R (2006) A phenomenological high strain rate model with failure for ice. International Journal of Solids and Structures 43(25–26): 7820–7839 https://doi.org/10.1016/j.ijsolstr.2006.04.005 Cui P, Zhang AM, Wang S, Khoo BC (2018) Ice breaking by a collapsing bubble. Journal of Fluid Mechanics 841: 287–309 https://doi.org/10.1017/jfm.2018.63 DerradjiAouat A (2003) Multisurface failure criterion for saline ice in the brittle regime. Cold Regions Science and Technology 36 (1–3): 47–70 https://doi.org/10.1016/S0165232X(02)000939 DerradjiAouat A (2000) A unified failure envelope for isotropic fresh water ice and iceberg ice. Proceedings of ETCE/OMAE Joint Conference, New Orleans Drucker DC, Prager W (1952) Soil mechanics and plastic analysis or limit design. Quarterly of Applied Mathematics 10(2): 157–165 https://doi.org/10.1090/qam/48291 Fan H, Bergel GL, Li S (2016) A hybrid peridynamicssph simulation of soil fragmentation by blast loads of buried explosive. International Journal of Impact Engineering 87: 14–27 https://doi.org/10.1016/j.ijimpeng.2015.08.006 Fan H, Li S (2017a) A peridynamicssph modeling and simulation of blast fragmentation of soil under buried explosive loads. Computer Methods in Applied Mechanics and Engineering 318: 349–381 https://doi.org/10.1016/j.cma.2017.01.026 Fan H, Li S (2017b) Parallel peridynamicsSPH simulation of explosion induced soil fragmentation by using openmp. Computational Particle Mechanics 4(2): 199–211 https://doi.org/10.1007/s4057101601165 Fan H, Ren B, Li S (2015) An adhesive contact mechanics formulation based on atomistically induced surface traction. Journal of Computational Physics 302: 420–438 https://doi.org/10.1016/j.jcp.2015.08.035 Gold LW (1988) On the elasticity of ice plates. Canadian Journal of Civil Engineering 15(6): 1080–1084 https://doi.org/10.1139/l88140 Han R, Zhang AM, Tan S, Li S (2022) Interaction of cavitation bubbles with the interface of two immiscible fluids on multiple time scales. Journal of Fluid Mechanics 932: A8 https://doi.org/10.1017/jfm.2021.976 Hu Y, Feng G, Li S, Sheng W, Zhang C (2020) Numerical modelling of ductile fracture in steel plates with nonordinary statebased peridynamics. Engineering Fracture Mechanics 225: 106446 https://doi.org/10.1016/j.engfracmech.2019.04.020 Hughes T, Winget J (1980) Finite rotation effects in numerical integration of rate constitutive equations arising in large deformation analysis. International Journal for Numerical Methods in Engineering 15(12): 1862–1867 https://doi.org/10.1002/nme.1620151210 Jia B, Ju L, Wang Q (2019) Numerical simulation of dynamic interaction between ice and wide vertical structure based on peridynamics. Computer Modeling in Engineering & Sciences 121(2): 501–522 Johnson GR, Holmquist TJ (1994) An improved computational constitutive model for brittle materials. AIP Conference Proceedings, 309: 981–984 https://doi.org/10.1063/1.46199 Jones SJ (1982) The confined compressive strength of polycrystalline ice. Journal of Glaciology 28(98): 171–178 https://doi.org/10.3189/S0022143000011874 Jones SJ (1997) High strainrate compression tests on ice. The Journal of Physical Chemistry B 101(32): 6099–6101 https://doi.org/10.1021/jp963162j Kolsky H (1949) An investigation of the mechanical properties of materials at very high rates of loading. Proceedings of the Physical Society, Section B 62(11): 676 https://doi.org/10.1088/03701301/62/11/302 Lai X, Liu L, Li S, Zeleke M, Liu Q, Wang Z (2018) A nonordinary statebased peridynamics modeling of fractures in quasibrittle materials. International Journal of Impact Engineering 111: 130–146 https://doi.org/10.1016/j.ijimpeng.2017.08.008 Li T, Zhang AM, Wang SP, Li S, Liu WT (2019) Bubble interactions and bursting behaviors near a free surface. Physics of Fluids 31(4): 042104 https://doi.org/10.1063/1.5088528 Liu M, Wang Q, Lu W (2017) Peridynamic simulation of brittleice crushed by a vertical structure. International Journal of Naval Architecture and Ocean Engineering 9(2): 209–218 https://doi.org/10.1016/j.ijnaoe.2016.10.003 Liu NN, Zhang AM, Cui P, Wang SP, Li S (2021) Interaction of two outofphase underwater explosion bubbles. Physics of Fluids 33(10): 106103 https://doi.org/10.1063/5.0064164 Lu W, Li M, Vazic B, Oterkus S, Oterkus E, Wang Q (2020) Peridynamic modelling of fracture in polycrystalline ice. Journal of Mechanics 36(2): 223–234 https://doi.org/10.1017/jmech.2019.61 Madenci E, Oterkus E (2014) Peridynamic theory. In: Peridynamic theory and its applications, 19–43 Mellor M, Cole DM (1982) Deformation and failure of ice under constant stress or constant strainrate. Cold Regions Science and Technology 5(3): 201–219 https://doi.org/10.1016/0165232X(82)900155 Palmer A, Dempsey J, Masterson D (2009) A revised ice pressurearea curve and a fracture mechanics explanation. Cold Regions Science and Technology 56(2–3): 73–76 https://doi.org/10.1016/j.coldregions.2008.11.009 PernasSánchez J, ArteroGuerrero JA, Varas D, LópezPuente J (2015) Analysis of ice impact process at high velocity. Experimental Mechanics 55(9): 1669–1679 https://doi.org/10.1007/s1134001500674 PernasSánchez J, Pedroche DA, Varas D, LópezPuente J, Zaera R (2012) Numerical modeling of ice behavior under high velocity impacts. International Journal of Solids and Structures 49(14): 1919–1927 https://doi.org/10.1016/j.ijsolstr.2012.03.038 Riska K, Frederking R (1987) Ice load penetration modelling. Proceedings of the Ninth Port and Ocean Engineering Under Arctic Conditions Conference, Fairbanks, 1: 317–327 Rubinstein R, Atluri S (1983) Objectivity of incremental constitutive relations over finite time steps in computational finite deformation analyses. Computer Methods in Applied Mechanics and Engineering 36(3): 277–290 https://doi.org/10.1016/00457825(83)901251 Sain T, Narasimhan R (2011) Constitutive modeling of ice in the high strain rate regime. International Journal of Solids and Structures 48(5): 817–827 https://doi.org/10.1016/j.ijsolstr.2010.11.016 Schulson E, Buck S (1995) The ductiletobrittle transition and ductile failure envelopes of orthotropic ice under biaxial compression. Acta Metallurgica et Materialia 43(10): 3661–3668 https://doi.org/10.1016/09567151(95)901493 Schulson EM (2001) Brittle failure of ice. Engineering Fracture Mechanics 68(17–18): 1839–1887 https://doi.org/10.1016/S00137944(01)000376 Silling SA (2000) Reformulation of elasticity theory for discontinuities and longrange forces. Journal of the Mechanics and Physics of Solids 48(1): 175–209 https://doi.org/10.1016/S00225096(99)000290 Snyder SA, Schulson EM, Renshaw CE (2016) Effects of prestrain on the ductiletobrittle transition of ice. Acta Materialia 108: 110–127 https://doi.org/10.1016/j.actamat.2016.01.062 Sun PN, Colagrossi A, Zhang AM (2018) Numerical simulation of the selfpropulsive motion of a fishlike swimming foil using the d^{+}SPH model. Theoretical and Applied Mechanics Letters 8(2): 115–125 https://doi.org/10.1016/j.taml.2018.02.007 Sun PN, Le Touze D, Oger G, Zhang AM (2021) An accurate FSISPH modeling of challenging fluidstructure interaction problems in two and three dimensions. Ocean Engineering 221, 108552 Sun PN, Luo M, Le Touzé D, Zhang AM (2019) The suction effect during freak wave slamming on a fixed platform deck: Smoothed particle hydrodynamics simulation and experimental study. Physics of Fluids 31(11): 117108 https://doi.org/10.1063/1.5124613 Vazic B, Oterkus E, Oterkus S (2020) Inplane and outof plane failure of an ice sheet using peridynamics. Journal of Mechanics 36(2): 265–271 https://doi.org/10.1017/jmech.2019.65 Wang G, Ji SY, Lv HX, Yue QJ (2006) Druckerprager yield criteria in viscoelasticplastic constitutive model for the study of sea ice dynamics. Journal of Hydrodynamics 18(6): 714–722 https://doi.org/10.1016/S10016058(07)600110 Wang Q, Wang Y, Zan Y, Lu W, Bai X, Guo J (2018) Peridynamics simulation of the fragmentation of ice cover by blast loads of an underwater explosion. Journal of Marine Science and Technology 23(1): 52–66 https://doi.org/10.1007/s007730170454x Xie Y, Li S (2021a) A stressdriven computational homogenization method based on complementary potential energy variational principle for elastic composites. Computational Mechanics 67: 637–652 https://doi.org/10.1007/s00466020019538 Xie Y, Li S (2021b) Finite temperature atomisticinformed crystal plasticity finite element modeling of single crystal tantalum (ata) at micron scale. International Journal for Numerical Methods in Engineering 122(17): 4660–4697 https://doi.org/10.1002/nme.6741 Xie Y, Li S, Hu X, Bishara D (2022b) An adhesive gurtinmurdoch surface hydrodynamics theory of moving contact line and modeling of droplet wettability on soft substrates. Journal of Computational Physics 456: 111074 https://doi.org/10.1016/j.jcp.2022.111074 Xie Y, Li S, Wu C, Lyu D, Wang C, Zeng D (2022a) A generalized bayesian regularization network approach on characterization of geometric defects in lattice structures for topology optimization in preliminary design of 3D printing. Computational Mechanics 69(5): 1191–1212 https://doi.org/10.1007/s00466021021378 Xu Y, Kujala P, Hu ZQ, Li F, Chen G (2020) Numerical simulation of level ice impact on landing craft bow considering the transverse isotropy of Baltic Sea ice based on XFEM. Marine Structures 71: 102735 https://doi.org/10.1016/j.marstruc.2020.102735 Zhang AM, Li SM, Cui P, Li S, Liu YL (2023) A unified theory for bubble dynamics. Physics of Fluids 35(3): 033323 https://doi.org/10.1063/5.0145415 Zhang AM, Sun PN, Ming FR, Colagrossi A (2017) Smoothed particle hydrodynamics and its applications in fluidstructure interactions. Journal of Hydrodynamics, Ser. B 29(2): 187–216 https://doi.org/10.1016/S10016058(16)607308 Zhang LW, Xie Y, Lyu D, Li S (2019) Multiscale modeling of dislocation patterns and simulation of nanoscale plasticity in bodycentered cubic (BCC) single crystals. Journal of the Mechanics and Physics of Solids 130: 297–319 https://doi.org/10.1016/j.jmps.2019.06.006