求解三維流固耦合問題的一種全隱全耦合區域分解并行算法1)

鄧小毛 廖子菊

*(廣東外語外貿大學數學與統計學院,廣州 510006)

?(暨南大學信息科學技術學院數學系,廣州 510632)

引言

流固耦合現象廣泛存在于自然界及各類工程問題中,如昆蟲與鳥類的飛行[1]、血流和血管壁的相互作用[2]、高層建筑、高聳結構及橋梁結構的風致振動[3-4]、飛行器的顫振與運動穩定性[5]以及醫療器械的設計與優化[6-7]等.在工程問題中,由于流場與彈性結構相互作用會誘發許多不良現象(如靜穩定性發散、顫振、極限環振蕩及渦激振動等),這些現象可能引起結構破壞或疲勞損傷,因此在工程設計與計算中考慮流固耦合作用對于結構安全具有重要意義.

流固耦合涉及到不同物理場的多尺度(時間和空間)非線性耦合,其快速準確求解一直以來都是力學與工程領域的一個難題[8].目前工程上流固耦合問題的求解策略主要分為兩類:分區算法和整體算法[9].在分區算法中,流體和固體區域分別獨立地進行求解并通過交界面交換位移和作用力等數據,這類方法具備良好的程序模塊性,可以直接利用現有的計算流體力學(CFD)和計算固體力學(CSD) 應用軟件及程序,且對計算機內存要求較低,是目前流固耦合問題的主流求解方法[10].

分區算法又分為分區顯式(弱耦合)算法[11-12]和分區隱式(強耦合)算法[13-14],其中分區顯式算法在每個時間步對流體和固體方程只分別求解一次,由于存在時間滯后效應,流固界面上的守恒律無法滿足,因此時間步長必須取得非常小,并且數值模擬結果有可能收斂到錯誤的解上.分區顯式算法主要適用于氣動彈性等較大質量比的情形[15].分區隱式算法則在每時間步內交替迭代求解各單場方程直至收斂以滿足界面守恒,常用的求解方法主要有固定點法[16]、Newton 類方法[17]和優化方法[18].分區隱式算法的這種交替迭代求解方式本質上是一種非線性Gauss-Seidel 型算法,它的收斂速度較慢,并且容易出現不收斂的情況,因此通常需要采用一些加速收斂技術如Aiken 松弛算法[19].此外,當流體和固體的密度非常接近時,迭代耦合方法會出現由附加質量效應導致的數值不穩定性[20].

為了避免分區算法的這些問題,近些年來整體算法逐漸受到越來越多的關注[21-25].整體算法又稱全耦合算法,其對流固控制方程以及界面守恒條件同時進行求解,可較好地適應強附加質量效應,以及不會惡化整體系統和各子系統的求解精度.此外,由于比分區算法少了一層迭代,收斂速度相對較快,且更適合并行計算,但其缺點是方程具有極強的非線性,且已失去了子問題控制方程系統矩陣的特點和可利用之處,不能直接利用現有的流體、固體算法和軟件,需要重新研究求解算法.文獻[26-27]對分區算法和整體算法的計算性能進行了詳細比較.

流固耦合問題整體求解的計算規模通常非常巨大,且離散矩陣具有病態性,因而對方程迭代求解技術要求很高,目前已經發展了一些有效的方法(主要是Newton 類方法)用于其離散非線性方程組及相應線性系統的求解[28-29],不過由于三維問題的計算量巨大,并且現有的算法大多缺乏較好的并行可擴展性,難以處理大規模問題的數值模擬,因此目前關于三維問題的數值結果尚不多.

本文將針對三維非定常流固耦合問題,構造一種適合大規模并行計算的整體求解并行算法.對于流固運動界面,我們采用任意拉格朗日-歐拉(arbitrary Lagrangian-Eulerian,ALE) 方法進行描述[30].對于流體、固體及動網格方程,采用非結構有限元進行統一離散,然后構造一種結合Newton 類非線性求解器,Krylov 子空間算法類線性求解器以及基于區域分解的Schwarz 預條件子[31]的Newton-Krylov-Schwarz(NKS) 算法進行高效求解.NKS 算法在求解大規模線性及非線性方程組時非常有效,已經成功地應用于各類大規模科學計算問題,我們在之前的工作中利用這類方法來構造求解流體動力學方程的并行算法[32-33],在數千核的計算規模中取得了良好的收斂性和并行可擴展性,本文將之拓展到流固耦合問題,通過標準測試算例驗證其計算效果,并對算法的收斂性、穩定性、魯棒性和并行可擴展性等數值性能進行測試分析.

1 控制方程

本文考慮不可壓縮流體與彈性體的相互作用問題.下面我們基于ALE 方法給出要求解的流固耦合控制方程及相應的邊界條件,對于流體、固體和網格等不同場的變量,我們分別用上標f,s,m來表示區分.

1.1 ALE 描述及動網格方程

其中dm(X,t) 表示網格位移,通過求解給定的網格運動方程得到.為了方便,這里我們采用調和方程來描述網格位移,即有

其中ds表示交界面處的固體位移.除了調和方程,也可以采用其他方程來描述網格位移,例如彈性方程,它對大變形的情況具有更好的適應性,不過與此同時也引入了一些附加的待定參數.我們的測試結果顯示,使用調和方程對于一般的變形已經可以得到良好的計算結果.

1.2 流體方程

ALE 坐標下,不可壓縮Navier-Stokes 方程可寫為

1.3 固體方程

對于固體變形,我們采用線彈性模型,在Lagrange 描述下的彈性動力學方程為

為了方便求解,將式(9)化為如下一階導數形式

式(1),式(4),式(5),式(10)統一構成流固耦合控制方程組,式(2)～式(3),式(6)～式(8),式(11)～式(12)組成邊界條件.

2 數值方法

2.1 有限元離散

針對流固耦合方程,空間方向我們采用非結構有限元進行離散,為此我們首先給出控制方程的Galerkin 變分形式.

對于動網格方程,分別定義如下試探函數空間和測試函數空間

注意流固界面處的位移連續性條件在Vm中作了限制.對于固體方程,類似地定義相應函數空間

注意上式左邊第三項利用到了流固界面處的應力連續性條件(12).

最后,對于流體方程,定義函數空間

對于變分式(13)～式(15),采用非結構有限元進行離散,其中網格方程和固體方程采用標準的線性有限元,流體方程采用P1-P1 穩定化有限元[34].通過定義P等函數空間相應的有限元空間,然后將式(13)～式(15)中的測試函數取為線性有限元的基函數,即可得到有限元離散方程.需要注意的是對于不可壓縮流動問題,當速度和壓力采用等階的P1-P1 有限單元時,其離散方程不滿足Ladyzhenskaya-Babu?ka-Breezi(LBB) inf-sup 條件,需要添加穩定化項.本文采用流線迎風Petrov-Galerkin(SUPG)方法[35],設流體區域的網格為={K},則流體方程的有限元變分形式改為

通過有限元離散后,我們得到如下半離散形式的非線性常微分方程組

其中y(t)為整個計算區域的所有待求未知量,即由網格位移,固體位移及速度,流體速度及壓力一起組成的向量,M為質量矩陣,L(y(t))為除時間導數項外的非線性函數項.注意由于流體區域的網格是隨時間變化的,因此這里質量矩陣也為y(t)的函數.

在時間方向,采用二階BDF 全隱格式進行離散,得

因此在每一時間步,需要求解如下大規模非線性方程組

在啟動步中,先采用半步長的隱式Euler 法計算出BDF2 格式所需要的初值.

2.2 Newton-Krylov-Schwarz 并行算法

在每一時間步,需要求解大規模非線性方程組(18).為了高效并行地進行求解,本文構造一種基于區域分解法的Newton-Krylov-Schwarz(NKS) 并行求解算法,主要包括非線性方程求解器、線性求解器和預條件子三部分.為了方便推導,略去上標,將方程組寫為一般形式

采用帶線搜索的非精確Newton 法來求解非線性方程組(19),假設第k迭代步的值為y(k),則第k+1 迭代步的值更新為

其中步長 α(k)采用線搜索方法[36]來計算,Newton 修正量 δy(k)通過求解如下線性方程組得到

其中 ηk為相對收斂誤差,當 ηk→0時,則上述算法退化為精確Newton 法.

對于Jacobian 矩陣Jk,一般可以通過無矩陣化方法[37]或者有限差分方法進行近似求解,也可以通過解析方法精確計算.在我們之前關于流體問題的測試中,采用精確的Jacobian 矩陣可以使算法具有更好的收斂性和魯棒性[32],因此,本文對Jacobian 矩陣也采用解析表達式進行計算.Jacobian 矩陣的主要計算量為流體方程的非線性對流項及穩定化項,我們注意到在這些項中對網格位移及速度的導數項的量級較小,因此進行了相應的簡化處理,即忽略對網格位移及速度的導數,這可以有效降低算法的計算量以及編程的工作量,并且數值測試表明這種處理幾乎不影響算法的收斂性.

最后,為了加快收斂速度,對線性方程組(21)進行預條件處理

而限制型加性Schwarz 預條件子則定義為

整個NKS 算法的計算流程如下:

(1) 采用上一個時間步的解作為初始值

(2) 對k=0,1,2,···,執行以下步驟,直至收斂:

①構造F(y(k)) 的Jacobian 矩陣Jk;

② 采用帶有預處理的Krylov 子空間迭代法GMRES 求解如下線性方程組

③通過線搜索獲得步長 αk;

④ 更新迭代解

3 數值算例與算法性能測試

本節首先采用標準測試算例對本文算法進行驗證,然后對算法的收斂性、魯棒性及并行可擴展性等數值性能進行測試研究.計算平臺為國家超級計算廣州中心的天河二號集群,其計算節點的硬件配置為雙路12 核Intel(R) Xeon(R) E5 CPU@ 2.4 GHz,64 G 內存.算法程序基于PETSc 軟件[38]編寫,計算網格采用非結構四面體網格,利用CUBIT 軟件[39]生成,網格并行分區采用ParMETIS 軟件[40]實現.計算過程中,Newton 法的收斂誤差限rtol設定為1.0×10-6,最大迭代次數設定為30;線性求解器采用重啟型GMRES,相對誤差限設為1.0×10-4,重啟次數設為400,最大迭代次數設為2000;默認情況下,區域分解的重疊層數設為2,子區域問題采用ILU(2)近似求解.

3.1 算法驗證

采用文獻[22,28]的彈性障礙物繞流標準測試計算模型,如圖2 所示,計算區域為(0,L)×(0,H)×,其中各幾何參數分別為:槽道的長L=1.5 m、寬W=0.8 m、高H=0.4 m,內部彈性阻礙物的長l=0.2 m、寬w=0.4 m、高h=0.2 m,彈性體與入口的距離為e=0.4 m ;流體的物理參數為ρf=1.0×103kg/m3,動力黏性系數為 μf=1.0Pa·s ;彈性體的物理參數為 ρs=1.0k g/m3,楊氏模量Es=1.4 MPa,Poisson 比 νs=0.4,相應的Lamé系數為 λs=2 MPa,μs=0.5 MPa.

圖2 計算區域Fig.2 Computational domain

入口處給定入流速度

其中s(t) 是一個時間方向的光滑因子,定義為

為平均流速,這里取=1.壁面采用無滑移邊界條件,出口處采用無應力條件n·σf=0;彈性體的底座處給定Dirichlet 邊界條件ds=us=0.

采用四個不同規模的網格進行計算,網格單元數與自由度總數如表1 所示.對于這四個不同網格,分別采用24,48,96,192 個處理器核進行計算,程序的運行內存,Newton 法的平均迭代次數,每Newton 步的平均GMRES 迭代次數,以及平均每時間步的計算時間也同時在表1 給出.

表1 計算網格Table 1 Computational meshes

為了驗證數值結果的準確性,我們計算位于彈性體頂點處的兩個觀測點P1(0.4,0.2,-0.2)及P2(0.5,0.2,-0.2)的位移,并與文獻[22]的計算結果進行比較,如表2,3 所示.可以看到本文計算結果與文獻結果基本一致,表明本文算法的求解結果是有效和可靠的.

表2 觀測點 P 1(0.4,0.2,-0.2) 處的位移Table 2 Displacements at pointP1(0.4,0.2,-0.2)

圖3 給出了t=3 s 時的彈性變形及流場速度分布圖,為了驗證較大變形的情況,我們同時也給出了彈性模量Es=350kPa 及Es=87.5 kPa 的結果.而當進一步降低彈性體的彈性模量時,由于彈性體變形過大,網格出現奇性,從而導致計算發散.注意這種發散不是數值格式不穩定性帶來,不能通過減小時間步長來消除.

圖3 不同楊氏模量下的彈性變形及速度分布圖Fig.3 Deformation and velocity distribution with different Young's modulus

表3 觀測點 P 2(0.5,0.2,-0.2) 處的位移Table 3 Displacements at pointP2(0.5,0.2,-0.2)

3.2 算法性能測試

首先對算法的收斂性能進行測試.影響算法收斂性能的參數主要有兩個:一是區域分解的重疊層數 δ;二是構造子區域局部預條件子的線性求解器.分別取 δ=1,2,然后采用不同填充級的子區域求解器ILU(k) 進行計算,測試結果如表4 所示,其中Newton 表示平均每個時間步的非線性迭代次數,GMRES 表示平均每個Newton 步的線性迭代次數,Time 表示平均每個時間步的計算時間(單位為秒).可以看到,Newton 迭代次數基本不受影響,表明不同填充級的ILU 線性求解器均可作為子區域問題的求解器.另外,由于填充級增加時,ILU對Jacobian 矩陣的逆有更好的逼近,因此可以看到線性迭代次數顯著隨之減少.不過每時間步的計算時間并沒有相應減少,原因是當填充級增加時,雖然線性迭代次數顯著減少了,但是花在子區域預條件子構造的時間相應增加了,因此有個折衷的選擇,在本算例中,填充級取為k≤2 均可,當k=3 時計算時間則顯著增加.另外,由于本算例中Newton法的收斂速度很快,因此子區域的重疊層數對計算時間影響不大.

表4 子區域的重疊層數及子問題求解器對算法性能的影響Table 4 Impact of the overlapping size and the subsolver on the NKS algorithm

接下來我們考察NKS 算法的穩定性.表5 給出取不同時間步長時算法的收斂性能.由于時間方向采用了全隱離散格式,因此算法具有很好的穩定性,當時間步長Δt從0.001 增大至0.128 時,雖然Newton迭代步數、線性迭代步數和每步的計算時間都略有增加,但總的來說算法仍能快速收斂,表明算法可以在精度允許的范圍內自由地選擇時間步長進行數值模擬.

表5 不同時間步長對NKS 算法收斂性能的影響Table 5 Performance of NKS with respect toΔt

隨后對NKS 算法的魯棒性進行測試,即考察取不同的物理參數對算法收斂性能的影響.表6給出了僅改變流體黏性系數,而其他參數保持不變時算法收斂性能的測試結果,表7 給出僅改變固體密度時的測試結果,表8 則給出僅改變固體的楊氏模量和Poisson 比的測試結果.可以看到在不同的物理參數下,算法的收斂性能基本保持不變,表明了算法具有良好的魯棒性.

表6 關于流體黏性系數的算法魯棒性測試Table 6 Robustness of the algorithm with respect to the fluid viscosity

表7 關于固體密度的算法魯棒性測試Table 7 Robustness of the algorithm with respect to the solid density

表8 關于楊氏模量和Poisson 比的算法魯棒性測試Table 8 Robustness of the algorithm with respect to the Young’s modulus and Poisson’s ratio

最后,對算法的并行性能進行測試,如表9所示,其中np表示計算的進程數(等于處理器的核數),speedup為加速比,ideal 表示理想加速比,efficiency 表示并行效率.可以看到,隨著進程數的增加,非線性迭代次數保持不變,線性迭代次數則略有增加,每迭代步的計算時間則顯著減少.當進程數從192 增加到1536 時,算法取得了超線性的加速比,當進程數增加到3072 時,并行效率為91%.算法的加速比曲線如圖4 所示,測試結果表明了本算法具有良好的并行效率和可擴展性.

圖4 算法加速比Fig.4 Parallel speedup of the algorithm

表9 算法并行可擴展性測試Table 9 Parallel performance and scalability of the algorithm

3.3 與其他算法性能的比較

本節將本文方法與文獻中已有的一些算法進行比較.對不同算法的性能直接進行對比是一項復雜的工作,一方面由于不同文獻通常采用了不同的測試算例,另一方面算例的計算規模、算法的參數設置以及進行測試的機器也不相同.本文選取平均每時間步的非線性迭代次數作為算法收斂性能的考量,另外,基于不同文獻測試算例的問題規模不同,參照文獻[29]的做法,我們定義一個時間步內“每CPU 核在每秒鐘內所求解的自由度數”這一指標來比較不同算法的計算效率.最后,為了考察算法處理大規模問題的能力,對算法的并行性能也進行對比.關于全耦合算法,目前文獻中的測試算例大多都是關于二維問題的,三維算例較少,在我們的檢索范圍內,找到了文獻[22,28-29]這幾種不同的數值方法,這幾種方法均采用Newton 法作為非線性求解器,Krylov 子空間迭代法作為線性求解器,它們的主要區別在于預條件算子及其求解算法不同,表10給出了本文方法與它們在收斂性能與并行性能方面的比較結果.可以看到,在不考慮算法收斂參數設定以及測試機器的情況下,本文方法每時間步的Newton 迭代次數最少,每CPU 核在每秒鐘內所求解的自由度數僅次于文獻[28]的方法,在并行可擴展性與并行效率上,本文方法遠遠高于其他方法,表明在處理大規模問題上具有較好的潛能.

表10 本文算法與其他文獻算法的性能比較Table 10 Comparison of different numerical approaches found in literature

4 結論

本文針對三維流固耦合問題,提出了一種有限元數值求解的全隱全耦合可擴展并行算法.為方便處理復雜的計算區域,采用了三維非結構網格進行空間離散,為使算法具有更好穩定性和魯棒性,采用二階BDF 全隱格式進行時間離散.對于時空離散后得到的大規模非線性代數系統,構造了基于區域分解的Newton-Krylov-Schwarz 并行算法對流體、固體與動網格方程進行一次性整體求解.采用彈性障礙物繞流問題的標準測試算例對算法的準確性進行了驗證.算法性能測試結果表明本文提出的算法在數千處理器核條件下具有良好的并行可擴展性能.通過測試算法在不同計算時間步長、不同物理參數下的收斂性能,驗證了算法具有良好的穩定性和魯棒性.另外,本文對動網格方程采用了較為簡單的調和方程模型,數值結果表明其在大變形的情況下計算會發散,為了避免這一問題,可以考慮更復雜的彈性方程模型,這些問題有待進一步深入研究.