Crank Nicholson method for cylindrical co-ordinates

18 views (last 30 days)
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
J. Alex Lee
J. Alex Lee on 22 Jan 2020
Edited: J. Alex Lee on 22 Jan 2020
Maybe there is something funky about your BCs...your final condition seems to imply that you want the temperature difference with the outside where outside has 0 temperature, but your intial temperature is 0. Don't you want to either start with a non-zero temperature field, or modify your outside boundary condition as
or something like that, where T_0 is a non-zero constant?
It doesn't solve your problem, but maybe helpful to understand the steady state solution, and see if your FD scheme will suffer oscillations with the steady state solution.
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!