基于空間非迭代Oseen格式的時間解耦局部并行方法

張 蕊 張 翀 李 劍*

1(陜西科技大學文理學院 陜西 西安 710021)2(寶雞文理學院陜西省災害監測與機理模擬重點實驗室 陜西 寶雞 721013)

0 引 言

Navier-Stokes方程是描述水、大氣和其他流體運動的一般方程，是粘性不可壓縮流體的經典模型。而Navier-Stokes/Navier-Stokes耦合方程[1-2]是兩個Navier-Stokes方程在接口界面上通過耦合條件得到的復雜方程組，主要出現在氣象學和海洋學等理論研究中。復雜的自然現象不僅僅是單一的Navier-Stokes方程，因此耦合問題相關方法的研究也是必不可少的。

目前，Navier-Stokes/Navier-Stokes耦合問題的串行求解方法已有很多。例如：Connors等[3]運用傳統方法解決Navier-Stokes/Navier-Stokes耦合問題，并且分析了其穩定性與收斂性，在很大程度上保證了計算的準確性；文獻[5]通過運用牛頓格式的解耦方法分析流體耦合問題的穩定性與誤差估計，證明該解耦方法不僅無條件穩定而且長時間穩定；文獻[6]提出空間非迭代Oseen格式的歐拉時間推進方法，通過空間一步校正非線性項，僅用三層時間推進求解耦合問題。

耦合問題的數值模擬面臨兩大困難：較大的數據存儲空間和較長的計算時間。因此，需要借助并行計算技術實現快速求解。目前，國內外大批的科研人員與學者致力于單一Navier-Stokes方程并行有限元方法的研究，并取得了相關成果[7-13]。其中，單位劃分與自適應網格相結合的統一分區并行方法被稱作是一種靈活可控的區域分解方法[14-15]。Zheng等[7]提出了一種求解定常Navier-Stokes方程的局部并行有限元算法，通過劃分粗細網格，結合粗網格上的全局解與細網格上的局部解，最終得到有限元數值解。文獻[8]設計了兩種求解非定常Navier-Stokes方程的并行數值算法，一種是將Picard迭代法和區域分解算法相結合的并行Oseen線性化方法，另一種是基于兩重網格離散方法的有限元并行算法。而文獻[12]通過結合有限元方法與統一劃分方法，提出了一種求解時間依賴對流擴散方程的局部并行方法。

相比于單一Navier-Stokes方程的并行方法，耦合方程的并行有限元方法研究較少。因此，根據空間非迭代Oseen格式的歐拉時間推進方法與單元劃分的局部并行方法[7]，提出Navier-Stokes/Navier-Stokes耦合方程空間非迭代時間解耦的局部并行方法。在FreeFem++軟件的MPI并行功能[16]下，通過解耦復雜區域、處理邊界條件以及劃分并行區域，最終得到耦合模型的數值解。最后通過數值實驗驗證歐拉時間推進方法的精確性以及時間解耦局部并行方法的高效性。

1 預備知識

1.1 Navier-Stokes/Navier-Stokes耦合方程

本文主要考慮的Navier-Stokes/Navier-Stokes耦合方程為：

ui=0 onΓi=?ΩiI

(6)

式中：耦合區域Ω=Ω1∪Ω2由子區域Ω1、Ω2組成，其中Ωi?Rd,d=2，3，ni為單位外法向量,i=1,2;I為接口界面；粘性系數νi>0,i=1,2；κ為給定的參數；Δ表示拉普拉斯算子；速度ui:Ωi×[0,T]→Rd，其中T表示時間；壓力pi:Ωi×[0,T]→R；體積力fi:[0,T]→[H1(Ωi)]d。耦合區域如圖1所示。

圖1 耦合區域圖

為了解決耦合問題，提出如下兩個空間：

Xi={vi∈[H1(Ωi)]d:vi=0 onΓi;vi·ni=0 onI}

式(1)與式(4)的兩端分別乘以vi∈Xi與qi∈Mi，運用格林公式，得到式(1)-式(6)的變分形式，求解(ui,pi)∈Xi×Mi，使得(vi,qi)∈Xi×Mi，滿足：

b(ui,ui,vi)-d(vi,pi)=(fi,vi)

(7)

d(ui,qi)=0

(8)

雙線性型a(·,·)與d(·,·)定義在Xi×Xi與Xi×Mi上：

并且滿足連續性和inf-sup條件：

式中：C是常數；β>0。

三線性項b(·,·,·)定義在Xi×Xi×Xi上：

并且滿足如下反對稱性：

b(ui,vi,wi)=-b(ui,wi,vi)ui,vi,wi∈Xi

以及連續性：

一個實際的模型包括很多復雜項，本文僅討論耦合模型的算法問題。

1.2 歐拉時間推進方法

空間非迭代Oseen格式的歐拉時間推進方法主要運用Oseen格式一步校正非線性項，而耦合邊界式(2)則用如下的標準幾何平均跳躍[3]進行處理：

式中：Span表示擴張的空間。

并且滿足性質：

式中：c表示常數。

下面介紹空間非迭代Oseen格式的歐拉時間推進算法。

算法1空間非迭代的歐拉時間推進算法

(9)

(10)

2 時間解耦的局部并行方法

基于空間非迭代歐拉時間推進方法，給出Navier-Stokes/Navier-Stokes耦合方程時間解耦的局部并行方法。

首先，復雜區域Ω解耦為兩個單一子區域Ωi，解耦區域如圖2所示。因此，區域Ω內耦合方程的求解轉換為兩個子區域Ωi內單一方程的求解。

圖2 區域解耦圖

為了方便，定義:

(11)

根據耦合區域Ω的邊界條件式(2)、式(3)、式(6)，設置并行區域Ωi,k的不同邊界條件：

(1) 若并行區域為耦合邊界，設置式(2)為并行區域的邊界條件。

(2) 若并行區域為人工邊界，則由如下的第一類邊界條件代替：

(12)

(3) 若并行區域為混合邊界，結合耦合邊界式(2)與人工邊界式(12)進行并行計算。

因此，通過將耦合區域Ω上的方程轉換為并行區域Ωi,ki上方程的求解，根據并行區域邊界條件的設置，本文提出時間解耦的局部并行算法，從而獲得Navier-Stokes/Navier-Stokes耦合方程的近似解。

算法2時間解耦的局部并行方法

步驟3得到耦合方程的近似解：

3 數值實驗

為了驗證時間解耦局部并行方法的準確性，提出如下三個數值實驗。所有實驗使用的計算機CPU型號為Inter酷睿i7-8700 3.20 GHz，8 GB內存。算法1由軟件FREEFEM++完成，算法2的實現以FREEFEM++的MPI為平臺，網格采用三角形單元，有限元為P1b-P1。耦合方程的精確解詳見文獻[3,5]。

假設區域Ω1=[0,1]×[0,1],Ω2=[0,1]×[-1,0]外法向量n1=[0,-1]T,n2=[0,1]T，邊界I是x軸上0到1的位置，給定參數ν1=0.5,ν2=0.05,α=1,κ=100,T=1。并且速度ui與壓力pi的誤差表示為：

1) 給出串行方法的誤差與時間結果。固定時間步長Δt=0.01，改變子區域Ωi網格剖分的尺寸h，運用算法1求解Navier-Stokes/Navier-Stokes耦合方程，觀察數值解的誤差精度與求解時間，結果如表1-表2所示。

表1 歐拉時間推進方法的速度誤差結果

表2 歐拉時間推進方法的壓力誤差和時間結果

表1-表2驗證了Navier-Stokes/Navier-Stokes耦合方程的數值解中，速度解和壓力解的收斂階均達到了一階，但隨著網格剖分的增大，運行時間越來越長。

2) 為了說明空間非迭代Oseen格式的優越性，比較空間非迭代的Oseen格式與Oseen迭代格式、牛頓迭代格式、簡單迭代格式四種格式的求解時間。固定子區域Ωi的網格剖分h=1/16，時間步長Δt=0.01得到時間對比結果，如表3所示。

表3 四種格式的誤差和時間對比結果

續表3

從表3可以看出，在誤差精度方面，四種格式相差不大，但從運行時間來看，我們所選擇的空間非迭代Oseen格式的求解時間相對最短。

3) 為了評估時間解耦局部并行方法的高效性，分別運用算法1與算法2來計算耦合方程的有限元近似解，從而比較兩種方法的求解時間。

首先，分別給定網格剖分h=1/12,1/48，改變時間步長Δt=0.1,0.05,0.025,0.012 5，得到時間結果如表4-表5所示。

表5 網格剖分h=1/48的時間結果

通過表4-表5可以發現：與空間非迭代Oseen格式的歐拉推進方法相比，時間解耦局部并行算法節省的時間均在45%以上。這說明當網格剖分h不變時，改變時間步長Δt，并行算法的求解時間相對較短，即收斂速度較快。

然后，分別固定時間步長Δt=0.1,0.01，改變網格剖分h=1/12,1/24,1/48，得到時間結果如表6-表7所示。

表6 時間步長Δt=0.1的時間對比結果

表7 網格剖分Δt=0.01的時間對比結果

從表6-表7可以看出：與歐拉時間推進方法相比，當給定時間步長Δt時，時間解耦的局部并行方法至少節省38%以上的時間，而最大節省的時間達到70%。即改變區域網格剖分h，時間解耦局部并行方法的運行時間相對較短。

4 結 語

基于空間非迭代Oseen格式的歐拉時間推進方法，提出了求解Navier-Stokes/Navier-Stokes耦合方程時間解耦的局部并行方法。通過對復雜區域進行解耦，將耦合方程轉換為并行區域上的單一方程進行求解。改變子區域網格剖分的尺寸，從而確定并行區域的大小，在多個并行區域上同時計算。可以發現：不論是改變網格剖分的尺寸，還是時間步長的大小，時間解耦的局部并行算法在取得與原有串行算法相同精度的情況下，均有相對較快的收斂速度，即減少了Navier-Stokes/Navier-Stokes耦合方程的求解時間，在一定程度上體現了并行算法的高效性。