Ismail B. Celik, Allen E. Badeau Jr., Andrew Burt and Sherif Kandil

Mechanical and Aerospace Engineering Department

West Virginia University, Morgantown, WV 26506-6106

E-mail: icelik@wvu.edu

Abstract: The computation of mixing in two component liquid-liquid flows using two-fluid models is difficult for large-scale industrial applications. This is due to two-fold increase in the number of differential equations to be solved with semi-empirical phase interaction terms. This task can be simplified if only the equations for the mixture are solved while the relative is prescribed empirically. However, this requires some information about the size of particles, e.g. fuel droplets in water. In this study, a set of viable transport equations are formulated for the mixture along with a droplet formation model (DFM) suitable for liquid-liquid flows. The DFM is based on local flow parameters such as turbulent Richardson numbers involving local turbulence length and velocity scales. Parallel experimental study was conducted at Johns Hopkins University (JHU) by Katz et. al. The model constants in DFM are calibrated using data from JHU in conjunction with the CFX and FLOW-3D (commercial CFD codes) calculations. In various simulated cases, the model predicted fuel droplet size distributions very similar to experimental observations.

 

Keywords: droplet formation, shear flow, impinging jet, immiscible stratified liquids

The present study is an element of a larger work, which is being conducted at West Virginia University and at the Naval Surface Warfare Center – Carderock Division (NSWC-CD) in Bethesda, MD, whose goal is to study the transient flow phenomena that occur during the refueling process of compensated fuel/ballast tanks (CFBT’s), which are used in U.S. naval surface ships. Of primary interest in these studies is the location and extent of fuel/water mixing, the amount of water trapped (water hideout) in the tanks after refueling is complete, and the estimated flow-through time of the fuel. To aid in this study, numerical sub-models have been developed which will predict the dispersed phase (i.e. fuel) volume fraction, the relative velocity between the two phases, as well as droplet formation.

When using a commercial code to perform multiphase flow simulations, incorporating the formulation of the transport equations via user subroutines into the generic scalar equation is necessary. The continuity equation for one of the phases is converted to a transport equation by via an empirical slip (or relative) velocity relationship. The slip velocity relation involves local flow parameters and a locally determined droplet size, which requires a droplet formation model to predict the size of fuel droplets formed within the mixture layer. Details are provided in Section 1.0 below. This model is explicitly dependent on the dispersed phase volume fraction and relative velocity between the two phases (or components). For further details concerning these other numerical sub-models, the reader is referred to Celik, et al (1998), Celik, et al (1999), Umbel (1998), and Wilson (1999). Friedman et. al. (1999) and Wu et. al. (1999) at Johns Hopkins University performed relevant experiments involving a diesel fuel-water interface. In these studies, global parameters such as the channel height, inlet water and fuel velocities, were used to characterize the flows. These parameters are not suitable for general three-dimensional transient calculations using CFD. The focus of the present study deals with developing a new computational model in light of experiments performed at J.H.U, pertaining directly to the study of liquid-liquid, diesel fuel and water, stratified flow phenomena. This model will be governed by the local flow parameters, i.e. the Weber number, the Reynolds number, and the Richardson number.

The cases we are concerned with in this study involve two incompressible fluid components (both liquids). The equations of motion for the mixture can be written as

Continuity:

                                            (1)

Momentum:

                         (2)

Scalar transport equation:

                                   (3)

where u is the velocity, r is the density, u_s is the slip velocity (see Eq. 8), r is the phase volume fraction, and the subscript a and b refer to the two different components of the mixture. The mixture properties are defined by

                                      (4)

                                      (5)

                                     (6)

                                    (7)

                                    (8)

where n _eff is the effective viscosity, i.e. n _eff=n _laminar +n _turbulent, and Einstein notation is used in Eqs. (1-3, 4 &8).

In deriving the above equations it is assumed that dynamic pressures in both phases, are equal. In addition, some minor assumptions are made in the construction of an effective diffusion term. The details can be found in Umbel (1998). The turbulent eddy viscosity was calculated using the variants of the well-known k-e model (see e.g. Rodi, 1980) with corrections to account for buoyancy production (or damping) and destruction. The variants are used in which case will be mentioned earlier in the text.

Equations (1) and (2) can be solved for the three mixture velocity components and the mixture pressure, P. An additional equation is needed to determine the volume fraction of one of the two components. This equation is derived from the continuity equation as follows. In case of no phase change, the continuity equation for the fuel component, a , is given by

                                         (9)

It is convenient to write the above equation in the form of the general scalar transport equation, Eq. (3) especially, when commercial software is being used. This is done via a variable transformation, f = (r a ra )/r _m, making use, at the same time, of the definitions of the mixture properties given by equations 4-8. Though Eq. (9) is already volume averaged (see Drew, 1983) it is customary to apply time averaging as a means of filtering (or smoothing) the un-resolvable small scales. This results in some additional terms involving higher order correlation between the fluctuation component of volume fraction and the velocity. Here, all these terms are lumped together and included in the diffusion term (Eq. 3) (see Celik and Wilson, 2000, for more details).

The resulting scalar transport equation used in this study is

            (10)

where R = r a /r b = r _f/r _w represents the ratio of the two fluid densities. The volume fraction can then be calculated from

                           (11)

If one takes advantage of the simplifying assumption that the ratio of the densities is nearly one, it can be shown that Eq. (10) reduces to the more familiar form (Halbronn, 1949).

             (12)

Equations 1, 2, and 10 form the basis of our model, which will be referred to as the Single Fluid Scalar Transport (SFST) model.

A sub-model for the SFST model was developed to introduce the effects of slip velocity. Celik et al (1999) modeled the slip velocity with a constitutive equation obtained from Ishii and Zuber (1979) given by

             (13)

where u¥ is the terminal velocity of a discrete particle or droplet. Our assumption of zero slip velocity in the horizontal and tangential directions is a limitation of the droplet formation model. However, the fuel facility that we considered resembles sedimentation tanks and Eq. (13) should provide adequate accuracy. The terminal velocity is calculated from (Umbel, 1998)

                           (14)

where d_p is the droplet diameter and g is the acceleration of gravity.

Figure 1 shows some characteristic droplets forming in the vicinity of the fuel-water interface. It is the size of these droplets that needs to be predicted by the droplet formation model. Parameters based on locally calculable flow properties should be used to capture the droplet behavior and formation. Here we report on the development of such a model.

The new droplet formation model is based on a previous model (Celik and Wilson, 2000) which used a gradient Richardson number, Ri_g defined by

                                             (15)

where U is the horizontal velocity of the fuel in the mixed layer. However, when the gradient Richardson number was used as the local flow parameter, singularities developed (primarily due to the velocity derivative in Eq. (15)) in areas of flow separation and recirculation zones. To prevent the occurrence of these singularities, a turbulent Richardson number expression is formulated using scale analysis on Eq. (15), which leads to

                                    (16)

where l_I is the integral length scale of turbulence and u_rms is root mean square velocity fluctuations. Making the assumptions that the turbulence is isotropic (i.e. u²_rms = (2/3)k), and using l_I ~ k^3/2/e , along with Eq. (5) and (6) to eliminate the mixture density, Eq. (16) can be written as

                                  (17)

where k is the turbulent kinetic energy, e is the rate of dissipation of k, and r _w / r _f is the water to fuel density ratio. It is further assumed that the maximum droplet size formed at the interface cannot be larger than roughly half of the mixed layer thickness, d _m, with this in mind, we suggest the characteristic length scale which determines the droplet size be given by

                                   (18)

where Ri_cr is the critical Richardson number [see Eq. (20)], which is defined as the value of the turbulent Richardson number beyond which no droplet is formed, (i.e. where the flow is no longer turbulent).

An expression for the critical Richardson, Ri_cr, number is derived from a force balance between the inertial and surface tension forces,

                                       (19)

where s is the surface tension. If du/dy @ D U/d _v =u_rms/d _v, and l @ d _v , Eq. (19) yields

                               (20)

In Eq. (20), d _v represents a momentum layer thickness. From the literature (Atsavapranee 1999), an emperical relationship is obtained that relates d _v to _{d m,} which is given by

                                            (21a&b)

The surface tension effects are thus accounted for through the use of the Weber number. However, other mechanisms such as pinch off, which may influence the formulation of these droplets, are not considered. It is further assumed that the droplet size formed is in equilibrium with the surrounding turbulence, hence, it should be limited by a lower bound proportional to the Taylor microscale. The resulting equation is given by.

                               (22 a&b)

The constant c₁ is selected somewhat arbitrarily as 0.5. Experimental observations (Katz, 2000) indicate that it could be as low as 0.1. The constants c₂ and m, present in Eq. (27), are calibrated using only two data points from the shear flow experiments (Wu and Katz, 1999) along with other necessary flow information from the CFX and FLOW-3D calculations. These values c₂ = 0.0365, m = 2.2, and c₃ = 1.0.

The default model for multiphase flows in CFX-4 is the homogeneous multiphase (HMP) model and application of this model was very problematic (see Chang et al, 1996; Chang, 1998, Umbel, 1998) in terms of convergence, accuracy, and computational efficiency. The HMP model showed too much mixing, which also causes the fuel to reach the outlet prematurely. For larger, more complex geometries, the HMP model was also limited in that a turbulence model could not be included because of extreme convergence difficulties and buoyancy effects could not be modeled using the HMP model. Because of these problems, the SFST model was implemented via the user subroutines provided in CFX - 4.

Multiphase flows are dealt with in FLOW-3D (trademark of Flow Sciences, Inc.) also through a variant drift flux approximation, allowing the slip velocity to be calculated from

                         (23)

where Ñ p_i is the pressure gradient. D_f, the particle relaxation time, (referred to as drift flux coefficient in FLOW-3D), is calculated from

                                   (24)

where r₀=D_p /2 is the droplet radius, and is the mixture kinematic viscosity.

Here too, the droplet size D_p is calculated as part of the solution using our droplet formation model DFM-V3.1.

The droplet formation model, DFM - V3.1 was applied to three different cases, which were also investigated experimentally. The first case deals with a two-dimensional shear flow geometry (Figure 2). The impinging jet and the two-compartment tank are the other applications, which are shown in Figures 3 and 4, respectively.

Figure 2 shows a two-dimensional shear flow test case in which the lighter fuel enters from the top left and flows right, and the heavier water enters from the bottom right and flows left. The inlet flow rate of the fuel was a kept at 15 gpm, while the water flow rate was varied.

The goal of this study was to predict the mixed layer thickness as well as the average droplet size. Figure 5 shows the predicted mixed layer thickness for two inlet flow rates, and shows that DFM – V3.1 accurately predicts the mixed layer thickness. In the experiments, the fuel droplets were contained within a 5 cm window. Figure 6 shows the predicted average axial fuel droplet size. It is evident that predicted droplet sizes in good agreement with experiments.

The impinging jet simulation was modeled after the experiments performed by Friedman and Katz (1999). Figure 3 shows the 2D (axisymmetric) model used in the simulations. Here the main interest was in the flow phenomena near the interface, and in droplet formation.

Figure 7 shows contours of the droplet sizes predicted by the DFM-V3.1 for an interface Richardson number of 0.57 and a jet Reynolds number of 12,000. This case corresponds to flow regime 3 as defined by Friedman et. al. 1999. Figure 8 shows the predicted average droplet size, over the mixed layer region, as a function of the normalized radial distance. The average of the all droplets in the region turned out to be 1.21mm, which is very close to the experimental results of 1.24mm.

The half-scale tank model was designed to enable the study of the buoyant jets common during the refueling process of CFBT's. Fuel is introduced through a vertical inlet pipe located in the first compartment. When the fuel fills the top of the first compartment and reaches the manhole it begins to spill into the second compartment forming a buoyant jet. The two flow rates considered were 23.5 gpm and 45 gpm.

The experimental results shown on the left in Figure 9a compare well with the pseudo-steady state solution obtained using the SFST model in CFX (Fig 9b). When Figure 9a is compared to Figure 9c1, it is clear that FLOW-3D has more numerical diffusion in the solution when the first order advection scheme is used. This is not the case when the second order advection scheme is applied, Fig. 9c2, which shows good agreement with experimental data. It should be noted however that the picture in Fig 9a was taken from outside of the tank, thus fuel in the foreground might have blocked the events happening in the center at the top of the tank.

Figure 10 shows the predicted fuel droplet size distribution within the tank using DFM - V3.1. Applying the current droplet formation model, DFM-V3.1 to the half-scale tank using FLOW-3D showed that the fuel droplet sizes ranged from 0.5 mm to 6.0 mm, where as in the experiments they typically ranged from 0.25 mm to 5.5 mm (Atsavapranee and Katz, 2000).

A droplet formation model (DFM-V3.1) has been developed to predict the local fuel droplet size distribution within the mixed layer in flows of two immiscible fluids dominated by shear and buoyancy forces. This model uses the local flow parameters and dimensionless numbers, such as the Weber number, Richardson number, and Reynolds number, to determine the size of the droplets, which is also dependent upon the size of the small eddies present at that location. The needed local flow quantities, such as the turbulent kinetic energy and rate of dissipation are obtained from the solution of a set of transport equations formulated for the mixture itself. A new scalar transport equation is formulated which is linked to the volume fraction of the fuel implicitly. These equations are solved with the help of commercial CFD codes, CFX – 4 and FLOW-3D, utilizing the user subroutines.

The new droplet formation model (DFM-V3.1) was applied to three different test cases, namely the 2-D shear flow, the 2-D axisymmetric impinging jet, and the 3-D two compartment half-scale tank. The predicted results are in good agreement with experimental data for all three cases. Accurate predictions of the mixed layer thickness, as well as the local droplet size distributions, are vital in prediction of the fuel present in the effluent discharge. The new model can in principle be used for tracking of small fuel droplets, which disperse much faster than the bulk fuel itself, and is transported to the outlet with the water. This topic is currently under investigation.

This work was mostly funded by the Office of Naval Research, under Award No. N00014-99-1-0475. The technical monitor at ONR was Dr. E. P. Rood. We thank P. Chang , P. Atsavapranee, Wesley Wilson, J Katz, P. Friedman, X. Wu and R. Schmidt for his administrative help.

[1]  Atsavapranee, P. and Katz, J. (2000) Private Communication.

[2]  Celik, I., Umbel, M., and Wilson, W. (1999) “Computations of turbulent mixing at the interface of a density stratified, shear layer.” Proceedings of the 3^rd ASME/JSME Joint Fluids Engineering Conference, July 18-23, San Francisco, CA.

[3]  Celik, I. and Wilson, W. (2000) “A single fluid mixture model for two-phase liquid-liquid flows.” Proceedings of the ASME 2000 Fluids Engineering Division Summer Meeting, June 11-15, Boston, MA.

[4]  CFX-F3D Flow Solver User Guide, Version 4.3. (1997) Computational Fluid Dynamics Services, Building 8.19, Harwell Laboratory, Oxfordshire OX11, ORA UK, December 1997.

[5]  Chang, P. A., Percival, S., and Hill, B. (1996) “Computations of two-fluid flows through two compartments of a compensated fuel/ballast tank.” Report No.: CRDKNSWC/HD-1370-01, Hydromechanics Directorate Research and Development Report, October 1996.

[6]  Drew, D.A. (1983) “Continuum modeling of two-phase flows.” Theory of dispersed multiphase flow; Proceedings of an Advanced Seminar, University of Wisconsin-Madison, May 26-28, 1982, ed. by Richard E. Meyer. pg. 173.

[7]  FLOW-3D User Manual, Version 7.5. (1999) Flow Sciences, Inc., Los Alamos, NM, USA.

[8]  Friedman, P. D., Katz, J. (1999) "The flow and mixing mechanisms caused by the impingement of an immiscible interface with a vertical jet." Physics of Fluids 11, 2598-2606 (September 1999).

[9]  Halbronn, G., “Remarque sur la theorie de l’Austausch appliquee’ au transprt des materiauz en suspension,” Proceedings, International Association of Hydraulic Research, 3^rd Meeting, Grenoble, France, 1949, 1-6.

[10]  Ishii, M. and Zuber, N. (1979) “Drag coefficident and realtive velocity in bubbly, droplet or particulate flows.” AIChE Journal, 25: 843-855.

[11]  Rodi, W. (1980) “Turbulence models and their application in hydraulics: A state-of-the-art review.” An IAHR (International Association of Hydraulic Research) Publication, Delft, The Netherlands.

[12]  mbel, M. (1998) “Predictions of turbulent mixing at the interface of density stratified, shear flows using CFD.” MS Thesis, West Virginia University.

[13]  Wilson, W. (1999) “The development of a droplet formation and entrainment model for simulations of immiscible liquid-liquid flows.” MS Thesis, West Virginia University.

[14]  Wu, X. and Katz, J. (1999) “On the flow structure and mixing phenomena in a fuel- water stratified shear flow.” Proceedings of the 3^rd ASME/JSME Joint Fluids Engineering Conference, July 18-23, San Francisco, C.A.