黨雷寧, 白智勇, 柳 森

(中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所, 四川 綿陽 621000)

0 引 言

近幾十年來，計算流體力學(xué)(CFD)越來越廣泛地應(yīng)用于航空航天，可以說CFD改變了航空航天的設(shè)計過程。計算機性能的提高和CFD技術(shù)的發(fā)展，是CFD大規(guī)模應(yīng)用的“助推器”。據(jù)統(tǒng)計，高性能計算機(HPC)計算速度每十年將提高三個量級[1]。然而，現(xiàn)今的CFD方法在建立過程中很多都沒有考慮到并行計算的特點，并行水平較低，已經(jīng)不能適應(yīng)HPC發(fā)展的要求。NASA在2030年的CFD展望[2]中指出，目前魯棒性的CFD流動解算器非常缺乏可擴展性，極少有可以有效利用超過O(1000)個內(nèi)核的應(yīng)用能力，為此需要開發(fā)并行程度很高的CFD算法，用以平衡計算與通訊。

在CFD并行計算中，普遍采用基于網(wǎng)格分區(qū)的并行策略。在CFD時間推進中，為了不受穩(wěn)定性條件的限制以及出于加速流動模擬收斂的目的，一般采用各種隱式計算方法。然而，傳統(tǒng)的以LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法[3]為代表的隱式近似因子分解(AF)類方法應(yīng)用于并行計算時，存在兩個問題。一方面，隱式算法中相鄰網(wǎng)格單元之間較強的數(shù)據(jù)依賴使得并行求解較為困難[4]。另一方面，如果直接在各網(wǎng)格分區(qū)中采用隱式算法，則隨著分區(qū)數(shù)的增加算法收斂性能將惡化[5]。

線松弛方法是一類隱式求解的迭代算法。MacCormack(1985)[6]提出求解Navier-Stokes方程的GSLR(Gauss-Seidel Line Relaxation)方法。Candler(1988)[7]將其應(yīng)用到熱化學(xué)非平衡流動的數(shù)值模擬中，它的特點是收斂很快。GSLR方法的思想是：在黏性流動中，相較于流向和周向，流場梯度在物面法向更大，流動在物面法向耦合更緊密，數(shù)值求解時在物面法向把隱式部分的對角項和非對角項耦合求解才能反映流動的物理本質(zhì)。GSLR方法在流向和周向采用Gauss-Seidel掃描，不利于并行計算。NASA Ames研究中心的Wright(1997)[8]，以GSLR為出發(fā)點，把Gauss-Seidel掃描改為松弛過程，消除相鄰物面法向網(wǎng)格線流場變量之間的數(shù)據(jù)依賴，既保留了GSLR收斂快的特點，又使它有利于并行計算，這就是DPLR(Data-Parallel Line Relaxation)方法。

Wright發(fā)明并在單塊結(jié)構(gòu)網(wǎng)格上應(yīng)用DPLR方法[8]，后以此為基礎(chǔ)研制著名的氣動熱力學(xué)軟件DPLR[9-10]，這是一款基于多塊結(jié)構(gòu)網(wǎng)格的CFD軟件。但是DPLR方法在多塊網(wǎng)格上的性能如何，網(wǎng)格分區(qū)對算法的收斂性能會產(chǎn)生多大影響，這些問題未在公開文獻中報道。由于DPLR需要結(jié)合隱式邊界條件使用以及多塊對接結(jié)構(gòu)網(wǎng)格中存在大量的對接邊界，本文首先提出一種隱式對接邊界條件的處理方法，為把DPLR應(yīng)用于多塊結(jié)構(gòu)網(wǎng)格奠定基礎(chǔ)。然后，以二維和三維球頭繞流為算例，分析比較了DPLR與LU-SGS的收斂速度，研究了網(wǎng)格分塊方式對算法穩(wěn)定性和收斂性能的影響。最后進行了返回艙外形繞流的數(shù)值模擬，給出返回艙外形在馬赫數(shù)9.96、雷諾數(shù)2.04×105、迎角10°～30°的氣動力特性，并與風(fēng)洞試驗結(jié)果對比，考察算法在復(fù)雜外形飛行器繞流數(shù)值模擬中的性能。

1 數(shù)值計算方法

1.1 控制方程

為敘述基于多塊對接結(jié)構(gòu)網(wǎng)格的DPLR隱格式實現(xiàn)途徑，以二維無量綱Navier-Stokes方程為例，說明本文的計算方法，三維控制方程的計算方法與之類似?？刂品匠蘙11]可寫為：

式中，Q為守恒變量矢量，E、F為無黏通量矢量，Ev、Fv為黏性通量矢量，它們的表達式詳見文獻[11]。

1.2 空間項離散

應(yīng)用有限體積法，以網(wǎng)格單元為控制體，對控制方程進行離散。在有限體積法中，流場變量儲存在單元格心，在單元內(nèi)部均勻分布。本文記流向為ξ，下標索引為i；法向為η，下標索引為j。把控制方程(1)在控制體i,j積分，并應(yīng)用Gauss-Green定理，則有：

式中，RHS表示殘差，E、F分別為通過垂直于ξ、η方向單元面的無黏通量，Ev、Fv分別為相應(yīng)的黏性通量。V是控制體單元的體積，S是單元面的面積。

在本文中，無黏通量采用Steger-Warming矢通量分裂[12]以及NND格式[13]計算；黏性通量采用中心差分格式計算。

1.3 時間項離散

對方程(2)取n+1時間層，左側(cè)的時間導(dǎo)數(shù)進行后向差分，在離散中需要對無黏和黏性通量線性化。無黏通量E的線性化為:

En+1≈En+AnδQn+1(4)

黏性通量Ev根據(jù)Tysinger和Caughey的方法[14]線性化為:

其中:

A±、B±分別是ξ、η方向的經(jīng)過分裂的無黏Jacobian矩陣，L和N分別是ξ、η方向黏性通量的線性化矩陣，Δt是時間步長。

1.4 DPLR格式

設(shè)η方向為物面法向，則DPLR格式[8]可寫為:

對于m=1,mmax

1.5 邊界條件

DPLR需要和隱式邊界條件結(jié)合起來使用[8-9]。本文采用虛擬網(wǎng)格方法[6-7,16]處理邊界條件。虛擬網(wǎng)格不由網(wǎng)格軟件生成，而是在代碼中根據(jù)不同邊界類型邏輯生成。圖1和圖2分別是物面邊界和對接邊界虛擬網(wǎng)格示意圖。

圖1 物面邊界虛擬網(wǎng)格示意圖Fig.1 Sketch of ghost cell at wall

圖2 對接邊界虛擬網(wǎng)格示意圖Fig.2 Sketch of ghost cell at interface boundary

文獻[16]認為，虛擬網(wǎng)格上流場變量設(shè)置的原則是邊界面上的無黏和黏性通量能夠正確計算。在GSLR方法中，MacCormack認為[6]可以在隱式計算中把邊界條件包含進來，這就是隱式邊界條件。因此，DPLR格式的邊界條件分為顯式和隱式兩部分，顯式部分是為了正確計算邊界面上的通量，隱式部分則體現(xiàn)了邊界條件對隱式求解的影響。這和LU-SGS方法是不同的，LU-SGS隱式部分通常不包含邊界條件。

參照圖1，設(shè)邊界面兩側(cè)網(wǎng)格中心的某流場變量分別為φ1和φ0,設(shè)邊界面上的物理量為φb，下標“1”表示內(nèi)點，“0”表示虛擬點，“b”表示邊界面。顯式邊界條件可寫為通式:

φ0=2φb-φ1(10)

對于整個求解區(qū)域的外邊界，顯示邊界條件可應(yīng)用式(10)。例如，對于超聲速入口，直接給出φ0=φ∞，其中“∞”表示自由來流；對于超聲速出口，φb=φ1；對于絕熱壁無滑移物面，Pb=P1，ub=0，vb=0，Tb=T1；對于等溫壁無滑移物面，Pb=P1，ub=0，vb=0，Tb=Tw。其中，P、u、v、T分別表示流場壓力、x方向速度分量、y方向速度分量和溫度。對于相鄰網(wǎng)格塊之間的對接邊界(見圖2)，虛擬網(wǎng)格上的流場變量φ0從與之對接的網(wǎng)格塊內(nèi)點直接得到，φ0和φ1之間沒有類似式(10)的函數(shù)關(guān)系。

對式(10)微分，得到dφ0和dφ1的函數(shù)關(guān)系，經(jīng)過進一步推導(dǎo)可得到邊界兩側(cè)守恒變量微分之間的函數(shù)關(guān)系:

dQ0=HdQ1(11)

其中H在文獻[16]中稱為折疊矩陣，式(11)就是隱式邊界條件。把式(11)代入離散方程(9)，就可得到靠近邊界內(nèi)點的離散方程。例如，設(shè)某邊界位于j=1/2位置，則在控制體(i,1)處的離散方程為:

從式(12)可看出，隱式邊界條件把邊界條件的影響包含進離散方程的隱式部分。

但是對于對接邊界，邊界面兩側(cè)的流場變量φ0和φ1不存在顯式的函數(shù)關(guān)系式(10)，也無法推導(dǎo)出隱式邊界條件式(11)，那么在對接邊界處怎樣把邊界條件對隱式求解的影響考慮進來呢？本文提出一種處理方法，還是以j=1/2邊界(假定是對接邊界)為例：

表1 算例驗證的計算條件Table 1 Computational condition of two cases

2 算例驗證及討論

2.1 計算條件及網(wǎng)格

為了驗證本文的計算方法和程序、比較DPLR和LU-SGS的收斂性能以及研究網(wǎng)格分區(qū)方式對算法穩(wěn)定性和收斂性能的影響，本文計算兩個算例。算例1是二維球頭繞流，算例2是三維球頭繞流，計算條件如表1所示。算例1的網(wǎng)格如圖3(a)所示，網(wǎng)格數(shù)為121×81(流向×法向)。算例2的網(wǎng)格如圖3(b)所示，總共3個網(wǎng)格塊，每個塊的網(wǎng)格數(shù)都是16×101×16(流向×法向×周向)。兩個算例法向最小網(wǎng)格間距都設(shè)置為使得網(wǎng)格雷諾數(shù)ReΔn=2。

(a) 二維 (b) 三維

2.2 計算結(jié)果正確性驗證

對于二維球頭繞流，本文把計算得到的物面無量綱壓力和熱流(與駐點值的比值)分布與Wieting試驗結(jié)果[17]比較(見圖4)，兩者吻合較好。對于三維球頭繞流，本文把計算結(jié)果和Coleman的試驗結(jié)果[18]比較，并如圖5所示。由圖可見，本文計算的三維球頭物面壓力(與自由來流值的比值)、熱流分布與試驗吻合較好。

(a) 壓力曲線 (b) 熱流曲線

(a) 壓力曲線 (b)熱流曲線

2.3 DPLR穩(wěn)定性和收斂性能以及網(wǎng)格分區(qū)方式影響的分析與討論

收斂快是DPLR算法的特點之一，Wright在這方面已經(jīng)做了詳盡的研究[8]。對于二維球頭和三維球頭算例，本文把DPLR與LU-SGS算法的收斂曲線進行對比，分別如圖6和圖7所示，其中橫坐標是CPU時間(Intel Core i7-4790 單核)。由于DPLR不是無條件穩(wěn)定格式[8]，在應(yīng)用DPLR格式中這兩個算例的CFL都取為能夠使計算正常進行并保證收斂的最大CFL數(shù)，其中二維球頭繞流算例CFL數(shù)取為3000，三維球頭繞流算例CFL數(shù)取為3500。LU-SGS是無條件穩(wěn)定格式，其CFL在兩個算例中都取為1.0×106。此外，基于Wright的經(jīng)驗[8,10]在DPLR計算時采用全局時間步長法；而在LU-SGS計算時應(yīng)用當(dāng)?shù)貢r間步長法以加速收斂。從DPLR和LU-SGS收斂性能的對比看，DPLR格式收斂更快，二維球頭繞流收斂到機器零的時間相當(dāng)于LU-SGS的35%，三維球頭繞流收斂到機器零的時間相當(dāng)于LU-SGS的22%。從殘差曲線的變化歷程看，在迭代的初始階段，LU-SGS方法由于單步計算量小而收斂較快，然而在迭代后期，DPLR方法就顯示出優(yōu)勢。這啟發(fā)我們，可以在迭代初期使用LU-SGS而在后期使用DPLR方法，這樣既能利用LU-SGS方法單步計算量小以及穩(wěn)定性好的特點，在時間推進的初始階段更好地避免強烈非定常特性給收斂帶來的不利影響并使殘差快速下降，又能發(fā)揮DPLR格式在迭代后期的優(yōu)勢。

圖6 二維球頭繞流DPLR與LU-SGS收斂曲線的比較Fig.6 Comparison between DPLR and LU-SGS residual curve for flowfield computation of 2D circular cylinder

圖7 三維球頭繞流DPLR與LU-SGS收斂曲線的比較Fig.7 Comparison between DPLR and LU-SGS residual curve for flowfield computation of 3D sphere

為了考察多塊網(wǎng)格下DPLR方法的性能。本文把二維球頭繞流網(wǎng)格和三維球頭繞流網(wǎng)格進行分塊。有兩種分塊方式，一是在流向和周向分塊；二是在法向分塊。

圖8是對初始二維球頭繞流網(wǎng)格在流向進行等間距分塊、得到的不同網(wǎng)格分塊數(shù)流場的殘差收斂曲線。由圖可見，不同網(wǎng)格分塊數(shù)流場的殘差收斂曲線基本重合，在流向進行網(wǎng)格分塊對DPLR算法的收斂性能幾乎沒有影響。圖9是對三維球頭繞流網(wǎng)格在流向和周向進行等間距分塊、得到的不同網(wǎng)格分塊數(shù)流場的殘差收斂曲線，由圖可以得到和二維球頭繞流同樣的結(jié)論。這很容易理解，參看式(9)，在第m松弛步，DPLR算法在每條法向網(wǎng)格線上求解，而流向網(wǎng)格上的δQ已經(jīng)在m-1松弛步得到計算，這就是說每條法向網(wǎng)格線上的隱式求解都是獨立的、不會相互影響，因此在流向和周向進行網(wǎng)格分塊不會影響算法的收斂過程。

下面考察在法向進行網(wǎng)格分塊對算法性能的影響。圖10是對二維球頭繞流網(wǎng)格在法向進行分塊的分塊位置示意圖，本文分別考察在兩個分塊位置分塊后形成的多塊網(wǎng)格流場的收斂性能，其中分塊位置1(splitting 1)遠離邊界層，分塊位置2(splitting 2)位于邊界層內(nèi)。圖11展示了二維球頭繞流網(wǎng)格在法向分塊后得到的多塊網(wǎng)格流場收斂曲線，由圖可見，在分塊位置1進行網(wǎng)格分塊對收斂性能影響較小，但在分塊位置2進行網(wǎng)格分塊降低了算法的收斂性能。分析原因，相較于分塊位置1，分塊位置2在邊界層內(nèi)，流場法向梯度大，相鄰法向網(wǎng)格單元之間流場變量耦合緊密，而DPLR算法在分塊后的兩個塊上分別求解使得兩個網(wǎng)格塊流場變量聯(lián)系得不那么緊密，收斂速度下降是必然的。需要說明的是，DPLR隱格式不是無條件穩(wěn)定的[8]，在初始網(wǎng)格求解時CFL數(shù)可以取到3000(圖6，圖8)，在分塊位置2分塊后CFL數(shù)最大只能取到1800(圖11)，這說明分塊后格式穩(wěn)定性下降。圖11中網(wǎng)格分塊后流場的DPLR計算，采用本文提出的隱式對接邊界條件(見式(13)，圖中標注為“implicit interface boundary condition”)，而不采用隱式對接邊界條件即直接舍棄與虛擬網(wǎng)格點有關(guān)的項，算法CFL數(shù)最大可取到100，收斂很慢，在圖中未給出。

圖8 二維球頭繞流網(wǎng)格在流向分塊時的殘差收斂曲線(CFL=3000)Fig.8 Residual history of flow around 2D circular cylinder when splitting grid blocks in streamwise direction(CFL=3000)

圖9 三維球頭繞流網(wǎng)格在流向和周向分塊時的殘差收斂曲線(CFL=3500)Fig.9 Residual history of flow around 3D sphere when splitting grid blocks in streamwise and circumferential direction(CFL=3500)

圖10 二維球頭繞流網(wǎng)格法向分塊位置示意圖Fig.10 Sketch of block splitting position in the grid of 2D circular cylinder

圖11 二維球頭繞流網(wǎng)格在法向分塊時的殘差收斂曲線(CFL=1800)Fig.11 Residual history of flow around 2D circular cylinder when splitting grid blocks in normal direction(CFL=1800)

圖12是三維球頭繞流網(wǎng)格在法向分塊后得到的多塊網(wǎng)格流場收斂曲線，圖中的兩個分塊位置(splitting 1和splitting 2)與二維球頭繞流網(wǎng)格法向分塊位置類似，即分塊位置1遠離邊界層，分塊位置2在邊界層內(nèi)。圖12中在分塊位置2分塊得到的多塊網(wǎng)格的流場在計算時CFL數(shù)取不到初始網(wǎng)格和分塊位置1計算時的3500，最大可取到2000，說明三維計算時在法向進行網(wǎng)格分塊導(dǎo)致格式穩(wěn)定性下降。從是否采用隱式對接邊界條件的收斂曲線對比可以看出，采用隱式對接邊界條件能夠改善DPLR算法在法向分塊后收斂速率下降的問題。

圖12 三維球頭繞流網(wǎng)格在法向分塊時的殘差收斂曲線(CFL=2000)Fig.12 Residual history of flow around 3D sphere when splitting grid blocks in normal direction(CFL=2000)

3 返回艙繞流計算分析

返回艙外形復(fù)雜，包含凸起物、耳片等局部結(jié)構(gòu)。計算網(wǎng)格分區(qū)如圖13所示，共77個分區(qū)，約116萬網(wǎng)格數(shù)。在這樣的網(wǎng)格上計算，DPLR算法將會出現(xiàn)相鄰網(wǎng)格求解方向交叉的現(xiàn)象，這是與上文球頭繞流計算的不同之處，也是對本文提出的隱式對接邊界條件模型的重要檢驗。為了與風(fēng)洞試驗結(jié)果對比，本文選擇風(fēng)洞試驗狀態(tài)[19]作為計算狀態(tài)：模型縮比6%，馬赫數(shù)9.96，基于模型長度的雷諾數(shù)2.04×105，總溫1005 K，總壓1.976 MPa，迎角10°～30°。風(fēng)洞試驗在FD-17AΦ1 m高超聲速低密度風(fēng)洞[20]進行，這是一座是高壓下吹-真空抽吸暫沖式風(fēng)洞。分析試驗條件，自由來流努森數(shù)為7.08×10-5，流動處于連續(xù)流范圍；雷諾數(shù)較低，流動保持層流的可能性大；總溫較低，流動的真實氣體效應(yīng)很弱。因而在計算中，本文選擇層流和完全氣體模型。

圖13 返回艙物面和對稱面網(wǎng)格拓撲結(jié)構(gòu)Fig.13 Grid topology at wall and symmetry boundary of capsule module

圖14給出迎角22°狀態(tài)的殘差收斂曲線，圖中包含了DPLR和LU-SGS兩種方法的對比，其中DPLR算法的CFL數(shù)為5000，而LU-SGS算法的CFL數(shù)取到了1.0×106，圖中的橫坐標是CPU時間。兩種隱式算法都基于原始網(wǎng)格分區(qū)在Intel(R) Xeon(R) CPU E5-2680共8核上進行并行計算，未考慮負載均衡。由圖可見，LU-SGS算法殘差最多只能下降3個量級，而DPLR算法的殘差可以下降到機器零。這不僅說明了DPLR算法和本文提出的隱式對接邊界條件能夠成功地應(yīng)用于復(fù)雜外形，也再次說明了DPLR相對于LU-SGS在收斂速度上的優(yōu)勢。

圖15給出返回艙在迎角22°對稱面和物面的壓力、馬赫數(shù)等值線云圖以及對稱面的流線結(jié)構(gòu)，可以看到返回艙頭部的弓形激波和表面凸起物附近模擬得到的流場結(jié)構(gòu)較為清晰。從圖15(b)的對稱面流線結(jié)構(gòu)可以看到，在返回艙背部的凸起物前部，由于凸起物的遮擋，流動被局部滯止，馬赫數(shù)減小并發(fā)生了流動分離；而在凸起物的后部，流動分離形成旋渦，部分分離的流線折向返回艙底部并與下部流線匯集。圖16給出返回艙的阻力、升力和俯仰力矩系數(shù)隨迎角變化關(guān)系?？煽闯觯枇蜕ο禂?shù)隨迎角的增大而減?。桓┭隽叵禂?shù)隨迎角的增大而減小，說明返回艙具有縱向靜穩(wěn)定性。計算和試驗符合較好，阻力系數(shù)相對偏差不超過為3%，升力系數(shù)相對偏差不超過2%；計算與試驗預(yù)測的配平迎角分別為22.05°和21.63°，相對偏差為1.9%。圖15、圖16較好地揭示了返回艙繞流面貌。

圖14 返回艙繞流DPLR與LU-SGS收斂曲線的比較Fig.14 Comparison between DPLR and LU-SGS residual curve for flowfield computation of reentry capsule

(a) 壓力云圖

(b) 馬赫數(shù)云圖

(a) 阻力和升力系數(shù)

(b) 俯仰力矩系數(shù)

4 結(jié)束語

由于DPLR算法需要和隱式邊界條件相結(jié)合，為把DPLR應(yīng)用到返回艙再入近空間連續(xù)流區(qū)氣動力/熱復(fù)雜繞流多塊對接結(jié)構(gòu)網(wǎng)格的計算，本文首先提出一種隱式對接邊界條件的處理方法，通過在對接邊界處應(yīng)用隱式邊界條件把邊界條件對隱式求解的影響考慮進來。然后研究了網(wǎng)格分塊方式對算法穩(wěn)定性和收斂性能的影響。最后進行了返回艙繞流的數(shù)值模擬，給出氣動力特性，并與風(fēng)洞試驗結(jié)果對比，考察算法在復(fù)雜外形數(shù)值模擬中的性能。可得出如下結(jié)論：

(1) DPLR算法采用松弛方法把求解方向(一般是物面法向)和其它方向解耦，在非求解方向上進行網(wǎng)格分塊不會對算法的收斂性能產(chǎn)生影響。而在求解方向上對網(wǎng)格分塊，特別是分塊位置在邊界層內(nèi)，會降低算法的穩(wěn)定性和收斂性能。

(2) 本文提出的隱式對接邊界條件處理方法，能夠改善在求解方向進行網(wǎng)格分塊時造成的算法穩(wěn)定性和收斂性能下降的問題。

(3) 如果應(yīng)用多塊結(jié)構(gòu)網(wǎng)格和DPLR算法進行數(shù)值模擬，在用網(wǎng)格軟件生成網(wǎng)格或基于并行計算的目的剖分網(wǎng)格時，本文的研究給出一些注意事項：可以在流向和周向進行網(wǎng)格分區(qū)，盡量不要在物面法向進行網(wǎng)格分區(qū)，更不要在邊界層內(nèi)沿物面法向進行網(wǎng)格分區(qū)。這為相關(guān)方面的研究提供指導(dǎo)原則。

(4) 返回艙連續(xù)流區(qū)的復(fù)雜繞流計算結(jié)果表明，DPLR算法結(jié)合隱式對接邊界條件能夠成功地應(yīng)用于復(fù)雜外形結(jié)構(gòu)網(wǎng)格數(shù)值模擬，并以較快的速度收斂，大大縮短了計算所需時間。

本文的研究是初步的階段性工作下，一步將把DPLR算法應(yīng)用到湍流以及熱化學(xué)非平衡流動的數(shù)值模擬。