Note

This notebook is already available in your BattMo installation. In Matlab, run

open runSiliconGraphiteBattery

Simulation of a composite active material

This example shows how to simulate a composite active material

Setup the properties of the battery

We load the property of a composite silicon graphite electrode.

jsonstruct_composite_material = parseBattmoJson('ParameterData/BatteryCellParameters/LithiumIonBatteryCell/composite_silicon_graphite.json');

In this structure, we have the material property of two active materials, ActiveMaterial1 and ActiveMaterial2.

flattenStruct(jsonstruct_composite_material);

$\begin{array} {|c|c|c|}\hline\text{ } & \text{parameter name} & \text{parameter value} \\ \hline\text{1} & \text{'NegativeElectrode.Coating.activeMaterialModelSetup.composite'} & \text{1} \\ \hline\text{2} & \text{'NegativeElectrode.Coating.ActiveMaterial1.specificHeatCapacity'} & \text{632} \\ \hline\text{3} & \text{'NegativeElectrode.Coating.ActiveMaterial1.heatCapacity'} & \text{632000} \\ \hline\text{4} & \text{'NegativeElectrode.Coating.ActiveMaterial1.thermalConductivity'} & \text{1.0400} \\ \hline\text{5} & \text{'NegativeElectrode.Coating.ActiveMaterial1.electronicConductivity'} & \text{100} \\ \hline\text{6} & \text{'NegativeElectrode.Coating.ActiveMaterial1.massFraction'} & \text{0.9200} \\ \hline\text{7} & \text{'NegativeElectrode.Coating.ActiveMaterial1.density'} & \text{2240} \\ \hline\text{8} & \text{'NegativeElectrode.Coating.ActiveMaterial1.Interface.saturationConcentration'} & \text{30555} \\ \hline\text{9} & \text{'NegativeElectrode.Coating.ActiveMaterial1.Interface.volumetricSurfaceArea'} & \text{723600} \\ \hline\text{10} & \text{'NegativeElectrode.Coating.ActiveMaterial1.Interface.numberOfElectronsTransferred'} & \text{1} \\ \hline\text{11} & \text{'NegativeElectrode.Coating.ActiveMaterial1.Interface.activationEnergyOfReaction'} & \text{5000} \\ \hline\text{12} & \text{'NegativeElectrode.Coating.ActiveMaterial1.Interface.reactionRateConstant'} & \text{5.0310e-11} \\ \hline\text{13} & \text{'NegativeElectrode.Coating.ActiveMaterial1.Interface.guestStoichiometry100'} & \text{0.8855} \\ \hline\text{14} & \text{'NegativeElectrode.Coating.ActiveMaterial1.Interface.guestStoichiometry0'} & \text{0.1429} \\ \hline\end{array}$

For the remaining properties, we load a standard data set

jsonstruct_cell = parseBattmoJson('ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_graphite.json');

We remove from this structure active material field. This step is not necessary but is cleaner and we avoid a warning.

jsonstruct_cell = removeStructField(jsonstruct_cell, {'NegativeElectrode', 'Coating', 'ActiveMaterial'});

we load a 1d geometry

jsonfilename = fullfile('Examples', 'JsonDataFiles', 'geometry1d.json');
jsonstruct_geometry = parseBattmoJson(jsonfilename);

We merge the json structures

jsonstruct = mergeStructs({jsonstruct_composite_material, ...
                               jsonstruct_cell              , ...
                               jsonstruct_geometry});

We do not consider the thermal model and remove the current collector. We also use a CV switch control.

jsonstruct.use_thermal                = false;
jsonstruct.include_current_collectors = false;

We define some shorcuts for the sub-models, for convenience

ne   = 'NegativeElectrode';
pe   = 'PositiveElectrode';
co   = 'Coating';
am1  = 'ActiveMaterial1';
am2  = 'ActiveMaterial2';
bd   = 'Binder';
ad   = 'ConductingAdditive';
sd   = 'SolidDiffusion';
itf  = 'Interface';
ctrl = 'Control';

We modify some parameters

We adjust the mass fractions parameters of the active material in the negative electrode

jsonstruct = setStructField(jsonstruct, {ne, co, am1, 'massFraction'}, 0.9, 'handleMisMatch', 'quiet');
jsonstruct = setStructField(jsonstruct, {ne, co, am2, 'massFraction'}, 0.08, 'handleMisMatch', 'quiet');
jsonstruct = setStructField(jsonstruct, {ne, co, bd , 'massFraction'}, 0.01, 'handleMisMatch', 'quiet');
jsonstruct = setStructField(jsonstruct, {ne, co, ad , 'massFraction'}, 0.01, 'handleMisMatch', 'quiet');

We run the simulations

output = runBatteryJson(jsonstruct);

Plotting

We extract the voltage, current and time from the simulation output

states = output.states;
model  = output.model;


E    = cellfun(@(x) x.Control.E, states);
I    = cellfun(@(x) x.Control.I, states);
time = cellfun(@(x) x.time, states);

We plot the voltage and current

figure
subplot(2, 1, 1);
plot(time/hour, E);
xlabel('Time / h');
ylabel('Voltage / V');
title('Voltage')
subplot(2, 1, 2);
plot(time/hour, I/milli);
xlabel('Time / h');
ylabel('Current / mA');
title('Current')

We compute and plot the state of charges in the different material

figure
hold on


for istate = 1 : numel(states)
    states{istate} = model.evalVarName(states{istate}, {ne, co, 'SOC'});
end


SOC  = cellfun(@(x) x.(ne).(co).SOC, states);
SOC1 = cellfun(@(x) x.(ne).(co).(am1).SOC, states);
SOC2 = cellfun(@(x) x.(ne).(co).(am2).SOC, states);


plot(time/hour, SOC, 'displayname', 'SOC - cumulated');
plot(time/hour, SOC1, 'displayname', 'SOC - Graphite');
plot(time/hour, SOC2, 'displayname', 'SOC - Silicon');


xlabel('Time / h');
ylabel('SOC / -');
title('SOCs')


legend show

plot of the particle concentration distribution in the particle at the end time

We recover the state at the last time step

state = states{end};

We iterate over the two active materials. The first one is the graphite and the second one the silicon

ams = {am1, am2};


for iam = 1 : numel(ams)


    am = ams{iam};


    model_sd = model.(ne).(co).(am).(sd);
    state_sd = state.(ne).(co).(am).(sd);

We recover the concentration as an array. The column index corresponds to the spatial direction (here x as we are considering a 1D model) and the row index corresponds to the particle radia direction

    c = model_sd.getParticleConcentrations(state_sd);
    r = model_sd.operators.radii;


    figure
    hold on

We plot the concentration distribution at the last point in the grid, which corresponds in this case to the closest to the positive electrode

    plot(r/(micro*meter), c(size(c, 1), :)/(mol/litre));
    title(sprintf('Particle concentration profile in %s', am));
    xlabel('radius / m');
    ylabel('concentration / mol/litre');

end

Charge step

initstate = states{end};
jsonstruct.(ctrl).CRate = 1;
jsonstruct = setStructField(jsonstruct, {'Control', 'controlPolicy'}, 'CCCharge', 'handleMisMatch', 'quiet');


jsonstruct.initializationSetup = 'given matlab object';


output = runBatteryJson(jsonstruct, 'initstate', initstate);

Visualisation

We concatenate the states we have computed

dischargeStates = states; % see previous assignement
chargeStates    = output.states; % assigned from last simulation output


allStates = vertcat(dischargeStates, chargeStates);

% We extract the voltage, current and time from the simulation output
E    = cellfun(@(x) x.Control.E, allStates);
I    = cellfun(@(x) x.Control.I, allStates);
time = cellfun(@(x) x.time, allStates);

We plot the voltage and current

figure
subplot(2, 1, 1);
plot(time/hour, E);
xlabel('Time / h');
ylabel('Voltage / V');
title('Voltage')
subplot(2, 1, 2);
plot(time/hour, I/milli);
xlabel('Time / h');
ylabel('Current / mA');
title('Current')

We compute and plot the state of charges in the different material

figure
hold on


for istate = 1 : numel(allStates)
    allStates{istate} = model.evalVarName(allStates{istate}, {ne, co, 'SOC'});
end


SOC  = cellfun(@(x) x.(ne).(co).SOC, allStates);
SOC1 = cellfun(@(x) x.(ne).(co).(am1).SOC, allStates);
SOC2 = cellfun(@(x) x.(ne).(co).(am2).SOC, allStates);


plot(time/hour, SOC, 'displayname', 'SOC - cumulated');
plot(time/hour, SOC1, 'displayname', 'SOC - Graphite');
plot(time/hour, SOC2, 'displayname', 'SOC - Silicon');


xlabel('Time / h');
ylabel('SOC / -');
title('SOCs')


legend show