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.
[1]:
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.
[2]:
flattenJsonStruct(jsonstruct_composite_material);
[2]:
parameter name parameter value
__________________________________________________________________________________________ _______________________________
{'NegativeElectrode.Coating.activeMaterialModelSetup.composite' } {[ 1]}
{'NegativeElectrode.Coating.ActiveMaterial1.specificHeatCapacity' } {[ 632]}
{'NegativeElectrode.Coating.ActiveMaterial1.heatCapacity' } {[ 632000]}
{'NegativeElectrode.Coating.ActiveMaterial1.thermalConductivity' } {[ 1.0400]}
{'NegativeElectrode.Coating.ActiveMaterial1.electronicConductivity' } {[ 100]}
{'NegativeElectrode.Coating.ActiveMaterial1.massFraction' } {[ 0.9200]}
{'NegativeElectrode.Coating.ActiveMaterial1.density' } {[ 2240]}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.saturationConcentration' } {[ 30555]}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.volumetricSurfaceArea' } {[ 723600]}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.numberOfElectronsTransferred' } {[ 1]}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.activationEnergyOfReaction' } {[ 5000]}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.reactionRateConstant' } {[ 5.0310e-11]}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.guestStoichiometry100' } {[ 0.8855]}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.guestStoichiometry0' } {[ 0.1429]}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.chargeTransferCoefficient' } {[ 0.5000]}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.openCircuitPotential.type' } {'function' }
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.openCircuitPotential.functionname' } {'computeOCP_Graphite_Torchio'}
{'NegativeElectrode.Coating.ActiveMaterial1.Interface.openCircuitPotential.argumentlist' } {3x1 cell }
{'NegativeElectrode.Coating.ActiveMaterial1.diffusionModelType' } {'full' }
{'NegativeElectrode.Coating.ActiveMaterial1.SolidDiffusion.activationEnergyOfDiffusion' } {[ 5000]}
{'NegativeElectrode.Coating.ActiveMaterial1.SolidDiffusion.referenceDiffusionCoefficient'} {[ 3.9000e-14]}
{'NegativeElectrode.Coating.ActiveMaterial1.SolidDiffusion.particleRadius' } {[ 1.0000e-06]}
{'NegativeElectrode.Coating.ActiveMaterial1.SolidDiffusion.N' } {[ 10]}
{'NegativeElectrode.Coating.ActiveMaterial2.specificHeatCapacity' } {[ 632]}
{'NegativeElectrode.Coating.ActiveMaterial2.heatCapacity' } {[ 632000]}
{'NegativeElectrode.Coating.ActiveMaterial2.thermalConductivity' } {[ 1.0400]}
{'NegativeElectrode.Coating.ActiveMaterial2.electronicConductivity' } {[ 100]}
{'NegativeElectrode.Coating.ActiveMaterial2.massFraction' } {[ 0.0800]}
{'NegativeElectrode.Coating.ActiveMaterial2.density' } {[ 2330]}
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.saturationConcentration' } {[ 3.1131e+05]}
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.volumetricSurfaceArea' } {[ 723600]}
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.numberOfElectronsTransferred' } {[ 1]}
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.activationEnergyOfReaction' } {[ 5000]}
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.reactionRateConstant' } {[ 5.0310e-11]}
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.guestStoichiometry100' } {[ 0.8855]}
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.guestStoichiometry0' } {[ 0.1429]}
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.chargeTransferCoefficient' } {[ 0.5000]}
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.openCircuitPotential.type' } {'function' }
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.openCircuitPotential.functionname' } {'computeOCP_Silicon_Li2012' }
{'NegativeElectrode.Coating.ActiveMaterial2.Interface.openCircuitPotential.argumentlist' } {3x1 cell }
{'NegativeElectrode.Coating.ActiveMaterial2.diffusionModelType' } {'full' }
{'NegativeElectrode.Coating.ActiveMaterial2.SolidDiffusion.activationEnergyOfDiffusion' } {[ 5000]}
{'NegativeElectrode.Coating.ActiveMaterial2.SolidDiffusion.referenceDiffusionCoefficient'} {[ 1.0000e-15]}
{'NegativeElectrode.Coating.ActiveMaterial2.SolidDiffusion.particleRadius' } {[ 6.0000e-08]}
{'NegativeElectrode.Coating.ActiveMaterial2.SolidDiffusion.N' } {[ 10]}
For the remaining properties, we load a standard data set
[3]:
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.
[4]:
jsonstruct_cell = removeJsonStructField(jsonstruct_cell, {'NegativeElectrode', 'Coating', 'ActiveMaterial'});
we load a 1d geometry
[5]:
jsonfilename = fullfile('Examples', 'JsonDataFiles', 'geometry1d.json');
jsonstruct_geometry = parseBattmoJson(jsonfilename);
We merge the json structures
[6]:
jsonstruct = mergeJsonStructs({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.
[7]:
jsonstruct.use_thermal = false;
jsonstruct.include_current_collectors = false;
We define some shorcuts for the sub-models, for convenience
[8]:
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
[9]:
jsonstruct = setJsonStructField(jsonstruct, {ne, co, am1, 'massFraction'}, 0.9, 'handleMisMatch', 'quiet');
jsonstruct = setJsonStructField(jsonstruct, {ne, co, am2, 'massFraction'}, 0.08, 'handleMisMatch', 'quiet');
jsonstruct = setJsonStructField(jsonstruct, {ne, co, bd , 'massFraction'}, 0.01, 'handleMisMatch', 'quiet');
jsonstruct = setJsonStructField(jsonstruct, {ne, co, ad , 'massFraction'}, 0.01, 'handleMisMatch', 'quiet');
We run the simulations
[10]:
output = runBatteryJson(jsonstruct);
Plotting
We extract the voltage, current and time from the simulation output
[11]:
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
[12]:
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')
[12]:
We compute and plot the state of charges in the different material
[13]:
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
[13]:
plot of the particle concentration distribution in the particle at the end time
We recover the state at the last time step
[14]:
state = states{end};
We iterate over the two active materials. The first one is the graphite and the second one the silicon
[15]:
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
[16]:
c = model_sd.getParticleConcentrations(state_sd);
r = model_sd.operators.radii;
figure
hold on
[16]:
We plot the concentration distribution at the last point in the grid, which corresponds in this case to the closest to the positive electrode
[17]:
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
[18]:
initstate = states{end};
jsonstruct.(ctrl).CRate = 1;
jsonstruct = setJsonStructField(jsonstruct, {'Control', 'controlPolicy'}, 'CCCharge', 'handleMisMatch', 'quiet');
jsonstruct.initializationSetup = 'given matlab object';
output = runBatteryJson(jsonstruct, 'initstate', initstate);
Visualisation
We concatenate the states we have computed
[19]:
dischargeStates = states; % see previous assignement
chargeStates = output.states; % assigned from last simulation output
allStates = vertcat(dischargeStates, chargeStates);
[20]:
% 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
[21]:
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')
[21]:
We compute and plot the state of charges in the different material
[22]:
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
[22]:
