第1章 CFD 基 础
计算流体动力学(computational fluid dynamics ，CFD)是流体力学的一个分支，它通过计算机模拟获得某种流体在特定条件下的有关信息，实现了用计算机代替试验装置完成“计算试验”，为工程技术人员提供了实际工况模拟仿真的操作平台，已广泛应用于航空航天、
热能动力、土木水利、汽车工程、铁道、船舶工业、化学工程、流体机械、环境工程等 领域。
本章介绍CFD 一些重要的基础知识，帮助读者熟悉CFD 的基本理论和基本概念，为计算时设置边界条件、对计算结果进行分析与整理提供参考。
1.1 流体力学的基本概念
1.1.1 流体的连续介质模型
流体质点(fluid particle)：几何尺寸同流动空间相比是极小量，又含有大量分子的微 元体。
连续介质(continuum/continuous medium)：质点连续地充满所占空间的流体或固体。 连续介质模型(continuum/continuous medium model)：把流体视为没有间隙地充满它所占据的整个空间的一种连续介质，且其所有的物理量都是空间坐标和时间的连续函数的一种假设模型：u =u (t ,x ,y ,z )。
1.1.2 流体的性质
1. 惯性
惯性(fluid inertia)指流体不受外力作用时，保持其原有运动状态的属性。惯性与质量有关，质量越大，惯性就越大。单位体积流体的质量称为密度(density)，以r 表示，单位为kg/m 3。对于均质流体，设其体积为V ，质量为m ，则其密度为
m
V
ρ= (1-1)
对于非均质流体，密度随点而异。若取包含某点在内的体积V ?，其中质量m ?，则该点密度需要用极限方式表示，即
0lim
V m
V
ρ?→?=? (1-2) 2. 压缩性
作用在流体上的压力变化可引起流体的体积变化或密度变化，这一现象称为流体的可压缩性。压缩性(compressibility)可用体积压缩率k 来量度
d /d /d d V V k p p
ρρ
=-
=
(1-3) 式中：p 为外部压强。
在研究流体流动过程中，若考虑到流体的压缩性，则称为可压缩流动，相应地称流体为可压缩流体，例如高速流动的气体。若不考虑流体的压缩性，则称为不可压缩流动，相应地称流体为不可压缩流体，如水、油等。
3. 粘性
粘性(viscosity)指在运动的状态下，流体所产生的抵抗剪切变形的性质。粘性大小由粘度来量度。流体的粘度是由流动流体的内聚力和分子的动量交换所引起的。粘度有动力粘度μ和运动粘度ν之分。动力粘度由牛顿内摩擦定律导出：
d d u
y
τμ= (1-4)
式中：τ为切应力，Pa ；μ为动力粘度，Pa ?s ；d /d u y 为流体的剪切变形速率。
运动粘度与动力粘度的关系为
μ
νρ
=
(1-5) 式中：ν为运动粘度，m 2/s 。
在研究流体流动过程中，考虑流体的粘性时，称为粘性流动，相应的流体称为粘性流体；当不考虑流体的粘性时，称为理想流体的流动，相应的流体称为理想流体。
根据流体是否满足牛顿内摩擦定律，将流体分为牛顿流体和非牛顿流体。牛顿流体严格满足牛顿内摩擦定律且μ保持为常数。非牛顿流体的切应力与速度梯度不成正比，一般又分为塑性流体、假塑性流体、胀塑性流体3种。
塑性流体，如牙膏等，它们有一个保持不产生剪切变形的初始应力0τ，只有克服了这个初始应力后，其切应力才与速度梯度成正比，即
0d d u
y ττμ=+ (1-6)
假塑性流体，如泥浆等，其切应力与速度梯度的关系是
d 1d n
u n y τμ??
=< ???
， (1-7)
胀塑性流体，如乳化液等，其切应力与速度梯度的关系是
d 1d n
u n y τμ??
=> ???
， (1-8)
1.1.3 流体力学中的力与压强
1. 质量力
与流体微团质量大小有关并且集中在微团质量中心的力称为质量力(body force)。在重力场中有重力mg ；直线运动时，有惯性力ma 。质量力是一个矢量，一般用单位质量所具
有的质量力来表示，其形式如下：
x y z f f i f j f k =++ (1-9)
式中：x f ，y f ，z f 为单位质量力在各轴上的投影。
2. 表面力大小与表面面积有关而且分布作用在流体表面上的力称为表面力(surface force)。表面力按其作用方向可以分为两种：一是沿表面内法线方向的压力，称为正压力；另一种是沿表面切向的摩擦力，称为切向力。
对于理想流体的流动，流体质点只受到正压力，没有切向力；对于粘性流体的流动，流体质点所受到的作用力既有正压力，也有切向力。
作用在静止流体上的表面力只有沿表面内法线方向的正压力。单位面积上所受到的表面力称为这一点处的静压强。静压强具有两个特征：①静压强的方向垂直指向作用面； ②流场内一点处静压强的大小与方向无关。
3. 表面张力
在液体表面，界面上液体间的相互作用力称为张力。在液体表面有自动收缩的趋势，收缩的液面存在相互作用的与该处液面相切的拉力，称为液体的表面张力(surface tension)。正是这种力的存在，引起弯曲液面内外出现压强差以及常见的毛细现象等。
试验表明，表面张力大小与液面的截线长度L 成正比，即
T L σ= (1-10)
式中：σ为表面张力系数，它表示液面上单位长度截线上的表面张力，其大小由物质种类决定，其单位为N/m 。
4. 绝对压强、相对压强及真空度
标准大气压的压强是101325Pa(760mm 汞柱)，通常用p atm 表示。若压强大于大气压，则以该压强为计算基准得到的压强称为相对压强(relative pressure)，也称为表压强，通常用p r 表示。若压强小于大气压，则压强低于大气压的值就称为真空度(vacuum)，通常用p v 表示。如以压强0Pa 为计算的基准，则这个压强就称为绝对压强(absolute pressure)，通常用p s 表示。这三者的关系如下：
r s atm p p p =- (1-11)
v atm s p p p =- (1-12)
在流体力学中，压强都用符号p 表示，但一般来说有一个约定：对于液体，压强用相对压强；对于气体，特别是马赫数大于0.1的流动，应视为可压缩流，压强用绝对压强。
压强的单位较多，一般用Pa ，也可用bar ，还可以用汞柱、水柱，这些单位换算如下：
1Pa=1N/m 2
1bar=105Pa
1p atm =760mmHg=10.33mH 2O=101325Pa
5. 静压、动压和总压
对于静止状态下的流体，只有静压强。对于流动状态的流体，有静压强(static pressure)、动压强(dynamic pressure)、测压管压强(manometric tube pressure)和总压强(total pressure)之
分。下面从伯努利(Bernoulli)方程(也有人称其为伯努里方程)中分析它们的意义。
伯努利方程阐述一条流线上流体质点的机械能守恒，对于理想流体的不可压缩流动其表达式如下：
2
2p v z H g g ρ++= (1-13) 式中：/p g ρ称为压强水头，也是压能项，为静压强；2/2v g 称为速度水头，也是动能项；z 称为位置水头，也是重力势能项，这三项之和就是流体质点的总的机械能；H 称为总的
水头高。
将式(1-13)两边同时乘以g ρ，则有
21
2p v gz gH ρρρ++= (1-14)
式中：p 称为静压强，简称静压；21
2
v ρ称为动压强，简称动压；gH ρ称为总压强，简称
总压。对于不考虑重力的流动，总压就是静压和动压之和。
1.1.4 流体运动的描述
1. 流体运动描述的方法
描述流体物理量有两种方法，一种是拉格朗日描述；一种是欧拉描述。
拉格朗日(Lagrange)描述也称随体描述，它着眼于流体质点，并将流体质点的物理量认为是随流体质点及时间变化的，即把流体质点的物理量表示为拉格朗日坐标及时间的函数。设拉格朗日坐标为(a ,b ,c )，以此坐标表示的流体质点的物理量，如矢径、速度、压强等等在任一时刻t 的值，便可以写为a 、b 、c 及t 的函数。
若以f 表示流体质点的某一物理量，其拉格朗日描述的数学表达式为
(,,,)f f a b c t = (1-15)
例如，设时刻t 流体质点的矢径即t 时刻流体质点的位置以r 表示，其拉格朗日描述为
(,,,)r r a b c t = (1-16) 同样，质点的速度的拉格朗日描述是
(,,,)v v a b c t = (1-17) 欧拉描述，也称空间描述，它着眼于空间点，认为流体的物理量随空间点及时间而变化，即把流体物理量表示为欧拉坐标及时间的函数。设欧拉坐标为(q 1,q 2,q 3)，用欧拉坐标表示的各空间点上的流体物理量如速度、压强等，在任一时刻t 的值，可写为q 1、q 2、q 3及t 的函数。从数学分析知道，当某时刻一个物理量在空间的分布一旦确定，该物理量在此空间形成一个场。因此，欧拉描述实际上描述了一个个物理量的场。
若以f 表示流体的一个物理量，其欧拉描述的数学表达式是(设空间坐标取用直角坐标)
(,,,)(,)f F x y z t F r t == (1-18) 如流体速度的欧拉描述是
(,,,)v v x y z t = (1-19)
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 )上，则应有
(,,,)(,,,)(,,,)x x a b c t y y a b c t z z a b c t =??
=??=?
(1-20) (,,,)(,,,)F x y z t f a b c t = (1-21)
事实上，将式(1-16)代入式(1-21)左端，即有
(,,,)[(,,,),(,,,),(,,,),]
(,,,)F x y z t F x a b c t y a b c t z a b c t t f a b c t == (1-22)
或者反解式(1-16)，得到
(,,,)(,,,)(,,,)a a x y z t b b x y z t c c x y z t =??
=??=?
(1-23) 将式(1-23)代入式(1-21)的右端，也应有
(,,,)[(,,,),(,,,),(,,,),]
(,,,)f a b c t f a x y z t b x y z t c x y z t t F x y z t == (1-24)
由此，可以通过拉格朗日描述推出欧拉描述，同样也可以由欧拉描述推出拉格朗日 描述。
3. 随体导数
流体质点物理量随时间的变化率称为随体导数(substantial derivative)，或物质导数、质点导数。
按拉格朗日描述，物理量f 表示为(,,,)f f a b c t =，f 的随体导数就是跟随质点(a ,b ,c )的物理量f 对时间t 的导数/f t ??。例如，速度(,,,)v a b c t 是矢径(,,,)r a b c t 对时间的偏导数，
(,,,)
(,,,)r a b c t v a b c t t
?=
? (1-25) 即随体导数就是偏导数。
按欧拉描述，物理量f 表示为(,,,)f F x y z t =，但/F t ??并不表示随体导数，它只表示物理量在空间点(,,,)x y z t 上的时间变化率。而随体导数必须跟随t 时刻位于(,,,)x y z t 空间点上的那个流体质点，其物理量f 的时间变化率。由于该流体质点是运动的，即x 、y 、z 是变的，若以a 、b 、c 表示该流体质点的拉格朗日坐标，则x 、y 、z 将依式(1-16)变化，从而f =F (x ,y ,z ,t )的变化依连锁法则处理。因此，物理量f =F (x ,y ,z ,t )的随体导数是
D (,,,)
D [(,,,),(,,,),(,,,),]
D ()F x y z t F x a b c t y a b c t z a b c t t t
F x F y F z F x t y t z t t
F F F F u v w x y z t
F
v F t
=???????=+++???????????=+++?????=??+
? (1-26)
式中：D/D t 表示随体导数。
从中可以看出，对于质点物理量的随体导数，欧拉描述与拉格朗日描述大不相同。前者是两者之和，而后者是直接的偏导数。
4. 定常流动与非定常流动
根据流体流动过程以及流动过程中的流体的物理参数是否与时间相关，可将流动分为定常流动(steady flow)与非定常流动(unsteady flow)。
定常流动：流体流动过程中各物理量均与时间无关，这种流动称为定常流动。
非定常流动：流体流动过程中某个或某些物理量与时间有关，则这种流动称为非定常流动。
5. 流线与迹线
常用流线和迹线来描述流体的流动。
迹线(track)：随着时间的变化，空间某一点处的流体质点在流动过程中所留下的痕迹称为迹线。在t =0时刻，位于空间坐标(a ,b ,c )处的流体质点，其迹线方程为
d (,,,)d d (,,,)d d (,,,)d x a b c t u t y a b c t v t z a b c t w t =??
=??=? (1-27) 式中：u 、v 、w 分别为流体质点速度的三个分量；x 、y 、z 为在t 时刻此流体质点的空间 位置。
流线(streamline)：在同一个时刻，由不同的无数多个流体质点组成的一条曲线，曲线上每一点处的切线与该质点处流体质点的运动方向平行。流场在某一时刻t 的流线方程为
d d d (,,,)(,,,)(,,,)x y z
u x y z t v x y z t w x y z t ==
(1-28) 对于定常流动，流线的形状不随时间变化，而且流体质点的迹线与流线重合。在实际流场中除驻点或奇点外，流线不能相交，不能突然转折。
6. 流量与净通量
流量(flux)：单位时间内流过某一控制面的流体体积称为该控制面的流量Q ，其单位为m 3/s 。若单位时间内流过的流体是以质量计算，则称为质量流量Q m ；不加说明时“流量”一词概指体积流量。在曲面控制面上有
d A
Q v n A =??? (1-29)
净通量(net flux)：在流场中取整个封闭曲面作为控制面A ，封闭曲面内的空间称为控制体。流体经一部分控制面流入控制体，同时也有流体经另一部分控制面从控制体中流出，此时流出的流体减去流入的流体，所得出的流量称为流过全部封闭控制面A 的净流量(或净通量)，通过式(1-30)计算：
d A
q v n A =??? (1-30)
对于不可压缩流体来说，流过任意封闭控制面的净通量等于0。 7. 有旋流动与有势流动
由速度分解定理，流体质点的运动可以分解为： (1) 随同其他质点的平动； (2) 自身的旋转运动；
(3) 自身的变形运动(拉伸变形和剪切变形)。
在流动过程中，若流体质点自身做无旋运动(irrotational flow)，则称流动是无旋的，也就是有势的，否则就称流动是有旋流动(rotational flow)。流体质点的旋度是一个矢量，通常用ω表示，其大小为
12ω???=???i j k x y z u v w (1-31)
若ω=0，则称流动为无旋流动，否则就是有旋流动。
ω与流体的流线或迹线形状无关；粘性流动一般为有旋流动；对于无旋流动，伯努利方程适用于流场中任意两点之间；无旋流动也称为有势流动(potential flow)，即存在一个势函数(,,,)x y z t ?，满足：
grad V ?= (1-32) 即
u v w x y z
???
???=
==
???，， (1-33) 8. 层流与湍流
流体的流动分为层流流动(laminar flow)和湍流流动(turbulent flow)。从试验的角度来看，层流流动就是流体层与层之间相互没有任何干扰，层与层之间既没有质量的传递也没有动量的传递；而湍流流动中层与层之间相互有干扰，而且干扰的力度还会随着流动而加大，层与层之间既有质量的传递又有动量的传递。
判断流动是层流还是湍流，是看其雷诺数是否超过临界雷诺数。雷诺数的定义如下：
VL
Re ν
=
(1-34) 式中：V 为截面的平均速度；L 为特征长度；ν为流体的运动粘度。
对于圆形管内流动，特征长度L 取圆管的直径d 。一般认为临界雷诺数为2320，即
vd
Re ν
=
(1-35)
当Re <2320时，管中是层流；当Re >2320时，管中是湍流。
对于异型管道内的流动，特征长度取水力直径d H ，则雷诺数的表达式为
H Vd
Re ν
= (1-36)
异型管道水力直径的定义如下：
H 4
A
d S
= (1-37) 式中：A 为过流断面的面积；S 为过流断面上流体与固体接触的周长。
临界雷诺数根据形状的不同而有所差别。根据试验几种异型管道的临界雷诺数如 表1-1所示。
表1-1 几种异型管道的临界雷诺数
Va
ν
3Va
ν
()V
D d ν
-
2070 1930 1000
对于平板的外部绕流，特征长度取沿流动方向的长度，其临界雷诺数为5×105～3×106。
1.2 CFD 基本模型
流体流动所遵循的物理定律，是建立流体运动基本方程组的依据。这些定律主要包括质量守恒、动量守恒、动量矩守恒、能量守恒、热力学第二定律，加上状态方程、本构方程。在实际计算时，还要考虑不同的流态，如层流与湍流。
1.2.1 基本控制方程
1. 系统与控制体
在流体力学中，系统是指某一确定流体质点集合的总体。系统以外的环境称为外界。分隔系统与外界的界面，称为系统的边界。系统通常是研究的对象，外界则用来区别于系统。系统将随系统内质点一起运动，系统内的质点始终包含在系统内，系统边界的形状和所围空间的大小可随运动而变化。系统与外界无质量交换，但可以有力的相互作用，及能
量(热和功)交换。
控制体是指在流体所在的空间中，以假想或真实流体边界包围，固定不动形状任意的空间体积。包围这个空间体积的边界面，称为控制面。控制体的形状与大小不变，并相对于某坐标系固定不动。控制体内的流体质点组成并非不变的。控制体既可通过控制面与外界有质量和能量交换，也可与控制体外的环境有力的相互作用。
2. 质量守恒方程(连续性方程)
在流场中，流体通过控制面A 1流入控制体，同时也会通过另一部分控制面A 2流出控制体，在这期间控制体内部的流体质量也会发生变化。按照质量守恒定律，流入的质量与流出的质量之差，应该等于控制体内部流体质量的增量，由此可导出流体流动连续性方程的积分形式为
d d d d 0V A x y z v n A t ρ?
+?=?????? (1-38) 式中：V 表示控制体，A 表示控制面。等式左边第一项表示控制体V 内部质量的增量；第
二项表示通过控制表面流入控制体的净通量。
根据数学中的奥-高公式，在直角坐标系下可将其化为微分形式：
()()()0u v w u v w t x y z
ρρρρ????+++=???? (1-39) 对于不可压缩均质流体，密度为常数，则有
0u v w x y z ???++=??? (1-40) 对于圆柱坐标系，其形式为
()()()0r r z v v v v t r r r z
θρρρρρθ????++++=???? (1-41) 对于不可压缩均质流体，密度为常数，则有
0r r z
v v v v r r r z θθ???+++=??? (1-42) 3. 动量守恒方程(运动方程)
动量守恒是流体运动时应遵循的另一个普遍定律，描述为：在一给定的流体系统，其动量的时间变化率等于作用于其上的外力总和，其数学表达式即为动量守恒方程，也称为运动方程，或N-S 方程，其微分形式表达如下：
d d d d d d yx xx zx bx
xy yy zy by
yz xz zz bz p p p u F t x y z p p p v F t x y z p p p w F t x y z ρρρρρρ????=+++??????????=+++?????
?????=+++
????? (1-43) 式中：bx F 、by F 、bz F 分别是单位质量流体上的质量力在三个方向上的分量；yx p 是流体内应力张量的分量。
动量守恒方程在实际应用中有许多表达形式，其中比较常见的有如下几种。 (1) 可压缩粘性流体的动量守恒方程
d 22d 3d 22d 3x y u p u u v w f t x x x x y z u v w u y y x z x z v p v u v w f t y y y x y z v w z z y ρρμμμρρμμ??????????????
=++-+++
???? ???????????????
??????????????+++?? ? ?????????????
??????????????????=++-+++
???? ???????????
??????????+ ???????d 22d 3z
u v x y x w p w u v w f t z z z x y z w u v w y x x z z z z μρρμμμ????
???????
???
??????++?????
??????????
?
???????????????=++-+++???? ????????????????????????????????+++? ? ???????????????????? (1-44) (2) 常粘性流体的动量守恒方程
2d grad grad(div )d 3
v F p v v t μ
ρ
ρμ=-++? (1-45) (3) 常密度常粘性流体的动量守恒方程
2d grad d v
F p v t
ρρμ=-+? (1-46) (4) 无粘性流体的动量守恒方程(欧拉方程)
d grad d v
F p t
ρρ=- (1-47) (5) 静力学方程
grad F p ρ= (1-48)
(6) 相对运动方程
在非惯性参考系中的相对运动方程是研究像大气、海洋及旋转系统中流体运动的所必须考虑的。由理论力学得知，绝对速度a v 为相对速度r v 及牵连速度e v 之和，即
a r e v v v =+ (1-49)
其中，e 0v v r Ω=+?，0v 为运动系中的平动速度，Ω是其转动角速度，r 为质点矢径。
而绝对加速度a a 为相对加速度r a 、牵连加速度e a 及科氏加速度c a 之和，即
a r e c a a a a =++ (1-50)
其中，0e d d ()d d v a r r t t
ΩΩΩ=
+?+??，c r 2a v Ω=?。 将绝对加速度代入运动方程，即得到流体的相对运动方程
r b c r d div 2d v
F P a v t
ρρΩ=+-- (1-51)
4. 能量守恒方程
将热力学第一定律应用于流体运动，把式(1-51)各项用有关的流体物理量表示出来，即是能量方程。如式(1-52)所示。
eff eff ()[()]()i j j j ij h j i i i T
E u E p k h J u S t x x x ρρτ'''??????++=-++????????∑ (1-52) 式中：2
2
i u p
E h ρ=-+；eff k 是有效热传导系数，eff t k k k =+，其中t k 是湍流热传导系数，
根据所使用的湍流模型来定义；j J '是组分j 的扩散流量；h S 包括了化学反应热以及其他用
户定义的体积热源项；方程右边的前3项分别描述了热传导、组分扩散和粘性耗散带来的能量输运。
1.2.2 湍流模型
湍流是自然界广泛存在的流动现象。大气、海洋环境的流动，飞行器和船舰的绕流，叶轮机械、化学反应器、核反应器中的流体运动都是湍流。湍流流动的核心特征是其在物理上近乎于无穷多的尺度和数学上强烈的非线性，这使得人们无论是通过理论分析、实验研究还是计算机模拟来彻底认识湍流都非常困难。回顾计算流体力学的发展，特别是活跃的20世纪80年代，不仅提出和发展了一大批高精度、高分辨率的计算格式，从主控方程看相当成功地解决了欧拉方程的数值模拟，可以说欧拉方程数值模拟方法的精度已接近于它有效使用范围的极限；同时还发展了一大批有效的网格生成技术及相应的软件，具体实现了工程计算所需要的复杂外形的计算网格；且随着计算机的发展，无论从计算时间还是从计算费用考虑，欧拉方程都已能适用于各种实践所需。在此基础上，20世纪80年代还进行了求解可压缩雷诺平均方程及其三维定态粘流流动的模拟。20世纪90年代又开始一个非定常粘流流场模拟的新局面，这里所说的粘流流场具有高雷诺数、非定常、不稳定、剧烈分离流动的特点，显然需要继续探求更高精度的计算方法和更实用可靠的网格生成技术。但更为重要的关键性的决策将是，研究湍流机理，建立相应的模式，并进行适当的模拟仍是解决湍流问题的重要途径。
1. 湍流模型分类
湍流流动模型很多，但大致可以归纳为以下3类。
第一类是湍流输运系数模型，即将速度脉动的二阶关联量表示成平均速度梯度与湍流粘性系数的乘积，用笛卡儿张量表示为
2
3j i i j t ij j
i u u u u k x x ρμρδ????''-=+- ? ????? (1-53) 模型的任务就是给出计算湍流粘性系数t μ的方法。根据建立模型所需要的微分方程的数目，可以分为零方程模型(代数方程模型)、单方程模型和双方程模型。
第二类是抛弃了湍流输运系数的概念，直接建立湍流应力和其他二阶关联量的输运 方程。
第三类是大涡模拟。前两类是以湍流的统计结构为基础，对所有涡旋进行统计平均。
大涡模拟把湍流分成大尺度湍流和小尺度湍流，通过求解三维经过修正的Navier-Stokes 方程(纳维-斯托克斯方程，简称N-S 方程)，得到大涡旋的运动特性，而对小涡旋运动还采用上述的模型。
实际求解中，选用什么模型要根据具体问题的特点来决定。选择的一般原则是精度要高，应用简单，节省计算时间，同时也具有通用性。
Fluent 提供的湍流模型包括：单方程(Spalart-Allmaras)模型、双方程模型(标准-k ε模型、重整化群-k ε模型、可实现-k ε模型)及雷诺应力模型和大涡模拟，如图1-1所示。
图1-1 湍流模型详解
2. 平均量输运方程
雷诺平均就是把Navier-Stokes 方程中的瞬时变量分解成平均量和脉动量两部分。对于速度，有
i i i u u u '=+ (1-54) 式中：i u 和i u '分别是平均速度和脉动速度(1,2,3i =)。
类似地，对于压力等其他标量，也有
φφφ'=+ (1-55)
式中：φ表示标量，如压力、能量、组分浓度等。
把上面的表达式代入瞬时的连续与动量方程，并取平均(去掉平均速度上的横线)，可以把连续与动量方程写成如下的笛卡儿坐标系下的张量形式：
()0i i
u t x ρρ??+=?? (1-56) d 2d 3j i i l ij i j i j
j i l j u u u u p u u t x x x x x x ρ
μδρ?????????
??
?''=++-+-?? ? ? ???????????????
(1-57) 上面两个方程称为雷诺平均的Navier-Stokes(RANS)方程。它们和瞬时Navier-Stokes 方
程有相同的形式，只是速度或其他求解变量变成了时间平均量。额外多出来的项i j u u ρ''-是雷诺应力，表示湍流的影响。
对于密度变化的流动过程，如燃烧问题，需要采用法夫雷(Favre)平均才可以求解。法夫雷平均就是除了压力和密度本身以外，所有变量都用密度加权平均。变量的密度加权平均定义如下：
/φρφρ= (1-58) 式中：符号~表示密度加权平均，对应于密度加权平均值的脉动值用φ''表示，有φφφ''=+。显然，这种脉动值的简单平均值不为零，但它的密度加权平均值等于零，即00φρφ''''≠=，。 为了求解方程(1-57)，必须模拟雷诺应力项以使方程封闭。通常的方法是应用Boussinesq 假设，认为雷诺应力与平均速度梯度成正比，表达式如下：
23j i i i j t t ij j i i u u u u u k x x x ρμρμδ???????''-=+-+ ? ? ???????? (1-59) Boussinesq 假设被用于单方程模型和-k ε双方程模型。这种近似方法好处是与求解湍
流粘性系数有关的计算时间比较少。例如，在Spalart-Allmaras 单方程模型中只多求解一个表示湍流粘性的输运方程；在-k ε双方程模型中只需多求解湍动能k 和耗散率ε两个方程，湍流粘性系数用湍动能k 和耗散率ε的函数来描述。Boussinesq 假设的不足之处是假设t μ是个各向同性标量，对于一些复杂流动，该条件并不是严格成立，所以具有其应用局限性。
另外的近似方法是求解雷诺应力各分量的输运方程。这也需要额外再求解一个标量方程，通常是耗散率ε方程。这就意味着对于二维湍流流动问题，需要多求解4个输运方程，而三维湍流问题需要多求解7个方程，需要较多的计算时间，要求更高的计算机内存。
在很多情况下基于Boussinesq 假设的模型很好用，而且计算量并不是很大。但是，如果湍流场各向异性很明显，如强旋流动以及应力取得的二次流等流动中，求解RSM 模型可以得到更好的结果。
3. 常用湍流模型简介
1) 单方程(Spalart-Allmaras)模型
单方程模型求解变量是v ，表征出了近壁(粘性影响)区域以外的湍流运动粘性系数。v 的输运方程为
2d 1()d v b v j j j v v v v G v C Y t x x x ρμρσ???????????
=+++-?? ??? ??????????????
(1-60) 式中：v G 是湍流粘性产生项；v Y 是由于壁面阻挡与粘性阻尼引起的湍流粘性的减少；v σ和2b C 是常数；v 是分子运动粘性系数。
湍流粘性系数1t v v f μρ=，其中，1v f 是粘性阻尼函数，定义为3133
1v v f C χχ=+，v
v
χ≡
。而湍流粘性产生项v G 模拟为1v b G C Sv ρ=，其中22
2
v v
S S f k d
≡+
，2111v v f f χχ=-+，1b C 和k
是常数，d
是计算点到壁面的距离；S 12j i ij i j u u x x Ω??
??=- ? ?????
。在Fluent 软件中，
考虑到平均应变率对湍流产生也起到很大作用，prod ||min(0,||||)ij ij ij S C S ΩΩ≡+-，其中，prod C =2.0
，||ij Ω
||ij S ≡，平均应变率12j i ij i j u
u S x x ????=+ ? ?????
。
在涡量超过应变率的计算区域计算出来的涡旋粘性系数变小。这适合涡流靠近涡旋中心的区域，那里只有“单纯”的旋转，湍流受到抑止。包含应变张量的影响更能体现旋转对湍流的影响。忽略了平均应变，估计的涡旋粘性系数产生项偏高。
湍流粘性系数减少项v Y 为2
1v w w v Y C f d ρ??= ???，其中，6
61/6
3
6
31w w w C f g g C ??+= ? ?+??，
62()w g r C r r =+-，22v
r Sk d
≡，1w C 、2w C 、3w C 是常数，在计算r 时用到的S 受平均应变
率的影响。
上面的模型常数在Fluent 软件中默认值为10.1335b C =，20.622b C =，2/3v σ=，
17.1v C =，2112/(1)/w b b v C C k C σ=++，20.3w C =，3 2.0w C =，0.41k =。
2) 标准-k ε模型
标准-k ε模型需要求解湍动能及其耗散率方程。湍动能输运方程是通过精确的方程推导得到的，但耗散率方程是通过物理推理，数学上模拟相似原形方程得到的。该模型假设流动为完全湍流，分子粘性的影响可以忽略。因此，标准-k ε模型只适合完全湍流的流动过程模拟。标准-k ε模型的湍动能k 和耗散率ε方程为如下形式：
d d t k b M i k i k k G G Y t x x μρμρεσ???
???=+++--?? ?????????
(1-61) 2
132d ()d t k b i i C G C G C t x x k k εεεεμεεεερμρσ??????=+++-?? ????????? (1-62) 式中：k G 表示由于平均速度梯度引起的湍动能产生，b G 表示由于浮力影响引起的湍动能产生；M Y 表示可压缩湍流脉动膨胀对总的耗散率的影响。湍流粘性系数2
t k C μμρε
=。
在Fluent 中，作为默认值常数，1C ε=1.44， 2C ε=1.92，3C ε = 0.09，湍动能k 与耗散
率ε的湍流普朗特数分别为k σ=1.0，εσ=1.3。
3) 重整化群-k ε模型
重整化群-k ε模型是对瞬时的Navier-Stokes 方程用重整化群的数学方法推导出来的模型。模型中的常数与标准-k ε模型不同，而且方程中也出现了新的函数或者项。其湍动能与耗散率方程与标准-k ε模型有相似的形式：
eff d ()d k k b M i i k k G G Y t x x ραμρε????=++--??????
(1-63) 2eff 132d ()()d k b i i C G C G C R t x x k k
εεεεεεεεραμρ????=++--?????? (1-64)
式中：k G 表示由于平均速度梯度引起的湍动能产生，b G 表示由于浮力影响引起的湍动能产生；M Y 表示可压缩湍流脉动膨胀对总的耗散率的影响，
这些参数与标准-k ε模型中相同。k α和εα分别是湍动能k 和耗散率ε的有效湍流普朗特数的倒数。湍流粘性系数计算公式
为2d d 1v v v Cv ??=--，其中，eff /v μμ=，100v C ≈。对于前面方程的积分，可以精确到有效雷诺数(涡旋尺度)对湍流输运的影响，这有助于处理低雷诺数和近壁流动问
题的模拟。对于高雷诺数，上面方程可以给出：2
t k C μμρε
=，0.0845C μ=。这个结果非常
有意思，和标准-k ε模型的半经验推导给出的常数0.09C μ=非常近似。在Fluent 中，如果是默认设置，用重整化群-k ε模型时是针对的高雷诺数流动问题。如果对低雷诺数问题进行数值模拟，必须进行相应的设置。
4) 可实现-k ε模型
可实现-k ε模型的湍动能及其耗散率输运方程为
d d t k b M i k i k k G G Y t x x μρμρεσ???
???=+++--?? ?????????
(1-65) 21213d d t b i i C S C C C G t x x k ε
εε
μεεε
ρμρερσ??????=++-+?? ?????????
(1-66) 式中：1max 0.43,5C ηη??
=??+??
，/Sk ηε=。 在上述方程中，k G 表示由于平均速度梯度引起的湍动能产生，b G 表示由于浮力影响
引起的湍动能产生；M Y 表示可压缩湍流脉动膨胀对总的耗散率的影响；2C 和1C ε是常数；
k σ和εσ分别是湍动能及其耗散率的湍流普朗特数。在Fluent 中，作为默认值常数，1C ε=1.44，2C =1.9，k σ=1.0，εσ=1.2。
该模型的湍流粘性系数与标准-k ε模型相同。不同的是，粘性系数中的C μ不是常数，
而是通过公式计算得到*
01s
C U K
A A με
=
+，其中， *ij ij U ΩΩ，2ij ij ij k k Ω=Ωεω-，
2ΩΩεω+ij ij ijk k =，ij Ω表示在角速度k ω
旋转参考系下的平均旋转张量率。模型常数
0 4.04A =，s
A φ=，1)3
φ=，式中i j j k k i
S S S W S =，ij ij S S S ≡，12j i ij i j u u S x x ??
??=+ ? ?????
。从这些式子中发现，C μ是平均应变率与旋度的函数。在平衡边界层惯性底层，可以得到0.09C μ=，与标准-k ε模型中采用的常数一样。
该模型适合的流动类型比较广泛，包括有旋均匀剪切流、自由流(射流和混合层)、腔道流动和边界层流动。对以上流动过程模拟结果都比标准-k ε模型的结果好，特别是可实现-k ε模型对圆口射流和平板射流模拟中，能给出较好的射流扩张角。
双方程模型中，无论是标准-k ε模型、重整化群-k ε模型还是可实现-k ε模型，三个模型有类似的形式，即都有k 和ε的输运方程，它们的区别在于：①计算湍流粘性的方法
不同；②控制湍流扩散的湍流普朗特数不同；③ε方程中的产生项和G k 关系不同。但都包含了相同的表示由于平均速度梯度引起的湍动能产生k G ，表示由于浮力影响引起的湍动能产生b G ；表示可压缩湍流脉动膨胀对总的耗散率的影响M Y 。
湍动能产生项
j k i j i
u G u u x ρ?''
=-? (1-67)
t b i
rt i
T
G g P x μβ?=? (1-68)
式中：rt P 是能量的湍流普特朗数，对于可实现-k ε模型，默认设置值为0.85；对于重整化群-k ε模型，1/rt P α=，1//rt p P k C αμ==。热膨胀系数1p
T ρβρ???
=- ????，对于理想气体，浮力引起的湍动能产生项变为
t b i
rt i
G g P x μρ
ρ?=-? (1-69)
5) 雷诺应力模型
雷诺应力模型(RSM)是求解雷诺应力张量的各个分量的输运方程。具体形式为
()()[()]()22()i j k i j i j k kj i ik j k k
j i i j i k j k i j j i k
k k k j j
i i k j m ikm i m jkm j
i k k u u U u u u u u p u u t x x U U u u u u u u g u g u x x x x u u u u p u u u u x x x x ρρρδδμρρβθθμρΩεε???+=-+++????????
???-+-++ ? ????????
???????+--+ ? ??????? (1-70) 式中：左边的第二项是对流项ij C ，右边第一项是湍流扩散项T ij D ，第二项是分子扩散项L ij D ，
第三项是应力产生项ij P ，第四项是浮力产生项ij G ，第五项是压力应变项ij φ，第六项是耗散项ij ε，第七项系统旋转产生项ij F 。
在式(1-69)中，ij C 、L ij D 、ij P 、ij F 不需要模拟，而T ij D 、ij G 、ij φ、ij ε需要模拟以封闭方程。下面简单对几个需要模拟项进行模拟。
T ij D 可以用Delay 和Harlow 的梯度扩散模型来模拟，但这个模型会导致数值不稳定，在Fluent 中是采用标量湍流扩散模型：
i j T t ij k
k k u u D x x μσ????=
? ?????
(1-71) 式中：湍流粘性系数用2
t k C μ
μρε
=来计算，根据Lien 和Leschziner ，0.82k σ=，这和标准
-k ε模型中选取1.0有所不同。
压力应变项ij φ可以分解为三项，即
,1,2w ij ij ij ij φφφφ=++ (1-72)
式中：,1ij φ、,2ij φ和w ij φ分别是慢速项、快速项和壁面反射项，具体表述可以参见文献[2]。
浮力引起的产生项ij G 模拟为
μβ?
?
??=+ ? ????
?t ij i j rt j i T T G g g P x x (1-73) 耗散张量ij ε模拟为
2
()3
ij ij M Y εδρε=+ (1-74)
式中：22M t Y M ρε=，t M 是马赫数；标量耗散率ε用标准-k ε模型中采用的耗散率输运方
程求解。
6) 大涡模拟
湍流中包含了不同时间与长度尺度的涡旋。最大长度尺度通常为平均流动的特征长度尺度。最小尺度为Komogrov 尺度。LES 的基本假设是：①动量、能量、质量及其他标量主要由大涡输运；②流动的几何和边界条件决定了大涡的特性，而流动特性主要在大涡中体现；③小尺度涡旋受几何和边界条件影响较小，并且各向同性，大涡模拟(LES)过程中，直接求解大涡，小尺度涡旋模拟，从而使得网格要求比DNS 低。
LES 的控制方程是对Navier-Stokes 方程在波数空间或者物理空间进行过滤得到的。过滤的过程是去掉比过滤宽度或者给定物理宽度小的涡旋，从而得到大涡旋的控制方程：
0i i
u u t x ρρ??+=?? (1-75) ()()()i ij i i j j j j j j
u p u u u t x x x x x τρρμ??????
+=--?????? (1-76) 式中：ij τ为亚网格应力，ij i j i j u u u u τρρ=-?。
很明显，上述方程与雷诺平均方程很相似，只不过大涡模拟中的变量是过滤过的量，
而非时间平均量，并且湍流应力也不同。
1.2.3 初始条件和边界条件
计算流体动力学(CFD)分析中，初始条件和边界条件的正确设置是关键的一步。现有的CFD 软件都提供了现成的各种类型的边界条件，这里对有关的初始条件和边界条件作一般讨论。
1. 初始条件
顾名思义，初始条件就是计算初始给定的参数，即0t t =时给出各未知量的函数分布，如 000000000000(,,,)(,,)(,,,)(,,)(,,,)(,,)
(,,,)(,,)(,,,)(,,)(,,,)(,,)
u u x y z t u x y z v v x y z t v x y z w w x y z t w x y z p p x y z t p x y z x y z t x y z T T x y z t T x y z ρρρ==??
==??==??
==?
?==?==?? (1-77) 很明显，当流体运动定常时，无初始条件问题。
2. 边界条件
所谓边界条件就是流体力学方程组在求解域的边界上，流体物理量应满足的条件。例如，流体被固壁所限，流体将不应有穿过固壁的速度分量；在水面这个边界上，大气压强认为是常数(一般在距离不大的范围内可如此)；在流体与外界无热传导的边界上，流体与边界之间无温差，如此等。由于各种具体问题不同，边界条件提法千差万别，一般要保持恰当：①保持在物理上是正确的；②要在数学上不多不少，刚好能用来确定积分微分方程中的积分常数，而不是矛盾的或有随意性。
通常流体边界分为流固交界面和流流(液液、液气)交界面，下面分别讨论。 1) 流固分界面边界条件
飞机、船舶在空气及水中运动时的流固分界面，水在岸边及底部的流固分界面，均属这一类。一般而言，流体在固体边界上的速度依流体有无粘性而定。对于粘性流体，流体将粘附于固体表面(无滑移)，即
F S ||v v = (1-78) 式中：F |v 是流体速度；S |v 是固壁面相应点的速度。式(1-78)表明，在流固边界面上，流体在一点的速度等于固体在该点的速度。对于无粘性流体，流体可沿界面滑移，即有速度的切向分量，但不能离开界面，也就是流体的法向速度分量等于固体的法向速度分量，即
n F S ||v v = (1-79)
另外，也可视所给条件，给出无温差条件：
F S ||T T = (1-80) 式中：F |T 是流体温度，S |T 是固壁面相应点的温度。
2) 液液分界面边界条件
密度不同的两种液体的分界面就属于这一类。一般而言，对分界面两侧的液体情况经常给出的条件是
121212,,v v T T p p === (1-81)
对应力及传导热情况给出的条件是
1
122||u u n n τμμ??==?? (1-82) 1122||T T Q k k n n
??==?? (1-83)
3) 液气分界面边界条件
液气分界面最典型的是水与大气的分界面，即自由面。由于自由面本身是运动和变形的，而且其形状常常也是一个需要求解的未知函数，因此就有一个自由面的运动学条件问题。设自由面方程为
(,,,)0F x y z t = (1-84) 并假定在自由面上的流体质点始终保持在自由面上，则流体质点在自由面上一点的法向速度，应该等于自由面本身在这一点的法向速度。经过一系列推导(参见文献[2])，得到自由液面运动学条件：
0F
v F t
?+??=? (1-85)
如果要考虑液气边界上的表面张力，则在界面两侧，两种介质的压强差与表面张力有如下关系：
211211p p R R σ??
-=+ ??? (1-86)
这就是自由面上的动力学条件。当不考虑表面张力时，有
a p p = (1-87) 式中：a p 为大气压强。
4) 无限远的条件
流体力学中的很多问题，流体域是无限远的。例如，飞机在空中飞行时，流体是无界的。如果将坐标系取在运动物体上，这时无限远处的边界条件为
当x →∞时，
u u p p ∞∞==， (1-88) 其中下标∞表示无穷远处的值。
1.3 CFD 模型的离散——有限体积法
1.3.1 CFD 模型的数值求解方法概述
从上面的分析看到，CFD 模型(控制方程)是一系列偏微分方程组，要得到解析解比较困难，目前，均采用数值方法得到其满足实际需要的近似解。
数值方法求解CFD 模型的基本思想是：把原来在空间与时间坐标中连续的物理量的场(如速度场、温度场、浓度场等)，用一系列有限个离散点(称为节点，node)上的值的集合来代替，通过一定的原则建立起这些离散点上变量值之间关系的代数方程(称为离散方程，discretization equation)，求解所建立起来的代数方程以获得所求解变量的近似解。在过去的几十年内已经发展了多种数值解法，其间的主要区别在于区域的离散方式、方程的离散方式及代数方程求解的方法这三个环节上。在CFD 求解计算中用得较多的数值方法有：有限差分法(finite difference method ，FDM)、有限体积法(finite volume method ，FVM)、有限元法(finite element method ，FEM)及有限分析法(finite analytic method ，FAM)。下面简要介绍，后面将着重介绍有限体积法。
1. 有限差分法
有限差分法是历史上采用最早的数值方法，对简单几何形状中的流动与换热问题也是一种最容易实施的数值方法。其基本点是：将求解区域用与坐标轴平行的一系列网格线的交点所组成的点的集合来代替，在每个节点上，将控制方程中每一个导数用相应的差分表达式来代替，从而在每个节点上形成一个代数方程，每个方程中包括了本节点及其附近一些节点上的未知值，求解这些代数方程就获得了所需的数值解。由于各阶导数的差分表达式可以从Taylor(泰勒)展开式来导出，这种方法又称建立离散方程的Taylor 展开法。
有限差分法软件一般研究者自己编写，很少看到商品的有限差分法软件。
2. 有限体积法
在有限体积法中将所计算的区域划分成一系列控制体积，每个控制体积都有一个节点作代表，通过将守恒型的控制方程对控制体积作积分来导出离散方程。在导出过程中，需要对界面上的被求函数本身及其一阶导数的构成作出假定，这种构成的方式就是有限体积法中的离散格式。用有限体积法导出的离散方程可以保证具有守恒特性，而且离散方程系数的物理意义明确，是目前流动与传热问题的数值计算中应用最广泛的一种方法。
Phoenics是最早投入市场的有限体积法软件，Fluent、STAR-CD和CFX都是常用的有限体积法软件，它们在流动、传热传质、燃烧和辐射等方面应用广泛。
3. 有限元法
在有限元法中把计算区域划分成一系列单元体(在二维情况下，单元体多为三角形或四边形)，在每个单元体上取数个点作为节点，然后通过对控制方程做积分来获得离散方程。它与有限体积法区别主要在于如下两点。
(1) 要选定一个形状函数(最简单的是线性函数)，并通过单元体中节点上的被求变量之值来表示该形状函数，在积分之前将该形状函数代入到控制方程中去。这一形状函数在建立离散方程及求解后结果的处理上都要应用。
(2) 控制方程在积分之前要乘上一个权函数，要求在整个计算区域上控制方程余量(即代入形状函数后使控制方程等号两端不相等的差值)的加权平均值等于零，从而得出一组关于节点上的被求变量的代数方程组。
有限元法的最大优点是对不规则区域的适应性好。但计算的工作量一般较有限体积法大，而且在求解流动与换热问题时，对流项的离散处理方法及不可压流体原始变量法求解方面没有有限体积法成熟。
Ansys、Sysweld和北京飞箭公司的FEPG(finite element programs generator)等有限元软件比较流行。
4. 有限分析法
有限分析法是由美籍华裔科学家陈景仁教授在1981年提出的。在这种方法中，也像有限差分法那样，用一系列网格线将区域离散，所不同的是每个节点与相邻的4个网格(二维)问题组成计算单元，即一个计算单元由一个中心节点与8个邻点组成。在计算单元中把控制方程中的非线性项(如Navier-Stokes方程中的对流项)局部线性化(即认为流速已知)，并对该单元上未知函数的变化型线作出假设，把所选定型线表达式中的系数和常数项用单元边界节点上未知的变量值来表示，这样该单元内的被求问题就转化为第一类边界条件下的一个定解问题，可以找出其分析解；然后利用这一分析解，得出该单元中点及边界上8个邻点上未知值间的代数方程，此即为单元中点的离散方程。有限分析法中的系数不像有限体积法中那样有明确的物理意义，对不规则区域的适应性也较差。
1.3.2 有限体积法
从上面的简介看到，有限体积法是一种分块近似的计算方法，其中比较重要步骤是计算区域的离散和控制方程的离散。