，

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

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

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

1 模型與變分形式

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

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

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

(2)

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

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

在給出耦合問題的弱形式之前，首先引入一些重要的記號(hào)?？紤]函數(shù)空間,

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

其中(· ,· )Ω為相應(yīng)區(qū)域Ω內(nèi)的內(nèi)積。Navier -Stokes/Darcy問題(1)—(3)弱形式如下：求w=(u,φ)∈W以及p∈Q，使得對(duì)?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 基于時(shí)間異步的穩(wěn)定化特征線法

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

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

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

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

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

則帶有穩(wěn)定化項(xiàng)的表示式為求wh=(uh,φh)∈Wh以及ph∈Qh，使得對(duì)于?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]假設(shè)(w(t),p(t))有足夠的光滑性，則有下列估計(jì)：

假定時(shí)間層是均勻的，對(duì)于多孔介質(zhì)流區(qū)域Ωp的每一層時(shí)間Sk都對(duì)應(yīng)自由流區(qū)域的相對(duì)應(yīng)的時(shí)間層tmk，即滿足如下對(duì)應(yīng)關(guān)系：

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進(jìn)行循環(huán)，直到k=M- 1結(jié)束。

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

3 穩(wěn)定性

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

(7)

(8)

此外,用逆不等式，?wh,zh∈Wh可以推導(dǎo)：

(9)

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

(10)

(11)

(12)

(13)

(14)

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

(15)

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

(16)

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

(17)

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

(18)

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

(19)

(20)

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

(21)

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

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

(22)

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

4 收斂性

首先假設(shè)方程(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的近似性質(zhì)，可以得到

(23)

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

結(jié)合式(4)—(6)，依次相減并整理，對(duì)于((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) ， 以及當(dāng)時(shí)間步長(zhǎng)Δt滿足限制條件對(duì)于l=0, 1, …,M-1, 有誤差估計(jì)

(26)

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

(27)

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

(28)

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

(29)

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

(30)

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

(31)

(32)

(33)

B5可以放縮為

(34)

(35)

對(duì)式(35)中的第2—5項(xiàng)估計(jì)，得到

(36)

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

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

下面證明Navier - Stokes方程在較小時(shí)間步長(zhǎng)上的誤差估計(jì)。

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

(37)

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

(38)

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

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

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

5 數(shù)值實(shí)驗(yàn)

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

在區(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 當(dāng)時(shí)間步長(zhǎng)、最終時(shí)刻固定為Δt=0.01 s ,T=1 s時(shí)網(wǎng)格尺寸不斷細(xì)化的數(shù)值結(jié)果

表2 當(dāng)網(wǎng)格大小、最終時(shí)刻固定為時(shí)時(shí)間步長(zhǎng)不斷細(xì)化的數(shù)值結(jié)果

6 結(jié)論

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

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