有限差分法

《有限差分法》由会员分享，可在线阅读，更多相关《有限差分法（23页珍藏版）》请在装配图网上搜索。

1、(11.1.51)(11.1.52)(11.1.53)于是 n (X)=CneIkxkk=s n+1 (x)=艺C n+1e Ikxkk=s n = g n (X ) =Cne叫iikk=s n+1 = n+1 (X ) =C n+1e Ikx.iikk =g(11.1.54)有限差分法傅立叶稳定性分析分析差分格式稳定性的方法很多，大部份应用于线性方程，这里只介绍其中最常用的 一种：傅立叶稳定性分析法。傅立叶稳定分析法由 Von Neumann 于 20 世纪 40 年代提出，所以又称为 Von Neumann 稳定性分析法。该方法的基本思想是，将解的误差作周期延拓并用傅立叶级数表示出来，然

2、后考察每一个傅立叶级数分量的增大和衰减情况。如果每一分量的强度(或振幅)是随时间 的推移而增大的，则所讨论的差分格式是不稳定的；反之，若每一分量的振幅是随时间的推 移而衰减或保持不变，则格式是稳定的。为了进行这种分析，可以把某一分量的表达式代入 到误差传播方程中，得出相邻二时间层间该分量的振幅比，通常称为放大因子。稳定性的条 件要求放大因子的绝对值(或模)小于或等于1。当放大因子等于1 时，称为中性稳定，在 这种情况下任何时刻引进的误差都不会衰减或放大。【例11.1E】讨论逼近以下一维对流方程的FTCS格式的稳定性:dudu+ a = 0 a 0dtdx该方程的 FTCS 格式为u n +1

3、u nu n u ni1 + a i+1i1 = 0At2AX将式(11.1.52)改写成易于递推计算的差分格式，有a .()Un+1 = un九 Un Un /ii 2i+1i1式中，九二At/(Ax)为网格比。相应于上式的误差传播方程为a .() n+1 = g n A g n n 丿ii 2i+1i1式中，是各节点上的离散量。如果对在正负方向上作周期延拓，即把看作是以某一定 值为周期的周期函数，则 n、 n+1可以展开为以下的傅立叶级数： n = n (X + Ax) =Cne Ik (x, 土Ax)i 1ik(11.1.55)k=-g n+1 = n+1 (x Ax)=Cn+1eIk(

4、xiAx)i 1ikk=-g其中，I = 1。将式(11.1.54)和(11.1.55)代入式(11.1.53)并整理，得Cn+1 C nkkk=-g1 _ Xt? IkAx e IkAr 丿 | ” Ikxi 二 0(11.1.56)式(11.1.56)相当于将零展开成傅立叶级数；式中 内相当于傅立叶系数，它对于所有的 k 都应该等于零：(11.1.57)(11.1.58)C n+1 C n 1 X (e IkAx e IkAx ) = 0G =1 C X (e 皿2e IkAx )k k 2则式(11.1.57)成为：(不失一般性，去掉式中的下标记号 k)Cn+1 = GCn(11.1.5

5、9)表示误差从第 n 层传播到第 n+1 层时，以傅立叶级数表示的每一误差分量的振幅放大或衰 减了 G倍。所以，称G为放大因子，傅立叶稳定性分析法就集中在对G的分析上。如果IGI 1,则误差的振幅随n的增大而增长，差分格式不稳定；如果IGI W 1,则误差的振幅随n 的增大而减小或不变，差分格式稳定。应用欧拉公式e土二cosZ 土 IsinZ，将式(11.1.58)改写为G=1 IX sin kAxIGI2=1+2X2sin2 kAx当sin kAx丰0时，不管怎样选取网格比久总有IGI 1,因此差分格式(11.1.52)是不稳定的。从上例的分析注意到，以傅立叶稳定性分析法判断格式稳定性时，是

6、从误差传播方程(11.1.53)出发、将计算节点的误差 n延拓为傅立叶级数、并通过分析式(11.1.56沖傅立叶级i数任一系数，来确定放大因子G进而确定差分格式的稳定性的。对于齐次线性微分方程，由于误差传播方程与其相应的差分方程形式相同，在傅立叶 稳定性分析中，只要令Un = C ne Ikx. iUn = Cne Ik (x, 土 Ax) I i 1Un1 二 Cn1e IkxiiUn 1 = C n 1 e Ik (x, Ax) i1(11.1.60)并将它们代入相应的差分格式中，同样可以得到与上例相同的放大因子G的表达式。为方 便叙述，在以后的傅立叶稳定性分析讨论中将采用式(11.1.6

7、0)的方式。【例11.1F】讨论逼近一维对流对方程(11.1.51)的FTBS格式的稳定性。 相应于方程(11.1.51 )的FTBS格式为+ =0AtAx(11.1.61)与上式等价的易于递推计算的差分格式为(11.1.62)1 t )Un+1 = un X Un Un /i iii 1将式(11.1.60)代入上式，有Cn+1eIkx, = CneIkx, X Cn e-IkAt eIkx,.消去公因子e叫后，得由此得放大因子为Cn+1 = 1 a九 1 e - IkAxCnG = 1 a九(一e -】kA )= 1 cos kAx) la九 sin kAxIGI2 = 1 -al(1 -

8、coskAxH2 +a2尢2 sin2 kAx=1 - 2a尢sin2 空丫1- sin2 空k2丿k2丿=1 - 4a尢G - a尢)sin 22如果00(11.1.63)dtdx2相应的 BTCS (Backward Time Central Space)隐式差分格式，为un+1 -unun+1 - 2un+1 + un+1ii- _ 卩一i= 0At(Ax)2将其变形为易于递推计算的差分格式，有(11.1.64)-pXu n+1 + (1 + 2pX)u n+1 J3Xu n+1 = u ni -1ii +1i式中，X = At/(Ax)2为网格比。将式(11.1.60)代入并消去公因子

9、e叫，易得放大因子=1/(1 + 4 PX sin 2kAx亍由于卩0,所以对任何网格比X都有|G| W 1,所以隐式差分格式(11.1. 64)是稳定的。从以上几个例子看到，差分格式(11.1. 61)在aX W1条件下是稳定的；差分格式(11.1. 64) 对任何网格比X = At/(Ax)2都是稳定的；差分格式(11.1. 52)对任何网格比X = At/ Ax都是 不稳定的。我们称第一种情况的差分格式是条件稳定的，第二种情况的差分格式是无条件或 绝对稳定的，第三种情况的差分格式则为绝对或无条件不稳定的。在实际应用中对于差分格式的选择，首先要排除不稳定的差分格式，还要寻求对稳定性 限制较

10、弱的差分格式，当然最好是选择无条件稳定的差分格式。1.1 流体力学的基本概念1.1.1流体的连续介质模型流体质点(fluid particle)：几何尺寸同流动空间相比是极小量，又含有大量分子的微元体。连续介质(continuum/continuous mediu m)：质点连续地充满所占空间的流体或固体。 连续介质模型(continuum/continuous medium model)：把流体视为没有间隙地充满它所占 据的整个空间的一种连续介质，且其所有的物理量都是空间坐标和时间的连续函数的一种假 设模型： u =u(t,x,y,z)。1.1.2 流体的性质1. 惯性惯性(fluid in

11、ertia)指流体不受外力作用时，保持其原有运动状态的属性。惯性与质量有 关，质量越大，惯性就越大。单位体积流体的质量称为密度(density),以r表示，单位为kg/m3。 对于均质流体，设其体积为V质量为m,则其密度为(1-1)mp= V对于非均质流体，密度随点而异。若取包含某点在内的体积AV，其中质量Am，则该 点密度需要用极限方式表示，即(1-2)Amp = limAV0 AV2. 压缩性作用在流体上的压力变化可引起流体的体积变化或密度变化，这一现象称为流体的可压缩性。压缩性(compressibility)可用体积压缩率k来量度k=dV / V dp / pdpdp(1-3)式中：

12、p 为外部压强。 在研究流体流动过程中，若考虑到流体的压缩性，则称为可压缩流动，相应地称流体为 可压缩流体，例如高速流动的气体。若不考虑流体的压缩性，则称为不可压缩流动，相应地 称流体为不可压缩流体，如水、油等。3. 粘性粘性(viscosity)指在运动的状态下，流体所产生的抵抗剪切变形的性质。粘性大小由粘 度来量度。流体的粘度是由流动流体的内聚力和分子的动量交换所引起的。粘度有动力粘度 卩和运动粘度v之分。动力粘度由牛顿内摩擦定律导出：duT = U - dy(1-4)式中：t为切应力，Pa；卩为动力粘度，Pa - s； du/dy为流体的剪切变形速率。运动粘度与动力粘度的关系为(1-5)

13、v=式中： v 为运动粘度， m2/s。 在研究流体流动过程中，考虑流体的粘性时，称为粘性流动，相应的流体称为粘性流体； 当不考虑流体的粘性时，称为理想流体的流动，相应的流体称为理想流体。根据流体是否满足牛顿内摩擦定律，将流体分为牛顿流体和非牛顿流体。牛顿流体严格 满足牛顿内摩擦定律且卩保持为常数。非牛顿流体的切应力与速度梯度不成正比，一般又分 为塑性流体、假塑性流体、胀塑性流体3 种。塑性流体，如牙膏等，它们有一个保持不产生剪切变形的初始应力T0,只有克服了这个 初始应力后，其切应力才与速度梯度成正比，即(1-6)(1-7)(1-8)(1-13)duT =T + H -假塑性流体，如泥浆等，

14、其切应力与速度梯度的关系是0dy胀塑性流体，如乳化液等，其切应力与速度梯度的关系是(du dy丿1.1.4 流体运动的描述1. 流体运动描述的方法描述流体物理量有两种方法，一种是拉格朗日描述；一种是欧拉描述。拉格朗日(Lagrange)描述也称随体描述，它着眼于流体质点，并将流体质点的物理量认 为是随流体质点及时间变化的，即把流体质点的物理量表示为拉格朗日坐标及时间的函数。 设拉格朗日坐标为(a,b,c)，以此坐标表示的流体质点的物理量，如矢径、速度、压强等等在 任一时刻t的值，便可以写为a、b、c及t的函数。若以f表示流体质点的某一物理量，其拉格朗日描述的数学表达式为f = f (a, b,

15、 c, t)(1-15)例如，设时刻t流体质点的矢径即t时刻流体质点的位置以r表示，其拉格朗日描述为r = r (a, b, c, t)(1-16)同样，质点的速度的拉格朗日描述是v = v(a, b, c, t)(1-17)欧拉描述，也称空间描述，它着眼于空间点，认为流体的物理量随空间点及时间而变化， 即把流体物理量表示为欧拉坐标及时间的函数。设欧拉坐标为(q,q2，q3)，用欧拉坐标表示的 各空间点上的流体物理量如速度、压强等，在任一时刻t的值，可写为qq2、q3及t的函 数。从数学分析知道，当某时刻一个物理量在空间的分布一旦确定，该物理量在此空间形成 一个场。因此，欧拉描述实际上描述了一

16、个个物理量的场。若以f表示流体的一个物理量，其欧拉描述的数学表达式是(设空间坐标取用直角坐标) f=F(x,y,z,t)=F(r,t)(1-18)如流体速度的欧拉描述是2. 拉格朗日描述与欧拉描述之间的关系拉格朗日描述着眼于流体质点，将物理量视为流体坐标与时间的函数；欧拉描述着眼于 空间点，将物理量视为空间坐标与时间的函数。它们可以描述同一物理量，必定互相相关 设表达式f = f (a,b,c,t)表示流体质点(a,b,c)在t时刻的物理量；表达式f = F(x,y,z,t)表示 空间点(x,y,z)在时刻t的同一物理量。如果流体质点(a,b,c)在t时刻恰好运动到空间点(x,y,z) 上，则

17、应有x= x(a,b,c,t)(1-20)(1-21)(1-22)(1-23) y = y(a, b, c, t)z = z (a, b, c, t)F(x,y,z,t)=f(a,b,c,t)事实上，将式(1-16)代入式(1-21)左端，即有F(x,y,z,t)= Fx(a,b,c,t),y(a,b,c,t),z(a,b,c,t),t= f (a, b, c, t)或者反解式(1-16)，得到a = a(x, y, z,t)b = b(x, y, z,t)c = c(x, y, z,t)将式(1-23)代入式(1-21)的右端，也应有f (a,b,c,t) = f a(x, y, z,t),

18、b(x, y, z,t),c(x, y, z,t),t(1-24)= F(x, y,z,t)由此，可以通过拉格朗日描述推出欧拉描述，同样也可以由欧拉描述推出拉格朗日 描 述。3. 随体导数流体质点物理量随时间的变化率称为随体导数(substantial derivative)，或物质导数、质 点导数。按拉格朗日描述，物理量f表示为f = f (a, b, c, t)，f的随体导数就是跟随质点(a,b,c)的 物理量f对时间t的导数df /&。例如，速度v(a,b,c,t)是矢径r(a,b,c,t)对时间的偏导数，(1-25)v(a, b, c, t)=业心即随体导数就是偏导数。按欧拉描述，物理

19、量f表示为f = F(x,y,z,t)，但dF/dt并不表示随体导数，它只表示 物理量在空间点(x,y,z,t)上的时间变化率。而随体导数必须跟随t时刻位于(x,y,z,t)空间点 上的那个流体质点，其物理量f的时间变化率。由于该流体质点是运动的，即x、y、z是变 的，若以a、b、c表示该流体质点的拉格朗日坐标，则x、z将依式(1-16)变化，从而f =F(x,y,z,t) 的变化依连锁法则处理。因此，物理量f=F(x,y,z,t)的随体导数是OF dx dF dy dF dzOFDF(x，y,*, = df x(a, b, c, t), y(a, b, c, t), z (a, b, c,

20、t), t Dt(1-26)+ + +Ox OtOy OtOz OtOtOFOFOFOF式中：D/Dt表示随体导数。从中可以看出，对于质点物理量的随体导数，欧拉描述与拉格朗日描述大不相同。前者 是两者之和，而后者是直接的偏导数。4. 定常流动与非定常流动 根据流体流动过程以及流动过程中的流体的物理参数是否与时间相关，可将流动分为定 常流动(steady flow)与非定常流动(unsteady flow)。定常流动：流体流动过程中各物理量均与时间无关，这种流动称为定常流动。 非定常流动：流体流动过程中某个或某些物理量与时间有关，则这种流动称为非定常流 动。5. 有旋流动与有势流动由速度分解定理

21、，流体质点的运动可以分解为：(1) 随同其他质点的平动；(2) 自身的旋转运动；(3) 自身的变形运动(拉伸变形和剪切变形)。(1-31)在流动过程中，若流体质点自身做无旋运动(irrotational flow),则称流动是无旋的，也 就是有势的，否则就称流动是有旋流动(rotational flow)。流体质点的旋度是一个矢量，通常 用表示，其大小为ijk1dddw 2dxdydzuvw若=0，则称流动为无旋流动，否则就是有旋流动。与流体的流线或迹线形状无关；粘性流动一般为有旋流动；对于无旋流动，伯努利方 程适用于流场中任意两点之间；无旋流动也称为有势流动(potential flow)，

22、即存在一个势函 数申(x, y, z, t)，满足:V = grad 申(1-32)U 9dx(1-33)w dz1.2 CFD 基本模型流体流动所遵循的物理定律，是建立流体运动基本方程组的依据。这些定律主要包括质 量守恒、动量守恒、动量矩守恒、能量守恒、热力学第二定律，加上状态方程、本构方程 在实际计算时，还要考虑不同的流态，如层流与湍流。1.2.1 基本控制方程1. 系统与控制体在流体力学中，系统是指某一确定流体质点集合的总体。系统以外的环境称为外界。分 隔系统与外界的界面，称为系统的边界。系统通常是研究的对象，外界则用来区别于系统 系统将随系统内质点一起运动，系统内的质点始终包含在系统内

23、，系统边界的形状和所围空 间的大小可随运动而变化。系统与外界无质量交换，但可以有力的相互作用，及能量(热和功)交换。 控制体是指在流体所在的空间中，以假想或真实流体边界包围，固定不动形状任意的空 间体积。包围这个空间体积的边界面，称为控制面。控制体的形状与大小不变，并相对于某 坐标系固定不动。控制体内的流体质点组成并非不变的。控制体既可通过控制面与外界有质 量和能量交换，也可与控制体外的环境有力的相互作用。2. 质量守恒方程(连续性方程)在流场中，流体通过控制面A流入控制体，同时也会通过另一部分控制面A2流出控制 体，在这期间控制体内部的流体质量也会发生变化。按照质量守恒定律，流入的质量与流出

24、 的质量之差，应该等于控制体内部流体质量的增量，由此可导出流体流动连续性方程的积分 形式为fffp dxdydz 花 M p v - ndA = 0(1-38)dtVA式中：V表示控制体，A表示控制面。等式左边第一项表示控制体V内部质量的增量；第二 项表示通过控制表面流入控制体的净通量。根据数学中的奥-高公式，在直角坐标系下可将其化为微分形式：(1-39)(1-40)(1-41)(1-42)Qp 丄 d(pu)丄 d(pv)丄d(pw)0+ u+ v+ w= 0QtQxQyQz对于不可压缩均质流体，密度为常数，则有Qu Qv Qw + + = 0Qx Qy Qz对于圆柱坐标系，其形式为Qp ,

25、 pv3(pv )丄 Q(PvJ ,(pv )0QtrQrrQOQz对于不可压缩均质流体，密度为常数，则有v Qv Qv Qv-丄-丄 8+- = 0r Qr rQ0 Qz3. 动量守恒方程(运动方程)动量守恒是流体运动时应遵循的另一个普遍定律，描述为：在一给定的流体系统，其动 量的时间变化率等于作用于其上的外力总和，其数学表达式即为动量守恒方程，也称为运动 方程，或 N-S 方程，其微分形式表达如下：du 厂 QpQPQpp = p F + xr + 丄zxdtbx QxQyQzdvQpQpQp(1-43)p = p F +xr +yy 丄:dtbyQxQyQzdwdtpFbzQpQp+ y

26、z丄Qp式中： F 、F 、F 分别是单位质量流体上的质量力在三个方向上的分量； p 是流体内应 bx by bzyx力张量的分量。动量守恒方程在实际应用中有许多表达形式，其中比较常见的有如下几种。(1) 可压缩粘性流体的动量守恒方程(2)(3)(4)(5)dudp d II Ip花=p fx +玉+玉I屮(du dv + ( dy dx 丿du 2(du dv dw 一 + 一 + 一 dx 3 (dx dy dz(dw du)I dx dz 丿d+ LXdzdv” dp 6 I I _p 応=p f w+iy f|2dv dw、+ ldz dy 丿ddzp也=pf +里+AL 2比 dtz

27、dz dz Idzdv 2( du dv dwdy 3、dx dy dz 丿(du dv )+ (dy dx 丿一dw 2 ( du dv dw-+ + 3 ( dx dy dz(1-44)dzdI(dw du)d(dv dw )一 + 一+ 一Xdxdx dz 丿dz(dzdz 丿_常粘性流体的动量守恒方程dv卩p = p F - gradp + grad(divv) + 凹2v dt3常密度常粘性流体的动量守恒方程dv p= p F - gradp +2 vdt无粘性流体的动量守恒方程(欧拉方程)p 芈=p F - gradpdt静力学方程(1-45)(1-46)(1-47)(1-48)相

28、对运动方程p F = gradp(6)在非惯性参考系中的相对运动方程是研究像大气、海洋及旋转系统中流体运动的所必须考虑的。由理论力学得知，绝对速度v为相对速度v及牵连速度v之和，即 arev = v + v(1-49)a r e其中，v = v +Ox r，v为运动系中的平动速度，。是其转动角速度，r为质点矢径。e 0 0(1-50)而绝对加速度a为相对加速度a、牵连加速度a及科氏加速度a之和，即 areca =a +a +aarec其中，a = 0 +x r + Q x (Ox r)， a = 20 x v。e dt dtcr将绝对加速度代入运动方程，即得到流体的相对运动方程(1-51)p=

29、 p F + divP - a - 2Q vdtbc r4. 能量守恒方程将热力学第一定律应用于流体运动，把式(1-51)各项用有关的流体物理量表示出来，即 是能量方程。如式(1-52)所示。ddd(PE) +u (pE + p)=-Ctex iexii式中：e = h -上+傑；k是有效热传导系数， p 2 effk 工 h J + u (t )eff Cxjj jj j j effi j(1-52)k = k + k，其中k是湍流热传导系数，根 eff t t据所使用的湍流模型来定义； J 是组分 j 的扩散流量； S 包括了化学反应热以及其他用户 jh定义的体积热源项；方程右边的前3 项

30、分别描述了热传导、组分扩散和粘性耗散带来的能量 输运。1. 湍流模型分类湍流流动模型很多，但大致可以归纳为以下3 类。第一类是湍流输运系数模型，即将速度脉动的二阶关联量表示成平均速度梯度与湍流粘 性系数的乘积，用笛卡儿张量表示为_pUu = |Lli jt(1-53)Cu Cu i + jCx Cxji 模型的任务就是给出计算湍流粘性系数 l 的方法。根据建立模型所需要的微分方程的 t 数目，可以分为零方程模型(代数方程模型)、单方程模型和双方程模型。第二类是抛弃了湍流输运系数的概念，直接建立湍流应力和其他二阶关联量的输运 方 程。第三类是大涡模拟。前两类是以湍流的统计结构为基础，对所有涡旋进

31、行统计平均。大 涡模拟把湍流分成大尺度湍流和小尺度湍流，通过求解三维经过修正的 Navier-Stokes 方程 (纳维-斯托克斯方程，简称 N-S 方程)，得到大涡旋的运动特性，而对小涡旋运动还采用上 述的模型。实际求解中，选用什么模型要根据具体问题的特点来决定。选择的一般原则是精度要高 应用简单，节省计算时间，同时也具有通用性。Fluent提供的湍流模型包括:单方程(Spalart-Allmaras)模型、双方程模型(标准k-&模型、 重整化群k - 8模型、可实现k - 8模型)及雷诺应力模型和大涡模拟，如图1-1所示。基于RAN$的模型/零方裡模型I快方程摸葩會多理理 包必物机取方程模

32、叽:标准肚模型!；重整化群肚模型!I可实现模型!雷诺说刀模塑宜接数值模祕Fluen礙供的计算摸唱图 1-1 湍流模型详解2. 平均量输运方程雷诺平均就是把 Navier-Stokes 方程中的瞬时变量分解成平均量和脉动量两部分。对于 速度，有u = u + u(1-54)i i i式中：u和u分别是平均速度和脉动速度(i = 1,2,3 )。ii类似地，对于压力等其他标量，也有0=$+炉(1-55)式中：Q表示标量，如压力、能量、组分浓度等。把上面的表达式代入瞬时的连续与动量方程，并取平均(去掉平均速度上的横线)，可以 把连续与动量方程写成如下的笛卡儿坐标系下的张量形式：dt+i(p u)=0

33、dugp gp r 二 + dtgxgxijidugu2。duLI r HjO iI gxgx3 j gxjilg+ 一gxj/ 、-pu u k i j丿(1-56)(1-57)上面两个方程称为雷诺平均的Navier-Stokes(RANS)方程。它们和瞬时Navier-Stokes方程有相同的形式，只是速度或其他求解变量变成了时间平均量。额外多出来的项-Puu 是 ij 雷诺应力，表示湍流的影响。对于密度变化的流动过程，如燃烧问题，需要采用法夫雷(Favre)平均才可以求解。法夫 雷平均就是除了压力和密度本身以外，所有变量都用密度加权平均。变量的密度加权平均定 义如下：& 二 p / p(

34、1-58)式中：符号表示密度加权平均，对应于密度加权平均值的脉动值用表示，有=& +0。显然，这种脉动值的简单平均值不为零，但它的密度加权平均值等于零，即帀工0,丽二0。为了求解方程(1-57),必须模拟雷诺应力项以使方程封闭。通常的方法是应用Boussinesq假设，认为雷诺应力与平均速度梯度成正比，表达式如下：Tgxgxj i 丿Boussinesq假设被用于单方程模型和k -双方程模型。这种近似方法好处是与求解湍流gugur +gut gxi 丿iOij(1-59)粘性系数有关的计算时间比较少。例如，在 Spalart-Allmaras 单方程模型中只多求解一个表 示湍流粘性的输运方程；

35、在k - 双方程模型中只需多求解湍动能k和耗散率两个方程，湍 流粘性系数用湍动能k和耗散率的函数来描述。Boussinesq假设的不足之处是假设L(是个 各向同性标量，对于一些复杂流动，该条件并不是严格成立，所以具有其应用局限性。另外的近似方法是求解雷诺应力各分量的输运方程。这也需要额外再求解一个标量方程，通常是耗散率方程。这就意味着对于二维湍流流动问题，需要多求解4个输运方程 而三维湍流问题需要多求解7个方程，需要较多的计算时间，要求更高的计算机内存。 在很多情况下基于 Boussinesq 假设的模型很好用，而且计算量并不是很大。但是，如 果湍流场各向异性很明显，如强旋流动以及应力取得的二

36、次流等流动中，求解RSM模型可 以得到更好的结果。3. 常用湍流模型简介1)单方程(Spalart-Allmaras)模型单方程模型求解变量是V ,表征出了近壁(粘性影响)区域以外的湍流运动粘性系数。V的 输运方程为dv1p = G + dtv -vddxj式中：G是湍流粘性产生项；Y是由于壁面阻挡与粘性阻尼引起的湍流粘性的减少；j和-dv (卩 + P v k dxj汁Cb 2 dx(1-60)vvvC 是常数； v 是分子运动粘性系数。b2湍流粘性系数p = pvf，其中，f是粘性阻尼函数，定义为f =tv1v1v1 X 3 + C 3v1v f，f = 1 -k2d 2 v2 v21+

37、X fv11 ( du du ji2 dx dxI i j丿v1而湍流粘性产生项G模拟为G = C pSv，其中/三S +vvb1是常数，d是计算点到壁面的距离；S三*2QQ，ij ijQ =ij， C 和 kb1。在 Fluent 软件中，考虑到平均应变率对湍流产生也起到很大作用，S三IQ I+C min(0,l S I -1 Q I),其中，ij prodijij1 (du du j + I2 dx dxij在涡量超过应变率的计算区域计算出来的涡旋粘性系数变小。这适合涡流靠近涡旋中心 的区域，那里只有“单纯”的旋转，湍流受到抑止。包含应变张量的影响更能体现旋转对湍 流的影响。忽略了平均应变

38、，估计的涡旋粘性系数产生项偏高。( Yd丿C =2.0，IQ I三、2Q Q ，IS I三 2S S，平均应变率S =4prodijj jjj jj湍流粘性系数减少项Y =C p fvw1w，其中，vr三Sk 2 d 2g = r + C (r 6 - r) ，w2率的影响。上面的模型常数在 FluentC、w1C 、 C 是常数，在计算 r1 + C 6w3、g + C 66w3时用到的 Sv 受平均应变1/ 6w 2w3软件中默认值为C = 0.1335， Cb1b 2C = 2.0， k = 0.41 。w3二 0.622，- = 2/3，vC =7.1， C =C /k2 +(1+C

39、)/- ， C =0.3v1w1b1b 2vvw22)标准k - s模型标准k -s模型需要求解湍动能及其耗散率方程。湍动能输运方程是通过精确的方程推 导得到的，但耗散率方程是通过物理推理，数学上模拟相似原形方程得到的。该模型假设流 动为完全湍流，分子粘性的影响可以忽略。因此，标准k -s模型只适合完全湍流的流动过 程模拟。标准k -s模型的湍动能k和耗散率s方程为如下形式：( 叮I 3丿P + -T -丿_ dkdp =dt dxidsdp =一dtdxidkdxi处 + C 二(G + C G ) - C p dxi+ G + G -ps- YkbMss21s kk3s b2sk(1-61

40、)(1-62)式中：G表示由于平均速度梯度引起的湍动能产生，G表示由于浮力影响引起的湍动能产 kbk2生；Y表示可压缩湍流脉动膨胀对总的耗散率的影响。湍流粘性系数卩=pC。t 卩s= 0.09，湍动能 k 与耗散M在 Fluent 中，作为默认值常数率8的湍流普朗特数分别为-=1.0，k3)重整化群k - 模型C =1.44，C =1.92， C1s2s3s- =1.3。s重整化群k-s模型是对瞬时的Navier-Stokes方程用重整化群的数学方法推导出来的模 型。模型中的常数与标准k -s模型不同，而且方程中也出现了新的函数或者项。其湍动能与耗散率方程与标准k - 模型有相似的形式：Qkd

41、k Q、Qk 厂 厂yp = （a p ） + G + G p Y dt Qxid Q p =一 dt Qxi式中： G 表示由于平均速度梯度引起的湍动能产生， G 表示由于浮力影响引起的湍动能产 kb生；Y表示可压缩湍流脉动膨胀对总的耗散率的影响，这些参数与标准k -模型中相同。aMk 和a分别是湍动能k和耗散率的有效湍流普朗特数的倒数。湍流粘性系数计算公式为 k effQxi、Qkb(1-63)(a 卩)+ C (G + C G ) - C p ReffQxi1 kk3 b2k(1-64)d=1.72 dV，其中，v = p /p， C 100。对于前面方程的积分，可以3 -1 - Cve

42、ffv精确到有效雷诺数（涡旋尺度）对湍流输运的影响，这有助于处理低雷诺数和近壁流动问题的k2模拟。对于高雷诺数，上面方程可以给出：p =pC ，C = 0.0845。这个结果非常有意t p p思，和标准k -模型的半经验推导给出的常数C = 0.09非常近似。在Fluent中，如果是默 p认设置，用重整化群k - 模型时是针对的高雷诺数流动问题。如果对低雷诺数问题进行数 值模拟，必须进行相应的设置。4）可实现k - 模型可实现k - 模型的湍动能及其耗散率输运方程为dkQp =dtQxidQp = dtQxi卩)Qk P + r I I 3丿 ( 叮 I Q丿 + G + G - p - Yk

43、bMQxi+ pC S pC Qxi2(1-65)(1-66)式中：C = max 0.43,在上述方程中， G 表示由于平均速度梯度引起的湍动能产生， G 表示由于浮力影响引kb起的湍动能产生；Y表示可压缩湍流脉动膨胀对总的耗散率的影响；C和C是常数；Q和M21kQ 分别是湍动能及其耗散率的湍流普朗特数。在 Fluent 中，作为默认值常数， C =1.44，1C =1.9 ， Q =1.0 ， Q =1.2 。2k该模型的湍流粘性系数与标准k -模型相同。不同的是，粘性系数中的C不是常数，P而是通过公式计算得到C =1 ，其中，U* =、：SS + Q Q，Q = Q - 2 W，PU *

44、 Kij ijij ijijijijk kA + A0 s 0 = G + 2 3 , d表示在角速度W旋转参考系下的平均旋转张量率。模型常数 ij；6cosQ ，Wijijijk kA= 4.04A=0sS1QuQu、J+ 1ij2QxQx Ik1S S S0 = arccosG 6W）， 式中 W = i k，S 三 jS S ，3Sj j。从这些式子中发现， C 是平均应变率与旋度的函数。在平衡边界层惯V ij丿性底层，可以得到C = 0.09，与标准k-模型中采用的常数一样。P 该模型适合的流动类型比较广泛，包括有旋均匀剪切流、自由流（射流和混合层）、腔道 流动和边界层流动。对以上流动

45、过程模拟结果都比标准k-模型的结果好，特别是可实现 k- 模型对圆口射流和平板射流模拟中，能给出较好的射流扩张角。双方程模型中，无论是标准k - 模型、重整化群k - 模型还是可实现k - 模型，三个模型有类似的形式，即都有k和的输运方程，它们的区别在于：计算湍流粘性的方法不同；控制湍流扩散的湍流普朗特数不同；方程中的产生项和Gk关系不同。但都包含 了相同的表示由于平均速度梯度引起的湍动能产生G，表示由于浮力影响引起的湍动能产 k生G ;表示可压缩湍流脉动膨胀对总的耗散率的影响Y。 bM 湍动能产生项 QuG = p u u jk i j QxiG =P g卑空b i P Qxrt i(1-6

46、7)(1-68)式中：P是能量的湍流普特朗数，对于可实现k -模型，默认设置值为0.85；对于重整化rt1群 k - 模型，P = 1/a， a = 1/ P = k / |LiC。热膨胀系数 p=rtrtp浮力引起的湍动能产生项变为，对于理想气体pc 卩QpG = - g 厂b i p PQxrt i(1-69)5) 雷诺应力模型雷诺应力模型(RSM )是求解雷诺应力张量的各个分量的输运方程。具体形式为(puu ) +(pU uu ) = -puu u + p(5 u +5 u ) +Qti jQxk i jQxi j k kj i ik jkk(Q(QUQU )口一u u-pu u j+

47、u uiJ Qxi j 丿JikQxjkQx丿QQxk“QuQu )-pp (g u 0 + g u 0) +(1-70)i j j iQu Qui j - 2(u u + u u )Qx Qxk j m ikm i m jkmkkr + j 2uQxQx 丿ji式中：左边的第二项是对流项C，右边第一项是湍流扩散项D T，第二项是分子扩散项DL， ijijij第三项是应力产生项P，第四项是浮力产生项G，第五项是压力应变项，第六项是耗散 ijijij项，第七项系统旋转产生项F。ijij在式(1-69)中，C、Dl、P、F不需要模拟，而Dt、G、需要模拟以封闭ijijijijijijijij方程。

48、下面简单对几个需要模拟项进行模拟。D T可以用Delay和Harlow的梯度扩散模型来模拟，但这个模型会导致数值不稳定，ij在 Fluent 中是采用标量湍流扩散模型：DtijQxkI Qu u、tijQQxJ k k丿(1-71)k2式中：湍流粘性系数用卩=pc 来计算，根据Lien和Leschziner, q = 0.82，这和标准k - t卩k模型中选取1.0有所不同。压力应变项 可以分解为三项，即ij = + + w(1-72)ijij ,1ij ,2ij式中： 、和 w分别是慢速项、快速项和壁面反射项，具体表述可以参见文献2。ij ,1ij ,2ij浮力引起的产生项G模拟为ij耗散张

49、量 模拟为ijGj =咋rtQTQTg+ g亍i Qxj QxJ ji丿(1-73)2 = 8 (p + Y )(1-74)ij 3 ij M式中：Y = 2pM 2, M是马赫数；标量耗散率用标准k -模型中采用的耗散率输运方 Mtt程求解。6) 大涡模拟 湍流中包含了不同时间与长度尺度的涡旋。最大长度尺度通常为平均流动的特征长度尺 度。最小尺度为Komogrov尺度。LES的基本假设是：动量、能量、质量及其他标量主要 由大涡输运；流动的几何和边界条件决定了大涡的特性，而流动特性主要在大涡中体现； 小尺度涡旋受几何和边界条件影响较小，并且各向同性，大涡模拟(LES)过程中，直接求 解大涡，小

50、尺度涡旋模拟，从而使得网格要求比DNS低。LES的控制方程是对Navier-Stokes方程在波数空间或者物理空间进行过滤得到的。过 滤的过程是去掉比过滤宽度或者给定物理宽度小的涡旋，从而得到大涡旋的控制方程：即+ u泗=0QtQx(1-75)QiQQQ uQ pQt(pu u ) =(U l)jQxi jQxQxQxQxjjjjj(1-76)Qt(pu) +式中：T为亚网格应力，T = pu u - pu - u。ij ij i j i j很明显，上述方程与雷诺平均方程很相似，只不过大涡模拟中的变量是过滤过的量，而 非时间平均量，并且湍流应力也不同。1.2.3 初始条件和边界条件计算流体动力

51、学(CFD)分析中，初始条件和边界条件的正确设置是关键的一步。现有的 CFD 软件都提供了现成的各种类型的边界条件，这里对有关的初始条件和边界条件作一般 讨论。1. 初始条件顾名思义，初始条件就是计算初始给定的参数，即=t0时给出各未知量的函数分布，如u = u(x, y, z, t ) = u (x, y, z)0 0v = v(x, y, z, t ) = v (x, y, z)00w = w( x, y, z, t ) = w (x, y, z)00(1-77)p = p(x, y, z,t ) = p (x, y, z)00P = P(x, y, z, t ) = p (x, y, z

52、)0 0T = T(x, y, z, t ) = T (x, y, z)00很明显，当流体运动定常时，无初始条件问题。2. 边界条件所谓边界条件就是流体力学方程组在求解域的边界上，流体物理量应满足的条件。例如， 流体被固壁所限，流体将不应有穿过固壁的速度分量；在水面这个边界上，大气压强认为是 常数(一般在距离不大的范围内可如此)；在流体与外界无热传导的边界上，流体与边界之间 无温差，如此等。由于各种具体问题不同，边界条件提法千差万别，一般要保持恰当：保 持在物理上是正确的;要在数学上不多不少，刚好能用来确定积分微分方程中的积分常数， 而不是矛盾的或有随意性。通常流体边界分为流固交界面和流流(液

53、液、液气)交界面，下面分别讨论。1) 流固分界面边界条件飞机、船舶在空气及水中运动时的流固分界面，水在岸边及底部的流固分界面，均属这 一类。一般而言，流体在固体边界上的速度依流体有无粘性而定。对于粘性流体，流体将粘 附于固体表面(无滑移)，即v I = v I(1-78)FS式中：vI是流体速度；vI是固壁面相应点的速度。式(1-78)表明，在流固边界面上，流体 FS在一点的速度等于固体在该点的速度。对于无粘性流体，流体可沿界面滑移，即有速度的切 向分量，但不能离开界面，也就是流体的法向速度分量等于固体的法向速度分量，即v I = vI(1-79)n F S另外，也可视所给条件，给出无温差条件

54、：TI = TI(1-80)FS式中： TI 是流体温度， TI 是固壁面相应点的温度。FS2) 液液分界面边界条件密度不同的两种液体的分界面就属于这一类。一般而言，对分界面两侧的液体情况经常给出的条件是v = v , T = T , p = p1 2 1 2 1 2 对应力及传导热情况给出的条件是=严I*1 dn 1du I2 dn 2Q =哼1 dndT=k 2 dn(1-81)(1-82)(1-83)3) 液气分界面边界条件 液气分界面最典型的是水与大气的分界面，即自由面。由于自由面本身是运动和变形的 而且其形状常常也是一个需要求解的未知函数，因此就有一个自由面的运动学条件问题。设 自由

55、面方程为F(x,y,z,t) =0(1-84)并假定在自由面上的流体质点始终保持在自由面上，则流体质点在自由面上一点的法向速度，应该等于自由面本身在这一点的法向速度。经过一系列推导(参见文献2)，得到自由液 面运动学条件：dF+ v -V F = 0(1-85)如果要考虑液气边界上的表面张力，则在界面两侧，两种介质的压强差与表面张力有如下关系：这就是自由面上的动力学条件r 1 1)IR R丿 12 当不考虑表面张力时，有p p =b21p = pa 式中： p 为大气压强。a4) 无限远的条件(1-86)(1-87)流体力学中的很多问题，流体域是无限远的。例如，飞机在空中飞行时，流体是无界的。如果将坐标系取在运动物体上，这时无限远处的边界条件为(1-88)u=u，p=pgg其中下标g表示无穷远处的值。