Simulate SEI Growth During Battery Cycling
R2026bThis example shows how to set up a pseudo-2D (P2D) battery model with a solid electrolyte interphase (SEI) aging model and simulate repeated cycling to observe SEI layer growth over time. The model uses parameters for a 12 Ah nickel cobalt aluminum (NCA) pouch lithium-ion cell from Jin et al. [1] with typical SEI parameters for a lithium-ion cell with an EC-based electrolyte.
Define Active Materials
Create the anode and cathode active materials. The open-circuit potential functions use polynomial fits from [1], evaluated with normalized stoichiometry.
anodeMaterial = batteryActiveMaterial( ... ParticleRadius=0.834e-6, ... MaximumSolidConcentration=16.1e3, ... VolumeFraction=0.58, ... DiffusionCoefficient=2e-14, ... ReactionRate=4.5229e-11, ... OpenCircuitPotential=@anodeOCP, ... StoichiometricLimits=[0.126 0.676]); cathodeMaterial = batteryActiveMaterial( ... ParticleRadius=0.834e-6, ... MaximumSolidConcentration=23.9e3, ... VolumeFraction=0.5, ... DiffusionCoefficient=8.38e-16, ... ReactionRate=2.0976e-11, ... OpenCircuitPotential=@cathodeOCP, ... StoichiometricLimits=[0.442 0.936]);
Define SEI Model
Create the SEI aging model. The SEI layer grows on the anode surface during cycling, consuming lithium and increasing impedance. The parameters represent typical values for an SEI layer composed primarily of lithium carbonate (LiCO) in an EC-based electrolyte.
seiModel = batterySEIModel( ... IonicResistivity=5e4, ... MolecularWeight=0.07369, ... Density=2110, ... KineticRateConstant=5e-13, ... SolventConcentration=4541, ... CathodicTransferCoeff=0.5, ... StoichiometricCoefficient=2, ... SolventDiffusivity=2.5e-22, ... EquilibriumPotential=0.4);
Define Electrodes and Separator
Create the electrodes and separator. The anode includes the SEI aging model.
anodeElectrode = batteryElectrode( ... Thickness=5e-5, ... Porosity=0.21, ... BruggemanCoefficient=1.5, ... ElectricalConductivity=1000, ... ActiveMaterial=anodeMaterial, ... Aging=seiModel); cathodeElectrode = batteryElectrode( ... Thickness=3.64e-5, ... Porosity=0.25, ... BruggemanCoefficient=1.5, ... ElectricalConductivity=0.002, ... ActiveMaterial=cathodeMaterial); separator = batterySeparator( ... Thickness=2.54e-5, ... Porosity=0.5, ... BruggemanCoefficient=1.5);
Define Electrolyte
Create the electrolyte with constant transport properties.
electrolyte = batteryElectrolyte( ... DiffusionCoefficient=1.66e-11, ... IonicConductivity=0.29, ... TransferenceNumber=0.35);
Set Initial Conditions and Assemble the Model
Define initial conditions including an initial SEI thickness of 5 nm (representing a formed cell). Assemble all components into a P2D battery model.
ic = batteryInitialConditions( ... ElectrolyteConcentration=1200, ... StateOfCharge=1, ... InitialSEIThickness=5e-9, ... Temperature=298.15); model = batteryP2DModel( ... Anode=anodeElectrode, ... Separator=separator, ... Cathode=cathodeElectrode, ... Electrolyte=electrolyte, ... InitialConditions=ic);
Define Cycling Protocol
Set up a cycling protocol designed to observe SEI growth:
Initial CV charge: Hold at 4.2 V until C/100 to ensure full charge
Rest: 1 hour rest period
Capacity check: 0.5C discharge to 2.5 V
Aging cycles (repeated 10 times):
- CC charge at 0.5C to 4.2 V
- CV hold at 4.2 V until C/20
- Discharge at 0.5C to 2.5 V
Final capacity check: 0.5C discharge after charging
% Initial CV hold to ensure full charge cvInit = batteryCyclingStep; cvInit.HoldVoltage = 4.2; cvInit.CutoffNormalizedCurrent = 0.01; cvInit.OutputTimeStep = 10; % Rest rest = batteryCyclingStep; rest.NormalizedCurrent = 0; rest.CutoffTime = 3600; rest.OutputTimeStep = 60; % Initial capacity check discharge (0.5C) capCheck = batteryCyclingStep; capCheck.NormalizedCurrent = -0.5; capCheck.CutoffVoltageLower = 2.5; capCheck.OutputTimeStep = 10; % Aging cycle: CC charge at 0.5C ccCharge = batteryCyclingStep; ccCharge.NormalizedCurrent = 0.5; ccCharge.CutoffVoltageUpper = 4.2; ccCharge.OutputTimeStep = 10; % Aging cycle: CV hold at 4.2 V cvHold = batteryCyclingStep; cvHold.HoldVoltage = 4.2; cvHold.CutoffNormalizedCurrent = 0.05; cvHold.OutputTimeStep = 10; % Aging cycle: Discharge at 0.5C discharge = batteryCyclingStep; discharge.NormalizedCurrent = -0.5; discharge.CutoffVoltageLower = 2.5; discharge.OutputTimeStep = 10; % Final charge before capacity check ccChargeFinal = batteryCyclingStep; ccChargeFinal.NormalizedCurrent = 0.5; ccChargeFinal.CutoffVoltageUpper = 4.2; ccChargeFinal.OutputTimeStep = 10; cvHoldFinal = batteryCyclingStep; cvHoldFinal.HoldVoltage = 4.2; cvHoldFinal.CutoffNormalizedCurrent = 0.01; cvHoldFinal.OutputTimeStep = 10; % Final capacity check discharge (0.5C) capCheckFinal = batteryCyclingStep; capCheckFinal.NormalizedCurrent = -0.5; capCheckFinal.CutoffVoltageLower = 2.5; capCheckFinal.OutputTimeStep = 10;
Assemble the complete cycling protocol with 10 aging cycles.
numCycles = 10; model.CyclingStep = [cvInit, rest, capCheck, ... repmat([ccCharge, cvHold, discharge], 1, numCycles), ... ccChargeFinal, cvHoldFinal, capCheckFinal];
Solve the Model
Run the simulation for the complete cycling protocol.
R = solve(model);
Performing cycling step: 1 of 36 Performing cycling step: 2 of 36 Performing cycling step: 3 of 36 Performing cycling step: 4 of 36 Performing cycling step: 5 of 36 Performing cycling step: 6 of 36 Performing cycling step: 7 of 36 Performing cycling step: 8 of 36 Performing cycling step: 9 of 36 Performing cycling step: 10 of 36 Performing cycling step: 11 of 36 Performing cycling step: 12 of 36 Performing cycling step: 13 of 36 Performing cycling step: 14 of 36 Performing cycling step: 15 of 36 Performing cycling step: 16 of 36 Performing cycling step: 17 of 36 Performing cycling step: 18 of 36 Performing cycling step: 19 of 36 Performing cycling step: 20 of 36 Performing cycling step: 21 of 36 Performing cycling step: 22 of 36 Performing cycling step: 23 of 36 Performing cycling step: 24 of 36 Performing cycling step: 25 of 36 Performing cycling step: 26 of 36 Performing cycling step: 27 of 36 Performing cycling step: 28 of 36 Performing cycling step: 29 of 36 Performing cycling step: 30 of 36 Performing cycling step: 31 of 36 Performing cycling step: 32 of 36 Performing cycling step: 33 of 36 Performing cycling step: 34 of 36 Performing cycling step: 35 of 36 Performing cycling step: 36 of 36
Plot Results
Plot a summary of all quantities of interest.
plotSummary(R)
![Figure contains 6 axes objects. Axes object 1 with title Terminal Voltage [V], xlabel Time [s] contains an object of type line. Axes object 2 with title Normalized Current, xlabel Time [s] contains an object of type line. Axes object 3 with title Liquid Concentration [mol/m Cubed baseline ], xlabel Thickness [m] contains 7 objects of type line, rectangle. Axes object 4 with title Normalized Average Solid Concentration, xlabel Thickness [m] contains 7 objects of type line, rectangle. Axes object 5 with title Liquid Potential [V], xlabel Thickness [m] contains 7 objects of type line, rectangle. Axes object 6 with title Solid Potential [V], xlabel Thickness [m] contains 7 objects of type line, rectangle. These objects represent 0, 51697.49489, 100379.1482, 149044.1849.](../../examples/pde/SimulateSEIGrowthDuringBatteryCyclingExample_01.png)
Visualize terminal voltage, current, and SEI thickness over time.
figure tiledlayout(3,1) nexttile plot(R.SolutionTimes/3600, R.TerminalVoltage, LineWidth=1.5) xlabel("Time (hours)") ylabel("Voltage (V)") title("Terminal Voltage") grid on nexttile plot(R.SolutionTimes/3600, R.NormalizedCurrent, LineWidth=1.5) xlabel("Time (hours)") ylabel("C-rate") title("Normalized Current") grid on nexttile plot(R.SolutionTimes/3600, R.SEIThickness*1e9, LineWidth=1.5) xlabel("Time (hours)") ylabel("SEI Thickness (nm)") title("SEI Layer Growth") grid on
References
[1] Jin N, Danilov DL, Van den Hof PMJ, Donkers MCF. Parameter estimation of an electrochemistry-based lithium-ion battery model using a two-step procedure and a parameter sensitivity analysis. Int J Energy Res. 2018;42:2417-2430. https://doi.org/10.1002/er.4022
Local Functions
Open-circuit potential functions for the anode and cathode, defined as polynomial fits from [1]. The functions convert absolute stoichiometry to normalized stoichiometry before evaluating the polynomials.
function Un = anodeOCP(theta) y = (theta - 0.126) ./ (0.676 - 0.126); gamma = [0.0004, -0.0145, 0.1115, -0.6830, 0.8020, ... -0.3611, 0.1115, 0.0171, 0.1115, 0.0171]; Un = gamma(1).*y.^(-1) ... + gamma(2).*y.^(-0.5) ... + gamma(3) ... + gamma(4).*y.^(0.5) ... + gamma(5).*y.^(1.0) ... + gamma(6).*y.^(1.5) ... + gamma(7).*exp(gamma(8).*y) ... + gamma(9).*exp(gamma(10).*y); end function Up = cathodeOCP(theta) y = (0.936 - theta) ./ (0.936 - 0.442); gamma = [-2.2049e3, 11.2250, -112.9966, 647.4250, -2.083e3, ... 3.7868e3, -3.6869e3, 1.5388e3, 2.2079e3, -0.0464]; Up = gamma(1) ... + gamma(2).*y.^1 ... + gamma(3).*y.^2 ... + gamma(4).*y.^3 ... + gamma(5).*y.^4 ... + gamma(6).*y.^5 ... + gamma(7).*y.^6 ... + gamma(8).*y.^7 ... + gamma(9).*exp(gamma(10).*y.^10); end
