火箭整流罩脈動壓力環(huán)境數(shù)值模擬與優(yōu)化設計

樊宇翔，于 洋，席 柯，趙 瑞，*，任 方

(1.北京理工大學宇航學院，北京 100081；2.北京強度環(huán)境研究所可靠性與環(huán)境工程技術(shù)重點實驗室，北京 100076；3.中國兵器工業(yè)導航與控制技術(shù)研究所，北京 100089)

0 引言

飛行器在主動飛行階段，整流罩外脈動壓力會對壁面產(chǎn)生隨時間變化的激勵力，這是飛行器結(jié)構(gòu)設計中需要考慮的關(guān)鍵問題之一。大量的實驗研究表明，氣動噪聲是一種強噪聲，聲壓總級高達130～170 dB，低中高頻都有所覆蓋，其不但會破壞飛行器結(jié)構(gòu)，還可能以聲透射方式進入飛行器內(nèi)部，影響內(nèi)部儀器正常工作并對人員造成傷害，甚至引發(fā)嚴重的飛行事故[1-3]。因此，整流罩作為保護飛行器的重要結(jié)構(gòu)，對改善其氣動噪聲環(huán)境具有重要意義。整流罩最惡劣的噪聲環(huán)境主要有兩個階段：第一階段是發(fā)射階段的噴流噪聲；第二階段則是在飛行器主動飛行時，整流罩表面由于非連續(xù)過渡的存在，流動容易出現(xiàn)膨脹、分離、壓縮等結(jié)構(gòu)，這些流動結(jié)構(gòu)的相互影響，尤其是在跨聲速階段，會使得整流罩受到極大的脈動壓力，形成極為惡劣的噪聲環(huán)境[4-6]。

目前，對飛行器脈動壓力環(huán)境預測的方法主要有數(shù)值模擬法、經(jīng)驗/半經(jīng)驗工程算法和風洞試驗法。對于數(shù)值模擬法，Tsutsumi等[7]采用改進的延遲脫體渦模擬方法（improved delayed detached-eddy simulation,IDDES）[8]獲取了典型整流罩結(jié)構(gòu)的脈動壓力等特征參數(shù)。趙瑞等[9]基于五階加權(quán)本質(zhì)無振蕩（weighted essentially non-oscillatory,WENO）格式構(gòu)造隱式大渦模擬方法（implicit large eddy simulation,ILES），對跨聲速來流條件下火箭整流罩外頭錐部位噪聲環(huán)境進行了數(shù)值模擬計算，驗證了該計算方法的可行性，并通過瞬時流場分析，改進了脈動壓力經(jīng)驗預測公式[10]。此外，他們還對比了整流罩頭錐不同母線形狀對脈動壓力環(huán)境的影響效果，指出馮·卡門母線形狀有利于減緩頭錐脈動壓力環(huán)境[11]。石小潘等[12]采用脫體渦方法（detached eddy simulation,DES）預示了火星進入器在跨聲速階段時的脈動壓力環(huán)境。之后他們又采用動網(wǎng)格技術(shù)，耦合了剛體運動方程，對火星進入器強迫震蕩運動開展了非定常數(shù)值模擬，指出強迫振蕩會誘導大底脫體激波進行同頻率振動，進而導致進入器大底的脈動壓力環(huán)境極為惡劣[13]。

綜上所述，目前研究飛行器壁面脈動壓力環(huán)境的數(shù)值模擬方法主要依靠DES和LES。但以上兩種方法主要針對分離流動，而實際飛行器繞流往往以附體流動為主導，DES和LES需要在邊界層內(nèi)布置較密的網(wǎng)格才能捕捉小尺度湍流脈動結(jié)構(gòu)，這樣會耗費大量的計算資源[14]。為解決上述問題，Morris等[15]提出使用RANS模型先求解流場，再求解非線性擾動方程以獲得聲場，該方法準確地計算了三維射流聲場問題。Batten等[14]將上述方法進行了改進，提出一種非線性噪聲求解器（non-linear acoustic solver, NLAS）方法。這種方法采用非線性湍流模型（nonlinear reynolds averaged Navier-Stokes simulation, Nonlinear RANS）進行流場求解，求解效率較高；而無法求解的小尺度量則進行湍流重構(gòu)[16-17]，兼顧了亞格子尺度聲源的影響。因此相比于LES與DES方法，NLAS在保證低耗散的同時降低了網(wǎng)格量。近年來，NLAS方法被成功應用于空腔噪聲模擬[18-20]、高鐵噪聲模擬[21]、飛機機翼噪聲模擬[22-23]、射流噪聲模擬[24]等。

雖然目前關(guān)于NLAS方法的應用已經(jīng)取得了一些進展，但將其應用于整流罩脈動壓力環(huán)境的研究還較少。任淑杰等[25]使用RANS/NLAS方法研究了跨聲速階段火箭表面典型位置處的脈動壓力，給出了不同馬赫數(shù)下的聲壓級與均方根脈動壓力系數(shù)，但是并未給出直觀的流動現(xiàn)象，也沒有分析迎角的影響。因此，本文使用RANS/NLAS方法，繼續(xù)深入研究整流罩的壁面脈動壓力環(huán)境產(chǎn)生機理并進行優(yōu)化設計。文章首先對典型的整流罩頭錐外形進行方法驗證；然后研究頭錐-柱-倒錐外形整流罩的脈動壓力環(huán)境，分析迎角對脈動壓力的影響規(guī)律；最后采用優(yōu)化倒錐型線的方式抑制脈動壓力環(huán)境。

1 控制方程

1.1 湍流模型

為準確描述NLAS方程中的非線性項，需采用一種非線性的湍流模型獲取雷諾應力分布。本文用兩方程cubick-ε湍流模型[26-29]，同時引入Favre平均[30-33]。即在雷諾N-S（RANS）方程組的基礎(chǔ)上提出密度加權(quán)平均概念，對壓強和密度采用時間平均，其他變量采用密度加權(quán)平均，存在以下Favre分解：f=。Favre平均解決了RANS方程組處理可壓縮流動較為復雜的問題，通適于不可壓縮流動與可壓縮流動，其輸運方程如下：

式中：常量Cε1= 1.44，Cε2= 1.92，σε= 1.3，σk= 1.0，

系數(shù)分別為：

其中：Pk為湍動能生成項，ρ為密度，t為時間，ui為i方向上的速度分量，μ為黏性系數(shù)，μt為渦黏性系數(shù)，Tt為時間尺度，常量AEτ= 0.15，Aμ= 0.01，Cτ=

1.2 NLAS控制方程

NLAS將湍流模型獲取的雷諾應力和重構(gòu)的小尺度亞格子源項共同作為聲學擾動方程的源項，通過求解聲學擾動方程，模擬脈動壓力環(huán)境。

將流體變量分為統(tǒng)計平均項和擾動項：將其代入Navier-Stokes方程組并進行重新整理，等式左邊為擾動項，右邊為平均項，可得非線性的擾動方程：

式中：

其中：p為壓力；θi為傳熱項；δij為克羅內(nèi)克符號；τ為剪切應力；e為單位體積能，表達式為：

忽略式（4）中的密度擾動并進行時間平均處理，得：

式中：LHS與RHS分別代表式（4）的左項與右項，Ri是標準雷諾應力張量和湍流熱通量相關(guān)項，由左式擾動項化解而來，表達式為：

式中cp為定壓比熱容。Ri中各項由之前的cubick-ε求解得出。

對于無法求解的小尺度脈動量，需要使用湍流統(tǒng)計結(jié)果重構(gòu)作為亞格子源項。湍流重構(gòu)指的是在cubick-ε求解得到的平均流場基礎(chǔ)上產(chǎn)生對應的擾動速度場[23]，再疊加到非線性擾動方程求解過程中。在LES和聲學相關(guān)領(lǐng)域，湍流重構(gòu)可以顯著降低計算成本[17]。Kraichnan[16]在1970年就提出了一種詳細的湍流重構(gòu)方法，用傅里葉模式疊加描述時間與空間關(guān)系，不過該方法只適用于各向同性湍流。之后Smirnov等[17]提出了基于雷諾應力張量相似變換的張量尺度概念，對Kraichnan的方法進行了改進，將其成功應用到了各向異性湍流。最后Batten等[14]對張量尺度進行了簡化處理，避免了計算相似變換的過程，其湍流擾動速度的傅里葉模式方程如下：

式中：

其中：為網(wǎng)格中特定位置的速度擾動；xj為笛卡爾坐標系中的位置；aik為當?shù)乩字Z應力張量的Cholesky分解；L為湍流長度尺度；τ為湍流時間尺度；與為余弦和正弦函數(shù)的系數(shù)；cn為速度縮放尺度；向量積運算的置換張量； 是高斯正態(tài)分布函數(shù)，其平均值為，標準偏差為β。上述湍流重構(gòu)過程都是以cubick-ε計算的相關(guān)張量和湍流長度尺度與時間尺度作為輸入，重構(gòu)的結(jié)果將在式（4）時間迭代過程中作為擾動項加入。

此外，NLAS可以使用原來的RANS計算網(wǎng)格來模擬壓力脈動的傳播，也可以通過插值在新的計算網(wǎng)格上進行求解。相對于RANS，NLAS對網(wǎng)格的要求較低，可大量降低網(wǎng)格量，節(jié)省計算成本，更好地適用于工程復雜外形?？臻g離散采用耦合TVD限制器的二階迎風格式，其限制器為Min-mod，聲場計算的時間積分采用雙時間步長隱式方法。迭代時，采用多重網(wǎng)格法加速收斂的同時，使用MPI并行計算，減小計算時間，提高計算效率。

2 NLAS方法驗證

2.1 驗證模型

采用文獻[7,9]研究的整流罩外形作為驗證算例，模型如圖1、圖2所示。計算迎角為0°，馬赫數(shù)為0.8，基于圓柱段直徑的雷諾數(shù)為2.66×106，計算網(wǎng)格參照文獻[9]，采用1/4模型，第一層網(wǎng)格高度y+＜1。文獻[7]采用1/72模型，網(wǎng)格量為1100萬。在RANS計算中，遠場邊界為特征邊界條件，對于NLAS聲學計算，遠場為吸收層邊界，以消除遠場邊界反射，時間步長2×10-7s。為避免初始計算數(shù)值擾動，舍棄掉初始3000步數(shù)據(jù)。

圖1 整流罩模型的幾何尺寸[7]Fig.1 Rocket fairing model geometry[7]

圖2 風洞試驗模型[7]Fig.2 Wind tunnel model of the rocket fairing[7]

2.2 計算結(jié)果對比

圖3 壁面壓力系數(shù)（時均）軸向分布對比Fig.3 Comparison of time-averaged surface pressure coefficient Cp profilesalong the axis-direction

圖4 再附點附近的流向速度型（時均）分布對比Fig.4 Comparison of time-averaged streamwise velocity profiles around the reattachment point

圖5均方根脈動壓力系數(shù)分布對比Fig.5 Comparison of root-mean-square fluctuating pressure coefficient Cp_rms profilesalong the axis-direction

圖3 ～圖5為RANS/NLAS聲場計算結(jié)果與文獻結(jié)果的對比。其中，文獻[7]和文獻[9]分別使用IDDES方法和基于五階WENO格式的ILES方法開展了數(shù)值計算。圖3為時均后壁面壓力系數(shù)軸向分布圖，通過與試驗值和文獻計算結(jié)果進行對比發(fā)現(xiàn)，使用NLAS方法能以較少的網(wǎng)格量獲得更優(yōu)的結(jié)果，與試驗值也更為接近。圖4為時均后再附點附近流向速度型分布對比圖，從圖中同樣可以看到，本文使用的NLAS方法能夠更加準確地預測出速度分布。圖5為整流罩壁面均方根脈動壓力系數(shù)分布對比圖。均方根脈動壓力系數(shù)Cp_rms的計算公式為：Cp_rms=是表征噪聲環(huán)境的重要參數(shù)；其中，是均方根脈動壓力，代表壁面上某點脈動壓力的總強度，q∞=0.5ρ∞是來流動壓頭?？梢钥吹?，均方根脈動壓力系數(shù)在折角區(qū)內(nèi)顯著增加，過折角之后開始下降，附體之后與實驗值更吻合。頭錐折角區(qū)域為激波/邊界層干擾的主要區(qū)域，圖中試驗值的差異是由于Kulite傳感器布局的限制，使其并沒有捕捉到折角區(qū)域的峰值[7]。不同數(shù)值模擬方法都有預測到這一段區(qū)域的峰值，峰值的具體大小與位置的差異應該與網(wǎng)格量有關(guān)[7]。相比于文獻[9]的基于五階WENO格式構(gòu)造的ILES方法，以及文獻[7]的使用1100萬網(wǎng)格量（1/72模型）的IDDES方法，本文采用更少的網(wǎng)格量以低階計算格式進行計算，計算效率更高，計算精度較好。

3 整流罩外流場數(shù)值模擬分析與型線優(yōu)化設計

3.1 幾何模型和計算網(wǎng)格

計算模型為頭錐-柱-倒錐的整流罩風洞模型，模型已做特定縮放處理，圖6為外形示意圖。

圖6 整流罩外形Fig.6 Rocket fairing shape

計算模型為1/2模型，圖7為RANS計算與NLAS計算的網(wǎng)格對比，RANS計算網(wǎng)格量約為620萬，第一層網(wǎng)格高度y+＜1，NLAS計算網(wǎng)格量約為370萬。邊界條件與算例驗證相同，自由來流馬赫數(shù)為0.8，單位雷諾數(shù)為1.58×107/m，迎角為0°、4°、8°，時間步長為 5×10-6s，共計算18000步，舍棄初始3000步數(shù)據(jù)。

圖7 RANS與NLAS網(wǎng)格比較Fig.7 Comparison of meshes used by RANSand NLAS

3.2 計算結(jié)果分析

3.2.1 定常RANS流場分析

圖8為迎角8°下的整流罩流場，對稱面為馬赫數(shù)云圖，壁面為無量綱壓力云圖。亞聲速來流沿整流罩頭部加速，使得壁面壓力不斷降低，在頭錐和圓柱段交界處即肩部附近顯著降低達到最小值，過折角后膨脹加速形成一道斜激波，過激波后速度降低為亞聲速，壁面壓力增大。在肩部區(qū)域因為激波強逆壓梯度發(fā)生分離，形成分離區(qū)，產(chǎn)生正激波，兩道激波構(gòu)成“λ”型，如圖9的肩部放大圖所示，對稱面改用無量綱密度云圖著色。而在倒錐處，由于外形突然的收縮產(chǎn)生大分離區(qū)，如圖10所示。此外由于有迎角的存在，在肩部附近，背風面的氣流迅速加速，流速更快，形成更大的高速區(qū)域；在倒錐處，迎風面的分離區(qū)更大，再附點靠后，背風面的分離區(qū)減小。

圖8 壁面壓力與對稱面馬赫數(shù)云圖Fig.8 Contours for the pressure on the wall and the Mach number in the symmetry plane

圖9 背風面肩部密度云圖Fig.9 Density contour in theleeward side of the cone-cylinder

圖10 背風面倒錐部位分離區(qū)Fig.10 Separation zone in the leeward side of the inverted cone

3.2.2 NLAS脈動壓力計算分析

迎角的存在導致流動在背風面肩部產(chǎn)生了更強的逆壓梯度，附體邊界層發(fā)生明顯分離，脈動壓力環(huán)境惡劣，因此本次研究主要針對背風面。

圖11為NLAS計算所得的整流罩肩部瞬時數(shù)值紋影圖，可以觀察到其上游的斜激波主要在圓柱段與頭錐連接的拐角處附近；分離剪切層不穩(wěn)定、極易破碎，并向上游傳播，其與激波相互作用，使得分離區(qū)的范圍和激波位置前后震蕩；在激波/邊界層的干擾下，形成了極其惡劣的脈動壓力環(huán)境。

圖11肩部瞬時數(shù)值紋影圖Fig.11 Instantaneous numerical schlieren images of the cone-cylinder

圖12 為整流罩背風面均方根脈動壓力系數(shù)沿軸向分布的風洞試驗與數(shù)值計算結(jié)果對比圖，各狀態(tài)下計算結(jié)果和試驗數(shù)據(jù)均已經(jīng)過特定參考值無量綱化處理（下文同）。肩部處由于激波/邊界層的干擾，產(chǎn)生了較大的峰值。倒錐處，分離區(qū)的不斷擾動使得脈動壓力增大，導致了倒錐部峰值的出現(xiàn)。NLAS預測結(jié)果與試驗值較為吻合。

圖12背風面均方根脈動壓力系數(shù)分布(Ma = 0.8，α = 8°)Fig.12 Root-mean-square fluctuating pressure coefficient Cp_rms distributionsin the leeward side(Ma = 0.8，α = 8°)

圖13 為不同迎角下背風面肩部附近的均方根脈動壓力系數(shù)分布圖。脈動壓力系數(shù)隨迎角的變化比較復雜，整體有隨迎角增大而增大的趨勢，這主要是因為隨著迎角的增大，肩部激波/邊界層的影響范圍增大，導致了更加惡劣的脈動壓力環(huán)境。

圖14為不同迎角下背風面倒錐的均方根脈動壓力系數(shù)分布圖，可以看出在倒錐部位，隨著迎角的增大，均方根脈動壓力系數(shù)會減小，這主要是因為背風面倒錐部位的分離區(qū)會因迎角的增大而減小，如圖15所示，分離區(qū)產(chǎn)生的擾動作用減弱，從而使得均方根脈動壓力系數(shù)值降低。

圖13 不同迎角下肩部均方根脈動壓力系數(shù)分布(Ma = 0.8)Fig.13 Root-mean-square fluctuating pressure coefficient Cp_rms distributions on the cone-cylinder at different angles of attack(Ma = 0.8)

圖14 不同迎角下倒錐處均方根脈動壓力系數(shù)(Ma = 0.8)Fig.14 Root-mean-squarepressure coefficient Cp_rms distributionson the inverted cone at different anglesof attack(Ma = 0.8)

圖15 不同迎角下倒錐處分離區(qū)大小Fig. 15 Separation zone on the inverted cone at different angles of attack

3.3 倒錐型線優(yōu)化設計

如圖16所示，為抑制整流罩倒錐部位的脈動壓力，改善其噪聲環(huán)境，現(xiàn)設計三種不同的方案進行優(yōu)化。方案一：采用一段直線連接兩端；方案二：采用分別與兩端相切的正弦曲線連接兩端；方案三：采用“切線弧+圓弧”將兩端連接，尾端的圓弧半徑為1.8×10-2m。來流工況見3.1節(jié)，迎角為8°。

圖17為不同外形背風面時均摩擦力系數(shù)Cf軸向分布圖，其計算公式為：Cf=其中τwall是壁面剪切應力。通過對比可以看出，“切線弧+圓弧”外形的摩擦力系數(shù)下降緩慢，到達零值后小幅下降，然后再迅速回升到零值，流動緩慢減速分離，分離點靠后，分離效果較弱，分離區(qū)最小；“正弦”外形相較于“切線弧+圓弧”外形效果較差；原始外形和“直線”外形因在起始端的非光滑過渡，導致摩擦力系數(shù)迅速下降，而原始外形由于中間段外形的改變，使得摩擦力系數(shù)有所回升，兩種外形的分離區(qū)都較大，分離效果較強。各方案的分離區(qū)情況如圖18所示。

圖16 倒錐部位輪廓線Fig.16 Outline profiles of the inverted cone

圖17 不同外形背風面倒錐處時均摩擦力系數(shù)分布Fig.17 Time-averaged friction coefficients Cf distributions in the leeward side of different inverted cones

圖18 不同倒錐部位分離區(qū)Fig.18 Separation zonesof different inverted cones

圖19為不同外形背風面均方根脈動壓力系數(shù)軸向分布對比圖。對于“直線”外形與原始外形，在起始端，強烈的流動分離導致該區(qū)域壓力脈動劇烈，均方根脈動壓力系數(shù)增大；“直線”外形在肩部之后為直線過渡，脈動壓力減弱，而原始外形由于中間依舊存在折角的非光滑過渡，數(shù)值波動上升，達到峰值；在尾端，由于非光滑過渡和分離再附的影響兩者均再次出現(xiàn)峰值。對于“正弦”外形與“切線弧+圓弧”外形，在起始端，均方根脈動壓力系數(shù)增長相對緩慢且幅值較??；但在中間段，“正弦”外形斜率較大，會引起較強的脈動壓力，達到峰值，而“切線弧+圓弧”外形的連接更為平緩，數(shù)值更低；在尾端，“正弦”外形的相切連接減弱了分離再附的劇烈脈動，均方根脈動壓力系數(shù)快速下降，而“切線弧+圓弧”外形在尾端出現(xiàn)峰值。尾端由于分離再附的脈動影響較為劇烈，數(shù)值都相對較高。綜合來看，“切線弧+圓弧”外形方案的優(yōu)化效果更好。

圖19 不同外形背風面均方根脈動壓力系數(shù)分布Fig.19 Root-mean-squarepressure coefficient Cp_rms distributions in the leeward side of different inverted cones

4 結(jié)論

本文使用RANS/NLAS的非線性噪聲求解方法對跨聲速條件下火箭整流罩外部脈動壓力環(huán)境進行了數(shù)值模擬計算，研究了迎角對脈動壓力的影響，并采用優(yōu)化倒錐型線的方式對倒錐部位的脈動壓力環(huán)境進行了改善。通過算例驗證，與試驗結(jié)果和文獻結(jié)果對比分析，得出以下結(jié)論：

1）RANS/NLAS可以準確預測跨聲速條件下整流罩外的非定常流動現(xiàn)象；并且由于NLAS方法對網(wǎng)格要求較低，可減少聲場計算的網(wǎng)格數(shù)量，降低計算成本，所以其適合于工程復雜外形脈動壓力環(huán)境計算。

2）在跨聲速來流條件下，氣流經(jīng)過火箭整流罩的頭錐肩部后會形成“λ”型激波和分離泡；在倒錐部位會產(chǎn)生大范圍的分離區(qū)，導致兩區(qū)域的脈動壓力環(huán)境極為惡劣，均方根脈動壓力系數(shù)達到峰值。

3）隨著迎角的增大，背風面肩部激波/邊界層影響的范圍增大，均方根脈動壓力系數(shù)整體有增大的趨勢；倒錐處，背風面的分離區(qū)減小，由此產(chǎn)生的擾動作用減弱，均方根脈動壓力系數(shù)有所降低。

4）采用優(yōu)化倒錐型線的設計方式來抑制倒錐部位的脈動壓力，分析了四種外形下的時均摩擦力系數(shù)、分離區(qū)、均方根脈動壓力系數(shù)，結(jié)果表明，“切線弧+圓弧”外形設計最有利于優(yōu)化其脈動壓力環(huán)境。