火星進入器作強迫震蕩運動壁面脈動壓力數(shù)值模擬

石小潘，趙 瑞，榮吉利，袁 武，李 齊

(1.北京理工大學飛行器動力學與控制教育部重點實驗室，北京 100081；2. 中國航空救生研究所，襄陽 441003；3. 中國科學院計算機網(wǎng)絡(luò)信息中心，北京 100190；4. 北京空間飛行器總體設(shè)計部，北京 100094)

0 引 言

長期以來，世界各國先后展開了多次的深空探測任務(wù)。其中包括火星、金星、木星、土衛(wèi)六等有大氣天體的無人進入或軟著陸[1]?；鹦亲鳛榫嚯x地球最近的地外行星，和其他行星相比其自然環(huán)境更接近地球，因此，火星成為人類探索宇宙重要的一部分。由于火星大氣比較稀薄，密度只有地球大氣的1%，幾乎所有的火星進入器均采用鈍體設(shè)計減速，沿雙曲軌道直接進入火星大氣[2]。進入器進入火星大氣層之后，由于進入器外形幾何剖面不連續(xù)過渡，流動繞過肩部時出現(xiàn)流動分離，在進入器尾段形成分離區(qū)，分離區(qū)內(nèi)旋渦流動的不穩(wěn)定性及強脫體渦的動力學特性使得作用于進入器上氣動力的脈動劇烈[3-4]，致使進入器的動態(tài)穩(wěn)定性遭到破壞，進入器可能會出現(xiàn)極限環(huán)震蕩的問題。進入器的非定常運動將會加劇進入器壁面的脈動壓力環(huán)境(或稱氣動噪聲環(huán)境)，使得進入器表面出現(xiàn)較大的局部載荷，引起進入器結(jié)構(gòu)的抖振響應(yīng)，極大地縮短材料的使用壽命。另外，脈動壓力信號以噪聲的形式通過透射及結(jié)構(gòu)共振轉(zhuǎn)變?yōu)榕搩?nèi)噪聲，極易引起進入器內(nèi)部儀器失效，導致飛行任務(wù)的失敗[5-6]。因此，研究火星進入器作非定常運動時壁面所受到的脈動壓力環(huán)境，對火星進入器的結(jié)構(gòu)設(shè)計及艙內(nèi)聲振環(huán)境的預測都具有重要的意義。

脈動壓力來源于飛行器邊界層的轉(zhuǎn)捩/湍流、分離/再附及激波震蕩，具有顯著的隨機特性[7]。對于脈動壓力的研究，目前常用的預示方法包括風洞實驗、經(jīng)驗公式和數(shù)值計算。自20世紀50年代以來，國內(nèi)外研究機構(gòu)對飛行器壁面脈動壓力作了諸多的研究。Plotkin等[8]通過對前人的研究成果進行分析歸納，整理了一套適用于預測航天飛機壁面脈動壓力的經(jīng)驗算法，給出了均方根脈動壓力系數(shù)、功率譜及空間相關(guān)性的計算方法。Laganelli等[9-10]分析研究了光滑/粗糙表面上附體和分離的邊界層流動，歸納了一套預測壁面脈動壓力分布的工程算法。徐立功等[11]對機動再入飛行器壁面所受的空氣動力脈動環(huán)境進行了分析，用統(tǒng)一的參變量給出了高超聲速范圍內(nèi)飛行器壁面脈動壓力統(tǒng)計特性的工程計算公式。張志成等[12]使用工程估算公式分析了來流馬赫數(shù)和來流攻角對球頭雙錐體表面脈動壓力分布的影響。王娜等[13]通過實驗對彈體模型表面的脈動壓力進行研究，分析了彈體表面上脈動壓力與馬赫數(shù)、迎角的關(guān)系以及脈動壓力的頻譜特征。趙瑞等[5,14]使用隱式大渦模擬方法(ILES)對跨聲速旋成體火箭外壁面的脈動壓力環(huán)境進行數(shù)值模擬，驗證了該方法的可行性，并對外壁面的流場結(jié)構(gòu)進行分區(qū)歸納，整理改進了各分區(qū)脈動壓力計算的經(jīng)驗公式。

對于返回艙這類飛行器，為了再入大氣層時減速和防熱的需要，一般采用大鈍頭倒錐外形[15]，其再入過程的穩(wěn)定性一直以來都是國內(nèi)外研究者首要關(guān)注的問題。根據(jù)已有的研究成果[16]，返回艙在亞跨超聲速階段可能會出現(xiàn)極限環(huán)震蕩。因此，對于這類飛行器壁面脈動壓力環(huán)境的研究，需要考慮飛行器非定常運動帶來的影響。Ross等[17]通過風洞實驗分析了雷諾數(shù)、來流馬赫數(shù)、攻角等對標準返回艙外形脈動壓力的影響。石小潘等[18]針對火星進入器在開傘之前的跨超聲速階段，采用DES方法分析了馬赫數(shù)、攻角及配平翼展開角對壁面脈動壓力環(huán)境的影響。以上對于標準返回艙及火星再入器壁面脈動壓力環(huán)境的研究均是基于飛行器保持相對靜止，由非定常流動誘導產(chǎn)生的脈動壓力，沒有考慮飛行器非定常運動帶來的影響。

綜上所述，國內(nèi)外對飛行器壁面脈動壓力環(huán)境作了較多的研究，并取得了顯著的成果。但對于鈍體外形的飛行器由于自身非定常運動對壁面脈動壓力環(huán)境的影響研究相對較少。本文基于DES方法，耦合進入器運動方程和流體力學方程，利用剛性動網(wǎng)格技術(shù)，對火星大氣下進入器作強迫俯仰震蕩運動進行非定常數(shù)值模擬，綜合分析進入器在開傘之前的超聲速氣動減速階段壁面的脈動壓力環(huán)境。

1 計算方法和模型

1.1 控制方程

流動控制方程為三維非定?？蓧嚎sN-S方程組，其在直角坐標系下的守恒積分形式可表示為：

(1)

式中：U為控制體，?U為控制面，Q為守恒變矢量，F(xiàn)為矢通量，n為邊界的外法向量。

采用有限體積法對控制方程進行數(shù)值離散，在數(shù)值離散過程中，空間對流項離散使用5階的WENO格式，黏性項使用4階的中心差分方法[19]。時間推進采用二階隱式雙時間法進行迭代[20]，內(nèi)迭代選擇常用的LU-SGS(lower-upper symmetric-Gauss-Seidel)方法[21]。定常計算采用一方程S-A湍流模型[22]并將其計算結(jié)果作為非定常DES計算的初始流場，DES方法的詳細介紹參考文獻[18]。

目前對于火星大氣介質(zhì)的模擬，有兩種方法：一種方法是采用真實氣體模型[23]，考慮氣體介質(zhì)中的各種化學組分；另一種方法是等效比熱比模型[24]。本文使用的是基于等效比熱比方法進行數(shù)值計算，其等效比熱比取為1.17[24]。另外，對于火星大氣環(huán)境，以CO2為主，其黏性系數(shù)隨溫度而改變，使用Sutherland公式來表征火星大氣以CO2為主，其黏性系數(shù)隨溫度而改變，使用Sutherland公式來表征

(2)

式中：參考黏性系數(shù)μ0=1.48×10-5kg/m·s，參考溫度T0=293.15 K。和黏性力一樣，氣體微團之間的熱傳導也是由分子運動造成，導熱系數(shù)和黏性系數(shù)之間存在相似關(guān)系，在處理高速黏性流動問題，導熱系數(shù)由Pr=μcp/k確定，Pr=0.71。

進入器強迫震蕩運動可認為是6自由度運動的簡化，忽略質(zhì)心的平動，僅考慮某一方向繞質(zhì)心的轉(zhuǎn)動。對于單自由度強迫俯仰震蕩，給定簡諧振動形式為：

α=αm+α0sin(kt)=αm+θ

(3)

式中：α為攻角，αm為初始攻角，α0為攻角振幅，θ為俯仰角，k=ωLref/2V∞為縮減頻率。本文計算和風洞實驗條件保持一致，強迫振動的縮減頻率為0.5 Hz，最大振幅2°。

1.2 動網(wǎng)格技術(shù)

由于進入器作強迫俯仰震蕩運動，網(wǎng)格模型將隨著改變，不同時刻必須重新生成新的網(wǎng)格，這將極大地增加計算的時間。采用剛性動網(wǎng)格技術(shù)，即令計算網(wǎng)格與物面剛性固連，整個非定常計算過程中，網(wǎng)格無須重新生成，這種方法動網(wǎng)格的質(zhì)量與靜網(wǎng)格質(zhì)量完全一致，精確地滿足幾何守恒定律。實現(xiàn)該方法較為簡單，在計算過程中對位置矢量乘以一個坐標轉(zhuǎn)化矩陣，n+1時刻的網(wǎng)格坐標：

Gn+1=Ln+1G0

(4)

式中：Ln+1為由n+1時刻姿態(tài)角確定的轉(zhuǎn)換矩陣；G0為姿態(tài)角全為0時的網(wǎng)格坐標。

1.3 計算模型和網(wǎng)格

本文研究進入器進入火星大氣層之后由于自身非定常運動誘發(fā)的壁面脈動壓力，同時作為對比，本文計算了進入器保持相對靜止僅由非定常分離流動引起的脈動壓力。計算工況如表1所示。

表1 進入器計算條件Table 1 Computation condition for the capsule

計算模型為類似文獻[25]中所示的具有配平翼的進入器外形，基本尺寸如圖1所示。配平翼展開后計算網(wǎng)格模型如圖2所示，網(wǎng)格量總計600萬，并在分離區(qū)附近(配平翼及艙體后方)進行加密，用于捕捉脫體渦結(jié)構(gòu)。壁面第一層網(wǎng)格保證y+<1，用于捕捉黏性邊界層流動。

圖1 進入器模型Fig.1 Physical model of the capsule

圖2 進入器計算網(wǎng)格示意圖Fig.2 Schematic of the capsule computational grid

對進入器z=0截面進行測點采樣，采樣點共計29個，艙體迎風面測點為1～14，艙體背風面測點為15～29。壁面脈動壓力的采樣頻率為200 Hz，從瞬時流場中提取壁面特征點壓力信號隨時間變化歷程曲線，獲得均方根脈動壓力系數(shù)，并采用快速傅里葉變換(FFT)，研究頻域空間中脈動壓力的功率譜密度分布。

1.4 數(shù)值方法驗證

為驗證本文使用的計算方法，以具有風洞實驗數(shù)據(jù)和計算數(shù)據(jù)的再入返回艙模型[26]作為驗證算例。風洞實驗模型及壁面測點如圖3所示，風洞實驗測試和計算條件見表2。針對該再入返回艙模型，基于SA-DES方法進行數(shù)值計算。

圖3 風洞實驗模型測點布置Fig.3 Observation points of wind tunnel model

MaU∞/(m·s-1)ρ∞/(kg·m3)T∞Re1.05341.3080.659263.031.35×1071.1326.080.0636218.671.21×1061.2380.3850.5752501.37×1071.4392.190.0428221.131.03×1062600.930.0246224.628.45×105

圖4和圖5分別給出了計算所得的壁面測點L1和W1的均方根脈動壓力系數(shù)、風洞實驗結(jié)果及Fujimoto等[26]的計算結(jié)果，通過對比可以看出數(shù)值計算所得的壁面測點的均方根脈動壓力系數(shù)與風洞實驗測量值及文獻中Fujimoto等[26]的計算值誤差相對較小，可以驗證本文所采用的數(shù)值計算方法具有一定的準確性。

圖4 L1測點計算與實驗均方根脈動壓力系數(shù)對比Fig.4 Surface pressure coefficient comparison with experiments at L1 position

圖5 W1測點計算與實驗均方根脈動壓力系數(shù)對比Fig.5 Surface fluctuating pressure coefficient comparison with experiments at W1 position

2 結(jié)果與分析

2.1 進入器非定常運動脈動壓力特征

圖6 來流馬赫數(shù)1.2不同狀態(tài)的各測點脈動壓力系數(shù)對比Fig.6 The coefficient of fluctuating pressure of different states under inflow Mach number is 1.2

在配平翼迎風面，進入器作強迫震蕩運動時，攻角震蕩加劇了分離區(qū)再附點的脈動，相比進入器保持相對靜止時，配平翼迎風面分離區(qū)的脈動壓力環(huán)境顯著增強。圖7和圖8分別為來流馬赫數(shù)為1.2時進入器保持相對靜止和強迫震蕩運動時不同時刻流場結(jié)構(gòu)示意圖。對于進入器背風面的脈動壓力環(huán)境，分離區(qū)內(nèi)渦流脈動是壁面壓力脈動的主要原因[18]。從圖7可以看出，進入器保持相對靜止時，分離區(qū)剪切層位置基本一致，流場內(nèi)大尺度結(jié)構(gòu)變化較小，脈動壓力主要由分離區(qū)內(nèi)的小尺度渦流脈動所致。而進入器作強迫震蕩運動時，流動分離區(qū)內(nèi)大尺度結(jié)構(gòu)非定常效應(yīng)顯著、變化劇烈(如圖8所示)，致使渦流脈動更加劇烈。因此，進入器作強迫震蕩運動時加劇了進入器背風壁面的脈動壓力環(huán)境，使得均方根脈動壓力系數(shù)為0.04～0.07，增加了近一倍。

圖7 Ma=1.2進入器相對靜止時不同時刻流場結(jié)構(gòu)示意圖Fig.7 The sketch of flowfield structure of different times at inflow Mach number 1.2 under relatively static condition

圖8 Ma=1.2進入器強迫震蕩時不同時刻流場結(jié)構(gòu)示意圖Fig.8 The sketch of flowfield structure of different times at inflow Mach number 1.2 under forced oscillation

圖9為來流馬赫數(shù)為3時進入器作強迫震蕩運動與進入器保持相對靜止時壁面測點均方根脈動壓力系數(shù)對比曲線。從圖9可以看出，進入器作強迫震蕩運動使得艙體迎風面的均方根脈動壓力系數(shù)明顯增加。

圖9 來流馬赫數(shù)3不同狀態(tài)的各測點脈動壓力系數(shù)對比Fig.9 The coefficient of fluctuating pressure of different states under inflow Mach number is 3

圖10和圖11分別為來流馬赫數(shù)為3時進入器保持相對靜止和強迫震蕩運動時不同時刻流場結(jié)構(gòu)示意圖。通過對比可以發(fā)現(xiàn)，進入器保持相對靜止脫體激波形態(tài)和位置保持穩(wěn)定時，對迎風面的壓力脈動起到抑制作用，均方根脈動壓力系數(shù)約為0。進入器作強迫震蕩運動時，迎風面脫體激波形態(tài)和位置發(fā)生改變，脫體激波震蕩將會在艙體迎風面誘導產(chǎn)生均方根脈動壓力系數(shù)在0.03～0.09的脈動壓力環(huán)境。對于進入器背風面，相比進入器保持相對靜止，進入器強迫震蕩運動時迎風面脫體激波震蕩使得背風面高速區(qū)范圍增大，低速分離區(qū)范圍減小(如圖11所示)。分離區(qū)范圍減小使得渦流脈動減弱，但背風面分離區(qū)脈動壓力環(huán)境同樣受到迎風面脫體激波震蕩的影響，兩者共同作用下，使得進入器背風面脈動壓力環(huán)境和進入器保持相對靜止時差異不大，均方根脈動壓力系數(shù)在0.01以下。

圖10 Ma=3進入器相對靜止時不同時刻流場結(jié)構(gòu)示意圖Fig.10 The sketch of flowfield structure of different times at inflow Mach number 3 under relatively static condition

圖11 Ma=3進入器強迫震蕩時不同時刻流場結(jié)構(gòu)示意圖Fig.11 The sketch of flowfield structure of different times at inflow Mach number 3 under forced oscillation

圖12為來流馬赫數(shù)為1.2和3時進入器艙翼連接段的局部流場圖，可以看出在配平翼迎風面靠近翼根區(qū)存在小分離區(qū)，進入器保持相對靜止時小分離區(qū)非定常效應(yīng)顯著，再附點劇烈脈動會在配平翼迎風面誘導產(chǎn)生均方根脈動壓力系數(shù)為0.022(Ma=1.2)和0.002(Ma=3)的脈動壓力環(huán)境[18]，而進入器作強迫震蕩運動時，在配平翼迎風面的誘導產(chǎn)生均方根脈動壓力系數(shù)高達0.06(Ma=1.2)和0.12(Ma=3)的脈動壓力環(huán)境。通過對比可以看出，低馬赫數(shù)時，進入器強迫震蕩運動使得均方根脈動壓力系數(shù)增大約3倍；高馬赫數(shù)時，進入器作強迫震蕩運動使得均方根脈動壓力系數(shù)增大約60倍。進入器強迫震蕩運動時的攻角震蕩和激波震蕩加劇了翼根分離區(qū)的非定常脈動，造成更為惡劣的脈動壓力環(huán)境，在結(jié)構(gòu)設(shè)計中需要尤為關(guān)注。

圖12 進入器作強迫震蕩時翼根區(qū)局部流場圖Fig.12 The local flowfield of the flap-root under forced oscillation

圖13 進入器強迫震蕩時迎風面三個測點的脈動壓力變化圖Fig.13 The fluctuating pressure at three points around the windward of the capsule under forced oscillation

圖14 進入器強迫震蕩時配平翼迎風面測點脈動壓力變化圖(測點11,Ma=3)Fig.14 The fluctuating pressure for the windward side of the capsule under forced oscillation (Point 11,Ma=3)

2.2 進入器非定常運動脈動壓力的功率譜密度特征

功率譜密度表征脈動能量隨頻率的分布特性。通過對火星進入器脈動壓力的頻譜特性分析，可以確定脈動能量集中的頻率區(qū)間，進而對進入器進行針對性的結(jié)構(gòu)設(shè)計。將進入器分為三個區(qū)域：艙體迎風面，配平翼迎風面，進入器背風面，選取各區(qū)域特征點進行功率譜密度分析。圖15為上述3個區(qū)域特征測點的功率譜密度頻譜在不同來流馬赫數(shù)下的變化曲線。測點5位置(圖15(a))位于艙體迎風面大底中心，來流馬赫數(shù)較小時，脫體激波較弱且距離進入器較遠，脈動能量集中在0.5 Hz，和攻角震蕩頻率一致；來流馬赫數(shù)較大時，進入器迎風面脈動壓力環(huán)境受攻角震蕩和激波震蕩共同影響，且激波震蕩的影響要大于攻角震蕩的影響，脈動能量密度在30 Hz左右達到峰值，激波震蕩誘導的脈動壓力在頻域表現(xiàn)出低頻高幅的特性，在結(jié)構(gòu)設(shè)計中要著重注意。測點11位置(圖15(b))位于配平翼迎風面分離區(qū)，來流馬赫數(shù)較小時，脈動能量主要集中在0.5 Hz，但相比大底迎風面，脈動能量幅值要高出幾個量級。來流馬赫數(shù)較大時，功率譜密度和大底迎風面類似，但由于該區(qū)域存在激波與分離區(qū)再附點相互作用，導致脈動能量幅值比大底迎風面要大。測點17位置(圖15(c))位于進入器背風面分離區(qū)，脈動壓力主要受艙體分離區(qū)內(nèi)渦流脈動和迎風面激波震蕩雙重影響，低馬赫數(shù)時背風面分離區(qū)較大，渦流脈動較強，致使10 Hz左右的脈動能量密度幅值較高；高馬赫數(shù)時，進入器背風面分離區(qū)減小，渦流脈動較弱，致使10 Hz左右的脈動能量密度有所減弱，但是激波震蕩對于背風面分離區(qū)影響增強，使得脈動能量密度在30～40 Hz范圍有所增強。

圖15 進入器特征點的頻譜分布Fig.15 Power spectral density distributions at the capsule

3 結(jié) 論

本文以火星進入器為研究對象，使用強迫震蕩方法分析了進入器自身非定常運動時壁面脈動壓力環(huán)境的變化規(guī)律，通過與進入器保持相對靜止時壁面脈動壓力環(huán)境作了對比，得出以下結(jié)論：

1)與靜止狀態(tài)不同，進入器作強迫震蕩運動時，由于攻角震蕩和脫體激波震蕩，會引起進入器迎風面。

2) 當來流馬赫數(shù)較小時(Ma=1.2)，進入器迎風面脈動壓力環(huán)境主要受攻角震蕩的影響，脈動能量集中頻率和攻角震蕩頻率保持一致。由于攻角震蕩的影響，背風面分離區(qū)流場結(jié)構(gòu)變化劇烈，致使渦流脈動更加劇烈。進入器強迫震蕩運動加劇了進入器迎風面和背風面的脈動壓力環(huán)境。

3) 當來流馬赫數(shù)較大時(Ma=3)，進入器迎風面脫體激波強度較強，由于運動脫體激波震蕩將在配平翼迎風面產(chǎn)生更為惡劣的脈動壓力環(huán)境，均方根脈動壓力系數(shù)高達0.12；對于進入器背風面，脫體激波震蕩使得背風面高速區(qū)范圍增大，低速分離區(qū)范圍減小，抑制了背風面渦流脈動，使得該區(qū)域脈動壓力環(huán)境減緩。