h2o_n2_system.hh 19.1 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
// $Id$
2
3
4
5
6
7
/*****************************************************************************
 *   Copyright (C) 2009 by Andreas Lauser
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
8
 *   This program is free software: you can redistribute it and/or modify    *
9
 *   it under the terms of the GNU General Public License as published by    *
10
11
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
12
 *                                                                           *
13
14
15
16
17
18
19
 *   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/>.   *
20
 *****************************************************************************/
21
/*!
Bernd Flemisch's avatar
Bernd Flemisch committed
22
 * \file
23
 *
24
 * \ingroup Fluidsystems
Bernd Flemisch's avatar
Bernd Flemisch committed
25
 * \brief A fluid system with a liquid and a gaseous phase and \f$H_2O\f$ and \f$N_2\f$
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
 *        as components.
 */
#ifndef DUMUX_H2O_N2_SYSTEM_HH
#define DUMUX_H2O_N2_SYSTEM_HH

#include <dumux/material/idealgas.hh>
#include <dumux/material/settablephase.hh>

#include <dumux/material/binarycoefficients/h2o_n2.hh>

#include "defaultcomponents.hh"

#include <dumux/common/propertysystem.hh>

namespace Dumux
{
Andreas Lauser's avatar
Andreas Lauser committed
42
43
44
45
46
47
// forward defintions of the property tags
namespace Properties {
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(Components);
};

48
/*!
49
50
 * \ingroup Fluidsystems
 *
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
 * \brief A compositional fluid with water and molecular nitrogen as
 *        components in both, the liquid and the gas phase.
 */
template <class TypeTag>
class H2O_N2_System
{
    typedef H2O_N2_System<TypeTag> ThisType;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;

    typedef Dumux::IdealGas<Scalar> IdealGas;
    typedef Dumux::SettablePhase<Scalar, ThisType> SettablePhase;

    typedef typename GET_PROP(TypeTag, PTAG(Components)) Components;

public:
    typedef typename Components::H2O H2O;
    typedef typename Components::N2 N2;

    static const int numComponents = 2;
    static const int numPhases = 2;

    static const int lPhaseIdx = 0; // index of the liquid phase
    static const int gPhaseIdx = 1; // index of the gas phase

    static const int wPhaseIdx = lPhaseIdx; // index of the wetting phase
    static const int nPhaseIdx = gPhaseIdx; // index of the non-wetting phase

    static const int H2OIdx = 0;
    static const int N2Idx = 1;

    static void init()
    { Components::init(); }

84
85
    /*!
     * \brief Return the human readable name of a phase
86
87
     *
     * \param phaseIdx index of the phase
88
89
90
91
92
93
94
95
96
97
     */
    static const char *phaseName(int phaseIdx)
    {
        switch (phaseIdx) {
        case lPhaseIdx: return "l";
        case gPhaseIdx: return "g";
        };
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }

98
99
    /*!
     * \brief Return the human readable name of a component
100
101
     *
     * \param compIdx index of the component
102
103
104
105
106
107
108
109
110
111
112
113
     */
    static const char *componentName(int compIdx)
    {
        switch (compIdx) {
        case H2OIdx: return H2O::name();
        case N2Idx: return N2::name();
        };
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }

    /*!
     * \brief Return the molar mass of a component [kg/mol].
114
115
     *
     * \param compIdx index of the component
116
117
118
119
120
121
122
123
124
125
126
127
128
129
     */
    static Scalar molarMass(int compIdx)
    {
        switch (compIdx) {
        case H2OIdx: return H2O::molarMass();
        case N2Idx: return N2::molarMass();
        };
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }

    /*!
     * \brief Given a phase's composition, temperature, pressure, and
     *        the partial pressures of all components, return its
     *        density [kg/m^3].
130
131
132
133
134
135
136
     *
     * \param phaseIdx index of the phase
     * \param temperature phase temperature in [K]
     * \param pressure phase pressure in [Pa]
     * \param fluidState the fluid state
     *
     * \tparam FluidState the fluid state class
137
138
139
140
141
142
143
     */
    template <class FluidState>
    static Scalar phaseDensity(int phaseIdx,
                               Scalar temperature,
                               Scalar pressure,
                               const FluidState &fluidState)
    {
Andreas Lauser's avatar
Andreas Lauser committed
144
145
146
147
148
        if (phaseIdx == lPhaseIdx) {
            // See: Ochs 2008
            // \todo: proper citation
            Scalar rholH2O = H2O::liquidDensity(temperature, pressure);
            Scalar clH2O = rholH2O/H2O::molarMass();
149

Andreas Lauser's avatar
Andreas Lauser committed
150
151
152
153
            // this assumes each nitrogen molecule displaces exactly one
            // water molecule in the liquid
            return
                clH2O*(H2O::molarMass()*fluidState.moleFrac(lPhaseIdx, H2OIdx)
154
                       +
Andreas Lauser's avatar
Andreas Lauser committed
155
156
157
                       N2::molarMass()*fluidState.moleFrac(lPhaseIdx, N2Idx));
        }
        else if (phaseIdx == gPhaseIdx) {
158
            Scalar fugH2O =
159
160
                fluidState.moleFrac(gPhaseIdx, H2OIdx)  *
                fluidState.phasePressure(gPhaseIdx);
161
            Scalar fugN2 =
162
163
                fluidState.moleFrac(gPhaseIdx, N2Idx)  *
                fluidState.phasePressure(gPhaseIdx);
164
            return
165
                H2O::gasDensity(temperature, fugH2O) +
Andreas Lauser's avatar
Andreas Lauser committed
166
167
                N2::gasDensity(temperature, fugN2);
        }
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }

    /*!
     * \brief Given a phase's composition, temperature and pressure,
     *        return its viscosity.
     */
    template <class FluidState>
    static Scalar phaseViscosity(int phaseIdx,
                                 Scalar temperature,
                                 Scalar pressure,
                                 const FluidState &fluidState)
    {
        if (phaseIdx == lPhaseIdx) {
            // assume pure water for the liquid phase
            // TODO: viscosity of mixture
            return H2O::liquidViscosity(temperature,
                                        pressure);
        }
        else {
            return N2::gasViscosity(temperature,
                                    pressure);
            /* Wilke method. See:
             *
             * S.O.Ochs: "Development of a multiphase multicomponent
             * model for PEMFC - Technical report: IRTG-NUPUS",
             * University of Stuttgart, 2008
             *
             * and:
             *
             * See: R. Reid, et al.: The Properties of Gases and Liquids, 4th
             * edition, McGraw-Hill, 1987, 407-410
             */
            Scalar muResult = 0;
            const Scalar mu[numComponents] = {
                H2O::gasViscosity(temperature,
                                  H2O::vaporPressure(temperature)),
                N2::gasViscosity(temperature,
                                 pressure)
            };
            // molar masses
            const Scalar M[numComponents] = {
                H2O::molarMass(),
                N2::molarMass()
            };

            for (int i = 0; i < numComponents; ++i) {
                Scalar divisor = 0;
                for (int j = 0; j < numComponents; ++j) {
                    Scalar phiIJ = 1 + sqrt(mu[i]/mu[j] *
                                            pow(M[i]/M[j], 1/4.0));
                    phiIJ *= phiIJ;
                    phiIJ /= sqrt(8*(1 + M[i]/M[j]));
                    divisor += fluidState.moleFrac(phaseIdx, j)*phiIJ;
                }
                muResult += fluidState.moleFrac(phaseIdx, i)*mu[i] / divisor;
            }

            return muResult;
        }
    }


    /*!
     * \brief Assuming the composition of a single phase and the
     *        pressure of all phases is known or that all phases are
     *        present, compute the thermodynamic equilibrium from the
     *        temperature and phase pressures. If the known phase
     *        index
     *
     */
    template <class FluidState>
    static void computeEquilibrium(FluidState &fluidState,
                                   int knownPhaseIdx = -1)
    {
        const Scalar T = fluidState.temperature();
        const Scalar pg = fluidState.phasePressure(gPhaseIdx);
        const Scalar pl = fluidState.phasePressure(lPhaseIdx);

        const Scalar betaH2O = H2O::vaporPressure(T);
        const Scalar betaN2 = BinaryCoeff::H2O_N2::henry(T);

        if (knownPhaseIdx < 0)
        {
            // we only have all phase pressures and temperature and
            // know that all phases are present
            Scalar xlH2O = (pg - betaN2)/(betaH2O - betaN2);
            Scalar xlN2 = 1 - xlH2O;
            Scalar rhol = liquidPhaseDensity_(T, pl, xlH2O, xlN2);

            Scalar cgH2O = H2O::gasDensity(T, betaH2O*xlH2O)/H2O::molarMass();
            Scalar cgN2 = N2::gasDensity(T, betaN2*xlN2)/N2::molarMass();

            Scalar xgH2O = cgH2O/(cgH2O + cgN2);
            Scalar xgN2 = cgN2/(cgH2O + cgN2);

            // set the liquid phase composition
            SettablePhase liquid;
            liquid.moleFrac_[H2OIdx] = xlH2O;
            liquid.moleFrac_[N2Idx] = xlN2;
            liquid.pressure_ = pl;
            liquid.density_ = rhol;
            liquid.xToX(); // compute mass fractions from mole fractions
            fluidState.assignPhase(lPhaseIdx, liquid);

            // set the gas phase composition
            SettablePhase gas;
            gas.moleFrac_[H2OIdx] = xgH2O;
            gas.moleFrac_[N2Idx] = xgN2;
            gas.pressure_ = pg;
            gas.density_ = cgH2O*H2O::molarMass() + cgN2*N2::molarMass();
            gas.xToX(); // compute mass fractions from mole fractions
            fluidState.assignPhase(gPhaseIdx, gas);
        }
        else if (knownPhaseIdx == lPhaseIdx) {
            // the composition of the liquid phase is given

            // retrieve the known mole fractions from the fluid state
            Scalar xlH2O = fluidState.moleFrac(lPhaseIdx, H2OIdx);
            Scalar xlN2 = fluidState.moleFrac(lPhaseIdx, N2Idx);

            // calculate the component contentrations in the gas phase
            Scalar pH2O = betaH2O*xlH2O; // fugacity of water
            Scalar pN2 = betaN2*xlN2; // fugacity of nitrogen
            Scalar cgH2O = H2O::gasDensity(T, pH2O)/H2O::molarMass();
            Scalar cgN2 = N2::gasDensity(T, pN2)/N2::molarMass();

            // convert concentrations to mole fractions
            Scalar xgH2O = cgH2O/(cgH2O + cgN2) * (pH2O + pN2)/pg;
            Scalar xgN2 = cgN2/(cgH2O + cgN2) * (pH2O + pN2)/pg;

            // set gas phase composition
            SettablePhase gas;
            gas.moleFrac_[H2OIdx] = xgH2O;
            gas.moleFrac_[N2Idx] = xgN2;
            gas.pressure_ = pg;
            gas.density_ = cgH2O*H2O::molarMass() + cgN2*N2::molarMass();
            gas.xToX(); // update mass fractions from mole fractions
            fluidState.assignPhase(gPhaseIdx, gas);
        }
        else if (knownPhaseIdx == gPhaseIdx) {
            // the composition of the gas phase is given

            Scalar xgH2O = fluidState.moleFrac(gPhaseIdx, H2OIdx);
            Scalar xgN2 = fluidState.moleFrac(gPhaseIdx, N2Idx);
            Scalar pgH2O = pg*xgH2O;
            Scalar pgN2 = pg*xgN2;

            Scalar xlH2O = pgH2O/betaH2O;
            Scalar xlN2 = pgN2/betaN2;

            SettablePhase liquid;
            liquid.moleFrac_[H2OIdx] = xlH2O;
            liquid.moleFrac_[N2Idx] = xlN2;
            liquid.pressure_ = pl;
            liquid.density_ = liquidPhaseDensity_(T, pl, xlH2O, xlN2);
            liquid.xToX(); // update mass fractions from mole fractions
            fluidState.assignPhase(lPhaseIdx, liquid);
        }
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
330
     * \brief Returns the activity coefficient of a component in a
331
332
     *        phase.
     *
Andreas Lauser's avatar
Andreas Lauser committed
333
334
335
336
     * We define the activity coefficent \f$\gamma_{\alpha,\kappa}\f$
     * of component \f$\kappa\f$ by the following equation:
     *  \f[ f_\kappa = p_\alpha \gamma_{\alpha,\kappa} \f]
     * where \f$f_\kappa\f$  is the component's fugacity and \f$p_\alpha\f$
337
338
339
     * is the phase' pressure. To yield a phase-independent name for this
     * method, it also treats gaseous phases although the term "activity"
     * might be misleading in that context.
Andreas Lauser's avatar
Andreas Lauser committed
340
341
342
343
344
345
346
347
     *
     * For liquids with very low miscibility this boils down to the
     * inverse Henry constant for the solutes and the partial pressure
     * for the solvent.
     *
     * For ideal gases this is equivalent to the gas pressure, in real
     * gases it is the gas pressure times the component's fugacity
     * coefficient.
348
349
     */
    template <class FluidState>
Andreas Lauser's avatar
Andreas Lauser committed
350
    static Scalar activityCoeff(int phaseIdx,
351
352
353
354
355
                                int compIdx,
                                Scalar temperature,
                                Scalar pressure,
                                const FluidState &state)
    {
Andreas Lauser's avatar
Andreas Lauser committed
356
357
358
359
360
361
        if (phaseIdx == gPhaseIdx) {
            return pressure;
            Scalar fugH2O = std::max(1e-3, state.fugacity(H2OIdx));
            Scalar fugN2 = std::max(1e-3, state.fugacity(N2Idx));
            Scalar cH2O = H2O::gasDensity(temperature, fugH2O) / H2O::molarMass();
            Scalar cN2 = N2::gasDensity(temperature, fugN2) / N2::molarMass();
362

Andreas Lauser's avatar
Andreas Lauser committed
363
364
365
366
367
368
            Scalar alpha = (fugH2O + fugN2)/pressure;

            if (compIdx == H2OIdx)
                return fugH2O/(alpha*cH2O/(cH2O + cN2));
            else if (compIdx == N2Idx)
                return fugN2/(alpha*cN2/(cH2O + cN2));
369

Andreas Lauser's avatar
Andreas Lauser committed
370
371
            DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
        }
372
373
374
375
376
377
378
379
380
381

        switch (compIdx) {
        case H2OIdx: return H2O::vaporPressure(temperature);
        case N2Idx: return BinaryCoeff::H2O_N2::henry(temperature);
        };
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }

    /*!
     * \brief Given a phase's composition, temperature and pressure,
382
     *        return the binary diffusion coefficient for components
383
     *        \f$i\f$ and \f$j\f$ in this phase.
384
     *
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
     */
    template <class FluidState>
    static Scalar diffCoeff(int phaseIdx,
                            int compIIdx,
                            int compJIdx,
                            Scalar temperature,
                            Scalar pressure,
                            const FluidState &fluidState)
    {
        if (compIIdx > compJIdx)
            std::swap(compIIdx, compJIdx);

#ifndef NDEBUG
        if (compIIdx == compJIdx ||
            phaseIdx > numPhases - 1 ||
            compJIdx > numComponents - 1)
        {
            DUNE_THROW(Dune::InvalidStateException,
                       "Binary diffusion coefficient of components "
                       << compIIdx << " and " << compJIdx
                       << " in phase " << phaseIdx << " is undefined!\n");
        }
#endif


        switch (phaseIdx) {
        case lPhaseIdx:
            switch (compIIdx) {
            case H2OIdx:
                switch (compJIdx) {
                case N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(temperature,
                                                                        pressure);
                }
            default:
                DUNE_THROW(Dune::InvalidStateException,
                           "Binary diffusion coefficients of trace "
                           "substances in liquid phase is undefined!\n");
            }
        case gPhaseIdx:
            switch (compIIdx) {
            case H2OIdx:
                switch (compJIdx) {
                case N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(temperature,
                                                                     pressure);
                }
            }
        }

        DUNE_THROW(Dune::InvalidStateException,
                   "Binary diffusion coefficient of components "
                   << compIIdx << " and " << compJIdx
                   << " in phase " << phaseIdx << " is undefined!\n");
    };

    /*!
     * \brief Given a phase's composition, temperature and pressure,
     *        return its specific enthalpy [J/kg].
     */
    template <class FluidState>
    static Scalar phaseEnthalpy(int phaseIdx,
                                Scalar temperature,
                                Scalar pressure,
                                const FluidState &fluidState)
    {
449
        if (phaseIdx == lPhaseIdx) {
450
451
452
453
454
455
456
            Scalar cN2 = fluidState.concentration(lPhaseIdx, N2Idx);
            Scalar pN2 = N2::gasPressure(temperature, cN2*N2::molarMass());

            // TODO: correct way to deal with the solutes???
            return
                fluidState.massFrac(lPhaseIdx, H2OIdx)*
                H2O::liquidEnthalpy(temperature, pressure)
457
                +
458
459
460
461
                fluidState.massFrac(lPhaseIdx, N2Idx)*
                N2::gasEnthalpy(temperature, pN2);
        }
        else {
462
463
            Scalar cH2O = fluidState.concentration(gPhaseIdx, H2OIdx);
            Scalar cN2 = fluidState.concentration(gPhaseIdx, N2Idx);
464

465
466
            Scalar pH2O = H2O::gasPressure(temperature, cH2O*H2O::molarMass());
            Scalar pN2 = N2::gasPressure(temperature, cN2*N2::molarMass());
467
468
469

            Scalar result = 0;
            result +=
470
                H2O::gasEnthalpy(temperature, pH2O) *
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
                fluidState.massFrac(gPhaseIdx, H2OIdx);
            result +=
                N2::gasEnthalpy(temperature, pN2) *
                fluidState.massFrac(gPhaseIdx, N2Idx);

            return result;
        }
    }

    /*!
     * \brief Given a phase's composition, temperature and pressure,
     *        return its specific internal energy [J/kg].
     */
    template <class FluidState>
    static Scalar phaseInternalEnergy(int phaseIdx,
                                      Scalar temperature,
                                      Scalar pressure,
                                      const FluidState &fluidState)
    {
        return
            phaseEnthalpy(phaseIdx, temperature, pressure, fluidState) -
            pressure/phaseDensity(phaseIdx, temperature, pressure, fluidState);
    }

private:
    static Scalar liquidPhaseDensity_(Scalar T, Scalar pl, Scalar xlH2O, Scalar xlN2)
    {
        // See: Ochs 2008
        // \todo: proper citation
        Scalar rholH2O = H2O::liquidDensity(T, pl);
        Scalar clH2O = rholH2O/H2O::molarMass();
502

503
504
505
506
        // this assumes each nitrogen molecule displaces exactly one
        // water molecule in the liquid
        return
            clH2O*(xlH2O*H2O::molarMass()
507
                   +
508
509
510
511
512
                   xlN2*N2::molarMass());
    }

    static Scalar gasPhaseDensity_(Scalar T, Scalar pg, Scalar xgH2O, Scalar xgN2)
    {
Andreas Lauser's avatar
Andreas Lauser committed
513
        return H2O::gasDensity(T, pg*xgH2O) + N2::gasDensity(T, pg*xgN2);
514
515
516
517
518
519
    };
};

} // end namepace

#endif