一種基于新型動態混合重疊網格的數值模擬方法

劉 超，劉仙名，王立強

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

一種基于新型動態混合重疊網格的數值模擬方法

劉 超，劉仙名，王立強

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

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

動態混合重疊網格； 非定常； 網格動態自適應； 三線性插值； 六自由度

0 引 言

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

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

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

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

1 網格技術

1.1 新型重疊網格

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

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

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

圖1 新型重疊網格

1.2 網格動態自適應

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

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

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

自適應加密算法步驟如下：

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

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

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

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

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

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

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

a.在格心處插入新節點；

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

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

d.將新面與原有的面連接形成新網格，并記錄父子關系。

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

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

1.3 三線性插值

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

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

圖2 六面體控制單元

假設頂點為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，型函數可以寫為

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

(1)

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

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

(2)

同理，y,z可求得。

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

2 數值模擬方法

2.1 流動控制方程及求解器

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

(3)

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

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

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

2.2 六自由度剛體運動方程

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

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

(4)

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

(5)

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

3 驗證算例

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

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

圖3 模型比例示意圖

表1 計算條件

外掛物質心位置和線速度變化曲線見圖4～5。總的來說，在所顯示的時間范圍內，計算得到的質心線速度和位移與實驗值都比較吻合。

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

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

圖4 質心位置曲線 圖5 質心速度曲線 圖6 俯仰、偏航、滾轉角曲線

圖7 角速度曲線 圖8 力系數曲線 圖9 力矩系數曲線

4 結 論

本文采用新型混合重疊網格生成方法，運用耦合六自由度方程的非定常數值模擬，以標準投彈算例進行驗證。可得出如下結論：

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

(2) 網格動態自適應可保證每個計算步內主體附近的網格密度，適用于運動狀態的計算；

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

[1] 達興亞, 陶洋, 趙忠良. 基于預估校正和嵌套網格的虛擬飛行數值模擬[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] 楊愛明, 喬志德. 基于運動潛逃網格的前飛旋翼繞流N-S方程數值計算[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] 蔣躍文. 基于廣義網格的CFD方法及其應用[D]. 西安: 西北工業大學，2012.

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

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

TJ760.11

A

1673-5048(2016)06-0044-05