Front. Mar. Sci. Frontiers in Marine Science Front. Mar. Sci. 2296-7745 Frontiers Media S.A. 10.3389/fmars.2022.956032 Marine Science Original Research Dynamic analysis of offshore wind turbines subjected to the combined wind and ice loads based on the cohesive element method Wang Bin 1 2 † Liu Yingzhou 3 4 † Zhang Jianhua 5 Shi Wei 3 6 * Li Xin 3 4 Li Ying 7 1 Key Laboratory of Far-shore Wind Power Technology of Zhejiang Province, Hangzhou, China 2 Powerchina Huadong Engineering Corporation Limited, Hangzhou, China 3 State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian, China 4 Institute of Earthquake Engineering, Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian, China 5 College of Aerospace and Civil Engineering, Harbin Engineering University, Harbin, China 6 Deepwater Engineering Research Center, Dalian University of Technology, Dalian, China 7 Chinese-German Institute of Engineering, Zhejiang University of Science and Technology, Hangzhou, China

Edited by: Yong Cheng, Jiangsu University of Science and Technology, China

Reviewed by: Zhenju Chuang, Dalian Maritime University, China; Bang-Fuh Chen, National Sun Yat-sen University, Taiwan

*Correspondence: Wei Shi, weishi@dlut.edu.cn

†These authors have contributed equally to this work

This article was submitted to Ocean Solutions, a section of the journal Frontiers in Marine Science

 06 09 2022 2022 9 956032 29 05 2022 15 08 2022 Copyright © 2022 Wang, Liu, Zhang, Shi, Li and Li 2022 Wang, Liu, Zhang, Shi, Li and Li

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Ice loads are an important and decisive factor for the safe operation of offshore wind turbines (OWTs). In severe environment load cases, it shall lead to prominent ice-induced vibration and ice-induced fatigue failure of OWT structures. Based on the cohesive element method (CEM) and considering the pile–soil interaction used by nonlinear distributed springs, the full interaction model of the ice and monopile OWT structure with an ice-breaking cone in a cold sea region is established in this study. Furthermore, the Tsai-Wu failure criterion and the empirical failure formula of maximum plastic failure strain are used to describe the mechanical behavior of ice bending failure in the collision simulation tool LS-DYNA, and the dynamic ice loads under different ice velocities and cone angles are statistically analyzed. Finally, according to the interaction process between sea ice and OWT containing the ice-breaking cone, the dynamic response of OWT under the combined wind and ice loads is studied, and the most reasonable ice-breaking cone angle is determined. The results show that the method adopted in this paper can well simulate the bending failure process of sea ice. Concurrently, the cone angle has a significant impact on the dynamic response and damage of the OWT, and the recommended optimal cone angle is 60.

 ice-structure interaction model offshore wind turbine cohesive element method bending failure dynamic analysis 1 Introduction

As a renewable and clean energy, wind energy plays an important role in the field of power generation, which has led to the rapid development of offshore wind turbines (OWTs) (Zhang et al., 2020; Zhao et al., 2021). Since the Europe 2020 agreement set the target of at least 20% of the electricity production coming from renewable energy sources by 2020, the cumulative installed capacity of OWTs in Europe has more than quadruplicated. In the OWT research, the dynamic analysis of OWTs (Wang et al., 2022; Ren et al., 2022) was investigated under combined environment load conditions recommended by design standards of IEC and DNV (IEC, 2012; DNV GL, 2016). However, as the construction of offshore wind farms gradually advances to the cold sea, the OWT shall be continuously affected by large-scale floating ice in the cold sea area. Therefore, ice-induced structural vibration has become a major factor affecting the safe operation of the OWT structure. In order to avoid ice-induced structural vibrations resulting in excessive motion response of OWTs, positive and inverted ice-breaking cones are extensively installed in the collision position of OWT substructure subjected to ice floe. It can reduce the ice loads by transforming the failure mode of sea ice from crushing failure mode to bending failure mode.

Currently, most studies on the dynamic response of OWTs under sea ice were carried out through ice force time histories or simplified ice force models. Based on the bending failure model of ice, Shi et al. (2016) developed a calculation program for the interaction between ice and OWTs in the coupled numerical analysis software, HAWC 2. The coupled dynamic response characteristics of the OWT were analyzed under different ice thicknesses and speeds. A forced ice-induced vibration model was used to simulate the motion response of OWTs under the action of turbulent wind and ice by (Ye et al., 2019). Based on the FAST numerical calculation software, a simplified ice force model was adopted to compare the dynamic response of monopile OWTs under the floating ice by Heinonen and Rissanen (2017). The Karna ice force spectrum model was implemented in analyzing the influence of ice loads and wind loads on the OWT fatigue by Wang (2015), and found that the influence of ice loads on the fatigue damage of OWT substructure was greater than that caused by wind loads.

Involving the influence of factors such as the physical and mechanical properties of sea ice, the sea ice dynamics, and the structural dynamic characteristics of OWTs, the coupling effect between sea ice and OWT structures is extremely complex. Furthermore, due to the complexity of the interaction coupling effect, OWTs are generally simplified as rigid structures in the ice force studies. On this basis, in order to better simulate the dynamic ice force, by establishing a three-dimensional sea ice-marine structure model, the nonlinear finite element method, discrete element method, and cohesive element method (CEM) are used to investigate the ice-structure interaction research. The finite element method was applied to simulate the collision process between floating ice and a simple rigid OWT substructure by (Song et al., 2019), and ensured that when a mesh size of the ice model is 1.5 times the ice thickness, the simulations could present more accurate estimations in terms of maximum ice loads for all static ice tests. Meanwhile, (Liu et al., 2022) considered the fluid–structure interaction and simulated the ice–OWT interaction process using the nonlinear finite element method, and discussed the crushing failure mode of sea ice. The discrete element model with a bonding-fragmentation effect was established to calculate the sea ice action process of the offshore platform structure with an ice-breaking cone in Bohai Sea, focusing on the ice loads and stress distribution of the cone by (Di et al., 2017). (Jang and Kim, 2021) estimated the dynamic ice loads on a conical marine structure on the basis of the discrete element method. The findings revealed that the comparisons between the simulated ice force and the experimental results were relatively consistent. For the CEM, the dynamic ice force of vertical structures based on the CEM was carried out by (Statoilhydro et al., 2009) and compared them with the field test results. The results showed that the proposed numerical method captured a large number of qualitative observations and quantitatively validates the ice force time history. (Zhan et al., 2021) used the CEM to establish a numerical model of layer ice, and performed a numerical simulation of the interaction between layer ice and rigid cone structure to study the effects of ice resistance cone angle, ice thickness, and other parameters on the ice loads.

From the above research results, it can be seen that the crushing failure and bending failure caused by the interaction between sea ice and the OWT structure, as well as the cracking and accumulation of crushed ice, have a significant impact on the ice loads and the dynamic response of the OWT structure. However, most of the studies about ice loads used the time-history fitting input method of ice force, ignoring the coupling effects of ice loads and structures. In addition, in the simulation of ice loads, in order to simplify the coupled calculation, the OWT structure is usually regarded as a rigid body, and the structural motion response and pile–soil interaction are ignored. In this paper, based on the CEM and considering the pile–soil interaction, a coupled model of ice–OWT with cone interaction is established. Furthermore, based on the established interaction coupled model, the dynamic response and cone damage analysis of the OWT under operating conditions are studied. The effects of ice speeds and cone angles on the dynamic ice loads and dynamic response of OWTs are discussed, and the optimal cone angle is determined.

 2 Theoretical model 2.1 Mathematical model of collision analysis

According to the finite element method of the dynamic nonlinear collision problem, the motion governing equation of the ice–OWT structure interaction subjected to wind and ice loads can be expressed as:

[ M ] { x ¨ ( t ) } + [ C ] { x ˙ ( t ) } + [ K ] { x ( t ) } = { f w i n d ( t ) } + { f i c e ( t ) } + { H ( t ) }

where [M], [C], and [K] are the mass, damping, and stiffness matrices of the OWT structure, respectively; { x ¨ ( t ) } , { x ˙ ( t ) } , and {x(t)} are the nodal acceleration, velocity, and displacement vectors, respectively; {f wind (t)} and {f ice (t)} are the external wind and ice load vectors, respectively. When considering the hourglass energy effect, it is necessary to add the hourglass resistance vector {H(t)}.

In LS-DYNA, the central difference method of explicit integration is mainly used to first obtain the time history of node displacement, and then calculate the time history of acceleration, collision force (rcforce), and internal force through the time history of displacement (Hallquist, 2006). The relationship between acceleration and force under ice is calculated by:

{ x ¨ ice } = [ M ] − 1 ( { f i c e ( t ) } − [ ∑ ( ∫ Ω B T σ n d Ω + F h g ) + F c o n t ] )

where x ¨ ice is the node acceleration under ice, F hg is the hourglass resistance, and F cont is the constant force.

For the discretized dynamic balance equation above, LS-DYNA calculates the structural motion response based on the explicit central difference method (Hallquist, 2006). The specific solution process is as follows:

{ x ¨ ( t n ) = M − 1 [ f wind ( t n ) + f ice ( t n ) − f int ( t n ) + H ( t n ) − C x ˙ ( ( t n + t n − 1 ) / 2 ) ] x ˙ ( ( t n + 1 + t n ) / 2 ) = x ˙ ( ( t n + t n − 1 ) / 2 ) + x ¨ ( t n ) ( Δ t n − 1 + Δ t n ) / 2 x ( t n + 1 ) = x ( t n ) + x ˙ ( ( t n + 1 + t n ) / 2 ) Δ t n

where {f int (t n )} is the internal force vector for any time tn.

In Eq. (3), the solution by the explicit central difference method does not need to calculate the overall matrix, nor does it need to perform complex equilibrium iteration, which effectively avoids the problem of nonlinear convergence. However, the solution stability is conditional; in order to ensure the stability of the numerical calculation, LS-DYNA adopts a variable integration step size smaller than the critical value. Critical values applicable to various element types can be expressed as:

Δ t e = α ( l e / c )

where Δt e represents the critical time step size of the element e; α represents the time step scaling factor, and the default value is 0.9; l e represents the feature size of the element; and c represents propagation speed of longitudinal waves.

 2.2 Soil–structure interaction

In this paper, nonlinear distributed springs are adopted in simulating the pile–soil interaction of OWT foundations under the sea ice and wind. Based on the p–y curve method, the stiffness and damping of the spring element are determined, among which the spring damping cs adopts frequency-independent radiation damping per unit length (Wang et al., 1980):

c s = 4 ρ v s D

where ρ is the soil density, vs is the shear wave velocity of soil, and D is the diameter of the monopile foundation.

According to the measured geological and environmental conditions of the offshore wind farm located in the southeastern offshore regions in China, the empirical formula of soil shear wave velocity is proposed in Eq. (6) (Lan et al., 2012).

v s = 150 H 0.243

where H is the buried depth.

Moreover, considering the dynamic effect of the nonlinear spring elements, it is assumed that the distributed spring dynamic force can be represented by the amplification factor and the static force (Gladman, 2013):

F d = ( 1 + k d v v ice ) · F s

where Fd and Fs are distributed spring dynamic force and static force, respectively; kd is the dynamic amplification factor; v is the absolute value of the relative velocity at both ends of the distribution spring; and v ice is the sea ice drift speed.

 2.3 Material model 2.3.1 Cohesive element model

The cohesive element model is a complex numerical algorithm extended under the finite element framework of the cohesion model (Hillervborg et al., 1976). This method divides the ice mesh into traditional solid elements and cohesive elements. Firstly, the solid ice elements are discretized, and then cohesive elements are inserted between the solid ice elements in the horizontal and vertical directions, as shown in Figure 1 . Sequentially, solid ice elements share nodes with cohesive elements to transfer deformations and stresses. Through the action of external loads, the solid ice element is deformed, which, in turn, leads to the corresponding displacement of the cohesive element nodes. During the interaction simulation, once the cohesive element reaches the maximum separation displacement, it is considered to be completely destroyed and removed, resulting in the formation of obvious cracks on the ice surface and inside. Meanwhile, some fractured solid ice elements shall adhere to the surface of the ice sheet to form the accumulation of ice fragments. It is worth noting that the failure of cohesive elements is theoretically based on fracture mechanics; however, the actual CEM is based on damage mechanics.

The cohesive element method (Hillervborg et al., 1976).

The cohesive element adopts the cohesive force constitutive model, and the failure of the element follows the traction force–separation displacement failure criterion. Three commonly used traction force–separation displacement curves, namely, bilinear softening, exponential softening, and trapezoidal softening, are shown in Figure 2 . The traction force–separation displacement curve is mainly determined by fracture energy, fracture strength, and curve forms above. The area formed by the curves and the coordinate axis is the maximum energy required for the failure of the cohesive element, that is, the fracture energy.

Cohesive element failure criterion models.

The relationship between the traction force and the separation displacement of the linear softening model can be expressed as:

T = { T 0 δ n 1 δ n , 0 ≤ δ n δ n 1 T 0 ( δ n f − δ n 1 ) ( δ n f − δ n ) , δ n 1 ≤ δ n δ n f 0 , δ n ≥ δ n f S = { S 0 δ n 1 δ τ , 0 ≤ δ τ δ τ 1 S 0 ( δ τ f − δ τ 1 ) ( δ τ f − δ τ ) , δ τ 1 ≤ δ τ δ τ f 0 , δ τ ≥ δ τ f

where T and S are tensile and shear stresses, respectively; T 0 and S 0 are the tensile and shear critical stresses, respectively; δ n and δ τ are separation displacements for tensile and shear fracture, respectively; δ n 1 and δ τ 1 are the corresponding tensile and shear separation displacements at this time step, respectively; and δ n f and δ τ f are the maximum separation displacements for tensile and shear fracture, respectively.

 2.3.2 Sea ice constitutive model

For the solid ice elements, in order to reflect the generation and expansion of ice micro-cracks in the interaction between ice and OWTs, a homogenized elastic–plastic linear softening constitutive model is adopted in the presented study. The constitutive relationship of ice is shown in Figure 3 . The material behaves elastically until the sea ice reaches yield strength. After reaching the yield point σ c , the ice material begins to enter the linear softening stage, which means that with the increase of plastic strain, the effective stress is greatly reduced, and the sea ice is more easily crushed. Furthermore, when the failure strain is reached, ϵ f , the sea ice elements are completely broken. The failure strain is used to control the failure deletion of solid ice elements in LS-DYNA and can be expressed as:

ϵ f = ϵ 0 + ( P 10 8 − 0.6 ) 2

where ϵ 0 is the initial strain, which can be adjusted according to the ice test data. Based on the research results of the Norwegian University of Science and Technology, such parameter is selected as 0.01; P is the hydrostatic pressure.

Ice solid element constitutive model.

Concurrently, considering the temperature effect and the sensitivity of ice strength to hydrostatic pressure, the Tsai-Wu failure criterion (Derradji-Aouat, 2003) is an effective criterion for describing the bending failure of sea ice, which can reflect the tension, compression, and shear of the ice element in the three-dimensional stress space. Meanwhile, the Tsai-Wu failure criterion is adopted as the criterion for ice crack initiation and propagation, which can be re-developed through the LS-DYNA user-defined subroutine. The Tsai-Wu failure criterion can be defined as:

f = J 2 − ( b 0 + b 1 P + b 2 P 2 ) = 0

where, b 0, b 1, and b 2 are the constant coefficients obtained for fitting data from a triaxial experiment; and J 2 is the deviatoric stress, which is written as:

J 2 = 1 6 [ ( σ x − σ y ) 2 + ( σ x − σ z ) 2 + ( σ y − σ z ) 2 ]

where σ x , σ y , and σ z are the first, second, and third principal stresses, respectively.

 3 Numerical model of the interaction between ice and OWTs 3.1 Structural parameters of OWTs

Based on the DTU 10 MW baseline OWT (Velarde and Bachynski, 2017) proposed by the Technical University of Denmark (DTU) in cooperation with Vestas Company, an ice–OWT interaction model is established to conduct the dynamic response of OWT under wind and ice. As shown in Figure 4 , the DTU 10 MW OWT is mainly composed of a rotor nacelle assembly (RNA), a supporting structure, and a pile foundation. The height from the hub to the water surface is 119 m, and the design water depth is selected as 20 m. The main parameters are shown in Table 1 .

Schematic diagrams of the DTU 10 MW OWT (Velarde and Bachynski, 2017).

Main specifications of the DTU 10 MW OWT (Velarde and Bachynski, 2017).

OWT Parameters Values
Rated power 10 MW
Blade numbers 3
Hub height 119 m
Rotor diameter 178.3 m
Cut-in/rated/cut-out wind speeds 4, 11.4, 25 m/s
Cut-in/rated rotor speeds 6.9, 12.1 rpm
Rotor mass 227,962 kg
Nacelle mass 446,036 kg
Design depth 20 m
3.2 Sea ice model parameters

Traditional ice solid elements are modeled with eight-node hexahedral elements. However, for the bending failure of ice, the hexahedral element mesh shall result in breaking paths larger than the actual ice crack length, which can lead to additional dissipation of fracture energy with a non-negligible effect on the simulation results. Therefore, in this study, the six-node regular triangular prism meshing method is used to simulate the ice fracture and accumulation phenomena, as shown in Figure 5 .

Sea ice finite element-cohesive element model.

The constitutive parameters of sea ice materials are affected by many factors. In this paper, the parameters of the elasto-plastic constitutive model of sea ice are determined based on the relevant research on the physical and mechanical properties of ice materials at home and abroad, and the effects of temperature, salinity, and porosity of sea ice on the properties of sea ice are considered, as shown in Table 2 . Furthermore, Timco and Weeks (2010) obtained the tensile and shear rupture strengths of cohesive elements based on experimental data. In addition, Dempsey et al. (1999) carried out experimental studies on thick ice layers, and the fracture toughness values obtained showed that the fracture energy of ice ranges from 1 to 40 N/m. Based on the above research results, the parameters of the ice cohesive element are selected, as tabulated in Table 2 .

Main material parameters of ice (Timco and Weeks, 2010).

Parameters Values Parameters Values
Elastic modulus E Density ρ Poisson’s ratio ν Yield stress σ y Flexural strength σ f Structure-ice coefficient of friction 5 GPa910 kg/m3 0.332.0 MPa1.0 MPa0.2 Transverse tensile strength T x Longitudinal tensile strength T z Transverse shear strength S x Longitudinal shear strength S z Transverse fracture energy G І Longitudinal fracture energy G II 0.6 MPa0.5 MPa0.8 MPa0.7 MPa30 N/m30 N/m

In order to avoid the influence of the size of the ice sheet on the bending failure of the ice element in the contact area, and considering the calculation efficiency, the size (length × width: 36 m × 20 m) of the ice sheet is determined based on the elastic foundation beam theory. The side that collided with the ice sheet is set as a free boundary condition, and the other three sides only constrain the y-direction degree of freedom and set a non-reflection boundary condition to exclude the influence of the reflection of stress waves from the boundary on the numerical simulation results. Meanwhile, an elastic base plate is established to simulate seawater buoyancy, as shown in Figure 6 .

Ice–OWT interaction coupling model.

 3.3 Pile–soil interaction parameters

According to the soil parameters in Ref (Yang et al., 2020)., the site soil conditions shown in Figure 7 are adopted in the monopile OWT. Soil conditions are designed as sandy soil with a depth of 45 m. The internal friction angle is 36°, and the saturation weight is selected as 20 kN/m3.

Pile–soil interaction model based on nonlinear distributed springs.

During the interaction between sea ice and OWTs, the pile foundation of OWT is mainly affected by the horizontal direction, and the displacement variations caused by the vertical pile side friction can be ignored. Sequentially, based on the p–y curve method, a nonlinear distributed spring in the horizontal direction is used to simulate the pile–soil interaction under the combined action of wind and ice, as shown in Figure 7 . At different depth positions in the horizontal x, y direction, one node of the distribution spring is connected to the pile foundation, and the other node is set as the fixed end.

 3.4 Contact settings and hourglass control

In order to ensure that the initial contact is eliminated between the OWT substructure and ice, an initial collision distance of 0.1 m is set in the fully coupling interaction model. The contact for eroding surface to surface is applied to the interaction of ice and OWT. However, in order to reduce the computation time, the contact between ice fragments adopts eroding single surface contact. Simultaneously, the contact stiffness is increased by four to five times to reduce unnecessary penetration between ice and the OWT substructure. It is found that such contact stiffness can obtain better numerical simulation results through tests (Liu et al., 2022).

Thereafter, the hourglass mode is a non-physical zero-energy deformation mode that results in zero strain and stress, which can cause errors in the calculation results in LS-DYNA. When the hourglass energy exceeds 10% of the total internal energy of the structure, the analysis results can be considered invalid. Therefore, the hourglass phenomenon needs to be effectively controlled in the ice–OWT interaction. This paper found that the rigid hourglass control can effectively suppress the hourglass phenomenon compared with the viscous hourglass control by conducting a comparative study on the hourglass control.

 4 Calculation results and discussion 4.1 Design of load cases

Considering the pile–soil interaction, this paper focuses on the effects of different ice speeds and cone angles on the motion response of the DTU 10 MW OWT under the wind and ice loads. Based on the relevant sea conditions in the Bohai Sea and the requirements of the IEC-61400-1 specification (IEC, 2012), and combined with the cut-in and cut-out wind speeds of the DTU 10 MW wind turbine, the typical combined wind and ice conditions are selected as shown in Table 3 .

Designed combined load conditions.

Load conditions Wind speeds Turbulence intensity Ice thicknesses Ice speeds Cone angles
Ice speed conditions LC 1.1 11.4 m/s 0.146 0.4 m 0.2 m/s 60°
LC 1.2 11.4 m/s 0.146 0.4 m 0.4 m/s 60°
LC 1.3 11.4 m/s 0.146 0.4 m 0.6 m/s 60°
Cone angle conditions LC 2.1 11.4 m/s 0.146 0.4 m 0.4 m/s 45°
LC 2.2 11.4 m/s 0.146 0.4 m 0.4 m/s 60°
LC 2.3 11.4 m/s 0.146 0.4 m 0.4 m/s 75°

The average reference wind speed at the hub height of the OWT is a rated wind speed of 11.4 m/s and a specified turbulence intensity of 0.146 is selected in the presented study. Using the IEC Kaimal turbulence model and exponential coherence model, a turbulent wind field with a time step of 0.02 s is generated by TurbSim (Kilcher, 2012). The time-domain variation curve of wind speed at the hub height is shown in Figure 8 . In addition, the cone angle conditions (LC 2.1–LC 2.3) only changed the cone angle by keeping the identical diameter of the water surface of the cone, that is, ensuring a uniform collision position of sea ice, as shown in Figure 9 .

Synthesized stochastic wind speeds based on IEC Kaimal turbulence model.

Schematic diagram of ice-breaking cone angle based on same collision location.

 4.2 Sea ice bending failure and ice force verification 4.2.1 Bending failure mode of ice

Figure 10 shows the simulation results of ice bending failure and ice fragment accumulation effect under colliding with a constant velocity of 0.4 m/s. When the sea ice makes initial contact with the cone, the front of sea ice crawls a short distance along the cone. As the bending stress of the sea ice reaches its bending strength, tiny radial cracks first appear on the ice surface. This is due to the deletion of cohesive elements at the front of ice, resulting in radial failure exhibiting a non-simultaneous failure mode subjected to non-uniform forces. Furthermore, the crushed ice broken from the ice floe shall climb up and down slightly on the cone surface and remain at the height of the water level under the combinations of gravity and the buoyancy provided by the elastic base in Figure 6 . Consequently, the accumulation phenomenon of ice fragments can be observed. Finally, with the continuous occurrence of collisions and the generation of broken annular cracks, the ice fragments gradually climb, slide, accumulate, and clear on the ice front until equilibrium. At this point, the cone completely invades the ice sheet. Therefore, in the process of the fully ice–OWT interaction, the dynamic ice loads on the ice-breaking cone can be divided into the direct action of sea ice on the cone surface and the contact action of accumulated ice fragments.

Bending failure mode of the sea ice model.

Furthermore, according to the degree of intrusion into the ice sheet by the OWT structure, the effective stress distribution of ice is illustrated in Figure 11 . As indicated in the figure, in the early stage of contact collision, cracks are initiated and the fractured ice elements begin to separate from the ice sheet. Then, as the OWT structure intrudes, the stress wave exhibits a circular spread.

Ice stress distribution at different time steps. (A) at 7.4 s. (B) at 17.7 s. (C) at 28.1 s.

 4.2.2 Ice force verification

The numerical calculation results of the ice force in the presented study are compared with the three-dimensional horizontal and vertical static ice force calculation models of the cone structure of Ralston (1980) and Croasdale and Cammaert (1994) to verify the accuracy of the numerical simulation. The static ice force calculation models are expressed as

F H R = A 4 [ A 1 σ f h 2 + A 2 ρ w g h D 2 + A 3 ρ w g h ( D 2 − D T 2 ) ] F V R =B 1 F H + B 2 ρ w g h ( D 2 − D T 2 ) F H C = C 1 D σ f ( ρ w g h 5 / E ) 0.25 + C 2 h R h D ρ i g F V C = F H / ς

where the superscripts R and C represent the Ralston and Croasdale ice force models, respectively; FH and FV are the horizontal and vertical loads, respectively; D is the diameter of the cone at the waterline; DT is the diameter of the top of the cone; ρi is the sea ice density; h is the sea ice thickness; hR is the ice climbing height; E is the elastic modulus of the sea ice; A 1, A 2 are dimensionless coefficients depending on the ice thickness and flexural strength; and A 3, A 4, B 1, B 2, C 1, and C 2 are the dimensionless coefficients that depend on the cone angle and the ice-cone friction coefficient, respectively.

Figure 12 shows the dynamic ice force time history of the collision between an ice sheet with an ice velocity of 0.4 m/s, an ice thickness of 0.4 m, and an ice-breaking cone with a cone angle of the 60°. The dynamic ice force with obvious fluctuations can be distinctly observed from the figure. This is because the ice floes have been going through a cyclic process of "contact deformation–fracture failure–accumulated climbing–ice fragment removal–contact deformation". Under the compressive stress and bending stress, the amplitudes of ice force gradually increases. Subsequently, the ice loads are unloaded due to the failure deletion of the cohesive elements. Furthermore, it can be seen from Figure 12A that the simulation results of the presented study and the horizontal maximum ice force values calculated by the Ralston model and the Croasdale model are 1.360 MN, 1.427 MN, and 0.829 MN, respectively. The certain difference in the calculation results of the two model formulas can be discovered from the figure, and the Croasdale model has the smallest calculation result for ice force. Moreover, simulation results of the horizontal ice force amplitude in the presented study are between the Ralston model and the Croasdale model, and the relative error with the Ralston model is the smallest, which is only 15.35%. However, for LC 2.3 of the combined load case with a cone angle of 75°, the maximum error of the horizontal ice force reaches 92.86%, which indicates that the static ice force model used in the study of ice-induced vibration of OWTs with large cone angles shall significantly overestimate the ice force results. This difference is mainly due to the discrepancies between the selection of physical and mechanical parameters of sea ice in the coupling interaction model of ice–OWT and the theoretical models for ice force. In addition, the pile–soil interaction is a factor not considered in the theoretical models, and the influence of the pile–soil interaction cannot be ignored in the dynamic response analysis of OWTs. The calculation of relative error Δ is shown in Eq. (17).

Δ = ( | X n u m e r i c a l − X t h e o r e t i c a l | X n u m e r i c a l ) × 100 %

where X numerical is the numerical simulation result in the presented study and X theoretical is the result calculated by ice force theoretical models.

Comparison of dynamic ice load calculation with CEM, DEM, and static ice force models.

The amplitude of dynamic ice force is not sensitive to the variations of ice speeds observed from the statistics of horizontal and vertical dynamic ice loads in Tables 4 , 5 . In the LC 1.1 load case, the simulated dynamic ice force amplitude is 1.226 MN, and the relative error with the Ralston model is 16.39%. In the load cases of LC 1.2 and LC 1.3, the above relative errors are 15.35% and 4.93%, respectively. On the other hand, as the cone angle increases, the ice force amplitude increases significantly. Taking the 60° cone angle as an example, the horizontal ice force amplitude when the 75° cone angle is used can obtain 1.681 MN, which is significantly larger than the 60° cone angle condition of LC 2.2. Furthermore, as noted from Table 5 , consistent with the variation law of the horizontal dynamic ice loads, the vertical dynamic ice loads increase slightly with the increase of the ice speed. However, with the increase in cone angles, the dynamic ice load amplitude decreases, which is consistent with the results calculated by the Croasdale model. It is demonstrated that the more obvious reduction of the horizontal dynamic ice force of the ice-breaking cone can be caused by the smaller cone angles. However, it shall lead to the increase of the vertical ice dynamic force, which is prominently unfavorable for the ice-breaking cone.

Comparisons of horizontal ice load statistics.

Load cases Maximum value/MN Relative error Δ
Simulation results Ralston model Croasdale model Ralston model Croasdale model
LC 1.1 1.226 1.427 0.829 16.39% 32.38%
LC 1.2 1.237 1.427 0.829 15.35% 32.98%
LC 1.3 1.360 1.427 0.829 4.93% 39.04%
LC 2.1 1.064 0.692 0.762 34.96% 28.38%
LC 2.2 1.237 1.427 0.829 15.35% 32.98%
LC 2.3 1.681 3.242 2.976 92.86% 77.03%

Comparisons of vertical ice load statistics.

Load cases Maximum value/MN Relative error Δ
Simulation results Ralston model Croasdale model Ralston model Croasdale model
LC 1.1 0.596 0.650 0.280 9.06% 53.02%
LC 1.2 0.642 0.650 0.280 1.24% 56.38%
LC 1.3 0.669 0.650 0.280 2.84% 58.14%
LC 2.1 0.999 0.500 0.507 50.45% 49.24%
LC 2.2 0.642 0.650 0.280 1.24% 56.38%
LC 2.3 0.442 0.653 0.193 47.73% 56.33%

Furthermore, the dynamic ice force comparison results on the basis of DEM can be found in Figure 12 . For the bending failure of ice, when applying the DEM to dynamic ice load calculations, the magnitude of dynamic ice force is less than that of CEM. For example, the maximum ice force based on DEM and CEM is 0.753 MN and 1.237 MN, respectively. In addition, by comparing the DEM results, static ice force calculation models are more conservative. This is mainly because there are many limitations and differences in the parameter selection of sea ice and the OWT structure that affect the basic assumptions of sea ice models and the results of numerical simulations. Meanwhile, the calculation parameters of DEM, such as the number of particle layers, stiffness, and bond strength, also need to be further determined. Moreover, since sea ice failure relies on small-sized cohesive spherical particles in DEM, differences in ice failure paths and crack lengths shall also lead to discrepancies of dynamic ice forces.

 4.3 Ice-induced vibration of OWT structures with cones

Figures 13A, B show the time-history curves of the displacement response at the tower top under different ice velocities and cone angles, respectively. From Figure 13A , it can be found that the ice speed has little effect on the structural motion response of the OWT, especially for the operating status of the OWT in this study. In addition, relative to the ice velocity, the obvious differences are depicted in the displacement time-history curves under the three cone angle conditions, and the difference in the above displacement amplitudes reaches 0.13 m, as shown in Figure 13B . In order to further evaluate OWT vibrations under ice and wind loads, tower top acceleration statistics are listed in Table 6 . As noted from the table, the same results can also be observed; that is, ice speeds have little effect on acceleration amplitude. However, as the cone angle increases, the acceleration amplitude at the tower top increases. For example, the acceleration amplitude of the OWT with 75° cone can exceed 2.4 m/s2. Moreover, the acceleration fluctuations are more consistent under different ice speeds and cone angles by comparing the standard deviation.

Histories of tower top displacements under selected wind-ice cases. (A) Ice speeds (B) Cone angles.

Comparisons of tower top acceleration statistics.

Load cases Maximum value (m/s2) Standard deviation (m/s2)
LC 1.1 2.32 1.33
LC 1.2 2.36 1.31
LC 1.3 2.39 1.32
LC 2.1 2.28 1.31
LC 2.2 2.36 1.31
LC 2.3 2.47 1.30

To further illustrate the discrepancies in the vibration characteristics at the tower top of the OWT, the frequency domain displacement responses at the tower top under the corresponding sea ice conditions are shown in Figure 14 . In the frequency domain, the tower top displacement is mainly affected by the fundamental frequency of the overall structure of 0.246 Hz. It should be noted that since the quality variations of the ice-breaking cone are induced by the cone angles in Figure 9 , the fundamental frequencies of the OWT with three ice-breaking cones are slightly different, which are 0.247 Hz, 0.246 Hz, and 0.245 Hz, respectively, as shown in Figure 14B . In addition, the increase in the weight of the ice-breaking cone shall cause the fundamental frequency of the OWT to decrease, which shall lead to the OWT structure resonance more easily. Meanwhile, considering the influence of the cone angle on the dynamic ice loads in the horizontal and vertical directions, it is recommended to select the cone angle as 60°.

Fourier amplitudes of tower top displacements under selected wind-ice cases. (A) Ice speeds (B) Cone angles.

 4.4 Damage analysis of the ice-breaking cone

Based on the CEM, the Von-Mises effective stress and pressure of the conical OWT structure under the sea ice can be obtained through the established ice–OWT interaction coupling model, as shown in Figure 15 . The annular spread of stress at the collision position can be observed. Moreover, the remarkable increasing trend of damage distribution with the increasing of the cone angles is indicated in the figure. The main reason for such result is that the increase of the cone angle leads to a larger contact area between the ice and the cone surface. Concurrently, it can be seen from Figure 16 that the stress distribution corresponds to the pressure distribution, and the maximum pressure is the region with the maximum Von-Mises effective stress.

Conical structure Von-Mises effective stress distribution. (A) 45° cone angle (B) 60° cone angle (C) 75° cone angle.

Pressure distribution of cone structure. (A) 45° cone angle (B) 60° cone angle (C) 75° cone angle.

The variation trends of Von-Mises effective stress and pressure with different cone angles are shown in Figure 17 . It can be seen from the comparison of the figures that the cone angle has a great influence on the stress. With the increase of the cone angle, the Von-Mises effective stress and pressure show an increasing trend, which is also consistent with the dynamic ice loads and the dynamic response of the OWT structure. Furthermore, the comparison results show that the maximum Von-Mises effective stress is 186 MPa, which does not reach the yield strength of the steel material of 355 MPa; that is, the OWT structure is in the elastic deformation stage during the ice–OWT interaction.

Maximum Von-Mises effective stress of pyramidal structures under different taper angles.

 Conclusion

In this paper, nonlinear distributed springs are used to simulate the pile–soil interaction of OWTs, and based on the CEM, an ice–OWT interaction model with ice-breaking cones is established. Furthermore, based on the secondary development of LS-DYNA, the Tsai-Wu yield criterion and empirical failure formula are used to describe the failure and failure modes of sea ice, and the bending failure process of sea ice is simulated. According to the established interaction coupling numerical simulation model, the analysis on the motion response and strength of the OWT with cone under the combined wind and ice are carried out. The main conclusions of the research are as follows:

Based on the CEM, the cohesive element model is combined with the finite element constitutive model to carry out numerical simulation research on the coupled interaction between ice and OWTs. Compared with the Ralston ice force model and the Croasdale ice force model, the result shows that the adopted method in the presented study can simulate the bending failure behavior of sea ice well, especially regarding the failure of ice elements and the accumulation of ice fragments consistent with actual field conditions.

With the increase in ice velocity, the motion response of OWTs is increasing. However, the ice velocity has little effect on the dynamic ice loads and the structural motion response of the OWTs, especially for the operating OWT with the pile–soil interaction.

The ice-breaking cone angle has a significant impact on the dynamic ice loads and the structural motion response of OWTs. As the cone angle increases, the magnitude of the horizontal ice load increases significantly; however, for the vertical ice load, the magnitude shows a decreasing trend. This is mainly because variations in the cone angle change the pattern of sea ice destruction. In addition, the increased weight of the ice-breaking cone may lead to the resonance of OWT. Therefore, considering the pile–soil interaction, it is suggested that the optimal ice-breaking cone angle of the OWT under the operating state is 60°.

The remarkable discrepancies in the damage of the ice-breaking cone surface with different cone angles can be directly observed. The Von-Mises effective stress distribution is consistent with pressure distribution. With the increase of the cone angles, the stress damage area and effective stress amplitude are also positively correlated. However, the maximum Von-Mises effective stress is much smaller than the yield strength of the steel material. This indicates that the linear analysis of the OWT structure in the elastic range is also reasonable for the ice-conical OWT structure interaction process.

 Data availability statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

 Author contributions

Conceptualization, WS, JZ and XL. Methodology, BW, YZL, and YL. Investigation, BW and YZL. Writing—first draft preparation, BW, YZL, JZ, and WS. Review and editing, WS and XL. Supervision, WS and XL. Project administration, JZ. Funding acquisition, WS and XL. All authors contributed to the article and approved the submitted version.

 Funding

This research was supported by Key Laboratory of far-shore wind power technology of Zhejiang Province (ZOE2020003), Zhejiang Provincial Natural Science Foundation of China (LHY21E090001) and the National Natural Science Foundation of China (Grant Nos. 52001052, 51939002, 52071058, 52071301 and 51909238). And partially supported by Ministry of Industry and Information Technology (2020) g-275.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Croasdale K. R. Cammaert A. B. (1984). An improved method for the calculation of ice loads on sloping structures in first-year ice. Hydrotech Constr. 28 (3), 174–179. doi: 10.1007/BF01545935 Dempsey J. P. Adamson R. M. Mulmule S. V. (1999). Scale effects on the in situ tensile strength and fracture of ice. part II: first-year sea ice at resolute, NWT. Intentional J. Fract 95 (1-4), 347–366. doi: 10.1023/a:1018650303385 Derradji-Aouat A. (2003). Multi-surface failure criterion for saline ice in the brittle regime. Cold Reg. Sci. Tech 36 (1-3), 470–470. doi: 10.1016/S0165-232X(02)00093-9 Di S. Xue Y. Wang Q. (2017). Discrete element simulation of ice loads on narrow conical structures. Ocean Eng. 146 (1), 282–297. doi: 10.1016/j.oceaneng.2017.09.033 DNV GL (2016). DNV GL-ST-0437: Loads and site conditions for wind turbines (Oslo, Norway: DNV GL). Gladman B. (2013). LS-DYNA keyword user's manual, volume Ⅱ, material models (Livermore, California: Livermore Software Technology Corporation). Hallquist J. (2006). LS-DYNA® theory manual. Livermore (Livermore, California: Livermore Software Technology Corporation), 3. Heinonen J. Rissanen S. (2017). Coupled-crushing analysis of a sea ice-wind turbine interaction-feasibility study of FAST simulation software. Ships Offshore Struct. 12 (7-8), 1–8. doi: 10.1080/17445302.2017.1308782 Hillervborg A. Mod E. M. Petersson P. E. (1976). Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res. 6 (6), 773–781. doi: 10.1016/0008-8846(76)90007-7 IEC (2012). IEC 61400-1: Wind turbine-part 1: Design requirements (Geneva, Switzerland: IEC). Jang H. K. Kim M. H. (2021). Dynamic ice force estimation on a conical structure by discrete element method. Int. J. Nav Archit Ocean Eng. 13 (1), 136–146. doi: 10.1016/j.ijnaoe.2021.01.003 Kilcher J. (2012). TurbSim user's guide (Draft version) (Golden, USA: National Renewable Energy Laboratory, Golden). Lan J. Y. Liu H. S. Lv Y. J. (2012). Dynamic nonliner parameters of soil in the bohai sea and their rationality. J. Harbin Eng. Univ. 33 (9), 1079–1085. doi: 10.3969/j.issn.1007-7043.201202036 Liu Y. Z. Shi W. Wang W. H. (2022). Dynamic analysis of monopile-type offshore wind turbine under Sea ice coupling with fluid-structure interaction. Front. Mar. Sci. 9, 839897. doi: 10.3389/fmars.2022.839897 Ralston T. D. (1980). Plastic limit analysis of sheet ice loads on conical structures, physics and mechanics of ice, 289–308. Ren Y. Vengatesan V. Shi W. (2022). Dynamic analysis of a multi-column TLP floating offshore wind turbine with tendon failure scenarios. Ocean Eng. 245, 110472. doi: 10.1016/j.oceaneng.2021.110472 Shi W. Tan X. Gao Z. (2016). Numerical study of ice-induced loads and responses of a monopile-type offshore wind turbine in parked and operating conditions. Cold Reg. Sci. Tech 123, 121–139. doi: 10.1016/j.coldregions.2015.12.007 Song M. Shi W. Ren A. R. Zhou L. (2019). Numerical study of the interaction between level ice and wind turbine tower for estimation of ice crushing loads on structure. J. Mar. Sci. Eng. 7 (12), 439. doi: 10.3390/jmse7120439 Statoilhydro A. G. Bjerks M. Kühnlein W. (2009). In: Proceedings of the ASME 28th International Conference on Ocean, Offshore and Arctic Engineering, Honolulu, Hawaii, USA, May 31 - June 5. doi: 10.1115/OMAE2009-80164 Timco G. W. Weeks W. F. (2010). A review of the engineering properties of sea ice. Cold Reg. Sci. Tech 60 (2), 107–129. doi: 10.1016/j.coldregions.2009.10.003 Velarde J. Bachynski E. E. (2017). Design and fatigue analysis of monopile foundations to support the DTU 10 MW offshore wind turbine. Energy Proc. 137, 3–13. doi: 10.1016/j.egypro.2017.10.330 Wang Q. (2015). Ice-inducted vibrations under continuous brittle crushing for an offshore wind turbine (Trondheim, Norway: Norwegian University of Science and Technology). Wang Y. P. Shi W. Michailides C. Wan L. Kim H. Li X. (2022). WEC shape effect on the motion response and power performance of a combined wind-wave energy converter. Ocean Eng. 250, 111038. doi: 10.1016/j.oceaneng.2022.111038 Wang Z. L. Wang Y. G. Han Q. Y. (1980). Visco-elastoplastic soil model for irregular shear cyclic dynamic loadings. Chin. J. Geotech Eng. 2 (3), 10–20. doi: CNKI:SUN:YTGC.0.1980-03-001 Yang Y. Bashir M. Li C. (2020). Mitigation of coupled wind-wave-earthquake responses of a 10 MW fixed-bottom offshore wind turbine. Renew Energy 157, 1171–1184. doi: 10.1016/j.renene.2020.05.077 Ye K. H. Li C. Yang Y. (2019). Research on influence of ice-induced vibration on offshore wind turbines. J. Renew Sustain Energy 11 (3), 033301. doi: 10.1063/1.5079302 Zhan K. Y. Cao L. S. Wan D. C. (2021). Cohesive element method for ice load on cortical structures. Ocean Eng. 39 (4), 8. doi: 10.16483/j.issn.1005-9865.2021.04.007 Zhang L. X. Shi W. Karimirad M. Michailides C. Jiang Z. Y. (2020). Second-order hydrodynamic effects on the response of three semisubmersible floating offshore wind turbines. Ocean Eng. 207, 107371. doi: 10.1016/j.oceaneng.2020.107371 Zhao Z. X. Shi W. Wang W. H. Qi S. Li X. (2021). Dynamic analysis of a novel semi-submersible platform for a 10 MW wind turbine in intermediate water depth. Ocean Eng. 237, 109688. doi: 10.1016/j.oceaneng.2021.109688