Crank-Nicholson solution to the cylindrical heat equation

[hunt_mat]

TL;DR 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.

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:

%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?

 

[Chestermiller]

I think the issue here is how you handle the boundary conditions, in particular, the boundary condition at r = 0.

I would approach the formulation of the difference equations a little differently from the way you have done it. It would be a two step process, first using the method of lines to discretize the differential equation spatially into a coupled set of 1st order ODEs in time, and then secondly to apply the Crank Nicholson formalism to the set of ODEs.

I will only show the first step here. The 2nd step is then obvious. So, applying the 1st step to the general region, I get $$\frac{dT_j}{dt}=\alpha \frac{T_{J+1}(1+1/(2j))-2T_j+T_{j-1}(1-1/(2j))}{(\Delta r)^2}+q$$where $\alpha=\kappa/\theta$ and $q=Q/\theta$

Now for the BC at r = 0. Since dT/dr is zero there, we can approximate the temperature profile in the vicinity of r = 0 by a series of quadratic accuracy $$T=A(t)+B(t)r^2\tag{1}$$, with $T=T_0$ at r = 0 and $T=T_1$ at $r=\Delta r$. So $A=T_0$ and $B=\frac{T_1-T_0}{(\Delta r)^2}$, From Eqn. 1, it follows that $$\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right)=4B(t)$$. So, for the BC at r = 0, the spatially discretized differential equation becomes: $$\frac{dT_0}{dt}=\alpha \frac{4(T_1-T_0)}{(\Delta r)^2}+q$$This is how the 2nd order accurate spatially discretized boundary condition should be applied at r = 0.

 

[hunt_mat]

Hmmm, I was hoping that I could do the same thing as I did with the spherical diffusion equation. I am aware of the method of lines but that's all I am. I've seen the method for dealing with the r=0 term before and I applied it here, so why does it work for spherical and not cylindrical? I've always viewed numerics as a bit of a dark art.

 

[Chestermiller]

Hmmm, I was hoping that I could do the same thing as I did with the spherical diffusion equation. I am aware of the method of lines but that's all I am. I've seen the method for dealing with the r=0 term before and I applied it here, so why does it work for spherical and not cylindrical? I've always viewed numerics as a bit of a dark art.

It's. not clear what you are saying. Are you saying that you tried the discretization method I suggested in spherical coordinates (where the 4 becomes a 6), or are you saying that you tried something entirely different. I've used the method I've suggested a huge number of times for both spherical coordinates and cylindrical coordinates, and it's worked flawlessly in all cases. I don't really view numerics as a dark art. It seems pretty straightforward to me.

If you haven't tried the method I suggested, please do. It should be easy enough to set up within your current framework.

 

[hunt_mat]

No, I computed the limit as $r\rightarrow 0$ as:
$$\lim_{r\rightarrow 0}\frac{1}{r}\frac{\partial u}{\partial r}=\frac{\partial^{2}u}{\partial r^{2}}$$

So for the first point, I would simply discretise a cartesian diffusion equation and I treat everything else normally.

Also for the method of lines, you basically have the equation:

$$\frac{d\mathbf{T}}{dt}=\mathbf{A}\mathbf{T}$$
Where A is some matrix, could I simply write $\mathbf{A}=\mathbf{P}^{T}\mathbf{D}\mathbf{P}$, where $\mathbf{D}$ is a diagonal matrix, then if I just multiply my equation through by $\mathbf{P}$ I obtain:
$$\frac{d}{dt}(\mathbf{P}\mathbf{T})=\mathbf{D}(\mathbf{P}\mathbf{T})$$
Which results in a set of decoupled equations which are simple to integrate.

I have no doubt that you can get your code to work perfectly, I, on the other hand, have to fight to get something to work.

 

[Chestermiller]

No, I computed the limit as $r\rightarrow 0$ as:
$$\lim_{r\rightarrow 0}\frac{1}{r}\frac{\partial u}{\partial r}=\frac{\partial^{2}u}{\partial r^{2}}$$

So for the first point, I would simply discretise a cartesian diffusion equation and I treat everything else normally.

. This is not a valid discretization of the spatial derivative. The expression I gave, on the other hand actually is a valid discretization. I can write the corresponding Crank-Nicholson version of the equation for the grid point at r = 0 so that you can incorporate it into your code if that is what you need. But I was hoping you would be able to do this on your own. It should take less than 5 minutes for you to change to code. Until you do that, in my judgment, you will not be able to get a satisfactory solution.

Also for the method of lines, you basically have the equation:

$$\frac{d\mathbf{T}}{dt}=\mathbf{A}\mathbf{T}$$
Where A is some matrix, could I simply write $\mathbf{A}=\mathbf{P}^{T}\mathbf{D}\mathbf{P}$, where $\mathbf{D}$ is a diagonal matrix, then if I just multiply my equation through by $\mathbf{P}$ I obtain:
$$\frac{d}{dt}(\mathbf{P}\mathbf{T})=\mathbf{D}(\mathbf{P}\mathbf{T})$$
Which results in a set of decoupled equations which are simple to integrate.

I have no doubt that you can get your code to work perfectly, I, on the other hand, have to fight to get something to work.

No. My equation using the method of lines is of the form:
$$\frac{dT}{dt}=AT+q$$where q is the vector of heat source at the grid points. The Crank-Nicholson approximation to this set of ODEs is:
$$\frac{T^{i+1}-T^i}{\Delta t}=\frac{1}{2}(AT^{i+1}+q^{i+1}+AT^i+q^i)$$

 

[hunt_mat]

Why isn't it a valid discretization? I saw this in various peer-reviewed research articles, so I naturally assumed it was perfectly fine.

So basically I would only need to change a row of the $\mathbf{A}$ matrix and that should be fine? I'm battling with a finite element code at the moment which isn't working but did in the spherical case.

 

[Chestermiller]

Why isn't it a valid discretization? I saw this in various peer-reviewed research articles, so I naturally assumed it was perfectly fine.

I don't want to have to go into that. Just please try the approximation that I showed you.

So basically I would only need to change a row of the $\mathbf{A}$ matrix and that should be fine? I'm battling with a finite element code at the moment which isn't working but did in the spherical case.

Yes, that's all you would have to do.

 

[hunt_mat]

Using the notation from previous posts, my new stencil will now become:
$$-\alpha A_{j,j-1}u_{i+1,j+1}+(\theta-\alpha A_{j,j})-\alpha A_{j,j+1}u_{i+1,j+1}=\alpha A_{j,j-1}u_{i,j-1}+(\theta+\alpha A_{j,j})u_{i,j}+\alpha A_{j,j+1}+\frac{\delta t}{2}(Q_{i+1}+Q_{i})$$

Certain rows and columns will have to be adjusted for boundary conditions. This is basically the code?

 

[Chestermiller]

Using the notation from previous posts, my new stencil will now become:
$$-\alpha A_{j,j-1}u_{i+1,j+1}+(\theta-\alpha A_{j,j})-\alpha A_{j,j+1}u_{i+1,j+1}=\alpha A_{j,j-1}u_{i,j-1}+(\theta+\alpha A_{j,j})u_{i,j}+\alpha A_{j,j+1}+\frac{\delta t}{2}(Q_{i+1}+Q_{i})$$

Certain rows and columns will have to be adjusted for boundary conditions. This is basically the code?

I think you made so typo's in the equation you wrote. It should be exactly the same as you previously obtained, except for row j = 0 (i.e., the $T_0$ equation). Did you not understand that?

 

[hunt_mat]

No, the code will be different. The way it was discretised was different and hence the stencil will be different. Otherwise, what's the difference between the method of lines and normal finite difference?

 

[Chestermiller]

No, the code will be different. The way it was discretised was different and hence the stencil will be different. Otherwise, what's the difference between the method of lines and normal finite difference?

The method of lines is just a process for deriving the finite difference discretization in a two step-process, by first discretizing with respect to the spatial parameter and then discretizing with respect to the time parameter. It should deliver the exact same Crank-Nicholson finite difference equation(s). I will derive the Crank-Nicholson finite difference equation using your parameters and demonstrate that it is exactly the same.

Your starting differential equation is

$$θ\frac{\partial T}{\partial t}=\frac{\kappa}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right)+Q$$​

The spatial finite difference discretization of this is:
$$\theta \frac{dT_j}{dt}=\kappa\left[\frac{T_{j+1}-T_{j-1}}{2r_j\Delta r}+\frac{(T_{j-1}-2T_j+T_{j+1}}{(\Delta r)^2}\right]+Q=
\kappa \frac{T_{j+1}(1+1/(2j))-2T_j+T_{j-1}(1-1/(2j))}{(\Delta r)^2}+Q$$where, in the second term of the rhs, we have made use of the relationship that rj=jΔrrj=jΔr.

I can now see that you have an error in Eqn. 4 of your write-up. There should be a 2 in the denominator of the 1st- and the 4th terms in brackets. Do you see that? So the value of your parameter B should be $B=\frac{\kappa \delta t}{4\delta r}$, not $B=\frac{\kappa \delta t}{2\delta r}$.

I will continue after we get agreement on this.

 

[Chestermiller]

Are you labeling the grid point at r = 0 as j = 1 and the grid point at r = R as j = N, so that there are (N - 1) grid intervals and $\Delta r=R/(N-1)$?

 

[hunt_mat]

I start from $j=1$ which is in-line with the MATLAB convention. So I do

r=linspace(0,R,N)

, this makes $r_{1}=0$.

 

[Chestermiller]

I found some sign errors in your algebra in going from Eqn. 12 to Eqn. 14. Please check.

 

[hunt_mat]

Using the notation before my stencil is:
$$(\theta\delta_{jk} -\alpha A_{j,k})T_{i+1,k}=(\alpha A_{j,k}+\theta\delta_{jk})u_{i,k}+\frac{\delta t}{2}(Q_{i+1}+Q_{i})$$

The top row of $\mathbf{A}$ is $(-4,4,0,\ldots,0)$, the bottom row is $(0,\dots,0,2,(h\delta r/\kappa-2))$

That's the code in a nutshell I believe.

 

[Chestermiller]

Using the notation before my stencil is:
$$(\theta\delta_{jk} -\alpha A_{j,k})T_{i+1,k}=(\alpha A_{j,k}+\theta\delta_{jk})u_{i,k}+\frac{\delta t}{2}(Q_{i+1}+Q_{i})$$

The top row of $\mathbf{A}$ is $(-4,4,0,\ldots,0)$, the bottom row is $(0,\dots,0,2,(h\delta r/\kappa-2))$

That's the code in a nutshell I believe.

I concur with the top row if alpha is what you previously called A in your write-up.

And what is your comment regarding my assertion about your parameter B being high by a factor of 2?

 

[Chestermiller]

It looks like you have the A and B parameter definitions switched in your code.

 

[hunt_mat]

This is the code I'm working with now:

%This solves the heat equation with source term in cylindrical
%co-ordinates.

%T_t=(k/r)*(rT_r)_r+Q(t,r)
N_r=800; %Radial steps
N_t=1000; %Time steps

r=linspace(0,1,N_r); %Radial co-ordinate
t=linspace(0,3,N_t)'; %Time co-ordinate
dr=r(2);dt=t(2); %Increments

%Solution matrix
T=zeros(N_t,N_r);

%Physical characteristics
k=1; %Thermal conductivity
h=1; %Thermal exchange constant
theta=1;
alpha=k*dt/(2*dr*2);  %Useful constants
Y_1=sech(t)*(r.^2.*(1-r).^2);
Y_2=-2*tanh(t)*(8*r.^2-9*r+2);
Q=Y_1+Y_2;
a_1=-2*ones(1,N_r);
a_2=1+dr./r(1:N_r-1);
a_3=1-dr./r(2:N_r);
A=diag(a_1)+diag(a_2,1)+diag(a_3,-1);
A(1,1)=-4;A(1,2)=4;
A(N_r,N_r-1)=2;A(N_r,N_r)=h*dr/kappa-2;
A_plus=theta*eye(N_r)-alpha*A;
A_minus=theta*eye(N_r)+alpha*A;

%Constructing the solution:
b=zeros(1,N_r);
analytical=tanh(t)*(r.^2.*(1-r).^2);
for i=2:N_t
    b=0.5*(Q(i,:)'+Q(i-1,:)')*dt;
    u=T(i-1,:)';
    C=A_minus*T(i-1,:)'+b;
    x=A_plus\C;
    T(i,:)=x';
end

for i=1:N_t
    plot(r,T(i,:));
    hold on
    plot(r,analytical(i,:));
    pause(0.1);
    hold off
end

 

[Chestermiller]

This is the code I'm working with now:

%This solves the heat equation with source term in cylindrical
%co-ordinates.

%T_t=(k/r)*(rT_r)_r+Q(t,r)
N_r=800; %Radial steps
N_t=1000; %Time steps

r=linspace(0,1,N_r); %Radial co-ordinate
t=linspace(0,3,N_t)'; %Time co-ordinate
dr=r(2);dt=t(2); %Increments

%Solution matrix
T=zeros(N_t,N_r);

%Physical characteristics
k=1; %Thermal conductivity
h=1; %Thermal exchange constant
theta=1;
alpha=k*dt/(2*dr*2);  %Useful constants
Y_1=sech(t)*(r.^2.*(1-r).^2);
Y_2=-2*tanh(t)*(8*r.^2-9*r+2);
Q=Y_1+Y_2;
a_1=-2*ones(1,N_r);
a_2=1+dr./r(1:N_r-1);
a_3=1-dr./r(2:N_r);
A=diag(a_1)+diag(a_2,1)+diag(a_3,-1);
A(1,1)=-4;A(1,2)=4;
A(N_r,N_r-1)=2;A(N_r,N_r)=h*dr/kappa-2;
A_plus=theta*eye(N_r)-alpha*A;
A_minus=theta*eye(N_r)+alpha*A;

%Constructing the solution:
b=zeros(1,N_r);
analytical=tanh(t)*(r.^2.*(1-r).^2);
for i=2:N_t
    b=0.5*(Q(i,:)'+Q(i-1,:)')*dt;
    u=T(i-1,:)';
    C=A_minus*T(i-1,:)'+b;
    x=A_plus\C;
    T(i,:)=x';
end

for i=1:N_t
    plot(r,T(i,:));
    hold on
    plot(r,analytical(i,:));
    pause(0.1);
    hold off
end

Y'know, it would have been very considerate of my valuable time to, from the onset, have given me a correct version of the writeup and a corresponding latest version of your code. I've wasted a lot of time looking over your flawed versions.

 

[hunt_mat]

That was my new code which I wrote yesterday.

Thank you for your time. I will use pdepe, as it seems that doing anything in cylindrical co-ordinates is a complete nightmare which I cannot get to work.

 

[hunt_mat]

These are my calculations for the method. I haven't made any mistakes and it simply doesn't work.

 

[Chestermiller]

These are my calculations for the method. I haven't made any mistakes and it simply doesn't work.

If you didn't make any mistakes, it would work. Try a calculation in which you use only 21 grid points and choose delta t such that it makes alpha close to 1.

 

[Chestermiller]

I have some ideas on how you can proceed to get your CN program to work properly.

Phase 1: To test out your spatial finite difference scheme, solve the problem using forward Euler (rather than CN ) with a small enough step size that the solution remains stable.

Phase 2: Use the CN equation, but solve the matrix using predictor-corrector (successive substitution) at each time step, rather than matrix inversion. The predictor for the first time step would be forward Euler. The corrector equation would involve only the diagonal elements, while the other values at the end of the time interval would be those from the previous iteration. Iterate at each time step until the corrector has converged.

I also recommend for debugging purposes that you do all these tests with only 3 spatial grid points. The CN solution at each time step would then require the solution of only 3 simultaneous linear equations, which could be done by hand (to check that the matrix was set up and solved properly).