運(yùn)載火箭推進(jìn)劑復(fù)雜流動傳熱問題數(shù)值模擬中的模型簡化方法

，,，，

(1.中國科學(xué)院國家空間科學(xué)中心, 北京 100190；2.中國科學(xué)院大學(xué), 北京 100190；3.北京宇航系統(tǒng)工程研究所，北京 100076)

0 引言

在運(yùn)載火箭動力系統(tǒng)總體設(shè)計中，越來越多地采用數(shù)值模擬方法。通過數(shù)值模擬，不僅可以對動力系統(tǒng)設(shè)計參數(shù)進(jìn)行優(yōu)化，還可以與理論設(shè)計、試驗(yàn)測量進(jìn)行交叉對比驗(yàn)證，提高系統(tǒng)的可靠性[1]。由于運(yùn)載火箭中的流動傳熱屬于經(jīng)典的流體力學(xué)和傳熱問題，相關(guān)的數(shù)值計算方法比較成熟，采用現(xiàn)成的商業(yè)軟件，就可以進(jìn)行計算，因而得到廣泛使用。工程型號設(shè)計人員采用數(shù)值模擬方法，解決了大量動力系統(tǒng)設(shè)計中的問題[2-4]。

但是，在運(yùn)載火箭改進(jìn)優(yōu)化或者是新型運(yùn)載火箭設(shè)計過程中，往往會遇到各種復(fù)雜的流動與傳熱問題。直接建立的數(shù)值計算模型過于復(fù)雜，多種不同材料和結(jié)構(gòu)耦合，同時出現(xiàn)在需要解決的問題之中。經(jīng)常出現(xiàn)計算時間太長(數(shù)年)、收斂困難、程序調(diào)試?yán)щy等問題。此時對數(shù)值計算模型進(jìn)行簡化，縮短計算時間的同時，保證結(jié)果的正確性、實(shí)用性，就是一個關(guān)鍵問題。美國Space X公司的Falcon 1火箭在初始的3次試射中均失敗，教訓(xùn)慘重。從公布的資料分析，F(xiàn)alcon 1號在設(shè)計上的經(jīng)驗(yàn)教訓(xùn)之一是流體管理方面，推進(jìn)劑貯箱沒有安裝防晃擋板。這可能是前兩次發(fā)射失敗的主要原因。之所以不安裝防晃擋板，是由于該公司數(shù)值模擬研究結(jié)果表明沒有必要安裝[1]。發(fā)射失敗后，他們找到了計算給出錯誤結(jié)果的原因是沒有考慮分離過程中的干擾，數(shù)值計算的初始條件設(shè)置錯誤。從第3次發(fā)射開始，增加了貯箱防晃擋板，克服了推進(jìn)劑晃動問題。由此可見，推進(jìn)劑晃動的數(shù)值模擬工作需要相關(guān)研究經(jīng)驗(yàn)作為基礎(chǔ)，模型簡化需要慎重。否則的話，很可能數(shù)值模擬得到與事實(shí)相反的錯誤結(jié)論，引起重大的設(shè)計錯誤。

早期的貯箱自生增壓設(shè)計，都采用集總參數(shù)法，大大簡化了計算。但缺點(diǎn)也十分明顯，這種方法精度非常低，也無法預(yù)測飛行過程中具體的壓力變化過程[5]。Ariane 5火箭上面級滑行期間壓力計算和流動計算分別進(jìn)行，使用了不同的計算模型[6]。其中，壓力變化采用集總參數(shù)法進(jìn)行計算，推進(jìn)劑分布狀態(tài)和流動則使用FLOW3D計算。這種簡化方法避開了計算量大的難點(diǎn)，但是壓力計算精度不夠，也無法應(yīng)用到更加復(fù)雜、要求更高的場合。

針對初始條件復(fù)雜、飛行過程復(fù)雜、邊界條件復(fù)雜3類具體問題，結(jié)合我們十幾年的研究經(jīng)驗(yàn)[7-8]，介紹3種相應(yīng)的模型簡化方法——工程經(jīng)驗(yàn)法、極限參數(shù)法、低維近似法，為推進(jìn)劑復(fù)雜流動與傳熱的數(shù)值模擬提供參考。具體的數(shù)值模擬計算方法，可以參考Kassemi等[9]的研究。

1 工程經(jīng)驗(yàn)法

液體火箭在飛行過程中，上面級可能關(guān)閉主發(fā)動機(jī)作空間慣性飛行(簡稱微重力滑行)，一段時間之后主發(fā)動機(jī)二次啟動。通過微重力滑行，可以提高運(yùn)載能力。為了保證主發(fā)動機(jī)二次啟動成功，一般在微重力滑行期間使用較小的沉底發(fā)動機(jī)工作，讓液體推進(jìn)劑始終保持在沉底位置[1]。我國的空間探測工程要求優(yōu)化微重力滑行段工作狀況,即延長微重力滑行時間，同時盡可能減小沉底發(fā)動機(jī)推進(jìn)劑用量。

針對微重力滑行后主發(fā)動機(jī)二次啟動問題，需要通過數(shù)值模擬計算微重力滑行期間液體推進(jìn)劑在貯箱中的分布、氣液界面的運(yùn)動情況。

推進(jìn)劑流動狀態(tài)高度依賴于初始條件。飛行過程中運(yùn)載火箭的振動、姿態(tài)控制干擾、主發(fā)動機(jī)瞬間的關(guān)機(jī)過程都會影響該時刻的流動狀態(tài)，該初始流動狀態(tài)的選擇和設(shè)置就成為一個關(guān)鍵問題。由于這些因素具有很強(qiáng)的不確定性，而且依賴于具體的工程型號，根據(jù)工程經(jīng)驗(yàn)選擇計算的初始條件是一個好的選擇。參考工程遙測數(shù)據(jù)和Berglund等[1]的研究，在某火箭的設(shè)計過程中，選用側(cè)向干擾力主推力10%大小作為初始干擾簡化計算條件。計算得到推進(jìn)劑晃動幅值為貯箱半徑10%的初始干擾運(yùn)動，作為微重力滑行期間流動計算的初始條件。

如圖1所示，初始干擾在軸向過載作用下形成近似一階晃動，晃動幅值為0.15m(最高點(diǎn)和最低點(diǎn)落差為兩倍幅值，即0.30m，貯箱直徑為3m)，t1時刻液體大部分區(qū)域速度接近于0，即流體近似處于動能極小，勢能極大的狀態(tài)。其左右兩側(cè)自由面附近液體位于極低和極高位置，t1時刻即最大勢能時刻。在t2時刻，自由面附近液體豎直方向速度分量達(dá)到最大值，而且它們分別位于箱體左右兩側(cè)。此時波動的幅值最小，但動能最大；這時刻如果主發(fā)動機(jī)關(guān)機(jī)，在慣性作用下將于下一時段形成最大晃幅的波動，t2時刻為最大動能時刻。

(a) t1 (b) t2

(c) t3 (d) t4圖1 液氫貯箱主發(fā)動機(jī)關(guān)機(jī)前液面附近的晃動干擾簡化，時間t1、t2、t3、t4分別是一個晃動周期內(nèi)的4個不同時刻(深色為液氫，淺色為增壓氣體，矢量表示中心截面上流體的速度)Fig.1 Simplification of sloshing near the interface, in the hydrogen tank, before main engine cut-off, t1,t2,t3, and t4 are four separate moments in a sloshing cycle

通過工程經(jīng)驗(yàn)法，快速得到了后續(xù)計算所需要的初始條件，省去了上面級火箭一次飛行段的直接三維非定常數(shù)值模擬計算。通常情況下，省去的部分采用10個CPU核計算，需要計算1m左右。

2 極限參數(shù)(狀態(tài))法

如前所述，主發(fā)動機(jī)關(guān)機(jī)過程也給推進(jìn)劑流動帶來了不確定性和復(fù)雜性。發(fā)動機(jī)關(guān)機(jī)過程中，1s～2s內(nèi)，上面級過載從幾個g快速過渡到10-3g的微重力狀態(tài)。進(jìn)入微重力時刻的流動狀態(tài)是未知的，可以是前面所述晃動周期中的任意時刻，具有不確定性。此時，可以用極限參數(shù)(狀態(tài))法進(jìn)行簡化，分別有最大動能關(guān)機(jī)、最大勢能關(guān)機(jī)兩種極限狀態(tài)。在保證主發(fā)動二次啟動成功問題中，最大勢能關(guān)機(jī)為最危險工況。此時推進(jìn)劑在微重力滑行中運(yùn)動最為劇烈，沉底困難。選用最大勢能時刻作為開始微重力滑行數(shù)值模擬計算的初始時刻進(jìn)行簡化，滿足主發(fā)動二次啟動設(shè)計要求。

盡管如此，由于微重力滑行期間流動的復(fù)雜性，工程中通常需要進(jìn)行保守設(shè)計，此時可以按照可能的最危險狀態(tài)，進(jìn)行沉底方案設(shè)計。假設(shè)經(jīng)過長時間微重力滑行后，推進(jìn)劑全部位于貯箱頂部，如圖2所示，然后施加沉底力作用。通過三維非定常直接數(shù)值模擬，得到保守的沉底力大小和沉底時間。

圖2 液氫貯箱推進(jìn)劑沉底前的極限狀態(tài)，液氫全部位于貯箱頂部Fig.2 Extreme status of the hydrogen propellant in the tank, liquid hydrogen all occupying the upper side of the propellant tank

3 低維近似法

上面級液氫貯箱在微重力滑行期間，在多種因素聯(lián)合作用下，貯箱壓力逐漸下降。引起貯箱壓力下降的因素包括氣枕溫度/密度重新分層，氣液界面?zhèn)鳠?、相變，貯箱壁面?zhèn)鳠帷Ｓ捎谕七M(jìn)劑干擾流動、非對稱壁面?zhèn)鳠?、非軸對稱調(diào)姿干擾的存在，貯箱壓力對這些因素都很敏感，只能進(jìn)行三維非定常直接數(shù)值模擬計算，才能夠獲得準(zhǔn)確的壓力下降預(yù)測，從而指導(dǎo)增壓輸送系統(tǒng)的優(yōu)化設(shè)計。

最直接的三維非定常數(shù)值計算模型，幾何模型上包括推進(jìn)劑、金屬貯箱、貯箱外表面的泡沫隔熱層，時間上從上面級發(fā)動機(jī)一次工作段開始計算數(shù)百秒，然后主發(fā)動機(jī)關(guān)機(jī)，進(jìn)入微重力滑行段。上述完整的幾何和時間模型會遇到多方面的困難。1)幾何模型復(fù)雜，推進(jìn)劑、金屬貯箱、泡沫分別是m、mm、cm 這3個尺度，網(wǎng)格數(shù)量大，網(wǎng)格質(zhì)量差，程序很難收斂。2)即使程序收斂，計算時間也會很長，預(yù)計比只考慮推進(jìn)劑建模、只計算微重力滑行期間的簡化模型長10倍。如果采用10個CPU核計算，需要計算1a左右。這個時間對于運(yùn)載火箭的任務(wù)適應(yīng)性是很不利的。

采用只對推進(jìn)劑建模，只計算微重力滑行期間的簡化模型，可以大大節(jié)省計算資源。由此帶來兩個新的難點(diǎn)，一個是如何獲得準(zhǔn)確的初始溫度分布，另一個是如何確定合適的熱流邊界條件?？梢杂玫途S度模型解決這兩個難點(diǎn)。

首先，上面級主發(fā)動機(jī)一次工作段，流動主要是氣枕區(qū)注入氫氣，液氫從下面流出，氣枕溫度分布主要由氫氣注入和推進(jìn)劑液位下降共同決定，可以簡化成軸對稱流動與傳熱問題。因此，可以建立二維軸對稱簡化模型，初始條件由任務(wù)加注情況給定，快速計算得到主動飛行段末了時刻氣枕的溫度分布。然后以這個溫度分布，作為微重力滑行段三維直接數(shù)值計算模型的初始條件，計算得到相應(yīng)的壓力下降過程。采用上述模型簡化方法，計算得到的某飛行工況微重力滑行段開始時刻的氣枕溫度如圖3所示。

圖3 軸對稱簡化模型計算得到的某火箭液氫貯箱，微重力滑行開始時刻的溫度分布(底部為液體推進(jìn)劑，上部分為氫氣)Fig.3 Temperature distribution in a launcher hydrogen propellant tank, at the beginning of microgravity flight, calcul-ated from an axial symmetry simplified model

貯箱和泡沫層的傳熱效應(yīng)可以通過對流傳熱邊界條件進(jìn)行描述。為了獲得較為準(zhǔn)確的傳熱邊界條件，與上面模型簡化計算的方法類似，可以用低維近似法計算得到推進(jìn)劑在貯箱壁面處的傳熱邊界條件，單獨(dú)對貯箱和絕熱泡沫層進(jìn)行建模。采用二維平板模型，分別對與載荷艙接觸的頂部貯箱，內(nèi)部與氣枕接觸、外部與地球或太空輻射換熱的中上部貯箱等各部分貯箱單獨(dú)計算，得到貯箱外表面泡沫平衡溫度。然后假設(shè)泡沫溫度沿貯箱厚度方向呈線性分布，計算得到熱流密度，將該熱流密度通過等效對流換熱邊界條件。

邊界溫度計算結(jié)果如圖4所示，對與載荷倉接觸的頂部貯箱，內(nèi)側(cè)采用氣枕等溫邊界條件，外側(cè)載荷艙300K輻射條件，計算得到泡沫外表面溫度變化過程，最終平衡溫度在270K附近。這樣就得到了簡化三維非定常計算模型所需要的熱邊界條件。

圖4 二維平板簡化模型計算得到的某火箭液氫貯箱在微重力滑行開始前的泡沫外表面溫度變化Fig.4 Temperature variation at the outer side of the tank foam, before microgravity flight, calculated from a two-dimension simplified plane model

4 結(jié)論與討論

結(jié)合運(yùn)載火箭上面級微重力滑行期間的復(fù)雜流動和傳熱問題，總結(jié)了3種數(shù)值計算模型簡化方法。這些模型簡化方法可以加速數(shù)值模擬計算，同時，是否保持了工程型號設(shè)計需要的計算精度，需要與具體工程數(shù)據(jù)進(jìn)行對比分析，進(jìn)一步確認(rèn)。

工程經(jīng)驗(yàn)法適用于初始條件影響因素復(fù)雜，不確定因素多，直接數(shù)值模擬困難的情況。這時，通過火箭飛行數(shù)據(jù)和理論分析，結(jié)合工程人員的經(jīng)驗(yàn)，給定一個合適的初始條件，不需要對長時間飛行進(jìn)行數(shù)值計算，大幅減少了計算時間。

極限參數(shù)法應(yīng)用于系統(tǒng)存在短時間的非線性變化，其變化過程(如主發(fā)動機(jī)關(guān)機(jī)等)復(fù)雜，變化的最后狀態(tài)不確定的情況。通過分析該時刻狀態(tài)的極限情況，設(shè)置相應(yīng)的極限狀態(tài)參數(shù)，進(jìn)行保守設(shè)計計算，解決工程型號設(shè)計需求。

低維近似法適用于推進(jìn)劑流動和傳熱，雖然復(fù)雜，但是具有對稱性特征，可以通過低維模型簡化，計算得到需要的初始條件和邊界條件，大大簡化核心計算模型。

必須指出，以上3種模型簡化方法具有很強(qiáng)的經(jīng)驗(yàn)性，運(yùn)用時須十分慎重。理論計算科研人員需要加強(qiáng)和工程設(shè)計人員的討論和溝通，盡可能對簡化模型計算的結(jié)果與地面試驗(yàn)、相關(guān)飛行遙測數(shù)據(jù)進(jìn)行比較驗(yàn)證，以保證簡化模型的正確性。