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

劉 超，劉仙名，王立強

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

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

劉 超，劉仙名，王立強

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

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

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

0 引 言

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

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

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

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

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

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

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

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

如圖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)格不足以滿足氣動力計算的要求，需有遠場網(wǎng)格。選用笛卡爾網(wǎng)格作為背景網(wǎng)格具有生成簡單、便于挖洞操作、易于做網(wǎng)格自適應(yīng)等優(yōu)點。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 三線性插值

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

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

圖2 六面體控制單元

假設(shè)頂點為Ai(i=1, 2, …, 8)，按逆時針進行排序。對其進行局部坐標變換，變換到(ξ,η,ζ)上的單元立方體。變換后頂點Ai(i=1, 2, …, 8)在(ξ,η,ζ)中的坐標為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)。因為頂點Ai(i=1, 2, …, 8)的局部坐標系(ξ,η,ζ)取值為0或1，型函數(shù)可以寫為

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

(1)

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

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

(2)

同理，y,z可求得。

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

2 數(shù)值模擬方法

2.1 流動控制方程及求解器

運用三維非定常可壓縮Navier-Stokes方程為流動控制方程。一般曲線坐標系中，無量綱化的方程守恒形式為[10]

(3)

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

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

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

2.2 六自由度剛體運動方程

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

剛體動力學方程如下[11]：

(4)

剛體運動學方程如下[11]：

(5)

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

3 驗證算例

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

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

圖3 模型比例示意圖

表1 計算條件

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

[2] 陶洋, 范召林, 吳繼飛. 基于CFD的方形截面導彈縱向虛擬飛行模擬[J]. 力學學報, 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] 楊愛明, 喬志德. 基于運動潛逃網(wǎng)格的前飛旋翼繞流N-S方程數(shù)值計算[J]. 航空學報, 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è)大學，2012.

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

[11] 李新國，方群.有翼導彈飛行動力學[M]. 西安：西北工業(yè)大學出版社，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-)，男，黑龍江牡丹江人，碩士，研究方向為飛行器設(shè)計。

TJ760.11

A

1673-5048(2016)06-0044-05