[数据结构与算法]用Matlab求解一维非稳态导热问题（显式离散）

章熙民的第六版《传热学》里，较为简单的介绍了非稳态导热的数值计算，本文根据此书，以计算一个可视为无限大平壁的复合墙体传热过程为例，讨论一维非稳态导热问题数值求解的问题。

初始条件：一个由四层材料构成的墙体（可视为无限大平壁），已知材料热工性质如下，该墙体的内表面换热系数为8.7W/(m^2*K)，外表面换热系数为23W/(m^2*K)，室内温度为20℃恒定不变，室外温度随昼夜变化。

现给出一星期的室外温度数据（每小时），求墙体内部温度场。

• 前期准备

在《传热学》第四章第三节中，对什么是非稳态导热的给出了如下解释：

因此根据上述场景将该问题提炼为：一个两边是第三类边界条件的无限大平壁的一维非稳态导热问题。这里我们选用显示离散格式（求解上更为简单，但是对步长要求较为严苛）。

内部节点方程式：

稳定条件：

边界节点方程式（第三类边界条件）：

?稳定条件：

• ?问题分析

在经过多次试验之后，将该墙体按间距3mm划分为132层，取时间步长为1s，此时满足稳定条件。由于这里是复合墙体，因此后续计算时也要注意节点的热物性参数会相应改变。

还有一个重要的问题是墙体的初始温度场，显式离散格式的逻辑就是从初始温度出发逐个计算出各个时刻的温度。这里没有给，所以就用初始温度条件（室内20℃，室外-25.4℃）计算这一条件下的稳态温度。

?根据上式求得墙体内初始温度场。

• ?Matlab求解

在完成上述工作之后，便可以在Matlab中求解了。

clc;
clf;
format short e;
global N delta_x delta_t TN;

h1 = 8.7; %墙体内表面传热系数 W/(m^2*K)
h4 = 23; %墙体外表面传热系数 W/(m^2*K)
D = 0.332; %墙体厚度/m
N = 124; %厚度划分数目
delta_x = 0.003; %划分间距/m
delta_t = 1; %时间步长/s
day = 7; %计算天数
TN = 3600*24*day; %计算时间节点数目
Tn = 20; %室内温度
T = zeros(TN,N+1); %预设数组大小
Tf = zeros(1,TN);

%设置室外空气温度 这一步是因为给的是每小时的温度 换成每秒均等增加的温度值
for hour = 1:24*day %取值范围
    %将第hour小时和hour+1小时的温度等间距分成3600个
    %自定义数组outdoor_air为室外温度
    temp = linspace(outdoor_air(hour),outdoor_air(hour+1),3600); 
    for t = 1:3600
        Tf(3600*(hour-1)+t) = temp(t);%每秒均匀变化的温度
    end
end

%设置初始墙体内温度
for I = 1:N+1
    %自定义数组wall为墙体初始温度
    T(1,I) = wall(I);
end

%第一层：石膏板
a1 = 3.265e-7; %热扩散率 m^2/s
lamb1 = 0.36; %热导率 W/(m*K)
Fo1 = a1*delta_t/(delta_x)^2;
Bi1 = h1*delta_x/lamb1;
c1 = Fo1 - 1/(2*Bi1+2);

%第二层：钢筋混凝土
a2 = 8.870e-7; %热扩散率 m^2/s
lamb2 = 2.04; %热导率 W/(m*K)
Fo2 = a2*delta_t/(delta_x)^2;

%第三层：聚苯乙烯泡沫板
a3 = 1.667e-6; %热扩散率 m^2/s
lamb3 = 0.046; %热导率 W/(m*K)
Fo3 = a3*delta_t/(delta_x)^2;

%第四层：陶瓷红砖
a4 = 4.419e-7; %热扩散率 m^2/s
lamb4 = 0.7; %热导率 W/(m*K)
Fo4 = a4*delta_t/(delta_x)^2;
Bi4 = h4*delta_x/lamb4;
c4 = Fo4 - 1/(2*Bi4+2);

while (Fo1>0.5 || Fo2>0.5 || Fo3>0.5 || Fo4>0.5 || c4>0  || c1>0)
    error('无法收敛');
end

?第一部分主要是输入参数，设置初始条件，判断选择的步长是否满足稳定性条件。

for t = 1:TN
    for i = 1:N+1      
        if i == 1
            T(t+1,i) = 2*Fo1*(T(t,i+1) + Bi1*Tn) + (1 - 2*Bi1*Fo1 - 2*Fo1)*T(t,i);
        end
        
        if (i>=2 && i<= 5)
            T(t+1,i) =  Fo1*(T(t,i-1) + T(t,i+1)) + (1 - 2*Fo1)*T(t,i);
        end
        
        if (i>=6 && i<= 45)
            T(t+1,i) =  Fo2*(T(t,i-1) + T(t,i+1)) + (1 - 2*Fo2)*T(t,i);
        end
        
        if (i>=46 && i<= 85)
            T(t+1,i) =  Fo3*(T(t,i-1) + T(t,i+1)) + (1 - 2*Fo3)*T(t,i);
        end
        
        if (i>=86 && i<= 124)
            T(t+1,i) =  Fo4*(T(t,i-1) + T(t,i+1)) + (1 - 2*Fo4)*T(t,i);
        end
        
        if i == N+1
            T(t+1,i) = 2*Fo4*(T(t,i-1) + Bi4*Tf(t)) + (1 - 2*Bi4*Fo4 - 2*Fo4)*T(t,i);
            %这里给了室外空气温度Tf一个时间函数关系
        end
    end
end

?第二部分就是计算显式离散节点方程了，没有什么好说的，就用了最简单的迭代法，甚至不需要求解方程组，唯一需要注意的就是各个材料层的节点方程式使用的参数不一样。

X = zeros(1,24*day*2);
Y = zeros(1,N+1);
TT = zeros(24*day*2,N+1);

J = 1;
for I = 1:24*day*2
    X(I) = 1800*I; %30分钟一个数
    for J = 1:N+1
        Y(J) = J;
        TT(I,J) = T(X(I),Y(J));
    end
end
 
figure(1)
hold on;
 
% Make a contour with the presence of isolines & texts
[XX,YY] = meshgrid(Y,X);
surfc(YY,XX,TT)

xlabel('Time - AXIS')
ylabel('Thickness - AXIS')
zlabel('T (^{\circ}C)')
title(colorbar,'^{\circ}C')

?第三部分画图，半个小时取一个数，因为一秒钟一个温度数太多了没有必要，画图也不好看。

最后呈现出来的温度-时间-节点的图形就是这样：