鄧小毛 廖子菊

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

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

引言

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

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

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

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

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

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

1 控制方程

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

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

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

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

1.2 流體方程

ALE 坐標(biāo)下,不可壓縮Navier-Stokes 方程可寫(xiě)為

1.3 固體方程

對(duì)于固體變形,我們采用線彈性模型,在Lagrange 描述下的彈性動(dòng)力學(xué)方程為

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

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

2 數(shù)值方法

2.1 有限元離散

針對(duì)流固耦合方程,空間方向我們采用非結(jié)構(gòu)有限元進(jìn)行離散,為此我們首先給出控制方程的Galerkin 變分形式.

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

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

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

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

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

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

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

在時(shí)間方向,采用二階BDF 全隱格式進(jìn)行離散,得

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

在啟動(dòng)步中,先采用半步長(zhǎng)的隱式Euler 法計(jì)算出BDF2 格式所需要的初值.

2.2 Newton-Krylov-Schwarz 并行算法

在每一時(shí)間步,需要求解大規(guī)模非線性方程組(18).為了高效并行地進(jìn)行求解,本文構(gòu)造一種基于區(qū)域分解法的Newton-Krylov-Schwarz(NKS) 并行求解算法,主要包括非線性方程求解器、線性求解器和預(yù)條件子三部分.為了方便推導(dǎo),略去上標(biāo),將方程組寫(xiě)為一般形式

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

其中步長(zhǎng) α(k)采用線搜索方法[36]來(lái)計(jì)算,Newton 修正量 δy(k)通過(guò)求解如下線性方程組得到

其中 ηk為相對(duì)收斂誤差,當(dāng) ηk→0時(shí),則上述算法退化為精確Newton 法.

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

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

而限制型加性Schwarz 預(yù)條件子則定義為

整個(gè)NKS 算法的計(jì)算流程如下:

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

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

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

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

③通過(guò)線搜索獲得步長(zhǎng) αk;

④ 更新迭代解

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

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

3.1 算法驗(yàn)證

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

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

入口處給定入流速度

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

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

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

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

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

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

圖3 給出了t=3 s 時(shí)的彈性變形及流場(chǎng)速度分布圖,為了驗(yàn)證較大變形的情況,我們同時(shí)也給出了彈性模量Es=350kPa 及Es=87.5 kPa 的結(jié)果.而當(dāng)進(jìn)一步降低彈性體的彈性模量時(shí),由于彈性體變形過(guò)大,網(wǎng)格出現(xiàn)奇性,從而導(dǎo)致計(jì)算發(fā)散.注意這種發(fā)散不是數(shù)值格式不穩(wěn)定性帶來(lái),不能通過(guò)減小時(shí)間步長(zhǎng)來(lái)消除.

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

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

3.2 算法性能測(cè)試

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

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

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

表5 不同時(shí)間步長(zhǎng)對(duì)NKS 算法收斂性能的影響Table 5 Performance of NKS with respect toΔt

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

表6 關(guān)于流體黏性系數(shù)的算法魯棒性測(cè)試Table 6 Robustness of the algorithm with respect to the fluid viscosity

表7 關(guān)于固體密度的算法魯棒性測(cè)試Table 7 Robustness of the algorithm with respect to the solid density

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

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

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

表9 算法并行可擴(kuò)展性測(cè)試Table 9 Parallel performance and scalability of the algorithm

3.3 與其他算法性能的比較

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

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

4 結(jié)論

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