Efficient discrete element simulation of managed ice actions on moored floating platforms

Ocean Engineering 190 (2019) 106483

Contents lists available at ScienceDirect

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

Efficient discrete element simulation of managed ice actions on moored floating platforms Mochamad Raditya Pradana a, Xudong Qian a, *, Aziz Ahmed b a b

Department of Civil and Environmental Engineering, National University of Singapore, 1 Engineering Drive 2, Singapore, 117576, Singapore School of Civil Mining and Environmental Engineering, University of Wollongong, Australia

A R T I C L E I N F O

A B S T R A C T

Keywords: Discrete element method Ice-structure interaction Managed ice Floating platforms Arctic

The current design standard lacks the detailed approach to estimate managed ice actions, which in turn forces the industry to rely on ice basin model tests and calibrated numerical simulations. To address the need for a reliable assessment of managed ice actions on large-scale floaters, this paper proposes an efficient two-dimensional discrete element simulation with polygonal elements, which integrates smoothly with a mooring system modeling and a hydrodynamic force representation. The study aims at establishing a range of numerical pa­ rameters suitable for managed ice, as well as the boundaries within which the two-dimensional assumption remains acceptable. The study includes validations against the empirical data on isolated ice impacts and the towing model test conducted in the Hamburg Ship Model Basin (HSVA), for which the discrete element simu­ lations demonstrate reasonable agreements. Using this validated tool, this paper conducts a case study using a hypothetical floater based on the Kulluk, a heavily studied conical drilling platform with years of deployment in the Arctic. The study examines the effect of different ice and mooring characteristics on the simulated load, which provides insights on the design of floating platforms and its mooring systems for managed ice conditions.

1. Introduction According to the U.S. Geological Survey (Bird et al., 2008), the Arctic region potentially holds about 22% of the world’s undiscovered oil and gas, of which more than 80% is expected to be offshore. Hamilton (2011) projected that approximately half of the untouched oil reserve exists in water depths greater than 100 m, where bottom-founded structures are no longer cost-effective for the Arctic environments (Paulin et al., 2008). Floating platforms with reliable station-keeping systems become the only realistic option for drilling and other key offshore operations in the Arctic’s deep-water. As unmanaged sea ice imposes load beyond the capacity of conventional mooring systems, which leads to failures at a low ambient temperature (Zhang and Qian, 2017; Feng and Qian, 2018), floating platforms are usually deployed with the support from ice management (Hamilton et al., 2011). Fig. 1 illustrates an ice management technique where two icebreakers move in a circular pattern to reduce the size of incoming ice floes to a safe level. Although the ice loads are reduced significantly, the floating platform and its mooring system still have to endure the managed ice actions while maintaining an acceptable offset for drilling or other critical operations.

The assessment of loads imposed by managed ice remains a chal­ lenging issue due to its complex and uncertain nature. For floating platforms in managed ice, the current design standard (ISO19906, 2010) recommends an empirical approach based on the field data collected on the Kulluk (Wright, 2000). The approach relies on data correction (proposed by Keinonen et al., 1996) to reflect the influence of different vessel dimensions, hull forms, hull surface conditions, ice strengths, and ambient temperatures. In the more rigorous attempt to assess the managed ice load, numerous ice basin model tests have investigated a wide range of cases, including managed ice load on ships (Zhou et al., 2019; Wang et al., 2016; Jenssen et al., 2012; Karulin and Karulina, 2011; Løset et al., 1998), semi-submersibles (Karulin et al., 2012; Kar­ ulin and Karulina, 2014; Williams, 1989), conical platforms (Cole, 2005; Abdelnour et al., 1987), and iceberg towing operations (Eik and Mar­ chenko, 2010). Comfort et al. (1999) presented a discussion on the trend observed in the previous model test data and concluded on the relative importance of the ice concentration, ice thickness, ice velocity and the correlation between the model tests. Due to the resource-intensive na­ ture of ice basin tests, various researches have invested in well-calibrated numerical models to extend the investigation beyond the scope of the basin tests. From the large scale perspective, ice floes in nature bear resemblance

* Corresponding author. E-mail address: [email protected] (X. Qian). https://doi.org/10.1016/j.oceaneng.2019.106483 Received 14 March 2019; Received in revised form 29 August 2019; Accepted 20 September 2019 Available online 25 September 2019 0029-8018/© 2019 Elsevier Ltd. All rights reserved.

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

Nomenclature Aoverlap

e eix ; eiy

Planar area of the overlapping region

h l m* Δt

Dfloe Average floe diameter E Elastic modulus of the discrete elements/particles Fdrag Water drag force Fmoor; damp Mooring damping force Fmoor Mooring restoring force Normal component of the contact force Fn Ft Tangential component of the contact force MCD Mean caliper diameter Horizontal tension of the ith mooring line Ti c Initial ice concentration

α ζn ζt

μ ψ

to a system of discrete particles. This inspired many scholars to inves­ tigate the behavior of ice floes using the Discrete Element Method (DEM), where the individual ice floe is represented using a single,

coefficient of restitution Eccentricity of the ith mooring line tension in x- and ydirections Thickness of the discrete elements/particles Length scale of the two contacting particles Effective mass of the two contacting particles Numerical time step Scaling factor for the Voronoi polygons Normal damping ratio Tangential damping ratio Coefficient of friction Normalized domain width

€rvi, 2018). The earlier at­ unbreakable particle (Tuhkuri and Poloja tempts utilized 2D disk elements to simulate the ice floe fields in the horizontal plane (Babic, 1990; Hopkins and Hibler, 1991; Løset, 1994a,

Ice Drift

Fig. 1. Ice management operation, adapted from Hamilton et al. (2011). 2

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

1994b; Hansen and Løset, 1999a, 1999b). One notable study by Løset (1994a, 1994b) demonstrated that even with simplified floe geometry and restriction to the horizontal plane, the load predicted by DEM aligned well with the model test results. The 2D DEM developments later incorporated polygonal elements to mimic the planar geometry of real ice floes (Zhan et al., 2010; Rabatel et al., 2015; Yulmetov et al., 2016; Yulmetov and Løset, 2017). Despite the promising results, 2D simulation inherently restricts the possibilities for out-of-plane motions, such as rafting and overturning of ice floes, which are inevitable in the case of high ice concentration with significant lateral confinement. DEM studies naturally progressed to 3D formulation utilizing disk elements (Hopkins and Tuhkuri, 1999; Hopkins and Tuthill, 2002; Ji et al., 2013; Sun and Shen, 2012), polyhedral elements (Lau et al., 2011; Lau, 2006; Lau and Simoes R� e, 2006; Zhan and Molyneux, 2012), dilated polyhedral ele­ ments (Hopkins, 2004; Liu and Ji, 2018) and spherical elements (Pra­ dana and Qian, 2019). Majority of the aforementioned 3D simulations demonstrated a close match to their benchmarking ice basin tests. Despite the advancements in 3D DEM, a number of recent studies still advocate the use of 2D DEM for managed ice problems, which evidently yielded good or at least qualitative agreement with the experimental results when the ice concentration and the lateral confinement is not significant (Yulmetov et al., 2016; Yulmetov and Løset, 2017; Karulin et al., 2012; Karulin and Karulina, 2011). In support of the previous statement, the 3D simulations in Hopkins and Tuhkuri (1999) revealed that a confined field of floating ice floes began showing 3D motions (rafting and overturning) only when the surface concentration reaches 79%. This suggests that there exists a limiting condition, possibly in terms of initial ice concentration or relative floe size, below which 2D simulations remain reasonably accurate. Past 2D investigations, how­ ever, have not clearly established these limiting conditions. Addressing the limit is therefore the key to maximize the potential of 2D solutions as compliments to the design standard and the ice basin tests. This study aims to address the assessment of loads generated by managed ice, specifically on moored floating platforms. In line with this aim, this study proposes a 2D discrete element framework where the ice floes are represented using random convex polygons created using the Voronoi Tessellation procedure. Via comparison against ice basin tests, this study attempts to establish a validity range where the 2D assump­ tion remains reasonably accurate. Despite the inherent limitations, the current 2D framework consumes considerably less computational re­ sources and execution time as compared to a 3D simulation. Thus, once the validity range is identified, the 2D framework holds significant po­ tential as a preliminary analysis tool for engineering applications (Ghoshal et al., 2018). The simplicity of the 2D formulation also allows for modeling various dynamic systems with relative ease. Hence, the framework developed herein is easily extendable to other relevant topics, such as vessels with dynamic positioning under managed ice (Kjerstad et al., 2015). The numerical framework in this paper is an extension of that in Matuttis and Chen (2014), which is originally designed to simulate granular materials. In Section 2, this paper describes the underlying formulation in the contact forces, the Voronoi-based ice floe creation, the moored floating platform modeling, and the simplified hydrody­ namic forces. Next, Section 3 validates the current approach against the field data on isolated ice floe impacts, as well as the iceberg towing model test presented in Eik and Marchenko (2010). As a demonstration of the proposed approach, Section 5 presents a case study on a hypo­ thetical floating platform based on the Kulluk (Wright, 2000), where the effect of the ice and the mooring system characteristics is examined. Finally, Section 6 summarizes the conclusions supported by the findings of this study.

2. Discrete element formulation 2.1. Overview The current framework computes the interaction between each discrete element according to the soft particle concept (Cundall and Strack, 1979), which establishes the contact force as a function of the overlap between non-deformable elements/particles. The particles may also experience other external actions, such as hydrodynamic or wind forces. By applying time-integration to the equation of motion, the framework calculates the position and the velocity of each particle at the next time-step. Due to the 2D assumption, the framework computes strictly the water-plane motion (two translational and one rotational degrees of freedom), and assumes that all out-of-plane actions are in equilibrium. In the current approach, the typical simulation comprises the simu­ lation walls, the ice floes, the floating platform, and the representation of the mooring system. The first three components are discrete elements with specific geometry and assignments. The simulation walls are un­ movable rectangular particles placed at the extremes of the domain. The ice floes are arbitrarily-shaped polygonal particles that move in response to the external forces. On the other hand, the floating platform is a particle that possesses the platform’s water plane geometry, mass, and mass moment of inertia. Motions of the platform elicit restoring forces from the mooring system, which is represented using a set of springs. Successive time-integration over a pre-determined duration provides the time-series of: 1) ice actions on the platform, 2) mooring restoring forces, and 3) motions of the platform. 2.2. Contact formulation Fig. 2 illustrates a contact between two polygonal particles. As illustrated, the particles in contact experience normal and tangential forces that are generated by the spring-dashpot system in the corre­ sponding direction. These forces act on the centroid of the overlapping region (see Fig .2a), and, due to the eccentricity with respect to the particle centroids, they produce in-plane moments as well. The geom­ etry of the overlap region (Fig. 2c) defines the normal and the tangential direction. The vector connecting the two intersection points represents the tangential direction, while the normal direction follows accordingly (see Fig. 2c). Details regarding the contact detection and the computa­ tion of the overlap geometry are available in Matuttis and Chen (2014). Following the traditional formulation for the soft particle contact, the current contact forces include, Fn ¼ Fns

Fnd � 0; ðNormal componentÞ

(1a)

Ft ¼ Fts

Ftd � μFn ; ðTangential componentÞ

(1b)

where the subscripts s and d designate the elastic and the damping components, respectively. The formulation restricts the negative normal force (signifying unphysical attractive force), which may occur when the damping term is dominant. On the other hand, the tangential force imitates friction, and as such, its magnitude is limited by the coefficient of friction (μ). At the end of every contact computation, the framework resolves the forces and moment according to the global coordinate system, resulting in Fcontact;x ; Fcontact;y and Mcontact at the centroid of each particle. 2.2.1. Normal forces The normal elastic force depends on the extent of overlap and the assigned material stiffness. At any given time-step, its magnitude follows, Fns ¼ Eh

3

Aoverlap ; l

(2)

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

(a)

Ft r2

r1

Ft M2 Particle 2

M1 Particle 1

Fts

Ftd

Fn Fnd

Ft Intersection Points

Fns

Particle 2

Fn

(b)

Particle 1

Aoverlap

Fn

Ft <

(c)

n

lcontact

Fig. 2. The particle-particle contact model.

of friction. To avoid non-decaying tangential oscillations, the current contact model includes a tangential damping force similar to the viscous damper in the normal direction, pffiffiffiffiffiffiffiffiffiffiffiffiffi (6) Ftd ¼ 2ζt ηEhm* ðvt Þ:

where, l¼

4r1 r2 : r1 þ r2

(3)

The parameters r1 ; r2 represent the distance between the corre­ sponding particle centroid and the center of the overlap area (see Fig. 2). The constants E and h denote the particle elastic modulus and the out-ofplane thickness, respectively. Together, they characterize the normal contact stiffness. Meanwhile, the area of the overlap region (Aoverlap ) and the particle length scale (l) quantify the extent of overlap. The normal damping force is modeled after the viscous damper in harmonic oscillators, and expressed as, ! 1 pffiffiffiffiffiffiffiffiffiffiffi Aioverlap Aioverlap 1 * Fnd ¼ 2ζn Ehm (4) ⋅ ; l Δt

2.2.3. Upper limit for the normal contact force The normal contact force in Eq. (1a) may exceed the strength of the simulated material and produces unrealistically high loads. This study proposes a force limit that reflects the compressive strength of the ma­ terial. The limit ensures that the maximum pressure on the contact area (i.e., Ac ¼ h � lcontact , see Fig. 2c) does not exceed the material’s compressive strength. Based on the ice pressure data from ship ramming tests, ISO19906 (2010) proposed the pressure-area curve that follows,

σ ¼ CP ADc P ;

where ζn represents the normal damping ratio, and m* denotes the effective mass of the two particles in contact. The term inside the pa­ rentheses is an approximation of the relative velocity based on the change of the overlap area over a single time-step (Δt). During separa­ tion, i.e., when the relative velocity is positive, the damping force may lead to a negative (attractive) total force in Eq. (1a) if the damping force exceeds the elastic force, which may occur when ζn or the initial impact velocity is high. The limit Fn � 0 (see Eq. [1a]) prevents this unphysical behavior.

(7)

where the coefficients CP and DP equal 3.0 and 0.4, respectively. Consistent with the focus on ice, the current formulation limits the normal contact pressure (Fn =Ac ) to that in Eq. (7). This study therefore does not simulate the physical crushing failure in the ice floes for confined ice fields. A 3D DEM approach will be able to predict such failures in the ice floes with enhanced accuracy (Pradana and Qian, 2019). 2.3. Ice field creation with Voronoi Tessellation

2.2.2. Tangential forces The tangential elastic force adheres to the Cundall-Strack friction model (Cundall and Strack, 1979), which computes the force incre­ mentally as follows, � � F its ¼ Fits 1 þ ηEh vit Δt � μFn : (5)

Realistic assessment of managed ice loads requires the simulation of the ice field’s randomness along with its essential characteristics, namely the ice concentration and the average floe size. The ice con­ centration represents the fraction of the water plane area that is occu­ pied by ice floes, while the ice floe size is often quantified using the mean caliper diameter (MCD). For any given polygon, this study approximates the MCD as twice the average centroid-to-vertex distance. A number of publications have generated artificial ice fields using existing or self-developed algorithm. Yulmetov et al. (2016) propose an algorithm that generates, scales, and places arbitrary polygons in such a way that the target size distribution and the desired ice concentration are achieved. Other scholars utilize the 2D Voronoi Tessellation to

The multiplier η represents the ratio between the tangential and the normal spring stiffness. This paper assumes that η ¼ 0.38, based on the shear-to-normal modulus correlation with a Poisson ratio of 0.3. The bracketed term approximates the relative tangential motion based on the current relative tangential velocity (vit ). At any time-step, the tangential elastic force must not exceed the ceiling set by the coefficient 4

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

achieve similar result (Liu and Ji, 2018; Lau and Simoes Re�, 2006). Considering its efficiency, this study adopts the Voronoi Tessellation algorithm to generate the managed ice fields. In the 2D context, the Voronoi Tessellation divides a planar domain into n convex polygons (see Okabe et al. (2009) for the detailed algo­ rithm). The variable n equals the number of randomly placed Voronoi generators (also referred to as seed points) that are present within the domain. The edges of the polygons are equidistant lines between one Voronoi generator and its adjacent neighbors, such that every polygon encloses precisely one Voronoi generator. Without further modification, the Voronoi polygons form a 100% concentration. Nevertheless, any desired concentration is achievable by scaling down the polygons, while the average floe size (in terms of MCD) can be regulated by adjusting n according to the domain size. In a case where every polygon is scaled uniformly, the concentration (c) and the scaling factor (α) obey the pffiffi following relationship, α ¼ c. Fig. 3a shows an example ice field with c ¼ 90% and a mean floe size of Dfloe ¼ 28 m. The drawback of a Voronoi-Tessellation-based approach is the lack of a direct control over the floe shape and size distribution. Nonetheless, indirect controls over these distributions are possible through the stra­ tegic positioning of Voronoi generators (Sotomayor and Tippur, 2014). Another method to control the shape distribution involves the Cen­ troidal Voronoi Tessellation (CVT), where an iterative tessellation is performed until the polygon centroids and the Voronoi generators converge (Du et al., 1999). This method requires that for the ith iteration, the Voronoi generators are the centroid of the polygons from the

floes. On the other hand, if the CVT is stopped before convergence, the resulting polygons appear to be more natural and thus closer to the real ice floes (see Fig. 3b). For this reason, the current study adopts the partial CVT, where the algorithm is stopped after the second iteration. 2.4. Modelling the mooring system In the current approach, the individual mooring line is a spring that follows a user-defined horizontal tension-offset curve (Fig. 4). This method enables the modeling of different mooring types (taut or cate­ nary, of any size and pre-tension) by simply providing the appropriate tension-offset curve. Furthermore, the layout of the mooring system (spread or grouped, in any orientation) is easily modifiable by adjusting the coordinates of the fairlead and the anchor points. Fig. 4b exemplifies a typical tension-offset curve required by the current framework. As shown, the curve should cover both the positive and the negative offset that are expected to occur during the simulation. The example curve in Fig. 4b represents the response of a single mooring line to a horizontal offset, which, for better accuracy, can be generated using a specialized mooring analysis tool. Throughout the simulation, the framework as­ sumes that the mooring offset equals the change in the projected mooring length. The mooring system provides restoring forces and moment propor­ tional to the platform’s motion. The magnitudes of these restoring ac­ tions equal the resultant of all mooring line tension, and are expressed as,

ði 1Þth iteration (except for the 1st iteration, where the Voronoi gen­ erators are randomly placed). Fig. 3a–c illustrates the polygon shape evolution as the iteration is performed, while Fig. 3d shows a field of ice floes in nature. As shown, the CVT algorithm eliminates the elongated polygons (Fig. 3c), yet they produce unnaturally rounded and uniform

k X

Fmoor; x ¼

Ti cosθi ;

(8a)

i¼1

Fig. 3. The ice field generated using the Centroidal Voronoi Tessellation: (a) after 1 iteration; (b) after 2 iteration; (c) after 10 iteration; and (d) the real ice field for comparison purpose, adapted from https://discoveringthearctic.org.uk. 5

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

T1

T8

(a)

T2

e2,x

T7 y

e2,y

4000

Horizontal Tension (kN)

3000

2

(b)

Assigned to line 2, 3, 6, 7 Assigned to line 1, 4, 5, 8

2000 Pre-tension

x

1000

T6

T3 T5

0 -3

T4

-2

-1

0

Offset (m)

1

2

3

Fig. 4. The modeling of floating platform in the current 2D DEM: (a) the rigid body representing the hull and the springs simulating the mooring system; and (b) the horizontal tension-offset curve assigned to the individual springs. k X

Fmoor; y ¼

Similar to other approximations of the fluid drag, Eq. (10) is a function of the fluid density (ρw ) as well as the relative velocity between the fluid (Vw ) and the ice floe (v). The current framework computes the drag force in the x- and y-directions independently, and applies them to the centroid of each particle. As suggested in Lu et al. (2011), the value of drag coefficients are Cf ¼ 1.0 and Cs ¼ 0.002. The formulation in Eq. (10) lacks the fluid resistance against in-plane rotations, which can lead to unrealistic floe movements (e.g., unre­ strained rigid-body rotation). To alleviate this, the current study im­ poses a drag moment similar to that in Ji et al. (2013), � � Cr D 2 D Mdrag ¼ ρ jωjω; (11) ðh⋅DÞ 2 w 2 4

(8b)

Ti sinθi ; i¼1

k X

Mmoor ¼

Ti cosθi *ei;y þ

i¼1

k X

Ti sinθi *ei;x :

(8c)

i¼1

where Ti represents the tension of line- i according to the assigned tension-offset curve, θi denotes the angle between the mooring line and the horizontal x-axis (see Fig. 4a), and the parameters ei;x and ei;y define the eccentricity between the mooring force and the centroid of the floater in the x and y directions, respectively. In addition to these forces, the motion of the submerged mooring lines introduces a water drag that exerts a damping effect on the system. Currently, the energy dissipation contributed by the mooring system is simulated using viscous dampers, which generate the following forces and moment, pffiffiffiffiffiffiffiffiffiffiffi Fmoor;damp; x ¼ 2ζmoor Kx mflt ⋅vflt;x ; (9a) Fmoor;damp; y ¼ Mmoor;damp ¼

2ζmoor 2ζmoor

pffiffiffiffiffiffiffiffiffiffiffi Ky mflt ⋅vflt;y ;

pffiffiffiffiffiffiffiffiffiffiffiffi Krot Jflt ⋅ωflt :

where ω is the particle’s angular velocity. The rotational drag coefficient (Cr ) assumes a value of 1.0 following the form drag coefficient. In addition to the drag forces, the current framework considers the added mass effect that occurs when an ice floe accelerates or decelerates in water. Meylan et al. (2015) demonstrated experimentally that the added mass on an ice floe constitutes about 10% of the actual mass. Following this finding, the ice floe mass considered in the simulation equals,

(9b) (9c)

Kx ; Ky ; and Krot represent the translational and the rotational stiffness of the mooring system, respectively. These stiffness parameters derive from separate static offset analyses, performed before the ice-structure simulation. The constant ζmoor represents the mooring system’s damp­ ing ratio, while mflt and Jflt denote the mass and the mass moment of inertia of the floating platform. The total mooring action applied to the platform’s centroid equals the summation of Eqs. (8) and (9).

m ¼ m0 ð1 þ Cm Þ;

(12)

where m0 is the actual mass of the ice floe, while the added mass coef­ ficient (Cm ) equals 0.1 in all simulations. 2.6. Equation of motion and time integration The x and y translational motions of a discrete element adheres to the Newton’s equation of motion,

2.5. Hydrodynamic forces

m

Lu et al. (2011) identified that the water drag force acting on a relatively flat ice floe equals the sum of the form drag and the skin drag, and follows, ! " �1=2 #2 � Cf c h Fdrag ¼ h⋅D þ Cs A ρw jVw vjðVw vÞ; (10) 1 1 cD 2

� dv ¼ Fcontact þ Fdrag þ Fmoor þ Fmoor;damp ; dt

(13)

where v designates the particle’s velocity in the corresponding direction. The right-hand side of Eq. (13) comprises the time-varying external forces that may act on any given discrete element. Clearly, the mooring forces are zero when the particular discrete element represents an ice floe. Parallel to the translational motion, the centroidal rotation of a discrete element follows,

where Cf and Cs represent the form and the skin drag coefficient, respectively. The form drag refers to the normal force acting on the vertical face of an ice floe as it moves through the water. Accordingly, it depends on the floe thickness (h) and the floe diameter (D). However, the wake effects in a field of ice floes attenuate this force in proportion to the ice concentration (c), as implied by Eq. (10). The skin drag, on the other hand, describes the tangential force acting on the bottom face of an ice floe, and as such, is governed by the planar area of the ice floe (A).

J

� dω ¼ Mcontact þ Mdrag þ Mmoor þ Mmoor;damp ; dt

(14)

where ω represents the particle’s angular velocity while J denotes the particle’s mass moment of inertia. By integrating Eqs. 13 and 14 with respect to time, the DEM predicts the elements’ position and velocity throughout the simulation. The time6

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

ice

Ek = 12 mvx 2 Cf = 0

Cs = 0 Cr = 0

Cm = 0

3. Validation of the proposed approach

�

" ln e

π

1þ

� �2 # ln e

π

1 2

:

=0

x Fig. 5. Setup for the isolated ice floe impact simulations.

reported in the empirical database. The simulations assigned three different values for the elastic modulus, namely E ¼ 1, 10, and 100 MPa. For every elastic modulus considered, the simulation comprises 530 unique cases. Fig. 6 contrasts the impact force-energy data in Timco (2011) against those obtained from the simulations. As illustrated, both the measured and the simulated impact forces increase with the kinetic energy, yet there are significant amounts of scatter. For the measured data, the scatter originates from many different factors. Some of the prominent factors are arguably the mechanical properties of ice (including strength and modulus), since the ice floes are reported to indicate varying extent of failure after collision. For the DEM results, since the ice strength is uniformly based on Eq. (7), the scatter occurs mainly due to eccentric collisions, which transform parts of the kinetic energy into rotational motion (Matskevitch, 1997). Thus, depending on the exact floe geom­ etry, the same kinetic energy may translate into a different impact force. As demonstrated through Fig. 6a–c, the DEM impact force increases with the elastic modulus. For the elastic modulus of 1–100 MPa, the DEM results occupy the margin between the mean and the upper-bound curve established by Timco (2011). However, the DEM and the measurements achieved the best agreement when E ¼ 1 MPa (Fig. 6a), which is several orders lower than the elastic modulus of sea ice (2–6 GPa according to ISO19906 [2011]). Earlier ice floe simulation with 2D disks similarly used a low elastic modulus of 0.1 MPa, which was justified by the close comparison to the ice basin tests (Løset, 1994a, 1994b; Hansen and Løset, 1999a, 1999b). This appears to be the consequence of the soft-particle approach, where the particle overlap is merely a crude estimation of the actual deformation. The elastic modulus under this framework, therefore, may not adhere to the physical value and should be determined through calibration against the relevant experimental data. Based on the calibration herein, the subsequent simulations adopt the elastic modulus of E ¼ 1 MPa.

(16b)

With the appropriate the contact parameters, the current DEM is ex­ pected to reproduce the empirical trend. In DEM, the contact behavior highly depends on the material prop­ erties assigned to the individual element. These properties include the elastic modulus (E), the damping ratio (ζn ;ζt ), the friction coefficient (μ), and the limiting normal stress (σ ). As mentioned earlier, the limiting normal stress predicates on the empirical pressure-area curve in Eq. (7). The coefficient for ice-ice friction assumes a value of 0.3 based on the experimental data reported by Sukhorukov and Løset (2013), while the coefficient for ice-steel friction equals 0.1 according to the recommen­ dation in ISO19906 (2010). For impact problems, the normal damping ratio has a strong link to the coefficient of restitution (e), i.e., the ratio between the velocity after and before impact. Yulmetov et al. (2016) demonstrated that, by taking into account the near-field hydrodynamic damping together with the energy dissipated during local crushing of ice floes, the coefficient of restitution for ice floe impacts should not exceed 0.1. According to Nagurka and Huang (2004), the damping ratio and the coefficient of restitution adhere to the following relationship, ζ¼

vy = 0

y

The impact of an isolated ice floe represents the simplest form of icestructure interaction. In such a case, the contact parameters are the main factors that dictate the impact force for a given floe motion. Calibration against isolated impact data is therefore a crucial first step in deter­ mining the appropriate contact parameters. In relation to that, Timco (2011) gathered and examined numerous isolated ice floe impacts from the field measurements as well as the ice basin tests. The data comprises descriptions of the ice floe features (size, mass, and velocity) and the impact force on the structures. The study discovered a strong correlation between the impact force (Fimpact , in MN) and the floe’s kinetic energy (Ek , in GJ), which can be described using a power law, � � Fimpact ¼ 61:7ðEk Þ0:532 ; Mean (16a) bound

vx ≠ 0

Ice Floe

3.1. Simulation of isolated ice floe impacts

� Fimpact ¼ 388ðEk Þ0:532 : Upper

= 920 kg/m3

Immoveable wall

integration in the current framework adopts the Gear’s predictorcorrector technique following the original work by Matuttis and Chen (2014). To ensure stability, the simulation time-step (Δt) is equal to 10% of the estimated contact duration of the smallest discrete element, which follows, rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi minðmÞ Δt ¼ 0:1 π: (15) Eh

(17)

The parameter e in Eq. (17) defines the coefficient of restitution. Thus for e ¼ 0.1, the normal damping ratio is ζn ¼ 0.59. The tangential damping ratio, on the other hand, assumes a value of ζt ¼ ηζn , where η represents the tangential-to-normal stiffness ratio. The remaining ma­ terial property, namely the elastic modulus, is the subject of the current investigation. To validate the current approach against the impact data in Timco (2011), this study simulates numerous collisions between an arbitrarily shaped ice floe and an immovable wall. In these simulations, the water drag and the added mass effect are neglected to avoid any form of energy dissipation prior to collision. The simulation begins with the Voronoi Tessellation to produce an array of randomly shaped ice floes, which are then isolated and driven one-by-one into the target wall (see Fig. 5). By adjusting the floe size, the out-of-plane thickness, and the drift velocity, the simulation covers the range of impact force and kinetic energy

3.2. Simulation of iceberg towing model tests The validation progresses to the complex interactions within a field of ice floes. As a benchmark, the current study adopts the ice basin test simulating the iceberg towing operations in managed ice (Eik and Marchenko, 2010) under the geometric scale factor of 40. Updated de­ tails regarding the tests are available in Yulmetov and Løset (2017). The following discussions refer to the full-scale dimensions. The model test was conducted in the HSVA testing facility where the ice basin is 78 m long and 10 m wide. The experimentalists carried out the towing of a cylindrical iceberg which corresponds to the full-scale iceberg with a diameter of 76 m, a thickness of 26 m, and a density of 887 kg/m3. The towing gear consists of a line looped around the iceberg, which is then connected to a steel hawser attached to the overhead 7

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

103

Impact Force ( MN )

(a)

103

Impact Force ( MN ) DEM, E = 10 MPa

DEM, E = 1 MPa 101

101 Fimpact = 388 ( Ek )

0.532

Fimpact = 388 ( Ek )

0.532

Fimpact = 75.8 ( Ek )

0.532

10-1

Fimpact = 165.7 ( Ek )

0.532

10-1

10-3

10-3 Fimpact = 61.7 ( Ek )

0.532

10-5 -12 10 103

(b)

Fimpact = 61.7 ( Ek )

0.532

10-6 10-3 Kinetic Energy (GJ )

10-9

10-5 -12 10

100

Impact Force ( MN )

10-9

Newmans Cove Grappling Island Hans Island Cushioned Hans Island Direct Rideau Bridge Pier Molikpaq – First Year ice Molikpaq – First Year + Old ice Hans Lab Impacts NRC-CHC Lab Hondo Bridge Pier

DEM, E = 100 MPa Fimpact = 388 ( Ek )

0.532

Fimpact = 357.3 ( Ek )

0.532

10-1 10-3 Fimpact = 61.7 ( Ek )

0.532

10-9

100

(c)

101

10-5 -12 10

10-6 10-3 Kinetic Energy (GJ )

10-6 10-3 Kinetic Energy (GJ )

100

Fig. 6. Results from the isolated ice floe impact simulations versus the measured data from Timco (2011): (a) Particle elastic modulus ¼ 1 MPa; (b) particle elastic modulus ¼ 10 MPa; (c) particle elastic modulus ¼ 100 MPa.

carriage. The line, converted to full scale, is 920 m long with a diameter of 160 mm, while the steel hawser is 80 m long with a diameter of 120 mm. For each tow, the carriage travelled in a straight-line trajectory while the load cell installed on the steel hawser measured the total towing force, which includes the water and the ice resistance. The tests examined the towing force under 40%, 71%, and 86% ice concentration, while the ice thickness is fixed at 1.16 m at full scale. During the prep­ aration stage, the intact ice is cut into triangles with an approximate size of 800 m2 at full scale (assuming equilateral triangle, the sides are 43 m long). However, since the same ice field is reused for different ice con­ centrations, the floe size is not preserved throughout the tests. There­ fore, the average floe size for c ¼ 71% and 40% is somewhat smaller (Yulmetov and Løset, 2017). The benchmark test may appear unrelated to floating platforms. However, it still provides a realistic measurement of managed ice ac­ tions within a controlled setting, in which complications due to other environmental loads do not occur. Moreover, the cylindrical iceberg suits the nature of the 2D simulation, where the influence of nonprismatic geometry cannot be considered directly. Yet, there are as­ pects unique to the towing test that need to be addressed: 1) the modeling of the towing gear; and 2) the water drag acting on the iceberg.

The following paragraphs describe the modeling of the aforementioned aspects. Fig. 7a illustrates the modeling of the towing gear, which comprises two linear springs representing the towing line looped around the iceberg. As shown, the simulation excludes the steel hawser due to the lack of information on the material properties. The initial geometry depicted in Fig. 7a is an approximation based on the reported line length and the iceberg diameter. One end of the spring resides at the side of the iceberg while the other is attached to a virtual carriage, which moves according to the sequence of velocity (vtow ) shown in Fig. 7b. The movement of the carriage stretches the spring, which in turn generates tension (Ttow ) proportional to the spring stiffness (ktow ). The spring stiffness and the corresponding tension follow, ktow ¼

1.0

vtow ( m/s ) Phase 0

Phase 3 vtow = 0.82 m/s

Phase 4

Phase 2

vtow

(b) Phase 1 vtow = 0.7 m/s

0.8 84.6°

(18b)

where Et and At denote the elastic modulus and the cross-sectional area of the towing line, respectively, while lt represents the un-stretched spring length. The current framework computes the towing line

ktow Iceberg

(18a)

Ttow ¼ ktow Δlt � 0;

(a)

76 m

Et At 95 � 109 *π � 0:082 ¼ ¼ 4:78 � 106 N=m; lt 400

0.6 0.4 0.2 0

0

10

20 30 Time ( min )

40

Fig. 7. Setup for of the iceberg towing simulations: (a) the modeling of the towing line, shown at the initial condition; and (b) the sequence of the towing velocity. 8

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

tension at every time-step, and subsequently applies the resultant tow­ ing forces and moment onto the iceberg, similar to the procedure for the mooring restoring forces (see Eq. [8]). Prior to the main test, the experiment measured the water resistance on the iceberg through an open-water towing. Yulmetov and Løset (2017) showed that the water resistance oscillates due to the complex flow around the cylindrical geometry. Assuming that the water resis­ tance is dictated mainly by the form drag mechanism, Yulmetov and Løset (2017) proposed the following expression to describe the water drag on the iceberg, � � Cf Fdrag; iceberg ¼ (19) h⋅D ρw jVw vjðVw vÞ 2

causing a lateral confinement effect that intensifies the ice forces. The time-series in Fig. 10 demonstrate that the simulated towing force in high concentration (c ¼ 71% and 86%) consists of sustained peaks that are considerably higher than the transient peaks observed under lower ice concentration (c ¼ 40%). However, compared to the test, the sim­ ulations for c ¼ 71% and 86% underestimate the towing force signifi­ cantly. The underestimation becomes more severe as the ice concentration increases, as depicted in Fig. 10b–c. These comparisons, once again, demonstrate the limitation of a 2D approach. The nature of a densely packed ice field, aggravated by the narrow and confining basin, unavoidably leads to 3D motions as well as ice floe fragmentations. These phenomena ultimately result in the accumulation of ice floes, not only in the horizontal plane, but also in the vertical direction. As the ice starts to accumulate at the interface between the ice floe and the structure, the effective thickness of the ice floe increases due to the in­ crease in the contact area between the ice floe and the structure. This leads to an increase in the ice load, which cannot be represent-ed in the 2D simulation. The 2D framework simply lacks the means to simulate the three-dimensional ice build-up, and thus leads to the force under­ estimation in Fig. 10b–c. The comparisons in Fig. 10 prove that the current 2D DEM works reasonably well under limited circumstances, i.e., when the ice con­ centration throughout the simulation remains low enough for the 2D assumption to stay valid, which was also confirmed in Hopkins and Tuhkuri (1999). The current results suggest that two factors dictate this condition: 1) the initial ice concentration, and 2) the relative domain width (which affects the lateral confinement). Nonetheless, defining the exact limits for these two parameters is challenging due to the complexity of the problem and the limited empirical data. From the current tests alone, the initial ice concentration of c ¼ 40% and the domain of width 400 m appear to be the limiting conditions. As an attempt to further investigate the effect of domain width, this study simulates several hypothetical cases where the basin width varies while the mean floe size and the initial concentration are maintained. A newly introduced parameter, namely the normalized domain width (ψ ), quantifies the net domain width relative to the average floe size, and follows, �� ψ ¼ Wdomain Diceberg Dfloe : (20)

where h and D correspond to the thickness and the diameter of the iceberg, respectively. The experimental data indicates that the form drag coefficient fluctuates around Cf ¼ 0.37 throughout the tow. The current simulation adopts this coefficient under the assumption that the average drag force represents the water resistance adequately. Table 1 summarizes the parameters assigned for the simulation, while Fig. 8 presents the ice floe size distribution for each ice concen­ tration considered. Since it is not possible to precisely reproduce the ice floes from the tests, the simulated ice field enforces only the key char­ acteristics, namely the average floe size and the ice concentration. The floe size reduction reflected in Fig. 8 is merely an estimate, since the actual value is not reported. Fig. 9 shows snapshots from the simulation, while Fig. 10 compares the simulated versus the measured towing force. In a lower ice concentration (c ¼ 40%), the ice floes have plenty of space to disperse around the advancing iceberg (see Fig. 9a). As a result, the ice actions come mostly in the form of isolated or simultaneous collisions with little to no floe accumulation ahead of the iceberg. The simulated towing force (Fig. 10a) reflects these collisions as transient peaks throughout the towing. However, the water resistance on the iceberg still constitutes a significant portion of the total towing action, as indicated by the relatively high force in between the force peaks. Compared to the measured towing force (Fig. 10a), the simulation in­ dicates a reasonable agreement although it is unable to reproduce the exact towing force evolution. The transient peaks are somewhat com­ parable, yet the simulation underestimates the average towing force due to the simplification in the iceberg water resistance. This comparison therefore suggests that the current DEM is able to estimate the ice collision forces, and with improved water resistance formulation, may reproduce the measured the towing force more accurately. As the ice concentration increases, the chain of contacts between ice floes gradually becomes the dominant driving factor in the total ice action. The snapshots from the case with c ¼ 71% (Fig. 9b) and c ¼ 86% (Fig. 9c) clearly portray this phenomenon, which further suggest a force build-up as opposed to the collision-induced force peaks. In addition to that, the contact network extends laterally towards the basin walls,

where Wdomain represents the domain width, Diceberg denotes the iceberg

diameter, and Dfloe corresponds to the average floe size. Fig. 11 presents the evolution of the maximum towing force with respect to the normalized domain width (ψ ) at c ¼ 40%, 71%, and 86%. As shown, the maximum towing force converges when the parameter ψ is greater than a certain critical value, depending on the initial ice concentration. When c ¼ 40%, the towing force converges at ψ ¼ 8.0, which coincidentally matches the actual ice basin width. At higher initial concentrations (c ¼ 71% and 86%), the towing force converges at ψ values larger than that of the actual basin. Conversely, when ψ is less than the critical value, the maximum towing force increases rapidly as a consequence of the lateral confinement effect, which promotes severe accumulation of ice floes. Nonetheless, due to the limitation of the 2D assumption, the maximum towing force at a small ψ may not be accurate, as is the case for c ¼ 71% and 86% (see Fig. 11). Overall, the results in Fig. 11 confirm that the lateral confinement effect becomes negligible past a certain ψ value. If the initial ice con­ centration is sufficiently low, this critical ψ can be a useful predictor whether the ice concentration will remain low enough for a valid 2D assumption. Conversely, if the initial ice concentration is high, there can be a critical ψ value, yet the 2D assumption is questionable from the start of the simulation. Based on the current benchmark tests, the proposed 2D discrete element framework is reasonably accurate up to c ¼ 40% with ψ > 8.0. The exact validity range may correspond to a slightly larger c and ψ value, yet they will remain unknown due to the lack of different ice concentrations in the test program. In conclusion, beyond

Table 1 Parameters for the iceberg towing simulation. Parameter

Description

Values

Unit

hice Vw E

Ice floe thickness Water current velocity Particle elastic modulus Tangential-to-normal property ratio Normal damping ratio Tangential damping ratio Form drag coefficient Skin drag coefficient Rotational drag coefficient Form drag coefficient for the iceberg Friction coefficient for ice-ice contact Friction coefficient for ice-boundary contact Ice floe density Iceberg density Water density

1.16 0.0 1.0 0.38 0.59 η ζn 1.0 0.002 1.0 0.37 0.3 0.05 920 887 1025

m m/s MPa – – – – – – – – – kg/m3 kg/m3 kg/m3

η ζn ζt Cf Cs Cr Cf, iceberg

μice-ice μice-boundary ρice ρiceberg ρw

9

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

1.2

Normalized Frequency

(a)

1.2

c = 40% Dfloe = 40.3 m

1.0

0.8

0.6

0.6

0.4

0.4

0.2

0.2

1.2

0

10

20

30

40

50

Dfloe ( m )

60

70

80

Normalized Frequency

0

90

(b) c = 71%

1.0

0.8

0

Normalized Frequency

Dfloe = 41.8 m

0

10

20

30

40 50 Dfloe ( m )

60

70

80

90

(c) c = 86% Dfloe = 43.7 m

1.0 0.8 0.6 0.4 0.2 0

0

10

20

30

40 50 Dfloe ( m )

60

70

80

90

Fig. 8. The distribution of the ice floe size in the iceberg towing simulation: (a) for 40% initial ice concentration; (b) for 71% initial ice concentration; and (c) for 86% initial ice concentration.

designed specifically to endure harsher ice condition, such that it can extend the annual drilling season up to the early winter. The platform’s instrumented mooring lines have amassed invaluable ice load data during the years of deployment in the Beaufort Sea. To date, the collected data is regarded as the most comprehensive reference since no other field measurements have covered the range of ice conditions documented by the Kulluk. Fig. 12a illustrates the Kulluk during one of its deployments. As reported in Wright (2000), Kulluk has a waterline diameter of 70 m, an operating draft of 11.5 m, and a displacement of 28,000 tonnes. Its mooring system consists of 12 catenary steel wires that are evenly-spread, each with a diameter of 90 mm. During opera­ tion, the mooring lines’ length and pretension vary according to the site water depth. An example data in Wright (2000) showed that, for a water depth of 32 m, the mooring lines are 530–732 m long while the pre­ tension ranges from 600 to 2200 kN. Fig. 12b illustrates the modeling of the hypothetical platform while Table 2 summarizes all of the relevant properties. Following the Kulluk, the hypothetical platform has a circular water plane geometry with a diameter of 70 m. The springs illustrated in Fig. 12c represent the hor­ izontal projection of the mooring system. Each spring obeys the tensionoffset relationship shown in Fig. 12c, which is obtained from a separate 3D static-offset analysis performed in HARP (Hull And Riser/mooring Program) – an industry-accepted coupled hull-mooring analysis package (Kim et al., 2001; Tahar and Kim, 2003). The HARP model assumes that all 12 mooring lines are of equal length and pretension, while the ma­ terial properties are taken from a steel wire catalogue in the current market (see Table 2 for the values). As a simple check, Fig. 13a contrasts the static-offset properties of the mooring system modeled in the discrete element framework against that in HARP. The comparison shows a close agreement, which verifies the mooring system representation in the 2D DEM. The case study examines the effect of various managed ice conditions on the ice load as well as the mooring load. As shown in Table 3, the

Fig. 9. Snapshots from the iceberg towing simulation: (a) for 40% initial ice concentration; (b) for 71% initial ice concentration; and (c) for 86% initial ice concentration.

the verified conditions, the 2D DEM may only provide qualitative pre­ diction of the actual ice action. 4. Managed ice actions on moored floating platforms 4.1. General description This section presents a case study on managed ice load on a moored floating platform as predicted by the proposed discrete element frame­ work. The subject of the study is a hypothetical floater inspired by the Kulluk. As background information, Kulluk is a conical floating structure 10

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

Fig. 10. The simulated versus the measured towing forces: (a) for 40% initial ice concentration; (b) for 71% initial ice concentration; and (c) for 86% initial ice concentration.

9

DEM, c = 40% DEM, c = 71% DEM, c = 86%

while others (Molyneux et al., 2011; Aksnes, 2011) showed otherwise. Accordingly, in addition to the “default” mooring stiffness in Fig. 13a, the case study simulates a “fixed” and a drastically more flexible system to examine the significance of structural compliance. In the upcoming section, the flexible system corresponds to 5% of the “default” stiffness.

6

Test, c = 40% Test, c = 71% Test, c = 86%

4.2. Simulation setup

12

Peak Towing Force ( MN )

For every case in Table 3, this study assigned a 1560 � 1180 m2 ice domain in which the floes are generated. Based on previous experiences from the narrow ice basin tests, the current simulation area is set to be considerably larger than the ice domain in order to avoid the effect of artificial lateral confinement (ψ > 50). Nonetheless, the case study purposely adopts a relatively wide ice domain (�17 times floater diameter) to simulate the natural confining effect from the mass of densely packed ice floes. Table 4 summarizes the parameters used in the simulation, which mostly remain the same as those in Section 3. In the case study, the ice floes are riding on the constant water current that moves towards the floating platform. The ice drift velocity therefore depends on the current velocity (Vw ) as implied by the drag force formulation. On the other hand, the case study ignores the water drag acting on the floater and focus exclusively on the ice-induced loads. In support of this treatment, a recent DEM study (Karulin and Karulina, 2013) has demonstrated that for a 50% ice concentration and a 0.5 m/s current velocity, the mooring

3 0 0

2

4

6

8

10 12 14 16 18 20 22

, Normalized Domain Width Fig. 11. The effect of domain width on the simulated ice action.

examined parameters consist of the ice concentration, the ice thickness, the ice drift velocity, and the ice floe size. The values for these param­ eters predicate on the ranges and descriptions given in Kulluk’s field report (Wright, 2000). The case study also examines the effect of the global mooring stiffness. Past studies have shown that the compliance afforded by the mooring system affects the ice actions to some extent. Some researches (Lawrence, 2009; Sayed and Barker, 2011) demon­ strated that an increasing structural compliance reduces the ice action, 11

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

Fig. 12. The hypothetical floating platform in the case study: (a) the basis for the hypothetical floater, the Kulluk; (b) the modeling of the floating platform; and (c) the tension-offset curve assigned to the springs.

motion comes exclusively from the mooring system damping (see Eq. [9]), yet with adjustments on the damping ratio to compensate the omitted hydrodynamic damping. The chosen damping ratio of ζmoor ¼ 4% derives from calibration against the free-decay analysis in HARP, in which the hydrodynamic damping on the hull and the mooring lines are taken into account. As confirmed in Fig. 13b, the 2D floater-mooring system with a 4% damping produces comparable mo­ tion decay to that from the 3D model in HARP.

Table 2 Parameters of the hypothetical floater and its mooring system. Parameters

Values

Unit

Water depth Floater waterline diameter Floater mass (mflt) Floater mass moment of inertia (Jflt) Mooring line diameter Mooring line length Mooring line pretension Mooring line axial rigidity (EA) Mooring line dry weight Mooring line wet weight Translational stiffness of mooring system (Kx, Ky) Rotational stiffness of mooring system (Krot)

32 70 2.80 � 107 1.64 � 109 90 650 2000 5.0 � 105 235 205 3.13 � 106 6.35 � 108

m m kg kg.m2 mm m kN kN kg/m kg/m N/m Nm/rad

4.3. Results and discussions 4.3.1. Simulation snapshots and typical load time series The snapshots in Fig. 14 exemplify the ice floe motion in two dras­ tically different ice concentrations, namely c ¼ 20% (Fig. 14a–b) and c ¼ 90% (Fig. 14c–d), subjected to the same ice velocity. The same ice velocity in the two different ice concentration implies the different waves, currents and winds experienced by the floater. Fig. 14 aims to examine the effect of ice concentration on the ice load, while

load due to water drag is insignificant compared to that due to the ice actions. Excluding the water drag on the floater also implies the omis­ sion of the hydrodynamic damping. Currently, the decay of the floater

15

Global Mooring Load ( MN )

(a)

4

DEM HARP

Offset ( m )

(b) DEM HARP

2

10

0 5

0

-2

Static-offset analysis Water depth = 32 m 0

1

2

Offset (m)

3

4

-4

5

Free-decay analysis Water depth = 32 m 0

200

400

600

Time (sec)

800

1000

Fig. 13. Response of the floating platform model in the DEM framework and in HARP: (a) from a static-offset test; and (b) from a free-decay test. 12

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

difficulties for a meaningful quantitative comparison. The field data in Wright (2000) consists of numerous load events, during which the mooring system configuration varies. The report, however, contains no description on the mooring system for the corresponding ice load events. Moreover, the overall ice conditions are described qualitatively using the “ice regime” parameter, which is proportional to the ice concen­ tration yet does not reflect the specific ice thickness, drift velocity, or floe size. From the simulation side, the current 2D assumption restricts the floater geometry to prismatic shapes, which means that the effect of the hull’s downward slope is not present in the simulations. Due to the discrepancy in the floater geometry, in addition to the uncertain ice and mooring characteristics, the current comparison can only verify the correctness of the simulated load in a qualitative manner, i.e., in terms of the overall trend and the order of magnitude. Fig. 16 contrasts the simulated global load against that from the field measurements. As shown, the field data scatters significantly, yet it in­ dicates a roughly exponential trend with respect to the ice regime. In this comparison, the simulation adopts the typical ice encountered by the Kulluk (Wright, 2000), which is characterized by a thickness of 1.0 m, a drift velocity of 0.5 m/s, an average floe size of 50 m, and an ice con­ centration that varies from 20% to 90%. Furthermore, the comparison predicates on the assumption that simulated ice concentration is equal to the ice regime. Fig. 16 shows that the discrete element simulation is able to reproduce the exponentially increasing mooring load. The simulation, with all of the assumed ice condition, also produces global loads that are in reasonable agreement with the mean curve of the field data. This qualitative agreement, at the very least, establishes some credibility for the simulated load in the upcoming case study. However, as demon­ strated through the towing simulation in Section 3, the 2D framework carries intrinsic limitations that compromise the accuracy in higher ice concentration. Although results from the higher ice concentrations may not be accurate, they are still included for the sake of completeness.

Table 3 Conditions examined in the hypothetical floater case study. hice (m)

Vice (m/s)

Mean Floe Size (m)

Concentration (%)

Mooring System

Ice concentration effect Floe size effect

1.0

0.5

50

20, 30, 40, 50, 60, 70, 80, 90

"Default"

1.0

0.5

50

"Default"

Ice thickness effect

1.0, 2.0, 4.0 1.0

0.5

30, 50, 80 50

50

"Default"

50

50

"Default"

50 50 50

50, 90 50, 90 50, 90

"Default" "Fixed" "5% Default"

Ice drift velocity effect Mooring stiffness effect

1.0 1.0 1.0

0.25, 0.5, 1.0 0.5 0.5 0.5

Table 4 DEM parameters for the hypothetical floater case study. Parameter

Description

Values

Unit

E

Particle elastic modulus Tangential-to-normal property ratio Normal damping ratio Tangential damping ratio Form drag coefficient Skin drag coefficient Rotational drag coefficient Friction coefficient for ice-ice contact Friction coefficient for ice-floater contact Ice floe density Water density Mooring damping ratio

1.0 0.38 0.59 η ζn 1.0 0.002 1.0 0.3 0.05 920 1025 0.04

MPa – – – – – – – – kg/m3 kg/m3 –

η ζn ζt Cf Cs Cr

μice-ice μice-floater ρice ρw ζmoor

4.3.3. Effect of managed ice characteristics Fig. 17 illustrates the effect of various ice characteristics on the mooring and the ice load. Since the typical global load consists of distinct peaks (see Fig. 15), the case study places emphasis on the maximum instead of the average load. In general, the global mooring load demonstrates the same trend as the global ice load, although their magnitudes are not necessarily equal. Fig. 17a illustrates the global load evolution under varying ice con­ centrations. While both the mooring and ice loads increase exponen­ tially with ice concentration, the mooring load noticeably exceeds the ice load at c ¼ 40% and 50%, thus indicating a critical dynamic amplification. As briefly discussed earlier, this phenomenon occurs due to the interaction between the ice load and the floater’s dynamic properties, which coincidentally resonates under the aforementioned ice concentrations. As also illustrated in Fig. 17a, the ice load surpasses the mooring load at higher ice concentrations (c � 80%). Again, the reason behind this trend lies on the dynamics of the global loads at the considered ice concentration. Fig. 17b presents the effect of the average floe size (Dfloe ) on the global load while the ice thickness, the ice drift velocity, and the ice concentration are maintained. The simulated loads indicate a linearly increasing trend with respect to the average floe size. Considering the proportionality between the impact load and the floe kinetic energy (see Fig. 6), this trend is expected since the energy level rises with the ice floe size. However, the results in Fig. 17b demonstrate that the global load is not very sensitive to the average floe size. As the floe size increases from 30 m to 80 m (�167% increase), the global ice load and the global mooring load escalate by only 83% and 82%, respectively. Despite the seemingly conclusive results, the simulation assumes unbreakable ice floes, which is no longer realistic when the floe is large relative to its thickness. The current framework by itself is incapable of establishing this floe size limit. As observed from the recent literatures, calibrated 3D

maintaining the other ice load parameters. Meanwhile, Fig. 15 presents the global ice load, the global mooring load, as well as the offset timeseries from the aforementioned cases. The global load in this context refers to the resultant of the x and y force components. In a low ice concentration, the global ice load consists of sparsely distributed impulsive forces, which cause the floater to undergo free vibration in between impacts. This is evident from the oscillating mooring load in between the ice load peaks (Fig. 15a) as well as the oscillating floater offset (Fig. 15b). Conversely, the higher ice concentration generates greater and longer-sustained ice load (Fig. 15c), which in turn lessens the floater’s free vibration response (Fig. 15d). The load time-series also demonstrate that, at times, the mooring load is greater than the corresponding ice impact load. These events indicate the presence of dynamic amplification, which depends heavily on the dynamic properties of both the ice load and the floating platform. The forces in Fig. 15a suggest that the dynamic amplification tend to occur when the ice load duration is comparable to the platforms’ oscillation period. Further post-processing of Fig. 15a reveals that the average ice load duration is indeed close to the natural period of the floater. The ice load duration fluctuates around 19.9 s while the floater’s natural period is 18.78 s (based on the simple estimation, Tn ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2π mflt =Kx ). Despite the dynamic amplifications, the overall maximum of the ice and the mooring load in Fig. 15a remains comparable. This may not be the case when the ice condition or the mooring stiffness changes. However, the ice-induced dynamic amplification is a complex topic that warrants a separate investigation. The current study, due to the inherent restrictions of the 2D framework, will not address this phenomenon in a detailed manner. 4.3.2. Comparison to field data on Kulluk Limitations on both the field data and the simulation bring 13

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

Fig. 14. Snapshots from the floating platform simulations: (a) 20% initial concentration, at the start of simulation; (b) 20% initial concentration, at mid-point of simulation; (c) 90% initial concentration, at the start of simulation; and (d) 90% initial concentration, at mid-point of simulation.

Global Load ( MN )

0.75

c = 20% Dfloe = 50 m Vw = 0.5 m/s hice = 1.0 m

Floater Offset ( m )

(a)

(b)

0.2

Ice Load Mooring Load

0.1

0.50

0.0 0.25

0.00

-0.1

0

1000

Time (sec )

2000

-0.2

3000

Global Load ( MN )

1000

Time (sec )

0.4

c = 90%

3000

(d)

0.2

Vw = 0.5 m/s

Ux Uy

2000

Ux Uy

Ice Load Mooring Load

1.5

c = 90% Vw = 0.5 m/s Dfloe = 50 m hice = 1.0 m

Dfloe = 50 m

1.0

hice = 1.0 m

0.0

0.5 0.0

0

Floater Offset ( m )

(c)

2.0

c = 20% Dfloe = 50 m Vw = 0.5 m/s hice = 1.0 m

0

1000

Time (sec )

2000

-0.2

3000

0

1000

Time (sec )

2000

3000

Fig. 15. The example time series from the floating platform simulations: (a) the global load and (b) the floater offset for c ¼ 20%; (c) the global load and (d) the floater offset for c ¼ 90%.

14

M. Raditya Pradana et al.

8

Ocean Engineering 190 (2019) 106483

Global Load ( MN ) Kulluk, mooring load data DEM mooring load, max. DEM ice load, max.

6

the velocity term is squared. The results in Fig. 17 confirm that the ice-induced actions demon­ strate stable and clear trends with respect to the managed ice charac­ teristics. Among the investigated parameters, the ice concentration appears to be the most complex determining factor. The increase in ice concentration leads to ice floe accumulation/jamming, which not only produces higher loads, but also changes the dynamics of the load as exemplified earlier in Fig. 15. On the other hand, the effect of the remaining parameters, namely the ice thickness, the floe size, and the ice drift velocity, are predictable through their contribution in the kinetic energy of the moving ice floe. Overall, under the current unbreakable floe assumption, the global load indicates the least sensitivity to the average ice floe size.

DEM

Vw = 0.5 m/s Dfloe = 50 m hice = 1.0 m

4

Field data upper-bound Field data mean

2 0

2

3

4

5

6

7

8

Ice Regime

9

4.3.4. Effect of mooring stiffness Fig. 18 presents the maximum global loads simulated under three different mooring systems. The “fixed” system is equivalent to a system with an infinite stiffness, under which the global mooring load is equal to the global ice load. The “default” system refers to that described in Fig. 12, while the “5% default” system maintains the same geometrical configuration, yet with tension-offset curve decreased to 5% of that in Fig. 12c. The drastic loss in stiffness undoubtedly leads to a much higher offset, which likely exceeds the maximum allowable limit (5% of water depth for drilling). Since the ice concentration changes the ice load dynamics, the case study also examines the mooring stiffness effect under two different ice concentrations, i.e., c ¼ 50% and 90%. In terms of the ice load magnitude, the variation in the mooring stiffness leads to some changes. Under c ¼ 50% (Fig. 18a), the difference in the ice load is practically negligible as compared to the dramatic change in the mooring stiffness. As the stiffness decreases from “fixed” to “default”, the ice load increases by 13.3%. Further stiffness reduction from “default” to “5% default”, on the contrary, lowers the ice load by 9.6%. When c ¼ 90% (Fig. 18b), the mooring stiffness effect is some­ what more substantial. The ice load indicates a mere 5% increase when

10

Fig. 16. The simulated ice and mooring load compared to the field data on the Kulluk, adapted from Wright (2000).

DEM with bonded particles (e.g., in Di et al., 2017) is arguably the more suitable tool to address larger ice floes. The global load increases with the ice thickness (Fig. 17c) and the ice drift velocity (Fig. 17d), as expected from their influence on the kinetic energy. The result demonstrates that increasing the ice thickness from 1.0 m to 4.0 m (300% increase) leads to a 239% and a 185% increase in the global ice load and the global mooring load, respectively. In contrast, as the ice drift velocity increases from 0.25 m/s to 1.0 m/s (300% in­ crease), the global ice load and the global mooring load sharply increase by 750% and 605%, respectively. The higher sensitivity towards the ice drift velocity is consistent with the definition of kinetic energy, in which

3.00

Global Load ( MN )

(a)

3.00

DEM mooring load, max. DEM ice load, max.

2.25

0.75

0

20

30

40

50

60

70

80

Ice concentration (%)

0

90 100

Global Load ( MN )

(c)

20

40

60

Mean floe size (m)

80

Global Load ( MN )

2.25

Drift velocity = 0.5 m/s Concentration = 50% Mean floe size = 50 m

1.50

0

3.00

DEM mooring load, max. DEM ice load, max.

2.25

1.50

100

(d)

DEM mooring load, max. DEM ice load, max. Mean floe size = 50 m Ice thickness = 1.0 m Concentration = 50%

0.75

0.75 0

Drift velocity = 0.5 m/s Ice thickness = 1.0 m Concentration = 50%

1.50

0.75

3.00

(b)

DEM mooring load, max. DEM ice load, max.

2.25

Drift velocity = 0.5 m/s Ice thickness = 1.0 m Mean floe size = 50 m

1.50

Global Load ( MN )

0

1

2

3

Ice thickness (m)

4

0

5

0

0.3

0.6

Ice drift velocity (m/s)

0.9

1.2

Fig. 17. The effect of managed ice characteristics on the global ice and mooring load: (a) ice concentration; (b) average ice floe size; (c) ice thickness; and (d) ice drift velocity. 15

M. Raditya Pradana et al.

1.5

Ocean Engineering 190 (2019) 106483

Global Load ( MN ) DEM mooring load, max. DEM ice load, max.

1.0

(a)

3

Concentration = 50% Drift velocity = 0.5 m/s Ice thickness = 1.0 m Mean floe size = 50 m

DEM mooring load, max. DEM ice load, max.

2

(b) Concentration = 90% Drift velocity = 0.5 m/s Ice thickness = 1.0 m Mean floe size = 50 m

1

0.5

0.0

Global Load ( MN )

Fixed

Default

Mooring Stiffness

0

5% Default

Fixed

Default

Mooring Stiffness

5% Default

Fig. 18. The effect of mooring stiffness on the global ice and mooring load: (a) at c ¼ 50%; and (b) at c ¼ 90%.

the stiffness is reduced from “fixed” to “default”, while stiffness reduc­ tion from “default” to “5% default” lowers the ice load by 25.5%. Overall, Fig. 18 demonstrates that the mooring stiffness influence on the ice impact load is relatively minor. The results suggest that, for a given managed ice characteristics, there exists a mooring stiffness that will maximize ice load, which should be avoided in design. On the other hand, considerable mooring stiffness reduction alleviates only a small portion of the ice load, and is unfeasible for drilling or other critical operations. In a broader context, Fig. 18 suggests that including the floater motion in an ice-structure interaction analysis do not alter the computed ice load significantly. In practice, this implies that the dy­ namic mooring-floater analysis may be conducted separately using specialized tools, by taking the DEM ice load (with fixed structure assumption) as input. The change in the mooring stiffness undoubtedly influences the global mooring load. However, the influence seemingly depends on the ice load dynamics as well. The ice concentration of c ¼ 50% have pre­ viously demonstrated the tendency to trigger dynamic amplification, which explains the 36% increase in the global mooring load as the mooring stiffness changes from “fixed” to “default” (see Fig. 18a). On the contrary, further stiffness reduction from “default” to “5% default” alters the dynamic properties of the floater relative to that of the ice load, which ultimately results in 48% reduction of the global mooring load. The same reasoning explains the trend under c ¼ 90% (Fig. 18b), except that the corresponding ice load do not have a tendency to induce dy­ namic amplifications.

isolated ice floe impact database (Timco, 2011) verifies that the most suitable elastic modulus for the current formulation is 1 MPa. Despite being several orders lower than the actual sea ice modulus (2–6 GPa), the chosen modulus yielded the best agreement with the average of the impact database. The seemingly unphysical modulus is likely a conse­ quence of the soft-particle principle, where the particle-particle overlap is a crude estimation of the actual deformation. 2) The 2D discrete element framework simulates a reasonably ac­ curate ice action under a certain condition, i.e., when the ice concen­ tration throughout the simulation remains sufficiently low. The validation against the iceberg towing model test reveals that agreement is achieved only when the initial ice concentration is 40%. Meanwhile, for the 71% and 86% initial concentrations, the discrete element anal­ ysis underestimates the measured ice load significantly. The underesti­ mation occurs since the 2D framework is incapable of simulating the three-dimensional ice buildup that tends to occur in high initial con­ centration, or under confined environment (e.g. a narrow ice basin). Based on the currently available results, the validity of the 2D formu­ lation is verified up to the initial ice concentration of c ¼ 40%, with a relative domain width of ψ ¼ 8.0. 3) The case study on the hypothetical floater presents the effect of managed ice characteristics on the ice load as well as the mooring load. The overall results demonstrate that the mooring load follows the same trend as the ice load. However, depending on the dynamics of the ice load relative to that of the mooring system, the mooring load may exceed the ice action as a result of dynamic amplification. The global load (refers to both the ice and the mooring load) increases in an approximately exponential trend with respect to the initial ice concen­ tration. With respect to the other ice characteristics, namely ice thick­ ness, floe size, and ice drift velocity, the global load indicates a roughly linear trend over the examined parameters. Among all ice characteris­ tics, the global load is the least sensitive to the average floe size. 4) In addition, the case study shows that the structural compliance offered by the mooring system leads to relatively insignificant variations in the ice load. The ice load slightly increases as the mooring stiffness shifts from infinite (“fixed”) to the “default” value assumed in the case study. In contrast, an additional stiffness reduction from “default” to “5% default” lowers the ice load by relatively small margin. The observed trend suggests that the ice load reaches a maximum under a certain mooring stiffness, which should be avoided in design.

5. Summary and conclusions This paper proposes an efficient 2D discrete element approach as a tool to assess managed ice load on moored floating platforms. The discrete element framework adheres to the soft-particle principle, yet with modifications to facilitate non-disk elements. In the proposed tool, the ice floes are arbitrarily-shaped polygons created using the Voronoi Tessellation algorithm, which enables direct control over the ice con­ centration and the average floe size. The floating platform is represented using a polygon that possesses the platform’s waterline geometry, mass, and mass moment of inertia. Meanwhile, the individual mooring line is a spring that follows a tension-offset relationship, generated using a specialized mooring analysis package. This paper validates the discrete element simulation against relevant experiments and field measure­ ments, and subsequently demonstrates the application through a case study on a hypothetical floating platform based on the Kulluk. The re­ sults presented throughout the paper support the following conclusions: 1) The elastic modulus assigned to the particles/discrete elements is critical to the estimated impact force. The comparison against the

Acknowledgement The authors thank the National Research Foundation, Keppel Cor­ poration and National University of Singapore (R-261-507-010-281) for supporting this work done in the Keppel-NUS Corporate Laboratory. The 16

M. Raditya Pradana et al.

Ocean Engineering 190 (2019) 106483

conclusions put forward reflect the views of the authors alone, and not necessarily those of the institutions within the Corporate Laboratory.

Lau, M., 2006. Discrete element modeling of ship manoeuvring in ice. In: 18th International Symposium on Ice, Sapporo, JAPAN (IAHR Ice Symposium 2006), Aug. Lau, M., Simoes R� e, A., 2006. Performance of survival craft in ice environments. In: 7th Int. Conf. And Exhibition on Performance of Ships and Structures in Ice, ICETECH 2006, Banff, Canada, pp. 51–58. Alexandria, VA. Lau, M., Lawrence, K.P., Rothenburg, L., 2011. Discrete element analysis of ice loads on ships and structures. Ships Offshore Struct. 6 (3), 211–221. Lawrence, P., 2009. Load Prediction for a Moored Conical Drillship in Level Unbroken Ice: a Discrete Element and Experimental Investigation. PhD thesis. University of Waterloo. Liu, L., Ji, S., 2018. Ice load on floating structure simulated with dilated polyhedral discrete element method in broken ice field. Appl. Ocean Res. 75, 53–65. Løset, S., 1994. Discrete element modelling of a broken ice field—Part I: model development. Cold Rons Sci Tech 22 (4), 339–347. Løset, S., 1994. Discrete element modelling of a broken ice field—Part II: simulation of ice loads on a boom. Cold Rons Sci Tech 22 (4), 349–360. Løset, S., Kanestrom, O., Pytte, T., 1998. Model tests of a submerged turret loading concept in level ice, broken ice and pressure ridges. Cold Rons Sci Tech 271, 57–73. Lu, P., Li, Z., Cheng, B., Lepp€ aranta, M., 2011. A parameterization of the ice-ocean drag coefficient. J Geophys Res-Oceans 116 (C7). Matskevitch, D.G., 1997. Eccentric impact of an ice feature: linearized model. Cold Rons Sci Tech 25 (3), 159–171. Matuttis, H.G., Chen, J., 2014. Understanding the Discrete Element Method: Simulation of Non-spherical Particles for Granular and Multi-Body Systems. John Wiley & Sons. Meylan, M.H., Yiew, L.J., Bennetts, L.G., French, B.J., Thomas, G.A., 2015. Surge motion of an ice floe in waves: comparison of a theoretical and an experimental model. Ann. Glaciol. 56 (69), 155–159. Molyneux, D., Derradji-Aouat, A., Cholley, J.M., 2011. First-year ridge loads on moored offshore structures. In: OTC Arctic Technology Conference. Offshore Technology Conference, Houston, TX, OTC-22099-MS. https://doi.org/10.4043/22099-MS. Nagurka, M., Huang, S., 2004. A mass-spring-damper model of a bouncing ball. In: Proceedings of the 2004 American Control Conference, vol. 1. IEEE, pp. 499–504. Okabe, A., Boots, B., Sugihara, K., Chiu, S.N., 2009. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams, vol. 501. John Wiley & Sons. Paulin, M., 2008. Arctic Offshore Technology Assessment of Exploration and Production Options for Cold Rons of the US Outer Continental Shelf. In: Prepared by IMV for American Bureau of Ocean Energy Management, Regulation and Enforcement, IMVPA Project No C-0506-15. Pradana, M.R., Qian, X., 2019. Bridging local parameters with global mechanical properties in bonded discrete elements for ice load prediction on conical structures. Cold Rons Sci Tech (submitted for publication). Rabatel, M., Labb� e, S., Weiss, J., 2015. Dynamics of an assembly of rigid ice floes. J. Geophys. Res.: Oceans 120 (9), 5887–5909. Sayed, M., Barker, A., 2011. Numerical simulations of ice interaction with a moored structure. In: OTC Arctic Technology Conference. Offshore Technology Conference, Houston, TX, OTC-22101-MS. https://doi.org/10.4043/22101-MS. Sotomayor, O.E., Tippur, H.V., 2014. Role of cell regularity and relative density on elasto-plastic compression response of random honeycombs generated using Voronoi diagrams. Int. J. Solids Struct. 51 (21–22), 3776–3786. Sukhorukov, S., Løset, S., 2013. Friction of sea ice on sea ice. Cold Rons Sci Tech 94, 1–12. Sun, S., Shen, H.H., 2012. Simulation of pancake ice load on a circular cylinder in a wave and current field. Cold Rons Sci Tech 78, 31–39. Tahar, A., Kim, M.H., 2003. Hull/mooring/riser coupled dynamic analysis and sensitivity study of a tanker-based FPSO. Appl. Ocean Res. 25 (6), 367–382. Timco, G.W., 2011. Isolated ice floe impacts. Cold Rons Sci Tech 68 (1–2), 35–48. Tuhkuri, J., Poloj€ arvi, A., 2018. A review of discrete element simulation of ice–structure interaction. Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 376 (2129). Wang, J., Sayeed, T., Millan, D., Gash, R., Islam, M., Millan, J., 2016. Ice model tests for dynamic positioning vessel in managed ice. In: Arctic Technology Conference. Offshore Technology Conference, Houston, TX, OTC-27430-MS. https://doi. org/10.4043/27430-MS. Williams, M., 1989. Nekton 8000 Model in Ice: Fixed Structure and Moored Structure. Technical report TR-1989-07 prepared by the Institute of Marine Dynamics. Wright, B., 2000. Full Scale Experience with Kulluk Station Keeping Operations in Pack Ice (With Reference to Grand Banks Developments), pp. 25–44. Report no. PERD/ CHC Report. https://doi.org/10.4224/12327366. Yulmetov, R., Løset, S., 2017. Validation of a numerical model for iceberg towing in broken ice. Cold Rons Sci Tech 138, 36–45. Yulmetov, R., Lubbad, R., Løset, S., 2016. Planar multi-body model of iceberg free drift and towing in broken ice. Cold Rons Sci Tech 121, 154–166. Zhan, D., Agar, D., He, M., Spencer, D., Molyneux, D., 2010. Numerical simulation of ship maneuvering in pack ice. In: ASME 2010 29th International Conference on Ocean, Offshore and Arctic Engineering. American Society of Mechanical Engineers, pp. 855–862. OMAE2010-21109. https://doi.org/10.1115/OMAE2010-21109. Zhan, D., Molyneux, D., 2012. 3-Dimensional numerical simulation of ship motion in pack ice. In: ASME 2012 31st International Conference on Ocean, Offshore and Arctic Engineering, OMAE2012–83105. American Society of Mechanical Engineers, pp. 407–414. https://doi.org/10.1115/OMAE2012-83105. Zhang, Z., Qian, X., 2017. Effect of experimental sample size on local Weibull assessment of cleavage fracture for steel. Fatigue Fract. Eng. Mater. Struct. 40, 1128-11142. Zhou, L., Ling, H., Chen, L., 2019. Model tests of an icebreaking tanker in broken ice. Int J Nav Arch Ocean Eng 11, 422–434.

References Abdelnour, R., Comfort, G.M., Pilkington, R., Wright, B.L., 1987. Ice forces on offshore structures; model and full scale comparison and future improvements. Oceans 87, 24–29. https://doi.org/10.1109/OCEANS.1987.1160869. IEEE. Aksnes, V., 2011. A panel method for modelling level ice actions on moored ships. Part 1: local ice force formulation. Cold Rons Sci Tech 65 (2), 128–136. Babic, M., 1990. Discrete element simulations of river ice transport. In: Proc. IAHR Ice Symp, pp. 564–574. Bird, K.J., Charpentier, R.R., Gautier, D.L., Houseknecht, D.W., Klett, T.R., Pitman, J.K., et al., 2008. Circum-Arctic Resource Appraisal: Estimates Of Undiscovered Oil and Gas North of the Arctic Circle (No. 2008-3049). Geological Survey (US). Cole, T., 2005. Model Tests on “Kulluk” in Ice Conditions. NRC Report SR-2005-27 IOT Project No. 2019. National Research Council Canada. https://doi.org/10.4224/ 8895895. Comfort, G., Singh, S., Spencer, D., 1999. Evaluation of Ice Model Test Data for Moored Structures. National Research Council Canada. NRC Publications Archive, PERD/ CHC Report 26-195. https://doi.org/10.4224/12328046. Cundall, P.A., Strack, O.D., 1979. A discrete numerical model for granular assemblies. Geotechnique 29 (1), 47–65. Di, S., Xue, Y., Wang, Q., Bai, X., 2017. Discrete element simulation of ice loads on narrow conical structures. Ocean. Eng. 146, 282–297. Du, Q., Faber, V., Gunzburger, M., 1999. Centroidal Voronoi tessellations: applications and algorithms. SIAM Rev. 41 (4), 637–676. Eik, K., Marchenko, A., 2010. Model tests of iceberg towing. Cold Rons Sci Tech 61 (1), 13–28. Feng, L., Qian, X., 2018. Size effect and life estimation for welded plate joints under low cycle actions at room and low ambient temperatures. Thin-Walled Struct. 132, 195–207. Ghoshal, R., Yenduri, A., Ahmed, A., Qian, X., Jaiman, R.K., 2018. Coupled nonlinear instability of cable subjected to combined hydrodynamic and ice loads. Ocean. Eng. 148, 486–499. Hamilton, J.M., 2011. The challenges of deep water Arctic development. In: The TwentyFirst International Offshore and Polar Engineering Conference. International Society of Offshore and Polar Engineers. Hamilton, J., Holub, C., Blunt, J., Mitchell, D., Kokkinis, T., 2011. Ice management for support of arctic floating operations. In: OTC Arctic Technology Conference. Offshore Technology Conference, Houston, TX, OTC-22105-MS. https://doi. org/10.4043/22105-MS. Hansen, E.H., Løset, S., 1999. Modelling floating offshore units moored in broken ice: comparing simulations with ice tank tests. Cold Rons Sci Tech 29 (2), 107–119. Hansen, E.H., Løset, S., 1999. Modelling floating offshore units moored in broken ice: model description. Cold Rons Sci Tech 29 (2), 97–106. Hopkins, M.A., 2004. Discrete element modeling with dilated particles. Eng. Comput. 21 (2/3/4), 422–430. Hopkins, M.A., Hibler, W., 1991. Numerical simulations of a compact convergent system of ice floes. Ann. Glaciol. 15, 26–30. Hopkins, M.A., Tuhkuri, J., 1999. Compression of floating ice fields. J Geophys ResOceans 104 (C7), 15815–15825. Hopkins, M.A., Tuthill, A.M., 2002. Ice boom simulations and experiments. J Cold Rons Eng 16 (3), 138–155. ISO, I, 2010. 19906: Petroleum and Natural Gas Industries–Arctic Offshore Structures. ISO, Geneva. Jenssen, N.A., Hals, T., Jochmann, P., Haase, A., Dal Santo, X., Kerkeni, S., et al., 2012. DYPIC–A multi-national R&D project on DP technology in ice. In: Dynamic Positioning Conference. Marine Technology Society, Houston, TX. Ji, S., Li, Z., Li, C., Shang, J., 2013. Discrete element modeling of ice loads on ship hulls in broken ice fields. Acta Oceanol. Sin. 32 (11), 50–58. Karulin, E.B., Karulina, M.M., 2011. Numerical and physical simulations of moored tanker behaviour. Ships Offshore Struct. 6 (3), 179–184. Karulin, E.B., Karulina, M.M., 2013. Determination of loads on mooring system during the semisubmersible interaction with ice. In: The Twenty-Third International Offshore and Polar Engineering Conference. International Society of Offshore and Polar Engineers, Anchorage, Alaska, 30 June-5 July, 2013. Karulin, E., Karulina, M., 2014. Peculiarities of multi-legged platform operation in ice condition. In: ASME 2014 33rd International Conference on Ocean, Offshore and Arctic Engineering, OMAE2014–23203, V010T07A012. American Society of Mechanical Engineers, San Francisco, CA. June 8-13, 2014. https://doi.org/10. 1115/OMAE2014-23203. Karulin, E., Karulina, M., Toropov, E., Yemelyanov, D., 2012. Influence of ice parameters on managed ice interaction with multi-legged structure. In: Proc. Of 21st IAHR International Symposium on Ice. “Ice Research for a Sustainable Environment”, Dalian, China. Keinonen, A., Browne, R., Revill, C., 1996. Icebreaker Characteristics Synthesis. Report for Transportation Development Centre, Transport Canada. TP 12812E. Kim, M.H., Tahar, A., Kim, Y.B., 2001. Variability of TLP motion analysis against various design methodologies/parameters. In: The Eleventh International Offshore and Polar Engineering Conference. International Society of Offshore and Polar Engineers. Kjerstad, Ø.K., Metrikin, I., Løset, S., Skjetne, R., 2015. Experimental and phenomenological investigation of dynamic positioning in managed ice. Cold Rons Sci Tech 111, 67–79.

17