Bao-yu Ni,Duan-feng Han,Shao-cheng Di,Yan-zhuo Xue

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001,China

Abstract:Ice-water-structure interaction (IWSI)is a novel extension of the fluid-structure interaction (FSI),which is significant for design and operating of polar ship and offshore structures. It involves multi-media and multi-interfacesand thus is quite complicated to solve,no matter from mathematical or mechanical perspectives.Although IWSI is complex and still very new,researchers try to develop various methods to deal with it and relevant literature starts to bloom.This paper aims to provide concise descriptions of typical analytical numerical and experimental methods to solve IWSI, together with a review of their major applications to date. Lastly,we succinctly highlight some development tendencies and some pieces of work to be investigated for each method.

Key words:Ice-water-structure interaction (IWSI),ice,polar ships,ocean engineering,fluid-structure interaction (FSI)

Biography of First author:

Bao-yu Ni(1986-),Professor and Doctoral Supervisor of College of Shipbuilding Engineering,Harbin Engineering University (HEU).He received his bachelor and doctor degree in Naval Architecture and Ocean Engineering from HEU in 2008 and 2012,respectively,during which he visited University College London (UCL)in UK for one year between 2010 and 2011 as a visiting PhD candidate.After graduation,he worked as a teacher in HEU and visited UCL again for two years between 2015 and 2017 as a post-doctor.His research fields of interest include ice-water-ship interaction,motion characteristics of polar ships in ice,bubble dynamics,water impact,etc.He has undertaken over 10 projects such as supported from EU Horizon 2020 Project,the National Natural Science Foundation of China,and National Key Research and Development Program of China,etc..He has published 2 monographs and over 60 journal/conference papers,authorized 7 invention patents and 7 software copyrights.He has been as one of the editorial board members for the Journal of Hydrodynamics(JHD),Shock and Vibration,Journal of Marine Science and Application and Chinese JHD.He has been selected into “Youth Talent Promotion Project”of China and “Postdoctoral International Exchange Program”of China,etc..

Introduction

With the global warming,the sea ice around the polar region is reducing gradually[1].The Arctic shipping route and polar resource development become more realizable,so more ships and offshore structures operate in polar region[1-2].Arctic shipping route provides a shortcut between Pacific and Atlantic ports,which saves time and fuel cost compared with the traditional routes.Take the European route from Rotterdam to Yokohama as an example,the distance and time through Northern Sea Route are 6 930 nautical miles (nm)and round 18 d,but those through traditional Suez Route are 11 580 nm and round 30 d[1].As a result,around 4 650 nm and 12 d are saved.According to the Northern Sea Route Administration(NSRA),Russia,a total of 250 vessels crossed the Arctic routes between Europe and East Asia from 2011 to 2017.Although it was still a small number in the entire shipping interests,it presented an increasing tendency.On the other hand,the Arctic consists of numerous natural resources including petroleum,gas,gas hydrates and minerals.For example, it has two of the world’s major oil and gas producing areas and reserves about 20%-32%of total estimated amounts of gas hydrates in the world ocean[3]. Accordingly,more offshore and ships start to work in the Arctic,especially before the decline of international oil pricefrom 2014.

Different from open water area,ships and offshore structures would encounter not only water loads but also iceloadsin the polar region.Asa result,it is important to understand ice-water-structure interaction (IWSI),which is one of the core problems for design and operating of polar ship and offshore structures.IWSI can be seen as an extension of the fluid-structure interaction (FSI),which has been one of the hot topics in hydrodynamics for many years[4-6].Due to various mechanical properties of the ice,IWSI is more complex and difficult,which involves multi-media and multi-interfaces[7].

To solve this problem,various analytical,numerical and experimental methods are developed based on different simplifications.Analytical method including semi-analytical method has a relative longer history,along with ice model from simple to complex.A distinctive method is to solve the interaction between ice sheets and waves by using hydroelasticity theory[8-10].In this area,the researchers,especially mathematicians,from institutes of Russia,New Zealand,Australia,UK and USA,etc.contributed a lot,for example Lavrentyev Institute of Hydrodynamics SB RAS[11-14](Russia),University of Otago[15-17](New Zealand),The University of Newcastle[18-20](Australia),University College London[21-25]and University of East Anglia[26-27](UK)and Clarkson University[28-29](USA).By contrast,experimental studies,especially model tests,are more concentrated in the institutes with ice tanks,for example,Krylov State Research Centre(Russia),The Hamburg Ship Model Basin(HSVA,Germany),Institute for Ocean Technology,National Research Council(NRC-IOT,Canada),The Maritime and Ocean Engineering Research Institute (MOERI,South Korea),National Maritime Research Institute(NMRI,Japan),Aalto University(Finland),Tianjin University(China).It can be seen that ice tanks played significant roles in the development of IWSI,considering the complex mechanical properties of ice.On the other hand,numerical methods have been developing quickly in recent years with the rapid development of computer capacity.They involve various methods, including mesh methods(e.g.,finite element method,boundary element method),meshless methods(e.g.,discrete element method,smooth particle hydrodynamics)and coupled methods(e.g.,FEM-DEM,CFD-DEM).New methods(e.g.,peridynamics)start to be used in simulating IWSI related problems also.The boom in numerical simulation of IWSI forms the prime motivation of the present work.

This paper presents a review of the work on the IWSI based on numerical methods mainly.Some analytical and experimental methods are also briefly involved.The focuses are on the key techniques for solving the complex interaction among ice,water and structure, for example treatment on interfaces,and the latest applications in ship and ocean engineering.Considering the particularity of ice,we shall classify the methods according to ice simulation and analyze the applicability of different methods to different types of ice.Some work on the ice-water interaction and ice-structure collision,which provide the basis for IWSI,is also introduced very briefly.Lastly,the potential development tendencies of these methods on IWSI are predicted and analyzed.

1.Analytical/semi-analytical and boundary element methods

For some simple cases,analytical or semianalytical solutions for IWSI can be found with some assumptions.One of the most common assumptions is assuming the ice as a thin elastic plate floating on the fluid surface[16,18,30-32],as shown in Fig.1,which usually corresponds to large ice sheet or continuous icecover.

Fig.1 (Color online)Sketch of interaction between fluid and ice plate

1.1 Model formulation

The fluid is assumed to be inviscid and incompressible,and the flow of this fluid is irrotational.In this way,a velocity potential φ( x, y , z )can be introduced,which satisfies Laplace’sequation in the fluid domain

where ρWis the density of water,and g is the gravitational acceleration.For the linearized Bernoulli’s equation,the second item on the right hand of Eq.(3)can be ignored and Eq.(3)is simplified as

In addition,on the interface of lower surface of ice sheet and the water surface needs to satisfy the kinematic boundary condition and the second and third items on the right hand of Eq.(4)can be ignored for the linearized kinematic boundary condition as?φ / ?z = ? w /?t .

If the fluid depth is finite,the flat bottom should satisfy the impenetrable condition

If the ice sheet is finite,corresponding boundary conditions should also be satisfied at the ends of the ice sheet,such asfree edge condition

Equations(1)-(4)with proper boundary and initial conditions constitute a typical FSI problem,which can be solved in various methods[15,17,33-34].When a body is taken into account,the problem becomes more complicated.

1.2 Applications

Here we divide this problem into three categories,according to three different typical types of body.

1.2.1 Bottom-mounted rigid body

First category of the body corresponds to the leg of the offshore platform mainly,which can be simplified into a bottom-mounted vertical cylinder in an ice sheet,as shown in Fig.2.The cylinder can either be frozen or non-frozen in the ice sheet,and the boundary condition at the contact line between the ice sheet and cylinder surface is either clamped edge condition Eq.(6b)or free edge condition Eq.(6a)correspondingly.

Fig.2 (Color online)Sketch of a bottom-mounted cylinder in the ice sheet (revised from Ref.[26])

On the surface of the rigid body,the impenetrable condition below needsto besatisfied

where n denotes the normal vector of the body surface.

In this area,the work of Prof.Korobkin from University of East Anglia[26-27]has provided many insights.For the problem of circular vertical cylinder frozen in an infinite ice sheet,Brocklehurst et al.[26]used the Weber transform to solve the problem considering the linear boundary conditions.It was found that the amplitude of vertical force on cylinder is an order of magnitude larger than that of horizontal force.The vertical and horizontal forces reach quite large magnitude even for small-amplitude long waves.Ren et al.[23]studied multiple bottom-mounted circular cylinders in ice sheet,by using the eigenfunction expansions and the Green’s second identity in an artificial domain for each single cylinder.Both clamped and free edge conditions between the ice sheet and cylinders are considered.It was found that the horizontal force acting on the middle cylinder in an array reached a local peak in the situation that cylinders with free edge conditions at relatively thicker ice sheet.By contrast,this phenomenon was obvious in the clamped edge condition at a far thinner ice sheet.In particular,each array of cylinders can be approximated as a wall for two arrays of cylinder arranged side by side,when the distance between two arrays is much larger than the radius of cylinders.Korobkin et al.[27]used both vertical mode methods and Weber integral transform to solve the same problem as Brocklehurst et al.[26]did.They found that these two methods can obtain identical results,but the former method is simpler to solve.At the end of their work,they proposed a more practical problem to be solved,that is a cylinder is frozen in a semi-infinite ice sheet at a distance from the free edge of the ice sheet.The propagation of wave from open water to ice sheet and then to body will be solved analytically,as well asthe force induced on the body.

There is another situation,that is a body is fixed in the water area which is surrounded by infinite ice sheet.It isusually called a body standing in a polynya.Only the free edge condition Eq.(6a)needs to be satisfied at the edge of the ice sheet.Ren et al.[24]simplified the body into a cylinder and treated the polynya as a circular lake,considering the linear boundary conditions.This problem was solved by using the matched eigenfunction expansions method.It was found that the hydrodynamic coefficient oscillates much along with the radius of the polynya and the wave number.In particular,even though the radiusof polynya wasextended to be much larger than the radius of cylinder,the results did not tend to those in open water case.

1.2.2 Movable rigid body

Second category of the body corresponds to the ships in the ice lead or submersibles in the water covered by ice[35].In this category,the structure moves either in a prescribed motion or a coupling motion under hydrodynamic forces propagated from ice sheet via water,as shown in Fig.3.They are usually named forced motion and free motion,respectively[36-37].The boundary condition on the moving body satisfies

where V ( t )is known for forced motion and unknown for free motion.For the free motion,V ( t )can be solved by decoupling the fluid-body interaction motion[37-38].

Fig.3 (Color online)Sketch of a structure floating in an ice lead[22]

Due to the complex ice-water and water-body coupling,this kind of work has not been much studied and only been limited to two-dimensional problem with linear boundary conditions,to our best knowledge.Sturova[12]solved the problem of a submerged cylinder in the water covered by a finite ice sheet or two semi-infinite ice sheets.This problem was solved through the Greens function obtained by using matched eigenfunction expansions.It was found that the added mass coefficient and damping coefficient both oscillate with the wave frequency in these two situations.Recent representative work comes from the group of Prof.Wu from University College London[21-25,39-40].For example,Ren et al.[25]solved a two-dimensional floating rectangular body in a polynya confined by two semi-infinite ice sheets on two sides.They divided the fluid domain into five sub-regions,and solved this problem by using matched eigenfunction expansion method.It was also found that the body’s hydrodynamic coefficients oscillate with the wave frequency,which is caused by the radiating waves generated by the body motion being reflected back by the ice sheets.However,this matched eigenfunction expansion method was limited to the rectangular floating box.To solve thislimitation,Li et al.[21]combined the research of a floating body in open water without ice sheets and ice sheets floating on water without floating bodies.In this way,the shape of floating body was no longer restricted to be rectangular,but this method was based on wide spacing approximate,which is an approximation solution.Recently, this restriction has been fully solved by Liet al.[22].They proposed a hybrid method utilizing a simple source function and eigenfunction matching,which was used in the lead and ice sheets,respectively.

1.2.3 Moving pressure over ice/water surface

The third category corresponds to a hovercraft moving over ice surface or water surface.According to working principle of the hovercraft,it is usually simplified to a pressure load moving over ice/water surface. Thus, the pressure of the hovercraft is exerted on the corresponding equation of ice or water as an external force.

Moving loads over an infinite ice sheet have been studied for a long time. Squire et al.[30]summarized the researches of different loads moving over an infinite ice sheet, including a point load, line load and that of a hovercraft. However, the theoretical results based on Eq. (2) cannot agree well with the experimental results quantitatively. To solve this problem, different models have been developed to describe the ice sheet. One of the most commonly used is the Kelvin-Vogit model[41]

Based on this model, work of Prof. Kozin from Institute of Machine Science and Metallurgy FEB RAS[34-35,42-43]and Prof. Sturova from Lavrentyev Institute of Hydrodynamics SB RAS[11-13], et al. has contribute a lot to this problem. Kozin and Pogorelova[42]studied the resistance of a hovercraft moving over an infinite viscoelastic ice sheet steadily,by using Fourier integral method. By using the same method, Kozin and Pogorelova[43]further compared the deflection of ice sheet of different viscoelastic model, including Kelvin-Vogit model, Maxwell model and Maxwell-Kelvin model. It was found that the Maxwell-Kelvin model has better agreement with the experiment results than other two models. With respect to numerical method, boundary element method (BEM) together with finite difference method(FDM) is a popular method to solve this coupled Eqs.(1) and (2) or (9). Parau and Vanden-Broeck[44]adopted this method to solve the steady motion of loads moving on ice sheets with nonlinear boundary conditions. Instead of using Eq. (9), they added an extra item reflecting the viscosity of ice on the left hand of Eq. (2) to make sure the results convergent. Li et al.[45]also adopted this method to solve similar problem, but they solved Eq. (9) with linear boundary conditions. It was found that linear numerical results still coincided with theoretical and experiments results well.

The deflection of ice sheet under moving loads relies on the speed of the loads and exhibits different characteristics at different velocity. In general, the velocity of load is divided into three regions with respect to the critical velocity of an infinite elastic ice sheet[46-47]

where k is the wave number, H is the water depth.They are subcritical velocity, critical velocity and supercritical velocity regions.

Once the deflection of the ice sheets is obtained,the stress of ice sheet and wave-making resistance of the hovercraft can be calculated. Here we provide the equations for σ and wave-making resistance coefficient Cwas below:

where P is the pressure of the hovercraft, Ω is the area of the pressure distribution.

Fig. 4 (Color online) Deflection of ice sheet at different velocity based on BEM simulation

Figure 5 shows the stress distribution of σxxon the ice sheet at the different speed regions. According to the color bar, one can still see that the induced stress in Fig. 5(b) is larger and wider than those in Figs. 5(a), 5(c). It is easy to observe that distributions of stress in Fig. 5 are similar to those of ice deflection in Fig. 4. Once the largest of stress is over the allowable stress of the ice sheet, cracks will be generated and the ice sheets will be broken by the waves induced by the moving pressure.

Fig. 5 (Color online) Stress distribution of σxx on the ice sheet at different velocity based on BEM simulation

Fig. 6 (Color online) Wave-making resistance coefficient with different velocity, where theoretical results come from Kozin and Pogorelova[42]

On the basis of infinite ice sheet, further work studies moving loads on semi-infinite ice sheet beside semi-infinite water surface, semi-infinite water surface beside semi-infinite ice sheet, or a confined lead in two semi-infinite ice sheets, etc. For the hovercraft moving on the semi-infinite ice sheet beside semi-infinite water surface, Tkacheva[48]studied this problem by using analytical method with elastic plate model for ice sheet.The deflection of ice was investigated and similar results with those for infinite ice sheet were observed.However,different from infinite ice sheet,the wave-making resistance was found tend to zero in the subcritical velocity region.Liet al.[45]studied this problem by using BEM.However,they did not consider the edge condition Eq.(6)of the ice sheet,which was proved to be improper.For the pressure moving on the water surface,either on semi-infinite water surface or on a lead confined by two semi-infinite ice sheets,Sturova and Tkacheva[13-14]studied these two cases,respectively.It was found that the wave-making resistance oscillates near the critical velocity.Recently,Ni and Zeng[49]preliminarily studied the moving loads in an ice lead and on an infinite ice sheet with a crack by using BEM and further more work will follow on.

2.SPH-based methods

Smoothed particle hydrodynamics(SPH)is a meshless adaptive Lagrangian particle method for hydrodynamic problems.SPH and its related variants have been developing as a new generation of computational fluid mechanics solvers for decades and have been widely applied in a wide range of engineering problems,including FSI problems[50-58].More detailed developments of SPH can refer to the reviewsof Zhang et al.[59]and Gotoh and Khayyer[60].

Although SPH developsfast in hydrodynamicsas its name indicates,it can still be effectively used in solid mechanics,especially in simulating large deformation and failure behavior of solids,because of its Lagrangian nature.Recently,SPH has also been used in ice mechanics,including the damage and cracks growth process of ice[61-63].Given that SPH can be both used in simulating ice and water,it provides a good alternative method for IWSI.Besides,with the rapid development of SPH,it has been added in many commercial software already,such as LS-DYNA.

2.1 Model formulation

2.1.1 Basic equation

Without regard to energy conservation,the governing equations in SPH method are the mass and momentum conservation equations written in the Lagrangian forms[64]

where superscriptsαand βindicate the Cartesian components in x,y and z directions,ρ ,u and σare the particle density,velocity vector and stress tensor,respectively,and D/D t denotes the particle derivative.

The stress tensor can be decomposed into two parts,the isotropic pressure and viscousstress

To close the equation set,state equation is still needed for water and ice,respectively.For water,it can be assumed as a weakly compressible fluid for example and thusit follows

2.1.2 Spatial derivativesand particle approximation in SPH

In SPH method,the quantitiesof particle i can be approximated by the summation of the relevant quantities of its neighboring particles j.There are many particle approximation methods in SPH[64]and here we just introduce some common ones.The continuity Eq.(13)can be approximated asfollows

where subscripts i and j indicate the particle number,m is the mass of the particle,and Wij= W ( Ri,j, l) is the smooth function of particle j to particle i,in whichRi,jis the distance between particles i and j,and l is the smooth length.The specific expressions of Wijcan refer to Liu and Liu[64].

The momentum Eq.(14)can be approximated as follows

When solving the problem involving dissipation,one can improve the numerical stability by introducing artificial viscosityΠij[67].As a result,the SPH approximation of the momentum equation for the fluid can be written as

When SPH is applied to solid,an artificial stress method proposed by Monaghan[68],Gray et al.[69]is usually used to remove numerical instability caused by the clumping of the SPH particles.As a result, the SPH approximation of the momentum equation for ice model can be written as

2.2 Applications

Considering it is still quite new for SPH in this area,we try to review the recent developments of SPH in ice-related problems before introducing its typical applicationsin IWSIin ocean engineering.

2.2.1 Ice-structure interaction

Firstly,ice and body interaction,such as extrusion and collision,has been reviewed.There are mainly two bodies of this kind of work,one of which adopts the SPH model in the commercial software,such as LS-DYNA.In this body, the SPH method is relatively fundamental and easy to combine with other methods,but has lower accuracy.By using SPH module in LS-DYNA,Liu et al.[70]simulated an ice cone colliding a rigid plate and compared it with the experimental results to validate the ice material model.On this basis,they simulated the collision between ice sheet and the bow of an icebreaker.The deformation and damage of the ice sheet were obtained. To study the response of the structure,Gui and Hu[71]used the coupled SPH-FEM method to simulate the collision between ice and propeller.In their work,SPH was adopted to simulate ice,while FEM was adopted to simulate propeller.The contact between SPH and FEM used the point-to-face contact algorithm provided by LS-DYNA.The damage of ice under the cutting of propeller and the structural response such as stress and strain variation of propeller were both obtained.The influences of collision speed and position on ice-propeller collision were also investigated preliminarily.

Because SPH is particle-based method,the simulation is slow and takes longer time than mesh-based method.When the structural response is considered,the simulation will be much slower.To promote simulation efficiency,Bian[72]adopted SPH-FEM method to simulate the collision between a rigid conical pile leg structure and ice sheet.Different from Gui and Hu[71],Bian[72]divided the ice sheet into near-region and far-region relative to the diameter of the pile leg,as shown in Fig.7.SPH and FEM were adopted to simulate ice in near-region and far-region,respectively.The displacement and stress was continuous at the interface of near and far regions.The horizontal and vertical ice loads on the conical pile leg were obtained under various ice thickness and collision speed.This coupled SPH-FEM method has both advantagesof SPH and FEM.

Fig.7 (Color online) Interaction between ice sheet and a conical structure[72]

Fig.8 (Color online)Comparisons of the failure progress of level iceagainst sloping structure

The other body of this work chooses selfprogramming SPH method.In this body,the accuracy of SPH simulation can be improved by updated techniques, and mechanical behaviors of the ice can be simulated better.Recent representative work on this topic comes from Prof.Zheng of Harbin Engineering University[62,74-75].For example,Zhang et al.[62]simulated the bending and compression failure processes of ice with the SPH_SFDI model,which is an updated SPH method based on the simplified finite difference interpolation (SFDI)scheme.It was shown that the SPH_SFDI significantly improved the capability and accuracy for simulating ice bending and compression failures,by comparing with the original SPH scheme.On this basis,Zhang et al.[74]adopted elastic-plastic cohesion softening Drucker-Prager failure model as the ice failure model to simulate the impact process of an ice plate on the rigid inclined plane.As shown in Fig.8,the numerical results were also compared with the experimental results.The initial state as shown in Fig.8(a)including ice properties was the same as that of experiment.It can be seen that the moving tendency of the broken ice was similar to that of experimental results,but the details cannot coincide exactly.The ignorance of water effects may be one of the main reasons.Furthermore,they simulated the icebreaking process of ship in level ice and compared the average ice resistance with the experimental data from Lau et al.[76].Similarly,although the trend of ice resistance versus moving velocity was in good agreement with that in the model test,the quantities differed to some degree,as the water effects were not included.As a result,it was suggested to consider water effect in the ice-structure interaction.

2.2.2 Ice-water interaction

Secondly,the interaction of ice and water is reviewed,such as the drift and accumulation of ice floe and deformation of ice sheet in waves.In early stage,most work concerned about the accumulation of river ice in the waterway.For example,Shen et al.[77]described a two-dimensional numerical model for dynamic transport and jamming of surface ice in rivers[78].The hydrodynamic component of the model used an Eulerian finite-element method,while the ice dynamic component used SPH method.

Recently,Zhang et al.[75]have studied two typical cases of ice-wave interaction by using twodimensional SPH model.One is the motion of a small ice floe in waves and the other is wave-induced flexure of a large ice sheet,as shown in Figs.9 and 10.For the former,three degree of freedom movement of the ice floe was obtained aswell asthe drift of the ice.For the latter,the deflection and stress of the ice plate were investigated and over-wash phenomena[19]was also observed.It was found that when the incident wave was large enough,the ice plate would be broken into pieces of small ice floes.Both two cases were compared with the experiments,and good agreements achieved.At present,the work on this area is still lacking and needed to be extended.

Fig.9 (Color online)Response of small sea ice floes in regular waves[75]

2.2.3 IWSI

Fig. 10 (Color online) Wave-induced flexure of an ice floe[75]

Along with the development of ice-structure and ice-water interactions, researchers are also trying to study IWSI by using SPH. Similar to ice-structure interaction, commercial software involving SPH module and self-programming SPH algorithm are two common practices. For the former, Hu and Zhan[79]simulated the interaction among ice, water and ship by using commercial software based on FEM-SPH coupled algorithm. In their work, ship and ice were simulated by FEM shell elements and solid elements,respectively, while water was simulated by SPH. All the contact interfaces, including ship-ice interface,ship-water interface, ice-water interface and sidewall of tank were treated by using face-to-face contact algorithm. The icebreaking process of a KVLCC2 hull form was simulated preliminarily and it was found the water resistance was larger than empirical value. As they described, the complex fluid-structure interaction was simplified by using the contact algorithm, whose accuracy was much lower than that of professional CFD software.

On the other hand, by using self-programming SPH, group of Prof. Han from Harbin Engineering University[80]opened the work on IWSI, to our best knowledge. For example, Liu et al.[80]simulated the continuous icebreaking process of an icebreaker involving water effects. They adopted corrected smooth particle method (CSPM)[81]to simulate ice,water and ship and improved CSPM to extend its applications in low-speed collision based on the low-speed collision fracture material model with density correction and artificial stress correction. As water is included, they found the ventilation phenomenon during icebreaking process, which means that the fluid has no time to backfill after the broken ice is turned over, thus creating a gap between the ice and the ship[82]. Ventilation phenomenon affects the ice buoyancy and ice resistance, which is a typical difference between ice-structure interaction with and without water effects. The results are shown in Fig. 11, where Vxand Vyare velocity of fluid in x and y directions. It can be seen that numerical simulation of IWSI captured the characteristics of the continuous icebreaking process, such as the rotation and sliding of the ice fragments along the bow. Water resistance and ice resistance were calculated as shown in Fig. 12 and the total resistance was also compared with experimental results from Puntigliano[83]. It was found that the water resistance of ship in ice-covered water domain was lower than that in open water, but the ice resistance with water effects was higher than that without water effects (namely water was ignored).It denoted that the support of fluid to sea ice cannot be easily ignored when calculating the total ice resistance.

Fig. 11 (Color online) Ice-water-ship coupling calculation model based on 3-D CSPM[80]

3. DEM-based methods

Discrete element method (DEM) has become a popular method to simulate the dynamics of ice since its first application in simulation of broken ice field[84].It has been extensively used in simulating various types of ice based on various granular elements and the interaction between ice and structure[85-89]. Recent representative work on DEM for simulating ice mechanics comes from group of Prof. Ji from Dalian University of Technology[87,79-92]. However, water effects in the interaction between ice and structure are either neglected or just simplified into buoyance force and drag force exerted on ice elements[91,93].

To consider the interaction between structure and fluid, researchers start to develop DEM-based coupling methods such as DEM-computational fluid dynamics (CFD), DEM-SPH and DEM-lattice Boltzmann method (LBM), etc.[90]. Among them,DEM-CFD is more general and has been used in the simulation of IWSI already.Thus,we discuss this method mainly in thischapter.

Fig.12 Resistances of water and ice during ice-water-ship coupling processsimulation based on CSPM[80]

3.1 Model formulation

In DEM-CFD coupling method,the Euler method is used to describe fluid flow and the Lagrange method is used to describe particle movement.Fluid motion is governed by the Navier-Stokes(N-S)equation,which takes into account the effects of solid particles,and adds the void rate and momentum exchange term.On the other hand,movement of particles depends on two forces:one is the force of collision from other particles and the other is fluid force.

3.1.1 Fluid phase equation

Here we consider incompressible viscous fluid with free surface.The continuity equation and N-S equation are used as the governing equations in the whole flow field[94]

where αis the volume fraction of the fluid,if the cell is full of water,one has=1α ,u is the fluid velocity,νWis the kinematic viscous coefficient of the fluid and Rprepresents the momentum exchange between the fluid phase and the particle phase.

There are many methods to calculate the momentum exchange term Rpand a common method is introduced here[95]

where Fpis the interaction force between the fluid and particle and upis the particle velocity.More detailson Fpcan refer to Klosset al.[95].

3.1.2 Particle phase equation

where subscript i,j,k denote the particle number,m is the particle mass, Fcand Flrare the contact and non-contact forces between particles, respectively,Fgis the gravity force of the particle.I is the moment of inertia of the particle,ωis angular velocity of the particle, Mtand Mrare the momentsof the sliding friction and the rolling friction,respectively.

3.2 Coupling scheme

It is challenging to solve FSI by using different methods.For DEM-CFD,there are also many numerical difficulties.For example,time steps of two iterations in DEM and CFD are quite different,and the unification on the time scale needs to be solved.Besides,the relative size of fluid mesh is usually needed to be slightly larger than that of solid particle.With the development of coupling model, these difficulties has been gradually overcome[95,97].The coupling model is generally divided into one-way coupling and two-way coupling,where the former denotes that particles have no effect on fluid and thus saves much calculation time.The flow of two-way coupling is briefly described below.The principle of one-way coupling is similar but there is no feedback of fluid motion to particle movement.

Figure 13 shows the general framework of coupling in the DEM-CFD model[98].To start with,all components of simulation,including DEM,CFD and coupling parts,are initialized.The coupling startswith calculating the fluid porosity in each fluid cell based on the positions of particles and fluid mesh information.Thereafter,velocity of particles and fluid as well as the pressure and stress tensor at the current fluid time step are used to calculate the fluid-particle interaction force f(acting on each particle)and F(acting on fluid).The next step is the iteration loop of the DEM.The obtained data in the coupling step is used in the equations of motion of each particle Eqs.(26),(27).After the DEM loop is completed,the new position and translational and rotational velocities of all particles in the next fluid time step are obtained.On this basis,the fluid phase equations Eqs.(23),(24)of the fluid are solved.The loop ends before a new loop starts until the whole simulation terminates.

Fig.13 Framework of coupling in the DEM-CFD model[98]

3.3 Applications

DEM-CFD coupling method has been applied in many fields, such as chemical industry[99],geomechanics[100],etc.and starts to be used in ice-related area.In this section,we will try to provide the latest applications of DEM-CFD coupling method in IWSI.

Wang et al.[101]adopted DEM-CFD coupling method to simulate a ship moving in brash-ice region based on commercial software STAR-CCM.By using the DEM module,they bonded spherical particles to build two different shapes of the brash-ice model.The brash-ice elements were released into the water region and moved along the water flow at a given speed before collided with the ship hull fixed at the water surface.The streamlines of the flow and movement of the brash ice were observed,and the induced ice resistance was concerned under different speed and iceconcentration,asshown in Fig.14.

Fig.14 (Color online)Wave making at the stem (top)and its influence on the motion of crushed ice(bottom)[101]

On this basis,Li[102]used similar methods to study on the resistance characteristics of a ship moving in a brash-ice channel with flat walls on two sides,which simulates a ship moving in an ice-breaking channel.Both water resistance and ice resistance on the ship were considered,with different sizes of the ice elements.The “force chain”found in the experiment[103]was also observed in the numerical simulation,as shown in Fig.15.The characteristics of ice resistance in the channel was compared with that in an unrestricted domain.The results showed that the presence of flat walls on both sides of the channel increased the ice resistance,relative to that in an unrestricted domain.With the reduction of the ice size,the ice resistance showed an increasing trend.

Fig.15 (Color online) Comparison of a ship moving in brashice channel[102]

Fig.16 (Color online)Comparison of propeller-ice interaction between one-way(a) and two-way (b) coupling memethods[104]

To reveal the effects of one-way and two-way coupling methods,Xu et al.[104]compared them when they explored the interaction among crushed-ice,water and propeller.In two-way coupling method,it could simulate the suction effect of the propeller on the crushed ice particles better and was closer to that in the actual situation.This is because the ice particles and water affect each other iteratively in two-way coupling and the motion states of ice particles are significantly changed than those in the one-way coupling,as shown in Fig.16.As a result,the propeller-ice contact force was found about 10 times larger than that in the one-way coupling.On the other hand,the thrust and torque of the propeller were calculated in open water,one-way coupling and two-way coupling,respectively,as shown in Fig.17.There is no surprise that the results in open water and one-way coupling coincided,as the straight line in Fig.17 shows,considering that ice has no effects on fluid motion in one-way coupling method and thrust and torque both come from fluid motion.It can be seen in Fig.17 that thrust and torque oscillated fiercely in two-way coupling methods which involves the influence of ice cubes on fluid motion.The movement of fluid in front of the propeller was significantly changed by ice cubes,so the hydrodynamic performances of the propellers were different.The averaged thrust and torque were found approximately 2 times larger than those in one-way coupling and open water.

Fig.17 (Color online)Comparison of propeller thrust (a)and torque (b) between one-way and two-way coupling methods[104]

To our best knowledge,there has been little work published on the interaction between level ice and ship by using DEM-CFD coupling method.In other words,the break of ice sheet and crack development in ice sheet are very challenging in this method.Recently,this paper has preliminarily solved this problem by adding bonds between DEM particles and establishing a simple failure criterion.As shown in Fig.18,the continuous ice-breaking process of a ship moving in ice-covered water domain was simulated successfully.The numerical ice resistance was further compared with the previous model test results from Huang et al.[105].As shown in Fig.19,it can be seen that the resistance coincides well with experimental results after the whole hull enters into ice-covered water region,which validated the accuracy of the interaction among level ice,water and ship.More detailed work will follow on.

4.LB M-based methods

The LBM is a mesoscopic simulation method between macroscopic continuous simulation and microscopic molecular dynamics simulation. Since its advent in 1988[106], it has been rapidly developed into a powerful numerical simulation method in computational fluid dynamics. Now LBM is widely used in the field of two-phase flow and multi-phase flow[107,108-111]. On the other hand, LBM has also been used in the FSI problem[112-113], by solving the FSI force combined with high-precision boundary methods, such as immersed boundary method (IBM).All of these developments provide LBM a possible attempt to solve the complex IWSI, which involves abundant and complex fluid-structure interfaces.

4.1 Model formulation

4.1.1 Continuous Boltzmann equation

LBM is considered as a mesoscopic numerical model for simulating fluid motion and also a special discrete scheme for continuous Boltzmann equation.continuous Boltzmann equation describes the evolution behavior of molecular velocity distribution function f(x ,ξ,t), where x is molecular space position vector, ξ is molecular velocity vector and t is time. Under the action of external forces,molecules migrate and collide. This process is expressed as[114]

in which a is acceleration and provided by external force F=ma, the right term Ω(f) represents the change of number of molecules caused by the collision between molecules, which is also called collision term. The difficulty of solving Boltzmann equation lies in the collision term. In order to solve this problem, many methods have been proposed, such as linearization of collision term, Chapman-Enskog expansion method[114], etc.. Here we introduce a common solution method below.

4.1.2 BGK equation

In 1954, Bhatnagar, Gross and Krook (BGK)[115]proposed the Boltzmann BGK equation and introduced a simple collision operator model, as shown in the following Eq. (29). They believed that the main physical meaning of collision term is to make the distribution of molecules approach Maxwell equilibrium distribution[116]

4.1.3 Lattice Boltzmann equation

Lattice Boltzmann equation is a special discrete form of BGK-Boltzmann equation. The Boltzmann equation is discretized in both time and space.Therefore, finite dimensional velocity set{e0,e1,…,em} replaces infinite dimensional velocity ξ and fireplaces f, where i is discrete velocity direction. One has[117]

where Δt is a time step. Equation (31) is the famous Lattice Boltzmann equation. In order to obtain the macroscopical governing equation, it is necessary to combine different mesh structures and corresponding equilibrium distribution functions to derive the macroscopical governing equation.

Fig. 18 (Color online) Continuous ice-breaking process of a ship moving in ice-covered water domain

Fig. 19 (Color online) Time history of ice resistance and its comparison with model test results from Huang et al.[105]

4.2 Applications

The advantages of LBM in simulating moving bodies floating on free surface prompt researchers to use LBM to solve the complex IWSI.Jan?en et al.[118]adopted LBM to develop an efficient numerical ice tank for the simulation of interaction among rigid ice floes,water and moving ship.The fully viscous and turbulent flow fields were calculated by LBM method,and the motions of ice floes were simulated by coupling open dynamics engine (ODE)and LBM.D3Q9 velocity model[116]and single-relaxation time(SRT)model[114]were used in the collision.The volume force was added directly to the distribution function in Eq.(31).In order to improve the solution speed,they designed a set of Efficient LBM (ELBM),which included different collision operators,boundary conditions,turbulence model,interface capture method and mesh refinement technology, etc..They adopted two different ELBM models to simulate a ship moving in a crushed ice-covered water region,as shown in Fig.20.For the first one,as shown in Fig.20(a),the boundary conditions and force calculation method were updated on each lattice in each time step.While for the second one,as shown in Fig.20(b),the boundary conditions were only updated on the lattice at the interface in each time step,as the state of unchanged lattice(i.e.inside the fluid or inside the inactive gas phase)did not need to be updated.By comparison,the second one improves the computing efficiency by around 25%.

Fig.20 (Color online) ELBM applied to (a)all lattice nodes of the computational domain and applied only to (b)interface lattices[118]

Based on LBM,more and more coupled methods are also being developed to solve complex FSI.For example,DEM-LBM coupled methods is proposed to solve the interaction between solid particles and fluid[119-120].Although it hasnot been used in IWSIyet,it would be a method to be attempted.

5.Other numerical methods

Besides the methods above,we will introduce another two methods in this chapter.One is finite element method (FEM)and the other is peridynamics(PD)method. Both have been used in simulating ice mechanics and ice-structure interaction successfully,but the involvement of water has not been considered or just started.However,they both have great potentials to be extended to IWSI.Thus,we review their developmentsvery briefly in thissection.

5.1 Finite element method

FEM is a mature method to solve continuum mechanics problem.Based on energy equation or weighted residual equation,FEM solves the equilibrium equations of force and moment at any volume of structure and ice by using variational principle.Due to its advantages of convenience,practicability and effectiveness,it has been widely used in various mechanical problems,including ice-structure interaction[121-122].Work from Norwegian University of Science and Technology[122-123], etc.has provided many insightson ice-structure interaction.

Because a considerable part of ice resistance of ships comes from submerging the broken ice into water[124],it is important to consider water effects in the interaction of ice and structure.Although in FEM simulation,water is very hard to be involved directly,researchers still try to combine water effects by using simplified models.One common method is to establish a water domain in commercial software based on FEM,such as LS-DYNA.Asfor the material properties of fluids,the empty materials are used to describe the constitutive models in the software and the relationships between pressure and density are described by defining the equations of state of fluids.It should be stated that as the fluid model is not solved by Navier-Stokes equation,the hydrodynamic force cannot be obtained accurately.However,the buoyancy force of the ice can be considered through this model,which contributes to the submersion resistance as a result.Solids are coupled with fluid by defining keywords of contact,which is also an approximation method in FSI.

By using this method,Kim et al.[126]simulated the resistance performance of a cargo ship under the conditions of broken ice.The numerical results were compared with the results of the non-frozen model ice test in a water tank in South Korea and the cut ice test in an ice tank in Canada, and agreements were both achieved.Guo et al.[125]adopted the arbitrary Lagrange-Euler algorithm (ALE)in the software LS-DYNA to consider the fluid-structure coupling and calculated the resistance of an ice-going container ship under brash ice,as shown in Fig.21.The result of the simulation was compared with the experimental data of Guo et al.[127]and good agreements achieved,including the movement of brash ice and the ice resistance.On this basis,the effects of speed and ice concentration were investigated.Recently,Ni et al.[128]adopted similar method to study the maneuverability of an icebreaker in level ice.They compared the ice resistance of the ship with and without water effects and found that the existence of water increased the ice load of ship hull on all directions whether in straight or rotational motion.Besides,when water effects was considered,the oscillation degree of the ice load curve decreased,which indicated that water cushioned the collision between ship and ice and depressed the randomness of ice load.

Fig.21 (Color online)Comparisons of the ship-ice interaction between the experiments (left) and numerical simulations(right)[125]

Besides improving FEM itself,combining FEM with other methods[129-130]is also a good alternative to extend the application of FEM in ice-related problems,such as FEM-SPH for IWSI mentioned in Section 2.2.3 and FEM-DEM for structure-ice interaction[92].It has an advantage that different method can play its strength in solving different medium.The difficulty lies in the exchange of information at the interfaces of different media.

5.2 Peridynamics

Compared with other methods mentioned above,PD is a new-born numerical method which was first proposed by Prof.Silling from Sandia National Laboratories,USA,in 2000[131].PD is based on non-local formulation in an integral form,instead of solving differential equations as classical continuum mechanics does.PD applies to both continuous and discontinuous problems and is good at solving the fracture problem especially.Silling and Askari[132]further refined PD and proposed a meshfree method that discretizes the model into particles.Particles interact with particles within a finite distance.PD has great advantages in solving the spontaneous generation and propagation of cracks in materials,and does not require any additional criteria to control the direction and length of crack propagation.

Due to its unique advantages in dealing with fracture damage,PD has been developing fast and attracted more and more attentions in fracture mechanics[133-136].Although PD is new,its applicationson ice-related problemsdevelop rapidly.

Xue et al.[138]adopted bond-based PD to establish the elastic-brittle constitutive model of ice,by considering compressive and tensile strengths of the ice.They used this model to simulate the quasi-static loading problems including cantilever beam bending and three-point bending.It was found numerical results coincided with experimental data well and proved that PD was suitable for damaging process of ice under low-speed loads.On the other hand,some other researchers exerted high-speed loads on the PD ice model.Wang et al.[137]applied bond-based PD approach to investigate the icebreaking process under blast loads of an underwater explosion.As shown in Fig.22,the damage and crack process of an ice plate were captured.Numerical results were compared with field experiments and qualitative agreement was achieved.On this basis,factors affecting the radius of breaking ice were investigated,such as the depth and massof the underwater explosive,and the thicknessof the ice.Similarly,Niet al.[139]used PD to simulate the damage process of an ice plate under high-pressure bubble jet load.BEM was used to simulate bubble dynamics and jet load while PD was used to simulate ice mechanics.Furthermore,PD has also been used in the interaction of ice and with various structure,including vertical structures[140],a sloped wall[141],and a moving ship[142].

Fig.22 (Color online) Crack propagation of an ice plate under blast loads[137]

Similar to other meshless method like DEM and SPH mentioned above,the computational efficiency of PD is still incomparable with mesh method like FEM,although the PD algorithm is being optimized.To overcome this problem,Lu[143]combined PD in the crack development region (near-field of structure)with FEM in the non-crack development region(far-field of structure),as shown in Fig.23.In this way,the generation and development of ice crack has been simulated with high accuracy and the total calculation time has been saved.On the other hand,if water effect needs to be considered,it is an alternative to combine different methods.Recently,Liu et al.[144]combined the bond-based PD for solids and the updated Lagrangian particle hydrodynamics[145-146]model of fluids to simulate ice-water interaction.They studied the impact of a rigid ball on an ice plate floating in water by using this combined method.It found that the main characteristics of the dynamic ice-breaking process were captured successfully,which proved the feasibility of this combined methods in applying to IWSI problems.More work will follow on to develop specific techniques.

Fig.23 (Color online)The coupling between PD and FEM to simulate ice breaking of a ship[143]

6.Experimental methods

Experimental method is the most direct method to study IWSI as we know.However,considering the particularity of ice,it needs relatively harsh environment to do the experiments,either for in-suit test or model test.In this section,we will introduce the development of experimental methods for IWSI very briefly,given that the emphasis of this paper is numerical methods.

In-suit test is believed as the most reliable method and a series of empirical formulas of ice loads were put forward based on the in-suit test results[147-151],which plays important roles in preliminary design and evaluation of polar structures.However,considering the various ice conditions,the data acquisition methods,the structural geometry and other factors,the in-suit results from different literatures may have discreteness and applicability of some empirical ice-load formulas is quite limited.For example,Bjerk?s[152]collected 31 public available full scale ice loads to fixed structures reported from Europe,America and Asia.He grouped all the data into 6 groups and plotted the upper bound of ice pressure versus structural width.On this basis,an approximate upper limit of ice forces could be predicted quickly in the early phase of estimation.Kellner et al.[153]examined the applicability and quality of prediction of empirical formulas of ice loads.They found that most formulas overestimated the ice loads and applicability of some empirical formulas was questionable.To avoid these deficiencies of deterministic methods,probabilistic approaches to analyze the in-suit data have been developed[154-157].Correspondingly,risk-based design and assessment of polar structures are under developing[158-160].

On the other hand,model test is an effective experiment method.Different from model test of other problems,ice should be simulated appropriately first for model test of IWSI.It is extremely hard to scale real ice to model ice because of the very complex material properties of ice[161].There are usually two kinds of model ice used widely[162]:non-frozen model ice (NMI)and frozen model ice (FMI).NMI does not need extra freezing equipment,which is mainly used in open water tank,as shown in Fig.24.Various materials and their mixture were used to prepare NMI,such as paraffin[163],parget[164], polyethylene[165],and polypropylene[166],etc..On the contrast,FMI is prepared and used in ice tanks[167],asshown in Fig.25,some representatives of which have been listed in the introduction.It is generally recognized that there are three generations of FMI, according to the time of advent:saline-doped ice[168],carbamide (urea)-doped ice[169]and EG/AD/S ice[170].All these three generations of FMI are still used in different ice tanks.Different ice tank has different recipes,even for the same generation of FMI,which usually becomes the secret or patent of the ice tank.Until now,none of FMIcan fully satisfy the scaled mechanical properties and new preparation methods of FMI with better mechanical properties are always under researched,such asmixing micro-bubbles[171].

Fig.24 (Color online) A model ship moving in the NMI in the outdoor ice tank of Harbin Engineering University,China[172]

Fig.25 (Color online)A model ship moving in the FMI generated by the Ice tank of Tianjin University,China[105]

With various NMI and FMI,researches have studied the IWSI for a long time.From the perspective of structural style,we can roughly divide two categories.First category is offshore platform.Because spud leg is the structure which encounters with ice directly,most researchers studied spud leg structure,which is also the preventive structure form for offshore wind turbine,bridge,lighthouse,etc.Ice-induced vibration is a major concern.Three regimesof ice-induced vibrations have been observed:intermittent crushing,frequency lock-in and continuous brittle crushing[173].However,the mechanical mechanism of ice-induced vibrations has not reached a consensus yet[174]. Some experts thought it as a compulsive vibration[175-176]but others thought it as an ice-structure interaction process[177-178]. Model experiments in ice tank tried to explore the mechanism of ice-induced vibration and have obtained some insights[174,179-181].Optimized designs to reduce the ice loads are always under research.For example,a most practical way is to add ice-breaking cones on conventional vertical legs[182-184],which can convert crushing failure mode to bending failure mode.

Second category is ship structure,which involves various structural styles.There are mainly two bodies of model experimental work on this category.One is adopting NMI in open water to study the interaction of ship and broken ice.Ice resistance except of icebreaking component can be obtained as well as the motions of broken ice under ship collisions or wave motions[126-127,185-186].The other is using FMI in ice tanks.Froude number and Cauchy number proposed by Timco[187]were usually adopted as similarity criterions.It can simulate various ice conditions,including level ice,crushed ice, ice floe,snow-covered ice,ice ridge and iceberg[105-190].Local ice pressure and total ice resistance have been obtained by using discrete pressure sensors, distributed tactile sensor,and/or force sensors[155,191]. On this basis,icebreaking modes[188],maneuverability[192]and propulsion performance[193]of the ship have also been studied.The results have played important roles in the design of polar ships.One challenge is how to return the model results to prototype results as accurately as possible,considering the complexity of ice properties.

7.Discussionsand development tendencies

With rapid development of numerical algorithms,IWSI has been developing fast recently.Different method is applicable for different situations.In this section,their applicability is discussed and their development tendenciesare predicted briefly.

For analytical or semi-analytical solutions,the ice is usually taken as large ice sheet or continuous ice cover,which cannot be broken,and the structures are simplified into simple and regular rigid bodies.Until now,there has been little work related to structural response of deforming body.Waves can be involved as the driving force.One of the tendencies is to develop more complex ice model based on elastic plate and viscoelastic plate,such as porous elastic plate[20]and porous viscoelastic plate[194].These ice models close to real ice more and thus will provide more insights into the ice-wave interaction as well as IWSI in the future.

For SPH,it can not only simulate level ice but also brash ice.Because SPH isgood at simulating FSI,it is easiest to extend IWSI in all the numerical methods mentioned here.At present,SPH still needs to be improved in simulating ice dynamics.On the other hand,as it is particle-based method,its computational efficiency is still low and limits its engineering applications.An alternative approach isto combine them with mesh-based method,such as FEM.

For DEM,it is one of the mainstream approaches in ice mechanics,which is suitable for various types of ice,such as level ice,brash ice,ice ridge and iceberg,etc..However,DEM is incapable of simulating fluid dynamics,which limits its extension to ice-water interaction.A common practice is combining DEM with other methods capable of simulating hydrodynamics,such as CFD,SPH,LBM,etc..Among them,DEM-CFD ismost popular and has been used in engineering problems already.Furthermore,DEM can also be combined with FEM to solve structural responseof deforming body.

For LBM,it is a computational fluid dynamics method in essence.Thus,it can just be applied to simulate brash ice without considering break.Its advantage lies on its treatment on the interaction between fluid and multi-bodies.Without solving the interface between fluid and bodies directly,it is good at solving interaction between large-scale fluid and large-number floating bodies.It provides a possible method for numerical ice tank in the future.

For FEM,similar to DEM,it is also one of the mainstream approaches in ice mechanics.However,it is mesh-based method and not quite suitable for the calculation of discontinuous medium and the generation of cracks and subsequent cracks,which limits its application in icebreaking problems.One alternative approach is to improve FEM based on fracture mechanics,such as cohesion element method(CEM)[195]and extended finite element method(XFEM)[196],etc.In order to considering water effect,a common approach is coupling with hydrodynamic methods,such as CFD,SPH,etc..

For PD,it is a new generation of numerical method for fracture mechanics.It has been used in icebreaking simulation and ice-structure interaction.As it is new,it still needs time to improve its algorithm and to apply in engineering problem.For example,similar to SPH,to combine with FEM in space can enhance its computational efficiency.Recently,an Eulerian PD model for compressible fluids has been developed[134], which may provide a possible way to simulate ice-water interaction via a framework combining Lagrangian and Eulerian PD models.

Lastly,for experimental studies,it is still the most reliable method to study IWSI.Experimental data have a relatively large discreteness due to several factors.Probabilistic approaches to solve experimental data have been developed and risk-based design and assessment of polar structures are recommended accordingly.New preparation methods of model ice with mechanical properties closer to real ice are alwaysneeded.Based on model experimental methods,new icebreaking or ice-load reducing methods are developing rapidly, such as high-pressure bubble icebreaking method[197-200],air bubbling system[172,201],etc..