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

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

[J. Alex Lee]

Finite difference is ill suited for non rectangular meshes, unless you can solve in cylindrical coordinates and make a "rectangular" mesh in r,theta coordinates

[Vignesh Ramakrishnan]

On proceeding that way, I am getting a groove along theta=0

[J. Alex Lee]

your original question was about generating a mesh. did you succeed there? results?

[Vignesh Ramakrishnan]

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]

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]

I have used circular coordinates(for cylindrical, we have a height term, which I have not considered).

I don't have an issue with the mesh circle, but there is an unnecessary groove along theta=0.

The cartesian conversion is done, just prior to plotting.