，,，，

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

0 引言

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

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

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

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

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

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

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

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

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

(a) t1 (b) t2

(c) t3 (d) t4圖1 液氫貯箱主發(fā)動(dòng)機(jī)關(guān)機(jī)前液面附近的晃動(dòng)干擾簡(jiǎn)化，時(shí)間t1、t2、t3、t4分別是一個(gè)晃動(dòng)周期內(nèi)的4個(gè)不同時(shí)刻(深色為液氫，淺色為增壓氣體，矢量表示中心截面上流體的速度)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

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

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

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

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

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

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

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

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

圖3 軸對(duì)稱簡(jiǎn)化模型計(jì)算得到的某火箭液氫貯箱，微重力滑行開(kāi)始時(shí)刻的溫度分布(底部為液體推進(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)可以通過(guò)對(duì)流傳熱邊界條件進(jìn)行描述。為了獲得較為準(zhǔn)確的傳熱邊界條件，與上面模型簡(jiǎn)化計(jì)算的方法類似，可以用低維近似法計(jì)算得到推進(jìn)劑在貯箱壁面處的傳熱邊界條件，單獨(dú)對(duì)貯箱和絕熱泡沫層進(jìn)行建模。采用二維平板模型，分別對(duì)與載荷艙接觸的頂部貯箱，內(nèi)部與氣枕接觸、外部與地球或太空輻射換熱的中上部貯箱等各部分貯箱單獨(dú)計(jì)算，得到貯箱外表面泡沫平衡溫度。然后假設(shè)泡沫溫度沿貯箱厚度方向呈線性分布，計(jì)算得到熱流密度，將該熱流密度通過(guò)等效對(duì)流換熱邊界條件。

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

圖4 二維平板簡(jiǎn)化模型計(jì)算得到的某火箭液氫貯箱在微重力滑行開(kāi)始前的泡沫外表面溫度變化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)載火箭上面級(jí)微重力滑行期間的復(fù)雜流動(dòng)和傳熱問(wèn)題，總結(jié)了3種數(shù)值計(jì)算模型簡(jiǎn)化方法。這些模型簡(jiǎn)化方法可以加速數(shù)值模擬計(jì)算，同時(shí)，是否保持了工程型號(hào)設(shè)計(jì)需要的計(jì)算精度，需要與具體工程數(shù)據(jù)進(jìn)行對(duì)比分析，進(jìn)一步確認(rèn)。

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

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

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

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