李凰立，蘇 虹，沈 丹

（北京宇航系統(tǒng)工程研究所，北京，100076）

0 引 言

脈動(dòng)壓力是一種隨時(shí)間變化的壓力，所以當(dāng)其作用于火箭上時(shí)，脈動(dòng)壓力是不均勻的，很可能對某個(gè)部位產(chǎn)生較集中或較大的壓力。同時(shí)，非定常流體因脈動(dòng)也會(huì)產(chǎn)生氣動(dòng)噪聲，在航空航天工程中，大量試驗(yàn)研究表明：脈動(dòng)壓力激起的氣動(dòng)噪聲的幅值一般在130～170 dB，頻率覆蓋了低、中、高頻。通常將5～100 Hz的振動(dòng)視為低頻段，將20～2000 Hz的振動(dòng)稱為高頻段（高、低頻不只是以頻率大小分類，還與振動(dòng)機(jī)理有關(guān)）?；鸺龤んw的頻率一般較低，低頻段的脈動(dòng)壓力往往會(huì)激起結(jié)構(gòu)的強(qiáng)烈振動(dòng)，造成設(shè)備破壞。非定常脈動(dòng)引起的噪聲在向外傳播的同時(shí)，也會(huì)傳入飛行器內(nèi)部，以寬帶噪聲的形式影響艙內(nèi)儀器或者工作人員。如果噪聲頻率與箭體特定模態(tài)的頻率范圍接近時(shí)，就有可能導(dǎo)致箭體結(jié)構(gòu)的疲勞甚至破壞，因此火箭結(jié)構(gòu)強(qiáng)度設(shè)計(jì)必須考慮脈動(dòng)壓力的影響。

飛行試驗(yàn)經(jīng)驗(yàn)表明，脈動(dòng)壓力在火箭跨聲速段至最大動(dòng)壓時(shí)間段較大，隨著火箭高度增加（空氣密度減?。┲饾u減小。因此，脈動(dòng)壓力研究的重點(diǎn)在跨聲速時(shí)段。但是，無論是開展風(fēng)洞試驗(yàn)?zāi)M還是工程預(yù)示、數(shù)值模擬，跨聲速脈動(dòng)壓力的研究難度都很大，尤其是數(shù)值模擬，對分析方法和數(shù)值算法都提出了較高的要求?？缏曀倜}動(dòng)壓力主要因強(qiáng)激波/邊界層相互干擾而產(chǎn)生，體現(xiàn)在激波位置的前后移動(dòng)和強(qiáng)激波后的非定常流動(dòng)分離等多方面。當(dāng)前學(xué)術(shù)界對脈動(dòng)壓力的非定常數(shù)值模擬大多為直接數(shù)值模擬（Direct Numerical Simulation，DNS），但是DNS在現(xiàn)有計(jì)算條件下只能分析小雷諾數(shù)情況下的流動(dòng)，對于高雷諾數(shù)DNS計(jì)算量計(jì)算時(shí)間過多；大渦模擬法（Large Eddy Simulation，LES）也是非定常模擬方法，它可精確計(jì)算大尺度非定常運(yùn)動(dòng)及分離流動(dòng)，然而模擬邊界層流動(dòng)時(shí)，LES所需計(jì)算資源也很大，且近壁模式不成熟，還需發(fā)展和完善；雷諾平均法（Reynolds Average Navier-Stokes，RANS）求解Reynolds平均的Navier-Stokes方程組，可以較好地模擬附體流動(dòng)，且對計(jì)算資源需求相對較低，但缺點(diǎn)是依賴于經(jīng)驗(yàn)參數(shù)，湍流信息少，高頻小尺度運(yùn)動(dòng)模擬不準(zhǔn)，即便獲得部分非定常效應(yīng)，其可信度也不高。

鑒于LES和RANS都存在各自的優(yōu)勢和局限性，因此若能有機(jī)結(jié)合兩種或多種方法的優(yōu)點(diǎn)并利用自身的優(yōu)點(diǎn)克服對方的不足，就有可能實(shí)現(xiàn)計(jì)算精度和效率的統(tǒng)一。RANS/LES的混合方法即具備上述特征。本文主要采用了RANS/LES方法分析脈動(dòng)壓力，實(shí)踐證明算法有效、數(shù)據(jù)可信、規(guī)律合理。

1 計(jì)算方法

為了準(zhǔn)確并且比較快速地預(yù)測跨聲速飛行器的壓力脈動(dòng)，如下的關(guān)鍵技術(shù)需要首先解決：預(yù)測非定常特性的時(shí)間推進(jìn)方法，預(yù)測非定常湍流方法。

1.1 預(yù)測非定常流動(dòng)的高精度非定常時(shí)間推進(jìn)方法

首先采用有限體積方法，按照某種空間格式離散Navier-Stokes方程組，將其轉(zhuǎn)化成為每個(gè)網(wǎng)格單元上關(guān)于時(shí)間的一階常微分方程組。對常微分方程組進(jìn)行時(shí)間推進(jìn)求解時(shí)，一般采用顯式格式或隱式格式。

顯式時(shí)間推進(jìn)簡單，易于實(shí)現(xiàn)，占用內(nèi)存較少，工程中較為實(shí)用且應(yīng)用廣泛。但是顯式推進(jìn)格式有其自身固有的缺陷：受穩(wěn)定性限制，CFL數(shù)不可取得太大，時(shí)間步長較小。在模擬非定常流動(dòng)時(shí)，往往關(guān)心流動(dòng)隨時(shí)間的變化過程，如果使用顯式方法，時(shí)間步長可能取得很小，導(dǎo)致計(jì)算效率降低。

相比之下隱式方法可取較大的全流場統(tǒng)一時(shí)間步長，計(jì)算效率較高，常用的隱式方法有LU-SGS分解方法，基本思想是將塊對角矩陣分解為上、下兩個(gè)三角矩陣，避免了繁雜的矩陣運(yùn)算，節(jié)約了計(jì)算時(shí)間，可用于定常或非定常流動(dòng)計(jì)算。但是 LU-SGS分解方法是一階時(shí)間精度，為了充分發(fā)揮LU-SGS隱式方法的優(yōu)勢，提高模擬非定常流動(dòng)隱式推進(jìn)過程的時(shí)間精度，現(xiàn)采用偽時(shí)間子迭代技術(shù)（Pseudo Time Sub-iteration，τTS），也常被稱作雙時(shí)間推進(jìn)技術(shù)。通過在原控制方程中引入偽時(shí)間導(dǎo)數(shù)項(xiàng)，即采用雙時(shí)間方法，借助偽時(shí)間方向的“子迭代”技術(shù)，對LU-SGS隱式方法進(jìn)行改進(jìn)，使其達(dá)到二階時(shí)間精度（本文簡稱為 LU-SGS-τTS）。

有限體積法離散后的N-S方程式為

式中 Q為守恒變量向量；S為源項(xiàng)，如無源項(xiàng)則S=0；I·n為凈通量，且，

式中e·In為無粘通量；ν·In為粘性通量。

其半離散格式為

時(shí)間方向采用向后差分離散，得到：

當(dāng)φ=0.5時(shí)，時(shí)間方向上可以實(shí)現(xiàn)二階精度。

引入了虛擬時(shí)間項(xiàng)之后可以計(jì)算非定常流場。下文即采用二階精度 LU-SGS-τTS方法，根據(jù)經(jīng)驗(yàn)、計(jì)算精度和計(jì)算效率，子迭代步數(shù)選取為3步。

1.2 預(yù)測非定常流動(dòng)的湍流模擬方法

RANS在時(shí)間域上對流場物理量進(jìn)行雷諾平均化處理，然后求解所得到的時(shí)均化控制方程。比較常用的模型包括S-A模型、k-ω模型、k-ε模型。雷諾時(shí)均模擬方法可以較好地模擬附體流動(dòng)，精度基本可以滿足工程需要，且計(jì)算效率高，對計(jì)算資源的需求相對較低；但RANS的缺點(diǎn)是只能提供湍流的平均信息，對于火箭關(guān)鍵部位的非定常脈動(dòng)壓力預(yù)測，給出的結(jié)果不夠準(zhǔn)確，不能很好地捕捉跨聲速激波/邊界層相互干擾的非定常特性。

LES屬于湍流計(jì)算方法中的尺度解析模擬方法，主要是對流場中一部分湍流進(jìn)行直接求解，其余部分通過數(shù)學(xué)模型來計(jì)算。LES方法在求解瞬態(tài)性和分離性比較強(qiáng)的流動(dòng)中，相對RANS具有優(yōu)勢，LES可以準(zhǔn)確地模擬分離流動(dòng)；但LES方法對流場計(jì)算網(wǎng)格要求較高，特別是近壁面區(qū)域網(wǎng)格密度要求大于 RANS方法，在計(jì)算邊界層流動(dòng)時(shí)卻需要耗費(fèi)極大的計(jì)算資源。

RANS/LES混合方法結(jié)合RANS、LES兩種或多種方法的優(yōu)點(diǎn)并利用自身的優(yōu)點(diǎn)克服對方的不足，可有效實(shí)現(xiàn)計(jì)算精度和效率的統(tǒng)一，即RANS用湍流模式高效且比較可靠地模擬高頻小尺度運(yùn)動(dòng)占主導(dǎo)地位的近壁區(qū)域，LES方法可比較準(zhǔn)確地計(jì)算低頻大尺度運(yùn)動(dòng)占優(yōu)的非定常分離流動(dòng)區(qū)域。采用RANS/LES混合方法是當(dāng)前計(jì)算資源有限情況下的合理選擇，從工程應(yīng)用角度分析，它也被認(rèn)為是結(jié)合LES方法處理高雷諾時(shí)數(shù)下非定常流動(dòng)的一種切實(shí)可行方法。

此外，在湍流模型上，本文采用了SST k-ω模型。SST k-ω模型混合了k-ω模型和k-ε模型的優(yōu)勢，在近壁面處使用k-ω模型，而在邊界層外采用k-ε模型。SST k-ω模型包含了修正的湍流粘性公式，考慮了湍流剪切應(yīng)力的效應(yīng)，一般能更精確地模擬逆壓梯度引起的分離點(diǎn)和分離區(qū)大小。

本文研究的主要途徑就是采用基于 SST k-ω模式的RANS/LES方法，數(shù)值分析了跨聲速飛行器頭部區(qū)域流動(dòng)中產(chǎn)生的非定常激波/邊界層相互干擾，包括流動(dòng)分離、激波位置振動(dòng)等問題。

2 計(jì)算模型

以圖 1所示的火箭頭部為例，該外形頭部采用結(jié)構(gòu)化網(wǎng)格，網(wǎng)格形式見圖2，采用的是O - C型結(jié)構(gòu)網(wǎng)格。該網(wǎng)格的優(yōu)點(diǎn)在于此種拓?fù)浣Y(jié)構(gòu)既保證了物面附近較高的網(wǎng)格密度，又可以使整個(gè)流場空間內(nèi)保持合理而又均勻的網(wǎng)格分布。該網(wǎng)格體流向布置 330～350網(wǎng)格點(diǎn)，在外形變化劇烈的位置加密網(wǎng)格點(diǎn)；法向布置90～120網(wǎng)格點(diǎn)，貼近壁面處加密到10-3m量級；周向120～150網(wǎng)格點(diǎn)，背風(fēng)面較密。據(jù)此計(jì)算，總網(wǎng)格點(diǎn)數(shù)目大約在520萬左右。

圖1 火箭頭部外形示意Fig.1 Launch Vehicle Fairing Contour

圖2 計(jì)算網(wǎng)格示意Fig.2 Grid Topology On Fairing

另外，為了便于下文將計(jì)算數(shù)據(jù)和試驗(yàn)數(shù)據(jù)對比，此處附上風(fēng)洞試驗(yàn)的外形設(shè)計(jì)，見圖3。模型設(shè)計(jì)中，在火箭頭部柱段和倒錐段各選取一個(gè)測量截面。

圖3 風(fēng)洞模型測點(diǎn)和截面位置示意Fig.3 Wind-tunnel Model and the Station of the Fairing

3 計(jì)算狀態(tài)和計(jì)算結(jié)果

以工程經(jīng)驗(yàn)判斷，火箭在飛行過程中強(qiáng)烈的脈動(dòng)壓力環(huán)境出現(xiàn)在跨聲速階段。對圖1所示火箭頭部外形，開展了若干跨聲速狀態(tài)下的脈動(dòng)壓力預(yù)示。表 1列出幾個(gè)典型的馬赫數(shù)和攻角組合狀態(tài)。表1中同時(shí)還說明了對應(yīng)狀態(tài)下的流動(dòng)特征。

表1 典型計(jì)算狀態(tài)Tab.1 Typical States of Computation

下文分析各狀態(tài)計(jì)算結(jié)果，其中脈動(dòng)壓力結(jié)果均已經(jīng)過特定參考值無量綱化處理，包括后文試驗(yàn)數(shù)據(jù)。

本狀態(tài)計(jì)算所得最大均方根脈動(dòng)壓力系數(shù)大致在1.5%左右，且其位置在流向截面Ⅰ附近，與真實(shí)流動(dòng)物理特性相符。脈動(dòng)壓力系數(shù)反映了脈動(dòng)壓力的強(qiáng)弱。圖 4給出瞬時(shí)流場圖，分別是表面流線、物面壓力分布、速度梯度二階不變量Q等值面，可見流動(dòng)具有強(qiáng)非定常特性，尤其是通過湍動(dòng)能變化，可看出最大脈動(dòng)壓力區(qū)域明顯在截面Ⅰ之后。

該狀態(tài)攻角同樣是 0°，但相比狀態(tài)Ma=0.70，α=0°，馬赫數(shù)為0.95更接近1，此時(shí)發(fā)生大振幅脈動(dòng)壓力的位置在倒錐段即流向截面Ⅱ之后，見圖5。

圖4 Ma=0.70，α=0°流場圖Fig.4 Flow Field of Ma=0.70，α=0°

圖5 Ma=0.95，α=0°流場圖Fig.5 Flow Field of Ma=0.95，α=0°

通過壓力等值線圖可以看出，頭部前錐-柱段的激波位置往復(fù)移動(dòng)，因此作用在火箭頭部的力也來回變化，它與火箭結(jié)構(gòu)相互作用會(huì)產(chǎn)生結(jié)構(gòu)噪聲；同時(shí)會(huì)形成非定常的俯仰力矩，導(dǎo)致火箭壓心位置變化，影響火箭的穩(wěn)定性和操縱性。

另外，圖5顯示前錐-柱段流動(dòng)變化不大，存在很多類似小波紋的渦流結(jié)構(gòu)，但是在倒錐段的湍動(dòng)能則變化劇烈。在倒錐段，由于收縮外形導(dǎo)致流動(dòng)擴(kuò)張產(chǎn)生更強(qiáng)的激波/邊界層干擾，強(qiáng)激波導(dǎo)致流動(dòng)大范圍分離，所以該部位的壓力脈動(dòng)最大。經(jīng)統(tǒng)計(jì)了解到，截面Ⅰ均方根脈動(dòng)壓力系數(shù)大致為0.3%～0.5%，截面Ⅱ部分測點(diǎn)均方根脈動(dòng)壓力系數(shù)達(dá)到1.6%以上，這說明在接近跨聲速M(fèi)a=1時(shí)，脈動(dòng)壓力較高能量已經(jīng)后移至箭體柱段-倒錐段（截面Ⅱ之后）。

相比狀態(tài) Ma=0.70，α=0°，此狀態(tài)有4°攻角，火箭整流罩頭部的流動(dòng)會(huì)出現(xiàn)上下不對稱現(xiàn)象，即上下對稱面的激波/邊界層干擾會(huì)不同，見圖6。

通過壓力脈動(dòng)幅值可見，背風(fēng)面的漩渦比較集中，背風(fēng)面分離區(qū)域較大，湍動(dòng)能占優(yōu)；但是激波強(qiáng)度并不大；因?yàn)楣ソ窃斐闪肆鲃?dòng)在背風(fēng)面分離，所以從截面Ⅰ開始，脈動(dòng)壓力就明顯增大，經(jīng)統(tǒng)計(jì)，截面Ⅰ最大均方根脈動(dòng)壓力系數(shù)約為1.3%；但是到截面Ⅱ之后，脈動(dòng)壓力的能量已經(jīng)大幅減弱，均方根脈動(dòng)壓力系數(shù)逐漸變小，截面Ⅱ之后均方根脈動(dòng)壓力系數(shù)大致只能達(dá)到0.2%。

圖6 Ma=0.71，α=4°流場圖Fig.6 Flow Field of Ma=0.71，α=4°

該狀態(tài)Ma=0.96，和狀態(tài)Ma=0.95接近，因此流場特性差異主要由攻角大小差異決定。該狀態(tài)流場如圖7所示。

對比圖 5（Ma=0.95，α=0°）、圖 7（Ma=0.96，α=4°）流線特征，可知有攻角時(shí)，流動(dòng)分離會(huì)發(fā)生在火箭整流罩頭部較前區(qū)域；而對比之前計(jì)算狀態(tài)，0°攻角時(shí)流動(dòng)分離會(huì)發(fā)生在倒錐段之后。

圖7 Ma=0.96，α=4°流場圖Fig.7 Flow Field of Ma=0.96，α=4°

狀態(tài) Ma=0.96，α=4°的最大均方根脈動(dòng)壓力系數(shù)出現(xiàn)在截面Ⅰ附近；對比之前狀態(tài) Ma=0.95，α=0°最大均方根脈動(dòng)壓力系數(shù)出現(xiàn)在截面Ⅱ附近。

另外，由于存在攻角，所以狀態(tài)Ma=0.96，α=4°在整個(gè)背風(fēng)面的區(qū)域內(nèi)，脈動(dòng)壓力的作用都比較明顯。

圖8 Ma=0.72，α=6°流場圖Fig.8 Flow Field of Ma=0.72，α=6°

本狀態(tài)攻角已經(jīng)達(dá)到 6°，相對 0°攻角、4°攻角，本狀態(tài)的迎風(fēng)面激波更強(qiáng)、背風(fēng)面出現(xiàn)漩渦。當(dāng)攻角為 6°時(shí)，漩渦已經(jīng)清晰可辨，從圖 8的流線分布中可以觀察到一個(gè)漩渦中心。另外，攻角6°時(shí)，湍動(dòng)能在漩渦區(qū)域內(nèi)的形態(tài)也更加豐富，和流動(dòng)速度相關(guān)的Q值也充分反映了流動(dòng)中漩渦的發(fā)展和演化。

該狀態(tài)脈動(dòng)壓力較大區(qū)域產(chǎn)生在流向截面Ⅰ處，此時(shí)最大均方根脈動(dòng)壓力系數(shù)約為1.3%，基本在迎風(fēng)面和背風(fēng)面交界處。對稱背風(fēng)面上壓力樣本點(diǎn)的均方根脈動(dòng)壓力系數(shù)值大約為0.5%，由此體現(xiàn)攻角對壓力的影響，反映了激波/邊界層/漩渦的相互作用。

該狀態(tài)可以和狀態(tài)Ma=0.72，α=6°進(jìn)行不同馬赫數(shù)、相同攻角的對比；可以和狀態(tài)Ma=0.95，α=0°、狀態(tài)Ma=0.96，α=4°進(jìn)行相同（相近）馬赫數(shù)、不同攻角的對比。流場圖參考圖9。

相對狀態(tài)Ma=0.72，α=6°，本狀態(tài)均方根脈動(dòng)壓力最大的截面是截面Ⅱ，尤其是在背風(fēng)面，部分點(diǎn)均方根脈動(dòng)壓力系數(shù)達(dá)到0.7%，較大的壓力脈動(dòng)來源于激波邊界層相互作用；相對狀態(tài)Ma=0.95，α=0°、狀態(tài) Ma=0.96，α=4°，本狀態(tài)攻角6°背風(fēng)面的旋渦比、4°情況更加集中，背風(fēng)面分離區(qū)域范圍更遠(yuǎn)，影響面積也更大；同時(shí)背風(fēng)面湍動(dòng)能占優(yōu)。在壓力場的圖中可以看到，截面Ⅱ之后，伴隨著較大的攻角，在背風(fēng)面擴(kuò)張角較大的地方容易產(chǎn)生空化、渦脫落現(xiàn)象，最終導(dǎo)致了背風(fēng)面壓力場的不穩(wěn)定，并且使這一區(qū)域保持了較大的壓力脈動(dòng)。

圖9 Ma=0.96，α=6°流場圖Fig.9 Flow Field of Ma=0.96，α=6°

4 數(shù)值預(yù)示果和風(fēng)洞試驗(yàn)的對比

風(fēng)洞試驗(yàn)和數(shù)值計(jì)算相輔相成，同時(shí)也是驗(yàn)證計(jì)算方法可靠性的重要手段。針對上述計(jì)算外形開展風(fēng)洞實(shí)驗(yàn)，下文將數(shù)值計(jì)算統(tǒng)計(jì)的均方根脈動(dòng)壓力系數(shù)峰值結(jié)果和試驗(yàn)結(jié)果進(jìn)行了比較，見表2。經(jīng)比較可知：0°攻角時(shí)計(jì)算結(jié)果和試驗(yàn)結(jié)果符合較好，有攻角時(shí)的計(jì)算結(jié)果與試驗(yàn)結(jié)果差異比較明顯；這說明流動(dòng)分離過程復(fù)雜，模擬比較困難，目前建議分離區(qū)的計(jì)算需要在背風(fēng)面加密網(wǎng)格，改進(jìn)背風(fēng)面近壁網(wǎng)格質(zhì)量，提高數(shù)值算法的精度，以辨識(shí)較大脈動(dòng)壓力信號。

表2 數(shù)值計(jì)算結(jié)果與風(fēng)洞試驗(yàn)數(shù)據(jù)對比Tab.2 Computation Results vs Wind-tunnel Data

5 結(jié)束語

針對跨聲速脈動(dòng)壓力提出了一種非定常計(jì)算方法，含雙時(shí)間推進(jìn)技術(shù)和RANS/LES混合湍流模擬方法。通過幾種典型狀態(tài)的數(shù)值模擬，了解到：外形變化劇烈的表面更容易引起較大的脈動(dòng)壓力；這些變形劇烈的區(qū)域或位置在火箭高速上升的同時(shí)，將引起強(qiáng)烈的流動(dòng)分離，在一定的馬赫數(shù)和攻角的狀態(tài)下，部分?jǐn)U張段在下游流場將產(chǎn)生或強(qiáng)或弱的激波，且隨著時(shí)間的推移前后振動(dòng)，從而產(chǎn)生或大或小的脈動(dòng)壓力。

通過上述計(jì)算結(jié)論可知，當(dāng)火箭以跨聲速飛行時(shí)，流動(dòng)是極不穩(wěn)定的；尤其是在火箭頭部外形變化的部段，例如錐-柱、柱-倒錐等外形變化部段，流場在強(qiáng)激波/邊界層相互作用下，流動(dòng)極易分離，分離流再引發(fā)激波前后不規(guī)則移動(dòng)，由此將形成非定常的激振力，嚴(yán)重時(shí)誘發(fā)箭體的非定常結(jié)構(gòu)響應(yīng)即抖振效應(yīng)。在馬赫數(shù)從亞跨聲速增加至跨聲速的過程中，抖振會(huì)越加明顯，當(dāng)攻角不為0°時(shí)，抖振的影響還將更加復(fù)雜。抖振存在諸多負(fù)面影響，如影響火箭的氣動(dòng)、穩(wěn)定和操縱性能、加速火箭結(jié)構(gòu)疲勞，嚴(yán)重時(shí)會(huì)引起全箭結(jié)構(gòu)振動(dòng)，還有前文提及的氣動(dòng)噪聲等不利影響，這些都應(yīng)該引起火箭設(shè)計(jì)者的重視，在研制階段需要落實(shí)相應(yīng)的解決方案。

最后需要說明的是，本文對跨聲速脈動(dòng)壓力的數(shù)值研究方法，還有一些待改進(jìn)的空間：a）當(dāng)開展有攻角計(jì)算時(shí)，背風(fēng)面的旋渦比較集中，背風(fēng)面的分離區(qū)域更大，因此背風(fēng)面需要加密背風(fēng)面的網(wǎng)格數(shù)量、改進(jìn)近壁面網(wǎng)格質(zhì)量；b）需要從湍流數(shù)值算法、網(wǎng)格設(shè)計(jì)方法等方面進(jìn)行改進(jìn)；c）結(jié)合相關(guān)外形的跨聲速脈動(dòng)壓力風(fēng)洞試驗(yàn)、或者飛行試驗(yàn)數(shù)據(jù)做出進(jìn)一步的修訂。