biomin.hh 21 KB
Newer Older
Simon Scholz's avatar
Simon Scholz committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
// -*- 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
 * \ingroup Fluidsystems
 * \brief A fluid system with water and gas as phases and brine and CO2
 *        as components.
 */
#ifndef DUMUX_BIOMIN_FLUID_SYSTEM_HH
#define DUMUX_BIOMIN_FLUID_SYSTEM_HH

#include <dumux/common/parameters.hh>
#include <dumux/material/idealgas.hh>
#include <dumux/material/fluidsystems/base.hh>

#include <dumux/material/components/brine.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/co2.hh>
#include <dumux/material/components/tabulatedcomponent.hh>

#include <dumux/material/components/calciumion.hh>
#include "../components/urea.hh"

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

namespace Dumux {
namespace FluidSystems {
44

Simon Scholz's avatar
Simon Scholz committed
45
46
47
48
/*!
 * \ingroup Fluidsystems
 * \brief A compositional fluid with brine and carbon dioxide as
 *        components in both, the liquid and the gas (supercritical) phase,
49
 *        additional biomineralisation components (Ca and Urea) in the liquid phase
Simon Scholz's avatar
Simon Scholz committed
50
51
52
53
54
55
56
57
58
59
60
61
62
63
 *
 * This class provides acess to the Bio fluid system when no property system is used.
 * For Dumux users, using BioMinFluid<TypeTag> and the documentation therein is
 * recommended.
 *
 *  The user can provide their own material table for co2 properties.
 *  This fluidsystem is initialized as default with the tabulated version of
 *  water of the IAPWS-formulation, and the tabularized adapter to transfer
 *  this into brine.
 *  In the non-TypeTagged version, salinity information has to be provided with
 *  the init() methods.
 */
template <class Scalar,
          class CO2Table,
64
          class H2OType = Dumux::Components::TabulatedComponent<Dumux::Components::H2O<Scalar>> >
Simon Scholz's avatar
Simon Scholz committed
65
class BioMin
66
: public Base<Scalar, BioMin<Scalar, CO2Table, H2OType> >
Simon Scholz's avatar
Simon Scholz committed
67
68
69
{

    using ThisType = BioMin<Scalar, H2OType>;
70
    using Base = Dumux::FluidSystems::Base<Scalar, ThisType>;
Simon Scholz's avatar
Simon Scholz committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
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
330
331
332
    using Brine_CO2 = BinaryCoeff::Brine_CO2<Scalar, CO2Table>;

    using IdealGas = Dumux::IdealGas<Scalar>;

public:
    using CO2 = Components::CO2<Scalar, CO2Table>;
    using H2O = H2OType;
    using Ca = Components::CalciumIon<Scalar>;
    using Urea = Components::Urea<Scalar>;
    using Brine = Components::Brine<Scalar, H2O>;

    // the type of parameter cache objects. this fluid system does not
    // cache anything, so it uses Dumux::NullParameterCache
    using ParameterCache = Dumux::NullParameterCache;

    /****************************************
     * Fluid phase related static parameters
     ****************************************/
    static constexpr int numPhases = 2; // liquid and gas phases
    static constexpr int liquidPhaseIdx = 0; // index of the liquid phase
    static constexpr int gasPhaseIdx = 1; // index of the gas phase

    static constexpr int phase0Idx = liquidPhaseIdx; // index of the first phase
    static constexpr int phase1Idx = gasPhaseIdx; // index of the second phase
    /*!
     * \brief Return the human readable name of a fluid phase
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static std::string phaseName(int phaseIdx)
    {
        static std::string name[] =
        {
            "liquid",
            "gas"
        };

        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return name[phaseIdx];
    }

    /*!
     * \brief Returns whether the fluids are miscible
     */
    static constexpr bool isMiscible()
    { return true; }

    /*!
     * \brief Return whether a phase is liquid
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static constexpr bool isLiquid(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return phaseIdx != gasPhaseIdx;
    }

    /*!
     * \brief Returns true if and only if a fluid phase is assumed to
     *        be an ideal mixture.
     *
     * We define an ideal mixture as a fluid phase where the fugacity
     * coefficients of all components times the pressure of the phase
     * are independent on the fluid composition. This assumption is true
     * if Henry's law and Raoult's law apply. If you are unsure what
     * this function should return, it is safe to return false. The
     * only damage done will be (slightly) increased computation times
     * in some cases.
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static constexpr bool isIdealMixture(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return true;
    }

    /*!
     * \brief Returns true if and only if a fluid phase is assumed to
     *        be compressible.
     *
     * Compressible means that the partial derivative of the density
     * to the fluid pressure is always larger than zero.
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static constexpr bool isCompressible(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        return true;
    }

    /*!
     * \brief Returns true if and only if a fluid phase is assumed to
     *        be an ideal gas.
     *
     * \param phaseIdx The index of the fluid phase to consider
     */
    static bool isIdealGas(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

        // let the fluids decide
        if (phaseIdx == gasPhaseIdx)
            return H2O::gasIsIdeal() && CO2::gasIsIdeal();
        return false; // not a gas
    }

    /****************************************
     * Component related static parameters
     ****************************************/
    static constexpr int numComponents = 4; // H2O/brine, CO2, Ca, urea
    static constexpr int H2OIdx = 0;
    static constexpr int CO2Idx = 1;
    static constexpr int CaIdx = 2;
    static constexpr int UreaIdx = 3;

    static constexpr int BrineIdx = H2OIdx;
    static constexpr int comp0Idx = BrineIdx;
    static constexpr int comp1Idx = CO2Idx;

    /*!
     * \brief Return the human readable name of a component
     *
     * \param compIdx The index of the component to consider
     */
    static std::string componentName(int compIdx)
    {
        switch (compIdx) {
        case BrineIdx: return Brine::name();
        case CO2Idx: return "CO2";
        case CaIdx: return Ca::name();
        case UreaIdx: return Urea::name();
        default: DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
        };
    }

    /*!
     * \brief Return the molar mass of a component in \f$\mathrm{[kg/mol]}\f$.
     *
     * \param compIdx The index of the component to consider
     */
    static Scalar molarMass(int compIdx)
    {
        switch (compIdx) {
        case H2OIdx: return H2O::molarMass();
        // actually, the molar mass of brine is only needed for diffusion
        // but since solutes are accounted for seperately
        // only the molar mass of water is returned.
        case CO2Idx: return CO2::molarMass();
        case CaIdx: return Ca::molarMass();
        case UreaIdx: return Urea::molarMass();
        default: DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
        };
    }

    /****************************************
     * thermodynamic relations
     ****************************************/

    static void init()
    {
        init(/*startTemp=*/295.15, /*endTemp=*/305.15, /*tempSteps=*/10,
             /*startPressure=*/1e4, /*endPressure=*/1e6, /*pressureSteps=*/200);

    }

    static void init(Scalar startTemp, Scalar endTemp, int tempSteps,
                     Scalar startPressure, Scalar endPressure, int pressureSteps)
    {
        std::cout << "Initializing tables for the pure-water properties.\n";
        H2O::init(startTemp, endTemp, tempSteps,
                  startPressure, endPressure, pressureSteps);
    }

    /*!
     * \brief Given a phase's composition, temperature, pressure, and
     *        the partial pressures of all components, return its
     *        density \f$\mathrm{[kg/m^3]}\f$.
     *
     * \param fluidState The fluid state
     * \param phaseIdx The index of the phase
     */
    template <class FluidState>
    static Scalar density(const FluidState &fluidState,
                          const ParameterCache &paramCache,
                          int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);

        switch (phaseIdx) {
            // assume pure brine for the liquid phase.
            case liquidPhaseIdx:
                return liquidDensity_(temperature,
                                      pressure,
                                      fluidState.moleFraction(liquidPhaseIdx, CO2Idx),
                                      fluidState.moleFraction(liquidPhaseIdx, H2OIdx),
                                      fluidState.massFraction(liquidPhaseIdx, CaIdx)); //consider density effect of dissolved calcium

            // assume pure CO2 for the gas phase.
            case gasPhaseIdx:
                return CO2::gasDensity(temperature, pressure);

            default:
                DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); break;
        }
    }

    /*!
     * \brief The molar density \f$\rho_{mol,\alpha}\f$
     *   of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
     *
     * The molar density is defined by the
     * mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
     *
     * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
     */
    template <class FluidState>
    static Scalar molarDensity(const FluidState &fluidState,
                               const ParameterCache &paramCache,
                               int phaseIdx)
    {
        const Scalar temperature = fluidState.temperature(phaseIdx);
        const Scalar pressure = fluidState.pressure(phaseIdx);
        if (phaseIdx == liquidPhaseIdx)
        {
            return density(fluidState, paramCache, phaseIdx)
                   / fluidState.averageMolarMass(phaseIdx);
        }
        else if (phaseIdx == gasPhaseIdx)
        {
            // for the gas phase assume an ideal gas
            return CO2::gasMolarDensity(temperature, pressure);
        }
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }

    /*!
     * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$.
     *
     * \param temperature temperature of component in \f$\mathrm{[K]}\f$
     * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
     *
     * Equation given in: - Batzle & Wang (1992)
     *                    - cited by: Bachu & Adams (2002)
     *                    "Equations of State for basin geofluids"
     */
    template <class FluidState>
    static Scalar viscosity(const FluidState &fluidState,
                            const ParameterCache &paramCache,
                            int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);

333
        // assume brine with viscosity effect of Ca for the liquid phase.
Simon Scholz's avatar
Simon Scholz committed
334
        if (phaseIdx == liquidPhaseIdx)
335
336
337
        {
            return Brine::liquidViscosity(temperature, pressure);
        }
Simon Scholz's avatar
Simon Scholz committed
338

339
        // assume pure CO2 for the gas phase.
Simon Scholz's avatar
Simon Scholz committed
340
        else if (phaseIdx == gasPhaseIdx)
341
        {
342
            return CO2::gasViscosity(temperature, pressure);
343
        }
344

Simon Scholz's avatar
Simon Scholz committed
345
        else
346
        {
Simon Scholz's avatar
Simon Scholz committed
347
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
348
        }
Simon Scholz's avatar
Simon Scholz committed
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388

    }

    /*!
     * \brief Returns the fugacity coefficient [Pa] of a component in a
     *        phase.
     *
     * The fugacity coefficient \f$\phi^\kappa_\alpha\f$ of a
     * component \f$\kappa\f$ for a fluid phase \f$\alpha\f$ defines
     * the fugacity \f$f^\kappa_\alpha\f$ by the equation
     *
     * \f[
     f^\kappa_\alpha := \phi^\kappa_\alpha x^\kappa_\alpha p_\alpha\;.
     \f]
     *
     * The fugacity itself is just an other way to express the
     * chemical potential \f$\zeta^\kappa_\alpha\f$ of the component:
     *
     * \f[
     f^\kappa_\alpha := \exp\left\{\frac{\zeta^\kappa_\alpha}{k_B T_\alpha} \right\}
     \f]
     * where \f$k_B\f$ is Boltzmann's constant.
     */
    template <class FluidState>
    static Scalar fugacityCoefficient(const FluidState &fluidState,
                                      const ParameterCache &paramCache,
                                      int phaseIdx,
                                      int compIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        assert(0 <= compIdx && compIdx < numComponents);

        if (phaseIdx == gasPhaseIdx)
            // use the fugacity coefficients of an ideal gas. the
            // actual value of the fugacity is not relevant, as long
            // as the relative fluid compositions are observed,
            return 1.0;

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);
389
        Scalar salinity = Brine::salinity(); // 0.1; //TODO major assumption in favor of runtime!
Simon Scholz's avatar
Simon Scholz committed
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
449
450
451
452
453
        //function is actually designed for use with NaCl not Ca.
        //Theoretically it should be: fluidState.massFraction(liquidPhaseIdx, CaIdx);

        assert(temperature > 0);
        assert(pressure > 0);
        // calulate the equilibrium composition for the given
        // temperature and pressure.
        Scalar xwH2O, xnH2O;
        Scalar xwCO2, xnCO2;
        Brine_CO2::calculateMoleFractions(temperature,
                                          pressure,
                                          salinity,
                                          /*knowgasPhaseIdx=*/-1,
                                          xwCO2,
                                          xnH2O);

        // normalize the phase compositions
        using std::min;
        using std::max;
        xwCO2 = max(0.0, min(1.0, xwCO2));
        xnH2O = max(0.0, min(1.0, xnH2O));

        xwH2O = 1.0 - xwCO2;
        xnCO2 = 1.0 - xnH2O;

        if (compIdx == BrineIdx)
        {
            Scalar phigH2O = 1.0;
            return phigH2O * xnH2O / xwH2O;
        }
        if (compIdx == CO2Idx)
        {
            Scalar phigCO2 = 1.0;
            return phigCO2 * xnCO2 / xwCO2;
        }
        else
            return 1/pressure; //all other components stay in the liquid phase
    }

    /*!
     * \brief Given the phase compositions, return the binary
     *        diffusion coefficient \f$\mathrm{[m^2/s]}\f$ of two components in a phase.
     * \param fluidState An arbitrary fluid state
     * \param phaseIdx The index of the fluid phase to consider
     * \param compIIdx Index of the component i
     * \param compJIdx Index of the component j
     */
    template <class FluidState>
    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
                                             const ParameterCache &paramCache,
                                             int phaseIdx,
                                             int compIIdx,
                                             int compJIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        assert(0 <= compIIdx && compIIdx < numComponents);
        assert(0 <= compJIdx && compJIdx < numComponents);

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);

        if (phaseIdx == liquidPhaseIdx)
        {
            assert(compIIdx == H2OIdx);
454
455
            if (compJIdx == CO2Idx)
                return Brine_CO2::liquidDiffCoeff(temperature, pressure);
Simon Scholz's avatar
Simon Scholz committed
456
            else if (compJIdx < numComponents) //Calcium and urea
457
                return 1.587e-9;  //[m²/s]  educated guess, value for NaCl from J. Phys. D: Appl. Phys. 40 (2007) 2769-2776
Simon Scholz's avatar
Simon Scholz committed
458
459
460
461
462
463
464
            else
                DUNE_THROW(Dune::NotImplemented, "Binary difussion coefficient : Incorrect compIdx");
        }
        else
        {
            assert(phaseIdx == gasPhaseIdx);
            assert(compIIdx == CO2Idx);
465
466
            if (compJIdx == H2OIdx)
                return Brine_CO2::gasDiffCoeff(temperature, pressure);
Simon Scholz's avatar
Simon Scholz committed
467
            else if (compJIdx <numComponents) //Calcium and urea will stay in brine, no gaseous calcium or urea!
468
                return 0.0;
Simon Scholz's avatar
Simon Scholz committed
469
470
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
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
            else
                DUNE_THROW(Dune::NotImplemented, "Binary difussion coefficient : Incorrect compIdx");
        }
    };

    /*!
     * \brief Given a phase's composition, temperature and pressure,
     *        return its specific enthalpy \f$\mathrm{[J/kg]}\f$.
     * \param fluidState An arbitrary fluid state
     * \param phaseIdx The index of the fluid phase to consider
     *
     * See:
     * Class 2001:
     * Theorie und numerische Modellierung nichtisothermer Mehrphasenprozesse in NAPL-kontaminierten porösen Medien
     * Chapter 2.1.13 Innere Energie, Wäremekapazität, Enthalpie \cite A3:class:2001 <BR>
     *
     * Formula (2.42):
     * the specific enthalpy of a gasphase result from the sum of (enthalpies*mass fraction) of the components
     *
     */
    template <class FluidState>
    static Scalar enthalpy(const FluidState &fluidState,
                           const ParameterCache &paramCache,
                           int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

        Scalar temperature = fluidState.temperature(phaseIdx);
        Scalar pressure = fluidState.pressure(phaseIdx);

        if (phaseIdx == liquidPhaseIdx)
        {
            // assume pure brine for the liquid phase.
            return Brine::liquidEnthalpy(temperature, pressure);
        }
        else
        {
            // assume pure CO2 for the gas phase.
            return CO2::gasEnthalpy(temperature, pressure);
        }
    };

private:
    //! calculate liquid density with respect to Water, CO2 and salt
    static Scalar liquidDensity_(Scalar T,
                                 Scalar pl,
                                 Scalar xwCO2,
                                 Scalar xwH2O,
                                 Scalar XlSal)
    {
        if(T < 273.15)
        {
            DUNE_THROW(NumericalProblem,
                       "Liquid density for Brine and CO2 is only "
                       "defined above 273.15K (is" << T << ")");
        }
        if(pl >= 2.5e8)
        {
            DUNE_THROW(NumericalProblem,
                       "Liquid density for Brine and CO2 is only "
                       "defined below 250MPa (is" << pl << ")");
        }

532
533
534
535
        const Scalar rho_brine = Brine::liquidDensity(T, pl);
        const Scalar rho_pure = H2O::liquidDensity(T, pl);
        const Scalar rho_lCO2 = liquidDensityWaterCO2_(T, pl, xwH2O, xwCO2);
        const Scalar contribCO2 = rho_lCO2 - rho_pure;
Simon Scholz's avatar
Simon Scholz committed
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
        return rho_brine + contribCO2;
    }

    //! calculate liquid Density of water and CO2
    static Scalar liquidDensityWaterCO2_(Scalar temperature,
                                         Scalar pl,
                                         Scalar xwH2O,
                                         Scalar xwCO2)
    {
        const Scalar M_CO2 = CO2::molarMass();
        const Scalar M_H2O = H2O::molarMass();

        const Scalar tempC = temperature - 273.15;        /* tempC : temperature in °C */
        const Scalar rho_pure = H2O::liquidDensity(temperature, pl);
        xwH2O = 1.0 - xwCO2; // xwH2O is available, but in case of a pure gas phase
                             // the value of M_T for the virtual liquid phase can become very large
        const Scalar M_T = M_H2O * xwH2O + M_CO2 * xwCO2;
        const Scalar V_phi =
            (37.51 +
             tempC*(-9.585e-2 +
                    tempC*(8.74e-4 -
                           tempC*5.044e-7))) / 1.0e6;
        return 1 / (xwCO2 * V_phi/M_T + M_H2O * xwH2O / (rho_pure * M_T));
    }

};

} // end namespace FluidSystems
} // end namespace Dumux

#endif