对流-扩散方程的离散格式
第五章对流扩散方程的离散格式
aP = aE + aW
aE = De – Fe / 2 aW = Dw + Fw / 2
在流场的实际求解过程中， 每一个迭代层次上，即使速度 场尚未收敛，也要保证连续方 程是满足的。
3. 对流项的中心差分与迎风差分
3.2 对流项的中心差分
三点说明：
系数 aE , aW 包含了对流 F 与扩散 D的作用的影响；
对均分网格：
2. 对流项离散格式的重要性及两种离散方式
2.2 构造对流项离散格式的两种方式 （2）控制容积积分
给出界面上被求函数的插值方式
2. 对流项离散格式的重要性及两种离散方式
2.2 构造对流项离散格式的两种方式 （3）两种定义之间的关系
对某种对流项离散格式，都可以用两种方法给出其相应 的定义；
3.1 一维对流-扩散问题模型方程的精确解
上游优势
3. 对流项的中心差分与迎风差分
3.1 一维对流-扩散问题模型方程的精确解
希望所构建的离散方程形式也具有这样的物理特性
3. 对流项的中心差分与迎风差分
3.2 对流项的中心差分 Central Scheme （CS）
分段线性
均分网格
令
对流项
----界面上的流量
1. 简 介
对流与扩散作用在物理本质上的区别
从物理过程来看，扩散作用与对流作用在传递信息或扰动方面 的特点有很大区别：
扩散是由分子的不规则热运动所致，分子不规则热运动对空间不同方向
的几率都是一样的，因而扩散过程可以把发生在某一地点上的扰动的影响 向各个方向传递。
对流是流体微团宏观的定向运动，带有强烈的方向性。在对流的作用下，
两种定义方式的截断误差阶数是一致的，均为二阶截差 （中心差分，分段线性）；
CFD大作业解析(对流扩散方程)
求解迁移方程
求解域为 初始条件为
求解当 t = 0.7 时，u 的值。
∂u ∂u ∂2u ∂t + ∂x = α ∂x2
x ∈ [0.0, 1.0]
{
u=
2.0, 0.1 < x < 0.2; 1.0, others.
2 迁移方程 α = 0 差分格式的构造
2.1 显式 Euler 法 + 空间后差
空间采取 Euler 后差。构造的差分格式如下
uni +1 − uni = − uni − uni−1
∆t
∆x
将其化为可求解的格式
uni +1
=
(1
−
∆t ∆x
)uni
+
∆t ∆x
uni−1
2.2 显式 Euler 法 + 空间前差
空间采取 Euler 前差。构造的差分格式如下
uni +1 − uni = − uni+1 − uni
∂u + ∂u = 2π cos(2πx) ∂t ∂x
边界条件为 x = 0, u = 0; x = 1, u = 0。 粘性项的验证，令 a = 0α = −1，令理论解为 u = x2 − x，则控制方程为
∂u ∂2u ∂t + ∂x2 = 2x
边界条件为 x = 0, u = 0; x = 1, u = 0。
1 N
∑精Ni=度1(测ui试−的Ui方)，法以为l空n N间, l离n e散为点坐数标N作图取，为计5算0、两1点00间、的20斜0，率进。行时间推进，直到达到理论解，计算误差
e
=
6 作业要求
• 编写程序，实现 2 中列出的四种方法，完成空间二阶导数的离散
哈尔滨工业大学 计算传热学 第五章 对流-扩散方程的离散格式-2013
aPP aEE aWW
Fe Fw exp( Pw ) aE , aW exp( Pe ) 1 exp( Pw ) 1
(D)
aP aE aW (Fe Fw )
区别就在函数 aE和aW
aE De
Pe aE De exp( Pe ) 1
aE Pe De
该格式计算量比指数小，且指数格式的解差别很小。
§ 5-3
为了在讨论中引入 PE 记
通用表达式
x
i
J*
i+1 i+1/2
x
1 界面i+ 上的值可以用界面两侧节点值表示 2
J * Bi Ai 1 （y）
系数A和B的性质的讨论 (1)当 i i 1 时，扩散量=0， J *完全由对流造成，即
即
aPP aEE aW W
显然不论那种格式，仅仅是 A(| P |) 表达式的区别。
A( P )
A(|P |)
中心 1 0.5 | P | 迎风 1 混合 [| 0,1 0.5 | P | |] 指数 | P | [exp(| P |) 1]
1.0
迎风
指数 乘方
乘方 | 0, (1 0.1| P |)5 |
中心
混合
P
§ 5-4
原始的假扩散概念
关于假扩散的讨论
一维非稳态对流方程（纯对流,没有扩散）
u t x
显示迎风差分格式
in1 in
t
u
in in 1
x
, o(x, t )
将上式在（i，n）点做Taylar级数展开，保留二阶。
上述若对任何成立，必得
B( P ) A( P ) A( P ) B( P )
对流扩散方程解析解
对流扩散方程解析解对流扩散方程（Convection-DiffusionEquation，CDE）是描述物理系统中物质扩散和热对流运动的方程。
它源于20世纪30年代真空磁体理论中发现的电子运动方程，在50年代被普及应用于各种工程、物理学和化学领域，如电子、热传输、水力学等，具有不可缺少的重要意义。
一般来说，对流扩散方程可以被描述为：$$frac{partial y}{partial t}=afrac{partial^2 y}{partial x^2}+bfrac{partial y}{partial x}+cfrac{partial y}{partial y}+d$$其中，a、b、c和d是常数，t和x分别代表时间和物理位置。
若把空间坐标投射到它们的平面上，则可以用更具体的形式表述为： $$frac{partial y}{partial t}=afrac{partial^2 y}{partial x^2}+bfrac{partial y}{partial x}+cfrac{partial y}{partial y}+d+frac{partial y}{partial z}$$其中，z是投射后的空间坐标，a、b、c和d也可以改变以适合不同的实际应用场景。
对于对流扩散方程的解析解，有两种基本方法：一种是用不定积分法；另一种是用微分平面法，也称作渐进分析方法。
从一般的原理上来看，不定积分法是把对流扩散方程拆解成多个简单的可求解的微分方程，然后分别求解它们，最后再综合求得总解。
此外，它还可以运用标准积分法来近似求解，特别有利于解复杂的多变量方程。
而渐进分析（Perturbation Analysis）是把复杂的问题划分成几个渐进步骤，每一步把问题简化为可以近似解决的状态，依此不断迭代，最终求得近似解。
这种技术通常用来求解非线性方程，对于对流扩散方程求解也非常有效，能有效地提高准确度和计算速度。
此外，还有其他一些求解方法，比如拉格朗日法（Lagrange Method）、拉普拉斯正则化（Laplace Regularization）以及偏微分方程的泛函理论方法（Functional Theory of Partial Differential Equations）等。
数值传热第六章作业
6-3 试在直角坐标系的交错网格上，写出动量离散方程式（6-5）、（6-6）中的系数nb a （即S N W E a a a a ,,,）,n n e e A a A a ,,,的表达式。
为简便起见，设（1）流体物性为常数；（2）在x, y 方向上网格各自均匀划分。
速度e u 的邻点可参阅图6-5， 速度n υ的邻点参见图6-32.对流、扩散项的离散可采用五种三点格式之一。
解：根据课本P145式(5-13)、(5-16)、(5-18)，对流、扩散项采用指数格式计算本题 在二维直角坐标系中，对流—扩散方程的通用形式为：()()()φφφφφρυφφρρφS y y x x y u x t +⎪⎪⎭⎫⎝⎛∂∂Γ∂∂+⎪⎭⎫ ⎝⎛∂∂Γ∂∂=∂∂+∂∂+∂∂ 对于动量方程，把压力梯度项放到源项中了。
引入在x 及y 方向的对流—扩散总通量密度，上式可改写为：()⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-=⎪⎪⎭⎫⎝⎛∂∂Γ-∂∂+⎪⎭⎫ ⎝⎛∂∂Γ-∂∂+∂∂y p x p S y y x u x t φρυφφφρρφφφ 即：()⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-=∂∂+∂∂+∂∂y p x p S y J x J t yx ρφ （1） 其中：yJ xu J y x ∂∂Γ-=∂∂Γ-=φρυφφφρφφ将（1）式对P 控制容积做时间与空间上的积分得：e E P p P c s n w e pp A P P V S S J J J J V t)()()()()()(0-+∆+=-+-+∆∆-φρφρφ将通用变量φ换成速度u ，相应的其控制容积变为：所以上式可改写为：e E P p p c s n w e ee A P P V S S J J J J V t u u )()()()()()(0-+∆+=-+-+∆∆-φρρ （2）式（6-5）为：()e E P nb nbe e A p p b u au a -++=∑对上式用界面总通量表达式为：ee E e e EE e u a u F a J -+=)( （3）e w W w W w u F a u a J )(--= （4）n N e n N n u a u F a J -+=)( （5）e s S s S s u F a u a J )(--= （6） 把以上方程代入方程（2）得：e E p e p c e s S s S n N en N w w e w W ee E e e EE ee A P P V u S S u F a u a u a u F a u a u F a u a u F a V tu u )()()()()()(0-+∆+=-+--++--+-++∆∆-ρρ整理得：eE p ec s S n N w W ee E ep s S n N w W e EE A P P u tVV S u a u a u a u a u V S F a F a F a F a tV)(])()()()([0-+∆∆+∆++++=∆--+++-+++∆∆ρρ当对流、扩散项的离散采用指数格式时， 则上式中的系数分别为：1)ex p()(-==∆∆e ee e EE P F P A D a 1)e x p ()e x p ()(-==∆∆∆w w w w w W P P F P B D a1)ex p()(-==∆∆n n n n N P F P A D a 1)e x p ()e x p ()(-==∆∆∆s s s s s S P P F P B D a tVa e ∆∆=ρ0V S a F F F F a a a a a p e s n w e S N W EE e ∆-+-+-++++=00e e c u a V S b +∆=y A e ∆=同理对（6－6）()n N P nb nbn n A p p b aa -++=∑υυ，类似地有：1)ex p()(-==∆∆n n n n NN P F P A D a 1)e x p ()e x p ()(-==∆∆∆s s s s s S P P F P B D a 1)ex p()(-==∆∆e e e e E P F P A D a 1)e x p ()e x p ()(-==∆∆∆ww w w w W P P F P B D atVa n ∆∆=ρ0V S a F F F F a a a a a p n s e w n S E W NN n ∆-+-+-++++=000n n c u a V S b +∆=x A n ∆=6-4 对图6-11所示的二维流动情形，已知：10,0,20,50====E N s w p p v u 流动是稳态的，且密度为常数。
2.4常用的离散格式
低阶格式的假扩散特性
迎风格式,指数格式,混合格式及乘方格式等 一阶格式应用于实际问题时都可能引起较严 重的假扩散,这在HVAC领域的高大空间流体 流动及传热计算中尤为明显. 因此,为了有效地克服或减轻假扩散所带来的 计算误差,空间导数应当采用二阶或更高阶的 格式(如QUICK格式,二阶迎风差分格式等).
离散格式
假设速度场已知，则为求解离散方程，需计算广义未 知量在边界e和w处的值。
为完成这一任务，必须决定界面物理量如何通过节点 物理量的插值表示。
各种不同的插值方法就构成了不同的离散格式。
中心差分格式
一阶迎风格式
混合格式
指数格式
乘方格式
1
2.4.1术语的约定
对离散格式的讨论以一维稳态对流扩散方程为例，不 涉及瞬态项。
3
Central differencing scheme 中心差分格式
(x) P
P
interpolated value
e E
eE
We determine the value of at the face by linear
interpolation between the cell centered values.就是界 面上的物理量采用线性插值公式来计算。
基于此限制，中心差分格式不能作为对于一般 流动问题的离散格式，需创建其它更合适的格 式（对纯扩散稳态，如热传导是适用的）。
5
对流扩散方程的精确解
6
精确解随Pe数的变化
（Pe＝0纯扩散，Pe增大对流增强）
7
具体算例
（不同计算工况意味着不同Pe数）
8
第一种工况Pe=0.2
尽管网格粗糙，但数值解与精确解非常接近。
第5章-对流-扩散方程的离散格式
uL
0
Pe表示对流与扩散作用 的相对大小。
0
4/59
传热与流体流动的数值计算
二、对流项的中心差分
d d d u 采用控制容积积分法 对方程 dx dx dx e u e w u w P 2 2 x w x e
aE De Fe ,0 , aW Dw Fw ,0
对流项一阶迎风：
aW i 1 aE i 1 P ,0 1 P ,0 P D D
12/59
传热与流体流动的数值计算
A P P
B P A P P A P P ,0 P B P A P P ,0
24/59
传热与流体流动的数值计算
四、aE、aW的通用表达式
* Je B Pe P A Pe E
J d J P D d x x
*
18/59
传热与流体流动的数值计算
一、通量密度及其离散表达式（续）
J*的离散表达式：
J * Bi Ai 1
Behind Ahead 界面后的项 界面前的项 以坐标轴正方向为依据的“前”、“后”。
19/59
传热与流体流动的数值计算
负系数会导致物理上不真实的解。
7/59
传热与流体流动的数值计算
三、对流项的迎风格式
Taylor展开法
d i i 1 , ui 0 dx i x
i 1 i , ui 0 x
控制容积积分法 e界面 ue 0 , P ; ue 0 , E w界面 uw 0 , W ; uw 0 , P
第五章对流扩散方程
• 比混合格式复杂，计算量增加，但准确性 提高
• 稳定性：若用中心差分格式不能体现对流 项的物理本质，常会引起数值解的振荡
• 经济性：若用高阶格式，无数值振荡，但 格式复杂，求解相对困难，机时消耗较多
5.2 一维稳态对流扩散问题
d (u)
dx
d dx
d
dx
5.2.1 模型方程的精确解
d (u)
dx
d dx
d
dx
边界条件：
x 0, 0；x L, L
采用迎风思想：从来流上游方向找依赖区
在界面e上,若 ue 0,则e P ;若 ue 0,则e E ; 在界面w上,若 uw 0,则w W ;若 uw 0,则w P
界面流量
• 引入符号 • 对流：
a1, a2 max(a1, a2 )
(u)e Fee P Fe,0 E Fe,0 , (u)w Fww W Fw,0 P Fw,0
• 控制方程变为： dJ 0; 或 J const dx
J
F
0
0 L
exp(Pe) 1
界面上通量
Jw
Fw
W
W P
exp(Pw )
1Hale Waihona Puke JeFeP
P
exp(
E
Pe )
1
Fe exp( exp(Pe
Pe )
) 1
Fw exp(Pw )
1 P
Fe exp(Pe )
1E
Fw exp(Pw exp(Pw )
) 1
W
合并整理结果
aPP aEE aWW
• 系数
aE
第五章 对流-扩散方程的离散格式
见下页表格：
5.3.5 5种三点格式系数计算式的汇总 不同格式离散方程的形式相同，但 系数不同。具体见下表5－1：
5.4 对流-扩散方程5种3点格式系 数特性的分析
5.4.1 通量密度及其离散表达式
J J ( / x )
*
由于 所以
d d J u [ P ] dx x d ( x / x)
i i 1 d ， ui 0 dx i x i 1 i x
ui 0
对多维问题，用此方法构造的对流 项的离散格式，只有在求解区域内 流速不发生逆向时，所形成的离散 方程才具有守恒性。
2、控制容积积分法定义
规定界面上的未知量恒取上游节点的值 e界面上： ue 0 , p ； ue 0 , E
把式（2）用于计算界面总通量密度Je, Jw: 对Je： ， , L x
0 P L E e
P E J e Fe [ P ] exp( Pe ) 1
对Jw：
0 W ， L P , L xw
W P J w Fw [W ] exp( Pw ) 1
对于坐标系I，C位于界面之后，而D位 于界面之前，于是： J * B( P )C A( P ) D 对于坐标系II，D位于界面之后，而C 位于界面之前，于是：
J B( P ) D A( P )C
*
由于
J J
*
*'
C [ B( P ) A( P )] D [ A( P ) B( P )]
exp( Pe ) 1
Fe ；
Fw exp( Pw ) aW exp( Pw ) 1
第六章对流与扩散
该格式计算量比指数小，且与指数格式的解差别很小。
§ 6-3
通用表达式
为了在讨论中引入 PE J* J u d * 记 J x x ( ) d( ) d i i+1 P d x i+1/2 d( ) 1 界面i+ 上的值可以有界面两侧节点值表示
第六章
对流扩散方程的差分格式
导热型方程：（原始或经过变换的）
二阶导数项（扩散），源项
对流扩散方程：（动量或能量）
二阶导数项（扩散），源项 一阶导数项（对流），压力梯度。 一维稳态无内热源的对流扩散方程:
d d d ( u ) ( ) 密度， 扩散系数。 dx dx dx
对流热能量方程
aE Pe De
aE 1 1 Pe De 2
指数
aE 0 De
二．混合格式
虽然指数格式是精确解,但计算过繁,通过对 随 Pe 变化及其三条切线 aE Pe 0 De aE Pe Pe De aE 1 Pe 0 1 Pe De 2 斯帕尔丁提出 aPP aEE aww
F u J * 而 P ,J D ( ) D x * 根据通量守恒 Je J De Je D J* 0
P{De B(P e ) D A(P )} De A(Pe )E e D BP W
aE De A(Pe ) De{A(| Pe |) [| Pe ,0 |]}
Pe 10
aE Pe DE aE (1 0.1Pe )5 Pe DE aE (1 0.1Pe )5 DE aE 0 DE
10 Pe 0
0 Pe 10
Pe 10
(f)
3 对流扩散方程的离散化(讲义)
令
dφ dx ρu dφ J* = φ− Γ d ( xδ ) δ J = ρ uφ − Γ
15 第三章 对流扩散方程的离散化 16
if Fw < 0 φ w = φ P ( ρu ) w φw = ( ρ u ) wφ P else φ w = φW ( ρ u ) w φw = ( ρ u ) w φW
第三章 对流扩散方程的离散化
3.1 一维对流扩散问题
• 指数格式
dφ dJ ⇒ = 0 ⇒ Je − J w = 0 dx dx • 将上面的精确解应用于P点和E点之间，得到指数格式 J = ρ uφ − Γ
第三章 对流扩散方程的离散化 17
aP = aE + aW + Fe − Fw
3
3.1 一维对流扩散问题
• 混合格式 由指数格式可知
6
3.1 一维对流扩散问题
• 采用分段近似
− Pe aE = 1 − Pe / 2 De 0 0 aW = 1 + Pw / 2 Dw Pw Pe<-2 -2<Pe<2 Pe>2 Pw<-2 -2<Pw<2 Pw>2
第三章 对流扩散方程的离散化
用图形表示的精确解为
1 1
0.8
0.8
( φ -φ 0 )/(φ (Px/L)-φ 0 )
0.6
0.6 图例
(φL − φ0 ) exp(
ρ uL ) Γ )
0.4
0.4
0.2
0.2
P=1 P=2 P=4 P=10 P=-10 P=-4 P=-2 p=-1 P=0
0 0 0.2 0.4 0.6 0.8 1
• 为了获得关于P点的离散化方程，必须将控制面上的控制变量用 节点上的值来表示。算术平均是最直接的
计算热物理复习提纲
第一章1、什么是计算热物理？《计算热物理》(也称《计算传热学》或《数值传热学》)是指用计算机通过数值方法来求解各类热量传递问题的一门《传热学》的分支学科。
2、做科研研究有哪几种方法？实验研究、数值研究、理论研究3、偏微分方程数值计算方法的本质？离散化、代数化第二章1、什么叫平衡问题？所谓平衡问题是一类单一的边值问题，其对应偏微分方程的定义域是一个封闭区域；定义域上每一点的解，依赖于封闭边界每一点上的边界条件。
注意：平衡问题不一定没有时间变量；平衡问题的判断是基于是否定义在封闭边界内，或者说是基于控制方程是否为椭圆型方程2、偏微分方程的物理分类及数学分类物理分类：平衡问题；行进问题（包括对流问题、扩散问题及对流扩散问题）；数学分类：椭圆形控制方程，双曲线控制方程，抛物性控制方程、联系物理实际：数学上的椭圆型方程对应物理上的平衡问题（且椭圆型方程的依赖区是整个封闭边界，而影响区是整个定义域）；抛物型方程对应物理上的扩散传播问题；双曲型方程对应物理上的对流传播问题；3、偏微分方程的适定性：是指偏微分方程的定解条件，能使方程解存在、解唯一、解连续的依赖它的初始或边界条件（即解稳定）4、定解条件（包括初始条件和边界条件）1)平衡问题（椭圆型方程）：必须提封闭边界上每一点的边界条件，要特别小心“无穷远边界”上的条件2）行进问题（抛物型、双曲型方程）：必须提初始条件（二阶偏微分方程需要有函数值、一阶导数值条件）；空间无界定义域问题，有的可以不提无界边界上的条件，但有界定义域问题，一般需要规定一定的边界条件。
初始条件：时间或类时间导数为一阶：只要给出函数在全域的初始值；时间或类时间导数为二阶：函数在全域的初始值，以及函数在全域对时间或类时间变量的一阶导数值。
边界条件：第一类边界条件：直接规定边界上的函数值（可以随时间或类时间变化）；第二类边界条件：直接规定在边界上的函数导数值，热物理问题一般规定：热流值在流入边界内部方向时为正；第三类边界条件：规定边界上的函数值与它的法向导数之间的某个关系；第三章1、离散格式的迁移性是指对离散方程中的对流项的某种离散格式，仅能使其扰动顺着流动方向传递的特性，又称为格式的迎风性或传输性。
对流扩散问题的有限体积法
流体仿真与应用第八讲二、对流-扩散问题的有限体积法◆中心差分格式（例子）节点增加到20个结果◆离散格式的性质在数学上，一个离散格式必须要引起很小的误差（包括离散误差和舍入误差）才能收敛于精确解，即要求离散格式必须要稳定或网格必须满足稳定性条件。
在物理上，离散格式所计算出的解必须要有物理意义，对于得到物理上不真实的解的离散方程，其数学上精度再高也没有价值。
通常，离散方程的误差都是因离散而引起，当网格步长无限小时，各种误差都会消失。
然而，在实际计算中，考虑到经济性（计算时间和所占的内存）都只能用有限个控制容积进行离散。
因此，格式需要满足一定的物理性质，计算结果才能令人满意。
主要的物理性质包括：守恒性、有界性和迁移性。
◆离散格式的性质——守恒性满足守恒性的离散方程不仅使计算结果与原问题在物理上保持一致，而且还可以使对任意体积（由许多个控制容积构成的计算区域）的计算结果具有对计算区域取单个控制容积上的格式所估计的误差。
◆离散格式的性质——迁移性③当Pe 为有限大小时，对流和扩散同时影响一个节点的上、下游相邻节点。
随着Pe 的增加，下游受的影响逐渐增大，而上游受的影响逐渐变小。
①，即纯扩散，无对流。
②，即纯对流，无扩散。
0=Pe ∞=Pe◆迎风格式迎风格式（Upwind Differencing Scheme ）在确定控制容积界面上的值时就考虑了流动的方向性，其思想为：在控制容积界面上对流项的取上游节点处的值，称之为第二类迎风格式。
中心差分格式的缺点是，它不能识别流动的方向，控制容积界面上的值取相邻上、下游节点的平均值。
当对流作用较强时，这样的处理就与其物理特征（某点的值受上游的影响，而不受下游的影响）不一致了。
φφφ◆迎风格式◆迎风格式在控制容积界面上对流项的取其上游节点处的值EW →φWw φφ=Pe φφ=()()W P w P E e W w P e D D F F φφφφφφ−−−=−()()[]()Ee W w w P w e e w w D F D F F D F D φφφ++=−+++WE →Pw φφ=Ee φφ=()()[]()Ee e W w Pw e e e w F D D F F F D D φφφ−+=−+−+◆迎风格式通用形式WW E E P P a a a φφφ+=()w e E W P F F a a a −++=EW →ww W F D a +=eE D a =W E →w W D a =ee E F D a −=◆迎风格式的特点迎风格式满足守恒性。
对流扩散方程背景
对流扩散方程背景提出一种隐格式用于解决二维时间依赖的Burgers型方程。
迎风单边差分格式被用于对流项离散；对扩散项用二阶中心差分格式离散。
我们建立了全隐的数值有限差分格式，分析了无条件稳定性和严格推导了收敛性，在空间是二阶收敛的和时间一阶收敛的。
给出数值结果验证理论正确性。
关键词：隐格式，单边差分逼近，Burgers方程，稳定性，收敛阶对流扩散方程背景对流扩散方程描述黏性流体的动力学行为，这在许多工程应用中发挥了重要作用。
对流占优型扩散方程一般具有对流比扩散的系数大得多的特点，通常数值模拟具有一定难度，因为一方面，扩散系数比传输速度小，并且在另一方面，由于数值扰动容易出现边界层现象。
许多格式已用于这些问题的模拟，并有大量成功的数值方法[1-3]。
通过离散方法来解决对流扩散问题时，一般运用标准Galerkin有限元方法求解，但此方法会导致非物理特征扰动。
为了解决这类缺陷，几种稳定的有限元方法已经在[4]中被提出了。
我们感兴趣的是建立非耗散方法来克服数值扰动，并有鲁棒性和二阶精度，尤其是对Burgers问题。
Burgers问题通常被认为非线性流体的流动和扰动的经典模型。
在二维非线性的情况下，可以描述对流和扩散的现象，Burgers方程代表一种最基本的非线性模型方程。
从一个数值格式出发研究是相当有趣的，因为Burgers已出现在众多的流体方程中[5-7]。
并已经由霍普夫-科尔计算出多种组合的初始条件和边界条件下的结果[10,11]。
此外，对于非线性Burgers方程的解析解也可以通过Homotop Perturbation法[12]得到。
众所周知，单独的选择一种基本差分格式如中心差分或者迎风格式，来计算纯对流式的方程，扩散项通常只是中心近似。
而解决问题的关键在于对流方面构造稳定的离散结构来克服数值扰动。
虽然单边差分近似格式已经被提出了30年之久[13]，人们却很少关注他们在计算流动问题。
一阶或者二阶单边迎风有限差分近似已经在[14]中用来分析双曲型偏微分方程的数值解，尤其在一维空间非线性守恒定律上。
数值传热第五章课件2陶文铨
主讲陶文铨西安交通大学能源与动力工程学院热流中心CFD-NHT-EHT CENTER2010年10月18日, 西安数值传热学第五章对流扩散方程的离散格式(2)对流项离散格式的重要性及两种离散方式5.5.1假扩散的含义与成因5.5.2一阶截差格式引起严重假扩散举例1.本来的含义2.扩充的含义3.Taylor 展开法的分析5.5关于假扩散的讨论5.5.3网格倾斜交叉引起的计算误差5.5.4 非常数源项引起的假扩散5.5.5 两个名例以一维非稳态纯对流过程为例俩分析，其中有两n nφφ2(,O x φΔΔ其中关于时间的二阶导数项可做如下变化：时才没有这部分的计算误差。
2. 扩充的含义现有文献中常常将较大的计算误差都称为假扩散，大致有以下几项原因：(1) 一阶导数的一阶截差格式；(2) 流动方向与网格线呈倾斜交叉；(3) 离散格式未计及非常数源项的影响。
5.5.2一阶截差格式引起严重假扩散举例1.一维稳态对流扩散问题对流项用FUD，扩散项用CD，当Pe较大时，数值计算结果严重偏离精确解。
Physically plausible solution纯对流传递纯对流传递由离散方程：1n−1此时只有对流，没有扩散！时则有严重假扩散！0.8C =0.8C =当时，产生了严重的扩散作此种误差称为流向假扩散Γ≠Γ气流01. 设UE对P 控制容积，有2. 设控制容积，此时：计算误差纯对流传递三个对流问题的归纳这就是假扩散纯对流传递3)网格倾斜交叉引起的计算误差E冷热流体之间产生了温度均匀化的过程，即交叉5.5.5 已知流场计算温度场232(1),2(1)u y x v x y =−=−−参考解xT严重假扩散2) Leonard细高方腔中的自然对流换热5.6.1采用高阶格式克服流向假扩散5.6可以克服或减轻假扩散的格式与方法5.2.2 克服、减轻交叉假扩散的方法1. 采用二阶迎风2.采用三阶迎风3. 采用QUICK 格式1. 采用有效扩散系数2.采用自适应网格4. 采用SGSD 格式可以克服或减轻假扩散的格式与方法相当于界面上的中心差分)W WWxφ+Δ如型线上凹，则(2) FVM向上游取两点定义界面插值2.采用三阶迎风展开定义－一阶导数的三阶偏差分格式3. 采用定义－界面的插值在中心差分基础上考虑曲中心差分插值率修正？需要满足两个条件：插值的正确修正：相邻(2)0W PE φφφ−+<型线下凹8Cur −对e-界面u e 小于零时，取,,W P φφφu e 大于零时，取怎样相邻的三点?QUICK(2)e φφ=1/2w i φφ−=有：4. 采用CD条件稳定，但没有二阶假扩散；二阶迎风绝对稳定，组合起来，但是：如何确定值，特别是如何由计算结果来5. 高阶格式实施中的问题f u f计算边界：固o2) 代数方程的求解：等时，5.6.2用减小扩散系采用自适应网格（以减轻流5.7 对流-扩散方程离散形式稳定性分析5.7.1 数值计算中常见的三种不稳定性5.7.2 分析对流项格式不稳定性的“符号不变原则”5.7.3 稳定性分析结果讨论5.7.4 对流项格式问题讨论小结2.“符号不变”原则的基本思想3. “符号不变”原则的实施步骤4. “符号不变”原则的实施例子1. 研究背景扩散方程离散形式稳定性分析也会产生振荡的解，称为对流项离散格式的不稳态定性，研究目的是，找出产生振荡的临界Peclet 数。
对流扩散方程数值方法
4
情况
E100,W50
E50,W100
P 25 P 125
适用条件有限
➢ 当P 2 时，系数变成负数，会产生振荡
➢ 小适，用或条只件能P处理2 流速这比意较味小着的：问网题格步长必须很
5.2.3 一阶迎风格式
➢ 1.积分离散和格式推演
d(u)
dx
ddx
d
dx
(u)e(u)w d d x e d d x w
5.2.4 指数格式
➢ 1．指数格式推演，根据精确解表达式
L 0 0 e e x x p p ( ( u u L x// ) ) 1 1 e x p e ( x P p e (P x e ) / L ) 1 1
对流-扩散总通量密度
J u d
dx
L 0 0 e e x x p p ( ( u u L x// ) ) 1 1 e x p e ( x P p e (P x e ) / L ) 1 1
1Pe 2
,
0
当 P e2 aED e0
同理可得
aWDwPw,
1Pw, 2
0
2．混合格式评述
➢ 当Pe 2 时，与中心差分格式相同 ➢ 当Pe 2 时，为扩散项=零的迎风格式 ➢ 形式简单，兼有中心差分和迎风格式的优点，
较为适用
5.2.6 乘方律格式
➢ 1．乘方律格式提出和构成[Patankar] 四条包络线
a P a E a W (F e F w )
2. 一阶迎风格式评述
➢ 符合对流项物理意义，离散系数恒>零，无数 值振荡。得到广泛的应用。
➢ 迎风思想：二阶迎风、三阶迎风和QUICK格式 都很好地吸取了迎风思想：让信息主要来自流 动上游方向
数值分析-对流项的离散
数值分析-对流项的离散对流项的数值离散扩散项的离散在先前已经介绍，对于动量⽅程，更为棘⼿的问题在于N-S⽅程中对流项和压⼒项的离散处理。
在不可压流场的数值求解过程中，最重要的问题均是这两个⼀阶偏导项所引起的。
这⾥主要讨论关于对流-扩散项的离散⽅法。
对流项的离散⽅法⾸先是关于对流项的两种离散思路的介绍Taylor展开离散Taylor展开的⽅法，就是直接套⽤泰勒公式将偏导转化为线性离散控制容积积分通过将对流项的⼀阶导数在控制容积P中进⾏积分，利⽤相邻两个结点的值来获得控制容积边界e，w处的插值两种思路所给出的截断误差的阶数是相同的⼀维稳态对流-扩散的离散由于整个N-S⽅程还是过于复杂，这⾥我们对整个物理问题可以进⾏相应的简化，从仅含有对流项和扩散的场景出发，来讨论如何解决对流项中的⼀阶导离散问题。
因此，⾸先对于⼀维稳态⽆内热源的对流扩散问题的数学描述，我们可以很轻松的给出⽅程的守恒形式：ddx(ρuϕ)=ddx(Γdϕdx)边界条件写为：ϕ(0)=ϕ0；ϕ(L)=ϕL值得注意的是，以下内容表⽰的是对对流项（convection term）采⽤不同的差分格式进⾏描述，⽽公式中的扩散项仍然保持中⼼差分格式。
中⼼差分格式采⽤控制容积法，对流项的中⼼差分格式相当于在控制容积P的界⾯上取分段的线性函数，即(ρuϕ)e−(ρuϕ)w=(ρu)e ϕE+ϕP2−(ρu)wϕP+ϕW2最终的离散结果可以写为：ϕP[(ρu)e2−(ρu)w2+Γe(δx)e+Γw(δx)w]=ϕE(Γe(δx)e−(ρu)e2)+ϕW(Γw(δx)w+(ρu)w2)这⾥可以将界⾯上的流量记为$ F = \rho u，界⾯上单位扩散阻⼒的倒数记为D = \frac{\Gamma}{\delta x}$最终式(2)可以化简为：a PϕP=a EϕE+a WϕWpython代码可以轻松实现中⼼差分算法的描述：# 中⼼差分格式（x结点数，y结点数，x⽅向扩散，y⽅向扩散，x界⾯流量，y界⾯流量，linux下debug符号）def CDiff(x_nodes,y_nodes,D_x, D_y, F_x, F_y, Boundary, debug="n"):"""Requires:x_nodes = number of nodes along xy_nodes = number of nodes along yD_x = Diffusivity along x, type = floatD_y = Diffusivity along y, type = floatF_x = Force along X, type = matrix (x by y)F_y = Force along Y, type = matrix (x by y)Optional:debug, "y" if you want debug, off by default"""## setting up stuff (simple math)dim = x_nodes * y_nodes # specified to reduce calculationsw_bound = [0] * y_nodes # creating 0-matricies to finde_bound = [0] * y_nodes # all boundary nodess_bound = [0] * x_nodesn_bound = [0] * x_nodesif debug == 'y': # Debugger notifies progressprint("finding boundaries")for i in range (x_nodes): # boundary nodes are populateds_bound[i] = i # used to populate matriciesn_bound[i] = (dim - x_nodes + i)for i in range (y_nodes):w_bound[i] = (i * x_nodes)e_bound[i] = (i+1) * x_nodes -1Processing math: 100%if debug == 'y': # Debugger prints boundary nodesxB = "n:" + str(n_bound) + "\n"+ " s:" + str(s_bound)yB = "e:" + str(e_bound) + "\n"+ " w:" + str(w_bound)print("boundaries: \n", xB, "\n", yB)# Declaring variablesa_w = [0] * dim # Creating empty matricies for all nodesa_e = [0] * dima_s = [0] * dima_n = [0] * dimsp_p = 0su_p = 0B = [0] * dimsp_W = [0] * dimsu_W = [0] * dimsp_E = [0] * dimsu_E = [0] * dimsp_S = [0] * dimsu_S = [0] * dimsp_N = [0] * dimsu_N = [0] * dim#Boundariesphi_W = Boundary[0]phi_E = Boundary[1]phi_S = Boundary[2]phi_N = Boundary[3]## answer matricies ##a = [0] * dim # make afor i in range(dim): # 0-arraya[i] = [0] * dimB = [0] * dim# POPULATING USING CENTRAL DIFFERENCEprint("populating central difference")pbar = tqdm(range(dim)) # load progress barfor i in range(dim):pbar.update(1) # progress indicatorif i not in w_bound: # populating non-boundary conditionsa[i][i - 1] = a[i][i - 1] - (D_x + (F_x[i - 1] / 2))a[i][i] = a[i][i] + (D_x + (F_x[i - 1] / 2))if i not in e_bound:a[i][i + 1] = a[i][i + 1] - (D_x - (F_x[i + 1] / 2))a[i][i] = a[i][i] + (D_x - (F_x[i + 1] / 2))if i not in s_bound:a[i][i - x_nodes] = a[i][i - x_nodes] - (D_y + (F_y[i - x_nodes] / 2))a[i][i] = a[i][i] + (D_y + (F_y[i - x_nodes] / 2))if i not in n_bound:a[i][i + x_nodes] = a[i][i + x_nodes] - (D_y - (F_y[i + x_nodes] / 2))a[i][i] = a[i][i] + (D_y - (F_y[i + x_nodes] / 2))if i in w_bound: # populating boundary conditionssp_W[i] = -(2*D_x + F_x[i])su_W[i] = (2*D_x + F_x[i])* phi_Wa[i][i] = a[i][i] - sp_W[i]B[i] = B[i] + su_W[i]if i in e_bound:sp_E[i] = -(2*D_x - F_x[i])su_E[i] = (2*D_x - F_x[i])* phi_Ea[i][i] = a[i][i] - sp_E[i]B[i] = B[i] + su_E[i]if i in s_bound:sp_S[i] = -(2*D_y + F_y[i])su_S[i] = (2*D_y + F_y[i])* phi_Sa[i][i] = a[i][i] - sp_S[i]B[i] = B[i] + su_S[i]if i in n_bound:sp_N[i] = -(2*D_y - F_y[i])su_N[i] = (2*D_y - F_y[i]) * phi_Na[i][i] = a[i][i] - sp_N[i]B[i] = B[i] + su_N[i]return a, B迎风差分格式利⽤控制容积法，对对流项的⼀阶迎风差分格式进⾏描述：在e界⾯上u e>0,ϕ=ϕP u e<0,ϕ=ϕE 在w界⾯上u w>0,ϕ=ϕW u w<0,ϕ=ϕP 因此对流迎风差分格式可以表达为：(ρuϕ)e=F eϕe=ϕP MAX[F e,0]−ϕE MAX[−F e,0]最终的离散结果可写为：ϕP[MAX[F e,0]+MAX[−F w,0]+Γe(δx)e+Γw(δx)w]=ϕE(Γe(δx)e+MAX[−F e,0])+ϕW(Γw(δx)w+MAX[F w,0])需要注意的是：在数值解不出现振荡的范围内，中⼼差分要⽐迎风差分结果的误差⼩⼀阶迎风格式在应⽤上误差较⼤，应⽤需由限制，但新发展出的⼆阶迎风、三阶迎风都从中吸取了思想python代码中同样也可以实现⼀阶迎风差分算法的描述：# ⼀阶迎风差分格式（x结点数，y结点数，x⽅向扩散，y⽅向扩散，x界⾯流量，y界⾯流量，linux下debug符号）def UDiff(x_nodes,y_nodes,D_x, D_y, F_x, F_y, Boundary, debug="n"):"""Requires:x_nodes = number of nodes along xy_nodes = number of nodes along yD_x = Diffusivity along x, type = floatD_y = Diffusivity along y, type = floatF_x = Force along X, type = matrix (x by y)F_y = Force along Y, type = matrix (x by y)Optional:debug, "y" if you want debug, off by default"""## setting up stuff (simple math)dim = x_nodes * y_nodes # specified to reduce calculationsw_bound = [0] * y_nodes # creating 0-matricies to finde_bound = [0] * y_nodes # all boundary nodess_bound = [0] * x_nodesn_bound = [0] * x_nodesif debug == 'y': # Debugger notifies progressprint("finding boundaries")for i in range (x_nodes): # boundary nodes are populateds_bound[i] = i # used to populate matriciesn_bound[i] = (dim - x_nodes + i)for i in range (y_nodes):w_bound[i] = (i * x_nodes)e_bound[i] = (i+1) * x_nodes -1if debug == 'y': # Debugger prints boundary nodesxB = "n:" + str(n_bound) + "\n"+ " s:" + str(s_bound)yB = "e:" + str(e_bound) + "\n"+ " w:" + str(w_bound)print("boundaries: \n", xB, "\n", yB)dim = x_nodes* y_nodesa_w = [0] * dim # Creating empty matricies for all nodesa_e = [0] * dima_s = [0] * dima_n = [0] * dimsp_p = 0su_p = 0B = [0] * dimsp_W = [0] * dimsu_W = [0] * dimsp_E = [0] * dimsu_E = [0] * dimsp_S = [0] * dimsu_S = [0] * dimsp_N = [0] * dimsu_N = [0] * dim#Boundariesphi_W = Boundary[0]phi_E = Boundary[1]phi_S = Boundary[2]phi_N = Boundary[3]## answer matricies ##a = [0] * dim # make afor i in range(dim): # 0-arraya[i] = [0] * dimB = [0] * dim# POPULATING USING CENTRAL DIFFERENCEprint("populating upwind difference")pbar = tqdm(range(dim))for i in range(dim):pbar.update(1)if i not in w_bound: # populating non-boundary nodesa[i][i - 1] = a[i][i - 1] - (D_x + max((F_x[i - 1]),0))a[i][i] = a[i][i] + (D_x + max((F_x[i - 1]),0))if i not in e_bound:a[i][i + 1] = a[i][i + 1] - (D_x + max(-(F_x[i + 1]),0))a[i][i] = a[i][i] + (D_x + max(-(F_x[i + 1]),0))if i not in s_bound:a[i][i - x_nodes] = a[i][i - x_nodes] - (D_y + max((F_y[i - x_nodes]),0))a[i][i] = a[i][i] + (D_y + max((F_y[i - x_nodes]),0))if i not in n_bound:a[i][i + x_nodes] = a[i][i + x_nodes] - (D_y + max(-(F_y[i + x_nodes]),0)) a[i][i] = a[i][i] + (D_y + max(-(F_y[i + x_nodes]),0))if i in w_bound: # populating boundary nodessp_W[i] = -(2*D_x + max((F_x[i]),0))su_W[i] = (2*D_x + max((F_x[i]),0))* phi_Wa[i][i] = a[i][i] - sp_W[i]B[i] = B[i] + su_W[i]if i in e_bound:sp_E[i] = -(2*D_x + max(-(F_x[i]),0))su_E[i] = (2*D_x + max(-(F_x[i]),0))* phi_Ea[i][i] = a[i][i] - sp_E[i]B[i] = B[i] + su_E[i]if i in s_bound:sp_S[i] = -(2*D_y + max((F_y[i]),0))su_S[i] = (2*D_y + max((F_y[i]),0))* phi_Sa[i][i] = a[i][i] - sp_S[i]B[i] = B[i] + su_S[i]if i in n_bound:sp_N[i] = -(2*D_y + max(-(F_y[i]),0))su_N[i] = (2*D_x + max(-(F_y[i]),0))* phi_Na[i][i] = a[i][i] - sp_N[i]B[i] = B[i] + su_N[i]return a, B混合格式乘⽅格式指数格式假扩散。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第5章 对流-扩散方程的离散格式
2009年3月13日
1/59
传热与流体流动的数值计算
§5.1 对流项离散格式的重要性 及两种离散方式
一、对流项离散格式的重要性
1、数值解的准确性（假扩散） 2、数值解的稳定性 3、数值解的经济性
二、构造离散格式的两种方式
1、Taylor展开法 2、控制容积积分法
u e
Fee
P
max Fe,0
E
max Fe,0
w界面
P Fe ,0 E Fe ,0
uw 0 , W ; uw 0 , P
u w
Fww
W
max Fw,0
P
max Fw,0
W Fw,0 P Fw,0
9/59
传热与流体流动的数值计算
三、对流项的迎风格式（续）
3、一阶迎风格式截差阶数低，除非采用相当密的网格， 否则计算结果的误差较大。
4、一阶迎风格式的启示：应当在迎风方向取更多的信 息构造格式，更好地反映对流过程的物理本质。
5、在调试程序或计算的中间过程仍可以采用一阶迎风 格式。
11/59
传热与流体流动的数值计算
§5.3 对流-扩散方程的混合格式及乘方格式
一、通量密度及其离散表达式
d dx
u
d dx
d
dx
总通量密度J：单位时间内、单位面积上由扩散
及对流作用而引起的某一物理量的总转移量。
J
u
d
dx
x
P
d
d
x
x
J*
J D
P
d
d
x
x
18/59
传热与流体流动的数值计算
一、通量密度及其离散表达式（续）
一、系数aE与aW 之间的内在联系
aE(i)与aW (i+1)共享同一个界面。
对流项中心差分：
aE
De
Fe 2
, aW
Dw
Fw 2
aW
i 1
D
aE i
D
1
P 2
1
P 2
P
对流项一阶迎风：
aE De Fe,0 , aW Dw Fw,0
aW i 1 aE i
D
D
1 P,0
Fe 2
Dw
Fw 2
P
De
Fe 2
E
Dw
Fw 2
W
aPP aEE aWW
aP aE aW Fe Fw
,
aE
De
Fe 2
,
aW
Dw
Fw 2
在数值计算过程中，如果连续性方程始终得到满足，
则： aP aE aW
在求解过程中，始终保持连续性方程满足非常重要。
常物性条件下均分网格：
五、5种3点格式系数汇总只需给出
aE De
定义式
格式
定义
中心差分
1 Pe 2
迎风格式
1 Pe,0
混合格式 乘方格式 指数格式
Pe , 1 0.5Pe , 0
0
,
1 0.1 Pe
5
+
0
,
Pe
Pe
exp Pe 1
17/59
传热与流体流动的数值计算
§5.4 对流-扩散方程5种3点格式系数特性的分析
Fe exp exp Pe
Pe
1
exp
Fw Pw
1
P
exp
Fe Pe
1E
Fw exp Pw exp Pw
1
W
aPP aEE aWW
aE
Fe
exp Pe 1
,
aW
Fw exp Pw exp Pw 1
aP aE aW Fe Fw
14/59
传热与流体流动的数值计算
d
dx
采用控制容积积分法
e
x e
u e 2
w
x w
u w 2
P
e
x e
u e 2
E
w
x
界面的流量。
D= x 界面上单位面积扩散阻力的倒数(扩导)。
F D
=
u x
u
x
P
5/59
传热与流体流动的数值计算
二、对流项的中心差分（续）
De
迎风格式离散形式： aPP aEE aWW
aE De Fe, 0 aW Dw Fw, 0
aP aE aW Fe Fw
10/59
传热与流体流动的数值计算
四、中心差分与一阶迎风格式的讨论
1、对流项中心差分在不发生振荡的参数范围内，比一 阶迎风格式的误差更小。
2、一阶迎风格式离散方程系数永远大于零，不会引起 解的振荡，得到物理上看似合理的解。
1 P,0
P
12/59
传热与流体流动的数值计算
二、混合格式(Spalding,1971)
0
aE De
1 0.5Pe Pe
, Pe 2 , 2 Pe 2 , Pe 2
aE De
Pe
,
1 0.5Pe
,
0
13/59
传热与流体流动的数值计算
三、指数格式
利用精确解得到相邻节点间符合精确解的关系式。
三、指数格式（续）
15/59
传热与流体流动的数值计算
四、乘方格式(Patankar,1979)
0
, Pe 10
aE De
1 0.1Pe 5
1
0.1Pe
5
Pe
, ,
0 Pe 10 10 Pe 0
Pe
, Pe 10
aE
De
0
,
1 0.1 Pe
5
+
0
,
Pe
16/59
传热与流体流动的数值计算
3/59
传热与流体流动的数值计算
一、一维对流-扩散问题模型方程的精确解（续）
0 L 0
eux euL
1 1
ePex L 1 ePe 1
Peclet数：
Pe uL
0
Pe表示对流与扩散作用 的相对大小。
0
4/59
传热与流体流动的数值计算
二、对流项的中心差分
对方程
d u
dx
d dx
P
1
0.5P
E
2
1
0.5P
W
6/59
传热与流体流动的数值计算
二、对流项的中心差分（续）
例：在一维模型方程离散求解的
均分网格中，已知W =100， E =200。试对P =0，1，2及4
四种情况按中心差分格式计算
P之值。
负系数会导致物理上不真实的解。
7/59
传热与流体流动的数值计算
三、对流项的迎风格式
两种定义截差阶数一致，但截差首项系数有所不同。
2/59
传热与流体流动的数值计算
§5.2 对流项的中心差分与迎风格式
一、一维对流-扩散问题模型方程的精确解
d dx
u
d dx
d
dx
边界条件： x 0 , 0 ; x L , L
d udx C
ln
C1
ux
C2
C1eux C2
Taylor展开法
d i i1 dx i x
,
ui 0
i1 i x
, ui 0
控制容积积分法 e界面 ue 0 , P ; ue 0 , E w界面 uw 0 , W ; uw 0 , P
8/59
传热与流体流动的数值计算
三、对流项的迎风格式（续）
e界面
ue 0 , P ; ue 0 , E