matlab温度场程序

clearechoond1=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;fori=3:d1-1;x(i)=x(i-1)+dx;endx(d1)=(d1-2)*dx;y(1)=0;y(2)=dy/2;fori=3:d2-1;y(i)=y(i-1)+dy;endy(d2)=(d2-2)*dy;t1=20*ones(d1,d2);t=zeros(d1);dt=0.1;ttt=30;%nnn=tt/dt;echooff%foriii=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;fori=1:42;ifx(i)v*tttt1(i,1)=tin+(300-tout)*x(i)/(v*ttt);elseifx(i)v*ttt+28*dt1(i,1)=300-(300-tout)*(x(i)-v*ttt-28*d)/(l-28*d-v*ttt);elsezz=-0.123*(x(i)-v*ttt)/d-3.52*exp(-0.123*(x(i)-v*ttt)/d);t1(i,1)=10060*exp(zz);endendfori=1:41;forj=2:6t1(i,j)=t1(i,1)-10*(j-1);endendforiii=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);fori=1:10;fncp=an(i)*cos(i*pi/9*t1)+bn(i)*sin(i*pi/9*t1)+fncp;endfncp=fncp/4.1868;b=1.3e-6;c=1.5e-9;fnk=(a+b*t1+c*t1.^2)*360;fori=2:42;forj=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));endendfnap0=cd*fncp*dx*dy/dt;fnbb=fnap0.*t1;%fori=2:41%forj=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;forj=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;endkk1=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;fori=3:40;forj=2:6;ifj==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;elseifj==6kk=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;elsefff=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;endendendforj=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;endkk1=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;ifabs(abs(t1)-abs(t0))20;break;%else%t1=t0+0.5*(t1-t0);endend%ifrem(iii,10)==0figuresubplot(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,b1a2=round((a1+3)/4);b2=round((b1+3)/4);disp('计算出网格加密前坐标点的温度为')Q1=(t111(b2,a2))display('计算出网格加密后坐标点的温度为')Q2=Z2(b1,a1)