基于動態網格的多體分離計算技術研究

徐琳 宋萬強

摘要：多體分離問題是航空航天等工業設計部門非常關心的問題之一。目前，多體分離軌跡問題主要依托風洞試驗和數值模擬兩類。本文基于網格變形與局部重構相結合的動態網格方法，在單步小位移的情況下用Laplace算法進行網格變形.當總的位移量逐漸增大至網格質量低于設定標準時，進行局部網格重新生成，采用非結構網格、Euler方程、準定常的方法開展外掛投放數值求解，模擬了典型的多體分離標模（機翼/外掛物分離），并與AEDC風洞試驗結果進行了對比，驗證了本文多體分離計算技術的準確性。

關健詞：動態網格;多體分離;數值模擬;網格變形;Laplace

中圖分類號：V224.5 文獻標識碼：A

多體分離過程相當復雜，不同物體之間存在相互干擾，為確保飛行器的飛行安全，必須研究他們之間的干擾與分離，因此，多體分離問題是航空航天等工業設計部門非常關心的力學問題之一[1]。目前，多體分離軌跡問題主要依托風洞試驗和數值模擬兩類。風洞試驗采用的技術方法主要有捕獲軌跡系統（CTS）、自由投放試驗、網格掃描方法和流向角方法等?？傮w來說，地面設備不足使得采用風洞試驗手段研究多體分離問題相對比較困難，尤其是在高速機動飛行等非常規情況下，試驗研究的局限性越發明顯[2，3]。近年來，數值方法的進步使得仿真計算越來越能逼近現實狀況，采用數值方法模擬多體分離軌跡也逐漸成為一種常用的技術手段[4]。

多個物體分離時各模型邊界之間的相對運動幅度較大，這對于采用數值模擬方法進行分離軌跡求解的計算網格提出了更高的要求，也是其所需解決的關鍵問題。目前國內外比較成熟的動態網格求解技術大致可以分為三種。

（1）嵌套網格技術

嵌套網格技術由運動物體網格和背景網格組成，通過對重疊區插值運算完成網格間的數據傳遞計算。此類方法為國外常用的運動問題求解算法，相對成熟。其缺點是頻繁的插值運算會帶來解的精度損失。

（2）網格變形技術

隨著模型相對位置的變化，計算網格會按照設定的變形算法進行壓縮或拉伸變形，以保證網格與模型邊界的貼合。常用算法有彈簧法、彈性體方法Delaunay圖法和RBF網格變形方法等。此類方法多用于非結構網格，其缺點是當相對運動較大時，網格質量無法保證。

（3）網格重生技術

當模型相對位置發生變化時，網格依據新的模型重新生成，并通過插值將原網格的物理量傳遞給新的網格。此方法能保證求解過程中的網格質量，其缺點是頻繁的插值運算會帶來解的精度損失，且網格重生會耗時較長，影響效率[5]。

本文采用網格變形與局部重構相結合的動態網格方法：在單步小位移的情況下用Laplace[6，7]算法進行網格變形，當總的位移量逐漸增大至網格質量低于設定標準時，進行局部網格重新生成。采用這種動態網格方法可以計算一個或多個運動物體的軌跡。本文采用非結構網格、Euler方程、準定常的方法開展外掛投放數值求解，模擬了典型的多體分離標模（機翼/外掛物分離），并與阿諾德工程發展復雜綜合體（AFDC）風洞試驗結果183進行了對比。

1 軌跡計算方法

本文運用運動網格變形/局部重構計算的動態網格求解方法，對機翼外掛物分離軌跡進行數值模擬。具體的計算過程如圖1所示。

具體為：（1）首先利用非結構網格生成技術生成背景網格和分離部件子網格，并設置網格變形及軌跡計算相關參數;（2）然后利用“挖洞”技術進行局部網格重構，生成背景網格與分離部件一體的初始計算網格，并計算定常流場作為初始流場;（3）由于計算的流場力為定常氣動力，需增加氣動阻尼進行運動物體受力修正;（4）求解六自由度方程，計算分離物體的下一位置和姿態;（5）移動分離物體，進行網格變形，進行網格質量檢查，若網格質量不能滿足計算要求，則利用“挖洞”技術進行局部網格重構，生成新的網格;（6）采用歐拉方程進行數值模擬，求取分離物體的受力分布;（7）判斷氣動力是否收斂，如果未收斂則繼續進行內迭代，如果收斂則進行下一步;（8）重復（3）～（7）步，直至軌跡模擬求解結束。

2飛行力學模型

本文中使用到的飛行力學模型是建立在剛體六自由度運動基礎上[9，10]，涉及到的運動方程如下所示：

L'_C=ω'I'_C+ω'×（ω'I'_C）

Ω'=ω'（4）式中：F為運動物體上的合力;m為指運動物體的質量;r_C（t）和r_C（t）為指速度及質心的位置;L_C為運動物體上相對于物體質心的合力矩;ω為角速度;I_C相對于運動物體質心的慣性張量;t為當前時間步;t⁰為上一個時間步;Ω為旋轉角速度。上面的公式可以寫為：

q=f（q）（6）式中：f是相對于右手坐標系。

A（5）、A（6）是運用五階Runge-Kutta方法進行求解，這種方法在多體分離的求解中可以提供一個較高精度的解。在準定常求解中，每一個時刻的流場求解運動物體的運動都是定常的。為了把運動物體的運動影響加入計算結果中，本文在計算運動物體受力中增加氣動阻尼：式中：C_lp，C_mq和C_nr為氣動阻尼系數;P，Q和R為運動參考面內的滾轉、俯仰和偏航速度;q_dynr為運動物體上的氣動壓力;S_r，c_r和b_r分別為運動物體的參考面積、參考弦長和參考展長。

3 網格變形/局部，構技術基本原理

本文采用網格變形與局部重構相結合的方法：在單步小位移的情況下用Laplace算法進行網格變形，當總的位移量逐漸增大至網格質量低于設定標準時，進行局部網格重新生成。在小位移的情況下，網格變形就是通過移動網格點的位置來貼合運動物體的邊界。最常見的方法就是拉普拉斯光順，使如式（8）所示的二次方程最小化：式中：δ為某一坐標方向的位移;n_i，j^ref為邊矢量;I_i，j^ref為與該邊矢量相關的對偶網格曲面矢量。這個二次方程最小化的求解是對位移場的每一項單獨求解，也就是說忽略了不同坐標矢量之間的耦合。

式（9）為網格拉伸率P的定義。當網格沒有負體積，且網格拉伸率ρ滿足：ρ_min<ρ<ρ_max。其中，最小網格拉伸率ρ_max與最大網格拉伸率ρ_max為用戶設定。當總的位移量逐漸增大至使網格質量低于設定標準時，進行局部網格重新生成。網格拉伸率P的定義如圖2所示。

當通過網格變形得到新的網格質量無法滿足要求時，則會進行局部網格重構，生成新的網格，再進行氣動計算。局部網格重構是指根據用戶定義的具體范圍，對運動背景網格（多指機體網格）進行挖洞，將運動物體及其周圍網格（多指彈體網格）根據運動物體運動位置填入機體網格，并將機體網格與彈體網格之間局部重構新的網格。具體流程為：（1）首先根據用戶自定義的一個或者多個outer box對背景網格進行挖洞（outer box可以穿越壁面），保留用戶定義的box外的網格;（2）根據用戶定義的一個或多個innerbox，保留box內運動物體周圍的網格，讓其隨運動物體一起運動（inner box不可以穿越壁面，需將運動物體全部包含在內，如果運動物體周圍有加密網格，需將運動物體周圍的加密網格包含在內）;（3）outer box之內inner box之外的部分生成新的網格;（4）最后將三部分網格合成新的網格，完成局部網格重構。如圖3所示。

4 算例驗證

4.1 二維雙翼型算例

為了高效地驗證本文軌跡計算的能力，方便直觀地展示計算過程中的網格變形情況，本文首先采用二維雙翼型算例。本次計算的二維雙翼型算例是由兩個NACA0012翼型構成，初始位置如圖4（a）所示，上邊的翼型充當飛機，下面的翼型充當導彈。本次計算來流馬赫數為0.8，雷諾數為1.8×10⁷，下翼型重量（質量）5000kg，當地重力加速度為9.8m/s²，橫搖阻尼系數取-4.0/rad，縱搖阻尼系數取-40.0/rad，偏航阻尼系數取-40.0/rad。流場求解采用歐拉方程準定常計算。

圖4為雙翼型軌跡模擬過程網格變形/重構圖。圖4（a）為采用非結構網格生成技術生成的背景網格以及運動物體子網格;圖4（b）為初始位置網格，它是由背景網格和運動物體子網格通過上文所講的局部網格重構技術生成的;圖4（c）為小位移情況下，通過拉普拉斯光順的方法進行網格變形，從圖4（c）可以清楚地看出，下面的翼型上翼面網格被拉伸，下翼面網格被壓縮;圖4（d）為大變形情況下，網格變形已經不能滿足計算對網格質量的要求，通過局部網格重構生成的該位置新的網格。圖5為雙翼型軌跡模擬過程壓力分布云圖。

4.2 AEDC算例計算

4.2.1 模型幾何參數

AEDC外掛物標模算例的捕獲軌跡系統（CTS）試驗是由AEDC于1990年完成的。其模型，即機翼/掛架/帶舵外掛物模型（Wing/Pylon/Finned-Store，WPSF）是美國發起第一次分離投放計算流體力學（CFD）驗證時使用的一個機翼和導彈的簡化模型[6]。如圖6所示。很多研究都選用AEDC模型做為運動軌跡數值模擬的算例，是因為這個模型的幾何數據和試驗數據都比較詳細。

AEDC模型由一個帶有雙掛架的三角翼和一個簡易導彈組成。其中，三角翼的展向截面為NACA-64A010翼型，機翼傾斜角為45°;機翼的掛架由一個對稱的弧面一平面一弧面組合而成，掛架的頂端與三角翼前緣的距離為0.61m，其對稱線與三角翼對稱線的展向距離為3.3m;導彈由一個旋轉體與4個對稱尾翼組合而成，其旋轉體的半徑為0.25m，總長度為3.017m。導彈和掛架之間的初始距離為3.66cm。在導彈分離的前0.054s內，導彈質心前后分別受到兩個大小為10.7kN和42.7kN的彈射力作用，直至0.054s后導彈離開掛架，作用力消失。表1為AEDC外掛物標模的其他具體參數。AEDC外掛物標模受力示意圖如圖7所示。

4.2.2 初始網格和流場參數

本次計算共準備兩套網格：背景網格wing.grid和運動物體網格store.grid。這兩套網格采用Pointwise繪制，其中wing.grid共包含131487個節點和787942個單元，store.grid共包含27891個節點和143945個單元，由于本次軌跡計算采用歐拉方程求解，故并沒有繪制棱柱層進行加密，僅對導彈軌跡區域附近進行了網格局部加密。具體網格示意圖如圖8所示。

本次計算來流馬赫數為0.95，單位雷諾數為7.874×10⁶，迎角為0°，飛行高度是8000m。流場求解采用歐拉方程準定常計算，時間步長△t為0.002s，求解總時長1s。

4.2.3 模擬結果

通過網格局部重構的方法生成的網格如圖9所示，為了保證導彈分離過程中對導彈流體計算的準確性，需保證計算過程中導彈周圍網格的疏密度，因此，在設置時，將導彈周圍加密網格區域在inner box的范圍內，如圖9所示導彈周圍的網格為導彈初始網格的保留。

初始時刻定常流場計算的壓力（壓強）分布如圖10所示，由于流動是跨聲速的，因此，在外掛物中部和尾部區域出現激波。由于吊艙和掛架很近，二者之間有較強的氣動干擾，使這一區域的流動比較復雜。圖11為質心位移與試驗值的對比圖。從圖11可以看出，采用準定常流動計算方法，實現了AEDC機翼/導彈分離標準模型軌跡模擬，Is時間內導彈投放軌跡計算結果與風洞試驗數據基本吻合。

5 結論

本文通過采用網格變形與局部重構相結合的網格變形技術，以NACA0012和AFDC為例，耦合求解Euler方程與六自由度方程，采用準定常的方法實現了機翼外掛物分離軌跡數值模擬，證明本文網格變形/局部網格重構方法有很好的實用性。AEDC仿真結果與試驗值吻合較好，驗證了本文多體分離計算方法的準確性。

參考文獻

[1]唐志共，李彬，鄭鳴，等.飛行器外掛投放數值模擬[J].空氣動力學學報，2009，27（5）：592-596.

[2]曾錚，王剛，葉正寅.RBF整體網格變形技術與多體軌跡仿真[J]，空氣動力學學報，2015，33（2）：170-177.

[3]Torsten B，Lars T.Time-accurate CFD approach to numericalsimulation of store separation trajectory prediction[R].AIAA2011-3958，2011.

[4]Mahmood T，Aizud M N，Zahir S.Aerodynamic effects of thestore release onthe roll attitude of a wing configuration[C]//Transonic Flight Proceedings of International BhurbanConference on Applied Sciences&Technology，Islamabad，Pakistan，2011.

[5]聶璐，向錦武，飛機外掛物投放安全性的參數影響分析陰.飛行力學，2011，29（2）：25-28.

[6]黃鵬飛.拉普拉斯加權聚類算法的研究[D].南京：南京航空航天大學，2009.

[7]胡勇全.基于自動生成控制點的拉普拉斯變形[D].大連：大連理工大學，2012.

[8]Elias E P，Spyridon D K.CID transonic store separationtrajectory predictions with comparison to wind tunnelinvestigations[J].International Journal of Engineering（UE），2010，3（6）：538-553.

[9]范立欽.飛機飛行力學研究工作的回顧與展望[J].飛行力學，1987（01）：18-25.

[10]方振平，陳萬春，張曙光，等.航空飛行器飛行動力學[M].1版.北京：北京航空航天大學出版社，2005.