張明星，蔡國飆，王永海

（1.北京航天長征飛行器研究所，北京100076；2.北京航空航天大學(xué)宇航學(xué)院，北京100191）

1 引言

姿軌控發(fā)動機(jī)在高空工作時，其噴流自由膨脹形成羽流［1-2］。羽流的氣動力、熱效應(yīng)會給航天器帶來干擾力、力矩和熱載荷等負(fù)面影響，輕則降低航天器的使用效率，重則使航天器無法正常工作，因此必須對發(fā)動機(jī)真空羽流進(jìn)行深入研究。

真空羽流包含了連續(xù)流、過渡流和自由分子流，流動狀態(tài)復(fù)雜，因此人們往往需要開展試驗(yàn)來研究發(fā)動機(jī)羽流。然而在很多情況下，限于經(jīng)費(fèi)、試驗(yàn)周期或設(shè)備能力等的制約，使用全尺寸模型開展試驗(yàn)是很困難的，甚至是無法實(shí)現(xiàn)的。此時就需要借助縮尺技術(shù)，利用縮尺技術(shù)“以小見大”的特點(diǎn)來解決這個問題。

縮尺技術(shù)的理論基礎(chǔ)是相似準(zhǔn)則，在過去幾十年里，研究人員對發(fā)動機(jī)流場的相似性開展了大量研究，并提出了很多相似準(zhǔn)則。早在1955年，S.S.Penner［3］針對液體火箭發(fā)動機(jī)燃燒室的縮尺技術(shù)進(jìn)行了理論研究，提出以Re、Sc、Pr數(shù)等為相似準(zhǔn)則，但為了滿足這些無量綱參數(shù)研究忽略了Ma數(shù)。隨后兩年，航空航天研究與發(fā)展咨詢組（AGARD）針對火箭發(fā)動機(jī)縮尺技術(shù)開展了大量的研究［4-7］，這些研究與S.S.Penner的研究方法類似，只是在某個物性參數(shù)的使用上進(jìn)行了改進(jìn)。 1957 年，S.S.Penner和 A.E.Fuhs［8］通過將液滴直徑與Re數(shù)和We數(shù)建立關(guān)系，對以前的相似準(zhǔn)則進(jìn)行了簡化。1963年，R.B.Lawhead等［9］提出模型發(fā)動機(jī)的燃燒室壓強(qiáng)和推進(jìn)劑噴注密度應(yīng)與原型發(fā)動機(jī)相同，若要研究發(fā)動機(jī)的破壞聲波模式，還應(yīng)保證頻率相同。隨后針對發(fā)動機(jī)的相似研究中斷了一段時間，直到2004年，C.E.Dexter等［10］對前人的工作進(jìn)行了總結(jié)，綜述了液體發(fā)動機(jī)的縮比技術(shù)。2006年，R.J.Kenny等［11］發(fā)現(xiàn)以前的縮尺技術(shù)沒有關(guān)注噴注面板，噴注器的尺寸并不隨模型尺寸變化，因此他用1∶2和1∶4的模型，研究了噴注面板尺寸對發(fā)動機(jī)縮尺技術(shù)的影響。 2008 年，J.R.Hulka［12］討論了液體發(fā)動機(jī)的縮尺技術(shù)。2010年，汪小衛(wèi)等［13-14］研究了氣氫氣氧發(fā)動機(jī)的相似準(zhǔn)則，研究從控制方程和π定理兩個角度出發(fā)，推出了適用于氣氫氣氧發(fā)動機(jī)的相似準(zhǔn)則。

雖然已經(jīng)有很多學(xué)者提出了相似準(zhǔn)則，但這些相似準(zhǔn)則尚有兩個不足：一是只針對連續(xù)流或只針對稀薄流，沒有一個針對全流域的研究；另一個是這些相似準(zhǔn)則往往只能保證部分流場參數(shù)的相似，比如只是保證氣動力的相似，而不能讓所有流場參數(shù)相似。這兩方面的原因限制了這些相似準(zhǔn)則在羽流試驗(yàn)研究上的應(yīng)用，需要提出一種新的相似準(zhǔn)則來滿足羽流研究的要求。

本文結(jié)合氣體物理屬性，從方程分析法和量綱分析法推導(dǎo)出一種適用于發(fā)動機(jī)羽流的相似準(zhǔn)則，并通過理論分析和仿真試驗(yàn)對該相似準(zhǔn)則進(jìn)行驗(yàn)證。

2 相似準(zhǔn)則理論推導(dǎo)

發(fā)動機(jī)羽流包含連續(xù)流、過渡流和自由分子流，其中連續(xù)流用Navier-Stokes（NS）方程來描述，過渡流和自由分子流用Boltzmann方程描述。由于流域描述方程不同，因此對羽流相似準(zhǔn)則的推導(dǎo)將包含兩部分，分別針對連續(xù)流域和稀薄流域。

一個完整的相似準(zhǔn)則包括幾何相似、運(yùn)動相似、動力相似、熱相似和邊界條件相似。幾何相似是指兩個流場對應(yīng)尺寸比都相同；運(yùn)動相似是指流體速度和速度梯度比在相似位置處相同；動力相似是指兩流場相似位置處力方面的比值相同；熱相似是指溫度和熱流的相似。幾何相似是容易實(shí)現(xiàn)的，因此這里主要關(guān)注其它三個相似要求。為了保證理論的正確，本文將采用方程分析法和無量綱分析法兩種方式來開展相似理論的推導(dǎo)。

2.1 連續(xù)流域相似理論

2.1.1 方程分析

對于真空羽流，氣體質(zhì)量力一般不予考慮，因此不考慮質(zhì)量力的無量綱Navier-Stokes方程可以寫成式（1）：

式中：p是壓強(qiáng)（單位：Pa）；t是時間 （單位：s）；cp是定壓比熱（單位：J／（kg·K））；K 是傳熱系數(shù)（單位：W／（m·k））；ρ是密度（單位：kg／m3）；μ是粘性系數(shù)（單位：Pa·s）；δij是 Kronecker張量。腳標(biāo)∞表示來流，?表示無量綱數(shù)。

根據(jù)動力學(xué)理論，粘性系數(shù)［15-16］可以寫作式（2）［17］：

式中：μ0是溫度T0時的粘性系數(shù)；ω是粘性溫度指數(shù)。

根據(jù) Chapman-Enskog 理論，有式（3）［17］：

式中：k是波爾茲曼常數(shù)（1.38 ×10－23J／k）；m是分子質(zhì)量（單位：kg）。

當(dāng)流體介質(zhì)可以被看作理想氣體時，即理想氣體公式 p＝ρkT／m 成立時，將公式（2）和（3）代入（1）中，可以得到式（4）：

若假設(shè)兩流場的cp值相同，當(dāng)幾何相似滿足且流體介質(zhì)相同時，NS方程中無量綱參數(shù)相同的條件是：

需要注意的是，上述推導(dǎo)過程中做了一個假設(shè)，即兩個流場相似位置處的cp值相同，因此我們需要明確這個假設(shè)何時才能成立。

對于理想的多原子氣體，cp是溫度的函數(shù)，可用式（5）表示［17］：

式中：ξr是轉(zhuǎn)動自由度，θvi是振動特征溫度。

由于發(fā)動機(jī)噴流總溫一般都較高，與振動特征溫度相當(dāng)，因此需要兩流場相似位置處溫度相同才能使cp值相同。

綜上所述，若要兩個氣流場嚴(yán)格相似，除了幾何相似和流體介質(zhì)相同外，還需要滿足以下條件：T1∞＝ T2∞， u1∞＝ u2∞且 L1ρ1∞＝ L2ρ2∞。

壁面邊界條件應(yīng)當(dāng)滿足以下條件以保證流場在近壁面區(qū)域的相似性：兩個流場壁面帶量綱參數(shù)的比值要與氣體流場相應(yīng)參數(shù)的比值相同，且兩流場壁面無量綱參數(shù)相等。例如，壁面溫度比（T1w／T2w）應(yīng)當(dāng)與流體來流溫度比（T1∞／T2∞）相同，而壁面熱適應(yīng)系數(shù)應(yīng)相同（α1＝α2）。

2.1.2 量綱分析

無量綱組合的提取分析已經(jīng)在Buckingham定理（π定理［18］）中標(biāo)準(zhǔn)化了，本文將采用這個理論來推導(dǎo)相似準(zhǔn)則。

熱對流是極為復(fù)雜的物理過程，幾乎所有的流場參數(shù)都會對其有影響，因此選擇熱對流作為研究目標(biāo)。與熱對流有關(guān)的參數(shù)包括：研究目標(biāo)尺寸l、單位力f、特征長度L、氣流速度u、氣流壓強(qiáng)p、氣流密度ρ、氣體溫度T、氣流粘性系數(shù)μ、氣流比熱 cp、熱傳導(dǎo)系數(shù) K、壁面溫度 Tw以及其它壁面條件。這里作者假設(shè)兩個流場相似位置處的比熱cp相同，同時其它壁面條件也都滿足相似條件，那么氣動熱流的計算公式可寫作式（6）：

q ＝G（l， f， L， u， ρ， p， T， μ， K， Tw） （6）

這些參數(shù)的基本單位包括：長度（m）、時間（s）、溫度（K）和質(zhì)量（kg）。 給定這些單位，本文選擇特征長度L、來流速度u∞、來流溫度T∞和來流密度ρ∞作為基本參數(shù)，可以得到8個無量綱參

公式（6）可寫作式（7）：

從公式（7）可以看出，只要 π2，…π8相同，π1就是不變的。當(dāng)不考慮質(zhì)量力，滿足以下條件可實(shí)現(xiàn)這個目標(biāo)：1）幾何相似；2）流體介質(zhì)相同；3）相同的 T∞／u2∞和 Lρ∞／Tω－0.5∞；4）壁面溫度比與流場溫度比相同?？紤]到發(fā)動機(jī)流場溫度很高，則相似條件為： T1∞＝ T2∞，u1∞＝ u2∞且L1ρ1∞＝L2ρ2∞?？梢钥闯鲞@個結(jié)論與前面方程分析法得到的結(jié)論是一致的。

就熱相關(guān)的問題而言，Normand［19］詳細(xì)描述了與該復(fù)雜物理過程相關(guān)的參數(shù)，這里我們使用這些參數(shù)來驗(yàn)證相似準(zhǔn)則的正確性。Normand認(rèn)為對熱傳導(dǎo)有重要影響的無量綱參數(shù)是馬赫數(shù)M、普朗特數(shù)Pr、雷諾數(shù) Re、格拉曉夫數(shù) Gr、瑞利數(shù)Ra、耗散數(shù)Di和努塞爾數(shù)Nu。當(dāng)以上這些無量綱數(shù)相同時，兩流場的氣動熱是相似的?，F(xiàn)在我們檢驗(yàn)前面推導(dǎo)的相似準(zhǔn)則能否使這些無量綱數(shù)相同。

根據(jù)推導(dǎo)的相似準(zhǔn)則，可得到以下方程：

對于自然對流，努塞爾數(shù)Nu是瑞利數(shù)Ra和普朗特數(shù)Pr的函數(shù)；而對于強(qiáng)迫對流，努塞爾數(shù)Nu是雷諾數(shù)Re和普朗特數(shù)Pr的函數(shù)。無論對哪一種對流，由于Ra、Re和Pr是相同的，兩流場的Nu也自動滿足相同。綜上所述，基于π定理推導(dǎo)得到的相似準(zhǔn)則是正確的。

2.2 稀薄流域相似理論

對于稀薄流的相似，沈青在《稀薄氣體動力學(xué)》一書中介紹了雙尺度律［17］，我們可以發(fā)現(xiàn)，雙尺度律在本質(zhì)上與連續(xù)流的相似要求是一致的，這說明只要上游的連續(xù)流相似，那么下游的稀薄流也是相似的。

3 相似準(zhǔn)則仿真驗(yàn)證

以我國航天器上使用的某種3000 N發(fā)動機(jī)為例，其推進(jìn)劑為一甲基肼／四氧化二氮（MMH／NTO），總壓0.8 MPa，燃?xì)鉁囟葹?040 K。選擇縮尺比例為0.2，根據(jù)相似準(zhǔn)則的要求，兩個發(fā)動機(jī)的參數(shù)設(shè)計如表1所示。

表1 原型和縮比發(fā)動機(jī)的相關(guān)參數(shù)Table 1 Parameters of the prototype and the scaling engines

根據(jù)相似準(zhǔn)則，兩流場參數(shù)理論上應(yīng)符合表2中的關(guān)系：

表2 兩發(fā)動機(jī)流場參數(shù)的理論關(guān)系Table 2 Theoretical parameter ratios of two similar flow fields

發(fā)動機(jī)羽流在真空中自由膨脹時會有明顯的膨脹波，但不一定有明顯的壓縮波。為了在流場中形成比較顯著的壓縮波，在距離發(fā)動機(jī)出口平面0.3 m的位置處放置一個圓盤，圓盤軸線與發(fā)動機(jī)軸線重合，如圖1所示?？s尺流場的圓盤尺寸和離發(fā)動機(jī)距離相應(yīng)調(diào)整為原流場的1／5。

發(fā)動機(jī)羽流包括了連續(xù)流和稀薄流，本文用CFD／DSMC弱耦合的方式對兩個流域分別進(jìn)行仿真，即先用CFD方法求解發(fā)動機(jī)內(nèi)流場及出口附近流場，然后選取CFD計算結(jié)果的某一截面參數(shù)作為DSMC方法計算的入口參數(shù)，最后用DSMC方法求解稀薄流域。

圖1 羽流場計算區(qū)域劃分示意圖Fig．1 Diagram of computation area division of flow field

3.1 連續(xù)流仿真

表3所示為燃?xì)獾臒崃τ嬎憬Y(jié)果，可以看出燃?xì)庵饕傻獨(dú)?、水蒸汽、一氧化碳、二氧化碳和氫氣組成，這里將噴流視為凍結(jié)流，那么燃?xì)獾钠骄镄詤?shù)為：平均摩爾質(zhì)量20.9 kg／mol、比熱1759 J／（kg·K）、導(dǎo)熱系數(shù)382 T0.75×e－6W／（m2·k）、粘性系數(shù) 0.256 T0.75× e－6Pa·s。壁面設(shè)為絕熱、無滑移，湍流模型選擇SST k－ω模型，通量格式使用Roe格式。由于燃?xì)饪倝翰桓?，計算中將氣體視作理想氣體。

表3 發(fā)動機(jī)燃?xì)鈪?shù)Table 3 Parameters of the engine plume

圖2 流場云圖對比Fig．2 Comparison of the two flow fields

流場最小網(wǎng)格的尺寸為0.002 mm，位于發(fā)動機(jī)喉部壁面附近；最大網(wǎng)格的尺寸為3 mm，位于自由邊界處；圓盤擋板附近使用等比網(wǎng)格進(jìn)行局部加密，網(wǎng)格最小尺寸為0.005 mm，在壁面附近。

將縮尺流場參數(shù)按表2所示的比例進(jìn)行相應(yīng)的縮放，并與原尺寸流場進(jìn)行對比，流場云圖對比如圖2所示。從圖中我們可以明顯地看出，經(jīng)過縮放后的縮尺流場計算結(jié)果與原流場結(jié)果高度一致：噴管內(nèi)外主流區(qū)流場形態(tài)基本相同，邊界層結(jié)構(gòu)基本相同，膨脹波和壓縮波的形態(tài)及位置基本相同。從云圖對比結(jié)果來看，縮尺流場與原流場是相似的。

需要注意的是，在發(fā)動機(jī)燃燒室附近區(qū)域兩流場的仿真結(jié)果有一定的差別（見圖2c、d），但這并不能說明兩流場不相似。因?yàn)閷?shí)際上這個區(qū)域是回流區(qū)，已經(jīng)屬于稀薄流區(qū)，不適合用CFD方法求解，因此該處的計算結(jié)果不作為分析依據(jù)。回流區(qū)的計算我們將在下一小節(jié)中介紹。

為了進(jìn)一步比較兩流場各參數(shù)，作兩流場的等高線圖，如圖3所示。從圖中可以清晰地看出，在發(fā)動機(jī)內(nèi)部及外部的連續(xù)流區(qū)，兩流場的等高線完全重合，只有在回流區(qū)有一些差異。表明本文的相似準(zhǔn)則在發(fā)動機(jī)連續(xù)流區(qū)是適用的。

圖3 流場各參數(shù)對比圖（連續(xù)流）Fig．3 Comparison of the contours （continuous flow field）

3.2 稀薄流仿真

在原流場中，以發(fā)動機(jī)出口邊緣為起點(diǎn)，用與發(fā)動機(jī)軸線60°夾角的錐截取流場，作為DSMC的入口參數(shù)；縮比的計算區(qū)域劃分與原流場相似，但對流場截面的參數(shù)進(jìn)行了稍許改動，使入口參數(shù)嚴(yán)格滿足相似關(guān)系。

使用北航自主研發(fā)的基于DSMC方法的羽流計算軟件PWS［20-22］對稀薄流域進(jìn)行仿真，在分子碰撞的模擬中選取了不計時間的方式，分子勢能的計算中選取了可變硬球模型（VHS），壁面設(shè)為粗糙壁面，壁面反射模型為漫反射模型。計算中采用5種組分的氣體，氣體種類分別為N2、H2O、CO、CO2、H2，其比例見表 3。 流場網(wǎng)格尺寸均小于當(dāng)?shù)胤肿幼杂沙痰囊话?，時間步長小于分子平均碰撞時間的1／3，每個網(wǎng)格內(nèi)粒子數(shù)不少于15。

將縮尺流場參數(shù)按表2所示的比例進(jìn)行相應(yīng)的縮放，并與原尺寸流場進(jìn)行對比，流場對比如圖4所示。可以看出兩流場各參數(shù)是基本一致的，說明對于稀薄流區(qū)，本文的相似準(zhǔn)則也是適用的。同時也證明了2.2中關(guān)于稀薄流相似推導(dǎo)的正確性。

4 結(jié)論

本文結(jié)合氣體物理屬性，從方程分析法和量綱分析法推導(dǎo)出了一種適用于發(fā)動機(jī)羽流的相似準(zhǔn)則，相似條件為條件：T1∞＝ T2∞，u1∞＝ u2∞且L1ρ1∞＝ L2ρ2∞。 隨后通過理論方法驗(yàn)證了相似準(zhǔn)則的正確性，該準(zhǔn)則能夠使全流域流場參數(shù)滿足相似關(guān)系。最后通過CFD／DSMC弱耦合的方法對所推導(dǎo)的相似準(zhǔn)則開展仿真驗(yàn)證。仿真結(jié)果與理論預(yù)測完全一致，進(jìn)一步驗(yàn)證了相似準(zhǔn)則的正確性。本文的相似準(zhǔn)則對于發(fā)動機(jī)羽流的試驗(yàn)研究具有一定的參考價值。

圖4 稀薄流區(qū)各參數(shù)對比圖Fig．4 Comparison of the contours （rarefied flow field）

參考文獻(xiàn)（References）

［1］ 蔡國飆.真空羽流問題的數(shù)值模擬研究［D］.北京：北京航空航天大學(xué)，1995.Cai G B.Numerical Investigation of the Vacuum Plume［D］.Beijing： Beihang University， 1995.（in Chinese）

［2］ 蔡國飆，祖國君，王慧玉，等.羽流問題研究概況［J］.推進(jìn)技術(shù)，1995，16（4）：50-54.Cai G B，Zu G J，Wang H Y，et al.The sketchy situation of study on plume ［J］.Journal of Propulsion Technology，1995，16（4）：50-54.（in Chinese）

［3］ Penner S S.Similarity analysis for chemical reactors and the scaling of liquid fuel rocket engines［R］.Combustion Research and Reviews，AGARD Publications by Butterworths Scientific Publication， 1955：1-41.

［4］ Weller A E.Similarities in combustion［C］／／Selected Combustion Problems，2nd AGARD Combustion Colloquium，Belgium，1956：371-383.

［ 5 ］ Stewart D G.Scaling of gas turbine combustion systems［C］／／Selected Combustion Problems，2nd AGARD Combustion Colloquium，Belgium，1956：384-413.

［6］ Crocco L.Considerations on the problem of scaling rocket engines［C］／／Selected Combustion Problems， 2nd AGARD Combustion Colloquium，Belgium，1956：457-468.

［7］ Ross C C.Scaling of liquid fuel rocket combustion chambers［C］／／Selected Combustion Problems， 2nd AGARD Combustion Colloquium，Belgium，1956：444-456.

［8］ Penner S S，F(xiàn)uhs A E.On generalized scaling procedures for liquid-fuel rocket engines［J］.Combustion and Flame， 1957，1（2）： 229-240.

［9］ Lawhead R B， Combs L P.Modeling techniques for liquid propellant rocket combustion processes［C］ ／／9th International Symposium on Combustion， The Combustion Institute，1963：973-981.

［10］ Yang V，Habiballah M，Hulka J，et al.Liquid Rocket Thrust Chambers： Aspects of Modeling， Analysis， and Design［M］.Washington DC：American Institute of Aeronautics＆ Astronautics， 2004： 553-600.

［11］ Kenny R J， Moser M D， Hulka J， et al.Cold flow testing for liquid propellant rocket injector scaling and throttling［C］／／42nd Aiaa／asme／sae／asee Joint Propulsion Conference ＆ Exhibit， Sacramento， California， 2006： 4705.

［12］ Hulka J R.Scaling of performance in liquid propellant rocket engine combustion devices［ C］ ／／44th Aiaa／asme／sae／asee Joint Propulsion Conference＆ Exhibit.Hartford，CT，2008：5113.

［13］ Wang X W， Cai G B， Jin P.Scaling of the flowfield in a combustion chamber with a gas-gas injector［ J］.Chinese Physics B，2010，19（1）：019401.

［14］ Wang X W， Jin P， Cai G B.Scaling of gas-gas cold flow fields［J］.Acta Aeronautica et AstronauticaSinica， 2010， 31（7）： 1305-1311.

［15］ Sun Q H， Boyd I D， Candler G V.A hybrid continuum／particle approach for modeling subsonic rarefied gas flows ［J］.Journal of Computational Physics， 2004， 194（1）： 256-277.

［16］ Schwartzentruber T E， Scalabrin L C，Boyd I D.A modular particle-continuum numerical method for hypersonic non-equilibrium gas flows ［J］.Journal of Computational Physics，2007， 225（1）： 1159-1174.

［17］ Shen C.Rarefied Gas Dynamics［M］.Germany： Springer-Verlag Berlin Heidelberg， 2005：51-129.

［18］ Buckingham E.On physically similar systems， illustrations of the use of dimensional equations ［J］.Physical Review，1914， 4（4）： 345-376.

［19］ Normand C， Pomeau Y， and Velarde M G.Convective instability： a physicist’s approach ［J］.Reviews of Modern Physics， 1977， 49（3）： 581.

［20］ He X Y， He B J， Cai G B.Simulation of rocket plume and lunar dust using DSMC method ［J］.Acta Astronautica，2012， 70（1）： 100-111.

［21］ He B J， Zhang J H， Cai G B.Research on vacuum plume and its effects［J］.Chinese Journal of Aeronautics， 2013， 26（1）： 27-36.

［22］ He B J， He X Y， Zhang M X， et al.Plume aerodynamic effects of cushion engine in lunar landing ［J］.Chinese Journal of Aeronautics， 2013， 26（2）： 269-278.