How to solve 2D heat equation for a sector of a circular disc using finite difference method?

7 次查看（过去 30 天）

显示更早的评论

Vignesh Ramakrishnan 2020-9-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/589483-how-to-solve-2d-heat-equation-for-a-sector-of-a-circular-disc-using-finite-difference-method

评论： Vignesh Ramakrishnan 2020-9-29

I am starting off with Dirichlet boundaries. How do I generate a sector of a circular mesh ?

6 个评论
显示 4更早的评论隐藏 4更早的评论

Vignesh Ramakrishnan 2020-9-28

编辑：Vignesh Ramakrishnan 2020-9-28

在 MATLAB Online 中打开

I get a good mesh, but there are unnecessary peaks along theta=0 radial line.

Here's the code

clear all;
clc;
 close all;
 alpha=0.001;
%grid size
 Nr = 100;
 Ntheta = 100;
 %time division
 tmax=1;
 dt=tmax/50;
 t=0:dt:tmax;
 Nt=length(t);
 %space discretization
 Rmax=10;
 %gridposition
 dr=Rmax/(Nr-1);
 dtheta= (2*pi)/(Ntheta-1);
 r=0:dr:Rmax-dr;
 theta=0:dtheta:((2*pi)-dtheta);
 Nr=length(r);
 Ntheta=length(theta);
 x= r.* cos(theta)';
 y= r.* sin(theta)';
 T= ones((Nr),(Ntheta),Nt);
%Initialisation with Dirichlet boundaries
 for i=1:1:(Nr)
 for j=1:1:(Ntheta)
 T(i,j,1)=45;
 if(r(j)==Rmax )
 T(i,j,:)=50;
 end
 end
 end
 T0=T(:,:,1);
 %Calculate A matrix
 d=(ones(Nr,1));
 dr1=d;
 drminus1=d;
 for i=2:length(d)
 dr1(i)= (d(i-1)*2)+(1/(i-1));
 drminus1(i)=(d(i-1)*2)+(1/(i-1));
 end
  d(1)=1;
 Ar=(alpha/2)*(1/(dr^2))* spdiags([dr1 -4*d drminus1],-1:1,Nr,Nr);
 Ar=full(Ar);
  d=(ones(Ntheta,1));
 Aphi= (alpha)*(1/(dtheta^2))*spdiags([d -2*d d],-1:1,Ntheta,Ntheta);
  Aphi=full(Aphi);
 N=Nr;
 if Nr>Ntheta
 Ar1=Ar;
 Aphi1=[horzcat(Aphi,zeros(length(Nr-Ntheta),Ntheta));zeros(Nr)];
 elseif Ntheta>Nr
 Aphi1=Aphi;
 Ar1=[horzcat(Ar,zeros(length(Ntheta-Nr),Nr));zeros(Ntheta)];
 N=Ntheta;
 else
 Ar1=Ar;
 Aphi1=Aphi;
 end
 A=Ar1+Aphi1;
 A(end,1)=A(end,end)*(dr^2);
 A(1,end)=A(1,1)*(dr^2);
 M=zeros(size(A));
 q=0;
 for p=t
 q=q+1;
 M=expm((A).*p);
  T(1,1,q)=mean(T0(2,:));
  disp(['t=','0'+(q-1)]);
  disp(T(1,1,q));
    for j=2:1:(Nr)
 for phi=1:1:(Ntheta)
 if(abs(complex(x(phi),y(phi)) )==Rmax)
 T(j,phi,:)=50;
 disp(T);
 end
 if(abs(r(j))==0 )
 T(j,phi,q)=mean(T0(:,1,q));
 end
  end
 end
 T(:,:,q)=M*(T0)';
 T(1,1,q)=mean(T0(2,:));
  for j=2:1:(Nr)
 for phi=1:1:(Ntheta)
 if(abs(complex(x(phi),y(phi)) )==Rmax)
 T(j,phi,:)=50;
 disp(T);
 end
 if(abs(r(j))==0 )
 T(j,phi,q)=mean(T0(:,1,q));
 end
  end
 end
 T0=T(:,:,q);
 end
 T0=T(:,:,1);
 Tfinal=T(:,:,q);
  %plotting
 figure(1);
  for i=1:Nt
    surf(x(1:5:Nr,1:5:Ntheta),y(1:5:Nr,1:5:Ntheta),T(1:5:Nr,1:5:Ntheta,i));
    shg
 end

J. Alex Lee 2020-9-28

It's hard to read this code, but are you formulating finite differences in cartesian space, or in cylindrical space?

My suggestion was to formulate the heat equation in cylindrical coordinates and just use a "cartesian" meshgrid in rho-theta space. This may not be so great, because you'll get tiny mesh near center of circle.

Vignesh Ramakrishnan 2020-9-29