计算流体力学基础简介

计算流体力学基础简介本硕111班王鹏主要内容：•流体力学的控制方程组•适合CFD使用的控制方程组•偏微分方程的离散化——有限差分•计算流体力学的基本方法流体力学的控制方程组一.粘性流动的N-S方程1.连续性方程非守恒形式守恒形式2.动量方程非守恒形式流体力学的控制方程组守恒形式3.能量方程非守恒形式流体力学的控制方程组守恒形式流体力学的控制方程组一.无粘流Euler方程1.连续性方程非守恒形式守恒形式2.动量方程非守恒形式流体力学的控制方程组守恒形式3.能量方程非守恒形式守恒形式流体力学的控制方程组对于牛顿流体，斯托克斯在1845年得到，在上面的方程组中：其中，𝝁是分子粘性系数，𝝀是第二粘性系数。斯托克斯提出假设，认为：根据傅里叶热传导定律，热传导产生的热流与当地的温度梯度成正比:其中k为热导率。适合CFD使用的控制方程组守恒形式的连续性方程、动量方程和能量方程可以使用同一个通用方程来表达，这有助于计算程序的简化和结构组织，为算法设计和编程提供了方便。守恒形式的通用方程（1）如果将上式中的U、F、G、H和J看成列向量，则上式就可代表整个守恒形式的控制方程。这些列向量为：适合CFD使用的控制方程组（1）式中的F、G和H称为通量项(通量向量），J代表源项，U称为解向量。对于一般的时间推进方法，我们可将（1）式变换一下以方便求解利用上一个时间步的结果计算出方程右边的这些项，再解出列向量U，这样就得到了𝝆、𝒖𝝆、v𝝆、w𝝆和𝝆(e+𝒗𝟐𝟐)，这些量被称为守恒变量。适合CFD使用的控制方程组根据下式，即可得到𝝆、𝒖、v、w和e，这些量被称为原始变量。针对Euler方程也可进行类似分析。偏微分方程的离散化——有限差分偏微分方程的解析解是封闭形式的表达式，它给出了未知函数在区域内的连续变化。偏微分方程的离散化是用另外一个类似的表达式来近似这些偏微分方程，但是这个近似表达式只在去域内有限多个离散点或控制体上规定了取值。利用泰勒级数展开可以推导出各阶导数的有限差分的一般形式。（2）（3）利用（2）、（3）式可以得到以下一些差分公式：偏微分方程的离散化——有限差分一阶导数的一阶向前差分一阶导数的二阶中心差分一阶导数的一阶向后差分二阶导数的二阶中心差分一阶混合导数的二阶中心差分……..偏微分方程的离散化——有限差分利用多项式也可以实现有限差分。例如，如右图所示的边界上，u可以表示成多项式：在网格点1，y=0,由上式给出，在网格点2，y=∆𝑦,在网格点3，y=2∆𝑦（4）（5）（6）（7）偏微分方程的离散化——有限差分解出b,再对（4）式求导，边界上y=0,代入上式有，最终得到上式就是单侧差分表达式，之所以被称为单侧，是因为它只用到了边界一侧的网格点上的信息。计算流体力学的基本方法计算流体力学包括很多方法，其中较基础的有拉克斯-温德罗夫方法、麦考马克方法、空间推进方法、松弛法、交替方向隐式方法、压力修正法等。1.拉克斯-温德罗夫方法以非定常二维无粘流为例，控制方程如下式中，P=𝜌𝑅𝑇,e=𝐶vT（8）（9）（10）（11）计算流体力学的基本方法如右图所示，在二维矩形网格中，用𝜌𝑖,𝑗𝑡表示在t时刻网格点(i,j)的密度，用𝜌𝑖,𝑗𝑡+∆t表示在t+∆𝑡时刻网格点(i,j)的密度。其他流场变量亦是同理。（12）（13）（14）（15）已知t时刻的𝜌𝑖,𝑗𝑡，u𝑖,𝑗𝑡，𝑣𝑖,𝑗𝑡，e𝑖,𝑗𝑡，只要能求出它们在t时刻的时间导数，根据上面四式即可将每个网格点处的流场变量推进到下一时刻。计算流体力学的基本方法比如，用（12）式计算t+∆𝑡时刻的密度。将方程（8）中的空间导数由二阶中心差分给出：（17）再将方程（8）对时间求导：上式中出现二阶混合导数，可由方程式（8）到（11）对相应的空间自变量求导得到。例如（16）将右边所有的项还是由t时刻的二阶中心差分表示，即计算流体力学的基本方法可用同样的方法处理（17）式中其他的二阶导数项，然后把得到的所有方程式都代入式（12）中，即可得到t+∆𝑡时刻的密度。另外，对于式（16）中出现的一阶空间导数仍由二阶中心差分代替，例如按照同样的方法可得到在t+∆𝑡时刻网格点(i,j)的其他流场变量。右图给出了计算过程中用到的网格，图中显示了t时刻和t+∆𝑡时刻两个时间平面上的空间网格。从右图中可以看到，利用t时刻流场变量在(i,j)(i+1,j)(i-1,j)(i,j+1)(i,j-1)处的已知值，拉克斯-温德罗夫方法可以显示地求出(i,j)处流场变量在t+∆𝑡时刻的值。流场中所有网格点处的流场变量在在t+∆𝑡时刻的值都可以通过这样的方法得到。拉克斯-温德罗夫方法在时间和空间上都有二阶精度，其优点是思路直观简单，缺点是代数表达式冗长。计算流体力学的基本方法二.麦考马克（Macormack）方法麦考马克方法是从拉克斯-温德罗夫方法中变化过来的一种方法，但更简单，在1969年一经提出，就成为解决流动问题最流行的显示有限差分方法，一直流行了15年，现在它已被更先进的方法取代了。和拉克斯-温德罗夫方法一样，我们的目的还是由t时刻网格点处的所有流场变量求出t+∆𝑡时刻的所有网格点变量。以密度为例，在麦考马克方法中，式中，(𝜕𝜌𝜕𝜌)𝑎𝑣表示𝜕𝜌𝜕𝑡在t时刻和t+∆𝑡时刻之间的平均值。其他流场变量也有类似关系（18）计算流体力学的基本方法预估步：在连续方程（8）中，用向前差分代替方程右边的空间导数：取泰勒级数前两项，即可求密度的估计值（𝜌）𝑡+∆t（19）同理可求出u、v和e的预估值，即计算流体力学的基本方法（20）校正步：先用向后差分替代连续方程右边的空间导数，然后用𝜌、u、v和e的预估值进行计算，就得到了时间导数在t+∆𝑡时刻的预估值式（18）中密度的时间导数的平均值把式（19）（20）的值代入即可。得到了这个平均值之后，就可以从式（18）中得到t+∆𝑡时刻密度的最终校正值，即上面描述的预估-校正方法，给出了t+∆𝑡时刻网格点（i,j）的密度值。在所有的网格点中重复这一过程，可得到t+∆𝑡时刻整个流场的密度值。同理可得到t+∆𝑡时刻整个流场的u、v和e的值。在使用预估-校正方法时，预估步用向前差分，校正步用向后差分，使得麦考马克方法具有二阶精度（？）。该方法避免了二阶导数的计算，从而简化了计算。计算流体力学的基本方法（21）松弛法是一种适合于求解椭圆型偏微分方程的迭代法。例如，亚声速无粘流动由椭圆型偏微分方程控制松弛法经常被用来求解亚声速的低速流动。松弛法可以是显式的也可以是隐式的，这里介绍显式松弛法。考虑无粘不可压流体的二维无旋流动，用∅表示速度势，则控制方程为我们将在右图所示的网格上数值求解方程（21），利用二阶导数的二阶中心差分式代替上式的偏导数，得仔细考察右图可知，网格点1~20构成了流场的边界条件。为了使椭圆型方程成为适定的，必须在包围区域的整个边界上给定边界条件。这意味着∅1~∅20都是已知的。从式（22）中解出∅𝑖,𝑗，有上式中，上标n和n+1表示迭代次数。（22）（23）三.松弛法计算流体力学的基本方法（24）在对所有内部网格点都使用了方程（23）后，就完成了第一次迭代，n=1。接着进行第二次迭代，n=2。这个迭代过程反复不断，直至收敛到解。（为什么一定是朝解的方向收敛？）再具体些，将方程（23）应用到网格点21，进行第n+1次迭代，有式中，∅21𝑛+1是未知量，∅22𝑛和∅24𝑛是上一次迭代中已经求出来的，∅20和∅2是由边界条件给定的。也可以将更新后的∅值更快的用到方程（23）中，例如，从方程（24）中解出∅21𝑛+1后，转移到网格点22中，从方程（23）得到在这种方式中，可沿一条给定的水平线或竖直线逐点第n+1次迭代时的未知量∅。（这种迭代法叫高斯-赛德尔迭代法）。进行多次迭代后，当所有的网格点处的变得小于预定值时，就认为迭代收敛了。通常，运用逐次超松弛法可加快收敛过程，思路如下：把方程（23）的结果看成是∅𝑖,𝑗的中间值，记作∅𝑖,𝑗𝑛+1，即（25）计算流体力学的基本方法最后，用上一次迭代得到的∅𝑖,𝑗𝑛和由式（25）得到的∅𝑖,𝑗𝑛+1，外推出∅𝑖,𝑗𝑛+1的值：式中，𝜔是松弛因子。若𝜔1，上面的过程就是逐次超松弛法；若𝜔1，就叫做逐次低松弛法。当收敛过程表现为在某个值附近来回摆动时，通常会用到地松弛法。对超松弛法，通常1𝜔2。选取适当的值，可以减少迭代次数。计算流体力学的基本方法（26）四.交替方向隐式（ADI）方法讲这个方法之前，有必要先介绍一下显示方法和隐式方法。所有的有限差分方法都可以归结为这两种不同的通用方法——显示方法或隐式方法。显示方法。以一位热传导方程为例等式左边用向前差分，等是右边用二阶中心差分，即整理此式后可写成方程（26）是抛物线型偏微分方程，可推进求解。此时，推进变量是时间t，为更简明起见，考虑右图所示的有限差分网格。假设在第n个时间层上的每个网格点上温度T是已知的，时间推进就意味着第n+1个时间层每个网格点上的T值都要用第n个时间层上的已知量计算出来。（28）（27）计算流体力学的基本方法（29）例如，在网格点2处，式（28）将写成按照相同方法，可相继求出𝑇3𝑛+1、𝑇4𝑛+1等。于是，我们可以给出定义：显示方法中每个差分方程只包含一个未知数，从而这个未知数可以用直接计算的方式显示地求解。隐示方法我们可以再大胆些，把式（27）右边的空间差分写成第n个时间层的量与第n+1个时间层上的量的平均值，即式（29）中所用的这种特定形式的差分称为克兰克-尼克尔森格式（克兰克-尼克尔森格式广泛被用于求解抛物型方程问题）。式（29）中包含三个未知量𝑇𝑖−1𝑛+1、𝑇𝑖𝑛+1、𝑇𝑖+1𝑛+1，暂时求解不了。需要把排列在同一时间层上的同一条直线上的未知量联立起来形成代数方程组才能解出这些未知量。计算流体力学的基本方法（30）考虑右图所示的有限差分网格，以七点空间网格点作为具体例子。整理式（29），将未知量放在等式左边，已知量放到等式右边，得为简化书写，用A,B和𝐾𝑖（均是已知量）表示上式的系数和右端项于是，式（30）可写成在网格点1和7上给定了边界条件，所以𝑇1也是已知的，关于𝑇1的项也应放到等式右边，即用𝐾2‘表示𝐾2−𝐴𝑇1，则𝐾2‘也是已知的，即（31）计算流体力学的基本方法把方程（31）~（35）写成矩阵形式：（32）（33）（34）（35）这是一个三对角矩阵，可用托马斯算法（国内称为追赶法）求解。解出之后即得到第n+1个时间层上的温度T值。重复以上过程，即可得到最终的数值解。计算流体力学的基本方法（37）现在，介绍交替方向隐式（AlternatingDirectionImplicit,ADI）方法。以二维热传导方程为例，对上式运用克兰克-尼克尔森方法，有上式有五个未知量，即𝑇𝑖−1,𝑗𝑛+1、𝑇𝑖,𝑗𝑛+1、𝑇𝑖+1,𝑗𝑛+1、、𝑇𝑖,𝑗−1𝑛+1和𝑇𝑖,𝑗+1𝑛+1，这让我们得不到三对角方程组。但运用交替方向隐式方法，可使得用三对角形式就可以求解方程（37）。方程（36）是通过推进方法求解的，即T(t+∆𝑡)是根据已知量T(t)用某种方式得到的。现在让我们分两步得出T(t+∆𝑡)。第一步的时间步长为∆𝒕/𝟐，用中心差分代替方程（36）中的空间导数，但只对x的导数采用隐式处理，即（36）（38）计算流体力学的基本方法（39）方程（38）可化为三对角形式式中如右图所示，对每一个固定的j，我们在x方向扫描，对该排所有的i写出方程，联立成方程组，用托马斯算法求解方程组，得到𝑇𝑖,𝑗𝑛+1/2。扫描一次就用一次托马斯算法。然后在标号j+1的下一排网格点上扫描……当这一步结束时，中间时刻t+∆𝑡/2所有网格点上的𝑇𝑛+1/2就都求出来了。ADI方法的第二步是利用t+∆𝑡/2时刻的已知值求