Simulate Battery Discharge-Rest-CC-CV Charge Cycle
R2026bThis example shows how to set up a pseudo-2D (P2D) battery model and simulate a full discharge-rest-charge cycle with a constant-current constant-voltage (CC-CV) charging protocol. The model uses parameters for a 12 Ah nickel cobalt aluminum (NCA) pouch lithium-ion cell from Jin et al. [1].
Define Active Materials
Create the anode and cathode active materials by specifying particle geometry, concentration limits, diffusion coefficients, reaction kinetics, and open-circuit potential (OCP) functions.
The OCP functions are polynomial fits from the parameterization study. The stoichiometric limits define the range of lithium intercalation used during cycling: the first value is the stoichiometry at 0% SOC and the second is the stoichiometry at 100% SOC.
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 Electrodes and Separator
Create the anode and cathode electrodes with thickness, porosity, Bruggeman coefficient, and electrical conductivity. The Bruggeman coefficient accounts for the tortuous path ions take through the porous electrode.
anodeElectrode = batteryElectrode( ... Thickness=5e-5, ... Porosity=0.21, ... BruggemanCoefficient=1.5, ... ElectricalConductivity=1000, ... ActiveMaterial=anodeMaterial); 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 a constant diffusion coefficient, ionic conductivity, and transference number. In this parameterization, these properties are not concentration-dependent.
electrolyte = batteryElectrolyte( ... DiffusionCoefficient=1.66e-11, ... IonicConductivity=0.29, ... TransferenceNumber=0.35);
Set Initial Conditions and Assemble the Model
Define the initial electrolyte concentration, state of charge (SOC), and temperature. Assemble all components into a P2D battery model.
ic = batteryInitialConditions( ... ElectrolyteConcentration=1200, ... StateOfCharge=1, ... Temperature=298.15); model = batteryP2DModel( ... Anode=anodeElectrode, ... Separator=separator, ... Cathode=cathodeElectrode, ... Electrolyte=electrolyte, ... InitialConditions=ic);
Define Cycling Steps
Set up a four-step cycling protocol:
Discharge: 1C constant current until voltage reaches 2.5 V
Rest: No current for 600 seconds
CC Charge: 1C constant current until voltage reaches 4.2 V
CV Charge: Hold voltage at 4.2 V until current drops to C/20
discharge = batteryCyclingStep; discharge.NormalizedCurrent = -1; discharge.CutoffVoltageLower = 2.5; discharge.OutputTimeStep = 10; rest = batteryCyclingStep; rest.NormalizedCurrent = 0; rest.CutoffTime = 600; rest.OutputTimeStep = 10; ccCharge = batteryCyclingStep; ccCharge.NormalizedCurrent = 1; ccCharge.CutoffVoltageUpper = 4.2; ccCharge.OutputTimeStep = 10; cvCharge = batteryCyclingStep; cvCharge.HoldVoltage = 4.2; cvCharge.CutoffNormalizedCurrent = 0.05; cvCharge.OutputTimeStep = 10; model.CyclingStep = [discharge, rest, ccCharge, cvCharge];
Solve the Model
Run the simulation for the complete cycling sequence.
R = solve(model);
Performing cycling step: 1 of 4 Performing cycling step: 2 of 4 Performing cycling step: 3 of 4 Performing cycling step: 4 of 4
Plot Results
Plot a summary of all quantities of interest from the simulation results.
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, 2380, 4766.5925, 7140.6833.](../../examples/pde/SimulateBatteryDischargeRestCCCVChargeCycleExample_01.png)
Visualize the terminal voltage and normalized current (C-rate) over time.
figure tiledlayout(2,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
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].
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
