 #1
hunt_mat
Homework Helper
 1,720
 16
Summary:
 I used a CrankNicholson 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 CrankNicholson which I think I've got right. I plotted a crosssection 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:
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 CrankNicholson which I think I've got right. I plotted a crosssection 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
%coordinates.
%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 coordinate
t=linspace(0,20,N_t)'; %Time coordinate
dr=r(2);dt=t(2); %Increments
T=zeros(N_t,N_r);
a_2=zeros(N_r1,1);
a_3=zeros(N_r1,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*Atheta);
a_1(N_r)=2*h*dr*(A+B/r(N_r))2*Atheta;
a_2(2:N_r1)=A+B./r(2:N_r1);
a_2(1)=4*A;
a_3(1:N_r2)=AB./r(2:N_r1);
a_3(N_r1)=2*A;
alpha=diag(a_1)+diag(a_2,1)+diag(a_3,1);
%ith matrix
b_1=(2*Atheta)*ones(N_r,1);
b_1(1)=(theta+4*A);
b_1(N_r)=2*Atheta+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(i1,:)')*dt;
u=T(i1,:)';
C=beta*T(i1,:)'+b;
x=alpha\C;
T(i,:)=x';
end
for i=1:N_t
plot(r,T(i,:));
pause(0.1);
end
Attachments

91.5 KB Views: 15

112.2 KB Views: 15