RBF整體網(wǎng)格變形技術與多體軌跡仿真

曾 錚，王 剛，葉正寅

（西北工業(yè)大學航空學院，陜西西安 710072）

RBF整體網(wǎng)格變形技術與多體軌跡仿真

曾 錚＊，王 剛，葉正寅

（西北工業(yè)大學航空學院，陜西西安 710072）

在多體分離的CFD動態(tài)軌跡仿真中，計算網(wǎng)格需要準確刻畫固壁邊界之間幅度很大的相對運動。理想的狀況是采用整體網(wǎng)格，在不增刪網(wǎng)格節(jié)點和不改變原網(wǎng)格拓撲結構的條件下，通過網(wǎng)格變形來刻畫多體邊界間的相對運動。為此，以先前發(fā)展的徑向基函數(shù)（RBF）網(wǎng)格變形技術為基礎，通過引入一種改進的拉普拉斯光順算法，大幅改善了單獨RBF網(wǎng)格變形技術在處理大幅度邊界位移時的網(wǎng)格質量下降和網(wǎng)格拓撲破壞問題，成功實現(xiàn)了多體分離中的超大幅度網(wǎng)格變形。在改進的拉普拉斯算法中，針對四面體、三棱柱、六面體、金字塔四種網(wǎng)格單元類型提出了一套統(tǒng)一的網(wǎng)格質量評估標準，并將其用于非結構混合網(wǎng)格下的拉普拉斯光順。通過耦合求解非定常Navier-Stokes方程和剛體六自由度運動方程，利用典型跨聲速下機翼與外掛物分離問題對以上網(wǎng)格變形和光順技術進行驗證，將計算獲得的外掛物投放過程的空間位置、姿態(tài)角和氣動力與風洞實驗結果進行了對比分析，結果顯示：在不采用網(wǎng)格嵌套和局部重生的前提下，與改進的拉普拉斯光順相結合的RBF網(wǎng)格變形技術可以有效地解決多體分離過程中的網(wǎng)格大幅度變形問題，且網(wǎng)格質量能夠滿足精確刻畫外掛物投放飛行軌跡的要求。

多體分離；網(wǎng)格變形；網(wǎng)格質量；拉普拉斯光順；徑向基函數(shù)；軌跡仿真

0 引 言

多體分離是實際飛行器設計中需要經(jīng)常面對的技術問題。常見多體分離問題包括外掛物與載機分離過程、座艙蓋/座椅的彈射過程、多級火箭的級間分離、多彈頭再入時的展開、子母彈拋撒過程等。此類問題包含的分離體數(shù)目多、相對位移大、氣動力非定常效應顯著、氣動干擾嚴重［1］。目前，多體分離軌跡問題主要依托風洞實驗和數(shù)值模擬兩類考察手段。風洞實驗采用的技術方法主要有捕獲軌跡系統(tǒng)（CTS），自由投放方法（又稱動態(tài)投放實驗），網(wǎng)格掃描方法和流向角方法等?？傮w來說，地面設備不足使得采用風洞實驗手段研究多體分離問題相對比較困難，尤其是在高速機動飛行等非常規(guī)情況下，實驗研究的局限性愈發(fā)明顯。自由投放實驗費用高，危險性強，而且不易測量。CTS往往費用高昂，實驗狀態(tài)極多［2］。近年來，數(shù)值方法的進步使得仿真計算越來越能逼近現(xiàn)實狀況，采用數(shù)值方法模擬多體分離軌跡也逐漸成為一種常用的技術手段。在CFD工作者的不懈努力下，多體分離問題的數(shù)值模擬方法已經(jīng)從最初的準定常計算發(fā)展為完全非定常模擬［3］。采用非定常數(shù)值方法進行多體分離軌跡仿真時必須解決的核心問題是：如何采用動態(tài)計算網(wǎng)格刻畫各固壁邊界之間幅度很大的相對運動。針對這一問題，國內(nèi)外采用的動態(tài)網(wǎng)格技術大致分為以下幾類：

第一類是嵌套網(wǎng)格技術，這類技術已經(jīng)發(fā)展的比較成熟，是一種國內(nèi)外很常用的求解多體分離問題的方法。但嵌套網(wǎng)格由于運動子塊之間及其與背景網(wǎng)格之間的相對位置隨時間變化，所以交疊區(qū)內(nèi)的人工邊界即所謂“洞”邊界必須每個時間步都重新建立，重新定義計算網(wǎng)格。在計算中需要在重疊區(qū)頻繁地進行插值運算，實現(xiàn)運動子塊的物理量傳遞，這會帶來解的精度損失［4］。

第二類是網(wǎng)格變形方法，這類方法在非結構網(wǎng)格中運用較多，代表性算法包括彈簧法、彈性體方法、Delaunay圖法和徑向基函數(shù)網(wǎng)格變形方法等。彈簧法需要構建和求解大型矩陣，計算量大，而且主要是針對四面體非結構網(wǎng)格。彈性體方法對任意網(wǎng)格類型都適用，并能很好地保證變形后的網(wǎng)格質量，但其計算量非常大，難以在非定常流動計算中多次反復使用。Delaunay圖法具有較高的效率和魯棒性，但難以處理較大變形。近年來，采用徑向基函數(shù)（RBF）插值的網(wǎng)格變形方法逐漸興起。該方法先運用徑向基函數(shù)對物面邊界網(wǎng)格節(jié)點的位移進行插值，再利用構造出來的徑向基函數(shù)序列將物面的位移效應光滑地分散到整個網(wǎng)格區(qū)域的節(jié)點上［5］。結合數(shù)據(jù)精簡算法的RBF網(wǎng)格變形技術在處理小幅變形問題時具有很高的效率，且網(wǎng)格質量可以得到保證。但是，單純的RBF網(wǎng)格變形技術存在極值聚集現(xiàn)象，對于多體分離問題中的超大幅度網(wǎng)格變形，在極值的積累效應下，計算網(wǎng)格易被破壞。

第三類是計算網(wǎng)格重新生成技術。對于結構化網(wǎng)格，由于其拓撲結構相對固定，因此可以在計算過程中根據(jù)邊界位置運用代數(shù)方法重新生成新的網(wǎng)格。而非結構網(wǎng)格下的網(wǎng)格局部重生辦法則要復雜得多。由于每一次網(wǎng)格重生都可能會改變原先的網(wǎng)格節(jié)點和網(wǎng)格單元的對應關系，這給流場信息的傳遞帶來了不便。對于非結構網(wǎng)格下的大幅度變形問題，人們往往采用網(wǎng)格變形和網(wǎng)格局部重生相結合的策略。在單步小位移的情況下利用彈簧法或RBF方法進行網(wǎng)格變形，當總位移逐漸增大至網(wǎng)格質量低于設定標準時，進行局部網(wǎng)格重新生成［6-7］。這種方法能夠保證整個分離過程中的網(wǎng)格質量，但網(wǎng)格重生過程影響效率，且在物理信息傳遞中會產(chǎn)生插值精度損失。

為了克服多體分離大變形中網(wǎng)格質量易下降、流場信息插值精度易損失的問題，本文發(fā)展了一套RBF和改進的拉普拉斯光順相結合的網(wǎng)格變形技術。采用整體網(wǎng)格，在單步小幅變形時選用RBF方法實現(xiàn)網(wǎng)格變形，當網(wǎng)格質量下降到設定標準以下時引入改進的拉普拉斯光順，在不破壞網(wǎng)格節(jié)點之間相鄰關系的前提下進行網(wǎng)格節(jié)點的局部光順調整，使網(wǎng)格質量得到改善。本文將這一網(wǎng)格變形方法應用于多體分離問題的數(shù)值模擬，耦合求解非定常NS方程和六自由度運動方程，形成了一套基于動態(tài)混合網(wǎng)格的多體分離問題數(shù)值模擬方法。利用該方法數(shù)值模擬了典型的多體分離標模問題（機翼/外掛物分離），并與AEDC風洞實驗結果［8］進行了對比。

1 理論背景

1.1 非定常N-S方程求解

本文使用的流場求解器是作者課題組自主開發(fā)的一套基于非結構混合網(wǎng)格的雷諾平均N-S方程求解程序，簡稱HUNS3D［9］。為了考慮固體邊界的六自由度運動，流動控制方程采用ALE（Arbitrary Lagrangeian Eulerian）方法描述的非定常雷諾平均NS方程。在計算網(wǎng)格的任意運動和變形情形下，該方程在直角坐標系下的積分形式為：

運用求解定常流動的時間推進算法對方程（4）進行虛時間迭代求解，當?shù)_到收斂標準，即≈0時，相應地有。顯然，此時的Q（k＋1）可以被視為方程（3）的近似解，也就是真實時間層第n＋1時間步的流動參數(shù)值Qn＋1。具體方法可參見文獻［10-11］。

1.2 六自由度運動方程求解

多體分離問題中，被考察的物體一般有六個自由度。相應的有六個動力學方程，其中三個描述平動，三個描述轉動。另外，還有六個相對應的運動學方程，分別確定物體的空間位置和姿態(tài)［12］?？紤]圖1中給定的地面坐標系［13］，X軸平行于地平面指向某任意選定方向，Z軸鉛垂向下，Y軸垂直XZ平面，按右手定則確定方向。地面坐標系與固連于飛行器并隨飛行器運動的隨體坐標系之間的夾角構成了飛行器的姿態(tài)角。隨體坐標系中X軸在飛行器對稱平面內(nèi)，平行于機身軸線指向前；Z軸亦在對稱平面內(nèi)，垂直于X軸；Y軸垂直對稱平面指向右。在圖1中，為滾轉角，θ為俯仰角，ψ為偏航角。

圖1 外掛物歐拉角示意圖Fig.1 Eular angles of store

（1）圖1所示的參考系中，六個動力學方程描述為：

（2）六個運動學方程包括：

公式（5）～（8）中，m代表質量，V為速度，F(xiàn)為力，ω為角速度，M為外力矩，I為慣性矩；上標i和b分別表示慣性坐標系和體軸系；下標a、e、g分別代表氣動力，外力（發(fā)動機推力），體積力。RB-I為體軸系到慣性系轉換矩陣，詳細表達式參見文獻［12］。

剛體運動方程的求解屬于常微分方程初值問題，令l＝［Vx，Vy，Vz，ωx，ωy，ωz，x，y，z，θ，ψ，］，則剛體運動方程可以寫為：

本文采用改進的Adams預估校正法［14］求解剛體運動方程。

1.3 徑向基函數(shù)（RBF）網(wǎng)格變形技術

徑向基函數(shù)插值的一般形式是：

這里的F（r）是插值函數(shù)，N代表插值問題所使用的徑向基函數(shù)的總數(shù)目，（‖r-ri‖）是所采用的徑向基函數(shù)的通用形式，‖r-ri‖是位置矢量r到ri的距離，ri代表第i號徑向基函數(shù)的支撐點的位置，wi是與第i號徑向基函數(shù)相對應的權重系數(shù)。徑向基函數(shù)的類型很多，本文選取的是比較適合網(wǎng)格變形插值使用的Wendland C2函數(shù)，其具體形式如下：

這里下標S表示物面屬性，分別表示N個物面邊界節(jié)點在x、y、z方向上的位移分量；是待定的權重系數(shù)序列；矩陣Φ的具體形式為：

通過方程（12-14）求解出WX、WY和WZ后就可以獲得物面位移的徑向基函數(shù)插值。這樣計算域內(nèi)任意網(wǎng)格節(jié)點j的位移就可以通過公式（10）確定。

上述徑向基函數(shù)網(wǎng)格變形方法的關鍵環(huán)節(jié)是將物面變形位移通過徑向基函數(shù)插值來近似代替的過程，也就是建立和求解方程（12）～（14）的過程。為了敘述方便，以下將方程（12-14）描述的插值條件統(tǒng)一表述為下面的形式：

為了提高網(wǎng)格變形效率，采用文獻［5］的改進的貪心算法在物面節(jié)點位移的擬合過程中進行數(shù)據(jù)精簡，在滿足網(wǎng)格變形精度要求的前提下，可以盡量降低條件矩陣Φ的維數(shù)。

1.4 一種改進的拉普拉斯光順和網(wǎng)格質量評估函數(shù)

徑向基函數(shù)網(wǎng)格變形算法雖然具有良好的魯棒性和很高的計算效率，但是在單獨使用它處理多體分離帶來的超大幅度網(wǎng)格變形時局部網(wǎng)格質量會下降，嚴重時甚至導致網(wǎng)格拓撲破壞。所以當局部網(wǎng)格單元質量下降到一定程度時，需要運用網(wǎng)格光順，改善變形后的網(wǎng)格質量。經(jīng)典拉普拉斯光順的原理是：調整每個網(wǎng)格頂點至其相鄰頂點的幾何中心。本文對其進行改進［15］，移動頂點時考慮網(wǎng)格質量，僅當該頂點連接的網(wǎng)格單元質量獲得有效提高時，才對該頂點進行移動。

當確定光順區(qū)域后，節(jié)點新的位置按如下方法計算：

其中，β為初始松弛因子，一般取0.5左右。

本文中拉普拉斯光順啟動的判據(jù)是：調整節(jié)點后，與該節(jié)點相鄰的所有網(wǎng)格單元的最差質量系數(shù)被提高且平均質量系數(shù)保持在原來80%以上。對于非結構混合網(wǎng)格，通常有四面體、金字塔、三棱柱、六面體四種網(wǎng)格單元類型。針對上述這四種網(wǎng)格單元，本文提出一套統(tǒng)一的網(wǎng)格質量評估函數(shù)。下面簡要介紹這種網(wǎng)格質量評估函數(shù)。

對于四面體單元，令T作為任意一個包含四個節(jié)點的四面體單元，節(jié)點用Pn表示，n＝0，1，2，3。定義邊矢量ek，n＝xk-xn（k≠n，k＝0，1，2，3），節(jié)點n處的四面體單元雅克比矩陣為：

令αn代表An的行列式。四面體體積，所以，An的構成應使αn＞0。

定義正四面體Te，其四個頂點坐標為（0，0，0），。

讓Wn代表Te第n個頂點的雅克比矩陣。如：

令T＋代表任意具有正體積的四面體，其雅克比矩陣An可逆，定義雅克比矩陣An的加權條件數(shù)為：

其中W＝W0。令Sn＝AnW-1，Sn實際反映的是由正四面體向所求四面體轉換的線性關系。kω（An）即線性轉換矩陣Sn的條件數(shù)。據(jù)此，定義網(wǎng)格質量評估函數(shù)為：

具體推導過程可參見文獻［16］。該函數(shù)具有優(yōu)良網(wǎng)格質量評估函數(shù)所應具備的五大特點，分別是：

①它是關于節(jié)點坐標的連續(xù)函數(shù)。

②網(wǎng)格單元剛性旋轉其函數(shù)值不變。

③函數(shù)取值范圍在0到1之間，取值越大代表網(wǎng)格質量越好。

④正四面體剛性放大或旋轉時其函數(shù)值始終為1。

⑤對于一個質量較好的網(wǎng)格，等比例縮小使其體積趨近于0，其值不變。只有當網(wǎng)格體積趨近于0且存在某相鄰棱邊的長度不相應趨近于0時其函數(shù)值才趨近于0。

本文將該函數(shù)擴展應用到金字塔、三棱柱和六面體三種網(wǎng)格單元類型中。把單元每個頂點及與其相連的三邊視作一個四面體，分別寫出各自的雅克比矩陣，求解各自的質量系數(shù)，即3/kω函數(shù)值，所有頂點的函數(shù)值的平均值即為該類型網(wǎng)格單元的質量系數(shù)，此時的理想四面體雅克比矩陣W應根據(jù)具體網(wǎng)格類型下具體頂點處理想的四面體形式相應給出。例如，六面體下任一頂點的理想四面體是：與頂點相鄰的三條棱邊為單位長度且互相垂直的四面體。對于金字塔單元，有一個頂點與四條邊相鄰，不能夠給出該頂點處的三階雅克比矩陣。此時應將其由對角線剖開成兩個四面體，分別計算這兩個四面體的質量系數(shù)，再取平均值，該值即為該頂點在金字塔網(wǎng)格單元中的質量系數(shù)。如圖2，對于頂點5，其質量系數(shù)為頂點5、1、4、2和頂點5、3、2、4所構成的兩個四面體質量系數(shù)的平均值。

理論和實踐均證明，上述網(wǎng)格質量評估函數(shù)可以對四種網(wǎng)格的網(wǎng)格單元的質量進行統(tǒng)一的刻畫，且具有良好的取值特性。

圖2 金字塔單元點序示意圖Fig.2 Node list of pyramid unit

2 算例驗證

2.1 模型幾何參數(shù)

選擇AEDC外掛物投放標模作為算例考察本文方法在多體分離問題中的魯棒性和精度。AEDC模型有較詳細的幾何數(shù)據(jù)和實驗數(shù)據(jù)［8］，很多研究者都選用它做多體分離問題的研究算例［17-18］。AEDC模型的機翼為一個具有45°傾斜角的三角翼，機翼展向截面為NACA-64A010翼型。掛架的構造為一個對稱的弧面-平面-弧面組合體，對稱線沿展向距翼根弦線的距離為3.3m，掛架頂點距機翼前緣的距離為0.61m。外掛物為帶有四個對稱尾翼的旋成體，其中心圓柱體直徑為0.5m，總長度為3.017m。外掛物和掛架之間的距離只有3.66cm。外掛物質心前后作用兩個彈射力分別為10.7kN和42.7kN，作用時間0.054s。其他具體參數(shù)見圖3和文獻［17］。

圖3 機翼外掛物模型幾何參數(shù)示意圖（mm）Fig.3 Geometry parameters of wing-store model（mm）

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

由于外掛物和掛架之間的距離太小，一般的網(wǎng)格生成軟件沒有辦法在附面層內(nèi)生成粘性網(wǎng)格。本文采用的方法是將機翼和外掛物拉開一段距離生成附面層內(nèi)的粘性網(wǎng)格后，再采用RBF網(wǎng)格變形技術將機翼外掛物還原至它們的初始位置。模型表面網(wǎng)格點數(shù)為53 135，表面單元數(shù)為75 992。整體網(wǎng)格共有1 531 886個節(jié)點和3 275 111個單元。表面網(wǎng)格及掛架外掛物之間的空間網(wǎng)格如圖4所示。來流馬赫數(shù)為0.95，單位雷諾數(shù)為7.874×106，攻角0°，飛行高度是8 000m。流場求解采用Roe格式進行有限體積的空間離散，雙時間方法進行時間推進，隱式LU-SGS格式進行虛時間層迭代，湍流模型采用Menter的k-ωSST兩方程模型。RBF網(wǎng)格變形的收斂標準為1×10-5，采用精簡算法時，徑向基函數(shù)的支撐點的個數(shù)限制在1500內(nèi)。若在網(wǎng)格變形過程中選擇1500個點還未降到收斂標準以下，不再繼續(xù)加入更多的支撐點。

圖4 表面網(wǎng)格和初始網(wǎng)格示意圖Fig.4 Surface mesh and initial volume mesh

2.3 模擬結果

初始時刻定常流場計算的壓力分布如圖5所示。由于流動是跨聲速的，因此在外掛物中部和尾部區(qū)域出現(xiàn)激波。由于吊艙和掛架很近，二者之間有較強的氣動干擾，使這一塊區(qū)域的流動比較復雜。

圖5 初始時刻壓強分布云圖Fig.5 Pressure contours at initial time

為了考察網(wǎng)格光順效果，本文選取單獨采用RBF網(wǎng)格變形技術進行100個非定常時間步（約0.2s）后的瞬態(tài)網(wǎng)格，對其執(zhí)行了150次光順，從圖6可以看出，網(wǎng)格質量改善效果非常顯著。

由圖6可見，光順次數(shù)越多，網(wǎng)格節(jié)點的整體分布越趨于均勻，意味著網(wǎng)格的質量越好，但耗費的計算時間也相應地增多。算例采用的參數(shù)設置是：從第50個非定常時間步開始啟動光順，每個非定常時間步執(zhí)行10次光順，初始松弛因子為0.5。單步動態(tài)混合網(wǎng)格的生成時間為459.856s，其中RBF網(wǎng)格變形耗時178.626s，10次光順耗時281.23s。相比起一次非定常時間步的計算時間，一次動態(tài)混合網(wǎng)格的生成時間只占了很小的比例。

圖6 光順效果圖Fig.6 Effect of smooth

圖7給出外掛物投放過程中六個不同時刻的網(wǎng)格截面示意圖。當外掛物與機翼之間的分離距離達到較大程度時，網(wǎng)格依然保持著良好的拓撲外形。

圖7 分離過程網(wǎng)格變化示意圖Fig.7 Mesh slices of separation

計算得到的投放過程中外掛物質心位移和姿態(tài)角與風洞實驗的比較見圖8和圖9，總體上計算結果和實驗值吻合很好。在圖8中，Xg、Yg、Zg表示質心位移，Z方向的質心位移匹配度最高，這是因為彈射力和重力在這個方向上遠大于氣動力，占了主要支配的地位。由于阻力的作用，外掛物會持續(xù)向后運動，CFD/RBD的仿真結果準確地反映了這一現(xiàn)象。在橫側向上，外掛物的質心會先略微向左，之后會緩慢地向右移動，這一點也與實驗相符。

圖8 質心位移與實驗值的對比圖Fig.8 Calculated centeroid displacement vs.experiment

圖9 姿態(tài)角與實驗值的對比圖Fig.9 Calculated orientation angle vs.experiment

圖9給出了三個姿態(tài)角隨時間變化的結果。由圖可知，俯仰角和偏航角和實驗結果匹配良好。由于彈射力矩的作用，外掛物一開始會有一個抬頭的效應，0.054s之后彈射力消失，俯仰力矩會逐漸使外掛物低頭。外掛物在前0.5s左右先向右偏航，而后會逐漸左偏。本文的CFD/RBD軌跡仿真真實地反應了這些現(xiàn)象。外掛物在前0.5s一直是順時針滾轉的趨勢，這一結果和實驗一致。滾轉角在0.3s后與實驗值在數(shù)值上出現(xiàn)了一定的偏差，這是由于滾轉慣性矩遠遠小于俯仰慣性矩和偏航慣性矩，這導致滾轉角對于氣動力的變化特別敏感，也使關于它的仿真變得尤其困難。

圖10和圖11是軌跡仿真過程中氣動力和氣動力矩隨時間變化的示意圖。從圖中可以看出整個仿真過程中氣動力和氣動力矩與實驗值的整體變化趨勢吻合。圖10中X方向氣動力略大，主要原因是實驗中外掛物后部用尾撐桿支撐，略微抵消了氣動阻尼的作用，使實驗測得的氣動力略小。

圖10 氣動力與實驗值的對比圖Fig.10 Calculated aeroforces vs.experiment

圖11 氣動力矩與實驗值的對比圖Fig.11 Calculated aeromoments vs.experiment

圖12給出外掛物在0.8s內(nèi)與機翼和外掛物分離過程的示意圖，可以更加直觀形象地看到這個過程。此時的外掛物到機翼大約有一個展長的距離，邊界之間的干擾已經(jīng)基本可以被忽略，應可以對外掛物進行單體軌跡仿真。

圖12 機翼與外掛物分離過程示意圖Fig.12 Separation process of wing-store model

3 結 論

本文以非結構混合網(wǎng)格為背景，針對四面體、三棱柱、六面體、金字塔四種網(wǎng)格單元提出了一套統(tǒng)一的網(wǎng)格質量評估標準，該質量評估標準可以準確統(tǒng)一地反映網(wǎng)格的質量。將上述質量標準作為光順的執(zhí)行判據(jù)，提出了一套改進的拉普拉斯網(wǎng)格光順方法。

將上述網(wǎng)格光順方法與徑向基函數(shù)網(wǎng)格變形結合，形成了一種改進的RBF整體網(wǎng)格變形方法。在不進行網(wǎng)格重構和網(wǎng)格嵌套的前提下，該方法改善了單純RBF網(wǎng)格變形技術所帶來的網(wǎng)格質量下降和網(wǎng)格拓撲破壞問題，成功實現(xiàn)了多體分離中的超大幅度網(wǎng)格變形。

通過耦合求解非定常雷諾平均NS方程與六自由度方程，并采用以上改進的網(wǎng)格變形技術，實現(xiàn)了機翼外掛物分離軌跡仿真，證明本文網(wǎng)格變形方法具有很好的實用性。仿真結果與實驗值對比吻合良好，驗證了本文多體分離軌跡仿真方法的魯棒性和精度。

［1］ 郭正.包含運動邊界的多體非定常流場數(shù)值模擬方法研究［D］.［博士學位論文］.長沙：國防科學技術大學，2002.

［2］ Wang Wei.Numerical simulation technique research and experiment verification for unsteady multi-body flowfield involving relative movement［D］.［PhD.Thesis］.Changsha：National University of Defense Technology，2002.（in Chinese）王巍.有相對運動的多體分離過程非定常數(shù)值算法研究及實驗驗證［D］.長沙：國防科學技術大學，2008.

［3］ Duan Xupeng，Chang Xinghua，Zhang Laiping.A CFD-and-6DOF-coupled solver for multiple moving object problems based on dynamic bybrid grids［J］.Acta Aerodynamica Sinica，2011，29（4）：447-452.（in Chinese）段旭鵬，常興華，張來平.基于動態(tài)混合網(wǎng)格的多體分離數(shù)值模擬方法［J］.空氣動力學學報，2011，29（4）：447-452.

［4］ Tian Shuling，Wu Yizhan，Xia Jian.3Dstore separation simulation using unstructured overlapping grid［J］.Acta Aerodynamica Sinica，2007，25（2）：245-249.（in Chinese）田書玲，伍貽兆，夏健.基于非結構重疊網(wǎng)格的二維外掛物投放模擬［J］.空氣動力學學報，2007，25（2）：245-249.

［5］ Gang Wang，Haris Hameed Mian，Zheng-Yin Ye，Jen-Der Lee.Improved point selection method for hybrid-unstructured mesh deformation using radial basis functions［J］.AIAA Journal，2015，53（4）：1016-1025.

［6］ Baum J D，Lohner R.Three-dimensional store separation using a finite element solver and adaptive remeshing［C］//29th Aerospace Sciences Meeting，1991-0602.

［7］ Zhang laiping，Deng Xiaogang，Zhang Hanxin.Reviews of moving grid generation techniques and numberical methods for unsteady flow［J］.Advances in Mechanics，2010，40（4）：424-447.（in Chinese）張來平，鄧小剛，張涵信.動網(wǎng)格生成技術及非定常計算方法進展綜述［J］.力學進展，2010，40（4）：424-447.

［8］ Heim E R.CFD wing/pylon/finned store mutual interference wind tunnel experiment［R］.Arnold Engineering Development Center，Arnold Afs Tn，1991.

［9］ Wang G，Ye Z.Generation of three dimensional mixed and unstructured grids and its application in solving Navier-Stokes equations［J］.Acta Aeronautica et Astronautica Sinica，2003，24（5）：385-390.（in Chinese）王剛，葉正寅.三維非結構混合網(wǎng)格生成與N-S方程求解［J］.航空學報，2003，24（5）：385-390.

［10］Jameson A.Time dependent calculations using multigrid，with applications to unsteady flows past airfoils and wings［C］//10th Computational Fluid Dynamics Conference，1991-1596.

［11］Gaitonde A L.A dual time method for the solution of the unsteady Euler equations［J］.Aeronautical Journal，1994，98（978）：283-291.

［12］Gao Hao，Zhu Peishen，Gao Zhenghong.The advanced flight dynamics［M］.Beijing：National Defence Industry Press，2003.（in Chinese）高浩，朱培申，高正紅.高等飛行動力學［M］.北京：國防工業(yè)出版社，2003.

［13］Sahu J.Time-accurate numerical prediction of free-flight aerodynamics of a finned projectile［J］.Journal of Spacecraft and Rockets，2008，45（5）：946-954.

［14］Jiang Yuewen，Ye Zhengyin，Zhang Weiwei.Semi-implicit solution of aeroelastic equations in time domain［J］.Engineering Mechanics，2012，29（4）：66-71.（in Chinese）蔣躍文，葉正寅，張偉偉.一種半隱式的氣動彈性時域求解方法［J］.工程力學，2012，29（4）：66-71.

［15］Field D A.Laplacian smoothing and Delaunay triangulations［J］.Communications in Applied Numerical Methods，1988，4（6）：709-712.

［16］Freitag L A，Knupp P M.Tetrahedral element shape optimization via the Jacobian determinant and condition number［C］//8th International Meshing Round Table，Lake Tahoe，1999：247-258.

［17］Lijewski L E，Suhs N E.Time-accurate computational fluid dynamics approach to transonic store separation trajectory prediction［J］.Journal of Aircraft，1994，31（4）：886-891.

［18］Panagiotopoulos E E，Kyparissis S D.CFD transonic store separation trajectory predictions with comparison to wind tunnel Investigations［J］.International Journal of Engineering，2010，3（6）：538-553.

Enhanced RBF mesh deformation method and multi-body trajectory simulation

Zeng Zheng，Wang Gang，Ye Zhengyin
（SchoolofAeronautics，NorthwesternPolytechnicalUniversity，Xi’an710072，China）

Multi-body separation is a kind of common mechanical problem which exists widely in the field of Aeronautics and Astronautics.Because of its complexity，studying this kind of problems is very difficult either using flight or wind tunnel test.Therefore，the establishment of an unified，efficient and accurate numerical method for simulating multi-body separation has great application value.In multi-body trajectory simulation，computational grids need to accurately depict the large magnitude of relative motion among all solid boundaries.The ideal approach is using a single set of mesh，without adding or deleting the grid nodes or changing the original grid topology，just through mesh deformation to depict the relative movement.On the basis of radial basis function（RBF）mesh deformation technology，the intelligent Laplacian smoothing algorithm is introduced，which greatly improve the ability of pure RBF mesh deformation technology to avoid quality decreasing and topology damaging，and successfully simulate the process of multi-body separation.The flow solver includes four types of grid cell，which respectively are tetrahedron，triangular prism，hexahedron and pyramid.A new unified mesh quality assessment standard is presented as the execute criterion of intelligent Laplacian smoothing.This mesh deformation and smoothing technology is fully integrated into Navier-Stokes equations and 6-DOF moving equations to validate its feasibility by a typical store separation problem.The computational results are compared with experimental data，which demonstrate the accuracy and ability of presented method indealing with the large magnitude of relative motion among multiple moving boundaries.

multi-body separation；mesh deform；mesh quality；Laplacian smoothing；radial basic function；trajectory simulation

V211.3

：Adoi：10.7638/kqdlxxb-2014.0143

0258-1825（2015）02-0170-08

2015-01-04；

：2015-03-02

國家自然科學基金（91216202）

曾錚＊（1989-），男，重慶人，研究生，研究方向：計算流體力學.E-mail：zzviolin＠m(xù)ail.nwpu.edu.cn

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

10.7638/kqdlxxb-2014.0143 Zeng Z，Wang G，Ye Z Y.Enhanced RBF mesh deformation method and multi-body trajectory simulation［J］.Acta Aerodynamica Sinica，2015，33（2）：170-177.