基于雅克比迭代算法室内温度分布问题的研究——封闭空间内点热源的散热情况

Research on Indoor Temperature Distribution Based on Jacobi Iterative Algorithm—Heat Dissipation of Point Heat Source in Closed Space

摘要: 本文章在给定环境下，求出室内的稳态、动态温度分布。通过傅立叶导热定律，采用热平衡法建立区域节点离散方程，利用雅克比迭代算法进行分析求解。通过围护结构耗热量公式计算出房间的基本耗热量，推导以房间温度为变量的能量守恒方程，求得出水口水温，然后确定定解条件，建立稳态导热内节点离散方程，并利用雅克比迭代的方法表达连续介质各点在同一时刻的温度分布。把进水口流量分成等差值情况进行分析讨论，在不同流量下，依据能量守恒方程，便可通过假设出水口温度求出室内平均温度，进而分别求出每种流量下的室内温度分布。

Abstract: In this paper, the steady and dynamic temperature distribution in a given environment is solved. According to Fourier heat conduction law, the discrete equation of regional nodes is established by heat balance method, and the Jacobian iterative algorithm is used for analysis and solution. The basic heat consumption of the room is calculated by the heat consumption formula of the enclosure structure, the energy conservation equation with the room temperature as the variable is derived, and the water temperature of the nozzle is obtained. Then, the definite solution conditions are determined, and the discrete equation of the internal nodes of the steady-state heat conduction is established. The temperature distribution of each point in the continuous medium at the same time is expressed by the Jacobian iteration method. According to the energy conservation equation, the indoor average temperature can be calculated by assuming the outlet temperature under different flow rates, and then the indoor temperature distribution under each flow rate can be calculated respectively. The numerical solution and temperature contrast chart of some temperature regions under different flow rates are given.

1. 引言

房间内的温度对于普通百姓而言，温度大小的变化对他们的影响不是很大，没有强制性要求。但对于某些特殊行业(比如医药行业)而言，室温是有严格规定的，必须按照行业准则进行调控室内温度。在给定时刻的房间内温度随空间位置的变化，这是热量传递研究中要解决的主要问题之一。在换热器、反应器等化工设备中，由于两流体存在温差，或因反应具有热效应，形成温度分布。

许多城市冬季供暖采用热水供热的方式。每一户有一个进水口和出水口，热水通过进水口进入房间，由出水口离开房间。一个基本问题是，如何调整进水流量，使得在室外天气不断变化时，室内的温度维持在一个相对稳定的范围内。

2. 建立模型

假设房间大小为50平方米，房间高为2.8米，具体规格如图所示。房间南北为外墙，与室外接触，长为5米，每面墙上有一扇窗户，规格为2米 * 1.5米；东西为内墙，不与室外接触，长为10米。房间采用暖气片供热方式，共两扇，南北墙窗户正下方各一扇，每扇暖气片长2米，高1米。

房间进水口的流量介于每小时0.1到1立方米之间。进水口的水温为40度，出水口的水温未知。房间的室内温度标准为每天的平均温度为21度。

本文研究当室外温度为0度，进水口流量为每小时0.5立方米时，房间的空间温度分布是怎样的？当室外温度为0度，进水口流量在0.1~1之间变化时，房间的动态温度分布是怎样的？

2.1. 模型的研究

2.1.1. 温度场

物体各部分之间不发生相对位移时，依靠分子、原子及自由电子等微观粒子的热运动进行的热量传递称为热传导。在紧密的不透明物体内部，热量只能依靠热传导方式传递。温度场是指某时刻空间内所有各点温度的总称，若温度场随着时间变化，则称作稳定温度场，由此产生的导热为稳定导热；若温度场随时间改变，则称作不稳定温度场，不稳定温度场的导热为不稳定导热。温度场是温度有关时间和空间的函数，即：

三维非稳态温度场： $t=f\left(x,y,z,\tau \right)$

三维稳态温度场： $t=f\left(x,y,z\right)$

二维稳态温度场： $t=f\left(x,y\right)$

一维稳态温度场： $t=f(\; x\; )$

2.1.2. 基于傅立叶定律下温度分布的研究

傅立叶导热定律，理论上讲，它给出的只是在已知温度梯度的某一个坐标向梯度方向传递热流的大小，是一个普适规律，一个物体内部任何一点的导热都要满足傅立叶定律，但整体上的能量传递效应，则需要热量守恒定律来管控。那么，对于三维直角坐标系下的某个坐标具有三个方向，即

$q\left(x\right)=-\lambda \frac{\partial t}{\partial x}$

$q\left(y\right)=-\lambda \frac{\partial t}{\partial y}$

$q\left(z\right)=-\lambda \frac{\partial t}{\partial z}$

温度梯度如图1所示。

Figure 1. Temperature gradient diagram

图1. 温度梯度图

三个坐标方向的单位矢量与该方向上温度梯度(热流密度)分量的乘积合成梯度矢量(热流密度矢量)

$gradt=\Delta t=\frac{\partial t}{\partial x}i+\frac{\partial t}{\partial y}j+\frac{\partial t}{\partial z}k$

$q=-\lambda \left(\frac{\partial t}{\partial x}i+\frac{\partial t}{\partial y}j+\frac{\partial t}{\partial z}k\right)=-\lambda \Delta t=-\lambda \frac{\partial t}{\partial n}n$

这个热流密度矢量就是一般形式的傅立叶定律。

傅立叶导热微分方程的一般形式：

$\rho c\frac{\partial t}{\partial \tau }=\frac{\partial }{\partial x}\left(\lambda \frac{\partial t}{\partial x}\right)+\frac{\partial t}{\partial }\left(\lambda \frac{\partial t}{\partial y}\right)+\frac{\partial }{\partial z}\left(\lambda \frac{\partial t}{\partial z}\right)+\stackrel{\dot{}}{\Phi }$

可知，运用该方程可求解房间内温度分布的通解，无法得到房间的具体位置的温度，方程也可表示为：

$\frac{\partial t}{\partial \tau }=\frac{\lambda }{\rho c}\left(\frac{{\partial }^{2}t}{\partial x^2}+\frac{{\partial }^{2}t}{\partial y^2}+\frac{{\partial }^{2}t}{\partial z^2}\right)+\frac{\stackrel{\dot{}}{\Phi }}{\rho c}$

其中 $\frac{\Phi }{\rho c}$ 是内热源项， $\frac{\lambda }{\rho c}$ 称为热扩散系数， $\frac{\partial t}{\partial \tau }$ 称作非稳态项。

2.2. 模型的建立

1) 结合实际生活可知，靠近暖气片的位置温度较高，随着距离的增加温度呈阶梯状减小。在分析暖气片散热对室内温度的影响情况时，假设房间是一个无内热源的二维空腔。由于本题要求解决房间的温度分布问题，考虑到其室内温度不仅仅与暖气片散热有关，还与房间内墙、外墙、外窗等围护结构的基本耗热情况有关。因围护结构的耗热消耗了室内的热量，故要把围护结构的基本耗热量考虑到解决吸收暖气片散热问题中去，具体的表达式为：

$暖气片的散热量\Delta Q=房间围护结构的基本耗热量{\displaystyle \sum Q}$

${\displaystyle \sum Q}=Q_{南外墙}+Q_{北外墙}+Q_{东西内墙}+Q_{南窗}+Q_{北窗}$

已知入水口温度为40℃，则可计算出水口水温：

查找围护结构传热系数表得：

外墙传热系数： $K_1=2.08\text{W}/\left({\text{m}}^2\cdot ˚\text{C}\right)$

外窗传热系数： $K_2=6.40\text{W}/\left({\text{m}}^2\cdot ˚\text{C}\right)$

内墙传热系数： $K_3=1.72\text{W}/\left({\text{m}}^2\cdot ˚\text{C}\right)$

南北朝向的外墙及外窗的温差修正系数： ${\alpha }_1=1$

东西朝向内墙温差修正系数： ${\alpha }_2=0.7$

取室外温度 $t_{wn}=0˚\text{C}$ ，室内计算温度 $t_n=21˚\text{C}$

南墙朝向修正率 ${\zeta }_1=1.5\text{\%}$

北墙朝向修正率 ${\zeta }_2=-2.5\%$

计算得：

$Q_{南外墙}=\left(1+1.5\%\right){\alpha }_1{K}_1{F}_1\left(t_n-t_{wn}\right)=\left(1+1.5\%\right)\times 1\times 2.08\times \left(2.8\times 5-2\times 1.5\right)\times \left(t_n-t_{wn}\right)$

$Q_{南外窗}=\left(1+1.5\%\right){\alpha }_1{K}_2{F}_2\left(t_n-t_{wn}\right)=\left(1+1.5\%\right)\times 1\times 6.4\times 2\times 1.5\times \left(t_n-t_{wn}\right)$

$Q_{北外墙}=\left(1-2.5\%\right){\alpha }_1{K}_1{F}_1\left(t_n-t_{wn}\right)=\left(1-2.5\%\right)\times 1\times 2.08\times \left(2.8\times 5-2\times 1.5\right)\times \left(t_n-t_{wn}\right)$

$Q_{北外窗}=\left(1-2.5\%\right){\alpha }_1{K}_2{F}_2\left(t_n-t_{wn}\right)=\left(1-2.5\%\right)\times 1\times 6.40\times 2\times 1.5\times \left(t_n-t_{wn}\right)$

$Q_{东西内墙}={\alpha }_2{K}_3{F}_3\left(t_n-t_{wn}\right)=0.7\times 1.72\times 10\times 2.5\times 2\times \left(t_n-t_{wn}\right)$

求得房间基本耗热量 ${\displaystyle \sum Q}=0.12\left(t_n-t_{wn}\right)\text{J}/\text{s}$

暖气片的散热量 $\Delta Q={\displaystyle \sum Q}$

根据水的散热量与温降关系，求得出水口水温和室内温差的热量守恒方程为：

$0.12\left(t_n-t_{wn}\right)=cm\Delta T$

整理得到热量守恒方程：

$\frac{0.12\left(t_n-t_{wn}\right)}{G}=40-T_{出}$

在问题中已确定入水口温度为40℃，先假设出水口温度为35℃并代入确定的热量守恒方程中便可得出房内的平均温度为20.8℃，将入水口温度为40℃出水口温度为35℃。再代入雅克比迭代中，并将迭代终止条件之一设定为计算平均温度与迭代的平均温度保持误差在足够小的范围内，便可求出在该定解条件下的室内平均温度以及各点的温度分布的数值解。

2) 坐标离散化

把房间区域进行坐标化离散化，以便于求解温度分布的数值解，对房间内的各点采用热平衡法建立离散方程并给出边界条件进行求解。

以二维矩形区域为例给出区域离散化过程，所谓区域离散化，就是将研究区域分解成有限数量的小区域，单元的顶点作为节点，每个节点都有自己的控制区域，称作控制体(控制容积)，控制体内所有特性都是均匀的，节点的温度代表其所在控制体的温度。节点之间的距离称为空间步长。节点之间的连线称为网格线，控制容积的分界面称为界面。对于本题采用的内节点法网格划分方式，单元的中心点即为节点，相邻单元中心线所围起来的区域就是控制体。

3) 离散方程的建立

空间区域离散化后，就要针对每个节点来建立其导热微分方程，也被称之为离散方程。对于本题采用控制容积平衡的方法(也称为热平衡法)建立离散方程。对每个有限大小的控制容积应用热量守恒分析，从而获得温度场的代数方程组，从基本的物理现象和基本定律出发，依据热量守恒和傅立叶导热定律即可获得离散方程。

对于任意节点的热量守恒都可以写成下式：

流入控制体的总热流量 + 控制体内热源生成热 = 控制体内能的增量

${\Phi }_i+{\Phi }_g={\Phi }_{st}$ 单位：[W]

流入控制体的总热流量 ${\Phi }_i$ ，是指从所有方向流入该控制体的热流量。以二维、稳态、有内热源的导热问题为例导出离散方程，假设导热是稳态的，所以从所有方向流入控制体的总热流量：

${\Phi }_i+{\Phi }_g=0$

4) 节点的选取

以房间地面的东南角为O点；东墙与地面的交线为y轴；南墙为地面的交线为x轴；东墙与南墙的交线为y轴建立空间直角坐标系，如图2所示。

Figure 2. Room coordinate system

图2. 房间坐标系

将这个长为十米宽为五米的房间划分成51 × 101的网格，即网格所划分的区域为50 × 100。对于每个区域的边长分别为 $\Delta x,\Delta y$ 其中 $\Delta x=\Delta y=0.1\text{m}$ 。

划分网格所产生了三种节点它们的离散方程分别为：

① 内节点

$t\left(m,n\right)=\left(t_0\times \left(m-1,n\right)+t_0\left(m+1,n\right)+t_0\left(m,n-1\right)+t_0\left(m,n+1\right)\right)/4$

② 四个角点

$t\left(1,1\right)=\left(t_0\times \left(2,1\right)+t_0\left(1,2\right)+2\times t_0\left(2,2\right)\right)/4$

$t\left(1,51\right)=\left(t_0\times \left(1,51\right)+t_0\left(2,51\right)+2\times t_0\left(2,50\right)\right)/4$

$t\left(101,1\right)=\left(t_0\times \left(100,1\right)+t_0\left(101,2\right)+2\times t_0\left(100,2\right)\right)/4$

$t\left(101,51\right)=\left(t_0\times \left(101,51\right)+t_0\left(101,50\right)+2\times t_0\left(100,50\right)\right)/4$

③ 边界节点

左边界节点： $t\left(a,b\right)=\left(t_0\times \left(a-1,b\right)+t_0\left(a+1,b\right)+2\times t_0\left(a,b+1\right)\right)/4$

右边界节点： $t\left(a,b\right)=\left(t_0\times \left(a-1,b\right)+t_0\left(a+1,b\right)+2\times t_0\left(a,b-1\right)\right)/4$

上边界节点： $t\left(a,b\right)=\left(t_0\times \left(a,b-1\right)+t_0\left(a,b+1\right)+2\times t_0\left(a+1,b\right)\right)/4$

下边界节点： $t\left(a,b\right)=\left(t_0\times \left(a,b-1\right)+t_0\left(a,b+1\right)+2\times t_0\left(a-1,b\right)\right)/4$

5) 节点列写具体方法：

对于 $\left(m,n\right)$ 节点来讲，与其相邻并具有传热面积的单元分别为 $\left(m-1,n\right)$ 、 $\left(m+1,n\right)$ 、 $\left(m,n-1\right)$ 和 $\left(m,n+1\right)$ 四个，温度节点分布情况如图3所示。

Figure 3. Temperature node distribution

图3. 温度节点分布情况

${\Phi }_i={\Phi }_{\left(m-1,n\right)}+{\Phi }_{\left(m+1,n\right)}+{\Phi }_{\left(m,n-1\right)}+{\Phi }_{\left(m,n+1\right)}$

故方程 ${\Phi }_i+{\Phi }_g=0$ 可以写成如下形式：

${\Phi }_i={\Phi }_{\left(m-1,n\right)}+{\Phi }_{\left(m+1,n\right)}+{\Phi }_{\left(m,n-1\right)}+{\Phi }_{\left(m,n+1\right)}+\Phi {}_g$

第一项 ${\Phi }_{\left(m-1,n\right)}$ 在单元 $\left(m,n\right)$ 的左侧，通过接触面积向右侧传递热流量时，用傅立叶定律计算：

${\Phi }_{\left(m-1,n\right)}=-\lambda A\frac{\text{d}t}{\text{d}x}=-\lambda \Delta y\frac{\text{d}t}{\text{d}x}$

当温度场不确定时， $\frac{\text{d}t}{\text{d}x}$ 未知，故先假设相邻节点的温度分布形式。并且当网格划分足够密集的时候，

两个节点间的距离、空间的步长足够小时，空间温度的分布形式对结果没有影响。假设温度呈分段线性分布，每个节点代表一个温度，数值计算得到的是离散的点的温度，把所有的点都用折线连接起来。当节点越多，假设的分段线性分布越接近真实的温度分布。傅立叶定律的表达式就可以表示成采用相邻节点的温差与坐标差的比值来代替温度的一阶导数。

同理，可得其他三个方向上流入控制体的热流量：

${\Phi }_{\left(m-1,n\right)}=-\lambda \Delta y\frac{\text{d}t}{\text{d}x}=-\lambda \Delta y\frac{t_{\left(m-1,n\right)}-t_{\left(m,n\right)}}{x_{\left(m-1,n\right)}-x_{\left(m,n\right)}}=\lambda \Delta y\frac{t_{\left(m-1,n\right)}-t_{\left(m,n\right)}}{\Delta x}$

${\Phi }_{\left(m+1,n\right)}=\lambda \Delta y\frac{t_{\left(m+1,n\right)}-t_{\left(m,n\right)}}{\Delta x}$

${\Phi }_{\left(m,n+1\right)}=\lambda \Delta x\frac{t_{\left(m,n+1\right)}-t_{\left(m,n\right)}}{\Delta y}$

${\Phi }_{\left(m,n-1\right)}=\lambda \Delta x\frac{t_{\left(m,n-1\right)}-t_{\left(m,n\right)}}{\Delta y}$

以上四项是从上下左右四个方向相邻节点导入到 $\left(m,n\right)$ 里面的热流量，内热源：

${\Phi }_g=\stackrel{\dot{}}{\Phi }\cdot V=\stackrel{\dot{}}{\Phi }\cdot \Delta x\Delta y$

把以上五项分别代入二维、稳态、有内热源的导热问题节点热量守恒方程：

${\Phi }_i={\Phi }_{\left(m-1,n\right)}+{\Phi }_{\left(m+1,n\right)}+{\Phi }_{\left(m,n-1\right)}+{\Phi }_{\left(m,n+1\right)}+{\Phi }_g$

即可得到针对节点 $\left(m,n\right)$ 的离散方程：

$\lambda \Delta y\frac{t_{\left(m-1,n\right)}-t_{\left(m,n\right)}}{\Delta x}+\lambda \Delta y\frac{t_{\left(m+1,n\right)}-t_{\left(m,n\right)}}{\Delta x}+\lambda \Delta x\frac{t_{\left(m,n+1\right)}-t_{\left(m,n\right)}}{\Delta y}+\lambda \Delta x\frac{t_{\left(m,n-1\right)}-t_{\left(m,n\right)}}{\Delta y}+\stackrel{\dot{}}{\Phi }\Delta x\Delta y=0$

该离散方程的每一项都有其明确的物理意义，当物性是常数且在x方向和y方向的空间步长均匀且相等的时候 $\Delta x=\Delta y$ 时，离散方程可以化简成：

$t_{\left(m-1,n\right)}+t_{\left(m+1,n\right)}+t_{\left(m,n+1\right)}+t_{\left(m,n-1\right)}-4t_{\left(m,n\right)}+\frac{\Delta x^2}{\lambda }\stackrel{\dot{}}{\Phi }=0$

$4t_{\left(m,n\right)}=t_{\left(m-1,n\right)}+t_{\left(m+1,n\right)}+t_{\left(m,n+1\right)}+t_{\left(m,n-1\right)}+\frac{\Delta x^2}{\lambda }\stackrel{\dot{}}{\Phi }$

无内热源时，上式可继续简化为：

$4t_{\left(m,n\right)}=t_{\left(m-1,n\right)}+t_{\left(m+1,n\right)}+t_{\left(m,n+1\right)}+t_{\left(m,n-1\right)}$

$t\left(m,n\right)=\left(t_0\times \left(m-1,n\right)+t_0\left(m+1,n\right)+t_0\left(m,n-1\right)+t_0\left(m,n+1\right)\right)/4$

采用热平衡的方法建立起来的离散方程，所求节点的温度前的系数等于其他所有相邻节点温度前的系数之和。这一结论也适用于边界节点，但这里不包括热流或热流密度和内热源之前的系数，这说明所建立的离散方程是正确的。

通过该分析方法便可得出其他三种节点，此处不再赘述。

6) 雅克比迭代法

将传统意义上的雅克比迭代通过改写约束条件，使迭代精度保持在0.005之内，使迭代结果中的平均温度与假设给定出水口温度情况下的计算平均温度二者误差在0.001之内，进而反证假设出水口温度时计算平均温度的合理性。

该算法可以得出室内各节点的温度的数值解以及室内平均温度和房间温度分布云图。

3. 模型的结果及理论分析

(一) 当室外温度为0度，进水口流量为每小时0.5立方米时，依据房间的长宽比，我们把房间分成了100 × 50个小正方形，并且利用热量守恒方程：

$\frac{0.12\left(t_n-t_{wn}\right)}{G}=40-T_{出}$

节点离散方程条件下的温度呈现点状发散分布，同时进水口温度为40℃，然后假设出水口温度为35℃，作为边界条件进行雅克比迭代演算(详细程序见《附录一》)，得到当室外温度为0℃，进水口流量为每小时0.5立方米时的演算结果如下。

结果一：ans = 20.7920；

结果二：Num = 3780；

结果三：

房间节点温度的数值解。(房间左半部分同右半部分结果呈对称分布，左半部分详细结果见《附录三》)；

结果四：

房间内的空间温度分布云图(见下图4和图5)。

Figure 4. West and south

图4. 西面和南面

Figure 5. East

图5. 东面

结果分析：

1) 通过热量守恒方程计算得到的室内平均温度为约等于20.8℃；而通过雅克比迭代计算出的室内平均温度为20.7920℃，由此证得热量守恒方程计算的室内平均温度的合理性。

2) 程序迭代次数为3780次。迭代次数越高，各节点温度分布精确度越高，但保真性较差，因而设置迭代次数为5000次以下比较合理。

3) 得到了房间各节点的温度分布的数值解，由于数据太大，我们每隔一列截取一列的节点温度，做成表格以供参考。左半部分房间节点温度见《附录三》。得到结果是靠近暖气片的温度最高，随着节点与暖气片距离的增加，温度逐渐降低。次结果仅适用于室内空气未强制对流情况下的结果。若室内有人员活动，各点平均温度将趋于室内的平均温度值。

4) 得到室内温度分布云图，直观的表达了室内的温度分布情况，可以看出温度分布与距离是呈指数函数关系递减。

(二) 当室外温度为0℃，进水口流量在0.1~1之间变化时，以0.1 m³/h的差值建立水流量间隔分配表(表1)，将每个给定的流量分别代入热量守恒方程中，计算出满足热量守恒得室内平均温度与出水口温度。再使用MATLAB进行雅克比迭代算法，基于问题的雅克比算法进行条件程序的改写，使得最终迭代的室内平均温度和热量守恒方程得出的室内平均温度在非常小的误差范围之内，这就保证了热量守恒方程所计算的出水口温度和室内平均温度的合理性。利用plot函数绘制房间重要区域温度分布集中曲线图来分析研究。

Table 1. Water flow interval distribution table

表1. 水流量间隔分配表

结果一：

分别对出水口温度和室内温度进行计算，并给出相应的结果与理由，并且对变流量下房间平均温度和出水口温度的变化情况建立了以下表2和折线图6。

Table 2. Indoor temperature verification

表2. 室内温度验证

结果二：

十种流量情况下房间内节点温度分布的数值解(部分)。

Figure 6. Dynamic average temperature of variable flow room

图6. 变流量房间动态平均温度

结果三：

并且对流量在0.1~1之间变化时，房间东墙、西墙以及中间部分建立了以下三个温度分布曲线图(见图7)。

(a) 东墙 (b)中间(c) 西墙

Figure 7. Temperature distribution curves under different flow rates

图7. 不同流量下温度分布曲线图

4. 结论

由不同流量下温度分布数值解汇总表以及房间重要区域温度分布集中曲线图可知：在不同流量下房间内的温度分布趋势基本一致，且只与出水口温度相关。

在同一流量下，出水口温度越高，房间内的平均温度偏高；出水口温度越低，房间的平均温度越低。各节点的温度变化趋势与平均温度变化基本一致。

在不同流量下，出水口温度在一个稳定区域内变化时，若流量越大，则房间平均温度越高；流量越小，房间平均温度越低。各节点的温度变化趋势与平均温度变化基本一致。

由此验证：所假设的热量守恒方程中，流量、出水口温度和室内温度存在一定的关系和相关性。出水口温度和室内平均温度是负相关的；室内平均温度与流量是正相关的。即验证了热量守恒方程假设的合理性。

所构建的室内温度分布模型具有通用性。此模型可以根据温度源的温度，计算任何封闭空间内各个点的温度大小，并且能确定封闭空间下的平均温度，并且做出房间温度分布云图直观地了解温度分布。可以被类似于植物温室和植物工厂等多种需要考虑房间温度的封闭空间利用。将该模型及模拟算法与控制理论相结合，便可控制房间温度，维持稳态。

附录一

MATLAB雅克比迭代程序代码

t=zeros(101,51);

t0=zeros(101,51);

e=0.005;

Num=0;

for j=16:36

t(1,j)=40;

end%上边界40度节点

for j=16:36

t(101,j)=35;

end%下边界35度节点

t0=t;

for k=1:4000%迭代次数上限为4000

for j=2:15

t(1,j)=(t0(1,j+1)+t0(1,j-1)+2*t0(2,j))/4;

end%上左边界剩余节点

for j=37:50

t(1,j)=(t0(1,j+1)+t0(1,j-1)+2*t0(2,j))/4;

end%上右边界剩余节点

for i=2:100

for j=2:50

t(i,j)=(t0(i-1,j)+t0(i+1,j)+t0(i,j-1)+t0(i,j+1))/4;

end

end%内部节点

for j=2:15

t(101,j)=(t0(101,j+1)+t0(101,j-1)+2*t0(100,j))/4;

end%下左边界剩余节点

for j=37:50

t(101,j)=(t0(101,j+1)+t0(101,j-1)+2*t0(100,j))/4;

end%下右边界剩余节点

for i=2:100

t(i,1)=(t0(i-1,1)+t0(i+1,1)+2*t0(i,2))/4;

end%左边界节点

for i=2:100

t(i,51)=(t0(i-1,51)+t0(i+1,51)+2*t0(i,50))/4;

end%右边界节点

t(1,1)=(t0(2,1)+t0(1,2)+2*t0(2,2))/4;

t(1,51)=(t0(1,50)+t0(2,51)+2*t0(2,50))/4;

t(101,1)=(t0(100,1)+t0(101,2)+2*t0(100,2))/4;

t(101,51)=(t0(100,51)+t0(101,50)+2*t0(100,50))/4;

%四个边界角节点

dtmax=0;

for i=1:101

for j=1:51

dtmax=max(abs(t(i,j)-t0(i,j)),dtmax);

end

end

t0=t;

Num=Num+1;

if dtmax

break;

%使两次迭代结果之差满足标准误差

end

end

subplot(1,2,1);

surf(t);%绘制图1

subplot(1,2,2);

contour(t',1000);%绘制图二

mean(mean(t))%求迭代算法得出的数值解的平均值

t%记录各点温度数值解

Num%计算最终迭代次数

附录二

水流量间隔分配表

附录三

左半部分房间节点温度

房间内坐标(x,y)；x = 0.2i、y = 0.1j (单位：m)