matlab温度场程序

clear
echo on
d1=43;d2=7;dx=0.15;dy=0.1;xy=dx/dy;yx=dy/dx;
t=zeros(d1,d2);t1=ones(d1,d2);t0=zeros(d1,d2);
x=zeros(d1);y=zeros(d2);
x(1)=0;x(2)=dx/2;
for i=3:d1-1;
x(i)=x(i-1)+dx;
end
x(d1)=(d1-2)*dx;
y(1)=0;y(2)=dy/2;
for i=3:d2-1;
y(i)=y(i-1)+dy;
end
y(d2)=(d2-2)*dy;
t1=20*ones(d1,d2);
t=zeros(d1);
dt=0.1;ttt=30;
%nnn=tt/dt;
echo off
%for iii=1:nnn
% ttt=iii*dt;
tf=30;af=6.6;af=1/af;bta=-6.6;v=0.0625;
tin=100;tout=100;d=0.05;l=6;
for i=1:42;
if x(i)t1(i,1)=tin+(300-tout)*x(i)/(v*ttt);
elseif x(i)>v*ttt+28*d
t1(i,1)=300-(300-tout)*(x(i)-v*ttt-28*d)/(l-28*d-v*ttt);
else
zz=-0.123*(x(i)-v*ttt)/d-3.52*exp(-0.123*(x(i)-v*ttt)/d);
t1(i,1)=10060*exp(zz);
end
end
for i=1:41;
for j=2:6
t1(i,j)=t1(i,1)-10*(j-1);
end
end
for iii=1:500;
t0=t1;
cd=1300;a=0.0003;
an=[1.69,-0.594,0.401,-0.168,0.027,-0.037,0.046,-0.05, 0.039,-0.012];
bn=[0, 0.333, 0.017,-0.131,0.054,0.0003,0.007,-0.012,0.026,0];
fncp=zeros(d1,d2);
for i=1:10;
fncp=an(i)*cos(i*pi/9*t1)+bn(i)*sin(i*pi/9*t1)+fncp;
end
fncp=fncp/4.1868;
b=1.3e-6;c=1.5e-9;
fnk=(a+b*t1+c*t1.^2)*360;
for i=2:42;
for j=2:6;
fnae(i,j)=2*yx*fnk(i,j)*fnk(i+1,j)/(fnk(i,j)+fnk(i+1,j));
fnaw(i,j)=2*yx*fnk(i,j)*fnk(i-1,j)/(fnk(i,j)+fnk(i-1,j));
fnan(i,j)=2*xy*fnk(i,j)*fnk(i,j+1)/(fnk(i,j)+fnk(i,j+1));
fnas(i,j)=2*xy*fnk(i,j)*fnk(i,j-1)/(fnk(i,j)+fnk(i,j-1));
end
end
fnap0=cd*fncp*dx*dy/dt;
fnbb=fnap0.*t1;
%for i=2:41
% for j=1:5
% t1(i,j)=t(i,1)-10*(j-1);
% end
% end
% t0=t1;
kk=af+0.5*dx/fnk(1,2);
bb=fnbb(2,2)+tf*dy/kk;
ap=fnae(2,2)+fnan(2,2)+fnap0(2,2)+dy/kk+2*xy*fnk(2,1)-bta*dx*dy;
fff=fnae(2,2)*t1(3,2)+fnan(2,2)*t1(2,3)+2*xy*fnk(2,1)*t1(2,1)+bb;
t1(2,2)=fff/ap;
for j=3:5;
kk=af+0.5*dx/fnk(1,j);
bb=fnbb(2,j)+tf*dy/kk;
ap=fnae(2,j)+fnas(2,j)+fnan(2,j)+fnap0(2,j)+dy/kk-bta*dx*dy;
fff=fnae(2,j)*t1(3,j)+fnan(2,j)*t1(2,j+1)+fnas(2,j)*t1(2,j-1)+bb;
t1(2,j)=fff/ap;
end
kk1=af+0.5*dx/fnk(1,6);
kk2=af+0.5*dy/fnk(1,6);
bb=fnbb(2,6)+tf*dy/kk1+tf*dx/kk2;
ap=fnae(2,6)+fnas(2,6)+fnap0(2,6)+dy/kk1+dx/kk2-bta*dx*dy;
fff=fnae(2,6)*t1(3,6)+fnas(2,6)*t1(2,5)+bb;
t1(2,6)=fff/ap;
for i=3:40;
for j=2:6;
if j==2;
as=2*xy*fnk(i,1);
fff=fnae(i,2)*t1(i+1,2)+fnaw(i,2)*t1(i-1,2)+fnan(i,2)*t1(i,3)+as*t1(i,1)+fnbb(i,2);
ap=fnae(i,2)+fnaw(i,2)+fnan(i,2)+as+fnap0(i,2)+dy/kk-bta*dx*dy;
t1(i,2)=fff/ap;
elseif j==6
kk=af+0.5*dy/fnk(i,6);
bb=fnbb(i,6)+tf*dx/kk;
fff=fnae(i,6)*t1(i+1,6)+fnaw(i,6)*t1(i-1,6)+fnas(i,6)*t1(i

,4)+bb;
ap=fnae(i,6)+fnaw(i,6)+fnas(i,6)+fnap0(i,6)+dx/kk-bta*dx*dy;
t1(i,5)=fff/ap;
else
fff=fnae(i,j)*t1(i+1,j)+fnaw(i,j)*t1(i-1,j)+fnan(i,j)*t1(i,j+1)+fnas(i,j)*t1(i,j-1)+fnbb(i,j);
ap=fnae(i,j)+fnaw(i,j)+fnan(i,j)+fnas(i,j)+fnap0(i,j)-bta*dx*dy;
t1(i,j)=fff/ap;
end
end
end
for j=3:5;
kk=af+0.5*dx/fnk(42,j);
bb=fnbb(41,j)+tf*dy/kk;
ap=dy/kk+fnaw(41,j)+fnan(41,j)+fnas(41,j)+fnap0(41,j)-bta*dx*dy;
fff=fnaw(41,j)*t1(40,j)+fnan(41,j)*t1(41,j+1)+fnas(41,j)*t1(41,j-1)+bb;
t1(41,j)=fff/ap;
end
kk1=af+0.5*dx/fnk(42,6);
kk2=af+0.5*dy/fnk(41,7);
bb=fnbb(41,6)+tf*dy/kk1+tf*dx/kk2;
ap=fnaw(41,6)+fnas(41,6)+fnap0(41,6)+dy/kk1+dx/kk2-bta*dx*dy;
fff=fnaw(41,6)*t1(40,6)+fnas(41,6)*t1(41,5)+bb;
t1(41,6)=fff/ap;
kk=af+0.5*dx/fnk(42,2);
bb=fnbb(41,2)+tf*dy/kk;
ap=fnaw(41,2)+fnan(41,2)+fnap0(41,2)+dy/kk+2*xy*fnk(41,1)-bta*dx*dy;
fff=fnaw(41,2)*t1(40,2)+fnan(41,2)*t1(41,3)+2*xy*fnk(41,1)*t1(41,1)+bb;
t1(41,2)=fff/ap;

if abs(abs(t1)-abs(t0))<20;
break;
%else
%t1=t0+0.5*(t1-t0);
end
end
%if rem(iii,10)==0
figure
subplot(3,2,1)
t111=t1';
%mesh(t111);
%view(2);
pcolor(t111);
subplot(3,2,5)
ccc=contour(t111,[200,400,600,800,1000]);
%ccc=contour(t111,5);
clabel(ccc);
%end
% 将平面网格映射到三维网格面上
x1=1:7;
y1=1:43;
[X1,Y1]=meshgrid(y1,x1);
%绘制加密前三维网格
subplot(3,2,3);
mesh(X1,Y1,t111);
% 用二元差值公式将网格加密并绘出温度场的三维图
x2=1:.25:7;
y2=1:.25:43;
[X2,Y2]=meshgrid(y2,x2);
Z2=griddata(X1,Y1,t111,X2,Y2);
t111=round(t111);
subplot(3,2,4);
mesh(X2,Y2,Z2)
% 绘制加密后的平面网格
subplot(3,2,2);
pcolor(Z2)
% 绘制加密后的温度场等高线图
subplot(3,2,6)
CCC=contour(Z2,[100,200,300,400,500,600,700,800,900,1000]);
clabel(CCC);
% end
%下面输入平板中某点坐标以查看此点温度
a1=input('请输入你要查看温度的点的X轴坐标值（范围为1到169）');
b1=input('请输入你要查看温度的点的Y轴坐标值（范围为1到25）');
disp('输入坐标值为'),a1,b1
a2=round((a1+3)/4);
b2=round((b1+3)/4);
disp('计算出网格加密前坐标点的温度为')
Q1=(t111(b2,a2))
display('计算出网格加密后坐标点的温度为')
Q2= Z2(b1,a1)