基于MATLAB的温度场模拟

clear
echo on
d1=43;d2=7;dx=0.15;dy=0.1;xy=dx/dy;yx=dy/dx;
t=zeros（d1d2）;t1=ones（d1d2）;t0=zeros（d1d2）;
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（d1d2）;
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（i1）=tin+（300-tout）*x（i）/（v*ttt）;
        elseif x（i）>v*ttt+28*d
            t1（i1）=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（i1）=10060*exp（zz）;
        end
    end
    for i=1:41
        for j=2:6
            t1（ij）=t1（i1）-10*（j-1）;
        end
    end
for iii=1:500
    t0=t1;
    cd=1300;a=0.0003;
    an=[1.69-0.5940.401-0.1680.027-0.0370.046-0.05 0.039-0.012];
    bn=[0   0.333 0.017-0.1310.0540.00030.007-0.0120.0260];
    fncp=zeros（d1d2）;
    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（ij）=2*yx*fnk（ij）*fnk（i+1j）/（fnk（ij）+fnk（i+1j））;
            fnaw（ij）=2*yx*fnk（ij）*fnk（i-1j）/（fnk（ij）+fnk（i-1j））;
            fnan（ij）=2*xy*fnk（ij）*fnk（ij+1）/（fnk（ij）+fnk（ij+1））;
            fnas（ij）=2*xy*fnk（ij）*fnk（ij-1）/（fnk（ij）+fnk（ij-1））;
        end
    end
    fnap0=cd*fncp*dx*dy/dt;
    fnbb=fnap0.*t1;
    %for i=2:41
    %    for j=1:5
    %        t1（ij）=t（i1）-10*（j-1）;
    %    end
    % end
   % t0=t1;
    kk=af+0.5*dx/fnk（12）;
    bb=fnbb（22）+tf*dy/kk;
    ap=fnae（22）+fnan（22）+fnap0（22）+dy/kk+2*xy*fnk（21）-bta*dx*dy;
    fff=fnae（22）*t1（32）+fnan（22）*t1（23）+2*xy*fnk（21）*t1（21）+bb;
    t1（22）=fff/ap;
    for j=3;5
        kk=af+0.5*dx/fnk（1j）;
        bb=fnbb（2j）+tf*dy/kk;
        ap=fnae（2j）+fnas（2j）+fnan（2j）+fnap0（2j）+dy/kk-bta*dx*dy;
        fff=fnae（2j）*t1（3j）+fnan（2j）*t1（2j+1）+fnas（2j）*t1（2j-1）+bb;
        t1（2j）=fff/ap;