对流扩散方程的离散格式

流体仿真与应用第八讲二、对流-扩散问题的有限体积法◆中心差分格式（例子）节点增加到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 −=◆迎风格式的特点迎风格式满足守恒性。

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是扩散系数，表示物质扩散的速率与浓度梯度的乘积。

对流扩散方程的物理意义是描述了物质在流体中传输的速率与物质浓度梯度之间的关系。

通过对流项，方程能够描述物质随着流体的运动快速传输的现象；而通过扩散项，方程能够描述物质由浓度高处向浓度低处传输的现象。

综合考虑对流和扩散的作用，对流扩散方程能够比较准确地描述物质在流体中的传输过程。

对流扩散方程在科学和工程领域有广泛的应用。

例如，在污染物传输和扩散模拟中，对流扩散方程可用于描述污染物由源区到周围空气或水体的传输过程。

在热传导模拟中，对流扩散方程可用于描述热量由高温区域到低温区域的传导过程。

在物质传递过程中，对流扩散方程也被广泛应用于描绘物质的传输行为。

总结起来，对流扩散方程是一种常见的输运方程，它能够描述物质由流体传输并扩散的过程。

通过对流项和扩散项的综合作用，对流扩散方程能够比较准确地描述物质在流体中的传输行为，所以在科学和工程领域有着广泛的应用。

数值分析-对流项的离散对流项的数值离散扩散项的离散在先前已经介绍，对于动量⽅程，更为棘⼿的问题在于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混合格式乘⽅格式指数格式假扩散。