CANDU fuel bundle deformation modelling with COMSOL multiphysics

CANDU fuel bundle deformation modelling with COMSOL multiphysics

Nuclear Engineering and Design 250 (2012) 134–141 Contents lists available at SciVerse ScienceDirect Nuclear Engineering and Design journal homepage...

932KB Sizes 2 Downloads 72 Views

Nuclear Engineering and Design 250 (2012) 134–141

Contents lists available at SciVerse ScienceDirect

Nuclear Engineering and Design journal homepage:

CANDU fuel bundle deformation modelling with COMSOL multiphysics J.S. Bell, B.J. Lewis ∗ Department of Chemistry and Chemical Engineering, Royal Military College of Canada, 17000 Station Forces, Kingston, Ontario, Canada K7K 7B4

h i g h l i g h t s    

The deformation behaviour of a CANDU fuel bundle was modelled. The model has been developed on a commercial finite-element platform. Pellet/sheath interaction and end-plate restraint effects were considered. The model was benchmarked against the BOW code and a variable-load experiment.

a r t i c l e

i n f o

Article history: Received 19 January 2012 Received in revised form 4 May 2012 Accepted 8 May 2012

a b s t r a c t A model to describe deformation behaviour of a CANDU 37-element bundle has been developed under the COMSOL Multiphysics finite-element platform. Beam elements were applied to the fuel elements (composed of fuel sheaths and pellets) and endplates in order to calculate the bowing behaviour of the fuel elements. This model is important to help assess bundle-deformation phenomena, which may lead to more restrictive coolant flow through the sub-channels of the horizontally oriented bundle. The bundle model was compared to the BOW code for the occurrence of a dry-out patch, and benchmarked against an out-reactor experiment with a variable load on an outer fuel element. Crown Copyright © 2012 Published by Elsevier B.V. All rights reserved.

1. Introduction The modelling of nuclear fuel behaviour involves the coupling of many inter-related complex phenomena (Van Uffelen, 2006). A number of computer codes have been developed to model fuel rod performance during normal and transient conditions (Godesar et al., 1970; Nakajima et al., 1994; Bonnaud et al., 1997; Berna et al., 1997; Cunningham et al., 2001; Lassmann, 1992; Rashid et al., 2004). Typically, separate codes are used to simulate the steady-state and transient behaviour; for instance, in Canada, the ELESTRES code has been developed for the prediction of fuel performance for normal (steady-state) conditions (Tayal, 1987; Sim et al., 2001), and the ELOCA code for upset/transient conditions (Sills, 1979; Williams, 2005). Most codes are principally restricted to a one-dimensional (1D) radial approximation or a 1.5D analysis (i.e., a radial calculation involving an axially segmented fuel rod). However, with improved computational capability, current efforts are able to better couple the multiphysics behavior of heat/mass transfer and mechanical phenomena for multi-dimensional analysis (Williamson, 2011; Newman et al., 2009; Gaston et al., 2009).

∗ Corresponding author. E-mail address: [email protected] (B.J. Lewis).

Modelling efforts have typically focused on the individual fuelrod entity and not on the entire fuel-assembly construct. For instance, Fig. 1 shows a fuel bundle containing 37-fuel elements as used in the CANDU reactor design (Canteach). As shown in Fig. 1, the bundle is comprised of concentric rings of fuel elements (i.e., an inner central element, and an inner, intermediate and outer ring of fuel elements). Coolant flow occurs around the bundle between the outer-ring elements and the pressure tube, and through the inner sub-channels of the bundle between the individual elements. The inter-element spacing is maintained by slit spacers and the bundle is separated from the pressure tube by three-planes of bearing pads. The element is hermitically sealed with the welding of endcaps at each end of the element. In addition, the end caps are attached to end-plates that maintain the bundle structure and provide further restraint to any element movement due to external forces. In particular, the deformation of fuel elements in a fuel bundle is important in order to assess ageing effects during normal operation, which may inhibit the cooling of the bundle with coolant flow through and around the sub-channels of the bundle. In addition, during anticipated operational occurrences with, for instance, lowflow conditions, which can lead to a dryout patch on the sheath, the fuel element could bow as a result of a circumferential temperature gradient around the sheath. Also under high-temperature accident conditions, the bundle could potentially sag and/or slump due to its own weight with a reduction in the structural strength of the materials at high temperature in the horizontal fuel-channel

0029-5493/$ – see front matter. Crown Copyright © 2012 Published by Elsevier B.V. All rights reserved.

J.S. Bell, B.J. Lewis / Nuclear Engineering and Design 250 (2012) 134–141


Fig. 1. 37-Element CANDU fuel bundle and fuel channel assembly in the CANDU reactor design.

configuration (Lewis et al., 2009). These latter situations could give rise to coolability restrictions with the possible occurrence of coolant bypass effects. The individual fuel elements in a CANDU bundle can deform (i.e., deflect or bow) during their residence in the fuel channel by thermally induced phenomena (Veeder and Schankula, 1974). An analytical beam approach has in fact been used to explain bowing phenomena during normal operation in the WR-1 research reactor, where element failure occurred as a result of localized over heating/coolant-flow restrictions due to element bow and contact with the channel wall (Veeder and Schankula, 1974). Element bowing was also a factor during re-fuelling operations with the “sticking” of a fuel assembly after irradiation. Element bending can result from differential thermal expansion of the fuel/sheath materials which arise from: (i) non-uniform coolant temperatures due to imperfect coolant mixing, (ii) non-uniform heat transfer between the sheath and coolant due to local coolant flow conditions, (iii) asymmetric heat production as a result of a neutron flux gradient across an element, and (iv) the interaction of the expanding fuel pellets against the outer sheath wall (Veeder and Schankula, 1974). Bowing phenomena has also been modelled numerically using a finite-element approach as a composite beam problem with the development of the BOW code (Tayal, 1989). This code also accounts for pellet/sheath interaction effects where the fuel pellets may provide some flexural rigidity to the sheath bending, as well as “gripping” of the pellet stack by the thin collapsible sheathing during in-reactor operation with a pressurized coolant. In place of a hinged (or pinned) boundary condition for the analytical treatment of Veeder and Schankula (1974) for the single-element representation, a full bundle model was explicitly considered for the BOW code in order to account for the influence of the endplates, which can deform and/or offer resistance to the bending of the individual elements. This interaction was originally modelled with a Fourier series solution with torsional spring constants assuming equilibrium conditions for the force and moments acting on an infinitesimal segment of a given endplate ring (Tayal and Choo, 1984). The endplate model was later developed as a series of interconnected straight beam segments (Xu et al., 2005). In addition, contact effects were considered between the outer elements and pressure tube, and between adjacent fuel elements, employing “gap” elements for any clearances, where an additional force was applied in the opposing direction if overlap between two bodies occurred (Xu et al., 2005). An inelastic creep model was also added for the sheath (Yu et al., 1995). The code was successfully validated against a number of analytical solutions for single-element loadings, an ANSYS bundle simulation as well as a series of outand in-reactor experiments (Xu et al., 2004, 2005; Yu et al., 1995, 1997). The current work describes the development of a full-bundle model accounting for coupled thermo-mechanical behaviour using the commercial COMSOL Multiphysics software package.

2. Model development As detailed in Section 2.1, since the fuel element is much longer than its diameter, beam theory can be used to describe element deformation behaviour for the fuel stack and sheath. The deformation model needs to be coupled to a heat conduction analysis in order to account for temperature gradients and differential expansion of various materials within the element. The individual elements can be restrained somewhat by the end-plate webbing that is modelled as a boundary condition (see Fig. 1 and Section 2.2). The thermophysical properties required for the fuel pellets and sheathing materials are summarized in Appendix A. The complete details of the model are given in (Bell, 2011). 2.1. Fuel element model 2.1.1. Deformation model The deformation of a beam can be calculated as a “weak form” of the virtual work principle within the structural mechanics module in COMSOL Multiphysics: (COSMOL, 2008)

(−εel  + uT F)dV = 0

ıW = ı



Here F (N m−3 ) is a force per unit volume vector which is known from the external loads, u (m) is the displacement vector,  (Pa) is the internal stress acting within the deformed body, εel is the elastic strain and ı is the virtual displacement. This principle describes the equilibrium state of a deformed material, where a solid has deformed due to external load forces with the energy stored within the material as a stress. In this approach, an infinitesimally small displacement field (ı) is applied with a deformation equal to ∇ ı (Ibrahimbegovic, 2009), which yields a partial differential  equation uT FdV .

with boundary conditions dependent on the value of ı V

For a beam element, a simplification is possible for εel  considering the total strain in the solid (COSMOL, 2008): ε = εi + εth + εel


Here εi and εth refer to the initial strain in the solid before the application of an external load and thermal strain, respectively. With a beam approximation, the total strain can be given by (COSMOL, 2008) ε=z

∂y ∂z ∂uaxi +y + ∂s ∂s ∂s


The parameters z and y (radians) depicted in Fig. 2 are the angle of rotation of the cross section with respect to the z- and yaxes, respectively, s (m) is the length of the body along the neutral axis and uaxi is the displacement of the cross section in the direction of the neutral axis. In this formulation, the strain depends on the


J.S. Bell, B.J. Lewis / Nuclear Engineering and Design 250 (2012) 134–141

Fig. 2. 3D meshing of the endplate and sheath geometry of a 37-element bundle modelled by beam theory. The bundle is oriented along the x-axis.

deformation vector u (where uaxi is a displacement). The thermal strain in Eq. (2) is:

εth = ˛

T + Ty

y z + Tz − Tref hy hz


where ˛ (K−1 ) is the thermal expansion coefficient of the material, T (K) is the temperature of the material, Tref is the reference temperature from which the thermal expansion is measured, T (K) is the temperature gradient across the cross section of the material in the y- and z-axes, and h (m) is the cross section of the total length in the given specified dimension. Thus, Eq. (2) can be solved for εel by substituting in expressions for ε and εth from Eqs. (3) and (4), respectively, and accounting for any initial strain εi from any external loads in a given direction. The strain terms ∂y /∂s, ∂z /∂s and ∂uaxi /∂s are dependent on the displacement at equilibrium, while the remaining terms follow from the initial conditions (εi ), known material properties (˛ and Tref ), or from a coupled heat conduction analysis in item (ii) (T and T ). The force (N) acting normal to the cross section of the body is:


 · dA



The energy stored by the body as it rotates around the neutral axis can be evaluated from the torsional moment Mx : ∂x Mx = GJ ∂s


(m4 )

Here J is the torsional constant of the cross section for the given geometry and G (Pa) is the shear modulus of the body, which depends on the geometry and material properties of the solid, i.e., for a cylindrical body (Beer and Johnson, 1981): G=


where E (Pa) is the elastic modulus of the material and  is Poisson’s ratio. As follows from the volume integral in Eq. (1), an equivalent line integral can be given for the internal stress component: ıW = ı

My L

∂y ∂u ∂z ∂x + Mz + Mx + N axi ∂s ∂s ∂s ∂s

uT FdV

dx + ı V


where the moments My and Mz are given by: Mz = E · Iz ·

∂y ∂s

∂z My = E · Iy · ∂s

MTot = Ms + Gr · Mp



If there is no pellet/sheath interaction Gr = 0. On the other hand, without the presence of a fuel-to-sheath gap as a consequence of hard contact between the fuel pellets and sheath, Gr = 1. Thus, the total bending moment in Eq. (10) provides a means to transfer stress/strain from the pellets to the sheath which, in turn, would influence the fuel element bow, i.e., with fuel/sheath contact, the pellets can add some flexural rigidity to the fuel element to resist in the bending of the sheath itself. The “height” of the cross section in both the y and z directions for the pellet and sheath is, respectively: hy = hz = 2 · rpel








= 2 · rso

where rpel and rso are the radius of the fuel pellet (=0.006081 m) and outer sheath radius (=0.006486 m), respectively, for the CANDU fuel design. In addition, this analysis also requires the moments of inertia for the fuel pellets and sheath, as well as the torsional constants for these components (Beer and Johnson, 1981): Iy = Iz =

 4 r 4 pel


Iys = Izs =

 4 (r − rsi4 ) 4 so



E 2(1 + )


Here I is the moment of inertia for the given body (see below). Hence, Eq. (8) can be expressed in terms of known constants and the deformation vector u. The overall bending of a fuel element will be affected by the interaction between the pellet and sheathing (see Fig. 3) (Veeder and Schankula, 1974). In the current treatment, independent beam geometries are considered for the cylindrical fuel-pellet stack (p) and sheath (s), which are coupled via a “gripping factor” Gr through their individual bending moments in order to account for the possibility of pellet/sheath interaction:


Jp =

 4 r 2 pel


Js =

 4 (r − rsi4 ) 2 so


where rsi is the inner radius of the sheath (=0.006082 m). The external loads acting on the bundle components are the force of gravity (i.e., with an acceleration of 9.81 m s−2 in the negative z-direction), and the hydraulic drag force from coolant flow in the axial direction. A hydraulic drag force was not included in this analysis (since the deformation model was validated principally against an out-reactor experiment in Section 3). Any additional

J.S. Bell, B.J. Lewis / Nuclear Engineering and Design 250 (2012) 134–141

Fuel-to-Sheath Gap




rsi rso

Sheath Fuel Fig. 3. Schematic of a portion of a CANDU fuel element taken from for the cross section (left figure) and axial view (right figure) (not to scale).

loads on a fuel element (such as the drag force or additional crossflows) can be included in the model forces, which would act in the corresponding direction and in the given bundle component. As mentioned, the dependent variables are the displacement and angle  with respect to the x, y, and z-axis coordinate system. If the normal of the cross section initially lies in a direction other than the global x-axis, a local coordinate system is established for the beam in order to perform a calculation around the local x-axis. The displacement is then transformed into the global reference frame. This results in a total of 12 degrees of freedom (DOF) at each node within the finite-element mesh for the beam (i.e., 6 DOF for displacement and rotation in the local x-axis, and an additional 6 DOF for the displacements and rotations that are transformed into the global x-axis) (COSMOL, 2008; Cook, 1995).

2.1.2. Heat conduction analysis A thermal solution must be coupled to the deformation model of Section 2.1.1 since bundle deformation is a thermally driven process (Veeder and Schankula, 1974). The temperature distribution in the fuel and sheath can be calculated with a solution of a time-dependent heat conduction equation (Olander, 1976):

For a time-varying element linear power history, the fuel burnup as a function of time can be evaluated from the solution of the first order differential equation (Morgan, 2007): Plin MUO2 dBu = 2 dt  · rpel UO2 MU


where MUO2 is the molecular weight of the uranium dioxide fuel (=270.03 g mol−1 ), MU is the atomic weight of natural uranium (=238.03 g mol−1 ) and UO2 is the density of the non-irradiated fuel (=10,650 kg m−3 ). To account for heat transfer across the thin gap, an overall heat transfer coefficient from the fuel to the sheath hfs (W m−2 K−1 ) was estimated with the ELESTRES fuel performance code as a function of burnup at three different element linear power ratings of 25, 40 and 55 kW m−1 , respectively (Bell, 2011): hfs = 54.9 + 0.404Bu − 0.00682Bu2 + 6.51 × 10−5 Bu3 − 3.32 × 10−7 Bu4 + 8.64 × 10−10 Bu5 − 9.07 × 10−13 Bu6 (21a) hfs = 70.8 − 0.0922Bu + 0.00156Bu2 − 3.95 × 10−6 Bu3 − 4.30 × 10−8 Bu4 + 2.71 × 10−10 Bu5 − 4.51 × 10−13 Bu6


∂T = ∇ · (k∇ T ) + Q ∂t

hfs = 78.5 + 0.347Bu − 0.0410Bu2 + 6.18 × 10−4 Bu3 (kg m−3 )

is the Here t is the time (s), T (K) is the temperature, density of the material, Cp (J kg−1 K−1 ) is the heat capacity of the material at constant pressure, and k (W m−1 K−1 ) is the thermal conductivity of the material. The heat generation profile due to fission Q (W m−3 ) at the radial position r from the centre of the pellet, accounting for flux depression effects through the pellet of radius rpel (m), is given by Olander (1976):

Q =

Plin 2 rpel


2I1 ( rpel )

I0 ( r)


where Plin (W m−1 ) is the linear power rating for the element,

(m−1 ) is the inverse diffusion length through the pellet, and I0 and I1 are the zeroth and first order modified Bessel function. The diffusion length depends on the fuel burnup Bu (MW h kg−1 U) (Morgan, 2007):

= 89.89 + 0.60Bu0.68




− 4.10 × 10−6 Bu4 + 1.29 × 10−8 Bu5 − 1.57 × 10−11 Bu6 (21c) This heat transfer coefficient therefore takes into account fission gas release into the gap as a function of burnup. Hence, assuming sheath collapse with an effective gap size of tg ∼ 1 ␮m, a thermal conductivity for the fuel-to-sheath gap of kg = hfs ·tg can be estimated. Consequently, steady-state heat transfer across the gap can be evaluated by assuming simple heat conduction in accordance with Eq. (17) (with no heat source). Given the thermophysical properties in Appendix A, with no heat source term in the sheath, Eq. (17) can be solved for both the fuel pellets and sheath. These solutions require an outer sheath temperature Tso (K) as a boundary condition: Tso = Tc +

Plin 2 · rso · hsc


where Tc (K) is the temperature of the coolant (∼553 K for normal operating conditions) and hsc is the heat transfer coefficient between the sheath and coolant (=50,000 W m−2 K−1 ).


J.S. Bell, B.J. Lewis / Nuclear Engineering and Design 250 (2012) 134–141

IzR3 =

1 3 th hR3 3 end


IyR3 =

1 th h3 3 end R3


J R3 =

1 th hR3 (h2R3 + th2end ) 12 end


Here the superscripts associated with the second moments of inertia and the torsional constants refer to the inner ring 1 (R1), the intermediate ring 2 (R2), and the outer ring 3 (R3), along with the end plate web connectors between the rings (W) (see Fig. 2). The variable hw12 is the height (or width) of the web connector between rings 1 and 2 (=3.4 mm), thend is the thickness of the endplate (=2.4 mm) and hR3 is the height of ring 3 (=4.5 mm) (Bell, 2011). Fig. 4. Meshed cross section of a fuel element for the heat conduction solution.

The 2D meshing for the heat conduction model is given in Fig. 4. The thin gap was meshed with 100 quadratic quadrilateral segments for a single fuel element cross section, and the pellets and sheath contained a total of 1020 quadratic triangular and quadrilateral elements. 2.2. Endplate model The endplates are relatively thin with a thickness of just a few millimeters (2.4 mm). Therefore, the deformation of all endplate components could be approximated by a series of beam segments (Tayal and Choo, 1984). Hence, the deformation of the endplates can be calculated using the same beam equations as from Section 2.1. By assuming that the endcap-to-endplate weld acts as a continuous solid, the deformation of the endplates can therefore be used as boundary conditions for the ends of the elements within the bundle (Fig. 2). The endplates are made of the same material as the fuel sheath. The second moments of inertia and torsional constants for each ring and web section of the endplate is given by Bell (2011): IyR1 = IyR2 = IyW =

1 th h3 3 end w12


IzR1 = IzR2 = IzW =

1 3 th hW 12 3 end


J R1 = J R2 = J W =

1 th hw12 (h2w12 + th2end ) 12 end


3. COMSOL bundle model validation The temperature profile of the fuel elements is further required for a complete analysis of the coupled thermal expansion/fuel deformation behaviour. Material temperatures are represented by average temperatures evaluated with the heat conduction model as described in Section 2.1. This calculation can be compared to temperature predictions with the ELESTRES-IST industry-standard fuel performance code. Fig. 5(a) and (b) shows a comparison of the average fuel and sheath temperatures for the two codes at three representative linear fuel ratings of 25, 40 and 55 kW m−1 , which cover the typical range of CANDU fuel operation. The average fuel temperatures in Fig. 5(a) with the COMSOL model are in reasonable agreement (i.e., within ∼60 ◦ C) of the ELESTRES code calculations. The discrepancy between the two codes increases with linear power and burnup since the COMSOL model does not consider irradiation phenomena associated with fission gas generation and fuel structural changes (i.e., fuel densification/swelling and grain growth). In particular, the COMSOL model only includes heat conduction with no feedback effects of irradiation on the fuel thermal conductivity, although a correlation is implemented to capture degradation effects of fission gas release on the heat transfer coefficient in the fuel-to-sheath gap (i.e., Eq. (21)). The average temperature in the fuel sheath in Fig. 5(b) is predicted within ∼1% for all cases. The COMSOL bundle deformation model was also benchmarked against analytical formulations for several simple-beam loadings, an out-reactor fuel-bundle deformation experiment and a dryoutpatch simulation with the BOW code.

Fig. 5. Comparison of the average (a) fuel and (b) sheath temperature calculated with the COMSOL and ELESTRES codes for three element linear powers as a function of burnup.

J.S. Bell, B.J. Lewis / Nuclear Engineering and Design 250 (2012) 134–141


Fig. 6. Comparison of a COMSOL simulation of a fuel bundle with a variable loading on a single outer element without (Gr = 0) and with (Gr = 0.5) pellet/sheath interaction to experimental results.

The analytical formulations considered for the BOW code validation, included: (i) a single fuel-element geometry with one-end of the element fixed, and 10 and 6 Nm moments applied horizontally and vertically along the element, respectively, and (ii) a single-element geometry simply supported at both ends, with a temperature gradient of 30 and 60 K applied horizontally and vertically along the element, respectively. The COMSOL simulations are typically within 2% of the analytical evaluations and BOW code simulations for these cases (Bell, 2011). The complete COMSOL bundle model (including endplate interaction) was further benchmarked against an out-of-reactor experiment for bundle deformation as conducted by Atomic Energy of Canada Limited (Yu et al., 1995). The experiment consisted of taking an unirradiated CANDU 6 fuel bundle at room temperature, and displacing an outer-ring fuel element away from the fuel bundle at the mid-plane location, where an applied force was varied up to 45 N. Fig. 6 compares the results of the simulation to experiment. The model predictions are in excellent agreement with the measurements when no pellet/sheath interaction is considered (i.e., Gr = 0). This result is to be expected since there is no sheath creep down in the out-reactor experiment, i.e., creepdown of the sheath onto the pellets (and hence fuel/sheath contact) would only occur in-reactor in the presence of a pressurized coolant. The sensitivity of the model to the effect of pellet/sheath gripping is further shown in Fig. 6 with the inclusion of a simulation showing partial pellet/sheath interaction (i.e., Gr = 0.5). Fig. 7 also shows the results of a COMSOL simulation of a dryout patch on a single fuel sheath compared to an original simulation with the BOW code that was used to demonstrate its capability (Tayal, 1989; Bell, 2011). In both sets of results, the maximum displacement occurs roughly at the center of the fuel element. The BOW code simulation indicates a maximum displacement of 0.33 mm compared to a prediction of 0.39 mm for the COMSOL analysis. Both codes reveal an asymmetric displacement caused by the presence of a second dry-out patch. Aside from the differences in the thermal expansion coefficients used for the two models, the BOW model also includes a torsional spring treatment to account for the interaction of the elements with the endplates, whereas a simple “hinged boundary” condition is applied in the COMSOL simulation to account for the endplate interaction. The use of this simpler boundary condition in the COMSOL simulation does not provide any additional resistance for restriction of the ends of the fuel element to rotate. Hence, the use of this different boundary condition would account for the slightly greater displacement caused by the second drypatch in the COMOL simulation.

Fig. 7. Comparison of the predicted element bowing with the (a) BOW code and (b) COMSOL model for the simulation of a dry-out patch with a circumferential sheath temperature gradient.

4. Conclusions 1. The current analysis demonstrates that bundle deformation phenomena can be simulated on the commercial COMSOL Multiphysics platform, in which the elastic stresses and strains produced within the fuel elements can be calculated with a beam approximation. The individual bending moments of the pellets and sheath are coupled via a phenomenological “gripping” factor to account for pellet/sheath interaction effects, although this empirical approach could be replaced with a more advanced model with physically based contact or with a solidelement mechanics model. Short straight beam approximations were also considered for the endplate webbing to model the restraint offered by the endplates. The current model conservatively neglects for any restraining effects offered by the element appendages in the bundle design (i.e., element-to-element spacing as maintained with split spacers and element-to-pressure tube restraint by several planes of bearing pads), which will become important with larger deformations of the elements. 2. The heat conduction solution for a fuel element in the COMSOL model was successfully benchmarked against the ELESTRES-IST fuel performance code for several power ratings. The bundle model was also successfully compared against analytical formulae of stress and strain for the deformation of a simple beam, and an out-reactor experiment, where the deflection of an outer ring fuel element in a fuel bundle was measured with a load applied to the mid-plane of the element. The COMSOL model was also successfully benchmarked against a BOW code simulation of a dry-out patch analysis with the presence of a circumferential temperature gradient.


J.S. Bell, B.J. Lewis / Nuclear Engineering and Design 250 (2012) 134–141


The parameter Cp in Eq. (17) for UO2 is given as a function of temperature (Fink, 2000):

The authors would like to acknowledge funding from the CANDU Owners Group (COG), University Network of Excellence in Nuclear Engineering (UNENE) and the Natural Sciences and Engineering Research Council of Canada (NSERC).

Cp,UO2 = 1.93238 × 102 + 3.25744 × 10−1 T − 3.12004 × 10−4 T 2 + 1.16822 × 10−7 T 3 − 9.75333 × 10−12 T 4 − 2.644 × 106 T −2

Appendix A. Thermophysical properties of the fuel and sheath

The density of uranium dioxide UO2 can be evaluated from:

The thermal expansion coefficient for uranium dioxide as a function of temperature T (K) is given by Fink (2000) and the corresponding correlation for Zircaloy-4 is obtained from the materials library, MATPRO (Hagerman, 1997):

˛UO2 =


UO2 = 0.316


zirc = 0.30


The thermal conductivity kUO2 (W m−1 K−1 ) in Eq. (17) for the pellet is given by the MATPRO relation (Hagerman, 1997):

+ 5.2997 × 10

D 1 + (1 − D)6.5 − 0.00469T  )


T ×e


Cv (A + BT ")






(A.4) The functions T and T are given by:

T  =


T < 1834 K


T ≥ 1834 K

⎧ ⎨ T, T < 1800 K ⎩

0.5T + 900, 2050,

1800 K ≤ T ≤ 2300 K

T > 2300 K

and D is the fraction of theoretical density. The parameter A takes into account effects of point defects in the lattice on the mean free path of a phonon, and B considers phonon self-scattering effects. Finally, εth is the thermal strain due to thermal expansion of the fuel. Since the fuel is initially stoichiometric, with no additional plutonium content, the parameters A and B are taken as constants equal to 0.339 m s kg−1 K−1 and 0.06867 m s kg−1 K−1 , respectively. The parameter Cv ( J kg−1 K−1 ) in Eq. (A.4) is the specific heat of uranium dioxide at constant volume, which accounts for heat transfer through lattice vibrations (Hagerman, 1997), Cv =

297.7(535.285)2 T 2 [e(535.285/T ) − 1]


[e535.285/T ]




where D is the fraction of theoretical density of the uranium dioxide. The Poisson’s ratio for these materials is further given by Morgan (2007) and Hagerman (1997):

1 (1 + 3εth )



for 923 K ≤ T ≤ 3120 K

Ezirc = 9.2 × 1010 − 4.05 × 107 T

T =


(0.9672 + 1.179 × 10−5 T − 2.492 × 10−9 T 2 + 1.219 × 10−12 T 3 ),

EUO2 = 2.334 × 1011 {1 − 2.752(1 − D)}{1 − 1.0915 × 10−4 T } (A.2a)

kUO2 =

th =


for 273 K ≤ T ≤ 923 K

The Young’s modulus for the fuel and sheath materials is given by Hagerman (1997):

UO2 = th D

(0.9973 + 9.082 × 10−6 T − 2.705 × 10−10 T 2 + 4.391 × 10−13 T 3 ),

˛zirc = 4.44 × 10−6 T − 1.485 × 10−3



where th,273 is the theoretical density at 273 K and ˛UO2 is determined from Eq. (A.1a). The thermal conductivity for the Zircaloy-4 sheath is given by the MATPRO correlation as a function of the T (K) (Hagerman, 1997): kzirc = 7.51 + 2.09 × 10−2 T − 1.45 × 10−5 T 2 + 7.67 × 10−9 T 3 (A.8) The heat capacity, Cp (J kg−1 K−1 ), for Zircaloy-4 in the temperature range between 400 and 1090 K is given by Hagerman (1997): Cp,zirc = 244.946 + 0.156T − 3.341 × 10−5 T 2


The density of the Zircaloy sheathing zirc can be taken as 6551 kg m−3 . References Beer, F.P., Johnson, E.R., 1981. Mechanics of Materials. McGraw-Hill, Inc, New York, NY, USA. Bell, J.S., 2011. Assessment of COMSOL Multiphysics for the Use of CANDU Fuel Bundle Deformation Modeling. MASc Thesis, Royal Military College of Canada. Berna, G.A., Beyer, C.E., Davis, K.L., Lanning, D.D., 1997 December. FRAPCON-3: A Computer Code for the Calculation of Steady-state, Thermal–Mechanical Behaviour of Oxide Fuel Rods for High Burnup, NUREG/CR-6534, PNNL-11513. Bonnaud, E., Bernard, C., Van Schel, E., 1997. COPERNIC: a state-of-the-art fuel rod performance code. Trans. Am. Nucl. Soc. 77 (December). Cook, R.D., 1995. Finite Element Modeling for Stress Analysis. John Wiley & Sons Ltd, Hoboken, NJ, USA. COMSOL Multiphysics Structural Mechanics Module User’s Guide, Version 3.5, COMSOL AB, 2008. Cunningham, M.E., Beyer, C.E., Medvedev, P.G., Berna, G.A., 2001. FRAPTRAN: A Computer Code for the Transient Analysis of Oxide Fuel Rods, vol. 1, NUREG/CR-6739, PNNL-13576. Fink, J.K., 2000. Thermophysical properties of uranium dioxide. J. Nucl. Mater. 279, 1–18. Gaston, D., Newman, C., Hansen, G., Lebrun-Grandie., D., 2009. MOOSE: a parallel computational framework for coupled systems of nonlinear equations. Nucl. Eng. Des. 239, 1768–1778. Godesar, R., Guyette, M., Hoppe, N., 1970. COMETHE II: a computer code for predicting the mechanical and thermal behavior of a fuel pin. Nucl. Appl. Technol. 9, 205–217 (Aug 1970). Hagerman, D.T. (Ed.), 1997. MATPRO-A Library of Material Properties for LightWater-Reactor Accident Analysis, SCDAP/RELAP5/Mod 3.2 Code manual, NUREG/CR-6150, vol. 4. Ibrahimbegovic, A., 2009. Nonlinear Solid Mechanics. Springer Science and Business Media, New York. Lassmann, K., 1992. TRANSURANUS: a fuel rod analysis code ready for use. J. Nucl. Mater. 188, 295–302. Lewis, B.J., Iglesias, F.C., Dickson, R.S., Williams, A., 2009. Overview of high temperature fuel behaviour with relevance to CANDU fuel. J. Nucl. Mater. 394, 67–86. Morgan, D., 2007. A Thermomechanical Model of CANDU Fuel, MASc Thesis, Royal Military College of Canada. Nakajima, T., Saito, H., Osaka, T., 1994. FEMAXI-IV: a computer code for the analysis of thermal and mechanical behavior of light water reactor fuel rods. Nucl. Eng. Des. 148, 41–52.

J.S. Bell, B.J. Lewis / Nuclear Engineering and Design 250 (2012) 134–141 Newman, C., Hansen, G., Gaston, D., 2009. Three dimensional coupled simulation of thermomechanics heat, and oxygen diffusion in UO2 nuclear fuel rods. J. Nucl. Mater. 392, 6–15. Olander, D.R., 1976. Fundamental Aspects of Nuclear Reactor Fuel Elements, TID26711-p1, U.S. Department of Energy. Rashid, Y., Dunham, R., Montgomery, R., 2004. Fuel analysis and licensing code: FALCON MOD01, Technical Report EPRI 1011308, Electric Power Research Institute, December. Sills, H.E., 1979. ELOCA Fuel Element Behaviour during High-Temperature Transients, Atomic Energy of Canada Limited Report AECL-6357. Sim, K.S., Chassie, G.G., Xu, Z., Tayal, M., Westbye, C., 2001. Progress in qualifying ELESTRES-IST 1.0 Code: verification and interim results of validation. In: 7th CNS International Conference on CANDU Fuel, Kingston, Ontario, Canada, 2001 September 23–27, pp. 5B 21–33. Tayal, M., 1987. Modelling CANDU Fuel under Normal Operating Conditions, ELESTRES Code Description, Atomic Energy of Canada Limited Report AECL9331, February. Tayal, M., 1989. Modeling the bending/bowing of composite beams such as nuclear fuel: the BOW code. Nucl. Eng. Des. 116, 149–159. Tayal, M., Choo, C.K., 1984. Fatigue analysis of CANDU nuclear fuel subjected to flow-induced vibrations, Atomic Energy of Canada Limited Report AECL-8331. In: Proceedings of the 25th Structures, Structural Dynamics, and Materials Conference, American Society of Mechanical Engineers, Palm Springs, California, USA, May 14–16.


Van Uffelen, P., 2006. Modelling of Nuclear Fuel Behaviour. European Commission Directorate-General Joint Research Centre Institute for Transuranium Elements, EUR 22321 EN, ISSN 1018-5593. Veeder, J., Schankula, M.H., 1974. Bowing of pelletized fuel elements: theory and in-reactor experiments. Nucl. Eng. Des. 29, 167–179. Williams, A.F., 2005. The ELOCA fuel modelling code: past, present and future. In: 9th International CNS Conference on CANDU Fuel, Belleville, Ontario, Canada, September 18–21. Williamson, R.L., 2011. Enhancing the ABAQUS thermomechanics code to simulate multipellet steady and transient LWR fuel rod behavior. J. Nucl. Mater. 415, 74–83. Xu, S.G., Yu, S.D., Tayal, M., Xu, Z., 2004. BOW: Comparisons of Creep Calculations with Analytical Solutions. Presented at the 6th CNS International Conference on Simulation Methods in Nuclear Engineering, Montreal, Canada, October 12–15. Xu, S.G., Xu, Z., Yu, S.D., Tayal, M., Lai, L., Wong, B., 2005. BOW code developmentmodeling of in-reactor bundle deformation. Presented at the 9th International Conference on CANDU Fuel, Belleville, Canada, September 19–21. Yu, S.D., Tayal, M., Singh, P.N., 1995. Improvements verifications and validations of the BOW code. Presented at the 4th International Conference on CANDU Fuel, Pembroke, Canada, October 1–4, pp. 4B.52–4B.69. Yu, S.D., Tayal, M., Xu, Z., 1997. Creep Bowing in CANDU Fuel: Modelling and Applications. Presented at the 5th International Conference on CANDU Fuel, Toronto, Canada, September 21–25, pp. 364–375.