- propName of the material property to provide the multiplier
C++ Type:MaterialPropertyName
Controllable:No
Description:Name of the material property to provide the multiplier
- variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Controllable:No
Description:The name of the variable that this residual object operates on
MaterialSource
Source term defined by the material property
Overview
This kernel implements the following weak form.
where is the test function, and the source defined by a material property.
Example Input File Syntax
[energy_balance_2]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
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
- coefficient1Coefficient to be multiplied to the source
Default:1
C++ Type:double
Controllable:No
Description:Coefficient to be multiplied to the source
- displacementsThe displacements
C++ Type:std::vector<VariableName>
Controllable:No
Description:The displacements
- 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
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.
Optional Parameters
- absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
- extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
- extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
- matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
- vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill
Tagging 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.
- diag_save_inThe name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Controllable:No
Description:The name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- 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
- save_inThe name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Controllable:No
Description:The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- 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
Input Files
- (examples/LiB/CC_discharging.i)
- (examples/SSB/CC_charging.i)
- (examples/SSB_3D/charging.i)
- (examples/PCM-GF/discharging.i)
- (examples/SSB_3D/CC_charging.i)
- (test/tests/chemical-electrical-thermal/thermal_effects.i)
- (test/tests/chemical-electrical-thermal-mechanical/pressure.i)
- (examples/arc-welding/welding.i)
- (examples/PCM-GF/equilibriate.i)
- (examples/SSB/charging.i)
- (examples/SSB/CV_charging.i)
- (examples/SSB_3D/CC_discharging.i)
- (examples/creep-diffusion/bulk/creep-diffusion.i)
- (examples/cable/cable.i)
- (examples/SSB_3D/CV_charging.i)
- (examples/LiB/CC_charging.i)
- (examples/LiB/CV_charging.i)
- (test/tests/electrical-thermal/jh.i)
- (examples/resin-printing/print.i)
- (examples/PCM-GF/charging_idling.i)
- (examples/SSB/CC_discharging.i)
(test/tests/electrical-thermal/jh.i)
n = 32
[GlobalParams]
energy_densities = 'E H'
[]
[Mesh]
[battery]
type = GeneratedMeshGenerator
dim = 1
xmin = 0
xmax = 1
nx = ${n}
[]
[]
[Variables]
[Phi]
[]
[T]
[]
[]
[Functions]
[Phi]
type = ParsedFunction
expression = '3.565*exp(x)-4.71*x^2+7.1*x-3.4'
[]
[sigma]
type = ParsedFunction
expression = '1+x'
[]
[q]
type = ParsedFunction
expression = '2.32+exp(x)*(-7.13-3.565*x)+18.84*x'
[]
[kappa]
type = ParsedFunction
expression = '-2.6*x+3'
[]
[]
[Kernels]
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge]
type = BodyForce
variable = Phi
function = q
[]
[energy_balance_1]
type = RankOneDivergence
variable = T
vector = h
[]
[energy_balance_2]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[BCs]
[fix]
type = FunctionDirichletBC
variable = Phi
boundary = 'left right'
function = Phi
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = right
value = -1
boundary_material = qconv
[]
[]
[Materials]
[electric_constants]
type = ADGenericFunctionMaterial
prop_names = 'sigma'
prop_values = 'sigma'
[]
[charge_trasport]
type = BulkChargeTransport
electrical_energy_density = E
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
[thermal_constants]
type = ADGenericFunctionMaterial
prop_names = 'kappa'
prop_values = 'kappa'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = H
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
[qconv]
type = ADParsedMaterial
property_name = qconv
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '1.35 300'
boundary = right
[]
[]
[Postprocessors]
[error]
type = ElementL2Error
variable = Phi
function = Phi
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
num_steps = 1
[]
[Outputs]
exodus = true
[]
(examples/LiB/CC_discharging.i)
I = 3e-3 #mA
width = 0.03 #mm
in = '${fparse I/width}'
t0 = '${fparse 1e-2/in}'
dt = '${fparse t0/100}'
sigma_a = 1e0 #mS/mm
sigma_e = 1e-1 #mS/mm
sigma_c = 1e-2 #mS/mm
l0 = 0
l1 = 0.04
l2 = 0.07
l3 = 0.12
cmin = 1e-4 #mmol/mm^3
cmax = 1e-3 #mmol/mm^3
D_a = 1e-3 #mm^2/s
D_e = 1e-4 #mm^2/s
D_c = 5e-5 #mm^2/s
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_c = 1e5
E_e = 1e4
E_a = 2e5
nu_c = 0.3
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 60
beta = 1e-3
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
T_penalty = 2e-1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q zeta m'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
[]
[Problem]
restart_file_base = 'CV_charging_I_${I}_cp/LATEST'
[]
[Mesh]
[battery]
type = GeneratedMeshGenerator
dim = 2
xmin = ${l0}
xmax = ${l3}
ymin = 0
ymax = ${width}
nx = 60
ny = 15
[]
[anode]
type = SubdomainBoundingBoxGenerator
input = battery
block_id = 1
block_name = anode
bottom_left = '${l0} 0 0'
top_right = '${l1} ${width} 0'
[]
[elyte]
type = SubdomainBoundingBoxGenerator
input = anode
block_id = 2
block_name = elyte
bottom_left = '${l1} 0 0'
top_right = '${l2} ${width} 0'
[]
[cathode]
type = SubdomainBoundingBoxGenerator
input = elyte
block_id = 3
block_name = cathode
bottom_left = '${l2} 0 0'
top_right = '${l3} ${width} 0'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = cathode
add_interface_on_two_sides = true
split_interface = true
[]
[]
[Variables]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[T]
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[Phi0]
[]
[]
[ICs]
[c_ref_min]
type = ConstantIC
variable = c_ref
value = ${cmin}
block = 'anode'
[]
[c_ref_mid]
type = ConstantIC
variable = c_ref
value = '${fparse (cmax+cmin)/2}'
block = 'elyte'
[]
[c_ref_max]
type = ConstantIC
variable = c_ref
value = ${cmax}
block = 'cathode'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'elyte_anode cathode_elyte'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'anode_elyte elyte_cathode'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'anode_elyte elyte_cathode elyte_anode cathode_elyte'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[left]
type = FunctionNeumannBC
variable = Phi
boundary = left
function = in
[]
[right]
type = DirichletBC
variable = Phi
boundary = right
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[open]
type = OpenBC
variable = c
flux = jm
boundary = 'left right'
[]
[]
[Constraints]
[y]
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'anode ${sigma_a} elyte ${sigma_e} cathode ${sigma_c}'
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'anode ${D_a} elyte ${D_e} cathode ${D_c}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T
reference_concentration = c_ref
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
[]
# Migration
[migration]
type = Migration
electrochemical_energy_density = m
electric_potential = Phi
chemical_potential = mu
electric_conductivity = sigma
faraday_constant = ${F}
[]
[migration_flux]
type = MassFlux
mass_flux = jm
energy_densities = 'm'
chemical_potential = mu
[]
# Redox
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)'
args = c
block = 'anode'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)'
args = c
block = 'cathode'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'anode_elyte'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'elyte_anode'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cathode_elyte'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'elyte_cathode'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
# Mechanical
[stiffness_c]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_c*nu_c/(1+nu_c)/(1-2*nu_c)} ${fparse E_c/2/(1+nu_c)}'
block = cathode
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = elyte
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = anode
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
displacements = 'disp_x disp_y'
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_r - V_l'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[I]
type = ADSideIntegralMaterialProperty
property = i
component = 0
boundary = right
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = '-dt*I'
pp_names = 'dt I'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_a_max]
type = NodalExtremeValue
variable = c
value_type = max
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_c_min]
type = NodalExtremeValue
variable = c
value_type = min
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_a_min]
type = NodalExtremeValue
variable = c
value_type = min
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_c_max]
type = NodalExtremeValue
variable = c
value_type = max
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_a]
type = ElementIntegralVariablePostprocessor
variable = c
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_e]
type = ElementIntegralVariablePostprocessor
variable = c
block = elyte
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_c]
type = ElementIntegralVariablePostprocessor
variable = c
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill_a]
type = Terminator
expression = 'c_a_min <= ${cmin}'
message = 'Concentration in anode is below the minimum allowable value.'
[]
[kill_c]
type = Terminator
expression = 'c_c_max >= ${cmax}'
message = 'Concentration in cathode exceeds the maximum allowable value.'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
nl_rel_tol = 1e-6
nl_abs_tol = 1e-10
nl_max_its = 20
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
optimal_iterations = 7
iteration_window = 2
growth_factor = 1.2
cutback_factor = 0.5
cutback_factor_at_failure = 0.2
linear_iteration_ratio = 1000000
[]
start_time = 0
end_time = 100000
[]
[Outputs]
file_base = 'CC_discharging_I_${I}'
csv = true
exodus = true
print_linear_residuals = false
checkpoint = true
[]
(examples/SSB/CC_charging.i)
I = 0.005 #mA
width = 0.05 #mm
in = '${fparse -I/width}'
t0 = '${fparse -1e-2/in}'
sigma_a = 0.2 #mS/mm
sigma_e = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.2 #mS/mm
sigma_cm = 0.05 #mS/mm
Phi_penalty = 10
cmin_a = 1e-4 #mmol/mm^3
cmax_a = 1e-3 #mmol/mm^3
c_e = 5e-4 #mmol/mm^3
cmax_c = 1e-3 #mmol/mm^3
c_ref_entropy = 5e-5
D_cp = 5e-5 #mm^2/s
D_cm = 1e-4 #mm^2/s
D_a = 5e-4 #mm^2/s
D_e = 1e-4 #mm^2/s
c_penalty = 1
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_cp = 6e4
E_cm = 5e4
E_e = 5e4
E_a = 1e5
nu_cp = 0.3
nu_cm = 0.25
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 140
beta = 1e-4
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
htc = 9.5e-3
T_penalty = 1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
displacements = 'disp_x disp_y'
[]
[Mesh]
[battery]
type = FileMeshGenerator
file = 'gold/ssb.msh'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = battery
add_interface_on_two_sides = true
split_interface = true
[]
use_displaced_mesh = false
[]
[Variables]
[Phi_ca]
block = cm
[]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[T]
initial_condition = ${T0}
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[stress]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = pk1
scalar_type = Hydrostatic
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Js]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Fs
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Jt]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Ft
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Phi0]
[]
[]
[ICs]
[c_a]
type = ConstantIC
variable = c
value = ${cmin_a}
block = 'a'
[]
[c_e]
type = ConstantIC
variable = c
value = ${c_e}
block = 'cm e'
[]
[c_c]
type = ConstantIC
variable = c
value = ${cmax_c}
block = 'cp'
[]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${cmin_a}
block = 'a'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c_e}
block = 'cm e'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${cmax_c}
block = 'cp'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge_balance_ca]
type = RankOneDivergence
variable = Phi_ca
vector = i_ca
block = cm
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'a_e cm_cp'
[]
[negative_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = 1
boundary = 'a_e cm_cp'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'a_e cm_cp e_a cp_cm'
[]
[continuity_c]
type = InterfaceContinuity
variable = c
neighbor_var = c
penalty = ${c_penalty}
boundary = 'cm_e'
[]
[continuity_Phi_ca]
type = InterfaceContinuity
variable = Phi_ca
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_cp'
[]
[continuity_Phi]
type = InterfaceContinuity
variable = Phi
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_e'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[current]
type = FunctionNeumannBC
variable = Phi
boundary = right
function = in
[]
[potential]
type = DirichletBC
variable = Phi_ca
boundary = left
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'left right'
value = -1
boundary_material = qconv
[]
[]
[Constraints]
[ev_y]
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'a ${sigma_a} e ${sigma_e} cm ${sigma_cm} cp ${sigma_cp}'
[]
[conductivity_ca]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma_ca'
subdomain_to_prop_value = 'cm ${sigma_ca}'
block = cm
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[charge_transport_ca]
type = BulkChargeTransport
electrical_energy_density = q_ca
electric_potential = Phi_ca
electric_conductivity = sigma_ca
temperature = T
block = cm
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
output_properties = i
outputs = exodus
[]
[current_density_ca]
type = CurrentDensity
current_density = i_ca
electric_potential = Phi_ca
output_properties = i_ca
outputs = exodus
block = cm
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'a ${D_a} e ${D_e} cm ${D_cm} cp ${D_cp}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T_ref
reference_concentration = ${c_ref_entropy}
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta m'
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
outputs = exodus
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_a}; 2.77e-4*x^2-0.0069*x+0.0785'
# function = 'x:=c/${cmax_a}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
args = c
material_property_names = 'ramp'
block = 'a'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_c}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
args = c
material_property_names = 'ramp'
block = 'cp'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'a_e'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'e_a'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cp_cm'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cm_cp'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
output_properties = h
outputs = exodus
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
output_properties = r
outputs = exodus
[]
[conv]
type = ADParsedMaterial
f_name = qconv
function = '${htc}*(T-T_ref)'
args = 'T T_ref'
boundary = 'left right'
[]
# Mechanical
[stiffness_cp]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cp*nu_cp/(1+nu_cp)/(1-2*nu_cp)} ${fparse E_cp/2/(1+nu_cp)}'
block = cp
[]
[stiffness_cm]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cm*nu_cm/(1+nu_cm)/(1-2*nu_cm)} ${fparse E_cm/2/(1+nu_cm)}'
block = cm
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = e
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = a
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
non_swelling_pressure = p
output_properties = 'p'
outputs = exodus
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi_ca
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_l - V_r'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[I]
type = ADSideIntegralMaterialProperty
property = i
component = 0
boundary = right
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = 'dt*I'
pp_names = 'dt I'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill_V]
type = Terminator
expression = 'V >= 4.6'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-ksp_converged_reason'
# petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_truncfactor'
# petsc_options_value = 'hypre boomeramg 301 0.25 ext+i PMIS 4 2 0.4'
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
ignore_variables_for_autoscaling = 'T'
verbose = true
line_search = none
l_max_its = 300
l_tol = 1e-6
nl_rel_tol = 1e-6
nl_abs_tol = 1e-9
nl_max_its = 12
[Predictor]
type = SimplePredictor
scale = 1
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse t0/50}'
optimal_iterations = 6
iteration_window = 1
growth_factor = 1.2
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
linear_iteration_ratio = 100
[]
end_time = 10000
[]
[Outputs]
exodus = true
csv = true
print_linear_residuals = false
checkpoint = true
[]
(examples/SSB_3D/charging.i)
I = 2.5e-4 #mA
width = 0.05 #mm
in = '${fparse -I/width/width}'
t0 = '${fparse -1e-2/in}'
dt = '${fparse t0/100}'
sigma_a = 0.2 #mS/mm
sigma_e = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.5 #mS/mm
sigma_cm = 0.05 #mS/mm
Phi_penalty = 100
cmin_a = 1e-4 #mmol/mm^3
cmax_a = 1e-3 #mmol/mm^3
c_e = 5e-4 #mmol/mm^3
cmin_c = 1e-4 #mmol/mm^3
cmax_c = 1e-3 #mmol/mm^3
c_ref_entropy = 5e-5
D_cp = 5e-5 #mm^2/s
D_cm = 1e-4 #mm^2/s
D_a = 5e-4 #mm^2/s
D_e = 1e-4 #mm^2/s
c_penalty = 5e-1
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-4 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_cp = 6e4
E_cm = 5e4
E_e = 5e4
E_a = 1e5
nu_cp = 0.3
nu_cm = 0.25
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 60
beta = 1e-4
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
T_penalty = 2
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
[battery]
type = FileMeshGenerator
file = 'coarse.e'
[]
[scale]
type = TransformGenerator
input = battery
transform = SCALE
vector_value = '1e-3 1e-3 1e-3' #um to mm
[]
[cathode_particle]
type = RenameBlockGenerator
input = scale
old_block = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44'
new_block = 'cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp'
[]
[cathode_matrix]
type = RenameBlockGenerator
input = cathode_particle
old_block = '48'
new_block = 'cm'
[]
[elyte]
type = RenameBlockGenerator
input = cathode_matrix
old_block = 49
new_block = 'e'
[]
[anode]
type = RenameBlockGenerator
input = elyte
old_block = 50
new_block = 'a'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = anode
add_interface_on_two_sides = true
split_interface = true
[]
[sidesets]
type = SideSetsFromNormalsGenerator
input = interfaces
normals = '-1 0 0 1 0 0'
new_boundary = 'left right'
[]
use_displaced_mesh = false
[]
[Variables]
[Phi_ca]
block = cm
[]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[T]
initial_condition = ${T0}
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[stress]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = pk1
scalar_type = VonMisesStress
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[]
[ICs]
[c_a]
type = ConstantIC
variable = c
value = ${cmin_a}
block = 'a'
[]
[c_e]
type = ConstantIC
variable = c
value = ${c_e}
block = 'cm e'
[]
[c_c]
type = ConstantIC
variable = c
value = ${cmax_c}
block = 'cp'
[]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${cmin_a}
block = 'a'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c_e}
block = 'cm e'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${cmax_c}
block = 'cp'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge_balance_ca]
type = RankOneDivergence
variable = Phi_ca
vector = i_ca
block = cm
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
[momentum_balance_z]
type = RankTwoDivergence
variable = disp_z
component = 2
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'a_e cm_cp'
[]
[negative_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = 1
boundary = 'a_e cm_cp'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'a_e cm_cp e_a cp_cm'
[]
[continuity_c]
type = InterfaceContinuity
variable = c
neighbor_var = c
penalty = ${c_penalty}
boundary = 'cm_e'
[]
[continuity_Phi_ca]
type = InterfaceContinuity
variable = Phi_ca
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_cp'
[]
[continuity_Phi]
type = InterfaceContinuity
variable = Phi
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_e'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_z]
type = InterfaceContinuity
variable = disp_z
neighbor_var = disp_z
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[current]
type = FunctionNeumannBC
variable = Phi
boundary = right
function = in
[]
[potential]
type = DirichletBC
variable = Phi_ca
boundary = left
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'left right'
[]
[fix_z]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'left right'
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'a ${sigma_a} e ${sigma_e} cm ${sigma_cm} cp ${sigma_cp}'
[]
[conductivity_ca]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma_ca'
subdomain_to_prop_value = 'cm ${sigma_ca}'
block = cm
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[charge_transport_ca]
type = BulkChargeTransport
electrical_energy_density = q_ca
electric_potential = Phi_ca
electric_conductivity = sigma_ca
temperature = T
block = cm
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
[current_density_ca]
type = CurrentDensity
current_density = i_ca
electric_potential = Phi_ca
block = cm
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'a ${D_a} e ${D_e} cm ${D_cm} cp ${D_cp}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T_ref
reference_concentration = ${c_ref_entropy}
[]
[diffusion]
type = CondensedMassDiffusion
mass_flux = j
mobility = M
concentration = c
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_a}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
args = c
material_property_names = 'ramp'
block = 'a'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_c}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
args = c
material_property_names = 'ramp'
block = 'cp'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'a_e'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'e_a'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cp_cm'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cm_cp'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
# Mechanical
[stiffness_cp]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cp*nu_cp/(1+nu_cp)/(1-2*nu_cp)} ${fparse E_cp/2/(1+nu_cp)}'
block = cp
[]
[stiffness_cm]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cm*nu_cm/(1+nu_cm)/(1-2*nu_cm)} ${fparse E_cm/2/(1+nu_cm)}'
block = cm
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = e
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = a
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi_ca
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_l - V_r'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[in]
type = FunctionValuePostprocessor
function = in
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = '-dt*in*${width}'
pp_names = 'dt in'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[cmin_c]
type = NodalExtremeValue
variable = c
value_type = min
block = 'cp'
[]
[cmax_a]
type = NodalExtremeValue
variable = c
value_type = max
block = 'a'
[]
[]
[UserObjects]
[kill_a]
type = Terminator
expression = 'cmax_a >= ${cmax_a}'
message = 'Concentration in anode exceeds the maximum allowable value.'
[]
[kill_cp]
type = Terminator
expression = 'cmin_c <= ${cmin_c}'
message = 'Concentration in cathode particle is below the minimum allowable value.'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-ksp_converged_reason'
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_truncfactor'
petsc_options_value = 'hypre boomeramg 601 0.6 ext+i PMIS 4 2 0.4'
automatic_scaling = true
ignore_variables_for_autoscaling = 'T'
verbose = true
l_max_its = 600
l_tol = 1e-6
nl_rel_tol = 1e-6
nl_abs_tol = 1e-9
nl_max_its = 12
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
optimal_iterations = 6
iteration_window = 1
growth_factor = 1.2
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
linear_iteration_ratio = 300
[]
end_time = 10000
[]
[Outputs]
[exo]
type = Exodus
interval = 5
file_base = '${outname}_I_${I}'
[]
[csv]
type = CSV
file_base = '${outname}_I_${I}'
[]
print_linear_residuals = false
[]
(examples/PCM-GF/discharging.i)
# units are in meter kelvin second (m,kg,s)
end_time = 14400 # 4 hrs
dtmax = 5
dt = 1
sigma_PCM = 5 # (from Wen's measurement of Gfoam+PCM in radial direction) (from Cfoam 70% dense foam = 28571.43) S/m (1/electrical resistivity (0.000035 ohm-m))
kappa_PCM = 10 #18.8 # (average of Kxy = 14 W/m-K, Kz = 23.6 W/mK at T=700C) #from Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
rho_PCM = 2050 # kg/m^3 #from Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
cp_PCM = 1074 # J/kg-K #from Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
sigma_pipe = 750750.75 # S/m (resistivity 1.332e-6 ohm-m at T = 700C) #Special metal data sheet
kappa_pipe = 23.9 # W/m-K (at 700C) #Special metal datasheet
rho_pipe = 8359.33 #kg/m^3
cp_pipe = 419 # J/kg-K
sigma_gas = 1e-12
kappa_gas = 0.03 #file:///C:/Users/barua/Downloads/PDS-FOAMGLAS%20ONE-US-en.pdf
rho_gas = 1.29 #file:///C:/Users/barua/Downloads/PDS-FOAMGLAS%20ONE-US-en.pdf
cp_gas = 1000 #file:///C:/Users/barua/Downloads/PDS-FOAMGLAS%20ONE-US-en.pdf
htc_gas = 0.1
T_inf_gas = 500
htc_insul = 5
T_inf_insul = 300
T_target = 800
[GlobalParams]
energy_densities = 'E H'
[]
[Mesh]
[fmg]
type = FileMeshGenerator
file = 'equilibriate_out.e'
use_for_exodus_restart = true
[]
coord_type = RZ
[]
[Variables]
[Phi]
[]
[T]
initial_from_file_var = T
[]
[v]
family = LAGRANGE_VEC
block = gas
[]
[p]
block = gas
[]
[]
[ICs]
[vel]
type = VectorConstantIC
variable = v
x_value = 1e-15
y_value = 1e-15
block = gas
[]
[]
[AuxVariables]
[T_old]
[AuxKernel]
type = ParsedAux
expression = 'T'
coupled_variables = 'T'
execute_on = 'INITIAL TIMESTEP_BEGIN'
[]
[]
[]
[Kernels]
[mass]
type = INSADMass
variable = p
block = gas
[]
[pspg]
type = INSADMassPSPG
variable = p
block = gas
[]
[momentum_convection]
type = INSADMomentumAdvection
variable = v
block = gas
[]
[momentum_viscous]
type = INSADMomentumViscous
variable = v
block = gas
[]
[momentum_pressure]
type = INSADMomentumPressure
variable = v
pressure = p
integrate_p_by_parts = true
block = gas
[]
[momentum_supg]
type = INSADMomentumSUPG
variable = v
velocity = v
block = gas
[]
[temperature_advection]
type = INSADEnergyAdvection
variable = T
block = gas
[]
[temperature_supg]
type = INSADEnergySUPG
variable = T
velocity = v
block = gas
[]
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cp
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[energy_balance_3]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[BCs]
[ground]
type = DirichletBC
variable = Phi
boundary = 'PCM_left'
value = 0
[]
[voltage]
type = TargetFeedbackDirichletBC
variable = Phi
monitor = T_outlet
target = ${T_target}
window = '${fparse 0.02*T_target}'
idle_value = 0
maintain_value = 12
compensate_value = 20
boundary = 'PCM_right'
[]
[T_inlet]
type = DirichletBC
variable = T
value = 300
boundary = 'inlet'
[]
[velocity_inlet]
type = VectorFunctionDirichletBC
variable = v
function_y = 0.1
boundary = 'inlet'
[]
[wall]
type = VectorFunctionDirichletBC
variable = v
boundary = 'wall'
[]
[hconv_outlet]
type = ADMatNeumannBC
variable = T
boundary = 'outlet'
value = -1
boundary_material = qconv_outlet
[]
[hconv_insul]
type = ADMatNeumannBC
variable = T
boundary = 'insul PCM_right'
value = -1
boundary_material = qconv_insul
[]
[]
[Materials]
[constant]
type = ADGenericConstantMaterial
prop_names = 'mu'
prop_values = '1.8e-5'
[]
[ins]
type = INSADStabilized3Eqn
pressure = p
velocity = v
temperature = T
k_name = kappa
block = gas
[]
[electrical_conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = sigma
subdomain_to_prop_value = 'PCM ${sigma_PCM} pipe ${sigma_pipe} gas ${sigma_gas}'
[]
[charge_trasport]
type = BulkChargeTransport
electrical_energy_density = E
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
[thermal_conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = kappa
subdomain_to_prop_value = 'PCM ${kappa_PCM} pipe ${kappa_pipe} gas ${kappa_gas}'
[]
[density]
type = ADPiecewiseConstantByBlockMaterial
prop_name = rho
subdomain_to_prop_value = 'PCM ${rho_PCM} pipe ${rho_pipe} gas ${rho_gas}'
[]
[specific_heat]
type = ADPiecewiseConstantByBlockMaterial
prop_name = cp
subdomain_to_prop_value = 'PCM ${cp_PCM} pipe ${cp_pipe} gas ${cp_gas}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = H
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
[qconv_outlet]
type = ADParsedMaterial
property_name = qconv_outlet
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '${htc_gas} ${T_inf_gas}'
boundary = 'outlet'
[]
[qconv_insul]
type = ADParsedMaterial
property_name = qconv_insul
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '${htc_insul} ${T_inf_insul}'
boundary = 'insul PCM_right'
[]
[delta_enthalpy]
type = ADParsedMaterial
property_name = delta_enthalpy
expression = 'rho*cp*(T-T_old)/2'
material_property_names = 'rho cp'
coupled_variables = 'T T_old'
[]
[]
[Postprocessors]
[T_outlet]
type = SideAverageValue
variable = T
boundary = 'outlet'
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_factor_shift_type'
petsc_options_value = 'lu NONZERO'
automatic_scaling = true
end_time = ${end_time}
dtmax = ${dtmax}
dtmin = 0.01
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
growth_factor = 1.2
optimal_iterations = 7
iteration_window = 2
linear_iteration_ratio = 100000
[]
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
nl_abs_tol = 1e-8
nl_rel_tol = 1e-6
nl_max_its = 12
[]
[Outputs]
file_base = 'T_target_${T_target}'
exodus = true
csv = true
[]
(examples/SSB_3D/CC_charging.i)
I = 0.005 #mA
width = 0.05 #mm
in = '${fparse -I/width/width}'
t0 = '${fparse -1e-2/in}'
sigma_a = 0.2 #mS/mm
sigma_e = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.2 #mS/mm
sigma_cm = 0.05 #mS/mm
Phi_penalty = 10
cmin_a = 1e-4 #mmol/mm^3
cmax_a = 1e-3 #mmol/mm^3
c_e = 5e-4 #mmol/mm^3
cmax_c = 1e-3 #mmol/mm^3
c_ref_entropy = 5e-5
D_cp = 5e-5 #mm^2/s
D_cm = 1e-4 #mm^2/s
D_a = 5e-4 #mm^2/s
D_e = 1e-4 #mm^2/s
c_penalty = 1
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_cp = 6e4
E_cm = 5e4
E_e = 5e4
E_a = 1e5
nu_cp = 0.3
nu_cm = 0.25
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 140
beta = 1e-4
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
htc = 9.5e-3
T_penalty = 1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
displacements = 'disp_x disp_y disp_z'
[]
[Mesh]
[battery]
type = FileMeshGenerator
file = '../coarse.e'
[]
[scale]
type = TransformGenerator
input = battery
transform = SCALE
vector_value = '1e-3 1e-3 1e-3' #um to mm
[]
[cathode_particle]
type = RenameBlockGenerator
input = scale
old_block = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44'
new_block = 'cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp'
[]
[cathode_matrix]
type = RenameBlockGenerator
input = cathode_particle
old_block = '48'
new_block = 'cm'
[]
[elyte]
type = RenameBlockGenerator
input = cathode_matrix
old_block = 49
new_block = 'e'
[]
[anode]
type = RenameBlockGenerator
input = elyte
old_block = 50
new_block = 'a'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = anode
add_interface_on_two_sides = true
split_interface = true
[]
[sidesets]
type = SideSetsFromNormalsGenerator
input = interfaces
normals = '-1 0 0 1 0 0 0 -1 0 0 1 0 0 0 -1 0 0 1'
new_boundary = 'left right bottom top back front'
[]
use_displaced_mesh = false
[]
[Variables]
[Phi_ca]
block = cm
[]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[T]
initial_condition = ${T0}
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[stress]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = pk1
scalar_type = Hydrostatic
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Js]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Fs
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Jt]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Ft
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Phi0]
[]
[]
[ICs]
[c_a]
type = ConstantIC
variable = c
value = ${cmin_a}
block = 'a'
[]
[c_e]
type = ConstantIC
variable = c
value = ${c_e}
block = 'cm e'
[]
[c_c]
type = ConstantIC
variable = c
value = ${cmax_c}
block = 'cp'
[]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${cmin_a}
block = 'a'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c_e}
block = 'cm e'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${cmax_c}
block = 'cp'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge_balance_ca]
type = RankOneDivergence
variable = Phi_ca
vector = i_ca
block = cm
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
[momentum_balance_z]
type = RankTwoDivergence
variable = disp_z
component = 2
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'a_e cm_cp'
[]
[negative_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = 1
boundary = 'a_e cm_cp'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'a_e cm_cp e_a cp_cm'
[]
[continuity_c]
type = InterfaceContinuity
variable = c
neighbor_var = c
penalty = ${c_penalty}
boundary = 'cm_e'
[]
[continuity_Phi_ca]
type = InterfaceContinuity
variable = Phi_ca
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_cp'
[]
[continuity_Phi]
type = InterfaceContinuity
variable = Phi
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_e'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_z]
type = InterfaceContinuity
variable = disp_z
neighbor_var = disp_z
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[current]
type = FunctionNeumannBC
variable = Phi
boundary = right
function = in
[]
[potential]
type = DirichletBC
variable = Phi_ca
boundary = left
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[fix_z]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back'
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'left right'
value = -1
boundary_material = qconv
[]
[]
[Constraints]
[ev_y]
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[ev_z]
type = EqualValueBoundaryConstraint
variable = disp_z
penalty = ${u_penalty}
secondary = front
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'a ${sigma_a} e ${sigma_e} cm ${sigma_cm} cp ${sigma_cp}'
[]
[conductivity_ca]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma_ca'
subdomain_to_prop_value = 'cm ${sigma_ca}'
block = cm
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[charge_transport_ca]
type = BulkChargeTransport
electrical_energy_density = q_ca
electric_potential = Phi_ca
electric_conductivity = sigma_ca
temperature = T
block = cm
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
output_properties = i
outputs = exodus
[]
[current_density_ca]
type = CurrentDensity
current_density = i_ca
electric_potential = Phi_ca
output_properties = i_ca
outputs = exodus
block = cm
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'a ${D_a} e ${D_e} cm ${D_cm} cp ${D_cp}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T_ref
reference_concentration = ${c_ref_entropy}
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta m'
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
outputs = exodus
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_a}; 2.77e-4*x^2-0.0069*x+0.0785'
# function = 'x:=c/${cmax_a}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
args = c
material_property_names = 'ramp'
block = 'a'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_c}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
args = c
material_property_names = 'ramp'
block = 'cp'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'a_e'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'e_a'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cp_cm'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cm_cp'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
output_properties = h
outputs = exodus
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
output_properties = r
outputs = exodus
[]
[conv]
type = ADParsedMaterial
f_name = qconv
function = '${htc}*(T-T_ref)'
args = 'T T_ref'
boundary = 'left right'
[]
# Mechanical
[stiffness_cp]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cp*nu_cp/(1+nu_cp)/(1-2*nu_cp)} ${fparse E_cp/2/(1+nu_cp)}'
block = cp
[]
[stiffness_cm]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cm*nu_cm/(1+nu_cm)/(1-2*nu_cm)} ${fparse E_cm/2/(1+nu_cm)}'
block = cm
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = e
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = a
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
non_swelling_pressure = p
output_properties = 'p'
outputs = exodus
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi_ca
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_l - V_r'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[I]
type = ADSideIntegralMaterialProperty
property = i
component = 0
boundary = right
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = 'dt*I'
pp_names = 'dt I'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill_V]
type = Terminator
expression = 'V >= 4.6'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-ksp_converged_reason'
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_truncfactor'
petsc_options_value = 'hypre boomeramg 301 0.7 ext+i PMIS 4 2 0.4'
automatic_scaling = true
ignore_variables_for_autoscaling = 'T'
verbose = true
line_search = none
l_max_its = 300
l_tol = 1e-6
nl_rel_tol = 1e-6
nl_abs_tol = 1e-9
nl_max_its = 12
[Predictor]
type = SimplePredictor
scale = 1
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse t0/50}'
optimal_iterations = 6
iteration_window = 1
growth_factor = 1.2
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
linear_iteration_ratio = 100
[]
end_time = 10000
[]
[Outputs]
file_base = '${outname}_CC_charging'
exodus = true
csv = true
print_linear_residuals = false
checkpoint = true
[]
(test/tests/chemical-electrical-thermal/thermal_effects.i)
R = 8.3145 #mJ/mmol/K
F = 96485 #mC/mmol
I = 3e-3 #mA
width = 0.03 #mm
in = '${fparse -I/width}'
t0 = '${fparse -1e-2/in}'
dt = '${fparse t0/100}'
Vmax = 4.3 #V
vf_se = 0.3
vf_cp = 0.5
vf_ca = 0.2
sigma_a = 0.2 #mS/mm
sigma_se = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.2 #mS/mm
sigma_e = ${sigma_se}
sigma_c = '${fparse vf_se*sigma_se + vf_cp*sigma_cp + vf_ca*sigma_ca}'
l0 = 0
l1 = 0.04
l2 = 0.07
l3 = 0.12
cmax = 1e-3 #mmol/mm^3
c0_a = 1e-4
c0_e = 5e-4
c0_c = 1e-3
M_a = 8e-11
M_se = 1e-11
M_cp = 4e-14
M_ca = 1e-13
M_e = ${M_se}
M_c = '${fparse vf_se*M_se + vf_cp*M_cp + vf_ca*M_ca}'
T0 = 300 #K
i0_a = 0.1 #mA/mm^2
i0_c = 0.1 #mA/mm^2
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
hconv = 9.5e-3 #mJ/mm^2/K/s
T_penalty = 2e-1
[GlobalParams]
energy_densities = 'dot(psi_c) q zeta chi m'
[]
[Mesh]
[battery]
type = GeneratedMeshGenerator
dim = 2
xmin = ${l0}
xmax = ${l3}
ymin = 0
ymax = ${width}
nx = 60
ny = 15
[]
[anode]
type = SubdomainBoundingBoxGenerator
input = battery
block_id = 1
block_name = anode
bottom_left = '${l0} 0 0'
top_right = '${l1} ${width} 0'
[]
[elyte]
type = SubdomainBoundingBoxGenerator
input = anode
block_id = 2
block_name = elyte
bottom_left = '${l1} 0 0'
top_right = '${l2} ${width} 0'
[]
[cathode]
type = SubdomainBoundingBoxGenerator
input = elyte
block_id = 3
block_name = cathode
bottom_left = '${l2} 0 0'
top_right = '${l3} ${width} 0'
[]
[anode_elyte]
type = BreakMeshByBlockGenerator
input = cathode
block_pairs = '1 2'
add_interface_on_two_sides = true
split_interface = true
[]
[cathode_elyte]
type = BreakMeshByBlockGenerator
input = anode_elyte
block_pairs = '2 3'
add_interface_on_two_sides = true
split_interface = true
[]
[]
[Variables]
[Phi]
[]
[c]
[]
[T]
initial_condition = ${T0}
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[]
[ICs]
[c_a]
type = ConstantIC
variable = c
value = ${c0_a}
block = 'anode'
[]
[c_e]
type = ConstantIC
variable = c
value = ${c0_e}
block = 'elyte'
[]
[c_c]
type = ConstantIC
variable = c
value = ${c0_c}
block = 'cathode'
[]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${c0_a}
block = 'anode'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c0_e}
block = 'elyte'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${c0_c}
block = 'cathode'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
factor_neighbor = 1
boundary = 'cathode_elyte'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'anode_elyte'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
factor_neighbor = 1
boundary = 'anode_elyte cathode_elyte'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'anode_elyte cathode_elyte'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[left]
type = FunctionNeumannBC
variable = Phi
boundary = left
function = in
[]
[right]
type = DirichletBC
variable = Phi
boundary = right
value = 0
[]
[open]
type = OpenBC
variable = c
flux = jm
boundary = 'left right'
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'left right'
value = -1
boundary_material = qconv
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'anode ${sigma_a} elyte ${sigma_e} cathode ${sigma_c}'
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
# Migration
[migration]
type = Migration
electrochemical_energy_density = m
electric_potential = Phi
chemical_potential = mu
electric_conductivity = sigma
faraday_constant = ${F}
[]
[migration_flux]
type = MassFlux
mass_flux = jm
energy_densities = 'm'
chemical_potential = mu
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'M'
subdomain_to_prop_value = 'anode ${M_a} elyte ${M_e} cathode ${M_c}'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T
reference_concentration = c_ref
reference_chemical_potential = 1e3
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
property_name = U
expression = 'x:=c/${cmax}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
coupled_variables = c
material_property_names = 'ramp'
boundary = 'anode_elyte'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
property_name = U
expression = 'x:=c/${cmax}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
coupled_variables = c
material_property_names = 'ramp'
boundary = 'cathode_elyte'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'anode_elyte'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cathode_elyte'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
[qconv]
type = ADParsedMaterial
property_name = qconv
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '${hconv} ${T0}'
boundary = 'left right'
[]
[enthalpy]
type = ADParsedMaterial
f_name = enthalpy
function = 'rho*cv*(T-T_ref)'
args = 'T T_ref'
material_property_names = 'rho cv'
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_r - V_l'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[in]
type = FunctionValuePostprocessor
function = in
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = '-dt*in*${width}'
pp_names = 'dt in'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_a_max]
type = NodalExtremeValue
variable = c
value_type = max
block = anode
[]
[c_c_min]
type = NodalExtremeValue
variable = c
value_type = min
block = cathode
[]
[mass_a]
type = ElementIntegralVariablePostprocessor
variable = c
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_e]
type = ElementIntegralVariablePostprocessor
variable = c
block = elyte
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_c]
type = ElementIntegralVariablePostprocessor
variable = c
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[H_a]
type = ADElementIntegralMaterialProperty
mat_prop = enthalpy
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[H_e]
type = ADElementIntegralMaterialProperty
mat_prop = enthalpy
block = elyte
execute_on = 'INITIAL TIMESTEP_END'
[]
[H_c]
type = ADElementIntegralMaterialProperty
mat_prop = enthalpy
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill]
type = Terminator
expression = 'V >= ${Vmax}'
message = 'Voltage reached Vmax'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
line_search = none
nl_rel_tol = 1e-6
nl_abs_tol = 1e-10
nl_max_its = 20
l_max_its = 150
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
optimal_iterations = 7
iteration_window = 2
growth_factor = 1.2
cutback_factor = 0.5
cutback_factor_at_failure = 0.2
linear_iteration_ratio = 1000000
[]
dtmax = 1
end_time = 100000
[]
[Outputs]
exodus = true
print_linear_residuals = false
[]
(test/tests/chemical-electrical-thermal-mechanical/pressure.i)
R = 8.3145 #mJ/mmol/K
F = 96485 #mC/mmol
I = 3e-3 #mA
width = 0.03 #mm
in = '${fparse -I/width}'
t0 = '${fparse -1e-2/in}'
dt = '${fparse t0/100}'
Vmax = 4.3 #V
vf_se = 0.3
vf_cp = 0.5
vf_ca = 0.2
sigma_a = 0.2 #mS/mm
sigma_se = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.2 #mS/mm
sigma_e = ${sigma_se}
sigma_c = '${fparse vf_se*sigma_se + vf_cp*sigma_cp + vf_ca*sigma_ca}'
l0 = 0
l1 = 0.04
l2 = 0.07
l3 = 0.12
cmax = 1e-3 #mmol/mm^3
c0_a = 1e-4
c0_e = 5e-4
c0_c = 1e-3
M_a = 8e-11
M_se = 1e-11
M_cp = 4e-14
M_ca = 1e-13
M_e = ${M_se}
M_c = '${fparse vf_se*M_se + vf_cp*M_cp + vf_ca*M_ca}'
T0 = 300 #K
i0_a = 0.1 #mA/mm^2
i0_c = 0.1 #mA/mm^2
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
hconv = 9.5e-3 #mJ/mm^2/K/s
T_penalty = 2e-1
E_e = 1e4
E_a = 2e5
nu_e = 0.25
nu_a = 0.3
nu_c = 0.3
d = 0
E_e_bar = '${fparse E_e/(1-nu_e*nu_e)}'
E_c_bar = '${fparse (1+d)/(1-d)*E_e_bar}'
E_c = '${fparse (1-nu_c*nu_c)*E_c_bar}'
Omega = 140
beta = 1
CTE = 1e-5
u_penalty = 1e8
P = 10
center = 0.12
spread = 0.01
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) q zeta chi m'
displacements = 'disp_x disp_y'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
[]
[Mesh]
[battery]
type = GeneratedMeshGenerator
dim = 2
xmin = ${l0}
xmax = ${l3}
ymin = 0
ymax = ${width}
nx = 60
ny = 15
[]
[anode]
type = SubdomainBoundingBoxGenerator
input = battery
block_id = 1
block_name = anode
bottom_left = '${l0} 0 0'
top_right = '${l1} ${width} 0'
[]
[elyte]
type = SubdomainBoundingBoxGenerator
input = anode
block_id = 2
block_name = elyte
bottom_left = '${l1} 0 0'
top_right = '${l2} ${width} 0'
[]
[cathode]
type = SubdomainBoundingBoxGenerator
input = elyte
block_id = 3
block_name = cathode
bottom_left = '${l2} 0 0'
top_right = '${l3} ${width} 0'
[]
[anode_elyte]
type = BreakMeshByBlockGenerator
input = cathode
block_pairs = '1 2'
add_interface_on_two_sides = true
split_interface = true
[]
[cathode_elyte]
type = BreakMeshByBlockGenerator
input = anode_elyte
block_pairs = '2 3'
add_interface_on_two_sides = true
split_interface = true
[]
[]
[Variables]
[Phi]
[]
[c]
[]
[T]
initial_condition = ${T0}
[]
[disp_x]
[]
[disp_y]
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[j]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADMaterialRealVectorValueAux
property = j
component = 0
[]
[]
[]
[ICs]
[c_a]
type = ConstantIC
variable = c
value = ${c0_a}
block = 'anode'
[]
[c_e]
type = ConstantIC
variable = c
value = ${c0_e}
block = 'elyte'
[]
[c_c]
type = ConstantIC
variable = c
value = ${c0_c}
block = 'cathode'
[]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${c0_a}
block = 'anode'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c0_e}
block = 'elyte'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${c0_c}
block = 'cathode'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
use_displaced_mesh = true
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
use_displaced_mesh = true
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
use_displaced_mesh = true
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
factor_neighbor = 1
boundary = 'cathode_elyte'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'anode_elyte'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
factor_neighbor = 1
boundary = 'anode_elyte cathode_elyte'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[left]
type = FunctionNeumannBC
variable = Phi
boundary = left
function = in
[]
[right]
type = DirichletBC
variable = Phi
boundary = right
value = 0
[]
[open]
type = OpenBC
variable = c
flux = jm
boundary = 'left right'
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'left right'
value = -1
boundary_material = qconv
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[traction]
enable = false
type = FunctionNeumannBC
variable = disp_y
boundary = 'top'
function = '-if(t<${t0},t/${t0}*${P},${P})/${spread}/sqrt(2*pi)*exp(-0.5*(x-${center})^2/${spread}^2)'
[]
[]
[Constraints]
[y]
enable = false
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'anode ${sigma_a} elyte ${sigma_e} cathode ${sigma_c}'
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
# Migration
[migration]
type = Migration
electrochemical_energy_density = m
electric_potential = Phi
chemical_potential = mu
electric_conductivity = sigma
faraday_constant = ${F}
[]
[migration_flux]
type = MassFlux
mass_flux = jm
energy_densities = 'm'
chemical_potential = mu
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'M'
subdomain_to_prop_value = 'anode ${M_a} elyte ${M_e} cathode ${M_c}'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T
reference_concentration = c_ref
reference_chemical_potential = 1e3
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
property_name = U
expression = 'x:=c/${cmax}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
coupled_variables = c
material_property_names = 'ramp'
boundary = 'anode_elyte'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
property_name = U
expression = 'x:=c/${cmax}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
coupled_variables = c
material_property_names = 'ramp'
boundary = 'cathode_elyte'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'anode_elyte'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cathode_elyte'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
[qconv]
type = ADParsedMaterial
property_name = qconv
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '${hconv} ${T0}'
boundary = 'left right'
[]
# Mechanical
[stiffness_c]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_c*nu_c/(1+nu_c)/(1-2*nu_c)} ${fparse E_c/2/(1+nu_c)}'
block = cathode
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = elyte
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = anode
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
displacements = 'disp_x disp_y'
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
non_swelling_pressure = p
output_properties = 'p'
outputs = exodus
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_r - V_l'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[in]
type = FunctionValuePostprocessor
function = in
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = '-dt*in*${width}'
pp_names = 'dt in'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_a_max]
type = NodalExtremeValue
variable = c
value_type = max
block = anode
[]
[c_c_min]
type = NodalExtremeValue
variable = c
value_type = min
block = cathode
[]
[mass_a]
type = ElementIntegralVariablePostprocessor
variable = c
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_e]
type = ElementIntegralVariablePostprocessor
variable = c
block = elyte
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_c]
type = ElementIntegralVariablePostprocessor
variable = c
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[E_c_bar]
type = ParsedPostprocessor
function = '${E_c_bar}'
execute_on = INITIAL
pp_names = ''
[]
[]
[UserObjects]
[kill]
type = Terminator
expression = 'V >= ${Vmax}'
message = 'Voltage reached Vmax'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
line_search = none
nl_rel_tol = 1e-6
nl_abs_tol = 1e-10
nl_max_its = 20
l_max_its = 150
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
optimal_iterations = 7
iteration_window = 2
growth_factor = 1.2
cutback_factor = 0.5
cutback_factor_at_failure = 0.2
linear_iteration_ratio = 1000000
[]
dtmax = 1
end_time = 100000
[]
[Outputs]
exodus = true
print_linear_residuals = false
[]
(examples/arc-welding/welding.i)
t0 = 10
tm = 30
V = 0.1
end_time = 600
dtmax = 10
dt = 1
W = 40
H = 40
R = 2
q0 = 2.2e4
T_melting = '${fparse 1400+273.15}'
delta_T_pc = 50
L = 2.7e11
kappa = 25
cp = 5e8
rho = 7.8e-9
alpha = 1e-5
E = 1.9e5
nu = 0.27
htc = 5
T_inf = 300
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2
xmax = ${W}
ymax = ${H}
nx = 40
ny = 40
[]
[left]
type = SubdomainBoundingBoxGenerator
input = gmg
block_id = 0
block_name = left
bottom_left = '0 0 0'
top_right = '${fparse W/2} ${H} 0'
[]
[right]
type = SubdomainBoundingBoxGenerator
input = left
block_id = 1
block_name = right
bottom_left = '${fparse W/2} 0 0'
top_right = '${W} ${H} 0'
[]
[break]
type = BreakMeshByBlockGenerator
input = right
[]
use_displaced_mesh = false
[]
[Variables]
[T]
initial_condition = 300
[]
[disp_x]
[]
[disp_y]
[]
[]
[AuxVariables]
[T_ref]
initial_condition = 300
[]
[arc_x]
[AuxKernel]
type = ConstantAux
value = '${fparse W/2}'
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[arc_y]
[AuxKernel]
type = ParsedAux
expression = 'if(t<${tm}, ${H}, ${H}-${V}*(t-${tm}))'
use_xyzt = true
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[r]
[AuxKernel]
type = ParsedAux
coupled_variables = 'arc_x arc_y'
expression = 'sqrt((x-arc_x)^2+(y-arc_y)^2)'
use_xyzt = true
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[]
[Modules]
[TensorMechanics]
[Master]
[all]
strain = FINITE
incremental = true
temperature = T
eigenstrain_names = thermal_eigenstrain
use_automatic_differentiation = true
[]
[]
[CohesiveZoneMaster]
[interface]
strain = FINITE
boundary = interface
use_automatic_differentiation = true
[]
[]
[]
[]
[Kernels]
[htime]
type = ADHeatConductionTimeDerivative
variable = T
density_name = rho
specific_heat = cp
[]
[hcond]
type = ADHeatConduction
variable = T
thermal_conductivity = kappa
[]
[hsource]
type = MaterialSource
variable = T
coefficient = -1
prop = q
[]
[]
[BCs]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = left
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'top bottom left right'
value = -1
boundary_material = qconv
[]
[]
[Functions]
[q_ramp]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${q0}'
[]
[yield]
type = PiecewiseLinear
x = '300 1000 2000'
y = '205 180 150'
[]
[]
[Materials]
[q_ramp]
type = ADGenericFunctionMaterial
prop_names = 'q_ramp'
prop_values = 'q_ramp'
[]
[q]
type = ADParsedMaterial
property_name = q
expression = 'q_ramp*exp(-0.5*(r/R)^2)/R/sqrt(2*3.1415926)'
material_property_names = 'q_ramp'
coupled_variables = 'r'
constant_names = 'R q0'
constant_expressions = '${R} ${q0}'
outputs = exodus
[]
[bulk]
type = ADGenericConstantMaterial
prop_names = 'rho kappa'
prop_values = '${rho} ${kappa}'
[]
[gaussian_function]
type = ADParsedMaterial
property_name = D
expression = 'exp(-T*(T-Tm)^2/dT^2)/sqrt(3.1415926*dT^2)'
coupled_variables = 'T'
constant_names = 'Tm dT'
constant_expressions = '${T_melting} ${delta_T_pc}'
[]
[phi]
type = ADParsedMaterial
property_name = phi
expression = 'if(T<Tm, 0, if(T<Tm+dT, (T-Tm)/dT, 1))'
coupled_variables = 'T'
constant_names = 'Tm dT'
constant_expressions = '${T_melting} ${delta_T_pc}'
outputs = exodus
[]
[cp]
type = ADParsedMaterial
property_name = cp
expression = '${cp} + ${L} * D'
material_property_names = 'D'
[]
[eigenstrain]
type = ADComputeThermalExpansionEigenstrain
eigenstrain_name = thermal_eigenstrain
thermal_expansion_coeff = ${alpha}
stress_free_temperature = T_ref
temperature = T
[]
[C]
type = ADComputeIsotropicElasticityTensor
youngs_modulus = ${E}
poissons_ratio = ${nu}
[]
[plasticity]
type = ADComputeMultipleInelasticStress
inelastic_models = 'hardening'
[]
[hardening]
type = ADIsotropicPlasticityStressUpdate
hardening_constant = 600
yield_stress_function = yield
temperature = T
[]
[interface_stress]
type = WeldedInterfaceTraction
normal_stiffness = '${fparse 100*E}'
tangential_stiffness = '${fparse 60*E}'
phase = phi
phase_history_maximum = phi_max
residual_stiffness = 1e-6
boundary = interface
[]
[qconv]
type = ADParsedMaterial
property_name = qconv
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '${htc} ${T_inf}'
boundary = 'top bottom left right'
[]
[hydrostatic]
type = ADRankTwoInvariant
property_name = hydrostatic_stress
rank_two_tensor = stress
invariant = Hydrostatic
outputs = exodus
[]
[vonmises]
type = ADRankTwoInvariant
property_name = vonmises_stress
rank_two_tensor = stress
invariant = VonMisesStress
outputs = exodus
[]
[plastic_strain]
type = ADRankTwoInvariant
property_name = ep
rank_two_tensor = combined_inelastic_strain
invariant = EffectiveStrain
outputs = exodus
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
end_time = ${end_time}
dtmax = ${dtmax}
dtmin = 0.01
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
growth_factor = 1.2
optimal_iterations = 7
iteration_window = 2
linear_iteration_ratio = 100000
[]
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
nl_abs_tol = 1e-8
nl_rel_tol = 1e-6
nl_max_its = 12
[]
[Outputs]
exodus = true
[]
(examples/PCM-GF/equilibriate.i)
# units are in meter kelvin second (m,kg,s)
end_time = 36000
dtmax = 100
dt = 1
sigma_PCM = 5 # (from Wen's measurement of Gfoam+PCM in radial direction) (from Cfoam 70% dense foam = 28571.43) S/m (1/electrical resistivity (0.000035 ohm-m))
kappa_PCM = 10 #18.8 # (average of Kxy = 14 W/m-K, Kz = 23.6 W/mK at T=700C) #from Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
rho_PCM = 2050 # kg/m^3 #from Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
cp_PCM = 1074 # J/kg-K #from Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
sigma_pipe = 750750.75 # S/m (resistivity 1.332e-6 ohm-m at T = 700C) #Special metal data sheet
kappa_pipe = 23.9 # W/m-K (at 700C) #Special metal datasheet
rho_pipe = 8359.33 #kg/m^3
cp_pipe = 419 # J/kg-K
sigma_gas = 1e-12
kappa_gas = 0.03 #file:///C:/Users/barua/Downloads/PDS-FOAMGLAS%20ONE-US-en.pdf
rho_gas = 1.29 #file:///C:/Users/barua/Downloads/PDS-FOAMGLAS%20ONE-US-en.pdf
cp_gas = 1000 #file:///C:/Users/barua/Downloads/PDS-FOAMGLAS%20ONE-US-en.pdf
[GlobalParams]
energy_densities = 'E H'
[]
[Mesh]
[fmg]
type = FileMeshGenerator
file = 'gold/partial.msh'
[]
coord_type = RZ
[]
[Variables]
[Phi]
[]
[T]
initial_condition = 900
[]
[]
[AuxVariables]
[T_old]
[AuxKernel]
type = ParsedAux
expression = 'T'
coupled_variables = 'T'
execute_on = 'INITIAL TIMESTEP_BEGIN'
[]
[]
[]
[Kernels]
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cp
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[energy_balance_3]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[BCs]
[ground]
type = DirichletBC
variable = Phi
boundary = 'PCM_left'
value = 0
[]
[CV]
type = FunctionDirichletBC
variable = Phi
boundary = 'PCM_right'
function = 0
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'inlet'
value = -1
boundary_material = qconv
[]
[]
[Materials]
[electrical_conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = sigma
subdomain_to_prop_value = 'PCM ${sigma_PCM} pipe ${sigma_pipe} gas ${sigma_gas}'
[]
[charge_trasport]
type = BulkChargeTransport
electrical_energy_density = E
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
[thermal_conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = kappa
subdomain_to_prop_value = 'PCM ${kappa_PCM} pipe ${kappa_pipe} gas ${kappa_gas}'
[]
[density]
type = ADPiecewiseConstantByBlockMaterial
prop_name = rho
subdomain_to_prop_value = 'PCM ${rho_PCM} pipe ${rho_pipe} gas ${rho_gas}'
[]
[specific_heat]
type = ADPiecewiseConstantByBlockMaterial
prop_name = cp
subdomain_to_prop_value = 'PCM ${cp_PCM} pipe ${cp_pipe} gas ${cp_gas}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = H
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
[qconv]
type = ADParsedMaterial
property_name = qconv
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '100 300'
boundary = 'inlet'
[]
[delta_enthalpy]
type = ADParsedMaterial
property_name = delta_enthalpy
expression = 'rho*cp*(T-T_old)/2'
material_property_names = 'rho cp'
coupled_variables = 'T T_old'
[]
[]
[Postprocessors]
[Tin]
type = SideExtremeValue
variable = T
value_type = min
boundary = inlet
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
end_time = ${end_time}
dtmax = ${dtmax}
dtmin = 0.01
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
growth_factor = 1.2
optimal_iterations = 7
iteration_window = 2
linear_iteration_ratio = 100000
[]
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
nl_abs_tol = 1e-8
nl_rel_tol = 1e-6
nl_max_its = 12
[]
[Outputs]
exodus = true
[]
(examples/SSB/charging.i)
I = 1.2e-3 #mA
width = 0.03 #mm
in = '${fparse -I/width}'
t0 = '${fparse -1e-2/in}'
sigma_a = 0.2 #mS/mm
sigma_e = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.5 #mS/mm
sigma_cm = 0.05 #mS/mm
Phi_penalty = 1
cmin_a = 1e-4 #mmol/mm^3
cmax_a = 1e-3 #mmol/mm^3
c_e = 5e-4 #mmol/mm^3
cmin_c = 1e-4 #mmol/mm^3
cmax_c = 1e-3 #mmol/mm^3
c_ref_entropy = 5e-5
D_cp = 5e-5 #mm^2/s
D_cm = 1e-4 #mm^2/s
D_a = 5e-4 #mm^2/s
D_e = 1e-4 #mm^2/s
c_penalty = 1
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-4 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_cp = 6e4
E_cm = 5e4
E_e = 5e4
E_a = 1e5
nu_cp = 0.3
nu_cm = 0.25
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 60
beta = 1e-4
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
T_penalty = 1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
[]
[Mesh]
[battery]
type = FileMeshGenerator
file = 'gold/ssb.msh'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = battery
add_interface_on_two_sides = true
split_interface = true
[]
[]
[Variables]
[Phi_ca]
block = cm
[]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[T]
initial_condition = ${T0}
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[stress]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = pk1
scalar_type = VonMisesStress
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[]
[ICs]
[c_a]
type = ConstantIC
variable = c
value = ${cmin_a}
block = 'a'
[]
[c_e]
type = ConstantIC
variable = c
value = ${c_e}
block = 'cm e'
[]
[c_c]
type = ConstantIC
variable = c
value = ${cmax_c}
block = 'cp'
[]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${cmin_a}
block = 'a'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c_e}
block = 'cm e'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${cmax_c}
block = 'cp'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge_balance_ca]
type = RankOneDivergence
variable = Phi_ca
vector = i_ca
block = cm
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'a_e cm_cp'
[]
[negative_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = 1
boundary = 'a_e cm_cp'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'a_e cm_cp e_a cp_cm'
[]
[continuity_c]
type = InterfaceContinuity
variable = c
neighbor_var = c
penalty = ${c_penalty}
boundary = 'cm_e'
[]
[continuity_Phi_ca]
type = InterfaceContinuity
variable = Phi_ca
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_cp'
[]
[continuity_Phi]
type = InterfaceContinuity
variable = Phi
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_e'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[current]
type = FunctionNeumannBC
variable = Phi
boundary = right
function = in
[]
[electric_potential]
type = DirichletBC
variable = Phi_ca
boundary = left
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'left right'
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'a ${sigma_a} e ${sigma_e} cm ${sigma_cm} cp ${sigma_cp}'
[]
[conductivity_ca]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma_ca'
subdomain_to_prop_value = 'cm ${sigma_ca}'
block = cm
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[charge_transport_ca]
type = BulkChargeTransport
electrical_energy_density = q_ca
electric_potential = Phi_ca
electric_conductivity = sigma_ca
temperature = T
block = cm
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
[current_density_ca]
type = CurrentDensity
current_density = i_ca
electric_potential = Phi_ca
block = cm
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'a ${D_a} e ${D_e} cm ${D_cm} cp ${D_cp}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T_ref
reference_concentration = ${c_ref_entropy}
[]
[diffusion]
type = CondensedMassDiffusion
mass_flux = j
mobility = M
concentration = c
output_properties = 'j'
outputs = exodus
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_a}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
args = c
material_property_names = 'ramp'
block = 'a'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_c}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
args = c
material_property_names = 'ramp'
block = 'cp'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'a_e'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'e_a'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cp_cm'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cm_cp'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
# Mechanical
[stiffness_cp]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cp*nu_cp/(1+nu_cp)/(1-2*nu_cp)} ${fparse E_cp/2/(1+nu_cp)}'
block = cp
[]
[stiffness_cm]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cm*nu_cm/(1+nu_cm)/(1-2*nu_cm)} ${fparse E_cm/2/(1+nu_cm)}'
block = cm
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = e
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = a
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
displacements = 'disp_x disp_y'
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi_ca
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_l - V_r'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[in]
type = FunctionValuePostprocessor
function = in
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = '-dt*in*${width}'
pp_names = 'dt in'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[cmin_c]
type = NodalExtremeValue
variable = c
value_type = min
block = 'cp'
[]
[cmax_a]
type = NodalExtremeValue
variable = c
value_type = max
block = 'a'
[]
[]
[UserObjects]
[kill_a]
type = Terminator
expression = 'cmax_a >= ${cmax_a}'
message = 'Concentration in anode exceeds the maximum allowable value.'
[]
[kill_cp]
type = Terminator
expression = 'cmin_c <= ${cmin_c}'
message = 'Concentration in cathode particle is below the minimum allowable value.'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-ksp_converged_reason'
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_truncfactor'
petsc_options_value = 'hypre boomeramg 151 0.25 ext+i PMIS 4 2 0.4'
# petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
# petsc_options_value = 'lu superlu_dist'
automatic_scaling = true
# line_search = none
ignore_variables_for_autoscaling = 'T'
verbose = true
l_tol = 1e-06
nl_rel_tol = 1e-6
nl_abs_tol = 1e-10
nl_max_its = 12
l_max_its = 150
# [TimeStepper]
# type = IterationAdaptiveDT
# dt = ${dt}
# optimal_iterations = 7
# iteration_window = 2
# growth_factor = 1.2
# cutback_factor = 0.2
# cutback_factor_at_failure = 0.1
# linear_iteration_ratio = 100000
# []
[TimeStepper]
type = FunctionDT
function = 'if(t<${t0}, ${fparse t0/50}, ${fparse -1e-2/in})'
growth_factor = 1.2
cutback_factor_at_failure = 0.1
[]
end_time = 100000
[]
[Outputs]
file_base = 'beta_${beta}'
csv = true
exodus = true
print_linear_residuals = false
[]
(examples/SSB/CV_charging.i)
I = 0.005 #mA
width = 0.05 #mm
in = '${fparse -I/width}'
t0 = '${fparse -1e-2/in}'
sigma_a = 0.2 #mS/mm
sigma_e = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.2 #mS/mm
sigma_cm = 0.05 #mS/mm
Phi_penalty = 10
cmin_a = 1e-4 #mmol/mm^3
cmax_a = 1e-3 #mmol/mm^3
c_e = 5e-4 #mmol/mm^3
cmax_c = 1e-3 #mmol/mm^3
c_ref_entropy = 5e-5
D_cp = 5e-5 #mm^2/s
D_cm = 1e-4 #mm^2/s
D_a = 5e-4 #mm^2/s
D_e = 1e-4 #mm^2/s
c_penalty = 1
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_cp = 6e4
E_cm = 5e4
E_e = 5e4
E_a = 1e5
nu_cp = 0.3
nu_cm = 0.25
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 140
beta = 1e-4
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
htc = 9.5e-3
T_penalty = 1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
displacements = 'disp_x disp_y'
[]
[Problem]
restart_file_base = 'CC_charging_out_cp/LATEST'
[]
[Mesh]
[battery]
type = FileMeshGenerator
file = 'gold/ssb.msh'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = battery
add_interface_on_two_sides = true
split_interface = true
[]
use_displaced_mesh = false
[]
[Variables]
[Phi_ca]
block = cm
[]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[T]
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[stress]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = pk1
scalar_type = Hydrostatic
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Js]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Fs
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Jt]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Ft
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Phi0]
[AuxKernel]
type = ParsedAux
function = 'Phi'
args = 'Phi'
execute_on = 'INITIAL'
[]
[]
[]
[ICs]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${cmin_a}
block = 'a'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c_e}
block = 'cm e'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${cmax_c}
block = 'cp'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge_balance_ca]
type = RankOneDivergence
variable = Phi_ca
vector = i_ca
block = cm
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'a_e cm_cp'
[]
[negative_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = 1
boundary = 'a_e cm_cp'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'a_e cm_cp e_a cp_cm'
[]
[continuity_c]
type = InterfaceContinuity
variable = c
neighbor_var = c
penalty = ${c_penalty}
boundary = 'cm_e'
[]
[continuity_Phi_ca]
type = InterfaceContinuity
variable = Phi_ca
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_cp'
[]
[continuity_Phi]
type = InterfaceContinuity
variable = Phi
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_e'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[current]
type = CoupledVarDirichletBC
variable = Phi
boundary = right
value = Phi0
[]
[potential]
type = DirichletBC
variable = Phi_ca
boundary = left
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'left right'
value = -1
boundary_material = qconv
[]
[]
[Constraints]
[ev_y]
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'a ${sigma_a} e ${sigma_e} cm ${sigma_cm} cp ${sigma_cp}'
[]
[conductivity_ca]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma_ca'
subdomain_to_prop_value = 'cm ${sigma_ca}'
block = cm
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[charge_transport_ca]
type = BulkChargeTransport
electrical_energy_density = q_ca
electric_potential = Phi_ca
electric_conductivity = sigma_ca
temperature = T
block = cm
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
output_properties = i
outputs = exodus
[]
[current_density_ca]
type = CurrentDensity
current_density = i_ca
electric_potential = Phi_ca
output_properties = i_ca
outputs = exodus
block = cm
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'a ${D_a} e ${D_e} cm ${D_cm} cp ${D_cp}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T_ref
reference_concentration = ${c_ref_entropy}
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta m'
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
outputs = exodus
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_a}; 2.77e-4*x^2-0.0069*x+0.0785'
# function = 'x:=c/${cmax_a}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
args = c
material_property_names = 'ramp'
block = 'a'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_c}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
args = c
material_property_names = 'ramp'
block = 'cp'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'a_e'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'e_a'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cp_cm'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cm_cp'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
output_properties = h
outputs = exodus
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
output_properties = r
outputs = exodus
[]
[conv]
type = ADParsedMaterial
f_name = qconv
function = '${htc}*(T-T_ref)'
args = 'T T_ref'
boundary = 'left right'
[]
# Mechanical
[stiffness_cp]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cp*nu_cp/(1+nu_cp)/(1-2*nu_cp)} ${fparse E_cp/2/(1+nu_cp)}'
block = cp
[]
[stiffness_cm]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cm*nu_cm/(1+nu_cm)/(1-2*nu_cm)} ${fparse E_cm/2/(1+nu_cm)}'
block = cm
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = e
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = a
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
non_swelling_pressure = p
output_properties = 'p'
outputs = exodus
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi_ca
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_l - V_r'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[I]
type = ADSideIntegralMaterialProperty
property = i
component = 0
boundary = right
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = 'dt*I'
pp_names = 'dt I'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill_V]
type = Terminator
expression = 'I <= 1e-6'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-ksp_converged_reason'
# petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_truncfactor'
# petsc_options_value = 'hypre boomeramg 301 0.25 ext+i PMIS 4 2 0.4'
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
ignore_variables_for_autoscaling = 'T'
verbose = true
line_search = none
l_max_its = 300
l_tol = 1e-6
nl_rel_tol = 1e-6
nl_abs_tol = 1e-9
nl_max_its = 12
[Predictor]
type = SimplePredictor
scale = 1
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse t0/50}'
optimal_iterations = 6
iteration_window = 1
growth_factor = 1.2
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
linear_iteration_ratio = 100
[]
end_time = 10000
[]
[Outputs]
exodus = true
csv = true
print_linear_residuals = false
checkpoint = true
[]
(examples/SSB_3D/CC_discharging.i)
I = 0.005 #mA
width = 0.05 #mm
in = '${fparse -I/width/width}'
t0 = '${fparse -1e-2/in}'
sigma_a = 0.2 #mS/mm
sigma_e = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.2 #mS/mm
sigma_cm = 0.05 #mS/mm
Phi_penalty = 10
cmin_a = 1e-4 #mmol/mm^3
cmax_a = 1e-3 #mmol/mm^3
c_e = 5e-4 #mmol/mm^3
cmax_c = 1e-3 #mmol/mm^3
c_ref_entropy = 5e-5
D_cp = 5e-5 #mm^2/s
D_cm = 1e-4 #mm^2/s
D_a = 5e-4 #mm^2/s
D_e = 1e-4 #mm^2/s
c_penalty = 1
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_cp = 6e4
E_cm = 5e4
E_e = 5e4
E_a = 1e5
nu_cp = 0.3
nu_cm = 0.25
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 140
beta = 1e-4
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
htc = 9.5e-3
T_penalty = 1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
displacements = 'disp_x disp_y disp_z'
[]
[Problem]
restart_file_base = '${outname}_CV_charging_cp/LATEST'
[]
[Mesh]
[battery]
type = FileMeshGenerator
file = '../coarse.e'
[]
[scale]
type = TransformGenerator
input = battery
transform = SCALE
vector_value = '1e-3 1e-3 1e-3' #um to mm
[]
[cathode_particle]
type = RenameBlockGenerator
input = scale
old_block = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44'
new_block = 'cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp'
[]
[cathode_matrix]
type = RenameBlockGenerator
input = cathode_particle
old_block = '48'
new_block = 'cm'
[]
[elyte]
type = RenameBlockGenerator
input = cathode_matrix
old_block = 49
new_block = 'e'
[]
[anode]
type = RenameBlockGenerator
input = elyte
old_block = 50
new_block = 'a'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = anode
add_interface_on_two_sides = true
split_interface = true
[]
[sidesets]
type = SideSetsFromNormalsGenerator
input = interfaces
normals = '-1 0 0 1 0 0 0 -1 0 0 1 0 0 0 -1 0 0 1'
new_boundary = 'left right bottom top back front'
[]
use_displaced_mesh = false
[]
[Variables]
[Phi_ca]
block = cm
[]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[T]
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[stress]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = pk1
scalar_type = Hydrostatic
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Js]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Fs
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Jt]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Ft
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Phi0]
[]
[]
[ICs]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${cmin_a}
block = 'a'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c_e}
block = 'cm e'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${cmax_c}
block = 'cp'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge_balance_ca]
type = RankOneDivergence
variable = Phi_ca
vector = i_ca
block = cm
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
[momentum_balance_z]
type = RankTwoDivergence
variable = disp_z
component = 2
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'a_e cm_cp'
[]
[negative_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = 1
boundary = 'a_e cm_cp'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'a_e cm_cp e_a cp_cm'
[]
[continuity_c]
type = InterfaceContinuity
variable = c
neighbor_var = c
penalty = ${c_penalty}
boundary = 'cm_e'
[]
[continuity_Phi_ca]
type = InterfaceContinuity
variable = Phi_ca
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_cp'
[]
[continuity_Phi]
type = InterfaceContinuity
variable = Phi
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_e'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_z]
type = InterfaceContinuity
variable = disp_z
neighbor_var = disp_z
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${fparse -in}'
[]
[]
[BCs]
[current]
type = FunctionNeumannBC
variable = Phi
boundary = right
function = in
[]
[potential]
type = DirichletBC
variable = Phi_ca
boundary = left
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[fix_z]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back'
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'left right'
value = -1
boundary_material = qconv
[]
[]
[Constraints]
[ev_y]
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[ev_z]
type = EqualValueBoundaryConstraint
variable = disp_z
penalty = ${u_penalty}
secondary = front
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'a ${sigma_a} e ${sigma_e} cm ${sigma_cm} cp ${sigma_cp}'
[]
[conductivity_ca]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma_ca'
subdomain_to_prop_value = 'cm ${sigma_ca}'
block = cm
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[charge_transport_ca]
type = BulkChargeTransport
electrical_energy_density = q_ca
electric_potential = Phi_ca
electric_conductivity = sigma_ca
temperature = T
block = cm
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
output_properties = i
outputs = exodus
[]
[current_density_ca]
type = CurrentDensity
current_density = i_ca
electric_potential = Phi_ca
output_properties = i_ca
outputs = exodus
block = cm
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'a ${D_a} e ${D_e} cm ${D_cm} cp ${D_cp}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T_ref
reference_concentration = ${c_ref_entropy}
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta m'
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
outputs = exodus
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_a}; 2.77e-4*x^2-0.0069*x+0.0785'
# function = 'x:=c/${cmax_a}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
args = c
material_property_names = 'ramp'
block = 'a'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_c}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
args = c
material_property_names = 'ramp'
block = 'cp'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'a_e'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'e_a'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cp_cm'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cm_cp'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
output_properties = h
outputs = exodus
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
output_properties = r
outputs = exodus
[]
[conv]
type = ADParsedMaterial
f_name = qconv
function = '${htc}*(T-T_ref)'
args = 'T T_ref'
boundary = 'left right'
[]
# Mechanical
[stiffness_cp]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cp*nu_cp/(1+nu_cp)/(1-2*nu_cp)} ${fparse E_cp/2/(1+nu_cp)}'
block = cp
[]
[stiffness_cm]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cm*nu_cm/(1+nu_cm)/(1-2*nu_cm)} ${fparse E_cm/2/(1+nu_cm)}'
block = cm
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = e
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = a
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
non_swelling_pressure = p
output_properties = 'p'
outputs = exodus
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi_ca
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_l - V_r'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[I]
type = ADSideIntegralMaterialProperty
property = i
component = 0
boundary = right
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = 'dt*I'
pp_names = 'dt I'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill_V]
type = Terminator
expression = 'V <= 2.5'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-ksp_converged_reason'
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_truncfactor'
petsc_options_value = 'hypre boomeramg 301 0.25 ext+i PMIS 4 2 0.4'
automatic_scaling = true
ignore_variables_for_autoscaling = 'T'
verbose = true
line_search = none
l_max_its = 300
l_tol = 1e-6
nl_rel_tol = 1e-6
nl_abs_tol = 1e-9
nl_max_its = 12
[Predictor]
type = SimplePredictor
scale = 1
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse t0/50}'
optimal_iterations = 6
iteration_window = 1
growth_factor = 1.2
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
linear_iteration_ratio = 100
[]
end_time = 10000
[]
[Outputs]
file_base = '${outname}_CC_discharging'
exodus = true
csv = true
print_linear_residuals = false
[]
(examples/creep-diffusion/bulk/creep-diffusion.i)
lambda = 1e5
G = 8e4
alpha = 1
Omega = 100
sigma_y = 300
n = 5
A = 1e-6
ad = 1
Nr = 5e-12
Qv = 1e4
c0 = 1e-6
T0 = 800
M = 1e-8
mu0 = 1e3
R = 8.3145
load = 200
t0 = '${fparse load*60}'
dt = '${fparse t0/100}'
tf = 1e9
dtmax = '${fparse tf/100}'
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
strain = E
swelling_strain = Es
mechanical_strain = Em
elastic_strain = Ee
plastic_strain = Ep
equivalent_plastic_strain = ep
eigen_strain = Eg
energy_densities = 'dot(psi_m) dot(psi_c) Delta_p zeta'
volumetric_locking_correction = true
[]
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 3
nx = 3
ny = 3
nz = 3
[]
use_displaced_mesh = false
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[c]
initial_condition = ${c0}
[]
[]
[AuxVariables]
[c_ref]
initial_condition = ${c0}
[]
[T]
initial_condition = ${T0}
[]
[mu_old]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADMaterialRealAux
property = mu
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[]
[Kernels]
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = cauchy
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = cauchy
factor = -1
[]
[momentum_balance_z]
type = RankTwoDivergence
variable = disp_z
component = 2
tensor = cauchy
factor = -1
[]
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
[mass_source]
type = MaterialSource
variable = c
prop = m
coefficient = -1
[]
[]
[Functions]
[load]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${load}'
[]
[]
[BCs]
[fix_x]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0
[]
[fix_y]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0
[]
[fix_z]
type = DirichletBC
variable = disp_z
boundary = 'back'
value = 0
[]
[force_y]
type = FunctionNeumannBC
variable = disp_y
boundary = 'top'
function = load
[]
[]
[Constraints]
[ev_x]
type = EqualValueBoundaryConstraint
variable = disp_x
secondary = right
penalty = 1e8
[]
[ev_y]
type = EqualValueBoundaryConstraint
variable = disp_y
secondary = top
penalty = 1e8
[]
[ev_z]
type = EqualValueBoundaryConstraint
variable = disp_z
secondary = front
penalty = 1e8
[]
[]
[Materials]
[chemical]
type = ADGenericConstantMaterial
prop_names = 'M mu0'
prop_values = '${M} ${mu0}'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T
reference_concentration = c_ref
reference_chemical_potential = mu0
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
outputs = exodus
[]
[mass_source]
type = MassSource
mass_source = m
chemical_potential = mu
outputs = exodus
[]
[stiffness]
type = ADGenericConstantMaterial
prop_names = 'lambda G sigma_y'
prop_values = '${lambda} ${G} ${sigma_y}'
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'alpha Nr'
prop_values = '${alpha} ${Nr}'
[]
[swelling]
type = SwellingStrain
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = alpha
[]
[strain]
type = Strain
[]
[mechanical_strain]
type = MechanicalStrain
[]
[envelope]
type = ADDerivativeParsedMaterial
property_name = Delta_p
expression = 'sigma_y * ep_dot * (1 - exp(-${ad}*mu_old/${R}/T))'
coupled_variables = 'mu_old T'
material_property_names = 'sigma_y ep_dot mu0'
additional_derivative_symbols = 'ep_dot'
derivative_order = 2
compute = false
[]
[vacancy_nucleation]
type = ADParsedMaterial
property_name = dDelta_p/dmu
# expression = 'sigma_y * ep_dot * exp(-${ad}*mu_old/${R}/T) * ${ad} / ${R} / T + if(p>0,1,0) * p * Nr * exp(-${av}*mu_old/${R}/T) * ${av} / ${R} / T'
expression = 'sigma_y * ep_dot * exp(-${ad}*mu_old/${R}/T) * ${ad} / ${R} / T + if(p>0,1,0) * p * Nr * exp(- ${Qv} / ${R} / T)'
coupled_variables = 'mu_old T'
material_property_names = 'sigma_y ep_dot Nr p'
[]
[elastic_energy_density]
type = SDElasticEnergyDensity
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
creep_exponent = ${n}
creep_coefficient = ${A}
plastic_dissipation_material = envelope
plastic_power_density = Delta_p
output_properties = 'ep ddot(psi_m)/ddot(c)'
outputs = 'exodus'
[]
[stress]
type = SDStress
cauchy_stress = cauchy
strain_rate = 'dot(E)'
outputs = exodus
[]
[hydrostatic_stress]
type = ADRankTwoInvariant
property_name = p
rank_two_tensor = cauchy
invariant = Hydrostatic
[]
[strain_yy]
type = ADRankTwoCartesianComponent
property_name = Eyy
rank_two_tensor = E
index_i = 1
index_j = 1
[]
[strain_rate_yy]
type = ADRankTwoCartesianComponent
property_name = dot(Eyy)
rank_two_tensor = dot(E)
index_i = 1
index_j = 1
[]
[creep_strain_yy]
type = ADRankTwoCartesianComponent
property_name = Epyy
rank_two_tensor = Ep
index_i = 1
index_j = 1
[]
[swelling_strain_yy]
type = ADRankTwoCartesianComponent
property_name = Esyy
rank_two_tensor = Es
index_i = 1
index_j = 1
[]
[]
[Postprocessors]
[Eyy_dot]
type = ADElementAverageMaterialProperty
mat_prop = dot(Eyy)
execute_on = 'INITIAL TIMESTEP_END'
[]
[Esyy]
type = ADElementAverageMaterialProperty
mat_prop = Esyy
execute_on = 'INITIAL TIMESTEP_END'
outputs = none
[]
[Epyy]
type = ADElementAverageMaterialProperty
mat_prop = Epyy
execute_on = 'INITIAL TIMESTEP_END'
outputs = none
[]
[dt]
type = TimestepSize
execute_on = 'INITIAL TIMESTEP_END'
outputs = none
[]
[delta_Esyy]
type = ChangeOverTimestepPostprocessor
postprocessor = Esyy
execute_on = 'INITIAL TIMESTEP_END'
outputs = none
[]
[delta_Epyy]
type = ChangeOverTimestepPostprocessor
postprocessor = Epyy
execute_on = 'INITIAL TIMESTEP_END'
outputs = none
[]
[Esyy_dot]
type = ParsedPostprocessor
function = 'delta_Esyy/dt'
pp_names = 'delta_Esyy dt'
execute_on = 'INITIAL TIMESTEP_END'
[]
[Epyy_dot]
type = ParsedPostprocessor
function = 'delta_Epyy/dt'
pp_names = 'delta_Epyy dt'
execute_on = 'INITIAL TIMESTEP_END'
[]
[cavity_density]
type = ElementAverageValue
variable = c
execute_on = 'INITIAL TIMESTEP_END'
[]
[cavity_potential]
type = ADElementAverageMaterialProperty
mat_prop = mu
execute_on = 'INITIAL TIMESTEP_END'
[]
[p]
type = ADElementAverageMaterialProperty
mat_prop = p
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
line_search = none
nl_rel_tol = 1e-6
nl_abs_tol = 1e-10
nl_max_its = 20
nl_forced_its = 1
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
optimal_iterations = 6
iteration_window = 1
growth_factor = 2
cutback_factor = 0.5
cutback_factor_at_failure = 0.2
linear_iteration_ratio = 100000
[]
end_time = ${tf}
dtmax = ${dtmax}
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
[]
[Outputs]
file_base = 'load_${load}'
csv = true
exodus = true
[]
(examples/cable/cable.i)
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
in = '${fparse 1e3/pi/34^2}'
sigma_con = 1
c_ref = 1e-3
c_ref_ent = 1e-6
M_con = 1e-1
kappa_con = 398 #mJ/mm/K/s
kappa_air = 3.98 #mJ/mm/K/s
kappa_ins = 0.29 #mJ/mm/K/s
kappa_jac = 0.39 #mJ/mm/K/s
E_con = 1.1e5
E_air = 1e2
E_ins = 6e2
E_jac = 3e2
nu_con = 0.32
nu_air = 0.25
nu_ins = 0.4
nu_jac = 0.3
CTE = 1e-5
htc = 1e1
[GlobalParams]
energy_densities = 'q dot(psi_c) zeta m chi dot(psi_m)'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
thermal_deformation_gradient = Ft
[]
[Problem]
kernel_coverage_check = false
material_coverage_check = false
[]
[Mesh]
[fmg]
type = FileMeshGenerator
file = 'gold/cable.msh'
[]
use_displaced_mesh = false
[]
[Variables]
[Phi]
block = 'conductor'
[]
[c]
initial_condition = ${c_ref}
block = 'conductor'
[]
[T]
initial_condition = ${T0}
[]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]
[AuxVariables]
[stress]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = pk1
scalar_type = Hydrostatic
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[]
[Functions]
[current_density]
type = PiecewiseLinear
x = '0 100'
y = '0 ${in}'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
block = 'conductor'
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
block = 'conductor'
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
block = 'conductor'
[]
# Energy balance
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
[momentum_balance_z]
type = RankTwoDivergence
variable = disp_z
component = 2
tensor = pk1
factor = -1
[]
[]
[BCs]
[current]
type = FunctionNeumannBC
variable = Phi
boundary = 'conductor_top'
function = current_density
[]
[Phi_ref]
type = DirichletBC
variable = Phi
boundary = 'conductor_bottom'
value = 0
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'outer'
value = -1
boundary_material = qconv
[]
[fix_x]
type = DirichletBC
variable = disp_x
boundary = 'center piny'
value = 0
[]
[fix_y]
type = DirichletBC
variable = disp_y
boundary = 'center pinx'
value = 0
[]
[fix_z]
type = DirichletBC
variable = disp_z
boundary = 'conductor_bottom air_bottom insulator_bottom jacket_bottom'
value = 0
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'conductor ${sigma_con}'
block = 'conductor'
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
block = 'conductor'
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
block = 'conductor'
[]
# Chemistry
[mobility]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'M'
subdomain_to_prop_value = 'conductor ${M_con}'
block = 'conductor'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T
reference_concentration = ${c_ref_ent}
block = 'conductor'
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
block = 'conductor'
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
block = 'conductor'
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
block = 'conductor'
[]
# Migration
[migration]
type = Migration
electrochemical_energy_density = m
electric_potential = Phi
chemical_potential = mu
electric_conductivity = sigma
faraday_constant = ${F}
block = 'conductor'
[]
[migration_flux]
type = MassFlux
mass_flux = jm
energy_densities = 'm'
chemical_potential = mu
block = 'conductor'
[]
# Thermal
[heat_conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'kappa'
subdomain_to_prop_value = 'conductor ${kappa_con} air ${kappa_air} insulator ${kappa_ins} jacket ${kappa_jac}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
[qconv]
type = ADParsedMaterial
f_name = qconv
function = 'htc*(T-T_inf)'
args = 'T'
constant_names = 'htc T_inf'
constant_expressions = '${htc} ${T0}'
boundary = 'outer'
[]
# Mechanics
[youngs_modulus]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'lambda'
subdomain_to_prop_value = 'conductor ${fparse E_con*nu_con/(1+nu_con)/(1-2*nu_con)} air ${fparse E_air*nu_air/(1+nu_air)/(1-2*nu_air)} insulator ${fparse E_ins*nu_ins/(1+nu_ins)/(1-2*nu_ins)} jacket ${fparse E_jac*nu_jac/(1+nu_jac)/(1-2*nu_jac)}'
[]
[poissons_ratio]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'G'
subdomain_to_prop_value = 'conductor ${fparse E_con/2/(1+nu_con)} air ${fparse E_air/2/(1+nu_air)} insulator ${fparse E_ins/2/(1+nu_ins)} jacket ${fparse E_jac/2/(1+nu_jac)}'
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = ${T0}
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
displacements = 'disp_x disp_y disp_z'
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[T_max]
type = NodalExtremeValue
variable = T
block = 'conductor'
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-ksp_converged_reason'
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
# petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_truncfactor'
# petsc_options_value = 'hypre boomeramg 301 0.25 ext+i PMIS 4 2 0.4'
automatic_scaling = true
ignore_variables_for_autoscaling = 'T c'
verbose = true
line_search = none
l_max_its = 300
l_tol = 1e-6
nl_rel_tol = 1e-6
nl_abs_tol = 1e-8
nl_max_its = 12
[Predictor]
type = SimplePredictor
scale = 1
[]
[TimeStepper]
type = FunctionDT
function = 'if(t<100,10,100)'
[]
end_time = 100
[]
[Outputs]
exodus = true
print_linear_residuals = false
[]
(examples/SSB_3D/CV_charging.i)
I = 0.005 #mA
width = 0.05 #mm
in = '${fparse -I/width/width}'
t0 = '${fparse -1e-2/in}'
sigma_a = 0.2 #mS/mm
sigma_e = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.2 #mS/mm
sigma_cm = 0.05 #mS/mm
Phi_penalty = 10
cmin_a = 1e-4 #mmol/mm^3
cmax_a = 1e-3 #mmol/mm^3
c_e = 5e-4 #mmol/mm^3
cmax_c = 1e-3 #mmol/mm^3
c_ref_entropy = 5e-5
D_cp = 5e-5 #mm^2/s
D_cm = 1e-4 #mm^2/s
D_a = 5e-4 #mm^2/s
D_e = 1e-4 #mm^2/s
c_penalty = 1
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_cp = 6e4
E_cm = 5e4
E_e = 5e4
E_a = 1e5
nu_cp = 0.3
nu_cm = 0.25
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 140
beta = 1e-4
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
htc = 9.5e-3
T_penalty = 1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
displacements = 'disp_x disp_y disp_z'
[]
[Problem]
restart_file_base = '${outname}_CC_charging_cp/LATEST'
[]
[Mesh]
[battery]
type = FileMeshGenerator
file = '../coarse.e'
[]
[scale]
type = TransformGenerator
input = battery
transform = SCALE
vector_value = '1e-3 1e-3 1e-3' #um to mm
[]
[cathode_particle]
type = RenameBlockGenerator
input = scale
old_block = '1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44'
new_block = 'cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp cp'
[]
[cathode_matrix]
type = RenameBlockGenerator
input = cathode_particle
old_block = '48'
new_block = 'cm'
[]
[elyte]
type = RenameBlockGenerator
input = cathode_matrix
old_block = 49
new_block = 'e'
[]
[anode]
type = RenameBlockGenerator
input = elyte
old_block = 50
new_block = 'a'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = anode
add_interface_on_two_sides = true
split_interface = true
[]
[sidesets]
type = SideSetsFromNormalsGenerator
input = interfaces
normals = '-1 0 0 1 0 0 0 -1 0 0 1 0 0 0 -1 0 0 1'
new_boundary = 'left right bottom top back front'
[]
use_displaced_mesh = false
[]
[Variables]
[Phi_ca]
block = cm
[]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[T]
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[stress]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = pk1
scalar_type = Hydrostatic
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Js]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Fs
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Jt]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Ft
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Phi0]
[AuxKernel]
type = ParsedAux
function = 'Phi'
args = 'Phi'
execute_on = 'INITIAL'
[]
[]
[]
[ICs]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${cmin_a}
block = 'a'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c_e}
block = 'cm e'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${cmax_c}
block = 'cp'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge_balance_ca]
type = RankOneDivergence
variable = Phi_ca
vector = i_ca
block = cm
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
[momentum_balance_z]
type = RankTwoDivergence
variable = disp_z
component = 2
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'a_e cm_cp'
[]
[negative_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = 1
boundary = 'a_e cm_cp'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'a_e cm_cp e_a cp_cm'
[]
[continuity_c]
type = InterfaceContinuity
variable = c
neighbor_var = c
penalty = ${c_penalty}
boundary = 'cm_e'
[]
[continuity_Phi_ca]
type = InterfaceContinuity
variable = Phi_ca
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_cp'
[]
[continuity_Phi]
type = InterfaceContinuity
variable = Phi
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_e'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_z]
type = InterfaceContinuity
variable = disp_z
neighbor_var = disp_z
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[current]
type = CoupledVarDirichletBC
variable = Phi
boundary = right
value = Phi0
[]
[potential]
type = DirichletBC
variable = Phi_ca
boundary = left
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[fix_z]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back'
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'left right'
value = -1
boundary_material = qconv
[]
[]
[Constraints]
[ev_y]
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[ev_z]
type = EqualValueBoundaryConstraint
variable = disp_z
penalty = ${u_penalty}
secondary = front
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'a ${sigma_a} e ${sigma_e} cm ${sigma_cm} cp ${sigma_cp}'
[]
[conductivity_ca]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma_ca'
subdomain_to_prop_value = 'cm ${sigma_ca}'
block = cm
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[charge_transport_ca]
type = BulkChargeTransport
electrical_energy_density = q_ca
electric_potential = Phi_ca
electric_conductivity = sigma_ca
temperature = T
block = cm
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
output_properties = i
outputs = exodus
[]
[current_density_ca]
type = CurrentDensity
current_density = i_ca
electric_potential = Phi_ca
output_properties = i_ca
outputs = exodus
block = cm
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'a ${D_a} e ${D_e} cm ${D_cm} cp ${D_cp}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T_ref
reference_concentration = ${c_ref_entropy}
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta m'
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
outputs = exodus
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_a}; 2.77e-4*x^2-0.0069*x+0.0785'
# function = 'x:=c/${cmax_a}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
args = c
material_property_names = 'ramp'
block = 'a'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_c}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
args = c
material_property_names = 'ramp'
block = 'cp'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'a_e'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'e_a'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cp_cm'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cm_cp'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
output_properties = h
outputs = exodus
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
output_properties = r
outputs = exodus
[]
[conv]
type = ADParsedMaterial
f_name = qconv
function = '${htc}*(T-T_ref)'
args = 'T T_ref'
boundary = 'left right'
[]
# Mechanical
[stiffness_cp]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cp*nu_cp/(1+nu_cp)/(1-2*nu_cp)} ${fparse E_cp/2/(1+nu_cp)}'
block = cp
[]
[stiffness_cm]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cm*nu_cm/(1+nu_cm)/(1-2*nu_cm)} ${fparse E_cm/2/(1+nu_cm)}'
block = cm
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = e
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = a
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
non_swelling_pressure = p
output_properties = 'p'
outputs = exodus
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi_ca
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_l - V_r'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[I]
type = ADSideIntegralMaterialProperty
property = i
component = 0
boundary = right
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = 'dt*I'
pp_names = 'dt I'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill_V]
type = Terminator
expression = 'I <= 1e-6'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-ksp_converged_reason'
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_truncfactor'
petsc_options_value = 'hypre boomeramg 301 0.25 ext+i PMIS 4 2 0.4'
automatic_scaling = true
ignore_variables_for_autoscaling = 'T'
verbose = true
line_search = none
l_max_its = 300
l_tol = 1e-6
nl_rel_tol = 1e-6
nl_abs_tol = 1e-9
nl_max_its = 12
[Predictor]
type = SimplePredictor
scale = 1
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse t0/50}'
optimal_iterations = 6
iteration_window = 1
growth_factor = 1.2
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
linear_iteration_ratio = 100
[]
end_time = 10000
[]
[Outputs]
file_base = '${outname}_CV_charging'
exodus = true
csv = true
print_linear_residuals = false
checkpoint = true
[]
(examples/LiB/CC_charging.i)
I = 3e-3 #mA
width = 0.03 #mm
in = '${fparse -I/width}'
t0 = '${fparse -1e-2/in}'
dt = '${fparse t0/100}'
sigma_a = 1e0 #mS/mm
sigma_e = 1e-1 #mS/mm
sigma_c = 1e-2 #mS/mm
l0 = 0
l1 = 0.04
l2 = 0.07
l3 = 0.12
cmin = 1e-4 #mmol/mm^3
cmax = 1e-3 #mmol/mm^3
D_a = 1e-3 #mm^2/s
D_e = 1e-4 #mm^2/s
D_c = 5e-5 #mm^2/s
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_c = 1e5
E_e = 1e4
E_a = 2e5
nu_c = 0.3
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 60
beta = 1e-3
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
T_penalty = 2e-1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q zeta m'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
[]
[Mesh]
[battery]
type = GeneratedMeshGenerator
dim = 2
xmin = ${l0}
xmax = ${l3}
ymin = 0
ymax = ${width}
nx = 60
ny = 15
[]
[anode]
type = SubdomainBoundingBoxGenerator
input = battery
block_id = 1
block_name = anode
bottom_left = '${l0} 0 0'
top_right = '${l1} ${width} 0'
[]
[elyte]
type = SubdomainBoundingBoxGenerator
input = anode
block_id = 2
block_name = elyte
bottom_left = '${l1} 0 0'
top_right = '${l2} ${width} 0'
[]
[cathode]
type = SubdomainBoundingBoxGenerator
input = elyte
block_id = 3
block_name = cathode
bottom_left = '${l2} 0 0'
top_right = '${l3} ${width} 0'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = cathode
add_interface_on_two_sides = true
split_interface = true
[]
[]
[Variables]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[T]
initial_condition = ${T0}
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[]
[ICs]
[c_min]
type = ConstantIC
variable = c
value = ${cmin}
block = 'anode'
[]
[c_mid]
type = ConstantIC
variable = c
value = '${fparse (cmax+cmin)/2}'
block = 'elyte'
[]
[c_max]
type = ConstantIC
variable = c
value = ${cmax}
block = 'cathode'
[]
[c_ref_min]
type = ConstantIC
variable = c_ref
value = ${cmin}
block = 'anode'
[]
[c_ref_mid]
type = ConstantIC
variable = c_ref
value = '${fparse (cmax+cmin)/2}'
block = 'elyte'
[]
[c_ref_max]
type = ConstantIC
variable = c_ref
value = ${cmax}
block = 'cathode'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'elyte_anode cathode_elyte'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'anode_elyte elyte_cathode'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'anode_elyte elyte_cathode elyte_anode cathode_elyte'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[left]
type = FunctionNeumannBC
variable = Phi
boundary = left
function = in
[]
[right]
type = DirichletBC
variable = Phi
boundary = right
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[open]
type = OpenBC
variable = c
flux = jm
boundary = 'left right'
[]
[]
[Constraints]
[y]
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'anode ${sigma_a} elyte ${sigma_e} cathode ${sigma_c}'
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'anode ${D_a} elyte ${D_e} cathode ${D_c}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T
reference_concentration = c_ref
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
[]
# Migration
[migration]
type = Migration
electrochemical_energy_density = m
electric_potential = Phi
chemical_potential = mu
electric_conductivity = sigma
faraday_constant = ${F}
[]
[migration_flux]
type = MassFlux
mass_flux = jm
energy_densities = 'm'
chemical_potential = mu
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
args = c
material_property_names = 'ramp'
block = 'anode'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
args = c
material_property_names = 'ramp'
block = 'cathode'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'anode_elyte'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'elyte_anode'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cathode_elyte'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'elyte_cathode'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
# Mechanical
[stiffness_c]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_c*nu_c/(1+nu_c)/(1-2*nu_c)} ${fparse E_c/2/(1+nu_c)}'
block = cathode
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = elyte
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = anode
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
displacements = 'disp_x disp_y'
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_r - V_l'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[I]
type = ADSideIntegralMaterialProperty
property = i
component = 0
boundary = right
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = '-dt*I'
pp_names = 'dt I'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_a_max]
type = NodalExtremeValue
variable = c
value_type = max
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_c_min]
type = NodalExtremeValue
variable = c
value_type = min
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_a_min]
type = NodalExtremeValue
variable = c
value_type = min
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_c_max]
type = NodalExtremeValue
variable = c
value_type = max
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_a]
type = ElementIntegralVariablePostprocessor
variable = c
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_e]
type = ElementIntegralVariablePostprocessor
variable = c
block = elyte
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_c]
type = ElementIntegralVariablePostprocessor
variable = c
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill_a]
type = Terminator
expression = 'c_a_max >= ${cmax}'
message = 'Concentration in anode exceeds the maximum allowable value.'
[]
[kill_c]
type = Terminator
expression = 'c_c_min <= ${cmin}'
message = 'Concentration in cathode is below the minimum allowable value.'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
line_search = none
nl_rel_tol = 1e-6
nl_abs_tol = 1e-10
nl_max_its = 20
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
optimal_iterations = 7
iteration_window = 2
growth_factor = 1.2
cutback_factor = 0.5
cutback_factor_at_failure = 0.2
linear_iteration_ratio = 1000000
[]
end_time = 100000
[]
[Outputs]
file_base = 'CC_charging_I_${I}'
csv = true
exodus = true
print_linear_residuals = false
checkpoint = true
[]
(examples/LiB/CV_charging.i)
I = 3e-3 #mA
width = 0.03 #mm
in = '${fparse -I/width}'
t0 = '${fparse -1e-2/in}'
dt = '${fparse t0/100}'
sigma_a = 1e0 #mS/mm
sigma_e = 1e-1 #mS/mm
sigma_c = 1e-2 #mS/mm
l0 = 0
l1 = 0.04
l2 = 0.07
l3 = 0.12
cmin = 1e-4 #mmol/mm^3
cmax = 1e-3 #mmol/mm^3
D_a = 1e-3 #mm^2/s
D_e = 1e-4 #mm^2/s
D_c = 5e-5 #mm^2/s
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_c = 1e5
E_e = 1e4
E_a = 2e5
nu_c = 0.3
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 60
beta = 1e-3
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
T_penalty = 2e-1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q zeta m'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
[]
[Problem]
restart_file_base = 'CC_charging_I_${I}_cp/LATEST'
[]
[Mesh]
[battery]
type = GeneratedMeshGenerator
dim = 2
xmin = ${l0}
xmax = ${l3}
ymin = 0
ymax = ${width}
nx = 60
ny = 15
[]
[anode]
type = SubdomainBoundingBoxGenerator
input = battery
block_id = 1
block_name = anode
bottom_left = '${l0} 0 0'
top_right = '${l1} ${width} 0'
[]
[elyte]
type = SubdomainBoundingBoxGenerator
input = anode
block_id = 2
block_name = elyte
bottom_left = '${l1} 0 0'
top_right = '${l2} ${width} 0'
[]
[cathode]
type = SubdomainBoundingBoxGenerator
input = elyte
block_id = 3
block_name = cathode
bottom_left = '${l2} 0 0'
top_right = '${l3} ${width} 0'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = cathode
add_interface_on_two_sides = true
split_interface = true
[]
[]
[Variables]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[T]
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[Phi0]
[]
[]
[ICs]
[c_ref_min]
type = ConstantIC
variable = c_ref
value = ${cmin}
block = 'anode'
[]
[c_ref_mid]
type = ConstantIC
variable = c_ref
value = '${fparse (cmax+cmin)/2}'
block = 'elyte'
[]
[c_ref_max]
type = ConstantIC
variable = c_ref
value = ${cmax}
block = 'cathode'
[]
[]
[AuxKernels]
[Phi0]
type = ParsedAux
variable = Phi0
function = 'Phi'
args = 'Phi'
execute_on = 'INITIAL'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'elyte_anode cathode_elyte'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'anode_elyte elyte_cathode'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'anode_elyte elyte_cathode elyte_anode cathode_elyte'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'anode_elyte elyte_cathode'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[]
[BCs]
[Phi]
type = CoupledVarDirichletBC
variable = Phi
boundary = 'left right'
value = Phi0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[open]
type = OpenBC
variable = c
flux = jm
boundary = 'left right'
[]
[]
[Constraints]
[y]
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'anode ${sigma_a} elyte ${sigma_e} cathode ${sigma_c}'
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'anode ${D_a} elyte ${D_e} cathode ${D_c}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T
reference_concentration = c_ref
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
[]
# Migration
[migration]
type = Migration
electrochemical_energy_density = m
electric_potential = Phi
chemical_potential = mu
electric_conductivity = sigma
faraday_constant = ${F}
[]
[migration_flux]
type = MassFlux
mass_flux = jm
energy_densities = 'm'
chemical_potential = mu
[]
# Redox
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)'
args = c
block = 'anode'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)'
args = c
block = 'cathode'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'anode_elyte'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'elyte_anode'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cathode_elyte'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'elyte_cathode'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
# Mechanical
[stiffness_c]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_c*nu_c/(1+nu_c)/(1-2*nu_c)} ${fparse E_c/2/(1+nu_c)}'
block = cathode
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = elyte
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = anode
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
displacements = 'disp_x disp_y'
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_r - V_l'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[I]
type = ADSideIntegralMaterialProperty
property = i
component = 0
boundary = right
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = '-dt*I'
pp_names = 'dt I'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_a_max]
type = NodalExtremeValue
variable = c
value_type = max
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_c_min]
type = NodalExtremeValue
variable = c
value_type = min
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_a_min]
type = NodalExtremeValue
variable = c
value_type = min
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_c_max]
type = NodalExtremeValue
variable = c
value_type = max
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_a]
type = ElementIntegralVariablePostprocessor
variable = c
block = anode
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_e]
type = ElementIntegralVariablePostprocessor
variable = c
block = elyte
execute_on = 'INITIAL TIMESTEP_END'
[]
[mass_c]
type = ElementIntegralVariablePostprocessor
variable = c
block = cathode
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill_i]
type = Terminator
expression = '-I <= 1e-6'
message = 'No current.'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
nl_rel_tol = 1e-6
nl_abs_tol = 1e-10
nl_max_its = 20
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
optimal_iterations = 7
iteration_window = 2
growth_factor = 1.2
cutback_factor = 0.5
cutback_factor_at_failure = 0.2
linear_iteration_ratio = 1000000
[]
start_time = 0
end_time = 100000
steady_state_start_time = 10
steady_state_detection = true
[]
[Outputs]
file_base = 'CV_charging_I_${I}'
csv = true
exodus = true
print_linear_residuals = false
checkpoint = true
[]
(test/tests/electrical-thermal/jh.i)
n = 32
[GlobalParams]
energy_densities = 'E H'
[]
[Mesh]
[battery]
type = GeneratedMeshGenerator
dim = 1
xmin = 0
xmax = 1
nx = ${n}
[]
[]
[Variables]
[Phi]
[]
[T]
[]
[]
[Functions]
[Phi]
type = ParsedFunction
expression = '3.565*exp(x)-4.71*x^2+7.1*x-3.4'
[]
[sigma]
type = ParsedFunction
expression = '1+x'
[]
[q]
type = ParsedFunction
expression = '2.32+exp(x)*(-7.13-3.565*x)+18.84*x'
[]
[kappa]
type = ParsedFunction
expression = '-2.6*x+3'
[]
[]
[Kernels]
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge]
type = BodyForce
variable = Phi
function = q
[]
[energy_balance_1]
type = RankOneDivergence
variable = T
vector = h
[]
[energy_balance_2]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[BCs]
[fix]
type = FunctionDirichletBC
variable = Phi
boundary = 'left right'
function = Phi
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = right
value = -1
boundary_material = qconv
[]
[]
[Materials]
[electric_constants]
type = ADGenericFunctionMaterial
prop_names = 'sigma'
prop_values = 'sigma'
[]
[charge_trasport]
type = BulkChargeTransport
electrical_energy_density = E
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
[thermal_constants]
type = ADGenericFunctionMaterial
prop_names = 'kappa'
prop_values = 'kappa'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = H
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
[qconv]
type = ADParsedMaterial
property_name = qconv
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '1.35 300'
boundary = right
[]
[]
[Postprocessors]
[error]
type = ElementL2Error
variable = Phi
function = Phi
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
num_steps = 1
[]
[Outputs]
exodus = true
[]
(examples/resin-printing/print.i)
W = 100
H = 80
Ht = 70
nw = 100
nh = 80
nlayer = 70
hlayer = '${fparse Ht/nlayer}'
tlayer = 10
end_time = '${fparse nlayer*tlayer}'
dtmax = 2
dt = 2
# FEP properties: I am using those of stainless steel which are apparently wrong, but whatever
kappa_fep = 25
cp_fep = 5e8
rho_fep = 7.8e-9
alpha_fep = 1e-5
E_fep = 1.9e5
nu_fep = 0.27
# liquid properties: Polyurethane resin
kappa_liquid = 0.2
cp_liquid = 1.8e8
rho_liquid = 1e-9
alpha_liquid = 0
E_liquid = 100
nu_liquid = 0.27
# solid properties: Polyurethane
kappa_solid = 0.02
cp_solid = 2.4e9
rho_solid = 5e-10
alpha_solid = 2e-4
E_solid = 100
nu_solid = 0.27
q = 0.2
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 2
xmax = ${W}
ymax = ${H}
nx = ${nw}
ny = ${nh}
[]
[FEP]
type = SubdomainBoundingBoxGenerator
input = gmg
block_id = 0
block_name = FEP
bottom_left = '0 ${Ht} 0'
top_right = '${W} ${H} 0'
[]
[solid]
type = SubdomainBoundingBoxGenerator
input = FEP
block_id = 1
block_name = solid
bottom_left = '0 ${fparse Ht/2} 0'
top_right = '${W} ${Ht} 0'
[]
[liquid]
type = SubdomainBoundingBoxGenerator
input = solid
block_id = 2
block_name = liquid
bottom_left = '0 0 0'
top_right = '${W} ${fparse Ht/2} 0'
[]
use_displaced_mesh = false
[]
[Variables]
[T]
initial_condition = 300
[]
[disp_x]
[]
[disp_y]
[]
[]
[AuxVariables]
[T_ref]
initial_condition = 300
[]
[bunny]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = FunctionAux
function = bunny_current
block = 'liquid solid'
execute_on = 'INITIAL TIMESTEP_BEGIN'
[]
[]
[bunny_layer]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = FunctionAux
function = bunny_layer
block = 'liquid solid'
execute_on = 'INITIAL TIMESTEP_BEGIN'
[]
[]
[]
[Functions]
[bunny]
type = ImageFunction
file = gold/bunny.png
dimensions = '${W} ${Ht} 0'
flip_y = true
lower_value = 1
upper_value = 0
threshold = 100
[]
[current]
type = ParsedFunction
expression = 'i:=ceil(t/${tlayer}); if(y>${Ht}-${hlayer}*i, 1, 0)'
[]
[layer]
type = ParsedFunction
expression = 'i:=ceil(t/${tlayer}); if(y>${Ht}-${hlayer}*i, if(y<${Ht}-${hlayer}*(i-1), 1, 0), 0)'
[]
[bunny_current]
type = CompositeFunction
functions = 'bunny current'
[]
[bunny_layer]
type = CompositeFunction
functions = 'bunny layer'
[]
[]
[UserObjects]
[esm]
type = CoupledVarThresholdElementSubdomainModifier
coupled_var = bunny
criterion_type = ABOVE
threshold = 1e-3
subdomain_id = 1
complement_subdomain_id = 2
apply_initial_conditions = false
block = 'liquid solid'
[]
[]
[Modules]
[TensorMechanics]
[Master]
[all]
strain = FINITE
incremental = true
temperature = T
eigenstrain_names = thermal_eigenstrain
use_automatic_differentiation = true
[]
[]
[]
[]
[Kernels]
[htime]
type = ADHeatConductionTimeDerivative
variable = T
density_name = rho
specific_heat = cp
[]
[hcond]
type = ADHeatConduction
variable = T
thermal_conductivity = kappa
[]
[hsource]
type = MaterialSource
variable = T
coefficient = -1
prop = q
block = 'liquid solid'
[]
[]
[BCs]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = top
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = top
[]
# [hconv]
# type = ADMatNeumannBC
# variable = T
# boundary = 'top bottom left right'
# value = -1
# boundary_material = qconv
# []
[]
[Materials]
[q]
type = ADParsedMaterial
property_name = q
coupled_variables = 'bunny_layer'
expression = 'bunny_layer * ${q}'
block = 'liquid solid'
[]
[kappa]
type = ADPiecewiseConstantByBlockMaterial
prop_name = kappa
subdomain_to_prop_value = 'FEP ${kappa_fep} liquid ${kappa_liquid} solid ${kappa_solid}'
[]
[cp]
type = ADPiecewiseConstantByBlockMaterial
prop_name = cp
subdomain_to_prop_value = 'FEP ${cp_fep} liquid ${cp_liquid} solid ${cp_solid}'
[]
[rho]
type = ADPiecewiseConstantByBlockMaterial
prop_name = rho
subdomain_to_prop_value = 'FEP ${rho_fep} liquid ${rho_liquid} solid ${rho_solid}'
[]
[eigenstrain_FEP]
type = ADComputeThermalExpansionEigenstrain
eigenstrain_name = thermal_eigenstrain
thermal_expansion_coeff = ${alpha_fep}
stress_free_temperature = T_ref
temperature = T
block = FEP
[]
[eigenstrain_liquid]
type = ADComputeThermalExpansionEigenstrain
eigenstrain_name = thermal_eigenstrain
thermal_expansion_coeff = ${alpha_liquid}
stress_free_temperature = T_ref
temperature = T
block = liquid
[]
[eigenstrain_solid]
type = ADComputeThermalExpansionEigenstrain
eigenstrain_name = thermal_eigenstrain
thermal_expansion_coeff = ${alpha_solid}
stress_free_temperature = T_ref
temperature = T
block = solid
[]
[C_FEP]
type = ADComputeIsotropicElasticityTensor
youngs_modulus = ${E_fep}
poissons_ratio = ${nu_fep}
block = FEP
[]
[C_liquid]
type = ADComputeIsotropicElasticityTensor
youngs_modulus = ${E_liquid}
poissons_ratio = ${nu_liquid}
block = liquid
[]
[C_solid]
type = ADComputeIsotropicElasticityTensor
youngs_modulus = ${E_solid}
poissons_ratio = ${nu_solid}
block = solid
[]
[stress]
type = ADComputeLinearElasticStress
[]
[hydrostatic]
type = ADRankTwoInvariant
property_name = hydrostatic_stress
rank_two_tensor = stress
invariant = Hydrostatic
block = solid
outputs = exodus
[]
[vonmises]
type = ADRankTwoInvariant
property_name = vonmises_stress
rank_two_tensor = stress
invariant = VonMisesStress
block = solid
outputs = exodus
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
end_time = ${end_time}
dtmax = ${dtmax}
dtmin = 0.01
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
growth_factor = 1.2
optimal_iterations = 7
iteration_window = 2
linear_iteration_ratio = 100000
[]
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
nl_abs_tol = 1e-8
nl_rel_tol = 1e-6
nl_max_its = 12
[]
[Outputs]
file_base = 'out/print'
exodus = true
[]
(examples/PCM-GF/charging_idling.i)
# units are in meter kelvin second (m,kg,s)
tramp = 10
tcharge = 8000 # 3hr*3600
tidle = 86400 #24hr*3600
end_time = '${fparse tcharge+tidle}'
dtmax = 600
dt = 1
T_melting = '${fparse 718+273.15}' # Temperature at which the melting begins, Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
delta_T_pc = 8 # The temperature range of the melting/solidification process
L = 373.9e3 # Latent heat, from Singh et al. (2015)
sigma_foam_PCM = 5 # (from Wen's measurement of Gfoam+PCM in radial direction) (from Cfoam 70% dense foam = 28571.43) S/m (1/electrical resistivity (0.000035 ohm-m))
kappa_foam_PCM = 10 #18.8 # (average of Kxy = 14 W/m-K, Kz = 23.6 W/mK at T=700C) #from Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
rho_foam_PCM = 2050 # kg/m^3 #from Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
cp_foam_PCM = 1074 # J/kg-K #from Singh et al. Solar energy 159(2018) 270-282 (Prototype 1)
sigma_htf_pipe = 750750.75 # S/m (resistivity 1.332e-6 ohm-m at T = 700C) #Special metal data sheet
kappa_htf_pipe = 23.9 # W/m-K (at 700C) #Special metal datasheet
rho_htf_pipe = 8359.33 #kg/m^3
cp_htf_pipe = 419 # J/kg-K
sigma_insul_ht = 1e-9
kappa_insul_ht = 0.12 # W/m-K(Durablanket S from UNIFRAX) Wen emailed on 2023-03-31
rho_insul_ht = 2730 # kg/m^3(Durablanket S from UNIFRAX) Wen emailed on 2023-03-31
cp_insul_ht = 1130 # J/kg-K(Durablanket S from UNIFRAX) Wen emailed on 2023-03-31
sigma_air = 1e-12
kappa_air = 0.03 #file:///C:/Users/barua/Downloads/PDS-FOAMGLAS%20ONE-US-en.pdf
rho_air = 1.29 #file:///C:/Users/barua/Downloads/PDS-FOAMGLAS%20ONE-US-en.pdf
cp_air = 1000 #file:///C:/Users/barua/Downloads/PDS-FOAMGLAS%20ONE-US-en.pdf
htc_insul = 5
T_inf_insul = 300
htc_pipe = 0.01
T_inf_pipe = 300
T0 = 300
# i = 1 # This is the maximum current in constant-current charging
V = 21 # This is the maximum voltage in constant-voltage charging
[GlobalParams]
energy_densities = 'E H'
[]
[Mesh]
[fmg]
type = FileMeshGenerator
file = 'gold/geo.e'
[]
coord_type = RZ
uniform_refine = 1
[]
[Variables]
[Phi]
[]
[T]
initial_condition = ${T0}
[]
[]
[AuxVariables]
[T_old]
[AuxKernel]
type = ParsedAux
expression = 'T'
coupled_variables = 'T'
execute_on = 'INITIAL TIMESTEP_BEGIN'
[]
[]
[]
[Kernels]
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cp
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[energy_balance_3]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[Functions]
# [ramp_current]
# type = PiecewiseLinear
# x = '0 ${t0}'
# y = '0 ${i}'
# []
[ramp_voltage]
type = PiecewiseLinear
x = '0 ${tramp} ${tcharge} ${fparse tcharge+tramp} ${fparse tcharge+tidle}'
y = '0 ${V} ${V} 0 0'
[]
[]
[BCs]
[ground]
type = DirichletBC
variable = Phi
boundary = 'foam_id'
value = 0
[]
# [current]
# type = FunctionNeumannBC
# variable = Phi
# boundary = 'foam_od'
# function = ramp_current
# []
[CV]
type = FunctionDirichletBC
variable = Phi
boundary = 'foam_od'
function = ramp_voltage
[]
[hconv_insul]
type = ADMatNeumannBC
variable = T
boundary = 'insul_surf'
value = -1
boundary_material = qconv_insul
[]
[hconv_pipe]
type = ADMatNeumannBC
variable = T
boundary = 'pipe_id'
value = -1
boundary_material = qconv_pipe
[]
[]
[Materials]
[electrical_conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = sigma
subdomain_to_prop_value = 'foam_PCM ${sigma_foam_PCM} htf_pipe ${sigma_htf_pipe} insul_ht ${sigma_insul_ht} air ${sigma_air}'
[]
[charge_trasport]
type = BulkChargeTransport
electrical_energy_density = E
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[current]
type = CurrentDensity
current_density = i
electric_potential = Phi
[]
[thermal_conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = kappa
subdomain_to_prop_value = 'foam_PCM ${kappa_foam_PCM} htf_pipe ${kappa_htf_pipe} insul_ht ${kappa_insul_ht} air ${kappa_air}'
[]
[density]
type = ADPiecewiseConstantByBlockMaterial
prop_name = rho
subdomain_to_prop_value = 'foam_PCM ${rho_foam_PCM} htf_pipe ${rho_htf_pipe} insul_ht ${rho_insul_ht} air ${rho_air}'
[]
[gaussian_function]
type = ADParsedMaterial
property_name = D
expression = 'exp(-T*(T-Tm)^2/dT^2)/sqrt(3.1415926*dT^2)'
coupled_variables = 'T'
constant_names = 'Tm dT'
constant_expressions = '${T_melting} ${delta_T_pc}'
[]
[specific_heat_foam_PCM]
type = ADParsedMaterial
property_name = cp
expression = '${cp_foam_PCM} + ${L} * D'
material_property_names = 'D'
block = foam_PCM
outputs = exodus
[]
[specific_heat]
type = ADPiecewiseConstantByBlockMaterial
prop_name = cp
subdomain_to_prop_value = 'htf_pipe ${cp_htf_pipe} insul_ht ${cp_insul_ht} air ${cp_air}'
block = 'htf_pipe insul_ht air'
outputs = exodus
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = H
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
[]
[qconv_insul]
type = ADParsedMaterial
property_name = qconv_insul
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '${htc_insul} ${T_inf_insul}'
boundary = 'insul_surf'
[]
[qconv_pipe]
type = ADParsedMaterial
property_name = qconv_pipe
expression = 'htc*(T-T_inf)'
coupled_variables = 'T'
constant_names = 'htc T_inf'
constant_expressions = '${htc_pipe} ${T_inf_pipe}'
boundary = 'pipe_id'
[]
[delta_enthalpy]
type = ADParsedMaterial
property_name = delta_enthalpy
expression = 'rho*cp*(T-T_old)/2'
material_property_names = 'rho cp'
coupled_variables = 'T T_old'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
end_time = ${end_time}
dtmax = ${dtmax}
dtmin = 0.01
[TimeStepper]
type = IterationAdaptiveDT
dt = ${dt}
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
growth_factor = 1.2
optimal_iterations = 7
iteration_window = 2
linear_iteration_ratio = 100000
[]
[Predictor]
type = SimplePredictor
scale = 1
skip_after_failed_timestep = true
[]
nl_abs_tol = 1e-8
nl_rel_tol = 1e-6
nl_max_its = 12
# steady_state_detection = true
[]
[Postprocessors]
[volume_PCM]
type = VolumePostprocessor
block = 'foam_PCM'
execute_on = 'INITIAL TIMESTEP_END'
# outputs = none
[]
[delta_energy_absorbed_by_PCM]
type = ADElementIntegralMaterialProperty
mat_prop = delta_enthalpy
execute_on = 'INITIAL TIMESTEP_END'
block = 'foam_PCM'
outputs = none
[]
[delta_energy_absorbed_by_insul_ht]
type = ADElementIntegralMaterialProperty
mat_prop = delta_enthalpy
execute_on = 'INITIAL TIMESTEP_END'
block = 'insul_ht'
outputs = none
[]
[delta_energy_absorbed_by_air]
type = ADElementIntegralMaterialProperty
mat_prop = delta_enthalpy
execute_on = 'INITIAL TIMESTEP_END'
block = 'air'
outputs = none
[]
[delta_energy_absorbed_by_pipe]
type = ADElementIntegralMaterialProperty
mat_prop = delta_enthalpy
execute_on = 'INITIAL TIMESTEP_END'
block = 'htf_pipe'
outputs = none
[]
[energy_absorbed_by_PCM]
type = CumulativeValuePostprocessor
postprocessor = delta_energy_absorbed_by_PCM
execute_on = 'INITIAL TIMESTEP_END'
[]
[energy_absorbed_by_insul_ht]
type = CumulativeValuePostprocessor
postprocessor = delta_energy_absorbed_by_insul_ht
execute_on = 'INITIAL TIMESTEP_END'
[]
[energy_absorbed_by_air]
type = CumulativeValuePostprocessor
postprocessor = delta_energy_absorbed_by_air
execute_on = 'INITIAL TIMESTEP_END'
[]
[energy_absorbed_by_pipe]
type = CumulativeValuePostprocessor
postprocessor = delta_energy_absorbed_by_pipe
execute_on = 'INITIAL TIMESTEP_END'
[]
[total_energy]
type = ParsedPostprocessor
pp_names = 'energy_absorbed_by_PCM energy_absorbed_by_insul_ht energy_absorbed_by_air energy_absorbed_by_pipe'
function = 'energy_absorbed_by_PCM+energy_absorbed_by_insul_ht+energy_absorbed_by_air+energy_absorbed_by_pipe'
execute_on = 'INITIAL TIMESTEP_END'
[]
[max_energy]
type = TimeExtremeValue
postprocessor = total_energy
execute_on = 'INITIAL TIMESTEP_END'
outputs = none
[]
[energy_percentage]
type = ParsedPostprocessor
pp_names = 'total_energy max_energy'
function = 'if(total_energy < max_energy, total_energy/max_energy*100, 100)'
execute_on = 'INITIAL TIMESTEP_END'
[]
[power_input]
type = ADElementIntegralMaterialProperty
mat_prop = E
execute_on = 'INITIAL TIMESTEP_END'
# outputs = none
[]
[voltage]
type = FunctionValuePostprocessor
function = ramp_voltage
execute_on = 'INITIAL TIMESTEP_END'
# outputs = none
[]
[current_input]
type = ParsedPostprocessor
pp_names = 'voltage power_input'
function = 'power_input / voltage'
execute_on = 'TIMESTEP_END'
[]
[dt]
type = TimestepSize
execute_on = 'INITIAL TIMESTEP_END'
outputs = none
[]
[delta_energy_input]
type = ParsedPostprocessor
pp_names = 'dt power_input'
function = 'power_input * dt'
execute_on = 'INITIAL TIMESTEP_END'
outputs = none
[]
[energy_input]
type = CumulativeValuePostprocessor
postprocessor = delta_energy_input
execute_on = 'INITIAL TIMESTEP_END'
[]
[PCM_max_temperature]
type = NodalExtremeValue
variable = T
value_type = max
block = 'foam_PCM'
execute_on = 'INITIAL TIMESTEP_END'
[]
[PCM_min_temperature]
type = NodalExtremeValue
variable = T
value_type = min
block = 'foam_PCM'
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Outputs]
exodus = true
[]
(examples/SSB/CC_discharging.i)
I = 0.005 #mA
width = 0.05 #mm
in = '${fparse -I/width}'
t0 = '${fparse -1e-2/in}'
sigma_a = 0.2 #mS/mm
sigma_e = 0.1 #mS/mm
sigma_cp = 0.05 #mS/mm
sigma_ca = 0.2 #mS/mm
sigma_cm = 0.05 #mS/mm
Phi_penalty = 10
cmin_a = 1e-4 #mmol/mm^3
cmax_a = 1e-3 #mmol/mm^3
c_e = 5e-4 #mmol/mm^3
cmax_c = 1e-3 #mmol/mm^3
c_ref_entropy = 5e-5
D_cp = 5e-5 #mm^2/s
D_cm = 1e-4 #mm^2/s
D_a = 5e-4 #mm^2/s
D_e = 1e-4 #mm^2/s
c_penalty = 1
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
E_cp = 6e4
E_cm = 5e4
E_e = 5e4
E_a = 1e5
nu_cp = 0.3
nu_cm = 0.25
nu_e = 0.25
nu_a = 0.3
u_penalty = 1e8
Omega = 140
beta = 1e-4
CTE = 1e-5
rho = 2.5e-9 #Mg/mm^3
cv = 2.7e8 #mJ/Mg/K
kappa = 2e-4 #mJ/mm/K/s
htc = 9.5e-3
T_penalty = 1
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta'
deformation_gradient = F
mechanical_deformation_gradient = Fm
eigen_deformation_gradient = Fg
swelling_deformation_gradient = Fs
thermal_deformation_gradient = Ft
displacements = 'disp_x disp_y'
[]
[Problem]
restart_file_base = 'CV_charging_out_cp/LATEST'
[]
[Mesh]
[battery]
type = FileMeshGenerator
file = 'gold/ssb.msh'
[]
[interfaces]
type = BreakMeshByBlockGenerator
input = battery
add_interface_on_two_sides = true
split_interface = true
[]
use_displaced_mesh = false
[]
[Variables]
[Phi_ca]
block = cm
[]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[T]
[]
[]
[AuxVariables]
[c_ref]
[]
[T_ref]
initial_condition = ${T0}
[]
[stress]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = pk1
scalar_type = Hydrostatic
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Js]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Fs
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Jt]
order = CONSTANT
family = MONOMIAL
[AuxKernel]
type = ADRankTwoScalarAux
rank_two_tensor = Ft
scalar_type = ThirdInvariant
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[Phi0]
[]
[]
[ICs]
[c_ref_a]
type = ConstantIC
variable = c_ref
value = ${cmin_a}
block = 'a'
[]
[c_ref_e]
type = ConstantIC
variable = c_ref
value = ${c_e}
block = 'cm e'
[]
[c_ref_c]
type = ConstantIC
variable = c_ref
value = ${cmax_c}
block = 'cp'
[]
[]
[Kernels]
# Charge balance
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[charge_balance_ca]
type = RankOneDivergence
variable = Phi_ca
vector = i_ca
block = cm
[]
# Mass balance
[mass_balance_1]
type = TimeDerivative
variable = c
[]
[mass_balance_2]
type = RankOneDivergence
variable = c
vector = j
[]
# Momentum balance
[momentum_balance_x]
type = RankTwoDivergence
variable = disp_x
component = 0
tensor = pk1
factor = -1
[]
[momentum_balance_y]
type = RankTwoDivergence
variable = disp_y
component = 1
tensor = pk1
factor = -1
[]
# Energy balance
[energy_balance_1]
type = EnergyBalanceTimeDerivative
variable = T
density = rho
specific_heat = cv
[]
[energy_balance_2]
type = RankOneDivergence
variable = T
vector = h
[]
[heat_source]
type = MaterialSource
variable = T
prop = r
coefficient = -1
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'a_e cm_cp'
[]
[negative_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = -1
boundary = 'e_a cp_cm'
[]
[positive_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = 1
boundary = 'a_e cm_cp'
[]
[heat]
type = MaterialInterfaceNeumannBC
variable = T
neighbor_var = T
prop = he
factor = 1
boundary = 'a_e cm_cp e_a cp_cm'
[]
[continuity_c]
type = InterfaceContinuity
variable = c
neighbor_var = c
penalty = ${c_penalty}
boundary = 'cm_e'
[]
[continuity_Phi_ca]
type = InterfaceContinuity
variable = Phi_ca
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_cp'
[]
[continuity_Phi]
type = InterfaceContinuity
variable = Phi
neighbor_var = Phi
penalty = ${Phi_penalty}
boundary = 'cm_e'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[continuity_T]
type = InterfaceContinuity
variable = T
neighbor_var = T
penalty = ${T_penalty}
boundary = 'cp_cm cm_e e_a'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${fparse -in}'
[]
[]
[BCs]
[current]
type = FunctionNeumannBC
variable = Phi
boundary = right
function = in
[]
[potential]
type = DirichletBC
variable = Phi_ca
boundary = left
value = 0
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[hconv]
type = ADMatNeumannBC
variable = T
boundary = 'left right'
value = -1
boundary_material = qconv
[]
[]
[Constraints]
[ev_y]
type = EqualValueBoundaryConstraint
variable = disp_y
penalty = ${u_penalty}
secondary = top
[]
[]
[Materials]
# Electrodynamics
[conductivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma'
subdomain_to_prop_value = 'a ${sigma_a} e ${sigma_e} cm ${sigma_cm} cp ${sigma_cp}'
[]
[conductivity_ca]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'sigma_ca'
subdomain_to_prop_value = 'cm ${sigma_ca}'
block = cm
[]
[charge_transport]
type = BulkChargeTransport
electrical_energy_density = q
electric_potential = Phi
electric_conductivity = sigma
temperature = T
[]
[charge_transport_ca]
type = BulkChargeTransport
electrical_energy_density = q_ca
electric_potential = Phi_ca
electric_conductivity = sigma_ca
temperature = T
block = cm
[]
[current_density]
type = CurrentDensity
current_density = i
electric_potential = Phi
output_properties = i
outputs = exodus
[]
[current_density_ca]
type = CurrentDensity
current_density = i_ca
electric_potential = Phi_ca
output_properties = i_ca
outputs = exodus
block = cm
[]
# Chemical reactions
[diffusivity]
type = ADPiecewiseConstantByBlockMaterial
prop_name = 'D'
subdomain_to_prop_value = 'a ${D_a} e ${D_e} cm ${D_cm} cp ${D_cp}'
[]
[mobility]
type = ADParsedMaterial
f_name = M
args = 'c_ref T_ref'
material_property_names = 'D'
function = 'D*c_ref/${R}/T_ref'
[]
[chemical_energy]
type = EntropicChemicalEnergyDensity
chemical_energy_density = psi_c
concentration = c
ideal_gas_constant = ${R}
temperature = T_ref
reference_concentration = ${c_ref_entropy}
[]
[chemical_potential]
type = ChemicalPotential
chemical_potential = mu
concentration = c
energy_densities = 'dot(psi_m) dot(psi_c) chi q q_ca zeta m'
[]
[diffusion]
type = MassDiffusion
dual_chemical_energy_density = zeta
chemical_potential = mu
mobility = M
[]
[mass_flux]
type = MassFlux
mass_flux = j
chemical_potential = mu
outputs = exodus
[]
# Redox
[ramp]
type = ADGenericFunctionMaterial
prop_names = 'ramp'
prop_values = 'if(t<${t0},t/${t0},1)'
[]
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_a}; 2.77e-4*x^2-0.0069*x+0.0785'
# function = 'x:=c/${cmax_a}; -(122.12*x^6-321.81*x^5+315.59*x^4-141.26*x^3+28.218*x^2-1.9057*x+0.0785)*ramp'
args = c
material_property_names = 'ramp'
block = 'a'
[]
[OCP_cathode_NMC111]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax_c}; (6.0826-6.9922*x+7.1062*x^2-5.4549e-5*exp(124.23*x-114.2593)-2.5947*x^3)*ramp'
args = c
material_property_names = 'ramp'
block = 'cp'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'a_e'
[]
[charge_transfer_elyte_anode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_a}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'e_a'
[]
[charge_transfer_cathode_elyte]
type = ChargeTransferReaction
electrode = true
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cp_cm'
[]
[charge_transfer_elyte_cathode]
type = ChargeTransferReaction
electrode = false
charge_transfer_current_density = ie
charge_transfer_mass_flux = je
charge_transfer_heat_flux = he
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = 0.5
exchange_current_density = ${i0_c}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'cm_cp'
[]
# Thermal
[thermal_properties]
type = ADGenericConstantMaterial
prop_names = 'rho cv kappa'
prop_values = '${rho} ${cv} ${kappa}'
[]
[heat_conduction]
type = FourierPotential
thermal_energy_density = chi
thermal_conductivity = kappa
temperature = T
[]
[heat_flux]
type = HeatFlux
heat_flux = h
temperature = T
output_properties = h
outputs = exodus
[]
[heat_source]
type = VariationalHeatSource
heat_source = r
temperature = T
output_properties = r
outputs = exodus
[]
[conv]
type = ADParsedMaterial
f_name = qconv
function = '${htc}*(T-T_ref)'
args = 'T T_ref'
boundary = 'left right'
[]
# Mechanical
[stiffness_cp]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cp*nu_cp/(1+nu_cp)/(1-2*nu_cp)} ${fparse E_cp/2/(1+nu_cp)}'
block = cp
[]
[stiffness_cm]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_cm*nu_cm/(1+nu_cm)/(1-2*nu_cm)} ${fparse E_cm/2/(1+nu_cm)}'
block = cm
[]
[stiffness_e]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_e*nu_e/(1+nu_e)/(1-2*nu_e)} ${fparse E_e/2/(1+nu_e)}'
block = e
[]
[stiffness_a]
type = ADGenericConstantMaterial
prop_names = 'lambda G'
prop_values = '${fparse E_a*nu_a/(1+nu_a)/(1-2*nu_a)} ${fparse E_a/2/(1+nu_a)}'
block = a
[]
[swelling_coefficient]
type = ADGenericConstantMaterial
prop_names = 'beta'
prop_values = '${beta}'
[]
[swelling]
type = SwellingDeformationGradient
concentration = c
reference_concentration = c_ref
molar_volume = ${Omega}
swelling_coefficient = beta
[]
[thermal_expansion]
type = ThermalDeformationGradient
temperature = T
reference_temperature = T_ref
CTE = ${CTE}
[]
[defgrad]
type = MechanicalDeformationGradient
[]
[neohookean]
type = NeoHookeanSolid
elastic_energy_density = psi_m
lambda = lambda
shear_modulus = G
concentration = c
temperature = T
non_swelling_pressure = p
output_properties = 'p'
outputs = exodus
[]
[pk1]
type = FirstPiolaKirchhoffStress
first_piola_kirchhoff_stress = pk1
deformation_gradient_rate = dot(F)
[]
[]
[Postprocessors]
[V_l]
type = SideAverageValue
variable = Phi_ca
boundary = left
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V_r]
type = SideAverageValue
variable = Phi
boundary = right
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[V]
type = ParsedPostprocessor
function = 'V_l - V_r'
pp_names = 'V_l V_r'
execute_on = 'INITIAL TIMESTEP_END'
[]
[I]
type = ADSideIntegralMaterialProperty
property = i
component = 0
boundary = right
execute_on = 'INITIAL TIMESTEP_END'
[]
[dt]
type = TimestepSize
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[dC]
type = ParsedPostprocessor
function = 'dt*I'
pp_names = 'dt I'
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[C]
type = CumulativeValuePostprocessor
postprocessor = dC
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
[UserObjects]
[kill_V]
type = Terminator
expression = 'V <= 2.5'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-ksp_converged_reason'
# petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_strong_threshold -pc_hypre_boomeramg_interp_type -pc_hypre_boomeramg_coarsen_type -pc_hypre_boomeramg_agg_nl -pc_hypre_boomeramg_agg_num_paths -pc_hypre_boomeramg_truncfactor'
# petsc_options_value = 'hypre boomeramg 301 0.25 ext+i PMIS 4 2 0.4'
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
ignore_variables_for_autoscaling = 'T'
verbose = true
line_search = none
l_max_its = 300
l_tol = 1e-6
nl_rel_tol = 1e-6
nl_abs_tol = 1e-9
nl_max_its = 12
[Predictor]
type = SimplePredictor
scale = 1
[]
[TimeStepper]
type = IterationAdaptiveDT
dt = '${fparse t0/50}'
optimal_iterations = 6
iteration_window = 1
growth_factor = 1.2
cutback_factor = 0.2
cutback_factor_at_failure = 0.1
linear_iteration_ratio = 100
[]
end_time = 10000
[]
[Outputs]
exodus = true
csv = true
print_linear_residuals = false
[]