Numerical Simulations of the Ice Load of a Ship Navigating in Level Ice Using Peridynamics

In this study, a numerical method was developed based on peridynamics to determine the ice loads for a ship navigating in level ice. Convergence analysis of threedimensional ice specimen with tensile and compression loading are carried out first. The effects of ice thickness, sailing speed, and ice properties on the mean ice loads were also investigated. It is observed that the ice fragments resulting from the icebreaking process will interact with one another as well as with the water and ship hull. The ice fragments may rotate, collide, or slide along the ship hull, and these ice fragments will eventually drift away from the ship. The key characteristics of the icebreaking process can be obtained using the peridynamic model such as the dynamic generation of cracks in the ice sheet, propagation and accumulation of ice fragments, as well as collision, rotation, and sliding of the ice fragments along the ship hull. The simulation results obtained for the ice loads and icebreaking process were validated against those determined from the Lindqvist empirical formula and there is good agreement between the results.


Introduction
As the demand for oil continues to increase over the years, oil companies are constantly searching for new oil reserves. It is estimated that 25% of the world's hydrocarbons are located in the Arctic. However, oil and gas exploration in the Arctic poses a significant challenge due to the extremely low temperatures and ice conditions in this region. With global warming on the rise, the thinning of ice sheets in the Arctic shortens the season in which activities can be carried out. However, the melting of sea ice driven by global warming will make it possible for commercial non-ice-strengthened ships to navigate along the edges of the Arctic, which will significantly reduce the distance between trading ports and thus, shipping times [Kwok and Cunningham (2012)]. Ship navigation in level ice leads to continuous breakup of the ice sheet into fragments. These ice fragments occupy a certain space within proximity of the ship, resulting in periodic and impulsive ice loads. The maximum ice load typically occurs within a relatively short period and most of the time, the ice load remains at a low level. Hence, it is crucial to fulfill the design requirements of ice-strengthened vessels such as bow structural strength, fatigue strength, ice-induced vibrations, and ice resistance. Numerous models have been developed to model ice-induced loads on the shell structure of ships. However, it is particularly challenging to develop reliable models that can accurately simulate the prevailing loads acting on the ship hull. Various methods have also been developed to predict the ice resistance of ship hulls. The reader may refer to previous studies [Daley (1991); Enkvist (1972); Kheisin and Popov (1973); Kujala (1994); Kujala and Vuorio (1986); Kwok and Cunningham (2012); Leira, Børsheim and Amdahl (2009); Varsta (1983) ;Vuorio, Riska and Varsta (1979)] for a detailed treatment on the statistical analysis of ice loads acting on ship hulls. Models have also been developed based on the finite element method (FEM) to predict ship-ice interactions and ice loads [Han, Feng and Yue (2007); Kolari, Kuutti and Kurkela (2009) ;Nylandsted, Jäättelä, Hoffmann et al. (2003); Premachandran and Horii (1994); Ranta, Polojärvi and Tuhkuri (2018); Xue (2016)]. The discrete element method (DEM) has been used to model various ice-related applications including real-time simulations of ship-ice interactions [Lubbad and Løset (2011)], modeling the ship behavior in broken ice [Karulin and Karulina (2011)], dynamic positioning in ice [Metrikin, Løset, Jenssen et al. (2013)], and ice floes jamming in rivers [Stockstill, Daly and Hopkins (2009)] as well as modeling ice ridges and ice rubble in conjunction with FEM [Arttu and Jukka (2009) ;Paavilainen and Tuhkuri (2013)]. Metrikin et al. ] provided an excellent review of different types of DEMs and the applications of DEMs in solving ice-related problems. However, due to differences in the discretization schemes between FEM and DEM, guaranteeing continuity at the solid boundaries using FEM-DEM coupled methods is not possible. Some of the disadvantages of using FEM in modeling the progressive failure of ice include mesh tangling and the use of erosion models to predict the failure patterns of ice [Das (2017)]. According to Di et al. [Di, Xue, Bai et al. (2017)] even though DEM is capable of capturing the discontinuous phase of the fracturing process, the method is less accurate when it comes to modeling the whole fracturing process from the continuous to discontinuous phase. Deb et al. [Deb and Pramanik (2013)] highlighted that when DEM is used, it is essential to conduct extensive calibration to identify the parameters for deformability and strength. Still model tests are the best method to measure the ice resistance of new ship hull designs, which have significantly improved over the last few years. Zhang et al. [Zhang, Shang and Zhang (2013)] conducted still model tests to determine the minimum engine power output for an oil tanker built for ice class 1A. The ice loads of an icebreaking vessel navigating in the Baltic Sea predicted by empirical formula and computational fluid dynamics were compared with those obtained from still model tests [Kim, Kim, Sim et al. (2015)]. A medium-scale still model test system was also built by Konno et al. [Konno, Nakane and Kanamori (2013)] to study ice-structure interactions, where the system was used to examine the different modes of dynamic ice loads and structural responses of narrow compliant structures such as the drilling platforms in Bohai Bay. Huang et al. [Huang, Ma and Tian (2013)] observed the ice failure modes and structural responses for different ice conditions and structural stiffnesses. Zhang et al. [Zhang, Wang, Gu et al. (2012)] established a speed measurement and positioning system for moving air cushion loads and they measured the mechanical properties of ice-like materials and the pressure distributions of air cushion loads. Luo et al. [Luo, Guo, Wu et al. (2018)] conducted ice resistance tests of a ship model under the combined effects of waves and ice floes and the results showed that the motion of the ship model was more unstable in the marginal ice zone than in ice floes. Even though still model tests in ice tanks are typically used to investigate the ice loads of ships navigating in Arctic areas, these tests are time-consuming and costly, especially in the design and optimization of hull forms [Kim, Lee, Lee et al. (2013)]. Meshfree methods are certainly advantageous to model the fracturing process and crack growth of ice sheets because these methods are independent of the element grids during the computational process [Han, Zhang, Wang et al. (2018); Li and Liu (2002); Wang, Wang, Zan et al. (2017)]. Libersky et al. [Libersky and Petschek (1991)] first applied the Smoothed Particle Hydrodynamics (SPH) approach to solid mechanics. Benz et al. [Benz and Asphaug (1994)] extended their work to simulate the fracturing process in brittle solids. Zhang et al. [Zhang, Zheng and Ma (2017)] embedded the Drucker-Prager model into the SPH code in order to simulate the four-point bending and uniaxial compression failure of ice. The cohesion softening elastic-plastic model has also been used within the SPH framework. Even though SPH is advantageous to numerically solve large deformation problems, its convergence has not been proven theoretically and therefore, there is a need to improve the accuracy of the solution [Li and Liu (2002)]. In order to overcome this limitation, Silling in 1998 coined peridynamics, which is reformulation of the classical continuum mechanics. In this approach, a force density function (also known as the pairwise force function) is used for each pair of particles and it is postulated that the particles interact with one another across a finite distance [Silling (2000); Silling and Askari (2005)]. Unlike molecular dynamics, peridynamic modeling has constitutive properties that are dependent on the reference configuration. In addition, unlike continuum mechanics, peridynamic modeling does not require a continuous response field [Gerstle (2015)]. Peridynamics is superior to solve discontinuous problems such as fracture in quasi-brittle solids [Ren, Zhuang, Cai et al. (2016) [Cheng, Zhang, Wang et al. (2015)]. A couple of studies have been carried out to simulate the ice loads and ice-structure interactions using peridynamics [Liu, Wang and Lu (2016) ; Ye, Wang, Chang et al. (2017); Wang, Wang, Zan et al. (2017);Liu, Xue, Lu et al. (2018)]. Even though there are a few studies concerning the use of peridynamics to simulate ice loads and ice-structure interactions, none of these studies are focused on the investigating the key characteristics of ship-ice interactions in continuous icebreaking mode, which forms the objective of this study. The effects of ice thickness, sailing speed, and ice properties on the mean ice loads during the continuous icebreaking process were also investigated. This paper is organized as follows. The kinematics of bond-based peridynamics is described briefly is in Section 2. Convergence analysis of three-dimensional ice specimen with tensile and compression loading and the determination of the micro-parameters of ice are presented in Section 3 while the simplifying assumptions (specifically the gravitational and buoyancy forces) and development of the ship-ice interaction model are presented in Section 4. The results obtained from the peridynamic model were compared with those determined from the Lindqvist empirical formula, which will be presented and discussed in Section 5. The effects of ice thickness, sailing speed, and ice properties on the mean ice loads are also presented in this section. Finally, the conclusions drawn based on the findings of this study are presented in Section 6.

Kinematics of bond-based peridynamics
In peridynamic modeling, the internal force acting at a particular point i in a continuum solid is defined as a function of the reference location Xj and the deformed location xj of the neighboring point j with respect to the reference location Xi and the deformed location xi of point i. The peridynamic equation of motion for point x in the reference configuration at time t is given by: where V ′ x is the volume of material point ′ x , Hx is called the family of x, u is the displacement vector field, b is the body force density field, ρ(x) is the density field, and ′′ u (x, t) is the acceleration field. The vector-valued function f is assumed to be a function of the undeformed bond (relative position), which is given by ξ=x'-x, and the deformed bond (relative displacement) of the reference points ′ x and x, which is given by η=u( ′ x , t)-u(x, t). The vector ξ+η describes the relative position of the reference points x′ and x in the current configuration. In peridynamic modeling, it is assumed that for a given material, the material points x and ′ x can only influence one another within a domain called horizon δ, as shown in Fig. 1.
The criterion for the undeformed bond ξ is given by: The damage and failure of the bond are contained by assuming that the bond force is solely dependent on the stretch (strain) of the bond, which is given by: It shall be noted that the value for the bond stretch s can be positive or negative, depending whether the bond is subjected to tensile or compressive loads. The material is assumed to be isotropic such that the bond stretch is not dependent on the direction of ξ. Failure occurs X X＇ H x , family of X if the bond stretch is greater than its predetermined value and the bond will no longer bear any force. The bond force for a prototype micro-elastic brittle material [Silling and Askari (2005)] is defined by: (4) where g(s)=cs and it is a linear scalar-valued function. Here, c is the bond constant, which is also known as the micro-modulus. Without going into the detailed derivation (Equating the expressions for strain energy density obtained by classical continuum mechanics and peridynamics in terms of both isotropic expansion and simple shear, see Madenci et al. [Madenci and Oterkus (2014)], one can obtain, where κ is bulk modulus, h equates size of particle in 2D, and A is cross-sectional area.
The parameter μ(t, ξ) is a history-dependent variable, which is defined for each bond that changes suddenly and irreversibly from 1 or 0 when the bond breaks. The parameter μ(t, ξ) is defined as: where s0 is a constant known as the critical bond strain, which is the strain when the bond breaks. In the prototype micro-elastic brittle material model, s0 can be calibrated based on the measured critical energy release rate G0 for a brittle material [Madenci and Oterkus (2014)]. The critical bond strain is given by: In peridynamic modeling, the crack growth is an autonomous process and it can be monitored by defining the local damage field φ(x, t), which is a fraction of the broken bonds connected to a point. The local damage field can be determined by:

Parameters determination 3.1 Peridynamic parameters for ice
Peridynamic body is discretized in a number of particles each describing some amount of volume and the side of a particle is commonly called a lattice. The results depend on particle positions, spacing dx, and material's horizon δ. Horizon mainly refers to the domain of influence which defines the range of interactions between material points. It can also be considered as a length-scale parameter which gives peridynamics a "non-local" character. The size of the horizon changes depending on the nature of the problem [Bobaru, Yang, Alves et al. (2009); Freimanis and Paeglitis (2017); Han, Lubineau, Yan et al. (2016)]. For three-dimensional analysis the most common are cubic lattices with uniform spacing in all directions, and the horizon may be chosen according to convergence, since for any value of δ, the parameters in the peridynamic material model can be chosen to match the measured Young's modulus of the material, and physical behaviour can be accurately represented.   [Schulson and Duval (2009)]. In this special case, the problem shows less "non-local" physical behavior, because the accuracy no longer increases with increasing horizon. Therefore, in practice, a sufficiently large horizon value can be chosen which does not show any significant non-local behavior, to make the computational time significantly decrease. M=3.015 is used in present model. One should note that in some cases, horizon can be used as a length scale parameter and can be adjusted in such a way that the required physical behavior can be accurately represented [Madenci and Oterkus (2014)]. Further study on the effect of the lattice size on critical strength as the horizon size stays constant M=3.015 is shown in Fig. 3. From Fig. 3, both compression and tension strengths are close to a stable value, except the lattice size is 17.5 mm (with values of 27501 and 5462 N showing more accurate). Therefore, the ratio of the characteristic scale to the lattice size should be no less than 4 (70/17.5) for getting reasonable results, however, additional studies may be required for this conclusion. One should note that the results obtained above are only applicable to three dimensional situations. For two-dimensional and one-dimensional problems, further studied need to be performed according to Eqs. (5) and (7).

Micro-parameters of ice
A large number of parameters such as the elastic modulus, Poisson's ratio, viscosity, ultimate bending strength, and compressive strength are required in order to ensure that the constitutive model adequately captures the complex physics of ice. These parameters are essential to reflect the physical and mechanical properties of ice in the numerical simulations. It is crucial to accurately determine the values of these parameters by carrying out a series of tests. However, in order to demonstrate the viability of the peridynamic model in the absence of such test data, these parameters can be estimated by using the test data published in Ji et al. [Ji, Wang, Jie et al. (2011); Wang, Ning and Ji (2014)] for icecovered waters in the Bohai Sea in winter. These test data were obtained from two typical model tests: (1) uniaxial compression test and (2) three-point bending test. Both of these tests were simulated in this study using the peridynamic model based on the assumption that ice is an isotropic prototype micro-elastic brittle material. The physical and mechanical properties of ice are defined by micro-parameters such as the micro-modulus c and critical bond strain s0. Still model test is an efficient approach to determine the micro-parameters of ice, where the micro-parameters are varied to cater for the macro properties of the specimens, which will improve the accuracy of the numerical model. The readers may refer to Liu et al. [Liu, Xue, Lu et al. (2018)] for details on the simulation methods.

Simplifying assumptions 4.1 Gravitational and buoyancy forces
A few assumptions need to be made in order to simplify the ship-ice interaction model. Firstly, it was assumed that the gravitational force is balanced by the buoyancy force of the ice sheet in calm waters. When the ship (icebreaker) is present, the ice sheet deviates from its equilibrium position due to imbalance of the gravitational and buoyancy forces, which will affect the motion of the ice sheet. It was assumed that the buoyancy force acting on the ice block is approximately linear to the change in position of the material point. Consequently, Δf, which is a function of the gravitational force density g, as well as the volume force density b of a material point will vary linearly with changes in the position of the material point. For example, assuming that the center coordinates of point x in the vertical direction is x with the waterline taken as the datum, Δf can be determined by: where ρi is the density of ice, ρw is the density of water, lw is the length of the material point immersed in the water, and dx is the scale of the material point.

Ship-ice contact model
The contact force model used in this study is defined by the repelling short-range force [Silling (2000)]. The force between two material points can be expressed as: where csh is the short-range force constant and rsh is the critical distance. The short-range force is non-existent if the critical distance is exceeded. Here, csh=5⋅c while rsh=dx/2 [Madenci and Oterkus (2014)]. In this study, the ice material model ( An indentation is produced in the ice sheet during the initial stage of the icebreaking process, as shown in Fig. 7. Force is generated when the ship and ice come into contact. Since the contact area is small in this stage, the force is not sufficiently large to induce crack growth in areas of the ice sheet far from the contact area. Compressive damage and failure occur in the contact area. The indentation gradually enlarges and cracks begin to form within vicinity of the contact area. The color scale is the contour plot of damage variable φ(x, t) defined by Eq. (8), and is the same in other figures.

Figure 7: Indentation and crack initiation within the contact area
Radial cracks begin to appear as the icebreaking process continues, as shown in Fig. 8. In this stage, the tensile strength of ice is lower than its compressive strength and thus, cracks begin to form at the lower surface of the ice sheet and the radial cracks propagate along with the ship motion. As the radial cracks continue to develop in the ice sheet, the ice surface near the contact area is split into separate areas, forming wedge-shaped cantilevers. The force acting on the ice gradually increases until the ice sheet can no longer support the continuous bending of the cantilevers. The upper surface of the ice sheet is subjected to tensile loads whereas the lower surface of the ice sheet is subjected to compressive loads and therefore, circumferential cracks form on the upper surface. These circumferential cracks eventually penetrate through the ice sheet across its thickness, as shown in Fig. 9. One should note that this simulation of ship-ice interaction in the level ice with fixed boundaries, and the boundary regions of the ice plate contains three layers of material points on each edge which are fixed and free of deformation. These fixed boundaries have a substantial influence on the interaction process shown in Figs. 6, 8 and 9. However, for calculating ice load, it is reasonable to study the level ice with the size of the channel. It can be observed from Figs. 10 and 11 that the level ice is broken into fragments and these ice fragments then drift away from the ship bow. In reality, the channel will be opened over a continuous period, but in this study, the ship-ice interactions are only considered within a small fraction of the period in which the icebreaking process takes place. Most of the time, the ship sails through the ice rubble or ice-free conditions, which results in impulsive ice loads. The simulation results shown in Figs. 10 and 11 conform well with the characteristics of the icebreaking process described above. It is evident that the cracks are fully developed at 1.74 s in the initial stage and the icebreaking process enters the following stage at 8.68 s. The ice loads imposed on the ship bow are rather low during this process.  Kujala (1994)], indicating the reliability of the peridynamic model. The magnitudes of the ice loads are relatively low in the x-axis due to the rigid motion and symmetry of the ship bow model. Hence, the ice loads are symmetrical at both sides of the ship in the simulations, but in reality, the non-simultaneous failure of ice results in asymmetry and fluctuations of the ice loads. The magnitudes of the ice loads are higher at 10 s, which is indeed expected as the fracture forms an asymmetrical contact surface. Based on the load-time history curve in the z-direction, the magnitudes of the vertical ice loads are indeed high, which is sufficient to lift the bow.

Effect of ice thickness on the ice load
Simulations were conducted using the peridynamic model in order to investigate the effect of ice thickness (0.50 m, 0.75 m, and 1.00 m) on the ice loads acting on the ship bow. The Ice load [N] ship was navigating in level ice at a sailing speed of 3 kn. The load-time history curves are shown in Fig. 13. The mean ice load was plotted against the ice thickness, as shown in Fig. 14, and it can be observed that the ice thickness has a significant effect on the mean ice load. The mean ice load acting on the ship bow in continuous icebreaking mode are 386801.37 N, 710446.18 N, and 1107258.38 N for an ice thickness of 0.50 m, 0.75 m, and 1.00 m, respectively. The ice loads were then calculated using the Lindqvist empirical formula [Kujala (1991); Lindqvist (1989); Matlock, Dawkins and Panak (1969)]. The breaking resistance Rc is defined as: where φ=tan -1 (tan β/sin α), α is the rake angle, β is the broadside rake angle, and μ is the ship-ice kinetic friction coefficient. The bending resistance Rb can be expressed in terms of the ice thickness hi as follows: where υ is the Poisson's ratio, σf is the bending strength of ice, B is the breadth of ship, and ρw is the density of sea water. The submersion resistance Rs is defined as: where ρi is the density of sea ice, htot is the thickness of ice and snow, T is the vertical distance between the water surface and lowest point of the ship hull (otherwise known as draft), and L is the length of the ship. The total resistance R is defined as: where V is velocity of the ship. In current study, much less ice sliding force is present in such a short simulation. Therefore, Rs should not be included in next calculation. Because Lindqvist assumed that Rs is calculated when 70% of the ship is covered by broken ice, that is why there is the parameter 0.7 in Eq. (15). The values of the parameters used in the Lindqvist empirical formula are presented in Tab. 5.  The mean ice loads determined from the peridynamic model and Lindqvist empirical formula are summarized in Tab. 6. Fig. 14 shows the mean ice loads determined from the peridynamic model and Lindqvist empirical formula plotted against the ice thickness. It can be seen that the mean ice load is almost linearly depended on the square of the ice thickness when the ship is traveling at a speed of 3 kn. The mean ice load at an ice thickness of 1.00 m is 3.13 times higher than that at an ice thickness of 0.50 m, indicating that ice thickness is indeed a factor that significantly affects the ice load. The mean ice loads determined from the peridynamic model are consistent with those determined from the Lindqvist empirical formula, which proves the viability of the peridynamic model. The difference between the results obtained from the peridynamic model and Lindqvist empirical formula may be due to the lack of consideration on the role of water and frictional effects between the ship and ice in the peridynamic model. Figure 14: Effect of ice thickness on the mean ice load determined from the peridynamic model and Lindqvist empirical formula

Effect of sailing speed on the ice load
In order to investigate the effect of sailing speed on the ice loads imposed on the ship bow, the sailing speed of the ship was varied in the simulations (3, 4, 5, 6, 7, and 8 kn) while the thickness of the ice sheet was fixed at 0.50 m. The load-time history curves are shown in Fig. 15.  16 shows the mean ice load plotted against the sailing speed and it can be observed that the mean ice load slightly increases when the sailing speed of the ship is increased from 3 to 5 kn. The mean ice load significantly increases when the sailing speed is increased beyond 5 kn.

Effect of ice properties on the ice load
The compressive strength and tensile strength of the materials were determined based on the critical stretch of the bond in tensile and compression, respectively. It is known that the elastic modulus of ice is related to its compressive strength and tensile strength. Hence, if there are changes in the elastic modulus of ice, there will be changes in the bond critical stretch. For this reason, the effects of elastic modulus and compressive strength of ice on the ice loads were investigated in this study and three cases were considered: Case I, Case II, and Case III. The ship was navigating in level ice with a thickness of 0.5 m and a sailing speed of 3 kn. In Case I, the elastic modulus of ice was varied at 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, and 3.0 GPa and its effect on the mean ice load acting on the ship bow are shown in Fig. 17. In Case II, the compressive strength of ice was fixed at 5.3 MPa whereas its elastic modulus was varied at 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, and 3.0 GPa by adjusting the critical stretch of the bond. The results are shown in Fig. 18.  for Case II (where only the elastic modulus of ice is varied while its compressive strength is kept constant) show that there are no significant variations in the mean ice load as the elastic modulus of ice increases. This is also confirmed by the results obtained from the Lindqvist empirical formula. The results for Case III (where the compressive strength is varied while the elastic modulus of ice is kept constant) show that the mean ice load increases with an increase in the compressive strength. The following correlations were obtained by applying linear fitting to the data for each case, which can be used to predict the mean ice load acting on the ship bow as a function of the elastic modulus or compressive strength of ice.
For Case II, the value of the slope is -0.029611 with the unit of meganewtons (MN). This indicates that the elastic modulus of ice affects the mean ice load by ~3.0% if the compressive strength of ice is constant. For Case III, the value of the slope is 0.076331, with ~2.6 times of the effect of elastic modulus. The results presented in this section demonstrate the effect of elastic modulus and compressive strength of ice on the ice loads imposed on the ship bow. In general, the compressive strength of ice has a higher contribution to the ice load.

Conclusions
In this study, a model was developed based on peridynamic theory in order to simulate the ship-ice interactions of a ship navigating in level ice in continuous icebreaking mode. The ice loads and key characteristics of the icebreaking process were captured at different ice thicknesses and sailing speeds. The effect of the physical and mechanical properties of ice on the ice loads was also investigated. The following conclusions were drawn based on the simulation results: Convergence study showed the different ratios of the horizon to lattice size show similar trend of change in both compression and tension. The horizon size of 3.015 lattice size give the most accurate result in compression; however, the horizon size of 4.015 lattices shows the most accurate result in tension. Additionally, given the horizon size of 3.015 lattice size, the ratio of the size of model to lattice size no less than four gives more accurate results. It can be generalized that the number of particles on any dimension should be more than the ratio of the horizon to lattice size. The mean ice load increases almost linearly depended on the square of the ice thickness. The mean ice load at an ice thickness of 1.00 m is 3.13 times higher than that at an ice thickness of 0.50 m.
The mean ice load slightly increases as the sailing speed of the ship is increased from 3 to 5 kn. However, the ice load increases significantly when the sailing speed is increased beyond 5 kn. Within the peridynamic framework, it is found that the compressive strength of ice has a higher contribution to the mean ice load about 2.6 times of the elastic modulus of ice.