不同辐射下反应羽流三维直接数值模拟

袁步平,邱 榕,蒋 勇

(中国科学技术大学火灾科学国家重点实验室,安徽合肥,230026)

不同辐射下反应羽流三维直接数值模拟

袁步平,邱 榕＊,蒋 勇

(中国科学技术大学火灾科学国家重点实验室,安徽合肥,230026)

运用直接数值模拟方法对低速流动下的浮力反应羽流进行了三维数值模拟,分析了不同辐射通量下无量纲流场中近场浮力羽流的结构特性。在有重力场和燃烧释热的情况下,反应羽流自然演化生成周期性的大涡结构,随着辐射强度的增加,不仅对温度场有显著影响,同时还影响着羽流涡旋的发展及其下游三维涡旋结构之间的相互作用。同时,辐射并不破坏流场的对称性,有无辐射的流场涡旋结构均关于几何中心对称。

直接数值模拟;反应羽流;热辐射;涡旋结构

0 引言

火灾科学基础研究的任务是探索和认识火灾现象及其过程的机理和规律,而对浮力反应羽流的研究是理解火灾现象非常重要的环节。实验和理论分析很难全面研究浮力反应羽流近场的动力学详细信息。最近发展起来的大涡模拟方法在分析反应羽流的精细结构时并不适合,传统的雷诺平均方法也不能很好反映真实流场的涡结构和不稳定性。而不依赖于任何模型参数的直接数值模拟可以提供所有时间和空间尺度的详细流场信息,通过直接求解Navier-Stokes方程,给出真实的随时空变化的流场信息。

由于涉及大范围的空间和时间尺度,用直接数值模拟方法研究流动现象本身就是一个巨大的挑战,而对耦合有化学反应和复杂流体运动的浮力反应羽流难度更大。近十几年,随着计算机的迅速发展,极大的促进了湍流燃烧数值模拟的发展。Mell和 Mcgrattan[1]、Poinsot和 Candel[2]、Vervisch 和Poinsot[3]等研究者就开始利用DNS能精确求解各个尺度运动的优势研究了受迫羽流近场动力学和化学释热的影响,但是这些研究都局限于在二维平面上。三维空间发展的浮力羽流研究受到限制,主要是因为三维直接数值模拟极大的数值难度以及相关的边界条件。Luo[4]、Bedat[5]等人在1999年对直接数值模拟方法计算耦合有流场动力学和燃烧的物理问题提供了可能,Wilson和Demuren[6]同时用数值模拟方法研究了不考虑浮力影响的非圆形射流。不久前,Jiang和Luo采用直接数值模拟技术较为系统地研究了二维和三维的浮力反应羽流[7-10],揭示了热羽流和浮力反应羽流大涡结构输运和演变的机理和规律。但这些国外学者都是对羽流的不稳定性和涡旋演变规律进行研究,由于数值计算上的难度,均没有考虑到热辐射的问题,而热辐射在火灾的热量传递中起着不可忽视的作用。因此,本文用直接数值模拟技术来研究不同辐射情况下三维反应羽流的特征变化。所研究的具体物理问题是用圆形喷口将燃料垂直射入空气环境中进行燃烧化学反应产生低速流动的非预混火焰。

1 数学模型

本文以美国国家标准技术研究所开发的模拟软件FDS为研究基础,利用其源代码的开放性,应用其直接数值模拟(DNS)方法来进行研究。与大涡模拟方法相比,直接数值模拟不需要对流动建立模型,而是直接求解流动控制方程—Navier-Stokes方程,因此能够获得所有尺度上的流动信息,准确地预测各种物理量的时空变化。Beji[11]等人曾用 FDS软件的DNS方法研究了二维的甲烷燃烧反应,并与Zhou和 Gore[12]的甲烷浮力羽流实验数据进行对比,模拟结果与实验数据得到了很好的吻合,说明了用DNS方法研究浮力反应羽流的模拟结果具有较高的真实性和可靠性。

本文研究对象为三维、非定常、考虑浮力的粘性反应流,化学反应采用有限速率的Arrhennius动力学方程求解。在数值方法方面,国内外学者在进行相关的直接数值模拟时,一般采用六阶空间离散高精度格式,因此必须采用非反射的边界条件去克服伪物理波反射。这种边界条件计算量及其巨大,在实际中难以推广。蒋勇[13]等在空间离散高精度与低精格式的对比中发现,对羽流的直接数值模拟而言,时间离散的精度对计算结果可靠性的影响极大,而二阶空间离散精度是足够的。本文DNS采用隐式迭代的方法,其时间精度可以达到三阶,因此空间精度采用二阶精度就足够了,而不需要采用伪物理边界。

本文对圆形喷口轴对称浮力反应羽流进行三维直接数值模拟,分别从流动参数和燃烧参量方面分析羽流的结构特性,讨论了自由羽流在三种不同辐射分数下(0%,10%,35%)对反应羽流结构的影响。甲烷燃料从底部圆形出口喷出,垂直流入静止的空气中,水平速度为零,在来流边界上没有扰动。燃料的质量分数 Ya=1,环境温度为 T0,喷口处边界采用无滑移的边界条件,四周及顶部为自然开口。为了研究的通用性,利用喷口直径D、环境温度 T0、出口速度U0、环境密度ρ∞等参数对流场进行无量纲化。在笛卡尔坐标系下,喷口平面为 X-Y面,流动方向为 Z轴方向,g是竖直向下的重力加速度。为了使边界对模拟结果的影响将至最低,以及考虑计算量的大小,通过对比分析确定无量纲后的计算区域尺寸为Lx×Ly×Lz=4×4×8。在模拟过程中,考虑到三维直接数值模拟计算量巨大,为了减少网格尺度划分难度(DNS要求最小网格尺度足以分辨最小涡的运动,即与 Kolmogorov长度尺度[14]在一个量级),故只模拟低雷诺数和低弗鲁德数的流场。根据入口参考量确定雷诺数 Re=500,Froude数Fr=0.375。计算模型如图1所示。

图1 计算模型图Fig.1 The simulation diagram

为了能够保证计算精度,又能尽可能节省计算时间,需要进行网格独立性分析。本文以比例约为1/4让网格数目逐步增加,分别对四种网格进行计算分析,结果如表1所示(对比量为中心线 Z=2位置处轴向速率值)。

表1 不同网格密度下的时均速度及其与中心线 Z=2位置处轴向速率的比值Table 1 Averaged velocity and ratio of different mesh grid

从计算结果看,比100×100×200网格点更密 的网格划分并不能明显改变计算结果,误差在1%以下,因此本文计算时采用的网格数为100×100×200,三个方向的网格密度相等且等于0.04D。为了保证数值稳定,时间步长由 Courant-Friedrichs-Lewy(CFL)条件限制。

2 计算结果及讨论

浮力羽流表现出来的闪烁或膨胀等周期性脉动现象是由于近场浮力引起大涡结构的形成和输运造成的[15-17]。图2给出了三个算例在 T=20时 Y=Ly/2平面内的瞬时等温线。从图中可以看出,在重力和反应放热的影响下流场产生浮力,由于绝对浮力的不稳定性,在计算区域内自然演化生成大涡结构。大涡结构的流向输运导致流场产生周期性的脉动变化。从图2中a,b,c三种不同辐射分数对比可以看出,辐射不仅使得羽流的最高温度有明显下降(辐射35%的最高等温线下降1.7T0),而且对涡旋的发展也有一定的影响。这是因为辐射导致燃烧产生的热量有部分损失于周围介质中,从而直接导致温度极值的大幅下降,另温度差减少使得羽流与环境的密度差变少,浮力卷吸作用减弱,影响涡旋发展。从无辐射的a图和辐射35%的c图可以明显的看出,在流动发展的下游,无辐射的流场涡旋结构更小更细。这是由于无辐射的流场三维涡旋相互作用更加强烈,在下游涡旋破碎所致。从图中还可以看出,由于没有外界的扰动影响,辐射并没有破坏流场的对称性,涡旋结构关于几何中心对称。

图2 三种辐射情况下浮力反应羽流在 T=20时Y=Ly/2平面的瞬时等温线(极值之间共10条等值线)Fig.2 Temperature contours onY=Ly/2 plane of buoyant reactive plume atT=20 for different radiation(10 contours between the minimum and maximum as indicated)

图3和图4分别给出了三个算例在 Z=3和 Z=4截面内的瞬时等温线。在两图中可以看出由于旋涡变形形成的复杂三维结构非常明显,越往下游相互作用越强烈,结构越复杂,其切向的羽流结构不再保持原有的圆形形状,是由于卷吸混合及涡旋结构发展变形所致。从图3,图4不同的辐射分数对比图可以看出,辐射可以改变流场的涡旋结构,由于辐射量的增加,导致流场温度明显降低,温度的下降使浮力和旋涡扩张作用受到影响。从两图的a,b,c对比还可以看到,随着辐射量的增加,火焰结构有所收缩,膨胀性变小。除此之外,由于辐射的作用,涡旋结构也变得简单,涡旋破碎程度减少。可能是由于辐射的增加导致温度降低,流场密度差异性变小,造成卷吸和涡旋的相互作用变弱。由于没有受到不对称的外加扰动,因此,羽流结构还是以几何中心对称。从图2,3,4还可以细微的发现,辐射的增加会使涡旋发展变缓慢,其大涡出现在更接近下游的位置。

图5为三种不同辐射分数下中心截面径向涡量ωy的瞬时等值线。与图2相对应,可以看出随着辐射量的增加,涡量极值也有明显的降低,说明辐射使得流场的漩涡强度减少。除此之外,同样可以得出在有35%辐射的情况下,涡量产生较为滞后,大涡结构迟于无辐射的流场。这可能是因为随着辐射的增加,燃烧产生的热量通过辐射而损失,羽流温度相应降低,羽流与环境的密度差异变少,浮力减少,从而羽流与空气的混合卷吸等相互作用变弱,导致涡旋强度减少,大涡结构发展变慢。

图3 三种辐射情况下浮力反应羽流在 T=20时 Z=3平面的瞬时等温线(极值之间共10条等值线)Fig.3 Temperature contours onZ=3 plane of buoyant reactive plume atT=20 for different radiation(10 contours between the minimum and maximum as indicated)

图4 三种辐射情况下浮力反应羽流在 T=20时 Z=4平面的瞬时等温线(极值之间共10条等值线)Fig.4 Temperature contours onZ=4 plane of buoyant reactive plume atT=20 for different radiation(10 contours between the minimum and maximum as indicated)

从图6中水平截面Z=3的瞬时涡量图可以看出,辐射不仅使得涡量值变小,还使得羽流对周围的影响范围变小,辐射的增加抑制了反应羽流涡旋的发展和膨胀。

图7分别给出了三个算例在不同高度处中心线流向速率随时间的变化曲线。由(a)可以看到,三个算例在 Z=1处的中心线流向速率的变化均比较规则,这与喷口附近浮力驱动的大涡结构的形成与对流有很大关系,这就是所谓的浮力羽流的膨胀或闪烁现象。差异仅仅是由于辐射的作用使得密度差变小,从而浮力引起的流向加速减少。但是,在靠近下游的位置 Z=3和 Z=5位置处(见图(b),(c)),中心线流向速率不是很连贯,有较大的波动。这可能是随着流向涡量的发展,三维旋涡相互作用导致涡旋破碎的结果。但是,在含有35%辐射的情况下,这种不连贯的波动却不明显,这正是由于随着辐射的增加,阻碍了涡量的发展,漩涡之间的相互作用变弱,其中心线流速的变化仍较为规则。

图8给出不同辐射下不同高度处平均流向速度的径向分布。在 Z=1高度处,三种情况下的速度还基本保持帽型分布,但是流向速度由于浮力的影响有明显的加速。在下游 Z=3处,羽流轴向速度向四周扩张,这是由于燃烧引起的三维涡旋结构相互作用加强混合的结果。对比有无辐射的情况可以发现,在整个流域,辐射的增加不仅使浮力引起的流向平均速率降低,而且减少了三维涡旋的相互作用,收缩了羽流涡旋的影响范围。

图5 三种辐射情况下浮力反应羽流在 T=20时Y=Ly/2平面的涡量图(极值之间共20条等值线)Fig.5 Vorticity contours onY=Ly/2 plane of buoyant reactive plume atT=20 for different radiation(20 contours between the minimum and maximum as indicated)

图6 两种辐射情况下浮力反应羽流在 T=20时 Z=3平面的涡量图(极值之间共20条等值线)Fig.6 Vorticity contours onZ=3 plane of buoyant reactive plume atT=20 for different radiation(20 contours between the minimum and maximum as indicated)

3 结论

本文对浮力反应羽流进行了三维直接数值模拟,用DNS的方法分析不同辐射情况对羽流结构的影响。计算结果表明:

(1)辐射不仅使反应羽流的温度有明显下降,对涡旋的发展也有一定的影响。其使流场火焰结构有所收缩,膨胀性变小,影响范围变窄。无辐射的流场三维涡旋相互作用更加强烈,在下游涡旋破碎更小更细。辐射的增加使涡量极值也有明显的降低,流场的漩涡强度减少,涡量产生较为滞后,迟于无辐射的流场;

(2)从中心线流速的变化可以看出,喷口附近流场周期性明显,变化均较为规则,辐射仅仅减缓其流向速度的增加。但在靠近下游的位置,无辐射的流场中心线流率不再连贯,有较大的波动,而辐射35%的流场这种波动却不明显,其能阻碍下游漩涡之间的相互作用,混合卷吸强度变弱;

(3)辐射并不能破坏流场的对称性,有无辐射的流场涡旋结构均关于几何中心对称。

图7 不同辐射下浮力反应羽流中心线流速随时间的变化Fig.7 Comparison of the instantaneous centerline streamwise velocity variations atZ=1,Z=3and Z=5 in different radiation

图8 不同辐射下浮力反应羽流中心线流速径向时均值Fig.8 Time-averaged sreamwise velocity profiles on the central plane of buoyant reactive plume at different radiation

[1]Mell WE,Mcgrattan KB and Baum HR.Numerical simulation of combustion in fire plumes[A],In:Proceedings of the Twenty-Sixth Symposium(International)on Combustion[C].The Combustion Institute,Pittsburgh.1996:1523-1530.

[2]Poinsot T,Candel S and Trouve A.Applications of direct numerical simulation to premixed turbulent combustion[J].Prog.Energy Combust.1996,21:531-576.

[3]Vervisch L and Poinsot T.Direct numerical simulation of non-premixed turbulent flames[J].Annual Rev.Fluid Mech,1998,30:655-691.

[4]Luo KH.Combustion effects on turbulence in a partially premixed supersonic diffusion flame[J].Combust and Flame,1999,119:417-435.

[5]Bedat B,Egolfopoulos FN and Poinsot T.Direct numerical simulation of heat release and NOxformation in turbulent nonpremixed flames[J].Combust and Flame,1999,119:69-83.

[6]Wilson RV,Demuren AO.Numerical simulation of turbulent jets with rectangular cross-section[J].ASME J.Fluids Eng.1998,120:285-290.

[7]Jang X,Luo KH.Spatial direct numerical simulation of the large vortical structures in forced plumes[J].Flow,Turbulence and Combustion,2000,64:43-69.

[8]Jang X,Luo KH.Combustion-induced buoyancy effects of an axisymmetric reactive plume[J].Proceeding of the combustion institute,2000,28:1986-1995.

[9]Jang X,Luo KH.Direct numerical simulation of the near field dynamics of a rectangular reactive plume[J].Heat and Fluid Flow,2001,22:633-642.

[10]Jang X,Luo KH.Direct numerical simulation of transitional noncircular buoyant reactive jets[J].Theoretical and ComputationalFluid Dynamics,2001,15(3):183-198.

[11]Beji T,Zhang J P,Yao W,Delichatsios M.On the limitations of constant Prandtl and Schmidt numbers assumption in LES simulations of reacting buoyant plumes[A],in:Fourth European Combustion Meeting[C].Vienna University of Technology, Vienna,Austria,2009.

[12]Zhou XC,Gore J P.Experimental Estimation of Thermal Expansion and Vorticity Distribution in a Buoyant Diffusion Flame[J].Proc.Combust Inst,1998,27:2767-2773.

[13]蒋勇,邱榕,宋崇林.基于多重条件矩封闭模型的反应羽流大涡模拟[J].燃烧科学与技术,2009,15(1):16-21.

[14]Frisch U.Turbulence:The legacy of A.N.Kolmogorov[J].Cambridge University Press,1995.

[15]Cetegen BM,Ahmed TA.Experiments on the periodic instability of buoyant plumes and pool fires[J].Combust.Flame,1993,93:157-184.

[16]Lingens A,Reeker M,Schreiber M.Instability of buoyant diffusion flames[J].Exp.Fluids,1996,20:241-248.

[17]Maxworthy T.The flickering candle:transition to a global oscillation in a thermal plume[J].J.Fluid Mech,1999,390:297-323.

Spatial direct numerical simulation of reactive plumes under different radiations

YUAN Bu-ping,QIU Rong,J IANG Yong

(State Key Laboratory of Fire Science,University of Science and Technology of China,Hefei 230026,China)

Spatial direct numerical simulation(DNS)is used to study the near field dynamics of low-speed buoyant reactive plumes.The structure characteristics of non-dimensional flow field of reactive plumes at different radiation magnitudes are investigated.Reactive plumes display periodic large vortical structures naturally in the flow field,due to the heat release and gravitational effects.As the radiation fraction is increased,the temperature of the flow field,the vorticity of the plume and the interactions of its 3D vortices at the downstream of the flow filed will all be significantly affected.It is also shown that radiation does not destroy the symmetry of the flow filed,and all the vortical structures are symmetrical as to the geometric center.

Direct numerical simulation;reactive plumes;radiation;vortical structure.

TK16

A

1004-5309(2011)-0087-06

2011-02-27;修改日期:2011-03-26

国家自然科学基金资助(No.50876097)

袁步平(1985-),男,汉族,中国科学技术大学火灾科学国家重点实验室硕士研究生,主要从事火灾动力学数值模拟研究。

邱榕,副教授,E-mail:rqh@ustc.eu.cn.