S abili y and con e gence o a comple e model o mass di usion
R.C. Cab ales a,1, F. Guillén-González b,∗,2, J.V. Gu ié ez-San ac eu c,2
aDp o. de Ciencias Básicas, Uni e sidad del Bío-Bío, Facul ad de Ciencias, Campus Fe nando May, Casilla 447, Chillán, Chile
bDp o. E.D.A.N., Uni e si y o Se illa, Ap do. 1160, 41080 Se illa, Spain
cDp o. de Ma emá ica Aplicada I, Uni e si y o Se illa, E. T. S. I. In o má ica, A da. Reina Me cedes, s/n, 41012 Se illa, Spain
a icle in o abs ac
A icle his o y:
Recei ed 15 No embe 2007
Recei ed in e ised o m 1 Feb ua y 2011
Accep ed 6 June 2011
A ailable online 24 Augus 2011
Keywo ds:
Th ee-dimensional Kazhikho –Smagulo
model
Densi y-dependen Na ie –S okes p oblem
Fini e elemen s
S abili y
Con e gence
We p opose a ully disc e e scheme o app oxima ing a h ee-dimensional, s ongly
nonlinea model o mass di usion, also called he comple e Kazhikho –Smagulo model.
TheschemeusesaC0fini e-elemen app oxima ion o all unknowns (densi y, eloci y
and p essu e), e en hough he densi y limi , solu ion o he con inuous p oblem, belongs
o H2. A fi s -o de ime disc e iza ion is used such ha , a each ime s ep, one only needs
o sol e wo decoupled linea p oblems o he disc e e densi y and he eloci y–p essu e,
sepa a ely.
We ex end o he comple e model, some s abili y and con e gence esul s al eady ob ained
by he las wo au ho s o a simplified model whe e λ2- e ms a e no conside ed, λbeing
he mass di usion coefficien . Now, di e en a gumen s mus be in oduced, based mainly
on an induc ion p ocess wi h espec o he ime s ep, ob aining a he same ime he
h ee main p ope ies o he scheme: an app oxima e disc e e maximum p inciple o he
densi y, weak es ima es o he eloci y and s ong ones o he densi y. Fu he mo e, he
con e gence owa ds a weak solu ion o he densi y-dependen Na ie –S okes p oblem is
also ob ained as λ→0 (join ly wi h he space and ime pa ame e s).
Finally, some nume ical compu a ions p o e he p ac ical use ulness o he scheme.
1. In oduc ion
1.1. The model
The Kazhikho –Smagulo equa ions desc ibe he mo ion o a iscous, incomp essible fluid ha ing wo di e en densi ies
o be subjec o a di usion e ec which is modeled by Fick’s law. Assume ha he fluid unde conside a ion fills Ω⊆R3a
bounded domain wi h (sufficien ly egula ) bounda y Γ. The e olu ion o such a fluid is ollowed du ing he ime in e al
[0,T] o 0 <T<∞. We adop he con ec ion ha bold- ace le e s deno e ec o ial elemen s and use he no a ion Q=
Ω×(0,T)and Σ=Γ×(0,T). Then he comple e Kazhikho –Smagulo is gi en in conse a i e o m by:
⎧
⎪
⎪
⎨
⎪
⎪
⎩
(ρu) +∇·(ρu−λ∇ρ)⊗u−λu⊗∇ρ−μu+λ2∇·1
ρ∇ρ⊗∇ρ+∇P=ρ in Q,
∇·u=0in Q,
ρ +∇·(ρu−λ∇ρ)=0inQ,
(1)
*Co esponding au ho .
E-mail add esses: obe ocab [email protected] (R.C. Cab ales), [email p o ec ed] (F. Guillén-González), [email p o ec ed] (J.V. Gu ié ez-San ac eu).
1This wo k was made while he fi s au ho was a Pos doc o al s uden a Dp o. E.D.A.N.
2This au ho ’s wo k was pa ially suppo ed by p ojec MTM2009-12927, Spain.
1162 R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185
whe e ρ:Q→R+is he fluid densi y, u:Q→R3is he incomp essible (a e aged) eloci y field, and P:Q→Ris
he fluid p essu e. Mo eo e , ep esen s olume ex e nal o ces ha a e applied o he fluid, and μ>0 and λ>0
s and o he kinema ic iscosi y and mass di usion coefficien s, espec i ely. The enso ial p oduc ma ix o wo ec o s
a=(ai)n
i=1,b=(bi)n
i=1is deno ed by a⊗bwi h coefficien s (a⊗b)i,j=aibj.
Using in (1) he equali ies
(ρu) +∇·(ρu−λ∇ρ)⊗u=ρu +(ρu−λ∇ρ)·∇u,
and
−λ∇·(u⊗∇ρ)=−λ(u·∇)∇ρ=−λ∇(u·∇ρ)+λ∇·ρ(∇u) ,
one a i es a he ollowing (non-conse a i e) o mula ion:
⎧
⎪
⎪
⎨
⎪
⎪
⎩
ρu +(ρu−λ∇ρ)·∇u−∇·(μ−λρ)(∇u) +μ∇u+λ2∇·1
ρ∇ρ⊗∇ρ+∇p=ρ in Q,
∇·u=0in Q,
ρ +u·∇ρ−λρ=0in Q,
(2)
whe e p=P−λu·∇ρis a new po en ial unc ion.
We comple e he model wi h he bounda y condi ions
u(x, )=0,∂ρ
∂n(x, )=0 o (x, )∈Σ, (3)
and he ini ial condi ions
ρ(0)=ρ0(x), u(0,x)=u0(x) o x∈Ω, (4)
whe e n(x) he ou wa ds uni no mal ec o o Ωa he poin x∈Γ. Th oughou his pape , he ini ial densi y will be
assumed o sa is y:
0<m⩽ρ0(x)⩽M o all x∈Ω. (5)
The densi y-dependen Na ie –S okes p oblem is o mally ob ained om sys em (1) by jus choosing λ=0:
⎧
⎨
⎩
ρu +(u·∇)u−μu+∇p=ρ in Q,
∇·u=0inQ,
ρ +u·∇ρ=0inQ,
(6)
join ly wi h ini ial condi ions (4) and only Di ichle bounda y condi ions o he eloci y. I is impo an o obse e ha
he Neumann bounda y condi ion o he densi y does no now make sense. The Kazhikho –Smagulo equa ions (1) can be
seen as a egula iza ion o he densi y-dependen Na ie –S okes equa ion (6).
1.2. Known esul s
The fi s au ho s who deal wi h he ma hema ical analysis o p oblem (1) in i s simplified e sion, i.e., d opping he λ2-
e m, we e Kazhikho and Smagulo [13]. They p o ed, ia a semi-Gale kin me hod, he exis ence o global-in- ime weak
solu ions, unde he hypo hesis on he coefficien s: λ<2μ/(M−m), and he exis ence o local-in- ime s ong solu ions
(which is global-in- ime o wo dimensions).
The comple e model (1) we e fi s s udied by Bei ão da Veiga [2] who es ablished he exis ence o local-in- ime s ong
solu ions using linea iza ion and a fixed poin a gumen . La e , unde he same me hod, Secchi [15] p o ed he exis ence
o global-in- ime s ong solu ions in he wo-dimensional domains i λ/μis small enough. Mo eo e , he es ablished he
asymp o ic beha io , as λ→0, owa ds global-in- ime s ong solu ions o he densi y-dependen Na ie –S okes p oblem.
Assuming nonnega i e ini ial densi y, Guillén-González [7] p o ed he exis ence o global-in- ime weak solu ions as well as
he asymp o ic beha io , as λ→0, owa ds global-in- ime weak solu ions o he densi y-dependen Na ie –S okes p oblem.
In [8], an i e a i e me hod is used o p o e he exis ence and egula i y o s ong solu ions o (1), ob aining mo eo e ,
some con e gence a es.
Wi h ega d o p oblem (6), exis ence o global-in- ime weak solu ions was p o en by Lions in [14]. The egula i y o
solu ions was ob ained by An on se , Kazhikho , and Monakho [1].
The e a e no many nume ical schemes o app oxima e p oblems (1) and (6). In [9], an uncondi ionally s able, con e gen ,
linea nume ical scheme o he wo-dimensional simplified model (1), wi hou λ2 e ms, is s udied. This scheme consis s o
C0fini e-elemen spa ial app oxima ion combined wi h he linea ized backwa d-in- ime Eule me hod. I is pa icula ly in-
e es ing ha his ime disc e iza ion decouples he compu a ion o he disc e e densi y om he eloci y–p essu e pai . Fo
he h ee-dimensional case, a condi ionally s able, con e gen ully disc e e scheme is s udied in [10]. The main di e ences
R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185 1163
gi en in [10] wi h espec o [9] a e o p o e an app oxima e disc e e maximum p inciple o he disc e e densi y and he
asymp o ic beha io as λ→0 owa ds a weak solu ion o he densi y-dependen Na ie –S okes p oblem (6). In [11] a linea
ully disc e e scheme is analyzed om he poin o iew o e o es ima es. Fi s -o de e o es ima es a e p o en assuming
a compa ibili y condi ion be ween he eloci y and densi y spaces ela ed o hypo hesis (H4) below. Mo eo e , he use o
egula i y which equi es an ex a compa ibili y condi ion a =0 o solu ions o be app oxima ed is a oided.
1.3. Ou line
The es o he pape is di ided as ollows. In Sec ion 2, we gi e he main ideas o he ma hema ical analysis o
p oblem (1). In Sec ion 3, we desc ibe he nume ical scheme and announce he wo main esul s. In Sec ion 4, by means
o an induc ion a gumen , we fi s p o e poin wise es ima es o he densi y, and hen ene gy es ima es o he eloci y
and s ong es ima es o he densi y a e ob ained, using he disc e e Laplacian o he densi y. A e wa ds, we show he
compac ness o he densi y and eloci y in Sec ion 5 and he passage o he limi in Sec ion 6, concluding he p oo o
Theo em 4. In Sec ion 7, we s udy he asymp o ic beha io as he di usion pa ame e goes o ze o, p o ing Theo em 5.
Finally, Sec ion 8 is de o ed o p esen ing some nume ical expe iences.
2. Ma hema ical analysis o he comple e Kazhikho –Smagulo model
We lis he e some s anda d no a ion used h oughou he pape . By Lp(Ω) and Hs(Ω),1⩽p⩽∞,s=0,1,...,we
deno e he classical Lebesgue and Sobole spaces, espec i ely. The usual no m in Lp(Ω) and Hs(Ω) is deno ed by ·Lp(Ω)
and ·Hs(Ω). The no m and inne p oduc in L2(Ω) will be deno ed by |·|and (·,·), espec i ely.
To define he concep o weak solu ion, we in oduce he ollowing space o unc ions:
H=u∈L2(Ω):∇·u=0inΩ,u·n=0onΓ,
V=u∈H1
0(Ω):∇·u=0inΩ,
L2
0(Ω) =p∈L2(Ω):
Ω
p(x)dx=0,
H2
N(Ω) =ρ∈H2(Ω):∂ρ
∂n=0on∂Ω, 
Ω
ρ(x)dx=
Ω
ρ0(x)dx.
I is known by Poinca é’s inequali y ha uH1(Ω) and |∇u|a e equi alen no ms on H1
0(Ω). On he o he hand, H2
N(Ω)
is an affine space, and ∇ρH1(Ω) and |ρ|a e equi alen semi-no ms on H2
N(Ω).
Defini ion 1. Apai (ρ,u)is said o be a weak solu ion o p oblem (1)–(3)–(4) on (0,T)i :
(a) u∈L∞(0,T;H)∩L2(0,T;V),ρ∈L∞(0,T;H1(Ω)) ∩L2(0,T;H2
N(Ω)),0<m⩽ρ(x, )⩽M,∀(x, )∈Q.
((b) ∀φ∈C1([0,T];V)such ha φ(T)=0,
T

0−u,ρφ +(ρu−λ∇ρ)·∇φ+μ(∇u,∇φ) −λρ(∇u) ,∇φd −λ2
T

01
ρ∇ρ⊗∇ρ,∇φd
=
T

0
(ρ ,φ)d +ρ0u0,φ(0).
(c) The equa ion o mass di usion (1)cis sa isfied almos e e ywhe e in Q.
Defini ion 2. Apai (ρ,u)is said o be a weak solu ion o p oblem (6)on (0,T)i i e ifies:
(a) u∈L∞(0,T;L2(Ω)) ∩L2(0,T;V),ρ∈L∞(Q)wi h 0 <m⩽ρ(x, )⩽Ma.e. (x, )∈Q.
(b) Fo all φ∈C1([0,T];V)wi h φ(T)=0,
T

0−ρu,φ
+(u·∇)φ+μ(∇u,∇φ)d =
T

0
(ρ ,φ)d +ρ0u0,φ(0).
(c) Fo all ϕ∈C1([0,T];H1(Ω)) wi h ϕ(T)=0,
1164 R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185
−
T

0
(ρ,ϕ )d −
T

0
(ρu,∇ϕ)d =ρ0,ϕ(0).
He e, we only gi e an ou line o he p oo o he exis ence o weak solu ions o (1) o he eade ’s con enience [15].
Some ideas will be used la e .
Theo em 3. Le u0∈H,ρ0∈H1(Ω) and ∈L2(0,T;L6/5(Ω)).I λ/μis sufficien ly small, hen he e exis s a weak solu ion o
p oblem (1),(3)–(4) on (0,T).
P oo . We p oceed o mally, assuming (ρ,u) o be a sufficien ly egula solu ion o (1), (3)–(4). Then he maximum p inci-
ple applied o equa ion (1)c, oge he wi h (5), leads o
0<m⩽ρ(x, )⩽Ma.e. (x, )∈Q.(7)
Mul iplying (1)cby −λρ, in eg a ing o e Ω, in eg a ing by pa s in he con ec i e e m, and aking in o accoun he
in e pola ion inequali y
∇ρL4(Ω) ⩽CΩρ1/2
L∞(Ω)|ρ|1/2,(8)
one a i es a
λd
d |∇ρ|2+λ2|ρ|2⩽C0|∇u|2.(9)
Adding he momen um sys em (1)aby u o he densi y equa ion (1)cby 1
2u·u,in eg a edo e Ω, he ollowing ene gy
equali y holds:
1
2
d
d 
Ω
ρ|u|2dx+μ|∇u|2=λ
Ω
ρ(∇u) :∇udx+λ21
ρ∇ρ⊗∇ρ,∇u+(ρ ,u). (10)
Reasoning as in he simplified model (see [13]), one has
λ
Ω
ρ(∇u) :∇udx=λ
Ωρ−M+m
2(∇u) :∇udx⩽λM−m
2|∇u|2.
To es ima e he second e m on he igh -hand side o (10), we use he in e pola ion inequali y (8), he poin wise es ima e
o he densi y m⩽ρ⩽Mand Young’s inequali y, ge ing
λ21
ρ∇ρ⊗∇ρ,∇u
⩽ε1μλ2|ρ|2+C2
ε1
λ2
μ|∇u|2,
wi h ε1>0 o be chosen la e on.
The las e m o (10) is easily bounded by
(ρ ,u)⩽μ
2|∇u|2+C 2
L6/5(Ω).
Compiling he abo e es ima es in o (10), we a i e a
d
d 
Ω
ρ|u|2dx+2μ|∇u|2⩽1+C1
λ
μ+C2
ε1
λ2
μ2μ|∇u|2+ε1μλ2|ρ|2+C 2
L6/5(Ω).(11)
Adding up (11) o (9) mul iplied by με2wi h ε2>0 o be chosen la e on, his gi es us
ε2μλd
d |∇ρ|2+d
d |√ρu|2+(ε2−ε1)μλ2|ρ|2+1−C1
λ
μ−C2
ε1
λ2
μ2−ε2C0μ|∇u|2⩽C 2
L6/5(Ω).
Selec ing ε1=ε2/2, and ε2and λ/μsmall enough such ha C1λ
μ+C2
ε1(λ
μ)2+ε2C0⩽1
2,wege
ε2μλd
d |∇ρ|2+d
d |√ρu|2+ε2
2λ2μ|ρ|2+μ
2|∇u|2⩽C 2
L6/5(Ω).
In eg a ing o ∈(0,T)and bounding om below he densi y, we ob ain
R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185 1165
ε2μλ∇ρ( )
2+mu( )
2+

0ε2
2μλ2ρ(s)
2+μ∇u(s)
2ds ⩽C

0
 (s)

2
L6/5(Ω) ds,(12)
On he o he hand, he ollowing “ ac ional in ime es ima e” holds [1]:
T−δ

0u( +δ) −u( )
2d ⩽Cδ1/2∀δ∈(0,T).
This es ima e implies compac ness o he eloci y uin L2(0,T;L2(Ω)) [16]. Then he exis ence o weak solu ions can be
deduced in a s anda d o m [1] by using he Faedo–Gale kin me hod. 2
3. The fini e-elemen app oxima ion
In his sec ion we se ou ou assump ions on he disc e iza ion o he Kazhikho –Smagulo p oblem. Then we s a e ou
main esul s, Theo ems 4 and 5.
3.1. Hypo heses
F om now on, we assume ha Ωis a bounded domain o R3wi h a polyhed al bounda y and ha he e exis s a
amily o iangula ions {Th}h>0o Ωmade up o e ahed a o hexahed a in h ee dimensions, so ha Ω=K∈ThK.Le
Wh⊂H1(Ω),Vh,˜
Vh⊂H1
0(Ω) and Mh,˜
Mh⊂L2
0(Ω) be fini e-elemen s spaces associa ed o densi y, eloci y and p essu e
espec i ely.
Th oughou his wo k we will suppose he ollowing hypo heses:
(H0) Regula i y o he da a:
u0∈V,ρ0∈H2
N(Ω) wi h 0 <m⩽ρ0⩽Min Ωand ∈L2(0,T;L6/5(Ω)).
Assume λ/μsufficien ly small.
(H1) Assume Ωan open, bounded se o R3, whose bounda y is polyhed al and such ha he con inuous dependencies in
H2-no m o he Poisson–Neumann p oblem and in H2×H1-no m o he S okes hold (see (25) o a Poisson–Neumann
p oblem). This is e ified o example i Ωis con ex [5].
(H2) The iangula ion o Ωand he disc e e spaces e i y
• he in e se inequali ies:
|∇ ¯
ρh|⩽Ch−1|¯
ρh|∀
¯
ρh∈Wh,
|∇ ¯
ρh|L3(Ω) ⩽Ch−1/2|∇ ¯
ρh|∀
¯
ρh∈Wh,
¯
ρhL∞(Ω) ⩽Ch−1/2¯
ρhH1(Ω) ∀¯
ρh∈Wh,
∇ ¯
ρhL4(Ω) ⩽Ch−3/4|∇ ¯
ρh|∀
¯
ρh∈Wh,
•and he in e pola ion e o s:
¯
u−˜Jh¯
uH1(Ω) +¯
u−Jh¯
uH1(Ω) ⩽Ch|¯
u|H2(Ω) ∀¯
u∈H2(Ω) ∩H1
0(Ω),
|¯
p−˜
Kh¯
p|+|¯
p−Kh¯
p|⩽Ch|¯
p|H1(Ω) ∀¯
p∈H1(Ω) ∩L2
0(Ω),
|¯
ρ−Ih¯
ρ|+h¯
ρ−Ih¯
ρH1(Ω) ⩽Ch2|¯
ρ|H2(Ω) ∀¯
ρ∈H2(Ω),
¯
ρ−Ih¯
ρW1,3(Ω)∩L∞(Ω) ⩽Ch1/2|¯
ρ|H2(Ω) ∀¯
ρ∈H2(Ω),
¯
ρ−Ih¯
ρW1,4(Ω) ⩽Ch1/4|¯
ρ|H2(Ω) ∀¯
ρ∈H2(Ω),
whe e Jh,˜Jh,Kh,˜
Khand Iha e in e pola ion ope a o s om H2(Ω) ∩H1
0(Ω) in o Vh,H2(Ω) ∩H1
0(Ω) in o ˜
Vh,
H1(Ω) ∩L2
0(Ω) in o Mh,H1(Ω) ∩L2
0(Ω) in o ˜
Mh, and H2(Ω) in o Wh, espec i ely.
(H3) In –sup condi ions. The e exis β>0 and ˜
β>0 (independen o h) such ha , ∀¯
ph∈Mhand ∀¯
qh∈˜
Mh,
¯
phL2
0(Ω) ⩽βsup
¯
uh∈Vh {0}
(¯
ph,∇·¯
uh)
|∇¯
uh|,
¯
qhL2
0(Ω) ⩽˜
βsup
¯
wh∈˜
Vh {0}
(¯
qh,∇·¯
wh)
|∇ ¯
wh|.

1166 R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185
(H4) Compa ibili y condi ion be ween ˜
Mhand Wh:
(Wh·Wh)∩L2
0(Ω) ⊂˜
Mh,
ha is,
∀¯
ρ1
h,¯
ρ2
h∈Wh,¯
ρ1
h¯
ρ2
h−1
|Ω|
Ω
¯
ρ1
h(x)¯
ρ2
h(x)dx∈˜
Mh.
(H5) Compa ibili y condi ion be ween (Mh,˜
Mh):
Mh⊂˜
Mh.
(H6) S abili y p ope ies
|Jhu|⩽C|u|∀u∈L2(Ω),
|∇ Jhu|⩽C|∇u|∀u∈H1
0(Ω),
IhρH1(Ω) ⩽CρH1(Ω) ∀ρ∈H1(Ω).
Fo ins ance, a way o defining he disc e e spaces (Wh,Vh,Mh,˜
Vh,˜
Mh) e i ying (H2)–(H6) is he ollowing. Le {Th}h>0
be a egula , quasi-uni o m amily o iangula ions o Ω,wi hh=maxK∈ThhK(hK=diame e o K), and
Xl
h=xh∈C0(Ω) such ha xh|K∈Pl(K), ∀K∈Th.
Then we define Wh=X1
h. The e a e se e al possibili ies o define (Vh,Mh)[5], by using he Taylo –Hood elemen (P2×P1)
o he mini-elemen (P1+bubble ×P1), o ins ance. Fo he spaces (˜
Vh,˜
Mh)we choose ˜
Vh=X3
h∩H1
0(Ω) and ˜
Mh=
X2
h∩L2
0(Ω).
No e ha i Vh=˜
Vhand Mh=˜
Mha e chosen, we do no need o conside he p ojec ion p oblem (13).
3.2. Main esul s
The aim o his wo k is o p o e he exis ence o a weak solu ion o p oblem (1) in a ully disc e e se ing by using
fini e elemen s. The essen ial pa o such a p oo lies in ob aining ene gy es ima es o he scheme independen o he
disc e e pa ame e s om which one can in e hei weak con e gence. Then a compac ness a gumen p o ides he s ong
con e gence. The fi s wo ks in ha line o he simplified p oblem (1) we e de eloped in [9] and [10]; un o una ely,
a gumen s he e do no ca y o e o he p esen case due o he oublesome e m λ2∇·(1
ρ∇ρ⊗∇ρ)in (1)1.
To cons uc fini e-elemen app oxima ions o weak solu ions o (1), we mus ace wo difficul ies. The fi s hing is ha
he maximum p inciple (7) has o be sa isfied by he fini e-elemen app oxima ion o he densi y. The second hing is how
o es ablish he disc e e e sion o he ene gy es ima e (12).
To o e come hese difficul ies, we p opose he ollowing scheme. Fo simplici y, we choose a uni o m pa i ion o [0,T]
wi h ime s ep k=T/N:( n=nk)n=N
n=0. We hus conside a backwa d Eule ype scheme which is implici wi h espec o
he di usion e ms and semi-implici wi h espec o he con ec i e e ms. The λ2- e m equi es special a en ion since de-
pending on how we in eg a e i in ime we will be able o no o p o e he s abili y o he scheme. Then he app oxima ions
o (1) emain as ollows.
Ini ializa ion: Le (u0
h,ρ0
h)∈Vh×Whbe app oxima ions o (u0,ρ0)as h→0.
Time s ep n+1:Gi en (un
h,pn
h,ρn
h)∈Vh×Mh×Wh.
•Find (wn
h,qn
h)∈˜
Vhט
Mhsuch ha , o each (¯
wh,¯
qh)∈˜
Vhט
Mh,
∇wn
h,∇¯
wh−qn
h,∇·¯
wh=∇un
h,∇¯
wh,
∇·wn
h,¯
qh=0.(13)
•Find ρn+1
h∈Whsuch ha , o each ¯
ρh∈Wh:
ρn+1
h−ρn
h
k,¯
ρh+wn
h·∇ρn+1
h,¯
ρh+λ∇ρn+1
h,∇¯
ρh=0.(14)
R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185 1167
•Find (un+1
h,pn+1
h)∈Vh×Mhsuch ha , o each (¯
uh,¯
ph)∈Vh×Mh:
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎩
ρn
h
un+1
h−un
h
k,¯
uh+1
2ρn+1
h−ρn
h
kun+1
h,¯
uh+cρn+1
hun
h−λ∇ρn+1
h,un+1
h,¯
uh
+aρn+1
h,un+1
h,¯
uh−λ21
ρn+1
h∇ρn+1
h⊗∇ρn
h,∇¯
uh
=ρn+1
h n+1,¯
uh+pn+1
h,∇·¯
uh,
(15)
∇·un+1
h,qh=0,(16)
whe e we ha e used he ollowing sho -hand no a ion:
n+1=1
k
n+1

n
( )d ,
a(ρ,u, )=μ(∇u,∇ )+
Ω
λ˜
M+˜
m
2−ρ(∇u) :∇ dx,
wi h ˜
M>Mand 0 <˜
m<m, and
c(w,u, )=1
2(w·∇)u, −(w·∇) ,u.
We lis he e some esul s ela ed o coe ci i y and con inui y p ope ies o he ilinea o m defined abo e which will
be used la e :
a(ρ,u,u)⩾μ−λ˜
M−˜
m
2|∇u|2i ˜
m⩽ρ⩽˜
M,(17)
a(ρ,u, )⩽CuH1 H1i ρL∞(Ω) ⩽C,
c(w,u,u)=0,(18)
c(w,u, )⩽CwL3uH1 H1.(19)
He e and below, we deno e by C, wi h o wi hou subsc ip , di e en posi i e cons an s, always independen o he
disc e e pa ame e s (k,h)and, e en ually, depending on he di usion pa ame e λ. In his las case, we will deno e he
cons an s by Cλ.
Scheme (13)–(16) has he ollowing main ea u es. A each ime s ep, h ee linea sys ems which need o be sol ed
sepa a ely o compu e (ρn+1
h,un+1
h,pn+1). This would be s a ed wi h wn
has an H1o hogonal p ojec ion o un
hon o he
disc e e ee-di e gence space ela ed o (˜
Vh,˜
Mh)which in u n depends on ˜
Mh( his p ojec ion will gua an ee an app oxi-
ma e disc e e maximum p inciple o he densi y), second ρn+1
has a fini e-elemen app oxima ion o a con ec ion–di usion
equa ion wi h wn
hbeing he con ec i e eloci y, and hi d (un+1
h,pn+1
h)as a mixed fini e-elemen app oxima ion o he
Na ie –S okes-like equa ions.
In [10] scheme (13)–(16) was p oposed wi hou he e m λ2(1
ρn+1
h∇ρn+1
h⊗∇ρn
h,∇¯
uh). The e we ollow he ollowing
s a egy o es ablish s abili y and con e gence. Fi s ly, s a ing om a sligh a ian o he unca ed scheme de eloped in
[9], we ob ain weak es ima es o he eloci y; a e wa ds an app oxima e disc e e maximum p inciple is es ablished o
he densi y; as a consequence o his, he unca ion ope a o is unnecessa y. Secondly, s ong es ima es o he densi y a e
a ained based on he disc e e in e pola ion (23) ins ead o he disc e e e sion o he Gaglia do–Ni enbe g in e pola ion o
wo dimensions used in [9]. Finally, he con e gence is es ablished by compac ness a gumen s.
The inco po a ion o his λ2 e m in oduces new difficul ies, obliga ing o change he s a egy o he main p oo s. Now,
he ela ion wi h a unca ed scheme does no wo k due o he weak es ima es o he eloci y and he s ong ones o he
densi y a e no ob ained in an independen way as seen in he p oo o Theo em 3. We now make an induc ion p ocess wi h
espec o he ime s ep o ob ain, he h ee main p ope ies o he scheme: an app oxima e disc e e maximum p inciple,
and weak and s ong es ima es o eloci y and densi y, espec i ely. Namely, fixed a ime s ep, we ob ain fi s ly poin wise
es ima es o he densi y, and hen weak es ima es o he eloci y and s ong es ima es o he densi y a he same ime
(see Lemma 8 below). We will see ha scheme (13)–(16) is condi ionally s able and con e gen by imposing he cons ain
(S) below (which was al eady used in [10]). We will also ob ain, as in [10], he con e gence owa ds a weak solu ion o
he densi y-dependen Na ie –S okes equa ions, as he di usion pa ame e goes o ze o join ly wi h he space and ime
pa ame e s acco ding o cons ain (S) below. Since his a gumen is a he simila o ha o [10], we only explici he e he
pass o he limi o he λ2- e ms.
1168 R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185
By uh,kand ρh,kwe deno e he piecewise cons an unc ions aking alues uh,k=un
hand ρh,k=ρn
hon ( n−1, n], espec-
i ely (which we will deno e by uh,k,λ,ρh,k,λ when s udied he asymp o ic beha io wi h espec o λ). Thus he ollowing
wo main esul s will be p o ed in his pape .
Theo em 4. Unde he hypo heses (H0)–(H6), and he cons ain
lim
(h,k)→0
h
k=0,(S)
he e exis s a con e gen subsequence o (uh,k,ρh,k)(deno ed in he same way)as (h,k)→0 owa ds a weak solu ion (u,ρ)o
p oblem (1),(3)–(4) in he sense o Defini ion 1.
Theo em 5. Unde he hypo heses o Theo em 4and ex ending (H2) by he addi ional app oxima ion hypo hesis
|¯
ρ−Ph¯
ρ|⩽Ch2/3¯
ρW1,3/2(Ω) ∀¯
ρ∈W1,3/2(Ω) (H2)
(he e Phis he L2-p ojec o on Wh)and changing cons ain (S) by
lim
(λ,h,k)→0
1
λh
k=0,(S)
hen he e exis s a con e gen subsequence o (uh,k,λ,ρh,k,λ),as(h,k,λ) →0, owa ds a weak solu ion (u,ρ)o he densi y-
dependen Na ie –S okes p oblem in he sense o Defini ion 2.
Finally, we show some nume ical compu a ions in o de o p o e he p ac ical use ulness o scheme (13)–(16). In ac ,
we use scheme (13)–(16), wi h λ=0, in o de o compu e some a es o con e gence and he ins abili y o Rayleigh–Taylo
o he densi y-dependen Na ie –S okes p oblem (6), and a sligh adap a ion o scheme (13)–(16) in o de o simula e a
powde -snow a alanche.
4. A p io i es ima es
Since (13)–(16) a e a sequence o h ee squa e linea sys ems, o p o e exis ence and uniqueness i suffices o p o e only
uniqueness, which will be a consequence o he s abili y o scheme (14)–(16) gi en in his sec ion.
4.1. App oxima e disc e e maximum p inciple
The p oo o he ollowing lemma can be ound in [10].
Lemma 6. Fixed n:0⩽n⩽N−1, i he disc e e eloci ies (wl
h)n
l=0sa is y k n
l=0|∇wl
h|2⩽Cd,wi hC
d>0independen o (h,k,λ),
and n, hen he disc e e solu ion ρn+1
ho (14) sa isfies he poin wise es ima es
0<˜
m⩽ρn+1
h⩽˜
MinΩ
p o ided k and h sufficien ly small sa is ying (S).
4.2. Weak es ima es o he eloci y and s ong ones o he densi y
Conside he linea ope a o h:Wh→Whdefined as:
−(hρh,¯
ρh)=(∇ρh,∇¯
ρh)∀¯
ρh∈Wh.(20)
Then he disc e e densi y equa ion (14) can be ew i en as:
ρn+1
h−ρn
h
k,¯
ρh+wn
h·∇ρn+1
h,¯
ρh−λhρn+1
h,¯
ρh=0.(21)
Be o e we p oceed any u he , we need o es ablish he disc e e e sion o he Gaglia do–Ni enbe g in e pola ion and (8).
Lemma 7. The e exis s C =C(Ω) > 0such ha , o any ρh∈Wh,onehas:
∇ρhL3(Ω) ⩽C|∇ρh|1/2|hρh|1/2,(22)
∇ρhL4(Ω) ⩽Ch1/4|hρh|+ρh1/2
L∞(Ω)|hρh|1/2.(23)
R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185 1169
P oo . Es ima e (22) is ob ained in [10]. To p o e (23), one fi s ly ob ains
∇ρhL4(Ω) ⩽Ch1/4|hρh|+
ρ(h)

1/2
L∞(Ω)|hρh|1/2,(24)
whe e ρ(h)∈H2(Ω) is he solu ion o he ollowing ellip ic p oblem
−ρ(h)=−hρhin Ω, ∂ρ(h)
∂n∂Ω =0,
Ω
ρ(h)(x)dx=0.(25)
Indeed (24) is based on he in e se inequali y ∇ρhL4(Ω) ⩽Ch−3/4|∇ρh|, he app oxima ion p ope y ∇ρ−∇IhρL4(Ω) ⩽
Ch1/4ρH2(Ω), and he in e pola ion inequali y (8).
Nex , o ob ain (23) om (24), i is necessa y o eason as in [10] by compa ing ρ(h)wi h ρh, and using he in e se
inequali y ρhL∞(Ω) ⩽Ch−1/2ρhH1and he in e pola ion e o ρ−IhρL∞(Ω) ⩽Ch1/2ρH2(Ω).
Now, we a e in a posi ion o p o e some ecu si e inequali ies o scheme (13)–(16).
Lemma 8. Fixed n:0
⩽n⩽N−1, assume
0<˜
m⩽ρn
h,ρn+1
h⩽˜
MinΩ(26)
and
1
16C2
λ2μkhρn
h
2+1
2μk∇un
h
2⩽Cd,(27)
wi h Cd>0independen o (h,k,λ), and n. Then, p o ided ha h and k a e sufficien ly small sa is ying (S), he e exis s a unique
solu ion (ρn+1
h,un+1
h,pn+1
h)o scheme (13)–(16) which sa isfies:
⎧
⎪
⎨
⎪
⎩ρn+1
hun+1
h
2−ρn
hun
h
2+ρn
hun+1
h−un
h
2+3μ
4k∇un+1
h
2
⩽kC1
 n+1

2
L6/5(Ω) +εμλ2khρn+1
h
2+hρn
h
2,
(28)
λ∇ρn+1
h
2−λ∇ρn
h
2+λ∇ρn+1
h−ρn
h
2+λ2
2khρn+1
h
2⩽C2k∇un
h
2,(29)
whe e C1,C2,ε>0a e cons an s independen o (h,k,λ),andn,wi hεbeing a bi a ily small.
P oo . Taking ¯
uh=2kun+1
hand ¯
ph=pn+1
hin (15)–(16), and using he iden i y (a−b,2a)=|a|2−|b|2+|a−b|2and p op-
e ies (17) and (18) gi es:
⎧
⎨
⎩ρn+1
hun+1
h
2−ρn
hun
h
2+ρn
hun+1
h−un
h
2+μ−λ( ˜
M−˜
m)k∇un+1
h
2
⩽kC1
 n+1

2
L6/5(Ω) +A,
(30)
whe e A=2kλ2
˜
m∇ρn+1
hL4(Ω)∇ρn
hL4(Ω)|∇un+1
h|and C1=C1(μ)>0 a cons an independen o (h,k,λ), and n.
In iew o inequali y (23) and hypo hesis (26), he e m Acan be bounded by
A⩽Ckλ2h1/2hρn+1
hhρn
h+˜
M1/2h1/4hρn+1
hhρn
h
1/2+hρn
hhρn+1
h
1/2
+˜
Mhρn+1
h
1/2hρn
h
1/2∇un+1
h
:= A1+A2+A3+A4.
By hypo hesis (27), |hρn
h|⩽4C1/2
2C1/2
d/(λμ1/2k1/2), hence he e m A1is bounded as ollows:
A1⩽CC
1/2
d
μ1/2h1/2k1/2λhρn+1
h∇un+1
h⩽kC0
ε1
h
kλ2μhρn+1
h
2+ε1μk∇un+1
h
2,
o any ε1>0 and C0=C0(μ). Thus, by s abili y cons ain (S), we a e allowed o ake hand ksmall enough such ha
C0
ε1
h
k⩽ε2,wi hε2>0 being a bi a y, o a i e a
A1⩽ε2λ2μkhρn+1
h
2+ε1μk∇un+1
h
2.
The e m A2has a simila ea men as A1,
1176 R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
−
T

0ˆ
ρh,k,λ,d
d ˜
ηh,kd +λ
T

0
(∇ρh,k,λ,∇ηh,k)d
−
T

0
(ˆ
Hh,k,λρh,k,λ,∇ηh,k)d =ρ0h,η0
h,
(39)
whe e h,k,˜
h,k( espec i ely, ηh,k,˜
ηh,k) a e sui able app oxima ions o a ee-di e gence es unc ions ∈C1([0,T];
C∞
c(Ω)) ( espec i ely, C1([0,T];C∞
c(Ω)))wi h (T)=0( espec i ely,η(T)=0) such ha he sequence { h,k,λ}h,k,λ is
bounded in L∞(0,T;W1,3(Ω) ∩L∞(Ω)).
He e, we only pass o he limi in he λ2- e m, because he con e gence o he es o e ms has been made in [10]. We
mus p o e ha
J:=λ2
T

01
ρh,k,λ ∇ρh,k,λ ⊗∇ˆ
ρh,k,λ,∇ h,k→0as(h,k,λ)→0.
Indeed, since in pa icula h,kis bounded in L2(0,T;W1,3(Ω)), one has
J⩽Cλ1/2λ3
T

0
∇ρh,k,λ4
L3(Ω)1/4λ3
T

0
∇ ˆ
ρh,k,λ4
L3(Ω)1/4T

0
∇ h,k2
L3(Ω)1/2
⩽Cλ1/2→0.
The e o e, he p oo o Theo em 5 is concluded.
Rema k 17. Replacing he semi-implici app oxima ion o he λ2- e m
λ21
ρn+1
h∇ρn+1
h⊗∇ρn
h,∇¯
uh
by he ully explici app oxima ion
λ21
ρn
h∇ρn
h⊗∇ρn
h,∇¯
uh,
one may es ablish he same s abili y, compac ness and con e gence esul s ob ained p e iously in his wo k.
On he con a y, i we conside he ully implici app oxima ion
λ21
ρn+1
h∇ρn+1
h⊗∇ρn+1
h,∇¯
uh,
hen he e m λ2(1
ρn+1
h∇ρn+1
h⊗∇ρn+1
h,∇un+1
h)is es ima ed using he disc e e in e pola ion inequali y (23) by
CΩkλ2
˜
mh1/2hρn+1
h
2+2˜
M1/2h1/4hρn+1
h
3/2+˜
Mhρn+1
h∇un+1
h.
Then, i is no clea how o con ol he e ms
CΩkλ2
˜
mh1/2hρn+1
h
2+2˜
M1/2h1/4hρn+1
h
3/2∇un+1
h
in o de o ob ain he s abili y es ima es.
8. Nume ical esul s
In his sec ion we p esen h ee ype o nume ical esul s ha show ha he scheme p esen ed be o e does p o ide
good app oxima ions o fluids wi h a iable densi y. In ac , we a e going o app oxima e bo h he mass di usion p oblem
and he limi case λ=0 i.e. he densi y-dependen Na ie –S okes p oblem. The scheme was implemen ed by using he ee
so wa e F eeFem++ [12]. In all simula ions, he app oxima ing spaces a e Vh=˜
Vh=P2and Mh=˜
Mh=P1 o he eloci y
and p essu e, espec i ely, hence he p ojec ion s ep (i.e. he fi s s ep o he scheme) is no necessa y. We conside wo
cases o he densi y: ei he Wh=P2o Wh=P1. No e ha , wi h his choice, hypo hesis (H4) is no e ified.

R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185 1177
Fig. 1. E o s in he L∞(L2)and L∞(H1)no ms o he densi y, eloci y and p essu e, by using P1elemen s o he densi y. We conside a fixed mesh size
h=0.0234458.
The fi s nume ical esul is ela ed o a con e gence es ha illus a es ha he scheme o λ=0 leads o op imal
a es o con e gence o smoo h solu ions, and nume ical a es o con e gence wi h espec o λ→0+a e also p esen ed.
The second nume ical esul is conce ned o he so-called Rayleigh–Taylo ins abili y o he densi y-dependen Na ie –
S okes p oblem, wi h Reynolds numbe 5000 and 20 000, whe e has been needed o in oduce a unca ion p ocedu e o
he densi y in o de o p e en nume ical oscilla ions and o p ese e s abili y in calcula ions ( aking in o accoun ha
s abili y cons ain s imposed in he p e ious nume ical analysis (S)is no e ified in ou nume ical compu a ions, because
we choose ko he same o de han ha ending o a compu a ional cos and λ=0). This unca ion consis s in edefining
he densi y ρn+1
hcalcula ed by he con inui y equa ion (14) as ρn+1
h=χ(ρn+1
h), whe e χ ep esen s he cha ac e is ic
unc ion on he in e al [ρmin,ρmax]defined by he physics o he p oblem. The las nume ical compu a ions a e ela ed
o a powde snow a alanche modeled by a mass di usion fluids wi h densi y-dependen iscosi y (see (43) below). In his
case, whe e Reynolds numbe is aken o o de 106, we also unca e he densi y and mo eo e , an adap i e mesh p ocedu e
is implemen ed, based on he esidual o he disc e e densi y equa ion.
8.1. Ra es o con e gence
Le Ωbe he uni disk. We conside he densi y-dependen Na ie –S okes equa ions (6) in Qwi h homogeneous slip
bounda y condi ions and he ex e nal o ce and ini ial alues such ha (ρex,uex,pex)defined by
ρex =2+xcossin( )+ysinsin( ),(40)
uex =−ycos( ), xcos( ),(41)
pex =sin(x)sin(y)sin( ). (42)
is an analy ical exac solu ion (see [6]). We will use i as a es o s udying he e o beha io o he scheme o he
p oblem wi h λ=0, wi h espec o he mesh size hand ime s ep k, and when λ→0+.
8.1.1. S udy o λ=0
Ou aim is o e i y he capabili y o he scheme o sol e nume ically he densi y-dependen Na ie –S okes equa ions.
We show he e o be ween he nume ical solu ions ob ained by he scheme, o λ=0, and he exac solu ion gi en
in (40)–(42) in he L∞(L2)and L∞(H1)no ms. Fi s ly, in Figs. 1 and 2, we measu e he e o in ime o he eloci y,
p essu e, and densi y ( his la e is app oxima ed by P1and P2fini e elemen s). In bo h cases, we conside he mesh size
h=0.0234458 and ake ime s eps ksmalle and smalle . The e o in ime o all unknowns is o o de 1. Quan i a i ely,
he bes app oxima ion is o he eloci y and he wo s is o he densi y in he L∞(L2)no m and o he p essu e in
he L∞(H1)no m. We poin ou ha he e o in ime in he L∞(L2)and L∞(H1)no ms is compa able wi h he e o
measu ed in he L2(L2)and L2(H1)no ms.
Secondly, in Fig. 3, we fix he ime s ep k=10−5and we conside he mesh sizes h1=0.4484541, h2=0.3068574, and
h3=0.2410898. The eloci y e o in space is o o de 2 in bo h no ms, and he p essu e e o is o o de 2 in he L∞(L2)
no m and o o de 1 in he L∞(H1)no m.
1178 R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185
Fig. 2. E o s in he L∞(L2)and L∞(H1)no ms o he densi y, eloci y, and p essu e, by using P2elemen s o he densi y. We conside a fixed mesh size
o h=0.0234458.
Fig. 3. E o s in he L∞(L2)and L∞(H1)no ms o he densi y, eloci y and p essu e, by using P2elemen s o he densi y wi h a fixed ime s ep k=10−5.
The conside ed mesh sizes h1=0.4484541,h2=0.3068574, and h3=0.2410898.
8.1.2. S udy o λ→0+
Now we in es iga e he con e gence beha io o he nume ical solu ions o he comple e Kazhikho –Smagulo sys em
when λ→0+. In his case, we show he e o be ween he nume ical solu ion gi en by he scheme, o λ>0, and he
exac solu ion (40)–(42) in he L∞(L2)and L∞(H1)no ms. These esul s a e p esen ed in Figs. 4 and 5.
The e o wi h espec o λis o o de 1 o he eloci y and p essu e in he L∞(L2)no m as shown in Fig. 4. As can
also be seen in Fig. 4, his e o is quan i a i ely bes o he eloci y ollowed by he p essu e, and wo s o he densi y.
In Fig. 5, we depic he e o wi h espec o λin he L∞(H1)no m being o o de 1/2 o he densi y and o o de 1 o
he eloci y and p essu e.
8.2. The Rayleigh–Taylo ins abili y
As a complemen o he s udy o λ=0 o ou scheme, we use i o sol e he p oblem o he iscous Rayleigh–Taylo
ins abili y. Nume ical esul s o his p oblem has been epo ed by se e al au ho s (see [6] and [3], and he e e ences
he ein). Since he solu ion has symme ies, we show he esul s o Ω=(0,d/2)×(−2d,2d), he hal o he o iginal
domain gi en by (−d/2,d/2)×(−2d,2d). The fluid, subjec o g a i y, is ini ially a es and i s ini ial densi y ρ0is gi en
by
ρ0=ρM+ρm
2+ρM−ρm
2 anhy+ηcos(2πx/d)
0.01d,
R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185 1179
Fig. 4. E o s in he L∞(L2)no m o he densi y, eloci y and p essu e.
wi h ρM>ρm>0 and η>0. The p oblem depends on he A wood numbe A defined as
A =ρM−ρm
ρM+ρm
,
and he Reynolds numbe Re defined as
Re =ρmd3/2G1/2
μ,
whe e Gis he accele a ion o he g a i y.
We ha e conside ed he ollowing bounda y condi ions o he eloci y: no slip on ho izon al bounda ies and slip on
e ical bounda ies. As in [3], he ime is scaled as = d/A G. We selec d=1, η=0.1, ρm=1 and ρM=3. Thus,
A =0.5. We compu e solu ions o mesh size h=0.01, ime s ep k=0.01, Re =5000 and Re =20000 by using P1and P2
elemen s o app oxima e he densi y. The esul s a e p esen ed in Figs. 6, 7, 8, and 9.
Fo P1elemen s app oxima ing he densi y, ou esul s a e in good quali a i e ag eemen wi h hose epo ed in [3] and
[6]. On he o he hand, i we use P2 o he densi y, ou esul s di e o ⩾2 in he uppe igh pa o he in e ace whe e
he densi y changes. In ou case, we obse e ex a olls up o ming. In [3] i is poin ed ou ha o p e en ing pe u ba ions
in he eloci y, i is ecommended o use a iangula mesh wi h al e na e di ec ions. Bu , as we see in plo s o Figs. 6,
whe e we use P1 o he densi y, his is no ou case. Howe e , in hose calcula ions whe e we use a P2 o he densi y (see
Fig. 7), we can see i . We a e no su e wha is he o igin o hese ex a olls up; we conjec ha he unca ion p ocedu e
and he mesh adap a ion (in he case o he a alanches) o P1elemen s in oduce mo e nume ical dissipa ion han o P2.
1180 R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185
Fig. 5. E o s in he L∞(H1)no m o he densi y, eloci y and p essu e.
8.3. Powde snow a alanche
We use ou scheme o simula e powde snow a alanches. Some nume ical simula ions can be seen o example in [4]
and he e e ences he ein. In pa icula , we adop he es p oblem in oduced in [4]. In his case, he domain Ωis he
ec angle (−√3,7)×(0,3/2)minus an obs acle O=(5,5+0.05h0)×(0,0.25h0)(h0>0). Ini ially we conside a wo-fluid
a es wi h densi y gi en by
ρ=ρ+1ω++ρ−1ω−,
whe e ω+=(−√3,0)×(0,1/2),ω−=Ω ω+,1ωis he cha ac e is ic unc ion o he se ω, and ρ+>ρ−>0a e wo
di e en cons an s ( he densi y o he hea y and ligh fluid, espec i ely). The bounda y condi ions o he eloci y a e: no
slip condi ion on he la e al su aces Γl={−√3,7}×(0,3/2)and on he op o he domain Γ =(−√3,7)×{3/2}and slip
on he bo om ∂Ω (Γl∪Γ ).
The iscosi y μdepends on he densi y. In ac , by ollowing [4], i we ha e a wo-fluid flow, we can conside he olume
ac ion o each fluid, deno ed φi,i=1,2. These a iables a e defined in he ollowing way. Le dΩbe an elemen a y fluid
olume su ounding an in e io poin x∈dΩ, filled wi h wo fluids o olume dΩ1and dΩ2, such ha
|dΩ|=|dΩ1|+|dΩ2|.
Then, he olume ac ion o each fluid in he poin xis defined as
R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185 1181
Fig. 6. Time cap ion o he Rayleigh–Taylo ins abili y o A =0.5 (densi y a io 3) and Re =5000. The densi y is app oxima ed by using P1elemen s.
Fig. 7. Time cap ion o he Rayleigh–Taylo ins abili y o A =0.5 (densi y a io 3) and Re =5000. The densi y is app oxima ed by using P2elemen s.
φi(x)=lim
|dΩ|→0
x∈dΩ
|dΩi|
|dΩ|.
We e ain only he hea y fluid olume ac ion, deno ed as φ(and 1 −φwill be he ligh fluid olume ac ion). By using
φ,weha e

1182 R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185
Fig. 8. Time cap ion o he Rayleigh–Taylo ins abili y o A =0.5 (densi y a io 3) and Re =20 000. The densi y is app oxima ed by using P1elemen s.
Fig. 9. Time cap ion o he Rayleigh–Taylo ins abili y o A =0.5 (densi y a io 3) and Re =20 000. The densi y is app oxima ed by using P2elemen s.
ρ=φρ++(1−φ)ρ−,μ=φρ+ν++(1−φ)ρ−ν−,
whe e ν+,ν−a e he co esponding kinema ic iscosi ies o each fluid. In his case, he model, in non-conse a i e o m, is
R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185 1183
Fig. 10. Time cap ion o he densi y in powde snow a alanche o =0.2, 0.4, 0.8, 1.0, 1.2, 1.4, 1.6, 1.7 and 1.8sec, o P1elemen s (le ) and o P2
elemen s ( igh ).
⎧
⎨
⎩
ρ(u +u·∇u−λ∇log(ρ)·∇u)−∇·μ∇u+(∇u) +∇P=ρ in Q,
∇·u=0inQ,
ρ +∇·(ρu−λ∇ρ)=0inQ,
(43)
whe e λ>0 is a mass di usion coefficien , Pis he p essu e, u he eloci y, and =g(sin θ,cos θ) is he ex e nal o ce
(depending on he g a i y gand he angle θo he slope o he domain wi h espec o he ho izon al di ec ion). No e ha
he λ2- e m in (15) is omi ed o sol e he co esponding momen um equa ions (43) wi h densi y-dependen iscosi y.
We conside he go e ning equa ions in hei dimensionless o m as in [4]. This p ocedu e in oduces wo scaling pa-
ame e s: he Reynolds numbe Re and he Schmid numbe Sc defined as
1184 R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185
Fig. 11. Time cap ion o he densi y in powde snow a alanche o =2.0, 2.2 and 2.4sec, o P1elemen s (le ) and o P2elemen s ( igh ).
Table 1
Values o he pa ame e s used o nume ical simula ions (see [4]).
Pa ame e Value
G a i y accele a ion, g9.8 m/s2
Slope, θ32 ◦C
Reynolds numbe , Re 106
Schmid numbe , Sc 0.3
Ini ial heigh , h01m
Obs acle heigh , hs0.25h0m
Obs acle hickness 0.05h0m
Pa ame e s o ligh and hea y fluids Hea y fluid Ligh fluid
Densi y, ρ±4kg/m
31kg/m
3
Kinema ic iscosi y, ν±10−5m2/s10
−5m2/s
Re := h0gh0
ν+,Sc := ν+
λ.
Fo his sys em we conside a li le a ian o he scheme p esen ed in Sec ion 3.2 o ollow he complex dynamics o
he phenomenon. The main ideas o his modified scheme a e:
1. The equa ion o calcula ing he app oxima e densi y ρn+1
his he same,
2. As we saw be o e, we unca e he app oxima e densi y as ρn+1
h=χ(ρn+1
h), whe e χis he cha ac e is ic unc ion on
he in e al [ρ−,ρ+](ρ−=1 and ρ+=4, see Table 1), and we use an adap i e mesh p ocedu e based on he esidual
o he densi y equa ion (see [12] o de ails).
We pu in Table 1 he alues conside ed o he simula ions. The ime s ep is fixed o k=0.005. We p esen he esul s
in Figs. 10 and 11 o p oblem (43) by using P1and P2app oxima ion o he densi y, espec i ely. The esul s o P1a e
quali a i ely simila o hose epo ed in [4]. These di e ences could be again caused by di e en nume ical dissipa ion o
bo h app oxima ions wi h espec o he densi y unca ion and he adap i e p ocedu e.
Re e ences
[1] S.N. An on se , A.V. Kazhikho , V.N. Monakho , Bounda y Value P oblems in Mechanics o Nonhomogeneous Fluids, S udies in Ma hema ical and I s
Applica ions, ol. 22, No h-Holland Publishing Co., Ams e dam, 1990.
[2] H. Bei ão da Veiga, Di usion on iscous fluids, exis ence and asymp o ic p ope ies o solu ions, Ann. Sc. No m. Sup. Pisa 10 (1983) 341–355.
[3] C. Calga o, E. C ausé, Th. Goudon, An hyb id fini e olume-fini e elemen me hod o a iable densi y incomp essible flow, J. Compu . Phys. 227 (2008)
4671–4696.
[4] D. Du ykh, C. Aca y-Robe , D. B esch, Nume ical simula ion o powde -snow a alanche in e ac ion wi h an obs acle, Applied Ma hema ical Modelling,
elec onically a ailable a h p://hal.a chi es-ou e es. /docs/00/35/88/80/PDF/A alSim_Du ykh-AR-B esch.pd .
[5] V. Gi aul , P.A. Ra ia , Fini e Elemen Me hods o Na ie –S okes Equa ions: Theo y and Algo i hms, Sp inge -Ve lag, Be lin, 1986.
[6] J.-L. Gue mond, L. Qua apelle, A p ojec ion FEM o a iable densi y incomp essible flows, J. Compu . Phys. 165 (1) (2000) 167–188.
[7] F. Guillén-González, Sob e un modelo asin ó ico de di usión de masa pa a fluidos incomp esibles, iscoso y no homogéneos, in: P oceedings o he
Thi d Ca alan Days on Applied Ma hema ics, ISBN 84-87029-87-6, 1996, pp. 103–114.
[8] F. Guillén-González, P. Damázio, M.A. Rojas-Meda , App oach o egula solu ions o incomp essible fluids wi h mass di usion by an in e a i e me hod,
J. Ma h. Anal. Appl. 326 (1) (2007) 468–487.
[9] F. Guillén-González, J.V. Gu ié ez-San ac eu, Uncondi ional s abili y and con e gence o a ully disc e e scheme o 2D iscous fluids models wi h
mass di usion, Ma h. Comp. 77 (263) (2008) 1495–1524.
R.C. Cab ales e al. / Applied Nume ical Ma hema ics 61 (2011) 1161–1185 1185
[10] F. Guillén-González, J.V. Gu ié ez-San ac eu, Condi ional s abili y and con e gence o a ully disc e e scheme o 3DNa ie –S okes equa ions wi h
mass di usion, SIAM J. Num. Anal. 46 (5) (2008) 2276–2308.
[11] F. Guillén-González, J.V. Gu ié ez-San ac eu, E o es ima es o a linea decoupled Eule –FEM scheme o a mass di usion model, Nume . Ma h. 117
(2011) 333–371.
[12] F. Hech , F eeFem++ manual, a ailable a h p://www. ee em.o g/ ++/.
[13] A. Kazhikho , Sh. Smagulo , The co ec ness o bounda y alue p oblems in a di usion model o an inhomogeneous fluid, So . Phys. Dokl. 22 (1)
(1977) 249–252.
[14] P.-L. Lions, Ma hema ical Topics in Fluid Dynamics, ol. 2, Incomp essible Models, Ox o d Uni e si y P ess, Uni ed Kingdom, 1996.
[15] P. Secchi, On he mo ion o iscous fluids in he p esence o di usion, SIAM J. Ma h. Anal. 19 (1988) 22–31.
[16] J. Simon, Compac se s in he space Lp(0,T;B), Ann. Ma . Pu a Appl. 146 (1987) 65–97.