Stability and convergence for a complete model of mass diffusion

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.