GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/chemistry/electrochemistry/electrochemistry.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 66 67 98.5%
Functions: 16 16 100.0%
Branches: 73 124 58.9%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Chemistry
10 * \brief Electrochemical model for a fuel cell application.
11 */
12 #ifndef DUMUX_ELECTROCHEMISTRY_HH
13 #define DUMUX_ELECTROCHEMISTRY_HH
14
15 #include <cmath>
16
17 #include <dumux/common/parameters.hh>
18 #include <dumux/common/exceptions.hh>
19 #include <dumux/discretization/method.hh>
20 #include <dumux/material/constants.hh>
21
22 #include <dumux/material/fluidsystems/h2on2o2.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup Chemistry
28 * \brief The type of electrochemistry models implemented here (Ochs (2008) \cite ochs2008 or Acosta et al. (2006) \cite A3:acosta:2006 )
29 */
30 enum ElectroChemistryModel { Ochs, Acosta };
31
32 /*!
33 * \ingroup Chemistry
34 * \brief This class calculates source terms and current densities for fuel cells
35 * with the electrochemical models suggested by Ochs (2008) \cite ochs2008 or Acosta et al. (2006) \cite A3:acosta:2006
36 * \todo TODO: Scalar type should be extracted from VolumeVariables!
37 * \todo TODO: This shouldn't depend on grid and discretization!!
38 */
39 template <class Scalar, class Indices, class FluidSystem, class GridGeometry, ElectroChemistryModel electroChemistryModel>
40 class ElectroChemistry
41 {
42
43 using Constant = Dumux::Constants<Scalar>;
44
45 enum
46 {
47 //indices of the primary variables
48 pressureIdx = Indices::pressureIdx, // gas-phase pressure
49 switchIdx = Indices::switchIdx, // liquid saturation or mole fraction
50 };
51 enum
52 {
53 //equation indices
54 contiH2OEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx,
55 contiO2EqIdx = Indices::conti0EqIdx + FluidSystem::O2Idx,
56 };
57
58 enum
59 {
60 // phase indices
61 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
62 gasPhaseIdx = FluidSystem::gasPhaseIdx
63 };
64
65 using GridView = typename GridGeometry::GridView;
66 static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethods::box;
67 using GlobalPosition = typename Dune::FieldVector<typename GridView::ctype, GridView::dimensionworld>;
68
69 public:
70 /*!
71 * \brief Calculates reaction sources with an electrochemical model approach.
72 *
73 * \param values The primary variable vector
74 * \param currentDensity The current density
75 * \param paramGroup The group containing the needed parameters
76 *
77 * For this method, the \a values parameter stores source values
78 */
79 template<class SourceValues>
80 40950 static void reactionSource(SourceValues &values, Scalar currentDensity,
81 const std::string& paramGroup = "")
82 {
83 // correction to account for actually relevant reaction area
84 // current density has to be divided by the half length of the box
85 // \todo Do we have multiply with the electrochemically active surface area (ECSA) here instead?
86
4/6
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 20471 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 40950 static Scalar gridYMax = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.UpperRight")[1];
87
4/6
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 20471 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 40950 static Scalar nCellsY = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.Cells")[1];
88
89 // Warning: This assumes the reaction layer is always just one cell (cell-centered) or half a box (box) thick
90 40950 const auto lengthBox = gridYMax/nCellsY;
91 if (isBox)
92 35112 currentDensity *= 2.0/lengthBox;
93 else
94 5838 currentDensity *= 1.0/lengthBox;
95
96
4/6
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 20471 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 40950 static Scalar transportNumberH2O = getParam<Scalar>("ElectroChemistry.TransportNumberH20");
97
98 //calculation of flux terms with Faraday equation
99 40950 values[contiH2OEqIdx] = currentDensity/(2*Constant::F); //reaction term in reaction layer
100 40950 values[contiH2OEqIdx] += currentDensity/Constant::F*transportNumberH2O; //osmotic term in membrane
101 40950 values[contiO2EqIdx] = -currentDensity/(4*Constant::F); //O2-equation
102 40950 }
103
104 /*!
105 * \brief Newton solver for calculation of the current density.
106 * \param volVars The volume variables
107 * \returns The current density in A/m^2
108 */
109 template<class VolumeVariables>
110 44667 static Scalar calculateCurrentDensity(const VolumeVariables &volVars)
111 {
112
113
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 44664 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 44667 static int maxIter = getParam<int>("ElectroChemistry.MaxIterations");
114
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 44664 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 44667 static Scalar specificResistance = getParam<Scalar>("ElectroChemistry.SpecificResistance");
115
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 44664 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 44667 static Scalar reversibleVoltage = getParam<Scalar>("ElectroChemistry.ReversibleVoltage");
116
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 44664 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 44667 static Scalar cellVoltage = getParam<Scalar>("ElectroChemistry.CellVoltage");
117
118 //initial guess for the current density and initial newton solver parameters
119 44667 Scalar currentDensity = reversibleVoltage - cellVoltage - 0.5;
120 44667 Scalar increment = 1e-4;
121 44667 Scalar deltaCurrentDensity = currentDensity*increment;
122 44667 Scalar deltaVoltage = 1.0;
123 44667 int iterations = 0;
124
125 //Newton Solver for current Density
126 using std::abs;
127
4/4
✓ Branch 0 taken 134001 times.
✓ Branch 1 taken 44667 times.
✓ Branch 2 taken 134001 times.
✓ Branch 3 taken 44667 times.
 357336 while (abs(deltaVoltage) > 1e-6)
128 {
129
130 134001 Scalar activationLosses = calculateActivationLosses_(volVars, currentDensity);
131 134001 Scalar activationLossesNext = calculateActivationLosses_(volVars, currentDensity+deltaCurrentDensity);
132 134001 Scalar concentrationLosses = calculateConcentrationLosses_(volVars);
133 134001 Scalar activationLossesDiff = activationLossesNext - activationLosses;
134 134001 Scalar sw = volVars.saturation(liquidPhaseIdx);
135
136 if(electroChemistryModel == ElectroChemistryModel::Acosta)
137 {
138 // Acosta calculation
139 deltaVoltage = currentDensity*specificResistance - reversibleVoltage + cellVoltage + activationLosses + concentrationLosses;
140
141 currentDensity = currentDensity - (deltaVoltage*deltaCurrentDensity)/(deltaCurrentDensity*specificResistance + activationLossesDiff);
142
143 activationLosses = calculateActivationLosses_(volVars, currentDensity);
144
145 deltaVoltage = currentDensity*specificResistance - reversibleVoltage + cellVoltage + activationLosses + concentrationLosses;
146
147 iterations++;
148
149 }
150 else
151 {
152 // Ochs & other calculation
153 134001 deltaVoltage = currentDensity*specificResistance/(1-sw) - reversibleVoltage + cellVoltage + activationLosses + concentrationLosses;
154
155 134001 currentDensity = currentDensity - (deltaVoltage*deltaCurrentDensity)/(deltaCurrentDensity*specificResistance/(1-sw) + activationLossesDiff);
156
157 134001 activationLosses = calculateActivationLosses_(volVars, currentDensity);
158
159 134001 deltaVoltage = currentDensity*specificResistance/(1-sw) - reversibleVoltage + cellVoltage + activationLosses + concentrationLosses;
160
161 134001 iterations++;
162 }
163
164
1/2
✓ Branch 0 taken 134001 times.
✗ Branch 1 not taken.
 134001 if(iterations >= maxIter)
165 {
166 ✗ DUNE_THROW(NumericalProblem, "Newton solver for electrochemistry didn't converge");
167 }
168 }
169
170 //conversion from [A/cm^2] to [A/m^2]
171 44667 return currentDensity*10000;
172 }
173
174 private:
175
176 /*!
177 * \brief Calculation of the activation losses
178 * \param volVars The volume variables
179 * \param currentDensity The current density
180 */
181 template<class VolumeVariables>
182 402003 static Scalar calculateActivationLosses_(const VolumeVariables &volVars, const Scalar currentDensity)
183 {
184
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 402000 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 402003 static Scalar refO2PartialPressure = getParam<Scalar>("ElectroChemistry.RefO2PartialPressure");
185
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 402000 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 402003 static Scalar numElectrons = getParam<Scalar>("ElectroChemistry.NumElectrons");
186
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 402000 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 402003 static Scalar transferCoefficient = getParam<Scalar>("ElectroChemistry.TransferCoefficient");
187
188 //Saturation sw for Acosta calculation
189 402003 Scalar sw = volVars.saturation(liquidPhaseIdx);
190 //Calculate prefactor
191 804006 Scalar preFactor = Constant::R*volVars.fluidState().temperature()/transferCoefficient/Constant::F/numElectrons;
192 //Get partial pressure of O2 in the gas phase
193 1206009 Scalar pO2 = volVars.pressure(gasPhaseIdx) * volVars.fluidState().moleFraction(gasPhaseIdx, FluidSystem::O2Idx);
194
195 402003 Scalar losses = 0.0;
196 //Calculate activation losses
197 using std::log;
198 using std::abs;
199 if(electroChemistryModel == ElectroChemistryModel::Acosta)
200 {
201 losses = preFactor
202 *( log(abs(currentDensity)/abs(exchangeCurrentDensity_(volVars)))
203 - log(pO2/refO2PartialPressure)
204 - log(1 - sw)
205 );
206 }
207 else
208 {
209 402003 losses = preFactor
210 804006 *( log(abs(currentDensity)/abs(exchangeCurrentDensity_(volVars)))
211 402003 - log(pO2/refO2PartialPressure)
212 );
213 }
214 402003 return losses;
215 }
216
217
218 /*!
219 * \brief Calculation of concentration losses.
220 * \param volVars The volume variables
221 */
222 template<class VolumeVariables>
223 134001 static Scalar calculateConcentrationLosses_(const VolumeVariables &volVars)
224 {
225
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 133998 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 134001 static Scalar pO2Inlet = getParam<Scalar>("ElectroChemistry.pO2Inlet");
226
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 133998 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 134001 static Scalar numElectrons = getParam<Scalar>("ElectroChemistry.NumElectrons");
227
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 133998 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 134001 static Scalar transferCoefficient = getParam<Scalar>("ElectroChemistry.TransferCoefficient");
228
229 //Calculate preFactor
230 134001 Scalar preFactor = Constant::R*volVars.temperature()/transferCoefficient/Constant::F/numElectrons;
231 //Get partial pressure of O2 in the gas phase
232 402003 Scalar pO2 = volVars.pressure(gasPhaseIdx) * volVars.fluidState().moleFraction(gasPhaseIdx, FluidSystem::O2Idx);
233
234 134001 Scalar losses = 0.0;
235 //Calculate concentration losses
236 using std::log;
237 if(electroChemistryModel == ElectroChemistryModel::Acosta)
238 {
239 losses = -1.0*preFactor*(transferCoefficient/2)*log(pO2/pO2Inlet);
240 }else
241 {
242 // +1 is the Nernst part of the equation
243 134001 losses = -1.0*preFactor*(transferCoefficient/2+1)*log(pO2/pO2Inlet);
244 }
245
246 134001 return losses;
247 }
248
249
250 /*!
251 * \brief Calculation of the exchange current density.
252 * \param volVars The volume variables
253 */
254 template<class VolumeVariables>
255 402003 static Scalar exchangeCurrentDensity_(const VolumeVariables &volVars)
256 {
257 using std::exp;
258
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 402000 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 402003 static Scalar activationBarrier = getParam<Scalar>("ElectroChemistry.ActivationBarrier");
259
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 402000 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 402003 static Scalar surfaceIncreasingFactor = getParam<Scalar>("ElectroChemistry.SurfaceIncreasingFactor");
260
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 402000 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 402003 static Scalar refTemperature = getParam<Scalar>("ElectroChemistry.RefTemperature");
261
4/6
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 402000 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 402003 static Scalar refCurrentDensity = getParam<Scalar>("ElectroChemistry.RefCurrentDensity");
262
263 804006 Scalar T = volVars.fluidState().temperature();
264 402003 Scalar refExchangeCurrentDensity = -1.0
265 402003 * refCurrentDensity
266 402003 * surfaceIncreasingFactor
267 402003 * exp(-1.0 * activationBarrier / Constant::R * (1/T-1/refTemperature));
268
269 402003 return refExchangeCurrentDensity;
270 }
271 };
272
273 } // end namespace Dumux
274
275 #endif
276