计算流体力学的基本思想

三 、求解控制方程的数值方法
（4） （5a） （5b） （5c） （6）
对于无法用解析方法求解的微分方程可以用数值方法求解， 所谓数值方法求解就是用近 似的数值解逼近微分方程的精确解。流动控制方程的精确解是流场计算域内流动参数（如速 度、压力、温度等）的连续分布，而数值解则是流场计算域内离散的点上的近似解对连续精 确解的逼近，换句话说，我们可以把连续的流场离散为一定数目的不连续的点，在这些离散 点上，守恒方程被近似满足，如果离散点之间的距离为无穷小，则近似解将无限趋近于精确 解，因此我们可以用近似解代替精确解。这就是流动微分方程数值求解的基本思想。 以数值方法求解流动微分方程，首先要把需要求解的流场的几何空间（或称为计算域） 离散为孤立的不连续的点，或者说用一定数量的点覆盖或代表要求解的连续的流场，然后将 流动控制方程的偏导数用离散点之间的有限变化来代替， 例如， 表示速度梯度的导数 ∂u / ∂x 用差商 Δu / Δx 来代替，其中 Δu 和 Δx 分别是 x 坐标方向的两个相邻的点的速度差和坐标 x 的增量。 可以想象， 如果控制微分方程中的所有导数或偏导数都被类似于差商的量代替的话， 偏微分方程将有可能变成一个线性方程，一个只包含离散点的坐标和待求函数值（如上述的 u）的线性方程。事实上，我们可以把流动控制方程组的每一个偏微分方程在每一个离散点 上转变为一个线性方程。假如我们用 100 个点离散一个计算域，那么对每个偏微分方程我们 将得到 100 个线性方程。至此，偏微分方程的求解已经转化为线性方程组的求解，如果得到 线性方程组的解，我们就得到了偏微分方程组的近似数值解。因此，我们也可以说，CFD 模 拟的过程本质上是在计算域上构建线性方程组并求解线性方程组的过程。 从上面的论述可以看出，数值方法求解流动微分方程至少包括三个步骤：首先，离散计 算域；其次，在离散后的计算域上离散控制方程；其三，求解离散得到的线性方程组。需要 补充的是，并不是所有的线性方程都需要求解，实际上有些特殊点上的流动变量值或其梯度 是已知的，这些特殊的点就是计算域边界上的点。通常为了限定微分方程的解，我们需要给
CFD 数值模拟实验指导书
第一部分
计算流体力学（CFD）的基本思想
一、什么是计算流体力学（CFD）？ 计算流体力学（Computational Fluid Dynamics）是流体力学的一个新兴的分支，是一 个采用数值方法利用计算机来求解流体流动的控制偏微分方程组， 并通过得到的流场和其它 物理场来研究流体流动现象以及相关的物理或化学过程的学科。事实上，研究流动现象就是 研究流动参数如速度、压力、温度等的空间分布和时间变化，而流动现象是由一些基本的守 恒方程（质量、动量、能量等）控制的，因此，通过求解这些流动控制方程，我们就可以得 到流动参数在流场中的分布以及随时间的变化，这听起来似乎十分简单。但遗憾的是，常见 的流动控制方程如纳维-斯托克斯（Navier-Stokes）方程或欧拉（Euler）方程都是复杂的 非线性的偏微分方程组，以解析方法求解在大多数情况下是不可能的。实际上，对于绝大多 数有实际意义的流动，其控制方程的求解通常都只能采用数值方法的求解。因此，采用 CFD 方法在计算机上模拟流体流动现象本质上是流动控制方程（多数情况下是纳维-斯托克斯方 程或欧拉方程）的数值求解，而 CFD 软件本质上就是一些求解流动控制方程的计算机程序。 二、计算流体力学的控制方程 计算流体力学的控制方程就是流体流动的质量、动量和能量守恒方程。守恒方程的常见 的推导方法是基于流体微元的质量、动量和能量衡算。通过质量衡算可以得到连续性方程， 通过动量守恒可以得到动量方程，通过能量衡算可以得到能量方程。式（1）-（3）是未经 任何简化的流动守恒微分方程，即纳维-斯托克斯方程（N-S 方程） 。
N-S 方程可以表示成许多不同形式，上面的 N-S 方程是所谓的守恒形式，之所以称为守恒形 式，是因为这种形式的 N-S 方程求解的变量 ρ 、 ρu 、 ρv 、 ρw 、 ρE 是守恒型的，是质量、
Fra Baidu bibliotek
1
CFD 数值模拟实验指导书
动量和能量的守恒变量。事实上也可以直接求解 u、v、w、T 等原始变量，这种形式的 N-S 方程被称为非守恒形式，因为这些变量并不守恒。也可以根据具体的流动状况进行简化，比 如对于无粘流动 N-S 方程可以简化为欧拉方程（粘性项被去掉） ，如式（4）-(6) 所示；对 于不可压缩流动（液体的流动，马赫数小于 0.3 的气体流动） ，N-S 方程可以简化为不可压 缩的 N-S 方程（密度恒定，因此被消去） ；对于定常流动，N-S 方程可以去掉时间导数项， 简化为稳态的 N-S 方程；流体流动往往具有三维性质，但是也常常可以简化为二维流动甚至 一维流动。对于 CFD 的计算来说三维简化为二维或一维意味着运算量的大幅度降低。
2
CFD 数值模拟实验指导书
出定解条件，在这里就是所谓边界条件。同样的道理，对于包含时间导数的微分方程，我们 需要给定初始条件。 上面我们用差商取代导数的方法介绍了离散 （把连续空间里的微分方程转化为该连续空 间内的不连续的点上的近似的线性方程的过程叫做离散化）微分方程的思想。但是应该注意 的是，流动控制微分方程的离散化需要严谨的数学推导、证明和分析。离散化方法的研究是 CFD 最重要的部分，也是 CFD 中的数值方法的基础。计算流体力学中有三大类主要离散化方 法，即：有限差分方法（FDM） ，有限体积方法（FVM）和有限元方法（FEM） 。三者的区别主 要在于它们处理最基本的离散单元的方法，其中有限差分和有限体积法更为常用。有限差分 法通常在离散点上直接以差分替代微分（即差商替代导数） ，差分可以分为向前、向后和中 心差分；有限体积法则首先对构造在离散点周围的控制体进行积分，将一阶导数项转换为代 数项，然后在控制体界面插值来实现离散化。对于不同的控制方程，每一类方法又有许多具 体的实施办法，这些实施方法被称为格式（scheme） 。 1. 计算网格的生成 在计算流体力学术语中，计算域的离散被称之计算网格生成，所谓网格实际上就是用上 述的离散点以某种方式连接而成的“网络”。最直观的网格是二维网格，例如，我们可以将 一个矩形计算域用一定间隔的 x 方向的若干条直线和类似的 y 方向的若干条直线划分为一个 个小的矩形单元组成的网状结构，这个网状结构就是一个最简单的二维网格。前述的用于离 散控制方程的点可以是网格线的交叉点，也可以是矩形单元的中心，这取决于离散控制方程 所采用的方法。 实际上，划分网格有很多方法，网格线可以是直线或曲线、正交的或非正交的，网格线 的间隔可以是均匀的或非均匀的。而有些网格并不存在有意义的网格线，或者说网格线没有 规则的结构，如用小的三角形单元构成的二维网格（类似于有限元网格） ，这样的网格被称 为非结构网格（unstructured grid） ，相对应的是前面所说的具有直线或曲线网格线的网格 被称为结构网格（structured grid） 。二维网格是最据直观意义的网格，而一维网格的划分 实际上是将一个有限长度的直线或曲线分割成长度一定数量的均匀或不均匀的小的线段， 控 制方程将在这些小线段的端点或中心离散。 三维网格则可以看作二维网格在第三维方向的延 伸，例如三维结构网格的网格单元常见的是长方体或扭曲的长方体（视直线网格或曲线网格 而定） ，三维非结构网格的网格单元多为四面体。 网格生成是 CFD 模拟的一个十分重要的部分，为了确保计算精度，网格必须足够密集， 事实上我们并不要求网格的密度在整个网格范围均匀一致， 通常对流动参数梯度大的地方要 采用较为密集的网格（例如激波的位置，边界层附近） ，梯度小的地方则可以适当采用疏松 的网格（比较开阔的空间、流动被扰动较少的地方） 。一个高质量的网格是 CFD 模拟成功的 关键因素，不合适的网格可能直接导致计算的失败。因此，人们在生成网格上花费的时间常 常超过全部 CFD 工作时间的 50％，对于复杂的几何形状网格生成所花费的时间甚至达到 70 ％。由于网格生成的复杂性和巨大的工作量，许多专业的网格生成工具应运而生，例如 ICEM CFD, Gambit 等。 2. 边界条件与初始条件 对于 CFD 模拟要求解的问题，计算域的几何边界定义了流场的范围，或者说计算域是由
3
CFD 数值模拟实验指导书
几何边界确定的，而边界的物理特性则定义了问题本身。如前所述，边界点的流动参数值常 常是给定的，因此是已知的，这就是边界条件。从给定方式来看，边界条件有三种形式：其 一， Dirichlet 边界条件， 直接给定流动参数的值， 如给定边界的速度、 温度； 其二， Neumann 边界条件，给定一阶导数，如给定压力梯度；其三，混合边界条件，是 Dirichlet 边界条件 和 Neumann 边界条件的混合。以上三种边界条件也被称为第一、第二、第三类边界条件。从 边界的物理性质来看，边界条件又可分为：固壁边界条件、入口边界条件、出口边界条件等 等。给定正确或合适的边界条件对于 CFD 计算也是十分重要的，实际上流场的特性很大程度 上是由边界条件决定的。 相对而言，初始条件的设定比较简单，我们需要给定的是一个初始时刻已知的流场。事 实上并不是所有的 CFD 计算都需要初始条件， 初始条件仅对于随时间变化的流场的求解才是 必不可少的。 3. 计算结果的后处理 一个成功的 CFD 计算环节完成之后，CFD 程序或软件将计算结果写入一个或多个特定格 式（因特定的软件而异）的数据文件，这些数据文件通常包括计算网格点的坐标，每个网格 点上的流动参数值（如速度，压力，温度，密度等） ，对于这些数据的分析还需要专门的工 具软件，这些工具软件将网格的结构、流动参数的分布等显示出来。常见的基本的显示方法 包括标量等值线分布（如温度、压力分布的云图） ，向量分布（例如用带箭头的线段表示速 度的大小和方向） ，X-Y 曲线图等。这些后处理方法将计算结果清晰地显示出来，供人们方 便地分析和评价计算结果。 4. 举例：一维激波管内流动的 CFD 模拟 一维激波管内压力驱动的气体可压缩流动是一个气体动力学的经典问题。如图 1 所示， 考虑一个两端封闭、中间用隔膜（Diaphragm）隔开的管子，隔膜左右的封闭空间里分别充 满不同压力的气体，其中左面为高压 p 4 ，右面为低压 p1 , 有 p 4 > p1 ，而且压差较大。假 设在某一瞬间隔膜破裂，隔膜右侧的空间产生一个向右运动的激波，而在隔膜的左面将产生 一个向左运动的膨胀波，该激波和膨胀波随时间的推移向两个不同方向传播，如图 2 所示。 实际上，在隔膜破裂、可压缩流动建立以后， 激波管内存在着四个不同的区域， 如图 2 所示， 从右到左，区 1 是激波前面的尚未受到扰动的低压区域（压力为初始压力 p1 ，密度为初始 ；区 2 是位于激波之后，是被激波扰动过的低压区域，在这个区中，压力和密度均 密度 ρ1 ） 大于区 1 中的压力和密度，即： p 2 > p1 ， ρ 2 > ρ 1 ；区 3 位于膨胀波之后，是已经受到膨 胀波扰动的高压区，其压力和密度均低于初始的压力和密度，即： p3 < p 4 ， ρ 3 < 区 2 相比较，尽管 ρ 3 >
（1） （2a）
（2b）
v ∂ ( ρw) ∂p ∂τ xz ∂τ yz ∂τ zz （2c） + ∇ ⋅ ( ρwV ) = − + + + + ρf z ∂t ∂z ∂x ∂y ∂z v ∂ ( ρE ) ∂ ∂T ∂ ∂T ∂ ∂T & + (k + ∇ ⋅ ( ρEV ) = ρq ) + (k ) + (k ) ∂t ∂x ∂x ∂y ∂y ∂z ∂z ∂ (up) ∂ (vp) ∂ ( wp ) ∂ ∂ ∂ − − − + (uτ xx ) + (uτ yx ) + (uτ zx ) ∂x ∂y ∂z ∂x ∂y ∂z （3） ∂ ∂ ∂ ∂ ∂ + (vτ xy ) + (vτ yy ) + (vτ zy ) + ( wτ xz ) + ( wτ yz ) ∂x ∂y ∂z ∂x ∂y v v ∂ + ( wτ zz ) + ρf ⋅ V ∂z
v ∂ρ + ∇ ⋅ ( ρV ) = 0 ∂t v ∂ ( ρu ) ∂p + ∇ ⋅ ( ρuV ) = − + ρf x ∂t ∂x v ∂ ( ρv) ∂p + ∇ ⋅ ( ρvV ) = − + ρf y ∂t ∂y v ∂ ( ρw) ∂p + ∇ ⋅ ( ρwV ) = − + ρf z ∂t ∂z v v v ∂ ( ρE ) ∂ (up) ∂ (vp) ∂ ( wp ) &− + ∇ ⋅ ( ρEV ) = ρq − − + ρf ⋅ V ∂t ∂x ∂y ∂z
v ∂ρ + ∇ ⋅ ( ρV ) = 0 ∂t v ∂ ( ρu ) ∂p ∂τ xx ∂τ xy ∂τ zx + ∇ ⋅ ( ρuV ) = − + + + + ρf x ∂t ∂x ∂x ∂y ∂z
v ∂ ( ρv ) ∂p ∂τ xy ∂τ yy ∂τ yz + ∇ ⋅ ( ρvV ) = − + + + + ρf y ∂t ∂y ∂x ∂y ∂z