diff --git a/_images/battMoTutorial_01.png b/_images/battMoTutorial_01.png index 736664e0..6bc56f83 100644 Binary files a/_images/battMoTutorial_01.png and b/_images/battMoTutorial_01.png differ diff --git a/_images/battMoTutorial_02.png b/_images/battMoTutorial_02.png index 747fc431..d58215d6 100644 Binary files a/_images/battMoTutorial_02.png and b/_images/battMoTutorial_02.png differ diff --git a/_images/runBatteryP2D_01.png b/_images/runBatteryP2D_01.png index c36027b9..4239be97 100644 Binary files a/_images/runBatteryP2D_01.png and b/_images/runBatteryP2D_01.png differ diff --git a/_images/runElectrolyser_01.png b/_images/runElectrolyser_01.png index 44bcbb3b..396ce403 100644 Binary files a/_images/runElectrolyser_01.png and b/_images/runElectrolyser_01.png differ diff --git a/_images/runElectrolyser_02.png b/_images/runElectrolyser_02.png deleted file mode 100644 index 94251d73..00000000 Binary files a/_images/runElectrolyser_02.png and /dev/null differ diff --git a/_images/runSEIActiveMaterial_01.png b/_images/runSEIActiveMaterial_01.png index 33e677a8..666d056c 100644 Binary files a/_images/runSEIActiveMaterial_01.png and b/_images/runSEIActiveMaterial_01.png differ diff --git a/_images/runSEIActiveMaterial_02.png b/_images/runSEIActiveMaterial_02.png index 867bab50..602793e5 100644 Binary files a/_images/runSEIActiveMaterial_02.png and b/_images/runSEIActiveMaterial_02.png differ diff --git a/_images/runSEIActiveMaterial_03.png b/_images/runSEIActiveMaterial_03.png index 26dab1b5..68c623e9 100644 Binary files a/_images/runSEIActiveMaterial_03.png and b/_images/runSEIActiveMaterial_03.png differ diff --git a/_images/runSEIActiveMaterial_04.png b/_images/runSEIActiveMaterial_04.png index dd9aeb37..c8a378a4 100644 Binary files a/_images/runSEIActiveMaterial_04.png and b/_images/runSEIActiveMaterial_04.png differ diff --git a/_images/runSEIActiveMaterial_05.png b/_images/runSEIActiveMaterial_05.png index f26f3609..3d011d0b 100644 Binary files a/_images/runSEIActiveMaterial_05.png and b/_images/runSEIActiveMaterial_05.png differ diff --git a/_images/runSEIActiveMaterial_06.png b/_images/runSEIActiveMaterial_06.png index ee7e6c06..911a33e3 100644 Binary files a/_images/runSEIActiveMaterial_06.png and b/_images/runSEIActiveMaterial_06.png differ diff --git a/_sources/compositeElectrode.rst.txt b/_sources/compositeElectrode.rst.txt index 6cb21217..e94606e1 100644 --- a/_sources/compositeElectrode.rst.txt +++ b/_sources/compositeElectrode.rst.txt @@ -8,7 +8,7 @@ Composite electrode json input file Here is json input for the composite electrode part, in :ref:`runSiliconGraphiteBattery` -.. literalinclude:: ../ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_silicon_graphite.json +.. literalinclude:: ../ParameterData/BatteryCellParameters/LithiumIonBatteryCell/composite_silicon_graphite.json :language: json diff --git a/_sources/intermediate.rst.txt b/_sources/intermediate.rst.txt index 865f901a..df9dada6 100644 --- a/_sources/intermediate.rst.txt +++ b/_sources/intermediate.rst.txt @@ -173,7 +173,7 @@ the :code:`jsonstruct` that is obtained is equivalent to the one where we would "guestStoichiometry0": 0.1429, "chargeTransferCoefficient": 0.5, "openCircuitPotential" : {"type": "function", - "functionname" : "computeOCP_graphite", + "functionname" : "computeOCP_Graphite_Torchio", "argumentlist" : ["cElectrode", "T", "cmax"] }}}, diff --git a/_sources/modelinitialisation.rst.txt b/_sources/modelinitialisation.rst.txt index 7792836d..bbae5806 100644 --- a/_sources/modelinitialisation.rst.txt +++ b/_sources/modelinitialisation.rst.txt @@ -261,12 +261,12 @@ Then, using the :code:`printSpecifications` method, we get Energy Density : 851.122 [Wh/L] Initial Voltage : 4.17686 [V] -We can also mention here the utility function :battmo:`computeCellEnergyGivenCrate`, even it is not a *static* property. +We can also mention here the utility function :battmo:`computeCellEnergyGivenDrate`, even it is not a *static* property. The function computes the energy produced by a cell for a given CRate. .. code:: matlab - >> output = computeCellEnergyGivenCrate(model, 2); + >> output = computeCellEnergyGivenDrate(model, 2); >> fprintf('Energy at Crate=2 : %g [Wh]', output.energy / (watt*hour)); Energy at Crate = 2 : 0.0110781 [Wh] diff --git a/_sources/parsets.rst.txt b/_sources/parsets.rst.txt index 31078240..062769a0 100644 --- a/_sources/parsets.rst.txt +++ b/_sources/parsets.rst.txt @@ -20,10 +20,10 @@ listed in the directory :battmofile:`ParameterSets` - 2009 - :cite:`Safari_2009` - :battmofile:`Safari2009` - * - Lin et al + * - Xu et al - 2015 - :cite:`LIN2015633` - - :battmofile:`Lin2015` + - :battmofile:`Xu2015` Otherwise, separate input json files corresponding to different part of a battery model are included in the directory diff --git a/_sources/publishedExamples/battMoTutorial.rst.txt b/_sources/publishedExamples/battMoTutorial.rst.txt index 4b4eabc9..9c0f4d2b 100644 --- a/_sources/publishedExamples/battMoTutorial.rst.txt +++ b/_sources/publishedExamples/battMoTutorial.rst.txt @@ -1,238 +1,227 @@ - -.. _battMoTutorial: - -BattMo Tutorial ------------------------------------- -*Generated from battMoTutorial.m* - - -This tutorial explains how to setup and run a simulation in BattMo - -.. code-block:: matlab - - % First clear variables stored in memory and close all figures - clear; - close all; - - % Set thicker line widths and larger font sizes - set(0, 'defaultLineLineWidth', 3); - set(0, 'defaultAxesFontSize', 16); - set(0, 'defaultTextFontSize', 18); - - -Setting up the environment -^^^^^^^^^^^^^^^^^^^^^^^^^^ -BattMo uses functionality from :mod:`MRST `. This functionality is collected into modules where each module contains code for doing specific things. To use this functionality we must add these modules to the matlab path by running: - -.. code-block:: matlab - - mrstModule add ad-core mrst-gui mpfa agmg linearsolvers - - -Specifying the physical model -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -In this tutorial we will simulate a lithium-ion battery consisting of a negative electrode, a positive electrode and an electrolyte. BattMo comes with some pre-defined models which can be loaded from JSON files. Here we will load the basic lithium-ion model JSON file which comes with Battmo. We use :battmo:`parseBattmoJson` to parse the file. - -.. code-block:: matlab - - fname = fullfile('ParameterData','BatteryCellParameters',... - 'LithiumIonBatteryCell','lithium_ion_battery_nmc_graphite.json'); - jsonstruct = parseBattmoJson(fname); - -The parseBattmoJson function parses the JSON input and creates a matlab structure containing the same fields as the JSON input. This structure can be changed to setup the model in the way that we want. -In this instance we will exclude temperature effects by setting use_thermal to false. - -.. code-block:: matlab - - jsonstruct.use_thermal = false; - -We will also not use current collectors in this example: - -.. code-block:: matlab - - jsonstruct.include_current_collectors = false; - -The structure created in the jsonstruct follows the same hierarchy as the fields in the JSON input file. These can be referenced by name in the jsonstruct. To make life easier for ourselves we define some shorthand names for various parts of the structure. - -.. code-block:: matlab - - ne = 'NegativeElectrode'; - pe = 'PositiveElectrode'; - elyte = 'Electrolyte'; - thermal = 'ThermalModel'; - co = 'Coating'; - am = 'ActiveMaterial'; - itf = 'Interface'; - sd = 'SolidDiffusion'; - ctrl = 'Control'; - cc = 'CurrentCollector'; - -Now we can set the diffusion model type for the active material (am) in the positive (pe) and negative (ne) electrodes to 'full'. - -.. code-block:: matlab - - jsonstruct.(pe).(am).diffusionModelType = 'full'; - jsonstruct.(ne).(am).diffusionModelType = 'full'; - -To see which other types of diffusion model are available one can view :battmo:`ActiveMaterialInputParams`. When running a simulation, BattMo requires that all model parameters are stored in an instance of :battmo:`BatteryInputParams`. This class is used to initialize the simulation and is accessed by various parts of the simulator during the simulation. This class is instantiated using the jsonstruct we just created: - -.. code-block:: matlab - - inputparams = BatteryInputParams(jsonstruct); - -It is also possible to update the properties of this inputparams in a similar way to updating the jsonstruct. Here we set the discretisation level for the diffusion model. Other input parameters for the full diffusion model can be found here: :battmo:`FullSolidDiffusionModelInputParams`. - -.. code-block:: matlab - - inputparams.(ne).(co).(am).(sd).N = 5; - inputparams.(pe).(co).(am).(sd).N = 5; - - % We can also change how the battery is operated, for example setting - % the cut off voltage. - inputparams.(ctrl).lowerCutoffVoltage = 2.5; - - -Setting up the geometry -^^^^^^^^^^^^^^^^^^^^^^^ -Here, we setup the 1D computational grid that will be used for the simulation. The required discretization parameters are already included in the class :battmo:`BatteryGeneratorP2D`. Classes for generating other geometries can be found in the BattMo/Battery/BatteryGeometry folder. - -.. code-block:: matlab - - gen = BatteryGeneratorP2D(); - -Now, we update the inputparams with the properties of the grid. This function will update relevent parameters in the inputparams object and make sure we have all the required parameters for the model geometry chosen. - -.. code-block:: matlab - - inputparams = gen.updateBatteryInputParams(inputparams); - - -Initialising the battery model object -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The battery model is initialized by sending inputparams to the Battery class constructor. see :battmo:`Battery`. -In BattMo a battery model is actually a collection of submodels: Electrolyte, Negative Electrode, Positive Electrode, Thermal Model and Control Model. The battery class contains all of these submodels and various other parameters necessary to run the simulation. - -.. code-block:: matlab - - model = Battery(inputparams); - - -Plotting the OCP curves against state of charge -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -We can inspect the model object to find out which parameters are being used. For instance the information we need to plot the OCP curves for the positive and negative electrodes can be found in the interface structure of each electrode. - -.. code-block:: matlab - - T = 298.15; - eldes = {ne, pe}; - - figure - hold on - - for ielde = 1:numel(eldes) - el_itf = model.(eldes{ielde}).(co).(am).(itf); - - theta100 = el_itf.guestStoichiometry100; - theta0 = el_itf.guestStoichiometry0; - cmax = el_itf.saturationConcentration; - - soc = linspace(0, 1); - theta = soc*theta100 + (1 - soc)*theta0; - c = theta.*cmax; - OCP = el_itf.computeOCPFunc(c, T, cmax); - - plot(soc, OCP) - end - - xlabel('SOC / -') - ylabel('OCP / V') - title('OCP for both electrodes'); - ylim([0, 5.5]) - legend(eldes, 'location', 'nw') - -.. figure:: battMoTutorial_01.png - :figwidth: 100% - - -Controlling the simulation -^^^^^^^^^^^^^^^^^^^^^^^^^^ -The control model specifies how the battery is operated, i.e., how the simulation is controlled. -In the first instance we use CCDischarge control policy. We set the total time scaled by the CRate in the model. The CRate has been set by the json file. We can access it here: - -.. code-block:: matlab - - CRate = model.Control.CRate; - total = 1.1*hour/CRate; - -We want to break this total time into 100 timesteps. To begin with we will use equal values for each timestep. -We create a structure containing the length of each step in seconds ('val') and also which control to use for each step ('control'). -In this case we use control 1 for all steps. This means that the functions used to setup the control values are the same at each step. - -.. code-block:: matlab - - n = 100; - dt = total/n; - step = struct('val', dt*ones(n, 1), 'control', ones(n, 1)); - -We create a control structure containing the source function and and a stopping criteria. The control parameters have been given in the json file :battmofile:`ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_graphite.json` -The :code:`setupScheduleControl` method contains the code to setup the control structure that is used in the schedule structure setup below. - -.. code-block:: matlab - - control = model.Control.setupScheduleControl(); - -Finally we collect the control and step structures together in a schedule struct which is the schedule which the simulation will follow: - -.. code-block:: matlab - - schedule = struct('control', control, 'step', step); - - -Setting the initial state of the battery -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -To run simulation we need to know the starting point which we will run it from, in terms of the value of the primary variables being modelled at the start of the simulation. The initial state of the model is setup using model.setupInitialState() Here we take the state of charge (SOC) given in the input and calculate equilibrium concentration based on theta0, theta100 and cmax. - -.. code-block:: matlab - - initstate = model.setupInitialState(); - - -Running the simulation -^^^^^^^^^^^^^^^^^^^^^^ -Once we have the initial state, the model and the schedule, we can call the simulateScheduleAD function which will actually run the simulation. -The outputs from the simulation are: - sols: which provides the current and voltage of the battery at each timestep. - states: which contains the values of the primary variables in the model at each timestep. - reports: which contains technical information about the steps used in the numerical solvers. - -.. code-block:: matlab - - [sols, states, report] = simulateScheduleAD(initstate, model, schedule); - - -Plotting the results -^^^^^^^^^^^^^^^^^^^^ -To get the results we use the matlab cellfun function to extract the values Control.E, Control.I and time from each timestep (cell in the cell array) in states. We can then plot the vectors. - -.. code-block:: matlab - - E = cellfun(@(x) x.Control.E, states); - I = cellfun(@(x) x.Control.I, states); - time = cellfun(@(x) x.time, states); - - figure() - - subplot(1,2,1) - plot(time/hour, E) - xlabel('time / h') - ylabel('Cell Voltage / V') - - subplot(1,2,2) - plot(time/hour, I) - ylim([0, 0.02]) - xlabel('time / h') - ylabel('Cell Current / A') - -.. figure:: battMoTutorial_02.png - :figwidth: 100% - - - -complete source code can be found :ref:`here` + +.. _battMoTutorial: + +BattMo Tutorial +------------------------------------ +*Generated from battMoTutorial.m* + + +This tutorial explains how to setup and run a simulation in BattMo + +.. code-block:: matlab + + % First clear variables stored in memory and close all figures + clear; + close all; + + % Set thicker line widths and larger font sizes + set(0, 'defaultLineLineWidth', 3); + set(0, 'defaultAxesFontSize', 16); + set(0, 'defaultTextFontSize', 18); + + +Setting up the environment +^^^^^^^^^^^^^^^^^^^^^^^^^^ +BattMo uses functionality from :mod:`MRST `. This functionality is collected into modules where each module contains code for doing specific things. To use this functionality we must add these modules to the matlab path by running: + +.. code-block:: matlab + + mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + + +Specifying the physical model +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +In this tutorial we will simulate a lithium-ion battery consisting of a negative electrode, a positive electrode and an electrolyte. BattMo comes with some pre-defined models which can be loaded from JSON files. Here we will load the basic lithium-ion model JSON file which comes with Battmo. We use :battmo:`parseBattmoJson` to parse the file. + +.. code-block:: matlab + + fname = fullfile('ParameterData','BatteryCellParameters',... + 'LithiumIonBatteryCell','lithium_ion_battery_nmc_graphite.json'); + jsonstruct = parseBattmoJson(fname); + +The parseBattmoJson function parses the JSON input and creates a matlab structure containing the same fields as the JSON input. This structure can be changed to setup the model in the way that we want. +In this instance we will exclude temperature effects by setting use_thermal to false. + +.. code-block:: matlab + + jsonstruct.use_thermal = false; + +We will also not use current collectors in this example: + +.. code-block:: matlab + + jsonstruct.include_current_collectors = false; + +The structure created in the jsonstruct follows the same hierarchy as the fields in the JSON input file. These can be referenced by name in the jsonstruct. To make life easier for ourselves we define some shorthand names for various parts of the structure. + +.. code-block:: matlab + + ne = 'NegativeElectrode'; + pe = 'PositiveElectrode'; + elyte = 'Electrolyte'; + thermal = 'ThermalModel'; + co = 'Coating'; + am = 'ActiveMaterial'; + itf = 'Interface'; + sd = 'SolidDiffusion'; + ctrl = 'Control'; + cc = 'CurrentCollector'; + +Now we can set the diffusion model type for the active material (am) in the positive (pe) and negative (ne) electrodes to 'full'. + +.. code-block:: matlab + + jsonstruct.(pe).(am).diffusionModelType = 'full'; + jsonstruct.(ne).(am).diffusionModelType = 'full'; + +To see which other types of diffusion model are available one can view :battmo:`ActiveMaterialInputParams`. When running a simulation, BattMo requires that all model parameters are stored in an instance of :battmo:`BatteryInputParams`. This class is used to initialize the simulation and is accessed by various parts of the simulator during the simulation. This class is instantiated using the jsonstruct we just created: + +.. code-block:: matlab + + inputparams = BatteryInputParams(jsonstruct); + +It is also possible to update the properties of this inputparams in a similar way to updating the jsonstruct. Here we set the discretisation level for the diffusion model. Other input parameters for the full diffusion model can be found here: :battmo:`FullSolidDiffusionModelInputParams`. + +.. code-block:: matlab + + inputparams.(ne).(co).(am).(sd).N = 5; + inputparams.(pe).(co).(am).(sd).N = 5; + + % We can also change how the battery is operated, for example setting + % the cut off voltage. + inputparams.(ctrl).lowerCutoffVoltage = 2.5; + + +Setting up the geometry +^^^^^^^^^^^^^^^^^^^^^^^ +Here, we setup the 1D computational grid that will be used for the simulation. The required discretization parameters are already included in the class :battmo:`BatteryGeneratorP2D`. Classes for generating other geometries can be found in the BattMo/Battery/BatteryGeometry folder. + +.. code-block:: matlab + + gen = BatteryGeneratorP2D(); + +Now, we update the inputparams with the properties of the grid. This function will update relevent parameters in the inputparams object and make sure we have all the required parameters for the model geometry chosen. + +.. code-block:: matlab + + inputparams = gen.updateBatteryInputParams(inputparams); + + +Initialising the battery model object +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The battery model is initialized by sending inputparams to the Battery class constructor. see :battmo:`Battery`. +In BattMo a battery model is actually a collection of submodels: Electrolyte, Negative Electrode, Positive Electrode, Thermal Model and Control Model. The battery class contains all of these submodels and various other parameters necessary to run the simulation. + +.. code-block:: matlab + + model = Battery(inputparams); + + +Plotting the OCP curves against state of charge +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +We can inspect the model object to find out which parameters are being used. For instance the information we need to plot the OCP curves for the positive and negative electrodes can be found in the interface structure of each electrode. + +.. code-block:: matlab + + T = 298.15; + eldes = {ne, pe}; + + figure + hold on + + for ielde = 1:numel(eldes) + el_itf = model.(eldes{ielde}).(co).(am).(itf); + + theta100 = el_itf.guestStoichiometry100; + theta0 = el_itf.guestStoichiometry0; + cmax = el_itf.saturationConcentration; + + soc = linspace(0, 1); + theta = soc*theta100 + (1 - soc)*theta0; + c = theta.*cmax; + OCP = el_itf.computeOCPFunc(c, T, cmax); + + plot(soc, OCP) + end + + xlabel('SOC / -') + ylabel('OCP / V') + title('OCP for both electrodes'); + ylim([0, 5.5]) + legend(eldes, 'location', 'nw') + +.. figure:: battMoTutorial_01.png + :figwidth: 100% + + +Controlling the simulation +^^^^^^^^^^^^^^^^^^^^^^^^^^ +The control model specifies how the battery is operated, i.e., how the simulation is controlled. +The input parameters for the control have been given as part of the json structure :battmofile:`ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_graphite.json`. The total simulation time is setup for us, computed from the CRate value. We use the method :code:`setupScheduleStep` in :battmo:`ControlModel` to setup the :code:`step` structure. + +.. code-block:: matlab + + step = model.Control.setupScheduleStep(); + +We create a control structure containing the source function and and a stopping criteria. The control parameters have been given in the json file :battmofile:`ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_graphite.json` +The :code:`setupScheduleControl` method contains the code to setup the control structure that is used in the schedule structure setup below. + +.. code-block:: matlab + + control = model.Control.setupScheduleControl(); + +Finally we collect the control and step structures together in a schedule struct which is the schedule which the simulation will follow: + +.. code-block:: matlab + + schedule = struct('control', control, 'step', step); + + +Setting the initial state of the battery +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +To run simulation we need to know the starting point which we will run it from, in terms of the value of the primary variables being modelled at the start of the simulation. The initial state of the model is setup using model.setupInitialState() Here we take the state of charge (SOC) given in the input and calculate equilibrium concentration based on theta0, theta100 and cmax. + +.. code-block:: matlab + + initstate = model.setupInitialState(); + + +Running the simulation +^^^^^^^^^^^^^^^^^^^^^^ +Once we have the initial state, the model and the schedule, we can call the simulateScheduleAD function which will actually run the simulation. +The outputs from the simulation are: - sols: which provides the current and voltage of the battery at each timestep. - states: which contains the values of the primary variables in the model at each timestep. - reports: which contains technical information about the steps used in the numerical solvers. + +.. code-block:: matlab + + [sols, states, report] = simulateScheduleAD(initstate, model, schedule); + + +Plotting the results +^^^^^^^^^^^^^^^^^^^^ +To get the results we use the matlab cellfun function to extract the values Control.E, Control.I and time from each timestep (cell in the cell array) in states. We can then plot the vectors. + +.. code-block:: matlab + + E = cellfun(@(x) x.Control.E, states); + I = cellfun(@(x) x.Control.I, states); + time = cellfun(@(x) x.time, states); + + figure() + + subplot(1,2,1) + plot(time/hour, E) + xlabel('time / h') + ylabel('Cell Voltage / V') + + subplot(1,2,2) + plot(time/hour, I) + ylim([0, 0.02]) + xlabel('time / h') + ylabel('Cell Current / A') + +.. figure:: battMoTutorial_02.png + :figwidth: 100% + + + +complete source code can be found :ref:`here` diff --git a/_sources/publishedExamples/battMoTutorial_source.rst.txt b/_sources/publishedExamples/battMoTutorial_source.rst.txt index 99549cd6..3272af5b 100644 --- a/_sources/publishedExamples/battMoTutorial_source.rst.txt +++ b/_sources/publishedExamples/battMoTutorial_source.rst.txt @@ -1,273 +1,258 @@ -:orphan: - -.. _battMoTutorial_source: - -Source code for battMoTutorial ------------------------------- - -.. code:: matlab - - - %% BattMo Tutorial - % This tutorial explains how to setup and run a simulation in BattMo - - % First clear variables stored in memory and close all figures - clear; - close all; - - % Set thicker line widths and larger font sizes - set(0, 'defaultLineLineWidth', 3); - set(0, 'defaultAxesFontSize', 16); - set(0, 'defaultTextFontSize', 18); - - %% Setting up the environment - % BattMo uses functionality from :mod:`MRST `. This functionality - % is collected into modules where each module contains code for doing - % specific things. To use this functionality we must add these modules to - % the matlab path by running: - - mrstModule add ad-core mrst-gui mpfa agmg linearsolvers - - %% Specifying the physical model - % In this tutorial we will simulate a lithium-ion battery consisting of a - % negative electrode, a positive electrode and an electrolyte. *BattMo* - % comes with some pre-defined models which can be loaded from JSON files. - % Here we will load the basic lithium-ion model JSON file which comes with - % Battmo. We use :battmo:`parseBattmoJson` to parse the file. - - fname = fullfile('ParameterData','BatteryCellParameters',... - 'LithiumIonBatteryCell','lithium_ion_battery_nmc_graphite.json'); - jsonstruct = parseBattmoJson(fname); - - %% - % The parseBattmoJson function parses the JSON input and creates a matlab - % structure containing the same fields as the JSON input. This structure - % can be changed to setup the model in the way that we want. - % - % In this instance we will exclude temperature effects by setting - % use_thermal to false. - - jsonstruct.use_thermal = false; - - %% - % We will also not use current collectors in this example: - - jsonstruct.include_current_collectors = false; - - %% - % The structure created in the jsonstruct follows the same hierarchy as the - % fields in the JSON input file. These can be referenced by name in the - % jsonstruct. To make life easier for ourselves we define some shorthand - % names for various parts of the structure. - - ne = 'NegativeElectrode'; - pe = 'PositiveElectrode'; - elyte = 'Electrolyte'; - thermal = 'ThermalModel'; - co = 'Coating'; - am = 'ActiveMaterial'; - itf = 'Interface'; - sd = 'SolidDiffusion'; - ctrl = 'Control'; - cc = 'CurrentCollector'; - - %% - % Now we can set the diffusion model type for the active material (am) in the - % positive (pe) and negative (ne) electrodes to 'full'. - - jsonstruct.(pe).(am).diffusionModelType = 'full'; - jsonstruct.(ne).(am).diffusionModelType = 'full'; - - %% - % To see which other types of diffusion model are available one can view - % :battmo:`ActiveMaterialInputParams`. When running a simulation, *BattMo* requires that all model parameters - % are stored in an instance of :battmo:`BatteryInputParams`. - % This class is used to initialize the simulation and is accessed by - % various parts of the simulator during the simulation. This class is - % instantiated using the jsonstruct we just created: - - inputparams = BatteryInputParams(jsonstruct); - - %% - % It is also possible to update the properties of this inputparams in a - % similar way to updating the jsonstruct. Here we set the discretisation - % level for the diffusion model. Other input parameters for the full diffusion - % model can be found here: - % :battmo:`FullSolidDiffusionModelInputParams`. - - inputparams.(ne).(co).(am).(sd).N = 5; - inputparams.(pe).(co).(am).(sd).N = 5; - - % We can also change how the battery is operated, for example setting - % the cut off voltage. - inputparams.(ctrl).lowerCutoffVoltage = 2.5; - - %% Setting up the geometry - % Here, we setup the 1D computational grid that will be used for the - % simulation. The required discretization parameters are already included - % in the class :battmo:`BatteryGeneratorP2D`. Classes for generating other geometries can - % be found in the BattMo/Battery/BatteryGeometry folder. - - gen = BatteryGeneratorP2D(); - - %% - % Now, we update the inputparams with the properties of the grid. This function - % will update relevent parameters in the inputparams object and make sure we have - % all the required parameters for the model geometry chosen. - - inputparams = gen.updateBatteryInputParams(inputparams); - - %% Initialising the battery model object - % The battery model is initialized by sending inputparams to the Battery class - % constructor. see :battmo:`Battery`. - % - % In BattMo a battery model is actually a collection of submodels: - % Electrolyte, Negative Electrode, Positive Electrode, Thermal Model and Control - % Model. The battery class contains all of these submodels and various other - % parameters necessary to run the simulation. - - model = Battery(inputparams); - - %% Plotting the OCP curves against state of charge - % We can inspect the model object to find out which parameters are being - % used. For instance the information we need to plot the OCP curves for the - % positive and negative electrodes can be found in the interface structure - % of each electrode. - - T = 298.15; - eldes = {ne, pe}; - - figure - hold on - - for ielde = 1:numel(eldes) - el_itf = model.(eldes{ielde}).(co).(am).(itf); - - theta100 = el_itf.guestStoichiometry100; - theta0 = el_itf.guestStoichiometry0; - cmax = el_itf.saturationConcentration; - - soc = linspace(0, 1); - theta = soc*theta100 + (1 - soc)*theta0; - c = theta.*cmax; - OCP = el_itf.computeOCPFunc(c, T, cmax); - - plot(soc, OCP) - end - - xlabel('SOC / -') - ylabel('OCP / V') - title('OCP for both electrodes'); - ylim([0, 5.5]) - legend(eldes, 'location', 'nw') - - %% Controlling the simulation - % The control model specifies how the battery is operated, i.e., how - % the simulation is controlled. - % - % In the first instance we use CCDischarge control policy. - % We set the total time scaled by the CRate in the model. - % The CRate has been set by the json file. We can access it here: - - CRate = model.Control.CRate; - total = 1.1*hour/CRate; - - %% - % We want to break this total time into 100 timesteps. To begin with we - % will use equal values for each timestep. - % - % We create a structure containing the length of each step in seconds - % ('val') and also which control to use for each step ('control'). - % - % In this case we use control 1 for all steps. This means that the functions - % used to setup the control values are the same at each step. - - n = 100; - dt = total/n; - step = struct('val', dt*ones(n, 1), 'control', ones(n, 1)); - - %% - % We create a control structure containing the source function and - % and a stopping criteria. The control parameters have been given in the json file - % :battmofile:`ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_graphite.json` - % - % The :code:`setupScheduleControl` method contains the code to setup the control structure that is used in the schedule structure setup below. - - control = model.Control.setupScheduleControl(); - - %% - % Finally we collect the control and step structures together in a schedule - % struct which is the schedule which the simulation will follow: - - schedule = struct('control', control, 'step', step); - - - %% Setting the initial state of the battery - % To run simulation we need to know the starting point which we will run it - % from, in terms of the value of the primary variables being modelled at - % the start of the simulation. - % The initial state of the model is setup using model.setupInitialState() - % Here we take the state of charge (SOC) given in the input and calculate - % equilibrium concentration based on theta0, theta100 and cmax. - - initstate = model.setupInitialState(); - - - %% Running the simulation - % Once we have the initial state, the model and the schedule, we can call - % the simulateScheduleAD function which will actually run the simulation. - % - % The outputs from the simulation are: - % - sols: which provides the current and voltage of the battery at each - % timestep. - % - states: which contains the values of the primary variables in the model - % at each timestep. - % - reports: which contains technical information about the steps used in - % the numerical solvers. - - [sols, states, report] = simulateScheduleAD(initstate, model, schedule); - - - %% Plotting the results - % To get the results we use the matlab cellfun function to extract the - % values Control.E, Control.I and time from each timestep (cell in the cell - % array) in states. We can then plot the vectors. - - E = cellfun(@(x) x.Control.E, states); - I = cellfun(@(x) x.Control.I, states); - time = cellfun(@(x) x.time, states); - - figure() - - subplot(1,2,1) - plot(time/hour, E) - xlabel('time / h') - ylabel('Cell Voltage / V') - - subplot(1,2,2) - plot(time/hour, I) - ylim([0, 0.02]) - xlabel('time / h') - ylabel('Cell Current / A') - - - %{ - Copyright 2021-2023 SINTEF Industry, Sustainable Energy Technology - and SINTEF Digital, Mathematics & Cybernetics. - - This file is part of The Battery Modeling Toolbox BattMo - - BattMo is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - BattMo is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with BattMo. If not, see . - %} - +:orphan: + +.. _battMoTutorial_source: + +Source code for battMoTutorial +------------------------------ + +.. code:: matlab + + + %% BattMo Tutorial + % This tutorial explains how to setup and run a simulation in BattMo + + % First clear variables stored in memory and close all figures + clear; + close all; + + % Set thicker line widths and larger font sizes + set(0, 'defaultLineLineWidth', 3); + set(0, 'defaultAxesFontSize', 16); + set(0, 'defaultTextFontSize', 18); + + %% Setting up the environment + % BattMo uses functionality from :mod:`MRST `. This functionality + % is collected into modules where each module contains code for doing + % specific things. To use this functionality we must add these modules to + % the matlab path by running: + + mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + + %% Specifying the physical model + % In this tutorial we will simulate a lithium-ion battery consisting of a + % negative electrode, a positive electrode and an electrolyte. *BattMo* + % comes with some pre-defined models which can be loaded from JSON files. + % Here we will load the basic lithium-ion model JSON file which comes with + % Battmo. We use :battmo:`parseBattmoJson` to parse the file. + + fname = fullfile('ParameterData','BatteryCellParameters',... + 'LithiumIonBatteryCell','lithium_ion_battery_nmc_graphite.json'); + jsonstruct = parseBattmoJson(fname); + + %% + % The parseBattmoJson function parses the JSON input and creates a matlab + % structure containing the same fields as the JSON input. This structure + % can be changed to setup the model in the way that we want. + % + % In this instance we will exclude temperature effects by setting + % use_thermal to false. + + jsonstruct.use_thermal = false; + + %% + % We will also not use current collectors in this example: + + jsonstruct.include_current_collectors = false; + + %% + % The structure created in the jsonstruct follows the same hierarchy as the + % fields in the JSON input file. These can be referenced by name in the + % jsonstruct. To make life easier for ourselves we define some shorthand + % names for various parts of the structure. + + ne = 'NegativeElectrode'; + pe = 'PositiveElectrode'; + elyte = 'Electrolyte'; + thermal = 'ThermalModel'; + co = 'Coating'; + am = 'ActiveMaterial'; + itf = 'Interface'; + sd = 'SolidDiffusion'; + ctrl = 'Control'; + cc = 'CurrentCollector'; + + %% + % Now we can set the diffusion model type for the active material (am) in the + % positive (pe) and negative (ne) electrodes to 'full'. + + jsonstruct.(pe).(am).diffusionModelType = 'full'; + jsonstruct.(ne).(am).diffusionModelType = 'full'; + + %% + % To see which other types of diffusion model are available one can view + % :battmo:`ActiveMaterialInputParams`. When running a simulation, *BattMo* requires that all model parameters + % are stored in an instance of :battmo:`BatteryInputParams`. + % This class is used to initialize the simulation and is accessed by + % various parts of the simulator during the simulation. This class is + % instantiated using the jsonstruct we just created: + + inputparams = BatteryInputParams(jsonstruct); + + %% + % It is also possible to update the properties of this inputparams in a + % similar way to updating the jsonstruct. Here we set the discretisation + % level for the diffusion model. Other input parameters for the full diffusion + % model can be found here: + % :battmo:`FullSolidDiffusionModelInputParams`. + + inputparams.(ne).(co).(am).(sd).N = 5; + inputparams.(pe).(co).(am).(sd).N = 5; + + % We can also change how the battery is operated, for example setting + % the cut off voltage. + inputparams.(ctrl).lowerCutoffVoltage = 2.5; + + %% Setting up the geometry + % Here, we setup the 1D computational grid that will be used for the + % simulation. The required discretization parameters are already included + % in the class :battmo:`BatteryGeneratorP2D`. Classes for generating other geometries can + % be found in the BattMo/Battery/BatteryGeometry folder. + + gen = BatteryGeneratorP2D(); + + %% + % Now, we update the inputparams with the properties of the grid. This function + % will update relevent parameters in the inputparams object and make sure we have + % all the required parameters for the model geometry chosen. + + inputparams = gen.updateBatteryInputParams(inputparams); + + %% Initialising the battery model object + % The battery model is initialized by sending inputparams to the Battery class + % constructor. see :battmo:`Battery`. + % + % In BattMo a battery model is actually a collection of submodels: + % Electrolyte, Negative Electrode, Positive Electrode, Thermal Model and Control + % Model. The battery class contains all of these submodels and various other + % parameters necessary to run the simulation. + + model = Battery(inputparams); + + %% Plotting the OCP curves against state of charge + % We can inspect the model object to find out which parameters are being + % used. For instance the information we need to plot the OCP curves for the + % positive and negative electrodes can be found in the interface structure + % of each electrode. + + T = 298.15; + eldes = {ne, pe}; + + figure + hold on + + for ielde = 1:numel(eldes) + el_itf = model.(eldes{ielde}).(co).(am).(itf); + + theta100 = el_itf.guestStoichiometry100; + theta0 = el_itf.guestStoichiometry0; + cmax = el_itf.saturationConcentration; + + soc = linspace(0, 1); + theta = soc*theta100 + (1 - soc)*theta0; + c = theta.*cmax; + OCP = el_itf.computeOCPFunc(c, T, cmax); + + plot(soc, OCP) + end + + xlabel('SOC / -') + ylabel('OCP / V') + title('OCP for both electrodes'); + ylim([0, 5.5]) + legend(eldes, 'location', 'nw') + + %% Controlling the simulation + % The control model specifies how the battery is operated, i.e., how the simulation is controlled. + % + % The input parameters for the control have been given as part of the json structure + % :battmofile:`ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_graphite.json`. The + % total simulation time is setup for us, computed from the CRate value. We use the method :code:`setupScheduleStep` in + % :battmo:`ControlModel` to setup the :code:`step` structure. + + step = model.Control.setupScheduleStep(); + + %% + % We create a control structure containing the source function and + % and a stopping criteria. The control parameters have been given in the json file + % :battmofile:`ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_graphite.json` + % + % The :code:`setupScheduleControl` method contains the code to setup the control structure that is used in the schedule + % structure setup below. + + control = model.Control.setupScheduleControl(); + + %% + % Finally we collect the control and step structures together in a schedule + % struct which is the schedule which the simulation will follow: + + schedule = struct('control', control, 'step', step); + + + %% Setting the initial state of the battery + % To run simulation we need to know the starting point which we will run it + % from, in terms of the value of the primary variables being modelled at + % the start of the simulation. + % The initial state of the model is setup using model.setupInitialState() + % Here we take the state of charge (SOC) given in the input and calculate + % equilibrium concentration based on theta0, theta100 and cmax. + + initstate = model.setupInitialState(); + + %% Running the simulation + % Once we have the initial state, the model and the schedule, we can call + % the simulateScheduleAD function which will actually run the simulation. + % + % The outputs from the simulation are: + % - sols: which provides the current and voltage of the battery at each + % timestep. + % - states: which contains the values of the primary variables in the model + % at each timestep. + % - reports: which contains technical information about the steps used in + % the numerical solvers. + + [sols, states, report] = simulateScheduleAD(initstate, model, schedule); + + + %% Plotting the results + % To get the results we use the matlab cellfun function to extract the + % values Control.E, Control.I and time from each timestep (cell in the cell + % array) in states. We can then plot the vectors. + + E = cellfun(@(x) x.Control.E, states); + I = cellfun(@(x) x.Control.I, states); + time = cellfun(@(x) x.time, states); + + figure() + + subplot(1,2,1) + plot(time/hour, E) + xlabel('time / h') + ylabel('Cell Voltage / V') + + subplot(1,2,2) + plot(time/hour, I) + ylim([0, 0.02]) + xlabel('time / h') + ylabel('Cell Current / A') + + + %{ + Copyright 2021-2024 SINTEF Industry, Sustainable Energy Technology + and SINTEF Digital, Mathematics & Cybernetics. + + This file is part of The Battery Modeling Toolbox BattMo + + BattMo is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + BattMo is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with BattMo. If not, see . + %} + diff --git a/_sources/publishedExamples/runBatteryP2D.rst.txt b/_sources/publishedExamples/runBatteryP2D.rst.txt index 353e04c3..a492395a 100644 --- a/_sources/publishedExamples/runBatteryP2D.rst.txt +++ b/_sources/publishedExamples/runBatteryP2D.rst.txt @@ -1,217 +1,203 @@ - -.. _runBatteryP2D: - -Pseudo-Two-Dimensional (P2D) Lithium-Ion Battery Model --------------------------------------------------------------------------- -*Generated from runBatteryP2D.m* - - -This example demonstrates how to setup a P2D model of a Li-ion battery and run a simple simulation. - -.. code-block:: matlab - - % Clear the workspace and close open figures - clear - close all - - -Import the required modules from MRST -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -load MRST modules - -.. code-block:: matlab - - mrstModule add ad-core mrst-gui mpfa agmg linearsolvers - - -Setup the properties of Li-ion battery materials and cell design -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The properties and parameters of the battery cell, including the architecture and materials, are set using an instance of :class:`BatteryInputParams `. This class is used to initialize the simulation and it propagates all the parameters throughout the submodels. The input parameters can be set manually or provided in json format. All the parameters for the model are stored in the inputparams object. - -.. code-block:: matlab - - jsonstruct = parseBattmoJson(fullfile('ParameterData','BatteryCellParameters','LithiumIonBatteryCell','lithium_ion_battery_nmc_graphite.json')); - - % We define some shorthand names for simplicity. - ne = 'NegativeElectrode'; - pe = 'PositiveElectrode'; - elyte = 'Electrolyte'; - thermal = 'ThermalModel'; - co = 'Coating'; - am = 'ActiveMaterial'; - itf = 'Interface'; - sd = 'SolidDiffusion'; - ctrl = 'Control'; - cc = 'CurrentCollector'; - - jsonstruct.use_thermal = false; - jsonstruct.include_current_collectors = false; - - inputparams = BatteryInputParams(jsonstruct); - - use_cccv = false; - if use_cccv - cccvstruct = struct( 'controlPolicy' , 'CCCV', ... - 'initialControl' , 'discharging', ... - 'CRate' , 1 , ... - 'lowerCutoffVoltage', 2.4 , ... - 'upperCutoffVoltage', 4.1 , ... - 'dIdtLimit' , 0.01 , ... - 'dEdtLimit' , 0.01); - cccvparamobj = CcCvControlModelInputParams(cccvstruct); - inputparams.Control = cccvinputparams; - end - - -Setup the geometry and computational grid -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Here, we setup the 1D computational grid that will be used for the simulation. The required discretization parameters are already included in the class BatteryGeneratorP2D. - -.. code-block:: matlab - - gen = BatteryGeneratorP2D(); - - % Now, we update the inputparams with the properties of the grid. - inputparams = gen.updateBatteryInputParams(inputparams); - - -Initialize the battery model. -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The battery model is initialized by sending inputparams to the Battery class constructor. see :class:`Battery `. - -.. code-block:: matlab - - model = Battery(inputparams); - - model.AutoDiffBackend= AutoDiffBackend(); - - inspectgraph = false; - if inspectgraph - cgt = model.computationalGraph; - return - end - - -Compute the nominal cell capacity and choose a C-Rate -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The nominal capacity of the cell is calculated from the active materials. This value is then combined with the user-defined C-Rate to set the cell operational current. - -.. code-block:: matlab - - CRate = model.Control.CRate; - - -Setup the time step schedule -^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Smaller time steps are used to ramp up the current from zero to its operational value. Larger time steps are then used for the normal operation. - -.. code-block:: matlab - - switch model.(ctrl).controlPolicy - case 'CCCV' - total = 3.5*hour/CRate; - case 'CCDischarge' - total = 1.4*hour/CRate; - otherwise - error('control policy not recognized'); - end - - n = 100; - dt = total/n; - step = struct('val', dt*ones(n, 1), 'control', ones(n, 1)); - - % we setup the control by assigning a source and stop function. - - control = model.Control.setupScheduleControl(); - - % This control is used to set up the schedule - schedule = struct('control', control, 'step', step); - - -Setup the initial state of the model -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The initial state of the model is setup using the model.setupInitialState() method. - -.. code-block:: matlab - - initstate = model.setupInitialState(); - - -Setup the properties of the nonlinear solver -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. code-block:: matlab - - nls = NonLinearSolver(); - - linearsolver = 'direct'; - switch linearsolver - case 'amgcl' - nls.LinearSolver = AMGCLSolverAD('verbose', true, 'reduceToCell', false); - nls.LinearSolver.tolerance = 1e-4; - nls.LinearSolver.maxIterations = 30; - nls.maxIterations = 10; - nls.verbose = 10; - case 'battery' - nls.LinearSolver = LinearSolverBatteryExtra('verbose' , false, ... - 'reduceToCell', true, ... - 'verbosity' , 3 , ... - 'reuse_setup' , false, ... - 'method' , 'direct'); - nls.LinearSolver.tolerance = 1e-4; - case 'direct' - disp('standard direct solver') - otherwise - error('Unknown solver %s', linearsolver); - end - - % Change default maximum iteration number in nonlinear solver - nls.maxIterations = 10; - % Change default behavior of nonlinear solver, in case of error - nls.errorOnFailure = false; - nls.timeStepSelector = StateChangeTimeStepSelector('TargetProps', {{'Control','E'}}, 'targetChangeAbs', 0.03); - % Change default tolerance for nonlinear solver - model.nonlinearTolerance = 1e-3*model.Control.Imax; - % Set verbosity - model.verbose = true; - - -Run the simulation -^^^^^^^^^^^^^^^^^^ - -.. code-block:: matlab - - [~, states, report] = simulateScheduleAD(initstate, model, schedule, 'OutputMinisteps', true, 'NonLinearSolver', nls); - - -Process output and recover the output voltage and current from the output states. -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. code-block:: matlab - - ind = cellfun(@(x) not(isempty(x)), states); - states = states(ind); - E = cellfun(@(x) x.Control.E, states); - I = cellfun(@(x) x.Control.I, states); - T = cellfun(@(x) max(x.(thermal).T), states); - Tmax = cellfun(@(x) max(x.ThermalModel.T), states); - % [SOCN, SOCP] = cellfun(@(x) model.calculateSOC(x), states); - time = cellfun(@(x) x.time, states); - - figure - plot(time/hour, E); - grid on - xlabel 'time / h'; - ylabel 'potential / V'; - - writeh5 = false; - if writeh5 - writeOutput(model, states, 'output.h5'); - end - -.. figure:: runBatteryP2D_01.png - :figwidth: 100% - - - -complete source code can be found :ref:`here` + +.. _runBatteryP2D: + +Pseudo-Two-Dimensional (P2D) Lithium-Ion Battery Model +-------------------------------------------------------------------------- +*Generated from runBatteryP2D.m* + + +This example demonstrates how to setup a P2D model of a Li-ion battery and run a simple simulation. + +.. code-block:: matlab + + % Clear the workspace and close open figures + clear + close all + + +Import the required modules from MRST +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +load MRST modules + +.. code-block:: matlab + + mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + + +Setup the properties of Li-ion battery materials and cell design +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The properties and parameters of the battery cell, including the architecture and materials, are set using an instance of :class:`BatteryInputParams `. This class is used to initialize the simulation and it propagates all the parameters throughout the submodels. The input parameters can be set manually or provided in json format. All the parameters for the model are stored in the inputparams object. + +.. code-block:: matlab + + jsonstruct = parseBattmoJson(fullfile('ParameterData','BatteryCellParameters','LithiumIonBatteryCell','lithium_ion_battery_nmc_graphite.json')); + + % We define some shorthand names for simplicity. + ne = 'NegativeElectrode'; + pe = 'PositiveElectrode'; + elyte = 'Electrolyte'; + thermal = 'ThermalModel'; + co = 'Coating'; + am = 'ActiveMaterial'; + itf = 'Interface'; + sd = 'SolidDiffusion'; + ctrl = 'Control'; + cc = 'CurrentCollector'; + + jsonstruct.use_thermal = false; + jsonstruct.include_current_collectors = false; + + inputparams = BatteryInputParams(jsonstruct); + + use_cccv = false; + if use_cccv + cccvstruct = struct( 'controlPolicy' , 'CCCV', ... + 'initialControl' , 'discharging', ... + 'numberOfCycles' , 1 , ... + 'CRate' , 1 , ... + 'lowerCutoffVoltage', 2.4 , ... + 'upperCutoffVoltage', 4.1 , ... + 'dIdtLimit' , 0.01 , ... + 'dEdtLimit' , 0.01); + cccvinputparams = CcCvControlModelInputParams(cccvstruct); + inputparams.Control = cccvinputparams; + end + + +Setup the geometry and computational grid +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Here, we setup the 1D computational grid that will be used for the simulation. The required discretization parameters are already included in the class BatteryGeneratorP2D. + +.. code-block:: matlab + + gen = BatteryGeneratorP2D(); + + % Now, we update the inputparams with the properties of the grid. + inputparams = gen.updateBatteryInputParams(inputparams); + + +Initialize the battery model. +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The battery model is initialized by sending inputparams to the Battery class constructor. see :class:`Battery `. + +.. code-block:: matlab + + model = Battery(inputparams); + + inspectgraph = false; + if inspectgraph + cgt = model.computationalGraph; + return + end + + +Compute the nominal cell capacity and choose a C-Rate +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The nominal capacity of the cell is calculated from the active materials. This value is then combined with the user-defined C-Rate to set the cell operational current. + +.. code-block:: matlab + + CRate = model.Control.CRate; + + +Setup the schedule +^^^^^^^^^^^^^^^^^^ + +.. code-block:: matlab + + timestep.numberOfTimeSteps = 20; + + step = model.Control.setupScheduleStep(timestep); + control = model.Control.setupScheduleControl(); + + % This control is used to set up the schedule + schedule = struct('control', control, 'step', step); + + +Setup the initial state of the model +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The initial state of the model is setup using the model.setupInitialState() method. + +.. code-block:: matlab + + initstate = model.setupInitialState(); + + +Setup the properties of the nonlinear solver +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: matlab + + nls = NonLinearSolver(); + + linearsolver = 'direct'; + switch linearsolver + case 'amgcl' + nls.LinearSolver = AMGCLSolverAD('verbose', true, 'reduceToCell', false); + nls.LinearSolver.tolerance = 1e-4; + nls.LinearSolver.maxIterations = 30; + nls.maxIterations = 10; + nls.verbose = 10; + case 'battery' + nls.LinearSolver = LinearSolverBatteryExtra('verbose' , false, ... + 'reduceToCell', true, ... + 'verbosity' , 3 , ... + 'reuse_setup' , false, ... + 'method' , 'direct'); + nls.LinearSolver.tolerance = 1e-4; + case 'direct' + disp('standard direct solver') + otherwise + error('Unknown solver %s', linearsolver); + end + + % Change default maximum iteration number in nonlinear solver + nls.maxIterations = 10; + % Change default behavior of nonlinear solver, in case of error + nls.errorOnFailure = false; + nls.timeStepSelector = StateChangeTimeStepSelector('TargetProps', {{'Control','E'}}, 'targetChangeAbs', 0.03); + % Change default tolerance for nonlinear solver + model.nonlinearTolerance = 1e-3*model.Control.Imax; + % Set verbosity + model.verbose = true; + + +Run the simulation +^^^^^^^^^^^^^^^^^^ + +.. code-block:: matlab + + [~, states, report] = simulateScheduleAD(initstate, model, schedule, 'OutputMinisteps', true, 'NonLinearSolver', nls); + + +Process output and recover the output voltage and current from the output states. +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: matlab + + ind = cellfun(@(x) not(isempty(x)), states); + states = states(ind); + E = cellfun(@(x) x.Control.E, states); + I = cellfun(@(x) x.Control.I, states); + T = cellfun(@(x) max(x.(thermal).T), states); + Tmax = cellfun(@(x) max(x.ThermalModel.T), states); + % [SOCN, SOCP] = cellfun(@(x) model.calculateSOC(x), states); + time = cellfun(@(x) x.time, states); + + figure + plot(time/hour, E); + grid on + xlabel 'time / h'; + ylabel 'potential / V'; + + writeh5 = false; + if writeh5 + writeOutput(model, states, 'output.h5'); + end + +.. figure:: runBatteryP2D_01.png + :figwidth: 100% + + + +complete source code can be found :ref:`here` diff --git a/_sources/publishedExamples/runBatteryP2D_source.rst.txt b/_sources/publishedExamples/runBatteryP2D_source.rst.txt index 98e2d250..eeef0714 100644 --- a/_sources/publishedExamples/runBatteryP2D_source.rst.txt +++ b/_sources/publishedExamples/runBatteryP2D_source.rst.txt @@ -1,203 +1,189 @@ -:orphan: - -.. _runBatteryP2D_source: - -Source code for runBatteryP2D ------------------------------ - -.. code:: matlab - - - %% Pseudo-Two-Dimensional (P2D) Lithium-Ion Battery Model - % This example demonstrates how to setup a P2D model of a Li-ion battery - % and run a simple simulation. - - % Clear the workspace and close open figures - clear - close all - - %% Import the required modules from MRST - % load MRST modules - - mrstModule add ad-core mrst-gui mpfa agmg linearsolvers - - %% Setup the properties of Li-ion battery materials and cell design - % The properties and parameters of the battery cell, including the - % architecture and materials, are set using an instance of - % :class:`BatteryInputParams `. This class is - % used to initialize the simulation and it propagates all the parameters - % throughout the submodels. The input parameters can be set manually or - % provided in json format. All the parameters for the model are stored in - % the inputparams object. - - jsonstruct = parseBattmoJson(fullfile('ParameterData','BatteryCellParameters','LithiumIonBatteryCell','lithium_ion_battery_nmc_graphite.json')); - - % We define some shorthand names for simplicity. - ne = 'NegativeElectrode'; - pe = 'PositiveElectrode'; - elyte = 'Electrolyte'; - thermal = 'ThermalModel'; - co = 'Coating'; - am = 'ActiveMaterial'; - itf = 'Interface'; - sd = 'SolidDiffusion'; - ctrl = 'Control'; - cc = 'CurrentCollector'; - - jsonstruct.use_thermal = false; - jsonstruct.include_current_collectors = false; - - inputparams = BatteryInputParams(jsonstruct); - - use_cccv = false; - if use_cccv - cccvstruct = struct( 'controlPolicy' , 'CCCV', ... - 'initialControl' , 'discharging', ... - 'CRate' , 1 , ... - 'lowerCutoffVoltage', 2.4 , ... - 'upperCutoffVoltage', 4.1 , ... - 'dIdtLimit' , 0.01 , ... - 'dEdtLimit' , 0.01); - cccvparamobj = CcCvControlModelInputParams(cccvstruct); - inputparams.Control = cccvinputparams; - end - - - %% Setup the geometry and computational grid - % Here, we setup the 1D computational grid that will be used for the - % simulation. The required discretization parameters are already included - % in the class BatteryGeneratorP2D. - gen = BatteryGeneratorP2D(); - - % Now, we update the inputparams with the properties of the grid. - inputparams = gen.updateBatteryInputParams(inputparams); - - - %% Initialize the battery model. - % The battery model is initialized by sending inputparams to the Battery class - % constructor. see :class:`Battery `. - model = Battery(inputparams); - - model.AutoDiffBackend= AutoDiffBackend(); - - inspectgraph = false; - if inspectgraph - cgt = model.computationalGraph; - return - end - - %% Compute the nominal cell capacity and choose a C-Rate - % The nominal capacity of the cell is calculated from the active materials. - % This value is then combined with the user-defined C-Rate to set the cell - % operational current. - - CRate = model.Control.CRate; - - %% Setup the time step schedule - % Smaller time steps are used to ramp up the current from zero to its - % operational value. Larger time steps are then used for the normal - % operation. - switch model.(ctrl).controlPolicy - case 'CCCV' - total = 3.5*hour/CRate; - case 'CCDischarge' - total = 1.4*hour/CRate; - otherwise - error('control policy not recognized'); - end - - n = 100; - dt = total/n; - step = struct('val', dt*ones(n, 1), 'control', ones(n, 1)); - - % we setup the control by assigning a source and stop function. - - control = model.Control.setupScheduleControl(); - - % This control is used to set up the schedule - schedule = struct('control', control, 'step', step); - - %% Setup the initial state of the model - % The initial state of the model is setup using the model.setupInitialState() method. - - initstate = model.setupInitialState(); - - %% Setup the properties of the nonlinear solver - nls = NonLinearSolver(); - - linearsolver = 'direct'; - switch linearsolver - case 'amgcl' - nls.LinearSolver = AMGCLSolverAD('verbose', true, 'reduceToCell', false); - nls.LinearSolver.tolerance = 1e-4; - nls.LinearSolver.maxIterations = 30; - nls.maxIterations = 10; - nls.verbose = 10; - case 'battery' - nls.LinearSolver = LinearSolverBatteryExtra('verbose' , false, ... - 'reduceToCell', true, ... - 'verbosity' , 3 , ... - 'reuse_setup' , false, ... - 'method' , 'direct'); - nls.LinearSolver.tolerance = 1e-4; - case 'direct' - disp('standard direct solver') - otherwise - error('Unknown solver %s', linearsolver); - end - - % Change default maximum iteration number in nonlinear solver - nls.maxIterations = 10; - % Change default behavior of nonlinear solver, in case of error - nls.errorOnFailure = false; - nls.timeStepSelector = StateChangeTimeStepSelector('TargetProps', {{'Control','E'}}, 'targetChangeAbs', 0.03); - % Change default tolerance for nonlinear solver - model.nonlinearTolerance = 1e-3*model.Control.Imax; - % Set verbosity - model.verbose = true; - - %% Run the simulation - [~, states, report] = simulateScheduleAD(initstate, model, schedule, 'OutputMinisteps', true, 'NonLinearSolver', nls); - - %% Process output and recover the output voltage and current from the output states. - ind = cellfun(@(x) not(isempty(x)), states); - states = states(ind); - E = cellfun(@(x) x.Control.E, states); - I = cellfun(@(x) x.Control.I, states); - T = cellfun(@(x) max(x.(thermal).T), states); - Tmax = cellfun(@(x) max(x.ThermalModel.T), states); - % [SOCN, SOCP] = cellfun(@(x) model.calculateSOC(x), states); - time = cellfun(@(x) x.time, states); - - figure - plot(time/hour, E); - grid on - xlabel 'time / h'; - ylabel 'potential / V'; - - writeh5 = false; - if writeh5 - writeOutput(model, states, 'output.h5'); - end - - - %{ - Copyright 2021-2023 SINTEF Industry, Sustainable Energy Technology - and SINTEF Digital, Mathematics & Cybernetics. - - This file is part of The Battery Modeling Toolbox BattMo - - BattMo is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - BattMo is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with BattMo. If not, see . - %} - +:orphan: + +.. _runBatteryP2D_source: + +Source code for runBatteryP2D +----------------------------- + +.. code:: matlab + + + %% Pseudo-Two-Dimensional (P2D) Lithium-Ion Battery Model + % This example demonstrates how to setup a P2D model of a Li-ion battery + % and run a simple simulation. + + % Clear the workspace and close open figures + clear + close all + + %% Import the required modules from MRST + % load MRST modules + + mrstModule add ad-core mrst-gui mpfa agmg linearsolvers + + %% Setup the properties of Li-ion battery materials and cell design + % The properties and parameters of the battery cell, including the + % architecture and materials, are set using an instance of + % :class:`BatteryInputParams `. This class is + % used to initialize the simulation and it propagates all the parameters + % throughout the submodels. The input parameters can be set manually or + % provided in json format. All the parameters for the model are stored in + % the inputparams object. + + jsonstruct = parseBattmoJson(fullfile('ParameterData','BatteryCellParameters','LithiumIonBatteryCell','lithium_ion_battery_nmc_graphite.json')); + + % We define some shorthand names for simplicity. + ne = 'NegativeElectrode'; + pe = 'PositiveElectrode'; + elyte = 'Electrolyte'; + thermal = 'ThermalModel'; + co = 'Coating'; + am = 'ActiveMaterial'; + itf = 'Interface'; + sd = 'SolidDiffusion'; + ctrl = 'Control'; + cc = 'CurrentCollector'; + + jsonstruct.use_thermal = false; + jsonstruct.include_current_collectors = false; + + inputparams = BatteryInputParams(jsonstruct); + + use_cccv = false; + if use_cccv + cccvstruct = struct( 'controlPolicy' , 'CCCV', ... + 'initialControl' , 'discharging', ... + 'numberOfCycles' , 1 , ... + 'CRate' , 1 , ... + 'lowerCutoffVoltage', 2.4 , ... + 'upperCutoffVoltage', 4.1 , ... + 'dIdtLimit' , 0.01 , ... + 'dEdtLimit' , 0.01); + cccvinputparams = CcCvControlModelInputParams(cccvstruct); + inputparams.Control = cccvinputparams; + end + + + %% Setup the geometry and computational grid + % Here, we setup the 1D computational grid that will be used for the + % simulation. The required discretization parameters are already included + % in the class BatteryGeneratorP2D. + gen = BatteryGeneratorP2D(); + + % Now, we update the inputparams with the properties of the grid. + inputparams = gen.updateBatteryInputParams(inputparams); + + + %% Initialize the battery model. + % The battery model is initialized by sending inputparams to the Battery class + % constructor. see :class:`Battery `. + model = Battery(inputparams); + + inspectgraph = false; + if inspectgraph + cgt = model.computationalGraph; + return + end + + %% Compute the nominal cell capacity and choose a C-Rate + % The nominal capacity of the cell is calculated from the active materials. + % This value is then combined with the user-defined C-Rate to set the cell + % operational current. + + CRate = model.Control.CRate; + + %% Setup the schedule + % + + timestep.numberOfTimeSteps = 20; + + step = model.Control.setupScheduleStep(timestep); + control = model.Control.setupScheduleControl(); + + % This control is used to set up the schedule + schedule = struct('control', control, 'step', step); + + %% Setup the initial state of the model + % The initial state of the model is setup using the model.setupInitialState() method. + + initstate = model.setupInitialState(); + + %% Setup the properties of the nonlinear solver + nls = NonLinearSolver(); + + linearsolver = 'direct'; + switch linearsolver + case 'amgcl' + nls.LinearSolver = AMGCLSolverAD('verbose', true, 'reduceToCell', false); + nls.LinearSolver.tolerance = 1e-4; + nls.LinearSolver.maxIterations = 30; + nls.maxIterations = 10; + nls.verbose = 10; + case 'battery' + nls.LinearSolver = LinearSolverBatteryExtra('verbose' , false, ... + 'reduceToCell', true, ... + 'verbosity' , 3 , ... + 'reuse_setup' , false, ... + 'method' , 'direct'); + nls.LinearSolver.tolerance = 1e-4; + case 'direct' + disp('standard direct solver') + otherwise + error('Unknown solver %s', linearsolver); + end + + % Change default maximum iteration number in nonlinear solver + nls.maxIterations = 10; + % Change default behavior of nonlinear solver, in case of error + nls.errorOnFailure = false; + nls.timeStepSelector = StateChangeTimeStepSelector('TargetProps', {{'Control','E'}}, 'targetChangeAbs', 0.03); + % Change default tolerance for nonlinear solver + model.nonlinearTolerance = 1e-3*model.Control.Imax; + % Set verbosity + model.verbose = true; + + %% Run the simulation + [~, states, report] = simulateScheduleAD(initstate, model, schedule, 'OutputMinisteps', true, 'NonLinearSolver', nls); + + %% Process output and recover the output voltage and current from the output states. + ind = cellfun(@(x) not(isempty(x)), states); + states = states(ind); + E = cellfun(@(x) x.Control.E, states); + I = cellfun(@(x) x.Control.I, states); + T = cellfun(@(x) max(x.(thermal).T), states); + Tmax = cellfun(@(x) max(x.ThermalModel.T), states); + % [SOCN, SOCP] = cellfun(@(x) model.calculateSOC(x), states); + time = cellfun(@(x) x.time, states); + + figure + plot(time/hour, E); + grid on + xlabel 'time / h'; + ylabel 'potential / V'; + + writeh5 = false; + if writeh5 + writeOutput(model, states, 'output.h5'); + end + + + %{ + Copyright 2021-2024 SINTEF Industry, Sustainable Energy Technology + and SINTEF Digital, Mathematics & Cybernetics. + + This file is part of The Battery Modeling Toolbox BattMo + + BattMo is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + BattMo is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with BattMo. If not, see . + %} + diff --git a/_sources/publishedExamples/runElectrolyser.rst.txt b/_sources/publishedExamples/runElectrolyser.rst.txt index 0e0a554b..454a78f4 100644 --- a/_sources/publishedExamples/runElectrolyser.rst.txt +++ b/_sources/publishedExamples/runElectrolyser.rst.txt @@ -1,210 +1,212 @@ - -.. _runElectrolyser: - -Alkaline Membrane Electrolyser ----------------------------------------------------- -*Generated from runElectrolyser.m* - -.. include:: runElectrolyserPreamble.rst - - -Load MRST modules -^^^^^^^^^^^^^^^^^ - -.. code-block:: matlab - - mrstModule add ad-core matlab_bgl - - -Setup input -^^^^^^^^^^^ -Setup the physical properties for the electrolyser using json input file :battmofile:`alkalineElectrolyser.json` - -.. code-block:: matlab - - jsonstruct= parseBattmoJson('Electrolyser/Parameters/alkalineElectrolyser.json'); - inputparams = ElectrolyserInputParams(jsonstruct); - -Setup the grids. We consider a 1D model and the specifications can be read from a json input file, here :battmofile:`electrolysergeometry1d.json`, using :battmo:`setupElectrolyserGridFromJson`. - -.. code-block:: matlab - - jsonstruct= parseBattmoJson('Electrolyser/Parameters/electrolysergeometry1d.json'); - inputparams = setupElectrolyserGridFromJson(inputparams, jsonstruct); - -We define shortcuts for the different submodels. - -.. code-block:: matlab - - inm = 'IonomerMembrane'; - her = 'HydrogenEvolutionElectrode'; - oer = 'OxygenEvolutionElectrode'; - ptl = 'PorousTransportLayer'; - exr = 'ExchangeReaction'; - ctl = 'CatalystLayer'; - - -Setup model -^^^^^^^^^^^ - -.. code-block:: matlab - - model = Electrolyser(inputparams); - - -Setup the initial condition -^^^^^^^^^^^^^^^^^^^^^^^^^^^ -We use the default initial setup implemented in the model - -.. code-block:: matlab - - [model, initstate] = model.setupBcAndInitialState(); - - -Setup the schedule with the time discretization -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -We run the simulation over 10 hours, increasing the input current linearly in time. - -.. code-block:: matlab - - total = 10*hour; - - n = 100; - dt = total/n; - dts = rampupTimesteps(total, dt, 5); - -We use the function :battmo:`rampupControl` to increase the current linearly in time - -.. code-block:: matlab - - controlI = -3*ampere/(centi*meter)^2; % if negative, O2 and H2 are produced - tup = total; - srcfunc = @(time) rampupControl(time, tup, controlI, 'rampupcase', 'linear'); - control = struct('src', srcfunc); - - step = struct('val', dts, 'control', ones(numel(dts), 1)); - schedule = struct('control', control, 'step', step); - - -Setup the non-linear solver -^^^^^^^^^^^^^^^^^^^^^^^^^^^ -We do only minor modifications here from the standard solver - -.. code-block:: matlab - - nls = NonLinearSolver(); - nls.verbose = false; - nls.errorOnFailure = false; - - model.verbose = false; - - -Run the simulation -^^^^^^^^^^^^^^^^^^ - -.. code-block:: matlab - - [~, states, report] = simulateScheduleAD(initstate, model, schedule, 'NonLinearSolver', nls, 'OutputMiniSteps', true); - - -Visualize the results -^^^^^^^^^^^^^^^^^^^^^ -The results contain only the primary variables of the system (the unknwons that descrive the state of the system). We use the method :code:`addVariables` to add all the intermediate quantities that are computed to solve the equations but not stored automatically in the result. - -.. code-block:: matlab - - for istate = 1 : numel(states) - states{istate} = model.addVariables(states{istate}); - end - -We extract the time, voltage and current values for each time step - -.. code-block:: matlab - - time = cellfun(@(state) state.time, states); - E = cellfun(@(state) state.(oer).(ptl).E, states); - I = cellfun(@(state) state.(oer).(ctl).I, states); - -We plot the results for the voltage and current - -.. code-block:: matlab - - set(0, 'defaultlinelinewidth', 3) - set(0, 'defaultaxesfontsize', 15) - - figure - subplot(2, 1, 1) - plot(time/hour, E) - xlabel('time [hour]'); - ylabel('voltage'); - title('Polarisation curve'); - - subplot(2, 1, 2) - plot(time/hour, -I/(1/(centi*meter)^2)); - xlabel('time [hour]'); - ylabel('Current [A/cm^2]'); - title('Input current') - -.. figure:: runElectrolyser_01.png - :figwidth: 100% - - -pH distribution plot -^^^^^^^^^^^^^^^^^^^^ -We consider the three domains and plot the pH in each of those. We setup the helper structures to iterate over each domain for the plot. - -.. code-block:: matlab - - models = {model.(oer).(ptl), ... - model.(her).(ptl), ... - model.(inm)}; - - fields = {{'OxygenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2} , ... - {'HydrogenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2}, ... - {'IonomerMembrane', 'cOH'}}; - - h = figure(); - set(h, 'position', [10, 10, 800, 450]); - hold on - - ntime = numel(time); - times = linspace(1, ntime, 10); - cmap = cmocean('deep', 10); - - for ifield = 1 : numel(fields) - - fd = fields{ifield}; - submodel = models{ifield}; - - x = submodel.grid.cells.centroids; - - for itimes = 1 : numel(times); - - itime = floor(times(itimes)); - % The method :code:`getProp` is used to recover the value from the state structure - val = model.getProp(states{itime}, fd); - pH = 14 + log10(val/(mol/litre)); - - % plot of pH for the current submodel. - plot(x/(milli*meter), pH, 'color', cmap(itimes, :)); - - end - - end - - xlabel('x / mm'); - ylabel('pH'); - title('pH distribition in electrolyser') - - colormap(cmap) - hColorbar = colorbar; - caxis([0 3]); - hTitle = get(hColorbar, 'Title'); - set(hTitle, 'string', 'J (A/cm^2)'); - -.. figure:: runElectrolyser_02.png - :figwidth: 100% - - - -complete source code can be found :ref:`here` + +.. _runElectrolyser: + +Alkaline Membrane Electrolyser +---------------------------------------------------- +*Generated from runElectrolyser.m* + +.. include:: runElectrolyserPreamble.rst + + +Load MRST modules +^^^^^^^^^^^^^^^^^ + +.. code-block:: matlab + + mrstModule add ad-core + + +Setup input +^^^^^^^^^^^ +Setup the physical properties for the electrolyser using json input file :battmofile:`alkalineElectrolyser.json` + +.. code-block:: matlab + + jsonstruct= parseBattmoJson('Electrolyser/Parameters/alkalineElectrolyser.json'); + inputparams = ElectrolyserInputParams(jsonstruct); + +Setup the grids. We consider a 1D model and the specifications can be read from a json input file, here :battmofile:`electrolysergeometry1d.json`, using :battmo:`setupElectrolyserGridFromJson`. + +.. code-block:: matlab + + jsonstruct= parseBattmoJson('Electrolyser/Parameters/electrolysergeometry1d.json'); + inputparams = setupElectrolyserGridFromJson(inputparams, jsonstruct); + +We define shortcuts for the different submodels. + +.. code-block:: matlab + + inm = 'IonomerMembrane'; + her = 'HydrogenEvolutionElectrode'; + oer = 'OxygenEvolutionElectrode'; + ptl = 'PorousTransportLayer'; + exr = 'ExchangeReaction'; + ctl = 'CatalystLayer'; + + +Setup model +^^^^^^^^^^^ + +.. code-block:: matlab + + model = Electrolyser(inputparams); + + +Setup the initial condition +^^^^^^^^^^^^^^^^^^^^^^^^^^^ +We use the default initial setup implemented in the model + +.. code-block:: matlab + + [model, initstate] = model.setupBcAndInitialState(); + + +Setup the schedule with the time discretization +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +We run the simulation over 10 hours, increasing the input current linearly in time. + +.. code-block:: matlab + + total = 10*hour; + + n = 100; + dt = total/n; + dts = rampupTimesteps(total, dt, 5); + +We use the function :battmo:`rampupControl` to increase the current linearly in time + +.. code-block:: matlab + + controlI = -3*ampere/(centi*meter)^2; % if negative, O2 and H2 are produced + tup = total; + srcfunc = @(time) rampupControl(time, tup, controlI, 'rampupcase', 'linear'); + control = struct('src', srcfunc); + + step = struct('val', dts, 'control', ones(numel(dts), 1)); + schedule = struct('control', control, 'step', step); + + +Setup the non-linear solver +^^^^^^^^^^^^^^^^^^^^^^^^^^^ +We do only minor modifications here from the standard solver + +.. code-block:: matlab + + nls = NonLinearSolver(); + nls.verbose = false; + nls.errorOnFailure = false; + + model.verbose = false; + + +Run the simulation +^^^^^^^^^^^^^^^^^^ + +.. code-block:: matlab + + [~, states, report] = simulateScheduleAD(initstate, model, schedule, 'NonLinearSolver', nls, 'OutputMiniSteps', true); + + +Visualize the results +^^^^^^^^^^^^^^^^^^^^^ +The results contain only the primary variables of the system (the unknwons that descrive the state of the system). We use the method :code:`addVariables` to add all the intermediate quantities that are computed to solve the equations but not stored automatically in the result. + +.. code-block:: matlab + + for istate = 1 : numel(states) + states{istate} = model.addVariables(states{istate}); + end + +We extract the time, voltage and current values for each time step + +.. code-block:: matlab + + time = cellfun(@(state) state.time, states); + E = cellfun(@(state) state.(oer).(ptl).E, states); + I = cellfun(@(state) state.(oer).(ctl).I, states); + +We plot the results for the voltage and current + +.. code-block:: matlab + + set(0, 'defaultlinelinewidth', 3) + set(0, 'defaultaxesfontsize', 15) + + figure + subplot(2, 1, 1) + plot(time/hour, E) + xlabel('time [hour]'); + ylabel('voltage'); + title('Polarisation curve'); + + subplot(2, 1, 2) + plot(time/hour, -I/(1/(centi*meter)^2)); + xlabel('time [hour]'); + ylabel('Current [A/cm^2]'); + title('Input current') + +.. figure:: runElectrolyser_01.png + :figwidth: 100% + + +pH distribution plot +^^^^^^^^^^^^^^^^^^^^ +We consider the three domains and plot the pH in each of those. We setup the helper structures to iterate over each domain for the plot. + +.. code-block:: matlab + + doplot = false; + + if doplot + + models = {model.(oer).(ptl), ... + model.(her).(ptl), ... + model.(inm)}; + + fields = {{'OxygenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2} , ... + {'HydrogenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2}, ... + {'IonomerMembrane', 'cOH'}}; + + h = figure(); + set(h, 'position', [10, 10, 800, 450]); + hold on + + ntime = numel(time); + times = linspace(1, ntime, 10); + cmap = cmocean('deep', 10); + + for ifield = 1 : numel(fields) + + fd = fields{ifield}; + submodel = models{ifield}; + + x = submodel.G.cells.centroids; + + for itimes = 1 : numel(times); + + itime = floor(times(itimes)); + % The method :code:`getProp` is used to recover the value from the state structure + val = model.getProp(states{itime}, fd); + pH = 14 + log10(val/(mol/litre)); + + % plot of pH for the current submodel. + plot(x/(milli*meter), pH, 'color', cmap(itimes, :)); + + end + + end + + xlabel('x / mm'); + ylabel('pH'); + title('pH distribition in electrolyser') + + colormap(cmap) + hColorbar = colorbar; + caxis([0 3]); + hTitle = get(hColorbar, 'Title'); + set(hTitle, 'string', 'J (A/cm^2)'); + end + + + +complete source code can be found :ref:`here` diff --git a/_sources/publishedExamples/runElectrolyser_source.rst.txt b/_sources/publishedExamples/runElectrolyser_source.rst.txt index eea8e505..00460a2e 100644 --- a/_sources/publishedExamples/runElectrolyser_source.rst.txt +++ b/_sources/publishedExamples/runElectrolyser_source.rst.txt @@ -1,190 +1,194 @@ -:orphan: - -.. _runElectrolyser_source: - -Source code for runElectrolyser -------------------------------- - -.. code:: matlab - - - %% Alkaline Membrane Electrolyser - - %% Load MRST modules - % - - mrstModule add ad-core matlab_bgl - - %% Setup input - % Setup the physical properties for the electrolyser using json input file :battmofile:`alkalineElectrolyser.json` - - jsonstruct= parseBattmoJson('Electrolyser/Parameters/alkalineElectrolyser.json'); - inputparams = ElectrolyserInputParams(jsonstruct); - - %% - % Setup the grids. We consider a 1D model and the specifications can be read from a json input file, here :battmofile:`electrolysergeometry1d.json`, using - % :battmo:`setupElectrolyserGridFromJson`. - - jsonstruct= parseBattmoJson('Electrolyser/Parameters/electrolysergeometry1d.json'); - inputparams = setupElectrolyserGridFromJson(inputparams, jsonstruct); - - %% - % We define shortcuts for the different submodels. - inm = 'IonomerMembrane'; - her = 'HydrogenEvolutionElectrode'; - oer = 'OxygenEvolutionElectrode'; - ptl = 'PorousTransportLayer'; - exr = 'ExchangeReaction'; - ctl = 'CatalystLayer'; - - %% Setup model - - model = Electrolyser(inputparams); - - %% Setup the initial condition - % We use the default initial setup implemented in the model - - [model, initstate] = model.setupBcAndInitialState(); - - %% Setup the schedule with the time discretization - % We run the simulation over 10 hours, increasing the input current linearly in time. - - total = 10*hour; - - n = 100; - dt = total/n; - dts = rampupTimesteps(total, dt, 5); - - %% - % We use the function :code:`rampupControl` to increase the current linearly in time - - controlI = -3*ampere/(centi*meter)^2; % if negative, O2 and H2 are produced - tup = total; - srcfunc = @(time) rampupControl(time, tup, controlI, 'rampupcase', 'linear'); - control = struct('src', srcfunc); - - step = struct('val', dts, 'control', ones(numel(dts), 1)); - schedule = struct('control', control, 'step', step); - - %% Setup the non-linear solver - % We do only minor modifications here from the standard solver - - nls = NonLinearSolver(); - nls.verbose = false; - nls.errorOnFailure = false; - - model.verbose = false; - - %% Run the simulation - - [~, states, report] = simulateScheduleAD(initstate, model, schedule, 'NonLinearSolver', nls, 'OutputMiniSteps', true); - - %% Visualize the results - % - % The results contain only the primary variables of the system (the unknwons that descrive the state of the system). We - % use the method :code:`addVariables` to add all the intermediate quantities that are computed to solve the equations - % but not stored automatically in the result. - - for istate = 1 : numel(states) - states{istate} = model.addVariables(states{istate}); - end - - %% - % We extract the time, voltage and current values for each time step - - time = cellfun(@(state) state.time, states); - E = cellfun(@(state) state.(oer).(ptl).E, states); - I = cellfun(@(state) state.(oer).(ctl).I, states); - - %% - % We plot the results for the voltage and current - - set(0, 'defaultlinelinewidth', 3) - set(0, 'defaultaxesfontsize', 15) - - figure - subplot(2, 1, 1) - plot(time/hour, E) - xlabel('time [hour]'); - ylabel('voltage'); - title('Polarisation curve'); - - subplot(2, 1, 2) - plot(time/hour, -I/(1/(centi*meter)^2)); - xlabel('time [hour]'); - ylabel('Current [A/cm^2]'); - title('Input current') - - %% pH distribution plot - % - % We consider the three domains and plot the pH in each of those. We setup the helper structures to iterate over each - % domain for the plot. - - models = {model.(oer).(ptl), ... - model.(her).(ptl), ... - model.(inm)}; - - fields = {{'OxygenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2} , ... - {'HydrogenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2}, ... - {'IonomerMembrane', 'cOH'}}; - - h = figure(); - set(h, 'position', [10, 10, 800, 450]); - hold on - - ntime = numel(time); - times = linspace(1, ntime, 10); - cmap = cmocean('deep', 10); - - for ifield = 1 : numel(fields) - - fd = fields{ifield}; - submodel = models{ifield}; - - x = submodel.grid.cells.centroids; - - for itimes = 1 : numel(times); - - itime = floor(times(itimes)); - % The method :code:`getProp` is used to recover the value from the state structure - val = model.getProp(states{itime}, fd); - pH = 14 + log10(val/(mol/litre)); - - % plot of pH for the current submodel. - plot(x/(milli*meter), pH, 'color', cmap(itimes, :)); - - end - - end - - xlabel('x / mm'); - ylabel('pH'); - title('pH distribition in electrolyser') - - colormap(cmap) - hColorbar = colorbar; - caxis([0 3]); - hTitle = get(hColorbar, 'Title'); - set(hTitle, 'string', 'J (A/cm^2)'); - - - %{ - Copyright 2021-2023 SINTEF Industry, Sustainable Energy Technology - and SINTEF Digital, Mathematics & Cybernetics. - - This file is part of The Battery Modeling Toolbox BattMo - - BattMo is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - BattMo is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with BattMo. If not, see . - %} - +:orphan: + +.. _runElectrolyser_source: + +Source code for runElectrolyser +------------------------------- + +.. code:: matlab + + + %% Alkaline Membrane Electrolyser + + %% Load MRST modules + % + + mrstModule add ad-core + + %% Setup input + % Setup the physical properties for the electrolyser using json input file :battmofile:`alkalineElectrolyser.json` + + jsonstruct= parseBattmoJson('Electrolyser/Parameters/alkalineElectrolyser.json'); + inputparams = ElectrolyserInputParams(jsonstruct); + + %% + % Setup the grids. We consider a 1D model and the specifications can be read from a json input file, here :battmofile:`electrolysergeometry1d.json`, using + % :battmo:`setupElectrolyserGridFromJson`. + + jsonstruct= parseBattmoJson('Electrolyser/Parameters/electrolysergeometry1d.json'); + inputparams = setupElectrolyserGridFromJson(inputparams, jsonstruct); + + %% + % We define shortcuts for the different submodels. + inm = 'IonomerMembrane'; + her = 'HydrogenEvolutionElectrode'; + oer = 'OxygenEvolutionElectrode'; + ptl = 'PorousTransportLayer'; + exr = 'ExchangeReaction'; + ctl = 'CatalystLayer'; + + %% Setup model + + model = Electrolyser(inputparams); + + %% Setup the initial condition + % We use the default initial setup implemented in the model + + [model, initstate] = model.setupBcAndInitialState(); + + %% Setup the schedule with the time discretization + % We run the simulation over 10 hours, increasing the input current linearly in time. + + total = 10*hour; + + n = 100; + dt = total/n; + dts = rampupTimesteps(total, dt, 5); + + %% + % We use the function :battmo:`rampupControl` to increase the current linearly in time + + controlI = -3*ampere/(centi*meter)^2; % if negative, O2 and H2 are produced + tup = total; + srcfunc = @(time) rampupControl(time, tup, controlI, 'rampupcase', 'linear'); + control = struct('src', srcfunc); + + step = struct('val', dts, 'control', ones(numel(dts), 1)); + schedule = struct('control', control, 'step', step); + + %% Setup the non-linear solver + % We do only minor modifications here from the standard solver + + nls = NonLinearSolver(); + nls.verbose = false; + nls.errorOnFailure = false; + + model.verbose = false; + + %% Run the simulation + + [~, states, report] = simulateScheduleAD(initstate, model, schedule, 'NonLinearSolver', nls, 'OutputMiniSteps', true); + + %% Visualize the results + % + % The results contain only the primary variables of the system (the unknwons that descrive the state of the system). We + % use the method :code:`addVariables` to add all the intermediate quantities that are computed to solve the equations + % but not stored automatically in the result. + + for istate = 1 : numel(states) + states{istate} = model.addVariables(states{istate}); + end + + %% + % We extract the time, voltage and current values for each time step + + time = cellfun(@(state) state.time, states); + E = cellfun(@(state) state.(oer).(ptl).E, states); + I = cellfun(@(state) state.(oer).(ctl).I, states); + + %% + % We plot the results for the voltage and current + + set(0, 'defaultlinelinewidth', 3) + set(0, 'defaultaxesfontsize', 15) + + figure + subplot(2, 1, 1) + plot(time/hour, E) + xlabel('time [hour]'); + ylabel('voltage'); + title('Polarisation curve'); + + subplot(2, 1, 2) + plot(time/hour, -I/(1/(centi*meter)^2)); + xlabel('time [hour]'); + ylabel('Current [A/cm^2]'); + title('Input current') + + %% pH distribution plot + % + % We consider the three domains and plot the pH in each of those. We setup the helper structures to iterate over each + % domain for the plot. + + doplot = false; + + if doplot + + models = {model.(oer).(ptl), ... + model.(her).(ptl), ... + model.(inm)}; + + fields = {{'OxygenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2} , ... + {'HydrogenEvolutionElectrode', 'PorousTransportLayer', 'concentrations', 2}, ... + {'IonomerMembrane', 'cOH'}}; + + h = figure(); + set(h, 'position', [10, 10, 800, 450]); + hold on + + ntime = numel(time); + times = linspace(1, ntime, 10); + cmap = cmocean('deep', 10); + + for ifield = 1 : numel(fields) + + fd = fields{ifield}; + submodel = models{ifield}; + + x = submodel.G.cells.centroids; + + for itimes = 1 : numel(times); + + itime = floor(times(itimes)); + % The method :code:`getProp` is used to recover the value from the state structure + val = model.getProp(states{itime}, fd); + pH = 14 + log10(val/(mol/litre)); + + % plot of pH for the current submodel. + plot(x/(milli*meter), pH, 'color', cmap(itimes, :)); + + end + + end + + xlabel('x / mm'); + ylabel('pH'); + title('pH distribition in electrolyser') + + colormap(cmap) + hColorbar = colorbar; + caxis([0 3]); + hTitle = get(hColorbar, 'Title'); + set(hTitle, 'string', 'J (A/cm^2)'); + end + + %{ + Copyright 2021-2024 SINTEF Industry, Sustainable Energy Technology + and SINTEF Digital, Mathematics & Cybernetics. + + This file is part of The Battery Modeling Toolbox BattMo + + BattMo is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + BattMo is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with BattMo. If not, see . + %} + diff --git a/_sources/publishedExamples/runSEIActiveMaterial.rst.txt b/_sources/publishedExamples/runSEIActiveMaterial.rst.txt index 84ba86d0..fc1aa6b0 100644 --- a/_sources/publishedExamples/runSEIActiveMaterial.rst.txt +++ b/_sources/publishedExamples/runSEIActiveMaterial.rst.txt @@ -52,10 +52,11 @@ Setup the model .. code-block:: matlab % We use a stand alone model for the particle - inputparams.standAlone = true; + inputparams.isRootSimulationModel = true; % We initiate the model model = SEIActiveMaterial(inputparams); + model = model.setupForSimulation(); Setup initial state diff --git a/_sources/publishedExamples/runSEIActiveMaterial_source.rst.txt b/_sources/publishedExamples/runSEIActiveMaterial_source.rst.txt index 1ac3e82a..c989b6aa 100644 --- a/_sources/publishedExamples/runSEIActiveMaterial_source.rst.txt +++ b/_sources/publishedExamples/runSEIActiveMaterial_source.rst.txt @@ -1,218 +1,219 @@ -:orphan: - -.. _runSEIActiveMaterial_source: - -Source code for runSEIActiveMaterial ------------------------------------- - -.. code:: matlab - - - %% Particle simulation with SEI layer growth - - % clear the workspace and close open figures - clear - close all - - %% Import the required modules from MRST - % load MRST modules - mrstModule add ad-core mrst-gui mpfa - - %% Setup the properties of Li-ion battery materials and cell design - jsonstruct = parseBattmoJson(fullfile('ParameterData','ParameterSets','Safari2009','anode_sei.json')); - - % Some shorthands used for the sub-models - ne = 'NegativeElectrode'; - pe = 'PositiveElectrode'; - am = 'ActiveMaterial'; - sd = 'SolidDiffusion'; - itf = 'Interface'; - sei = 'SolidElectrodeInterface'; - sr = 'SideReaction'; - elyte = 'Electrolyte'; - - inputparams = SEIActiveMaterialInputParams(jsonstruct); - - inputparams.(sd).N = 10; - inputparams.(sei).N = 10; - - %% Setup the model - - % We use a stand alone model for the particle - inputparams.standAlone = true; - - % We initiate the model - model = SEIActiveMaterial(inputparams); - - %% Setup initial state - - Nsd = model.(sd).N; - Nsei = model.(sei).N; - - % Initial concentration value at the electrode - cElectrodeInit = 0.75*model.(itf).saturationConcentration; - % Initial value of the potential at the electrode - phiElectrodeInit = 0; - % Initial concentration value in the electrolyte - cElectrolyte = 5e-1*mol/litre; - % Temperature - T = 298.15; - - % The following datas come from :cite:`Safari_2009` - % Porosity of the SEI film - epsiSEI = 0.05; - % Solvent concentration in the bulk of the electrolyte - cECsolution = 4.541*mol/litre; - % Solvent concentration in the SEI film - cECexternal = epsiSEI*cECsolution; - - % We compute the OCP from the given data and use it to assign electrical potential in electrolyte - initState.T = T; - initState.(sd).cSurface = cElectrodeInit; - initState = model.evalVarName(initState, {itf, 'OCP'}); - - OCP = initState.(itf).OCP; - phiElectrolyte = phiElectrodeInit - OCP; - - % From the values computed above we set the values of the initial state - initState.E = phiElectrodeInit; - initState.I = 0; - initState.(sd).c = cElectrodeInit*ones(Nsd, 1); - initState.(sei).c = cECexternal*ones(Nsei, 1); - initState.(sei).cInterface = cECexternal; - initState.(sei).delta = 5*nano*meter; - initState.R = 0; - - % we set also static variable fields - initState.T = T; - initState.(itf).cElectrolyte = cElectrolyte; - initState.(itf).phiElectrolyte = phiElectrolyte; - initState.(sr).phiElectrolyte = phiElectrolyte; - initState.(sei).cExternal = cECexternal; - - - %% Setup schedule - - % Reference rate which roughly corresponds to 1 hour for the data of this example - Iref = 1.3e-4*ampere/(1*centi*meter)^2; - - Imax = 1e1*Iref; - - total = 1*hour*(Iref/Imax); - n = 100; - dt = total/n; - step = struct('val', dt*ones(n, 1), 'control', ones(n, 1)); - - % rampup value for the current function, see rampupSwitchControl - tup = dt; - srcfunc = @(time) rampupControl(time, tup, Imax); - - cmin = (model.(itf).guestStoichiometry0)*(model.(itf).saturationConcentration); - control.stopFunction = @(model, state, state0_inner) (state.(sd).cSurface <= cmin); - control.src = srcfunc; - - schedule = struct('control', control, 'step', step); - - %% Setup non-linear solver - - nls = NonLinearSolver(); - nls.errorOnFailure = false; - - model.nonlinearTolerance = 1e-5; - - %% Run simulation - - model.verbose = true; - [~, states, report] = simulateScheduleAD(initState, model, schedule, 'OutputMinisteps', true, 'NonLinearSolver', nls); - - %% Plotting - - set(0, 'defaulttextfontsize', 15); - set(0, 'defaultaxesfontsize', 15); - set(0, 'defaultlinelinewidth', 3); - set(0, 'defaultfigureposition', [10, 10, 800, 400]); - - ind = cellfun(@(state) ~isempty(state), states); - states = states(ind); - - time = cellfun(@(state) state.time, states); - - cSurface = cellfun(@(state) state.(sd).cSurface, states); - figure - plot(time/hour, cSurface/(1/litre)); - xlabel('time / h'); - ylabel('Surface concentration / mol/L'); - title('Surface concentration'); - - E = cellfun(@(state) state.E, states); - figure - plot(time/hour, E); - xlabel('time / h'); - ylabel('Potential / V'); - title('Potential'); - - - cmin = cellfun(@(state) min(state.(sd).c), states); - cmax = cellfun(@(state) max(state.(sd).c), states); - - for istate = 1 : numel(states) - states{istate} = model.evalVarName(states{istate}, {sd, 'cAverage'}); - end - - caver = cellfun(@(state) max(state.(sd).cAverage), states); - - figure - hold on - plot(time/hour, cmin /(mol/litre), 'displayname', 'cmin'); - plot(time/hour, cmax /(mol/litre), 'displayname', 'cmax'); - plot(time/hour, caver/(mol/litre), 'displayname', 'total concentration'); - title('Concentration in particle / mol/L') - legend show - - delta = cellfun(@(state) state.(sei).delta, states); - figure - plot(time/hour, delta/(nano*meter)); - xlabel('time [hour]'); - ylabel('thickness / nm'); - title('SEI thickness') - - c = states{end}.(sd).c; - r = linspace(0, model.(sd).particleRadius, model.(sd).N); - - figure - plot(r, c/(mol/litre)); - xlabel('radius / m') - ylabel('concentration / mol/L') - title('Particle concentration profile (last time step)') - - r = states{end}.(sei).delta; - r = linspace(0, r, model.(sei).N); - c = states{end}.(sei).c; - - figure - plot(r/(nano*meter), c/(mol/litre)); - xlabel('x / mm') - ylabel('concentration / mol/L'); - title('Concentration profile in SEI layer (last time step)'); - - - %{ - Copyright 2021-2023 SINTEF Industry, Sustainable Energy Technology - and SINTEF Digital, Mathematics & Cybernetics. - - This file is part of The Battery Modeling Toolbox BattMo - - BattMo is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - BattMo is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with BattMo. If not, see . - %} - +:orphan: + +.. _runSEIActiveMaterial_source: + +Source code for runSEIActiveMaterial +------------------------------------ + +.. code:: matlab + + + %% Particle simulation with SEI layer growth + + % clear the workspace and close open figures + clear + close all + + %% Import the required modules from MRST + % load MRST modules + mrstModule add ad-core mrst-gui mpfa + + %% Setup the properties of Li-ion battery materials and cell design + jsonstruct = parseBattmoJson(fullfile('ParameterData','ParameterSets','Safari2009','anode_sei.json')); + + % Some shorthands used for the sub-models + ne = 'NegativeElectrode'; + pe = 'PositiveElectrode'; + am = 'ActiveMaterial'; + sd = 'SolidDiffusion'; + itf = 'Interface'; + sei = 'SolidElectrodeInterface'; + sr = 'SideReaction'; + elyte = 'Electrolyte'; + + inputparams = SEIActiveMaterialInputParams(jsonstruct); + + inputparams.(sd).N = 10; + inputparams.(sei).N = 10; + + %% Setup the model + + % We use a stand alone model for the particle + inputparams.isRootSimulationModel = true; + + % We initiate the model + model = SEIActiveMaterial(inputparams); + model = model.setupForSimulation(); + + %% Setup initial state + + Nsd = model.(sd).N; + Nsei = model.(sei).N; + + % Initial concentration value at the electrode + cElectrodeInit = 0.75*model.(itf).saturationConcentration; + % Initial value of the potential at the electrode + phiElectrodeInit = 0; + % Initial concentration value in the electrolyte + cElectrolyte = 5e-1*mol/litre; + % Temperature + T = 298.15; + + % The following datas come from :cite:`Safari_2009` + % Porosity of the SEI film + epsiSEI = 0.05; + % Solvent concentration in the bulk of the electrolyte + cECsolution = 4.541*mol/litre; + % Solvent concentration in the SEI film + cECexternal = epsiSEI*cECsolution; + + % We compute the OCP from the given data and use it to assign electrical potential in electrolyte + initState.T = T; + initState.(sd).cSurface = cElectrodeInit; + initState = model.evalVarName(initState, {itf, 'OCP'}); + + OCP = initState.(itf).OCP; + phiElectrolyte = phiElectrodeInit - OCP; + + % From the values computed above we set the values of the initial state + initState.E = phiElectrodeInit; + initState.I = 0; + initState.(sd).c = cElectrodeInit*ones(Nsd, 1); + initState.(sei).c = cECexternal*ones(Nsei, 1); + initState.(sei).cInterface = cECexternal; + initState.(sei).delta = 5*nano*meter; + initState.R = 0; + + % we set also static variable fields + initState.T = T; + initState.(itf).cElectrolyte = cElectrolyte; + initState.(itf).phiElectrolyte = phiElectrolyte; + initState.(sr).phiElectrolyte = phiElectrolyte; + initState.(sei).cExternal = cECexternal; + + + %% Setup schedule + + % Reference rate which roughly corresponds to 1 hour for the data of this example + Iref = 1.3e-4*ampere/(1*centi*meter)^2; + + Imax = 1e1*Iref; + + total = 1*hour*(Iref/Imax); + n = 100; + dt = total/n; + step = struct('val', dt*ones(n, 1), 'control', ones(n, 1)); + + % rampup value for the current function, see rampupSwitchControl + tup = dt; + srcfunc = @(time) rampupControl(time, tup, Imax); + + cmin = (model.(itf).guestStoichiometry0)*(model.(itf).saturationConcentration); + control.stopFunction = @(model, state, state0_inner) (state.(sd).cSurface <= cmin); + control.src = srcfunc; + + schedule = struct('control', control, 'step', step); + + %% Setup non-linear solver + + nls = NonLinearSolver(); + nls.errorOnFailure = false; + + model.nonlinearTolerance = 1e-5; + + %% Run simulation + + model.verbose = true; + [~, states, report] = simulateScheduleAD(initState, model, schedule, 'OutputMinisteps', true, 'NonLinearSolver', nls); + + %% Plotting + + set(0, 'defaulttextfontsize', 15); + set(0, 'defaultaxesfontsize', 15); + set(0, 'defaultlinelinewidth', 3); + set(0, 'defaultfigureposition', [10, 10, 800, 400]); + + ind = cellfun(@(state) ~isempty(state), states); + states = states(ind); + + time = cellfun(@(state) state.time, states); + + cSurface = cellfun(@(state) state.(sd).cSurface, states); + figure + plot(time/hour, cSurface/(1/litre)); + xlabel('time / h'); + ylabel('Surface concentration / mol/L'); + title('Surface concentration'); + + E = cellfun(@(state) state.E, states); + figure + plot(time/hour, E); + xlabel('time / h'); + ylabel('Potential / V'); + title('Potential'); + + + cmin = cellfun(@(state) min(state.(sd).c), states); + cmax = cellfun(@(state) max(state.(sd).c), states); + + for istate = 1 : numel(states) + states{istate} = model.evalVarName(states{istate}, {sd, 'cAverage'}); + end + + caver = cellfun(@(state) max(state.(sd).cAverage), states); + + figure + hold on + plot(time/hour, cmin /(mol/litre), 'displayname', 'cmin'); + plot(time/hour, cmax /(mol/litre), 'displayname', 'cmax'); + plot(time/hour, caver/(mol/litre), 'displayname', 'total concentration'); + title('Concentration in particle / mol/L') + legend show + + delta = cellfun(@(state) state.(sei).delta, states); + figure + plot(time/hour, delta/(nano*meter)); + xlabel('time [hour]'); + ylabel('thickness / nm'); + title('SEI thickness') + + c = states{end}.(sd).c; + r = linspace(0, model.(sd).particleRadius, model.(sd).N); + + figure + plot(r, c/(mol/litre)); + xlabel('radius / m') + ylabel('concentration / mol/L') + title('Particle concentration profile (last time step)') + + r = states{end}.(sei).delta; + r = linspace(0, r, model.(sei).N); + c = states{end}.(sei).c; + + figure + plot(r/(nano*meter), c/(mol/litre)); + xlabel('x / mm') + ylabel('concentration / mol/L'); + title('Concentration profile in SEI layer (last time step)'); + + + %{ + Copyright 2021-2024 SINTEF Industry, Sustainable Energy Technology + and SINTEF Digital, Mathematics & Cybernetics. + + This file is part of The Battery Modeling Toolbox BattMo + + BattMo is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + BattMo is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with BattMo. If not, see . + %} + diff --git a/_sources/publishedExamples/runSiliconGraphiteBattery.rst.txt b/_sources/publishedExamples/runSiliconGraphiteBattery.rst.txt index bf09f74f..4653b625 100644 --- a/_sources/publishedExamples/runSiliconGraphiteBattery.rst.txt +++ b/_sources/publishedExamples/runSiliconGraphiteBattery.rst.txt @@ -42,7 +42,7 @@ We load the property of a composite silicon graphite electrode, see :ref:`compos .. code-block:: matlab - jsonstruct_composite_material = parseBattmoJson('ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_silicon_graphite.json'); + jsonstruct_composite_material = parseBattmoJson('ParameterData/BatteryCellParameters/LithiumIonBatteryCell/composite_silicon_graphite.json'); For the remaining properties, we consider a standard data set diff --git a/_sources/publishedExamples/runSiliconGraphiteBattery_source.rst.txt b/_sources/publishedExamples/runSiliconGraphiteBattery_source.rst.txt index 8849c36a..a5000e38 100644 --- a/_sources/publishedExamples/runSiliconGraphiteBattery_source.rst.txt +++ b/_sources/publishedExamples/runSiliconGraphiteBattery_source.rst.txt @@ -35,7 +35,7 @@ Source code for runSiliconGraphiteBattery % We load the property of a composite silicon graphite electrode, see :ref:`compositeElectrode` % - jsonstruct_composite_material = parseBattmoJson('ParameterData/BatteryCellParameters/LithiumIonBatteryCell/lithium_ion_battery_nmc_silicon_graphite.json'); + jsonstruct_composite_material = parseBattmoJson('ParameterData/BatteryCellParameters/LithiumIonBatteryCell/composite_silicon_graphite.json'); %% % For the remaining properties, we consider a standard data set diff --git a/_static/notebooks/part_1_battery_modeling_guide.html b/_static/notebooks/part_1_battery_modeling_guide.html new file mode 100644 index 00000000..3efbaf95 --- /dev/null +++ b/_static/notebooks/part_1_battery_modeling_guide.html @@ -0,0 +1,1022 @@ + +Untitled

Introduction to battery modelling with BATTMO

(I have added a few questions in Italics)
Rechargeable Li-ion batteries have been transforming our lives over the past two decades. With the current focus on electrifying the mobility sector and the growing need for large energy storage devices to support intermittent renewable energy sources, there is much to learn about batteries and the materials that make them. The goal is to develop batteries with higher energy density and improved safety, manufactured in a sustainable manner.
Battmo is a powerful tool for conducting battery simulations, aiding in the understanding of battery performance, lifecycle, and safety concerns. It utilizes continuum scale P2D (pseudo-two-dimensional) simulations based on the Doyle-Fuller-Newman (DFN) model. By inputting material parameters for the various components of the battery, Battmo can simulate capacity losses that occur during charge-discharge cycles. It supports simulations from simple 1D systems to complex multi-layer pouch cells and jelly rolls, incorporating thermal effects and degradation mechanisms. Its modular design facilitates easy switching between different cell chemistries, making Battmo versatile for applications in battery research and development.
For a detailed overview of the modelling approach and the DFN model, refer to part 2 of the modelling guide.

Basics

battery.png
A schematic representation of a typical Li-ion cell during charging: The cell is a stack of a graphite anode, a LiCoO2 cathode and a Li+ containing electrolyte. The external voltage drives out Li+ and electrons from the cathode into the anode. Electrons are released from the oxidation of cobalt in the cathode (Co3+ -> Co4+)[1].
Let's begin by explaining the fundamental principles and operation of a lithium-ion (Li-ion) battery. The battery cell consists of three main components: the anode, typically made of a graphitic material; the cathode, usually composed of a transition metal oxide; and an electrolyte, which is an organic solvent containing a lithium salt. We will look at the exact types of materials used for these components in below. A rechargeable Li-ion battery operates on the principles of intercalation and deintercalation of lithium ions. During charging and discharging cycles, lithium ions move between the anode and cathode through the electrolyte, via a separator, the electrons are driven through the external circuit, providing electrical energy to connected devices. The energy stored is calculated as the charge that is transported times the potential difference between the electrodes.
Should I add more to this section, for instance, the images of the oxide materials itself?
Electrodes:
Although we tend to think of electrodes as solid materials, they are actually porous, consisting of a mixture of active material, binder, and solvent, all immersed in the electrolyte, as can be seen in the schematic representation above.
The electrode materials are applied as a coating on the current collectors. This coating contains the active material of the electrode, such as graphite or metal oxides, along with other inactive materials like binders, conducting additives and solvents.
Cathode:
Standard active materials for Li-ion batteries are the layered Lithium transition metal oxides, these materials come in different compositions and crystal structures.
The most common oxides include NMC (Nickel Manganese Cobalt), LFP (Lithium Iron Phosphate), LMO (Lithium Manganese Oxide), LCO (Lithium Cobalt Oxide) and NCA (Nickel Cobalt Aluminum). Because these materials differ in energy density, cycling stability, commercial availability, and environmental impact, selecting the appropriate one for each application is crucial. For a thorough review of cathode materials, the reader can refer here.
Anode:
The anode in early rechargeable lithium batteries was pure lithium metal, but its tendency to form dendrites or lithium plating, which could lead to short circuits during charging, limited its use. Today, graphite is the standard anode material for Li-ion cells. Compared to cathode materials, graphite is more affordable and widely available. For a review on on carbon based anode materials, one can refer here.
Electrolyte:
The electrolyte serves as the medium for Li-ion conduction and can be either solid or liquid. Liquid electrolytes typically consist of lithium salts dissolved in an organic solvent, with common examples being LiPF6, LiBF4, and LiClO4. Commonly used solvents are ethylene carbonate, dimethyl carbonate, and diethyl carbonate. A good electrolyte, whether solid or liquid must possess high ionic conductivity, low electronic conductivity, support high voltages, and remain stable across a wide temperature range. For a review on liquid electrolytes, one can refer here.
Working:
As an example, Let's take NMC as the cathode material, During a typical discharge cycle,
The reaction taking place at the cathode (postive electrode):
At the anode (negative electrode):
The overall cell reaction:
The system is at its lowest energy when the Li-ions are in the cathode, during discharge, the Li ions migrate from the anode to the cathode and the potential difference is released as energy.

Modelling

Battmo comes with several default material and cell parameter sets; for more details, refer here. Users can also input their own parameterization data from experiments. It uses the JSON format to manage input parameters. This allows you to easily save, document, and share complete parameter sets from specific simulations. If you are new to JSON, you can learn more about it here.
For the rest of this guide, we will use a starter JSON file, sample_input.json that describes a Lithium nickel manganese cobalt oxide (NMC 811) - Graphite cell [2] with LiPF6 [3] as the electrolyte, in a 1D cell according to BattMo's JSON input specification.
This sample_input.json file contains all the information needed to run the simulation. These parameters are neatly arranged into structured schemas. Unless otherwise specified, BattMo uses SI base units for physical quantities.
  • A schema for the geometry, contains the dimensions of the whole cell and the different components within the cell.
  • A schema for the battery cell material parameters, contains material parameters for the NMC - Graphite cell sourced from the literature.
  • A schema for the initialization of the state of the battery, specifies the initial state of the battery, including the initial state of charge (SOC).
  • A schema for control model, defines the charging and discharging protocols,
  • A schema for the time stepping,
  • A schema for the solver parameters, specifies the parameters for the numerical solver used in the simulation.
  • A schema for the output specification
As a first step, let us look at how a battery looks in BattMo. There are convenient built-in functions that help in parsing and visualizing data.
Parse and Visualize your battery grid:
  • Load the default sample_input.json file using parseBattmoJson function.
  • Visualize the 1D system using plotBatteryGrid.
  • Merge the default file with a 3D Geometry File with mergeJsonStructs and visualise a 3D battery cell.
jsonstruct_default = parseBattmoJson('Examples/JsonDataFiles/sample_input.json');
model_1d = setupModelFromJson(jsonstruct_default);
% use the plotBatteryGrid function to show the grid
plotBatteryGrid(model_1d)
axis tight
view(45,45)
jsonstruct_3d_geometry = parseBattmoJson('Examples/JsonDataFiles/geometry3d.json');
jsonstruct_modified_3d = mergeJsonStructs({jsonstruct_3d_geometry , ...
jsonstruct_default});
parameter include_current_collectors is assigned twice with different values. we use the value from first jsonstruct. +parameter Geometry.case is assigned twice with different values. we use the value from first jsonstruct. +parameter NegativeElectrode.Coating.thickness is assigned twice with different values. we use the value from first jsonstruct. +parameter NegativeElectrode.Coating.N is assigned twice with different values. we use the value from first jsonstruct. +parameter NegativeElectrode.CurrentCollector.N is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.Coating.thickness is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.Coating.N is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.CurrentCollector.N is assigned twice with different values. we use the value from first jsonstruct. +parameter Separator.thickness is assigned twice with different values. we use the value from first jsonstruct. +parameter Separator.N is assigned twice with different values. we use the value from first jsonstruct. +parameter ThermalModel.externalHeatTransferCoefficient is assigned twice with different values. we use the value from first jsonstruct.
model_3d = setupModelFromJson(jsonstruct_modified_3d);
plotBatteryGrid(model_3d)
% make the axis tight and set the camera viewing angle
axis tight
view(45,45)

Characterising a battery

A good battery should posess high energy density, ensuring it can store a significant amount of energy relative to its size and weight. Additionally, it should maintain a long cycle life, allowing it to endure numerous charge and discharge cycles without significant degradation. Minimal energy loss during these cycles is crucial, ensuring efficient energy usage and prolonged operational efficiency.
Let's delve into essential terms and parameters crucial for understanding batteries and their simulation in Battmo. These help in evaluating and comparing battery performance and characteristics.
Definitions:
  • Lower (upper) cut off voltage:The minimum (maximum) allowable voltage. It is this voltage that generally defines the “empty” ("full") state of the battery, this is set as an input parameter in the simulation.The value depends on the cell chemistry.
  • C-rate (cycling rate): this is the rate at which a battery is discharged in relation to its maximum capacity. Specifically, a 1C rate signifies that the discharge current will exhaust the total battery capacity in 1 hour. The CRate chosen depends on your application, this is set as an input parameter in the simulation
  • State of Charge (SOC %): ratio of the current battery capacity to its maximum capacity, this is set as an input parameter in the simulation
  • OpenCircuitVoltage (OCV): The voltage between the cathode and the anode with no external load, to calculate this in Battmo,one can use a very small charge/discharge current, a CRate of C/20. OCV curves are typically plotted as voltage versus state of charge (SOC), the curves depend on the chemistry of the cell.
Characterization of a cell:
  • Capacity (Ah): electric charge which a cell or battery can deliver under specified discharge conditions, calculated as a product of discharge current and discharge time. Increasing the CRate decreases the capacity due to increasing internal losses.
  • Nominal Capacity (Ah): the rated capacity that is representative of the typical value associated with a cell as stated by the manufacturer
  • Theoretical capacity: maximum possible charge, a cell can store, it is a fixed value determined by the amount and type of active material used in the electrodes.
  • Energy(Wh):calculated as the product of capacity and the average discharge voltage when discharged at a specific C-rate,Increasing the CRate decreases the energy.
  • Energy density (Wh/l): calculated as the energy divided by the volume of the cell
  • Specific energy (Wh/kg): calculated as the energy divided by the mass of the cell
  • Efficiency: ratio of the energy released during discharge to the energy stored during charge. This parameter indicates how effectively a battery converts stored energy into usable energy.
Charging/discharging protocols:
In BattMo, we set the protocols for the charging and discharging of the cell using the control model schema. We discuss the basic protocols here,
Charging-CCCV
A Li-ion battery is typically charged using a constant current/constant voltage protocol (CC-CV). This means the battery is initially charged at a constant current until it reaches a maximum voltage (upperCutoffVoltage). Once this voltage is reached, the protocol switches to maintaining a constant voltage between the terminals while the charging current gradually decreases.
Discharging-CC
During discharge, a constant current (CC) based on the chosen C-rate is drawn from the cell until the voltage reaches the minimum cut-off voltage (lowerCutoffVoltage) of the cell.

CC discharge of NMC-Graphite cell:

Now we are ready to run our first simulation. Refer to the `sample_input.json` for the exact input values. In summary, we will simulate the constant current (CC) discharge of an NMC-Graphite system, starting from a 100% state of charge (SOC) at a C-rate of 1. Before running the actual simulation, the function `computeCellCapacity.m` processes the material parameters for the active materials in the electrodes, along with their mass fractions, to compute the maximum cell capacity (theoretical capacity). Using this capacity, the program sets the discharge current accordingly()
Can we show the theoretical capacity value, 'cell specification summary'? For this system, the typical lower and upper cutoff voltages are 2.4V and 4.1V, respectively. The simulation stops when the lower cut-off voltage is reached.
jsonstruct = parseBattmoJson('Examples/JsonDataFiles/sample_input.json');% parse the sample_input.json
disp(jsonstruct.Control);% see the default values set in the control policy
controlPolicy: 'CCDischarge' + CRate: 1 + lowerCutoffVoltage: 2.4000 + upperCutoffVoltage: 4.1000 + dIdtLimit: 0.0100 + dEdtLimit: 0.0100 + rampupTime: 0.1000
The simulation can then be executed using the `runBatteryJson` function.
jsonstruct.TimeStepping.totalTime = (1./jsonstruct.Control.CRate) .* 3600 .* 1.2;
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 3 Seconds, 375 Milliseconds +Solving timestep 02/45: 3 Seconds, 375 Milliseconds -> 6 Seconds, 750 Milliseconds +Solving timestep 03/45: 6 Seconds, 750 Milliseconds -> 13 Seconds, 500 Milliseconds +Solving timestep 04/45: 13 Seconds, 500 Milliseconds -> 27 Seconds +Solving timestep 05/45: 27 Seconds -> 54 Seconds +Solving timestep 06/45: 54 Seconds -> 108 Seconds +Solving timestep 07/45: 108 Seconds -> 216 Seconds +Solving timestep 08/45: 216 Seconds -> 324 Seconds +Solving timestep 09/45: 324 Seconds -> 432 Seconds +Solving timestep 10/45: 432 Seconds -> 540 Seconds +Solving timestep 11/45: 540 Seconds -> 648 Seconds +Solving timestep 12/45: 648 Seconds -> 756 Seconds +Solving timestep 13/45: 756 Seconds -> 864 Seconds +Solving timestep 14/45: 864 Seconds -> 972 Seconds +Solving timestep 15/45: 972 Seconds -> 1080 Seconds +Solving timestep 16/45: 1080 Seconds -> 1188 Seconds +Solving timestep 17/45: 1188 Seconds -> 1296 Seconds +Solving timestep 18/45: 1296 Seconds -> 1404 Seconds +Solving timestep 19/45: 1404 Seconds -> 1512 Seconds +Solving timestep 20/45: 1512 Seconds -> 1620 Seconds +Solving timestep 21/45: 1620 Seconds -> 1728 Seconds +Solving timestep 22/45: 1728 Seconds -> 1836 Seconds +Solving timestep 23/45: 1836 Seconds -> 1944 Seconds +Solving timestep 24/45: 1944 Seconds -> 2052 Seconds +Solving timestep 25/45: 2052 Seconds -> 2160 Seconds +Solving timestep 26/45: 2160 Seconds -> 2268 Seconds +Solving timestep 27/45: 2268 Seconds -> 2376 Seconds +Solving timestep 28/45: 2376 Seconds -> 2484 Seconds +Solving timestep 29/45: 2484 Seconds -> 2592 Seconds +Solving timestep 30/45: 2592 Seconds -> 2700 Seconds +Solving timestep 31/45: 2700 Seconds -> 2808 Seconds +Solving timestep 32/45: 2808 Seconds -> 2916 Seconds +Solving timestep 33/45: 2916 Seconds -> 3024 Seconds +Solving timestep 34/45: 3024 Seconds -> 3132 Seconds +Solving timestep 35/45: 3132 Seconds -> 3240 Seconds +Solving timestep 36/45: 3240 Seconds -> 3348 Seconds +Solving timestep 37/45: 3348 Seconds -> 3456 Seconds +Solving timestep 38/45: 3456 Seconds -> 3564 Seconds +Solving timestep 39/45: 3564 Seconds -> 1 Hour, 72 Seconds +*** Simulation complete. Solved 39 control steps in 2 Seconds, 933 Milliseconds (termination triggered by stopFunction) ***
The output matrix returns among other things the model and the states. The model contains information about the setup of the model and initial conditions, while states contains the results of the simulation at each timestep. The model can also output the energy, energy density and specific energy of the system.
disp(output);
model: [1×1 Battery] + inputparams: [1×1 BatteryInputParams] + schedule: [1×1 struct] + initstate: [1×1 struct] + time: [62×1 double] + E: [62×1 double] + I: [62×1 double] + states: {62×1 cell} + energy: [61×1 double] + energyDensity: [61×1 double] + specificEnergy: []
The `plotDashboard` function displays the Li-ion concentration profile in both the solids and electrolytes, as well as the electrode potentials at a specific time step.
The concentration profiles help understand how efficiently the intercalation (and deintercalation) processes have occurred. They can identify performance issues such as lithium plating or irreversible lithium loss.
plotDashboard(output.model, output.states, 'step', 40);
CC discharge voltage profile of NMC-Graphite cell:
Let's examine the voltage versus capacity, which represents the cell discharge profile. The voltage curve typically exhibits three distinct segments:
  • An initial drop, caused by internal resistances and the activation of electrochemical processes.
  • A plateau region where the voltage remains relatively flat. A longer plateau region in a battery's discharge profile is beneficial as it signifies that the voltage remains stable with minimal fluctuations across a substantial portion of its capacity, ensuring consistent and efficient energy delivery. The length and voltage of this plateau depend on the cell chemistry.
  • A final drop as the cell approaches its capacity limit.
% Extract the states from the output
states = output.states;
 
% Extract the time, voltage, and current quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
 
% Calculate the capacity
capacity = time .* current;
 
% Instantiate an empty figure
figure()
 
% Plot the discharge curve in the figure
plot((capacity / (hour * milli)), voltage, '-', 'linewidth', 3)
 
% Label the axes
xlabel('Capacity / mA \cdot h')
ylabel('Voltage / V')
 
% Adjust the axis limits
axis tight
Should we add the influence of C-rates here, copy from notebook example?

Influence of structural and material parameters on the capacity

For a given cell chemistry, the capacity of a cell is constrained by the number of lithium ions that can be safely intercalated and deintercalated from the electrodes. As cells age, they experience irreversible loss of lithium ions, which reduces their capacity over time. We will look at ageing mechanisms in a later chapter. Now, let's examine how material parameters affect cell capacity. One method to increase a cell's capacity, or the amount of energy stored, is by increasing the mass loading, which is the amount of active material per unit area of the electrode. This can be achieved by reducing the number of pores to increase the effective density or by increasing the mass fraction of the coating relative to other components, such as binders and additives. Alternatively, the thickness of the electrode can be increased while maintaining the same ratio of components, thereby adding more active material.
(i) Increasing the thickness of the anode / cathode:
Please set the electrode type (electrodeType) to be either 'NegativeElectrode' or 'PositiveElectrode' in the code below.
Anode:
We maintain the NMC cathode at a constant thickness while increasing the graphite electrode thickness from 32 micrometers to 72 micrometers to observe its impact on the capacity of our 1D system. Additionally, we switch to the CC-CV protocol, where we first charge the cell with constant current-constant voltage (CC-CV) and then discharge it at constant current.
We plot the discharge curve and observe that the capacity increases for 32 < 45 < 64 micrometers. However, beyond 64 micrometers, further increasing the thickness has no effect on capacity. This is because the amount of lithium in this scenario is limited by the NMC cathode. Although more lithium can initially intercalate into the thicker graphite, beyond a certain thickness, there is no additional lithium available to intercalate into the graphite that can be extracted during discharge. Therefore, the capacity remains unchanged.
Cathode:
We increase the thickness of the NMC cathode from 32 micrometers to 72 micrometers. In this scenario, the negative electrode is intentionally oversized relative to the capacity of the positive electrode by approximately 1.3 times. This design ensures that the capacity of the negative electrode exceeds that of the positive electrode, thereby helping to prevent dangerous lithium plating [4].
The thicker electrode indeed shows a notable increase in capacity. However, this comes with its own set of challenges. The longer diffusion pathways for lithium ions within the electrode, higher overpotentials, and extended paths for electronic transport all contribute to a reduced rate capacity. Despite the increase in energy density, these factors cause the capacity to diminish more rapidly at higher discharge rates, this makes it unsuitable for high power applications.
I cant see the effect of high C-rates, should I change some other parameter as well to see it in this toy example?
% create a vector of diffent thickness values
thickness = [32,45,64,72,].*micro;
markers={'-',"-",'o',"x"};
jsonstruct_default = parseBattmoJson('Examples/JsonDataFiles/sample_input.json');
cccv_control_protocol = parseBattmoJson('cccv_control_charge.json');%loads the CCCV protocol
jsonstruct_modified = mergeJsonStructs({cccv_control_protocol, jsonstruct_default});
parameter Control.controlPolicy is assigned twice with different values. we use the value from first jsonstruct. +parameter Control.lowerCutoffVoltage is assigned twice with different values. we use the value from first jsonstruct. +parameter SOC is assigned twice with different values. we use the value from first jsonstruct.
% instantiate and empty cell array to store the outputs of the simulations
output = cell(size(thickness));
 
% instantiate an empty figure
figure()
 
% Specify the type of electrode to modify
%electrodeType = 'NegativeElectrode'; % or 'PositiveElectrode'
electrodeType = 'PositiveElectrode'; % or 'PositiveElectrode'
dischargeRate ='high'; % or 'default'
 
% use a for-loop to iterate through the vector of thickness
for i = 1 : numel(thickness)
% Modify the value for the coating thickness for the specified electrode type
if strcmp(electrodeType, 'NegativeElectrode')
jsonstruct_modified.NegativeElectrode.Coating.thickness = thickness(i);
elseif strcmp(electrodeType, 'PositiveElectrode')
if strcmp(dischargeRate, 'high')
jsonstruct_modified.CRate=10;
end
jsonstruct_modified.PositiveElectrode.Coating.thickness = thickness(i);
jsonstruct_modified.NegativeElectrode.Coating.thickness = thickness(i)*1.3;
else
error('Unknown electrode type specified');
end
% run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct_modified);
% retrieve the states from the simulation result
states = output{i}.states;
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
% calculate the capacity
capacity = time .* current;
% plot only the discharge part of the curve in the figure to observe the capacity
start_index = find(capacity > 0, 1, 'first');
end_index = find(capacity > 0, 1, 'last');
plot((capacity(start_index:end_index)/(hour*milli)), voltage(start_index:end_index), markers{i}, 'linewidth', 3);
hold on
end
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds +Solving timestep 02/45: 6 Seconds, 187 Milliseconds -> 12 Seconds, 375 Milliseconds +Solving timestep 03/45: 12 Seconds, 375 Milliseconds -> 24 Seconds, 750 Milliseconds +Solving timestep 04/45: 24 Seconds, 750 Milliseconds -> 49 Seconds, 500 Milliseconds +Solving timestep 05/45: 49 Seconds, 500 Milliseconds -> 99 Seconds +Solving timestep 06/45: 99 Seconds -> 198 Seconds +Solving timestep 07/45: 198 Seconds -> 396 Seconds +Solving timestep 08/45: 396 Seconds -> 594 Seconds +Solving timestep 09/45: 594 Seconds -> 792 Seconds +Solving timestep 10/45: 792 Seconds -> 990 Seconds +Solving timestep 11/45: 990 Seconds -> 1188 Seconds +Solving timestep 12/45: 1188 Seconds -> 1386 Seconds +Solving timestep 13/45: 1386 Seconds -> 1584 Seconds +Solving timestep 14/45: 1584 Seconds -> 1782 Seconds +Solving timestep 15/45: 1782 Seconds -> 1980 Seconds +Switch control type from CV_charge2 to CC_discharge1 +Solving timestep 16/45: 1980 Seconds -> 2178 Seconds +Solving timestep 17/45: 2178 Seconds -> 2376 Seconds +Solving timestep 18/45: 2376 Seconds -> 2574 Seconds +Solving timestep 19/45: 2574 Seconds -> 2772 Seconds +Solving timestep 20/45: 2772 Seconds -> 2970 Seconds +Solving timestep 21/45: 2970 Seconds -> 3168 Seconds +Solving timestep 22/45: 3168 Seconds -> 3366 Seconds +Solving timestep 23/45: 3366 Seconds -> 3564 Seconds +Solving timestep 24/45: 3564 Seconds -> 1 Hour, 162 Seconds +Solving timestep 25/45: 1 Hour, 162 Seconds -> 1 Hour, 360 Seconds +Solving timestep 26/45: 1 Hour, 360 Seconds -> 1 Hour, 558 Seconds +Solving timestep 27/45: 1 Hour, 558 Seconds -> 1 Hour, 756 Seconds +Solving timestep 28/45: 1 Hour, 756 Seconds -> 1 Hour, 954 Seconds +Solving timestep 29/45: 1 Hour, 954 Seconds -> 1 Hour, 1152 Seconds +Solving timestep 30/45: 1 Hour, 1152 Seconds -> 1 Hour, 1350 Seconds +Solving timestep 31/45: 1 Hour, 1350 Seconds -> 1 Hour, 1548 Seconds +Solving timestep 32/45: 1 Hour, 1548 Seconds -> 1 Hour, 1746 Seconds +switch control from CC_discharge1 to CC_discharge2 +Solving timestep 33/45: 1 Hour, 1746 Seconds -> 1 Hour, 1944 Seconds +Switch control type from CC_discharge2 to CC_charge1 +*** Simulation complete. Solved 33 control steps in 2 Seconds, 504 Milliseconds (termination triggered by stopFunction) ***
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds +Solving timestep 02/45: 6 Seconds, 187 Milliseconds -> 12 Seconds, 375 Milliseconds +Solving timestep 03/45: 12 Seconds, 375 Milliseconds -> 24 Seconds, 750 Milliseconds +Solving timestep 04/45: 24 Seconds, 750 Milliseconds -> 49 Seconds, 500 Milliseconds +Solving timestep 05/45: 49 Seconds, 500 Milliseconds -> 99 Seconds +Solving timestep 06/45: 99 Seconds -> 198 Seconds +Solving timestep 07/45: 198 Seconds -> 396 Seconds +Solving timestep 08/45: 396 Seconds -> 594 Seconds +Solving timestep 09/45: 594 Seconds -> 792 Seconds +Solving timestep 10/45: 792 Seconds -> 990 Seconds +Solving timestep 11/45: 990 Seconds -> 1188 Seconds +Solving timestep 12/45: 1188 Seconds -> 1386 Seconds +Solving timestep 13/45: 1386 Seconds -> 1584 Seconds +Solving timestep 14/45: 1584 Seconds -> 1782 Seconds +Switch control type from CV_charge2 to CC_discharge1 +Solving timestep 15/45: 1782 Seconds -> 1980 Seconds +Solving timestep 16/45: 1980 Seconds -> 2178 Seconds +Solving timestep 17/45: 2178 Seconds -> 2376 Seconds +Solving timestep 18/45: 2376 Seconds -> 2574 Seconds +Solving timestep 19/45: 2574 Seconds -> 2772 Seconds +Solving timestep 20/45: 2772 Seconds -> 2970 Seconds +Solving timestep 21/45: 2970 Seconds -> 3168 Seconds +Solving timestep 22/45: 3168 Seconds -> 3366 Seconds +Solving timestep 23/45: 3366 Seconds -> 3564 Seconds +Solving timestep 24/45: 3564 Seconds -> 1 Hour, 162 Seconds +Solving timestep 25/45: 1 Hour, 162 Seconds -> 1 Hour, 360 Seconds +Solving timestep 26/45: 1 Hour, 360 Seconds -> 1 Hour, 558 Seconds +Solving timestep 27/45: 1 Hour, 558 Seconds -> 1 Hour, 756 Seconds +Solving timestep 28/45: 1 Hour, 756 Seconds -> 1 Hour, 954 Seconds +Solving timestep 29/45: 1 Hour, 954 Seconds -> 1 Hour, 1152 Seconds +Solving timestep 30/45: 1 Hour, 1152 Seconds -> 1 Hour, 1350 Seconds +Solving timestep 31/45: 1 Hour, 1350 Seconds -> 1 Hour, 1548 Seconds +switch control from CC_discharge1 to CC_discharge2 +Switch control type from CC_discharge2 to CC_charge1 +*** Simulation complete. Solved 31 control steps in 2 Seconds, 384 Milliseconds (termination triggered by stopFunction) ***
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds +Solving timestep 02/45: 6 Seconds, 187 Milliseconds -> 12 Seconds, 375 Milliseconds +Solving timestep 03/45: 12 Seconds, 375 Milliseconds -> 24 Seconds, 750 Milliseconds +Solving timestep 04/45: 24 Seconds, 750 Milliseconds -> 49 Seconds, 500 Milliseconds +Solving timestep 05/45: 49 Seconds, 500 Milliseconds -> 99 Seconds +Solving timestep 06/45: 99 Seconds -> 198 Seconds +Solving timestep 07/45: 198 Seconds -> 396 Seconds +Solving timestep 08/45: 396 Seconds -> 594 Seconds +Solving timestep 09/45: 594 Seconds -> 792 Seconds +Solving timestep 10/45: 792 Seconds -> 990 Seconds +Solving timestep 11/45: 990 Seconds -> 1188 Seconds +Solving timestep 12/45: 1188 Seconds -> 1386 Seconds +Solving timestep 13/45: 1386 Seconds -> 1584 Seconds +Solving timestep 14/45: 1584 Seconds -> 1782 Seconds +Switch control type from CV_charge2 to CC_discharge1 +Solving timestep 15/45: 1782 Seconds -> 1980 Seconds +Solving timestep 16/45: 1980 Seconds -> 2178 Seconds +Solving timestep 17/45: 2178 Seconds -> 2376 Seconds +Solving timestep 18/45: 2376 Seconds -> 2574 Seconds +Solving timestep 19/45: 2574 Seconds -> 2772 Seconds +Solving timestep 20/45: 2772 Seconds -> 2970 Seconds +Solving timestep 21/45: 2970 Seconds -> 3168 Seconds +Solving timestep 22/45: 3168 Seconds -> 3366 Seconds +Solving timestep 23/45: 3366 Seconds -> 3564 Seconds +Solving timestep 24/45: 3564 Seconds -> 1 Hour, 162 Seconds +Solving timestep 25/45: 1 Hour, 162 Seconds -> 1 Hour, 360 Seconds +Solving timestep 26/45: 1 Hour, 360 Seconds -> 1 Hour, 558 Seconds +Solving timestep 27/45: 1 Hour, 558 Seconds -> 1 Hour, 756 Seconds +Solving timestep 28/45: 1 Hour, 756 Seconds -> 1 Hour, 954 Seconds +Solving timestep 29/45: 1 Hour, 954 Seconds -> 1 Hour, 1152 Seconds +Solving timestep 30/45: 1 Hour, 1152 Seconds -> 1 Hour, 1350 Seconds +Solving timestep 31/45: 1 Hour, 1350 Seconds -> 1 Hour, 1548 Seconds +switch control from CC_discharge1 to CC_discharge2 +Switch control type from CC_discharge2 to CC_charge1 +*** Simulation complete. Solved 31 control steps in 2 Seconds, 634 Milliseconds (termination triggered by stopFunction) ***
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds +Solving timestep 02/45: 6 Seconds, 187 Milliseconds -> 12 Seconds, 375 Milliseconds +Solving timestep 03/45: 12 Seconds, 375 Milliseconds -> 24 Seconds, 750 Milliseconds +Solving timestep 04/45: 24 Seconds, 750 Milliseconds -> 49 Seconds, 500 Milliseconds +Solving timestep 05/45: 49 Seconds, 500 Milliseconds -> 99 Seconds +Solving timestep 06/45: 99 Seconds -> 198 Seconds +Solving timestep 07/45: 198 Seconds -> 396 Seconds +Solving timestep 08/45: 396 Seconds -> 594 Seconds +Solving timestep 09/45: 594 Seconds -> 792 Seconds +Solving timestep 10/45: 792 Seconds -> 990 Seconds +Solving timestep 11/45: 990 Seconds -> 1188 Seconds +Solving timestep 12/45: 1188 Seconds -> 1386 Seconds +Solving timestep 13/45: 1386 Seconds -> 1584 Seconds +Switch control type from CV_charge2 to CC_discharge1 +Solving timestep 14/45: 1584 Seconds -> 1782 Seconds +Solving timestep 15/45: 1782 Seconds -> 1980 Seconds +Solving timestep 16/45: 1980 Seconds -> 2178 Seconds +Solving timestep 17/45: 2178 Seconds -> 2376 Seconds +Solving timestep 18/45: 2376 Seconds -> 2574 Seconds +Solving timestep 19/45: 2574 Seconds -> 2772 Seconds +Solving timestep 20/45: 2772 Seconds -> 2970 Seconds +Solving timestep 21/45: 2970 Seconds -> 3168 Seconds +Solving timestep 22/45: 3168 Seconds -> 3366 Seconds +Solving timestep 23/45: 3366 Seconds -> 3564 Seconds +Solving timestep 24/45: 3564 Seconds -> 1 Hour, 162 Seconds +Solving timestep 25/45: 1 Hour, 162 Seconds -> 1 Hour, 360 Seconds +Solving timestep 26/45: 1 Hour, 360 Seconds -> 1 Hour, 558 Seconds +Solving timestep 27/45: 1 Hour, 558 Seconds -> 1 Hour, 756 Seconds +Solving timestep 28/45: 1 Hour, 756 Seconds -> 1 Hour, 954 Seconds +Solving timestep 29/45: 1 Hour, 954 Seconds -> 1 Hour, 1152 Seconds +switch control from CC_discharge1 to CC_discharge2 +Switch control type from CC_discharge2 to CC_charge1 +*** Simulation complete. Solved 29 control steps in 2 Seconds, 342 Milliseconds (termination triggered by stopFunction) ***
hold off
xlabel('Capacity / mA \cdot h')
ylabel('Voltage / V')
legend( 't_{elde} = 32 µm', 't_{elde} = 45 µm','t_{elde} = 64 µm','t_{elde} = 72 µm')
(ii) Increasing the effective density of electrode coating
The effective density of the electrode coating comprises both the active material and the pores. A crucial aspect of cell design is optimizing the ratio of pores to active material. While increasing the amount of active material enhances capacity and energy density by providing more lithium atoms, a less porous electrode reduces the electrolyte volume fraction. This reduction can impair the rate of ion transfer, leading to lower power densities and potential issues with the electrode’s mechanical stability. The density of the cathode, NMC is typically around 4650 kg/m³. To assess the impact of effective density on capacity, we vary it from 3500 to 3900 kg/m³. As discussed earlier, the electrolyte fills the pores between the electrode active particles. Therefore, increasing the effective density indefinitely is not possible, as it would reduce the electrolyte volume and hinder the transport of Li-ions.
jsonstruct = parseBattmoJson('Examples/JsonDataFiles/sample_input.json');% parse the sample_input.json for default values.
effective_density = [3500,3700,3900,];
markers={'-',"-",'-'};
output = cell(size(effective_density));
% instantiate an empty figure
figure()
 
% use a for-loop to iterate through the vector of effective density
for i = 1 : numel(effective_density)
% Modify the value for the coating thickness for the positiv electrode
jsonstruct.PositiveElectrode.Coating.effectiveDensity = effective_density(i);
% run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct);
% retrieve the states from the simulation result
states = output{i}.states;
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
% calculate the capacity
capacity = time .* current;
% plot capacity vs voltage
plot((capacity/(hour*milli)), voltage, markers{i}, 'linewidth', 3);
hold on
end
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds +Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds +Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds +Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds +Solving timestep 05/45: 31 Seconds, 500 Milliseconds -> 63 Seconds +Solving timestep 06/45: 63 Seconds -> 126 Seconds +Solving timestep 07/45: 126 Seconds -> 252 Seconds +Solving timestep 08/45: 252 Seconds -> 378 Seconds +Solving timestep 09/45: 378 Seconds -> 504 Seconds +Solving timestep 10/45: 504 Seconds -> 630 Seconds +Solving timestep 11/45: 630 Seconds -> 756 Seconds +Solving timestep 12/45: 756 Seconds -> 882 Seconds +Solving timestep 13/45: 882 Seconds -> 1008 Seconds +Solving timestep 14/45: 1008 Seconds -> 1134 Seconds +Solving timestep 15/45: 1134 Seconds -> 1260 Seconds +Solving timestep 16/45: 1260 Seconds -> 1386 Seconds +Solving timestep 17/45: 1386 Seconds -> 1512 Seconds +Solving timestep 18/45: 1512 Seconds -> 1638 Seconds +Solving timestep 19/45: 1638 Seconds -> 1764 Seconds +Solving timestep 20/45: 1764 Seconds -> 1890 Seconds +Solving timestep 21/45: 1890 Seconds -> 2016 Seconds +Solving timestep 22/45: 2016 Seconds -> 2142 Seconds +Solving timestep 23/45: 2142 Seconds -> 2268 Seconds +Solving timestep 24/45: 2268 Seconds -> 2394 Seconds +Solving timestep 25/45: 2394 Seconds -> 2520 Seconds +Solving timestep 26/45: 2520 Seconds -> 2646 Seconds +Solving timestep 27/45: 2646 Seconds -> 2772 Seconds +Solving timestep 28/45: 2772 Seconds -> 2898 Seconds +Solving timestep 29/45: 2898 Seconds -> 3024 Seconds +Solving timestep 30/45: 3024 Seconds -> 3150 Seconds +Solving timestep 31/45: 3150 Seconds -> 3276 Seconds +Solving timestep 32/45: 3276 Seconds -> 3402 Seconds +Solving timestep 33/45: 3402 Seconds -> 3528 Seconds +Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds +*** Simulation complete. Solved 34 control steps in 2 Seconds, 371 Milliseconds (termination triggered by stopFunction) *** +Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds +Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds +Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds +Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds +Solving timestep 05/45: 31 Seconds, 500 Milliseconds -> 63 Seconds +Solving timestep 06/45: 63 Seconds -> 126 Seconds +Solving timestep 07/45: 126 Seconds -> 252 Seconds +Solving timestep 08/45: 252 Seconds -> 378 Seconds +Solving timestep 09/45: 378 Seconds -> 504 Seconds +Solving timestep 10/45: 504 Seconds -> 630 Seconds +Solving timestep 11/45: 630 Seconds -> 756 Seconds +Solving timestep 12/45: 756 Seconds -> 882 Seconds +Solving timestep 13/45: 882 Seconds -> 1008 Seconds +Solving timestep 14/45: 1008 Seconds -> 1134 Seconds +Solving timestep 15/45: 1134 Seconds -> 1260 Seconds +Solving timestep 16/45: 1260 Seconds -> 1386 Seconds +Solving timestep 17/45: 1386 Seconds -> 1512 Seconds +Solving timestep 18/45: 1512 Seconds -> 1638 Seconds +Solving timestep 19/45: 1638 Seconds -> 1764 Seconds +Solving timestep 20/45: 1764 Seconds -> 1890 Seconds +Solving timestep 21/45: 1890 Seconds -> 2016 Seconds +Solving timestep 22/45: 2016 Seconds -> 2142 Seconds +Solving timestep 23/45: 2142 Seconds -> 2268 Seconds +Solving timestep 24/45: 2268 Seconds -> 2394 Seconds +Solving timestep 25/45: 2394 Seconds -> 2520 Seconds +Solving timestep 26/45: 2520 Seconds -> 2646 Seconds +Solving timestep 27/45: 2646 Seconds -> 2772 Seconds +Solving timestep 28/45: 2772 Seconds -> 2898 Seconds +Solving timestep 29/45: 2898 Seconds -> 3024 Seconds +Solving timestep 30/45: 3024 Seconds -> 3150 Seconds +Solving timestep 31/45: 3150 Seconds -> 3276 Seconds +Solving timestep 32/45: 3276 Seconds -> 3402 Seconds +Solving timestep 33/45: 3402 Seconds -> 3528 Seconds +Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds +Solving timestep 35/45: 1 Hour, 54 Seconds -> 1 Hour, 180 Seconds +*** Simulation complete. Solved 35 control steps in 2 Seconds, 358 Milliseconds (termination triggered by stopFunction) *** +Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds +Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds +Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds +Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds +Solving timestep 05/45: 31 Seconds, 500 Milliseconds -> 63 Seconds +Solving timestep 06/45: 63 Seconds -> 126 Seconds +Solving timestep 07/45: 126 Seconds -> 252 Seconds +Solving timestep 08/45: 252 Seconds -> 378 Seconds +Solving timestep 09/45: 378 Seconds -> 504 Seconds +Solving timestep 10/45: 504 Seconds -> 630 Seconds +Solving timestep 11/45: 630 Seconds -> 756 Seconds +Solving timestep 12/45: 756 Seconds -> 882 Seconds +Solving timestep 13/45: 882 Seconds -> 1008 Seconds +Solving timestep 14/45: 1008 Seconds -> 1134 Seconds +Solving timestep 15/45: 1134 Seconds -> 1260 Seconds +Solving timestep 16/45: 1260 Seconds -> 1386 Seconds +Solving timestep 17/45: 1386 Seconds -> 1512 Seconds +Solving timestep 18/45: 1512 Seconds -> 1638 Seconds +Solving timestep 19/45: 1638 Seconds -> 1764 Seconds +Solving timestep 20/45: 1764 Seconds -> 1890 Seconds +Solving timestep 21/45: 1890 Seconds -> 2016 Seconds +Solving timestep 22/45: 2016 Seconds -> 2142 Seconds +Solving timestep 23/45: 2142 Seconds -> 2268 Seconds +Solving timestep 24/45: 2268 Seconds -> 2394 Seconds +Solving timestep 25/45: 2394 Seconds -> 2520 Seconds +Solving timestep 26/45: 2520 Seconds -> 2646 Seconds +Solving timestep 27/45: 2646 Seconds -> 2772 Seconds +Solving timestep 28/45: 2772 Seconds -> 2898 Seconds +Solving timestep 29/45: 2898 Seconds -> 3024 Seconds +Solving timestep 30/45: 3024 Seconds -> 3150 Seconds +Solving timestep 31/45: 3150 Seconds -> 3276 Seconds +Solving timestep 32/45: 3276 Seconds -> 3402 Seconds +Solving timestep 33/45: 3402 Seconds -> 3528 Seconds +Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds +Solving timestep 35/45: 1 Hour, 54 Seconds -> 1 Hour, 180 Seconds +Solving timestep 36/45: 1 Hour, 180 Seconds -> 1 Hour, 306 Seconds +*** Simulation complete. Solved 36 control steps in 2 Seconds, 366 Milliseconds (termination triggered by stopFunction) ***
hold off
xlabel('Capacity / mA \cdot h')
ylabel('Voltage / V')
legend( 'eff-dens = 3500', 'eff-dens = 3700','eff-dens = 3900')
(iii) Influence of material chemistry on capacity:
Let's now replace the cathode material NMC 811 with LiFePO4 (LFP) [2]. Different cathode materials have varying lithium saturation concentrations, reaction rate constants, and densities. LFP has a lower density and a lower lithium ion saturation concentration compared to NMC 811. As a result, we expect to observe a lower capacity with LFP. Additionally, the shape of the discharge curve will differ, reflecting the distinct electrochemical behavior of LFP compared to NMC 811.
cathode_materials = ["NMC","LFP"];
markers={'-',"-",};
output = cell(size(cathode_materials));
% instantiate an empty figure
figure()
% use a for-loop to iterate through the vector of cathode materials
for i = 1:numel(cathode_materials)
% Parse the sample input JSON for default values, not the default
% material used here is NMC
jsonstruct = parseBattmoJson('Examples/JsonDataFiles/sample_input.json');
if cathode_materials(i) == "LFP"
% Load the LFP json parameters
lfp = parseBattmoJson('ParameterData/MaterialProperties/LFP/LFP.json');
jsonstruct_lfp.PositiveElectrode.Coating.ActiveMaterial.Interface = lfp;
jsonstruct = mergeJsonStructs({jsonstruct_lfp, jsonstruct});
jsonstruct.PositiveElectrode.Coating.ActiveMaterial.density = jsonstruct.PositiveElectrode.Coating.ActiveMaterial.Interface.density;
% the effective density also needs to be changed here for
% the new material
jsonstruct.PositiveElectrode.Coating.effectiveDensity = 2709;
end
% Run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct);
% Retrieve the states from the simulation result
states = output{i}.states;
% Extract time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
% Calculate the capacity
capacity = time .* current;
% plot capacity vs voltage
plot((capacity/(hour*milli)), voltage, markers{i}, 'linewidth', 3);
hold on
end
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds +Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds +Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds +Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds +Solving timestep 05/45: 31 Seconds, 500 Milliseconds -> 63 Seconds +Solving timestep 06/45: 63 Seconds -> 126 Seconds +Solving timestep 07/45: 126 Seconds -> 252 Seconds +Solving timestep 08/45: 252 Seconds -> 378 Seconds +Solving timestep 09/45: 378 Seconds -> 504 Seconds +Solving timestep 10/45: 504 Seconds -> 630 Seconds +Solving timestep 11/45: 630 Seconds -> 756 Seconds +Solving timestep 12/45: 756 Seconds -> 882 Seconds +Solving timestep 13/45: 882 Seconds -> 1008 Seconds +Solving timestep 14/45: 1008 Seconds -> 1134 Seconds +Solving timestep 15/45: 1134 Seconds -> 1260 Seconds +Solving timestep 16/45: 1260 Seconds -> 1386 Seconds +Solving timestep 17/45: 1386 Seconds -> 1512 Seconds +Solving timestep 18/45: 1512 Seconds -> 1638 Seconds +Solving timestep 19/45: 1638 Seconds -> 1764 Seconds +Solving timestep 20/45: 1764 Seconds -> 1890 Seconds +Solving timestep 21/45: 1890 Seconds -> 2016 Seconds +Solving timestep 22/45: 2016 Seconds -> 2142 Seconds +Solving timestep 23/45: 2142 Seconds -> 2268 Seconds +Solving timestep 24/45: 2268 Seconds -> 2394 Seconds +Solving timestep 25/45: 2394 Seconds -> 2520 Seconds +Solving timestep 26/45: 2520 Seconds -> 2646 Seconds +Solving timestep 27/45: 2646 Seconds -> 2772 Seconds +Solving timestep 28/45: 2772 Seconds -> 2898 Seconds +Solving timestep 29/45: 2898 Seconds -> 3024 Seconds +Solving timestep 30/45: 3024 Seconds -> 3150 Seconds +Solving timestep 31/45: 3150 Seconds -> 3276 Seconds +Solving timestep 32/45: 3276 Seconds -> 3402 Seconds +Solving timestep 33/45: 3402 Seconds -> 3528 Seconds +Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds +*** Simulation complete. Solved 34 control steps in 2 Seconds, 451 Milliseconds (termination triggered by stopFunction) ***
parameter PositiveElectrode.Coating.ActiveMaterial.Interface.saturationConcentration is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.Coating.ActiveMaterial.Interface.volumetricSurfaceArea is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.Coating.ActiveMaterial.Interface.density is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.Coating.ActiveMaterial.Interface.activationEnergyOfReaction is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.Coating.ActiveMaterial.Interface.reactionRateConstant is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.Coating.ActiveMaterial.Interface.guestStoichiometry100 is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.Coating.ActiveMaterial.Interface.guestStoichiometry0 is assigned twice with different values. we use the value from first jsonstruct. +parameter PositiveElectrode.Coating.ActiveMaterial.Interface.openCircuitPotential.functionname is assigned twice with different values. we use the value from first jsonstruct.
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds +Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds +Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds +Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds +Solving timestep 05/45: 31 Seconds, 500 Milliseconds -> 63 Seconds +Solving timestep 06/45: 63 Seconds -> 126 Seconds +Solving timestep 07/45: 126 Seconds -> 252 Seconds +Solving timestep 08/45: 252 Seconds -> 378 Seconds +Solving timestep 09/45: 378 Seconds -> 504 Seconds +Solving timestep 10/45: 504 Seconds -> 630 Seconds +Solving timestep 11/45: 630 Seconds -> 756 Seconds +Solving timestep 12/45: 756 Seconds -> 882 Seconds +Solving timestep 13/45: 882 Seconds -> 1008 Seconds +Solving timestep 14/45: 1008 Seconds -> 1134 Seconds +Solving timestep 15/45: 1134 Seconds -> 1260 Seconds +Solving timestep 16/45: 1260 Seconds -> 1386 Seconds +Solving timestep 17/45: 1386 Seconds -> 1512 Seconds +Solving timestep 18/45: 1512 Seconds -> 1638 Seconds +Solving timestep 19/45: 1638 Seconds -> 1764 Seconds +Solving timestep 20/45: 1764 Seconds -> 1890 Seconds +Solving timestep 21/45: 1890 Seconds -> 2016 Seconds +Solving timestep 22/45: 2016 Seconds -> 2142 Seconds +Solving timestep 23/45: 2142 Seconds -> 2268 Seconds +Solving timestep 24/45: 2268 Seconds -> 2394 Seconds +Solving timestep 25/45: 2394 Seconds -> 2520 Seconds +Solving timestep 26/45: 2520 Seconds -> 2646 Seconds +Solving timestep 27/45: 2646 Seconds -> 2772 Seconds +Solving timestep 28/45: 2772 Seconds -> 2898 Seconds +Solving timestep 29/45: 2898 Seconds -> 3024 Seconds +Solving timestep 30/45: 3024 Seconds -> 3150 Seconds +Solving timestep 31/45: 3150 Seconds -> 3276 Seconds +Solver did not converge in 10 iterations for timestep of length 126 Seconds. Cutting timestep. +Solving timestep 32/45: 3276 Seconds -> 3402 Seconds +Solving timestep 33/45: 3402 Seconds -> 3528 Seconds +Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds +*** Simulation complete. Solved 34 control steps in 2 Seconds, 358 Milliseconds (termination triggered by stopFunction) ***
hold off
xlabel('Capacity / mA \cdot h')
ylabel('Voltage / V')
legend( 'NMC', 'LFP')
As expected, we observe that the capacity of LiFePO4 (LFP) is lower compared to NMC, and the maximum open-circuit voltage (OCP) for the LFP-graphite system is also less than that of the NMC cathode. Additionally, the discharge curve for LFP is flatter, indicating that the voltage remains more constant throughout the discharge cycle compared to NMC.
Should we add a comment explaining when a 1D model is sufficient, which parameters can be tested in 1D, and when it is necessary to use a full 3D model? Can you let me know?
References
[1] Flores Cedeño, E. J. (2020). Development of operando diagnostics for Li-ion cathodes by Raman spectroscopy (Doctoral thesis, ETH Zurich).
[2] Safari, M., et al. "Multimodal physics–based aging model for life prediction of Li–ion batteries." Journal of The Electrochemical Society 156.3 (2008): A145.
[3] Xu, Meng, et al. "A pseudo three–dimensional electrochemical–thermal model of a prismatic LiFePO4 battery during discharge process." Energy 80 (2015): 303–317.
[4]Karin Kleiner, Peter Jakes, Sebastian Scharner, Verena Liebau, Helmut Ehrenberg, "Changes of the balancing between anode and cathode due to fatigue in commercial lithium-ion cells," Journal of Power Sources, Volume 317, 2016, Pages 25-34, ISSN 0378-7753, https://doi.org/10.1016/j.jpowsour.2016.03.049.
+
+ +
\ No newline at end of file diff --git a/_static/notebooks/tutorial_1_a_simple_p2d_model_live.html b/_static/notebooks/tutorial_1_a_simple_p2d_model_live.html index bbeed4ac..25e5fabc 100644 --- a/_static/notebooks/tutorial_1_a_simple_p2d_model_live.html +++ b/_static/notebooks/tutorial_1_a_simple_p2d_model_live.html @@ -2,10 +2,13 @@ Tutorial 1 - Your First BattMo Model

Tutorial 1 - Your First BattMo Model

Introduction

Let’s make your first BattMo model! ?⚡?
In this tutorial, we will build, run, and visualize a simple pseudo-two-dimensional (P2D) simulation for a Li-ion battery cell. We will first introduce how BattMo handles model parameters, then run the simulation, dashboard the results, and explore the details of the simulation output.

Define Parameters

BattMo uses JSON to manage parameters. This allows you to easily save, document, and share complete parameter sets from specific simulations. We have used long and explicit key names for good readability. If you are new to JSON, you can learn more about it here. Details on the BattMo specification are available in the JSON input specification.
For this example, we provide a sample JSON file sample_input.json that describes a simple NMC-Graphite cell.
We load and parse the JSON input file into BattMo using the command:
jsonstruct = parseBattmoJson('Examples/JsonDataFiles/sample_input.json');
This transforms the parameter data as a MATLAB structure jsonstruct that is used to setup the simulation. We can explore the structure within the MATLAB Command Window by navigating the different levels of the structure.
disp(jsonstruct)
use_thermal: 0 +.S9 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: 700; text-align: left; } +.S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S11 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S12 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S13 { margin: 15px 10px 5px 4px; padding: 0px; line-height: 18px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 17px; font-weight: 700; text-align: left; }

Tutorial 1 - Your First BattMo Model

Introduction

Let’s make your first BattMo model! ?⚡?
In this tutorial, we will build, run, and visualize a simple pseudo-two-dimensional (P2D) simulation for a Li-ion battery cell. We will first introduce how BattMo handles model parameters, then run the simulation, dashboard the results, and explore the details of the simulation output.

Define Parameters

BattMo uses JSON to manage input parameters. This allows you to easily save, document, and share complete parameter sets from specific simulations. We have used long and explicit key names for good readability. If you are new to JSON, you can learn more about it here.
For this example, we provide a sample JSON file sample_input.json describing a a simple NMC-Graphite cell [1] with LiPF6 [2] as the electrolyte in a 1D geometry according to BattMo's JSON input specification.
This sample_input.json contains the following schemas, needed to run the simulation.
In this tutorial, we simulate the constant current discharge profile for the cell with a CRate of 1 and a SoC=1.
We load and parse the JSON input file into BattMo using the command:
jsonstruct = parseBattmoJson('Examples/JsonDataFiles/sample_input.json');
This transforms the parameter data as a MATLAB structure jsonstruct that is used to setup the simulation. We can explore the structure within the MATLAB Command Window by navigating the different levels of the structure.
disp(jsonstruct)
use_thermal: 0 include_current_collectors: 0 - Geometry: [1x1 struct] - NegativeElectrode: [1x1 struct] - PositiveElectrode: [1x1 struct] - Separator: [1x1 struct] - Control: [1x1 struct] - Electrolyte: [1x1 struct] - ThermalModel: [1x1 struct] - TimeStepping: [1x1 struct] - Output: [1x1 struct] + Geometry: [1×1 struct] + NegativeElectrode: [1×1 struct] + PositiveElectrode: [1×1 struct] + Separator: [1×1 struct] + Control: [1×1 struct] + Electrolyte: [1×1 struct] + ThermalModel: [1×1 struct] + TimeStepping: [1×1 struct] + Output: [1×1 struct] SOC: 0.9900 - initT: 298.1500
We can see that the structure contains various attributes for the different aspects of the model. For example, if we want to know the thickness of the negative electrode coating, we can navigate the hierarchy to find it:
disp(jsonstruct.NegativeElectrode.Coating.thickness)
6.4000e-05
Unless otherwise specified, BattMo uses SI base units for physical quantities.

Run Simulation

We can run the simulation with the command:
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds + initT: 298.1500
We can see that the structure contains various attributes for the different schema. For example, if we want to know the thickness of the negative electrode coating, we can navigate the hierarchy to find it:
disp(jsonstruct.NegativeElectrode.Coating.thickness)
6.4000e-05
Let's also check the time step parameter, we see that the maximum time for the simulation is set at ~5040 seconds with 40 timesteps.
disp(jsonstruct.TimeStepping)
totalTime: 5040 + numberOfTimeSteps: 40 + useRampup: 1
Unless otherwise specified, BattMo uses SI base units for physical quantities.

Run Simulation

we can now use the runBatteryJson function to simulate the constant current current discharge (half cycle) of the cell. At each time step in the simulation, BattMo computes the Li-ion concentrations in the electrodes and the electrolyte as well as the electrode and electrolyte potentials.
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds @@ -79,27 +85,27 @@ Solving timestep 32/45: 3276 Seconds -> 3402 Seconds Solving timestep 33/45: 3402 Seconds -> 3528 Seconds Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds -*** Simulation complete. Solved 34 control steps in 8 Seconds, 822 Milliseconds (termination triggered by stopFunction) ***

Show the Dashboard

We can dashboard the main results using plotDashboard. Here, for example at time step 10,
plotDashboard(output.model, output.states, 'step', 10);
The left 3 columns of the dashboard shows the profiles for the main state quantities (concentration and electric potential) in the negative electrode, electrolyte, and positive electrode. The rightmost column shows the calculated cell current and voltage. In the following subsections, we will explore how to access and plot this data from the simulation output.

Explore the Output

The output structure returns among other thing the model and the states.
disp(output)
model: [1x1 Battery] - inputparams: [1x1 BatteryInputParams] - schedule: [1x1 struct] - initstate: [1x1 struct] - time: [57x1 double] - E: [57x1 double] - I: [57x1 double] - states: {57x1 cell} - energy: [56x1 double] - energyDensity: [56x1 double] - specificEnergy: []
The model contains information about the setup of the model and initial conditions, while states contains the results of the simulation at each timestep. Plotting the simulation results requires information about the grid (i.e. what is the position where the quantity is calculated?) and the state (i.e. what is the value of the quantity in that position at a given time?).

Explore the Grid

The grid (or mesh) is one of the most used properties of the model, which can be accessed with the command:
disp(output.model.grid)
cells: [1x1 struct] - faces: [1x1 struct] - nodes: [1x1 struct] +*** Simulation complete. Solved 34 control steps in 2 Seconds, 358 Milliseconds (termination triggered by stopFunction) ***

Show the Dashboard

We can dashboard the main results using plotDashboard. The left 3 columns of the dashboard shows the profiles for the main state quantities (concentration and electric potential) in the negative electrode, electrolyte, and positive electrode. The rightmost column shows the calculated cell current and voltage. Here, for example at time step 40, we can see that the Li has fully de-intercalated from the negative electrode (graphite) and intercalated almost uniformly into the positive electrode (NMC). The electric potentials at the electrodes also show oxidation at the anode and reduction at the cathode. We also see the cell voltage drop from 4.1 V to 2.4 V at the end of the discharge cycle.
plotDashboard(output.model, output.states, 'step', 40);
In the following subsections, we will explore how to access and plot this data from the simulation output.

Explore the Output

The output structure returns among other thing the model and the states.
disp(output)
model: [1×1 Battery] + inputparams: [1×1 BatteryInputParams] + schedule: [1×1 struct] + initstate: [1×1 struct] + time: [58×1 double] + E: [58×1 double] + I: [58×1 double] + states: {58×1 cell} + energy: [57×1 double] + energyDensity: [57×1 double] + specificEnergy: []
The model contains information about the setup of the model and initial conditions, while states contains the results of the simulation at each timestep. Plotting the simulation results requires information about the grid (i.e. what is the position where the quantity is calculated?) and the state (i.e. what is the value of the quantity in that position at a given time?).

Explore the Grid

The grid (or mesh) is one of the most used properties of the model, which can be accessed with the command:
disp(output.model.grid)
cells: [1×1 struct] + faces: [1×1 struct] + nodes: [1×1 struct] griddim: 1 - type: 'generic'
We can see that the grid is stored as a structure with information about the cells, faces, nodes, etc. The values of the state quantities (e.g. concentration and electric potential) are calculated at the centroids of the cells. To plot the positions of the centroids, we can use the following commands:
% get the centroid locations for the full grid
x = output.model.grid.cells.centroids;
 
% plot them in a figure
figure();
grid_axes = axes();
plot(grid_axes, x, zeros(size(x)), '+', 'LineWidth', 2)
 
% remove the y-axis ticks
yticks([]);
 
% add plot annotations
title('Grid Centroids');
xlabel('Position / m');
This shows the overall grid that is used for the model. However, BattMo models use a modular hierarchy where the overall cell model is composed of smaller submodels for electrodes, electrolyte, and current collectors. Each of these submodels has its own grid.
For example, if we want to plot the grid associated with the different submodels in different colors, we can use the following commands:
% get the centroid locations for the negative electrode (x_ne), separator
% (x_sep), and positive electrode (x_pe)
x_ne = output.model.NegativeElectrode.grid.cells.centroids;
x_sep = output.model.Separator.grid.cells.centroids;
x_pe = output.model.PositiveElectrode.grid.cells.centroids;
 
% plot them together on one figure
figure()
plot(x_ne, zeros(size(x_ne)), '+', 'LineWidth', 2)
hold on
plot(x_sep, zeros(size(x_sep)), '+k', 'LineWidth', 2)
plot(x_pe, zeros(size(x_pe)), '+r', 'LineWidth', 2)
 
% remove the y-axis tick marks
yticks([]);
 
% add plot annotations
title('Grid Centroids by Component')
xlabel('Position / m')
legend('Negative Electrode', 'Separator', 'Positive Electrode')
If you would like more information about the BattMo model hierarchy, please see BattMo Model Architecture.

Explore the States

The values of the state quantities at each time step are stored in the states cell array. Each entry in the array describes the state of the simulation at a given timestep.
For example, we can look at the state of the simulation at timestep 10 (shown in the dashboard plot above) using the command:
timestep = 10;
disp(output.states{timestep})
Electrolyte: [1x1 struct] - NegativeElectrode: [1x1 struct] - PositiveElectrode: [1x1 struct] - Control: [1x1 struct] - time: 504 - ThermalModel: [1x1 struct]
We see that the time of the state is 504 seconds and there are other structures containing the states of the electrodes and electrolyte. We can look into the state of the electrolyte using the command:
disp(output.states{timestep}.Electrolyte)
c: [30x1 double] - phi: [30x1 double]
which shows that there are two quantities there for concentration, c, and electric potential, phi.

Plot a Result

Let’s plot the concentration in the electrolyte at timestep 10. We can plot the results using basic MATLAB commands this way:
% get the centroid locations and the values for the electrolyte
% concentration at the given timestep
x = output.model.grid.cells.centroids;
c = output.states{timestep}.Electrolyte.c;
 
% plot the concentration values at the given grid centroid locations
figure()
plot(x, c, 'LineWidth', 3)
 
% add plot annotations
xlabel('Position / m')
ylabel('Concentration / mol \cdot m^{-3}')
BattMo also includes dedicated plotting functions that will come in handy when we start working with more complex systems (e.g. P4D grids). To demonstrate how it works, we can plot the electrolyte concentration profile at every timestep using the BattMo function plotCellData:
% create a figure and use hold on to show multiple lines in one figure
figure();
hold on
 
% define a colormap for the line colors
cm = cmocean('matter', length(output.states));
 
% iterate through each of the output states and plot the calculated
% electrolyte concentration at each timestep
for timestep = 1:length(output.states)
plotCellData(output.model.grid, output.states{timestep}.Electrolyte.c, 'LineWidth', 3, 'Color', cm(timestep,:))
end
 
% add plot annotations and turn hold off
xlabel('Position / m')
ylabel('Concentration / mol \cdot m^{-3}')
hold off
You did it! ? You have run an post-processed your first BattMo simulation! But there are still a lot of exciting features to discover. Let’s keep going and explore making chages to the model. ?
+ type: 'generic'
We can see that the grid is stored as a structure with information about the cells, faces, nodes, etc. The values of the state quantities (e.g. Li-ion concentration and electric potential) are calculated at the centroids of the cells. To plot the positions of the centroids, we can use the following commands:
% get the centroid locations for the full grid
x = output.model.grid.cells.centroids;
 
% plot them in a figure
figure();
grid_axes = axes();
plot(grid_axes, x, zeros(size(x)), '+', 'LineWidth', 2)
 
% remove the y-axis ticks
yticks([]);
 
% add plot annotations
title('Grid Centroids');
xlabel('Position / m');
This shows the overall grid that is used for the model. However, BattMo models use a modular hierarchy where the overall cell model is composed of smaller submodels for electrodes, electrolyte, and current collectors. Each of these submodels has its own grid.If you would like more information about the BattMo model hierarchy, please see BattMo Model Architecture.
For example, if we want to plot the grid associated with the different submodels in different colors, we can use the following commands:
% get the centroid locations for the negative electrode (x_ne), separator
% (x_sep), and positive electrode (x_pe)
x_ne = output.model.NegativeElectrode.grid.cells.centroids;
x_sep = output.model.Separator.grid.cells.centroids;
x_pe = output.model.PositiveElectrode.grid.cells.centroids;
 
% plot them together on one figure
figure()
hold on
plot(x_ne, zeros(size(x_ne)), '+', 'LineWidth', 2)
plot(x_sep, zeros(size(x_sep)), '+k', 'LineWidth', 2)
plot(x_pe, zeros(size(x_pe)), '+r', 'LineWidth', 2)
% remove the y-axis tick marks
yticks([]);
 
% add plot annotations
title('Grid Centroids by Component')
xlabel('Position / m')
legend('Negative Electrode', 'Separator', 'Positive Electrode')
Go ahead and add the centroids for the electrolyte to the above figure to see its grid distribution in the cell. The electrolyte overlaps with the electrodes throughout the cell.

Explore the States

The values of the state quantities at each time step are stored in the states cell array. Each entry in the array describes the state of the simulation at a given timestep.
For example, we can look at the state of the simulation at timestep 10 using the command:
timestep = 40;
disp(output.states{timestep})
Electrolyte: [1×1 struct] + NegativeElectrode: [1×1 struct] + PositiveElectrode: [1×1 struct] + Control: [1×1 struct] + time: 3.5692e+03 + ThermalModel: [1×1 struct]
We see that the time of the state is 3569 seconds and there are other structures containing the states of the electrodes and electrolyte. We can look into the state of the electrolyte using the command:
disp(output.states{timestep}.Electrolyte)
c: [30×1 double] + phi: [30×1 double]
which shows that there are two quantities there for concentration, c, and electric potential, phi.

Plot a Result

Let’s plot the concentration in the electrolyte at timestep 40. We can plot the results using basic MATLAB commands this way:
% get the centroid locations and the values for the electrolyte
% concentration at the given timestep
x = output.model.grid.cells.centroids;
c = output.states{timestep}.Electrolyte.c;
 
% plot the concentration values at the given grid centroid locations
figure()
plot(x, c, 'LineWidth', 3)
 
% add plot annotations
xlabel('Position / m')
ylabel('Concentration / mol \cdot m^{-3}')
BattMo also includes dedicated plotting functions that will come in handy when we start working with more complex systems (e.g. P4D grids). To demonstrate how it works, we can plot the electrolyte concentration profile at every timestep using the BattMo function plotCellData:
% create a figure and use hold on to show multiple lines in one figure
figure();
hold on
 
% define a colormap for the line colors
cm = cmocean('matter', length(output.states));
 
% iterate through each of the output states and plot the calculated
% electrolyte concentration at each timestep
for timestep = 1:length(output.states)
plotCellData(output.model.grid, output.states{timestep}.Electrolyte.c, 'LineWidth', 3, 'Color', cm(timestep,:))
end
 
% add plot annotations and turn hold off
xlabel('Position / m')
ylabel('Concentration / mol \cdot m^{-3}')
hold off
You did it! ? You have run and post-processed your first BattMo simulation! But there are still a lot of exciting features to discover. Let’s keep going and explore making changes to the model. ?

References

[1] Safari, M., et al. "Multimodal physics–based aging model for life prediction of Li–ion batteries." Journal of The Electrochemical Society 156.3 (2008): A145.
[2] Xu, Meng, et al. "A pseudo three–dimensional electrochemical–thermal model of a prismatic LiFePO4 battery during discharge process." Energy 80 (2015): 303–317.

\ No newline at end of file diff --git a/_static/notebooks/tutorial_2_changing_control_protocol_live.html b/_static/notebooks/tutorial_2_changing_control_protocol_live.html index 219bdfe6..a2076f3b 100644 --- a/_static/notebooks/tutorial_2_changing_control_protocol_live.html +++ b/_static/notebooks/tutorial_2_changing_control_protocol_live.html @@ -37,15 +37,15 @@ .S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } .S11 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } .S12 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S13 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

Tutorial 2 - Change the Control Protocol

Introduction

In this tutorial, we will use a P2D model to simulate the discharge of an NMC-Graphite cell at different rates. After completing this tutorial, you should have a working knowledge of:
  • Basics of the control protocol definitions in BattMo
  • The link between battery discharge rate and simulation time
  • How to setup and execute a parameter sweep
We'll use the same model from Tutorial 1.
jsonstruct = parseBattmoJson('Examples/jsondatafiles/sample_input.json');

Explore the Control Protocol

The control protocol is defined in the JSON parameter file and parsed into the MATLAB structure. Once the JSON parameter file has been read into MATLAB as a jsonstruct, its properties can be modified programmatically.
Let's begin by exploring the control protocol definition with the following command:
disp(jsonstruct.Control)
controlPolicy: 'CCDischarge' +.S13 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

Tutorial 2 - Change the Control Protocol

Introduction

In this tutorial, we will use a P2D model to simulate the discharge of an NMC-Graphite cell at different rates. After completing this tutorial, you should have a working knowledge of:
  • Basics of the control protocol definitions in BattMo
  • The link between battery discharge rate and simulation time
  • How to setup and execute a parameter sweep
We'll use the same input parameter set from Tutorial 1.
jsonstruct = parseBattmoJson('Examples/jsondatafiles/sample_input.json');

Explore the Control Protocol

The control protocol is defined in the JSON parameter file and parsed into the MATLAB structure. Once the JSON parameter file has been read into MATLAB as a jsonstruct, its properties can be modified programmatically.
Let's begin by exploring the control protocol definition with the following command:
disp(jsonstruct.Control)
controlPolicy: 'CCDischarge' CRate: 1 lowerCutoffVoltage: 2.4000 upperCutoffVoltage: 4.1000 dIdtLimit: 0.0100 dEdtLimit: 0.0100 - rampupTime: 0.1000
Here, we can see that the control policy follows the schema for a constant current discharge (CCDischarge). The three main control parameters in this schema are the CRate, lowerCutoffVoltage, and upperCutoffVoltage. Let's first try changing the protocol to discharge at a rate of C/10 to a lower cutoff voltage of 3 V.
jsonstruct.Control.CRate = 0.1;
jsonstruct.Control.lowerCutoffVoltage = 3.0;
Remember that if we change the rate of the discharge then we also need to change the duration of the simulation. To do this we can explore the TimeStepping property of the structure.
disp(jsonstruct.TimeStepping)
totalTime: 5040 + rampupTime: 0.1000
Here, we can see that the control policy follows the schema for a constant current discharge (CCDischarge). The three main control parameters in this schema are the CRate, lowerCutoffVoltage, and upperCutoffVoltage. Let's first try changing the protocol to discharge at a rate of C/10 to a lower cutoff voltage of 3 V. Read more about battery C-rates here
jsonstruct.Control.CRate = 0.1;
jsonstruct.Control.lowerCutoffVoltage = 3.0;
Remember that if we change the rate of the discharge then we also need to change the duration of the simulation. To do this we can explore the TimeStepping property of the structure.
disp(jsonstruct.TimeStepping)
totalTime: 5040 numberOfTimeSteps: 40 - useRampup: 1
The total duration of the simulation is given by the totalTime property. We can adjust it to an appropriate duration using the following command.
jsonstruct.TimeStepping.totalTime = (1./jsonstruct.Control.CRate) .* 3600 .* 1.1;
This sets the total time of the simulation to be 10% longer than the expected duration based on the C-Rate. We can then run the simulation and plot the discharge curve.
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 30 Seconds, 937 Milliseconds + useRampup: 1
The total duration of the simulation is given by the totalTime property. We can adjust it to an appropriate duration using the following command.
jsonstruct.TimeStepping.totalTime = (1./jsonstruct.Control.CRate) .* 3600 .* 1.1;
This sets the total time of the simulation to be 10% longer than the expected duration based on the C-Rate. We can then run the simulation using the same command before and plot the discharge curve.
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 30 Seconds, 937 Milliseconds Solving timestep 02/45: 30 Seconds, 937 Milliseconds -> 61 Seconds, 875 Milliseconds Solving timestep 03/45: 61 Seconds, 875 Milliseconds -> 123 Seconds, 750 Milliseconds Solving timestep 04/45: 123 Seconds, 750 Milliseconds -> 247 Seconds, 500 Milliseconds @@ -87,7 +87,7 @@ Solving timestep 40/45: 9 Hours, 1260 Seconds -> 9 Hours, 2250 Seconds Solving timestep 41/45: 9 Hours, 2250 Seconds -> 9 Hours, 3240 Seconds Solving timestep 42/45: 9 Hours, 3240 Seconds -> 10 Hours, 630 Seconds -*** Simulation complete. Solved 42 control steps in 5 Seconds, 966 Milliseconds (termination triggered by stopFunction) ***
retrieve the states from the simulation result
states = output.states;
 
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
 
% plot the discharge curve in a new figure
figure();
plot((time/hour), voltage, '-', 'linewidth', 3)
xlabel('Time / h')
ylabel('Cell Voltage / V')
title('Cell Discharge at C/10')

Setup and Run a Parameter Sweep

Now, let's setup and run a paramter sweep that simulates the performance of the cell at many different rates. To do this, we can define the different rates that we want to simulate in a vector and then use a for-loop to simulate the cell discharge at each value from the vector. We can also store the outputs of the various simulations as elements in a MATLAB cell array so that they are accessible at the end of the sweep.
% create a vector of different c-rates for the parameter sweep
CRates = [1, 2, 3];
 
% instantiate and empty cell array to store the outputs of the simulations
output = cell(size(CRates));
 
% instantiate an empty figure
figure()
 
% use a for-loop to iterate through the vector of c-rates
for i = 1 : numel(CRates)
% modify the value for the c-rate in the control definition and update
% the total duration of the simulation accordingly
jsonstruct.Control.CRate = CRates(i);
jsonstruct.TimeStepping.totalTime = (1./jsonstruct.Control.CRate) .* 3600 .* 1.2;
% run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct);
% retrieve the states from the simulation result
states = output{i}.states;
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
% plot the discharge curve in the figure
plot((time/hour), voltage, '-', 'linewidth', 3)
hold on
end
Solving timestep 01/45: -> 3 Seconds, 375 Milliseconds +*** Simulation complete. Solved 42 control steps in 2 Seconds, 18 Milliseconds (termination triggered by stopFunction) ***
Retrieve the states from the simulation result
states = output.states;
 
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
 
% plot the discharge curve in a new figure
figure();
plot((time/hour), voltage, '-', 'linewidth', 3)
xlabel('Time / h')
ylabel('Cell Voltage / V')
title('Cell Discharge at C/10')

Setup and Run a Parameter Sweep

Now, let's setup and run a paramter sweep that simulates the performance of the cell at many different rates. To do this, we can define the different rates that we want to simulate in a vector and then use a for-loop to simulate the cell discharge at each value from the vector. We can also store the outputs of the various simulations as elements in a MATLAB cell array so that they are accessible at the end of the sweep. We will keep the lower cut-off voltage at 3 V for all CRates.
% create a vector of different c-rates for the parameter sweep
CRates = [1, 2, 3];
 
% instantiate and empty cell array to store the outputs of the simulations
output = cell(size(CRates));
 
% instantiate an empty figure
figure()
 
% use a for-loop to iterate through the vector of c-rates
for i = 1 : numel(CRates)
% modify the value for the c-rate in the control definition and update
% the total duration of the simulation accordingly
jsonstruct.Control.CRate = CRates(i);
jsonstruct.TimeStepping.totalTime = (1./jsonstruct.Control.CRate) .* 3600 .* 1.2;
% run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct);
% retrieve the states from the simulation result
states = output{i}.states;
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
% plot the discharge curve in the figure
plot((time/hour), voltage, '-', 'linewidth', 3)
hold on
end
Solving timestep 01/45: -> 3 Seconds, 375 Milliseconds Solving timestep 02/45: 3 Seconds, 375 Milliseconds -> 6 Seconds, 750 Milliseconds Solving timestep 03/45: 6 Seconds, 750 Milliseconds -> 13 Seconds, 500 Milliseconds Solving timestep 04/45: 13 Seconds, 500 Milliseconds -> 27 Seconds @@ -126,7 +126,7 @@ Solving timestep 37/45: 3348 Seconds -> 3456 Seconds Solving timestep 38/45: 3456 Seconds -> 3564 Seconds Solving timestep 39/45: 3564 Seconds -> 1 Hour, 72 Seconds -*** Simulation complete. Solved 39 control steps in 5 Seconds, 314 Milliseconds (termination triggered by stopFunction) *** +*** Simulation complete. Solved 39 control steps in 1 Second, 778 Milliseconds (termination triggered by stopFunction) *** Solving timestep 01/45: -> 1 Second, 687 Milliseconds Solving timestep 02/45: 1 Second, 687 Milliseconds -> 3 Seconds, 375 Milliseconds Solving timestep 03/45: 3 Seconds, 375 Milliseconds -> 6 Seconds, 750 Milliseconds @@ -165,7 +165,7 @@ Solving timestep 36/45: 1620 Seconds -> 1674 Seconds Solving timestep 37/45: 1674 Seconds -> 1728 Seconds Solving timestep 38/45: 1728 Seconds -> 1782 Seconds -*** Simulation complete. Solved 38 control steps in 5 Seconds, 631 Milliseconds (termination triggered by stopFunction) *** +*** Simulation complete. Solved 38 control steps in 1 Second, 985 Milliseconds (termination triggered by stopFunction) *** Solving timestep 01/45: -> 1 Second, 125 Milliseconds Solving timestep 02/45: 1 Second, 125 Milliseconds -> 2 Seconds, 250 Milliseconds Solving timestep 03/45: 2 Seconds, 250 Milliseconds -> 4 Seconds, 500 Milliseconds @@ -201,7 +201,7 @@ Solving timestep 33/45: 972 Seconds -> 1008 Seconds Solving timestep 34/45: 1008 Seconds -> 1044 Seconds Solving timestep 35/45: 1044 Seconds -> 1080 Seconds -*** Simulation complete. Solved 35 control steps in 4 Seconds, 767 Milliseconds (termination triggered by stopFunction) ***
hold off
xlabel('Time / h')
ylabel('Voltage / V')
legend('1C', '2C', '3C')
Often, it is a more direct comparison of the discharge curves to plot against an x-axis of capacity rather than time. We can use the following commands to make such a plot.
% instantiate an empty figure
figure()
 
% use a for-loop to iterate through the vector of c-rates
for i = 1 : numel(output)
% modify the value for the c-rate in the control definition and update
% the total duration of the simulation accordingly
jsonstruct.Control.CRate = CRates(i);
jsonstruct.TimeStepping.totalTime = (1./jsonstruct.Control.CRate) .* 3600 .* 1.2;
% run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct);
% retrieve the states from the simulation result
states = output{i}.states;
% extract the time, voltage, and current quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
% calculate the capacity
capacity = time .* current;
% plot the discharge curve in the figure
plot((capacity/(hour*milli)), voltage, '-', 'linewidth', 3)
hold on
end
Solving timestep 01/45: -> 3 Seconds, 375 Milliseconds +*** Simulation complete. Solved 35 control steps in 1 Second, 807 Milliseconds (termination triggered by stopFunction) ***
hold off
xlabel('Time / h')
ylabel('Voltage / V')
legend('1C', '2C', '3C')
Often, it is a more direct comparison of the discharge curves to plot against an x-axis of capacity rather than time. The cell capacity is calculated as discharge time times the cell curent. We can use the following commands to make such a plot. We will
% instantiate an empty figure
figure()
 
% use a for-loop to iterate through the vector of c-rates
for i = 1 : numel(output)
% modify the value for the c-rate in the control definition and update
% the total duration of the simulation accordingly
jsonstruct.Control.CRate = CRates(i);
jsonstruct.TimeStepping.totalTime = (1./jsonstruct.Control.CRate) .* 3600 .* 1.2;
% run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct);
% retrieve the states from the simulation result
states = output{i}.states;
% extract the time, voltage, and current quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
% calculate the capacity
capacity = time .* current;
% plot the discharge curve in the figure
plot((capacity/(hour*milli)), voltage, '-', 'linewidth', 3)
hold on
end
Solving timestep 01/45: -> 3 Seconds, 375 Milliseconds Solving timestep 02/45: 3 Seconds, 375 Milliseconds -> 6 Seconds, 750 Milliseconds Solving timestep 03/45: 6 Seconds, 750 Milliseconds -> 13 Seconds, 500 Milliseconds Solving timestep 04/45: 13 Seconds, 500 Milliseconds -> 27 Seconds @@ -240,7 +240,7 @@ Solving timestep 37/45: 3348 Seconds -> 3456 Seconds Solving timestep 38/45: 3456 Seconds -> 3564 Seconds Solving timestep 39/45: 3564 Seconds -> 1 Hour, 72 Seconds -*** Simulation complete. Solved 39 control steps in 4 Seconds, 694 Milliseconds (termination triggered by stopFunction) *** +*** Simulation complete. Solved 39 control steps in 1 Second, 782 Milliseconds (termination triggered by stopFunction) *** Solving timestep 01/45: -> 1 Second, 687 Milliseconds Solving timestep 02/45: 1 Second, 687 Milliseconds -> 3 Seconds, 375 Milliseconds Solving timestep 03/45: 3 Seconds, 375 Milliseconds -> 6 Seconds, 750 Milliseconds @@ -279,7 +279,7 @@ Solving timestep 36/45: 1620 Seconds -> 1674 Seconds Solving timestep 37/45: 1674 Seconds -> 1728 Seconds Solving timestep 38/45: 1728 Seconds -> 1782 Seconds -*** Simulation complete. Solved 38 control steps in 5 Seconds, 519 Milliseconds (termination triggered by stopFunction) *** +*** Simulation complete. Solved 38 control steps in 1 Second, 924 Milliseconds (termination triggered by stopFunction) *** Solving timestep 01/45: -> 1 Second, 125 Milliseconds Solving timestep 02/45: 1 Second, 125 Milliseconds -> 2 Seconds, 250 Milliseconds Solving timestep 03/45: 2 Seconds, 250 Milliseconds -> 4 Seconds, 500 Milliseconds @@ -315,7 +315,7 @@ Solving timestep 33/45: 972 Seconds -> 1008 Seconds Solving timestep 34/45: 1008 Seconds -> 1044 Seconds Solving timestep 35/45: 1044 Seconds -> 1080 Seconds -*** Simulation complete. Solved 35 control steps in 5 Seconds, 390 Milliseconds (termination triggered by stopFunction) ***
hold off
xlabel('Capacity / mA \cdot h')
ylabel('Voltage / V')
legend('1C', '2C', '3C')

Summary

In this tutorial, we explored the basic structure of the control protocol definition and made changes to simulate cell discharge at a variety of rates. We showed that the control protocol is defined as part of the overall BattMo parameter structure. It can be modified by changing the properties of the Control structure. We learned that the total duration of the simulation is dependent on the rate of the battery discharge, and should be updated accordingly. Finally, we showed how to setup and execute a basic parameter sweep by defining a vector with different C-rates and using a for-loop to simulate the cell discharge at each rate.
+*** Simulation complete. Solved 35 control steps in 1 Second, 839 Milliseconds (termination triggered by stopFunction) ***
hold off
xlabel('Capacity / mA \cdot h')
ylabel('Voltage / V')
legend('1C', '2C', '3C')

Summary

In this tutorial, we explored the basic structure of the control protocol definition and made changes to simulate cell discharge at a variety of rates. We showed that the control protocol is defined as part of the overall BattMo parameter structure. It can be modified by changing the properties of the Control structure. We learned that the total duration of the simulation is dependent on the rate of the battery discharge, and should be updated accordingly. Finally, we showed how to setup and execute a basic parameter sweep by defining a vector with different C-rates and using a for-loop to simulate the cell discharge at each rate.

\ No newline at end of file diff --git a/_static/notebooks/tutorial_4_modify_material_parameters_live.html b/_static/notebooks/tutorial_4_modify_material_parameters_live.html index bd83d349..7b9bf41b 100644 --- a/_static/notebooks/tutorial_4_modify_material_parameters_live.html +++ b/_static/notebooks/tutorial_4_modify_material_parameters_live.html @@ -32,19 +32,20 @@ .inlineElement .textElement {} .embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement { min-height: 16px;} .rightPaneElement .textElement { padding-top: 2px; padding-left: 9px;} -.S8 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S9 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S11 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S12 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: 700; text-align: left; } -.S13 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

Tutorial 4 - Modify Material Parameters

Introduction

In this tutorial, we will use a P2D model to simulate the effects of changing material parameters on cell discharge. After completing this tutorial, you should have a working knowledge of:
  • Basics of how material properties are defined in BattMo
  • How to make some simple changes to material property definitions
  • How to change material models from the BattMo material library
We'll use the same model from Tutorial 1.
jsonstruct = parseBattmoJson('Examples/jsondatafiles/sample_input.json');
Parameters are defined in the JSON parameter file and parsed into the MATLAB structure. Once the JSON file has been read into MATLAB as a jsonstruct, its properties can be modified programmatically.

Explore the Material Parameters

Let's begin by exploring the active material of the positive electrode coating with the following command:
disp(jsonstruct.PositiveElectrode.Coating.ActiveMaterial)
massFraction: 0.9500 +.S8 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: 400; text-align: left; } +.S9 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S11 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S12 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S13 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S14 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: 700; text-align: left; }

Tutorial 4 - Modify Material Parameters

Introduction

In this tutorial, we will use a P2D model to simulate the effects of changing material parameters on cell discharge. After completing this tutorial, you should have a working knowledge of:
  • Basics of how material properties are defined in BattMo
  • How to make some simple changes to material property definitions
  • How to change material models from the BattMo material library
We'll use the same input parameter set from Tutorial 1.
jsonstruct = parseBattmoJson('Examples/jsondatafiles/sample_input.json');
Parameters are defined in the JSON parameter file and parsed into the MATLAB structure. Once the JSON file has been read into MATLAB as a jsonstruct, its properties can be modified programmatically.

Explore the Material Parameters

Let's begin by exploring the active material of the positive electrode (NMC) coating with the following command:
disp(jsonstruct.PositiveElectrode.Coating.ActiveMaterial)
massFraction: 0.9500 density: 4650 specificHeatCapacity: 700 thermalConductivity: 2.1000 electronicConductivity: 100 - Interface: [1x1 struct] + Interface: [1×1 struct] diffusionModelType: 'full' - SolidDiffusion: [1x1 struct]
We can see that a variety of parameters are defined that determine the physicochemical properties of the material. Many of the parameters that describe the electrochemical reactions that take place at the active material - electrolyte interface are contained in the Interface structure, which we can explore with the command:
disp(jsonstruct.PositiveElectrode.Coating.ActiveMaterial.Interface)
saturationConcentration: 55554 + SolidDiffusion: [1×1 struct]
We can see that a variety of parameters are defined that determine the physicochemical properties of the material. Many of the parameters that describe the electrochemical reactions that take place at the active material - electrolyte interface are contained in the Interface structure, which we can explore with the command:
disp(jsonstruct.PositiveElectrode.Coating.ActiveMaterial.Interface)
saturationConcentration: 55554 volumetricSurfaceArea: 885000 density: 4650 numberOfElectronsTransferred: 1 @@ -53,10 +54,10 @@ guestStoichiometry100: 0.4955 guestStoichiometry0: 0.9917 chargeTransferCoefficient: 0.5000 - openCircuitPotential: [1x1 struct]
Furthermore, it is well known that the diffusion of lithium in the active materials can be a limiting factor in Li-ion battery performance. The diffusion model used in this simulation can be explored with the following command:
disp(jsonstruct.PositiveElectrode.Coating.ActiveMaterial.SolidDiffusion)
activationEnergyOfDiffusion: 5000 + openCircuitPotential: [1×1 struct]
Furthermore, it is well known that the diffusion of lithium in the active materials can be a limiting factor in Li-ion battery performance. The diffusion model used in this simulation can be explored with the following command:
disp(jsonstruct.PositiveElectrode.Coating.ActiveMaterial.SolidDiffusion)
activationEnergyOfDiffusion: 5000 referenceDiffusionCoefficient: 1.0000e-14 particleRadius: 1.0000e-06 - N: 10
We can see there that the solid diffusion model includes a reference diffusion coefficient, activation energy of diffusion, and the effective particle radius.

Make a Simple Change

Let's make a simple change to the positive electrode active material properties in our model and see how it affects the results. Let's compare how c
% create a vector of diffent thickness values
radius = [1, 5, 10].*micro;
 
% instantiate and empty cell array to store the outputs of the simulations
output = cell(size(radius));
 
% instantiate an empty figure
figure()
 
% use a for-loop to iterate through the vector of c-rates
for i = 1 : numel(radius)
% modify the value for the c-rate in the control definition and update
% the total duration of the simulation accordingly
jsonstruct.PositiveElectrode.Coating.ActiveMaterial.SolidDiffusion.particleRadius = radius(i);
% run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct);
% retrieve the states from the simulation result
states = output{i}.states;
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
% calculate the capacity
capacity = time .* current;
% plot the discharge curve in the figure
plot((capacity/(milli*hour)), voltage, '-', 'linewidth', 3)
hold on
end
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds + N: 10
We can see there that the solid diffusion model includes a reference diffusion coefficient, activation energy of diffusion, and the effective particle radius.

Make a Simple Change

Let's make a simple change to the positive electrode active material properties in our model and see how it affects the results. Let's compare how cell capacity is affected by changes in the active particle radius. With a larger particle radius and the same amount of active material, the Li ions have a longer path to traverse, hence one would expect the capacity to actually drop with increase in the radius. Let's check that with our parameter sweep.
% create a vector of diffent thickness values
radius = [1, 5, 10].*micro;
 
% instantiate and empty cell array to store the outputs of the simulations
output = cell(size(radius));
 
% instantiate an empty figure
figure()
 
% use a for-loop to iterate through the vector of c-rates
for i = 1 : numel(radius)
% modify the value for the radius in the active material definition
jsonstruct.PositiveElectrode.Coating.ActiveMaterial.SolidDiffusion.particleRadius = radius(i);
% run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct);
% retrieve the states from the simulation result
states = output{i}.states;
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
% calculate the capacity
capacity = time .* current;
% plot the discharge curve in the figure
plot((capacity/(milli*hour)), voltage, '-', 'linewidth', 3)
hold on
end
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds @@ -90,7 +91,7 @@ Solving timestep 32/45: 3276 Seconds -> 3402 Seconds Solving timestep 33/45: 3402 Seconds -> 3528 Seconds Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds -*** Simulation complete. Solved 34 control steps in 7 Seconds, 380 Milliseconds (termination triggered by stopFunction) *** +*** Simulation complete. Solved 34 control steps in 2 Seconds, 319 Milliseconds (termination triggered by stopFunction) *** Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds @@ -124,7 +125,7 @@ Solving timestep 31/45: 3150 Seconds -> 3276 Seconds Solving timestep 32/45: 3276 Seconds -> 3402 Seconds Solving timestep 33/45: 3402 Seconds -> 3528 Seconds -*** Simulation complete. Solved 33 control steps in 6 Seconds, 954 Milliseconds (termination triggered by stopFunction) *** +*** Simulation complete. Solved 33 control steps in 2 Seconds, 313 Milliseconds (termination triggered by stopFunction) *** Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds @@ -154,14 +155,14 @@ Solving timestep 27/45: 2646 Seconds -> 2772 Seconds Solving timestep 28/45: 2772 Seconds -> 2898 Seconds Solving timestep 29/45: 2898 Seconds -> 3024 Seconds -*** Simulation complete. Solved 29 control steps in 7 Seconds, 773 Milliseconds (termination triggered by stopFunction) ***
hold off
 
% add plot annotations
xlabel('Capacity / mA \cdot h')
ylabel('Cell Voltage / V')
 
% we can pull the radius values directly from the vector to make the plot
% annotations more resilient to changes in the code
legend(strcat('r=', num2str(radius(1)/micro), ' µm'), strcat('r=', num2str(radius(2)/micro), ' µm'), strcat('r=', num2str(radius(3)/micro), ' µm'))

Swap Active Materials

Now let's try to replace the active material with a different material from the BattMo library. In this example, we will replace the NMC material with LFP. First, let's re-load our baseline example parameter set to get rid of the changes from the previous sections.
clear
jsonstruct = parseBattmoJson('Examples/jsondatafiles/sample_input.json');
Now we load and parse the LFP material parameters from the BattMo library and move it to the right place in the model hierarchy:
lfp = parseBattmoJson('ParameterData/MaterialProperties/LFP/LFP.json');
jsonstruct_lfp.PositiveElectrode.Coating.ActiveMaterial.Interface = lfp;
To merge new parameter data into our existing model, we can use the BattMo function mergeJsonStructs.
jsonstruct = mergeJsonStructs({jsonstruct_lfp, ...
jsonstruct});
parameter PositiveElectrode.Coating.ActiveMaterial.Interface.saturationConcentration is assigned twice with different values. we use the value from first jsonstruct. +*** Simulation complete. Solved 29 control steps in 2 Seconds, 606 Milliseconds (termination triggered by stopFunction) ***
hold off
% add plot annotations
xlabel('Capacity / mA \cdot h')
ylabel('Cell Voltage / V')
% we can pull the radius values directly from the vector to make the plot
% annotations more resilient to changes in the code
legend(strcat('r=', num2str(radius(1)/micro), ' µm'), strcat('r=', num2str(radius(2)/micro), ' µm'), strcat('r=', num2str(radius(3)/micro), ' µm'))

Swap Active Materials

Now let's try to replace the active material with a different material from the BattMo library. In this example, we will replace the NMC material with LFP [1]. First, let's re-load our baseline example parameter set to get rid of the changes from the previous sections.
clear
jsonstruct = parseBattmoJson('Examples/jsondatafiles/sample_input.json');
Now we load and parse the LFP material parameters from the BattMo library and move it to the right place in the model hierarchy:
lfp = parseBattmoJson('ParameterData/MaterialProperties/LFP/LFP.json');
jsonstruct_lfp.PositiveElectrode.Coating.ActiveMaterial.Interface = lfp;
 
To merge new parameter data into our existing input model, we can use the BattMo function mergeJsonStructs.
jsonstruct = mergeJsonStructs({jsonstruct_lfp, ...
jsonstruct});
parameter PositiveElectrode.Coating.ActiveMaterial.Interface.saturationConcentration is assigned twice with different values. we use the value from first jsonstruct. parameter PositiveElectrode.Coating.ActiveMaterial.Interface.volumetricSurfaceArea is assigned twice with different values. we use the value from first jsonstruct. parameter PositiveElectrode.Coating.ActiveMaterial.Interface.density is assigned twice with different values. we use the value from first jsonstruct. parameter PositiveElectrode.Coating.ActiveMaterial.Interface.activationEnergyOfReaction is assigned twice with different values. we use the value from first jsonstruct. parameter PositiveElectrode.Coating.ActiveMaterial.Interface.reactionRateConstant is assigned twice with different values. we use the value from first jsonstruct. parameter PositiveElectrode.Coating.ActiveMaterial.Interface.guestStoichiometry100 is assigned twice with different values. we use the value from first jsonstruct. parameter PositiveElectrode.Coating.ActiveMaterial.Interface.guestStoichiometry0 is assigned twice with different values. we use the value from first jsonstruct. -parameter PositiveElectrode.Coating.ActiveMaterial.Interface.openCircuitPotential.functionname is assigned twice with different values. we use the value from first jsonstruct.
We need to be sure that the parameters are consistent across the hierarchy:
jsonstruct.PositiveElectrode.Coating.ActiveMaterial.density = jsonstruct.PositiveElectrode.Coating.ActiveMaterial.Interface.density;
jsonstruct.PositiveElectrode.Coating.effectiveDensity = 900;
And now, we can run the simulation and plot the discharge curve:
% run the simulation
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds +parameter PositiveElectrode.Coating.ActiveMaterial.Interface.openCircuitPotential.functionname is assigned twice with different values. we use the value from first jsonstruct.
We need to be sure that the parameters are consistent across the hierarchy: The effective density (the mass density of the material (either in wet or calendared state). This density is computed with respect to the total volume (i.e. including the empty pores)) of this material, assuming a porosity of 0.2 is calculated to be 2807 kg/m3.
disp(jsonstruct.PositiveElectrode.Coating.ActiveMaterial.Interface.density);
3600
jsonstruct.PositiveElectrode.Coating.ActiveMaterial.density = jsonstruct.PositiveElectrode.Coating.ActiveMaterial.Interface.density;
jsonstruct.PositiveElectrode.Coating.effectiveDensity = 2807;
And now, we can run the simulation and plot the discharge curve:
% run the simulation
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds @@ -192,10 +193,11 @@ Solving timestep 29/45: 2898 Seconds -> 3024 Seconds Solving timestep 30/45: 3024 Seconds -> 3150 Seconds Solving timestep 31/45: 3150 Seconds -> 3276 Seconds +Solver did not converge in 10 iterations for timestep of length 126 Seconds. Cutting timestep. Solving timestep 32/45: 3276 Seconds -> 3402 Seconds Solving timestep 33/45: 3402 Seconds -> 3528 Seconds Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds -*** Simulation complete. Solved 34 control steps in 6 Seconds, 809 Milliseconds (termination triggered by stopFunction) ***
get the states
states = output.states;
 
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
 
% calculate the capacity
capacity = time .* current;
 
% plot the discharge curve in the figure
plot((capacity/(milli*hour)), voltage, '-', 'linewidth', 3)
 
% add plot annotations
xlabel('Capacity / mA \cdot h')
ylabel('Cell Voltage / V')

Summary

In this tutorial, we explored how to modify material parameters in BattMo. We first explored the material parameter structure. Then we modified the effective particle radius and simulated its affect on the cell discharge curve. Finally we swapped out the NMC active material for a LFP material and simulated the new performance of the cell.
+*** Simulation complete. Solved 34 control steps in 2 Seconds, 274 Milliseconds (termination triggered by stopFunction) ***
We now get the states,
states = output.states;
 
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
 
% calculate the capacity
capacity = time .* current;
 
% plot the discharge curve in the figure
plot((capacity/(milli*hour)), voltage, '-', 'linewidth', 3)
 
% add plot annotations
xlabel('Capacity / mA \cdot h')
ylabel('Cell Voltage / V')

Summary

In this tutorial, we explored how to modify material parameters in BattMo. We first explored the material parameter structure. Then we modified the effective particle radius and simulated its affect on the cell discharge curve. Finally we swapped out the NMC active material for a LFP material and simulated the new performance of the cell, we saw that the LFP-graphite cell has reduced capacity compared to the NMC cell.

References

[1] Xu, Meng, et al. "A pseudo three–dimensional electrochemical–thermal model of a prismatic LiFePO4 battery during discharge process." Energy 80 (2015): 303–317.

\ No newline at end of file diff --git a/_static/notebooks/tutorial_5_simulate_CCCV_cycling_live.html b/_static/notebooks/tutorial_5_simulate_CCCV_cycling_live.html index 17743096..86c3b84e 100644 --- a/_static/notebooks/tutorial_5_simulate_CCCV_cycling_live.html +++ b/_static/notebooks/tutorial_5_simulate_CCCV_cycling_live.html @@ -34,22 +34,24 @@ .rightPaneElement .textElement { padding-top: 2px; padding-left: 9px;} .S8 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } .S9 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S11 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: 700; text-align: left; }

Tutorial 5 - Simulate CC-CV Cycling

Introduction

In this tutorial, we will use a P2D model to simulate CC-CV cycling. After completing this tutorial, you should have a working knowledge of:
  • How to define and modify cycling protocols in BattMo
We'll use the same model from Tutorial 1.
jsonstruct = parseBattmoJson('Examples/jsondatafiles/sample_input.json');
Parameters are defined in the JSON parameter file and parsed into the MATLAB structure. Once the JSON file has been read into MATLAB as a jsonstruct, its properties can be modified programmatically.

Explore the Control Definition

Let's begin by reviewing the control protocol in BattMo, with the command:
disp(jsonstruct.Control)
controlPolicy: 'CCDischarge' +.S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S11 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S12 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S13 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: 400; text-align: left; }

Tutorial 5 - Simulate CC-CV Cycling

Introduction

In this tutorial, we will use a P2D model to simulate CC-CV cycling. After completing this tutorial, you should have a working knowledge of:
  • How to define and modify cycling protocols in BattMo
In the CC-CV cycling, the cell is discharged at constant current (CC) and then the charging is also done at constant current till we get close to the upper cut-off voltage and then it is switched to constant voltage mode, where the cell potential is kept constant and the charging current falls off.
We'll use the same model from Tutorial 1.
jsonstruct = parseBattmoJson('Examples/jsondatafiles/sample_input.json');
Parameters are defined in the JSON parameter file and parsed into the MATLAB structure. Once the JSON file has been read into MATLAB as a jsonstruct, its properties can be modified programmatically.

Explore the Control Definition

Let's begin by reviewing the control protocol in BattMo, with the command:
disp(jsonstruct.Control)
controlPolicy: 'CCDischarge' CRate: 1 lowerCutoffVoltage: 2.4000 upperCutoffVoltage: 4.1000 dIdtLimit: 0.0100 dEdtLimit: 0.0100 - rampupTime: 0.1000
We see that the default control protocol is set to a constant current (galvanostatic) discharge. To change to a CC-CV cycling protocol, we can use the command:
cccv_control_protocol = parseBattmoJson('cccv_control.json');
jsonstruct_modified = mergeJsonStructs({cccv_control_protocol, jsonstruct});
parameter Control.controlPolicy is assigned twice with different values. we use the value from first jsonstruct. -parameter Control.lowerCutoffVoltage is assigned twice with different values. we use the value from first jsonstruct.
Now we can explore the modified control protocol definition with the command:
disp(jsonstruct_modified.Control)
controlPolicy: 'CCCV' + rampupTime: 0.1000
We see that the default control protocol is set to a constant current (galvanostatic) discharge. To change to a CC-CV cycling protocol, we can use the command:
cccv_control_protocol = parseBattmoJson('cccv_control.json');
jsonstruct_modified = mergeJsonStructs({cccv_control_protocol, jsonstruct});
parameter Control.controlPolicy is assigned twice with different values. we use the value from first jsonstruct. +parameter Control.lowerCutoffVoltage is assigned twice with different values. we use the value from first jsonstruct.
Now we can explore the modified control protocol definition with the command:
disp(jsonstruct_modified.Control)
controlPolicy: 'CCCV' initialControl: 'discharging' CRate: 1 lowerCutoffVoltage: 2.6000 upperCutoffVoltage: 4.1000 dIdtLimit: 0.0100 dEdtLimit: 0.0100 - rampupTime: 0.1000
Let's run the simulation and plot the cell voltage curve.
% run the simulation
output = runBatteryJson(jsonstruct_modified);
Warning: Number of cycles has not been given in CCCV control. We use numberOfCycles = 1.
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds + rampupTime: 0.1000
 
jsonstruct_modified.NegativeElectrode.Coating.thickness = 32.*micro;
Let's run the simulation and plot the cell voltage curve.
% run the simulation
output_01 = runBatteryJson(jsonstruct_modified);
Warning: Number of cycles has not been given in CCCV control. We use numberOfCycles = 1.
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds Solving timestep 02/45: 6 Seconds, 187 Milliseconds -> 12 Seconds, 375 Milliseconds Solving timestep 03/45: 12 Seconds, 375 Milliseconds -> 24 Seconds, 750 Milliseconds Solving timestep 04/45: 24 Seconds, 750 Milliseconds -> 49 Seconds, 500 Milliseconds @@ -73,11 +75,106 @@ Solving timestep 22/45: 3168 Seconds -> 3366 Seconds Solving timestep 23/45: 3366 Seconds -> 3564 Seconds Solving timestep 24/45: 3564 Seconds -> 1 Hour, 162 Seconds +Solving timestep 25/45: 1 Hour, 162 Seconds -> 1 Hour, 360 Seconds +Solving timestep 26/45: 1 Hour, 360 Seconds -> 1 Hour, 558 Seconds +Solving timestep 27/45: 1 Hour, 558 Seconds -> 1 Hour, 756 Seconds +switch control from CC_discharge1 to CC_discharge2 +Switch control type from CC_discharge2 to CC_charge1 +Solving timestep 28/45: 1 Hour, 756 Seconds -> 1 Hour, 954 Seconds +Solving timestep 29/45: 1 Hour, 954 Seconds -> 1 Hour, 1152 Seconds +Solving timestep 30/45: 1 Hour, 1152 Seconds -> 1 Hour, 1350 Seconds +Solving timestep 31/45: 1 Hour, 1350 Seconds -> 1 Hour, 1548 Seconds +Solving timestep 32/45: 1 Hour, 1548 Seconds -> 1 Hour, 1746 Seconds +Solving timestep 33/45: 1 Hour, 1746 Seconds -> 1 Hour, 1944 Seconds +Solving timestep 34/45: 1 Hour, 1944 Seconds -> 1 Hour, 2142 Seconds +Solving timestep 35/45: 1 Hour, 2142 Seconds -> 1 Hour, 2340 Seconds +Solving timestep 36/45: 1 Hour, 2340 Seconds -> 1 Hour, 2538 Seconds +Solving timestep 37/45: 1 Hour, 2538 Seconds -> 1 Hour, 2736 Seconds +Solving timestep 38/45: 1 Hour, 2736 Seconds -> 1 Hour, 2934 Seconds +Solving timestep 39/45: 1 Hour, 2934 Seconds -> 1 Hour, 3132 Seconds +Solving timestep 40/45: 1 Hour, 3132 Seconds -> 1 Hour, 3330 Seconds +Solving timestep 41/45: 1 Hour, 3330 Seconds -> 1 Hour, 3528 Seconds +Solving timestep 42/45: 1 Hour, 3528 Seconds -> 2 Hours, 126 Seconds +Solving timestep 43/45: 2 Hours, 126 Seconds -> 2 Hours, 324 Seconds +Solving timestep 44/45: 2 Hours, 324 Seconds -> 2 Hours, 522 Seconds +Solving timestep 45/45: 2 Hours, 522 Seconds -> 2 Hours, 720 Seconds +Switch control type from CV_charge2 to CC_discharge1 +*** Simulation complete. Solved 45 control steps in 4 Seconds, 497 Milliseconds (termination triggered by stopFunction) ***
get the states
states_01 = output_01.states;
% instantiate an empty figure
figure()
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states_01);
voltage = cellfun(@(state) state.('Control').E, states_01);
current = cellfun(@(state) state.('Control').I, states_01);
 
% calculate the capacity
capacity = time .* current;
index=find(capacity == 0);
% plot the discharge curve in the figure
 
plot(time/(hour), voltage, '-', 'linewidth', 3)
 
% add plot annotations
xlabel('Time / h')
ylabel('Cell Voltage / V')
Now let us create electrodes of different thickeness' and run through the CC-CV cycling to see how thickness affects the cell's capacity.
% create a vector of diffent thickness values
thickness = [16,32,64,128].*micro;
markers={'-','-','-','o'};
% instantiate and empty cell array to store the outputs of the simulations
output = cell(size(thickness));
 
% instantiate an empty figure
figure()
 
% use a for-loop to iterate through the vector of c-rates
for i = 1 : numel(thickness)
% modify the value for the coating thickness for the negative electrode
jsonstruct_modified.NegativeElectrode.Coating.thickness = thickness(i);
% run the simulation and store the results in the output cell array
output{i} = runBatteryJson(jsonstruct_modified);
% retrieve the states from the simulation result
states = output{i}.states;
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
% calculate the capacity
capacity = time .* current;
index=find(capacity == 0);
% plot only the discharge part of the curve in the figure
plot((capacity(1:index-1)/(hour*milli)), voltage(1:index-1), markers{i}, 'linewidth', 3)
hold on
end
Warning: Number of cycles has not been given in CCCV control. We use numberOfCycles = 1.
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds +Solving timestep 02/45: 6 Seconds, 187 Milliseconds -> 12 Seconds, 375 Milliseconds +Solving timestep 03/45: 12 Seconds, 375 Milliseconds -> 24 Seconds, 750 Milliseconds +Solving timestep 04/45: 24 Seconds, 750 Milliseconds -> 49 Seconds, 500 Milliseconds +Solving timestep 05/45: 49 Seconds, 500 Milliseconds -> 99 Seconds +Solving timestep 06/45: 99 Seconds -> 198 Seconds +Solving timestep 07/45: 198 Seconds -> 396 Seconds +Solving timestep 08/45: 396 Seconds -> 594 Seconds +Solving timestep 09/45: 594 Seconds -> 792 Seconds +Solving timestep 10/45: 792 Seconds -> 990 Seconds +Solving timestep 11/45: 990 Seconds -> 1188 Seconds +Solving timestep 12/45: 1188 Seconds -> 1386 Seconds +Solving timestep 13/45: 1386 Seconds -> 1584 Seconds +Solving timestep 14/45: 1584 Seconds -> 1782 Seconds +Solving timestep 15/45: 1782 Seconds -> 1980 Seconds +Solving timestep 16/45: 1980 Seconds -> 2178 Seconds +Solving timestep 17/45: 2178 Seconds -> 2376 Seconds +Solving timestep 18/45: 2376 Seconds -> 2574 Seconds +Solving timestep 19/45: 2574 Seconds -> 2772 Seconds +Solving timestep 20/45: 2772 Seconds -> 2970 Seconds +Solving timestep 21/45: 2970 Seconds -> 3168 Seconds +Solving timestep 22/45: 3168 Seconds -> 3366 Seconds +Solving timestep 23/45: 3366 Seconds -> 3564 Seconds +Solving timestep 24/45: 3564 Seconds -> 1 Hour, 162 Seconds +Solving timestep 25/45: 1 Hour, 162 Seconds -> 1 Hour, 360 Seconds +Solving timestep 26/45: 1 Hour, 360 Seconds -> 1 Hour, 558 Seconds +Solving timestep 27/45: 1 Hour, 558 Seconds -> 1 Hour, 756 Seconds switch control from CC_discharge1 to CC_discharge2 Switch control type from CC_discharge2 to CC_charge1 +Solving timestep 28/45: 1 Hour, 756 Seconds -> 1 Hour, 954 Seconds +Solving timestep 29/45: 1 Hour, 954 Seconds -> 1 Hour, 1152 Seconds +Solving timestep 30/45: 1 Hour, 1152 Seconds -> 1 Hour, 1350 Seconds +Solving timestep 31/45: 1 Hour, 1350 Seconds -> 1 Hour, 1548 Seconds +Solving timestep 32/45: 1 Hour, 1548 Seconds -> 1 Hour, 1746 Seconds +Solving timestep 33/45: 1 Hour, 1746 Seconds -> 1 Hour, 1944 Seconds +Solving timestep 34/45: 1 Hour, 1944 Seconds -> 1 Hour, 2142 Seconds +Solving timestep 35/45: 1 Hour, 2142 Seconds -> 1 Hour, 2340 Seconds +Solving timestep 36/45: 1 Hour, 2340 Seconds -> 1 Hour, 2538 Seconds +Solving timestep 37/45: 1 Hour, 2538 Seconds -> 1 Hour, 2736 Seconds +Solving timestep 38/45: 1 Hour, 2736 Seconds -> 1 Hour, 2934 Seconds +Solving timestep 39/45: 1 Hour, 2934 Seconds -> 1 Hour, 3132 Seconds +Solving timestep 40/45: 1 Hour, 3132 Seconds -> 1 Hour, 3330 Seconds +Solving timestep 41/45: 1 Hour, 3330 Seconds -> 1 Hour, 3528 Seconds +Solving timestep 42/45: 1 Hour, 3528 Seconds -> 2 Hours, 126 Seconds +Solving timestep 43/45: 2 Hours, 126 Seconds -> 2 Hours, 324 Seconds +Solving timestep 44/45: 2 Hours, 324 Seconds -> 2 Hours, 522 Seconds +Switch control type from CV_charge2 to CC_discharge1 +*** Simulation complete. Solved 44 control steps in 4 Seconds, 765 Milliseconds (termination triggered by stopFunction) ***
Warning: Number of cycles has not been given in CCCV control. We use numberOfCycles = 1.
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds +Solving timestep 02/45: 6 Seconds, 187 Milliseconds -> 12 Seconds, 375 Milliseconds +Solving timestep 03/45: 12 Seconds, 375 Milliseconds -> 24 Seconds, 750 Milliseconds +Solving timestep 04/45: 24 Seconds, 750 Milliseconds -> 49 Seconds, 500 Milliseconds +Solving timestep 05/45: 49 Seconds, 500 Milliseconds -> 99 Seconds +Solving timestep 06/45: 99 Seconds -> 198 Seconds +Solving timestep 07/45: 198 Seconds -> 396 Seconds +Solving timestep 08/45: 396 Seconds -> 594 Seconds +Solving timestep 09/45: 594 Seconds -> 792 Seconds +Solving timestep 10/45: 792 Seconds -> 990 Seconds +Solving timestep 11/45: 990 Seconds -> 1188 Seconds +Solving timestep 12/45: 1188 Seconds -> 1386 Seconds +Solving timestep 13/45: 1386 Seconds -> 1584 Seconds +Solving timestep 14/45: 1584 Seconds -> 1782 Seconds +Solving timestep 15/45: 1782 Seconds -> 1980 Seconds +Solving timestep 16/45: 1980 Seconds -> 2178 Seconds +Solving timestep 17/45: 2178 Seconds -> 2376 Seconds +Solving timestep 18/45: 2376 Seconds -> 2574 Seconds +Solving timestep 19/45: 2574 Seconds -> 2772 Seconds +Solving timestep 20/45: 2772 Seconds -> 2970 Seconds +Solving timestep 21/45: 2970 Seconds -> 3168 Seconds +Solving timestep 22/45: 3168 Seconds -> 3366 Seconds +Solving timestep 23/45: 3366 Seconds -> 3564 Seconds +Solving timestep 24/45: 3564 Seconds -> 1 Hour, 162 Seconds Solving timestep 25/45: 1 Hour, 162 Seconds -> 1 Hour, 360 Seconds Solving timestep 26/45: 1 Hour, 360 Seconds -> 1 Hour, 558 Seconds Solving timestep 27/45: 1 Hour, 558 Seconds -> 1 Hour, 756 Seconds +switch control from CC_discharge1 to CC_discharge2 +Switch control type from CC_discharge2 to CC_charge1 Solving timestep 28/45: 1 Hour, 756 Seconds -> 1 Hour, 954 Seconds Solving timestep 29/45: 1 Hour, 954 Seconds -> 1 Hour, 1152 Seconds Solving timestep 30/45: 1 Hour, 1152 Seconds -> 1 Hour, 1350 Seconds @@ -90,8 +187,101 @@ Solving timestep 37/45: 1 Hour, 2538 Seconds -> 1 Hour, 2736 Seconds Solving timestep 38/45: 1 Hour, 2736 Seconds -> 1 Hour, 2934 Seconds Solving timestep 39/45: 1 Hour, 2934 Seconds -> 1 Hour, 3132 Seconds +Solving timestep 40/45: 1 Hour, 3132 Seconds -> 1 Hour, 3330 Seconds +Solving timestep 41/45: 1 Hour, 3330 Seconds -> 1 Hour, 3528 Seconds +Solving timestep 42/45: 1 Hour, 3528 Seconds -> 2 Hours, 126 Seconds +Solving timestep 43/45: 2 Hours, 126 Seconds -> 2 Hours, 324 Seconds +Solving timestep 44/45: 2 Hours, 324 Seconds -> 2 Hours, 522 Seconds +Solving timestep 45/45: 2 Hours, 522 Seconds -> 2 Hours, 720 Seconds Switch control type from CV_charge2 to CC_discharge1 -*** Simulation complete. Solved 39 control steps in 10 Seconds, 117 Milliseconds (termination triggered by stopFunction) ***
get the states
states = output.states;
 
% extract the time and voltage quantities
time = cellfun(@(state) state.time, states);
voltage = cellfun(@(state) state.('Control').E, states);
current = cellfun(@(state) state.('Control').I, states);
 
% calculate the capacity
capacity = time .* current;
 
% plot the discharge curve in the figure
plot(time/hour, voltage, '-', 'linewidth', 3)
 
% add plot annotations
xlabel('Capacity / mA \cdot h')
ylabel('Cell Voltage / V')

Summary

In this tutorial, we explored how to modify material parameters in BattMo.
+*** Simulation complete. Solved 45 control steps in 4 Seconds, 525 Milliseconds (termination triggered by stopFunction) ***
Warning: Number of cycles has not been given in CCCV control. We use numberOfCycles = 1.
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds +Solving timestep 02/45: 6 Seconds, 187 Milliseconds -> 12 Seconds, 375 Milliseconds +Solving timestep 03/45: 12 Seconds, 375 Milliseconds -> 24 Seconds, 750 Milliseconds +Solving timestep 04/45: 24 Seconds, 750 Milliseconds -> 49 Seconds, 500 Milliseconds +Solving timestep 05/45: 49 Seconds, 500 Milliseconds -> 99 Seconds +Solving timestep 06/45: 99 Seconds -> 198 Seconds +Solving timestep 07/45: 198 Seconds -> 396 Seconds +Solving timestep 08/45: 396 Seconds -> 594 Seconds +Solving timestep 09/45: 594 Seconds -> 792 Seconds +Solving timestep 10/45: 792 Seconds -> 990 Seconds +Solving timestep 11/45: 990 Seconds -> 1188 Seconds +Solving timestep 12/45: 1188 Seconds -> 1386 Seconds +Solving timestep 13/45: 1386 Seconds -> 1584 Seconds +Solving timestep 14/45: 1584 Seconds -> 1782 Seconds +Solving timestep 15/45: 1782 Seconds -> 1980 Seconds +Solving timestep 16/45: 1980 Seconds -> 2178 Seconds +Solving timestep 17/45: 2178 Seconds -> 2376 Seconds +Solving timestep 18/45: 2376 Seconds -> 2574 Seconds +Solving timestep 19/45: 2574 Seconds -> 2772 Seconds +Solving timestep 20/45: 2772 Seconds -> 2970 Seconds +Solving timestep 21/45: 2970 Seconds -> 3168 Seconds +Solving timestep 22/45: 3168 Seconds -> 3366 Seconds +Solving timestep 23/45: 3366 Seconds -> 3564 Seconds +Solving timestep 24/45: 3564 Seconds -> 1 Hour, 162 Seconds +switch control from CC_discharge1 to CC_discharge2 +Switch control type from CC_discharge2 to CC_charge1 +Solving timestep 25/45: 1 Hour, 162 Seconds -> 1 Hour, 360 Seconds +Solving timestep 26/45: 1 Hour, 360 Seconds -> 1 Hour, 558 Seconds +Solving timestep 27/45: 1 Hour, 558 Seconds -> 1 Hour, 756 Seconds +Solving timestep 28/45: 1 Hour, 756 Seconds -> 1 Hour, 954 Seconds +Solving timestep 29/45: 1 Hour, 954 Seconds -> 1 Hour, 1152 Seconds +Solving timestep 30/45: 1 Hour, 1152 Seconds -> 1 Hour, 1350 Seconds +Solving timestep 31/45: 1 Hour, 1350 Seconds -> 1 Hour, 1548 Seconds +Solving timestep 32/45: 1 Hour, 1548 Seconds -> 1 Hour, 1746 Seconds +Solving timestep 33/45: 1 Hour, 1746 Seconds -> 1 Hour, 1944 Seconds +Solving timestep 34/45: 1 Hour, 1944 Seconds -> 1 Hour, 2142 Seconds +Solving timestep 35/45: 1 Hour, 2142 Seconds -> 1 Hour, 2340 Seconds +Solving timestep 36/45: 1 Hour, 2340 Seconds -> 1 Hour, 2538 Seconds +Solving timestep 37/45: 1 Hour, 2538 Seconds -> 1 Hour, 2736 Seconds +Solving timestep 38/45: 1 Hour, 2736 Seconds -> 1 Hour, 2934 Seconds +Solving timestep 39/45: 1 Hour, 2934 Seconds -> 1 Hour, 3132 Seconds +Switch control type from CV_charge2 to CC_discharge1 +*** Simulation complete. Solved 39 control steps in 4 Seconds, 661 Milliseconds (termination triggered by stopFunction) ***
Warning: Number of cycles has not been given in CCCV control. We use numberOfCycles = 1.
Warning: Both the total time and the number of cycles are given.\nWe do not use the given total time value but compute it instead from the number of cycles.
Solving timestep 01/45: -> 6 Seconds, 187 Milliseconds +Solving timestep 02/45: 6 Seconds, 187 Milliseconds -> 12 Seconds, 375 Milliseconds +Solving timestep 03/45: 12 Seconds, 375 Milliseconds -> 24 Seconds, 750 Milliseconds +Solving timestep 04/45: 24 Seconds, 750 Milliseconds -> 49 Seconds, 500 Milliseconds +Solving timestep 05/45: 49 Seconds, 500 Milliseconds -> 99 Seconds +Solving timestep 06/45: 99 Seconds -> 198 Seconds +Solving timestep 07/45: 198 Seconds -> 396 Seconds +Solving timestep 08/45: 396 Seconds -> 594 Seconds +Solving timestep 09/45: 594 Seconds -> 792 Seconds +Solving timestep 10/45: 792 Seconds -> 990 Seconds +Solving timestep 11/45: 990 Seconds -> 1188 Seconds +Solving timestep 12/45: 1188 Seconds -> 1386 Seconds +Solving timestep 13/45: 1386 Seconds -> 1584 Seconds +Solving timestep 14/45: 1584 Seconds -> 1782 Seconds +Solving timestep 15/45: 1782 Seconds -> 1980 Seconds +Solving timestep 16/45: 1980 Seconds -> 2178 Seconds +Solving timestep 17/45: 2178 Seconds -> 2376 Seconds +Solving timestep 18/45: 2376 Seconds -> 2574 Seconds +Solving timestep 19/45: 2574 Seconds -> 2772 Seconds +Solving timestep 20/45: 2772 Seconds -> 2970 Seconds +Solving timestep 21/45: 2970 Seconds -> 3168 Seconds +Solving timestep 22/45: 3168 Seconds -> 3366 Seconds +Solving timestep 23/45: 3366 Seconds -> 3564 Seconds +Solving timestep 24/45: 3564 Seconds -> 1 Hour, 162 Seconds +switch control from CC_discharge1 to CC_discharge2 +Switch control type from CC_discharge2 to CC_charge1 +Solving timestep 25/45: 1 Hour, 162 Seconds -> 1 Hour, 360 Seconds +Solving timestep 26/45: 1 Hour, 360 Seconds -> 1 Hour, 558 Seconds +Solving timestep 27/45: 1 Hour, 558 Seconds -> 1 Hour, 756 Seconds +Solving timestep 28/45: 1 Hour, 756 Seconds -> 1 Hour, 954 Seconds +Solving timestep 29/45: 1 Hour, 954 Seconds -> 1 Hour, 1152 Seconds +Solving timestep 30/45: 1 Hour, 1152 Seconds -> 1 Hour, 1350 Seconds +Solving timestep 31/45: 1 Hour, 1350 Seconds -> 1 Hour, 1548 Seconds +Solving timestep 32/45: 1 Hour, 1548 Seconds -> 1 Hour, 1746 Seconds +Solving timestep 33/45: 1 Hour, 1746 Seconds -> 1 Hour, 1944 Seconds +Solving timestep 34/45: 1 Hour, 1944 Seconds -> 1 Hour, 2142 Seconds +Solving timestep 35/45: 1 Hour, 2142 Seconds -> 1 Hour, 2340 Seconds +Solving timestep 36/45: 1 Hour, 2340 Seconds -> 1 Hour, 2538 Seconds +Solving timestep 37/45: 1 Hour, 2538 Seconds -> 1 Hour, 2736 Seconds +Solving timestep 38/45: 1 Hour, 2736 Seconds -> 1 Hour, 2934 Seconds +Solving timestep 39/45: 1 Hour, 2934 Seconds -> 1 Hour, 3132 Seconds +Solver did not converge in 10 iterations for timestep of length 198 Seconds. Cutting timestep. +Solver did not converge in 10 iterations for timestep of length 99 Seconds. Cutting timestep. +Solver did not converge in 10 iterations for timestep of length 49 Seconds, 500 Milliseconds. Cutting timestep. +Switch control type from CV_charge2 to CC_discharge1 +*** Simulation complete. Solved 39 control steps in 22 Seconds, 707 Milliseconds (termination triggered by stopFunction) ***
hold off
xlabel('Capacity / mA \cdot h')
ylabel('Voltage / V')
legend('t_{NE} = 16 µm', 't_{NE} = 32 µm', 't_{NE} = 64 µm','t_{NE} = 128 µm')
We see that increasing the thickness of the negative electrode beyond a certain value has no effect on the capacity of the cell, since the Li source (at the cathode and its thickness) isn't changing.

Summary

In this tutorial, we explored how to modify control parameters in BattMo and simulate CC-CV cycling of NMC-Graphite cell at a CRate of 1 and also saw the effect of thickness on cell capacity.

\ No newline at end of file diff --git a/_static/notebooks/tutorial_6_simulate_thermal_performance_live.html b/_static/notebooks/tutorial_6_simulate_thermal_performance_live.html index 0545a9a8..d4569f3b 100644 --- a/_static/notebooks/tutorial_6_simulate_thermal_performance_live.html +++ b/_static/notebooks/tutorial_6_simulate_thermal_performance_live.html @@ -32,15 +32,18 @@ .inlineElement .textElement {} .embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement { min-height: 16px;} .rightPaneElement .textElement { padding-top: 2px; padding-left: 9px;} -.S8 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S9 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } -.S11 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

Tutorial 6 - Simulate Thermal Performance

Introduction

In this tutorial, we will use a P4D model to simulate the thermal performance. We'll use the same model from Tutorial 1

Setup material properties

jsonstruct_material = parseBattmoJson('Examples/jsondatafiles/sample_input.json');
jsonstruct_material.include_current_collectors = true;

Setup the geometry

We use a 3D geometry where it is easier to visualize the thermal effects.
jsonstruct_geometry = parseBattmoJson('Examples/JsonDataFiles/geometry3d.json');
disp(jsonstruct_geometry)
include_current_collectors: 1 - Geometry: [1x1 struct] - NegativeElectrode: [1x1 struct] - PositiveElectrode: [1x1 struct] - Separator: [1x1 struct] - ThermalModel: [1x1 struct]
We merge the json structure. We get a warning for each field that gets different values from the given inputs. The rule is that the first input takes precedence, the warning can be switched off be setting the 'warn' option to false in the call to mergeJsonStructs.
jsonstruct = mergeJsonStructs({jsonstruct_geometry , ...
jsonstruct_material});
parameter Geometry.case is assigned twice with different values. we use the value from first jsonstruct. +.S8 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: 400; text-align: left; } +.S9 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px; } +.S10 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: pre-wrap; } +.S11 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S12 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S13 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S14 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

Tutorial 6 - Simulate Thermal Performance

Introduction

In this tutorial, we will use a P4D model to simulate the thermal performance. We'll use the same input json from Tutorial 1. Note: You need the MinGW-w64 compiler to build MEX files, a MATLAB® interface to a C++ library for running this notebook. Follow the instructions here to install it.

Setup material properties

jsonstruct_material = parseBattmoJson('Examples/jsondatafiles/sample_input.json');
jsonstruct_material.include_current_collectors = true;

Setup the geometry

We use a 3D geometry where it is easier to visualize the thermal effects. This 3d geometry file also contains the parameters externalHeatTransferCoefficientTab and externalHeatTransferCoefficient needed as input for the thermal model.
jsonstruct_geometry = parseBattmoJson('Examples/JsonDataFiles/geometry3d.json');
disp(jsonstruct_geometry)
include_current_collectors: 1 + Geometry: [1×1 struct] + NegativeElectrode: [1×1 struct] + PositiveElectrode: [1×1 struct] + Separator: [1×1 struct] + ThermalModel: [1×1 struct]
We merge the json structure. We get a warning for each field that gets different values from the given inputs. The rule is that the first input takes precedence, the warning can be switched off be setting the 'warn' option to false in the call to mergeJsonStructs.
jsonstruct = mergeJsonStructs({jsonstruct_geometry , ...
jsonstruct_material});
parameter Geometry.case is assigned twice with different values. we use the value from first jsonstruct. parameter NegativeElectrode.Coating.thickness is assigned twice with different values. we use the value from first jsonstruct. parameter NegativeElectrode.Coating.N is assigned twice with different values. we use the value from first jsonstruct. parameter NegativeElectrode.CurrentCollector.N is assigned twice with different values. we use the value from first jsonstruct. @@ -49,7 +52,7 @@ parameter PositiveElectrode.CurrentCollector.N is assigned twice with different values. we use the value from first jsonstruct. parameter Separator.thickness is assigned twice with different values. we use the value from first jsonstruct. parameter Separator.N is assigned twice with different values. we use the value from first jsonstruct. -parameter ThermalModel.externalHeatTransferCoefficient is assigned twice with different values. we use the value from first jsonstruct.
 
jsonstruct.use_thermal = true;
We setup the model using the json structure.
model = setupModelFromJson(jsonstruct);

Run the simulation

We run the simulation. We have used here a standard discharge control at C-Rate=1.
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds +parameter ThermalModel.externalHeatTransferCoefficient is assigned twice with different values. we use the value from first jsonstruct.
We now set the use_thermal flag in the input to true. The thermal parameters are given for each component in the input file. They typically consists of
  • Specific Heat Capacity
  • Thermal Conducitivty
From those, we compute effective quantities. For the coating, it is done by processing the corresponding properties of the constituents. The effective thermal conductivity takes into account the volume fraction. The parameters needed to compute the heat exchange with the exterior are
  • The external temperature (K) (given in the input)
  • The external heat transfer coefficient (W/m²/s) (given in the input)
The external heat transfer coefficient can take different values for each external face of the model. It will then depend on the chosen geometrical domain: A default value used for all external faces which is overwritten by the value given for the tabs. In our example we use 100 W/m²/s and 1000 W/m²/s, respectively.
jsonstruct.use_thermal = true;

Setup the model for inspection

When we run the simulation using function runBatteryJson.m, the model is setup. In the case where we want to setup the model for inspection, prior to simulation, we can use the function setupModelFromJson.m
model = setupModelFromJson(jsonstruct);

Run the simulation

We run the simulation. We have used here a standard discharge control at C-Rate=1.
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds @@ -83,60 +86,32 @@ Solving timestep 32/45: 3276 Seconds -> 3402 Seconds Solving timestep 33/45: 3402 Seconds -> 3528 Seconds Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds -*** Simulation complete. Solved 34 control steps in 44 Seconds, 813 Milliseconds (termination triggered by stopFunction) ***

Visualisation of the results

We plot the voltage versus time for the simulation.
time = output.time;
E = output.E;
 
set(0, 'defaulttextfontsize', 15);
set(0, 'defaultaxesfontsize', 15);
set(0, 'defaultlinelinewidth', 3);
 
figure
plot(time/hour, E)
title('Voltage / V');
xlabel('time / h');
ylabel('voltage / V');
We plot the minimum and maximum values of the temperature in the model as a function of time. The temperature for each grid cell, is stored in state.ThermalModel.T. For each computed step, we extract the minimum and maximum value of the temperature
T0 = PhysicalConstants.absoluteTemperature;
 
states = output.states;
Tmin = cellfun(@(state) min(state.ThermalModel.T + T0), states);
Tmax = cellfun(@(state) max(state.ThermalModel.T + T0), states);
 
figure
hold on
plot(time / hour, Tmin, 'displayname', 'min T');
plot(time / hour, Tmax, 'displayname', 'max T');
title('Temperature / C')
xlabel('time / h');
ylabel('Temperature / C');
 
legend show

Temperature distribution at final time

We have access to the temperature distribution for all the time steps. Here, we plot the temperature on the 3D model. We find a extensive set of plotting functions in MRST. You may be interested to have a look at the Visualization Tutorial. We use plotCellData to plot the temperature. The function plotToolbar.m is also convenient to visualize values inside the domain in a interactive way.
state = states{end}
state = struct with fields:
Electrolyte: [1x1 struct] - NegativeElectrode: [1x1 struct] - PositiveElectrode: [1x1 struct] - ThermalModel: [1x1 struct] - Control: [1x1 struct] +*** Simulation complete. Solved 34 control steps in 25 Seconds, 370 Milliseconds (termination triggered by stopFunction) ***

Visualisation of the results

We plot the voltage versus time for the simulation.
time = output.time;
E = output.E;
 
set(0, 'defaulttextfontsize', 15);
set(0, 'defaultaxesfontsize', 15);
set(0, 'defaultlinelinewidth', 3);
 
figure
plot(time/hour, E)
title('Voltage / V');
xlabel('time / h');
ylabel('voltage / V');
We plot the minimum and maximum values of the temperature in the model as a function of time. The temperature for each grid cell, is stored in state.ThermalModel.T. For each computed step, we extract the minimum and maximum value of the temperature. We notice that there is very little temperature variation. The reason is that we have a small cell which is very thin, with a lot of external contact where heat can be released.
T0 = PhysicalConstants.absoluteTemperature;
 
states = output.states;
Tmin = cellfun(@(state) min(state.ThermalModel.T + T0), states);
Tmax = cellfun(@(state) max(state.ThermalModel.T + T0), states);
 
figure
hold on
plot(time / hour, Tmin, 'displayname', 'min T');
plot(time / hour, Tmax, 'displayname', 'max T');
title('Temperature / C')
xlabel('time / h');
ylabel('Temperature / C');
 
legend show

Temperature distribution at final time

We have access to the temperature distribution for all the time steps. Here, we plot the temperature on the 3D model. We find a extensive set of plotting functions in MRST. You may be interested to have a look at the Visualization Tutorial. We use plotCellData to plot the temperature. The function plotToolbar.m is also convenient to visualize values inside the domain in a interactive way.
state = states{end}
state = struct with fields:
Electrolyte: [1×1 struct] + NegativeElectrode: [1×1 struct] + PositiveElectrode: [1×1 struct] + ThermalModel: [1×1 struct] + Control: [1×1 struct] time: 3.6017e+03 -
figure
plotCellData(model.ThermalModel.grid, ...
state.ThermalModel.T + T0);
colorbar
title('Temperature / C');
view([50, 50]);

The external heat transfer coefficient

The external heat transfer coefficient and the external temperature have a strong influence on the value of the temperature inside the model. The value of the external temperature has been given in the json input and passed to the model.
disp(model.ThermalModel.externalTemperature);
298.1500
The values of the heat transfer coefficient have been processed from the json structure. For this 3D geometry, the external heat transfer coefficient is given by a value that is used at the tab
disp(jsonstruct.ThermalModel.externalHeatTransferCoefficientTab);
1000
and an other that is used for the remaining external faces,
disp(jsonstruct.ThermalModel.externalHeatTransferCoefficient);
100
From those values, when the model is setup, a value is assigned to all the external faces of the model. This value is stored in model.ThermalModel.externalHeatTransferCoefficient. Let us plot its value the 3D model. In the coupling term of the thermal model (coupling to the exterior model), we get the list of the face indices that are coupled thermally, which are in fact allthe external faces.
extfaceind = model.ThermalModel.couplingTerm.couplingfaces;
G = model.ThermalModel.grid;
nf = G.faces.num;
We create a vector with one value per face and we set this value equal to the external heat transfer coefficient for the corresponding external face.
val = NaN(nf, 1);
val(extfaceind) = model.ThermalModel.externalHeatTransferCoefficient;
 
figure
plotFaceData(G, val, 'edgecolor', 'black');
axis equal
view([50, 20]);
title('External Heat Transfer Coefficient / W/s/m^2')
colorbar
Let us modify the heat transfer coefficients and rerun the simulation to observe the effects.
jsonstruct.ThermalModel.externalHeatTransferCoefficientTab = 100;
jsonstruct.ThermalModel.externalHeatTransferCoefficient = 0;
We run the model for the new values
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds -Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds -Solving timestep 03/45: 7 Seconds, 875 Milliseconds -> 15 Seconds, 750 Milliseconds -Solving timestep 04/45: 15 Seconds, 750 Milliseconds -> 31 Seconds, 500 Milliseconds -Solving timestep 05/45: 31 Seconds, 500 Milliseconds -> 63 Seconds -Solving timestep 06/45: 63 Seconds -> 126 Seconds -Solving timestep 07/45: 126 Seconds -> 252 Seconds -Solving timestep 08/45: 252 Seconds -> 378 Seconds -Solving timestep 09/45: 378 Seconds -> 504 Seconds -Solving timestep 10/45: 504 Seconds -> 630 Seconds -Solving timestep 11/45: 630 Seconds -> 756 Seconds -Solving timestep 12/45: 756 Seconds -> 882 Seconds -Solving timestep 13/45: 882 Seconds -> 1008 Seconds -Solving timestep 14/45: 1008 Seconds -> 1134 Seconds -Solving timestep 15/45: 1134 Seconds -> 1260 Seconds -Solving timestep 16/45: 1260 Seconds -> 1386 Seconds -Solving timestep 17/45: 1386 Seconds -> 1512 Seconds -Solving timestep 18/45: 1512 Seconds -> 1638 Seconds -Solving timestep 19/45: 1638 Seconds -> 1764 Seconds -Solving timestep 20/45: 1764 Seconds -> 1890 Seconds -Solving timestep 21/45: 1890 Seconds -> 2016 Seconds -Solving timestep 22/45: 2016 Seconds -> 2142 Seconds -Solving timestep 23/45: 2142 Seconds -> 2268 Seconds -Solving timestep 24/45: 2268 Seconds -> 2394 Seconds -Solving timestep 25/45: 2394 Seconds -> 2520 Seconds -Solving timestep 26/45: 2520 Seconds -> 2646 Seconds -Solving timestep 27/45: 2646 Seconds -> 2772 Seconds -Solving timestep 28/45: 2772 Seconds -> 2898 Seconds -Solving timestep 29/45: 2898 Seconds -> 3024 Seconds -Solving timestep 30/45: 3024 Seconds -> 3150 Seconds -Solving timestep 31/45: 3150 Seconds -> 3276 Seconds -Solving timestep 32/45: 3276 Seconds -> 3402 Seconds -Solving timestep 33/45: 3402 Seconds -> 3528 Seconds -Solving timestep 34/45: 3528 Seconds -> 1 Hour, 54 Seconds -*** Simulation complete. Solved 34 control steps in 41 Seconds, 805 Milliseconds (termination triggered by stopFunction) ***
We plot the results
states = output.states;
 
Tmin = cellfun(@(state) min(state.ThermalModel.T + T0), states);
Tmax = cellfun(@(state) max(state.ThermalModel.T + T0), states);
time = output.time;
 
figure
hold on
plot(time / hour, Tmin, 'displayname', 'min T');
plot(time / hour, Tmax, 'displayname', 'max T');
title('Temperature / C')
xlabel('time / h');
ylabel('Temperature / C');
 
legend show
+
figure
plotCellData(model.ThermalModel.grid, ...
state.ThermalModel.T + T0);
colorbar
title('Temperature / C');
view([50, 50]);

The external heat transfer coefficient

The external heat transfer coefficient and the external temperature have a strong influence on the value of the temperature inside the model.
disp(model.ThermalModel.externalTemperature);
298.1500
The values of the heat transfer coefficient have been processed from the json structure. For this 3D geometry, the external heat transfer coefficient is given by a value that is used at the tab
disp(jsonstruct.ThermalModel.externalHeatTransferCoefficientTab);
100
and an other that is used for the remaining external faces,
disp(jsonstruct.ThermalModel.externalHeatTransferCoefficient);
0
From those values, when the model is setup, a value is assigned to all the external faces of the model. This value is stored in model.ThermalModel.externalHeatTransferCoefficient. Let us plot its value the 3D model. In the coupling term of the thermal model (coupling to the exterior model), we get the list of the face indices that are coupled thermally, which are in fact allthe external faces.
extfaceind = model.ThermalModel.couplingTerm.couplingfaces;
G = model.ThermalModel.grid;
nf = G.faces.num;
We create a vector with one value per face and we set this value equal to the external heat transfer coefficient for the corresponding external face.
val = NaN(nf, 1);
val(extfaceind) = model.ThermalModel.externalHeatTransferCoefficient;
 
figure
plotFaceData(G, val, 'edgecolor', 'black');
axis equal
view([50, 20]);
title('External Heat Transfer Coefficient / W/s/m^2')
colorbar
Let us modify the heat transfer coefficients and rerun the simulation to observe the effects. We set the default heat exchange to zero so that heat exchange can only occur through the tabs.
jsonstruct.ThermalModel.externalHeatTransferCoefficientTab = 100;
jsonstruct.ThermalModel.externalHeatTransferCoefficient = 0;
We run the model for the new values
output = runBatteryJson(jsonstruct);
Solving timestep 01/45: -> 3 Seconds, 937 Milliseconds +Solving timestep 02/45: 3 Seconds, 937 Milliseconds -> 7 Seconds, 875 Milliseconds
We plot the results, we see we obtain higher temperature during the operation when the heat release happens only via the tabs.
states = output.states;
 
Tmin = cellfun(@(state) min(state.ThermalModel.T + T0), states);
Tmax = cellfun(@(state) max(state.ThermalModel.T + T0), states);
time = output.time;
 
figure
hold on
plot(time / hour, Tmin, 'displayname', 'min T');
plot(time / hour, Tmax, 'displayname', 'max T');
title('Temperature / C')
xlabel('time / h');
ylabel('Temperature / C');
 
legend show