有限体积法（4）——一维扩散方程数值求解（第二类边界条件）

物理模型：

有对流换热的细棒导热（如下图），固定端为恒温边界，自由端为绝热边界，棒的侧壁面与环境之间存在对流换热。固定端温度$T_B=100℃$，环境温度$T_{\infty }=20℃$。

假设棒的截面内温度均匀，问题简化为一维导热模型，其控制方程为
$$\frac{d}{dx}\left(kA\frac{dT}{dx}\right)-hP(T-T_{\infty })=0\text{\hspace{1em}(1)}$$
其中 $h$ 是对流换热系数，$P$ 是侧面周长， $k$ 是棒的导热系数。方程左边第二项在这里相当于源项。

该方程的解析解为，
$$\frac{T-T_{\infty }}{T_B-T_{\infty }}=\frac{cosh[n(L-x)]}{cosh(nL)}\text{\hspace{1em}(2)}$$
其中$n^2=hP/(kA)$， $A$是截面积， $L$是棒的长度。$L=1m，n^2=25/m^2$, $kA$是常数。

方程离散

网格单元数为5，因为$kA$是常数，所以方程式$(1)$可以简化为，
$$\frac{d}{dx}\left(\frac{dT}{dx}\right)-n^2(T-T_{\infty })=0\text{\hspace{1em}(3)}$$
其中，$n^2=hP/(kA)$。

按一维扩散方程的离散推导中方程式$(4)$的形式，这里相当于
$$\begin{array}{rl}\Gamma =& 1\\ & \\ S=& n^2(T-T_{\infty })\end{array}\}\text{\hspace{1em}(4)}$$

则方程离散为，
$$\left[{\left(A\frac{dT}{dx}\right)}_e-{\left(A\frac{dT}{dx}\right)}_w\right]-[n^2(T_P-T_{\infty })A\delta x]=0\text{\hspace{1em}(5)}$$

其中，$A\delta x=\Delta V$。由于$n^2$是常数，所以$n^2(T-T_{\infty })=\bar{S}$。

方程两边消去$A$，有
$$\left[\left(\frac{T_E-T_P}{\delta x}\right)-\left(\frac{T_P-T_W}{\delta x}\right)\right]-[n^2(T_P-T_{\infty })\delta x]=0\text{\hspace{1em}(6)}$$
整理之，
$$\left(\frac{1}{\delta x}+\frac{1}{\delta x}+n^2\delta x\right)T_P=\frac{1}{\delta x}{T}_W+\frac{1}{\delta x}{T}_E+n^2\delta x{T}_{\infty }\text{\hspace{1em}(7)}$$
即，
$$a_P{T}_P=a_W{T}_W+a_E{T}_E+S_u\text{\hspace{1em}(8)}$$
$$\begin{array}{rl}{a}_W=& \frac{1}{\delta x}\\ & \\ a_E=& \frac{1}{\delta x}\\ & \\ S_P=& -n^2\delta x\\ & \\ S_u=& n^2\delta x{T}_{\infty }\\ & \\ a_P=& a_W+a_E-S_P\end{array}\}\text{\hspace{1em}(9)}$$

式$(8)$和$(9)$是针对于内部节点的，如网格节点2,3,4

节点1

节点1包含了恒温边界，这是第一类边界条件，与一维扩散方程数值求解（第一类边界条件）中的边界条件相同，方程离散与其中的方程式$22、23、24$类似，
$$\left[\left(\frac{T_E-T_p}{\delta x}\right)-\left(\frac{T_P-T_B}{\delta x/2}\right)\right]-[n^2(T_P-T_{\infty })\delta x]=0\text{\hspace{1em}(10)}$$
即，
$$\begin{array}{rl}{a}_W=& 0\\ & \\ a_E=& \frac{1}{\delta x}\\ & \\ S_P=& -n^2\delta x-\frac{2}{\delta x}\\ & \\ S_u=& n^2\delta x{T}_{\infty }+\frac{2}{\delta x}{T}_B\\ & \\ a_P=& a_W+a_E-S_P\end{array}\}\text{\hspace{1em}(11)}$$

节点 5

自由端为绝热边界，所以截面的热流密度为零
$$q=k\frac{dT}{dx}=0\text{\hspace{1em}(12)}$$

即，
$${\left(\frac{dT}{dx}\right)}_e=0\text{\hspace{1em}(13)}$$
那么方程离散为
$$\left[0-\left(\frac{T_P-T_W}{\delta x}\right)\right]-[n^2(T_P-T_{\infty })\delta x]=0\text{\hspace{1em}(14)}$$
即，
$$\begin{array}{rl}{a}_W=& \frac{1}{\delta x}\\ & \\ a_E=& 0\\ & \\ S_P=& -n^2\delta x\\ & \\ S_u=& n^2\delta x{T}_{\infty }\\ & \\ a_P=& a_W+a_E-S_P\end{array}\}\text{\hspace{1em}(15)}$$

把全部节点的离散方程换成$-a_W+a_P-a_E=S_u$ 的形式，然后组成方程组，带入数值，计算求解，数值结果与解析解的对比如下，

网格数：5

网格数：10

网格数：20

计算程序

import numpy as np 
import matplotlib.pyplot as plt
import math

#== parameters ===
nx = 20  # 网格单元数
Nnodes = nx + 1 # 节点数，含左边界点
L = 1 # 长度，m
n2 = 25 # /m2, hP/kA
Tb = 100 # 边界B的温度值 ,℃
q = 0 # 自由端截面的热流密度值, W/m2
Tc = 20 # 环境温度，℃

#==  x grid ==
dx = L / nx  # 网格间距
print('dx = ',dx)
x = np.zeros(Nnodes)  # x网格
x[1:] = np.linspace(dx/2, L-dx/2, nx) # 以边界A为原点创建网格点的坐标值
print('x grid = ', x, '\n')

#==  solution array ==
tt = np.zeros(Nnodes)  # 解向量
tt[0] = Tb # 边界值

#== matrix ==
A = np.zeros((nx, nx))
b = np.zeros(nx)

su = n2 * dx * Tc
sp = -n2 * dx
for i in range(1, nx-1): #内部节点
    A[i][i-1] = -1/dx
    A[i][i+1] = -1/dx
    A[i][i] = -(A[i][i-1] + A[i][i+1]) - sp
    b[i] = su

# for B boundary
i = 0  
A[i][i+1] = -1/dx
su = 2*Tb/dx + n2*dx*Tc
sp = -2/dx - n2*dx
A[i][i] = -A[i][i+1] - sp
b[i] = su

# for insulated boundary
i = nx-1  
A[i][i-1] = -1/dx
su = n2 * dx * Tc
sp = -n2 * dx
A[i][i] = -A[i][i-1] - sp
b[i] = su

print('A = \n', A, '\n')
print('b = \n', np.matrix(b).T ,'\n')

t_temp = np.linalg.solve(A, b)
print('solution = \n', np.matrix(t_temp).T, '\n')
tt[1:] = t_temp

xx = np.linspace(0, L, 50, endpoint=True)
exact_tt = np.zeros(50)
n = math.sqrt(n2)
for i in range(50):
    exact_tt[i] = Tc + (Tb-Tc) * (math.cosh(n * (L-xx[i]))) / math.cosh(n*L)

plt.xlabel('DisTbnce (m)')
plt.ylabel('Temperature (C)')
plt.plot(x,tt ,'bs', label='Numerical')
plt.plot(xx,exact_tt,'k', label='Exct')
plt.legend()
plt.show()

参考资料

Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.