1. Introduction

2. State of the Art

3. Numerical Scheme

3.1. Analysis Process

Our study sought to develop a reasonable numerical model to predict the behavior of drift ice and estimate the wear caused by the collision between drift ice and a vessel’s hull. The environmental load of the Arctic route can be modeled with the CFD scheme, after which the computed load is transferred to the ice floes. The motion of many ice particles colliding on a ship can be modeled with the DEM scheme. The impact computed with the DEM algorithm can be used to construct models to describe the scattering of ice floes approaching an object, as well as to assess the contact forces on the hull surface. Afterward, the wear magnitude was evaluated based on the collision mechanism between the hull surface and the ice floe, and was simulated with three-dimensional DEM. In this study, we sought to predict hull wear caused by a set of ice floe collisions by taking advantage of the unique strengths of the CFD and DEM methods. The flow-induced load was computed using the continuum method and the behavior of drift ice was simulated with the discrete method, whereas the continuum and discrete methods were implemented using the Ansys Fluent code and Rocky-DEM code, respectively. Ansys Fluent is a VOF (volume of fluid)-based CFD code that simulates the hydrodynamic load of seawater on ice, whereas the Rocky-DEM is a particle behavior analysis solver that can accurately predict the load and behavior of particles colliding with other particles, as well as particles colliding with hulls. From the perspective of particle geometry, the Rocky-DEM approach can be used to model the free shape of the particles whereas other DEM schemes (Cleary et al., [34]; Morrison and Cleary [35]) exclusively assume that the particles are spherical, resulting in unrealistic results.

In the coupling analysis method, hydrodynamic loads such as drag force and buoyancy caused by seawater are calculated through CFD analysis and then transferred to the DEM analysis. In the DEM analysis, the load calculated from the CFD analysis is applied to evaluate the behavior and collision load of the particles. In this approach, the fluid phase and particle phase are analyzed through each specialized solver, which greatly enhances the accuracy of the analysis. Coupling analysis can be classified as one-way and two-way according to the interrelationships between the fluid phase and the particle phase. One-way coupling analysis is suitable when particles are exposed to the flow environment but the particle behavior does not have a substantial influence on the flow. In contrast, two-way coupling analysis may be suitable when the behavior and flow of particles closely affect each other. The flow of the route acts as the dominant load on the drift ice. However, the behavior of the drift ice can slightly affect the flow of the route but has no effect on the overall flow. Therefore, the one-way coupling analysis is reasonably well suited as a basis for the wear evaluation model due to the collision between the drift ice and the hull.

Figure 1. Schematic diagram of numerical codes to simulate the CFD – DEM coupling behavior.

Figure 1. Schematic diagram of numerical codes to simulate the CFD – DEM coupling behavior.

3.2. Theoretical Background of DEM

DEM is a numerical method that is often used to compute the motion and behavior of a large number of particles. In the DEM analysis, the load applied to the particle is defined as the combination of the contact load (surface force) caused by the collision and the non-contact load (body force). Non-contact load can be defined as a load that affects particles regardless of whether collision occurs, such as gravity and ambient fluid forces. The inertia force due to particle motion and the effects of the marine environment are applied as a non-contact load to the center of gravity of the particles. The speed and direction of motion are changed by the collision of particles, and some kinetic energy is dissipated. DEM analysis can be defined as the repeated process of calculating particle behavior (body force) and collision (surface force) during the analysis [36].

Although some particles are modeled as flexible in the DEM algorithm, particles are generally modeled as rigid objects to make the analysis more efficient. When a particle collides with other particles and structures, deformation occurs due to the collision energy. Particles modeled as rigid objects do not deform even when a collision occurs. The deformation caused by the collision was replaced by an overlap of particles because the DEM algorithm generally uses rigid particles. The relationship between load and displacement is defined as a function of either $F=kx$ or $F=k{s}_n$ when describing the behavior of a flexible body or a rigid particle, respectively. Where F denotes force, $k$ is stiffness, $x$ is deformation, and $s_n$ represents the overlap distance.

Figure 2. Comparison of the behavior of flexible and rigid particles.

Figure 2. Comparison of the behavior of flexible and rigid particles.

Collision load is defined as the component of normal force and tangential force, and each load is defined as a function of stiffness and overlap. In this study, the normal force was defined using the hysteretic linear spring model (HLS). The HLS model is an elastoplastic model that is used by separately defining the loading stiffness, in which the direction of motion is maintained in the direction of the initial collision, and the un-loading stiffness, where the direction of motion is reversed after the initial collision. In Figure 3, the slope of $\overline{AB}$ is the loading stiffness, whereas $\overline{BC}$ represents the un-loading stiffness. $\overline{CA}$, which is an overlap that remains after being completely numerically un-loaded, is plastic deformation, and the area of the triangle is the energy lost due to noise, heat, and other factors after collision.

The linear spring Coulomb limit (LSCL) model was used to evaluate the tangential force that directly affects the wear. The Coulomb friction model, which defines the tangential force as the product of the friction coefficient ($\mu$) and the normal force ($F_n$), is widely used to calculate the tangential forces between contact surfaces. However, unlike the LSCL model, it cannot account for contact area and sliding stiffness (Figure 5). Therefore, the tangential force calculated through the Coulomb friction model was applied as a limitation of the LSCL model. The normal force model and tangential force model are shown in Figure 4 and Figure 5, respectively.

In Figure 4, the loading stiffness $K_{nl}$ and the unloading stiffness $K_{nu}$ are respectively defined as follows:

$$\frac{1}{K_{nl}}=\frac{1}{K_{nl,p1}}+\frac{1}{K_{nl,p2}}$$

(1)

$$K_{nu}=\frac{K_{nl}}{{\epsilon }^2}$$

(2)

$$K_{nl,p}=E_{p}L$$

(3)

where $K_{nl,p1}$ and $K_{nl,p2}$ represent the stiffness of the particle, L is the particle size, and $E_p$ represents the bulk Young's modulus of the particles. Moreover, ε denotes the restitution coefficient, which is not an individual property of particles but a relative property of a collision pair.

where $K_{\tau }$ represents the sliding stiffness and is defined as the product of $r_K$ (tangential stiffness ratio) and $K_{nl}$, as shown in the following equation:

$$K_T=r_K{K}_{nl}$$

(4)

4. Evaluation of Hull Wear Due to Collision with Ice Floes

4.1. Evaluation Conditions

In the early stages of operation, the painted hull surface is damaged by friction with the floating ice. Since the shape of the hull does not change, there is no need to account for the shape change due to wear. If the operation is prolonged, friction with the ice floe accumulates and the material of the hull is lost. Since the shape of the hull changes due to material loss, the shape change due to wear must now be accounted for. Evaluation models suitable for each situation were thus presented. The evaluation was conducted assuming a scenario in which a vessel sails 200 m on a route with ice floes. To calculate the fluid force of the route, an evaluation model was constructed with Ansys Fluent. It was assumed that there was no seawater flow, and the buoyancy force was considered by implementing the pressure gradient of the fluid force. To simulate buoyancy and conditions in which there is no fluid flow, an arbitrary speed was applied upward from the lower part of the CFD model and analyzed, after which the speeds in all directions were patched to 0. The velocity in all directions was fixed to 0 to prevent particle movement due to seawater flow, and the depth direction pressure gradient was implemented to define the buoyancy force. The results of CFD analysis were applied as the load of DEM analysis. For the coupling of CFD and DEM, the Ganser drag law was applied to calculate drag force. The density of water was defined as 1000 kg/m³ and the viscosity as 0.001003 kg/m^-s.

Transport Canada [42], a Canadian company that operates through the Arctic Ocean, suggests a range of safe operating speeds to avoid potential accidents from ice collision. Therefore, our evaluations were carried out based on the safe speeds of 4, 6, and 10 knots suggested by Transport Canada (AMNS) [42].

According to the Arctic Ice Regime Shipping System (AIRSS) [43] developed by Transport Canada, one of the main factors defining the ice regime that affects ships operating in polar regions is ice concentration, which is defined according to the percentage of space occupied by ice in the route. In this study, ice concentrations of 60% and 80% were evaluated, as illustrated in Figure 7. Ice concentration was defined as the area ratio of particles in the evaluation area. The method for the evaluation of ice concentration is shown in Figure 8.

The shape of the ice particles was defined as described by Zhang et al. [23]. The ice particle modeling results are shown in Figure 9. The size distribution of the ice particles was determined by Liu and Ji [22]: 25% of 1.0 m particles, 50% of particles larger than 1 m and smaller than 1.5 m, and 25% of particles larger than 1.5 m and smaller than 2.0 m. This size distribution was used for the three particle shapes. The sizes of the three different particle shapes were defined as the sizes of spherical particles having the same volume. Figure 9 illustrates the ice particle shape modeling.

An oil tanker with a length of 43 m, a width of 6.5 m, a depth of 3.6 m, and a draft of 2.4 m was evaluated. The material properties applied to the analysis are shown in Table 1 and Table 2.

Figure 10. Ship modeling.

Figure 10. Ship modeling.

The DEM analysis conditions reflecting the above conditions are shown in Figure 11. The evaluation models were constructed by initially placing particles according to the ice concentration on the surface of the 200 m route and then moving the hull. The blue area in Figure 11 is the seawater area defined by the CFD analysis result. To perform the evaluation in the same ice arrangement for each ice concentration, the particle arrangement was stored and used as an initial condition for each evaluation case. The initial arrangements of drift ice according to ice concentration are shown in Figure 12 and Figure 13.

4.2. Coating Material Wear Assessment

The evaluations were performed assuming a scenario in which a coating material such as paint was separated from the hull by friction between the ice and the hull. Evaluations were performed at ice concentration of 60% and 80% and ship speeds of 4, 6, and 10 knots.

Given that the coating material is very thin compared to the thickness of the hull, the effect of material loss on the flow and behavior of ice is negligible. Therefore, the analysis was performed without reflecting the shape change, and the wear was evaluated based on the shear work according to Archard's wear law. Given that wear is a type of damage that occurs over a long period, it was not evaluated as a result of a specific point in time but was evaluated as an accumulated value on the hull after passing a 200 m route. Shear work, tangential force, and sliding distance were evaluated in the form of cumulative values for each evaluation condition according to the Archard wear law. The hull was divided into three regions to conduct the evaluations: forward (FWD; i.e., the bow), midship, and after (AFT; i.e., the stern) of port.

Figure 14. Evaluation area.

Figure 14. Evaluation area.

Six analysis scenarios were simulated by combining three cases of ship speed and two cases of ice concentration. Among the six simulated scenarios, Figure 15 and Figure 16 show the simulation results of the scenario with the slowest ship speed and the lowest ice concentration (4 knots, 60%), whereas Figure 17 and Figure 18 illustrate the scenario with the fastest ship speed and the highest ice concentration (10 knots, 80%).

The average, maximum, and total cumulative shear work causing wear on the FWD, midship, and AFT are shown in Figure 19, Figure 20 and Figure 21.

The average, maximum, and total cumulative tangential force, which are components of shear work, at the FWD, midship, and AFT are shown in Figure 22, Figure 23 and Figure 24.

The average, maximum, and total cumulative sliding distance, which are components of shear work, at the FWD, midship, and AFT are shown in Figure 25, Figure 26 and Figure 27. The total area affected by wear is shown in Figure 28.

Our findings indicated that the accumulated total shear work, tangential force, and sliding distance in all areas at 80% ice concentration was greater than 60%. As the ice concentration increases, the cumulative results of each item increase because the number of collisions also increases.

Our results also demonstrated that as the ship speed increases, the cumulative shear work, tangential force, sliding distance of the drift ice per unit time, and wear affected area at the FWD also increased. These results indicate that the wear occurs in a deep and wide area as the ship speed and the ice concentration increases. Since the direct collision between drift ice and the hull is dominant in the FWD, the tangential force and sliding distance generated in each collision increases as the ship speed increases. Therefore, the cumulative average, maximum, and total amounts of each evaluation item tended to increase.

In the midship and AFT, as the ship speed increases, the cumulative shear work, tangential force, sliding distance, and the area affected by wear tend to decrease inversely. This is because the drift ice, which first collided with FWD, rubs against the hull while flowing along the outer wall. When the ship speed is relatively slow, the drift ice that first collided with FWD moves smoothly along the outer wall of the hull. However, as the ship speed increases, it tends to move away from the hull rather than flowing along the outer wall.

This trend is clearly shown in Table 3, which summarizes the total cumulative shear work per unit area as a ratio of FWD. Given that the total volume of wear according to the Archard's wear law is proportional to the shear work, the results are summarized based on the shear work. The total shear work accumulated per unit area of the midship and AFT relative to the FWD region was evaluated to be as high as 4.99%, and the ratio tends to decrease as the speed increases. As expected, our results demonstrated that FWD directly collides with the drift ice, and therefore the amount of wear and the ship speed are directly related. In contrast, midship and AFT wear was not related to ship speed because the particles that first collided with the FWD induced wear while moving along the outer wall of the hull. Therefore, to prevent abrasion in the FWD, the ship speed must be reduced, or sufficient reinforcement must be provided according to the ship speed. Moreover, to prevent midship and AFT wear, the outer shape of the hull could be redesigned to minimize the damage caused the drift ice that first collides with the FWD.

Models without shape changes due to wear-induced material loss can also evaluate the wear depth by applying the accumulated shear work to the Archard wear law. As shown in Equation (27), the wear depth was evaluated by multiplying the accumulated shear work by the volume loss ratio C according to the Archard wear law. The wear depth was evaluated by dividing the evaluated wear volume by the individual cell area of the hull surface. The volume loss ratio C per unit shear work was defined as 5e-7m³/J, which is the same as the value applied in the following section. The wear depth according to ship speed and ice concentration for each hull position was evaluated as shown in Table 4.

$$\mathrm{W}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{D}\mathrm{e}\mathrm{p}\mathrm{t}\mathrm{h}=\frac{C\ast CumulativeShearWork}{CellArea}$$

(27)

5. Conclusions

References

1. Ardhuin, F.; Sutherland, P.; Doble, M.; Wadhams, P. Ocean waves across the Arctic: Attenuation due to dissipation dominates over scattering for periods longer than 19s. Geophys. Res. Lett. 2016, 43, 5775–5783. [Google Scholar] [CrossRef]
2. Riska, K.; Bridges, R. Limit state design and methodologies in ice class rules for ships and standards for Arctic offshore structures. Mar. Struct. 2019, 63, 462–479. [Google Scholar] [CrossRef]
3. Adumene, S.; Ikue-John, H. Offshore system safety and operational challenges in harsh Arctic operations. J. Saf. Sci. Resil. 2022, 3, 153–168. [Google Scholar] [CrossRef]
4. International Association of Classification Societies (IACS). Polar ship rules, URI. 2016.
5. International Association of Classification Societies (IACS). Technical background (TB) documents for URs concerning polar class (UR I). 2016.
6. International Maritime Organization (IMO). Descriptions of polar classes, submitted by the international association of classification Societies (IACS), development of a mandatory code for ships operating in polar waters, DE 55/12/10. 2011.
7. International Maritime Organization (IMO). International code for ships operating in polar waters (polar code), resolution MEPC.264(68). 2014. adopted on. 15 May.
8. Gao, Y.; Hu, Z.; Wang, J. Sensitivity analysis for iceberg geometry shape in ship–iceberg collision in view of different material models. Math. Probl. Eng. 2014, 2014, 11. [Google Scholar] [CrossRef]
9. Kim, H.; Daley, C.; Colbourne, B. A numerical model for ice crushing on concave surfaces. Ocean. Eng. 2015, 106, 289–297. [Google Scholar] [CrossRef]
10. Liu, L.; Sun, S.; Ji, S. Interaction between Floater and Sea Ice Simulated with Dilated Polyhedral DEM, DEM 2016: Proceedings of the 7th International Conference on Discrete Element Methods, 1-. 4 August.
11. van den Berg, M. Discrete Numerical Modelling of the Interaction Between Broken Ice Fields and Structures, Norwegian University of Science and Technology, 2019.
12. Sun, S.; Shen, H. Simulation of pancake ice load in a circular cylinder in a wave and current field. Cold Reg. Sci. Technol. 2012, 78, 31–39. [Google Scholar] [CrossRef]
13. Kwon, J.; Jeon, B.; Kim, J.; Bo Wang, H.Y.; Basu, R.; Lee, H.; Kendrick, A. Structural Integrity Assessment of Cargo Containment Systems in Arctic LNG Carriers under Ice Loads. In Fourth Annual Arctic Shipping St. Petersburg, Russia.
14. Nho, I.S.; Lim, S.J.; Kang, K.J. Safe Speed Estimation of Arctic Ships considering Structural Safety. J. Soc. Nav. Archit. Korea 2018, 55, 236–242. [Google Scholar] [CrossRef]
15. Han, Y.; Sawamura, J. Fatigue Damage Calculation for Ship Hulls Operating in Pack Ice, Proceeding of the 24th International Conference on Port and Ocean Engineering under Arctic Conditions, 11–16 June.
16. Cai, W.; Zhu, L.; Yu, T.X.; Li, Y. Numerical simulations for plates under ice impact based on a concrete constitutive ice model. Int. J. Impact Eng. 2020, 143, 103594. [Google Scholar] [CrossRef]
17. Zhu, L.; Qiu, X.; Chen, M.; Yu, T.X. Simplified ship-ice collision numerical simulations, In The 26th International Ocean and Polar Engineering Conference, Rhodes, Greece, 16. 20 June.
18. Kim, M.; Lee, S.; Lee, W.; Wang, J. Numerical and experimental investigation of the resistance performance of an icebreaking cargo vessel in pack ice conditions. Int. J. Nav. Archit. Ocean. Eng. 2013, 5, 116–131. [Google Scholar] [CrossRef]
19. Robb, D.M.; Gaskin, S.J.; Marongiu, J.C. SPH-DEM model for free-surface flows containing solids applied to river ice jams. J. Hydraulic Research 2016, 54, 27–40. [Google Scholar] [CrossRef]
20. Huang, L.; Tuhkuri, L.; Igrec, B.; Li, M.; Stagonas, D.; Toffoli, A.; Cardiff, P.; Thomas, G. Ship resistance when operating in floating ice floes: A combined CFD&DEM approach. Mar. Struct. 2020, 74, 102817. [Google Scholar]
21. Luo, W.; Jiang, D.; Wu, T.; Guo, C.; Wang, C.; Deng, R.; Dai, S. Numerical simulation of an ice-strengthened bulk carrier in brash ice channel. Ocean Eng. 2020, 196, 106830. [Google Scholar] [CrossRef]
22. Liu, L.; Ji, S. Dilated-polyhedron-based DEM analysis of the ice resistance on ship hulls in escort operations in level ice. Mar. Struct. 2021, 80, 103092. [Google Scholar] [CrossRef]
23. Zhang, J.; Zhang, Y.; Shang, Y.; Jin, Q.; Zhang, L. CFD-DEM based full-scale ship-ice interaction research under FSICR ice condition in restricted brash ice channel. Cold Reg. Sci. Technol. 2022, 194, 103454. [Google Scholar] [CrossRef]
24. Suyuthi, A.; Leira, B.J.; Riska, K. Fatigue damage of ship hulls due to local ice-induced stresses. Appl. Ocean. Res. 2013, 42, 87–104. [Google Scholar] [CrossRef]
25. Cho, S.; Chun, E.; Yoo, C.; Jeong, S.; Lee, C. The measuring methodology of friction coefficient between ice and ship hull. J. Soc. Nav. Archit. Korea 2011, 48, 363–367. [Google Scholar] [CrossRef]
26. Kietzig, A.M.; Hatzikiriakos, S.G.; Englezos, P. Physics of ice friction. J. Appl. Phys. 2010, 107, 081101. [Google Scholar] [CrossRef]
27. Shimizu, K.; Noguchi, T.; Seitoh, H.; Okada, M.; Matsubara, Y. FEM Analysis of Erosive Wear. Wear 2001, 250, 779–784. [Google Scholar] [CrossRef]
28. Xie, L.-J.; Schmidt, J.; Schmidt, C.; Biesinger, F. 2D FEM estimate of tool wear in turning operation. Wear 2005, 258, 1479–1490. [Google Scholar] [CrossRef]
29. Chen, J.; Wang, Y.; Li, X.; He, R.; Han, S.; Chen, Y. Erosion prediction of liquid-particle two-phase flow in pipeline elbows via CFD–DEM coupling method. Powder Technol. 2015, 282, 25–31. [Google Scholar] [CrossRef]
30. Xu, L.; Zhang, Q.; Zheng, J.; Zhao, Y. Numerical prediction of erosion in elbow based on CFD-DEM simulation. Powder Technol. 2016, 302, 236–246. [Google Scholar] [CrossRef]
31. Zhang, D.M.; Gao, C.P.; Yin, Z.Y. CFD-DEM modeling of seepage erosion around shield tunnels. Tunn. Undergr. Space Technol. 2019, 83, 60–72. [Google Scholar] [CrossRef]
32. Walker, C.I.; Hambe, M. Influence of particle shape on slurry wear of white iron. Wear 2015, 332, 1021–1027. [Google Scholar] [CrossRef]
33. Shunying, J.I.; Zilin, L.I.; Chunhua, L.I.; Shang, J. Discrete element modeling of ice loads on ship hulls in broken ice fields. Acta Oceanol. Sin. 2013, 32, 50–58. [Google Scholar]
34. Cleary, P.W.; Morrisson, R.; Morrell, S. Comparison of DEM and experiment for a scale model SAG mill. Int. J. Miner. Process. 2003, 68, 129–165. [Google Scholar] [CrossRef]
35. Morrison, R.D.; Cleary, P.W. Using DEM to model ore breakage within a pilot scale SAG mill. Miner. Eng. 2004, 17, 1117–1124. [Google Scholar] [CrossRef]
36. ESSS Rocky, Release 2022 R1.1, DEM Technical Manual, ESSS Rocky DEM, S.R.L.
37. Lee, S.; Lee, J.; Hwang, S. Suggestion of Practical Application of Discrete Element Method for Long-Term Wear of Metallic Materials. Appl. Sci. 2022, 12, 10423. [Google Scholar] [CrossRef]
38. Ansys Inc., Fluent Theory Guide Release 2020, Ansys Inc.
39. ESSS Rocky, Release 2022 R1.1, CFD Coupling Technical Manual, ESSS Rocky DEM, S.R.L.
40. Ganser, G.H. A rational approach to drag prediction of spherical and nonspherical particles. Powder Technol. 1993, 77, 143–152. [Google Scholar] [CrossRef]
41. Archard, J.F. Wear theory and mechanisms. In Wear control handbook; American Society of Mechanical Engineers, 1980. [Google Scholar]
42. Transport Canada, Safe speed in ice: an analysis of transit speed and ice decision numerals. Ship Safety Northern (AMNS), 1996.
43. Transport Canada, Arctic Ice Regime Shipping System (AIRSS), Pictorial Guide, 03. 20 January.