孫維佳,鐘 偉,王平陽(yáng),蔡 健,Micheal. F. MODEST

(1. 上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，上海 200240; 2. 上海宇航系統(tǒng)工程研究所，上海 201109；3. 懷俄明大學(xué) 機(jī)械學(xué)院，懷俄明州 82071; 4. 加利福尼亞大學(xué) 默塞德分校工學(xué)院，加利福尼亞州 95344)

0 引言

超燃沖壓發(fā)動(dòng)機(jī)被視為航空航天領(lǐng)域最有希望的革命性技術(shù)[1]，部分歸功于其不攜帶氧化劑的推力方式使得飛行器步入太空變得更加可行。典型scramjet的結(jié)構(gòu)中，進(jìn)氣道負(fù)責(zé)壓縮來(lái)流空氣，而尾噴管將高溫燃?xì)馀蛎浖铀俸笈懦觯瑥亩@得推力[2]。劉興洲[3]介紹了高超聲速流動(dòng)和超燃沖壓發(fā)動(dòng)機(jī)的基礎(chǔ)原理和進(jìn)展。由scramjet的工作原理可知，燃燒室是scramjet設(shè)計(jì)中最重要的部分[4]。在scramjet的發(fā)展中遇到諸多挑戰(zhàn)，比如復(fù)雜的超聲速氣流特性、穩(wěn)焰、毫秒級(jí)燃燒過(guò)程中的摻混及熱負(fù)荷的準(zhǔn)確計(jì)算等。楊事民等[5]就曾為提高穩(wěn)焰效果研究過(guò)一種臺(tái)階和凹腔組合的燃燒室內(nèi)的燃燒和流動(dòng)情況。為了克服這些挑戰(zhàn)，對(duì)燃燒室內(nèi)熱輻射的作用的研究不可小覷。

本文利用文獻(xiàn)[9]的燃燒室?guī)缀文Ｐ筒⒒贠penFOAM平臺(tái)，首先采用基于k-ω湍流模型的剪切應(yīng)力輸運(yùn)(SST)模型及隨機(jī)歐拉解(SEF)模型計(jì)算流場(chǎng)和燃燒，以解耦方式采用P1 模型計(jì)算輻射熱效應(yīng)。模擬結(jié)果與文獻(xiàn)[9]結(jié)果進(jìn)行對(duì)比以驗(yàn)證所用求解器的正確性。在此基礎(chǔ)上，采用將熱輻射與燃燒和流動(dòng)耦合計(jì)算的方式，考察了燃?xì)鉁囟取⒈诿孑椛錈崃鞯淖兓?，進(jìn)而分析與解耦計(jì)算的不同。通過(guò)本文及后續(xù)工作，希望能為采用基于OpenFOAM平臺(tái)的自定義求解器以耦合方式求解超聲速燃燒及熱輻射提供新的方式。

1 燃燒流動(dòng)及熱輻射模型

1.1 OpenFOAM及求解器簡(jiǎn)介

OpenFOAM(http://www.openfoam.com/)是基于C++語(yǔ)言開(kāi)發(fā)的面向?qū)ο蟮拈_(kāi)源CFD軟件。對(duì)大多數(shù)工程和自然科學(xué)領(lǐng)域，OpenFOAM都自帶大量的求解器；同時(shí)，用戶(hù)可根據(jù)自身需要對(duì)標(biāo)準(zhǔn)求解器進(jìn)行修改，如更改邊界條件、湍流模型等；其開(kāi)源特性更是允許用戶(hù)按其語(yǔ)法規(guī)則創(chuàng)建自定義求解器，這些自定義求解器可任意調(diào)用OpenFOAM內(nèi)置的庫(kù)及頭文件。這些特性賦予了OpenFOAM極強(qiáng)的擴(kuò)展性、實(shí)用性和可移植性。另外，其數(shù)據(jù)封裝方式也較好地保護(hù)了求解器中的私有類(lèi)被惡意修改和非法使用，從而保證了求解器的穩(wěn)定性和安全性。

本文使用的求解器是在OpenFOAM中基于壓力的有限體積反應(yīng)流求解器基礎(chǔ)上建立的，具有分析超聲速氣流燃燒和換熱情況的功能。該求解器可以指定計(jì)算域被劃分的個(gè)數(shù)，每一小塊計(jì)算域中都包括對(duì)氣體質(zhì)量分?jǐn)?shù)及焓的求解，從而實(shí)現(xiàn)并行計(jì)算，大大縮減求解時(shí)間。在改進(jìn)后的求解器能量方程中，采用顯焓的方式對(duì)原方程進(jìn)行改進(jìn)，即

(1)

1.2 燃燒模型

本文算例選用的燃料為H2，其當(dāng)量比為0.033；選用的燃燒模型由9種反應(yīng)物、19種反應(yīng)組成(見(jiàn)表1)，燃燒室所有壁面均設(shè)為非滑移壁面。由于和流動(dòng)的時(shí)間尺度相比，H2燃燒的時(shí)間非常短，因此對(duì)湍流-化學(xué)反應(yīng)模型進(jìn)行了簡(jiǎn)化。選用了SEF法作為湍流燃燒模型，這是一種耦合了RANS/LES解法的求解高速湍流流動(dòng)的歐拉輸運(yùn)偏微分解法。根據(jù)之前驗(yàn)證模擬的結(jié)果可知，與單純有限體積法計(jì)算所得的結(jié)果相比，SEF方法得到的溫升規(guī)律及主要物質(zhì)的濃度和實(shí)驗(yàn)結(jié)果更為吻合[19]。表1中A為指前因子，n為反應(yīng)級(jí)數(shù)，E為活化能。

另外，相關(guān)的施工人員通過(guò)借鑒其他建筑項(xiàng)目案例，引進(jìn)先進(jìn)的技術(shù)設(shè)備，避免對(duì)周?chē)沫h(huán)境造成影響。比如，我們可以利用永久性圍墻來(lái)替代臨時(shí)性的圍墻，這樣不僅減少了資源的消耗，減少了項(xiàng)目的成本，而且也降低了對(duì)周?chē)h(huán)境的負(fù)面影響。我們還可以增加相應(yīng)的灑水裝置，以免由于大量的灰塵影響周?chē)用裆畹沫h(huán)境。最為常用的就是懸掛式的噴灑裝置，施用起來(lái)比較方便，撒水量也會(huì)比較均勻，避免那種由于噴灑量過(guò)多形成水坑。當(dāng)然，施工環(huán)節(jié)的噪聲也是需要我們密切注意的問(wèn)題，設(shè)置噪聲預(yù)警裝置，如果超出規(guī)定范圍就會(huì)報(bào)警題型[4]。

表1 H2-空氣的燃燒動(dòng)力學(xué)模型[19]

1.3 輻射光譜模型

由于光譜吸收系數(shù)是每個(gè)波長(zhǎng)的吸收系數(shù)的總和，而準(zhǔn)確地將所有波長(zhǎng)的吸收系數(shù)求解和整合是項(xiàng)龐大而艱巨的任務(wù)，因此對(duì)光譜模型進(jìn)行求解具有重要意義。選取窄譜帶 k-distribution模型作為光譜模型[20-21]，該模型通過(guò)重組吸收系數(shù)將其轉(zhuǎn)換為光滑函數(shù), 進(jìn)而進(jìn)行求解。

(2)

式中:η表示波數(shù)。由于窄譜帶中任何僅取決于吸收系數(shù)的平均值均可由此公式表示，在任意兩個(gè)光譜間積分可得

(3)

比較以上兩個(gè)方程，當(dāng)κη=k時(shí)可得

(4)

(5)

式中:k(g)是g(k)的反函數(shù)。將式(4)代入g(k)的表達(dá)式,得

(6)

式中:H(k)函數(shù)的定義為

(7)

至此，g(k)表示k值以下波長(zhǎng)的光譜吸收系數(shù)；重排列后的吸收系數(shù)k(g)就是一個(gè)光滑、單調(diào)遞增的函數(shù)，其最大最小值與κη(η)保持一致。

當(dāng)k在一個(gè)較小的范圍內(nèi)變化時(shí)，

dg(kj)=f(kj)δkj=

(8)

式(8)即本文求解器采用的光譜模型公式。

2 幾何模型及計(jì)算條件

本文選取的幾何模型為三維矩形，主要的幾何結(jié)構(gòu)參數(shù)和來(lái)流參數(shù)均取自文獻(xiàn)[8]和文獻(xiàn)[9]，如圖1及表2所示。在該結(jié)構(gòu)中，H2由燃燒室側(cè)面的噴射器噴入。該模型的寬度為10 cm；總長(zhǎng)為140 cm，其中100 cm為燃燒室長(zhǎng)度。采用的矩形結(jié)構(gòu)網(wǎng)格總數(shù)為88 900，并且在燃燒室壁面處加密(見(jiàn)圖2)。

表2 模擬計(jì)算采用參數(shù)

圖1 所用模型幾何參數(shù)Fig.1 Geometric parameters of the model

圖2 模擬計(jì)算采用網(wǎng)格Fig.2 Grids used in simulated calculation

3 結(jié)果分析

3.1 非耦合計(jì)算結(jié)果分析

采用上述與文獻(xiàn)[9]一致的幾何參數(shù)及來(lái)流條件進(jìn)行燃燒/流動(dòng)與熱輻射的解耦計(jì)算，分別從Ma、溫度、壓力及主要物質(zhì)濃度的分布等與文獻(xiàn)[9]的結(jié)果比較，以驗(yàn)證本文改進(jìn)后的求解器和模擬結(jié)果的正確性。

圖3給出了本文算例中Ma的分布情況，其變化規(guī)律和文獻(xiàn)[9]保持一致。由圖3和圖4所示，激波在x≈58 cm處產(chǎn)生；由于氣體混合燃燒反應(yīng)的進(jìn)行，激波在x≈64 cm處首次從燃燒室壁面反射，緊接著的第二次反射發(fā)生在x≈128 cm處。這兩個(gè)位置均與文獻(xiàn)[9]基本一致。但由于文獻(xiàn)[9]并未給出所采用的燃燒模型及流體物性設(shè)置，因此，本文采用的燃燒模型及流體黏度設(shè)置等可能與之不同，致使結(jié)果數(shù)值上出現(xiàn)一些不同，如激波的頂點(diǎn)比文獻(xiàn)[9]中的要靠后等。

圖4分別為在不考慮輻射和考慮輻射兩種情況下的溫度分布。其中，考慮輻射的溫度分布圖是在后處理時(shí)通過(guò)計(jì)算將溫度的影響加入不考慮輻射情況下的溫度分布圖得到的。隨燃燒反應(yīng)的開(kāi)始，燃燒室內(nèi)溫度迅速升高。由圖4可明顯看出，由于輻射的作用，原來(lái)集中于壁面的熱量向燃燒室內(nèi)部轉(zhuǎn)移，由于燃?xì)獍l(fā)射熱輻射，圖4(b)燃?xì)庾罡邷囟冉档土?20 K，燃?xì)庠诎l(fā)射熱輻射的同時(shí)，也吸收熱輻射，使燃燒室內(nèi)高溫區(qū)域的范圍遠(yuǎn)大于圖4(a)的，表明考慮輻射的燃燒和流場(chǎng)參數(shù)更符合實(shí)際情況。

圖4 溫度分布云圖Fig.4 Cloud image of the temperature distribution

圖5 壓力分布云圖Fig.5 Cloud image of the pressure distribution

燃燒室內(nèi)的壓力分布如圖5所示。由圖5可知，激波從壁面反射時(shí)損失相當(dāng)一部分能量，且圖5中激波反射的位置與圖3吻合，并在出口處有明顯的壓力增加。由圖4和圖5可知，空氣與H2較好地混合點(diǎn)位于x≈96 cm處，與文獻(xiàn)[9]中的位置(x≈90 cm)吻合較好。

圖6～8給出了H2，H2O和OH的摩爾分?jǐn)?shù)分布，這三種物質(zhì)的摩爾分?jǐn)?shù)及分布反映了燃燒的程度及點(diǎn)火位置。由圖6～8可知，點(diǎn)火發(fā)生在x≈64 cm處，即激波在燃燒室壁面反射附近。H2O的摩爾分?jǐn)?shù)達(dá)23%，OH和H2O的分布情況表明燃燒在噴射器附近開(kāi)始后，隨燃?xì)獠粩嘞蛉紵液蠖瓮七M(jìn)，且燃燒的程度比較令人滿(mǎn)意。而H2的分布主要集中在噴射器附近，從另一方面表明燃燒程度較高。

圖6 OH摩爾分?jǐn)?shù)分布云圖Fig.6 Cloud image of the mole fraction distribution of OH

圖7 H2O摩爾分?jǐn)?shù)分布云圖Fig.7 Cloud image of the mole fraction distribution of H2O

圖8 H2摩爾分?jǐn)?shù)分布云圖Fig.8 Cloud image of the mole fraction distribution of H2

圖9 壁面輻射熱流密度Fig.9 Wall radiative heat flux density

壁面輻射熱流密度(Qr)和對(duì)流熱流密度(Qc)如圖9和圖10所示。對(duì)于二維模型，熱輻射可以從燃?xì)獾母邷夭糠謧鬟f到其x方向兩側(cè)的壁面，因此，相對(duì)文獻(xiàn)[9]來(lái)講，兩段偏低，中間偏高；對(duì)流熱流密度則趨勢(shì)一致，但由于前面所述模型等方面的原因，具體數(shù)值上有些差別。本文輻射熱流密度的趨勢(shì)與文獻(xiàn)[12]和文獻(xiàn)[21]的計(jì)算結(jié)果也是一致的，結(jié)合之前對(duì)溫度、壓力等云圖的分析比較，本文所用求解器和所得到的結(jié)果具有一定的可信度。由圖9可知，壁面輻射熱流密度的平均值為45.6 W/cm2，最大值約為55 W/cm2；這個(gè)結(jié)果雖然比文獻(xiàn)[9]中的70 W/cm2要小，卻遠(yuǎn)大于文獻(xiàn)[8]中的24 W/cm2。文獻(xiàn)[8]中結(jié)果較小可能由于其將燃燒室內(nèi)氣體假設(shè)為常物性。計(jì)算輻射熱流密度時(shí)所設(shè)定的出口黑度、壁面吸收率及反射率、氣體黏度等模型所需參數(shù)的不同，也是導(dǎo)致當(dāng)前結(jié)果與文獻(xiàn)[9]結(jié)果差異的原因。但是這些算例及曲線(xiàn)圖均表明，對(duì)燃燒室內(nèi)熱輻射在壁面某些部位的作用應(yīng)當(dāng)給予充分的重視。

圖10 壁面對(duì)流熱流密度Fig.10 Wall convective heat flux density

圖11給出了燃燒室內(nèi)壁面對(duì)流熱流密度和輻射熱流密度的對(duì)比。由圖可知，在壁面的不同位置處，輻射熱流密度占對(duì)流熱流密度的比值有很大區(qū)別。在x≈120 cm處，考慮計(jì)算偏差后，所得輻射熱流密度占對(duì)流熱流密度值的30%以上。圖11進(jìn)一步表明熱輻射對(duì)燃燒室壁面熱防護(hù)的重要性。

圖11 壁面輻射熱流密度和對(duì)流熱流密度對(duì)比圖Fig.11 Comparison between wall radiative heat flux density and wall convective heat flux density

3.2 耦合計(jì)算結(jié)果分析

解耦算法具有計(jì)算效率高、操作方便等優(yōu)點(diǎn)，但燃燒/流動(dòng)與熱輻射本質(zhì)上是耦合的過(guò)程，而由于超聲速流場(chǎng)中大量激波和渦系等的存在，現(xiàn)有商業(yè)軟件在超聲速燃燒耦合方面計(jì)算能力有限，目前的國(guó)內(nèi)外文獻(xiàn)中也鮮有耦合計(jì)算的先例。為研究解耦計(jì)算與實(shí)際的不同情況，本文也對(duì)耦合計(jì)算開(kāi)展了初步研究。圖12給出了耦合計(jì)算的溫度分布，與圖4(b)對(duì)比可知，燃燒室后半段內(nèi)溫度分布更為均勻，燃?xì)獾淖罡邷囟扰c解耦計(jì)算時(shí)相比進(jìn)一步降低了78 K，與不加輻射時(shí)相比更是下降了近200 K。壁面高溫區(qū)面積明顯縮小。比較圖12與圖4(b)中x=40～60 cm處點(diǎn)火區(qū)域的圖像可知，耦合計(jì)算在燃燒室中心點(diǎn)火區(qū)域附近的溫度有明顯升高，該結(jié)論表明在scramjet的點(diǎn)火及穩(wěn)焰研究中熱輻射的作用值得深入探究。圖13給出了耦合計(jì)算與解耦計(jì)算中壁面輻射熱流密度的對(duì)比，該圖表明耦合計(jì)算壁面輻射熱流的分布趨向更平穩(wěn)，所得燃燒室內(nèi)溫度分布更均勻。

圖12 耦合計(jì)算時(shí)溫度分布云圖Fig.12 Cloud image of the temperature distribution in couple calculation

圖13 耦合及解耦計(jì)算壁面輻射熱流密度對(duì)比圖Fig.13 Comparison of the radiative heat flux density in couple and decouple calculation

4 結(jié)束語(yǔ)

為計(jì)算熱輻射對(duì)scramjet燃燒室壁面熱流的作用，基于OpenFOAM平臺(tái)，采用了三維矩形H2燃料scramjet模型，結(jié)合基于k-ω湍流模型的SST模型和SEF湍流燃燒模型進(jìn)行流場(chǎng)計(jì)算。燃?xì)夤庾V特性取自k-distribution模型，二維P1模型用于計(jì)算輻射輸運(yùn)過(guò)程。在燃燒/流動(dòng)和熱輻射解耦計(jì)算模型基礎(chǔ)上，嘗試開(kāi)展耦合計(jì)算，得到的主要結(jié)論如下：

1) 通過(guò)與文獻(xiàn)[9]結(jié)果比較可知，本文所得壁面輻射熱流密度曲線(xiàn)的變化趨勢(shì)和該文獻(xiàn)結(jié)果基本一致，圖3～8給出的燃燒室內(nèi)幾種主要參數(shù)的云圖與文獻(xiàn)[9]的基本一致，本文模型和改進(jìn)求解器的正確性得到驗(yàn)證。

2) 進(jìn)行燃燒/流動(dòng)和熱輻射的解耦計(jì)算可知，考慮輻射后，原來(lái)集中于壁面處的熱量向燃燒室內(nèi)部擴(kuò)散，且燃?xì)庾罡邷囟认陆颠_(dá)120 K (圖4)；H2、OH和H2O的摩爾分?jǐn)?shù)分布(見(jiàn)圖6～8)表明，本文求解器得到的燃燒情況良好。壁面輻射熱流密度最高可達(dá)55 W/cm2，且可占對(duì)流熱流密度的40%以上，因此燃燒室內(nèi)熱輻射的影響不應(yīng)該被忽略。

3) 由燃燒/流動(dòng)與熱輻射的耦合計(jì)算結(jié)果可知，燃燒室中段及后段溫度分布更加均勻，燃?xì)庾罡邷囟缺炔豢紤]輻射時(shí)相比進(jìn)一步降低，降幅達(dá)到200 K左右；并且耦合計(jì)算中點(diǎn)火區(qū)域溫度升高。將耦合計(jì)算的壁面輻射熱流密度曲線(xiàn)與解耦計(jì)算所得曲線(xiàn)對(duì)比可知(見(jiàn)圖13)，耦合計(jì)算中壁面輻射熱流的分布趨于平穩(wěn)。由此可知，在燃燒室中對(duì)燃燒/流動(dòng)和熱輻射的耦合計(jì)算有必要進(jìn)行深入探討。

[1] FRY R S. A century of ramjet propulsion technology evolution[J]. Journal of Propulsion and Power, 2004, 20(1): 27-58.

[2] BERGLUND M, FEDINA E, FUREBY C, et al. Finite rate chemistry large-eddy simulation of self-ignition in supersonic combustion ramjet[J]. AIAA Journal, 2010, 48(3): 540-550.

[3] 劉興洲. 超燃沖壓發(fā)動(dòng)機(jī)性能初步研究[J]. 航空發(fā)動(dòng)機(jī), 2007, 33(2): 1-4.

[4] CHAN J, SISLIAN J P, ALEXANDER D. Numerically simulated comparative performance of a scramjet and shcramjet at Mach 11[J]. Journal of Propulsion and Power, 2010, 26(5): 1125-1134.

[5] 楊事民, 張建良. 超燃沖壓發(fā)動(dòng)機(jī)燃燒室流場(chǎng)數(shù)值模擬[J]. 航空發(fā)動(dòng)機(jī), 2009, 35(4): 25-28.

[6] TIWARI S N, PIDUGU S B, MOHIELDIN T O. Radiative interactions in supersonic flows of premixed hydrogen in expanding nozzles[C]∥37th Aerospace Sciences Meeting and Exhibit. Reno:[s.n.], 1999: 1052.

[7] MANI M, TIWARI S N, DRUMMOND J P. Investigation of chemically reacting and radiating supersonic flow in channels[J]. Applied Scientific Research, 1993, 50(1): 43-68.

[8] NELSON H F. Radiative heating in scramjet combustors[J]. Journal of Thermophysics and Heat Transfer, 1997, 11(1): 59-64.

[9] KOTOV D V, SURZHIKOV S T. Computation of hypersonic flow and radiation of viscous chemically reacting gas in a channel modeling a section of a scramjet[J]. High Temperature, 2012, 50(1): 120-130.

[10] SURZHIKOV S, SHANG J. Radiative heat exchange in a hydrogen-fueled scramjet combustion chambers[C]∥51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. Grapevine:[s.n.], 2013: 448.

[11] SURZHIKOV S, SHANG J. Numerical prediction of convective and radiative heating of scramjet combustion chamber with hydrocarbon fuels[J]. 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 2013.

[12] CROW A J, BOYD I D, TERRAPON V E. Radiation modeling of a hydrogen fueled scramjet[J]. Journal of Thermophysics and Heat Transfer, 2013, 27(1): 11-21.

[13] 黃志偉, 何國(guó)強(qiáng), 秦飛, 等. 支板增強(qiáng)混合超聲速燃燒的大渦模擬研究[J]. 固體火箭技術(shù), 2015, 38(5): 664-670.

[14] 鞏偉杰, 唐碩, 李世珍. 高超聲速飛行器氣動(dòng)加熱三維數(shù)值分析方法研究[J]. 飛行力學(xué), 2011, 29(2): 78-81.

[15] 鮑文, 張聰, 秦江, 等. 超燃沖壓發(fā)動(dòng)機(jī)主被動(dòng)復(fù)合熱防護(hù)系統(tǒng)方案設(shè)計(jì)思考[J]. 推進(jìn)技術(shù), 2013, 34(12): 1659-1663.

[16] WANG J Y, GAO Z X, LEE C H, et al. A decoupled procedure for convection-radiation simulation in scramjets[J]. Science China: Technological Sciences, 2014, 57(12): 2551-2566.

[17] 魏鑫, 孫冰, 鄭力銘. 超燃沖壓發(fā)動(dòng)機(jī)前體/進(jìn)氣道熱結(jié)構(gòu)分析[C]∥中國(guó)宇航學(xué)會(huì)2005年固體火箭推進(jìn)年會(huì)第22屆年會(huì)論文集(發(fā)動(dòng)機(jī)分冊(cè)). [S.l.: s.n.], 2005.

[18] 艾青, 夏新林, 孫鳳賢. 壁面熱特性對(duì)超聲速燃燒室熱環(huán)境的影響[J]. 工程熱物理學(xué)報(bào), 2009, 30(8): 1373-1375.

[19] NASA. Eulerian transported PDF framework for scramjet flowpath analysis: Phase I project[OL]. (2012-02-01). https:∥techport.nasa.gov/view/9451

[20] LIU J, BROWN M. Radiative heating in hydrocarbon-fueled scramjet engines[C]∥48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit. Atlanta: [s.n.], 2012: 3775.

[21] CAI J, MARQUEZ R, MODEST M F. Comparisons of radiative heat transfer calculations in a jet diffusion flame using spherical harmonics and k-distributions[J]. Journal of Heat Transfer, 2014, 136(11): 112702.