Numerical Simulation of an Air-Bubble System for Ice Resistance Reduction

: Ships sailing through cold regions frequently encounter ﬂoe ice ﬁelds. An air-bubble system that reduces friction between the hull and ice ﬂoes is thus considered useful for the reduction of ice-induced resistance. In this study, a numerical analysis procedure based on coupled ﬁnite volume method (FVM) and discrete element method (DEM) is proposed to simulate complicated hull-water-gas-ice interactions for ice-going ships installed with air-bubble systems. The simulations reveal that after turning on the air-bubble system ice ﬂoes in contact with the hull side wall are pushed away from the hull by the gas-water mixture, resulting in an ice-free zone close to the side hull. It is found that the drag reduction rate increases with the increase of ventilation, while the bow ventilation plays a deciding role in the overall ice-resistance reduction. The proposed procedure is expected to facilitate design of new generations of ice-going ships.


Introduction
As a consequence of global warming, the polar marginal ice zone has been observed to become wider [1].The polar marginal ice zone is characterized by floe ice fields that cover 15-80% sea surface [2].The existence of ice in water induces extra resistance for ships sailing in floe ice fields [3][4][5].During the hull-ice interaction, significant kinetic energy from ship propulsion is dissipated, resulting in speed loss or additional fuel consumption [6,7].An air-bubble system that reduces the friction between the hull and the ice floes is thus considered useful for the reduction of ice-induced resistance in floe ice fields.
The air-bubble system discussed in this article has something in common with but differs from the air lubrication technology that has been utilized for drag reduction.The ship's air lubrication technology can be subdivided into two main categories.The first group is termed bubble-induced skin-friction drag reduction (BDR), for which a large number of micro-bubbles are injected into the boundary layer.The second group is termed air layer drag reduction (ALDR), which forms a continuous gas layer on the hull surface [8].Both groups utilize air as a lubricant, which has been proven to decrease the friction between the ship and the seawater, see, e.g., in [9,10].The air-bubble system of this study is more similar to the former group of air lubrication technology, but aims to reduce hullice friction instead.This is achieved by injecting air from a series of nozzles at the bow and bilge.When the air bubbles arise along the hull, the mixed air and water create a strong current, forming a layer between the hull and the ice floes, consequently reducing ice resistance.
The idea of reducing ice resistance by air bubbles originated in the late 1960s [11].An air-bubble system was studied mainly through model tests.Till the early 1990s, several icebreakers and ice-going vessels that mainly operate in the Baltic Sea were installed with air-bubble systems [12].Little progress in air-bubble systems for ice resistance reduction has been reported in recent decades.Nevertheless, as more ships sail into the polar floe ice fields, an air-bubble system may bring added value to ice-going ships' efficiency and safety, and thus deserves further investigation.Furthermore, the fast development of computational fluid dynamic (CFD) methods makes it possible to investigate complicated hull-water-air-ice interaction processes with numerical simulations, instead of depending on high-cost model tests.In this work, the authors aim to make use the state-of-the-art numerical methods to quantify the ice resistance characteristics of ice-going ships installed with air-bubble systems.
There are two major numerical methods for simulating air-bubbles.The first one is termed the interface tracking method [13].With this approach, the fluid interface can be accurately defined.However, this method requires a complicated process of mesh reconstruction, and it has the disadvantage of mass and energy loss of the bubbles.The second method is called the interface capturing method [14], in which the fluid interface does not have to be accurately defined.The different liquids are instead distinguished through additional fluid variables such as the mass fraction.The interface capturing method requires a large number of grid cells to keep the accuracy.This method is represented by the method of the volume of fluid (VOF) approach, which is often employed in the simulation of large bubble motions and free surface flow in liquids [15].Some recent research making use of the VOF approach are as follows: Zhu et al. [16] used the VOF method to investigate the effects of gas velocity, liquid velocity and other factors on the bubble detachment diameter.Based on the VOF method, Li et al. [17] described the deformation during the ascent of a single bubble in gas-liquid, gas-liquid-solid multiphase flow under high pressure.Tsui et al. combined the VOF method with the Level Set method to simulate rising bubbles in still water and got agreeable results [18].
Numerical simulations of ship resistance characteristics in ice-infested waters have been presented by many researchers.A recent review paper by Li and Huang [19] indicates that the discrete element method (DEM) predicates reasonable resistance induced by broken ice.Hansen and Loset [20] applied a two-dimensional DEM model to simulate the ice force on ships in broken ice.Ji et al. [21,22] used the GPU parallel algorithm to accelerate the DEM calculation, making it possible to use the DEM method to calculate sea ice structure interaction in the large-scale calculation domain.Luo et al. [23] applied combined CFD-DEM to study the coupling characteristics of ship-ice-water in the brash ice channel and analyzed the difference between one-way coupling and two-way coupling.In addition to the DEM approach, other numerical methods were also employed for ship-ice interaction simulations.Kim et al. [24] simulated ice resistance in ice channels through the finite element method (FEM) and found the simulation results were in good agreement with the model test results.Lubbad and Loset [25] simulated ships in ice through the physics engine PhysX, and compared it with the full-scale measurements.Furthermore, there are other emerging numerical methodologies, such as the Peridynamics (PD) method, the Smooth Particle Hydrodynamics (SPH) method, and the Extended Finite Element (XFEM) method.All those numerical methods have the potential to simulate ship-ice interactions with reasonable accuracy [19,26].
In this paper, the authors utilized the combined CFD-DEM approach to simulate ship resistance in floe ice fields with the air-bubble system installed.The air-water interface was simulated by using the VOF method.An icebreaker was chosen as the case study vessel.By this means, we aim to find out how effective the air bubbles are for drag reduction in ice-infested waters.

The Numerical Models
In this section, the features of the numerical models of this study are summarized and the key theoretical formulations are presented as follows.

The Governing Equations of the Incompressible Fluid
In this study, the finite volume method (FVM) is used to discretize the fluid domain.The governing equations are: ∂ ∂t where t is time, V is the volume of the fluid element, a is the area vector, v is the velocity vector of the fluid element, S u is the source term of the continuity equation, p is the pressure, T is the viscous stress tensor, f b is the resultant force of the body force, s u is the source term of the momentum conservation equation.The viscous stress tensor can be expressed as: where µ is the dynamic viscosity coefficient, I is the unit tensor.

The Governing Equations of the Discrete Phase in Numerical Simulation
The governing equations of the discrete phase follow the Lagrangian framework.The surface force and physical force acting on the particle jointly determine the change of particle momentum, and its momentum conservation equation is: In this equation, m p is particle mass, v p is particle velocity, F s is surface force, F b is body force, F d is the drag force, F p is pressure gradient force, F vm is virtual mass force, F g is the gravity, F con is the contact force.The two most critical items are the calculation of the drag force and the contact force.The former involves the treatment of the gas-liquid mixed phase, and the latter depends on the choice of the contact model.The calculation of both will be introduced in subsequent sub-sections.
The conservation of angular momentum of the particle can be expressed as: where I p is the moment of inertia of the particle, ω p is the angular velocity of the particle, M b is the resistance moment, r c is the vector from the contact point to the center of gravity, F ci is the contact force between particle c and particle i, and M ci is the moment of rolling resistance acting on the particle.

The Turbulence Model and the Free Surface Treatment
The turbulence model of the RANS equation used in the numerical simulation of this paper is the standard k-ε model.The existence of the turbulence model is to make the Reynolds-averaged Navier-Stokes equation closed.For the RANS equation, the average value can be regarded as the time average of the steady-state situation and the overall average of repeatable transient situations.Inserting the decomposed solution variables into the Navier-Stokes equations yields equations of average quantities.The conservation equations of average mass and average momentum can be expressed as: where ρ is the density, v is the mean velocity, p is the mean pressure, T is the mean viscous stress tensor, and f b is the resultant force of body forces (such as gravity, centrifugal force, etc.).This equation is essentially the same as the original N-S equation, with an extra term T RANS added to the momentum equation, which is the stress tensor.
The standard k-ε model is a two-equation model that determines the turbulent length and time scale by solving two independent transport equations.This model assumes a fully turbulent fluid flow and does not take into account the effects of molecular viscosity.The transport equation corresponding to the turbulent kinetic energy and dissipation rate of the standard k-ε model is of the form: where G k is the term from the turbulent kinetic energy, k due to the average velocity gradient, G b is the term from the turbulent kinetic energy caused by the buoyancy effect, C 1ε , C 2ε and C 3ε are empirical constants.σ k and σ ε are the Prandtl numbers corresponding to the turbulent kinetic energy and dissipation rate, respectively.S k and S ε are user-defined source terms.
The relationship between the turbulent kinetic viscosity and turbulent kinetic energy and the dissipation rate can be expressed as: In this equation, C µ is the empirical constant.The standard k-ε model is a semiempirical formula derived from physical experiments combined with theory.
In order to capture the water-air interface better to simulate the effect of the bubble assist system, this paper adopts the VOF (Volume of Fluids) method and HRIC (High-Resolution Interface Capture) format.The VOF method is used to capture incompatible terms and assumes that the mesh resolution is sufficient to resolve the position and shape of the interface between the different phases.Therefore, we should pay attention to the mesh size during the numerical simulation.Figure 1 shows the unsuitable mesh and the suitable mesh:  The VOF model describes the phase distribution and position of the interface through the field of the phase volume fraction α i , the volume fraction α i = V i V of the phase i, where V i is the volume of phase i in the grid cell, V is the volume of the grid cell, And the sum of the volume fractions of all phases within each grid cell is 1.When the grid contains only a single fluid, the material properties of the fluid are used in the calculation.If there are multiple fluid phases in the grid, it is regarded as a mixture, and its material properties use the weighted average of each phase.
The distribution of fluid phase i is determined by the mass conservation equation: where a is the surface area vector, v is the velocity of the mixed fluid, v d,i is the diffusion velocity, S αi is the source term of the phase i, and Dρ i /Dt is the Lagrangian derivative of the phase density ρ i .When only two phases, water and air, are present in the simulation, the mass conservation equation is solved for the first term only, and the volume fraction of the second phase is adjusted in each grid cell so that the sum of the volume fractions equals 1.
The momentum equation of the fluid can be expressed as: where p is the pressure, I is the unit tensor, T is the stress tensor, f b is the vector of the body force, S i α is the momentum source term of the phase, and g is the gravity acceleration.The drag force F d provided by the mixed fluid can be calculated as: where C d is the drag coefficient, ρ is the density of the continuous phase (mixing density for multiphase flow), v s = v − v p , v is the instantaneous velocity of the continuous phase, v p is the particle slip velocity, and A p is the projected area of the particle.The drag coefficient C d in this equation is determined by the Schiller-Naumann correlation, which applies to fluids with bubbles, which is set as: Re p (1 + 0.15Re p 0.687 ), Re p ≤ 1000 0.44, Re p > 1000 (14) where Re p is the particle Reynolds number, which is defined as Re p ≡ , where D p is the particle equivalent diameter and µ is the kinematic viscosity.

The Ice Model
Ice is modeled by using DEM, and the governing equation is Newton's law which has been described in detail in Section 2.2.The contact force in the governing equation needs to be calculated by the contact model.The contact model of DEM will be introduced below.
We employ the DEM method in this study to model the ice floes.Two major assumptions are taken to simplify the computation.Firstly, we assume that the ice floes will be pushed away but not broken during the ship-ice interaction process.Secondly, the contacts of ship-ice and ice-ice are assumed as elastic.Based on these assumptions, the Hertz-Mindlin model can be implemented.In this model, the spring simulates the elastic part of the collision process, and the damper reflects the energy dissipation of the collision process.The contact force between two DEM particles is described by the following equations: where F con is the contact force, F n is the normal force, F t is the tangential force.The normal force is expressed as: where K n is the normal spring stiffness, d n is the overlap of the local normal directions of the contact between the two particles, N n is the Normal damping, v n is the normal velocity of the particle.
The normal spring stiffness is: The normal damping is: where E eq is the equivalent Young's modulus, R eq is the equivalent radius, M eq is the equivalent particle mass and The tangential direction is defined by: where K t is the tangential spring stiffness, N t is the tangential damping, d t is the overlap of the local tangential directions of the contact between the two particles, C fs is the coefficient of static friction.The tangential spring stiffness is: The tangential damping is: where The equivalent radius is: The equivalent particle mass is: The equivalent Young's modulus is: The equivalent shear modulus is: where M A and M B are the masses of spheres A and B; R A and R B are the radii of the sphere; E A and E B are Young's modulus of the sphere; ν A and ν B represent Poisson's ratio of A and B, respectively.For collisions between particles and walls, the above formula remains the same but assumes that the wall radius and mass are summed R wall = ∞ and M wall = ∞, so the equivalent radius decreases to R eq = R partical , and the equivalent mass decreases to M wall = M partical .

The Computational Domain Settings
The main ship particulars of the case study vessel are listed in Table 1. Figure 2 illustrates the three-dimensional hull model as well as the location of the nozzles of the air-bubble system.It is noticeable that the nozzles are placed on the bow and along the side instead of in the bottom and the keel areas.This is because conventional air-bubble systems aim at reducing the water resistance of the ship, which requires the air-bubble system to cover the wet surface of the hull as much as possible.For that purpose, the bottom and the keel areas need to be covered by air-bubbles.The air-bubble system in this study, in contrast, is supposed to reduce ice resistance instead.It would be sufficient to use a smaller volume of air from the bow/sides to push the crushed ice away from the hull.The numerical modelling and simulations are carried out in the commercial CFD software STAR-CCM+.In the calculation domain, the stern bottom is considered the coordinate origin.The positive X-direction is defined as the direction from the stern to the bow; the positive Y-direction is to the port side; the positive Z-direction is upwards.In order to minimize the influence of the boundary on the flow field, the fluid domain should be set as large as possible.In this study, following the experience of a previous study [27], the distances between the boundary from the stern and the bow are set to be twice the ship's length.The water depth is also set to be twice the ship's length.The width of the ice field is set to be three times the ship's breadth.The ice field is modelled with ice floes made of DEM elements.The size of the ice floe is assumed as 4 m × 4 m × 1 m.The ice concentration is controlled by setting the distance of the injecting point as the injecting speed of the DEM element.The computational domain with a concentration of 60% is illustrated in Figure 3.
software STAR-CCM+.In the calculation domain, the stern bottom is considered the coordinate origin.The positive X-direction is defined as the direction from the stern to the bow; the positive Y-direction is to the port side; the positive Z-direction is upwards.In order to minimize the influence of the boundary on the flow field, the fluid domain should be set as large as possible.In this study, following the experience of a previous study [27], the distances between the boundary from the stern and the bow are set to be twice the ship's length.The water depth is also set to be twice the ship's length.The width of the ice field is set to be three times the ship's breadth.The ice field is modelled with ice floes made of DEM elements.The size of the ice floe is assumed as 4 m × 4 m × 1 m.The ice concentration is controlled by setting the distance of the injecting point as the injecting speed of the DEM element.The computational domain with a concentration of 60% is illustrated in Figure 3.  Table 2 lists the boundary conditions of the computational domain.The boundary to the right is set as the velocity inlet, while the boundary to the right is the outlet.The other surrounding boundary surfaces are set as slip wall conditions.The hull surface and the sides of the ice field are set as non-slip wall conditions.The air inlet is set as the velocity inlet boundary.For the computational domain and ship model, the fluid domain meshes with the trimmed meshing model and the boundary layer grids are divided around the ship's surface.The y+ value is in the range of 30-60.The meshes on the ship's surface are refined near the waterline, the stem, the stern, and around the nozzles as shown in Figure 4.The total mesh number is about 3.9 million.
The numerical simulation is carried out by discrete solution of N-S equation based on finite volume method, and the multiphase flow model adopts the VOF method to realize interface tracking [28].In this paper, there are two kinds of fluids in the computational domain, α 0 represents the volume function of the air phase and α 1 represents the volume function of the water phase.In the computational domain, the sum of the volume fractions of the two phases is 1 (α 0 + α 1 = 1).The turbulence model adopts the standard k-ε Model [29].
For the computational domain and ship model, the fluid domain meshes with the trimmed meshing model and the boundary layer grids are divided around the ship's surface.The y+ value is in the range of 30-60.The meshes on the ship's surface are refined near the waterline, the stem, the stern, and around the nozzles as shown in Figure 4.The total mesh number is about 3.9 million.The numerical simulation is carried out by discrete solution of N-S equation based on finite volume method, and the multiphase flow model adopts the VOF method to realize interface tracking [28].In this paper, there are two kinds of fluids in the computational domain, 0  represents the volume function of the air phase and 1  represents The nozzle on the hull surface is set as the velocity inlet of the air.In this study, the gas velocity and pressure of the air-bubble system are small that the compressibility is ignored.The separation flow model is employed to describe the liquid phase and gas phase, which solves the momentum equation corresponding to each dimension, and associates the momentum equation with the continuity equation through the prediction correction method.The second-order discretization of the convective flux is in use, which is deemed particularly suitable for constant-density fluids.Each nozzle contains at least 25 complete meshes to avoid excessive numerical loss when the computational domains are generated.Figure 5 illustrates the numerical losses of the computational domains under different mesh numbers.The color in the figure represents the volume fraction of air, red represents the gas volume ratio of 100%, yellow represents the gas-liquid interface (gas volume ratio is 50%), and green color represents the gas attached to the surface of the hull.ar.Sci.Eng.2022, 10, x FOR PEER REVIEW 10 of the volume function of the water phase.In the computational domain, the sum of the v ume fractions of the two phases is 1 ( 0 ).The turbulence model adopts standard k-ε Model [29].
The nozzle on the hull surface is set as the velocity inlet of the air.In this study, gas velocity and pressure of the air-bubble system are small that the compressibility ignored.The separation flow model is employed to describe the liquid phase and g phase, which solves the momentum equation corresponding to each dimension, and sociates the momentum equation with the continuity equation through the prediction c rection method.The second-order discretization of the convective flux is in use, which deemed particularly suitable for constant-density fluids.Each nozzle contains at least complete meshes to avoid excessive numerical loss when the computational domains generated.Figure 5 illustrates the numerical losses of the computational domains und different mesh numbers.The color in the figure represents the volume fraction of air, r represents the gas volume ratio of 100%, yellow represents the gas-liquid interface (g volume ratio is 50%), and green color represents the gas attached to the surface of the hu The ice floes are modelled using the DEM element, following the theoretical mod as described in Section 2.4.In this work, the ice density is set as 900 kg/m 3 , Young's mo ulus is assumed as 1 GPa and Poisson's ratio is assumed as 0.3.

Comparision of Simulations
Simulations with the air-bubble system activated are compared with those witho the air-bubble system under identical operational conditions.For the air-bubble syste the air inlet velocity is 2.5 m/s, and the air jet direction is perpendicular to the hull surfa The ice floes are modelled using the DEM element, following the theoretical models as described in Section 2.4.In this work, the ice density is set as 900 kg/m 3 , Young's modulus is assumed as 1 GPa and Poisson's ratio is assumed as 0.3.

Comparision of Simulations
Simulations with the air-bubble system activated are compared with those without the air-bubble system under identical operational conditions.For the air-bubble system, the air inlet velocity is 2.5 m/s, and the air jet direction is perpendicular to the hull surface.The ship speed is 6 knots.Figure 6 illustrated the wave patterns and streamlines around the hull for the cases of when the air-bubble system is deactivated and when the air-bubble system is activated, respectively.When the air-bubble system is turned on, the wave pattern around the hull is found to differ obviously from that when the air-bubble system is turned off.After the gas is pumped out from the nozzles, the gas-water mixture rises along the side of the hull and then the gas escapes from the free surface, resulting in a more distorted free surface around the hull.It is also observed that when the gas-water mixture exists the streamlines differ from those when there is no gas in the fluid flow.When ice exists in the water, additional resistance is induced by the interaction between the hull and the ice floes.In this study, the ice fields are assumed to be composed of ice floes of identical size and shape of 4 m × 4 m × 1 m.The ice concentration is assumed to be 60%.Figure 7 illustrates the snapshot at 120 s of the simulation in the ice field with a ship speed of 6 knots.The air-bubble system has not been activated.When a ship enters the floe ice fields, the speeds of the ice floes around the bow drop rapidly due to wavemaking and the collision with the bow.Consequently, the ice floes accumulate around the bow and some of them slide then along the bow area to the ship's sides and the bottom, as shown in Figure 7.It is observed that during the hull-water-ice interaction, the ice floes are overturned by the wave system rising from the bow/shoulder area.Some of the ice floes collide with the bow and also with other ice floes.Then some ice floes move along the ship's side or the bottom, resulting in friction forces on the hull.An ice-free channel slightly narrower than the width of the ship is formed behind the ship.When ice exists in the water, additional resistance is induced by the interaction between the hull and the ice floes.In this study, the ice fields are assumed to be composed of ice floes of identical size and shape of 4 m × 4 m × 1 m.The ice concentration is assumed to be 60%.Figure 7 illustrates the snapshot at 120 s of the simulation in the ice field with a ship speed of 6 knots.The air-bubble system has not been activated.When a ship enters the floe ice fields, the speeds of the ice floes around the bow drop rapidly due to wave-making and the collision with the bow.Consequently, the ice floes accumulate around the bow and some of them slide then along the bow area to the ship's sides and the bottom, as shown in Figure 7.It is observed that during the hull-water-ice interaction, the ice floes are overturned by the wave system rising from the bow/shoulder area.Some of the ice floes collide with the bow and also with other ice floes.Then some ice floes move along the ship's side or the bottom, resulting in friction forces on the hull.An ice-free channel slightly narrower than the width of the ship is formed behind the ship.
In contrast, when the air-bubble system is turned on, the hull-water-ice interaction becomes significantly different, which is illustrated in Figure 8.It is observed that when the air-bubble system is on, despite ice accumulation remaining unchanged around the bow, much fewer ice floes become in contact with slide through the shoulder due to the gas-water mixture.Most of the ice floes are overturned before passing through the ship's shoulder and drifting away from the hull, which greatly reduces the hull-ice contact occurrence at the sides and the bottom.
a ship speed of 6 knots.The air-bubble system has not been activated.When a ship enters the floe ice fields, the speeds of the ice floes around the bow drop rapidly due to wavemaking and the collision with the bow.Consequently, the ice floes accumulate around the bow and some of them slide then along the bow area to the ship's sides and the bottom, as shown in Figure 7.It is observed that during the hull-water-ice interaction, the ice floes are overturned by the wave system rising from the bow/shoulder area.Some of the ice floes collide with the bow and also with other ice floes.Then some ice floes move along the ship's side or the bottom, resulting in friction forces on the hull.An ice-free channel slightly narrower than the width of the ship is formed behind the ship.In contrast, when the air-bubble system is turned on, the hull-water-ice interaction becomes significantly different, which is illustrated in Figure 8.It is observed that when the air-bubble system is on, despite ice accumulation remaining unchanged around the bow, much fewer ice floes become in contact with slide through the shoulder due to the gas-water mixture.Most of the ice floes are overturned before passing through the ship's shoulder and drifting away from the hull, which greatly reduces the hull-ice contact occurrence at the sides and the bottom.

Ice Resistance Calculation
In this subsection, we quantify ice-induced resistance for the cases with and without the air-bubble system.The top subplot of Figure 9 shows the time history of the total ice resistance when the air-bubble system is on, for which t = 11.9 s is the timestep when the ship bow reaches the ice field.It is observed that the ice resistance gradually increases and becomes stable around t = 20 s, which corresponds to the timestep when the entire hull has entered the ice field.The middle and bottom sub-plots of Figure 9, illustrate the time series of the ice resistance components from the bow and the ship's sides, respectively.The ice resistance values of the stable stage, i.e., 20-120 s are listed in Table 3.The bow area accounts for a major part of the total ice resistance.

Ice Resistance Calculation
In this subsection, we quantify ice-induced resistance for the cases with and without the air-bubble system.The top subplot of Figure 9 shows the time history of the total ice resistance when the air-bubble system is on, for which t = 11.9 s is the timestep when the ship bow reaches the ice field.It is observed that the ice resistance gradually increases and becomes stable around t = 20 s, which corresponds to the timestep when the entire hull has entered the ice field.The middle and bottom sub-plots of Figure 9, illustrate the time series of the ice resistance components from the bow and the ship's sides, respectively.The ice resistance values of the stable stage, i.e., 20-120 s are listed in Table 3.The bow area accounts for a major part of the total ice resistance.
As mentioned previously, when the air-bubble system is turned on with the air injection rate of 2.5 m/s, the hull-ice contact on the ship side is greatly reduced.Figure 10 illustrates the time series of the ice resistances for this case.In comparison with the ice resistances in Figure 9, the total resistance as well as the components from the bow area and the sides are found to be smaller.The resistance values are listed in Table 3, together with the resistance when the air-bubble system is off.The resistance reductions are also included in Table 3.It is seen that when the air-bubble system is off, the total ice resistance is reduced by 15.3%.When it comes to the ice resistance components, the bow resistance remains almost unchanged, with a reduction rate of 10.3%.In contrast, the ice resistance from the ship's sides is greatly reduced, with a drag reduction rate of 70.8%.
ship bow reaches the ice field.It is observed that the ice resistance gradually increases and becomes stable around t = 20 s, which corresponds to the timestep when the entire hull has entered the ice field.The middle and bottom sub-plots of Figure 9, illustrate the time series of the ice resistance components from the bow and the ship's sides, respectively.The ice resistance values of the stable stage, i.e., 20-120 s are listed in Table 3.The bow area accounts for a major part of the total ice resistance.The effect of the ship speed on the drag reduction rate was also investigated.In addition to the abovementioned speed 6 knots, ship-ice interactions under four other speeds were simulated and ice resistances with and without the air-bubble system were compared.Table 4 listed the resistances as well as the drag reduction rates, which are also plotted in Figure 11.It is observed that the drag reduction rate decreases with the speed increase.This can be explained by the fact that with the increase of the ship's speed, the location where the bubbles reach the free surface moves backwards due to the drag effect of the fluid.This implies the area covered by the gas-water mixture moves backwards at a higher speed, resulting in a larger hull surface in the front in contact with ice.The drag reduction rate is thus reduced.This is however a tentative explanation of this interesting phenomenon.The effect of ship speed on the drag reduction rate of an air-bubble system requires systematic investigation, in a combination of other factors such as the injected air volume, which is included in the authors' future work.
In addition to the drag force in the longitudinal direction, the hull-ice interaction forces in the transverse direction regarding the air-bubble system are also analyzed.Figure 12 shows the time series of ice-induced drift force with and without the air-bubble system.It is seen from the figure that the ice-induced drift forces increase gradually in the first stage of the ice-going voyage for both cases.This is similar to the drag force, which can be explained by the fact that the air-bubble system has not been utilized when the ship enters the ice field.After the entrance stage up to t = 25 s, the ice-induced drift force is found to be smaller in both magnitude and variation.The mean value and the standard deviation of the ice-induced drift force without the air-bubble system are 14.0 kN and 94.3 kN, respectively.For comparison, when the air-bubble system has been activated, the mean and the standard deviation of the ice-induced drift force become 8.4 kN and 64.3 kN, respectively, which indicates the ice-induced drift force has also been reduced significantly.It is also noticeable that the ice-induced drift force has a large standard deviation.This is because the ice forces on the ship's sides are asymmetric.
As mentioned previously, when the air-bubble system is turned on with the air injection rate of 2.5 m/s, the hull-ice contact on the ship side is greatly reduced.Figure 10 illustrates the time series of the ice resistances for this case.In comparison with the ice resistances in Figure 9, the total resistance as well as the components from the bow area and the sides are found to be smaller.The resistance values are listed in Table 3, together with the resistance when the air-bubble system is off.The resistance reductions are also included in Table 3.It is seen that when the air-bubble system is off, the total ice resistance is reduced by 15.3%.When it comes to the ice resistance components, the bow resistance remains almost unchanged, with a reduction rate of 10.3%.In contrast, the ice resistance from the ship's sides is greatly reduced, with a drag reduction rate of 70.8%.The effect of the ship speed on the drag reduction rate was also investigated.In addition to the abovementioned ship speed of 6 knots, ship-ice interactions under four other speeds were simulated and ice resistances with and without the air-bubble system were compared.Table 4 listed the resistances as well as the drag reduction rates, which are also plotted in Figure 11.It is observed that the drag reduction rate decreases with the speed increase.This can be explained by the fact that with the increase of the ship's speed, the location where the bubbles reach the free surface moves backwards due to the drag effect of the fluid.This implies the area covered by the gas-water mixture moves backwards at a higher speed, resulting in a larger hull surface in the front in contact with ice.The drag reduction rate is thus reduced.This is however a tentative explanation of this interesting phenomenon.The effect of ship speed on the drag reduction rate of an air-bubble system requires systematic investigation, in a combination of other factors such as the injected air volume, which is included in the authors' future work.In addition to the drag force in the longitudinal direction, the hull-ice interaction forces in the transverse direction regarding the air-bubble system are also analyzed.Figure 12 shows the time series of ice-induced drift force with and without the air-bubble system.It is seen from the figure that the ice-induced drift forces increase gradually in the first stage of the ice-going voyage for both cases.This is similar to the drag force, which

Effects of Ventilation Rate
In this sub-section, a sensitivity study about the ventilation rate at the nozzle regarding drag reduction rate was carried out.Different air velocities from the nozzles were investigated under the condition of an ice concentration of 60% and a ship speed of 6 kn.The nozzles are divided into two groups: the bow nozzles are the first 3 pairs of nozzles at the bow area, and the side nozzles are the rest of 5 pairs at the sides that are near the bilge; see Figure 1 for the exact locations of the nozzles.The ventilation rate in m/s represents the gas flow rate at the nozzle.Four air velocities, i.e., 1, 2.5, 5 and 10 m/s, were investigated.A reference case with zero air velocity represents the condition when the airbubble system is turned off.The drag reduction for a total of eight cases with different air velocities was calculated.The ice resistance data and drag reduction rate under for these cases are listed in Table 5.In this table, Case D represents the air-bubble system with a ventilation rate of 2.5 m/s, which has been mentioned in the previous sub-sections.It is observed from Table 5, that with the increased ventilation rate the ice resistance is reduced.Let us look closer at the specific cases.Case A is when the bow-nozzles are turned off and the side-nozzle ventilation rate is set as 1 m/s.Figure 13 shows the simulation of the ice resistance time series of this case.When the side nozzles are turned on, air bubbles go up along the hull wall to the free surface.The ice floes at the ship's side become

Effects of Ventilation Rate
In this sub-section, a sensitivity study about the ventilation rate at the nozzle regarding drag reduction rate was carried out.Different air velocities from the nozzles were investigated under the condition of an ice concentration of 60% and a ship speed of 6 kn.The nozzles are divided into two groups: the bow nozzles are the first 3 pairs of nozzles at the bow area, and the side nozzles are the rest of 5 pairs at the sides that are near the bilge; see Figure 1 for the exact locations of the nozzles.The ventilation rate in m/s represents the gas flow rate at the nozzle.Four air velocities, i.e., 1, 2.5, 5 and 10 m/s, were investigated.A reference case with zero air velocity represents the condition when the air-bubble system is turned off.The drag reduction for a total of eight cases with different air velocities was calculated.The ice resistance data and drag reduction rate under for these cases are listed in Table 5.In this table, Case D represents the air-bubble system with a ventilation rate of 2.5 m/s, which has been mentioned in the previous sub-sections.It is observed from Table 5, that with the increased ventilation rate the ice resistance is reduced.Let us look closer at the specific cases.Case A is when the bow-nozzles are turned off and the side-nozzle ventilation rate is set as 1 m/s.Figure 13 shows the simulation of the ice resistance time series of this case.When the side nozzles are turned on, air bubbles go up along the hull wall to the free surface.The ice floes at the ship's side become thus overturned and move away from the hull.As a result, the occurrence of hull-ice interaction on both sides of the ship side is reduced.A further step is to turn on the bow nozzles.This is Case B, for which the simulation and ice resistance time series are illustrated in Figure 14.
For this case, the air velocities are set to 1 m/s for all the nozzles.It is seen from Figure 14 that observed the ice-free zone moves forward to the vicinity of the ship's shoulder, and at the same time, the width of the ice-free zone increases.However, the ice accumulation at the bow area remains almost unchanged.Comparing Case B with Case A, the drag reduction effect is not significant.
thus overturned and move away from the hull.As a result, the occurrence of hull-ice interaction on both sides of the ship side is reduced.A further step is to turn on the bow nozzles.This is Case B, for which the simulation and ice resistance time series are illustrated in Figure 14.For this case, the air velocities are set to 1 m/s for all the nozzles.It is seen from Figure 14 that observed the ice-free zone moves forward to the vicinity of the ship's shoulder, and at the same time, the width of the ice-free zone increases.However, the ice accumulation at the bow area remains almost unchanged.Comparing Case B with Case A, the drag reduction effect is not significant.To better interpret the drag reduction rates in Table 5, the ice resistance values under the various air velocities are paired for comparison, as shown in Figure 15.The cases in Table 5 are put into two categories: one is featured by the side-nozzle ventilation rate kept as 1 m/s; the other is characterized by the nozzles having the same ventilation rate.Comparing the two groups, it is found that the ice resistance of the two groups is quite close under the same bow-nozzle ventilation rate.This implies that the drag reduction effect is more sensitive to the bow ventilation volume.If the bow ventilation volume is sufficiently large, the side ventilation volume has a marginal contribution to drag reduction.This is To better interpret the drag reduction rates in Table 5, the ice resistance values under various air velocities are paired for comparison, as shown in Figure 15.The cases in Table 5 are put into two categories: one is featured by the side-nozzle ventilation rate kept as 1 m/s; the other is characterized by the nozzles having the same ventilation rate.Comparing the two groups, it is found that the ice resistance of the two groups is quite close under the same bow-nozzle ventilation rate.This implies that the drag reduction effect is more sensitive to the bow ventilation volume.If the bow ventilation volume is sufficiently large, the side ventilation volume has a marginal contribution to drag reduction.This is in line with the fact that for this vessel, the ice resistance on the side of the ship accounts for less than 30% of the total ice resistance.thus overturned and move away from the hull.As a result, the occurrence of hull-ice interaction on both sides of the ship side is reduced.A further step is to turn on the bow nozzles.This is Case B, for which the simulation and ice resistance time series are illustrated in Figure 14.For this case, the air velocities are set to 1 m/s for all the nozzles.It is seen from Figure 14 that observed the ice-free zone moves forward to the vicinity of the ship's shoulder, and at the same time, the width of the ice-free zone increases.However, the ice accumulation at the bow area remains almost unchanged.Comparing Case B with Case A, the drag reduction effect is not significant.To better interpret the drag reduction rates in Table 5, the ice resistance values under the various air velocities are paired for comparison, as shown in Figure 15.The cases in Table 5 are put into two categories: one is featured by the side-nozzle ventilation rate kept as 1 m/s; the other is characterized by the nozzles having the same ventilation rate.Comparing the two groups, it is found that the ice resistance of the two groups is quite close in line with the fact that for this vessel, the ice resistance on the side of the ship accounts for less than 30% of the total ice resistance.The work presented in this article is one of the first investigations on numerical simulation of air-bubble systems regarding drag reduction in floe ice fields.The ice conditions were simplified to ice floes of identical size and shape, which are the delimitations of the current work.An icebreaker was employed as the case study vessel for the demonstration of the proposed procedure.Other ship types with different hull forms need to be modelled to verify the robustness of the proposed procedure.Ice model tests are also required for validation of the numerical results, which are included in the authors' ongoing work.Despite the limitations, the proposed procedure is expected to facilitate design of new generations of ice-going ships.

Conclusions
In this paper, the authors made use of a coupling CFD-DEM approach in combination with the VOF method to simulate resistance in floe ice fields, aiming to establish a numerical analysis procedure for ice-going ships installed with air-bubble systems.From the simulations and analyses, a more distorted wave making around the hull is observed after turning on the air-bubble system.Ice floes in contact with the hull side wall are pushed away from the hull by the gas-water mixture, resulting in an ice-free zone close to the side hull.The ventilation rate of the air-bubble system is also studied.It is found that the drag reduction rate increases with the increase of ventilation but decreases somewhat at higher speeds.Side ventilation only contributes to reducing the side friction resistance, and the side friction resistance can be eliminated under low ventilation.In general, the bow ventilation plays a deciding role in the overall drag reduction.The work presented in this article is one of the first investigations on numerical simulation of air-bubble systems regarding drag reduction in floe ice fields.The ice conditions were simplified to ice floes of identical size and shape, which are the delimitations of the current work.An icebreaker was employed as the case study vessel for the demonstration of the proposed procedure.Other ship types with different hull forms need to be modelled to verify the robustness of the proposed procedure.Ice model tests are also required for validation of the numerical results, which are included in the authors' ongoing work.Despite the limitations, the proposed procedure is expected to facilitate design of new generations of ice-going ships.

Conclusions
In this paper, the authors made use of a coupling CFD-DEM approach in combination with the VOF method to simulate resistance in floe ice fields, aiming to establish a numerical analysis procedure for ice-going ships installed with air-bubble systems.From the simulations and analyses, a more distorted wave making around the hull is observed after turning on the air-bubble system.Ice floes in contact with the hull side wall are pushed away from the hull by the gas-water mixture, resulting in an ice-free zone close to the side hull.The ventilation rate of the air-bubble system is also studied.It is found that the drag reduction rate increases with the increase of ventilation but decreases somewhat at higher speeds.Side ventilation only contributes to reducing the side friction resistance, and the side friction resistance can be eliminated under low ventilation.In general, the bow ventilation plays a deciding role in the overall drag reduction.

Figure 1 .
Figure 1.(a) unsuitable mesh size; (b) suitable mesh size.The VOF model describes the phase distribution and position of the interface through the field of the phase volume fraction αi, the volume fraction   =    of the phase i, where Vi is the volume of phase i in the grid cell, V is the volume of the grid cell, And the sum of the volume fractions of all phases within each grid cell is 1.When the grid contains only a single fluid, the material properties of the fluid are used in the calculation.If there are

22 ,
10, x FOR PEER REVIEW 8 of 1 contrast, is supposed to reduce ice resistance instead.It would be sufficient to use a smaller volume of air from the bow/sides to push the crushed ice away from the hull.

Figure 2 .
Figure 2. the three-dimensional hull model and the nozzle locations of the air-bubble system.

Figure 2 .
Figure 2. The three-dimensional hull model and the nozzle locations of the air-bubble system.

Figure 3 .
Figure 3. Schematic diagram of the computational domain; (a) Schematic diagram of the computational domain at the initial stage; (b) the layout of the ice field with an ice concentration of 60%.

Figure 3 .
Figure 3. Schematic diagram of the computational domain; (a) Schematic diagram of the computational domain at the initial stage; (b) the layout of the ice field with an ice concentration of 60%.

Figure 4 .
Figure 4.The meshing of the numerical model.(a) Ship hull mesh; (b) Front view of the flow domain; (c) Side view of the flow domain.

Figure 4 .
Figure 4.The meshing of the numerical model.(a) Ship hull mesh; (b) Front view of the flow domain; (c) Side view of the flow domain.

JFigure 6 .
Figure 6.Wave patterns and streamlines around the hull for: (a) when the air-bubble system is turned off; (b) when the air-bubble system is activated.

Figure 6 .
Figure 6.Wave patterns and streamlines around the hull for: (a) when the air-bubble system is turned off; (b) when the air-bubble system is activated.

Figure 7 .
Figure 7.The snapshot at 120 s of the simulation when the air-bubble system has not been activated: (a) isometric view; (b) bottom view.

Figure 7 .Figure 8 .
Figure 7.The snapshot at 120 s of the simulation when the air-bubble system has not been activated: (a) isometric view; (b) bottom view.J. Mar.Sci.Eng.2022, 10, x FOR PEER REVIEW 12 of 19

Figure 8 .
Figure 8.The snapshot at 120 s of the simulation when the air-bubble system has been activated: (a) isometric view; (b) bottom view.

Figure 9 .
Figure 9.Time series of ice resistance when the air-bubble system is off: (top) total ice resistance; (middle) ice resistance from the bow area; (bottom) ice resistance from the sides.

Figure 10 .
Figure 10.Time series of ice resistance when the air-bubble system is on: (top) total ice resistance; (middle) ice resistance from the bow area; (bottom) ice resistance from the sides.

Figure 11 .
Figure 11.The resistance drag reduction rate at different speeds.

Figure 12 .
Figure 12.Comparison of time history curves of ice contact force along ship width.

Figure 12 .
Figure 12.Comparison of time history curves of ice contact force along ship width.

Figure 13 .
Figure 13.simulation and ice resistance of Case A: the side-nozzle ventilation rate is 1 m/s; the bownozzles are turned off.

Figure 14 .
Figure 14.simulation and ice resistance of Case B: the side-nozzle ventilation rate is kept as 1 m/s; the bow nozzles are turned on as 1 m/s.

Figure 13 .
Figure 13.Simulation and ice resistance Case A: the side-nozzle ventilation rate is 1 m/s; the bow-nozzles are turned off.

Figure 13 .
Figure 13.simulation and ice resistance of Case A: the side-nozzle ventilation rate is 1 m/s; the bownozzles are turned off.

Figure 14 .
Figure 14.simulation and ice resistance of Case B: the side-nozzle ventilation rate is kept as 1 m/s; the bow nozzles are turned on as 1 m/s.

Figure
Figure Simulation and ice of Case B: the side-nozzle ventilation rate is kept as 1 m/s; the bow nozzles are turned on as 1 m/s.

Figure 15 .
Figure 15.Ice resistance comparison for the various air velocities.

Author
Contributions: Conceptualization, B.-Y.N. and Z.L.; methodology, B.-Y.N. and H.W.; software, Z.L.; validation, H.W.; writing-original draft preparation, H.W.; writing-review and editing, B.-Y.N., Y.X. and Z.L.; project administration, Y.X. and B.F.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by National Natural Science Foundation of China, grant numbers 52192690, 52192693, 51979051, 51979056, U20A20327; and by National Key Research and Development Program of China, grant number 2021YFC2803400.And the APC was funded by National Natural Science Foundation of China, grant number 52192690.

Figure 15 .
Figure 15.Ice resistance comparison for the various air velocities.

Table 1 .
The main particulars of the case study vessel.

Table 1 .
The main particulars of the case study vessel.

Table 2 .
The boundary conditions of the computational domain.

Table 3 .
Ice resistance of the stable stage.
Figure 9.Time series of ice resistance when the air-bubble system is off: (top) total ice resistance; (middle) ice resistance from the bow area; (bottom) ice resistance from the sides.

Table 3 .
Ice resistance of the stable stage.

Table 4 .
Ice resistance at different speeds.

Bubble System off) Ice Resistance (kN) (Air-Bubble System on)
Figure 10.Time series of ice resistance when the air-bubble system is on: (top) total ice resistance; (middle) ice resistance from the bow area; (bottom) ice resistance from the sides.

Table 4 .
Ice resistance at different speeds.

Bubble System off) Ice Resistance (kN) (Air-Bubble System on)
Figure 11.The ice resistance and drag reduction rate at different speeds.

Table 5 .
Ice resistance under different air velocities.

Table 5 .
Ice resistance under different air velocities.