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.
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.
 
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
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 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]
                 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