Bubble effect on Kelvin-Helmholtz' instability

0.27 MB·

by  Sergey L. Gavrilyuk

Bubble effect on Kelvin-Helmholtz’ instability
8
0
0
S. L. Gavrilyuk∗, H. Gouin∗ and V.M.Teshukov†
2

n
February 2, 2008
a
J

6
1
Abstract

]
n We derive boundary conditions at interfaces (contact discontinuities)
y for a class of Lagrangian models describing, in particular, bubbly flows.
d Weusetheseconditions tostudyKelvin-Helmholtz’instability which de-
- velops in theflow of two superposed layers of a pureincompressible fluid
u
andafluidcontaininggasbubbles,co-flowingwithdifferentvelocities. We
l
f show thatthepresenceofbubblesin onelayerstabilizes the flowin some
.
s intervals of wavelengths.
c
i
s
y 1 Introduction
h
p
Many mathematical models of fluid mechanics are derivedthrough the approx-
[
imation of the solution of boundary value problems for Euler equations. For

1 example, equations for bubbly flows are derived as an approximation of the
v
solution of a complex free-boundary problem describing the motion of the mix-
7
ture of water and gas bubbles (Iordanski (1960), Kogarko(1961), Wijngaarden
5
4 (1968)). Green-Naghdi’s model describing wave motions of a liquid layer of fi-
2 nite depth in a shallow water approximation with account of dispersion effects
1. is another example (Green et al (1974), Green & Naghdi (1976)). Due to the
0 fact that averaging procedures and asymptotic expansions have been used in
8 thederivation,itisnotobvioustodecidewhatboundaryconditionsarenatural
0
for these systems of equations. In this paper, a Hamiltonian formulation of the
:
v problemin(t,x)-spaceisproposedwhichallowsonetofindboundaryconditions
i for a general class of models. This class includes a model of bubbly fluid and
X
dispersiveshallowwater. Weusetheboundaryconditionstostudythestability
r
a of two co-flowing layers of bubbly and pure fluids. In absence of gravity and
capillarity, it is well known that Kelvin-Helmholtz’ instability develops for any
wave lengths in the flow of two superposed layers of pure incompressible fluids
∗LaboratoiredeMod´elisationenM´ecaniqueetThermodynamique,Facult´edesScienceset
TechniquesdeSaint-J´eroˆme,Case322,AvenueEscadrilleNormandie-Niemen,13397Marseille
Cedex20,France,e-mail: [email protected];[email protected]
†LavrentyevInstituteofHydrodynamics,Lavrentyevprospect15,Novosibirsk630090,Rus-
sia, also INRIA, SMASH, 2004 route de Lucioles, 06902 Sophia Antipolis, France, e-mail:
[email protected]
1
(see, for example, Drazin & Reid (1981)). We prove that the presence of bub-
bles in one layer can produce a stabilizing effect: the flow becomes stable with
respect to perturbations of some wave lengths.
2 Governing equations and conditions on mov-
ing interfaces
2.1 Variation of Hamilton’s action
Here we calculate the variation of Hamilton’s action for a special class of La-
grangians. Usually this procedureisusedonlyfor derivationofgoverningequa-
tions. We obtain below not only governing equations but also boundary condi-
tions at inner contact surfaces.
Let us consider the Lagrangianof the form:
∂J
L=L(J, ,z) (2.1)
∂z
t
where z= zi , (i = 0,1,2,3), t is the time, x is the space variable,
x ≡
(cid:18) (cid:19)
ρ (cid:0) (cid:1)
J= , ρ is the fluid density and u is the velocity field. Let us calculate
ρu
(cid:18) (cid:19)
thevariationofHamilton’sactioninthecasewhenthe4-DmomentumJverifies
the equation of continuity:
∂ρ
Div J +div(ρu)=0 (2.2)
≡ ∂t
Hamilton’s action is defined by
a= Ldz (2.3)
Ω
Z
where Ω is a material domain. Considering a smooth one-parameter family of
virtual motions
z=Φ(Z,ε), Φ(Z,0)=ϕ(Z)
(Z stands forthe Lagrangiancoordinates,ε is a smallparameterat the vicinity
of zero and z = ϕ(Z) is the real motion), we define the virtual displacements
ζ(Z) and the Lagrangianvariations δJ(Z) by the formulae:
∂Φ(Z,ε) ∂J(Z,ε)
ζ(Z)= , δJ(Z)= (2.4)
∂ε ∂ε
(cid:12)ε=0 (cid:12)ε=0
(cid:12) (cid:12)
Due to the fact that Z = ϕ−1(z(cid:12)(cid:12)), we consider the variat(cid:12)(cid:12)ions as functions of
Eulerian coordinates and use the same notations ζ(z) and δJ(z) in Eulerian
variables. Hamilton’s principle assumes that ζ(z) = 0 on the boundary ∂Ω of
Ω.
2
In the following, the transposition is denoted by T. For any vectors a and b
we use the notation aTb for their scalar product a b and abT for their tensor
·
producta b. ThedivergenceofthesecondordertensorAisacovectordefined
⊗
by
Div(Ah)=Div(A)h
whereh is anyconstantvectorfield. Div anddiv, Grad and arerespectively
∇
divergence and gradient operators in the 4-D and 3-D space. The identity
matrix and the zero matrix of dimension n are denoted by I and O .
n n
In calculations we shall use the equality
∂ζ
δJ= (Div ζ) I4 J (2.5)
∂z −
(cid:18) (cid:19)
which was proved in Gavrilyuk & Gouin (1999). Variation of the Hamilton
action is
da ∂ζ
δa= = δL+L tr dz (2.6)
dε(cid:12)ε=0 ZΩ(cid:18) (cid:18)∂z(cid:19)(cid:19)
(cid:12)
Since (cid:12)
(cid:12)
∂L ∂L ∂J ∂L
δL= δJ+tr δ + ζ
∂J ∂J ∂z ∂z
∂ (cid:18) (cid:19)
 ∂z 
 (cid:18) (cid:19) 
 
and
∂J ∂δJ ∂J∂ζ
δ =
∂z ∂z − ∂z ∂z
(cid:18) (cid:19)
we get from (2.5) and (2.6):
∂L ∂ζ ∂δJ ∂J∂ζ ∂L
δa= (Div ζ)I4 J+tr AT + ζ+L Div ζ dz
Ω ∂J ∂z − ∂z − ∂z ∂z ∂z
Z (cid:18) (cid:18) (cid:19) (cid:18) (cid:18) (cid:19)(cid:19) (cid:19)
with
AT = ∂L or AT j = ∂L (2.7)
∂J i ∂Ji
∂
∂z (cid:0) (cid:1) ∂zj
(cid:18) (cid:19) (cid:18) (cid:19)
Forthesakeofsimplicitythemeasureofintegrationwillnotbeindicated. Since
for any linear transformation A and vector field v
∂v
Div(Av)=(Div A)v+tr(A )
∂z
we get
∂L ∂J ∂ζ ∂L ∂δJ ∂L
δa= tr J AT + L J Divζ+tr AT + ζ =
Ω ∂J − ∂z ∂z − ∂J ∂z ∂z
Z (cid:18)(cid:18) (cid:19) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)
∂L ∂J ∂ζ ∂L
= tr J AT + L J Div ζ+
Ω ∂J − ∂z ∂z − ∂J
Z (cid:18)(cid:18) (cid:19) (cid:19) (cid:18) (cid:19)
3
∂L
+Div ATδJ Div AT δJ+ ζ ,
− ∂z
and finally (cid:0) (cid:1) (cid:0) (cid:1)
∂L ∂J ∂ζ ∂L
δa= tr J AT + L J Div ζ+Div ATδJ
Ω ∂J − ∂z ∂z − ∂J −
Z (cid:18)(cid:18) (cid:19) (cid:19) (cid:18) (cid:19)
(cid:0) (cid:1)
∂ζ ∂L
Div AT (Div ζ)I4 J+ ζ =
− ∂z − ∂z
(cid:18) (cid:19)
(cid:0) (cid:1)
∂L ∂J ∂ζ
= tr J AT JDiv AT +
Ω ∂J − ∂z − ∂z
Z (cid:18)(cid:18) (cid:19) (cid:19)
(cid:0) (cid:1)
∂L ∂L
+ L J+ DivAT J Div ζ+Div ATδJ + ζ =
− ∂J ∂z
(cid:18) (cid:19)
(cid:0) (cid:1) (cid:0) (cid:1)
∂L ∂J
= Div J AT J Div AT ζ+ATδJ
Ω ∂J − ∂z − −
Z (cid:18)(cid:18) (cid:19) (cid:19)
(cid:0) (cid:1)
∂L ∂J
Div J AT J Div AT ζ+
− ∂J − ∂z −
(cid:18) (cid:19)
(cid:0) (cid:1)
∂L ∂L T ∂L
+Div L J+ DivAT J ζ Grad L J+ DivAT J ζ+ ζ
− ∂J − − ∂J ∂z
(cid:18)(cid:18) (cid:19) (cid:19) (cid:18) (cid:19)
(cid:0) (cid:1) (cid:0) (cid:1)
Let us denote
δL ∂L
Div AT =KT.
δJ ≡ ∂J −
Then (cid:0) (cid:1)
∂L ∂J
δa= Div J KT AT + L KT J I4 ζ+ (2.8)
Ω ∂z − − ∂z −
Z (cid:18) (cid:18) (cid:19)(cid:19)
(cid:0) (cid:1)
∂J
+ NT J KT AT + L KT J I4 ζ+NTATδJ
Z∂Ω (cid:18) − ∂z − (cid:19)
(cid:0) (cid:1)
with NT = ( D ,nT), D denotes the surface velocity and n is the space
n n
−
unit normal vector. Virtual displacements vanish at the boundary ∂Ω and the
surface integral is zero. The volume integral yields the equations of motion in
conservative form as in Gavrilyuk & Gouin (1999).
Inthecasewhenfluidtensorialquantitiesarediscontinuousattheinnerinterface
Σ, expression (2.8) becomes
∂L ∂J
δa= Div J KT AT + L KT J I4 ζ+ (2.9)
Ω ∂z − − ∂z −
Z (cid:18) (cid:18) (cid:19)(cid:19)
(cid:0) (cid:1)
∂J
+ NT J KT AT + L KT J I4 ζ+NTATδJ
Σ − ∂z −
Z (cid:20) (cid:18) (cid:19) (cid:21)
(cid:0) (cid:1)
where the jump through Σ is denoted by [ ].
4
2.2 Governing equations and boundary conditions
We explicit the governing equations and the inner boundary conditions for the
Lagrangian
L= 1ρ u2 W ρ,ρ· , where (·) = d = ∂() +uT∇() (2.10)
2 | | − dt ∂t
(cid:16) (cid:17)
Such a Lagrangian appears in the study of wave propagation in both shallow
water flows with dispersion and bubbly flows (a complete discussion of these
models is given in Gavrilyuk & Teshukov (2001)). We get
∂ρ ∂ρ
∂J ∂t ∂x ∂L u2 ∂W 1∂W 1∂W
= , = | | + ( ρ)T u, uT ( ρ)T
∂z ∂j ∂j ∂J − 2 − ∂ρ ρ ∂ρ· ∇ − ρ ∂ρ· ∇ !
 
 ∂t ∂x 
 
 
Matrix (2.7) becomes
∂W
0T
·
− ∂ρ
∂L  
AT = = .
∂J
∂  ∂W 
∂z  · u O3 
(cid:18) (cid:19)  − ∂ρ 
 
 
We get
∂ ∂W ∂W
Div(AT)= div u , 0T , (2.11)
· ·
−∂t ∂ρ !− ∂ρ ! !
2
u δW 1∂W 1∂W
KT= | | + ( ρ)T u, uT ( ρ)T ,
· ·
− 2 − δρ ρ ∂ρ ∇ − ρ ∂ρ ∇ !
∂W
NT J=ρ nT u D , NTAT = nT u D ,0T .
n · n
− − ∂ρ − !
(cid:0) (cid:1) (cid:0) (cid:1)
We study the case when Σ is a contact surface and, consequently, NT J = 0,
NTAT =0T. We see that the surface integralin (2.9) vanishes if the boundary
conditions :
L KT J [p]=0 (2.12)
− ≡
are fulfilled with p L K(cid:2)T J. Vani(cid:3)shing the volume integral in (2.9) and
≡ −
using relations (2.11), we obtain the governing equations in the form
∂ρ
+div(ρu)=0, (2.13)
∂t
∂ρu
+div ρuuT+pI =0,
∂t
(cid:0) (cid:1)
5
∂W ∂ ∂W ∂W δW
p ρ div u W =ρ W.
· ·
≡ ∂ρ − ∂t ∂ρ !− ∂ρ !!− δρ −
A complete set of boundary conditions at the contact interfaces for the class of
models considered is:
nTu D =0, [p]=0 (2.14)
n
−
Equations (2.13) have been obtained earlier in Gavrilyuk & Shugrin (1996) by
meansofthe classicalmethod ofLagrangemultipliers. Notice thatthis method
did not give jump conditions (2.12). System (2.13) is reminiscent of Euler
equations for compressible fluids. The term p in equations of motion (2.13)
andboundaryconditions(2.14)standsforthe pressure. Nevertheless,pisnota
functionofdensityasinthecaseofbarotropicfluids: itdependsalsoonmaterial
time derivatives of the density. Thus, p is not usual thermodynamic pressure.
In general, it implies that the boundary conditions are not just consequences
of conservation laws. For example, for Korteweg-type fluids where the stress
tensordependsonthedensitygradient,theboundaryconditionscontainnormal
derivatives of the density at interfaces (Seppecher (1989), Gouin & Gavrilyuk
(1999)). As we have shownhere, the boundary conditions at inner interfaces in
bubbly fluids do not depend on normal derivatives of the density: they involve
only tangential derivatives.
If the flow domain is sharedbetween two domains consisting of pure ideal fluid
andbubblyfluid,allpreviouscalculationsarestillvalidandboundaryconditions
(2.14) can be extended to such a case.
2.3 Exact statement of the problem
Now,weproposetouseconditions(2.14)tostudyKelvin-Helmholtz’instability
of a parallel flow between two rigid walls of a layer of bubbly fluid in contact
with a layer of pure incompressible fluid.
Considerthe 2-Dflowoftwosuperposedfinite layersofinviscidfluids inthe
channel: < x < + , 0 < y < H0. The first one is governed by the Euler
−∞ ∞
equations :
div(u1)=0, ∂u1 + uT1 u1+ ∇p1 =0, <x< , 0<y <h(t,x)
∂t ∇ ρ10 −∞ ∞
(cid:0) (cid:1) (2.15)
Hereρ10 isthe constantdensityofpurefluid, u1 =(u1,v1)T isthe velocityfield
and p1 is the pressure.
In the domain h(t,x)<y <H0 governing equations (2.13) for the second fluid
are
∂ρ2
+div(ρ2u2)=0, (2.16)
∂t
∂u2 + uT2 u2+ ∇p2 =0,
∂t ∇ ρ2
(cid:0) (cid:1)
∂W ∂ ∂W ∂W ·
p2 =ρ2 · div · u2 W, W =W(ρ2,ρ2)
∂ρ2 − ∂t ∂ρ2!− ∂ρ2 !!−
6
Here ρ2, u2 = (u2,v2)T and p2 are the average density, velocity and mixture
·
pressure, respectively; W = W(ρ2,ρ2) is a given potential. For the flow of
compressible bubbles having the same radius, potential W has the form (see,
for example, Gavrilyuk (1994) and Gavrilyuk & Teshukov (2001)):
· 3 ·2
W(ρ2,ρ2)=ρ2 cgεg(ρg) 2πnρlR R (2.17)
−
(cid:18) (cid:19)
where ε is the internalenergyof the gasin bubbles, R is the bubble radius,ρ
g g
is the gas density, ρ =const is the density of carrying phase, c =const is the
l g
bubblemassconcentration,nisthenumberofbubblesperunitmass. Potential
W is the difference between the internal energy of gas and the kinetic energy
of fluid due to radial bubble oscillations. The bubble radius and the bubble
density are functions of the average density ρ2 :
−1
4 1 1 1 c 1 1 c
3 g g
πR = − , ρ =c −
g g
3 n(cid:18)ρ2 − ρl (cid:19) (cid:18)ρ2 − ρl (cid:19)
It can be shown that governing equations (2.16) with potential (2.17) coincide
exactly with classical equations of bubbly fluids (Kogarko (1961) and van Wi-
jngaarden (1968)). In particular, the role of ”equation of state”
δW ·
p2 =ρ2 W, W =W(ρ2,ρ2) (2.18)
δρ2 −
plays the Rayleigh-Lamb equation which governs the radial oscillations of a
spherical bubble:
·· 3 ·2 1 dε
2 g
RR+ 2R = ρ (pg−p2), pg =ρgdρ (2.19)
l g
The equivalence between (2.19) and (2.18) with W given by (2.17) has been
proved, for example, in Gavrilyuk (1994). Dissipation-free model (2.16)-(2.17)
assumesthattheslidingbetweencomponentsisnegligible. Moreover,thismodel
is valid only when the volume fraction of bubbles is very small. Notice that a
detailed description of the bubble interaction has been recently done by Russo
&Smereka(1996)andHerrero,Lucquin-Desreux&Perthame(1999)inthecase
of rigidbubbles, and Smereka (2002),Teshukov& Gavrilyuk (2002)in the case
of compressible bubbles, by using a kinetic approach.
At the rigid walls y = 0 and y = H0, the vertical components of the velocity
are equal to zero
v1(t,x,0)=v2(t,x,H0)=0.
In accordance with the previous Section, we prescribe the following boundary
conditions at the contact interface y =h(t,x):
ht+u1hx =v1, ht+u2hx =v2, p1 =p2 (2.20)
7
3 Linear stability problem
3.1 Linearization
Consider the following main parallel flow of two fluids:
u1 =u10 =const, v1 =0, p1 =p0 =const, 0<y <h0, h0 =const
(3.1)
u2 =u20 =const, v2 =0, p2 =p0 =const, ρ2 =ρ20 =const, h0 <y <H0
(3.2)
′
Small perturbations denoted by :
′ ′ ′
u1 =u10+u1, v1 =v1, p1 =p0+p1,
′ ′ ′ ′
ρ2 =ρ20+ρ2, u2 =u20+u2, v2 =v2, p2 =p0+p2
satisfy the linearized system of equations
′ ′ ′ ′ ′ ′
u1x+v1y =0, ρ10D1u1+p1x =0, ρ10D1v1+p1y =0, 0≤y ≤h0 (3.3)
′ ′ ′ ′ ′ ′ ′
D2ρ2+ρ20(u2x+v2y)=0, ρ20D2u2+p2x =0, ρ20D2v +p2y =0, h0 ≤y ≤H0.
(3.4)
We use the notations
∂ ∂
Di = +u0i , i=1,2
∂t ∂x
In (3.4)
′ 2 ′ 2 2 ′
p =a ρ +b D ρ
2 2 2 2
with
∂2W ∂2W
2 2
a =ρ20 ∂ρ2 (ρ20,0), b =−ρ20 ·2 (ρ20,0).
2 ∂ρ
2
Here a is the equilibrium sound velocity, and b is a characteristic wave length
depending on the bubble radius and gas volume fraction. We suppose that
∂2W
· (ρ20,0)=0
∂ρ2∂ρ2
This condition is obviously fulfilled for potential (2.17). The coefficients a and
b calculated in equilibrium R=R0 and p2 =p0 are :
dp 1 1 ρ
2 g 2 l
a = (τ0) , b = (3.5)
− dτ N0ρ20 4πN0R0ρ20
For dilute mixtures we can simplify these expressions by replacing ρ20 by ρl:
dp 1 1 R2
2 g 2 0
a = (τ0) , b = =
− dτ N0ρl 4πN0R0 3α0
8
Here N0 = ρ20n is the number of bubbles per unit volume, R = R0 is the
4
equilibrium bubble radius, α0 =τ0N0 is the volume fraction of gas, τ0 = πR03
3
isthebubblevolume andp =p (τ)isthegaspressureinbubblesexpressedasa
g g
functionofbubblevolume. Theequilibriumsoundspeedaisusuallysmallwith
respect to the gas sound velocity. Expressions (3.5) can be obtained directly
after linearization of Rayleigh-Lamb’s equation for bubbles (2.19).
Boundary conditions (2.20) at y =h0 are
′ ′ ′ ′ ′ ′ ′ ′
ht+u10hx =v1, ht+u20hx =v2, p1 =p2 (3.6)
Let us consider the normal modes of the linear problem (3.3), (3.4), (3.6):
′ ′
u1 =U1(y)exp(i(kx ωt)), v1 =ikV1(y)exp(i(kx ωt)), (3.7)
− −
′ ′
p1 =P1(y)exp(i(kx ωt)), h =Hexp(i(kx ωt))
− −
′ ′
ρ2 =R2(y)exp(i(kx ωt)), u2 =U2(y)exp(i(kx ωt)),
− −
′ ′
v2 =ikV2(y)exp(i(kx ωt)), p2 =P2(y)exp(i(kx ωt))
− −
By substituting into the linearized system we get the system of equations for
unknown amplitudes
U1+V1y =0, ρ10(u10 c)U1+P1 =0, (3.8)
−
2
ρ10k (u10 c)V1+P1y =0,
− −
(u20 c)R2+ρ20(U2+V2y)=0, ρ20(u20 c)U2+P2 =0,
− −
2 2 2 2 2
ρ20k (u20 c)V2+P2y =0, P2 = a b k (u20 c) R2
− − − −
with the following boundary conditions (cid:0) (cid:1)
V1(h0) V2(h0)
V1(0)=0, V2(H0)=0, P1(H0)=P2(H0), = =H
u10 c u20 c
− −
ω
Here c= is the phase velocity.
k
Fromequations(3.8)weobtaintheeigenvalueproblemforpressureamplitudes:
2
P1yy k P1 =0, 0<y <h0, (3.9)
−
2 2
P2yy λ k P2 =0, h0 <y <H0,
−
P1y(h0) P2y(h0)
P1(h0)=P2(h0), ρ10(u10 c)2 = ρ20(u20 c)2, P1y(0)=P2y(H0)=0
− −
where
2 (u20 c)2
λ =1 − .
− a2 b2k2(u20 c)2
− −
9
It follows from (3.9) that the eigenvalues c are solutions of the equation
th(kh0) λ th(λk(H0 h0))
+ − =0 (3.10)
ρ10(u10 c)2 ρ20(u20 c)2
− −
Next we assume that both layers are thin: k(H0 h0) 1, kh0 1, and
− ≪ ≪
dispersion relation (3.10) is simplified into
(H0 h0) 1 1 h0
− + =0 (3.11)
ρ20 (u20 c)2 − a2 b2k2(u20 c)2 ρ10(u10 c)2
(cid:18) − − − (cid:19) −
3.2 Study of dispersion relation
In case u20 =u10, equation (3.11) reduces to the quadratic equation
2 2 2 2 h0ρ20
(1+A)a (1+A)b k +1 (u20 c) =0, with A= .
− − (H0 h0)ρ10
−
(cid:0) (cid:1)
which has only real roots. It means that the flows with equal velocities are
stable.
Letusconsiderthe generalcasewhenu20 =u10.Equation(3.11)canberewrit-
6
ten as a polynomialof fourthdegree. The stability needs that allfour roots are
real. To study this problem, it is convenient to rewrite the dispersion relation
(3.11) in the form
F(z)+A 0 (3.12)
≡
with
d2
2
F(z)=(1 Nz) 1 , (3.13)
− − z2 1
(cid:18) − (cid:19)
ad 1 M u20 u10
z = , d= , N = , M = −
u20 c bk d a
−
Here M is similar to the Machnumber, d is the dimensionless length of pertur-
bationwave. Obviously,z isrealifandonlyifcisreal. The derivativeof(3.13)
is
′ 2
F (z)= 2Nd (1 Nz)ϕ(z) (3.14)
− −
where
1 1 N−1z
ϕ(z)= + − . (3.15)
d2 (z2 1)2
−
Four cases must be considered
◦ 2 −2
1 0<N <1, d N 1 (3.16)
≤ −
◦ 2 −2
2 0<N <1, d >N 1 (3.17)
−
16N2 3N−2
◦ 2 2 2
3 N >1, d > 8N 6 1 +8N 9 (3.18)
27 − r − 4 − !
(cid:0) (cid:1)
10