基于動(dòng)態(tài)網(wǎng)格的多體分離計(jì)算技術(shù)研究

徐琳 宋萬強(qiáng)

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

關(guān)健詞：動(dòng)態(tài)網(wǎng)格;多體分離;數(shù)值模擬;網(wǎng)格變形;Laplace

中圖分類號(hào)：V224.5 文獻(xiàn)標(biāo)識(shí)碼：A

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

多個(gè)物體分離時(shí)各模型邊界之間的相對(duì)運(yùn)動(dòng)幅度較大，這對(duì)于采用數(shù)值模擬方法進(jìn)行分離軌跡求解的計(jì)算網(wǎng)格提出了更高的要求，也是其所需解決的關(guān)鍵問題。目前國(guó)內(nèi)外比較成熟的動(dòng)態(tài)網(wǎng)格求解技術(shù)大致可以分為三種。

（1）嵌套網(wǎng)格技術(shù)

嵌套網(wǎng)格技術(shù)由運(yùn)動(dòng)物體網(wǎng)格和背景網(wǎng)格組成，通過對(duì)重疊區(qū)插值運(yùn)算完成網(wǎng)格間的數(shù)據(jù)傳遞計(jì)算。此類方法為國(guó)外常用的運(yùn)動(dòng)問題求解算法，相對(duì)成熟。其缺點(diǎn)是頻繁的插值運(yùn)算會(huì)帶來解的精度損失。

（2）網(wǎng)格變形技術(shù)

隨著模型相對(duì)位置的變化，計(jì)算網(wǎng)格會(huì)按照設(shè)定的變形算法進(jìn)行壓縮或拉伸變形，以保證網(wǎng)格與模型邊界的貼合。常用算法有彈簧法、彈性體方法Delaunay圖法和RBF網(wǎng)格變形方法等。此類方法多用于非結(jié)構(gòu)網(wǎng)格，其缺點(diǎn)是當(dāng)相對(duì)運(yùn)動(dòng)較大時(shí)，網(wǎng)格質(zhì)量無法保證。

（3）網(wǎng)格重生技術(shù)

當(dāng)模型相對(duì)位置發(fā)生變化時(shí)，網(wǎng)格依據(jù)新的模型重新生成，并通過插值將原網(wǎng)格的物理量傳遞給新的網(wǎng)格。此方法能保證求解過程中的網(wǎng)格質(zhì)量，其缺點(diǎn)是頻繁的插值運(yùn)算會(huì)帶來解的精度損失，且網(wǎng)格重生會(huì)耗時(shí)較長(zhǎng)，影響效率[5]。

本文采用網(wǎng)格變形與局部重構(gòu)相結(jié)合的動(dòng)態(tài)網(wǎng)格方法：在單步小位移的情況下用Laplace[6，7]算法進(jìn)行網(wǎng)格變形，當(dāng)總的位移量逐漸增大至網(wǎng)格質(zhì)量低于設(shè)定標(biāo)準(zhǔn)時(shí)，進(jìn)行局部網(wǎng)格重新生成。采用這種動(dòng)態(tài)網(wǎng)格方法可以計(jì)算一個(gè)或多個(gè)運(yùn)動(dòng)物體的軌跡。本文采用非結(jié)構(gòu)網(wǎng)格、Euler方程、準(zhǔn)定常的方法開展外掛投放數(shù)值求解，模擬了典型的多體分離標(biāo)模（機(jī)翼/外掛物分離），并與阿諾德工程發(fā)展復(fù)雜綜合體（AFDC）風(fēng)洞試驗(yàn)結(jié)果183進(jìn)行了對(duì)比。

1 軌跡計(jì)算方法

本文運(yùn)用運(yùn)動(dòng)網(wǎng)格變形/局部重構(gòu)計(jì)算的動(dòng)態(tài)網(wǎng)格求解方法，對(duì)機(jī)翼外掛物分離軌跡進(jìn)行數(shù)值模擬。具體的計(jì)算過程如圖1所示。

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

2飛行力學(xué)模型

本文中使用到的飛行力學(xué)模型是建立在剛體六自由度運(yùn)動(dòng)基礎(chǔ)上[9，10]，涉及到的運(yùn)動(dòng)方程如下所示：

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

Ω'=ω'（4）式中：F為運(yùn)動(dòng)物體上的合力;m為指運(yùn)動(dòng)物體的質(zhì)量;r_C（t）和r_C（t）為指速度及質(zhì)心的位置;L_C為運(yùn)動(dòng)物體上相對(duì)于物體質(zhì)心的合力矩;ω為角速度;I_C相對(duì)于運(yùn)動(dòng)物體質(zhì)心的慣性張量;t為當(dāng)前時(shí)間步;t⁰為上一個(gè)時(shí)間步;Ω為旋轉(zhuǎn)角速度。上面的公式可以寫為：

q=f（q）（6）式中：f是相對(duì)于右手坐標(biāo)系。

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

3 網(wǎng)格變形/局部，構(gòu)技術(shù)基本原理

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

式（9）為網(wǎng)格拉伸率P的定義。當(dāng)網(wǎng)格沒有負(fù)體積，且網(wǎng)格拉伸率ρ滿足：ρ_min<ρ<ρ_max。其中，最小網(wǎng)格拉伸率ρ_max與最大網(wǎng)格拉伸率ρ_max為用戶設(shè)定。當(dāng)總的位移量逐漸增大至使網(wǎng)格質(zhì)量低于設(shè)定標(biāo)準(zhǔn)時(shí)，進(jìn)行局部網(wǎng)格重新生成。網(wǎng)格拉伸率P的定義如圖2所示。

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

4 算例驗(yàn)證

4.1 二維雙翼型算例

為了高效地驗(yàn)證本文軌跡計(jì)算的能力，方便直觀地展示計(jì)算過程中的網(wǎng)格變形情況，本文首先采用二維雙翼型算例。本次計(jì)算的二維雙翼型算例是由兩個(gè)NACA0012翼型構(gòu)成，初始位置如圖4（a）所示，上邊的翼型充當(dāng)飛機(jī)，下面的翼型充當(dāng)導(dǎo)彈。本次計(jì)算來流馬赫數(shù)為0.8，雷諾數(shù)為1.8×10⁷，下翼型重量（質(zhì)量）5000kg，當(dāng)?shù)刂亓铀俣葹?.8m/s²，橫搖阻尼系數(shù)取-4.0/rad，縱搖阻尼系數(shù)取-40.0/rad，偏航阻尼系數(shù)取-40.0/rad。流場(chǎng)求解采用歐拉方程準(zhǔn)定常計(jì)算。

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

4.2 AEDC算例計(jì)算

4.2.1 模型幾何參數(shù)

AEDC外掛物標(biāo)模算例的捕獲軌跡系統(tǒng)（CTS）試驗(yàn)是由AEDC于1990年完成的。其模型，即機(jī)翼/掛架/帶舵外掛物模型（Wing/Pylon/Finned-Store，WPSF）是美國(guó)發(fā)起第一次分離投放計(jì)算流體力學(xué)（CFD）驗(yàn)證時(shí)使用的一個(gè)機(jī)翼和導(dǎo)彈的簡(jiǎn)化模型[6]。如圖6所示。很多研究都選用AEDC模型做為運(yùn)動(dòng)軌跡數(shù)值模擬的算例，是因?yàn)檫@個(gè)模型的幾何數(shù)據(jù)和試驗(yàn)數(shù)據(jù)都比較詳細(xì)。

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

4.2.2 初始網(wǎng)格和流場(chǎng)參數(shù)

本次計(jì)算共準(zhǔn)備兩套網(wǎng)格：背景網(wǎng)格wing.grid和運(yùn)動(dòng)物體網(wǎng)格store.grid。這兩套網(wǎng)格采用Pointwise繪制，其中wing.grid共包含131487個(gè)節(jié)點(diǎn)和787942個(gè)單元，store.grid共包含27891個(gè)節(jié)點(diǎn)和143945個(gè)單元，由于本次軌跡計(jì)算采用歐拉方程求解，故并沒有繪制棱柱層進(jìn)行加密，僅對(duì)導(dǎo)彈軌跡區(qū)域附近進(jìn)行了網(wǎng)格局部加密。具體網(wǎng)格示意圖如圖8所示。

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

4.2.3 模擬結(jié)果

通過網(wǎng)格局部重構(gòu)的方法生成的網(wǎng)格如圖9所示，為了保證導(dǎo)彈分離過程中對(duì)導(dǎo)彈流體計(jì)算的準(zhǔn)確性，需保證計(jì)算過程中導(dǎo)彈周圍網(wǎng)格的疏密度，因此，在設(shè)置時(shí)，將導(dǎo)彈周圍加密網(wǎng)格區(qū)域在inner box的范圍內(nèi)，如圖9所示導(dǎo)彈周圍的網(wǎng)格為導(dǎo)彈初始網(wǎng)格的保留。

初始時(shí)刻定常流場(chǎng)計(jì)算的壓力（壓強(qiáng)）分布如圖10所示，由于流動(dòng)是跨聲速的，因此，在外掛物中部和尾部區(qū)域出現(xiàn)激波。由于吊艙和掛架很近，二者之間有較強(qiáng)的氣動(dòng)干擾，使這一區(qū)域的流動(dòng)比較復(fù)雜。圖11為質(zhì)心位移與試驗(yàn)值的對(duì)比圖。從圖11可以看出，采用準(zhǔn)定常流動(dòng)計(jì)算方法，實(shí)現(xiàn)了AEDC機(jī)翼/導(dǎo)彈分離標(biāo)準(zhǔn)模型軌跡模擬，Is時(shí)間內(nèi)導(dǎo)彈投放軌跡計(jì)算結(jié)果與風(fēng)洞試驗(yàn)數(shù)據(jù)基本吻合。

5 結(jié)論

本文通過采用網(wǎng)格變形與局部重構(gòu)相結(jié)合的網(wǎng)格變形技術(shù)，以NACA0012和AFDC為例，耦合求解Euler方程與六自由度方程，采用準(zhǔn)定常的方法實(shí)現(xiàn)了機(jī)翼外掛物分離軌跡數(shù)值模擬，證明本文網(wǎng)格變形/局部網(wǎng)格重構(gòu)方法有很好的實(shí)用性。AEDC仿真結(jié)果與試驗(yàn)值吻合較好，驗(yàn)證了本文多體分離計(jì)算方法的準(zhǔn)確性。

參考文獻(xiàn)

[1]唐志共，李彬，鄭鳴，等.飛行器外掛投放數(shù)值模擬[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2009，27（5）：592-596.

[2]曾錚，王剛，葉正寅.RBF整體網(wǎng)格變形技術(shù)與多體軌跡仿真[J]，空氣動(dòng)力學(xué)學(xué)報(bào)，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]聶璐，向錦武，飛機(jī)外掛物投放安全性的參數(shù)影響分析陰.飛行力學(xué)，2011，29（2）：25-28.

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

[7]胡勇全.基于自動(dòng)生成控制點(diǎn)的拉普拉斯變形[D].大連：大連理工大學(xué)，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ī)飛行力學(xué)研究工作的回顧與展望[J].飛行力學(xué)，1987（01）：18-25.

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