首页/文章/ 详情

二维热导方程Matlab数值解案例

过冷水

3年前浏览3157

过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，群号：927550334

过冷水之前有和大家分享一维热传导方程的有限差分法求解的具体实现过程。本次就和大家一起来看看如何根据一维热传导有限差分法的思想求解二维热传导方程进行求解。形式如下的二维热传导方程：

边界条件是：

初值为：

差微分方法思路：

变形为迭代式如下：

根据该迭代过程就可求出任意时刻的温度T随空间的分布情况。现在回到我们之前和大家分享的二维微分方程的具体案例中看一下怎么解：

clear all
a=0.5;
u_xy0=inline('0','x','y');
u_xyt=inline('x^2*cos(pi*y)-y^2*cos(pi*x)','x','y','t');
f=inline('0','x','y','t');
D=[0,2*pi,0,2*pi];
T=500;
Mx=100;
My=100;
N=100;
ox=(D(2)-D(1))/Mx;
x=D(1) [0:Mx]*ox;
oy=(D(4)-D(3))/My;
y=D(3) [0:My]*oy;
ot=T/N;
t=[0:N]*ot;
%初始化 u
for j=1:Mx 1;
    for i=1:My 1;
        u(i,j)=u_xy0(x(j),y(i));
    end
end
rx=a*ot/(ox*ox);
ry=a*ot/(oy*oy);
rx1=1 2*rx;
rx2=1-2*rx;
ry1=1 2*ry;
ry2=1-2*ry;
for j=1:Mx-1;
    A(j,j)=ry1;
    if j>1
        A(j-1,j)=-ry;
        A(j,j-1)=-ry;
    end
end
%A为y方向隐式时的系数矩阵
for i=1:My-1;
    B(i,i)=rx1;
    if i>1;
        B(i-1,i)=-rx;
        B(i,i-1)=-rx;
    end
end
%B为x方向隐式时的系数矩阵
for k=1:N
    t=k*ot;u_1=u;
    %for m=1:Mx 1;
        %for n=1:My 1;
            %v(n,m)=feval(f,x(m),y(n),t);
        %end
    %end
    u_1=u;
    for i=1:My 1
        u(i,1)=feval(u_xyt,x(1),y(i),t);
        u(i,Mx 1)=feval(u_xyt,x(Mx 1),y(i),t);
    end
    for j=1:My 1;
        u(1,j)=feval(u_xyt,x(j),y(1),t);
        u(My 1,j)=feval(u_xyt,x(j),y(My 1),t);
    end
    if mod(k,2)==0;
        for i=2:My;
            jj=2:Mx;
            bx=[ry*u(i,1),zeros(1,Mx-3),ry*u(i,My 1)] rx*(u_1(i-1,jj) u_1(i 1,jj)) rx2*u_1(i,jj) ot*u_1(i,jj);
            opts.UT = true; opts.TRANSA = true;
            u(i,jj)=linsolve(A,bx',opts)';
        end
    else
        for j=2:Mx
            ii=2:My;
            by=[rx*u(1,j);zeros(My-3,1);rx*u(Mx 1,j)] ry*(u_1(ii,j-1) u_1(ii,j-1)) ry2*u_1(ii,j) ot*u_1(ii,j);
            opts.UT = true; opts.TRANSA = true;
            u(ii,j)=linsolve(B,by,opts);
        end
    end
end
mesh(x,y,u);
xlabel('x');
zlabel('温度T')

过冷水本期要和大家分享的二维热传导的问题就是这些了，热传导属于比较复杂的问题，物理模型不容易理解，解决了物理模型其实就是偏微分方程的求解问题，