2.1 Peridynamic model for elastic sea ice

PD is a non-local formulation of continuum mechanics that is oriented towards deformations with discontinuities, especially fractures. Unlike spatial derivatives common in classic continuum mechanics (e.g. the concept of stress and strain), PD employs integral operators to represent these phenomena. The discretization of the domain in PD involves particles that interact with their neighbouring particles within a specific distance known as the horizon (H _x). Figure 1 provides a visual representation of the fundamental principle of PD theory, where Particle i at position ${\bf x}$ [m] interacts with Particle j at position ${{\bf x}^{\prime}}$[m]. Under the influence of external forces, the particle i experiences displacement ${\bf u}$ [m] while the particle j undergoes displacement ${{\bf u}^{\prime}}$ [m]. Consequently, this deformation induces a force state term ${\bf t}$, which has a unit of [N/m⁶] acting on the particle i and another force state term ${{\bf t}^{\prime}}$ [N/m⁶] acting on the particle j. In PD theory, the equation of motion for a material point is defined in Eqn (1) (Madenci and Oterkus, Reference Madenci and Oterkus2014). The previously introduced ‘force state’ terms ${\bf t}$ and ${{\bf t}^{\prime}}$, after a volume integration, yield a concept of force density [N/m³] on both sides of Eqn (1).

(1)$$\eqalignno{\rho ( {\bf x} ) \ddot{{\bf u}}( {{\bf x}, \;t} ) & = \mathop \int \limits_{H_x}^{} ( {{\bf t}( {{{\bf u}^{\prime}}-{\bf u}, \;{{\bf x}^{\prime}}-{\bf x}, \;t} ) -{\bf t}{\rm^{\prime}}( {{\bf u}-{\bf u}{\rm^{\prime}}, \;{\bf x}-{\bf x}{\rm^{\prime}}, \;t} ) } ) dV{\rm ^{\prime}} \cr& \quad+ {\bf b}( {{\bf x}, \;t} ) , \;}$$

Figure 1. Discretization and particle interactions in PD theory (Zhang and others, Reference Zhang, Tao, Wang, Ye and Sun2021a, Reference Zhang, Tao, Wang, Ye and Guo2021b).

In Eqn (1), $\rho ( {\bf x})$ [kg/m³] represents the density of the sea ice, $\ddot{{\bf u}}( {{\bf x}, \;t} )$ [m/s²] represents the acceleration of the discretized ice particle, and ${\bf b}( {\bf x}, \;t)$ [N/m³] is a body force. In the present work, the fracture problem of 2-D ice floes is studied, and the Poisson's ratio of sea ice equals to 0.33. Consequently, we employ the bond-based PD equation rather than the general state-based PD for numerical simulation. The rationale behind this choice lies in the fact that the force density vectors ${\bf t}$ and ${{\bf t}^{\prime}}$ in the bond-based PD are equal in magnitude and parallel to the relative position vector. Under this assumption, the Poisson's ratio (ν) is fixed at 1/3 for 2D simulations, which suits our requirements for simulating sea ice. In contrast, in the general state-based PD, the force density vector is unconstrained, allowing for the free setting of Poisson's ratio. However, state-based PD is computationally expensive.

The expression of the force state ${\bf t}$ [N/m⁶] for an elastic and isotropic ice material in Eqn (1) is

(2)$${\bf t} = {-}{{\bf t}^{\prime}} = 2\delta sb\displaystyle{{{{\bf y}^{\prime}}-{\bf y}} \over {\vert {{{\bf y}^{\prime}}-{\bf y}} \vert }}.$$

In the above formula, b represents a parameter of PD (presented in Eqn (4)), δ [m] represents the size of the horizon (H _x). ${\bf y}$ [m] and ${{\bf y}^{\prime}}$ [m] are the position vectors of particle i and particle j after deformation, respectively. Generally, the shape of the horizon is a circular shape in 2D (disc). Therefore, horizon size δ refers to the radius of a disc, as shown in Figure 1. s is a unitless scalar, called the bond stretch, representing the deformation between two particles (see Eqn (3)). The relation between bond stretch and the ‘force state’ in Eqn (2) demonstrates the constitutive model of linear elastic ice material.

(3)$$s = \displaystyle{{\vert {{{\bf y}^{\prime}-y}} \vert -\vert {{{\bf x}^{\prime}-x}} \vert } \over {\vert {{{\bf x}^{\prime}-x}} \vert }}.$$

b [Pa/m⁵] is a scalar related to a bond constant c [Pa/m⁴] which usually is expressed by ice properties (elastic modulus, E [Pa]), and horizon size δ [m]. The relations are:

(4)$$b = \displaystyle{c \over {4\delta }}\;{\rm with}\;c = \displaystyle{{9E} \over {\pi h\delta ^3}}{\rm for\ 2D\ with\ }\upsilon = \displaystyle{1 \over 3}, \;$$

where h [m] is the thickness of the ice floe.

2.2 Material failure model

In PD theory, material failure simulation can be viewed at two levels (see Fig. 2). The first level is the particle–particle interaction elimination through a binary damage function (i.e. 0 or 1). For example, in the ‘crack initiation phase’ in Figure 2, for particle i = 1, its bonds 1–11 and 1–12 are assigned a value of 0 and are thus eliminated. The failure at this level does not necessarily entail a complete failure for a material point (e.g. for particle i = 1 in Fig. 2, there are still 10 particles remaining connected to it). The second level failure reminiscent a domain damage concept, in which, a percentage of particle–particle interactions within the horizon of a material point are eliminated. This percentage (from 0 to 1) is termed as a damage variable and is a continuous function represents the level of ‘damage’, with ‘1’ representing the complete break-off of a material point (e.g. see the total elimination of particle i = 1 in the last phase of Fig. 2). Any other value in between 0–1 gives us the possibility to characterize the location of a macroscopic crack, for example with a damage variable of 0.5 (e.g. see the crack formation phase in Fig. 2), it represents a crack cutting through the material point, whose 50% particle–particle connections in the horizon have been eliminated.

Figure 2. Failure process and its mathematical expression in PD theory.

Mathematically, the above damaging processes are achieved through the introduction of two additional functions. At the first level failure, particle–particle interaction can be eliminated through the concept of bond rupture. Bond rupture can be determined using a variety of failure criteria. We choose to use the critical stretch criterion. It is a straightforward and widely used criterion to predict material failure; and is particularly derived for tensile failures. This criterion involves comparing the stretch s (in Eqn (3)) between two material particles with a critical value, which is called the critical stretch (s _C).

(5)$$\Omega ( t, \;{{\bf x}^{\prime}}-{\bf x}) = \left\{{\matrix{ {1, \;s < s_C{\rm , \;\ unbroken\ or\ visible\ bond\ }} \cr {0, \;s \ge s_C, \;{\rm broken\ or\ invisible\ bond\ }} \cr } } \right..$$

A binary failure function $\Omega ( t, \;{{\bf x}^{\prime}}-{\bf x})$ is introduced in accordance with the critical stretch criterion in Eqn (5). This function is incorporated in Eqn (1) to characterize material failure in Eqn (6) as (Silling and Askari, Reference Silling and Askari2005; Madenci and others, Reference Madenci, Roy and Behera2022)

(6)$$\eqalignno{\rho ( {\bf x} ) \ddot{{\bf u}}( {{\bf x}, \;t} ) \ =& \mathop \int \limits_{H_x}^{} ( {{\bf t}( {{{\bf u}^{\prime}}-{\bf u}, \;{{\bf x}^{\prime}}-{\bf x}, \;t} ) -{\bf t}^{\prime}( {{\bf u}-{\bf u}^{\prime}, \;{\bf x}-{\bf x}^{\prime}, \;t} ) } ) \Omega ( {t, \;{\bf x}^{\prime}-{\bf x}} ) dV^{\prime} \cr& \quad +\ {\bf b}( {{\bf x}, \;t} ) .}$$

$\Omega ( t, \;{{\bf x}^{\prime}}-{\bf x})$ defines the magnitude of an irreversible failure, thereby changing the load distribution within the body and allowing crack initiation and propagation. With this criterion, crack propagation occurs spontaneously without requiring a predefined direction.

In Eqn (5), the unitless critical stretch s _C is derived from the critical energy release rate G _c [N/m]. The 2D bond-based version of s _C is expressed as (Madenci and others, Reference Madenci, Roy and Behera2022)

(7)$$s_C = \sqrt {\displaystyle{{\pi G_c} \over {3K_{2D}\delta }}} , \;$$

in which K _2D [Pa] is bulk modulus in 2D.

Failure function in Eqn (5) only defines the interaction elimination between two particles. In the second level, the determination of crack initiation and propagation involves the introduction of a PD damage variable (denoted as $\varphi ( {\bf x}, \;t)$, see Fig. 2), which integrates all the failure function in the horizon of particle of interest. Therefore, $\varphi ( {\bf x}, \;t)$ is simply a weight ratio between the broken bonds and the total number of initial bonds connecting a material particle within the horizon. At crack initiation, a particle engages in interactions with all ice particles within its horizon, resulting in a local damage variable equalling to 0 (i.e. no damage or no crack initiation). On the other hand, the formation of a crack surface leads to the elimination of half of the interactions within the horizon, yielding a local damage variable equalling to 0.5. The damage variable is presented in Eqn (8).

(8)$$\varphi ( {\bf x}, \;t) = 1-\displaystyle{{\int_{H_x} {\Omega ( t, \;{{\bf x}^{\prime}}-{\bf x}) } {\rm d}{V}^{\prime}} \over {\int_{H_x} {} {\rm d}{V}^{\prime}}}.$$

2.3 Discretization and solvers

Our approach primarily relies on the open-source PD software Peridigm (Littlewood and others, Reference Littlewood, Parks, Foster, Mitchell and Diehl2023). We also utilize a finite-element mesh generator to construct a hexahedral mesh of the ice domain, which is then saved as an external file for PD model construction purposes. As mentioned before, the 2D bond-based explicit PD solver is employed. This solver evaluates the force states (denoted as ${\bf t}$ and ${{\bf t}^{\prime}}$) in Eqn (1) at each time step and applies them to each bond in the discrete model. The time steps are determined based on the criterion in Bobaru and others (Reference Bobaru, Foster, Geubelle and Silling2016), that is

(9)$$\Delta t_{critical} = \sqrt {\displaystyle{{2\rho } \over {\sum\limits_p {\Delta V_p( 18K_{2D}/( {{\bf x}^{\prime}}-{\bf x}) \pi \delta ^4) } }}} , \;$$

in which ρ is the density, p iterates over all the neighbours of the given material point, ΔV _p is the volume associated with neighbours p. Then, we utilize Paraview for post-processing.

3.1 Model set-up

Figure 3 illustrates the model set-up, in which each mesh element represents a spatial volume occupied by PD particle. The solution of the governing equation (Eqn (6)) is obtained by performing a volume integral over all hexahedron volumes in the horizon. To compare with existing analytical solutions, the ECRP is modelled as homogeneous, isotropic and linear elastic in the PD model. This ice plate is also assumed to fail elastically, where LEFM can be applied. However, it should be noted that PD is not limited to these idealizations (Zhang and others, Reference Zhang, Tao, Wang, Ye and Guo2021b).

Figure 3. Model set-up for the in-plane splitting of an ECRP.

Following the approach in the study conducted by Vazic and others (Reference Vazic, Oterkus and Oterkus2020), a constant particle body force is assumed as the loading condition in the current study and is applied to the three particles in the volume they occupy along the pre-crack, as shown in Figure 3. The magnitude of this constant load is expected to be small (to minimize dynamic effects), yet large enough to propagate the central crack. In post processing, we extract the splitting force F _y for comparison with analytical solutions. Since we have performed a load-controlled simulation here, we took an unconventional but equivalent approach to extract the splitting force. In detail, we integrate bond–bond forces along the remaining ligament of the central crack. This gives us the ‘fracture resistance’ in force term. In this way, we manage to extract the splitting force F _y, which is in equilibrium of the ‘fracture resistance’.

The particle spacing (dx) and the horizon size (δ) are critical factors influencing the computational process, establishing the optimal values for these parameters is crucial to obtain accurate results. Before directly presenting the results and comparison, a series of convergence studies (δ-convergence and dx-convergence) are performed in the following.

3.2 Convergence analysis

The δ-convergence and dx-convergence are first studied to determine the optimal parameters for the current PD simulations. We choose dx as 0.1, 0.15, 0.2, 0.3, 0.5 m for the dx-convergence study. The horizon size is δ = mdx, where m denotes a multiplier factor. We chose δ = mdx = (2.015, 3.015, 4.015)dx for the δ-convergence study. Discretization and critical stretches calculated by Eqn (7) are shown in Table 2.

Table 2. Information for discretization and the critical stretch in ice tensile failure

Convergence study results are presented in Figures 4a and 4b. In these figures, we presented the normalized splitting force (F _y/(hK _ICL ^1/2), in which $K_{IC} = \sqrt {EG_c}$ represents fracture toughness of sea ice vs the normalized crack length (A/L). Simulation outputs from PD are generated in the time domain. Normalized splitting force vs normalized crack length in Figure 4 is constructed by synchronizing the simulated force and crack length history.

Figure 4. Normalized splitting force vs normalized crack length for the splitting of an ECRP: (a) dx -convergence and (b) δ -convergence study with PD method.

Based on Figure 4a, it can be observed that varying particle spacings has minimal impact on the splitting force value. All curves reach their maximum value when A/L is between 0.14 and 0.19. However, as the particle spacing increases, the oscillation amplitude of the splitting force also increases. Notably, when the particle distances are set at 0.1 m and 0.15 m, the vibration amplitude of the dimensionless splitting force is the smallest, and all curves converge at around these values. Consequently, a particle distance of 0.15 m is chosen as the optimal setting for accuracy considerations.

δ-Convergence analysis plays a crucial role in determining the appropriate horizon size, which greatly affects the simulation results. As depicted in Figure 4b, it is evident that the results are in close agreement when the multiplier factors m for the horizon size are set to 3.015 and 4.015. Notably, the multiplier factor of 3.015 is widely used in various other research fields (Madenci and Oterkus, Reference Madenci and Oterkus2014), further affirming its reliability.

The convergence results are encouraging. Especially, it has revealed an optimal setting, that is, m = 3.015 and dx = 0.15 m for our benchmark test (to be discussed in Section 5.1).

3.3 Simulation results of the splitting of ECRP

Fixing the multiplier factor m = 3.015 and dx = 0.15 m, the simulated results of crack path for ice tensile failure are depicted in Figure 5. The central edge crack propagates straight ahead, as theoretically expected.

Figure 5. Crack path of the splitting of an ECRP: initial snapshot (left) and crack propagation snapshot (right).

Then, the normalized splitting force of PD simulation is compared with analytical solution found in Lu and others (Reference Lu, Lubbad, Høyland and Løset2014a), numerical results obtained by Lu and others (Reference Lu, Lubbad and Løset2015c), and FEM results obtained by Bhat (Reference Bhat1988) are shown in Figure 6.

Figure 6. Comparison of the 2D numerical PD scheme with the analytical solution and other numerical methods for the benchmark test: normalized splitting force vs dimensionless crack length.

The normalized splitting force initially increases and then decreases as the normalized crack length increases, reaching a peak when the normalized crack length is between 0.14 and 0.2. This trend is consistent with the results from the compression of an edge-cracked rectangular plate conducted by Hallam and others (Reference Hallam, Jones and Howard1989). Moreover, we compare the calculated results with those calculated by other methods for the same case. The results depicted in Figure 6 indicate that the calculated values for the ice plate's splitting force from the 2D bond-based PD exhibit a reasonably satisfactory agreement with those methods. Further discussions will be presented in Section 5.2. Next, we will explore PD's capability in a more realistic tensile fracture scenario, that is, multiple fracturing of a heterogeneous ice floe.

4.1 Model set-up

Following the work of Dinh and others (Reference Dinh, Giannakis, Slawinska and Stadler2023), we built a same square ice floe (1 km in length and 1 m in thickness) in the PD discretization domain (see Fig. 7). All the input parameters for simulation are listed in Table 3. A prescribed boundary condition with a displacement of 5 mm is applied on the left edge of ice floe. The right edge is fixed (see Fig. 7a). Note that the boundary conditions are imposed by three layers of fictional boundary particles as required by the PD theory. These layers consist of particles that occupy a volume of space and have the same parameters as other particles in the ice body. However, these particles are not within the actual geometry of the ice and are added additionally, as shown in Figure 7(a).

Figure 7. Ice floe model for fracture simulation (Each ice floe has ten different weak zones (#1–10) with different thickness.): (a) illustration of boundary conditions of the heterogenous ice floe fracture; (b) the heterogenous ice floe for case 1; (c) the heterogenous ice floe for case 2.

Table 3. Input parameters for numerical experiment of heterogenous ice floe fracture

Varying ice thickness values are assigned to different zones matching those of literature (Dinh and others, Reference Dinh, Giannakis, Slawinska and Stadler2023). Two ice floes are modelled in this study (see Figs 7b and 7c). Corresponding ice thickness information is presented in Table 4.

Table 4. Thickness distribution of the two ice floe models

4.2 Simulation results on the fracture of a heterogeneous ice floe

Figure 8 illustrates the displacement contour and damage variable contour (i.e. φ value in Eqn (8)) computed using the PD method, as well as the displacement and damage result obtained by Dinh and others (Reference Dinh, Giannakis, Slawinska and Stadler2023) through PFM modelling.

Figure 8. Comparison of the simulation results of the fracture of a heterogeneous ice floe by PD and PFM: (a-1) case 1, PD displacement result, unit [m]; (a-2) case 1, PD damage result, unit [-]; (a-3) case 1, PFM displacement and damage result obtained by Dinh and others (Reference Dinh, Giannakis, Slawinska and Stadler2023) (black arrows represent the displacement field over the ice floe, while red line indicates crack); (b-1) case 2, PD displacement result; (b-2) case 2, PD damage result; (b-3) case 2, PFM displacement and damage result obtained by Dinh and others (Reference Dinh, Giannakis, Slawinska and Stadler2023) (black arrows represent the displacement field over the ice floe, while red line indicates crack).

As Figure 8 demonstrates, our simulated crack paths agree well with the results from PFM: (1) both ice floes were fractured into two pieces. (2) The cracks predominantly initiate from weak zones (i.e. thin ice) and propagate in parallel to the boundary, aligning with any existing defects. The area of non-zero motion in Figure 8(b-1) appears because this area is almost completely separated from the rest of the ice plate. Only a small portion of the particles in this area interact with its horizon particles, which are located inside the more intact part of the ice body. Consequently, the velocity and displacement of particles in this area are subject to minimal constraints, leading to rapid changes in speed and displacement during numerical calculations. Please also note that positive displacement is to the right and negative displacement is to the left. As you can see, there is a white area without a damage value in Figures 8(b-1) and (b-2). This is because zero thickness is represented by an absence of material, which is why there is no damage value in this gap area.

Additionally, PD simulates crack branches (as shown in case 1, Fig. 8(b-2)). Even no additional complete fragmentations are formed, PD reflects the location of potential cracks about to occur through the damage variable (Eqn (7)), which can be observed in Figure 8(b-2). In this scenario, a damage value equal to or greater than 0.5 is designated as indicative of a clear crack path. Conversely, when the damage falls below 0.5, it signifies a weakened area that may serve as a potential site for crack propagation. These demonstrate an advantage of PD in solving failure problems.