- Email: [email protected]
The deformation and kinematic mixing of a collapsing interface
Erratum An error occurred in Volume 21, Number 2 of Applied Mathematical Modelling. The first page of “The deformation and kinematic mixing of a collapsing interface,” by Olusegun 1. Ilegbusi, Mahmut Mat, and Malcolm Andrews was transposed with the first page of “Boundaries of safe operation of mechanical oscillators,” by Jesse R. de Souza and Steven R. Bishop. The correct and complete articles now appear on the following pages.
Appl. Math. Modelling 1997, 21:177, March 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
0307-904x/97/$17.00
The deformation and kinematic mixing of a collapsing interface Olusegun
J. Ilegbusi and Mahmut
Department of Mechanical, Boston, MA
Industrial
Mat and Manufacturing
Engineering,
Northeastern
University,
Malcolm Andrews Department
of Mechanical
Engineeting,
Texas A&M
University, College Station, TX
Interface morphology and flow characteristics of buoyancy driven flow inside a closed cavity are investigated numerically with a two-fluid model. A van Leer scheme is used to resolve the sharp property gradient across the interface. Transport equations are solved for zone-averaged variables of each fluid with allowance for interface transfer of momentum. The volume fraction of one of the fluids is used to track the interface. At a small Atwood number the model agrees with the results based on the Boussinesq approximation. Three flow regimes (chaotic, conoective, and diffusive) are formed, depending on the parametric range of the Reynolds number, which is shown to be related to the Grashof number. The mixing eficiency is found to be highest in the chaotic regime and lowest in the diffusive regime. 0 1997 by Elsevier Science Inc. Keywords: buoyancy driven flow, mixing, two-fluids model, Rayleigh-Taylor
1. Introduction Buoyancy-induced mixing inside a cavity is of considerable interest in many materials processing applications, including crystal growth, steel, and alloy making. This paper examines the effects of a Reynolds number defined relative to a Grashof number on the mixing of two fluids in a rectangular cavity. The problem is similar to the early stage of Rayleigh-Taylor mixing, but the initial density jump is here perpendicular to gravity forces. The mixing phenomenon involves a reduction of mixing length scale accomplished by stretching and folding of material lines or surfaces. In Rayleigh-Taylor mixing a heavier fluid initially is placed over a lighter one. Driven by buoyancy, the interface mixing region expands in proportion to gt2 according to Youngs.’ Andrews and Spalding2 conducted a numerical and experimental investigation of Rayleigh-Taylor instability in a tilted cavity. They showed that in such a system a large-scale two-dimensional overturning motion is superimposed on the
Address reprint requests to Dr. Ilegbusi chanical, Industrial, and Manufacturing University, Boston, MA 0211.5. Received 21 March October 1996.
1996; revised
at the Department of MeEngineering, Northeastern
6 September
Appl. Math. Modelling 1997, Vol. 21, February 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
1996; accepted
16
instability
Rayleigh-Taylor mixing process. The experiment showed that the mixing region between the different density fluids contracts in width as the overturning motion dominates the flow. Interface motion and the effect of vortex distribution on the mixing process has been investigated by Aref et al.” In this work two kinds of interfaces, active and passive, were defined by considering their effect on the flow field. Ottino 4 Ottino et al.,‘*” Khakar et a1.,7 and a number of researdhers have studied mixing phenomena theoretically by means of motion of fluid particles from Eulerian and Lagragian view points. Chien, Rising, and Ottino’ have studied laminar and chaotic mixing in several two-dimensional cavity flows by means of material line and blob deformation. They have shown that alternate (periodic) cavity flow is more efficient for mixing than steady flow. At a low Reynolds number, the dimensionless frequency of oscillation is an important parameter to characterize the stretching of the interface. Mixing at relatively large scale has been studied by Turner” and McEwan.i”~” Turner conducted an extensive review of mixing phenomena in the ocean and the atmosphere. McEwan i”.” has explained the formation of internal wave and mixing by a wave-breaking mechanism. In a recent study of mixing in a rectangular cavity Duvali2 has shown that three basic mixing regimes can be identified, depending on the parametric range of Grashof
0307-904x/97/$17.00 PII SO307-904X(96)00147-3
Deformation
and kinematic mixing of a collapsing interface: 0. J. llegbusi et al.
number. These regimes are classified as chaotic (at large Grashof numbers), convective (at intermediate Grashof numbers), and diffusive (at low Grashof numbers). While this study has provided useful insight into flow structure and interface behavior, it was subject to the limitations of the Boussinesq approximation. The objective of this paper is to present an alternative model that is valid beyond the theoretical limit of Boussinesq approximation. To achieve this objective, a novel two-fluid model is employed to calculate zone-averaged velocities and volume fractions of the fluids. This model is based on a proposal of Spalding’” and derives from the theory of two-phase flow. The two fluids are assumed to share occupancy of space but retain their identity throughout the mixing process. This model has been applied to a range of problems including Rayleigh-Taylor mking,Z 14315 plasma synthesis,” material processing in general,” equations are and shear flows. ‘s Transport solved for the zone-averaged quantities, and allowance is made for interface transfer of momentum. The calculated volume fraction (or existence probability of each fluid) coupled with the interface evolution provides a measure of the degree of mixing in the cavity. It is shown in the paper that the three regimes identified by Duval12 can be characterized by a Reynolds number that is shown to be directly related to the Grashof number. The paper is divided into four sections, of which this introduction is the first. In the next section details of mathematical methodology are presented. In Section 3 the results are presented and discussed. Section 4 contains the concluding remarks.
2. Mathematical
Formulation
Consideration is given to two incompressible and miscible fluids of different densities and/or viscosities in a cavity, as shown in Figure 1. The lighter fluid is on the right, and the heavier fluid is on the left. The interface between these fluids can be active or passive, depending on the Reynolds number of the flow field. The system will be assumed to be isothermal such that thermal buoyancy and cross-flux effects are negligible. It is further assumed that no initial perturbation of the interface occurs due to removal of a plate separating the fluids. After the plate is removed a pressure gradient develops as a result of unbal-
T 1 T9
anced body forces at the interface. Three independent dimensionless numbers describe the mixing characteristics of this system. These numbers are: Reynolds number Re, = HpV/p (ratio of inertia to viscous forces), which is shown later to be related to the Grashof number Gr = ApgH3/(( p, + p2)v2) (ratio of buoyancy to viscous forces), Atwood number A, = ( p1 - pz)/( p1 + p2) (ratio of density difference to average density), and aspect ratio Ar = H/L (ratio of height of cavity to its length). Depending on the magnitude of the Reynolds number, the flow field may stretch and deform the interface into an internal wave. As the Reynolds number increases, nonlinear convective terms in the momentum equations become important and the flow field exhibits chaotic behavior. In the diffusive limit (Re + 0 or Gr + 0) these nonlinear terms are negligible, and the flow field cannot deform the interface. To represent mixing in this system, a two-fluid model has been employed. This model is based on the concept of mixing by a combined shear and “sifting” mechanism.13 Sifting occurs when fragments of fluid subjected to larger body forces move through those subjected to smaller body forces in a pressure gradient field. This phenomenon is closely related to the Rayleigh-Taylor mixing.” The two fluids are regarded as two intermingled phases that are separated by sharp, flexible boundaries and interact with each other through the sharing of space and exchange of momentum. Thus at any location there are two flow quantities such as velocity components and volume fractions. The governing transport equations expressing the conservation of mass and momentum for the two-fluid model can be expressed as follows: 2.1 Continuity
at
1.
The initial
positions
+
dX
a(h Pi’i)
=
o
(1)
JY
2.2 U momentum Pi’,)
a(fi
+
d(fj piuf) + a(fi Pi”iui)
’
at
ax
dY
+ABi,x (2) 2.3 V momentum J(fi
H
Figure
+ d(fi Pi"i)
d(fipi>
PiUi)
+ d(fi
dt
=
Pi”iuj)
+ ‘(.fi
Pi”f) dY
dX
-fi$ +F~(u, -
~2) +fiBi,y (3)
of the two fluids.
ApPl.
Math.
Modelling,
1997,
Vol. 21, February
67
Deformation
and kinematic
mixing
of a collapsing
interface:
In equations (l)-(3) i = 1,2 represents the first and second fluid respectively, u and u represent velocity components in the x and y coordinate directions, p is a shared pressure, and I_Lis the viscosity. The second term on the right side of the momentum equations is an interfluid exchange term due to friction at the interface of the fluids, F is interface friction coefficient, and B is the body force. The interface friction coefficient is here expressed as: F, =
-F*=cciPfifz
(4)
where cd is an empirical constant that has been established by Andrews” to be approximately 20.0 for a range of Rayleigh-Taylor Instabilities. The mixture density 5 is calculated from the relation: P
=f,
Pl
ff2
(5)
P2
2.4 Total conservation
(6)
fl +f2 = 1 2.5 Initial and bounda y conditions
The two fluids are at rest at t = 0, and the plate separating them is assumed to be removed instantaneously without perturbation of the interface. Thus the initial conditions may be expressed mathematically as (see Figure I) u1 = u* = 0
(7)
u1 = u2 = 0
(8)
x L/2:
f, = 1.0,f2 = 0.0 f, = 0.0,fi = 1.0
(9)
(10)
The boundary walls are fixed and a no-slip imposed on all velocities such that:
condition
V=o 2.6 Numerical
is
(11) details
A highly accurate van Leer numerical scheme that can handle a sharp density interface is employed to solve the coupled differential equations (l)-(4). Numerical methods with upwind schemes generally suffer from numerical diffusion, which can act as an additional source term to produce unphysical results. l5 A high-order scheme can reduce this problem. However, due to the sharp density
Table 1.
68
Appt.
Math.
Grid independence
results
0. J. llegbusi
et al.
gradient in the present problem, a conventional high-order scheme may cause numerical oscillations. For example the volume fraction for each fluid should lie between 0.0 and 1.0 at every time step. Undershoots or overshoots of these values give unrealistic results and adversely affect the momentum equation since the interface friction term involves a product of volume fractions. The van Leer method is used to prevent such undershoots and overshoots while maintaining high-order numerical accuracy. The superiority of the van Leer over the upwind scheme for such a situation has been extensively dealt with by Andrews.” A summary of this numerical method is given here, and details can be found in Andrews.i5 The momentum equations are solved in two steps such that an intermediate value for velocities is first calculated from convective terms and later updated with Lagrangian source terms. At this stage these updated velocities do not necessarily satisfy the continuity equation. A Poisson equation for pressure is obtained by substituting velocity and pressure correction into the continuity equation. This Poisson equation is then solved using a Gauss-Seidel iteration method that ensures that new velocities and pressures satisfy the continuity and momentum equations. The calculations are carried out until the sum of the absolute mass residuals over all cell faces is less than a prescribed small value of 10p6. This procedure is similar to the SIMPLE method of Patankar and Spalding,” but has been adapted to two-phase flow. A numerically accurate result is obtained with a 40 x 40 grid system in all cases considered. This grid system was selected from a systematic grid-refinement test performed on 20 X 20, 40 X 40, and 60 X 60 grid systems. The results of this test at three points after 30 set are presented in Table 1. It is seen that the values of the volume fraction of fluid 1 at these points do not change significantly beyond the 40 x 40 grid system. The observed trends in the fi values are typical of other locations in the system. A time step of 0.01 set was chosen to ensure the Courant number E = U,t/Ay of less than 0.5. A typical calculation for 110 set of real time in the chaotic regime took approximately 4 hr of cpu time on a Spare 5 Workstation. The van Leer method is used for calculation of the convective terms. For example, the van Leer limiter in the calculation of volume fraction value on the east cell face for the first fluid can be expressed in the form:
fi,, =fi,upwind + skn(~,,.) ’ -rl’e’ Ax&, i 1 (12)
at t= 30 sec.
Grid
20 x 20
30 x 30
40x40
50x50
60X60
f,L/8,3H/16) f&/2. H/2) f,(L/2,3H/4)
0.4631 0.3505 0.2773
0.4421 0.3328 0.2569
0.4378 0.3282 0.2517
0.4355 0.3263 0.2496
0.4340 0.3252 0.2482
Modelling,
1997,
Vol. 21, February
Deformation
and kinematic
where,
(13)
Aw =f;,p -f;,w,
Ae=fCe-f;,,
v,,
El, e = 6 e/l/;lupwmd SV,,,: Volume flux on east cell face of fluid 1 V”1,upwind : Volume of upwind cell occupied by fluid 1 at time n 1 S=
-1 0
i
if Ae and Aw > 0 ifAeandAw
(14)
The IDI value in the van ‘Leer limiter can be chosen as D = 0.0 for first-order accuracy, D = (Ae + Aw)/2Ax for second-order accuracy and,
1 + IE,,,IAw +
2
D=
1+ i
lc:,,lAe 3
2 - I~l,plAe 2 - JE;,~lAw
+
3
of a collapsing
interface:
0. J. llegbusi
I,e
a0
for E~,~ < 0
It is useful to consider this problem from an energy point of view. The difference between the initial potential energy of the system in Figzue 1 and its lowest potential (when the heavy fluid lies horizontally under the light) is given by: APE=-
gH 2L 4
(19)
(P,-Pz)
Assuming no energy losses, the potential energy released is converted into kinetic energy, given by: KE=
1 HL y+p,
+p2w2
(20)
where V is a mean velocity in each of the fluids. Thus the maximum value (with no losses) for the velocity is readily obtained by equating equations (19) and (20) to give:
This velocity provides a useful kinematic problem as:
time scale for the
(15)
for third-order accuracy. The latter scheme employed in the computations presented here.
Two parameters are used to quantify the state of mixing inside the cavity. These parameters are12: Dimensionless length, which is defined as the ratio of the interface length 1 at any time t to the initial length 1,. For convenience, this length has been expressed here as: (16)
The ratio l/T, is the Brunt-VlisalC frequency” that plays a key role in describing oceanographic fluid oscillations. Tk provides an estimate of the slopping time scale in the problem at early times. It is interesting to define a Reynolds number based on the velocity in equation (21) thus: Re, = H~v-~{~=Gr1/2
(23)
Table 2 presents the Gr and equivalent Re, considered here and indicates that the present problems are not turbulent. Hence our need to classify the conditions. The result also implies a fourth regime, “turbulent,” can be defined for Re, > 2000 or Gr > 4 x 106.
from the expressior?:
_/*M(l, t>: n(Z,, t,> dt 0
(17)
where A4 is the symmetric part of the viscous stress tensor and n is the unit vector. Dimensionless width, similarly defined as the ratio of the interface width w, which is the average distance between the 0.05 and 0.95 contours of the volume fraction of one of the fluids, at time t( >O) to a reference width wO, expressed as: w* = (w-w,)
(22)
has been
2.7 Analysis of mixing
and 1 is calculated
et al.
(21)
for E
2
mixing
(181
2.8 Cases considered
Three different flow regimes are considered, as identified by Duval12: chaotic, convective, and diffusive. These regimes are characterized by the Reynolds number defined in equation (23). One case is considered in each regime, details of which are presented in Table 2.
3. Results 3.1 Chaotic regime (Re, = 608)
w0 where w, is here taken to be equal to the separation between the nodes adjacent to the interface at t = 0.
The simulated evolution of the interface and corresponding velocity vectors and streamlines for the chaotic regime (Re, = 608) are presented in Figures 2 to 4. Due to the
Appl. Math.
Modelling,
1997,
Vol. 21, February
69
Deformation
and kinematic Table 2. Case # 1 2 3
mixing
of a collapsing
interface:
Regime
Re,
Gr
Ar
At
Chaotic Convective Diffusive
608 120 0.6
3.7 x 105 1.45 x 104 0.37
1 .o 1.3 1 .o
5x10-5 5x10-5 5x 10-b
.. t=1u set
^ ^^, . ..-7
yrmu=u3mx1u~
1y,,,,+.136x10~
t=20 see
70
Appl.
point is reached. The mixing region contracts in width as the overturning motion dominates the flow (t < 40 set). This width reduction is associated with the elongation of the interface to preserve mass conservation, since the fluids are incompressible and there is very little time for diffusion. Kelvin-Helmhotz instability develops at the late stage of wave formation because of the perturbation introduced by the pressure gradient across the interface. A
^
*“A
I&“=-U.4 IIX
yr&=o.463X
.,.a
Iu
lo4
Evolution of interface and flow pattern in the chaotic regime (Re,=608,
Math.
Modelling,
et al.
Summary of cases considered.
steep density stratification across the interface a horizontal pressure gradient develops that drives the flow as observed in the experiment of Andrews and Spalding’ and the work of Duval.” Initially, this interface is pushed symmetrically from top and bottom of the cavity as shown in Figure 2 at 10 sec. The buoyancy force stretches and deforms the interface to form an internal wave (t = 20 set). This interface continually elongates until a maximum
Figure 2.
0. J. llegbusi
1997,
Vol. 21, February
V-=3.36x10-”
Ar=
1 .O).
Deformation
and kinematic mixing of a collapsing interface: 0. J. llegbusi et al. V,=3.02x1t13
-t=80set
Figure 2.
II
I
I
1&,+0.148x10-’
Continued
strong vortex at the center causes this instability to grow, leading to the formation of tendril structure (t = 30 set), typical of chaotic mixing systems.” This vortex is, however, too weak to produce a spiral-like structure called whorl, which is another common feature of chaotic mixing.’ Figure 2 shows that a single vortex at the center is displaced by two adjacent vortices at about t = 34 sec. These two vortices play an important role in determining
the shape of the folding interface. The interface is rotated clockwise by these vortices, and it is advected at the folding region because of the strength of the vortices. The interface subsequently breaks up (t = 40 set), producing many vortices that serve to stir and homogenize the mixing region. Complete homogenization of the two miscible fluids on a molecular scale can occur only by means of molecular diffusion, which has a substantially long time scale compared to the present macroscopic mixing.
Appl.
Math.
Modelling,
1997,
Vol. 21, February
71
Deformation
and kinematic
mixing
of a collapsing
interface:
et al.
gation can be explained by means of the Liapunov exponent, which expresses the stretching of a material line from an initial state. If l&l represents the length at the initial condition and l&l the length at a later time f, the m can be Liapunov exponent for a given orientation expressed as:
0
OOM 0030
a(m)
= lim fin 1-+p
If the Liapunov
4
Figure 3. ber.
0. J. llegbusi
exponent
j$jj i
(24) ii
is nonzero
and positive, then,
Liapunov exponent as a function of Reynolds num-
IdSI = IdSlexpCat)
Figure 2 shows that two fragments are formed at a later time ( > 30 set) near the top and bottom walls. The frontal area of the fragments grows, while the necks get thinner. The fragments finally separate from the necks due to their weight. This result is similar to droplet behavior in Rayleigh-Taylor instabilities. One part of the fragment moves into the flow field and finally decays. The second part, which is closer to the wall, changes direction and subsequently diffuses into the flow-field until about t = 110 set when a stably stratified configuration forms. This behavior resembles “fossil” turbulence2’ in the sense that the fragment mixes as it falls. It is seen that the flow field has been transformed from multiple vortices to a single vortex in the core region (t > 80 set). The nonlinear flow field positioned the fluids such that the lighter fluid lies above the heavier fluid. The interface region diminishes as the body forces balance the pressure gradient and the flow field losses kinetic energy due to viscous dissipation. The mixing process finally degenerates through viscous diffusion until all density stratification disappears. To confirm that this flow does indeed display chaotic behavior, it should satisfy at least one of the necessary conditions outlined by Ottino.” The most important condition is the presence of a horseshoe function, which is the signature of chaos. Such a function involves stretching and folding of a material line onto itself with a forward and backward transformation. Despite the difficulty of detecting the existence of a horseshoe function in the present system, it can be demonstrated that it satisfies another important condition for chaos, which is the exponential growth of the material interface. Exponential elon-
(25)
and the material line elongates exponentially with time.20 The value obtained for the Liapunov exponent in this regime is 0.0341, confirming that the flow is weakly chaotic. The variation of Liapunov exponent as a function of Reynolds number is shown in Figure 3. It is seen that the Liapunov exponent increases significantly at large Reynolds numbers. Figure 46~) presents the simulated elongation of this system in the regime considered here. In this figure, the characteristic time t” has been scaled with the kinematic time scale T introduced in equation (22). It is seen that the material interface grows exponentially, confirming the chaotic nature of the mixing process in this parametric range of Grashof number. The corresponding width is presented in Figure 4(b). It is seen that the interface width increases slowly at the early stages of flow development and contracts during the wave-breaking period. The mixing region expands rapidly until stable stratification occurs. The interface width subsequently displays a slow and constant increment at later stages. 3.2 Convective regime (Re,
= 120)
Results for the convective regime (Re, = 120) are presented in Figures 5 and 6. It is seen that the initial flow behavior is similar to the chaotic regime with interface stretching and deformation into an internal wave. However, the interface does not break up but oscillates and decays to stable stratification of the two fluids. After wave formation, the flow field transforms from single to multiple waves.
A- Wave FormatIon 5
waveBreaking
C-stabte stratlfmtlcfl
Figure 4.
72
Interface
Appl.
Math.
elongation and width as a function of time in the chaotic regime: (a) elongation; (b) width (Re,,=
Modelling,
1997,
Vol. 21, February
608. Ar=
1 .O).
Deformation t=!m set
and kinematic mixing of a collapsing interface: 0. J. Ilegbusi et al. 0X10-’
r,_ =0.63x I
0
V,=4.18x104
Q
I =0.324x
yr,”
E-0.l4x10-’
t=240 set
t=360 set
-
t=910
Figure 5. Deformation (Re,= 120, Ar= 1.3).
Figure 6.
of interface,
wave formation,
oscillation,
wave breaking,
see
and stable stratification
in the convective regime
Interface elongation and width as a function of time in the convective regime: (a) elongation; (b) width (Re,=
Appl.
Math. Modelling,
120, Arc
1997, Vol. 21, February
1.3).
73
Deformation
and kinematic
mixing
of a collapsing
interface:
Figure 5 shows that the interface elongates until the internal wave formation is completed. At this stage the convective transport of momentum of the flow field is balanced by the gravitational and hydrodynamic pressure forces. The interface length reaches a maximum value at the end of this cycle at t = 180 sec. There still exists a residual vortex located in the central region of the flowfield. It is seen that the material interface decreases in length in the oscillatory regime (t = 240 to 910 set). At the end of this regime, the interface forms a horizontal line at the center of the cavity. Kelvin-Helmholtz instability develops due to the shear stresses at the interface. The growth of this instability is limited by the strength of the flow field. At the end of the oscillatory regime, the flow field again transforms from multiple vortices to a single vortex. The wave subsequently dissipates after a series of oscillatory motions due to viscous damping. However, the perturbation at the interface developed during the oscillatory regime still persists even after the wave has been completely dissipated. The stretching and deformation of the interface coupled with the formation and decay of internal wave effectively mix the fluids, leading to stable stratification with the lighter fluid overlying the heavier fluid, as seen in Figure 5 at t = 910 sec. Subsequently the fluids mix diffusively until a final state of uniform volume fraction field is achieved. Figure 5 shows that plumes, which are not evident in the classical Rayleigh-Taylor instabilities, are formed at the top and bottom walls. These plumes continuously diffuse and dissipate until a stable stratification is reached. It is interesting to note that the independent fragments that were observed at the very early stage of the chaotic regime do not appear in this convective regime until the internal wave has decayed. The diffusion of these fragments at the bottom and top walls precede the decay of the plumes as seen in this figure. At the last stage of the mixing process (t > 510 secl, the length of the interface region remains essentially constant. Thus the mixing width presented in Figure 6 represents the most appropriate characteristic parameter for quantifying the degree of mixing in the system. It is seen that the width does not decrease significantly in the wave formation region. This behavior may be attributed to the relatively mild stretching of the interface in this regime. In the oscillatory region the interface width increases as a result of both molecular diffusion and the rapid decrease of interface length. The wave looses all potential energy due to this oscillation. It should be remarked that the final and
t=30
L Figure
74
7.
Appl.
Interface
Math.
behavior
II
in the diffusive
Modelling,
1997,
regime
Vol. 21,
(Re,=0.6,
February
et al.
initial lengths of material interface are not equal in the specific example presented here because the aspect ratio is greater than 1.0. The Liapunov exponent in this convective regime is 0.0011 (based on equation [24]), which is too small to produce exponential growth in elongation. This result is further confirmed by the relatively low value of Liapunov exponent presented in Figure 3 for the range of Reynolds numbers within the convective regime.
3.3 Di@sive
regime (Re,,
= 0.6)
The results for the lowest Reynolds number considered = 0.6) are presented in Figures 7 and 8. It is seen that due to relatively small nonlinear convective forces in this regime, mixing takes place diffusively. The material interface remains essentially undistorted, as shown in Figure 7. Thus the change in interface width is here used to characterize mixing as presented in Figure 8. As in the previous cases the characteristic time is here scaled with the kinematic time scale. These figures show that diffusive mixing occurs over a much longer time scale compared to the two mixing mechanisms described above. The calculated value of the Liapunov exponent is zero, which indicates that there is no evidence of chaotic behavior in this regime.
(Re,
4. Conclusion
The mechanism of mixing of two fluids in a cavity have been studied with a non-Boussinesq model. This model is based on a two-fluid formulation that involves the solution of transport equations for zone-averaged properties of each fluid with allowance for interface transfer of mass and momentum. Three different flow regimes characterized by the Reynolds number have been considered, namely: chaotic (Re - 6081, convective (Re, = 1201, and diffusive (Re -- 0 .6jH . The majorHfindings of this work may be summarized as follows: 1. The chaotic regime is very efficient for mixing fluids characterized by an active flow field that develops due to unbalanced body forces. The material interface is stretched and deformed to form an internal wave. The material interface elongates exponentially before breaking up, thus promoting the mixing.
t=450
set
0. J. llegbusi
Ar=
set
1 .O).
t=1050
set
Deformation
0.019
and kinematic
f t
0.016
-
0.013
-
F Gr H I L P Re, t T u
‘k 0.0100.007
mixing
-
V
V
0.094
W 0.001
4.3
-3.3
-2.3
Figure
8.
regime
(Re,=0.6,
Interface
width Ar=
as a function
x
Y
In t’ of time
2. In the
convective regime the flow field is relatively weaker than the chaotic regime. Internal waves are formed but rather than break up, they oscillate and subsequently decay to stably stratified mixing. 3. In the diffusive regime nonlinear convection is negligible. The length of material interface remains constant, while mixing occurs diffusively. The two limiting cases of chaotic and diffusive mixing considered in this study have wide application in materials processing. In situations requiring rapid and homogeneous mixing, as in ladle metallurgy, alloying processes, and rheocasting, the chaotic regime would be the preferred mode. On the other hand diffusive mixing would be ideal for crystal growth processes to prevent segregation resulting from excessive convection. The results have also shown that the two-fluid mode1 performs as well as models based on the Boussinesq approximation for the relatively small Atwood number considered in this study. The relative advantage of the two-fluid model lies in its ability to represent systems with high Atwood numbers beyond the theoretical limit of Boussinesq approximation. This will be the subject of a future investigation.
Acknowledgment Partial financial support for this work was provided by NASA/ASEE and the Microgravity Processing group at NASA Lewis through a summer fellowship for one of the authors (O.J.I.). The authors would like to acknowledge the contribution of Dr. W. Duval, who suggested the problem and provided useful advice.
Nomenclature aspect ratio Atwood number A, B body force
interface:
0. J. llegbusi
et al.
volume fraction inter-fluid friction coefficient Grashof number height of the cavity length of the interface length of the cavity static pressure Reynolds number time kinematic time scale x direction velocity y direction velocity generalized velocity mixing width horizontal coordinate vertical coordinate
in the diffusive
1 .O).
Greek Symbols
EL
Ar
of a collapsing
P
ij u V
viscosity fluid density mixture density liapunov exponent kinematic viscosity
References 1. Youngs, D. L. Numerical simulation of turbulent mixing by Rayleigh-Taylor instability. Physica 1984, 12D, 32-44 M. J. and Spalding, D. B. A simple experiment to 2. Andrews, investigate two-dimensional mixing by Rayleigh-Taylor instability. Phys. Fluids A 1990, 2(6), 922-924 G. Vortex dynamics of passive and 3. Aref, H. and Tryggvason, active interfaces. Physica 1984, 12D, 59-70 of mixing with diffusion and reaction in 4. Ottino, J. M. Description terms of the concept of material surfaces. .I. Fluid M. 1982, 114, 83-103 5. Ottino, J. M., Ranz, W. E., and Macosko, C. W. A lamellar model for analysis of liquid-liquid mixing. Chem. Eng. Sci. 1979, 34, 877-890 for 6. Ottino, J. M., Ranz, W. E., and Macosko, C. W. A framework description of mechanical mixing of fluids. AIChE J. 1981,27(4), 565-571 D. V., Rising, H., and Ottino, J. M. A comparative I. Khakhar, computational and experimental study of chaotic mixing of viscous fluids. J. Fluid Mech. 1990, 213, 227-249 mixing and 8. Chien, W. L., Rising, H., and Ottino, J. M. Laminar chaotic mixing in several cavity flows. J. Fluid Mech. 1986, 170, 355-377 J. S., Buoyancy Eflecfs in Fluids. Cambridge University 9. Turner, Press, Cambridge, MA, 1973 A. D. The kinematics of stratified mixing through 10. McEwan, internal wavebreaking. J. Fluid Mech. 1983a, 128, 47-57 11. McEwan, A. D. Internal mixing in stratified fluids. J. Fluid Mech. 1983b, 128, 59-80 study of mixing of two fluids under 12. Duval, W. M. B. Numerical low gravity. NASA TM 105865, 1992 NASA, Longley 13. Spalding, D. B. Two-fluid models of turbulence. Workshop of Theoretical approaches to turbulence, Hampton, VA, 1984 D. L. Modelling turbulent mixing by Rayleigh-Taylor 14. Youngs, instability. Physica D 1989, 37, 270-287 computation of convective transport in 15. Andrews, M. J. Accurate transient two-phase flow. Int. J. Numer. M&h. Fluids 1995,21(3), 205-222
Appl.
Math.
Modelling,
1997,
Vol. 21, February
75
Deformation
and kinematic
mixing
of a collapsing
interface:
16. Ilegbusi, 0. J. Mixing and unmixedness in plasma jets 1: near-field analysis. Proceedings of the Fifth Annual Thermal and Fluids Analysis Workshop, NASA Lewis/Ohio Aerospace Inslitute, August 16-20, 1993, pp. 307-328 17. Ilegbusi, 0. J. The two-fluid model of turbulence and Its application in metals processing. J. Mat. Proc. Manufact. Sci. 1994, 3, 143-157 18. Ilegbusi, 0. J. and Spalding, D. B. Prediction of fluid flow and heat transfer characteristics of turbulent shear flows with a two-
76
Appl.
Math.
Modelling,
1997,
Vol. 21, February
0. J. llegbusi
et al.
fluid model of turbulence. Inr. J. Heat Mass Trans. 1989, 32(4), 767-774 19. Patankar, S. V. and Spalding, D. B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. J. Heat Mass Trans. 1972, 15, 1787-1806 20. Ottino, J. M.. The Kinematics of Mixing: Stretching, Chaos, and Transport, Cambridge University Press, New York, 1989 21. Gibson, C. H. Fossil turbulence and internal waves. American Institute of Physics, Conference Proceedings No. 76, Nonlinear Properties of Internal Waues, 1981, pp. 159-179
Appl. Math. Modelling 1997, 21:177, March 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
0307-904x/97/$17.00
The deformation and kinematic mixing of a collapsing interface Olusegun
J. Ilegbusi and Mahmut
Department of Mechanical, Boston, MA
Industrial
Mat and Manufacturing
Engineering,
Northeastern
University,
Malcolm Andrews Department
of Mechanical
Engineeting,
Texas A&M
University, College Station, TX
Interface morphology and flow characteristics of buoyancy driven flow inside a closed cavity are investigated numerically with a two-fluid model. A van Leer scheme is used to resolve the sharp property gradient across the interface. Transport equations are solved for zone-averaged variables of each fluid with allowance for interface transfer of momentum. The volume fraction of one of the fluids is used to track the interface. At a small Atwood number the model agrees with the results based on the Boussinesq approximation. Three flow regimes (chaotic, conoective, and diffusive) are formed, depending on the parametric range of the Reynolds number, which is shown to be related to the Grashof number. The mixing eficiency is found to be highest in the chaotic regime and lowest in the diffusive regime. 0 1997 by Elsevier Science Inc. Keywords: buoyancy driven flow, mixing, two-fluids model, Rayleigh-Taylor
1. Introduction Buoyancy-induced mixing inside a cavity is of considerable interest in many materials processing applications, including crystal growth, steel, and alloy making. This paper examines the effects of a Reynolds number defined relative to a Grashof number on the mixing of two fluids in a rectangular cavity. The problem is similar to the early stage of Rayleigh-Taylor mixing, but the initial density jump is here perpendicular to gravity forces. The mixing phenomenon involves a reduction of mixing length scale accomplished by stretching and folding of material lines or surfaces. In Rayleigh-Taylor mixing a heavier fluid initially is placed over a lighter one. Driven by buoyancy, the interface mixing region expands in proportion to gt2 according to Youngs.’ Andrews and Spalding2 conducted a numerical and experimental investigation of Rayleigh-Taylor instability in a tilted cavity. They showed that in such a system a large-scale two-dimensional overturning motion is superimposed on the
Address reprint requests to Dr. Ilegbusi chanical, Industrial, and Manufacturing University, Boston, MA 0211.5. Received 21 March October 1996.
1996; revised
at the Department of MeEngineering, Northeastern
6 September
Appl. Math. Modelling 1997, Vol. 21, February 0 1997 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010
1996; accepted
16
instability
Rayleigh-Taylor mixing process. The experiment showed that the mixing region between the different density fluids contracts in width as the overturning motion dominates the flow. Interface motion and the effect of vortex distribution on the mixing process has been investigated by Aref et al.” In this work two kinds of interfaces, active and passive, were defined by considering their effect on the flow field. Ottino 4 Ottino et al.,‘*” Khakar et a1.,7 and a number of researdhers have studied mixing phenomena theoretically by means of motion of fluid particles from Eulerian and Lagragian view points. Chien, Rising, and Ottino’ have studied laminar and chaotic mixing in several two-dimensional cavity flows by means of material line and blob deformation. They have shown that alternate (periodic) cavity flow is more efficient for mixing than steady flow. At a low Reynolds number, the dimensionless frequency of oscillation is an important parameter to characterize the stretching of the interface. Mixing at relatively large scale has been studied by Turner” and McEwan.i”~” Turner conducted an extensive review of mixing phenomena in the ocean and the atmosphere. McEwan i”.” has explained the formation of internal wave and mixing by a wave-breaking mechanism. In a recent study of mixing in a rectangular cavity Duvali2 has shown that three basic mixing regimes can be identified, depending on the parametric range of Grashof
0307-904x/97/$17.00 PII SO307-904X(96)00147-3
Deformation
and kinematic mixing of a collapsing interface: 0. J. llegbusi et al.
number. These regimes are classified as chaotic (at large Grashof numbers), convective (at intermediate Grashof numbers), and diffusive (at low Grashof numbers). While this study has provided useful insight into flow structure and interface behavior, it was subject to the limitations of the Boussinesq approximation. The objective of this paper is to present an alternative model that is valid beyond the theoretical limit of Boussinesq approximation. To achieve this objective, a novel two-fluid model is employed to calculate zone-averaged velocities and volume fractions of the fluids. This model is based on a proposal of Spalding’” and derives from the theory of two-phase flow. The two fluids are assumed to share occupancy of space but retain their identity throughout the mixing process. This model has been applied to a range of problems including Rayleigh-Taylor mking,Z 14315 plasma synthesis,” material processing in general,” equations are and shear flows. ‘s Transport solved for the zone-averaged quantities, and allowance is made for interface transfer of momentum. The calculated volume fraction (or existence probability of each fluid) coupled with the interface evolution provides a measure of the degree of mixing in the cavity. It is shown in the paper that the three regimes identified by Duval12 can be characterized by a Reynolds number that is shown to be directly related to the Grashof number. The paper is divided into four sections, of which this introduction is the first. In the next section details of mathematical methodology are presented. In Section 3 the results are presented and discussed. Section 4 contains the concluding remarks.
2. Mathematical
Formulation
Consideration is given to two incompressible and miscible fluids of different densities and/or viscosities in a cavity, as shown in Figure 1. The lighter fluid is on the right, and the heavier fluid is on the left. The interface between these fluids can be active or passive, depending on the Reynolds number of the flow field. The system will be assumed to be isothermal such that thermal buoyancy and cross-flux effects are negligible. It is further assumed that no initial perturbation of the interface occurs due to removal of a plate separating the fluids. After the plate is removed a pressure gradient develops as a result of unbal-
T 1 T9
anced body forces at the interface. Three independent dimensionless numbers describe the mixing characteristics of this system. These numbers are: Reynolds number Re, = HpV/p (ratio of inertia to viscous forces), which is shown later to be related to the Grashof number Gr = ApgH3/(( p, + p2)v2) (ratio of buoyancy to viscous forces), Atwood number A, = ( p1 - pz)/( p1 + p2) (ratio of density difference to average density), and aspect ratio Ar = H/L (ratio of height of cavity to its length). Depending on the magnitude of the Reynolds number, the flow field may stretch and deform the interface into an internal wave. As the Reynolds number increases, nonlinear convective terms in the momentum equations become important and the flow field exhibits chaotic behavior. In the diffusive limit (Re + 0 or Gr + 0) these nonlinear terms are negligible, and the flow field cannot deform the interface. To represent mixing in this system, a two-fluid model has been employed. This model is based on the concept of mixing by a combined shear and “sifting” mechanism.13 Sifting occurs when fragments of fluid subjected to larger body forces move through those subjected to smaller body forces in a pressure gradient field. This phenomenon is closely related to the Rayleigh-Taylor mixing.” The two fluids are regarded as two intermingled phases that are separated by sharp, flexible boundaries and interact with each other through the sharing of space and exchange of momentum. Thus at any location there are two flow quantities such as velocity components and volume fractions. The governing transport equations expressing the conservation of mass and momentum for the two-fluid model can be expressed as follows: 2.1 Continuity
at
1.
The initial
positions
+
dX
a(h Pi’i)
=
o
(1)
JY
2.2 U momentum Pi’,)
a(fi
+
d(fj piuf) + a(fi Pi”iui)
’
at
ax
dY
+ABi,x (2) 2.3 V momentum J(fi
H
Figure
+ d(fi Pi"i)
d(fipi>
PiUi)
+ d(fi
dt
=
Pi”iuj)
+ ‘(.fi
Pi”f) dY
dX
-fi$ +F~(u, -
~2) +fiBi,y (3)
of the two fluids.
ApPl.
Math.
Modelling,
1997,
Vol. 21, February
67
Deformation
and kinematic
mixing
of a collapsing
interface:
In equations (l)-(3) i = 1,2 represents the first and second fluid respectively, u and u represent velocity components in the x and y coordinate directions, p is a shared pressure, and I_Lis the viscosity. The second term on the right side of the momentum equations is an interfluid exchange term due to friction at the interface of the fluids, F is interface friction coefficient, and B is the body force. The interface friction coefficient is here expressed as: F, =
-F*=cciPfifz
(4)
where cd is an empirical constant that has been established by Andrews” to be approximately 20.0 for a range of Rayleigh-Taylor Instabilities. The mixture density 5 is calculated from the relation: P
=f,
Pl
ff2
(5)
P2
2.4 Total conservation
(6)
fl +f2 = 1 2.5 Initial and bounda y conditions
The two fluids are at rest at t = 0, and the plate separating them is assumed to be removed instantaneously without perturbation of the interface. Thus the initial conditions may be expressed mathematically as (see Figure I) u1 = u* = 0
(7)
u1 = u2 = 0
(8)
x
f, = 1.0,f2 = 0.0 f, = 0.0,fi = 1.0
(9)
(10)
The boundary walls are fixed and a no-slip imposed on all velocities such that:
condition
V=o 2.6 Numerical
is
(11) details
A highly accurate van Leer numerical scheme that can handle a sharp density interface is employed to solve the coupled differential equations (l)-(4). Numerical methods with upwind schemes generally suffer from numerical diffusion, which can act as an additional source term to produce unphysical results. l5 A high-order scheme can reduce this problem. However, due to the sharp density
Table 1.
68
Appt.
Math.
Grid independence
results
0. J. llegbusi
et al.
gradient in the present problem, a conventional high-order scheme may cause numerical oscillations. For example the volume fraction for each fluid should lie between 0.0 and 1.0 at every time step. Undershoots or overshoots of these values give unrealistic results and adversely affect the momentum equation since the interface friction term involves a product of volume fractions. The van Leer method is used to prevent such undershoots and overshoots while maintaining high-order numerical accuracy. The superiority of the van Leer over the upwind scheme for such a situation has been extensively dealt with by Andrews.” A summary of this numerical method is given here, and details can be found in Andrews.i5 The momentum equations are solved in two steps such that an intermediate value for velocities is first calculated from convective terms and later updated with Lagrangian source terms. At this stage these updated velocities do not necessarily satisfy the continuity equation. A Poisson equation for pressure is obtained by substituting velocity and pressure correction into the continuity equation. This Poisson equation is then solved using a Gauss-Seidel iteration method that ensures that new velocities and pressures satisfy the continuity and momentum equations. The calculations are carried out until the sum of the absolute mass residuals over all cell faces is less than a prescribed small value of 10p6. This procedure is similar to the SIMPLE method of Patankar and Spalding,” but has been adapted to two-phase flow. A numerically accurate result is obtained with a 40 x 40 grid system in all cases considered. This grid system was selected from a systematic grid-refinement test performed on 20 X 20, 40 X 40, and 60 X 60 grid systems. The results of this test at three points after 30 set are presented in Table 1. It is seen that the values of the volume fraction of fluid 1 at these points do not change significantly beyond the 40 x 40 grid system. The observed trends in the fi values are typical of other locations in the system. A time step of 0.01 set was chosen to ensure the Courant number E = U,t/Ay of less than 0.5. A typical calculation for 110 set of real time in the chaotic regime took approximately 4 hr of cpu time on a Spare 5 Workstation. The van Leer method is used for calculation of the convective terms. For example, the van Leer limiter in the calculation of volume fraction value on the east cell face for the first fluid can be expressed in the form:
fi,, =fi,upwind + skn(~,,.) ’ -rl’e’ Ax&, i 1 (12)
at t= 30 sec.
Grid
20 x 20
30 x 30
40x40
50x50
60X60
f,L/8,3H/16) f&/2. H/2) f,(L/2,3H/4)
0.4631 0.3505 0.2773
0.4421 0.3328 0.2569
0.4378 0.3282 0.2517
0.4355 0.3263 0.2496
0.4340 0.3252 0.2482
Modelling,
1997,
Vol. 21, February
Deformation
and kinematic
where,
(13)
Aw =f;,p -f;,w,
Ae=fCe-f;,,
v,,
El, e = 6 e/l/;lupwmd SV,,,: Volume flux on east cell face of fluid 1 V”1,upwind : Volume of upwind cell occupied by fluid 1 at time n 1 S=
-1 0
i
if Ae and Aw > 0 ifAeandAw
(14)
The IDI value in the van ‘Leer limiter can be chosen as D = 0.0 for first-order accuracy, D = (Ae + Aw)/2Ax for second-order accuracy and,
1 + IE,,,IAw +
2
D=
1+ i
lc:,,lAe 3
2 - I~l,plAe 2 - JE;,~lAw
+
3
of a collapsing
interface:
0. J. llegbusi
I,e
a0
for E~,~ < 0
It is useful to consider this problem from an energy point of view. The difference between the initial potential energy of the system in Figzue 1 and its lowest potential (when the heavy fluid lies horizontally under the light) is given by: APE=-
gH 2L 4
(19)
(P,-Pz)
Assuming no energy losses, the potential energy released is converted into kinetic energy, given by: KE=
1 HL y+p,
+p2w2
(20)
where V is a mean velocity in each of the fluids. Thus the maximum value (with no losses) for the velocity is readily obtained by equating equations (19) and (20) to give:
This velocity provides a useful kinematic problem as:
time scale for the
(15)
for third-order accuracy. The latter scheme employed in the computations presented here.
Two parameters are used to quantify the state of mixing inside the cavity. These parameters are12: Dimensionless length, which is defined as the ratio of the interface length 1 at any time t to the initial length 1,. For convenience, this length has been expressed here as: (16)
The ratio l/T, is the Brunt-VlisalC frequency” that plays a key role in describing oceanographic fluid oscillations. Tk provides an estimate of the slopping time scale in the problem at early times. It is interesting to define a Reynolds number based on the velocity in equation (21) thus: Re, = H~v-~{~=Gr1/2
(23)
Table 2 presents the Gr and equivalent Re, considered here and indicates that the present problems are not turbulent. Hence our need to classify the conditions. The result also implies a fourth regime, “turbulent,” can be defined for Re, > 2000 or Gr > 4 x 106.
from the expressior?:
_/*M(l, t>: n(Z,, t,> dt 0
(17)
where A4 is the symmetric part of the viscous stress tensor and n is the unit vector. Dimensionless width, similarly defined as the ratio of the interface width w, which is the average distance between the 0.05 and 0.95 contours of the volume fraction of one of the fluids, at time t( >O) to a reference width wO, expressed as: w* = (w-w,)
(22)
has been
2.7 Analysis of mixing
and 1 is calculated
et al.
(21)
for E
2
mixing
(181
2.8 Cases considered
Three different flow regimes are considered, as identified by Duval12: chaotic, convective, and diffusive. These regimes are characterized by the Reynolds number defined in equation (23). One case is considered in each regime, details of which are presented in Table 2.
3. Results 3.1 Chaotic regime (Re, = 608)
w0 where w, is here taken to be equal to the separation between the nodes adjacent to the interface at t = 0.
The simulated evolution of the interface and corresponding velocity vectors and streamlines for the chaotic regime (Re, = 608) are presented in Figures 2 to 4. Due to the
Appl. Math.
Modelling,
1997,
Vol. 21, February
69
Deformation
and kinematic Table 2. Case # 1 2 3
mixing
of a collapsing
interface:
Regime
Re,
Gr
Ar
At
Chaotic Convective Diffusive
608 120 0.6
3.7 x 105 1.45 x 104 0.37
1 .o 1.3 1 .o
5x10-5 5x10-5 5x 10-b
.. t=1u set
^ ^^, . ..-7
yrmu=u3mx1u~
1y,,,,+.136x10~
t=20 see
70
Appl.
point is reached. The mixing region contracts in width as the overturning motion dominates the flow (t < 40 set). This width reduction is associated with the elongation of the interface to preserve mass conservation, since the fluids are incompressible and there is very little time for diffusion. Kelvin-Helmhotz instability develops at the late stage of wave formation because of the perturbation introduced by the pressure gradient across the interface. A
^
*“A
I&“=-U.4 IIX
yr&=o.463X
.,.a
Iu
lo4
Evolution of interface and flow pattern in the chaotic regime (Re,=608,
Math.
Modelling,
et al.
Summary of cases considered.
steep density stratification across the interface a horizontal pressure gradient develops that drives the flow as observed in the experiment of Andrews and Spalding’ and the work of Duval.” Initially, this interface is pushed symmetrically from top and bottom of the cavity as shown in Figure 2 at 10 sec. The buoyancy force stretches and deforms the interface to form an internal wave (t = 20 set). This interface continually elongates until a maximum
Figure 2.
0. J. llegbusi
1997,
Vol. 21, February
V-=3.36x10-”
Ar=
1 .O).
Deformation
and kinematic mixing of a collapsing interface: 0. J. llegbusi et al. V,=3.02x1t13
-t=80set
Figure 2.
II
I
I
1&,+0.148x10-’
Continued
strong vortex at the center causes this instability to grow, leading to the formation of tendril structure (t = 30 set), typical of chaotic mixing systems.” This vortex is, however, too weak to produce a spiral-like structure called whorl, which is another common feature of chaotic mixing.’ Figure 2 shows that a single vortex at the center is displaced by two adjacent vortices at about t = 34 sec. These two vortices play an important role in determining
the shape of the folding interface. The interface is rotated clockwise by these vortices, and it is advected at the folding region because of the strength of the vortices. The interface subsequently breaks up (t = 40 set), producing many vortices that serve to stir and homogenize the mixing region. Complete homogenization of the two miscible fluids on a molecular scale can occur only by means of molecular diffusion, which has a substantially long time scale compared to the present macroscopic mixing.
Appl.
Math.
Modelling,
1997,
Vol. 21, February
71
Deformation
and kinematic
mixing
of a collapsing
interface:
et al.
gation can be explained by means of the Liapunov exponent, which expresses the stretching of a material line from an initial state. If l&l represents the length at the initial condition and l&l the length at a later time f, the m can be Liapunov exponent for a given orientation expressed as:
0
OOM 0030
a(m)
= lim fin 1-+p
If the Liapunov
4
Figure 3. ber.
0. J. llegbusi
exponent
j$jj i
(24) ii
is nonzero
and positive, then,
Liapunov exponent as a function of Reynolds num-
IdSI = IdSlexpCat)
Figure 2 shows that two fragments are formed at a later time ( > 30 set) near the top and bottom walls. The frontal area of the fragments grows, while the necks get thinner. The fragments finally separate from the necks due to their weight. This result is similar to droplet behavior in Rayleigh-Taylor instabilities. One part of the fragment moves into the flow field and finally decays. The second part, which is closer to the wall, changes direction and subsequently diffuses into the flow-field until about t = 110 set when a stably stratified configuration forms. This behavior resembles “fossil” turbulence2’ in the sense that the fragment mixes as it falls. It is seen that the flow field has been transformed from multiple vortices to a single vortex in the core region (t > 80 set). The nonlinear flow field positioned the fluids such that the lighter fluid lies above the heavier fluid. The interface region diminishes as the body forces balance the pressure gradient and the flow field losses kinetic energy due to viscous dissipation. The mixing process finally degenerates through viscous diffusion until all density stratification disappears. To confirm that this flow does indeed display chaotic behavior, it should satisfy at least one of the necessary conditions outlined by Ottino.” The most important condition is the presence of a horseshoe function, which is the signature of chaos. Such a function involves stretching and folding of a material line onto itself with a forward and backward transformation. Despite the difficulty of detecting the existence of a horseshoe function in the present system, it can be demonstrated that it satisfies another important condition for chaos, which is the exponential growth of the material interface. Exponential elon-
(25)
and the material line elongates exponentially with time.20 The value obtained for the Liapunov exponent in this regime is 0.0341, confirming that the flow is weakly chaotic. The variation of Liapunov exponent as a function of Reynolds number is shown in Figure 3. It is seen that the Liapunov exponent increases significantly at large Reynolds numbers. Figure 46~) presents the simulated elongation of this system in the regime considered here. In this figure, the characteristic time t” has been scaled with the kinematic time scale T introduced in equation (22). It is seen that the material interface grows exponentially, confirming the chaotic nature of the mixing process in this parametric range of Grashof number. The corresponding width is presented in Figure 4(b). It is seen that the interface width increases slowly at the early stages of flow development and contracts during the wave-breaking period. The mixing region expands rapidly until stable stratification occurs. The interface width subsequently displays a slow and constant increment at later stages. 3.2 Convective regime (Re,
= 120)
Results for the convective regime (Re, = 120) are presented in Figures 5 and 6. It is seen that the initial flow behavior is similar to the chaotic regime with interface stretching and deformation into an internal wave. However, the interface does not break up but oscillates and decays to stable stratification of the two fluids. After wave formation, the flow field transforms from single to multiple waves.
A- Wave FormatIon 5
waveBreaking
C-stabte stratlfmtlcfl
Figure 4.
72
Interface
Appl.
Math.
elongation and width as a function of time in the chaotic regime: (a) elongation; (b) width (Re,,=
Modelling,
1997,
Vol. 21, February
608. Ar=
1 .O).
Deformation t=!m set
and kinematic mixing of a collapsing interface: 0. J. Ilegbusi et al. 0X10-’
r,_ =0.63x I
0
V,=4.18x104
Q
I =0.324x
yr,”
E-0.l4x10-’
t=240 set
t=360 set
-
t=910
Figure 5. Deformation (Re,= 120, Ar= 1.3).
Figure 6.
of interface,
wave formation,
oscillation,
wave breaking,
see
and stable stratification
in the convective regime
Interface elongation and width as a function of time in the convective regime: (a) elongation; (b) width (Re,=
Appl.
Math. Modelling,
120, Arc
1997, Vol. 21, February
1.3).
73
Deformation
and kinematic
mixing
of a collapsing
interface:
Figure 5 shows that the interface elongates until the internal wave formation is completed. At this stage the convective transport of momentum of the flow field is balanced by the gravitational and hydrodynamic pressure forces. The interface length reaches a maximum value at the end of this cycle at t = 180 sec. There still exists a residual vortex located in the central region of the flowfield. It is seen that the material interface decreases in length in the oscillatory regime (t = 240 to 910 set). At the end of this regime, the interface forms a horizontal line at the center of the cavity. Kelvin-Helmholtz instability develops due to the shear stresses at the interface. The growth of this instability is limited by the strength of the flow field. At the end of the oscillatory regime, the flow field again transforms from multiple vortices to a single vortex. The wave subsequently dissipates after a series of oscillatory motions due to viscous damping. However, the perturbation at the interface developed during the oscillatory regime still persists even after the wave has been completely dissipated. The stretching and deformation of the interface coupled with the formation and decay of internal wave effectively mix the fluids, leading to stable stratification with the lighter fluid overlying the heavier fluid, as seen in Figure 5 at t = 910 sec. Subsequently the fluids mix diffusively until a final state of uniform volume fraction field is achieved. Figure 5 shows that plumes, which are not evident in the classical Rayleigh-Taylor instabilities, are formed at the top and bottom walls. These plumes continuously diffuse and dissipate until a stable stratification is reached. It is interesting to note that the independent fragments that were observed at the very early stage of the chaotic regime do not appear in this convective regime until the internal wave has decayed. The diffusion of these fragments at the bottom and top walls precede the decay of the plumes as seen in this figure. At the last stage of the mixing process (t > 510 secl, the length of the interface region remains essentially constant. Thus the mixing width presented in Figure 6 represents the most appropriate characteristic parameter for quantifying the degree of mixing in the system. It is seen that the width does not decrease significantly in the wave formation region. This behavior may be attributed to the relatively mild stretching of the interface in this regime. In the oscillatory region the interface width increases as a result of both molecular diffusion and the rapid decrease of interface length. The wave looses all potential energy due to this oscillation. It should be remarked that the final and
t=30
L Figure
74
7.
Appl.
Interface
Math.
behavior
II
in the diffusive
Modelling,
1997,
regime
Vol. 21,
(Re,=0.6,
February
et al.
initial lengths of material interface are not equal in the specific example presented here because the aspect ratio is greater than 1.0. The Liapunov exponent in this convective regime is 0.0011 (based on equation [24]), which is too small to produce exponential growth in elongation. This result is further confirmed by the relatively low value of Liapunov exponent presented in Figure 3 for the range of Reynolds numbers within the convective regime.
3.3 Di@sive
regime (Re,,
= 0.6)
The results for the lowest Reynolds number considered = 0.6) are presented in Figures 7 and 8. It is seen that due to relatively small nonlinear convective forces in this regime, mixing takes place diffusively. The material interface remains essentially undistorted, as shown in Figure 7. Thus the change in interface width is here used to characterize mixing as presented in Figure 8. As in the previous cases the characteristic time is here scaled with the kinematic time scale. These figures show that diffusive mixing occurs over a much longer time scale compared to the two mixing mechanisms described above. The calculated value of the Liapunov exponent is zero, which indicates that there is no evidence of chaotic behavior in this regime.
(Re,
4. Conclusion
The mechanism of mixing of two fluids in a cavity have been studied with a non-Boussinesq model. This model is based on a two-fluid formulation that involves the solution of transport equations for zone-averaged properties of each fluid with allowance for interface transfer of mass and momentum. Three different flow regimes characterized by the Reynolds number have been considered, namely: chaotic (Re - 6081, convective (Re, = 1201, and diffusive (Re -- 0 .6jH . The majorHfindings of this work may be summarized as follows: 1. The chaotic regime is very efficient for mixing fluids characterized by an active flow field that develops due to unbalanced body forces. The material interface is stretched and deformed to form an internal wave. The material interface elongates exponentially before breaking up, thus promoting the mixing.
t=450
set
0. J. llegbusi
Ar=
set
1 .O).
t=1050
set
Deformation
0.019
and kinematic
f t
0.016
-
0.013
-
F Gr H I L P Re, t T u
‘k 0.0100.007
mixing
-
V
V
0.094
W 0.001
4.3
-3.3
-2.3
Figure
8.
regime
(Re,=0.6,
Interface
width Ar=
as a function
x
Y
In t’ of time
2. In the
convective regime the flow field is relatively weaker than the chaotic regime. Internal waves are formed but rather than break up, they oscillate and subsequently decay to stably stratified mixing. 3. In the diffusive regime nonlinear convection is negligible. The length of material interface remains constant, while mixing occurs diffusively. The two limiting cases of chaotic and diffusive mixing considered in this study have wide application in materials processing. In situations requiring rapid and homogeneous mixing, as in ladle metallurgy, alloying processes, and rheocasting, the chaotic regime would be the preferred mode. On the other hand diffusive mixing would be ideal for crystal growth processes to prevent segregation resulting from excessive convection. The results have also shown that the two-fluid mode1 performs as well as models based on the Boussinesq approximation for the relatively small Atwood number considered in this study. The relative advantage of the two-fluid model lies in its ability to represent systems with high Atwood numbers beyond the theoretical limit of Boussinesq approximation. This will be the subject of a future investigation.
Acknowledgment Partial financial support for this work was provided by NASA/ASEE and the Microgravity Processing group at NASA Lewis through a summer fellowship for one of the authors (O.J.I.). The authors would like to acknowledge the contribution of Dr. W. Duval, who suggested the problem and provided useful advice.
Nomenclature aspect ratio Atwood number A, B body force
interface:
0. J. llegbusi
et al.
volume fraction inter-fluid friction coefficient Grashof number height of the cavity length of the interface length of the cavity static pressure Reynolds number time kinematic time scale x direction velocity y direction velocity generalized velocity mixing width horizontal coordinate vertical coordinate
in the diffusive
1 .O).
Greek Symbols
EL
Ar
of a collapsing
P
ij u V
viscosity fluid density mixture density liapunov exponent kinematic viscosity
References 1. Youngs, D. L. Numerical simulation of turbulent mixing by Rayleigh-Taylor instability. Physica 1984, 12D, 32-44 M. J. and Spalding, D. B. A simple experiment to 2. Andrews, investigate two-dimensional mixing by Rayleigh-Taylor instability. Phys. Fluids A 1990, 2(6), 922-924 G. Vortex dynamics of passive and 3. Aref, H. and Tryggvason, active interfaces. Physica 1984, 12D, 59-70 of mixing with diffusion and reaction in 4. Ottino, J. M. Description terms of the concept of material surfaces. .I. Fluid M. 1982, 114, 83-103 5. Ottino, J. M., Ranz, W. E., and Macosko, C. W. A lamellar model for analysis of liquid-liquid mixing. Chem. Eng. Sci. 1979, 34, 877-890 for 6. Ottino, J. M., Ranz, W. E., and Macosko, C. W. A framework description of mechanical mixing of fluids. AIChE J. 1981,27(4), 565-571 D. V., Rising, H., and Ottino, J. M. A comparative I. Khakhar, computational and experimental study of chaotic mixing of viscous fluids. J. Fluid Mech. 1990, 213, 227-249 mixing and 8. Chien, W. L., Rising, H., and Ottino, J. M. Laminar chaotic mixing in several cavity flows. J. Fluid Mech. 1986, 170, 355-377 J. S., Buoyancy Eflecfs in Fluids. Cambridge University 9. Turner, Press, Cambridge, MA, 1973 A. D. The kinematics of stratified mixing through 10. McEwan, internal wavebreaking. J. Fluid Mech. 1983a, 128, 47-57 11. McEwan, A. D. Internal mixing in stratified fluids. J. Fluid Mech. 1983b, 128, 59-80 study of mixing of two fluids under 12. Duval, W. M. B. Numerical low gravity. NASA TM 105865, 1992 NASA, Longley 13. Spalding, D. B. Two-fluid models of turbulence. Workshop of Theoretical approaches to turbulence, Hampton, VA, 1984 D. L. Modelling turbulent mixing by Rayleigh-Taylor 14. Youngs, instability. Physica D 1989, 37, 270-287 computation of convective transport in 15. Andrews, M. J. Accurate transient two-phase flow. Int. J. Numer. M&h. Fluids 1995,21(3), 205-222
Appl.
Math.
Modelling,
1997,
Vol. 21, February
75
Deformation
and kinematic
mixing
of a collapsing
interface:
16. Ilegbusi, 0. J. Mixing and unmixedness in plasma jets 1: near-field analysis. Proceedings of the Fifth Annual Thermal and Fluids Analysis Workshop, NASA Lewis/Ohio Aerospace Inslitute, August 16-20, 1993, pp. 307-328 17. Ilegbusi, 0. J. The two-fluid model of turbulence and Its application in metals processing. J. Mat. Proc. Manufact. Sci. 1994, 3, 143-157 18. Ilegbusi, 0. J. and Spalding, D. B. Prediction of fluid flow and heat transfer characteristics of turbulent shear flows with a two-
76
Appl.
Math.
Modelling,
1997,
Vol. 21, February
0. J. llegbusi
et al.
fluid model of turbulence. Inr. J. Heat Mass Trans. 1989, 32(4), 767-774 19. Patankar, S. V. and Spalding, D. B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. J. Heat Mass Trans. 1972, 15, 1787-1806 20. Ottino, J. M.. The Kinematics of Mixing: Stretching, Chaos, and Transport, Cambridge University Press, New York, 1989 21. Gibson, C. H. Fossil turbulence and internal waves. American Institute of Physics, Conference Proceedings No. 76, Nonlinear Properties of Internal Waues, 1981, pp. 159-179