Numerical Analysis of Ice-Structure Impact: Validating Material Models and Yield Criteria for Enhanced Prediction of Impact Pressure

1. Introduction

Vessels and offshore structures operating in Arctic regions constantly face the risk of structural damage from collisions with ice. A notable incident took place in January 1994, when the tanker Overseas Ohio, weighing 10,000 tons and traveling at approximately 5 m/s, sustained bow damage due to an ice collision off the coast of Alaska, resulting in severe marine pollution [1]. In the design phase of Arctic ships and offshore structures, predicting ice collision force is crucial for identifying structures with impact-resistant stiffness to minimize the risk of hull damage. This study introduces a numerical procedure and methodology for estimating ice collision force and compares the results with experimental findings to assess the validity of the load prediction process. The determination of collision forces can be categorized into two main analytical frameworks: dynamic structural analysis and fluid-structure interaction analysis. Several studies [2,3,4] have employed dynamic structural analysis to calculate ice collision force. Additionally, methods utilizing fluid-structure interaction to simulate ice-structure collisions have been explored [5], along with dynamic structural analyses focusing on the damage to both ice and structures [7,8]. Various constitutive equations and damage models for ice have been applied to predict collision force [3,9,10]. Despite the numerous approaches to ice collision force analysis discussed in the literature, this topic remains a subject of ongoing investigation. This continued interest is due to the significant variability in material properties, including elastic and plastic deformation and ice fracture, as well as the irregular changes in contact area during collisions, which pose significant challenges in accurately simulating ice damage and structural response. Different models have been applied to characterize the interactions between ice and other materials, including the Drucker-Prager model [11,12] and the crushable foam model [7,8]. However, variations in the ice-induced resistance faced by ships have been observed, highlighting the need to evaluate load predictions based on different yield criterion assumptions. Furthermore, studies using the discrete element method (DEM) have simultaneously considered the fragmentation of level ice during ship collisions, the formation and movement of ice blocks, and the interaction between ice blocks and ship structures [13,14,15]. Importantly, the DEM accounts for both ship and ice motion, including contact loading between individual ice blocks and the ship’s hull. Moreover, this approach represents the shape of ice as polygons or polyhedra and can simulate the mutual motion between ice and ship structure, as well as the ice fragmentation process. Models evaluating the impact of local ice collision force on ship structures, including various material properties and ice shapes, continue to attract considerable attention. Previous studies have conducted ice collision tests using various impact apparatuses. For example, Zhu et al. [7] and Cai et al. [8] conducted collision experiments with wedge-shaped ice specimens at speeds of 4 to 7 knots to measure the loads applied to the impacted bodies. Additionally, Jang et al. [16] presented the results of collision experiments between hemispherical ice specimens and steel plates using a collision pendulum at controlled speeds (5, 7, 10 knots) in a low-temperature chamber. These experiments accurately predict loads under experimental conditions, but they require scaling down the size of the ship and ice, and the experimental conditions are limited. In contrast, numerical simulation models do not account for size and shape limitations for ships and ice or collision speeds, but accuracy verification is essential. Therefore, numerical simulations must be performed to predict the ice collision force and to conduct analytical and experimental work for accuracy verification. The purpose of this study is to propose a procedure for numerically evaluating the local ice collision force and the response of steel plates by simulating simplified experiments. This procedure considers various factors, including yield criteria, material behavior, and strain rate effects in the material model of ice. To achieve this, several numerical simulations were conducted to reproduce the local collision force and structural deformation at the collision speeds realized in experiments. These simulations vary the material model to propose results and evaluate the impact of material properties on ice collision force predictions. This study also sought to expand the range of environmental variables associated with ice collision, (e.g., collision speed) in existing empirical design formulas to demonstrate the applicability of numerical analysis for load prediction. By systematically examining the performance of proposed material models against experimentally observed responses, this research aims to develop the appropriate numerical model for ice collision force prediction.

2. Material Properties of Ice

The strength and behavior of ice are strongly influenced by its crystal structure, salinity, temperature, and strain rate [17,18]. The determination of material properties and constants associated with ice deformation is uniquely challenging due to the complexity of the variables that influence this phenomenon. Numerous studies [3,5,9,19,20] have focused on the mechanical properties of ice, such as its constitutive relationships, fracture properties, and strain rate dependence, rather than the subtleties of its microstructure or thermal variations. Building on these findings, the present study incorporates the elastoplastic behavior, fracture characteristics, and strain rate sensitivity of ice. The referenced literature describes the mechanical properties of ice in detail, which are used as the basis for deriving the material properties and constitutive equations listed in Table 1. Therefore, this study sought to apply these data in an analytical framework where the results of different combinations of material properties can be compared and analyzed.

2.1. Constitutive Equations for Ice Elasto-Plastic Behavior

In the initial phase of this study, two different constitutive models were evaluated to characterize the elastic-plastic behavior of ice. The first constitutive model considered herein was the linear elastic-perfectly plastic with damage model (LEPP), dominated by elastic behavior. The second material model was the linear elastic-plastic with damage model (LEP), assuming plastic deformation as the primary response. The load-displacement curves for these models are depicted in Figure 1, illustrating a linear increase in behavior with applied load.

The LEPP model, which emphasizes the elastic domain, suggests that when the load is removed, the strain and stress are fully recoverable, indicating that no permanent deformation occurs until the yield strength (${\sigma }_y$)) is reached. At the yield point, the material begins to suffer damage (Damage, D=0). This constitutes a crucial aspect of the model, as it marks the onset of fracture or failure in the ice structure. In contrast, the LEP model assumes a low yield stress threshold. Therefore, yielding occurs immediately upon loading, resulting in a negligible elastic phase. Consequently, after yielding, the ice does not return to its original state when the load is removed but continues to suffer damage until it reaches its ultimate strength (${\sigma }_y$), at which point damage propagation becomes more pronounced. The adequacy of these two material models is discussed by comparing their predicted behavior with experimental observations. Finally, this study aims to compare the predictive capabilities of these two material models with experimental results to determine the most appropriate description of ice behavior.

2.2. Yield Criteria of Ice

(1)

Crushable foam model

The crushable foam model effectively represents the hardening characteristics associated with the volumetric strain of ice under compressive loading. This model has been successfully applied by researchers such as Gagnon & Derradji-Aouat [21], Kim et al. [4], Gagnon [3], and Gao et al. [22], and was also utilized by Obisesan & Sriramula [23] in their analyses. A key feature of the crushable foam model is the assumption that the material does not exhibit elastic recovery after yielding, resulting in permanent deformation. This model not only accounts for the failure stress associated with tension and compression but also considers the hardening characteristics. The yield surface ($F$) of the crushable foam model is defined by Equations (1) to (3) and is represented as an elliptical shape, as shown in Figure 2. Due to these characteristics, yielding occurs under the influence of hydrostatic compression pressures.

$$\mathit{F}=\sqrt{{\mathit{q}}^2+{\mathit{\alpha }}^2{\left(\mathit{p}-{\mathit{p}}_0\right)}^2}-\mathit{B}=0$$

(1)

$$A=\frac{\left|p_c\right|+p_{\mathrm{t}}}{2}, B=q_{\mathrm{m}\mathrm{a}\mathrm{x}}, p_0=\frac{p_c-p_{\mathrm{t}}}{2} , \alpha =\frac{B}{A}=\frac{3k}{\sqrt{\left(3k_t+k\right)\left(3-k\right)}}$$

(2)

$$k=\frac{{\sigma }_c^0}{p_c^0}, k_t=\frac{p_c}{p_c^0}$$

(3)

where $A$ and $B$ represent the major and minor axes of the elliptical yield surface in the crushable foam model, respectively, with $p_0$ denoting the center. $p_c$ and $p_{\mathrm{t}}$ denote the yield strengths under hydrostatic compression and hydrostatic tension conditions, respectively. The factor $k$ represents the ratio of yield in compression, whereas $k_t$ represents the ratio of yield under hydrostatic conditions. The initial yield stress in hydrostatic conditions is indicated by $p_c^0$, whereas ${\sigma }_c^0$ represents the uniaxial compressive yield stress, which is shown in Figure 2 as the point where the yield surface intersects the $q$-axis. As the deformation of the ice enters the plastic regime, the size and shape of the yield surface change; this change can be expressed by the relationship between the uniaxial compressive yield stress and the volumetric strain. These changes are described by a strain hardening model where the yield surface varies based on a line connecting $p_c$ and $p_{\mathrm{t}}$ intersecting with the $q$axis. The evolution of the yield surface can be represented as a function of the volumetric plastic hardening, $p_c\left({\epsilon }_{vol}^{pl}\right)$, as expressed in Equation (4), which expresses the relationship between the volumetric plasticity and uniaxial stress. Volumetric hardening assumes that uniaxial plastic deformation (${\epsilon }_{axial}^{pl}$) is equivalent to volumetric deformation (${\epsilon }_{vol}^{pl}$). To implement the hardening of the crushable foam plastic model in ABAQUS, the change in the yield surface is reflected by the uniaxial compressive yield stress and plastic strain rate.

$$p_c\left({\epsilon }_{vol}^{pl}\right)=\frac{{\sigma }_c\left({\epsilon }_{axial}^{pl}\right)\left[{\sigma }_c\left({\epsilon }_{axial}^{pl}\right)\left(\frac{1}{{\alpha }^2}+\frac{1}{9}\right)+\frac{p_t}{3}\right]}{p_t+\frac{{\sigma }_c\left({\epsilon }_{axial}^{pl}\right)}{3}}$$

(4)

(2)

Drucker-Prager model

The Drucker-Prager model is known to be suitable for materials with internal friction, such as rock and concrete. This model effectively captures the strengthening of certain materials under compressive loads and is appropriate for describing the behavior of materials where the compressive strength is higher than the tensile strength. Noting the similarity between the behavior of ice and concrete, Kajaste-Rudnitski & Kujala [11], as well as Zhang et al. [12], have applied the Drucker-Prager approach to ice collision models. A key feature of the Drucker-Prager model is its use of a non-circular yield surface to adjust the yield stress values in both tensile and compressive scenarios. This model allows for the independent setting of the dilation and friction angles, and as shown in Figure 3, it defines the yield surface through three invariants. The yield surface ($F$) of the Drucker-Prager model can be expressed as described in Equations (5), (6), and (7).

$$F=t-p\bullet tan\beta -d=0$$

(5)

$$d=\left(1-\frac{1}{3}tan\beta \right){\sigma }_c$$

(6)

$$t=\frac{1}{2}\left[1+\frac{1}{K}-(1-\frac{1}{K}){\left(\frac{r}{q}\right)}^3\right]$$

(7)

where $\beta$ represents the material’s angle of friction in the $p$-$t$ plane, $d$ denotes the material’s cohesion, $K$ is the ratio of the triaxial tensile yield stress to the triaxial compressive yield stress, $p$ is the equivalent pressure, and $q$ represents the Von-Mises stress. In other words, the yield condition of the Drucker-Prager model typically takes on a non-circular form, with the yield surface determined by the material’s angle of friction ($\beta$), cohesion ($d$), equivalent pressure ($p$), and Von-Mises stress ($q$). The study by Kajaste-Rudnitski & Kujala [11] provides a detailed application of such a model.

(3)

Ductile damage model

Material failure or damage is defined as the gradual deterioration of a material’s stiffness, resulting in a complete loss of load-carrying capacity. To address this behavior in ice, a ductile damage model was adopted in this study. Fracture-based theory illustrates the stress-strain behavior of the material, as depicted in Figure 4.

• Damage Initiation (D=0): The fracture model identifies the commencement of damage when the material reaches its ultimate strength. Damage initiation is typically associated with the equivalent plastic strain (${\epsilon }_{eqv}$).
• Damage Evolution (0<D<1): In Figure 4, the dotted line represents the undamaged response, indicating the material’s pure plastic behavior without the influence of a damage model. The solid line, labeled as “softening,” depicts the behavior of the material once it has sustained damage. Following damage initiation (D=0), the model simulates the progression of damage. As the damage evolves, the material’s stiffness starts to decrease, and the damage variable (D) increases.
• Fracture (D=1): When the damage variable (D) reaches 1, signifying the completion of the damage progression, the material’s internal cohesive strength is completely lost. At this point, material fracture occurs, and bonded elements may separate or be removed (i.e., element deletion).

2.3. Strain Rate Dependency of Ice: Yield Strength Ratio

This section examines the complex relationship between yield strength ratio and strain rate. Notably, higher strain rates typically lead to an increase in the yield strength of ice, as highlighted by Tippmann [24]. Moreover, Barette and Jordaan [17] reported that lower temperatures also contribute to an increase in the intrinsic strength of ice. It is also important to consider that both strain rate and yield stress tend to increase simultaneously as collision energy increases. Tippmann [24] explored these dynamics by investigating interactions between ice and structural elements, particularly in scenarios involving spherical ice specimens colliding with carbon fiber-reinforced plastic (CFRP) plates. Furthermore, in the context of high-velocity collisions, it is important to account for the strain rate dependence of the materials involved. This aspect has been addressed by adopting the model proposed by Tippmann [24]. Notably, the range of strain rates outlined by Tippmann [24] extends beyond the strain rate of this study, enabling the analysis of a broader spectrum of collision velocities. Based on observations from collision experiments with ice, it is evident that yield strength is dependent on strain rate. Tippmann [24] and Tippmann et al. [25] conducted studies focusing on this dependency, analyzing the relationship between strain rate and yield strength. Their findings demonstrated that as the strain rate increases, the compressive strength also shows an upward trend. During the plastic deformation of ice, yield occurs at a specific strain rate, after which the stress remains relatively constant. This phenomenon is particularly observed during high-speed ice collisions, making a perfect plasticity model appropriate, assuming a yield hardening coefficient of 0 [26,27]. Tippmann [24] and Tippmann et al. [25] compiled various research data and represented the relationship between strain rate and yield strength using a linear-logarithmic curve. In this study, this relationship, specifically the ‘strain rate-yield ratio,’ was incorporated into the analysis. Figure 5 illustrates the relationship between compressive plastic strength and the yield strength ration within a strain rate range of 0.1/s to 106/s, providing a visual representation of this curve and the relevant values. The yield strength ratio is normalized to the initial compressive plastic strength of 5.2 MPa.

3. Ice material Calibration

In the previous section, various material properties of ice were described, focusing on plastic deformation relationships, yield conditions, and the strain rate dependency of yield strength. As reported in previous studies, different models are employed for different material properties, reflecting the complex and uncertain nature of the deformation properties of ice. In this section, considering the diversity of material models mentioned in Section 2, the selection process for material models and material constants that can account for the dynamic compressive deformation characteristics of ice is presented.

• Plastic deformation relationships: Models such as the linear elastic-perfectly plastic (LEPP) and linear elastic-plastic (LEP) models define how materials behave when they exceed the elastic limit. While the LEPP model assumes that no deformation occurs beyond the yield point, the LEP model describes post-yield hardening or softening. These material plastic properties are applied to the analyses to verify their results.
• Yield criteria: The crushable foam (CF) and Drucker-Prager (DP) models define the conditions under which materials reach yield. The CF model simulates the collapse of cellular materials under compression, whereas the DP model models the yield of pressure-dependent materials. Given the sensitivity of ice to pressure and density, the appropriateness of these conditions was investigated.
• Strain rate dependence: In dynamic events such as collisions, strain rate affects the yield stress and failure modes of the material. Results are analyzed based on whether or not strain rate dependence is considered.

In this study, the combination of these three models was applied to material modeling efforts and compared with experimental data to select the most appropriate model. This approach helps to reduce the complexity and uncertainty associated with material modeling and allows for more accurate predictions of ice behavior in real-world scenarios. In this study, four models were adopted by combining (1) two plastic constitutive equations and (2) two yield conditions. Additionally, (3) the effects of strain rate are included in all eight material model combinations, as summarized in Table 2. Determining material constants for the ICE material models, including elasticity coefficients, yield stresses, hardening coefficients, and yield strength ratios, typically requires performing tensile tests, bending tests, or triaxial compression tests. In some cases, these constants can be indirectly estimated through theoretical analysis. In this research, however, these values were derived by performing quasi-static compression experiments on ice and simulating these experiments using finite element analysis. A detailed description of this process is provided in Section 3.1 and Section 3.2. Additionally, the compression strength ratio as a function of strain rate, as measured by Tippmann [24] in high-speed collisions with CFRP, was used to assess the dynamic deformation properties of ice with strain rate dependence. The values required to define the strain rate dependence are provided in Figure 5 of the previous section, with the lower limit of the compression strength ratio curve being applied.

3.1. Determination of Material Model Constants

This section discusses the process of determining the material constants required for the elastic constitutive equations (i.e., linear elastic deformation relations) for ice, including the LEPP and LEP models, as well as the yield conditions required for the CF and DP models. The constants that must be defined for each material model are summarized below:

• LEPP model: This model requires the determination of two fundamental material constants, Young’s modulus ($E$) and the yield stress (${\sigma }_Y$). Young’s modulus is conventionally obtained from tensile tests, whereas the yield stress is the stress point at which the material begins to undergo plastic deformation.
• LEP model: Similar to LEPP, the LEP model requires the determination of Young’s modulus ($E$) and yield stress (${\sigma }_Y$). In addition, it requires consideration of additional constants associated with the hardening model and the hardening coefficient or plastic flow rule.
• CF model: This model requires the definition of material constants that govern the compressive behavior of foams. Key constants include density, initial strength (${\sigma }_i$), and hardening coefficient. The determination of these constants is facilitated by the analysis of stress-strain curves derived from compression tests on specimens.
• DP model: The DP model requires a set of constants related to the angle of internal friction (β), material cohesion (d), drag parameters such as the dilatancy angle (ψ), and constants related to the hardening rule. Typically, these constants are obtained from triaxial compression tests and should be considered in conjunction with stress paths.

For material constants that are difficult to measure or calculate directly, assumptions were made based on existing literature. Specifically, values from Tippmann [24], Kajaste-Rudnitski & Kujala [11], Zhang et al. [12], Kim et al. [4], and Han et al. [9] were used to set the density and Poisson’s ratio of ice to 900 kg/m³ and 0.3, respectively. Undetermined material constants for the CF and DP models were hypothetically varied and then incorporated into finite element analyses, ultimately yielding material constants consistent with experimental results. It should be noted, however, that the material models applied to quasi-static compression simulations, which primarily represent low rates of deformation, are tailored to exhibit characteristics specific to quasi-static compression. The constants that define the yield criteria, including parameters associated with damage initiation and stiffness reduction after initiation, are taken from the works of Obisesan & Sriramula [23] and Jeon & Kim [20]. The determined parameters for the yielding models are summarized in Table 3. Additionally, the undetermined material constants specific to each material model were refined to best match the experimental results through simulation analyses simulating ice compression tests, as described in Section 3.2. The material constants defined in the ice compression simulations are as follows:

• LEPP model: The apparent elastic modulus was assumed to be 0.14 GPa. The yield stresses for the CF and DP models were set at 4.8 MPa and 3 MPa respectively. The material damage model parameters included a failure strain of 1.0 × 10^-4 and an energy release rate of 15.0 J.
• LEP model: The elastic modulus, influenced by stress and strain rate, was set to 9 GPa. Hardening constants for this model are detailed in Table 4 for both the CF and DP models, and failure strains were designated as 0.019 and 0.025, respectively, with a damage evolution energy of 15.0 J.

These material constants were carefully tuned to accurately simulate the ice compression behavior under the chosen material models.

3.2. Finite Element Analysis Simulating Ice Compression Experiments

This study successfully conducted quasi-static compression tests on ice specimens using compression testing equipment, which enabled the validation of specific material constants. Building on these experiments, finite element analysis (FEA) was then performed. Throughout this process, multiple FEA iterations were undertaken to determine previously undefined material constants. According to the quasi-static compression test results reported by Jang et al. [16], it was observed that ice specimens, configured as square columns with cross-sectional dimensions of 100 mm × 100 mm and a height of 250 mm, failed under loads ranging from approximately 25 to 32 kN when compressed by 5 mm (Figure 6).

For the finite element modeling, these specimens were represented using 39,375 solid elements (C3D8R), with each element having a grid size of 4 mm. The compression plate in the model was assumed to be a rigid body (Figure 7). The setup for the analysis model was as follows:

• Contact conditions: The lower surface of the ice specimen was fixed in place. A hard node-to-surface contact condition was applied to the upper surface, which came into contact with the rigid compression plate. This setup was crucial to prevent any penetration between the specimen and the plate.
• Compression conditions: The rigid compression plate was set to move at a speed of 1.0 mm/s, mirroring the experimental conditions. This movement applied a controlled displacement to the specimen, thereby replicating the actual test conditions in the numerical model.

The quasi-static compression experiments on ice were reproduced using finite element analysis models to calculate the load and displacement values leading to failure. Verification was achieved through a comparison with experimental results. One of the primary objectives of this study was to optimize the material constants for each finite element model to accurately reproduce the average failure load of 28 kN, as observed in the experiments. This involved finely tuning various material constants, including yield stress, hardening coefficient, and parameters of the damage model. The aim was to ensure that the predicted load-displacement relationship derived from the FEA closely matched the experimental data, maintaining an error margin of less than 1.5%.

In the experimental setup, the reaction force corresponding to the displacement was measured. In contrast, the total compression load acting on the bottom surface of the ice was calculated via FEA. Through iterative processes, the material constants were precisely adjusted to align with experimental observations. The refined load-displacement relationships for each material model are depicted in Figure 8. Table 5 summarizes the compression loads at specific points, confirming that the material constants deduced in this study effectively replicate the ice compression behavior observed in the experiments.

4. Simulation Model for the Ice-Structure Collision Experiment

In this study, the analysis of ice-structure collisions is divided into two aspects: ice damage and the dynamic response of the structure. Explicit finite element analysis was conducted to investigate these aspects. The commercial software ABAQUS/Explicit was used to simulate the nonlinear dynamic response during the collision process between the two objects. Solid elements (C3D8R) were chosen to represent the yielding and failure behavior of the ice, whereas the structure was represented by shell elements (S4R). The main features of this finite element analysis model are explained below:

• Ice yield criteria: The CF and DP models were used to model the yielding and failure behavior of ice, as described in Section 2.1.
• Contact conditions: Contact conditions were established to account for the interaction between the two objects and to calculate the load transferred from the ice to the structure. For the contact between ice elements and plates, hard contact and tie conditions are applied.
• Element with reduced integration: Reduced integration was used to establish the stiffness matrix to enhance the computational efficiency, along with the selection of the hourglass control.
• Time increment: In the explicit finite element analysis, the determination of the time increment (△t) is crucial for achieving both accuracy and computational efficiency. In our study, the time increment for explicit integration was calculated using the following formula: $△t=L_{min}/c_d$, where $L_{min}$ represents the size of the smallest element within mesh. Moreover, $c_d$ represents the material’s dilatational wave speed, which is a function of the material’s mechanical properties and is calculated using the formula: $\sqrt{E/\rho }$, where E denotes the elastic modulus of the material, and ρ denotes the density.

Both quasi-static material models and dynamic material models, which account for strain rate effects were discussed in the previous section, were applied in the analysis. Analyses were conducted for each yield model considered, and the obtained strain and ice load on the plate were compared with the experimental values. This study encompasses four distinct material models: the LEPP, LEP, CF, and DP models. These models are integrated into the FEA to assess and compare the fidelity of ice collision simulations across each material model. This comparative analysis aims to identify the most accurate model in simulating the behavior of ice-structure impact.

4.1. Overview of the Ice-Structural Collision Simulation Model

In this study, the collision pendulum experiment conducted by Jang et al. [16] was chosen for numerical simulation. A detailed discussion on the material models associated this simulation is provided in Section 2, Section 3.1, and Section 3.2. Figure 9 (left) shows the ice specimen, the impactor, and the specimen holder used in the experiment. Correspondingly, the finite element model (FEM) mirrors the geometry of the ice and steel plate specimens, as well as the specimen holder. Figure 9 (right) illustrates the analysis model, which incorporates the experimental apparatus. The dynamic behavior of collision and the steel plate’s response are simulated using an explicit scheme of FEM. The mesh of the FEM is presented on the Figure 9 (b). The collision velocities tested in the experiment, specifically 2.5, 3.5, and 5 knots, were applied to the FEA, which spans a duration of 0.05 seconds.

Figure 10 and Figure 11 illustrate the configurations of the ice specimen and holder, as well as the steel plate and its corresponding jig, respectively. The ice specimen is modeled as a hemisphere with a 250 mm height and a 300 mm diameter, represented by solid elements in the simulation. The total weight of the ice holder, including the pendulum’s mass, was defined as 150 kg. The steel plate, with dimensions of 1,000 × 500 mm and a thickness of 20 mm, was modeled using shell elements. Additionally, the jig for supporting the steel plate was assumed as a rigid body, measuring 200 × 500 × 50 mm. At the interfaces where the ice specimen contacts the holder and the steel plate contacts the jig, node-to-surface contact with hard contact and frictionless conditions were implemented. In these contact conditions, the steel plate specimen, which possesses higher rigidity, was designated as the master surface. Additionally, tie contact conditions were established between the ice specimen and its holder, and between the steel plate specimen and its jig. Table 6 summarizes the shapes and conditions of each component employed in the analysis. A critical aspect of our analysis involved assessing mesh convergence. During the analysis process, mesh convergence was examined, and the effect of mesh size on the collision force was investigated by varying the ice grid size from 2 to 7 mm. The results indicated that a biased mesh size of 2–5 mm was the most appropriate.

4.2. Collision Analysis Reflecting Quasi-Static Compressed Material Model

The four material models introduced in Section 3.1 were applied to the analyses for specific collision velocities. However, this section did not include the dynamic strain rate effect of ice. Four separate analyses corresponding to a collision velocity of 10 knots were conducted, with the predicted shapes and stresses of the specimen deformation presented in Figure 12, Figure 13, and Figure 14. The top portions of each figure display the results using the LEPP model coupled with the CF and DP yield models. The bottom portions of each figure represent the results calculated with the LEP model of the ice material coupled with the CF and DP yield models, respectively. Interesting results were obtained when the CF yield model was applied to the LEPP model. In this case, the analysis results reveal the generation of large fragments due to spalling following the collision. On the other hand, the remaining three analysis models produce smaller fragments, which were more consistent with crushing phenomena. Significant differences in the observed loads during collision were noted depending on the characteristics of the analyzed fragments. Particularly, the model where large fragments were generated due to spalling exhibited significant deviations from the actual experimental results. These results underscore the importance of material model selection and its influence on collision analysis.

• Predicted strains on steel plate

Figure 15 and Figure 16 illustrate the strain rate history in the width direction $\left({\epsilon }_x\right)$ and height direction $\left({\epsilon }_y\right)$ of the plate specimen for various collision velocities. The maximum strain rates for ice properties based on collision velocity are summarized in Table 7, Table 8, and Figure 17. Both the LEPP and LEP models, in combination with the DP yield model, performed better than the combination with the CF model. Among these, the LEP model exhibited results closest to the experimental data, with error rates within 15%.

• Ice collision force prediction results

Figure 18 presents the collision force history applied to the plate specimen for various collision velocities. Figure 19 and Table 9 summarize the maximum collision force in conjunction with the strain rate for different material model combinations. Similar to the strain rate prediction results, the DP yield model outperformed the CF model in terms of accuracy. Particularly, the combination involving the LEP model exhibited a deviation of 22.1% from the experimental results. Among the material properties obtained from quasi-static compression experiments (W/O strain rate effect), this combination proved to be the closest to the experimental results. These results emphasize the significance of material model selection in collision-related analyses and highlight the importance of choosing a model that accurately predicts responses under various conditions.

4.3. Collision Analysis Reflecting Yield Ratio According to Strain Rate

This section presents the results of collision analysis, including the dynamic material properties of ice as a function of strain rate, as well as the models introduced in Section 4.1. Initially, the four material models introduced earlier were applied to calculate collision force and strain rates. The analysis was conducted for a collision velocity of 10 knots, and the calculated deformation shapes of the specimens are presented in Figure 21, Figure 21 and Figure 22. Figure 20 illustrates the failure pattern of ice predicted by combining the LEPP ice material model with the CF and DP yield models. Figure 21 illustrates the plastic strain distribution on the ice surface after the collision, whereas Figure 22 displays the Mises equivalent stress in the plate specimen. The analysis that combines LEP with the CF yield model does not consider the yield ratio in relation to strain rate. However, the analysis in Section 4.1 confirmed that the strain rate on the collision surface falls within a 20–30 mm/mm range. Based on the relationship between strain rate and yield strength ratio presented in Figure 5, the yield stress ratio could be estimated as 2.5. Additionally, the analysis that combined the LEPP model with the DP yield model illustrated the formation of spalling phenomena and the generation of large, fragmented particles.

• Strain of steel plate prediction results

The strain rate histories in the width $\left({\epsilon }_x\right)$ and height $\left({\epsilon }_y\right)$ directions of the plate specimen for various collision velocities are presented in Figure 23 and Figure 24, respectively. Figure 25 illustrates the maximum strain rates as a function of collision velocity for each material model. Detailed strain rate values are summarized in Table 10 and Table 11. Consistent with the results of the previous section, the DP yield model exhibited higher accuracy than the CF model. Among these models, the material model associated with the LPE model, which emphasizes plastic behavior, exhibited a lower deviation from the experimental results, with differences of 8.7% and 12.4%.

• Ice collision force prediction results

Figure 26 illustrates the history of collision force applied to the plate specimen as a function of collision velocity. The maximum strain rates for each material model are summarized in Figure 27 and Table 12. The DP yield model accurately predicted the ice loads applied to the plate specimen, outperforming the CF model. Furthermore, the material model associated with the LEP model, which accentuates plastic behavior, showed a deviation of 12.2% from the experimental results. This indicates a better fit compared to models that do not account for strain rate effects.

4.4. Evaluation of Suitability for each Combination of Ice Material Models

In this section, we compare the analysis results from Sections 3.3 and 3.4 with experimental data, summarizing the suitability of the analysis for each combination of material models. Below, the accuracies of the strain rate and collision force predictions for each ice material model are summarized.

• The strain rates for each material model were compared to identify models that closely match the experimental results. Figure 28 presents a comparison of the accuracy of the maximum strain rates in the width and height directions of the plate specimen as a function of collision velocity. Table 13 and Table 14 provide a comparison of the calculated and experimental values for the horizontal strain $\left({\epsilon }_x\right)$ and vertical strain $\left({\epsilon }_y\right)$ for each analysis combination. Here, we observed that dynamic material properties combined with the DP model (LEP-DP) accurately simulated ${\epsilon }_x$, with an 8.7% deviation, and ${\epsilon }_y$, with a 12.4% deviation. In summary, our findings demonstrated that considering dynamic effects yields better strain rate predictions, and when considering both width and height directions, the combination of the dynamic material properties of ice with the DP model is deemed the most appropriate model.

Table 13. Comparison between the experimental and FEA-derived maximum strain x (${\mathit{\epsilon }}_{\mathbf{x}}$).

Collision velocity	LEP-CF	LEP-DP	LEP-CF + Strain rate effect	LEP-DP + Strain rate effect
5 kts	6.90 × 10−4 (7.1%)	6.43 × 10−4 (0.2%)	5.59 × 10−4 (13.2%)	7.55 × 10−4 (−17.3%)
7 kts	8.06 × 10−4 (8.9%)	8.63 × 10−4 (2.3%)	7.22 × 10−4 (18.4%)	8.26 × 10−4 (6.7%)
10 kts	1.17 × 10−3 (8.3%)	1.11 × 10−3 (2.8%)	8.26 × 10−4 (6.7%)	1.06 × 10−3 (2.0%)
Avg. Difference	8.1%	1.7%	15.0%	8.7%

Table 13. Comparison between the experimental and FEA-derived maximum strain x (${\mathit{\epsilon }}_{\mathbf{x}}$).

Collision velocity	LEP-CF	LEP-DP	LEP-CF + Strain rate effect	LEP-DP + Strain rate effect
5 kts	6.90 × 10−4 (7.1%)	6.43 × 10−4 (0.2%)	5.59 × 10−4 (13.2%)	7.55 × 10−4 (−17.3%)
7 kts	8.06 × 10−4 (8.9%)	8.63 × 10−4 (2.3%)	7.22 × 10−4 (18.4%)	8.26 × 10−4 (6.7%)
10 kts	1.17 × 10−3 (8.3%)	1.11 × 10−3 (2.8%)	8.26 × 10−4 (6.7%)	1.06 × 10−3 (2.0%)
Avg. Difference	8.1%	1.7%	15.0%	8.7%

Table 14. Comparison of maximum strain y (${\mathit{\epsilon }}_{\mathbf{y}}$) from experiment and FEA.

Collision Velocity	LEP-CF	LEP-DP	LEP-CF + Strain rate effect	LEP-DP + Strain rate effect
5 kts	2.80 × 10−4 (31.7%)	2.40 × 10−4 (41.5%)	2.96 × 10−4 (27.9%)	3.32 × 10−4 (19.0%)
7 kts	3.12 × 10−4 (21.2%)	2.99 × 10−4 (24.5%)	2.90 × 10−4 (26.8%)	3.92 × 10−4 (0.90%)
10 kts	5.53 × 10−3 (47.9%)	3.76 × 10−3 (0.5%)	3.80 × 10−4 (−1.5%)	4.39 × 10−3 (−17.4%)
Avg. Difference	33.6%	22.1%	17.7%	12.4%

Table 14. Comparison of maximum strain y (${\mathit{\epsilon }}_{\mathbf{y}}$) from experiment and FEA.

Collision Velocity	LEP-CF	LEP-DP	LEP-CF + Strain rate effect	LEP-DP + Strain rate effect
5 kts	2.80 × 10−4 (31.7%)	2.40 × 10−4 (41.5%)	2.96 × 10−4 (27.9%)	3.32 × 10−4 (19.0%)
7 kts	3.12 × 10−4 (21.2%)	2.99 × 10−4 (24.5%)	2.90 × 10−4 (26.8%)	3.92 × 10−4 (0.90%)
10 kts	5.53 × 10−3 (47.9%)	3.76 × 10−3 (0.5%)	3.80 × 10−4 (−1.5%)	4.39 × 10−3 (−17.4%)
Avg. Difference	33.6%	22.1%	17.7%	12.4%

• Subsequently, we turned our attention to the maximum collision force exerted by the ice on the plate specimen. The maximum collision force applied by the ice to the plate specimen was compared for different ice material property models. Figure 29 compares the predicted collision force for each material model combination against the experimental data.
• Table 15 provides a comprehensive summary of the differences between experimental and analytical collision force for each material model combination. The inclusion of dynamic strain rate effects to the LEP material model with the DP yield condition (LEP-DP with strain rate effect) accurately predicts collision force with a 12.2% deviation, demonstrating high accuracy in both strain rates and collision force. This material model combination shows high accuracy in both strain rates and collision force.

Overall, upon comparing the different numerical analysis configurations analyzed in this study, our findings demonstrated that the LEP material model, combining the DP yield condition and strain rate effects, best matches the experimental results. Table 16 summarizes the comparison between both strain rates and collision force values obtained through experiments and numerical modeling. Notably, the strain rate effects in the analysis results were very similar to the experimental results, with deviations of 8.7% for ${\epsilon }_x$, 12.4% for ${\epsilon }_y$, and 12.2% for collision force. Therefore, our findings confirmed that the ice material properties formulated by combining the DP yield condition and strain rate-yield strength ratio within the LEP material model best approximate the results of the ice collision experiments.

Table 15. Comparison between the experimental and FEA-derived maximum force.

Collision Velocity	LEP-CF	LEP-DP	LEP-CF + Strain rate effect	LEP-DP + Strain rate effect
5 kts	31.5 (2.9%)	31.1 (4.2%)	27.5 (25.0%)	36.5 (-12.6%)
7 kts	44.7 (1.3%)	42.7 (5.7%)	36.9 (18.5%)	39.9 (11.8%)
10 kts	68.8 (3.7%)	63.9 (10.5%)	52.8 (25.0%)	61.8 (12.2%)
Avg. Difference	2.6%	6.8%	19.5%	12.2%

Table 15. Comparison between the experimental and FEA-derived maximum force.

Collision Velocity	LEP-CF	LEP-DP	LEP-CF + Strain rate effect	LEP-DP + Strain rate effect
5 kts	31.5 (2.9%)	31.1 (4.2%)	27.5 (25.0%)	36.5 (-12.6%)
7 kts	44.7 (1.3%)	42.7 (5.7%)	36.9 (18.5%)	39.9 (11.8%)
10 kts	68.8 (3.7%)	63.9 (10.5%)	52.8 (25.0%)	61.8 (12.2%)
Avg. Difference	2.6%	6.8%	19.5%	12.2%

Table 16. Comparative summary of the experimental and FEA-derived results.

FEA Result	Material model	5 kts	7 kts	10 kts	Mean Diff.
$\mathrm{Max}. \mathrm{strain} \mathrm{x} ({\mathsf{\epsilon }}_{\mathrm{x}}$) (Difference with experiment)	LEP -DP	$5.59{\times 10}^{-4}$ (−13.2%)	$7.22{\times 10}^{-4}$ (−18.4%)	$8.26{\times 10}^{-4}$ (−6.7%)	15.0%
	LEP -DP + Strain rate effect	$7.55{\times 10}^{-4}$ (−17.3%)	$8.26{\times 10}^{-4}$ (−6.7%)	$1.06{\times 10}^{-3}$ (−2.0%)	8.7%
$\mathrm{Max}. \mathrm{strain} \mathrm{y} ({\mathsf{\epsilon }}_{\mathrm{y}}$) (Difference with experiment)	LEP -DP	$2.96{\times 10}^{-4}$ (−27.9%)	$2.90{\times 10}^{-4}$ (−26.8%)	$3.80{\times 10}^{-4}$ (−1.5%)	17.7%
	LEP -DP + Strain rate effect	$3.32{\times 10}^{-4}$ (−19.0%)	$3.92{\times 10}^{-4}$ (−0.9%)	$4.39{\times 10}^{-4}$ (−17.4%)	12.4%
Max. force (Difference with experiment)	LEP -DP	27.5 (−25.0%)	36.9 (−18.5%)	52.8 (−25.0%)	19.5%
	LEP -DP + Strain rate effect	36.5 (−12.6%)	39.9 (−11.8%)	61.8 (−12.2%)	12.2%

Table 16. Comparative summary of the experimental and FEA-derived results.

FEA Result	Material model	5 kts	7 kts	10 kts	Mean Diff.
$\mathrm{Max}. \mathrm{strain} \mathrm{x} ({\mathsf{\epsilon }}_{\mathrm{x}}$) (Difference with experiment)	LEP -DP	$5.59{\times 10}^{-4}$ (−13.2%)	$7.22{\times 10}^{-4}$ (−18.4%)	$8.26{\times 10}^{-4}$ (−6.7%)	15.0%
	LEP -DP + Strain rate effect	$7.55{\times 10}^{-4}$ (−17.3%)	$8.26{\times 10}^{-4}$ (−6.7%)	$1.06{\times 10}^{-3}$ (−2.0%)	8.7%
$\mathrm{Max}. \mathrm{strain} \mathrm{y} ({\mathsf{\epsilon }}_{\mathrm{y}}$) (Difference with experiment)	LEP -DP	$2.96{\times 10}^{-4}$ (−27.9%)	$2.90{\times 10}^{-4}$ (−26.8%)	$3.80{\times 10}^{-4}$ (−1.5%)	17.7%
	LEP -DP + Strain rate effect	$3.32{\times 10}^{-4}$ (−19.0%)	$3.92{\times 10}^{-4}$ (−0.9%)	$4.39{\times 10}^{-4}$ (−17.4%)	12.4%
Max. force (Difference with experiment)	LEP -DP	27.5 (−25.0%)	36.9 (−18.5%)	52.8 (−25.0%)	19.5%
	LEP -DP + Strain rate effect	36.5 (−12.6%)	39.9 (−11.8%)	61.8 (−12.2%)	12.2%

5. Conclusion

This study reviews various material models and numerical simulation methodologies relevant to the impact of ice on structures, with a focus on developing a sophisticated numerical simulation model capable of accurately capturing the pressure characteristics induced by ice impacts. Using dynamic explicit finite element analysis, this study incorporated the elastic-plastic material properties of ice specimens, the effect of strain rate, yield condition, and the fracture behavior of materials. Initial analyses were conducted to validate necessary material constants for each ice material model, serving as a foundation for subsequent finite element simulations of ice impact scenarios. These simulations were iteratively refined to align closely with empirical impact load. The analysis process entailed exploring combinations of various material models. These included the linear elastic-perfectly plastic relationship, the linear elastic-plastic relationship, the crushable foam model, and the Drucker-Prager model, along with considering the strain rate dependency of ice. Collision simulations were instrumental in computing the strain and reaction forces of plate specimens, which were then compared with values measured in collision experiments. Our findings suggest that a material model that incorporates the Drucker-Prager yield criterion and a dynamic yield ratio within a Linear Elastic-Plastic framework most accurately reflects the experimental results. This study sought to validate the feasibility of simulating experiments with different material models of ice, which are essential for the analysis of ice-structure interactions. The validation of this numerical model enhances the predictive capabilities for offshore structural loads and allows for the derivation of more realistic ice material models through empirical comparisons.