第13章_计算流体力学CFD(5)

1、 Lax-Wendroff方法方法 是一种显式有限差是一种显式有限差 分方法，适合于推分方法，适合于推 进求解。进求解。 二维时间推进网格二维时间推进网格 Lax-Wendroff方法方法 在时间和空间上都在时间和空间上都 具有二阶精度。具有二阶精度。 二维时间推进网格二维时间推进网格 非定常二维无粘流（欧拉方程非守恒形式）：非定常二维无粘流（欧拉方程非守恒形式）： Lax-Wendroff显式推进求解显式推进求解 (沿时间方向进行泰勒级沿时间方向进行泰勒级 数展开数展开)： 空间导数采用中心差分：空间导数采用中心差分： ( ) 求对时间求对时间t的二阶导数：的二阶导数： Lax-Wendro

2、ff显式推进求解显式推进求解 ： MacCormack方法在时间和空间上都方法在时间和空间上都 具有二阶精度。具有二阶精度。 MacCormack方法是一种显式有限差分方法是一种显式有限差分 方法，适合于推进求解。方法，适合于推进求解。 MacCormack方法比方法比Lax-Wendroff方方 法应用起来更简单。法应用起来更简单。 校正步校正步 预估步预估步 预估步：空间导数用向前差分计算。预估步：空间导数用向前差分计算。 预估步：空间导数用向前差分计算。预估步：空间导数用向前差分计算。 预估值：预估值： 校正步：空间导数用向后差分计算。校正步：空间导数用向后差分计算。 在在MacCorm

3、ack方法中，预估步用向前差分，方法中，预估步用向前差分， 校正步用向后差分；也可以预估步用向后差分，校正步用向后差分；也可以预估步用向后差分， 校正步用向前差分。或者在时间推进解法的相校正步用向前差分。或者在时间推进解法的相 继两个时间步中轮流使用这两种办法。继两个时间步中轮流使用这两种办法。 粘性流动的控制方程是粘性流动的控制方程是N-S方程。方程。 对定常流动，对定常流动，N-S方程的数学性质更多地表现为方程的数学性质更多地表现为 椭圆型的，不能采用椭圆型的，不能采用Lax-Wendroff方法和方法和 MacCormack方法求解。方法求解。 对非定常流动，可以采用对非定常流动，可以采

4、用Lax-Wendroff方法或方法或 MacCormack方法求解方法求解N-S方程。方程。 2 2 v e v u U 可以采用可以采用Lax-Wendroff方法或方法或 MacCormack方法求解方法求解U的分的分 量在各时间步的值。量在各时间步的值。 非定常守恒形式欧拉方程（二维）：非定常守恒形式欧拉方程（二维）： 定常守恒型二维欧拉方程：定常守恒型二维欧拉方程： 对于亚声速流动，上述对于亚声速流动，上述 方程是椭圆型的，所有方程是椭圆型的，所有 空间推进方法都不适用，空间推进方法都不适用， MacCormack方法也不方法也不 适用。适用。 对于超声速流动，上述方对于超声速流动，

5、上述方 程是双曲型的，空间推进程是双曲型的，空间推进 方法适用，方法适用，MacCormack 方法也适用。方法也适用。 定常守恒型二维欧拉方程：定常守恒型二维欧拉方程： MacCormack方法：方法： 定常守恒型二维欧拉方程：定常守恒型二维欧拉方程： 预测步预测步：（向前差分）：（向前差分） 预估值：预估值： 预估值：预估值： 校正步校正步：（向后差分）：（向后差分） 松弛法特别适合于求解椭圆型偏微分方程，松弛法特别适合于求解椭圆型偏微分方程， 常被用来求解无粘亚声速的低速流动。常被用来求解无粘亚声速的低速流动。 考虑无粘不可压流体的二维无旋流动，控考虑无粘不可压流体的二维无旋流动，控 制

6、方程为制方程为Laplace方程：方程： 松弛法是一种迭代法松弛法是一种迭代法 上标上标n和和n+1表示迭代次数表示迭代次数 松弛法是一种迭代法松弛法是一种迭代法 松弛法是一种迭代法松弛法是一种迭代法 松弛法是一种迭代法松弛法是一种迭代法 从左至右扫描从左至右扫描 松弛法是一种迭代法松弛法是一种迭代法 当所有网格点处的当所有网格点处的 都小于一个预定的值时，迭代都小于一个预定的值时，迭代 收敛。收敛。 n ji n ji, 1 , 运用逐次松弛法可加运用逐次松弛法可加 快收敛的过程。快收敛的过程。 从左至右扫描从左至右扫描从下至上扫描从下至上扫描 运用逐次松弛法可加运用逐次松弛法可加 快收敛的

7、过程。快收敛的过程。 是松弛因子，如果是松弛因子，如果 1，叫做逐次超松弛法；，叫做逐次超松弛法； 如果如果 1，叫做逐次低松弛法。，叫做逐次低松弛法。 运用逐次松弛法可加运用逐次松弛法可加 快收敛的过程。快收敛的过程。 选取合适的选取合适的 值，可以减少迭代次数，从而减少计算值，可以减少迭代次数，从而减少计算 时间。在某些问题中，迭代次数可减少到原来的时间。在某些问题中，迭代次数可减少到原来的1/30 一维波动方程：一维波动方程： 差分方程：差分方程： 截断误差：截断误差： 差分方程：差分方程： 泰勒级数展开：泰勒级数展开： 差分方程：差分方程： 将泰勒级数展开代入差分方程得：将泰勒级数展开

8、代入差分方程得： 差分方程：差分方程： 将泰勒级数展开代入差分方程得：将泰勒级数展开代入差分方程得： 差分方程：差分方程： 等号右边将对等号右边将对t的偏导数转化为对的偏导数转化为对x的偏导数得：的偏导数得： 差分方程：差分方程： 偏微分方程（修正方程）：偏微分方程（修正方程）： 一维波动方程（偏微分方程）：一维波动方程（偏微分方程）： 差分方程：差分方程： 一维波动方程（偏微分方程）：一维波动方程（偏微分方程）： 差分方程的精确解是上述一维波动方程的数值解（含误差）差分方程的精确解是上述一维波动方程的数值解（含误差） 差分方程：差分方程： 差分方程的精确解是上述修正方程的精确解（不含误差）差

9、分方程的精确解是上述修正方程的精确解（不含误差） 偏微分方程（修正方程）：偏微分方程（修正方程）： 差分方程：差分方程： 偏微分方程（修正方程）：偏微分方程（修正方程）： 一维波动方程（偏微分方程）：一维波动方程（偏微分方程）： 修正方程等号右端的项是截断误差，如果截断误差的主项修正方程等号右端的项是截断误差，如果截断误差的主项 是偶数阶导数，数值解将主要表现出耗散行为；如果主项是偶数阶导数，数值解将主要表现出耗散行为；如果主项 是奇数阶导数，数值解将主要表现出色散行为。是奇数阶导数，数值解将主要表现出色散行为。 偏微分方程（修正方程）：偏微分方程（修正方程）： 等号右端的偶数阶导数项起数值耗

10、散的作用，奇数阶导数等号右端的偶数阶导数项起数值耗散的作用，奇数阶导数 项起数值色散的作用。项起数值色散的作用。 偏微分方程（修正方程）：偏微分方程（修正方程）： 数值耗散的作用很象物理粘性，二阶导数项前的系数被称数值耗散的作用很象物理粘性，二阶导数项前的系数被称 为人工粘性。为人工粘性。 偏微分方程（修正方程）：偏微分方程（修正方程）： 数值耗散的影响会将波抹平数值耗散的影响会将波抹平 色散导致波的不同相位在传播中产生畸变，色散导致波的不同相位在传播中产生畸变， 表现为波前和波后出现振荡。表现为波前和波后出现振荡。 尽管人工粘性降低了解的精度，但通常有助于提高解的稳尽管人工粘性降低了解的精度

11、，但通常有助于提高解的稳 定性。定性。 偏微分方程（修正方程）：偏微分方程（修正方程）： 考虑二维热传导方程：考虑二维热传导方程： 等号右端有五个未知量，不能得到三对角方程组，不能采等号右端有五个未知量，不能得到三对角方程组，不能采 用托马斯算法（追赶法）求解。用托马斯算法（追赶法）求解。 采用采用Crank-Nicolson方法（隐式）：方法（隐式）： 考虑二维热传导方程：考虑二维热传导方程： 第一步：第一步： 时间步长为时间步长为 ，空间导数采用中心差分，只对，空间导数采用中心差分，只对x的的 导数采用隐式处理。导数采用隐式处理。 第一步：第一步： 简化为三对角形式简化为三对角形式 第一步

12、：第一步： 对每一个固定的对每一个固定的j，对所有，对所有 的的i联立形成方程组。联立形成方程组。 对不同的对不同的j，重复上述过程。，重复上述过程。 考虑二维热传导方程：考虑二维热传导方程： 第二步：第二步： 时间步长为时间步长为 ，空间导数采用中心差分，只对，空间导数采用中心差分，只对y的的 导数采用隐式处理。导数采用隐式处理。 第二步：第二步： 简化为三对角形式简化为三对角形式 第二步：第二步： 对每一个固定的对每一个固定的i，对所有，对所有 的的j联立形成方程组。联立形成方程组。 对不同的对不同的i，重复上述过程。，重复上述过程。 两步结束之后，两步结束之后，T在时间方向上推进了一个时

13、间步长在时间方向上推进了一个时间步长 t. 考虑二维热传导方程：考虑二维热传导方程： 推进过程只涉及三对角方程组。推进过程只涉及三对角方程组。 第一步，差分方程的第一步，差分方程的x方向是隐式的。方向是隐式的。 考虑二维热传导方程：考虑二维热传导方程： 所以这种方法叫交替方向隐式方法所以这种方法叫交替方向隐式方法(Alternating Direction Implicit, ADI) 第二步，差分方程的第二步，差分方程的y方向是隐式的。方向是隐式的。 考虑二维热传导方程：考虑二维热传导方程： ADI格式对格式对t,x,y都是二阶精度的都是二阶精度的 截断误差为：截断误差为： 不可压无粘流动受

14、椭圆型偏微分方程控制（不可压欧不可压无粘流动受椭圆型偏微分方程控制（不可压欧 拉方程），松弛法是求解椭圆型问题经典的数值方法，拉方程），松弛法是求解椭圆型问题经典的数值方法， 本质上是一个迭代过程。本质上是一个迭代过程。 不可压粘性流动的控制方程是不可压的不可压粘性流动的控制方程是不可压的NS方程，这方程，这 个方程具有椭圆型和抛物型的混合特性，松弛法不是个方程具有椭圆型和抛物型的混合特性，松弛法不是 特别适用。特别适用。 压力修正法也是一种迭代过程，在不可压压力修正法也是一种迭代过程，在不可压NS方程的方程的 数值求解中得到了广泛的应用。数值求解中得到了广泛的应用。 假设假设 =常数，常数，

15、 =常数，可压缩常数，可压缩NS方程转化为不可方程转化为不可 压压NS方程：方程： 上述四个方程封闭，含上述四个方程封闭，含 四个未知数。四个未知数。 二维不可压流体的连续二维不可压流体的连续 性方程为：性方程为： 中心差分格式为：中心差分格式为： 右上角是右上角是u的值，的值， 左下角是左下角是v的值的值 速度会出现右图的速度会出现右图的 棋盘式分布棋盘式分布 右上角是右上角是u的值，的值， 左下角是左下角是v的值的值 可压流动中不会发生右可压流动中不会发生右 图的问题，因为连续性图的问题，因为连续性 方程中包含了密度对时方程中包含了密度对时 间和空间的变化。间和空间的变化。 0V t 在可

16、压缩流动中，右图在可压缩流动中，右图 速度的棋盘分布经过一速度的棋盘分布经过一 个时间步就会被抹平。个时间步就会被抹平。 二维不可压流体压力梯二维不可压流体压力梯 度采用中心差分：度采用中心差分： 棋盘式的离散压力分棋盘式的离散压力分 布布 压力会出现右图的压力会出现右图的 棋盘式分布棋盘式分布 在交错网格上使用中心在交错网格上使用中心 差分就不会出现速度和差分就不会出现速度和 压力的棋盘式分布问题。压力的棋盘式分布问题。 交错网格交错网格 在在(i-1,j), (i,j), (i+1,j), (i,j+1),(i,j-1)等图中的实等图中的实 心原点上计算压力心原点上计算压力 交错网格交错网

17、格 在在(i-1/2,j), (i+1/2,j)等图等图 中的空心原点上计算中的空心原点上计算u 交错网格交错网格 在在(i,j-1/2), (i,j+1/2)等图等图 中的空心原点上计算中的空心原点上计算v 连续性方程在网格点连续性方程在网格点(i,j) 的中心差分表达式为：的中心差分表达式为： 交错网格交错网格 压力修正法本质上是一种迭代法，思路如下：压力修正法本质上是一种迭代法，思路如下： 1）迭代开始时，先给定压力的初始近似）迭代开始时，先给定压力的初始近似p* 2）用）用p*的值从动量方程中求解的值从动量方程中求解u,v,w,得到与得到与p*有关有关 的的u*,v*,w* 压力修正法

18、本质上是一种迭代法，思路如下：压力修正法本质上是一种迭代法，思路如下： 修正后的压力为修正后的压力为 3）将）将u*,v*,w*代入连续性方程，它们不一定满足连代入连续性方程，它们不一定满足连 续性方程。用连续性方程构造压力的修正量续性方程。用连续性方程构造压力的修正量 ,加加 到到p*上，使速度场满足连续性方程。上，使速度场满足连续性方程。 修正后的速度为修正后的速度为 速度修正量速度修正量 可以从可以从 得到。得到。 压力修正法本质上是一种迭代法，思路如下：压力修正法本质上是一种迭代法，思路如下： 4） 用步骤用步骤3)中修正后的压力做为新的中修正后的压力做为新的p*，回到步，回到步 骤骤

19、2)。重复这个过程，直到速度场满足连续性方程。重复这个过程，直到速度场满足连续性方程 为止。为止。 这样就得到修正好了的流场。这样就得到修正好了的流场。 压力修正公式为：压力修正公式为： 压力修正公式为：压力修正公式为： 上述压力修正公式具有椭圆型的性质，可以用松弛上述压力修正公式具有椭圆型的性质，可以用松弛 法数值求解。法数值求解。 在不可压流场中，压力的扰动将会传遍整个流场，在不可压流场中，压力的扰动将会传遍整个流场， 这与上述方程的椭圆型性质相吻合。这与上述方程的椭圆型性质相吻合。 压力修正公式为：压力修正公式为： 压力修正公式是压力修正压力修正公式是压力修正 的泊松方程的中心差的泊松方

20、程的中心差 分表达式。分表达式。 上述泊松方程中的二阶偏导数用中心差分替代。上述泊松方程中的二阶偏导数用中心差分替代。 tdQ/式中：式中： 压力修正压力修正 的泊松方程的泊松方程(为椭圆型为椭圆型)： tdQ/ d相当于一个质量源项。相当于一个质量源项。 SIMPLE是是Semi-implicit method for pressure-linked equation (压力耦合方程的半隐式算法压力耦合方程的半隐式算法)的缩写。的缩写。 SIMPLE算法的步算法的步 骤如下：骤如下： 1）在右图所示的交）在右图所示的交 错网格上分别给出错网格上分别给出 * , n p * , n u * n

21、 v SIMPLE算法的步算法的步 骤如下：骤如下： 2）求出）求出 1 * , n u 1 * n v 采用动量方程求解。采用动量方程求解。 2） 的求法：的求法： 1 * n u X方向的动量方程：方向的动量方程： 2） 的求法：的求法： 1 * n u 在在a点：点： 在在b点：点： X方向的动量方程：方向的动量方程： 差分方程：差分方程： 差分方程：差分方程： X方向的动量方程：方向的动量方程： 2） 的求法：的求法： 1 * n u 2） 的求法：的求法： 1 * n u 2） 的求法：的求法： 1 * n v Y方向的动量方程：方向的动量方程： 2） 的求法：的求法： 1 * n v 在在c点：点： 在在d点：点： Y方向的动量方程：方向的动量方程： 差分方程：差分方程： 2） 的求法：的求法： 1 * n v 3）将）将 和和 代入压力修正公式，在所有内代入压力修正公式，在所有内 部网格点上求解部网格点上求解 1 * n u 1 * n v SIMPLE算法