对流扩散方程的离散格式
对流扩散问题有限体积法
流体仿真与应用第八讲二、对流-扩散问题的有限体积法◆中心差分格式（例子）节点增加到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 −=◆迎风格式的特点迎风格式满足守恒性。
第五章对流扩散方程的离散格式
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 )
数值传热第六章作业
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 流动是稳态的，且密度为常数。
输运方程对流扩散方程
输运方程对流扩散方程输运方程是描述物质传输过程的数学模型，常见的有对流扩散方程。
对流扩散方程是由对流和扩散两种机制共同产生的输运过程来描述的，它的一般形式为：∂c/∂t+∇·(v*c)=∇·(D*∇c)其中，c表示物质的浓度或者响应变量，t表示时间，v表示流体的速度场，D表示物质的扩散系数，∇表示梯度运算符。
对流项描述了物质的对流运动，即物质随着流体的移动而移动。
对于三维坐标系来说，对流项可以表示为∇·(v*c)。
具体来说，对流项的每一项分别表示了物质在x、y和z方向上的携带速度与浓度梯度的乘积。
扩散项描述了物质由浓度高处至浓度低处的扩散现象，即物质自发性地从高浓度区域向低浓度区域传播。
扩散项可以表示为∇·(D*∇c)，其中D是扩散系数，表示物质扩散的速率与浓度梯度的乘积。
对流扩散方程的物理意义是描述了物质在流体中传输的速率与物质浓度梯度之间的关系。
通过对流项，方程能够描述物质随着流体的运动快速传输的现象；而通过扩散项，方程能够描述物质由浓度高处向浓度低处传输的现象。
综合考虑对流和扩散的作用，对流扩散方程能够比较准确地描述物质在流体中的传输过程。
对流扩散方程在科学和工程领域有广泛的应用。
例如，在污染物传输和扩散模拟中，对流扩散方程可用于描述污染物由源区到周围空气或水体的传输过程。
在热传导模拟中，对流扩散方程可用于描述热量由高温区域到低温区域的传导过程。
在物质传递过程中，对流扩散方程也被广泛应用于描绘物质的传输行为。
总结起来，对流扩散方程是一种常见的输运方程，它能够描述物质由流体传输并扩散的过程。
通过对流项和扩散项的综合作用，对流扩散方程能够比较准确地描述物质在流体中的传输行为，所以在科学和工程领域有着广泛的应用。
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章 对流-扩散方程的离散格式
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
传热与流体流动的数值计算
数值分析-对流项的离散
数值分析-对流项的离散对流项的数值离散扩散项的离散在先前已经介绍，对于动量⽅程，更为棘⼿的问题在于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混合格式乘⽅格式指数格式假扩散。
《对流扩散方程》课件
环境科学
描述污染物在大气、水体等环境 介质中的扩散、输移和归宿。
在环境科学中，对流扩散方程用 于模拟污染物在大气、水体等环 境介质中的扩散、输移和归宿过
程。
在环境保护、污染治理等领域， 对流扩散方程具有重要的应用价
值。
化学反应动力学
描述化学反应在流体或固定床 反应器中的传递和反应过程。
在化学反应动力学中，对流 扩散方程用于模拟化学反应 在流体或固定床反应器中的
初始条件
指定在求解开始时刻的解的性质，如 常数、函数等。
03 对流扩散方程的应用
流体动力学
01
描述流体在运动状态下的物质传递和扩散现象。
02
在流体动力学中，对流扩散方程用于模拟流体中的物质传递过
程，如温度、浓度、速度等。
在航空航天、船舶、汽车等领域的流体动力学分析中，对流扩
03
散方程被广泛应用。
应用于多尺度问题
研究对流扩散方程在多尺度问题中的应用，如 微纳尺度流动、大气污染扩散等。
探索新的应用领域
将该方程应用于其他领域，如生物医学、环境科学等。
与其他领域的交叉研究
与流体动力学结合
研究对流扩散方程与流体动力学之间的相互 作用和影响，探索更深入的物理机制。
与偏微分方程理论的交叉
将对流扩散方程的研究与偏微分方程理论相 结合，推动数学理论的发展。
02
03
有限体积法
将连续的求解域离散化为有限个小的 体积，在每个体积上近似函数，将微 分方程转化为代方程进行求解。
有限差分法
向前差分法
将微分方程中的导数项用前一步的函数值近似代替，得到向前差 分方程。
向后差分法
将微分方程中的导数项用后一步的函数值近似代替，得到向后差 分方程。
第五章 对流-扩散方程的离散格式
见下页表格：
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、离散格式的迁移性是指对离散方程中的对流项的某种离散格式，仅能使其扰动顺着流动方向传递的特性，又称为格式的迎风性或传输性。
对流扩散方程背景
对流扩散方程背景提出一种隐格式用于解决二维时间依赖的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 数。
03 第三章 对流与扩散
( ) ( u ) ( v ) ( w ) ( ) ( ) ( )S t x y z x x y y z z
• 离散方程为：
aPP aEE aW W a N N aSS aT T aBB b
六、通用表达式
• 不同化方程
• 一、二维问题的离散化方程 • 二、三维问题的离散化方程
一、二维问题的离散化方程
一、二维问题的离散化方程
• 二维问题的微分方程
( ) ( u ) ( v ) ( ) ( )S t x y x x y y
二、逆风格式
• 逆风格式假定交界面上的 值等于该界面逆风 边节点的 值，故：
w W ( Fw 0)
e P ( Fe 0) e E ( Fe 0)
w P ( Fw 0)
• 所以
Fww W [ Fw ,0] P [ FW ,0] Fee P [ Fe ,0] E [ FE ,0]
Fe aE exp( Fe / De ) 1
aW
Fw exp( Fw / Dw ) exp( Fw / Dw ) 1
aP aE aW ( Fe Fw )
Fe Fw
aP aE aW
• 此时 aP aE aW 都大于0，满足离散化的条件。
四、混合格式
• 混合格式由Spalding提出。
0.1 Fw 5 aW Dw [0, (1 ) ] [0, Fw ] Dw
aP aE aW ( Fe Fw ) aE aW
六、通用表达式
• 引入函数 A( P ) • 方程离散为 aPP aEE aWW • 其中
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
§5.1 对流项离散格式的重要性 及两种离散方式 一、对流项离散格式的重要性
1、数值解的准确性（假扩散） 2、数值解的稳定性 3、数值解的经济性
二、构造离散格式的两种方式
1、Taylor展开法 2、控制容积积分法
1/59
两种定义截差阶数一致，但截差首项系数有所不同。
传热与流体流动的数值计算
2/59
传热与流体流动的数值计算
一、一维对流-扩散问题模型方程的精确解（续）
0 eux 1 e Pex L 1 uL Pe L 0 e 1 e 1
Peclet数：
Pe
uL
0
Pe表示对流与扩散作用 的相对大小。
0
3/59
传热与流体流动的数值计算
三、对流项的迎风格式（续）
e界面
ue 0 , P ; ue 0 , E
u e Fee P max Fe ,0 E max Fe ,0
w界面
P Fe ,0 E Fe ,0
uw 0 , W ; uw 0 , P
aP aE aW Fe Fw Fe Fw , aE De , aW Dw 2 2
在数值计算过程中，如果连续性方程始终得到满足， 则： aP aE aW 在求解过程中，始终保持连续性方程满足非常重要。 常物性条件下均分网格： 1 0.5P E 1 0.5P W P
aE De Fe ,0 , aW Dw Fw ,0
对流项一阶迎风：
aW i 1 aE i 1 P ,0 1 P ,0 P D D
11/59
传热与流体流动的数值计算
2
5/59
传热与流体流动的数值计算
二、对流项的中心差分（续）
例：在一维模型方程离散求解的 均分网格中，已知W =100， E =200。试对P =0，1，2及4 四种情况按中心差分格式计算 P之值。
负系数会导致物理上不真实的解。
6/59
传热与流体流动的数值计算
三、对流项的迎风格式
Taylor展开法
§5.2 对流项的中心差分与迎风格式
一、一维对流-扩散问题模型方程的精确解
d d d u dx dx dx
边界条件： x 0 , 0 ; x L , L
d
udx
C
C2
ln C1
ux
C1eux C2
二、混合格式(Spalding,1971)
0 aE 1 0.5Pe De Pe , Pe 2 , , 2 Pe 2 Pe 2
aE Pe , 1 0.5Pe , 0 De10/59来自传热与流体流动的数值计算
§5.3 对流-扩散方程的混合格式及乘方格式
一、系数aE与aW 之间的内在联系
aE(i)与aW (i+1)共享同一个界面。 对流项中心差分：
Fe Fw aE De , aW Dw 2 2 aW i 1 aE i P P 1 1 P D D 2 2
u w Fww W max Fw ,0 P max Fw ,0
W Fw ,0 P Fw ,0
8/59
传热与流体流动的数值计算
三、对流项的迎风格式（续）
迎风格式离散形式：
aPP aEE aWW
二、对流项的中心差分
d d d u 采用控制容积积分法 对方程 dx dx dx e u e w u w P 2 2 x w x e
e w u e u w E W 2 2 x e x w
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
7/59
传热与流体流动的数值计算
记：F = u 通过界面的流量。 D= x 界面上单位面积扩散阻力的倒数(扩导)。 F u u x = P D x 4/59
传热与流体流动的数值计算
二、对流项的中心差分（续）
Fe Fw Fe Fw De Dw P De E Dw W 2 2 2 2 aPP aEE aWW
a E De Fe , 0 aW Dw Fw , 0 a P a E aW Fe Fw
9/59
传热与流体流动的数值计算
四、中心差分与一阶迎风格式的讨论
1、对流项中心差分在不发生振荡的参数范围内，比一 阶迎风格式的误差更小。 2、一阶迎风格式离散方程系数永远大于零，不会引起 解的振荡，得到物理上看似合理的解。 3、一阶迎风格式截差阶数低，除非采用相当密的网格， 否则计算结果的误差较大。 4、一阶迎风格式的启示：应当在迎风方向取更多的信 息构造格式，更好地反映对流过程的物理本质。 5、在调试程序或计算的中间过程仍可以采用一阶迎风 格式。