electrochemistry.hh 13.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program 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 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
21
 * \ingroup Chemistry
22
23
 * \brief Electrochemical model for a fuel cell application.
 */
24
25
#ifndef DUMUX_ELECTROCHEMISTRY_HH
#define DUMUX_ELECTROCHEMISTRY_HH
26

Timo Koch's avatar
Timo Koch committed
27
28
#include <cmath>

29
#include <dumux/common/properties.hh>
Timo Koch's avatar
Timo Koch committed
30
#include <dumux/common/parameters.hh>
31
#include <dumux/common/exceptions.hh>
Timo Koch's avatar
Timo Koch committed
32
#include <dumux/discretization/methods.hh>
33
34
#include <dumux/material/constants.hh>
#include <dumux/material/components/component.hh>
35
#include <dumux/material/fluidsystems/h2on2o2.hh>
36
37
38
39

namespace Dumux
{

Timo Koch's avatar
Timo Koch committed
40
/*!
41
42
 * \ingroup Chemistry
 * \brief The type of electrochemistry models implemented here (Ochs (2008) \cite ochs2008 or Acosta et al. (2006) \cite A3:acosta:2006 )
Timo Koch's avatar
Timo Koch committed
43
44
45
46
 */
enum ElectroChemistryModel { Ochs, Acosta };

/*!
47
 * \brief This class calculates source terms and current densities for fuel cells
48
 * with the electrochemical models suggested by Ochs (2008) \cite ochs2008 or Acosta et al. (2006) \cite A3:acosta:2006
Timo Koch's avatar
Timo Koch committed
49
 */
50
template <class TypeTag, ElectroChemistryModel electroChemistryModel>
Timo Koch's avatar
Timo Koch committed
51
52
class ElectroChemistry
{
53
54
55
56
57
58
59
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Element = typename GridView::template Codim<0>::Entity;
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
Timo Koch's avatar
Timo Koch committed
60

61
    using Constant = Constants<Scalar>;
Timo Koch's avatar
Timo Koch committed
62

63
    using ThisType = ElectroChemistry<TypeTag, electroChemistryModel>;
Timo Koch's avatar
Timo Koch committed
64

65
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
Timo Koch's avatar
Timo Koch committed
66
67
68
69

    enum {
        //indices of the phases
        wPhaseIdx = Indices::wPhaseIdx,
70
71
72
        nPhaseIdx = Indices::nPhaseIdx
    };
    enum {
Timo Koch's avatar
Timo Koch committed
73
74
75
        //indices of the components
        wCompIdx = FluidSystem::wCompIdx, //major component of the liquid phase
        nCompIdx = FluidSystem::nCompIdx, //major component of the gas phase
76
77
78
        O2Idx = wCompIdx + 2
    };
    enum {
Timo Koch's avatar
Timo Koch committed
79
80
81
        //indices of the primary variables
        pressureIdx = Indices::pressureIdx, //gas-phase pressure
        switchIdx = Indices::switchIdx, //liquid saturation or mole fraction
82
83
84
        temperatureIdx = FluidSystem::numComponents //temperature
    };
    enum {
Timo Koch's avatar
Timo Koch committed
85
86
87
88
        //equation indices
        conti0EqIdx = Indices::conti0EqIdx,
        contiH2OEqIdx = conti0EqIdx + wCompIdx,
        contiO2EqIdx = conti0EqIdx + wCompIdx + 2,
89
        energyEqIdx = FluidSystem::numComponents //energy equation
Timo Koch's avatar
Timo Koch committed
90
91
    };

92
    static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
93
94
95
96
97
    enum { dofCodim = isBox ? GridView::dimension : 0 };

    using GlobalPosition = typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
    using CellVector = typename Dune::FieldVector<Scalar, GridView::dimension>;

Timo Koch's avatar
Timo Koch committed
98
public:
99
    /*!
100
101
102
103
104
105
106
     * \brief Calculates reaction sources with an electrochemical model approach.
     *
     * \param values The primary variable vector
     * \param currentDensity The current density
     *
     * For this method, the \a values parameter stores source values
     */
107
    static void reactionSource(SourceValues &values,
108
                               Scalar currentDensity)
109
    {
Timo Koch's avatar
Timo Koch committed
110
111
        //correction to account for actually relevant reaction area
        //current density has to be devided by the half length of the box
112
        //\todo Do we have multiply with the electrochemically active surface area (ECSA) here instead?
113
114
        static Scalar gridYMax =getParamFromGroup<GlobalPosition>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Grid.UpperRight")[1];
        static Scalar nCellsY = getParamFromGroup<GlobalPosition>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Grid.Cells")[1];
115
116
117
118
119

        // Warning: This assumes the reaction layer is always just one cell (cell-centered) or half a box (box) thick
        const auto lengthBox = gridYMax/nCellsY;
        if (isBox)
            currentDensity *= 2.0/lengthBox;
Timo Koch's avatar
Timo Koch committed
120
        else
121
            currentDensity *= 1.0/lengthBox;
Timo Koch's avatar
Timo Koch committed
122

123
        static Scalar transportNumberH2O = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.TransportNumberH20");
Timo Koch's avatar
Timo Koch committed
124
125
126
127
128
129
130

        //calculation of flux terms with faraday equation
        values[contiH2OEqIdx] = currentDensity/(2*Constant::F);                  //reaction term in reaction layer
        values[contiH2OEqIdx] += currentDensity/Constant::F*transportNumberH2O;  //osmotic term in membrane
        values[contiO2EqIdx]  = -currentDensity/(4*Constant::F);                 //O2-equation
    }

131
    /*!
132
133
134
135
     * \brief Newton solver for calculation of the current density.
     * \param volVars The volume variables
     * \returns The current density in A/m^2
     */
136
    static Scalar calculateCurrentDensity(const VolumeVariables &volVars)
137
    {
138
        static Scalar maxIter = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.MaxIterations");
139

140
141
142
        static Scalar specificResistance = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.SpecificResistance");
        static Scalar reversibleVoltage = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.ReversibleVoltage");
        static Scalar cellVoltage = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.CellVoltage");
Timo Koch's avatar
Timo Koch committed
143

144
        //initial guess for the current density and initial newton solver parameters
Timo Koch's avatar
Timo Koch committed
145
        Scalar currentDensity = reversibleVoltage - cellVoltage - 0.5;
146
147
148
149
        Scalar increment = 1e-4;
        Scalar deltaCurrentDensity = currentDensity*increment;
        Scalar deltaVoltage = 1.0;
        int iterations = 0;
Timo Koch's avatar
Timo Koch committed
150

151
        //Newton Solver for current Density
152
153
        using std::abs;
        while (abs(deltaVoltage) > 1e-6)
154
        {
Timo Koch's avatar
Timo Koch committed
155
156
157
158
159
160
161
162

            Scalar activationLosses        = calculateActivationLosses_(volVars, currentDensity);
            Scalar activationLossesNext    = calculateActivationLosses_(volVars, currentDensity+deltaCurrentDensity);
            Scalar concentrationLosses     = calculateConcentrationLosses_(volVars);
            Scalar activationLossesDiff    = activationLossesNext - activationLosses;
            Scalar sw                      = volVars.saturation(wPhaseIdx);

            if(electroChemistryModel == ElectroChemistryModel::Acosta)
163
            {
Timo Koch's avatar
Timo Koch committed
164
165
166
167
168
                // Acosta calculation
                deltaVoltage = currentDensity*specificResistance - reversibleVoltage + cellVoltage + activationLosses + concentrationLosses;

                currentDensity = currentDensity - (deltaVoltage*deltaCurrentDensity)/(deltaCurrentDensity*specificResistance + activationLossesDiff);

169
                activationLosses = calculateActivationLosses_(volVars, currentDensity);
Timo Koch's avatar
Timo Koch committed
170
171
172

                deltaVoltage = currentDensity*specificResistance - reversibleVoltage + cellVoltage + activationLosses + concentrationLosses;

173
                iterations++;
Timo Koch's avatar
Timo Koch committed
174

175
176
177
            }
            else
            {
Timo Koch's avatar
Timo Koch committed
178
179
180
181
                // Ochs & other calculation
                deltaVoltage = currentDensity*specificResistance/(1-sw) - reversibleVoltage + cellVoltage + activationLosses + concentrationLosses;

                currentDensity = currentDensity - (deltaVoltage*deltaCurrentDensity)/(deltaCurrentDensity*specificResistance/(1-sw) + activationLossesDiff);
182
183
184

                activationLosses = calculateActivationLosses_(volVars, currentDensity);

Timo Koch's avatar
Timo Koch committed
185
                deltaVoltage = currentDensity*specificResistance/(1-sw) - reversibleVoltage + cellVoltage + activationLosses + concentrationLosses;
186
187
188

                iterations++;
            }
Timo Koch's avatar
Timo Koch committed
189
190

            if(iterations >= maxIter)
191
            {
192
                DUNE_THROW(NumericalProblem, "Newton solver for electrochemistry didn't converge");
193
194
            }
        }
Timo Koch's avatar
Timo Koch committed
195

196
197
        //conversion from [A/cm^2] to [A/m^2]
        return currentDensity*10000;
198
    }
Timo Koch's avatar
Timo Koch committed
199
200
201

private:

202
    /*!
203
204
205
206
     * \brief Calculation of the activation losses
     * \param volVars The volume variables
     * \param currentDensity The current density
     */
Timo Koch's avatar
Timo Koch committed
207
    static Scalar calculateActivationLosses_(const VolumeVariables &volVars, const Scalar currentDensity)
208
    {
209
210
211
        static Scalar refO2PartialPressure = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.RefO2PartialPressure");
        static Scalar numElectrons = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.NumElectrons");
        static Scalar transferCoefficient = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.TransferCoefficient");
Timo Koch's avatar
Timo Koch committed
212

213
        //Saturation sw for Acosta calculation
Timo Koch's avatar
Timo Koch committed
214
        Scalar sw = volVars.saturation(wPhaseIdx);
215
        //Calculate prefactor
Timo Koch's avatar
Timo Koch committed
216
        Scalar preFactor = Constant::R*volVars.fluidState().temperature()/transferCoefficient/Constant::F/numElectrons;
217
        //Get partial pressure of O2 in the gas phase
Timo Koch's avatar
Timo Koch committed
218
219
        Scalar pO2 = volVars.pressure(nPhaseIdx) * volVars.fluidState().moleFraction(nPhaseIdx, O2Idx);

220
221
        Scalar losses = 0.0;
        //Calculate activation losses
222
223
        using std::log;
        using std::abs;
Timo Koch's avatar
Timo Koch committed
224
        if(electroChemistryModel == ElectroChemistryModel::Acosta)
225
226
        {
            losses = preFactor
227
228
229
                            *(  log(abs(currentDensity)/abs(exchangeCurrentDensity_(volVars)))
                            - log(pO2/refO2PartialPressure)
                            - log(1 - sw)
230
                            );
Timo Koch's avatar
Timo Koch committed
231
232
        }
        else
233
234
        {
            losses = preFactor
235
236
            *(  log(abs(currentDensity)/abs(exchangeCurrentDensity_(volVars)))
                - log(pO2/refO2PartialPressure)
237
238
239
240
                );
        }
        return losses;
    }
Timo Koch's avatar
Timo Koch committed
241
242


243
    /*!
244
245
246
     * \brief Calculation of concentration losses.
     * \param volVars The volume variables
     */
Timo Koch's avatar
Timo Koch committed
247
    static Scalar calculateConcentrationLosses_(const VolumeVariables &volVars)
248
    {
249
250
251
        static Scalar pO2Inlet = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.pO2Inlet");
        static Scalar numElectrons = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.NumElectrons");
        static Scalar transferCoefficient =getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.TransferCoefficient");
Timo Koch's avatar
Timo Koch committed
252

253
        //Calculate preFactor
Timo Koch's avatar
Timo Koch committed
254
        Scalar preFactor = Constant::R*volVars.temperature()/transferCoefficient/Constant::F/numElectrons;
255
        //Get partial pressure of O2 in the gas phase
Timo Koch's avatar
Timo Koch committed
256
257
        Scalar pO2 = volVars.pressure(nPhaseIdx) * volVars.fluidState().moleFraction(nPhaseIdx, O2Idx);

258
259
        Scalar losses = 0.0;
        //Calculate concentration losses
260
        using std::log;
Timo Koch's avatar
Timo Koch committed
261
        if(electroChemistryModel == ElectroChemistryModel::Acosta)
262
        {
263
            losses = -1.0*preFactor*(transferCoefficient/2)*log(pO2/pO2Inlet);
264
265
        }else
        {
Timo Koch's avatar
Timo Koch committed
266
            // +1 is the Nernst part of the equation
267
            losses = -1.0*preFactor*(transferCoefficient/2+1)*log(pO2/pO2Inlet);
268
269
270
271
        }

        return losses;
    }
Timo Koch's avatar
Timo Koch committed
272
273


274
    /*!
275
276
277
     * \brief Calculation of the exchange current density.
     * \param volVars The volume variables
     */
Timo Koch's avatar
Timo Koch committed
278
    static Scalar exchangeCurrentDensity_(const VolumeVariables &volVars)
279
    {
280
        using std::exp;
281
282
283
284
        static Scalar activationBarrier =getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.ActivationBarrier");
        static Scalar surfaceIncreasingFactor = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.SurfaceIncreasingFactor");
        static Scalar refTemperature = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.RefTemperature");
        static Scalar refCurrentDensity = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "ElectroChemistry.RefCurrentDensity");
Timo Koch's avatar
Timo Koch committed
285

286
        Scalar T = volVars.fluidState().temperature();
Timo Koch's avatar
Timo Koch committed
287
288
289
        Scalar refExchangeCurrentDensity = -1.0
                            * refCurrentDensity
                            * surfaceIncreasingFactor
290
                            * exp(-1.0 * activationBarrier / Constant::R * (1/T-1/refTemperature));
Timo Koch's avatar
Timo Koch committed
291

292
293
        return refExchangeCurrentDensity;
    }
Timo Koch's avatar
Timo Koch committed
294
};
295
296
297

}// end namespace

Timo Koch's avatar
Timo Koch committed
298
#endif