DPLR隱格式在多塊結(jié)構(gòu)網(wǎng)格的計(jì)算實(shí)現(xiàn)

黨雷寧, 白智勇, 柳 森

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 超高速空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)

0 引 言

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

在CFD并行計(jì)算中，普遍采用基于網(wǎng)格分區(qū)的并行策略。在CFD時(shí)間推進(jìn)中，為了不受穩(wěn)定性條件的限制以及出于加速流動(dòng)模擬收斂的目的，一般采用各種隱式計(jì)算方法。然而，傳統(tǒng)的以LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法[3]為代表的隱式近似因子分解(AF)類方法應(yīng)用于并行計(jì)算時(shí)，存在兩個(gè)問(wèn)題。一方面，隱式算法中相鄰網(wǎng)格單元之間較強(qiá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é)非平衡流動(dòng)的數(shù)值模擬中，它的特點(diǎn)是收斂很快。GSLR方法的思想是：在黏性流動(dòng)中，相較于流向和周向，流場(chǎng)梯度在物面法向更大，流動(dòng)在物面法向耦合更緊密，數(shù)值求解時(shí)在物面法向把隱式部分的對(duì)角項(xiàng)和非對(duì)角項(xiàng)耦合求解才能反映流動(dòng)的物理本質(zhì)。GSLR方法在流向和周向采用Gauss-Seidel掃描，不利于并行計(jì)算。NASA Ames研究中心的Wright(1997)[8]，以GSLR為出發(fā)點(diǎn)，把Gauss-Seidel掃描改為松弛過(guò)程，消除相鄰物面法向網(wǎng)格線流場(chǎng)變量之間的數(shù)據(jù)依賴，既保留了GSLR收斂快的特點(diǎn)，又使它有利于并行計(jì)算，這就是DPLR(Data-Parallel Line Relaxation)方法。

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

1 數(shù)值計(jì)算方法

1.1 控制方程

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

式中，Q為守恒變量矢量，E、F為無(wú)黏通量矢量，Ev、Fv為黏性通量矢量，它們的表達(dá)式詳見(jiàn)文獻(xiàn)[11]。

1.2 空間項(xiàng)離散

應(yīng)用有限體積法，以網(wǎng)格單元為控制體，對(duì)控制方程進(jìn)行離散。在有限體積法中，流場(chǎng)變量?jī)?chǔ)存在單元格心，在單元內(nèi)部均勻分布。本文記流向?yàn)棣危聵?biāo)索引為i；法向?yàn)棣牵聵?biāo)索引為j。把控制方程(1)在控制體i,j積分，并應(yīng)用Gauss-Green定理，則有：

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

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

1.3 時(shí)間項(xiàng)離散

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

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

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

其中:

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

1.4 DPLR格式

設(shè)η方向?yàn)槲锩娣ㄏ?，則DPLR格式[8]可寫為:

對(duì)于m=1,mmax

1.5 邊界條件

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

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

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

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

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

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

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

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

dQ0=HdQ1(11)

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

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

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

表1 算例驗(yàn)證的計(jì)算條件Table 1 Computational condition of two cases

2 算例驗(yàn)證及討論

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

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

(a) 二維 (b) 三維

2.2 計(jì)算結(jié)果正確性驗(yàn)證

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

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

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

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

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

圖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)格進(jìn)行分塊。有兩種分塊方式，一是在流向和周向分塊；二是在法向分塊。

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

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

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

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

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

3 返回艙繞流計(jì)算分析

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

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

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

圖15給出返回艙在迎角22°對(duì)稱面和物面的壓力、馬赫數(shù)等值線云圖以及對(duì)稱面的流線結(jié)構(gòu)，可以看到返回艙頭部的弓形激波和表面凸起物附近模擬得到的流場(chǎng)結(jié)構(gòu)較為清晰。從圖15(b)的對(duì)稱面流線結(jié)構(gòu)可以看到，在返回艙背部的凸起物前部，由于凸起物的遮擋，流動(dòng)被局部滯止，馬赫數(shù)減小并發(fā)生了流動(dòng)分離；而在凸起物的后部，流動(dòng)分離形成旋渦，部分分離的流線折向返回艙底部并與下部流線匯集。圖16給出返回艙的阻力、升力和俯仰力矩系數(shù)隨迎角變化關(guān)系。可看出，阻力和升力系數(shù)隨迎角的增大而減??；俯仰力矩系數(shù)隨迎角的增大而減小，說(shuō)明返回艙具有縱向靜穩(wěn)定性。計(jì)算和試驗(yàn)符合較好，阻力系數(shù)相對(duì)偏差不超過(guò)為3%，升力系數(shù)相對(duì)偏差不超過(guò)2%；計(jì)算與試驗(yàn)預(yù)測(cè)的配平迎角分別為22.05°和21.63°，相對(duì)偏差為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é)束語(yǔ)

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

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

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

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

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

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