Governing Equations

The Bison governing relations consist of fully-coupled partial differential equations for energy, species, and momentum conservation. The energy balance is given in terms of the heat conduction equation (1) where , and are the temperature, density and specific heat, respectively, is the energy released in a single fission event, and is the volumetric fission rate. can be prescribed as a function of time and space, or input from a separate neutronics calculation. The heat flux is given as (2) where denotes the thermal conductivity of the material.

Species conservation is given by (3) where , , and are the concentration, radioactive decay constant, and source rate of a given species, respectively. The mass flux is specified as (4) where is the diffusion coefficient; this definition has been used to simulate fission product transport within the fuel. Also implemented in Bison is a hyperstoichiometric model for oxygen diffusion in fuel as described in Newman et al. (2009). In this case denotes the oxygen flux in the hyperstoichiometric regime with, (5) where is diffusivity, is the thermodynamic factor of oxygen, is the heat of transport of oxygen, and is the universal gas constant.

Momentum conservation is prescribed assuming static equilibrium at each time increment using Cauchy's equation, (6) where is the Cauchy stress tensor and is the body force per unit mass (e.g. gravity). The displacement field , which is the primary solution variable, is connected to the stress field via the strain, through a constitutive relation.

Element Kinematics

For geometrically linear analysis, the strain is defined as . Furthermore, with a linear elastic constitutive model, the stress is simply . We now outline our approach for nonlinear analysis. We follow the approach in Rashid (1993) and the software package Key (2011).

We begin with a complete set of data for step and seek the displacements and stresses at step . We first compute an incremental deformation gradient, (7) With , we next compute a strain increment that represents the rotation-free deformation from the configuration at to the configuration at . Following Rashid (1993), we seek the stretching rate : (8) Here, is the incremental stretch tensor, and is the incremental Green deformation tensor. Through a Taylor series expansion, this can be determined in a straightforward, efficient manner. is passed to the constitutive model as an input for computing at .

The next step is computing the incremental rotation, where . Like for , an efficient algorithm exists for computing . It is also possible to compute these quantities using an eigenvalue/eigenvector routine. With and , we rotate the stress to the current configuration.

Axisymmetric Equations

For the axisymmetric case (RZ), the nonlinear strains are derived starting from the Green-Lagrange strain: (9) This leads to: (10) We can recover the linear strain by ignoring the higher-order terms.

Spherically Symmetric Equations

For the spherically symmetric case, the nonlinear strains are derived starting from the Green-Lagrange strain: (11) This leads to: (12) We can recover the linear strain by ignoring the higher-order terms.


For elastic behavior, a hypoelastic formulation is used, specifically, (13) where is the elasticity tensor. For isotropic elasticity, this becomes (14) with as Lame's first parameter and as the shear modulus. This stress update occurs in the configuration at . Thus as a final step, the stress must be rotated to the configuration at .

Nonlinear Materials

Fuel materials often exhibit nonlinear mechanical behavior. As a first step to modeling this behavior, von Mises linear isotropic strain hardening via an implicit radial return method was implemented in Bison. A summary of this implementation is described in the following steps:

First, an elastic trial stress is calculated using the previous stress state and a total strain increment (15) where is the linear isotropic elasticity tensor, is the total strain increment tensor, and is the stress from the previous time step. Second, a yield function is evaluated (16) where is the yield function, is the effective trial stress based on the deviatoric trial stress, is the hardening variable, and is the yield stress. If the yield function is greater than zero, then permanent deformation has occurred and a the plastic strain increment must be calculated. Otherwise, the trial stress is the new stress.

Third, the hardening variable, , and the plastic strain increment are solved via Newton iteration. (17) In the third step, is the hardening variable from the previous time step, is the hardening constant, which defines the slope of the linear strain hardening section of the stress vs. strain plot, and are the plastic strain increment for the current and previous time steps respectively, and is the shear modulus. In this Newton iteration, the residual is driven to some predefined small number as the hardening variable and the plastic strain increment are updated to achieve such a small residual.

Fourth, when the residual is sufficiently small, the new plastic strain increment is used to update a plastic strain increment tensor () that is used to calculate an elastic strain () from the total strain () and a new stress increment is calculated using this new elastic strain. (18)

Fifth, now, the stress and the plastic strain are updated (19) (20)

Sixth, in conventional nonlinear solvers, the material Jacobian is calculated and used to solve the nonlinear problem. Note however, that the material Jacobian is NOT required using the JFNK method. It can be used as a preconditioner and it is therefore presented here. (21) where (22) and (23) The source used for guidance in implementing this plasticity model into Bison was Dunne and Petrinic (2005).


  1. Fionn Dunne and Nik Petrinic. Introduction to Computational Plasticity. Oxford University Press, Oxford, 2005.[BibTeX]
  2. Sam Key. Fma-3d theoretical manual. Technical Report Ver 32, FMA Development, LLC, Great Falls, Montanta, 2011.[BibTeX]
  3. C. Newman, G. Hansen, and D. Gaston. Three dimensional coupled simulation of thermomechanics, heat, and oxygen diffusion in $\mathrm UO_2$ nuclear fuel rods. Journal of Nuclear Materials, 392:6–15, 2009.[BibTeX]
  4. M. M. Rashid. Incremental kinematics for finite element applications. International Journal for Numerical Methods in Engineering, 36:3937–3956, 1993.[BibTeX]