一種用于TSTO 級間分離CFD 計算的網(wǎng)格動態(tài)優(yōu)化技術(shù)

唐 靜，張 健，張耀冰，周乃春，劉 剛

(中國空氣動力研究與發(fā)展中心，綿陽 621000)

0 引言

在天地往返空天運輸系統(tǒng)中，兩級入軌（twostage-to-orbit,TSTO）系統(tǒng)相比單級入軌具有較低的技術(shù)實現(xiàn)難度，可通過運載級的重復(fù)使用降低運輸成本，還可實現(xiàn)軌道器和洲際導(dǎo)彈同平臺發(fā)射[1]。兩級飛行器的安全分離是TSTO 系統(tǒng)設(shè)計的一項關(guān)鍵技術(shù)，尤其是升力體并聯(lián)分離時，級間流動相互干擾強烈，分離過程的預(yù)測和安全性評估極其重要。

計算流體力學（CFD）已廣泛應(yīng)用于多體分離軌跡和姿態(tài)的數(shù)值預(yù)測[2]。數(shù)值格式與計算網(wǎng)格直接影響著預(yù)測結(jié)果的精準度。相比數(shù)值格式，對計算網(wǎng)格影響的研究很少，且常規(guī)網(wǎng)格技術(shù)難以實現(xiàn)網(wǎng)格與流動的實時匹配。近年來，耦合流動特征的網(wǎng)格動態(tài)優(yōu)化技術(shù)被用于提高數(shù)值模擬的精準度[3-4]，得到了國內(nèi)外專家學者的廣泛關(guān)注[5-6]。網(wǎng)格自適應(yīng)技術(shù)可明顯提升大攻角分離流動的模擬能力[7-8]、精細模擬激波和剪切層流動[9-10]，有助于提高氣動力/熱的模擬精度[11-12]。目前網(wǎng)格自適應(yīng)技術(shù)主要應(yīng)用于定常流動的模擬，在非定常流動上的應(yīng)用較少且計算模型較為簡單。Tang 和Zhang 等將自適應(yīng)技術(shù)應(yīng)用于二維三角楔繞流的非定常模擬[13]，Luo 和Fidkowski 結(jié)合自適應(yīng)技術(shù)開展了翼型俯仰運動流動模擬[14]。

本文為了更好地捕捉和精細模擬TSTO 級間分離過程中的級間運動及反射激波，結(jié)合非結(jié)構(gòu)混合網(wǎng)格加密/稀疏優(yōu)化方法和基于壓力比值的激波識別方法，建立了用于飛行器分離過程的并行網(wǎng)格動態(tài)優(yōu)化技術(shù)。最后，基于NNW-FlowStar 軟件和建立的網(wǎng)格動態(tài)優(yōu)化技術(shù)，開展了TSTO 并聯(lián)分過程數(shù)值仿真。結(jié)果表明，通過網(wǎng)格動態(tài)優(yōu)化技術(shù)，可以實時地精細捕捉運動激波及其反射激波，精確預(yù)測分離軌跡和姿態(tài)的演化歷程。

1 網(wǎng)格動態(tài)優(yōu)化技術(shù)

對于真實飛行器相對運動過程的非定常流動數(shù)值模擬，首先要求數(shù)值方法應(yīng)具備極高的魯棒性，以實時地求取飛行器的氣動特性；其次還需要數(shù)值方法具備大規(guī)模并行計算能力，盡可能減少模擬時間開銷。

1.1 魯棒的網(wǎng)格分布優(yōu)化技術(shù)

對于TSTO 并聯(lián)分離過程的非定常流動，網(wǎng)格優(yōu)化不僅包含網(wǎng)格實時加密，還必須包含網(wǎng)格的稀疏，以適應(yīng)兩級之間激波運動的實時捕獲，提高網(wǎng)格的使用效率。

1.1.1 網(wǎng)格加密技術(shù)

為了盡可能提高加密方法的魯棒性，采用各向同性的方式加密網(wǎng)格單元，每類單元加密模式唯一，大大降低了算法的復(fù)雜度。通常情況下，非結(jié)構(gòu)混合網(wǎng)格包括三角形和四邊形2 種二維單元，包括四面體、三棱柱、金字塔和六面體4 種三維單元，各類單元加密方法如圖1 所示。

圖1 非結(jié)構(gòu)混合網(wǎng)格單元加密模式Fig.1 Refinement patterns for different cell types of unstructured mixed grid

在加密網(wǎng)格單元和未加密網(wǎng)格單元之間，將出現(xiàn)懸空節(jié)點，如圖2 所示。常規(guī)的流場模擬軟件不支持懸空節(jié)點，因此，為了增加網(wǎng)格自適應(yīng)方法對流場模擬軟件的通用性，采用多面體轉(zhuǎn)換方法消除懸空節(jié)點，如圖3 所示。

圖2 網(wǎng)格單元間懸空節(jié)點示意圖Fig.2 Suspending nodes between two cells

圖3 消除懸空節(jié)點的多面體網(wǎng)格單元Fig.3 Eliminating suspending nodes with polyhedron cell

1.1.2 網(wǎng)格稀疏技術(shù)

網(wǎng)格稀疏方法采用“回退”方法，如圖4 所示，即網(wǎng)格稀疏是網(wǎng)格加密的逆向過程。由加密-稀疏過程可以看出，網(wǎng)格單元間形成多級多路的樹形關(guān)系結(jié)構(gòu)，可以采用K-D 樹的數(shù)據(jù)結(jié)構(gòu)進行網(wǎng)格單元的管理。更加詳細的非結(jié)構(gòu)混合網(wǎng)格優(yōu)化方法可參考文獻[15]。

圖4 單元加密/稀疏過程及相應(yīng)樹形關(guān)系結(jié)構(gòu)Fig.4 The refinement and coarsening procedure and corresponding tree-based data structure

1.2 高效的網(wǎng)格優(yōu)化并行技術(shù)

計算流體力學求解流場通常采用區(qū)域分解策略實現(xiàn)并行計算，即在負載平衡的約束下進行網(wǎng)格分區(qū)。為了減少計算機內(nèi)存的開銷，求解時每個并行進程只保存局部區(qū)域的網(wǎng)格，各并行進程在分區(qū)交界面附近進行流場數(shù)據(jù)的通信，并采用虛擬網(wǎng)格來存儲并行通信接收的流場數(shù)據(jù)。

為了準確建立進程間的網(wǎng)格并行對應(yīng)關(guān)系，新增的網(wǎng)格對象（包括網(wǎng)格單元和網(wǎng)格點）需滿足兩個條件：

1）唯一性：真實物理網(wǎng)格對象的編號和并行進程編號在所有進程中唯一；

2）同一性：虛擬網(wǎng)格對象的編號和并行進程編號與存在于其他進程上的真實對象相同。

1.2.1 網(wǎng)格加密并行技術(shù)

由網(wǎng)格稀疏方法可知，稀疏過程不涉及并行通信問題，各并行進程獨立開展即可。網(wǎng)格加密涉及新增網(wǎng)格對象，本文采用“先唯一性，后同一性”兩步法策略實現(xiàn)運算的并行。

第一步：為每個并行進程分配互不重疊的編號區(qū)間，新增對象編號只需在編號區(qū)間內(nèi)取值即可滿足唯一性要求。只有在某個并行進程取值超出可用區(qū)間時，才需要通過并行通信確定新的可用區(qū)間。并行進程p的可用編號區(qū)間由下式確定：

式中：Nmax為當前所有進程的計數(shù)器的最大值，n為并行進程總數(shù)。NT為計數(shù)器區(qū)間跨度，此處取512。

第二步：為虛擬網(wǎng)格對象設(shè)置特征數(shù)據(jù)，采用并行通信傳遞特征數(shù)據(jù)到存儲真實網(wǎng)格對象的并行進程，通過特征數(shù)據(jù)識別出相應(yīng)真實網(wǎng)格對象后返回其真實編號。

圖5 給出了局部的網(wǎng)格對象經(jīng)過兩步法滿足兩個并行相容條件的示意圖。經(jīng)過“第一步”后，新增網(wǎng)格單元的編號實現(xiàn)全并行域的唯一，但虛擬網(wǎng)格對象的編號與對應(yīng)的真實對象編號不同，通過“第二步”后實現(xiàn)了二者的同一性。

圖5 兩步法策略實現(xiàn)新增網(wǎng)格對象的并行相容Fig.5 Parallel compatibility recovery of new grid object with two-steps strategy

1.2.2 負載平衡技術(shù)

由于增添或刪除網(wǎng)格單元的操作通常集中在少數(shù)幾個并行進程中，造成后續(xù)流場求解的負載嚴重不平衡，顯著降低并行效率。此處采用“并行重分區(qū)-網(wǎng)格數(shù)據(jù)遷移”的兩步法實現(xiàn)負載的再平衡。并行重分區(qū)采用ParMetis 程序庫[16]實現(xiàn)，然后通過并行通信將不再屬于當前進程的網(wǎng)格單元（將作為虛擬網(wǎng)格單元的臨近單元）傳送到目標進程。

如圖6 所示，初始網(wǎng)格采用4 個并行分區(qū)（相同顏色為同一分區(qū)）進行計算，由于在特征區(qū)域加密網(wǎng)格，造成了負載不平衡，通過“并行重分區(qū)-網(wǎng)格數(shù)據(jù)遷移”的兩步法實現(xiàn)負載再平衡后，并行分區(qū)的交界面也發(fā)生了移動。更詳細的非結(jié)構(gòu)混合網(wǎng)格優(yōu)化并行技術(shù)可參考文獻[17]。

圖6 網(wǎng)格自適應(yīng)加密后負載再平衡Fig.6 Load rebalancing after mesh refinement

2 網(wǎng)格單元加密/稀疏探測器

流場求解是網(wǎng)格動態(tài)優(yōu)化的基礎(chǔ)，本文流場求解采用NNW-FlowStar 軟件[18]。該軟件基于積分型式的雷諾平均Navier-Stokes 方程和采用格心格式的有限體積方法[19]。文中采用Roe 格式[20]計算對流通量，采用中心差分格式計算黏性通量，采用LU-SGS方法[21]進行方程的迭代求解，采用標準SA 模型[22]進行湍流模擬。

TSTO 并聯(lián)分離過程在兩級之間將會出現(xiàn)激波的多次反射，精確模擬并聯(lián)分離過程的關(guān)鍵是精細模擬激波的運動過程。本文采用激波前后壓力比值關(guān)系捕捉激波，即：

式中：下標1 和2 分別表示激波前后的流場參數(shù)；p為壓力；Man為法向馬赫數(shù)；v為速度矢量；c為聲速；η為誤差消除因子，取0.95。

當相鄰網(wǎng)格壓力比值滿足式（2）時，表示網(wǎng)格間存在激波，即網(wǎng)格需要加密，否則網(wǎng)格需要粗化。

3 兩級入軌飛行器并聯(lián)分離模擬

圖7 給出了兩級入軌飛行器的構(gòu)型幾何模型[23]，相應(yīng)的初始基準網(wǎng)格見圖8，其中一級飛行器為下面尺寸更大的運載級，二級為上面尺寸較小的載荷級。并聯(lián)分離過程模擬采用的網(wǎng)格為非結(jié)構(gòu)混合網(wǎng)格，共包含約510 萬個網(wǎng)格單元，其中一級約270 萬，二級約240 萬。為了對比計算結(jié)果，采用2 800 萬的密網(wǎng)格作為參考。圖9 給出了與更密網(wǎng)格（約3 400 萬個網(wǎng)格，加密網(wǎng)格集中在級間干擾區(qū)域）計算俯仰角的對比，可以看出，參考網(wǎng)格基本達到了網(wǎng)格無關(guān)性要求。級間分離高度為20 km，來流馬赫數(shù)為6.0，飛行器級間分離僅由氣動作用力提供驅(qū)動力。

圖7 TSTO 幾何構(gòu)型Fig.7 The geometry configuration of TSTO

圖8 TSTO 初始網(wǎng)格Fig.8 The base mesh of TSTO

圖9 不同網(wǎng)格規(guī)模計算的俯仰角Fig.9 Pitch angles computed by different number of grid cells

圖10 給出了分離過程中兩級飛行器的法向位移（差量反映飛行器的間距）隨時間的變化，其中實線和虛線分別表示一級飛行器和二級飛行器的法向位移。可以看出，采用初始基準網(wǎng)格計算的法向位移與參考值有較大的區(qū)別，在采用網(wǎng)格動態(tài)優(yōu)化后，計算的法向位移與參考值符合得很好。

圖10 分離過程法向位移比較Fig.10 Comparison of vertical displacement of separation

通常情況下，氣動力計算精度比氣動力矩更高，因而位移的預(yù)測相對容易。分離過程的姿態(tài)角變化預(yù)測難度更大，對分離安全性的影響也更顯著。圖11給出了俯仰姿態(tài)角隨時間的演化過程，可以看出，1 s后初始基準網(wǎng)格計算的二級俯仰姿態(tài)角與參考值有近1°的偏差。通過網(wǎng)格的動態(tài)優(yōu)化，計算的俯仰角與參考值差別很小。

圖11 分離過程俯仰角比較Fig.11 Comparison of pitch angle of separation

圖12 給出了分離前期飛行器表面和空間對稱面流場圖，此時二級頭部激波在一級上表面發(fā)生反射，反射激波到達二級下表面后發(fā)生再次反射；然后反射激波向下與一級相遇發(fā)生反射，最終到達二級尾部發(fā)生反射后逐漸減弱。從網(wǎng)格圖可以看出，在多次反射的激波面上，網(wǎng)格都進行了加密優(yōu)化，以精細捕捉激波。

圖12 分離前期流場圖與相應(yīng)網(wǎng)格圖Fig.12 The flow field and mesh slice at early separation

圖13 給出了分離中期流場圖和相應(yīng)的網(wǎng)格圖，可以看出，此時二級頭部激波僅在一級上表面發(fā)生一次反射后便進入空間流場，激波結(jié)構(gòu)相比分離前期簡單。對比圖12 和圖13 可以看出，隨著激波運動和結(jié)構(gòu)不斷變化，加密網(wǎng)格的區(qū)域與動態(tài)激波面始終一致。激波面上的網(wǎng)格實現(xiàn)了實時的動態(tài)優(yōu)化，以準確捕捉復(fù)雜的波系結(jié)構(gòu)，提高飛行器表面壓力分布計算精度，準確模擬了氣動力矩，實現(xiàn)了姿態(tài)角的準確預(yù)測。

圖13 分離中期流場圖與相應(yīng)網(wǎng)格圖Fig.13 The flow field and mesh slice at middle separation

4 結(jié)論

本文針對TSTO 級間分離過程中涉及的級間復(fù)雜激波運動數(shù)值精細模擬問題，建立了非結(jié)構(gòu)混合網(wǎng)格動態(tài)優(yōu)化技術(shù)，實現(xiàn)了飛行器并聯(lián)級間分離過程中網(wǎng)格實時優(yōu)化和運動規(guī)律預(yù)測，結(jié)果表明：

1）本文采用的動態(tài)網(wǎng)格優(yōu)化技術(shù)能精確捕捉到結(jié)構(gòu)復(fù)雜的運動激波及其反射激波，實現(xiàn)了激波面附近網(wǎng)格分辨率的實時優(yōu)化，精細模擬了激波的強度和壓縮效應(yīng)。

2）通過網(wǎng)格動態(tài)優(yōu)化技術(shù)，降低了網(wǎng)格生成對經(jīng)驗的依賴，提高了級間分離過程的運動軌跡和姿態(tài)角的預(yù)測精準度。

本文動態(tài)網(wǎng)格優(yōu)化技術(shù)是以網(wǎng)格本身為優(yōu)化對象，不僅適用于計算流體力學，還可以推廣到其他基于網(wǎng)格的數(shù)值模擬算法中，以提高關(guān)鍵區(qū)域的網(wǎng)格分辨率，進而提高數(shù)值模擬的精準度。