支板增强混合超声速燃烧的大涡模拟研究

黄志伟，何国强，秦 飞，魏祥庚

(西北工业大学 燃烧、热结构与内流场重点实验室，西安 710072)



支板增强混合超声速燃烧的大涡模拟研究

黄志伟，何国强，秦 飞，魏祥庚

(西北工业大学 燃烧、热结构与内流场重点实验室，西安 710072)

基于开放源代码软件OpenFOAM，建立了三维可压缩反应流动大涡模拟求解器，采用了PaSR亚格子燃烧模型和27步的氢气-空气反应动力学机理，开展了支板增强混合的超声速燃烧大涡模拟研究，对比了滑移和粘性2种不同壁面边界条件的影响。计算结果表明，不同截面上的平均轴向速度和温度与实验数据吻合良好，较好捕捉了超声速扩散燃烧的火焰空间发展过程。详细讨论了剪切层增长、发展和破碎对燃烧过程的影响，揭示了支板后旋涡脱落与燃烧过程的耦合作用，区分了支板下游亚声速区和超声速区内不同的掺混模式。利用化学爆炸模式分析方法，获取了爆炸化学过程及其特征时间尺度，得到了详细的火焰结构及其稳定机制。

超声速燃烧；化学爆炸模式分析；大涡模拟；PaSR燃烧模型；OpenFOAM

0 引言

双模态超燃冲压发动机和组合循环发动机是天地往返输运系统和临近空间飞行器的潜在高效动力系统。这类吸气式发动机工作过程中涉及诸多基础性科学问题。例如，高超声速和高温气体动力学，超声速条件下的点火与火焰稳定，复合防热材料与主动热防护，激波/附面层相互作用和吸热型碳氢燃料技术等[1]。在这些前沿技术挑战中，最基本的问题之一是燃烧室内超声速流动条件下的湍流混合与高效燃烧。随着计算机和数值计算研究的发展，具有较高时/空分辨率的湍流燃烧大涡模拟分析，成为超声速燃烧和吸气式发动机研究的有效手段[2]。

国内，清华大学的周建兴等采用RANS方法和单步反应模型，通过隐式耦合求解可压缩N-S方程，对氢气的超声速燃烧进行了三维数值模拟[3]。国防科技大学的范周琴等利用Flamelet燃烧模型和混合LES/RANS方法对支板构型超燃冲压发动机进行模拟，氢气动力学模型为9组分19方程，对比了三维与二维计算的区别[4]。中国科学院力学研究所的李晓鹏等分析了超声速燃烧中的特征尺度及其影响因素，讨论了Flamelet模型的适用性[5]。国外，Genin和Menon等使用LES-LEM和LES-EBU模型结合两步氢气动力学机理，对比研究了德国宇航院模型超燃冲压发动机燃烧室，验证了大涡模拟在超声速流动与燃烧方面的潜在适用性[6]。Fureby等利用OpenFOAM平台开展了超声速燃烧大涡模拟方面的研究[7-9]，通过与大量试验对比，验证了该软件在超燃冲压发动机设计和超声速燃烧机理研究方面的适用性。

国内外针对氢燃料超声速燃烧的大涡模拟，大多采用Flamelet或EBU亚格子燃烧模型结合单步或简化的动力学模型。本文基于近年来极具吸引力的开放源代码计算平台OpenFOAM，采用PaSR亚格子燃烧模型结合氢气-空气的详细动力学机理，分析了湍流与详细动力学的耦合作用，讨论了超声速燃烧流场火焰结构的演化，分析了壁面粘性对流场参数分布的影响。

1 数值方法与计算实现

1.1 OpenFOAM计算平台与数值方法

OpenFOAM是由伦敦帝国理工学院开发的求解连续介质力学问题的开放源代码软件平台[10]。其本质是基于面向对象程序设计的C++库，具有优异的可移植性，包含丰富的物理和数学模型，可根据需要进行开发扩展，与其它软件具有良好接口。近年来，已经在基础研究和工程领域得到了广泛的校验和应用[7-9]。

本文的数值模拟基于OpenFOAM平台，建立了适用于任意非结构化网格的可压缩反应流动大涡模拟求解器。对所有对流通量的重建采用了基于非线性vanLeer限制器的二阶精度TVD格式，对粘性和亚格子通量的空间离散化采用了中心差分格式。计算的最大CFL数不超过0.3，对应本文计算的物理时间步长为10-8s量级。当二阶统计矩收敛时，认为求解收敛。

1.2 大涡模拟反应流的控制方程

1.3 亚格子流动与燃烧模型

(1)

(2)

该燃烧模型同时考虑了Kolmogorov尺度与求解尺度的综合效应，避免了将湍流与火焰的解耦处理，较深刻地刻画了湍流-化学相互作用过程，并得到了DNS和实验数据的广泛校验，显示了对于高速流动燃烧计算较优的预测能力[7-9,13]。

1.4 化学爆炸模式分析方法

典型反应流的微分方程可表达成如下离散形式[14]：

(3)

其中，D/Dt为随流导数(物质导数);y为依赖变量矢量，如温度、组分浓度等,对于空间离散化的流动方程，不同网格点上的化学组分浓度对应着不同的y入口值;ω表示化学源项;s包含所有的非化学源项，如火焰中的扩散、搅拌反应器中的均质混合等。当地化学信息完全包含在化学源项ω的Jacobian矩阵Jω中，化学模式即定义为Jω的特征模式，即一个特征值及其相应的左、右特征向量。如果化学模式特征值的实部λexp为正，则可进一步定义其为化学爆炸模式(CEM)[15]。因此，当非化学反应源项s在方程(3)中可忽略时，CEM的存在表征了混合物的本质是爆炸性的，即化学反应速率沿着其特征向量的方向呈指数型增长趋势。

1.5 化学动力学模型

在湍流燃烧的大涡模拟中，选择恰当的化学反应机理是关键的问题之一。为了更深入地研究反应机理与湍流的耦合作用，本文使用了根据文献[16]中的氢气详细动力学机理得到的9组分、27步反应的氢气-空气化学动力学模型，在不过多增加计算量的同时，较好地反映氢气燃烧的主要动力学过程，为获得湍流火焰稳定机制提供了保证。表1给出了该动力学模型每个反应速率常数表达式中的参数。

表1 氢气9组分、27步反应机理的速率常数Table1 Rate constants for 9 species 27-step chemical kinetics of hydrogen (cm3 ·mol·s·cal·K)

2 计算构型与边界条件

图1给出了本文所模拟的三维模型超声速燃烧室尺寸示意图[7]。计算网格包含590万结构化的六面体单元，沿燃烧室上下壁面、支板壁面、支板前缘和尾迹区及剪切层区域进行了局部加密，以保证大涡模拟在这些区域内具有较好的网格分辨率。支板底端开有15个直径为1 mm的小孔，每个小孔间隔2.4 mm。由于在宽度方向具有对称性，同时为了节约计算成本，本文的计算区域只包含3个燃料喷孔。氢气-空气的等值比为0.034。

本文采用不同类型的壁面边界条件，计算了2种工况。在空气和燃料入口，由于流动是超声速或声速的，根据特征理论，所有的变量应用Dirichlet边界条件；在燃烧室出口，采用Neumann边界；在宽度方向使用周期性边界条件。

图1 模型超燃冲压发动机燃烧室示意图Fig.1 Schematic of the model scramjet combustor

在燃烧室的上下壁面及支板壁面上，第一个工况应用了无粘、绝热边界条件，第二个工况则应用了粘性壁面。

表2给出了空气和燃料入口的状态参数，来流空气的湍流强度为0.5%，氢气为5%。

表2 空气来流与氢气射流的入口条件Table2 Inflow conditions for air stream and hydrogen jet

3 结果与讨论

3.1 LES结果与实验数据的校验

图2给出了不同截面上轴向速度的分布曲线。

(a)x=0.120 m

(b)x=0.249 m

图2(a)的位置非常靠近支板底端，由于燃烧释热，使得回流区内压力升高，同时温度升高，导致气体粘性增大。因此，氢气自支板喷出后受到十分明显的减速作用。但由于氢气出口速度很高，导致射流核心区相比于外围气体仍具有略高的速度，因而形成双峰状的速度分布。大涡模拟的结果在回流区内总体偏差较大，Menon和Fureby等指出，由于该位置恰好处于高度湍流的尾迹区，实验中很难对速度进行精确地测量[6-7]，是造成偏差较大的原因之一。图2(b)在距离支板很远的下游位置，大尺度旋涡结构逐步转化为小尺度漩涡，在燃烧释热膨胀和外层空气流剪切力的双重作用下，核心流逐步加速至接近主流的速度，同一截面上速度分布较为平缓。

图3给出了不同位置上时间平均的温度与实验数据的对比曲线。图3(a)中，反应较为剧烈，燃烧比较充分，释放出大量的热，形成高温燃烧区，射流中心温度达到局部最大值。图3(b)中，由于氢气和中间活性物质逐渐消耗完毕，反应强度减弱，流场中心温度也有所降低。同时，由于分子的湍流输运和辐射传热等作用，流场高温区向上下壁面扩展，表现为曲线在高度方向上拓宽而峰值降低。无粘壁面的结果在第二个位置气流核心区内偏差较大，可能由于假设了支板壁面无粘而导致燃料与空气混合更充分，因而在燃烧室后段仍具有较强的热释放。

(a)x=0.167 m

(b)x=0.275 m

通过对比显示，本文的大涡模拟计算结果与实验数据测量结果吻合较好，有效地捕捉到了超声速/亚声速混合燃烧流场的演化过程。粘性壁面得到的速度和温度分布准确性高于无粘壁面假设的结果。

3.2 无粘壁面与粘性壁面的结果对比

为了研究不同壁面边界条件对流场结构发展、两股气流掺混及燃烧过程的影响，定义了标量场λ，它表征了流场不同的流动类型[17]：

(4)

其中

(5)

当λ=-1时，表示有旋流动；λ=0表示简单剪切流动；λ=1表示平面拉伸流动。图4(a)和(b)分别给出了粘性和无粘壁面条件下，燃烧室中心截面上λ的分布。考虑壁面粘性时，附面层内仍保持为简单剪切流动，而附面层外迅速转变成平面拉伸流动。同时，由于流体往下游发展，并不存在完全意义上的“纯旋转流动”，流体发生旋转运动的同时，总伴随着往下游的平动，因而并不存在λ=-1的情况。支板后的简单剪切流动区域加宽，往下游发展时，表现出明显向外扩展的趋势。图4(a)与图4(b)的对比，清楚地表明了壁面粘性效应的影响范围十分有限，当主要关心核心流的发展时，假设壁面无粘具有一定的合理性。

图4 粘性与无粘壁面下λ分布的对比Fig.4 Distribution of λ for comparison of viscous and non-viscous wall

3.3 支板掺混增强与火焰稳定的机理分析

图5给出了标量耗散率(Scalar Dissipation Rate)的等值面图和中间活性物质OH质量分数的切面，以区分支板下游亚声速和超声速不同区域内的混合与燃烧模式。标量耗散率χ表征了燃料与氧化剂的混合速率，而OH基团可表征火焰面位置[18]。在近支板底端区域及初始反应混合层外边界，由于两股气流之间存在较大参数梯度，混合主要通过扩散实现；而在更远的下游，大尺度旋涡结构的卷吸与拉伸对掺混过程起决定性作用。在燃烧室末端，燃料与氧化剂已经达到了充分的混合，并实现了火焰的保持。从图5可看到，该部分存在着大量的OH自由基。燃烧首先发生在支板尾迹区较薄的剪切层内，该区域聚集了大量的富燃燃气和活性自由基，并对空气和燃料进行预加热，从而有利于反应剪切层内的燃烧与火焰稳定。

图5 标量耗散率等值面与OH质量分数切面图Fig.5 An iso-surface of scalar dissipation rate and contours of OH

来流空气与氢气射流作用形成反应剪切层，发生剧烈的湍流燃烧，并实现质量、动量和能量的交换。图6给出了表征旋涡强度的速度梯度张量第二不变量 (Second Invariant of Velocity Gradient Tensor)等值面图[19]。可看到，剪切层内含有丰富的旋涡结构。支板底部射流-尾迹区内大的旋涡结构往下游发展时，破碎成很多小旋涡。一方面，促进了反应物之间的掺混及随后的燃烧；另一方面，又给热释放带来较大的扰动，使火焰在三维空间内产生褶皱，并发生扭曲。支板前缘形成的两道压缩波以及支板末端形成的两道膨胀波与燃烧室的上下壁面相互作用，形成复杂的反射波结构。由于超声速反应剪切层较为稳定，且厚度较大，反射波无法直接穿透中心射流，从而在壁面与剪切层之间形成新的反射波系。燃烧释热还增加了回流区内的压力，使得剪切层开始向外扩展，回流区扩大，有利于火焰的保持和稳定。

图6 速度梯度张量第二不变量等值面与温度切面图Fig.6 An iso-surface of second invariant of velocity gradient tensor and contours of temperature

为进一步研究燃烧室内的湍流-化学相互作用过程，图7给出了燃烧产物H2O作为混合物分数函数的散点图分布，揭示燃烧室不同位置上掺混好坏与反应进行程度之间的关系。可看到，H2O的质量分数与混合分数随着位置的后移表现出愈加明显的线性关系，且分布更趋于集中化。在x=0.167 m位置处，相对较为分散的分布表明，该位置主要的反应活性强烈地受到湍流的影响，空气与燃料通过大的湍流脉动强制进入到混合区的反应混合层内。掺混越充分，相应的燃烧反应越剧烈，产物H2O的质量分数也越大，但由于湍流脉动自身的随机性，导致水的质量分数与混合物分数之间并不是严格的线性函数。而在更远的下游位置，如x=0.199 m和x=0.275 m，随着大尺度旋涡结构的破碎及湍动能的耗散，湍流的作用相对减弱，随着位置的后移趋势更加明显。对于H2-O2燃烧系统，当量比混合物分数为Zst=0.028 5，在燃烧室后段，混合物分数已经低于该值，说明主燃烧反应区结束，在支板下游，实现了高效燃烧。

(a)x=0.167 m

(b)x=0.199 m

(c)x=0.275 m

图8给出了放热率作为混合物分数函数的散点图分布。在x=0.167 m处，随着混合物分数增加，放热率首先缓慢上升，超过某一临界值后迅速下降。更高的混合物分数意味着更充分的掺混，进而预示着更有效的燃烧。恰当的混合对于高效燃烧是必要的，但过大的混合速率会导致局部热量和活性基团流失，化学反应反而受到抑制，甚至导致火焰不能自持稳定。在x=0.180 m处，放热率随混合物分数增加只保持缓慢增加，且分布极为分散，表明该位置放热已趋于饱和，掺混对放热率不起决定性作用。放热率表现出强烈的脉动表明该位置的燃烧反应主要受湍流脉动影响。在x=0.199 m处，放热率普遍较低，主放热区结束。较低的混合物分数表明，H2已基本消耗完毕，中间活性物质对该位置的化学反应进程起主导作用。

(a)x=0.167 m

(b)x=0.180 m

(c)x=0.199 m

3.4 化学爆炸模式分析

利用化学爆炸模式分析方法，对每个网格节点进行了处理，图9给出了化学源项的特征值，即爆炸模式时间尺度倒数的分布(对数坐标)。由于主流的流通时间不超过1 ms量级，故对时间尺度大于1 s的模式进行了截断(该时间尺度相比于流通时间几乎为“静止的”，认为是非爆炸性的)。从图9中可看到，整个燃烧室可划分成两个不同类型的区域：即蓝色背景的非爆炸区域和彩色标识的核心流爆炸区域。其中，非爆炸区域又包含3个子区域：氢气核心射流外围的“贫燃”来流空气流动区，氢气核心射流紧邻出口下游的低温“富燃”湍流区，以及富含平衡燃烧产物的上、下混合层区域内。在前两个子区域内主要是因为混合物超过了其贫/富可燃极限。

在近支板底端区域，由于主要通过扩散来实现掺混，因此在非爆炸和爆炸混合物之间存在较为显著的界限，CEM特征时间尺度表现出较大的梯度。与图6中相对低温的高温混合物相对应的两个混合层外边界将爆炸混合物与非爆炸混合物分开，即在该层内形成了所谓的预混燃烧模式。在燃烧室的下游位置，主要通过大尺度旋涡的卷曲拉伸来实现混合。相对而言，爆炸/非爆炸混合物的混合更为均匀，二者间形成较厚的中间过渡剪切层。

图9 当地化学爆炸模式时间尺度lgλexp的空间分布Fig.9 Spatial distribution of the time scale of local chemical explosive mode

4 结论

(1)跨越支板回流区后大尺度旋涡的破碎对燃烧过程具有双重作用：一方面，由于形成大量的小旋涡结构，有效促进了燃料-氧化剂的掺混，进而促进了化学反应进程；另一方面，也给热释放率带来很大扰动。PaSR亚格子燃烧模型能捕获超声速燃烧流场的复杂物理化学过程。

(2)在近支板底端区域内的掺混主要通过扩散来实现；而在较远的下游位置，湍流混合主要是由大尺度旋涡结构的卷吸与拉伸作用来完成。通过化学爆炸模式分析，获得了燃烧室火焰结构及其稳定机理。

(3)在支板下游不同位置处，放热速率存在不同的影响因素：随着离开支板轴向距离的增加，首先是燃料与氧化剂的掺混，然后是湍流脉动，再者是中间活性物质，分别对放热率起主导作用。

(4)粘性壁面相比于无粘壁面最主要的影响在壁面附面层很窄的区域内，对核心流场的发展影响很小。在主要关心核心流的参数变化时，假定壁面无粘具有一定的合理性。

(5)大涡模拟的平均轴向速度和温度均与试验数据吻合较好，表明OpenFOAM求解器有较高的数值计算精度。其开放源代码的特性为后续开展燃烧流场计算与发动机设计、多学科优化和不确定性分析方法的耦合提供了充分的拓展空间。

[1] 蔡国飚, 徐大军. 高超声速飞行器技术[M]. 北京:科学出版社, 2012.

[2] Clark Ryan J, Shrestha SO Bade. A review of numerical simulation and modeling of combustion in scramjets[J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2014:1-23.

[3] 周建兴, 朴英, 岂兴明, 等. 超声速燃烧室内氢气燃烧的三维数值研究[J]. 航空动力学报, 2007, 22(12):1984-1988.

[4] 范周琴, 孙明波, 刘卫东. 支板喷射超声速燃烧流场三维大涡模拟[J]. 国防科技大学学报, 2011, 33(1):1-6.

[5] 李晓鹏, 张泰昌, 齐力, 等. 超声速燃烧中的特征尺度及影响因素[J]. 航空动力学报, 2013, 28(7):1458-1466.

[6] Genin Franklin, Menon Suresh. Simulation of turbulent mixing behind a strut injector in supersonic flow[J]. AIAA Journal, 2009, 48(3):526-539

[7] Fureby C, Fedina E, Teger J. A computational study of supersonic combustion behind a wedge-shaped flameholder[J]. Shock Waves, 2014, 24(1):41-50.

[8] Fureby C, Chapuis M, Fedina E, et al. CFD analysis of the HyShot II scramjet combustor[J]. Proceedings of the Combustion Institute, 2011, 33(2):2399-2405.

[9] Berglund M, Fedina E, Fureby C, et al. Finite rate chemistry large-eddy simulation of self-ignition in a supersonic combustion ramjet[J]. AIAA Journal, 2010, 48(3):540-550.

[10] www.openfoam.com[DB/OL].

[11] Poinsot T J, Veynante D. Theoretical and Numerical Combustion, Third Edition[M]. US: R. T. Edwards, 2012.

[12] Bensow R E, Fureby C. On the justification and extension of mixed models in LES[J]. Journal of Turbulence, 2008, 54(8):1-17.

[13] Sabelnikov V, Fureby C. Extended LES-PaSR model for simulation of turbulent combustion[J]. Progress in Propulsion Physics, 2013, 4(1):539-568.

[14] Luo Z Y, Yoo C S, Richardson E S, et al. Chemical explosive mode analysis for a turbulent lifted ethylene jet flame in highly-heated coflow[J]. Combustion and Flame, 2012, 159(1):265-274.

[15] Lu T F, Yoo C S, Chen J H, et al. Three-dimensional direct numerical simulation of a turbulent lifted hydrogen jet flame in heated coflow: a chemical explosive mode analysis[J]. Journal of Fluid Mechanics, 2010, 652(1):45-64.

[16] Marinov N M, Westbrook C K, Pitz W J. Detailed and global chemical kinetics model for hydrogen[C]//8th International Symposium on Transport Processes, 1995.

[17] http://openfoamwiki.net/index.php/FlowType[DB/OL].

[18] Hawkes, Evatt R, Sankaran Ramanan, et al. Direct numerical simulation of turbulent combustion: fundamental insights towards predictive models[J]. Journal of Physics: Conference Series, 2005, 16(1):65-79.

[19] Khashehchi M, Elsinga G E, Bornstein H, et al. Studying the invariants of the velocity-gradient tensor of a round turbulent jet using Tomo-PIV[C]//8th International Symposium on Particle Image Velocimetry, 2009.

(编辑：崔贤彬)

Large eddy simulation of strut enhanced mixing for supersonic combustion

HUANG Zhi-wei, HE Guo-qiang, QIN Fei, WEI Xiang-geng

(Science and Technology on Combustion, Internal Flow and Thermal-Structure Laboratory, Northwestern Polytechnical University, Xi'an 710072, China)

Large Eddy Simulation (LES)of supersonic combustion in a model scramjet combustor based on an Open Source Field Operation and Manipulation (OpenOAM)computing platform was established, with two different wall boundary conditions, i.e. slip and viscous walls applied. The three-dimensional LES solver, which adopts a Partially Stirred Reactor (PaSR)sub-grid combustion model along with a skeleton 27 steps hydrogen chemical kinetics, was used to study strut-enhanced mixing and combustion. LES results show that mean axial velocity and temperature at different cross sections match well with experimental data, and spatial evolution of the supersonic diffusion flame is well captured. Effects of shear layers growth, development and breaking down on combustion processes were discussed in detail, and the coupling effects with vortex shedding at the strut base were revealed. Different mixing modes were recognized after the strut where subsonic and supersonic flows co-exist. Explosive chemical processes and their characteristic time scales were acquired by the Chemical Explosive Mode Analysis (CEMA)method, and the detailed structure and stabilization mechanism of the flame was identified.

supersonic combustion；chemical explosive mode；large eddy simulation；PaSR combustion model；OpenFOAM

2014-10-17；

：2015-03-26。

黄志伟(1989—)，男，博士生，研究方向为航空宇航推进理论与工程。E-mail:huangzhiwei504@mail.nwpu. edu.cn

V435

A

1006-2793(2015)05-0664-07

10.7673/j.issn.1006-2793.2015.05.012