- charge_transfer_coefficientThe dimensionless charge transfer coefficient
C++ Type:double
Controllable:No
Description:The dimensionless charge transfer coefficient
- charge_transfer_current_densityGive the current density a name
C++ Type:MaterialPropertyName
Controllable:No
Description:Give the current density a name
- charge_transfer_heat_fluxGive the heat flux a name
C++ Type:MaterialPropertyName
Controllable:No
Description:Give the heat flux a name
- charge_transfer_mass_fluxGive the mass flux a name
C++ Type:MaterialPropertyName
Controllable:No
Description:Give the mass flux a name
- electric_potentialElectric potential
C++ Type:std::vector<VariableName>
Controllable:No
Description:Electric potential
- exchange_current_densityThe exchange current density (normal to the interface) for the electrode/electrolyte interface
C++ Type:double
Controllable:No
Description:The exchange current density (normal to the interface) for the electrode/electrolyte interface
- faraday_constantThe Faraday's constant
C++ Type:double
Controllable:No
Description:The Faraday's constant
- ideal_gas_constantThe ideal gas constant
C++ Type:double
Controllable:No
Description:The ideal gas constant
- neighbor_electric_potentialNeighbor electric potential
C++ Type:std::vector<VariableName>
Controllable:No
Description:Neighbor electric potential
- open_circuit_potentialThe open-circuit potential
C++ Type:MaterialPropertyName
Controllable:No
Description:The open-circuit potential
- temperatureThe temperature
C++ Type:std::vector<VariableName>
Controllable:No
Description:The temperature
ChargeTransferReaction
buildconstruction:Undocumented Class
The ChargeTransferReaction has not been documented. The content listed below should be used as a starting point for documenting the class, which includes the typical automatic documentation associated with a MooseObject; however, what is contained is ultimately determined by what is necessary to make the documentation clear for users.
The Butler-Volmer condition for current density across the electrode/electrolyte interface.
Overview
Example Input File Syntax
Input Parameters
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
- declare_suffixAn optional suffix parameter that can be appended to any declared 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 declared properties. The suffix will be prepended with a '_' character.
- degradationInterface degradation
C++ Type:MaterialPropertyName
Controllable:No
Description:Interface degradation
- interface_resistanceInterface resistance. The resistance is assumed to be independent of the transfer current, for now.
C++ Type:MaterialPropertyName
Controllable:No
Description:Interface resistance. The resistance is assumed to be independent of the transfer current, for now.
- 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
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
- outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
Outputs Parameters
Input Files
- (test/tests/chemical-electrical/diffusion_vs_migration.i)
- (examples/LiB/CC_discharging.i)
- (examples/SSB/charging.i)
- (examples/SSB/CC_charging.i)
- (examples/SSB_3D/charging.i)
- (examples/degradation/base.i)
- (examples/SSB_3D/CC_charging.i)
- (examples/SSB/CV_charging.i)
- (examples/SSB_3D/CC_discharging.i)
- (test/tests/chemical-electrical-thermal/thermal_effects.i)
- (examples/SSB_3D/CV_charging.i)
- (test/tests/chemical-electrical-thermal-mechanical/pressure.i)
- (examples/SSB/CC_discharging.i)
- (examples/LiB/CV_charging.i)
- (examples/LiB/CC_charging.i)
- (test/tests/electrical/bv.i)
(test/tests/chemical-electrical/diffusion_vs_migration.i)
m = 1 # larger m means more migration compared to diffusion
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 = '${fparse sigma_se/F/F/m}'
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 = 1e-1 #mA/mm^2
i0_c = 1e-1 #mA/mm^2
[GlobalParams]
energy_densities = 'dot(psi_c) q zeta 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]
[]
[]
[AuxVariables]
[c_ref]
[]
[T]
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
[]
[]
[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'
[]
[]
[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'
[]
[]
[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'
[]
[]
[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'
[]
[]
[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]
file_base = 'diffusion_vs_migration'
exodus = true
print_linear_residuals = false
[]
(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/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/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/degradation/base.i)
I = 5e-4 #mA
width = 0.03 #mm
in = '${fparse -I/width}'
t0 = '${fparse -1e-2/in}'
dt = '${fparse t0/50}'
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
mu0 = 1e3
h0 = 1e-9
hc = 1e-3
Omega_sei = 100
A_sei = 0.8
Q_sei = 5e3
rho0_sei = 5e2
R = 8.3145 #mJ/mmol/K
T0 = 300 #K
F = 96485 #mC/mmol
i0_a = 1e-4 #mA/mm^2
i0_c = 1e-4 #mA/mm^2
E_c = 1e5
E_e = 1e4
E_a = 2e5
nu_c = 0.3
nu_e = 0.25
nu_a = 0.3
Ei = 1e7
Gi = 5e6
Acr = 1e-5
Qcr = 4e3
ncr = 5
Tn0 = 10
Gc = 2.06e-5
u_penalty = 1e8
Omega = 60
beta = 1e-3
CTE = 1e-5
[GlobalParams]
energy_densities = 'dot(psi_m) dot(psi_c) q zeta'
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 = 1
[]
[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
[]
use_displaced_mesh = false
[]
[Variables]
[Phi]
[]
[c]
[]
[disp_x]
[]
[disp_y]
[]
[]
[AuxVariables]
[c_ref]
[]
[Phi0]
[]
[T]
initial_condition = ${T0}
[]
[T_ref]
initial_condition = ${T0}
[]
[]
[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
expression = 'Phi'
coupled_variables = '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
[]
[]
[Modules]
[TensorMechanics]
[CohesiveZoneMaster]
[interface]
boundary = 'anode_elyte'
strain = SMALL
use_automatic_differentiation = true
[]
[]
[]
[]
[InterfaceKernels]
[negative_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
factor = -1
boundary = 'cathode_elyte'
[]
[positive_current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ie
boundary = 'anode_elyte'
[]
[negative_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
factor = -1
boundary = 'cathode_elyte'
[]
[positive_mass]
type = MaterialInterfaceNeumannBC
variable = c
neighbor_var = c
prop = je
boundary = 'anode_elyte'
[]
[continuity_disp_x]
type = InterfaceContinuity
variable = disp_x
neighbor_var = disp_x
penalty = ${u_penalty}
boundary = 'cathode_elyte'
[]
[continuity_disp_y]
type = InterfaceContinuity
variable = disp_y
neighbor_var = disp_y
penalty = ${u_penalty}
boundary = 'cathode_elyte'
[]
[]
[BCs]
[pull_x]
type = NeumannBC
variable = disp_x
value = 1
boundary = 'right'
[]
[fix_x]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left'
[]
[fix_y]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom'
[]
[]
[Functions]
[in]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${in}'
[]
[in_discharging]
type = PiecewiseLinear
x = '0 ${t0}'
y = '0 ${fparse -in}'
[]
[]
[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
property_name = M
coupled_variables = 'c_ref T_ref'
material_property_names = 'D'
expression = '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
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
[]
# Redox
[OCP_anode_graphite]
type = ADParsedMaterial
f_name = U
function = 'x:=c/${cmax}; (2.77e-4*x^2-0.0069*x+0.0785)'
args = c
material_property_names = 'ramp'
boundary = 'anode_elyte'
[]
[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'
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
interface_resistance = rho_sei
degradation = g
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'
[]
[SEI]
type = SEIGrowth
thickness = h_sei
initial_thickness = ${h0}
charge_transfer_mass_flux = je
characteristic_thickness = ${hc}
molar_volume = ${Omega_sei}
correction = ${A_sei}
activation_energy = ${Q_sei}
ideal_gas_constant = ${R}
temperature = T
boundary = 'anode_elyte'
[]
[SEI_resistance]
type = ADParsedMaterial
property_name = rho_sei
expression = '${rho0_sei} * h_sei'
material_property_names = 'h_sei'
boundary = 'anode_elyte'
[]
# 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
[]
[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)
[]
[traction]
type = InterfaceTractionWithCreepDegradation
normal_stiffness = ${Ei}
tangential_stiffness = ${Gi}
normal_traction = Tn
activation_energy = ${Qcr}
ideal_gas_constant = ${R}
temperature = T
damage = Di
degradation = g
creep_displacement_jump = juc
fracture_driving_energy = psi_f
energy_release_rate = ${Gc}
reference_normal_traction = ${Tn0}
creep_coefficient = ${Acr}
creep_exponent = ${ncr}
residual_stiffness = 1e-4
boundary = 'anode_elyte'
[]
[]
[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
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_c_min]
type = NodalExtremeValue
variable = c
value_type = min
block = cathode
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_a_min]
type = NodalExtremeValue
variable = c
value_type = min
block = anode
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[c_c_max]
type = NodalExtremeValue
variable = c
value_type = max
block = cathode
outputs = none
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'
[]
[SEI_thickness]
type = ADSideAverageMaterialProperty
property = h_sei
boundary = 'anode_elyte'
execute_on = 'INITIAL TIMESTEP_END'
[]
[Di]
type = ADSideAverageMaterialProperty
property = Di
boundary = 'anode_elyte'
execute_on = 'INITIAL TIMESTEP_END'
[]
[psi_f]
type = ADSideAverageMaterialProperty
property = psi_f
boundary = 'anode_elyte'
execute_on = 'INITIAL TIMESTEP_END'
[]
[dV]
type = ChangeOverTimePostprocessor
postprocessor = V
outputs = none
execute_on = 'INITIAL TIMESTEP_END'
[]
[target_V_dt]
type = ParsedPostprocessor
function = 'if(abs(dV)>0.01,${dt},1e6)'
pp_names = 'dV'
outputs = none
execute_on = 'TIMESTEP_END'
[]
[]
(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
[]
(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
[]
(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
[]
(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
[]
(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/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
[]
(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
[]
(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
[]
(test/tests/electrical/bv.i)
n = 100
sigma = 10
alpha = 0.5
i0 = 1
F = 96485
R = 8.3145
T = 300
U = 0.1
[Mesh]
[left]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = 1
ymin = 0
ymax = 1
nx = ${n}
ny = 1
boundary_name_prefix = anode
[]
[right]
type = GeneratedMeshGenerator
dim = 2
xmin = 1
xmax = 2
ymin = 0
ymax = 1
nx = ${n}
ny = 1
boundary_name_prefix = elyte
boundary_id_offset = 100
[]
[stitch]
type = StitchedMeshGenerator
inputs = 'left right'
stitch_boundaries_pairs = 'anode_right elyte_left'
[]
[rename]
type = RenameBoundaryGenerator
input = stitch
old_boundary = 'anode_left elyte_right'
new_boundary = 'left right'
[]
[anode]
type = SubdomainBoundingBoxGenerator
input = rename
bottom_left = '0 0 0'
top_right = '1 1 0'
block_id = 0
block_name = anode
[]
[elyte]
type = SubdomainBoundingBoxGenerator
input = anode
bottom_left = '1 0 0'
top_right = '2 1 0'
block_id = 1
block_name = elyte
[]
[break]
type = BreakMeshByBlockGenerator
input = elyte
split_interface = true
add_interface_on_two_sides = true
[]
[]
[Variables]
[Phi]
[]
[]
[AuxVariables]
[T]
initial_condition = ${T}
[]
[]
[Kernels]
[charge_balance]
type = RankOneDivergence
variable = Phi
vector = i
[]
[]
[InterfaceKernels]
[current]
type = MaterialInterfaceNeumannBC
variable = Phi
neighbor_var = Phi
prop = ibv
boundary = 'anode_elyte'
[]
[]
[Functions]
[in]
type = ParsedFunction
expression = '-t'
[]
[]
[BCs]
[left]
type = DirichletBC
variable = Phi
boundary = left
value = 0
[]
[right]
type = FunctionNeumannBC
variable = Phi
boundary = right
function = in
[]
[]
[Materials]
[electric_constants]
type = ADGenericConstantMaterial
prop_names = 'sigma'
prop_values = '${sigma}'
[]
[charge_trasport]
type = BulkChargeTransport
electrical_energy_density = E
electric_potential = Phi
electric_conductivity = sigma
[]
[current]
type = CurrentDensity
current_density = i
energy_densities = 'E'
electric_potential = Phi
[]
[OCP_anode]
type = ADParsedMaterial
property_name = U
expression = '${U}'
boundary = 'anode_elyte'
[]
[charge_transfer_anode_elyte]
type = ChargeTransferReaction
charge_transfer_current_density = ibv
charge_transfer_mass_flux = jbv
charge_transfer_heat_flux = hbv
electric_potential = Phi
neighbor_electric_potential = Phi
charge_transfer_coefficient = ${alpha}
exchange_current_density = ${i0}
faraday_constant = ${F}
ideal_gas_constant = ${R}
temperature = T
open_circuit_potential = U
boundary = 'anode_elyte'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
automatic_scaling = true
dt = 1e-3
end_time = 1
timestep_tolerance = 1e-10
[]
[Postprocessors]
[Phi_r]
type = PointValue
variable = Phi
point = '2 0 0'
outputs = none
[]
[delta_Phi]
type = ParsedPostprocessor
function = '-Phi_r'
pp_names = 'Phi_r'
[]
[in]
type = FunctionValuePostprocessor
function = in
[]
[]
[Outputs]
exodus = true
[]