Crank Nicholson method for cylindrical co-ordinates

32 views (last 30 days)
Matthew Hunt
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_r=300; %Radial steps
N_t=200; %Time steps
r=linspace(0,1,N_r); %Radial co-ordinate
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?
  5 Comments
Matthew Hunt
Matthew Hunt on 24 Jan 2020
There should be a heat source term in the equation.

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!