扩散方程的数值解法
一维稳态导热的离散
一维稳态导热问题的数学描述可表示为：
\[\begin{equation}
\frac{1}{A(x)} \frac{d}{dx} [ \lambda A(x) \frac{d T}{dx}] + S = 0
\end{equation}
\tag{1}
\]
其中\(S\)表示为源项，假设源项可表示为温度的线性函数，即\(S = S_C + S_P T_P\)
对方程\((1)\)式两边同乘\(A(x)\) ，并对控制体\(P\) 积分，可得
\[\begin{equation}
T_P(\frac{\lambda_e A_e}{(\delta x)_e} + \frac{\lambda_w A_w}{(\delta x)_w} -S_P A_P\Delta x) = T_W (\frac{\lambda_w A_w}{(\delta x)_w}) +T_E(\frac{\lambda_e A_e}{(\delta x)_e}) + S_C A_P \Delta x
\end{equation}
\tag{2}
\]
最终可将其简化为：
\[a_P T_P = a_W T_W + a_E T_E + b
\tag{3}
\]
\((3)\)式即为一维稳态导热问题的最终离散形式，这时的问题关键在于表示控制体上的导热系数\(\lambda_w , \lambda_e\)，确定界面上的导热系数，方法包括算术平均(arithmetic mean) 和调和平均( harmonic mean)两种，经验证，调和平均法更接近物理事实。
一维非稳态导热的离散
一维非稳态导热问题的数学描述为：
\[\begin{equation}
\rho c \frac{\part T}{\part \tau} = \frac{1}{A(x)} \frac{d}{dx} [ \lambda A(x) \frac{d T}{dx}] + S
\end{equation}
\tag{4}
\]
对控制体在\([\tau , \tau + \Delta \tau]\)其进行积分，
\[(\rho c)_P A_P \Delta x (T_P^{\tau + \Delta \tau} - T_P^\tau) = \int_\tau^{\tau + \Delta \tau} [\frac{\lambda_e A_e (T_E - T_P)}{(\delta x )_e} - \frac{\lambda_w A_w(T_P - T_W)}{(\delta x)_w}]d\tau \\ + \int_\tau^{\tau + \Delta \tau}(S_C + S_PT_P) A_P \Delta x d\tau
\tag{5}
\]
对于积分部分的离散，引入加权因子\(f\)，将温度在时间段\(\Delta\tau\)上的积分表示为：
\(\int_\tau^{\tau + \Delta \tau} T d\tau = [f T - (1-f)T^0]\Delta \tau\)
\((5)\)式的积分式可表示为：
\[(\rho c)_P \frac{A_P \Delta x}{\Delta \tau} (T_P - T_P^0) = f[\frac{\lambda_e A_e (T_E - T_P)}{(\delta x )_e} - \frac{\lambda_w A_w(T_P - T_W)}{(\delta x)_w}] \\ +(1-f)[\frac{\lambda_e A_e (T_E^0 - T_P^0)}{(\delta x )_e} - \frac{\lambda_w A_w(T_P^0 - T_W^0)}{(\delta x)_w}]\\ + [f(S_C + S_PT_P)- (1-f)(S_C + S_PT_P^0)] A_P \Delta x
\tag{6}
\]
进一步简化为：
\[a_P T_P = a_E[f T_E + (1-f)T_E^0] + a_W[fT_W + (1-f)T_W^0] \\+ T_P^0[a_P^0 - (1-f)a_W +(1-f)S_P A_P \Delta x] + S_C A_P \Delta x
\tag{7}
\]
其中，
\[a_P = f a_E + f a_W + a_P^0 - f S_P A_P \Delta x \tag{7.1}
\]
\[a_P^0 = \frac{(\rho c)_P A_P \Delta x}{\Delta x}
\tag{7.2}
\]
\[a_E = \frac{\lambda _e A_e}{(\delta x)e} \tag{7.3}
\]
\[a_W = \frac{\lambda _w A_w}{(\delta x)w} \tag{7.4}
\]
通过改变加权因子\(f\) 的值，0，1，0.5可以分别得到显式(explicit)，隐式(implicit)，Crank-Nicolson格式。
多维导热问题
在直角坐标系下，二维非稳态导热问题的数学描述为：
\[\begin{equation}
\rho c\frac{\part T}{\part t} = \frac{\part}{\part x}(\lambda \frac{\part T}{\part x})+ \frac{\part}{\part y}(\lambda \frac{\part T}{\part y}) +S
\end{equation}
\tag{8}
\]
分别对控制方程的各项进行积分离散，
对于非稳态项，对控制体P积分后离散可得：
\[\int _t ^{t + \Delta t} \int _x ^ n \int_w ^e \rho c \frac{\part T}{\part t}dx dy dt = (\rho c) _P (T _P - T_P^0)\Delta x \Delta y
\tag{8.1}
\]
对于扩散项，对控制体P积分后离散可得：
\[\int _t ^{t + \Delta t} \int _x ^ n \int_w ^e \frac{\part}{\part x}(\lambda \frac{\part T}{ \part x}) dx dy dt + \int _t ^{t + \Delta t} \int_w ^e \int _x ^ n\frac{\part}{\part y}(\lambda \frac{\part T}{ \part y}) dy dx dt
\\ = [\lambda _e \frac{T_E - T_P}{(\delta x)_e}- \lambda _w\frac{T_P - T_W}{(\delta x)_w}]\Delta y \Delta t + [\lambda _n \frac{T_N - T_P}{(\delta y)_n}- \lambda_s\frac{T_P - T_S}{(\delta y)_s}] \Delta x \Delta t
\tag{8.2}
\]
对于源项，对控制体P积分后离散可得：
\[\int _t ^{t + \Delta t} \int _x ^ n \int_w ^e S dx dy dt = (S_C +S_P T_P)\Delta x \Delta y \Delta t
\tag{8.3}
\]
对于上式，可以整理成以下形式：
\[a_P T_P = a_E T_E + a_WT_W + a_N T_N +a_ST_S + b
\tag{9}
\]
其中，
\[a_P = a_W + a_E +a_N + a_S +a_P^0 - S_P \Delta x \Delta y \tag{9.1}
\]
\[a_W = \frac{\Delta y \lambda_w}{(\delta x)_w}, a_E = \frac{\Delta y \lambda_e}{(\delta x)_e}, a_N = \frac{\Delta x \lambda_n}{(\delta y)_n},a_S = \frac{\Delta x \lambda_s}{(\delta x)_s} \tag{9.2}
\]
\[a_P^0 = \frac{(\rho c)_P \Delta x \Delta y}{\Delta t} \tag{9.3}
\]
\[b = S_C \Delta x \Delta y + a_P^0 T_P^0 \tag{9.4}
\]
式\((9) \ -(9.4)\) 即为二维直角系导热问题的离散方程，界面上的导热系数可通过上述的调和平均法( harmonic mean )求得。
源项处理
上式中的源项包括所有不能代表控制方程里非稳态项，对流项，扩散项的项。对于源项的处理十分重要，上述案例中通过将源项局部线性化，即\(S = S_C +S_P T_P\)，相比较与将源项看为常数的方法，且每次$S $ 都能随着\(T_P\) 变化更新值，因此局部线性化更能保证源项更接近物理事实。
此外，若需要保证代数方程的求解不收敛，需要保证\(S_P < 0\)
边界处理
对于第二类和第三类边界条件上的处理方法主要有两种，附加源项( additional source term method )和补充边界结点代数方程的方法，经过实验验证，附加源项法相比于补充结点法，更能明确的表现物理意义，且实施方法简单。因此以下主要介绍附加源项法的实施过程。
附加源项法，简而言之是将第二类和第三类边界条件所规定的热量条件，转化为同控制容积中的源项具有同等地位的附加源项。
假设结点P在x方向上的前一个结点W，正好是外边界
此时，对上述多维非稳态的离散方程\((9)\) 进行变化，消去不可知的温度\(T_W\)，可得：
\[(a_P - a_W)T_P = a_ET_E + a_N T_N+a_S T_S + a_W(T_P - T_W) + b
\tag{10}
\]
由Fourier定律，将\((T_P - T_W)\)转化为热流量\(q_B \Delta y\)
因此，对于第二类边界条件来说，\(q_B\) 是已知量，最终可将离散方程转化为：
\[a^* _P T_P =a_ET_E + a_N T_N+a_S T_S + (S_C + \frac{q_B \Delta y}{\Delta x \Delta y})\Delta x \Delta y
\\=a_ET_E + a_N T_N+a_S T_S + (S_C +S_{C,ad})\Delta x \Delta y
\tag{10.1}
\]
而关于第三类边界条件\(q_B = (T_f - T_W)h\)，需要联立对于\(w\) 界面上的Fourier定律$ q_B = \lambda_B (T_W - T_P)/(\delta x)_w $
也即是说，外边界温度\(T_f\) 和最外层的结点温度\(T_P\) 之间的热阻包括原先假想的结点边界\(w\)处的导热热阻加上外面的对流热阻，$1/h + (\delta x)_w/\lambda_B $
最终对于\((10)\) 式，化简可得：
\[(a^*_P + \frac{A}{1/h + (\delta x)_w/\lambda_B })T_P = (a* + S_{C,ad}) \\ = a_ET_E + a_N T_N+a_S T_S +\{ S_C +\frac{A T_f}{\Delta V [1/h + (\delta x)_w/\lambda_B]} \}\Delta V
\\ =a_ET_E + a_N T_N+a_S T_S + (S_C + S_{P,ad}) \Delta V
\tag{10.2}
\]
对于附加源项法的具体实施步骤：
- 计算边界上的附加源项 \(S_{C,ad}, S_{P,ad}\)，并将其分别带入原有的\(S_C, S_P\) 中
- 令边界上的导热系数\(\lambda _B = 0\) ，使\(a^* _W = 0\)
- 对内结点建立离散方程组，求解代数方程组
- 获得收敛后的解后，按照Fourier 定律或Newton冷却定律解出边界层上的温度\(T_f\)
Reference:
Patankar S. Numerical heat transfer and fluid flow[M]. Taylor & Francis, 2018.
浙公网安备 33010602011771号