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