超燃沖壓發(fā)動機燃燒室內對流與輻射加熱基于OpenFOAM的數(shù)值模擬

孫維佳,鐘 偉,王平陽,蔡 健,Micheal. F. MODEST

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

0 引言

超燃沖壓發(fā)動機被視為航空航天領域最有希望的革命性技術[1]，部分歸功于其不攜帶氧化劑的推力方式使得飛行器步入太空變得更加可行。典型scramjet的結構中，進氣道負責壓縮來流空氣，而尾噴管將高溫燃氣膨脹加速后排出，從而獲得推力[2]。劉興洲[3]介紹了高超聲速流動和超燃沖壓發(fā)動機的基礎原理和進展。由scramjet的工作原理可知，燃燒室是scramjet設計中最重要的部分[4]。在scramjet的發(fā)展中遇到諸多挑戰(zhàn)，比如復雜的超聲速氣流特性、穩(wěn)焰、毫秒級燃燒過程中的摻混及熱負荷的準確計算等。楊事民等[5]就曾為提高穩(wěn)焰效果研究過一種臺階和凹腔組合的燃燒室內的燃燒和流動情況。為了克服這些挑戰(zhàn)，對燃燒室內熱輻射的作用的研究不可小覷。

本文利用文獻[9]的燃燒室?guī)缀文Ｐ筒⒒贠penFOAM平臺，首先采用基于k-ω湍流模型的剪切應力輸運(SST)模型及隨機歐拉解(SEF)模型計算流場和燃燒，以解耦方式采用P1 模型計算輻射熱效應。模擬結果與文獻[9]結果進行對比以驗證所用求解器的正確性。在此基礎上，采用將熱輻射與燃燒和流動耦合計算的方式，考察了燃氣溫度、壁面輻射熱流的變化，進而分析與解耦計算的不同。通過本文及后續(xù)工作，希望能為采用基于OpenFOAM平臺的自定義求解器以耦合方式求解超聲速燃燒及熱輻射提供新的方式。

1 燃燒流動及熱輻射模型

1.1 OpenFOAM及求解器簡介

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

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

(1)

1.2 燃燒模型

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

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

表1 H2-空氣的燃燒動力學模型[19]

1.3 輻射光譜模型

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

(2)

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

(3)

比較以上兩個方程，當κη=k時可得

(4)

(5)

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

(6)

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

(7)

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

當k在一個較小的范圍內變化時，

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

(8)

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

2 幾何模型及計算條件

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

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

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

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

3 結果分析

3.1 非耦合計算結果分析

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

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

圖4分別為在不考慮輻射和考慮輻射兩種情況下的溫度分布。其中，考慮輻射的溫度分布圖是在后處理時通過計算將溫度的影響加入不考慮輻射情況下的溫度分布圖得到的。隨燃燒反應的開始，燃燒室內溫度迅速升高。由圖4可明顯看出，由于輻射的作用，原來集中于壁面的熱量向燃燒室內部轉移，由于燃氣發(fā)射熱輻射，圖4(b)燃氣最高溫度降低了120 K，燃氣在發(fā)射熱輻射的同時，也吸收熱輻射，使燃燒室內高溫區(qū)域的范圍遠大于圖4(a)的，表明考慮輻射的燃燒和流場參數(shù)更符合實際情況。

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

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

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

圖6～8給出了H2，H2O和OH的摩爾分數(shù)分布，這三種物質的摩爾分數(shù)及分布反映了燃燒的程度及點火位置。由圖6～8可知，點火發(fā)生在x≈64 cm處，即激波在燃燒室壁面反射附近。H2O的摩爾分數(shù)達23%，OH和H2O的分布情況表明燃燒在噴射器附近開始后，隨燃氣不斷向燃燒室后段推進，且燃燒的程度比較令人滿意。而H2的分布主要集中在噴射器附近，從另一方面表明燃燒程度較高。

圖6 OH摩爾分數(shù)分布云圖Fig.6 Cloud image of the mole fraction distribution of OH

圖7 H2O摩爾分數(shù)分布云圖Fig.7 Cloud image of the mole fraction distribution of H2O

圖8 H2摩爾分數(shù)分布云圖Fig.8 Cloud image of the mole fraction distribution of H2

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

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

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

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

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

3.2 耦合計算結果分析

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

圖12 耦合計算時溫度分布云圖Fig.12 Cloud image of the temperature distribution in couple calculation

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

4 結束語

為計算熱輻射對scramjet燃燒室壁面熱流的作用，基于OpenFOAM平臺，采用了三維矩形H2燃料scramjet模型，結合基于k-ω湍流模型的SST模型和SEF湍流燃燒模型進行流場計算。燃氣光譜特性取自k-distribution模型，二維P1模型用于計算輻射輸運過程。在燃燒/流動和熱輻射解耦計算模型基礎上，嘗試開展耦合計算，得到的主要結論如下：

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

2) 進行燃燒/流動和熱輻射的解耦計算可知，考慮輻射后，原來集中于壁面處的熱量向燃燒室內部擴散，且燃氣最高溫度下降達120 K (圖4)；H2、OH和H2O的摩爾分數(shù)分布(見圖6～8)表明，本文求解器得到的燃燒情況良好。壁面輻射熱流密度最高可達55 W/cm2，且可占對流熱流密度的40%以上，因此燃燒室內熱輻射的影響不應該被忽略。

3) 由燃燒/流動與熱輻射的耦合計算結果可知，燃燒室中段及后段溫度分布更加均勻，燃氣最高溫度比不考慮輻射時相比進一步降低，降幅達到200 K左右；并且耦合計算中點火區(qū)域溫度升高。將耦合計算的壁面輻射熱流密度曲線與解耦計算所得曲線對比可知(見圖13)，耦合計算中壁面輻射熱流的分布趨于平穩(wě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ā)動機性能初步研究[J]. 航空發(fā)動機, 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ā)動機燃燒室流場數(shù)值模擬[J]. 航空發(fā)動機, 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] 黃志偉, 何國強, 秦飛, 等. 支板增強混合超聲速燃燒的大渦模擬研究[J]. 固體火箭技術, 2015, 38(5): 664-670.

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

[15] 鮑文, 張聰, 秦江, 等. 超燃沖壓發(fā)動機主被動復合熱防護系統(tǒng)方案設計思考[J]. 推進技術, 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ā)動機前體/進氣道熱結構分析[C]∥中國宇航學會2005年固體火箭推進年會第22屆年會論文集(發(fā)動機分冊). [S.l.: s.n.], 2005.

[18] 艾青, 夏新林, 孫鳳賢. 壁面熱特性對超聲速燃燒室熱環(huán)境的影響[J]. 工程熱物理學報, 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.