第二章控制方程离散化之低阶格式（一）计算流体力学

好久没有更新了，每次敲完公式，都要自己添加\nonumber才不会显示标号，对CSDn也是服了。
最近也不知道忙什么，课也不多，琐事倒是不少，正所谓春水碧于天，划船听雨眠，好了，进入正题
目标

通过本节，了解有限体积法基本原理，通过中心差分、一阶迎风、混合迎风了解这些格式的优缺点

2.1离散化概述在对CFD计算之前，需要对计算区域离散化，也就是说要对空间连续的计算区域进行划分，把它划分为许多离散子区域，通过确定每个区域中的节点，从而生成网格，说人话就是把偏微分方程转化为代数方程组进行求解，对于瞬态，还要对时间进行离散化，那我们先从空间离散进行讨论吧
2.1.1离散化目的目的我就不说了，还不是因为偏微分方程那魔性，至今还没有办法求出真解，但是，还好，计算机最不怕的就是迭代，我们可以通过逼近，得到最接近真解的数值解。
2.1.2离散化使用的网格网格的质量对计算结果影响非常大，至于原因，参考积分怎么求的吧，在二维时，常见的网格有三角形，四边形单元，在三维问题，有四面体、六面体，棱锥， 混合网格等单元
2.1.3 常见离散化方法 1.有限差分法
这个方法在教程中经常提到，通过将求解区域划分为差分网格，用有限个网格节点代替连续的求解域，用差商代替偏导数，这个方法可以用于双曲型和抛物型问题，对于复杂的问题就呵呵了
2.有限元法
有限元法在固体力学占据着主导地位，原因呢，知乎找找就有了，由于其在求解CFD时较慢，在CFD领域用的并不广泛，有限元就是将一个连续的求解区域任意划分为任意形状的许多微小单元，对微小单元构成的方程组进行求解。
3.有限体积法
有限体积法（Finite Volume Method）在流体力学领域用的可就非常广泛了，下面我们开始介绍它
2.2有限体积法及其网格简介 2.2.1 基本思想
有限体积法又叫控制体积法，基础思路为：把计算区域划分为网格，并使每个网格点周围有一个互不重复的控制体积；将待求解微分方程对每一个控制体积积分，从而得到一组离散的方程。其中的未知数是网格点上的因变量 ?? 。为了求出控制体积的积分，必须假设 ?? 值在网格点之间的变化规律。从积分区域的选取方法看来，有限体积属于加权余量法中的子域法，从未知解的近似方法看来，有限体积法属于局部近似的离散方法。简言之，子域法加离散，就是有限体积的基本方法。
有限体积法得出的离散方程，要求因变量的积分守恒对任意一组控制体积都得到满足，对整个计算区域，自然也得到满足。这就是其优点，对于其他的离散方法，如差分，必须网格及其密时，离散方程才满足积分守恒，而有限体积法即使在粗网格情况下，也显示出准确的积分守恒
2.3一维稳态问题的有限体积法
控制方程

（dρu?）dx=ddx(Γd?dx)+S（ d ρ u ? ） d x = d d x ( Γ d ? d x ) + S
该方程又称为一维模型方程，方程中包含对流项、扩散项及源项，方程中的 ??? ?是广义变量，可以为速度。温度或浓度等一些待求的物理量， Γ?Γ ?是相应于 ??? ?的广义扩散系数， S?S ?是广义源项，变量 ??? ?在端点A和B的边界值为已知。
【第二章控制方程离散化之低阶格式（一）】第一步：生成计算网格 δδ

建立离散方程，对控制体积P进行积分，

∫ΔVdρu?dxdV=∫ΔVddx(Γd?dx)dV+∫ΔVSdV∫ Δ V d ρ u ? d x d V = ∫ Δ V d d x ( Γ d ? d x ) d V + ∫ Δ V S d V
式中， ΔVΔ V是控制体积的体积值，当控制体积很微小时， ΔVΔ V可以表示为 ΔV+AΔ V + A，这里 AA是控制体积界面的体积
这样积分就可以进行计算了，这个简化过程很重要

(ρu?A)e?(ρu?A)w=(ΓAd?dx)e?(ΓAd?dx)w+SΔV(1)(1) ( ρ u ? A ) e ? ( ρ u ? A ) w = ( Γ A d ? d x ) e ? ( Γ A d ? d x ) w + S Δ V
目前我们的位置数为界面 ww和 ee处的 ρρ、 uu、 ΓΓ、 ??d?dxd ? d x。为了计算界面出的未知物理量，需要知道物理参数在节点间的近似分布。简单的说，如果是线性分布的话，我们采用的方法就叫做 中心差分，对于均匀网格，各个物理量可以用下面的公式进行求解
ΓeΓw?e?w(d?dx)e(d?dx)w=ΓP+ΓE2=ΓW+ΓP2=?P+?E2=?W+?P2=?E??P(δx)e=?P??W(δx)eΓ e = Γ P + Γ E 2 Γ w = Γ W + Γ P 2 ? e = ? P + ? E 2 ? w = ? W + ? P 2 ( d ? d x ) e = ? E ? ? P ( δ x ) e ( d ? d x ) w = ? P ? ? W ( δ x ) e
S?S ?通常为时间和 ??? ?的函数，为了简化，可以做如下线性处理：
S=Sc+Sp?pS = S c + S p ? p
式中， ScS c是常数， SpS p随时间和物理量 ??变化而变化
这样式(1)就可以进行化简了，带入！

(ρu)eAe?P+?E2?(ρu)wAw?W+?P2=ΓeAe[?E??P(δx)e]?ΓwAw[?P??W(δx)w]+(Sc+Sp?p)ΔV( ρ u ) e A e ? P + ? E 2 ? ( ρ u ) w A w ? W + ? P 2 = Γ e A e [ ? E ? ? P ( δ x ) e ] ? Γ w A w [ ? P ? ? W ( δ x ) w ] + ( S c + S p ? p ) Δ V
整理得：

(Γe(δx)eAe+Γw(δw)wAw?SpΔV)?P=(Γw(δx)wAw+(ρu)w2Aw)?W+(Γe(δe)eAe+(ρu)e2Ae)?E+ScΔV( Γ e ( δ x ) e A e + Γ w ( δ w ) w A w ? S p Δ V ) ? P = ( Γ w ( δ x ) w A w + ( ρ u ) w 2 A w ) ? W + ( Γ e ( δ e ) e A e + ( ρ u ) e 2 A e ) ? E + S c Δ V
可以记为：

aP?P=aW?W+aE?E+b(2)(2) a P ? P = a W ? W + a E ? E + b
因此只要把整个计算域的每个节点建立式(2)的离散方程，组成一个大的含有节点未知量的线性方程组，通过求解方程组，就可以得到 ??? ?在各节点处的值，至于如何求解方程组，将在后面介绍
2.3常见离散格式离散格式又叫插值方式，比如我们上一节介绍的线性插值(中心差分)就是一种很常见的离散格式
一维稳态无源项问题的控制方程

（dρu?）dx=ddx(Γd?dx)（ d ρ u ? ） d x = d d x ( Γ d ? d x )
该流动必须满足连续方程：
（dρu）dx=0（ d ρ u ） d x = 0

在控制体积P上积分控制方程和连续方程：

(ρuA?)e?(ρuA?)w=(ΓAd?dx)e?(ΓAd?dx)w(ρuA)e?(ρuA)w=0( ρ u A ? ) e ? ( ρ u A ? ) w = ( Γ A d ? d x ) e ? ( Γ A d ? d x ) w ( ρ u A ) e ? ( ρ u A ) w = 0
为了方便，定义两个新的物理量， FF代表通过界面处单位面积的对流质量通量， DD表示界面的扩散传导性(diffusion conductance)， F=ρuF = ρ u， D=Γ/(δx)D = Γ / ( δ x )
故，界面处的值可以表示为：
Fw=(ρu)w,Fe=(ρu)eF w = ( ρ u ) w , F e = ( ρ u ) e
Dw=Γ/(δx)w，De=Γ/(δx)eD w = Γ / ( δ x ) w ， D e = Γ / ( δ x ) e
同时定义一维单元的Peelet数 PeP e :
Pe=FD=ρUΓ/δxP e = F D = ρ U Γ / δ x
PeP e 表示对流和扩散强度之比，当 PeP e 为0时，问题就变为了纯扩散问题了，即流场中没有流动，只有扩散； Pe>0P e > 0 ，流体沿着正 xx 方向流动，当 PeP e 很大时，问题变为纯对流，扩散可以忽略
同时，可以假设 Ae=Aw=AA e = A w = A
利用中心差分，可以得到：
Fe?e?Fw?w=De(?E??P)?Dw(?P??W)F e ? e ? F w ? w = D e ( ? E ? ? P ) ? D w ( ? P ? ? W )
Fe?Fw=0F e ? F w = 0
这样我们需要求解的变量有 FF 和广义未知数 ?? 在界面 ee 和 ww 处的值， FF 的求解或者叫速度场的求解，我们放到以后讲，现在为了求解 ?? 不得不利用插值来表示
2.4中心差分格式利用中心差分
?e=?P+?E2? e = ? P + ? E 2
?w=?W+?P2? w = ? W + ? P 2
带入：

Fe2(?P+?E)?Fw2(?P+?W)=De(?E??P)?Dw(?P??W)F e 2 ( ? P + ? E ) ? F w 2 ( ? P + ? W ) = D e ( ? E ? ? P ) ? D w ( ? P ? ? W )
通过化简，依旧可以写成：
aP?P=aW?W+aE?E(3)(3) a P ? P = a W ? W + a E ? E
式中，
???????????aWaEaP=DW+Fw2=De?Fe2=aW+aE+(Fe?Fw){ a W = D W + F w 2 a E = D e ? F e 2 a P = a W + a E + ( F e ? F w )
注：

当 Pe<2P e < 2 时，中心差分格式的计算结果与精确解基本吻合
当 Pe<2P e < 2 时，中心差分格式失去物理意义，因为这时 aE<0a E < 0 ，会对结果产生振荡影响
系数 aEa E 和 aWa W 代表了邻点 EE 和 WW 处的物理量通过对流及扩散作用对 PP 点产生影响的大小，为了满足式(3),要求每一个系数都大于零，(正系数是考虑到后面的迭代)
要满足 Pe<2P e < 2 ,只能是速度非常小或者网格间距很小，
Pe=FDP e = F D ,也就是说要么F很小，要么D很大，对于一般的流动中心差分已经不能满足了，因此，必须引入新的离散格式，

2.5 一阶迎风格式经过之前的分析，我们知道 ?? 同时受到 ?W? W 和 ?P? P 的影响。在一个对流占据主导地位的由西向东的流动中， ww 界面受到来自节点 WW 应该比节点 PP 更强烈，举例子就是不能脚踏两只船，迎风格式在确定界面的 ?? 时考虑了流动方向
2.5.1一阶迎风数学描述

定义：因对流造成的界面上的 ?? 值被认为等于上游(迎风测节点)的 ?? 值，即当流动沿着正方向时( uw>0,ue>0,Fw>0,Fe>0u w > 0 , u e > 0 , F w > 0 , F e > 0 )，有：
?w=?W,?we=?P? w = ? W , ? w e = ? P

得：

Fe?P?Fw?W=De(?E??P)?Dw(?P??W)F e ? P ? F w ? W = D e ( ? E ? ? P ) ? D w ( ? P ? ? W )
由 Fe?Fw=0F e ? F w = 0，带入化简（这里对左边进行了技术处理，主要为了后面的计算）：
[(Dw+Fw)+De+(Fe?Fw)]?p=(Dw+Fw)?W+De?E[ ( D w + F w ) + D e + ( F e ? F w ) ] ? p = ( D w + F w ) ? W + D e ? E

定义：因对流造成的界面上的 ?? 值被认为等于上游(迎风测节点)的 ?? 值，即当流动沿着负方向时( uw<0,ue<0,Fw<0,Fe<0u w < 0 , u e < 0 , F w < 0 , F e < 0 )，有：
?w=?P,?we=?E? w = ? P , ? w e = ? E

得：

Fe?E?Fw?P=De(?E??P)?Dw(?P??W)F e ? E ? F w ? P = D e ( ? E ? ? P ) ? D w ( ? P ? ? W )
由 Fe?Fw=0F e ? F w = 0，带入化简（这里对左边进行了技术处理，主要为了后面的计算）：
[Dw+(De?Fe)+(Fe?Fw)]?p=Dw?W+(De?Fe)?E[ D w + ( D e ? F e ) + ( F e ? F w ) ] ? p = D w ? W + ( D e ? F e ) ? E
结合两个过程，可以得出以下结果
aP?P=aW?W+aE?E(4)(4) a P ? P = a W ? W + a E ? E
式中，
?????aWaEaP=DW+max(Fw,0)=De+max(0,?Fe)=aW+aE+(Fe?Fw){ a W = D W + max ( F w , 0 ) a E = D e + max ( 0 , ? F e ) a P = a W + a E + ( F e ? F w )
注：

1.一阶迎风格式考虑了流动方向的影响，其系数永远大于零，故不会引起解的振荡
2.缺点就是按界面上流动大于或者小于0而决定其取值，但实际上精确解还和 PeP e 数的大小有关
3.迎风格式不管 PeP e 数的大小，扩散项永远按中心差分计算，可是，当 |Pe|| P e | 巨大时，扩散作用接近于零，
此时迎风格式夸大可扩散项的影响，也就是说迎风格式在大的 |Pe|| P e | 条件下过高地估计了扩散值。
4.由于只有一阶精度，限制了其解的精度，除非有非常细密的网格，否则计算结果的误差较大
5.目前，一阶迎风常被二阶迎风所代替

2.6 混合格式综合中心差分和迎风的优点，在 |Pe|<2| P e | < 2 采用具有二阶精度的中心差分格式；当 |Pe|<≥2| P e | <≥ 2 采用一阶迎风
离散方程如下：

aP?P=aW?W+aE?E(4)(4) a P ? P = a W ? W + a E ? E
式中，
???????????????????aWaEaP=max[(Fw,(Dw+Fw2),0]=max[(Fw,(De?Fe2),0]=aW+aE+(Fe?Fw){ a W = max [ ( F w , ( D w + F w 2 ) , 0 ] a E = max [ ( F w , ( D e ? F e 2 ) , 0 ] a P = a W + a E + ( F e ? F w )
注：

1.混合格式根据 PeP e 数选择合适的格式，综合了我们之前讲的优点，因为其系数总是正的，故是无条件满的，
2.混合格式在CFD中广泛采纳，是非常实用的离散格式。
3.缺点还是精度了，但是！现在计算机那么好，加密网格怕什么！

2.7 指数形式

这个就不介绍了，对一维该方法是对的，到了高维度及源项不为0时，指数函数失去了意义，且指数计算相当费时

好了，写了这么多，总结一下，中心差分、一阶迎风和混合格式均属于低阶离散格式，这样难免会带来误差。由于一阶的截断误差小于二阶引起的误差，这个现象我们叫假扩散。有一点，必须指出，尤其一阶迎风人为的引入中心差分，在计算过程中可能使得扩散作用被认为放大，相当于引入了人工粘度或数值粘度
就物理过程本身而言，扩散使得物理量的变化率减小，使得流场处于均匀化。但是，在离散格式中，假扩散的存在却带来了计算结果偏离真实解了。