Navier-Stokes/Darcy模型的BDF2 模塊化梯度散度穩(wěn)定格式的數(shù)值分析

楊翠平, 王江珊, 賈宏恩

(1. 太原學(xué)院數(shù)學(xué)系，太原 030012; 2. 太原理工大學(xué)數(shù)學(xué)學(xué)院，晉中 030600)

0 引言

一直以來(lái)，對(duì)Navier-Stokes/Darcy 和Stokes/Darcy 模型數(shù)值方法的研究得到了廣泛的關(guān)注。迄今為止，大量的數(shù)值格式被提出用來(lái)求解這一耦合問(wèn)題。例如：有限元方法[1]、兩層網(wǎng)格法[2–4]、區(qū)域分解法[5–7]、邊界積分法[8–9]、特征線穩(wěn)定化有限元法[10]、梯度散度穩(wěn)定方法[11–12]。然而，梯度散度穩(wěn)定方法由于穩(wěn)定項(xiàng)導(dǎo)出的矩陣是奇異的，且對(duì)于大的穩(wěn)定項(xiàng)參數(shù)γ會(huì)導(dǎo)致求解器故障。為了減弱或克服這些困難，在文獻(xiàn)[13]中模塊化梯度散度穩(wěn)定方法被提出用于求解Navier-Stokes 方程，文獻(xiàn)[14]中對(duì)Navier-Stokes 方程提出了BDF2 模塊化梯度散度穩(wěn)定格式。本文對(duì)Navier-Stokes/Darcy 模型提出了BDF2 模塊化梯度散度穩(wěn)定格式，這種格式在保留梯度散度穩(wěn)定格式優(yōu)點(diǎn)的同時(shí)，更易于操作，且可以有效地避免大的穩(wěn)定化參數(shù)γ對(duì)解的非正常影響。

本文的結(jié)構(gòu)如下：第1 節(jié)，引入了Navier-Stokes/Darcy 模型和一些相關(guān)的符號(hào)、函數(shù)空間等；第2 節(jié)，給出了求解Navier-Stokes/Darcy 模型的BDF2 模塊化梯度散度穩(wěn)定格式并證明了該格式的穩(wěn)定性；第3 節(jié)，給出了誤差分析；第4 節(jié)，通過(guò)數(shù)值算例對(duì)于所提出方法的有效性進(jìn)行了驗(yàn)證。

1 Navier-Stokes/Darcy 模型及其相關(guān)符號(hào)

設(shè)在有界區(qū)域? ∈Rd(d= 2,3)，這里?包含兩個(gè)不重疊的流體區(qū)域，分別是自由流體區(qū)域?f和多孔介質(zhì)區(qū)域?p，見(jiàn)圖1。假定?f和?p的邊界滿足Lipschitz 條件，且

圖1 問(wèn)題區(qū)域? 及其邊界和交界面示意圖

自由流體區(qū)域?f由Navier-Stokes 方程控制

其中u=u(x,t)表示流體的速率，p=p(x,t)表示區(qū)域?f所受到的壓力，f1表示區(qū)域?f所受到的外力，ν>0 是?f流體的粘性系數(shù)。

多孔介質(zhì)區(qū)域?p由Darcy 方程控制

其中S0是物質(zhì)的釋水系數(shù)，up是區(qū)域?p中的流速，f2是源項(xiàng)，?表示測(cè)壓水頭函數(shù)，K 是滲透率張量。假設(shè)K=diag(K,··· ,K)，且滿足K>0。

通過(guò)使用Darcy 定律，(2)式可以整理為如下形式：

在交界面Γ上，我們考慮特定的交界面條件：質(zhì)量守恒條件、力的平衡條件和Beavers-Joseph-Saffman(BJS)條件

其中nf和np分別是?f和?p的單位外法向量，τi(i=1,2,··· ,d ?1)是交界面Γ上的單位正切向量，g是重力加速度，α取決于多孔介質(zhì)的性質(zhì)并由實(shí)驗(yàn)決定。

為便于以后的分析，我們?cè)谶吔绂、Γp上引入如下的邊界條件

2 全離散格式與穩(wěn)定性分析

估計(jì)不等式(22)右端的每一項(xiàng)如下：

最后，使用引理1，可得到(20)式成立。

3 誤差分析

定義投影算子[20]：Ph:(u(t),p(t),?(t))∈Hf×Q×Hp →(Phu(t),Php(t),Ph?(t))∈Hfh×Qh×Hph。假定(Phu(t),Php(t),Ph?(t))是(u(t),p(t),?(t))的近似，并且有以下性質(zhì)

記un、pn、?n為時(shí)間tn=n?t處的精確解，并假設(shè)精確解滿足如下的正則性

引入以下誤差

定義1 定義下述的誤差

成立。

定理2 在正則性假設(shè)條件(25)下，假定β>0，則存在C>0，滿足

(27)式、(28)式分別與(16)式、(17)式相減，有

合并兩式并使用引理5，有

接下來(lái)，利用Cauchy-Schwarz-Young 不等式和引理2，對(duì)(31)式右端的每一項(xiàng)進(jìn)行估計(jì)

使用引理2，有

對(duì)邊界項(xiàng)的估計(jì)如下

對(duì)非線性項(xiàng)的估計(jì)如下

4 數(shù)值算例

本節(jié)將利用數(shù)值算例來(lái)驗(yàn)證理論分析，區(qū)域?由?f= (0,1)×(1,2), ?p= (0,1)×(0,1)和邊界Γ=(0,1)×{1}組成，精確解為

實(shí)驗(yàn)中參數(shù)n、ρ、g、ν、K、S0和α都設(shè)為1，梯度散度參數(shù)γ= 1, β= 0.2，初始條件、邊界條件及穩(wěn)定項(xiàng)都滿足精確解。為了確保實(shí)驗(yàn)的準(zhǔn)確性，設(shè)定h= ?t，對(duì)Navier-Stokes 方程和Darcy 方程分別使用P2—P1 元和P1 元，數(shù)值實(shí)驗(yàn)結(jié)果參見(jiàn)表1 至表3。

表1 在T =1 時(shí)刻，具有L2 范數(shù)的標(biāo)準(zhǔn)格式的數(shù)值結(jié)果

表3 在T =1 時(shí)刻，具有L2 范數(shù)的BDF2 模塊化梯度散度穩(wěn)定格式的數(shù)值結(jié)果

來(lái)計(jì)算收斂階。通過(guò)觀察可以發(fā)現(xiàn)，梯度散度格式的速度的散度誤差要比標(biāo)準(zhǔn)格式的小，特別是BDF2 穩(wěn)定格式要比標(biāo)準(zhǔn)穩(wěn)定格式的更精確，所有的數(shù)值實(shí)驗(yàn)都證實(shí)了理論分析的正確性。

表2 在T =1 時(shí)刻，具有L2 范數(shù)的標(biāo)準(zhǔn)穩(wěn)定格式的數(shù)值結(jié)果