GPU-Based DEM Simulations of Global Ice Resistance on Ship Hull During Navigation in Level Ice

The ice resistance on a ship hull affects the safety of the hull structure and the ship maneuvering performance in ice-covered regions. In this paper, the discrete element method (DEM) is adopted to simulate the interaction between level ice and ship hull. The level ice is modeled with 3D bonded spherical elements considering the buoyancy and drag force of the water. The parallel bonding approach and the de-bonding criterion are adopted to model the freezing and breakage of level ice. The ship hull is constructed with rigid triangle elements. To improve computational efficiency, the GPU-based parallel computational algorithm was developed for the DEM simulations. During the interaction between the ship hull and level ice, the ice cover is broken into small blocks when the inter-particle stress approaches the bonding strength. The global ice resistance on the hull is calculated through the contacts between ice elements and hull elements during the navigation process. The influences of the ice thickness and navigation speed on the dynamic ice force are analyzed considering the breakage mechanism of ice cover. The Lindqvist and Riska formulas for the determination of ice resistance on ship hull are employed to validate the DEM simulation. The comparison of results of DEM, Lindqvist, and Riska formula show that the DEM result is between those the Lindqvist formula and Riska formula. Therefore the proposed DEM is an effective approach to determine the ice resistance on the ship hull. This work can be aided in the hull structure design and the navigation operation in ice-covered fields.


Introduction
With the increasing scientific investigations, oil and gas exploitations as well as marine transportations in polar regions, more attention has been paid to shipping performance in ice-covered waters. Icebreaking ships need the ability to maintain a reasonable speed in level ice of certain ice thickness, and the turning ability under specified ice conditions (Kim et al., 2005). The ice load, an integrated effect over the hull area of the ship, is the dominant factor that governs the ship maneuvering performance. The interaction between ship hull and sea ice cover is a complex process, and the ice load is affected severely by many factors, such as ice conditions (level ice, pancake ice, ice ridge, and iceberg), the hull geometry, and the relative velocity between the ship and ice. As the level ice is fractured by the ship's gravity, the flexural failure of ice cover occurs when the bending stress is beyond the ice strength at a distance from the ship. This process depends on the ice thickness and the shipping speed (Valanto, 2001;Kotovirta et al., 2011).
Studies on ship maneuverability in ice have been conducted by field measurements and laboratory experiments. Riska et al. (2001) observed the ice performance of the Swedish Multi-Purpose icebreaker for Viking II. Kotisalo and Kujala (1999) and Hänninen (2003) measured the ice load onboard MT Uikku during the ARCDEV Voyage. A series of tests were also performed with a laboratory-scale model ship to simulate the effect of ice load parameters on an icebreaking ship (Zhou et al., 2013). Von Bock und Polach and Ehler (2011) carried out model tests in both constrained and free modes. The impact of heave and pitch mo-tions on peak values of the ice resistance and the icebreaking pattern were discussed. However, the systems required field measurements and model tests are always very complex and expensive as well as the long study period. Therefore, numerical simulations have been applied widely to determine the global ice resistance on ship hull under different ice conditions and shipping performances. The finite element model (FEM) has been applied to model the ice pressure and global ice resistance on ship hull of the level ice, ice ridge, and iceberg (Aksnes, 2010;Su et al., 2011;Sayed and Kubat, 2011;Zhou et al., 2012). Concerning the break-up feature of ice during the ice-hull interaction, the discrete element method (DEM) shows intrinsic advantages and has been applied widely to simulate the ice resistance on the ship hull. Hansen and Løset (1999) used a 2D disk DEM to calculate ice resistance on ship hull in the broken ice field (Hansen and Løset, 1999). Ji et al. (2013b) adopted the 3D disk element to simulate the ice floes and obtained the global ice resistance on the ship hull. The DEM software, DECICE, was successfully applied to simulate the ice breaking and maneuvering ability of icebreakers (Zhan et al., 2010;Lau et al., 2011). In the DEM simulations of ice resistance on the ship hull, the discrete distribution of ice floes and the breakage of level ice can be described following their physical behaviors (Ji and Liu, 2012). In the DEM simulation of level ice, the ice cover can be constructed with bonded spheres. The dynamic ice resistance on the ship hull is generated when the ice cover fractures under the collision of the ship hull. The magnitude and frequency of the ice resistance are dominated by failure modes. The contact model of bonded ice elements is the key factor in DEM simulation of ice resistance on ship hull. In the DEM simulation of the interaction between ice cover and offshore structure, the bonded model and relative parameters have been validated with the flexural and compressive strength tests (Ji et al., 2012(Ji et al., , 2013b. This model is also applied in the numerical simulation of the interaction between the ice cover and the ship hull. The graphic processing unit (GPU), a lightweight parallel computing device, can be used to improve the computational efficiency of the DEM simulation. The GPU-based parallel algorithm for the neighbor list of the contact detection between elements, such as the uniform grid method and the tree-based grid method, can achieve high computational efficiency (Kureck et al., 2019). The acceleration ratio of the GPU-based parallel algorithm of spherical DEM can be larger than 30 as the spherical DEM is suitable for the lightweight computing device (Gan et al., 2020). The computational efficiency of the polyhedral DEM can also be improved with an appropriate data structure, considering the polyhedral element has complex geometric data (Govender et al., 2015). The GPU-based parallel technology is also applied in the DEM simulation of the interaction between sea ice and marine structures. Small scaled elements can be used in the compression and bending test of sea ice for the parameter calibration based on the GPU (Di et al., 2017). Accelerating by the GPU device, the DEM simulations of the ice bending failure and the ice load during the interaction between sea ice and conical structures illustrate good agreements with the field observation data (Long et al., 2020). Hence, the GPU-based parallel technology is an effective approach to accelerate the DEM simulation of the ice load on marine structures.
In this paper, the discrete element method is developed with a GPU-based parallel algorithm to simulate the interaction between level ice and ship hull. The global ice resistance on the ship hull is calculated with the collision between the ice element and the hull structure. The influences of ice thickness and ship speed on global ice resistance are compared with those obtained from the Lindqvist and Riska formula.

Discrete element method for level ice
In the DEM simulation of the interaction between level ice and ship hull, the ice cover is constructed with 3D bonded spherical elements considering the buoyancy, gravity, and water drag force. For the boundary conditions of the computational domain, a series of springs are set on each boundary element in the vertical direction. The boundary elements are fixed on a rigid wall in the horizontal direction considering the resistance of far-field ice cover under natural conditions, as shown in Fig. 1. Here, a small area of ice cover is zoomed out to show the packing pattern of ice elements.  All of the ice elements are bonded together to model the mechanical behavior of level ice. The parallel bonding model and linear contact force model are adopted to determine the inter-particle contact forces, as shown in Fig. 2. The parallel-bond disk is set over a circular cross-section between the particles as shown in Fig. 2a, and can transmit both force and moments (Potyondy and Cundall, 2004). Here, X A and X B are the position vectors of elements A and B, , and , are the normal and shear component vectors of force and moment, respectively.
The maximum normal and shear stresses in the bonding section are determined with the inter-particle force and moment, and can be written as: (1) where A, R, J, and I are the area, radius, polar inertia moment, and inertia moment of the bonding disk, respectively. Here , , , where R is the radius of the bonding section and is set as the bonded particle radius here.
In the parallel bonded elements of sea ice, two major criteria are introduced to determine the failure of bonded ice elements, as shown in Fig. 3. One is the bonding strength σ max and the other is the maximum strain ε max . The bonded disk will be damaged when the inter-particle stress σ approaches the bonding strength σ max , and will be broken when the inter-particle strain ε is larger than the maximum strain ε max . When , the sea ice performs as strain soften behavior. Here is the inter-particle strain when the stress equals the bonding strength σ max, i.e., . This failure model had been applied in the breakage of rock material with DEM simulation under uniaxial compression (Ji and Di, 2013). In the traditional parallel bonding model, the strain soften phase is ignored when . Here we set to consider the strain soften behavior of sea ice.
With this strain soften model, the elastic potential energy of bonded ice elements can be dissipated to avoid the elements splash when the bonded elements are broken.
The contact force between two ice elements is calculated with an elastic-viscous contact model based on the Mohr-Coulomb shear friction law, as shown in Fig. 2b, where M A and M B are the mass of ice elements A and B, K n and K s are the normal and tangential stiffness, C n and C s are the normal and tangential damping coefficients, and μ is the inter-particle friction coefficient. The normal and the tangential forces are calculated as: where , and , are the relative displacement and velocity of the two contacting particles in normal and tangential directions, respectively.
is the unit vector in the tangential direction. The damping coefficient , the dimensionless normal damping coefficient , where e is the coefficient of restitution. The particle stiffness K n is proportional to the elastic modulus of sea ice, E, and can be written as K n = πDE/4, where D is the particle diameter. The normal and tangential stiffnesses have the relationship of K s = 0.5K n and C s = 0.5C n (Ji et al., 2012).
In the linear viscous-elastic contact force model, the binary contact time of two free elements can be calculated as: where M is the effective mass to two collision elements. This binary contact time is determined by the element size and material properties. In this DEM simulation, the computational time step is set at 1/30 of the binary contact time (Ji et al., 2013a, 2013c).
2.2 Structure construction of ship hull A cargo ship, Tian You ship owned by COSCO Shipping, is used in the DEM simulation, as shown in Fig. 4a. The length and width in the water level of Tian You ship are 186.4 m and 28.5 m, and the draught is 11.0 m. The ship  hull is constructed with triangular elements without considering the deformation of the hull structure. Here, the ship hull is generated with 13 582 triangular elements and 7 282 nodes, as shown in Fig. 4b. The interaction between ice cover and ship hull is divided into contacts between triangular hull elements and ice elements. The contact modes between ice spheres and hull triangular mainly include three different patterns, i.e., sphere-corner contact, sphere-edge contact, and sphere-plane contact, as shown in Fig. 5. The contact force model between ice particle and hull element is similar to that of inter-particles of sea ice. Here, the hull element is simplified as a rigid surface. A contact detection scheme is presented for hull boundary surfaces and ice spheres based on a series of vector projections to reduce the problem dimensionally. 2.3 GPU-based parallel algorithm for DEM simulations of ice resistance on ship hull In DEM simulations, the neighbor list and contact detec-tion of particles are the most cumbersome part which seriously affects the computational efficiency. According to the characteristics of GPU-based parallel technology and DEM, the cell-list method is employed to detect contact between spherical particles (Di et al., 2017). The simulation domain is divided into grids whose size is slightly larger than the particle diameter so that each element would only contact the particles in neighbor cells. Here the cell-list is implemented within the 3×3×3 cells which enclose the objective element. This method can achieve a high speed-up ratio on a multi-core processor with shared memory. Fig. 6 shows the cell-list method for contact detection between particles. According to the space and grid of each element in Fig. 6a, the list of new particle ID can be established by sorting the cell ID, as shown in Fig. 6b. Nishiura and Sakaguchi (2011) used the shared memory to determine the first and the last element ID in each cell, i.e., the cell ID which each element belongs to, and every element ID in each cell can be obtained. Accordingly, the contact detection between adjacent particles in the neighbor list can be implemented. The parallel algorithm of this method on a GPU device shows the high computational efficiency of the DEM simulation. 3 Ice resistance on ship hull simulated with DEM

Simulation setup
In the DEM simulation of the interaction between level ice and ship hull, the ship speed is constant. Here we ignore the ice-induced vibration and deformation of hull structures, and the ice cover is static. The main computational parameters are listed in Table 1. The inter-element bonding strength of level ice is a key parameter to affect the ice resistance on the ship hull. Based on the previous investigation, the flex- HU Bing et al. China Ocean Eng., 2021, Vol. 35, No. 2, P. 228-237 231 ural strength of sea ice is determined as σ f = 0.7 MPa in this study. The field and laboratory tests showed that the compressive, flexural and shear strengths of sea ice are the function of ice temperature, salinity and stress rate (Timco et al., 2010;Ji et al., 2011). The compressive and flexural strengths of sea ice have been simulated with DEM in the investigations of the interaction between sea ice and offshore structures (Ji et al., 2013c;Di et al., 2017). In the DEM simulation, the interactions between ship hull and ice elements are characterized as forces (unit: N). Therefore, the ice resistance can be obtained as the resultant force on the ship hull from ice elements on each time step. The size effect in the parallel bonding model based DEM is an important issue for the breakage simulation of sea ice. The main problem is how to establish the relationship of the macro-micro parameters (Di et al., 2017), where the macro parameters are the strengths of sea ice materials such as flexural strength or compressive strength, and the micro parameters are the bonding strength, element size, etc., in the DEM. The dimensionless bending or compressive test can be simulated to establish the relationship of the macro-micro parameters (Long et al., 2019), so that the bonding strength can be determined with the ice bending strength, and the element size which can be characterized as the layer number of the sea ice. In the following simulations, the element size (at least 2 layers) and the bonding strength may be different, but the ice bending strength keeps the same based on the relationship by Long et al. (2019). The inter-element bonding strength has been calibrated with the strength of sea ice obtained with physical experiments. By considering the particle size effect in the DEM, the interelement bonding strength is determined according to different particle sizes. The specific element number used in each case will be introduced. Fig. 7 shows the simulated ice cover-hull interaction process at time t = 200 s and t = 515 s, respectively. A water channel is generated behind the ship. The global ice resistance on the ship is determined by summing the contact forces between ice element and hull in different directions. Fig. 8 plots the ice resistance histories in the three directions on the ship hull during the ship running through the entire ice cover. Here the ship navigating direction is the xdirection, the vertical direction is the z-direction, and the other one is the y-direction. Normally the ice resistance in the x-direction is mainly focused. When the ship navigates from the open water into the level ice, the ice resistance in the x direction increases from zero, then reach steady-state with a certain fluctuation accompanying the breakage of ice cover. The global ice resistance on the ship hull mainly includes two portions: one is bending and crushing from the bow, and the other is from the friction of two broadsides. During the shipping from open water into ice, the ice resistance on the bow generates and increases along the crushing length into the ice cover. Meanwhile, the broadsides frictions increase obviously along the crushing length. The mean ice resistance is determined as the mean value of the peak loads (Liu and Ji, 2018). The maximum and mean ice resistances are 30.47 MN and 15.58 MN in the x-direction. In the y-and z-direction, the ice resistance always fluctuates around zero during the ship running through ice cover. The sea ice interacts with the two sides of the ship hull, which causes the fluctuation of the ice resistance in the ydirection. But it is not usual of the fluctuation around zero of the ice resistance in the z-direction, although the magnitude is smaller than that in the x and y directions. The bending failure is the main failure mode of sea ice interact- ing with the ice-going ships, especially the icebreaker, owing to the large inclined geometric structure at the bow. Thus, the global ice resistance of ship hull in the z-direction will be vertically upward, i.e., positive in this instance, for ice-going ships. However, the ship hull used in this study does not have a large inclined geometric structure at the bow. The bow structure is generally straight in fact. Hence, the friction between the ship and sea ice can be vertically downward or upward, which leads to the fluctuation around zero of the ice resistance in the z-direction. The failure mode of sea ice directly influences the ice resistance on the ship hull. The snapshot of the DEM simulation at t = 200 s is shown in Fig. 9 for the analysis of shipice interaction. The color represents the speed variation. We can find that the ice cover is broken into a couple of large blocks in front of the bow. In the DEM simulation of iceship interaction, the breakage of ice cover dominates the ice resistance period and magnitude. In the analysis of bending failure during the interaction between sea ice and inclined structure, as the broken ice size is larger, the ice resistance becomes larger. With the analysis of the broken size of ice cover, the ice resistance characteristics can also be determined. However, the geometric shape of the ship bow used in this study is sharp in the horizontal plane. The ice cover is breaking like paper cutting during interacting with the ship bow. Meanwhile, the geometric shape of the ship bow is concave for the inner side of the hull but nearly straight in the vertical direction. Thus, the sea ice appears bending failure at the two sides of the ship and shapes a large broken size. Besides, the crashing failure of sea ice also occurs because of the vertical shape of the ship bow, i.e., the sea ice is broken into small pieces.

Simulation results
In most circumstances, the cargo ship which is not an ice-class ship would never navigate in level ice. In fact, an icebreaker will be deployed to escort the cargo ship. However, the ice resistance on the cargo ship in level ice is still required in specific circumstances. For example, it is required in the evaluation of the load reduction rate considering the cost. Besides, the stem angle is an important parameter for the ice load and icebreaking performance. Normally, the cargo ship has a larger stem angle, but the icebreaker has a small stem angle. The ice resistance on the cargo ship in level ice can be used to investigate the change rate of the ice resistance to the stem angle. Therefore, it is still necessary to study the ice resistance on the cargo ship in level ice.

Influences of ice thickness and navigation speed on global ice resistance
The ice resistance on the ship hull is dominated by ice conditions, ship structure, and shipping performance. The breaking lengths of ice pieces also depend on the ice thickness and shipping speed. The influences of ice thickness and shipping speed are analyzed on ice resistance here. The simulation results are validated with empirical formulas, including the Lindqvist and Riska formulas. Lindqvist (1989) proposed a formula for global ice resistance on ship hull based on field measurement and theoretical analysis (Su et al., 2010). Based on the fundamental idea of the Lindqvist formula, Riska et al. (1997) proposed a correction formula with the field monitoring data in the  HU Bing et al. China Ocean Eng., 2021, Vol. 35, No. 2, P. 228-237 233 Baltic Sea (Tan et al., 2013). The two formulas are both important approaches to quickly evaluate the ice resistance on ships navigating in level ice (Hu and Zhou, 2016).

Empirical formulas: Lindqvist and Riska
The ice resistance on ships defined by the Lindqvist formula, , can be expressed as: ; where , , , and are total ice resistance, crushing resistance, bending resistance and resistance due to submersion, respectively; , , and are the stem angle, waterline entrance angle and flare angle respectively, and ; is the gravity acceleration; , , and are the ship waterline length, breadth, and draught, respectively. The ice resistance on ships defined by the Riska formula, , can be expressed as: where and are the lengths of the bow and parallel sides section, respectively; is the stem angle in degrees; , , , , , , and are all constant parameters that can refer to Hu and Zhou (2016). The parameters of the ship hull for the Lindqvist and Riska formulas are listed in Table 2. Other parameters can refer to Table 1. 4.2 Influence of ice thickness on global ice resistance The element sizes used in each case are different according to different ice thicknesses, so that the element numbers and the computational time are also different, as shown in Table 3. The minimum computational time is 5.05 hours for 324220 elements, and the maximum computational time is 14.92 hours for 728252 elements. To analyze the influence of ice thickness on the ice resistance, the DEM simulation is implemented when the ice thickness H ice varies from 0.5 m to 2.5 m at the constant ship speed V ship = 2.0 m/s. The ice resistance histories in the x-direction under different sea ice thicknesses during t = 245−305 s are plotted in Fig. 10. The ice resistance histories show obvious fluctuation with different ice thicknesses while the amplitude of the ice resistance increases with the increase of the ice thickness. Fig. 11 is the comparison of DEM, Lindqvist, and Riska formula results with different ice thicknesses. The comparison shows that all the ice resistance increases in a near-linear trend. It is reasonable that the  Lindqvist formula result is larger than the Riska formula result because the Riska formula considers the specific field data in the Baltic Sea. The maximum ice resistance of the DEM is larger than that of the Lindqvist formula, while the mean ice resistance of the DEM is quite close to that of the Lindqvist formula. This comparison indicates that the DEM simulation of the ice resistance on the cargo ship can obtain reasonable results.

Influence of navigation speed on global ice resistance
The ship speed is also a key factor affecting the ice resistance on ship hull through the failure mode of ice cover. Here the ship speed varies from 0.5 m/s to 4.0 m/s to analyze the influence of navigation speed. Because the ice thickness is the same, the element size and element number are also the same. The computational time for each case is shown in Fig. 12. The computational time grows with the decrease of the ship speed, because the decrease of the ship speed caused the growth of the simulation time.  Fig. 13 shows the ice resistance histories versus navigating distance at different ship speeds. The ice resistance histories at different ship speeds fluctuate and illustrate impulse shape. When the ship speed is large, such as V ship = 3.0 m/s and 4.0 m/s, large amplitudes appear in the ice resistance history. With the maximum and mean ice resistances, the Lindqvist and Riska formula results are shown in Fig. 14. It shows that all of them increase with the increas of the ship speed. The growth of ice resistance is not sensitive compared with the factor of sea ice thickness, especially when the ship speed is high. This is reasonable because the de-bonding criterion between elements including the bonding strength and stiffness is not related to the loading rate. Hence, the growth of the ice resistance is mainly caused by the increasing damping force related to the ship speed. Normally, the ice strength increases with the increase of the loading rate (Timco et al., 2010). The effect of ship speed will be more obvious if the influence of the stress rate on the sea ice strength is considered (Ji et al., 2011). On the other hand, the maximum ice resistance of the DEM result is the largest, while the mean ice resistance of the DEM result is quite close but a bit smaller than the Lindqvist formula result.

Discussions
The comparison between numerical results and empirical formula results obtains good agreements. But there is still a difference. We believe that the main reasons that caused the difference include the following reasons, although it is difficult to be completely explained.
(1) The Lindqvist formula is based on several ideal models, while the Riska formula is a correction of the Lindqvist formula based on the field data in the Baltic Sea. These two formulas have a lot of differences with the proposed DEM in the ice model and ice-ship interaction process.
(2) The empirical formulas calculate the static force on the ship, which neglects the dynamics of the ice-ship interaction.
(3) The empirical mainly considers the first breaking of level ice during the ice-ship interactions but neglects that the broken fragments may be compressed closely on the ship hull and broken again. It can be even more complex.

Conclusions
The sphere-based DEM with the parallel bonding model is built to simulate the breakage of the level ice. The GPU-based parallel algorithm is developed to accelerate the computational efficiency. Accordingly, the ice cover is modelled with spherical elements considering the buoyancy  and gravity, while the hull structure of a cargo ship is constructed with triangle elements. This proposed approach is used in the analysis of the interaction between the ship hull and the level ice. The ice resistance on the ship hull, calculated as the resultant force on the ship hull from spherical elements, is analyzed and validated with the Lindqvist and Riska formula with different ice thicknesses and ship speeds. The comparisons of both the maximum and mean ice forces to the empirical formulas show good agreements.
In the next study, the six-degree motion of the ship hull will be considered to investigate the icebreaking performance. The distribution of the ice pressure on the ship hull will be obtained for the strength analysis of the ice-going ship. Besides the level ice, ice ridge and rubble ice will be considered in the analysis of ice loads on ship hulls.