Xing Zheng,Zhizong Tian,Zhigang Xie and Ningbo Zhang

Received:28 June 2021/Accepted:23 January 2022

?Harbin Engineering University and Springer-Verlag GmbH Germany,part of Springer Nature 2022

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.

Keywords Icebreaker;Smoothed-particle hydrodynamics;Ice breaking resistance;the Yellow River;Ice thickness

1 Introduction

Given its high latitude and severe cold in winter,the Yel‐low 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 phenome‐non 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 in‐clude ice breaking by blasting (Xie et al. 2021), icebreak‐ers (Cheng 2013), and so on. Compared with other types of ice breaking, icebreakers have the advantages of stron‐ger mobility and lower cost. Thus, icebreakers are impor‐tant 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 dis‐crete 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), Peridyn‐amics(PD)(Xue et al.2020),and smoothed particle hydro‐dynamics(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 non‐continuous problem for mesh-based methods. Ice resis‐tance 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 up‐dating 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 uni‐axial 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 ap‐plied for ice dynamic simulations. Xue et al. (2018) estab‐lished an elastic brittle failure model of ice material by us‐ing 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 dynam‐ics 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 me‐chanics,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 float‐ing 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 consti‐tutive 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 dy‐namic 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 mod‐el 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 re‐duced 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 nu‐merical model by using SPH to predict and analyze ice re‐sistance.The ice–water coupling model was based on a sim‐ple repulsive force model without considering the broken ice fields. Zhang et al. (2019b) and Zhang (2020) estab‐lished 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 resis‐tance. 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 icebreak‐er in the Yellow River is predicted based on SPH.For SPH,it can be easily applied for fluid and ice dynamics.Further‐more,on the basis of the previous work of ice–water inter‐actions, SPH has been applied for ice cracking and break‐ing,which is unsuitable for mesh-based methods.The con‐struction of the method includes the elastic–plastic consti‐tutive 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. Previ‐ous 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 interac‐tion 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.Accord‐ing to the practical application of this SPH for ice–ship in‐teraction, 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 in‐teraction force and a new ideal for SPH.The inlet and out‐let 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:

where D/Dtis the particle derivative,αandβare the Carte‐sian component in thex,y,zdirections,ρis the particle density,νis the particle velocity,gis the acceleration of gravity,andσis the stress tensor.In this paper,an artificial diffusion term proposed by Antuono et al. (2010) is intro‐duced into the continuity equation to eliminate the unrea‐sonable high-frequency oscillation of pressure field in flu‐id simulation. This method is called delta-SPH. Sun et al.(2017) reported the use of delta-SPH. Only the main for‐mulas 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:

Figure 1 Schematic of fluid and solid particles at interface boundary

whereρiandmiare the density and mass of particleiwith velocity componentvi, respectively, andρjandmjare the density and mass of particlejwith the velocity componentvj.The discrete form of momentum equation of ice particles is

In this paper, the simplified finite difference interpola‐tion method is used to calculate the strain rate of ice parti‐cles.In addition,the cubic spline kernel function is applied to simulate the smooth kernel function,and the virtual par‐ticle 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 inter‐action zone(Figure 1)is highly important.A simple and ef‐fective fluid–ice interface scheme is proposed to deal with the interaction between the fluid and ice particles in the in‐teraction 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 particleialso act as the dummy particles for imposing boundary conditions.For example,when cal‐culating the momentum equations of fluid particleiwith fluid and solid particles in its neighbor supporting domain,solid particlejin the computational domain of fluid parti‐clei(Figure 1(a)) is considered in the calculation of the momentum equations:

When acting as dummy particles, the corresponding pressure of ice particlejin the above equation can be inter‐polated using the neighboring fluid particles by Eq. (4b),which can be referred to Adami et al.(2012).

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:

whereδαβis the Kronecker sign,Gis the shear modulus,λis the plastic product factor,Qis the plastic potential func‐tion determining the development direction of plastic strain,e?αβ=ε?αβ?1/3ε?γγδαβis the partial shear strain rate tensor,K=E/(3(1 ?2υ)) is the elastic bulk modulus,Eis Young’s modulus,G=E/(2(1 +υ))is the shear modulus,J2is the second invariant of the stress tensor, andηˉis the function of dilatancy angle.

During the simulation of ice failure, the stress tensor must be modified, and the cohesive force softening meth‐od must be used to reduce the cohesive force of ice in the plastic flow stage. The damage model of ice is not an im‐portant research content of this paper; thus, it will not be carried out in detail, and the specific process can be re‐ferred 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 flu‐id and ice phases, which can be shown by the black curve in Figure 2. First, the geometry of the icebreaker is gener‐ated in SolidWorks. Then, tetrahedral meshes are generat‐ed 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 parti‐cles in Figure 2(a).The boundary particles on the hull sur‐face 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 fol‐lows. When solving the governing equation of the ice phase, these ship boundary particles in the computing do‐main of the ice particles can be regarded as virtual parti‐cles to impose boundary conditions.These boundary parti‐cles participate in the continuity and momentum equa‐tions 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 comput‐ing 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 ve‐locity boundary conditions of the fluid phase.The calcula‐tion process is shown in Figure 3.

Figure 2 Sketch of the ship–ice–water interaction model

In 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 real‐ized. 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

Figure 3 SPH framework of the ship–ice–water program including the elastic–plastic ice model

The ice breaking resistance under different conditions can be predicted using the above numerical model.Accord‐ing to the characteristics of the ice condition on the Yellow River,the calculation results of ice breaking resistance un‐der the condition of the level ice are provided as a refer‐ence 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 par‐ticles 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 estimat‐ed using the surface integral of the stress tensor of adja‐cent ice particles. The force of ice particles on the bound‐ary particles can be expressed as

whereFisthe interaction force betweentheice particlesand the ship boundary particles,σαis the stress t ensor,Wijis the kernel function,andmandρare the mass and density ofeach 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 lengthL=12 m,widthB=9 m,and thicknessH.In this part,one case is consideredH=0.06 m. The level ice is fixed with no motion in six directions. Some characteristic pa‐rameters 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/m3and cohesionc=17.0 kPa for ice thicknessH=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 an‐gleφ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 ARAON

Figure 5 Ice breaking resistance time histories with different ship velocities when H=0.06 m

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 break‐ing 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 averag‐ing the corresponding ice resistances in Figure 5 aftert= 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.Before presenting the results of the Yellow River Ice‐breaker, 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 re‐sistance, and bending failure resistance (Lindqvist 1989).The main characteristics of the Lindqvist formula are list‐ed in Table 2. The formula of ice breaking resistance is as follows:

Figure 6 Comparison of ice breaking resistances with different ship velocities of SPH results with experimental data when H=0.06 m

whereRCis the crushing resistance,RBis the bending resis‐tance,Rsis the immersion resistance,Lis the ship length,hiis the thickness of ice, andVis the ship speed. The cal‐culation formula of crushing resistanceRCis as follows:

whereσfis the bending strength,φ1is the bow column an‐gle,ψ= arctan(tanφ1/sinα)is the drift angle,andαis the waterline angle. The expression of bending resistanceRBis as follows:

whereBis the ship width.The expression of immersion re‐sistanceRsis as follows:

whereρwis the density of sea water,ρiis the density of sea ice,Tis the draft, andμis the friction coefficient between the ship and the ice.

Table 2 Characteristics of the Lindqvist formula

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 ob‐tained by SPH has some noise, but it can be maintained at a certain range after it is stable. When the stable ice resis‐tance 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 m

The results of SPH and the Lindqvist empirical formula are compared when the ice thickness isH=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 cal‐culation results in Figure 7 are averaged aftert=10 s,and the calculation formula isF?=∑F/nto extract the deter‐mined 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.

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 his‐tories 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 ten‐dency 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 oscil‐lation also increases.The SPH results of Table 4 show that the case of ice thicknessH=0.6 m is approximately 3.01 times compared with the case of ice thicknessH= 0.3 m. The SPH results in Table 5 indicate that the case of ship veloci‐tyV=4 kn is 1.76 times compared with the case of ship ve‐locityV=1 kn.All the results of SPH with different veloci‐ties and ice thicknesses agree with the results of the Lindqvist formula.

Figure 10 shows the comparison of the ice breaking re‐sistance 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 break‐ing 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 in‐creases by approximately 2.5 times. This finding suggests that the ice thickness is the most important factor influenc‐ing ship–ice breaking resistance.

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

Table 3 Comparison of ice breaking resistance with different bending strengths

Figure 11 shows the results of ice breaking patterns according to different ice thicknesses whent= 20.0 s.Figure 5 considers the level ice case, where the ice thick‐ness isH= 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 increas‐ing 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 af‐ter the ice breaking.

Table 4 Comparisons of ice breaking resistance with different ice thicknesses

Table 5 Comparison of ice breaking resistance with different velocities

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 differ‐ent ice thicknesses when t=20.0 s

4 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 ice‐breakers 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 model

4.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 concentrationC, and the average size of the floating ice floe. The ice concentrationCis 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 aver‐age ratio of the total ice covered area against the number of ice floes.Then,the same number of random points are gen‐erated in the river model area, and the Voronoi polygon is generated according to the random points to ensure the dis‐tance between any point in the polygon and the correspond‐ing 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 parti‐cles.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 fields

The flow field distribution of the river must be ob‐tained 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 solv‐er 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 flowU= 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)

whereAis the contact area between ice particles and the flow field,andVis the flow velocity.

The framework of the SPH code clearly is introduced in Figure 14. First, the SPH code is based on the origi‐nal SPH formulas. Second, the ice dynamic model is considered. The velocity distributions of the river chan‐nel are given on the basis of the results of iceFoam.The current flow force is determined according to the veloci‐ty distribution, which can be added to the moment equa‐tion 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 be‐tween the ice and the ship is calculated, and results are outputted.

4.3 Numerical simulation

Figure 15 shows the simulation results of icebreaker operation progress in the Yellow River with the ice con‐centrationC= 70% at the simulation timest= 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 atU= 4 kn. Figure 15 shows that SPH can ef‐fectively 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 re‐sistance of the icebreaker under different ice concentra‐tions 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 concentra‐tion is small withC= 50% and obviously increases when the concentration reachesC= 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 rela‐tively small. With the increase in concentration, the ac‐tion time of ice resistance becomes continuous and the amplitude of ice force increases. The ice concentration significantly affects the time distribution but not the am‐plitude of ice resistance.×106

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 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 in‐creases, the amplitudes of high-frequency parts are larger than the ones of the higher concentrations.

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

The effect of ice thickness on the ice breaking resistance of the icebreaker in the case of broken ice fields is inves‐tigated.Figure 18 provides a comparison of the time histo‐ry of ice resistance with different ice thicknesses calculat‐ed by SPH when the ship velocity is 4 kn and the ice densi‐ty is 70%.The simulated time history is 120 s in Figure 18.Comparison of the results in Figure 18 shows that the peri‐odic trend of ice resistance exhibits no obvious change,but the amplitude of ice resistance displays some signifi‐cant changes under different ice thicknesses. The maxi‐mum amplitudes of ice breaking resistance are approxi‐mately 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 thick‐nessH= 0.6 m is 2.7 times compared with the case of the thicknessH= 0.3 m. This result indicates that ice thick‐ness is the main factor influencing ice breaking resistance.

Figure 19 presents the FFT results against the frequen‐cies and amplitudes for the case on Figure 18 with dif‐ferent ice thicknesses.As displayed in Figure 19, the av‐erage ice resistance increases with ice thickness.The fre‐quencies of different ice thicknesses do not change obvi‐ously. 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 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%

5 Conclusion

A ship–ice–water coupling numerical model is estab‐lished using the SPH method. According to the compari‐son 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 fol‐lowing 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 for‐mula;

2) For the case of breaking ice simulation, generating the breaking ice on the river by the Voronoi diagram meth‐od is a good method.The ice floe model can realize the ini‐tial 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 resis‐tance in the real river.The thickness of ice floe is the most important factor affecting ice breaking resistance.

Founding InfromationSupported by the National Key Research and Development Program of China (No. 2018YFC1508405),National Natural Science Foundation of China (Nos.51879051 and 51739001),the Open Fund of the Key Laboratory of Far-shore Wind Power Technology of Zhejiang Province (ZOE20200007), Natural Science Foundation of Heilongjiang Province in China (LH2020 E071).