Numerical Research on Global Ice Loads of Maneuvering Captive Motion in Level Ice

: In level ice, the maneuvering motion of icebreakers has a major inﬂuence on the global ice loads of the hull. This study researched the inﬂuences of the drift angle and turning radius on the ice loads of the icebreaker Xue Long through a partial numerical method based on the linear superposition theory of ice loads. First, with reference to the Araon model tests performed by the Korea Research Institute of Ships and Ocean Engineering (KRISO), numerical simulations of Araon ’s direct motion were carried out at different speeds, and the average deviation between numerical results and model test results was about 13.8%. Meanwhile, the icebreaking process and modes were analyzed and discussed, compared with a model test and a full-scale ship trial. Next, the maneuvering captive motions of oblique and constant radius were simulated to study the characteristics of ice loads under different drift angles and turning radii. Compared with the maneuvering motion model tests in the ice tank of Tianjin University and the Institute for Ocean Technology of the National Research Council of Canada (NRC/IOT), the numerical results had good agreement with the model test results in terms of the variation trend of ice loads and ice–hull interaction, and the inﬂuences of drift angle and turning radius on ice resistance and transverse force, which have a certain reference value for sailing performance research and the design of the hull form of icebreaker ships, are discussed.


Introduction
The continuous shrinking of sea ice makes the summer navigation of Arctic routes a reality and confers a natural advantage in terms of shortening transportation distances, which has led to extensive research on icebreakers [1,2]. However, compared with open water, maneuvering motion is different, and much more difficult, in ice. The maneuvering motion of an icebreaker makes the ship-ice interaction more complicated, which has a severe impact on the ice loads. To summarize, it is challenging to calculate ice loads under maneuvering motion in level ice.
Early research on global ice loads mainly focused on ice resistance, and many theories of ice-ship interaction were developed in this period [3][4][5][6]. Among them, the linear superposition theory and empirical formulas of ice resistance were proposed; these have been widely accepted by many researchers and have provided a theoretical basis for the study of ice loads. Lindqvist [6] divided the ice resistance formula into an icebreaking resistance term, a submersion resistance term, and speed-related terms, and introduced stem-shaped parameters and sea ice properties. In addition, several analytical methods, including empirical, semi-empirical, and theoretical, were proposed to predict ship maneuvering in ice. The analytical method is a cost-effective manner with rapid estimating capability. Lindstrom [7,8] adopted a dynamic equation for an elastic plate on an elastic foundation to lating contact collision, deformation, damage, etc. Lau [37] used DECICE, a commercial discrete element code, to simulate captive motion in level ice, and the numerical results, ice forces, and yaw moments were in good agreement with the experimental data. For these simulations, a hybrid, simple, deformable, three-dimensional plate bending element of the Mindlin type (SDFE) was used, and the ice plate was divided into three zones and meshed based on experience. Meanwhile, a water foundation was resting under the level ice plate, and the buoyancy forces and moments were calculated by the wetted surface.
Zhan [38] simulated the turning circle and Zig-Zag maneuvers using DECICE2D combined with a ship maneuvering code in two-dimensional pack ice conditions, and the code was developed by the Japanese mathematical maneuvering group. For these simulations, mass damping was used to simulate the drag or viscous force applied to the element, and all the ice elements are unbreakable. Meanwhile, ice loads were joined in MMG as external forces, and the Runge-Kutta algorithm was implemented to solve governing equations of ship maneuvering and obtain the trajectories of ship maneuvering. The simulation results between cases with and without ice were compared and discussed. An expansion of this method was used to simulate ship motions with different drift angles, turning circles, and Zig-Zag maneuvers in three-dimensional pack ice, achieving good agreement [39]. According to the observation results during full-scale ship trials, Li and Kujula [40] proposed a novel ice failure model based on elliptical ice cusps. Then, they established a partial numerical framework: Lindqvist formulas were used to calculate submersion forces, and the finite difference method (FDM) simulated the performance of a ship in level ice [41].
Generally, the research of ship maneuvering in ice includes two aspects. First is the prediction of channel and ice loads under ship maneuvering motion; the second is to simulate the maneuvering motion based on a mathematical model. At present, using a numerical method to predict ship maneuverability in ice, the numerical method of the partial empirical formula is usually adopted [41]. Meanwhile, the research on ice-ship interactions is mainly about ice resistance, while there has been little research on ice loads under maneuvering motion. The characteristics of maneuvering ice loads in level ice are a matter of considerable importance, but the knowledge is still limited. In open water, captive model tests are a common method to establish a hydrodynamic model for predicted ship maneuverability [42]. Therefore, this paper researched the influences of drift angle and turning radius on ice loads by simulating the maneuvering captive motion in level ice. Firstly, the numerical method was verified by the mean ice resistance and icebreaking process of the model icebreaker Araon test. Then, the maneuvering captive motions of oblique and constant radius were simulated, and the influences of drift angle and turning radius on ice loads were discussed.

Calculation Method of Ice Loads
In this paper, the ice loads are divided into two types, icebreaking loads and submersion loads, based on the commonly used ice load linear superposition theory, which can be expressed as follows [18]: where F ice is the ice loads, F bre is the component of ice loads due to crushing and bending, and F sub is the component of ice loads due to the motion of broken ice, such as rotation, sliding, or accumulation. In order to simulate the level icebreaking process and reduce the influence of empirical factors in the ice load calculation process, LS-DYNA software was used to calculate the icebreaking loads. The calculation of submersion loads involves broken ice, water, and hull; the interactions between them are complex and have not been thoroughly studied. The Lindqvist empirical formula has been tested by many engineering practices and has a certain calculation accuracy. In the present ship performance model, the Lindqvist submersion resistance formula was extended to account for the calculation of submersion loads in maneuvering motion [18,19,24,41].

Calculation of Icebreaking Forces F bre
The icebreaking process of an icebreaker ship is a typical dynamic contact collision process. LS-DYNA is explicit, dynamic software mainly used to simulate and calculate the impact and collision of structures. The main points of the explicit solution method are as follows: The concentrated mass method is used to solve the dynamic response equation directly, which avoids the transformation of complex stiffness matrix and greatly improves the solving efficiency. The time is discretized by the solution format of central difference, and the acceleration is assumed to be constant within a time step. The midpoint velocity of the current time step is obtained by adding the midpoint velocity of the previous time step and the change of velocity. The displacement at the end of each time step can be obtained by adding the displacement at the beginning of each time step to the time integral of velocity. At the beginning of the time step, the acceleration is provided to meet the dynamic equilibrium conditions. Through the explicit solution, the node displacement and velocity can be further obtained.
The dynamic response equation of the structure after finite element discretization is expressed as follows: M ..
where M, C, K, x, and F are the mass matrix, damping matrix, stiffness matrix, displacement matrix, and external load matrix, which includes the interactions of ice-ice and ice-ship, respectively. Assuming that the acceleration remains constant in a time step and the explicit central difference method is used in the time dimension, the Formula (3) is solved in the form of: where, x(t n+1 ) is the displacement vector at t n+1 , .
x(t n−1/2 ) is the velocity vector at t n−1/2 . As shown in formula 1, the displacement of the node x(t n+1 ) can be obtained by the displacement at t n and the velocity at t n+1/2 . Thus, the nodes' displacement, velocity, and acceleration at all discrete time points can be obtained through the above recursive formulas. The above is the main solution process of the explicit integration method. The linear equations used in this method are non-coupled and can be directly calculated, which saves a lot of storage space and solution time for numerical calculation. The finite element software LS-DYNA used in this paper uses this method to solve the Lagrange problem. The calculation process of a single time step in the solution process is shown in Figure 1.
The time step is determined by the minimum grid of the numerical model: where ∆t i is the time step of the unit, n is the number of units, α is the scaling parameter, and the value of α is 0.9 to ensure the stability of the calculation. The time steps of numerical simulations are all on the magnitude of 10 −5 .
In order to ensure the transmission of the contact force between a hull and level ice, the choice of contact model is essential. When a hull interacts with level ice, it will continuously erode and destroy the level ice as the hull advances. This is a typical surface-to-surface contact type with erosion, so the *CONTACT_ERODING_ SURFACE_TO_SURFACE contact model was adopted. In this contact model, the hull is the master surface and the level ice is the slave surface. The determination of the master and slave surfaces mainly depends on the grid size and other factors. lution process. The contact between the ship and ice is judged and the interaction fo calculated. Meanwhile, the stress and strain of ice are calculated. If the pressure or p strain reaches the value of failure, the failure element will be deleted; otherwise, t eration continues. When the ice breaks off level ice, the edge of level ice is updat calculate the new contact force in the next iterative step. The ice-breaking process an ice loads are obtained through such iterative calculation.  The solution process of ship-ice interaction is shown in Figure 2. The physical parameters of ice and the information of ship motion are input at the beginning of the solution process. The contact between the ship and ice is judged and the interaction force is calculated. Meanwhile, the stress and strain of ice are calculated. If the pressure or plastic strain reaches the value of failure, the failure element will be deleted; otherwise, the iteration continues. When the ice breaks off level ice, the edge of level ice is updated to calculate the new contact force in the next iterative step. The ice-breaking process and the ice loads are obtained through such iterative calculation.

Calculation of Submersion Forces Fsub
The submersion forces are caused by the rotation or accumulation of broken ic its sliding along the hull, as shown in Figure 3.

Calculation of Submersion Forces F sub
The submersion forces are caused by the rotation or accumulation of broken ice and its sliding along the hull, as shown in Figure 3.

Calculation of Submersion Forces Fsub
The submersion forces are caused by the rotation or accumulation of broken ice and its sliding along the hull, as shown in Figure 3. The Lindqvist empirical formula is: where R c is the crushing force; R b is the bending force, and R s was assumed as the sum o p R , f R ,and R μ ; p R is simplified as the potential energy of ice; f R I s the frictiona force of between ice and bow; and R μ is the frictional force between ice and the ship' bottom.
where δ ρ is the density difference of ice water, g is the acceleration of gravity, H ice is the ic thickness, B is the beam, T is the draft, μ is the coefficient of friction, α is the waterlin The Lindqvist empirical formula is: where R c is the crushing force; R b is the bending force, and R s was assumed as the sum of R p , R f ,and R µ ; R p is simplified as the potential energy of ice; R f is the frictional force of between ice and bow; and R µ is the frictional force between ice and the ship's bottom.
where δ ρ is the density difference of ice water, g is the acceleration of gravity, H ice is the ice thickness, B is the beam, T is the draft, µ is the coefficient of friction, α is the waterline entrance angle, ϕ is the stem angle, and ψ is the angle between the normal surface and a vertical vector. When the icebreaker is maneuvering with small and medium drift angles, the submersion forces generated by broken ice are dominated by friction, and the contact area between the broken ice and the hull is close. To calculate the global ice loads in the transversal direction, this paper adopted the approximation in Equation (10) when calculating the submersion forces of the small drift angle maneuvering motion in level ice [18,19,24,41]. Figure 4 gives the schematic diagram of the submersion forces. [18,19,24,41]. Figure 4 gives the schematic diagram of the submersion forces. x sub x sub x y y sub y sub

Finite Element Numerical Model
This study focused on the global ice loads of the icebreakers, without considering the influences of other environmental loads, and the hull was assumed to be a rigid body. According to the assumption that the icebreakers are located in an infinite navigation field, the boundary of the calculation domain adopted rigid constraints, as shown in Figure 5. The shell and hexahedral solid mesh are used to discretize the hull and level ice, respectively. When the marine structures interact with the level ice, the level ice is usually simplified as a cuboid with uniform thickness, and the same simplification was applied to level ice in our article. For those numerical simulations, the high-quality hexahedral meshes were used to divide the level ice, and all the mesh sizes of ice were the same, as shown in Figure 5. According to the convergence test of mesh size by Subiao [43], the sensitivity analysis of mesh and ice parameters by Long Xue [44], and our [45][46][47] and other researchers' [48,49] experience in the calculation of global ice loads, to achieve a balance between the computation time and the convergence, the level ice was divided into two layers in thickness, and the mesh sizes were 300 × 500 × 500 mm, 376 × 500 × 500 mm, and 400 × 500 × 500 mm in different cases. The number of meshes and the convergence of the results met the requirements.
Due to the complex material properties of sea ice, a constitutive model of sea ice has not yet been developed. Based on the research of Wang [50] and Liu [51], the isotropic material model with failure of MAT-13 was adopted in the LS-DYNA material library, ignoring the parameters that weaken the strength of sea ice materials, such as porosity, salt degree, etc. The values of density, shear modulus, yield stress, bulk modulus, plastic failure strain, and failure pressure were determined to establish the constitutive model of sea ice ( Table 1). The values of shear modulus and bulk modulus were calculated from

Finite Element Numerical Model
This study focused on the global ice loads of the icebreakers, without considering the influences of other environmental loads, and the hull was assumed to be a rigid body. According to the assumption that the icebreakers are located in an infinite navigation field, the boundary of the calculation domain adopted rigid constraints, as shown in Figure 5. The shell and hexahedral solid mesh are used to discretize the hull and level ice, respectively. When the marine structures interact with the level ice, the level ice is usually simplified as a cuboid with uniform thickness, and the same simplification was applied to level ice in our article. For those numerical simulations, the high-quality hexahedral meshes were used to divide the level ice, and all the mesh sizes of ice were the same, as shown in Figure 5. According to the convergence test of mesh size by Subiao [43], the sensitivity analysis of mesh and ice parameters by Long Xue [44], and our [45][46][47] and other researchers' [48,49] experience in the calculation of global ice loads, to achieve a balance between the computation time and the convergence, the level ice was divided into two layers in thickness, and the mesh sizes were 300 × 500 × 500 mm, 376 × 500 × 500 mm, and 400 × 500 × 500 mm in different cases. The number of meshes and the convergence of the results met the requirements. the elastic modulus and Poisson's ratio given in the Araon model test ( Table 2). Modulus is the ratio of stress and strain under the stress state of the material, and strength is the ability of the material to resist fracture and excessive deformation. Due to the limitation of the numerical method, the modulus parameters of sea ice in numerical simulation and experiment were consistent, but the flexural and compressive strengths were not used as the failure criterion of sea ice, and the minimum pressure and maximum plastic strain were used as the failure criteria of level ice (Equation (11)). When the node satisfies the failure criterion, it will be deleted from FEM to simulate the ice breaking.
1 min (11) where p n+1 is the pressure of the n + 1 step, p min is the pressure cutoff, p eff ε is the effective plastic strain, and m ax p ε is the plastic failure strain. There are two coordinate systems, the earth coordinate system and the ship coordinate system, to describe the ship's motion, as shown in Figure 6.   Due to the complex material properties of sea ice, a constitutive model of sea ice has not yet been developed. Based on the research of Wang [50] and Liu [51], the isotropic material model with failure of MAT-13 was adopted in the LS-DYNA material library, ignoring the parameters that weaken the strength of sea ice materials, such as porosity, salt degree, etc. The values of density, shear modulus, yield stress, bulk modulus, plastic failure strain, and failure pressure were determined to establish the constitutive model of sea ice ( Table 1). The values of shear modulus and bulk modulus were calculated from the elastic modulus and Poisson's ratio given in the Araon model test (Table 2). Modulus is the ratio of stress and strain under the stress state of the material, and strength is the ability of the material to resist fracture and excessive deformation. Due to the limitation of the numerical method, the modulus parameters of sea ice in numerical simulation and experiment were consistent, but the flexural and compressive strengths were not used as the failure criterion of sea ice, and the minimum pressure and maximum plastic strain were used as the failure criteria of level ice (Equation (11)). When the node satisfies the failure criterion, it will be deleted from FEM to simulate the ice breaking. (11) where p n+1 is the pressure of the n + 1 step, p min is the pressure cutoff, ε p e f f is the effective plastic strain, and ε p max is the plastic failure strain. There are two coordinate systems, the earth coordinate system and the ship coordinate system, to describe the ship's motion, as shown in Figure 6.  (Table 2). Modulus is the ratio of stress and strain under the stress state of the material, and strength is the ability of the material to resist fracture and excessive deformation. Due to the limitation of the numerical method, the modulus parameters of sea ice in numerical simulation and experiment were consistent, but the flexural and compressive strengths were not used as the failure criterion of sea ice, and the minimum pressure and maximum plastic strain were used as the failure criteria of level ice (Equation (11)). When the node satisfies the failure criterion, it will be deleted from FEM to simulate the ice breaking.  (11) where p n+1 is the pressure of the n + 1 step, p min is the pressure cutoff, p eff ε is the effective plastic strain, and m ax p ε is the plastic failure strain. There are two coordinate systems, the earth coordinate system and the ship coordinate system, to describe the ship's motion, as shown in Figure 6.

Verification of the Numerical Method
Due to the lack of test data, the applicability and accuracy of this method were verified through the literature data on Araon and a multipurpose icebreaker. The main particulars of icebreakers are summarized in Table 3, and Figure 7 shows the geometry of Araon and a multipurpose icebreaker.

Verification of the Numerical Method
Due to the lack of test data, the applicability and accuracy of this method were verified through the literature data on Araon and a multipurpose icebreaker. The main particulars of icebreakers are summarized in Table 3, and Figure 7 shows the geometry of Araon and a multipurpose icebreaker.    Table 1. Ice tank model tests generally follow Cauchy, Ch, and Froude, Fr, similarity with a full-scall ship (Equation (12)) [53,54]. Therefore, the relationship between full-scale ship ice resistance and model ice resistance is shown in Equation (13), where Fres is the ice resistance of a numerical method and Fres-model is the model test's ice resistance, which is total resistance minus the open water resistance; the open water resistance, Fow, of Araon is given in Equation (14) [52]:   [52]. The parameters of the sea ice material are shown in Table 1. Ice tank model tests generally follow Cauchy, Ch, and Froude, Fr, similarity with a full-scall ship (Equation (12)) [53,54]. Therefore, the relationship between full-scale ship ice resistance and model ice resistance is shown in Equation (13), where F res is the ice resistance of a numerical method and F res-model is the model test's ice resistance, which is total resistance minus the open water resistance; the open water resistance, F ow , of Araon is given in Equation (14) [52]: Figure 8 shows an underwater camera view of the model test and the numerical ice-hull interaction. The model test shows how the bow breaks the level ice into pack ice, which forms accumulation and slides along the bottom of the ship; it should be noted that the numerical simulation is used to simulate ice breaking and F bre calculation. Figure 8a shows the interaction between level ice and the bow; it can be clearly observed that the level ice was broken into pack ice under the action of the bow. Figure 8b is the bottom view of the interaction phenomena between level ice and the stern. The icebreaking phenomena can still be observed in the stern, but much less than in the bow. Comparing Figure 8a,b, we see that the bow is the main icebreaking structure and the numerical simulation agrees well with the model test phenomena in terms of the breaking of level ice and the width of the channel in direct motion. It is ideal to directly compare the resistance curves obtained from the test and numerical simulation; but due to the limited test data, only the average value of resistance can be used for the verification of the method. The ice resistance predicted by the numerical method and the Araon model test results are shown in Figure 9 and Table 4. The numerical results overestimate the mean ice resistance from the model test: The maximum deviation is 17.6%, the minimum deviation is 9.4%, and the average deviation is 13.8%, which shows good agreement with the Araon model test. The main reason is that the constitutive relation and failure mode of sea ice used in the numerical simulation cannot comprehensively describe the complex material properties of actual sea ice.

Verification of Calculation Method by Icebreaker Araon Model Test
ice and the width of the channel in direct motion. It is ideal to directly compare the resistance curves obtained from the test and numerical simulation; but due to the limited test data, only the average value of resistance can be used for the verification of the method. The ice resistance predicted by the numerical method and the Araon model test results are shown in Figure 9 and Table 4. The numerical results overestimate the mean ice resistance from the model test: The maximum deviation is 17.6%, the minimum deviation is 9.4%, and the average deviation is 13.8%, which shows good agreement with the Araon model test. The main reason is that the constitutive relation and failure mode of sea ice used in the numerical simulation cannot comprehensively describe the complex material properties of actual sea ice.

Analysis of a Typical Direct Motion
The direct motion of the multipurpose icebreaker was simulated with reference to the literature [24] at 11 knots and 0.6-m ice thickness. The main parameters of the icebreaker and ship geometry model are given in Table 3 and Figure 7, respectively. The direct motion can be analyzed from the icebreaking process and mode. Observing the icebreaking resistance curve in Figure 10, there are three peak intervals, whose period is about 1.5 s. The resistance increases rapidly at the beginning of the direct motion and then fluctuates within a certain range.
The specific method of analysis is as follows: (1) At the beginning of the motion, the level ice was squeezed as it collided with the stem, which increased the icebreaking resistance. The first peak of the resistance curve appeared at 0.8 s, and bending of the level ice was observed. (2) After the broken ice fell off, the icebreaking resistance decreased until 1.5 s. (3) As the icebreaker advanced, the stem continued to collide with the level ice and squeeze it, causing the icebreaking resistance to increase again and reach a second peak at 2.3 s. (4) Broken ice started falling off; the icebreaking resistance decreased from 2.5 to 2.9 s. This process repeats continuously during direct sailing, which makes the icebreaking resistance fluctuate within a certain range. The average value of the icebreaking resistance in the peak interval of II and III is 323 KN. According to Equation (10), the submersion resistance is 577 KN, so the ice resistance under this operating condition was about 900 KN. The mean ice resistance was 823 KN in the test.

Analysis of a Typical Direct Motion
The direct motion of the multipurpose icebreaker was simulated with reference to the literature [24] at 11 knots and 0.6-m ice thickness. The main parameters of the icebreaker and ship geometry model are given in Table 3 and Figure 7, respectively. The direct motion can be analyzed from the icebreaking process and mode. Observing the icebreaking resistance curve in Figure 10, there are three peak intervals, whose period is about 1.5 s. The resistance increases rapidly at the beginning of the direct motion and then fluctuates within a certain range.
2.5 to 2.9 s. This process repeats continuously during direct sailing, which makes the icebreaking resistance fluctuate within a certain range. The average value of the icebreaking resistance in the peak interval of II and III is 323 KN. According to Equation (10), the submersion resistance is 577 KN, so the ice resistance under this operating condition was about 900 KN. The mean ice resistance was 823 KN in the test. The study of the icebreaking mode is helpful to understand the icebreaking process under different types of motion. Crushing and bending are the two main modes of icebreaking [24]. The icebreaking modes are mainly related to the shape of the structure. When the level ice collides with a slope structure (such as the stem), the level ice will bend and deform until the stress reaches the flexural strength and the level ice breaks and The specific method of analysis is as follows: (1) At the beginning of the motion, the level ice was squeezed as it collided with the stem, which increased the icebreaking resistance. The first peak of the resistance curve appeared at 0.8 s, and bending of the level ice was observed. (2) After the broken ice fell off, the icebreaking resistance decreased until 1.5 s. (3) As the icebreaker advanced, the stem continued to collide with the level ice and squeeze it, causing the icebreaking resistance to increase again and reach a second peak at 2.3 s. (4) Broken ice started falling off; the icebreaking resistance decreased from 2.5 to 2.9 s. This process repeats continuously during direct sailing, which makes the icebreaking resistance fluctuate within a certain range. The average value of the icebreaking resistance in the peak interval of II and III is 323 KN. According to Equation (10), the submersion resistance is 577 KN, so the ice resistance under this operating condition was about 900 KN. The mean ice resistance was 823 KN in the test.
The study of the icebreaking mode is helpful to understand the icebreaking process under different types of motion. Crushing and bending are the two main modes of icebreaking [24]. The icebreaking modes are mainly related to the shape of the structure. When the level ice collides with a slope structure (such as the stem), the level ice will bend and deform until the stress reaches the flexural strength and the level ice breaks and falls off, as shown in Figure 11a. When the level ice collides with a vertical structure (such as a parallel midship), or the contact area between level ice and a slope structure is too small to bend and break level ice, the breaking of the level ice is mainly via crushing, as shown in Figure 11b. falls off, as shown in Figure 11a. When the level ice collides with a vertical structure (such as a parallel midship), or the contact area between level ice and a slope structure is too small to bend and break level ice, the breaking of the level ice is mainly via crushing, as shown in Figure 11b. Figure 12 shows the process of level ice breaking. At the beginning of direct motion, the stem near the waterline began to make contact with the level ice. Meanwhile, the area of contact between the stem and the level ice was too small to break the level ice. Therefore, only crushing was observed in the time interval of 0-0.5 s, as shown in Figure 12a. With the icebreaker's advance, the contact area between the stem and level ice increased, and the contact force increased as well. At 0.6 s, the level ice began to bend downward, as shown in Figure 12b. Figure 12c shows the bending degree of level ice continuously intensifying from 0.6 s to 1s, and level ice breakup appeared at 1 s, which can be clearly observed in the red field. The level ice completed bending and fracturing at 1.1-1.5 s, and the broken ice began to fall off, as shown in Figure 12d. In summary, the breaking mode at 0-0.5 s was crushing, and the breaking mode at 0.6-1.5 s was bending, which is in good agreement with the actual icebreaking process shown in Figure 12.   Figure 12 shows the process of level ice breaking. At the beginning of direct motion, the stem near the waterline began to make contact with the level ice. Meanwhile, the area of contact between the stem and the level ice was too small to break the level ice. Therefore, only crushing was observed in the time interval of 0-0.5 s, as shown in Figure 12a. With the icebreaker's advance, the contact area between the stem and level ice increased, and the contact force increased as well. At 0.6 s, the level ice began to bend downward, as shown in Figure 12b. Figure 12c shows the bending degree of level ice continuously intensifying from 0.6 s to 1 s, and level ice breakup appeared at 1 s, which can be clearly observed in the red field. The level ice completed bending and fracturing at 1.1-1.5 s, and the broken ice began to fall off, as shown in Figure 12d. In summary, the breaking mode at 0-0.5 s was crushing, and the breaking mode at 0.6-1.5 s was bending, which is in good agreement with the actual icebreaking process shown in Figure 12.

Numerical Calculation of Ice Loads in Oblique Motion
Drift angle directly affects the ship-ice interaction area, which determines the distribution of ice loads near the waterline and ultimately determines the global forces on the hull. Therefore, it is worthwhile to research the influence of the drift angle on the ice loads on the hull, which are the factors that directly affect the ship's maneuverability in the ice [22]. The research on the characteristics of ice loads in maneuvering motion took the icebreaker Xue Long as the research object. The main particulars and geometry of the

Numerical Calculation of Ice Loads in Oblique Motion
Drift angle directly affects the ship-ice interaction area, which determines the distribution of ice loads near the waterline and ultimately determines the global forces on the hull. Therefore, it is worthwhile to research the influence of the drift angle on the ice loads on the hull, which are the factors that directly affect the ship's maneuverability in the ice [22]. The research on the characteristics of ice loads in maneuvering motion took the icebreaker Xue Long as the research object. The main particulars and geometry of the icebreaker Xue Long are shown in Table 5 and Figure 13. Considering the effects of small drift angle on ice loads, the oblique motion was simulated at 6 knots of speed, 0.8 m of ice thickness, and five drift angles (from 0 deg to 8 deg). In addition, the computational domain of the numerical model was broadened to 120 m.  Figure 13. The geometry of the icebreaker Xue Long. Figure 13. The geometry of the icebreaker Xue Long.
Box plots of ice resistance and transverse force are presented in Figures 14 and 15. With the increase in drift angle, the trend of the median line and mean value of ice resistance were not obvious, while those of transverse force increased obviously. This is closely related to the shape and velocity of the icebreaking structure. Ice resistance mainly comes from the icebreaking activity of the stem and the friction of the hull-ice interaction. With the change in drift angle, the velocity in the X-axis and the bow-ice and hull-ice contact area experienced little change, which means the median line and mean value of the box plot showed no obvious change. During oblique sailing, the transverse force mainly comes from vertical structures such as the midship and stern, and the icebreaking ability of the vertical structure is not as good as that of slope structures such as the stem. As the drift angle increases, the midship and stern participate in the icebreaking process, making the transverse force increase continuously.      Figure 16 is a comparison of icebreaking under different drift angles. It shows that the stem is always the main icebreaking structure, but the midship and stern are involved in the process as well. When comparing the level icebreaking of direct and oblique motion, we see that level ice splitting during oblique motion weakens, which is consistent with the model test under different drift angles at the Laboratory of Ice Mechanics and Ice Engineering, Tianjin University [36].

Numerical Calculation of Ice Loads in Constant Radius Motion
In a constant-radius maneuver, both drift angle and turning radius will affect the ice-hull interaction area. Hence, a constant radius maneuver was simulated to study the influence of drift angle and turning radii on ice loads at a speed of 6 knots, with 0.8 m of

Numerical Calculation of Ice Loads in Constant Radius Motion
In a constant-radius maneuver, both drift angle and turning radius will affect the ice-hull interaction area. Hence, a constant radius maneuver was simulated to study the influence of drift angle and turning radii on ice loads at a speed of 6 knots, with 0.8 m of ice thickness. The computational domain width was about 7 B, and the length was 1/8 c round.
Three turning radii (from 2 L to 3.5 L) and five drift angles (from 0 deg to 8 deg) were used in this simulation. Figures 17 and 18 are box plots of ice resistance and transverse force, respectively. When the drift angle increased, the mean value and median line of ice resistance and transverse force increased obviously; the increasing trend of the ice resistance was weaker than that of the transverse force, however. The icebreaker can open up a wider channel with a drift angle, and the stern and midship are involved in the icebreaking process, which leads to the increasing trend of transverse force being greater than that of ice resistance. The effects of drift angle and turning radius are illustrated in Figures 19 and 20. When the turning radius increased, the ice resistance and transverse force had a slight downward trend, but not a very obvious one. This may have been due to the decrease in the angular velocity of the icebreaker, which slowed down the instantaneous ship-ice interaction. Moreover, due to the limited computing resources, the small radius range may be one of the reasons why the trend was not obvious. However, with increasing drift angle, the ice resistance and transverse force both showed an obvious upward trend. It can be seen that the drift angle had a greater impact on the ice loads, which may be related to the effect of drift angle on channel width. National Research Council-Ocean, Coastal and River Engineering (NRC-OCRE) (former National Research Council of Canada-Institute for Ocean Technology, NRC-IOT) carried out a series of planar motion mechanism tests using a 1:21.8-scaled model of Terry Fox, a Canadian icebreaker, in level ice [22]. The numerical results of transverse force showed a trend similar to the model Terry Fox's test results.         Figure 21 shows the icebreaking process of 5 s, 10 s, 30 s, and 50 s when the turning radius and the drift angle were 3.5 L and 8 deg, respectively. During turning motions, the stern, midship, and stem began to break ice continuously and there was a gap on the starboard side. Figure 22 shows the ice breaking under the drift angles of 0 • , 4 • , 6 • , and 8 • . As the drift angle increased, the midship and stern were involved more in the icebreaking process, while the stem was always the main icebreaking structure.

Conclusions
The influences of the drift angle and turning radius on the global ice loads of the icebreaker Xue Long were studied based on the widely used combination of numerical simulation and empirical formulas, in which the numerical simulation computes the icebreaking processes and the Lindqvist formula estimates the underwater process. The results are summarized as follows.
First, the accuracy and applicability of the ice load calculation method were verified by the Araon model test of the KRISO ice tank in level ice. The numerical results, mean ice resistance, and icebreaking process were in good agreement with the Araon model test. In addition, the icebreaking process of a multipurpose icebreaker at a higher speed was studied, and the relationship between bending and slope structure, crushing, and vertical structure was discussed. The main mode for slope structures such as the stem is bending; however, the level ice will be crushed under the action of vertical structures, such as ship midships.
Secondly, the influences of drift angle and turning radius on the global ice loads of the icebreaker Xue Long under the maneuvering of the captive oblique and constant radius motion were studied. When the drift angle increased, the variation trend of ice resistance was not obvious, but the ice resistance of the cantilever circular motion increased obviously, and the transverse force of both the constant radius and oblique motion increased. Meanwhile, the transverse force variation of constant-radius motion at different drift angles can be compared with the Terry Fox model test at the NRC/IOT, with which the variation trend of transverse force had a good agreement. The difference in the influence of drift angle on the ice resistance of oblique and constant-radius motion should be related to the rotational angular velocity of the hull. In addition, the splitting of the level ice during oblique motion was reduced compared with the direct motion, which was consistent with the results of the model test in the ice tank of Tianjin University.
There is still limited knowledge about ice load characteristics in level ice, however. Further research should be conducted on establishing more comprehensive sea ice failure criteria and the effects of drift angle and turning radius on local ice loads in terms of maneuvering motion and the effects of high-frequency ice loads on low-frequency maneuvering motions.

Conclusions
The influences of the drift angle and turning radius on the global ice loads of the icebreaker Xue Long were studied based on the widely used combination of numerical simulation and empirical formulas, in which the numerical simulation computes the icebreaking processes and the Lindqvist formula estimates the underwater process. The results are summarized as follows.
First, the accuracy and applicability of the ice load calculation method were verified by the Araon model test of the KRISO ice tank in level ice. The numerical results, mean ice resistance, and icebreaking process were in good agreement with the Araon model test. In addition, the icebreaking process of a multipurpose icebreaker at a higher speed was studied, and the relationship between bending and slope structure, crushing, and vertical structure was discussed. The main mode for slope structures such as the stem is bending; however, the level ice will be crushed under the action of vertical structures, such as ship midships.
Secondly, the influences of drift angle and turning radius on the global ice loads of the icebreaker Xue Long under the maneuvering of the captive oblique and constant radius motion were studied. When the drift angle increased, the variation trend of ice resistance was not obvious, but the ice resistance of the cantilever circular motion increased obviously, and the transverse force of both the constant radius and oblique motion increased. Meanwhile, the transverse force variation of constant-radius motion at different drift angles can be compared with the Terry Fox model test at the NRC/IOT, with which the variation trend of transverse force had a good agreement. The difference in the influence of drift angle on the ice resistance of oblique and constant-radius motion should be related to the rotational angular velocity of the hull. In addition, the splitting of the level ice during oblique motion was reduced compared with the direct motion, which was consistent with the results of the model test in the ice tank of Tianjin University.
There is still limited knowledge about ice load characteristics in level ice, however. Further research should be conducted on establishing more comprehensive sea ice failure criteria and the effects of drift angle and turning radius on local ice loads in terms of maneuvering motion and the effects of high-frequency ice loads on low-frequency maneuvering motions.