steamn2cao2h2.hh 22 KB
Newer Older
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
// -*- 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
 *
 * \brief @copybrief Dumux::FluidSystems::SteamN2CaO2H2
 */
#ifndef DUMUX_STEAM_N2_CAO2H2_SYSTEM_HH
#define DUMUX_STEAM_N2_CAO2H2_SYSTEM_HH

#include <cassert>

#include <dune/common/exceptions.hh>

#include <dumux/common/valgrind.hh>

#include <dumux/material/idealgas.hh>
#include <dumux/material/fluidsystems/base.hh>
#include <dumux/material/components/n2.hh>
#include <dumux/material/components/h2o.hh>
// #include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/CaO2H2.hh>
39
#include <dumux/material/components/CaOtest.hh>
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include <dumux/material/binarycoefficients/h2o_n2.hh>
#include <dumux/material/components/tabulatedcomponent.hh>

namespace Dumux
{
namespace FluidSystems
{
/*!
 * \ingroup Fluidsystems
 *
 * \brief A compositional one-phase fluid system with \f$H_2O\f$ \f$N_2\f$ as gaseous components
 *              and \f$CaO\f$  and \f$Ca(OH)_2/f$ as solid components drawn for thermo-chemical
 *              heat storage.
 *
54
55
56
 *  This fluidsystem is applied by default with the tabulated version of
 *  water of the IAPWS-formulation. However, the IAPWS-formulation has to be
 *  adapted, if to higher temperatures and higher pressures occur.
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
 */

template <class Scalar,
           class H2Otype = Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar>>,
          bool useComplexRelations=true>
class SteamN2CaO2H2
: public BaseFluidSystem<Scalar, SteamN2CaO2H2<Scalar, H2Otype, useComplexRelations> >
{
    typedef SteamN2CaO2H2<Scalar, H2Otype, useComplexRelations> ThisType;
    typedef BaseFluidSystem <Scalar, ThisType> Base;

    typedef Dumux::IdealGas<Scalar> IdealGas;

public:
//     typedef Dumux::SimpleH2O<Scalar> H2O;
    typedef H2Otype H2O;
    typedef Dumux::BinaryCoeff::H2O_N2 H2O_N2;
    typedef Dumux::N2<Scalar> N2;

76
    typedef Dumux::CaOTest<Scalar> CaO;
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
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
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
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
449
    typedef Dumux::CaO2H2<Scalar> CaO2H2;

    // the type of parameter cache objects. this fluid system does not
    typedef Dumux::NullParameterCache ParameterCache;

    /****************************************
     * Fluid phase related static parameters
     ****************************************/

    //! Number of phases in the fluid system
    static constexpr int numPhases = 1; //gas phase: N2 and steam
    static const int numSPhases = 2;//  solid phases CaO and CaO2H2

    static constexpr int gPhaseIdx = 0;
//     static constexpr int gPhaseIdx = phaseIdx;
    static const int nPhaseIdx = gPhaseIdx; // index of the gas phase

    static constexpr int cPhaseIdx = 1; // CaO-phaseIdx
    static constexpr int hPhaseIdx = 2; // CaO2H2-phaseIdx


    /*!
     * \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)
    {
        switch (phaseIdx) {
        case nPhaseIdx: return "gas";
        case cPhaseIdx : return "CaO";
        case hPhaseIdx : return "CaOH2";
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }

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

    /*!
     * \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 Rault'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 bool isIdealMixture(int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);
        // we assume no interaction between gas molecules of different
        // components

        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 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 fluid decide
        return H2O::gasIsIdeal() && N2::gasIsIdeal();
    }

     /****************************************
     * Component related static parameters
     ****************************************/

    static const int numComponents = 2;//3; // H2O, Air
    static const int numMajorComponents = 2;// H2O, Air
    static const int numSComponents = 2;// CaO2H2, CaO


    static const int N2Idx = 0;
    static const int H2OIdx = 1;
    static const int CaOIdx  = 2;
    static const int CaO2H2Idx  = 3;


     /*!
     * \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 H2OIdx: return "H2O";
        case N2Idx: return "N2";
        case CaOIdx:return "CaO";
        case CaO2H2Idx:return "CaO2H2";

        }
        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();
        case N2Idx: return N2::molarMass();
        case CaOIdx: return CaO::molarMass();
        case CaO2H2Idx: return CaO2H2::molarMass();
        }
        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
    }

        /*!
     * \brief Return the mass density of the solid \f$\mathrm{[kg/m^3]}\f$.
     *
     * \param phaseIdx The index of the solid phase to consider
     */
    static Scalar precipitateDensity(int phaseIdx)
    {
      if(phaseIdx==cPhaseIdx)
            return CaO::density();
      if(phaseIdx==hPhaseIdx)
            return CaO2H2::density();
      else
        DUNE_THROW(Dune::InvalidStateException, "Invalid solid index " << phaseIdx);
        return 1;
    }

     /*!
     * \brief Return the salt specific heat capacity \f$\mathrm{[J/molK]}\f$.
     *
     * \param phaseIdx The index of the solid phase to consider
     */

     static Scalar precipitateHeatCapacity(int phaseIdx)
     {
      if(phaseIdx==cPhaseIdx)
            return CaO::heatCapacity();
      if(phaseIdx==hPhaseIdx)
            return CaO2H2::heatCapacity();
      else
      DUNE_THROW(Dune::InvalidStateException, "Invalid solid index " << phaseIdx);
     return 1;
    }

    /*!
     * \brief Return the molar density of the solid \f$\mathrm{[mol/m^3]}\f$.
     *
     * \param phaseIdx The index of the solid phase to consider
     */
    static Scalar precipitateMolarDensity(int phaseIdx)
     {
     if(phaseIdx==1){
         return precipitateDensity(phaseIdx)/ molarMass(CaOIdx);
     }
      if(phaseIdx==2){
            return precipitateDensity(phaseIdx)/molarMass(CaO2H2Idx);  }
      else
         DUNE_THROW(Dune::InvalidStateException, "Invalid solid phase index " << phaseIdx);
     }

    /****************************************
     * thermodynamic relations
     ****************************************/
    /*!
     *\brief Initialize the fluid system's static parameters generically
     *
     * If a tabulated H2O component is used, we do our best to create
     * tables that always work.
     */
    static void init()
    {
        init(/*tempMin=*/473.15,
             /*tempMax=*/723.0,
             /*numTemptempSteps=*/25,
             /*startPressure=*/0,
             /*endPressure=*/9e6,
             /*pressureSteps=*/200);
    }

   /*!
    * \brief Initialize the fluid system's static parameters using
    *        problem specific temperature and pressure ranges
    *
    * \param tempMin The minimum temperature used for tabulation of water \f$\mathrm{[K]}\f$
    * \param tempMax The maximum temperature used for tabulation of water \f$\mathrm{[K]}\f$
    * \param nTemp The number of ticks on the temperature axis of the  table of water
    * \param pressMin The minimum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
    * \param pressMax The maximum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$
    * \param nPress The number of ticks on the pressure axis of the  table of water
    */

    static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
                      Scalar pressMin, Scalar pressMax, unsigned nPress)
    {
        if (useComplexRelations)
            std::cout << "Using complex H2O-N2-CaO2H2 fluid system\n";
        else
            std::cout << "Using fast H2O-N2-CaO2H2 fluid system\n";

        if (H2O::isTabulated) {
            std::cout << "Initializing tables for the H2O fluid properties ("
                        << nTemp*nPress
                        << " entries).\n";

            H2O::init(tempMin, tempMax, nTemp,
                                pressMin, pressMax, nPress);
        }
    }

    /*!
     * \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 phaseIdx index of the phase
     * \param temperature phase temperature in \f$\mathrm{[K]}\f$
     * \param pressure phase pressure in \f$\mathrm{[Pa]}\f$
     * \param fluidState the fluid state
     *
     * Equation given in:
     * - Batzle & Wang (1992) \cite batzle1992
     * - cited by: Bachu & Adams (2002)
     *   "Equations of State for basin geofluids" \cite adams2002
     */
    using Base::density;
    template <class FluidState>
    static Scalar density(const FluidState &fluidState,
                          int phaseIdx)
    {
        assert(0 <= phaseIdx && phaseIdx < numPhases);

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

        Scalar sumMoleFrac = 0;
        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
            sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);

        if (!useComplexRelations)
            // for the gas phase assume an ideal gas
            return
                IdealGas::molarDensity(T, p)
                * fluidState.averageMolarMass(nPhaseIdx)
                / std::max(1e-5, sumMoleFrac);
        else
        {
            return
                (H2O::gasDensity(T, fluidState.partialPressure(nPhaseIdx, H2OIdx)) +
                N2::gasDensity(T, fluidState.partialPressure(nPhaseIdx, N2Idx)));
        }
      }

    /*!
     * \brief Calculate the dynamic viscosity of a fluid phase \f$\mathrm{[Pa*s]}\f$
     *
     * \param fluidState An arbitrary fluid state
     * \param phaseIdx The index of the fluid phase to consider
     *
     * \note For the viscosity of the phases the contribution of the minor
     *       component is neglected. This contribution is probably not big, but somebody
     *       would have to find out its influence.
     */
    using Base::viscosity;
    template <class FluidState>
    static Scalar viscosity(const FluidState &fluidState,
                            int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);

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

        // Wilke method (Reid et al.):
        Scalar muResult = 0;
        const Scalar mu[numComponents] = {
            H2O::gasViscosity(T, H2O::vaporPressure(T)),
            N2::gasViscosity(T, p)
        };

        const Scalar M[numComponents] = {
            H2O::molarMass(),
            N2::molarMass()
        };

        Scalar sumx = 0.0;
        for (int compIdx = 0; compIdx < 2; ++compIdx)
            sumx += fluidState.moleFraction(phaseIdx, compIdx);

        sumx = std::max(1e-10, sumx);

        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[j]/M[i], 1/4.0);
                phiIJ *= phiIJ;
                phiIJ /= sqrt(8*(1 + M[i]/M[j]));
                divisor += fluidState.moleFraction(phaseIdx, j)/sumx * phiIJ;
            }
            muResult += fluidState.moleFraction(phaseIdx, i)/sumx * mu[i] / divisor;
        }
        return muResult;
    }

    /*!
     * \brief Given a phase's composition, temperature and pressure,
     *        return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components
     *        \f$\mathrm{i}\f$ and \f$\mathrm{j}\f$ in this phase.
     *
     * \param fluidState An arbitrary fluid state
     * \param phaseIdx The index of the fluid phase to consider
     * \param compIIdx The index of the first component to consider
     * \param compJIdx The index of the second component to consider
     */
    using Base::binaryDiffusionCoefficient;
    template <class FluidState>
    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
                                             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);

        assert(phaseIdx == gPhaseIdx);

        if (compIIdx != N2Idx)
            std::swap(compIIdx, compJIdx);

        Scalar result = 0.0;
        if(compJIdx == H2OIdx)
        result = H2O_N2::gasDiffCoeff(temperature, pressure);
450

451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
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
532
533
534
535
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
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
        else
            DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components "
                                                 << compIIdx << " and " << compJIdx
                                                 << " in phase " << phaseIdx);
        Valgrind::CheckDefined(result);
        return result;
    }

    /*!
     * \brief Given a phase's composition, temperature and pressure,
     *        return its specific enthalpy \f$\mathrm{[J/kg]}\f$.
     * \param fluidState The fluid state
     * \param phaseIdx The index of the phase
     *
     * See:
     * Class 2000
     * 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
     *
     * Formula (2.42):
     * the specific enthalpy of a gas phase result from the sum of (enthalpies*mass fraction) of the components
     */
    using Base::enthalpy;
    template <class FluidState>
    static Scalar enthalpy(const FluidState &fluidState,
                           int phaseIdx)
    {
        assert(phaseIdx == gPhaseIdx);

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

        Scalar XN2 = fluidState.massFraction(gPhaseIdx, N2Idx);
        Scalar XH2O = fluidState.massFraction(gPhaseIdx, H2OIdx);

        Scalar result = 0;
        result += XH2O * H2O::gasEnthalpy(T, p);
        result += XN2 * N2::gasEnthalpy(T, p);
        Valgrind::CheckDefined(result);

        return result;
    }

    /*!
    * \brief Returns the specific enthalpy \f$\mathrm{[J/kg]}\f$ of a component in a specific phase
    * \param fluidState The fluid state
    * \param phaseIdx The index of the phase
    * \param componentIdx The index of the component
    */
    template <class FluidState>
    static Scalar componentEnthalpy(const FluidState &fluidState,
                                    int phaseIdx,
                                    int componentIdx)
    {
        Scalar T = fluidState.temperature(nPhaseIdx);
        Scalar p = fluidState.pressure(nPhaseIdx);
        Valgrind::CheckDefined(T);
        Valgrind::CheckDefined(p);

        if (componentIdx ==  H2OIdx)
        {
            return H2O::gasEnthalpy(T, p);
        }
        else if (componentIdx == N2Idx)
        {
            return N2::gasEnthalpy(T, p);
        }
    }

   //for the boundary condition T = 573.15 K;
   template <class FluidState>
    static Scalar componentEnthalpyBorder(const FluidState &fluidState,
                                    int phaseIdx,
                                    int componentIdx)
    {
        Scalar T = 573.15;
        Scalar p = fluidState.pressure(nPhaseIdx);
        Valgrind::CheckDefined(T);
        Valgrind::CheckDefined(p);

        if (componentIdx ==  H2OIdx)
        {
            return H2O::gasEnthalpy(T, p);
        }
        else if (componentIdx == N2Idx)
        {
            return N2::gasEnthalpy(T, p);
        }
    }

    /*!
     * \brief Thermal conductivity of a fluid phase \f$\mathrm{[W/(m K)]}\f$.
     * \param phaseIdx The index of the fluid phase to consider
     *
     * \note For the thermal conductivity of the phases the contribution of the minor
     *       component is neglected. This contribution is probably not big, but somebody
     *       would have to find out its influence.
     */
    using Base::thermalConductivity;
    template <class FluidState>
    static Scalar thermalConductivity(const FluidState &fluidState,
                                      int phaseIdx)
    {
        assert(0 <= phaseIdx  && phaseIdx < numPhases);
        Scalar temperature  = fluidState.temperature(phaseIdx) ;
        Scalar pressure = fluidState.pressure(phaseIdx);

        // Isobaric Properties for Air and Carbondioxide in: NIST Standard
        // Reference Database Number 69, Eds. P.J. Linstrom and
        // W.G. Mallard evaluated at p=.1 MPa, does not
        // change dramatically with p
        // and can be interpolated linearly with temperature
//         Scalar lambdaPureN2 = 6.525e-5 * temperature + 0.024031;
        Scalar lambdaPureN2 = N2::gasThermalConductivity(temperature, pressure);

        if (useComplexRelations){
            Scalar xN2 = fluidState.moleFraction(phaseIdx, N2Idx);
            Scalar xH2O = fluidState.moleFraction(phaseIdx, H2OIdx);
            Scalar lambdaN2 = xN2 * lambdaPureN2;
                // Assuming Raoult's, Daltons law and ideal gas
                // in order to obtain the partial density of water in the air phase
            if(xH2O <= 0+ 1e-6) return lambdaN2;

                Scalar partialPressure  = pressure * xH2O;
                Scalar lambdaH2O = xH2O * H2O::gasThermalConductivity(temperature, partialPressure);

                return lambdaN2 + lambdaH2O;
            }
            else
                return lambdaPureN2; // conductivity of Air [W / (m K ) ]
    }

    /*!
     * \brief Specific isobaric heat capacity of a fluid phase.
     *        \f$\mathrm{[J/(kg*K)}\f$.
     * \param fluidState An abitrary fluid state
     * \param phaseIdx The index of the fluid phase to consider
     *
     * \note The calculation of the isobaric heat capacity is preliminary. A better
     *       description of the influence of the composition on the phase property
     *       has to be found.
     */
    using Base::heatCapacity;
    template <class FluidState>
    static Scalar heatCapacity(const FluidState &fluidState,
                               int phaseIdx)
    {
        const Scalar temperature  = fluidState.temperature(phaseIdx);
        const Scalar pressure = fluidState.pressure(phaseIdx);

        Scalar c_pN2;
        Scalar c_pH2O;
        // let the water and air components do things their own way
            c_pN2= N2::gasHeatCapacity(fluidState.temperature(phaseIdx),
                                        fluidState.pressure(phaseIdx)
                                        * fluidState.moleFraction(phaseIdx, N2Idx));

            c_pH2O = H2O::gasHeatCapacity(fluidState.temperature(phaseIdx),
                                          fluidState.pressure(phaseIdx)
                                          * fluidState.moleFraction(phaseIdx, H2OIdx));

        return
            c_pH2O*fluidState.moleFraction(nPhaseIdx, H2OIdx) + c_pN2*fluidState.moleFraction(nPhaseIdx, N2Idx);
    }

};

} // end namespace
} // end namespace

#endif