1D Generic PWR Tutorial
Note
This case can be found in the Cases/testCases/generic_pwr_1D/ folder of the OFFBEAT repository.
This tutorial serves as a demonstration on how to setup an OFFBEAT case for the 1.5D simulation of a generic PWR fuel rod. The tutorial case consist in a one-year long irradiation of standard UO2 fuel with the following characteristics:
The axially averaged power prescribed to the rod as a function of time is rather simple, and consists of an initial 1-hour ramp to power that is then kept constant throughout the irradiation, as depicted in Figure 1.
Linear heat generation rate is prescribed with an axial profile that varies over time. For visualization purposes, Figure 2 shows the normalized axial power distribution at the simulation start.
Mesh
OFFBEAT provides a Python script to facilitate geometry creation and meshing that can be found at ./offbeat/tools/rodMaker. In the input file named rodDict, users can specify the parameters for the geometry. For this case, the mesh consists of a single block for the fuel and two vertically stacked blocks that forms the cladding: the first as tall as the fuel column, and the second constituting the portion of the cladding that bounds the upper plenum volume. The mesh used for this tutorial is shown in Figure 3, reduced 100x in the z-direction for visualization.
Tip
The cladding mesh is here divided in two blocks to facilitate the construction of a computational grid with conformal faces between fuel outer and cladding inner surface, which is optimal for the mapping algorithm of OFFBEAT.
Note
The plena volumes are automatically computed by OFFEBAT from the geometry of the computational domain. The user can however specify an additional reservoir volume without explicitely modeling it (see gap gas model options).
Warning
The axial discretization of the mesh determines the maximum amount of axial slices for the simulation, as they cannot be larger than the number of axial divisions of the computational domain.
Setting the Boundary Conditions
According to the OpenFOAM standard, initial and boundary conditions are set in the 0/ folder for each of the fields solved by the code. For this case, three sub-physics solvers are enabled:
-
A small-strain mechanical sub-solver, for which
Dboundary conditions are needed, -
A solid diffusion thermal sub-solver, for which
Tboundary conditions are needed, -
A neutron diffusion sub-solver, for which
neutronFlux0boundary conditions are needed.
According to the 1.5D approach, all the governing equations are solved radially (i.e. along x in this case) for each of the axial slices. To disable the solution of the governing equations along y and z:
-
Boundary conditions for surfaces with a normal aligned along z are set of type
empty(or not specified at all), -
Boundary conditions for surfaces with a normal perpendicular to the radial direction are set of type
wedge, anOpenFOAMboundary condition that imposes axi-symmetry.
Warning
The surfaces to which empty (or wedge) boundary condition is imposed need to be set as an empty (or wedge) patch type at the stage of building the mesh. This can be done by modifying the entry type in the sub-dictionary of the boundary list in the constant/blockMeshDict dictionary.
Displacement BCs
Displacement boundary conditions are specified in the 0/D file.
Contact boundary conditions are imposed to the fuel outer surface and cladding inner surface according to the following dictionary:
"fuelOuter|cladInner"
{
type gapContact;
patchType regionCoupledOFFBEAT;
penaltyFactor 0.1;//(1)
frictionCoefficient 0; //(2)
relax 1.0;
relaxInterfacePressure 0.1; //(3)
value $internalField;
}
A penalty factor lower than 1 allows slight inter-penetration of bodies but facilitates convergence.- Friction interaction between fuel and cladding is here neglected.
- Contact pressure between fuel and cladding is under-relaxed to facilitate convergence.
On the cladding outer surface, the effect of coolant is modeled by imposing the coolantPressure boundary condition. This allows the user to specify a time-dependent list of coolant pressures:
cladOuter
{
type coolantPressure;
coolantPressureList
{
file "$FOAM_CASE/constant/systemPressure";//(1)
outOfBounds clamp;//(2)
};
relax 1.0;
value $internalField;
}
- In this case the time-dependend pressure is read from a file.
- If simulation time goes beyong the maximum time specified in the time-dependent list, the value of pressure is kept constant to the last available one in the list.
Temperature BCs
Thermal boundary conditions are specified in the 0/T file.
The thermal coupling between the fuel outer surface and cladding inner surface is enforced with the fuelRodGap boundary condition:
fuelOuter
{
type fuelRodGap;
patchType regionCoupledOFFBEAT;
kappa k;
coupled true;
roughness uniform 2.2e-6; //(1)
value $internalField;
relax 1;
}
cladInner
{
type fuelRodGap;
patchType regionCoupledOFFBEAT;
kappa k;
coupled true;
roughness uniform 0.5e-6;
value $internalField;
relax 1;
}
- The roughness of the surface plays a role in the computation of the solid conduction thermal resistance when the gap closes.
On the cladding outer surface, temperature is prescribed with a time-dependent Dirichlet boundary condition, called uniformFixedValue:
Neutron Flux BCs
Neutron flux boundary conditions are specified in the 0/neutronFlux0 file.
The solution of a diffusion-based neutron transport equations is done in OFFBEAT to feed the burnup model with the right radial flux distribution, that in turn allows to determine the radial power distribution. What is therefore relevant in OFFBEAT is the calculation of the flux shape, rather than determining it in terms of absolute values.
Therefore, the strategy commonly used in OFFBEAT calculations to have a well posed mathematical problem is to set Dirichlet boundary conditions on the fuel outer surface (e.g. = 1) and to the cladding surfaces( = 0, as flux will not contribute to power production within the cladding).
"cladInner|cladOuter"
{
type fixedValue;
value uniform 0;
}
fuelOuter
{
type fixedValue;
value uniform 1;
}
Setting the solverDict
The constant/solverDict dictionary is the main input dictionary of OFFBEAT, where physics sub-solvers, materials and models are set by the user.
Physics Sub-Solvers
The user can enable or disable the solution of the different physics from the constant/solverDict file:
thermalSolver solidConduction;
mechanicsSolver smallStrain;
neutronicsSolver diffusion;
elementTransport fromLatestTime; //(1)
- Transport of porosity/actinides is switched off here, as not relevant for LWR fuels.
Warning
Enabling the neutronicsSolver is necessary if:
- The
burnupLassmanntype is enabled for burnup; - The
fromBurnuptype is selected under theheatSourceOptions\radialProfilesub-dictionary.
A dedicated sub-dictionary is available for each of the sub-physics solvers, where the user can input additional parameters/settings. For this case for example, the multi-material correction is enabled from the mechanical solver settings as shown hereafter:
mechanicsSolverOptions
{
forceSummary off; //(1)
cylindricalStress on; //(2)
multiMaterialCorrection
{
type uniform;
defaultWeights 1;
defaultWeightsGrad 0;
}
}
- If enabled, it prints in terminal the force balance across all the boundaries of the domain.
- Enable the writing to disk of stress tensor in cylindrical coordinates.
Other Fuel Behavior Models
Besides enabling the specific sub-solvers, OFFBEAT also requires the user to select additional fuel performance models in the constant/solverDict dictionary, as shown in the following code snippet:
materialProperties byZone; //(1)
rheology byMaterial; //(2)
burnup Lassmann; //(3)
heatSource timeDependentLhgr; //(4)
fastFlux timeDependentAxialProfile; //(5)
gapGas FRAPCON; //(6)
fgr SCIANTIX; //(7)
sliceMapper autoAxialSlices; //(8)
corrosion fromLatestTime; //(9)
- Material models are assigned for each mesh
cellZone. This is an OpenFOAM entity that simply groups ensembles of cells, and it is assigned in thesystem/blockMeshDict. - Rheology models (e.g., elastic, visco-plastic material laws) are assigned for each material.
-
Lassmannmodel is chosen for nuclides depletion and for computing the fuel radial power density profile. - This heat source model allows to define a time-dependent Linear Heat Generation Rate, whose settings are defined in the
heatSourceOptionssub-dictionary. - This fast flux model allows to define a time-dependent Fast Flux, whose settings are defined in the
fastFluxOptionssub-dictionary. - This model takes care of tracking gas pressure and composition along irradiation by monitoring changes in the rod inner free volume.
- Fission gas behavior is simulated through internal online coupling with the
SCIANTIXcode. - Axial slices are needed by some OFFBEAT models (e.g.,
burnup,relocation).autoAxialSlicesis thought for use with structured meshes, building an axial slice per axial cell layer. - Corrosion model is disabled here.
For most of these models sub-dictionaries are available to specify additional options. For example, the rheologyOptions sub-dictionary allows to specify important settings for 1.5D simulations, reported hereafter:
rheologyOptions
{
thermalExpansion on; //(1)
modifiedPlaneStrain on; //(2)
springModulus 3500; //(3)
coolantPressureList
{
file "$FOAM_CASE/constant/systemPressure"; //(4)
outOfBounds clamp;
}
planeStress off; //(5)
}
- Switch to enable/disable thermal deformation.
- If activated, the strain tensor z-component is modified to consider the typical modified-plane strain approximation of standard 1.5D fuel performance calculations.
- Modulus in N/m for the plenum spring.
- Time-dependent list for coolant pressure that acts on the top cladding cap.
- If activated, the strain tensor z-component is modified to consider the plane stress approximation.
To accurately compute the rod inner free volume, the code needs to be informed by the user about the names of fuel bottom and top surfaces, together with names of inner (if any), outer fuel and inner cladding surfaces. This is specified in the gapGasOptions dictionary:
gapGasOptions
{
gapPatches ( fuelOuter cladInner );
holePatches (); //(3)
topFuelPatches ( fuelTop);
bottomFuelPatches ( fuelBottom);
gapVolumeOffset 0.0; //(1)
gasReserveVolume 0.0; //(2)
gasReserveTemperature 290;
}
- It is possible to define an additional gap volume, i.e. an additional volume that will have the same temperature and pressure conditions of the ones computed for the rod gap.
- The user can also specify an additional reserve gas volume with a constant temperature.
- In case of annular fuel pellets, also the surfaces that bound the hole must be specified.
Settings related to the imposed power linear density are collected in the heatSourceOptions dictionary, where it is possible to set the axial and radial profiles:
heatSourceOptions
{
timePoints ( 0 3600 1.26E+08 );
lhgr ( 0 15000 15000 ); //(1)
timeInterpolationMethod linear; //(2)
axialProfile
{
type timeDependentTabulated; //(3)
#include "axialProfile";
axialInterpolationMethod linear;
burnupInterpolationMethod linear;
}
radialProfile
{
type fromBurnup;
}
materials ( fuel );
}
- This is the axially averaged linear power variation over time.
- The power is linearly interpolated for the simulated times that fall between the defined
timePoints. - This option allows to define a time-varying axial profile.
The timeDependentTabulated option requires the user to define:
timePoints: times at which the axial power profile is defined. They can differ from thetimePointsdefined in the generalheatSourceOptionsdictionary.axialLocations: normalized axial coordinates at which the profile is defined, they must cover the inteval [0,1].data: a list of profile values, i.e. one list per time point. They are multiplicators for time-dependent values defined inlhgr.
For this case, these data are collected in the axialProfile file:
timePoints ( 0.00E+00 1.58E+07 ... 1.26E+08 );
axialLocations ( 0 0.0909 0.1818 ... 1 );
data
(
( 0.4697 0.9138 ... 0.2714 )
( 0.6259 1.0305 ... 0.4570 )
( 0.8070 1.1132 ... 0.6295 )
);
Warning
The profile values defined in data must be normalized by their integral average.
The settings for SCIANTIX, the fission gas behavior code coupled to OFFBEAT, are found in the solverDict as well:
fgrOptions
{
nFrequency 1; //(1)
relax 1; //(2)
addToRegistry on; //(3)
SCIANTIX
{
iverification 0; // (0= no verification)
igrain_growth 1; // (1 = ainscough)
iinert_gas_behavior 1; // (1= do IGB)
igas_diffusion_coefficient 1; // (1= Turnbull et al., 1988)
iintra_bubble_evolution 1; // (1=Pizzocri et al., 2018)
ibubble_radius 1; // (1= Olander&Wongy, 2006)
iresolution_rate 1; // (1=Turnbull 1971)
itrapping_rate 1; // (1= Ham 1958)
inucleation_rate 1; // (1= Baker 1971)
isolver 1; // (1= SDA, Pizzocri et al., 2019)
iformat_output 1; // (1 = output.txt, values separated by tabs)
igrain_boundary_vacancy_diffusion_coefficient 1; // (1= Reynolds and Burton, 1979)
igrain_boundary_behaviour 1; // (1= do InterGranularGasBehavior - Pastore et al., 2013; Barani et al., 2017)
igrain_boundary_micro_cracking 1; // (1 = Barani et al., 2017)
igrain_recrystallization 0; // (0 = non active)
ifuel_reactor_type 0; // (0=UO2/PWR)
igas_effective_coefficient 1; // gas effective (=0) or single atom (=1)
igas_sweeping 1; // add boundary gas sweeping
}
}
- SCIANTIX is run after
nFrequencyouter iterations are completed. - Under relaxation factor for inter- and intra-granular swelling computation.
- Switch to enable the storage on object registry of SCIANTIX fields.
Materials Definition
Materials are defined in OFFBEAT in the dedicate solverDict subdictionary called materials:
materials
{
fuel // Name should match the corresponding cellZone name in the mesh
{
material UO2;
...
}
cladding // Name should match the corresponding cellZone name in the mesh
{
material zircaloy;
...
}
}
A different sub-dictionary for each of the different cellZones is required by the code. The sub-dictionary name should match the cellZone name (fuel and cladding in this case). The material model is then selected for each material, in this case UO2 for fuel and zircaloy as cladding. The fuel materials require the user to specify a larger number of information than cladding materials, such as enrichment, grain size, theorethical density and density fraction. In addition, based on the specific models that are switched on, additional information might be required. As an example, the activation of the radial relocation model requires the user to specify the diametral gap and the fuel pellet diameter in cold conditions, together with the name of the fuel outer patch.
The rheology model (pure elasticity, visco-plasticity, etc.) to be used needs to be specified for each material too.
The following code snippet summarizes the choices for the case at hand:
materials
{
fuel
{
material UO2;
Tref Tref [ 0 0 0 1 0 ] 293; // Ref. T for thermal expansion
densificationModel UO2FRAPCON;
swellingModel UO2FRAPCON;
relocationModel UO2FRAPCON;
enrichment 0.045;
rGrain 5e-06; // Fuel grain radius
GdContent 0.0; // Gadolinium content
theoreticalDensity 10960; // Density UO2
densityFraction 0.95; // 1-fractional porosity
dishFraction 0.0;
resinteringDensityChange 1; // Allowed % shrinkage due to densification
GapCold 0.13e-3; // Cold diamteral gap size
DiamCold 9.0e-3; // Cold pellet diameter
recoveryFraction 0.5; // Fractional relocation strain allowed for recovery
outerPatch "fuelOuter";
isotropicCracking on; // Allow fuel softening to mimic cracking effect
nCracksMax 12; // Max. number of allowed cracks
rheologyModel elasticity; // Purely elastic behavior
}
cladding
{
material zircaloy;
Tref Tref [ 0 0 0 1 0 ] 293; // Ref. T for thermal expansion
PoissonRatioModel ZyConstant; // Constant poisson ratio
rheologyModel misesPlasticCreep; // Visco-plastic behavior
rheologyModelOptions
{
plasticStrainVsYieldStress table
(
(0 250e6)
);
creepModel LimbackCreepModel;
relax 1.0;
}
}
}
Time-Step Control
OFFBEAT provides a variety of options for the user to select criteria that eventually determine the adaptive time-step selection. The code selects the next time-step according to the criteria set in the system/controlDict dictionary:
// Enable adaptive time-step for the simulation
adjustableTimeStep true;
// Minimum and Maximum allowed time-step size
maxDeltaT 6.048e5;
minDeltaT 0.001;
// Relative allowed increase/decrease of time-step size between two adjacent time-steps
maxRelativeDeltaTIncrease 1e9;
minRelativeDeltaTDecrease 1e9;
// Relative allowed increase/decrease of avg. power between two adjacent time-steps
maxRelativePowerIncrease 1e9;
maxRelativePowerDecrease 1e9;
// Maximum allowed increase of Burnup in MWd/kg between two adjacent time-steps
maxBurnupIncrease 0.1;
// Maximum allowed increase of max. and avg. creep between two adjacent time-steps
maxMaximumCreep 1e-4;
maxAverageCreep 1e-4;
// Maximum allowed number of moles released during time-step
maxFGR 1e-8;
OFFBEAT also features a custom option to write to disk specific time-step folders defined by the user in a list form. If, for instance, the user desires to print only the final time step folder, adjustableWriteTimeStep can be enabled:
Warning
For adjustableWriteTimeStep to work as desired, the user should set the general printing option of OpenFOAM to writeControl timeStep; and writeInterval 1;.
Postprocessing
A simple probe postProcessor is configured for this simulation. This allow to 'measure' the value of one or more fields in a specific point of the simulation domain over time. In the case at hand, temperature and burnup evolution at the location (0 0 1.5) is recorded. The postprocessor is enabled in system/controlDict:
where probes is the name of a file located in system that reads:
points
(
(0.0 0.0 1.5) // Probe location
);
fields (T Bu); // Names of fields to be probed
#includeEtc "caseDicts/postProcessing/probes/probes.cfg"
// ************************************************************************* //
Upon run of the OFFBEAT code, this function produces output files in the postProcessing/probes/0/ folder with the time-evolution of the specified fields. For temperature for instance, the produced file reads:
Results
The plot in Figure 4 shows the results for the time evolution of temperature and burnup at the point (0 0 1.5) along irradiation, as extracted by the probes postProcessor.
To have an overall visualization of the fields computed by OFFBEAT, the code's results can be visualized through paraview. This is done by launching the command:
from the case folder. Through the ParaView GUI, the user can select and visualize the available fields mapped on the computational domain. Examples are provided in Figure 5, where temperature distribution in the fuel rod at the end of the irradiation is visualized, and in Figure 6, where the equivalent creep strain is mapped onto a rotational extrusion of the cladding geometry.
