计算流体力学

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