劉 超，劉仙名，王立強(qiáng)

(中國空空導(dǎo)彈研究院，河南 洛陽 471009)

一種基于新型動態(tài)混合重疊網(wǎng)格的數(shù)值模擬方法

劉 超，劉仙名，王立強(qiáng)

(中國空空導(dǎo)彈研究院，河南 洛陽 471009)

外掛導(dǎo)彈發(fā)射時導(dǎo)彈處于載機(jī)復(fù)雜干擾流場中，會導(dǎo)致空氣動力的非定常、非線性特性。此時運(yùn)用常規(guī)風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)插值作為工程氣動數(shù)據(jù)基礎(chǔ)在一定程度上存在誤差。本文提出一種基于計算流體力學(xué)(CFD)的非定常數(shù)值模擬方法，將飛行力學(xué)方程與空氣動力學(xué)方程耦合求解，更精確地模擬真實(shí)的飛行狀態(tài)。所采用的新型動態(tài)混合重疊網(wǎng)格相比于以往重疊網(wǎng)格具有生成簡便、針對運(yùn)動過程可自適應(yīng)調(diào)整等優(yōu)點(diǎn)。以標(biāo)準(zhǔn)投彈算例作為驗(yàn)證算例，結(jié)果體現(xiàn)了本文網(wǎng)格和數(shù)值模擬方法的工程應(yīng)用價值。

動態(tài)混合重疊網(wǎng)格； 非定常； 網(wǎng)格動態(tài)自適應(yīng)； 三線性插值； 六自由度

0 引 言

現(xiàn)代戰(zhàn)斗機(jī)的外掛導(dǎo)彈武器，在發(fā)射初始階段對載機(jī)的安全性有重要影響。發(fā)射導(dǎo)彈時，導(dǎo)彈處于載機(jī)復(fù)雜干擾流場中，會出現(xiàn)空氣動力的非定常、非線性甚至非對稱特性，引起氣動和運(yùn)動的交叉耦合[1]，此時往往出現(xiàn)非正常分離，如若導(dǎo)彈與載機(jī)發(fā)生碰撞，則會嚴(yán)重威脅載機(jī)與飛行人員安全。導(dǎo)彈氣動設(shè)計師必須對外掛導(dǎo)彈發(fā)射時的氣動特性加以研究，以確保發(fā)射條件不受太大影響，保證載機(jī)安全。世界各航空大國均投入大量人力和物力，從計算和實(shí)驗(yàn)兩方面研究多體干擾及多體分離問題。此外，為了有效擊中目標(biāo)，需要掌握投放初始階段導(dǎo)彈姿態(tài)和運(yùn)動軌跡數(shù)據(jù)。因此，對導(dǎo)彈與載機(jī)分離進(jìn)行數(shù)值模擬，研究導(dǎo)彈發(fā)射過程中的機(jī)彈氣動力干擾，考慮其對導(dǎo)彈發(fā)射后初始彈道的影響十分必要。

由一般風(fēng)洞實(shí)驗(yàn)和計算流體力學(xué)(Computational Fluid Dynamics，CFD)方法得到的局部線性化氣動力數(shù)據(jù)，在運(yùn)動狀態(tài)較為復(fù)雜的情況下不能準(zhǔn)確描述飛行器的非定常、非線性氣動力和運(yùn)動規(guī)律。隨著CFD方法的發(fā)展，綜合運(yùn)用空氣動力學(xué)、飛行力學(xué)耦合一體化數(shù)值模擬技術(shù)是研究和解決這類問題的一個新途徑。其特殊的研究方式，也能在一定程度上填補(bǔ)靜態(tài)風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)以及靜態(tài)CFD計算數(shù)據(jù)和飛行實(shí)驗(yàn)之間的數(shù)據(jù)空白[2]。

重疊網(wǎng)格(overlapping grids)也叫Chimera網(wǎng)格[3]或者嵌套網(wǎng)格(overset grids)，該類網(wǎng)格的實(shí)質(zhì)是一種網(wǎng)格的區(qū)域分割和組合策略[4]，其主要優(yōu)點(diǎn)是降低了網(wǎng)格生成的難度，提高了網(wǎng)格生成的靈活性，保證了原始網(wǎng)格的質(zhì)量等。而重疊網(wǎng)格方法特別適用于復(fù)雜外形和多體相對運(yùn)動，特別是20世紀(jì)90年代以來，動態(tài)重疊網(wǎng)格技術(shù)的提出和發(fā)展[5-6]，拓展了外掛物分離等狀態(tài)的研究途徑，國外在這方面已經(jīng)取得很大成就[7-8]。

本文綜合運(yùn)用笛卡爾網(wǎng)格和結(jié)構(gòu)重疊網(wǎng)格，提出一種新型重疊網(wǎng)格。運(yùn)用基于CFD的非定常數(shù)值模擬方法，將飛行力學(xué)方程、空氣動力學(xué)方程耦合求解，以標(biāo)準(zhǔn)投彈算例為驗(yàn)證算例。

1 網(wǎng)格技術(shù)

1.1 新型重疊網(wǎng)格

現(xiàn)階段的重疊網(wǎng)格多以部件作為基本單元，例如在標(biāo)準(zhǔn)投彈算例中，對導(dǎo)彈做一套體網(wǎng)格，對機(jī)翼做一套體網(wǎng)格，兩套體網(wǎng)格之間存在相互運(yùn)動。這種形式的體網(wǎng)格有兩個不足之處：(1)若部件外形較為復(fù)雜，則體網(wǎng)格的生成較為困難。(2)若部件之間相對運(yùn)動時間較長，距離較遠(yuǎn)，部件之間可能會出現(xiàn)網(wǎng)格不匹配，導(dǎo)致計算無法繼續(xù)。

為解決上述問題，提出一種新型重疊網(wǎng)格。這種重疊網(wǎng)格在生成線和面時就進(jìn)行重疊，把復(fù)雜幾何外形結(jié)構(gòu)的整體拆分成簡單幾何外形結(jié)構(gòu)的個體，分別生成面網(wǎng)格，每片面網(wǎng)格單獨(dú)外推生成附面層網(wǎng)格(物面附近的體網(wǎng)格)。相鄰面網(wǎng)格、體網(wǎng)格之間存在的重疊部分用作信息交換。用此方法對幾何外形進(jìn)行剖分，可以簡單生成任意復(fù)雜外形結(jié)構(gòu)的網(wǎng)格拓?fù)洹?/p>

如圖1所示，(a)為面網(wǎng)格，由彈翼翼梢面網(wǎng)格(有一部分“長”到彈翼翼面上形成重疊)、彈翼翼面面網(wǎng)格(有一部分“長”到彈體面網(wǎng)格上形成重疊)、彈體面網(wǎng)格三部分組成。(b)(c)(d)為三部分分別外推形成的體網(wǎng)格，(e)為(b)(c)(d)的組合體，(f)為全彈的體網(wǎng)格。

圖1 新型重疊網(wǎng)格

1.2 網(wǎng)格動態(tài)自適應(yīng)

僅有物面附近的網(wǎng)格不足以滿足氣動力計算的要求，需有遠(yuǎn)場網(wǎng)格。選用笛卡爾網(wǎng)格作為背景網(wǎng)格具有生成簡單、便于挖洞操作、易于做網(wǎng)格自適應(yīng)等優(yōu)點(diǎn)。

在CFD中，流動控制方程在離散過程中產(chǎn)生誤差是不可避免的。提高離散格式的精度和加密網(wǎng)格將有助于減小誤差。然而提高格式的精度和全區(qū)域的加密網(wǎng)格往往會帶來存儲量和計算量的較大增加，而自適應(yīng)技術(shù)可有效地解決這一問題[9]。尤其在遇到多體相對運(yùn)動的問題時，隨著物體相對位置的改變，網(wǎng)格加密部分也需要進(jìn)行相應(yīng)的改變，以保證計算的效率和準(zhǔn)確性。

網(wǎng)格的動態(tài)自適應(yīng)是根據(jù)流場的計算結(jié)果和計算對象的位置信息，實(shí)時實(shí)地或人為指定每n真實(shí)時間步對網(wǎng)格進(jìn)行加密或稀疏。本文的背景網(wǎng)格為笛卡爾網(wǎng)格，借鑒笛卡爾網(wǎng)格的剖分思想引入父子關(guān)系。網(wǎng)格自適應(yīng)分為自適應(yīng)加密算法和自適應(yīng)稀疏算法。

自適應(yīng)加密算法步驟如下：

(1) 根據(jù)梯度大小(壓力、馬赫數(shù)、溫度等)對流場區(qū)域進(jìn)行劃分，若梯度值符合要求，便將該類網(wǎng)格標(biāo)記為A，對標(biāo)記為A的網(wǎng)格，若相鄰有多于兩個A網(wǎng)格則對此網(wǎng)格標(biāo)記為“加密”。

(2) 將標(biāo)記為“加密”網(wǎng)格的所有面標(biāo)記為“加密”，處理被加密的面，將標(biāo)記為“加密”面的所有邊標(biāo)記為“加密”，在邊的中點(diǎn)插入新節(jié)點(diǎn)，形成邊上的懸點(diǎn)，生成兩條邊，記錄邊的父子關(guān)系。

(3) 對每個面進(jìn)行循環(huán)，如果該面標(biāo)記為“加密”，則：

a.在面心處插入新節(jié)點(diǎn)，形成面上懸點(diǎn)；

b.將新節(jié)點(diǎn)與每條邊上懸點(diǎn)連接，形成新邊；

c.將新邊與原有的邊連接形成新面，并記錄父子關(guān)系。

(4) 對每個網(wǎng)格進(jìn)行循環(huán)，如果該網(wǎng)格標(biāo)記為“加密”，則：

a.在格心處插入新節(jié)點(diǎn)；

b.將新節(jié)點(diǎn)與每個面上懸點(diǎn)連接，形成新邊；

c.將新邊與原有邊連接形成新面；

d.將新面與原有的面連接形成新網(wǎng)格，并記錄父子關(guān)系。

(5) 處理未標(biāo)記為“加密”但有懸點(diǎn)的邊的面，生成新的子面，處理未標(biāo)記為“加密”但有懸點(diǎn)的面的網(wǎng)格，生成新的子網(wǎng)格。

自適應(yīng)稀疏算法與自適應(yīng)加密算法是相交的過程，不再詳述。經(jīng)過自適應(yīng)加密以及稀疏之后的笛卡爾網(wǎng)格能很好地適應(yīng)計算條件，優(yōu)化計算速度以及計算準(zhǔn)確性。

1.3 三線性插值

重疊網(wǎng)格須通過插值的方式實(shí)現(xiàn)網(wǎng)格之間的數(shù)據(jù)交換。交換數(shù)據(jù)的重要形式是質(zhì)量、能量和動量通量。如果流動中存在激波或大梯度物理區(qū)跨越插值洞邊界時，須使用守恒型插值方法來正確捕捉激波。對基于密度求解法的可壓縮流動來說，國外研究者[4]認(rèn)為三線性插值方法具有保證正常精度的能力，這一結(jié)論已經(jīng)得到許多數(shù)值實(shí)驗(yàn)的證實(shí)。

本文使用三線性插值來進(jìn)行重疊網(wǎng)格之間的數(shù)據(jù)交換。引入型函數(shù)思想來解決重疊區(qū)域三線性插值問題。精度為二階精度，能適應(yīng)任何形狀的插值問題。該方法計算簡單，對重疊網(wǎng)格的計算結(jié)果較好。三維問題的控制體為類似于六面體的單元，如圖2所示。

圖2 六面體控制單元

假設(shè)頂點(diǎn)為Ai(i=1, 2, …, 8)，按逆時針進(jìn)行排序。對其進(jìn)行局部坐標(biāo)變換，變換到(ξ,η,ζ)上的單元立方體。變換后頂點(diǎn)Ai(i=1, 2, …, 8)在(ξ,η,ζ)中的坐標(biāo)為A1(0, 0, 0)；A2(1, 0, 0)；A3(1, 1, 0)；A4(0, 1, 0)；A5(0, 0, 1)；A6(1, 0, 1)；A7(1, 1, 1)；A8(0, 1, 1)。因?yàn)轫旤c(diǎn)Ai(i=1, 2, …, 8)的局部坐標(biāo)系(ξ,η,ζ)取值為0或1，型函數(shù)可以寫為

φi=(-1)1+ξi+ηi+ζi(ξ+ξi-1)(η+ηi-1)(ζ+ζi-1) (i=1, 2, …, 8)

(1)

上式為三線性，即對于ξ,η,ζ都是線性的。

然后進(jìn)行坐標(biāo)變換，目的是由網(wǎng)格點(diǎn)(x,y,z)求當(dāng)前面網(wǎng)格局部坐標(biāo)(ξ,η,ζ), 其中：

(2)

同理，y,z可求得。

這三組數(shù)據(jù)組成三個方程、三個未知變量的非線性方程組，采用牛頓迭代法求解該非線性方程組, 即可得到(ξ,η,ζ)。

2 數(shù)值模擬方法

2.1 流動控制方程及求解器

運(yùn)用三維非定?？蓧嚎sNavier-Stokes方程為流動控制方程。一般曲線坐標(biāo)系中，無量綱化的方程守恒形式為[10]

(3)

式中:Q為守恒矢量；F,G,H為對流通量；Fv,Gv,Hv為粘性通量；ξ,η,ζ為三個貼體坐標(biāo)系方向；t為時間；Re∞為自由來流雷諾數(shù)。

流場求解采用基于結(jié)構(gòu)網(wǎng)格的Roe Upwind格式和van Albada限制器。湍流模型選擇Spalart-Allmaras一方程模型，該模型是從量綱分析和經(jīng)驗(yàn)出發(fā)，針對簡單流動逐步發(fā)展起來的模型，其容錯功能好，處理復(fù)雜流動能力強(qiáng)，同時相對于兩方程模型計算量小、穩(wěn)定性好，有較高精度。

時間推進(jìn)采用雙時間步法，為了加速收斂，使用網(wǎng)格分塊算法和網(wǎng)格分塊并行計算技術(shù)。

2.2 六自由度剛體運(yùn)動方程

本文算例的變形相對飛行器的尺寸較小，其運(yùn)動可以近似采用剛體來描述。而剛體的運(yùn)動可以分為兩部分：質(zhì)心的運(yùn)動和物體繞質(zhì)心的轉(zhuǎn)動。

剛體動力學(xué)方程如下[11]：

(4)

剛體運(yùn)動學(xué)方程如下[11]：

(5)

式中：m為剛體質(zhì)量；V為速度；P為推力；X，Y，Z分別為阻力、升力、側(cè)向力；α,β分別為攻角、側(cè)滑角； ?為彈道傾角；Ψv為彈道偏角；θ,Ψ,φ分別為俯仰角、偏航角、滾轉(zhuǎn)角；rv為速度滾轉(zhuǎn)角；J為轉(zhuǎn)動慣量；ω為角速度；M為力矩。

3 驗(yàn)證算例

為了驗(yàn)證本文方法，以實(shí)驗(yàn)數(shù)據(jù)較為全面的標(biāo)準(zhǔn)彈翼分離作為驗(yàn)證算例。該算例來源于美國空軍實(shí)驗(yàn)室進(jìn)行的外掛物分離實(shí)驗(yàn)研究，風(fēng)洞實(shí)驗(yàn)在阿諾德工程發(fā)展中心(AEDC)4英尺跨音速空氣動力風(fēng)洞進(jìn)行，具有比較完備的實(shí)驗(yàn)數(shù)據(jù)[3]。機(jī)翼翼型為NACA 64A010。

模型的尺寸如圖3所示，計算條件如表1所示。為了直觀顯示結(jié)果數(shù)據(jù)，彈體坐標(biāo)系定義為：坐標(biāo)原點(diǎn)為導(dǎo)彈質(zhì)心，X軸沿導(dǎo)彈縱軸指向前方，Z軸為0時刻重力方向，Y軸符合右手定則，0時刻地軸系與彈體坐標(biāo)系重合。

圖3 模型比例示意圖

表1 計算條件

外掛物質(zhì)心位置和線速度變化曲線見圖4～5?？偟膩碚f，在所顯示的時間范圍內(nèi)，計算得到的質(zhì)心線速度和位移與實(shí)驗(yàn)值都比較吻合。

外掛物的姿態(tài)和角速度的時間歷程見圖6～7，p，q，r分別為彈體系XYZ軸的三個角速度分量，φ，θ和Ψ分別為滾轉(zhuǎn)角、俯仰角和偏航角，用來表示從地軸系到彈體系的變換。對于俯仰運(yùn)動來說，在分離初始階段，彈射力產(chǎn)生較大的抬頭力矩，彈射力維持的時間到t=0.055 s，之后彈射力作用消失，作用在外掛物上的氣動力和力矩開始起主導(dǎo)作用，這時候氣動力產(chǎn)生的是低頭力矩，在t=0.2 s以后，外掛物開始下俯運(yùn)動，在偏航方向和滾轉(zhuǎn)方向，計算的偏航角和滾轉(zhuǎn)角與實(shí)驗(yàn)值吻合良好。

圖8～9分別給出了氣動力和氣動力矩的時間歷程比較。從氣動力方面來看，除了初始的Cy計算值略小于實(shí)驗(yàn)值之外，其他方面都吻合良好； 對氣動力矩來說，滾轉(zhuǎn)力矩和偏航力矩吻合較好，初始俯仰力矩Cmy略小于實(shí)驗(yàn)值。俯仰力矩的偏差直接導(dǎo)致俯仰角θ出現(xiàn)誤差。

圖4 質(zhì)心位置曲線 圖5 質(zhì)心速度曲線 圖6 俯仰、偏航、滾轉(zhuǎn)角曲線

圖7 角速度曲線 圖8 力系數(shù)曲線 圖9 力矩系數(shù)曲線

4 結(jié) 論

本文采用新型混合重疊網(wǎng)格生成方法，運(yùn)用耦合六自由度方程的非定常數(shù)值模擬，以標(biāo)準(zhǔn)投彈算例進(jìn)行驗(yàn)證?？傻贸鋈缦陆Y(jié)論：

(1) 在復(fù)雜繞流的CFD計算中，新型混合重疊網(wǎng)格能夠降低網(wǎng)格生成難度，保證主體每一部分網(wǎng)格質(zhì)量；

(2) 網(wǎng)格動態(tài)自適應(yīng)可保證每個計算步內(nèi)主體附近的網(wǎng)格密度，適用于運(yùn)動狀態(tài)的計算；

(3) 外掛物投放過程數(shù)值模擬結(jié)果表明，本文提出的新型混合重疊網(wǎng)格算法為外掛物投放分離安全分析預(yù)測和優(yōu)化設(shè)計提供了一種高效的解決方案，具有較高的工程應(yīng)用價值。

[1] 達(dá)興亞, 陶洋, 趙忠良. 基于預(yù)估校正和嵌套網(wǎng)格的虛擬飛行數(shù)值模擬[J]. 航空學(xué)報, 2012, 33(6): 977-983.

[2] 陶洋, 范召林, 吳繼飛. 基于CFD的方形截面導(dǎo)彈縱向虛擬飛行模擬[J]. 力學(xué)學(xué)報, 2010, 42(2): 169-176

[3] Noack R W, Slotnick J P. A Summary of the 2004 Overset Symposium on Composite Grids and Solution Technology[C]∥43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 2005.

[4] Meakin R L. On the Spatial and Temporal Accuracy of Overset Grid Methods for Moving Body Problems[C]∥12th Applied Aerodynamics Conference, Fluid Dynamics and Co-Located Conferences, Colorado Springs, Colorado, 1994.

[5] Lee S, Park M, Cho K W, et al. A New Automated Chimera Method for Prediction of Store Trajectory[C]∥17th Applied Aerodynamics Conference, Fluid Dynamics and Co-Located Conferences, Norfolk, Virgina, 1999.

[6] 楊愛明, 喬志德. 基于運(yùn)動潛逃網(wǎng)格的前飛旋翼繞流N-S方程數(shù)值計算[J]. 航空學(xué)報, 2001, 22(5): 434-436.

[7] Dougherty F C, Kuan J H. Transonic Store Separation Using a Three-Dimensional Chimera Grid Scheme[C]∥27th Aerospace Sciences Meeting, Reno, Nevada, 1989.

[8] Snyder D O, Koutsavdis E K, Anttonen J S R. Transonic Store Separation Using Unstructured CFD with Dynamic Meshing[C]∥33rd AIAA Fluid Dynamics Conference and Exhibit, Fluid Dynamics and Co-Located Conferences, Orlando, Florida, 2003.

[9] 蔣躍文. 基于廣義網(wǎng)格的CFD方法及其應(yīng)用[D]. 西安: 西北工業(yè)大學(xué)，2012.

[10] 閆超. 計算流體力學(xué)方法及應(yīng)用[M]. 北京：北京航空航天大學(xué)出版社，2006: 18-25.

[11] 李新國，方群.有翼導(dǎo)彈飛行動力學(xué)[M]. 西安：西北工業(yè)大學(xué)出版社，2008: 48-49.

A Numerical Simulation Method Based on the New Dynamic Mixing Overlapping Grid

Liu Chao, Liu Xianming, Wang Liqiang

(China Airborne Missile Academy，Luoyang 471009，China)

During launching of the external-carried missile, the missile is in a complex flow field of airborne, which will cause the unsteady and nonlinear characteristics of aerodynamic. At this time, using common method of wind tunnel test data interpolation as the aerodynamic data base will create error to a certain extent. An unsteady numerical simulation method based on the computational fluid dynamics (CFD) is proposed, and the flight mechanics equations and aerodynamic equations are coupled and solved, so as to simulate the real flight status more accurately. Compared with the past overlapping grids, the used new dynamic mixing overlapping grid has advantages such as easy generation, selfadapting adjustment for moving process. Taking the standard bomb-dropping example as a verification example, the results indicate the engineering importance of the dynamic mixing overlapping grid and the numerical simulation method.

dynamic mixing overlapping grid； unsteady； dynamic self-adapting； tri-linear interpolation； six-degree of freedom

10.19297/j.cnki.41-1228/tj.2016.06.010

2016-01-05

劉超(1991-)，男，黑龍江牡丹江人，碩士，研究方向?yàn)轱w行器設(shè)計。

TJ760.11

A

1673-5048(2016)06-0044-05