Crank-Nicholson solution to cylindrical heat equation

  • A
  • Thread starter hunt_mat
  • Start date
  • #1
hunt_mat
Homework Helper
1,720
16

Summary:

I used a Crank-Nicholson method to solve a radially symmetric heat equation. I get some odd results. I'm not sure if it's my 1) Numerical scheme, 2) Parameters chosen or 3) Code which is wrong.

Main Question or Discussion Point

Hi,

I am solving the radially symmetric heat equation with an internal heat source(this is meant to model the heating of a cylindrical battery). It's meant to model heat in a cylinder with conduction to the environment, so my outer boundary condition is Newton's law of cooling. The $T$ in the equation is excess temperature from the environment. I used a Crank-Nicholson which I think I've got right. I plotted a cross-section to see if it looked okay and it didn't basically and now I'm not sure what's going on. I've checked my code, and it matches my numerical algorithm. I think I have my physics correct but again I could be wrong. My code is:
Code:
%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,20,N_t)'; %Time co-ordinate
dr=r(2);dt=t(2); %Increments

T=zeros(N_t,N_r);
a_2=zeros(N_r-1,1);
a_3=zeros(N_r-1,1);

%Physical characteristics
k=3; %Thermal conductivity
h=2.5; %Thermal exchange constant
theta=1;
A=k*dt/(4*dr); B=k*dt/(2*dr^2); %Useful constants
Q=10*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(2:N_r-1)=A+B./r(2:N_r-1);
a_2(1)=4*A;
a_3(1:N_r-2)=A-B./r(2:N_r-1);
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)=-(theta+4*A);
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;
    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
Can anyone see where I have made a glaring error?
 

Attachments

Answers and Replies

Related Threads for: Crank-Nicholson solution to cylindrical heat equation

  • Last Post
Replies
0
Views
2K
Replies
4
Views
4K
Replies
6
Views
1K
  • Last Post
Replies
5
Views
4K
  • Last Post
Replies
1
Views
4K
Replies
1
Views
4K
  • Last Post
Replies
2
Views
9K
Replies
4
Views
4K
Top