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

物理模型：

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

假设棒的截面内温度均匀，问题简化为一维导热模型，其控制方程为
ddx(kAdTdx)−hP(T−T∞)=0(1)\frac{d}{dx} \left( kA\frac{dT}{dx} \right) - hP(T-T_\infin)=0 \tag{1} dxd(kAdxdT)−hP(T−T∞)=0(1)
其中 hhh 是对流换热系数，PPP 是侧面周长， kkk 是棒的导热系数。方程左边第二项在这里相当于源项。

该方程的解析解为，
T−T∞TB−T∞=cosh[n(L−x)]cosh(nL)(2)\frac{T-T_\infin}{T_B-T_\infin}=\frac{cosh[n(L-x)]}{cosh(nL)} \tag{2} TB−T∞T−T∞=cosh(nL)cosh[n(L−x)](2)
其中n2=hP/(kA)n^2=hP/(kA)n2=hP/(kA)， AAA是截面积， LLL是棒的长度。L=1m，n2=25/m2L=1m，n^2=25 /m^2L=1m，n2=25/m2, kAkAkA是常数。

方程离散

网格单元数为5，因为kAkAkA是常数，所以方程式(1)(1)(1)可以简化为，
ddx(dTdx)−n2(T−T∞)=0(3)\frac{d}{dx} \left( \frac{dT}{dx} \right) - n^2(T-T_\infin)=0 \tag{3} dxd(dxdT)−n2(T−T∞)=0(3)
其中，n2=hP/(kA)n^2=hP/(kA)n2=hP/(kA)。

按一维扩散方程的离散推导中方程式(4)(4)(4)的形式，这里相当于
Γ=1S=n2(T−T∞)}(4)\left. \begin{aligned} \Gamma =&1 \\ \\ S =& n^2(T-T_\infin) \end{aligned} \right \} \tag{4}Γ=S=1n2(T−T∞)⎭⎪⎬⎪⎫(4)

则方程离散为，
[(AdTdx)e−(AdTdx)w]−[n2(TP−T∞)Aδx]=0(5)\left[ \left( A\frac{dT}{dx} \right)_e - \left( A\frac{dT}{dx} \right)_w \right] - [n^2 (T_P-T_\infin)A\delta x] = 0 \tag{5}[(AdxdT)e−(AdxdT)w]−[n2(TP−T∞)Aδx]=0(5)

其中，Aδx=ΔVA \delta x =\Delta VAδx=ΔV。由于n2n^2n2是常数，所以n2(T−T∞)=Sˉn^2(T-T_\infin)=\bar Sn2(T−T∞)=Sˉ。

方程两边消去AAA，有
[(TE−TPδx)−(TP−TWδx)]−[n2(TP−T∞)δx]=0(6)\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_\infin)\delta x] = 0 \tag{6} [(δxTE−TP)−(δxTP−TW)]−[n2(TP−T∞)δx]=0(6)
整理之，
(1δx+1δx+n2δx)TP=1δxTW+1δxTE+n2δxT∞(7)\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_\infin \tag{7} (δx1+δx1+n2δx)TP=δx1TW+δx1TE+n2δxT∞(7)
即，
aPTP=aWTW+aETE+Su(8)a_P T_P = a_W T_W + a_E T_E + S_u \tag{8} aPTP=aWTW+aETE+Su(8)
aW=1δxaE=1δxSP=−n2δxSu=n2δxT∞aP=aW+aE−SP}(9)\left. \begin{aligned} 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_\infin \\ \\ a_P = &a_W + a_E - S_P \end{aligned} \right \} \tag{9} aW=aE=SP=Su=aP=δx1δx1−n2δxn2δxT∞aW+aE−SP⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(9)

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

节点1

节点1包含了恒温边界，这是第一类边界条件，与一维扩散方程数值求解（第一类边界条件）中的边界条件相同，方程离散与其中的方程式22、23、2422、23、2422、23、24类似，
[(TE−Tpδx)−(TP−TBδx/2)]−[n2(TP−T∞)δx]=0(10)\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_\infin)\delta x]=0 \tag{10} [(δxTE−Tp)−(δx/2TP−TB)]−[n2(TP−T∞)δx]=0(10)
即，
aW=0aE=1δxSP=−n2δx−2δxSu=n2δxT∞+2δxTBaP=aW+aE−SP}(11)\left. \begin{aligned} 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_\infin + \frac{2}{\delta x}T_B \\ \\ a_P = &a_W + a_E - S_P \end{aligned} \right \} \tag{11} aW=aE=SP=Su=aP=0δx1−n2δx−δx2n2δxT∞+δx2TBaW+aE−SP⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(11)

节点 5

自由端为绝热边界，所以截面的热流密度为零
q=kdTdx=0(12)q=k \frac{dT}{dx}=0 \tag{12} q=kdxdT=0(12)

即，
(dTdx)e=0(13)\left( \frac{dT}{dx} \right)_e = 0 \tag{13} (dxdT)e=0(13)
那么方程离散为
[0−(TP−TWδx)]−[n2(TP−T∞)δx]=0(14)\left[0-\left( \frac{T_P -T_W}{\delta x} \right) \right] - [n^2(T_P-T_\infin)\delta x]=0 \tag{14} [0−(δxTP−TW)]−[n2(TP−T∞)δx]=0(14)
即，
aW=1δxaE=0SP=−n2δxSu=n2δxT∞aP=aW+aE−SP}(15)\left. \begin{aligned} a_W =& \frac{1}{\delta x} \\ \\ a_E = & 0 \\ \\ S_P = & -n^2 \delta x \\ \\ S_u = &n^2 \delta x T_\infin \\ \\ a_P = &a_W + a_E - S_P \end{aligned} \right \} \tag{15} aW=aE=SP=Su=aP=δx10−n2δxn2δxT∞aW+aE−SP⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(15)

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

网格数：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/dxA[i][i+1] = -1/dxA[i][i] = -(A[i][i-1] + A[i][i+1]) - spb[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] = suprint('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_tempxx = 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.