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

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

文章图片

假设棒的截面内温度均匀，问题简化为一维导热模型，其控制方程为
d d x ( k A d T d x ) ? h P ( 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)
其中h h h 是对流换热系数， P P P 是侧面周长，k k k 是棒的导热系数。方程左边第二项在这里相当于源项。
【有限体积法（4）——一维扩散方程数值求解（第二类边界条件）】该方程的解析解为，
T ? T ∞ T B ? T ∞ = c o s h [ n ( L ? x ) ] c o s h ( n L ) (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)
其中 n 2 = h P / ( k A ) n^2=hP/(kA) n2=hP/(kA)，A A A是截面积，L L L是棒的长度。 L = 1 m ， n 2 = 25 / m 2 L=1m，n^2=25 /m^2 L=1m，n2=25/m2,k A kA kA是常数。
方程离散
文章图片

网格单元数为5，因为 k A kA kA是常数，所以方程式 ( 1 ) (1) (1)可以简化为，
d d x ( d T d x ) ? n 2 ( 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)
其中， n 2 = h P / ( k A ) n^2=hP/(kA) n2=hP/(kA)。
按一维扩散方程的离散推导中方程式 ( 4 ) (4) (4)的形式，这里相当于
Γ = 1 S = n 2 ( T ? T ∞ ) } (4) \left. \begin{aligned} \Gamma =&1 \\ \\ S =& n^2(T-T_\infin) \end{aligned} \right \} \tag{4} Γ=S=?1n2(T?T∞?)???????(4)
则方程离散为，
[ ( A d T d x ) e ? ( A d T d x ) w ] ? [ n 2 ( T P ? 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 = Δ V A \delta x =\Delta V Aδx=ΔV。由于 n 2 n^2 n2是常数，所以 n 2 ( T ? T ∞ ) = S ˉ n^2(T-T_\infin)=\bar S n2(T?T∞?)=Sˉ。
方程两边消去 A A A，有
[ ( T E ? T P δ x ) ? ( T P ? T W δ x ) ] ? [ n 2 ( T P ? 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 + n 2 δ x ) T P = 1 δ x T W + 1 δ x T E + n 2 δ x T ∞ (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?=δx1?TW?+δx1?TE?+n2δxT∞?(7)
即，
a P T P = a W T W + a E T E + S u (8) a_P T_P = a_W T_W + a_E T_E + S_u \tag{8} aP?TP?=aW?TW?+aE?TE?+Su?(8)
a W = 1 δ x a E = 1 δ x S P = ? n 2 δ x S u = n 2 δ x T ∞ a P = a W + a E ? S P } (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 、 24 22、23、24 22、23、24类似，
[ ( T E ? T p δ x ) ? ( T P ? T B δ x / 2 ) ] ? [ n 2 ( T P ? 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)
即，
a W = 0 a E = 1 δ x S P = ? n 2 δ x ? 2 δ x S u = n 2 δ x T ∞ + 2 δ x T B a P = a W + a E ? S P } (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?δx2?n2δxT∞?+δx2?TB?aW?+aE??SP??????????????????????????????????????????????(11)
节点 5
自由端为绝热边界，所以截面的热流密度为零
q = k d T d x = 0 (12) q=k \frac{dT}{dx}=0 \tag{12} q=kdxdT?=0(12)
即，
( d T d x ) e = 0 (13) \left( \frac{dT}{dx} \right)_e = 0 \tag{13} (dxdT?)e?=0(13)
那么方程离散为
[ 0 ? ( T P ? T W δ x ) ] ? [ n 2 ( T P ? 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)
即，
a W = 1 δ x a E = 0 S P = ? n 2 δ x S u = n 2 δ x T ∞ a P = a W + a E ? S P } (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?=?δx1?0?n2δxn2δxT∞?aW?+aE??SP??????????????????????????????????????????(15)
把全部节点的离散方程换成 ? a W + a P ? a E = S u -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/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] = 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.