湍流预混射流火焰直接模拟中入口条件的确定
吴韶华, 张健
清华大学 工程力学系, 北京 100084
张健, 教授, E-mail:jianzhang@mail.tsinghua.edu.cn

作者简介: 吴韶华(1982—), 男(汉), 河南, 博士研究生。

摘要

该文对甲烷-空气湍流平面射流预混火焰进行了直接数值模拟。射流入口的扰动速度数据依据给定的湍流能谱生成,利用入口处的湍流积分长度尺度和湍动能确定了能谱的峰值波数。计算结果给出了气体温度、质量分数和涡量模的瞬态分布,表明在剪切层内随着旋涡尺度的增大出现了拟序结构。化学反应受到湍流的作用,瞬时反应面出现了明显的皱折,反应面积增大。沿射流中心线,湍动能逐渐衰减,温度脉动和甲烷质量分数脉动均方根值则逐渐增大。

关键词: 湍流平面射流火焰; 预混火焰; 直接数值模拟; 入口条件; 速度扰动
中图分类号:TK16 文献标志码:A 文章编号:1000-0054(2014)06-0834-05
Inlet conditions for direct numerical simulations of turbulent premixed jet flames
Shaohua WU, Jian ZHANG
Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China
Abstract

Methane-air turbulent premixed planar jet flames are simulated using direct numerical simulations. A fluctuating velocity field at the jet inflow boundary is generated based on a prescribed turbulent energy spectrum. The peak wavenumber in the energy spectrum is determined from the turbulent integral length scale and turbulent kinetic energy at the jet inlet. The model also gives the instantaneous distributions of the gas temperature, species concentrations, and vorticity. The results show that coherent structures in the shear layer gradually appear as the eddy sizes increase. The chemical reactions are affected by the turbulence and the instantaneous reaction surface is quite wrinkled with its area increase. The turbulent kinetic energy gradually decreases, while the root mean squares of the temperature and methane concentration fluctuations increase along the jet centerline.

Keyword: turbulent planar jet flame; premixed flame; direct numerical simulation; inlet condition; fluctuating velocity

湍流燃烧的直接数值模拟(direct numerical simulation, DNS)可揭示出火焰的瞬态多尺度结构,为深入认识湍流燃烧的特点和规律提供了大量丰富的信息,因而受到了越来越多的重视。在湍流燃烧的直接模拟中,入口扰动数据的生成是关键和难点。较简单的方法是在入口平均速度分布上添加随机扰动[1]。较复杂的方式是将计算区域扩展到入口的上游,或者预先进行入口上游区域的DNS计算,以获取入口扰动数据[2]。还有一种方法是采用数字过滤的方式,生成满足给定的二阶统计矩与局部自相关函数的入口扰动数据[1,2,3]。而另外一种常用的方法是依据真实的湍流能谱生成入口扰动数据[4,5,6,7]。为了确定能谱的形状,该方法需要指定能谱峰值对应的波数。

然而,在入口平均速度分布上添加随机扰动的方法在高波数区引入了过多的能量,可能会导致扰动快速地衰减。将计算区伸展至入口上游的方法会增加额外的计算量,给定局部自相关函数也有一定的难度。因此,本文采用了按照各向同性衰变湍流的能谱生成入口扰动数据的方法,通过入口处的湍流积分长度尺度和湍流动能确定能谱峰值对应的波数,实现了对甲烷-空气湍流预混射流火焰的直接模拟。模拟得到了火焰的瞬时温度场、质量分数场和速度场,并进一步得到了火焰的统计平均特性。

1 控制方程组与数值模拟方法

可压缩湍流反应流的瞬时控制方程组包括连续、动量、能量与组分质量分数方程,它们分别为

ρt+xi(ρui)=0,(1)ρuit+xj(ρuiuj)+pxi=τijxj+ρk=1NsYkfk,i,(2)ett+xi(et+p)ui=xj(uiτij)-qixi+ρk=1NsYkfk,i(ui+Vk,i)+qR,(3)(ρYk)t+xiρYk(ui+Vk,i)=-Wk.(4)

其中:

et=12ρui2+ρh-p,(5)qi=-λTxi+ρk=1NshkYkVk,i,(6)τij=μ(uixj+ujxi)-23μukxkδij,(7)YkVk,i=-DkYkxi.(8)

压力通过理想气体状态方程求出,

p=ρRTk=1NsYkMk.(9)

以上各式中: et代表气体总能; qi τij分别为导热通量与粘性应力; Yk, hk, Vk, i Wk分别为气体第 k组分的质量分数、焓、扩散速度分量和反应率; fk, i为作用在 k组分上的体积力分量。

为便于求解,对上述方程组作了如下的简化: 忽略热辐射源项 qR; 各个组分的Lewis数为常数且都等于1; 各个组分的比热容为常数且相等; 粘性系数的计算采用Sutherland公式; Prandtl数为常数,并取为0.75。

上述瞬时控制方程组的求解采用可压缩流动的数值求解方法,将密度和速度作为待求变量,分别通过连续方程和动量方程求出。在各方程的离散化中,空间离散采用6阶精度的中心有限差分格式。以 x方向的导数为例,在网格均分的条件下, 1阶和2阶空间导数的离散格式分别为:

ϕxi,j=34ϕi+1,j-ϕi-1,jΔx-320ϕi+2,j-ϕi-2,jΔx+160ϕi+3,j-ϕi-3,jΔx,(10)2ϕx2i,j=13600Δx2-4214ϕi,j+828(ϕi+1,j+ϕi-1,j)+1935(ϕi+2,j+ϕi-2,j)-810(ϕi+3,j+ϕi-3,j)+171(ϕi+4,j+ϕi-4,j)-18(ϕi+5,j+ϕi-5,j)+ϕi+6,j+ϕi-6,j.(11)

其中, Δ x x方向的网格尺寸。利用Fourier变换进行分析可以看到,上述有限差分格式的误差随波数的增加而增大。应用8阶精度的空间过滤来抑制离散格式导致的这类误差[8],变量 ϕ的空间过滤格式为

ϕ^i=ϕi+1256-70ϕi+56(ϕi+1+ϕi-1)-28(ϕi+2+ϕi-2)+8(ϕi+3+ϕi-3)-(ϕi+4+ϕi-4).(12)

过滤格式采用的结点都在同一个空间维度上。在2维的情况下,首先沿一个空间维度进行过滤,之后再沿另一个空间维度进行过滤。对各方程的时间离散采用4阶精度的低存储Runge-Kutta格式[9],该格式包含5个子步。设控制方程可表示为 ∂ϕ/∂t=f( ϕ),取时间步长Δ t=tn+1 -tn,在 t=tn时,有 ϕ=ϕn。Runge-Kutta各子步的计算式依次为:

ϕ1=ϕn+a1f(ϕn)Δt,(13)ϕ2=ϕn+b1f(ϕn)Δt+a2f(ϕ1)Δt,(14)ϕ3=ϕn+b1f(ϕn)Δt+b2f(ϕ1)Δt+a3f(ϕ2)Δt,(15)ϕ4=ϕn+b1f(ϕn)Δt+b2f(ϕ1)Δt+b3f(ϕ2)Δt+a4f(ϕ3)Δt,(16)ϕn+1=ϕn+b1f(ϕn)Δt+b2f(ϕ1)Δt+b3f(ϕ2)Δt+b4f(ϕ3)Δt+b5f(ϕ4)Δt.(17)

其中,计算参数 ai bi( i=1,2,…,5)的取值可见文[9]。

2 模拟对象与网格划分

计算模拟的对象为2维湍流平面射流预混火焰。在入口处有中心射流和伴流。中心射流为甲烷 空气混合物,燃空当量比为0.8, 温度为300 K, 平均速度为36 m/s。射流出口宽度 H=0.8 mm, 基于平均速度的Reynolds数为1 765。伴流为中心射流完全燃烧后的高温产物,速度为9 m/s。

化学反应采用甲烷燃烧的单步反应机理为

CH4+2O2CO2+2H2O.(18)

燃料的瞬时反应速率为[10]

WCH4=Aρ1.5YCH4YO20.5exp-ERT.(19)

其中: A=6.04×107 m1.5/(kg0.5s), E=83 600 J/mol。应用该反应机理进行了1维甲烷层流预混火焰的计算,取燃空当量比为0.8。计算得到的层流火焰传播速度 SL=0.24 m/s, 实验值则为0.259 m/s[11],计算值与实验值相符合。计算的甲烷层流火焰厚度为4.9´10-4 m。火焰厚度的定义为[10]

δL=(T2-T1)/maxTx.(20)

其中: T1为未燃气体混合物的温度, T2为燃烧产物的温度。

x y方向的计算区域尺寸取为 Lx×Ly=1012 H x方向为流向,采用均匀网格划分。横向即 y方向的中心区域划分为均匀网格,两侧采用拉伸网格。由于流场中有燃烧化学反应,为正确求解与分辨反应区结构,在火焰内部需要划分足够多的网格。取 x y方向均匀网格划分区域内的网格尺寸都为40 μm, 约为层流火焰厚度的1/12[10],而湍流火焰厚度一般要大于层流火焰厚度。 y方向边界处的网格尺寸为中心区域网格尺寸的3倍。总的结点数目为201×171。考虑流体受到的重力作用,重力方向与中心射流的流向相反。

3 入口条件的确定

对计算域的各边界均采用特征边界条件[12]。给定中心射流入口处的瞬时速度为 ui= u̅i + u'i。其中: u̅i为射流的平均速度, u'i为扰动速度。射流的流向平均速度及反应标量分布利用双曲正切函数得到[13],横向平均速度取为0, 反应标量分布上不添加扰动。在伴流入口处的各变量分布上也不添加扰动。

中心射流入口扰动速度的生成利用各向同性衰变湍流的能谱。该能谱为[14]

其中: K为湍流动能, kp为能谱峰值对应的波数。依据Rogallo的方法[15],在2维的情况下,得到谱空间中速度矢量 u^的表达式为

其中: e1 e2为谱空间中的两个基矢量; k1与k2分别为 e1 e2方向上的波数,k= ; α(k)为由能谱确定的复数系数,其表达式为

θ为在[0,2π]内均匀分布的随机数。通过式(22)得到的谱空间速度场的能谱仅与给定的能谱形状相同,进一步变换之前需要对该速度场进行重构,以确保在任意的波数处计算的能谱与给定的能谱精确匹配。对重构后的速度场进行Fourier变换得到物理空间的速度场。该速度场仅在谱空间满足散度为0的条件,在物理空间中为使采用前述的有限差分格式离散后的连续方程得到满足,需要进行求解计算。首先求解涡量场,并基于涡量场求解流函数,然后再根据流函数得到速度场。这样获得的速度场满足在给定的离散格式下散度为0的条件。由于速度场作了修正,且有限差分运算在高波数区引入了离散误差,改变了高波数区能谱的形状。为确保在任意波数处计算的能谱与给定的能谱精确匹配,需要在谱空间再次对速度场进行重构,之后进行Fourier变换得到最终的物理空间的速度场。

能谱的形状由峰值波数 kp决定, kp对应的峰值波长为 lp=2π/ kp。随着 kp的增大, lp减小,能量的分布趋于小尺度; 随着 kp的降低, lp增大,能量的分布趋于大尺度。这里引入湍流积分长度尺度 lt=K3 /2。其中, ε=2ν0k2E(k)dk为湍动能耗散率。结合式(21)可以看到, lt K kp的函数,即峰值波数 kp可由湍动能 K和湍流积分长度尺度 lt确定。本文的计算取射流入口处湍流动能 =0.067 u̅c,其中 u̅c为射流入口中心的流向平均速度,并取射流入口处湍流积分长度尺度 lt=0 .297 H[7]。由此确定能谱的峰值波长 lp=0 .5 H

4 计算结果及讨论

DNS计算结果显示, Kolmogorov尺度的最小值出现在射流入口处,为1.543´10-5 m, 而该处网格尺寸是它的2.6倍。随着射流向下游发展, Kolmogorov 尺度逐渐增大,但仍与网格尺寸较为接近。沿射流中心线,湍流积分长度尺度也是逐渐增大的。在中心线上的 x/H=0、 3、 6和9处,湍流积分长度尺度分别为 ltL=0.509 7、 0.705 9、 0.857 3和1.157 2。文[2]在相近条件下的DNS模拟给出的结果则分别为 ltL=0.50、 0.63、 0.82和1.13, 二者的变化趋势是一致的,这说明了本文DNS模拟结果的合理性。

图1a—1c分别给出了 t=8 .02 tu时气体温度、甲烷质量分数和甲烷反应速率的瞬时分布云图,图1d给出了涡量模的瞬时等值线图,这里流动时间尺度 tu=Lx/ u̅c。由图1a的瞬时温度场可以看到,随着射流的发展,剪切层中涡的尺度逐渐增大,并出现了大尺度的拟序结构。从图1b可以看到,甲烷质量分数分布与温度分布有相似的瞬态结构。由图1c的瞬时甲烷反应速率场可以看到,反应速率较大的区域位于剪切层靠近中心射流的内侧。随着射流的发展,瞬时反应速率受到的湍流作用增强,瞬时反应面的折皱变大,面积增多。图1d的瞬时涡量场中给出的是涡量模数值较高的区域,可以看到,近喷口处涡量值较高的湍流旋涡主要发生在剪切层区域。随着射流的发展,剪切层中旋涡的结构增大,涡量值有所降低,湍流旋涡逐渐穿透到射流中心线处。

图1 气体温度、甲烷质量分数、反应速率与涡量模的瞬时分布

图2a—2c分别给出了射流中心线上湍动能、温度脉动均方根值和甲烷质量分数脉动均方根值的分布,展现了从射流入口处即 x/H=0至下游 x/H=9的范围内各物理量的变化。从图2a可以看到,湍动能沿射流中心线逐渐衰减。尽管湍动能呈逐渐衰减的趋势,但它仍保持了一定的数值。随着射流的发展,湍动能沿中心线的衰减逐渐趋于平缓。这说明本文采用的射流入口扰动速度生成方法并未导致扰动快速地衰减。依据能谱生成的入口扰动速度场中包含大小不同尺度的旋涡结构,而能谱的峰值波数是由射流入口处湍动能与湍流积分长度尺度确定的。这就使得入口扰动的湍动能主要分布在尺度较大的旋涡中,这些旋涡的衰减速率相对较慢,因此扰动的衰减也相对较慢。由图2b的温度脉动均方根值分布可以看到,随着射流的发展,温度脉动均方根值沿中心线是逐渐增大的。从图2c可以看到,沿射流中心线,甲烷质量分数脉动均方根值也是逐渐增加的。在射流入口处,仅速度场添加了扰动,而温度与甲烷质量分数的脉动值则均为0。温度与甲烷质量分数脉动均方根值沿中心线逐渐升高,一方面是由于湍流的作用,另一方面则是由于反应区逐渐向中心线处扩展,导致温度脉动和质量分数脉动逐渐增强。而温度和质量分数脉动的增强,又必然会对燃烧化学反应产生影响。

图2 气体湍动能和温度与甲烷质量分数脉动均方根值沿射流中心线的分布

5 结 论

本文利用各向同性衰变湍流的能谱生成了射流入口扰动速度的数据,对甲烷-空气的2维平面湍流预混射流火焰进行了直接数值模拟,得到以下结论:

1) 对气体温度场和甲烷质量分数场的瞬态模拟结果表明,随着射流的发展,剪切层内涡的尺度不断增大,并出现了大尺度的拟序结构。

2) 甲烷反应率分布的瞬态结果表明,剪切层内的瞬时反应速率较高。随着射流的发展和燃烧的进行,反应速率受到的湍流作用增强,导致瞬时反应面出现明显的皱折,反应面积增大。

3) 沿射流中心线,气体湍动能逐渐衰减,但仍然保持了一定的数值,而温度脉动均方根值和甲烷质量分数脉动均方根值则逐渐增长。本文采用的射流入口扰动速度生成方法并未导致出现扰动的快速衰减。

The authors have declared that no competing interests exist.

参考文献
[1] Klein M, Sadiki A, Janicka J. A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations[J]. Journal of Computational Physics, 2003, 186: 652-665. [本文引用:1] [JCR: 2.485]
[2] Sankaran R, Hawkes E R, Chen J H, et al. Structure of a spatially developing turbulent lean methane-air Bunsen flame[J]. Proceedings of the Combustion Institute, 2007, 31: 1291-1298. [本文引用:1] [JCR: 3.828]
[3] Klein M, Sadiki A, Janicka J. Investigation of the influence of the Reynolds number on a plane jet using direct numerical simulation[J]. International Journal of Heat and Fluid Flow, 2003, 24: 785-794. [本文引用:1] [JCR: 1.777]
[4] Lee S, Lele S K, Moin P. Simulation of spatially evolving turbulence and the applicability of Taylor's hypothesis in compressible flow[J]. Physics of Fluids A, 1992, 4: 1521-1530. [本文引用:1]
[5] Le H, Moin P, Kim J. Direct numerical simulation of turbulent flow over a backward-facing step[J]. Journal of Fluid Mechanics, 1997, 330: 349-374. [本文引用:1] [JCR: 2.294]
[6] Stanley S A, Sarkar S, Mellado J P. A study of the flow-field evolution and mixing in a planar turbulent jet using direct numerical simulation[J]. Journal of Fluid Mechanics, 2002, 450: 377-407. [本文引用:1] [JCR: 2.294]
[7] Hawkes E R, Sankaran R, Sutherland J C, et al. Scalar mixing in direct numerical simulations of temporally evolving plane jet flames with skeletal CO/H2kinetics[J]. Proceedings of the Combustion Institute, 2007, 31: 1633-1640. [本文引用:1] [JCR: 3.828]
[8] Kennedy C A, Carpenter M H. Several new numerical methods for compressible shear-layer simulations[J]. Applied Numerical Mathematics, 1994, 14: 397-433. [本文引用:1] [JCR: 1.036]
[9] Kennedy C A, Carpenter M H, Lewis R H. Low-storage, explicit Runge-Kutta scheme for the compressible Navier-Stokes equations[J]. Applied Numerical Mathematics, 2000, 35: 177-219. [本文引用:1] [JCR: 1.036]
[10] Poinsot T, Veynante D. Theoretical and Numerical Combustion[M]. 2 Ed. Philadelphia, PA, USA: Edwards, 2005. [本文引用:3]
[11] Gu X J, Hao M Z, Lawes M, et al. Laminar burning velocity and Markstein lengths of methane-air mixtures[J]. Combustion and Flame, 2000, 121: 41-58. [本文引用:1] [JCR: 3.708]
[12] Sutherland J C, Kennedy C A. Improved boundary conditons for viscous, reacting compressible flows[J]. Journal of Computational Physics, 2003, 191: 502-524. [本文引用:1] [JCR: 2.485]
[13] Sutherland J C. Evaluation of Mixing and Reacting Models for Large-Eddy Simulation of Nonpremixed Combustion Using Direct Numerical Simulation [D]. Salt Lake City, UT, USA: University of Utah, 2004. [本文引用:1]
[14] Haworth D C, Poinsot T J. Numerical simulations of Lewis number effects in turbulent premixed flames[J]. Journal of Fluid Mechanics, 1992, 244: 405-436. [本文引用:1] [JCR: 2.294]
[15] Rogallo R S. Numerical Experiments in Homogeneous Turbulence, NASA-TM-81315 [R]. Moffett Field, CA, USA: NASA Ames Research Center, 1981. [本文引用:1]