A Finite Element Model for Compressive Ice Loads Based on a Mohr-Coulomb Material and the Node Splitting Technique

Ice loads are crucial for marine structures in arctic environments. The Finnish-Swedish Ice Class Rules [1] or the POLAR CODE [2] for ships and the ISO19906 [3] for offshore structures are well established as industrial standards for the design against service ice loads. These standards are to a large extent based on empirical data and idealization. For example, pressure-area considerations play an important role in the ISO19906 and the Polar Code. For unclassified field data, this correlation is undoubtedly correct [4,5]. However, under controlled conditions, the ice pressures seem to depend much more on the relative speed and degree of confinement as on the loaded area [6]. In addition, it is expected that the ice type and temperature are also of great importance for the ice loading.

Shipbuilding industry follows the general trend from experience-based to model and data-driven design processes [7]. For structural applications, the finite element method is mostly applied [8]. However, generally applicable, and extensively validated numerical ice models for fully coupled ice–structure interaction simulations are currently missing [9]. Numerical models are particularly important because full-instrumented and repeatable full-scale ship–ship or ship–ice collision tests are extremely scarce and limited to single scenarios.

A challenge for numerical ice–structure interaction simulations is the multitude of ice models that have been developed for specific purposes. However, these models are typically not transferable to other applications. Therefore, this paper presents a physical-based ice model for brittle and compressive dominated ice loads, based on the Mohr-Coulomb theory.

The mechanical behavior of ice is extremely versatile. The failure mechanism depends strongly on the strain rate. At small strain rates, ice fails in a ductile manner and is dominated by creep. At higher strain rates, ice fails in brittle mode [10–12], which is typical for ship–ice interaction. The transition from ductile failure to brittle failure takes place at strain rates in the order of $∼10−3s−1$⁠.

Compression test allows the determination of physical meaningful strain rates in experiments. However, for real ice–structure interaction scenarios, strain rates are difficult to determine because of the brittleness of the ice. As an indication, the specification of a transition speed is helpful. The transition from ductile to brittle behavior commonly takes place in the velocity range of 0.1 and 3 mm/s [5,6].

The material behavior of ice is characterized by a strong anisotropy of the failure envelope in tension and compression [11]. For low confinement, the failure of ice is limited by Coulombic faults [13]. Under higher confinement, frictional sliding is suppressed on the crack surface. Renshaw et al. [13] describe this type of failure as plastic faulting.

Compressive dominated ice–structure interaction in a brittle failure mode can be analyzed with approaches using continuum and fracture mechanics. In brittle ice–structure interaction, the ice fails by spalling on the free boundaries [14,15]. In the spalled zone, no significant mechanical forces are transferred to the structure. Virtually all mechanical loads are transferred into the structure through a damaged layer, also called High-Pressure Zone (HPZ). In the brittle mode, the pressure distribution of the HPZ varies significantly [6,12,16,17]. The failure processes of crushing and spalling can be identified in Fig. 1.

Three different types of fully coupled finite element models for maritime ice–structure interaction problems are often used in the literature:

1. Global ice load models: The aim of these models is to compute the correct ice-induced drag forces. Such simulations often include large calculation areas and empirical or semi-empirical low fidelity methods. The correct contact pressure distribution is not the objective of these simulation approaches. An example of these types of models is the coarse Cohesive Zone Model of Gürtner [18]. These coarse models are not appropriate to represent ice pressures correctly in the spatial domain on the ship structure.

2. Local load or design load models: The aim of these models is to represent the correct ice forces and pressures in full scale. For this purpose, the ice–interaction process is implemented numerically in a simplified way. Despite the weakness that finite element method (FEM) can only represent discontinuities with difficulty, it is still the most widely used approach for ice–structure interaction simulations.

   The simplest models assume elastic or rigid ice behavior. These models are applicable to simple ice geometries, e.g., a large sphere, where the primary concern is a conservative estimation of the structural deformation. However, the assumption of a rigid or elastic material does not allow to simulate the complex nonlinear material behavior of ice in general, as spalling and crushing are not considered.

   Plastic approaches are another possibility. First plastic ice models are given by, e.g., Ralston [19] and Derradji-Aouat [20]. Application examples for these types of models are the iceberg design load model of Liu et al. [21] or the constitutive model of Ince et al. [22]. Iceberg material models are not applicable to general ice collision scenarios because spalling of the contact interface is not represented.

   Another approach in ice modeling is the application of crushable foam models for the description of ice [23,24]. The crushable foam models achieve reasonable results with respect to global loads. Nevertheless, the transferability of the model and the determination of the input curves are challenging.

   A drawback of many FEM-based ice models is the violation of the principle of mass conservation. The failure of ice is often investigated with the element erosion technique. However, if ice under compression is simulated, element erosion leads to unphysical results because actually the existing ability to transmit loads lacks. Furthermore, the elemental erosion technique violates mass and energy conservation [18,25].

3. Physical correct models: Complex physically based models already exist for selected and specific academic problems. Examples are Representative Volume Models [26,27], Phase Field Models [28], Full Field Modeling [29,30], or in certain cases Cohesive Zone Models [31]. Due to the high computational effort and the high degree of specialization, these models are not yet suitable for the simulation of full-scale ice–structure interaction problems.

As an alternative to the FEM-based approaches, a variety of meshless methods for simulating brittle problems are available in the literature. In this paper, we used the finite-element code LS-Dyna [32] that includes routines for brittle and semi-brittle problems, as well as different meshless methods, e.g., Bounded Discrete Element method (BDEM) [33], Peridynamic [34], Smoothed Particle Galerkin method (SPG) [35], or Smooth Particle Hydrodynamics (SPH) [36]. However, no meshless approach has yet been proven to be generally applicable for brittle problems [37]. According to the author's experience, there are major difficulties in simulating the ice self-contact.

The development of a versatile and robust ice load model, which can represent spalling, material transport, and the continuous nature of the ice failure processes is the subject of this paper. The model is intended to be applicable for the simulation of various crushing-dominated ice–structure interaction problems in the range of low (level ice and ice floes) to high confinement (iceberg collisions) conditions.

The objective is the presentation of an effective finite element model for the full-scale application, which considers cracking in a simplified form. For the validation of the Mohr-Coulomb Nodal Split (MCNS) model, simulated maximum ice forces and contact pressures are compared with experimental results of an ice extrusion test series [6].

Despite a variety of alternative mesh-free methods, FEM is considered the most well-developed choice for simulations involving ice. A big advantage of FEM is many contact algorithms, which allow the coupling of the ice- and structural model. Therefore, the model developed in this paper is implemented into the explicit finite element (FE) solver LS-DYNA R11.1.0.

Kellner et al. [38] collected the results of a large number of uniaxial and triaxle ice compression experiments available in the literature. Brittle tests with fresh or tab water ice in the temperature range of −10 °C + −1 °C are selected in Fig. 2. For the selected data sets, the peak shear stress is plotted against the simultaneously acting peak pressure. The fundamental trend that is shown in Fig. 3 can be represented by a best-fit curve according to the Mohr-Coulomb theory.

The Mohr-Coulomb failure citation is often used for the description of soils, rocks, and concrete materials [39,40]. The brittle failure of rocks and ice are in many respects similar [41]. Furthermore, Mohr-Coulomb is also applicable as a plastic flow condition, for cohesive-frictional materials. The Mohr-Coulomb theory has been already used for ice [42].

Based on Fig. 2, the novel MCNS ice model is developed. The following idealizations are made in the MCNS model with respect to ice–structure interaction problems:

• Creep is neglected.

• A thermodynamic adiabatic process with constant material properties is assumed.

• The effects of the grain structure on the stress distribution are neglected.

• Broken ice remains broken; pressure-induced healing processes of cracks and dynamic recrystallization are neglected.

In the model, an elastic-ideal plastic Mohr-Coulomb Material model is combined with a node splitting approach. During node splitting, the failed element is detached from its neighbors, if the critical plastic failure strain is exceeded. The plastic failure strain is selected as small as possible to achieve a brittle model behavior. Since the elements remain in the model after failure mass, volume, and energy are preserved [43]. For the LS-Dyna keyword deck of the MCNS model, the following cards are used:

• Solver: Explicit finite element solver LS-Dyna R11.1.0.

• Element type: Fully integrated solid elements, ELFORM = −1.

• Mesh: A hexahedron mesh is created as uniformly as possible via splitting of a tetrahedron mesh once. The mesh allows the elements to slide and the fracture paths to be as random as possible. Unphysical stacking of detached elements is prevented. An exemplary mesh of a cylindrical specimen with a diameter of 200 mm and an element size of 12.5 mm is shown in Fig. 3. If elements must pass through a gap, the size of the element should not be larger than half the size of the gap, so that the elements can pass through the constriction.

• Material model: *MAT_MOHR_COULOMB—The parameters for freshwater ice at approximately −10 °C are presented in Table 1. The material parameters are defined with a best fit on the selected data in Fig. 2. The analytical tensile stress is 1.72 MPa. Which is in line with Timco and O’Brien [44]. The corresponding uniaxial compression failure stress is 5.23 MPa.

  To model pressure melting, a maximum pressure according to the equation of state of Feistel [45] is defined with *ADD_EROSION. The melting pressure of −10 °C cold freshwater ice is about 100 MPa. To represent pulverization and extrusion phenomena, as observed during crushing at high confinement [6,46,47], the maximum plastic strain is limited. A critical value of 1 is assumed and implemented with *ADD_EROSION.

• Failure model: The failure model is realized by the simple node splitting approach of the LS-Dyna keyword *CONSTRAINED_TIED_NODES_FAILURE. A critical plastic strain at failure of 0.002 was used for all simulations.

• Ice self-contact: *CONTACT_SINGLE_SURFACE (SOFT = 2, SBOPT = 5, and DEPTH = 5) with a constant friction coefficient of 0.05 is utilized.

• Time-step size: The critical time-step size is reduced precautionary with TSSFAC = 0.5.

• In general: The best practice recommendations of Refs. [48,49] are used.

For further information regarding the LS-Dyna inputs is referred to the LS-Dyna documentation [32,50].

The parameters of the model are derived sequentially as presented in Fig. 4. The sequence was chosen so that the parameters could be determined as unique as possible.

Contact models have a crucial function for the proposed model. During the collision simulation, ice elements become detached if the equivalent plastic strain exceeds 0.2%. For these elements, it is necessary to consider contact, and a proper contact formulation is required to transfer mechanical loads.

To simulate the contact problems in an optimal way, only surface-based formulations are used. The ice self-contact (contact between broken ice elements) is represented by a single surface contact SOFT = 2 with an edge-to-edge treatment. The surface-based SOFT = 2 contact is recommended for complex problems in LS-Dyna [48]. During the first MCNS ice–structure interaction simulations with different element sizes, it became apparent that the contact stiffness has a major influence on the stability of the model. A too high contact stiffness leads to unstable results.

For the investigation of the contact stiffness, numerical tests were conducted. A cylindrical ice model, according to Fig. 4 is considered. The applied boundary conditions are as follows: (i) a prescribed displacement at the top of the cylinder, and (ii) a fixed boundary condition in normal direction for the remaining boundary surfaces. The simulations consider two extreme cases. At first, a conventionally assembled finite element mesh without contact is used. Second, a mesh where all elements are detached is studied. In this second case, the forces of adjacent elements are transferred only via the contact algorithm.

As expected, the stiffness of the model decreases in the second case considering the detached elements. For the element size of 12.5 mm, the stiffness reduction is 27.6%. This effect is reasonable for ice. During experiments, broken ice is softer than intact ice, due to the cracks that develop on the micro and macro scale [51–53].

The recommended contact stiffness scale factors are given in Table 2. With larger elements, the simulations initially were unstable. Therefore, the contact stiffness for the element size of 50 mm has to be reduced to 0.75. With the adaption of the contact stiffness scale factor, the stiffness reduction of the broken ice changes only slightly in comparison to the element size of 12.5 mm.

This section presents the results of the MCNS model based on the parameter set according to Table 1, Table 2. It is noted that the models of the uniaxial tension and compression test as well as the ice extrusion test were also used for parameter identification based on the process in Fig. 4.

In order to setup and test the MCNS model in well-defined applications, uniaxial tensile and compression tests was simulated. A cylindrical specimen (diameter = 200 mm and length = 200 mm) similar to Fig. 3 is considered. The element size is 12.5 mm. In the simulation, a fixed boundary condition for the bottom in the axial direction is applied. The top surface of the model is subjected to a controlled translational velocity of 10 mm/s. The lateral boundaries are all unconstrained.

The simulation results of the model are given in Fig. 5. The corresponding stress–strain curves are given in Fig. 6. In particular, the stress shown in Fig. 6 is calculated by determining the total reaction force of the top nodes divided by the constant top surface area. As expected, the stress–strain curves depend on the loading direction due to the used Mohr-Coulomb yield condition. Figure 6 shows that response of the overall structure can be attributed to the ideal plastic material behavior. The constant plateau according to Fig. 6 corresponds to the compressive failure stress σ_c and tensile failure stress σ_t of the used plasticity model. In tension, the simulation achieves a maximum stress of 1.69 MPa in comparison to analytical tensile failure stress of 1.73 MPa according to Eq. (2). The corresponding deviation is 2.3%. In the compressive case, the simulation reaches 5.12 MPa. The analytical value of Eq. (3) is 5.21 MPa. The deviation for the compressive case is specified as 1.73%.

According to Fig. 6, depending on the loading direction, the post-fracture behavior ranges from virtually no transmitted forces in the case of tension to the transmission of forces during compression.

Finally, this section addresses the post-fracture behavior of the used Mohr-Coulomb material model in case of compressive failure. For the study, an idealized failure was simulated of the previously considered compression test, in which at 0.08 s, the upper boundary condition was deleted. In Fig. 7, the z-stresses (in the axial direction) of the MAT173 material model used in the MCNS Model and the common elastic material model MAT1 are presented. In the elastic case, tensile and compressive oscillations occur with an amplitude equal to the previous compressive stress. In the case of the Mohr-Coulomb model, vibration energy dissipates into plastic deformation and tensile vibration therefore only reach values in the range of maximum tensile strength. This characteristic is important for the MCNS model to obtain a stable numerical behavior under compressive failure and considering the anisotropic behavior of ice.

Subsequently, the ice extrusion tests with specimens having a diameter of 100 mm and 200 mm are simulated. In the ice extrusion test, a cylindrical ice specimen is pushed out of a confining pipe against a test structure.

All experiments were performed with the same type of ice. The ice was isotropic fine-grained freshwater ice. The average temperature of the ice specimens during the experiment was −10 °C, and a grainsize of approximately 2.6 mm was determined. The test setup of the ice extrusion test is shown in Fig. 8. During the tests, the specimen diameter, the speed of the cylinder, the gap height between the pipe and the specimen, and the cone angle are varied. For more information regarding the experiments, the reader is referred to Herrnring et al. [6]. The experimental results of the novel ice extrusion test are especially suited for the development of numerical models since it directly represents a well-documented continuous ice–structure interaction process.

A comparison of the experimental and numerical results of the MCNS model for two different gap heights of two brittle ice extrusion tests with a specimen diameter of 200 mm and a cone angle of 20 deg against a quasi-ridged structure is given in Fig. 9. An assumed and constant coefficient of friction of 0.03 is used for the contact between the structure and the ice. The small coefficient of friction is chosen according to Timco and Weeks [54] since the steel structure was covered with a very smooth Kapton HN500 foil, which protects the TekScan sensor.

According to Fig. 10, the MCNS model can reproduce satisfactorily the maximum force of both experiments. The increase from approximately 60 kN to over 160 kN can be attributed to the reduction of the gap height from 50 mm (Fig. 9(a)) to 25 mm (Fig. 9(b)). These values correspond well to the experimental peak forces. Due to the reduction of the gap height, the confinement in the crushing zone is increased [6]. In other words, it becomes difficult to push the ice out of the interaction area within the pipe diameter. The force increases accordingly.

A significant force peak at the beginning characterizes the force curve of the 50 mm experiment. Large spalls are extruded when the maximum force is reached. Until the spalling progress starts, the pressure distribution of the simulation is circular shaped.

As shown in Fig. 9(b), the ice behavior is different for both gap heights. In the case of the smaller gap height, the maximum forces occur at the beginning of the extrusion phase. This effect is very well represented by the MCNS model in comparison to the experimental results. The loads depend in particular on the selected friction coefficient for the ice–ice contact. Due to the high confinement, the ice elements were more plastically deformed, and a significant part of the elements is eroded according to the pulverization criteria.

The contact pressures on the structure corresponding to the maximal reached reaction force are shown in Fig. 9. The 25-mm simulation shows a much larger contact area compared to the 50-mm simulation. A distinct HPZ can be seen in the middle of the contact area.

Moreover, a detailed evaluation of the loaded contact area for the same simulations is given in Fig. 10. In the plot, the measured and simulated process–area as well as force–area curves are presented. The measured data are collected with a TekScan pressure measurement system [6]. In this evaluation, the TekScan data are only used to determine the real transient contact area. For the following comparison, the resolution of the TekScan measurement was coarsened according to the structural FE model during post-processing.

The MCNS model reflects the increase of the contact area with increasing confinement. Thus, the loaded area increases from about 60% for a gap of 50 mm to 95% for a gap of 25 mm. Nevertheless, the loaded area is overestimated compared to the experimental results by the MCNS model by about 33% in both cases. Moreover, it can be observed that the maximum forces of the simulation occur at comparable normalized loaded areas. The pressures depend strongly on the selected element size of the structural mesh, which is however subject to further investigations.

A comparison of the maximum nominal pressures for extrusion tests and the corresponding simulations with 100 mm and 200 mm diameter and different gap heights is given in Fig. 11. For comparison, the nominal gap height, which is defined as gap height divided by specimen diameter (G/D) is used. The simulation results of the 100 mm and 200 mm series are in very good agreement with the experimental results. The results of the simulations reflect significantly increasing load capacity in line with increased confinement.

Finally, the influence of the cone angle on the ice forces was investigated. To enable fracturing due to low confinement, a large gap height of 100 mm was chosen for the study. The experimental and numerical results are compared in Fig. 12. Except for the simulation with a cone angle of 30 deg, the simulations match the failure forces well.

To show the applicability and transferability to larger and energy-limited problems, a double pendulum experiment of Gagnon et al. [16] was simulated. During the test, two counter-rotating pendulums of equal speed and weight collide. An ice sample is attached to one of the pendulums. The second pendulum is equipped with a flat acrylic plate, which is part of a pressure measuring device. The cone diameter was 1 m and the cone angle 30 deg. The effective impact mass of both colliding bodies is given at approximately 4330 kg each. The impact speed of the simulated test “May22_2014” was given with 4.1 m/s.

A larger element size of max. 50 mm is used. To achieve a stable ice self-contact, the SFS of the *CONTACT_SINGLE_SURFACE is reduced according to Table 2 to 0.75. Apart from that, the model is adopted unchanged.

As observed in the experiments, the sample fails crushing dominated. Large ice pieces were not detached. A plot of the results for t = 0.025 s is given in Fig. 13. To compare the experimental force–time curve with the simulation results, the contact reaction force is computed. The simulated and measured curves are given in Fig. 14. Both the force level as the qualitative behavior was reproduced by the MCNS model.

The MCNS model was able to simulate crushing-dominated ice loads. Forces were well represented with the same material parameters for various tests. The method provided valid results for force- and energy-limited problems shown in this study. One drawback for low confined processes of the model is that the forces of the second force peak are too high. This phenomenon was observed for example during the ice extrusion test as shown in Fig. 9(a). In addition, the maximum contact areas in the test cases shown are 33% larger than those in the experiments. A possible explanation could be the relatively coarse mesh of the ice model used or the simplified mechanical continuums behavior.

Most model parameters have a physical meaning and are associated with specific physical processes. Spalling, continuum mechanical deformation of the detached ice elements, pressure melting, the ice self-contact, and pulverization were considered in the MCNS model. All these processes are essential for the ice–structure interaction process. It must be clearly stated that the description of the continuum behavior of the crushed material, the friction coefficient for the ice–ice contact, and pulverization have hardly been provided in the literature. Accordingly, the parameters given in the paper provide a basis for further modeling and analyses.

As the crack tip is not represented directly in the MCNS model due to the node splitting technic and the plastic failure citation, which is mostly implemented for numerical reasons, for simulations of tensile-dominated single crack-driven problems, other methods such as cohesive zone method (CZM) are better suited [31,55]. However, these methods are difficult to apply for arbitrary and interacting crack paths for compressive 3D problems [56].

Due to the complex contact conditions, the MCNS model is numerically more expansive than, conventional FE models. Nevertheless, the time-step sizes are within a normal range. The time-step of the small-scale simulations with an edge length of the elements of 12.5 mm was 2.45E-07 s, and in the case of the double pendulum test, the time-step with an element size of 50 mm was 1.28E-06 s. With an appropriate element size, the application of the MCNS model for ship–ice interactions seems feasible.

The MCNS ice model using Mohr-Coulomb plasticity and a simple node splitting failure technique was successfully implemented and validated for brittle ice–structure interaction problems. For a variety of complex crushing-dominated problems, such as the ice extrusion test and the double pendulum test, the numerical results are in particularly good agreement with the observed values.

The number of required model parameters is relatively small, and most parameters have a direct physical meaning. Input parameters of the ice model were partly determined using experimental databases or were calibrated during simulations of the ice extrusion test. The model enables geometric change due to fragmentation, material transport and is at low confinement rations almost energy and mass-preserving. First evaluations of the contact pressures showed slightly too large contact areas. However, basic effects and pressure values are reproduced comprehensibly.

The interested reader may contact the authors for access to exemplary input files of the MCNS model.

We highly appreciate the trustworthy and inspiring collaboration and support of Dr. Paul Hess of ONR and Douglas Lesar.

This work was supported by the US Office of Naval Research Global (ONRG) under NICOP Grant No. N62909-18-1-2127. We also acknowledge the financial support of the Lloyd’s Register Foundation in the research project “Recommended practice of scenario-based risk management for polar waters.” Lloyd’s Register Foundation helps to protect life and property by supporting engineering-related education, public engagement, and the application of research. It is stated that all funders are not responsible for any of the content of this publication.