- ATemperature-independent coefficient of the intercept term in the concentration-pressure proportionality
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Temperature-independent coefficient of the intercept term in the concentration-pressure proportionality
- BTemperature-dependent coefficient of the intercept term in the concentration-pressure proportionality
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Temperature-dependent coefficient of the intercept term in the concentration-pressure proportionality
- DTemperature-independent coefficient of the order in the concentration-pressure proportionality
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Temperature-independent coefficient of the order in the concentration-pressure proportionality
- ETemperature-dependent coefficient of the order in the concentration-pressure proportionality
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Temperature-dependent coefficient of the order in the concentration-pressure proportionality
- concentrationThe coupled concentration variable
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The coupled concentration variable
- d1Temperature-independent coefficient of the transition concentration
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Temperature-independent coefficient of the transition concentration
- d2Temperature-dependent coefficient of the transition concentration
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Temperature-dependent coefficient of the transition concentration
- densityMass density of the material
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Mass density of the material
- temperatureThe coupled temperature variable
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The coupled temperature variable
Sorption Partial Pressure
Computes the sorption partial pressure of a solute in a material.
Description
The SorptionPartialPressure class is used to calculate the partial pressure due to sorption of a solute at a material surface. The relationship between the partial pressure and the concentration is linear at low concentrations and transitions to power-law behavior at higher concentrations. The two regimes are referred to as Hernian and Freundlich, respectively. The partial pressure is given by (1) where is the solute concentration, and is the temperature. Here also, , , , , and are sorption constants specific to a particular solute and material.
Sorption constants in the literature may only be compatible with concentrations expressed in certain units, and those units may differ from those used elsewhere in the simulation, e.g., versus . The model accepts material mass densities and concentration unit conversion factors to accomodate data from different sources. For example, densities in and unit conversion factors of can be supplied to accomodate sets of sorption constants fit for concentrations in . The model assumes that the units associated with the partial pressures are equal and that all sorption constants are compatible with absolute temperature.
An optional base name can be prepended to the calculated partial pressure material property, allowing calculation of multiple partial pressures on the same block, i.e., for multiple sorping species. The names of the partial pressure material properties can be passed to SorptionIsothermGapInterface or MassSorptionConstraint to model sorption conditions.
Example Input Syntax
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[partial_pressure]
type = SorptionPartialPressure<<<{"description": "Computes the sorption partial pressure of a solute in a material.", "href": "SorptionPartialPressure.html"}>>>
A<<<{"description": "Temperature-independent coefficient of the intercept term in the concentration-pressure proportionality"}>>> = 19.3
B<<<{"description": "Temperature-dependent coefficient of the intercept term in the concentration-pressure proportionality"}>>> = -47300
D<<<{"description": "Temperature-independent coefficient of the order in the concentration-pressure proportionality"}>>> = 1.51
E<<<{"description": "Temperature-dependent coefficient of the order in the concentration-pressure proportionality"}>>> = 4340
d1<<<{"description": "Temperature-independent coefficient of the transition concentration"}>>> = 3.4
d2<<<{"description": "Temperature-dependent coefficient of the transition concentration"}>>> = 6.15e-4
unit_scale<<<{"description": "Unit conversion factor used to scale the concentration"}>>> = 1
density<<<{"description": "Mass density of the material"}>>> = 1
concentration<<<{"description": "The coupled concentration variable"}>>> = c
temperature<<<{"description": "The coupled temperature variable"}>>> = T
partial_pressure_base_name<<<{"description": "Optional parameter that allows the user to define multiple partial pressures on the same block, i.e. for multiple concentrations"}>>> = test
[]
[]Input Parameters
- 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.
- partial_pressure_base_nameOptional parameter that allows the user to define multiple partial pressures on the same block, i.e. for multiple concentrations
C++ Type:std::string
Controllable:No
Description:Optional parameter that allows the user to define multiple partial pressures on the same block, i.e. for multiple concentrations
- unit_scale1Unit conversion factor used to scale the concentration
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Unit conversion factor used to scale the concentration
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/sorption_constraint/test.i)
- (test/tests/sorption_interface/sorption.i)
- (test/tests/sorption_interface/desorption.i)
- (examples/TRISO/accident_simulation/triso1D_accident.i)
- (assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Sr.i)
- (test/tests/sorption_interface/test.i)
- (test/tests/sorption_interface/adsorption.i)
- (assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Ag.i)
- (examples/TRISO/accident_simulation/triso2D_accident_ad.i)
- (assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Cs.i)
- (test/tests/sorption_partial_pressure/test.i)
- (test/tests/sorption_interface/sorption_ad.i)
References
No citations exist within this document.(test/tests/sorption_partial_pressure/test.i)
# This test is to verify the implementation of SorptionPartialPressure material.
# It compares calculated partial pressure parameters to analytical solutions.
[Problem]
solve = false
[]
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[Functions]
[partial_pressure_avg_exact]
type = ParsedFunction
symbol_names = 'c T A B D E d1 d2'
symbol_values = 'c_avg T_avg 19.3 -47300 1.51 4340 3.4 6.15e-4'
expression = 'exp(A + B / T) * c ^ (D + E / T) + exp(A + B / T + (D - 1 + E / T) * (d1 - d2 * T)) * c'
[]
[partial_pressure_dc_avg_exact]
type = ParsedFunction
symbol_names = 'c T A B D E d1 d2'
symbol_values = 'c_avg T_avg 19.3 -47300 1.51 4340 3.4 6.15e-4'
expression = '(D + E / T) * exp(A + B / T) * c ^ (D + E / T - 1) + exp(A + B / T + (D - 1 + E / T) * (d1 - d2 * T))'
[]
[]
[AuxVariables]
[c]
order = FIRST
family = LAGRANGE
initial_condition = 1
[]
[T]
order = FIRST
family = LAGRANGE
initial_condition = 1000
[]
[]
[Materials]
[partial_pressure]
type = SorptionPartialPressure
A = 19.3
B = -47300
D = 1.51
E = 4340
d1 = 3.4
d2 = 6.15e-4
unit_scale = 1
density = 1
concentration = c
temperature = T
partial_pressure_base_name = test
[]
[]
[Executioner]
type = Steady
solve_type = NEWTON
[]
[Postprocessors]
[c_avg]
type = ElementAverageValue
variable = c
[]
[T_avg]
type = ElementAverageValue
variable = T
[]
[partial_pressure_avg]
type = ElementAverageMaterialProperty
mat_prop = test_partial_pressure
[]
[partial_pressure_dc_avg]
type = ElementAverageMaterialProperty
mat_prop = test_partial_pressure_dc
[]
[partial_pressure_avg_exact]
type = FunctionValuePostprocessor
function = partial_pressure_avg_exact
[]
[partial_pressure_dc_avg_exact]
type = FunctionValuePostprocessor
function = partial_pressure_dc_avg_exact
[]
[]
[Outputs]
csv = true
[]
(test/tests/sorption_constraint/test.i)
# This test is to verify the implementation of MassSorptionConstraint constraint.
# It contains two 2D blocks separated by an unmeshed flat gap.
# 1D mortar subdomains are set up to enforce sorption and preserve flux between the blocks.
# Checks are performed to verify concentration conservation, sorption behavior, and flux preservation.
[GlobalParams]
order = FIRST
family = LAGRANGE
[]
[Mesh]
coord_type = XYZ
patch_update_strategy = auto
[inner] # create inner block
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = 0.5
ymin = 0
ymax = 1
nx = 10
ny = 10
boundary_name_prefix = inner
elem_type = QUAD4
[]
[inner_id]
type = SubdomainIDGenerator
input = inner
subdomain_id = 1
[]
[outer] # create outer block
type = GeneratedMeshGenerator
dim = 2
xmin = 0.51
xmax = 1.01
ymin = 0
ymax = 1
nx = 10
ny = 10
boundary_name_prefix = outer
boundary_id_offset = 10
elem_type = QUAD4
[]
[outer_id]
type = SubdomainIDGenerator
input = outer
subdomain_id = 2
[]
[combined]
type = MeshCollectionGenerator
inputs = 'inner_id outer_id'
[]
[block_rename]
type = RenameBlockGenerator
input = combined
old_block = '1 2'
new_block = 'inner outer'
[]
[primary] # create left mortar subdomain between the inner and outer blocks
type = LowerDBlockFromSidesetGenerator
input = block_rename
sidesets = 1
new_block_id = 100
new_block_name = primary
[]
[secondary] # create right mortar subdomain between the inner and outer blocks
type = LowerDBlockFromSidesetGenerator
input = primary
sidesets = 13
new_block_id = 200
new_block_name = secondary
[]
[]
[Variables]
[u] # concentration
block = 'inner outer'
initial_condition = 2
[]
[temperature]
block = 'inner outer'
initial_condition = 550
[]
[lm_u] # constraint on the concentration in the gap
block = secondary
[]
[lm_dx] # constraint on the flux in the x direction in the gap
block = secondary
[]
[lm_dy] # constraint on the flux in the y direction in the gap
block = secondary
[]
[]
[BCs]
[u_left]
type = ADDirichletBC
value = 1
variable = u
boundary = inner_left
[]
[u_right]
type = ADDirichletBC
value = 1
variable = u
boundary = outer_right
[]
[temperature_left]
type = ADDirichletBC
value = 1100
variable = temperature
boundary = inner_left
[]
[temperature_primary]
type = ADDirichletBC
value = 1000
variable = temperature
boundary = inner_right
[]
[temperature_secondary]
type = ADDirichletBC
value = 900
variable = temperature
boundary = outer_left
[]
[temperature_right]
type = ADDirichletBC
value = 0
variable = temperature
boundary = outer_right
[]
[]
[Kernels]
[u]
type = ADMatDiffusion
variable = u
diffusivity = diffusivity
block = 'inner outer'
[]
[temperature]
type = ADHeatConduction
variable = temperature
thermal_conductivity = thermal_conductivity
block = 'inner outer'
[]
[]
[Constraints]
[u]
type = MassSorptionConstraint
variable = lm_u
primary_variable = u
primary_boundary = inner_right
primary_subdomain = primary
secondary_variable = u
secondary_boundary = outer_left
secondary_subdomain = secondary
partial_pressure_name = partial_pressure
epsilon = 1e-9
correct_edge_dropping = true
[]
[dx]
type = MassFluxConstraint
variable = lm_dx
primary_variable = u
diffusivity_primary = diffusivity
primary_boundary = inner_right
primary_subdomain = primary
secondary_variable = u
diffusivity_secondary = diffusivity
secondary_boundary = outer_left
secondary_subdomain = secondary
component = 0
epsilon = 1e-5
correct_edge_dropping = true
[]
[dy]
type = MassFluxConstraint
variable = lm_dy
primary_variable = u
diffusivity_primary = diffusivity
primary_boundary = inner_right
primary_subdomain = primary
secondary_variable = u
diffusivity_secondary = diffusivity
secondary_boundary = outer_left
secondary_subdomain = secondary
component = 1
epsilon = 1e-5
correct_edge_dropping = true
[]
[]
[Materials]
[properties_inner]
type = ADGenericConstantMaterial
prop_names = 'thermal_conductivity diffusivity density'
prop_values = '1 1 1'
block = inner
outputs = all
[]
[partial_pressure_inner]
type = ADSorptionPartialPressure
A = 19.3
B = -47300
D = 1.51
E = 4340
d1 = 3.4
d2 = 6.15e-4
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = inner
outputs = 'all'
output_properties = partial_pressure
[]
[properties_outer]
type = ADGenericConstantMaterial
prop_names = 'thermal_conductivity diffusivity density'
prop_values = '2 2 1'
block = outer
outputs = all
[]
[partial_pressure_outer]
type = ADSorptionPartialPressure
A = 24
B = -35700
D = -1.56
E = 6120
d1 = 2.04
d2 = 1.79e-3
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = outer
outputs = 'all'
output_properties = partial_pressure
[]
[]
[Functions]
[u_mid_diff]
type = ParsedFunction
symbol_names = 'u_mid_inner u_mid_outer'
symbol_values = 'u_mid_inner u_mid_outer'
expression = '(abs(u_mid_outer) - abs(u_mid_inner)) / abs(u_mid_inner)'
[]
[pressure_mid_error]
type = ParsedFunction
symbol_names = 'pressure_mid_inner pressure_mid_outer'
symbol_values = 'pressure_mid_inner pressure_mid_outer'
expression = '(abs(pressure_mid_outer) - abs(pressure_mid_inner)) / abs(pressure_mid_inner)'
[]
[flux_error]
type = ParsedFunction
symbol_names = 'flux_inner flux_outer'
symbol_values = 'flux_inner flux_outer'
expression = '(abs(flux_outer) - abs(flux_inner)) / abs(flux_inner)'
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Steady
solve_type = NEWTON
automatic_scaling = true
scaling_group_variables = 'u; temperature; lm_u; lm_dx; lm_dy'
petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
petsc_options_value = 'lu NONZERO 1e-15' # provide nonzero on-diagonals for lagrange multiplier
line_search = none
nl_rel_tol = 1e-15
nl_abs_tol = 1e-9
l_tol = 1e-3
l_max_its = 50
nl_max_its = 50
[]
[Postprocessors]
[u_mid_inner] # verify sorption behavior
type = PointValue
variable = u
point = '0.5 0.5 0'
outputs = csv
[]
[u_mid_outer]
type = PointValue
variable = u
point = '0.51 0.5 0'
outputs = csv
[]
[u_mid_diff]
type = FunctionValuePostprocessor
function = u_mid_diff
outputs = console
[]
[temperature_mid_inner]
type = PointValue
variable = temperature
point = '0.5 0.5 0'
outputs = csv
[]
[temperature_mid_outer]
type = PointValue
variable = temperature
point = '0.51 0.5 0'
outputs = csv
[]
[pressure_mid_inner]
type = PointValue
variable = partial_pressure
point = '0.51 0.5 0'
outputs = csv
[]
[pressure_mid_outer]
type = PointValue
variable = partial_pressure
point = '0.51 0.5 0'
outputs = csv
[]
[pressure_mid_error]
type = FunctionValuePostprocessor
function = pressure_mid_error
outputs = console
[]
[flux_inner] # verify flux preservation
type = ADSideDiffusiveFluxIntegral
variable = u
boundary = inner_right
diffusivity = diffusivity
outputs = csv
[]
[flux_outer]
type = ADSideDiffusiveFluxIntegral
variable = u
boundary = outer_left
diffusivity = diffusivity
outputs = csv
[]
[flux_error]
type = FunctionValuePostprocessor
function = flux_error
outputs = console
[]
[]
[Outputs]
exodus = true
csv = true
[]
[Debug]
show_var_residual_norms = true
[]
(test/tests/sorption_interface/sorption.i)
# This test is verifies the implementation of SorptionIsothermGapInterface with
# interface_configuration = solid_gas. It contains two solids separated by a gas.
# One solid acts as a source of u. u desorbs into the gap and then onto the
# other solid, which acts as a sink. Temperature increases over time to verify
# sorption behavior and species conservation. For SorptionIsothermGapInterface,
# the solid is always the element and the gas is always the neighbor.
# _________________________
# | source | gap | sink |
# --------------------------
[GlobalParams]
order = SECOND
[]
[Mesh]
coord_type = RZ
[gen]
type = GeneratedMeshGenerator
nx = 150
xmax = 1.5
dim = 1
elem_type = EDGE3
[]
[source_gap]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 0'
top_right = '0.75 0 0'
input = gen
[]
[gap_sink]
type = SubdomainBoundingBoxGenerator
block_id = 2
bottom_left = '0.75 0 0'
top_right = '1.5 0 0'
input = source_gap
[]
[source]
type = SubdomainBoundingBoxGenerator
block_id = 3
bottom_left = '0 0 0'
top_right = '0.5 0 0'
input = gap_sink
[]
[sink]
type = SubdomainBoundingBoxGenerator
block_id = 4
bottom_left = '1 0 0'
top_right = '1.5 0 0'
input = source
[]
[break]
type = BreakMeshByBlockGenerator
block_pairs = '1 3; 2 4'
split_interface = true
add_interface_on_two_sides = true
input = sink
[]
[]
[Variables]
[u]
[]
[]
[AuxVariables]
[temperature]
[]
[]
[Kernels]
[u_dt]
type = TimeDerivative
variable = u
[]
[u]
type = MatDiffusion
variable = u
diffusivity = diffusivity
[]
[]
[ICs]
[u_source]
type = FunctionIC
variable = u
function = 'if(x < 0.75, 1, 0)'
block = 3
[]
[]
[AuxKernels]
[temperature]
type = FunctionAux
variable = temperature
function = '(4000 - 2000 * x) / 2.0 * (t + 0.1)'
execute_on = 'INITIAL TIMESTEP_BEGIN LINEAR NONLINEAR TIMESTEP_END'
[]
[]
[InterfaceKernels]
[u_solid_gas]
type = SorptionIsothermGapInterface
variable = u
neighbor_var = u
interface_configuration = solid_gas
temperature = temperature
partial_pressure_name = partial_pressure
sorption_penalty = 1e10
boundary = 'Block3_Block1 Block4_Block2'
[]
[]
[Materials]
[properties_solid]
type = GenericConstantMaterial
prop_names = 'density diffusivity'
prop_values = '1 0.1'
block = '3 4'
[]
[properties_gas]
type = GenericConstantMaterial
prop_names = diffusivity
prop_values = 1
block = '1 2'
[]
[partial_pressure_source]
type = SorptionPartialPressure
A = 19.3
B = -47300
D = 1.51
E = 4340
d1 = 3.4
d2 = 6.15e-4
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = 3
outputs = 'all'
output_properties = partial_pressure
[]
[partial_pressure_sink]
type = SorptionPartialPressure
A = 24
B = -35700
D = -1.56
E = 6120
d1 = 2.04
d2 = 1.79e-3
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = 4
outputs = 'all'
output_properties = partial_pressure
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
compute_scaling_once = false
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = none
nl_rel_tol = 1e-15
nl_abs_tol = 1e-15
l_tol = 1e-3
l_max_its = 50
nl_max_its = 50
dt = 0.1
end_time = 2
[]
[Functions]
[source_surface_u_exact]
type = ParsedFunction
symbol_names = 'P T R'
symbol_values = 'source_surface_pressure source_surface_temperature 8.31446261815324'
expression = '1.0 / R / T * P'
[]
[sink_surface_u_exact]
type = ParsedFunction
symbol_names = 'P T R'
symbol_values = 'sink_surface_pressure sink_surface_temperature 8.31446261815324'
expression = '1.0 / R / T * P'
[]
[]
[Postprocessors]
[source_surface_pressure]
type = SideAverageMaterialProperty
property = partial_pressure
boundary = Block3_Block1
[]
[source_surface_temperature]
type = SideAverageValue
variable = temperature
boundary = Block3_Block1
[]
[source_surface_u]
type = SideAverageValue
variable = u
boundary = Block1_Block3
[]
[source_surface_u_exact]
type = FunctionValuePostprocessor
function = source_surface_u_exact
[]
[sink_surface_pressure]
type = SideAverageMaterialProperty
property = partial_pressure
boundary = Block4_Block2
[]
[sink_surface_temperature]
type = SideAverageValue
variable = temperature
boundary = Block4_Block2
[]
[sink_surface_u]
type = SideAverageValue
variable = u
boundary = Block2_Block4
[]
[sink_surface_u_exact]
type = FunctionValuePostprocessor
function = sink_surface_u_exact
[]
[u_source_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = 3
[]
[u_gap_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = '1 2'
[]
[u_sink_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = 4
[]
[u_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = '1 2 3 4'
[]
[]
[Outputs]
exodus = true
csv = true
[]
[Debug]
show_var_residual_norms = true
[]
(test/tests/sorption_interface/desorption.i)
# This test is verifies the implementation of SorptionIsothermGapInterface with
# interface_configuration = solid_gas. It contains a solid having a single
# interface with a gas. The solid acts as a source of u. u desorbs into the gap.
# Temperature increases over time to verify desorption behavior and species
# conservation. For SorptionIsothermGapInterface, the solid is always the
# element and the gas is always the neighbor.
# _________________
# | source | gap |
# -----------------
[GlobalParams]
order = SECOND
[]
[Mesh]
[gen]
type = GeneratedMeshGenerator
nx = 10
ny = 2
ymax = 0.2
dim = 2
elem_type = QUAD8
[]
[source]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 0'
top_right = '0.5 0.2 0'
input = gen
[]
[gap]
type = SubdomainBoundingBoxGenerator
block_id = 2
bottom_left = '0.5 0 0'
top_right = '1 0.2 0'
input = source
[]
[break]
type = BreakMeshByBlockGenerator
block_pairs = '1 2'
split_interface = true
add_interface_on_two_sides = true
input = gap
[]
[]
[Variables]
[u]
[]
[]
[AuxVariables]
[temperature]
[]
[]
[Kernels]
[u_dt]
type = TimeDerivative
variable = u
[]
[u]
type = MatDiffusion
variable = u
diffusivity = diffusivity
[]
[]
[ICs]
[u_source]
type = ConstantIC
variable = u
value = 1
block = 1
[]
[]
[AuxKernels]
[temperature]
type = FunctionAux
variable = temperature
function = '(4000 - 1000 * x) / 2.0 * (t + 0.1)'
execute_on = 'INITIAL TIMESTEP_BEGIN LINEAR NONLINEAR TIMESTEP_END'
[]
[]
[InterfaceKernels]
[u_source_gap]
type = SorptionIsothermGapInterface
variable = u
neighbor_var = u
interface_configuration = solid_gas
temperature = temperature
partial_pressure_name = partial_pressure
sorption_penalty = 1e10
boundary = Block1_Block2
[]
[]
[Materials]
[properties_solid]
type = GenericConstantMaterial
prop_names = 'density diffusivity'
prop_values = '1 0.1'
block = 1
[]
[properties_gas]
type = GenericConstantMaterial
prop_names = diffusivity
prop_values = 1
block = 2
[]
[partial_pressure]
type = SorptionPartialPressure
A = 19.3
B = -47300
D = 1.51
E = 4340
d1 = 3.4
d2 = 6.15e-4
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = 1
outputs = 'all'
output_properties = partial_pressure
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
compute_scaling_once = false
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = none
nl_rel_tol = 1e-15
nl_abs_tol = 1e-15
l_tol = 1e-3
l_max_its = 50
nl_max_its = 50
dt = 0.1
end_time = 2
[]
[Functions]
[surface_u_exact]
type = ParsedFunction
symbol_names = 'P T R'
symbol_values = 'surface_pressure surface_temperature 8.31446261815324'
expression = '1.0 / R / T * P'
[]
[]
[Postprocessors]
[surface_pressure]
type = SideAverageMaterialProperty
property = partial_pressure
boundary = Block1_Block2
[]
[surface_temperature]
type = SideAverageValue
variable = temperature
boundary = Block1_Block2
[]
[surface_u]
type = SideAverageValue
variable = u
boundary = Block2_Block1
[]
[surface_u_exact]
type = FunctionValuePostprocessor
function = surface_u_exact
[]
[u_source_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = 1
[]
[u_gap_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = 2
[]
[u_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = '1 2'
[]
[]
[Outputs]
exodus = true
csv = true
[]
[Debug]
show_var_residual_norms = true
[]
(examples/TRISO/accident_simulation/triso1D_accident.i)
# This example is 1D spherical analysis of a TRISO fuel particle. Fully coupled
# heat transfer and solid mechanics, plus diffusion of the fission product
# species cesium (Cs) are simulated. The mesh includes contact surfaces
# between the buffer and IPyC layers to facilitate a gap opening between
# these layers. These surfaces are initially in mechanical contact but
# are assumed to have no strength in tension. A coarse mesh is used to
# provide a short run time.
# The calculation simulates fuel-life in three steps. The first step is an
# irradiation period, where constant power and a fixed particle surface
# temperature (1500 K) are assumed over a lifetime of 76 Ms (2.4 yrs).
# For the second step, fuel removal and storage are simulated by setting
# the reactor power and Cs source terms to zero, reducing the particle
# surface temperature to ambient (300 K), and then holding it
# for 100 days. A third and final step simulates accident
# behavior by increasing the particle surface temperature from ambient
# to 2073 K over 2 hrs, and then holding it at this elevated temperature
# for an additional 200 hrs. At the particle outer boundary, the Cs
# concentration is held at zero and the pressure at ambient during the
# entire simulation. The particle is assumed to be stress-free at an
# initial temperature of 1500 K.
#
# Details about this simulation are given in Section 4 of the following
# article: J. D. Hales, R. L. Williamson, S. R. Novascone, D. M. Perez,
# B. W. Spencer and G. Pastore, "Multidimensional multiphysics simulation
# of TRISO particle fuel", Journal of Nuclear Materials, Vol. 443, p. 531,
# 2013.
# This is a version using an interface kernel to model gap mass transfer.
# Sorption constants are given in Table 1 of the following article: A.
# Londono-Hurtado, I. Szlufarska, R. Bratton and D. Morgan, "A review of
# fission product sorption in carbon structures", Journal of Nuclear
# Materials, Vol. 426, p. 254, 2012.
initial_fuel_density = 11000
[GlobalParams]
order = SECOND
family = LAGRANGE
displacements = disp_x
flux_conversion_factor = 0.85
[]
[Mesh]
coord_type = RSPHERICAL
[gen] # exclude gap to establish buffer-IPyC neighbor relationships for the sorption interface kernel
type = TRISO1DMeshGenerator
elem_type = EDGE3
coordinates = '0 2.125e-4 3.125e-4 3.525e-4 3.875e-4 4.275e-4'
mesh_density = '10 5 2 2 2'
block_names = 'fuel buffer IPyC SiC OPyC'
[]
[break] # create gap between buffer and IPyC to model mechanical and thermal contact
type = BreakMeshByBlockGenerator
input = gen
block_pairs = 'buffer IPyC'
split_interface = true
add_interface_on_two_sides = true
[]
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
[]
[Variables]
[disp_x]
[]
[temp]
initial_condition = 1500.0
[]
[conc]
initial_condition = 0.0
[]
[]
[AuxVariables]
[fission_rate]
block = fuel
order = CONSTANT
family = MONOMIAL
[]
[fluence]
order = CONSTANT
family = MONOMIAL
[]
[burnup]
block = fuel
order = CONSTANT
family = MONOMIAL
[]
[gap_cond]
order = CONSTANT
family = MONOMIAL
[]
[creep_xx]
order = CONSTANT
family = MONOMIAL
[]
[creep_yy]
order = CONSTANT
family = MONOMIAL
[]
[creep_zz]
order = CONSTANT
family = MONOMIAL
[]
[gap_HTC]
order = CONSTANT
family = MONOMIAL
[]
[gap_distance]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear
x = '0 76e6 76.001e6'
y = '1 1 0'
[]
[temp_bc]
type = PiecewiseLinear
x = '0 76e6 76.001e6 84.641e6 84.6482e6'
y = '1500 1500 300 300 2073'
[]
[k_function]
type = PiecewiseLinear
x = '0 200e6'
y = '4e-37 4e-37'
[]
[d1_function]
type = ParsedFunction
expression = 'exp(t/4.5e25)'
[]
[integral_flux_error]
type = ParsedFunction
symbol_names = 'buffer_integral_flux IPyC_integral_flux'
symbol_values = 'buffer_integral_flux IPyC_integral_flux'
expression = 'IPyC_integral_flux + buffer_integral_flux'
[]
[partial_pressure_error]
type = ParsedFunction
symbol_names = 'buffer_partial_pressure IPyC_partial_pressure'
symbol_values = 'buffer_partial_pressure IPyC_partial_pressure'
expression = 'IPyC_partial_pressure - buffer_partial_pressure'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
generate_output = 'stress_xx stress_yy stress_zz stress_xy stress_yz stress_zx hydrostatic_stress'
strain = FINITE
incremental = true
add_variables = false
[default]
block = 'fuel buffer IPyC OPyC'
eigenstrain_names = 'thermal_strain swelling_strain'
extra_vector_tags = 'ref'
[]
[SiC]
block = 'SiC'
eigenstrain_names = 'thermal_strain'
extra_vector_tags = 'ref'
[]
[]
[Kernels]
[heat_ie]
type = HeatConductionTimeDerivative
variable = temp
extra_vector_tags = 'ref'
[]
[heat]
type = HeatConduction
variable = temp
extra_vector_tags = 'ref'
[]
[heat_source]
type = NeutronHeatSource
variable = temp
block = fuel
energy_per_fission = 3.2e-11 # units of J/fission
fission_rate = fission_rate
extra_vector_tags = 'ref'
[]
[mass_ie]
type = TimeDerivative
variable = conc
extra_vector_tags = 'ref'
[]
[mass]
type = ArrheniusDiffusion
variable = conc
extra_vector_tags = 'ref'
[]
[mass_source]
type = BodyForce
variable = conc
function = power_history
value = 1.22e-5 # units of moles/m**3-s
block = fuel
extra_vector_tags = 'ref'
[]
[mass_decay]
type = Decay
variable = conc
radioactive_decay_constant = 7.297e-10 # units:(1/sec) The constant for Cesium
extra_vector_tags = 'ref'
[]
[]
[AuxKernels]
[fission_rate]
type = FissionRateGeneral
fission_rate_formulation = GENERIC
variable = fission_rate
block = fuel
fission_rate_function = power_history
value = 3.89e19
execute_on = timestep_begin
[]
[fluence]
type = MaterialRealAux
property = fast_neutron_fluence
variable = fluence
[]
[burnup]
type = BurnupAux
variable = burnup
block = fuel
fission_rate = fission_rate
molecular_weight = 0.270 # units of kg/mole
execute_on = timestep_begin
density = ${initial_fuel_density}
[]
[creep_xx]
type = RankTwoAux
rank_two_tensor = creep_strain
variable = creep_xx
index_i = 0
index_j = 0
block = 'buffer IPyC SiC OPyC'
execute_on = timestep_end
[]
[creep_yy]
type = RankTwoAux
rank_two_tensor = creep_strain
variable = creep_yy
index_i = 1
index_j = 1
block = 'buffer IPyC SiC OPyC'
execute_on = timestep_end
[]
[creep_zz]
type = RankTwoAux
rank_two_tensor = creep_strain
variable = creep_zz
index_i = 2
index_j = 2
block = 'buffer IPyC SiC OPyC'
execute_on = timestep_end
[]
[conductance]
type = MaterialRealAux
property = gap_conductance
variable = gap_cond
boundary = buffer_IPyC
execute_on = linear
[]
[]
[Contact]
[pellet_clad_mechanical]
primary = IPyC_buffer
secondary = buffer_IPyC
penalty = 1e5
model = frictionless
formulation = penalty
[]
[]
[ThermalContact]
[thermal_contact]
type = GasGapHeatTransfer
variable = temp
primary = IPyC_buffer
secondary = buffer_IPyC
initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
gas_released = 'fis_gas_released co_production' # coupling to postprocessors which supply the fission gas addition, co addition
released_gas_types = 'Kr Xe;
CO'
released_fractions = '0.153 0.847;
1'
gap_geometry_type = SPHERE
tangential_tolerance = 1e-6
roughness_coef = 0.0
quadrature = true
[]
[]
[InterfaceKernels]
[cesium_gap]
type = SorptionIsothermGapInterface
variable = conc
neighbor_var = conc
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef
use_flux_penalty = true
flux_penalty = 1e3
boundary = buffer_IPyC
extra_vector_tags = 'ref'
[]
[]
[BCs]
# pin particle along symmetry planes
[no_disp_x]
type = DirichletBC
variable = disp_x
boundary = xzero
value = 0.0
extra_vector_tags = 'ref'
[]
# fix temperature on free surface
[freesurf_temp]
type = FunctionDirichletBC
variable = temp
boundary = exterior
function = temp_bc
extra_vector_tags = 'ref'
[]
# fix concentration on free surface
[freesurf_conc]
type = DirichletBC
variable = conc
boundary = exterior
value = 0.0
extra_vector_tags = 'ref'
[]
[PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
[plenumPressure]
boundary = 'buffer_IPyC IPyC_buffer'
initial_pressure = 0
startup_time = 1.0e4
R = 8.3145
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 = volumeGas # coupling to post processor to get gas volume
material_input = 'fis_gas_released co_production' # coupling to post processor to get fission gas added, co added
output = plenum_pressure # coupling to post processor to output plenum/gap pressure
[]
[]
[]
[Materials]
[flux]
type = FastNeutronFlux
calculate_fluence = true
factor = 5e17
[]
[fission_gas_release] # Sifgrs fission gas release mode
type = UO2Sifgrs
block = fuel
temperature = temp
fission_rate = fission_rate # coupling to fission_rate aux variable
grain_radius_const = 5.0e-6
[]
[fuel_thermal]
type = UO2Thermal
thermal_conductivity_model = FINK_LUCUTA
block = fuel
temperature = temp
burnup = burnup
initial_porosity = 0.0
[]
[fuel_swelling]
type = UO2VolumetricSwellingEigenstrain
gas_swelling_model_type = MATPRO
block = fuel
temperature = temp
burnup = burnup
eigenstrain_name = 'swelling_strain'
initial_fuel_density = ${initial_fuel_density}
[]
[fuel_stress]
type = ComputeFiniteStrainElasticStress
block = 'fuel'
[]
[fuel_elasticity]
type = ComputeIsotropicElasticityTensor
block = fuel
youngs_modulus = 2.2e11
poissons_ratio = .345
[]
[fuel_thermal_strain]
type = ComputeThermalExpansionEigenstrain
block = fuel
thermal_expansion_coeff = 10e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[fuel_den]
type = StrainAdjustedDensity
block = fuel
strain_free_density = ${initial_fuel_density} # kg/m^3
[]
[fuel_conc]
type = ArrheniusDiffusionCoef
block = fuel
d1 = 5.6e-8 # m^2/s
q1 = 209.0e+3 # J/mol
d2 = 5.2e-4 # m^2/s
q2 = 362.0e+3 # J/mol
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[buffer_eigenstrain]
type = PyCIrradiationEigenstrain
block = buffer
pyc_type = buffer
eigenstrain_name = 'swelling_strain'
[]
[buffer_thermal_strain]
type = ComputeThermalExpansionEigenstrain
block = buffer
thermal_expansion_coeff = 5.65e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[buffer_elasticity]
type = ComputeIsotropicElasticityTensor
block = buffer
youngs_modulus = 2e10
poissons_ratio = .23
[]
[buffer_stress]
type = PyCCreep
block = buffer
temperature = temp
[]
[buffer_temp]
type = HeatConductionMaterial
block = buffer
thermal_conductivity = 0.5 # J/m-s-K
specific_heat = 720.0 # J/kg-K
[]
[buffer_den]
type = StrainAdjustedDensity
strain_free_density = 1000.0 #kg/m^3
block = buffer
[]
[buffer_conc]
type = ArrheniusDiffusionCoef
block = buffer
d1 = 1.0e-12 # m^2/s
q1 = 0.0
d2 = 0.0
q2 = 0.0
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[buffer_partial_pressure]
type = SorptionPartialPressure
A = 19.33
B = -47290
D = 1.518
E = 4338
d1 = 3.397
d2 = 6.15e-4
unit_scale = 1e3 # convert from mol to mmol
density = density # convert from mmol/m^3 to mmol/kg
concentration = conc
temperature = temp
block = 'buffer IPyC'
outputs = 'all'
output_properties = partial_pressure
[]
[normal_vectors_triso]
type = NormalVectorsTRISO
block = 'IPyC OPyC buffer'
[]
[IPyC_eigenstrain]
type = PyCIrradiationEigenstrain
block = IPyC
pyc_type = dense
eigenstrain_name = 'swelling_strain'
[]
[IPyC_thermal_strain]
type = ComputeThermalExpansionEigenstrain
block = IPyC
thermal_expansion_coeff = 5.65e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[IPyC_elasticity]
type = ComputeIsotropicElasticityTensor
block = IPyC
youngs_modulus = 4.74e10
poissons_ratio = .23
[]
[IPyC_disp]
type = PyCCreep
block = 'IPyC OPyC'
temperature = temp
[]
[IPyC_temp]
type = HeatConductionMaterial
block = 'IPyC OPyC'
thermal_conductivity = 4.0
specific_heat = 720.0
[]
[IPyC_den]
type = StrainAdjustedDensity
block = 'IPyC OPyC'
strain_free_density = 1900.0
[]
[IPyC_conc]
type = ArrheniusDiffusionCoef
block = IPyC
d1 = 6.3e-8
q1 = 222.0e+3
d2 = 0.0
q2 = 0.0
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[SiC_thermal_strain]
type = ComputeThermalExpansionEigenstrain
block = SiC
thermal_expansion_coeff = 4.9e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[SiC_elasticity]
type = ComputeIsotropicElasticityTensor
block = SiC
youngs_modulus = 3.4e11
poissons_ratio = .13
[]
[SiC_creep]
type = MonolithicSiCCreepUpdate
block = SiC
temperature = temp
k_function = k_function
[]
[SiC_stress]
type = ComputeMultipleInelasticStress
block = SiC
tangent_operator = elastic
inelastic_models = 'SiC_creep'
[]
[SiC_temp]
type = HeatConductionMaterial
block = SiC
thermal_conductivity = 13.9 # J/m-s-K
specific_heat = 620.0 # J/kg-K
[]
[SiC_den]
type = StrainAdjustedDensity
strain_free_density = 3180.0 # kg/m^3
block = SiC
[]
[SiC_conc]
type = ArrheniusDiffusionCoef
block = SiC
d1 = 5.5e-14 # m^2/s
d1_function = d1_function
d1_function_variable = fluence
q1 = 125.0e+3 # J/mol
d2 = 1.6e-2 # m^2/s
q2 = 514.0e+3 # J/mol
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[OPyC_eigenstrain]
type = PyCIrradiationEigenstrain
block = OPyC
pyc_type = dense
eigenstrain_name = 'swelling_strain'
[]
[OPyC_thermal_strain]
type = ComputeThermalExpansionEigenstrain
block = OPyC
thermal_expansion_coeff = 5.65e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[OPyC_elasticity]
type = ComputeIsotropicElasticityTensor
block = OPyC
youngs_modulus = 4.74e10
poissons_ratio = .23
[]
[OPyC_conc]
type = ArrheniusDiffusionCoef
block = OPyC
d1 = 6.3e-8 # m^2/s
q1 = 222.0e+3 # J/mol
d2 = 0.0
q2 = 0.0
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[]
[Dampers]
[temp]
type = MaxIncrement
variable = temp
max_increment = 50
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = 'none'
nl_rel_tol = 5e-4
nl_abs_tol = 1e-9
nl_max_its = 15
l_tol = 1e-3
l_max_its = 50
start_time = 0.0
end_time = 85.3682e6
dt = 100
dtmax = 2e6
dtmin = 1
[TimeStepper]
type = IterationAdaptiveDT
dt = 100
optimal_iterations = 6
growth_factor = 1.5
linear_iteration_ratio = 100
time_t = '0 76e6 76.001e6 84.641e6 84.6482e6'
time_dt = '20 20 20 20 20'
[]
[Predictor]
type = SimplePredictor
scale = 1
skip_times_old = '0 76e6 76.001e6 84.641e6 84.6482e6'
[]
[]
[Outputs]
perf_graph = true
exodus = true
[console]
type = Console
max_rows = 25
[]
[csv]
type = CSV
sync_times = '100 6308007 75696087'
sync_only = true
[]
[]
[Postprocessors]
[Cs_release]
type = SideIntegralMassFlux
variable = conc
boundary = exterior
execute_on = timestep_end
[]
[dt]
type = TimestepSize
execute_on = timestep_end
[]
[fis_gas_produced] # fission gas produced (moles)
type = ElementIntegralFisGasGeneratedSifgrs
block = fuel
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[fis_gas_released] # fission gas released to plenum (moles)
type = ElementIntegralFisGasReleasedSifgrs
block = fuel
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[volumeTotal]
type = InternalVolume
boundary = exterior
scale_factor = -1
execute_on = 'initial timestep_end'
[]
[volumeFuel]
type = InternalVolume
boundary = fuel_outer_boundary
scale_factor = -1
execute_on = 'initial timestep_end'
[]
[volumeGas]
type = InternalVolume
boundary = 'buffer_IPyC IPyC_buffer'
# ro = 3.125e-4
# ri = 2.125e-4
# vb = 4/3*pi*(ro^3-ri^3) = 8.76e-11
# buffer density = 1000
# PyC density = 1900
# fill ratio = 10/19
# vb*10/19 = 4.6e-11
# Must remove 4.6e-11 m^3 from the volume
addition = -4.6e-11
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[volumeBufferShell]
type = InternalVolume
boundary = 'buffer_IPyC IPyC_buffer'
execute_on = 'initial timestep_end'
[]
[ave_temp_interior]
type = SideAverageValue
boundary = 'buffer_IPyC IPyC_buffer'
variable = temp
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
# Postprocessors for CO production
[total_fission_rate]
type = ElementIntegralPower
variable = temp
fission_rate = fission_rate
block = fuel
energy_per_fission = 1.0
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[total_fissions]
type = TimeIntegratedPostprocessor
value = total_fission_rate
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[avg_surface_temp]
type = SideAverageValue
variable = temp
boundary = exterior
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[time_int_surf_temp]
type = TimeIntegratedPostprocessor
value = avg_surface_temp
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[co_production]
type = CarbonMonoxideProduction
total_fissions = total_fissions
time_integrated_triso_temperature = time_int_surf_temp
initial_enrichment = 0.14029
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[num_lin_it]
type = NumLinearIterations
[]
[num_nonlin_it]
type = NumNonlinearIterations
[]
[tot_lin_it]
type = CumulativeValuePostprocessor
postprocessor = num_lin_it
[]
[tot_nonlin_it]
type = CumulativeValuePostprocessor
postprocessor = num_nonlin_it
[]
[alive_time]
type = PerfGraphData
section_name = Root
data_type = TOTAL
[]
[buffer_integral_flux]
type = SideDiffusiveFluxIntegral
variable = conc
boundary = buffer_IPyC
diffusivity = arrhenius_diffusion_coef
[]
[IPyC_integral_flux]
type = SideDiffusiveFluxIntegral
variable = conc
boundary = IPyC_buffer
diffusivity = arrhenius_diffusion_coef
[]
[buffer_partial_pressure]
type = SideAverageMaterialProperty
property = partial_pressure
boundary = buffer_IPyC
[]
[IPyC_partial_pressure]
type = SideAverageMaterialProperty
property = partial_pressure
boundary = IPyC_buffer
[]
[integral_flux_error]
type = FunctionValuePostprocessor
function = integral_flux_error
[]
[partial_pressure_error]
type = FunctionValuePostprocessor
function = partial_pressure_error
[]
[integral_Cs_release]
type = TimeIntegratedPostprocessor
value = Cs_release
[]
[Cs_production]
type = FunctionValuePostprocessor
function = power_history
scale_factor = 1.22e-5 # units of moles/m**3-s
[]
[time_integral_Cs_production]
type = TimeIntegratedPostprocessor
value = Cs_production
[]
[volumeFuel_initial]
type = InternalVolume
boundary = fuel_outer_boundary
scale_factor = -1
execute_on = initial
[]
[integral_Cs_production]
type = ParsedPostprocessor
pp_names = 'time_integral_Cs_production volumeFuel_initial'
expression = 'time_integral_Cs_production * volumeFuel_initial'
[]
[Cs_release_fraction]
type = ParsedPostprocessor
pp_names = 'integral_Cs_release integral_Cs_production'
expression = 'integral_Cs_release / integral_Cs_production'
[]
[]
[VectorPostprocessors]
[temperaturevpp]
type = SideValueSampler
boundary = 11
variable = temp
sort_by = x
outputs = 'csv'
use_displaced_mesh = true
[]
[]
(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Sr.i)
[InterfaceKernels]
[fuel_matrix]
type = SorptionIsothermGapInterface
variable = conc_Sr
neighbor_var = conc_Sr
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef_Sr
use_flux_penalty = true
flux_penalty = 1e3
boundary = 'DRV_IR'
[]
[matrix_graphite]
type = SorptionIsothermGapInterface
variable = conc_Sr
neighbor_var = conc_Sr
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef_Sr
use_flux_penalty = true
flux_penalty = 1e3
boundary = 'IR_OR'
[]
[graphite_sink]
type = SorptionIsothermGapInterface
variable = conc_Sr
neighbor_var = conc_Sr
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef_Sr
use_flux_penalty = true
flux_penalty = 1e3
boundary = 'OR_SR'
[]
[]
[Variables]
[conc_Sr]
scaling = 1e12 #1e18
[]
[]
[AuxVariables]
[Sr_diff_coef]
order = CONSTANT
family = MONOMIAL
[]
[]
[Kernels]
[mass_Sr_dt]
type = TimeDerivative
variable = conc_Sr
extra_vector_tags = 'ref'
[]
[mass_Sr]
type = ArrheniusDiffusion
variable = conc_Sr
arrhenius_prpty_name = arrhenius_diffusion_coef_Sr
extra_vector_tags = 'ref'
[]
[conc_source_dtf]
type = BodyForce
variable = conc_Sr
block = 'DTF'
function = DTF_MassIsotope
extra_vector_tags = 'ref'
[]
[conc_source_drv]
type = BodyForce
variable = conc_Sr
block = 'DRV'
function = DRV_MassIsotope
extra_vector_tags = 'ref'
[]
[]
[AuxKernels]
[Sr_diff_coef]
type = MaterialRealAux
variable = Sr_diff_coef
property = arrhenius_diffusion_coef_Sr
execute_on = timestep_end
[]
[]
[BCs]
[freesurf_conc_Sr]
type = DirichletBC
variable = conc_Sr
boundary = 'right'
value = 0.0
extra_vector_tags = 'ref'
[]
[]
[Materials]
[matrix_conc_Sr]
type = ArrheniusDiffusionCoef
d1 = 1.00e-2 # m^2/s
q1 = 303000 # J/mol
temperature = temperature
arrhenius_prpty_name = arrhenius_diffusion_coef_Sr
[]
[graphite_conc_Sr]
type = ArrheniusDiffusionCoef
d1 = 1.70e-2 # m^2/s
q1 = 268000 # J/mol
temperature = temperature
arrhenius_prpty_name = arrhenius_diffusion_coef_Sr
[]
[Sr_matrix_partial_pressure]
type = SorptionPartialPressure
A = 54.3
B = -149000
D = -8.52
E = 28500
d1 = 3.13
d2 = 0
unit_scale = 1e3 # convert from mol to mmol
density = density # convert from mmol/m^3 to mmol/kg
concentration = conc_Sr
temperature = temperature
output_properties = partial_pressure
[]
[Sr_graphite_partial_pressure]
type = SorptionPartialPressure
A = 19.4
B = -40100
D = -0.32
E = 4090
d1 = -2.12
d2 = 0
unit_scale = 1e3 # convert from mol to mmol
density = density # convert from mmol/m^3 to mmol/kg
concentration = conc_Sr
temperature = temperature
output_properties = partial_pressure
[]
[]
[Postprocessors]
[release_Sr_inc_pebble]
type = SideIntegralMassFlux
variable = conc_Sr
boundary = right
arrhenius_prpty_name = arrhenius_diffusion_coef_Sr
execute_on = 'initial timestep_end'
[]
[released_Sr]
type = TimeIntegratedPostprocessor
value = release_Sr_inc_pebble
execute_on = 'initial timestep_end'
[]
[retained_Sr]
type = ElementIntegralVariablePostprocessor
variable = conc_Sr
[]
# check partial pressure equality and integral flux conservation from DRV to IR
[flux_DRV_outer]
type = SideDiffusiveFluxIntegral
boundary = 'DRV_IR'
variable = conc_Sr
diffusivity = arrhenius_diffusion_coef_Sr
[]
[flux_IR_inner]
type = SideDiffusiveFluxIntegral
boundary = 'IR_DRV'
variable = conc_Sr
diffusivity = arrhenius_diffusion_coef_Sr
[]
[flux_DRV_IR_error]
type = FunctionValuePostprocessor
function = flux_DRV_IR_error
[]
# check partial pressure equality and integral flux conservation from IR to OR
[flux_IR_outer]
type = SideDiffusiveFluxIntegral
boundary = 'IR_OR'
variable = conc_Sr
diffusivity = arrhenius_diffusion_coef_Sr
[]
[flux_OR_inner]
type = SideDiffusiveFluxIntegral
boundary = 'OR_IR'
variable = conc_Sr
diffusivity = arrhenius_diffusion_coef_Sr
[]
[flux_IR_OR_error]
type = FunctionValuePostprocessor
function = flux_IR_OR_error
[]
# check partial pressure equality and integral flux conservation from OR to SR
[flux_OR_outer]
type = SideDiffusiveFluxIntegral
boundary = 'OR_SR'
variable = conc_Sr
diffusivity = arrhenius_diffusion_coef_Sr
[]
[flux_SR_inner]
type = SideDiffusiveFluxIntegral
boundary = 'SR_OR'
variable = conc_Sr
diffusivity = arrhenius_diffusion_coef_Sr
[]
[flux_OR_SR_error]
type = FunctionValuePostprocessor
function = flux_OR_SR_error
[]
[]
[VectorPostprocessors]
[Sr_conc]
type = LineValueSampler
variable = conc_Sr
start_point = '0.0 0 0.0'
num_points = 1000
use_displaced_mesh = false
execute_on = final
sort_by = x
outputs = line_plot_Sr
[]
[]
[Outputs]
[line_plot_Sr]
type = CSV
execute_on = 'FINAL'
time_step_interval = 1
create_final_symlink = true
file_base = capsule_radial
[]
[]
(test/tests/sorption_interface/test.i)
# This test is to verify the implementation of SorptionIsothermGapInterface and its AD counterpart.
# It contains two 2D blocks separated by a continuous interface.
# SorptionIsothermGapInterface is used to enforce sorption and preserve flux between the blocks.
# Checks are performed to verify concentration conservation, sorption behavior, and flux preservation.
[GlobalParams]
order = FIRST
family = LAGRANGE
[]
[Mesh]
[gen]
type = GeneratedMeshGenerator
nx = 10
ny = 10
dim = 2
[]
[block1]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 0'
top_right = '0.5 1 0'
input = gen
[]
[block2]
type = SubdomainBoundingBoxGenerator
block_id = 2
bottom_left = '0.5 0 0'
top_right = '1 1 0'
input = block1
[]
[breakmesh]
input = block2
type = BreakMeshByBlockGenerator
block_pairs = '1 2'
split_interface = true
add_interface_on_two_sides = true
[]
[]
[Variables]
[u]
[]
[temperature]
[]
[]
[Kernels]
[u]
type = MatDiffusion
variable = u
diffusivity = diffusivity
[]
[temp]
type = HeatConduction
variable = temperature
[]
[]
[BCs]
[left_u]
type = DirichletBC
value = 1
variable = u
boundary = left
[]
[right_u]
type = DirichletBC
value = 1
variable = u
boundary = right
[]
[left_temp]
type = DirichletBC
value = 1100
variable = temperature
boundary = left
[]
[right_temp]
type = DirichletBC
value = 0
variable = temperature
boundary = right
[]
[block1_2_temp]
type = DirichletBC
value = 1000
variable = temperature
boundary = Block1_Block2
[]
[block2_1_temp]
type = DirichletBC
value = 900
variable = temperature
boundary = Block2_Block1
[]
[]
[InterfaceKernels]
[interface]
type = SorptionIsothermGapInterface
variable = u
neighbor_var = u
interface_configuration = solid_solid
partial_pressure_name = partial_pressure
sorption_penalty = 1e10
boundary = Block1_Block2
[]
[]
[Materials]
[properties_1]
type = GenericConstantMaterial
prop_names = 'thermal_conductivity diffusivity density'
prop_values = '1 1 1'
block = 1
[]
[partial_pressure_1]
type = SorptionPartialPressure
A = 19.3
B = -47300
D = 1.51
E = 4340
d1 = 3.4
d2 = 6.15e-4
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = 1
outputs = 'all'
output_properties = partial_pressure
[]
[properties_2]
type = GenericConstantMaterial
prop_names = 'thermal_conductivity diffusivity density'
prop_values = '2 2 1'
block = 2
[]
[partial_pressure_2]
type = SorptionPartialPressure
A = 24
B = -35700
D = -1.56
E = 6120
d1 = 2.04
d2 = 1.79e-3
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = 2
outputs = 'all'
output_properties = partial_pressure
[]
[]
[Functions]
[u_mid_diff]
type = ParsedFunction
symbol_names = 'u_mid_inner u_mid_outer'
symbol_values = 'u_mid_inner u_mid_outer'
expression = '(abs(u_mid_outer) - abs(u_mid_inner)) / abs(u_mid_inner)'
[]
[pressure_mid_error]
type = ParsedFunction
symbol_names = 'pressure_mid_inner pressure_mid_outer'
symbol_values = 'pressure_mid_inner pressure_mid_outer'
expression = '(abs(pressure_mid_outer) - abs(pressure_mid_inner)) / abs(pressure_mid_inner)'
[]
[flux_error]
type = ParsedFunction
symbol_names = 'flux_inner flux_outer'
symbol_values = 'flux_inner flux_outer'
expression = '(abs(flux_outer) - abs(flux_inner)) / abs(flux_inner)'
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Steady
solve_type = NEWTON
automatic_scaling = true
scaling_group_variables = 'u; temperature'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = none
nl_rel_tol = 1e-15
nl_abs_tol = 1e-9
l_tol = 1e-3
l_max_its = 50
nl_max_its = 50
[]
[Postprocessors]
[u_mid_inner] # verify sorption behavior
type = PointValue
variable = u
point = '0.4999 0.5 0'
outputs = csv
[]
[u_mid_outer]
type = PointValue
variable = u
point = '0.5001 0.5 0'
outputs = csv
[]
[u_mid_diff]
type = FunctionValuePostprocessor
function = u_mid_diff
outputs = console
[]
[temperature_mid_inner]
type = PointValue
variable = temperature
point = '0.4999 0.5 0'
outputs = csv
[]
[temperature_mid_outer]
type = PointValue
variable = temperature
point = '0.5001 0.5 0'
outputs = csv
[]
[pressure_mid_inner]
type = PointValue
variable = partial_pressure
point = '0.5001 0.5 0'
outputs = csv
[]
[pressure_mid_outer]
type = PointValue
variable = partial_pressure
point = '0.5001 0.5 0'
outputs = csv
[]
[pressure_mid_error]
type = FunctionValuePostprocessor
function = pressure_mid_error
outputs = console
[]
[flux_inner] # verify flux preservation
type = SideDiffusiveFluxIntegral
variable = u
boundary = Block1_Block2
diffusivity = diffusivity
outputs = csv
[]
[flux_outer]
type = SideDiffusiveFluxIntegral
variable = u
boundary = Block2_Block1
diffusivity = diffusivity
outputs = csv
[]
[flux_error]
type = FunctionValuePostprocessor
function = flux_error
outputs = console
[]
[]
[Outputs]
exodus = true
csv = true
[]
[Debug]
show_var_residual_norms = true
[]
(test/tests/sorption_interface/adsorption.i)
# This test is verifies the implementation of SorptionIsothermGapInterface with
# interface_configuration = solid_gas. It contains a solid having a single
# interface with a gas. The gas acts as a source of u. u adsorbs onto the sink.
# Temperature increases over time to verify adsorption behavior and species
# conservation. For SorptionIsothermGapInterface, the solid is always the
# element and the gas is always the neighbor.
# _________________
# | sink | gap |
# -----------------
[GlobalParams]
order = SECOND
[]
[Mesh]
[gen]
type = GeneratedMeshGenerator
nx = 10
ny = 2
ymax = 0.2
dim = 2
elem_type = QUAD8
[]
[sink]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 0'
top_right = '0.5 0.2 0'
input = gen
[]
[gap]
type = SubdomainBoundingBoxGenerator
block_id = 2
bottom_left = '0.5 0 0'
top_right = '1 0.2 0'
input = sink
[]
[break]
type = BreakMeshByBlockGenerator
block_pairs = '1 2'
split_interface = true
add_interface_on_two_sides = true
input = gap
[]
[]
[Variables]
[u]
[]
[]
[AuxVariables]
[temperature]
[]
[]
[Kernels]
[u_dt]
type = TimeDerivative
variable = u
[]
[u]
type = MatDiffusion
variable = u
diffusivity = diffusivity
[]
[]
[ICs]
[u_gap]
type = ConstantIC
variable = u
value = 1
block = 2
[]
[]
[AuxKernels]
[temperature]
type = FunctionAux
variable = temperature
function = '(4000 - 1000 * x) / 2.0 * (2.0 - t + 0.1)'
execute_on = 'INITIAL TIMESTEP_BEGIN LINEAR NONLINEAR TIMESTEP_END'
[]
[]
[InterfaceKernels]
[u_sink_gap]
type = SorptionIsothermGapInterface
variable = u
neighbor_var = u
interface_configuration = solid_gas
temperature = temperature
partial_pressure_name = partial_pressure
sorption_penalty = 1e10
boundary = Block1_Block2
[]
[]
[Materials]
[properties_solid]
type = GenericConstantMaterial
prop_names = 'density diffusivity'
prop_values = '1 0.1'
block = 1
[]
[properties_gas]
type = GenericConstantMaterial
prop_names = diffusivity
prop_values = 1
block = 2
[]
[partial_pressure]
type = SorptionPartialPressure
A = 19.3
B = -47300
D = 1.51
E = 4340
d1 = 3.4
d2 = 6.15e-4
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = 1
outputs = 'all'
output_properties = partial_pressure
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
compute_scaling_once = false
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = none
nl_rel_tol = 1e-15
nl_abs_tol = 1e-15
l_tol = 1e-3
l_max_its = 50
nl_max_its = 50
dt = 0.1
end_time = 2
[]
[Functions]
[surface_u_exact]
type = ParsedFunction
symbol_names = 'P T R'
symbol_values = 'surface_pressure surface_temperature 8.31446261815324'
expression = '1.0 / R / T * P'
[]
[]
[Postprocessors]
[surface_pressure]
type = SideAverageMaterialProperty
property = partial_pressure
boundary = Block1_Block2
[]
[surface_temperature]
type = SideAverageValue
variable = temperature
boundary = Block1_Block2
[]
[surface_u]
type = SideAverageValue
variable = u
boundary = Block2_Block1
[]
[surface_u_exact]
type = FunctionValuePostprocessor
function = surface_u_exact
[]
[u_sink_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = 1
[]
[u_gap_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = 2
[]
[u_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = '1 2'
[]
[]
[Outputs]
exodus = true
csv = true
[]
[Debug]
show_var_residual_norms = true
[]
(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Ag.i)
[InterfaceKernels]
[fuel_matrix]
type = SorptionIsothermGapInterface
variable = conc_Ag
neighbor_var = conc_Ag
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef_Ag
use_flux_penalty = true
flux_penalty = 1e3
boundary = 'DRV_IR'
[]
[matrix_graphite]
type = SorptionIsothermGapInterface
variable = conc_Ag
neighbor_var = conc_Ag
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef_Ag
use_flux_penalty = true
flux_penalty = 1e3
boundary = 'IR_OR'
[]
[graphite_sink]
type = SorptionIsothermGapInterface
variable = conc_Ag
neighbor_var = conc_Ag
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef_Ag
use_flux_penalty = true
flux_penalty = 1e3
boundary = 'OR_SR'
[]
[]
[Variables]
[conc_Ag]
scaling = 1e12 #1e18
[]
[]
[AuxVariables]
[Ag_diff_coef]
order = CONSTANT
family = MONOMIAL
[]
[]
[Kernels]
[mass_Ag_dt]
type = TimeDerivative
variable = conc_Ag
extra_vector_tags = 'ref'
[]
[mass_Ag]
type = ArrheniusDiffusion
variable = conc_Ag
arrhenius_prpty_name = arrhenius_diffusion_coef_Ag
extra_vector_tags = 'ref'
[]
[conc_source_dtf]
type = BodyForce
variable = conc_Ag
block = 'DTF'
function = DTF_MassIsotope
extra_vector_tags = 'ref'
[]
[conc_source_drv]
type = BodyForce
variable = conc_Ag
block = 'DRV'
function = DRV_MassIsotope
extra_vector_tags = 'ref'
[]
[]
[AuxKernels]
[Ag_diff_coef]
type = MaterialRealAux
variable = Ag_diff_coef
property = arrhenius_diffusion_coef_Ag
execute_on = timestep_end
[]
[]
[BCs]
[freesurf_conc_Ag]
type = DirichletBC
variable = conc_Ag
boundary = 'right'
value = 0.0
extra_vector_tags = 'ref'
[]
[]
[Materials]
[matrix_conc_Ag]
type = ArrheniusDiffusionCoef
d1 = 1.6 # m^2/s
q1 = 258000 # J/mol
temperature = temperature
arrhenius_prpty_name = arrhenius_diffusion_coef_Ag
[]
[graphite_conc_Ag]
type = ArrheniusDiffusionCoef
d1 = 1.38e-2 # m^2/s
q1 = 226000 # J/mol
temperature = temperature
arrhenius_prpty_name = arrhenius_diffusion_coef_Ag
[]
[Ag_matrix_partial_pressure]
type = SorptionPartialPressure
A = 19.33
B = -47290
D = 1.518
E = 4338
d1 = 3.397
d2 = 6.15e-4
unit_scale = 1e3 # convert from mol to mmol
density = density # convert from mmol/m^3 to mmol/kg
concentration = conc_Ag
temperature = temperature
output_properties = partial_pressure
[]
[Ag_graphite_partial_pressure]
type = SorptionPartialPressure
A = 24
B = -35700
D = -1.56
E = 6120
d1 = 2.04
d2 = 1.79e-3
unit_scale = 1e3 # convert from mol to mmol
density = density # convert from mmol/m^3 to mmol/kg
concentration = conc_Ag
temperature = temperature
output_properties = partial_pressure
[]
[]
[Postprocessors]
[release_Ag_inc_pebble]
type = SideIntegralMassFlux
variable = conc_Ag
boundary = right
arrhenius_prpty_name = arrhenius_diffusion_coef_Ag
execute_on = 'initial timestep_end'
[]
[released_Ag]
type = TimeIntegratedPostprocessor
value = release_Ag_inc_pebble
execute_on = 'initial timestep_end'
[]
[retained_Ag]
type = ElementIntegralVariablePostprocessor
variable = conc_Ag
[]
# check partial pressure equality and integral flux conservation from DRV to IR
[flux_DRV_outer]
type = SideDiffusiveFluxIntegral
boundary = 'DRV_IR'
variable = conc_Ag
diffusivity = arrhenius_diffusion_coef_Ag
[]
[flux_IR_inner]
type = SideDiffusiveFluxIntegral
boundary = 'IR_DRV'
variable = conc_Ag
diffusivity = arrhenius_diffusion_coef_Ag
[]
[flux_DRV_IR_error]
type = FunctionValuePostprocessor
function = flux_DRV_IR_error
[]
# check partial pressure equality and integral flux conservation from IR to OR
[flux_IR_outer]
type = SideDiffusiveFluxIntegral
boundary = 'IR_OR'
variable = conc_Ag
diffusivity = arrhenius_diffusion_coef_Ag
[]
[flux_OR_inner]
type = SideDiffusiveFluxIntegral
boundary = 'OR_IR'
variable = conc_Ag
diffusivity = arrhenius_diffusion_coef_Ag
[]
[flux_IR_OR_error]
type = FunctionValuePostprocessor
function = flux_IR_OR_error
[]
# check partial pressure equality and integral flux conservation from OR to SR
[flux_OR_outer]
type = SideDiffusiveFluxIntegral
boundary = 'OR_SR'
variable = conc_Ag
diffusivity = arrhenius_diffusion_coef_Ag
[]
[flux_SR_inner]
type = SideDiffusiveFluxIntegral
boundary = 'SR_OR'
variable = conc_Ag
diffusivity = arrhenius_diffusion_coef_Ag
[]
[flux_OR_SR_error]
type = FunctionValuePostprocessor
function = flux_OR_SR_error
[]
[]
[VectorPostprocessors]
[Ag_conc]
type = LineValueSampler
variable = conc_Ag
start_point = '0.0 0 0.0'
num_points = 1000
use_displaced_mesh = false
execute_on = final
sort_by = x
outputs = line_plot_Ag
[]
[]
[Outputs]
[line_plot_Ag]
type = CSV
execute_on = 'FINAL'
time_step_interval = 1
create_final_symlink = true
file_base = capsule_radial
[]
[]
(examples/TRISO/accident_simulation/triso2D_accident_ad.i)
# This example is 2D-RZ analysis of a TRISO fuel particle. Fully coupled
# heat transfer and solid mechanics, plus diffusion of the fission product
# species cesium (Cs) are simulated. The mesh includes contact surfaces
# between the buffer and IPyC layers to facilitate a gap opening between
# these layers. These surfaces are initially in mechanical contact but
# are assumed to have no strength in tension. A coarse mesh is used to
# provide a short run time.
# The calculation simulates fuel-life in three steps. The first step is an
# irradiation period, where constant power and a fixed particle surface
# temperature (1500 K) are assumed over a lifetime of 76 Ms (2.4 yrs).
# For the second step, fuel removal and storage are simulated by setting
# the reactor power and Cs source terms to zero, reducing the particle
# surface temperature to ambient (300 K), and then holding it
# for 100 days. A third and final step simulates accident
# behavior by increasing the particle surface temperature from ambient
# to 2073 K over 2 hrs, and then holding it at this elevated temperature
# for an additional 200 hrs. At the particle outer boundary, the Cs
# concentration is held at zero and the pressure at ambient during the
# entire simulation. The particle is assumed to be stress-free at an
# initial temperature of 1500 K.
#
# Details about this simulation are given in Section 4 of the following
# article: J. D. Hales, R. L. Williamson, S. R. Novascone, D. M. Perez,
# B. W. Spencer and G. Pastore, "Multidimensional multiphysics simulation
# of TRISO particle fuel", Journal of Nuclear Materials, Vol. 443, p. 531,
# 2013.
# This is a version using a thermomechanical mortar approach. It uses
# Automatic Differentiation classes and models gap mass transfer using
# flux preserving and sorption mortar constraints. Sorption constants are
# given in Table 1 of the following article: A. Londono-Hurtado, I.
# Szlufarska, R. Bratton and D. Morgan, "A review of fission product
# sorption in carbon structures", Journal of Nuclear Materials, Vol. 426,
# p. 254, 2012.
initial_fuel_density = 11000.0
[GlobalParams]
order = SECOND
family = LAGRANGE
displacements = 'disp_x disp_y'
flux_conversion_factor = 0.85
use_automatic_differentiation = true
[]
[Mesh]
coord_type = RZ
[file]
type = FileMeshGenerator
file = triso2Dmed.e
[]
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
converge_on = 'disp_x disp_y temp conc'
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[temp]
initial_condition = 1500.0
[]
[conc]
initial_condition = 0.0
[]
[conc_lm]
block = pellet_clad_mechanical_secondary_subdomain
[]
[conc_dx_lm]
block = pellet_clad_mechanical_secondary_subdomain
[]
[conc_dy_lm]
block = pellet_clad_mechanical_secondary_subdomain
[]
[]
[AuxVariables]
[fission_rate]
block = fuel
order = CONSTANT
family = MONOMIAL
[]
[fluence]
order = CONSTANT
family = MONOMIAL
[]
[burnup]
block = fuel
order = CONSTANT
family = MONOMIAL
[]
[creep_xx]
order = CONSTANT
family = MONOMIAL
[]
[creep_yy]
order = CONSTANT
family = MONOMIAL
[]
[creep_zz]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[power_history]
type = PiecewiseLinear
x = '0 76e6 76.001e6'
y = '1 1 0'
[]
[temp_bc]
type = PiecewiseLinear
x = '0 76e6 76.001e6 84.641e6 84.6482e6'
y = '1500 1500 300 300 2073'
[]
[k_function]
type = PiecewiseLinear
x = '0 200e6'
y = '4e-37 4e-37'
[]
[d1_function]
type = ParsedFunction
expression = 'exp(t/4.5e25)'
[]
[integral_flux_error]
type = ParsedFunction
symbol_names = 'buffer_integral_flux IPyC_integral_flux'
symbol_values = 'buffer_integral_flux IPyC_integral_flux'
expression = 'IPyC_integral_flux + buffer_integral_flux'
[]
[partial_pressure_error]
type = ParsedFunction
symbol_names = 'buffer_partial_pressure IPyC_partial_pressure'
symbol_values = 'buffer_partial_pressure IPyC_partial_pressure'
expression = 'IPyC_partial_pressure - buffer_partial_pressure'
[]
[]
[Physics/SolidMechanics/QuasiStatic]
generate_output = 'stress_xx stress_yy stress_zz stress_xy stress_yz stress_zx hydrostatic_stress'
strain = FINITE
incremental = true
add_variables = false
[default]
block = 'fuel buffer IPyC OPyC'
eigenstrain_names = 'thermal_strain swelling_strain'
extra_vector_tags = 'ref'
[]
[SiC]
block = 'SiC'
eigenstrain_names = 'thermal_strain'
extra_vector_tags = 'ref'
[]
[]
[Kernels]
[heat_ie]
type = ADHeatConductionTimeDerivative
variable = temp
extra_vector_tags = 'ref'
block = 'fuel buffer IPyC SiC OPyC'
[]
[heat]
type = ADHeatConduction
variable = temp
extra_vector_tags = 'ref'
block = 'fuel buffer IPyC SiC OPyC'
[]
[heat_source]
type = ADNeutronHeatSource
variable = temp
block = fuel
energy_per_fission = 3.2e-11 # units of J/fission
fission_rate = fission_rate
extra_vector_tags = 'ref'
[]
[mass_ie]
type = ADTimeDerivative
variable = conc
extra_vector_tags = 'ref'
block = 'fuel buffer IPyC SiC OPyC'
[]
[mass]
type = ADArrheniusDiffusion
variable = conc
extra_vector_tags = 'ref'
block = 'fuel buffer IPyC SiC OPyC'
[]
[mass_source]
type = ADBodyForce
variable = conc
function = power_history
value = 1.22e-5 # units of moles/m**3-s
block = fuel
extra_vector_tags = 'ref'
[]
[mass_decay]
type = Decay
variable = conc
radioactive_decay_constant = 7.297e-10 # units:(1/sec) The constant for Cesium
block = 'fuel buffer IPyC SiC OPyC'
extra_vector_tags = 'ref'
[]
[]
[AuxKernels]
[fission_rate]
type = FissionRateGeneral
fission_rate_formulation = GENERIC
variable = fission_rate
block = fuel
fission_rate_function = power_history
value = 3.89e19
execute_on = timestep_begin
[]
[fluence]
type = ADMaterialRealAux
property = fast_neutron_fluence
variable = fluence
[]
[burnup]
type = ADBurnupAux
variable = burnup
block = fuel
fission_rate = fission_rate
molecular_weight = 0.270 # units of kg/mole
execute_on = timestep_begin
density = ${initial_fuel_density}
[]
[creep_xx]
type = ADRankTwoAux
rank_two_tensor = creep_strain
variable = creep_xx
index_i = 0
index_j = 0
block = 'buffer IPyC SiC OPyC'
execute_on = timestep_end
[]
[creep_yy]
type = ADRankTwoAux
rank_two_tensor = creep_strain
variable = creep_yy
index_i = 1
index_j = 1
block = 'buffer IPyC SiC OPyC'
execute_on = timestep_end
[]
[creep_zz]
type = ADRankTwoAux
rank_two_tensor = creep_strain
variable = creep_zz
index_i = 2
index_j = 2
block = 'buffer IPyC SiC OPyC'
execute_on = timestep_end
[]
[]
[ThermalContactMortar]
[thermal]
secondary_variable = temp
primary_boundary = 15
secondary_boundary = 17
initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
gas_released = 'fis_gas_released co_production' # coupling to postprocessors which supply the fission gas addition, co addition
released_gas_types = 'Kr Xe;
CO'
released_fractions = '0.153 0.847;
1'
gap_geometry_type = CYLINDER
min_gap = 1e-7
max_gap = 50e-6
roughness_coef = 0.0
correct_edge_dropping = true
[]
[]
[Contact]
[pellet_clad_mechanical]
primary = 15
secondary = 17
model = frictionless
formulation = mortar
c_normal = 1.0e8
correct_edge_dropping = true
[]
[]
[Constraints]
[cesium_gap_value]
type = MassSorptionConstraint
variable = conc_lm
primary_variable = conc
primary_boundary = 15
primary_subdomain = pellet_clad_mechanical_primary_subdomain
secondary_variable = conc
secondary_boundary = 17
secondary_subdomain = pellet_clad_mechanical_secondary_subdomain
partial_pressure_name = partial_pressure
epsilon = 1e-4
correct_edge_dropping = true
[]
[cesium_gap_flux_x]
type = MassFluxConstraint
variable = conc_dx_lm
primary_variable = conc
diffusivity_primary = arrhenius_diffusion_coef
primary_boundary = 15
primary_subdomain = pellet_clad_mechanical_primary_subdomain
secondary_variable = conc
diffusivity_secondary = arrhenius_diffusion_coef
secondary_boundary = 17
secondary_subdomain = pellet_clad_mechanical_secondary_subdomain
component = 0
epsilon = 1e-5
correct_edge_dropping = true
[]
[cesium_gap_flux_y]
type = MassFluxConstraint
variable = conc_dy_lm
primary_variable = conc
diffusivity_primary = arrhenius_diffusion_coef
primary_boundary = 15
primary_subdomain = pellet_clad_mechanical_primary_subdomain
secondary_variable = conc
diffusivity_secondary = arrhenius_diffusion_coef
secondary_boundary = 17
secondary_subdomain = pellet_clad_mechanical_secondary_subdomain
component = 1
epsilon = 1e-5
correct_edge_dropping = true
[]
[]
[BCs]
# pin particle along symmetry planes
[no_disp_x]
type = ADDirichletBC
variable = disp_x
boundary = xzero
value = 0.0
extra_vector_tags = 'ref'
[]
[no_disp_y]
type = ADDirichletBC
variable = disp_y
boundary = yzero
value = 0.0
extra_vector_tags = 'ref'
[]
# fix temperature on free surface
[freesurf_temp]
type = ADFunctionDirichletBC
variable = temp
boundary = exterior
function = temp_bc
extra_vector_tags = 'ref'
[]
# fix concentration on free surface
[freesurf_conc]
type = ADDirichletBC
variable = conc
boundary = exterior
value = 0.0
extra_vector_tags = 'ref'
[]
[PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
[plenumPressure]
boundary = BufferGapVol
initial_pressure = 0
startup_time = 1.0e4
R = 8.3145
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 = volumeGas # coupling to post processor to get gas volume
material_input = 'fis_gas_released co_production' # coupling to post processor to get fission gas added, co added
output = plenum_pressure # coupling to post processor to output plenum/gap pressure
[]
[]
[]
[Materials]
[flux]
type = ADFastNeutronFlux
calculate_fluence = true
factor = 5e17
[]
[fission_gas_release] # Sifgrs fission gas release mode
type = ADUO2Sifgrs
block = fuel
temperature = temp
fission_rate = fission_rate # coupling to fission_rate aux variable
grain_radius_const = 5.0e-6
[]
[fuel_thermal]
type = ADUO2Thermal
thermal_conductivity_model = FINK_LUCUTA
block = fuel
temperature = temp
burnup = burnup
initial_porosity = 0.0
[]
[fuel_swelling]
type = ADUO2VolumetricSwellingEigenstrain
gas_swelling_model_type = MATPRO
block = fuel
temperature = temp
burnup = burnup
eigenstrain_name = 'swelling_strain'
initial_fuel_density = ${initial_fuel_density}
[]
[fuel_stress]
type = ADComputeFiniteStrainElasticStress
block = 'fuel'
[]
[fuel_elasticity]
type = ADComputeIsotropicElasticityTensor
block = fuel
youngs_modulus = 2.2e11
poissons_ratio = .345
[]
[fuel_thermal_strain]
type = ADComputeThermalExpansionEigenstrain
block = fuel
thermal_expansion_coeff = 10e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[fuel_den]
type = ADStrainAdjustedDensity
block = fuel
strain_free_density = ${initial_fuel_density} # kg/m^3
[]
[fuel_conc]
type = ADArrheniusDiffusionCoef
block = fuel
d1 = 5.6e-8 # m^2/s
q1 = 209.0e+3 # J/mol
d2 = 5.2e-4 # m^2/s
q2 = 362.0e+3 # J/mol
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[buffer_eigenstrain]
type = ADPyCIrradiationEigenstrain
block = buffer
pyc_type = buffer
eigenstrain_name = 'swelling_strain'
[]
[buffer_thermal_strain]
type = ADComputeThermalExpansionEigenstrain
block = buffer
thermal_expansion_coeff = 5.65e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[buffer_elasticity]
type = ADComputeIsotropicElasticityTensor
block = buffer
youngs_modulus = 2e10
poissons_ratio = .23
[]
[buffer_stress]
type = ADPyCCreep
block = buffer
temperature = temp
[]
[buffer_temp]
type = ADHeatConductionMaterial
block = buffer
thermal_conductivity = 0.5 # J/m-s-K
specific_heat = 720.0 # J/kg-K
[]
[buffer_den]
type = ADStrainAdjustedDensity
strain_free_density = 1000.0 #kg/m^3
block = buffer
[]
[buffer_conc]
type = ADArrheniusDiffusionCoef
block = buffer
d1 = 1.0e-12 # m^2/s
q1 = 0.0
d2 = 0.0
q2 = 0.0
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[buffer_partial_pressure]
type = ADSorptionPartialPressure
A = 19.33
B = -47290
D = 1.518
E = 4338
d1 = 3.397
d2 = 6.15e-4
unit_scale = 1e3 # convert from mol to mmol
density = 1000 # convert from mmol/m^3 to mmol/kg, using constant for compatibility with default AD derivative container size
concentration = conc
temperature = temp
block = buffer
outputs = 'all'
output_properties = partial_pressure
[]
[normal_vectors_triso]
type = NormalVectorsTRISO
block = 'IPyC OPyC buffer'
[]
[IPyC_eigenstrain]
type = ADPyCIrradiationEigenstrain
block = IPyC
pyc_type = dense
eigenstrain_name = 'swelling_strain'
[]
[IPyC_thermal_strain]
type = ADComputeThermalExpansionEigenstrain
block = IPyC
thermal_expansion_coeff = 5.65e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[IPyC_elasticity]
type = ADComputeIsotropicElasticityTensor
block = IPyC
youngs_modulus = 4.74e10
poissons_ratio = .23
[]
[IPyC_disp]
type = ADPyCCreep
block = 'IPyC OPyC'
temperature = temp
[]
[IPyC_temp]
type = ADHeatConductionMaterial
block = 'IPyC OPyC'
thermal_conductivity = 4.0
specific_heat = 720.0
[]
[IPyC_den]
type = ADStrainAdjustedDensity
block = 'IPyC OPyC'
strain_free_density = 1900.0
[]
[IPyC_conc]
type = ADArrheniusDiffusionCoef
block = IPyC
d1 = 6.3e-8
q1 = 222.0e+3
d2 = 0.0
q2 = 0.0
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[IPyC_partial_pressure]
type = ADSorptionPartialPressure
A = 19.33
B = -47290
D = 1.518
E = 4338
d1 = 3.397
d2 = 6.15e-4
unit_scale = 1e3 # convert from mol to mmol
density = 1900 # convert from mmol/m^3 to mmol/kg, using constant for compatibility with default AD derivative container size
concentration = conc
temperature = temp
block = IPyC
outputs = 'all'
output_properties = partial_pressure
[]
[SiC_thermal_strain]
type = ADComputeThermalExpansionEigenstrain
block = SiC
thermal_expansion_coeff = 4.9e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[SiC_elasticity]
type = ADComputeIsotropicElasticityTensor
block = SiC
youngs_modulus = 3.4e11
poissons_ratio = .13
[]
[SiC_creep]
type = ADMonolithicSiCCreepUpdate
block = SiC
temperature = temp
k_function = k_function
[]
[SiC_stress]
type = ADComputeMultipleInelasticStress
block = SiC
inelastic_models = 'SiC_creep'
[]
[SiC_temp]
type = ADHeatConductionMaterial
block = SiC
thermal_conductivity = 13.9 # J/m-s-K
specific_heat = 620.0 # J/kg-K
[]
[SiC_den]
type = ADStrainAdjustedDensity
strain_free_density = 3180.0 # kg/m^3
block = SiC
[]
[SiC_conc]
type = ADArrheniusDiffusionCoef
block = SiC
d1 = 5.5e-14 # m^2/s
d1_function = d1_function
d1_function_variable = fluence
q1 = 125.0e+3 # J/mol
d2 = 1.6e-2 # m^2/s
q2 = 514.0e+3 # J/mol
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[OPyC_eigenstrain]
type = ADPyCIrradiationEigenstrain
block = OPyC
pyc_type = dense
eigenstrain_name = 'swelling_strain'
[]
[OPyC_thermal_strain]
type = ADComputeThermalExpansionEigenstrain
block = OPyC
thermal_expansion_coeff = 5.65e-6
stress_free_temperature = 1500.0
eigenstrain_name = thermal_strain
temperature = temp
[]
[OPyC_elasticity]
type = ADComputeIsotropicElasticityTensor
block = OPyC
youngs_modulus = 4.74e10
poissons_ratio = .23
[]
[OPyC_conc]
type = ADArrheniusDiffusionCoef
block = OPyC
d1 = 6.3e-8 # m^2/s
q1 = 222.0e+3 # J/mol
d2 = 0.0
q2 = 0.0
gas_constant = 8.3143 # J/K-mol
temperature = temp
[]
[]
[Dampers]
[temp]
type = MaxIncrement
variable = temp
max_increment = 50
[]
[]
[Debug]
show_var_residual_norms = true
[]
[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_ksp_ew'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
petsc_options_value = 'lu superlu_dist 1e-5 NONZERO 1e-14'
snesmf_reuse_base = false
line_search = 'none'
nl_rel_tol = 5e-4
nl_abs_tol = 1e-10
nl_max_its = 20
l_max_its = 8
start_time = 0.0
end_time = 85.3682e6
dt = 100
dtmax = 2e6
dtmin = 1
[TimeStepper]
type = IterationAdaptiveDT
dt = 100
optimal_iterations = 10
growth_factor = 1.5
linear_iteration_ratio = 100
time_t = '0 76e6 76.001e6 84.641e6 84.6482e6'
time_dt = '20 20 20 20 20'
[]
[Predictor]
type = SimplePredictor
scale = 0.5
skip_times_old = '0 76e6 76.001e6 84.641e6 84.6482e6'
[]
[]
[Outputs]
perf_graph = true
exodus = true
[console]
type = Console
max_rows = 25
[]
[csv]
type = CSV
sync_times = '100 6308007 75696087'
sync_only = true
[]
[]
[Postprocessors]
[Cs_release]
type = ADSideDiffusiveFluxIntegral
variable = conc
diffusivity = arrhenius_diffusion_coef
boundary = exterior
execute_on = timestep_end
[]
[dt]
type = TimestepSize
execute_on = timestep_end
[]
[fis_gas_produced] # fission gas produced (moles)
type = ADElementIntegralFisGasGeneratedSifgrs
block = fuel
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[fis_gas_released] # fission gas released to plenum (moles)
type = ADElementIntegralFisGasReleasedSifgrs
block = fuel
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[volumeTotal]
type = InternalVolume
boundary = exterior
execute_on = 'initial timestep_end'
[]
[volumeFuel]
type = InternalVolume
boundary = fuel
execute_on = 'initial timestep_end'
[]
[volumeGas]
type = InternalVolume
boundary = BufferGapVol
# ro = 3.125e-4
# ri = 2.125e-4
# vb = 4/3*pi*(ro^3-ri^3) = 8.76e-11
# buffer density = 1000
# PyC density = 1900
# fill ratio = 10/19
# vb*10/19 = 4.6e-11
# Must remove 4.6e-11 m^3 from the volume
addition = -4.6e-11
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[volumeBufferShell]
type = InternalVolume
boundary = BufferGapVol
execute_on = 'initial timestep_end'
[]
[ave_temp_interior]
type = SideAverageValue
boundary = BufferGapVol
variable = temp
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
# Postprocessors for CO production
[total_fission_rate]
type = ElementIntegralPower
variable = temp
fission_rate = fission_rate
block = fuel
energy_per_fission = 1.0
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[total_fissions]
type = TimeIntegratedPostprocessor
value = total_fission_rate
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[avg_surface_temp]
type = SideAverageValue
variable = temp
boundary = exterior
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[time_int_surf_temp]
type = TimeIntegratedPostprocessor
value = avg_surface_temp
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[co_production]
type = CarbonMonoxideProduction
total_fissions = total_fissions
time_integrated_triso_temperature = time_int_surf_temp
initial_enrichment = 0.14029
execute_on = 'initial linear nonlinear timestep_begin timestep_end'
[]
[num_lin_it]
type = NumLinearIterations
[]
[num_nonlin_it]
type = NumNonlinearIterations
[]
[tot_lin_it]
type = CumulativeValuePostprocessor
postprocessor = num_lin_it
[]
[tot_nonlin_it]
type = CumulativeValuePostprocessor
postprocessor = num_nonlin_it
[]
[alive_time]
type = PerfGraphData
section_name = Root
data_type = TOTAL
[]
[buffer_integral_flux]
type = ADSideDiffusiveFluxIntegral
variable = conc
boundary = 17
diffusivity = arrhenius_diffusion_coef
[]
[IPyC_integral_flux]
type = ADSideDiffusiveFluxIntegral
variable = conc
boundary = 15
diffusivity = arrhenius_diffusion_coef
[]
[buffer_partial_pressure]
type = ADSideAverageMaterialProperty
property = partial_pressure
boundary = 17
[]
[IPyC_partial_pressure]
type = ADSideAverageMaterialProperty
property = partial_pressure
boundary = 15
[]
[integral_flux_error]
type = FunctionValuePostprocessor
function = integral_flux_error
[]
[partial_pressure_error]
type = FunctionValuePostprocessor
function = partial_pressure_error
[]
[integral_Cs_release]
type = TimeIntegratedPostprocessor
value = Cs_release
[]
[Cs_production]
type = FunctionValuePostprocessor
function = power_history
scale_factor = 1.22e-5 # units of moles/m**3-s
[]
[time_integral_Cs_production]
type = TimeIntegratedPostprocessor
value = Cs_production
[]
[volumeFuel_initial]
type = InternalVolume
boundary = fuel
execute_on = initial
[]
[integral_Cs_production]
type = ParsedPostprocessor
pp_names = 'time_integral_Cs_production volumeFuel_initial'
expression = 'time_integral_Cs_production * volumeFuel_initial'
[]
[Cs_release_fraction]
type = ParsedPostprocessor
pp_names = 'integral_Cs_release integral_Cs_production'
expression = 'integral_Cs_release / integral_Cs_production'
[]
[]
[VectorPostprocessors]
[temperaturevpp]
type = SideValueSampler
boundary = 11
variable = temp
sort_by = x
outputs = 'csv'
use_displaced_mesh = true
[]
[]
(assessment/TRISO/validation/AGR-34/SharedFiles/capsule_Cs.i)
[InterfaceKernels]
[fuel_matrix]
type = SorptionIsothermGapInterface
variable = conc_Cs
neighbor_var = conc_Cs
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef_Cs
use_flux_penalty = true
flux_penalty = 1e3
boundary = 'DRV_IR'
[]
[matrix_graphite]
type = SorptionIsothermGapInterface
variable = conc_Cs
neighbor_var = conc_Cs
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef_Cs
use_flux_penalty = true
flux_penalty = 1e3
boundary = 'IR_OR'
[]
[graphite_sink]
type = SorptionIsothermGapInterface
variable = conc_Cs
neighbor_var = conc_Cs
partial_pressure_name = partial_pressure
sorption_penalty = 1e5
diffusivity = arrhenius_diffusion_coef_Cs
use_flux_penalty = true
flux_penalty = 1e3
boundary = 'OR_SR'
[]
[]
[Variables]
[conc_Cs]
scaling = 1e12 #1e18
[]
[]
[AuxVariables]
[Cs_diff_coef]
order = CONSTANT
family = MONOMIAL
[]
[]
[Kernels]
[mass_Cs_dt]
type = TimeDerivative
variable = conc_Cs
extra_vector_tags = 'ref'
[]
[mass_Cs]
type = ArrheniusDiffusion
variable = conc_Cs
arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
extra_vector_tags = 'ref'
[]
[conc_source_dtf]
type = BodyForce
variable = conc_Cs
block = 'DTF'
function = DTF_MassIsotope
extra_vector_tags = 'ref'
[]
[conc_source_drv]
type = BodyForce
variable = conc_Cs
block = 'DRV'
function = DRV_MassIsotope
extra_vector_tags = 'ref'
[]
[]
[AuxKernels]
[Cs_diff_coef]
type = MaterialRealAux
variable = Cs_diff_coef
property = arrhenius_diffusion_coef_Cs
execute_on = timestep_end
[]
[]
[BCs]
[freesurf_conc_Cs]
type = DirichletBC
variable = conc_Cs
boundary = 'right'
value = 0.0
extra_vector_tags = 'ref'
[]
[]
[Materials]
[matrix_conc_Cs]
type = ArrheniusDiffusionCoef
d1 = 3.6e-4 #3.6e-4m^2/s
q1 = 189000 # J/mol
temperature = temperature
arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
[]
[graphite_conc_Cs]
type = ArrheniusDiffusionCoef
d1 = 1.7e-6 # 1.7e-6 m^2/s
q1 = 149000 # J/mol
temperature = temperature
arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
[]
[Cs_matrix_partial_pressure]
type = SorptionPartialPressure
A = 19.33
B = -47290
D = 1.518
E = 4338
d1 = 3.397
d2 = 6.15e-4
unit_scale = 1e3 # convert from mol to mmol
density = density # convert from mmol/m^3 to mmol/kg
concentration = conc_Cs
temperature = temperature
output_properties = partial_pressure
[]
[Cs_graphite_partial_pressure]
type = SorptionPartialPressure
A = 24.0
B = -35730.
D = -1.561
E = 6123.
d1 = -2.035
d2 = 1.79e-3
unit_scale = 1e3 # convert from mol to mmol
density = density # convert from mmol/m^3 to mmol/kg
concentration = conc_Cs
temperature = temperature
output_properties = partial_pressure
[]
[]
[Postprocessors]
[release_Cs_inc_pebble]
type = SideIntegralMassFlux
variable = conc_Cs
boundary = right
arrhenius_prpty_name = arrhenius_diffusion_coef_Cs
execute_on = 'initial timestep_end'
[]
[released_Cs]
type = TimeIntegratedPostprocessor
value = release_Cs_inc_pebble
execute_on = 'initial timestep_end'
[]
[retained_Cs]
type = ElementIntegralVariablePostprocessor
variable = conc_Cs
[]
# check partial pressure equality and integral flux conservation from DRV to IR
[flux_DRV_outer]
type = SideDiffusiveFluxIntegral
boundary = 'DRV_IR'
variable = conc_Cs
diffusivity = arrhenius_diffusion_coef_Cs
[]
[flux_IR_inner]
type = SideDiffusiveFluxIntegral
boundary = 'IR_DRV'
variable = conc_Cs
diffusivity = arrhenius_diffusion_coef_Cs
[]
[flux_DRV_IR_error]
type = FunctionValuePostprocessor
function = flux_DRV_IR_error
[]
# check partial pressure equality and integral flux conservation from IR to OR
[flux_IR_outer]
type = SideDiffusiveFluxIntegral
boundary = 'IR_OR'
variable = conc_Cs
diffusivity = arrhenius_diffusion_coef_Cs
[]
[flux_OR_inner]
type = SideDiffusiveFluxIntegral
boundary = 'OR_IR'
variable = conc_Cs
diffusivity = arrhenius_diffusion_coef_Cs
[]
[flux_IR_OR_error]
type = FunctionValuePostprocessor
function = flux_IR_OR_error
[]
# check partial pressure equality and integral flux conservation from OR to SR
[flux_OR_outer]
type = SideDiffusiveFluxIntegral
boundary = 'OR_SR'
variable = conc_Cs
diffusivity = arrhenius_diffusion_coef_Cs
[]
[flux_SR_inner]
type = SideDiffusiveFluxIntegral
boundary = 'SR_OR'
variable = conc_Cs
diffusivity = arrhenius_diffusion_coef_Cs
[]
[flux_OR_SR_error]
type = FunctionValuePostprocessor
function = flux_OR_SR_error
[]
[]
[VectorPostprocessors]
[Cs_conc]
type = LineValueSampler
variable = conc_Cs
start_point = '0.0 0 0.0'
num_points = 1000
use_displaced_mesh = false
execute_on = final
sort_by = x
outputs = line_plot_Cs
[]
[]
[Outputs]
[line_plot_Cs]
type = CSV
execute_on = 'FINAL'
time_step_interval = 1
create_final_symlink = true
file_base = capsule_radial
[]
[]
(test/tests/sorption_partial_pressure/test.i)
# This test is to verify the implementation of SorptionPartialPressure material.
# It compares calculated partial pressure parameters to analytical solutions.
[Problem]
solve = false
[]
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 1
[]
[]
[Functions]
[partial_pressure_avg_exact]
type = ParsedFunction
symbol_names = 'c T A B D E d1 d2'
symbol_values = 'c_avg T_avg 19.3 -47300 1.51 4340 3.4 6.15e-4'
expression = 'exp(A + B / T) * c ^ (D + E / T) + exp(A + B / T + (D - 1 + E / T) * (d1 - d2 * T)) * c'
[]
[partial_pressure_dc_avg_exact]
type = ParsedFunction
symbol_names = 'c T A B D E d1 d2'
symbol_values = 'c_avg T_avg 19.3 -47300 1.51 4340 3.4 6.15e-4'
expression = '(D + E / T) * exp(A + B / T) * c ^ (D + E / T - 1) + exp(A + B / T + (D - 1 + E / T) * (d1 - d2 * T))'
[]
[]
[AuxVariables]
[c]
order = FIRST
family = LAGRANGE
initial_condition = 1
[]
[T]
order = FIRST
family = LAGRANGE
initial_condition = 1000
[]
[]
[Materials]
[partial_pressure]
type = SorptionPartialPressure
A = 19.3
B = -47300
D = 1.51
E = 4340
d1 = 3.4
d2 = 6.15e-4
unit_scale = 1
density = 1
concentration = c
temperature = T
partial_pressure_base_name = test
[]
[]
[Executioner]
type = Steady
solve_type = NEWTON
[]
[Postprocessors]
[c_avg]
type = ElementAverageValue
variable = c
[]
[T_avg]
type = ElementAverageValue
variable = T
[]
[partial_pressure_avg]
type = ElementAverageMaterialProperty
mat_prop = test_partial_pressure
[]
[partial_pressure_dc_avg]
type = ElementAverageMaterialProperty
mat_prop = test_partial_pressure_dc
[]
[partial_pressure_avg_exact]
type = FunctionValuePostprocessor
function = partial_pressure_avg_exact
[]
[partial_pressure_dc_avg_exact]
type = FunctionValuePostprocessor
function = partial_pressure_dc_avg_exact
[]
[]
[Outputs]
csv = true
[]
(test/tests/sorption_interface/sorption_ad.i)
# This test is verifies the implementation of SorptionIsothermGapInterface with
# interface_configuration = solid_gas. It contains two solids separated by a gas.
# One solid acts as a source of u. u desorbs into the gap and then onto the
# other solid, which acts as a sink. Temperature increases over time to verify
# sorption behavior and species conservation. For SorptionIsothermGapInterface,
# the solid is always the element and the gas is always the neighbor.
# _________________________
# | source | gap | sink |
# --------------------------
[GlobalParams]
order = SECOND
[]
[Mesh]
coord_type = RZ
[gen]
type = GeneratedMeshGenerator
nx = 150
xmax = 1.5
dim = 1
elem_type = EDGE3
[]
[source_gap]
type = SubdomainBoundingBoxGenerator
block_id = 1
bottom_left = '0 0 0'
top_right = '0.75 0 0'
input = gen
[]
[gap_sink]
type = SubdomainBoundingBoxGenerator
block_id = 2
bottom_left = '0.75 0 0'
top_right = '1.5 0 0'
input = source_gap
[]
[source]
type = SubdomainBoundingBoxGenerator
block_id = 3
bottom_left = '0 0 0'
top_right = '0.5 0 0'
input = gap_sink
[]
[sink]
type = SubdomainBoundingBoxGenerator
block_id = 4
bottom_left = '1 0 0'
top_right = '1.5 0 0'
input = source
[]
[break]
type = BreakMeshByBlockGenerator
block_pairs = '1 3; 2 4'
split_interface = true
add_interface_on_two_sides = true
input = sink
[]
[]
[Variables]
[u]
[]
[]
[AuxVariables]
[temperature]
[]
[]
[Kernels]
[u_dt]
type = ADTimeDerivative
variable = u
[]
[u]
type = ADMatDiffusion
variable = u
diffusivity = diffusivity
[]
[]
[ICs]
[u_source]
type = FunctionIC
variable = u
function = 'if(x < 0.75, 1, 0)'
block = 3
[]
[]
[AuxKernels]
[temperature]
type = FunctionAux
variable = temperature
function = '(4000 - 2000 * x) / 2.0 * (t + 0.1)'
execute_on = 'INITIAL TIMESTEP_BEGIN LINEAR NONLINEAR TIMESTEP_END'
[]
[]
[InterfaceKernels]
[u_solid_gas]
type = ADSorptionIsothermGapInterface
variable = u
neighbor_var = u
interface_configuration = solid_gas
temperature = temperature
partial_pressure_name = partial_pressure
sorption_penalty = 1e10
boundary = 'Block3_Block1 Block4_Block2'
[]
[]
[Materials]
[properties_solid]
type = ADGenericConstantMaterial
prop_names = 'density diffusivity'
prop_values = '1 0.1'
block = '3 4'
[]
[properties_gas]
type = ADGenericConstantMaterial
prop_names = diffusivity
prop_values = 1
block = '1 2'
[]
[partial_pressure_source]
type = ADSorptionPartialPressure
A = 19.3
B = -47300
D = 1.51
E = 4340
d1 = 3.4
d2 = 6.15e-4
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = 3
outputs = 'all'
output_properties = partial_pressure
[]
[partial_pressure_sink]
type = ADSorptionPartialPressure
A = 24
B = -35700
D = -1.56
E = 6120
d1 = 2.04
d2 = 1.79e-3
unit_scale = 1
density = density
concentration = u
temperature = temperature
block = 4
outputs = 'all'
output_properties = partial_pressure
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
compute_scaling_once = false
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
line_search = none
nl_rel_tol = 1e-15
nl_abs_tol = 1e-15
l_tol = 1e-3
l_max_its = 50
nl_max_its = 50
dt = 0.1
end_time = 2
[]
[Functions]
[source_surface_u_exact]
type = ParsedFunction
symbol_names = 'P T R'
symbol_values = 'source_surface_pressure source_surface_temperature 8.31446261815324'
expression = '1.0 / R / T * P'
[]
[sink_surface_u_exact]
type = ParsedFunction
symbol_names = 'P T R'
symbol_values = 'sink_surface_pressure sink_surface_temperature 8.31446261815324'
expression = '1.0 / R / T * P'
[]
[]
[Postprocessors]
[source_surface_pressure]
type = ADSideAverageMaterialProperty
property = partial_pressure
boundary = Block3_Block1
[]
[source_surface_temperature]
type = SideAverageValue
variable = temperature
boundary = Block3_Block1
[]
[source_surface_u]
type = SideAverageValue
variable = u
boundary = Block1_Block3
[]
[source_surface_u_exact]
type = FunctionValuePostprocessor
function = source_surface_u_exact
[]
[sink_surface_pressure]
type = ADSideAverageMaterialProperty
property = partial_pressure
boundary = Block4_Block2
[]
[sink_surface_temperature]
type = SideAverageValue
variable = temperature
boundary = Block4_Block2
[]
[sink_surface_u]
type = SideAverageValue
variable = u
boundary = Block2_Block4
[]
[sink_surface_u_exact]
type = FunctionValuePostprocessor
function = sink_surface_u_exact
[]
[u_source_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = 3
[]
[u_gap_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = '1 2'
[]
[u_sink_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = 4
[]
[u_tot]
type = ElementIntegralVariablePostprocessor
variable = u
block = '1 2 3 4'
[]
[]
[Outputs]
exodus = true
csv = true
[]
[Debug]
show_var_residual_norms = true
[]