超燃冲压发动机燃烧室内对流与辐射加热基于OpenFOAM的数值模拟

孙维佳,钟 伟,王平阳,蔡 健,Micheal. F. MODEST

(1. 上海交通大学 机械与动力工程学院，上海 200240; 2. 上海宇航系统工程研究所，上海 201109；3. 怀俄明大学 机械学院，怀俄明州 82071; 4. 加利福尼亚大学 默塞德分校工学院，加利福尼亚州 95344)

0 引言

超燃冲压发动机被视为航空航天领域最有希望的革命性技术[1]，部分归功于其不携带氧化剂的推力方式使得飞行器步入太空变得更加可行。典型scramjet的结构中，进气道负责压缩来流空气，而尾喷管将高温燃气膨胀加速后排出，从而获得推力[2]。刘兴洲[3]介绍了高超声速流动和超燃冲压发动机的基础原理和进展。由scramjet的工作原理可知，燃烧室是scramjet设计中最重要的部分[4]。在scramjet的发展中遇到诸多挑战，比如复杂的超声速气流特性、稳焰、毫秒级燃烧过程中的掺混及热负荷的准确计算等。杨事民等[5]就曾为提高稳焰效果研究过一种台阶和凹腔组合的燃烧室内的燃烧和流动情况。为了克服这些挑战，对燃烧室内热辐射的作用的研究不可小觑。

本文利用文献[9]的燃烧室几何模型并基于OpenFOAM平台，首先采用基于k-ω湍流模型的剪切应力输运(SST)模型及随机欧拉解(SEF)模型计算流场和燃烧，以解耦方式采用P1 模型计算辐射热效应。模拟结果与文献[9]结果进行对比以验证所用求解器的正确性。在此基础上，采用将热辐射与燃烧和流动耦合计算的方式，考察了燃气温度、壁面辐射热流的变化，进而分析与解耦计算的不同。通过本文及后续工作，希望能为采用基于OpenFOAM平台的自定义求解器以耦合方式求解超声速燃烧及热辐射提供新的方式。

1 燃烧流动及热辐射模型

1.1 OpenFOAM及求解器简介

OpenFOAM(http://www.openfoam.com/)是基于C++语言开发的面向对象的开源CFD软件。对大多数工程和自然科学领域，OpenFOAM都自带大量的求解器；同时，用户可根据自身需要对标准求解器进行修改，如更改边界条件、湍流模型等；其开源特性更是允许用户按其语法规则创建自定义求解器，这些自定义求解器可任意调用OpenFOAM内置的库及头文件。这些特性赋予了OpenFOAM极强的扩展性、实用性和可移植性。另外，其数据封装方式也较好地保护了求解器中的私有类被恶意修改和非法使用，从而保证了求解器的稳定性和安全性。

本文使用的求解器是在OpenFOAM中基于压力的有限体积反应流求解器基础上建立的，具有分析超声速气流燃烧和换热情况的功能。该求解器可以指定计算域被划分的个数，每一小块计算域中都包括对气体质量分数及焓的求解，从而实现并行计算，大大缩减求解时间。在改进后的求解器能量方程中，采用显焓的方式对原方程进行改进，即

(1)

1.2 燃烧模型

本文算例选用的燃料为H2，其当量比为0.033；选用的燃烧模型由9种反应物、19种反应组成(见表1)，燃烧室所有壁面均设为非滑移壁面。由于和流动的时间尺度相比，H2燃烧的时间非常短，因此对湍流-化学反应模型进行了简化。选用了SEF法作为湍流燃烧模型，这是一种耦合了RANS/LES解法的求解高速湍流流动的欧拉输运偏微分解法。根据之前验证模拟的结果可知，与单纯有限体积法计算所得的结果相比，SEF方法得到的温升规律及主要物质的浓度和实验结果更为吻合[19]。表1中A为指前因子，n为反应级数，E为活化能。

另外，相关的施工人员通过借鉴其他建筑项目案例，引进先进的技术设备，避免对周围的环境造成影响。比如，我们可以利用永久性围墙来替代临时性的围墙，这样不仅减少了资源的消耗，减少了项目的成本，而且也降低了对周围环境的负面影响。我们还可以增加相应的洒水装置，以免由于大量的灰尘影响周围居民生活的环境。最为常用的就是悬挂式的喷洒装置，施用起来比较方便，撒水量也会比较均匀，避免那种由于喷洒量过多形成水坑。当然，施工环节的噪声也是需要我们密切注意的问题，设置噪声预警装置，如果超出规定范围就会报警题型[4]。

表1 H2-空气的燃烧动力学模型[19]

1.3 辐射光谱模型

由于光谱吸收系数是每个波长的吸收系数的总和，而准确地将所有波长的吸收系数求解和整合是项庞大而艰巨的任务，因此对光谱模型进行求解具有重要意义。选取窄谱带 k-distribution模型作为光谱模型[20-21]，该模型通过重组吸收系数将其转换为光滑函数, 进而进行求解。

(2)

式中:η表示波数。由于窄谱带中任何仅取决于吸收系数的平均值均可由此公式表示，在任意两个光谱间积分可得

(3)

比较以上两个方程，当κη=k时可得

(4)

(5)

式中:k(g)是g(k)的反函数。将式(4)代入g(k)的表达式,得

(6)

式中:H(k)函数的定义为

(7)

至此，g(k)表示k值以下波长的光谱吸收系数；重排列后的吸收系数k(g)就是一个光滑、单调递增的函数，其最大最小值与κη(η)保持一致。

当k在一个较小的范围内变化时，

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

(8)

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

2 几何模型及计算条件

本文选取的几何模型为三维矩形，主要的几何结构参数和来流参数均取自文献[8]和文献[9]，如图1及表2所示。在该结构中，H2由燃烧室侧面的喷射器喷入。该模型的宽度为10 cm；总长为140 cm，其中100 cm为燃烧室长度。采用的矩形结构网格总数为88 900，并且在燃烧室壁面处加密(见图2)。

表2 模拟计算采用参数

图1 所用模型几何参数Fig.1 Geometric parameters of the model

图2 模拟计算采用网格Fig.2 Grids used in simulated calculation

3 结果分析

3.1 非耦合计算结果分析

采用上述与文献[9]一致的几何参数及来流条件进行燃烧/流动与热辐射的解耦计算，分别从Ma、温度、压力及主要物质浓度的分布等与文献[9]的结果比较，以验证本文改进后的求解器和模拟结果的正确性。

图3给出了本文算例中Ma的分布情况，其变化规律和文献[9]保持一致。由图3和图4所示，激波在x≈58 cm处产生；由于气体混合燃烧反应的进行，激波在x≈64 cm处首次从燃烧室壁面反射，紧接着的第二次反射发生在x≈128 cm处。这两个位置均与文献[9]基本一致。但由于文献[9]并未给出所采用的燃烧模型及流体物性设置，因此，本文采用的燃烧模型及流体黏度设置等可能与之不同，致使结果数值上出现一些不同，如激波的顶点比文献[9]中的要靠后等。

图4分别为在不考虑辐射和考虑辐射两种情况下的温度分布。其中，考虑辐射的温度分布图是在后处理时通过计算将温度的影响加入不考虑辐射情况下的温度分布图得到的。随燃烧反应的开始，燃烧室内温度迅速升高。由图4可明显看出，由于辐射的作用，原来集中于壁面的热量向燃烧室内部转移，由于燃气发射热辐射，图4(b)燃气最高温度降低了120 K，燃气在发射热辐射的同时，也吸收热辐射，使燃烧室内高温区域的范围远大于图4(a)的，表明考虑辐射的燃烧和流场参数更符合实际情况。

图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的摩尔分数分布，这三种物质的摩尔分数及分布反映了燃烧的程度及点火位置。由图6～8可知，点火发生在x≈64 cm处，即激波在燃烧室壁面反射附近。H2O的摩尔分数达23%，OH和H2O的分布情况表明燃烧在喷射器附近开始后，随燃气不断向燃烧室后段推进，且燃烧的程度比较令人满意。而H2的分布主要集中在喷射器附近，从另一方面表明燃烧程度较高。

图6 OH摩尔分数分布云图Fig.6 Cloud image of the mole fraction distribution of OH

图7 H2O摩尔分数分布云图Fig.7 Cloud image of the mole fraction distribution of H2O

图8 H2摩尔分数分布云图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]来讲，两段偏低，中间偏高；对流热流密度则趋势一致，但由于前面所述模型等方面的原因，具体数值上有些差别。本文辐射热流密度的趋势与文献[12]和文献[21]的计算结果也是一致的，结合之前对温度、压力等云图的分析比较，本文所用求解器和所得到的结果具有一定的可信度。由图9可知，壁面辐射热流密度的平均值为45.6 W/cm2，最大值约为55 W/cm2；这个结果虽然比文献[9]中的70 W/cm2要小，却远大于文献[8]中的24 W/cm2。文献[8]中结果较小可能由于其将燃烧室内气体假设为常物性。计算辐射热流密度时所设定的出口黑度、壁面吸收率及反射率、气体黏度等模型所需参数的不同，也是导致当前结果与文献[9]结果差异的原因。但是这些算例及曲线图均表明，对燃烧室内热辐射在壁面某些部位的作用应当给予充分的重视。

图10 壁面对流热流密度Fig.10 Wall convective heat flux density

图11给出了燃烧室内壁面对流热流密度和辐射热流密度的对比。由图可知，在壁面的不同位置处，辐射热流密度占对流热流密度的比值有很大区别。在x≈120 cm处，考虑计算偏差后，所得辐射热流密度占对流热流密度值的30%以上。图11进一步表明热辐射对燃烧室壁面热防护的重要性。

图11 壁面辐射热流密度和对流热流密度对比图Fig.11 Comparison between wall radiative heat flux density and wall convective heat flux density

3.2 耦合计算结果分析

解耦算法具有计算效率高、操作方便等优点，但燃烧/流动与热辐射本质上是耦合的过程，而由于超声速流场中大量激波和涡系等的存在，现有商业软件在超声速燃烧耦合方面计算能力有限，目前的国内外文献中也鲜有耦合计算的先例。为研究解耦计算与实际的不同情况，本文也对耦合计算开展了初步研究。图12给出了耦合计算的温度分布，与图4(b)对比可知，燃烧室后半段内温度分布更为均匀，燃气的最高温度与解耦计算时相比进一步降低了78 K，与不加辐射时相比更是下降了近200 K。壁面高温区面积明显缩小。比较图12与图4(b)中x=40～60 cm处点火区域的图像可知，耦合计算在燃烧室中心点火区域附近的温度有明显升高，该结论表明在scramjet的点火及稳焰研究中热辐射的作用值得深入探究。图13给出了耦合计算与解耦计算中壁面辐射热流密度的对比，该图表明耦合计算壁面辐射热流的分布趋向更平稳，所得燃烧室内温度分布更均匀。

图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给出的燃烧室内几种主要参数的云图与文献[9]的基本一致，本文模型和改进求解器的正确性得到验证。

2) 进行燃烧/流动和热辐射的解耦计算可知，考虑辐射后，原来集中于壁面处的热量向燃烧室内部扩散，且燃气最高温度下降达120 K (图4)；H2、OH和H2O的摩尔分数分布(见图6～8)表明，本文求解器得到的燃烧情况良好。壁面辐射热流密度最高可达55 W/cm2，且可占对流热流密度的40%以上，因此燃烧室内热辐射的影响不应该被忽略。

3) 由燃烧/流动与热辐射的耦合计算结果可知，燃烧室中段及后段温度分布更加均匀，燃气最高温度比不考虑辐射时相比进一步降低，降幅达到200 K左右；并且耦合计算中点火区域温度升高。将耦合计算的壁面辐射热流密度曲线与解耦计算所得曲线对比可知(见图13)，耦合计算中壁面辐射热流的分布趋于平稳。由此可知，在燃烧室中对燃烧/流动和热辐射的耦合计算有必要进行深入探讨。

[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] 刘兴洲. 超燃冲压发动机性能初步研究[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] 杨事民, 张建良. 超燃冲压发动机燃烧室流场数值模拟[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] 黄志伟, 何国强, 秦飞, 等. 支板增强混合超声速燃烧的大涡模拟研究[J]. 固体火箭技术, 2015, 38(5): 664-670.

[14] 巩伟杰, 唐硕, 李世珍. 高超声速飞行器气动加热三维数值分析方法研究[J]. 飞行力学, 2011, 29(2): 78-81.

[15] 鲍文, 张聪, 秦江, 等. 超燃冲压发动机主被动复合热防护系统方案设计思考[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] 魏鑫, 孙冰, 郑力铭. 超燃冲压发动机前体/进气道热结构分析[C]∥中国宇航学会2005年固体火箭推进年会第22届年会论文集(发动机分册). [S.l.: s.n.], 2005.

[18] 艾青, 夏新林, 孙凤贤. 壁面热特性对超声速燃烧室热环境的影响[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.