- eigenstrain_nameMaterial property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator.
C++ Type:std::string
Controllable:No
Description:Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator.
- fast_neutron_fluencefast neutron fluence
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:fast neutron fluence
- temperatureCoupled Temperature in Kelvin
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled Temperature in Kelvin
CompositeSiCVolumetricSwellingEigenstrain
Computes volumetric swelling of SiC as a function of temperature and fluence.
Description
Swelling due to irradiation of composite silicon carbide is modeled in CompositeSiCVolumetricSwellingEigenstrain. The same model is also used for monolithic silicon carbide. Two correlations are available: KATOH and MIELOSZYK.
For the KATOH option, the differential equation for the swelling strain is given by Stone et al. (2015): where is the fast neutron fluence in dpa (displacements per atom), is a temperature dependent rate constant, and is a temperature dependent characteristic fast neutron fluence in dpa. A conversion factor of 1 dpa = 1n/m is used to convert the fast neutron fluence (Katoh et al., 2014). The temperature dependent rate constant is given by: where is the temperature in K. The temperature dependent characteristic fast neutron fluence is given by: where is the temperature in K.
warningwarning:Timestep Size Dependence
The KATOH swelling model cannot be analytically integrated and is particularly nonlinear at fluences less than 1.0 dpa. An incremental approach is used to calculate the swelling strain and internal sub cycling is applied at fluences less than 1.0 dpa to help mitigate the error due to time step size. The number of sub cycles is controlled by the _number_of_substeps parameter.
For the MIELOSZYK option, the increment of volumetric swelling strain is given by Mieloszyk (2015): where is the fluence at the current time, is the fluence at the previous time, and , which represents the fluence in DPA at which the swelling is saturated. Also, where is the temperature (K). and are algorithmic parameters. See Mieloszyk (2015) for details.
An incremental approach is used to solve for the total volumetric strain.
Example Input Syntax
[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain<<<{"description": "Computes volumetric swelling of SiC as a function of temperature and fluence.", "href": "CompositeSiCVolumetricSwellingEigenstrain.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 0
temperature<<<{"description": "Coupled Temperature in Kelvin"}>>> = temp
swelling_model<<<{"description": "Options for the correlation used to calculate thermal conductivity"}>>> = KATOH
number_of_substeps<<<{"description": "The number of substeps to take when solving for the swelling increment in the Katoh formulation at low fluences."}>>> = 100
fast_neutron_fluence<<<{"description": "fast neutron fluence"}>>> = fast_neutron_fluence
eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = swelling_strain
[]
[][Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain<<<{"description": "Computes volumetric swelling of SiC as a function of temperature and fluence.", "href": "CompositeSiCVolumetricSwellingEigenstrain.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 0
temperature<<<{"description": "Coupled Temperature in Kelvin"}>>> = temp
swelling_model<<<{"description": "Options for the correlation used to calculate thermal conductivity"}>>> = MIELOSZYK
fast_neutron_fluence<<<{"description": "fast neutron fluence"}>>> = fast_neutron_fluence
eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = swelling_strain
[]
[]The eigenstrain name must also be passed to the strain calculator, and an example parameter setting in the Solid Mechanics QuasiStatic Action is shown below:
[Physics<<<{"href": "../../../syntax/Physics/index.html"}>>>]
[SolidMechanics<<<{"href": "../../../syntax/Physics/SolidMechanics/index.html"}>>>]
[QuasiStatic<<<{"href": "../../../syntax/Physics/SolidMechanics/QuasiStatic/index.html"}>>>]
[square]
block<<<{"description": "The list of ids of the blocks (subdomain) that the stress divergence kernels will be applied to"}>>> = 0
add_variables<<<{"description": "Add the displacement variables"}>>> = false
strain<<<{"description": "Strain formulation"}>>> = FINITE
eigenstrain_names<<<{"description": "List of eigenstrains to be applied in this strain calculation"}>>> = 'swelling_strain'
[]
[]
[]
[]Input Parameters
- base_nameOptional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
C++ Type:std::string
Controllable:No
Description:Optional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
- 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
- computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
- 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
Options:NONE, ELEMENT, SUBDOMAIN
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.
- number_of_substeps100The number of substeps to take when solving for the swelling increment in the Katoh formulation at low fluences.
Default:100
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The number of substeps to take when solving for the swelling increment in the Katoh formulation at low fluences.
- swelling_modelMIELOSZYKOptions for the correlation used to calculate thermal conductivity
Default:MIELOSZYK
C++ Type:MooseEnum
Options:KATOH, MIELOSZYK
Controllable:No
Description:Options for the correlation used to calculate thermal conductivity
Optional Parameters
- 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.
- 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
- 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
- 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
- (test/tests/solid_mechanics/sic_creep/composite_sic_creep.i)
- (test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_mieloszyk_ad.i)
- (test/tests/monolithicSiCThermal/thermal_stone_unirradiated.i)
- (test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh_ad.i)
- (test/tests/solid_mechanics/sic_creep/ad_composite_sic_creep.i)
- (test/tests/monolithicSiCThermal/thermal_stone_irradiated.i)
- (assessment/LWR/validation/ATF_SiC_HFIR/analysis/base.i)
- (test/tests/thermalCompositeSiC/thermal_stone_irradiated.i)
- (test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_mieloszyk.i)
- (examples/accident_tolerant_fuel/u3si2_sic/u3si2_outer_monolith_1.5D.i)
- (test/tests/thermalCompositeSiC/thermal_koyanagi_stone_combined.i)
- (test/tests/thermalCompositeSiC/thermal_stone_unirradiated.i)
- (test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh.i)
References
- Y. Katoh, K. Ozawa, C. Shih, T. Nozawa, R. J. Shinavski, A. Hasegawa, and L. L. Snead.
Continuous SiC fiber, CVI SiC matrix composites for nuclear applications: properties and irradiation effects.
Journal of Nuclear Materials, 448:448–476, 2014.[BibTeX]
@article{katoh2014, author = "Katoh, Y. and Ozawa, K. and Shih, C. and Nozawa, T. and Shinavski, R. J. and Hasegawa, A. and Snead, L. L.", title = "Continuous {SiC} fiber, {CVI} {SiC} matrix composites for nuclear applications: Properties and irradiation effects", journal = "Journal of Nuclear Materials", volume = "448", pages = "448-476", year = "2014", publisher = "Elsevier" } - Alexander Mieloszyk.
Assessing thermo-mechanical performance of ThO2 and SiC clad light water reactor fuel rods with a modular simulation tool.
PhD thesis, Massachusetts Institute of Technology, September 2015.
http://hdl.handle.net/1721.1/103660.[BibTeX]
@phdthesis{mieloszyk2015, author = "Mieloszyk, Alexander", title = "Assessing thermo-mechanical performance of {ThO2} and {SiC} clad light water reactor fuel rods with a modular simulation tool", school = "Massachusetts Institute of Technology", year = "2015", month = "September", note = "http://hdl.handle.net/1721.1/103660" } - J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, and C. A. Back.
Stress analysis and probabilistic assessment of multi-layer sic-based accident tolerant nuclear fuel cladding.
Journal of Nuclear Materials, 466:682–697, 2015.[BibTeX]
@article{stone2015, author = "Stone, J. G. and Schleicher, R. and Deck, C. P. and Jacobsen, G. M. and H. E. Khalifa and Back, C. A.", title = "Stress analysis and probabilistic assessment of multi-layer SiC-based accident tolerant nuclear fuel cladding", journal = "Journal of Nuclear Materials", volume = "466", pages = "682-697", year = "2015", publisher = "Elsevier" }
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh.i)
# This test tests the Katoh swelling model used for composite and monolithic
# silicon carbide. The model is obtained from:
#
# Y. Katoh, K. Ozawa, C. Shih, T. Nozawa, R. J. Shinavski, A. Hasegawa, and
# L. L. Snead. Continuous SiC fiber, CVI SiC matrix composites for nuclear
# applications: properties and irradiation effects. Journal of Nuclear Materials,
# 448, pp. 448-476, 2014.
#
# The swelling rate is a function of the fast neutron fluence:
#
# dS/dgamma = k(T) * gamma^(-1/3) * exp(-gamma/gamma_sc(T))
#
# where gamma is the fluence in dpa, S is the swelling strain, and k(T) and
# gamma_sc(T) are temperature dependent constants. The above equation cannot
# be analytically integrated and an incremental approach is used to calculate the
# swelling strain. The "analytical" solution is taken as the results of the
# numerical integral calculated using Wolfram Alpha or Mathematica. Depending
# upon the number of substeps used the BISON simulations approach this analytical
# solution.
#
# Here, a single 1x1 m 2D element is subjected to a constant temperature of
# 1073.15 K with a varying fast neutron fluence from 0 to 1e25 n/m^2 (0 to 5 dpa)
# over 1 second. Below BISON results are included for 100, 1000 and 10000 substeps to
# illustrate that the results approach that of the analytical solution. This test
# is for the 100 substep case. Tests for the 1000 and 10000 substeps at high
# temperature as well as low and mid range temperatures of 473.15 K
# and 773.15 K (with 100 substeps) are included using cli_args in the tests file.
#
# Fluence BISON Swelling Strain (-) Analytical
# (dpa) 100 substeps 1000 substeps 10000 substeps Swelling
# 1.0 5.91955e-3 5.83353e-3 5.83899e-3 5.84045e-3
# 2.0 7.20081e-3 7.22617e-3 7.23163e-3 7.23270e-3
# 3.0 7.70534e-3 7.73070e-3 7.73617e-3 7.73713e-3
# 4.0 7.90061e-3 7.92597e-3 7.93143e-3 7.93235e-3
# 5.0 7.97852e-3 8.00388e-3 8.00934e-3 8.01025e-3
#
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
order = FIRST
family = LAGRANGE
[]
[disp_y]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1073.15
[]
[]
[AuxVariables]
[fast_neutron_fluence]
order = FIRST
family = LAGRANGE
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[Functions]
[fast_neutron_fluence_func]
type = PiecewiseLinear
x = '0 1.0'
y = '0 5e25'
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence_func
execute_on = timestep_begin
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[]
[BCs]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[all_temp]
type = DirichletBC
variable = temp
boundary = 'bottom left right top'
value = 1073.15
[]
[]
[Materials]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temp
swelling_model = KATOH
number_of_substeps = 100
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
block = 0
[]
[thermalCompositeSiC]
type = CompositeSiCThermal
temperature = temp
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 3100.0
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-snes_ksp_ew '
petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type'
petsc_options_value = '70 hypre boomeramg'
l_max_its = 60
nl_rel_tol = 1e-5
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
dt = 0.01
end_time = 1.0
[]
[Postprocessors]
[avgFuelTemp]
type = ElementAverageValue
variable = temp
execute_on = 'initial timestep_end'
block = 0
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[average_fluence]
type = ElementAverageValue
variable = fast_neutron_fluence
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_mieloszyk.i)
# The purpose of this test is to verify the implementation of the Mieloszyk
# swelling model for composite and monolithic silicon carbide. The algorithm
# used is from:
#
# A. J. Mieloszyk, "Assessing Thermo-Mechanical Performance of ThO2 and SiC Clad
# Light Water Reactor Fuel Rods with a Modular Simulation Tool," PhD Dissertation,
# Massachusetts Institute of Technology, 2015.
#
# In the above reference an example is provided that applies a temperature and
# fluence profile over 300 days (here we simplify that to 300 seconds with the
# same profiles). The reference also provides an expected profile of the
# volumetric swelling as a function of time given those profiles. The BISON
# results from this test are compared to the results provided in Mieloszyk's
# dissertation in the swelling_mieloszyk_testing.xlsx Excel file contained
# within this directory.
#
# The small discrepancies between the BISON and Mieloszyk curves can be
# attributed to the number of points used in the digitization of
# Mieloszyk's data.
#
# The results of predicted swelling at selected times used in the temperature
# function of the test are tabulated below:
#
# Time (s) Temperature (K) Fluence (dpa) Swelling Strain (-)
# 0 650.0 0 0
# 100 650.0 2.0 0.014267
# 120 350.0 2.0 0.014267
# 200 750.0 4.2 0.011403
# 250 650.0 4.9 0.014604
# 300 550.0 5.6 0.017932
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
order = FIRST
family = LAGRANGE
[]
[disp_y]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 650.0
[]
[]
[AuxVariables]
[fast_neutron_fluence]
order = FIRST
family = LAGRANGE
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[]
[Functions]
[fast_neutron_fluence_func]
type = PiecewiseLinear
x = '0 100.0 120.0 200.0 300.0'
y = '0 2e25 2e25 4.2e25 5.6e25'
[]
[temperature]
type = PiecewiseLinear
x = '0 100 100.01 120 120.01 200.0 200.01 250 250.01 300'
y = '650 650 350 350 675 750 600 650 625 550'
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence_func
execute_on = timestep_begin
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[]
[BCs]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[rightX]
type = FunctionDirichletBC
variable = disp_x
boundary = right
function = '1e-6*t'
[]
[topY]
type = FunctionDirichletBC
variable = disp_y
boundary = top
function = '1e-6*t'
[]
[all_temp]
type = FunctionDirichletBC
variable = temp
boundary = 'bottom right top left'
function = temperature
[]
[]
[Materials]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temp
swelling_model = MIELOSZYK
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
block = 0
[]
[thermalCompositeSiC]
type = CompositeSiCThermal
temperature = temp
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 3100.0
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-snes_ksp_ew '
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = none
l_max_its = 60
nl_rel_tol = 1e-8
nl_abs_tol = 1e-3
l_tol = 1e-5
start_time = 0.0
dt = 1
end_time = 300.0
[]
[Postprocessors]
[avgFuelTemp]
type = ElementAverageValue
variable = temp
execute_on = 'initial timestep_end'
block = 0
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[average_fluence]
type = ElementAverageValue
variable = fast_neutron_fluence
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh.i)
# This test tests the Katoh swelling model used for composite and monolithic
# silicon carbide. The model is obtained from:
#
# Y. Katoh, K. Ozawa, C. Shih, T. Nozawa, R. J. Shinavski, A. Hasegawa, and
# L. L. Snead. Continuous SiC fiber, CVI SiC matrix composites for nuclear
# applications: properties and irradiation effects. Journal of Nuclear Materials,
# 448, pp. 448-476, 2014.
#
# The swelling rate is a function of the fast neutron fluence:
#
# dS/dgamma = k(T) * gamma^(-1/3) * exp(-gamma/gamma_sc(T))
#
# where gamma is the fluence in dpa, S is the swelling strain, and k(T) and
# gamma_sc(T) are temperature dependent constants. The above equation cannot
# be analytically integrated and an incremental approach is used to calculate the
# swelling strain. The "analytical" solution is taken as the results of the
# numerical integral calculated using Wolfram Alpha or Mathematica. Depending
# upon the number of substeps used the BISON simulations approach this analytical
# solution.
#
# Here, a single 1x1 m 2D element is subjected to a constant temperature of
# 1073.15 K with a varying fast neutron fluence from 0 to 1e25 n/m^2 (0 to 5 dpa)
# over 1 second. Below BISON results are included for 100, 1000 and 10000 substeps to
# illustrate that the results approach that of the analytical solution. This test
# is for the 100 substep case. Tests for the 1000 and 10000 substeps at high
# temperature as well as low and mid range temperatures of 473.15 K
# and 773.15 K (with 100 substeps) are included using cli_args in the tests file.
#
# Fluence BISON Swelling Strain (-) Analytical
# (dpa) 100 substeps 1000 substeps 10000 substeps Swelling
# 1.0 5.91955e-3 5.83353e-3 5.83899e-3 5.84045e-3
# 2.0 7.20081e-3 7.22617e-3 7.23163e-3 7.23270e-3
# 3.0 7.70534e-3 7.73070e-3 7.73617e-3 7.73713e-3
# 4.0 7.90061e-3 7.92597e-3 7.93143e-3 7.93235e-3
# 5.0 7.97852e-3 8.00388e-3 8.00934e-3 8.01025e-3
#
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
order = FIRST
family = LAGRANGE
[]
[disp_y]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1073.15
[]
[]
[AuxVariables]
[fast_neutron_fluence]
order = FIRST
family = LAGRANGE
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[Functions]
[fast_neutron_fluence_func]
type = PiecewiseLinear
x = '0 1.0'
y = '0 5e25'
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence_func
execute_on = timestep_begin
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[]
[BCs]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[all_temp]
type = DirichletBC
variable = temp
boundary = 'bottom left right top'
value = 1073.15
[]
[]
[Materials]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temp
swelling_model = KATOH
number_of_substeps = 100
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
block = 0
[]
[thermalCompositeSiC]
type = CompositeSiCThermal
temperature = temp
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 3100.0
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-snes_ksp_ew '
petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type'
petsc_options_value = '70 hypre boomeramg'
l_max_its = 60
nl_rel_tol = 1e-5
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
dt = 0.01
end_time = 1.0
[]
[Postprocessors]
[avgFuelTemp]
type = ElementAverageValue
variable = temp
execute_on = 'initial timestep_end'
block = 0
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[average_fluence]
type = ElementAverageValue
variable = fast_neutron_fluence
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(test/tests/solid_mechanics/sic_creep/composite_sic_creep.i)
# The mesh is a 1x1 square (single element) assuming RZ axisymmetry.
# The temperature is held constant at 1073.15 K with a varying fast fluence
# from 0 to 1e25 n/m^2 over 1e6 second. A timestep of 1e4 s is taken and
# swelling is calculated by the Katoh volumetric swelling model with 100 substeps.
# An axial pressure of -3MPa is applied (i.e., pulling outward from the surface).
#
# The creep rate is calculated from:
#
# e_dot = e_dot_pri+e_dot_sec
# e_dot_pri = K_pri * sigma * e_dot_swell
# e_dot_sec = K_sec * sigma * phi
#
# where K_pri and K_sec are conversion factors, sigma is the Von Mises stress in TPa, phi is the fast neutron flux, and
# e_dot_swell is the swelling strain rate. The conversion factors are given by
#
# K_pri = 3.5626e-4 * T^2 - 0.41704 * T + 156.8507 (TPa)^-1
# K_sec = 1e-7 (MPa*dpa)^-1
#
# where T is the temperature in K. The swelling strain rate is calculated from
# swelling increment for the timestep divided by the timestep size.
#
# Results for selected timesteps corresponding to a particular fluence value
# are provided in the table below
# BISON Analytical
# Fluence Swelling Swelling Strain Creep Strain Creep Strain
# (dpa) Strain (-) Rate (s^-1) Rate (s^-1) Rate (s^-1)
# 0.25 2.86074e-3 3.59778e-8 1.44079e-11 1.44079e-11
# 0.5 4.21743e-3 2.27542e-8 9.66361e-12 9.66362e-12
# 1.0 5.80817e-3 1.18073e-8 5.73616e-12 5.73615e-12
# 2.0 7.20081e-3 4.05640e-9 2.95533e-12 2.95533e-12
# 3.0 7.70534e-3 1.53843e-9 2.05195e-12 2.05195e-12
# 4.0 7.90061e-3 6.07271e-10 1.71787e-12 1.71787e-12
#
# Command line arguments are used to test the SINGH model at low temperatures as well as the KOYANAGI model, which
# has a slightly different form than described above. See the documentation for more details.
[GlobalParams]
displacements = 'disp_x disp_y'
order = FIRST
family = LAGRANGE
[]
[Mesh]
coord_type = RZ
[cylinder]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[temperature]
initial_condition = 1073.15
[]
[]
[AuxVariables]
[fast_neutron_fluence]
[]
[fast_neutron_flux]
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[creep_rate_aux]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
generate_output = vonmises_stress
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temperature
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temperature
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FastNeutronFluenceAux
fast_neutron_flux = fast_neutron_flux
variable = fast_neutron_fluence
block = 0
execute_on = 'timestep_begin'
[]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
block = 0
factor = 50e18
execute_on = timestep_begin
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[creep_rate]
type = MaterialRealAux
variable = creep_rate_aux
property = creep_rate
block = 0
execute_on = 'initial linear'
[]
[]
[BCs]
[temperature_all]
type = DirichletBC
boundary = 'left right bottom top'
variable = temperature
value = 1073.15
[]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[topY]
type = Pressure
variable = disp_y
boundary = top
factor = -3e6
[]
[]
[Materials]
[thermal]
type = HeatConductionMaterial
thermal_conductivity = 10.0
specific_heat = 1.0
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 2700.0
[]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temperature
swelling_model = KATOH
number_of_substeps = 100
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ComputeMultipleInelasticStress
block = 0
inelastic_models = composite_sic_creep
[]
[composite_sic_creep]
type = CompositeSiCCreepUpdate
fast_neutron_flux = fast_neutron_flux
block = 0
temperature = temperature
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart'
petsc_options_value = '101'
line_search = 'none'
nl_rel_tol = 1e-8
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
dt = 1e4
end_time = 1e6
[]
[Postprocessors]
[avg_temperature]
type = SideAverageValue
boundary = 'left right top bottom'
execute_on = 'initial timestep_end'
variable = temperature
[]
[creep_rate]
type = ElementAverageValue
block = 0
variable = creep_rate_aux
execute_on = 'initial timestep_end'
[]
[von_mises]
type = ElementAverageValue
variable = vonmises_stress
execute_on = 'initial timestep_end'
block = 0
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[fast_neutron_fluence]
type = AverageNodalVariableValue
variable = fast_neutron_fluence
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_mieloszyk_ad.i)
# The purpose of this test is to verify the implementation of the Mieloszyk
# swelling model for composite and monolithic silicon carbide. The algorithm
# used is from:
#
# A. J. Mieloszyk, "Assessing Thermo-Mechanical Performance of ThO2 and SiC Clad
# Light Water Reactor Fuel Rods with a Modular Simulation Tool," PhD Dissertation,
# Massachusetts Institute of Technology, 2015.
#
# In the above reference an example is provided that applies a temperature and
# fluence profile over 300 days (here we simplify that to 300 seconds with the
# same profiles). The reference also provides an expected profile of the
# volumetric swelling as a function of time given those profiles. The BISON
# results from this test are compared to the results provided in Mieloszyk's
# dissertation in the swelling_mieloszyk_testing.xlsx Excel file contained
# within this directory.
#
# The small discrepancies between the BISON and Mieloszyk curves can be
# attributed to the number of points used in the digitization of
# Mieloszyk's data.
#
# The results of predicted swelling at selected times used in the temperature
# function of the test are tabulated below:
#
# Time (s) Temperature (K) Fluence (dpa) Swelling Strain (-)
# 0 650.0 0 0
# 100 650.0 2.0 0.014267
# 120 350.0 2.0 0.014267
# 200 750.0 4.2 0.011403
# 250 650.0 4.9 0.014604
# 300 550.0 5.6 0.017932
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
order = FIRST
family = LAGRANGE
[]
[disp_y]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 650.0
[]
[]
[AuxVariables]
[fast_neutron_fluence]
order = FIRST
family = LAGRANGE
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
use_automatic_differentiation = true
[]
[]
[Kernels]
[heat]
type = ADHeatConduction
variable = temp
[]
[]
[Functions]
[fast_neutron_fluence_func]
type = PiecewiseLinear
x = '0 100.0 120.0 200.0 300.0'
y = '0 2e25 2e25 4.2e25 5.6e25'
[]
[temperature]
type = PiecewiseLinear
x = '0 100 100.01 120 120.01 200.0 200.01 250 250.01 300'
y = '650 650 350 350 675 750 600 650 625 550'
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence_func
execute_on = timestep_begin
[]
[total_swelling]
type = ADMaterialRealAux
variable = swell
property = swelling
block = 0
[]
[]
[BCs]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[rightX]
type = FunctionDirichletBC
variable = disp_x
boundary = right
function = '1e-6*t'
[]
[topY]
type = FunctionDirichletBC
variable = disp_y
boundary = top
function = '1e-6*t'
[]
[all_temp]
type = FunctionDirichletBC
variable = temp
boundary = 'bottom right top left'
function = temperature
[]
[]
[Materials]
[swelling]
type = ADCompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temp
swelling_model = MIELOSZYK
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ADComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ADComputeFiniteStrainElasticStress
block = 0
[]
[thermalCompositeSiC]
type = ADCompositeSiCThermal
temperature = temp
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = ADStrainAdjustedDensity
block = 0
strain_free_density = 3100.0
[]
[]
[Executioner]
type = Transient
solve_type = Newton
petsc_options = '-snes_ksp_ew '
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = none
l_max_its = 60
nl_rel_tol = 1e-8
nl_abs_tol = 1e-3
l_tol = 1e-5
start_time = 0.0
dt = 1
end_time = 300.0
[]
[Postprocessors]
[avgFuelTemp]
type = ElementAverageValue
variable = temp
execute_on = 'initial timestep_end'
block = 0
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[average_fluence]
type = ElementAverageValue
variable = fast_neutron_fluence
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(test/tests/monolithicSiCThermal/thermal_stone_unirradiated.i)
# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500 K
# over 12.0 s.
#
# Specific heat capacity is computed using the correlation from:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# and thermal conductivity is computed using the thermal conductivity correlation
# from:
#
# J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, C. A.
# Back, "Stress analysis and probablistic assessment of multi-layer SiC-based
# accident tolerant nuclear fuel cladding," Journal of Nuclear Materials, 466,
# 2015, p. 682-697.
#
# The specific heat correlation is tested elsewhere. The results of the BISON
# simulation for unirradiatedthermal conductivity is compared to the analytical
# solution below.
#
# SiC/SiC k SiC/SiC k
# Temp(K) BISON(W/m-K) analytical(W/m-K)
# 300 100.951 100.951
# 600 68.9080 68.9080
# 900 50.9770 50.9770
# 1200 41.1640 41.1640
# 1500 33.4750 33.4750
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 3
[]
[]
[Variables]
[temperature]
order = FIRST
family = LAGRANGE
initial_condition = 300.0
[]
[]
[AuxVariables]
[specific_heat]
order = CONSTANT
family = MONOMIAL
[]
[thermal_conductivity]
order = CONSTANT
family = MONOMIAL
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temperature
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temperature
[]
[]
[AuxKernels]
[thermal_conductivity]
type = MaterialRealAux
variable = thermal_conductivity
property = thermal_conductivity
block = 0
execute_on = 'initial linear'
[]
[specific_heat]
type = MaterialRealAux
variable = specific_heat
property = specific_heat
block = 0
execute_on = 'initial linear'
[]
[]
[Functions]
[temp_ramp]
type = PiecewiseLinear
x = '0.0 12.0'
y = '300 1500'
[]
[]
[BCs]
[temperature_all]
type = FunctionDirichletBC
boundary = 'left right front back bottom top' # All cube faces
variable = temperature
function = temp_ramp
[]
[]
[Materials]
[MonolithicSiCThermal]
type = MonolithicSiCThermal
temperature = temperature
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = ParsedMaterial
block = 0
property_name = density
expression = 3210.0
[]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temperature
fast_neutron_fluence = 0.0
eigenstrain_name = swelling_strain
[]
[]
[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-6
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
num_steps = 12
dt = 1.0
[]
[Postprocessors]
[avg_temperature]
type = SideAverageValue
boundary = 'left right top bottom front back'
execute_on = 'initial timestep_end'
variable = temperature
[]
[avg_th_cond]
type = ElementAverageValue
block = 0
variable = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[avg_specific_heat]
type = ElementAverageValue
block = 0
variable = specific_heat
execute_on = 'initial timestep_end'
[]
[]
[Outputs]
exodus = true
[]
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh_ad.i)
# This test tests the Katoh swelling model used for composite and monolithic
# silicon carbide. The model is obtained from:
#
# Y. Katoh, K. Ozawa, C. Shih, T. Nozawa, R. J. Shinavski, A. Hasegawa, and
# L. L. Snead. Continuous SiC fiber, CVI SiC matrix composites for nuclear
# applications: properties and irradiation effects. Journal of Nuclear Materials,
# 448, pp. 448-476, 2014.
#
# The swelling rate is a function of the fast neutron fluence:
#
# dS/dgamma = k(T) * gamma^(-1/3) * exp(-gamma/gamma_sc(T))
#
# where gamma is the fluence in dpa, S is the swelling strain, and k(T) and
# gamma_sc(T) are temperature dependent constants. The above equation cannot
# be analytically integrated and an incremental approach is used to calculate the
# swelling strain. The "analytical" solution is taken as the results of the
# numerical integral calculated using Wolfram Alpha or Mathematica. Depending
# upon the number of substeps used the BISON simulations approach this analytical
# solution.
#
# Here, a single 1x1 m 2D element is subjected to a constant temperature of
# 1073.15 K with a varying fast neutron fluence from 0 to 1e25 n/m^2 (0 to 5 dpa)
# over 1 second. Below BISON results are included for 100, 1000 and 10000 substeps to
# illustrate that the results approach that of the analytical solution. This test
# is for the 100 substep case. Tests for the 1000 and 10000 substeps at high
# temperature as well as low and mid range temperatures of 473.15 K
# and 773.15 K (with 100 substeps) are included using cli_args in the tests file.
#
# Fluence BISON Swelling Strain (-) Analytical
# (dpa) 100 substeps 1000 substeps 10000 substeps Swelling
# 1.0 5.91955e-3 5.83353e-3 5.83899e-3 5.84045e-3
# 2.0 7.20081e-3 7.22617e-3 7.23163e-3 7.23270e-3
# 3.0 7.70534e-3 7.73070e-3 7.73617e-3 7.73713e-3
# 4.0 7.90061e-3 7.92597e-3 7.93143e-3 7.93235e-3
# 5.0 7.97852e-3 8.00388e-3 8.00934e-3 8.01025e-3
#
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
order = FIRST
family = LAGRANGE
[]
[disp_y]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1073.15
[]
[]
[AuxVariables]
[fast_neutron_fluence]
order = FIRST
family = LAGRANGE
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
use_automatic_differentiation = true
[]
[]
[Kernels]
[heat]
type = ADHeatConduction
variable = temp
[]
[heat_ie]
type = ADHeatConductionTimeDerivative
variable = temp
[]
[]
[Functions]
[fast_neutron_fluence_func]
type = PiecewiseLinear
x = '0 1.0'
y = '0 5e25'
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence_func
execute_on = timestep_begin
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[]
[BCs]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[all_temp]
type = DirichletBC
variable = temp
boundary = 'bottom left right top'
value = 1073.15
[]
[]
[Materials]
[swelling]
type = ADCompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temp
swelling_model = KATOH
number_of_substeps = 100
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ADComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ADComputeFiniteStrainElasticStress
block = 0
[]
[thermalCompositeSiC]
type = ADCompositeSiCThermal
temperature = temp
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = ADStrainAdjustedDensity
block = 0
density = 3100.0
[]
[]
[Executioner]
type = Transient
solve_type = Newton
petsc_options = '-snes_ksp_ew '
petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type'
petsc_options_value = '70 hypre boomeramg'
l_max_its = 60
nl_rel_tol = 1e-5
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
dt = 0.01
end_time = 1.0
[]
[Postprocessors]
[avgFuelTemp]
type = ElementAverageValue
variable = temp
execute_on = 'initial timestep_end'
block = 0
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[average_fluence]
type = ElementAverageValue
variable = fast_neutron_fluence
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(test/tests/solid_mechanics/sic_creep/ad_composite_sic_creep.i)
# The mesh is a 1x1 square (single element) assuming RZ axisymmetry.
# The temperature is held constant at 1073.15 K with a varying fast fluence
# from 0 to 1e25 n/m^2 over 1e6 second. A timestep of 1e4 s is taken and
# swelling is calculated by the Katoh volumetric swelling model with 100 substeps.
# An axial pressure of -3MPa is applied (i.e., pulling outward from the surface).
#
# The creep rate is calculated from:
#
# e_dot = e_dot_pri+e_dot_sec
# e_dot_pri = K_pri * sigma * e_dot_swell
# e_dot_sec = K_sec * sigma * phi
#
# where K_pri and K_sec are conversion factors, sigma is the Von Mises stress in TPa, phi is the fast neutron flux, and
# e_dot_swell is the swelling strain rate. The conversion factors are given by
#
# K_pri = 3.5626e-4 * T^2 - 0.41704 * T + 156.8507 (TPa)^-1
# K_sec = 1e-7 (MPa*dpa)^-1
#
# where T is the temperature in K. The swelling strain rate is calculated from
# swelling increment for the timestep divided by the timestep size.
#
# Results for selected timesteps corresponding to a particular fluence value
# are provided in the table below
# BISON Analytical
# Fluence Swelling Swelling Strain Creep Strain Creep Strain
# (dpa) Strain (-) Rate (s^-1) Rate (s^-1) Rate (s^-1)
# 0.25 2.86074e-3 3.59778e-8 1.44079e-11 1.44079e-11
# 0.5 4.21743e-3 2.27542e-8 9.66361e-12 9.66362e-12
# 1.0 5.80817e-3 1.18073e-8 5.73616e-12 5.73615e-12
# 2.0 7.20081e-3 4.05640e-9 2.95533e-12 2.95533e-12
# 3.0 7.70534e-3 1.53843e-9 2.05195e-12 2.05195e-12
# 4.0 7.90061e-3 6.07271e-10 1.71787e-12 1.71787e-12
#
# Command line arguments are used to test the SINGH model at low temperatures as well as the KOYANAGI model, which
# has a slightly different form than described above. See the documentation for more details.
[GlobalParams]
displacements = 'disp_x disp_y'
order = FIRST
family = LAGRANGE
[]
[Mesh]
coord_type = RZ
[cylinder]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[temperature]
initial_condition = 1073.15
[]
[]
[AuxVariables]
[fast_neutron_fluence]
[]
[fast_neutron_flux]
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[creep_rate_aux]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
generate_output = vonmises_stress
use_automatic_differentiation = true
[]
[]
[Kernels]
[heat]
type = ADHeatConduction
variable = temperature
[]
[heat_ie]
type = ADHeatConductionTimeDerivative
variable = temperature
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FastNeutronFluenceAux
fast_neutron_flux = fast_neutron_flux
variable = fast_neutron_fluence
block = 0
execute_on = 'timestep_begin'
[]
[fast_neutron_flux]
type = FastNeutronFluxAux
variable = fast_neutron_flux
block = 0
factor = 50e18
execute_on = timestep_begin
[]
[total_swelling]
type = ADMaterialRealAux
variable = swell
property = swelling
block = 0
[]
[creep_rate]
type = ADMaterialRealAux
variable = creep_rate_aux
property = creep_rate
block = 0
execute_on = 'initial linear'
[]
[]
[BCs]
[temperature_all]
type = ADDirichletBC
boundary = 'left right bottom top'
variable = temperature
value = 1073.15
[]
[leftX]
type = ADDirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = ADDirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[topY]
type = ADPressure
variable = disp_y
boundary = top
factor = -3e6
[]
[]
[Materials]
[thermal]
type = ADHeatConductionMaterial
thermal_conductivity = 10.0
specific_heat = 1.0
[]
[density]
type = ADStrainAdjustedDensity
block = 0
strain_free_density = 2700.0
[]
[swelling]
type = ADCompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temperature
swelling_model = KATOH
number_of_substeps = 100
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor]
type = ADComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ADComputeMultipleInelasticStress
block = 0
inelastic_models = composite_sic_creep
[]
[composite_sic_creep]
type = ADCompositeSiCCreepUpdate
fast_neutron_flux = fast_neutron_flux
block = 0
temperature = temperature
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart'
petsc_options_value = '101'
line_search = 'none'
nl_rel_tol = 1e-8
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
dt = 1e4
end_time = 1e6
[]
[Postprocessors]
[avg_temperature]
type = SideAverageValue
boundary = 'left right top bottom'
execute_on = 'initial timestep_end'
variable = temperature
[]
[creep_rate]
type = ElementAverageValue
block = 0
variable = creep_rate_aux
execute_on = 'initial timestep_end'
[]
[von_mises]
type = ElementAverageValue
variable = vonmises_stress
execute_on = 'initial timestep_end'
block = 0
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[fast_neutron_fluence]
type = AverageNodalVariableValue
variable = fast_neutron_fluence
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(test/tests/monolithicSiCThermal/thermal_stone_irradiated.i)
# The mesh is a 1x1 square (single element).
# The temperature is held constant at 1073.15 K with a varying fast fluence
# from 0 to 1e25 n/m^2 over 1 second. The swelling is calculated by the Katoh
# volumetric swelling model. At 1073.15 K the calculated unirradiated thermal
# conductivity using the model of:
#
# J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, C. A.
# Back, "Stress analysis and probablistic assessment of multi-layer SiC-based
# accident tolerant nuclear fuel cladding," Journal of Nuclear Materials, 466,
# 2015, p. 682-697.
#
# is 44.7072 W/m-K. Then, the degradation due to swelling for monolithic silicon
# carbide is given by:
#
# 1 / k_irradiated = R_0 + R_irr
#
# where R_0 = 1 / k_unirradiated and R_irr = 6.08 * S where S is the swelling
# strain. The BISON results are compared to the analytical using the BISON
# values for swelling.
#
# Fluence Swelling Monolithic SiC k Monolithic SiC k
# (dpa) Strain (-) BISON(W/m-K) analytical(W/m-K)
# 1.0 5.91955e-3 17.1354 17.1354
# 2.0 7.20081e-3 15.1175 15.1175
# 3.0 7.70534e-3 14.4475 14.4475
# 4.0 7.90061e-3 14.2039 14.2039
# 5.0 7.97852e-3 14.1089 14.1089
[GlobalParams]
displacements = 'disp_x disp_y'
order = FIRST
family = LAGRANGE
[]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[temperature]
initial_condition = 1073.15
[]
[]
[AuxVariables]
[fast_neutron_fluence]
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[thermal_conductivity]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temperature
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temperature
[]
[]
[AuxKernels]
[thermal_conductivity]
type = MaterialRealAux
variable = thermal_conductivity
property = thermal_conductivity
block = 0
execute_on = 'initial linear'
[]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence
block = 0
execute_on = 'timestep_begin'
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[]
[Functions]
[fast_neutron_fluence]
type = PiecewiseLinear
x = '0 1.0'
y = '0 5e25'
[]
[]
[BCs]
[temperature_all]
type = DirichletBC
boundary = 'left right bottom top'
variable = temperature
value = 1073.15
[]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[]
[Materials]
[MonolithicSiCThermal]
type = MonolithicSiCThermal
temperature = temperature
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = ParsedMaterial
block = 0
property_name = density
expression = 3210.0
[]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temperature
swelling_model = KATOH
number_of_substeps = 100
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
block = 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'
nl_rel_tol = 1e-5
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
dt = 0.01
end_time = 1.0
[]
[Postprocessors]
[avg_temperature]
type = SideAverageValue
boundary = 'left right top bottom'
execute_on = 'initial timestep_end'
variable = temperature
[]
[avg_th_cond]
type = ElementAverageValue
block = 0
variable = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(assessment/LWR/validation/ATF_SiC_HFIR/analysis/base.i)
[GlobalParams]
temperature = temperature
displacements = 'disp_x disp_y'
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
[]
[Mesh]
coord_type = RZ
[clad_mesh]
type = GeneratedMeshGenerator
dim = 2
nx = 10
ny = 100
xmin = 3.51e-3
xmax = 4.27e-3
ymin = 0.0
ymax = 16.04e-3
[]
[]
[Variables]
[temperature]
initial_condition = 295.0 # set initial temperature to coolant inlet
[]
[]
[AuxVariables]
[fast_neutron_flux]
order = CONSTANT
family = MONOMIAL
[]
[fast_neutron_fluence]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear
x = '0 1e4 2.16e6 2.3328e6'
y = '0 333333.0 333333.0 0'
scale_factor = 1
[]
[axial_peaking_factors]
type = ParsedFunction
expression = 1
[]
[pressure_ramp]
type = PiecewiseLinear
x = '-200 0 1e8'
y = '6.537e-3 1 1'
[]
[max_time_step_size_func]
type = ParsedFunction
expression = 'if(t < 2.16e6, 5e4, 1e4)'
[]
[heat_flux_func]
type = ParsedFunction
expression = 'if(t < 2.16e6, 6e5, 0)'
[]
[]
[Constraints]
[top]
type = EqualValueBoundaryConstraint
variable = disp_y
secondary = 'top'
penalty = 1e+11
formulation = penalty
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[clad]
add_variables = true
strain = FINITE
eigenstrain_names = 'composite_thermal_strain composite_swelling_strain'
automatic_eigenstrain_names = true
generate_output = 'hoop_stress hoop_strain stress_xx stress_yy stress_zz vonmises_stress strain_zz strain_yy'
cylindrical_axis_point1 = '0 0 0'
cylindrical_axis_point2 = '0 1 0'
extra_vector_tags = 'ref'
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temperature
extra_vector_tags = 'ref'
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temperature
extra_vector_tags = 'ref'
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = MaterialRealAux
variable = fast_neutron_flux
property = fast_neutron_flux
execute_on = timestep_begin
[]
[fast_neutron_fluence]
type = MaterialRealAux
variable = fast_neutron_fluence
property = fast_neutron_fluence
execute_on = timestep_begin
[]
[]
[BCs]
[no_y]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[]
[heat_flux_clad_inner_surface]
type = FunctionNeumannBC
boundary = left
variable = temperature
function = heat_flux_func
[]
[right]
type = ConvectiveHeatFluxBC
variable = temperature
boundary = 'right'
T_infinity = 295.0
heat_transfer_coefficient = 12000
heat_transfer_coefficient_dT = 0
[]
[]
[Materials]
[flux]
type = FastNeutronFlux
calculate_fluence = true
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
factor = 3e13
[]
[composite_thermal]
type = CompositeSiCThermal
thermal_conductivity_model = STONE
temperature = temperature
outputs = 'all'
[]
[composite_density]
type = StrainAdjustedDensity
strain_free_density = 2700.0
[]
[composite_elasticity_tensor]
type = CompositeSiCElasticityTensor
elasticity_model = SINGH
[]
[composite_stress]
type = ComputeStrainIncrementBasedStress
[]
[composite_thermal_expansion]
type = CompositeSiCThermalExpansionEigenstrain
stress_free_temperature = 295.0
temperature = temperature
eigenstrain_name = composite_thermal_strain
outputs = 'all'
[]
[composite_irradiation_swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
temperature = temperature
fast_neutron_fluence = fast_neutron_fluence
swelling_model = KATOH
number_of_substeps = 1000
eigenstrain_name = composite_swelling_strain
outputs = 'all'
[]
[]
[Dampers]
[limitT]
type = MaxIncrement
max_increment = 50.0
variable = temperature
[]
[limitX]
type = MaxIncrement
max_increment = 1e-5
variable = disp_x
[]
[]
[PerformanceMetricOutputs]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON' #PJFNK'
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = 'none'
l_max_its = 50
l_tol = 8e-3
nl_max_its = 50
nl_rel_tol = 1e-4
nl_abs_tol = 1e-7
start_time = -200
n_startup_steps = 1
end_time = 2.3328e6
dtmin = 1
[TimeStepper]
type = IterationAdaptiveDT
dt = 2e2
optimal_iterations = 25
iteration_window = 5
growth_factor = 2
cutback_factor = .5
timestep_limiting_postprocessor = max_time_step_size_pp
[]
[]
[Postprocessors]
[hoop_stress_outer]
type = ElementalVariableValue
variable = hoop_stress
elementid = 499
[]
[axial_stress_outer]
type = ElementalVariableValue
variable = stress_yy
elementid = 499
[]
[hoop_stress_inner]
type = ElementalVariableValue
variable = hoop_stress
elementid = 490
[]
[axial_stress_inner]
type = ElementalVariableValue
variable = stress_yy
elementid = 490
[]
[hoop_strain_inner]
type = ElementalVariableValue
variable = strain_zz
elementid = 490
[]
[axial_strain_inner]
type = ElementalVariableValue
variable = strain_yy
elementid = 490
[]
[hoop_strain_outer]
type = ElementalVariableValue
variable = strain_zz
elementid = 499
[]
[axial_strain_outer]
type = ElementalVariableValue
variable = strain_yy
elementid = 499
[]
[max_time_step_size_pp]
type = FunctionValuePostprocessor
function = max_time_step_size_func
execute_on = 'INITIAL TIMESTEP_END'
[]
[outlet_heat_flux]
type = SideDiffusiveFluxIntegral
variable = temperature
boundary = left
diffusivity = thermal_conductivity
[]
[swelling_strain_outer]
type = ElementalVariableValue
variable = composite_swelling_strain_22
elementid = 499
[]
[swelling_strain_inner]
type = ElementalVariableValue
variable = composite_swelling_strain_22
elementid = 490
[]
[thermal_strain_outer]
type = ElementalVariableValue
variable = composite_thermal_strain_22
elementid = 499
[]
[thermal_strain_inner]
type = ElementalVariableValue
variable = composite_thermal_strain_22
elementid = 490
[]
[]
[Outputs]
perf_graph = true
exodus = true
[csv]
type = CSV
execute_on = final
[]
color = false
[]
(test/tests/thermalCompositeSiC/thermal_stone_irradiated.i)
# The mesh is a 1x1 square (single element).
# The temperature is held constant at 1073.15 K with a varying fast fluence
# from 0 to 1e25 n/m^2 over 1 second. The swelling is calculated by the Katoh
# volumetric swelling model. At 1073.15 K the calculated unirradiated thermal
# conductivity using the model of:
#
# J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, C. A.
# Back, "Stress analysis and probablistic assessment of multi-layer SiC-based
# accident tolerant nuclear fuel cladding," Journal of Nuclear Materials, 466,
# 2015, p. 682-697.
#
# is 14.9090 W/m-K. Then, the degradation due to swelling for composite silicon
# carbide is given by:
#
# 1 / k_irradiated = R_0 + R_irr
#
# where R_0 = 1 / k_unirradiated and R_irr = 15.11 * S where S is the swelling
# strain. The BISON results are compared to the analytical using the BISON
# values for swelling.
#
# Fluence Swelling SiC/SiC k SiC/SiC k
# (dpa) Strain (-) BISON(W/m-K) analytical(W/m-K)
# 1.0 5.91955e-3 6.38904 6.38904
# 2.0 7.20081e-3 5.68577 5.68577
# 3.0 7.70534e-3 5.44955 5.44955
# 4.0 7.90061e-3 5.36332 5.36332
# 5.0 7.97852e-3 5.32967 5.32967
[GlobalParams]
displacements = 'disp_x disp_y'
order = FIRST
family = LAGRANGE
[]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[temperature]
initial_condition = 1073.15
[]
[]
[AuxVariables]
[fast_neutron_fluence]
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[thermal_conductivity]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temperature
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temperature
[]
[]
[AuxKernels]
[thermal_conductivity]
type = MaterialRealAux
variable = thermal_conductivity
property = thermal_conductivity
block = 0
execute_on = 'initial linear'
[]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence
block = 0
execute_on = 'timestep_begin'
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[]
[Functions]
[fast_neutron_fluence]
type = PiecewiseLinear
x = '0 1.0'
y = '0 5e25'
[]
[]
[BCs]
[temperature_all]
type = DirichletBC
boundary = 'left right bottom top'
variable = temperature
value = 1073.15
[]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[]
[Materials]
[thermalCompositeSiC]
type = CompositeSiCThermal
temperature = temperature
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 2700.0
[]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temperature
swelling_model = KATOH
number_of_substeps = 100
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
block = 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'
nl_rel_tol = 1e-5
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
dt = 0.01
end_time = 1.0
[]
[Postprocessors]
[avg_temperature]
type = SideAverageValue
boundary = 'left right top bottom'
execute_on = 'initial timestep_end'
variable = temperature
[]
[avg_th_cond]
type = ElementAverageValue
block = 0
variable = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_mieloszyk.i)
# The purpose of this test is to verify the implementation of the Mieloszyk
# swelling model for composite and monolithic silicon carbide. The algorithm
# used is from:
#
# A. J. Mieloszyk, "Assessing Thermo-Mechanical Performance of ThO2 and SiC Clad
# Light Water Reactor Fuel Rods with a Modular Simulation Tool," PhD Dissertation,
# Massachusetts Institute of Technology, 2015.
#
# In the above reference an example is provided that applies a temperature and
# fluence profile over 300 days (here we simplify that to 300 seconds with the
# same profiles). The reference also provides an expected profile of the
# volumetric swelling as a function of time given those profiles. The BISON
# results from this test are compared to the results provided in Mieloszyk's
# dissertation in the swelling_mieloszyk_testing.xlsx Excel file contained
# within this directory.
#
# The small discrepancies between the BISON and Mieloszyk curves can be
# attributed to the number of points used in the digitization of
# Mieloszyk's data.
#
# The results of predicted swelling at selected times used in the temperature
# function of the test are tabulated below:
#
# Time (s) Temperature (K) Fluence (dpa) Swelling Strain (-)
# 0 650.0 0 0
# 100 650.0 2.0 0.014267
# 120 350.0 2.0 0.014267
# 200 750.0 4.2 0.011403
# 250 650.0 4.9 0.014604
# 300 550.0 5.6 0.017932
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
order = FIRST
family = LAGRANGE
[]
[disp_y]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 650.0
[]
[]
[AuxVariables]
[fast_neutron_fluence]
order = FIRST
family = LAGRANGE
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[]
[Functions]
[fast_neutron_fluence_func]
type = PiecewiseLinear
x = '0 100.0 120.0 200.0 300.0'
y = '0 2e25 2e25 4.2e25 5.6e25'
[]
[temperature]
type = PiecewiseLinear
x = '0 100 100.01 120 120.01 200.0 200.01 250 250.01 300'
y = '650 650 350 350 675 750 600 650 625 550'
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence_func
execute_on = timestep_begin
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[]
[BCs]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[rightX]
type = FunctionDirichletBC
variable = disp_x
boundary = right
function = '1e-6*t'
[]
[topY]
type = FunctionDirichletBC
variable = disp_y
boundary = top
function = '1e-6*t'
[]
[all_temp]
type = FunctionDirichletBC
variable = temp
boundary = 'bottom right top left'
function = temperature
[]
[]
[Materials]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temp
swelling_model = MIELOSZYK
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
block = 0
[]
[thermalCompositeSiC]
type = CompositeSiCThermal
temperature = temp
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 3100.0
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-snes_ksp_ew '
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = none
l_max_its = 60
nl_rel_tol = 1e-8
nl_abs_tol = 1e-3
l_tol = 1e-5
start_time = 0.0
dt = 1
end_time = 300.0
[]
[Postprocessors]
[avgFuelTemp]
type = ElementAverageValue
variable = temp
execute_on = 'initial timestep_end'
block = 0
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[average_fluence]
type = ElementAverageValue
variable = fast_neutron_fluence
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(examples/accident_tolerant_fuel/u3si2_sic/u3si2_outer_monolith_1.5D.i)
# Model is of a 10 pellet fuel rodlet modeled in 1.5D. The rodlet contains
# U3Si2 fuel and a multilayer silicon carbide cladding (an inner composite
# winding layer) and an outer monolithic layer. The inner composite layer is
# 0.75 mm thick and the outer monolithic layer is 0.25 mm thick. The internal
# layered1D mesh generator can model a clad an arbitrary number of additional blocks.
# Therefore, to create the multilayer SiC clad the composite layer is assigned to the
# clad block and the monolithic_layer is assigned to the monolithic_layer block
# as specified in the additional_block_names parameter in the Mesh block.
initial_fuel_density = 11590.0
[GlobalParams]
density = ${initial_fuel_density}
initial_porosity = 0.05
order = SECOND
family = LAGRANGE
energy_per_fission = 3.2e-11 # J/fission
displacements = disp_x
temperature = temperature
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
[]
[Mesh]
coord_type = RZ
[layered1D_mesh]
type = Layered1DMeshGenerator
slices_per_block = 10
clad_gap_width = 8.0e-5
clad_mesh_density = customize
clad_thickness = 0.00075
nx_c = 5
additional_block_names = 'monolithic_layer'
additional_elements_per_ring = '3'
additional_ring_thicknesses = '0.00025'
fuel_height = 0.1186
plenum_height = 0.027
[]
patch_update_strategy = auto
partitioner = centroid
centroid_partitioner_direction = y
[]
[UserObjects]
[pin_geometry]
type = Layered1DFuelPinGeometry
mesh_generator = layered1D_mesh
[]
[]
[Variables]
[temperature]
initial_condition = 580.0 # set initial temperature to coolant inlet
[]
[]
[AuxVariables]
[disp_y] ## Required for easier visualization in Paraview
[]
[disp_z] ## Required for easier visualization in Paraview
[]
[fast_neutron_flux]
block = 'clad monolithic_layer'
order = CONSTANT
family = MONOMIAL
[]
[fast_neutron_fluence]
block = 'clad monolithic_layer'
order = CONSTANT
family = MONOMIAL
[]
[grain_radius]
block = fuel
initial_condition = 10e-6
[]
[solid_swell]
order = CONSTANT
family = MONOMIAL
block = fuel
[]
[gaseous_swelling]
order = CONSTANT
family = MONOMIAL
block = fuel
[]
[densification]
order = CONSTANT
family = MONOMIAL
block = fuel
[]
[volumetric_swelling_strain]
order = CONSTANT
family = MONOMIAL
block = fuel
[]
[relocation]
order = CONSTANT
family = MONOMIAL
block = fuel
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear
x = '0 1e4 1e8'
y = '0 25000 25000'
scale_factor = 1
[]
[axial_peaking_factors]
type = ParsedFunction
expression = 1
[]
[pressure_ramp]
type = PiecewiseLinear
x = '-200 0 1e8'
y = '6.537e-3 1 1'
[]
[q]
type = CompositeFunction
functions = 'power_history axial_peaking_factors'
[]
[clad_axial_pressure]
type = CladdingAxialPressureFunction
plenum_pressure = plenum_pressure
coolant_pressure = pressure_ramp
coolant_pressure_scaling_factor = 15.5e6
fuel_pin_geometry = pin_geometry
[]
[fuel_axial_pressure]
type = ParsedFunction
expression = plenum_pressure
symbol_names = plenum_pressure
symbol_values = plenum_pressure
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temperature
extra_vector_tags = 'ref'
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temperature
extra_vector_tags = 'ref'
[]
[heat_source]
type = NeutronHeatSource
variable = temperature
block = fuel
burnup_function = burnup
extra_vector_tags = 'ref'
[]
[]
[Physics]
[SolidMechanics]
[Layered1D]
[fuel]
block = fuel
add_variables = true
strain = SMALL
incremental = true
add_scalar_variables = true
out_of_plane_strain_name = strain_yy
fuel_pin_geometry = pin_geometry
out_of_plane_pressure_function = fuel_axial_pressure
eigenstrain_names = 'fuel_thermal_strain fuel_swelling_strain'
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
extra_vector_tags = 'ref'
group_scalar_vars_in_reference_residual = true
mesh_generator = layered1D_mesh
[]
[composite]
block = clad
add_variables = true
strain = SMALL
incremental = true
add_scalar_variables = true
out_of_plane_strain_name = strain_yy
fuel_pin_geometry = pin_geometry
out_of_plane_pressure_function = clad_axial_pressure
eigenstrain_names = 'composite_thermal_strain composite_swelling_strain'
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
extra_vector_tags = 'ref'
group_scalar_vars_in_reference_residual = true
mesh_generator = layered1D_mesh
[]
[monolith]
block = monolithic_layer
add_variables = true
strain = SMALL
incremental = true
add_scalar_variables = true
out_of_plane_strain_name = strain_yy
fuel_pin_geometry = pin_geometry
out_of_plane_pressure_function = clad_axial_pressure
eigenstrain_names = 'monolith_thermal_strain monolith_swelling_strain'
generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
extra_vector_tags = 'ref'
group_scalar_vars_in_reference_residual = true
mesh_generator = layered1D_mesh
[]
[]
[]
[]
[Burnup]
[burnup]
block = fuel
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
num_radial = 80
num_axial = 11
order = CONSTANT
family = MONOMIAL
fuel_pin_geometry = pin_geometry
fuel_volume_ratio = 1.0
RPF = RPF
fuel_type = U3Si2
[]
[]
[AuxKernels]
[fast_neutron_flux]
type = MaterialRealAux
variable = fast_neutron_flux
property = fast_neutron_flux
block = 'clad monolithic_layer'
execute_on = timestep_begin
[]
[fast_neutron_fluence]
type = MaterialRealAux
variable = fast_neutron_fluence
property = fast_neutron_fluence
block = 'clad monolithic_layer'
execute_on = timestep_begin
[]
[grain_radius]
type = GrainRadiusAux
block = fuel
variable = grain_radius
temperature = temperature
execute_on = linear
[]
[solid_swell]
type = MaterialRealAux
variable = solid_swell
property = solid_swelling
execute_on = timestep_end
block = fuel
[]
[gas_swell]
type = MaterialRealAux
variable = gaseous_swelling
property = gaseous_swelling
execute_on = timestep_end
block = fuel
[]
[densification]
type = MaterialRealAux
variable = densification
property = densification
execute_on = timestep_end
block = fuel
[]
[volumetric_swelling_strain]
type = MaterialRealAux
variable = volumetric_swelling_strain
property = volumetric_swelling_strain
execute_on = timestep_end
block = fuel
[]
[]
[Contact]
[pellet_clad_mechanical]
primary = 5
secondary = 10
formulation = kinematic
model = frictionless
penalty = 1e7
[]
[]
[ThermalContact]
[thermal_contact]
type = GasGapHeatTransfer
variable = temperature
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
[]
[Pressure] # apply coolant pressure on clad outer walls
[coolantPressure]
use_displaced_mesh = false
boundary = 2
function = pressure_ramp # use the pressure_ramp function defined above
factor = 15.5e6
[]
[]
[PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
[plenumPressure]
use_displaced_mesh = false
boundary = 9
initial_pressure = 2.0e6
startup_time = 0
R = 8.314
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
[]
[]
[]
[CoolantChannel]
[convective_clad_surface] # apply convective boundary to clad outer surface
variable = temperature
boundary = 2
inlet_temperature = 580 # K
inlet_pressure = 15.5e6 # Pa
inlet_massflux = 3800 # kg/m^2-sec
rod_diameter = 10.368e-3 # m
rod_pitch = 1.26e-2 # m
linear_heat_rate = power_history
axial_power_profile = axial_peaking_factors
[]
[]
[Materials]
[flux]
type = FastNeutronFlux
calculate_fluence = true
block = 'clad monolithic_layer'
rod_ave_lin_pow = power_history
axial_power_profile = axial_peaking_factors
factor = 3e13
[]
### U3Si2 Fuel
[fuel_thermal]
type = SilicideFuelThermal
block = fuel
thermal_conductivity_model = WHITE
temperature = temperature
[]
[fuel_elasticity_tensor]
type = U3Si2ElasticityTensor
block = fuel
[]
[fuel_stress]
type = ComputeMultipleInelasticStress
block = fuel
tangent_operator = elastic
inelastic_models = 'fuel_creep'
[]
[fuel_creep]
type = U3Si2CreepUpdate
block = fuel
temperature = temperature
[]
[fuel_thermal_expansion]
type = U3Si2ThermalExpansionEigenstrain
block = fuel
temperature = temperature
stress_free_temperature = 295.0
eigenstrain_name = fuel_thermal_strain
[]
[fuel_volumetric_swelling]
type = U3Si2VolumetricSwellingEigenstrain
block = fuel
gaseous_swelling_type = FINLAY
temperature = temperature
burnup_function = burnup
eigenstrain_name = fuel_swelling_strain
[]
[fission_gas_release]
type = U3Si2Sifgrs
block = fuel
temperature = temperature
burnup_function = burnup
grain_radius_const = 2.5e-05
[]
[fuel_density]
type = StrainAdjustedDensity
block = fuel
strain_free_density = ${initial_fuel_density}
[]
### Composite SiC
[composite_thermal]
type = CompositeSiCThermal
thermal_conductivity_model = STONE
temperature = temperature
block = clad
[]
[composite_density]
type = StrainAdjustedDensity
block = clad
strain_free_density = 2700.0
[]
[composite_elasticity_tensor]
type = CompositeSiCElasticityTensor
block = clad
[]
[composite_stress]
type = ComputeStrainIncrementBasedStress
block = clad
[]
[composite_thermal_expansion]
type = CompositeSiCThermalExpansionEigenstrain
block = clad
stress_free_temperature = 295.0
temperature = temperature
eigenstrain_name = composite_thermal_strain
[]
[composite_irradiation_swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = clad
temperature = temperature
fast_neutron_fluence = fast_neutron_fluence
swelling_model = KATOH
number_of_substeps = 1000
eigenstrain_name = composite_swelling_strain
[]
### Monolithic SiC
[monolith_thermal]
type = MonolithicSiCThermal
temperature = temperature
thermal_conductivity_model = STONE
block = monolithic_layer
[]
[monolith_density]
type = StrainAdjustedDensity
block = monolithic_layer
strain_free_density = 3120.0
[]
[monolith_elasticity_tensor]
type = MonolithicSiCElasticityTensor
block = monolithic_layer
[]
[monolith_stress]
type = ComputeMultipleInelasticStress
block = monolithic_layer
tangent_operator = elastic
inelastic_models = 'monolith_creep'
[]
[monolith_creep]
type = MonolithicSiCCreepUpdate
block = monolithic_layer
fast_neutron_flux = fast_neutron_flux
temperature = temperature
k_function = 2e-37
[]
[monolith_thermal_expansion]
type = MonolithicSiCThermalExpansionEigenstrain
block = monolithic_layer
stress_free_temperature = 295.0
temperature = temperature
eigenstrain_name = monolith_thermal_strain
[]
[monolith_irradiation_swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = monolithic_layer
temperature = temperature
fast_neutron_fluence = fast_neutron_fluence
swelling_model = KATOH
number_of_substeps = 1000
eigenstrain_name = monolith_swelling_strain
[]
[]
[Dampers]
[limitT]
type = MaxIncrement
max_increment = 100.0
variable = temperature
[]
[limitX]
type = MaxIncrement
max_increment = 1e-5
variable = disp_x
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = 'none'
l_max_its = 50
l_tol = 8e-3
nl_max_its = 50
nl_rel_tol = 1e-4
nl_abs_tol = 1e-7
start_time = -200
n_startup_steps = 1
end_time = 8.0e7
dtmax = 1e6
dtmin = 1
[TimeStepper]
type = IterationAdaptiveDT
dt = 2e2
optimal_iterations = 25
iteration_window = 5
growth_factor = 2
cutback_factor = .5
[]
[]
[Postprocessors]
[ave_temp_interior] # average temperature of the cladding interior and all pellet exteriors
type = LayeredSideAverageValuePostprocessor
boundary = 9
variable = temperature
execute_on = 'initial linear'
fuel_pin_geometry = pin_geometry
[]
[clad_inner_vol] # volume inside of cladding
type = LayeredInternalVolumePostprocessor
boundary = 7
component = 0
fuel_pin_geometry = pin_geometry
out_of_plane_strain = strain_yy
execute_on = 'initial linear'
[]
[pellet_volume] # fuel pellet total volume
type = LayeredInternalVolumePostprocessor
boundary = 8
component = 0
fuel_pin_geometry = pin_geometry
out_of_plane_strain = strain_yy
execute_on = 'initial linear'
[]
[avg_clad_temp] # average temperature of cladding interior
type = LayeredSideAverageValuePostprocessor
boundary = 7
variable = temperature
fuel_pin_geometry = pin_geometry
execute_on = 'initial linear'
[]
[fis_gas_produced] # fission gas produced (moles)
type = LayeredElementIntegralFisGasGeneratedSifgrsPostprocessor
block = fuel
fuel_pin_geometry = pin_geometry
[]
[fis_gas_released] # fission gas released to plenum (moles)
type = LayeredElementIntegralFisGasReleasedSifgrsPostprocessor
block = fuel
fuel_pin_geometry = pin_geometry
[]
[fis_gas_grain]
type = LayeredElementIntegralFisGasGrainSifgrsPostprocessor
block = fuel
outputs = exodus
fuel_pin_geometry = pin_geometry
[]
[fis_gas_boundary]
type = LayeredElementIntegralFisGasBoundarySifgrsPostprocessor
block = fuel
outputs = exodus
fuel_pin_geometry = pin_geometry
[]
[fission_gas_release]
type = FGRPercent
fission_gas_released = fis_gas_released
fission_gas_generated = fis_gas_produced
[]
[gas_volume]
type = LayeredInternalVolumePostprocessor
boundary = 9
execute_on = 'initial linear'
component = 0
out_of_plane_strain = strain_yy
fuel_pin_geometry = pin_geometry
[]
[flux_from_clad] # area integrated heat flux from the cladding
type = LayeredSideFluxIntegralPostprocessor
variable = temperature
boundary = 5
diffusivity = thermal_conductivity
fuel_pin_geometry = pin_geometry
[]
[flux_from_fuel] # area integrated heat flux from the fuel
type = LayeredSideFluxIntegralPostprocessor
variable = temperature
boundary = 10
diffusivity = thermal_conductivity
fuel_pin_geometry = pin_geometry
[]
[_dt] # time step
type = TimestepSize
[]
[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 = LayeredElementIntegralPowerPostprocessor
variable = temperature
burnup_function = burnup
block = fuel
fuel_pin_geometry = pin_geometry
[]
[rod_input_power]
type = FunctionValuePostprocessor
function = power_history
scale_factor = 0.1186 # rod height
[]
[solid_swelling]
type = ElementAverageValue
variable = solid_swell
block = fuel
[]
[gaseous_swelling]
type = ElementAverageValue
variable = gaseous_swelling
block = fuel
[]
[densification]
type = ElementAverageValue
variable = densification
block = fuel
[]
[volumetric_swelling]
type = ElementAverageValue
variable = volumetric_swelling_strain
block = fuel
[]
[]
[Outputs]
perf_graph = true
exodus = true
csv = true
color = false
[]
(test/tests/thermalCompositeSiC/thermal_koyanagi_stone_combined.i)
# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1800 K
# over 15.0 s.
#
#
# Unirradiated thermal conductivity is computed using the thermal diffusivity data from:
#
# T. Koyanagi, Y. Katoh, G. Singh, M. Snead, "SiC/SiC Cladding Materials
# Properties Handbook," Technical Report ORNL/TM-2017/385, 2017.
#
# Irradiation resistance is computed according to:
#
# J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, C. A.
# Back, "Stress analysis and probablistic assessment of multi-layer SiC-based
# accident tolerant nuclear fuel cladding," Journal of Nuclear Materials, 466,
# 2015, p. 682-697.
#
# 1 / k_irradiated = R_0 + R_irr
#
# where R_0 = 1 / k_unirradiated and R_irr = 15.11 * S where S is the swelling
# strain.
#
# The thermal conductivity and specific heat computed by BISON for the
# temperatures shown below and compared with analytical solutions.
# The results are as follows:
#
# SiC/SiC k SiC/SiC k SiC/SiC Cp SiC/SiC Cp
# Temp(K) BISON(W/m-K) analytical(W/m-K) BISON(J/kg-K) analytical(J/kg-K)
# 300 8.20058 8.20058 676.721 676.721
# 600 1.67342 1.66730 1034.698 1034.697
# 900 1.62824 1.63774 1161.491 1161.491
# 1200 1.59471 1.60440 1241.972 1241.972
# 1800 1.44797 1.45651 1337.951 1337.9509
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
[cube]
type = GeneratedMeshGenerator
dim = 3
[]
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[temperature]
order = FIRST
family = LAGRANGE
initial_condition = 300.0
[]
[]
[AuxVariables]
[fast_neutron_fluence]
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[specific_heat]
order = CONSTANT
family = MONOMIAL
[]
[thermal_conductivity]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[block]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temperature
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temperature
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence
block = 0
execute_on = 'timestep_begin'
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[thermal_conductivity]
type = MaterialRealAux
variable = thermal_conductivity
property = thermal_conductivity
block = 0
execute_on = 'initial linear'
[]
[specific_heat]
type = MaterialRealAux
variable = specific_heat
property = specific_heat
block = 0
execute_on = 'initial linear'
[]
[]
[Functions]
[fast_neutron_fluence]
type = PiecewiseLinear
x = '0 5.0'
y = '0 5e25'
[]
[temp_ramp]
type = PiecewiseLinear
x = '0.0 5.0'
y = '300 1800'
[]
[]
[BCs]
[temperature_all]
type = FunctionDirichletBC
boundary = 'left right front back bottom top' # All cube faces
variable = temperature
function = temp_ramp
[]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[frontY]
type = DirichletBC
variable = disp_y
boundary = 'front'
value = 0.0
[]
[bottomZ]
type = DirichletBC
variable = disp_z
boundary = 'bottom'
value = 0.0
[]
[]
[Materials]
[thermalCompositeSiC]
type = CompositeSiCThermal
temperature = temperature
thermal_conductivity_model = COMBINED
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 2700.0
[]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temperature
swelling_model = KATOH
number_of_substeps = 100
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
block = 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-5
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
num_steps = 10
dt = 0.5
[]
[Postprocessors]
[avg_temperature]
type = SideAverageValue
boundary = 'left right top bottom front back'
execute_on = 'initial timestep_end'
variable = temperature
[]
[avg_th_cond]
type = ElementAverageValue
block = 0
variable = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[avg_specific_heat]
type = ElementAverageValue
block = 0
variable = specific_heat
execute_on = 'initial timestep_end'
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]
(test/tests/thermalCompositeSiC/thermal_stone_unirradiated.i)
# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500 K
# over 12.0 s.
#
# Specific heat capacity is computed using the correlation from:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# and thermal conductivity is computed using the thermal conductivity correlation
# from:
#
# J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, C. A.
# Back, "Stress analysis and probablistic assessment of multi-layer SiC-based
# accident tolerant nuclear fuel cladding," Journal of Nuclear Materials, 466,
# 2015, p. 682-697.
#
# The specific heat correlation is tested elsewhere. The results of the BISON
# simulation for unirradiatedthermal conductivity is compared to the analytical
# solution below.
#
# SiC/SiC k SiC/SiC k
# Temp(K) BISON(W/m-K) analytical(W/m-K)
# 300 18.21599 18.21599
# 600 18.62984 18.62984
# 900 16.13219 16.13219
# 1200 14.31944 14.31944
# 1500 13.46375 13.46375
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 3
[]
[]
[Variables]
[temperature]
order = FIRST
family = LAGRANGE
initial_condition = 300.0
[]
[]
[AuxVariables]
[specific_heat]
order = CONSTANT
family = MONOMIAL
[]
[thermal_conductivity]
order = CONSTANT
family = MONOMIAL
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temperature
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temperature
[]
[]
[AuxKernels]
[thermal_conductivity]
type = MaterialRealAux
variable = thermal_conductivity
property = thermal_conductivity
block = 0
execute_on = 'initial linear'
[]
[specific_heat]
type = MaterialRealAux
variable = specific_heat
property = specific_heat
block = 0
execute_on = 'initial linear'
[]
[]
[Functions]
[temp_ramp]
type = PiecewiseLinear
x = '0.0 12.0'
y = '300 1500'
[]
[]
[BCs]
[temperature_all]
type = FunctionDirichletBC
boundary = 'left right front back bottom top' # All cube faces
variable = temperature
function = temp_ramp
[]
[]
[Materials]
[thermalCompositeSiC]
type = CompositeSiCThermal
temperature = temperature
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = ParsedMaterial
block = 0
property_name = density
expression = 2700.0
[]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temperature
fast_neutron_fluence = 0.0
eigenstrain_name = swelling_strain
[]
[]
[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-6
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
num_steps = 12
dt = 1.0
[]
[Postprocessors]
[avg_temperature]
type = SideAverageValue
boundary = 'left right top bottom front back'
execute_on = 'initial timestep_end'
variable = temperature
[]
[avg_th_cond]
type = ElementAverageValue
block = 0
variable = thermal_conductivity
execute_on = 'initial timestep_end'
[]
[avg_specific_heat]
type = ElementAverageValue
block = 0
variable = specific_heat
execute_on = 'initial timestep_end'
[]
[]
[Outputs]
exodus = true
[]
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh.i)
# This test tests the Katoh swelling model used for composite and monolithic
# silicon carbide. The model is obtained from:
#
# Y. Katoh, K. Ozawa, C. Shih, T. Nozawa, R. J. Shinavski, A. Hasegawa, and
# L. L. Snead. Continuous SiC fiber, CVI SiC matrix composites for nuclear
# applications: properties and irradiation effects. Journal of Nuclear Materials,
# 448, pp. 448-476, 2014.
#
# The swelling rate is a function of the fast neutron fluence:
#
# dS/dgamma = k(T) * gamma^(-1/3) * exp(-gamma/gamma_sc(T))
#
# where gamma is the fluence in dpa, S is the swelling strain, and k(T) and
# gamma_sc(T) are temperature dependent constants. The above equation cannot
# be analytically integrated and an incremental approach is used to calculate the
# swelling strain. The "analytical" solution is taken as the results of the
# numerical integral calculated using Wolfram Alpha or Mathematica. Depending
# upon the number of substeps used the BISON simulations approach this analytical
# solution.
#
# Here, a single 1x1 m 2D element is subjected to a constant temperature of
# 1073.15 K with a varying fast neutron fluence from 0 to 1e25 n/m^2 (0 to 5 dpa)
# over 1 second. Below BISON results are included for 100, 1000 and 10000 substeps to
# illustrate that the results approach that of the analytical solution. This test
# is for the 100 substep case. Tests for the 1000 and 10000 substeps at high
# temperature as well as low and mid range temperatures of 473.15 K
# and 773.15 K (with 100 substeps) are included using cli_args in the tests file.
#
# Fluence BISON Swelling Strain (-) Analytical
# (dpa) 100 substeps 1000 substeps 10000 substeps Swelling
# 1.0 5.91955e-3 5.83353e-3 5.83899e-3 5.84045e-3
# 2.0 7.20081e-3 7.22617e-3 7.23163e-3 7.23270e-3
# 3.0 7.70534e-3 7.73070e-3 7.73617e-3 7.73713e-3
# 4.0 7.90061e-3 7.92597e-3 7.93143e-3 7.93235e-3
# 5.0 7.97852e-3 8.00388e-3 8.00934e-3 8.01025e-3
#
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = GeneratedMeshGenerator
dim = 2
[]
[]
[Variables]
[disp_x]
order = FIRST
family = LAGRANGE
[]
[disp_y]
order = FIRST
family = LAGRANGE
[]
[temp]
order = FIRST
family = LAGRANGE
initial_condition = 1073.15
[]
[]
[AuxVariables]
[fast_neutron_fluence]
order = FIRST
family = LAGRANGE
[]
[swell]
order = CONSTANT
family = MONOMIAL
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[square]
block = 0
add_variables = false
strain = FINITE
eigenstrain_names = 'swelling_strain'
[]
[]
[Kernels]
[heat]
type = HeatConduction
variable = temp
[]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
[]
[]
[Functions]
[fast_neutron_fluence_func]
type = PiecewiseLinear
x = '0 1.0'
y = '0 5e25'
[]
[]
[AuxKernels]
[fast_neutron_fluence]
type = FunctionAux
variable = fast_neutron_fluence
function = fast_neutron_fluence_func
execute_on = timestep_begin
[]
[total_swelling]
type = MaterialRealAux
variable = swell
property = swelling
block = 0
[]
[]
[BCs]
[leftX]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[bottomY]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[all_temp]
type = DirichletBC
variable = temp
boundary = 'bottom left right top'
value = 1073.15
[]
[]
[Materials]
[swelling]
type = CompositeSiCVolumetricSwellingEigenstrain
block = 0
temperature = temp
swelling_model = KATOH
number_of_substeps = 100
fast_neutron_fluence = fast_neutron_fluence
eigenstrain_name = swelling_strain
[]
[elasticity_tensor] # isotropic elasticity tensor
type = ComputeIsotropicElasticityTensor
block = 0
youngs_modulus = 3e11
poissons_ratio = 0.17
[]
[elastic_stress]
type = ComputeFiniteStrainElasticStress
block = 0
[]
[thermalCompositeSiC]
type = CompositeSiCThermal
temperature = temp
thermal_conductivity_model = STONE # This is the default
[]
[density]
type = StrainAdjustedDensity
block = 0
strain_free_density = 3100.0
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options = '-snes_ksp_ew '
petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type'
petsc_options_value = '70 hypre boomeramg'
l_max_its = 60
nl_rel_tol = 1e-5
nl_abs_tol = 1e-10
l_tol = 1e-5
start_time = 0.0
dt = 0.01
end_time = 1.0
[]
[Postprocessors]
[avgFuelTemp]
type = ElementAverageValue
variable = temp
execute_on = 'initial timestep_end'
block = 0
[]
[swelling_out]
type = ElementAverageValue
variable = swell
execute_on = 'initial timestep_end'
block = 0
[]
[average_fluence]
type = ElementAverageValue
variable = fast_neutron_fluence
execute_on = 'initial timestep_end'
block = 0
[]
[]
[Outputs]
exodus = true
[]