1.1　CFD基础

1.1.1　CFD的基本方程

CFD是计算流体力学的英文缩写，即利用数值解算方法求解流体力学的基本控制方程，包括连续性方程、动量方程和能量方程等。

一、流体力学的连续性方程

连续性方程即质量守恒方程，任何流动问题都必须满足质量守恒定律。按照质量守恒定律，单位时间内流出控制体的流体净质量总和应等于同时间间隔控制体内因密度变化而减少的质量，由此可导出流体流动连续性方程的微分形式为

式中，u_x、u_y、u_z分别为x、y、z三个方向的速度分量，m/s；t为时间，s；ρ为密度，kg/m³。

若用哈密顿微分算子表示，可以写成散度的形式

二、流体力学的动量方程

动量方程的本质是满足牛顿第二定律。该定律可描述为：对于一给定的流体微元体，其动量对时间的变化率等于外界作用在该微元体上的各种力之和。依据这一定律，可导出x、y和z三个方向的动量方程为

式中，p为流体微元体上的压强，Pa；τ_xx、τ_yy、τ_zz等是因分子黏性作用而产生的作用在微元体表面上的黏性应力τ的分量，Pa；f_x、f_y、f_z为三个方向的单位质量力，m/s²。

动量方程在实际应用中有许多表达形式，其中比较常见的有如下几种。

（1）可压缩黏性流体的动量方程。

（2）常黏性流体的动量方程。

（3）常密度黏性流体的动量守恒方程。

（4）无黏性流体的动量守恒方程。

（5）静力学方程。

实际流体的连续性方程与动量方程组成的方程组就称为纳维-斯托克斯（N-S）方程。

三、流体力学的能量方程

能量守恒定律是包含有热交换的流动系统必须满足的基本定律，其本质是热力学第一定律。依据能量守恒定律，微元体中能量的增加率等于进入微元体的净热流通量加上质量力与表面力对微元体所做的功，可得其表达式为

式中，E为流体微团的总能，J/kg，包含内能、动能和势能之和，；h为焓，J/kg；h_j为组分j的焓，J/kg，定义为，其中T_ref=298.15K；k_eff为有效热传导系数，W/(m·K)，，k_t为湍流热传导系数，根据所用的湍流模型来确定；J_j为组分j的扩散通量；S_h为包括了化学反应热及其他用户定义的体积热源项。

1.1.2　CFD的数值解法

有限体积法（FVM）又称有限容积法，是近年来发展迅速的一种离散化方法，其特点是计算效率高。有限体积法在CFD领域得到了广泛的应用，大多数商用CFD软件都采用了这种方法，如FLUENT、STAR-CD和CFX等。

有限体积法将所计算的区域划分成一系列控制体积，每个控制体积都有一个节点作代表，通过将控制方程对控制体积作积分来导出离散方程。在积分的过程中，需要对控制体积界面上的被求函数本身（对流通量）及其一阶导数（扩散通量）的构成做出假定，这就形成了不同的格式。由于扩散项多是采用相当于二阶精度的线性插值，因而格式的区别主要体现在对流通量项上。

用有限体积法导出的离散方程可以保证守恒特性，而且离散方程系数的物理意义明确，是目前流动与传热问题的数值计算中应用最普遍的一种方法。

由于流体控制方程较多，为研究方便，我们将控制方程写成通用变量方程的形式，即

式中的、Γ、的取值见表1-1。

表1-1　通用变量方程中的各参量取值

写成通用变量的形式，是为了下一步离散做准备的，这样只需离散一个通用变量方程即可。通用变量方程的意义为

采用有限体积法对方程进行离散，即将微分方程式（1-10）在控制容积内进行积分，则

利用奥式公式或高斯定律，式（1-12）化为

对于稳态问题，由于时间相关项等于零，式（1-13）变为

对于非稳态问题，还需在时间间隔内积分，则式（1-13）变为

下面以三维对流扩散问题为例，介绍一下有限体积法的离散过程。三维对流扩散问题的控制微分方程为

采用图1-1所示的离散网格系统，对方程式（1-16）在控制容积内积分，有

图1-1　三维问题的有限体积法结构化网格

针对内部节点（P点）研究其控制方程的离散形式，W、E分别为P节点左右侧邻近节点（相当于西、东两侧）；N、S分别为P节点前后侧邻近节点（相当于北、南两侧）；T、B分别为P节点上下侧邻近节点，由图1-1可知，对于结构化网格，控制容积的边界面积，，，于是由奥式公式，方程（1-17）可写成

采用线性插值（中央差分），有

源项也用线性化处理，有

令，，，，，，，，，，，，代入式（1-19）后，代入式（1-18），可得

整理节点变量，有

这时，需要在系数中引入，则式（1-21）变为

设，，，，，，，，那么（1-22）式可化为

式（1-23）就是三维对流扩散问题的离散方程格式，适用于三维问题的计算区域所有内部节点的离散方程构造。

有限体积法除了采用上述中心差分格式离散外，还有多种离散格式，如统一的离散方程为

式（1-24）中，ap取值取决于问题的阶数，对于一阶问题，，对于二阶问题，，其中a_w、a_ww、a_E、a_EE取决于所使用的离散格式，如表1-2所示。

表1-2　不同离散格式下系数a_w、a_ww、a_E、a_EE的计算公式

对于任何一种离散格式，读者都希望其既具有稳定性，又具有较高的精度，同时又能适应于不同的流动形式，但是实际上这种理想的离散格式是不存在的。表1-3所示为常见离散格式的性能粗略对比，方便读者实际计算时选用。

表1-3　常见离散格式的性能对比

流场计算的基本过程是在空间上用有限体积法或其他类似方法将计算区域离散成许多小的体积单元，在每个体积单元上对离散后的控制方程组进行求解。上面已经建立了与控制方程相应的离散方程，即代数方程组。但是，除了已知速度场求温度场这类简单的问题之外，所生成的离散方程不能直接用来求解，还必须对离散方程进行某种调整，并且对各未知量（压力、速度、温度等）的求解顺序及方式进行特殊处理。流场数值计算是针对常规解法存在的主要问题进行改善而形成的一系列方法集，其本质是对离散后的控制方程进行求解，可分为分离式解法和耦合式解法两大类，如图1-2所示。

图1-2　流场数值计算方法的分类

分离式解法不直接解联立方程组，而是顺序地、逐个地求解各变量代数方程组。

依据是否直接求解原始变量，分离式解法可分为非原始变量法和原始变量法。

涡量-速度法和涡量-流函数法是两种典型的非原始变量法。涡量-速度法不直接求解流场的原始变量p，而是求解旋度和速度u。涡量-流函数法不直接求解原始变量u和p，而是求解旋度和流函数。

常用的原始变量法有压力修正法、解压力泊松方程法和人为压缩法。解压力泊松方程法需要采用对方程取散度等方法将动量方程转变为泊松方程，然后求解泊松方程。人为压缩法主要是受可压气体可以通过联立求解速度分量和密度的方法来求解的启发，引入人为压缩性和人为状态方程，以此对不可压缩流体的连续性方程施加干扰。目前，工程上使用最广泛的流场计算方法是压力修正法，其实质是迭代法。压力修正法有很多实现方式，其中，压力耦合方程组的半隐式方法（SIMPLE算法）应用最为广泛，也是各商用CFD软件普遍采用的方法。

耦合式解法则同时求解离散化的控制方程组，联立求解各变量。耦合式解法可以分为所有变量全场联立求解（隐式解法）、部分变量全场联立求解（显隐式解法）、在局部地区（如一个单元上）对所有变量联立求解（显式解法）。总体而言，耦合式解法计算效率低、内存消耗大。

下面详细介绍几种数值计算常用的流场迭代求解方法。

一、SIMPLE算法

SIMPLE算法基于交错网格，所谓交错网格是指将速度分量和压力在不同的网格系统上离散，如压力、温度和密度等在正常的网格节点上存储和计算；而将速度的分量分别在错位后的网格上存储和计算，错位后网格的中心位于原控制体积的界面上。二维流动计算的交错网格如图1-3所示，有3套不同的网格系统。对于三维问题有4套网格系统，在此不再详述。

图1-3　二维流动计算的交错网格

SIMPLE算法由Patankar和Spalding于1972年提出，是一种主要用于求解不可压缩流场的数值方法（也可用于求解可压流动）。它的核心是采用“猜测-修正”的过程，在交错网格的基础上来计算压力场，从而达到求解动量方程的目的。其基本思路可描述如下：对于给定的压力场（它可以是假定的值，或是上一次迭代计算所得到的结果），求解离散形式的动量方程，得出速度场；因为压力场是假定的或不精确的，由此得到的速度场一般不满足连续性方程，因此，必须对给定的压力场加以修正；修正的原则是，与修正后的压力场相应的速度场能满足这一迭代层次上的连续性方程。据此原则，我们把由动量方程的离散形式所规定的压力与速度的关系代入连续方程的离散格式，从而得到压力修正方程，由压力修正方程得出压力修正值；接着，根据修正后的压力场，求得新的速度场；然后检查速度场是否收敛；若不收敛，用修正后的压力值作给定的压力场，开始下一层次的计算；如此反复，直到获得收敛的解。

在上述求解过程中，如何获得压力修正值（即如何构造压力修正方程），以及如何根据压力修正值确定“正确”的速度（即如何构造速度修正方程），是SIMPLE算法的两个关键问题。

（1）速度修正方程。

针对直角坐标系下的二维层流稳态问题进行介绍，设有初始的猜测压力场p*，则动量方程的离散方程可借助该压力场得以求解，从而求出相应的速度分量和。

动量方程的离散方程为

式中，A为控制体积的界面面积，对于二维问题，即为△x或△y；b为动量方程的源项。

将猜测压力场代入得

定义修正的压力值p’为正确压力场p与猜测的压力场p*之差，有

同样，定义速度修正值、，有

将正确的压力场p代入动量离散方程，得到正确的速度场。将方程（1-25）减去（12-6），并假定源项b不变，有

引入压力修正值和速度修正值，则式（1-28）化为

由上式，可以通过压力修正值求出速度修正值。式中表明，任一点上速度的修正值由两部分组成，一部分是与该速度在同一方向上的相邻两节点间压力修正值之差，这是产生速度修正值的直接动力；另一部分是由邻点速度的修正值所引起的，这又可以视为四周压力的修正值对所讨论位置上速度改进的间接影响。

为简化式（1-29）的求解过程，可引入如下近似处理，即略去方程中与速度修正值相关的和，该近似是SIMPLE算法的重要特征，于是有

其中

将式（1-30）所示的速度修正值代入式（1-27），有

对于和，存在类似的表达式为

其中

式（1-32）与式（1-33）对比表明，如果已知压力修正值，便可对猜测的速度场作出相应的速度修正，从而得到正确的速度场。

（2）压力修正方程。

在速度修正方程的推导中，只考虑了动量方程，其实，如前所述，速度场还受到连续性方程的约束。这里针对二维稳态问题进行分析，其连续性方程为

对图1-4所示的标量控制体积，连续性方程满足如下离散形式

图1-4　用于离散连续性方程的标量控制体积

将正确的速度值代入式（1-36），有

该式可简记为

其中

式（1-38）为连续性方程的离散方程，也是压力修正值的离散方程。方程中的源项是由于不正确的速度场所导致的“连续性”不平衡量。通过求解方程（1-38），可得到空间所有位置的压力修正值。

式（1-39）中的密度是标量控制体积界面上的密度值，需要通过插值得到。无论采用何种插值方法，对于交界面所属的两个控制体积，必须采用同样的密度值。

为了求解式（1-38），还必须对压力修正值的边界条件进行说明。因为压力修正方程是动量方程和连续性方程的派生物，不是基本方程，所以其边界条件也与动量方程的边界条件有关。在一般的流场计算中，动量方程的边界条件有两种：一种是已知边界上的压力；另一种是已知沿边界法向的速度分量。若为第一种边界，已知边界压力，可以令，则压力修正值应为零；若为第二种边界，已知边界上的法向速度，在网格设计时令控制体积的界面与边界一致，这样，控制体积界面上的速度就为已知条件了。

（3）SIMPLE算法的计算流程如图1-5所示。

图1-5　SIMPLE算法流程图

二、SIMPLEC算法

SIMPLEC是英文SIMPLE Consisent的缩写，意为协调一致的SIMPLE算法。它是SIMPLE的改进算法之一，由Van Doormal和Raithby提出。

在SIMPLE算法中，为求解的方便，略去了速度修正方程中的项，从而把速度的修正完全归结于压差项的直接作用。这一作法虽然并不影响收敛解的值，但加重了修正值的负担，使得整个速度场迭代收敛速度降低。实际上，当我们在略去时，犯了一个“不协调一致”的错误。为了能略去而同时又能使方程基本协调，试在方程（1-29）的等号两端同时减去，有

可以预期，与其邻点的修正值具有相同的数量级，因而略去所产生的影响远比在方程（1-29）中不计所产生的影响要小得多。于是有

其中

同理，有

其中

得出修正后的速度计算公式为

式（1-45）与式（1-32）形式上一致，但系数d的计算公式不一样，这里需要用式（142）和式（1-44）计算。

这就是SIMPLEC算法，它与SIMPLE算法的计算步骤相同，只是速度修正方程中的系数项d的计算公式有所区别。由于SIMPLEC算法没有忽略项，所以得到的压力修正值一般比较合适，因此可不对进行欠松弛处理。SIMPLEC算法的流程图如图1-6所示。

图1-6　SIMPLEC算法流程图

三、PISO算法

PISO算法是Pressure Implicit with Splitting of Operators的缩写，意为压力的隐式算子分割算法，由Issa于1986年提出，起初是针对非稳态可压流动的无迭代计算所建立的一种压力速度计算程序，后来在稳态问题的迭代计算中也较广泛地使用了该方法。

SIMPLE算法和SIMPLEC算法是两步算法，即一步预测和一步修正，而PISO算法增加了一个修正步，在完成了第一步修正后寻求第二次改进值，目的是使它们更好地同时满足动量方程和连续性方程。由于PISO算法使用了预测-修正-再修正三步，从而加快了单个迭代步的收敛速度，其三个步骤介绍如下。

（1）预测步。

利用猜测的压力场p^*，求解动量离散方程，得到速度分量、。

（2）第一修正步。

引入对SIMPLE算法的第一个修正步，给出一个速度场（，），使其满足连续性方程。此处的修正公式与SIMPLE算法中的式（1-30）完全一致，但考虑到PISO算法还有第二个修正步，因此使用不同的记法

这组公式用于定义修正后的速度和

将式（1-47）代入连续性方程（1-36），产生与式（1-38）具有相同系数与源项的压力修正方程。求解该方程，产生第一个压力修正值。一旦压力修正值已知，可通过方程（1-47）获得速度分量和。

（3）第二修正步。

为强化SIMPLE算法的计算，PISO算法要进行第二步修正。和的动量离散方程是

再次求解动量方程，得出两次修正的速度场

注意修正步中的求和项是用速度分量和来计算的。

将式（1-49）减去（1-48），得

其中

将（1-50）式代入连续性方程（1-36），得到二次压力修正方程

其中

求解方程（1-52），就可以得到二次压力修正值，从而得到二次修正的压力场

最后，求解方程（1-50），得到二次修正的速度场。

在瞬态问题的非迭代计算中，压力场与速度场（，）被认为是准确的。对于稳态流动的迭代计算，PISO算法的实施过程如图1-7所示。

图1-7　PISO算法流程图

PISO算法要两次求解压力修正方程，因此，它需要额外的存储空间来计算二次压力修正方程中的源项。尽管该方法涉及较多的计算，但对比发现，它的计算速度很快，总体效率比较高。FLUENT的用户手册推荐，对于瞬态问题，PISO算法有明显的优势；而对于稳态问题，可能选SIMPLE或SIMPLEC算法更合适。

SIMPLE算法是压力-速度耦合求解系列算法的基础。目前在各种CFD软件中均提供了这种算法。SIMPLE的各种改进算法主要是提高了计算的收敛性，从而缩短计算时间。

在SIMPLE算法中，压力修正值能够很好地满足速度修正的要求，但对压力修正不是很理想。SIMPLEC和PISO算法总体上计算量要高出SIMPLE算法30%左右，收敛速度较SIMPLE算法快30%～50%。但对于不同类型的问题而言，每种算法有各自的优势。一般而言，动量方程与标量方程如果不是耦合在一起的，用PISO算法能够较好地收敛，且体现出较高的计算效率；而在动量方程与标量方程耦合非常密切时，SIMPLEC算法的效果更好些。

1.1.3　CFD的求解流程

CFD的求解流程主要分为3个部分：前处理、解算过程以及后处理。每个部分包括一些具体的操作，如图1-8所示。

图1-8　CFD软件求解流程

1.1.4　CFD商业软件

CFD技术的发展依托于计算机技术与数值计算方法的发展，在高性能计算机硬件平台的支撑下，CFD通用软件包得以出现并迅速商业化，对CFD技术在工程应用中的进一步推广起到了巨大的促进作用。

CFD软件于20世纪70年代诞生于美国，但近10年才真正得到较广泛的应用。为了完成CFD计算，早期需要用户自己编写计算程序，但由于CFD的复杂性及计算机软硬件的多样性，使得用户各自的应用程序往往缺乏通用性，而CFD本身又有其鲜明的系统性与规律性，因此，比较适合于制成通用的商用软件。自1981年以来，出现了如PHOENICS、CFX、STAR-CD、FIDAP、FLUENT等多个商用CFD软件，且随着计算机技术的快速发展，这些商用软件在工程界正发挥着越来越大的作用。

一、PHOENICS

这是世界上第一个投放市场的CFD商业软件，可以算是CFD商用软件的鼻祖。这一软件中所采用的一些基本算法，如SIMPLE方法、混合格式等，正是由该软件创始人D.B.Spalding及其合作者S.V.Patankar等提出的，对以后开发的商业软件有较大的影响。近年来，PHOENICS软件在功能上与方法上做了较大的改进，包括纳入拼片式多网格及细密网格嵌入技术，同位网格及非结构化网格技术；在湍流模型方面开发了通用的零方程、低Reynolds k-ε模型、RNG k-ε模型等。应用这一软件可计算大量的实际工作问题，包括城市污染预测、叶轮中的流动、管道流动等。

二、CFX

该软件采用有限体积法、拼片式块结构化网络，在非正交曲线坐标（适体坐标）系上进行离散，变量的布置采用同位网格方式。对流项的离散格式包括一阶迎风、混合格式、QUICK、CONDIF、MUSCI及高阶迎风格式。压力与速度的耦合关系采用SIMPLE系列算法（SIMPLEC），代数方程求解的方法中包括线迭代、代数多重网格、ICCG、STONE强隐方法及块隐式（BIM）。软件可计算不可压缩及可压缩流动、耦合传热、多相流、化学反应、气体燃烧等问题。

CFX属于ANSYS软件大家庭的一员，一直跟随ANSYS的版本号升级，至今也已经发展到CFX 15.0。

三、STAR-CD

STAR-CD软件名称是英文Simulation of Turbulent Flow in Arbitrary Region的缩写，连字符后的CD是开发商Computational Dynamics Ltd的缩写。这是基于有限容积法的一个通用软件。在网格生成方面，采用非结构化网格，单元的形态可以有六面体、四面体、三角形截面的棱柱体、金字塔形的锥体及六种形状的其他多面体。应用这一软件可以计算稳态与非稳态流动、牛顿流体及非牛顿流体的流动、多孔介质中的流动、亚音速及超音速流动，并且这一软件在汽车工业中的应用十分广泛。

四、FIDAP

FIDAP是英文Fluid Dynamics Analysis Package的缩写，美国Fluid Dynamics International Inc.于1983年推出，是世界上第一个使用有限元法（FEM）的CFD软件。可以接受如I-DEAS、PATRAN、ANSYS和ICEMCFD等著名生成网格的软件所产生的网格。该软件可以计算可压缩及不可压缩流、层流与湍流、单相与两相流、牛顿流体及非牛顿流体的流动问题。

五、FLUENT

FLUENT软件由美国FLUENT公司于1983年推出，是继PHOENICS软件之后第二个投放市场的基于有限体积法的软件。它包含结构化及非结构化网格两个版本。在结构化网格版本中有适体坐标的前处理软件，同时也可以纳入I-DEAS、PATRAN、ANSYS和ICEMCFD等著名软件生成的网格。速度与压力耦合采用同位网格上的SIMPLEC算法。对流项差分格式纳入了一阶迎风、中心差分及QUICK等格式。软件能计算可压缩及不可压缩流动、含有粒子的蒸发、燃烧过程、多组分介质的化学反应过程等问题。

2006年，FLUENT软件被ANSYS公司收购，作为ANSYS旗下的流体分析计算软件之一，FLUENT版本编号沿用ANSYS主软件版本编号，即FLUENT软件从FLUENT 6.3直接跨越为FLUENT 12，至今又经历了FLUENT 13、FLUENT 14、FLUENT 15.0三次版本升级。在ANSYS大家庭下，FLUENT的功能不断完善，且与ANSYS其他软件模块之间建立了数据交换的桥梁，如通过Workbench平台实现流固耦合模拟等。因此，FLUENT逐渐超越了其他流体分析软件，成为当下工程应用领域使用最为广泛的商用CFD软件。