有限差分法求解二维热传导方程

ziguayida

方程：

边界条件：

我想用有限差分法进行求解，方程的离散过程以及编程如下：

编程实现：

NX=100;

NY=40;

NT=73;

ox=50/100;

ot=1;

P0=75000;

Pa=101325;

A1=(1-0.873/800)*220*0.15*0.1386;

A2=(1-0.873/800)*220*0.55*0.0231;

A3=(1-0.873/800)*220*0.3*0.017328;

%初始条件

for j=1:NX+1

    for l=1:NY+1

        P(j,l,1)=P0;

        P(j,l,2)=P0;

    end

end

%顶部边界条件

for m=2:NT

    for j=2:NX-1

        P(j,1,m)=Pa;

    end

end

%左边边界条件

for m=2:NT

    for l=1:NY

        P(1,l,m)=-2000;

    end

end

%底部边界条件

for m=2:NT

    for j=2:NX-1

        P(j,NY+1,m)=P(j,NY-1,m);

    end

end

%右边边界条件

for m=2:NT

    for l=2:NY-1

       P(NX+1,l,m)=P(NX-1,l,m);

    end

end

%内部网格

for m=2:NT

    for i=1:10

        for j=2:NX-1

            for l=2:NY-1

                t(l,m-1)=(m-1)*ot+l*ox/2;

                B(l,m-1)=0.873*1100*9.8*l*ox*0.365*10^(-6)*exp(-10^(-3)*(m-1)*ot*365);

                f(l,m-1)=1.145*(A1*exp(-0.1386*t(l,m-1))+A2*exp(-0.0231*t(l,m-1))+A3*exp(-0.017328*t(l,m-1))+B(l,m-1));

                P(j,l,m)=4.3*365*24*3600*((3*P(j+1,l,m)+3*P(j-1,l,m)+P(j,l+1,m)+P(j,l-1,m)-8*P(j,l,m))/(ox*ox))*ot+f(l,m-1)*ot+P(j,l,m-1)

            end

        end

            for j=2:NX-1

                t(NY,m-1)=(m-1)*ot+NY*ox/2;

                B(NY,m-1)=0.873*1100*9.8*NY*ox*0.365*10^(-6)*exp(-10^(-3)*(m-1)*ot*365);

                f(NY,m-1)=1.145*(A1*exp(-0.1386*t(NY,m-1))+A2*exp(-0.0231*t(NY,m-1))+A3*exp(-0.017328*t(NY,m-1))+B(NY,m-1));        

                P(j,NY,m)=4.3*365*24*3600*((3*P(j+1,NY,m)+3*P(j-1,NY,m)+P(j,NY+1,m)+P(j,NY-1,m)-8*P(j,l,m))/(ox*ox))*ot+f(NY,m-1)*ot+P(j,NY,m-1);

            end

            for l=2:NY-1

                t(l,m-1)=(m-1)*ot+l*ox/2;

                B(l,m-1)=0.873*1100*9.8*l*ox*0.365*10^(-6)*exp(-10^(-3)*(m-1)*ot*365);

                f(l,m-1)=1.145*(A1*exp(-0.1386*t(l,m-1))+A2*exp(-0.0231*t(l,m-1))+A3*exp(-0.017328*t(l,m-1))+B(NY,m-1));

                P(NX,l,m)=4.3*365*24*3600*((3*P(NX+1,l,m)+3*P(NX-1,l,m)+P(NX,l+1,m)+P(NX,l-1,m)-8*P(j,l,m))/(ox*ox))*ot+f(l,m-1)*ot+P(j,l,m-1);

            end

    end

end

根据我自己编的程序在matlab中计算，结果十分不理想，不知道怎么回事，跪求各位大神指导！

ziguayida

跪求大神们的帮助啊，小弟在此赶紧不尽啊！！！

ziguayida

跪等各位大神的援手啊！不要沉呀！

ziguayida

求好心人帮忙:'(

ziguayida

:Q:'(跪等大神的帮助与指点啊！！！

徐全

也想看看哪位牛人可以解决啊:o