定常不可壓Navier-Stokes方程的并行有限元算法

尚月強，鄭 波，周康瑞，丁 琪

(西南大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，重慶 北碚 400715)

0 引言

Navier-Stokes(N-S)方程組是流體力學(xué)的基本方程組，出現(xiàn)于許多理論與應(yīng)用研究領(lǐng)域，如地球物理科學(xué)、水動力學(xué)、航天航空學(xué)、石油工業(yè)等，其研究對我國的經(jīng)濟與國防建設(shè)、工業(yè)設(shè)計具有非常重要的現(xiàn)實意義。然而，盡管對N-S方程的理論研究由來已久，數(shù)學(xué)上，三維N-S方程解的存在性與光滑性還沒解決，是千禧年七大數(shù)學(xué)難題之一[1]；物理上，許多基本問題也還未弄清楚(如：湍流的慣性區(qū)統(tǒng)計的間歇性修正問題，粘性與大尺度效應(yīng)等)[2]。這使得數(shù)值模擬成為一種十分重要的研究手段，N-S方程的數(shù)值模擬方法更是當(dāng)前科學(xué)與工程計算研究中的前沿?zé)衢T課題之一。

不可壓N-S方程的有限元數(shù)值模擬存在諸多困難：1)非線性，N-S方程是一種典型的非線性方程，數(shù)值求解中需采用合適的迭代格式進行線性化處理；2)不可壓縮性，N-S方程中的不可壓縮條件(連續(xù)方程)將速度變量與壓力變量耦合在一起，使得速度和壓力的有限元空間需滿足離散的inf-sup條件；3)大雷諾數(shù)問題，當(dāng)雷諾數(shù)較大時，N-S方程對流占優(yōu)，數(shù)值逼近中會出現(xiàn)偽震蕩；4)大規(guī)模計算問題，存在巨大的求解規(guī)模和有限的計算資源之間的矛盾。特別是高維復(fù)雜區(qū)域上N-S方程的數(shù)值模擬，不論是計算規(guī)模還是存儲需求，單機難以滿足其需求，因此，需借助并行計算技術(shù)，實現(xiàn)N-S方程的快速大規(guī)模并行數(shù)值模擬。然而，要實現(xiàn)并行計算，除了高性能并行機、一套完善的并行編程系統(tǒng)軟件外，高效的并行算法更是不可缺少的重要組成部分，它無論是對N-S方程的高效數(shù)值模擬，還是對提高并行機的運行效率，都起著至關(guān)重要的作用。

近年來，盡管并行計算技術(shù)取得了長足進步，但并行編程軟件遠(yuǎn)遠(yuǎn)落后于并行計算機硬件的發(fā)展，使得并行編程的難度遠(yuǎn)遠(yuǎn)大于串行編程，制約了并行計算的發(fā)展。因此，一個并行算法除了應(yīng)具有一個數(shù)值方法所應(yīng)具有的數(shù)學(xué)特征(如穩(wěn)定性、收斂性、有效性)外，還需實現(xiàn)簡單、通信需求少、具有良好的可擴展性，才具有廣闊的應(yīng)用前景，這在多核處理器普及的今天，顯得尤為重要。正是基于此認(rèn)識，一類簡單而高效的局部與并行有限元離散方法受到業(yè)界的普遍關(guān)注。該類方法最先由許進超教授和周愛輝教授針對非線性橢圓邊值問題提出[3]，其基本思想是基于對有限元解的局部與整體性質(zhì)的認(rèn)識(即解的整體性質(zhì)主要由低頻分量控制，可用粗網(wǎng)格較好地逼近, 其局部性質(zhì)則由高頻分量控制，可在細(xì)網(wǎng)格上局部和并行計算)，首先在一粗網(wǎng)格上求解一個非線性問題，得粗網(wǎng)格的有限元解，然后在相互重疊的細(xì)網(wǎng)格子區(qū)域上局部與并行求解一個線性化的殘差問題，以校正粗網(wǎng)格的有限元解。由于該類方法實現(xiàn)簡單、通信需求少、具有良好的并行性能，被廣泛應(yīng)用于求解各種問題，如特征值問題[4-6]、Stokes問題[7-11]、N-S問題[12-17]、磁動流體力學(xué)問題[18-20]、Stokes/Darcy流問題[21]，N-S/Darcy流問題[22]、四階微分方程[23]等；特別地，將該局部與并行有限元離散方法與區(qū)域分解算法和使用同階元的壓力穩(wěn)定化方法相結(jié)合，我們分別研究了Stokes方程的Schwarz區(qū)域分解方法[24]和并行穩(wěn)定化方法[25]，同時與完全重疊型區(qū)域分解技巧和求解大雷諾數(shù)問題的穩(wěn)定化方法相結(jié)合，研究了從小雷諾數(shù)到大雷諾數(shù)N-S方程的系列并行算法[26-36]。理論分析與數(shù)值試驗驗證了這類算法的高效性。

基于局部與并行有限元離散技巧，結(jié)合區(qū)域分解方法，本文主要討論N-S方程的簡單而高效的并行有限元算法。這些算法實現(xiàn)簡單，稍加修改現(xiàn)有的串行程序即可實現(xiàn)并行計算,可充分利用現(xiàn)有的串行軟件進行并行編程，通信需求少，能快速有效地模擬復(fù)雜的流體流動行為。我們將給出這些算法的描述、相關(guān)理論結(jié)果和數(shù)值算例，驗證算法的有效性。

本文內(nèi)容安排如下：第1節(jié)，給出N-S方程及其混合有限元數(shù)值逼近的相關(guān)知識；第2節(jié)，給出兩種并行策略；第3節(jié)，給出Dirichlet邊界條件N-S方程的完全重疊型區(qū)域分解并行算法；第4節(jié)，給出Dirichlet邊界條件N-S方程的并行兩水平有限元算法；第5節(jié)，討論非線性滑移邊界條件N-S方程并行有限元算法；第6節(jié)，結(jié)束語。

1 N-S方程及其有限元逼近

1.1 不可壓N-S方程

設(shè)Ω是Rn(n=2,3)中具有Lipschitz連續(xù)邊界的有界開集，我們考慮下列定常不可壓N-S方程:

(1)

我們引入Hilbert空間：

并定義a(u,v),b(u,v,w),d(v,q)如下：

a(u,v)=ν(▽u,▽v),d(v,q)=(▽·v,q),

本文用(·,·)表示L2(Ω)n(n=1,2,3)上的標(biāo)準(zhǔn)內(nèi)積,用‖·‖m表示Hm(Ω)n(m≥0)的范數(shù)。方程(1)的變分形式為:求(u,p)∈X×M, 滿足

a(u,v)+b(u,u,v)-d(v,p)+d(u,q)=(f,v), ?(v,q)∈X×M。

(2)

1.2 混合有限元逼近

設(shè)Tμ(Ω)={K}是Ω的一個尺寸為μ(0<μ<1,μ=H,h,H>h)的正則網(wǎng)格剖分，Xμ(Ω)?X,Mμ(Ω)?M是相應(yīng)于網(wǎng)格Tμ(Ω)的兩個有限元空間，

Xμ0(Ω) =Xμ(Ω)∩X,Mμ0(Ω) =Mμ(Ω)∩M。

給定G??Ω0?Ω(這里G??Ω0?Ω表示dist(?G?Ω0, ?Ω0?Ω)>0)，我們定義Xμ(G),Mμ(G)和Tμ(G)分別是Xμ(Ω),Mμ(Ω)和Tμ(Ω)在G上的限制，

關(guān)于混合有限元空間的假設(shè)如下：對任何(u,p)∈Hk+1(G)×Hk(G)(k≥1),存在(πμu,ρμp)∈Xμ(G)×Mμ(G),使得

‖μ-1(u-πμu)‖0,G+‖u-πμu‖1,G≤

cμs‖u‖s+1,G0≤s≤k,

‖μ-1(p-ρμp)‖-1,G+‖p-ρμp‖0,G≤

cμs‖p‖s,G0≤s≤k。

(3)

是標(biāo)準(zhǔn)的L2-投影：

1) 亞格子模型方法：

2) 基于兩個Gauss積分的變分多尺度方法：

2 并行策略

我們首先對求解區(qū)域Ω進行區(qū)域分解，得到互不重疊的子區(qū)域D1,D2,…,Dm，然后將Dj向外擴展一定尺寸得到相互重疊的子區(qū)域Ωj，滿足Dj??Ωj?Ω(j=1,2,…,m)。

2.1 完全重疊型區(qū)域分解的并行策略

在該并行策略中，每臺處理器根據(jù)流體在其所負(fù)責(zé)子區(qū)域上的物理特征，相互獨立地生成一個局部加密的全局多尺度網(wǎng)格(該網(wǎng)格絕大部分自由度來自其所負(fù)責(zé)的子區(qū)域)，然后在該全局多尺度網(wǎng)格上使用有限元方法計算其所負(fù)責(zé)子區(qū)域上的局部有限元解。具體說來，處理器j與其他處理器相互獨立地對其所負(fù)責(zé)的子區(qū)域Ωj生成尺度為h的細(xì)網(wǎng)格，對其余區(qū)域生成尺度為H>>h的粗網(wǎng)格，記所得的全局多尺度網(wǎng)格為Th,H(Ωj), 在這樣的全局多尺度網(wǎng)格Th,H(Ωj)上使用有限元方法求解N-S方程，計算結(jié)束時取子區(qū)域Dj上的計算結(jié)果，得到子區(qū)域Dj上的局部有限元解。由于每個子問題的網(wǎng)格都覆蓋整個區(qū)域，我們稱這種并行策略為完全重疊型區(qū)域分解。圖 1展示了兩個子區(qū)域情形的完全重疊型區(qū)域分解。

2.2 基于兩重網(wǎng)格離散的并行策略

在該并行策略中，我們分三步進行計算：1)在尺度為H的粗網(wǎng)格TH(Ω)上使用有限元方法求解非線性的N-S方程；2)使用區(qū)域分解技巧，在尺度為h的細(xì)網(wǎng)格Th(Ωj)上使用有限元方法并行求解線性化的局部殘差方程；3)在細(xì)網(wǎng)格互不重疊的子區(qū)域Dj上并行校正有限元解。如圖2所示：

圖1 完全重疊型區(qū)域分解Fig.1 A fully overlapping domain decomposition

圖2 基于兩重網(wǎng)格離散的并行策略Fig.2 A parallel strategy based on two-grid discretizations

3 完全重疊型區(qū)域分解并行算法

本節(jié)考慮Dirichlet邊界條件的N-S方程基于完全重疊型區(qū)域分解的并行算法，并分別考慮大雷諾數(shù)問題和使用同階元的并行穩(wěn)定化方法。

3.1 算法描述及理論結(jié)果

算法A 完全重疊型區(qū)域分解并行算法

備注3.1：

1)算法A的每一個子問題事實上是一個全局問題，但它用來求解在特定子區(qū)域的局部有限元解；

2)上述算法中，各處理器可自適應(yīng)地根據(jù)所模擬的流體在各自所負(fù)責(zé)子區(qū)域上的物理特征，相互獨立地生成恰當(dāng)?shù)木W(wǎng)格，選取恰當(dāng)?shù)挠邢拊瘮?shù)空間、穩(wěn)定化方法及代數(shù)方程組的求解器進行計算，具有很大的靈活性；

3)上述算法中，各處理器相互獨立地完成所有計算，無需交換數(shù)據(jù)，使得算法實現(xiàn)簡單，稍加修改現(xiàn)有的串行程序即可實現(xiàn)并行計算，具有良好的并行性能；

4)上述算法得到的有限元解時分片連續(xù)的，在整個區(qū)域上是間斷的。為了得到整體連續(xù)的有限元解，可在原算法的基礎(chǔ)上增加一后處理步，如文[3]及[37-38]的方法。

對于上述算法，使用inf-sup穩(wěn)定的混合有限元(即速度和壓力的有限元空間滿足離散的inf-sup條件，此時壓力穩(wěn)定項=0)，在一定的準(zhǔn)確解的正則條件下，我們有下列收斂性結(jié)果[34]：

|||▽(u-uh)|||0,Ω+|||p-ph|||0,Ω≤c(hs+Hs+1+α),1≤s≤t。

其中，α是速度穩(wěn)定項的系數(shù)，|||·|||0,Ω=

3.2 數(shù)值算例

算例1 已知準(zhǔn)確解[34]

針對已知準(zhǔn)確解的數(shù)值算例：

u1=10x2(x-1)2y(y-1)(2y-1),

u2=-10y2(y-1)2x(x-1)(2x-1),

p=x2-y2。

我們將求解區(qū)域Ω=[0,1]×[0,1]分解成4個子區(qū)域

D1=(0,1/2)×(0,1/2),

D2=(1/2,1)×(0,1/2),

D3=(1/2,1)×(1/2,1),

D4=(0,1/2)×(1/2,1),

并使用二階Taylor-Hood元進行了數(shù)值試驗，驗證了理論分析的正確性，數(shù)值結(jié)果如表 1所示。

表1 算法A所得近似解的誤差Tab.1 Errors of approximate solutions by Algorithm A

同時，針對不同的粘性系數(shù)，我們比較了使用不同的穩(wěn)定化方法的計算結(jié)果，即變分多尺度方法(Variational multiscale method), 亞格子穩(wěn)定化方法(Subgrid stabilized method)和不使用穩(wěn)定化方法(Method without stabilization)。表2表明，對于大粘性系數(shù)(即小雷諾數(shù))問題，相對于不使用穩(wěn)定化方法，算法中所使用的穩(wěn)定項沒有影響近似解的精度；而對于小粘性系數(shù)(即大雷諾數(shù))問題，穩(wěn)定化方法的使用非常必要(如表2中，當(dāng)ν=0.000 1時，不使用穩(wěn)定化方法根本得不到有限元解)。

表2 基于不同穩(wěn)定化方法的并行算法比較Tab.2 Comparison of numerical results by different parallel algorithms based on different stabilization methods

算例2后臺階流[34]

我們還將算法應(yīng)用于后臺階流的數(shù)值模擬，圖 3比較了雷諾數(shù)Re=800時我們的計算結(jié)果與文獻中的數(shù)據(jù)比較，圖4和圖5分別描繪了Re=800、1 000時計算所得的流線和壓力等值線。

算例3 方腔驅(qū)動流[36]

算法A中的問題是非線性問題，實際計算中需要采用迭代方法進行線性化處理。使用同階的P1-P1元和基于兩個Gauss積分的壓力穩(wěn)定化方法，我們考查了3種不同的迭代方法，即Newton迭代法、Oseen迭代法和Stokes迭代法, 得到了不同的穩(wěn)定性條件和相應(yīng)的收斂性結(jié)果，同時進行了系列數(shù)值試驗，進一步驗證了算法A的有效性。

使用不同的并行壓力穩(wěn)定化方法，即基于兩個Gauss積分的壓力穩(wěn)定化方法(Parallel Gauss method)、壓力投影的Petrov-Galerkin方法(Parallel PSPG method)、流線迎風(fēng)Petrov-Galerkin方法(Parallel SUPG method)、Galerkin最小二乘法(Parallel GLS method),我們將算法用于模擬二維方腔驅(qū)動流，圖6給出了雷諾數(shù)Re=400時所得數(shù)值結(jié)果與文獻中已有結(jié)果的比較，驗證了并行算法的有效性。

圖3 Re=800時后臺階流數(shù)值結(jié)果比較Fig.3 Comparison of numerical results for backward-facing step flow at Re=800

圖4 后臺階流的流線Fig.4 Streamlines of backward-facing step flow

圖5 后臺階流的壓力等值線Fig.5 Pressure contours of backward-facing step flow

圖6 方腔驅(qū)動流的并行壓力穩(wěn)定化方法比較Fig.6 Comparison of parallel pressure stabilization methods for lid-driven cavity flow

算例4 圓柱繞流[39]

我們使用同階的P2-P2元和基于壓力梯度投影的穩(wěn)定化方法，將算法A應(yīng)用于圓柱繞流的數(shù)值模擬。圖7和圖8顯示了我們的算法將求解區(qū)域分成3子區(qū)域時的計算結(jié)果與全局串行算法所得數(shù)值結(jié)果的比較，從中可以看出我們并行算法的有效性。

圖7 圓柱繞流的流線：(左)并行算法，(右)全局串行算法Fig.7 Streamlines of flow around a circular cylinder:(left) parallel algorithm,(right) global serial algorithm

圖8 圓柱繞流的壓力等值線：(左)并行算法，(右)全局串行算法Fig.8 Pressure contours of flow around a circular cylinder:(left) parallel algorithm,(right) global serial algorithm

4 并行兩水平有限元方法

4.1 算法描述與理論結(jié)果

將區(qū)域分解方法和兩重網(wǎng)格離散方法相結(jié)合，我們的并行兩水平有限元算法如下：

算法B1并行Oseen兩水平有限元方法

3)Dj(j=1,2,…,m)中并行校正：(uh,ph)=(uH,pH)+(ej,ηj)。

上述并行算法的第二步是基于Picard迭代線性化的殘差方程，我們也可采用其他線性化策略，得到類似算法：

算法B2 并行Newton兩水平有限元方法

1)同算法B1。

3)同算法B1。

算法B3 并行Stokes兩水平有限元方法

1)同算法B1。

3)同算法B1。

備注3.2：

2)并行算法B1-B3的第一步是在粗網(wǎng)格上求解，計算量很少，各處理器可相互獨立地進行，因此，在整個計算過程中，處理器可相互獨立地完成所有計算，從而可相互獨立地選取不同的有限元空間、不同的穩(wěn)定化方法、不同的線性化策略等，具有極大的靈活性；

3)與算法A類似，并行算法B1-B3中各處理器可相互獨立地完成所有計算，無需交換數(shù)據(jù)，使得算法實現(xiàn)簡單，稍加修改現(xiàn)有的串行程序即可實現(xiàn)并行計算，具有良好的并行性能。

使用inf-sup穩(wěn)定的混合有限元，對于算法B1-B3所得近似解，我們有下列誤差估計[35]

其中αH,αh中分別是粗細(xì)網(wǎng)格上關(guān)于速度穩(wěn)定項的系數(shù)。上述結(jié)果表明，當(dāng)

時，我們的并行算法能得到最優(yōu)的收斂階。

4.2 數(shù)值算例

算例1 解析解情形[35]

圖9給出了將算法B1-B3用于求解一個已知數(shù)值解算例，不同子區(qū)域個數(shù)時的收斂階，與理論結(jié)果一致。表3給出了與串行單水平相比，我們的并行算法使用四個子區(qū)域時所節(jié)約的計算時間。

表3 與串行單水平算法相比，并行算法所節(jié)約的時間Tab.3 Saved time by the parallel algorithms compared with the serial one-level algorithm

圖9 不同區(qū)域數(shù)的并行算法收斂階Fig.9 Convergence rates of parallel algorithms with different numbers of domains

圖10子區(qū)域之間的距離對近似解精度的影響，由此可看出并行算法的健壯性。

算例2 方腔動流[35]

圖11描述了將并行算法應(yīng)用于雷諾數(shù)Re=10 000時的二維方腔動流且子區(qū)域個數(shù)分別為4和16個時所得的流線。圖12則給出了Re=15 000和20 000時計算所得流線，這驗證了我們的并行算法在大雷諾數(shù)不可壓縮流數(shù)值模擬中的有效性。

圖10 子區(qū)域之間的距離對近似解精度的影響Fig.10 Influence of the distance between subdomains on the accuracy of approximate solutions

圖11 方腔驅(qū)動流Re=10 000時的流線：(左)4個子區(qū)域情形，(右)16個子區(qū)域情形 Fig.11 Streamlines of lid-driven cavity flow at Re=10 000:(left) 4 subdomains,(right) 16 subdomains

圖12 方腔驅(qū)動流的流線：(左)Re=15 000，(右)Re=20 000Fig.12 Streamlines of lid-driven cavity flow:(left) Re=15 000,(right) Re=20 000

5 非線性滑移邊界條件N-S方程的并行有限元算法

本節(jié)考慮帶非線性滑移邊界條件的N-S方程

(4)

非線性滑移邊界條件的N-S方程刻畫了醫(yī)學(xué)和環(huán)境學(xué)領(lǐng)域的某些流體流動現(xiàn)象，如動脈硬化患者的血液流動、油在沙礫上或沙礫下的流動等。我們分別研究了非線性滑移邊界條件的N-S方程并行兩水平有限元算法和完全重疊型區(qū)域分解并行算法，這里僅給出完全重疊型區(qū)域分解并行算法：

算法C 非線性滑移邊界條件N-S方程的完全重疊型區(qū)域分解并行算法

(Ω)(j=1,2,…,m), 使得

圖13 S=0.25時分叉血液流數(shù)值結(jié)果比較：(a,b)串行有限元方法，(c,d)并行算法[40]Fig.13 Comparison of numerical results for bifurcated blood flow with S=0.25:(a,b) serial finite element method,(c,d) parallel algorithm[40]

圖14 S=0.25時分叉血液流數(shù)值結(jié)果比較:(a,b) 串行有限元方法，(c,d)并行算法Fig.14 Comparison of numerical results for bifurcated blood flow with S=0.5:(a,b) serial finite element method,(c,d) parallel algorithm

6 結(jié)束語

隨著高性能并行機的發(fā)展和多核處理器的普及，高性能并行計算成為當(dāng)前科學(xué)計算的主流方向。并行算法在高性能并行計算扮演著非常重要的基礎(chǔ)作用，使得高性能并行算法的研究成為當(dāng)前科學(xué)計算領(lǐng)域的最前沿?zé)衢T課題之一。本文給出了數(shù)值求解定常不可壓N-S方程的若干并行算法，這些算法分別采用完全重疊型區(qū)域分解和并行兩重網(wǎng)格離散技巧，實現(xiàn)簡單，稍加修改現(xiàn)有的串行程序即可實現(xiàn)并行計算，通信需求少，具有良好的并行性能，能快速有效地模擬復(fù)雜的流體流動行為，具有較廣闊的應(yīng)用前景。下一步的工作是并行后處理方法及非定常問題的高效并行算法。