Numerical investigations of micro bubble drag reduction effect for container ships

In the course of the most recent decades reduction of ship resistance and saving fuel consumption to accomplish higher speed with reduction of pollutants has been the significant subject for researchers. Micro bubble drag reduction technique is one of the most interesting thoughts in this field owing to its great advantages, such as considerable potential drag reduction, easy operations, environmental friendliness and low costs. In this study a 3-D numerical investigation into frictional drag reduction by air micro bubbles is applied on KRISO container ship model. The objective is to understand the mechanism of resistance reduction through micro bubbles injection under model ship at different Froude numbers, injection rate and of course volume fractions. The numerical simulations are performed using a commercial CFD code solving Reynolds averaged Navier–Stokes (RANS) equations. A large number of simulations has been performed to investigate the effect of injection of micro bubble under ship model hull to estimate the local coefficient of friction values along ship hull model. The results show that at all of the examined Froude’s numbers, frictional resistance reduction attained at different rates and a maximum drag reduction of 27.6% was obtained at 0.282 Froude number with 4.8% air volume fraction.


Introduction
Ships are considered one of the most important vehicles worldwide as they transport 90 % of the world cargo and they are commonly one of the most economical method contrasted with other methods of transportation [1].
Despite the fact that they are very productive as far as energy utilization per ton of (cargo x mile) and therefore low emissions, their worldwide fuel consumption and environmental pollution are important parameters because of their broad use [2].
Marine engineers and maritime companies are striving to design and build ships with the lowest fuel consumption, in order to achieve the maximization of benefits and reduction of pollutants.
Moreover, rising fuel prices and restricted rules in emissions intensify the need for lower ship resistance and requisite propulsion power [3].
As ship travels in the water, total resistance on the submerged hull includes wave, pressure and skin friction resistance components.
Wave and pressure resistances are unavoidable and can be handled by the detailed design of the hull. However, the reduction in the skin friction resistance remains proportional to the viscosity of water, wetted surface area and the ship speed [4].
Frictional resistance causes 70-90% of total resistance for low-speed ships such as tankers and bulk carriers and less than 40% for high-speed ships like container carriers [5].
Active and passive methods are applied to achieve drag reduction in marine vehicles. The active method is accomplished through applying anti fouling materials and coatings, air lubrication techniques, and the use of riblets [6].
Passive methods include improving the shape of the vessel by application of modern hull forms or hull form optimization techniques [7].
Air lubrication presents one of the most efficient methods for the reduction of hull friction resistance. The economic and environmental impact of successfully implemented air lubrication could be significant, as a ship's fuel consumption can be reduced by 5 % to 20 % [8].
Air lubrication is a relatively young technology initially developed in 1961 at the Krylov Shipbuilding Institute; ALDR began with artificially inflated air cavities applied to the surface of low-speed river vessels and barges [12].
One of the first drag reduction techniques was the application of electrolysis-induced micro bubbles, reported by McCormick 1973 [13]. Bogdevich [12] mapped the peak drag reduction, which they found to occur immediately downstream of the point of air injection, while skin-friction returned to the conventional levels further downstream. Madavan [10] similarly showed high local drag reductions of up to 80%, which were later verified by Sanders [14].
Enthusiasm for these technologies waned as the oil crisis of the 1970s came to an end.
Since Madavan [15] accomplished a noteworthy 80% drag reduction in a water tunnel test, numerous analysts dedicated themselves to the study of the micro bubble drag reduction technique [16].
The significant exploration headings of the micro bubble drag reduction technique include: (1) Modeling and validating the mechanism of the micro bubble drag reduction technique by numerical methods [17]. (2) Validating the drag reduction effect in a water tunnel by applying different injection techniques, such as electrolysis method [13], the porous method slot and discrete hole injection [14], or the winged air induction pipe method [18] which is fitted on MV Olivia Maersk a full-scale 3000 TEU Maersk Line container carriers. (3) Investigating the effect of micro bubble drag reduction technique to the model ship in a towing tank [19].
Recently, the air lubrication system was tried in full scale on a new build ship from Mitsubishi Heavy Industries [20], called Till-Deymann. (4) Analyzing the mechanisms of drag reduction experimentally and studying the various factors affecting of the micro bubble drag reduction technique, such as the bubble concentration effect [21], the bubble size effect [14], the effect of micro bubbles on the structure of turbulence [21], location of micro-bubble injection [21], distance of micro-bubbles away from the solid surface and void fraction [22].
Lately, some significant exploration results have been confirmed for the micro bubble drag reduction techniques.
First, total of 18 full-scale ship performance tests between 2002 and 2015 with net energy savings reported between 4% and 10% [23].
Second, Ceccio [24] observsed that air bubble diameter mainly depends on main flow velocity at the place where air is injected and the air flow rate, while the micro bubble diameter has no effect on the reduction of frictional resistance. (Not by the size of the hole) [2].
Murai 2014 [25] reported that the best size of bubbles is an order of magnitude larger than the viscous sublayer and an order of magnitude smaller than the boundary layer thickness, while, bubbles with the size of same order as boundary layer size can increase the drag [2].
However, Azlina [26] reported that till now, our knowledge of the influence of bubble size on micro-bubble drag reduction is incomplete.
Third, most of researchers fitted the micro bubbles injectors close to the bow of the model in their studies. [16] and [20] However, in 2013 [27] it was recommended that a design standard for air injection location outlets based on ship type to cover the hull surface by air as widely as possible.
So, it was that recommended for the blunt ship, which has a wide flat area such as bulk carrier, only one outlet is placed near the bow and, for slender ship such as ferries, the flat area of which is narrow on the bow and aft, three outlets are placed separately.
Researches in [28] have established a variety of experiment with the application of air lubrication method on a tanker model ship and they obtained a maximum drag reduction of 18.13% near the design speed 15 kts and they found out that it was more effective to distribute the air injection at multiple locations than single injection.
Fourth, researchers in [22,25] reported that the bubble drag reduction has its main impact just near the injection point but continuing downstream from the injection point, the bubbles leave the wall boundary layer and consequently their impact reduces significantly.
Fifth, the most important parameter for microbubble drag reduction is the volume fraction of the injected air in the water boundary layer [16].
Ceccio [24] found that the effective gas phase volumetric flow rate is the most important parameter in achieving drag reduction by gas injection.
Authors in [22] reported that the most important parameter in deciding the MBDR effect is the effective void fraction at the point of analysis, which is influenced by both the injection rate and the static pressure under the test conditions.
They mention that reduction in the frictional drag increased with increasing mean void fraction, but at higher injection rate, excessive micro bubble injection destroys the favorable turbulent boundary layer and in doing so decreases the drag reduction effect.
Sixth, Azlina [26] believed that many factors contributed to the positive impact on the drag reduction, which among them is density effect, turbulence suppression, and Reynold's stress and near-wall void fraction.
Authors in [17] carried out a numerical simulation in order to study the relationship between drag reduction effect and the combination of jet position for a low-speed ship with a long mid-plane section and the obtained results are used as the reference for the towing experiment and practical application.
In addition, they reach with friction drag reduction rate up to 22% and found that good distribution of gas film on the parallel middle body of the bottom is an immediate reason for drag reduction; hence they concluded that numerical results are in good agreement with those of the model experiment, and the relative error is between 2% and 6%.
Sudhir et al. [22] conducted a numerical investigation for a flat plate in a 3D turbulent flow channel, and they studied the effect of micro bubbles on frictional drag reduction with varying flow speed, void fraction and geometry at the injection point. They compared numerical results with experimental study conducted by Moriguchi Y, Kato H [29] and they found that numerical results are in good agreement with those of the model experiment.
Tsai [16] derived a simple boundary layer mixture model to predict the drag reduction effect of the micro bubble drag reduction technique for a flat plate. The flat plate resistance measuring system with a porous material micro bubble injecting system is set up to measure the total resistances of the plate without and with injecting micro bubbles. The tests are conducted in a water tunnel to verify the drag reduction effect predicted by the boundary layer mixture model. In addition, they reported that: the drag reduction effect predicted by the boundary mixture model is almost directly proportional to the density ratio of mixture and water and the frictional resistance coefficient is in good agreement with the value predicted by the boundary mixture model.
In this study, a 3-D numerical scheme is applied to a model scale of KCS with and without micro bubbles air injection using ANSYS FLUENT software for different Froude's numbers and air flow rates. In the numerical setup, air bubbles were injected through a slot at the upstream upper surface, generating flow of air bubble mixture. A boundary mixture model stated by Jing-Fa Tsai [16] is used to estimate air injection flow rate. More than 36 simulations performed to investigate the effect of injection of micro bubble under ship's hull for different Froude numbers and void fraction. The numerical framework consists of Reynolds averaged Navier-Stokes (RANS) equation and the standard k-ɛ with enhanced wall function treatment, because of two-phase flow, mixture model is chosen for multi-phase flow.

Ship model
A physical model of KRISO container ship was numerically tested in Fig. 1.
The model is developed in the Korea Research Institute for Ships and Ocean Engineering (KRISO)-No full-scale ship ever exists-and it is used as a benchmark for flow computations around ships [30].
Recently Elsherbiny et. al [31] performed numerical and experimental study on KRISO model to examine a hydrodynamic performance of ships advancing through different canal cross sections.
In our study the model tests performed based on Froude similarity with a scale of 1/31.6. Table 1 gives the ship and model particulars [32].

Parameter definitions and investigation of model conditions without air micro bubbles injection
The flow around model hull is modeled using RANS equations employing the Finite Volume Method (FVM). Water is assumed incompressible in calculations, thus the volume of water entering computational cells in vicinity of model hull will match with an equal volume of water flowing out, leads to continuity equation.
Together with Navier-Stokes equations, conservation of momentum of the flow can be defined and flow field around model hull can be characterized.
Detailed description of the mathematical model and numerical methods are presented in [33].

Mathematical model
The governing equations applied in this study are conservation of mass (continuity equations) and momentum (RANS) for incompressible turbulence flow. These equations are expressed in cartesian tensor form as follows: where u i is the time averaged velocity components in Cartesian Coordinates x i (i = 1, 2, 3), p, ρ and μ are the static pressure, fluid density and fluid viscosity respectively, ij is the Kronecker delta and −ù iùj is the Reynolds stress, where To close the RANS equations, the Standard k-ε was used, where: k , , u i , , , t are turbulent kinetic energy, dissipation rate, time-averaged velocity, density, dynamic viscous coefficient and turbulent viscosity of the mixed fluid respectively. t is time and k , are Prandtl numbers corresponding to the mixed fluid turbulent kinetic energy and dissipation rate. P k is the generation term of the turbulent kinetic energy due to the average velocity gradient. c 1 , c 1 are constants where the in this study have taken the following value:

Computational domain and boundary conditions
The numerical setup used for the present study is shown in Fig. 2 and it is according to ITTC Procedures [34].
Only the wetted surface area at the design draft will be studied, the control volume was selected to be of box shape.
Calm water conditions and air velocity are considered the same, however the aerodynamic forces are neglected.
The model fixed in its floating position and it is not free to sink and trim. Six Froude's numbers are studied: 0.108, 0.152, 0.195, 0.227, 0.260 and 0.282.
The computational domain height is 1.15 LBP and its width is 1.5 LBP because of the symmetry of the model.
In order to save computational resources in the simulation time the geometry is split into halves and symmetry boundary conditions are applied for the physical geometry of interest.
The inlet boundary located at 2.0 LBP ahead of the model, while the outlet boundary is at distance of 5.0 LBP from the model stern.
As shown in Fig. 2, uniform flow velocity condition is applied at the inlet of the boundary. A hydrostatic pressure outlet is applied for the boundary outlet.
The mid-plane is modelled as a symmetry plane and the model hull surface specified as a no slip condition wall. The top surface modeled as free slip wall and considered as static water surface with no free surface effect in addition, the bottom and far field surface modeled as free slip walls.
The objective for this case is to realize a verification and validation study. The validation is done through comparison to both experimental Tokyo 2015 workshop [32] and ITTC 57 results.

Grid generation
The computational domains are discretized with tetrahedral unstructured grid system, ANSYS mesh software, The preprocessor of FLUENT is used to generate the three-dimensional grid around the hull under present study.
Grid generation of the computational fluid domain followed the Mesh-Independence solution method as the density of meshing affects the flow field and the results. Therefore, the mesh is intensified to a limited extent where the skin friction coefficient results varied not any more. In our study grid convergence study is performed based the ITTC uncertainty analysis recommendation [34].
At least three meshes are required to estimate convergence performance which were categorized into coarse, medium and fine mesh.
These should be refined with same refinement ratio ( r G ) in all coordinate directions I, J, K.
Based on ITTC, it's recommended to implement refinement ratio with (r G = √ 2) . The mesh was varied by modification of the control volume and model hull element size.
In order to resolve the boundary layer accurately, 20 inflation layers are placed along the hull surface. The stretching factor (Gap factor) of inflation layers is constant throughout the analysis. Table 2 shows the details of meshes used for the convergence study.
Typical computational fine mesh used for simulation is shown in Figs. 3, 4.

Computational setup and numerical simulation
The SIMPLE algorithm is used to couple the momentum and continuity equations. Iterative technique is used for solving k-ɛ turbulence quantities as computational boundary conditions contains certain approximations.
The solution is converged when the normalized residuals of mass, momentum and turbulence values are less than 10 −3 .
The simulation conditions are carried out for viscous flow field around KCS model at six different Froude's numbers ranging from 0.l088 to 0.28287.
The simulations carried out in a single node of an Intel Core I7 CPU with 4 cores, clock speed 2.4 GHz and 6 GB of RAM memory. Solution setting for the cases are summarized in Table 3.

Grid independence study
In order to perform grid independence study a simulation has been processed to calculate total resistance at 0.108 Froude's number.    Table 4 shows the grid convergence study for the total resistance coefficient at 0.108 Froude's number.
It is noted that as number of grid elements increase, the total resistance results agree fairly well with result obtained from experiment.
Based on ITTC, grid convergence is calculated using ratio of changes in solution where R G , 21 and 32 is convergence ratio, changes of solution between the medium-fine and coarse-medium grids respectively.
From Table 4, we estimate R G = 0.00032 and according to the following cases which can occur: The solution can be considered as converged as estimated R G meets the convergence condition and the fine mesh setting will be applied for all different Froude's numbers.

Resistance results
The simulation data has been processed to predict drag coefficient of the ship model and a total of six cases were simulated for different Froude's and Reynolds numbers.
The total resistance acting on the ship is divided into two parts, frictional R F and residuary resistance R R where it is equal to sum of viscous pressure resistance R pf and wave making resistance R w .
The resistance component resulted from CFD simulation is divided to frictional resistance and viscous pressure resistance component neglecting wave making resistance.
The total resistance expressed in a coefficient form is given by: where is the water density, S is the wetted surface area of the model and V is the model velocity.
The frictional resistance acting on model can be calculated by performing an integral of the wall shear stresses over the wetted surface.
Where the frictional resistance coefficient is given by: To validate the resistance results for the present solver, the results were compared with the experimental data presented in the Tokyo 2015 workshop [30] and with ITTC 57 formula as:

Relative deviations are calculated according the following equation is
where C f CFD is the frictional resistance coefficient obtained numerically and C f EFD is the frictional resistance obtained experimentally. Table 5 shows the relative deviation RD for all test results comparing with ITTC 57 and EFD results.
It is observed that mean of relative deviation in case of frictional resistance prediction is found to be around 4.6 %. However, the mean of relative deviation for total resistance is around − 11.67 %, where the relative deviation increases as Froude's number gets higher. This is mainly because of the assumption that the water surface as static surface with no wave making resistance.
The prediction of total resistance coefficient can be improved by applying free surface effect in top surface boundary condition. However, in our study we will focus only on the frictional resistance for all cases. Figures. 5, 6, show the comparison of the total resistance coefficients and fractional resistance coefficient obtained numerically in the present study with the experimental data obtained in Tokyo 2015 workshop [32] and with frictional resistance coefficients computed using ITTC 57.

Parameter definitions and investigation of model conditions with air micro bubble injection
In our case, air micro bubbles are injected through a slot at the upstream upper surface, generating flow of air bubble mixture.
The injection slot is 0.25 m long, 0.02 m wide and 2.5 m from the bow. Fig. 7 shows the micro bubble air injection location.

Mathematical model
When the air micro bubbles are injected, the flow around the model becomes two phase. For simulating two phase flow, the 'mixture model' is implemented in FLUENT. This can be used to model multiphase flow where the phases move at different velocities.
The mixture model can model 'n' phases (fluid or particulate) by solving the momentum and continuity for the mixture, the volume fraction equations for the secondary phases, and the algebraic expressions for the relative velocities.An extensive discussion on this model may be found in Manninen et al [35].
The governing equations: The continuity equation, the momentum equation and the equation of the gas-liquid two phase mixed flow are as follows: The continuity equation for the mixture is are the density of mixed fluid and the mass averaged velocity respectively. w and b are the volume fraction of fluid and gas respectively. b + w = 1 subscripts b, w represents bubbles and water. The momentum equation of the mixture is obtained by summing the individual momentum equations for all phases. It can be expressed as  where P is the pressure and m = b m + w w is the viscosity of the mixed fluid.
⇀ v m are the drift velocity of water and gas respectively.
The volume Fraction Equation for the Secondary Phases: From the continuity equation for secondary phase b, the volume fraction equation for secondary phase can be obtained: Slip between water and gas is reflected in the interaction between the two phases. The relative velocities between gas and water is defined as the slip velocity: The relationship between gas drift velocity and slip velocity is

Numerical setting
The model is simulated with the fine mesh produced in section 3.3, boundary condition is the same as previous case with the difference that for the air-injection section, the air mass flow rate boundary used where the micro bubble injection is perpendicular to slot surface and the air volume fraction α is estimated according to (13) 3. Solution settings for the case are summarized in Table 6.

Estimation of air injection rate Q a
There is no immediate set up technique to assess micro bubble air injection flow rate.
However, this parameter can be estimated as function of boundary water flow, which can be determined just by considering some common turbulent fluid relations [36].
The air volume fraction C v is defined as the ratio of the injected air flow rate divided by the summation of the air flow rate and the water flow rate within the boundary layer.
where Q a the air injection flow rate, and Q w is the water flow rate within the.boundary layer of model.
In this study the micro bubble injection flow rate is measured in terms of the equivalent air film thickness t l ranging from 2-9 mm as follows: where Q a , B and V m are the micro bubble air injection flow rate, the wide of air injection slot and model ship speed, respectively.
A survey on [28] regarding patents [37,38] for the air lubrication tells that the suitable value for t l is recommended to be from 4 to 8 mm.

Estimation of water flow rate Q w
The water flow rate within the boundary layer of the KRISO ship model can be predicted based on the turbulent boundary theory [16] using the following equation: where b is the width of the model ship, U 0 is the model velocity, L wl is the length of the water line for the KRISO model and R el is Reynolds number. Figure 8 shows water injection flow rate versus Froude's number.

Estimation of air volume fraction C v
The air volume fraction C v in Eq. (16) can subsequently be calculated with. measured injected air flow rate and water flow rate which was estimated in section 4.3.2.
As shown in Table 7 simulations were performed at six different Froude's numbers ranging from 0.108 to 0.282 with different air mass flow rate (air volume fraction) More than 36 simulations are performed to investigate the effect of injection of micro bubble under ship's hull where at each simulation local coefficient of friction values measured along ship hull model for 11 longitudinal positions, and at each longitudinal position, five in number at transverse positions.
As shown in Fig. 9 the location points for estimation the local coefficient of fraction in longitudinal and transverse direction on one side with 100, 200, 300, 400 mm a way from center line of the model. Table 8 lists the values of skin friction resistance and the drag reduction percentage for all cases considered, and Fig. 10 displays the relationship between drag reduction percentage and air volume fraction for different Froude's numbers.  It can be seen that at all examined Froude's numbers, drag reduction is attained. The maximum drag reduction of 27.6 % was found at 0.282 Froude's number with 4.8% air void fraction.

Results and discussion
Moreover, it is observed for 0.108 Froude's number with lower void fraction, the total drag reduction was 8.0 % as air bubbles cover the area below the hull fully against 5.4% of drag reduction at the higher void fraction where the air bubbles escape from the sides of the hull.
Further investigation is carried out for lower injection flow rate with 0.0046 air volume fraction, the percentage of drag reduction was 9.5% for this case. Figure 11 shows the comparison of local air/water mixture percentage distributed below model hull for two different volume fraction 0.0046 and 0.04 at 0.108 Froude's number. This indicates that micro bubble drag reduction effect decreases with increasing airflow rate.
It is noticeable here that for lower ship speed, the drag reduction increases with decreasing injection flow rate.  At higher injection, the bubbles escape from the sides of the hull and increases the turbulence of the flow, which change the flow parameters and increases the local frictional drag coefficient these results agree with Sudhir et al. [23] finding. Here, w refers to local shear stress with the injection of micro air bubbles and w0 refers to local shear stress without the injection of bubbles.
From these Figures, it's observed that with injection of air bubbles the local shear stress decreases gradually in the downstream direction.
It was also concluded that with increasing both air volume fraction and Froude's number the total drag reduction increased from 8% in the first case to 27.0% in the second case.
It is observed in Fig. 14 that for most of the transverse locations increasing distance from the injection point in the longitudinal direction, values of local C F initially decreased due to the formation of air bubbles mixture and with the further increase downstream, the values increased.
This reconfirms that, micro bubble effect reduces with the distance from the injection point both in longitudinal and in the transverse directions.
At 0.195 Froude's number, Comparison of local air/water mixture percentage distributed below the hull at two different air volume fractions 0.01 and 0.045 is shown in Figs. 15, 16.    As seen with lower injection of volume fraction of 0.01, injected air bubbles do not cover the area below the hull fully.
However, at the higher volume fraction, a continuous layer of air has been formed displacing the water in contact with hull.
With the void fraction of 0.01 the reduction in the total drag of 14.5 % was observed as against 20.7 % of reduction at the injection void fraction of 0.045.

Conclusions
Based on the results obtained from the present simulations the following conclusion can be made: (1) At all examined Froude number, resistance reduction occurred at different rates with a maximum drag reduction of 27.6% obtained at 0.282 Froude's number with 4.8% air void fraction. (2) At higher Froude's number, the air injection flow rate must be increased to get higher drag reduction effect because at higher Froude's number bubbles are pushed out of the boundary layer and do not cover the hull bottom fully. (3) At lower Froude's number, the air injection flow rate must be decreased to get higher drag reduction because the flow stream can move air bubbles to whole region of the boundary layer and cover the hull bottom fully. (4) Micro bubble drag reduction effect decreases with increasing airflow rate. At higher injection flow rate, the air bubbles escape from side of the hull and increases the turbulence of the flow, which change the flow parameters and increase the local frictional drag coefficient. (5) The distance from injection point is a paramount factor in achieving the reduction in frictional drag as micro bubble drag reduction effect reduces slowly in streamwise direction and quickly reduces in span-wise direction.

Recommendation for future research
To the best knowledge of the authors, this study is the first numerical investigation on the effect of air micro bubbles injection for KRISO ship model and future experimental investigations are necessary to validate numerical results.