湍流燃烧反应数值模拟的研究与实现

孙兆国，刘志勤，刘 涛，刘敏贤，吴颖川

（1.西南科技大学 计算机科学与技术学院，四川 绵阳621010；2.中国空气动力研究与发展中心，四川绵阳621000）

0 引 言

燃烧流动的数值模拟是指应用计算机为工具，将流体力学、传热学、化学反应动力学和数值计算方法相结合所得到的求解化学流体力学基本方程的理论和方法。数值模拟的基本思想是：把空间与时间坐标中连续的物理量的场（如速度场，温度场，浓度场等），用一系列有限个离散点上的值的集合来代替，通过一定的原则建立起这些离散点上变量值之间关系的代数方程，并根据所建立起来的代数方程求解这些变量的近似值。数值模拟的研究方法对燃烧室性能评估有一定的参考价值。

本文的研究重点是基于RANS［1］求解的定长可压缩湍流燃烧流动的数值模拟。湍流是工程技术领域中常见的流动形式，在湍流流动中流体的各种物理参数，如速度、压力、温度等都随时间与空间发生随机变化。流体做湍流燃烧时涉及到流体流动、对流换热以及化学反应等过程，数值模拟比较繁琐。在燃烧流动数值模拟中需要借助计算流体力学平台。当前计算流体力学数值模拟领域多是商业软件，如Fluent、CFX等。由于商业机密等原因，其源代码难以开放，当计算结果不理想时，针对特定问题的判断和解决较困难。本文燃烧流动数值模拟采用C＋＋语言开发的开源计算流体力学软件OpenFOAM［2］，这对于开发新的针对特定问题的解算器有重要的帮助作用。

在OpenFOAM平台上做燃烧数值模拟的研究尚处于初级阶段，Aalborg University的 Christian Andersen等［3］使用OpenFOAM中的dieselFoam解算器对基于Burner Flow Reactor（BFR）的柴油喷射进行数值模拟，中国科学技术大学的陆阳、武文等［4－6］改进了OpenFOAM平台自带的基于RANS求解的非定常燃烧流动解算器reactingFoam对湍流燃烧进行数值模拟，这些都是对OpenFOAM中自带的解算器进行改进且是对非定常可压缩燃烧过程做数值模拟。

本文根据定常可压缩湍流燃烧流动机理和化学动力反应模型库Cantera设计出湍流燃烧流动的解算器，然后使用该解算器对Sydney钝体驻定火焰［7－8］进行数值模拟。模拟得到了燃烧流动组分浓度分布图和温度曲线变化图等。通过将模拟的温度变化与Sydney钝体驻定火焰的实验数据进行对比验证表明：模拟效果满足燃烧组分变化的研究要求，流场温度较好的符合Sydney钝体驻定火焰HM1的实验数据，这验证了解算器的合理性和可行性。

1 数学模型

在湍流燃烧反应中，反应物湍流流动过程和化学反应过程相互关联又相互影响，湍流流动组分浓度及温度的脉动可以强化组分的混合与传热，因而影响时平均化学反应速率，化学反应放热过程迅速放热而引起密度变化，同时使流体输运系数变化，进而影响着湍流流动过程。这种相互关联和相互影响使得湍流燃烧反应极为复杂。湍流燃烧流动数值模拟的基本思想是分别独立描述湍流流动和化学反应过程，然后考虑湍流流动和化学反应相互作用，基于这个思想，在数值模拟过程中使用的燃烧模型包括湍流模型，湍流燃烧相互作用模型和化学反应模型。

湍流模型主要包括层流（laminar），直接数值模拟（direct numerical simulation，DNS），大涡模拟（large eddy simulation，LES）和雷诺时均数值模拟（RANS）［9］。雷诺时均数值模拟（RANS）是将瞬态N—S方程的瞬时量分解为时均值和脉动值之和（雷诺分解），再取时间平均，得到雷诺时均方程，然后利用某些模拟假设，将方程中的高阶的未知关联项用低阶项或时均量来表达，从而使雷诺（reynolds）时均方程封闭。

本文中使用的模型湍流是基于RANS的k－ε双方程模型和Cantera化学动力反应模型库。

1.1 湍流模型

在反应物进行湍流混合时，使用k－ε双方程模型。kε双方程模型对于一般的流动情况，如平壁边界层、无力平面射流、管流等流动双方程模型能给出相当满意的计算结果，而且计算工作量小。k－ε双方程模型［10］，模型如下

在上述方程中，Gk表示由于平均速度梯度引起的湍动能产生，Gb表示用于浮力影响引起的湍动能产生，YM表示可压速湍流脉动膨胀对总的耗散率的影响。

k－ε双方程模型的湍动能k和耗散率ε方程求解式为式中：Ux′，Uy′，Uz′——x，y，z这3个方向的速度矢量，l——特征长度，Cμ设定——0.09。湍流粘性系数求解式为

1.2 化学反应模型

化学反应模型使用Cantera化学动力反应模型库［11］。Cantera是跨平台的面向对象的开源软件，主要应用于化学反应，热力学和传输过程等问题的求解。Cantera化学反应模型是基于美国气体研究所和加州大学的研究成果GRIMech3.0详细化学反应机理［12］。GRI－Mech3.0详细化学反应机理总共涉及到53种组分和325个基元反应，其中包括C1反应、C2反应、C3反应和NOX的生成机理。

对于燃料CmHn，其和氧气的总包反应可以写成［13］

对于CO，其和氧气的总包反应可以写成

2 仿真平台的建立

2.1 OpenFOAM平台

本文解算器设计使用近几年流行起来的计算流体力学软件OpenFOAM作为平台。OpenFOAM的全称是Open Field Operation and Manipulation，它的核心是使用 C＋＋编写的面向对象的计算流体力学类库，有较高效率的偏微分方程求解模块。OpenFOAM是一款开源软件，其源代码是对外公开的，这就可以方便用户编写解算器。Open－FOAM软件可以模拟复杂的流体流动、化学反应、传热传质等问题，还可以进行结构动力学分析、电磁场分析等的数值仿真。

从程序实现功能的角度来看，OpenFOAM软件包括核心解算器、前处理和后处理3大模块。前处理主要用网格设计的工具来进行，如OpenFOAM本身自带的BlockMesh工具，也可以使用其他商业网格设计软件，如GridGen等。在本文中，燃烧室网格设计采用GridGen，然后转换成OpenFOAM可以识别的网格文件。解算器主要是求解计算流体力学问题常用的一些解算器，OpenFOAM的解算器支持的流动模型有层流（laminar），基于雷诺时均湍流模型（RANS），大涡模拟（LES）以及直接数值模拟（DNS）等。后处理主要是将计算结果转化成可视化图像，如云图，燃烧等值面图等。OpenFOAM的组件构成如图1所示。

图1 OpenFOAM 的组件构成［14－15］

2.2 解算器

解算器位于/opt/application/solvers/combustion/turbulentFoam目录下。解算器包括Make目录，头文件和主程序3部分。其中Make目录包括options文件和files文件，options文件定义编译选项，用于指定编译用到的头文件位置及其动态库，files文件用于指定当前要编译的文件。头文件是主程序编译过程中调用用到的程序，解算器主程序是turbulentFoam.C。解算器结构如图2所示。

图2 解算器结构

3 算例

3.1 物理模型

计算区域如图3所示，钝体半径为50mm，燃料入口半径为3.6mm。钝体网格结构简图如图4所示，使用结构化网格离散计算区域，网格数为5580，并在轴向划分5°，以便在OpenFOAM中模拟轴对称问题。

网格有上下两个进气口，上口为伴随流，注入成分为空气，为了简化，使用3组分空气构成，分别为N2，O2和Ar，伴随流的温度为300K。射流入口注入燃料煤气，煤气中包含CH4，CO，H2，CO2，H2O和N2。反应工况如表1所示。反应开始时，钝体内充满空气。燃烧反应的边界条件见表1。

表1 燃烧反应的边界条件和初始条件

3.2 离散方法

对控制方程采用有限体积法进行离散，对控制方程中的瞬态项采用LU－SGS隐式迭代法求解，对流项差分采用二阶迎风的高斯差分格式，扩散项差分采用二阶线性修正的高斯差分格式。由于是定常可压缩流动计算，本文采用SIMPLE（semi－implicit method for pressure－linked equations）［16］半隐式压力校正方程算法，SIMPLE算法是流体流场数值模拟中最重要的模拟方法。SIMPLE算法的求解流程如图5所示。

图5 SIMPLE算法求解流程

3.3 求解过程及分析

计算平台使用基于MPI的并行计算系统［17］。并行计算原理即在计算前将燃烧室网格剖分为4个子区域分配给各CPU内核，把子区域的初始相关信息分别装载入各子区域对应的CPU内存，每一个CPU内核分别启动一个计算进程，这样同时有四个进程并行运算，可以加快计算时间。

在数值模拟阶段，将Coalgas和Air分别从fuelInlet和airInlet以118m/s和40m/s的速度注入燃烧室。反应物注入后在燃烧室内进行燃烧反应，经过21600s的计算，燃烧反应稳定。图6，图7和图8分别为燃烧反应稳定后产物二氧化碳，水，速度分布图。从图中可以看出燃烧反应从喷口到狂涨段都富含燃烧反应产物二氧化碳和水，并在输出口富集。分布云图较好的反映了反应流场在反应稳定时的一些基本特征。

图8 速度分布（m/s）

燃烧反应的参数温度是与实验结果对比的重要参数。图9是燃烧反应稳定时燃烧室截面中取一条平行于X轴的端点坐标为（0，0）和（0，0.4）的直线并根据此直线上坐标对应的温度得到的温度分布曲线图，从曲线图可以得出最高温度超过2200K。

Sydney钝体驻定火焰是国际湍流扩散火焰测量与计算研讨会（TNF workshop）的标准系列火焰之一，常用于与仿真结果做对比。实验温度分布与Sydney钝体驻定火焰的HM1实验结果基本吻合，尤其是燃烧的最高温度等变化情况的计算结果很好。

图9 燃烧反应稳定时温度分布（T/K）

从图6，图7和图8可以看出，火焰上游没有出现由燃料射流扩散而形成的回流区，这使得模拟出的效果偏向层流燃烧，因此在数值模拟时，k－ε双方程模型中k和ε的设置还得进一步探讨，以便得到更好的模拟效果。另外，由于在数值模拟时考虑涉及详细化学反应机理，求解的规模巨大，为了加快计算时间，将网格规模设计较小，这使得有些区域网格不够密集，尤其是火焰上游区域。这对数值模拟精度也有一定的影响。

4 结束语

本文在开源计算流体力学平台OpenFOAM和Cantera化学动力反应模型库的基础上，设计出基于定常可压缩的湍流燃烧流动反应解算器，并通过该解算器对Sydney钝体驻定湍流扩散火焰HM1进行数值模拟，在Coalas/Air详细反应机理的条件下得到的计算结果并得出了燃烧流场中各组分浓度和温度变化曲线图，仿真效果较好的符合燃烧组分变化的研究要求，尤其是温度曲线变化效果比较理想。模拟效果的具体情况基本验证了湍流燃烧流动解算器设计的合理性，同时，对OpenFOAM平台和Cantera反应库的研究对于将开源软件应用在燃烧流动数值模拟领域有重要的开拓意义。

［1］ZHANG Zhaoshun，CUI Guixiang，XU Chunxiao.Theory and modeling of turbulence ［M］.Beijing：Tsinghua University Press，2005（in Chinese）.［张兆顺，崔桂香，许春晓.湍流理论与模拟 ［M］.北京：清华大学出版社，2005.］

［2］Hrvoje Jasak.OpenFOAM：Open source CFD in research and industry ［J］.International Journal of Naval Architecture and Ocean Engineering，，2009，1（2）：89－94.

［3］Niels E L Nielsen.Numerical investigation of a BFR using OpenFOAM ［R］.Fluids and Combustion Engineering，2008.

［4］LU Yang.A study on flamelet model in combustion simulation［D］.Hefei：University of Science and Technology of China，2009：93－117（in Chinese）.［陆阳.燃烧计算中火焰面模型的研究［D］.合肥：中国科学技术大学博士论文，2009：93－117.］

［5］WU Wen，HUANG Wei，ZHAO Pinghui，et al.Numerical simulation of turbulent CH4/Air of jet diffusion ［C］.Hefei：The Chinese Society of Engineering Thermophysics Conference，2009（in Chinese）.［武文，黄威，赵平辉，等.湍流CH4/Air射流扩散燃烧数值模拟 ［C］.合肥：中国工程热物理学会学术会议论文，2009.］

［6］WU Wen，HUANG Wei，ZHAO Pinghui，et al.Study of RANS based flamelet/progress variable turbulent combustion model ［J］.Journal of University of Science and Technology of China，2010，40（10）：1016－1022（in Chinese）.［武文，黄威，赵平辉，等.基于RANS求解的火焰面/反应进度变量湍流燃烧模型研究 ［J］.中国科学技术大学学报，2010，40（10）：1016－1022.］

［7］Dally B B，Masri A R，Barlow R S，et al.Instantaneous and mean compositional steucture of bluff body stabilized non－premixed flames［J］.Combustion F1ame，1998，114（1－2）：119－148.

［8］HUANG Qing，ZHU Minming，YE Taohong，et al.PDF simulation of bluff－body stabilized turbulent non－premixed flame［J］.Chinese Journal of Comutational Physics，2008，25（6）：733－743（in Chinese）.［黄庆，朱旻明，叶桃红，等.PDF方法模拟钝体驻定的湍流扩散火焰 ［J］.计算物理，2008，25（6）：733－743.］

［9］LU Yu.Kinetic reduction for combustion chemistry and investigation of methane turbulent jet flame using direct numerical simulation［D］.Hangzhou：Zhejiang University，2011（in Chinese）.［吕钰.燃烧化学机理简化及甲烷湍流射流火焰的直接数值模拟研究［D］.杭州：浙江大学硕士论文，2011.］

［10］SONG Baojun.On tow－equation turbulence model for complex turbulent flow by using DI approximation ［C］.Kyoto，Japan：The 9th International Conference of Turbulent Shear Flow，1993.

［11］Moffat H K.CADS：Cantera aerosol dynamics simulator ［M］.Albuquerque：Sandia National Laboratories，2007：7－9.

［12］Michael Frenklach.GRI－Mech30 ［EB/OL］.http：//www.me.berkeley.edu/gri＿mech/version30/text30.html，2011.

［13］WANG Hui，HOU Lingyun.Numerical simulation of hydrocarbon gas combustion and supersonic shear flow ［J］.Journal of Aerospace Power，2007，22（4）：559－564（in Chinese）.［王慧，侯凌云.碳氢燃气超声速剪切流动燃烧数值模拟［J］.航空动力学报，2007，22（4）：559－564.］

［14］OpenFOAM Company.OpenFOAM user guide ［M］.Berkshire：OpenCFD Company，2010.

［15］Mattijs janssens andy heather sergio ferraris graham macpherson helene blanchnet henry weller CFD toolbox－user guide［M］.Berkshire：OpenCFD Company，2010.

［16］WANG Fujun.The analysis of computational fluid dynamics－principles and applications of CFD software ［M］.Beijing：Tsinghua University Press，2005：74－75（in Chinese）.［王福军.计算流体动力学分析－CFD软件原理与应用 ［M］.北京：清华大学出版社，2005：74－75.］

［17］Edgar Gabriel.OpenMPI：Open source high performance computing ［EB/OL］.http：//www.open－mpi.org/，2011.