求解Navier - Stokes/Darcy方程的時間異步穩(wěn)定化特征線方法

，

(太原理工大學數(shù)學學院， 山西太原030024)

流- 固耦合問題一直是流體力學研究的熱點，特別是自由流與多孔介質流的耦合問題，在實際的工業(yè)與工程生產中和醫(yī)學、地球科學中都有很大的應用，例如，模擬平緩河床的流動過程、模擬工業(yè)中石油開采的過程、描述血液在經過血管壁的流體運動過程等。求解這類耦合問題的解耦方法[1-7]已經有很多，例如，文獻[4-5]中分別提出區(qū)域分解法、兩重網格法來求解耦合問題，文獻[6]中提出時間異步的方法來求解非穩(wěn)態(tài)的Stokes- Darcy耦合問題，文獻[7]中將時間異步拓展到非穩(wěn)態(tài)的Navier-Stokes/Darcy方程。

基于時間異步的優(yōu)勢以及特征線方法[8-9]在處理對流占優(yōu)問題上的優(yōu)點，本文中提出求解耦合問題的時間異步穩(wěn)定化特征線有限元法， 進行數(shù)值格式的穩(wěn)定性和收斂性分析，給出數(shù)值算例并得出相應的結論。

1 模型與變分形式

在自由流區(qū)域Ωf中的流體流動可以由如下方程控制，

?u?t-υ?2u+(u·?)u+?p=f1, Ωf×(0,T],?·u=0, Ωf×(0,T],u(x,0)=u0, Ωf ,u=0, Ω—f?！?0,T],ì?í????????(1)

式中：x為Ωf內的位置；t為時間；u(x,t)為區(qū)域Ωf中流體的流速；u0為u在t=0時刻的取值；2為拉普拉斯算子；為梯度算子；·u為u的散度； (u·)為u與的內積；p(x,t)為區(qū)域Ωf中流體的壓力；f1為外力；υ>0為流體的黏性系數(shù)；常數(shù)T>0為最終時刻。

對于多孔介質流體區(qū)域Ωp的流動過程由如下方程控制，

(2)

交界面條件包括質量守恒條件、力的平衡條件和Beavers - Joseph - Saffman(BJS)條件，即

式中：nf、np為相應的Ωf、Ωp的單位外法向量；τi,i=1,2,…,d-1為交界面Γ的單位正切向量；g為重力加速度；α為取決于多孔介質的性質并且通過實驗確定的正參數(shù)。

在給出耦合問題的弱形式之前，首先引入一些重要的記號。考慮函數(shù)空間,

Hf、Hp中元素的范數(shù)分別定義為

其中(· ,· )Ω為相應區(qū)域Ω內的內積。Navier -Stokes/Darcy問題(1)—(3)弱形式如下：求w=(u,φ)∈W以及p∈Q，使得對?t∈(0,T]都有?z=(v,ψ)∈W,?q∈Q，

其中

[wt,z]=(ut,v)Ωf+gS0(φt,ψ)Ωp,

a(w,z)=af(u,v)+ap(φ,ψ)+a1(w,z),

af(u,v)=v(u,v)Ωf+

ap(φ,ψ)=g(Kφ,ψ)Ωpb(z,p)=-(p,·v)Ωf,

B(u,v,w)=((u·)v,w)Ωf,

(f,z)=(f1,v)Ωf+g(f2,ψ)Ωp。

2 基于時間異步的穩(wěn)定化特征線法

設τh為區(qū)域Ωf∪Ωp依賴于正參數(shù)h>0的擬一致三角剖分，當d=2時，由三角形構成，當d=3時，由四面體構成。定義有限元子空間Wh=Hfh×Hph?W,Qh?Q[7]，S為三角單元，

為了克服最低階有限元對不滿足LBB(Ladyzhenskaya - Brezzi - Babuska)條件，引入標準L2- 投影算子Πh∶L2(Ω)→P0,P0={p∈Q∶pS∈P0(S),?S∈τh},穩(wěn)定化項表達式定義為Gh(· ,·)，

G(p,q)=(p-Πhp,q-Πhq),

且投影算子Πh具有下列性質[10]：

(p,qh)=(Πhp,qh), ?p∈Q,qh∈P0,

則帶有穩(wěn)定化項的表示式為求wh=(uh,φh)∈Wh以及ph∈Qh，使得對于?t∈(0,T]，?zh=(vh,ψh)∈Wh,qh∈Qh,都有

(wht,zh)+a(wh,zh) +b(zh,ph) +B(uh;uh,vh) +

βG(ph,qh) =(f,zh),

b(wh,qh) =0,

wh(0) =w0,

式中β為正的穩(wěn)定化參數(shù)。

其中

定義投影算子，滿足

?t∈[0,T],?zh∈Wh,qh∈Qh,滿足

a(Pwhw(t),zh)+b(zh,Pphp(t))+βGh(Pphp(t),q)=a(w(t),zh)+b(zh,p(t)),b(Pwhw(t),qh)=0。ì?í??????(4)

引理1[1]假設(w(t),p(t))有足夠的光滑性，則有下列估計：

假定時間層是均勻的，對于多孔介質流區(qū)域Ωp的每一層時間Sk都對應自由流區(qū)域的相對應的時間層tmk，即滿足如下對應關系：

Sk=kΔs,k=0,1,…,M, Δs=rΔt，

tm=mΔt,m=0,1,…,N,

其中

gS0?mk+1h-?mkhΔs,ψh()+ap(?mk+1h, ψh)=g(fm+12,ψh)+g∫ΓψhSmk·nf dx,?m0h=Pph?0 。ì?í??????(6)

令k=k+ 1進行循環(huán)，直到k=M- 1結束。

如何選取適當?shù)膔非常重要。r的選取是由流體的參數(shù)、邊界條件和初始條件決定的，這就導致r的選取比較復雜，因此對于r值的選取是通過數(shù)值實驗來確定的。

3 穩(wěn)定性

為了更方便地分析穩(wěn)定性,首先引入引理2、3。

(7)

(8)

此外,用逆不等式，?wh,zh∈Wh可以推導：

(9)

在穩(wěn)定性證明的這部分中，都假設滿足

(10)

(11)

(12)

(13)

(14)

合并式(13)、(14)，有

(15)

運用Young不等式、Holder’s不等式、Poincare不等式，式(15)右端的第1、2項可以放縮為

(16)

利用引理3中的式(8), 以及式(15)右端的第3項可以放縮為

(17)

利用引理2，式(15)的右端最后一項放縮為

(18)

合并式(15)—(18)，對k=0,1,…,l,0≤l≤M-1求和，得到

(19)

(20)

通過與式(15)相似的討論，有

(21)

考慮上式的特殊情形l=-1, 得到

合并式(19)、(21), 整理得

(22)

利用離散的Gronwall’s引理，可以得到時間步長[s1,T]內的穩(wěn)定性，即

4 收斂性

首先假設方程(1)—(3)的解(u,p,φ)滿足下列正則性條件[7]:

1)u∈L∞(0,T,H2(Ωf)2),φ∈L∞(0,T,H2(Ωp));

2)ut∈L2(0,T,H2(Ωf)2),φl∈L2(0,T,H2(Ωp));

3)p∈L∞(0,T,H1(Ωf)),pt∈(0,T,H1(Ωf))。

根據(jù)引理1的近似性質，可以得到

(23)

為了估計數(shù)值解與精確解的誤差， 對em進行分析。

結合式(4)—(6)，依次相減并整理，對于((vh,ψh),qh)∈(Wh,Qh)，

em+1-emΔt,vh()+af(em+1,vh)+b(vh,ηm+1)+βG(ηm+1,qh)=u^mh-umhΔt+u(tm+1)-u～m+1Δt-u—(xm)-u～mΔt,vh()-u(tm+1)-u—(xm)Δt-ut(tm+1)-u(tm+1)·?u(tm+1),vh()+g∫Γ(?～m+1-?～m)vh·nfdx+g∫Γ(?～m-?mh)vh·nfdx,b(em+1, qh)=0。ì?í??????????????????(24)

(25)

定理2 如果精確解(u,p,φ)是光滑的， 并且初始解的近似足夠精確(即e0=0,0=0) ， 以及當時間步長Δt滿足限制條件對于l=0, 1, …,M-1, 有誤差估計

(26)

證明：將vh=2Δtem+1和qh==2Δtηm+1代入式(24)， 使用無散度條件, 對m=mk,mk+1, …,mk+1-1, 求和得到

(27)

將ψh=2Δsmk+1代入式(25)，利用和得到

(28)

合并式(27)、(28)，有

(29)

利用文獻[12]中引理3.4中的式(3.15)，以及Holder’s不等式,B1可以放縮為

(30)

利用引理4以及Holder’s不等式,B2可以放縮為

(31)

(32)

(33)

B5可以放縮為

(34)

(35)

對式(35)中的第2—5項估計，得到

(36)

將式(36)代入式(35)，利用離散的Gronwall’s引理，可以推導出定理2的結果。證畢。

注2 定理2分析uh在L∞(0,T;L2(Ωf))下大時間步長上的誤差估計。由定理2可知，uh在L2(0,T;Hf)下以及φh在L2(0,T;Hp)下的誤差估計是最優(yōu)的。

下面證明Navier - Stokes方程在較小時間步長上的誤差估計。

定理3 如果滿足定理2的條件，則對于J=0,1,…,r-1和k=0,1,…,M-1，有

(37)

證明：將vh= 2Δtem+1,qh= 2Δtηm+1代入式(24)，利用及無散度性質對m=mk,mk+1,…,mk+J求和，與定理2的證明過程相似，推導出

(38)

根據(jù)式(36)，以及定理2,可以得出

將以上估計式代入式(38)，應用離散的Gronwall’s引理，定理3即得證。

注3 定理3分析了在較小時間步長的條件下，uh在L∞(0,T;L2(Ωf))上的誤差估計。

5 數(shù)值實驗

以下給出數(shù)值算例來驗證本文中所提出的新方法的有效性。

在區(qū)域Ωf=[0,1]×[1,2]和Ωp=[0,1]×[0,1]交界面為Γ=(0,1)×{1}上來研究該耦合問題，給定問題的精確解為

u=(u1,u2)=([x2(y-1)2+y]cost,

[2-πsin(πx)]cost)，

p=[2-πsin(πx)]sin(0.5πy)cost,

φ=[2-πsin(πx)][1-y-cos(πy)]cost。

表1 當時間步長、最終時刻固定為Δt=0.01 s ,T=1 s時網格尺寸不斷細化的數(shù)值結果

表2 當網格大小、最終時刻固定為時時間步長不斷細化的數(shù)值結果

6 結論

在時間異步的基礎上，將穩(wěn)定化方法與特征線法相結合，對求解非穩(wěn)態(tài)的Navier - Stokes/Darcy方程的基于時間異步的穩(wěn)定化特征線有限元方法進行了研究。數(shù)值實驗證明了方法的有效性。新的方法具有更好的穩(wěn)定性以及最優(yōu)誤差估計。

對于求解非穩(wěn)態(tài)的耦合模型而言，最難處理是Navier - Stokes方程中非線性項。在未來的研究中，可以通過基于牛頓迭代法的兩重網格方法進行處理，進一步改善解的精度與求解的效率。