Numerical Simulation of a Polar Ship Moving in Level Ice Based on a One-Way Coupling Method

In most previous ice–ship interaction studies involving fluid effects, ice was taken as unbreakable. Building breakable level ice on water domain is still a big challenge in numerical simulation. This paper overcomes this difficulty and presents a numerical modeling of a ship moving in level ice on the water by using a one-way CFD-DEM (computational fluid dynamics-discrete element method) coupling method. The detailed numerical processes and techniques are introduced. The ice crack propagation process including radial and circular cracks have been observed. Numerical results are compared with previous experimental data and good agreement has been achieved. The results show that water resistance is an order of magnitude smaller than ice resistance during the ice-breaking process. Ice resistance shows strong oscillation along with ice failure process, which are affected by ship speed and ice thickness significantly.


Introduction
With global warming, the sea ice has melted gradually [1], making the Arctic shipping routes more navigable [2]. More and more ships are beginning to sail in the polar ice regions (NSRA) [3]. Different from open water areas, ships would encounter not only water loads but also ice loads in the polar ice region. In particular, the effect of ice on ships presents different character from that of water, which affects the motion response and structural safety of ships. Therefore, it is of significance to understand ice-water-ship interaction and predict the ice and water loads of a polar ship more accurately, which has become one of the core problems for the design and operation of polar ships [4,5].
The problem of ice-water-ship interactions is very complex [5]. To solve this problem, various experimental, analytical, and numerical methods have been developed. Experimental studies, including field tests [6][7][8] and model tests in ice tanks [9][10][11][12][13][14], are generally viewed as reliable; however, they are strict regarding the experimental facilities and methods used and are quite expensive so cannot be done easily. Analytical methods have a relatively long history from simple models to complex models [15][16][17][18][19][20]. However, they are usually based on many simplifications of body shape and ice model, which made them hard to extend to complex ship hulls and actual ice conditions. Numerical methods [21][22][23][24][25] have been developing quickly in recent years with the rapid development of computer capacity. In contrast to experimental and analytical methods, numerical methods are easy to extend to various body configurations.
From the perspective of ice mechanics, the finite element method (FEM) and discrete element method (DEM) are two mainstream numerical methods used to simulate ice-ship interactions. FEM is a relatively mature method to solve the continuum mechanics problem and has been used in ice-ship interaction recently, especially adopted by various commercial software. Valanto [26] adopted FEM to

Fluid Model (CFD)
The fluid domain is governed by the equation of continuity and Navier-Stokes equation for an incompressible fluid, as expressed in Equations (1) and (2) [46], where η is the volume fraction of the fluid term in the control volume and it satisfies η ϕ pi V pi )/V cell , in which ε is the volume fraction of the particle term in the control volume, or the commonly called 'void ratio', V cell is the total volume of the calculated cell, n is the number of the particles in the cell, ϕ pi is the weighting coefficient of the particle i according to the ratio of the particle volume in the cell to its total particle volume V pi . u f is the fluid velocity, p is the fluid pressure, ρ f is the fluid density, ν f is the kinematic viscous coefficient of the fluid, and R p f represents the momentum exchange between the fluid phase and the particle phase. Standard k − ε turbulence model is selected in this paper. For more details on turbulent models, refer to reference [46]. There are many methods to calculate the momentum exchange term R p f and a common method is introduced here [47], which is relatively accurate and easy to implement.
where F p f is the interaction force between the fluid and particle, which is obtained by integrating the pressure of the fluid on the particle, and u p is the particle velocity. The volume of fluid (VOF) method is used to deal with the interface between water and air. The VOF model defines α w and α a as the water volume fraction and the air volume fraction, respectively. Considering volume fraction of the fluid term η= (V w + V a )/V cell , where V w and V a are water volume and air volume in the cell, α i is defined as α i = V i / 2 i=1 V i , where phase i = 1 and 2 denote water and air, respectively. Thus, in a control volume, it satisfies If the cell is full of water, one has α w = 1. If there is no water in the cell, one has α w = 0. Otherwise, one has 0 < α w < 1. The interface is tracked by solving the volume fraction transport equation [48], as below: (5) where A is the surface of the control volume, n is the unit normal vector of the surface, S α i is user-defined source term for phase i, v dr,i is the diffusion velocity for phase i. Ship is viewed as a rigid body upright floating in the water at a given draught in this simulation, which can move at a given speed in just one direction, and the hull form comes from an icebreaker. The principal dimensions of the hull are shown in Table 1. In numerical simulations of a ship moving in a brash ice region, to make the numerical modeling easier, part of previous work [43] kept the ship stationary, set water moving in a constant velocity against the hull, and released the ice particles into the water with the same initial velocity. In this way, water carried ice particles towards the ship, which was sometimes used as an alternative in simulating a ship moving in a brash ice region. However, this does not work for level ice, given that the push of water on the level ice is not large enough to support the collision between ship and level ice, inducing results wildly inconsistent with reality. In order to solve this problem, it is necessary to set the ship to move in and break the level ice floating on the water.
The overset grid technique is used in this paper to achieve this process. The overset grid is composed of two sets of grids, the background grid, and the moving grid. Information is exchanged between these two sets of grids using interpolation methods. In this paper, the linear interpolation method is used. The grid accepting information is named as 'acceptor', while the grid providing information is named as 'provider'. In linear interpolation formulation, one has where Φ is the physical quantity, χ is weighting factor, and the subscript i denotes the number of the provider. As shown in Figure 1, the physical quantity of a grid C comes from those of the neighboring grids N 1 , N 2 , and N 3 of the same set of grids and the neighboring grids N 4 , N 5 , and N 6 of the other set of grids. The equation is solved iteratively until the residual error is small enough. During the construction, it should be noted that the overlapping area between two sets of grids should contain at least 4-5 grid cell layers, and the boundary mesh sizes should be similar, which reduces the error of interpolation as much as possible. More information about overset grid technique can refer to Hadzic [49]. The flow field is truncated as far as possible in order to reduce the influence of the outer boundary of the flow field. In this paper, the flow field is taken as 600 m × 600 m × 450 m. The fluid domain is meshed with a trimmed meshing model and boundary layer grids are divided around the ship surface, as shown in Figure 2. According to the ship speed, the Reynolds number in this paper is 1.94 × 10 8 ≤ Re ≤ 4.54 × 10 8 , which is a high-Reynolds-number problem. According to the user's manual of STAR-CCM+ [48], y+ for a high-Reynolds-number problem is recommended to be larger than 30 in standard k − ε turbulence model. Thus, in this paper, the thickness of the first layer grid of the boundary layer is 8 × 10 −4 m, and the y+ value is chosen around 60 based on our numerical experience. The meshes on the ship surface are refined near waterline, stem, and stern area, as shown The flow field is truncated as far as possible in order to reduce the influence of the outer boundary of the flow field. In this paper, the flow field is taken as 600 m × 600 m × 450 m. The fluid domain is meshed with a trimmed meshing model and boundary layer grids are divided around the ship surface, as shown in Figure 2. According to the ship speed, the Reynolds number in this paper is 8 8 1.94 10 Re 4.54 10 × ≤ ≤ × , which is a high-Reynolds-number problem. According to the user's manual of STAR-CCM+ [48], y+ for a high-Reynolds-number problem is recommended to be larger than 30 in standard k ε − turbulence model. Thus, in this paper, the thickness of the first layer grid of the boundary layer is    The flow field is truncated as far as possible in order to reduce the influence of the outer boundary of the flow field. In this paper, the flow field is taken as 600 m × 600 m × 450 m. The fluid domain is meshed with a trimmed meshing model and boundary layer grids are divided around the ship surface, as shown in Figure 2. According to the ship speed, the Reynolds number in this paper is 8 8 1.94 10 Re 4.54 10 × ≤ ≤ × , which is a high-Reynolds-number problem. According to the user's manual of STAR-CCM+ [48], y+ for a high-Reynolds-number problem is recommended to be larger than 30 in standard k ε − turbulence model. Thus, in this paper, the thickness of the first layer grid of the boundary layer is   The truncated surfaces of the flow domain are defined as front, rear, left, right, top, and bottom surfaces based on the ship, where the front surface is the surface that ship bow points. The boundary conditions of the flow domain include a velocity inlet for the rear, left, right, top, and bottom surfaces, and a pressure outlet for the front surface. Rigid wall boundary condition is exerted on the ship surface, which satisfies both impermeable condition and nonslip condition, namely, where u f and u s are the fluid velocity and ship velocity, respectively, and n and τ are unit normal and tangential vectors of the ship surface. On the ice surfaces, no boundary conditions are satisfied because the coupling method is based on force and moment exchange balance and not the interface track or capture method [50]. More details about the coupling method will be discussed in Section 2.3. The initial conditions of the flow field are stationary.
An appropriate number of grids should be chosen to ensure the calculation accuracy and save calculation time. After grid independence verification based on the water resistance of the ship in 5 kn, as shown in Table 2, the overall grid is taken about 1.63 million, namely mesh3 in Table 2.

Ice Model (DEM)
Ice is modeled by using DEM, and the governing equation is Newton's second law. Taking a particle element i as an example, one has where subscript i, j, and k denote the particle number, m is the particle mass, F c and F lr are the contact and non-contact forces between particles, respectively, F p f is the interaction force between the fluid and particle, same to which in Equation (3), and F g is the gravity force of the particle. I is the moment of inertia of the particle, ω is angular velocity of the particle, M t and M r are the moments of the sliding friction and the rolling friction, respectively. Considering F lr , such as electromagnetic force or van der Waals forces, is just involved in special cases, F lr is not included in this paper. Here the method to calculate the contact force F c is introduced briefly. There are many methods to calculate the contact force F c , and a common model 'Hertz-Mindlin' collision model is adopted here [47]. In the Hertz-Mindlin collision model, the contact force F c between two spheres, A and B, are described by the following set of equations: where subscript n and t denote normal and tangential components, respectively, and n and τ are unit vector in normal and tangential directions, respectively. Normal and tangential forces F n and F t are expressed as: where d and v are relative displacement and velocity, respectively. Figure 3 provides the normal and tangential relative displacements d n and d t , and they are obtained by the integral of normal and tangential relative velocities v n and v t , which are obtained by using Equation (9). K and N are spring stiffness and damping, respectively, and C f s is a static friction coefficient, which is taken as 0.15 according to the recommended value for sea ice [51]. K and N are not constant in Hertz-Mindlin collision model and they are calculated each step by using K n = 4 3 E eq d n R eq , K t = 8G eq d t R eq , N n = 5K n M eq N n damp , N t = 5K t M eq N t damp , in which equivalent Young modulus E eq = with ν as Poisson's ratio and subscript A and B denoting particles number, equivalent radius R eq = , equivalent shear , equivalent mass M eq = , normal damping coefficient π 2 +ln (C n rest ) 2 and tangential damping coefficient , in which C n rest and C t rest are the normal and tangential coefficients of restitution, which are taken as 0.5 [45]. , equivalent shear modulus  The generation of level ice is one of the main difficulties in the process of ship-level ice-water interaction. Care must be taken to build breakable level ice in DEM. In STAR-CCM+, there are a variety of particle models in the DEM module, and the composite particle model is widely used because it can be freely combined into the shape of the target object. However, the composite particle model cannot be broken, which is inconsistent with the properties of the level ice in this paper. On the other hand, the particle clumps model can realize the fracture between particles, but it requires The generation of level ice is one of the main difficulties in the process of ship-level ice-water interaction. Care must be taken to build breakable level ice in DEM. In STAR-CCM+, there are a variety of particle models in the DEM module, and the composite particle model is widely used because it can be freely combined into the shape of the target object. However, the composite particle model cannot be broken, which is inconsistent with the properties of the level ice in this paper. On the other hand, the particle clumps model can realize the fracture between particles, but it requires initial input of the position and size of each particle to define the initial shape of the object. It would be a huge amount of work to generate such a model of level ice with tens of thousands of particles. As an alternative, in this paper, a simple spherical particle model is adopted, and the bonding bond and fracture criterion between the particles are applied to establish the model of breakable level ice, which will be proved to be successful and relatively easy. The model is chosen as viscoelastic with ice density of 900 kg/m 3 and Young's modulus of 1 Gpa. The detailed processes are described. It includes two main steps: 'particle generation' and 'particle bonding'. For the first step, particles are generated based on seed points. For the success of bonding a level ice, it is necessary to ensure that the particles are relatively compact and no large gaps exist. Generating particles with a conventional 'part ejector' does not guarantee needed particle density. To solve this problem, a 'random ejector' with the maximum filling form is used. This 'random ejector' injects a specified number of particle seeds with small starting volume into the space in random positions, in which the number is obtained by the area of the space and the prescribed radius of the particle. Then, spherical particles start to grow from the seeds. As the maximum filling form is adopted, the spherical particles grow until there is no room left. For the 'random ejector', the space to fill must be real. A direct idea is to build a real space where the level ice locates, before filling it. However, building such a real space in the fluid domain will affect the computation of the fluid. To avoid this disturbance, a technique was developed to create a new space in a spare fluid domain. This space and this fluid domain are not used for calculation, but just for generating ice particles. After the particles are generated by 'random ejector' with the maximum filling form there, the positions and sizes of these particles can be derived to make a table, which consists of all the information of the particles. Then, the original fluid domain is returned, and another ejector, 'table ejector', is used to generate ice particles at the above locations and sizes in the original fluid domain. There is another problem with this method. The 'table ejector' does not generate ice particles simultaneously, which makes the early ejected ice particles deviate from initial positions when all the ice particles have been ejected, inducing the bonding step to fail. Therefore, another technique was put forward, which is to use a custom function to generate these particles simultaneously. In this way, all the particles are kept at their initial positions before bonding, laying a good foundation for a successful bonding step subsequently. These two techniques above are of significance to this step.
After all the ice particles have been released to their prescribed positions in the first step 'particle generation', they need to be bonded together as a level ice. Thus, the second step, 'parallel bonding', is implemented. As shown in Figure 3a, the parallel bonding model uses the concept of massless bars connecting a pair of bonded particles. The bars can transmit force and torque between particles. Alongside force Equation (11), torque equation is expressed as: where M n and M t denote normal and tangential torque components, respectively. One has where ∆ is increment, J = 1 2 πR 4 eq and L = 1 4 πR 4 eq , K is still spring stiffness and the value is the same to that in Equations (12) and (13), Ω is relative angular displacement.
Based on 'parallel bonding' model, a simple rupture model is used. It means that once the tensile or shear stresses between particles exceeds the maximum limits, the bar breaks and corresponding force and torque disappear. Then the particles separate, and cracks generate. In this way, a breakable ice model becomes available. It should be stated that the interparticle tensile and shear bonding strength of ice in the DEM model are affected by many factors, including particle shape, particle radius, number of layers, and arrangement modes (regular or random arrangement) [52]. According to previous studies [53], the interparticle tensile and shear bonding strength of ice in random arrangement mode should be larger than macroscopic measured tensile and shear strength of sea ice [54], so that the simulated ice behaviors including compression and bending agree with those of real ice. Based on the suggested values of reference [53] and our numerical experience, interparticle tensile and shear bonding strength were both taken as 3.0 Mpa. The final successful bonded level ice model with two layers of ice particles is shown in Figure 4, which is also adopted in this paper.  [54], so that the simulated ice behaviors including compression and bending agree with those of real ice. Based on the suggested values of reference [53] and our numerical experience, interparticle tensile and shear bonding strength were both taken as 3.0 Mpa. The final successful bonded level ice model with two layers of ice particles is shown in Figure 4, which is also adopted in this paper.

Coupling Scheme
As mentioned in Luo et al. [45] and Ni et al. [5], there are usually two kinds of coupling schemes in the coupling framework of CFD and DEM: two-way coupling (TWC) and one-way coupling (OWC). TWC considers both the full interaction between fluid (solved by CFD) and particles (solved by DEM), while OWC just considers the force of fluid on the particle but ignores the force of ice on the fluid. In other words, fluid exerts forces on the particles, but the particles do not affect the movement of the fluids. Luo et al. [45] adopted both TWC and OWC methods to calculate a ship moving in brash ice channel and compared the influences of these two schemes. It was found that

Coupling Scheme
As mentioned in Luo et al. [45] and Ni et al. [5], there are usually two kinds of coupling schemes in the coupling framework of CFD and DEM: two-way coupling (TWC) and one-way coupling (OWC).
TWC considers both the full interaction between fluid (solved by CFD) and particles (solved by DEM), while OWC just considers the force of fluid on the particle but ignores the force of ice on the fluid. In other words, fluid exerts forces on the particles, but the particles do not affect the movement of the fluids. Luo et al. [45] adopted both TWC and OWC methods to calculate a ship moving in brash ice channel and compared the influences of these two schemes. It was found that although the minimum error of TWC was slightly smaller than that of OWC relative to experimental resistance, the average error of them was very close. However, TWC occupied much more computing resources than OWC. The total solving CPU time of TWC rose more sharply and became much longer, up to five times or longer, than that of OWC along with the calculation [55]. As a result, Luo et al. [45] recommended OWC for the case of a ship moving in ice at a low speed. Considering the speed of a ship moving in level ice is much lower than that in brash ice, the OWC method is chosen in the simulation in this paper.
Here the framework of OWC is introduced briefly, as shown in Figure 5, and that of TWC can refer to [34]. To start with, all components of simulation, including DEM, CFD and coupling parts, are initialized. All the initial conditions and initial value, including stationary fluid with hydrostatic pressure, stationary level ice floating on the water surface and rigid ship moving forward in a given speed, are assigned. The OWC starts with calculating the fluid porosity, namely the ratio of the particle volume in each fluid cell, based on the particle information (position and size) and fluid mesh information (mesh size and node location). As mentioned in Section 2.1, the interface between water and ice particle is not tracked or captured, and the fluid-particle interaction force and momentum exchange are calculated based on fluid porosity, so it is important to calculate fluid porosity. Thereafter, velocity of particles and fluid as well as the pressure and stress tensor of fluid at the current time step are used to calculate the fluid-particle interaction force F p f ,i in Equation (9). The next step is the iteration loop of the DEM. After the DEM loop involving Equations (9) and (10) (1) and (2) are solved. In OWC, R p f in Equation (2) is taken as zero. In this way, the solving process becomes much easier and enormous computation time and storage are saved. Then the obtained particle and fluid information will be used in next loop, until the whole simulation terminates.

Validation
Before studying the motion of a ship in the level ice, the coupling method is firstly validated by comparing numerical results with experimental data. There have been previous studies of model ships moving in level ice in an ice tank as mentioned in Section 1. Considering the hull form

Validation
Before studying the motion of a ship in the level ice, the coupling method is firstly validated by comparing numerical results with experimental data. There have been previous studies of model ships moving in level ice in an ice tank as mentioned in Section 1. Considering the hull form information, the experimental study from ice tank of Tianjin University [14] was chosen. The corresponding prototype used by Huang et al. [14] is that same as that used in this paper, as shown in Table 1 and Figure 2, and the scale ratio is 37.5 for their model. Considering that they have predicted the resistance of the prototype from model resistance under the similarity laws, we calculated the prototype and compared our numerical results with their predicted resistance of prototype. The experimental picture and the state of the ship before moving into level ice region in numerical simulation are shown in Figure 6a,b, respectively.

Validation
Before studying the motion of a ship in the level ice, the coupling method is firstly validated by comparing numerical results with experimental data. There have been previous studies of model ships moving in level ice in an ice tank as mentioned in Section 1. Considering the hull form information, the experimental study from ice tank of Tianjin University [14] was chosen. The corresponding prototype used by Huang et al. [14] is that same as that used in this paper, as shown in Table 1 and Figure 2, and the scale ratio is 37.5 for their model. Considering that they have predicted the resistance of the prototype from model resistance under the similarity laws, we calculated the prototype and compared our numerical results with their predicted resistance of prototype. The experimental picture and the state of the ship before moving into level ice region in numerical simulation are shown in Figure 6a,b, respectively. Two main physical quantities are concerned. One is the damage of the level ice under impact of the ship, including the crack development and broken ice movement. The other is the total resistance. The former can just be compared qualitatively while the later can be compared quantitatively. In order to validate the former, the case with ice thickness 1.5 m and ship speed 5 kn in prototype is chosen as an example. The damage of ice obtained by numerical simulation in this paper is compared with that in the model test, as shown in Figure 7. Comparison of ice damage between model tests (a) Figure 6. Bird's eye view of experimental picture (a) (reproduced from [14], with permission from ELSEVIER, 4 September 2020), and the state of the ship before moving into level ice region in numerical simulation (b).
Two main physical quantities are concerned. One is the damage of the level ice under impact of the ship, including the crack development and broken ice movement. The other is the total resistance. The former can just be compared qualitatively while the later can be compared quantitatively. In order to validate the former, the case with ice thickness 1.5 m and ship speed 5 kn in prototype is chosen as an example. The damage of ice obtained by numerical simulation in this paper is compared with that in the model test, as shown in Figure 7. Comparison of ice damage between model tests (a) [14] and (b) [56] and numerical simulation (c) is examined, where the general outline of the ice crack is marked by lines.  [14] and (b) [56] and numerical simulation (c) is examined, where the general outline of the ice crack is marked by lines. In Figure 7, two cases of experimental data are chosen for comparison as shown in Figure 7a,b, respectively. Figure 7a provides the experimental data from [14], where the red lines are the outlines of the circular cracks. Experimental scene was taken by camera, while the model vessel was driven back to prepare the next test run. Figure 7c is the numerical result with the same parameters in Figure  7a. In Figure 7c, ice particles in the same fragment are shown in the same color, so different colors can be seen to distinguish the shapes of the broken ice. Considering the cracks in Figure 7a from [14] are not clear enough, although we tried to highlight them by red lines, the picture from another reference [56] is adopted as supplementary in Figure 7b. It is worth mentioning that there is no direct link between numerical simulation and model test in Figure 7b, as the ship forms are different. As a result, the characteristics of the ice breakup in numerical simulation Figure 7c are compared qualitatively with those in Figure 7a,b. It can be found that the simulated ice destruction area and the shape of the broken ice are similar to those of the model test. There are mainly two points of qualitative similarities. One is the ice-breaking areas extending outwards in a V-shape, as shown by the general outlines of the crack in Figure 7a,c. Circular cracks and radial cracks cross each other, forming crushed ice of various sizes. The other similarity is that crushed ice has a smaller size when it is closer to the center line of the ship, especially shown by the fragments in Figure 7b,c. Secondly, the ice resistance values of numerical simulation are compared with those transformed from model test, as plotted in Figure 8. In Figure 8, the squares denote the average resistance at ship velocity 2 kn, 3 kn, 4 kn, and 5 kn in prototype transformed from model test, while the dots denote the average resistance at ship velocity 3 kn, 5 kn, and 7 kn in numerical simulation. Considering the ship with lower speed taking much longer calculation time and resources, we did  [14], with permission from ELSEVIER, 4 September 2020), and (b) (reproduced from [56], with permission from ELSEVIER, 4 September 2020), and numerical simulation (c), where the general outline of the ice crack is marked by lines.
In Figure 7, two cases of experimental data are chosen for comparison as shown in Figure 7a,b, respectively. Figure 7a provides the experimental data from [14], where the red lines are the outlines of the circular cracks. Experimental scene was taken by camera, while the model vessel was driven back to prepare the next test run. Figure 7c is the numerical result with the same parameters in Figure 7a. In Figure 7c, ice particles in the same fragment are shown in the same color, so different colors can be seen to distinguish the shapes of the broken ice. Considering the cracks in Figure 7a from [14] are not clear enough, although we tried to highlight them by red lines, the picture from another reference [56] is adopted as supplementary in Figure 7b. It is worth mentioning that there is no direct link between numerical simulation and model test in Figure 7b, as the ship forms are different. As a result, the characteristics of the ice breakup in numerical simulation Figure 7c are compared qualitatively with those in Figure 7a,b. It can be found that the simulated ice destruction area and the shape of the broken ice are similar to those of the model test. There are mainly two points of qualitative similarities. One is the ice-breaking areas extending outwards in a V-shape, as shown by the general outlines of the crack in Figure 7a,c. Circular cracks and radial cracks cross each other, forming crushed ice of various sizes. The other similarity is that crushed ice has a smaller size when it is closer to the center line of the ship, especially shown by the fragments in Figure 7b,c.
Secondly, the ice resistance values of numerical simulation are compared with those transformed from model test, as plotted in Figure 8. In Figure 8, the squares denote the average resistance at ship velocity 2 kn, 3 kn, 4 kn, and 5 kn in prototype transformed from model test, while the dots denote the average resistance at ship velocity 3 kn, 5 kn, and 7 kn in numerical simulation. Considering the ship with lower speed taking much longer calculation time and resources, we did not calculate cases 2 kn and 4 kn and took cases 3 kn and 5 kn as examples for comparisons. As can be seen from Figure 8, when the speed is 5 kn, the mean value of numerical simulation resistance is about 1.9 MN, which is very close to the measured resistance 2 MN in ice tank, with the deviation just around 5%. The total ice resistance deviation is only 2% when the speed is 3 knots. The good agreement validates the accuracy of the coupling method to some extent. Based on the above analysis, it can be seen agreement has achieved between numerical results and experimental data both qualitatively and quantitatively. In this way, it is considered that the coupling numerical model can simulate ship navigation in level ice effectively.

Results and Discussions
The one-way CFD-DEM coupling model established in this paper is applied to the numerical simulation of a ship moving in the level ice region. The responses of ice in this process, including crack propagation, are obtained. In addition, the characteristics of the resistance of the ship in this Based on the above analysis, it can be seen agreement has achieved between numerical results and experimental data both qualitatively and quantitatively. In this way, it is considered that the coupling numerical model can simulate ship navigation in level ice effectively.

Results and Discussion
The one-way CFD-DEM coupling model established in this paper is applied to the numerical simulation of a ship moving in the level ice region. The responses of ice in this process, including crack propagation, are obtained. In addition, the characteristics of the resistance of the ship in this process are also one focus of this paper.

Ice-Breaking Process
The interaction between a ship and level ice is a very complicated problem. The process of ice breaking and crack propagation is accompanied by various failure modes. To study this problem, the case with ship speed of 5 kn and the ice thickness of 1.5 m is selected as a typical case study in this section. Figure 9 shows the ice evolution during the ice-breaking process. As mentioned before, different colors denote different ice fragments. Figure 8. Comparison of ice resistance between model test (the squares), transformed into prototype by Huang et al. [14] and numerical simulation (the dots) at various ship velocity. Based on the above analysis, it can be seen agreement has achieved between numerical results and experimental data both qualitatively and quantitatively. In this way, it is considered that the coupling numerical model can simulate ship navigation in level ice effectively.

Results and Discussions
The one-way CFD-DEM coupling model established in this paper is applied to the numerical simulation of a ship moving in the level ice region. The responses of ice in this process, including crack propagation, are obtained. In addition, the characteristics of the resistance of the ship in this process are also one focus of this paper.

Ice-Breaking Process
The interaction between a ship and level ice is a very complicated problem. The process of ice breaking and crack propagation is accompanied by various failure modes. To study this problem, the case with ship speed of 5 kn and the ice thickness of 1.5 m is selected as a typical case study in this section. Figure 9 shows the ice evolution during the ice-breaking process. As mentioned before, different colors denote different ice fragments. As mentioned in Section 2.2, ice is taken as viscoelastic model. When the ship comes into contact with the ice, the elastic deformation of the ice occurs first, and no cracks appear at this time. As the ship continues forward, the contact force between ship and ice increases further. When the stress inside the ice exceeds its ultimate stress, the particle bond of the ice in DEM model breaks. With the increasing of the broken bonds, the crack of the ice starts to expand. Radial cracks are first observed As mentioned in Section 2.2, ice is taken as viscoelastic model. When the ship comes into contact with the ice, the elastic deformation of the ice occurs first, and no cracks appear at this time. As the ship continues forward, the contact force between ship and ice increases further. When the stress inside the ice exceeds its ultimate stress, the particle bond of the ice in DEM model breaks. With the increasing of the broken bonds, the crack of the ice starts to expand. Radial cracks are first observed propagating from the contact points between ship and the ice sheet, as Figure 9a t = 9 s shows. Along with the forward movement of the ship, radial cracks extend outwards gradually. A first-order circular crack appears at the ship's shoulder, as Figure 9b t = 14 s shows, and a small amount of ice fragments appear at the bottom of the bow, forming the so-called local crushing zone [13]. As the ship continues to move forward, the interaction between the ship's shoulder and the ice sheet gets severer. The ice bends under the action of ships. The first-order circular crack continues to extend forward, inducing large pieces of ice separating from the ice sheet, as shown at Figure 9c t = 15s. These ice fragments overturn under the bow and may collide with ship bottom or other ice fragments and get broken again, resulting in several smaller pieces of ice. Then a second-order circular crack appears on the basis of the first-order circular crack around ship bow, as shown at t = 17 s in Figure 9d. Many pieces of crushed ice with various sizes are formed between the radial cracks and the circular cracks. As the ship continues to move further, more radial and circular cracks generate in a similar and repeated ways as described above, which have also been observed in experiments [13]. Finally, a V-shaped crushing region will be formed as shown in Figure 9b.

Total Resistance and Ice Resistance
The total resistance is one of the most significant problems for ship performance in ice regions. It is a common way to divide the total resistance into ice resistance and water resistance [57]. Although water resistance is usually small compared with ice resistance, especially for a ship moving in level ice, water effects cannot be easily ignored. On the one hand, the ice sheet is not fixed in the direction of gravity and needs the support of water to keep it afloat and stable. On the other hand, the broken ice fragments usually drift and overturn under the action of water, which affects the change of ice resistance, especially the so-called submersion component [58]. That is also the reason why we choose the coupling model, which considers the interaction between ship, ice and water.
In the numerical modeling, on the ship surface, there are two forces, one is the fluid force and the other is the ice force. Therefore, one can obtain two resistances from the integral of fluid and ice forces, respectively. They are defined as 'water resistance' and 'ice resistance' in this paper. Because OWC method is adopted in this paper, fluid motion and pressure are not affected by the ice, so the 'water resistance' is not affected by ice either. 'Ice resistance' is affected by the fluid motion under OWC method, so it can be seen as the ice resistance considering the influence of fluid. This division method is easy to compare the contribution of ice and water in a rough way. Therefore, water and ice resistances are checked separately. Figure 10a,b provides curves of water resistance and ice resistance, respectively. As shown in Figure 10a, one can see that water resistance when stable is about 0.1-0.2 MN, which is at least an order of magnitude smaller than the ice resistance shown in Figure 10b, which is in MN. This verifies again that the water resistance component is very small in the total resistance for a ship moving in level ice. Considering that the one-way coupling method is adopted in this paper, the water resistance is not affected by ice and remains the same under different ice conditions. Therefore, the following discussion will focus on the characteristics of ice resistance instead of water resistance.  Figure 10b shows the time-history curve of ice resistance. Corresponding to Figure 9, the influence of various ice failure patterns on the ice resistance can be vividly seen in the curve. Firstly, it can be seen from the figure that at the first 9 s, the ice resistance is rising gradually in a small extent. Before 9 s, the ice is bent under the moment of ship bow. Local large elastic deformations occur but not the whole cracks. At about 9 s, the first crack starts to form and propagate, as shown in Figure 9a, so the ice resistance reaches its first peak before it decreases sharply along with the unloading of the ice force. As the ship moves further and contacts with more ice, the ice resistance rises slowly with fluctuations due to cracks. Around 14 s, the ice resistance reaches another peak, when the first-order circular crack just generates, as shown in Figure 9b. Then the ice is completely broken and the circular crack has been completely formed around 15 s, as shown in Figure 9c, so the ice resistance reaches bottom sharply due to unloading of the ice force. The generation of secondary circular crack propagation in Figure 9d also induces a rise in the resistance curve but with a smaller magnitude than that induced by the first order circular crack. Then there exists a distinct rise of ice resistance  Figure 10b shows the time-history curve of ice resistance. Corresponding to Figure 9, the influence of various ice failure patterns on the ice resistance can be vividly seen in the curve. Firstly, it can be seen from the figure that at the first 9 s, the ice resistance is rising gradually in a small extent. Before 9 s, the ice is bent under the moment of ship bow. Local large elastic deformations occur but not the whole cracks. At about 9 s, the first crack starts to form and propagate, as shown in Figure 9a, so the ice resistance reaches its first peak before it decreases sharply along with the unloading of the ice force. As the ship moves further and contacts with more ice, the ice resistance rises slowly with fluctuations due to cracks. Around 14 s, the ice resistance reaches another peak, when the first-order circular crack just generates, as shown in Figure 9b. Then the ice is completely broken and the circular crack has been completely formed around 15 s, as shown in Figure 9c, so the ice resistance reaches bottom sharply due to unloading of the ice force. The generation of secondary circular crack propagation in Figure 9d also induces a rise in the resistance curve but with a smaller magnitude than that induced by the first order circular crack. Then there exists a distinct rise of ice resistance around 18 s, when the whole ship bow enters the ice region. After that moment, although the ice resistance still fluctuates strongly, its mean value tends to be stable. This indicates that the parallel middle body of the ship has a smaller contribution to the ice resistance than bow, because the force resulted from ice collision in this area is perpendicular to the motion direction. Although the frictional resistance between broken force and hull surface also has a contribution, it is small compared with ice-ship collision force. To compare with the mean ice resistance in model test mentioned above, the mean value after the whole ship bow enters the level ice region is defined as the average ice resistance here, as denoted by the red dotted line in Figure 10.

Effect of Ship Speed
Ship speed is a key factor that affects ice breaking and ice resistance. Other parameters are kept the same as those in Section 4.1, and ship speed is changed from 3, 5, to 7 kn to determine its effect. Figure 11 shows the ice conditions under different ship speed, in which the distances of the ship into the level ice are the same with different moments marked below each figure. When the ship speed is 3 kn, the damage range of the ice sheet in the lateral direction is largest. The ice far away from the ship side is also affected and damaged, with larger pieces of broken ice and longer cracks propagated. As the ship speed increases to 7 kn, it can be found that the ice-breaking channel is narrowest. The ice beyond ship breadth has been little affected by the ship motion. Furthermore, the difference in size between the pieces of broken ice becomes smaller and the very large ice fragments get fewer. speed is 3 kn, the damage range of the ice sheet in the lateral direction is largest. The ice far away from the ship side is also affected and damaged, with larger pieces of broken ice and longer cracks propagated. As the ship speed increases to 7 kn, it can be found that the ice-breaking channel is narrowest. The ice beyond ship breadth has been little affected by the ship motion. Furthermore, the difference in size between the pieces of broken ice becomes smaller and the very large ice fragments get fewer. The change of ship speed can be considered as the change of the loading rate on the ice. As one may know, the responses of the ice are affected by the loading rate significantly [54]. When ship speed increases, the loading rate of the ice increases and the ice sheet presents a greater brittle response [54]. As a result, the size of ice fragments gets small, the crack propagation becomes short, and the lateral damage area reduces, and the ice-breaking channel gets relatively narrow.
Figure12 shows the time histories and typical values of ice resistance at different ship speed. From time history curves, one can still see two stages for each velocity. Although the time when they enter the second stage differs, the distance they take is the same, which is just equal to the length of the ship bow. This verifies the analysis in Section 4.2 again. The variation of maximum and average ice resistance in the second stage of each case are further checked in Figure 12d. It can be seen that ship speed has a significant effect on the ice resistance. Both maximum and mean values of ice resistance increase with the increase of speed, almost in a linear relationship. Similar trends are also observed in the previous experimental studies [14,56]. Furthermore, the rising slopes of maximum and mean values are close. The change of ship speed can be considered as the change of the loading rate on the ice. As one may know, the responses of the ice are affected by the loading rate significantly [54]. When ship speed increases, the loading rate of the ice increases and the ice sheet presents a greater brittle response [54]. As a result, the size of ice fragments gets small, the crack propagation becomes short, and the lateral damage area reduces, and the ice-breaking channel gets relatively narrow. Figure 12 shows the time histories and typical values of ice resistance at different ship speed. From time history curves, one can still see two stages for each velocity. Although the time when they enter the second stage differs, the distance they take is the same, which is just equal to the length of the ship bow. This verifies the analysis in Section 4.2 again. The variation of maximum and average ice resistance in the second stage of each case are further checked in Figure 12d. It can be seen that ship speed has a significant effect on the ice resistance. Both maximum and mean values of ice resistance increase with the increase of speed, almost in a linear relationship. Similar trends are also observed in the previous experimental studies [14,56]. Furthermore, the rising slopes of maximum and mean values are close.

Effect of Ice Thickness
Ice thickness is also an important factor affecting ice breaking. The ice sheet in this paper is modeled by using DEM particles. The mechanical strength of the ice sheet is determined by the bonding bars between the particles. The values of the bonding parameters are closely related to the particle size. In order to eliminate this error, the particle size is kept unchanged, and the ice thickness is changed by increasing or decreasing the number of particle layers. Considering the diameter of the particle is 0.75 m, single-layer, double-layer, and three-layer ice can be modeled, so the ice thickness is 0.75 m, 1.5 m, and 2.25 m, respectively. Other parameters are kept the same with those in Section 4.1. Figure 13 shows the ice conditions under different ice thicknesses at the same time t = 20 s. It can be found that when the ice thickness is 0.75 m, there are more long radial cracks formed on the ice surface. The ice fragments are relatively larger in size and the total amount decreases instead. As the ice grows to 1.5 m in thickness, the size of the crushed ice varies. Both large and small ice fragments exist, and the amount of crushed ice increases compared with 0.75 m case. When the ice thickness rises to 2.25 m, the damage of ice sheet reduces distinctively. Many small-sized ice fragments are generated just around the ship-ice contact region.

Effect of Ice Thickness
Ice thickness is also an important factor affecting ice breaking. The ice sheet in this paper is modeled by using DEM particles. The mechanical strength of the ice sheet is determined by the bonding bars between the particles. The values of the bonding parameters are closely related to the particle size. In order to eliminate this error, the particle size is kept unchanged, and the ice thickness is changed by increasing or decreasing the number of particle layers. Considering the diameter of the particle is 0.75 m, single-layer, double-layer, and three-layer ice can be modeled, so the ice thickness is 0.75 m, 1.5 m, and 2.25 m, respectively. Other parameters are kept the same with those in Section 4.1. Figure 13 shows the ice conditions under different ice thicknesses at the same time t = 20 s. It can be found that when the ice thickness is 0.75 m, there are more long radial cracks formed on the ice surface. The ice fragments are relatively larger in size and the total amount decreases instead. As the ice grows to 1.5 m in thickness, the size of the crushed ice varies. Both large and small ice fragments exist, and the amount of crushed ice increases compared with 0.75 m case. When the ice thickness rises to 2.25 m, the damage of ice sheet reduces distinctively. Many small-sized ice fragments are generated just around the ship-ice contact region.
The ice resistance under different ice thickness is shown in Figure 14. It can be found that both maximum and mean ice resistances increase with ice thickness sharply. As the ice thickness rises, the difference between the maximum and average ice resistance increases also. It denotes that maximum ice resistance is more sensitive to ice thickness.
be found that when the ice thickness is 0.75 m, there are more long radial cracks formed on the ice surface. The ice fragments are relatively larger in size and the total amount decreases instead. As the ice grows to 1.5 m in thickness, the size of the crushed ice varies. Both large and small ice fragments exist, and the amount of crushed ice increases compared with 0.75 m case. When the ice thickness rises to 2.25 m, the damage of ice sheet reduces distinctively. Many small-sized ice fragments are generated just around the ship-ice contact region. The ice resistance under different ice thickness is shown in Figure 14. It can be found that both maximum and mean ice resistances increase with ice thickness sharply. As the ice thickness rises, the difference between the maximum and average ice resistance increases also. It denotes that maximum ice resistance is more sensitive to ice thickness.

Conclusions
In this paper, an icebreaker moving in level ice has been simulated by using the one-way CFD-DEM coupling method. Numerical results are compared with experimental data and good agreement has been achieved. On this basis, influences of various parameters on ice resistance and ice breaking conditions are further investigated. The following conclusions are drawn preliminarily: (1) One-way CFD-DEM coupling method is firstly used to simulate a ship moving in level ice. To solve the challenge of building level ice on the water domain, this paper puts forward two numerical techniques, which are crucial for the success of the modeling. (2) For a ship moving in level ice, water resistance is an order of magnitude smaller than ice resistance. However, this does not mean water can be easily ignored. Actually, ice resistance is affected by the fluid motion in the OWC method. In other words, the ice resistance in OWC method has included the influence of fluid already.
(3) Ice resistance shows strong oscillation in the process of ice breaking, presenting cycling of 'gradual rising-reaching maximum-sudden drop' generally. This is closely related to the ice failure process of 'elastic deformation-reaching stress, limit-crack generation, and unloading'. After the ship bow fully enters into the level ice region, the average ice resistance reaches stable broadly. Both ice resistance and ice destruction are affected by ship speed and/or ice thickness significantly.

Conclusions
In this paper, an icebreaker moving in level ice has been simulated by using the one-way CFD-DEM coupling method. Numerical results are compared with experimental data and good agreement has been achieved. On this basis, influences of various parameters on ice resistance and ice breaking conditions are further investigated. The following conclusions are drawn preliminarily: (1) One-way CFD-DEM coupling method is firstly used to simulate a ship moving in level ice.
To solve the challenge of building level ice on the water domain, this paper puts forward two numerical techniques, which are crucial for the success of the modeling. (2) For a ship moving in level ice, water resistance is an order of magnitude smaller than ice resistance.
However, this does not mean water can be easily ignored. Actually, ice resistance is affected by the fluid motion in the OWC method. In other words, the ice resistance in OWC method has included the influence of fluid already.
(3) Ice resistance shows strong oscillation in the process of ice breaking, presenting cycling of 'gradual rising-reaching maximum-sudden drop' generally. This is closely related to the ice failure process of 'elastic deformation-reaching stress, limit-crack generation, and unloading'. After the ship bow fully enters into the level ice region, the average ice resistance reaches stable broadly. Both ice resistance and ice destruction are affected by ship speed and/or ice thickness significantly.
In future work, TWC method will be studied further and used in the ship moving in level ice. Furthermore, the variation speed of the ship will be studied, as well as the motion of the ship in six degree of freedom.