Exp Fluids (2016) 57:178

DOI 10.1007/s00348-016-2268-8

RESEARCH ARTICLE

‌Void waves propagating in the bubbly two‑phase turbulent boundary layer beneath a flat‑bottom model ship during drag reduction

Hyun Jin Park1

· Yoshihiko Oishi2 · Yuji Tasaka1 · Yuichi Murai1

Received: 1 August 2016 / Revised: 12 October 2016 / Accepted: 21 October 2016

© Springer-Verlag Berlin Heidelberg 2016

‌Abstract The injection of bubbles into a turbulent bound- ary layer can reduce the skin friction of a wall. Conven- tionally, the drag reduction rate is evaluated using time- averaged quantities of the mean gas flow rate or mean void fraction. Actually, as bubbles are subject to strong shear stresses near the wall, void waves and local bubble clusters appear. For pipe and channel flows, such wave- like behavior of the dispersed phase has been investigated intensely as an internal two-phase flow problem. We inves- tigate how this wavy structure forms within the boundary layer as an external spatially developing two-phase flow along a horizontal flat plate. We describe how our model ship is designed to meet that purpose and report bubble- traveling behavior that accompanies unexpectedly strong wavy oscillations in the streamwise direction. A theoretical explanation based on a simplified two-fluid model is given to support this experimental fact, which suggests that void waves naturally stand out when drag reduction is enhanced through the local spatial gradient of the void fraction.

List of symbols

C Void wave propagation speed (m/s) Ca Cavitation number (dimensionless) Cf Frictional coefficient (dimensionless)

dcb Distance between the closest pair of bubbles (m)

de Equivalent diameter of bubble (m)

Fr Froude number (dimensionless)

1. Local instantaneous volume fraction of liquid phase (dimensionless)fG Local instantaneous void fraction (dimensionless)fvoid Frequency of voidage wave (Hz)G1 Impact factor of drag reduction to void fraction in the boundary layer (dimensionless)G2 Impact factor of drag reduction to local gradient of void fraction (m)
2. Acceleration due to gravity (m/s2)H Height of the model ship (m)
3. Thickness of the reflective index matching mate- rial (m)

L Length of the model ship (m)

lτ Friction length (m)

p Local pressure of the water (Pa) pv Vapor pressure of the water (Pa) Qg Injected gas flow rate (m3/s)

Ql Liquid flow rate in the boundary layer (m3/s) Rex Reynolds number on a flat plate (dimensionless) t Time (s)

tg Apparent air layer thickness (m)

Umain Main flow velocity (equivalent to towing speed) (m/s)

Hyun Jin Park

park@eng.hokudai.ac.jp

1 Laboratory for Flow Control, Division of Energy and Environmental Systems, Faculty of Engineering,

Hokkaido University, N13 W8, Kita-ku, Sapporo 060-8628, Japan

2 Fluid Engineering Laboratory, College of Design and Manufacturing Technology, Muroran Institute

of Technology, 27-1, Mizumoto-cho, Muroran 050-8585, Japan

Uδ Averaged flow velocity in the boundary layer (m/s)

u, v, w Velocity components in x, y, z directions (m/s)

ub Averaged advection velocity of bubbles (m/s)

uy Averaged streamwise velocity at each depth (m/s)

1. Averaged downward velocity of liquid phase on the border of the boundary layer (m/s)
2. Width of the model ship (m)

x, y, z Cartesian coordinates of the model ship (m)

αδ Void fraction in the boundary layer (dimensionless)

δ 99% thickness of the boundary layer (m)

δg Superficial air layer thickness (m)

μ Viscosity of water (kg/m s)

ν Kinematic viscosity of water (m2/s)

ρ Density of water (kg/m3)

τw Wall shear stress (Pa)

τw0 Wall shear stress in single-phase flow (Pa)

1. ‌Introduction

   Bubbly drag reduction (BDR) is a collective term for attempts to reduce drag by injecting bubbles into turbulent boundary layers. Over the past four decades, it has received attention as a means to decrease fuel consumption of large marine vessels. Several successes of BDR in sea trials have been recently reported in the literature (Mizokami et al. 2013; Jang et al. 2014; Kumagai et al. 2015). These groups have confirmed independently good reproducibility in fuel- saving performance ranging from 5 to 15% for different types of vessels. In academic fields, drag reduction per- formance and its parametric dependency have been inves- tigated experimentally with in-house flow geometries such as a horizontal channel flows. We have found hundreds of papers on BDR since its first reporting by McCormick and Bhattacharyya (1973). In idealized conditions, local drag reduction rates of 80% can be achieved (e.g., Madavan et al. 1984). Nevertheless, there is a gap between ideal and practical conditions; the recent success for ships relies on expertise obtained with fundamental two-phase flow exper- iments. For further improvements in drag reduction perfor- mance, the behavior of bubbles skimming along a ship’s hull during BDR operations needs clarification. There are two reasons for this requirement. One is the difficulty in extrapolative expectations on how bubbles travel beneath a hull of a real ship from simply laboratory-downsized results. Reynolds number are of much higher order in real ships than those in-house, whereas Froude number can be set by reducing the length scale. The other is that the domi- nant BDR mechanism changes undesirably depending on a combination of dimensionless numbers associated with bubbles, such as Weber number in turbulence. In the review article on BDR performance by Ceccio (2010), disparate results in different facilities were reported despite a similar volume fraction of bubbles treated. Murai (2014) stressed the complexity of parametric combinations in two-phase turbulence while classifying the drag reduction mechanism over a two-dimensional parameter space given by bubble size and main flow velocity.BDR uses dispersed bubbles in the turbulent boundarylayer, and therefore bubbles are distributed more or lessnon-uniformly along the wall. Even with a random distri- bution, bubbles inevitably exhibit local non-uniformity. In addition to the randomness, the Lagrangian motions of individual bubbles accompanying the slip velocity in the liquid phase will amplify this non-uniformity. For vertical pipe geometry, there are several reports dealing with bubble non-uniformity that enhances the flow transition (Lisseter and Fowler 1992; Lammers and Biesheuvel 1996). Bubble clusters can be generated as the slip velocity brings bub- bles closer together. This clustering was clearly observed in free-rising bubble flow by Kitagawa et al. (2004) and Mercado et al. (2010), in vertical bubbly channel flow by Takagi and Matsumoto (2011), and also in horizontal bub- bly channel flow by our group (Murai 2014).Viewed in the Eulerian frame, these are identified as void fluctuations in time and space. Here two questions arise when considering the relationship between void fluc- tuation and drag reduction performance. One is how such fluctuations are amplified in the main direction of flow as it interacts with the turbulent boundary layer. The other is how fluctuations affect the average drag reduction rate. Oishi et al. (2009) found in their channel flow experiments that a local void fluctuation, which was generated naturally within a boundary layer, contributes positively to a time- averaged drag reduction. They also detected a significant phase delay in the local void fraction with respect to the local drag reduction. Pursuing the mechanism hidden in this fact, Park et al. (2015a) generated artificial void waves by repetitive bubble injections at the upstream location of a fully developed horizontal channel flow. They confirmed that BDR performance was improved significantly by gen- erating these artificial void waves.The above-mentioned in-house experiments also alert marine vessel designers that they need to be careful of void waves, which may occur around ship hulls. There has been no attempt yet to quantitatively visualize the void fluctua- tion in an actual application of BDR to ships. The vast dif- ference between model ships of in-house experiments and marine vessels may be whether the system is closed or opened in terms of the two-phase turbulent boundary layer. While channel flow experiments inquire bubble–turbulence interaction in fully developed turbulence between two par- allel walls, marine vessels apply bubbles in spatially devel- oping boundary layers open to outer potential flows. The main question we try to solve experimentally is whether void waves emerge stronger in an open system than in bounded shear flows. To this end, we have designed a fully transparent model vessel, which is essentially equivalent to an experimental cabin cruising in water that allows various bubble behaviors to be quantified visually. In the following, we describe the design features of the model ship to meet this purpose. We then report on the drag reduction perfor- mance and its relation to measured statistics of the bubble

   ‌Table 1 Details of model ships reported in previous research studying bubbly drag reduction

   Titov (1975)	1	8.46	1.0–1.8
   Tokunaga (1987)	2	3.5	5–9
   Yim and Kim (1996)	3	2.6	1.1–1.3
   Watanabe et al. (1998)	4	20	5.0–7.0
   	5	40	5.0–7.0
   Fukuda et al. (2000)	6	7.267	0.7–1.9
   	7	12.0	1–6
   Hirayama et al. (2003)	8	16	1.0–7.0
   Latorre et al. (2003)	9	2.31	5.7–10.4
   Takahashi et al. (2003)	10	12	5–7
   	11	50	5–7
   Katsui et al. (2003)	12	7	1.0–1.4
   Foeth et al. (2010)	13	11.8	3.1–4.9
   Amromin et al. (2011)	14	4.55	12.3–16.5
   Mäkiharju et al. (2013)	15	12.9	5–7.5
   	16	0.912	1.3–2.5
   Jang et al. (2014)	17	8.3	6.7–8.2

   Researchers and reporting year Label in Fig. 5 Ship length (m) Ship speed (m/s)behavior. Finally, void waves, which stand out strongly in the present flow geometry, are analyzed precisely, sup- ported also by a theoretical discussion based on the simpli- fied two-fluid model equations.

2. ‌Design of experimental model ship

   In the study of BDR performance, details of several model ships have already been reported (see Table 1). All these model ships were designed with a flat bottom to avoid bub- bles escaping from the target wall. Each model ship suc- cessfully obtained a certain level of drag reduction, which depended on a combination of various parameters in opera- tion. Whereas total or pointwise drag reduction was con- firmed in these systems, bubble behavior was not carefully measured. We therefore designed a model ship consist- ing of a transparent flat acrylic plate for visualizing bub- bles traveling in the spatially developing boundary layer. Experiments with the model ship were performed in a tow- ing tank to establish conditions of no hydrostatic pressure gradient in the horizontal plane, i.e., the bottom plane of the ship was adjusted horizontally and towed at a constant advancing speed in the stationary water of a large under- ground pool, 100 m in total length. Details of the model ship and the towing system are explained in Sects. 2.1 and 2.2.

   1. ‌Model ship incorporating measuring devices‌

      To analyze the formation of void waves at various length scales, the behavior of bubbles and their relevance to dragreduction performance were investigated simultaneously. We designed a model ship that allows various measuring devices to be attached including ultrasonic bubble meas- urement systems, wall shear stress sensors, and optical bubble-imaging systems. Figure 1 presents schematic dia- grams of the model ship; Table 2 lists its dimensions and basic specifications. Except for aluminum rims, the model ship is for the most part made of transparent acrylic resin; its overall length (L) is 4000 mm, width (W) 600 mm, and height (H) 500 mm. The x, y, and z coordinates are defined respectively as the streamwise distance from the lead- ing edge, the vertical downward position from the bot- tom plate, and the spanwise position from the central axis of the ship. To avoid influences of bow-generated splash- ing waves to the boundary layer structure, two guide walls protruding 20 mm below the bottom plate are attached to both edges (Fig. 1c). To impose an ideal spatially devel- oping boundary layer on the flat bottom plate, the hull of the ship is completely flat with no ridges and has a smooth plate surface. The leading edge of the bottom plate has a 45° bevel to minimize downstream influences of the front- edge flow separation. For this model ship, we also calcu- lated its center of gravity, meta-center in the body-neutral floating state, and recovery moment vector from various inclined attitudes after the adopted strength of materials were determined in high-speed towing conditions. We omit the description of these design features to focus on the two- phase fluid dynamics of the ship.The air injector supplying the bubbles consists of a com-pressor, an airflow rate control system, a buffer chamber of total volume 5.0 × 10−3 m3, and a multi-hole bubble injec- tion plate having 42 open holes, each of 5 mm diameteraLength (L) 4000 mm Height (H) 500 mm Width (W) 600 mmbAir compressorMovable cameraxxyForwardMid-shipzyx = 0.7 mx = 1.1 mAftCamera 30°c65 mmdx = 2.3 mL = 4 mex = 3.5 mProwBottom plateAirRefractive indexShear stress sensor135°matching materialUltrasonic transduceryxBufferchamberFlowSide wallyxPlatewith holes42 × Φ5 mmFlowxy FlowInjected bubblesθ = 8°25 mmh10 mm20 mm90 mmH200 mm

      ‌Fig. 1 Schematic diagram of the experimental facility; a top view of the model ship, b side view of the model ship, c details of the leading edge,

      d schematic of the air injector and e combined mounted system of the shear stress and ultrasonic transducers embedded on the bottom plate

      ‌Table 2 Dimensions and weight of the model ship

      Material Transparent acrylic resinLength (L) 4000 mmHeight (H) 500 mmWidth (W) 600 mmWeight (unloaded) 149 kgDraft (unloaded) 68 mm(see Fig. 1d). The bubble injection plate is located 0.7 m downstream from the leading edge. A servo valve in the air- flow rate control system is operated and managed automati- cally by a PC to supply air with a constant stable volume flow rate (see Fig. 2). The control system for the airflow rate is designed similar to a device developed by Tokyo Gas Corporation (Takeuchi and Kagawa 2013) which is able to supply air at constant volume flow rate up to 2.5 × 10−3m3/s. Pressure logs in the buffer chamber for air injection are shown in Fig. 3, where the unit is in gauge pressure, i.e., pressure differential increase from atmospheric pressure. In our system, the flow rate of air supply is kept mostly con- stant in time and standard deviation of the temporal fluc- tuation relative to mean value was smaller than 1%. The power spectra of the corresponding condition are shown at right panels. The spectra do not include strong peaks (note that the scale of the ordinate is three digits different from the left panels) and can be regarded as white-noise pattern. From this, it can be denied that the void waves observed downstream come straightforward from the initial fluctua- tion in air injection flow rates.To record bubble motion, two cameras were installed on the model ship (Fig. 1b). One was a high-speed video camera (FASTCAM Mini UX 100, Photron, Ltd.) set above the hull. The camera was mounted on two parallel rails that enabled the photographing location to change in the mainCompressorTowing trainTriggerPC & Dater logger(50 MHz)Data logger(50 Hz)Flow Amp (20 Hz)controllerInjectorUltrasonic velocity profiler(Pulse generator)Isothermal chamberThermometer Pressure sensor Flow meter Servo valveShear Ultrasonictransducer transducer

      ‌Fig. 2 Schematic diagram of the flow control system syn- chronized and integrated with multiple measuring devices

      Pressure [kPa]a3.02.01.000 1.00 2.00 3.00 4.00 5.00 6.00 7.00Time (t) [s]Pressure [kPa]b3.02.01.000 1.00 2.00 3.00 4.00 5.00 6.00 7.00Time (t) [s]Pressure [kPa]c3.02.01.000 1.00 2.00 3.00 4.00 5.00 6.00 7.00Time (t) [s]20Intensity [Pa]1000 10.0 20.0 30.0 40.0 50.0Frequency [Hz]Intensity [Pa]201000 10.0 20.0 30.0 40.0 50.0Frequency [Hz]Intensity [Pa]201000 10.0 20.0 30.0 40.0 50.0Frequency [Hz]‌Fig. 3 Time series of pressure in buffer chamber (left side) and linear spectra (right side) at Qg = 1.67 × 10−3 m3/s; a Umain = 2.00 m/s, b2.50 m/s and c 3.00 m/sdirection of flow. The camera recorded a local top view image of the bubbles through the transparent bottom of the ship, from which size, shape, and velocity of individual bubbles were obtained. The other was a waterproof cam- era (HERO3, GoPro, Inc.) set near the rear edge of the ship bottom plate under the water line. This camera recorded all of the bubbles flowing beneath the bottom plate from an oblique direction. In imaging the bubble, three light sources were used and arranged after several trial-and-error attempts to optimize the field of illumination. Basically, the top viewing of bubbles relies on underwater lighting of six large white screens laid on the bottom of the water reser- voir. Each screen is 3 m in length and 8 m in width and is placed at 4-m intervals. Because the white screens reflect diffuse light upward, the images captured the bubbles as backlit shadows. However, this is insufficient for identify- ing film-state bubbles, and we implemented underwater lateral lighting as well. Furthermore, a metal halide lamp supplemented lighting inside the model ship to detect bub- bles smaller than 0.5 mm. The camera then captured strong light-scattering points of all spherical bubbles.A combined system of shear stress sensor (S10W-4,SSK Co., Ltd, Tokyo, Japan.) and ultrasonic transducer was installed at three points along the direction of flow, 1.1,2.3 and 3.5 m from the leading edge; these are hereafter referred to as forward, mid-ship, and aft, respectively. The same type of shear transducer was used in previous studies to monitor the local wall shear stress in bubbly two-phase flows (Kodama et al. 2000; Takahashi et al. 2003; Parket al. 2015a). The frequency response characteristics of the shear transducer were reported by Murai et al. (2007); the temporal resolution is about 20 Hz. Signals from the shear transducers, airflow rate control system, and the velocity of the towing train are recorded simultaneously into a data logger (NR-500, KEYENCE Co.) as shown in Fig. 2. This signal integration allows us to monitor the drag reduction performance in real time during the ship-towing opera- tion. To determine the vertical interfacial position of bub- bles away from the bottom plane, ultrasonic bubble ech- ography developed within our group was applied, details of which were reported by Park et al. (2015b). The ultra- sonic transducer was embedded inside the bottom wall at a tilt angle of θ = 8° to the vertical direction. The trans- ducer is located at approximately 28 mm away from the shear transducer, both of which are collocated in a disk- hold combined mount system (Fig. 1e). Furthermore, an ultrasonic velocity profiler (UVP-DUO MX, MET-FLOW S.A., Lausanne, Switzerland) was adopted not only for use in velocity profile measurements of the liquid phase but also as the ultrasonic pulse generator in the bubble ech- ography. Table 3 summarizes the parameter settings for all instrumentation.In this paper, we report on drag reduction performance of the present model ship and focus in particular on the bub- ble behavior that is observed using the two bubble-imaging techniques recorded by the top-viewing high-speed cam- era and the underwater camera. Liquid velocity and bubble echography information will be reported separately.

      ‌Table 3 Parameters of the measuring instruments

      Shear stress sensorMeasurement area 25π mm2Temporal resolution 20 HzRange of shear stress ±250 PaUltrasonic pulse generator and ultrasonic transducerUltrasonic basic frequency 4 MHzNumber of cycle 4Pulse repetition frequency 3.2 kHzVoltage for ultrasonic emission 150 VUltrasonic beam diameter 5 mmDivergence half-angle 2.2°Data logger for echographySampling frequency 50 MHzRange of voltage ±2 VResolution of voltage 4 mVCamera above the shipFrame rate 500 fpsResolution 0.08 mm/pixelUnderwater cameraFrame rate 120 fps

   2. ‌Towing test facility‌

      The present series of experiments was performed at a large towing test facility in Hiroshima University, Japan. The facility is certified as a standards institute for ship perfor- mance examination, which was registered by ISO 9001

      ‌Fig. 4 Photographs of the ship-towing test facility in Hiroshima Uni- versity; a towing water reservoir, b towing train and c model ship hooked up to the train

      distance in constant-speed testing. Under the maximum pos- sible speed, the duration for steady-state operation is 7 s or 21 m in distance.

   3. ‌Experimental conditions

      The two controllable parameters of the present experiment are towing speed (Umain) and air volume flow rate (Qg). Umain ranges from 2.00 to 3.00 m/s, at which the superficial air layer thickness (δg) defined by(International Standard Organization) in 2012. The towing test facility comprises a 100-m-long rectangular water res-δ =   Qg  .g UmainW(1)ervoir and a powerful towing train (Fig. 4a, b).The train is a large steel carriage that runs on two par- allel rails either side of reservoir and running along its full length. The motion of the train is managed by the operator at the cockpit above the train. Speeds of the train are con- trolled. Table 4 summarizes dimensions, performance, and conditions of the present towing test facility. Two rigid pil- lars inside the train support our model ship (Fig. 4c). The support system does not allow the model ship to incline because of hydrodynamic moments, i.e., pitching, yawing, and rolling do not occur in towing operations. The bottomis maintained below 2 mm; here W denotes the spanwise width of the ship, which is 600 mm in the present model ship. This air layer thickness is often used in experimen- tal studies for BDR because it roughly estimates the dis- placement thickness of the liquid velocity boundary layer agitated by bubbles. Figure 5 shows a comparison of pre- sent experimental conditions with those of past studies (see Table 1) mapped on the parameter space of Umain and L. Onthe map, the oblique lines indicate the contours of Froudenumber defined asUmainplane of the ship was fixed horizontally at 90 mm below the water surface. This depth is close to the natural draft of theFr =√gL .(2)ship determined by the balance between weight and buoy- ancy. To mitigate hydrostatic pressure gradient along at the bottom plane, we adjusted the water depth precisely to beThe horizontal line in the map indicates the unit contour of the cavitation number defined by2(p − pv)the same along the entire ship bottom. Because of safety issue associated with limits in train deceleration, the modelCa =2 .ρUmain(3)ship travels approximately 80 m. Within this distance, train acceleration and deceleration further reduce the availableHere, g, p and pv are the acceleration due to gravity, the local static pressure at the ship bottom, and the vapor

      ‌Table 4 Dimensions, performance and conditions of the towing test facility

      ‌Table 5 Experimental conditions of bubble injection into spatially developing boundary layer

      Water reservoir			Controlled parameters
      Length in towing direction	100 m		Towing speed	2.00–3.00 m/s
      Depth Width	3.5 m 8.0 m		Air volume flow rate (Qg) At the bubble injector, x = 0.7 m	0.42 × 10−3–2.50 × 10−3 m3/s
      Content in the tank Temperature of water	Water (clean city water) 29 °C		Reynolds number (Rex) Thickness of the boundary layer (δ)	1.6 × 106–2.5 × 106 13.6–14.8 mm
      Density of water (ρ)	996 kg/m3		Void fraction inside the boundary	1.9–16.1%

      Kinematic viscosity of water (ν) 0.847 × 10−6 m2/sSurface tension of water 71.2 × 10−3 N/mSpeed of sound in water 1507 m/sTowing trainlayer (αδ)At the forward, x = 1.1 mReynolds number (Rex) 2.6 × 106–3.9 × 106 Thickness of the boundary layer (δ) 19.6–21.2 mmLength 8 mWeight 20 × 103 kgVoid fraction inside the boundary layer (αδ)1.4–11.2%Total power of the towing train by four electric cells65 kWAt the mid-ship, x = 2.3 mReynolds number (Rex) 5.4 × 106–8.1 × 106Maximum speed 3.00 m/sMinimum speed 0.10 m/sVariable range of acceleration 0.1–0.8 m/s2 Period of constant towing speed at maxi- 7.0 smum speedThickness of the boundary layer (δ) 35.3–38.3 mm Void fraction inside the boundary 0.7–6.2%layer (αδ)At the aft, x = 3.5 mReynolds number (Rex) 8.2 × 106–1.2 × 107 Thickness of the boundary layer (δ) 49.4–53.6 mmVoid fraction inside the boundary layer (αδ)14Ca = 1.0921710 15137 845 1116312 16200.5–4.4%Ship speed (Umain) [m/s]10 and also the side walls separate the BDR-target bottom plane from the outer flow.On the target plane, the Reynolds number defined byRex= xUmain ,ν(4)10.50.5 110 60Length of ship (L) [m]‌where x and ν are the streamwise distance from the front edge of the bottom plate and kinematic viscosity of water, respectively, describes the dynamic similarity of the spa- tially developing boundary layer for single-phase flow. Generally, a turbulent boundary layer forms on a flat plate at Rex > 5 × 105. Thus, for the present model ship, a turbu- lent state is reached 0.25 m from the front edge, which is forward of the point of bubble injection at 0.70 m from the front edge.

      Fig. 5 Ellipse plots of the experimental parameter conditions exam-

      ined in the past relative to the present experimental condition (closed ellipse); numbers labeling open ellipses correspond to the experimen- tal ships listed in Table 1pressure of the water, respectively. For reference, a line is drawn assuming a water depth of 100 mm. In actual situa-The air volume flow rate (Qg) is regulated at fixed values taken from the range 0.42 × 10−3 to 2.50 × 10−3 m3/s. The void fraction inside the turbulent boundary layer is esti- mated by0 W∫ δ yα =     Qg    ≈     Qg     ≈         Qg        ,tions, large vessels such as container carriers and oil tank- ers sail at Fr ≤ 0.3 so as not to intensify their wave makingδ Ql + QgW ∫ δ uydy10 Umain δ 7 dy‌(5)resistance. Although Fr in our experiments is larger than 0.3, the wave does not affect the BDR phenomena at the bottom plane because the model ship is fixed in altitude,where we assume the initial liquid phase velocity distribu- tion (uy) is given to estimate the liquid volume flow rate inside the turbulent boundary layer. In single-phase flow,the boundary layer thickness (δ) spatially increases down- stream as estimated by   ν   1aShear stress (τw0) [Pa]15.0‌10.0ForwardMid-shipAft4δ = 0.37x 55.Umain(6)5.0Equations (5) and (6) imply that the mean void fraction inside the boundary layer, αδ, decreases downstream when Umain and Qg are fixed. However, when drag reduction occurs, the local instantaneous values of Qg and uy couplein the space–time domain along the wall. This aspect is ourfocus of attention in the rest of the paper. All other details of the experimental conditions are listed in Table 5.

3. ‌Performance of the model ship

   0.0Friction coefficient (Cf)b 10–20.00 1.00 2.00 3.00Towing speed (U) [m/s]ForwardMid-shipAftIn this section, we present the basic drag reduction per- formance obtained with the present model ship and the detailed conditions of bubbles streaming beneath the bot- tom plane.

   1. ‌Wall shear stress

      Figure 6a shows the relationship between the wall shear stress with no bubble injection (τw0) measured at the three locations and several Umain. Wall shear stress increases with ship speed more than linearly but less than quadratic. This trend agrees in general with a previous study using a tow- ing flat plate (e.g., Mori et al. 2009). Of the three locations, the forward location has a wall shear stress higher than the other two aft locations. This is explained by the expanding thickness of the boundary layer along the main direction of flow. Figure 6b shows the same trend plotted in the two- dimensionless parameter space, i.e., Rex and coefficient of friction which is defined by2τw0  10–3Turbulent flowTransitionLaminar flow105 106 107Reynolds number (Rex)‌Fig. 6 Frictional drag of the model ship in single-phase flow condi- tions, where error bars indicate standard deviations of the temporal fluctuation; a local wall shear stress as a function of ship speed, and b coefficient of friction as a function of Reynolds number; the two curves represent the Blasius friction law for laminar flow and the empirical coefficient of friction for turbulent flow (Schlichting 1979)In this measurement, nominal accuracy of the shear transducer is 0.05 Pa at a frequency lower than 20 Hz (tem- poral response limit). The error bars in the graph indicate the standard deviation of temporally fluctuating shear stress subject to wall turbulence. It is noted that the deviation for the forward position were obtained relatively high. This occurs because the boundary layer induced from the front edge of the flat bottom intermittently keeps laminar state due to smooth wall surface. In case of bubbling, randomlyρUCf = 2 .main(7)distributed bubbles play a role as trigger for more stable transition to turbulent boundary layer in the beginning partThe two curves refer to the Blasius friction law for lami- nar boundary layer flows and the empirical coefficient of friction for the turbulent state by Schlichting (1979). The corresponding formulae are given byof the boundary layer.Both the wall shear stress during bubble injection (τw) and that without bubbles (τw0) were measured and time- averaged over the 7 s at Umain = 3.00 m/s. From theirxCf = 1.328Re−1/2(8)ratio (Fig. 7), friction in the forward section intensifies with increasing αδ except for dilute fractions αδ < 2.0%.for laminar flow, and−0.558C−1/2 −1We deduce from this trend that bubbles at high flow ratesnear the injector suddenly perturb the boundary layer and enhance momentum transfer, resulting in an increase inCfef = Rex(9)drag. In contrast, drag reduction is maintained under dilutefor turbulent flow. The measured data points support the notion that the present flat bottom ship properly forms a spatially developing turbulent boundary layer in accord- ance with Schlichting’s formula.bubble injection. At mid-ship and aft, the wall shear stress produces a well-known decline with increasing void frac- tion. The maximum drag reduction rate that we confirmed within the tested range is 30% at void fraction αδ = 4.2%.aRatio of shear stress τ w/τw0)2.0)1.81.61.41.21.00.80.0 2.04.06.0 8.0Forwardbubble-mixed fluid density. We also found that the drag reduction curves have similar trends both at mid-ship and aft when parameterized by the boundary layer void frac- tion, αδ. This suggests that drag reduction occurs with a quasi-steady mechanism between mid-ship and aft. Hence, we can analyze the behavior of bubbles within this range to elucidate the bubble-to-drag influence.

   2. ‌Bubble distribution

      Ratio of shear stress τ w/τw0)b 1.1)1.00.90.80.70.6Void fraction (αδ) [%]0.0 2.0 4.0 6.0 8.0Mid-ship AftA question that has been a long-debated issue in the field of BDR research is: What kind of bubble contributes actively to drag reduction? With our fully transparent ship model, an answer does appear in a series of bubble snapshots (Fig. 8). The snapshots are sampled from images taken with the high-speed video camera recording top views of bubbles near the wall. The camera is located mid-ship (x ≈ 2.1 m) on the centerline close to the shear stress measurement point. From the images, the bubble sizes range widely from 1 to 50 mm. With larger Umain, sizes are seemingly smaller. This general trend is understandable given the bubble-Void fraction (αδ ) [%]‌Fig. 7 Wall shear stress modified by bubble injection with respect to void fraction at a ship speed of Umain = 3.00 m/s; a forward, and b mid-ship and aft, where error bars indicate standard deviationsFrom the experimental plots, the mean impact factor of the void fraction to the drag reduction rate, defined by (1 − τw/τw0)/αδ ranges from 4 to 7. An impact factor larger than unity proves that the present drag reduction is enhanced by the two-phase mutual interaction inside the boundary layer and not just because of a bulk decrease inshearing force, which is characterized by Weber number. For αδ > 2%, bubble size exceeds the mean bubble–bubble interval distance observed in the two-dimensionally pro- jected top views. For αδ > 3%, bubbles occupy more than a half of the imaged area in these top views, and are likely to be an air barrier against high-speed outer flow.To obtain the bubble size distribution, we produced a digital image from the high-speed video images. Figure 9 shows a video image of a single bubble taken under condi- tions Qg ≈ 0.42 × 10−3 m3/s at x ≈ 3.3 m. The bubble sur- face is measured by image binarization using an adequate threshold value close to the background brightness. The

      ‌Fig. 8 Snapshots of stream- ing bubbles near mid-ship,

      x ≈ 2.1 m; αδ below each panel z2.00is the boundary layer void x0.42㽢10−30.83㽢10−3Air flowrate (Qg) [m3/s]abcdefαδ = 1.1%αδ = 2.2%αδ = 3.3%αδ = 4.5%g h i jαδ = 5.6%kαδ = 6.7%lαδ = 0.9%αδ = 1.9%αδ = 2.8%αδ = 3.7%n o pαδ = 4.7%qαδ = 5.6%mr1.25㽢10−3 1.67㽢10−32.08㽢10−32.50㽢10−3fraction20 mmMain flow velocity (Umain) [m/s]2.503.00αδ = 0.8%αδ = 1.6%αδ = 2.4%αδ = 3.2%αδ = 4.0%αδ = 4.8%(de) [mm]θb [rad]a b Equivalent diameterCenter of massxz‌Fig. 9 Definition of equivalent diameter (d ) and phase (θ ) of a bub-(2010). This value suggests qualitatively that the bubbly two-phase layer slides downstream creating a difference in velocity between the solid wall and the outer flow. How- ever, an exact explanation is at present impossible to val- idate. From a careful look at the graphs, we find that the velocity ratio ranges between 0.45 < ub/Umain < 0.65 and rises with large Umain and high Qg towing operations. Such conditions correspond to conditions for high drag reduc- tion (see Fig. 7b). In the downstream region where the drag reduction rate relaxes, the bubble velocity ratio alsoble; ae b b relaxes. What implications does this have for a mechanism?raw image sampled from the high-speed video camera, and binary image: gray dashed line marks the average radius of the bub- ble surfacecoordinates of the surface are used to find the two-dimen- sional center of gravity of the bubble to compute the bub- ble velocity. The bubble size is measured from the angular average of the distance from the center to the surface coor- dinates; doubling the value gives the circle-equivalent bub- ble diameter (de). Figure 10 shows the probability distribu- tion of de as it changes with ship speed. Note that bubbles smaller than 1 mm and larger than 8 mm are considered out of range, because the quality of back-lit imaging of their surfaces is limited.The measured results indicate that the peak bubble diam- eter ranged from 3 to 4 mm as Umain increases from 2 to 3 m/s. This relationship is in very good agreement with the formula which Hinze (1955) proposed and Sanders et al. (2006) confirmed by experiment. Their formula describes a survivable bubble size against surrounding turbulent shear. Our present data involve both smaller and larger ones due to active fragmentation and coalescence among bubbles accumulated close to the wall. It is understandable that the theoretical survivable bubble size takes the peak population since fragmentation and coalescence rate are reversed at this criterion.

   3. ‌Bubble velocity

      Predicting theoretically the mean bubble velocity is very difficult as the bubbles are suspended inside the turbulent boundary layer during drag reduction. Bubble velocities depend sensitively on the normal distance from wall, i.e., where the bubbles accumulate in the boundary layer in which 99% of the spread in velocity is localized.Figure 11 shows the averaged advection velocity (ub) of the bubbles measured by particle image velocimetryapplied for bubbles imaged in Fig. 8. The error bar around each data point indicates the standard deviation from 300 sampled images. The values of ub are roughly half of Umian for all experimental conditions; this trend is basically con- sistent with a previous model ship study by Johansen et al.We have found a certain relevance to void waves that we disclose in the next section.

4. ‌Void wave measurement

   During towing operations of the present model ship, we observed visually the propagation of waves constituted by the sparseness and denseness of the local bubbles on the bottom plane. In this section, we analyze the wave quantita- tively. Such void waves were observed to form in spatially developing bubbly flows inside vertical pipes (Lisseter and Fowler 1992; Lammers and Biesheuvel 1996); being similar, the phenomenon we observed is termed likewise. Because these void waves can have broad spectra, we ana- lyze the wave on two different scales as reported separately below.

   1. ‌Behavior of void waves

      To analyze the void waves quantitatively, we used the video images taken by the underwater camera located at the stern (x ≈ 3.7 m). Figure 12 shows instantaneous images of the bubbles beneath the bottom plane of the model ship, in which bubbles are flowing upwards. A sampling area that is marked by the thin white rectangle in the figure was chosen for a line scan of the bubbles integrated over time. Figure 13 depicts the timeline scanning image under condi- tions Umain = 3.00 m/s for three different Qg values. Many lateral waves comprising bubble clusters can be seen in these timeline scanning images.To determine the peak frequency of the void wave, the scanning images were analyzed using Fourier analysis after a background subtraction was applied. The spectra obtained are presented in Fig. 14, where the ordinate of each panel represents the bubble brightness level normalized by the maximum bubble brightness. Note that brightness is not a physical quantity, and therefore each value gives only the magnitude of the void fluctuation. The abscissa rep- resents the frequency (fvoid) directly obtained by frame interval time; on top is the corresponding wavelength scale, estimated using ub/fvoid. All the spectra are subjecta 0.40PDF [(mm)-1]0.300.200.10b 0.40PDF [(mm)-1]0.300.200.10c 0.40PDF [(mm)-1]0.300.200.100.000.02.0 4.06.08.00.000.02.0 4.0 6.08.00.000.02.0 4.0 6.08.0Equivalent diameter (de) [mm]Equivalent diameter (de) [mm]Equivalent diameter (de) [mm]‌Fig. 10 Probability density distributions of the equivalent bubble diameter evaluated for Qg ≈ 0.42 × 10−3 m3/s at x ≈ 3.3 m, where small bub- bles, de < 1.0 mm, are regarded as out of range; a Umain = 2.00 m/s, b 2.50 m/s and c 3.00 m/saAveraged advection velocity of bubbles (ub) [m/s]3.002.001.000.00 0.500.00 1.00 1.50 2.00 2.50Air flowrate (Qg) [m3/s]bAveraged advection velocity of bubbles (ub) [m/s]3.002.001.00Umain2.00 m/s2.50 m/s3.00 m/sUmain2.00 m/s2.50 m/s3.00 m/sto broad background power because of the discrete nature of the dispersed phase. Against the background, clear peak frequencies can be seen in the range 3 < fvoid < 8 Hz for Qg = 1.67 × 10−3 and 2.50 × 10−3 m3/s. Within this range, there are two general trends seen in comparing the nine spectra. One is that the peak frequency rises with increas- ing ship speed Umain, particularly at high QG. Its rise is roughly linear with Umain inferring that the wavelength of the void wave is independent of ship speed. The other is that the peaks shift to lower frequencies with increasing Qg. This trend may be attributed to the coalescence of two neighboring void waves or collapse of a single void wave in consequence of an upper limit in the local bubble num- ber density sustainable inside a single wave. The upper limit and capacity are unclear quantitatively in this spectral analysis, and therefore we proceed to an analysis of bubble clustering.

   2. ‌Local bubble arrangement

      0.000.00 0.501.00 1.50 2.00 2.50Air flowrate (Qg) [m3/s]

      ‌Fig. 11 Averaged advection velocity of bubbles measured at a

      x ≈ 2.1 m and b 3.3 m, where error bars indicate standard deviationsRear regionzy‌Fig. 12 Snapshot taken by the underwater camera; the framed area is located at x ≈ 3.7 m and used to generate timeline scanning imagesTo find out how the void waves are initiated and amplified at the bottom plane, we here analyze local bubble images to quantify a network-like organization of bubble distribu- tion. Even if a peak is not observed in the void frequency spectra, many bubble clusters clearly emerge as evi- denced in Figs. 13a and 14h. We believe that these meso- scale structures in the bubble distribution are perturbative source triggering the void waves. To begin, accumulations of bubbles are evaluated using distances between the cent- ers of two bubbles (dcb). The results for three different shipspeeds are presented in Fig. 15, which shows probabilitydensity function (PDF) of the mutual distance. A vertical line in each panel indicates the reference distance when all the bubbles are distributed uniformly. The abscissa shows dimensional and non-dimensional scales based on the mean bubble diameter. In all three cases, the PDF has a profile weighted on the left-hand side of the reference value. For a perfectly random distribution, the reference value shifts to‌Fig. 13 Timeline scanning images taken at x ≈ 3.7 m for Umain = 3.00 m/s; aQg ≈ 0.83 × 10−3 m3/s,b 1.67 × 10−3 m3/s and c2.50 × 10−3 m3/sa 0.3Width (z/W)0.20.10.0-0.1-0.2-0.30.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00Time (t) [s]b 0.3Width (z/W)0.20.10.0-0.1-0.2-0.30.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00Time (t) [s]c 0.3Width (z/W)0.20.10.0-0.1-0.2-0.30.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00Time (t) [s]1/√2 times the original. Nonetheless, the peak of the PDF is located to the left from the random state. This proves that bubble clusters actively form more often than that occur- ring naturally in a random state. Moreover, we can confirm that the peak of the PDF shifts left, i.e., toward narrower bubble spacing, as the ship speed increases. This suggests a positive correlation between drag reduction and the void wave inside which the bubble spacing narrows. In addition, it is worth noting that the spacing can be narrower than the mean bubble diameter exceeding the contact limit of spher- ical bubbles. This occurs mainly due to the large deviation in bubble size and may play as a role of trigger for down- stream void wave generation.We also analyzed the bubble probability distribution in the azimuthal directions. The aim of the analysis is to find a source of void wave generation from microscopic point of view. From the distribution, we subtract perfect uni- form distribution in order to emphasize its anisotropy of the local bubble arrangement patterns. The result is shown in Fig. 16. The same approach was employed by Kitagawaclustering of bubbles because of significant bubble–bub- ble interaction in high viscous oil. Surprisingly, the pre- sent results also show a certain statistical heterogeneity although the target is distinct from a laminar boundary layer. At Umain = 2.00 m/s, bubbles are arranged mostly in the isotropic state (see Fig. 16a). It becomes heterogene- ous at Umain = 2.50 m/s and Umain = 3.00 m/s; the bub- bles exhibit a high probability density at angles θb = 0 and θb = π. This implies that bubbles aggregate mainly in the streamwise direction. This is consistent with the observa- tion (Fig. 8) that bubbles form many chains in the main direction of flow. Nierhaus et al. (2007) and Harleman et al. (2011) have reported a preferential concentration of bubbles in turbulent boundary layers, although the bubbles were smaller than coherent structures inside the boundary layer. Smith and Metzler (1983) and Zacksenhouse et al. (2001) found that the spanwise width of low-speed streaks in a single-phase turbulent boundary layer is around 100 times the friction length defined byw0νet al. (2004) for wall-sliding bubbles along a vertical platein a stationary liquid. They found a clear heterogeneouslτ = √τ /ρ .(10)Qg [m3/s]0.83㽢10−31.67㽢10−32.50㽢10−3Umain [m/s]a Wavelength (ub/fvoid) [m] bWavelength (ub/fvoid) [m] cWavelength (ub/fvoid) [m]2.00Arbitrary linear scale00 5 10 15 20Frequency (fvoid) [Hz]00 5 10 15 20Arbitrary linear scaleFrequency (fvoid) [Hz]00 5 10 15 20Arbitrary linear scaleFrequency (fvoid) [Hz]d Wavelength (ub/fvoid) [m] eWavelength (ub/fvoid) [m] fWavelength (ub/fvoid) [m]2.50Arbitrary linear scale00 5 10 15 20Frequency (fvoid) [Hz]00 5 10 15 20Arbitrary linear scaleFrequency (fvoid) [Hz]00 5 10 15 20Arbitrary linear scaleFrequency (fvoid) [Hz]g Wavelength (ub/fvoid) [m] hWavelength (ub/fvoid) [m] iWavelength (ub/fvoid) [m]3.00Arbitrary linear scale00 5 10 15 20Frequency (fvoid) [Hz]00 5 10 15 20Arbitrary linear scaleFrequency (fvoid) [Hz]00 5 10 15 20Arbitrary linear scaleFrequency (fvoid) [Hz]

      ‌Fig. 14 Linear spectra of the void wave naturally generated under- neath the flat bottom of a model ship; the spectra are obtained from timeline scanning images (e.g., Fig. 13), and given as normalized val-

      ues using the maximum brightness of the images with wavelengths calculated using ub at x ≈ 3.3 m (see Fig. 11b)a Normalized distance (dcb/de)
      1. Normalized distance (dcb/de)
      2. Normalized distance (dcb/de)
      0.00.40PDF [(mm)-1]0.300.200.101.02.03.00.00.40PDF [(mm)-1]0.300.200.101.0 2.03.00.00.40PDF [(mm)-1]0.300.200.101.0 2.03.00.000.05.010.015.00.000.05.0 10.015.00.000.05.0 10.015.0Distance (dcb) [mm]Distance (dcb) [mm]Distance (dcb) [mm]‌Fig. 15 Probability density distribution of the mutual distance between the closest pair of bubbles for Qg ≈ 0.42 × 10−3 m3/s at x ≈ 3.3 m: a Umain = 2.00 m/s, b 2.50 m/s and c 3.00 m/s. Small bub-bles, de < 1.0 mm, are regarded outside the analysis; the vertical solid and dashed lines indicate the distances between bubbles for perfect uniform distribution and random distributioncπ/2πaπ/2πbπ/2πzx0 0 03π/23π/23π/2‌Fig. 16 Deviation of probability density distribution in the azimuthal direction of the closest bubble obtained for Qg ≈ 0.42 × 10−3 m3/s at x ≈ 3.3 m from uniform distribution: a Umain = 2.00 m/s, b 2.50 m/s,and c 3.00 m/s. Small bubbles, de < 1.0 mm, are regarded as outside the analysis; gray circles indicate values of PDF, and positive values mean higher bubble density than the uniform situationThe width has a comparable scale with that of the streamwise vortices because of coupling in a process termed a self-sustaining cycle (Hamilton et al. 1995). Therefore, 100 lτ can be regarded as a representative length scale of the coherent structure compared with bub- ble size. In our experimental condition, the length scale is estimated to be approximately 0.7–1.0 mm for Umain in the range 2.00–3.00 m/s. Hence, the present bubble size is sev- eral times larger than the coherent structure. With such a condition, the bubble motion may be subject to a stochas- tic behavior that weakens the spatial structure of the bub- ble distribution. Nevertheless, our experimental visualiza-upon high velocity gradient (e.g., Oishi and Murai (2014)). Furthermore, 100% of bubbles exist inside the turbulent boundary layer, and thereby it is important to consider the bubble motion in relevance to the wall shear stress rather in unbounded space.

      1. ‌Mathematical derivation

         The conservation of momentum for a two-phase mixture (e.g., Murai and Matsumoto 2000) in the main direction of flow is described bytion showed a clear formation of many chained bubbles. We believe a synergy arises between the local high void∂fuρ∂t∂fu2+∂x∂fuv+∂y∂fuw+∂z‌(11)fraction and local drag reduction, which alters the original∂p ∂2u∂2u∂2ucoherent structure inside the void wave. As this is not com- pletely explained in the present study, we suggest this as an=− + µ∂x∂x2 +∂y2 +,∂z2open problem for future research, which would be solved using artificial void waves.

5. ‌Theoretical description of void wave

   To support theoretically the wave-outstanding phenome- non in the bubbly two-phase boundary layer, we attempted to derive a wave equation for the void fraction from con-where f, p, ρ, and μ are volume fraction, pressure, density, and viscosity of liquid phase, respectively. Here f is the same as 1 − fG, where fG is the local instantaneous void fraction. The velocity components, u, v, and w are those of the liq-uid phase in the streamwise (x), wall-perpendicular (y), and spanwise (z) directions, respectively. A spatial integration of Eq. (11) with respect to the boundary layer thickness (δ) givesservation laws. Since bubbles migrate in the horizontal direction without certain base slip velocity to liquid phase, drift flux model approach cannot be employed. Thus, thepresent observation of void wave is essentially different∂fu dy +δ∫∫∂t0 0∫δ2δδ∂fu ∫dy +∂x0∫δ∂fuv dy∂y‌(12)from kinematic void wave that is a wave phenomenon relying on void-to-slip correlation as reported by Pauchon and Banerjee (1988), and Lahey (1991). Even bubble’s1 ∂p=−ρ ∂x0dy + µρ0∂u2∂y2 dy.slip velocity is hardly defined since the relative velocity of the liquid phase to the bubble interface takes opposite sign between top and bottom surface of a single bubbleHere the terms including w and z in Eq. (11) disappear through the averaging inside the boundary layer. The third term in L.H.S. of Eq. (12) is estimated asδ∫ ∂fuv∫ ∂f (u − U    )∂(1 − α )(U − U    )0dy =∂y0fuv δ = fuv|— fuv|0= fv|δ Umain = VUmain,‌(13) 0mainδ∂tdy = δδ δ main∂t,‌(20)δwhere Umain stands for flow speed outside the boundary layer. Here we assume f = 1 at y = δ, i.e., no bubbles exist on theborder of the boundary layer. V denotes the mean downward∫ ∂fu(u − Umain) dy = δ (1 − αδ)Uδ(Uδ − Umain) ,velocity of liquid phase on the border, and can be a function of time in the boundary layer flow subject to unsteadiness. However, here we treat V as constant following to the similar theory for single-phase boundary layer flows.From volume conservation law for liquid phase, another∂xδ0and∂x(21)equation stands for the same boundary layer as‌1 ∫ ∂pµ ∫ ∂u2µ ∂u δ∂f ∂fu+∂fv+∂fw+= 0.(14)—dy +δδρ ∂x ρ0 0∂y2 dy = 0 + ρ∂y 0∂t ∂x ∂y ∂zµ ∂u τwSpatial integration of Eq. (14) inside the boundary layer gives= ρ 0 − ∂y=− .δρ‌(22)δ∫∫δ∂f dy +∂t0 0δ∫∂fu dy +∂x0∂fvdy = 0,∂y‌(15)Equation (15) satisfying Eq. (16) can be also rewrittenaswhere spanwise velocity w disappears due to averaging. The third term of Eq. (15) is estimated by∂(1 − αδ) + ∂(1 − αδ)Uδ = −V .∂t ∂x‌(23)δ Using Eqs. (20)–(23), Eq. (18) can be simplified as∫ ∂fv dy = fv δ = fv|= V .‌(16)∂(1 − α )U  ∂(1 − α )U2 τ  ∂y 0 δ0δ δ +∂tδ δ =− w + VUmain.∂x ρδ‌(24)Hence, Eq. (15) has the following relation. The derivation of Eq. (24) from Eq. (11) is in part simi-δV = − ∫0δ∫∂f dy −∂t0∂fudy∂x‌(17)lar to von Karman’s momentum integral approach for a sin- gle-phase spatially developing boundary layer. In our case, void fraction profile as functions of time and streamwise coordinate, αδ(x, t), is taken into account under the condi-Substituting Eq. (17) to Eq. (13) and then to Eq. (12), weobtaintion of αδ = 0 at y = δ.Equation (23) also produces following equation as being taking derivative in time,δδ∫ ∂f (u − Umain) dy + ∫ ∂fu(u − Umain) dy‌∂2(1 − α )  ∂ ∂(1 − α )U∂t ∂x0 0∫∫δ δ(18)δ +∂t2 ∂xδ δ = 0.∂t‌(25)1 ∂p=−ρ ∂x0dy + µρδ+0∂u2∂y2 dy.As the first term of Eq. (24) is identical to the time dif- ferential in the second term above, Eq. (25) is rewritten asHereafter, we use boundary-layer-averaged quantities,∂2αδ∂ ∂(1 − αδ)U2+ τw= 0.(26)which are defined by∫∫δ δ∂t2 ∂x ∂x ρδα = 1δ δ1(1 − f )dy, Uδ = δudy.(19)This equation expands to giveδ0 0 ∂2αδ+ (1 − αδ)∂2U22 ∂2αδ= Uδ—  ∂ τw .‌(27)Each term in Eq. (18) is therefore linearly approximated ∂t2as∂x2∂x2∂x ρδIn a slowly developing boundary layer, the second term which is rearranged to give2on the left-hand side in Eq. (27) is negligible compared with other terms, and thus we obtain∂2αδ = C2∂2αδ − G1τw0∂αδ , C2 = U2 − Gτw0 .∂2αδ2 ∂2αδ∂ τw∂t2∂x2ρδ ∂xδ ρδ‌(32)‌∂t2 = Uδ∂x2 − ∂x ρδ .(28)Equation (32) takes the form of a wave equation for void fraction. Hence, it supports mathematically that the genera-Noe that τw in Eqs. from (22) to (28) stands for the shear stress acting on fluid at the wall, but not on the wall. In order to avoid confusion due to action–reaction force law,let us rewrite Eq. (28) reversing the sign as∂t2 = Uδ∂x2 + ∂x,(29)∂2αδ 2 ∂2αδ  ∂ τw1 −ρδtion of the void wave is inherent to the system that accom- panies drag reduction. Equation (32) also indicates that the void fraction propagates with speed (C) that decreases with G2. This is consistent with the present experimental obser- vation, i.e., void waves propagate at a speed significantly slower than the ship speed. The value of G2 is estimatedG =from the present experimental result to be approximatelyNow this equation expresses the relationship among local void fraction, local flow speed, and local wall shear stress that acts on the wall.2δ2 CfC 2Uδ≈ 0.1 m0.0031 − (0.6)2≈ 21 m.(33)

   1. ‌Void wave equation

   Our problem is how to model the wall shear stress in Eq. (29). Let us assume a wall shear stress given as an experimentally correlated model,Here 99% thickness of boundary layer is adopted for δ for which Uδ is estimated by (7/8)Umain. The estimated value of G2 means that bubbles have a transient effect on drag reduction within 21 m in the main direction of flow, being much longer than the boundary layer thickness. It alsoτw = τw01 − G1αδ − G2∂αδ ,∂x‌(30)infers that the wavelength of the void wave is shorter than 21 m. From another perspective, void waves do not amplify themselves for wavelengths longer than 21 m in the rangewhere τw0 denotes wall shear stress in the absence of bub- bles, i.e., the original value before bubble injection. The parameter G1 represents an impact factor for drag reduc-tion to void fraction supplied inside the boundary layer. The performance of G1 for a variety of flow conditions was previously reported by Murai (2014). The other parameter G2 in Eq. (30) is a new term that describes the impact of the local gradient of the void fraction. This term is mod- eled based on a previous finding that local wall shear stress has a significant phase shift in the oscillatory fluctuation of the local void fraction (Oishi et al. 2009). The two param- eters, G1 and G2, can take both positive and negative values depending on bubble size and flow speed. We emphasize that there is no physical evidence why these two terms are linearly additive to describe local wall shear stress. That is, Eq. (30) simply expresses a mathematical constitution of the wall shear stress composed of the first two derivative terms with respect to void fraction, i.e., α(0) and α(1). Besidetested because transients no longer are present beyond this length. The second term in Eq. (32) plays the role of a source term in the wave equation. When the value of G1 is high, the wave is initiated.Note that the present attempt at a theoretical descrip- tion only proves the potential existence of void waves in the system, but does not provide a mechanism for genera- tion or a determination of its wavelength. We are presently studying this phenomenon using artificially generated void waves in the upstream location to observe its amplification along the main direction of flow. We shall report on this aspect in a separate contribution.

6. ‌Conclusions

We designed a 4-m-long model ship with a fully transpar- ent flat bottom plate to investigate frictional drag reduction

δ δ

it, the authors are not insisting that the term with G1 would be replaced with the term with G2 in interpretation of bub- ble drag reduction since both terms stand simultaneously and contribute to independent phenomena.

Substituting Eq. (30) into Eq. (29) obtains

following the injection of air bubbles. Various measure- ment instruments were incorporated into the model ship, including a high-speed video camera for imaging bubbles from above, an underwater camera to capture the behavior of bubbles, ultrasonic measurement systems, and wall fric-

tion measurement systems. The model ship was towed in

∂2αδ 2 ∂2αδ

τw0

∂αδ ∂2αδ

a 100-m-long water reservoir, to confirm drag reduction

∂t2 = Uδ

∂x2 + ρδ

−G1 ∂x − G2 ∂x2 ,

(31)

of up to 30% at an average void fraction of 5% inside the

turbulent boundary layer. In this paper, details of the ship design and measurement systems were presented, particu- larly regarding void waves generated entirely along the ship bottom. Through image processing of several sets of bub- ble behavior visualizations, we confirmed that the average bubble streaming velocity is approximately half the ship’s speed. We observed clearly formed void waves having a distinct wavelength irrespective of ship speed under the conditions tested. The void waves involve local heteroge- neous bubble clustering, which was confirmed statistically in a pattern analysis of bubble arrangement. Because void waves stand out when drag reduction becomes effective, we attempted mathematically to derive a wave equation for void fraction. The wave equation originates from a con- tributing term to drag reduction that depends on the spa- tial gradient of void fraction. This term also constrains the propagating speed of the void fraction wave relative to the ship speed.

Throughout this experiment, our detailed measure- ments of bubble behavior have suggested spin-off ques- tions regarding bubbly two-phase turbulent boundary layer dynamics during drag reduction. We leave these issues as open problems from which we expect further advances in scientific understanding as well as practical improvements in ship drag reduction.

‌Acknowledgements This work was supported by the Fundamen- tal Research Developing Association for Shipbuilding and Offshore (REDAS), Grant-in-Aid for JSPS Fellows No. 15J00147, JSPS KAK- ENHI Grant Nos. 24246033 and 23760143, and Grant-in-Aid for Young Scientists (B) No. 16K18006. The authors express their appre- ciation for all the support. Also, the authors express thanks to Prof. Yasukawa of Hiroshima University for his full support during the towing experiments at the Graduate School of Engineering.

‌References

‌Amromin E, Karafiath G, Metcalf B (2011) Ship drag reduction by air bottom ventilated cavitation in calm water and in waves. J Ship Res 55:196–207‌

‌Ceccio SL (2010) Frictional drag reduction of external flow with bub- ble and gas injection. Annu Rev Fluid Mech 42:183–203

Foeth EJ, Eggers R, Quadvlieg EHHA (2010) The efficiency of air- bubble lubrication for decreasing friction resistance, Paper no.

‌12. Prof Int Conf Ship Drag Reduction (SMOOTH-SHIP), Istan- bul, Turkey

Fukuda K, Tokunaga J, Nobunaga T, Nakatani T, Iwasaki T (2000) Frictional drag reduction with air lubricant over a super-water- repellent surface. J Mar Sci Technol 5:123–130

‌Hamilton JM, Kim J, Waleffe F (1995) Regeneration mechanisms of near-wall turbulence structures. J Fluid Mech 287:317–348

‌Harleman MJW, Delfos R, Terwisga TJCV, Westerweel J (2011) Dis- persion of bubbles in fully developed channel flow. J Phys Conf Ser 318:052007

‌Hinze JO (1955) Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. AIChE J 1:289–295

‌Hirayama A, Soejima S, Miyata H, Tatsui T, Kasahara Y, Okamoto Y, Iwasaki Y, Shimoyama N (2003) A study of air lubrication method to reduce frictional resistance of ship—an experimental study using flat plate and 16 m-model. In: West-Japan Society of Naval Architects meeting, vol 105, pp 1–9 (in Japanese)‌

‌Jang J, Choi SH, Ahn S, Kim B, Seo JS (2014) Experimental investi- gation of frictional resistance reduction with air layer on the hull bottom of a ship. Int J Nav Archit Ocean Eng 6:363–379

‌Johansen J, Castro AM, Carrica P (2010) Full-scale two-phase flow measurements on Athena research vessel. Int J Multiphase Flow 36:720–737

‌Katsui T, Okamoto Y, Kasahara Y, Shimoyama N, Iwasaki Y, Soejima S (2003) A study of air lubrication method to reduce frictional resistance of ship: experimental investigation by tanker form model ship and estimation of full scale ship performance. J Kan- sai Soc Nav Archit Jpn 239:45–53 (in Japanese)

‌Kitagawa A, Sugiyama K, Murai Y (2004) Experimental detection of bubble–bubble interactions in a wall-sliding bubble swarm. Int J Multiphase Flow 30:1213–1234

‌Kodama Y, Kakugawa A, Takahashi T, Kawashima H (2000) Experi- mental study on microbubbles and their applicability to ships for skin friction reduction. Int J Heat Fluid Flow 21:582–588

‌Kumagai I, Takahashi Y, Murai Y (2015) A new power-saving device for air bubble generation using a hydrofoil for reducing ship drag: the- ory, experiments, and applications to ships. Ocean Eng 95:183–194

‌Lahey RT Jr (1991) Void wave propagation phenomena in two-phase flow. AIChE J 37:123–135

‌Lammers JH, Biesheuvel A (1996) Concentration waves and the instability of bubbly flows. J Fluid Mech 328:67–93

‌Latorre R, Miller A, Philips R (2003) Micro-bubble resistance reduc- tion on a model SES catamaran. Ocean Eng 30:2297–2309

‌Lisseter PE, Fowler AC (1992) Bubbly flow—II: Modelling void frac- tion waves. Int J Multiphase Flow 18:205–215

‌Madavan NK, Deutsch S, Merkle CL (1984) Reduction of turbulent skin friction by microbubbles. Phys Fluids 27:356–363

‌Mäkiharju SA, Elbing BR, Wiggins A, Schinasi S, Vanden-Broeck JM, Perlin M, Dowling DR, Ceccio SL (2013) On the scaling of air entrainment from a ventilated partial cavity. J Fluid Mech 732:47–76

‌McCormick M, Bhattacharyya R (1973) Drag reduction of a submers- ible hull by electrolysis. Nav Eng J 85:11–16

‌Mercado JM, Gómez DC, van Gils D, Sun C, Lohse D (2010) On bubble clustering and energy spectra in pseudo-turbulence. J Fluid Mech 650:287–306

‌Mizokami S, Kawakado M, Kawano M, Hasegawa T, Hirakawa I (2013) Implementation of ship energy-saving operations with Mitsubishi air lubrication system. MHI Tech Rev 50:44–49

‌Mori K, Imanishi H, Tsuji Y, Hattori T, Matsubara M, Mochizuki S (2009) Direct total skin-friction measurement of a flat plate in zero-pressure-gradient boundary layers. Fluid Dyn Res 41:021406

‌Murai Y (2014) Frictional drag reduction by bubble injection. Exp Fluids 55:1733

‌Murai Y, Matsumoto Y (2000) Numerical study of the three-dimen- sional structure of a bubble plume. J Fluids Eng 122:754–760

‌Murai Y, Fukuda H, Oishi Y, Kodama Y, Yamamoto F (2007) Skin friction reduction by large air bubbles in a horizontal channel flow. Int J Multiphase Flow 33:147–163

‌Nierhaus T, Vanden Abeele D, Deconinck H (2007) Direct numerical simulation of bubbly flow in the turbulent boundary layer of a horizontal parallel plate electrochemical reactor. Int J Heat Fluid Flow 28:542–551

Oishi Y, Murai Y (2014) Horizontal turbulent channel flow interacted by a single large bubble. Exp Thermal Fluid Sci 55:128–139

‌Oishi Y, Murai Y, Tasaka Y, Takeda Y (2009) Frictional drag reduc- tion by wavy advection of deformation bubbles. J Phys Conf Ser 147:012020‌

‌Park HJ, Tasaka Y, Oishi Y, Murai Y (2015a) Drag reduction promoted by repetitive bubble injection in turbulent channel flows. Int J Multiphase Flow 75:12–25

‌Park HJ, Tasaka Y, Murai Y (2015b) Ultrasonic pulse echography for bubbles traveling in the proximity of a wall. Meas Sci Technol 26:125301

‌Pauchon C, Banerjee S (1988) Interphase momentum interaction effects in the averaged multifield model. Part II: Kinematic wave and interfacial drag in bubbly flows. Int J Multiphase Flow 14:253–264

‌Sanders WC, Winkel ES, Dowling DR, Perlin M, Ceccio SL (2006) Bubble friction drag reduction in a high-Reynolds-number flat- plate turbulent boundary layer. J Fluid Mech 552:353–380

‌Schlichting H (1979) Boundary-layer theory, 7th edn. McGraw-Hill Higher Education, New York

‌Smith CR, Metzler SP (1983) The characteristics of low-speed streaks in the near-wall region of a turbulent boundary layer. J Fluid Mech 129:27–54

Takagi S, Matsumoto Y (2011) Surfactant effects on bubble motion and bubbly flows. Annu Rev Fluid Mech 43:615–636

‌Takahashi T, Kakugawa A, Makino M, Kodama Y (2003) Experimen- tal study on scale effect of drag reduction by microbubbles using very large flat plate ships. J Kansai Soc Nav Archit Jpn 239:11– 20 (in Japanese)‌

‌Takeuchi T, Kagawa T (2013) Applicability of frequency response test for stability evaluation of gas pressure regulator. Trans Soc Instrum Control Eng 49:747–754 (in Japanese)

Titov I (ed) (1975) Practical problems in ship hydromechanics.

‌Sudostroeniye Publishing Hose, Leningrad (in Russian)

‌Tokunaga K (1987) Reduction of frictional resistance of a flat plate by microbubbles. In: West-Japan Society of Naval Architects meet- ing, vol 73, pp 79–82

‌Watanabe O, Masuko A, Shirose Y (1998) Measurements of drag reduction by microbubbles using very long ship models. J Soc Nav Archit Jpn 1998:53–63

‌Yim KT, Kim H (1996) On the variation of resistance components due to air bubble blowing on bulb surface of a ship. Trans SNAK 33:54–64 (in Korean)

Zacksenhouse M, Abramovich G, Hetsroni G (2001) Automatic spa- tial characterization of low-speed streaks from thermal images. Exp Fluids 31:229–239