Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
LES of Bubble Dynamics in Wake Flows T. Celik, A. Smirnov, S.Shi West Virginia University Department of Mechanical arid Aerospace Engineering Morgantown, WV26506-6106 Email: andrei @ smirnov.mae.wvu.edu Nomenclature DNS LES LPD NS RANS RFG SGS ~ Abstract Direct Numerical Simulation Large Eddy Simulation Lagrangian Particle Dynamics Navier Stokes Reynolds Averaged NS equation Random Flow Generation Sub-Grid Scale Length scale of turbulence Time scale of turbulence In this paper we present the results obtained using the Lagrangian particle dynamics (LPD) method and a ran- dom flow generation technique (RFG) developed by the authors in conjunction with large eddy simulations (LES) in application to turbulent bubbly wake flows. The hy- brid LES+LPD+RFG approach was applied to the case of a two-phase bubbly mixing layer and the high-Reynolds number bubbly ship-wake flows. The simulations were performed for the wake of a widely used ship Navy model 5415. We also present the results of a validation study where the same method is applied to the bubbly mixing layer experiments. The dynamics of the vortices and bub- ble concentrations reproduced in simulations are in a good agreement with experimental data. 2 Introduction Lagrangian particle dynamics (LPD) method is practi- cal and efficient in computing dispersed phases in multi- phase flow systems. It is usually favored over Eulerian two-fluid methods in the cases with dilute suspensions or large concentration variability of the dispersed phase (El- ghobashi, 1994; Crowe, 1998~. This is true for bubbly wake flows, such as ship-wakes, where bubbles experi- ence preferential concentration and clustering effects in the near wake region and are rather dilute in the far wake. The LED model is also conceptually simple and method- ologically robust providing a numerically stable statisti- cal method for evaluating dispersed phase statistics and particle/bubble distribution functions. The algorithm for particle tracking and population dynamics developed by the authors demonstrated the ability to efficiently simu- late large populations of particles including coalescence effects with even modest computer resources (Shi et al., 2000a; Smirnov et al., 2000; Smirnov and Celik, 2000~. Simulations of bubbles in turbulent shear layers re- quire accurate representation of the fluctuating flow-field that governs bubble dynamics. When conventional RANS (Reynolds Averaged Navier-Stokes) models are used this accuracy is often lost or comes at a cost of empirical phys- ical sub-models, e.g. turbulence. On the other hand, direct numerical simulations (DNS) or LES require no or a relatively simple subgrid-scale (SOS) model. LES has the advantage of being able to handle flows with higher Reynolds numbers and still accurately reproduce large-scale turbulent structures (eddies), which are im- portant in bubble dynamics. The majority of simula- tions of turbulent bubbly shear layers using Eulerian-
Lagrangian approach is done using RANS methods for (Ad) and lift (Al) given by the following expressions high Reynolds number flows (Joie et al., 1997; Mu- rai, 2000) or DNS-type simulations for relatively low Reynolds number flows (Elghobashi and Lasheras, 1996; Ruetsch and Meiburg, 1994; Okawa et al., 2001; Murai et al., 2001). An LES-based approach pursued in this study provides the basis for multi-phase LES of large- Reynolds number flows, like ship-wakes. A major disadvantage of the Lagrangian approach is the perception that an extremely large number of particle trajectories is needed to obtain sufficiently smooth statis- tics. This study elaborates that this is not necessary true as long as the flow can be categorized as "nominally dilute". 3 Method The LED algorithm for particle tracking and population dynamics developed by the authors demonstrated the abil- ity to efficiently simulate large populations of particles in- cluding coalescence effects with even modest computer resources (Shi et al., 2000b; Shi et al., 2000a; Smirnov et al., 2000; Smirnov and Celik, 2000~. In this study the LED algorithm was combined with the Large-Eddy sim- ulation approach (Smirnov et al., 2001a; Piomelli, 1999) to compute the distribution of bubbly phase in the near- wake flow of a ship-model. The LES technique (Piomelli, 1999) was enhanced with The RFG procedure (Smirnov et al., 2000; Shi et al., 2000b) was implemented into the LES method to enable the appropriate representation of initial/inlet turbulence conditions and as a subgrid scale for particle dynamics solver. The combined LES/RFG method was validated on the case of turbulent mixing layer, and used later to simu- late the wake of a Navy 5415 ship model as described in (Smirnov et al., 2001a). The fluctuating velocities at the inflow boundary provided by RFG were generated from the data obtained in prior RANS calculations (Larreteguy, 1999~. To compute bubble dynamics the following equations of motion were adapted from Shridhar and Katz (1995) ~Ub =Aa+Ab+Ad + A~ (1) A a 3 (aatf + (Uf V)Uf) Ab = -2g Ad = 4r Cd | Urel | Ure Al = 4 Cl | Urel I U ! (o ~ c') = 0.5*(VxUf) where rb is bubble radius, Uf is fluid velocity, Urel is the relative velocity of fluid and bubble, i.e. UfUb, and Cd ,Cl are drag and lift coefficients respectively. The above relations were obtained for spherical bubbles and by considering the limit of small air-to-water density ra- tios Pb/PI ~ 0. The coefficients of drag (Cd) and lift (cl) are themselves empirical functions of Urel and bubble di- ameter, which were obtained from experimental measure- ments of bubbles dynamics in turbulent vortices. The Bas- sett term involving the history integral was neglected in this case, following the conclusions of Shridhar and Katz (1995) on small contribution of this term compared to the buoyancy term. Equation (1) was discretized using a second-order Runge-Kutta time-stepping procedure. Bubble tracking throughout the computational domain was simplified by using a uniform Cartesian grid, so that the index of the cell containing the bubble could easily be obtained using sim- ple division by modulus operation. Fluid velocities were interpolated to the location of a bubble using tri-linear in- terpolation formula. Bubble dissolution effects were accounted for by the following equation (Hymen, 1994; Carrica et al., 1998) dm 6~C Co)k2/3 ~ u ~1/3 r4/3 (2) where r is the bubble radius, ~ u ~ is the relative air- bubble velocity and k is the diffusivity of air in water. Constants CoO7Csurf represent the concentrations of dis- solved gas at a distant point and at the bubble surface re- spectively and were derived using Henry's law (Carrica et al., 1998) for air solubility in water. where the A vectors on the right-hand side represent ac- In a joint LES+LPD+RFG approach the flow and celerations due to added mass (Aa), buoyancy (Ab), drag the particle solvers were using different time-stepping 2
schemes with independent selection of time-steps. Usu- 7c~is ally sub-cycling of particle iterations was required to reach a stable solution. 4 Results 4.l Validation on a turbulent mixing layer. The validation of the LES/LPD approach was done on ex- perimental data of a mixing layer (Rightley and Lasheras, 2000) (Fig.1~. The computational domain size is 0.55m x 0.2m x 0.2m in stream-wise, vertical, and span-wise di- rections respectively. The mixing layer is generated by two separate parallel flows with different incoming veloc- ities. A thin flat plate, which is O.l5m length, 0.003m height, and 0.2m width (the whole span-wise extent), is mounted in the middle of the inlet plane. The velocity of the lower half flow is 0.28m/s, while the upper half value is 0.07m/s. Bubbles of 40,u diameter are carried in from the lower half. As it was the case in the experiments, a sinusoidal perturbation with the amplitude of five per- cent of the mean flow, was added to the vertical velocity component of the flow in the lower half of the inlet plane during the simulations. The grid used in this simulation was a 194 x 66 x 42 uniform grid on all three directions. The cell size is 2.8mm x 3.0mm x 4.8mm. Free gradient boundary condition was used for the outflow boundary. Slip-wall boundary conditions were used at the top and bottom (y-direction). In the span-wise direction (z), pe- riodic boundary conditions were applied. At the surface of the flat plate, slip-wall boundary conditions were used. The reason for using slip-wall boundary conditions is that in this computation the exact resolution of the boundary layer was outside of the scope of the current study. More- The predicted rms velocity in axial direction is gener- over the mixing process is dominated by the effect of the ally lower than the measured values. This leads to under- shearlayer. estimation of turbulence intensities. Farther away from the region where the mixing occurs, the difference is even larger. The slip-wall boundary condition at the top and bottom boundaries may be the primary reason for the ob- served differences. Very little mixing occurs in these re- gions, so the turbulence intensity is almost zero. A simi- lar phenomenon is observed for the vertical mean veloci- ties and their fluctuations (Fig. 2(b)~. It should be noted that the rms values shown in Fig.2 represent only the re- solved part of velocity field reproduced by LES, which Figure 1: Mixing layer experimental setup. locity, u' = t(2/3)k]~l2, which is used as an input to the RFG method: u' = (2/3) ~l2(cs/c~u~s/\2/(a~m) (3) where Cs = 0.04, C,~' = 0.09 are the turbulence con- stants, S is the strain rate tensor computed at each point, /\ = 3 mm is the grid cell size, 6,'~ is the mixing layer thick- ness, and a = 0.035 is an empirical constant. Because the volume fraction of the bubbles was kept small (~ 5 10-4), the effect of bubbles on the carrier flow was not considered. This approximation is also supported by Elghobashi and Truesdell (1993) where they indicated that the influence of such low dispersed phase particles (bubbles) would be seen in the smallest scales, not in the large, energy-containing scales dominating the mean and rms velocity profiles. Figure 2(a) shows the stream-wise mean velocities and the fluctuations at stream-wise position x/X = 1.25, where ~ = 79.36mm is the length of Kelvin-Helmholtz instabil- ity wave, representing the characteristic flow length scale. The profile of the mean stream-wise velocity is in a good agreement with the measurements. Central differencing (CD) scheme was used to dis- cretize the convection term. The the Smagorinsky SGS model was used to represent the effects of subgrid scale motion on the carrier flow and the RFG technique (Smirnov et al., 2001 b) was applied to reproduce the effect of turbulence on the bubbles. In the latter case the standard relation for turbulent eddy viscosity, v, = C,Ukl/2 * Im''r' and the Smagorinsky model v, = CsA2S were used to provide the estimate of the fluctuating ve- 3
domain with time interval of 27 milliseconds for the du- the boundary layer. Some irregularity of computed data ration of the 12.7 seconds, resulting in 47,000 bubbles that further downstream is due to statistical uncertainty of the entered the domain during the course of the run. Fig- sample as bubble concentration becomes more dilute. ure 5 shows a typical bubble distribution computed by It should be noted that, being statistical in nature, the LES/LPD algorithm compared with those observed in ex- LED method enables the refinement of the histogram and periments. The bubble cloud can be seen to be entrained improvement of the accuracy by subsequent accumulation by the fluid entering the mixing region from the high of bubble statistics. Since the statistical error is propor- speed side into the cores of the coherent vertical flow tional to 1/ni/2, where n is the number of bubbles in a structures present in the mixing region. The pictures from histogram box, it would require four times as many bub- the simulations represent an instantaneous distribution, bles for a two-dimensional histogram to double the reso- which is different for any given time, but shows common lution along each axis, or to reduce the statistical error by statistics and similar dynamical features as seen in the ex- half with the same resolution. The choice of histogram periments. size is a trade-off between the spatial resolution and sta- Spatial histograms of bubble distributions were ob- tistical error. rained by counting all bubbles passing the cells of the 3D-histogram and accumulating the statistics. Normaliz- ing the histogram data by the total number of bubbles in- jected, gave bubble probability density functions and con- centrations. The statistical error, a, in the number of bub- ble counts, n, for each slot of the histogram can be esti- mated from binomial distribution as c, = En(~1n/N)] ~/2, where N is the maximum number of bubbles injected dur- ing the simulation. This will add to the uncertainty in the layer thickness calculations for larger X, and smaller bub- ble concentrations. Figured shows the surface of the histogram and the contour plot of the constant concentration levels. Fig- ure 7 shows the growth of the mixing layer thickness in the stream-wise direction (X) with and without the RFG model. Over-prediction of layer thickness at small X is due to the finite resolution of the histogram. The inclusion of the RFG model improves the predictions at greater X. Figure7 also shows predicted and measured develop- ment of the mixed layer thickness. The overall agreement is good except at large axial distances where a smaller rate of growth of the mixing layer in the stream-wise di- rection is observed when no subgrid scale model was ap- plied. Although the inclusion of the RFG model (Smirnov et al., 2001 b) improved the predictions of the boundary layer growth at higher axial asymptotics, it also increased the uncertainty in the boundary layer thickness because of the higher dispersion in bubble distribution. Increasing the bubble sample will improve the accuracy. Figure8 provides the comparison of predicted bub- ble concentrations against experimental data. The agree- ment is reasonably good especially at the beginning of s 4.2 Ship-wake simulations. The ship-wake simulations were performed for the wake of a widely used ship Navy model 5415 ~ (Smirnov et al., 2001a; Carrica et al., 1998~. In the simulations a com- putational 164 x 98 x 92 grid was used to represent the near-wake region of one ship-length L in the axial direc- tion, 0.3L in depth and 0.6L in the span-wise direction. The RFG procedure was used to initialize turbulent fluctu- ations in the inlet plane (Smirnov et al., 2001b). The free surface was assumed to be flat, invoking the low Froude number approximation. Bubbles with a uniform size of lOO,u were randomly injected at the inlet plane. The prob- ability distribution of the injection points was set propor- tional to the turbulent kinetic energy level obtained from the preliminary RANS solution. At the free surface the life-time of the bubbles was set to zero for the absence of more exact empirical data. The integration of Eq.~1), via a second-order Runge-Kutta scheme provides the new velocity, vitt), in the xi location for each particle as a func- tion of time. Figure 9 shows an instantaneous picture of the bubbles in the wake as observed from the rear-bottom corner of the computational domain (ship stern is not shown). It can be seen that the bubbles experience a trend to a preferential concentration regions. These regions are dictated by the dynamics of the flow in the near wake. Figure 10 provides a cross-sectional view of a instantaneous velocity distribu- tions in the wake, which is rather typical for ship wakes ~ htip://wwwSO.~.navy.mil/54 ~ S
(Hymen, 2000~. Two symmetrical vertical structures in the velocity distribution coincide with the bubble cluster- ing regions. With the higher vortex intensities, which will take place at a higher Reynolds number, bubble distribu- tion in the wake will be increasingly affected by these vor- 5 tices. This is especially true for smaller bubbles. Figurell shows typicalcross-sectionaldistributionsof The results of this study show the viability of joint bubbles in the wake. These distributions are not instanta- LES/LPD method for computing turbulent bubbly wakes, neous, but represent total counts of bubbles accumulated which can be applied to high Reynolds-number ship-wake over the simulation run. flows. Realistic bubble distributions for mixing layers and wakes were obtained. The effect of bubble influence on the flow field was not considered in this study. Nevertheless, the work by others (Elghobashi and Truesdell, 1993; Truesdell and Elghobashi, 1994; Elghobashi and Lasheras, 1996) indi- cates that when there is a large density ratio between the phases, such as that in this study, the influence of bubbles on the carrier phase may become important even though the void fraction remains small. Work on this topic are continuing under the current project. To further improve the accuracy of predictions for the flat-plate wake a more realistic no-slip wall boundary con- dition may be used together with the grid refinement at the wall so as to resolve the turbulent boundary layer. More- over a more realistic model of bubble/turbulence interac- tion near the free-surface needs to be implemented in the future. for smaller bubbles (~ 50,u) and at lower depths (~ 10m). This was also confirmed by LES simulations. Conclusions and future work PLANE COUNT RATIO 1 46998 0.01 2 9803 0.21 3 3326 0.07 4 1119 0.02 Table 1: Cumulative bubble distributions in different planes Table 1 shows the total bubble counts in different planes Ci=l,4 and the corresponding normalized ratios: Ri = Cinch. These distributions indicate a strong depletion of bubbles over the half ship-length distance, shown in Fig. 13. The depletion of bubbles is due to the buoyancy driven migration to the surface and the dissolution effects. Buoyancy forces affect mainly large bubbles, whereas the dissolution affects primarily small ones. Therefore, even the small bubbles entrained by the vortices will eventually References disappear from the domain. Figures 12 show the contours of probability density functions of bubble occurrence in the wake. These probabilities were computed from bubble-distribution his- tograms obtained in the way as for the mixing layer vali- dation case described above. The profiles shown in Fig.12 indicate a rapid bubble population decay and gradual spreading of the bubble cloud on the distance of one ship- length. Although the classified nature of bubble measure- ment data in the wakes of Navy ships prevents us from making a direct comparison, we believe that the predicted bubble distributions are similar to those observed in typi- cal ship wakes (Hymen, 1998; Hyman, 20001. Separate estimates of bubble dissolution rate have Carrica, P., Bonetto, F., Drew, D., and Lathey, J.: 1998, The interaction of background ocean air bubbles with a surface ship, Int. J. Numer. Meth. Fluids 28, 571 Crowe, C.: 1998, An assessment of multiphase flow models for industrial applications, in Proceed- ing of FEDSM'98, Vol. FEDSM-5093, Washing- tow,DC,USA Elghobashi, S.: 1994, On predicting particle-laden turbu- lent flows, Applied Scientific Research 52, 309 Elghobashi, S. and Lasheras, J.: 1996, Effects of Grav- ity on Sheared Turbulence Laden with Bubbles or Droplets, in 3-rd Microgravity Fluid Physics Confer- ence, Cleveland, OH shown that its effect can be neglected in the case of mix- Elghobashi, S. and Truesdell, G.: 1993, On the two- ing layer and can be noticeable in the case of ship-wake way interaction between homogeneous turbulence and 6
dispersed solid particles. i: Turbulence modification, Physics of Fluids A 5, 1790 Hyman, M.: 1994, Modeling Ship Microbubble Wakes, Technical Report CSS/TR-94/39, Naval Surface War- fare Center. Dahlgren Division. Hyman, M.: 1998, Computation of ship wake flows with free-surface/turbulence interaction, in 22nd Sympo- sium on Naval Hydrodynamics, pp 11-32, Washing- ton,D.C. Hyman, M.: 2000, Bubble Concentrations Near a Ship Surface, Private communication Vat de Reuil, France Smirnov, A. and Celik, I.: 2000, A Lagrangian particle dynamics model with an implicit four-way coupling scheme, in The 2000ASME International Mechanical Engineering Congress and Exposition. Fluids Engi- neering Division, Vol. FED-253, pp 93-100, Orlando, F1 Smirnov, A., Shi, S., and Celik, I.: 2000, Random Flow Simulations with a Bubble Dynamics Model, in ASME Fluids Engineering Division Summer Meeting, No. 11215 in FEDSM2000, Boston, MA Joia, I., Ushima, T., and Perking, R.: 1997, Numeri- Smirnov,A., Shi,S.,andCelik,I.: 2001a, Lesof a bubbly cat study of buble and particle motion in a turbulent boundary layer using propoer orthogonal decomposi- tion, Applied Scientific Research 57, 263 Larreteguy, A.: 1999, Ship-Wake simulations, Private communication Murai, Y.: 2000, 1 Y' . 1 ~ ~~ ~ Numerical study of the three- dimensional structure of a bubble plume, Transactions of ASME 122, 754 Murai, Y., Ohno, Y., Bae, D., Abdulmouti, H., Ishikawa, M., and Yamamoto, F.: 2001, Bubble generated con- vection in immiscible two-phase stratified liquids, in ASME FEDSM-01, New Orleans, LA Okawa, T., Nakazumi, M., Yoshida, K., Matsumoto, T., and Kataoka, I.: 2001, Interfacial forces acting on a bubble in verticla upflow, in ASME FEDSM-01, New Orleans, LA Piomelli, U.: 1999, Large-eddy simulation: achievements and challenges, Progress in Aerospace Sciences 35, 335 Rightley, P. and Lasheras, J.: 2000, Bubble dispersion and interphase coupling in a free-shear flow, Journal of Fluid Mechanics pp 21-59 Ruetsch, G. and Meiburg, E.: 1994, Two-way cou- pling in shear layers with dilute bubble concentrations, Physics of Fluids 6~8), 2656 Shi, S., Smirnov, A., and Celik, I.: 2000a, Large eddy simulations of particle-laden turbulent wakes using a random flow generation technique, in ONR 2000 Free Surface Turbulence and Bubbly Flows Work- shop, pp 13.1-13.7, California Institute of Technol- ogy, Pasadena, CA Shi, S., Smirnov, A., and Celik, I.: 2000b, Large-Eddy simulations of turbulent wake flows, in Twenty-Third Symposium on Naval Hydrodynamics, pp 203-209, sn~pwa~e now, ~n bympos~um on CFD Applications in Aerospace, 2001 ASME Fluids Engineering Divison Summer Meeting, No. 18013 in FEDSM-2001, New Orleans Smirnov, A., Shi, S., and Celik, I.: 2001b, Random flow generation technique for large eddy simulations and particle-dynamics modeling, Trans. ASME. Journal of Fluids Engineering 123, 359 Sridhar, G. and Katz, J.: 1995, Drag and lift forces on mi- croscopic bubbles entrained by a vortex, Phys. Fluids 7(2), 389 Truesdell, G. and Elghobashi, S.: 1994, On the two-way interaction between homogeneous turbulence and dis- persed solid particles. ii: Particle dispersion, Physics of Fluids 6, 1405 7
At: Ad := :" - ~ O . 03~ O . 0 1- _ O~ _.: i. _~. A: · ~ i_ ~ C~ t: ~ ~ ~ ~ ~ : ' -] ~~ -I q ~~ :~- ~ ~ ~' ~ :~: by, ~ ~ ~ ~ `~ ~ I: ~ ~ ~ ~' :` : ,, :~ ~~::J~ ::~;:~: 7~ ~ ~ ~ :` Or .' ~ ~ ~ ~ ~.~i hi ~ ~ ~ ~ ~ ~ by: Aft: t. ~ :§ ~ ~ ~ ~ ~ ~ i: .. ., : _- ~ g: : ~ : I: 0.04~ . A: ~ ~ ~ id. A.: :~ : L., ~~~~ i. - =~C~ : :: : : :::: :-: : as, :: :: :: ::::: ~ ::: ::::: : : ~ : : ad. : :: ::: ::~::: :: A::.. :::: : : ~ : : :: ::~: ' ~ . :: : I::: : :: ~ : ::::: : : : : : . : ::: ,: . ~ - o': : ~ ~: -~ ~ ~: ~ ~, ~ :~ ~ ~ ~ ~ ~:~ :~: ~ '~ : :~: '' ~ Ph~:~ ~ t r :~ ~ ': ~ :^ ~~: : :~ : ~ a. ~ ~ :Ir ~ : ~> ~L. ' AL 1 (a) Measurements by Rightley (1995) . . . . . . ~ . . _ _ ~ ~ ~ ~ ~ t t t - ' t ~ ~ ~ ~ ~ j ~ 6, . ~ t ~ ~ ~ ~ t ~ ~ ~ I ~ ~ ~ ~ ~ ~ ~ ~ t, ~ t ~ t t t t O_ ~ ~ 27 1 0~129 1 0~31 1 0~33 1 0~35 1 0~37 1 0~39 1 0~41 1 0~43 0 180 360 540 720 x (m) / Phase Angle (degree) (b) Simulations Figure 3: Velocity vectors at x/\ = 2.50 8
0.14 0.05] ~ at 0 . 05- ~9 Figure 4: Computed vertical velocity contours. (a) Computed (b) Rightley and Laseras (2000) Figure 5: Instantaneous bubble distribution 9 Level v 23 8.5 E-02 21 6.3 E-02 1 9 4.0E-02 1 7 3.0E-02 1 5 2.4E-02 1 3 -2.6E-02 1 1 -3.5E-02 9 -5.0E-02 7 -6.1 E-02 5 -8.3E-02 3 -1.1 E-01 1 -1 .3E-01
'PDF' 1 .2 0.8 0.6 0.4 0.2 o 4 X/L ~/~//A///~//// 3 0\ -1 ~ 0.8 ~0.6 j/~/~/ 60 40 \/e°rtical distance [mm] -60 Figure 6: Normalized histogram of probability density function of bubble distributions. 60 50 40 cot c,, 30 ._ cat .= 20 > 10 o ... ^-'3X<;'l"'"~""' 0 0.5 aid...+ Y.` - ~ at, ..~.., . I;;; ;! ~,..~. - ,,,,,...,~,,~ 1 ~ 1__;~ ..,d'.~, 1 1.5 ~ Streamwise distance [X/L] I. .--- ,.H ~ ~ ~ ............... =.,,)~..~.~ _ 2.5 3 Figure 7: Mixing layer thickness. 3 measured; * - computed without RFG; x - computed with RFG 10
0.8 0.6 0.4 0.2 o 1.2 , , . . , 1 . 1.2 .--\\ Doer , . . . , \ · · O -40 -30 -20 -10 0 10 20 30 Vertical distance [mm] (a) x/=0.31 40 -40 -30 -20 -10 0 10 20 30 40 Vertical distance [mm] (b) x/0.63 0.8 0.6 0.4 0.2 o A . ~ A \ ~ . ''. . ~ 'I.;, ~ .A ~ -40 -30 -20 -10 0 10 20 30 Vertical distance [mm] (c) x/=1.25 1.2 0.8 0.6 0.4 0.2 o 1 ~ 0.8 0.6 0.4 0.2 n , , . , , , . 1 .1 1 1 At' . 0.9 \~ 0.8 ~ \ it/\ o76 >- 0.3 , , , , , , me,, of o - 0 -30 -20 -1 0 0 1 0 20 30 40 Vertical distance [mm] (e) x/=1 .88 Figure 8: Bubble concentrations (~0.079365 m). -+- experiments, computations. 11 40 -40 -30 -20 -1 0 0 1 0 20 30 40 Vertical distance [mm] . . . . · ~ ~ ~ \ Hi\ A\ · iK (d) /1.56 .. ~ \ ;' -40 -30 -20 -1 0 0 1 0 20 30 40 Vertical distance [mm] (f) x/2.50
!: .j I.. Figure 9: Bubbles in a ship-wake. Rear bottom view. -o.o,~ -0.03 -0.04 -0.05 -0.06 -0.07 -0.08 ~ } ~ I, I, ~ ~ ~ ~ ~ ~ ~ , , , , , , ~ , ~ , . , , ~ . . ., it,, ~ -0.05 0 0.05 Figure 10: Velocity vectors at X/L = 0.2 12
-0.01 -n nob -0.01 -0.02 -0.025 -0.03 -0.04 -0.045 -0.05 -0.1 -0.05 n (a) X/L = 0.01 ~ I-:' -0.015 -n.025 -n no -0.045 -0.05 + n ns 0.1 -0.1 item ~ + ++++ + +++ *I+ +*+ +$+ + be+ ++ ti , 1 11$ -0.1 -0.05 0 0.05 (c) X/L = 0.50 t3 .,.~. ~ ~ -0.05 0 0.05 0.1 (b) X/L = 0.25 -0.015 -0.02 -t -0.025 .+ -0.03 -0.035 . -0.04 -0.045 ;+~ At + ++ +` +++ ++{ ++~++ l++ Be++ + $+~++~$+ + + ++++ + + ++++~+ + ++ +++ + ++ + +t~ +$+ t ++ + + $ + + + , , f+ 0.1 -0.1 -0.05 0 (d) X/L = 0.75 Figure 11: Bubble distributions at different cross-sections 13 0-05 0.1
O, ,.-~6;' , -0,005 my,,! 0.01 -0.015 0.02 0.025 -0.03 0.035 -0.04 0.045 -0.05 -0.1 0.05 0 0.05 0.1 0.02 --------- 0.015 0.01 ~ ohs . ~ . . (a) X/L = 0.01 O :~.-': ·00005 ooo2o53 hi;',. oooo4 0.045 . 0.05 1 , , -0.1 0.05 0 0.05 0.1 0.002 ------- ).0015 (b) X/L = 0.25 0.0008 ~ 1 ~ 0.00025 0.0007 0.0002 0005 0005 ooos 0015-- 55:, 0002 0OOO5~ : -O.M ~ '' <~ jOOO1 -- ·0 02 ~ ~ ~ 1 ooo45 ~ 1 .l ~ 1 004 t 'I-- ~ 0.05 -0.05 0.1 0.05 0 0.05 0.1 0.1 0.05 0 0.05 0.1 (c) X/L = 0.50 eO5 ~ (d) X/L = 0.75 Figure 12: Bubble probability density functions at different cross-sections 14
oc9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 1 1~' 1 1 1 1 1 1 - - - - - - 0.1 o - 1 , , , , _ ~ 0 0.1 0.2 0.3 0.4 X/L Figure 13: Bubble depletion in the wake 0.5 0.6 0.7 0.8 Figure shows the ratio of total number of bubbles in each plane to the number of bubbles in the first plane. 15
DISCUSSION P. Atsavapranee Naval Surface Warfare Center, Carderock, USA Do you have any physical experimental data to compare to your result of bubble distribution around 5415 hull? AUTHORS' REPLY The question reflects a legitimate concern that without proper validation the confidence level in the simulations will low. Unfortunately there is very little data available in the open literature on bubble distribution in ship wakes. The data known to us is not appropriate for comparison with the present simulations. Other sources of data seem to exist (Hymen, 2002) but they are classified. There is, however, an ongoing program under ONR to quantify bubble distribution around the ship model 5415, but not directly in he wake. It would be highly desirable to conduct some carefully planned experiments for measuring the bubble distribution in the wake of surface ships or ship models to supplement our computational work. References 1. Hyman, M. (2002) Private Communication. DISCUSSION R.E. Melnik Mississippi State University, USA The authors should be commended for their development of an effective approach for the numerical solution of the very difficult problem of predicting the bubble distribution in a turbulent wake, a problem of great current interest in the US Navy. The authors approach based on coupling their LED and RFG methods with an established numerical method for the numerical solution of the 3D unsteady Navier- Stokes equations with LES modeling, for dilute bubble concentrations, is potentially the most accurate approach actually implemented to-date. The authors provided some validation of their method through comparison with experimental data for a flat plate generated wake that showed good agreement for mean flow stream-wise velocity but otherwise relatively poor agreement for the rms velocity fluctuations. Although, these results provide encouragement in the overall accuracy of the method, much uncertainty remains on the question its quantitative accuracy. My discussion will focus on the uncertainty in the accuracy of the method. Inaccuracies are basically introduced from three sources: approximations employed in the physical models (LPD, RFG, LES), errors from the slip flow boundary condition, and numerical errors arising from the computational methods employed for the LED and LES solutions. Because these errors are highly intertwined it is really not possible to draw meaningful conclusions about the accuracy of the individual sources of error. This point was forcefully stated in the conclusions of the 1981 Stanford workshop on the validation of turbulence models in vogue in that timeframe. In that workshop, numerical solutions of the turbulent boundary layer and RANS equations, using many of the available turbulence models, were compared with experimental data for a wide variety of flow conditions and geometry. Since that time computing power has increased to the point that it has now become customary to check the accuracy of numerical solutions of the RANS equations through grid refinement studies. Indeed, some journals will not accept papers on numerical simulation unless accuracy assessment data are included in the paper. Today it is equally important to assess the level of numerical error before drawing conclusions on the adequacy of LES models. The authors of the present paper point to the need to add a no-slip condition and grid refinement to better resolve the turbulent boundary layer emanating from the flat plate. They also pointed to the need for a better model of free surface effects on bubble turbulence interaction. I think it is equally important to conduct grid refinement studies on the numerical solution of the LES equations, before drawing conclusions on the adequacy of the physical models employed in the computations. LES grid refinement studies are more difficult than those for the RANS equations, because of the need for much greater computer resources as well as the need to consider the two length scales that appear in numerical solutions involving LES models, namely the filter width, /\, employed in the LES model and the mesh spacing, h, used in the computational grid. This issue has been studied in great detail in Refs 1-4, listed below. These papers considered numerical solutions of various LES modeled equations, including the Smagorinsky model, for the relatively simple
problem of a temporal mixing layer. The Reynolds numbers employed in these computations was low enough to permit computation of highly resolved DNS solutions, which then served as a basis for assessing the accuracy of the LES solutions. Here, I briefly summarize the methodology and the main conclusions reached in those papers. The quality of the solution depends on the choices of the parameters, ~ and h. The papers considered various values of these parameters organized in terms of the ratio, r = i\ /h. If ~ is held fixed while h is reduced a grid resolved solution of the LES equations will be approached in the limit ho 0. An alternative approach is to consider the ratio, r, fixed, while decreasing h. This would yield a grid independent DNS solution as h approaches zero, although at considerably higher computational costs. The papers considered the convergence of both mean and fluctuating quantities and they also provided an analysis of the computational cost of the various choices of ~ and h. The principle conclusion of these studies is that the ratio, A/h, should greater than 2 ~ r >2) and for some other quantities be as large as 4 -5 in order to achieve grid converged solutions of LES equations. Most LES solutions in the literature, including those in the paper under discussion, employ a value of 1` /h = 1, which implies that these results are likely contaminated by numerical error. I realize that these conclusions, which are based on relatively low Reynolds number computations, and are for a different type of mixing layer and therefore, may not apply to the computations of the paper under discussion. Nevertheless the results clearly establish the need to demonstrate that numerical solutions with LES models are free of effects of numerical contamination. I will close my discussion of the paper with a question for the authors. The RFG procedure is employed to provide turbulence inlet and initial conditions as well as to provide a subgrid scale model in the particle dynamics solver. The subgrid scale model is used to provide a turbulent flow representation to account for the small scale dispersion of particles. It is possible to couple the RFG technique to a RANS solver over the entire flow domain and to thereby avoid the need to use an expensive LES solver, greatly reduce required computer resources. It would be very interesting to compare the accuracy and cost of such a procedure with that of the LES method used in the paper. Have the authors considered such an approach? References 1. Vreman B., Geurts B., & Kuerten, H., 1996 "COMPARISONS OF NUMERICAL SCHEMES IN LARGE-EDDY SIMULATION OF THE TEMPORAL MIXING LAYER" Intl J. Numer Meth. Fluids 22, pp. 297-31 1 2. Verman B., Geurts, B., and Kuerten H., 1997 " Large Eddy simulation of the turbulent mixing layer", J. Fluid Mech. 339, pp. 357-390 3. Geurts B. J. and Frolich, J., 2001 "Numerical Effects contaminating LES; a mixed story", in Modern Simulation Strategies for Turbulent Flow, edited by B. J. Geurts (Edwards, Ann Arbor, MI ), pp. 309-327 4. Geurts B., Frohlich J., 2002 " A framework for predicting accuracy limitations in large-eddy simulation", Physics of Fluids, 14, pp. L41-L44 AUTHORS' REPLY We appreciate the comments of this reviewer with regards to quantification of numerical uncertainty in LES (Large Eddy Simulation). The reviewer seem to appreciate that this is not a trivial issue. In LES calculations both the modeling errors (i.e. the inaccuracies resulting from the sub-gri-scale (sgs) model) and numerical (primarily the discretization) errors are directly proportional to the grid size. This is especially true in cases of implicit top-hat filters inherent to finite volume formulations such as that used in the present paper. As the grid resolution is refined a good LES should approach to direct numerical simulations (DNS). If very fine grid resolution is used in LES calculations such that the simulation results are close to those obtained from DNS it will not be economical. If on the other hand relatively coarse grids are used to take advantage of the LES technique, then the discretization errors are large and they pollute the sgs contribution. Hence the user is faced with a dilemma of what to use for an optimum accuracy and benefit from LES. Although the present authors are aware of this non-trivial problem, no attempt was made to determine the numerical uncertainty in the present LES calculations. However, work is underway by our group to develop a practical
methodology with which the quality of an LES can be measured. Preliminary results from this effort will be reported in an upcoming ASME Fluids Engineering Conference (see Celik et al, 2003). We thank the reviewer for making us aware of the relevant study by Geurts and Froehlich (2002~. We studied this work closely and found out that the method presented in that paper is not directly applicable to our case. Since no DNS results are available for our case, the calculation of the "exact" solution is not possible. Moreover, the activity parameter "s" suggested by Geurts and Froelich requires the calculation of volume averaged turbulent dissipation rate which inherently includes both the modeled dissipation and the numerical dissipation; segregation of he two is necessary but not easy. As mentioned above, this is the topic of a future work by our group. In this approach we attempt to estimate the percentage of resolved turbulent kinetic energy relative to the total; if more than 70% is captured we judge the LES as adequate. Preliminary calculations have revealed that our current LES results presented in this proceedings satisfy this criteria. Simulations with much finer grid resolution using parallel version of the same code will be presented in the future Symposia on Naval Hydrodynamics. As for the last question, no we did not consider running RAN S supplemented by RFG. There is a good reason for that. The reason for using LES, in the first place, was that RANS smeared out even the relatively large coherent vorticies which are important for bubble dynamics. RGF on the other hand deals essentially with fine scale turbulence. Hence, it would be counter productive to resort back to RANS. References 1. Celik, I.B., Cehreli, Z., Yavuz, I. (2003) "Index of quality for large eddy simulations." Proceedings of the 2003 ASME Fluids Engineering Division Summer Meeting, Honolulu, Hawaii, July 6-10. 2. Guests, B.J. and Froehlich, J. (2002) "A Framework for predicting accuracy limitations in large-eddy simulation," Physics of Fluids, Vol. 14, NO. 6, pp.L42-L44.
