張培紅， 張耀冰， 趙 煒， 周桂宇， 吳曉軍， 楊福軍

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000)

1 引 言

非結(jié)構(gòu)網(wǎng)格方法消除了網(wǎng)格節(jié)點(diǎn)的結(jié)構(gòu)性限制，易于控制網(wǎng)格單元的大小、形狀及網(wǎng)格點(diǎn)的位置，比結(jié)構(gòu)網(wǎng)格具有更大的靈活性，對(duì)復(fù)雜外形的適應(yīng)能力非常強(qiáng)，是最近十多年發(fā)展最快的一種CFD技術(shù)[1,2]。以AIAA阻力會(huì)議為例，第一屆會(huì)議18個(gè)程序中有7個(gè)為非結(jié)構(gòu)程序，第二屆會(huì)議30個(gè)結(jié)果中非結(jié)構(gòu)程序有11個(gè)，第三屆會(huì)議26個(gè)參加者中有11個(gè)使用非結(jié)構(gòu)程序，第四屆會(huì)議提交的28個(gè)結(jié)果中非結(jié)構(gòu)網(wǎng)格的計(jì)算結(jié)果達(dá)到了17個(gè)[3-7]。美國(guó)NASA的USM3Dns，NSU3D和FUN3D、美國(guó)空軍的Cobalt60、波音公司的BCFD、歐洲D(zhuǎn)LR的TAU、FOI的EDGE以及日本JAXA的TAS等著名CFD軟件都是基于非結(jié)構(gòu)網(wǎng)格[8-10]；商業(yè)的非結(jié)構(gòu)網(wǎng)格軟件也非常多，如Metacomp Technologies公司的CFD++、ANSYS公司的CFX和Fluent以及CFDRC公司的CFD -FASTRAN等；還有許多開(kāi)源軟件也是非結(jié)構(gòu)網(wǎng)格代碼，如OpenFlower和OpenFoam等[11]。

但是，由于非結(jié)構(gòu)網(wǎng)格在黏性區(qū)域很難使用類似結(jié)構(gòu)網(wǎng)格的大長(zhǎng)寬比網(wǎng)格，導(dǎo)致黏性預(yù)測(cè)精度降低，嚴(yán)重阻礙非結(jié)構(gòu)網(wǎng)格的應(yīng)用和發(fā)展[12-14]。近年來(lái)，CFD工作者針對(duì)非結(jié)構(gòu)網(wǎng)格的特點(diǎn)開(kāi)展了大量的網(wǎng)格生成、梯度求解和計(jì)算格式等研究，特別是非結(jié)構(gòu)混合網(wǎng)格生成和求解技術(shù)的發(fā)展，解決了非結(jié)構(gòu)網(wǎng)格氣動(dòng)特性預(yù)測(cè)精度低的缺陷，促進(jìn)了非結(jié)構(gòu)網(wǎng)格在飛行器氣動(dòng)特性預(yù)測(cè)模擬中的應(yīng)用[15-20]。非結(jié)構(gòu)混合網(wǎng)格一般由四面體、六面體、三棱柱或者金字塔組成，本文在生成混合網(wǎng)格時(shí)，在物面粘性作用區(qū)生成大伸展比三棱柱形和金字塔形網(wǎng)格；在其他流動(dòng)區(qū)域采用陣面推進(jìn)方法生成四面體網(wǎng)格。但對(duì)于高超聲速氣動(dòng)熱流計(jì)算，由于高超聲速流動(dòng)中含有激波等強(qiáng)間斷，流體黏性起主導(dǎo)作用，物理黏性摩擦決定了熱流的大小，對(duì)黏性計(jì)算精度提出了更高的要求，在計(jì)算中受到物理模型、數(shù)值格式、網(wǎng)格分布、收斂過(guò)程以及熱流后處理方法等諸多因素的影響和制約[21-24]。即使是CFD方法能夠給出較為精確的氣動(dòng)力結(jié)果，對(duì)于壁面熱流的預(yù)測(cè)往往也可能會(huì)有數(shù)量級(jí)的差異。高超聲速氣動(dòng)熱的數(shù)值模擬一直是CFD研究的難點(diǎn)和熱點(diǎn)問(wèn)題，尤其是采用非結(jié)構(gòu)網(wǎng)格計(jì)算熱流，還存在較大的困難，很難達(dá)到和結(jié)構(gòu)網(wǎng)格相當(dāng)?shù)挠?jì)算精度。主要有兩個(gè)原因，一是非結(jié)構(gòu)四面體網(wǎng)格單元無(wú)法保證對(duì)稱性；二是四面體網(wǎng)格至少有一個(gè)面是和激波斜交，無(wú)法保證網(wǎng)格與激波正交，影響激波捕捉精度。

針對(duì)非結(jié)構(gòu)四面體網(wǎng)格在計(jì)算熱流時(shí)的缺陷，Gnoffo[25]發(fā)展了多維通量重構(gòu)方法，通過(guò)改進(jìn)通量重構(gòu)方法，減少四面體網(wǎng)格的影響，在六面體剖分的純四面體網(wǎng)格上得到了較好的熱流計(jì)算結(jié)果。McCloud[26]在網(wǎng)格自適應(yīng)技術(shù)方面開(kāi)展了相關(guān)研究，發(fā)展了魯棒穩(wěn)定的激波區(qū)域網(wǎng)格自適應(yīng)技術(shù)，基于常規(guī)的混合網(wǎng)格在探測(cè)到的激波面附近重新生成三棱柱網(wǎng)格，其余區(qū)域使用四面體網(wǎng)格填充，較好地改善了非結(jié)構(gòu)/混合網(wǎng)格上的熱流計(jì)算結(jié)果。

近年來(lái)，國(guó)內(nèi)在結(jié)構(gòu)網(wǎng)格的熱流計(jì)算方面開(kāi)展了較多的研究工作，在非結(jié)構(gòu)/混合網(wǎng)格的熱流計(jì)算方法研究方面也取得了一定進(jìn)展，但開(kāi)展的工作相對(duì)較少。賀立新等[27]使用有限元/有限體積混合方法基于混合網(wǎng)格在雙橢球上得到了較好的熱流計(jì)算結(jié)果。李澤禹[28]基于非平衡氣體在混合網(wǎng)格上就球頭繞流算例進(jìn)行了數(shù)值計(jì)算，發(fā)現(xiàn)很難得到光滑連續(xù)的分布結(jié)果。許和勇等[29]利用自適應(yīng)網(wǎng)格技術(shù)開(kāi)展了典型高超聲速飛行器氣動(dòng)熱環(huán)境的數(shù)值模擬。萬(wàn)云博[30]應(yīng)用多維方法基于常規(guī)混合網(wǎng)格對(duì)包含鈍錐、鈍雙錐和雙橢球等多種外形進(jìn)行了數(shù)值模擬，得到了與實(shí)驗(yàn)值吻合較好的熱流計(jì)算結(jié)果。

相關(guān)研究表明[23,24]，網(wǎng)格分布是影響熱流計(jì)算的重要因素之一。本文基于非結(jié)構(gòu)混合網(wǎng)格，以典型鈍錐標(biāo)模外形的高超聲速繞流為研究對(duì)象，驗(yàn)證了計(jì)算方法的計(jì)算精度和可靠性，開(kāi)展了不同網(wǎng)格形式和第一層網(wǎng)格不同間距的影響研究，為采用非結(jié)構(gòu)混合網(wǎng)格氣動(dòng)熱計(jì)算時(shí)的網(wǎng)格生成形式和第一層間距等提供了參考建議和指導(dǎo)。

2 數(shù)值方法

本文計(jì)算采用國(guó)家數(shù)值風(fēng)洞工程(NNW)項(xiàng)目團(tuán)隊(duì)自主開(kāi)發(fā)的基于格心的非結(jié)構(gòu)混合網(wǎng)格流場(chǎng)解算器FlowStar。該解算器支持任意形狀的網(wǎng)格單元，針對(duì)非結(jié)構(gòu)混合網(wǎng)格的特點(diǎn)，通過(guò)改進(jìn)傳統(tǒng)的Green-Gauss梯度求解方法和傳統(tǒng)的Roe格式Harten-Yee熵修正方法，提出了一種新的節(jié)點(diǎn)型Green-Gauss梯度求解方法和Harten-Yee熵修正改進(jìn)方法，大大提高了非結(jié)構(gòu)混合網(wǎng)格粘性計(jì)算精度，軟件經(jīng)過(guò)了大量標(biāo)準(zhǔn)算例的考核驗(yàn)證，在航空航天領(lǐng)域得到了廣泛的應(yīng)用。

采用有限體積法對(duì)空間進(jìn)行離散，未知變量位于網(wǎng)格單元的體心?？刂品匠虨槭睾阈问降姆嵌ǔ？蓧嚎sN-S方程,

(1)

式中Ω為控制體的體積，?Ω為控制體封閉面的面積，W為守恒變量，F(xiàn)c為無(wú)粘通量，F(xiàn)v為粘性通量。

本文主控方程對(duì)流項(xiàng)采用二階迎風(fēng)Roe通量差分裂格式進(jìn)行離散，粘性項(xiàng)采用中心差分格式離散，并采用多重網(wǎng)格技術(shù)進(jìn)行收斂加速。湍流模型采用SST兩方程湍流模型，湍流控制方程空間離散采用一階迎風(fēng)格式。主控方程和湍流方程的時(shí)間迭代采用LU-SGS(Lower Upper-Symmetric Gauss -Seidel)方法 。

Roe通量差分分裂格式具有很高的分辨定態(tài)激波的能力，較其他的一些通量分裂格式(如Van Leer格式)，該格式的耗散較小。Roe格式在控制體單元面上的通量表達(dá)式為

(2)

式中Fc為無(wú)粘通量，W為守恒變量，ARoe為Roe平均矩陣。

Roe平均矩陣和左右狀態(tài)差的乘積的求法為

|ARoe|I J(WJ-WI)=|ΔF1|+|ΔF2,3,4|+|ΔF5|

(3)

為提高流場(chǎng)計(jì)算效率，采用了基于MPI的大規(guī)模并行計(jì)算技術(shù)、多重網(wǎng)格方法、低速預(yù)處理技術(shù)和當(dāng)?shù)貢r(shí)間步方法等加速收斂技術(shù)。

3 計(jì)算模型及網(wǎng)格

鈍錐標(biāo)模外形高超聲速繞流是典型的高超聲速三維粘性流動(dòng)。圖1為鈍錐標(biāo)模的幾何尺寸示意圖，圖2為鈍錐標(biāo)模三維計(jì)算模型。鈍錐球頭半徑為R=27.94 mm，半錐角為15°，長(zhǎng)度L=447.04 mm。計(jì)算來(lái)流條件為Ma=10.6，Re=3.973×106/m，T∞=47.3 K，Tw=294.44 K，分別計(jì)算了0°和20°兩個(gè)工況。

圖1 鈍錐幾何尺寸

圖2 鈍錐計(jì)算模型

本文采用的非結(jié)構(gòu)混合網(wǎng)格，均是以生成的一套結(jié)構(gòu)網(wǎng)格grid-str為基礎(chǔ)，通過(guò)剖分得到。結(jié)構(gòu)網(wǎng)格grid-str第一層間距5×10-6m，網(wǎng)格雷諾數(shù)約為20，網(wǎng)格量約為196.6萬(wàn)個(gè)單元，其中周向點(diǎn)數(shù)257個(gè)，法向點(diǎn)數(shù)81個(gè)。圖3給出了結(jié)構(gòu)網(wǎng)格grid-str鈍錐錐段和頭部表面網(wǎng)格。

圖3 鈍錐錐段和頭部表面網(wǎng)格(grid-str)

4 網(wǎng)格生成形式影響

為研究網(wǎng)格生成形式的影響，采用roe格式，對(duì)不同形式網(wǎng)格進(jìn)行了計(jì)算研究。在結(jié)構(gòu)網(wǎng)格grid-str的基礎(chǔ)上，通過(guò)對(duì)頭部網(wǎng)格采用不同剖分方式，共生成了7套非結(jié)構(gòu)混合網(wǎng)格，分別為grid-unstr1，grid-unstr2，grid-unstr3，grid-unstr4，grid-unstr5，grid-unstr6和grid-unstr7，保持第一層網(wǎng)格間距5×10-6m不變，即保持網(wǎng)格雷諾數(shù)20不變。其中，grid-unstr1，grid-unstr2，grid-unstr3，grid-unstr4，grid-unstr5和grid-unstr6的頭部網(wǎng)格在grid-str基礎(chǔ)上采用不同方式剖分，尾部網(wǎng)格采用Advancing Front剖分直接生成非結(jié)構(gòu)網(wǎng)格；grid-unstr7的頭部網(wǎng)格與grid-unstr6完全相同，尾部網(wǎng)格是在grid-str基礎(chǔ)上每個(gè)網(wǎng)格剖分成4個(gè)三角形單元。

grid-unstr1，grid-unstr2和grid-unstr3三套網(wǎng)格的頭部網(wǎng)格是在grid-str基礎(chǔ)上把每個(gè)四邊形網(wǎng)格單元剖分為2個(gè)三角形單元，尾部采用 Advancing Front剖分直接生成非結(jié)構(gòu)網(wǎng)格，不同在于三套網(wǎng)格頭部采用的剖分方式，grid-unstr1采用繞圈剖分，grid-unstr2采用交叉剖分，grid-unstr3 采用放射剖分。grid-unstr1網(wǎng)格量為4125951個(gè)單元，grid-unstr2和grid-unstr3的網(wǎng)格量與grid-unstr1相近。圖4給出了三種剖分方式得到的頭部中心附近網(wǎng)格比較。

grid-unstr4采用Delaunay剖分方式，網(wǎng)格分布具有明顯的任意性，grid-unstr5采用Advancing Front剖分方式，除頭部中心附近外，網(wǎng)格具有較好的排列。grid-unstr4和grid-unstr5的尾部也均采用Advancing Front剖分直接生成非結(jié)構(gòu)網(wǎng)格。

圖4 不同剖分方式頭部中心附近網(wǎng)格比較

grid-unstr6與grid-unstr2一樣，頭部網(wǎng)格是在grid-str網(wǎng)格的基礎(chǔ)上采用交叉剖分方式生成，不同在于每個(gè)網(wǎng)格剖分成4個(gè)三角形單元，相比grid-unstr2網(wǎng)格量增加一倍。同樣，grid-unstr6尾部也采用Advancing Front剖分直接生成非結(jié)構(gòu)網(wǎng)格。grid-unstr7的頭部網(wǎng)格與grid-unstr6完全相同，不同在于尾部網(wǎng)格也是在結(jié)構(gòu)網(wǎng)格grid-str的基礎(chǔ)上把每個(gè)網(wǎng)格剖分成4個(gè)三角形單元生成。圖5給出了grid-unstr1和grid-unstr7錐段表面網(wǎng)格分布的比較。圖6給出了grid-unstr1，grid-unstr2，grid-unstr3，grid-unstr4，grid-unstr5，grid-unstr6和grid-unstr7不同剖分方式頭部網(wǎng)格比較。

圖5 grid-unstr1和grid-unstr7錐段表面網(wǎng)格分布比較

圖6 頭部網(wǎng)格分布比較

圖7為Ma=10.6，α=0°時(shí)，不同網(wǎng)格頭部壁面的熱流分布?？梢钥闯觯煌W(wǎng)格頭部駐點(diǎn)附近壁面熱流分布有明顯的差異，相對(duì)來(lái)說(shuō)，grid-unstr1，grid-unstr2，grid-unstr6和grid-unstr7四套網(wǎng)格計(jì)算得到的壁面熱流周向分布最好，但grid-unstr1 和grid-unstr2兩套網(wǎng)格計(jì)算得到的駐點(diǎn)處熱流值不是最大，明顯不合理，grid-unstr6和grid-unstr7兩套網(wǎng)格計(jì)算得到的壁面熱流分布最好，兩者基本上完全相同。

圖8給出了Ma=10.6，α=20° 時(shí)，不同網(wǎng)格頭部壁面熱流不同剖面分布與實(shí)驗(yàn)比較。計(jì)算結(jié)果與試驗(yàn)值吻合較好，除六面體網(wǎng)格grid-str外，其他不同網(wǎng)格的計(jì)算結(jié)果相差很小。

圖8 不同網(wǎng)格壁面熱流與試驗(yàn)比較(Ma=10.6,α=20°)

綜合不同網(wǎng)格計(jì)算得到的壁面熱流云圖分布和不同剖面壁面熱流實(shí)驗(yàn)對(duì)比結(jié)果分析，要想較好地計(jì)算鈍錐標(biāo)模頭部的熱流分布，在研究的幾種網(wǎng)格形式中，頭部物面網(wǎng)格最好采用四邊形或四邊形交叉剖分得到的三角形網(wǎng)格，根據(jù)商業(yè)軟件生成網(wǎng)格的便利性，推薦采用四邊形交叉剖分得到的三角形網(wǎng)格。同時(shí)頭部激波附近網(wǎng)格最好沿激波排列，激波靠近附面層的情況可以采用附面層三棱柱網(wǎng)格適當(dāng)多向外推幾層，激波遠(yuǎn)離附面層時(shí)可以采用baffle面的方式生成。

5 網(wǎng)格加密及第一層間距影響

為研究網(wǎng)格加密對(duì)計(jì)算結(jié)果的影響，在計(jì)算結(jié)果較好的網(wǎng)格grid-unstr6的基礎(chǔ)上，將物面和空間網(wǎng)格尺度都加密一倍，得到grid-unstr8網(wǎng)格。圖9給出了α=20° 時(shí)，網(wǎng)格加密前(grid-unstr6)和網(wǎng)格加密后(grid-unstr8)計(jì)算得到的壁面熱流分布云圖對(duì)比。圖10給出了α=20°時(shí)，網(wǎng)格加密前(grid-unstr6)和網(wǎng)格加密后(grid-unstr8)計(jì)算得到的壁面熱流分布比較。可以看出，與grid-unstr6相比，加密后網(wǎng)格grid-unstr8計(jì)算的壁面熱流分布并沒(méi)有明顯的改進(jìn)，說(shuō)明針對(duì)目前的計(jì)算狀態(tài)，已經(jīng)無(wú)法通過(guò)網(wǎng)格加密明顯改進(jìn)計(jì)算結(jié)果。

圖9 網(wǎng)格加密前后壁面熱流分布云圖對(duì)比(Ma=10.6,α=20°)

圖10 網(wǎng)格加密前后不同位置壁面熱流分布與試驗(yàn)比較(Ma=10.6,α=20°)

為進(jìn)一步研究第一層網(wǎng)格間距的影響，在網(wǎng)格grid-unstr6的基礎(chǔ)上，保持表面網(wǎng)格等其他參數(shù)不變，將法向第一層減小為五分之一，即第一層間距變?yōu)?×10-6m得到grid-unstr9，對(duì)應(yīng)的網(wǎng)格雷諾數(shù)約為4；將法向第一層增大為五倍，即第一層間距變?yōu)?.5×10-5m得到grid-unstr10，對(duì)應(yīng)的網(wǎng)格雷諾數(shù)約為100。

圖11給出了α=20°時(shí)，不同第一層網(wǎng)格間距得到的壁面熱流分布云圖。圖12給出了α=20° 時(shí)，不同第一層網(wǎng)格間距得到的壁面熱流分布比較。三套網(wǎng)格得到的計(jì)算結(jié)果沒(méi)有明顯差別，相對(duì)來(lái)說(shuō)，法向最密的grid-unstr9網(wǎng)格得到的駐點(diǎn)壁面熱流分布軸對(duì)稱性要好于其他兩個(gè)，而第一層網(wǎng)格間距最大的grid-unstr10網(wǎng)格計(jì)算得到的壁面熱流分布與物面網(wǎng)格明顯相關(guān)。但是從計(jì)算收斂的角度，grid-unstr9由于網(wǎng)格太密，計(jì)算過(guò)程中會(huì)產(chǎn)生明顯的網(wǎng)格剛性現(xiàn)象，殘差下降異常慢，計(jì)算效率急劇下降。因此，法向間距不宜取得過(guò)密。

圖11 網(wǎng)格第一層間距對(duì)頭部壁面熱流計(jì)算結(jié)果的影響(Ma=10.6,α=20°)

圖12 網(wǎng)格第一層間距對(duì)不同位置壁面熱流計(jì)算結(jié)果影響(Ma=10.6,α=20°)

6 結(jié) 論

本文以典型鈍錐標(biāo)模外形的高超聲速繞流為研究對(duì)象，采用基于非結(jié)構(gòu)混合網(wǎng)格建立的熱流計(jì)算方法，開(kāi)展了不同網(wǎng)格形式和第一層網(wǎng)格不同間距的影響研究，可以得到以下結(jié)論。

(1) 建立的熱流計(jì)算方法，可以較好地模擬鈍錐標(biāo)模外形的高超聲速繞流壁面熱流分布。

(2) 熱流計(jì)算時(shí)，頭部物面網(wǎng)格最好采用四邊形或四邊形交叉剖分得到的三角形網(wǎng)格。

(3) 物面法向的網(wǎng)格雷諾數(shù)取20左右，過(guò)密的網(wǎng)格可能會(huì)產(chǎn)生網(wǎng)格剛性，降低計(jì)算效率。