计算流体力学基础

1、计算流体力学基础目录u 一、CFD概述u 二、流体力学的相关基本概念u 三、描述流体流动的基本控制方程u 四、实现CFD的基本步骤u 五、有限容积法离散控制方程的基本步骤u 六、CFD实现中的关键技术u 七、CFD的初始和边界条件u 八、商用CFD软件介绍一、CFD概述 什么是CFD？ 为什么要使用CFD？ CFD的应用发展简介 CFD的典型应用领域什么是CFD？ CFD-Computational Fluid Dynamics：计算流体动力学，是一种以流体力学和传热学为学科背景，基于计算机技术发展起来的模拟分析工具； 本质是采用计算机技术来获得流体力学/传热学等问题的数值解； 主要用来模拟流

2、体流动、换热以及可能伴随的化学反应等物理、化学过程； 现代CFD的商用软件已经扩展了CFD的求解功能， 不但能求解流固耦合的换热问题，还可进行流 固耦合的结构分析，甚至多物理场的分析等。为什么要使用CFD？不管对于科学研究，还是工业级别的应用，CFD都是一个很强大的工具。相对纯粹的实验手段，CFD具有以下优势： 缩短新产品的开发周期，降低开发成本； 对于实验难以进行，甚至根本不能进行的情况进行分析； 对工作于危险工况下的情况进行分析； 可以得到比实验方法丰富得多的系统细节信息； 可以很容易地实现参数化的分析，从而达到系统优化的目的.CFD的应用发展简介 CFD起于上世纪60年代。在飞行器和喷气

3、发动机的设计中使用了CFD技术； 直到80年代初，CFD的使用主要在高校和研究所等科研部门，使用者需具有相当的流体力学和数值分析的背景； 自1981年CHAM公司推出第一款商用CFD软件PHOENICS之后，CFD技术在工业领域的应用得到大 大推动。商用CFD软件对使用者的知识背景要求有所降低； 80年代之后的二十多年是商用CFD软件大力发展的阶段。现在出现的新一代的商用CFD软件，其自动化程度大大提高，对使用者的知识背景要求大大降低，甚至可较容易地被非CFD相关专业的工程师掌握，在工业领域具有广阔的应用前景。CFD的典型应用领域(1) 飞行器和车辆等的空气动力学分析 船舶行业的水力学分析 内

4、燃机和气轮机的燃烧过程分析 车辆发动机舱内的流动换热分析 车辆/飞机/房间等的舒适性分析 电子工程中的器件冷却分析 化工过程中的混合、分离及聚合物的形成分析CFD的典型应用领域(2) 建筑物的内部和外部环境分析 环境工程中污染物的扩散分析 海事工程中近海建筑物的负载分析 河流、港湾和海洋流动分析 透平机械内部流动分析 天气预测 生物医学工程领域的血管内血液流动分析 原子能工程的冷却(沸腾)换热分析二、流体力学的相关基本概念 定常流动和非定常流动 可压缩流体和不可压缩流体 粘性流体和无粘性流体 牛顿流体和非牛顿流体 层流和湍流定常流动和非定常流动 定常流动：也称为稳态流动，在流体流动过程中，所有

5、物理量的变化均和时间无关，而只是空间位置的函数； 非定常流动：也称为瞬态流动，在流体流动过程中，某个或某些物理量随时间的变化而变化； 在CFD中，非定常流动的模拟要耗费比定常流动多得多的计算资源。可压缩流体和不可压缩流体 流体流动过程中压力的变化可引起密度的相应变化； 如果在流动研究中需要考虑流体的这种密度变化，则称为可压缩流体，反之则称为不可压缩流体； 高速流动的情况，如超音速流和马赫数较高的亚音速流动等，一般按可压缩流体处理； 液体的流动，以及绝大多数的低速流动一般按不可压缩流体处理。粘性流体和无粘性流体 粘性是指流体在运动状态下所产生的抵抗剪切变形的能力； 如果在研究流体流动时考虑流体的

6、粘性，则流体称为粘性流体，否则称为无粘性流体， 也叫理想流体。无粘性流体粘性流体牛顿流体和非牛顿流体 在研究粘性流体的流动时，如果流体遵守牛顿内摩擦定律，即则此种流体称为牛顿流体，反之则称为非牛顿流体；依据剪切应力跟剪切变形速率的关系，非牛顿流体又可分为塑性流体、假塑性流体和胀塑性流体（不在本次讨论范围之内，具体可参考相关书籍）。t = m du (t为剪切应力，m为动力粘度，du 为流体的剪切变形速率 ) dydy层流和湍流 层流流动的特点是流体层与层之间无任何干扰（既无质量传递也无动量传递），流体微团互不掺混，从现象上来看，流动轨迹有条不紊； 层流流动的雷诺数超过临界雷诺数时，就会向湍流发

7、展； 层流和湍流两种形态的流动在物理现象上有很大 的区别。湍流流体内充满了漩涡，做不规则的， 混乱方式运动。其漩涡快速流经某一点时，局部 参数值的大小和方向随时间和空间变化特别频繁； 随机性是湍流的主要特性。湍流运动中的任何物理量随时间和空间随机地变化着。层流湍流三、描述流体流动的基本控制方程 散度和梯度 基本控制方程 控制方程的通用形式 通用控制方程的物理意义散度和梯度U =u,v速度的散度：速度的梯度：v2u2yu2u2d yd yxd xd xdivU = u + v u2 - u1 + v2 - v1j u2 - u1 i + u2 -u1 jgradu = u i + ud xd y

8、d xd yxyxy注意：散度的对象矢量，结果是标量；梯度的对象是标量，结果是矢量。vu1uPu1vu1uPv1基本控制方程质量守恒：动量守恒：能量守恒：组分守恒： sm _ j = smjYj =1j(rYj ) + div(rUY ) -div(G gradY ) = stjjjm _ j(rT ) + div(rUT ) -div( l gradT ) = sf + shtcpcp(rv) +div(rUv) -div(m grad v) = - p + s+ stya _ vv(rw) +div(rUw) -div(m grad w) = - p + s+ stza _ ww(ru)

9、+div(rUu) -div(m gradu) = - p + s+ stxa _uur + div(rU ) = stm控制方程的通用形式非稳态项对流项扩散项广义源项(rf)t+div(rUf)-div(Gf gradf)=sf通用控制方程的物理意义Ff2Ff 2Ff控制体积内随时间的变化率控制体积内 的净生成率流入控制体积的的净通量(rf) = -div+ stf (rf)dvt dV=- div Ff dvdV+ sf dvdVFf1dVFf1(rUf -Gf gradf)四、实现CFD的基本步骤 定义求解几何区域并进行离散 确定求解参数和边界条件 控制方程的离散化 求解代数方程组定义求

10、解几何区域并进行离散(1) 根据所要求解的具体问题，确定所要求解的几何区域； 根据所要选取的方程离散方法和计算区域几何形状进行计算区域的离散； 计算区域的离散，就是用若干离散的点来代替连续的空间。结果就是把整个求解区域划分成许多不重叠的子区域。控制方程就根据这些划分的子区域进行离散； 计算区域离散之后的网格分为两类：结构化网格和非结构化网格；定义求解几何区域并进行离散(2) 结构化网格结构化网格是指网格区域内所有的内部点都具有相同的毗邻单元数，给出一个节点的编号后，可立即得出其相邻节点的编号；结构化网格的生成技术比较成熟，它实现容易、生成速度快、网格质量好、数据结构简单，而且对曲面或空间的拟合

11、大多数采用参数化或样条插值的方法得到，区域光滑，与实际的模型更容易接近；结构化网格的典型缺点是适用范围窄，很难实现复杂边界区域的离散。因此不适用现在越来越复杂的几何模型需求；比较典型的结构化网格有二维的四边形网格和三维的六面体网格。定义求解几何区域并进行离散(3) 非结构化网格非结构化网格的内部节点分布不规则，各节点周围的网格数目不尽相同；非结构化网格的生成过程比较复杂，但却有极大的适应性，对复杂边界的流场计算问题特别有效；关于网格生成技术，请参考相关的资料。确定求解参数和边界条件 物性参数的确定； 初始条件的确定； 边界条件的确定（具体后面将会涉及）； 方程离散和求解等相关参数的确定。控制方

12、程的离散化(1) 流动问题的数学描述，归结为一系列的偏微分方程。除了一些简单的流动外，绝大多数的偏微分方程是没有理论解的，但我们可以通过数值计算的手段来得到这些流动问题的数值解。而方程离散就是获得数值解的最重要的一环； CFD偏微分方程离散的方法，主要包括有限差 分法、有限容积法和有限元法。其中有限元法 主要用于结构分析，在求解流动问题时，在某 些技术的处理上还不是很成熟，在此不做介绍。下面简单介绍一下有限差分法和有限容积法：控制方程的离散化(2) 有限差分法： 这是最早出现的一种离散方法。基本实现过程是将求解区域划分成一系列平行于坐标轴的节点，将偏微分方程在每个节点上进行离散，即用差分表达式

13、代替偏微分方程中的各阶导数，这样就在每个节点形成了一个代数方程式，而每个代数方程式又包含了本节点以及周围节点的相关信息。将这些代数方程式联立进行求解，就可得到整个流场的数值解。i-1ii+1xDxDx一阶导数二阶导数中心差分：向前差分：向后差分：中心差分： f fi+1 -fi-1 x 2Dxi f fi -fi-1 x Dxi f fi+1 -fi x Dxi 2f f- 2f +fi+1ii-1 x2 Dx2i控制方程的离散化(3) 有限容积法将整个求解区域划分成一系列的控制容积，每个控制容积内有一个节点，通过将控制方程在控制容积内进行积分而导出离散方程；采用有限容积法进行方程离散可以保证

14、离散方程的守恒性，这对于流动问题的求解是很重要的；目前主流的商用CFD软件都是采用的有限容积法。 关于有限容积法的实施方法，将在下面详细论述求解代数方程组 控制方程完成离散化之后，会生成若干代数方程组，要得到流场的数值解，必须对这些代数方程组进行求解； CFD生成的代数方程组，求解的未知量数目非常庞大（取决于网格数目），而且常常是非线性的，因此通常采用迭代的方式来求解； 关于代数方程的迭代求解方法，在下面的内容中将会进一步进行讨论。五、有限容积法离散控制方程的基本步骤 求解参数的空间型线和时间型线 控制方程在控制容积内的积分 控制方程的离散求解参数的空间型线和时间型线求解参数的两种空间型线：f

15、EfEfffPfPfWfWxxWPEWPE求解参数的两种时间型线：ffft +Dtft +Dt隐式f t显式f txxt + Dtt + Dttt控制方程在控制容积内的积分以一维方程的离散为例，介绍有限容积法的实现。图中阴影部分为控制容积P(我们假设网格为均分网格)：DxWExwe(d x)e(dw一维控制方程的通用形式为：将方程对上图所示的控制容积P在时间间隔Dt内进行积分可得(假设流体不可压)：e r(ft+Dt -ft )dx + t+Dt r(uf) -(uf) dt = t+Dt G f - f t +t+Dte s dxdtwtewtf x x dtw few (rf) + (ru

16、f) = G f + stxx f x fPx)控制方程的离散(1)1. 非稳态项-选定随x变化的型线为阶梯式，即控制容积P内各点值均相同，则非稳态项可离散成：2. 对流项-进行对流项的离散时，需要选定随时间t变化的型线。有三种可选的方式，即显式、隐式和C-N格式(半隐方式)。不同的选择方式会产生不同的离散方程。这里我们选择显式格式。则得对 流项的离散形式为：t+Dttr(uf)e -(uf)wdt = r(uf)e - (uf)w Dttte r(ft+Dt -ft )dx = r(f t +Dt -f t )DxwPP控制方程的离散(2)3. 扩散项-对于扩散项的离散，选取的一阶导数随时间

17、t变化的型线为显式格式，则得：4. 控制容积界面的处理-由于离散方程是针对控制容积的中心节点来求解的，因此在控制容积界面处的参数(对流项中界面处的通量 和扩散项中的一阶导数)需要确定其表达方式。在此我们选择随x的变化型线为分段线性。则：(uf)= (uf)W + (uf)P w2 f = fP -fW x (d x)ww(uf) = (uf)P + (uf)E e2 f fE -fP x = (d x)eet+Dt G f - f = f t f t tf x x dtGf x - x Dtew ew 控制方程的离散(3)5. 源项的处理-假设源项对t和x均呈阶梯状的变化，则源项可离散成(对于

18、源项是求解参数函数的情况，需要特殊处理，请参考相关文献)：6. 最终的离散方程-根据我们前面所做的网格均分的假设，Dx = d x，则将上面各项离散后的结果代入控制容积积分后的方程可得一维模 型通用方程的离散形式为：r ft+Dt -ft + r (uf)t -(uf)t= G ft - 2ft +ft + stPP EW EPW Dt2DxfDx2t+Dtetw sf dxdt = sf DxDtt六、CFD实现中的关键技术 对流项的差分格式 N-S方程中速度和压力的耦合算法 代数方程组的求解方法对流项的差分格式(1) 对流项强烈的方向性，导致其离散具有很大的难度Dx如何确定WEfe、fw？

19、xwe(d x)e(dw对流项的离散格式，就是指如何利用相邻节点的信息来获得 fe及的插fw 值方式e f dx = f - fw xewPx)对流项的差分格式(2) 对流项差分格式的影响：数值解的准确性。扩散项二阶截差的离散格式， 能很好反映扩散过程的特点，可以满足大多数实际问题的准确度要求。数值解误差的主要来源， 在于对流项的离散格式；数值解的稳定性。如对流项采用中心差分格式， 在某些情况下会引起数值解的振荡，即不稳定性；数值解的经济性。很多针对对流项发展的高阶差分格式，虽然能满足数值解的准确性和稳定性的要求，但是r 由d x 于其构造过程比较复杂，使得代数=解u占用较多的内存和时间。方程

20、的P求DG 网格Pe数对流项的差分格式(3)一些基本差分格式的简单介绍迎风差分格式(UD)：对控制容积界面上变量，永远取其上游节点的值。其优点是充分考虑了流动方向的影响，而且是绝对稳定的格式，不会带来解的振荡；缺点是其具有一阶截差， 误差较大；中心差分格式(CD)：采用控制容积界面上下游两个节点值的线性插值来计算界面上的变量值。此种格式具有二阶截差， 在解不发生振荡的情况下，其精度高于迎风格式，但由于其没有充分考虑流动的方向性，在某些情况下会导致解的振荡， 甚至得出无物理意义的解；混合格式：综合考虑了中心差分和迎风格式的特点，根据网格Pe数的大小来判断如何使用差分格式；QUICK(Quadra

21、tic Upstream Interpolation of Convective Kinetics，对流项的二次迎风插值 )格式：对于界面值采用三点（界面上游两点和下游一点）二次插值来计算，具有三阶精度 。可减小假扩散现象，精度比中心差分和混合格式高， 并且具有迎风权重的特点。但有时会导致在k-湍流模型的计算中出现负的湍动能。N-S方程中速度和压力的耦合算法(1) N-S方程求解的问题在N-S方程的广义源项中包含了压力梯度项。在不可压缩流体的求解中，压力梯度项的处理成了最主要的问题；采用常规网格及中心差分法离散压力梯度项时， 动量方程的离散形式可能会得到不合理的压力场(波形压力场)，导致速度和

22、压力失耦；根本原因在于采用中心差分时，压力梯度的离散不是以相邻两点的压差来表示的，而且压力梯度项包含在广义源项中，不存在独立的控制方程。 解决办法 压力梯度的离散以相邻两点的压差来表示； 将压力的求解同连续性方程关联，使压力值在求解过程中得到合理的更新N-S方程中速度和压力的耦合算法(2) 采用交叉网格的SIMPLE算法什么是交叉网格(staggered grid)？-将每个速度分量(u, v, w)及压力p分别存储于不同的网格系统中。以二维为例：NWEWEWEu的制容积v的制容积主控制容积（p及其它标量）NvnuwuePvsSNvnuwPvsueSvnuwuePvsSN-S方程中速度和压力的

23、耦合算法(3)SIMPLE(Semi-implicit Method for PressurelinkedEquations，求解压力耦合方程的半隐方法)由Patankar和Spalding在1972年提出，其基本思想如下：假定一个初始的速度分布和压力分布，并根据此假定求解动量方程，得出新的速度场；求解根据连续性方程导出的压力修正方程，得到压力场的修正值；根据压力修正值改进速度场和压力场；根据改进后的速度场求解其它变量的方程；依次迭代直到收敛。N-S方程中速度和压力的耦合算法(4) 采用同位网格的SIMPLE算法什么是同位网格(collocated grid)？-将每个速度分量(u, v, w

24、)及压力p都存储于同一套网格系统中。以二维为例：WEu, v和压力p均在P点离散NvnuwPvsueSN-S方程中速度和压力的耦合算法(5)为什么要采用同位网格？采用同位网格计算时在主节点上求得的流速和压力修正值还要插值到界面上去（这一点下面将述及），因此它的计算步骤要比交错网格复杂一些，单纯从这方面来讲， 或许对于二维的直角坐标来说，采用同位网格并不会比交错网格有什么优势。但是考虑到同位网格的生成要比交错网格简单，而且由于速度在主节点上离散，u, v动量离散方程的系数除了源项之外都是相同的，可一次计算得到。相比之下，还是同位网格可取。而当计算由二维扩展到三维，由规则区域发展到不规则区域，由正

25、交坐标发展到非正交坐标，由结构化网格发展到非结构化网格的时候，同位网格在编制程序方面的优点便非常突出了。N-S方程中速度和压力的耦合算法(6)在同位网格上实现SIMPLE算法的步骤：根据上一层次计算得到的界面流速和压力求解动量方程， 得到新的速度场；根据新的速度场按动量插值的方式求得界面流速，进而求解压力修正方程；根据压力修正值来计算界面上的速度修正值，进而得到节点上的速度修正值(其中要用线性插值的方式计算界面压力)；根据速度和压力的修正值更新流场和压力场最关键的问题是，如何将界面上的压力引入动量方程的求解当 中，从而保证速度和压力的耦合！N-S方程中速度和压力的耦合算法(7) 最初SIMPL

26、E算法的提出是为了解决不可压流体的速度和压力耦合问题，现在的SIMPLE算法也可以用来求解可压缩流体流动； 有许多针对SIMPLE算法的改进，以提高其收敛性，如SIMPLER、SIMPLREC、SIMPLEX 等算法；N-S方程中速度和压力的耦合算法(8) PISO(Pressure Implicit with Splitting of Operators，压力隐式分裂算子)算法属于SIMPLE算法族的一部分其主要思想是使用一个或更多的PISO循环来代替SIMPLE算法中压力校正阶段所需的重复计算，校正后的速度更接近满足连续性和动量方程。这一迭代过程被称为动量校正；另外，PISO还提供了一个偏

27、斜矫正的过程-极大的减少了计算高扭曲度网格所带来的收敛性问题， 可保证在基本相同的迭代步中，得到和正交性较好的网格同样好的解；对于所有的过渡流动问题，推荐使用PISO算法代数方程组的求解方法(1) TDMA算法是求解三对角阵型方程组的直接解法；一维流动问题的离散方程，其系数构成典型的三对角阵；对于一维问题，使用TDMA算法，其求解方法直接， 编程简单，计算量也不大；对于二维和三维的情况，可使用TDMA算法来求解迭代求解过程中的直接求解部分。代数方程组的求解方法(2) 点迭代法每一步只能改进求解区域中的一个节点之值Jacobi迭代：任一点上未知值的更新采用上一轮迭代中所获得的邻点之值来计算。其收

28、敛速度很慢， 一般很少使用。但可用来解决强烈非线性问题的 收敛问题；Gauss-Seidel迭代：每一步计算采用邻点的最新值进行；SOR(逐次超松弛)/SUR(逐次亚松弛)迭代法：对于非线性的求解一般采用亚松弛(松弛因子小于1)， 在一定范围内，收敛因子越大，收敛越快。但当收敛因子超过一定值时，解会出现振荡，甚至发散；代数方程组的求解方法(3) 块迭代法 把求解区域分成若干块，同一块中各节点之值采用代数方程的直接解法求得，从一块到另一块的推进采用迭代的方式进行； 获得收敛解的外迭代次数减小，但每一轮迭代中的代数运算次数增加。 交替方向隐式迭代法(ADI) 采用先逐行扫描，再逐列扫描的方法进行块

29、迭代，两次扫描组成一轮迭代； 可加快收敛速度，是采用较多的一种方法。 强隐迭代法(SIP) 比上面的ADI方法采用更高的直接解法的比例，收敛速度更快； 适合于结构化网格所形成的代数方程组的求解。 多重网格法 这是促进代数方程迭代求解收敛速度的有效方法，在最近二十余年得到广泛应用。具体请参考相关书籍。七、CFD的初始和边界条件 为什么要使用初始和边界条件？ 常用边界条件的简单介绍为什么要使用初始和边界条件 初始条件 对于瞬态过程的模拟，需要给定初始条件。这决定了物理过程基于何种状态进行发展，因此也就决定了数值求解结果的准确性。 边界条件上一章节中所讲的代数方程的求解，必须给定边界条件，这是为了保证方程组的封闭，得到唯一的解；边界条件的给定，需要根据具体的物理过程来确定，一定要保证边界条件的合理性，否则所得到的解就是不合理的，甚至代数方程的求解根本没法进行。常用边界条件的简单介绍(1) 给定流速或流量 一般用在进口处，但在某些情况下也可用于出口； 一般还要同时给定其它参数的分布情况(压力除外)； 对称边界 一般用于几何的对称面(但不表示所有的几何对称面都可应用此类边界条件)； 给定垂直于对称面的速度分量为零，对称面上其它物理