- absolute_tolerance1e-11Absolute convergence tolerance for Newton iteration
Default:1e-11
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Absolute convergence tolerance for Newton iteration
- acceptable_multiplier10Factor applied to relative and absolute tolerance for acceptable convergence if iterations are no longer making progress
Default:10
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Factor applied to relative and absolute tolerance for acceptable convergence if iterations are no longer making progress
- adaptive_substeppingFalseUse adaptive substepping, where the number of substeps is successively doubled until the return mapping model successfully converges or the maximum number of substeps is reached.
Default:False
C++ Type:bool
Controllable:No
Description:Use adaptive substepping, where the number of substeps is successively doubled until the return mapping model successfully converges or the maximum number of substeps is reached.
- automatic_differentiation_return_mappingFalseWhether to use automatic differentiation to compute the derivative.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use automatic differentiation to compute the derivative.
- base_nameOptional parameter that defines a prefix for all material properties related to this stress update model. This allows for multiple models of the same type to be used without naming conflicts.
C++ Type:std::string
Controllable:No
Description:Optional parameter that defines a prefix for all material properties related to this stress update model. This allows for multiple models of the same type to be used without naming conflicts.
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
- declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
- fast_neutron_fluxThe fast neutron flux
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The fast neutron flux
- max_creep_increment0.001Maximum creep strain increment allowed by accuracy time step criterion
Default:0.001
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Maximum creep strain increment allowed by accuracy time step criterion
- max_inelastic_increment0.0001The maximum inelastic strain increment allowed in a time step
Default:0.0001
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The maximum inelastic strain increment allowed in a time step
- maximum_number_substeps25The maximum number of substeps allowed before cutting the time step.
Default:25
C++ Type:unsigned int
Controllable:No
Description:The maximum number of substeps allowed before cutting the time step.
- model_irradiation_creepTrueSet true to activate irradiation induced creep
Default:True
C++ Type:bool
Controllable:No
Description:Set true to activate irradiation induced creep
- model_thermal_creepTrueSet true to activate steady state thermal creep
Default:True
C++ Type:bool
Controllable:No
Description:Set true to activate steady state thermal creep
- outputThe reporting postprocessor to use for the max_iterations value.
C++ Type:PostprocessorName
Unit:(no unit assumed)
Controllable:No
Description:The reporting postprocessor to use for the max_iterations value.
- relative_tolerance1e-08Relative convergence tolerance for Newton iteration
Default:1e-08
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Relative convergence tolerance for Newton iteration
- temperatureThe coupled temperature (K)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The coupled temperature (K)
- use_substep_integration_errorFalseIf true, it establishes a substep size that will yield, at most,the creep numerical integration error given by substep_strain_tolerance.
Default:False
C++ Type:bool
Controllable:No
Description:If true, it establishes a substep size that will yield, at most,the creep numerical integration error given by substep_strain_tolerance.
- use_substeppingNONEWhether and how to use substepping
Default:NONE
C++ Type:MooseEnum
Controllable:No
Description:Whether and how to use substepping
- zircaloy_material_typeSTRESS_RELIEF_ANNEALEDType of zircaloy material properties to use in calculating creep. Note: ESCORE_IRRADIATIONGROWTHZR4 is not valid.
Default:STRESS_RELIEF_ANNEALED
C++ Type:MooseEnum
Controllable:No
Description:Type of zircaloy material properties to use in calculating creep. Note: ESCORE_IRRADIATIONGROWTHZR4 is not valid.
ZryCreepHayesHoppeUpdate
Computes the secondary thermal Hayes and Kassner secondary creep and the Hoppe irradiation creep for Zircaloy cladding. This material must be run in conjunction with ComputeMultipleInelasticStress.
Description
Hayes-Kassner secondary thermal creep and Hoppe secondary irradiation creep are both calculated in this class, ZryCreepHayesHoppeUpdate. This material, which must be run in conjunction with ComputeMultipleInelasticStress calculates the inelastic creep strain, the elastic strain, and the resulting stress for zircaloy materials. The contributions to creep from irradiation, primary, and thermal secondary creep are summed at each iteration.
warningwarning:Not Use with High Temperature
This material is not suitable for use in high temperatures, such as under LOCA conditions, and primary creep is not integrated with the secondary creep.
Thermal Creep in Standard Operating Conditions
The standard operating temperature thermal creep model used to calculate the secondary thermal creep with a power-law model is based on the work of Hayes and Kassner (2006). Hayes and Kassner found that zircaloy creep rate can be described with a conventional five-power creep law, as described in the review by Kassner and Pérez-Prado (2000).
Hayes-Kassner Secondary Thermal Creep
The secondary thermal creep strain rate equation implemented in BISON is (1) where is the effective thermal creep rate (1/s), is the effective (Mises) stress (Pa), is the activation energy (J/mol), is the universal gas constant (J/mol-K), is the temperature (K), is the shear modulus (Pa), and and are material constants.
The shear modulus value is taken directly from the elasticity tensor within BISON. Hayes and Kassner (2006) specify a creep law power (n) of 5. A value for is not reported in Hayes and Kassner (2006); however, based on experimental data presented in that paper, an approximate value of = 3.14 10 (1/s) was computed based on the relationship for the diffusion given in Moon et al. (2006); this value is used in BISON.
Primary Thermal Creep
Primary creep is not implemented for use with the Hayes-Kassner thermal creep model.
Irradiation Secondary Creep
Irradiation-induced creep of cladding materials is based on an empirical model developed by Hoppe (1991) that relates the creep rate to the current fast neutron flux and stress. The specific relation implemented is: where is the effective irradiation creep rate (1/s), is the fast neutron flux (n/m-s), is the effective (Mises) stress (MPa), and , , and are material constants. The material constant values = 3.557 10, = 0.85, and = 1.0 are for stress relief annealed zircaloy.
commentnote
The original Hoppe formulation is given in terms of circumferential stress, whereas the relation implemented in BISON assumes an effective (Mises) stress.
Example Input Syntax
[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
[zrycreep]
type = ZryCreepHayesHoppeUpdate<<<{"description": "Computes the secondary thermal Hayes and Kassner secondary creep and the Hoppe irradiation creep for Zircaloy cladding. This material must be run in conjunction with ComputeMultipleInelasticStress.", "href": "ZryCreepHayesHoppeUpdate.html"}>>>
temperature<<<{"description": "The coupled temperature (K)"}>>> = temp
fast_neutron_flux<<<{"description": "The fast neutron flux"}>>> = fast_neutron_flux
model_irradiation_creep<<<{"description": "Set true to activate irradiation induced creep"}>>> = true
model_thermal_creep<<<{"description": "Set true to activate steady state thermal creep"}>>> = true
[]
[]ZryCreepHayesHoppeUpdate must be run in conjunction with the inelastic strain return mapping stress calculator as shown below:
[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
[stress]
type = ComputeMultipleInelasticStress<<<{"description": "Compute state (stress and internal parameters such as plastic strains and internal parameters) using an iterative process. Combinations of creep models and plastic models may be used.", "href": "../ComputeMultipleInelasticStress.html"}>>>
tangent_operator<<<{"description": "Type of tangent operator to return. 'elastic': return the elasticity tensor. 'nonlinear': return the full, general consistent tangent operator."}>>> = elastic
inelastic_models<<<{"description": "The material objects to use to calculate stress and inelastic strains. Note: specify creep models first and plasticity models second."}>>> = 'zrycreep'
[]
[]Input Parameters
- apply_strainTrueFlag to apply strain. Used for testing.
Default:True
C++ Type:bool
Controllable:No
Description:Flag to apply strain. Used for testing.
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- effective_inelastic_strain_nameeffective_creep_strainName of the material property that stores the effective inelastic strain
Default:effective_creep_strain
C++ Type:std::string
Controllable:No
Description:Name of the material property that stores the effective inelastic strain
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- substep_strain_tolerance0.1Maximum ratio of the initial elastic strain increment at start of the return mapping solve to the maximum inelastic strain allowable in a single substep. Reduce this value to increase the number of substeps
Default:0.1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Maximum ratio of the initial elastic strain increment at start of the return mapping solve to the maximum inelastic strain allowable in a single substep. Reduce this value to increase the number of substeps
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- creeprate_scale_factor1scaling factor for total creep rate. Used for calibration and sensitivity studies
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:scaling factor for total creep rate. Used for calibration and sensitivity studies
Advanced: Scaling Factors Parameters
- internal_solve_full_iteration_historyFalseSet true to output full internal Newton iteration history at times determined by `internal_solve_output_on`. If false, only a summary is output.
Default:False
C++ Type:bool
Controllable:No
Description:Set true to output full internal Newton iteration history at times determined by `internal_solve_output_on`. If false, only a summary is output.
- internal_solve_output_onon_errorWhen to output internal Newton solve information
Default:on_error
C++ Type:MooseEnum
Options:never, on_error, always
Controllable:No
Description:When to output internal Newton solve information
Debug Parameters
- output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
- outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
Outputs Parameters
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Material Property Retrieval Parameters
Input Files
- (examples/2D_plane_strain_rod/planestrain.i)
- (test/tests/fuelrodlinevaluesampler/example_problem_smeared_test2.i)
- (examples/2D_plane_strain_fretting_wear/fretting-wear-initial.i)
- (test/tests/example_problem_test/example_problem_test.i)
- (test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep_rz.i)
- (examples/2D_plane_strain_fretting_wear/fretting-wear-initial-dyn-exc.i)
- (test/tests/side_ave_incr_tensor_component/side_ave_incr_component.i)
- (test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep.i)
- (test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep_highflux.i)
- (test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep_rz_highflux.i)
- (test/tests/fuelrodlinevaluesampler/example_problem_smeared_test.i)
Child Objects
References
- T. A. Hayes and M. Kassner.
Creep of zirconium and zirconium alloys.
Metallurgical and Materials Transactions A, 37A:2389–2396, 2006.[BibTeX]
@ARTICLE{hayes2006, author = "Hayes, T. A. and Kassner, M.", title = "Creep of Zirconium and Zirconium Alloys", journal = "Metallurgical and Materials Transactions A", year = "2006", volume = "37A", pages = "2389-2396" } - N. E. Hoppe.
Engineering model for zircaloy creep and growth.
In Proceedings of the ANS-ENS International Topical Meeting on LWR Fuel Performance, 157–172. Avignon, France, April 21-24, 1991.[BibTeX]
@INPROCEEDINGS{hoppe91, author = "Hoppe, N. E.", title = "Engineering model for zircaloy creep and growth", booktitle = "Proceedings of the ANS-ENS International Topical Meeting on {LWR} Fuel Performance", year = "1991", pages = "157-172", address = "Avignon, France", month = "April 21-24," } - ME Kassner and M-T Pérez-Prado.
Five-power-law creep in single phase metals and alloys.
Progress in Materials Science, 45(1):1–102, 2000.[BibTeX]
@article{kassner2000, author = "Kassner, ME and P{\'e}rez-Prado, M-T", title = "Five-power-law creep in single phase metals and alloys", journal = "Progress in Materials Science", volume = "45", number = "1", pages = "1-102", year = "2000", publisher = "Pergamon" } - J.H. Moon, P.E. Cantonwine, K.R. Anderson, S. Karthikeyan, and M.J. Mills.
Characterization and modeling of creep mechanisms in zircaloy-4.
Journal of Nuclear Materials, 353(3):177 – 189, 2006.
doi:10.1016/j.jnucmat.2006.01.023.[BibTeX]
@ARTICLE{Moon2006177, author = "Moon, J.H. and Cantonwine, P.E. and Anderson, K.R. and Karthikeyan, S. and Mills, M.J.", title = "Characterization and modeling of creep mechanisms in Zircaloy-4", journal = "Journal of Nuclear Materials", year = "2006", volume = "353", pages = "177 - 189", number = "3", bdsk-url-1 = "http://www.sciencedirect.com/science/article/pii/S0022311506001759", bdsk-url-2 = "http://dx.doi.org/10.1016/j.jnucmat.2006.01.023", doi = "10.1016/j.jnucmat.2006.01.023", issn = "0022-3115" }
(test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep.i)
#
# The mesh is a 1x1x1 cube with a constant pressure of 10 MPa on the top face.
# Symmetry boundary conditions on three planes provide a uniaxial stress
# field. The temperature is held constant at 1000; only creep is modeled: no
# platicity material model is included in the material blocks. The solution is
# advanced through ten time steps of 1e4 for a total time of 1e5.
#
# The total strain at time 1e5 can be computed as:
#
# e_tot = e_elas +
# e_thermal_creep +
# e_irradiation_creep
#
# = P/E +
# A * (sigma/G)**n * exp(-Q/(RT)) * dt +
# C0 * (phi)**C1 * (sigma/10**6)**C2 * dt
#
# where P = pressure load
# E = Young's modulus
# A = thermal creep coefficient
# G = shear modulus = 0.5*E/(1+nu)
# nu = Poisson's ratio
# sigma = stress
# n = thermal creep exponent
# Q = activation energy
# R = gas constant
# T = temperature
# C0 = irradiation creep coefficient
# phi = fast neutron flux
# C1 = neutron flux exponent
# C2 = stress exponent
# dt = total time
#
# For this test, the analytical solutuon at t = 1e5 is:
#
# G = 0.5*2e11/(1+0.3) = 7.692308e10
#
# e_tot = 10e6/2e11 +
# 3.14e24 * (10e6/G)**5 * exp(-2.7e5/(8.3143*1000) * 1e5
# 9.881e-28 * (1e19)**0.85 * (10e6/1e6)**1 * 1e5
#
# = 5e-5 +
# 9.189527e-5
# 1.395728e-5
#
# = 1.558526e-4
#
# For both small strain (ComputeIncrementalSmallStrain) or finite strain (ComputeFiniteStrain)
# kinematics, BISON computes:
#
# e_elas = 5e-5
# e_creep = 1.05819e-4
# e_tot = 1.55819e-4
#
#
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 3
[]
[]
[GlobalParams]
displacements = 'x_disp y_disp z_disp'
[]
[Variables]
[x_disp]
order = FIRST
family = LAGRANGE
[]
[y_disp]
order = FIRST
family = LAGRANGE
[]
[z_disp]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1000.0
[]
[]
[AuxVariables]
[fast_neutron_flux]
order = FIRST
family = LAGRANGE
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[elastic_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[creep_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[top_pull]
type = PiecewiseLinear
x = '0 1'
y = '1 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = finite
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
factor = 1e19
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[elastic_strain_yy]
type = RankTwoAux
rank_two_tensor = elastic_strain
variable = elastic_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[creep_strain_yy]
type = RankTwoAux
rank_two_tensor = creep_strain
variable = creep_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[]
[BCs]
[u_top_pull]
type = Pressure
variable = y_disp
boundary = top
factor = -10.0e6
function = top_pull
[]
[u_bottom_fix]
type = DirichletBC
variable = y_disp
boundary = bottom
value = 0.0
[]
[u_yz_fix]
type = DirichletBC
variable = x_disp
boundary = left
value = 0.0
[]
[u_xy_fix]
type = DirichletBC
variable = z_disp
boundary = back
value = 0.0
[]
[temp_fix]
type = DirichletBC
variable = temp
boundary = 'bottom top'
value = 1000.0
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 2.0e11
poissons_ratio = 0.3
[]
[stress]
type = ComputeMultipleInelasticStress
tangent_operator = elastic
inelastic_models = 'zrycreep'
[]
[zrycreep]
type = ZryCreepHayesHoppeUpdate
temperature = temp
fast_neutron_flux = fast_neutron_flux
model_irradiation_creep = true
model_thermal_creep = true
[]
[thermal]
type = HeatConductionMaterial
block = 0
specific_heat = 1.0
thermal_conductivity = 100.
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 1.0
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart'
petsc_options_value = '101'
line_search = 'none'
l_max_its = 100
nl_max_its = 100
l_tol = 1e-5
start_time = 0.0
num_steps = 10
dt = 10000
[]
[Outputs]
show = 'x_disp y_disp z_disp temp stress_yy elastic_strain_yy creep_strain_yy'
[out]
type = Exodus
[]
[]
(test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep.i)
#
# The mesh is a 1x1x1 cube with a constant pressure of 10 MPa on the top face.
# Symmetry boundary conditions on three planes provide a uniaxial stress
# field. The temperature is held constant at 1000; only creep is modeled: no
# platicity material model is included in the material blocks. The solution is
# advanced through ten time steps of 1e4 for a total time of 1e5.
#
# The total strain at time 1e5 can be computed as:
#
# e_tot = e_elas +
# e_thermal_creep +
# e_irradiation_creep
#
# = P/E +
# A * (sigma/G)**n * exp(-Q/(RT)) * dt +
# C0 * (phi)**C1 * (sigma/10**6)**C2 * dt
#
# where P = pressure load
# E = Young's modulus
# A = thermal creep coefficient
# G = shear modulus = 0.5*E/(1+nu)
# nu = Poisson's ratio
# sigma = stress
# n = thermal creep exponent
# Q = activation energy
# R = gas constant
# T = temperature
# C0 = irradiation creep coefficient
# phi = fast neutron flux
# C1 = neutron flux exponent
# C2 = stress exponent
# dt = total time
#
# For this test, the analytical solutuon at t = 1e5 is:
#
# G = 0.5*2e11/(1+0.3) = 7.692308e10
#
# e_tot = 10e6/2e11 +
# 3.14e24 * (10e6/G)**5 * exp(-2.7e5/(8.3143*1000) * 1e5
# 9.881e-28 * (1e19)**0.85 * (10e6/1e6)**1 * 1e5
#
# = 5e-5 +
# 9.189527e-5
# 1.395728e-5
#
# = 1.558526e-4
#
# For both small strain (ComputeIncrementalSmallStrain) or finite strain (ComputeFiniteStrain)
# kinematics, BISON computes:
#
# e_elas = 5e-5
# e_creep = 1.05819e-4
# e_tot = 1.55819e-4
#
#
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 3
[]
[]
[GlobalParams]
displacements = 'x_disp y_disp z_disp'
[]
[Variables]
[x_disp]
order = FIRST
family = LAGRANGE
[]
[y_disp]
order = FIRST
family = LAGRANGE
[]
[z_disp]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1000.0
[]
[]
[AuxVariables]
[fast_neutron_flux]
order = FIRST
family = LAGRANGE
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[elastic_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[creep_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[top_pull]
type = PiecewiseLinear
x = '0 1'
y = '1 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = finite
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
factor = 1e19
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[elastic_strain_yy]
type = RankTwoAux
rank_two_tensor = elastic_strain
variable = elastic_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[creep_strain_yy]
type = RankTwoAux
rank_two_tensor = creep_strain
variable = creep_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[]
[BCs]
[u_top_pull]
type = Pressure
variable = y_disp
boundary = top
factor = -10.0e6
function = top_pull
[]
[u_bottom_fix]
type = DirichletBC
variable = y_disp
boundary = bottom
value = 0.0
[]
[u_yz_fix]
type = DirichletBC
variable = x_disp
boundary = left
value = 0.0
[]
[u_xy_fix]
type = DirichletBC
variable = z_disp
boundary = back
value = 0.0
[]
[temp_fix]
type = DirichletBC
variable = temp
boundary = 'bottom top'
value = 1000.0
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 2.0e11
poissons_ratio = 0.3
[]
[stress]
type = ComputeMultipleInelasticStress
tangent_operator = elastic
inelastic_models = 'zrycreep'
[]
[zrycreep]
type = ZryCreepHayesHoppeUpdate
temperature = temp
fast_neutron_flux = fast_neutron_flux
model_irradiation_creep = true
model_thermal_creep = true
[]
[thermal]
type = HeatConductionMaterial
block = 0
specific_heat = 1.0
thermal_conductivity = 100.
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 1.0
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart'
petsc_options_value = '101'
line_search = 'none'
l_max_its = 100
nl_max_its = 100
l_tol = 1e-5
start_time = 0.0
num_steps = 10
dt = 10000
[]
[Outputs]
show = 'x_disp y_disp z_disp temp stress_yy elastic_strain_yy creep_strain_yy'
[out]
type = Exodus
[]
[]
(examples/2D_plane_strain_rod/planestrain.i)
initial_fuel_density = 10431.0
[GlobalParams]
temperature = temp
displacements = 'disp_x disp_y'
order = SECOND
family = LAGRANGE
energy_per_fission = 3.2e-11 # J/fission
[]
[Mesh]
patch_size = 100 # For contact algorithm
[mesh]
type = FileMeshGenerator
file = planestrain.e
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[temp]
initial_condition = 580.0 # set initial temp to ambient
[]
[]
[AuxVariables]
[fission_rate]
block = pellet_type_1
[]
[burnup]
block = pellet_type_1
[]
[fast_neutron_flux]
block = clad
[]
[fast_neutron_fluence]
block = clad
[]
[relocation_strain]
order = CONSTANT
family = MONOMIAL
[]
[gap_cond]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear # reads and interpolates an input file containing rod average linear power vs time
data_file = powerhistory.csv
scale_factor = 1
[]
[axial_peaking_factors]
type = ConstantFunction
value = 1
[]
[pressure_ramp] # reads and interpolates input data defining amplitude curve for fill gas pressure
type = PiecewiseLinear
x = '0 1e4'
y = '0 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[pellets]
block = pellet_type_1
strain = FINITE
planar_formulation = PLANE_STRAIN
eigenstrain_names = 'fuel_relocation_eigenstrain fuel_thermal_eigenstrain
fuel_volumetric_eigenstrain'
generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
decomposition_method = EigenSolution
[]
[clad]
block = clad
strain = FINITE
planar_formulation = PLANE_STRAIN
eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_eigenstrain'
generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
decomposition_method = EigenSolution
[]
[]
[Kernels]
[heat] # gradient term in heat conduction equation
type = HeatConduction
variable = temp
[]
[heat_ie] # time term in heat conduction equation
type = HeatConductionTimeDerivative
variable = temp
[]
[heat_source] # source term in heat conduction equation
type = NeutronHeatSource
variable = temp
block = pellet_type_1
fission_rate = fission_rate
[]
[]
[Burnup]
[burnup]
block = pellet_type_1
rod_ave_lin_pow = power_history # using the power function defined above
axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
num_radial = 80
num_axial = 21
axial_direction = z
density = ${initial_fuel_density}
a_lower = -1e-3 # mesh dependent!
a_upper = 1e-3 # mesh dependent!
fuel_inner_radius = 0
fuel_outer_radius = .0041
fuel_volume_ratio = 0.987775 # for use with dished pellets (ratio of actual volume to cylinder volume)
#N235 = N235 # Activate to write N235 concentration to output file
#N238 = N238 # Activate to write N238 concentration to output file
#N239 = N239 # Activate to write N239 concentration to output file
#N240 = N240 # Activate to write N240 concentration to output file
#N241 = N241 # Activate to write N241 concentration to output file
#N242 = N242 # Activate to write N242 concentration to output file
RPF = RPF
[]
[]
[AuxKernels]
# Define auxilliary kernels for each of the aux variables
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
block = clad
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
factor = 3e13
execute_on = timestep_begin
[]
[fast_neutron_fluence]
type = FastNeutronFluenceAux
variable = fast_neutron_fluence
block = clad
fast_neutron_flux = fast_neutron_flux
execute_on = timestep_begin
[]
[relocation_strain]
type = MaterialRealAux
property = relocation_strain
variable = relocation_strain
block = pellet_type_1
execute_on = timestep_end
[]
[conductance]
type = MaterialRealAux
property = gap_conductance
variable = gap_cond
boundary = 8
execute_on = linear
[]
[]
[Contact]
# Define mechanical contact between the fuel (sideset=10) and the clad (sideset=5)
[pellet_clad_mechanical]
primary = 5
secondary = 10
penalty = 1e7
[]
[]
[ThermalContact]
# Define thermal contact between the fuel (sideset=10) and the clad (sideset=5)
[thermal_contact]
type = GasGapHeatTransfer
variable = temp
primary = 7
secondary = 8
initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
gas_released = fission_gas_released # coupling to a postprocessor which supplies the fission gas addition
[]
[]
[BCs]
# Define boundary conditions
[no_y_all] # pin pellets and clad along axis of symmetry (y)
type = DirichletBC
variable = disp_y
boundary = 15
value = 0.0
[]
[no_x_all] # pin pellets and clad along axis of symmetry (x)
type = DirichletBC
variable = disp_x
boundary = 16
value = 0.0
[]
[Pressure] # apply coolant pressure on clad outer walls
[coolantPressure]
boundary = '2'
factor = 15.5e6
function = pressure_ramp # use the pressure_ramp function defined above
[]
[]
[PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
[plenumPressure]
boundary = 9
initial_pressure = 2.0e6
R = 8.3143
output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
temperature = plenum_temperature # coupling to post processor to get gas temperature approximation
volume = plenum_volume # coupling to post processor to get gas volume
material_input = fission_gas_released # coupling to post processor to get fission gas added
output = plenum_pressure # coupling to post processor to output plenum/gap pressure
displacements = 'disp_x disp_y'
[]
[]
[convective_clad_surface] # apply convective boundary to clad outer surface
type = ConvectiveFluxBC
boundary = '2'
variable = temp
rate = 38200.0 #convection coefficient (h)
initial = 580.0
final = 580.0
duration = 1.0e4 #duration of initial power ramp
[]
[]
[Materials]
# Define material behavior models and input material property data
[fuel_thermal] # temperature and burnup dependent thermal properties of UO2 (BISON kernel)
type = UO2Thermal
thermal_conductivity_model = FINK_LUCUTA
block = pellet_type_1
temperature = temp
burnup = burnup
initial_porosity = 0.0
[]
[fuel_solid_mechanics_swelling] # free expansion strains (swelling and densification) for UO2 (BISON kernel)
type = UO2VolumetricSwellingEigenstrain
gas_swelling_model_type = MATPRO
block = pellet_type_1
burnup = burnup
initial_fuel_density = 10431.0
temperature = temp
eigenstrain_name = 'fuel_volumetric_eigenstrain'
[]
[fuel_creep]
type = UO2CreepUpdate
block = pellet_type_1
temperature = temp
fission_rate = fission_rate
density = 10431.0
initial_grain_radius = 10.0e-6
oxygen_to_metal_ratio = 2.0
[]
[fuel_elasticity_tensor]
type = ComputeIsotropicElasticityTensor
block = pellet_type_1
youngs_modulus = 2.0e11
poissons_ratio = 0.345
[]
[fuel_stress]
type = ComputeMultipleInelasticStress
block = pellet_type_1
inelastic_models = 'fuel_creep'
[]
[fuel_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = pellet_type_1
thermal_expansion_coeff = 10.0e-6
temperature = temp
stress_free_temperature = 580.0
eigenstrain_name = 'fuel_thermal_eigenstrain'
[]
[fuel_relocation]
type = UO2RelocationEigenstrain
block = pellet_type_1
burnup = burnup
diameter = 0.0082
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
diametral_gap =160e-6
burnup_relocation_stop = 1.e20
relocation_activation1 = 5000
axial_direction = z
eigenstrain_name = 'fuel_relocation_eigenstrain'
[]
[clad_thermal]
type = HeatConductionMaterial
block = clad
thermal_conductivity = 16.0
specific_heat = 330.0
[]
[clad_elasticity_tensor]
type = ZryElasticityTensor
block = clad
[]
[clad_creep_model]
type = ZryCreepHayesHoppeUpdate
block = clad
fast_neutron_flux = fast_neutron_flux
temperature = temp
zircaloy_material_type = stress_relief_annealed
model_irradiation_creep = true
model_thermal_creep = true
[]
[clad_stress]
type = ComputeMultipleInelasticStress
block = clad
tangent_operator = elastic
inelastic_models = 'clad_creep_model'
[]
[clad_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = clad
thermal_expansion_coeff = 5.0e-6
temperature = temp
stress_free_temperature = 580.0
eigenstrain_name = 'clad_thermal_eigenstrain'
[]
[clad_irrgrowth]
type = ZryIrradiationGrowthEigenstrain
block = clad
fast_neutron_fluence = fast_neutron_fluence
axial_direction = 2
zircaloy_material_type = ESCORE_IrradiationGrowthZr4
eigenstrain_name = 'clad_irradiation_eigenstrain'
[]
[fission_gas_release] # Forsberg-Massih fission gas release mode
type = UO2Sifgrs
block = pellet_type_1
temperature = temp
fission_rate = fission_rate # coupling to fission_rate aux variable
grain_radius = 10.0e-6
#external_pressure = 40e6
[]
[clad_density]
type = StrainAdjustedDensity
block = clad
strain_free_density = 6551.0
[]
[fuel_density]
type = StrainAdjustedDensity
block = pellet_type_1
strain_free_density = 10431.0
[]
[]
[Dampers]
[limitT]
type = MaxIncrement
max_increment = 100.0
variable = temp
[]
[limitX]
type = MaxIncrement
max_increment = 1e-5
variable = disp_x
[]
[]
[Executioner]
# PETSC options:
# petsc_options
# petsc_options_iname
# petsc_options_value
#
# controls for linear iterations
# l_max_its
# l_tol
#
# controls for nonlinear iterations
# nl_max_its
# nl_rel_tol
# nl_abs_tol
#
# time control
# start_time
# dt
# optimal_iterations
# iteration_window
# linear_iteration_ratio
type = Transient
solve_type = 'PJFNK'
petsc_options = '-ksp_gmres_modifiedgramschmidt'
petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_composite_pcs -sub_0_pc_hypre_type -sub_0_pc_hypre_boomeramg_max_iter -sub_0_pc_hypre_boomeramg_grid_sweeps_all -sub_1_sub_pc_type -pc_composite_type -ksp_type -mat_mffd_type'
petsc_options_value = '201 composite hypre,asm boomeramg 2 2 lu multiplicative fgmres ds'
line_search = 'none'
l_max_its = 100
l_tol = 8e-3
nl_max_its = 15
nl_rel_tol = 1e-4
nl_abs_tol = 1e-10
start_time = 0.0
end_time = 1.0e6
[TimeStepper]
type = IterationAdaptiveDT
dt = 2.0e2
time_t = '1e4 1e5 1e6'
time_dt = '1e3 1e4 1e5'
[]
dtmax = 2e6
dtmin = 1
# optimal_iterations = 6
# iteration_window = 2
# linear_iteration_ratio = 100
[Quadrature]
order = THIRD
[]
[]
[Postprocessors]
# Define postprocessors (some are required as specified above; others are optional; many others are available)
[average_interior_clad_temperature] # average temperature of cladding interior
type = SideAverageValue
boundary = 7
variable = temp
execute_on = 'initial timestep_end'
[]
[average_centerline_fuel_temperature] # average temperature of the cladding interior and all pellet exteriors
type = SideAverageValue
boundary = 9
variable = temp
execute_on = 'initial linear'
[]
[plenum_temperature]
type = SideAverageValue
boundary = 9
variable = temp
execute_on = 'initial timestep_end'
[]
[plenum_volume] # gas volume
type = InternalVolume
boundary = 9
addition = 1.3e-5 #rough guess of plenum volume/unit length of fuel
execute_on = 'initial linear'
[]
[pellet_volume] # fuel pellet total volume
type = InternalVolume
boundary = 8
execute_on = 'initial timestep_end'
[]
[clad_inner_vol] # volume inside of cladding
type = InternalVolume
boundary = 7
outputs = exodus
execute_on = 'initial timestep_end'
[]
[fission_gas_generated] # fission gas produced (moles)
type = ElementIntegralFisGasGeneratedSifgrs
block = pellet_type_1
execute_on = linear
[]
[fission_gas_released] # fission gas released to plenum (moles)
type = ElementIntegralFisGasReleasedSifgrs
block = pellet_type_1
execute_on = linear
[]
[flux_from_clad] # area integrated heat flux from the cladding
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 5
diffusivity = thermal_conductivity
execute_on = timestep_end
[]
[flux_from_fuel] # area integrated heat flux from the fuel
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 10
diffusivity = thermal_conductivity
execute_on = timestep_end
[]
[_dt] # time step
type = TimestepSize
execute_on = timestep_end
[]
[num_lin_it]
type = NumLinearIterations
[]
[num_nonlin_it]
type = NumNonlinearIterations
[]
[tot_lin_it]
type = CumulativeValuePostprocessor
postprocessor = num_lin_it
[]
[tot_nonlin_it]
type = CumulativeValuePostprocessor
postprocessor = num_nonlin_it
[]
[alive_time]
type = PerfGraphData
section_name = Root
data_type = TOTAL
[]
[rod_total_power]
type = ElementIntegralPower
variable = temp
fission_rate = fission_rate
block = pellet_type_1
execute_on = timestep_end
[]
[rod_input_power]
type = FunctionValuePostprocessor
function = power_history
scale_factor = 0.1186 # rod height
execute_on = timestep_end
[]
[fission_gas_released_percentage]
type = FGRPercent
fission_gas_released = fission_gas_released
fission_gas_generated = fission_gas_generated
[]
[]
[Outputs]
perf_graph = true
exodus = true
csv = true
[console]
type = Console
max_rows = 25
[]
[]
(test/tests/fuelrodlinevaluesampler/example_problem_smeared_test2.i)
[GlobalParams]
density = 10431.0
displacements = 'disp_x disp_y'
energy_per_fission = 3.2e-11 # J/fission
temperature = temp
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
[]
[Mesh]
coord_type = RZ
patch_update_strategy = auto
patch_size = 10
partitioner = centroid
centroid_partitioner_direction = y
[mesh]
type = FileMeshGenerator
file = SmearedTwoPelletOneType2D.e
[]
[]
[Variables]
[temp]
initial_condition = 580.0
[]
[]
[AuxVariables]
[fast_neutron_flux]
block = clad
[]
[fast_neutron_fluence]
block = clad
[]
[grain_radius]
block = pellet_type_1
initial_condition = 10e-6
[]
[gap_cond]
order = CONSTANT
family = MONOMIAL
[]
[coolant_htc]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear
data_file = powerhistory.csv
scale_factor = 1
[]
[axial_peaking_factors]
type = ParsedFunction
expression = 1
[]
[pressure_ramp]
type = PiecewiseLinear
x = '-200 0'
y = '0 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[fuel]
block = pellet_type_1
strain = FINITE
incremental = true
extra_vector_tags = 'ref'
add_variables = true
decomposition_method = EigenSolution
eigenstrain_names = 'fuel_volumetric_swelling_eigenstrain
fuel_relocation_eigenstrain fuel_thermal_eigenstrain'
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress'
[]
[clad]
block = clad
strain = FINITE
incremental = true
extra_vector_tags = 'ref'
add_variables = true
decomposition_method = EigenSolution
eigenstrain_names = 'clad_thermal_strain clad_irradiation_growth_eigenstrain'
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress'
[]
[]
[Kernels]
[gravity]
type = Gravity
variable = disp_y
value = -9.81
[]
[heat]
type = HeatConduction
variable = temp
extra_vector_tags = 'ref'
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
extra_vector_tags = 'ref'
[]
[heat_source]
type = NeutronHeatSource
variable = temp
extra_vector_tags = 'ref'
block = pellet_type_1
burnup_function = burnup
[]
[]
[UserObjects]
[pin_geometry]
type = FuelPinGeometry
[]
[]
[Burnup]
[burnup]
block = pellet_type_1
rod_ave_lin_pow = power_history # using the power function defined above
axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
num_radial = 80
num_axial = 11
fuel_pin_geometry = 'pin_geometry'
fuel_volume_ratio = 0.987775 # for use with dished pellets (ratio of actual volume to cylinder volume)
order = CONSTANT
family = MONOMIAL
RPF = RPF
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
block = clad
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
factor = 3e13
execute_on = timestep_begin
[]
[fast_neutron_fluence]
type = FastNeutronFluenceAux
variable = fast_neutron_fluence
block = clad
fast_neutron_flux = fast_neutron_flux
execute_on = timestep_begin
[]
[grain_radius]
type = GrainRadiusAux
block = pellet_type_1
variable = grain_radius
temperature = temp
execute_on = linear
[]
[conductance]
type = MaterialRealAux
property = gap_conductance
variable = gap_cond
boundary = 10
execute_on = 'initial timestep_end'
[]
[coolant_htc]
type = MaterialRealAux
property = coolant_channel_htc
variable = coolant_htc
boundary = 2
execute_on = 'initial timestep_end'
[]
[]
[Contact]
[pellet_clad_mechanical]
primary = 5
secondary = 10
formulation = KINEMATIC
model = frictionless
normalize_penalty = true
penalty = 1e14
[]
[]
[ThermalContact]
[thermal_contact]
type = GasGapHeatTransfer
variable = temp
primary = 5
secondary = 10
initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
gas_released = fis_gas_released # coupling to a postprocessor which supplies the fission gas addition
contact_pressure = contact_pressure
[]
[]
[BCs]
[no_x_all] # pin pellets and clad along axis of symmetry (y)
type = DirichletBC
variable = disp_x
boundary = 12
value = 0.0
[]
[no_y_clad_bottom] # pin clad bottom in the axial direction (y)
type = DirichletBC
variable = disp_y
boundary = '1'
value = 0.0
[]
[no_y_fuel_bottom] # pin fuel bottom in the axial direction (y)
type = DirichletBC
variable = disp_y
boundary = '1020'
value = 0.0
[]
[Pressure] # apply coolant pressure on clad outer walls
[coolantPressure]
boundary = '1 2 3'
factor = 15.5e6
function = pressure_ramp # use the pressure_ramp function defined above
[]
[]
[PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
[plenumPressure]
boundary = 9
initial_pressure = 2.0e6
startup_time = -200
R = 8.3143
output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
temperature = ave_temp_interior # coupling to post processor to get gas temperature approximation
volume = gas_volume # coupling to post processor to get gas volume
material_input = fis_gas_released # coupling to post processor to get fission gas added
output = plenum_pressure # coupling to post processor to output plenum/gap pressure
displacements = 'disp_x disp_y'
execute_on = 'initial linear'
[]
[]
[]
[CoolantChannel]
[convective_clad_surface] # apply convective boundary to clad outer surface
boundary = '1 2 3'
variable = temp
inlet_temperature = 580 # K
inlet_pressure = 15.5e6 # Pa
inlet_massflux = 3800 # kg/m^2-sec
rod_diameter = 0.948e-2 # m
rod_pitch = 1.26e-2 # m
linear_heat_rate = power_history
axial_power_profile = axial_peaking_factors
[]
[]
[Materials]
[fuel_thermal]
type = UO2Thermal
block = pellet_type_1
thermal_conductivity_model = NFIR
initial_porosity = 0.0
temperature = temp
burnup_function = burnup
[]
[fuel_volumetric_swelling]
type = UO2VolumetricSwellingEigenstrain
block = pellet_type_1
burnup_function = burnup
initial_fuel_density = 10431.0
eigenstrain_name = 'fuel_volumetric_swelling_eigenstrain'
[]
[fuel_elasticity_tensor]
type = UO2ElasticityTensor
block = pellet_type_1
[]
[fuel_thermal_expansion]
type = UO2ThermalExpansionMartinEigenstrain
block = pellet_type_1
stress_free_temperature = 295
eigenstrain_name = 'fuel_thermal_eigenstrain'
[]
[hotpressing]
type = UO2HotPressingCreepUpdate
block = pellet_type_1
burnup_function = burnup
initial_grain_radius = 10.0e-6
[]
[radial_return_stress]
type = ComputeMultipleInelasticStress
tangent_operator = elastic
inelastic_models = 'hotpressing'
block = pellet_type_1
[]
[fuel_relocation]
type = UO2RelocationEigenstrain
block = pellet_type_1
burnup_function = burnup
fuel_pin_geometry = 'pin_geometry'
relocation_activation1 = 5000 #TM default value
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
burnup_relocation_stop = 1.e20
eigenstrain_name = fuel_relocation_eigenstrain
[]
[clad_thermal]
type = HeatConductionMaterial
block = clad
thermal_conductivity = 16.0
specific_heat = 330.0
[]
[clad_elasticity_tensor]
type = ZryElasticityTensor
block = clad
[]
[clad_creep_model]
type = ZryCreepHayesHoppeUpdate
block = clad
fast_neutron_flux = fast_neutron_flux
model_irradiation_creep = true
model_thermal_creep = true
[]
[clad_inelastic_stress]
type = ComputeMultipleInelasticStress
block = clad
tangent_operator = elastic
inelastic_models = 'clad_creep_model'
[]
[clad_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = clad
thermal_expansion_coeff = 5.0e-6
stress_free_temperature = 295.0
eigenstrain_name = clad_thermal_strain
[]
[clad_irradiation_growth]
type = ZryIrradiationGrowthEigenstrain
block = clad
fast_neutron_fluence = fast_neutron_fluence
zircaloy_material_type = ESCORE_IrradiationGrowthZr4
eigenstrain_name = clad_irradiation_growth_eigenstrain
[]
[fission_gas_release]
type = UO2Sifgrs
block = pellet_type_1
temperature = temp
burnup_function = burnup
grain_radius = grain_radius
gbs_model = true
[]
[clad_density]
type = StrainAdjustedDensity
block = clad
strain_free_density = 6551.0
[]
[fuel_density]
type = StrainAdjustedDensity
block = pellet_type_1
strain_free_density = 10431.0
[]
[]
[Dampers]
[BoundingValueNodalDamper]
type = BoundingValueNodalDamper
variable = temp
max_value = 3200
min_value = 300
[]
[]
[Preconditioning]
[SMP]
type = SMP
coupled_groups = 'disp_x,disp_y'
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-pc_type_asm'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
petsc_options_value = 'lu superlu_dist 51'
line_search = 'none'
verbose = true
l_max_its = 100
l_tol = 1e-5 #8e-3
nl_max_its = 15
nl_rel_tol = 1e-10
nl_abs_tol = 1e-8
start_time = -200
num_steps = 2
dtmax = 2e6
dtmin = 1
[TimeStepper]
type = IterationAdaptiveDT
dt = 2.0e2
optimal_iterations = 6
iteration_window = 2
[]
[]
[Postprocessors]
[ave_temp_interior] # average temperature of the cladding interior and all pellet exteriors
type = SideAverageValue
boundary = 9
variable = temp
execute_on = 'initial linear'
[]
[clad_inner_vol] # volume inside of cladding
type = InternalVolume
boundary = 7
outputs = exodus
execute_on = 'initial timestep_end'
[]
[pellet_volume] # fuel pellet total volume
type = InternalVolume
boundary = 8
outputs = exodus
execute_on = 'initial timestep_end'
[]
[avg_clad_temp] # average temperature of cladding interior
type = SideAverageValue
boundary = 7
variable = temp
execute_on = 'initial timestep_end'
[]
[fis_gas_produced] # fission gas produced (moles)
type = ElementIntegralFisGasGeneratedSifgrs
block = pellet_type_1
execute_on = timestep_end
[]
[fis_gas_released] # fission gas released to plenum (moles)
type = ElementIntegralFisGasReleasedSifgrs
block = pellet_type_1
execute_on = timestep_end
[]
[fis_gas_grain]
type = ElementIntegralFisGasGrainSifgrs
block = pellet_type_1
outputs = exodus
execute_on = timestep_end
[]
[fis_gas_boundary]
type = ElementIntegralFisGasBoundarySifgrs
block = pellet_type_1
outputs = exodus
execute_on = timestep_end
[]
[gas_volume] # gas volume
type = InternalVolume
boundary = 9
component = 1
execute_on = 'initial linear'
[]
[flux_from_clad] # area integrated heat flux from the cladding
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 5
diffusivity = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[flux_from_fuel] # area integrated heat flux from the fuel
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 10
diffusivity = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[_dt] # time step
type = TimestepSize
execute_on = timestep_end
[]
[nonlinear_its] # number of nonlinear iterations at each timestep
type = NumNonlinearIterations
execute_on = timestep_end
[]
[rod_total_power]
type = ElementIntegralPower
variable = temp
burnup_function = burnup
block = pellet_type_1
execute_on = 'initial timestep_end'
[]
[rod_input_power]
type = FunctionValuePostprocessor
function = power_history
scale_factor = 0.02372 # rod height
execute_on = 'initial timestep_end'
[]
[]
[VectorPostprocessors]
[fuel_vonmises]
type = FuelRodLineValueSampler
variable = vonmises_stress
material = 'fuel'
fraction = 0.51
num_points = 20
orientation = 'vertical'
fuel_pin_geometry = 'pin_geometry'
outputs = chkfile
[]
[clad_vonmises]
type = FuelRodLineValueSampler
variable = vonmises_stress
material = 'clad'
fraction = 0.51
num_points = 20
orientation = 'vertical'
fuel_pin_geometry = 'pin_geometry'
outputs = chkfile
[]
[]
[Outputs]
exodus = true
color = false
csv = true
[console]
type = Console
output_linear = true
max_rows = 25
[]
[chkfile]
type = CSV
execute_on = 'FINAL'
[]
[]
(examples/2D_plane_strain_fretting_wear/fretting-wear-initial.i)
initial_fuel_density = 10431.0
[GlobalParams]
temperature = temp
displacements = 'disp_x disp_y'
order = FIRST
family = LAGRANGE
energy_per_fission = 3.2e-11 # J/fission
volumetric_locking_correction = true
[]
[Mesh]
[file]
type = FileMeshGenerator
file = refined_excitation_better_mesh.e
[]
construct_node_list_from_side_list = true
patch_size = 100 # For contact algorithm
[]
[Variables]
[temp]
initial_condition = 580.0 # set initial temp to ambient
[]
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
group_variables = 'disp_x disp_y '
converge_on = 'disp_x disp_y temp'
material_coverage_check = false
kernel_coverage_check = false
# restart_file_base = planestrain_grid_aux_vars_out_cp/LATEST
[]
[AuxVariables]
[fission_rate]
block = pellet_type_1
[]
[burnup]
block = pellet_type_1
[]
[fast_neutron_flux]
block = 'clad grid'
[]
[fast_neutron_fluence]
block = 'clad grid'
[]
[relocation_strain]
order = CONSTANT
family = MONOMIAL
[]
[worn_depth]
order = FIRST
family = LAGRANGE
block = 'spacer_clad_mechanical_secondary_subdomain'
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear # reads and interpolates an input file containing rod average linear power vs time
data_file = powerhistory.csv
scale_factor = 1
[]
[axial_peaking_factors]
type = ConstantFunction
value = 1
[]
[pressure_var] # reads and interpolates input data defining amplitude curve for fill gas pressure
type = PiecewiseLinear
x = '0 1e4'
y = '0 1'
[]
[pressure_var_variable] # reads and interpolates input data defining amplitude curve for fill gas pressure
type = ParsedFunction
expression = 'if(t < 1e4, 1, 1 + sin((t-1e4)*pi/10.0) * (t-1e4))'
[]
[]
[Physics/SolidMechanics/Dynamic]
[pellets]
add_variables = true
newmark_beta = 0.25
newmark_gamma = 0.5
block = pellet_type_1
strain = FINITE
planar_formulation = PLANE_STRAIN
eigenstrain_names = 'fuel_relocation_eigenstrain fuel_thermal_eigenstrain
fuel_volumetric_eigenstrain'
generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
decomposition_method = EigenSolution
temperature = temp
extra_vector_tags = 'ref'
[]
[clad]
add_variables = true
newmark_beta = 0.25
newmark_gamma = 0.5
block = clad
strain = FINITE
planar_formulation = PLANE_STRAIN
eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_eigenstrain'
generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
decomposition_method = EigenSolution
temperature = temp
extra_vector_tags = 'ref'
[]
[grid]
add_variables = true
newmark_beta = 0.25
newmark_gamma = 0.5
block = grid
strain = FINITE
planar_formulation = PLANE_STRAIN
eigenstrain_names = 'grid_thermal_eigenstrain grid_irradiation_eigenstrain'
generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
decomposition_method = EigenSolution
temperature = temp
extra_vector_tags = 'ref'
[]
[]
[Kernels]
[heat] # gradient term in heat conduction equation
type = HeatConduction
variable = temp
block = 'pellet_type_1 clad grid'
extra_vector_tags = 'ref'
[]
[heat_ie] # time term in heat conduction equation
type = HeatConductionTimeDerivative
variable = temp
block = 'pellet_type_1 clad'
extra_vector_tags = 'ref'
[]
[heat_source] # source term in heat conduction equation
type = NeutronHeatSource
variable = temp
block = pellet_type_1
fission_rate = fission_rate
extra_vector_tags = 'ref'
[]
[]
[Contact]
# Define mechanical contact between the fuel (sideset=10) and the clad (sideset=5)
[spacer_clad_mechanical]
formulation = mortar
model = coulomb
primary = 101
secondary = 102
c_normal = 1e+16 # 1e+7
c_tangential = 1e+20
friction_coefficient = 0.4
# Do not apply dynamic stabilization
newmark_beta = 0.0001
newmark_gamma = 0.5
capture_tolerance = 0.0
mortar_dynamics = true
interpolate_normals = false
generate_mortar_mesh = true
wear_depth = worn_depth
[]
[]
[Contact]
# Define mechanical contact between the fuel (sideset=10) and the clad (sideset=5)
[pellet_clad_mechanical_real]
formulation = mortar
model = frictionless
primary = 7
secondary = 8
c_normal = 1e+16 #
c_tangential = 1e+16
friction_coefficient = 0.4
# Do not apply dynamic stabilization
newmark_beta = 0.0001
newmark_gamma = 0.5
capture_tolerance = 0.0
mortar_dynamics = true
interpolate_normals = false
generate_mortar_mesh = true
[]
[]
[ThermalContactMortar]
[thermal_contact]
secondary_variable = temp
primary_boundary = 7
secondary_boundary = 8
initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
gas_released = fission_gas_released # coupling to a postprocessor which supplies the fission gas addition
[]
[]
[Burnup]
[burnup]
block = pellet_type_1
rod_ave_lin_pow = power_history # using the power function defined above
axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
num_radial = 80
num_axial = 21
axial_axis = 2
density = ${initial_fuel_density}
a_lower = -1e-3 # mesh dependent!
a_upper = 1e-3 # mesh dependent!
fuel_inner_radius = 0
fuel_outer_radius = .0041
fuel_volume_ratio = 0.987775 # for use with dished pellets (ratio of actual volume to cylinder volume)
#N235 = N235 # Activate to write N235 concentration to output file
#N238 = N238 # Activate to write N238 concentration to output file
#N239 = N239 # Activate to write N239 concentration to output file
#N240 = N240 # Activate to write N240 concentration to output file
#N241 = N241 # Activate to write N241 concentration to output file
#N242 = N242 # Activate to write N242 concentration to output file
RPF = RPF
[]
[]
[AuxKernels]
[worn_depth]
type = MortarArchardsLawAux
variable = worn_depth
primary_boundary = 101
secondary_boundary = 102
primary_subdomain = 'spacer_clad_mechanical_primary_subdomain'
secondary_subdomain = 'spacer_clad_mechanical_secondary_subdomain'
displacements = 'disp_x disp_y'
friction_coefficient = 0.5
energy_wear_coefficient = 0.1e-9
normal_pressure = spacer_clad_mechanical_normal_lm
execute_on = 'TIMESTEP_END'
[]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
block = clad
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
factor = 3e13
execute_on = timestep_begin
[]
[fast_neutron_fluence]
type = FastNeutronFluenceAux
variable = fast_neutron_fluence
block = clad
fast_neutron_flux = fast_neutron_flux
execute_on = timestep_begin
[]
[relocation_strain]
type = MaterialRealAux
property = relocation_strain
variable = relocation_strain
block = pellet_type_1
execute_on = timestep_end
[]
[]
[BCs]
# Define boundary conditions
[no_y_all] # pin pellets and clad along axis of symmetry (y)
type = DirichletBC
variable = disp_y
boundary = 15
value = 0.0
[]
[no_x_all] # pin pellets and clad along axis of symmetry (x)
type = DirichletBC
variable = disp_x
boundary = 16
value = 0.0
[]
[no_y_all_grid] # pin pellets and clad along axis of symmetry (y)
type = FunctionDirichletBC
variable = disp_y
boundary = '112'
function = 'if(t < 1.0e4,1.0e-4 * t/1.0e4 - 1.0e-5,0.9e-4)'
[]
[no_x_all_grid] # pin pellets and clad along axis of symmetry (x)
type = DirichletBC
variable = disp_x
boundary = '112'
value = 0.0
[]
[Pressure] # apply coolant pressure on clad outer walls
[coolantPressure]
boundary = '2'
factor = 15.5e6
function = pressure_var # use the pressure_ramp function defined above
[]
[]
[PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
[plenumPressure]
boundary = 9
initial_pressure = 2.0e6
R = 8.3143
output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
temperature = plenum_temperature # coupling to post processor to get gas temperature approximation
volume = plenum_volume # coupling to post processor to get gas volume
material_input = fission_gas_released # coupling to post processor to get fission gas added
output = plenum_pressure # coupling to post processor to output plenum/gap pressure
displacements = 'disp_x disp_y'
[]
[]
[convective_clad_surface] # apply convective boundary to clad outer surface
type = ConvectiveFluxBC
boundary = '2'
variable = temp
rate = 38200.0 #convection coefficient (h)
initial = 580.0
final = 580.0
duration = 1.0e4 #duration of initial power ramp
[]
[]
[Materials]
# Define material behavior models and input material property data
[fuel_thermal] # temperature and burnup dependent thermal properties of UO2 (BISON kernel)
type = UO2Thermal
thermal_conductivity_model = FINK_LUCUTA
block = pellet_type_1
temperature = temp
burnup = burnup
initial_porosity = 0.0
[]
[fuel_solid_mechanics_swelling] # free expansion strains (swelling and densification) for UO2 (BISON kernel)
type = UO2VolumetricSwellingEigenstrain
gas_swelling_model_type = MATPRO
block = pellet_type_1
burnup = burnup
initial_fuel_density = 10431.0
temperature = temp
eigenstrain_name = 'fuel_volumetric_eigenstrain'
[]
[fuel_creep]
type = UO2CreepUpdate
block = pellet_type_1
temperature = temp
fission_rate = fission_rate
density = 10431.0
initial_grain_radius = 10.0e-6
oxygen_to_metal_ratio = 2.0
[]
[fuel_elasticity_tensor]
type = ComputeIsotropicElasticityTensor
block = 'pellet_type_1'
youngs_modulus = 906e6
poissons_ratio = 0.345
[]
[fuel_stress]
type = ComputeMultipleInelasticStress
block = pellet_type_1
inelastic_models = 'fuel_creep'
[]
[fuel_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = pellet_type_1
thermal_expansion_coeff = 10.0e-6
temperature = temp
stress_free_temperature = 580.0
eigenstrain_name = 'fuel_thermal_eigenstrain'
[]
[fuel_relocation]
type = UO2RelocationEigenstrain
block = pellet_type_1
burnup = burnup
diameter = 0.0082
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
diametral_gap =160e-6
burnup_relocation_stop = 1.e20
relocation_activation1 = 5000
axial_axis = 2
eigenstrain_name = 'fuel_relocation_eigenstrain'
[]
[clad_thermal]
type = HeatConductionMaterial
block = 'clad'
thermal_conductivity = 16.0
specific_heat = 330.0
[]
[clad_elasticity_tensor]
type = ZryElasticityTensor
block = 'clad'
[]
[clad_creep_model]
type = ZryCreepHayesHoppeUpdate
block = 'clad'
fast_neutron_flux = fast_neutron_flux
temperature = temp
zircaloy_material_type = stress_relief_annealed
model_irradiation_creep = true
model_thermal_creep = true
[]
[clad_stress]
type = ComputeMultipleInelasticStress
block = 'clad'
tangent_operator = elastic
inelastic_models = 'clad_creep_model'
[]
[clad_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = 'clad'
thermal_expansion_coeff = 5.0e-6
temperature = temp
stress_free_temperature = 580.0
eigenstrain_name = 'clad_thermal_eigenstrain'
[]
[clad_irrgrowth]
type = ZryIrradiationGrowthEigenstrain
block = clad
fast_neutron_fluence = fast_neutron_fluence
axial_direction = 2
zircaloy_material_type = ESCORE_IrradiationGrowthZr4
eigenstrain_name = 'clad_irradiation_eigenstrain'
[]
[grid_thermal]
type = HeatConductionMaterial
block = 'grid'
thermal_conductivity = 16.0
specific_heat = 330.0
[]
[grid_elasticity_tensor]
type = ZryElasticityTensor
block = 'grid'
[]
[grid_creep_model]
type = ZryCreepHayesHoppeUpdate
block = 'grid'
fast_neutron_flux = fast_neutron_flux
temperature = temp
zircaloy_material_type = stress_relief_annealed
model_irradiation_creep = true
model_thermal_creep = true
[]
[grid_stress]
type = ComputeMultipleInelasticStress
block = 'grid'
tangent_operator = elastic
inelastic_models = 'grid_creep_model'
[]
[grid_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = 'grid'
thermal_expansion_coeff = 5.0e-6
temperature = temp
stress_free_temperature = 580.0
eigenstrain_name = 'grid_thermal_eigenstrain'
[]
[grid_irrgrowth]
type = ZryIrradiationGrowthEigenstrain
block = grid
fast_neutron_fluence = fast_neutron_fluence
axial_direction = 2
zircaloy_material_type = ESCORE_IrradiationGrowthZr4
eigenstrain_name = 'grid_irradiation_eigenstrain'
[]
[fission_gas_release] # Forsberg-Massih fission gas release mode
type = UO2Sifgrs
block = pellet_type_1
temperature = temp
fission_rate = fission_rate # coupling to fission_rate aux variable
grain_radius = 10.0e-6
#external_pressure = 40e6
[]
[clad_density]
type = StrainAdjustedDensity
block = 'clad'
density = 6551.0
[]
[fuel_density]
type = StrainAdjustedDensity
block = pellet
strain_free_density = 10431.0
[]
[grid]
type = StrainAdjustedDensity
block = grid
density = 6560
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason -ksp_converged_reason'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
petsc_options_value = 'lu superlu_dist 1e-6 NONZERO 1e-10'
snesmf_reuse_base = true
line_search = 'none'
l_max_its = 100
l_tol = 8e-3
nl_max_its = 45
nl_rel_tol = 1e-10 # was -7 and nl 25. Tightening tangential contact forces.
nl_abs_tol = 1e-12
[TimeIntegrator]
type = NewmarkBeta
beta = 0.25
gamma = 0.5
[]
start_time = 0.0
end_time = 1.0e5
[TimeStepper]
type = IterationAdaptiveDT
dt = 2.0e2
time_t = '1e4 1e5 1e6'
time_dt = '2e2 1e4 1e5'
growth_factor = 1.4
iteration_window = 5.0
optimal_iterations = 35
[]
dtmax = 2e5 # Larger causes instabilities 2e6
dtmin = 1
[]
[Postprocessors]
# Define postprocessors (some are required as specified above; others are optional; many others are available)
[average_interior_clad_temperature] # average temperature of cladding interior
type = SideAverageValue
boundary = 7
variable = temp
execute_on = 'initial timestep_end'
[]
[average_centerline_fuel_temperature] # average temperature of the cladding interior and all pellet exteriors
type = SideAverageValue
boundary = 9
variable = temp
execute_on = 'initial linear'
[]
[plenum_temperature]
type = SideAverageValue
boundary = 9
variable = temp
execute_on = 'initial timestep_end'
[]
[plenum_volume] # gas volume
type = InternalVolume
boundary = 9
addition = 1.3e-5 #rough guess of plenum volume/unit length of fuel
execute_on = 'initial linear'
[]
[pellet_volume] # fuel pellet total volume
type = InternalVolume
boundary = 8
execute_on = 'initial timestep_end'
[]
[clad_inner_vol] # volume inside of cladding
type = InternalVolume
boundary = 7
outputs = exodus
execute_on = 'initial timestep_end'
[]
[fission_gas_generated] # fission gas produced (moles)
type = ElementIntegralFisGasGeneratedSifgrs
block = pellet_type_1
execute_on = linear
[]
[fission_gas_released] # fission gas released to plenum (moles)
type = ElementIntegralFisGasReleasedSifgrs
block = pellet_type_1
execute_on = linear
[]
[flux_from_clad] # area integrated heat flux from the cladding
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 5
diffusivity = thermal_conductivity
execute_on = timestep_end
[]
[flux_from_fuel] # area integrated heat flux from the fuel
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 10
diffusivity = thermal_conductivity
execute_on = timestep_end
[]
[_dt] # time step
type = TimestepSize
execute_on = timestep_end
[]
[num_lin_it]
type = NumLinearIterations
[]
[num_nonlin_it]
type = NumNonlinearIterations
[]
[tot_lin_it]
type = CumulativeValuePostprocessor
postprocessor = num_lin_it
[]
[tot_nonlin_it]
type = CumulativeValuePostprocessor
postprocessor = num_nonlin_it
[]
[alive_time]
type = PerfGraphData
section_name = Root
data_type = TOTAL
[]
[rod_total_power]
type = ElementIntegralPower
variable = temp
fission_rate = fission_rate
block = pellet_type_1
execute_on = timestep_end
[]
[rod_input_power]
type = FunctionValuePostprocessor
function = power_history
scale_factor = 0.1186 # rod height
execute_on = timestep_end
[]
[fission_gas_released_percentage]
type = FGRPercent
fission_gas_released = fission_gas_released
fission_gas_generated = fission_gas_generated
[]
[]
[VectorPostprocessors]
[contact_pressure]
type = NodalValueSampler
sort_by = x
use_displaced_mesh = true
variable = spacer_clad_mechanical_normal_lm
boundary = 102
[]
[frictional_pressure]
type = NodalValueSampler
sort_by = x
use_displaced_mesh = true
variable = spacer_clad_mechanical_tangential_lm
boundary = 102
[]
[worn_depth]
type = NodalValueSampler
sort_by = x
use_displaced_mesh = true
variable = worn_depth
boundary = 102
execute_on = TIMESTEP_END
[]
[]
[Outputs]
perf_graph = true
exodus = true
csv = true
[console]
type = Console
max_rows = 25
[]
checkpoint = true
[]
(test/tests/example_problem_test/example_problem_test.i)
[GlobalParams]
density = 10431.0
displacements = 'disp_x disp_y'
energy_per_fission = 3.2e-11 # J/fission
temperature = temp
volumetric_locking_correction = false
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
[]
[Mesh]
coord_type = RZ
patch_update_strategy = auto
patch_size = 10 # For contact algorithm
partitioner = centroid
centroid_partitioner_direction = y
[mesh]
type = FileMeshGenerator
file = 2_pellet_discrete.e
[]
[]
[Variables]
[temp]
initial_condition = 580.0
[]
[]
[AuxVariables]
[fast_neutron_flux]
block = clad
[]
[fast_neutron_fluence]
block = clad
[]
[grain_radius]
block = pellet_type_1
initial_condition = 10e-6
[]
[gap_cond]
order = CONSTANT
family = MONOMIAL
[]
[coolant_htc]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear
data_file = powerhistory.csv
scale_factor = 1
[]
[axial_peaking_factors]
type = ParsedFunction
expression = 1
[]
[pressure_ramp]
type = PiecewiseLinear
x = '-200 0'
y = '0 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[fuel]
block = pellet_type_1
strain = FINITE
incremental = true
extra_vector_tags = 'ref'
add_variables = true
decomposition_method = EigenSolution
eigenstrain_names = 'fuel_volumetric_swelling_eigenstrain
fuel_relocation_eigenstrain fuel_thermal_eigenstrain'
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress'
[]
[clad]
block = clad
strain = FINITE
incremental = true
extra_vector_tags = 'ref'
add_variables = true
decomposition_method = EigenSolution
eigenstrain_names = 'clad_thermal_strain clad_irradiation_growth_eigenstrain'
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress'
[]
[]
[Kernels]
[gravity]
type = Gravity
variable = disp_y
value = -9.81
[]
[heat]
type = HeatConduction
variable = temp
extra_vector_tags = 'ref'
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
extra_vector_tags = 'ref'
[]
[heat_source]
type = NeutronHeatSource
variable = temp
extra_vector_tags = 'ref'
block = pellet_type_1
burnup_function = burnup
[]
[]
[UserObjects]
[pin_geometry]
type = FuelPinGeometry
[]
[]
[Burnup]
[burnup]
block = pellet_type_1
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
num_radial = 80
num_axial = 11
fuel_pin_geometry = 'pin_geometry'
fuel_volume_ratio = 0.987775
order = CONSTANT
family = MONOMIAL
RPF = RPF
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
block = clad
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
factor = 3e13
execute_on = timestep_begin
[]
[fast_neutron_fluence]
type = FastNeutronFluenceAux
variable = fast_neutron_fluence
block = clad
fast_neutron_flux = fast_neutron_flux
execute_on = timestep_begin
[]
[grain_radius]
type = GrainRadiusAux
block = pellet_type_1
variable = grain_radius
temperature = temp
execute_on = linear
[]
[conductance]
type = MaterialRealAux
property = gap_conductance
variable = gap_cond
boundary = 10
execute_on = 'initial timestep_end'
[]
[coolant_htc]
type = MaterialRealAux
property = coolant_channel_htc
variable = coolant_htc
boundary = 2
execute_on = 'initial timestep_end'
[]
[]
[Contact]
[pellet_clad_mechanical]
primary = 5
secondary = 10
formulation = KINEMATIC
model = frictionless
normalize_penalty = true
penalty = 1e14
[]
[]
[ThermalContact]
[thermal_contact]
type = GasGapHeatTransfer
variable = temp
primary = 5
secondary = 10
initial_moles = initial_moles
gas_released = fis_gas_released
contact_pressure = contact_pressure
[]
[]
[BCs]
[no_x_all]
type = DirichletBC
variable = disp_x
boundary = 12
value = 0.0
[]
[no_y_clad_bottom]
type = DirichletBC
variable = disp_y
boundary = '1'
value = 0.0
[]
[no_y_fuel_bottom]
type = DirichletBC
variable = disp_y
boundary = '1020'
value = 0.0
[]
[Pressure]
[coolantPressure]
boundary = '1 2 3'
factor = 15.5e6
function = pressure_ramp
[]
[]
[PlenumPressure]
[plenumPressure]
boundary = 9
initial_pressure = 2.0e6
startup_time = -200
R = 8.3143
output_initial_moles = initial_moles
temperature = ave_temp_interior
volume = gas_volume
material_input = fis_gas_released
output = plenum_pressure
displacements = 'disp_x disp_y'
execute_on = 'initial linear'
[]
[]
[]
[CoolantChannel]
[convective_clad_surface]
boundary = '1 2 3'
variable = temp
inlet_temperature = 580 # K
inlet_pressure = 15.5e6 # Pa
inlet_massflux = 3800 # kg/m^2-sec
rod_diameter = 0.948e-2 # m
rod_pitch = 1.26e-2 # m
linear_heat_rate = power_history
axial_power_profile = axial_peaking_factors
[]
[]
[Materials]
[fuel_thermal]
type = UO2Thermal
block = pellet_type_1
thermal_conductivity_model = NFIR
initial_porosity = 0.0
temperature = temp
burnup_function = burnup
[]
[fuel_volumetric_swelling]
type = UO2VolumetricSwellingEigenstrain
block = pellet_type_1
burnup = burnup
initial_fuel_density = 10431.0
eigenstrain_name = 'fuel_volumetric_swelling_eigenstrain'
[]
[fuel_elasticity_tensor]
type = UO2ElasticityTensor
block = pellet_type_1
[]
[fuel_thermal_expansion]
type = UO2ThermalExpansionMartinEigenstrain
block = pellet_type_1
stress_free_temperature = 295
eigenstrain_name = 'fuel_thermal_eigenstrain'
[]
[hotpressing]
type = UO2HotPressingCreepUpdate
block = pellet_type_1
burnup_function = burnup
initial_grain_radius = 10.0e-6
[]
[radial_return_stress]
type = ComputeMultipleInelasticStress
tangent_operator = elastic
inelastic_models = ' hotpressing'
block = pellet_type_1
[]
[fuel_relocation]
type = UO2RelocationEigenstrain
block = pellet_type_1
burnup_function = burnup
fuel_pin_geometry = 'pin_geometry'
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
relocation_activation1 = 5000 #TM default value
burnup_relocation_stop = 1.e20
eigenstrain_name = 'fuel_relocation_eigenstrain'
[]
[clad_thermal]
type = HeatConductionMaterial
block = clad
thermal_conductivity = 16.0
specific_heat = 330.0
[]
[clad_elasticity_tensor]
type = ZryElasticityTensor
block = clad
[]
[clad_creep_model]
type = ZryCreepHayesHoppeUpdate
block = clad
fast_neutron_flux = fast_neutron_flux
model_irradiation_creep = true
model_thermal_creep = true
[]
[clad_inelastic_stress]
type = ComputeMultipleInelasticStress
block = clad
tangent_operator = elastic
inelastic_models = 'clad_creep_model'
[]
[clad_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = clad
thermal_expansion_coeff = 5.0e-6
stress_free_temperature = 295.0
eigenstrain_name = clad_thermal_strain
[]
[clad_irradiation_growth]
type = ZryIrradiationGrowthEigenstrain
block = clad
fast_neutron_fluence = fast_neutron_fluence
zircaloy_material_type = ESCORE_IrradiationGrowthZr4
eigenstrain_name = clad_irradiation_growth_eigenstrain
[]
[fission_gas_release]
type = UO2Sifgrs
block = pellet_type_1
temperature = temp
burnup_function = burnup
grain_radius = grain_radius
gbs_model = true
[]
[clad_density]
type = StrainAdjustedDensity
block = clad
strain_free_density = 6551.0
[]
[fuel_density]
type = StrainAdjustedDensity
block = pellet_type_1
strain_free_density = 10431.0
[]
[]
[Dampers]
[BoundingValueNodalDamper]
type = BoundingValueNodalDamper
variable = temp
max_value = 3200
min_value = 300
[]
[]
[Preconditioning]
[SMP]
type = SMP
coupled_groups = 'disp_x,disp_y'
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-pc_type_asm'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
petsc_options_value = 'lu superlu_dist 51'
line_search = 'none'
verbose = true
l_max_its = 100
l_tol = 1e-5 #8e-3
nl_max_its = 15
nl_rel_tol = 1e-10
nl_abs_tol = 1e-8
start_time = -200
num_steps = 2
dtmax = 2e6
dtmin = 1
[TimeStepper]
type = IterationAdaptiveDT
dt = 2.0e2
optimal_iterations = 6
iteration_window = 2
[]
[]
[Postprocessors]
[ave_temp_interior] # average temperature of the cladding interior and all pellet exteriors
type = SideAverageValue
boundary = 9
variable = temp
execute_on = 'initial linear'
[]
[clad_inner_vol] # volume inside of cladding
type = InternalVolume
boundary = 7
outputs = exodus
execute_on = 'initial timestep_end'
[]
[pellet_volume] # fuel pellet total volume
type = InternalVolume
boundary = 8
outputs = exodus
execute_on = 'initial timestep_end'
[]
[avg_clad_temp] # average temperature of cladding interior
type = SideAverageValue
boundary = 7
variable = temp
execute_on = 'initial timestep_end'
[]
[fis_gas_produced] # fission gas produced (moles)
type = ElementIntegralFisGasGeneratedSifgrs
block = pellet_type_1
execute_on = timestep_end
[]
[fis_gas_released] # fission gas released to plenum (moles)
type = ElementIntegralFisGasReleasedSifgrs
block = pellet_type_1
execute_on = timestep_end
[]
[fis_gas_grain]
type = ElementIntegralFisGasGrainSifgrs
block = pellet_type_1
outputs = exodus
execute_on = timestep_end
[]
[fis_gas_boundary]
type = ElementIntegralFisGasBoundarySifgrs
block = pellet_type_1
outputs = exodus
execute_on = timestep_end
[]
[gas_volume] # gas volume
type = InternalVolume
boundary = 9
component = 1
execute_on = 'initial linear'
[]
[flux_from_clad] # area integrated heat flux from the cladding
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 5
diffusivity = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[flux_from_fuel] # area integrated heat flux from the fuel
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 10
diffusivity = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[_dt] # time step
type = TimestepSize
execute_on = timestep_end
[]
[nonlinear_its] # number of nonlinear iterations at each timestep
type = NumNonlinearIterations
execute_on = timestep_end
[]
[rod_total_power]
type = ElementIntegralPower
variable = temp
burnup_function = burnup
block = pellet_type_1
execute_on = 'initial timestep_end'
[]
[rod_input_power]
type = FunctionValuePostprocessor
function = power_history
scale_factor = 0.02372 # rod height
execute_on = 'initial timestep_end'
[]
[]
[Outputs]
exodus = true
color = false
[console]
type = Console
output_linear = true
max_rows = 25
[]
[]
(test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep_rz.i)
#
# The mesh is a 1x1 square with a constant pressure of 10 MPa on the top face.
# Symmetry boundary conditions on three planes provide a uniaxial stress
# field. The temperature is held constant at 1000; only creep is modeled: no
# platicity material model is included in the material blocks. The solution is
# advanced through ten time steps of 1e4 for a total time of 1e5.
#
# The total strain at time 1e5 can be computed as:
#
# e_tot = e_elas +
# e_thermal_creep +
# e_irradiation_creep
#
# = P/E +
# A * (sigma/G)**n * exp(-Q/(RT)) * dt +
# C0 * (phi)**C1 * (sigma/10**6)**C2 * dt
#
# where P = pressure load
# E = Young's modulus
# A = thermal creep coefficient
# G = shear modulus = 0.5*E/(1+nu)
# nu = Poisson's ratio
# sigma = stress
# n = thermal creep exponent
# Q = activation energy
# R = gas constant
# T = temperature
# C0 = irradiation creep coefficient
# phi = fast neutron flux
# C1 = neutron flux exponent
# C2 = stress exponent
# dt = total time
#
# For this test, the analytical solutuon at t = 1e5 is:
#
# G = 0.5*2e11/(1+0.3) = 7.692308e10
#
# e_tot = 10e6/2e11 +
# 3.14e24 * (10e6/G)**5 * exp(-2.7e5/(8.3143*1000) * 1e5
# 9.881e-28 * (1e19)**0.85 * (10e6/1e6)**1 * 1e5
#
# = 5e-5 +
# 9.189527e-5
# 1.395728e-5
#
# = 1.558526e-4
#
# For both small strain (ComputeIncrementalSmallStrain) or finite strain (ComputeFiniteStrain)
# kinematics, BISON computes:
#
# e_elas = 5e-5
# e_creep = 1.05819e-4
# e_tot = 1.55819e-4
#
#
[Mesh]
coord_type = RZ
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Variables]
[disp_x]
order = FIRST
family = LAGRANGE
[]
[disp_y]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1000.0
[]
[]
[AuxVariables]
[fast_neutron_flux]
order = FIRST
family = LAGRANGE
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[elastic_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[creep_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[top_pull]
type = PiecewiseLinear
x = '0 1'
y = '1 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = finite
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
factor = 1e19
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[elastic_strain_yy]
type = RankTwoAux
rank_two_tensor = elastic_strain
variable = elastic_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[creep_strain_yy]
type = RankTwoAux
rank_two_tensor = creep_strain
variable = creep_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[]
[BCs]
[Pressure]
[u_top_pull]
boundary = top
factor = -10.0e6
function = top_pull
[]
[]
[u_bottom_fix]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[]
[temp_fix]
type = DirichletBC
variable = temp
boundary = 'top bottom left right'
value = 1000.0
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 2.0e11
poissons_ratio = 0.3
[]
[stress]
type = ComputeMultipleInelasticStress
tangent_operator = elastic
inelastic_models = 'zrycreep'
[]
[zrycreep]
type = ZryCreepHayesHoppeUpdate
temperature = temp
fast_neutron_flux = fast_neutron_flux
model_irradiation_creep = true
model_thermal_creep = true
[]
[thermal]
type = HeatConductionMaterial
block = 0
specific_heat = 1.0
thermal_conductivity = 100.
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 1.0
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart'
petsc_options_value = '101'
line_search = 'none'
l_max_its = 100
nl_max_its = 100
nl_rel_tol = 1e-8
nl_abs_tol = 1e-6
l_tol = 1e-5
start_time = 0.0
num_steps = 10
dt = 10000
[]
[Outputs]
show = 'disp_x disp_y temp stress_yy elastic_strain_yy creep_strain_yy'
[out]
type = Exodus
[]
[]
(examples/2D_plane_strain_fretting_wear/fretting-wear-initial-dyn-exc.i)
user_start_time = 1.0e5
user_end_time = 1.000002e5
end_dynamic_excitation = 1.000002e5
time_step_dynamics = 2.0e-3
step_number = 1
initial_fuel_density = 10431.0
[GlobalParams]
temperature = temp
displacements = 'disp_x disp_y'
order = FIRST
family = LAGRANGE
energy_per_fission = 3.2e-11 # J/fission
volumetric_locking_correction = true
[]
[Mesh]
[file]
type = FileMeshGenerator
file = fretting-wear-initial_out_cp/LATEST
skip_partitioning = true
allow_renumbering = false
[]
patch_size = 100 # For contact algorithm
[]
[Variables]
[temp]
[]
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
group_variables = 'disp_x disp_y '
converge_on = 'disp_x disp_y temp'
restart_file_base = ./fretting-wear-initial_out_cp/LATEST
material_coverage_check = false
kernel_coverage_check = false
[]
[AuxVariables]
[fission_rate]
block = pellet_type_1
[]
[burnup]
block = pellet_type_1
[]
[fast_neutron_flux]
block = 'clad grid'
[]
[fast_neutron_fluence]
block = 'clad grid'
[]
[relocation_strain]
order = CONSTANT
family = MONOMIAL
[]
[worn_depth]
order = FIRST
family = LAGRANGE
block = 'spacer_clad_mechanical_secondary_subdomain'
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear # reads and interpolates an input file containing rod average linear power vs time
data_file = powerhistory.csv
scale_factor = 1
[]
[axial_peaking_factors]
type = ConstantFunction
value = 1
[]
[pressure_var] # reads and interpolates input data defining amplitude curve for fill gas pressure
type = PiecewiseLinear
x = '0 1e4'
y = '0 1'
[]
[pressure_var_variable] # reads and interpolates input data defining amplitude curve for fill gas pressure
type = ParsedFunction
expression = 'if(t < 1e4, 1, 1 + sin((t-1e4)*pi/10.0) * (t-1e4))'
[]
[]
[Physics/SolidMechanics/Dynamic]
[pellets]
add_variables = true
newmark_beta = 0.25
newmark_gamma = 0.5
block = pellet_type_1
strain = FINITE
planar_formulation = PLANE_STRAIN
eigenstrain_names = 'fuel_relocation_eigenstrain fuel_thermal_eigenstrain
fuel_volumetric_eigenstrain'
generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
decomposition_method = EigenSolution
temperature = temp
extra_vector_tags = 'ref'
[]
[clad]
add_variables = true
newmark_beta = 0.25
newmark_gamma = 0.5
block = clad
strain = FINITE
planar_formulation = PLANE_STRAIN
eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_eigenstrain'
generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
decomposition_method = EigenSolution
temperature = temp
extra_vector_tags = 'ref'
[]
[grid]
add_variables = true
newmark_beta = 0.25
newmark_gamma = 0.5
block = grid
strain = FINITE
planar_formulation = PLANE_STRAIN
eigenstrain_names = 'grid_thermal_eigenstrain grid_irradiation_eigenstrain'
generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
decomposition_method = EigenSolution
temperature = temp
extra_vector_tags = 'ref'
[]
[]
[Kernels]
[heat] # gradient term in heat conduction equation
type = HeatConduction
variable = temp
block = 'pellet_type_1 clad grid'
extra_vector_tags = 'ref'
[]
[heat_ie] # time term in heat conduction equation
type = HeatConductionTimeDerivative
variable = temp
block = 'pellet_type_1 clad'
extra_vector_tags = 'ref'
[]
[heat_source] # source term in heat conduction equation
type = NeutronHeatSource
variable = temp
block = pellet_type_1
fission_rate = fission_rate
extra_vector_tags = 'ref'
[]
[]
[Contact]
# Define mechanical contact between the fuel (sideset=10) and the clad (sideset=5)
[spacer_clad_mechanical]
formulation = mortar
model = coulomb
primary = 101
secondary = 102
c_normal = 1e+12 # 5e13
c_tangential = 1e+18
friction_coefficient = 0.4
# Do not apply dynamic stabilization
newmark_beta = 0.0001
newmark_gamma = 0.5
capture_tolerance = 0.0
mortar_dynamics = true
interpolate_normals = false
normal_lm_scaling = 1.0e-6
tangential_lm_scaling = 1.0e-6
generate_mortar_mesh = false
wear_depth = worn_depth
[]
[]
[Contact]
# Define mechanical contact between the fuel (sideset=10) and the clad (sideset=5)
[pellet_clad_mechanical_real]
formulation = mortar
model = frictionless
primary = 7
secondary = 8
c_normal = 1e+16 #
c_tangential = 1e+16
friction_coefficient = 0.4
# Do not apply dynamic stabilization
newmark_beta = 0.0001
newmark_gamma = 0.5
capture_tolerance = 0.0
mortar_dynamics = true
interpolate_normals = false
generate_mortar_mesh = false
[]
[]
[ThermalContactMortar]
[thermal_contact]
secondary_variable = temp
primary_boundary = 7
secondary_boundary = 8
initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
gas_released = fission_gas_released # coupling to a postprocessor which supplies the fission gas addition
primary_subdomain = 'pellet_clad_mechanical_real_primary_subdomain'
secondary_subdomain = 'pellet_clad_mechanical_real_secondary_subdomain'
[]
[]
[Burnup]
[burnup]
block = pellet_type_1
rod_ave_lin_pow = power_history # using the power function defined above
axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
num_radial = 80
num_axial = 21
axial_axis = 2
density = ${initial_fuel_density}
a_lower = -1e-3 # mesh dependent!
a_upper = 1e-3 # mesh dependent!
fuel_inner_radius = 0
fuel_outer_radius = .0041
fuel_volume_ratio = 0.987775 # for use with dished pellets (ratio of actual volume to cylinder volume)
#N235 = N235 # Activate to write N235 concentration to output file
#N238 = N238 # Activate to write N238 concentration to output file
#N239 = N239 # Activate to write N239 concentration to output file
#N240 = N240 # Activate to write N240 concentration to output file
#N241 = N241 # Activate to write N241 concentration to output file
#N242 = N242 # Activate to write N242 concentration to output file
RPF = RPF
[]
[]
[AuxKernels]
# Define auxilliary kernels for each of the aux variables
[worn_depth]
type = MortarArchardsLawAux
variable = worn_depth
primary_boundary = 101
secondary_boundary = 102
primary_subdomain = 'spacer_clad_mechanical_primary_subdomain'
secondary_subdomain = 'spacer_clad_mechanical_secondary_subdomain'
displacements = 'disp_x disp_y'
friction_coefficient = 0.5
energy_wear_coefficient = 0.1e-9
normal_pressure = spacer_clad_mechanical_normal_lm
execute_on = 'TIMESTEP_END'
[]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
block = clad
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
factor = 3e13
execute_on = timestep_begin
[]
[fast_neutron_fluence]
type = FastNeutronFluenceAux
variable = fast_neutron_fluence
block = clad
fast_neutron_flux = fast_neutron_flux
execute_on = timestep_begin
[]
[relocation_strain]
type = MaterialRealAux
property = relocation_strain
variable = relocation_strain
block = pellet_type_1
execute_on = timestep_end
[]
[]
[BCs]
# Define boundary conditions
[no_y_all] # pin pellets and clad along axis of symmetry (y)
type = DirichletBC
variable = disp_y
boundary = 15
value = 0.0
[]
[no_x_all] # pin pellets and clad along axis of symmetry (x)
type = DirichletBC
variable = disp_x
boundary = 16
value = 0.0
[]
# Flow-induced vibrations refined_excitation
[vibration_x] # pin pellets and clad along axis of symmetry (y)
type = FunctionDirichletBC
variable = disp_x
boundary = '112'
expression = 'if(t < ${end_dynamic_excitation}, 10.0*1.0e-6*sin(2*3.1415926535*20* (t - ${user_start_time})) + 2.0*1.0e-6*sin(2*3.1415926535*35*(t - ${user_start_time})), 0)'
#expression = '0'
[]
[vibration_y] # pin pellets and clad along axis of symmetry (y)
type = FunctionDirichletBC
variable = disp_y
boundary = '112'
expression = 'if(t < ${end_dynamic_excitation}, 10.0*1.0e-6*sin(2*3.1415926535*20*(t-${user_start_time})) + 2.0*1.0e-6*sin(2*3.1415926535*35*(t-${user_start_time})) + 0.9e-4, 0.9e-4)'
#expression = '5.9e-4'
[]
[Pressure] # apply coolant pressure on clad outer walls
[coolantPressure]
boundary = '2'
factor = 15.5e6
function = pressure_var # use the pressure_ramp function defined above
[]
[]
[PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
[plenumPressure]
boundary = 9
initial_pressure = 2.0e6
R = 8.3143
output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
temperature = plenum_temperature # coupling to post processor to get gas temperature approximation
volume = plenum_volume # coupling to post processor to get gas volume
material_input = fission_gas_released # coupling to post processor to get fission gas added
output = plenum_pressure # coupling to post processor to output plenum/gap pressure
displacements = 'disp_x disp_y'
[]
[]
[convective_clad_surface] # apply convective boundary to clad outer surface
type = ConvectiveFluxBC
boundary = '2'
variable = temp
rate = 38200.0 #convection coefficient (h)
initial = 580.0
final = 580.0
duration = 1.0e4 #duration of initial power ramp
[]
[]
[Materials]
# Define material behavior models and input material property data
[fuel_thermal] # temperature and burnup dependent thermal properties of UO2 (BISON kernel)
type = UO2Thermal
thermal_conductivity_model = FINK_LUCUTA
block = pellet_type_1
temperature = temp
burnup = burnup
initial_porosity = 0.0
[]
[fuel_solid_mechanics_swelling] # free expansion strains (swelling and densification) for UO2 (BISON kernel)
type = UO2VolumetricSwellingEigenstrain
gas_swelling_model_type = MATPRO
block = pellet_type_1
burnup = burnup
initial_fuel_density = 10431.0
temperature = temp
eigenstrain_name = 'fuel_volumetric_eigenstrain'
[]
[fuel_creep]
type = UO2CreepUpdate
block = pellet_type_1
temperature = temp
fission_rate = fission_rate
density = 10431.0
initial_grain_radius = 10.0e-6
oxygen_to_metal_ratio = 2.0
[]
[fuel_elasticity_tensor]
type = ComputeIsotropicElasticityTensor
block = 'pellet_type_1'
youngs_modulus = 906e6
poissons_ratio = 0.345
[]
[fuel_stress]
type = ComputeMultipleInelasticStress
block = pellet_type_1
inelastic_models = 'fuel_creep'
[]
[fuel_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = pellet_type_1
thermal_expansion_coeff = 10.0e-6
temperature = temp
stress_free_temperature = 580.0
eigenstrain_name = 'fuel_thermal_eigenstrain'
[]
[fuel_relocation]
type = UO2RelocationEigenstrain
block = pellet_type_1
burnup = burnup
diameter = 0.0082
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
diametral_gap =160e-6
burnup_relocation_stop = 1.e20
relocation_activation1 = 5000
axial_axis = 2
eigenstrain_name = 'fuel_relocation_eigenstrain'
[]
[clad_thermal]
type = HeatConductionMaterial
block = 'clad'
thermal_conductivity = 16.0
specific_heat = 330.0
[]
[clad_elasticity_tensor]
type = ZryElasticityTensor
block = clad
[]
[clad_creep_model]
type = ZryCreepHayesHoppeUpdate
block = clad
fast_neutron_flux = fast_neutron_flux
temperature = temp
zircaloy_material_type = stress_relief_annealed
model_irradiation_creep = true
model_thermal_creep = true
[]
[clad_stress]
type = ComputeMultipleInelasticStress
block = clad
tangent_operator = elastic
inelastic_models = 'clad_creep_model'
[]
[clad_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = clad
thermal_expansion_coeff = 5.0e-6
temperature = temp
stress_free_temperature = 580.0
eigenstrain_name = 'clad_thermal_eigenstrain'
[]
[clad_irrgrowth]
type = ZryIrradiationGrowthEigenstrain
block = clad
fast_neutron_fluence = fast_neutron_fluence
axial_direction = 2
zircaloy_material_type = ESCORE_IrradiationGrowthZr4
eigenstrain_name = 'clad_irradiation_eigenstrain'
[]
[grid_thermal]
type = HeatConductionMaterial
block = 'grid'
thermal_conductivity = 16.0
specific_heat = 330.0
[]
[grid_elasticity_tensor]
type = ZryElasticityTensor
block = 'grid'
[]
[grid_creep_model]
type = ZryCreepHayesHoppeUpdate
block = 'grid'
fast_neutron_flux = fast_neutron_flux
temperature = temp
zircaloy_material_type = stress_relief_annealed
model_irradiation_creep = true
model_thermal_creep = true
[]
[grid_stress]
type = ComputeMultipleInelasticStress
block = 'grid'
tangent_operator = elastic
inelastic_models = 'grid_creep_model'
[]
[grid_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = 'grid'
thermal_expansion_coeff = 5.0e-6
temperature = temp
stress_free_temperature = 580.0
eigenstrain_name = 'grid_thermal_eigenstrain'
[]
[grid_irrgrowth]
type = ZryIrradiationGrowthEigenstrain
block = grid
fast_neutron_fluence = fast_neutron_fluence
axial_direction = 2
zircaloy_material_type = ESCORE_IrradiationGrowthZr4
eigenstrain_name = 'grid_irradiation_eigenstrain'
[]
[fission_gas_release] # Forsberg-Massih fission gas release mode
type = UO2Sifgrs
block = pellet_type_1
temperature = temp
fission_rate = fission_rate # coupling to fission_rate aux variable
grain_radius = 10.0e-6
#external_pressure = 40e6
[]
[clad_density]
type = StrainAdjustedDensity
block = 'clad'
density = 6551.0
[]
[fuel_density]
type = StrainAdjustedDensity
block = pellet
strain_free_density = 10431.0
[]
[grid]
type = StrainAdjustedDensity
block = grid
density = 6560
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason -ksp_converged_reason'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
petsc_options_value = 'lu superlu_dist 1e-6 NONZERO 1e-14'
snesmf_reuse_base = true
line_search = 'basic'
l_max_its = 100
l_tol = 8e-3
nl_max_its = 35
nl_rel_tol = 1e-7
nl_abs_tol = 1e-11
[TimeIntegrator]
type = NewmarkBeta
beta = 0.25
gamma = 0.5
[]
start_time = '${user_start_time}'
end_time = '${user_end_time}'
timestep_tolerance = 1e-8
[TimeStepper]
type = IterationAdaptiveDT
dt = '${time_step_dynamics}'
time_t = '${end_dynamic_excitation}'
time_dt = '${time_step_dynamics}'
growth_factor = 1.2
cutback_factor = 0.75
[]
dtmax = 3e-2
[]
[Postprocessors]
# Define postprocessors (some are required as specified above; others are optional; many others are available)
[average_interior_clad_temperature] # average temperature of cladding interior
type = SideAverageValue
boundary = 7
variable = temp
execute_on = 'initial timestep_end'
[]
[average_centerline_fuel_temperature] # average temperature of the cladding interior and all pellet exteriors
type = SideAverageValue
boundary = 9
variable = temp
execute_on = 'initial linear'
[]
[plenum_temperature]
type = SideAverageValue
boundary = 9
variable = temp
execute_on = 'initial timestep_end'
[]
[plenum_volume] # gas volume
type = InternalVolume
boundary = 9
addition = 1.3e-5 #rough guess of plenum volume/unit length of fuel
execute_on = 'initial linear'
[]
[pellet_volume] # fuel pellet total volume
type = InternalVolume
boundary = 8
execute_on = 'initial timestep_end'
[]
[clad_inner_vol] # volume inside of cladding
type = InternalVolume
boundary = 7
outputs = exodus
execute_on = 'initial timestep_end'
[]
[fission_gas_generated] # fission gas produced (moles)
type = ElementIntegralFisGasGeneratedSifgrs
block = pellet_type_1
execute_on = linear
[]
[fission_gas_released] # fission gas released to plenum (moles)
type = ElementIntegralFisGasReleasedSifgrs
block = pellet_type_1
execute_on = linear
[]
[flux_from_clad] # area integrated heat flux from the cladding
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 5
diffusivity = thermal_conductivity
execute_on = timestep_end
[]
[flux_from_fuel] # area integrated heat flux from the fuel
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 10
diffusivity = thermal_conductivity
execute_on = timestep_end
[]
[_dt] # time step
type = TimestepSize
execute_on = timestep_end
[]
[num_lin_it]
type = NumLinearIterations
[]
[num_nonlin_it]
type = NumNonlinearIterations
[]
[tot_lin_it]
type = CumulativeValuePostprocessor
postprocessor = num_lin_it
[]
[tot_nonlin_it]
type = CumulativeValuePostprocessor
postprocessor = num_nonlin_it
[]
[alive_time]
type = PerfGraphData
section_name = Root
data_type = TOTAL
[]
[rod_total_power]
type = ElementIntegralPower
variable = temp
fission_rate = fission_rate
block = pellet_type_1
execute_on = timestep_end
[]
[rod_input_power]
type = FunctionValuePostprocessor
function = power_history
scale_factor = 0.1186 # rod height
execute_on = timestep_end
[]
[fission_gas_released_percentage]
type = FGRPercent
fission_gas_released = fission_gas_released
fission_gas_generated = fission_gas_generated
[]
[]
[VectorPostprocessors]
[contact_pressure]
type = NodalValueSampler
sort_by = x
use_displaced_mesh = true
variable = spacer_clad_mechanical_normal_lm
boundary = 102
[]
[frictional_pressure]
type = NodalValueSampler
sort_by = x
use_displaced_mesh = true
variable = spacer_clad_mechanical_tangential_lm
boundary = 102
[]
[worn_depth]
type = NodalValueSampler
sort_by = x
use_displaced_mesh = true
variable = worn_depth
boundary = 102
execute_on = TIMESTEP_END
[]
[]
[Outputs]
perf_graph = true
exodus = true
csv = true
execute_on = 'FINAL'
[console]
type = Console
max_rows = 25
[]
checkpoint = true
file_base = 'step_${step_number}'
[]
(test/tests/side_ave_incr_tensor_component/side_ave_incr_component.i)
#--------------------------------------------------------------------------------------------------------
# REGRESSION TEST FOR THE POSTPROCESSOR: SideAverageIncrementTensorComponent
# ====================================================
#
# time creep_zz incremental creep incremental creep
# (calculated) (SideAverageIncrementTensorComponent)
# 0 0.00000 0.00000 0.00000
# 100 1.5011e-07 1.5011e-07 1.5011e-07
# 200 3.0022e-07 1.5011e-07 1.5011e-07
# 300 4.5033e-07 1.5011e-07 1.5011e-07
# 400 6.0044e-07 1.5011e-07 1.5011e-07
# 500 7.5055e-07 1.5011e-07 1.5011e-07
# 600 9.0066e-07 1.5011e-07 1.5011e-07
# 700 1.0508e-06 1.5011e-07 1.5011e-07
# 800 1.2009e-06 1.5011e-07 1.5011e-07
# 900 1.3501e-06 1.5011e-07 1.5011e-07
# 1000 1.5011e-06 1.5011e-07 1.5011e-07
#--------------------------------------------------------------------------------------------------------
[GlobalParams]
displacements = 'disp_x disp_y'
temperature = temp
[]
[Mesh]
coord_type = RZ
[mesh]
type = FileMeshGenerator
file = clad_rz.e
[]
[]
[Variables]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 630.0
[]
[]
[AuxVariables]
[fast_neutron_flux]
order = FIRST
family = LAGRANGE
[]
[creep_strain_zz]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[pressure_ramp]
type = PiecewiseConstant
data_file = int_pressure.csv
format = columns
scale_factor = 1.0
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[clad]
block = 1
add_variables = true
strain = FINITE
eigenstrain_names = 'thermal_eigenstrain'
generate_output = 'creep_strain_zz'
decomposition_method = EigenSolution
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
factor = 3e17
[]
[]
[BCs]
[Pressure]
[internal_pressure]
boundary = 1
factor = 1.0
function = pressure_ramp
[]
[]
[u_bottom_fix]
type = DirichletBC
variable = disp_y
boundary = 4
value = 0.0
[]
[temp_bc]
type = DirichletBC
variable = temp
boundary = '1 2 3 4'
value = 630.0
[]
[]
[Materials]
[clad_thermal]
type = HeatConductionMaterial
block = 1
thermal_conductivity = 100.0
specific_heat = 1.0
[]
[clad_elasticity_tensor]
type = ComputeIsotropicElasticityTensor
block = 1
youngs_modulus = 8.0e10
poissons_ratio = 0.3
[]
[clad_creep_model]
type = ZryCreepHayesHoppeUpdate
block = 1
fast_neutron_flux = fast_neutron_flux
zircaloy_material_type = stress_relief_annealed
model_irradiation_creep = true
model_thermal_creep = false
[]
[clad_inelastic_stress]
type = ComputeMultipleInelasticStress
block = 1
tangent_operator = elastic
inelastic_models = 'clad_creep_model'
[]
[clad_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = 1
thermal_expansion_coeff = 5.0e-6
stress_free_temperature = 573.15
eigenstrain_name = 'thermal_eigenstrain'
[]
[density]
type = StrainAdjustedDensity
block = 1
strain_free_density = 6500
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-snes_ksp_ew'
# petsc_options = '-snes_mf -ksp_monitor -snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart'
petsc_options_value = '101'
line_search = 'none'
l_max_its = 100
nl_max_its = 100
nl_rel_tol = 1e-8
nl_abs_tol = 1e-6
l_tol = 1e-5
start_time = 0.0
end_time = 1000
dt = 100
[]
[Postprocessors]
[crp_zz]
type = SideAverageValue
boundary = 1
variable = creep_strain_zz
[]
[incr_cr]
type = SideAverageIncrementTensorComponent
boundary = 1
variable = creep_strain_zz
[]
[]
[Outputs]
file_base = out
exodus = true
[]
(test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep.i)
#
# The mesh is a 1x1x1 cube with a constant pressure of 10 MPa on the top face.
# Symmetry boundary conditions on three planes provide a uniaxial stress
# field. The temperature is held constant at 1000; only creep is modeled: no
# platicity material model is included in the material blocks. The solution is
# advanced through ten time steps of 1e4 for a total time of 1e5.
#
# The total strain at time 1e5 can be computed as:
#
# e_tot = e_elas +
# e_thermal_creep +
# e_irradiation_creep
#
# = P/E +
# A * (sigma/G)**n * exp(-Q/(RT)) * dt +
# C0 * (phi)**C1 * (sigma/10**6)**C2 * dt
#
# where P = pressure load
# E = Young's modulus
# A = thermal creep coefficient
# G = shear modulus = 0.5*E/(1+nu)
# nu = Poisson's ratio
# sigma = stress
# n = thermal creep exponent
# Q = activation energy
# R = gas constant
# T = temperature
# C0 = irradiation creep coefficient
# phi = fast neutron flux
# C1 = neutron flux exponent
# C2 = stress exponent
# dt = total time
#
# For this test, the analytical solutuon at t = 1e5 is:
#
# G = 0.5*2e11/(1+0.3) = 7.692308e10
#
# e_tot = 10e6/2e11 +
# 3.14e24 * (10e6/G)**5 * exp(-2.7e5/(8.3143*1000) * 1e5
# 9.881e-28 * (1e19)**0.85 * (10e6/1e6)**1 * 1e5
#
# = 5e-5 +
# 9.189527e-5
# 1.395728e-5
#
# = 1.558526e-4
#
# For both small strain (ComputeIncrementalSmallStrain) or finite strain (ComputeFiniteStrain)
# kinematics, BISON computes:
#
# e_elas = 5e-5
# e_creep = 1.05819e-4
# e_tot = 1.55819e-4
#
#
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 3
[]
[]
[GlobalParams]
displacements = 'x_disp y_disp z_disp'
[]
[Variables]
[x_disp]
order = FIRST
family = LAGRANGE
[]
[y_disp]
order = FIRST
family = LAGRANGE
[]
[z_disp]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1000.0
[]
[]
[AuxVariables]
[fast_neutron_flux]
order = FIRST
family = LAGRANGE
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[elastic_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[creep_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[top_pull]
type = PiecewiseLinear
x = '0 1'
y = '1 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = finite
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
factor = 1e19
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[elastic_strain_yy]
type = RankTwoAux
rank_two_tensor = elastic_strain
variable = elastic_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[creep_strain_yy]
type = RankTwoAux
rank_two_tensor = creep_strain
variable = creep_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[]
[BCs]
[u_top_pull]
type = Pressure
variable = y_disp
boundary = top
factor = -10.0e6
function = top_pull
[]
[u_bottom_fix]
type = DirichletBC
variable = y_disp
boundary = bottom
value = 0.0
[]
[u_yz_fix]
type = DirichletBC
variable = x_disp
boundary = left
value = 0.0
[]
[u_xy_fix]
type = DirichletBC
variable = z_disp
boundary = back
value = 0.0
[]
[temp_fix]
type = DirichletBC
variable = temp
boundary = 'bottom top'
value = 1000.0
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 2.0e11
poissons_ratio = 0.3
[]
[stress]
type = ComputeMultipleInelasticStress
tangent_operator = elastic
inelastic_models = 'zrycreep'
[]
[zrycreep]
type = ZryCreepHayesHoppeUpdate
temperature = temp
fast_neutron_flux = fast_neutron_flux
model_irradiation_creep = true
model_thermal_creep = true
[]
[thermal]
type = HeatConductionMaterial
block = 0
specific_heat = 1.0
thermal_conductivity = 100.
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 1.0
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart'
petsc_options_value = '101'
line_search = 'none'
l_max_its = 100
nl_max_its = 100
l_tol = 1e-5
start_time = 0.0
num_steps = 10
dt = 10000
[]
[Outputs]
show = 'x_disp y_disp z_disp temp stress_yy elastic_strain_yy creep_strain_yy'
[out]
type = Exodus
[]
[]
(test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep_highflux.i)
# Tests material model for Hayes-Kassner creep in the clad which computes the combined thermal
# and irradiation creep for Zr4 alloy.
#
# The mesh is a 1x1x1 cube with a constant pressure of 10 MPa on the top face.
# Symmetry boundary conditions on three planes provide a uniaxial stress
# field. The temperature is held constant at 1000; only creep is modeled: no
# platicity material model is included in the material blocks. The solution is
# advanced through ten time steps of 1e4 for a total time of 1e5.
#
# The total strain at time 1e5 can be computed as:
#
# e_tot = e_elas +
# e_thermal_creep +
# e_irradiation_creep
#
# = P/E +
# A * (sigma/G)**n * exp(-Q/(RT)) * dt +
# C0 * (phi)**C1 * (sigma/10**6)**C2 * dt
#
# where P = pressure load
# E = Young's modulus
# A = thermal creep coefficient
# G = shear modulus = 0.5*E/(1+nu)
# nu = Poisson's ratio
# sigma = stress
# n = thermal creep exponent
# Q = activation energy
# R = gas constant
# T = temperature
# C0 = irradiation creep coefficient
# phi = fast neutron flux
# C1 = neutron flux exponent
# C2 = stress exponent
# dt = total time
#
# For this test, the analytical solutuon at t = 1e5 is:
#
# G = 0.5*2e11/(1+0.3) = 7.692308e10
#
# e_tot = 10e6/2e11 +
# 3.14e24 * (10e6/G)**5 * exp(-2.7e5/(8.3143*1000) * 1e5 +
# 9.881e-28 * (5e19)**0.85 * (10e6/1e6)**1 * 1e5
#
# = 5e-5 +
# 9.189527e-5
# 5.481828e-5
#
# = 1.9467135e-4
#
# For both small strain (ComputeIncrementalSmallStrain) or finite strain (ComputeFiniteStrain)
# kinematics, BISON computes:
# e_elas = 5e-5
# e_creep = 1.4671e-4
# e_tot = 1.9671e-4
#
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 3
[]
[]
[GlobalParams]
displacements = 'x_disp y_disp z_disp'
[]
[Variables]
[x_disp]
order = FIRST
family = LAGRANGE
[]
[y_disp]
order = FIRST
family = LAGRANGE
[]
[z_disp]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1000.0
[]
[]
[AuxVariables]
[fast_neutron_flux]
order = FIRST
family = LAGRANGE
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[elastic_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[creep_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[top_pull]
type = PiecewiseLinear
x = '0 1'
y = '1 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = finite
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
factor = 5.0e19
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[elastic_strain_yy]
type = RankTwoAux
rank_two_tensor = elastic_strain
variable = elastic_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[creep_strain_yy]
type = RankTwoAux
rank_two_tensor = creep_strain
variable = creep_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[]
[BCs]
[u_top_pull]
type = Pressure
variable = y_disp
boundary = top
factor = -10.0e6
function = top_pull
[]
[u_bottom_fix]
type = DirichletBC
variable = y_disp
boundary = bottom
value = 0.0
[]
[u_yz_fix]
type = DirichletBC
variable = x_disp
boundary = left
value = 0.0
[]
[u_xy_fix]
type = DirichletBC
variable = z_disp
boundary = back
value = 0.0
[]
[temp_fix]
type = DirichletBC
variable = temp
boundary = 'top bottom'
value = 1000.0
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 2.0e11
poissons_ratio = 0.3
[]
[stress]
type = ComputeMultipleInelasticStress
tangent_operator = elastic
inelastic_models = 'zrycreep'
[]
[zrycreep]
type = ZryCreepHayesHoppeUpdate
temperature = temp
fast_neutron_flux = fast_neutron_flux
model_irradiation_creep = true
model_thermal_creep = true
[]
[thermal]
type = HeatConductionMaterial
block = 0
specific_heat = 1.0
thermal_conductivity = 100.
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 1.0
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart'
petsc_options_value = '101'
line_search = 'none'
l_max_its = 100
nl_max_its = 100
l_tol = 1e-5
start_time = 0.0
num_steps = 10
dt = 10000
[]
[Outputs]
show = 'x_disp y_disp z_disp temp stress_yy'
[out]
type = Exodus
[]
[]
(test/tests/solid_mechanics/zry_creep/hayes_hoppe/thermal_irradiation_creep_rz_highflux.i)
# Tests material model for Hayes-Kassner creep in the clad which computes the combined thermal
# and irradiation creep for Zr4 alloy.
#
# The mesh is a 1x1x1 cube with a constant pressure of 10 MPa on the top face.
# Symmetry boundary conditions on three planes provide a uniaxial stress
# field. The temperature is held constant at 1000; only creep is modeled: no
# platicity material model is included in the material blocks. The solution is
# advanced through ten time steps of 1e4 for a total time of 1e5.
#
# The total strain at time 1e5 can be computed as:
#
# e_tot = e_elas +
# e_thermal_creep +
# e_irradiation_creep
#
# = P/E +
# A * (sigma/G)**n * exp(-Q/(RT)) * dt +
# C0 * (phi)**C1 * (sigma/10**6)**C2 * dt
#
# where P = pressure load
# E = Young's modulus
# A = thermal creep coefficient
# G = shear modulus = 0.5*E/(1+nu)
# nu = Poisson's ratio
# sigma = stress
# n = thermal creep exponent
# Q = activation energy
# R = gas constant
# T = temperature
# C0 = irradiation creep coefficient
# phi = fast neutron flux
# C1 = neutron flux exponent
# C2 = stress exponent
# dt = total time
#
# For this test, the analytical solutuon at t = 1e5 is:
#
# G = 0.5*2e11/(1+0.3) = 7.692308e10
#
# e_tot = 10e6/2e11 +
# 3.14e24 * (10e6/G)**5 * exp(-2.7e5/(8.3143*1000) * 1e5 +
# 9.881e-28 * (5e19)**0.85 * (10e6/1e6)**1 * 1e5
#
# = 5e-5 +
# 9.189527e-5
# 5.481828e-5
#
# = 1.9467135e-4
#
# For both small strain (ComputeIncrementalSmallStrain) or finite strain (ComputeFiniteStrain)
# kinematics, BISON computes:
# e_elas = 5e-5
# e_creep = 1.4671e-4
# e_tot = 1.9671e-4
#
[Mesh]
coord_type = RZ
use_displaced_mesh = false
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Variables]
[disp_x]
order = FIRST
family = LAGRANGE
[]
[disp_y]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1000.0
[]
[]
[AuxVariables]
[fast_neutron_flux]
order = FIRST
family = LAGRANGE
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[elastic_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[creep_strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[top_pull]
type = PiecewiseLinear
x = '0 1'
y = '1 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = small
incremental = true
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
factor = 5e19
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[elastic_strain_yy]
type = RankTwoAux
rank_two_tensor = elastic_strain
variable = elastic_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[creep_strain_yy]
type = RankTwoAux
rank_two_tensor = creep_strain
variable = creep_strain_yy
index_j = 1
index_i = 1
execute_on = timestep_end
[]
[]
[BCs]
[Pressure]
[u_top_pull]
boundary = top
factor = -10.0e6
function = top_pull
[]
[]
[u_bottom_fix]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[]
[temp_fix]
type = DirichletBC
variable = temp
boundary = 'top bottom left right'
value = 1000.0
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 2.0e11
poissons_ratio = 0.3
[]
[stress]
type = ComputeMultipleInelasticStress
tangent_operator = elastic
inelastic_models = 'zrycreep'
[]
[zrycreep]
type = ZryCreepHayesHoppeUpdate
temperature = temp
fast_neutron_flux = fast_neutron_flux
model_irradiation_creep = true
model_thermal_creep = true
[]
[thermal]
type = HeatConductionMaterial
block = 0
specific_heat = 1.0
thermal_conductivity = 100.
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 1.0
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart'
petsc_options_value = '101'
line_search = 'none'
l_max_its = 100
nl_max_its = 100
nl_rel_tol = 1e-8
nl_abs_tol = 1e-6
l_tol = 1e-5
start_time = 0.0
num_steps = 10
dt = 10000
[]
[Outputs]
[out]
type = Exodus
[]
[]
(test/tests/fuelrodlinevaluesampler/example_problem_smeared_test.i)
initial_fuel_density = 10431.0
[GlobalParams]
density = ${initial_fuel_density}
displacements = 'disp_x disp_y'
energy_per_fission = 3.2e-11 # J/fission
temperature = temp
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
[]
[Mesh]
coord_type = RZ
displacements = 'disp_x disp_y'
patch_update_strategy = auto
patch_size = 10 # For contact algorithm
partitioner = centroid
centroid_partitioner_direction = y
[mesh]
type = FileMeshGenerator
file = SmearedTwoPelletOneType2D.e
[]
[]
[Variables]
[temp]
initial_condition = 580.0 # set initial temp to coolant inlet
[]
[]
[AuxVariables]
[fast_neutron_flux]
block = clad
[]
[fast_neutron_fluence]
block = clad
[]
[grain_radius]
block = pellet_type_1
initial_condition = 10e-6
[]
[gap_cond]
order = CONSTANT
family = MONOMIAL
[]
[coolant_htc]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear
data_file = powerhistory.csv
scale_factor = 1
[]
[axial_peaking_factors]
type = ParsedFunction
expression = 1
[]
[pressure_ramp]
type = PiecewiseLinear
x = '-200 0'
y = '0 1'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[fuel]
block = pellet_type_1
strain = FINITE
incremental = true
extra_vector_tags = 'ref'
add_variables = true
decomposition_method = EigenSolution
eigenstrain_names = 'fuel_volumetric_swelling_eigenstrain
fuel_relocation_eigenstrain fuel_thermal_eigenstrain'
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress'
[]
[clad]
block = clad
strain = FINITE
incremental = true
extra_vector_tags = 'ref'
add_variables = true
decomposition_method = EigenSolution
eigenstrain_names = 'clad_thermal_strain clad_irradiation_growth_eigenstrain'
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress'
[]
[]
[Kernels]
[gravity]
type = Gravity
variable = disp_y
value = -9.81
[]
[heat]
type = HeatConduction
variable = temp
extra_vector_tags = 'ref'
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
extra_vector_tags = 'ref'
[]
[heat_source]
type = NeutronHeatSource
variable = temp
extra_vector_tags = 'ref'
block = pellet_type_1
burnup_function = burnup
[]
[]
[UserObjects]
[pin_geometry]
type = FuelPinGeometry
[]
[]
[Burnup]
[burnup]
block = pellet_type_1
rod_ave_lin_pow = power_history # using the power function defined above
axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
num_radial = 80
num_axial = 11
fuel_pin_geometry = 'pin_geometry'
fuel_volume_ratio = 0.987775 # for use with dished pellets (ratio of actual volume to cylinder volume)
order = CONSTANT
family = MONOMIAL
RPF = RPF
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
block = clad
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
factor = 3e13
execute_on = timestep_begin
[]
[fast_neutron_fluence]
type = FastNeutronFluenceAux
variable = fast_neutron_fluence
block = clad
fast_neutron_flux = fast_neutron_flux
execute_on = timestep_begin
[]
[grain_radius]
type = GrainRadiusAux
block = pellet_type_1
variable = grain_radius
temperature = temp
execute_on = linear
[]
[conductance]
type = MaterialRealAux
property = gap_conductance
variable = gap_cond
boundary = 10
execute_on = 'initial timestep_end'
[]
[coolant_htc]
type = MaterialRealAux
property = coolant_channel_htc
variable = coolant_htc
boundary = 2
execute_on = 'initial timestep_end'
[]
[]
[Contact]
[pellet_clad_mechanical]
primary = 5
secondary = 10
formulation = KINEMATIC
model = frictionless
normalize_penalty = true
penalty = 1e14
[]
[]
[ThermalContact]
[thermal_contact]
type = GasGapHeatTransfer
variable = temp
primary = 5
secondary = 10
initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
gas_released = fis_gas_released # coupling to a postprocessor which supplies the fission gas addition
contact_pressure = contact_pressure
[]
[]
[BCs]
[no_x_all] # pin pellets and clad along axis of symmetry (y)
type = DirichletBC
variable = disp_x
boundary = 12
value = 0.0
[]
[no_y_clad_bottom] # pin clad bottom in the axial direction (y)
type = DirichletBC
variable = disp_y
boundary = '1'
value = 0.0
[]
[no_y_fuel_bottom] # pin fuel bottom in the axial direction (y)
type = DirichletBC
variable = disp_y
boundary = '1020'
value = 0.0
[]
[Pressure] # apply coolant pressure on clad outer walls
[coolantPressure]
boundary = '1 2 3'
factor = 15.5e6
function = pressure_ramp # use the pressure_ramp function defined above
[]
[]
[PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
[plenumPressure]
boundary = 9
initial_pressure = 2.0e6
startup_time = -200
R = 8.3143
output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
temperature = ave_temp_interior # coupling to post processor to get gas temperature approximation
volume = gas_volume # coupling to post processor to get gas volume
material_input = fis_gas_released # coupling to post processor to get fission gas added
output = plenum_pressure # coupling to post processor to output plenum/gap pressure
displacements = 'disp_x disp_y'
execute_on = 'initial linear'
[]
[]
[]
[CoolantChannel]
[convective_clad_surface] # apply convective boundary to clad outer surface
boundary = '1 2 3'
variable = temp
inlet_temperature = 580 # K
inlet_pressure = 15.5e6 # Pa
inlet_massflux = 3800 # kg/m^2-sec
rod_diameter = 0.948e-2 # m
rod_pitch = 1.26e-2 # m
linear_heat_rate = power_history
axial_power_profile = axial_peaking_factors
[]
[]
[Materials]
[fuel_thermal]
type = UO2Thermal
block = pellet_type_1
thermal_conductivity_model = NFIR
initial_porosity = 0.0
temperature = temp
burnup_function = burnup
[]
[fuel_volumetric_swelling]
type = UO2VolumetricSwellingEigenstrain
block = pellet_type_1
burnup_function = burnup
initial_fuel_density = 10431.0
eigenstrain_name = 'fuel_volumetric_swelling_eigenstrain'
[]
[fuel_elasticity_tensor]
type = UO2ElasticityTensor
block = pellet_type_1
[]
[fuel_thermal_expansion]
type = UO2ThermalExpansionMartinEigenstrain
block = pellet_type_1
stress_free_temperature = 295
eigenstrain_name = 'fuel_thermal_eigenstrain'
[]
[hotpressing]
type = UO2HotPressingCreepUpdate
block = pellet_type_1
burnup_function = burnup
initial_grain_radius = 10.0e-6
[]
[radial_return_stress]
type = ComputeMultipleInelasticStress
tangent_operator = elastic
inelastic_models = 'hotpressing'
block = pellet_type_1
[]
[fuel_relocation]
type = UO2RelocationEigenstrain
block = pellet_type_1
burnup_function = burnup
fuel_pin_geometry = 'pin_geometry'
relocation_activation1 = 5000 #TM default value
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
burnup_relocation_stop = 1.e20
eigenstrain_name = 'fuel_relocation_eigenstrain'
[]
[clad_thermal]
type = HeatConductionMaterial
block = clad
thermal_conductivity = 16.0
specific_heat = 330.0
[]
[clad_elasticity_tensor]
type = ZryElasticityTensor
block = clad
[]
[clad_creep_model]
type = ZryCreepHayesHoppeUpdate
block = clad
fast_neutron_flux = fast_neutron_flux
model_irradiation_creep = true
model_thermal_creep = true
[]
[clad_inelastic_stress]
type = ComputeMultipleInelasticStress
block = clad
tangent_operator = elastic
inelastic_models = 'clad_creep_model'
[]
[clad_thermal_expansion]
type = ComputeThermalExpansionEigenstrain
block = clad
thermal_expansion_coeff = 5.0e-6
stress_free_temperature = 295.0
eigenstrain_name = clad_thermal_strain
[]
[clad_irradiation_growth]
type = ZryIrradiationGrowthEigenstrain
block = clad
fast_neutron_fluence = fast_neutron_fluence
zircaloy_material_type = ESCORE_IrradiationGrowthZr4
eigenstrain_name = clad_irradiation_growth_eigenstrain
[]
[fission_gas_release]
type = UO2Sifgrs
block = pellet_type_1
temperature = temp
burnup_function = burnup
grain_radius = grain_radius
gbs_model = true
[]
[clad_density]
type = StrainAdjustedDensity
block = clad
strain_free_density = 6551.0
[]
[fuel_density]
type = StrainAdjustedDensity
block = pellet_type_1
strain_free_density = ${initial_fuel_density}
[]
[]
[Dampers]
[BoundingValueNodalDamper]
type = BoundingValueNodalDamper
variable = temp
max_value = 3200
min_value = 300
[]
[]
[Preconditioning]
[SMP]
type = SMP
coupled_groups = 'disp_x,disp_y'
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-pc_type_asm'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
petsc_options_value = 'lu superlu_dist 51'
line_search = 'none'
verbose = true
l_max_its = 100
l_tol = 1e-5 #8e-3
nl_max_its = 15
nl_rel_tol = 1e-10
nl_abs_tol = 1e-8
start_time = -200
num_steps = 2
dtmax = 2e6
dtmin = 1
[TimeStepper]
type = IterationAdaptiveDT
dt = 2.0e2
optimal_iterations = 6
iteration_window = 2
[]
[]
[Postprocessors]
[ave_temp_interior] # average temperature of the cladding interior and all pellet exteriors
type = SideAverageValue
boundary = 9
variable = temp
execute_on = 'initial linear'
[]
[clad_inner_vol] # volume inside of cladding
type = InternalVolume
boundary = 7
outputs = exodus
execute_on = 'initial timestep_end'
[]
[pellet_volume] # fuel pellet total volume
type = InternalVolume
boundary = 8
outputs = exodus
execute_on = 'initial timestep_end'
[]
[avg_clad_temp] # average temperature of cladding interior
type = SideAverageValue
boundary = 7
variable = temp
execute_on = 'initial timestep_end'
[]
[fis_gas_produced] # fission gas produced (moles)
type = ElementIntegralFisGasGeneratedSifgrs
block = pellet_type_1
execute_on = timestep_end
[]
[fis_gas_released] # fission gas released to plenum (moles)
type = ElementIntegralFisGasReleasedSifgrs
block = pellet_type_1
execute_on = timestep_end
[]
[fis_gas_grain]
type = ElementIntegralFisGasGrainSifgrs
block = pellet_type_1
outputs = exodus
execute_on = timestep_end
[]
[fis_gas_boundary]
type = ElementIntegralFisGasBoundarySifgrs
block = pellet_type_1
outputs = exodus
execute_on = timestep_end
[]
[gas_volume] # gas volume
type = InternalVolume
boundary = 9
component = 1
execute_on = 'initial linear'
[]
[flux_from_clad] # area integrated heat flux from the cladding
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 5
diffusivity = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[flux_from_fuel] # area integrated heat flux from the fuel
type = SideDiffusiveFluxIntegral
variable = temp
boundary = 10
diffusivity = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[_dt] # time step
type = TimestepSize
execute_on = timestep_end
[]
[nonlinear_its] # number of nonlinear iterations at each timestep
type = NumNonlinearIterations
execute_on = timestep_end
[]
[rod_total_power]
type = ElementIntegralPower
variable = temp
burnup_function = burnup
block = pellet_type_1
execute_on = 'initial timestep_end'
[]
[rod_input_power]
type = FunctionValuePostprocessor
function = power_history
scale_factor = 0.02372 # rod height
execute_on = 'initial timestep_end'
[]
[]
[VectorPostprocessors]
[fuel_vonmises]
type = FuelRodLineValueSampler
variable = vonmises_stress
material = 'fuel'
fraction = 0.51
num_points = 20
orientation = 'horizontal'
fuel_pin_geometry = 'pin_geometry'
outputs = chkfile
[]
[clad_vonmises]
type = FuelRodLineValueSampler
variable = vonmises_stress
material = 'clad'
fraction = 0.51
num_points = 9
orientation = 'horizontal'
fuel_pin_geometry = 'pin_geometry'
outputs = chkfile
[]
[]
[Outputs]
exodus = true
color = false
csv = true
[console]
type = Console
output_linear = true
max_rows = 25
[]
[chkfile]
type = CSV
execute_on = 'FINAL'
[]
[]
(include/materials/solid_mechanics/ZryCreepTulkkiHayesHoppeUpdate.h)
/*************************************************/
/* DO NOT MODIFY THIS HEADER */
/* */
/* BISON */
/* */
/* (c) 2015 Battelle Energy Alliance, LLC */
/* ALL RIGHTS RESERVED */
/* */
/* Prepared by Battelle Energy Alliance, LLC */
/* Under Contract No. DE-AC07-05ID14517 */
/* With the U. S. Department of Energy */
/* */
/* See COPYRIGHT for full restrictions */
/*************************************************/
#pragma once
#include "ZryCreepHayesHoppeUpdate.h"
class ZryCreepTulkkiHayesHoppeUpdate : public ZryCreepHayesHoppeUpdate
{
public:
static InputParameters validParams();
ZryCreepTulkkiHayesHoppeUpdate(const InputParameters & parameters);
protected:
/*
* Calculates the Primary Viscoelastic Creep from Tulkki (2014) model.
* The secondary Thermal Creep is computed using the Hayes and Kasner (2006)
* model, while the secondary Irradiation Creep is computed using the Hoppe
* model (1991)
*/
virtual void initQpStatefulProperties() override;
virtual void computePrimaryViscoElasticCreep();
virtual void computeStressFinalize(const RankTwoTensor & inelasticStrainIncrement) override;
/// Parameters for Visco Elastic Creep for sudden load drops and load reversals
MaterialProperty<Real> & _primary_viscoelastic_creep;
const MaterialProperty<Real> & _primary_viscoelastic_creep_old;
/// Stress Like variable Tulkki(2014)
MaterialProperty<Real> & _stress_like_variable;
const MaterialProperty<Real> & _stress_like_variable_old;
/// Creep Strain in the viscoelastic setting (primary creep is recoverable)
MaterialProperty<Real> & _creep_strain_viscoelastic;
const MaterialProperty<RankTwoTensor> & _stress;
const MaterialProperty<RankTwoTensor> & _stress_old;
const MaterialProperty<RankTwoTensor> & _creep_strain;
};