Crank Nicholson method for cylindrical co-ordinates

32 views (last 30 days)
Matthew Hunt on 17 Jan 2020
Commented: Matthew Hunt on 24 Jan 2020
I am trying to solve the heat equation in cylindrical co-ordinates using the Crank-Nicholson method, the basic equation along with boundary/initial conditions are:
The numerical algorithm is contained in the document. The code I wrote for it is:
%This solves the heat equation with source term in cylindrical
%co-ordinates.
%T_t=(k/r)*(rT_r)_r+Q(t,r)
N_t=200; %Time steps
t=linspace(0,2,N_t)'; %Time co-ordinate
dr=r(2);dt=t(2); %Increments
%Solution matrix
T=zeros(N_t,N_r);
%Physical characteristics
k=0.01; %Thermal conductivity
h=0.01; %Thermal exchange constant
theta=0.1;
A=k*dt/(2*dr*2); B=k*dt/(2*dr); %Useful constants
Q=-0.1*t.^2.*exp(-t);
%Construct the matrices
%(i+1)th matrix
a_1=-(theta+2*A)*ones(N_r,1);
a_1(1)=-(4*A-theta);
a_1(N_r)=-(2*h*dr*(A+B/r(N_r))+2*A+theta);
a_2=A+B./r(1:N_r-1);
a_2(1)=4*A;
a_3=A-B./r(2:N_r);
a_3(N_r-1)=2*A;
alpha=diag(a_1)+diag(a_2,1)+diag(a_3,-1);
%ith matrix
b_1=(2*A-theta)*ones(N_r,1);
b_1(1)=4*A-theta;
b_1(N_r)=2*A-theta+2*h*dr*(A+B/r(N_r));
b_2=-a_2;
b_3=-a_3;
beta=diag(b_1)+diag(b_2,1)+diag(b_3,-1);
%Constructing the solution:
b=zeros(1,N_r);
for i=2:N_t
b=-0.5*(Q(i,:)'+Q(i-1,:)')*dt;
b(N_r)=2*dr*(Q(i-1)+Q(i))*(A+B/r(N_r));
u=T(i-1,:)';
C=beta*T(i-1,:)'+b';
x=alpha\C;
T(i,:)=x';
end
for i=1:N_t
plot(r,T(i,:));
pause(0.1);
end
There is a great deal of oscillations in the region of r=0 and I'm not sure how to get rid of it. Any suggestions?
Matthew Hunt on 24 Jan 2020
There should be a heat source term in the equation.