clc;
clear;
%分为五层，假人为一层，各项参数近似与真人相同
Lambda = [0.082 0.37 0.045 0.028 0.432];%热导率
Ci = [1377 2100 1726 1005 4.2];%比热容
Rho = [300 862 74.2 1.18 1060];%密度
T = 3600;%时间分为T个单位
deltat = 3600 / T;%每个单位时间
deltax = 0.1 * 1e-3;%取衣物及人体厚度0.01mm为单位长度求解
Tmax = 65.00;%最外层温度
Tmin = 37.00;%最内层温度
n = 61;
%第一问算出假人有效吸热厚度为61个单位长度
Hb = 0;
Ub = zeros(T);
for hi = 25 : -0.1 : 0.6
H = [0.6 hi 3.6 5.5 n * deltax * 1e+3] * 1e-3;
%各层厚度
Hsum = zeros(1,6);
for i = 1 : 5
Hsum(i + 1) = Hsum(i) + H(i);
end
L = Hsum(6);
M = int32(L / deltax);%总单位长度数
m = M - n;
x = linspace(0,L,M);%将厚度分为M块
Lambdax = zeros(1,M);%得到各个单位长度lambda
Cix = zeros(1,M);%得到各个单位长度Ci
Rhox = zeros(1,M);%得到各个单位长度Rho
for i = 1 : M
for j = 1 : 5
if x(i) <= Hsum(j + 1)
Lambdax(i) = Lambda(j);
Cix(i) = Ci(j);
Rhox(i) = Rho(j);
break;
end
end
end
alphax = [0 Lambdax * deltat ./ Cix ./ Rhox / 2 / deltax^2 0];
U = zeros(T,M + 2);%温度分布
U(1,:) = zeros(1,M + 2) + Tmin;%初始化温度分布
U(1,1) = Tmax;
for i = 2 : T%微分方程求解
MatrixU = zeros(M + 2);
D_last = zeros(M + 2,1);
MatrixU(1,1) = 1 + 2 * alphax(1);MatrixU(1,2) = -alphax(1);
D_last(1) = Tmax;
for row = 2 : M + 1
MatrixU(row,row - 1) = -alphax(row);
MatrixU(row,row) = 2 * alphax(row);
MatrixU(row,row + 1) = -alphax(row);
D_last(row) = alphax(row) * U(i - 1,row + 1) + alphax(row) * U(i - 1,row - 1) + (1 - 2 * alphax(row)) * U(i - 1,row);
end
MatrixU(M + 2,M + 1) = -alphax(M + 2);
MatrixU(M + 2,M + 2) = 1 + 2 * alphax(M + 2);
D_last(M + 2) = Tmin;
U(i,:) = reshape(MatrixU \ D_last,1,M + 2);
for row = 1 : M + 2
if U(i,row) > 65.00
U(i,row) = 65.00;
end
end
end
temp = 0;
tcount = 0;
for i = 1 : T
if U(i,m + 1) > 47.00 || tcount > 5 * 60
temp = 1;
break;
end
if U(i,m + 1) > 44.00
tcount = tcount + 1;
end
end
if temp == 0
Hb = hi;
Ub = U(:,m + 1);
end
end
Hb %8.9
figure(1);
axis([0 5500 35 55]);
hold on;
plot(1:T,Ub);%最优解为Ub(m+1)
