Cyclic Ageing Simulation

In this example we run a cyclic ageing study for the Molicel P45B cell, with periodic reference performance tests (RPTs).

Note

This is a showcase example for demonstration only. We are happy to discuss how we develop and supply degradation models for your specific cell and use case. Contact us at hello@breathebatteries.com.

Objective

Run cyclic ageing experiments under defined conditions and track degradation state evolution using the Breathe PBM-Deg Simscape model.

The workflow includes:

SEI growth contribution
Lithium plating and stripping contribution
Mechanical damage (particle fracture) contribution
RPT-based tracking of capacity, resistance, and degradation states

Workflow

1. Set simulation settings

%Clear the workspace
bdclose all;
close all;
clear all;
clc;

numRPTs = 6;
numCyclesBtwnRPTs = 10;
crateCharge = 3.4;
crateDischarge = 1;
socLowerLim = 0;
socUpperLim = 1;

2. Setup the experiment

%% Load model parameters

% Get user profile directory
userProfile = getenv('USERPROFILE');

% Build the path relative to user profile
filePath = fullfile(userProfile, ...
    'Breathe Battery Technologies', ...
    'Breathe Model Projects - Documents', ...
    'Proj_BM_Mono', ...
    '07_Breathe_Model_Data', ...
    '03_PBMDegv1', ...
    '02_Parameters', ...
    'MoliP45_PBMDegv1_GD_Format');

% Load optimal parameter set
load(filePath);

% Create the ageing state object
ageState = Degradation_Experiments.CurrentAgeingState("myCyclicAgeingExp");
ageState.Update_State('porosityAnode', ones(6,1) * BrPBMDegv1.Components.porosityAnode);

% Create Cyclic Ageing Experiment object
cycObj = Degradation_Experiments.CyclicAgeingStep("CyclicAgeing");

% Create RPT object
rptObj = Degradation_Experiments.RPT("myRPT");

3. Set experimental conditions

% Set cyclic ageing experiment settings
cycObj.List_Exp_Settings;

normCap = BrPBMDegv1.Components.capCell_BOL;
cycObj.Update_Exp_Settings('numberOfCycles', numCyclesBtwnRPTs, ...
    'currCharge', crateCharge * normCap, ...
    'currDischarge', -crateDischarge * normCap, ...
    'socMinCyc', socLowerLim, ...
    'socMaxCyc', socUpperLim);

% List and update RPT settings
rptObj.List_Exp_Settings();
rptObj.Update_Exp_Settings('currTrkdCapMes', 4.5, 'currCharge', 4.5);

4. Build the model and run simulation

% Build the experiment with proper signal specification
potentials = ["voltModel", "voltAnode", "voltCathode", "tempSurfaceModel"];
cycler     = ["stepNumModel"];
states     = ["socModel"];
fluxes     = ["currModel"];

rptObj.Build_Experiment( ...
    "BrPBMDegv1", ...
    'outPotBusSignals', potentials, ...
    'outBusCyclerSignals', cycler, ...
    'outStateBusSignals', states, ...
    'outFluxBusSignals', fluxes ...
    );

cycObj.Build_Experiment( ...
    "BrPBMDegv1", ...
    'outPotBusSignals', ["voltModel"], ...
    'outFluxBusSignals', ["currModel"], ...
    'outStateBusSignals', ["lamneModel", "lliModel", "lliSeiModel", "lliPlatingModel"]);

% Run the experiment
rptObj.Run_Experiment(ageState);
for k = 1:numRPTs-1
    cycObj.Run_Experiment(ageState);
    rptObj.Run_Experiment(ageState);
end

5. Plot RPT summary results

cycleNumber = rptObj.results.cycleNumber;
capacity = rptObj.results.trackedCapacity_Ah;
pulseRes = rptObj.results.pulseResistance_Ohm * 1000;
lli = rptObj.results.lli * 100;
lamne = rptObj.results.lamne * 100;
lampe = rptObj.results.lampe * 100;

figure;
tiledlayout(2,3, 'Padding', 'compact', 'TileSpacing', 'compact');

% 1. Capacity
nexttile;
plot(cycleNumber, capacity / capacity(1) * 100, '-o', 'LineWidth', 1.5, 'MarkerSize', 6);
xlabel('Cycle Number'); ylabel('Capacity (Ah)');
title('Capacity vs Cycle Number');
grid on;

% 2. Pulse Resistance
nexttile;
plot(cycleNumber, pulseRes, '-o', 'LineWidth', 1.5, 'MarkerSize', 6);
xlabel('Cycle Number'); ylabel('Resistance (m\Omega)');
title('Pulse Resistance vs Cycle Number');
grid on;

% 3. LLI
nexttile;
plot(cycleNumber, lli, '-o', 'LineWidth', 1.5, 'MarkerSize', 6);
xlabel('Cycle Number'); ylabel('LLI (%)');
title('Loss of Lithium Inventory');
grid on;

% 4. LAM_NE
nexttile;
plot(cycleNumber, lamne, '-o', 'LineWidth', 1.5, 'MarkerSize', 6);
xlabel('Cycle Number'); ylabel('LAM-NE (%)');
title('Loss of Active Material - NE');
grid on;

% 5. LAM_PE
nexttile;
plot(cycleNumber, lampe, '-o', 'LineWidth', 1.5, 'MarkerSize', 6);
xlabel('Cycle Number'); ylabel('LAM-PE (%)');
title('Loss of Active Material - PE');
grid on;

Example output

RPT summary output

6. Plot electrochemical signals during RPTs

desiredStepNumber = 4;
signalNames = {'voltAnode', 'voltModel', 'tempSurfaceModel'};

logs = rptObj.results.logsTable;
nCyc = numel(logs);
colors = [linspace(0,1,nCyc)', zeros(nCyc,1), zeros(nCyc,1)];

for s = 1:numel(signalNames)
    sig = signalNames{s};
    anyPlotted = false;
    figure('Color','w');
    hold on; grid on; box on;

    for i = 1:nCyc
        T = logs{i};
        if ~istable(T) || ~all(ismember({'Time','stepNumModel',sig}, T.Properties.VariableNames))
            continue;
        end
        try
            [tSeg, ySeg] = getStepSegment(T.Time, T.(sig), T.stepNumModel, desiredStepNumber);
        catch
            continue;
        end
        plot(tSeg, ySeg, 'LineWidth', 1.2, 'Color', colors(i,:), 'DisplayName', sprintf('RPT No: %d', i));
        anyPlotted = true;
    end

    if ~anyPlotted
        close(gcf);
        warning('No data for "%s" at step %d.', sig, desiredStepNumber);
        continue;
    end

    xlabel('Time (s)');
    ylabel(sig);
    title(sprintf('%s vs Time across Cycles', sig));
    if nCyc <= 12
        legend('Location', 'best');
    end
    colormap(colors);
    caxis([1 nCyc]);
    cb = colorbar;
    set(cb, 'Ticks', [1 nCyc], 'TickLabels', {'First','Last'});
    hold off;
end

Example outputs

voltAnode:

voltAnode output

voltModel:

voltModel output

tempSurfaceModel:

tempSurfaceModel output

7. Utility function

function [tSeg, ySeg] = getStepSegment(t, y, step, stepNum)
idx = (step == stepNum);
if ~any(idx), error('Step %d not found.', stepNum); end
t0 = t(find(idx,1,'first'));
tSeg = t(idx) - t0;
ySeg = y(idx);
end

Version 73723266