基于結(jié)構(gòu)網(wǎng)格的柵格翼繞流數(shù)值模擬

周培培

（中國(guó)航天空氣動(dòng)力技術(shù)研究院，北京 100074）

0 引 言

柵格翼是由許多柵格壁連接起來(lái)的復(fù)雜的空間受力系統(tǒng)，常用的有框式柵格翼和蜂窩式柵格翼。相對(duì)于傳統(tǒng)舵面而言，柵格翼質(zhì)量輕，適用攻角范圍大、鉸鏈力矩小，因而舵機(jī)功率很小，允許采用較輕和較小的控制系統(tǒng)，導(dǎo)彈質(zhì)量明顯減輕，具有緊貼彈體折疊安裝減少空間等優(yōu)點(diǎn)。因此在各種飛行器的設(shè)計(jì)中得到廣泛應(yīng)用，如俄羅斯的R－77、AAM－AE導(dǎo)彈控制舵面，美國(guó)的MOAB炸彈，中國(guó)飛船逃逸系統(tǒng)的穩(wěn)定翼面、德國(guó)高超聲速導(dǎo)彈的控制舵面等。

鑒于柵格翼在亞、跨、超聲速流中的獨(dú)特氣動(dòng)特性，國(guó)內(nèi)外開展了很多研究工作。在數(shù)值模擬方面，國(guó)外的 Despirito，et al［1－3］采用非結(jié)構(gòu)網(wǎng)格對(duì)柵格翼導(dǎo)彈開展了粘性繞流的數(shù)值模擬研究，并將計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了詳細(xì)的對(duì)比分析，但得到的氣動(dòng)力系數(shù)存在一定的偏差；Chen和Khalid［4］等人基于歐拉方程，針對(duì)柵格翼的柵格壁厚度、前緣形狀等開展了詳細(xì)的研究工作；Lin和Huang［5］等人采用有效的分區(qū)方法，針對(duì)柵格翼導(dǎo)彈生成了分區(qū)結(jié)構(gòu)網(wǎng)格，結(jié)合Baldwin－Barth湍流模型求解N－S方程，對(duì)柵格翼導(dǎo)彈進(jìn)行數(shù)值模擬，并且與試驗(yàn)進(jìn)行了對(duì)比分析，吻合程度較高；Ravubdra［6］等人對(duì)二維的柵格框生成結(jié)構(gòu)網(wǎng)格，對(duì)單獨(dú)的柵格翼生成非結(jié)構(gòu)網(wǎng)格，以此對(duì)柵格框內(nèi)流動(dòng)機(jī)理及氣動(dòng)特性展開了研究；Yan［7］采用結(jié)構(gòu)非結(jié)構(gòu)混合網(wǎng)格的方法研究了后掠式柵格翼的減租特性。國(guó)內(nèi)的馬明生、吳曉軍、劉剛［9－12］等人采用結(jié)構(gòu)／非結(jié)構(gòu)混合網(wǎng)格對(duì)單個(gè)柵格框和柵格翼翼身組合體的氣動(dòng)特性進(jìn)行了數(shù)值模擬研究，在柵格翼區(qū)域采用非結(jié)構(gòu)網(wǎng)格，其它區(qū)域采用結(jié)構(gòu)網(wǎng)格；Cai［13］研究了柵格翼的減阻特性；鄧帆、陳少松［14－15］等人的研究主要集中于單個(gè)柵格框的氣動(dòng)性能上；王斌等［16］采用直角網(wǎng)格，王學(xué)占等［17］采用結(jié)構(gòu)重疊網(wǎng)格、周張等［18］采用面元法，分別對(duì)柵格翼導(dǎo)彈的流場(chǎng)進(jìn)行了數(shù)值模擬研究。

對(duì)柵格翼進(jìn)行數(shù)值模擬的困難在于生成高質(zhì)量的計(jì)算網(wǎng)格，以前的研究中大多采用的是直角網(wǎng)格、非結(jié)構(gòu)網(wǎng)格、混合網(wǎng)格等。本文的研究?jī)?nèi)容在于采取一種有效的分塊結(jié)構(gòu)網(wǎng)格生成方法，針對(duì)柵格翼這種復(fù)雜框式結(jié)構(gòu)生成高質(zhì)量的貼體結(jié)構(gòu)網(wǎng)格。以此為基礎(chǔ)開展柵格翼導(dǎo)彈亞跨超聲速流場(chǎng)的數(shù)值模擬，為柵格翼導(dǎo)彈繞流精確數(shù)值預(yù)測(cè)提供行之有效的方法。

1 數(shù)值模擬方法

1.1 控制方程

對(duì)于三維非定?？蓧嚎s粘性流動(dòng)，一般曲線坐標(biāo)系下的N－S方程可寫作：

無(wú)粘通量矢量為：

粘性通量矢量為：

對(duì)于湍流的模擬采用兩方程的k－ε渦粘性湍流模型。

1.2 數(shù)值方法及邊界條件

對(duì)于流場(chǎng)的求解，采用有限體積法進(jìn)行空間離散，隱式方法結(jié)合多重網(wǎng)格加速收斂技術(shù)進(jìn)行推進(jìn)求解?？臻g無(wú)粘通量采用ROE的通量差分分裂方法來(lái)離散，粘性通量采用二階中心差分格式離散，并且采用MIN MOD限制器以保證格式單調(diào)性。時(shí)間方向采用隱式LU－SGS方法推進(jìn)求解，并結(jié)合多重網(wǎng)格技術(shù)加速收斂。超聲速遠(yuǎn)場(chǎng)邊界條件入流邊界所有參數(shù)為自由來(lái)流值，出流邊界所有參數(shù)由內(nèi)流場(chǎng)外插得到；亞聲速遠(yuǎn)場(chǎng)邊界條件采用Riemann不變量處理的無(wú)反射邊界條件；柵格翼導(dǎo)彈對(duì)稱面按對(duì)稱條件處理，物面采用無(wú)滑移邊界條件。

1.3 多重網(wǎng)格技術(shù)

多重網(wǎng)格技術(shù)是數(shù)值算法中一種加速收斂的迭代技術(shù)。其基本思想就是細(xì)網(wǎng)格松弛、粗網(wǎng)格修正和套迭代技術(shù)。即先在細(xì)網(wǎng)格上迭代以消除高頻振蕩誤差，再在粗網(wǎng)格上迭代消除低頻光滑誤差，如此逐步使網(wǎng)格粗化從而把各種頻率誤差消除，再有套迭代技術(shù)將粗網(wǎng)格上的解依次返回到各級(jí)細(xì)網(wǎng)格上。計(jì)算過程中采用多重網(wǎng)格法中非線性方程的全近似格式（FAS格式）；網(wǎng)格粗化則采用標(biāo)準(zhǔn)粗化策略；套迭代技術(shù)采用W循環(huán)，保持收斂因子不隨網(wǎng)格層的變化而變化，魯棒性好；插值算子采用三線性插值算子；限制算子則為體積加權(quán)算子。

1.4 數(shù)值計(jì)算模型

數(shù)值計(jì)算模型采用在加拿大DREV風(fēng)洞中的進(jìn)行大量試驗(yàn)的柵格翼導(dǎo)彈外形，如下圖所示，具體尺寸參見文獻(xiàn)［19］。針對(duì)這種復(fù)雜的框架結(jié)構(gòu)，生成點(diǎn)對(duì)點(diǎn)的結(jié)構(gòu)網(wǎng)格，以此為基礎(chǔ)開展柵格翼導(dǎo)彈亞跨超氣動(dòng)特性的數(shù)值模擬。

圖1 DREV風(fēng)洞試驗(yàn)外形圖Fig.1 Schematic of the DREV grid fin missile

2 網(wǎng)格生成

隨著應(yīng)用環(huán)境的不同，對(duì)柵格翼的氣動(dòng)性能也提出了不同的要求。這就導(dǎo)致柵格翼的結(jié)構(gòu)形狀和包含的格框數(shù)也大不相同，這給結(jié)構(gòu)網(wǎng)格的生成帶來(lái)很大困難。本文采用多塊對(duì)接網(wǎng)格的思想，將柵格翼這種復(fù)雜的升力面系統(tǒng)分成一個(gè)個(gè)的柵格框（如圖2所示），然后分別在不同柵格框內(nèi)生成結(jié)構(gòu)網(wǎng)格，相鄰兩塊公共邊界上的網(wǎng)格點(diǎn)為兩塊公共的網(wǎng)格點(diǎn)。相鄰塊網(wǎng)格點(diǎn)數(shù)相同，為點(diǎn)對(duì)點(diǎn)的對(duì)接，不同子塊的網(wǎng)格共同充滿整個(gè)流場(chǎng)。各個(gè)子塊根據(jù)不同的物面形狀和流場(chǎng)特點(diǎn)采取不同的網(wǎng)格拓?fù)洌跂鸥窨騼?nèi)生成O型網(wǎng)格，其他區(qū)域內(nèi)生成H型網(wǎng)格。隨著來(lái)流速度的不同，柵格框內(nèi)的流態(tài)會(huì)發(fā)生很大變化，復(fù)雜的激波與邊界層干涉，對(duì)流場(chǎng)特性有重大影響。為了保證數(shù)值模擬結(jié)果的準(zhǔn)確性，單個(gè)柵格框內(nèi)要保證一定密度的網(wǎng)格，從而導(dǎo)致整個(gè)柵格翼導(dǎo)彈的網(wǎng)格量比較大。在本文的數(shù)值模擬中，柵格翼導(dǎo)彈半模型結(jié)構(gòu)網(wǎng)格點(diǎn)數(shù)為500萬(wàn)，彈體附近垂直方向網(wǎng)格尺度為4×10－5m，柵格翼上垂直方向網(wǎng)格尺度為2×10－5m。

圖2 柵格翼導(dǎo)彈結(jié)構(gòu)網(wǎng)格分塊示意圖Fig.2 Mesh topology of the grid fin missile

圖3 柵格翼導(dǎo)彈對(duì)稱面網(wǎng)格示意圖Fig.3 Structure grid system on the surface of symmetry

圖4 柵格翼局部網(wǎng)格圖Fig.4 Structure grid schematic of the fin

3 計(jì)算結(jié)果與分析

圖5 Ma＝3，α＝2°條件下殘值收斂史Fig.5 The residual history（Ma＝3，α＝2°）

圖6 Ma＝3，α＝2°條件下柵格翼表面壓力分布Fig.6 Pressure coefficient contours on the fin surface（Ma＝3，α＝2°）

圖7 Ma＝3條件下法向力系數(shù)與試驗(yàn)的比較Fig.7 Normal force coefficient on the grid fin missile（Ma＝3）

圖8 Ma＝3條件下俯仰力矩系數(shù)與試驗(yàn)的比較Fig.8 Pitching moment coefficient on the grid fin missile（Ma＝3）

我們以DREV風(fēng)洞中進(jìn)行大量試驗(yàn)的柵格翼導(dǎo)彈外形為基礎(chǔ)，生成了點(diǎn)對(duì)點(diǎn)的多塊結(jié)構(gòu)網(wǎng)格，并進(jìn)行了不同馬赫數(shù)（Ma＝0.5到Ma＝3）流場(chǎng)特性的數(shù)值模擬研究。其中，我們將部分?jǐn)?shù)值模擬結(jié)果與試驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比分析。圖5給出了Ma＝3、α＝2°狀態(tài)下殘值收斂曲線，從中看出殘值下降了6個(gè)量級(jí)，收斂性較好。圖6給出了圖5狀態(tài)下水平柵格翼表面壓力云圖。圖7給出了Ma＝3條件下法向力系數(shù)計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)的比較。圖8給出了Ma＝3條件下俯仰力矩系數(shù)與試驗(yàn)的比較。圖9給出了Ma＝0.5、1.15、1.5、3條件下軸向力系數(shù)與試驗(yàn)的比較。圖10給出了法向力系數(shù)導(dǎo)數(shù)隨馬赫數(shù)變化與試驗(yàn)結(jié)果的比較。通過對(duì)比分析，以結(jié)構(gòu)網(wǎng)格為基礎(chǔ)，求解N－S方程得到的計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)吻合較好，能夠比較準(zhǔn)確的捕捉柵格翼導(dǎo)彈的流場(chǎng)特性。其中，法向力系數(shù)導(dǎo)數(shù)隨馬赫數(shù)變化規(guī)律與試驗(yàn)結(jié)果一致，在亞聲速段與吻合較好，在跨聲速和超聲速段計(jì)算結(jié)果偏大。圖11和圖12給出了柵格翼內(nèi)部流態(tài)隨馬赫數(shù)增大的變化。當(dāng)來(lái)流馬赫數(shù)為Ma＝0.9時(shí)，在柵格翼通道后面形成了正激波，導(dǎo)致柵格翼升力減小阻力增加；當(dāng)來(lái)流馬赫數(shù)為Ma＝1.5時(shí)，在柵格翼前出現(xiàn)脫體正激波，阻力急劇增大（如圖11所示），波后的氣流是亞聲速的，隨著來(lái)流速度的繼續(xù)增加，激波會(huì)向后移動(dòng)貼在柵格翼的前緣上；當(dāng)來(lái)流馬赫數(shù)達(dá)到2時(shí)，在柵格翼前緣出現(xiàn)斜激波，并且在相鄰柵格壁之間出現(xiàn)了反射，在靠近彈體壁面的柵格框內(nèi)，由于受到壁面附面層的影響，前緣處形成的仍然是正激波；隨著來(lái)流速度的進(jìn)一步加大，達(dá)到馬赫數(shù)3時(shí)，柵格翼前緣形成的斜激波不再落在相鄰的柵格壁上，柵格壁之間不再出現(xiàn)相互干涉，阻力也會(huì)隨之降低，氣動(dòng)性能也有所提升。

圖9 不同馬赫數(shù)下軸向力系數(shù)與試驗(yàn)的比較Fig.9 Axial force coefficient for different Mach numbers

圖10 法向力系數(shù)斜率與試驗(yàn)的比較Fig.10 The slope of normal force coefficient

圖11 柵格翼內(nèi)部馬赫數(shù)云圖（α＝0°）Fig.11 Mach number coefficient contours on the surface of symmetry（α＝0°）

圖12 柵格翼內(nèi)部馬赫數(shù)云圖（α＝0°）Fig.12 Mach number coefficient contours on the surface of symmetry（α＝0°）

4 結(jié) 論

本文通過多塊對(duì)接網(wǎng)格的方法，對(duì)柵格翼導(dǎo)彈這種復(fù)雜構(gòu)型生成了點(diǎn)對(duì)點(diǎn)的結(jié)構(gòu)網(wǎng)格。并以此為基礎(chǔ)，通過求解N－S方程對(duì)柵格翼導(dǎo)彈亞跨超聲速流場(chǎng)進(jìn)行了數(shù)值模擬。通過與試驗(yàn)結(jié)果的對(duì)比分析，本文數(shù)值模擬的結(jié)果在局部流場(chǎng)細(xì)節(jié)和柵格翼導(dǎo)彈氣動(dòng)力方面與試驗(yàn)數(shù)據(jù)吻合較好。說(shuō)明本文采用的結(jié)構(gòu)網(wǎng)格可靠合理，質(zhì)量?jī)?yōu)良，能夠?qū)鸥褚盹w行器復(fù)雜流場(chǎng)進(jìn)行比較精確的數(shù)值模擬，提供準(zhǔn)確的氣動(dòng)力數(shù)據(jù)，為工程應(yīng)用評(píng)估提供了技術(shù)支撐。

［1］DESPIRITO J，EDGE H L，WEINACHT P，et al.Computational fluid dynamic（CFD）analysis of a generic missile with grid fins［R］.ADA 383925，2000.

［2］DESPIRITO J，EDGE H L，et al.CFD analysis of grid fins for maneuvering missiles［R］.AIAA 2000－0391.

［3］DESPIRITO J，SAHU J.Viscous CFD calculations of grid fin missile aerodynamic in the supersonic flow regine［R］.AIAA 2001－0257.

［4］CHEN S，KHALID M，XU H.A comprehensive CFD investigation of grid fins as efficient control surface devices［R］.AIAA 2000－0987.

［5］LIN H，HUANG J C.Navier－Stokes computation for body／cruciform grid fin configuration［R］.AIAA 2002－2722.

［6］RAVINDRA K，NIKHIL V，SHENDE N.Balakrishnan.CFD simulation of the grid fin flows［R］.AIAA 2013－3023.

［7］YAN Z.Drag reduction for sweptback grid fin with blunt and sharp leading edges［J］.Journal of Aircraft，2012，49（5）：1526－1531.

［8］HUGHSON M C，BLADES E L.Transonic aerodynamic analysis of lattice grid tail fin missiles［R］.AIAA 2006－3651.

［9］MA M S，DENG Y Q，ZHENG M，et al.Navier－Stokes computations for a grid fin missile［R］.AIAA 2005－4973.

［10］WU X J，MA M S，DENG Y Q，et al.Navier－Stokes computations of a grid fin missile on hybrid structured－unstructured grids［J］.ACTA Aerodynamica Sinica，2009，27（4）：419－423.（in Chinese）吳曉軍，馬明生，鄧有奇，等.結(jié)構(gòu)／非結(jié)構(gòu)混合網(wǎng)格數(shù)值模擬柵格翼［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2009，27（4）：419－423.

［11］WU X J，CHEN H Q，DENG Y Q，et al.CFD analysis of grid fins［J］.ACTA Aerodynamica Sinica，2009，27（1）：78－82.（in Chinese）吳曉軍，陳紅全，鄧有奇，等.柵格繞流數(shù)值模擬研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2009，27（1）：78－82.

［12］LIU G，XIAO Z Y，JIANG X，et al.Flow simulation of missile with grid fins using structured－unstructured hybrid grids［J］.ACTA Aerodynamica Sinica，2009，25（2）：277－280.（in Chinese）劉剛，肖中云，江雄，等.混合網(wǎng)格方法在柵格翼數(shù)值模擬中的應(yīng)用研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2007，25（2）：277－280.

［13］CAI J S.Numerical study on choked flow over grid－fin configurations［J］.Journal of Spacecraft and Rockets，2009，46（5）：949－956.

［14］DENG F，CHEN S S，YU Y G.Numerical simulation about lift and drag characteristics of grid at subsonic／transonic speeds［J］.Journal of Ballistics，2009，21（4）：34－37.（in Chinese）鄧帆，陳少松，余勇剛.亞跨聲速柵格的升阻特性數(shù)值研究［J］.彈道學(xué)報(bào)，2009，21（4）：34－37.

［15］DENG F，CHEN S S.The investigations of windward swept grid fin′s lift and drag characterization［J］.ACTA Aerodynamica Sinica，2012，30（5）：628－633.（in Chinese）鄧帆，陳少松.前緣后掠式柵格翼升阻特性研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2012，30（5）：628－633.

［16］WANG B，ZHANG W M.Numerical investigation into aerodynamic characteristics of grid fins［J］.Tactical Missile Technology，2011，（2）：31－33.（in Chinese）王斌，張衛(wèi)民.柵格翼氣動(dòng)特性數(shù)值研究［J］.戰(zhàn)術(shù)導(dǎo)彈技術(shù)，2011，（2）：31－33.

［17］WANG X Z，WANG L Q.Numerical simulation of missile with grid fins in supersonic flow based on the overlapping structured grids［J］.AERO Weaponry，2012，（2）：44－47.（in Chinese）王學(xué)占，王立強(qiáng).基于結(jié)構(gòu)重疊網(wǎng)格的柵格翼導(dǎo)彈流場(chǎng)數(shù)值仿真［J］.航空兵器，2012，（2）：44－47.

［18］ZHOU Z，HU F，YU Y.Hypersonic aerodynamic s research of grid fin－body configurations base on surfaces element modeling［J］.Journal of Solid Rocket Technology，2012，35（4）：442－445.（in Chinese）周張，胡凡，于勇.基于面元法的柵格翼翼身組合體高超聲速氣動(dòng)特性研究［J］.固體火箭技術(shù)，2012，35（4）：442－445.

［19］FOUNDER E Y.Wind tunnel investigation of grid fin and conventional planar control surfaces［R］.AIAA 2001－0256.