李 崢，向紅軍

(北京航空航天大學(xué)宇航學(xué)院，北京 100191)

0 引言

固體火箭發(fā)動(dòng)機(jī)燃燒室和噴管中的流場(chǎng)是不可分割的統(tǒng)一流場(chǎng)。燃燒室中流速較低，壓力變化幅度小，可按不可壓流處理，燃?xì)庠趪姽苤械募铀倥蛎涍^(guò)程則為可壓縮流動(dòng)［1］。目前，大多數(shù)固體火箭發(fā)動(dòng)機(jī)藥柱幾何形狀都呈現(xiàn)明顯的三維特征，使燃面邊界及相應(yīng)的內(nèi)流場(chǎng)亦具有三維性質(zhì)。與此同時(shí)，發(fā)動(dòng)機(jī)在工作過(guò)程中會(huì)出現(xiàn)過(guò)載情況。文獻(xiàn)［2－3］指出，在高過(guò)載情況下，燃燒室內(nèi)部Al2O3顆粒偏向過(guò)載方向，并在藥柱表面沉積，受到影響的藥柱燃燒部位燃速增加，藥柱燃面平行層燃燒規(guī)律被打破，藥柱不同部位的燃速受顆粒運(yùn)動(dòng)的影響差異很大，使得流場(chǎng)、燃面以及燃燒不再保持軸對(duì)稱性。因此，有必要對(duì)固體火箭發(fā)動(dòng)機(jī)燃燒室噴管進(jìn)行全三維數(shù)值模擬。

以往對(duì)于求解固體火箭發(fā)動(dòng)機(jī)內(nèi)部流動(dòng)問(wèn)題，大多采用基于壓力的方法，將連續(xù)方程改寫為壓力修正方程，解出壓力場(chǎng)與速度場(chǎng)，這就是著名的SIMPLE算法［4］。由于固體火箭發(fā)動(dòng)機(jī)內(nèi)部流動(dòng)包含不可壓和可壓縮問(wèn)題，壓力修正方程的迭代矩陣往往呈現(xiàn)病態(tài)，這增加了求解難度。早期工程上多采用直接求解法以及Krylov子空間迭代法求解壓力修正方程，受限于其計(jì)算時(shí)間復(fù)雜度［5］。此類方法大多應(yīng)用于二維或網(wǎng)格數(shù)較少的偽三維情況，難以求解全三維大規(guī)模實(shí)際問(wèn)題。代數(shù)多重網(wǎng)格方法(Algebraic Multigrid—AMG)具有最優(yōu)的時(shí)間復(fù)雜度，即迭代收斂速度與問(wèn)題規(guī)模無(wú)關(guān)，計(jì)算量線性依賴于未知量個(gè)數(shù)。目前，該算法已成功用于求解計(jì)算流體力學(xué)問(wèn)題［6－7］，但在固體火箭發(fā)動(dòng)機(jī)模擬方面的相關(guān)報(bào)道較少。

本文將AMG方法引入原有基于SIMPLE算法的求解器，用于求解壓力修正方程。在介紹控制方程和離散方法的基礎(chǔ)上，給出了AMG方法的主要內(nèi)容;然后，將Gauss消元法、不完全LU分解BiCGStab方法及AMG方法用于雙燃速星型藥柱固體發(fā)動(dòng)機(jī)三維兩相流場(chǎng)仿真，對(duì)比了3種方法在不同網(wǎng)格規(guī)模下的計(jì)算效率;最后，給出了AMG方法下發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)計(jì)算結(jié)果，分析了顆粒相與過(guò)載對(duì)發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)的影響。

1 數(shù)學(xué)模型

表1給出了SIMPLE算法求解步驟及單步循環(huán)各步所用時(shí)間。其中，求解壓力修正方程的時(shí)間占用了總時(shí)間的86%，為SIMPLE算法的核心步驟。

表1 SIMPLE算法單次循環(huán)時(shí)間分布(AMG，100×40×20)Table 1 Time distribution of a single SIMPLE cycle(AMG，100×40×20)

1.1 氣相控制方程及離散

穩(wěn)態(tài)氣相通用控制方程組

式中 φ表示氣相通用變量，在三維坐標(biāo)下分別為速度u、v、w及溫度T;φ為1時(shí)代表連續(xù)方程;Γ表示擴(kuò)散系數(shù);Sφ為源項(xiàng);U為速度矢量;ρ為密度，由完全氣體狀態(tài)方程 p=ρRT 確定［3］。

采用基于同位網(wǎng)格的有限體積法離散氣相通用控制方程(1)，最終得到線性代數(shù)方程組:

其中的系數(shù)依賴于離散化過(guò)程中所使用的具體格式。對(duì)于可壓縮流動(dòng)，壓力不僅出現(xiàn)在動(dòng)量方程中，而且通過(guò)狀態(tài)方程直接與密度以及能量方程相關(guān)聯(lián)。密度的變化也會(huì)影響壓力變化，所以必須建立壓力-速度-密度的耦合關(guān)系式，即可壓縮形式的壓力校正方程:

基于SIMPLE算法的可壓縮流動(dòng)求解步驟與不可壓縮幾乎完全相同。其中:

從中可看出，該方程的系數(shù)矩陣不再具有對(duì)角占優(yōu)特性(aP≠n=∑nban)，且為病態(tài)方程組，其求解出現(xiàn)

(i)了困難。

1.2 顆粒相求解方法

采用拉格朗日顆粒軌道模型模擬Al2O3顆粒與連續(xù)相之間的動(dòng)量、能量交換，不考慮顆粒的燃燒、蒸發(fā)、碰撞、聚合等，只考慮破碎效應(yīng)。采用Crowe C T［8］所提出的PSIC方法(單元中的顆粒源項(xiàng)法)，進(jìn)行氣粒兩相間的耦合計(jì)算。

(1)顆粒破碎模型

當(dāng)顆粒在氣相的牽引下加速時(shí)，會(huì)逐漸變形直至破碎，采用韋伯?dāng)?shù)(Weber)［9］來(lái)控制顆粒破碎是否發(fā)生。文中簡(jiǎn)單地認(rèn)為，當(dāng)顆粒保持為液態(tài)，且韋伯?dāng)?shù)超過(guò)臨界值時(shí)，顆粒破碎成直徑為原來(lái)一半的粒子。

式中 Dk為顆粒直徑;ρ為氣相密度;為氣相速度;為顆粒速度;σ為表面張力，Al2O3顆粒在2 300 K時(shí)的表面張力為0.69 N/m。

當(dāng)韋伯?dāng)?shù)超過(guò)4時(shí)，球形液滴開(kāi)始變形，在We=10～20范圍內(nèi)粒子發(fā)生破碎。

(2)顆粒尺寸分布

采用 Braithwaite［10］和 Salita 等［11］提出的顆粒分布函數(shù)，描述推進(jìn)劑中Al2O3顆粒的直徑分布，顆粒尺寸分布為對(duì)數(shù)正態(tài)分布。當(dāng)采用質(zhì)量平均直徑為5 μm、標(biāo)準(zhǔn)差為0.46的8組顆粒群計(jì)算時(shí)，平均直徑及所代表的尺寸范圍和質(zhì)量分?jǐn)?shù)見(jiàn)表2。

表2 顆粒群平均直徑Table 2 Average diameter of particle cloud

2 代數(shù)多重網(wǎng)格方法

多重網(wǎng)格方法的思想起源于20世紀(jì)60年代Brakhage H［12］、Fedorenko R［13］和 Bachvalov N［14］的工作，最早出現(xiàn)的是幾何多重網(wǎng)格(Geometric Multigrid-GMG)方法。GMG方法通過(guò)在不同疏密網(wǎng)格下求解原方程組，平順高、低頻誤差，以達(dá)到加速收斂的效果。

AMG方法是在GMG方法的思想上延伸建立起來(lái)的，它不需要利用求解問(wèn)題的物理和幾何背景，僅利用方程組稀疏矩陣中的代數(shù)信息，構(gòu)造多重網(wǎng)格組成部件。這一特點(diǎn)使得AMG方法既保持了GMG方法良好的算法可擴(kuò)展能力，又可作為具有即插即用性質(zhì)的黑盒子解法器，大大增強(qiáng)了其魯棒性［15］。

2.1 基本原理

對(duì)于式(3)，可寫成如下的線性系統(tǒng)［16］:

式中 A 為 n×n矩陣;u，f∈Rn。

插值算子P將粗網(wǎng)格向量空間轉(zhuǎn)換為細(xì)網(wǎng)格向量空間Rnc→Rn，為n×nc矩陣;限制算子R將細(xì)網(wǎng)格向量空間轉(zhuǎn)換為粗網(wǎng)格向量空間，為插值算子的轉(zhuǎn)置R=PT。

盡管實(shí)際的多重網(wǎng)格格式多采用幾層網(wǎng)格，為簡(jiǎn)便起見(jiàn)，這里僅討論二層網(wǎng)格(即細(xì)網(wǎng)格Ωh和粗網(wǎng)格Ω2h)的V循環(huán)，并先介紹最簡(jiǎn)單的粗網(wǎng)格修正格式:

(1)在細(xì)網(wǎng)格上求解

迭代μ1步后，計(jì)算所得近似值的殘余:

(2)在粗網(wǎng)格上求解誤差方程:

(3)進(jìn)行粗網(wǎng)格修正:

(4)以此為初值，在細(xì)網(wǎng)格上再迭代μ2步，然后回到(1)，開(kāi)始下一個(gè)V循環(huán)，直到達(dá)到一定的收斂標(biāo)準(zhǔn)為止。

經(jīng)典Ac構(gòu)建按Galerkin方法得到，即Ac=PTAP。該方法為能量范數(shù)形式，可更好地減少校正后的誤差。

2.2 粗網(wǎng)格的選取與插值算子［17］

在大多數(shù)AMG方法中，嵌套的網(wǎng)格層由細(xì)網(wǎng)格變量集合與粗網(wǎng)格變量集合組成，從細(xì)變量集合中篩選粗變量的過(guò)程為粗網(wǎng)格的選取。粗網(wǎng)格的選取對(duì)于AMG方法性能有很大影響，其他部件的構(gòu)造往往直接依賴于網(wǎng)格粗化的結(jié)果。本文以C-AMG方法為例，介紹網(wǎng)格粗化的一般方法。

假設(shè)aij為系數(shù)矩陣A中的元素，若

則稱i點(diǎn)與j點(diǎn)是強(qiáng)連接，粗網(wǎng)格點(diǎn)由強(qiáng)連接方向上選出，同時(shí)文獻(xiàn)［18］還給出了粗化過(guò)程應(yīng)滿足2個(gè)原則:(1)對(duì)于當(dāng)前網(wǎng)格上不屬于粗網(wǎng)格的節(jié)點(diǎn)，它的每個(gè)強(qiáng)連接點(diǎn)都應(yīng)在粗網(wǎng)格上，或與粗網(wǎng)格上的至少一點(diǎn)有強(qiáng)連接;(2)粗網(wǎng)格點(diǎn)集應(yīng)為細(xì)網(wǎng)格點(diǎn)集的極大獨(dú)立集，即細(xì)網(wǎng)格點(diǎn)集中的任意2個(gè)粗點(diǎn)之間不存在強(qiáng)連接關(guān)系，在此前提下使得粗網(wǎng)格點(diǎn)集網(wǎng)格結(jié)點(diǎn)盡量多。文獻(xiàn)［16］中詳細(xì)描述了選取粗網(wǎng)格結(jié)點(diǎn)的具體過(guò)程，本文不再贅述。

到目前為止，算法的核心問(wèn)題變?yōu)槎x插值算P:

假設(shè)j點(diǎn)為細(xì)網(wǎng)格上的節(jié)點(diǎn)，而i為粗網(wǎng)格上的結(jié)點(diǎn)。若j點(diǎn)與i點(diǎn)重合，則令j點(diǎn)的誤差校正量等于粗網(wǎng)格上i點(diǎn)的值，若j點(diǎn)與i點(diǎn)不重合，可由與j點(diǎn)相鄰的粗網(wǎng)格點(diǎn)插值得到。

3 計(jì)算結(jié)果

算例所模擬的發(fā)動(dòng)機(jī)由星形藥柱、圓柱藥柱、過(guò)渡段及長(zhǎng)尾噴管組成(圖1)，以兩種燃速推進(jìn)劑實(shí)現(xiàn)單室雙推力。其中，后部的第1級(jí)為星形藥柱;第2級(jí)包括前端的圓柱藥柱及過(guò)渡段藥柱。第1級(jí)裝高燃速推進(jìn)劑，第2級(jí)裝低燃速推進(jìn)劑。為滿足導(dǎo)彈總體需求，采用了長(zhǎng)尾噴管。圖2給出了圓柱藥柱和星形藥柱橫截面的網(wǎng)格劃分。

圖1 發(fā)動(dòng)機(jī)結(jié)構(gòu)Fig.1 Configuration of solid rocket motor

圖2 網(wǎng)格劃分Fig.2 Computational grid

推進(jìn)劑由 65%AP、8%HTPB、13%RDX、10%鋁粉及4%癸二酸二辛脂組成，燃燒室平衡壓強(qiáng)為6.0 MPa，由熱力計(jì)算得出燃燒室溫度為3 302 K，只計(jì)算發(fā)動(dòng)機(jī)初始時(shí)刻的流場(chǎng)情況。為考慮固體火箭發(fā)動(dòng)機(jī)三維特性，對(duì)其施加一徑向載荷30 g。

3.1 代數(shù)多重網(wǎng)格方法數(shù)值效率對(duì)比

首先，對(duì)AMG方法數(shù)值計(jì)算效率進(jìn)行測(cè)試，將上述算例在不同網(wǎng)格密度下展開(kāi)，分別對(duì)比了直接求解法(高斯消元法-GE)、不完全LU分解BiCGStab方法以及AMG方法的計(jì)算效率。

圖3給出了二維網(wǎng)格下3種方法單次循環(huán)的時(shí)間消耗，收斂精度為10－8。從中可看出，對(duì)于網(wǎng)格數(shù)較少的情況，AMG方法求解效率不明顯，3種方法耗時(shí)相當(dāng);隨著網(wǎng)格數(shù)增加，3種方法的耗時(shí)均出現(xiàn)了不同程度增加，AMG方法的時(shí)間增長(zhǎng)率遠(yuǎn)小于高斯消元法和BiCGStab方法;在200×80網(wǎng)格下，AMG方法取代直接求解法，成為耗時(shí)最短的算法。圖4分別給出了BiCGStab方法與AMG方法單次循環(huán)迭代收斂次數(shù)?？煽闯?，隨網(wǎng)格數(shù)增加，AMG方法收斂次數(shù)基本不變，而B(niǎo)iCGStab方法收斂次數(shù)大大增加，增加了該算法計(jì)算時(shí)間。

表3給出了三維網(wǎng)格下3種方法單次循環(huán)的計(jì)算效率，收斂精度為10－5。從中可看出，直接求解法并不適用于三維流場(chǎng)計(jì)算，而AMG方法不受計(jì)算規(guī)模的影響。在三維網(wǎng)格下，其收斂效率大大提高，收斂時(shí)間僅為BiCGStab方法的1/20，且仍能保持10次左右迭代收斂。

圖3 二維網(wǎng)格下單次循環(huán)計(jì)算時(shí)間對(duì)比Fig.3 Comparison of single cycle computational cost in 2D

圖4 二維網(wǎng)格下單次循環(huán)迭代次數(shù)對(duì)比Fig.4 Comparison of single cycle iteration in 2D

表3 三維情況下單次循環(huán)計(jì)算效率對(duì)比Table 3 Comparison of single cycle computational efficiency in 3D

3.2 雙燃速星型藥柱固體火箭發(fā)動(dòng)機(jī)內(nèi)部流場(chǎng)分析

該算例的網(wǎng)格劃分為200×40×20，采用AMG算法，可大大提高收斂效率，總計(jì)算時(shí)長(zhǎng)不到BiCGStab方法的1/20。

圖5給出了發(fā)動(dòng)機(jī)內(nèi)部流場(chǎng)縱截面的溫度分布，圖6給出了燃燒室后部及長(zhǎng)尾噴管內(nèi)馬赫數(shù)的分布。固體推進(jìn)劑藥柱燃燒生成的高溫高壓燃?xì)馓畛湔麄€(gè)燃燒室，隨后在長(zhǎng)尾噴管作用下，轉(zhuǎn)化為高速燃?xì)鈬姵觯瑴囟冉档?，馬赫數(shù)上升。在噴管擴(kuò)張段靠近壁面區(qū)域出現(xiàn)高速低溫區(qū)，這是顆粒相作用的結(jié)果。顆粒的物理特性決定了顆粒不能通過(guò)膨脹加速，將熱能轉(zhuǎn)換為動(dòng)能。受慣性的影響，噴管擴(kuò)張段對(duì)顆粒的作用降低，出現(xiàn)無(wú)顆粒區(qū)。在顆粒為主流的核心區(qū)域，顆粒對(duì)氣體傳熱，減緩氣體溫度降低的趨勢(shì)，并阻礙氣體流動(dòng)，同時(shí)將其擠壓到無(wú)顆粒區(qū);在以燃?xì)鉃橹鞯臒o(wú)顆粒區(qū)，氣體溫度較低，流速較高。

圖5 溫度分布Fig.5 Temperature distrbution

圖6 燃燒室后部及長(zhǎng)尾噴管內(nèi)馬赫數(shù)分布Fig.6 Mach contour at the rear parts of chamber and long-tail nozzle

圖7給出了4個(gè)不同位置上的橫截面云圖，分別位于二級(jí)圓柱藥柱(Z=0.9 m)、一級(jí)星型藥柱(Z=1.6 m)、長(zhǎng)尾噴管喉部(Z=1.88 m)及噴管出口位置(Z=2.0 m)。在星形藥柱及長(zhǎng)尾噴管內(nèi)，流動(dòng)均呈現(xiàn)明顯的三維特征，顆粒分布保持星角對(duì)稱性，這是由星邊燃面的周向加質(zhì)造成的，顆粒的分布使得流場(chǎng)內(nèi)溫度的分布亦成星角型;對(duì)于過(guò)載問(wèn)題，在圓柱藥柱燃燒室內(nèi)軸向速度與粒子分布明顯地偏向于過(guò)載作用方向(圖7(a))，而噴管內(nèi)的流場(chǎng)分布受過(guò)載影響不大。在燃燒室內(nèi)部，燃?xì)馀c顆粒相的流動(dòng)較慢，顆粒受過(guò)載作用較大，軸向速度較小，更容易發(fā)生偏轉(zhuǎn);在噴管內(nèi)部，燃?xì)饬鲃?dòng)較快，顆粒主要受氣相作用，軸向流動(dòng)速度較大，受過(guò)載的影響較少。

4 結(jié)論

(1)高斯消元法和BiCGStab方法在網(wǎng)格數(shù)較少的計(jì)算時(shí)具有簡(jiǎn)單高效的特性，但在三維計(jì)算中耗時(shí)較多;AMG方法計(jì)算收斂時(shí)間隨網(wǎng)格數(shù)增加線性增加，在數(shù)值模擬三維固體火箭發(fā)動(dòng)機(jī)內(nèi)部流動(dòng)問(wèn)題時(shí)，收斂時(shí)間遠(yuǎn)小于高斯消元法和BiCGStab方法;

(2)對(duì)于過(guò)載問(wèn)題，在圓柱藥柱燃燒室內(nèi)軸向速度與粒子分布明顯地偏向于過(guò)載作用方向，而噴管內(nèi)的流場(chǎng)分布受過(guò)載影響不大;

(3)顆粒的加入使得星形藥柱固體火箭發(fā)動(dòng)機(jī)噴管內(nèi)流動(dòng)呈星角對(duì)稱，顆粒對(duì)周圍燃?xì)庥袛D壓加熱作用，燃?xì)庠谠搮^(qū)域溫度升高，速度降低;在以燃?xì)鉃橹鞯膰姽鼙诿娓浇細(xì)鉁囟容^低，流速較高。

圖7 不同橫截面上的顆粒密度、溫度及燃?xì)廨S向速度分布Fig.7 Particle density distribution，temperature and gas axial velocity at different cross sections

［1］向紅軍，方國(guó)堯，崔濟(jì)亞.固體火箭發(fā)動(dòng)機(jī)燃燒室噴管統(tǒng)一流場(chǎng)計(jì)算［J］.固體火箭技術(shù)，2000，22(3).

［2］Langhenry M T，Marietta M.Acceleration effects in solid propellant rocket motors［R］.AIAA 86-1577.

［3］郭顏紅，梁曉庚，陳斌.雙燃速星孔藥柱長(zhǎng)尾噴管發(fā)動(dòng)機(jī)三維兩相流場(chǎng)數(shù)值模擬［J］.固體火箭技術(shù)，2007，30(3).

［4］蔡國(guó)彪，王慧玉.固體火箭發(fā)動(dòng)機(jī)燃燒室三維流場(chǎng)數(shù)值模擬［J］.推進(jìn)技術(shù)，1995(5).

［5］郭飛宵，楊力，劉榮，等.大型稀疏法方程組的代數(shù)多重網(wǎng)格解法［J］.測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2012，29(1).

［6］MacLachlan S P，Tang J M，Vuik C.Fast and robust solvers for pressure-correction in bubbly flow problems［J］.Journal of Computational Physics，2008.

［7］Darwish M，Sraj I，Moukalled F.A coupled finite volume solver for the solution of incompressible flows on unstructured grids［J］.Journal of Computational Physics，2009.

［8］周力行.湍流兩相流動(dòng)與燃燒的數(shù)值模擬［M］.北京:清華大學(xué)出版社，1991.

［9］Clift R，Grace J R.Bubbles，drops and particles［M］.New York:Academic Press，1978.

［10］Braithwaite P C.Quench bomb investigation of aluminum oxide formation from solid rocket propellants(part I):experimental methodology［C］//Proceeding of the 25th JANNAF Combustion Meeting，CPIA 498，1988.

［11］Salita.Quench bomb investigation of Al2O3formation from solid rocket propellants(part II):analysis of data［C］//25th JANNAF Combustion Meeting，CPIA 498，1988.

［12］Brakhage H.Uber die numerishe behandlung von integralgleichungen nach der quadraturformelmethode［J］.Numer.Math.，1960.

［13］Fedorenko R P.A relaxation method for solving elliptic difference equations［J］.USSR Comput.Math.Phys.，1961.

［14］Bachvalov N S.On the convergence of a relaxation method with natural constraints on the elliptic operator［J］.USSR Comput.Math.Phys.

［15］徐小文.可擴(kuò)展并行代數(shù)多重網(wǎng)格算法研究［D］.中國(guó)工程物理研究院，2007.

［16］Falgout R D.An introduction to algebraic multigrid［R］.UCRL-JRNL-220851.

［17］Ashby S F，F(xiàn)algout R D.A parallel multigrid preconditioned conjugate gradient algorithm for groundwater flow simulations［R］.UCRL-JC-122359，1995.

［18］Ruge J W，Stuben K.Algebraic multigrid［M］.In Multigrid Methods，F(xiàn)rontiers Appl.Math.3，McCormick S F editor，SIAM，PhiladelPhia，1987:73-130.