Discrete element simulation of ice loads on narrow conical structures

Ocean Engineering 146 (2017) 282–297

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Discrete element simulation of ice loads on narrow conical structures Shaocheng Di *, Yanzhuo Xue, Qing Wang, Xiaolong Bai Harbin Engineering University, College of Shipbuilding Engineering, 150001 Harbin, China

A R T I C L E I N F O

A B S T R A C T

Keywords: Ice loads Narrow conical structure Discrete element method Numerical model

This paper investigates a numerical approach to predict the ice loads on a conical structure through simulating the failure process of level ice sheet using the discrete element method (DEM). In the simulation, the level ice sheet is modeled with the bonded spherical particles, which can be de-bonded based on the bonding-breaking criteria. The ice sample compression and three-point bending tests are used to set up the relationship between the ice strengths, i.e. compression strength and flexural strength, and the bonding strengths in the DEM model; the sensitivity of particle size, normal and shear bonding strength, normal and shear contact stiffness are analyzed. The ice loads on narrow conical structures with different cone angles are simulated using the DEM model. Meanwhile, a comparison between numerical results and the ISO ice force standard is carried out. In addition, the interaction between level ice and multi-legs conical structure is simulated to study the shadowing effect.

1. Introduction The failure of sea ice cover against a conical structure is a fragmentation process where a continuous ice sheet is broken into discrete blocks mainly under the bending failure. For designing of offshore structures in ice covered waters, it is important to well understand this process and properly estimate ice loads. The ice load measurements on full-scale structures can obtain the significant accurate ice force data. There have been scarce published full-scale ice load data measured on conical structures around the world. Some well-known tests were the Finnish Kemi-I test in the Gulf of Bothnia (M€a€att€anen et al., 1996; Brown and M€ a€ att€ anen, 2009), the Confederation Bridge ice load measurements in the Southern Gulf of Lawrence (Brown et al., 2010) and the ice forces measurements on a jacket platform equipped with cones in the JZ20-2 field in the Bohai Sea of China (Yue et al., 2007). Compared with the full-scale tests of the ice loads on conical structures, the model scale tests conducted in ice tanks have some advantages of studying the influences of various factors (ice thickness, ice speed, water line diameter, etc.) on the ice loads of conical structures (Huang, 2010; Huang et al., 2013; Dalane, 2014). Meanwhile, there were some numerical models developed to simulate the failure process between ice sheet and conical structures and to calculate the ice loads acting on the conical structures. Matlock simulated the interaction between ice sheet and conical structure in 1971 (Matlock et al., 1971). In this model, he described structures and ice with a spring-mass-damper system and moving bars separately. The “Matlock-model” is able to simulate the ductile and

brittle ice crushing, corresponding to low and high velocities of moving ice, respectively (Withalm and Hoffmann, 2010). Although it is available for the structures with a high ratio between structure width and ice thickness, the ice sheet is divided into several zones by cracks instead of fails as one piece across the structure. Based on this, Eranti (1991) presented a zonal model of ice interacting with a cylindrical structure, and introduced two fracture lengths in ice advancing direction and radial direction to estimate the maximum ice force. This method is practical for modeling small-scale configurations and has a highly agreement with model-scale measurements. However, it is difficult to determine the fracture length in the Eranti's model which is only available for conical structures. Therefore, it is crucial to investigate some more effective and accurate numerical methods to simulate the interaction processes between sea ice and offshore structures with complex geometries. DEM method has been used in the simulation of pancake ice floes on cylindrical piles by Sun and Shen (2012) (Sun and Shen, 2012), and rubble pile formation process against a wide inclined structure by Hopkins (1997) and Paavilainen et al. (2011) (Hopkins, 1997; Paavilainen et al. 2011). The combined finite-discrete element method is also an effective approach to model ice sheet and its failure with FEM and subsequent pile-up process with DEM (Poloj€arvi and Tuhkuri, 2009). This paper examines ice sheet failure process against a narrow conical structure with the discrete element method by treating the ice as elasticbrittle material. The failure mode of ice is strongly dependent on the ice conditions, structure geometry and the relative velocity between the ice

* Corresponding author. E-mail address: [email protected] (S. Di). https://doi.org/10.1016/j.oceaneng.2017.09.033 Received 10 February 2017; Received in revised form 26 August 2017; Accepted 24 September 2017 0029-8018/© 2017 Elsevier Ltd. All rights reserved.

S. Di et al.

Ocean Engineering 146 (2017) 282–297

where F nb , F sb and M nb , M sb denote the normal- and tangential-direction forces and moments, respectively. Each subsequent relative displacement and rotation increment produces an increment of elastic force and moment. The increments of elastic force and moment acting on the bond can be written as follows:

and structure. The simulation results have high agreement with the measurements on a conical structure in the Bohai sea (Yue and Bi, 2000; Yue et al., 2007; Xu et al., 2011) as the purpose of deploying a conical structure is to reduce the ice loads on the flexible vertical structure (Qu et al., 2006). The peak ice loads are analyzed as well against conical structures with different slope angles and diameters. Firstly, this paper introduces the discrete element model for sea ice. Secondly, the parameter sensitivity analysis of this numerical method, mainly including the choice of micro-parameters, is presented. Numerical tests of uniaxial compressive and three-point bending of sea ice are performed to determine the stiffness and strength parameters in the failure criterion of the bonding model. Finally, the influence of the particle diameter on ice mechanical properties is investigated. With this model, simulations are conducted to study the failure process of level ice advancing against a conical structure. The results are verified with fullscale measurements from the JZ20-2 MUQ platform in Bohai Bay. Furthermore, the maximum ice loads on offshore structures with different slope angles and diameters are studied and compared with the ISO ice force standard.

The ice sheet is represented by a dense packing of uniform-sized spherical particles that are bonded together at their contact points with parallel bonds. Its mechanical behavior is simulated using the threedimensional discrete element method. In this method, the ice sheet and its fracture are modeled by the beam theory. A comprehensive presentation of the method can be found in the work of Potyondy and Cundall (2004) (Potyondy and Cundall, 2004). The mechanical behavior of an elastic-brittle bond joining the two bonded particles is approximated by a parallel bond. Fig. 1 shows the sketch of a parallel bond, where xA and xB are the position vectors of element A and B, respectively. The parallel bond can be envisioned as a set of elastic springs uniformly distributed over a circular cross-section lying on the contact plane and centered at the contact point. Parallel bond can transmit both force and moment between particles. The total force and moment carried by the parallel bond are denoted by F b and M b , respectively, which represent the action of the bond on particle B. The force and moment vectors can be resolved into normal and tangential components as

(1)

M b ¼ M nb þ M sb

(2)

(3)

ΔFbs ¼ k s ΔU s

(4)

ΔMbn ¼ k s JΔθn

(5)

ΔMbs ¼ k n IΔθs

(6)

where A (¼ πR2 ), I (¼ 14 πR4 ) and J (¼ 12 πR4 ) are the area, moment of inertia and polar moment of inertia of the parallel bond cross-section, respectively. kn (Nm3 ) and ks (Nm3 ) are the bond stiffness in the normal and tangential direction. ΔU n and ΔU s are the relative normal and tangential displacement increments. Δθn and Δθs are the relative normal and tangential rotation increments. The maximum tensile and shear stress acting on the parallel-bond periphery are derived from beam theory (Eqs. (7) and (8))

2. DEM model of sea ice

F b ¼ F nb þ F sb

ΔFbn ¼ k n ΔU n

σ max ¼

τmax ¼


n

F
b

A

þ


s

F
b

A


s

M
b

I

R


n

M
b

þ

J

R

(7)

(8)

The bonds break when either the maximum tensile stress exceeds the tensile strength (σ max > σ nb ) or the maximum shear stress exceeds the shear strength (τmax > σ sb ). σ nb is the inter-particle normal bonding strength, and σ sb the inter-particle shear bonding strength. The failure criterion of bond between two particles is indicated as follows (Eq. (9)) (Ji et al., 2016)

τb ¼ σ sb þ μb σ n

(9)

In this failure criterion, the shear strength τb is determined by the bonding strength σ sb in the shear direction and the friction induced by the normal stress σ n following the Mohr-Coulomb law, and a sliding friction coefficient μb between the two bonded particles also be introduced. The detail of the failure criterion and the determination of the values of σ nb , σ sb and μb can be found in the previous work (Ji et al., 2016). 3. Contact detection between ice particle and cone A narrow conical structure operating in Bohai Bay is selected for the calculation examples in this paper. The jacket leg is shown in Fig. 2. Fig. 3 shows the DEM model of one of the jacket's legs with ice-breaking cone.

Fig. 1. Parallel bonding model between two spherical particles.

Fig. 2. Ice-breaking cone on jacket leg in Bohai Bay. 283

S. Di et al.

Ocean Engineering 146 (2017) 282–297

Fig. 3. DEM model of the jacket leg equipped an ice-breaking cone.

3.1. Contact mode The jacket leg consists of the conical column members. The contacts between the spherical particle and the cone have three types, as shown in Fig. 4. (a) The interaction with the rotary surface of cone; (b) The interaction with the bottom/top surface; (c) The interaction with the edges. In order to judge the contact between ice particle and cone, some specific points are first defined: the circle centers of two base surfaces, A and B, radiuses, R1 and R2 , position and radius of spherical particle, P and Rp . The projection of the point P at the axis of the cone is Q, as shown in Fig. 4 (a and b). The below vectors connecting those points are then acquired, AP, BP, AB, BA, BQ, QP. Where BQ ¼ BP  QP, AB QP ¼ AP  AP⋅AB jABj ⋅jABj. These vectors are also acquired, AM, BN, MN, NM, QP QP , BN ¼ r2 jQPj , MN ¼ AB þ BN  AM, MP, NP, RP. Where AM ¼ r1 jQPj MN NM ¼ MN, RP ¼ MP  MP⋅MN jMNj ⋅jMNj.

The contact detection for the spherical particle and the cone include these ways: (a) The contact between the spherical particle and the rotary surface of cone The center of the spherical particle is first checked as to whether it is in the area A1, as shown in Fig. 4(a). If ðMP⋅NMÞðNP⋅NMÞ  0, then the particle is in the area A1 and the contact between the particle and the cone is further checked. If jRPj < Rp , then the particle contacts with the rotary surface of the cone. (b) The contact between the spherical particle and the base surfaces. If ðMP⋅NMÞðNP⋅NMÞ > 0, then the spherical particle is in the area A2, A3, A4, A5. If jQPj  r1 and jAQj < jBQj, then the spherical particle is in the area A2, as shown in Fig. 4(b), and the contact between particle and cone is further checked. For the base surface with radius r1, if jAQj < Rp , then the particle contacts with the bottom surface. If the spherical particle is not in the area A1, A2, A4, and jQPj  r2 and jBQj < jAQj, then the particle is in the area A3, as shown in Fig. 4(b). If jBQj < Rp , the spherical particle contacts with the top surface.

Fig. 4. Schematic of contact detection between sea ice particle and ice breaking cone; there are three contact modes: (a) Rotary surface contact; (b) Base surface contact; (c) Edge contact; P is the position of spherical particle, the shaded areas A1~A5 in subfigure (a)~(c) are the potential contact regions between spherical particle and cone.

contacts with the top edge. As shown in Fig. 4(c).

(c) The contact between the spherical particle and the edges.

(a) Rotary surface contact mode (b) Base surface contact mode (c) Edge contact mode

If jQPj > r1 or jAQj > jBQj, then the spherical particle is in the area A4, if jMPj < Rp , then the particle contacts with the bottom edge. If jQPj > r2 or jBQj > jAQj, then the particle is in the area A5. If jNPj < Rp , the particle

284

S. Di et al.

Ocean Engineering 146 (2017) 282–297

F s ←F s þ kws Δxs ;

3.2. Contact force

 

 F s ← F s μw F n
F s
;

Having identified the type and location of a potential contact point the contact detection algorithm finishes with calculation of the amount of overlap between spherical particle and the cone





ΔL ¼
PC P
 Rp

(10)





for sliding with
F s

μw F n


(16) (17)

4. Micro-parameter sensibility analysis in DEM simulation Parameter sensibility analysis is used to determine the correct numerical representation for a given specimen. The determination of parameters is vital for a numerical model. To establish the microparameters, a series of uniaxial compressive and three-point bending tests are performed on this model.

(11)

The force acting at the cone can be calculated by employing a contact law based on the relationship between the forces and incremental relative displacements of the spherical particle and the cone. The relative displacement vector Δx of the particle and cone is

  Δx ¼ vp  vw Δt



F s <
μ w F n


where ksw is the tangential contact stiffness, and μw is the Coulomb friction coefficient. The acceleration, velocity and position of each particle can be calculated by Newton's second law.

where PC is the contact point on the cone, and it denotes R (for rotary surface contact), Q (for base surface contact), and P (for edge contact), respectively. Finally, the contact unit vector nw pointing from the centroid of the spherical particle to the contact point on the cone is

PC P
nw ¼



PC P


for no sliding with

4.1. Numerical tests The sketches of the uniaxial compressive and flexural tests are shown in Figs. 5(a) and 6(a). Figs. 5(b) and 6(b) show that a packing pattern of regular hexagonal lattice was used to construct the compressive and flexural specimen where spherical particles represent the basic elements.

(12)

where, vp is the velocity of the spherical particle, and vw is the velocity of the contact point on the cone. The relative displacement vector can be decomposed into normal and tangential components using the contact unit vector nw

Δxn ¼ ðnw ΔxÞnw

(13)

Δxs ¼ Δx  Δxn

(14)

where Δxn is the normal component of the relative displacement and Δxs is the tangential component of the relative displacement between the spherical particle and the cone. A linear contact force model is employed to calculate the normal and tangential contact force, and the normal force can be expressed as

F n ¼ kwn ΔL

(15)

where knw is the normal contact stiffness which is function of radius and elastic modulus of the spherical particle. The tangential contact forces are

Fig. 6. Demonstration of the three-point bending test of sea ice: (a) schematic diagram of sample; (b) numerical sample in DEM simulation.

Fig. 5. Demonstration of the uniaxial compressive test of sea ice: (a) schematic diagram of sample; (b) numerical sample in DEM simulation. 285

S. Di et al.

Ocean Engineering 146 (2017) 282–297

Table 1 Computational parameters matrix. Definitions

Symbols

Values

Density of sea ice Diameter of particles Normal contact stiffness Shear contact stiffness Inter-particle friction coefficient Normal bonding strength Shear bonding strength Number of particles Loading rate

ρ D kn ks μb σ nb σ sb NP vL

920 kg=m3 10 mm 1:0  1011 Nm3 5:0  1010 Nm3 0.2 0.73 MPa 0.73 MPa 3300(Compression)/5040 (flexural) 0.02 m/s

The sample sizes of sea ice are a  a  h ¼ 100 mm  100 mm  250 mm for the compressive specimen and b  h  l ¼ 70 mm  70 mm  500 mm for the flexural specimen. The loading rate is 0.02 m/s in both compressive and flexural tests. Table 1 lists the main computational parameters in the simulation of compressive and flexural numerical tests. The failure processes of sea ice sample simulated with DEM under uniaxial compressive and three-point bending tests were indicated in Reference (Ji et al., 2016). From the stress-strain curve of the compressive test plotted in Fig. 7, The maximum normal stress can be obtained as the uniaxial compressive strength. For the three-point bending test, the maximum stress can be obtained as the flexural strength, as shown in Fig. 8. For the two tests simulated with DEM described above, the uniaxial compressive strength σ c ¼ 5.64 MPa, and the flexural strength σ f ¼ 1.25 MPa. In Table 1, the computational parameters were given roughly to gain the generic numerical results. The selection of the parameters reasonably will be indicated in the next section.

The scaling laws connecting the macro- and micro-models are determined by invoking dimensional analysis. Since there are eight parameters and three independent dimensions (e.g., [F], [T] and [L]), so the mechanical responses of the element assembly of sea ice are governed by five independent dimensionless parameters, and can be expressed as

4.2. Parameter sensibility analysis

Rs ¼ f

Fig. 8. Maximum flexural stress versus vertical deflection in three-point bending failure process of sea ice simulated with DEM; the maximum stress is the flexural strength.

  Rs ¼ f L; D; ρ; k n ; ks ; μb ; σ nb ; σ sb ; vL

(18)

In these parameters, the L is the specimen dimension of sea ice, and ρ is the density of sea ice. These two parameters can be obtained directly. So equation above can be changed to

  Rs ¼ f L; D; k n ; k s ; μb ; σ nb ; σ sb ; vL

A discrete element model of sea ice with parallel bonds is characterized by two groups of parameters at the micro-scale: geometric and physical parameters (the diameter of element D, the size of specimen L (The value here is just used to illustrate the influence factor, not referring to a certain size), the density of sea ice ρ ) and constitutive parameters for the contacts (the normal contact stiffness kn , the tangential contact stiffness ks , the normal bonding strength σ nb , the shear bonding strength σ sb , the friction coefficient μb ). At the macro-scale, however, a material such as sea ice is described by phenomenological constants such as the elastic constants E, the compressive and flexural strengths σ c and σ f . Since there is no complete theory to predict the macro-scale properties from the above micro-scale parameters, establishing the relationship between the parameters at the two scales is prerequisite for performing numerical simulations with DEM (Huang, 1999; Yang et al., 2006). The macro-scale response Rs can be expressed using these parameters as follows

! D ks σs vL ; n ; μb ; bn ; pffiffiffiffiffiffiffiffiffi L k σ b k n =ρ

(19)

(20)

If the numerical simulations are performed with enough degrees of freedom (D=L≪1) and under the condition of quasi-static loading pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi (vL = kn =ρ≪1), the ratios D=L and vL = kn =ρ can be dropped from the pffiffiffiffiffiffiffiffiffiffi expression above. Note that kn =ρ is dependent on the speed of the compressive wave propagating in an elastic medium. The loading velocity is very small compared with the speed of elastic wave, and it can maintain the loading process to be quasi-static. So the parameter pffiffiffiffiffiffiffiffiffiffi vL = kn =ρ is no longer considered in the expression (20). However, it is difficult to satisfy the condition that the diameter of spherical particle D is much smaller than the size of specimen L (i.e. D=L≪1). The influence of spherical particle D on macro-mechanical properties of sea ice is discussed in the next section. So far, the mechanical responses of sea ice can be reduced to the following expression:

 Rs ¼ f

D ks σs ; n ; μb ; bn L k σb

 (21)

In the following sections, the influence of these four dimensionless parameters D=L, ks =kn , μb and σ sb =σ nb on the macro-mechanical properties will be further discussed. 4.2.1. The influence of D=L To study the influence of spherical particle size on the mechanical properties of sea ice, the uniaxial compression and three-point bending Table 2 Specimen dimension matrix. Uniaxial compressive test (a  a  h) (mm3) 70  70  175, 100  100  250, 150  150  375 200  200  500, 1000  1000  2 500, 2000  2000  5000 Three-point bending test (b  h  l) (mm3)

Fig. 7. Axial stress versus axial strain in compressive failure of sea ice simulated with DEM; the slope of line segment is the elastic modulus of the sea ice, and the maximum stress is the uniaxial compressive strength.

37.5  37.5  250, 75  75  500, 112.5  112.5  750, 150  150  1000

286

S. Di et al.

Ocean Engineering 146 (2017) 282–297

tests are simulated with DEM using various specimen and particle sizes. Here, as listed in Table 2, we have six different specimen sizes for the uniaxial compressive tests, and four for the three-point bending tests. For each sea ice specimen, different particle sizes are adopted to determine the elastic modulus E, uniaxial compressive strength σ c and flexural strength σ f , respectively. The elastic modulus E, the uniaxial compressive strength σ c , and the flexural strength σ f of sea ice under various specimen and particle sizes can be obtained with DEM simulations. To analyze the combined influence of specimen size and particle size on the mechanical properties of sea ice comprehensively in DEM simulations, here, we introduce a relative particle size D=L, and set L ¼ a and b (or h) in the uniaxial compressive and three-point bending test, as shown in Figs. 9–11. We can find that the elastic modulus E, the uniaxial compressive strength σ c , and the flexural strength σ f decrease linearly with the increase of the relative particle size D=L, and the fitting equations are as follows

E ¼ 3:993ðD=LÞ þ 1:80

(22)

σ c ¼ 11:36ðD=LÞ þ 5:96

(23)

σ f ¼ 3:47ðD=LÞ þ 1:94

(24)

Fig. 10. Variation of the uniaxial compressive strength σ c with the ratio of D=L simulated with DEM under various specimen sizes.

From Figs. 9–11, we can find that the mechanical properties of sea ice simulated with DEM are dependent on the relative element size D=L and without the influence of specimen size. All of the E, σ c and σ f decrease perfectly linearly with the increase of D=L. 4.2.2. The influence of ks =kn When spherical particles are arranged in the pattern of Hexagonal Close Packing lattice, the elastic modulus and Poisson's ratio can be obtained from analytical solution as follows (Wang and Alonso-Marroquin, 2009)

E¼

ν¼

pffiffiffi 3 2k n ðk n þ k s Þ Rð5k n þ k s Þ   ðk n  k s Þ 6ðk n Þ2 þ 23ðks Þ2 þ 31k n k s 18ðk n Þ3 þ 119ðk n Þ2 k s þ 128k n ðk s Þ2 þ 23ðk s Þ3

Fig. 11. Variation of the flexural strength σ f with the ratio of D=L simulated with DEM under various specimen sizes.

(25)

(26)

where ks and kn are the tangential and normal contact stiffness in the DEM model, respectively. The analytical results of Eqs. (25) and (26) have been tested numerically using a 3D HCP(Hexagonal close-packing) block under uniaxial compression. We change the tangential stiffness ks , and plot the computed E and ν against ks =kn in Fig. 12. For the elastic modulus, the red solid dots correspond to the analytical result from Eq. (25), and the blue solid triangles correspond to the measurements from the numerical simulations. The close agreement between the numerical measurements and the analytical result validates the numerical results

Fig. 12. The variations of elastic modulus and Poisson's ratio with ks =kn .

under the condition of ks =kn  0:5. Similarly, for the Poisson's ratio, numerical measurements and the analytical result are also agreement with each other. The Poisson's ratio ν measured in the field test of the uniaixal compressive test is usually 0.29. So we choose the ratio of ks =kn of 0.1 as the input computational parameter in the subsequent DEM simulation of ice loads on narrow conical structures. 4.2.3. The influence of σ sb =σ nb and μb In the fracture criterion, the micro-scale parameters those can influence the uniaxial compressive and flexural strengths are the friction coefficient μb , normal bonding strength σ nb and shear bonding strength σ sb . The effect of the ratio of σ sb =σ nb on the uniaxial compressive strength σ c and flexural strength σ f are investigated by varying the values of σ nb and (σ sb . The bonding strengths σ nb and σ sb are chosen as: 0.1 MPa, 0.25 MPa, 0.5 MPa, 1.0 MPa, 1.5 MPa, and 2.0 MPa. Here, the friction coefficient μb

Fig. 9. Variation of the elastic modulus E with the ratio of D=L simulated with DEM under various specimen sizes. 287

S. Di et al.

Ocean Engineering 146 (2017) 282–297

varies from 0.0 to 0.3. The results of influence of normal and shear bonding strength on compressive and flexural strength were indicated in Reference (Ji et al., 2016). In general, the ratio of uniaxial compressive to flexural strength obtained from field tests is in the range from 2.5 to 3.0. Therefore, the micro-scale parameters σ sb =σ nb ¼ 1:0 and μb ¼ 0:2 will be used in the following DEM simulations. So far we have discussed the influence of micro-scale parameters (the normal contact stiffness kn , the shear contact stiffness ks , the particle size D, the normal bonding strength σ nb , shear bonding strength σ sb and the friction coefficient μb ) on the macro-scale mechanical properties (the elastic modulus E, uniaxial compressive strength σ c and flexural strength σ f ). In the next section, we will study the DEM numerical simulations of ice loads of level ice on the conical structures based on the parameters calibrated in this section.

Table 3 Input parameters matrix. Definitions

Symbols

Values

Density of sea ice Density of sea water Initial ice cover area Ice thickness Ice velocity Diameter of particles Normal bonding strength Shear bonding strength Number of particles

ρi ρw LW H Vi D σ nb σ sb NP

920 kg=m3 1035 kg=m3 40  18 m2 0.23 m 0.43 m/s 8.735 cm 0.5 MPa 0.5 MPa 326 304

zoomed out to clearly show the ice cover in the DEM simulation. It is shown that the surface of ice cover constructed with bonded particles is not as smooth as that in natural conditions. It is better to model the ice cover with smaller elements. However, due to the limitation of the computational resources, the larger elements are adopted in this study. A 90-s duration of the ice cover and Jacket interaction was simulated in the calculation examples, requiring 27 h running time on the K40 GPU Card workstation. The failure mode of the level ice against a narrow conical structure is simulated with DEM model as shown in Fig. 14. The observation shows that the bending failure mode dominating the simulation, which highly fits the practical scenario as shown in Fig. 2. Fig. 15 also clearly shows the ice fragments behind a breakage. Fig. 16 presents the comparisons of the simulation results and field measurements on ice load time histories. It can be seen that the simulation results agree fairly well with the measurements. The maximum ice load in the ice moving direction in field measurement is 117 kN, and the maximum ice load FH in the horizontal and FV in the vertical direction in the simulation is 109 kN and 70.0 kN, respectively. Here, the “up-down” property of the force-time curve is obvious and the simulation time is relatively short, so the peak ice loads were picked manually from the force-time curve. In the ice advancing direction, the minimum peak ice load is 70.6 kN, and the average value of peak ice loads is 91.2 kN. After statistical analysis, it is obvious that the peak ice loads mainly vary from 75 kN to 110 kN, and the maximum value of peak ice loads is 95 kN. The periods of the peak ice loads are in the range from 1.5 s to 4.5 s, and the period of 2.5s is the dominant. From the field observations, the periods of adjacent peak forces are in the range from 1.4 s to 1.7s. (Qu et al., 2006). Field measurements and DEM simulation results of the ice loads were also analyzed in the frequency domain. The power spectrum density curves of the ice loads obtained from the field data and DEM simulation are plotted in Fig. 17. It is shown that the simulation results agree well with the field test data in terms of the peak value of the power spectrum density. The slight difference of the frequency may be caused by the fragment size of broken ice block. Sea ice is a complex material which exhibits strong discreteness in the basic mechanical properties. We tested many sea ice samples collected from the Bohai Bay in China. The

5. DEM simulation of ice loads on conical structures 5.1. Numerical simulation The DEM model was validated using the field data measured on a single leg of a Jacket, JZ20-2 MUQ B2. This leg has an installed cone to break the level ice in flexural as shown in Fig. 2. The real geometries of the jacket cone were used in the simulation. Ice thickness, velocity and other environmental data were measured in field. In this model, finite level ice is used to model the infinite ice field, and the in-plane resistance on the modeled ice sheet edges from far-field level ice was represented using springs, as shown in Fig. 13(a). The stiffness of these springs is consistent with that between ice particles in DEM model. The geometry size of one leg of the Jacket is shown in Fig. 13(b), and the slope angle of icebreaking cone is 65 . Generally, a combination of uniaxial compressive and flexural strength of sea ice impacting the Jacket platform is needed to measure accurately to carry out the calibration of model parameters. Unfortunately, these two values were not measured. So we cannot determinate the contact stiffness and bonding strength in the DEM model directly. A trial computation was carried out until the ice loads time histories simulated with DEM are closer to the measurement results in field tests. Consequently, the model parameters used in the trial simulation are chosen and listed in Table 3. The ice cover is modeled with three layers of bonded particles with regular packing pattern, and subjected to drag force of the current. The drag force on the ice from the current was considered using the Morison theories (Sun and Shen, 2012). The particles at the far end and side boundaries move at a constant velocity in the ice advancing direction. Springs are set up in the horizontal and vertical direction to model the resistance from the far field ice cover in the DEM simulation. During the processes of the interaction between ice cover and the cone, fracture of the ice cover is represented by the breakage of bonds between associated particles. Dynamic ice loads on the cone are obtained by accumulating the impact load on the cone surface. In Fig. 13(a), one small portion is

Fig. 13. Construction of level ice with bonded spherical particles and the geometry size of one leg of JZ20-2 Jacket in ice cover-conical structure interaction simulation. 288

S. Di et al.

Ocean Engineering 146 (2017) 282–297

Fig. 14. Interaction between ice cover and conical structure simulated with DEM, here the magnitude of the velocity of each particle is presented in different colors. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Fig. 15. Detail views of the ice breaking process in the DEM simulation, the color of the particles represent the velocities of the particles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

compressive strengths range from 0.15 MPa to 5.9 MPa, and the flexural strengths range from 0.15 MPa to 2.4 MPa. So it is hard to obtain the compressive strength and flexural strength for a given block of ice

simultaneously. That is, the ratio of the compressive strength to flexural strength is just in an approximate range, rather than a certain value. In the parameter sensitivity analysis, the choice of σ sb =σ nb , ks =kn and μb had 289

S. Di et al.

Ocean Engineering 146 (2017) 282–297

Fig. 16. Dynamic ice loads on conical structure simulated with DEM and measured in the field.

5.2. Comparison of DEM results and ISO ice force standard The ISO 19906 (2010) provides empirical formulas for assessing the ice loads on conical structures based on physical model tests. There are two formulas based on plastic theory and elastic beam bending theory (refer to Appendix). The comparison of DEM simulation results and these two formulas will be presented here. (a) Plastic theory The computational parameters of ice load based on plastic theory in ISO standard are listed in Table 4.

Table 4 The parameters in ISO standard.

Fig. 17. Comparison of power spectral density of ice loads from field measurement and DEM simulation.

some slight deviation. This reason may lead to larger size of broken ice block than that in field measurement, and the time interval during two adjacent peak forces is longer than the field test result. Meanwhile, Due to the calculation capacity, three layers of particles were adopted to conduct the ice sheet. This may also be the reason inducing to the difference between two frequencies. The simulation will be fine if adopting more layers to conduct the ice sheet. Even so, we think that the model parameters stated above are reasonable, and can be used to calculate the ice loads on cones with different angles and diameters.

290

Definitions

Symbols

Values

Density of sea ice Diameter of cone at waterline Top diameter of the cone Ice thickness Flexural strength of sea ice Slope of cone Ice-structure friction coefficient Ice-ice friction coefficient Rubble height Porosity of ice rubble Angle of rubble Cohesion of rubble Friction angle of rubble

ρi w wT h σf α μ μi hr e θ c ϕ

920 kg=m3 2.85 m 1.708 m 0.23 m 0.54 MPa 65 0.02 0.15 0.3 m 0.3 10 0 kPa 15

S. Di et al.

Ocean Engineering 146 (2017) 282–297

In the table above, the flexural strength σ f is a key parameter in the calculation of ice load in the ISO standard. The determination of σ f is based on the relationship between the flexural strength of sea ice and the ratio of the particle size to specimen size, Eq. (24) can be rewritten as follows

σf ¼ 

3:47 þ 1:94 L=D

Table 5 Values of σ f and σ b in ice properties and ice-cone interaction simulations.

Parameter sensitivity analysis Ice-cone interaction simulation

(27)

We adopted three layers of elements in the simulation of ice load of ice cover on the conical structure, so L=D ¼ 3:0. The flexural strength σ f ¼ 0:78 MPa can be calculated using Eq. (27). The bonding strength between two particles corresponding to the value of σ f is σ b ¼ 0:73 MPa (Section 4.1). According to the linear relationship between σ f and σ b given by Reference (Ji et al., 2016) and the bonding strength of σ b ¼ 0:5 MPa used in the simulation (Section 5.1), the flexural strength of the numerical sea ice σ f ¼ 0:780:5 ¼ 0:54 MPa can be obtained finally. The 0:73 derivation of relationship between σ f and σ b is represented in Table 5. Substitute all computational parameters into Eqs. (A.1)–(A.12), we can get

Bonding strength between two spherical particle σ b

Flexural strength of ice sample under three-point bending σ f

0.73 MPa (Table 1)

0.78 MPa (Eq. (27))

0.5 MPa (Table 3)

0.54 MPa (Table 4)

against a conical structure with different slopes. In the simulation results, the breakage characteristics and the broken ice sizes of sea ice are shown clearly. Fig. 19 shows the maximum horizontal and vertical components of ice loads under various cone slopes, and their comparisons with the ice loads formulas in the ISO standard. The magnitude of ice load in the horizontal direction increases with the increase of cone slope, which has no obvious effect on the magnitude of ice load in the vertical direction. Compared with the ice loads of the plastic theory, the horizontal component load in the DEM simulation are much closer to those in the elastic theory, while the ice loads in the vertical direction in the DEM simulation agree well with those in the plastic theory. The locations of ice action on the conical structure are also important factors influencing the breaking mode and ice loads. Here, the effect of cone diameter is examined by changing the cone diameters D at waterline while keeping the cone angle as a constant of 60 . The cone diameters D are 1.81 m, 2.28 m, 2.74 m, 3.2 m and 3.66 m, respectively. The model parameters in the DEM simulation are listed in Table 3. Fig. 20 shows the failure process of level ice advancing against conical leg with diameters of 2.28 m and 3.2 m. In the simulation results, the breakage characteristics and the fragment sizes of sea ice are shown clearly. Fig. 21 shows the maximum ice loads in the horizontal and vertical components with various cone diameters at waterline, and their comparisons with the ice loads formulas in the ISO standard. The magnitude of ice loads in the horizontal and vertical direction increase with the increase of cone diameter. Compared to the ice loads of the plastic theory, the ice loads in the horizontal direction in the DEM simulation are much closer to those in the elastic theory, while the ice loads in the vertical direction in the DEM simulation agree well with those in the plastic theory.

FH ¼ 128:3 kN FV ¼ 74:6 kN

(b) Elastic theory The broken ice blocks accumulate around an ice breaking conical structure to form ice rubble. Besides failure induced load, rubble transportation also needs to be taken into count. There was no measurement for the geometrical and mechanical properties of ice rubble, including the height, porosity, angle, cohesion and the friction angle of ice rubble. Therefore, in the process of calibration with the elastic theory, some approximate ice rubble parameters were selected to calculate the ice load. The computational parameters of ice load based on elastic theory in ISO standard are also listed in Table 4. Substitute all computational parameters into the equations from Eq. (A.13)–(A.21), we can get

FH ¼ 107:6 kN FV ¼ 47:6 kN

5.4. Shielding effect of ice loads on multi-leg structure

From the comparison of the DEM simulations and the ISO standard of ice loads, we can find that the DEM model of ice load presented in this paper can give a comparatively accurate prediction in the simulation of the interaction between level ice and conical structure. In the simulation result, the horizontal component of ice load (109 kN) has a well agreement with the elastic theory (107.6 kN), while vertical component (70 kN) is close to the plastic theory (74.6 kN). The relationship between vertical ice load and horizontal ice load in Eq. (A.14) was obtained from a 2D theory of ice action on a sloping structure, which may be different from a 3D condition of ice action on a narrow conical structure.

To study the influence of the ice advancing direction on the ice actions on the multi-leg structure, a DEM numerical model of interaction between ice cover and JZ20-2 MUQ platform was presented. Fig. 22(a) shows a view of the full scale JZ20-2 platform, and (b) shows the DEM model of the platform constructed with cylindrical and conical elements. The cones are numbered in counter-clockwise order in which cone 1 and cone 4 are in the head-on direction of ice movement, and cone 2 and cone 3 are in the downstream direction. The span of two legs is 20 m. The size of ice cover was set as 50 m  40 m to reduce the boundary effect. The change of advancing angles is done by rotating the structure. The ice sheet maintain a consistent direction in the simulation, and the orientation of the structure θ was changed with respect to the direction of ice movement and was 0 , 5 , 10 , 15 , 20 , 25 , 30 , 35 , 40 , 45 . An interaction of 100-s duration between the Jacket platform and the ice cover of 904 860 particles was simulated in the calculation examples, requiring 88 h running time on the K40 GPU Card workstation. Fig. 23 shows the interaction process between ice cover and four-legs conical structure in the ice attacking direction of 0 , 15 , 30 , 45 . The channel forms at the rear of the two frontal cones, and influences the failure mode of ice sheet and ice load on the rear cones. For each numerical test, the ice load in the ice moving direction was obtained by

5.3. Influence of the cone slope and diameter During the interaction between ice cover and conical structures, the cone slope is also an important parameter, which influences the failure mode of ice sheet and the magnitude of the ice loads (Timco and Cornett, 1997; Xu et al., 2009). Here, the effect of cone slope is examined by changing the cone angles θ while keeping the cone diameter at waterline as a constant of 3.2 m. The cone angles θ are 40 , 45 , 50 , 55 , 60 and 65 , respectively. The model parameters in the DEM simulation are listed in Table 3. Fig. 18 shows the simulated failure process of level ice advances 291

S. Di et al.

Ocean Engineering 146 (2017) 282–297

Fig. 18. Interaction between sea ice and conical structure with different cone angles; particle color represents the velocity of particle. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

case, the frontal cones experienced the impact of intact ice and the two rear cones nearly suffered little affect. The rear cones are located in the channel formed by frontal cones. The rear cones contacted with the broken ice, and the ice loads on rear cones were lower than those on frontal ones contacted with intact ice sheet. In 5 ice angle case, cone 2 and cone 3 began to contact with the intact ice sheet, while most portion of these two cones still located in the channel formed by frontal cones. In this case, cone 2 contacted with the ice sheet which has two free boundaries, and cone 3 contacted with the ice sheet which have one free boundary. So the mean peak ice force on cone 2 was lower than cone 3. As the rotation angle of the jacket increases, cone 3 escaped from the channel formed by cone 4. It is obvious from Fig. 23 that the mean peak

recording and summarizing the contact force between the spherical particles and the cone. Fig. 24 shows the time-series behavior of the ice loads on the conical  structure in the ice moving direction for the numerical test of θ ¼ 15 . For this test, the ice sheet was approaching the cones at a constant speed. The zero level at the start of the time-series represents the time before the ice came in contact with the cone. The first cone to contact with the ice sheet is cone 4, the second one is cone 1, the third one is cone 3 and the final one is cone 2. The sequence of contact can also be obtained from timeseries of the ice loads on four cones. The peak ice forces picked from ice loads on all cones in the ice moving direction were averaged and plotted in Fig. 25. In 0 ice angle 292

S. Di et al.

Ocean Engineering 146 (2017) 282–297

Fig. 19. Comparison of the maximum ice loads on conical structures with different cone slopes in ISO standard and DEM simulation.

Fig. 20. Interaction between sea ice and conical structure with different cone diameters; particle color represents the velocity of particle. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Fig. 21. Comparison of the maximum ice loads on conical structures with different cone diameters in ISO standard and DEM simulation.

cone 2 was fully located in the channel formed by cone 4, so the mean peak ice load is much lower than other three cones.

force on cone 2 was lower than that on cone 3. At the process of the rotation of jacket from angle 0 –45 , the mean peak ice load on cone 3 increases and that on cone 2 increases at first then decreases with the increase of the rotation angle. In 45 ice attacking direction case, three cones (1, 2 and 3) directly faced the incoming ice sheet. Among these three cones, the mid-cone (cone 4) protrudes to be the forefront towards the ice impacting direction, and the ice failure process is not influenced by other cones, and the two back-side cones are also not influenced by cone 1. So their mean peak loads are very close to each other. In this case,

6. Conclusions In this paper, the simulations of the level ice advancing against a conical structure were conducted using the DEM model. The ice failure mode against the conical structure was discussed, and the ice loads were analyzed. Through qualitative comparison with full-scale measurements 293

S. Di et al.

Ocean Engineering 146 (2017) 282–297

Fig. 22. Four-legs conical platform JZ20-2 MUQ in the Bohai Sea and its DEM model.

Fig. 23. The interaction process between ice cover and four-legs conical structure in different ice directions; particle color represents the velocity of particle. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

technique, there is little possibility to conduct simulations either in realtime or with huge amount of elements. It is necessary to reduce the number of elements to equivalently simulate the mechanical properties of sea ice, so the size effect caused by particle sizes has to be considered. Some liner relationships between macro-mechanical properties of sea ice and relative element size (D=L) are obtained from the DEM numerical result. We can determine the bonding strength between two bonded particles using these linear relationships, and then, the same value of the

and industry experiences (ISO), it can be concluded that the simulation results using the DEM model provide reasonable information for the ice loads prediction of the offshore structure in ice regions. The determination of micro-scale parameters in the DEM model is an important part in the simulation. Firstly, the normal and shear contact stiffness, normal and shear bonding strength need to be determined according to the field test results of sea ice and the process of parameter determination presented in this paper. Base on the current hardware 294

S. Di et al.

Ocean Engineering 146 (2017) 282–297



Fig. 24. Dynamic ice loads on each cone in the ice attacking direction.θ ¼ 15 .

The numerical simulation results not only are validated with the field measurement but also provide more information about the continuous failure process and the effect of the cone slope and the acting position on the cone on the ice loads. It is found that the larger the value of the cone angle, the higher the ice load on cone. Apparently, a flat slope icebreaking cone should be equipped on the offshore structure. However, it is insufficient to only consider the cone angle; the influence of cone angle on the vibration of offshore structure should be taken into account. In general, the numerical simulation results can supplement the field data to some extent by providing more information about the offshore structure. However, limitations do exist in the numerical simulations. The ice-wind and ice-water interaction are not included in the numerical model, and the fracture-damage constitutive relation of sea ice is insufficient to model all the behaviors of sea ice. Regarding these issues, further studies and field- or model-scale measurements are required.

Fig. 25. The peak force averaged of each cone under various drifting directions.

bonding strength is used in the simulation of the interaction process between ice cover and offshore structure. The innovation of this paper is having established a numerical approach to predict the ice loads on conical structure using DEM method. For a level ice-covered region, as long as the mechanical properties of sea ice are known, the interaction between sea ice and structures can be simulated and the ice loads can be calculated.

Acknowledgements This work was supported financially by the National Natural Science Foundation of China (NSFC) [Grant Nos. 41606213, 51639004, 51579054], and the Fundamental Research Funds for the Central Universities [HEUCFP201701, HEUCFP201777].

Appendix The ISO 19906 (2010) provides empirical formulas for assessing the ice loads on conical structures based on physical model tests. There are two formulas based on plastic theory and elastic beam bending theory. (a) Plastic theory This method is based on a limit analysis solution for level ice actions on upward and downward breaking cones. The model considers ice loads due to the flexural failure of the ice sheet and the ride-up ice loads due to ice pieces. The total ice loads in the horizontal and vertical directions as follows, and some parameters can be expressed as

295

S. Di et al.

Ocean Engineering 146 (2017) 282–297

f ¼ sinα þ μE1 cosα

(A.1)

f cosα  μE2 2 sin α þ μαcosα 4

(A.2)

α sinα þ cosα 2 sin α þ 2μαcosα 2

(A.3)

hV ¼ π

gr ¼ π

W ¼ ρi ghr

w2  w2T 4cosα

(A.4)

where α is the slope of the concial structure measured from the horizontal; hr is the ice ride-up thickness (hr  h); μ is the ice-structure friction coefficient; w is the waterline diameter of the cone; wT is the top diameter of the cone. The effects of a rubble accumulation on the cone can be considered by using a value that exceeds the single sheet thickness for the ride-up thickness. The parameters E1 and E2 are the complete elliptical integrals of the first and second kind, defined as given as follows π=2 

E1 ¼ ∫ 0

π=2 

E2 ¼ ∫ 0

1  sin2 αsin2 η 1  sin2 αsin2 η

1=2 1=2

(A.5)

dη

(A.6)

dη

Assuming a single sheet thickness of ride-up ice, the horizontal ride-up HR ice load and the vertical ride-up ice load VR are obtained as follows

HR ¼ W

tanα þ μE2  μfgr cosα 1  μgr

VR ¼ Wcosα

(A.7)

π cosα  μα  fhV þ HR hV a 2

(A.8)

The horizontal breaking ice load HB and the vertical breaking ice load VB are given as follows

HB ¼

 σ f h2 tanα 1 þ Yxlnx þ Gðx  1Þðx þ 2Þ x1 3 1  μgr

(A.9)

VB ¼ HB hV

(A.10)

 1=2 where Y ¼ 2:711, G ¼ ðρi gw2 Þ=ð4σ f hÞ; x ¼ 1 þ 3G þ Y2 , ρi is the density of sea ice and σ f is the flexural strength of sea ice. The total ice load components in the horizontal and vertical directions are given by following equations

FH ¼ HB þ HR

(A.11)

FV ¼ VB þ VR

(A.12)

(b) Elastic theory In the ISO standard, a formulas of ice load based on elastic beam theory is also given. When an ice sheet acts on a cone, the flexural failure component can be evaluated considering the ice sheet as an elastic beam on elastic foundation.

FH ¼

HB þ HP þ HR þ HL þ HT 1  σHf lBc h

(A.13)

FV ¼

FH ξ

(A.14)

ξ¼

sinα þ μcosα cosα  μsinα

(A.15)

where HB is the breaking ice load; HP is the load component required to push the sheet ice through the ice rubble; HR is the load to push the ice blocks up the slope through the ice rubble; HL is the load required to lift the ice rubble on top of the advancing ice sheet prior to breaking it; HT is the load to turn the ice blok at the top of the slope.

296

S. Di et al.

Ocean Engineering 146 (2017) 282–297

The ice load component HB is given as follows

HB ¼ 0:68ξσ f

0:25    ρw gh5 π 2 Lc wþ E 4

(A.16) 

where E is the elastic modulus; ν is the Poisson ratio; ρw is the density of sea water and Lc ¼ The load component HP is expressed as follows

Eh3 12ρw gð1ν2 Þ

1=4 .

tanμ 2 1 HP ¼ wh2r μi ρi gð1  eÞ 1  tanα 2tanθ

(A.17)

where hr is the rubble height; μi is the ice-to-ice friction coefficient; e is the porosity of the ice rubble; θ is the angle the rubble makes with the horizontal. The load component HR is given as follows

HR ¼

wP cosα  μsinα

(A.18)

where

     1 1 tanθ cosα tanθ sinα þ μcosα  1 þ 0:5ðμi þ μÞρi gð1  eÞh2r 1 þ hr hρi g P ¼ 0:5μi ðμi þ μÞρi gð1  eÞh2r sinα tanθ tanα tanα tanα tanα sinα The load component HL is given as follows

  2     1 1 tanθ tanθ tanθ  1 þ 0:5wh2r ρi gð1  eÞξtanϕ 1  HL ¼ 0:5wh2r ρi gð1  eÞξ þ ξcwhr 1  tanθ tanα tanα tanα tanα

(A.19)

where c is the cohesion of the ice rubble; ϕ is the friction angle of the ice rubble. The load component HT is given by as follows

HT ¼ 1:5wh2 ρi g

cosα sinα  μcosα

(A.20)

Equation (A.13) need to be modified to account for compressive stress in ice sheet as a result of the applied horizontal load. The flexural strength is modified as follows ð1Þ

σf ¼

FH þ σf lc h

(A.21)

where lc is the total length of the circumferential crack, and estimated as lc ¼ w þ π4 Lc . 2

Poloj€arvi, A., Tuhkuri, J., 2009. 3D discrete numerical modelling of ridge keel punch through tests. Cold Reg. Sci. Technol. 56, 18–29. Potyondy, D.O., Cundall, P.A., 2004. A bonded-particle model for rock. Int. J. Rock. Mech. Min. 41, 1329–1364. Qu, Y., Yue, Q.J., Bi, X.J., K€arn€a, T., 2006. A random ice force model for narrow conical structures. Cold Reg. Sci. Technol. 45, 148–157. Sun, S.S., Shen, H.H., 2012. Simulation of pancake ice load on a circular cylinder in a wave and current field. Cold Reg. Sci. Technol. 78, 31–39. Timco, G.W., Cornett, A.M., 1997. The influence of variable-thickness ice on the loads exerted on sloping structures. Cold Reg. Sci. Technol. 26, 39–53. Wang, Y., Alonso-Marroquin, F., 2009. A finite deformation method for discrete modeling: particle rotation and parameter calibration. Granul. Matter 11, 331–343. Withalm, M., Hoffmann, N.P., 2010. Simulation of full-scale ice-structure-interaction by an extended Matlock-model. Cold Reg. Sci. Technol. 60, 130–136. Xu, N., Qu, Y., Yue, Q.J., Bi, X.J., Palmer, A.C., 2009. Main factors of ice sheet–conical structure interaction process based on field monitoring. In: Proceedings of the ASME 2009 28th International Conference on Ocean, Offshore and Arctic Engineering (OMAE2009), Honolulu, Hawaii, USA. May 31 – June 5, 2009, OMAE2009–79848. Xu, N., Yue, Q.J., Qu, Y., Bi, X.J., 2011. Results of field monitoring on ice actions on conical structures. J. Ocean. Mech. Eng. 144, 041502. Yang, B., Jiao, Y., Lei, S., 2006. A study on the effects of microparameters on macroproperties for specimens created by bonded particles. Eng. Comput. 23, 607–631. Yue, Q.J., Bi, X.J., 2000. Ice induced jacket structure vibration. J. Cold Reg. Eng. 14 (2), 81–92. Yue, Q.J., Qu, Y., Bi, X.J., K€arn€a, T., 2007. Ice force spectrum on narrow conical structures. Cold Reg. Sci. Technol. 49, 161–169.

References Brown, T.G., Tibbo, J.S., Tripathi, D., Obert, K., Shrestha, N., 2010. Extreme ice load events on the Confederation Bridge. Cold Reg. Sci. Technol. 60, 1–14. Brown, T.G., M€ a€ att€ anen, M., 2009. Comparison of Kemi-I and Confederation Bridge cone ice load measurement results. Cold Reg. Sci. Technol. 55, 3–13. Dalane, O., 2014. Influence of pitch motion on level ice actions. Cold Reg. Sci. Technol. 108, 18–27. Eranti, E., 1991. General theory of dynamic ice structure interaction with applications. In: Proceeding of the First International Offshore and Polar Engineering Conference (ISOPE), Edinburgh, UK, pp. 489–498. Hopkins, M.A., 1997. Onshore ice pile-up: a comparison between experiments and simulations. Cold Reg. Sci. Technol. 26, 205–214. Huang, H., 1999. Discrete Element Modeling of Tool-rock Interaction (PhD. thesis). University of Minnesota. Huang, Y., 2010. Model test study of the non-simultaneous failure of ice before wide conical structures. Cold Reg. Sci. Technol. 63, 87–96. Huang, Y., Ma, J.J., Tian, Y.F., 2013. Model tests of four-legged jacket platforms in ice: Part 1. Model tests and results. Cold Reg. Sci. Technol. 95, 74–85. ISO, 2010. Petroleum, and Natural Gas Industries –Arctic Offshore Structures. ISO 19906. Ji, S.Y., Di, S.C., Long, X., 2016. DEM simulation of uniaxial compressive and flexural strength of sea ice: parametric study. J. Eng. Mech.-ASCE C4016010–C4016011. M€ a€ att€ anen, M., Hoikkanen, A.N., Avis, J., 1996. Ice failure and ice loads on a conical structure Kemi-1 cone full scale ice force measurement data analysis. In: Proc 13th International Symposium on Ice (IAHR), Beijing, China, pp. 8–17. Matlock, H., Dawkins, W.P., Panak, J.J., 1971. Analytical model for ice-structure interaction. J. Eng. Mech. Div. 97, 1083–1092. Paavilainen, J., Tuhkuri, J., Poloj€arvi, A., 2011. 2D numerical simulations of ice rubble formation process against an inclined structure. Cold Reg. Sci. Technol. 68, 20–34.

297