求解Navier - Stokes/Darcy方程的时间异步稳定化特征线方法

2019-03-23，

济南大学学报（自然科学版） 2019年2期

，

(太原理工大学数学学院，山西太原030024)

流- 固耦合问题一直是流体力学研究的热点，特别是自由流与多孔介质流的耦合问题，在实际的工业与工程生产中和医学、地球科学中都有很大的应用，例如，模拟平缓河床的流动过程、模拟工业中石油开采的过程、描述血液在经过血管壁的流体运动过程等。求解这类耦合问题的解耦方法[1-7]已经有很多，例如，文献[4-5]中分别提出区域分解法、两重网格法来求解耦合问题，文献[6]中提出时间异步的方法来求解非稳态的Stokes- Darcy耦合问题，文献[7]中将时间异步拓展到非稳态的Navier-Stokes/Darcy方程。

基于时间异步的优势以及特征线方法[8-9]在处理对流占优问题上的优点，本文中提出求解耦合问题的时间异步稳定化特征线有限元法，进行数值格式的稳定性和收敛性分析，给出数值算例并得出相应的结论。

1 模型与变分形式

在自由流区域Ω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)为区域Ωf中流体的流速；u0为u在t=0时刻的取值；2为拉普拉斯算子；为梯度算子；·u为u的散度； (u·)为u与的内积；p(x,t)为区域Ωf中流体的压力；f1为外力；υ>0为流体的黏性系数；常数T>0为最终时刻。

对于多孔介质流体区域Ωp的流动过程由如下方程控制，

(2)

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

式中：nf、np为相应的Ωf、Ωp的单位外法向量；τi,i=1,2,…,d-1为交界面Γ的单位正切向量；g为重力加速度；α为取决于多孔介质的性质并且通过实验确定的正参数。

在给出耦合问题的弱形式之前，首先引入一些重要的记号。考虑函数空间,

Hf、Hp中元素的范数分别定义为

其中(· ,· )Ω为相应区域Ω内的内积。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 基于时间异步的稳定化特征线法

设τh为区域Ωf∪Ωp依赖于正参数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},稳定化项表达式定义为Gh(· ,·)，

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

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

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

则带有稳定化项的表示式为求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,

式中β为正的稳定化参数。

其中

定义投影算子，满足

∀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))有足够的光滑性，则有下列估计：

假定时间层是均匀的，对于多孔介质流区域Ωp的每一层时间Sk都对应自由流区域的相对应的时间层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进行循环，直到k=M- 1结束。

如何选取适当的r非常重要。r的选取是由流体的参数、边界条件和初始条件决定的，这就导致r的选取比较复杂，因此对于r值的选取是通过数值实验来确定的。

3 稳定性

为了更方便地分析稳定性,首先引入引理2、3。

(7)

(8)

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

(9)

在稳定性证明的这部分中，都假设满足

(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]内的稳定性，即

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))。

根据引理1的近似性质，可以得到

(23)

为了估计数值解与精确解的误差，对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)