- #1

- 1,782

- 32

- 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

Can anyone see where I have made a glaring error?

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?