The interaction between sea ice and structure is a complex non-linear dynamic process. According to the structural dynamics, the equation of motion of the sea ice collision problem can be expressed as:

where [M] is the mass matrix; [C] is the damping matrix; [K] is the stiffness matrix; _{x¨} is the acceleration vector; _{x˙} is the velocity vector; _{x} is the displacement vector and _{Fc} is the sea ice collision force vector. _Fc can be described as:

where _{Fs (h)} is the static ice force and _Fd is the dynamic ice force assumed to be a Gaussian process with a mean value of 0. The static ice forces can be described as:

where m is the collision surface shape factor (0.9 for circular cross section); σ_c is the compressive strength of sea ice; D is the structural outside diameter; I is the embedding coefficient and f_c is the contact coefficient. For cylindrical structures with a diameter of 2.5 to 10 m, the coefficient I⋅f_c is 0.4, according to the offshore standard IEC 61400-3 (International Electrotechnical Commission [IEC], 2009).

The Eq. 1 can be solved numerically to obtain the transient response of the sea ice collision problem. The explicit integration using the central differential method was used to calculate the accelerations, the collision forces and the structural internal forces.

Arbitrary Lagrange-Euler (ALE) and Euler algorithm are the commonly used methods in the simulation of FSI. According to the ALE algorithm, the meshed elements of the sea ice can be moved and co-existed with each other, rather than being fixed in one position. So, the consistency should be ensured by the continuities between the displacements and pressures on the contact surface of the fluid and ice, rather than the consistencies between the meshed fluid and solid elements. The velocity potential in the fluid domain can be introduced with Laplace’s equation in Eq. 4.

Further, the deformation continuum condition of the contact surface between fluid and sea ice can be expressed as:

In the ALE method, in addition to Lagrangian and Euler coordinates, another reference coordinate assumption needs to be defined. The derivative of any physical quantity f with respect to time t consists of two parts:

where _Xi and _xi are Lagrangian and Euler coordinates and w_i is the relative velocity. The governing equation of the fluid is derived based on the Navier–Stokes equation of incompressible fluid (Huang and Tuhkuri, 2019; Tang et al., 2020) and Eq. 6.

As a conservation of mass and momentum, it can be written as Eqs 7 and 8 using the ALE algorithm.

Thereafter, the conservation of energy can be expressed as:

Concurrently, for solid, the governing equation is required when moving solid elements are as follows:

where _ρ is the density; _v→, _a→, and _w→ are the velocity, acceleration and grid velocity of the material, respectively; σ¯ is the Cauchy stress; F_V is the volume force and e is the energy.

In LS-DYNA, the Euler description is used for the fluid solution and the Lagrangian description is used for the OWT structure and ice. The process of solving the FSI force is to solve the interaction between the both and the most common penalty method is also used for the solution method. Therefore, the coupling mechanism of fluids and solid can be achieved by solid-fluids penetration, in which a spring representing the pushing FSI force F is placed at a fluid node with the fluid penetrating depth _Δd into the solid. The FSI force at the node can be expressed as a product of stiffness and penetration depth:

where _pf is the coefficient of the penalty function; K is the volume modulus of the fluid element; V is the volume of the fluid element and A is the average area coupled with the fluid element.

The seawater is described by the Grüneisen state equation, which calculates the material compression force by defining the pressure volume. Based on the Grüneisen equation, the water pressure equation can be expressed as:

where _ρ0 and _ρ1 are the initial density and current density of the fluid, respectively; C is the velocity of sound shock wave; _γ0 is the Grüneisen constant; _α is the first-order correction factor of _γ0; E_C is the fluid-specific internal energy; _μ is the coefficient of fluid viscosity and S₁, S₂, and S₃ are the massless coefficients associated with the shock wave propagation.

The von Mises yield criterion is applied to determine the plastic phase of sea ice and the failure criterion of maximum failure strain and cut-off pressure are applied to determine the failure of sea ice. Meanwhile, the isotropic elastoplastic fracture characteristic is chosen as the constitutive model of ice. The sea ice presents an elastic behavior before reaching the yield state. After the first crack, the ice presents linear softening behavior with the developing of cracks. Finally, the ice exhibits a viscous fluid state when the crushing failure is completed (Derradji-Aouat, 2005). Based on the data from references (Song et al., 2019), the selected sea ice parameters are given in Table 1.

The NREL 5 MW monopile OWT is used in this study to investigate the interactions between the sea ice and vertical structures (Jonkman, 2005), as shown in Figure 1. The global coordinate system is located at the intersection of the wind turbine center and the mean sea level (MSL) with the x and y directions corresponding to the downwind and horizontal directions, respectively.

The rotor nacelle assembly (RNA) is supported by an extended tower with a height of 77.5 m placed on top of the monopile supporting structure. The water depth is selected as 20 m and the fluid parameters used in Eq. 12 are shown in Table 2. The material density for the tower is set to be 8,500 kg/m³ instead of the default value of 7,850 kg/m³ to account for the additional mass of flanges, bolts, internal equipment, paint, mass of secondary structures, etc.

The fully ice-OWT interaction analysis process is given in Figure 1. The IEC Kaimal turbulence model and the exponential coherence model are used to generate turbulent wind field inputs with a turbulence intensity in category A by TurbSim (Kaimal et al., 1972) and wind turbine loads can be calculated through FAST. Furthermore, the apparent fixity length (AFL) method is used herein to consider the soil-structure interaction (SSI) effect. In the AFL method, the properties of the fictive beam are tuned such that the mudline displacement and rotation for both the models would be the same, when loaded by a mudline shear force and bending moment that are representative. The fictitious length in the AFL method is 17.5 m, which is 2.9 times of the diameter of the monopile foundation.

In the analysis of ice-OWT interaction, the constitutive model of ice is the key for the simulation. When sea ice collides with an OWT, the interaction area of the ice sheet will be locally crushed and broken after cracks in the ice sheet. In order to verify the correctness of the sea ice model used in this study, a simple rigid column and sea ice collision numerical model are established and a theoretical pressure-area relationship proposed by Masterson and Frederking (1993) in international organization for standardization (ISO)-ALIE (International Organization for Standardization [ISO], 2010) is adopted to compare calculation results using such ice constitutive model. The calculated contact pressure is verified against ISO-ALIE (Figure 2). It is shown that the collision pressure can be expressed as an inversely proportional function of the collision area under the static ice forces and the comparison results are more consistent indicating the ice numerical model could be further used for FSI analysis.

In order to verify the influence of the main constitutive and failure parameters of sea ice on the ice-OWT interaction, numerical simulations of the interaction between sea ice and NREL 5 MW monopile OWT are carried out. In the simulations, different failure strains, failure pressure and plastic hardening modulus are considered. The detailed parameters are shown in Table 1. The plastic hardening modulus E_P is related to the elastic modulus E and the tangential modulus E_tan:

Figures 3–5 show the comparisons of the ice forces with the different ice parameters in Table 3. The extreme, mean and SD ice force all increase linearly with the increasing ice parameters. The fitted linear relations in the graphs show that the maximum load can generate the largest slopes and the fastest growth rates. The slopes of the maximum dynamic ice forces under the influence of failure strain and the failure pressure are 3.5356 and 3.4955, which are significantly greater than the influence of plastic hardening modulus, while the SD curves are not extremely different. The detailed comparisons of the mean, SD and extreme ice load values are shown in Tables 4–6.

In summary, the effect of plastic hardening modulus, failure pressure and failure strain of ice are significant in influencing ice force; therefore, these variations can play an important role in numerical simulations. Based on the data from references and the discussion of ice parameters above, the sea ice parameters selected are shown in Table 1. The failure strain value of 0.3 is used to simulate the crushing failure of sea ice.

The interaction numerical model with FSI consists of three main parts: the sea ice, NERL 5 MW monopile OWT and seawater. The fluid model is 80 m × 72 m (length × width) with 20 m depth of water, which is naturally modeled by null materials with non-reflecting boundary (infinite boundary condition), as shown in Figure 6A. Considering the large amount of computational effort FSI requires, the fluid-solid element mesh is no longer refined. Instead, a separate overlapping modeling approach is applied to divide the sea ice and the fluid model into Lagrangian meshes and Euler meshes, respectively, so that the two meshes can overlap and move spatial freely.

The three-dimensional effect and meshing of the interaction between sea ice and the OWT with FSI coupling are given in Figure 6. It is worth noting that F-A and S-S directions of OWTs are defined under the coordinate system (X, Y, and Z), as shown in Figure 6. The F-A direction is along the X direction, i.e., the sea ice drifting direction and the S-S direction is along the Y direction. As indicated in the figure, the ice model is immersed in a non-viscous, non-rotating and incompressible fluid and the sea ice is discretized using fully integrated solid elements. In order to prevent the ice boundary condition affecting the stress wave reflection, the outer boundary of the ice solid element adopts the non-reflecting boundary condition to ignore the ice size and shape effect, as shown in Figure 6B. Furthermore, to avoid an hourglass phenomenon caused by explicit integration calculations in ANSYS/LS-DYNA, the Flanagan–Belytschko stiffness method is applied to evaluate whether the calculation process produces an hourglass. The contact between the sea ice and the structure is set as the eroding surface-to-surface contact. This contact condition ensures that the remaining elements can still be considered after crushing and removal of the crushed elements.

Concurrently, the SHELL163 element is used to establish the plastic kinematic model of the steel and the RNA of the OWT is simplified as a concentrated mass, which is applied to the wind turbine loads calculated by FAST and wind turbine loads include shear force in F-A direction and bending moment around S-S direction. In order to ensure the accuracy of numerical simulation, while improving computational efficiency, the mesh elements are only encrypted at the location of the ice crushing and structure collision, which are 0.15 and 0.2 m for the ice and OWT, respectively.

During the fully ice-OWT with FSI interaction simulation, the OWT structure is defined as Lagrangian element, and the seawater are ALE elements. The ALE algorithm and contact collision algorithm in LS-DYNA are used to couple ice-seawater, OWT-seawater and ice-OWT collision settings. In addition, if the structural grid of ice overlaps with that of seawater, in order to avoid the hydrodynamic solution error caused by the overlapping of different material grids, the seawater material occupying the ice area should be removed.

The dynamic ice forces simulated from the numerical model with FSI in this study are compared with the dynamic ice forces obtained by the DEM (Ji et al., 2013); dynamic ice forces calculated by the semi-empirical formulation using ice force spectrum were proposed by Kärnä and Qu (2004) and the static ice force formulas in the IEC standard.

The simulated ice force from the proposed FEM, together with ice force spectrum fitting and the DEM are compared against the design values, was calculated by the Korzhavin formula, IEC standard, as given in Figure 7. For the static ice force, the ice force from present method provides better approximations with results closest to the IEC formula, the difference of maximum ice force is only 391 kN. Furthermore, the American petroleum institute (API) standard gives the more conservative upper bound of the maximum ice force design value, which can attain 3,600 kN. On the other hand, the mean ice force of three dynamic methods, namely, present FEM, DEM and ice force spectrum fitting method are 727.93, 374.29, and 833.57 kN, respectively, and the ice force amplitude obtained by DEM is the smallest related to other two methods. It can be seen that proposed FEM in this study basically consistent with the simulation result of ice force spectrum fitting method. On the basis of the above discussion, the interaction simulation of ice and OWT through proposed method in this article can be well verified.

The structural responses are in the form of forced vibrations. The displacement at the collision area of the OWT supporting structure under LC 1.1 is shown in Figure 8. In the F-A direction of the OWT, which is the primarily loaded direction, the displacement at the collision area on the supporting structure can reach the maximum of 0.086 m. In addition, the structural responses in the y direction are larger compared with z direction, although both of which are at a low level of structural responses.

The effective stress damage distribution of the structure at different moments is shown in Figure 9. As indicated in Figure 9A, the main damage to the OWT occurs only within the collision zone of contact, with little or no damage occurring at other locations. As the collision progresses in Figure 9B, the damage surface of the structure gradually spreads to the surrounding area. At 19.8 s, the effective stress continues to increase slowly and the damaged surface area gradually stabilized in Figure 9C. The collision mainly occurs in the elastic phase within the material failure strain value of 0.348; therefore, no failure occurs at the collision area.

On the other hand, the breaking of sea ice is also important in OWT structural design. It may cause ice-induced resonance generating a large ice force resulting in a severe vibration response of the structure. The main frequency of sea ice loads is between 0.2 and 0.84 Hz (Liu and Ji, 2018), which covers the fundamental frequency of the monopile OWT, therefore ice-induced resonance analysis, i.e., structure steady-state vibration is required. Numerical simulations with FSI are performed to provide the velocity range of the possible ice-induced resonance.

Zhang (2002) pointed out that the self-excited vibration of sea ice at ice velocities of 1–2 cm/s can lead to steady-state vibration of marine structures. The results shown in Figure 10 indicate that steady-state vibration of the structure is of high probability to occur when the sea ice velocity is in the range of 2.5–3.5 cm/s. The simulation results are consistent with the measured data in Zhang (2002) and Zhang et al. (2018).

The FSI effect includes the influence of seawater on OWTs and ice. To systematically explore the effect of fluid in full collision process, ice-OWT interaction numerical simulations are studied with the initial velocity of ice, as given in LC 1.2 of Table 7. Comparable cases of coupling and without coupling FSI are considered. Further, the collision force and energy variation are compared in terms of effect of FSI.

The ice-OWT interaction is performed at an initial velocity of 1.2 m/s in LC 1.3. The time history of the collision force with and without FSI is shown in Figure 11 where the general trend of the two cases mostly follow the same route, while the case with FSI shows an obvious hysteresis. The ice force with FSI fluctuates more frequently than that without FSI coupling because of the large number of element failures in the sea ice caused by FSI. The maximum ice force with and without coupling is 0.458 and 0.569 MN, respectively, indicating that the fluid can impair the load within the allowable error range. The ice force reaches its maximum between 0.8 and 1.2 s, which can be corresponded to the larger growth and the fluctuations in sea ice erosion energy. Thereafter, as the structural deformation and the collision force continue to increase, the sea ice slows down and finally bounces back with ice force gradually decreases until unloading completed. Therefore, the seawater can affect not only the ice force amplitude, but also the interaction duration, i.e., lowering the maximum ice force and lengthening the duration of the interaction.

The energy variation of the collision system is analyzed to understand the fully coupled collision process. The system energy shown in Figure 12A satisfies the law of energy conservation. The energy exchange process goes from a transient state to a steady state during the extrusion of the sea ice and the OWT structure. The system energy mainly consists of kinetic energy and internal energy. The kinetic energy drops rapidly to a minimum inflection point of 0.027 MJ at approximately 2 s after the ice floe coming into contact with the structure. It means that the ice velocity gradually decreases to zero at this moment. After the sea ice bounces back, the kinetic energy of the collision system becomes stable and gradually decreases to the minimum value. Additionally, the internal energy of the system increases slowly and gradually stabilized.

Except for the system energy of kinetic energy and internal energy, the interface slip energy and the hourglass energy generated by the explicit integral calculation are dissipated during the energy exchange process, as shown in Figure 12A. Interfacial slip energy should be kept at a small value to reduce the energy loss. Furthermore, the hourglass energy ratio derived in the numerical calculation is 0.037% (less than the requirement of 5%), indicating that the model, mesh density and calculation results are satisfactory.

The erosion energy of sea ice is shown in Figure 12B. With consideration of FSI, the erosion energy and the erosion failure effect of the sea ice are seen greater by comparing with the case without FSI. This indicates that fluids accelerate sea ice breaking and failure. Moreover, the ratio of the kinetic energy to the internal energy of the OWT supporting structure is shown in Figure 12C, in which smoother fluctuations and smaller energy are observed in the case with FSI coupling.