Numerical Study of the Ice Breaking Resistance of the Icebreaker in the Yellow River Through Smoothed-Particle Hydrodynamics
https://doi.org/10.1007/s11804-022-00259-w
-
Abstract
A ship – ice – water interaction model is established using smoothed-particle hydrodynamics (SPH) to predict the ice breaking resistance of the icebreaker in the Yellow River effectively. This method includes the numerical process of the constitutive equation, yield criterion, and the coupling model in SPH. The ice breaking resistance is determined under different conditions. The numerical results of the ice breaking resistance agree with the empirical formula results. Results show that the prediction accuracy of ice resistance is less than 17.6% compared with the empirical formula in the level ice. The method can also be extended to predict the floe motion and ice breaking resistance in actual river channels. The validation against the empirical formula indicates that the proposed ship – ice – water SPH method can predict the ice breaking resistance of icebreakers in actual rivers effectively. The predicted ice breaking resistance is analyzed under different conditions. The ice breaking resistance increases with increasing bending strength and ice thickness, and the latter is the most important factor influencing ice resistance.Article Highlights• Ice floes in a river are generated using the Voronoi polygon method based on SPH. Different ice concentrations can be considered.• Ice resistance characteristics of the level ice, including the bending strengths moments, ice thicknesses, and ship velocities, are obtained.• Smoothed-particle hydrodynamics (SPH) is applied to predict the ice breaking resistance of the icebreaker in the Yellow River. -
1 Introduction
Given its high latitude and severe cold in winter, the Yellow River always freezes in winter. In the Inner Mongolia basin, the northernmost part of the Yellow River, ice floes easily accumulate in the narrow and curved section of the river under the action of current in the early stage of river closure and the thawing period of the river. Accumulated ice floes can form ice rafts and ice dams, which can block the river channel and cause ice disasters. This phenomenon has caused huge losses to people's life and property along the coast (Gao et al. 2019). Ice breaking measures that dredge accumulated ice floes are necessary to deal with ice disasters. At present, ice breaking measures include ice breaking by blasting (Xie et al. 2021), icebreakers (Cheng 2013), and so on. Compared with other types of ice breaking, icebreakers have the advantages of stronger mobility and lower cost. Thus, icebreakers are important in ice breaking and dredging ice accumulation in some river sections. Effective numerical analysis can be used to estimate the ice breaking resistance of the Yellow River icebreaker. Results of such analysis are important in the design and operation of the icebreaker.
Computational fluid dynamics (CFD) is commonly used in fluid domain calculation but not in ice dynamics, which includes ice cracking and breaking. CFD is unsuitable for mesh-based methods, such as StarCCM+ coupled with discrete element method (DEM) (Hao et al. 2020; Huang et al. 2021) and with LS-DYNA simulation (Li et al. 2020). Meshless methods include DEM (Di et al. 2017), Peridynamics (PD) (Xue et al. 2020), and smoothed particle hydrodynamics (SPH) (Zhang et al. 2019a; Khayyer et al. 2021a; 2021b). SPH is a particle-based method used to simulate ice breaking and cracking with a suitable government equation that can generate the singular point and the noncontinuous problem for mesh-based methods. Ice resistance prediction based on SPH considers ice floe dynamic interaction. Ice floes can be broken by the ship–ice and ice–ice interactions. Moreover, the ice – water interaction can be simulated by SPH.
Kolari et al. (2009) proposed a finite element simulation method for continuous ice failure based on the model updating technology, in which the anisotropic continuum damage mechanics model is used to predict the crack growth direction, and the model updating technique is used to predict the crack growth in the finite element model. Carne et al. (2006) established a phenomenological failure model of ice and solved the plastic sensitive failure model of ice by using the finite element program of LS-DYNA. Pernas et al. (2012) established an elastic–plastic material model based on the Drucker–Prager (D–P) yield criterion and simulated the ice breaking behavior under a high stress rate. They solved the model with the finite element code LS-DYNA and integrated the Lagrangian, ALE, and SPH methods into the model to obtain the impact force on the slender cylindrical ice sample. Di et al. (2015) established a discrete element model to simulate the contact between sea ice and two elements and studied ice failure under uniaxial compression and three-point bending. The simulation results of this study agree with the experimental results. Kong et al. (2021) analyzed the ice load of a polar floating platform based on DEM. Furthermore, PD theory is applied for ice dynamic simulations. Xue et al. (2018) established an elastic brittle failure model of ice material by using the PD theory and simulated the three-point bending failure of ice beams. They compared the numerical results with the experimental data and found good consistency.
The aforementioned methods are focused on ice dynamics and do not consider the effect of ice–water interaction. SPH is a type of meshless Lagrangian particle algorithm. With the rapid development of SPH in fluid and solid mechanics, SPH has also been applied to the simulation of ice dynamics. For example, Gutfraind and Savage (1997) and Oger and Savage (1999) applied rheology based on the Mohr Coulomb yield criterion in SPH to simulate the floating and movement of broken ice floes on the water surface under the action of wind. Shen et al. (2000) proposed a 2D numerical model in which SPH was used to simulate the drift motion and accumulation blocking of river ice. Ji et al. (2005) proposed a new dynamic viscoelastic plastic constitutive model of sea ice in which SPH was used to simulate the ice motion in a rectangular basin. In addition, Ji et al. (2007) developed a hybrid Lagrangian Eulerian sea ice dynamic method in which the sea ice cover is represented by an SPH model with its own thickness and concentration. Pan et al. (2012) proposed a new SPH non-Newtonian model to study the coupled dynamics of ice sheets and ice shelves. Das (2017) used the SPH model to simulate the four-point bending failure of ice beam in LS-DYNA, in which the von Mises yield criterion was used to judge the failure of ice particles. In this method, once the ice particles reach the failure state, the deviatoric stress component is reduced to zero. Zhang et al. (2017) used an improved SPH method combined with the D–P yield criterion and elastic–plastic constitutive model of cohesive softening to simulate the bending and compression failure characteristics of the ice. Qiao (2018) established a ship–ice–water coupling numerical model by using SPH to predict and analyze ice resistance. The ice–water coupling model was based on a simple repulsive force model without considering the broken ice fields. Zhang et al. (2019b) and Zhang (2020) established the ship–ice and ship–ice–wave coupling numerical models by using SPH to simulate the ice failure, wave–ice interaction, and ice load on the hull. Studies showed that SPH is feasible and effective in studying ice breaking resistance. Xue et al. (2020) summarized different numerical methods for ice–ship interactions and found that SPH is an important numerical tool for this research.
In this study, the ice breaking resistance of the icebreaker in the Yellow River is predicted based on SPH. For SPH, it can be easily applied for fluid and ice dynamics. Furthermore, on the basis of the previous work of ice–water interactions, SPH has been applied for ice cracking and breaking, which is unsuitable for mesh-based methods. The construction of the method includes the elastic–plastic constitutive equation of the ice, the ship – ice – water numerical model, the broken ice fields, and the interaction between ice floes with the river flow field. The numerical method is validated against the Lindqvist empirical formula. Previous SPH works can include three aspects for this topic. First, SPH can provide the results of ice foe drift in a very large area, and it cannot consider the ice–ship interaction. Second, SPH is focused on ice dynamics. It considers the failure of ice impact and the crash development. Third, SPH can obtain reliable results for ship–ice–water interaction for some typical icebreakers, and these icebreakers are in an open area. The boundary of the ice cover and the inlet and outlet of the ice boundary are not considered. According to the practical application of this SPH for ice–ship interaction, the novelty of the proposed SPH code is mainly focused on the application for ice resistance prediction. The governing equation of the ice is a new application to SPH. The ice and river flow interaction is based on the interaction force and a new ideal for SPH. The inlet and outlet of ice floes are also the keys to solving these problems.
2 Establishment of the mathematical model
2.1 Basic equations of SPH
In SPH, the basic governing equations include the mass and momentum conservation equations in the Lagrangian form, which are expressed as follows:
$$ \frac{{{\rm{D}}\rho }}{{{\rm{D}}t}} = - \rho \frac{{\partial {v^\alpha }}}{{\partial {x^\alpha }}} $$ (1) $$ \frac{{{\rm{D}}{v^\alpha }}}{{{\rm{D}}t}} = \frac{1}{\rho }\frac{{\partial {\sigma ^{\alpha \beta }}}}{{\partial {x^\beta }}} + {g^\alpha } $$ (2) where D/Dt is the particle derivative, α and β are the Cartesian component in the x, y, z directions, ρ is the particle density, ν is the particle velocity, g is the acceleration of gravity, and σ is the stress tensor. In this paper, an artificial diffusion term proposed by Antuono et al. (2010) is introduced into the continuity equation to eliminate the unreasonable high-frequency oscillation of pressure field in fluid simulation. This method is called delta-SPH. Sun et al. (2017) reported the use of delta-SPH. Only the main formulas are included because delta-SPH is not the focus of this paper. The particle approximate expression of the mass conservation equation of solid phase is as follows:
$$ \frac{{{\rm{D}}{\rho _i}}}{{{\rm{D}}t}} = {\rho _i}\sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}} \left({v_i^\alpha - v_j^\alpha } \right)\frac{{\partial {W_{ij}}}}{{\partial x_i^\alpha }} $$ (3) where ρi and mi are the density and mass of particle i with velocity component vi, respectively, and ρj and mj are the density and mass of particle j with the velocity component vj. The discrete form of momentum equation of ice particles is
$$ \begin{array}{l} \frac{{{\rm{d}}v_i^\alpha }}{{{\rm{d}}t}} = \sum\limits_{j = 1}^N {{m_j}} \left({\frac{{\sigma _i^{\alpha \beta }}}{{\rho _i^2}} + \frac{{\sigma _j^{\alpha \beta }}}{{\rho _j^2}} - {\Pi _{ij}} \cdot {\delta ^{\alpha \beta }}} \right.\\ \;\;\;\;\;\;\;\;\;\;\left. { + f_{ij}^n\left({R_i^{\alpha \beta } + R_j^{\alpha \beta }} \right)} \right)\frac{{\partial {W_{ij}}}}{{\partial x_i^\beta }} + {g^\alpha } \end{array} $$ (4) where the artificial viscous term Πij was proposed by Monaghan (1994) to improve the stability of numerical calculation and reduce the unstable oscillation, Πij = (− acijμij + βμij2)/pij, μij = (huij·rij) /(rij2 + 0.01h2), h = Δx, α = 0.1, β = 0 fijn(Riαβ + Rjαβ) is the artificial repulsive term. fij is defined as fij = Wij /W (Δd, h), where Δd is the initial distance between neighbor particles. The Riαβ and Rjαβ in Eq. (4) are the artificial stress tensor of particles i and j, respectively with the correction parameter ε =0.2, Riαβ =− εσiαβ /ρ2 when σiαβ > 0. Additional details are provided in Monaghan (2000).
In this paper, the simplified finite difference interpolation method is used to calculate the strain rate of ice particles. In addition, the cubic spline kernel function is applied to simulate the smooth kernel function, and the virtual particle method is used to deal with the solid wall boundary.
During the wave–ice interaction simulation, the contact algorithm between the fluid and solid particles in the interaction zone (Figure 1) is highly important. A simple and effective fluid–ice interface scheme is proposed to deal with the interaction between the fluid and ice particles in the interaction domain. In this new method, the ice particles act as dummy particles to approximate the interface between the fluid phase and the solid body. When solving the fluid domain equations, the solid particles in the computational domain of fluid particle i also act as the dummy particles for imposing boundary conditions. For example, when calculating the momentum equations of fluid particle i with fluid and solid particles in its neighbor supporting domain, solid particle j in the computational domain of fluid particle i (Figure 1(a)) is considered in the calculation of the momentum equations:
$$ \begin{array}{l} \frac{{{\rm{d}}v_i^\alpha }}{{{\rm{d}}t}} = \sum\limits_{j \in {\rm{ ice }}}^N {{m_j}} \left({ - \frac{{{p_i} + {p_j}}}{{{\rho _i}{\rho _j}}}} \right)\frac{{\partial {W_{ij}}}}{{\partial x_i^\alpha }}\\ \; + \sum\limits_{j \in {\rm{ fluid }}}^N {{m_j}} \left({ - \frac{{{p_i} + {p_j}}}{{{\rho _i}{\rho _j}}} - {\Pi _{ij}}} \right)\frac{{\partial {W_{ij}}}}{{\partial x_i^\alpha }} + {g^\alpha } \end{array} $$ (4a) 
Figure 1 Schematic of fluid and solid particles at interface boundaryWhen acting as dummy particles, the corresponding pressure of ice particle j in the above equation can be interpolated using the neighboring fluid particles by Eq. (4b), which can be referred to Adami et al. (2012).
$$ {p_j} = \frac{{\sum\limits_{i \in {\rm{ fluid }}}^M {{p_i}} {W_{ji}} + \left({g - {a_j}} \right) \cdot \sum\limits_{i \in {\rm{ fluid }}}^M {{\rho _i}} {\mathit{\boldsymbol{r}}_{ji}}{W_{ji}}}}{{\sum\limits_{i \in {\rm{ fluid }}}^M {{W_{ji}}} }} $$ (4b) Additional details can be found in the study by Zhang (2020).
2.2 Damage model of ice
The elastic – plastic constitutive model is applied to SPH to simulate the failure of the ice. Combined with the D–P yield criterion, the stress–strain equation of the ice model with non-correlated flow law can be obtained as follows:
$$ {\dot \sigma ^{\alpha \beta }} = \left\{ \begin{array}{l} {\sigma ^{\alpha \gamma }}{{\dot \omega }^{\beta \gamma }} + {\sigma ^{\gamma \beta }}{{\dot \omega }^{\alpha \gamma }} + 2G{{\dot e}^{\alpha \beta }} + K{{\dot \varepsilon }^{\gamma \gamma }}{\delta ^{\alpha \beta }}\quad \;\;F\left({{\sigma ^{\alpha \beta }}, c} \right) < 0\\ {\sigma ^{\alpha \gamma }}{{\dot \omega }^{\beta \gamma }} + {\sigma ^{\gamma \beta }}{{\dot \omega }^{\alpha \gamma }} + 2G{{\dot e}^{\alpha \beta }} + K{{\dot \varepsilon }^{\gamma \gamma }}{\delta ^{\alpha \beta }}\\ - \dot \lambda \left({\bar \eta K{\delta ^{\alpha \beta }} + \frac{G}{{\sqrt {{J_2}} }}{s^{\alpha \beta }}} \right)\quad \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;F\left({{\sigma ^{\alpha \beta }}, c} \right) \ge 0 \end{array} \right. $$ (5) where δαβ is the Kronecker sign, G is the shear modulus, λ is the plastic product factor, Q is the plastic potential function determining the development direction of plastic strain, $ {\dot e^{\alpha \beta }} = {\dot \varepsilon ^{\alpha \beta }} - 1/3{\dot \varepsilon ^{\gamma \gamma }}{\delta ^{\alpha \beta }}$is the partial shear strain rate tensor, K=E/(3(1 − 2υ)) is the elastic bulk modulus, E is Young's modulus, G = E/ (2 (1 + υ))is the shear modulus, J2is the second invariant of the stress tensor, and $ \bar \eta $ is the function of dilatancy angle.
During the simulation of ice failure, the stress tensor must be modified, and the cohesive force softening method must be used to reduce the cohesive force of ice in the plastic flow stage. The damage model of ice is not an important research content of this paper; thus, it will not be carried out in detail, and the specific process can be referred to Zhang (2020).
2.3 Ship–ice–water interaction model
The SPH model of the ship – ice – water interaction scheme is provided in this section. The ship hull of the icebreaker is treated as the solid wall boundary for the fluid and ice phases, which can be shown by the black curve in Figure 2. First, the geometry of the icebreaker is generated in SolidWorks. Then, tetrahedral meshes are generated using Gambit, and the SPH particles are placed on these grid nodes to realize the arrangement of the ship's solid boundary, which can be shown by the black particles in Figure 2(a). The boundary particles on the hull surface should be sufficiently dense to prevent the ice and fluid particles from penetrating the ship boundary. These boundary particles in the solution of the ice and fluid phases act as dummy particles for imposing boundary conditions. The specific implementation process is as follows. When solving the governing equation of the ice phase, these ship boundary particles in the computing domain of the ice particles can be regarded as virtual particles to impose boundary conditions. These boundary particles participate in the continuity and momentum equations of the ice phase, so that the ice particles near the hull boundary satisfy the continuity of stress and velocity. Similarly, when solving the fluid governing equations, these ship boundary particles in the fluid particle computing domain also act as the virtual particles for the fluid phase. These boundary particles are also included in the solution of the continuity and momentum conservation equations of the fluid phase to satisfy the pressure and velocity boundary conditions of the fluid phase. The calculation process is shown in Figure 3.
Figure 2 Sketch of the ship–ice–water interaction model
Figure 3 SPH framework of the ship–ice–water program including the elastic–plastic ice modelIn addition, the fluid–ice interface scheme introduced by Zhang et al. (2019b) can be used to simulate the fluid– ice interaction. The slip boundary condition is imposed by neglecting the viscous interaction between virtual particles on the hull boundary and ice or fluid particles. Through the above method, the interaction of the fluid phase, the ice phase, and the ship boundary can be realized. The calculation process of this boundary treatment is simple; that is, the complex geometric information of the hull boundary and coupling interface is unnecessary to calculate.
3 Numerical simulation and verification of the level ice
The ice breaking resistance under different conditions can be predicted using the above numerical model. According to the characteristics of the ice condition on the Yellow River, the calculation results of ice breaking resistance under the condition of the level ice are provided as a reference for the ice breaking load in the actual river.
3.1 Ship model generation
The general layout of the light icebreaker in the Yellow River is shown in Figure 4(a). The total length of the ship is 28.0 m, the shape width is 7.0 m, the design draft is 1.0 m, and the propulsion power is 2×500 kW. The 3D model of the ship is shown in Figure 4(b).
Figure 4 Icebreaker model diagram (28 m)The uniform distribution of particles on the hull surface is established to represent the effect of the hull boundary to simulate the ship–ice interaction by SPH. Then, the ice load on the ship is obtained by solving the force of ice particles on the boundary particles of the ship. When these hull boundary particles act as virtual particles, the force of ice particles on the hull boundary particles can be estimated using the surface integral of the stress tensor of adjacent ice particles. The force of ice particles on the boundary particles can be expressed as
$$ {F_{{\rm{ice - shipboundary }}}} = \sum\limits_{i \in {\rm{ boundary }}}^N {{m_i}} \sum\limits_{j \in {\rm{ ice }}}^N {{m_j}} \left({\frac{{{{\rm{ \mathsf{ σ} }}_i}}}{{\rho _i^2}} + \frac{{{{\rm{ \mathsf{ σ} }}_j}}}{{\rho _j^2}}} \right){\nabla _i}{W_{ij}} $$ (6) where Fice−shipboundaryα is the interaction force between the ice particles and the ship boundary particles, σα is the stress tensor, Wij is the kernel function, and m and ρ are the mass and density of each particle, respectively.
3.2 Simulation and verification of ice resistance in the level ice
The results of the ship – ice interaction in the level ice are presented in this section. The ice breaking resistances on the icebreaker under different situations are predicted using SPH to validate the numerical results. The Korean icebreaker ARAON is used, and the scale factor is 1:20 in the numerical simulation (Lau and Akinturk 2011). The main parameters of icebreaker ARAON in the full and model scales are shown in Table 1. The level ice has the length L=12 m, width B=9 m, and thickness H. In this part, one case is considered H=0.06 m. The level ice is fixed with no motion in six directions. Some characteristic parameters are considered, such as the thickness of the ice sheet and the icebreaker velocities, which are as same as those in the study by Lau and Akinturk (2011). The ice density ρ =864.0 kg/m3 and cohesion c =17.0 kPa for ice thickness H=0.06 m are the same as those in the model test (Lau and Akinturk, 2011). The elastic modulus is 25.0 MPa, and the friction angle is 22.5°. The dilatancy angle φ in the non-associative plastic rule is set to be onethird of the friction angle φ = ϕ/3.
Table 1 Parameters for the full and model scales of icebreaker ARAONParameters (m) Full scale Model scale Overall length 107.2 5.36 Breadth 19.0 0.95 Draft 6.8 0.34 The ice breaking resistance time histories in the ice–ship interaction are shown in Figure 5. The resistance time histories computed by SPH are generally stable and reasonable, although some oscillations exist. The noise of the ice breaking resistance curve may be caused by the accumulation of crushed ice around the ice breaker. The general trend of ice breaking resistance increases with the ship speed and the level of ice thickness. Figure 6 illustrates the comparison of the time series of ice breaking resistance of ship against the ship velocities in the experimental data with the results calculated by SPH. The results of SPH in Figure 6 are obtained by averaging the corresponding ice resistances in Figure 5 after t = 4 s, where the ice breaker has basically sailed into the level ice, and the resistance value becomes stable in general. The trends of the ice resistance increasing with icebreaker speed agree with the results of the model test even though some errors exist. In Figure 6, the predicted results from the presented SPH model overestimate the ice breaking resistance for the model test data.
Figure 5 Ice breaking resistance time histories with different ship velocities when H = 0.06 m
Figure 6 Comparison of ice breaking resistances with different ship velocities of SPH results with experimental data when H = 0.06 mBefore presenting the results of the Yellow River Icebreaker, the Lindqvist empirical formula is introduced to verify the results of the ice breaking resistance. The Lindqvist empirical formula is based on the assumption that the ice resistance has a linear relationship with the ship velocity and divides the ice breaking resistance into three parts, including extrusion resistance, immersion resistance, and bending failure resistance (Lindqvist 1989). The main characteristics of the Lindqvist formula are listed in Table 2. The formula of ice breaking resistance is as follows:
$$ {R_{{\rm{ice}}}} = \left({{R_C} + {R_B}} \right)\left({1 + 1.4\frac{V}{{\sqrt {g{h_i}} }}} \right) + {R_S}\left({1 + 9.4\frac{V}{{\sqrt {gL} }}} \right) $$ (7) Table 2 Characteristics of the Lindqvist formulaIce density ρi (kg/m3) 920 Ship–ice friction Coefficient μ 0.15 Ice bending strength σi (MPa) 0.6/0.9/1.2 Ship length L 28 Bow column inclination φ1 (°) 40 Ice thickness hi (m) 0.3 / 0.4 / 0.5 / 0.6 Water density ρw (kg/m3) 998 Draft T (m) 1.0 Ship width B (m) 7.0 Waterline entry angleα (°) 33.5 Ship velocity V (kn) 1.0 / 2.0 / 3.0 / 4.0 where RC is the crushing resistance, RB is the bending resistance, Rs is the immersion resistance, L is the ship length, hi is the thickness of ice, and V is the ship speed. The calculation formula of crushing resistance RC is as follows:
$$ {R_C} = 0.5{\sigma _f}h_i^2\left({\tan {\varphi _1} + \mu \frac{{\cos {\varphi _1}}}{{\cos \psi }}} \right)/\left({1 - \mu \frac{{\sin {\varphi _1}}}{{\cos \psi }}} \right) $$ (8) where σf is the bending strength, φ1is the bow column angle, ψ = arctan (tan φ1 / sin α)is the drift angle, and α is the waterline angle. The expression of bending resistance RB is as follows:
$$ {R_B} = 0.003{\sigma _f}Bh_i^{1.5}\left({\frac{{\tan \psi + \mu \cos {\varphi _1}}}{{\cos \psi \sin \alpha }}} \right)\left({1 + \frac{1}{{\cos \psi }}} \right) $$ (9) where B is the ship width. The expression of immersion resistance Rs is as follows:
$$ {R_S} = \left({{\rho _w} - {\rho _i}} \right)g{h_i}B\left[ {\begin{array}{*{20}{l}} {\frac{{T(B + T)}}{{B + 2T}} + \mu \left({0.7L - \frac{T}{{\tan {\varphi _1}}} - \frac{B}{{4\tan \alpha }}} \right.}\\ {\left. {T\cos {\varphi _1}\cos \psi \sqrt {\frac{1}{{{{\sin }^2}{\varphi _1}}} + \frac{1}{{{{\tan }^2}\alpha }}} } \right)} \end{array}} \right] $$ (10) where ρw is the density of sea water, ρi is the density of sea ice, T is the draft, and μ is the friction coefficient between the ship and the ice.
Figure 7 shows the time history comparison between SPH and the Lindqvist empirical formula. As shown in Figure 7, the time history of the ice breaking resistance obtained by SPH has some noise, but it can be maintained at a certain range after it is stable. When the stable ice resistance results are achieved, the SPH results with different bending strengths oscillate near the results of the Lindqvist empirical formula.
Figure 7 Comparison of time histories of ice breaking resistance with different bending strengths when the ship velocity is 4 kn and ice thickness is 0.3 mThe results of SPH and the Lindqvist empirical formula are compared when the ice thickness is H = 0.3 m, the ship speed is 4 kn, and the ice bending strengths are 0.6, 0.9, and 1.2 MPa to verify the SPH coupling method. The calculation results in Figure 7 are averaged after t = 10 s, and the calculation formula is $ \tilde F = \sum F /n$ to extract the determined values for comparison. The comparison results are shown in Table 3. As shown in Table 3, the errors between SPH and the Lindqvist empirical formula are 15.7%, 14.0%, and 17.6%, when the bending strengths are 0.6, 0.9, and 1.2 MPa, respectively. This finding shows that the ice breaking resistance calculated by SPH agrees with that predicted by the Lindqvist empirical formula.
Table 3 Comparison of ice breaking resistance with different bending strengthsIce thickness (m) Bending strength (MPa) Ship velocity (kn) Lindqvist (kN) SPH (kN) Error (%) 0.3 0.6 4 1.79 1.51 15.7 0.3 0.9 4 2.67 2.30 14.0 0.3 1.2 4 3.55 2.93 17.6 Comparisons of the SPH results with different cases are performed. Time histories of the ice breaking resistance with different ice thicknesses when the ship velocity is 4 kn and the ice bending strength is 0.6 MPa are displayed in Figure 8 and Figure 9 displays the comparisons of time histories of the ice breaking resistance with different velocities when the ice thickness is 0.6 m and the ice bending strength is 1.2 MPa. Tables 3 and 4 illustrate the comparisons of the ice breaking resistance with different ice thicknesses and different velocities. As shown in Figures 8 and 9, the tendency of the ice resistance still has very large oscillations, but the average value is close to the results of the Lindqvist formula. With increasing ice thickness, the ice resistance increases quickly, and the amplitude of the oscillation also increases. The SPH results of Table 4 show that the case of ice thickness H = 0.6 m is approximately 3.01 times compared with the case of ice thickness H = 0.3 m. The SPH results in Table 5 indicate that the case of ship velocity V= 4 kn is 1.76 times compared with the case of ship velocity V=1 kn. All the results of SPH with different velocities and ice thicknesses agree with the results of the Lindqvist formula.
Figure 8 Comparisons of time histories of ice breaking resistance with different ice thicknesses when the ship velocity is 4 kn and ice bending strength is 0.6 MPa
Figure 9 Comparisons of time histories of ice breaking resistance with different velocities when the ice thickness is 0.6 m and ice bending strength is 1.2 MPaTable 4 Comparisons of ice breaking resistance with different ice thicknessesIce thickness (m) Bending strength (MPa) Ship velocity (kn) Lindqvist (105N) SPH (105N) Error (%) 0.3 0.6 4 1.79 1.51 15.7 0.4 0.6 4 2.72 2.53 7.0 0.5 0.6 4 3.80 3.25 14.5 0.6 0.6 4 5.00 4.56 8.8 Table 5 Comparison of ice breaking resistance with different velocitiesIce thickness (m) Bending strength (MPa) Ship velocity (kn) Lindqvist (105N) SPH (105N) Error (%) 0.6 1.2 1 5.90 5.07 14.1 0.6 1.2 2 7.25 6.95 4.2 0.6 1.2 3 8.60 8.13 5.5 0.6 1.2 4 9.95 8.91 10.5 Figure 10 shows the comparison of the ice breaking resistance calculated by SPH under different ice thicknesses and different bending strengths to study the ice breaking resistance with different ice thicknesses. Comparison of the results in Figure 10 shows that the change in ice breaking resistance increases linearly with the bending strength, whereas the ice breaking resistance changes dramatically with different ice thicknesses. As shown in Figure 10, when the ice thickness increases by one time and the other conditions remain the same, the ice breaking resistance increases by approximately 2.5 times. This finding suggests that the ice thickness is the most important factor influencing ship–ice breaking resistance.
Figure 10 Comparison of ice breaking resistance according to different ice thicknesses and different bending strengthsFigure 11 shows the results of ice breaking patterns according to different ice thicknesses when t = 20.0 s. Figure 5 considers the level ice case, where the ice thickness is H = 0.3, 0.4, 0.5, 0.6 m. The color bar represents the value of the accumulated plastic strain εp, which shows the ice breaking extent. As displayed in Figure 11, the crack propagation near the bow is not obvious when the ice thickness is small, and the boundary of the broken massive ice floe is unclear. However, the degree of ice fragmentation near the bow is more obvious with increasing ice thickness, and clear massive ice floes appear after ice breaking. This finding shows that SPH can predict not only the ice breaking resistance but also the ice cracks after the ice breaking.
Figure 11 Comparison of ice breaking patterns according to different ice thicknesses when t=20.0 s4 Simulation and analysis of ice breaking resistance in the river
4.1 River model building
The ice resistance in the river is predicted and analyzed in this section to study the ice breaking resistance of icebreakers in the Yellow River. First, numerical modeling of the river channel is introduced. On the basis of the satellite image of a section of the Yellow River, as shown in Figure 12(a), the edge model of the river is drawn, as shown in Figure 12(b). Two river sections at the outlet and the inlet of the river are established using SolidWorks, and then the profile shape at each different position is obtained by induced stretching through the terrain on both sides of the river, thereby creating a 3D river profile.
Figure 12 Setting up of river channel model4.2 Broken ice fields
After the river model is established, the SPH particle model must be established for broken ice floes in the river. In this paper, the Voronoi diagram is used to build the ice floe model. The number of floating ice floes in the broken ice fields is calculated according to the area of the river channel, the ice concentration C, and the average size of the floating ice floe. The ice concentration C is calculated according to the ratio of the actual broken ice floes and the coverage area. The total area is the corresponding river channel area. The average size of floating ice is the average ratio of the total ice covered area against the number of ice floes. Then, the same number of random points are generated in the river model area, and the Voronoi polygon is generated according to the random points to ensure the distance between any point in the polygon and the corresponding random points of the polygon is the closest. Finally, the Voronoi polygon is contracted according to the selected ice concentration to determine the filling area of SPH particles. With this method, the calculation model of broken ice fields in the river channel is shown in Figure 13(c).
Figure 13 Setting up of broken ice fieldsThe flow field distribution of the river must be obtained to determine the accumulation and movement of ice floes in the river. SPH is difficult to use to solve large-scale river flow velocity directly. Thus, icoFoam in OpenFOAM is applied to solve the problem. The solver is used to solve the unsteady incompressible Navier Stokes equations with the finite volume method. The hexahedral grid is generated by blockmesh, and then the STL file of the river surface is read by snappyhexmesh to draw the outline of the river and refine the grid. The velocity of flow U = 0.5, 0.7, 1.0 m/s can be set at the river inlet. After obtaining the velocity of the flow field on the grid of the river surface by icoFoam, the force of the flow field on the ice particles can be obtained using the following formula (Pan, 1986)
$$ F = 0.59A{V^2} $$ (11) where A is the contact area between ice particles and the flow field, and V is the flow velocity.
The framework of the SPH code clearly is introduced in Figure 14. First, the SPH code is based on the original SPH formulas. Second, the ice dynamic model is considered. The velocity distributions of the river channel are given on the basis of the results of iceFoam. The current flow force is determined according to the velocity distribution, which can be added to the moment equation of ice particles. Third, some damage models of the ice are added to calculate the interaction between the ice particles and ship boundary particles. Last, the force between the ice and the ship is calculated, and results are outputted.
Figure 14 Framework of the SPH code and the coupling of the icoFoam4.3 Numerical simulation
Figure 15 shows the simulation results of icebreaker operation progress in the Yellow River with the ice concentration C = 70% at the simulation times t = 40, 80, and 120 s. The ship model is shown in Figure 4. Starting from the entrance of the river, the trajectory coincides with the main channel of the river, and its speed remains constant at U = 4 kn. Figure 15 shows that SPH can effectively simulate the drift of ice floes under the action of the icebreaker in the river. In addition, obvious ice breaking traces appear at the tail of the ship. The ice concentration on the ice breaking resistance is analyzed. Figure 16 displays the time history comparison of ice resistance of the icebreaker under different ice concentrations when the ship speed is 4 kn and the ice thickness is 0.3 m. As demonstrated in Figure 16, the ice resistance of the icebreaker is relatively small when the concentration is small with C = 50% and obviously increases when the concentration reaches C = 90%. The ice load curve in Figure 16 illustrates that when the density is small, the action time of ice resistance is relatively short, and the maximum force of each period is also relatively small. With the increase in concentration, the action time of ice resistance becomes continuous and the amplitude of ice force increases. The ice concentration significantly affects the time distribution but not the amplitude of ice resistance.×106
Figure 15 Numerical results of icebreaker sailing in the river with the ice concentration C=70%
Figure 16 Comparison of time histories of ice breaking resistance with different concentrations when the ship velocity is 4 kn and the ice thickness is 0.3 mFigure 17 shows the FFT results against the frequencies and amplitudes for the case in Figure 16. As displayed in Figure 17, the average ice resistance increases with ice concentration. However, the frequencies of different ice concentrations are unclear. When ice concentration increases, the amplitudes of high-frequency parts are larger than the ones of the higher concentrations.
Figure 17 Comparisons of the results of FFT of time histories of ice breaking resistance when the ship velocity is 4 kn and the ice thickness is 0.3 mThe effect of ice thickness on the ice breaking resistance of the icebreaker in the case of broken ice fields is investigated. Figure 18 provides a comparison of the time history of ice resistance with different ice thicknesses calculated by SPH when the ship velocity is 4 kn and the ice density is 70%. The simulated time history is 120 s in Figure 18. Comparison of the results in Figure 18 shows that the periodic trend of ice resistance exhibits no obvious change, but the amplitude of ice resistance displays some significant changes under different ice thicknesses. The maximum amplitudes of ice breaking resistance are approximately 2.13×105, 3.11×105, 4.93×105, and 5.81×105N for ice thicknesses 0.3, 0.4, 0.5, and 0.6 m, respectively. Figure 18 shows that the maximum ice resistance of the ice thickness H = 0.6 m is 2.7 times compared with the case of the thickness H = 0.3 m. This result indicates that ice thickness is the main factor influencing ice breaking resistance.
Figure 18 Comparison of time histories of ice breaking resistance with different ice thicknesses when the ship velocity is 4 knots and the ice concentration is 70 %Figure 19 presents the FFT results against the frequencies and amplitudes for the case on Figure 18 with different ice thicknesses. As displayed in Figure 19, the average ice resistance increases with ice thickness. The frequencies of different ice thicknesses do not change obviously. In the case of different ice thicknesses of ice floes, the amplitudes of its resistances can increase obviously, and the parts of the frequency are not the main factor.
Figure 19 Comparisons of the FFT of time histories of ice breaking resistance when the ship velocity is 4 kn and the ice concentration is 70%5 Conclusion
A ship – ice – water coupling numerical model is established using the SPH method. According to the comparison with the Lindqvist empirical formula, the ice breaking resistance range of the icebreaker in the Yellow River is studied by the results under different conditions. The following conclusions can be drawn:
1) In accordance with the D–P yield criterion, the error of ice breaking resistance obtained by SPH is less than 17.6% at level ice cases with the Lindqvist empirical formula;
2) For the case of breaking ice simulation, generating the breaking ice on the river by the Voronoi diagram method is a good method. The ice floe model can realize the initial division of different broken ice fields by adjusting the ice concentration;
3) With the help of icoFoam solution, SPH can realize the calculation of ice floe motion and ice breaking resistance in the real river. The thickness of ice floe is the most important factor affecting ice breaking resistance.
-
Figure 1 Schematic of fluid and solid particles at interface boundary
Figure 2 Sketch of the ship–ice–water interaction model
Figure 3 SPH framework of the ship–ice–water program including the elastic–plastic ice model
Figure 4 Icebreaker model diagram (28 m)
Figure 5 Ice breaking resistance time histories with different ship velocities when H = 0.06 m
Figure 6 Comparison of ice breaking resistances with different ship velocities of SPH results with experimental data when H = 0.06 m
Figure 7 Comparison of time histories of ice breaking resistance with different bending strengths when the ship velocity is 4 kn and ice thickness is 0.3 m
Figure 8 Comparisons of time histories of ice breaking resistance with different ice thicknesses when the ship velocity is 4 kn and ice bending strength is 0.6 MPa
Figure 9 Comparisons of time histories of ice breaking resistance with different velocities when the ice thickness is 0.6 m and ice bending strength is 1.2 MPa
Figure 10 Comparison of ice breaking resistance according to different ice thicknesses and different bending strengths
Figure 11 Comparison of ice breaking patterns according to different ice thicknesses when t=20.0 s
Figure 12 Setting up of river channel model
Figure 13 Setting up of broken ice fields
Figure 14 Framework of the SPH code and the coupling of the icoFoam
Figure 15 Numerical results of icebreaker sailing in the river with the ice concentration C=70%
Figure 16 Comparison of time histories of ice breaking resistance with different concentrations when the ship velocity is 4 kn and the ice thickness is 0.3 m
Figure 17 Comparisons of the results of FFT of time histories of ice breaking resistance when the ship velocity is 4 kn and the ice thickness is 0.3 m
Figure 18 Comparison of time histories of ice breaking resistance with different ice thicknesses when the ship velocity is 4 knots and the ice concentration is 70 %
Figure 19 Comparisons of the FFT of time histories of ice breaking resistance when the ship velocity is 4 kn and the ice concentration is 70%
Table 1 Parameters for the full and model scales of icebreaker ARAON
Parameters (m) Full scale Model scale Overall length 107.2 5.36 Breadth 19.0 0.95 Draft 6.8 0.34 Table 2 Characteristics of the Lindqvist formula
Ice density ρi (kg/m3) 920 Ship–ice friction Coefficient μ 0.15 Ice bending strength σi (MPa) 0.6/0.9/1.2 Ship length L 28 Bow column inclination φ1 (°) 40 Ice thickness hi (m) 0.3 / 0.4 / 0.5 / 0.6 Water density ρw (kg/m3) 998 Draft T (m) 1.0 Ship width B (m) 7.0 Waterline entry angleα (°) 33.5 Ship velocity V (kn) 1.0 / 2.0 / 3.0 / 4.0 Table 3 Comparison of ice breaking resistance with different bending strengths
Ice thickness (m) Bending strength (MPa) Ship velocity (kn) Lindqvist (kN) SPH (kN) Error (%) 0.3 0.6 4 1.79 1.51 15.7 0.3 0.9 4 2.67 2.30 14.0 0.3 1.2 4 3.55 2.93 17.6 Table 4 Comparisons of ice breaking resistance with different ice thicknesses
Ice thickness (m) Bending strength (MPa) Ship velocity (kn) Lindqvist (105N) SPH (105N) Error (%) 0.3 0.6 4 1.79 1.51 15.7 0.4 0.6 4 2.72 2.53 7.0 0.5 0.6 4 3.80 3.25 14.5 0.6 0.6 4 5.00 4.56 8.8 Table 5 Comparison of ice breaking resistance with different velocities
Ice thickness (m) Bending strength (MPa) Ship velocity (kn) Lindqvist (105N) SPH (105N) Error (%) 0.6 1.2 1 5.90 5.07 14.1 0.6 1.2 2 7.25 6.95 4.2 0.6 1.2 3 8.60 8.13 5.5 0.6 1.2 4 9.95 8.91 10.5 -
Adami S, Hu XY, Adams NA (2012) A generalized wall boundary condition for smoothed particle hydrodynamics. Journal of Computational Physics 231(21): 7057-7075 https://doi.org/10.1016/j.jcp.2012.05.005 Antuono M, Colagrossi A, Marrone S, Molteni D (2010) Free-surface flows solved by means of SPH schemes with numerical diffusive terms. Computer Physics Communications 181(3): 532-549 https://doi.org/10.1016/j.cpc.2009.11.002 Cheng HT (2013) Research on a new method of ice prevention using icebreaker in Yellow River. Master thesis, North China University of Water Resources and Electric Power, Zhengzhou, 5-10. (in Chinese) Das J (2017) Modeling and validation of simulation results of an ice beam in four-point bending using smoothed particle hydro dynamics. International Journal of Offshore and Polar Engineering 27(1): 82-89 https://doi.org/10.17736/ijope.2017.ik01 Di SC (2015) Discrete element simulation of ice load on offshore platform and ship hull based on GPU parallel algorism. PhD thesis, Dalian University of Technology, Dalian, 1-5 (in Chinese) Di SC, Ji SY, Xue YZ (2017) Analysis of ship navigation in level ice covered regions with discrete element method. The Ocean Engineering 35(3): 59-69. (in Chinese) Gao GM, Deng Y, Tian ZZ, Li SX, Zhang BS (2019) Brief introduction and prospect of recent ice research in the Yellow River. Yellow River 41(10): 77-78. (in Chinese) Gutfraind R, Savage SB (1997). Smoothed particle hydrodynamics for the simulation of broken-ice fields: mohr-coulomb-type rheology and frictional boundary conditions. Journal of Computational Physics 134(2): 203-215 https://doi.org/10.1006/jcph.1997.5681 Hao YH, Zheng X, Ma QW (2020) Calculation of ice breaking resistance of level ice based on CFD-DEM coupling method. Proceedings of the International Offshore and Polar Engineering Conference, Shanghai, 590-595 Huang LF, Li MH, Romu TM, Dolatshah A, Thomas G (2021) Simulation of a ship operating in an open-water ice channel. Ships and Offshore Structures 16(4): 353-362 https://doi.org/10.1080/17445302.2020.1729595 Ji SY, Li H, Shen HT, Wang R, Yue Q (2007) A hybrid lagrangianeulerian numerical model for sea-ice dynamics. Acta Oceanologica Sinica 26(5): 12-24 Ji SY, Shen HT, Wang ZL, Shen HH, Yue Q (2005) A viscoelasticplastic constitutive model with mohr-coulomb yielding criterion for sea ice dynamics. Acta Oceanologica Sinica 24(4): 54-65 Khayyer A, Gotoh H, Shimizu Y, Nishijima Y (2021a) A 3D Lagrangian meshfree projection-based solver for hydroelastic Fluid-Structure Interactions. Journal of Fluids and Structures 105, 103342 https://doi.org/10.1016/j.jfluidstructs.2021.103342 Khayyer A, Shimizu Y, Gotoh H, NagashimaK (2021b) A coupled incompressible SPH-Hamiltonian SPH solver for hydroelastic FSI corresponding to composite structures. Applied Mathematical Modelling 94: 242-271 https://doi.org/10.1016/j.apm.2021.01.011 Kolari K, Kuutti J, Kurkela J (2009) FE-simulation of continuous ice failure based on model update technique. Proceedings of the 20th International Conference on Port and Ocean Engineering under Arctic Conditions, Lulea, Sweden, 100-115. DOI: 10.1109/APPEEC.2009.4918848 Kong S, Ji SY, Ji SP, Wang YY, Gang XH (2021) Numerical analysis of ice load on floating platform in polar region based on highperformance discrete element method. Chinese Journal of Ship Research 16(5), 1-7 (in Chinese) Lau M, Akinturk A (2011) KORDI Araon model tests in ice using the planar motion mechanism, report LM-2011-04. Canada: St. John's, Newfoundland and Labrador, National Research Council-Institute for Ocean Technology Li F, Kujala MKP, Goerlandt F (2020) Finite element based metamodeling of ship-ice interaction at shoulder and midship areas for ship performance simulation. Marine Structure 71, 102736 https://doi.org/10.1016/j.marstruc.2020.102736 Lindqvist GA (1989) Straightforward method for calculation of ice resistance of ships. Proceedings of the Port and Ocean Engineering under Arctic Conditions, Lulea, Sweden, 722-735 Monaghan JJ (1994) Simulating free surface flows with SPH. Journal of Computational Physics 110, 399-406 https://doi.org/10.1006/jcph.1994.1034 Monaghan JJ (2000) SPH without a tensile instability. Journal of Computational Physics 159(2), 290-311 https://doi.org/10.1006/jcph.2000.6439 Oger L, Savage SB (1999) Smoothed particle hydrodynamics for cohesive grains. Computer Methods in Applied Mechanics and Engineering 180(1-2): 169-183 https://doi.org/10.1016/S0045-7825(99)00054-7 Pan SH (1986) Building standards and codes of the State Construction Committee of the Council of Ministers of the Soviet Union. Chapter 57. Loads and actions of waves, ice and ships on hydraulic structures, Ocean Press. (in Chinese) Pan W, Tartakovsky AM, Monaghan JJ (2012) A smoothed-particle hydrodynamics model for ice-sheet and ice-shelf dynamics. Journal of Glaciology 58(208): 828-842 Pernas-Sánchez DA, Pedroche DV, López-Puente RZ (2012) Numerical modeling of ice behavior under high velocity impacts. International Journal of Solids and Structures 49: 1919-1927 https://doi.org/10.1016/j.ijsolstr.2012.03.038 Qiao Y (2018) Study on the coupled of ship-ice-water interaction and the prediction method of ice resistance. PhD thesis, Harbin Engineering University, Harbin, 1-6. (in Chinese) Shen HT, Su J, Liu L (2000) SPH simulation of river ice dynamics. Journal of Computational Physics, 165(2): 752-770 https://doi.org/10.1006/jcph.2000.6639 Sun PN, Colagrossi A, Marrone S, Zhang AM (2017) The δplus-SPH model: Simple procedures for a further improvement of the SPH scheme. Computer Methods in Applied Mechanics and Engineering 315: 25-49 https://doi.org/10.1016/j.cma.2016.10.028 Xie ZG, Zen H, Li G, Deng Y, Tian ZZ (2021) The timing and method of eliminating the ice dam with airborne bomb by unmanned aerial vehicle. Yellow River 43(2): 70-72. (in Chinese) Xue YZ, Liu RW, Li Z, Han D (2020) A review for numerical simulation methods of ship-ice interaction. Ocean Engineering 215, 107853 https://doi.org/10.1016/j.oceaneng.2020.107853 Xue YZ, Lu XK, Wang Q (2018) Simulation of three-point bending test of ice based on peridynamic. Journal of Harbin Engineering University 39(4), 607-613 Zhang NB (2020) Research on the ship-ice-wave coupling numerical model based on the SPH method. PhD thesis, Harbin Engineering University, Harbin, 56-65. (in Chinese) Zhang NB, Zheng X, Ma QW (2017) Updated Smoothed Particle Hydrodynamics for simulating bending and compression failure progress of ice. Water 9(11): 882 https://doi.org/10.3390/w9110882 Zhang NB, Zheng X, Ma QW, Hu ZH (2019a) A numerical study on ice failure process and ice-ship interactions by Smoothed Particle Hydrodynamics. International Journal of Naval Architecture and Ocean Engineering 11(2): 796-808 https://doi.org/10.1016/j.ijnaoe.2019.02.008 Zhang NB, Zheng X, Ma QW (2019b) Study on wave-induced kinematic responses and flexures of ice floe by Smoothed Particle Hydrodynamics. Computers & Fluids 189, 46-59
