CC BY-NC-ND
0ELSEVIER
Contents lists available at ScienceDirect
Cold Regions Science and Technology
journal homepage: www.elsevier.com/locate/coldregions
Numerical investigations of ship-ice interaction and maneuvering performance in level ice
Quan Zhou, Heather Peng *, Wei Qiu
Faculty of Engineering and Applied Science, Memorial University of Newfoundland, St John's, Canada
ARTICLE INFO ABSTRACT
This paper presents a solution to the ship-ice interaction problem solved in the time domain by a combined method involving numerical simulations and semi-empirical formula. The breaking process of an intact level ice by an advancing ice breaker has been investigated by a 2D numerical method. An interaction detection technique has been introduced. In order to compute the ice load, a pressure-area relationship for the ice was applied, and a flexural ice plate model was developed. A 3-degree-of-freedom (3DOF) model was adopted to simulate ship maneuvering in level ice. The proposed numerical method has been validated by using the experimental results of an 1:20 scaled CCG R-Class icebreaker model and the full-scale turning circle results.
© 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/4.0/).
CrossMark
Article history: Received 24 July 2014
Received in revised form 29 September 2015 Accepted 27 October 2015 Available online 14 November 2015
Keywords: Ship-ice interaction Ice breaking Contact detection Flexural ice plate Ice loads in 3DOF PMM tests
1. Introduction
The growing interests in oil and gas exploration in Arctic and subArctic regions and in the transportation through the Northern Sea Route have led to extensive research on better understanding of ship-ice interaction and vessels' maneuvering performance in ice covered water. For a ship advancing and turning in level ice, the assessment of ship maneuvering performance needs to be addressed. For a ship advancing in ice, the primary interest is on the prediction of the ice resistance. For maneuvering studies, the transverse force and turning moment are as important as the resistance when the ship is in turning.
It is challenging to predict the resistance that a ship is encountered in the intact ice field. One of the common assumptions, the principle of superposition to the total resistance is widely accepted by many researchers when the global ice load models were developed (for example, Enkvist, 1972; Kashteljan et al., 1969; Keinonen, 1996; Lewis and Edwards, 1970; Lindqvist, 1989; Milano, 1972; Riska et al., 1997; Spencer, 1992). In the past decades, efforts have been made to improve these models and implement them in numerical methods. For instance, Wang (2001) adopted the conceptual framework of the nested hierarchy of discrete events, proposed by Daley (1991, 1992) and Daley et al. (1998), and simplified the method by considering three continuum processes of crushing, bending, and rubble formation. A geometric grid method, which requires the discretization of the entire
* Corresponding author. E-mail address: hpeng@mun.ca (H. Peng).
computational domain, was proposed to simulate continuous contacts between the structure and the level ice. The mechanics of ice crushing-bending failure was applied and the ice loads were numerically computed. This approach has been adopted and modified by many other researchers, such as Nguyen et al. (2009), Sawamura et al. (2009), Su et al. (2010) and Zhou and Peng (2014). Tan et al. (2013) proposed a 6 degree-of-freedom (DOF) ship-ice interaction model by discretizing a belt area around the ship waterline. Their work provided a possibility to include ship heave, roll and pitch motions during the icebreaking process. A further study on the effect of dynamic bending of level ice during ship-ice interaction was carried out with this 6 DOF model by Tan et al. (2014).
Although the fundamental mechanism of icebreaking is not fully understood, many semi-empirical methods were proposed to estimate the resistance of a vessel in ice. However, limited effort has been made on the estimation of the yaw moment for a ship maneuvering in level ice. Lau et al. (2004) proposed a method to estimate the total yaw moment that is analogous to that of ice resistance in the work of Spencer (1992). The total yaw moment was divided into hydrodynamic, breaking, submergence, and ice clearing components. The formulas for the terms associated with breaking and submergence were presented. In their work, the ice-induced forces were considered as three concentrated loads among which two were acting at the bow and the other was on the midship section. The yaw moment was obtained by multiplying those forces and the corresponding lever arm lengths. Martio (2007) developed a numerical program to simulate the vessel's maneuvering performance in uniform level ice based on Lindqvist's ice resistance model (Lindqvist, 1989). The major contribution was to extend the
http ://dx.doi.org/10.1016/j.coldregions.2015.10.015
0165-232X/© 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
analytical formulas to 3 DOF. The effects of bending and submergence forces on sway and yaw motions were taken into account. A relationship between the resistance and the yaw moment, as well as the transverse force, was developed. Only crushing was calculated numerically, and the other two terms were obtained by using analytical formulas. Nguyen et al. (2009) and Su et al. (2010) developed 3DOF ship-ice interaction models. In their work, it was assumed that only the term associated with breaking had effect on sway and yaw motions. The remaining components of ice loads were considered as resistance that only affected the surge motion.
When a ship is maneuvred in level ice, the ship-ice interaction can be described as a continuous process. As the ship moves into the ice sheet, local shear fracture, crushing failure, and vertical deflection occur at the edge. These ice failure modes manifest themselves as a breaking force on the hull. When the stress somewhere inside the ice sheet exceeds its limit, a bending failure will occur which leads to the breaking of ice floes. The broken ice pieces are further accelerated and rotated by the ship until they are parallel to the hull. Subsequently, the broken pieces slide along the hull until they lose the contact with the ship. During this process, two force components are assumed: the buoyancy caused by the density difference between the ice and the water, and the clearing force that is attributed to the accelerating and rotating ice floes as well as the frictional force during ice sliding. The total ice resistance therefore consists of these three components discussed above. In this paper, the breaking process was simulated by a 2D numerical method involving the discretization of ice edges and ship's waterline. The contact loads were calculated, compared with the bending capacity of the ice sheet, and were used to determine the possibility of ice breaking. Equations of ship motions were then solved, and the ice edges were updated. The clearing and buoyancy forces were calculated using published empirical formulas. The effects on sway and yaw motions were taken into account in this paper. A flexural ice plate model was also developed to modify ice loads. Experimental results of the PMM tests of an 1:20 scale R-Class icebreaker ship model were used to validate the numerical method.
2. Modeling of ship-ice interaction
2.1. Coordinate systems and mathematical formulations
As shown in Fig. 1, three coordinate systems are employed in the computations. The earth-fixed coordinate system is denoted as Oe-XeY,Ze. The position and heading of a ship is described in this coordinate system. The ship-fixed coordinate system, o-xyz, is with ox-axis pointing to the bow, oy-axis pointing starboard and oz-axis positive downward. The origin, o, is located at the interaction point of the longitudinal central plane, the waterplane and the midship section. Equations of ship motions are solved in the ship-fixed coordinate system; The ice-fixed coordinate system, oi-xiyizi, is on the ice plate and parallel to Oe-XeYeZe. The coordinates of the discretized points on ice edges are defined in oi-xiyizi.
The equations of motion for maneuvering, including surge, sway and yaw, are given below:
x = cos(^)u— sin(^)v
y = sin(^)u + cos(^)v (1)
miu — mvr — mxGr2 = XH + XP + XR + Xice mV + mxcr + mur = YH + Yr + Yke (2)
¡z-r + mxG V + mxGur = NH + Nr + Nice
where m and ¡z are the mass of the ship and the mass moment of inertia, respectively; xG is the center of gravity defined in the ship-fixed coordinate system; X, Y and N denote the forces and moment with respect to ox, oy and oz-axis, respectively; x andy are the coordinates; tp is the orientation ; u, v and r denote the corresponding velocities in the ship-fixed coordinate system; the overhead dot represents the time derivative; variables with the subscripts, H, P, R and ice, represent forces and moment due to hydrodynamic force, propeller force, rudder force and ice load, respectively.
22. ¡ce-induced forces
The following assumptions were made in order to simplify the hull-ice interaction problem:
• The intact ice sheet is semi-infinite with uniform thickness. It is stationary with respect to Oe-XeYeZe;
• The hull-ice interaction is a continuous process which involves repeating cycles of contacting, crushing and bending and breaking;
• Superposition principle is applied to the total ice induced resistance;
• Vertical displacement of the ship is neglected;
• Contacting surfaces remain flat during crushing;
• The shape of the broken ice floe is considered circular for computing in the bending failure.
The total resistance due to ice can be written as:
Xice = Xbr + Xcl + Xbuoy
Yice = Ybr + Ycl + Ybuoy (3)
Nice = Nbr + Nc, + Nbuoy
The breaking force, denoted as Fbr = [Xbr, Ybr, Nbr]T, exists when the ice plate is crushed and bent downward until the breaking occurs. The dynamic clearing force, Fcl, acts on the hull as the broken ice floes rotates, accelerates, submerges and slides along the ship. The static force due to the buoyancy of ice is denoted as Fbuoy. The open-water resistance is considered as hydrodynamic force and is not part of ice loads. Each term in the equation above has effect on surge, sway and yaw motions. The breaking term is calculated numerically based on the detected contacts at each time step. The resistances due to clearing and buoyancy (Xcl and Xbuoy) are calculated using published empirical formulas. The
Fig. 1. Three coordinate systems.
ice edge at 100s
at 101s
ice edge at 102s
Fig. 2. Discretization of ship waterline and ice edges.
relationship between X and N in the clearing and buoyancy terms is derived in the following sections.
2.2.1. Ice breaking force
In the work of Zhou and Peng (2014), a numerical method was proposed to study ship maneuvering problems in level ice. A 2D interaction detection technique was introduced, and the calculation of ice breaking force was provided with details. This study is a continuation of the previous work and the same detection technique was applied. A brief summary of the method and its numerical implementation are presented below.
The waterline of the ship is discretized as a closed polygon while an ice edge is discretized as broken lines (Fig. 2). Note that this method is different from those proposed by Nguyen et al. (2009), Sawamura et al. (2009) and Wang (2001). Instead of discretizing the entire ice sheet, the present method only focuses on the ice edge where the ship-ice interaction may occur, which leads to improved computational efficiency.
Four contact detection scenarios as shown in Fig. 3 were considered, including (1) no ice node is inside the ship and no ship node is inside the ice, (2) ice nodes are inside the ship and ship nodes are inside the ice, (3) ice nodes are inside the ship but no ship node is inside the ice, (4) no ice node is inside the ship but ship nodes are inside the ice. They were distinguished by the Point-in-Polygon method which determines the relationship between a point and a closed polygon. For example, if an arbitrary ice node is inside the
polygon of the ship waterline, Scenario 2 or 3 will occur. Then all the nodes on the ship waterline need to be checked. If all ship nodes are inside the ice sheet, Scenario 2 needs to be considered; otherwise, Scenario 3 occurs. Scenario 4 could be identified similarly if no ice node is inside the ship waterline and the ship nodes are inside the ice sheet. As the high fidelity and the computational efficiency need to be satisfied, it is crucial to maintain the accuracy for Scenario 4 while using relatively coarse mesh to discretize the ice edge. Once a contact scenario was found, interpolations were applied to identify the intersections and apexes of the ideal ice wedge, which is illustrated by the blue dash lines in Fig. 3. The open angle of the ice wedge is denoted as 9. The breaking force was calculated based on the interaction area and the ideal ice wedge. If the force exceeds the bending limit of the ice sheet, a circumferential crack will occur and a floe will be broken off the ice sheet. The newly formed ice edge will be discretized into nodes and used to replace the nodes that are enclosed by the arc. This process is known as the ice edge update. The flow chart for the numerical implementation of the algorithm is given in Fig. 4.
As the hull contacts with ice edge, ice crushing starts at the contact point. It continues until the ice plate is broken by bending moment. More than one contact zone could exist at the same time. In each contact zone there will be a crushing force acting on the hull. To determine the magnitude of the crushing force in each contact zone, the following equation is adopted:
Fcr = Pa
where Fcr is the crushing force, which is normal to the contact surface and pointing inward to the ship; pave is the average pressure on the contact surface; Acr is the area of the contact surface that is obtained numerically at each time step, which is calculated by
Acr —
2 sin(/')
Lc-hi= tan(/')
Lh + Lh---
cos / 0
ifLc < hi ■ tan (/') ifLc > hi ■ tan(/')
where Lh is the characteristic width of the overlap area; Lc is the indentation; h{ is the ice thickness; Lh and Lc are numerically determined at
Fig. 3. Contact detection scenarios.
Fig. 4. Flow chart.
each contact point when the overlap area has been identified; and ¡3' is the normal frame angle of the ship hull that is determined by
tan(ß') = tan(ß) ■ cos(a) tan(ß) = tan(a)/tan(y)
where a and y are the waterline angle and the buttock angle, respectively, as shown in Fig. 1; and ß' is the normal section angle.
In the previous work, for example, Lau et al. (2004), Martio (2007), Nguyen et al. (2009), Sawamura et al. (2009), Su et al. (2010) and Zhou and Peng (2014), a constant contact pressure (crushing strength, ocr) was used when calculating the crushing force. However, the analysis of full-scale data shows that the contact pressure varies as the nominal contact area changes. This phenomenon can be expressed by the pressure-area relationship (also known as the P-A curve) and a form of power relation as below is widely accepted:
Pave — p0Acr
where p0 is assumed equal to the crushing strength and here ex is chosen as —0.5.
Based on the work of Zhou and Peng (2014), the three components of ice breaking force are calculated by:
Xbr — ( sin(ß') tan(a) + ^ 1 + tan2(a) cos2(ß')j ■ Fa
Ybr —
Zbr —
sin( ß') —n tan(a)-
cos(a)— cos2( ß')
' cos2 (a) + sin2 (a) cos2 (ß')
, ,) sin(a) sin(ß') cos(ß') 1
cos(ß') —^ ffi ffi ffi ffiffi ■ Fcr
y cos2 (a) + sin2 (a) cos2(ß') J
■ Fcr (9)
where Xbr and Ybr are the force components along x-axis and y-axis, respectively; Zbr is the vertical component that is further compared to the bearing capacity of the ice plate to determine if the breaking process would occur.
The yaw moment is obtained from:
Nbr = -Xbr ■ y + Yb,
where x and y are the coordinates of the contact location with respect to o-xyz.
It should be pointed out that Eqs. (8) to (11) give the breaking force at one contact location. The global breaking force can be calculated by summing the forces in different locations.
22.2. Crushing force adjustment
Valanto (1989) reported that a rapid flexural failure was observed in the experiment when the flexural strength ratio was large (E/cf = 6400), while considerably more time was required for flexural failure in the low ratio case (E/cf = 1400). It was also pointed out that the difference in two cases was important since it would significantly affect the average ice resistance. In the work of Nguyen et al. (2009), Su et al. (2010) and Wang (2001), rigid ice plate models were assumed and the vertical deflection by bending was not considered. Their models were thus unable to simulate what was observed in the experiment by Valanto (1989). In this study, the ice plate flexural deflection was taken into consideration, and an adjustment factor of contact surface was derived.
The actual ship-ice interaction is a dynamic process as both the ship and the ice plate have vertical motions. However, the process was simplified in this work by assuming that the ship does not move vertically. As the ship moves forward, the elastic ice plate will be bent. Since the ice deflection is small compared to the characteristic length ofthe ice, a parallel downward movement of the ice plate was further assumed and there is no bending involved (see Fig. 6). The contact force, Fcr, will lead to ice crushing and bending simultaneously. The water support was considered as an elastic foundation.
As shown in Fig. 6, from the initial contact to the present time step, the ship penetrated into ice by 8s, which was determined numerically, causing a vertical elastic deflection, Se, and a vertical crushing height 8c, which satisfy the following relationship:
Se + Sc = Ss / tan(/')
Both Se and Sc are unknown. Both are induced by the same force, the vertical component of the contact force, Fcr.
The relation between Se and Fcr is derived by considering the response of a homogeneous and isotropic elastic ice plate with an open angle of 9, which rests on a liquid and is subjected to a static vertical concentrated load P on the apex. Kerr (1976) proposed the
/////////////////// Fig. 6. Flexural deflection of ice plate.
following formula to calculate the deflection at the apex of a semiinfinite plate with an arbitrary open angle:
wo =2(0)
where P is the concentrated load; y = pwg is the weight per unit volume of the liquid; D = Eh3/(12(1 — v2)) is the flexural rigidity of ice plate; E is the ice Young's modulus; v is Poisson's ratio. Substituting P = Fcrv and w0 = 8e into Eq. (13), the following relationship can be achieved
Fcrv = ~nrVPwgD ■ Se
The relationship between Fcr,v and 8c is implicitly expressed by Eqs. (4), (5), and (10).
The binary search method was applied to obtain 8c. Once it is obtained, the contact area Acr in Eq. (4) can be modified by multiplying (8c/8v)2, because the flexural deflection of ice plate will reduce the contact area. Therefore, the adjusted crushing force is:
Fcr = Pave ■ Acr ■ (
The effect of the flexural deflection of the ice sheet can be considered by replacing Fcr with F£ in Eqs. (8) to (10).
2.2.3. Bending failure criterion
The vertical force component increases as the ship penetrates into the ice plate. When it exceeds the bearing capacity of the ice sheet, the bending failure would occur and a circular ice floe would be broken
3D view
side view
case 1
case 2
Fig. 5. Ideal ice wedge and contact surface.
offthe plate. The bearing capacity is calculated by Kashtelyan's method (Kerr, 1976).
Pflex = C/l-1 afh
where 9 is the open angle of the ideal wedge; of is the flexural strength of ice plate; h, is the ice thickness; and Cf is an empirical coefficient. It should be noted that Eq. (16) is derived from a quasi-static condition. It does not consider the dynamic effect due to varying ship velocities.
Kashtelyan suggested a small value (around 1.0) for Cf, while Nguyen et al. (2009) used a value of 4.5 with no explanation. In the work of Su et al. (2010), a study of Cf was carried out and a value of 3.1 was finally adopted. The value of Cf is adopted as 2.2 in this work and detailed discussions can be found in Section 3.2.
After the vertical force, Zbr, and the bearing capacity of determined, the next step is to investigate whether the ice wedge will be broken from the ice sheet. If the vertical force exceeds the capacity, i.e., Zbr > Pflex, the bending failure will occur and a circular ice floe will be broken off the ice plate; otherwise the ice edge will only be crushed. The radius of the broken ice floe is determined by using the formula in the work of Wang (2001):
R = Ci ■ l ■ (1 + Cvvn,2)
where R is the radius of circumferential crack; vn2 is the velocity component normal to the contact surface; Cl is the parameter that defines the broken ice size to its characteristic length; Cv is a tunable coefficients, and both are discussed in Section 3.2; and l is the characteristic length of ice plate given by:
( Eh3 )
\12(1-V2) pwg)
where E is Young's modulus of ice; v is Poisson's ratio; and pw is the density of sea water.
2.2.4. Buoyancy and clearing forces
The ice resistances due to buoyancy and ice clearing which includes ice rotating and sliding are obtained using the regression model proposed by Spencer and Jones (2001):
Xcl = 2.03 F-0971 piBhiu2 Xb = 2.67 A pghiBT
When a vessel starts turning, the forces on the port and starboard sides are not symmetrical anymore. A net transverse bow force exists. The net force results from the yaw rate of the vessel, i.e., the vessel will encounter larger transverse force if it turns rapidly.
u Sinffl!) Y LAI i,bow --i,bow
tan(a)
01 = atan^r (Lwl —Lbow), u)
where Lwl is the ship length at the design waterline, and Lbow is the length of the bow as shown in Fig. 7.
The force acting on the parallel midship section is determined by assuming the load is in the opposite direction of the velocity:
V -r Y _ 1 imid — Lmid ■ ^i ' ,,
where u and v are surge and sway velocities, respectively; and Cmid is equal to 2.0.
After the transverse components of buoyancy and clearing forces, AYitbow and Yifmid, are obtained, the yaw moment can be obtained by multiplying the transverse force by the corresponding lever arm:
Ni = AYi>bo„ ■ xb + Yimid ■ Xm where xm is calculated by
1 2 1 2
3 L?(L1r + v)-3 L2(L2r + v)
2 L1 (L1 r + v)+ 2 L2(L2r + v)
L1 = 2 LWL Lstern + x0 L2 = 1LWL Lbow x0
and Lstern is the length of the stern and x0 = — v/r.
It is assumed that AYiibow acts on the middle of the bow, i.e.,
where Fh = is the Froude number of the ice; u is the advancing
speed; B and Tare the ship breadth and draft, respectively; p, is the density of ice; Ap is the density difference between ice and sea water.
To calculate the y-component of buoyancy and clearing forces, it is assumed that the ice forces are acting on the bow part of the ship and the midship part. It is also assumed the ship is symmetric about the longitudinal centreplane and the ice force acts on the ship symmetrically when the ship moves straight. It is further considered that the bow is formed of two flat planes. Therefore, the directions of buoyancy and clearing forces should be the same as that of the breaking force, i.e.:
cos(a)- cos2( 3')
sin(/3') —¡1 tan(a
Xi Ybr \/cos2(«
v — —i — Yibow = 2 ■ Xbr =
2 (a) cos2 (3') Xi
sin (3') tan(a) + ¡J 1 + tan2 (a) cos2(3') 2
where the subscript i is either "cl" or "b", representing the clearing or buoyancy force, respectively. Note that the force is divided by 2 as the force is assumed to act on either side of the bow.
xb = 2 LWL — 22 Lbow
At each time step, the contact will be first detected by the numerical method. If a contact is found, the pressure-area relationship will be applied to determine the pressure. After that, the crushing force is calculated by Eq. (4). The breaking force can be obtained by using Eqs. (8) to (11) and further modified by considering the effect of flexural deflection of the ice plate. Finally, the buoyancy and clearing forces are determined by using Eqs. (19) to (20), (22), (24) and (25).
Ymidî
Fig. 7. Arms of rotation.
2.3. Hydrodynamic forces
The hydrodynamic forces are determined as follows according to the work of Hirano (1981):
Xh = -Xu u +X(u) + (my + Xvr)vr Yh = - Yv v-Yr r -mxur + Yho Nh = -Nr r -Nv v +Nho + Yho ■ xc where
YHO = YvV + Yrr + Yv|v|V|v| + Yr|r|r|r| + Yvrvr
Nho = Nvv + Nrr + Nvvrv2r + Nvrrvr2 + Nr|r|r|r|
the constant coefficients, referred to hydrodynamic derivatives, are
Xu= dù ' Yv = dv ' Nvvr = Bv2dr'
X(u) represents the resistance due to water and is calculated by using the prediction method proposed by Holtrop and Mennen (1982). In this paper, the hydrodynamic derivatives were obtained from the PMM tests with a 1:20 scale ship model.
2.4. Operational forces
The operational forces consist of two parts: the thrust from the main propeller and the lift and drag forces due to the rudder deflection. The model used by Gong (1993) was adopted in this work.
Xp = (1 -tp )p„n2D2pI<T (Jp ) Xr = -(1—r)Fn sin(S) Yr = -(1 + aH ) Fn cos(S) Nr = -(1 + aH )xrFn cos(S)
where tP and tR are the wake fractions due to propeller and rudder, respectively; n is the propeller revolution per second, (rps); DP is the propeller diameter; KT is the thrust coefficient; JP = u/(nDP) is the propeller advance ratio; aH is the coefficient of hydrodynamic force on ship hull induced by the rudder action that can be estimated from model test results; xR is the longitudinal location of the rudder; 8 is the rudder deflection angle; and FN is the rudder normal force determined by:
Fn = ~ P
2' wA + 2.25
ArUr sin(aR)
4th-order Runge-Kutta method according to the following vector form rewritten from Eq. (2):
Vfc+1 = Vk + £ (k + 2k2 + 2k3 + k4)
k1 = f (Vk, flc)dt k2 = f (vk + kr, fk )dt
k3 = f(Vk + f, fkjdt
= f (Vk + k3, fk)dt
Vk = f (Vk, fk)
= A-1 fk-A-1N(Vk) Vk = [ùk, vk, rk\T
XH + XP + XR + Xic
N(Vk)=
Y H + Yr + Yic Nh + Nr + NiC
m mxG mxc iz
mvkrk-mxcr2 mùkrk mxcùkrk
where the subscript k and k +1 represent the current and the subsequent time steps, respectively. The position and the orientation of the ship are calculated by using Eq. (38) along with Eq. (1). Coordinates of the nodes on the ice edge are then updated based on the new position and orientation.
3. Simulation results
The R-Class icebreaker, CCGS Sir John Franklin, was used in the present studies. It has a relatively simple underwater hull form. The vessel is one of most referred icebreaking hull forms in the literature since it was built in 1979. Both full-scale sea trials and ship model tests at various scales of the ship have been carried out for this icebreaker (Liu, 2009). To validate the proposed method, a series of Planar Motion Mechanism (PMM) test results were used. The numerical predictions were compared with the results of full-scale sea trials involving straight motion and turning.
The principal dimensions of the ship and the ice properties at model scale are listed in Tables 1 and 2, respectively. The waterline profile is presented in Fig. 8.
3.1. Convergence studies
where A is the aspect ratio of the rudder; AR is the projected rudder area; UR is the effective rudder inflow velocity; aR is the effective rudder inflow angle. Different models are available for the calculation of the effective rudder inflow velocity and angle (Gong, 1993; Hirano, 1981; Kobayashi et al., 2003). The one proposed by Hirano (1981) was employed in this paper.
2.5. Numerical implementation
At each time step, the external forces on the right hand side of Eq. (2) are obtained by the numerical methods described above. Eq. (2) is then solved to obtain the acceleration in the o-xyz frame. The velocity and the position of the ship are integrated by using the
Prior to simulating the PMM tests, convergence studies were carried out to investigate the effect of the discrete length (the distance between two adjacent points) and the time step on the mean ice load. The
Fig. 8. R-Class icebreaker ship model (model scale, unit: m).
Table 1
Principal dimensions (model scale).
Parameter Notation Value Unit
Scale K 20 -
Length lwl 4.65 m
Breadth i 0.96 m
Draft T 0.35 m
Block Coefficient Ci 0.611 -
Table 2
Ice mechanical properties (model scale).
Parameter
Notation
Flexural strength Crushing strength Ice thickness Young's modulus Poisson ratio Friction coefficient Ice density
oc 100
28,140
mm kPa
Model Scale 0.5m/s —*—
Full Scale 2.236mfe —»— Full Scale 3 130m/s _a_ ■
—S—O
0.1 1 10 100 Non-dimensional time
Fig. 10. Convergence of ice resistance with respect to the time step.
applicability of the present method to ships at different scales and speeds were also studied. The 1:20 scale R-Class icebreaker model and the full-scale ship were used. In the simulation of the straight tow tests, the chosen velocities were 0.5 m/s and 0.7 m/s in model scale.
The nondimensional discrete length, l', and the nondimensional time step, dt', were defined as follows:
From Figs. 9 and 10, it can be seen that the computations are converged in all the cases when the nondimensional discrete length is smaller than 0.01 and the nondimensional time step is greater than 5. This indicates that the scale and the ship speed will not affect the accuracy of simulation if proper values of nondimensional length and time step are selected, i.e., the proposed numerical method is suitable for various scales and ship speeds. This conclusion needs to be further verified by using experimental and/or numerical results in different scales. In the following simulations, the nondimensional discrete length and time step were chosen as 0.005 and 5, respectively.
3.2. Study of empirical parameters
where l is the discrete length and dt is the time step used in the simulations.
Fig. 9 presents the mean ice resistances computed with various discrete lengths, and Fig. 10 shows the resistances computed using different time steps. Good convergence in model scale is observed. The corresponding speeds for full scale are 2.236 m/s (4 knots) and 3.13 m/s (6 knots). Since the full-scale measurements are not available at these speeds, the simulated results were compared to the empirical solution based on the work of Spencer and Jones (2001). The resistance in full scale was estimated to be 389.6 kN and 442.5 kN, respectively. Good agreement is also found for full scale cases.
The parameters, Cf, Cl and Cv, need to be determined for different types of ice and vessels. Note that Cf is used to determine the peak force of a single contact and Cl and Cv are used to determine the size of a broken ice piece. The parameter, Q, significantly affects the time history of ice loads since smaller ice pieces will be broken off if a ship rapidly collides with the ice sheet.
The first step is to select Q according to physical observations. Wang (2001) used 0.32 in her study; Lau et al. (2004) used a different model but the crack to be 0.2 of the characteristic length. According to the observation from Liu et al. (2008), Q is selected to be 0.3 in this study. Once C is determined, the values of Cf and Cv are tuned by comparing the average resistance of simulation results against experimental measurements. Fig. 11 presents the mean resistance against ship velocity with different Cf values. Good agreement is observed between measurements
0.1 001
Non-dimensional length
Fig. 9. Convergence of ice resistance with respect to the discrete length.
Fig. 11. Resistance versus velocity with different C/.
I 87.45 m
1 80 <D
a. 73.68 60
v=0.5m/s —x— v=0.7m/s ---*■—-
0.3 0.4
CI value
Fig. 12. Resistance versus Cl with different velocities.
and simulations when Cf is 2.2. The parameter, Cv, is used to include the speed effect on the size of the broken ice floe. Wang (2001) used — 0.14 for Cv in her study; Lau et al. (2004) does not include the velocity effect when determining the size of the crack cusp. Cv is set as — 0.5 in model scale. It should be pointed out that Cv has unit of s/m, and Froude's law is applied to get the value for full scale (Cv = — 0.11s/m).
As shown in Fig. 11, the change of Cf will affect the relationship between the speed and the resistance (e.g., larger values of Cf result in a more steep increase trend of resistance). Once the value is determined,
3> 0.41 8
il I I ll
■v=0.3m/s I v=0.4m/s I lv=0.5m/s I lv=0.6m/s ■v=0.7m/s
Ratlo=1400 Ratio=6400 Rigid ice
2.85 2.9
Radius/thickness
0.3 0.35 0.4 0.45 0.5 0.55 0.6
Time, s
Fig. 15. Breaking force in a single crushing-breaking event (0.5 m/s).
it can be used for all the simulations. A larger Cf results in a larger peak value and hence a larger mean resistance.
Note that Cl denotes the proportion of broken ice radius and the characteristic length as shown in Eq. (17). With a smaller value of Q, smaller ice floe will be generated, which results in more frequent interaction between the vessel and the ice. This manifests into a larger mean ice resistance which is shown in Fig. 12.
Cv illustrates the size of ice floe varies according to normal velocity of collision. A higher velocity of the vessel results in a higher normal velocity of collision and hence a smaller ice floe. Fig. 13 presents the histogram of ratio of crack radius and ice thickness at varying velocities. The values of radius fall within a limited range. There is no obvious trend that radius distribution changes according to vessel velocities. Fig. 14 gives the relationship between mean resistance and velocity with different Cv values. Linear relationship is observed in all cases, but the slopes are different. Smaller value of Cv results in larger resistance. However, comparing to Cf and Q, the effect of Cv on resistance is limited.
All three parameters will affect the mean resistance. Cf determines the peak value of the breaking force and directly affects the mean resistance; Ci and Cv directly determine the size of ice floe and have indirect effects on the resistance. Further studies are recommended for those parameters if the measured peak value and ice floe sizes are available from field or model test data.
Fig. 13. Ratio of crack radius and ice thickness at varying velocities.
3.3. Effect of flexural ice plate model
A flexural ice plate model was introduced to solve the ship-level ice interaction problem. The effect and the significance of the model were
Table 3
Parameters for a single crushing-bending-breaking event.
Ice type fiv, 10-2m fie/fiv Duration, s fie/l Xbr, N
Rigid ice 1.295 0.00% 0.050 0.000 -7.609
High ratio 2.343 44.92% 0.090 0.032 -9.260
Low ratio 3.553 63.67% 0.135 0.069 -9.493
Fig. 14. Resistance versus velocities with different Cv
Fig. 16. Broken ice channel in a pure yaw model test (unit: m).
0 5 10 15 20 25 30 35 40 45 50
Time, s
Fig. 17. Time history of surge resistance.
studied in model scale. A single crushing-bending- breaking event was simulated, which started from the ship stem contacting with a semiinfinite ice plate until an ice floe was broken off the plate. Three different ice plates were used: a rigid ice plate, a flexural ice plate with a high flexural strength ratio (E/a/ = 6400), and an ice plate with a low flexural strength ratio (E/a/ = 1400). The strength ratios were determined based on the work of Valanto (1989). Fig. 15 presents the breaking force versus the simulation time. When the ship contacts with the ice plate, the breaking force increases from zero to a maximum value where the breaking failure occurs. After that, the ship loses contact with ice and the breaking force falls back to zero.
It is clear that the maximum values of breaking forces are identical (— 20.90N) in the three cases; however, the interaction durations and the mean forces vary during the interactions. Fig. 15
indicates that the lower the flexural strength ratio of the ice sheet is, the longer the duration of interaction will be. Table 3 presents a comparison between different types of ice sheets. The duration and the maximum vertical deformation increase greatly, which indicates bending is an component of this single event. However, a very small difference was observed in the mean breaking forces, especially for the ice floes with different flexural strength ratios.
The deflection ratios in the table represent the ratios between the flexural deflection and the characteristic length of the ice plate. Those values are smaller than 0.1, which verifies that the assumption of parallel downward movement is valid. The assumption that neglects the ship vertical movement can also be verified by comparing the maximum breaking force with the ship displacement.
Fig. 18. Time history of sway force.
0 5 10 15 20 25 30 35 40 45 50
Time, s
Fig. 19. Time history of yaw moment.
3.4. Model-scale PMM tests
PMM tests were simulated and compared with model test results to validate the numerical predictions. In the PMM tests, the amplitude of the lateral displacement was 2.5 m and the ice loads in surge, sway and yaw were measured. Average loads were recorded every 2 seconds. Each test lasted 50s. In order to compare the simulation results with the measured data, the identical test setup and data processing technique were applied in the simulation. A sample of ship model track and the broken ice channel is shown in Fig. 16.
Figs. 17,18 and 19 present the time histories of surge resistance, sway force, and yaw moment, respectively. The red lines represent the measured data and the blue scattered dots represent the simulation results. From Fig. 17, it can be observed that the measured surge resistance oscillates with a constant period (approximately 25 s) which is half of
the set value of the test. This is because the sway velocity was not strictly zero in the model tests. The average measured resistance is 73.680 N. A consistent average resistance (73.677 N) was obtained from the simulation. However, the explicit periodicity cannot be observed in the numerical results.
In Figs. 18 and 19, the consistent and explicit oscillatory motions can be seen in both the measured and numerical results. It is also clear that the amplitudes of computed sway force and yaw moment are slightly smaller than the measured data. This phenomenon is validated by the regression relationship as shown in Fig. 20. One possible factor that results in this phenomenon is the non-zero sway velocity. Due to the existence of sway velocity in the pure yaw test, the mean force acting on midship may be estimated incorrectly so that the predicted sway force and yaw moment are smaller than measured values. This should be further investigated with more experimental data.
Fig. 20. Relationship between yaw moment and yaw rate.
Sea trial data Simulation
Ship speed, m/s
Sea trial data Simulation Spencer and Jones
4 4.5 Ship speed, m/s
Fig. 21. Resistance in open water.
Fig. 24. Thrust in level ice covered water.
Ship speed, m/s
Fig. 22. Thrust in open water.
3.5. Full-scale straight running test
The full-scale results of free running tests were used to validate the total resistance and the propulsive performance of the ship. The full-scale free running tests were simulated in open water and in level ice.
Sea trial data Simulation Spencer and Jones
The computed propeller thrust and the total resistance were compared to the sea trial data from Keinonen (1996) and the analytical solutions using the formulas from Spencer and Jones (2001). The speed range of the tests in open water is from 6.34 m/s to 8.35 m/s, and that of the test in level ice is from 2.3 m/s to 6.3 m/s. The ice thicknesses vary from 0.489 m to 0.592 m. The simulations were conducted by restricting sway and yaw motions of the ship. Comparisons were made among the numerical results, empirical solutions and the sea-trial data.
Figs. 21 and 22 present the resistance and the propeller thrust in open water. It can be observed that the computed open water resistance and the thrust follow the same trend of the sea trial results, but are slightly smaller.
The total resistance and the thrust in ice are presented in Figs. 23 and 24. The numerical results by the present method are compared with the sea-trial data and the empirical solutions. The computed resistances are in good agreement with the empirical results, while the propeller thrusts based on numerical simulations and empirical formulas are not in good agreement. Neither the empirical results of the thrust nor the predicted thrusts show good agreement with the sea-trial measurement.
The relatively poor agreement between the numerical prediction and the measurement is not necessarily a weakness of this model since the thrust is difficult to measure accurately in the field trials. The fluctuated nature of the predicted ice resistance may be responsible for the discrepancy between the predicted and empirical thrusts. From an energy perspective, the computed ice load oscillates with a relatively small mean value and high peaks. This may result in more energy consumption in order to maintain the speed. Therefore, a larger mean thrust was obtained from the simulation.
3.6. Full-scale turning test
The turning performance was also investigated and compared with the results of sea trials in open water and in level ice (Keinonen, 1996). All trials were conducted with 35' starboard rudder. The comparison of the steady speed in turn and the turning diameter between sea trials and the simulation results is presented in Table 4. The turning
Table 4
Full-scale turning performance, Lwl = 92.31 m.
3 3.5 4 4.5 5 ! Ship speed, m/s
Fig. 23. Resistance in level ice covered water.
Condition Shaft resolution Speed in turn, m/s Turning diameter,
RPM Sea trial Simulation Sea trial Simulation
Open water 120 3.8 3.9 4.21 4.1
Level ice 170 5.3 5.14 13.12 12.80
1200 1000
g 800 c
5 600 400 200
0 200 400 600 800 1000 1200 1400 Transfer, m
Fig. 25. Turning track and ice channel.
circles and the ice breaking channel are presented in Fig. 25. The sample time histories of ship velocity and ice load in surge, sway, and yaw in the ice field are presented in Figs. 26 and 27.
In the simulations, the ship was kept straight for 600 m and then was turned using a full rudder angle (35 ). The initial speed was 6.17 m/s (12.0knots). The shaft speed was 170 rpm in ice and 120 rpm in water. Good agreement between the sea-trial and simulation results indicates that the numerical method is able to simulate the turning of the ship. The turning circles are shown in Fig. 25.
Fig. 26 presents the 3DOF ship velocities during the simulation period (700 sec). Fig. 27 shows the 3DOF ice loads in a time period of 40 s during turning in the level ice. The ship started turning at 90.48 s and obtained a stable speed (5.14 m/s) at 200 s. The steady sway speed is 0.58 m/s to the port side and the yaw rate is 0.5 deg/s clockwise. The
mean resistance during the steady turning is 0.566 MN. The mean value of yaw moment is 11.93 MNm. The velocities and the ice loads in sway and yaw each have non-zero value due to the turning. They oscillate around the mean value due to the ice breaking process.
4. Conclusions
A numerical model for simulating ship maneuvering performance in level ice has been developed and verified by comparing the numerical results to measured data and empirical solutions of the R-Class icebreaker, CCG Sir John Franklin. The numerical model has been developed considering similar strategies and assumptions made by previous researchers. It has been extended by considering the effect of buoyancy and clearing terms on sway and yaw motions. Moreover, the effect of the flexural deflection of ice plate is considered in the model. The water support of the ice plate is treated as an elastic foundation. The main conclusions can be drawn as follows:
1. A contact detecting method (the polygon-point algorithm) was proposed. As shown in Figs. 5 and 6, good convergence was obtained in a wide range of discrete lengths and time steps. By using nondimen-sional length and time step, the effects of scaling and velocity were avoided. This indicates that the identical values of nondimensional discrete length and time step could apply to different cases.
2. The effect of the elastic foundation for the ice plate on ship-ice interaction process was considered. The elastic foundation tends to increase the mean ice resistance and extends the duration of each single contact-bending-breaking process.
3. A method was developed to consider the effect of buoyancy and clearing forces on sway and yaw motions. It is based on an empirical calculation of resistance and takes the hull geometry and the ship speed into account. It is validated by using the model PMM tests and the full-scale sea trial results. Therefore, the proposed numerical method is suitable for both model-scale and full-scale simulations.
Fig. 26. Time history of velocities (full scale).
Fig. 27. Time history of ice loads and mean values (full scale).
Acknowledgements
The financial support from the NSERC CREATE Training Program for Offshore Technology Research (OTR) is highly appreciated. Mr. Don Spencer is gratefully acknowledged for his generous help on this research. The assistance by Dr. Claude Daley on the development of the flexural ice plate model is also acknowledged.
References
Daley, C., 1991. Ice Edge Contact-A Brittle Failure Process Model. D.Sc. thesis, Acta Polytechnica Scandinavica, Mechanical Engineering Series No. 100. Helsinki.
Daley, C., 1992. Ice edge contact and failure. Cold Reg. Sci. Technol. 21,1-23.
Daley, C., Tuhkuri, J., Riska, K., 1998. The role of discrete failures in local ice loads. Cold Reg. Sci. Technol. 27,197-211.
Enkvist, E., 1972. On the ice resistance encountered by ships operating in the continuous mode of icebreaking. Technical Report 24 The Swedish Academy of Engineering Science in Finland.
Gong, I.-Y., 1993. Development of maneuvering simulation program, part I: mathematical models for maneuvering motions in deep/shallow water with wind/current effects. Technical Report IR-1993-09. Institute of Marine Dynamics, NRC.
Hirano, M., 1981. A practical calculation method of ship maneuvering motion at initial design stage. Nav. Archit. Ocean Eng. 19, 68-80.
Holtrop, J., Mennen, G.G., 1982. An approximate power prediction method. Int. Shipbuild. Prog. 29,166-170.
Kashteljan, V.I., Poznjak I.I., Ryvlin, A.J., 1969. Ice Resistance To Motion Of A Ship. Marine Computer Application Corporation.
Keinonen, A., 1996. Icebreaker characteristics synthesis. Technical Report TP 12812E. AKAC Inc.
Kerr, A.D., 1976. The bearing capacity of floating ice plates subjected to static or quasi-static loads. J. Glaciol. 17, 229-268.
Kobayashi, H., Blok, J., Barr, R., Kim, Y., Nowicki, J., 2003. Specialist committee on Esso Osaka: final report and recommendations to the 23rd ITTC. Proceedings of the 23rd International Towing Tank Conference, pp. 8-14.
Lau, M., Liu, J.C., Derradji-Aouat, A., Williams, F.M., 2004. Preliminary results of ship maneuvering in ice experiments using a planar motion mechanism. Proceedings of the 17th International IAHR Symposium on Ice, pp. 479-487.
Lewis, J.W., Edwards Jr., R.Y., 1970. Methods For Predicting Icebreaking And Ice Resistance Characteristics Of Icebreakers. Transactions of Society of Naval Architects and Marine Engineers (SNAME), pp. 213-249.
Lindqvist, G., 1989. A straightforward method for calculation of ice resistance of ships. Proceedings of 10th International Conference on Port and Ocean Engineering under Arctic Conditions (POAC), pp. 722-735.
Liu, J., 2009. Mathematical Modeling Ice-Hull Interaction For Real Time Simulations Of Ship Manoeuvring In Level Ice (Ph.D. thesis) Memorial University of Newfoundland.
Liu, J., Lau, M., Williams, F., 2008. Numerical Implementation And Benchmark Of Ice-Hull Interaction Model For Ship Manoeuvring Simulations.
Martio, J., 2007. Numerical simulation of vessel's maneuvering performance in uniform ice. Technical Report Ship Laboratory. Helsinki University of Technology, Espoo, Finland.
Milano, V.R., 1972. Ship Resistance To Continuous Motion In Ice (Ph.D. thesis) Stevens Institute of Technology.
Nguyen, D.T., S0rb0, A.H., Sorensen, A.J., 2009. Modelling and control for dynamic positioned vessels in level ice. Proceedings of 8th Conference of Manoeuvring and Control of Marine Craft (MCMC'2009), pp. 16-18.
Riska, K., Wilhelmson, M., Englund, K., Leiviska, T., 1997. Performance of merchant vessels in ice in the Baltic. Technical Report 52. Helsinki University of Technology - Ship Laboratory.
Sawamura, J., Riska, K., Moan, T., 2009. Numerical simulation of breaking patterns in level ice at ship's bow. Proceedings of 19th International Offshore and Polar Engineering Conference (ISOPE), pp. 600-607.
Spencer, D., 1992. A standard method for the conduct and analysis of ice resistance model tests. Technical Report IR-1992-12. Institute for Ocean Technology (IOT), NRC.
Spencer, D., Jones, S.J., 2001. Model-scale/full-scale correlation in open water and ice for Canadian coast guard "r-class": icebreakers. J. Ship Res. 45,249-261.
Su, B., Riska, K., Moan, T., 2010. A numerical method for the prediction of ship performance in level ice. Cold Reg. Sci. Technol. 60,177-188.
Tan, X., Su, B., Riska, K., Moan, T., 2013. A six-degrees-of-freedom numerical model for level ice-ship interaction. Cold Reg. Sci. Technol. 92,1-16.
Tan, X., Riska, K., Moan, T., 2014. Effect of dynamic bending of level ice on ship's continuous-mode icebreaking. Cold Reg. Sci. Technol. 106,82-95.
Valanto, P., 1989. Experimental study of the icebreaking cycle in 2-d. Proceedings of the 8th International Conference on Ocean. Offshore and Arctic Engineering (OMAE).
Wang, S., 2001. A dynamic model for breaking pattern of level ice by conical structures. Acta Polytechnica Scandinavica Mechanical Engineering Me 156.
Zhou, Q., Peng, H., 2014. Numerical simulation of a dynamically controlled ship in level ice. Int. J. Offshore Polar Eng. 24,184-191.
