Navier-Stokes-Brinkman方程的后驗(yàn)誤差估計(jì)

高 亮， 李 劍

(陜西科技大學(xué) 文理學(xué)院， 陜西 西安 710021)

0 引言

多孔介質(zhì)中流動(dòng)的模擬有許多應(yīng)用，包括儲(chǔ)層模擬，核廢料處理和二氧化碳封存.由于各種原因，這種模擬具有挑戰(zhàn)性.首先，區(qū)域的不規(guī)則性，使模型的幾何構(gòu)建復(fù)雜化.其次，地質(zhì)構(gòu)造由許多不同的物質(zhì)組成，具有不同的地質(zhì)性質(zhì).區(qū)域內(nèi)經(jīng)常出現(xiàn)裂縫和孔洞，改變了有效滲透率.建模這類(lèi)問(wèn)題的標(biāo)準(zhǔn)方法是耦合Darcy 和Navier-Stokes，并沿著界面執(zhí)行BJS條件[1-3].自由流區(qū)(裂縫、溶洞)采用Navier-Stokes 模型，而多孔區(qū)采用Darcy定律模型[4].然而，在儲(chǔ)層中，這兩種類(lèi)型的區(qū)域并沒(méi)有很好地分離，可能很難確定沿界面施加的適當(dāng)條件.

本文使用Navier-Stokes-Brinkman方程對(duì)多孔介質(zhì)中的流動(dòng)進(jìn)行建模，該方程將Navier-Stokes方程和Darcy方程合并為一個(gè)方程系統(tǒng)，Navier-Stokes-Brinkman方程根據(jù)系數(shù)的不同而簡(jiǎn)化為Navier-Stokes方程或Darcy方程，并提出用它來(lái)代替耦合的Navier-Stokes-Darcy方程.通過(guò)對(duì)系數(shù)的選擇，該方程可以同時(shí)對(duì)自由流動(dòng)和多孔區(qū)域進(jìn)行建模，從而解決沿界面的問(wèn)題.

本文提出了一種基于殘差的后驗(yàn)誤差估計(jì)方法，并使用該方法對(duì)離散化的Navier-Stokes-Brinkman問(wèn)題進(jìn)行自適應(yīng)網(wǎng)格細(xì)分，文獻(xiàn)中提出了后驗(yàn)誤差估計(jì)策略[5].與本文工作最接近的是Burda等[6-9]對(duì)Navier-Stokes方程的估計(jì),特別是Burda在二維空間中推導(dǎo)的Stokes-Brinkman問(wèn)題.

1 Navier-Stokes-Brinkman方程

(1)

這里μ為流體的粘度，u為流體速度，p為壓力，f為外力.式(1)分別是由動(dòng)量守恒定律和質(zhì)量守恒定律導(dǎo)出.在多孔區(qū)域Ωp中，流體流動(dòng)情況受Darcy定律控制：

(2)

這里K表示對(duì)稱(chēng)正定的滲透率張量.需要注意的是, 速度u和壓力p與Navier-Stokes方程中的相同項(xiàng)不同.在Navier-Stokes方程中，它們表示流體的實(shí)際速度和壓力，而在Darcy定律中，它們表示一些代表性元素體積的平均值.

在Navier-Stokes-Brinkman方程中，Navier-Stokes方程和Darcy方程組合成統(tǒng)一的方程組：

(3)

這里μ*表示流體的有效粘度.

納維爾-斯托克斯流和達(dá)西流都是Navier-Stokes-Brinkman方程選擇適當(dāng)?shù)摩?和K的極限情況：

(1)μ*=0時(shí)，Navier-Stokes-Brinkman方程會(huì)退化為Darcy方程；

(2)K-1?0時(shí)，Navier-Stokes-Brinkman方程會(huì)退化為Navier-Stokes方程.

由于μ*在多孔區(qū)域中對(duì)達(dá)西定律只有一點(diǎn)微擾的結(jié)果，所以可以在整個(gè)領(lǐng)域中選擇μ*=μ.

Navier-Stokes-Brinkman方程的邊界條件選取Dirichlet邊界條件，形式如下：

u=0,on?Ω

(4)

為了尋求方程式(3)和式(4)的解，定義如下空間：

定義X×M上的雙線性形式如下：

定義X×X×X上的三線性形式如下：

對(duì)于u,v,w∈X，滿(mǎn)足

b(u,v,w)=-b(u,w,v),b(u,v,v)=0,

|b(u,v,w)|≤N‖u‖0‖v‖0‖w‖0.

2 非線性微分方程的后驗(yàn)誤差估計(jì)

記L(U,V)為U到V的連續(xù)線性映射，其范數(shù)記為‖·‖L(U,v)，I(U,V)?L(U,V)為同胚開(kāi)子集，V的對(duì)偶空間V*∶=L(V,R)，設(shè)F∈C(U,V*)是一個(gè)連續(xù)可微函數(shù)，DF(u0)表示F在u0處的微分值.

定義Rmin形式如下 ：

2γ-1‖DF(u0)‖L(U,V*)}

對(duì)于任意的u∈[B(u0,Rmin)∩UD]，均有下面的結(jié)論[10]

‖F(xiàn)(u)‖S*≤‖DF(u0)‖L(U,V*)‖u-u0‖U

引理3存在常數(shù)C1,C2僅取決于定量hk，對(duì)于K∈Th,E∈Th,1≤q≤∞(q為向量維數(shù))以下誤差估計(jì)是有效的[11].

這里Phv是v的插值逼近函數(shù)，Th是一個(gè)有限元剖分，K是有限元剖分的子單元內(nèi)部，E是有限元剖分的子單元公共交界.

3 Navier-Stokes-Brinkman方程后誤差估計(jì)

本節(jié)構(gòu)造了Navier-Stokes-Brinkman方程的誤差估計(jì).令

定義式(3)和式(4)的變分形式如下：對(duì)于任意(v,q)∈X×M

〈F([u,p]),[v,q]〉:=μ((u,v)+

K-1(u,v)+((u·)u,v))-

(p,·v)+(q,·u)-(f,v)=

a(u,v)-d(v,p)+d(u,q)+

b(u,u,v)-(f,v).

(5)

使用混合有限元離散式(5)，并且構(gòu)造離散的速度空間Xh?X和離散的壓力空間Mh?M如下：

其中P1(K)表示一次多項(xiàng)式，P2(K)表示二次多項(xiàng)式.P2-P1有限元對(duì)滿(mǎn)足LBB條件

則問(wèn)題的離散化是穩(wěn)定的.

(6)

3.1 能量范數(shù)的后驗(yàn)誤差估計(jì)

本小節(jié)主要介紹能量范數(shù)中速度和壓力的后驗(yàn)誤差估計(jì).定義局部誤差指示子如下：

定理1記(u0,p0)是式(2)和式(3)的弱解，(uh,ph)是利用P2-P1元所求的逼近解，且滿(mǎn)足

〈F([uh,ph]),[vh,ph]〉=0,

則有如下估計(jì)：

證明：首先，證明F的導(dǎo)數(shù)的存在性，以及它在(u0,p0)∈X×M的鄰域中是連續(xù)的.

(u0·)uv+(u·)u0v)-p·v+

q·u]dx=a(u,v)-d(v,p)+d(u,q)+

b(u0,u,v)+b(u,u0,v),

由上述定義可知

〈F([u,p])-F([u0,p0])-DF([u-u0,p-

因此

故F是可微的.接下來(lái)，證明DF是Lipschitz連續(xù)的.

記DF1(·)為DF在[u1,p1]處的導(dǎo)數(shù).對(duì)于任意(v,q)∈Q，(u,p)∈B([u0,p0],R0)有下列估計(jì)成立：

〈DF1([u,p])-DF0([u,p]),[v,q]〉=

2N‖(u1-u0)‖0‖u‖0‖v‖0≤

2N‖[u1,p1]-[u0,p0]‖P‖[u,p]‖P‖v,q‖Q

由上式可得

2N‖[u,p]‖P≤2N(‖[u0,p0]‖P+R0)≤γ

取(vh,qh)∈Xh×Mh，

〈F([uh,ph])-Fh([uh,ph]),[vh,qh]〉=0

因此，

將上述不等式結(jié)合得到所證結(jié)論.

3.2 L2范數(shù)的后驗(yàn)誤差估計(jì)

對(duì)于速度的L2范數(shù)，考慮如下線性化問(wèn)題.尋求(φ,ξ)∈X×M,(u0,p0)是原問(wèn)題的弱解，對(duì)于某些給定的g∈Y和任意(v,q)∈X×M，使得

a(φ,v)-d(φ,q)+b(u0,v,φ)+

b(v,u0,φ)+d(v,ξ)=(g,v)

(6)

可知雙線性形式a(φ,v)-d(φ,q)+d(v,ξ)是連續(xù)的.由于Necas定理，上式有唯一的對(duì)偶解(φ,ξ)∈X×M，且滿(mǎn)足

‖φ‖2,Ω+‖ξ‖1,Ω≤C‖g‖0,Ω,

這里，

δK=hK‖·uh‖0,K,

定理2令(uh,ph)是(u0,p0)的近似，在定理1的條件下，可得結(jié)論如下：

證明：取v=g=u0-uh,q=p0-ph代入式(6)，可得

d(u0-uh,ξ)+b(u0,u0-uh,φ)+b(u0-

uh,u0,φ),

另一方面，由式(5)及其離散形式作差可得

a(u0-uh,vh)-d(vh,p0-ph)+d(u0-uh,qh)+

b(uh,u0-uh,vh)+b(u0-uh,uh,vh)=0,

將上述兩式求和，并令(vh,qh)=(Phφ,Phξ),

Phφ,p0-ph)+d(u0-uh,ξ-Phξ)+b(u0,u0-

uh,φ-Phφ)+b(u0-uh,uh,φ-Phφ)+b(u0-

uh,u0-uh,φ)=(f+μ(△uh-(uh·)uh-

K-1uh)-ph,φ-Phφ)-(·uh,ξ-Phξ)+

uh)·)(u0-uh),φ)≤

由上式可得速度的L2范數(shù)的誤差估計(jì).

4 數(shù)值算例

本節(jié)中，使用兩個(gè)數(shù)值算例去證明本文理論結(jié)果的正確性，并且驗(yàn)證了通過(guò)后驗(yàn)誤差估計(jì)法所得到的誤差指示子對(duì)于Navier-Stokes-Brinkman方程自適應(yīng)過(guò)程的有效性.自適應(yīng)網(wǎng)格算法如下：

給定一個(gè)網(wǎng)格剖分Th，并計(jì)算出每一個(gè)子單元K所對(duì)應(yīng)的誤差指示子ηK，再給定一個(gè)常數(shù)θ∈(0,1).

(2)如果ηK≥θηT,max，則對(duì)單元K進(jìn)行網(wǎng)格加密.加密后得到一個(gè)新的網(wǎng)格剖分，返回步驟(1).

注：在自適應(yīng)過(guò)程中，有些單元需要進(jìn)行粗化處理，其基本思想是將誤差太小的單元聚集在一起.

4.1 圓形區(qū)域的數(shù)值算例

第一個(gè)算例中流體所在的區(qū)域?yàn)榘霃綖?的圓盤(pán)，圓心與邊界之間有一條裂紋.從初始的網(wǎng)格剖分圖1(a)開(kāi)始自適應(yīng)計(jì)算，經(jīng)過(guò)兩次迭代計(jì)算得到了如圖1(b)所示的網(wǎng)格剖分，然后再經(jīng)一次迭代計(jì)算，結(jié)果如圖1(c)所示.圖2所體現(xiàn)的是由于在裂紋處的自適應(yīng)迭代計(jì)算，壓力線的剖面變得更加光滑，這也意味著計(jì)算結(jié)果更加精確.

(a)初始剖分 (b)自適應(yīng)迭代兩次 (c)自適應(yīng)迭代三次圖1 P2-P1有限元對(duì)在圓形區(qū)域自適應(yīng)過(guò)程

(a)初始剖分 (b)自適應(yīng)迭代兩次 (c)自適應(yīng)迭代三次圖2 P2-P1有限元對(duì)在圓形區(qū)域自 適應(yīng)過(guò)程中壓力線變化

4.2 L形區(qū)域的數(shù)值算例

第二個(gè)數(shù)值算例中選取了L型區(qū)域Ω=(0,1)×(0,0.5)/(0,0.5)×(0.5,1)，并且邊界為Dirichlet邊界條件，進(jìn)行自適應(yīng)迭代計(jì)算，得到如下結(jié)果，從初始的網(wǎng)格剖分圖3(a)開(kāi)始自適應(yīng)計(jì)算，經(jīng)過(guò)兩次迭代計(jì)算得到了如圖3(b)所示的網(wǎng)格剖分，然后再經(jīng)一次迭代計(jì)算，結(jié)果如圖3(c)所示.通過(guò)在拐角處進(jìn)行自適應(yīng)迭代計(jì)算，網(wǎng)格剖分在拐角處進(jìn)行了加密，計(jì)算結(jié)果更加精確；同時(shí)也避免了因統(tǒng)一加密而導(dǎo)致的計(jì)算量的增加.圖4所體現(xiàn)的是由于在拐角處的自適應(yīng)迭代計(jì)算，壓力線的剖面變得更加光滑，這也意味著在計(jì)算量小的前提下，得到了更加準(zhǔn)確地表示.

(a)初始剖分 (b)自適應(yīng)迭代兩次 (c)自適應(yīng)迭代三次圖3 P2-P1有限元對(duì)在L型區(qū)域自適應(yīng)過(guò)程

(a)初始剖分 (b)自適應(yīng)迭代兩次 (c)自適應(yīng)迭代三次圖4 P2-P1有限元對(duì)在L型區(qū)域自 適應(yīng)過(guò)程中壓力線變化

5 結(jié)論

本文在混合有限元法的基礎(chǔ)上，推導(dǎo)出了Navier-Stokes-Brinkman問(wèn)題的殘差型后驗(yàn)誤差估計(jì)，得到了可計(jì)算的能量范數(shù)全局上界以及速度的L2范數(shù)的誤差估計(jì)值，并且利用計(jì)算所得的估計(jì)值驅(qū)動(dòng)自適應(yīng)算法去求解Navier-Stokes-Brinkman問(wèn)題，通過(guò)數(shù)值算例，證實(shí)了理論的正確性以及算法的高效性.