Project 2

Jeffrey A. Anderson

December 7, 2018
Jeffrey A. Anderson ME560: Project 2

Figure 1: Nuclear Fuel Rod Assembly [1].

Introduction
Uranium dioxide is a common radio-active ceramic material used to make fuel
pellets for nuclear reactors. As with other ceramic materials, it may undergo
brittle fracture if tensile stresses exceed it’s maximum tensile strength. Figure
1 shows a typical fuel rod assembly comprised of individual rods. Each rod has
an outer cladding that is filled with uranium dioxide pellets stacked one on top
of another. The cladding is designed to help prevent the spread of radiological
contamination from the fuel pellets.
Fuel rods function by producing heat. They produce some level of radiation
as a result of the natural radioactive decay process. As more radioactive ma-
terial is brought into close proximity, neutrons from the decay reaction become
more likely to strike a nucleus and split another atom, causing a chain reaction.
Control rods made of neutron absorbing material are inserted between the fuel
rods and are drawn in and out as necessary to control the chain reaction.
The radioactive energy from this decay process in released in the form of
heat, and is treated in heat transfer as internal heat generation, q, in the heat
conduction equation. Prior to being assembled, the heat generation of a single
pellet would be much lower than after it has been put together with a large
amount of other pellets in a reactor core and the control rods partially with-
drawn.
Due to this change in internal heat generation within the pellets, a thermal
gradient develops causing internal stresses to develop. If these stresses exceed
the fracture strength of the uranium dioxide, then failure occurs as shown in Fig-
ure 2a. Stresses and stress drivers may be analyzed to help in the development

1
Jeffrey A. Anderson ME560: Project 2

Figure 2: Left: (a) Failed Ceramic Fuel Pellet [2]. Right: (b) Fuel Rod
Schematic [3].

of more robust nuclear fuels.

This project will develop a simplified stress model based on constant material
properties. Additional predictive ability could be obtained if material properties
were considered to be a function of temperature.

Methods
Real World to Solid Mechanics
A detailed fuel rod schematic is shown in Figure 2b giving several representative
dimensions and identifying the various aspects of the general design. The ge-
ometry is cylindrical, and as such the polar coordinate system will be employed
in setting up the problem and developing the solution.
It will be assumed that the stress state is approximately planar in the r − θ
directions. This assumption is justified by the relatively low stress that will
be caused by the fuel stack hold-down spring. Any thermal expansion in the
z-direction will be compensated for by this spring with only minimal resulting
stress (assuming the spring has a very low stiffness compared to the ceramic).
It can also be seen from Figure 2, that the fuel pellets are not directly in
contact with the water flowing through the reactor core. For the purposes of
this project, a the cladding will be ignored to simplify geometry and the heat
transfer equations.
Further, material properties will be assumed to be constant. This avoids any
non-linearity that may be introduced by using properties that are a function of
temperature. While this simplification will not be justified, it will be accepted
as a necessity for reaching a solution within the scope of this project.

2
Jeffrey A. Anderson ME560: Project 2

Geometry and Kinematics

Since the geometry is cylindrical, and there is nothing special about any par-
ticular angle θ, the axi-symmetric case will be assumed. This results in several
very convenient simplifications of the governing equations.
The dimensions used for final solution evaluation and plotting will be taken
from Figure 2. Specifically, the cylinder radius will be assumed to be 0.005
meters. Since a planar stress state is assumed, the length dimension will be
ignored.

Constitutive Model
The material will be assumed to be linear elastic and isotropic. Therefore, the
generalized formulation of Hooke’s Law will be employed in the solution as
1
rr = (σrr − νσθθ ) + α∆T (1)
E
1
θθ = (σθθ − νσrr ) + α∆T (2)
E

Governing Equations
The static equilibrium equation for polar plane stress was used in the solution
technique. This is given by
∂σrr 1
+ (σrr − σθθ ) = 0 (3)
∂r r
The form of the strain compatibility equations that will be used is derived
from the relations
du
rr = (4)
dr
and
u
θθ = (5)
r
which leads to the compatibility equation
dθθ 1
= (rr − θθ ) (6)
dr r
The radial temperature distribution for a round rod of infinite length having
an internal energy generation is

q̇r2 r2
 
T (r) = o 1 − 2 + Ts (7)
4k ro
Since nuclear reactors are typically cooled using forced convection, the sur-
face temperature is given by

3
Jeffrey A. Anderson ME560: Project 2

q̇ro
Ts = T∞ + (8)
2h
Substituting equation 8 into equation 7 and simplifying yields
 
1 2 2 ro
T (r) = q̇ (r − r ) + + T∞ (9)
4k o 2h
Assuming there exists some reference temperature T0 at which the stress field
is identically zero (neglecting residual stresses created during manufacturing)
then
 
1 2 ro
∆T = q̇ (ro − r2 ) + + T∞ − T0 (10)
4k 2h
which can be substituted into equations 1 and 2.

Stress/Strain Deformations or Tractions

The traction on a given radial face defined by normal vector n is given by

t=σ·n (11)

The normal component of t is calculated using

tr = t · n (12)

And the transverse component of t is found by

tθ = t − tr (13)

Displacement Fields, Deformation Gradients, and Strain

Fields
The displacement gradient tensor is defined as
∂ui
u⊗∇= (14)
∂xj

and the deformation gradient tensor is defined as

∂ui
Fij = δij + (15)
∂xj

As for the strain, we will assume the strains to be small such that the infinites-
imal strain tensor provides a satisfactory approximation of the strain field. It
is given by  
1 ∂ui ∂uj
ij = + (16)
2 ∂xj ∂xi

4
Jeffrey A. Anderson ME560: Project 2

Solution
A summary of the solution is presented here. The full derivation of the solution
is given in the appendix.
First, Hooke’s Law (1 and 2) are substituted into the compatibility equations
(6). The equilibrium equation is then solved for σθθ and substituted to get the
whole thing in terms of σrr . After simplification, the resulting ODE is

d2 σrr 3 dσrr q̇αE

+ = (17)
dr2 r dr 2k
There are two boundary conditions, namely

σrr (ro ) = 0 (18)

and
dσrr
=0 (19)
dr r=0
which when solved, yields a solution of
q̇αE 2
σrr = (r − ro2 ) (20)
16k
Then, going back to solve for σθθ using the re-arranged equilibrium equation
we find that
q̇αE
σθθ = (3r2 − ro2 ) (21)
16k
The stress state is now defined by
 
σrr 0
σ= (22)
0 σθθ

The maximum stress values are found using the first derivative test. This shows
that the maximum stresses in both coordinates appear at r = 0 with values of
q̇αE 2
σrr (0) = − r (23)
16k o
q̇αE 2
σθθ (0) = r (24)
16k o
This result makes sense when looking at the failure fractures in Figure 2a.
Ceramics do not fail in compression as easily as they do in tension. When they
do fail, the failure is often brittle in nature. Since the radial stress is compressive
(positive), it makes sense that most of the fractures would be normal to the θ
stress direction where the stress is tensile (negative). However, it is curious that
there is a fracture ring normal to the radial direction as well. The failure of
this result to indicate the reason behind this may be due to the neglect of the
temperature dependence of the material properties which themselves also affect
the stress field.

5
Jeffrey A. Anderson ME560: Project 2

Results
To evaluate the solution results numerically, reasonable values for the constants
were obtained from the literature. The input q̇ had to be derived from available
data. Starting with rearranging equation 7 to give
4k(T (0) − Ts )
q̇ =
ro2
and inserting approximate temperatures of 400 o C at the surface and 1500 o C
at the center [4], ro = 0.005 m [3], and k = 7.6 W/mK [4] results in q̇ = 1.3376
GW/m3 . The elastic modulus used was E = 2.196 × 105 e−3.025P , where P is
porosity fraction, assumed to be about 0.02, giving E = 206.7 GP a [5]. Finally,
α had to be derived from an equation taken from an Atomic Energy Commission
report [6] that was given as

L = L0 (1 + 6.0 × 10−6 T + 2.0 × 10−9 T 2 + 1.7 × 10−12 T 3 ) (25)

with T in o C. A plot of equation 25 is shown in Figure 3 (left) in terms of

L/L0 . This equation can be rearranged to get

α∆T = 6.0 × 10−6 T + 2.0 × 10−9 T 2 + 1.7 × 10−12 T 3

if T0 is taken to be 0 o C. This gives us a result of

α = 6.0 × 10−6 + 2.0 × 10−9 T + 1.7 × 10−12 T 2 (26)

with T in o C. Plugging in an average expected temperature of 950 o C results

in α = 9.434 × 10−6 o C −1 .
While a constant coefficient of thermal expansion makes the math much
easier, Figure 3 (right) makes it abundantly clear that this is a poor assump-
tion since the thermal coefficient of expansion nearly doubles over the expected
temperature range (400-1500 o C) from edge to core.
The compressive fracture strength for uranium dioxide is given by

σcmax = 1.067 × 1 −3 e−3.505P (27)

Figure 3: Left: Thermal Expansion of Uranium Dioxide. Right: Thermal Ex-

pansion Coefficient of Uranium Dioxide

6
Jeffrey A. Anderson ME560: Project 2

Figure 4: Left: Radial Stress in Pa. Right: Hoop Stress in Pa

where P is again the porosity which is taken to be 0.02 [5]. This gives an ultimate
compressive strength of σcmax = 995 M P a. No value or equations were found
for tensile strength.
Plugging in all of the above constants into the solutions from equations 20
and 21 gives the overall stress distribution as a function of r. The Results are
plotted in Figure 4.

7
Bibliography

[1] Dave Lochbaum. Preventing fuel damage in nuclear

reactors. https://allthingsnuclear.org/dlochbaum/
preventing-fuel-damage-in-nuclear-reactors, January 2014.
[2] D. Banks. Some observations of density and porosity changes in uo2 fuel
irradiated in water-cooled reactors. 54:97–107.
[3] D. Olander. Nuclear fuels - present and future. 389:1–22.

[4] C. Ronchi. Thermal conductivity of uranium dioxide up to 2900 k from

simultaneous measurements of heat capacity and thermal diffusivity. 85:
776–789.
[5] N. Igata and K. Domoto. Fracture stress and elastic modulus of uranium
dioxide including excess oxygen. 45:317–322.

[6] H. G. Wohlers F. A. Halden and R. H. Reinhart. Thermal expansion of

uranium dioxide.

8
Jeffrey A. Anderson ME560: Project 2

Appendix
Substituting Hooke’s law (1 and 2) into compatibility (6) gives
    
d 1 1 2 2 ro
(σθθ − νσrr ) + α q̇ (r − r ) + + T∞ − T0 =
dr E 4k o 2h
 
1 1 1
(σrr − νσθθ ) − (σθθ − νσrr )
r E E
Solving equilibrium (3) for σθθ results in

dσrr
σθθ = r + σrr
dr
and  
dσθθ d dσrr
= r + σrr
dr dr dr
which combines to give
   2   
dσrr ν σrr ν ν dσrr q̇rα 1 d σrr dσrr σrr 1 dσrr
− − + r +σrr − + r +2 − + r +σrr = 0
dr E Er Er dr 2k E dr2 dr Er Er dr

r d2 σrr 1 dσrr 2
   
1 ν ν 1 1 ν ν q̇rαE
+ + + − + σ rr − + − − =0
E dr2 r dr E E E E Er Er Er Er 2k
which then simplifies to

d2 σrr 3 dσrr q̇αE

2
+ =
dr r dr 2k
This ODE will have a homogeneous and a particular solution given by Yh and
Yp , respectively. Yh is found by assuming a solution of the form rn such that

dσrr
= nrn−1
dr
d2 σrr
= n(n − 1)rn−2
dr2
which upon substitution into the ODE gives
3
n(n − 1)rn−2 + nrn−1 = 0
r
n(n − 1)rn−2 + 3nrn−2 = 0
rn−2 (n(n − 1) + 3n) = 0
Since rn−2 cannot be zero without yielding the trivial solution, we get

n(n − 1) + 3n = 0

9
Jeffrey A. Anderson ME560: Project 2

n2 − n + 3n = 0
n2 + 2n = 0
Solving for the roots gives
p
−2 ± 4 − 4(1)(0)
n=
2
n = −2, 0
Such that
C2
Yh = + C1
r2
Now, turning to Yp we try

C3 q̇αE 2
Yp = r
2k
which gives
dYp C3 q̇αE
= r
dr k
d2 Yp C3 q̇αE
2
=
dr k
Upon substitution back into the ODE we get

Solving for C3 results in

C3 q̇αE 3 C3 q̇αE 4C3 q̇αE
+ r=
k r k k
and therefore
1
C3 =
8
And so we reach the solution for Yp

q̇αE 2
Yp = r
16k
Which can be checked by substitution
q̇αE 3 q̇αE q̇αE
+ r=
8k r 8k 2k
Check!
Now, putting Yh and Yp together give us the whole solution

q̇αE 2 C2
σrr = r + 2 + C1
16k r

10
Jeffrey A. Anderson ME560: Project 2

The boundary conditions for the ODE are as follows:

σrr (ro ) = 0
dσrr
=0
dr r=0
Solving the BC for r = 0
dσrr 2C2 q̇αE
= 3 + r
dr r 8k
2C2 q̇αE
+ (0) = 0
(0)3 8k
we find that C2 must be zero or else the solution would not even be defined at
r = 0. This gives
q̇αE 2
σrr = C1 + r
16k
Now for the BC at r = ro
q̇αE 2
σrr (ro ) = C1 + r =0
16k o
q̇αE 2
C1 = − r
16k o
Which gives the final solution for σrr to be
q̇αE 2
σrr = (r − ro2 )
16k
Now we return to the equation
dσrr
σθθ = r + σrr
dr
to solve for σθθ by first taking the derivative of σrr
dσrr q̇αE
= r
dr 8k
and substituting
q̇αE 2 q̇αE 2
σθθ = r + (r − Ro2 )
8k 16k
q̇αE 2r2 r2 − ro2
 
σθθ = +
k 16 16
q̇αE 2r + r2 − ro2
 2 
σθθ =
k 16
and simplifying to find the final result for σθθ to be
q̇αE
σθθ = (3r2 − ro2 )
16k

11
Jeffrey A. Anderson ME560: Project 2

The stress state is now defined by

 
σ 0
σ = rr
0 σθθ

with maximum stresses found by conducting first derivative tests as follows

dσrr q̇αEr
=
dr 8k
dσθθ 3q̇αEr
=
dr 8k
This shows that the maximum stresses in both coordinates appear at r = 0 with
values of
q̇αE 2
σrr (0) = − r
16k o
q̇αE 2
σθθ (0) = r
16k o
While most of the values needed to evaluate these expressions are available
in the literature, the internal heat generation needs to be calculated from given
temperatures. This can be solved for starting with equation 7

q̇r2 r2
 
T (r) = o 1 − 2 + Ts
4k ro

q̇ 2
T (r) = (r − r2 ) + Ts
4k o
4k(T (0) − Ts )
q̇ =
ro2 − r2
to yield
4k(T (0) − Ts )
q̇ =
ro2
at the center temperature where r = 0.

12