Accelerating the pace of engineering and science

# Partial Differential Equation Toolbox

## Heat Distribution in Radioactive Rod

This example shows how a 3D axisymmetric model can be analyzed using Partial Differential Equation Toolbox™ which supports only 2D analysis. The model geometry, material properties, and boundary conditions must all be symmetric about a single axis for this simplification from 3D to 2D to be appropriate.

This particular example shows heat transfer in a rod with a circular cross section. There is a heat source at the left end of the rod and a fixed temperature at the right end. The outer surface of the rod exchanges heat with the environment due to convection. In addition, heat is generated within the rod due to radioactive decay.

We would like to calculate the temperature in the rod as a function of time. The parabolic equation describing heat transfer is

where are the density, specific heat, and thermal conductivity of the material, respectively, is the temperature, and is the heat generated in the rod.

Since the problem is axisymmetric, it is convenient to write this equation in a cylindrical coordinate system.

where , and are the three coordinate variables of the cylindrical system. Because the problem is axisymmetric, and after multiplying by the equation can be rewritten

The equation can be converted to the form supported by PDE Toolbox if is defined as and is defined as . Rewriting the above equation gives

Definition of PDE Coefficients

The rod is composed of a material with the following properties.

k = 40; % thermal conductivity, W/(m-degree C)
rho = 7800; % density, kg/m^3
cp = 500; % specific heat, W-s/(kg-degree C)
% PDE Toolbox allows the coefficients to be input as string expressions
c = sprintf('%g*y', k);
d = sprintf('%g*y', rho*cp);
q = 20000; % heat source, W/m^3
f = sprintf('%g*y', q);
a = 0;


Geometry and Mesh

The 2D model is a rectangular strip whose y-dimension extends from the axis of symmetry to the outer surface and x-dimension extends over the actual length of the rod (from -1.5 m to 1.5 m). The geometry and mesh for this rectangular section are easily defined by specifying the x and y locations of the four corners as shown below.

g = decsg([3 4 -1.5 1.5 1.5 -1.5 0 0 .2 .2]');


When we define boundary conditions below, it is necessary to know the edge numbers for the boundary edges of the geometry. A convenient way to obtain these edge numbers is to plot the geometry using pdegplot with option edgeLabels set to on.

figure;pdegplot(g,'edgeLabels','on');axis equal;xlim([-2 2]);
title 'Rod Section Geometry With Edge Labels Displayed';
hmax = .1; % element size
[p, e, t] = initmesh(g, 'Hmax', hmax);
figure; pdeplot(p,e,t); axis equal;
title 'Rod Section With Triangular Element Mesh'


Boundary Conditions

The boundary conditions are defined in the function below which is referred to in the PDE Toolbox documentation as a boundary file. The edge at equal zero is along the axis of symmetry so there is no heat transfered in the direction normal to this edge. This boundary is modeled by the default, zero Neumann boundary condition. Boundary conditions for the other three edges are described in the function shown at the end of this example.

b = @boundaryFileRadioactiveRod;


In transient problems of this type it is often useful to first compute the steady state solution-- the solution to the time-independent, elliptic equation. If the final time in the transient analysis is sufficiently large, the transient solution at the final time should be close to this steady state solution. This provides a valuable check on the accuracy of the transient analysis.

u = assempde(b,p,e,t,c,a,f);
figure; pdeplot(p, e, t, 'xydata', u, 'contour', 'on'); axis equal;


Transient Solution

The time-dependent solution is computed from to seconds using the parabolic function.

tfinal = 20000;
tlist = 0:100:tfinal;
u0 = 0; % initial temperature in the rod is zero
u = parabolic(u0, tlist, b, p,e,t,c,a,f, d);
figure; pdeplot(p, e, t, 'xydata', u(:,end), 'contour', 'on'); axis equal;
title(sprintf('Transient Temperature at Final Time (%g seconds)', tfinal));

96 successful steps
0 failed attempts
194 function evaluations
1 partial derivatives
20 LU decompositions
193 solutions of linear systems


As expected, the steady state solution and the transient solution at 20000 seconds are in close agreement. This can be seen by comparing the two figures.

The third figure below shows the temperature at the left end of the rod as a function of time. The outer surface of the rod is exposed to the environment with a constant temperature of 100 degrees-C. Consequently, when the surface temperature of the rod is less than 100, the rod is being heated by the environment and when greater than 100, cooled. When the rod temperature is less than 100 degrees, the outer surface is slightly warmer than the inner axis but when the temperature is around 100 degrees, the outer surface becomes cooler than the interior of the rod.

figure; plot(tlist, u(1,:)); hold all; plot(tlist, u(4,:), '--');
title 'Temperature at Left End as a Function of Time'
xlabel 'Time, seconds'
ylabel 'Temperature, degrees-C'
grid on;
legend('Center Axis', 'Outer Surface', 'Location', 'SouthEast');


function [ q, g, h, r ] = boundaryFileRadioactiveRod( p, e, u, time )
%   [ q, g, h, r ] = BOUNDARYFILERADIOACTIVEROD( p, e, u, time ) returns the
%   Neumann BC (q, g) and Dirichlet BC (h, r) matrices for the
%   p is the point matrix returned from INITMESH
%   e is the edge matrix returned from INITMESH
%   u is the solution vector (used only for nonlinear cases)
%   time (used only for parabolic and hyperbolic cases)
%

%   Copyright 2012 The MathWorks, Inc.

N = 1;
ne = size(e,2);
q = zeros(N^2, ne);
g = zeros(N, ne);
h = zeros(N^2, 2*ne);
r = zeros(N, 2*ne);

% Neumann boundary conditions are defined at the mid-points of
% the element edges so we calculate the coordinates of these points
% for use in the expressions below.
xy1 = p(:,e(1,:));
xy2 = p(:,e(2,:));
xyMidEdge = (xy1+xy2)/2;

for i=1:ne
switch(e(5,i))
case 2
% right edge, temperature is prescribed to be 100 degrees-C
h(1,i) = 1;
h(1,i+ne) = 1;
r(1,i) = 100;
r(1,i+ne) = 100;
case 3
% outer edge, convection to 100 degree environment
y = xyMidEdge(2,i);
q(1,i) = 50*y;
g(1,i) = 50*y*100;
case 4
% left edge, heat flux of 5000 W/m^2
y = xyMidEdge(2,i);
g(1, i) = 5000*y;
end

end