Parabolic Partial Differential Equations 
Crank-Nicolson Method

Sri Harsha Garapati, Autar Kaw
University of South Florida
United States of America
kaw@eng.usf.edu

Introduction

This worksheet demonstrates the use of Mathematica to illustrate the Crank-Nicolson method of solving parabolic partial differential equations.

Problem Statement

Use the implicit method to solve for the temperature distribution of a long, thin rod with a length of 0.05 metres subjected to a temperature of 100 degrees centigrade on the left end and 25 degress centigrade on the right end. Find the temperature distribution in the rod from t=0 secs to t=9 secs. The following values: K'=54 W/(m.k), C=490 J/(Kg-K) and density=7800 Kg/m^3, deltax=0.01 cm, deltat=3s. At t=0, the temperature of the rod is 20 degress centigrade worksheet demonstrates the use of Mathematica to illustrate the explicit method of solving parabolic partia; differential equations.

Input

The following are the inputs to this problem
      a) Length of the rod, l
      b) Thermal conductivity of the material of the rod, k
      c) Specific heat of the material of the rod, c
      d) density of the material of the rod, density
      f) Temperature on the left end of the rod, T1
      g) Temperature on the right end of the rod, T2
      h) Intial temperature of the rod, Ti
      i) number of divisions of the rod along it length(mesh size), n
      j) End time of the analysis, t
      k) Time step, dt

l=0.05;
k=54;
c=490;
density=7800;
T1=100;
T2=25;
Ti=20;
n=5;
t=9;
dt=3;
Print["Length of the rod=",l];
Print["Thermal conductivity of the material of the rod=",k];
Print["Specific heat of the material of the rod=",c];
Print["Density of the material of the rod=",density];
Print["Temperature on the left end of the rod=",T1];
Print["Temperature on the right end of the rod=",T2];
Print["Initial temperature of the rod=",Ti];
Print["Number divisions along the length of the rod ( mesh size)=",n];
Print["End time of the analysis=",t];
Print["Time step choosen for the analysis=",dt];

Calculations

a) Number of timesteps, runs
b) Length of the element, dx
c) Thermal diffusivity, alpha
d) lambda

runs=(t/dt)+1;
dx=(l/n);
alpha=k/(c*density);
lambda=alpha*dt/(dx)^2;
Print["Number of time steps=",runs];
Print["Edge length of the element=",dx];
Print["Thermal Diffusivity of the material of the rod=",N[alpha]];

Explicit Method

Creating an Array which stores the nodal locations of the nodes

L=Table[(i-1)*dx,{i,1,n+1,1}];

Creating a matrix to store the temperatures at each node at each time step
In the temperature matrix each column represents the node and the row represents the respective time step i.e. first row represents the nodal temperatures of the nodes at the intial time (t=0) and the final row reprsents the nodal temperatures at the end time of the analysis

T=Table[0,{i,1,runs},{j,1,n+1}];

Assigning the boundary temperatures to the Temperature matrix

Do[T[[i,1]]=T1;T[[i,n+1]]=T2,{i,1,runs,1}]

Assigning the the intial temperature (Intial conditions) to the temperature matrix

Do[T[[1,j]]=Ti,{j,2,n,1}]
Print["Temperature Matrix at initial time is, T="];Print[TableForm[T]];

Programming the explicit method

count=0;
For[j=1,j≤runs-1,j=j+1,
count=count+1;
temp=Table[0,{i,1,n-1},{j,1,n-1}];
cons=Table[0,{i,1,n-1},{j,1,1}];
For[i=2,i≤n,i=i+1,
k=i-1;
If[k==1,temp[[k,k]]=2*(1+lambda);temp[[k,k+1]]=-lambda;cons[[k,1]]=(lambda*T[[j,i-1]])+(2*(1-lambda)*T[[j,i]])+(lambda*T[[j,i+1]])+(lambda*T[[j+1,i-1]])];
If[k!=1 &&  k!=(n-1),temp[[k,k-1]]=-lambda;temp[[k,k]]=2*(1+lambda);temp[[k,k+1]]=-lambda;cons[[k,1]]=(lambda*T[[j,i-1]])+(2*(1-lambda)*T[[j,i]])+(lambda*T[[j,i+1]])];
If[k==(n-1),temp[[k,k-1]]=-lambda;temp[[k,k]]=2*(1+lambda);cons[[k,1]]=(lambda*T[[j,i-1]])+(2*(1-lambda)*T[[j,i]])+(lambda*T[[j,i+1]])+(lambda*T[[j+1,i+1]])];
]
TableForm[temp];
sol=Table[Inverse[temp].cons];
For[i=2,i≤n,i=i+1,
T[[j+1,i]]=sol[[i-1,1]];
]

]

Displaying the temperature matrix
Each row represents the time step and the column represents the node number.
i.e. First row represents the intial time (zero time step) and last row repsents the end time of analysis .
Ex: Element T(2,3) represents the nodal temperature of the third node at the first time step

TableForm[T]

Plotting

For[i=1,i≤runs,i=i+1,
B=Table[T[[i,j]],{j,1,n+1}];
dataToPlot=Transpose[{L,B}];

frames[i]=ListLinePlot[dataToPlot,AxesLabel->{length of the rod, Temperature},PlotRange->All];
]

ListAnimate[Table[frames[i],{i,1,runs}],{i,1,runs,1}]