（工程热物理专业论文）基于无结构网格的有限体积颗粒相结合的方法求解湍流联合的pdf方程.pdf

摘要 本文采用一种新的混合算法。对端流的脉动速度一频率 一标量联合p d f 输誊撞型左程进行数值求解采用一种基于 三角形无结构网格的有限体积法，求解质量、动量、能量平均 守恒方程，压力通过理想气体的状态方程求解，所求得的密度 场、速度场、压力场用于p d f 方程的计算求解；脉动速度一频 率一标量联合p d f 的输运模型方程用拉格朗日m o n t e - c a r l o ( 颗 粒) 方法求解，并为平均守恒方程提供雷诺应力、热通量、能 量源项和7 有限体积法和颗粒方法之间的数据交换用一基 于伽辽金算法的线性插值算法来实现本文对该算法的相容性 进行了深入的讨论，提出相应的修正算法和时间平均方法，以 减小颗粒方法的结果中的统计误差，以确保获得在数值解水平 上的相容利用该算法，对二维v 型槽后的湍流冷态流场和湍 流反应流场进行了数值模拟，最后对计算结果进行了分析和 讨论 关键词有限体积法，m o n t e - c a r l o 方法，几率密度函数，无 结构网格，湍流反应流，相容性 a b s t r a c t i nt h i sp a p e r ，an e w h y b r i da l g o r i t h mf o rj o i n tf l u c t u a t i n g - v e l o c i t y - f r e q u e n c y - c o m p o s i t i o np d fe q u a t i o n so ft u r b u l e n t ( r e a c t i v e ) f l o wi sp r e s e n t e d t h e 血e t h o d c o m b i e so faf i n i t e v o l u m em e t h o da n da m o n t e - c a r l o ( p a r t i c l e lm e t h o d t h e f i n i t e - v o l u m em e t h o di sc a r r i e do u to nt r i a n g l eu n s t r u c t u r e dm e s h a n di 8u s e dt os o l v e t h er e y n o l d sa v e r a g e de u l e re q u a t i o n s t h ep r e s s u r eo ft h ef i e l di sc a l c u l a g e db ya s t a t ee q u a t i o no fi d e a lg a s i tp r o v i d e st h em e a nd e n s i t y ，m e a nv e l o c i t ya n dp r e s s u r e f o rp d fe q u a t i o n s al a g r a n g i a nm o n t e - c a r l om e t h o di su s e dt os o l y et h ej o i n t f l u c t u a t i n g - v e l o c i t y - f r e q u e n c y - c o m p o s i t i o np d fe q u a t i o n s ，a n dp r o v i d e sr e y n o l d s s t r e s s ，h e a tf l u x ，甲( a n de n e r g ys o u r c et e r mf o rr e a c t i v ef l o w ) f o rt h ef i n i t e v o l u m e s c h e m e ag a l e r k i n b a s e dl i n e a r i n t e r p o l a t i o nt e c h n i q u e i se m p l o y e dt oe x c h a n g et h e d a t ab e t w e e nt h ef i n i t e - v o l u m em e t h o da n dm o n t e c a r l om e t h o d t h ec o n s i s t e n to f t h eh y b r i ds c h e m ei sd i s c u s s e d s o m ec o r r e c ta l g o r i t h ma n dt i m e - a v e r a g et e c h n i q u e i si n t r o d u c e dt or e d u c et h es t a t i s t i c a le r r o r n u m e r i c a le x p e r i m e n t sa l ec a r r i e do u t o i lt h et u r b u l e n tn o n r e a c t i n gf l o w sa sw e l la sp r e m i x e dt u r b u l e n tr e a c t i v ef l o w sp a s t t h eb i n i f - b o d vi nar e c t i l i n e a rc h a n n e l t h er e s u l t sa l ec o m p a r e dw i t he x p e r i m e n t a l m e a s u r e m e n t s l e sa n dk em o d e l i ts h o w st h a tt h er e s u l to ft h i ss c h e m ei s i n g o o da g r e e m e n tw i t he x p e r i m e n t a ld a t a k e y w o r d s f i n i t ev o l u m em e t h o d ，m o n t e - c a r l om e t h o d ，p d ff u n c t i o n ，u n s t r u c t u r e dm e s h ，t u r b u l e n tr e a c t i v ef l o w s ，c o n s i s t e n t i i 第1 章引言 1 1 湍流燃烧数值模拟的方法 湍流广泛存在于自然界中，工程实际中涉及的燃烧过程绝大多数是湍流 燃烧，燃烧是当前获取电力的重要手段，是交通运输的主要运行方式，但燃 烧过程同时产生大量的污染物，是全球温室气体的主要来源因此，湍流燃 烧的研究具有重要的学术价值和应用背景湍流燃烧的过程，涉及到流体力 学、传热传质学、化学反应动力学等诸多个学科，这些分过程相互耦合，相互 作用对该领域的研究，能有效地促进这些学科和交叉学科的发展，并为工 程应用提供积极的理论指导 随计算机和数值方法的发展，直接数值模拟( d i r e c tn u m e r i c a ls i m u l a t i o n ) 的 方法已经取得一定的进展，但离实际工程应用还很远如果用直接数值模拟 的方法求解三维非定常n a v i e r s t o k e s 方程，目前的计算机只能用于简单的低雷 诺数流场的计算大涡模拟( l a r g e - e d d ys i m u l a t i o n ) 方法采用较大的网格，用亚 网格模型来模拟小尺度湍流的运动，与直接数值模拟相比较，对计算机的要 求较低，但计算量仍然很大，现在也还只能用于一些简单的流动 目前在工程实践中广泛使用的是基于r e y n o l d s 平均的湍流宏观模拟方 法，其基本点是利用数学模型，将雷诺时均方程或湍流特征量输运方程中高 阶的未知关联项用低阶关联项或时均量来表示，从而使雷诺时均方程封闭 常用模型有混合长度模型、双方程模型、雷诺应力输运方程模型等用湍流 统计矩模型处理湍流反应流问题时，我们需要对平均化学反应速率项进行模 拟，现在常用的有涡旋破碎模型( e b u ) 、 一e 一9 等模型由于这些模型都是 建立在简单化学反应系统的假定基础上，无法全面考虑分子输运和化学动力 学等因素的作用，不能解决复杂的化学反应问题 几率密度函数( p r o b a b i l i t yd e n s i t yf u n c t i o n 输运方程的方法中，各个单点 相关都可以从联合几率密度函数求得，获得封闭所以p d f 输运方程的形式 与宏观模拟方程截然不同在变密度流动中，p d f 方程中与湍流输运和化学 反应速率有关的项都是以封闭的形式出现的，无须模拟 4 4 ，4 6 】因此p d f 方 法在湍流反应流的数值计算中有着独到的优势理论上，任意复杂的化学反 应机理都可以用p d f 方法精确计算同时，p d f 输运方程所能提供的信息也 比普通的湍流模型要多因此，p d f 输运方程很适合模拟那些必须考虑湍流 流动、复杂化学反应机理及其相合耦台的湍流反应流问题，如点火、熄火、 污染物的生成等当然，p d f 方程不是完全封闭的，其中粘性应力张量，脉 动压力梯度项和小尺度分子扩散的条件平均值尚需模拟在求解p d f 方程的 2 第1 章引言 时候。平均压力梯度无法从p d f 方程中得到，需要用其它的方法给出为了 获得压力项，我们需要求解p o i s s o n 方程，而p o i s s o n 方程比较难收敛，因此影 响了求解效率研究还表明，平均场的统计脉动，对计算结果的影响很大 提供给p d f 方程中模型的平均场数值是否准确，对整个p d f 方程的求解结果 和求解效率都有很大的影响因此，p d f 方法需要有一套有效的数值计算方 法 5 6 】，才能更好地应用于工程计算本文工作的目的就是试图发展一种新的 算法，把求解脉动速度一标量一频率联合p d f 方程和求解平均速度、平均能 量的统计矩方程结合起来，以降低计算误差和提高计算效率 1 2 p d f 方程求解方法的发展 几率密度函数方法始于1 9 6 9 年l u n d 昏 e n 3 3 1 的工作，他推导出速度的联合 p d f 的方程1 9 7 6 年d o p a z o 1 2 】、p o p e 4 2 】各自独立推导出化学热力学惨伤的 联合p d f 方程在以后p d f 、方法的研究中，许多学者致力于方程的模化和数 值模拟，而对p d f 方程求解算法的研究则相对比较少 p d f 输运方程是一个多维的标量输运方程，很不适合用有限体积法或者 有限差分法求解m o n t e - c a r l o ( 颗粒) 方法在求解多维问题时，计算量只随着维 敦的提高而线性增加，因此，m o n t e - c a r l o 方法在p d f 方程的研究中被广琵地 使用m o n t e - c a r l o 方法中，首先根据p d f 方程建立组描述颗粒变化的随机 微分方程，继而考虑大量的具有速度、化学组分以及温度等性质的随机颗粒 随着时间演化，最后根据颗粒达到统计定常时的性质，获得p d f 方程的解 m o n t e - c a r l o 方法分为欧拉( e u l e r i a n ) 方法和拉格朗日( l a g m l g i a n ) 方法欧 拉方法中，颗粒定位在物理空间的网格节点上田 ，而拉格朗丑方法中的颗粒 则是在物理空间中连续分布的两者相比较，后者更易建立模型，所获得的 结果也更加直观所以近十多年来人们把目光主要投向了拉格朗日方法， 拉格朗日方法中，通常有一个或者多个平均量是用有限体积或有限差分：- 方法求解偏微分方程的解给出的，其它平均量由颗粒方法给出( 为了区别于 从偏微分方程求得的平均量这类平均量称为“颗粒场量”) p d f 输运方程 不同的求解算法之间的重要区别，就在于利用什么样的平均场及在颗粒方程 中如何使用( 见表1 1 ) p d f 输运方程可以用单纯的m o n t e - c a r l o 方法求解w e l t o n 和p o p e 5 3 将 s p h ( s m o o t h e dp a r t i c l eh y d r o d y n a m i c s ) 的思想，应用于p d f 方程的求解该方 法就是完全用颗粒方法求解速度标量频率联合p d f 输运方程，不徭要从 阿格处获得统计平均量速度一频率一标量联合p d f 方程构成一个封闭的湍 流模型方程，除了流体的一些热物性和初始、边界条件之外，无需额外的信 息，这为建立一套自相容的颗粒方法提供了可能除了平均压强，颗粒方程 中所利用的平均参数都是颗粒场量【5 6 由于颗粒场量本身就包含着较大的 1 , 2p d f 方程求解方法的发展3 统计误差，颗粒场量在颗粒方程中的使用可能会引起更大的计算偏差为了 使偏差和统计误差降低到足够小，要求每个单元里有大量的颗粒，大大增加 了计算量在计算平均压力时采用了一种复杂的算法求解p o i s s o n 方程，计算 精度和效率都不高d e l a r u e 和p o p e 1 1 采用p d f 方法求解可压缩湍流时，平 均压力场改为由颗粒场的属性直接获取，从而使p d f 求解器不荐依赖于类似 于有限体积法的求解器 p o p e 在文献【4 3 中设计了一种有限差分与m o n t e - c a r l o 相结合的方法，用 于求解湍流反应流其主要思路是：用有限差分方法将化学热力学参数的联 合p d f 方程在物理空间和时间坐标上离散化，然后将每个节点上的p d f 转化 为反映这个分布的、由n 个随机样本或随机颗粒组成的一个系综用节点之 间随机颗粒的交换模拟物理空间上的对流和扩散过程，用同一节点上不同颗 粒的混合模拟湍流的小尺度混合，最后计算各颗粒由于化学反应所引起的参 数的变化 h a w o r t h 2 1 等人采用有限体积法和m o n t e - c a r l o 相结合的方法对湍流流场 进行了求解求解速度一频率联合p d f 输运模型方程，用基于s i m p l e 算法 的有限体积法求解雷诺方程，求解的量包括平均动量、平均压强、湍流动能、 湍流动能耗散率和平均内能，平均密度根据状态方程获得，为p d f 方程提供 平均速度、压力梯度，湍流耗散率( 湍流时间尺度) 等；用m o n t e r - c a r l o 方法求 解p d f 方程，为雷诺方程提供雷诺应力但该算法并没有解决两组方程的相 容性问题，其中的速度场被重复计算了另一方面，在平均场方程中由于颗 粒密度的脉动会引起一系列收敛困难的问题 陈义良【5 7 等人吸收了 一模型简单和p d f 方法能精确计算详细化学反 应机理的优点，发展了一个采用差分和m o n t e - c a r l o 求解相结合的方法具体 做法是采用差分方法计算速度场，用m o n t e - c a r l o 方法求解化学热力学参数， 应用于湍流燃烧的工程计算 在c o n e a 和p o p e 8 】开发的p d f 2 d s 程序中，速度和化学热力学参数的 p d f 方程用m o n t e - c a r l o 方法求解，m o n t ec a r l o 方法所需的平均速度，压强、 湍流动能和耗散率由有限体积法提供，并将平均密度场返回给有限体积法， 用以求解 一e 方程；压强是利用连续方程和动量方程获得的方程来求解的， 为了与平均压力场保持一致，在每个时间步后，对平均速度场进行了修正 其中平均速度和雷诺应力被重复计算由于颗粒方法所给出的结果中包含着 较大的统计误差，在平均场方程中使用颗粒方法所获得的密度，会引起一系 列收敛性的问题 。 h s u 等f 2 2 1 用m o n t e - c a r l o 方法求解组分和焓的联合p d f 方程，用有限体 积法求解超声速湍流反应流有限体积法为m o n t e - c a r l o 方法提供速度场和压 力场。并利用m o n t e c a r l o 方法所获得的组分的质量分数，求解平均场的温度 4 第1 章引言 和压强 j j a b e r i 2 3 ，7 】等人发展了大涡模拟和p d f 方法相结合的求解方法，用、 f d f ( f i l t e r e dd e a s 时f 慨t i o n ) ，f m d f ( f i l t e r e d m a s s d e n s i t yf u n c t i o n ) 方法推导出 大涡模拟的方程，大涡模拟为p d f 方程提供平均速度场，p d f 方程为l e s 提 供较为精确的s g s 量与前面各种方法不一样的是，有限体积法求解的是连 续方程、动量方程和能量方程，因为颗粒密度场的统计误差所造成的计算不 稳定现象得到很大改善 上面所提及的方法都存在相容性方面的问题，m o r a d o g l u 、p e n n y 和 p o p e 2 5 ，如】建立了一套在理论上完全相容的混合有限体积法和颗粒方法的算 法，并编制了p d f 2 d f 矿和h y b 2 d 程序具体所采用的算法是的混合算 法该算法中，某些量( 如平均密度、平均内能) 被重复求解了“理论上完 全相容的意思是，如果雷诺平均方程和联合p d f 输运模型方程精确求解的 话，被重复求解的量是相同的并对射流扩散火焰【2 5 】、钝体后回流【2 6 】进行 了计算 。 1 3 本文的工作 本文利用h y b 2 d 的思想，在无结构三角形网格上，用显式的有限体积法 求解二维空间里的雷诺平均e u l e r 方程用m o n t e - c a r l o 方法求解脉动速度一标 量一频率的联合p d f 输运方程；压强根据理想气体的状态方程来确定；有限 体积法为颗粒方法提供平均密度、平均速度和平均压强，颗粒方法则为有限 体积法提供雷诺应力、能量源项、热通量和比热比7 等信息两者之间的数 据交换采用一种基于伽辽金( g a r l e l d a ) 的插值算法来实现 本文的主要内容有如下几个方面； l 。首先介绍本文中所涉及的热力学状态量的求解； 2 求解的方程僵。诺平均e u l e r 方程、状态方程、压力方程、脉动速度标 量频率联合p d f 方程等) ，以及求解过程中所采用的湍流模型和燃烧模型； 3 数值计算方法和f v m c 2 d 程序的计算流程； 4 对整个算法的相容性进行了分析，并介绍若干修正算法( 包括程序中 所使用的和其它等价的修正算法) ； 。 。 5 对二维v 型槽后的湍流场和燃烧场进行了计算，并对计算结果进行分 析和讨论 、 7 1 3 本文的工作 5 颗粒方程平均方程 颗粒方颗粒方平均方重复 求解颗粒程中的程中的程中的计算 方法属性颗粒量平均量平均量颗粒量的量 标量联合m ，y y ，fu ，u ，( p )卢 p d f 【4 4 】 七 七，e 标量联合m 4 y ，= ( y + ) y u ，a6 - ，女，e ，i 。 z p d f 【2 3 】 饼，厉j( t ：，啪 速度一标m ，y ，u y ，fu ，御，u ，fu ，“q ) 量联合p d f 矗e怠、e ( p d f 2 d s ) 8 】 速度一标m ，y ，u y ，ku e ，u ，e ，厉 u 量联合妇) u ：u j l x ) p d f l 2 i 】 速度频m ，y ，u ，u ，寸，f如)( p u 率一标量 u 。0 n 联合p d f 似u ；i x ) ( p d f 2 d v ) 5 6 】 速度频m ，y ，u ，u ，y ， - 率标量u ，p o ，( p ) ， 联合p d f ( u ：q 盼 ( s p h ) 5 3 速度频m ，y ，u ，y ，au 劬u ) ，( u ：“；l x ) ， 一= 卢e 率一标量u ， 0 ，n ，己，卢( u ：f ：i x ) ， 联合p d f扣：u ：i x ) e ( p d f 2 d - f v ) c 2 5 ，如 表1 1 ：p d f 方程的不同求解算法的比较 第2 章湍流反应流的基本方程及模型 本章介绍计算湍流反应流所涉及的一些基本方程，以及所采用的具体模 型假定流体是理想气体 2 1 热力学状态方程 求解时需要给定热力学状态方程( 混合物的焓和温度之间的关系) 混合 气体的化学热力学状态参数包括密度p 、温度t 、压强p 、各组分的质量分数 y = m ，硷，y i ) r 和焓等设组分a 的分子量为，则该组分的气体常数 为： 、 孵 = 暑 ( 2 1 ) rr n 其中瞬= 8 3 1 4 5 1 j m o l - k 为通用气体常数 混合气体、各种组分的的焓值和比热由c h e m k i n 程序包提供 2 8 】该程序 包提供的热力学参数数据库t h e r m d a t 给出了每一种物质在参考温度t = 2 9 8 1 5 k 时的比焓破，以及以t 的多项式函数形式给出的定压比热口) 温度为t 时气体的比焓为： 。( t ) = 印。( t ) d t + = 弓。( t f ( 2 2 ) 其中平均定压比热由下式给出： ( t ) = 亍1 五t 铀( t 7 ) d r 比内能为： e s n ( t ) = ( b o ( r ) 一) t = 日时( ? ) t 其中平均定容比热的定义跟乙。类似： k ( t ) = 亍1 五t c e l 。( t ) d r 混合气体的比显焓： ( 2 3 ) j 。 ( 2 4 ) ( 2 5 ) k ( y ，t ) = 圪h 。( t ) ( 2 6 ) 下标s 表示显参数( s e n s i b l e ) ，不包含化学反应引起的变化混合气体的其他 参数，如： e 。( y ，t ) ，弓。( y ，t ) ，矗。( y ，刃，r ( y ) 的定义与k ( y ，t ) 类似 8 第2 章湍流反应流的基本方程及模型 混合气体的总焓和总能为： 1 ， 日。( y ，t ) = k ( y ，t ) + ；u i 巩( 2 7 ) 毋( y ，t ) = e 。( y ，t ) + ；仉阢( 2 8 l 其中u 为流动速度 理想气体状态方程： p = p r t( 2 9 ) 其中p 为密度利用( 2 4 ) 式，可将上面的状态方程表示成e ，的表达式： p = p 兄考= 警帕娟l 1 ) 盹 ( 2 1 0 ) 7 ”捌= 詈= 南 ( 2 1 1 ) 设由于化学反应引起的组分。的质量变化速率为& ，则对均匀的混合气 体，质量分数k 满足： 鲁= s o ( y , p , 研 ( 2 1 2 ) 化学反应放热项电为： - 0 = 一s o ( y ，p ，t ) h o a ( 2 1 3 ) 2 2 雷诺平均的统计矩方程 本文只考虑远离壁面的高r e 数流动，守恒方程中分子输运的影响可以忽 略不计这时的雷诺平均的统计矩方程为： 妄( p ) + 去( 如) 反) = 。 击( ( 鹏) + 茜( ( 西玩岛+ 酬= 一南( 泐蛹) ( 2 1 4 ) 击( 赢) + 鑫慨( ( 西直+ 列= - 鑫 泐( 诵+ ；蛹+ 吃t 再) l 式e e - 表示质量( 脚r e ) 加权平均量， ) 表示雷诺平均量；p ，u = ( u t ，醍，u 3 ) t ，p ， e 分别是混合气体的平璁密度、平均速度、压强和比总能，脉动速度u = u o 。 显焓的脉动值= k h 。、 ? 方程组( 2 1 4 ) 等号右边的项，由p d f 输运方程提供 2 3 联合p d f 输运方程 9 平均状态方程可表示为： ( p ) = ( 7 1 ) ( p ) 度一；( 玩政+ 面磁) 比热比 ：生 毛 2 3 联合p d f 输运方程 ( 2 1 5 ) ( 2 1 6 ) p d f 方法的基本思想就是：用质量加权的几率密度函数( v ，霍；x t ) 描述 在时刻、坐标x = ( z l ，z 2 ，黝) 处的流动状态样本空间v = ( u ，v 2 ，n ) t 对应于 速度相空间的u = ( 巩，巩，u 3 ) t ，霍= ( 妒，也，c n + 1 ) t 对应于化学热力学参 数相空间的壬= ( 币1 ，咖2 ，+ 1 ) r ，其中m 为所涉及组分的总数在一个系 综中，p d f 可以用样本几率密度函数的平均值表示： “二 ( p ) ，( v ，皿，眺x ，) = p ( 皿) ( d ( u v ) 6 ( 圣一雪) d 一日) ) ( 2 1 7 ) 从化学热力学参数的守恒方程 警+ 篱= 一；篝+ & a 吐，肌+ - 和统计矩方程出发，可推出，( v ，霍；x 】t ) 的输运方程 4 s 】； 薯+ t 力峙若一等蓦+ 去) s o ) = 者( ( 一鬻+ 筹m 雪) 力+ 去( ( 篆l v ，雪) 力 ( 2 - 8 ) 方程( 2 1 8 ) 是标量，的输运方程方程等号左边的项以封闭的形式出现 可以精确求解，等号右边的项以条件平均值的形式出现，需要采用适当的模 型进行模拟方程第一项是，的时间导数项，第二项是，在物理空间的对流 项，第三项是由平均压力梯度引起的，在速度相空间的输运项，第四项是化 学反应引起的，在标量相空间的输运项右边三项分别由粘性应力张量勺、 脉动压力梯度母一a x 5 引起的，在速度相空间的输运和分子扩散通量圩( 表示 标量a 在z t 方向上的通量) 引起的，在组分相空间的输运 在求解p d f 方程时，采用l a g r a _ n g i a n 的观点，即流体颗粒的概念将几率密 度函数离散成由大量样本( j 睡粒) 所组成的一个系综，每一个样本的状态用一组 随机微分方程来描述每一个样本( 颗粒) 包含的信息有： m 。，x + ，u ，圣，矿) ， 分别表示颗粒的质量、空间位置、速度、化学热力学参数矢量、湍流频率( 上 1 0第2 章湍流反应流的基本方程及模型 标”表示该变量是颗粒的一个性质) 。求解随机微分方程即可得到几率密度 函数在空间变化的信息但在几率密度函数方程或随机微分方程中。脉动压 力梯度、粘性项和分子扩散通量的条件平均值项需要模拟另外求解时还需 湍流时间尺度的信息在后面的章节里，我们将介绍颗粒速度变化的模型、 分子混合模型和湍流频率模型必须说明的是这些模型不是最佳的，我们 只是用简单的模型来说明混合算法的求解思路 2 4 颗粒速度变化的模型 在p d f 方法中，流体颗粒速度u ( t ) 的演化通常用随机扩散方程模拟 p o v a 4 发展的 , a n g e v i n 模型可以相当好地模拟颗粒在速度相空间的变化，为 简单起见，我们这里仅采用简化的l a n g e v i n 模型( s i m p l i s e dl a n g e v i nm 0 d 哪： d w ( t ) = 一击磐出一( i 1 + ；岛) n ( w ( t ) 一玩冲+ ( 岛 。) 1 2 d i ( 2 1 9 ) 其中n 为密度加权的湍流频率的条件平均值，一定程度上反映了湍流间歇性 的影响，其定义为 湍流动能 n = 掣 ( 2 2 0 ) = 警 ( 2 2 1 ) 岛和c h 为模型常数( 表2 1 ) w ( o ) 为各向同性的w e i n e r 扩散过程，其变化量d 暇( t ) = 髓( 件d t ) 一矾( t ) 满足 联合正态分布且平均值为零( ( a w e ( t ) ) = 0 ) ，协方差矩阵( d 眦( z ) d ( ) ) = d t 幻 w e i n e r 过程中随机量 的取值范围是( 一o o ，+ o o ) ，该模型的使用将导致雷诺应 力趋于各向同性、湍流动能不断衰减正确模拟了脉动压力和粘性项引起的 各向同性化和衰变过程 与雷诺应力方程相比，上述简化的l a n g e v i n 模型等价于r o t t a 模型。由于 该模型不包含快速压力项的模型，导致求出的雷诺应力不完全符合实际，其 中的不足将在对计算结果的分析中进行讨论详细的l a n g e v i n 模型则考虑了 快速压力项的影响e 4 4 2 5 湍流频率模型 在速度和化学热力学参数的联台p d f 方程中， 度或时间尺度的信息为了使p d f 方程完全封闭， 不包含湍流模型所需的长 j p o p e 和c h e n 4 5 提出采用 2 6 湍流混合模型 速度一耗散率联合p d f 的方法以后，s l o o t e n 、j a y e s h 和p o p e 5 0 又提出湍 流频率的随机模型，建立速度一频率一标量联合p d f 方程 并采用条件平均频率来确定湍流的时间尺度，这相当于考虑了内部间歇 性的影响对速度一频率一标量联合p d f 模型，湍流频率是飘粒的一种属 性n 的值由方程( 2 2 0 ) 给出，当流动充分发展为均匀湍流时，n 等价于西 u ( t ) 的随机模型1 5 0 】为： 血d + ( ) = 一c 3 ( u 一面) n d t 一5 0 n u 出+ ( 2 c 3 c 4 & p w ) 1 1 2 d w( 2 2 2 ) 其中伤和。均为模型常数( 表2 1 ) 这里的w e i n e r 过程w ( t ) 与速度模型中的 w e i n e r 过程无关揣流频率源项咒的模化形式为： 咒= 2 一瓯1 齐 ( 2 2 3 ) 其中q 1 和吼2 为模型常数( 表2 1 ) ，p 为湍流产生项： p = 一诵筹( 2 2 4 ) 2 6 湍流混合模型 i 常数数值用途 岛 2 1s l m 0 6 8 9 3 n 的定义 1 0 5 6 湍流频率模型 c 0 2 0 9 湍流频率模型 伤 1 o 湍流频率模型 g 4 0 2 5 湍流频率模型 2 o i e m 混合模型 表2 1 ：模型常数 p d f 方程中标量分子输运的条件平均值项需要模拟，该项反映了分子输运 使湍流脉动减小的耗散作用，在有化学反应的流动中。由于分子输运和化学反应 之间有很强的耦合关系使得该项的模拟相当困难，目前提出的模型有颗粒相互 作用模型、i e m 模型和映射封闭模型等各种混合模型在文 5 l ，5 9 中有详细的 描述、分析和比较本文采用最简单的i e m ( i n t e x a c t i o n - b y - e x c h a n g e - w i t h - t h e - m e a n ) 模型假定颗粒化学热力学参数的变化满足t 。 一击，+ 11 = = 0 1 = 一；o n ( 壬( t ) 一圣) + s ( o ( t ) ) ( 2 2 5 ) 第2 章湍流反应流的基本方程及模型 其中瓯是模型常数。s 是化学反应源项，对第a 种组分的质量分数，是 化学反应速率；而对焓( a = 儿+ 1 ) ，则是化学反应放热速率； ，) 轧+ t = 一品蠼 ( 2 2 6 ) 其中 ：是参考温度下第a 中组分的生成热该式暗含了化学反应放热与温度 无关的假定 2 7 模化的联合p d f 方程 根据前面所列的模型，我们可以写出速度一频率一标量的联合p d f ( p ) ，( v ，雪，比x ，t ) = p ( 皿) 扣( u v ) d ( 圣一皿) d ( u 一日) ) ( 2 2 7 ) 的转运方程，其中v 、屯、0 分别是u 、壬、u 的相空间 一 在实际计算过程中，我们求解的是脉动速度一频率一标量的联合几率密 度函数( v ，田，0 ；x ，) ，v 是密度加权脉动速度u 的相空间颗粒速度”( t ) 中 的脉动部分u + ( t ) 成为颗粒的一个特征参数根据文献【4 4 中的方法，很容易 给出脉动速度“；( t ) 的演化方程： 酬= 丽1 锷字出一哼筹州j 1 + 轴喇懈酬2 帆( 2 2 8 ) 函数，和是完全等价的，所不同的是，口不包含平均速度疗的信息根据 方程( 2 2 2 ) 、( 2 2 5 ) 和( 2 2 8 ) 可以写出( v ，田，0 ；x ，t ) 的输运方程： 鬈鹕州差= 一鬻挈鬟一嚣掣+ c i l + 轫n 掣 + ；岛枷磊亳+ 品n 品( 口) + 岛n 蔷茵( 口一。) 】+ 岛q 赢( 口o ) + ；o n 最p ( 札一孔) 】一彘( 西) ( 2 _ 2 9 ) 在三维流动中，若有1 种组分，该几率密度函数方程的维数为+ 8 ( 时 间t ，空间坐标x ，脉动速度v ，地种组分的质量分数，焓h ) 如此高的维 数，很难用有限差分或者有限体积法求解如果用m o n t ec a r l o 方法求解，计、 算量只近似与样本空间的维数成线性关系，使得p d f 方程的数值求解成为可 能 第3 章化学反应模型 在第6 章给出的算例中的燃料都是丙烷 不同的化学反应机理，对燃烧过程的计算结果影响很大采用详细的化 学反应机理能够精确地描述火焰的结构，获取更精确的信息但在二维甚至 三维燃烧室流场的数值模拟中，采用详细的化学反应机理计算量将大大增 加本文的最终目的是发展一个混合算法的计算程序，为了程序调试易于进 行，本文采用简化的单步和两步的化学反应机理如果需要，计算程序很容 易被推广到详细的化学反应机理 3 1 单步反应机理 、 ? 丙烷和空气燃烧的反应机理可以用一个单步不可逆反应方程式来描述： 岛风+ 5 ( 仉+ 3 7 6 n 2 ) _ 3 c 0 2 + 4 1 1 2 0 + 1 8 8 2( 3 1 ) 这是最简单的丙烷空气化学反应机理w e s t b r o o k 等人 5 4 ，55 】已经发展了详 细的丙烷，空气化学反应机理，但详细的化学反应机理会使得计算量急剧增 大在很多情况下，单步反应模型也已经能满足湍流反应流的计算要求了 燃料质量分数方程的源项可用经验关系来近似： ，配、 r = a c 3 h s 8 0 2 】6 e x p 【_ 寺) ( 3 2 ) 、 其中m ，、m o ：分别为丙烷和氧气的质量分数；根据文献【5 5 ，系数a = 8 6 1 0 1 1 m d f 3 8 ，a = 0 1 ，b = 1 6 5 ，e = 1 2 6x1 0 5 j 。 根据方程式( 3 1 ) ，0 2 、c 0 2 、h 2 0 的反应速率为： p , o lj5 r i ，p , c o ：= 一3 r i ，r 凰0 = 一4 r f( 3 3 ) 3 2 两步反应机理 为了研究不同化学反应机理对计算结果的影响，我们也采用丙烷和空气 燃烧的多步反应机理模型 5 5 】，对同一算例进行数值求解，并与单步反应机理 的结果进行比较 在单步反应机理中，由于假定产物只有g 0 2 和皿o ，反应放热被高估了 实际上，在绝热温度( 一2 0 0 0 k ) 下，碳氢燃料的燃烧产物在平衡状态下有一定 量的c o 、1 1 2 以及h 、0 、o h 等自由基，致使用方程式( 3 1 ) 计算出的火焰 绝热温度偏高 1 4 第3 章化学反应模型 与单步反应机理相比较，两步反应机理多了有关c o 浓度的信息丙烷 和空气燃烧的两步反应方程式为：) 白凰+ 3 5 ( d 2 + 3 7 6 n 2 ) _ + 3 c 04 - 4 1 1 2 04 - 1 3 1 6 2 c o 4 - 0 5 ( d 2 + 3 7 6 2 ) + + 口0 2 + l t 8 8 n 2( 3 4 ) 设方程式( 3 4 ) 给出的两个反应的反应速率分别为r 1 和岛，则各组分的反应 速率为： 如h a = r i ，a c o = 恐一3 r i ，硒2 = 5 r i + 9 5 恐( 3 5 ) 其中： r t = 川岛剐6 唧( 一嘉) ， 岛= 岛+ 一r 2 一 = 如+ h 2 0 0 2 ，( 一等) 一缸i c o n f e x p ( 一而e 2 - ) ( 3 6 ) 其中常数 a l = 1 0x1 0 1 2 m “c m 3 s ，a 2 + = 1 0 1 4 6 m d l c m 3 s ，a 2 一= 5 i c m o l c m 3 8 口= o i 抽= 1 6 5 ，c = 1 ，d = 0 , 5 ，e = 0 2 5 ，= 1 e l = 1 2 6 1 0 5 j ，e 2 + = 砀一= 1 6 8x 1 0 5 j 第4 章数值计算方法 湍流燃烧的数值模拟精确与否，在很大程度上决定于湍流模型、燃烧模 型( 简化化学反应机理) 和数值计算方法三个方面湍流模型是为了使控制方 程封闭，简化化学反应机理则是为了减少化学反应计算所需要的时间，而精 确、有效的数值算法，对获得方程的精确解至关重要 一种理想的算法必须具备如下几个性能： 稳定； 能真实反映物理量的特性； 。 精确； 计算效率高； 算法与单个工况的初始条件、边界条件无关 对单个算法而言，很难全部满足上述要求，比如高效与精确就是一对相互排 斥、很难调和的矛盾一般算法有颗粒方法，特征线方法、有限差分( 体积) 法、有限元方法、谱方法等每一种方法都有其优势和缺陷，因此每一种方 法都有各自的适用范围随着对算法研究的深入，计算机的飞速发展，以及 人们对数值模拟的要求提高，融合了多个方法的优点的混合算法应运而生 本文用来求解p d f 方程的方法就是一种有限体积和颗粒方法的混合算法，该 混合算法具有有限体积法和颗粒方法的双重优势，更为重要的是两者之间形 成一种互补关系，最终达到一种精确和高效的。双赢”效果当然，由此而带 来的新问题是两种方法之间的相容性问题，详细的讨论将在下一章中进行 本章的内容包括：p d f 输运方程的整个求解思路；无结构网格技术，包 括颗粒的定位和不同算法之间的数据交换；有限体积法，包括方程的离散、 各项的求解算法和龙格一库塔法；颗粒方法，包括随机微分方程的求解；时 间步长的确定和边界条件 4 1 基本思路 p d f 输运方程与传统的统计矩方程有着完全不同的形式和结构，是一个 高维的标量输运方程由于其维数过高，无法用有限体积或者有限差分方法 求解p d f 方法的发展历史表明，m o n t e - c a r l o ( 颗粒) 方法是求解联合p d f 方 程的一种行之有效的方法，因为m o n t e - c a r l o 方法的计算量只随维数的增多而 线性增加联合p d f 方法也已经成功地应用于许多实际流动闻题的计算由 于m o n t e - c a x l o 方法一般都需要大量的颗粒以降低数值偏差和统计误差。所以 相对来说，计算量都偏大人们采用各种不同的数值解法，来求解联合p d f 方程( 如图4 1 ) 1 6 第4 章数值计算方法 * 一e 模型 混合算法1 ；( p d f 2 d s ) ：混合算法2： i p d f 2 d - f v ，f v m c 2 d l l 一j 口i i 确 = = = ：掖拟 图4 1 ：各种不同的p d f 方程求解方法 单纯的颗粒方法 ；( p d f 2 d v ) 研究发现，当采用单纯的颗垃网格方法求解联合p d f 模型方程时，方 程( 2 1 9 ) 中的右边第二项是计算偏差的主要来源该项的作用是使颗粒的速 度w 向平均速度取( 从颗粒场获得) 松弛当从颗粒场获得平均速度场误差 较大时。往往导致计算很难收敛解决这一问题的办法是用脉动速度的随机 微分方程( 2 2 8 ) 代替平均速度的随机微分方程( 2 1 9 ) ，这时颗粒的脉动速度嵋 严 向0 松弛，而平均速度则通过有限体积法求解雷诺平均e u l e r 方程( 2 1 4 ) 来获 得，方程( 2 1 4 ) 右边的项由颗粒方法来提供 图( 4 2 ) 是本文发展的程序f v m c 2 d 的流程图首先，对有限体积法和 颗粒方法的数据进行初始化；然后有限体积法向前计算一个时间步长，并将 计算所得的平均速度场和平均压力场作为m o n t e - c a r l o 方法的输入；颗粒方法 中，颗粒在物理空问、速度相空间，频率相空间、标量相空间也向前计算一 个时间步长为了获得因化学反应引起的颗粒在标量相空问的变化，每个颗 粒都需调用计算化学反应速率的程序；从颗粒场可以获得雷诺应力、湍流通 量、能量源项和寺，；为了减小统计误差和偏差，对这些量进行了时间平均；经 过时间平均的量将作为有限体积法的输入，开始新的一个时间步长的计算 重复这一操作，直到计算结果收敛 4 1 基本思路 1 7 图4 2 ：f v m c 2 d 程序的流程图 1 8第4 章数值计算方法 4 2 无结构网格 本文所采用的是无结构网格在流体力学计算中，通常采用矩形的结构 网格所谓结构网格，就是区域内离某节点最近的网格能够通过简单增加或 减小序号来获得这种网格技术在遇到较复杂的几何形状时，就很难解决边 界问题了如果区域的边界对边界附近的流场影响很强，那么边界的处理方 法会对整个流场的计算结果产生很大的影响因此用矩形网格在处理不规则