ˇ sen´ı Navierov´ych-Stokesov´ych rovnic metodou Reˇ koneˇcn´ych prvk˚ u ˇ Libor Cerm´ ak prosinec 2009 Abstrakt Text obsahuje klasickou a variaˇcn´ı formulaci 2D-´ ulohy nestlaˇciteln´eho nestacion´arn´ıho proudˇen´ı, popis prostorov´e diskretizace metodou koneˇcn´ ych prvk˚ u, ˇcasovou diskretizaci implicitn´ı BDF1 a BDF2 formul´ı, linearizaci Oseenovou a Newtonovou metodou, stabilizaci metodou SUPG/PSPG/LSIC nebo GLS.
Oznaˇcen´ı: Ω [0, T ] x = (x1 , x2 ) n = (n1 , n2 ) t u = (u1 , u2 ), ui = ui (x, t) p(x, t) ε(u) = {εij (u)}2i,j=1 f = (f1 , f2 ), fi = fi (x, t) ν
... ... ... ... ... ... ... ... ... ...
rovinn´a oblast ˇcasov´ y interval bod rovinn´e oblasti Ω jednotkov´a vnˇejˇs´ı norm´ala hranice ∂Ω ˇcas rychlost kinematick´ y tlak tenzor rychlosti deformace kinematick´a setrvaˇcn´a s´ıla kinematick´a viskozita
Uˇz´ıv´ame Einsteinovu sumaˇcn´ı konvenci: pokud se v jednom ˇclenu vyskytuje index dvakr´at, sˇc´ıt´a se pˇres nˇej od 1 do 2. Navier-Stokesovy rovnice: ∂ui ∂ui ∂εij (u) ∂p + uj − 2ν + = fi , ∂t ∂xj ∂xj ∂xi
i = 1, 2,
(x, t) ∈ Ω × (0, T ).
(1)
Rovnice kontinuity: ∂ui = 0, ∂xi
(x, t) ∈ Ω × (0, T ).
(2)
Okrajov´a podm´ınka Dirichletova typu, pˇredeps´ana rychlost: ui = gi ,
i = 1, 2,
(x, t) ∈ Γ1 × (0, T ). 1
(3)
Okrajov´a podm´ınka Neumannova typu, pˇredeps´ano napˇet´ı: 2νεij (u)nj − pni = σi ,
(x, t) ∈ Γ2 × (0, T ).
i = 1, 2,
(4)
Poˇc´ateˇcn´ı podm´ınka: ui (x, 0) = u0i (x),
x ∈ Ω.
(5)
Sloˇzky tenzoru rychlosti deformace: 1 ∂ui ∂uj εij (u) = + , i, j = 1, 2. 2 ∂xj ∂xi
(6)
Hranice ∂Ω je sjednocen´ım ˇca´st´ı Γ1 a Γ2 : ¯1 ∪ Γ ¯ 2, ∂Ω = Γ
Γ1 ∩ Γ2 = ∅.
(7)
Zvol´ıme testovac´ı funkci v = (v1 , v2 ), vi = vi (x), s vlastnost´ı: vi = 0,
x ∈ Γ1 ,
i = 1, 2,
(8)
a testovac´ı funkci q. Rovnici (1) n´asob´ıme vi , integrujeme pˇres Ω, uˇzijeme Greenovu vˇetu, uplatn´ıme podm´ınku (8) a okrajovou podm´ınku (4). Rovnici (2) n´asob´ıme q a integrujeme pˇres Ω. Takto z´ıskan´e rovnice seˇcteme a dostaneme a(u, p, u; v, q) =0, kde Z ∂ui ∂vi ∂ui ∂ui vi + wj vi + 2νεij (u)εij (v) − p − fi vi + q dx− a(u, p, w; v, q) = ∂t ∂xj ∂xi ∂xi Ω Z − σi vi dS. (9) Γ2
¯ h = S e, kaˇzd´e dva r˚ Oblast Ω triangulujeme, prvky znaˇc´ıme e, Ω uzn´e prvky jsou bud’to disjunktn´ı nebo maj´ı spoleˇcnou jednu stranu nebo jeden vrchol. Necht’ ¯ 1h ∪ Γ ¯ 2h , ∂Ωh = Γ
Γ1h ∩ Γ2h = ∅,
¯1 ∩ Γ ¯2 = Γ ¯ 1h ∩ Γ ¯ 2h . Γ
Aproximace uh (x, t) rychlosti u je po prvc´ıch polynom stupnˇe dv jednoznaˇcnˇe urˇcen´ y sv´ ymi hodnotami v uzlech Pje prvku e, ¯ h ), wi ∈ Pdv (e), wi (P e ) = gi (P e , t) ∀P e ∈ Γ ¯ 1 } (10) uh (·, t) ∈ Vhg (t) = {w | wi ∈ C(Ω j j j e ¯ 1h ⇒ Pje ∈ Γ ¯ 1 ). Aproximace vh (x) testovac´ı funkce v: (pˇredpokl´ad´ame, ˇze Pje ∈ Γ ¯ 1 }. ¯ h ), wi ∈ Pdv (e), wi (P e ) = 0 ∀P e ∈ Γ vh ∈ Vh0 = {w | wi ∈ C(Ω j j e 2
(11)
Aproximace ph (x, t) tlaku p a aproximace qh (x) odpov´ıdaj´ıc´ı testovac´ı funkce q je po prvc´ıch polynom stupnˇe dp (obvykle dp ≤ dv ): ¯ h ), w ∈ Pdp (e)}. ph (·, t), qh ∈ Mh = {w ∈ C(Ω (12) e Definujme stabilizaˇcn´ı ˇclen: X Z ∂ui ∂ui ∂εij (u) ∂p as (u, p, w; v, q) = + wj − 2ν + − fi ψie (w; v, q) dx+ ∂t ∂x ∂x ∂x j j i e e (13) X Z ∂ui ∂vj e + δ dx. ∂xi ∂xj e e V nˇem jako testovac´ı funkci bereme: ∂vi ∂εi` (v) ∂q e e e ψi (w; v, q) = τu w` − τs 2ν + τpe . ∂x` ∂x` ∂xi
(14)
Vhodnou volbu stabilizuj´ıc´ıch ˇc´ıseln´ ych parametr˚ u τue , τpe , τse a δ e uvedeme pozdˇeji. Vˇsimnˇete si: je-li u a p klasick´e ˇreˇsen´ı rovnic (1), (2), pak as (u, p, u; v, q) = 0. Definujme formu A jako souˇcet forem a, as : A(u, p, w; v, q) = a(u, p, w; v, q) + as (u, p, w; v, q).
(15)
ˇ ad´ame, aby pro pˇribliˇzn´e ˇreˇsen´ı uh , ph platilo: Z´ (uh , ph )(·, t) ∈ Vhg (t) × Mh : A(uh , ph , uh ; vh , qh ) = 0 ∀(vh , qh ) ∈ Vh0 × Mh , t ∈ (0, T ). Na intervalu [0, T ] zvol´ıme dˇelen´ı 0 = t0 < t1 < · · · < tn−1 < tn < · · · < tm = T. ˇ Casov´ y krok ∆tn = tn − tn−1 ,
n = 1, 2, . . . , m.
V rovnici (16) zapsan´e pro ˇcas t = tn nahrad´ıme ˇcasovou derivaci zpˇetnou diferenc´ı: ∂uh (x, tn ) uh (x, tn ) − uh (x, tn−1 ) ≈ . ∂t ∆tn Pˇribliˇzn´e ˇreˇsen´ı v ˇcase tn znaˇc´ıme unh ≈ u(·, tn ),
pnh ≈ p(·, tn ). 3
(16)
Z formy a dostaneme dvˇe formy b a c. Biline´arn´ı forma b obsahuje diskretizaci ˇcasov´e derivace, oproti formˇe a v n´ı vˇsak chyb´ı konvekˇcn´ı ˇclen: Z ui − un−1 ∂ui ∂vi i n vi + 2νεij (u)εij (v) − p − fi vi + q dx− (17) b(u, p; v, q) = ∆tn ∂xi ∂xi Ωh Z − σin vi dS. Γ2h
Konvekˇcn´ı ˇclen tvoˇr´ı triline´arn´ı formu Z ∂ui c(u, w; v) = wj vi dx. ∂xj Ωh
(18)
Tak´e ve stabilizaˇcn´ım ˇclenu as nahrad´ıme ˇcasovou derivaci diferenc´ı a dostaneme formu bs (u, p, w; v, q) = X Z ui − un−1 ∂ui ∂εij (u) ∂p i n + wj − 2ν + − fi ψie (w; v, q) dx+ = ∆t ∂x ∂x ∂x n j j i e e Z ∂ui ∂vj δe + dx. ∂xi ∂xj Ωh
(19)
Pˇribliˇzn´e ˇreˇsen´ı tedy m´a splˇ novat: (unh , pnh ) ∈ Vhg (tn ) × Mh : B(unh , unh , pnh ; v, q) = 0 ∀(vh , qh ) ∈ Vh0 × Mh ,
n = 1, 2, . . . ,
kde
(20)
B(uh , wh , ph ; v, q) = b(uh , ph ; v, q) + c(uh , wh ; v) + bs (uh , ph , wh ; v, q). Pro n > 1 lze m´ısto implicitn´ı Eulerovy formule (nebo-li BDF1 formule ˇr´adu 1) pouˇz´ıt pˇresnˇejˇs´ı BDF2 formuli (ˇra´du 2) zaloˇzenou na aproximaci 3uh (x, tn ) − 4uh (x, tn−1 ) + uh (x, tn−2 ) ∂uh (x, tn ) ≈ . ∂t 2∆tn Ve form´ach b, bs v tom pˇr´ıpadˇe nahrad´ıme ˇclen
ui − un−1 i ∆tn
3ui − 4uin−1 + un−2 i . 2∆tn
ˇclenem
´ Uloha (20) je neline´arn´ı a to jak v konvekˇcn´ım ˇclenu c tak ve stabilizaˇcn´ım ˇclenu bs . Je proto tˇreba pouˇz´ıt iteraˇcn´ı metodu. Aproximace poˇc´ıtan´e iteraˇcn´ım procesem znaˇc´ıme n,0 n,k c´ateˇcn´ı aproximaci un,0 lze vz´ıt hodnoty horn´ım indexem k, tj. un,k h , ph h , ph . Jako poˇ n−1 , p z pˇ r edchoz´ ıho ˇ c asu t . un−1 n−1 h h Jednoduch´a linearizace, zn´am´a jako Oseenova metoda, poˇc´ıt´a un,k , pn,k z rovnice n,k−1 n,k B(un,k , ph ; v, q) = 0, h , uh
k = 1, 2, . . . . 4
(21)
D˚ umyslnˇejˇs´ı linearizace zaloˇzen´a na Newtonovˇe metodˇe aplikovan´e na konvekˇcn´ı ˇclen c vede na sch´ema n,k−1 n,k n,k−1 B(un,k , ph ; v, q) + c(uhn,k−1 , un,k , un,k−1 ; v) = 0, h , uh h ; v) − c(uh h
k = 1, 2, . . . . (22)
Parametry stabilizace τue , τpe , τse a δ e se nastavuj´ı uˇzit´ım jedn´e nebo souˇcasnˇe nˇekolika z n´asleduj´ıc´ıch metod, viz [6]: 1) SUPG (streamline upwind Petrov-Galerkin), 2) PSPG (pressure stabilizing Petrov-Galerkin), 3) LSIC (least squares on incompressibility constraint). Dalˇs´ı alternativou je pouˇzit´ı metody GLS (Galerkin least squares). N´asleduj´ı tˇri osvˇedˇcen´e volby stabilizaˇcn´ıch parametr˚ u: 1) SUPG+PSPG+LSIC pro pˇr´ıpad, ˇze rychlosti jsou aproximov´any polynomem stupnˇe dv vˇetˇs´ım neˇz je stupeˇ n dp aproximace tlaku. Popul´arn´ı Taylor-Hookovy prvky pouˇz´ıvaj´ı dv = dp + 1. Vol´ıme, viz [3]: τue = τpe =
1 ¯e 2 [h ] , d2v
τse = 0,
δe = 1 ,
(23)
¯ e je charakteristick´ ¯ e = he kde h y pr˚ umˇer troj´ uheln´ıka e, napˇr´ıklad h s´ı max je nejdelˇ strana e. 2) SUPG+PSPG+LSIC pro pˇr´ıpad, ˇze rychlosti i tlak jsou aproximov´any polynomem t´ehoˇz stupnˇe. Vol´ıme, viz [1]: ¯ e |¯ ¯ e ]2 ¯ e ]2 h ue,n,k−1 | [h [h e δ =ν 1+ + , τue = τpe = e , τse = 0, (24) ν ν∆t δ kde |¯ ue,n,k−1 | je charakteristick´a d´elka vektoru rychlosti ue,n,k−1 na troj´ uheln´ıku e, e,n,k−1 e napˇr´ıklad rychlosti uC v tˇeˇziˇsti C troj´ uheln´ıka e. 3) V pˇr´ıpadˇe metody GLS vol´ıme, viz [2]: 1 0 ≤ Ree < 1 , 4νλe , max τue = τpe = τse = 1 , Ree ≥ 1 , p e e,n,k−1 λmax |¯ u |
kde Ree =
δ e = τse |¯ ue,n,k−1 |2 .
|¯ ue,n,k−1 | p , 4ν λemax (25)
Pˇritom λemax je nejvˇetˇs´ı vlastn´ı ˇc´ıslo u ´lohy Z Z ∂εij (v) ∂εik (v) e dx = λe εij (v)εij (v) dx ∀v ∈ Vhe /Vh0 , ∂x ∂x j k e e (26) Z e Vh0 = v ∈ Vhe εij (v)εij (v) dx = 0 . Vhe = {v | vi ∈ Pdv (e), i = 1, 2}, e
5
e e . rozum´ıme faktorov´ y prostor prostoru Vhe podle podprostoru Vh0 Symbolem Vhe /Vh0 e e e ypoˇctu lze faktorov´ y prostor Vh /Vh0 Prostor Vh0 m´a dimenzi 3. Pˇri praktick´em v´ nahradit vhodn´ ym podprostorem Vˆhe prostoru Vhe dimenze 2dv − 3, napˇr´ıklad
Vˆhe = {v ∈ Vhe | v1 (P1e ) = v1 (P2e ) = v2 (P1e ) = 0}. Vlastn´ı ˇc´ısla λemax nez´avisej´ı na n ani na k, poˇc´ıtaj´ı se tedy jen jednou. ¯ e troj´ Jako charakteristick´ y pr˚ umˇer h uheln´ıka e lze v (24), (25) pouˇz´ıt rovnˇeˇz pr˚ umˇer heu ¯ e,n,k−1 , viz [5]: troj´ uheln´ıka e ve smˇeru rychlosti u 2|¯ ue,n,k−1 | , e,n,k−1 · ∇Le | |¯ u j j=1
heu = P3
(27)
kde {Lej }3j=1 jsou line´arn´ı b´azov´e funkce pˇr´ısluˇsn´e vrchol˚ um troj´ uheln´ıka e, viz [6].
Literatura [1] M. Braack, E. Burman, V. John, G. Lube: Stabilized finite element for the generalized Oseen problem, Comput. Methods Appl. Mech. Engrg. 196 (2007), 853-866. [2] L. P. Franca, A. L. Madureira: Element diameter free stability parameters for stabilized methods applied to fluids, Comput. Methods Appl. Mech. Engrg. 105 (1993), 395-403. [3] T. Gelhard, G. Lube, M.A. Olshanskii, J.H. Starcke: Stabilized finite element schemes with LBB-stable elements for incompressible flows, J. Comp. Appl. Math. 177 (2005), 243-267. [4] P.G. Gresho, R.S. Sani: Incompressible Flow and the Finite Element Method. Vol 2: Isothermal Laminar Flow, John Wiley & Sons, Chichester, 2000. [5] P. Sv´aˇcek: Numerical simulation of aeroelastic problems with consideration of nonlinear effects, Engineering MECHANICS, Vol. 16, No. 1, (2009), 13-28. [6] T.E. Tezduyar: Finite elements in fluids: Stabilized formulations and moving boundaries and interfaces, Computer & Fluids 36 (2007), 191-206.
6
