%%详细推导求解过程请移步 微信公众号 CFD入门指南 %%
%%二维稳态导热问题第三类边界%%
clc
clear
L1=0.6; %板长
L2=0.6; %板宽
T1=80; %边界温度
tf=20;
lan=50; %导热系数
m=101; %长度的节点数量
n=101; %宽度的节点数量
dx=L1/(m-1); %长度的网格尺寸
dy=L2/(n-1); %宽度的网格尺寸
h1=50;
h2=100;
h3=100;
Bi1=h1*dx/lan;
Bi2=h2*dx/lan;
Bi3=h3*dx/lan;
v1=zeros(m,n);
for i=1:m %温度边界条件
v1(i,n)=T1;
end
dt=1;
v2=v1;
k=0;
while dt>1e-3 %两次迭代误差小于10^-3时计算结束
k=k+1;
dt=0;
for i=2:m-1
for j=2:n-1
v2(i,1)=1/(2*Bi1+4)*(v1(i+1,1)+v1(i-1,1)+2*v1(i,2)+2*Bi1*tf); %对流换热边界
v2(1,j)=1/(2*Bi2+4)*(v1(1,j+1)+v1(1,j-1)+2*v1(2,j)+2*Bi2*tf);
v2(m,j)=1/(2*Bi3+4)*(v1(m,j+1)+v1(m,j-1)+2*v1(m-1,j)+2*Bi3*tf);
v2(1,1)=1/(2*(Bi1+Bi2)+2)*(v1(2,1)+v1(1,2)+(Bi1+Bi2)*tf);
v2(m,1)=1/(2*(Bi1+Bi3)+2)*(v1(m-1,1)+v1(m,2)+(Bi1+Bi3)*tf);
v2(i,j)=(v1(i,j+1)+v1(i,j-1))*(dx^2/2/(dx^2+dy^2))+(v1(i+1,j)+v1(i-1,j))*(dy^2/2/(dx^2+dy^2));
dt=max((v2(i,j)-v1(i,j)),dt);
end
end
v1=v2;
end
%%详细推导求解过程请移步 微信公众号 CFD入门指南 %%
figure
surf(v2);
a=xlabel('长度');
b=ylabel('宽度');
c=zlabel('温度');
set(a,'FontSize',15');
set(b,'FontSize',15');
set(c,'FontSize',15');
k
%%详细推导求解过程请移步 微信公众号 CFD入门指南 %%
