Skip to content

Commit

Permalink
Use primitive plotting
Browse files Browse the repository at this point in the history
  • Loading branch information
@augustjohansson
augustjohansson committed May 12, 2022
1 parent d546bda commit 48d9c17
Showing 1 changed file with 82 additions and 40 deletions.
122 changes: 82 additions & 40 deletions Examples/runBattery2D.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,71 +35,96 @@
%% Setup the geometry and computational mesh
% Here, we setup the 2D computational mesh that will be used for the
% simulation. The required discretization parameters are already included
% in the class BatteryGenerator2D.
% in the class BatteryGenerator2D.
gen = BatteryGenerator2D();

% Now, we update the paramobj with the properties of the mesh.
paramobj = gen.updateBatteryInputParams(paramobj);

% !!! REMOVE THIS. SET THE RIGHT VALUES IN THE JSON !!! In this case, we
% change some of the values of the paramaters that were given in the json
% file to other values. This is done directly on the object paramobj.
% file to other values. This is done directly on the object paramobj.
paramobj.(ne).(cc).EffectiveElectricalConductivity = 1e5;
paramobj.(pe).(cc).EffectiveElectricalConductivity = 1e5;

%% Initialize the battery model.
%% Initialize the battery model.
% The battery model is initialized by sending paramobj to the Battery class
% constructor. see :class:`Battery <Battery.Battery>`.
model = Battery(paramobj);

%% Plot the mesh
% The mesh is plotted using the plotGrid() function from MRST.
% The mesh is plotted using the plotGrid() function from MRST.
colors = crameri('vik', 5);
figure
plotGrid(model.(ne).(cc).G, 'facecolor', colors(1,:), 'edgealpha', 0.5, 'edgecolor', [1, 1, 1]);
plotGrid(model.(ne).(am).G, 'facecolor', colors(2,:), 'edgealpha', 0.5, 'edgecolor', [1, 1, 1]);
plotGrid(model.(elyte).(sep).G, 'facecolor', colors(3,:), 'edgealpha', 0.5, 'edgecolor', [1, 1, 1]);
plotGrid(model.(pe).(am).G, 'facecolor', colors(4,:), 'edgealpha', 0.5, 'edgecolor', [1, 1, 1]);
plotGrid(model.(pe).(cc).G, 'facecolor', colors(5,:), 'edgealpha', 0.5, 'edgecolor', [1, 1, 1]);
fig = figure; hold on
if mrstPlatform('octave')
edgealpha = 1;
else
edgealpha = 0.5;
end
plotGrid(model.(ne).(cc).G, 'facecolor', colors(1,:), 'edgealpha', edgealpha, 'edgecolor', [1, 1, 1]);
plotGrid(model.(ne).(am).G, 'facecolor', colors(2,:), 'edgealpha', edgealpha, 'edgecolor', [1, 1, 1]);
plotGrid(model.(elyte).(sep).G, 'facecolor', colors(3,:), 'edgealpha', edgealpha, 'edgecolor', [1, 1, 1]);
plotGrid(model.(pe).(am).G, 'facecolor', colors(4,:), 'edgealpha', edgealpha, 'edgecolor', [1, 1, 1]);
plotGrid(model.(pe).(cc).G, 'facecolor', colors(5,:), 'edgealpha', edgealpha, 'edgecolor', [1, 1, 1]);
axis tight;
legend({'negative electrode current collector' , ...
'negative electrode active material' , ...
'separator' , ...
'positive electrode active material' , ...
'positive electrode current collector'}, ...
'location', 'southwest'),
setFigureStyle('quantity', 'single');
'location', 'southwest')

if mrstPlatform('octave')
style.lineWidth = 5;
lines = findobj(gcf,'Type','Line');
for i = 1:numel(lines)
set(lines(i), 'LineWidth', style.lineWidth);
end

ax = gca;
style.fontName = 'Helvetica';
style.fontSize = 24;
style.fontColor = 'k';

set(ax, ...
'FontSize', style.fontSize, ...
'FontName', style.fontName, ...
'GridColor', style.fontColor);

else
setFigureStyle('quantity', 'single');
end
drawnow();
pause(0.1);

%% 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.
% operational current.
C = computeCellCapacity(model);
CRate = 1;
inputI = (C/hour)*CRate; % current
inputI = (C/hour)*CRate; % current

%% Setup the time step schedule
%% 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.
n = 25;
dt = [];
dt = [dt; repmat(0.5e-4, n, 1).*1.5.^[1:n]'];
% operation.
n = 25;
dt = [];
dt = [dt; repmat(0.5e-4, n, 1).*1.5.^[1:n]'];
totalTime = 1.4*hour/CRate;
n = 40;
dt = [dt; repmat(totalTime/n, n, 1)];
times = [0; cumsum(dt)];
tt = times(2 : end);
step = struct('val', diff(times), 'control', ones(numel(tt), 1));
n = 40;
dt = [dt; repmat(totalTime/n, n, 1)];
times = [0; cumsum(dt)];
tt = times(2 : end);
step = struct('val', diff(times), 'control', ones(numel(tt), 1));

%% Setup the operating limits for the cell
% The maximum and minimum voltage limits for the cell are defined using
% stopping and source functions. A stopping function is used to set the
% lower voltage cutoff limit. A source function is used to set the upper
% voltage cutoff limit.
stopFunc = @(model, state, state_prev) (state.(ctrl).E < 2.0);
% voltage cutoff limit.
stopFunc = @(model, state, state_prev) (state.(ctrl).E < 2.0);
tup = 0.1; % rampup value for the current function, see rampupSwitchControl
srcfunc = @(time, I, E) rampupSwitchControl(time, tup, I, E, ...
model.Control.Imax, ...
Expand All @@ -108,43 +133,60 @@
control = struct('src', srcfunc, 'IEswitch', true);

% This control is used to set up the schedule
schedule = struct('control', control, 'step', step);
schedule = struct('control', control, 'step', step);

%% Setup the initial state of the model
% The initial state of the model is dispatched using the
% model.setupInitialState()method.
initstate = model.setupInitialState();
% model.setupInitialState()method.
initstate = model.setupInitialState();

%% Setup the properties of the nonlinear solver
nls = NonLinearSolver();
%% Setup the properties of the nonlinear solver
nls = NonLinearSolver();
% Change default maximum iteration number in nonlinear solver
nls.maxIterations = 10;
nls.maxIterations = 10;
% Change default behavior of nonlinear solver, in case of error
nls.errorOnFailure = false;
nls.errorOnFailure = false;
% Timestep selector
nls.timeStepSelector = StateChangeTimeStepSelector('TargetProps', {{ctrl, 'E'}}, ...
'targetChangeAbs', 0.03);

% Change default tolerance for nonlinear solver
model.nonlinearTolerance = 1e-5;
model.nonlinearTolerance = 1e-5;
% Set verbosity of the solver (if true, value of the residuals for every equation is given)
model.verbose = true;

%% Run simulation
[wellSols, states, report] = simulateScheduleAD(initstate, model, schedule, ...
'OutputMinisteps', true, ...
'NonLinearSolver', nls);
'NonLinearSolver', nls);

%% Process output and recover the output voltage and current from the output states.
ind = cellfun(@(x) not(isempty(x)), states);
ind = cellfun(@(x) not(isempty(x)), states);
states = states(ind);
Enew = cellfun(@(x) x.(ctrl).E, states);
Enew = cellfun(@(x) x.(ctrl).E, states);
Inew = cellfun(@(x) x.(ctrl).I, states);
time = cellfun(@(x) x.time, states);
time = cellfun(@(x) x.time, states);

%% Plot an animated summary of the results

plotDashboard(model, states, 'step', 0);
if mrstPlatform('octave')
figure
plot(time, Enew, 'linewidth', style.lineWidth, ...
time, Inew, 'linewidth', style.lineWidth);
xlabel 'time'
legend('output voltage [V]', 'output current [A]')

ax = gca;
style.fontName = 'Helvetica';
style.fontSize = 24;

set(ax, ...
'FontSize', style.fontSize, ...
'FontName', style.fontName);

grid on
else
plotDashboard(model, states, 'step', 0);
end

%{
Copyright 2021-2022 SINTEF Industry, Sustainable Energy Technology
Expand Down

0 comments on commit 48d9c17

Please sign in to comment.