- eigenstrain_nameMaterial property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator.
C++ Type:std::string
Controllable:No
Description:Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator.
- stress_free_temperatureReference temperature at which there is no thermal expansion for thermal eigenstrain calculation
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation
GraphiteThermalExpansionEigenstrain
Calculates eigenstrain due to isotropic thermal expansion in AT 101 graphite
Description
This material computes the coefficient of thermal expansion, in 1/K, for AT 101 graphite based on the material properties given for a spark plasma sintering die in (Cincotti et al., 2007). The coefficient of thermal expansion value is computed solely as a function of temperature, T: (1) where the calibration range is 290.9K to 2383.0K.
The relationship for this thermal material property was set as a piecewise function with constant values based on the experimental data points (Cincotti et al., 2007), as shown in Figure 1.
Figure 1: Shown against the experimental data points, the coefficient of thermal expansion for AT 101 graphite is fit as a logarithmic function.
The local temperature is checked against the calibration range limits only once at the start of each timestep, using the value of the temperature from the previous converged timestep. The value of the coefficient of thermal expansion is calculated with the current value of the temperature.
Example Input File Syntax
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[graphite_thermal_expansion]
type = GraphiteThermalExpansionEigenstrain<<<{"description": "Calculates eigenstrain due to isotropic thermal expansion in AT 101 graphite", "href": "GraphiteThermalExpansionEigenstrain.html"}>>>
eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = graphite_thermal_expansion
stress_free_temperature<<<{"description": "Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation"}>>> = 300
temperature<<<{"description": "Coupled temperature"}>>> = temperature
[]
[]Input Parameters
- base_nameOptional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
C++ Type:std::string
Controllable:No
Description:Optional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- coeffient_thermal_expansion_scale_factor1The scaling factor for the coefficient of thermal expansion
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The scaling factor for the coefficient of thermal expansion
- 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.
- mean_thermal_expansion_coefficient_nameName of the mean coefficient of thermal expansion.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Name of the mean coefficient of thermal expansion.
- temperatureCoupled temperature
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled temperature
- use_old_temperatureFalseFlag to optionally use the temperature value from the previous timestep.
Default:False
C++ Type:bool
Controllable:No
Description:Flag to optionally use the temperature value from the previous timestep.
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
- (examples/sps/multiapp/diemodel/electrothermomechs/electrothermomechanical_plunger_powder_creep.i)
- (test/tests/materials/graphite/mechanical/graphite_thermal_expansion.i)
- (examples/sps/multiapp/diemodel/electrothermomechs/electrothermomechanical_plunger_powder_bothends_pressurebc.i)
- (examples/sps/multiapp/diemodel/electrothermomechs/mechanical_plunger_powder_plastic.i)
- (test/tests/materials/graphite/mechanical/ad_graphite_thermal_expansion.i)
- (examples/sps/multiapp/diemodel/electrothermomechs/simple_geometries/graphite_electrothermalmech_singleblock.i)
References
- A. Cincotti, A. M. Locci, R. Orrù, and G. Cao.
Modeling of SPS apparatus: temperature, current and strain distribution with no powders.
AIChE Journal, 53(3):703–719, 2007.
doi:10.1002/aic.11102.[BibTeX]
@article{cincotti2007sps, author = "Cincotti, A. and Locci, A. M. and Orrù, R. and Cao, G.", title = "Modeling of {SPS} apparatus: Temperature, current and strain distribution with no powders", journal = "AIChE Journal", volume = "53", number = "3", pages = "703-719", doi = "10.1002/aic.11102", year = "2007" }
(test/tests/materials/graphite/mechanical/graphite_thermal_expansion.i)
# This model includes only mechanics on a patch mesh and is intended
# verify the implementation of the piecewise function for the coefficient of
# thermal expansion from Cincotti et al (2007) AIChE Journal, Vol 53 No 3,
# page 710 Figure 7d. The piecewise function is given as a function within
# the input file (thexp_exact), and the value of that function is compared
# at each timestep against the computed strain_{xx} calculated value
#
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = ExamplePatchMeshGenerator
dim = 2
[]
coord_type = RZ
[]
[AuxVariables]
[temperature]
initial_condition = 300.0
[]
[]
[AuxKernels]
[temp_aux]
type = FunctionAux
variable = temperature
function = temp_ramp
[]
[]
[Functions]
[temp_ramp]
type = ParsedFunction
expression = '300 + 400.0/60.*t' #stand-in for a 400C/min heating rate
[]
[thermal_expansion]
type = ParsedFunction
symbol_names = T
symbol_values = temperature
expression = '(1.996e-6*log(T*4.799e-2)-4.041e-6)'
[]
[thexp_exact]
type = ParsedFunction
symbol_names = 'ref_temperature thermal_expansion temperature'
symbol_values = '300.0 thermal_expansion temperature'
expression = 'thermal_expansion * (temperature - ref_temperature)'
[]
[]
[Physics]
[SolidMechanics]
[QuasiStatic]
[all]
strain = SMALL
decomposition_method = EigenSolution
add_variables = true
eigenstrain_names = graphite_thermal_expansion
generate_output = 'strain_xx'
[]
[]
[]
[]
[BCs]
[left]
type = DirichletBC
variable = disp_x
value = 0
boundary = left
[]
[bottom]
type = DirichletBC
variable = disp_y
value = 0
boundary = bottom
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1
poissons_ratio = 0.3
[]
[graphite_thermal_expansion]
type = GraphiteThermalExpansionEigenstrain
eigenstrain_name = graphite_thermal_expansion
stress_free_temperature = 300
temperature = temperature
[]
[stress]
type = ComputeLinearElasticStress
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap'
petsc_options_value = 'asm ilu 1'
nl_max_its = 20
l_max_its = 50
dt = 25
dtmin = 1.0e-8
end_time = 300
[]
[Postprocessors]
[temperature]
type = ElementAverageValue
variable = temperature
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
[]
[thexp_exact]
type = FunctionValuePostprocessor
function = thexp_exact
[]
[thexp_diff]
type = DifferencePostprocessor
value1 = strain_xx
value2 = thexp_exact
outputs = none
[]
[thexp_max_diff]
type = TimeExtremeValue
postprocessor = thexp_diff
value_type = abs_max
[]
[]
[Outputs]
csv = true
perf_graph = true
[]
(examples/sps/multiapp/diemodel/electrothermomechs/electrothermomechanical_plunger_powder_creep.i)
## Units in the input file: m-Pa-s-K
initial_temperature = 300 #roughly 600C where the pyrometer kicks in
[GlobalParams]
displacements = 'disp_x disp_y'
volumetric_locking_correction = false
order = SECOND
[]
[Mesh]
file = stepped_plunger_powder_2d.e
coord_type = RZ
construct_side_list_from_node_list = true
patch_update_strategy = iteration
patch_size = 20
# ghosting_patch_size = 5*patch_size
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
[]
[Variables]
[temperature]
initial_condition = ${initial_temperature}
[]
[electric_potential]
[]
[]
[AuxVariables]
[vonmises_stress]
family = MONOMIAL
order = FIRST
[]
[specific_heat_capacity_va]
initial_condition = 842.2 # at 1500K #568.73 at 1000K #447.281 # at 293K
[]
[density_va]
initial_condition = 3106.2 ##5010.0*(1-${initial_porosity}) #in kg/m^3
[]
[heat_transfer_radiation]
[]
[thermal_conductivity]
family = MONOMIAL
order = FIRST
block = 'upper_plunger lower_plunger die_wall'
[]
[sigma_aeh]
initial_condition = 2.0e-10 #in units eV/((nV)^2-s-nm)
[]
# [electrical_conductivity]
# family = MONOMIAL
# order = FIRST
# block = 'upper_plunger lower_plunger die_wall'
# []
[microapp_potential]
#converted to microapp electronVolts units
[]
[E_x]
order = FIRST
family = MONOMIAL
[]
[E_y]
order = FIRST
family = MONOMIAL
[]
[yttria_current_density]
order = FIRST
family = MONOMIAL
[]
[graphite_current_density]
order = FIRST
family = MONOMIAL
[]
[current_density_J]
family = NEDELEC_ONE
order = FIRST
[]
[]
[Physics]
[SolidMechanics]
[QuasiStatic]
[graphite]
strain = FINITE
add_variables = true
use_automatic_differentiation = true
generate_output = 'strain_xx strain_xy strain_yy strain_zz stress_xx stress_xy stress_yy stress_zz'
extra_vector_tags = 'ref'
eigenstrain_names = 'graphite_thermal_expansion'
block = 'upper_plunger lower_plunger die_wall'
[]
[yttria]
strain = FINITE
add_variables = true
use_automatic_differentiation = true
generate_output = 'strain_xx strain_xy strain_yy strain_zz stress_xx stress_xy stress_yy stress_zz'
extra_vector_tags = 'ref'
eigenstrain_names = 'yttria_thermal_expansion'
block = 'powder_compact'
[]
[]
[]
[]
[Kernels]
[HeatDiff_graphite]
type = ADHeatConduction
variable = temperature
thermal_conductivity = thermal_conductivity
extra_vector_tags = 'ref'
block = 'upper_plunger lower_plunger die_wall'
[]
[HeatTdot_graphite]
type = ADHeatConductionTimeDerivative
variable = temperature
specific_heat = heat_capacity
density_name = graphite_density
extra_vector_tags = 'ref'
block = 'upper_plunger lower_plunger die_wall'
[]
[JouleHeating_graphite]
type = ADJouleHeatingSource
variable = temperature
elec = electric_potential
electrical_conductivity = electrical_conductivity
use_displaced_mesh = true
extra_vector_tags = 'ref'
block = 'upper_plunger lower_plunger die_wall'
[]
[electric_graphite]
type = ADMatDiffusion
variable = electric_potential
diffusivity = electrical_conductivity
extra_vector_tags = 'ref'
block = 'upper_plunger lower_plunger die_wall'
[]
[HeatDiff_yttria]
type = ADHeatConduction
variable = temperature
thermal_conductivity = yttria_thermal_conductivity #use parsed material property, hope it works
extra_vector_tags = 'ref'
block = powder_compact
[]
[HeatTdot_yttria]
type = ADHeatConductionTimeDerivative
variable = temperature
specific_heat = yttria_specific_heat_capacity #use parsed material property
density_name = yttria_density
extra_vector_tags = 'ref'
block = powder_compact
[]
[JouleHeating_yttria]
type = ADJouleHeatingSource
variable = temperature
elec = electric_potential
electrical_conductivity = electrical_conductivity
use_displaced_mesh = true
extra_vector_tags = 'ref'
block = powder_compact
[]
[electric_yttria]
type = ADMatDiffusion
variable = electric_potential
diffusivity = electrical_conductivity
use_displaced_mesh = true
extra_vector_tags = 'ref'
block = powder_compact
[]
[]
[AuxKernels]
[vonmises_stress]
type = ADRankTwoScalarAux
variable = vonmises_stress
rank_two_tensor = stress
scalar_type = VonMisesStress
[]
[heat_transfer_radiation]
type = ParsedAux
variable = heat_transfer_radiation
boundary = 'outer_die_wall'
coupled_variables = 'temperature'
constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
constant_expressions = '5.67e-8 0.85 300.0' #roughly room temperature, which is probably too cold
expression = '-boltzmann*epsilon*(temperature^4-temperature_farfield^4)'
[]
[thermal_conductivity]
type = ADMaterialRealAux
variable = thermal_conductivity
property = thermal_conductivity
block = 'upper_plunger lower_plunger die_wall'
[]
[microapp_potential]
type = ParsedAux
variable = microapp_potential
coupled_variables = electric_potential
expression = 'electric_potential*1e9' #convert from V to nV
block = 'powder_compact'
[]
[E_x]
type = VariableGradientComponent
variable = E_x
gradient_variable = electric_potential
component = x
# block = 'powder_compact'
[]
[E_y]
type = VariableGradientComponent
variable = E_y
gradient_variable = electric_potential
component = y
# block = 'powder_compact'
[]
[yttria_current_density]
type = ParsedAux
variable = yttria_current_density
coupled_variables = 'electrical_conductivity E_y'
expression = '-1.0*electrical_conductivity*E_y'
block = 'powder_compact'
[]
[graphite_current_density]
type = ParsedAux
variable = graphite_current_density
coupled_variables = 'electrical_conductivity E_y'
expression = '-1.0*electrical_conductivity*E_y'
block = 'upper_plunger lower_plunger die_wall'
[]
[current_density_J]
type = ADCurrentDensity
variable = current_density_J
potential = electric_potential
[]
[]
[BCs]
[centerline_no_dispx]
type = ADDirichletBC
preset = true
variable = disp_x
value = 0
boundary = 'centerline_powder_compact centerline_upper_plunger centerline_lower_plunger'
[]
[fixed_in_y]
type = ADDirichletBC
preset = true
variable = disp_y
value = 0
boundary = 'axial_centerpoint centerpoint_outer_die_wall' #centerpoint_inner_die_wall
[]
[bottom_pressure_ydirection]
type = ADPressure
variable = disp_y
boundary = 'bottom_lower_plunger'
function = 'if(t<1.0, (20.7e6/1.0)*t, 20.7e6)'
[]
[top_pressure_ydirection]
type = ADPressure
variable = disp_y
boundary = 'top_upper_plunger'
function = 'if(t<1.0, (20.7e6/1.0)*t, 20.7e6)'
[]
[external_surface]
type = CoupledVarNeumannBC
boundary = 'outer_die_wall'
variable = temperature
v = heat_transfer_radiation
[]
[electric_top]
type = ADFunctionDirichletBC
variable = electric_potential
boundary = 'top_upper_plunger'
function = 'if(t<1.0,0.0, if(t<21.0, 4.0e-2*t, 0.8))' #maximum value is 0.38V to maintain 10V/m drop
[]
[electric_bottom]
type = ADDirichletBC
variable = electric_potential
boundary = 'bottom_lower_plunger'
value = 0.0 #1.0 #0.0
[]
[]
[Functions]
[graphite_thermal_conductivity_fcn]
type = ParsedFunction
expression = '-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659'
[]
# [yttria_thermal_conductivity_fcn] #from the multiapp
# type = ParsedFunction
# expression = '3214.46/(t-147.73)'
# []
[harmonic_mean_thermal_conductivity]
type = ParsedFunction
expression = '2*(-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)*(3214.46/(t-147.73))/((-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)+(3214.46/(t-147.73)))'
# symbol_names = 'k_graphite k_yttria'
# symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
[]
[graphite_electrical_conductivity_fcn]
type = ParsedFunction
expression = '1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5)'
[]
# [electrical_conductivity_fcn]
# type = ParsedFunction
# # symbol_names = porosity
# # symbol_values = initial_porosity
# expression = '(1-0.62)*2.0025e4*exp(-1.61/8.617343e-5/t)'
# []
[harmonic_mean_electrical_conductivity]
type = ParsedFunction
expression = '2*(1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))*((1)*2.0025e4*exp(-1.61/8.617343e-5/t))/((1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))+((1)*2.0025e4*exp(-1.61/8.617343e-5/t)))'
# symbol_names = 'k_graphite k_yttria'
# symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
[]
[]
[Contact]
[upper_plunger_powder_mechanical]
primary = bottom_upper_plunger
secondary = top_powder_compact
formulation = kinematic
model = frictionless
penalty = 1e14 #1e4 is good in the mm-MPa system
normal_smoothing_distance = 0.1 #was 0.1 in mm-MPa system
normalize_penalty = true
[]
[upper_plunger_die_mechanical]
secondary = outer_upper_plunger
primary = inner_die_wall
formulation = kinematic
model = frictionless
penalty = 1e14 #isotropic yttria block takes 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[powder_die_mechanical]
primary = inner_die_wall
secondary = outer_powder_compact
formulation = kinematic
model = frictionless
penalty = 1e14 #isotropic yttria block takes 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[lower_plunger_die_mechanical]
primary = inner_die_wall
secondary = outer_lower_plunger
formulation = kinematic
model = frictionless
penalty = 1e14 #isotropic yttria block takes 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[powder_bottom_plunger_mechanical]
secondary = bottom_powder_compact
primary = top_lower_plunger
formulation = kinematic
model = frictionless
penalty = 1e14 #isotropic takes 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[]
[ThermalContact]
[lower_plunger_die_electric]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_lower_plunger
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'graphite_electrical_conductivity_fcn'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[top_plunger_die_electric]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_upper_plunger
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'graphite_electrical_conductivity_fcn'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[upper_plunger_powder_electric]
type = GapHeatTransfer
primary = bottom_upper_plunger
secondary = top_powder_compact
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85
gap_geometry_type = PLATE
gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[powder_die_electric]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_powder_compact
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[powder_bottom_plunger_electric]
type = GapHeatTransfer
primary = top_lower_plunger
secondary = bottom_powder_compact #expect more heat transfer from the die to the powder
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85 #assume no radiation heat transfer because no gap
gap_geometry_type = PLATE
gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[lower_plunger_die_thermal]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_lower_plunger
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[top_plunger_die_thermal]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_upper_plunger
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[upper_plunger_powder_thermal]
type = GapHeatTransfer
primary = bottom_upper_plunger
secondary = top_powder_compact
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85
gap_geometry_type = PLATE
gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[powder_die_thermal]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_powder_compact
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[powder_bottom_plunger_thermal]
type = GapHeatTransfer
primary = top_lower_plunger
secondary = bottom_powder_compact #expect more heat transfer from the die to the powder
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85 #assume no radiation heat transfer because no gap
gap_geometry_type = PLATE
gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[]
[Materials]
[graphite_elasticity_tensor]
type = ADComputeIsotropicElasticityTensor
youngs_modulus = 1.08e10 #in Pa, 10.8 GPa, Cincotti
poissons_ratio = 0.33
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_stress]
type = ADComputeFiniteStrainElasticStress
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_thermal_expansion]
type = ADGraphiteThermalExpansionEigenstrain
eigenstrain_name = graphite_thermal_expansion
stress_free_temperature = ${initial_temperature}
temperature = temperature
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_density]
type = ADGenericConstantMaterial
prop_names = 'graphite_density'
prop_values = 1.750e3 #in kg/m^3 from Cincotti et al 2007, Table 2, doi:10.1002/aic
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_thermal]
type = ADGraphiteThermal
temperature = temperature
block = 'upper_plunger lower_plunger die_wall'
[]
# [graphite_electrical]
# type = ADGraphiteElectricalConductivity
# temperature = temperature
# output_properties = all
# outputs = 'exodus'
# block = 'upper_plunger lower_plunger die_wall'
# []
# [graphite_electrical_resistivity]
# type = ADParsedMaterial
# property_name = electrical_resistivity
# material_property_names = electrical_conductivity
# expression = '1 / electrical_conductivity'
# output_properties = electrical_conductivity
# outputs = 'exodus'
# block = 'upper_plunger lower_plunger die_wall'
# []
[graphite_electrical_conductivity]
type = ADParsedMaterial
property_name = electrical_conductivity
coupled_variables = 'temperature'
expression = '1.0/(-2.705e-15*temperature^3+1.263e-11*temperature^2-1.836e-8*temperature+1.813e-5)'
output_properties = electrical_conductivity
outputs = 'exodus'
block = 'upper_plunger lower_plunger die_wall'
[]
[yttria_elasticity_tensor]
type = ADComputeIsotropicElasticityTensor
# youngs_modulus = 1.16e11 #in Pa for fully dense part
youngs_modulus = 1.38e10 #in Pa assuming 62% initial density and theoretical coeff. from Phani and Niyogi (1987)
poissons_ratio = 0.36
block = powder_compact
[]
[yttria_stress]
type = ADComputeMultipleInelasticStress
block = powder_compact
inelastic_models = yttria_creep_model
[]
[yttria_creep_model]
type = ADPowerLawCreepStressUpdate
block = powder_compact
coefficient = 3.75e-7 # from Al's work
n_exponent = 0.714 # from Al's work
activation_energy = 0.0
use_substepping = INCREMENT_BASED
absolute_tolerance = 5e-9
[]
[yttria_thermal_expansion]
type = ADComputeThermalExpansionEigenstrain
thermal_expansion_coeff = 9.3e-6 # from https://doi.org/10.1111/j.1151-2916.1957.tb12619.x
eigenstrain_name = yttria_thermal_expansion
stress_free_temperature = 300
temperature = temperature
[]
[yttria_thermal_conductivity]
type = ADParsedMaterial
coupled_variables = 'temperature'
expression = '3214.46 / (temperature - 147.73)' #in W/(m-K) #Given from Larry's curve fitting, data from Klein and Croft, JAP, v. 38, p. 1603 and UC report "For Computer Heat Conduction Calculations - A compilation of thermal properties data" by A.L. Edwards, UCRL-50589 (1969)
# coupled_variables = 'thermal_conductivity_aeh'
# expression = 'thermal_conductivity_aeh' #in W/(m-K) directly, for now
property_name = 'yttria_thermal_conductivity'
output_properties = yttria_thermal_conductivity
outputs = 'exodus'
block = powder_compact
[]
[yttria_specific_heat_capacity]
type = ADParsedMaterial
property_name = yttria_specific_heat_capacity
coupled_variables = 'specific_heat_capacity_va'
expression = 'specific_heat_capacity_va' #in J/(K-kg)
# output_properties = yttria_specific_heat_capacity
# outputs = 'exodus'
block = powder_compact
[]
[yttria_density]
type = ADParsedMaterial
property_name = 'yttria_density'
coupled_variables = 'density_va'
expression = 'density_va'
# output_properties = yttria_density
# outputs = 'exodus'
block = powder_compact
[]
[electrical_conductivity]
type = ADParsedMaterial
# coupled_variables = 'sigma_aeh'
# expression = 'sigma_aeh*1.602e8' #converts to units of J/(V^2-m-s)
property_name = 'electrical_conductivity'
output_properties = electrical_conductivity
outputs = 'exodus'
block = powder_compact
# type = ADDerivativeParsedMaterial
# property_name = electrical_conductivity
coupled_variables = 'temperature'
constant_names = 'Q_elec kB prefactor_solid initial_porosity'
constant_expressions = '1.61 8.617343e-5 1.25e-4 0.38'
expression = '(1-initial_porosity) * prefactor_solid * exp(-Q_elec/kB/temperature) * 1.602e8' # in eV/(nV^2 s nm) per chat with Larry, last term converts to units of J/(V^2-m-s)
[]
[]
[Dampers]
[disp_x_damp]
type = ElementJacobianDamper
# max_increment = 2.0e-4
use_displaced_mesh = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
line_search = 'none'
# compute_scaling_once = false
# force running options
# petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type -pc_factor_shift_amount'
# petsc_options_value = 'lu basic NONZERO 1e-15'
# #mechanical contact options
# petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu superlu_dist'
petsc_options = '-snes_converged_reason -ksp_converged_reason'
nl_forced_its = 1
nl_rel_tol = 2e-5 #was 1e-10, for temperature only
nl_abs_tol = 2e-12 #was 1e-10, before that 1e-12
nl_max_its = 20
l_max_its = 50
dtmin = 1.0e-4
end_time = 600 #900 #15 minutes, rule of thumb from Dennis is 10 minutes
[Quadrature]
order = FIFTH
side_order = SEVENTH
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = 0.05
optimal_iterations = 8
iteration_window = 2
[]
[]
[Postprocessors]
[stress_xx]
type = ElementAverageValue
variable = stress_xx
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
[]
[stress_yy]
type = ElementAverageValue
variable = stress_yy
[]
[strain_yy]
type = ElementAverageValue
variable = strain_yy
[]
[vonmises_stress]
type = ElementAverageValue
variable = vonmises_stress
[]
[temperature_pp]
type = AverageNodalVariableValue
variable = temperature
[]
[graphite_thermal_conductivity]
type = ElementAverageValue
variable = thermal_conductivity
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_sigma]
type = ElementAverageValue
variable = electrical_conductivity
block = 'upper_plunger lower_plunger die_wall'
[]
[yttria_thermal_conductivity]
type = ElementAverageValue
variable = yttria_thermal_conductivity
block = powder_compact
[]
[yttria_sigma]
type = ElementAverageValue
variable = electrical_conductivity
block = powder_compact
[]
[]
[Outputs]
csv = true
exodus = true
perf_graph = true
[ckpt]
type = Checkpoint
time_step_interval = 1
num_files = 2
[]
[]
(test/tests/materials/graphite/mechanical/graphite_thermal_expansion.i)
# This model includes only mechanics on a patch mesh and is intended
# verify the implementation of the piecewise function for the coefficient of
# thermal expansion from Cincotti et al (2007) AIChE Journal, Vol 53 No 3,
# page 710 Figure 7d. The piecewise function is given as a function within
# the input file (thexp_exact), and the value of that function is compared
# at each timestep against the computed strain_{xx} calculated value
#
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = ExamplePatchMeshGenerator
dim = 2
[]
coord_type = RZ
[]
[AuxVariables]
[temperature]
initial_condition = 300.0
[]
[]
[AuxKernels]
[temp_aux]
type = FunctionAux
variable = temperature
function = temp_ramp
[]
[]
[Functions]
[temp_ramp]
type = ParsedFunction
expression = '300 + 400.0/60.*t' #stand-in for a 400C/min heating rate
[]
[thermal_expansion]
type = ParsedFunction
symbol_names = T
symbol_values = temperature
expression = '(1.996e-6*log(T*4.799e-2)-4.041e-6)'
[]
[thexp_exact]
type = ParsedFunction
symbol_names = 'ref_temperature thermal_expansion temperature'
symbol_values = '300.0 thermal_expansion temperature'
expression = 'thermal_expansion * (temperature - ref_temperature)'
[]
[]
[Physics]
[SolidMechanics]
[QuasiStatic]
[all]
strain = SMALL
decomposition_method = EigenSolution
add_variables = true
eigenstrain_names = graphite_thermal_expansion
generate_output = 'strain_xx'
[]
[]
[]
[]
[BCs]
[left]
type = DirichletBC
variable = disp_x
value = 0
boundary = left
[]
[bottom]
type = DirichletBC
variable = disp_y
value = 0
boundary = bottom
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 1
poissons_ratio = 0.3
[]
[graphite_thermal_expansion]
type = GraphiteThermalExpansionEigenstrain
eigenstrain_name = graphite_thermal_expansion
stress_free_temperature = 300
temperature = temperature
[]
[stress]
type = ComputeLinearElasticStress
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap'
petsc_options_value = 'asm ilu 1'
nl_max_its = 20
l_max_its = 50
dt = 25
dtmin = 1.0e-8
end_time = 300
[]
[Postprocessors]
[temperature]
type = ElementAverageValue
variable = temperature
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
[]
[thexp_exact]
type = FunctionValuePostprocessor
function = thexp_exact
[]
[thexp_diff]
type = DifferencePostprocessor
value1 = strain_xx
value2 = thexp_exact
outputs = none
[]
[thexp_max_diff]
type = TimeExtremeValue
postprocessor = thexp_diff
value_type = abs_max
[]
[]
[Outputs]
csv = true
perf_graph = true
[]
(examples/sps/multiapp/diemodel/electrothermomechs/electrothermomechanical_plunger_powder_bothends_pressurebc.i)
## Units in the input file: m-Pa-s-K
initial_temperature = 300 #roughly 600C where the pyrometer kicks in
[GlobalParams]
displacements = 'disp_x disp_y'
volumetric_locking_correction = false
order = SECOND
[]
[Mesh]
file = stepped_plunger_powder_2d.e
coord_type = RZ
construct_side_list_from_node_list = true
patch_update_strategy = iteration
patch_size = 20
# ghosting_patch_size = 5*patch_size
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
[]
[Variables]
[temperature]
initial_condition = ${initial_temperature}
[]
[electric_potential]
[]
[]
[AuxVariables]
[vonmises_stress]
family = MONOMIAL
order = FIRST
[]
[specific_heat_capacity_va]
initial_condition = 842.2 # at 1500K #568.73 at 1000K #447.281 # at 293K
[]
[density_va]
initial_condition = 3106.2 ##5010.0*(1-${initial_porosity}) #in kg/m^3
[]
[heat_transfer_radiation]
[]
[thermal_conductivity]
family = MONOMIAL
order = FIRST
block = 'upper_plunger lower_plunger die_wall'
[]
[sigma_aeh]
initial_condition = 2.0e-10 #in units eV/((nV)^2-s-nm)
[]
# [electrical_conductivity]
# family = MONOMIAL
# order = FIRST
# block = 'upper_plunger lower_plunger die_wall'
# []
[microapp_potential]
#converted to microapp electronVolts units
[]
[E_x]
order = FIRST
family = MONOMIAL
[]
[E_y]
order = FIRST
family = MONOMIAL
[]
[yttria_current_density]
order = FIRST
family = MONOMIAL
[]
[graphite_current_density]
order = FIRST
family = MONOMIAL
[]
[current_density_J]
family = NEDELEC_ONE
order = FIRST
[]
[]
[Physics]
[SolidMechanics]
[QuasiStatic]
[graphite]
strain = FINITE
add_variables = true
use_automatic_differentiation = true
generate_output = 'strain_xx strain_xy strain_yy strain_zz stress_xx stress_xy stress_yy stress_zz'
extra_vector_tags = 'ref'
eigenstrain_names = 'graphite_thermal_expansion'
block = 'upper_plunger lower_plunger die_wall'
[]
[yttria]
strain = FINITE
add_variables = true
use_automatic_differentiation = true
generate_output = 'strain_xx strain_xy strain_yy strain_zz stress_xx stress_xy stress_yy stress_zz'
extra_vector_tags = 'ref'
block = 'powder_compact'
[]
[]
[]
[]
[Kernels]
[HeatDiff_graphite]
type = ADHeatConduction
variable = temperature
thermal_conductivity = thermal_conductivity
extra_vector_tags = 'ref'
block = 'upper_plunger lower_plunger die_wall'
[]
[HeatTdot_graphite]
type = ADHeatConductionTimeDerivative
variable = temperature
specific_heat = heat_capacity
density_name = graphite_density
extra_vector_tags = 'ref'
block = 'upper_plunger lower_plunger die_wall'
[]
[JouleHeating_graphite]
type = ADJouleHeatingSource
variable = temperature
elec = electric_potential
electrical_conductivity = electrical_conductivity
use_displaced_mesh = true
extra_vector_tags = 'ref'
block = 'upper_plunger lower_plunger die_wall'
[]
[electric_graphite]
type = ADMatDiffusion
variable = electric_potential
diffusivity = electrical_conductivity
extra_vector_tags = 'ref'
block = 'upper_plunger lower_plunger die_wall'
[]
[HeatDiff_yttria]
type = ADHeatConduction
variable = temperature
thermal_conductivity = yttria_thermal_conductivity #use parsed material property, hope it works
extra_vector_tags = 'ref'
block = powder_compact
[]
[HeatTdot_yttria]
type = ADHeatConductionTimeDerivative
variable = temperature
specific_heat = yttria_specific_heat_capacity #use parsed material property
density_name = yttria_density
extra_vector_tags = 'ref'
block = powder_compact
[]
[JouleHeating_yttria]
type = ADJouleHeatingSource
variable = temperature
elec = electric_potential
electrical_conductivity = electrical_conductivity
use_displaced_mesh = true
extra_vector_tags = 'ref'
block = powder_compact
[]
[electric_yttria]
type = ADMatDiffusion
variable = electric_potential
diffusivity = electrical_conductivity
use_displaced_mesh = true
extra_vector_tags = 'ref'
block = powder_compact
[]
[]
[AuxKernels]
[vonmises_stress]
type = ADRankTwoScalarAux
variable = vonmises_stress
rank_two_tensor = stress
scalar_type = VonMisesStress
[]
[heat_transfer_radiation]
type = ParsedAux
variable = heat_transfer_radiation
boundary = 'outer_die_wall'
coupled_variables = 'temperature'
constant_names = 'boltzmann epsilon temperature_farfield' #published emissivity for graphite is 0.85
constant_expressions = '5.67e-8 0.85 300.0' #roughly room temperature, which is probably too cold
expression = '-boltzmann*epsilon*(temperature^4-temperature_farfield^4)'
[]
[thermal_conductivity]
type = ADMaterialRealAux
variable = thermal_conductivity
property = thermal_conductivity
block = 'upper_plunger lower_plunger die_wall'
[]
[microapp_potential]
type = ParsedAux
variable = microapp_potential
coupled_variables = electric_potential
expression = 'electric_potential*1e9' #convert from V to nV
block = 'powder_compact'
[]
[E_x]
type = VariableGradientComponent
variable = E_x
gradient_variable = electric_potential
component = x
# block = 'powder_compact'
[]
[E_y]
type = VariableGradientComponent
variable = E_y
gradient_variable = electric_potential
component = y
# block = 'powder_compact'
[]
[yttria_current_density]
type = ParsedAux
variable = yttria_current_density
coupled_variables = 'electrical_conductivity E_y'
expression = '-1.0*electrical_conductivity*E_y'
block = 'powder_compact'
[]
[graphite_current_density]
type = ParsedAux
variable = graphite_current_density
coupled_variables = 'electrical_conductivity E_y'
expression = '-1.0*electrical_conductivity*E_y'
block = 'upper_plunger lower_plunger die_wall'
[]
[current_density_J]
type = ADCurrentDensity
variable = current_density_J
potential = electric_potential
[]
[]
[BCs]
[centerline_no_dispx]
type = ADDirichletBC
preset = true
variable = disp_x
value = 0
boundary = 'centerline_powder_compact centerline_upper_plunger centerline_lower_plunger'
[]
[fixed_in_y]
type = ADDirichletBC
preset = true
variable = disp_y
value = 0
boundary = 'axial_centerpoint'
[]
[bottom_pressure_ydirection]
type = ADPressure
variable = disp_y
boundary = 'bottom_lower_plunger'
function = 'if(t<1.0, (20.7e6/1.0)*t, 20.7e6)'
[]
[top_pressure_ydirection]
type = ADPressure
variable = disp_y
boundary = 'top_upper_plunger'
function = 'if(t<1.0, (20.7e6/1.0)*t, 20.7e6)'
[]
[external_surface]
type = CoupledVarNeumannBC
boundary = 'outer_die_wall'
variable = temperature
v = heat_transfer_radiation
[]
[electric_top]
type = ADFunctionDirichletBC
variable = electric_potential
boundary = 'top_upper_plunger'
function = 'if(t<1.0,0.0, if(t<21.0, 4.0e-2*t, 0.8))' #maximum value is 0.38V to maintain 10V/m drop
[]
[electric_bottom]
type = ADDirichletBC
variable = electric_potential
boundary = 'bottom_lower_plunger'
value = 0.0 #1.0 #0.0
[]
[]
[Functions]
[graphite_thermal_conductivity_fcn]
type = ParsedFunction
expression = '-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659'
[]
# [yttria_thermal_conductivity_fcn] #from the multiapp
# type = ParsedFunction
# expression = '3214.46/(t-147.73)'
# []
[harmonic_mean_thermal_conductivity]
type = ParsedFunction
expression = '2*(-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)*(3214.46/(t-147.73))/((-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)+(3214.46/(t-147.73)))'
# symbol_names = 'k_graphite k_yttria'
# symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
[]
[graphite_electrical_conductivity_fcn]
type = ParsedFunction
expression = '1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5)'
[]
# [electrical_conductivity_fcn]
# type = ParsedFunction
# # symbol_names = porosity
# # symbol_values = initial_porosity
# expression = '(1-0.62)*2.0025e4*exp(-1.61/8.617343e-5/t)'
# []
[harmonic_mean_electrical_conductivity]
type = ParsedFunction
expression = '2*(1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))*((1)*2.0025e4*exp(-1.61/8.617343e-5/t))/((1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))+((1)*2.0025e4*exp(-1.61/8.617343e-5/t)))'
# symbol_names = 'k_graphite k_yttria'
# symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
[]
[]
[Contact]
[upper_plunger_powder_mechanical]
primary = bottom_upper_plunger
secondary = top_powder_compact
formulation = kinematic
model = frictionless
penalty = 1e12 #1e4 is good in the mm-MPa system
normal_smoothing_distance = 0.1 #was 0.1 in mm-MPa system
normalize_penalty = true
[]
[upper_plunger_die_mechanical]
secondary = outer_upper_plunger
primary = inner_die_wall
formulation = kinematic
model = frictionless
penalty = 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[powder_die_mechanical]
primary = inner_die_wall
secondary = outer_powder_compact
formulation = kinematic
model = frictionless
penalty = 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[lower_plunger_die_mechanical]
primary = inner_die_wall
secondary = outer_lower_plunger
formulation = kinematic
model = frictionless
penalty = 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[powder_bottom_plunger_mechanical]
secondary = bottom_powder_compact
primary = top_lower_plunger
formulation = kinematic
model = frictionless
penalty = 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[]
[ThermalContact]
[lower_plunger_die_electric]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_lower_plunger
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'graphite_electrical_conductivity_fcn'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[top_plunger_die_electric]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_upper_plunger
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'graphite_electrical_conductivity_fcn'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[upper_plunger_powder_electric]
type = GapHeatTransfer
primary = bottom_upper_plunger
secondary = top_powder_compact
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85
gap_geometry_type = PLATE
gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[powder_die_electric]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_powder_compact
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[powder_bottom_plunger_electric]
type = GapHeatTransfer
primary = top_lower_plunger
secondary = bottom_powder_compact #expect more heat transfer from the die to the powder
variable = electric_potential
quadrature = true
emissivity_primary = 0.0 #0.85 #cincotti 2007, table 2
emissivity_secondary = 0.0 #0.85 #assume no radiation heat transfer because no gap
gap_geometry_type = PLATE
gap_conductivity_function = 'harmonic_mean_electrical_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[lower_plunger_die_thermal]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_lower_plunger
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[top_plunger_die_thermal]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_upper_plunger
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'graphite_thermal_conductivity_fcn'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[upper_plunger_powder_thermal]
type = GapHeatTransfer
primary = bottom_upper_plunger
secondary = top_powder_compact
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85
gap_geometry_type = PLATE
gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[powder_die_thermal]
type = GapHeatTransfer
primary = inner_die_wall
secondary = outer_powder_compact
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85
# gap_geometry_type = PLATE # Not for vertical surfaces
gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[powder_bottom_plunger_thermal]
type = GapHeatTransfer
primary = top_lower_plunger
secondary = bottom_powder_compact #expect more heat transfer from the die to the powder
variable = temperature
quadrature = true
emissivity_primary = 0.85 #cincotti 2007, table 2
emissivity_secondary = 0.85 #assume no radiation heat transfer because no gap
gap_geometry_type = PLATE
gap_conductivity_function = 'harmonic_mean_thermal_conductivity'
gap_conductivity_function_variable = temperature
normal_smoothing_distance = 0.1
[]
[]
[Materials]
[graphite_elasticity_tensor]
type = ADComputeIsotropicElasticityTensor
youngs_modulus = 1.08e10 #in Pa, 10.8 GPa, Cincotti
poissons_ratio = 0.33
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_stress]
type = ADComputeFiniteStrainElasticStress
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_thermal_expansion]
type = ADGraphiteThermalExpansionEigenstrain
eigenstrain_name = graphite_thermal_expansion
stress_free_temperature = ${initial_temperature}
temperature = temperature
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_density]
type = ADGenericConstantMaterial
prop_names = 'graphite_density'
prop_values = 1.750e3 #in kg/m^3 from Cincotti et al 2007, Table 2, doi:10.1002/aic
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_thermal]
type = ADGraphiteThermal
temperature = temperature
block = 'upper_plunger lower_plunger die_wall'
[]
# [graphite_electrical]
# type = ADGraphiteElectricalConductivity
# temperature = temperature
# output_properties = all
# outputs = 'exodus'
# block = 'upper_plunger lower_plunger die_wall'
# []
# [graphite_electrical_resistivity]
# type = ADParsedMaterial
# property_name = electrical_resistivity
# material_property_names = electrical_conductivity
# expression = '1 / electrical_conductivity'
# output_properties = electrical_conductivity
# outputs = 'exodus'
# block = 'upper_plunger lower_plunger die_wall'
# []
[graphite_electrical_conductivity]
type = ADParsedMaterial
property_name = electrical_conductivity
coupled_variables = 'temperature'
expression = '1.0/(-2.705e-15*temperature^3+1.263e-11*temperature^2-1.836e-8*temperature+1.813e-5)'
output_properties = electrical_conductivity
outputs = 'exodus'
block = 'upper_plunger lower_plunger die_wall'
[]
[yttria_elasticity_tensor]
type = ADComputeIsotropicElasticityTensor
# youngs_modulus = 1.16e5 #in MPa, 116 GPa, de Jong et al (2015) Scientific Data 2
youngs_modulus = 1.38e10 #in Pa assuming 62% initial density and theoretical coeff. from Phani and Niyogi (1987)
# youngs_modulus = 1.38e9
poissons_ratio = 0.36
block = powder_compact
[]
[yttria_stress]
type = ADComputeFiniteStrainElasticStress
block = powder_compact
[]
[yttria_thermal_conductivity]
type = ADParsedMaterial
coupled_variables = 'temperature'
expression = '3214.46 / (temperature - 147.73)' #in W/(m-K) #Given from Larry's curve fitting, data from Klein and Croft, JAP, v. 38, p. 1603 and UC report "For Computer Heat Conduction Calculations - A compilation of thermal properties data" by A.L. Edwards, UCRL-50589 (1969)
# coupled_variables = 'thermal_conductivity_aeh'
# expression = 'thermal_conductivity_aeh' #in W/(m-K) directly, for now
property_name = 'yttria_thermal_conductivity'
output_properties = yttria_thermal_conductivity
outputs = 'exodus'
block = powder_compact
[]
[yttria_specific_heat_capacity]
type = ADParsedMaterial
property_name = yttria_specific_heat_capacity
coupled_variables = 'specific_heat_capacity_va'
expression = 'specific_heat_capacity_va' #in J/(K-kg)
# output_properties = yttria_specific_heat_capacity
# outputs = 'exodus'
block = powder_compact
[]
[yttria_density]
type = ADParsedMaterial
property_name = 'yttria_density'
coupled_variables = 'density_va'
expression = 'density_va'
# output_properties = yttria_density
# outputs = 'exodus'
block = powder_compact
[]
[electrical_conductivity]
type = ADParsedMaterial
# coupled_variables = 'sigma_aeh'
# expression = 'sigma_aeh*1.602e8' #converts to units of J/(V^2-m-s)
property_name = 'electrical_conductivity'
output_properties = electrical_conductivity
outputs = 'exodus'
block = powder_compact
# type = ADDerivativeParsedMaterial
# property_name = electrical_conductivity
coupled_variables = 'temperature'
constant_names = 'Q_elec kB prefactor_solid initial_porosity'
constant_expressions = '1.61 8.617343e-5 1.25e-4 0.38'
expression = '(1-initial_porosity) * prefactor_solid * exp(-Q_elec/kB/temperature) * 1.602e8' # in eV/(nV^2 s nm) per chat with Larry, last term converts to units of J/(V^2-m-s)
[]
[]
[Dampers]
[disp_x_damp]
type = ElementJacobianDamper
# max_increment = 2.0e-4
use_displaced_mesh = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
line_search = 'none'
# compute_scaling_once = false
# force running options
# petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type -pc_factor_shift_amount'
# petsc_options_value = 'lu basic NONZERO 1e-15'
# #mechanical contact options
# petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu superlu_dist'
petsc_options = '-snes_converged_reason -ksp_converged_reason'
nl_forced_its = 1
nl_rel_tol = 2e-5 #was 1e-10, for temperature only
nl_abs_tol = 2e-12 #was 1e-10, before that 1e-12
nl_max_its = 20
l_max_its = 50
dtmin = 1.0e-3
end_time = 600 #900 #15 minutes, rule of thumb from Dennis is 10 minutes
[Quadrature]
order = FIFTH
side_order = SEVENTH
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = 0.05
optimal_iterations = 8
iteration_window = 2
[]
[]
[Postprocessors]
[stress_xx]
type = ElementAverageValue
variable = stress_xx
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
[]
[stress_yy]
type = ElementAverageValue
variable = stress_yy
[]
[strain_yy]
type = ElementAverageValue
variable = strain_yy
[]
[vonmises_stress]
type = ElementAverageValue
variable = vonmises_stress
[]
[temperature_pp]
type = AverageNodalVariableValue
variable = temperature
[]
[graphite_thermal_conductivity]
type = ElementAverageValue
variable = thermal_conductivity
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_sigma]
type = ElementAverageValue
variable = electrical_conductivity
block = 'upper_plunger lower_plunger die_wall'
[]
[yttria_thermal_conductivity]
type = ElementAverageValue
variable = yttria_thermal_conductivity
block = powder_compact
[]
[yttria_sigma]
type = ElementAverageValue
variable = electrical_conductivity
block = powder_compact
[]
[]
[Outputs]
csv = true
exodus = true
perf_graph = true
[ckpt]
type = Checkpoint
time_step_interval = 1
num_files = 2
[]
[]
(examples/sps/multiapp/diemodel/electrothermomechs/mechanical_plunger_powder_plastic.i)
## Units in the input file: m-Pa-s-K
initial_temperature = 300 #roughly 600C where the pyrometer kicks in
[GlobalParams]
displacements = 'disp_x disp_y'
volumetric_locking_correction = false
order = SECOND
[]
[Mesh]
[fmg]
type = FileMeshGenerator
file = stepped_plunger_powder_2d.e
[]
coord_type = RZ
construct_side_list_from_node_list = true
patch_update_strategy = iteration
patch_size = 20
# ghosting_patch_size = 5*patch_size
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
group_variables = 'disp_x disp_y'
[]
# [Variables]
# [temperature]
# initial_condition = ${initial_temperature}
# []
# []
[AuxVariables]
[vonmises_stress]
family = MONOMIAL
order = FIRST
[]
[density_va]
initial_condition = 3106.2 ##5010.0*(1-${initial_porosity}) #in kg/m^3
[]
[temperature]
order = FIRST
family = LAGRANGE
initial_condition = ${initial_temperature}
[]
[]
[Physics]
[SolidMechanics]
[QuasiStatic]
[graphite]
strain = FINITE
add_variables = true
use_automatic_differentiation = true
generate_output = 'strain_xx strain_xy strain_yy strain_zz stress_xx stress_xy stress_yy stress_zz'
extra_vector_tags = 'ref'
eigenstrain_names = 'graphite_thermal_expansion'
block = 'upper_plunger lower_plunger die_wall'
[]
[yttria]
strain = FINITE
add_variables = true
use_automatic_differentiation = true
generate_output = 'plastic_strain_yy strain_xx strain_xy strain_yy strain_zz stress_xx stress_xy stress_yy stress_zz'
extra_vector_tags = 'ref'
eigenstrain_names = 'yttria_thermal_expansion'
block = 'powder_compact'
[]
[]
[]
[]
[AuxKernels]
[vonmises_stress]
type = ADRankTwoScalarAux
variable = vonmises_stress
rank_two_tensor = stress
scalar_type = VonMisesStress
[]
[temp_aux]
type = FunctionAux
variable = temperature
function = temp_hist
[]
[]
[BCs]
[centerline_no_dispx]
type = ADDirichletBC
preset = true
variable = disp_x
value = 0
boundary = 'centerline_powder_compact centerline_upper_plunger centerline_lower_plunger'
[]
[die_no_disp]
type = ADDirichletBC
preset = true
variable = disp_y
value = 0
boundary = 'centerpoint_outer_die_wall'
[]
[bottom_disp]
type = ADFunctionDirichletBC
variable = disp_y
boundary = 'bottom_lower_plunger'
function = 'if(t<5.0, (125.0e-6/5.0)*t, 125.0e-6)'
[]
[top_disp]
type = ADFunctionDirichletBC
variable = disp_y
boundary = 'top_upper_plunger'
function = 'if(t<5.0, (-125.0e-6/5.0)*t, -125.0e-6)'
[]
[]
[Functions]
[yield]
type = PiecewiseLinear
x = '100 600 700 800 900 5000' #temperature
y = '15e6 15e6 14e6 13e6 10e6 10e6' #yield stress
[]
[temp_hist]
type = PiecewiseLinear
x = '0 5 10' #time
y = '300 600 900' #temperature
[]
[]
[Contact]
[upper_plunger_powder_mechanical]
primary = bottom_upper_plunger
secondary = top_powder_compact
formulation = penalty
model = frictionless
penalty = 1e14 #1e4 is good in the mm-MPa system
normal_smoothing_distance = 0.1 #was 0.1 in mm-MPa system
normalize_penalty = true
[]
[upper_plunger_die_mechanical]
secondary = outer_upper_plunger
primary = inner_die_wall
formulation = penalty
model = frictionless
penalty = 1e14 #isotropic yttria block takes 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[powder_die_mechanical]
primary = inner_die_wall
secondary = outer_powder_compact
formulation = penalty
model = frictionless
penalty = 1e14 #isotropic yttria block takes 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[lower_plunger_die_mechanical]
primary = inner_die_wall
secondary = outer_lower_plunger
formulation = penalty
model = frictionless
penalty = 1e14 #isotropic yttria block takes 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[powder_bottom_plunger_mechanical]
secondary = bottom_powder_compact
primary = top_lower_plunger
formulation = penalty
model = frictionless
penalty = 1e14 #isotropic takes 1e12
normal_smoothing_distance = 0.1
normalize_penalty = true
[]
[]
[Materials]
[graphite_elasticity_tensor]
type = ADComputeIsotropicElasticityTensor
youngs_modulus = 1.08e10 #in Pa, 10.8 GPa, Cincotti
poissons_ratio = 0.33
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_stress]
type = ADComputeFiniteStrainElasticStress
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_thermal_expansion]
type = ADGraphiteThermalExpansionEigenstrain
eigenstrain_name = graphite_thermal_expansion
stress_free_temperature = 300
temperature = temperature
block = 'upper_plunger lower_plunger die_wall'
[]
[graphite_density]
type = ADGenericConstantMaterial
prop_names = 'graphite_density'
prop_values = 1.750e3 #in kg/m^3 from Cincotti et al 2007, Table 2, doi:10.1002/aic
block = 'upper_plunger lower_plunger die_wall'
[]
[yttria_elasticity_tensor]
type = ADComputeIsotropicElasticityTensor
# youngs_modulus = 1.16e11 #in Pa for fully dense part
youngs_modulus = 1.38e10 #in Pa assuming 62% initial density and theoretical coeff. from Phani and Niyogi (1987)
poissons_ratio = 0.36
block = powder_compact
[]
[yttria_stress]
type = ADComputeMultipleInelasticStress
block = powder_compact
inelastic_models = yttria_plastic_model
[]
[yttria_plastic_model]
type = ADIsotropicPlasticityStressUpdate
block = powder_compact
hardening_constant = 0.1e10
# yield_stress = 15e6
yield_stress_function = yield
temperature = temperature
outputs = all
# relative_tolerance = 1e-20
# absolute_tolerance = 1e-8
# max_inelastic_increment = 0.000001
[]
[yttria_thermal_expansion]
type = ADComputeThermalExpansionEigenstrain
thermal_expansion_coeff = 9.3e-6 # from https://doi.org/10.1111/j.1151-2916.1957.tb12619.x
eigenstrain_name = yttria_thermal_expansion
stress_free_temperature = 300
temperature = temperature
[]
[yttria_density]
type = ADParsedMaterial
property_name = 'yttria_density'
coupled_variables = 'density_va'
expression = 'density_va'
# output_properties = yttria_density
# outputs = 'exodus'
block = powder_compact
[]
[]
[Dampers]
[disp_x_damp]
type = ElementJacobianDamper
# max_increment = 2.0e-4
use_displaced_mesh = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
line_search = 'none'
# compute_scaling_once = false
# force running options
# petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type -pc_factor_shift_amount'
# petsc_options_value = 'lu basic NONZERO 1e-15'
# #mechanical contact options
# petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu superlu_dist'
petsc_options = '-snes_converged_reason -ksp_converged_reason'
nl_forced_its = 1
nl_rel_tol = 2e-5 #was 1e-10, for temperature only
nl_abs_tol = 2e-12 #was 1e-10, before that 1e-12
nl_max_its = 20
l_max_its = 50
dtmin = 1.0e-4
dtmax = .4
dt = 0.4
end_time = 10 #900 #15 minutes, rule of thumb from Dennis is 10 minutes
[Quadrature]
order = FIFTH
side_order = SEVENTH
[]
# [TimeStepper]
# type = IterationAdaptiveDT
# dt = 0.1
# optimal_iterations = 8
# iteration_window = 2
# []
[]
[Postprocessors]
[powder_strain_yy]
type = ElementAverageValue
variable = strain_yy
block = powder_compact
[]
[powder_stress_yy]
type = ElementAverageValue
variable = stress_yy
block = powder_compact
[]
[powder_vonmises_stress]
type = ElementAverageValue
variable = vonmises_stress
block = powder_compact
[]
[powder_plastic_strain_yy]
type = ElementAverageValue
variable = plastic_strain_yy
block = powder_compact
[]
[powder_max_strain_yy]
type = ElementExtremeValue
variable = strain_yy
block = powder_compact
[]
[powder_max_stress_yy]
type = ElementExtremeValue
variable = stress_yy
block = powder_compact
[]
[powder_max_vonmises_stress]
type = ElementExtremeValue
variable = vonmises_stress
block = powder_compact
[]
[powder_max_plastic_strain_yy]
type = ElementExtremeValue
variable = plastic_strain_yy
block = powder_compact
[]
[powder_temperature_pp]
type = AverageNodalVariableValue
variable = temperature
block = powder_compact
[]
[plunger_strain_yy]
type = ElementAverageValue
variable = strain_yy
block = 'upper_plunger lower_plunger'
[]
[plunger_stress_yy]
type = ElementAverageValue
variable = stress_yy
block = 'upper_plunger lower_plunger'
[]
[plunger_vonmises_stress]
type = ElementAverageValue
variable = vonmises_stress
block = 'upper_plunger lower_plunger'
[]
[]
[Outputs]
csv = true
exodus = true
perf_graph = true
# [ckpt]
# type =Checkpoint
# time_step_interval = 1
# num_files = 2
# []
[]
(test/tests/materials/graphite/mechanical/ad_graphite_thermal_expansion.i)
# This model includes only mechanics on a patch mesh and is intended
# verify the implementation of the piecewise function for the coefficient of
# thermal expansion from Cincotti et al (2007) AIChE Journal, Vol 53 No 3,
# page 710 Figure 7d. The piecewise function is given as a function within
# the input file (thexp_exact), and the value of that function is compared
# at each timestep against the computed strain_{xx} calculated value
#
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[mesh]
type = ExamplePatchMeshGenerator
dim = 2
[]
coord_type = RZ
[]
[AuxVariables]
[temperature]
initial_condition = 300.0
[]
[]
[AuxKernels]
[temp_aux]
type = FunctionAux
variable = temperature
function = temp_ramp
[]
[]
[Functions]
[temp_ramp]
type = ParsedFunction
expression = '300 + 400.0/60.*t' #stand-in for a 400C/min heating rate
[]
[thermal_expansion]
type = ParsedFunction
symbol_names = T
symbol_values = temperature
expression = '(1.996e-6*log(T*4.799e-2)-4.041e-6)'
[]
[thexp_exact]
type = ParsedFunction
symbol_names = 'ref_temperature thermal_expansion temperature'
symbol_values = '300.0 thermal_expansion temperature'
expression = 'thermal_expansion * (temperature - ref_temperature)'
[]
[]
[Physics]
[SolidMechanics]
[QuasiStatic]
[all]
strain = SMALL
decomposition_method = EigenSolution
add_variables = true
use_automatic_differentiation = true
eigenstrain_names = graphite_thermal_expansion
generate_output = 'strain_xx'
[]
[]
[]
[]
[BCs]
[left]
type = ADDirichletBC
variable = disp_x
value = 0
boundary = left
[]
[bottom]
type = ADDirichletBC
variable = disp_y
value = 0
boundary = bottom
[]
[]
[Materials]
[elasticity_tensor]
type = ADComputeIsotropicElasticityTensor
youngs_modulus = 1
poissons_ratio = 0.3
[]
[graphite_thermal_expansion]
type = ADGraphiteThermalExpansionEigenstrain
eigenstrain_name = graphite_thermal_expansion
stress_free_temperature = 300
temperature = temperature
[]
[stress]
type = ADComputeLinearElasticStress
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
compute_scaling_once = true
petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap'
petsc_options_value = 'asm ilu 1'
nl_max_its = 20
l_max_its = 50
dt = 25
dtmin = 1.0e-8
end_time = 300
[]
[Postprocessors]
[temperature]
type = ElementAverageValue
variable = temperature
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
[]
[thexp_exact]
type = FunctionValuePostprocessor
function = thexp_exact
[]
[thexp_diff]
type = DifferencePostprocessor
value1 = strain_xx
value2 = thexp_exact
outputs = none
[]
[thexp_max_diff]
type = TimeExtremeValue
postprocessor = thexp_diff
value_type = abs_max
[]
[]
[Outputs]
csv = true
perf_graph = true
file_base = graphite_thermal_expansion_out
[]
(examples/sps/multiapp/diemodel/electrothermomechs/simple_geometries/graphite_electrothermalmech_singleblock.i)
## Units in the input file: m-Pa-s-K
#initial_porosity = 0.38
initial_temperature = 300 #roughly 600C where the pyrometer kicks in
[GlobalParams]
displacements = 'disp_x disp_y'
volumetric_locking_correction = false
order = SECOND
[]
[Mesh]
[top_square]
type = GeneratedMeshGenerator
dim = 2
nx = 20
ny = 20
xmax = 0.0005
ymin = 0.0005
ymax = 0.001
boundary_name_prefix = upper_plunger
elem_type = QUAD8
[]
[top_square_block]
type = SubdomainIDGenerator
input = top_square
subdomain_id = 1
[]
[two_blocks]
type = MeshCollectionGenerator
inputs = 'top_square_block'
[]
[block_rename]
type = RenameBlockGenerator
input = two_blocks
old_block = '1'
new_block = 'upper_plunger'
[]
[inner_centerpoint]
type = BoundingBoxNodeSetGenerator
input = block_rename
bottom_left = '0.0 0.00074 0'
top_right = '0.25e-5 0.000755 0'
new_boundary = 'centerpoint'
[]
[outer_centerpoint]
type = BoundingBoxNodeSetGenerator
input = inner_centerpoint
bottom_left = '0.00049 0.00074 0'
top_right = '5.25e-5 0.000755 0'
new_boundary = 'outer_centerpoint'
[]
coord_type = RZ
patch_update_strategy = iteration
patch_size = 10
[]
[Problem]
type = ReferenceResidualProblem
reference_vector = 'ref'
extra_tag_vectors = 'ref'
[]
[Variables]
[temperature]
initial_condition = ${initial_temperature}
[]
[electric_potential]
[]
[]
[AuxVariables]
[vonmises_stress]
family = MONOMIAL
order = FIRST
[]
# [thermal_conductivity_aeh]
# initial_condition = 0.4
# []
[specific_heat_capacity_va]
initial_condition = 842.2 # at 1500K #568.73 at 1000K #447.281 # at 293K
[]
[density_va]
initial_condition = 3106.2 ##5010.0*(1-${initial_porosity}) #in kg/m^3
[]
[heat_transfer_radiation]
[]
[thermal_conductivity]
family = MONOMIAL
order = FIRST
block = 'upper_plunger'
[]
[sigma_aeh]
initial_condition = 2.0e-10 #in units eV/((nV)^2-s-nm)
[]
# [electrical_conductivity]
# family = MONOMIAL
# order = FIRST
# block = 'upper_plunger'
# []
[microapp_potential]
#converted to microapp electronVolts units
[]
[E_x]
order = FIRST
family = MONOMIAL
[]
[E_y]
order = FIRST
family = MONOMIAL
[]
[graphite_current_density]
order = FIRST
family = MONOMIAL
[]
[current_density_J]
family = NEDELEC_ONE
order = FIRST
[]
[]
[Physics]
[SolidMechanics]
[QuasiStatic]
[graphite]
strain = FINITE
add_variables = true
use_automatic_differentiation = true
generate_output = 'strain_xx strain_xy strain_yy strain_zz stress_xx stress_xy stress_yy stress_zz'
extra_vector_tags = 'ref'
# eigenstrain_names = 'graphite_thermal_expansion'
block = 'upper_plunger'
[]
[]
[]
[]
[Kernels]
[HeatDiff_graphite]
type = ADHeatConduction
variable = temperature
thermal_conductivity = thermal_conductivity #use parsed material property, hope it works
extra_vector_tags = 'ref'
block = 'upper_plunger'
[]
[HeatTdot_graphite]
type = ADHeatConductionTimeDerivative
variable = temperature
specific_heat = heat_capacity #use parsed material property
density_name = graphite_density
extra_vector_tags = 'ref'
block = 'upper_plunger'
[]
[JouleHeating_graphite]
type = ADJouleHeatingSource
variable = temperature
elec = electric_potential
electrical_conductivity = electrical_conductivity
# use_displaced_mesh = true
extra_vector_tags = 'ref'
block = 'upper_plunger'
[]
[electric_graphite]
type = ADMatDiffusion
variable = electric_potential
diffusivity = electrical_conductivity
extra_vector_tags = 'ref'
block = 'upper_plunger'
[]
[]
[AuxKernels]
[vonmises_stress]
type = ADRankTwoScalarAux
variable = vonmises_stress
rank_two_tensor = stress
scalar_type = VonMisesStress
[]
[heat_transfer_radiation]
type = ParsedAux
variable = heat_transfer_radiation
boundary = 'upper_plunger_right'
coupled_variables = 'temperature'
constant_names = 'boltzmann epsilon temperature_farfield'
constant_expressions = '5.67e-8 0.85 300.0' #roughly room temperature, which is probably too cold
expression = '-boltzmann*epsilon*(temperature^4-temperature_farfield^4)'
[]
[electrical_conductivity]
type = ADMaterialRealAux
variable = electrical_conductivity
property = electrical_conductivity
block = 'upper_plunger'
[]
[thermal_conductivity]
type = ADMaterialRealAux
variable = thermal_conductivity
property = thermal_conductivity
block = 'upper_plunger'
[]
[E_x]
type = VariableGradientComponent
variable = E_x
gradient_variable = electric_potential
component = x
block = 'upper_plunger'
[]
[E_y]
type = VariableGradientComponent
variable = E_y
gradient_variable = electric_potential
component = y
block = 'upper_plunger'
[]
[graphite_current_density]
type = ParsedAux
variable = graphite_current_density
coupled_variables = 'electrical_conductivity E_y'
expression = '-1.0*electrical_conductivity*E_y'
block = 'upper_plunger'
[]
[current_density_J]
type = ADCurrentDensity
variable = current_density_J
potential = electric_potential
[]
[]
[BCs]
[centerline_no_dispx]
type = ADDirichletBC
preset = true
variable = disp_x
value = 0
boundary = 'upper_plunger_left'
[]
[fixed_in_y]
type = ADDirichletBC
preset = true
variable = disp_y
value = 0
boundary = 'centerpoint'
[]
[bottom_pressure_ydirection]
type = ADPressure
variable = disp_y
boundary = 'upper_plunger_bottom'
function = 20.7e6 #'if(t<1, 20.7e6*t, 20.7e6)'
[]
[top_pressure_ydirection]
type = ADPressure
variable = disp_y
boundary = 'upper_plunger_top'
function = 20.7e6 #'if(t<1, 20.7e6*t, 20.7e6)'
[]
# [top_temperature]
# type = ADFunctionDirichletBC
# variable = temperature
# function = 'if(t<1100, ${initial_temperature} + 50.0/60.0*t, 1800)' #stand-in for a 50C/min heating rate)'
# boundary = 'top_upper_plunger'
# []
[external_surface]
type = CoupledVarNeumannBC
boundary = 'upper_plunger_right'
variable = temperature
v = heat_transfer_radiation
[]
[electric_top]
type = ADFunctionDirichletBC
variable = electric_potential
boundary = 'upper_plunger_top'
function = 'if(t<20.0, 4.0e-3*t, 0.08)'
# value = 1.5 #2.5 #'3.5' ## V #+1.0e-6*t'
[]
[electric_bottom]
type = ADDirichletBC
variable = electric_potential
boundary = 'upper_plunger_bottom'
value = 0.0 #1.0 #0.0
[]
[]
[Functions]
[graphite_thermal_conductivity_fcn]
type = ParsedFunction
expression = 29.0 #'-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659'
[]
# [yttria_thermal_conductivity_fcn] #from the multiapp
# type = ParsedFunction
# expression = '3214.46/(t-147.73)'
# []
# [harmonic_mean_thermal_conductivity]
# type = ParsedFunction
# expression = '2*(-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)*(3214.46/(t-147.73))/((-4.418e-12*t^4+2.904e-8*t^3-4.688e-5*t^2-0.0316*t+119.659)+(3214.46/(t-147.73)))'
# # symbol_names = 'k_graphite k_yttria'
# # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
# []
[graphite_electrical_conductivity_fcn]
type = ParsedFunction
expression = 1.1e5 #'1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5)'
[]
# [yttria_electrical_conductivity_fcn]
# type = ParsedFunction
# # symbol_names = porosity
# # symbol_values = initial_porosity
# expression = '(1-0.62)*2.0025e4*exp(-1.61/8.617343e-5/t)'
# []
# [harmonic_mean_electrical_conductivity]
# type = ParsedFunction
# expression = '2*(1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))*((1)*2.0025e4*exp(-1.61/8.617343e-5/t))/((1.0/(-2.705e-15*t^3+1.263e-11*t^2-1.836e-8*t+1.813e-5))+((1)*2.0025e4*exp(-1.61/8.617343e-5/t)))'
# # symbol_names = 'k_graphite k_yttria'
# # symbol_values = 'graphite_thermal_conductivity_fcn yttria_thermal_conductivity_fcn'
# []
[]
[Materials]
[graphite_elasticity_tensor]
type = ADComputeIsotropicElasticityTensor
youngs_modulus = 1.08e10 #in Pa, 10.8 GPa, Cincotti
poissons_ratio = 0.33
block = 'upper_plunger'
[]
[graphite_stress]
type = ADComputeFiniteStrainElasticStress
block = 'upper_plunger'
[]
# [graphite_thermal_expansion]
# type = ADGraphiteThermalExpansionEigenstrain
# # type = ADComputeThermalExpansionEigenstrain
# # thermal_expansion_coeff = 3.5e-6
# eigenstrain_name = graphite_thermal_expansion
# stress_free_temperature = 300
# temperature = temperature
# block = 'upper_plunger'
# []
[graphite_density]
type = ADGenericConstantMaterial
prop_names = 'graphite_density'
prop_values = 1.750e3 #in kg/m^3 from Cincotti et al 2007, Table 2, doi:10.1002/aic
block = 'upper_plunger'
[]
[graphite_thermal]
type = ADGraphiteThermal
temperature = temperature
block = 'upper_plunger'
[]
[graphite_electrical_conductivity]
type = ADParsedMaterial
property_name = electrical_conductivity
coupled_variables = 'temperature'
expression = '1.0/(-2.705e-15*temperature^3+1.263e-11*temperature^2-1.836e-8*temperature+1.813e-5)'
output_properties = electrical_conductivity
outputs = 'exodus'
block = 'upper_plunger'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
line_search = 'none'
compute_scaling_once = false
# petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap'
# petsc_options_value = 'asm ilu 1'
# petsc_options_iname = '-ksp_max_it -ksp_gmres_restart -pc_type -snes_max_funcs -sub_pc_factor_levels'
# petsc_options_value = ' 100 100 asm 1000000 1'
# petsc_options_iname = '-pc_type -snes_linesearch_type -pc_factor_shift_type -pc_factor_shift_amount'
# petsc_options_value = 'lu basic NONZERO 1e-15'
# petsc_options_iname = '-pc_type -ksp_grmres_restart -sub_ksp_type -sub_pc_type -pc_asm_overlap'
# petsc_options_value = ' asm 101 preonly ilu 2'
# petsc_options = '-snes_converged_reason -ksp_converged_reason -options_left -ksp_monitor_singular_value'
#bison options
# petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu superlu_dist'
nl_rel_tol = 2e-5 #was 1e-10, for temperature only
nl_abs_tol = 2e-12 #was 1e-10, before that 1e-12
nl_max_its = 20
l_max_its = 50
dtmin = 1.0e-3
end_time = 1 #900 #15 minutes, rule of thumb from Dennis is 10 minutes
[Quadrature]
order = FIFTH
side_order = SEVENTH
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = 0.05
optimal_iterations = 8
iteration_window = 2
[]
[]
[Postprocessors]
[stress_xx]
type = ElementAverageValue
variable = stress_xx
[]
[strain_xx]
type = ElementAverageValue
variable = strain_xx
[]
[stress_yy]
type = ElementAverageValue
variable = stress_yy
[]
[strain_yy]
type = ElementAverageValue
variable = strain_yy
[]
[vonmises_stress]
type = ElementAverageValue
variable = vonmises_stress
[]
[temperature_pp]
type = AverageNodalVariableValue
variable = temperature
[]
[graphite_thermal_conductivity]
type = ElementAverageValue
variable = thermal_conductivity
block = 'upper_plunger'
[]
[graphite_sigma]
type = ElementAverageValue
variable = electrical_conductivity
block = 'upper_plunger'
[]
# [yttria_thermal_conductivity]
# type = ElementAverageValue
# variable = yttria_thermal_conductivity
# block = powder_compact
# []
# [temperature_powder_midpoint]
# type = PointValue
# variable = temperature
# point = '0.00025 0.00025 0'
# []
# [temperature_graphite_midpoint]
# type = PointValue
# variable = temperature
# point = '0.00025 0.00075 0'
# []
# [electric_potential_powder_midpoint]
# type = PointValue
# variable = electric_potential
# point = '0.00025 0.00025 0'
# []
# [electric_potential_graphite_midpoint]
# type = PointValue
# variable = electric_potential
# point = '0.00025 0.00075 0'
# []
[]
[Outputs]
csv = true
exodus = true
perf_graph = true
[]