1 数值方法A 续
下一步是用中心差分法近似\(\frac{dT}{dx}\)。
首先，对温度函数\(T(x)\)用泰勒展开，其中由于2阶以上项计算复杂、对结果影响小，故忽略。
假设在节点之间温度线性变化。
\[T(x)=T(x_i)+(x-x_i)\frac{dT}{dx}|_{x_i}+\frac{(x-x_i)^2}2\frac{d^2T}{dT^2}|_{x_i} \]\[T(x_{P})=T(x_{e})-(\delta x)_{e^{-}}\cdot(-\frac{dT}{dx})|_{x_{e}}+\frac{(\delta x)_{e^{-}}^2}{2}\cdot(\frac{d^2T}{dT^2})|_{x_{e}} \]\[T(x_E)=T(x_e)-(\delta x)_{e^+}\cdot (\frac{dT}{dx})|_{x_e}+\frac{(\delta x)_{e^+}^2}2 \cdot(\frac{d^2T}{dT^2})|_{x_e} \]令e为P、E中点，则有\((\delta x)_{e^-}=(\delta x)_{e^+}=\frac{1}{2} (\delta x)_e\)，故：
\[\frac{dT}{dx}|_{x_e}=\frac{T(x_E)-T(x_P)}{(\delta x)_e} \]回代入公式\((\lambda \frac{dT}{dx})_e-(\lambda \frac{dT}{dx})_w+\bar{S} \Delta x=0\)得：
\[(\lambda \frac{dT}{dx})_e=\lambda_e \cdot \frac{T_E-T_P}{(\delta x)_e} \]\[(\lambda \frac{dT}{dx})_w=\lambda_w \cdot \frac{T_P-T_W}{(\delta x)_w} \]另有：
\[\bar{S}=S_p \cdot T_P+S_c \]把源项在控制体积内的均值\(\bar{S}\)线性化，其中\(S_P\)为线性系数（或斜率），\(S_c\)为常数项（或截距），二者均应当为常数。
那么，将我们得到的这三个公式组合起来，可以得到最终的线性化方程：
\[a_w T_w +a_p T_p+a_ET_E=b \]其中，代表三个节点（当前控制节点左边一个、它本身、右边一个）的系数，以及右边一项b，分别是：
\[\quad a_W=-\frac{\lambda_w}{\delta x_w}\quad a_P=\frac{\lambda_w}{\delta x_w}+\frac{\lambda_e}{\delta x_e}-S_P\Delta x\quad a_E=-\frac{\lambda_e}{\delta x_e}\quad b=S_C\Delta x \]注意，这个式子只能用于内部控制体积，不适用于边界节点。
对于边界条件：
-
已知温度：（\(T(x=0)=T_a\text,~T(x=L)=T_b\)）
\[a_PT_P+a_ET_E=b,~a_P=1,~a_E=0,~b=T_a\\a_WT_W+a_PT_P=b,~a_P=1,~a_W=0,~b=T_b \] -
已知温度随空间变化的导数，描述热流边界条件：
\[(\lambda\frac{dT}{dx})_i-\left(\lambda\frac{dT}{dx}\right)_B+\bar{S}\Delta x_B=0 \]
如图，定义热通量\(q_B=-\lambda \frac{dT}{dx}\)。这里可以把\(-\left(\lambda\frac{dT}{dx}\right)_B\)这一项用\(q_B\)代替，于是可推出最后的线性化方程：
\[a_PT_P+a_ET_E=b,~a_P=\frac{\lambda_e}{\delta x_e}-S_p\Delta x_B,~a_E=-\frac{\lambda_e}{\delta x_e},~\mathrm{b}=S_p\Delta x_B+q_B \]
如图，最后得到了一个三对角矩阵。每个节点只影响其左侧一个和右侧一个节点。
求解线性方程组
判断解的稳定性
-
解是否有边界？
-
边界条件是否被遵守？
-
结果如何与理论解或实验数据比较，以验证数值解的正确性？
-
残差误差：
\[R=\frac{|B-A\cdot T|}{|diag(A)\cdot T|}<tol \]其中，\(|B-A\cdot T|\)是残差，表示实际计算的结果\(A\cdot T\)与理论上应当得到的结果B之间的差；
\(diag(A)\cdot T\)是矩阵A的对角线部分与T的乘积，用于对残差进行非量纲化，使得残差大小与系统的尺度保持一致。
若R小于某个预设的容差tol，可认为解的精读是能接受的。 -
能量平衡：能量守恒必须在每个控制体中保持。有时源项\(\bar{S}=0\)，则控制体内的流入和流出热量\(q_w - q_e = 0\)。
-
网格收敛性：若网格细化，节点间距减半，误差应减少4倍。这种收敛特性在已知精确解时表现最好，通常通过对数-对数图（log-log plot）来评估收敛行为。
直接求解：高斯消元法
对如下矩阵：
\[\begin{bmatrix}a_{1,1}&a_{1,2}&a_{1,3}&...&a_{1,n}\\a_{2,1}&a_{2,2}&a_{2,3}&...&a_{2,n}\\\vdots&\vdots&\vdots&\ddots&\vdots\\a_{n,1}&a_{n,2}&a_{n,3}&...&a_{n,n}\end{bmatrix}\times\begin{bmatrix}T_1\\T_2\\\vdots\\T_n\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\\vdots\\b_n\end{bmatrix} \]-
前向消元：将第i行减去第1行的某个倍数，使得\(a_{i,~1}=0\)，这样就能消去除了\(a_{1,~1}\)以外，矩阵A第1列的所有数；
重复直到矩阵A变成上三角矩阵。
-
回代：从最后一行开始，最后一行相当于解:
\[a_{n,~n}^{(N-1)}\cdot T_n = b_n^{(N-1)} \]倒数第2行相当于解：
\[a_{n-1,n-1}^{(N-1)}T_{n-1}+a_{n-1,n}^{(N-1)}T_n=b_{n-1}^{(N-1)} \]最后得到的通项公式为：
\[T_i=\frac{b_i^{(N-1)}-\sum_{k=i+1}^na_{i,k}^{(N-1)}T_k}{a_{i,i}^{(N-1)}} \]
迭代求解：高斯-赛德尔法
通项公式：
\[T_i^{(m)}=\frac{b_i}{a_{i,i}}-\sum_{j=1}^{i-1}\frac{a_{i,j}}{a_{i,i}}T_j^{(m)}-\sum_{j=i+1}^n\frac{a_{i,j}}{a_{i,i}}T_j^{(m-1)}\\ =\frac{1}{a_{i,~i}}(b_i-\sum_{j=1}^{i-1}a_{i,j}T_j^{(m)}-\sum_{j=i+1}^na_{i,j}T_j^{(m-1)}) \]其中：
-
\(T_i^{(m)}\)表示第m次迭代中，第i个未知数的近似值；
-
\(T_j^{(m)}\)表示第m次迭代中，第j个未知数近似值，j<i时，这些值在本次迭代中已更新；
-
\(T_j^{(m-1)}\)表示第m次迭代中，第j个未知数近似值，j>i时，这些值在本次迭代中未更新；
-
\(a_{i,~i}\)：矩阵A第i行i列元素；
-
\(b_i\)：矩阵b第i行元素；
-
\[\sum_{j=1}^{i-1}a_{i,j}T_j^{(m)}
\]
这一项对在第i行中，i之前的未知数项求和（已迭代）
-
\[\sum_{j=i+1}^na_{i,j}T_j^{(m-1)}
\]
这一项对在第i行中，i之后的未知数项求和（未迭代）。
该法的物理意义为：通过\(b_i\)减掉加权求和，意为扣除已知部分，剩余未知数项的影响；除以\(a_{i,~i}\)则解出未知数。
标签：frac,cdot,数值,dx,delta,dT,方法,lambda From： https://www.cnblogs.com/yukibvvd/p/18513261/1-numerical-method-a-continuation-n0s1t