syms i j n ;
t0=18;
tf=44;
H=1;
w=30; %导热系数
h1=10;h2=25;h3=15; %对流换热系数
n=20;%节点数
detax=H/(n-1);%步长
t=t0*ones(n,n);%设一个初始温度矩阵
t(:,1)=t0;
k=1; %迭代次数
h=t;
tic %迭代初始时间
while k<1e99
for i=2:1:n-1
for j=2:1:n-1
t(:,1)=t0; %左边界
t(i,j)=0.25*(t(i-1,j)+t(i+1,j)+t(i,j-1)+t(i,j+1)); %内部节点
end
end
for j=2:1:n-1
t(j,n)=1/(2*(h2*detax/w)+4)*(2*t(j,n-1)+t(j+1,n)+t(j-1,n)+2*h2*detax/w*tf); %右边界
end
for i=2:1:n-1
t(n,i)=1/(2*(h1*detax/w)+4)*(2*t(n-1,i)+t(n,i+1)+t(n,i-1)+2*h1*detax/w*tf);%下边界
t(1,i)=1/(2*(h3*detax/w)+4)*(2*t(2,i)+t(1,i-1)+t(1,i+1)+2*h3*detax/w*tf);%上边界
end
t(n,n)=1/((h1+h2)*detax/w+2)*(t(n,n-1)+t(n-1,n)+(h2+h1)*detax/w*tf); %C点
t(1,n)=1/((h3+h2)*detax/w+2)*(t(1,n-1)+t(2,n)+(h2+h3)*detax/w*tf); %D点
if max(max(t-h))<1e-6
break
else
k=k+1;
h=t;
end
end
z=toc; %迭代结束时间
%%%%%%——输出温度云图——%%%%%%
xx=linspace(0,n,n);
yy=linspace(0,n,n);
[x,y]=meshgrid(xx,yy);
surf(x,y,t) %输出温度图
%%%%%%——边界热流——%%%%%%
q0=(sum(t(2:n-1,1))-sum(t(2:n-1,2)))*w+(t(1,1)-t(1,2)+t(n,1)-t(n,2))*w*0.5; %左边界
q1=((n-2)*tf-sum(t(n,2:n-1)))*h1*detax+(2*tf-t(n,1)-t(n,n))*h1*0.5*detax; %下边界
q2=((n-2)*tf-sum(t(2:n-1,n)))*h2*detax+(2*tf-t(1,n)-t(n,n))*h2*0.5*detax; %右边界
q3=((n-2)*tf-sum(t(1,2:n-1)))*h3*detax+(2*tf-t(1,1)-t(1,n))*h3*0.5*detax; %上边界
%%%%%%——平均温度——%%%%%%
tm=mean(t(:));
评论0