ANSYS FLUENT 16.0超级学习手册

内容提要

本书以有限体积分析法（又称为控制容积法）为基础，结合作者多年的使用和开发经验，通过丰富的工程实例详细介绍ANSYS FLUENT 16.0在各个专业领域的应用。

全书分为基础和实例两个部分，共16章。基础部分详细介绍了流体力学的相关理论基础知识和ANSYS FLUENT 16.0软件，包括FLUENT软件、前处理、后处理、常用的边界条件等内容；实例部分包括导热问题、流体流动与传热、自然对流与辐射换热、凝固和融化过程、多相流模型、离散相、组分传输与气体燃烧、动网格问题、多孔介质内部流动与换热、UDF基础应用和燃料电池问题等的数值模拟。本书每个实例都有详细的说明和操作步骤，读者只需按书中的方法和步骤进行软件操作，即可完成一个具体问题的数值模拟和分析，进而逐步学会ANSYS FLUENT 16.0软件的使用。本书光盘配有书中实例的几何模型以及实例的网格模型，方便读者查阅。

本书内容翔实，既可以作为动力、能源、水利、航空、冶金、海洋、环境、气象、流体工程等专业领域的工程技术人员参考用书，也可以作为高等院校相关专业高年级本科生、研究生的学习用书。

前言

流体的流动规律以三大守恒定律为基础，即质量守恒定律、动量守恒定律和能量守恒定律。这些定律由数学方程组来描述，但由于这些方程组都是非线性的，对于一些复杂问题，传统的求解方法很难得到分析解。另一方面，随着计算机技术的不断发展和进步，计算流体动力学（CFD）逐渐在流体力学研究领域崭露头角，它通过计算机数值计算和图像显示方法，在时间和空间上定量描述流场的数值解，从而达到研究物理问题的目的。它兼具理论性和实践性，成为继理论流体力学和实验流体力学之后的又一种重要研究手段。

CFD软件最早于20世纪70年代诞生于美国，但其较广泛的应用是近十几年的事。目前，它已成为解决各种流体流动与传热问题的强有力工具，在水利、航运、海洋、环境、流体机械与流体工程等各种技术学科都有广泛的应用。

FLUENT是国际上流行的商用CFD软件包，包含基于压力的分离求解器、基于压力的耦合求解器、基于密度的隐式求解器、基于密度的显式求解器。它具有丰富的物理模型、先进的数值方法和强大的前后处理功能，可对高超音速流场、传热与相变、化学反应与燃烧、多相流、旋转机械、动/变形网格、噪声、材料加工复杂激励等流动问题进行精确的模拟，具有较高的可信度。

本书以ANSYS FLUENT 16.0作为软件平台，详尽地讲解了FLUENT软件的使用，全书共16章，各章的主要内容安排如下。

第1章：流体力学与计算流体力学基础。首先介绍了流体力学的基础理论知识，在此基础上介绍了计算流体力学的相关知识等内容。

第2章：FLUENT软件介绍。讲解FLUENT软件的特点、FLUENT与ANSYS Workbench之间的关系以及在Workbench中使用FLUENT的方法等内容，并在此基础上讲解了FLUENT的基本操作，最后通过一个“水流计算”的简单实例介绍了FLUENT的数值模拟方法。

第3章：前处理方法。首先简要介绍了主流前处理软件Gambit、ANSYS ICEM CFD、TGrid、GridPro和Gridgen的功能及特点，接着重点介绍了Gambit和ANSYS ICEM CFD 16.0的基本功能、基本用法和应用实例。通过实例介绍了用Gambit进行网格划分的基本步骤，详细介绍了用ANSYS ICEM CFD划分三维结构化网格的方法。

第4章：后处理方法。介绍了3种对FLUENT结果文件进行后处理的途径：FLUENT内置后处理、Workbench CFD-Post通用后处理器及Tecplot后处理软件，详细介绍了运用这些途径进行可视化图形处理、渲染以及图表、曲线和报告的生成方法。

第5章：FLUENT中常用的边界条件。首先对FLUENT中提供的各种边界条件进行了分类，接着阐述了FLUENT中流动入口和出口边界的各种参数确定方法，重点介绍了FLUENT中若干种常用边界条件的使用条件及方法。

第6章：导热问题的数值模拟。介绍了导热的基础理论，即傅里叶定律，然后通过两个实例对导热问题进行具体的数值模拟分析，包括有内热源的导热问题以及钢球非稳态冷却过程的数值模拟。

第7章：流体流动与传热的数值模拟。首先介绍了流体的两种流动状态—层流和湍流，然后介绍了FLUENT中的湍流模型，包括Spalart-Allmaras模型、k-ε模型、k-ω模型等。最后通过5个实例对其流场和温度场进行数值模拟。

其中引射器内流场、圆柱绕流和二维离心泵内流场的数值模拟属于流体流动的数值模拟；扇形教室空调通风的数值模拟属于流体流动与换热的数值模拟；地埋管流固耦合换热的数值模拟属于强制对流与导热耦合的数值模拟。

第8章：自然对流与辐射换热的数值模拟。首先介绍了自然对流与辐射换热的理论知识，然后通过3个实例分别对自然对流与辐射换热进行数值模拟。

两相连方腔内自然对流换热的数值模拟，左侧高温壁面以自然对流的形式通过中间壁面向右侧壁面传热，通过数值模拟可准确预测其内部温度场、压力场和速度场。

烟道内烟气对流辐射换热的数值模拟，主要是烟气中的三原子气体、非对称结构的双原子气体等对壁面有辐射换热，通过数值模拟可准确预测其内部的温度场、速度场。

室内通风问题模拟的是在英国菲尔德FLUENT欧洲办事处接待区的通风问题，考虑了不同材质墙体的传热和辐射问题，同时加载了夏季的太阳模型，得到室内温度分布情况和墙面太阳热流分布。

第9章：凝固和融化过程的数值模拟。介绍了凝固融化模型的基础理论。然后通过一个实例对其进行数值模拟，通过数值模拟可清晰地看到融化过程固液相的变化，并计算出冰块融化所需要的时间。

第10章：多相流模型的数值模拟。首先介绍了多相流的基础知识，然后介绍了FLUENT中的3种多相流模型，最后通过5个实例进行数值模拟。

其中孔口出流、水中气泡的上升和储油罐液面问题属于VOF模型，气穴现象的数值模拟属于Mixture模型，水流对沙滩冲刷过程的数值模拟属于Eulerian模型。

第11章：离散相的数值模拟。首先介绍了离散相模型的基础知识，然后通过两个实例进行详细的数值模拟分析。

引射器离散相流场的数值模拟，是在第7.2节引射器内流场的基础上添加离散相模型，用于模拟其内部烟灰的流动特性；喷淋过程的数值模拟是利用离散相的喷雾模型，对喷淋过程进行数值模拟。

第12章：组分传输与气体燃烧的数值模拟。首先介绍了基础理论知识，然后通过3个实例进行数值模拟分析。

室内甲醛污染物浓度的数值模拟，利用数值模拟方法准确预测室内甲醛的浓度；焦炉煤气燃烧和预混气体化学反应的模拟，利用数值模拟方法对多组分气体燃烧进行模拟，得到其温度场、速度场和各组分的浓度场。

第13章：动网格问题的数值模拟。首先介绍了FLUENT动网格的基础理论知识，然后通过4个实例进行数值模拟，包括两车交会过程、运动物体强制对流换热、双叶轮旋转流场和单级轴流涡轮机模型内部流场的数值模拟。

第14章：多孔介质内流动与换热的数值模拟。首先介绍了多孔介质的基础理论知识，然后介绍了FLUENT多孔介质模型，最后通过3个实例进行数值模拟分析。

第15章：UDF基础应用。首先介绍了UDF（用户自定义函数）的基本用法，然后用3个实例演示了UDF在定义物性参数、求解多孔介质和定义运动参数等方面的应用。

第16章：燃料电池问题模拟。主要向读者介绍如何使用FLUENT中的燃料电池附件模块来求解单通道逆流聚合物电解质膜（PEM）燃料电池问题。

通过本书的学习，读者可以在较短时间内把握FLUENT 16.0软件的学习要领，掌握FLUENT 16.0的详细操作步骤。各章所用到的实例均可从本书的配套光盘中找到。

本书以FLUENT 16.0版本为基础，其操作界面与老版本有较大不同，因此对新版本的操作界面进行了详细的说明，使读者能较快地掌握新版本的特点。

本书内容丰富、结构清晰，所有实例均经过精心设计与筛选，代表性强，并且每个实例都通过用户图形交互界面进行全过程操作。编写本书的主要目的不是求解多么复杂的物理问题，而是让读者学习FLUENT软件的求解思路和数值模拟软件的求解方法，强调软件的实用性，比如导热问题的数值模拟，其求解过程并不复杂，以往的书籍很少有涉及此问题的数值模拟，但实际工程中却有广泛的应用。

本书紧跟ANSYS软件发展的最前沿，对目前最新版FLUENT 16.0软件的部分新功能进行了详细的介绍与案例分析，希望对渴望入门的读者有所帮助。

本书由唐家鹏编著，虽然作者在本书的编写过程中力求叙述准确、完善，但由于水平有限，书中欠妥之处在所难免，希望读者能够及时指出，共同促进本书质量的提高。

读者在学习过程中遇到与本书有关的问题，可以发邮件到邮箱book_hai@126.com，或者访问博客http://blog.sina.com.cn/tecbook，作者会尽快给予解答。

第1章　流体力学与计算流体力学基础

流体力学是力学的一个重要分支，它主要研究流体本身的静止状态和运动状态，以及流体和固体界壁间有相对运动时的相互作用和流动的规律，在生活、环保、科学技术及工程中具有重要的应用价值。

计算流体力学或计算流体动力学（Computational Fluid Dynamics，CFD），是用电子计算机和离散化的数值方法对流体力学问题进行数值模拟和分析的一个分支。

本章先介绍流体力学中支配流体流动的基本物理定律，然后在此基础上介绍用数值方法求解流体力学问题的基本思想，进而阐述计算流体力学的相关基础知识，最后简要介绍常用的计算流体力学商业软件。

学习目标：

• 学习流体力学的基础知识，包括基本概念和重要理论；
• 学习计算流体力学的相关理论和方法；
• 了解CFD软件的构成；
• 了解常用的商业CFD软件。

1.1　流体力学基础

流体力学是连续介质力学的一个分支，是研究流体（包含气体及液体）现象以及相关力学行为的科学。

1.1.1　流体力学概述

1738年，伯努利在他的专著中首次采用了水动力学这个名词并作为书名；1880年前后出现了空气动力学这个名词；1935年以后，人们概括了这两方面的知识，建立了统一的体系，统称为流体力学。

在人们的生活和生产活动中随时随地都可遇到流体，因此流体力学是与人类日常生活和生产密切相关的。大气和水是最常见的两种流体，大气包围着整个地球，地球表面的70%是水面。大气运动、海水运动（包括波浪、潮汐、中尺度涡旋、环流等）乃至地球深处熔浆的流动都是流体力学的研究内容。

20世纪初，世界上第一架飞机出现以后，飞机和其他各种飞行器得到迅速发展。20世纪50年代开始的航天飞行，使人类的活动范围扩展到其他星球和银河系。航空航天事业的蓬勃发展是同流体力学的分支学科—空气动力学和气体动力学的发展紧密相联的。这些学科是流体力学中最活跃、最富有成果的领域。

石油和天然气的开采、地下水的开发利用，要求人们了解流体在多孔或缝隙介质中的运动，这是流体力学分支之一—渗流力学研究的主要对象。渗流力学还涉及土壤盐碱化的防治，化工中的浓缩、分离和多孔过滤，燃烧室的冷却等技术问题。

燃烧离不开气体，这是有化学反应和热能变化的流体力学问题，是物理化学流体动力学的内容之一。爆炸是猛烈的瞬间能量变化和传递过程，涉及气体动力学，从而形成了爆炸力学。

沙漠迁移、河流泥沙运动、管道中的煤粉输送、化工中气体催化剂的运动等，都涉及流体中带有固体颗粒或液体中带有气泡等问题，这类问题是多相流体力学研究的范围。

等离子体是自由电子、带等量正电荷的离子以及中性粒子的集合体。等离子体在磁场作用下有特殊的运动规律。研究等离子体运动规律的学科称为等离子体动力学和电磁流体力学，它们在受控热核反应、磁流体发电、宇宙气体运动等方面有广泛的应用。

风对建筑物、桥梁、电缆等的作用使它们承受载荷和激发振动，废气和废水的排放造成环境污染，河床冲刷迁移和海岸遭受侵蚀，研究这些流体本身的运动及其同人类、动植物间的相互作用的学科称为环境流体力学（其中包括环境空气动力学、建筑空气动力学）。这是一门涉及经典流体力学、气象学、海洋学和水力学、结构动力学等学科的新兴边缘学科。

生物流变学研究人体或其他动植物中有关的流体力学问题。例如血液在血管中的流动，心、肺、肾中的生理流体运动和植物中营养液的输送。此外，还研究鸟类在空中的飞翔，动物在水中的游动等。

因此，流体力学既包含自然科学的基础理论，又涉及工程技术科学方面的应用。

目前，研究流体力学问题的方法有理论分析研究、实验模拟研究和数值模拟方法研究3种。

流体力学理论分析的一般过程是：建立力学模型，用物理学基本定律推导流体力学数学方程，用数学方法求解方程，然后检验和解释求解结果。理论分析结果能揭示流动的内在规律，物理概念清晰，物理规律能公式化，具有普遍适用性，但分析范围有限，只能分析简单的流动。而且，线性问题能得到结果，非线性问题分析非常困难。

实验研究的一般过程是：在相似理论的指导下建立模拟实验系统，用流体测量技术测量流动参数，处理和分析实验数据。

典型的流体力学实验有风洞实验、水洞实验、水池实验等。测量技术有热线、激光测速，粒子图像、迹线测速，高速摄影，全息照相，压力密度测量等。现代测量技术在计算机、光学和图像技术配合下，在提高空间分辨率和实时测量方面已取得长足进步。

实验结果能反映工程中的实际流动规律，发现新现象，检验理论结果等，现象直观，测试结果可靠。但流体的实验研究对测试设备要求较高，设计制造周期长，且调试复杂。实验研究的方法只能得到有限的实验数据，真实模拟物理问题比较困难。

数值研究的一般过程是：对流体力学数学方程进行简化和数值离散化，编制程序进行数值计算，将计算结果与实验结果比较。

常用的数值模拟方法有有限差分法、有限元法、有限体积法、边界元法、谱分析法等。计算的内容包括飞机、汽车、河道、桥梁、涡轮机等流场的计算，湍流、流动稳定性、非线性流动等的数值模拟。大型工程计算软件已成为研究工程流动问题的有力武器。数值模拟方法的优点是能计算理论分析方法无法求解的数学方程（适用于线性和非线性问题），能处理各种复杂流动问题，比实验方法省时省钱。但毕竟是一种近似解方法，适用范围受数学模型的正确性和计算机的性能所限制。

流体力学的3种研究方法各有优缺点，在实际研究流体力学问题时，应结合实际问题，取长补短，互为补充和印证。

1.1.2　连续介质模型

如固体一样，流体也是由大量的分子组成的，而分子间都存在比分子本身尺度大得多的间隙，同时，由于每个分子都在不停地运动，因此，从微观的角度看，流体的物理量在空间分布上是不连续的，且随时间不断变化。

在流体力学中仅限于研究流体的宏观运动，其特征尺寸（如日常见到的是m、cm、mm那样的量级）比分子自由程大得多。描述宏观运动的物理参数，是大量分子的统计平均值，而不是个别分子的值。在这种情形下，流体可近似用连续介质模型处理。

连续介质模型认为，物质连续地分布于其所占有的整个空间，物质宏观运动的物理参数是空间及时间的可微连续函数。

根据连续介质模型假设，可以把流体介质的一切物理属性，如密度、速度、压强等都看作是空间的连续函数。因而，对于连续介质模型，微积分等现代数学工具可以加以应用。

连续介质模型假设成立的条件是建立在流体平均自由程远远小于物体特征尺寸的基础上的，即

（1-1）

在式1-1中，L为求解问题中物体或空间的特征尺寸；为流体分子的平均自由程。

在某些情况下，例如，在120 km的高空，如果空气分子的平均自由程和飞行器的特征尺寸在同一数量级，连续介质模型假设就不再成立。这时，必须把空气看成是不连续的介质，这个范围属于稀薄空气动力学范畴。

1.1.3　流体的基本概念及性质

1．密度

流体的密度定义为单位体积所含物质的多少，以ρ表示。密度是流体的一种固有物理属性，国际单位为。对于均质流体，设其体积为V，质量为m，则其密度为

（1-2）

对于非均质流体，不同位置的密度不同。若取包含某点的流体微团，其体积为，质量为，则该点的密度定义为

（1-3）

液体和气体的密度可以查询相关手册，这里为了方便读者，表 1-1 给出部分常见液体和气体的密度。

表1-1　常见气体和液体的密度（25℃常压情况下）

2．质量力和表面力

作用在流通微团上的力可分为质量力与表面力。

与流体微团质量大小有关并且集中作用在微团质量中心上的力称为质量力，如在重力场中的重力mg、直线运动的惯性力ma等。质量力是一个矢量，一般用单位质量所具有的质量力来表示，其形式为

（1-4）

式1-4中，、、为单位质量力在x、y、z轴上的投影，或简称为单位质量分力。

大小与表面面积有关而且分布作用在流体表面上的力称为表面力。表面力按其作用方向可分为两种：一种是沿表面内法线方向的压力，称为正压力；另一种是沿表面切向的摩擦力，称为切向力。

作用在静止流体上的表面力只有表面内法线方向的正压力，单位面积上所受到的表面力称为这一点处的静压强。静压强有两个特征：

• 静压强的方向垂直指向作用面。
• 流场内一点处静压强的大小与方向无关。

对于理想流体流动，流体质点只受到正压力，没有切向力。

对于黏性流体流动，流体质点所受到的表面力既有正压力，也有切向力。单位面积上所受到的切向力称为切应力。对于一元流动，切应力由牛顿内摩擦定律给出；对于多元流动，切应力可由广义牛顿内摩擦定律求得。

3．绝对压强、相对压强和真空度

单位面积上受到的压力叫作压强。在流体静力学中压强的定义是：作用在浸没于流体中物体表面上单位面积上的法向正应力。

压强在国际单位制中的单位是N/m²或者Pa（Pascal或帕斯卡）。另外压强也常用液柱高（汞柱、水柱等）、标准大气压（atm）和bar等单位进行度量，常用转换关系如下：

（1-5）

一个标准大气压的压强是760 mmHg，相当于101 325 Pa，通常用p_atm表示。

若压强大于标准大气压，则以标准大气压为计算基准得到的压强称为相对压强，也称为表压强（Gauge Pressure），通常用表示。

在FLUENT软件中求解器运算过程中实际使用的压强值都是表压强。FLUENT中默认的参考压力为一个标准大气压，用户可以指定参考压力。

若压强小于标准大气压，则压强低于大气压的值就称为真空度，通常用表示。

绝对压强、相对压强和真空度之间的关系为

（1-6）

所以在FLUENT中，如果用户将参考压力设置为0，表压强就等于绝对压强，有时候这么做可以方便边界条件的设置。

在流体力学中有如下约定：对于液体来说，压强用相对压强；对于气体，特别是马赫数大于0.3的流动，压强用绝对压强。

4．静压、动压和总压

物体在流体中运动时，在正对流体运动的方向的表面，流体完全受阻，此处的流体速度为0，其动能转变为压力能，压力增大，其压力称为全受阻压力（简称全压或总压），它与未受扰动处的压力（即静压）之差称为动压，即

（1-7）

式1-7中，为总压，p为静压，为动压，ρ为流体密度，为流体速度。也可以表达为：对于不考虑重力的流动，总压就是静压和动压之和。

根据伯努利方程的物理意义可知，对于一条理想流体，在一条流线上流体质点的机械能是守恒的，不可压缩流动的表达式为

（1-8）

式1-8中，称为压强水头，也是压能项，p为静压；为速度水头，也是动能项；第三项为位置水头，也是重力势能项。这三项之和就是流体质点的总机械能，等式右边的称为总水头。

将式1-8的左右两边同时乘以，则有

（1-9）

5．黏性

黏性是施加于流体的应力和由此产生的变形速率以一定的关系联系起来的流体的一种宏观属性，表现为流体的内摩擦。由于黏性的耗能作用，在无外界能量补充的情况下，运动的流体将逐渐停止下来。

黏性对物体表面附近的流体运动产生重要作用，使流速逐层减小并在物面上为零，在一定条件下也可使流体脱离物体表面。

黏性又称黏性系数、动力黏度，记为μ。牛顿黏性定律指出，在纯剪切流动中，相邻两流体层之间的剪应力（或黏性摩擦应力）为

（1-10）

式中，为剪切应力，为垂直流动方向的法向速度梯度。黏性数值上等于单位速度梯度下流体所受的剪应力。

速度梯度也表示流体运动中的角变形率，故黏性也表示剪应力与角变形率之间比值关系。按国际单位制，黏性的单位为kg/(m·s)。

黏性系数与密度之比称为运动黏性系数，常记作，则有

（1-11）

液体的黏性系数随温度的增加而减小，气体的黏性系数随温度的增加而变大。对于气体而言，黏性系数与温度的关系可以用萨瑟兰公式表示：

（1-12）

式1-12中，和分别为参考黏性系数和参考温度。

这里为了方便读者，表1-2列出了部分常见液体和气体的动力黏度和运动黏度。

表1-2 常见气体和液体的动力黏度和运动黏度（25℃常压情况下）

6．传热性

当气体中沿某一方向存在温度梯度时，热量就会由温度高的地方传向温度低的地方，这种性质称为气体的传热性。实验证明，单位时间内所传递的热量与传热面积成正比，与沿热流方向的温度梯度成正比，即

（1-13）

式中，q为单位时间通过单位面积的热量，负号表示热流量传递的方向永远和温度梯度的方向相反。

流体的导热系数λ随流体介质而不同，同一种流体介质的导热系数随温度变化而略有差异。

7．扩散性

当流体混合物中存在组元的浓度差时，浓度高的地方将向浓度低的地方输送该组元的物质，这种现象称为扩散。

流体的宏观性质，如扩散、黏性和热传导等，是分子输运性质的统计平均。由于分子的不规则运动，在各层流体间交换着质量、动量和能量，使不同流体层内的平均物理量均匀化。这种性质称为分子运动的输运性质。质量输运在宏观上表现为扩散现象，动量输运表现为黏性现象，能量输运则表现为热传导现象。

理想流体忽略了黏性，即忽略了分子运动的动量输运性质，因此，在理想流体中也不应考虑质量和能量输运性质—扩散和热传导，它们具有相同的微观机制。

8．流线和迹线

所谓流线，就是这样一种曲线，在某时刻，曲线上任意一点的切线方向正好与那一时刻该处的流速方向重合。可见，流线是由同一时刻的不同流点组成的曲线，它给出了该时刻不同流体质点的速度方向，是速度场的几何表示。

所谓迹线，就是流体质点在各时刻所行路径的轨迹线（或流体质点在空间运动时所描绘出来的曲线），如喷气式飞机飞过后留下的尾迹、台风的路径、纸船在小河中行走的路径等。

流线具有以下性质。

（1）同一时刻的不同流线，不能相交。

（2）流线不能是折线，而是一条光滑的曲线。

（3）流线族的疏密反映了速度的大小。

（4）实际流场中，除驻点或奇点外，流线不能突然转折。

流线的微分方程为

（1-14）

式1-14中，、和分别表示点在t时刻的速度在x、y和z方向上的分量。

迹线的微分方程为

（1-15）

式1-15中，t为关于时间的自变量，、、为t时刻此流体质点的空间位置。u、v和w分别为流体质点速度在x、y和z方向上的分量。

可见，对于定常流动，流线的形状不随时间变化，而且流体质点的迹线与流线重合。

9．流量和净通量

（1）流量。单位时间内流过某一控制面的流体体积称为该控制面的流量Q，其单位为m³/s。若单位时间内流过的流体是以质量计算的，则称为质量流体Q _m，不加说明时，“流量”一词概指体积流量。在曲面控制面上有

（1-16）

（2）净流量。在流场中取整个封闭曲面作为控制面A，封闭曲面内的空间称为控制体。流体经一部分控制面流入控制体，同时也有流体经另一部分控制面从控制体中流出。此时流出的流体减去流入的流体，所得出的流量称为流过全部封闭控制面A的净流量（或净通量），计算式为

（1-17）

对于不可压流体来说，流过任意封闭控制面的净通量等于0。

10．流速和音速

当把流体视为可压缩流体时，扰动波在流体中的传播速度是一个特征值，称为音速。音速方程式的微分形式为

（1-18）

音速在气体中的传播过程是一个等熵过程。将等熵方程式代入式1-18，并由理想气体状态方程，得到音速方程为

（1-19）

对于空气来说，k=1.4，R=287 J/(kg·K)，得到空气中的音速为

（1-20）

流速是流体流动的速度，而音速是扰动波的传播速度，两者之间的关系为

（1-21）

式1-21中，Ma称为马赫数。

11．马赫数和马赫锥

（1）马赫数。流体流动速度v与当地音速c之比称为马赫数，用Ma表示为

（1-22）

Ma<1的流动称为亚音速流动，Ma>1的流动称为超音速流动，Ma>3的流动称为高超音速流动。

（2）马赫锥。对于超音速流动，扰动波传播范围只能允许充满在一个锥形的空间内，这就是马赫锥，其半锥角θ称为马赫角，计算式如下：

（1-23）

12．滞止参数和等熵关系

（1）滞止参数。流场中速度为0的点上的各物理量称为滞止参数，用下标“0”表示，如滞止温度、滞止压强p ₀、滞止密度等。

（2）等熵流动基本关系式。流动参数与滞止参数及马赫数之间的基本关系为

（1-24）

（1-25）

（1-26）

式1-24的使用条件是绝热流动，式1-25和式1-26的使用条件是等熵流动。

13．正激波和斜激波

气流发生参数发生显著、突跃变化的地方，称为激波。激波常在超音速气流的特定条件下产生；激波的厚度非常小，约为10⁻⁴mm，因此一般不对激波内部的情况进行研究，所关心的是气流经过激波前后参数的变化。气流经过激波时受到激烈的压缩，其压缩过程是很迅速的，可以看作是绝热的压缩过程。

（1）正激波。

激波面与气流方向垂直，气流经过激波后方向不变，这称为正激波。

假设激波固定不动，激波前的气流速度、压强、温度和密度分别为v ₁、p ₁、T ₁和ρ ₁，经过激波后突跃地增加到v ₂、p ₂、T ₂和ρ ₂。设激波前气流马赫数为，则激波前后气流应满足如下公式。

连续性方程

（1-27）

动量方程

（1-28）

能量方程（绝热过程）

（1-29）

状态方程

（1-30）

由，可将式1-30改写成

（1-31）

在以上几个基本关系式的基础上，可导出以下的重要关系式。

（1-32）

（1-33）

（1-34）

（1-35）

（1-36）

（2）斜激波。

气流经过激波后方向发生改变，这种激波称为斜激波。v _1t、v _2t和v _1n、v _2n各表示斜激波前后速度v ₁和v ₂的切向分速度和法向分速度，α为气流折转角，β为激波角。

由于沿激波面没有切面压强的变化，所以气流经过斜激波后，沿激波面的分速度没有变化，即有

（1-37）

因为发生变化的只有法向分速度，所以斜激波相当于法向分速度的正激波。

由于斜激波是相当于法向分速度的正激波，所以只要把正激波关系式中各下标“1”、“2”换成“1n”、“2n”，则正激波有关方程式可以应用于斜激波。

对于斜激波前后气流参数之间的关系，有如下的关系式：

（1-38）

（1-39）

（1-40）

激波角β与气流折转角α之间满足如下的关系式：

（1-41）

1.1.4　流体流动分类

1．理想流体与黏性流体

如果忽略流动中流体黏性的影响，则可以近似地把流体看成是无黏的，称为无黏流体（inviscid liquid），也叫作理想流体。这时的流动称为理想流动，理想流体中没有摩擦，也就没有耗散损失。

事实上，真正的理想流体是不存在的。但是在一定的情形下，至少在特定的流动区域中，某些流体的流动非常接近于理想的流动条件，在分析处理中可以当作理想流体。

例如，在空气绕物体的流动（空气动力学）中，除去邻近与物体表面的薄层（称为边界层）之外，在其余的流动区域中，空气动力学中都处理成理想流动，此时所求解的控制方程组是不考虑黏性的欧拉方程组。

2．牛顿流体与非牛顿流体

根据内摩擦剪应力与剪应力变率的关系不同，黏性流体又可分为牛顿流体与非牛顿流体。

如果流体的剪应力与剪应力变率遵守牛顿内摩擦定律，即公式，则这种流体就称为牛顿流体。尽管这个线性的牛顿关系式只是一种近似，但是却很好地适用于一类范围很广的流体。水、空气和气体等绝大多数工业中常用的流体都是牛顿流体。

但是，对于某些物质而言，剪应力不只是速度梯度的函数（速度梯度和剪应变率是相同的），通常还可以是应变的函数，这种物质称为黏-弹性流体。剪应力只依赖于速度梯度的简单黏性流体，也可以不是牛顿流体。

事实上存在这样的流体，其剪应力与应变率之间有着相当复杂的非线性关系。如果流体的应力-应变关系还取决于实际的工况，即应变工况，则称为触变流体（如印刷油墨）。

非牛顿流体具有塑性行为，其特征是有一个表现的屈服应力，在达到表现的屈服应力之前，流体的性态像固体一样，一旦超过这个表现的屈服应力，则和黏性流体一样。

塑性流体的另一个极端情形是：在低应变率时，黏性系数很小，很容易流动，但随着应变率的增加，变得越来越像固体（如流沙）。这种流体称为膨胀流体。在图1-1中用曲线说明了牛顿流体和非牛顿流体的特征。

图1-1　牛顿流体与非牛顿流体

3．可压流体与不可压流体

流体的压缩性是指在外界条件变化时，其密度和体积发生了变化。这里的条件有两种，一种是外部压强发生了变化，另一种就是流体的温度发生了变化。描述流体的压缩性常用以下两个量。

（1）流体的等温压缩率。

当质量为m，体积为V的流体外部压强发生的变化时，相应其体积也发生了的变化，则定义流体的等温压缩率为

（1-42）

这里的负号是考虑到与总是符号相反的缘故；β的单位为1/Pa。流体等温压缩率的物理意义是，当温度不变时，每增加单位压强所产生的流体体积的相对变化率。

考虑到压缩前后流体的质量不变，式1-42还有另外一种表示形式，即

（1-43）

将理想气体状态方程代入式1-43，得到理想气体的等温压缩率为

（1-44）

（2）流体的体积膨胀系数α。

当质量为m，体积为V的物体温度发生ΔT的变化时，相应其体积也发生了ΔV的变化，则定义流体的体积膨胀系数为

（1-45）

考虑到膨胀前后流体的质量不变，式1-45还有另外一种表示形式，即

（1-46）

这里的负号是考虑到随着温度的升高，体积必然增大，则密度必然减少；α的单位为1/K。体积膨胀系数的物理意义是，当压强不变时，每增加单位温度所产生的流体体积的相对变化率。

对于物理气体，将气体状态方程代入式1-46，得到

（1-47）

在研究流体流动过程时，若考虑到流体的压缩性，则称为可压缩流动，相应地称流体为可压缩流体，如马赫数较高的气体流动。若不考虑流体的压缩性，则称为不可压缩流动，相应地称流体为不可压缩流体，如水、油等液体的流动。

4．定常流动与非定常流动

根据流体流动过程以及流动过程中的物理参数是否与时间相关，可将流动分为定常流动与非定常流动两种。

流体流动过程中，各物理量均与时间无关的流动称为定常流动。

流体流动过程中，某个或某些物理量与时间有关的流动称为非定常流动。

5．层流流动与湍流流动

流体的流动分为层流流动和湍流流动，从试验的角度来看，层流流动就是流体层与层之间相互没有任何干扰，层与层之间既没有质量的传递也没有动量的传递；而湍流流动中层与层之间相互有干扰，而且干扰的力度还会随着流动的加速而加大，层与层之间既有质量的传递，又有动量的传递。

判断流动是层流还是湍流，是看其雷诺数是否超过临界雷诺数。雷诺数的定义如下

（1-48）

式中，v为截面的平均速度，L为特征长度，为流体的运动黏度。

对于圆形管内流动，特征长度L取圆管的直径d，即

（1-49）

一般认为临界雷诺数为2 000，当Re<2 000时，管内流动是层流，否则为湍流。

对于异型管道内的流动，特征长度L取水力直径d _H，则雷诺数的计算式为

（1-50）

异型管道水力直径的定义为

（1-51）

式中，A为过流断面的面积；S为过流断面上流体与固体接触的周长，称为湿周。

1.1.5　流体流动描述的方法

描述流体物理量有两种方法，一种是拉格朗日描述，另一种是欧拉描述。

拉格朗日（Lagrange）描述也称随体描述，它着眼于流体质点，并将流体质点的物理量认为是随流体质点及时间变化的，即把流体质点的物理量表示为拉格朗日坐标及时间的函数。

设拉格朗日坐标为（a，b，c），以此坐标表示的流体质点的物理量，如矢径、速度、压强等在任意时刻t的值，便可以写为a、b、c及t的函数。

若以f表示流体质点的某一物理量，其拉格朗日描述的数学表达式为

（1-52）

例如，设时刻t流体质点的矢径，即t时刻流体质点的位置以r表示，其拉格朗日描述为

（1-53）

同样，质点速度的拉格朗日描述为

（1-54）

欧拉描述也称空间描述，它着眼于空间点，认为流体的物理量随空间点及时间而变化，即把流体物理量表示为欧拉坐标及时间的函数。设欧拉坐标为（q ₁，q ₂，q ₃），用欧拉坐标表示的各空间点上的流体物理量，如速度、压强等，在任意时刻t的值，可写为q ₁、q ₂、q ₃及t的函数。

从数学分析知道，当某时刻一个物理量在空间的分布一旦确定，该物理量在此空间形成一个场。因此，欧拉描述实际上描述了一个个物理量的场。

若以f表示流体的一个物理量，其欧拉描述的数学表达式为（设空间坐标取用直角坐标）

（1-55）

如流体速度的欧拉描述为

（1-56）

拉格朗日描述着眼于流体质点，将物理量视为流体坐标与时间的函数；欧拉描述着眼于空间点，将物理量视为空间坐标与时间的函数。它们可以描述同一物理量，必定互相相关。

设表达式表示流体质点（a，b，c）在t时刻的物理量；表达式表示空间点（x，y，z）在时刻t的同一物理量。如果流体质点（a，b，c）在t时刻恰好运动到空间点（x，y，z）上，则应有

（1-57）

（1-58）

事实上，将式1-57代入式1-58左端，即有

（1-59）

或者反解式1-57，得到

（1-60）

将式1-60代入式1-58的右端，也应有

（1-61）

由此，可以通过拉格朗日描述推出欧拉描述，同样也可以由欧拉描述推出拉格朗日描述。

1.1.6　流体力学基本方程组

流体力学基本方程组包括质量守恒方程、动量守恒方程、组分质量守恒方程、能量守恒方程、本构方程、状态方程及通用形式守恒方程。

下面对各种形式方程进行总结和对比，并分析它们之间的转化关系，以帮助读者理解流体力学基本方程组的数学物理意义，为离散计算这些方程组打下基础。

1．质量守恒方程

质量守恒方程可由4种方法得到，分别在拉格朗日法（L法）下对有限体积和体积元应用质量守恒定律、在欧拉法（E法）下对有限体积应用质量守恒定律及在直角坐标系中直接应用质量守恒定律。

（1）L法有限体积分析。

取体积为，质量为的一定流体质点团，则有

（1-62）

因为速度散度的物理意义是相对体积膨胀率及密度的随体导数，即

（1-63）

（1-64）

代入式1-62，得

（1-65）

运用奥高定理

（1-66）

得

（1-67）

式1-67即是连续性方程的积分形式。

假定被积函数连续，而且体积是任意选取的，由此可知被积函数必须等于0，即

（1-68）

或

（1-69）

在直角坐标系中，连续性方程为

（1-70）

或

（1-71）

式1-71表明，密度变化（随时间和位置）等于密度和体积变形的乘积。

（2）L法体积元分析。

考虑质量为dm的体积元，对其用拉格朗日观点，根据质量守恒定律有

（1-72）

（1-73）

两边同除以，得

（1-74）

或写成

（1-75）

式1-75表明，要维持质量守恒定律，相对体积变化率必须等于负的相对密度变化率。

（3）E法有限体积分析。

着眼坐标空间，取空间中以S面为界的有限体积，则称S面为控制面，为控制体。取外法线方向为法线的正方向，为外法线方向的单位矢量。考虑该体积内流体质量的变化，该变化主要由以下两方面原因引起。

① 通过表面S有流体流出或流入，单位时间内流出流入变化的总和为

（1-76）

② 由于密度场的不定常性（注意，欧拉观点下空间点是固定的，密度的变化只由场的不定常性刻画），单位时间内体积的质量将变化，变化量为

（1-77）

上述两者应相等，即

（1-78）

由于体积是任意的，且被积函数连续，则

（1-79）

（4）E法直角坐标系分析。

控制体体积如图1-2所示，单位时间内通过表面EFGH的通量为。

图1-2　控制体体积

通过表面ABCD的通量为

其他两对表面类似。另外，该控制体内质量的变化率为

则

（1-80）

特殊情况下的连续性方程为

定常态

（1-81）

不可压缩流体

（1-82）

2．动量守恒方程

任取一个体积为τ的流体，它的边界为S。根据动量定理，体积τ中流体动量的变化率等于作用在该体积上的质量力和面力（应力）之和。

单位面积上的面力，其中P是二阶对称应力张量，所以不是通常指的在（单位体积面元的法线方向）方向的分量。单位质量上的质量力为，则作用在该体积上的质量力和面力分别为

（1-83）

及

（1-84）

动量变化率为

（1-85）

上述动量变化率的表达式可采用如下两种处理方法。

（1）求解式1-85右边第二项内对体积元的随体导数，则

（1-86）

（2）对动量变化率表达式1-85右边第二项应用质量守恒定律

（1-87）

由上可得两种积分形式的动量方程，即

（1-88）

或

（1-89）

动量方程的微分形式为

（1-90）

或

（1-91）

微分方程中各项的物理意义为，表示单位体积上惯性力，为单位体积上的质量力，为单位体积上应力张量的散度，它是与面力等效的体力分布函数（由奥高公式转化而来）。

在直角坐标系下以应力表示的运动方程可采取下列形式。

（1-92）

或

（1-93）

这两种表达方式的等号左边实际只差了一个连续性方程，由基本微分公式

（1-94）

得

（1-95）

由连续性方程知

（1-96）

所以有

（1-97）

上述运动方程是以应力表示的黏性流体的运动方程，它们对任何黏性流体和任何运动状态都是适用的。但它们没有反映出不同属性的流体受力后的不同表现。

另外，方程数和未知量之数不等，运动方程有3个，加上连续性方程共4个，但未知量却有9个（6个应力张量分量（9个张量分量因对称关系减少为6个）和3个速度分量），所以该方程组不封闭。为使该方程组可解，必须考虑应力张量和变形速度张量之间的关系（将应力张量用速度分量表示出来），补足所需的方程。

3．本构方程

本构方程是表征流体宏观性质的一种微分方程，它用于表达流体黏性定律的应力张量和变形速度张量之间的关系。

最简单的应力与应变之间的关系是对牛顿流体作一维运动，即牛顿剪切定律可写为

（1-98）

要得到普遍意义上的广义牛顿定律需做一定的假设，但首先应理解流体速度分解定理和变形速度张量。

（1）流体速度分解定理。

刚体运动包括平动和转动两部分，一般可表为

（1-99）

其中是刚体中选定一点上的平动速度，是刚体绕点转动的瞬时角速度矢量，r是要确定速度那一点到点的矢径。转动角速度可用表示。因为

（1-100）

故

（1-101）

流体运动除平动、转动外，还有变形运动。设微团内点的速度为，邻域内任一点的速度为。将在点泰勒展开并略去二阶无穷小项，得

（1-102）

显然，是一个二阶张量（局部速度梯度张量），由张量分解定理可将该张量分解成对称张量S和反对称张量A之和，于是

（1-103）

所以

（1-104）

式1-104右边第二项、第三项可具体表示为

（1-105）

及

（1-106）

其中

（1-107）

另外

；

（1-108）

所以

（1-109）

式1-109表明流体运动可分为平动、转动和变形3种形式，S称为变形速度张量，该定理称为亥姆霍兹（Helmholtz）速度分解定理。

另外，流变学中常用应变速率张量来表示流体的变形和拉伸（或压缩），而用转动张量表示转动，它们与流体力学中的变形速度张量和转动张量的关系为

（1-110）

（1-111）

（1-112）

（2）变形速度张量的物理意义。

变形速度的表达式为

（1-113）

经分析可得

；

（1-114）

其中、及是角变形速率，亦称剪切应变速率（称拉伸应变速率）。

由上可知，变形速度张量的对角线分量、、的物理意义分别是x、y、z轴线上线段元dx、dy、dz的相对拉伸速度或相对压缩速度。

而非对角线分量的物理意义分别是y轴与z轴、z轴与x轴、x轴与y轴之间夹角的剪切速率的负值。

（3）广义牛顿定律及基本假设。

① 运动流体的应力张量在运动停止后应趋于静止流体的应力张量。据此将应力张量写成各向同性部分和各向异性部分P＇是方便的，因此

（1-115）

是除去后得到的张量，称为偏应力张量。当运动消失时它趋于0。可见，偏应力张量和应力张量一样也是对称张量。

② 偏应力张量的各分量是局部速度梯度张量各分量的线性齐次函数。当速度在空间均匀分布时，偏应力张量为0；当速度偏离均匀分布时，在黏性流体中产生了偏应力，它力图使速度恢复到均匀分布情形。

③ 流体是各向同性的，即流体性质不依赖于方向或坐标系的转换。

根据假设②，有

（1-116）

显然是一个四阶张量，它是表征流体黏性的常数，共个。根据假设③，是各向同性张量且对、对称，故

（1-117）

观察式1-117可知，对k、l也是对称的，物性常数减少至只有2个，即第二黏度和黏度。将式1-117代入偏应力表达式1-116（反对称项为零），得

（1-118）

则应力张量为

（1-119）

引进

（1-120）

则

（1-121）

根据式1-121，可知

（1-122）

（1-123）

（1-124）

将式1-122、式1-123和式1-124相加，得

（1-125）

对于可压缩流体，流体的体积在运动过程中发生膨胀或收缩，将引起平均法应力（由奥高公式可证某固定点处所有方向上法应力的平均值等于x、y、z三个方向上法应力的平均值，这是一个不随坐标系改变的不变量）的值发生的改变，称为第二黏性系数，亦称膨胀黏性系数。应用斯托克斯假定，即，则本构方程为

（1-126）

（1-127）

（1-128）

一般处理的是不可压缩流体，则

（1-129）

（1-130）

（1-131）

（1-132）

在直角坐标系下有

（1-133）

这里，有的文献中将应力张量用表示。将上述应力张量与变形速度张量的关系式代入运动方程，得

（1-134）

即

（1-135）

写成直角坐标系下的形式为

（1-136）

有些文献中常把式1-136等号右边表示分子黏性作用的3项做如下变化，以第一式为例。

（1-137）

其中

（1-138）

据此，有

（1-139）

及

（1-140）

其中广义源项定义为

（1-141）

当流体黏度不变且不可压缩时（牛顿流体），有

（1-142）

所以运动方程简化为

（1-143）

其中是运动黏度，亦是动量扩散系数，单位为m²/s。

本构方程和运动方程是紧密联系在一起的，通过本构方程可将应力张量用变形速度张量表示出来，即应力可用应变速率表示，而应变速率实际由速度分量决定，故使运动方程和连续性方程原则上封闭可解。

注意:

这里讨论的本构方程仅局限于牛顿流体。符合广义牛顿定律的流体称为牛顿流体，否则称为非牛顿流体。非牛顿流体的本构方程不能用广义牛顿定律描述，如对聚合物溶液等流体应该参考相关文献。

4．能量守恒方程

由能量守恒定律知，体积内流体的动能和内能的变换率等于单位时间内质量力和表面力所做的功加上单位时间内给予体积的热量。

体积内流体的动能和内能的总和为

（1-144）

其中是单位体积内的流体内能。质量力对体积τ内流体所做的功（单位时间内移动距离v，点积求做功）为

（1-145）

表面力对体积τ内流体所做的功（单位时间内移动距离v，点积求做功）为

（1-146）

单位时间内以热传导方式通过表面传给体积τ的热量为

（1-147）

式1-147中的被积函数实际就是傅里叶热传导定律，即热流密度矢量正比于传热面法向温度梯度。单位时间内由于辐射或其他原因（反应、蒸发等）传入τ内的总热量（q为单位时间内传入单位质量的热量分布函数）为

（1-148）

将式1-144～式1-148进行守恒计算，得

（1-149）

这是积分形式的能量守恒方程。求解体积分的随体导数并运用奥高公式把面积分转化为体积分，可得到微分形式的能量守恒方程，即

（1-150）

由质量守恒定律可得

（1-151）

所以

（1-152）

另外

（1-153）

（1-154）

则能量方程的微分形式为

（1-155）

或

（1-156）

或

（1-157）

式1-157各项的物理意义为：左边第一、第二项代表内能和动能的随体导数，右边第一项是单位体积内的质量力做的功，第二项是单位体积内面力所做的功，第三项是单位体积内热传导输入的热量，最后一项表示由于辐射或其他物理或化学原因的热量贡献。

能量守恒方程的另一种形式为

（1-158）

式1-158的物理意义为：单位体积内由于流体变形，面力所做的功加上热传导及辐射等其他原因传入的热量恰好等于单位体积内的内能在单位时间内的增加。将该式进一步简化，有

（1-159）

设

（1-160）

为由于黏性作用机械能转化为热能的部分，称为耗散函数。另外，在考虑液体流体时，比焓与内能值可看作相等，即，压力不做功，则

（1-161）

所以有

（1-162）

其中是单位体积内热源或由于辐射或其他物理或化学原因的热量贡献。一般较小，可以忽略。对液体及固体可以取，进一步取为常数，并把耗散函数纳入源项，于是

（1-163）

对于不可压缩流体，有

（1-164）

对于可以忽略黏性耗散作用的稳态低速流，能量方程可以简化为

（1-165）

及

（1-166）

取速度为0，则可得到稳态的热传导方程（对流项消失）为

（1-167）

5．状态方程

由连续性方程、运动方程、能量方程确定的未知量有u、v、w、p、T、 共6个，但方程数只有5个，为使方程组封闭，需补充一个联系、的状态方程：

（1-168）

6．组分质量守恒方程

在一个特定的系统中可能存在质的交换，或者存在多种化学组分，每一种组分都需要遵守组分质量守恒定律，即系统内某种化学组分对时间的变化率等于通过系统界面的净扩散流量与由反应产生的生成率之和，可表示为

（1-169）

其中代表单位体积内组分的质量变化率，是组分的对流流量密度。代表扩散流量密度，它由费克定律给出。是单位体积内组分的生成率。费克定律可表示为

（1-170）

其中为扩散系数。将扩散定律代入守恒方程，得

（1-171）

为便于以后引用，表1-3列出了本节的守恒型控制方程。

表1-3　常用流动与传热问题的守恒型控制方程

1.1.7　湍流模型

描述流体运动（层流）的流体力学基本方程组是封闭的，而描述湍流运动的方程组由于采用了某种平均（时间平均或网格平均等）而不封闭，必须对方程组中出现的新未知量采用模型而使其封闭，这就是CFD中的湍流模型。

湍流模型的主要作用是将新未知量和平均速度梯度联系起来。目前，工程应用中湍流的数值模拟主要分三大类：直接数值模拟（DNS）、大涡模拟（LES）和基于雷诺平均N-S方程组（RANS）的模型。

1．直接数值模拟（DNS）

直接数值模拟（DNS）方法直接求解湍流运动的N-S方程，得到湍流的瞬时流场，即各种尺度的随机运动，可以获得湍流的全部信息。

随着现代计算机的发展和先进数值方法的研究，DNS方法已经成为解决湍流的一种实际的方法。但由于计算机条件的约束，目前只能限于一些低Re数的简单流动，不能用于处理工程中的复杂流动问题。

目前国际上正在做的湍流直接数值模拟还只限于较低的雷诺数（Re约200）和非常简单的流动外形，如平板边界层、完全发展的槽道流以及后台阶流动等。

2．大涡模拟（LES）

大涡模拟（LES）方法即对湍流脉动部分的直接模拟，将N-S方程在一个小空间域内进行平均（或称之为滤波），以便从流场中去掉小尺度涡，导出大涡所满足的方程。

小涡对大涡的影响会出现在大涡方程中，再通过建立模型（亚格子尺度模型）来模拟小涡的影响。

由于湍流的大涡结构强烈地依赖于流场的边界形状和边界条件，难以找出普遍的湍流模型来描述具有不同的边界特征的大涡结构，所以宜做直接模拟。

相反地，由于小尺度涡对边界条件不存在直接依赖关系，而且一般具有各向同性性质，所以亚格子尺度模型具有更大的普适性，比较容易构造，这是它比雷诺平均方法要优越的地方。

LES方法已经成为计算湍流的最强有力的工具之一，应用的方向也在逐步扩展，但是仍然受计算机条件等的限制，使之成为解决大量工程问题的成熟方法仍有很长的路要走。

3．基于雷诺平均N-S方程组（RANS）的模型

目前能够用于工程计算的方法就是模式理论。所谓湍流模式理论，就是依据湍流的理论知识、实验数据或直接数值模拟结果，对雷诺应力做出各种假设，即假设各种经验的和半经验的本构关系，从而使湍流的平均雷诺方程封闭。

从对模式处理的出发点不同，可以将湍流模式理论分成两大类：一类称为二阶矩封闭模式（雷诺应为模式），另一类称为涡黏性封闭模式。

（1）雷诺应力模式。

雷诺应力模式即二阶矩封闭模式，是从雷诺应力满足的方程出发，将方程右端未知的项（生成项、扩散项、耗散项等）用平均流动的物理量和湍流的特征尺度表示出来。

典型的平均流动的变量是平均速度和平均温度的空间导数。这种模式理论由于保留了雷诺应力所满足的方程，如果模拟得好，可以较好地反映雷诺应力随空间和时间的变化规律，因而可以较好地反映湍流运动规律。

因此，二阶矩封闭模式是一种较高级的模式，但是，由于保留了雷诺应力的方程，加上平均运动的方程，整个方程组总计15个方程，是一个庞大的方程组，应用这样一个庞大的方程组来解决实际工程问题，计算量很大，这就极大地限制了二阶矩封闭模式在工程问题中的应用。

（2）涡黏性封闭模式。

在工程湍流问题中得到广泛应用的模式是涡黏性封闭模式。这是由Boussinesq仿照分子黏性的思路提出的，即设雷诺应力为

（1-172）

这里是湍动能，称为涡黏性系数，这是最早提出的基准涡黏性封闭模式，即假设雷诺应力与平均速度应变率成线性关系。当平均速度应变率确定后，6个雷诺应力只需要通过确定一个涡黏性系数就可完全确定，且涡黏性系数各向同性，可以通过附加的湍流量来模化，如湍动能k、耗散率、比耗散率ω以及其他湍流量、、。根据引入的湍流量的不同，可以得到不同的涡黏性模式，如常见的、k-ω模式，以及后来不断得到发展的、q-w、k-l等模式，涡黏性系数可以分别表示为

，， ，，

为了使控制方程封闭，引入多少个附加的湍流量，就要同时求解多少个附加的微分方程，根据求解的附加微分方程的数目，一般可将涡黏性模式划分为4类：零方程模式、半方程模型、一方程模式和两方程模式。

① 零方程模式。

所谓零方程模式，就是试图直接用平均流动物理量模化v_T，而不引入任何湍流量（如k、ε等）。

例如，Prandtl的混合长理论就是一种零方程模式：

（1-173）

式中，称为混合长。

在零方程模式的框架下，得到最为广泛应用的是Baldwin-Lomax模式（BL模式）。该模式是对湍流边界层的内层和外层采用不同的混合长假设。

这是因为靠近壁面处，湍流脉动受到很大的抑制，含能涡的尺度减小很多，因此长度尺度减小很多；另一方面，在边界层外缘，湍流呈间歇状，质量、动量和能量的输运能力大大下降，即湍流的扩散能力减小。

这样，应用混合长理论来确定涡黏性系数在这两个不同的区域应该有不同的形式。Baldwin-Lomax模式的具体数学描述如下。

（1-174）

这里y_c是 (v_T)_in=(v_T)_out的离壁面最小距离y值。

对于内层，即y≤yc，有

(v_T)_in=l ²Ω

（1-175）

Ω是涡量，，l是长度尺度。

（1-176）

式中，k=0.4是Karman常数，A⁺是模化常数，y⁺是无量纲法向距离。

（1-177）

其中是摩擦速度，其含义为

（1-178）

此处下标w表示壁面。

对于外层，即y>y_c，有

（1-179）

其中

（1-180）

F _max是函数F(y)的最大值。

（1-181）

而y _max是F(y)达到最大值的位置。F _kleb是Klebanoff间歇函数。

（1-182）

U _dif是平均速度分布中最大值和最小值之差。

几个模化常数的值如下。

由上述模化关系可以看出，雷诺应力完全由当时当地的平均流参数用代数关系式决定。

平均流场的任何变化立刻为当地的湍流所感知，这表明零方程模式是一个平衡态模式，假定湍流运动永远处于与平均运动的平衡之中。

实际上对于大多数湍流运动而言并非如此，特别是对于平均流空间和时间有剧烈变化的情形，再有因为坐标y显式地出现在湍流模式中，零方程模式不具有张量不变性，所以将它应用到复杂几何外形的流动的数值模拟会带来困难。

当流动发生分离时，Baldwin-Lomax模式会遇到困难，这是因为在分离点和再附点附近，摩擦速度为0，此时要引入一些人为的干涉来消除这些困难。

计算实践表明，只要流动是附体的，零方程模式一般都可以较好地确定压强分布，但是摩阻和传热率的估算不够准确，特别是当流动有分离和再附时。这是因为附体流压强分布对湍流应力不敏感。

总之，对于附体流动，如果只关心压强分布，应用零方程模式通常可以给出满意的结果，而且模式应用起来十分简便。但是对于计算摩阻的需求，零方程模式是不能满足要求的。对于有分离、再附等复杂流动，零方程模式是不适用的。

② 半方程模式。

为了能计算具有较强压强梯度，特别是较强逆压梯度的非平衡湍流边界层，Johnson-King于1985年提出了一个非平衡代数模式（JK模式），该模式仍采用涡黏性假设，把涡黏性的分布与最大剪切应力联系在一起，内层涡黏性与外层涡黏性分布用一个指数函数作为光滑拟合，外层涡黏性系数作为一个自由参数，由描述最大剪切应力沿流向变化的常微分方程来确定，此常微分方程是由湍流动能方程导出的，故此模式又称为半方程模式。

JK模式虽然仍采用涡黏性假设，却包含雷诺应力模式的特点。由于求解常微分方程比一方程、两方程模式中求解偏微分方程要简单、省时得多，故用JK模式的工作量只略高于通常平衡状态的零方程代数模式的工作量。

③ 一方程模式。

Baldwin-Barth（BB）模式是在两方程模型中，将某一导出的应变量作为基本物理量而得到的，应用此一方程模式可避免求解两方程时会遇到的某些数值困难。BB 一方程模式所选择的导出应变量为“湍流雷诺数”Re_t。

BB模式对计算网格的要求低，壁面的网格可以与采用BL代数模式的相当，而不像两方程k-ε模式那样要求壁面网格很细，这样就避免了在k-ε模式中流场求解的刚性问题。

Spalart-Allmaras（SA）模式与BB模式不同，它不是直接利用k-ε两方程模式加以简化而得，而是从经验和量纲分析出发，由针对简单流动再逐渐补充发展，进而适用于带有层流流动的固壁湍流流动的一方程模式。该模式中选用的应变量是与涡黏性相关的量，除在黏性次层外，与是相等的。

上述两种一方程模式具有相似的特点，它们不像代数模式那样需要分为内层模式、外层模式或壁面模式、尾流模式，同时也不需要沿法向网格寻找最大值，因此易于应用到非结构网格中；但由于在每个时间步长内，需要对整个流场求解一组偏微分方程，故比BL和JK模式更费机时。

④ 两方程模式。

常用的两方程模式有标准k-ε两方程模式、可实现型k-ε两方程模式、低雷诺数k-ε模式及k-ω两方程模式等。

● 标准k-ε两方程模式。

k-ε模式是最为人所知和应用最广泛的两方程涡黏性模式，为积分到壁面的不可压缩/可压缩湍流的两方程涡黏性模式。

雷诺应力的涡黏性模式为

（1-183）

其中μ_t为涡黏性（eddy viscosity），S_ij为平均速度应变率张量（mean-velocity strain-rate tensor），ρ为流体密度，k为湍动能，δ_ij为克罗内克算子（Kronecker delta）。涡黏性定义为湍动能k和湍流耗散率ε的函数。

（1-184）

基于量纲分析，涡黏性由流体密度ρ、湍流速度尺度k²和长度尺度k^3/2/ε来标度，衰减函数f由湍流雷诺数来模化。

湍流输运方程可表示成湍流能量输运方程1-185和能量耗散输运方程1-186。

（1-185）

（1-186）

其中右端项分别表示生成项（production term）、耗散项（dissipation term）和壁面项（wall term）。

模式中各常数的值如下。

=0.09　 　

近壁衰减函数

（1-187）

（1-188）

壁面项

（1-189）

（1-190）

其中u _s为平行于壁面的流动速度。

积分到壁面的无滑移边界条件为

k= 0　　　　ε = 0

● 可实现型k-ε 两方程模式。

上述标准k - ε模式，对于高平均切变率流动会出现非物理的结果（如当Sk/ε>3.7时，其中 ）。为了保证模式的可实现性，模式函数C_μ不应该是常数，而应当是平均切变率的函数。实验表明，对边界层流动和均匀切变流，C_μ的值是不同的。为此，根据可实现性对模式的约束条件，建议采用以下形式的C_μ。

（1-191）

式1-191中

（1-192）

是在以角速度ω_k旋转的旋转坐标系中得到的平均旋转速率。

（1-193）

关系式1-191中唯一未确定的系数是A ₀。为简单起见，可以设其为常数。对边界层流动，可以取A ₀= 4.0；对于其他流动，A ₀的数值可以调节。

● 低雷诺数k-ε 模式。

上述几种k-ε模式适用于高雷诺数情形。但是对于近壁区，湍流雷诺数很低，对湍流动力学而言，黏性效应非常重要，此时湍流雷诺数的效应必须加以考虑。研究摩阻的计算关注的恰恰是近壁区，因此低雷诺数k-ε模式的研究是十分重要的。

现将有关结果整理如下。

低雷诺数下的涡黏性和k-ε模式方程为

（1-194）

（1-195）

（1-196）

式中

（1-197）

所有模化常数如下。

（1-198）

其中

（1-199）

此处，f_u和f ₁、f ₂称为阻尼函数，是用来反映近壁区低雷诺数效应的一个经验公式，系数 、见表1-4。

表1-4　系数 、

● k-ω 两方程模式。

k-ω 模式是积分到壁面的不可压缩/可压缩湍流的两方程涡黏性模式，最主要的文献来自Wilcox。

下面求解湍动能k和它的的对流输运方程。

雷诺应力的涡黏性模型为

（1-200）

这里μ _t为涡黏性，S _ij为平均速度应变率张量，ρ为流体密度，k为湍动能，δ _ij为克罗内克算子。涡黏性定义为湍动能k和比耗散率ω的函数。

（1-201）

k和ω的输运方程为

（1-202）

（1-203）

模式中各常数的值如下。

对于边界层流动，壁面无滑移边界条件为

k= 0

（1-204）

这里y₁为离开壁面第一个点的距离，且<1。

除以上介绍的FLUENT中常用的湍流模型外，还有众多诸如SST两方程模式、k-τ模型、q-ω模型、双尺度两方程模式等湍流模型，在此不再一一详细介绍，感兴趣的读者可以查阅相关文献。

在实际求解中，选用什么模型要根据具体问题的特点来决定。选择的一般原则是精度高，应用简单，节省计算时间，同时也具有通用性。

不同软件所包含的湍流模型有略微区别，但常用的湍流模型一般CFD软件都包含。图1-3为常用的湍流模型及其计算量的变化趋势。

图1-3　CFD软件中常用湍流模型及其计算量变化趋势

1.2　计算流体力学（CFD）基础

计算流体动力学（Computational Fluid Dynamics，CFD）是近代流体力学、数值数学和计算机科学结合的产物，是一门具有强大生命力的边缘科学。

1.2.1　CFD概述

CFD以电子计算机为工具，应用各种离散化的数学方法，对流体力学的各类问题进行数值实验、计算机模拟和分析研究，以解决各种实际问题。

计算流体力学和相关的计算传热学、计算燃烧学的原理是用数值方法求解非线性联立的质量、能量、组分、动量和自定义的标量的微分方程组，求解结果能预报流动、传热、传质、燃烧等过程的细节，并成为过程装置优化和放大定量设计的有力工具。计算流体力学的基本特征是数值模拟和计算机实验，它从基本物理定理出发，在很大程度上替代了耗资巨大的流体动力学实验设备，在科学研究和工程技术中产生巨大的影响。

计算流体力学是目前国际上一个强有力的研究领域，是进行传热、传质、动量传递及燃烧、多相流和化学反应研究的核心和重要技术，广泛应用于航天设计、汽车设计、生物医学工业、化工处理工业、涡轮机设计、半导体设计、HAVC&R 等诸多工程领域，板翅式换热器 设计是CFD技术应用的重要领域之一。

CFD在最近20 年中得到飞速的发展，除了计算机硬件工业的发展给它提供了坚实的物质基础外，还主要因为无论分析的方法或实验的方法都有较大的限制。例如，由于问题的复杂性，既无法作分析解，也因费用昂贵而无力进行实验确定，而CFD 的方法正具有成本低和能模拟较复杂或较理想的过程等优点。

经过一定考核的CFD软件可以拓宽实验研究的范围，减少成本昂贵的实验工作量。在给定的参数下用计算机对现象进行一次数值模拟相当于进行一次数值实验，历史上也曾有过首先由CFD数值模拟发现新现象而后由实验予以证实的例子。

CFD软件一般都能推出多种优化的物理模型，如定常和非定常流动、层流、紊流、不可压缩和可压缩流动、传热、化学反应等。对每一种物理问题的流动特点，都有适合它的数值解法，用户可选择显式或隐式差分格式，以期在计算速度、稳定性和精度等方面达到最佳。

CFD 软件之间可以方便地进行数值交换，并采用统一的前、后处理工具，这就省去了科研工作者在计算机方法、编程、前后处理等方面投入的重复、低效的劳动，而可以将主要精力和智慧用于物理问题本身的探索上。

1.2.2　CFD求解力学问题的过程

所有CFD问题的求解过程都可用图1-4表示。如果所求解的问题是瞬态问题，则可将图1-4的过程理解为一个时间步的计算过程，循环这一过程求解下个时间步的解。下面对各求解步骤进行简单介绍。

图1-4　CFD求解流程框图

1．建立控制方程

建立控制方程是求解任何问题前都必须首先进行的。一般来讲，这一步是比较简单的。因为对于一般的流体流动而言，可直接写出其控制方程。假定没有热交换发生，则可直接将连续方程与动量方程作为控制方程使用。一般情况下，需要增加湍流方程。

2．确定边界条件和初始条件

初始条件与边界条件是控制方程有确定解的前提，控制方程与相应的初始条件、边界条件的组合构成对一个物理过程完整的数学描述。

初始条件是所研究对象在过程开始时刻各个求解变量的空间分布情况。对于瞬态问题，必须给定初始条件。对于稳态问题，不需要初始条件。

边界条件是在求解区域的边界上所求解的变量或其导数随地点和时间的变化规律。对于任何问题，都需要给定边界条件。

3．划分计算网格

采用数值方法求解控制方程时，都是想办法将控制方程在空间区域上进行离散，然后求解得到离散方程组。要想在空间域上离散控制方程，必须使用网格。现已发展出多种对各种区域进行离散以生成网格的方法，这些方法统称为网格生成技术。

不同的问题采用不同数值解法时，所需要的网格形式是有一定区别的，但生成网格的方法基本是一致的。目前网格分结构网格和非结构网格两大类。

简单地讲，结构网格在空间上比较规范，如对一个四边形区域，网格往往是成行成列分布的，行线和列线比较明显。而非结构网格在空间分布上没有明显的行线和列线。

对于二维问题，常用的网格单元有三角形和四边形等形式；对于三维问题，常用的网格单元有四面体、六面体、三菱体等形式。在整个计算域上，网格通过节点联系在一起。

目前各种CFD软件都配有专用的网格生成工具，如FLUENT 使用Gambit作为前处理软件。多数CFD软件可接收采用其他CAD或CFD/FEM软件产生的网格模型。例如，FLUENT可以接收ANSYS所生成的网格。

4．建立离散方程

对于在求解域内所建立的偏微分方程，理论上是有真解（或称精确解或解析解）的。但由于所处理问题自身的复杂性，一般很难获得方程的真解。因此就需要通过数值方法把计算域内有限数量位置（网格节点或网格中心点）上的因变量值当作基本未知量来处理，从而建立一组关于这些未知量的代数方程组，然后通过求解代数方程组来得到这些节点值，而计算域内其他位置上的值则根据节点位置上的值来确定。

由于所引入的应变量在节点之间的分数假设及推导离散化方程的方法不同，所以形成了有限差分法、有限元法、有限体积法等不同类型的离散化方法。

对于瞬态问题，除了在空间域上的离散外，还要涉及在时间域上的离散。离散后，将要涉及使用何种时间积分方案的问题。

5．离散初始条件和边界条件

前面所给定的初始条件和边界条件是连续性的，如在静止壁面上速度为0，现在需要针对所生成的网格，将连续型的初始条件和边界条件转化为特定节点上的值，如静止壁面上共有90个节点，则这些节点上的速度值应均设为0。

商用CFD软件往往在前处理阶段完成网格划分后，直接在边界上指定初始条件和边界条件，然后由前处理软件自动将这些初始条件和边界条件按离散的方式分配到相应的节点上。

6．给定求解控制参数

在离散空间上建立了离散化的代数方程组，并施加离散化的初始条件和边界条件后，还需要给定流体的物理参数和湍流模型的经验系数等。此外，还要给定迭代计算的控制精度、瞬态问题的时间步长和输出频率等。

7．求解离散方程

进行上述设置后，生成了具有定解条件的代数方程组。对于这些方程组，数学上已有相应的解法，如线性方程组可采用Gauss消去法或Gauss-Seidel迭代法求解，而对于非线性方程组，可采用Newton-Raphson方法。

商用CFD软件往往提供多种不同的解法，以适应不同类型的问题。这部分内容属于求解器设置的范畴。

8．显示计算结果

通过上述求解过程得出了各计算节点上的解后，需要通过适当的手段将整个计算域上的结果表示出来，这时，可采用线值图、矢量图、等值线图、流线图、云图等方式来表示计算结果。

• 线值图是指在二维或三维空间上，将横坐标取为空间长度或时间历程，将纵坐标取为某一物理量，然后用光滑曲线或曲面在坐标系内绘制出某一物理量沿空间或时间的变化情况。
• 矢量图是直接给出二维或三维空间里矢量（如速度）的方向及大小，一般用不同颜色和长度的箭头表示速度矢量。矢量图可以比较容易地让用户发现其中存在的漩涡区。
• 等值线图是用不同颜色的线条表示相等物理量（如温度）的一条线。
• 流线图是用不同颜色的线条表示质点运动轨迹。
• 云图是使用渲染的方式，将流场某个截面上的物理量（如压力或温度）用连续变化的颜色块表示其分布。

1.2.3　CFD数值模拟方法和分类

CFD的数值解法有很多分支，这些方法之间的区别主要在于对控制方程的离散方式。根据离散原理的不同，CFD大体上可以分为有限差分法（FDM）、有限元法（FEM）和有限体积法（FVM）。

1．有限差分法

有限差分法（FDM）是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。

有限差分法以Taylor级数展开的方法，把控制方程中的导数用网格节点上的函数值的差商代替，从而创建以网格节点上的值为未知数的代数方程组。

该方法是一种直接将微分问题变为代数问题，从而可以用近似数值解法求解，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。

从有限差分格式的精度来划分，有一阶格式、二阶格式和高阶格式；从差分的空间形式来考虑，可分为中心格式和逆风格式；考虑时间因子的影响，差分格式还可分为显格式、隐格式、显隐交替格式等。

目前常见的差分格式主要是上诉几种格式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际情况和柯郎稳定条件决定。

2．有限元法

有限元法（FEM）的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。

采用不同的权函数和插值函数形式，便于构成不同的有限元方法。有限元法最早应用于结构力学，后来随着计算机的发展逐渐用于流体力学的数值模拟。

在有限元法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线性组合来逼近单元中的真解，整个计算域上总体的基函数可以看作由每个单元基函数组成，而整个计算域内的解可以看作由所有单元上的近似解构成。

有限元方法的基本思路和解题步骤可归纳如下。

（1）建立积分方程。根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。

（2）区域单元剖分。根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，并列出自然边界和本质边界的节点序号和相应的边界值。

（3）确定单元基函数。根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元具有规则的几何形状，所以在选取基函数时可遵循一定的法则。

（4）单元分析。将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数（即单元中各节点的参数值）的代数方程组（称为单元有限元方程）。

（5）总体合成。在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进行累加，形成总体有限元方程。

（6）边界条件的处理。一般边界条件有3种形式，分别为本质边界条件（狄里克雷边界条件）、自然边界条件（黎曼边界条件）、混合边界条件（柯西边界条件）。对于自然边界条件，一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法则对总体有限元方程进行修正满足。

（7）解有限元方程。根据边界条件修正的总体有限元方程组，是含所有待定未知量的封闭方程组，采用适当的数值计算方法求解，可求得各节点的函数值。

3．有限体积法

有限体积法（Finite Volume Method，FVM）又称为控制体积法。其基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。

其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分，必须假定值在网格点之间的变化规律，即假设值分段分布的剖面。

从积分区域的选取方法看来，有限体积法属于加权剩余法中的子区域法；从未知解的近似方法看来，有限体积法属于采用局部近似的离散方法。简而言之，子区域法属于有限体积法的基本方法。

有限体积法的基本思路易于理解，并能得出直接的物理解释。离散方程的物理意义就是因变量在有限大小的控制体积中的守恒原理，如同微分方程表示因变量在无限小的控制体积都得到满足，在整个计算区域，自然也就得到满足一样，这是有限体积法吸引人的优点。

某些离散方法，如有限差分法，仅当网格极其细密时，离散方程才满足积分守恒；而有限体积法即使在粗网格情况下，也显示出准确的积分守恒。

就离散方法而言，有限体积法可视为有限单元法和有限差分法的中间物，有限单元法必须假定值符号网格点之间的变化规律（即插值函数），并将其作为近似解；有限差分法只考虑网格点上的数值而不考虑其在网格点之间如何变化；有限体积法只寻求节点值，这与有限差分法相类似，但有限体积法在寻求控制体积的积分时，必须假定值在网格点之间的分布，这又与有限单元法相类似。

在有限体积法中，插值函数只用于计算控制体积的积分，得出离散方程后，便可忘掉插值函数；如果需要的话，可以对微分方程中不同的项采取不同的插值函数。

1.2.4　有限体积法计算区域的离散

从前面的介绍中可以看出，有限体积法是一种分块近似的计算方法，因此其中比较重要的步骤是计算区域的离散和控制方程的离散。

所谓区域的离散化，实际上就是用一组有限个离散的点来代替原来的连续空间。一般的实施过程是：把所计算的区域划分成许多个互不重叠的子区域（sub-domain），确定每个子区域中的节点位置及该节点所代表的控制体积。区域离散后，得到以下4种几何要素。

• 节点：需要求解的未知物理量的几何位置。
• 控制体积：应用控制方程或守恒定律的最小几何单位。
• 界面：定义了与各节点相对应的控制体积的界面位置。
• 网格线：连接相邻两节点面形成的曲线簇。

一般把节点看成是控制体积的代表。在离散过程中，将一个控制体积上的物理量定义并存储在该节点处。一维问题的有限体积法计算网格如图1-5所示，二维问题的有限体积法计算网格如图1-6所示。

图1-5　一维控制体积网格

图1-6　二维控制体积网格

计算区域离散的网格有两类：结构化网格和非结构化网格。

结构化网格（structured grid）的节点排列有序，即给出了一个节点的编号后，立即可以得出其相邻节点的编号，所有内部节点周围的网格数目相同。

结构化网格具有实现容易、生成速度快、网格质量好、数据结构简单化的优点，但不能实现复杂边界区域的离散。

非结构化网格的内部节点以一种不规则的方式布置在流场中，各节点周围的网格数目不尽相同。这种网格虽然生成过程比较复杂，但却有极大的适应性，对复杂边界的流场计算问题特别有效。

1.2.5　有限体积法控制方程的离散

前面给出的流体流动问题的控制方程，无论是连续性方程、动量方程，还是能量方程，都可写成如式1-205所示的通用形式。

（1-205）

对于一维稳态问题，其控制方程如式1-206所示。

（1-206）

式1-206中，从左到右各项分别为对流项、扩散项和源项。方程中的是广义变量，可以是速度、温度或浓度等一些待求的物理量。是相应于的广义扩散系数，是广义源项。变量在端点A和B的边界值为已知。

有限体积法的关键一步是在控制体积上积分控制方程，在控制体积节点上产生离散的方程。对一维模型方程1-206，在图1-5所示的控制体积P上进行积分，有

（1-207）

式1-207中，是控制体积的体积值。当控制体积很微小时，可以表示为，这里A是控制体积界面的面积。从而有

（1-208）

从式1-208可以看到，对流项和扩散项均已转化为控制体积界面上的值。有限体积法最显著的特点之一就是离散方程中具有明确的物理插值，即界面的物理量要通过插值的方式由节点的物理量来表示。

为了建立所需形式的离散方程，需要找出如何表示式1-208中界面e和w处的ρ、u、Γ、和。在有限体积法中规定，ρ、u、Γ、__和等物理量均是在节点处定义和计算的。因此，为了计算界面上的这些物理参数（包括其导数），需要一个物理参数在节点间的近似分布。

可以想象，线性近似是可以用来计算界面物性值的最直接，也是最简单的方式。这种分布叫作中心差分。如果网格是均匀的，则单个物理参数（以扩散系数为例）的线性插值结果是

（1-209）

的线性插值结果是

（1-210）

与梯度项相关的扩散通量的线性插值结果是

（1-211）

对于源项S，它通常是时间和物理量的函数。为了简化处理，将S转化为如下线性　 方式：

（1-212）

式中，是常数，是随时间和物理量变化的项。将式1-210～式1-212代入式1-208，有

（1-213）

整理后得到

（1-214）

记为

式1-214中

（1-215）

对于一维问题，控制体积界面e和w处的面积和均为1，即单位面积。这样，式1-215各系数可转化为

（1-216）

根据，每个节点上都可建立此离散方程，通过求解方程组，就可得到各物理量在各节点处的值。

为了方便，定义两个新的物理量F和D，其中F表示通过界面上单位面积的对流质量通量（convective mass flux），简称对流质量流量，D表示界面的扩散传导性（diffusion con{\text{d}}uctance）。定义表达式为

（1-217）

这样，F和D在控制界面上的值分别为

（1-218）

在此基础上，定义一维单元的Peclet数Pe为

（1-219）

式1-219中，表示对流与扩散的强度之比。当为0时，对流扩散演变为纯扩散问题，即流场中没有流动，只有扩散；当>0时，流体沿x方向流动，当数很大时，对流扩散问题演变为纯对流问题。一般在中心差分格式中，有<2的要求。

将式1-218代入式1-216中，有

（1-220）

瞬态问题与稳态问题相似，主要是瞬态项的离散。其一维瞬态问题的通用控制方　 程为

（1-221）

该方程是一个包含瞬态及源项的对流扩散方程。从左到右，方程中的各项分别是瞬态项、对流项、扩散项及源项。方程中的是广义变量，如速度分量、温度、浓度等，为相应于的广义扩散系数，为广义源项。

对瞬态问题用有限体积法求解时，在将控制方程对控制体积进行空间积分的同时，还必须对时间间隔进行时间积分。对控制体积所作的空间积分与稳态问题相同，这里仅叙述对时间的积分。

将式1-221在一维计算网格上对时间及控制体积进行积分，有

（1-222）

整理后，有

（1-223）

式1-223中，A是图1-5中控制体积P的界面处的面积。

在处理瞬态项时，假定物理量在整个控制体积P上均具有节点处值，并用线性插值来表示。源项也分解为线性方式。对流项和扩散项的值按中心差分格式通过节点处的值来表示，则有

（1-224）

假定变量对时间的积分为

（1-225）

式1-225中，上标0代表t时刻；是时刻的值；f为0与1之间的加权因子，当f=0时，变量取老值进行时间积分，当f=1时，变量采用新值进行时间积分。将、、及进行时间积分，由式1-224可得到

（1-226）

整理后得

（1-227）

扩展F和D的定义，即乘以面积A，有

（1-228）

代入方程1-227，得

（1-229）

同样也向稳态问题引入、和，式1-229变为

（1-230）

式1-230中

（1-231）

根据f的取值，瞬态问题对时间的积分有几种方案，当f=0时，变量的初值出现在方程1-230的右端，从而可直接求出在现时刻的未知变量值，这种方案称为显式时间积分方案。

当0< f<1时，现时刻的未知变量出现在方程的两端，需要解由若干方程组成的方程组才能求出现时刻的变量值，这种方案称为隐式时间积分方案。

• 当f=1时，称为全隐式时间积分方案。
• 当f=1/2时，称为Crank-Nicolson时间积分方案。

进一步将一维问题扩展为二维与三维问题。在二维问题中，计算区域离散见图1-6。发现只是增加第二坐标y，控制体积增加的上下界面，分别用n（north）和s（south）表示，相应的两个邻点记为N和S。在全隐式时间积分方案下的二维瞬态对流-扩散问题的离散方程为

（1-232）

式1-232中

（1-233）

在三维问题中，计算区域离散如图1-7所示（两个方向的投影）。在二维的基础上增加第三坐标z，控制体积增加的前后界面，分别用t（top）和b（bottom）表示，相应的两个邻点记为T和B。在全隐式时间积分方案下的三维瞬态对流扩散问题的离散方程为

（1-234）

式1-234中

（1-235）

图1-7　三维计算区域离散网格在两个方向上的投影

有限体积法常用的离散格式有中心差分格式、一阶迎风格式、混合格式、指数格式、乘方格式、二阶迎风格式、QUICK格式。

各种离散格式对一维、稳态、无源项的对流-扩散问题的通用控制方程（如式1-236所示），均能得到如式1-237所示的形式，对于高阶情况如式1-238所示。

（1-236）

（1-237）

（1-238）

式1-238中，对于一阶情况，，对于二阶情况， ，其中系数和取决于所使用的离散格式（高阶还有和）。各种离散格式系数和的计算公式见表1-5。

表1-5　 不同离散格式下系数a_W和a_E的计算公式

离散格式	系数	系数
中心差分格式
一阶迎风格式
混合格式
指数格式
乘方格式
二阶迎风格式
QUICK格式

1.2.6　CFD常用算法

流场计算的基本过程是在空间上用有限体积法（或其他类似方法）将计算区域离散成许多小的体积单元，在每个体积单元上对离散后的控制方程组进行求解。

其本质是对离散方程进行求解，一般可以分为分离解法（segregated method）和耦合解法（coupled method）两大类，各自又根据实际情况扩展成具体的计算方法。

分离解法不直接求解联立方程组，而是顺序地、逐个地求解各变量代数方程组。分离解法中应用广泛的是压力修正法，其求解基本过程如下。

（1）假定初始压力场。

（2）利用压力场求解动量方程，得到速度场。

（3）利用速度场求解连续方程，使压力场得到修正。

（4）根据需要，求解湍流方程及其他标量方程。

（5）判断当前时间步上的计算是否收敛。若不收敛，返回第二步，迭代计算；若收敛，重复上述步骤，计算下一时间步的物理量。

压力修正法有很多实现方式，其中，压力耦合方程组的半隐式方法（SIMPLE算法）应用最为广泛，也是各种商用CFD软件普遍采纳的算法。

耦合解法同时求解离散方程组，联立求解出各变量（等），其求解过程如下。

（1）假定初始压力和速度等变量，确定离散方程的系数及常数项等。

（2）联立求解连续方程、动量方程、能量方程。

（3）求解湍流方程及其他标量方程。

（4）判断当前时间步上的计算是否收敛。若不收敛，返回第二步，迭代计算；若收敛，重复上述步骤，计算下一时间步的物理量。

1．SIMPLE算法

SIMPLE算法就是求解压力耦合方程的半隐方法（Semi-Implicit Method for Pressure Linked Equations），它是Patankar与Spalding在1972年提出的。

在常规离散方法中，压力梯度项的离散会遇到问题。

如图1-8所示，对P控制体积分后，的贡献为P_w−P_e，如w和e为单元中点，则

（1-239）

图1-8　一维控制体积网格

因此，动量方程将包含相间隔（而非相邻）节点间的压力差。

这样导致求解精度降低，且形成锯齿状的压力场，如图1-9所示。

图1-9　一维锯齿状压力场

这类锯齿状压力场对动量方程而言与均匀场相同（奇偶差），因此，高度不均匀的压力场将被动量方程的特殊离散化当作均匀的压力场处理。

连续方程的离散也会遇到问题，如对于一维问题。

对控制体积分后：

即

则

（跳开了P点）

这样也会导致奇偶差。

对以上出现的离散问题，用交错网格法能较好地解决，如图1-10所示。在此方法中，将速度变量u、v直接设置在P控制体的边界面上，即P控制体边界面上的u、v不再通过主节点上的值求得，而是直接解得。

图1-10　交错网格

SIMPLE算法的基本思想可描述为：对于给定的压力场（它可以是假定的值，或是上一次迭代计算所得到的结果），求解离散形式的动量方程，得出速度场。因为压力场是假定的或不精确的，这样，由此得到的速度场一般不满足连续方程，因此，必须对给定的压力场加以修正。

修正的原则是：与修正后的压力场相对应的速度场能满足这一迭代层次上的连续方程。据此原则，把由动量方程离散形式所规定的压力与速度的关系代入连续方程的离散形式，从而得到压力修正方程，由压力修正方程得出压力修正值，接着根据修正后的压力场，求得新的速度场，然后检查速度场是否收敛，若不收敛，用修正后的压力值作为给定的压力场，开始下一层次的计算；如此反复，直到获得收敛的解。

在上述求解过程中，如何获得压力修正值（即如何构造压力修正方程），以及如何根据压力修正值确定“正确”的速度（即如何构造速度修正方程），是SIMPLE算法的两个关键问题。为此，下面先解决这两个问题，然后给出SIMPLE算法的求解步骤。

（1）速度修正方程。

考察一个直角坐标系下的二维层流稳态问题。设有初始的猜测压力场p*，我们知道，动量方程的离散方程可借助该压力场得以求解，从而求出相应的速度分量 u* 和v*。

根据动量方程的离散方程，有

（1-240）

（1-241）

现在，定义压力修正值为正确的压力场与猜测的压力场之差，有

（1-242）

同样地，定义速度修正值和，以联系正确的速度场(u,v)与猜测的速度场，有

（1-243）

（1-244）

将正确的压力场p代入动量离散方程，得到正确的速度场(u,v)。现在，通过离散的动量方程与式1-240和式1-241，并假定源项b不变，有

（1-245）

（1-246）

引入压力修正值与速度修正值的表达式1-242～式1-244，方程1-245和方程1-246可写成

（1-247）

（1-248）

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

为了简化方程1-247和方程1-248的求解过程，在此，引入如下近似处理：略去方程中与速度修正值相关的和。该近似是SIMPLE算法的重要特征。略去后的影响将在后面介绍的SIMPLEC算法中讨论。于是有

（1-249）

（1-250）

以上两式中

（1-251）

将式1-249和式1-250所描述的速度修正值代入方程1-243和方程1-244，有

（1-252）

（1-253）

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

（1-254）

（1-255）

以上两式中

（1-256）

式1-252～式1-256表明，如果已知压力修正值，便可对猜测的速度场作出相应的速度修正，得到正确的速度场(u,v)。

（2）压力修正方程。

在上面的推导中，只考虑了动量方程，其实，如前所述，速度场还受连续方程的约束。这里暂不讨论瞬态问题。对于稳态问题，连续方程可写为

（1-257）

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

图1-11　离散连续方程的标量控制体积

（1-258）

将正确的速度值，即式1-252～式1-256，代入连续方程的离散方程1-258，有

（1-259）

整理后，得

（1-260）

式1-260可简化为

（1-261）

式1-261中

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

是标量控制体积界面上的密度值，同样需要通过插值得到，这是因为密度是在标量控制体积中的节点（即控制体积的中心）定义和存储的，在标量控制体积界面上不存在可直接引用的值。无论采用何种插值方法，对于交界面所属的两个控制体积，必须采用同样的值。

为了求解方程1-261，还必须对压力修正值的边界条件作出说明。实际上，压力修正方程是动量方程和连续方程的派生物，不是基本方程，故其边界条件也与动量方程的边界条件相联系。

在一般的流场计算中，动量方程的边界条件通常有两类。

第一类，已知边界上的压力（速度未知）。

第二类，已知沿边界法向的速度分量。

若已知边界压力，可在该段边界上令，则该段边界上的压力修正值应为零。这类边界条件类似于热传导问题中已知温度的边界条件。

若已知边界上的法向速度，在设计网格时，最好令控制体积的界面与边界相一致，这样，控制体积界面上的速度为已知。

（3）SIMPLE算法的计算步骤。

至此，已经得出了求解速度分量和压力所需要的所有方程。根据前面介绍的SIMPLE算法的基本思想，可给出SIMPLE算法的计算流程，如图1-12所示。

图1-12　SIMPLE算法流程图

2．其他算法

基于SIMPLE算法的改进算法包括SIMPLEC、SIMPLER和PISO。下面介绍这些改进算法，并对各算法进行对比。

（1）SIMPLER算法。

SIMPLER是英文SIMPLE revised的缩写，是SIMPLE算法的改进版本。它是由SIMPLER算法的创始人之一Patankar完成的。

在SIMPLER算法中，为了确定动量离散方程的系数，一开始就假定了一个速度分布，同时又独立地假定了一个压力分布，两者之间一般是不协调的，从而影响了迭代计算的收敛速度。

实际上，与假定的速度场相协调的压力场是可以通过动量方程求出的，因此不必在初始时刻单独假定一个压力场。

另外，在SIMPLER算法中对压力修正值采用了欠松弛处理，而松弛因子是比较难确定的，因此，速度场的改进与压力场的改进不能同步进行，最终影响收敛速度。于是，Patankar便提出了这样的想法：只用修正速度，压力场的改进则另谋更合适的方法。将上述两方面的思想结合起来，就构成了SIMPLER算法。

在SIMPLER算法中，经过离散后的连续方程1-258用于建立一个压力的离散方程，而不用来建立压力修正方程，从而可直接得到压力，而不需要修正。但是，速度仍需要通过SIMPLE算法中的方程1-252～方程1-256来修正。

将离散后的动量方程式重新改写后，有

（1-262）

（1-263）

在SIMPLER算法中，定义伪速度与为

（1-264）

（1-265）

这样，式1-262和式1-263可以变为

（1-266）

（1-267）

以上两式中的系数d，仍沿用式1-251。同样可写出与的表达式。然后将、、 与 的表达式代入离散后的连续方程1-258，有

（1-268）

整理后，得到离散后的压力方程为

（1-269）

式1-269中

方程1-269中的系数与前面的压力修正方程中的系数相同，差别仅在于源项b。这里的源项b是用伪速度来计算的。因此，离散后的动量方程式，可借助上面得到的压力场来直接求解。这样，可求出速度分量和。SIMPLER算法流程图如图 1-13所示。

在SIMPLER算法中，初始的压力场与速度场是协调的，且由SIMPLER算法算出的压力场不必作欠松弛处理，迭代计算时比较容易得到收敛解。但在SIMPLER的每一层迭代中，要比SIMPLE算法多解一个关于压力的方程组，一个迭代步内的计算量较大。总体而言，SIMPLER算法的计算效率要高于SIMPLE算法。

图1-13　SIMPLER算法流程图

（2）SIMPLEC算法。

SIMPLEC是英文SIMPLE consistent的缩写，意为协调一致的SIMPLE算法。它也是SIMPLE的改进算法之一。

在SIMPLE算法中，为求解的方便，略去了速度修正方程中的项，从而把速度的修正完全归结为由于压差项的直接作用。这一做法虽然不影响收敛解的值，但加重了修正值的负担，使得整个速度场迭代收敛速度降低。

在略去时，犯了一个“不协调一致”的错误。为了能略去而同时又能使方程基本协调，在方程1-247的等号两端同时减去，有

（1-270）

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

（1-271）

式1-271中

（1-272）

类似地，有

（1-273）

式1-273中

（1-274）

将式1-273和式1-274代入式1-252和式1-253，得到修正后的速度计算式：

（1-275）

（1-276）

这就是SIMPLEC算法。SIMPLEC算法与SIMPLE算法的计算步骤相同，只是速度修正方程中的系数项d的计算公式有所区别。

由于SIMPLEC算法没有像SIMPLE算法那样将项忽略，因此，得到的压力修正值一般是比较合适的，因此，在SIMPLEC算法中可不再对进行欠松弛处理。但据作者的经验，适当选取一个稍小于1的a_p对进行欠松弛处理，对加快迭代过程中解的收敛也是有效的。

（3）PISO算法。

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

PISO算法与SIMPLE、SIMPLEC算法的不同之处在于：SIMPLE和SIMPLEC算法是两步算法，即一步预测（图1-12中的步骤1）和一步修正（图1-12中的步骤2和步骤3）；而PISO算法增加了一个修正步，包含一个预测步和两个修正步，在完成了第一步修正得到(u,v,p)后寻求二次改进值，目的是使它们更好地同时满足动量方程和连续方程。PISO算法由于使用了预测—修正—再修正3步，从而可加快单个迭代步中的收敛速度。下面介绍这3个步骤。

① 预测步。

使用与SIMPLE算法相同的方法，利用猜测的压力场p*，求解动量离散方程式1-240与式1-241，得到速度分量u*与v*

② 第一步修正。

所得到的速度场（，）一般不满足连续方程，除非压力场p*是准确的。现引入对SIMPLE的第一个修正步，该修正步给出一个速度场（，），使其满足连续方程。此处的修正公式与SIMPLE算法中的完全一致，只不过考虑到在PISO算法还有第二个修正步，因此，使用不同的记法。

（1-277）

（1-278）

（1-279）

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

（1-280）

（1-281）

就像在SIMPLE算法中一样，将式1-280和式1-281代入连续方程1-258，得到压力修正方程。求解该方程，产生第一个压力修正值。一旦压力修正值已知，可通过式1-280和式1-281获得速度分量与。

③ 第二步修正。

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

（1-282）

（1-283）

再次求解动量方程，可以得到两次修正的速度场。

（1-284）

（1-285）

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

从式1-284中减去式1-282，从式1-285中减去式1-283，有

（1-286）

（1-287）

式1-286和式1-287中，记号是压力的二次修正值。有了该记号，可表示为

（1-288）

将和的表达式代入连续方程1-258，得到二次压力修正方程：

（1-289）

式1-289中，。可写出各系数如下。

（1-290）

（1-291）

（1-292）

（1-293）

（1-294）

现在，求解方程1-289，就可得到二次压力修正值，从而可得到二次修正的压力场。

（1-295）

最后，求解方程1-286与方程1-287，得到二次修正的速度场。

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

图1-14　PISO算法流程图

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

3．SIMPLE系列算法的比较

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

在SIMPLE算法中，压力修正值 能够很好地满足速度修正的要求，但压力修正不是十分理想。改进后的SIMPLER算法只用压力修正值 来修正速度，另外构建一个更加有效的压力方程来产生“正确”的压力场。

由于在推导SIMPLER算法的离散化压力方程时，没有任何项被忽略，因此所得到的压力场与速度场相适应。

在SIMPLER算法中，正确的速度场将导致正确的压力场，而在SIMPLE算法中则不是这样。因此SIMPLER算法是在很高的效率下正确计算压力场的，这一点在求解动量方程时有明显优势。

虽然SIMPLER算法的计算量比SIMPLE算法高出30%左右，但其较快的收敛速度使得计算时间减少30%～50%。

SIMPLEC算法和PISO算法总体上与SIMPLER算法具有同样的计算效率，相互之间很难区分谁高谁低，对于不同类型的问题每种算法都有自己的优势。

一般来讲，动量方程与标量方程（如温度方程）如果不是耦合在一起的，则PISO算法在收敛性方面显得很好，且效率较高。而在动量方程与标量方程耦合非常密切时，SIMPLEC和SIMPLER算法的效果可能更好些。

1.2.7　计算域网格生成技术

网格分布是流动控制方程数值离散的基础，因此，网格生成技术是CFD成功实现数值模拟的关键前提之一，网格质量的好坏直接影响到计算的敛散性及结果的精度。

网格生成的难度和耗费在整个模拟计算程序中占有较大的比重，“从某种角度看，自动生成绕复杂外形的理想网格甚至比编制一个三维的流场解程序更困难，即使在CFD高度发达的国家，网格生成仍占一个计算任务全部人力与时间的60%”，由此可见网格生成在实现流场解算功能过程中的重要性。

网格生成的实质是物理求解域与计算求解域的转换。一般而言，物理域与计算域间的转换应满足下述基本条件。

（1）生成的网格使物理求解域上的计算节点与计算求解域上的计算节点一一对应，不致于出现物理对应关系不确定的多重映射节点。

（2）生成的网格能够准确反映求解域的复杂几何边界形状变化，能够便于边界条件的处理。

（3）物理求解域上的网格应连续光滑求导，保证控制方程离散过程中一阶甚至多阶偏导数的存在性、连续性。网格中出现的尖点、突跃点都将导致算法发散。

（4）网格的疏密易于控制，能够在气动参数变化剧烈的位置，如击波面、壁面等处加密网格，而在气动参数变化平缓的位置拉疏网格。

计算机技术的发展为计算流体力学步入工程实用阶段提供了可能，如何有效地处理复杂的物面边界，生成高质量的计算网格，是目前计算流体力学一个重要的研究课题。

结构化网格在拓扑结构上相当于矩形域内的均匀网格，其节点定义在每一层的网格线上，因此对于复杂外形物体要生成贴体的结构网格是比较困难的。而非结构化网格节点的空间分布完全是随意的，没有任何结构特性，适应性强，因而适合于处理复杂几何外形，并且由于非结构化网格在其生成过程中都要采用一定的准则进行优化判定，因而生成的网格质量较高。

Wimslow最早在20世纪60年代利用有限面积法用三角形网格对Poisson方程进行了数值求解；而20世纪90年代以后，国外学者Dawes、Hah、Prekwas等采用了非结构化网格进行流场的数值求解，国外的著名商业CFD软件，如FLUENT、Star-CD等在20世纪90年代后都将结构化网格计算方法推广到非结构化网格上。

近年来，国内学者也开始对非结构化网格进行深入的研究与探讨，由此可见非结构化网格已成为目前计算流体力学学科中的一个重要方向。

非结构化网格的生成方法有很多，较常用到的有两种：Delaunay三角化方法和推进阵面法（advancing front method）。

（1）Delaunay三角化方法。

Delaunay三角化方法是按一定的方式在控制体内布置节点。定义一个凸多边形外壳，将所有的点都包含进去，并在外壳上进行三角化的初始化。

将节点逐个加入已有的三角化结构中，根据优化准则破坏原有的三角化结构，并建立新的三角化结构，对有关数据结构进行更新后，继续加点，直到所有的节点都加入其中，三角化过程结束。所生成的三角形网格，可以通过光顺技术进一步提高质量。

常用的一种网格光顺方法称为Laplacian 光顺方法，这种光顺方法是通过将节点向这个节点周围的三角形所构成的多边形的形心的移动来实现的。

（2）推进阵面法。

推进阵面法是网格和节点同时生成的一种生成方法，它的基本方法为：根据网格密度控制的需要，在平面上布置一些控制点，给每一个控制点定义一个尺度。根据这些控制点将平面划分成大块的三角形背景网格。

每一个背景网格中的所有点的尺度都可以根据其3个顶点的尺度插值得到。因此相当于布置了一个遍布整个平面的网格尺度函数。根据内边界定义初始阵面，按顺时针方向进行阵面初始化。初始化阵面上的每一条都称为活动边。

由初始阵面上的一条活动边开始推进，根据该活动边的中点所落入的背景网格插值确定该点的尺度，根据该点的尺度及有关的规则确定将要生成的节点位置。

判定该节点是否应被接纳，并根据情况生成新的三角形单元，更新阵面，并沿阵面的方向继续推进生成三角形，直至遇到外边界，网格生成结束。

1.3　CFD软件的构成

CFD软件一般由前处理器、求解器、后处理器3部分组成。这三大模块各有其独特的作用。

1.3.1　前处理器

前处理器（preprocessor）用于完成前处理工作。前处理环节是向CFD软件输　　入所求问题的相关数据，该过程一般是借助与求解器相对应的对话框等图形界面来完成的。流动问题的解是在单元内部的节点上定义的，解的精度由网格中单元的数量所决定。

一般来讲，单元越多，尺寸越小，所得到的解的精度越高，但所需要的计算机内存资源及CPU时间也相应增加。

为了提高计算精度，在物理量梯度较大的区域，以及我们感兴趣的区域，往往要加密计算网格。在前处理阶段生成计算网格的关键是把握好计算精度与计算成本之间的平衡。在前处理阶段需要用户进行以下工作。

（1）定义所求问题的几何计算域。

（2）将计算域划分成多个互不重叠的子区域，形成由单元组成的网格。

（3）对所要研究的物理和化学现象进行抽象，选择相应的控制方程。

（4）定义流体的属性参数。

（5）为计算域边界处的单元指定边界条件。

（6）对于瞬态问题，指定初始条件。

1.3.2　求解器

求解器（solver）的核心是数值求解算法。常用的数值求解方案如1.2.6节所述。总体上讲这些方法的求解过程大致相同，包括以下步骤。

（1）使用简单函数近似待求的流动变量。

（2）将该近似关系代入连续性的控制方程中，形成离散方程组。

（3）求解代数方程组。

各种数值求解方案的主要差别在于流动变量被近似的方式及相应的离散化过程。

1.3.3　后处理器

后处理的目的是有效地观察和分析流动计算结果。随着计算机图形处理功能的提高，目前的CFD软件均配备了后处理器（postprocessor），它提供了较为完善的后处理功能，具体包括以下几方面。

（1）计算域的几何模型及网格显示。

（2）矢量图（如速度矢量线）。

（3）等值线图。

（4）填充型的等值线图（云图）。

（5）XY散点图。

（6）粒子轨迹图。

（7）图像处理功能（平移、缩放、旋转等）。

借助后处理功能，可以动态模拟流动效果，直观地了解CFD的计算结果。

1.4　常用的商业CFD软件

下面介绍近30年来，出现的较为著名的商业CFD软件，包括Phoenics、CFX、STAR-CD和FLUENT等。

1.4.1　Phoenics软件

Phoenics是英国CHAM公司开发的模拟传热、流动、反应、燃烧过程的通用CFD软件，有30多年的历史。网格系统包括直角、圆柱、曲面（包括非正交和运动网格，但在其VR环境不可以）、多重网格、精密网格。

它可以对三维稳态或非稳态的可压缩流或不可压缩流进行模拟，包括非牛顿流、多孔介质中的流动，并且可以考虑黏度、密度、温度变化的影响。

在流体模型上面，Phoenics内置了22种适合于各种Re数场合的湍流模型，包括雷诺应力模型、多流体湍流模型和通量模型及k-e模型的各种变异，共计21个湍流模型、8个多相流模型和十多个差分格式。

Phoenics的VR（虚拟现实）彩色图形界面菜单系统是CFD软件中前处理最方便的一个，可以直接读入Pro/Engineer建立的模型（需转换成STL格式），使复杂几何体的生成更为方便，在边界条件的定义方面也极为简单，并且网格自动生成。但其缺点则是网格比较单一粗糙，不能细分复杂曲面或曲率小的地方的网格，即不能在VR环境里采用贴体网格。

Phoenics的VR后处理也不是很好。要进行更高级的分析则要采用命令格式进行，但这在易用性上比其他软件就要差了。

另外，Phoenics自带了1 000多个例题与验证题，附有完整的可读可改的输入文件。其中就有CHAM公司做的一个PDC钻头的流场分析。Phoenics的开放性很好，提供对软件现有模型进行修改、增加新模型的功能和接口，可以用FORTRAN语言进行二次开发。

1.4.2　STAR-CD软件

STAR-CD的创始人之一Gosman与Phoenics的创始人Spalding都是英国伦敦大学同一教研室的教授。STAR是Simulation of Turbulent flow in Arbitrary Region的缩写，CD是Computational Dynamics Ltd的缩写。

STAR-CD是基于有限容积法的通用流体计算软件，在网格生成方面，采用非结构化网格，单元体可为六面体、四面体、三角形界面的棱柱、金字塔形的锥体以及6种形状的多面体，还可与CAD、CAE软件接口，如ANSYS、IDEAS、NASTRAN、PATRAN、ICEMCFD、GRIDGEN等，这是STAR-CD在适应复杂区域方面的特别优势。

STAR-CD能处理移动网格，用于多级透平的计算，在差分格式方面，纳入了一阶迎风、二阶迎风、CDS、QUICK以及一阶迎风与CDS或QUICK的混合格式。

在压力耦合方面采用SIMPLE、PISO以及称为SIMPLO的算法。

在湍流模型方面，有k-e、RNK-ke、ke两层等模型，可计算稳态、非稳态、牛顿流体、非牛顿流体、多孔介质、亚音速、超音速和多项流等问题。STAR-CD的强项在于汽车工业、汽车发动机内的流动和传热。

1.4.3　ANSYS CFX软件

ANSYS CFX系列软件是拥有世界级先进算法的成熟商业流体计算软件。功能强大的前处理器、求解器和后处理模块使得ANSYS CFX系列软件的应用范围遍及航空、航天、船舶、能源、石油化工、机械制造、汽车、生物技术、水处理、火灾安全、冶金、环保等众多领域。

ANSYS CFX软件系列包括CFX Solver、BladeModeler、TurboPre、TurboPost、TurboGrid等软件，如图1-15所示。

图1-15　ANSYS CFX系列软件

DeSignModeler基于ANSYS的公共CAE平台—Workbench，提供了完全参数化的几何生成、几何修正、几何简化，以及概念模型的创建能力，它与Workbench支持的所有CAD软件能够直接双向关联。

CFX软件提供了从网格到流体计算以及后处理的整体解决方案。核心模块包括CFXMesh、CFXPre、CFXSolver和CFXPost 4个部分。其中CFXSolver是CFX软件的求解器，是CFX软件的内核，它的先进性和精确性主要体现在以下方面。

（1）不同于大多数CFD软件，CFXSolver采用基于有限元的有限体积法，在保证有限体积法守恒特性的基础上，吸收了有限元法的数值精确性。

（2）CFXSolver采用先进的全隐式耦合多网格线性求解，再加上自适应多网格技术，同等条件下比其他流体软件快1～2个数量级。

（3）CFXSolver支持真实流体、燃烧、化学反应和多相流等复杂的物理模型，这使得CFX软件在航空工业、化学及过程工业领域有着非常广泛的应用。

ANSYS CFX特别为旋转机械定制了完整的软件体系，向用户提供从设计到CFD分析的一体化解决方案，因此，ANSYS CFX被公认为是全球最好的旋转机械工程CFD软件，被旋转机械领域80%以上的企业选作动力分析和设计工具，包括GE、Pratt & Whitney、Rolls-Royce、Westinghouse、ABB、Siemens等企业界巨头。

ANSYS CFX包含的专用旋转机械设计分析工具有BladeModeler、TurboGrid、TurboPre和TurboPost。

BladeModeler：是交互式涡轮机械叶片设计软件，用户通过修改元件库参数或完全依靠BladeModeler中提供的工具设计各种旋转和静止叶片元件及新型叶片。软件简单实用，模块丰富，具有自动化程度高和叶片几何生成迅速的特点。

TurboGrid：是专业的涡轮叶栅通道网格划分软件，用户只需提供叶片数目、叶片及轮毂和外罩的外形数据文件。自动化程度高，网格生成迅速，生成网格质量高是它的优点。

TurboPre：包含于CFXPre中，是专业的旋转机械物理模型设置模块，以旋转机械的专业术语完成模型设置。

TurboPost：包含于CFXPost中，是专用的旋转机械问题模拟结果后处理模块，可以自动生成子午面等专业视图，同时提供效率、压比和扭矩等旋转机械性能参数。

1.4.4　ANSYS FLUENT软件

FLUENT自1983年问世以来，就一直是CFD软件技术的领先者，被广泛应用于航空航天、旋转机械、航海、石油化工、汽车、能源、计算机/电子、材料、冶金、生物、医药等领域，这使FLUENT公司成为占有最大市场份额的CFD软件供应商。作为通用的CFD软件，FLUENT可用于模拟从不可压缩到高度可压缩范围内的复杂流动。

其代表性客户包括美国国家航空航天局（NASA）、美国国防部（DOD）、美国能源部（DOE）等政府部门以及BMW-RR、波音、福特、GE、三菱等企业。

关于ANSYS FLUENT软件的相关介绍及详细应用，请查阅后面章节。

1.5　本章小结

为了方便初学者迅速进入学习计算流体动力学的大门，本章首先介绍流体力学中支配流体流动的基本物理定律，包括一些基本假设、概念及流体力学基本方程组等。

在此基础上介绍了用数值方法求解流体力学问题的基本思想，进而阐述计算流体力学（CFD）的相关基础知识，包括CFD数值模拟方法、空间和方程的离散方法及CFD常用算法等。最后简要介绍常用计算流体力学商业软件。学好本章能为读者更好地学习计算流体力学以及掌握CFD软件打下坚实的基础。.