Investigation of different flow parameters on air layer drag reduction (ALDR) performance using a hybrid stability analysis and numerical solution of the two-phase flow equations

Ocean Engineering xxx (xxxx) xxx

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Investigation of different flow parameters on air layer drag reduction (ALDR) performance using a hybrid stability analysis and numerical solution of the two-phase flow equations Mohammad Hossein Montazeri *, Mohammad Mehdi Alishahi Mechanical Engineering Department, Shiraz University, Shiraz, Iran

A R T I C L E I N F O

A B S T R A C T

Keywords: Drag reduction Computational fluid dynamics Stability Hydrophobic

Among different methods of drag reduction, Air Layer Drag Reduction (ALDR) has been shown to be promising, at least in laboratories. Although a few examples of its practical application is available, important questions regarding the physics of ALDR and consequently its design rules and implementation remain to be answered in order to enhance its reliable, efficient and widespread use. Montazeri and Alishahi (2019) introduced a hybrid numerical procedure for such a problem, capable of producing acceptable and feasible results in comparison with experiments. In the present study, this procedure is modified and reapplied to the problem and its results are studied in more depth and some aspects of ALDR that had been overlooked or even misinterpreted in experi­ mental tests, are discussed. The results show that the air layer develops and then disintegrates downstream in most cases. The ruptured air layer results in mild drag reduction for small air fluxes but high drag reduction happens for large flux values. Also, it is shown that the ALDR conduct is insensitive to bumps usually procured during construction. Results show that employing a superhydrophobic coating on a small surface prior to the injection point provides improved efficiency that would seem promising in practice.

1. Introduction Ships are one of the main transport vessels and large share of com­ mercial transportation belongs to them. The main advantage of marine transportation especially via ships over other transportation means is their lower price tag per mile per ton of freight loading. Large cargo ships bring about efficiency and cost effectiveness in spite of their huge fuel consumption. Large number of containers in a cargo ship can make the overall cost per loading lower than other transportation means. However, any reduction of fuel consumption in ships is still interesting as it may result in even cheaper transportation or faster ships. Reduction of fuel consumption even as small as 1 percent would bring about large amount of oil savings which is not only interesting financially, but also environmentally. Fuel savings can be categorized into different methods such as: increasing the engine or thruster efficiency, optimizing the geometry of the hull to reduce its form drag and reducing the hull skin friction drag. Although any of these choices can be more important than the others per case, the share of skin friction drag is usually more than form drag (pressure drag) in most cases of conventional commercial ships (Butterworth et al. (2015); Park et al. (2015)). This is more

pronounced for larger ships e.g. oil tankers and as a result, many re­ searches have been focused on reduction of the skin friction drag of these ship’s hull. Researchers have investigated different methods of friction reduc­ tion considering especially their performance and applicability. Using super-hydrophobic surfaces, injection of polymers, micro-bubbles or air layer, using partial cavity and traveling inside a supercavity bubble are the most important investigations. However, some of them (e.g. super­ cavity) cannot be applied to large ships at relatively low speeds. Although most of these drag reduction methods result in acceptable drag reduction, their practical application on commercial ships encompass difficulties. Implementation of super-hydrophobic coatings on wall surface decreases the frictional drag (Henoch et al. (2006); Daniello et al. (2009)) accompanied with an initial and maintenance cost as well. However, this method entails some disadvantages; the most important one is that the drag reduction by this method depends completely on micro/nano-bubbles of air entrapped inside the micro/nano surface roughness structures. Based on experimental data, Wang et al. (2014) stated that these bubbles may escape from their place when water flows over them for enough time period. The super-hydrophobic surface

* Corresponding author. E-mail address: [email protected] (M.H. Montazeri). https://doi.org/10.1016/j.oceaneng.2019.106779 Received 18 June 2019; Received in revised form 24 November 2019; Accepted 25 November 2019 0029-8018/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: Mohammad https://doi.org/10.1016/j.oceaneng.2019.106779

Hossein

Montazeri,

Mohammad

Mehdi

Alishahi,

Ocean

Engineering,

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

without its micro/nano air bubbles would result in no drag reduction. They claim that this problem can be solved by continuous air injection from the wall. In this method the hydrophobic surface should also have porosity and the injected air replaces the escaped air bubbles. Although the efficiency of this method seems to be high because of low air leakage rate from hydrophobic surface, construction and maintenance of a porous hydrophobic hull’s surface adjacent to sea water would be difficult. Another drag reduction method investigated in the literature is bubble injection method. In this method, micro/mini bubbles are injected from the wall. The presence of air bubbles inside the boundary layer with lower density and viscosity relative to water, changes the velocity and turbulence kinetic energy distribution in wall boundary layer (McCormick and Bhattacharyya (1973); Kanai and Miyata (2001)). This change finally results in an acceptable drag reduction (Sanders et al., 2006). Air is usually injected from a porous wall but injection from thin slot is also investigated by researchers (Maltzev (1995); Elbing et al. (2008)). Murai (2014) listed many experimental and numerical in­ vestigations on drag reduction using bubble injection. In this method, the ratio of the drag reduction to the required air flux (which is called “gain factor” by Murai (2014)) can be very high because of low air flux requirement. However, there are some issues on the applicability of this method; many researchers, e.g. Butterworth et al. (2015), studied the size of injected bubbles and reported that an acceptable drag reduction can be achieved using proper bubble size. Madavan, Deutsch and Merkle (1984) and Elbing et al. (2008) stated that the size of bubbles changes when they flow downstream from the injection point due to many breakups and coalescences. This change in bubble size may decrease the drag reduction such that Verschoof et al. (2016) experimentally showed 40% drag reduction using 1 mm air bubbles while 0.1 mm bubbles at the same conditions cannot produce drag reduction more than 5%. Also many researchers including Kodama et al. (2000) have shown that the injected bubbles leave the wall adjacency when flowing downstream due to boundary layer velocity distribution and corresponding effects. In this manner, the effect of these bubbles in drag reduction would be almost vanished. Watanabe, Masuko and Shirose (1998) showed that in their experimental test on a 40 m flat plate, the effect of injected bubbles decreased from 40% drag reduction at the injection point to about 10% at 20 m downstream. Also, Elbing et al. (2008) showed that in their experimental study on 13 m flat plate, the drag reduction due to bubble injection decreases from about 80% at 0.5 m downstream of the injec­ tion point to less than 20% at 4.5 m downstream. This rapid decrease in drag reduction shows the inefficiency of bubble injection method for large hulls. Although, multiple injection points may be used to overcome this problem, the required air flux would be higher and construction and maintenance of multiple porous injectors would be more involved. One of the interesting methods of drag reduction that has been investigated more recently, is making an air layer at the hulls surface. Required air flux of this method is higher than the micro/mini bubble injection method but the drag reduction percentage is also remarkable. DNS result reported by Kim and Moin (2010) showed about 100% skin friction reduction and experimental tests results reported by Elbing et al. (2008), Elbing et al. (2013) and Jang et al. (2014) conveyed more than 80% reduction in friction drag. This level of drag reduction is reported only for air layer and partial cavity methods. In partial cavity method, a cavity should be constructed in the hull with a backstep at upstream and a recess at the end part of the hull. This cavity should be filled with air and a continuous air injection inside the cavity would compensate the leakage of air from the recess. As the leakage flux is low, the required air flux of this method is lower than air layer method. Makiharju et al. (2013) reported the air flux requirement for the partial cavity method is about one third of the air layer method. Alongside the low required air flux and high level of drag reduction, partial cavity method entails some weaknesses. The most important weakness of this method is that the hull’s cavity with a specific depth and length should be designed only for a specific speed of ship.

Makiharju et al. (2013) results showed that the required air flux for partial cavity grows very rapidly when flow speed is lower or higher than the designed point. Also as Rotte et al. (2016) mentioned, if the air injection system is turned off, total drag of hull with cavity would be more than a simple flat hull geometry. These two problems of partial cavity method make it practically undesirable. Regarding different aspects of the above methods, the conclusion would be that the Air Layer Drag Reduction (ALDR) is capable of pro­ ducing high drag reduction ratio in a vast range of freestream velocities. The ALDR is mostly studied in experimental tests by researchers, e.g. through investigations by Elbing et al. (2008) and Elbing et al. (2013). In these studies, different air fluxes are injected onto a 2D flat plate wall surface from small values (about 0.01 m2/s) to large values (about 0.15 m2/s). Elbing et al. (2008) used wall skin friction sensors at 6 different positions downstream of the injection point and their drag reduction were reported for different air fluxes and different water speeds ranging from 6.7 to 18 m/s. Elbing et al. (2013) added LDV, pitot tube and Time-of-Flight (ToF) sensors to find the velocity and void fraction dis­ tributions. Elbing et al. (2008) and Elbing et al. (2013) showed that minor drag reduction (less than 30%) can be achieved with relatively low air flux. Increasing the air flux raise the drag reduction percentage almost linearly to values more than 80% which they choose as the air layer development threshold. They interpreted the results of the skin friction sensors such that there is a bubble drag reduction regime for low air fluxes and then for higher air flux, a transitional air layer would occur where the bubbles are large and air patches are formed in the wall ad­ jacency but the water stream remains in contact with some parts of the wall surface. At last, for a specific air flux (which depends on flow conditions such as water free-stream velocity, surface roughness and upstream perturbations) the air layer would be completely developed and the drag reduction would be more than 80% at the whole wall surface. They call this air flux as the “critical air flux”. Elbing et al. (2013) repeated the experimental test of Elbing et al. (2008) for free-stream velocity of 6.3 m/s. In this new test the air was injected just downstream of a backward step. Although the existence of the backstep is claimed to reduce the sensitivity of the air layer to up­ stream perturbations by Elbing et al. (2008), Jang et al. (2014) showed that implementation of side guards are necessary to prevent side air leakages when using injection behind a backstep. Elbing et al. (2013) investigated the effect of upstream perturbations by inserting vortex generators at different angles of attack and showed that the drag reduction does not depend on vortex generators angle of attack signifi­ cantly when air is injected behind a backstep. There are also some numerical simulations of the ALDR flow. Wang, Yang and Stern (2010) used a URANS turbulence model (Spallart-Al­ maras model) and tried to simulate the Elbing et al. (2008)’s experi­ mental test, numerically. They found that for very low air fluxes (about 0.007 m2/s) there is only a minor drag reduction ratio (lower than 10%) but at 0.007 m2/s air flux, a jump in drag reduction is predicted and for air fluxes higher than this value the drag reduction ratio is about 100%. This jump in drag reduction is not the case in the experimental results of neither Elbing et al. (2008) nor other researchers. Also the critical air flux predicted by Wang et al. (2010) is about 0.007 m2/s in 6.7 m/s water free-stream velocity where the experimental results show a value of about 0.04 m2/s. In another numerical simulation, Kim and Moin (2010) chose a small numerical domain (with the length of about 40 cm) and two-phase direct numerical simulation of air injection behind a backward facing step was carried out. Injection of air behind a back-step was previously investi­ gated by Elbing et al. (2008) and it was expressed that the insertion of a back-step just upstream of injection point would reduce the air layer sensitivity to upstream perturbations. This problem is also studied by Elbing et al. (2013) and Jang et al. (2014) afterwards. The DNS results of Kim and Moin (2010) were qualitatively compared to experimental re­ sults as there were no similar experimental case with the simulated geometry and dimensions. Kim and Moin (2010) showed that a low air 2

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

flux injection cannot result in a stable and continuous air layer but at high air flux, a continuous air layer on the wall would totally prevail. Although this is qualitatively in agreement with Elbing et al. (2008) results, an important difference can be seen in the DNS results (Fig. 6 of Kim and Moin (2010) paper) in comparison with Elbing et al. (2008) claims. The DNS results show that when the air flux is not sufficiently high, an air layer would be formed near the injection point, but there are some growing perturbations that finally cause disintegration of the air layer. Both the free surface shape and wall skin friction plot show that the air layer exists at near injection region for low air flux. This is not in agreement with Elbing et al. (2008) and Elbing et al. (2013) statements. They stated that for air fluxes lower than the critical value, the air layer could not be developed and only bubbles and air patches would be present. In another numerical study of ALDR, Montazeri and Alishahi (2019) investigated different URANS turbulence models and showed that the most physical results among the investigated models would be produced by the well-known SST model. They showed that this model predicts the overall attribute but the predicted flow is more stable and continuous air layer is present even at lower air fluxes than what was reported in experimental test of Elbing et al. (2008) (this was previously reported by Wang et al. (2010) paper using Spallart-Almaras model). They showed that the simulation weakness is due to incorrect stability prediction in URANS simulations. With the addition of proper perturbation to the inlet boundary condition, comparable results with experimental tests involving minor errors would be obtained. They performed many CFD simulations with different perturbation frequencies and found that there is a frequency which leads to the most unstable flow; this most unstable flow results in drag reduction similar to the experimental results re­ ported by Elbing et al. (2008). They also introduced a solution procedure which uses linear stability analysis to find the proper perturbation fre­ quency and showed that the most unstable frequency extracted from this linear stability analysis is in fair agreement with the most unstable case predicted by CFD simulations with different perturbation frequencies. Two CFD simulations with proper URANS model plus one linear stability analysis comprising in their method require much less com­ puter resources and time if compared with DNS solutions and hence, it is suitable for practical applications. It would be a valuable tool for studying the physical aspects of the flow field and specifically, deter­ mining the outcome of various parameters adjustment involved in this problem. Montazeri and Alishahi (2019) validated this procedure mainly by comparing the numerically predicted drag reduction values with those of experimental results. In this paper the same numerical procedure is used to investigate the physical phenomena of ALDR flow in more depth. Also, the sensitivity of the air layer to surface bumps is studied to find the robustness and stability of air layer encountering construction faults and appendices. Then the effects of some parameters is studied to find out if the efficiency of ALDR can be enhanced in practical situations.

boundary conditions should be performed. The most unstable fre­ quency from different modes of the flow should be extracted (usually the internal mode is dominant when both internal and interfacial modes are present). 4 The CFD simulation in the 1st step should be repeated with the addition of a harmonic perturbation with frequency of the unstable mode obtained in the 3rd step. As the first step, a simple 2-dimensional CFD simulation is carried out by commercial software, ANSYS Fluent, at two free-stream velocities of 6.7 and 11.1 m/s, each with different air fluxes. The geometry and flow conditions are the same as Elbing et al. (2008)’s tests. Water flows on a 13 m long flat plate and air is injected at 1.38 m location downstream of the flat plate leading edge. This step is implemented exactly in the same way as Montazeri and Alishahi (2019)’s investigation and hence, full descriptions of the grid independence study, turbulence model selection and other numerical aspects are omitted for the sake of briefness in this paper. As they mentioned, simple simulations with common turbulence models does not produce appropriate results and among different investigated URANS models, SST model is the only one that correctly predicts the overall behavior of the flow and the relevant physical phenomena as revealed in experimental tests. However, simulations using SST model result in an air layer with more stability and longer continuous air layer than experimental results. Fig. 1 shows the air layer thickness for these simulations. There is no additional perturbations in these simulations other than the inevitable numerical (discretization) error. In this figure, nominal air layer thick­ ness (ta) is used instead of air flux to make two figures with different free-stream velocities comparable. Nominal air layer thickness (ta) is a simple parameter that represents the ratio of 2D air flux to free-stream velocity, (Qa/U). This parameter has been used by other researchers including Elbing et al. (2008) and Jang et al. (2014). As it can be seen, minor perturbations happen in air layer but they are not strong enough to disintegrate the air layer except for the lowest nominal air layer thickness in both free-stream velocities. The disintegration occurs in the lowest nominal thickness for both free-stream velocities at locations of about 5.5–7 m downstream the flat plate leading edge which is not in agreement with experimental data. Jumps of the air layer interface position downstream of the disintegration point are due to the genera­ tion of bubbles and patches of air adjacent to and detached from the wall (pay attention to the scales of plot axes which make bubbles shapes similar to high frequency noises). The nominal thickness of air layer is not the same for simulated cases in these two free-stream velocities but two plots in this figure show qualitatively similar behavior, although the magnitude of free-stream velocities are quite different. Fig. 2 shows the time-averaged shape of air layer for these cases. It can be seen that the averaged shape of the air layer for cases with close nominal air layer thickness coincide with each other within acceptable error margin. It also shows that the shape of the air layer in similar flows depends on nominal air layer thickness and hence the air layer thickness is proportional to the air flux and the in­ verse of free stream water velocity. A close approximation of the required base flow data for the third step in the hybrid numerical algorithm, i.e. the stability analysis, can be obtained using these stable air layer flow simulations of the first step. Thus the air layer thickness and velocity distribution should be extracted from simulations of this step. The second step of the procedure is the curve fitting. Equation (1) is the curve fit template that Montazeri and Alishahi (2019) used. They compared the results of linear stability using this velocity distribution with those of linear velocity distribution (which is more common in literature) and showed that the correct estimation of base velocity dis­ tribution is very important as the stability spectrum using these two conditions is very different. In this equation um is a mean velocity of air layer flow, u0 is the interface velocity, x is the downstream position, y is the distance from wall and d is the air layer thickness. j-subscript

2. Hybrid numerical method Montazeri and Alishahi (2019) proposed a four-step numerical pro­ cedure for simulation of ALDR with an acceptable numerical cost and accuracy. These four steps are: 1 A simple 2-dimensional CFD simulation with appropriate two-phase model, URANS turbulence model and other numerical methods should be carried out. Then some quantities such as air layer thick­ ness and velocity distribution of flow in boundary layers should be extracted from CFD results. 2 A curve fit must be applied on the velocity distribution of CFD results with proper (and preferably physical) template as introduced by Montazeri and Alishahi (2019). 3 Using the ensued velocity distribution from the 2nd step as the base velocity, a two-phase linear stability analysis with appropriate 3

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

Fig. 1. Air layer shape predicted by CFD simulations without additional perturbations (first step) for free-stream velocities of 11.1 m/s (left) and 6.7 m/s (right).

upstream of the injection region. Montazeri and Alishahi (2019) used simple exponential decay because the air-water interface has just a minor effect on velocity distribution and it was supposed that an expo­ nential decay term can lead to accurate enough results. Although the exponential decay term succeeded in accurate curve fit, in this paper the water stream velocity distribution is modified to have a more physical shape by changing this exponential term to error function term similar to the other parts of the velocity distribution. The effect of this modifica­ tion on the analysis result is shown in the following. The new modified velocity distribution is given in equation (2). � �1 �1 0 � 0 d y y erf erf B B A1 C A2 C B B C � � A þ ðu0 um Þ@1 � �C ujet ¼ um @ d d A erf erf A1 A2 � � � �� � y yp y uenv ¼ U∞ erf þ ðu0 ub Þ 1 erf A4 A3 � � yp (2) ub ¼ U∞ erf A3 sffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi μj C1 x μj C 2 x A1 ¼ π ; A2 ¼ π ρj um ρj jum u0 j sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi μe C3 x þ xinj μe C4 x ; A4 ¼ π A3 ¼ π ρe U ∞ ρe ðu0 ub Þ

Fig. 2. Time-averaged air layer shape predicted by the first step CFD (no additional perturbation) for U ¼ 11.1 m/s (black) and U ¼ 6.7 m/s (grey) freestream velocities.

indicates the air layer properties and e-subscript indicates the water stream ones. C1 to C4 are constant coefficients that must be tuned for a good fit. � �1 �1 0 � 0 d y y erf erf B B C C A A 1 B � � C � 2 �C ujet ¼ um B @ A þ ðu0 um Þ@1 d d A erf erf A1 A2 � � �� (1) y uenv ¼ up y expðC4 yÞ þ U∞ þ ðu0 U∞ Þ 1 erf A3 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi μj C 1 x μj C 2 x μe C 3 x A1 ¼ 2 ; A2 ¼ 2 ; A3 ¼ 2 ρj u m ρj jum u0 j ρe ðu∞ u0 Þ

This template uses the Stokes-like profile for both pre-injection boundary layer and air-water interface in post-injection region, each with its appropriate downstream position (x). In this version of equa­ tion, a set of new parameters (ub, yp) has been introduced. The preinjection boundary layer is assumed to shift up in y-direction by a value of yp at locations downstream of injection point. ub is the velocity of water flow inside the boundary layer due to flat plate wall in noinjection condition at location of (x, yp). xinj is the length of flat plate prior to the injection point. In previous version of the equation used by Montazeri and Alishahi (2019), equation (1), the equality of shear stress at the interface sides determined up, but in this new version this equality fixes yp instead. In this modified version of the template, the first term in uenv is the boundary layer of the flat plate in no-injection condition. In this term the xinj is added to x to represent the effect of pre-injection boundary layer growth. Also, the y-coordinate is shifted by yp which means that the pre-injection boundary layer shifts by yp in y-direction in post-injection locations. The amount of this shift (yp) is determined by the condition of shear stress equality at the interface sides. The second term of uenv is the effect of air-water interface which also assumed to grow similar to Stokes profile. As the first term results a velocity at interface, the second term uses u0-ub as its governing velocity to have the correct interface velocity at y ¼ 0. Using the new modified template and tuning the constant co­ efficients, the best curve fit on CFD results was obtained. Fig. 3 shows

This special velocity distribution consists of a specific function for each phase; the air layer velocity distribution (ujet) is a combination of two boundary layer profiles similar to Stokes’ first problem solution. This type of velocity distribution was used by Soderberg (2003) and later Montazeri and Alishahi (2019) showed that adding a constant coeffi­ cient to the formula can also represent an accurate boundary layer profile for flat plate in compliance with the well-known Blasius velocity profile. The water stream velocity distribution (uenv) consists of two parts but only one of them is the Stokes-like profile. The other part in uenv is an exponential decay term corresponding to the effect of air-water interface imposed on the previously developed boundary layer on wall surface 4

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

this fitting of the results for three air fluxes in free-stream velocity of 11.1 m/s. As can be seen, the fitted curve matches well on CFD results points. As the previous template was also shown to match the CFD data by Montazeri and Alishahi (2019), it seems that this modification does not provide any sensible improvement. However, the main difference of these two versions is in the shape of boundary layer of air-water inter­ face in water stream and the new modified version is expected to be more physical than the previous exponential decay form as it applies a form of Stokes’ problem solution (a physical solution) but the previous version used a simple exponential form (just a mathematical fit) with no physical background. As mentioned before, both versions of velocity distribution templates compromise a proper curve fit and many other functions can also be fitted on velocity distribution data from CFD results with good accuracy. Also, introducing a complex velocity distribution to have an accurate curve fit may be questionable. The key point is that the first and second derivative of the velocity distribution are present in eigen system matrix equations and play important roles in the stability analysis. The presence of these derivatives can be seen in Orr-Sommerfeld equation (equation (3)) and its boundary conditions (equation (4)) at air-water interface location. In these equations U is the base velocity and j and e subscripts mean values corresponding to air jet and water stream environment sides of interface. ψ and φ are perturbation stream functions of air and water stream respectively. k and c are non-dimensional wave number and wave speed. m is the water to air viscosity ratio and Re and Fr are Reynolds and Froude numbers based on air properties and air layer thickness reference length. � k2 ψ ðU

ψ ’’

ψ U’’ ¼

cÞ

iν kνj Re

k4 ψ þ 2k2 ψ ’’

ψ ’’’’

�

When these equations are solved numerically using the Chebyshev polynomials, values of velocity derivatives are derived by direct differ­ entiation of the fitted velocity curve. Since in a normal curve fit process, one only tries to lower the difference between values of the curve fit and the velocity data points (not paying attention to the velocity de­ rivatives), the usual curve fit of the base velocity would not provide reliable stability analysis results. One of the proper methods (if do not call it the best) to provide correct velocity and its derivatives informa­ tion simultaneously, is to use a function that represents the physical aspects of the flow. In this manner, the derivatives accuracy only de­ pends on fitting accuracy of the velocity data, as the derivative infor­ mation is previously contained in the velocity template. To clarify the difference of these two versions of velocity distributions, first and sec­ ond derivatives of the fitted curve using both versions of template is compared in Fig. 4. As can be seen, although the curve of velocity distribution using two versions are not sensibly different, shapes of the derivatives are completely different. As the shear stress equality is solved for both versions, the first derivative of velocity at origin is identical for two distributions but there is one additional maxima in the new version derivative (the origin is put on the interface as was shown by Montazeri and Alishahi (2019).) As these derivatives are present in the governing differential equation and its boundary conditions, they are crucially important in finding correct stability spectrum. It is noteworthy that why this point has not received enough attention from previous re­ searchers. That might be due to simple boundary conditions that they employed straightforward velocity distribution (uniform or linear ve­ locity profile) for their specific stability analysis. The 3rd step in the above solution procedure is the linear stability analysis. Two-phase Orr-Sommerfeld equation was solved numerically using sufficient number of Chebyshev polynomials. Boundary conditions of perturbations was set at flat plate wall, free-stream virtual wall and instantaneous air-water interface location. This analysis was carried out just similar to Montazeri and Alishahi (2019) except the modification in the template formulation of base velocity distribution as explained above. Fig. 5 shows results of linear stability analysis using both versions of the template formulation (previous version at left and new modified version at middle). This figure shows that the most unstable frequency for different air fluxes are changed in the new analysis version. Also the magnitude of growth rate is changed to a large extent. It is noteworthy to mention that all of these differences are due to the minor change in formulation of velocity distribution with no significant change in curve fit accuracy. Fig. 5 also shows CFD results of Montazeri and Alishahi (2019). As mentioned before, they performed many CFD simulations with different perturbation frequencies. Different instabilities with different power and tendency of the air layer toward disintegration were found almost in all cases. They used continuous air layer length as the variable that reveals these different flow conditions and tendencies.

(3)

ψ ¼φ ψ’

U ’j c

φ’ þ ψ �

ψ ’’ þ ψ

’’

Uj c

U ’e ¼0 1

� ’’ � � �� Ue þ k2 ¼ m φ’’ þ φ þ k2 1 c 1

"

!

k3 φ’

2S U ’j

U ’e

!# þ ψ’

hU ’ Þ þ φ’ hU þ

U’ Þ

ψ’ U þ

2S U ’j

!

Fr U ’j

U ’e !#

Fr U ’j

U ’e

U ’e

" � � � �� � 3mi 3i þ k2 φ’ þ ψ’ þ k φð Re Re

ψð ¼ φ’ ðωhÞ

ψ ’ ðωÞ þ

i ’’’ ψ Re

mi ’’’ φ Re (4)

Fig. 3. Curve fitting of new modified version (equation (2)) of template on CFD data for 11.1 m/s free-stream velocity and three different air fluxes- Qa ¼ 0.04 m2/s (left), Qa ¼ 0.06 m2/s (middle) and Qa ¼ 0.08 m2/s (right). 5

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

Fig. 4. Effect of modification in curve fit template for velocity distribution and its derivatives – distribution of velocity (left), first derivative (middle) and second derivative (right).

Fig. 5. Stability analysis results using previous version of base velocity template (left) and new modified version (middle) and results of CFD simulations with the addition of perturbation with various frequencies (right).

Stability analysis results and CFD ones in Fig. 5 can be compared by comparing the most unstable frequencies. The agreement of the present linear stability analysis with Montazeri and Alishahi (2019)’s CFD simulations are better than the previous linear analysis. This can be seen specifically for 0.06 m2/s air flux where the CFD simulations show that the flow experiences its most unstable state when it is perturbed at 200 rad/s frequency and the linear stability analysis using the new velocity distribution predicts 240 rad/s as the most unstable frequency while the old velocity distribution forecast was 340 rad/s. This improvement again shows the fact that the base velocity should be set as accurate as possible by correct modelling of base flow field physics. The last step in the numerical simulation procedure is to repeat the first step CFD simulations with the addition of harmonic perturbation in boundary condition with the frequency of the most unstable mode of flow extracted from linear stability analysis. The perturbation is added to the inlet boundary condition of air as a fluctuation of velocity about its mean value. The air layer shape predicted in these simulations are plotted in Fig. 6. The air layer is completely different from what was shown in Fig. 1 due to the mentioned perturbations added in the air injection boundary condition. The effect of the fluctuation amplitude was also investigated by comparing 1, 3, 5 and 10% values. The conclusion is that there is no significant difference between these perturbation results. This implies that the fluctuations in air layer interface is mainly related to the flow instabilities and not to the magnitude of forced perturbations at the inlet boundary. The added velocity fluctuation at the boundary only triggers the unstable growth. Thus, the CFD simulations were carried out with

3% perturbation of air velocity about its mean value. In practice, a perturbation with 3% fluctuation amplitude seems to be physical due to the compressor and other mechanical and environmental vibrations and noises. As is shown in Fig. 6, the predicted unstable flow shows very large fluctuations compared to this initial amplitude which confirms that the fluctuation is due to the instability of the flow rather than an initial forced perturbation.

Fig. 6. Air layer shape predicted by the last step in numerical procedure (CFD simulation with the addition of appropriate perturbation) for different nominal air layer thicknesses in 11.1 m/s free-stream velocity. 6

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

3. Results and discussion 3.1. Investigation of ALDR flow field Drag reduction results of the hybrid numerical simulations at a sensor point located at 6 m downstream of the injection point are plotted in Fig. 7 alongside the experimental results from Elbing et al. (2008). Results of the mentioned procedure shows acceptable accuracy in comparison with experimental data. Pay attention to the difference between the first (simple numerical simulations with no additional perturbation) and forth step results where the first step simulations with no additional perturbation predict 100% drag reduction even for low air fluxes but the addition of the perturbations (step 4) brings the results to an acceptable range from the experiment. Another test case is also studied here with the injection of air behind a backstep (Elbing et al. (2013)’s test case geometry and conditions). Elbing et al. (2013) reported the average of three drag reduction ratio measured by three skin friction sensors located at three different downstream locations. Similar results were extracted from numerical simulations and are plotted in Fig. 8. The error of the drag reduction ratio predicted by the numerical simulations is in quite satisfactory range from experiment. Elbing et al. (2013) also reported velocity and void fraction distri­ butions. These data were reported as the average of three sensors output at three different positions (distributed along the total flat plate length) for different vortex generator configurations (including the case of no vortex generators). Similar results were extracted from numerical sim­ ulations without any vortex generators and are plotted in Fig. 9. As the computed flow field contains different time-dependent fluctuations, the numerical data are presented as time-averaged values in a relatively large time period of about 3 s. In this figure, the experimental data of void fraction are plotted for two different air fluxes. Elbing et al. (2013) stated that the drag reduction for these two air fluxes are above 80% and accordingly, the air layer development should be complete in these cases. Fig. 9 shows the velocity distribution of both experimental and nu­ merical results. Elbing et al. (2013) measured velocity of the water flow using LDV in regions of almost no void fraction. Velocity distribution in regions of higher void fraction where both air and water phases are present was measured by time-of-flight sensor. The time-of-flight sensor measures the velocity of two phase flow using its phase-interface ve­ locity but it does not provide any data when only one phase is present. Thus, as Elbing et al. (2013) also stated, the measurements of velocity in regions with high void fraction is not as accurate as other regions. Therefore, the anomaly in experimental void fraction close to the wall is somehow expectable. As Fig. 9 shows, the numerical water velocity distribution is in complete agreement with experiment but the two-phase velocity

Fig. 8. Drag reduction of air injection behind a backstep at 6.3 m/s freestream velocity (Elbing et al., 2013 test case).

distributions are different to some extent. Elbing et al. (2013) proposed a three layer velocity distribution in ALDR; the first part adjacent to the wall filled mainly with air and a constant velocity about 20% of free-stream velocity, the third part in water-only region with velocity distribution very similar to no-injection boundary layer distribution, and a buffer layer (second part) between these two parts with almost linear velocity variations. However, there is a shortcoming with this proposal as it doesn’t account for the value of air flux. It is noticeable that the Elbing et al. (2013) proposal is stated only for air fluxes that result in complete air layer development. On the other hand, the numerical ve­ locity distribution can be described as it consists of only two parts; one is the outer region where the velocity distribution is very similar to no-injection boundary layer (just like the experimental data acquired by LDV measurements), and the other is the inner region where the air layer is qualitatively similar to an unsteady Couette flow velocity profile in starting times of flow. The velocity distribution in this inner part consists of two boundary layers on the wall surface and air-water interface and may gradually change toward a fully developed profile if the air layer remains stable and persists without significant change in thickness. Later in this section in Fig. 10, variation of velocity distribution at different downstream positions is presented and discussed. Although the location of void fraction rise in numerical and experi­ mental results are close as shown in Fig. 9, there are two main differ­ ences in these results; one is a minor void fraction rise up to 20% in both experimental data while it is not present in the numerical output. The other difference is a drop in experimental void fraction at wall adjacency especially for air flux of 0.053 m2/s. The first difference may be due to some small free bubbles above the air layer in experimental test and seems not to be so important. The drop in void fraction at wall adjacency in experimental data is not in agreement with the physics of a contin­ uous full air layer. The existence of an air layer by itself implies a continuous stratified air layer adjacent to the wall and consequently a void fraction of unity at the wall. This strange behavior in Fig. 9 can be seen mainly in the lower air flux and the data for the higher air flux are more reasonable. It is noteworthy that according to Elbing et al. (2013), the error in void fraction measurement is up to 15% in regions of high void fraction (namely inside the air layer). Note that both air fluxes result in high drag reductions in the experiments and expressed as developed air layer by Elbing et al. (2013). To get a better insight of the flow field, these velocity and void fraction distributions are shown in Fig. 10 at each position separately (not the average of three positions) for three different air fluxes (these data are also the time-averaged numerical results). In this figure, the horizontal axis is the normal distance to the wall scaled by nominal air layer thickness instead of no-injection boundary layer thickness. As can be seen in the top row of Fig. 10, the void fraction is unity close to the wall and at the first station (closest to the injection point) for all three air fluxes. This shows the existence of an air layer even for low

Fig. 7. Drag reduction of air injection from simple slot at 11.1 m/s freestream velocity (Elbing et al., 2008 test case). 7

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

Fig. 9. Comparison of velocity and void fraction distribution predicted by numerical simulation and Elbing et al. (2013) experimental test measurements results.

Fig. 10. Normalized velocity (dotted line) and void fraction (solid line) distribution for different nominal air layer thicknesses: ta ¼ 3.6 mm (left), ta ¼ 5.4 mm (middle) and ta ¼ 6.3 mm (right) at positions with different distances from injection.

air fluxes in the vicinity of the injection point. At this station, a relatively sharp interface between air and water can be inferred from large slope of void fraction distribution. At the second station, i.e. the middle row, it seems that the air layer thickness is decreased and the interface of airwater is not as sharp as the first station. In other words, the air layer bears some time fluctuations and the air-water interface moves up and down resulting in a void fraction value between 0 and 1 from the time averaging process and hence a decrease in void fraction slope would be seen. Also this figure shows that there is almost no air layer in second position for the lowest air flux (ta ¼ 3.6 mm), i.e., the air layer is ruptured, as the void fraction is not unity at the wall. Further down­ stream, at the third station, the air layer is also ruptured for medium and high air fluxes (ta ¼ 5.4 and 6.3 mm). As can be seen in this figure, when the air layer disintegrates and the void fraction drops from unity to a lower value at the wall, there is a maximum in void fraction distribution near the wall but not at the wall in all cases. In other words, when the air layer disintegrates into bubbles

and air patches, they start to move away from the wall surface and flow mostly along a stream parallel to the wall but in a distance away from it. This is in agreement with the finding of some previous researchers including Kodama et al. (2000) and Merkle et al. (1990). They stated that air bubbles tend to escape from the wall surface in a boundary layer flow. This is also consistent with the earlier mentioned drop in void fraction at the wall in Elbing et al. (2013) measurements for the lower air flux in Fig. 9. They stated that at this air flux a continuous air layer would prevail throughout the flat plate and high drag reduction ratio is resulted. In contrast to their reasoning, if the present numerical results are accepted, the void fraction drop in Elbing et al. (2013)’s test data for 0.053 m2/s air flux should be due to a degenerated air layer and we are faced with a ruptured air layer instead of continuous one. In Fig. 8 the measured drag reduction ratio for the cases of Qa ¼ 0.06 and 0.08 m2/s (ta ¼ 0.01 and 0.013 m) is more than 80% and probably these high drag reduction ratios suggested the existence of continuous air layer in Elbing et al. (2008) and Elbing et al. (2013). However, the present numerical 8

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

simulations show that despite of the air layer rupture, void fraction at the wall remains high and consequently the drag reduction ratio is also high. It is important to note that for a wide range of air fluxes, the ex­ istence of air layer with different continuous length does not have any conflict with raw experimental measurements, in contrast to some re­ searchers’ misinterpretation of their measurement. High values of drag reduction without a continuous air layer is also reported in an experi­ mental test on flat plate by Park et al. (2015). They used pulsatile in­ jection instead of continuous steady injection of air and hence resulted in discontinuous air layer. Repetitive injection of air leads to transitional ALDR which can be assumed to be similar to the ruptured air layer as both contain large bubbles and air patches with high void fraction near the wall. In their experiment, they showed that this type of flow brings about high drag reduction ratios similar to continuous air layer. Their finding confirms that high drag reduction values can be achieved even when there is no continuous air layer. There is also another point in Fig. 10 void fraction plots; it can be seen that the air layer has an initial thickness (thickness of regions with more than 50% void fraction) and this thickness decreases downstream. The air layer disintegrates at a position which is different for different air fluxes. To further investigate the disintegration point of the air layer, the time-averaged air layer thickness is plotted in Fig. 11 for different nominal air layer thicknesses. The ordinate in this figure is the normalized air layer interface position which is the ratio of actual timeaveraged interface position to its nominal thickness. The figure repre­ sents results for different air fluxes, for both Elbing et al. (2008) and Elbing et al. (2013) test cases (without and with backstep implementa­ tion upstream of the injection point). Air layer thickness in cases containing the backstep starts from the step height and gradually decreases downstream. Note that the step height is the same in different cases but the air flux changes and hence the normalized step height is different for different cases. In cases without the backstep, air layer starts from zero thickness and a growth in thickness near-injection location can be seen which is due to the injec­ tion angle and air-water stream momentum balance. Following this increase in air layer thickness, the air layer continues downstream with a gradual reduction in thickness. The air layer fluc­ tuates before disintegration point but as the fluctuation is nearly sym­ metric, the time-averaged air layer thickness does not change significantly due to these fluctuations. For all cases that are plotted in this figure (including cases with and without backstep), the air layer thickness jumps to a lower value at some specific point which is the disintegration point of air layer. When the air layer disintegrates into bubbles and air patches, they flow inside the boundary layer and often touch the wall surface. These touches make the time-averaged air layer

thickness to have a non-zero value even when there is no actual continuous air layer. After the disintegration point, there are many bubbles and air patches which many of them do not touch the wall and so they have no effect in time-averaged air layer thickness. Hence, at the disintegration point the averaged thickness should suddenly drop to a small value which can be seen for almost all cases. However, in Fig. 11 it can be seen that higher air fluxes result in smoother change of air layer thickness at disintegration point which may be due to higher void fraction in ruptured air layer region and many contacts of bubbles and wall surface in that region. In Fig. 11, there is a remarkable topic; the horizontal thick light grey line shows normalized interface position of 1 which means that the actual time-averaged thickness of air layer is equal to its nominal value. Almost all of studied cases show a change in their interface position slope at this line. Hence, the point of disintegration for these different cases have interestingly almost identical normalized air layer thickness equal to unity. This means that the air layer thickness is equal to its nominal thickness at the disintegration point. In other words, the air layer thickness decreases gradually from an initial value and when the time-averaged thickness reaches its nominal thickness (Qa/U), the air layer disintegrates suddenly. This phenomenon can be seen in all cases including cases with backstep and cases without backstep. Fig. 11 brings about the hint that the air layer seems to disintegrate at a point where the air layer thickness equals its nominal value. This may suggest the simplification of numerical procedure, using the air layer thickness result from simple non-perturbed simulation to find the disintegration point without accomplishing the forth step CFD simula­ tion. Fig. 1 shows the air layer thickness for non-perturbed simulation (first step) while Fig. 11 shows the time-averaged of air layer thickness for perturbed simulation (forth step). As the added perturbations are harmonic and symmetric, the time-average of perturbed flow may be expected to be similar to the non-perturbed flow at least before the disintegration. Paying more attention reveals that the gradual thinning of the air layer in Fig. 11 is different from what can be seen in Fig. 1. The addition of perturbation in inflow velocity boundary condition seems to change the averaged air layer shape although the perturbation is har­ monic and symmetric (maybe due to nonlinearities of the flow). This shows that the numerical simulation of the flow without the additional perturbations does not predict even the averaged air layer thickness before disintegration and performing the perturbed simulation is essential. Fig. 11 can also be used to investigate the effect of backstep imple­ mentation on the performance of ALDR. Although the nominal air layer thicknesses are not exactly identical for two cases of with and without backstep, it can be seen that for cases with almost the same nominal

Fig. 11. Time-averaged air layer interface position normalized by nominal air layer thickness for different cases. 9

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

thickness, the existence of the backstep causes a fixed initial air layer thickness and the air layer thickness decreases gradually toward the thickness of no-backstep case. The rupture of air layer in both cases are almost at the same position except for low nominal thickness case. This means that the existence of the backstep cannot postpone the disinte­ gration point of the air layer to have a higher performance. This is in agreement with the experimental data of Elbing et al. (2008) and Elbing et al. (2013) as explained in the following from another point of view. As can be seen in Fig. 12, the drag reduction ratios reported by Elbing et al. (2008) without implementation of backstep at free-stream velocity of 6.7 m/s is almost the same as the data reported by Elbing et al. (2013) using backstep at free-stream velocity of 6.3 m/s. This shows that the backstep cannot increase the drag reduction performance of ALDR, although it may lower the sensitivity of air layer to upstream pertur­ bations. The disintegration point for low air fluxes is relatively near the injection point in no-backstep cases and the step height is much larger than the air layer thickness in identical case using backstep. Present numerical simulations show that large initial thickness due to the backstep in these cases may send the disintegration point further downstream (as can be seen for nominal thickness of 3 mm in Fig. 11). This is favorable if low air flux is used. More than one point of injection may be chosen in this situation to accomplish the air layer on the whole hull surface, although it is not so advantageous as it increases the no-injection drag and the overall air flux may also be high.

nominal air layer thickness would be unacceptable. For lower bumps on the wall surface, the velocity distribution of the air layer may change. New velocity distribution might increase or decrease the air layer thickness depending on the shape of the new distribution. This means that the change in velocity distribution (and whether to increase or decrease the air layer thickness), may have a stabilizing or destabilizing effect on the flow. As the first case of this study, circular curves with the length of 0.6 m were placed on the wall surface. The numerical simulation was carried out at 11.1 m/s free-stream velocity and 0.04 m2/s air flux. Two bumps were inserted on the wall; one at 1.5 m and the other at 3 m downstream of the injection point. Results for this flow without any bumps showed that the air layer degenerates at about 1.8 m downstream of the injection point. Destabilizing or stabilizing effect of bumps may cause the air layer to disintegrate before (destabilizing effect) or after the first bump or even the air layer may survive to the second bump location (stabilizing effect). As the nominal air layer thickness is about 3.6 mm for the studied case, three sizes of bumps were investigated with heights of 1, 2 and 3 mm. Instantaneous and time-averaged results of the simulation are shown in Fig. 13 and indicate a different conduct than what is expected. The air layer thickness change is insensible at bumps location and also the fluctuations in air layer do not increase or decrease due to bumps. It shows that this type of smooth bumps does not have any favorable or unfavorable (stabilizing or destabilizing) effects on the air layer. Another configuration of bumps is also studied in this investigation. A triangular section was considered with a small slope in upstream di­ rection and a backward facing step downstream, Fig. 14. Similar to the smooth bumps, three sizes of 1, 2 and 3 mm height and 0.6 m length of the triangular bumps were studied and they were placed at the same positions of previous bumps. To check if the position of the protrusions changes their effect, the triangular bumps with 1 mm height were also placed at new position of 0.9 and 2.4 m downstream of injection point in another simulation. It might be expected that the triangular bumps initiate a separated flow at their backstep base and thus the air layer thickness increases. The instantaneous and averaged air-water interface location is plotted in Fig. 14 for these four configurations. As can be seen again, the triangles do not have any noticeable effect on air layer (similar to smooth bumps).

3.2. ALDR sensitivity investigation As mentioned earlier, ALDR is capable of hydrodynamic frictional drag reduction to a great extent. In the previous section it was shown that a continuous and stable air layer may cancel out the frictional drag by almost 100%. Although high drag reductions are also present in cases of ruptured air layers for relatively high air fluxes, it is more favorable to maintain a continuous air layer instead of a ruptured layer which in­ cludes bubbles and air patches. Ruptured air layer may gradually loose the tendency of high drag reduction when bubbles move away from the wall surface and void fraction in boundary layer decreases. Thus, any disintegration of air layer should be assumed as non-favorable. One of the noticeable impacts that an air layer may encounter is the existence of wall surface bumps. During the construction of a ship’s hull, there would be many bumps on the surface due to manufacturing errors and tolerances. If the air layer is so sensitive to these bumps, disinte­ gration of the layer may occur which might construct the most impor­ tant weakness of the ALDR when applied practically. Thus, the air layer sensitivity to bumps on the flat plate surface is studied in this section. Referring to Fig. 11, it seems that the nominal air layer thickness would be the critical value for time-averaged air layer thickness such that the disintegration of air layer would occur if its time-averaged thickness reaches this critical value. Hence, bump heights more than

3.3. Enhancement of ALDR performance Performance of drag reduction methods based on air injection can be evaluated via gain factor (ratio of overall drag reduction to the required air flux). Any increase in drag reduction or lowering the required air flux improves the efficiency of the process. Regarding this, some parameters are changed to check whether the performance of ALDR can be enhanced. One of the parameters that might play a role in the instabilities of this flow (the Kelvin-Helmholtz instability), is the surface tension of water. In standard conditions, the surface tension coefficient (σ) of water-air interface is about 0.07 N/m but it can be changed by addition of some soluble liquids into the water. Injection of surfactant into the water flow is shown to change the surface tension coefficient by researchers, e.g. Elbing et al. (2008). Elbing et al. (2008) injected Triton-X-100 into water stream and their measurements showed that the surface tension coeffi­ cient was reduced from 0.07 to 0.05 N/m. They investigated the effect of this change on bubble drag reduction (BDR). They showed that the change in bubble drag reduction and bubble size due to the change in surface tension is not significant. They also measured the critical air flux for both normal water and water with additional surfactant. Their result shows just a minor difference between these two conditions. If this minor difference is assumed not to be the measurement error, the lower surface tension is showed to result in higher critical air flux by Elbing et al. (2008). In the present study, the effect of surface tension is studied for ALDR. Numerical simulation of the ALDR flow field for 0.04 m2/s air flux and 11.1 m/s free-stream velocity was reproduced with surface

Fig. 12. Comparison of Elbing et al. (2008) data of simple air injection from slot and Elbing et al. (2013) data of air injection behind a backstep. 10

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

Fig. 13. Air layer thickness for cases with different smooth bump heights – instantaneous shape (grey solid) and time-averaged (black dots): no-bump (top left), 1 mm bump (top right), 2 mm bump (bottom left) and 3 mm bump (bottom right).

Fig. 14. Air layer thickness for different triangular bumps – instantaneous shape (grey solid) and time-averaged (black dots):1 mm bump (top left), 2 mm bump (top right), 3 mm bump (bottom left) and 1 mm bump positioned 0.6 m upstream of others (bottom right). 11

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

tension coefficients of 0.05 and 0.1 N/m. Surface tension coefficient is present in linear stability analysis in the phase-interface boundary condition i.e. the last boundary condition in equation (4). On the other hand, surface tension is considered in numerical simulations as a force proportional to its coefficient and phase-interface surface curvature in momentum equation for cells which contain the interface. Fig. 15 shows the results of these two cases alongside the original simulation with surface tension coefficient of 0.07 N/m. As can be seen, the change in surface tension does not produce any visible enhancement in air layer. In order to make the air layer more stable or to increase its contin­ uous length without any rupture, the air layer thickness should be increased. If the air layer thickness increases, in any other means, the disintegration point might be pushed back downstream. The air layer thickness can be raised by increasing the air flux; but it is not desired as the efficiency of the ALDR decreases. The other way to expand the air layer is to change the air velocity distribution bounded by mass con­ servation of air. This change in velocity distribution may be brought up by a change in wall shear rate probably due to non-smoothness. Velocity distribution of water in contact with the wall can produce an effect on the interface velocity distribution and even on the air layer velocity distribution. All of the previous CFD simulations in this paper were carried out assuming smooth surface condition. Non-favorable effect of increasing the surface roughness on ALDR is reported by Elb­ ing et al. (2008) as the required air flux to establish the air layer in­ creases. Thus, the hydrophobicity of the wall surface is investigated in this study as it might be expected that the hydrophobicity has an opposite effect comparing with the wall roughness. A hydrophobic surface is a surface with a material that repels water molecules. Having such a material, the hydrophobicity of a surface can be enhanced by the use of micro/nano roughness. Numerical modelling of flow on hydro­ phobic surface is usually done in two steps; first, a very small numerical domain with the specific roughness on the wall is simulated using DNS or at least LES methods. Then, an equivalent Navier slip length is esti­ mated from such a simulation. Vast range of slip lengths is reported in experimental and numerical researches from very small values (some nanometers) for simpler hydrophobic surfaces to large values up to 0.4 mm for ultrahydrophobic surfaces with hierarchial wall roughness structures (Lee and Kim, 2009). Slip lengths of some tens of micrometers is more common in literature as large and yet practical values. Song, Daniello and Rothstein (2014) reported 0.02 mm and Li et al. (2009) reported 0.04 mm slip length in their measurements. Larger slip lengths are also reported in academic researches, i.e. 0.1 mm by Lee and Kim (2011), 0.05–0.25 mm by Maali and Bhushan (2012) and 0.4 mm by Lee and Choi (2008). Production of ultrahydrophobic surfaces with such a large slip length in practical applications has not been accomplished, yet; but perhaps it would be achieved in the future. Nowadays, super­ hydrophobic surfaces with slip lengths of tens of micrometers are con­ structed and used in practice. In this paper, a Navier slip length value is used and it is assumed that a surface with some specific material and

roughness structure and size can be constructed to have such a Navier slip length value. The Navier slip length in ALDR solution was implemented by setting up a slip velocity at the wall surface. This slip velocity was set to zero for cells with void fractions greater than 0.5 (assumed to be air-filled cells) that is usual no-slip condition. On the other hand, cells with void frac­ tions less than 0.5 are assumed to be water-filled and the hydrophobicity effect should be set. In these cells the slip velocity was set equal to

τw l μ

where τw is the wall shear stress and l is the slip length. As the wall shear stress depends on wall slip velocity, enough iterations should be per­ formed at each time step to ensure convergence of wall shear stress and slip velocity. Some slip lengths ranging from 0.02 mm to 2 mm were considered and ALDR flow was simulated for 11.1 m/s water free-stream velocity and 0.04 m2/s air flux. Fig. 16 shows the length of continuous air layer in these cases with different Navier slip lengths. Also, Fig. 17 shows the air layer shape for three cases among these. The hydropho­ bicity of the wall makes the air layer more stable and the continuous air layer length is increased from 1.6 m to about 3.5 m in cases with larger slip lengths but it does not provide major enhancement in cases with small slip lengths. As can be seen in Fig. 16, in the studied flow condi­ tions, 0.4 mm is the critical slip length and air layer enhancement is evident for larger slip lengths while minor improvements can be seen for smaller values. As mentioned before, achieving the slip lengths larger than this critical value is not straightforward, yet; but it seems not to be impossible. This critical slip length could be different for other flow conditions (water stream velocity and air flux). The next test case was carried out only for 1 mm slip length. This test case with 0.06 m2/s air flux was chosen to ensure if such an improve­ ment exists for other conditions. Fig. 18 shows the result of this case and similar improvement can be seen for 0.06 m2/s air flux case. In this combined application of air layer method and hydrophobic wall surface, if the air layer develops and continues to the end of flat plate (or practically, hull wall surface), almost no water-wall contact would be present. In such cases the hydrophobicity of wall would have no positive effect as the wall does not have any contact with water. Thus the hydrophobicity should affect the flow at regions with no air layer, i. e., regions prior to the injection point. Thus the difference of water boundary layer velocity distribution at injection point between the normal case and hydrophobic surface might be the origin of the stabi­ lizing effect if there is such an effect. To investigate this statement, previous simulation is repeated with new boundary conditions. Two new cases are studied; one with a hydrophobic surface at the leading edge of flat plate up to 0.5 m upstream of injection point (farfield hy­ drophobic case), and the other case is the nearfield hydrophobic case with a hydrophobic surface just upstream of injection point (from 0.5 m

Fig. 16. Variation of continuous air layer length with slip length of the hy­ drophobic wall surface.

Fig. 15. Air layer thickness for different surface tension coefficient values. 12

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

Fig. 19. Effect of implementation of hydrophobic surface on different regions (far and near upstream of injection point).

Fig. 17. Effect of hydrophobicity on air layer thickness with Navier slip lengths of 0.02, 0.5 and 1 mm – free stream velocity ¼ 11.1 m/s, Qa ¼ 0.04m2/s.

velocities, (low nominal air layer thickness), the air layer encounters instabilities that grow up and bring about disintegration of air layer. Numerical simulation of different cases showed that the air layer dis­ integrates when its time-averaged thickness reaches its nominal thick­ ness. This phenomenon was seen in almost all simulated cases in this paper. However, it can be imagined that for large enough nominal air layer thickness, the instabilities occur in such long distances down­ stream that no effect would be observed in practice. A drop in void fraction at wall adjacency is reported in experimental results of Elbing et al. (2013) for ALDR cases. Present numerical simu­ lations show that this means “ruptured air layer”. For large enough nominal air layer thickness, this ruptured air layer can still result in high drag reduction ratio which may be misinterpreted as continuous air layer drag reduction. For the enhancement of ALDR performance, different parameters effects were studied. It was shown that implementation of backstep just upstream of injection point cannot increase the performance of ALDR unless for small nominal thickness relative to the step height. Changing the surface tension was also non-effective. Covering the wall surface with hydrophobic coating was shown to be a successful option which resulted in a longer continuous air layer. Minor improvement was shown for small and more common slip lengths while major improvement was seen for large slip lengths which may be realized in future. The conclusion is that the main effect of hydrophobicity comes from an upstream region near the injection point. This is an important result as it shows that with a small amount of effort and money for implementation of hydrophobic surface on this small area, the performance of ALDR can be noticeably improved. As the construction and maintenance of a hy­ drophobic surface on small region would be cheap and cost effective, this method should be desirable. Implementation of bumps both in smooth circular and triangular step shapes showed no noticeable effect on the air layer. This means that the air layer sensitivity to minor construction errors and tolerances is negligible and no instabilities or rupture of the air layer would occur in practice. As a general conclusion, the numerical result of this study shows that the ALDR seems to be a robust enough method to be used practically and its combination with hydrophobic surface in a relatively small upstream area enhances its performance remarkably.

Fig. 18. Effect of hydrophobicity on air layer thickness - free stream velocity ¼ 11.1 m/s, Qa ¼ 0.06 m2/s.

upstream of injection point to the point of injection). Both cases have normal no-slip wall condition downstream of injection. Fig. 19 shows the result of these two cases in comparison with the case of totally hy­ drophobic flat plate. As can be seen, the farfield hydrophobic case does not show the total hydrophobic wall improvement but the nearfield hydrophobic case shows very similar result to the reference case. This shows that the main effect of hydrophobicity comes from a region just upstream of injection point. This is an important result as it means that in practice, there is no need to construct the hydrophobic coating on the whole ship’s hull, but only a relatively small region upstream of the injection point is enough to stabilize the air layer for a longer distance. 4. Conclusions Air layer drag reduction is said to be a successful and cost-efficient method for hydrodynamic skin friction reduction. It was studied experimentally in literature but just a few complete numerical (DNS or LES) investigations are available, due to huge computational resource requirements. Montazeri and Alishahi (2019) proposed a numerical procedure to provide acceptable numerical predictions of this flow with low computational cost. This procedure with minor modification was reapplied to various ALDR problems in this study. Numerical results showed that the air layer develops for almost all air fluxes in any free stream water velocity while the length of continuous air layer changes by conditions of each specific case. For low air fluxes in high free stream

Author contributions section Mohammad Hossein Montazeri: Conceptualization, Methodology, Software, Formal analysis, Investigation, Writing - Original Draft. Mohammad Mehdi Alishahi: Conceptualization, Methodology, Resources, Writing - Review & Editing. 13

M.H. Montazeri and M.M. Alishahi

Ocean Engineering xxx (xxxx) xxx

Declaration of competing interest

Li, J., Zhou, M., Cai, L., Ye, X., Yuan, R., 2009. On the measurement of slip length for liquid flow over super-hydrophobic surface. Chin. Sci. Bull. 54 (24), 4560–4565. Maali, A., Bhushan, B., 2012. Measurement of slip length on superhydrophobic surfaces. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 370 (1967), 2304–2320. Madavan, N., Deutsch, S., Merkle, C., 1984. Reduction of turbulent skin friction by microbubbles. Phys. Fluids 27, 356–363. Makiharju, S.A., Elbing, B.R., Wiggins, A., Schinasi, S., Vanden-Broeck, J.-M., Perlin, M., et al., 2013. On the scaling of air entrainment from a ventilated partial cavity. J. Fluid Mech. 732, 47–76. Maltzev, L.I., 1995. Jet methods of gas injection into fluid boundary layer for drag reduction. Appl. Sci. Res. 54, 281–291. McCormick, M.E., Bhattacharyya, R., 1973. Drag reduction of a submersible hull by electrolysis. Nav. Eng. J. 85, 11–16. Merkle, C.L., Deutsch, S., Bushnell, D.M., Hefner, J.N., 1990. Drag reduction in liquid boundary layers by gas injection. Viscous Drag Reduct. Bound. Layers 123, 351–412. Montazeri, M., Alishahi, M., 2019. An efficient method for numerical modeling of thin air layer drag reduction on flat plate and prediction of flow instabilities. Ocean. Eng. 179, 22–37. Murai, Y., 2014. Frictional drag reduction by bubble injection. Exp. Fluid 55, 1773. Park, H.J., Tasaka, Y., Oishi, Y., Murai, Y., 2015. Drag reduction promoted by repetitive bubble injection in turbulent channel flows. Int. J. Multiph. Flow 75, 12–25. Rotte, G., Zverkhovskyi, O., Kerkvliet, M., Terwisga, T., Rotte, G., Zverkhovskyi, O., et al., 2016. On the physical mechanisms for the numerical modelling of flows around air lubricated ships. In: International Conference on Hydrodynamics. Wijk aan Zee, The Netherlands. Sanders, W.C., Winkel, E.S., Dowling, D.R., Perlin, M., Ceccio, S.L., 2006. Bubble friction drag reduction in a high-Reynolds-number flat-plate turbulent boundary layer. J. Fluid Mech. 552, 353–380. Soderberg, L.D., 2003. Absolute and convective instability of a relaxational plane liquid jet, 10 J. Fluid Mech. 493, 89–119. https://doi.org/10.1017/s0022112003005871. Song, D., Daniello, R.J., Rothstein, J.P., 2014. Drag reduction using superhydrophobic sanded Teflon surfaces. Exp. Fluid 55 (8), 1783. Verschoof, R.A., Van Der Veen, R.C., Sun, C., Lohse, D., 2016. Bubble drag reduction requires large bubbles. Phys. Rev. Lett. 117, 104502. Wang, J., Wang, B., Chen, D., 2014. Underwater drag reduction by gas, 12 Friction 2, 295–309. https://doi.org/10.1007/s40544-014-0070-2. Wang, Z., Yang, J., Stern, F., 2010. URANS Study of Air-Layer Drag Reduction in a HighReynolds-Number Flat-Plate Turbulent Boundary Layer. 40th Fluid Dynamics Conference and Exhibit, p. 4276. Watanabe, O., Masuko, A., Shirose, Y., 1998. Measurements of drag reduction by microbubbles using very long ship models. J. Soc. Nav. Archit. Jpn. 53–63. https:// doi.org/10.2534/jjasnaoe1968.1998.53, 1998.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References Butterworth, J., Atlar, M., Shi, W., 2015. Experimental analysis of an air cavity concept applied on a ship hull to improve the hull resistance. Ocean. Eng. 110, 2–10. Daniello, R.J., Waterhouse, N.E., Rothstein, J.P., 2009. Drag reduction in turbulent flows over superhydrophobic surfaces. Phys. Fluids 8, 085103. https://doi.org/10.1063/ 1.3207885, 21. Elbing, B., Makiharju, S., Wiggins, A., Perlin, M., Dowling, D., Ceccio, S., 2013. On the scaling of air layer drag reduction. J. Fluid Mech. 717, 484–513. Elbing, B., Winkel, E., Lay, K., Ceccio, S., Dowling, D., Perlin, M., 2008. Bubble-induced skin-friction drag reduction and the abrupt transition to air-layer drag reduction. J. Fluid Mech. 612, 201–236. Henoch, C., Krupenkin, T., Kolodner, P., Taylor, J., Hodes, M., Lyons, A., et al., 2006, 6. Turbulent Drag Reduction Using Superhydrophobic Surfaces. 3rd AIAA Flow Control Conference. American Institute of Aeronautics and Astronautics. https://doi.org/ 10.2514/6.2006-3192. Jang, J., Choi, S., Ahn, S., Kim, B., Seo, J., 2014. Experimental investigation of frictional resistance reduction with air layer on the hull bottom of a ship. Int. J. Nav. Archit. Ocean. Eng. 6, 363–379. Kanai, A., Miyata, H., 2001. Direct numerical simulation of wall turbulent flows with microbubbles. Int. J. Numer. Methods Fluids 35, 595–615. Kim, D., Moin, P., 2010. Direct Numerical Study of Air Layer Drag Reduction Phenomenon Over A Backward-Facing Step. Tech. rep. Kodama, Y., Kakugawa, A., Takahashi, T., Kawashima, H., 2000. Experimental study on microbubbles and their applicability to ships for skin friction reduction. Int. J. Heat Fluid Flow 21, 582–588. Lee, C., Choi, C.H., 2008. Structured surfaces for a giant liquid slip. Phys. Rev. Lett. 101 (6), 064501. Lee, C., Kim, C.J.C., 2009. Maximizing the giant liquid slip on superhydrophobic microstructures by nanostructuring their sidewalls. Langmuir 25 (21), 12812–12818. Lee, C., Kim, C.J.C., 2011. Influence of surface hierarchy of superhydrophobic surfaces on liquid slip. Langmuir 27 (7), 4243–4248.

14