对流项的数值离散
扩散项的离散在先前已经介绍，对于动量方程，更为棘手的问题在于N-S方程中对流项和压力项的离散处理。
在不可压流场的数值求解过程中，最重要的问题均是这两个一阶偏导项所引起的。
这里主要讨论关于对流-扩散项的离散方法。
对流项的离散方法
首先是关于对流项的两种离散思路的介绍
-
Taylor展开离散
Taylor展开的方法，就是直接套用泰勒公式将偏导转化为线性离散
-
控制容积积分
通过将对流项的一阶导数在控制容积P中进行积分，利用相邻两个结点的值来获得控制容积边界e，w处的插值
两种思路所给出的截断误差的阶数是相同的
一维稳态对流-扩散的离散
由于整个N-S方程还是过于复杂，这里我们对整个物理问题可以进行相应的简化，从仅含有对流项和扩散的场景出发，来讨论如何解决对流项中的一阶导离散问题。
因此，首先对于一维稳态无内热源的对流扩散问题的数学描述，我们可以很轻松的给出方程的守恒形式：
\[\frac{d}{dx}(\rho u \phi) = \frac{d}{dx}(\Gamma \frac{d \phi}{dx}) \tag{1} \]边界条件写为：
\[\phi(0) = \phi_0； \phi(L) = \phi_L \tag{1.1} \]值得注意的是，以下内容表示的是对对流项（convection term）采用不同的差分格式进行描述，而公式中的扩散项仍然保持中心差分格式。
中心差分格式
采用控制容积法，对流项的中心差分格式相当于在控制容积P的界面上取分段的线性函数，即
\[(\rho u \phi)_e - (\rho u \phi)_w = (\rho u )_e \frac{\phi_E +\phi_P}{2}-(\rho u )_w \frac{\phi_P + \phi_W}{2} \tag{2.1} \]最终的离散结果可以写为：
\[\phi_P[ \frac{(\rho u)_e}{2} - \frac{(\rho u)_w}{2} + \frac{\Gamma_e}{(\delta x)e} + \frac{\Gamma_w}{(\delta x)w}] = \phi_E(\frac{\Gamma_e}{(\delta x)_e} - \frac{(\rho u)_e}{2}) + \phi_W(\frac{\Gamma_w}{(\delta x)_w} +\frac{(\rho u)_w}{2}) \tag{2} \]这里可以将界面上的流量记为$ F = \rho u\(，界面上单位扩散阻力的倒数记为\)D = \frac{\Gamma}{\delta x}$
最终式\((2)\)可以化简为：
\[a_P \phi_P = a_E \phi_E + a_W \phi_W \tag{2.2} \]python代码可以轻松实现中心差分算法的描述：
# 中心差分格式（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 x
y_nodes = number of nodes along y
D_x = Diffusivity along x, type = float
D_y = Diffusivity along y, type = float
F_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 calculations
w_bound = [0] * y_nodes # creating 0-matricies to find
e_bound = [0] * y_nodes # all boundary nodes
s_bound = [0] * x_nodes
n_bound = [0] * x_nodes
if debug == 'y': # Debugger notifies progress
print("finding boundaries")
for i in range (x_nodes): # boundary nodes are populated
s_bound[i] = i # used to populate matricies
n_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 -1
if debug == 'y': # Debugger prints boundary nodes
xB = "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 variables
a_w = [0] * dim # Creating empty matricies for all nodes
a_e = [0] * dim
a_s = [0] * dim
a_n = [0] * dim
sp_p = 0
su_p = 0
B = [0] * dim
sp_W = [0] * dim
su_W = [0] * dim
sp_E = [0] * dim
su_E = [0] * dim
sp_S = [0] * dim
su_S = [0] * dim
sp_N = [0] * dim
su_N = [0] * dim
#Boundaries
phi_W = Boundary[0]
phi_E = Boundary[1]
phi_S = Boundary[2]
phi_N = Boundary[3]
## answer matricies ##
a = [0] * dim # make a
for i in range(dim): # 0-array
a[i] = [0] * dim
B = [0] * dim
# POPULATING USING CENTRAL DIFFERENCE
print("populating central difference")
pbar = tqdm(range(dim)) # load progress bar
for i in range(dim):
pbar.update(1) # progress indicator
if i not in w_bound: # populating non-boundary conditions
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 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 conditions
sp_W[i] = -(2*D_x + F_x[i])
su_W[i] = (2*D_x + F_x[i])* phi_W
a[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_E
a[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_S
a[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_N
a[i][i] = a[i][i] - sp_N[i]
B[i] = B[i] + su_N[i]
return a, B
迎风差分格式
利用控制容积法，对对流项的一阶迎风差分格式进行描述：
| 在e界面上 | \(u_e > 0, \phi = \phi_P\) | \(u_e<0,\phi = \phi_E\) |
|---|---|---|
| 在w界面上 | \(u_w > 0, \phi = \phi_W\) | \(u_w < 0, \phi = \phi_P\) |
因此对流迎风差分格式可以表达为：
\[(\rho u \phi)_e = F_e \phi_e = \phi_P MAX[F_e,0] - \phi_E MAX[-F_e,0] \tag{3.1} \]最终的离散结果可写为：
\[\phi_P[ MAX[F_e,0] +MAX[-F_w,0] + \frac{\Gamma_e}{(\delta x)e} + \frac{\Gamma_w}{(\delta x)w}] \\ = \phi_E(\frac{\Gamma_e}{(\delta x)_e} + MAX[-F_e,0]) + \phi_W(\frac{\Gamma_w}{(\delta x)_w} +MAX[F_w,0]) \tag{3} \]需要注意的是：
- 在数值解不出现振荡的范围内，中心差分要比迎风差分结果的误差小
- 一阶迎风格式在应用上误差较大，应用需由限制，但新发展出的二阶迎风、三阶迎风都从中吸取了思想
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 x
y_nodes = number of nodes along y
D_x = Diffusivity along x, type = float
D_y = Diffusivity along y, type = float
F_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 calculations
w_bound = [0] * y_nodes # creating 0-matricies to find
e_bound = [0] * y_nodes # all boundary nodes
s_bound = [0] * x_nodes
n_bound = [0] * x_nodes
if debug == 'y': # Debugger notifies progress
print("finding boundaries")
for i in range (x_nodes): # boundary nodes are populated
s_bound[i] = i # used to populate matricies
n_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 -1
if debug == 'y': # Debugger prints boundary nodes
xB = "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_nodes
a_w = [0] * dim # Creating empty matricies for all nodes
a_e = [0] * dim
a_s = [0] * dim
a_n = [0] * dim
sp_p = 0
su_p = 0
B = [0] * dim
sp_W = [0] * dim
su_W = [0] * dim
sp_E = [0] * dim
su_E = [0] * dim
sp_S = [0] * dim
su_S = [0] * dim
sp_N = [0] * dim
su_N = [0] * dim
#Boundaries
phi_W = Boundary[0]
phi_E = Boundary[1]
phi_S = Boundary[2]
phi_N = Boundary[3]
## answer matricies ##
a = [0] * dim # make a
for i in range(dim): # 0-array
a[i] = [0] * dim
B = [0] * dim
# POPULATING USING CENTRAL DIFFERENCE
print("populating upwind difference")
pbar = tqdm(range(dim))
for i in range(dim):
pbar.update(1)
if i not in w_bound: # populating non-boundary nodes
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 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 nodes
sp_W[i] = -(2*D_x + max((F_x[i]),0))
su_W[i] = (2*D_x + max((F_x[i]),0))* phi_W
a[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_E
a[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_S
a[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_N
a[i][i] = a[i][i] - sp_N[i]
B[i] = B[i] + su_N[i]
return a, B