真空羽流智能化計算

蔡國飆，張百一，賀碧蛟，翁惠焱，劉立輝

1. 北京航空航天大學(xué) 宇航學(xué)院, 北京 100191 2. 航天器設(shè)計優(yōu)化與動態(tài)模擬教育部重點實驗室, 北京 100191

衛(wèi)星、飛船、空間站和地外天體探測器等航天器在太空中采用姿軌控發(fā)動機實現(xiàn)位置保持、姿態(tài)控制和軌道轉(zhuǎn)移等。真空環(huán)境中，發(fā)動機噴流向外部環(huán)境自由膨脹形成羽毛狀流場，稱為真空羽流[1]。真空羽流會對航天器產(chǎn)生氣動力、氣動熱、污染、電磁干擾和視場干擾等效應(yīng)，統(tǒng)稱為羽流效應(yīng)[1]。羽流效應(yīng)會干擾航天器正常工作狀態(tài)，甚至影響航天器壽命和任務(wù)成敗。因此，真空羽流及其效應(yīng)評估和防護是航天領(lǐng)域的重要科學(xué)和工程問題。真空羽流包括連續(xù)介質(zhì)流、過渡區(qū)域流和自由分子流等全流域狀態(tài)，流動機理極其復(fù)雜。常用的真空羽流及其效應(yīng)研究方法包括理論研究[2]、在軌實驗[3]、地面實驗[4-5]和數(shù)值模擬[1,6-7]等。北京航空航天大學(xué)研制了基于液氦/液氮組合雙層一體深冷泵的真空羽流效應(yīng)實驗系統(tǒng)[8]，開發(fā)了具有自主知識產(chǎn)權(quán)的、基于直接模擬蒙特卡洛方法(DSMC)[9-11]的通用并行仿真軟件羽流工作站(Plume Work Station, PWS)[1,6-7]，建立了航天器全狀態(tài)羽流效應(yīng)實驗仿真耦合評估體系。研究成果已成功應(yīng)用于長征、神舟、天宮、嫦娥和天問等系列航天器。

真空羽流控制方程為Boltzmann方程，理論研究非常困難；在軌實驗風(fēng)險大，且費效比巨大；地面實驗常受限于真空艙尺寸和動態(tài)真空度而無法進行全尺寸發(fā)動機羽流及其效應(yīng)研究。因此，數(shù)值模擬是當(dāng)前航天領(lǐng)域應(yīng)用最廣泛的研究方法。DSMC方法是稀薄流領(lǐng)域中較為成熟且精度最高的數(shù)值模擬方法，但DSMC方法本質(zhì)上是基于第一性原理的粒子模擬，需要同時跟蹤大量模擬粒子的運動，非常耗時[12]。如月球探測器著陸月面過程中，數(shù)值模擬全尺度(～30 m)發(fā)動機羽流及其與月壤作用過程，DSMC計算的時間可達數(shù)天甚至幾周，嚴(yán)重影響工程設(shè)計部分的迭代設(shè)計進度。此外，中國正在開展載人登月關(guān)鍵技術(shù)攻關(guān)，真空羽流與月壤相互作用評估及控制是其中一項關(guān)鍵技術(shù)。真空羽流與月面作用會激起月塵[6,12]，引發(fā)儀器讀數(shù)錯誤[13]和視野遮擋[14]等問題，因此實時預(yù)測羽流場及月塵可為保障航天員/航天器安全及任務(wù)成功提供關(guān)鍵數(shù)據(jù)支撐；但受限于計算效率，DSMC方法無法實現(xiàn)實時預(yù)測。由此可見，大幅提高真空羽流數(shù)值模擬效率十分必要。

近年來，人工智能[15]和深度學(xué)習(xí)[16]的發(fā)展為科學(xué)家和工程師們提供了新的技術(shù)手段，相關(guān)技術(shù)已廣泛應(yīng)用于計算機視覺[17]，自然語言處理[18]、語音識別[19]以及航空航天領(lǐng)域的相關(guān)問題[20-27]。研究發(fā)現(xiàn)[28-35]，深度學(xué)習(xí)技術(shù)可在連續(xù)流領(lǐng)域?qū)崿F(xiàn)快速計算。Sekar等[28]使用卷積神經(jīng)網(wǎng)絡(luò)(CNN)和多層感知機(MLP)預(yù)測了不同幾何結(jié)構(gòu)機翼的升力系數(shù)：首先使用CNN由多種機翼構(gòu)型抽象出16個幾何和氣動參數(shù)，然后將這些參數(shù)輸入MLP網(wǎng)絡(luò)預(yù)測升力系數(shù)。結(jié)果表明，升力系數(shù)預(yù)測精度可達97%，并且計算速度相比傳統(tǒng)的計算流體力學(xué)(CFD)方法快了約150倍。Hui等[29]同時使用前饋神經(jīng)網(wǎng)絡(luò)和CNN計算了機翼附近的壓力分布。研究中機翼的幾何形狀被抽象為符號距離函數(shù)，并和CFD方法計算出的機翼附近的壓力分布數(shù)據(jù)共同輸入神經(jīng)網(wǎng)絡(luò)進行訓(xùn)練，以此預(yù)測壓力分布；結(jié)果顯示預(yù)測的均方根誤差在2%左右，并且其計算速度相比CFD方法提升了約3個量級。Wu等[30]使用數(shù)據(jù)驅(qū)動的方法，提出了一種數(shù)據(jù)增強的生成對抗神經(jīng)網(wǎng)絡(luò)(GAN)，用于機翼層流壓力分布的快速精確預(yù)測。該計算方法分為流場預(yù)訓(xùn)練模塊和微調(diào)模塊，預(yù)訓(xùn)練模塊使用條件GAN估計訓(xùn)練集的分布；微調(diào)模塊用于增強模型的泛化性能。結(jié)果顯示上述方法可以準(zhǔn)確計算機翼附近的流場和壓力系數(shù)。總而言之，深度學(xué)習(xí)技術(shù)已成功應(yīng)用于連續(xù)流領(lǐng)域，且其計算效率顯著優(yōu)于傳統(tǒng)CFD方法，為提高稀薄流領(lǐng)域的計算效率提供了可能，但當(dāng)前尚未開展針對真空羽流快速計算的研究。

本文提出了一種基于卷積神經(jīng)網(wǎng)絡(luò)的直接模擬蒙特卡洛方法(Convolutional Neural Networks-based Direct Simulation Monte Carlo method, CNN-DSMC)，實現(xiàn)了月面探測器月面著陸過程中真空羽流的快速計算。該方法使用卷積神經(jīng)網(wǎng)絡(luò)從DSMC數(shù)值模擬結(jié)果中提取特征，以此訓(xùn)練真空羽流智能計算模型，然后將DSMC數(shù)值模擬中不同的幾何拓撲和邊界條件輸入計算模型完成真空羽流場智能化求解。結(jié)果表明，和傳統(tǒng)的DSMC方法相比，CNN-DSMC方法可以在保持較高計算精度的同時顯著降低計算時間，大幅提高了真空羽流評估效率。

1 理論方法

本文以月面探測器月面著陸過程中真空羽流流程為例，分別通過基于卷積神經(jīng)網(wǎng)絡(luò)的直接模擬蒙特卡洛方法(CNN-DSMC)和DSMC方法模擬月面探測器在不同懸停高度時的真空羽流流場。本節(jié)主要闡述CNN-DSMC和DSMC方法。

1.1 基于卷積神經(jīng)網(wǎng)絡(luò)的直接模擬蒙特卡洛方法

圖1為CNN-DSMC方法求解流程示意圖。在CNN-DSMC方法中，計算分為數(shù)據(jù)預(yù)處理和模型訓(xùn)練2個過程。在數(shù)據(jù)預(yù)處理中，真空羽流仿真模型中的幾何拓撲信息被抽象為符號距離函數(shù)(Signed Distance Function, SDF)，邊界條件信息被抽象為標(biāo)識符矩陣(Identifier Matrix, IM)；SDF和IM共同作為訓(xùn)練集的輸入；將由DSMC數(shù)值模擬得到的真空羽流速度場(3個)和密度場作為訓(xùn)練集的輸出；測試集為未參與訓(xùn)練的DSMC數(shù)值模擬算例，用于驗證CNN-DSMC方法的準(zhǔn)確性。在完成訓(xùn)練之后，就得到了真空羽流智能計算模型：

(V,ρ)=f(D,M)

(1)

式中：V和ρ分別為真空羽流速度矢量和密度；D和M分別代表符號距離函數(shù)和標(biāo)識符矩陣。接下來分別對CNN-DSMC方法中的真空羽流仿真模型、卷積神經(jīng)網(wǎng)絡(luò)和數(shù)據(jù)預(yù)處理方法進行介紹。

1.2 真空羽流數(shù)值模擬模型

本文中DSMC算例均是通過北京航空航天大學(xué)羽流工作站PWS[6-7]完成的。PWS軟件中航天器面網(wǎng)格和DSMC計算的體網(wǎng)格是解耦的，并且采用了自適應(yīng)網(wǎng)格加密策略，方便計算各種復(fù)雜工況。此外，PWS軟件可以進行多核并行計算，且實現(xiàn)了真空羽流熱效應(yīng)分析[5]、力效應(yīng)分析[36-37]和污染效應(yīng)分析[38-39]等方面的數(shù)值模擬。地面實驗結(jié)果表明，PWS軟件數(shù)值模擬結(jié)果與實驗符合較好，可以滿足本文數(shù)據(jù)的精度要求，詳見附錄A。

本文研究的月面著陸過程羽流仿真所采用的計算域如圖2(a)所示。除了月面，計算域的其余5個邊界均設(shè)置為開放邊界。月面和月球探測器表面設(shè)置為固體邊界，熱適應(yīng)系數(shù)設(shè)置為1.0[40]。所有邊界的溫度設(shè)置為300 K，粒子與邊界相互作用模型使用Maxwell模型[41]。DSMC數(shù)值模擬的時間步長設(shè)置為10-7s。

定義月球探測器足墊到月面的距離為h，且以h= 10 m為界限區(qū)分羽流智能計算的近場模型和遠場模型。當(dāng)h≤ 10 m時，羽流流場中的激波相互作用非常復(fù)雜，流場形態(tài)比遠場情況更加混亂。因此，在近場情況下，由于計算域范圍較小且流場更為復(fù)雜，應(yīng)用CNN-DSMC方法對全局羽流流場進行訓(xùn)練；在遠場情況下，由于流場結(jié)構(gòu)更簡單，且計算域范圍較大，因此截取月面上方6 m處的流場用于CNN-DSMC訓(xùn)練，如圖2(b)所示。近場和遠場的訓(xùn)練集與測試集設(shè)置如表1所示。訓(xùn)練集和測試集選取依據(jù)詳見附錄B。

表1 訓(xùn)練集和測試集設(shè)置Table 1 Configurations of training and test datasets

1.3 卷積神經(jīng)網(wǎng)絡(luò)

CNN-DSMC方法中使用的整體網(wǎng)絡(luò)結(jié)構(gòu)如圖3(a)所示。該網(wǎng)絡(luò)由1個編碼器和2個解碼器組成，其中每個編碼器(解碼器)都由7個(反)卷積塊組成。每個(反)卷積塊的結(jié)構(gòu)組成包括3個(反)卷積層和1個最大(反)池化層，如圖3(b)所示。單個(反)卷積層包括(反)卷積、激活函數(shù)和批量正則化3個過程。卷積本質(zhì)上是一種矩陣變換，對于給定的矩陣A，卷積操作定義為[42]

B=wA

(2)

式中：B為卷積后得到的矩陣；w為卷積核，也是一個矩陣，其矩陣元素將在神經(jīng)網(wǎng)絡(luò)訓(xùn)練中不斷進行優(yōu)化。反卷積則是上述計算的逆過程。文中使用的激活函數(shù)為Relu[43]，其定義為

Relu=max(0,x)

(3)

式中：x為(反)卷積的輸出矩陣。批量正則化[44]主要用于修正每一層輸入數(shù)據(jù)的期望和方差，有利于訓(xùn)練過程的效率與穩(wěn)定性。最大池化層的作用本質(zhì)上是下采樣，而最大反池化層用于上采樣。

除了前饋過程，編碼器中每個卷積塊的輸出都會輸入解碼器中對應(yīng)位置的反卷積塊，該設(shè)置參考了ResNet[17]，是訓(xùn)練中避免梯度消失和梯度爆炸的關(guān)鍵操作。表2給出了CNN-DSMC方法中的參數(shù)設(shè)置，其中的參數(shù)是大量數(shù)值實驗后得到的較優(yōu)數(shù)值。卷積層結(jié)構(gòu){8, 16, 32, 32, 64, 64, 128}對應(yīng)7個卷積塊，對于反卷積塊，上述結(jié)構(gòu)反序。卷積核大小為5，即卷積核矩陣的形狀為5×5×5。在訓(xùn)練過程中，使用的優(yōu)化方法為AdamW[45]，其存在權(quán)重衰減功能，有助于提高訓(xùn)練過程的穩(wěn)定性。訓(xùn)練中使用的損失函數(shù)為真空羽流速度場和密度場的均方根誤差，表達式為

表2 CNN-DSMC方法中的參數(shù)設(shè)置Table 2 Parameter setting in CNN-DSMC

(4)

1.4 數(shù)據(jù)預(yù)處理

在CNN-DSMC方法中，真空羽流仿真模型中的幾何拓撲信息被抽象為SDF[24,46-47]。定義Ω為度量空間Y的一個子空間，度量空間Y代表了DSMC羽流數(shù)值模型的幾何拓撲信息，其本質(zhì)是一個矩陣，該矩陣元素的值僅包括-1和1兩種情況，在航天器內(nèi)部為-1，航天器外部為1。SDF被定義為[48]

(5)

式中：d為歐幾里得距離；?Ω為Ω的邊界；Ωc為Ω的內(nèi)部。對于任意y∈Y，有

(6)

式中：inf代表下確界。SDF的物理意義為空間中的某個點到邊界的最小距離，其符號由該點是否在邊界內(nèi)決定。式(5)中的D可以通過快速行進算法[48-49]得到。本文中輸入CNN-DSMC網(wǎng)絡(luò)的SDF在基于式(5)的情況下還進行了歸一化處理，如圖4(a)所示，月球探測器內(nèi)部區(qū)域的SDF為負值，月球探測器外部的SDF為正值，且值的大小隨著接近月球探測器邊界而減小。

在CNN-DSMC中，邊界條件信息被抽象為IM。IM本質(zhì)上是一個三維的矩陣，其矩陣元素作為區(qū)分三維空間不同區(qū)域的標(biāo)識符。在本文中，共選取4種不同的標(biāo)識符：開放邊界、航天器邊界、月面邊界和真空羽流區(qū)域，這4種標(biāo)識符也與實際DSMC數(shù)值模擬中使用的邊界條件相對應(yīng)，如圖4(b)所示。在具體設(shè)置中，除了月面，最外面的5個邊界的標(biāo)識符均設(shè)置為開放邊界；航天器表面及內(nèi)部設(shè)置為航天器邊界；月面設(shè)置為月面邊界；真空羽流流場設(shè)置為真空羽流區(qū)域。

遠場情況下，由于流場區(qū)域很大，所以只截取了月面上6 m作為訓(xùn)練區(qū)域(圖2(b))。在此種情況下，輸入CNN-DSMC中的信息額外包括了所截區(qū)域第1層網(wǎng)格的流場(速度、密度)信息，信息從h=30 m的算例中截取。

2 結(jié)果與討論

通過DSMC數(shù)值模擬獲得落月過程中遠場和近場的真空羽流流場數(shù)據(jù)，將其輸入CNN-DSMC網(wǎng)絡(luò)中進行訓(xùn)練，共訓(xùn)練40 000步，分別得到的遠場和近場流場訓(xùn)練過程中的損失函數(shù)曲線如圖5所示。結(jié)果表明，訓(xùn)練過程中，近場的損失函數(shù)下降更慢，且最終達到穩(wěn)定的值(約為10)也比遠場的值要大；遠場的損失函數(shù)最終穩(wěn)定在0.1附近。盡管遠場的實際區(qū)域大于近場，但遠場用于訓(xùn)練的區(qū)域高度僅有6 m,小于近場的訓(xùn)練區(qū)域高度10 m,且近場流態(tài)復(fù)雜，訓(xùn)練出的計算模型精度要低于遠場情況的模型，因而遠場訓(xùn)練損失函數(shù)的穩(wěn)定值更小。接下來將分別展示遠場和近場的結(jié)果對比。

2.1 遠場(h > 10 m)結(jié)果分析

圖6為遠場驗證算例(h=18 m)真空羽流速度場的DSMC數(shù)值模擬結(jié)果和CNN-DSMC計算結(jié)果的三維切片。圖中的u代表真空羽流在x

方向的分速度。由于遠場情況只截取了月面上方6 m處的流場(圖2(b))用于CNN-DSMC訓(xùn)練，在截取區(qū)域上方，兩者的流場是完全相同的，所以圖6只給出了這一局部范圍流場的對比。從圖中可以看出，整體上兩者得到的流場是非常相似的，而且流場計算結(jié)果表明CNN-DSMC可以正確地計算激波以及多條激波相互作用后的復(fù)雜流場。

利用CNN-DSMC還能同時得到真空羽流密度場，圖7給出了月面上方0.1 m處的密度分布。相比于速度場的結(jié)果，密度場的CNN-DSMC計算結(jié)果與DSMC數(shù)值模擬結(jié)果的一致性稍差一些，但是不一致的區(qū)域基本都位于非核心流域(>10 m)，這些區(qū)域的密度相比核心區(qū)的密度低了4～5倍，體現(xiàn)為式(4)中的損失函數(shù)密度項權(quán)重更小，因此計算誤差會稍大一些。

為了定量分析CNN-DSMC計算結(jié)果和DSMC數(shù)值模擬結(jié)果之間的差距，截取了月球探測器軸線(y=0 m,z=0 m)上的真空羽流速度和密度數(shù)據(jù)，如圖8所示。圖8(a)和圖8(b)分別為x方向速度u和密度隨距離的變化曲線，需要說明的是，給出的曲線在x方向的范圍為-10～-4 m，對應(yīng)于所截取的月面上方6 m的真空羽流流場；其中，-10 m對應(yīng)于月面的位置，-4 m對應(yīng)于截取平面的位置。圖8所示結(jié)果顯示CNN-DSMC計算結(jié)果與DSMC數(shù)值模擬結(jié)果幾乎相符，特別是激波附近處的結(jié)果也相符。

計算結(jié)果表明，CNN-DSMC計算和DSMC數(shù)值模擬得到的速度和密度的平均相對誤差分別為4.1%和8.2%。進一步計算相對誤差發(fā)現(xiàn)，流場速度和密度在核心區(qū)域的最大誤差分別為3.2%和9.9%，但在激波附近的最大誤差會增加到14.9%和28.9%。受限于DSMC計算時間，本文遠場條件下的訓(xùn)練集只包含6個算例(見表1)，因此CNN-DSMC計算結(jié)果局部偏差較大；但從圖8所示結(jié)果仍能看到CNN-DSMC方法的可行性。

2.2 近場(h ≤ 10 m)結(jié)果分析

圖9為近場驗證算例(h=8.2 m)真空羽流速度場的DSMC數(shù)值模擬結(jié)果和CNN-DSMC計算結(jié)果的對比。與圖6不同，由于近場情況使用全局真空羽流流場數(shù)據(jù)進行訓(xùn)練，因此圖9展示的結(jié)果為全局的流場結(jié)果。可以看出，兩者得到的速度場幾乎完全一致，CNN-DSMC計算的激波形狀也與DSMC數(shù)值模擬結(jié)果一致。同時，考慮到真空羽流近場結(jié)構(gòu)的復(fù)雜性，圖9還給出了流線的分布，結(jié)果表明DSMC數(shù)值模擬結(jié)果和CNN-DSMC計算的粒子運動軌跡也基本相同。

圖10所示為月面上方0.1 m處的真空羽流密度分布。與圖7展示的結(jié)果不同，近場情況下的CNN-DSMC計算出的密度分布和DSMC數(shù)值模擬結(jié)果差別更小，這是因為近場情況下的整體羽流密度(10-4)比遠場情況(10-5)高了一個量級，使得式(4)中的損失函數(shù)密度項權(quán)重更大。

圖11(a)和圖11(b)分別近場情況下月球探測器x方向速度和密度的變化曲線，范圍為-9～-1 m，其中：-9 m對應(yīng)于月面位置，-1 m對應(yīng)于月球探測器正下方與足墊同一高度的位置。結(jié)果表明，和遠場情況類似，CNN-DSMC計算得到的速度和密度與DSMC數(shù)值模擬結(jié)果基本一致；其平均相對誤差分別為6.0%和8.8%，在核心區(qū)的最大誤差分別為1.3%和4.5%，在激波附近的最大誤差分別為10.9%和14.6%。相比于遠場結(jié)果(圖8)，雖然近場結(jié)果的最大誤差有所減小，但平均相對誤差大于遠場結(jié)果。如前所述，近場情況下，流場結(jié)構(gòu)較為復(fù)雜，因此CNN-DSMC計算精度相對較低。

2.3 計算效率評估

表3為遠場情況下DSMC和CNN-DSMC計算時間對比。DSMC數(shù)值模擬的計算總時間取訓(xùn)練集中6個算例運行的平均時間。PWS使用的CPU型號為Intel Xeon E5-2670 v3 (2.3 GHz)，使用240核并行計算。CNN-DSMC的計算基于GPU，型號為Nvidia Quadro RTX-A6000。GPU的架構(gòu)使得CNN-DSMC計算可以同時處理多個算例，表3中Batch size代表算例數(shù)目。結(jié)果顯示，隨著算例數(shù)目的增加，CNN-DSMC計算總時間增加不多，但是平均單個算例時間從2.216 8 s降低到了0.087 3 s，相比傳統(tǒng)DSMC的加速比在1.62×105～4.11×106之間變化，極大地提高了真空羽流的計算效率。表4給出了近場情況下的DSMC和CNN-DSMC計算時間對比。由于近場的DSMC數(shù)值模擬實際的計算域比遠場小，所以DSMC數(shù)值模擬時間有所下降；且近場訓(xùn)練的區(qū)域要比遠場訓(xùn)練的區(qū)域(6 m)要大，所以CNN-DSMC的計算時間稍長一些。具體表現(xiàn)為，CNN-DSMC相比于DSMC數(shù)值模擬的加速比在6.06×104～9.03×105之間變化，同樣達到了很好的加速效果。近場情況下只統(tǒng)計到Batch size=20，這是由于所用GPU的顯存存在限制(48 GB)。

表3 遠場(h>10 m)情況下DSMC和CNN-DSMC計算時間對比Table 3 Comparison of computation time between DSMC and CNN-DSMC in far-field (h > 10 m) case

表4 近場(h≤10 m)情況下DSMC和CNN-DSMC計算時間對比Table 4 Comparison of computation time between DSMC and CNN-DSMC in near-field (h≤10 m) case

從以上分析中可以看出，CNN-DSMC方法不依賴于固定的航天器，也不僅僅局限于在化學(xué)推進羽流計算中使用。對于任意航天器的化學(xué)推進(或電推進)羽流問題，在CNN-DSMC方法訓(xùn)練完成后，只要輸入航天器的幾何外形和計算邊界條件，即可完成化學(xué)推進(或電推進)流場的快速計算。

3 結(jié) 論

本文提出了一種真空羽流智能化計算方法CNN-DSMC。該方法首先將真空羽流數(shù)值模擬模型中的幾何拓撲信息和邊界條件信息分別抽象為符號距離函數(shù)SDF和標(biāo)識符矩陣IM；然后將SDF、IM和DSMC數(shù)值模擬的流場數(shù)據(jù)作為訓(xùn)練集輸入CNN-DSMC中的卷積神經(jīng)網(wǎng)絡(luò)，訓(xùn)練得到真空羽流智能化計算模型，該模型可以用于化學(xué)推進和電推進的真空羽流高精度、高效率計算。

1) 結(jié)果顯示，CNN-DSMC方法計算出的真空羽流速度場和密度場無論是在遠場情況還是在近場情況下都與傳統(tǒng)的DSMC數(shù)值模擬結(jié)果基本一致。其中，在遠場情況下，CNN-DSMC計算的速度和密度與DSMC數(shù)值模擬結(jié)果的平均相對誤差分別為4.1%和8.2%，在近場情況下分別為6.0%和8.8%。

2) 在保證較高計算精度的同時，CNN-DSMC方法相比DSMC方法具有顯著的加速效果。在遠場情況下，加速比范圍為1.62×105～4.11×106；在近場情況下，加速比范圍為6.06×104～9.03×105。卓越的加速性能和較高的流場計算精度表明，CNN-DSMC是一個非常有應(yīng)用潛力的真空羽流智能化計算方法。