李凰立，蘇 虹，沈 丹

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

0 引 言

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

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

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

1 計算方法

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

1.1 預測非定常流動的高精度非定常時間推進方法

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

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

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

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

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

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

其半離散格式為

時間方向采用向后差分離散，得到：

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

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

1.2 預測非定常流動的湍流模擬方法

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

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

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

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

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

2 計算模型

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

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

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

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

圖3 風洞模型測點和截面位置示意Fig.3 Wind-tunnel Model and the Station of the Fairing

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

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

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

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

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

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

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

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

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

通過壓力脈動幅值可見，背風面的漩渦比較集中，背風面分離區(qū)域較大，湍動能占優(yōu)；但是激波強度并不大；因為攻角造成了流動在背風面分離，所以從截面Ⅰ開始，脈動壓力就明顯增大，經(jīng)統(tǒng)計，截面Ⅰ最大均方根脈動壓力系數(shù)約為1.3%；但是到截面Ⅱ之后，脈動壓力的能量已經(jīng)大幅減弱，均方根脈動壓力系數(shù)逐漸變小，截面Ⅱ之后均方根脈動壓力系數(shù)大致只能達到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°）流線特征，可知有攻角時，流動分離會發(fā)生在火箭整流罩頭部較前區(qū)域；而對比之前計算狀態(tài)，0°攻角時流動分離會發(fā)生在倒錐段之后。

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

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

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

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

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

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

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

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

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

4 數(shù)值預示果和風洞試驗的對比

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

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

5 結(jié)束語

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

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

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