本文详细介绍了如何使用热平衡法和高斯-赛德尔迭代法解决二维稳态导热问题。问题背景为一墙角的温度分布，设定了等温和对流两种边界条件，给出了不同类型的节点方程，并通过FORTRAN程序实现了迭代计算。最后，讨论了如何用MATLAB进行数据可视化，展示温度分布情况。
目录
(1) 物理问题
已知一个建筑物墙角如图所示，其中墙体为二维、稳态、无内热源、常物性导热。现在要求用数值方法求出整个墙体的温度分布情况。已知：
（1）在等温边界条件下，墙角外表面的温度 t w , e = 30 ℃ t_{w,e}=30℃ tw,e=30℃，内表面温度 t w , i = 0 ℃ t_{w,i}=0℃ tw,i=0℃。
（2）在对流边界条件下，墙角外表面与周围流体的对流换热系数 h e = 10 W / ( m 2 ⋅ K ) h_{e}=10W/(m^{2}{\cdot}K) he=10W/(m2⋅K)，外表面流体温度 t f , e = 30 ℃ t_{f,e}=30℃ tf,e=30℃；墙角内表面与周围流体的对流换热系数 h i = 4 W / ( m 2 ⋅ K ) h_{i}=4W/(m^{2}{\cdot}K) hi=4W/(m2⋅K)，内表面流体温度 t f , i = 10 ℃ t_{f,i}=10℃ tf,i=10℃。
(2) 用热平衡法建立节点代数方程
用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点（下面将其称为节点）作为需要确定温度值的空间点位置。相邻节点之间的距离称为步长，在两个方向上分别为 Δ x \Delta{x} Δx、 Δ y \Delta{y} Δy（为了研究问题方便，这里取 Δ x = Δ y = 0.1 m \Delta{x}=\Delta{y}=0.1m Δx=Δy=0.1m）。节点的位置用数对 ( m , n ) (m,n) (m,n)表示，每一个节点可以代表以它为中心的一个小区域（由相邻节点连线的中垂线构成）。
(2.1) 内节点
λ t m − 1 , n − t m , n Δ x Δ y + λ t m + 1 , n − t m , n Δ x Δ y + λ t m , n + 1 − t m , n Δ y Δ x + λ t m , n − 1 − t m , n Δ y Δ x = 0 \lambda{\frac{t_{m-1,n}-t_{m,n}}{\Delta{x}}}{\Delta{y}}+\lambda{\frac{t_{m+1,n}-t_{m,n}}{\Delta{x}}}{\Delta{y}}+\lambda{\frac{t_{m,n+1}-t_{m,n}}{\Delta{y}}}{\Delta{x}}+\lambda{\frac{t_{m,n-1}-t_{m,n}}{\Delta{y}}}{\Delta{x}}=0 λΔxtm−1,n−tm,nΔy+λΔxtm+1,n−tm,nΔy+λΔytm,n+1−tm,nΔx+λΔytm,n−1−tm,nΔx=0
(2.2) 外部平直边界节点
等温边界条件：
t m , n = T w = c o n s t t_{m,n}=T_{w}=const tm,n=Tw=const
对流边界条件：
λ t m − 1 , n − t m , n Δ x Δ y + λ t m , n + 1 − t m , n 2 Δ y Δ x + λ t m , n − 1 − t m , n 2 Δ y Δ x + h Δ y ( t f − t m , n ) = 0 \lambda{\frac{t_{m-1,n}-t_{m,n}}{\Delta{x}}}{\Delta{y}}+\lambda{\frac{t_{m,n+1}-t_{m,n}}{2\Delta{y}}}{\Delta{x}}+\lambda{\frac{t_{m,n-1}-t_{m,n}}{2\Delta{y}}}{\Delta{x}}+h\Delta{y}(t_{f}-t_{m,n})=0 λΔxtm−1,n−tm,nΔy+λ2Δytm,n+1−tm,nΔx+λ2Δytm,n−1−tm,nΔx+hΔy(tf−tm,n)=0
绝热边界条件（此问题没有涉及）：
λ t m − 1 , n − t m , n Δ x Δ y + λ t m , n + 1 − t m , n 2 Δ y Δ x + λ t m , n − 1 − t m , n 2 Δ y Δ x = 0 \lambda{\frac{t_{m-1,n}-t_{m,n}}{\Delta{x}}}{\Delta{y}}+\lambda{\frac{t_{m,n+1}-t_{m,n}}{2\Delta{y}}}{\Delta{x}}+\lambda{\frac{t_{m,n-1}-t_{m,n}}{2\Delta{y}}}{\Delta{x}}=0 λΔxtm−1,n−tm,nΔy+λ2Δytm,n+1−tm,nΔx+λ2Δytm,n−1−tm,nΔx=0
(2.3) 外部角点
等温边界条件：
t m , n = T w 1 = T w 2 = c o n s t t_{m,n}=T_{w1}=T_{w2}=const tm,n=Tw1=Tw2=const
对流边界条件：
λ t m − 1 , n − t m , n 2 Δ x Δ y + λ t m , n + 1 − t m , n 2 Δ y Δ x + h ( Δ x 2 + Δ y 2 ) ( t f − t m , n ) = 0 \lambda{\frac{t_{m-1,n}-t_{m,n}}{2\Delta{x}}}{\Delta{y}}+\lambda{\frac{t_{m,n+1}-t_{m,n}}{2\Delta{y}}}{\Delta{x}}+h(\frac{\Delta{x}}{2}+\frac{\Delta{y}}{2})(t_{f}-t_{m,n})=0 λ2Δxtm−1,n−tm,nΔy+λ2Δytm,n+1−tm,nΔx+h(2Δx+2Δy)(tf−tm,n)=0
绝热边界条件（此问题没有涉及）：
λ t m − 1 , n − t m , n Δ x Δ y + λ t m , n + 1 − t m , n 2 Δ y Δ x = 0 \lambda{\frac{t_{m-1,n}-t_{m,n}}{\Delta{x}}}{\Delta{y}}+\lambda{\frac{t_{m,n+1}-t_{m,n}}{2\Delta{y}}}{\Delta{x}}=0 λΔxtm−1,n−tm,nΔy+λ2Δytm,n+1−tm,nΔx=0
(2.4) 内部角点
等温边界条件：
t m , n = T w = c o n s t t_{m,n}=T_{w}=const tm,n=Tw=const
对流边界条件：
λ t m − 1 , n − t m , n Δ x Δ y + λ t m , n + 1 − t m , n Δ y Δ x + λ t m , n − 1 − t m , n 2 Δ y Δ x + λ t m + 1 , n − t m , n 2 Δ x Δ y + h ( Δ x 2 + Δ y 2 ) ( t f − t m , n ) = 0 \lambda{\frac{t_{m-1,n}-t_{m,n}}{\Delta{x}}}{\Delta{y}}+\lambda{\frac{t_{m,n+1}-t_{m,n}}{\Delta{y}}}{\Delta{x}}+\lambda{\frac{t_{m,n-1}-t_{m,n}}{2\Delta{y}}}{\Delta{x}}+\lambda{\frac{t_{m+1,n}-t_{m,n}}{2\Delta{x}}}{\Delta{y}}+h(\frac{\Delta{x}}{2}+\frac{\Delta{y}}{2})(t_{f}-t_{m,n})=0 λΔxtm−1,n−tm,nΔy+λΔytm,n+1−tm,nΔx+λ2Δytm,n−1−tm,nΔx+λ2Δxtm+1,n−tm,nΔy+h(2Δx+2Δy)(tf−tm,n)=0
绝热边界条件（此问题没有涉及）：
λ t m − 1 , n − t m , n Δ x Δ y + λ t m , n + 1 − t m , n Δ y Δ x + λ t m , n − 1 − t m , n 2 Δ y Δ x + λ t m + 1 , n − t m , n 2 Δ x Δ y = 0 \lambda{\frac{t_{m-1,n}-t_{m,n}}{\Delta{x}}}{\Delta{y}}+\lambda{\frac{t_{m,n+1}-t_{m,n}}{\Delta{y}}}{\Delta{x}}+\lambda{\frac{t_{m,n-1}-t_{m,n}}{2\Delta{y}}}{\Delta{x}}+\lambda{\frac{t_{m+1,n}-t_{m,n}}{2\Delta{x}}}{\Delta{y}}=0 λΔxtm−1,n−tm,nΔy+λΔytm,n+1−tm,nΔx+λ2Δytm,n−1−tm,nΔx+λ2Δxtm+1,n−tm,nΔy=0
(3) 高斯-赛德尔(Gauss-Seidel)迭代法
现有待求解三元一次方程组：
{ a 11 t 1 + a 12 t 2 + a 13 t 3 = b 1 a 21 t 1 + a 22 t 2 + a 23 t 3 = b 2 a 31 t 1 + a 32 t 2 + a 33 t 3 = b 3 \begin{cases} a_{11}t_{1}+a_{12}t_{2}+a_{13}t_{3}=b_{1}\\ a_{21}t_{1}+a_{22}t_{2}+a_{23}t_{3}=b_{2}\\ a_{31}t_{1}+a_{32}t_{2}+a_{33}t_{3}=b_{3}\\ \end{cases} ⎩ ⎨ ⎧a11t1+a12t2+a13t3=b1a21t1+a22t2+a23t3=b2a31t1+a32t2+a33t3=b3
第一步：先将方程组化为显式形式：
{ t 1 = 1 a 11 ( b 1 − a 12 t 2 − a 13 t 3 ) t 2 = 1 a 22 ( b 2 − a 21 t 1 − a 23 t 3 ) t 3 = 1 a 33 ( b 3 − a 31 t 1 − a 32 t 2 ) \begin{cases} t_{1}=\frac{1}{a_{11}}(b_{1}-a_{12}t_{2}-a_{13}t_{3})\\ t_{2}=\frac{1}{a_{22}}(b_{2}-a_{21}t_{1}-a_{23}t_{3})\\ t_{3}=\frac{1}{a_{33}}(b_{3}-a_{31}t_{1}-a_{32}t_{2})\\ \end{cases} ⎩ ⎨ ⎧t1=a111(b1−a12t2−a13t3)t2=a221(b2−a21t1−a23t3)t3=a331(b3−a31t1−a32t2)
第二步：设定迭代初场(上标0表示第0次迭代)：
t 1 ( 0 ) 、 t 2 ( 0 ) 、 t 3 ( 0 ) t_{1}^{(0)}、t_{2}^{(0)}、t_{3}^{(0)} t1(0)、t2(0)、t3(0)
第一轮迭代：
{ t 2 ( 0 ) 、 t 3 ( 0 ) → t 1 ( 1 ) t 1 ( 1 ) 、 t 3 ( 0 ) → t 2 ( 1 ) t 1 ( 1 ) 、 t 2 ( 1 ) → t 3 ( 1 ) \begin{cases} t_{2}^{(0)}、t_{3}^{(0)}{\rightarrow}t_{1}^{(1)}\\ t_{1}^{(1)}、t_{3}^{(0)}{\rightarrow}t_{2}^{(1)}\\ t_{1}^{(1)}、t_{2}^{(1)}{\rightarrow}t_{3}^{(1)}\\ \end{cases} ⎩ ⎨ ⎧t2(0)、t3(0)→t1(1)t1(1)、t3(0)→t2(1)t1(1)、t2(1)→t3(1)
第二轮迭代：
{ t 2 ( 1 ) 、 t 3 ( 1 ) → t 1 ( 2 ) t 1 ( 2 ) 、 t 3 ( 1 ) → t 2 ( 2 ) t 1 ( 2 ) 、 t 2 ( 2 ) → t 3 ( 2 ) \begin{cases} t_{2}^{(1)}、t_{3}^{(1)}{\rightarrow}t_{1}^{(2)}\\ t_{1}^{(2)}、t_{3}^{(1)}{\rightarrow}t_{2}^{(2)}\\ t_{1}^{(2)}、t_{2}^{(2)}{\rightarrow}t_{3}^{(2)}\\ \end{cases} ⎩ ⎨ ⎧t2(1)、t3(1)→t1(2)t1(2)、t3(1)→t2(2)t1(2)、t2(2)→t3(2)
如此反复一直到前后两次迭代的差值在误差允许范围内，即收敛到误差允许范围内停止迭代，此时便得出了各个未知数的数值解（本次编程实现的算法中用到的是最简单的以前后两次迭代的差值的绝对值不小于给定的误差作为收敛依据）。
迭代过程是否可以收敛可根据以下的主对角线占优条件来判断：{ ∣ a 12 + a 13 ∣ ≤ ∣ a 11 ∣ ∣ a 21 + a 23 ∣ ≤ ∣ a 22 ∣ ∣ a 31 + a 32 ∣ ≤ ∣ a 33 ∣ \begin{cases} \vert{a_{12}}+a_{13}\vert\leq\vert{a_{11}}\vert\\ \vert{a_{21}}+a_{23}\vert\leq\vert{a_{22}}\vert\\ \vert{a_{31}}+a_{32}\vert\leq\vert{a_{33}}\vert\\ \end{cases} ⎩ ⎨ ⎧∣a12+a13∣≤∣a11∣∣a21+a23∣≤∣a22∣∣a31+a32∣≤∣a33∣
值得注意的是，在用热平衡法导出的差分方程，若每一个方程都选用导出该方程的中心节点的温度作为迭代变量，则上述条件必满足，迭代一定收敛。
(4) 显化节点方程
(4.1) 等温边界条件
由于具有对称性，只取墙角的四分之一进行研究。在等温边界条件下，对墙体等间隔分为一系列节点并建立坐标系如图所示（这里声明一下，由于节点数对(m,n)与二维数组的索引(i,j)对第二个参数的变化方向不一样，所以最终用MATLAB画出来的墙角温度分布图应该是整个墙的左上四分之一部分），有A、B、C一共3类节点，节点热平衡方程经过显化后分别为：
A ： t m , n = ( Δ y ) 2 ( t m + 1 , n + t m − 1 , n ) + 2 ( Δ x ) 2 t m , n − 1 2 [ Δ x ) 2 + Δ y ) 2 ] A：t_{m,n}=\frac{(\Delta{y})^{2}(t_{m+1,n}+t_{m-1,n})+2(\Delta{x})^{2}t_{m,n-1}}{2[\Delta{x})^{2}+\Delta{y})^{2}]} A：tm,n=2[Δx)2+Δy)2](Δy)2(tm+1,n+tm−1,n)+2(Δx)2tm,n−1
B ： t m , n = ( Δ y ) 2 ( t m + 1 , n + t m − 1 , n ) + ( Δ x ) 2 ( t m , n + 1 + t m , n − 1 ) 2 [ Δ x ) 2 + Δ y ) 2 ] B：t_{m,n}=\frac{(\Delta{y})^{2}(t_{m+1,n}+t_{m-1,n})+(\Delta{x})^{2}(t_{m,n+1}+t_{m,n-1})}{2[\Delta{x})^{2}+\Delta{y})^{2}]} B：tm,n=2[Δx)2+Δy)2](Δy)2(tm+1,n+tm−1,n)+(Δx)2(tm,n+1+tm,n−1)
C ： t m , n = ( Δ x ) 2 ( t m , n + 1 + t m , n − 1 ) + 2 ( Δ y ) 2 t m − 1 , n 2 [ Δ x ) 2 + Δ y ) 2 ] C：t_{m,n}=\frac{(\Delta{x})^{2}(t_{m,n+1}+t_{m,n-1})+2(\Delta{y})^{2}t_{m-1,n}}{2[\Delta{x})^{2}+\Delta{y})^{2}]} C：tm,n=2[Δx)2+Δy)2](Δx)2(tm,n+1+tm,n−1)+2(Δy)2tm−1,n
代入数据得：
A ： t m , n = t m + 1 , n + t m − 1 , n + 2 t m , n − 1 4 A：t_{m,n}=\frac{t_{m+1,n}+t_{m-1,n}+2t_{m,n-1}}{4} A：tm,n=4tm+1,n+tm−1,n+2tm,n−1
B ： t m , n = t m + 1 , n + t m − 1 , n + t m , n + 1 + t m , n − 1 4 B：t_{m,n}=\frac{t_{m+1,n}+t_{m-1,n}+t_{m,n+1}+t_{m,n-1}}{4} B：tm,n=4tm+1,n+tm−1,n+tm,n+1+tm,n−1
C ： t m , n = t m , n + 1 + t m , n − 1 + 2 t m − 1 , n 4 C：t_{m,n}=\frac{t_{m,n+1}+t_{m,n-1}+2t_{m-1,n}}{4} C：tm,n=4tm,n+1+tm,n−1+2tm−1,n
(4.2) 对流边界条件
在对流边界条件下，对墙体等间隔分为一系列节点并建立坐标系如图所示，有A、B、C、D、E、F、G、H、I、J、K、L、M一共13类节点，节点热平衡方程经过显化后分别为：：
A ： t m , n = h e Δ x ( Δ y ) 2 t f , e + λ ( Δ y ) 2 t m + 1 , n + λ ( Δ x ) 2 t m , n − 1 h e Δ x ( Δ y ) 2 + λ [ ( Δ x ) 2 + Δ y ) 2 ] A：t_{m,n}=\frac{h_{e}\Delta{x}(\Delta{y})^{2}t_{f,e}+\lambda(\Delta{y})^{2}t_{m+1,n}+\lambda(\Delta{x})^{2}t_{m,n-1}}{h_{e}\Delta{x}(\Delta{y})^{2}+\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} A：tm,n=heΔx(Δy)2+λ[(Δx)2+Δy)2]heΔx(Δy)2tf,e+λ(Δy)2tm+1,n+λ(Δx)2tm,n−1
B ： t m , n = h i Δ x ( Δ y ) 2 t f , i + λ ( Δ y ) 2 t m − 1 , n + λ ( Δ x ) 2 t m , n − 1 h i Δ x ( Δ y ) 2 + λ [ ( Δ x ) 2 + Δ y ) 2 ] B：t_{m,n}=\frac{h_{i}\Delta{x}(\Delta{y})^{2}t_{f,i}+\lambda(\Delta{y})^{2}t_{m-1,n}+\lambda(\Delta{x})^{2}t_{m,n-1}}{h_{i}\Delta{x}(\Delta{y})^{2}+\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} B：tm,n=hiΔx(Δy)2+λ[(Δx)2+Δy)2]hiΔx(Δy)2tf,i+λ(Δy)2tm−1,n+λ(Δx)2tm,n−1
C ： t m , n = h i Δ x Δ y ( Δ x + Δ y ) t f , i + 2 λ ( Δ x ) 2 t m , n − 1 + 2 λ ( Δ y ) 2 t m − 1 , n + λ ( Δ x ) 2 t m , n + 1 + λ ( Δ y ) 2 t m + 1 , n h i Δ x Δ y ( Δ x + Δ y ) + 3 λ [ ( Δ x ) 2 + Δ y ) 2 ] C：t_{m,n}=\frac{h_{i}\Delta{x}\Delta{y}(\Delta{x}+\Delta{y})t_{f,i}+2\lambda(\Delta{x})^{2}t_{m,n-1}+2\lambda(\Delta{y})^{2}t_{m-1,n}+\lambda(\Delta{x})^{2}t_{m,n+1}+\lambda(\Delta{y})^{2}t_{m+1,n}}{h_{i}\Delta{x}\Delta{y}(\Delta{x}+\Delta{y})+3\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} C：tm,n=hiΔxΔy(Δx+Δy)+3λ[(Δx)2+Δy)2]hiΔxΔy(Δx+Δy)tf,i+2λ(Δx)2tm,n−1+2λ(Δy)2tm−1,n+λ(Δx)2tm,n+1+λ(Δy)2tm+1,n
D ： t m , n = h i Δ y ( Δ x ) 2 t f , i + λ ( Δ y ) 2 t m − 1 , n + λ ( Δ x ) 2 t m , n − 1 h i Δ y ( Δ x ) 2 + λ [ ( Δ x ) 2 + Δ y ) 2 ] D：t_{m,n}=\frac{h_{i}\Delta{y}(\Delta{x})^{2}t_{f,i}+\lambda(\Delta{y})^{2}t_{m-1,n}+\lambda(\Delta{x})^{2}t_{m,n-1}}{h_{i}\Delta{y}(\Delta{x})^{2}+\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} D：tm,n=hiΔy(Δx)2+λ[(Δx)2+Δy)2]hiΔy(Δx)2tf,i+λ(Δy)2tm−1,n+λ(Δx)2tm,n−1
E ： t m , n = h e Δ y ( Δ x ) 2 t f , e + λ ( Δ y ) 2 t m − 1 , n + λ ( Δ x ) 2 t m , n + 1 h e Δ y ( Δ x ) 2 + λ [ ( Δ x ) 2 + Δ y ) 2 ] E：t_{m,n}=\frac{h_{e}\Delta{y}(\Delta{x})^{2}t_{f,e}+\lambda(\Delta{y})^{2}t_{m-1,n}+\lambda(\Delta{x})^{2}t_{m,n+1}}{h_{e}\Delta{y}(\Delta{x})^{2}+\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} E：tm,n=heΔy(Δx)2+λ[(Δx)2+Δy)2]heΔy(Δx)2tf,e+λ(Δy)2tm−1,n+λ(Δx)2tm,n+1
F ： t m , n = h e Δ x Δ y ( Δ x + Δ y ) t f , e + λ ( Δ x ) 2 t m , n + 1 + λ ( Δ y ) 2 t m + 1 , n h e Δ x Δ y ( Δ x + Δ y ) + λ [ ( Δ x ) 2 + Δ y ) 2 ] F：t_{m,n}=\frac{h_{e}\Delta{x}\Delta{y}(\Delta{x}+\Delta{y})t_{f,e}+\lambda(\Delta{x})^{2}t_{m,n+1}+\lambda(\Delta{y})^{2}t_{m+1,n}}{h_{e}\Delta{x}\Delta{y}(\Delta{x}+\Delta{y})+\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} F：tm,n=heΔxΔy(Δx+Δy)+λ[(Δx)2+Δy)2]heΔxΔy(Δx+Δy)tf,e+λ(Δx)2tm,n+1+λ(Δy)2tm+1,n
G ： t m , n = ( Δ y ) 2 ( t m + 1 , n + t m − 1 , n ) + 2 ( Δ x ) 2 t m , n − 1 2 [ Δ x ) 2 + Δ y ) 2 ] G：t_{m,n}=\frac{(\Delta{y})^{2}(t_{m+1,n}+t_{m-1,n})+2(\Delta{x})^{2}t_{m,n-1}}{2[\Delta{x})^{2}+\Delta{y})^{2}]} G：tm,n=2[Δx)2+Δy)2](Δy)2(tm+1,n+tm−1,n)+2(Δx)2tm,n−1
H ： t m , n = λ ( Δ x ) 2 ( t m , n + 1 + t m , n − 1 ) + 2 λ ( Δ y ) 2 t m − 1 , n + 2 h i Δ x ( Δ y ) 2 t f , i 2 h i Δ x ( Δ y ) 2 + 2 λ [ ( Δ x ) 2 + Δ y ) 2 ] H：t_{m,n}=\frac{\lambda(\Delta{x})^{2}(t_{m,n+1}+t_{m,n-1})+2\lambda(\Delta{y})^{2}t_{m-1,n}+2h_{i}\Delta{x}(\Delta{y})^{2}t_{f,i}}{2h_{i}\Delta{x}(\Delta{y})^{2}+2\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} H：tm,n=2hiΔx(Δy)2+2λ[(Δx)2+Δy)2]λ(Δx)2(tm,n+1+tm,n−1)+2λ(Δy)2tm−1,n+2hiΔx(Δy)2tf,i
I ： t m , n = λ ( Δ y ) 2 ( t m + 1 , n + t m − 1 , n ) + 2 λ ( Δ x ) 2 t m , n − 1 + 2 h i Δ y ( Δ x ) 2 t f , i 2 h i Δ y ( Δ x ) 2 + 2 λ [ ( Δ x ) 2 + Δ y ) 2 ] I：t_{m,n}=\frac{\lambda(\Delta{y})^{2}(t_{m+1,n}+t_{m-1,n})+2\lambda(\Delta{x})^{2}t_{m,n-1}+2h_{i}\Delta{y}(\Delta{x})^{2}t_{f,i}}{2h_{i}\Delta{y}(\Delta{x})^{2}+2\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} I：tm,n=2hiΔy(Δx)2+2λ[(Δx)2+Δy)2]λ(Δy)2(tm+1,n+tm−1,n)+2λ(Δx)2tm,n−1+2hiΔy(Δx)2tf,i
J ： t m , n = ( Δ x ) 2 ( t m , n + 1 + t m , n − 1 ) + 2 ( Δ y ) 2 t m − 1 , n 2 [ Δ x ) 2 + Δ y ) 2 ] J：t_{m,n}=\frac{(\Delta{x})^{2}(t_{m,n+1}+t_{m,n-1})+2(\Delta{y})^{2}t_{m-1,n}}{2[\Delta{x})^{2}+\Delta{y})^{2}]} J：tm,n=2[Δx)2+Δy)2](Δx)2(tm,n+1+tm,n−1)+2(Δy)2tm−1,n
K ： t m , n = λ ( Δ y ) 2 ( t m + 1 , n + t m − 1 , n ) + 2 λ ( Δ x ) 2 t m , n + 1 + 2 h e Δ y ( Δ x ) 2 t f , e 2 h e Δ y ( Δ x ) 2 + 2 λ [ ( Δ x ) 2 + Δ y ) 2 ] K：t_{m,n}=\frac{\lambda(\Delta{y})^{2}(t_{m+1,n}+t_{m-1,n})+2\lambda(\Delta{x})^{2}t_{m,n+1}+2h_{e}\Delta{y}(\Delta{x})^{2}t_{f,e}}{2h_{e}\Delta{y}(\Delta{x})^{2}+2\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} K：tm,n=2heΔy(Δx)2+2λ[(Δx)2+Δy)2]λ(Δy)2(tm+1,n+tm−1,n)+2λ(Δx)2tm,n+1+2heΔy(Δx)2tf,e
L ： t m , n = λ ( Δ x ) 2 ( t m , n + 1 + t m , n − 1 ) + 2 λ ( Δ y ) 2 t m + 1 , n + 2 h e Δ x ( Δ y ) 2 t f , e 2 h e Δ x ( Δ y ) 2 + 2 λ [ ( Δ x ) 2 + Δ y ) 2 ] L：t_{m,n}=\frac{\lambda(\Delta{x})^{2}(t_{m,n+1}+t_{m,n-1})+2\lambda(\Delta{y})^{2}t_{m+1,n}+2h_{e}\Delta{x}(\Delta{y})^{2}t_{f,e}}{2h_{e}\Delta{x}(\Delta{y})^{2}+2\lambda[(\Delta{x})^{2}+\Delta{y})^{2}]} L：tm,n=2heΔx(Δy)2+2λ[(Δx)2+Δy)2]λ(Δx)2(tm,n+1+tm,n−1)+2λ(Δy)2tm+1,n+2heΔx(Δy)2tf,e
M ： t m , n = ( Δ y ) 2 ( t m + 1 , n + t m − 1 , n ) + ( Δ x ) 2 ( t m , n + 1 + t m , n − 1 ) 2 [ Δ x ) 2 + Δ y ) 2 ] M：t_{m,n}=\frac{(\Delta{y})^{2}(t_{m+1,n}+t_{m-1,n})+(\Delta{x})^{2}(t_{m,n+1}+t_{m,n-1})}{2[\Delta{x})^{2}+\Delta{y})^{2}]} M：tm,n=2[Δx)2+Δy)2](Δy)2(tm+1,n+tm−1,n)+(Δx)2(tm,n+1+tm,n−1)
代入数据得：
A ： t m , n = 30 + 0.53 ( t m + 1 , n + t m , n − 1 ) 2.06 A：t_{m,n}=\frac{30+0.53(t_{m+1,n}+t_{m,n-1})}{2.06} A：tm,n=2.0630+0.53(tm+1,n+tm,n−1)
B ： t m , n = 4 + 0.53 ( t m − 1 , n + t m , n − 1 ) 1.46 B：t_{m,n}=\frac{4+0.53(t_{m-1,n}+t_{m,n-1})}{1.46} B：tm,n=1.464+0.53(tm−1,n+tm,n−1)
C ： t m , n = 8 + 1.06 ( t m , n − 1 + t m − 1 , n ) + 0.53 ( t m , n + 1 + t m + 1 , n ) 3.98 C：t_{m,n}=\frac{8+1.06(t_{m,n-1}+t_{m-1,n})+0.53(t_{m,n+1}+t_{m+1,n})}{3.98} C：tm,n=3.988+1.06(tm,n−1+tm−1,n)+0.53(tm,n+1+tm+1,n)
D ： t m , n = 4 + 0.53 ( t m − 1 , n + t m , n − 1 ) 1.46 D：t_{m,n}=\frac{4+0.53(t_{m-1,n}+t_{m,n-1})}{1.46} D：tm,n=1.464+0.53(tm−1,n+tm,n−1)
E ： t m , n = 30 + 0.53 ( t m − 1 , n + t m , n + 1 ) 2.06 E：t_{m,n}=\frac{30+0.53(t_{m-1,n}+t_{m,n+1})}{2.06} E：tm,n=2.0630+0.53(tm−1,n+tm,n+1)
F ： t m , n = 60 + 0.53 ( t m , n + 1 + t m + 1 , n ) 3.06 F：t_{m,n}=\frac{60+0.53(t_{m,n+1}+t_{m+1,n})}{3.06} F：tm,n=3.0660+0.53(tm,n+1+tm+1,n)
G ： t m , n = t m + 1 , n + t m − 1 , n + 2 t m , n − 1 4 G：t_{m,n}=\frac{t_{m+1,n}+t_{m-1,n}+2t_{m,n-1}}{4} G：tm,n=4tm+1,n+tm−1,n+2tm,n−1
H ： t m , n = 0.53 ( t m , n + 1 + t m , n − 1 ) + 1.06 t m − 1 , n + 8 2.92 H：t_{m,n}=\frac{0.53(t_{m,n+1}+t_{m,n-1})+1.06t_{m-1,n}+8}{2.92} H：tm,n=2.920.53(tm,n+1+tm,n−1)+1.06tm−1,n+8
I ： t m , n = 0.53 ( t m + 1 , n + t m − 1 , n ) + 1.06 t m , n − 1 + 8 2.92 I：t_{m,n}=\frac{0.53(t_{m+1,n}+t_{m-1,n})+1.06t_{m,n-1}+8}{2.92} I：tm,n=2.920.53(tm+1,n+tm−1,n)+1.06tm,n−1+8
J ： t m , n = t m , n + 1 + t m , n − 1 + 2 t m − 1 , n 4 J：t_{m,n}=\frac{t_{m,n+1}+t_{m,n-1}+2t_{m-1,n}}{4} J：tm,n=4tm,n+1+tm,n−1+2tm−1,n
K ： t m , n = 0.53 ( t m + 1 , n + t m − 1 , n ) + 1.06 t m , n + 1 + 60 4.12 K：t_{m,n}=\frac{0.53(t_{m+1,n}+t_{m-1,n})+1.06t_{m,n+1}+60}{4.12} K：tm,n=4.120.53(tm+1,n+tm−1,n)+1.06tm,n+1+60
L ： t m , n = 0.53 ( t m , n + 1 + t m , n − 1 ) + 1.06 t m + 1 , n + 60 4.12 L：t_{m,n}=\frac{0.53(t_{m,n+1}+t_{m,n-1})+1.06t_{m+1,n}+60}{4.12} L：tm,n=4.120.53(tm,n+1+tm,n−1)+1.06tm+1,n+60
M ： t m , n = t m + 1 , n + t m − 1 , n + t m , n + 1 + t m , n − 1 4 M：t_{m,n}=\frac{t_{m+1,n}+t_{m-1,n}+t_{m,n+1}+t_{m,n-1}}{4} M：tm,n=4tm+1,n+tm−1,n+tm,n+1+tm,n−1
(5) FORTRAN迭代
PROGRAM main
IMPLICIT NONE
TYPE :: node !定义节点对象
REAL :: x !节点横坐标(m)
REAL :: y !节点纵坐标(m)
REAL :: t !节点温度(℃)
LOGICAL :: is_in_wall !判断节点是否在墙体内
END TYPE node
REAL,PARAMETER :: dx = 0.1,dy = 0.1 !网格步长(m)
REAL,PARAMETER :: e = 1.e-3 !给定允许误差
INTEGER,PARAMETER :: m_len = 15,n_len = 11 !四分之一墙角节点的(m,n)数对的m、n范围
INTEGER,PARAMETER :: m_inner_corner = 5,n_inner_corner = 5 !四分之一墙角内角点的(m,n)数对
REAL,PARAMETER :: t_w_e = 30,t_w_i = 0 !内外表面温度(℃)
REAL,PARAMETER :: t_f_e = 30,t_f_i = 10 !内外流体温度(℃)
REAL,PARAMETER :: h_e = 10,h_i = 4 !对流换热系数(W/(m^2·K))
INTEGER :: i,j !循环变量
REAL :: dt !初场相邻等温线之间的温差(℃)
TYPE(node),DIMENSION(0:m_len,0:n_len) :: dengwen_nodes,duiliu_nodes,temp_nodes !节点数组
!节点初始化
DO j = 0,n_len
DO i = 0,m_len
dengwen_nodes(i,j)%t = -273.15
duiliu_nodes(i,j)%t = -273.15
IF(i > m_inner_corner .AND. j > n_inner_corner) THEN
dengwen_nodes(i,j)%x = dx*i
dengwen_nodes(i,j)%y = dy*j
dengwen_nodes(i,j)%is_in_wall = .FALSE.
duiliu_nodes(i,j)%x = dx*i
duiliu_nodes(i,j)%y = dy*j
duiliu_nodes(i,j)%is_in_wall = .FALSE.
ELSE
dengwen_nodes(i,j)%x = dx*i
dengwen_nodes(i,j)%y = dy*j
dengwen_nodes(i,j)%is_in_wall = .TRUE.
duiliu_nodes(i,j)%x = dx*i
duiliu_nodes(i,j)%y = dy*j
duiliu_nodes(i,j)%is_in_wall = .TRUE.
END IF
END DO
END DO
!等温边界条件
!设定初场
dt = ABS(t_w_e - t_w_i)/m_inner_corner !或者除以n_inner_corner
dengwen_nodes(0,:)%t = t_w_e
dengwen_nodes(1:,0)%t = t_w_e
DO i = 1,m_inner_corner
dengwen_nodes(i,1:)%t = t_w_e - i*dt
END DO
DO j = 1,n_inner_corner
dengwen_nodes(m_inner_corner + 1:,j)%t = t_w_e - j*dt
END DO
!开始迭代
100 DO WHILE(.TRUE.)
temp_nodes = dengwen_nodes
!A类节点
DO i = 1,m_inner_corner - 1
dengwen_nodes(i,n_len)%t = (dengwen_nodes(i+1,n_len)%t+dengwen_nodes(i-1,n_len)%t+2*dengwen_nodes(i,n_len-1)%t)/4
END DO
!B类节点第一部分
DO j = 1,n_inner_corner - 1
DO i = 1,m_len - 1
dengwen_nodes(i,j)%t = (dengwen_nodes(i+1,j)%t+dengwen_nodes(i-1,j)%t+dengwen_nodes(i,j+1)%t+dengwen_nodes(i,j-1)%t)/4
END DO
END DO
!B类节点第二部分
DO j = n_inner_corner,n_len - 1
DO i = 1,m_inner_corner - 1
dengwen_nodes(i,j)%t = (dengwen_nodes(i+1,j)%t+dengwen_nodes(i-1,j)%t+dengwen_nodes(i,j+1)%t+dengwen_nodes(i,j-1)%t)/4
END DO
END DO
!C类节点
DO j = 1,n_inner_corner - 1
dengwen_nodes(i,j)%t = (dengwen_nodes(m_len,j+1)%t+dengwen_nodes(m_len,j-1)%t+2*dengwen_nodes(m_len-1,j)%t)/4
END DO
DO j = 0,n_len
DO i = 0,m_len
IF(ABS(temp_nodes(i,j)%t-dengwen_nodes(i,j)%t)>e .AND. dengwen_nodes(i,j)%is_in_wall) GOTO 100
END DO
END DO
EXIT
END DO
!对流边界条件
!设定初场
dt = ABS(t_f_e - t_f_i)/m_inner_corner !或者除以n_inner_corner
duiliu_nodes(0,:)%t = t_f_e
duiliu_nodes(1:,0)%t = t_f_e
DO i = 1,m_inner_corner
duiliu_nodes(i,1:)%t = t_f_e - i*dt
END DO
DO j = 1,n_inner_corner
duiliu_nodes(m_inner_corner + 1:,j)%t = t_f_e - j*dt
END DO
!开始迭代
200 DO WHILE(.TRUE.)
temp_nodes = duiliu_nodes
!A类节点
duiliu_nodes(0,n_len)%t = (30+0.53*(duiliu_nodes(1,n_len)%t+duiliu_nodes(0,n_len-1)%t))/2.06
!B类节点
duiliu_nodes(m_inner_corner,n_len)%t = (4+0.53*(duiliu_nodes(m_inner_corner-1,n_len)%t+duiliu_nodes(m_inner_corner,n_len-1)%t))/1.46
!C类节点
duiliu_nodes(m_inner_corner,n_inner_corner)%t = (8+1.06*(duiliu_nodes(m_inner_corner,n_inner_corner-1)%t+duiliu_nodes(m_inner_corner-1,n_inner_corner)%t)+0.53*(duiliu_nodes(m_inner_corner,n_inner_corner+1)%t+duiliu_nodes(m_inner_corner+1,n_inner_corner)%t))/3.98
!D类节点
duiliu_nodes(m_len,n_inner_corner)%t = (4+0.53*(duiliu_nodes(m_len-1,n_inner_corner)%t+duiliu_nodes(m_len,n_inner_corner-1)%t))/1.46
!E类节点
duiliu_nodes(m_len,0)%t = (30+0.53*(duiliu_nodes(m_len-1,0)%t+duiliu_nodes(m_len,1)%t))/2.06
!F类节点
duiliu_nodes(0,0)%t = (60+0.53*(duiliu_nodes(0,1)%t+duiliu_nodes(1,0)%t))/3.06
!G类节点
do i = 1,m_inner_corner - 1
duiliu_nodes(i,n_len)%t = (duiliu_nodes(i+1,n_len)%t+duiliu_nodes(i-1,n_len)%t+2*duiliu_nodes(i,n_len-1)%t)/4
END DO
!H类节点
DO j = n_inner_corner+1,n_len - 1
duiliu_nodes(m_inner_corner,j)%t = (0.53*(duiliu_nodes(m_inner_corner,j+1)%t+duiliu_nodes(m_inner_corner,j-1)%t)+1.06*duiliu_nodes(m_inner_corner-1,j)%t+8)/2.92
END DO
!I类节点
DO i = m_inner_corner + 1,m_len - 1
duiliu_nodes(i,n_inner_corner)%t = (0.53*(duiliu_nodes(i+1,n_inner_corner)%t+duiliu_nodes(i-1,n_inner_corner)%t)+1.06*duiliu_nodes(i,n_inner_corner-1)%t+8)/2.92
END DO
!J类节点
DO j = 1,n_inner_corner - 1
duiliu_nodes(m_len,j)%t = (duiliu_nodes(m_len,j+1)%t+duiliu_nodes(m_len,j-1)%t+2*duiliu_nodes(m_len-1,j)%t)/4
END DO
!K类节点
DO i = 1,m_len - 1
duiliu_nodes(i,0)%t = (0.53*(duiliu_nodes(i+1,0)%t+duiliu_nodes(i-1,0)%t)+1.06*duiliu_nodes(i,1)%t+60)/4.12
END DO
!L类节点
DO j = 1,n_len - 1
duiliu_nodes(0,j)%t = (0.53*(duiliu_nodes(0,j+1)%t+duiliu_nodes(0,j-1)%t)+1.06*duiliu_nodes(1,j)%t+60)/4.12
END DO
!M类节点第一部分
DO j = 1,n_len - 1
DO i = 1,m_inner_corner - 1
duiliu_nodes(i,j)%t = (duiliu_nodes(i+1,j)%t+duiliu_nodes(i-1,j)%t+duiliu_nodes(i,j+1)%t+duiliu_nodes(i,j-1)%t)/4
END DO
END DO
!M类节点第二部分
DO j = 1,n_inner_corner - 1
DO i = m_inner_corner,m_len - 1
duiliu_nodes(i,j)%t = (duiliu_nodes(i+1,j)%t+duiliu_nodes(i-1,j)%t+duiliu_nodes(i,j+1)%t+duiliu_nodes(i,j-1)%t)/4
END DO
END DO
DO j = 0,n_len
DO i = 0,m_len
IF(ABS(temp_nodes(i,j)%t-duiliu_nodes(i,j)%t)>e .AND. duiliu_nodes(i,j)%is_in_wall) GOTO 200
END DO
END DO
EXIT
END DO
!在屏幕上输出计算结果
WRITE(*,*) '每个三维数对表示：（节点横坐标(m)，节点纵坐标(m)，节点温度(℃)）'
WRITE(*,*) '等温边界条件下墙角温度数值解：'
DO j = n_len,n_inner_corner + 1,-1
DO i = 0,m_inner_corner
WRITE(*,"('(',f4.1,',',f4.1,',',f5.2,')',4x,$)") dengwen_nodes(i,j)%x,dengwen_nodes(i,j)%y,dengwen_nodes(i,j)%t
END DO
WRITE(*,"(/)")
END DO
DO j = n_inner_corner,0,-1
DO i = 0,m_len
WRITE(*,"('(',f4.1,',',f4.1,',',f5.2,')',4x,$)") dengwen_nodes(i,j)%x,dengwen_nodes(i,j)%y,dengwen_nodes(i,j)%t
END DO
WRITE(*,"(/)")
END DO
WRITE(*,*) '对流边界条件下墙角温度数值解：'
DO j = n_len,n_inner_corner + 1,-1
DO i = 0,m_inner_corner
WRITE(*,"('(',f4.1,',',f4.1,',',f5.2,')',4x,$)") duiliu_nodes(i,j)%x,duiliu_nodes(i,j)%y,duiliu_nodes(i,j)%t
END DO
WRITE(*,"(/)")
END DO
DO j = n_inner_corner,0,-1
DO i = 0,m_len
WRITE(*,"('(',f4.1,',',f4.1,',',f5.2,')',4x,$)") duiliu_nodes(i,j)%x,duiliu_nodes(i,j)%y,duiliu_nodes(i,j)%t
END DO
WRITE(*,"(/)")
END DO
!将所得数据写入.csv文件，后续用MATLAB提取并进行数据可视化处理
OPEN(100,FILE = 'e:/coordinate_x.csv',STATUS = 'replace')
DO i = 0,m_len
DO j = 0,n_len
IF(dengwen_nodes(i,j)%is_in_wall) WRITE(100,"(f4.1)") dengwen_nodes(i,j)%x
END DO
END DO
CLOSE(100)
OPEN(200,FILE = 'e:/coordinate_y.csv',STATUS = 'replace')
DO i = 0,m_len
DO j = 0,n_len
IF(dengwen_nodes(i,j)%is_in_wall) WRITE(200,"(f4.1)") dengwen_nodes(i,j)%y
END DO
WRITE(200,*)
END DO
CLOSE(200)
OPEN(300,FILE = 'e:/dengwen_t.csv',STATUS = 'replace')
DO i = 0,m_len
DO j = 0,n_len
IF(dengwen_nodes(i,j)%is_in_wall) WRITE(300,"(f8.3)") dengwen_nodes(i,j)%t
END DO
WRITE(300,*)
END DO
CLOSE(300)
OPEN(400,FILE = 'e:/duiliu_t.csv',STATUS = 'replace')
DO i = 0,m_len
DO j = 0,n_len
IF(duiliu_nodes(i,j)%is_in_wall) WRITE(400,"(f8.3)") duiliu_nodes(i,j)%t
END DO
WRITE(400,*)
END DO
CLOSE(400)
END PROGRAM main
程序运行后会得到四个.scv文件，分别记录节点横坐标、节点纵坐标、等温边界条件下节点温度、对流边界条件下节点温度。
(6) MATLAB数据可视化
clc;clear;
dx = 0.1;dy = 0.1; %网格步长(m)
lamda = 0.53; %导热系数(W/(m·K))
t_f_i = 10;t_f_e = 30; %流体温度(℃)
h_i = 4;h_e = 10; %对流换热系数(W/(m^2·K))
dengwen_q_e = 0;dengwen_q_i = 0;dengwen_q = 0; %等温边界条件换热量(W/m)
duiliu_q_e = 0;duiliu_q_i = 0;duiliu_q = 0; %对流边界条件换热量(W/m)
%读取文件数据
x = csvread('e:\coordinate_x.csv');
y = csvread('e:\coordinate_y.csv');
dw_t = csvread('e:\dengwen_t.csv');
dl_t = csvread('e:\duiliu_t.csv');
[row,column] = size(x);
%数据可视化
colordef none; %个人比较喜欢暗黑风格
xi = 0:0.01:max(x); %网格加密
yi = 0:0.01:max(y);
[X,Y] = meshgrid(xi,yi); %生成矩阵网格
dw_T = griddata(x,y,dw_t,X,Y,'cubic'); %插值
dl_T = griddata(x,y,dl_t,X,Y,'cubic');
%用二维图形表示温度场
figure;
subplot(121);
contour(X,Y,dw_T);
colorbar('southoutside');
xlabel('x/m');
ylabel('y/m');
title('墙角等温线图(等温边界条件)');
subplot(122);
pcolor(X,Y,dw_T);
colorbar('southoutside');
xlabel('x/m');
ylabel('y/m');
title('墙角色温图(等温边界条件)');
figure;
subplot(121);
contour(X,Y,dl_T);
colorbar('southoutside');
xlabel('x/m');
ylabel('y/m');
title('墙角等温线图(对流边界条件)');
subplot(122);
pcolor(X,Y,dl_T);
colorbar('southoutside');
xlabel('x/m');
ylabel('y/m');
title('墙角色温图(对流边界条件)');
% 三维图形表示温度场
figure;
subplot(121);
surfc(X,Y,dw_T);
xlabel('x/m');
ylabel('y/m');
zlabel('temperature/℃');
title('三维墙角温度场分布(等温边界条件)');
colorbar('southoutside');
subplot(122);
surfc(X,Y,dl_T);
xlabel('x/m');
ylabel('y/m');
zlabel('temperature/℃');
title('三维墙角温度场分布(对流边界条件)');
colorbar('southoutside');
参考资料来源 参考资料来源 参考资料来源
- 《传热学(第五版)》.陶文铨.杨世铭.高等教育出版社.
- 《热与流体实验教程》.王小丹.孟倩.张可.吴青平.唐上朝.西安交通大学出版社.
源码作者： A i d e n L e e 源码作者：Aiden\ Lee 源码作者：Aiden Lee
博客创作： A i d e n L e e 博客创作：Aiden\ Lee 博客创作：Aiden Lee
特别声明：文章仅供学习参考，转载请注明出处，严禁盗用！
扫一扫
135
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
