鄧小毛 廖子菊

*(廣東外語外貿(mào)大學數(shù)學與統(tǒng)計學院,廣州 510006)

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

引言

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

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

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

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

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

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

1 控制方程

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

1.1 ALE 描述及動網(wǎng)格方程

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

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

1.2 流體方程

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

1.3 固體方程

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

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

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

2 數(shù)值方法

2.1 有限元離散

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

對于動網(wǎng)格方程,分別定義如下試探函數(shù)空間和測試函數(shù)空間

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

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

最后,對于流體方程,定義函數(shù)空間

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

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

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

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

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

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

2.2 Newton-Krylov-Schwarz 并行算法

在每一時間步,需要求解大規(guī)模非線性方程組(18).為了高效并行地進行求解,本文構造一種基于區(qū)域分解法的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 矩陣的主要計算量為流體方程的非線性對流項及穩(wěn)定化項,我們注意到在這些項中對網(wǎng)格位移及速度的導數(shù)項的量級較小,因此進行了相應的簡化處理,即忽略對網(wǎng)格位移及速度的導數(shù),這可以有效降低算法的計算量以及編程的工作量,并且數(shù)值測試表明這種處理幾乎不影響算法的收斂性.

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

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

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

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

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

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

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

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

④ 更新迭代解

3 數(shù)值算例與算法性能測試

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

3.1 算法驗證

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

圖2 計算區(qū)域Fig.2 Computational domain

入口處給定入流速度

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

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

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

表1 計算網(wǎng)格Table 1 Computational meshes

為了驗證數(shù)值結果的準確性,我們計算位于彈性體頂點處的兩個觀測點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 的結果.而當進一步降低彈性體的彈性模量時,由于彈性體變形過大,網(wǎng)格出現(xiàn)奇性,從而導致計算發(fā)散.注意這種發(fā)散不是數(shù)值格式不穩(wěn)定性帶來,不能通過減小時間步長來消除.

圖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 算法性能測試

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

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

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

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

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

表6 關于流體黏性系數(shù)的算法魯棒性測試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表示計算的進程數(shù)(等于處理器的核數(shù)),speedup為加速比,ideal 表示理想加速比,efficiency 表示并行效率.可以看到,隨著進程數(shù)的增加,非線性迭代次數(shù)保持不變,線性迭代次數(shù)則略有增加,每迭代步的計算時間則顯著減少.當進程數(shù)從192 增加到1536 時,算法取得了超線性的加速比,當進程數(shù)增加到3072 時,并行效率為91%.算法的加速比曲線如圖4 所示,測試結果表明了本算法具有良好的并行效率和可擴展性.

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

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

3.3 與其他算法性能的比較

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

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

4 結論

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