volumevariables.hh 40.8 KB
Newer Older
1
2
3
4
5
6
7
// -*- 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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
13
14
15
16
17
18
19
20
 *   (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 TwoPTwoCModel
22
23
24
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
 */
25

26
27
28
#ifndef DUMUX_2P2C_VOLUME_VARIABLES_HH
#define DUMUX_2P2C_VOLUME_VARIABLES_HH

29
#include <dumux/material/fluidstates/compositional.hh>
30
31
32
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>

33
#include <dumux/discretization/method.hh>
34
#include <dumux/porousmediumflow/volumevariables.hh>
35
#include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
36
#include <dumux/porousmediumflow/2p/formulation.hh>
Timo Koch's avatar
Timo Koch committed
37
#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
38
#include <dumux/porousmediumflow/2pnc/primaryvariableswitch.hh>
Timo Koch's avatar
Timo Koch committed
39
40

namespace Dumux {
41

42
43
44
45
46
// forward declaration
template <class Traits,  bool enableChemicalNonEquilibrium, bool useConstraintSolver>
class TwoPTwoCVolumeVariablesImplementation;


47
48
49
50
51
/*!
 * \ingroup TwoPTwoCModel
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
 */
52
template <class Traits, bool useConstraintSolver = true>
53
54
55
56
57
58
59
60
61
62
using TwoPTwoCVolumeVariables =  TwoPTwoCVolumeVariablesImplementation<Traits,  Traits::ModelTraits::enableChemicalNonEquilibrium(), useConstraintSolver>;

/*!
 * \ingroup TwoPTwoCModel
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
 *        This is the base class for a 2p2c model with and without chemical nonequilibrium
 */
template <class Traits, class Impl>
class TwoPTwoCVolumeVariablesBase
63
: public PorousMediumFlowVolumeVariables<Traits>
Katharina Heck's avatar
Katharina Heck committed
64
, public EnergyVolumeVariables<Traits, TwoPTwoCVolumeVariables<Traits> >
65
{
66
    using ParentType = PorousMediumFlowVolumeVariables< Traits>;
67
    using EnergyVolVars = EnergyVolumeVariables<Traits, TwoPTwoCVolumeVariables<Traits> >;
68

69
70
    using Scalar = typename Traits::PrimaryVariables::value_type;
    using ModelTraits = typename Traits::ModelTraits;
71

72
    static constexpr int numFluidComps = ParentType::numFluidComponents();
73
    // component indices
74
75
    enum
    {
76
77
78
79
        comp0Idx = Traits::FluidSystem::comp0Idx,
        comp1Idx = Traits::FluidSystem::comp1Idx,
        phase0Idx = Traits::FluidSystem::phase0Idx,
        phase1Idx = Traits::FluidSystem::phase1Idx
80
81
    };

82
    // phase presence indices
83
84
    enum
    {
85
86
        firstPhaseOnly = ModelTraits::Indices::firstPhaseOnly,
        secondPhaseOnly = ModelTraits::Indices::secondPhaseOnly,
87
        bothPhases = ModelTraits::Indices::bothPhases
88
89
90
    };

    // primary variable indices
91
92
    enum
    {
93
94
        switchIdx = ModelTraits::Indices::switchIdx,
        pressureIdx = ModelTraits::Indices::pressureIdx
95
96
    };

97
98
99
100
101
    // formulations
    static constexpr auto formulation = ModelTraits::priVarFormulation();

    using PermeabilityType = typename Traits::PermeabilityType;
    using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, typename Traits::FluidSystem>;
102
    using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition< Scalar, typename Traits::FluidSystem >;
103
public:
Timo Koch's avatar
Timo Koch committed
104
    //! The type of the object returned by the fluidState() method
105
106
107
    using FluidState = typename Traits::FluidState;
    //! The fluid system used here
    using FluidSystem = typename Traits::FluidSystem;
108
    //! Export type of solid state
109
    using SolidState = typename Traits::SolidState;
110
    //! Export type of solid system
111
    using SolidSystem = typename Traits::SolidSystem;
112
    //! Export the primary variable switch
113
    using PrimaryVariableSwitch = TwoPNCPrimaryVariableSwitch;
114

115
    //! Return whether moles or masses are balanced
116
    static constexpr bool useMoles() { return ModelTraits::useMoles(); }
117
    //! Return the two-phase formulation used here
118
    static constexpr TwoPFormulation priVarFormulation() { return formulation; }
119

120
    // check for permissive combinations
121
122
    static_assert(ModelTraits::numFluidPhases() == 2, "NumPhases set in the model is not two!");
    static_assert(ModelTraits::numFluidComponents() == 2, "NumComponents set in the model is not two!");
123
124
    static_assert((formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0), "Chosen TwoPFormulation not supported!");

125
    /*!
126
     * \brief Updates all quantities for a given control volume.
127
128
129
130
131
132
133
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param problem The object specifying the problem which ought to
     *                be simulated
     * \param element An element which contains part of the control volume
     * \param scv The sub control volume
    */
134
135
    template<class ElemSol, class Problem, class Element, class Scv>
    void update(const ElemSol& elemSol, const Problem& problem, const Element& element, const Scv& scv)
136
    {
137
        ParentType::update(elemSol, problem, element, scv);
138
        asImp_().completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
139

140
141
        // Second instance of a parameter cache. Could be avoided if
        // diffusion coefficients also became part of the fluid state.
142
143
144
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState_);

145
146
        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
        const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
147
148
149
150
151
152
153
154
155
156
157

        const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
        const int nPhaseIdx = 1 - wPhaseIdx;

        // relative permeabilities -> require wetting phase saturation as parameter!
        relativePermeability_[wPhaseIdx] = MaterialLaw::krw(matParams, saturation(wPhaseIdx));
        relativePermeability_[nPhaseIdx] = MaterialLaw::krn(matParams, saturation(wPhaseIdx));

        // binary diffusion coefficients
        diffCoeff_[phase0Idx] = FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phase0Idx, comp0Idx, comp1Idx);
        diffCoeff_[phase1Idx] = FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phase1Idx, comp0Idx, comp1Idx);
158

159
        // porosity & permeabilty
Katharina Heck's avatar
Katharina Heck committed
160
        updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
161
        EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
162
        permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
163
164
165
    }

    /*!
166
     * \brief Sets complete fluid state.
167
168
     *
     * \param elemSol A vector containing all primary variables connected to the element
169
170
171
172
173
174
     * \param problem The object specifying the problem which ought to
     *                be simulated
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     * \param fluidState A container with the current (physical) state of the fluid
     * \param solidState A container with the current (physical) state of the solid
175
176
     *
     * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies.
177
     */
178
    template<class ElemSol, class Problem, class Element, class Scv>
179
180
181
182
    void completeFluidState(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
183
184
                            FluidState& fluidState,
                            SolidState& solidState)
185
    {
186
        EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState);
187

188
        const auto& priVars = elemSol[scv.localDofIndex()];
189
        const auto phasePresence = priVars.state();
190

191
192
193
194
195
        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
        const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
        fluidState.setWettingPhase(wPhaseIdx);

196
        // set the saturations
197
198
199
200
201
202
203
204
205
206
        if (phasePresence == firstPhaseOnly)
        {
            fluidState.setSaturation(phase0Idx, 1.0);
            fluidState.setSaturation(phase1Idx, 0.0);
        }
        else if (phasePresence == secondPhaseOnly)
        {
            fluidState.setSaturation(phase0Idx, 0.0);
            fluidState.setSaturation(phase1Idx, 1.0);
        }
207
208
        else if (phasePresence == bothPhases)
        {
209
210
211
212
213
            if (formulation == TwoPFormulation::p0s1)
            {
                fluidState.setSaturation(phase1Idx, priVars[switchIdx]);
                fluidState.setSaturation(phase0Idx, 1 - priVars[switchIdx]);
            }
214
            else
215
216
217
218
            {
                fluidState.setSaturation(phase0Idx, priVars[switchIdx]);
                fluidState.setSaturation(phase1Idx, 1 - priVars[switchIdx]);
            }
219
        }
220
221
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase presence.");
222

223
224
225
226
227
228
229
        // set pressures of the fluid phases
        pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
        if (formulation == TwoPFormulation::p0s1)
        {
            fluidState.setPressure(phase0Idx, priVars[pressureIdx]);
            fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_
                                                                       : priVars[pressureIdx] - pc_);
230
        }
231
232
233
234
235
        else
        {
            fluidState.setPressure(phase1Idx, priVars[pressureIdx]);
            fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_
                                                                       : priVars[pressureIdx] + pc_);
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
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
   }

    /*!
     * \brief Returns the phase state within the control volume.
     */
    const FluidState &fluidState() const
    { return fluidState_; }

    /*!
     * \brief Returns the phase state for the control-volume.
     */
    const SolidState &solidState() const
    { return solidState_; }

    /*!
     * \brief Returns the saturation of a given phase within
     *        the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar saturation(const int phaseIdx) const
    { return fluidState_.saturation(phaseIdx); }

    /*!
     * \brief Returns the mass fraction of a given component in a
     *        given phase within the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     * \param compIdx The component index
     */
    Scalar massFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.massFraction(phaseIdx, compIdx); }

    /*!
     * \brief Returns the mole fraction of a given component in a
     *        given phase within the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     * \param compIdx The component index
     */
    Scalar moleFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.moleFraction(phaseIdx, compIdx); }

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume in \f$[kg/m^3]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar density(const int phaseIdx) const
    { return fluidState_.density(phaseIdx); }

    /*!
     * \brief Returns the dynamic viscosity of the fluid within the
     *        control volume in \f$\mathrm{[Pa s]}\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar viscosity(const int phaseIdx) const
    { return fluidState_.viscosity(phaseIdx); }

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume in \f$[mol/m^3]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar molarDensity(const int phaseIdx) const
    { return fluidState_.molarDensity(phaseIdx); }

    /*!
     * \brief Returns the effective pressure of a given phase within
     *        the control volume in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar pressure(const int phaseIdx) const
    { return fluidState_.pressure(phaseIdx); }

    /*!
     * \brief Returns temperature within the control volume in \f$[K]\f$.
     *
     * Note that we assume thermodynamic equilibrium, i.e. the
     * temperature of the rock matrix and of all fluid phases are
     * identical.
     */
    Scalar temperature() const
    { return fluidState_.temperature(/*phaseIdx=*/0); }

    /*!
     * \brief Returns the relative permeability of a given phase within
     *        the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar relativePermeability(const int phaseIdx) const
    { return relativePermeability_[phaseIdx]; }

    /*!
     * \brief Returns the effective mobility of a given phase within
     *        the control volume in \f$[s*m/kg]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar mobility(const int phaseIdx) const
    { return relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx); }

    /*!
     * \brief Returns the effective capillary pressure within the control volume
     *        in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
     */
    Scalar capillaryPressure() const
    { return pc_; }

    /*!
     * \brief Returns the average porosity within the control volume in \f$[-]\f$.
     */
    Scalar porosity() const
    { return solidState_.porosity(); }

    /*!
     * \brief Returns the average permeability within the control volume in \f$[m^2]\f$.
     */
    const PermeabilityType& permeability() const
    { return permeability_; }

    /*!
     * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
     */
    Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
    {
        if(phaseIdx == compIdx)
            DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
        else
            return diffCoeff_[phaseIdx];
    }

private:
    FluidState fluidState_;
    SolidState solidState_;

    Scalar pc_;                     // The capillary pressure
    PermeabilityType permeability_; // Effective permeability within the control volume

    // Relative permeability within the control volume
    std::array<Scalar, ModelTraits::numFluidPhases()> relativePermeability_;

    // Binary diffusion coefficients for the phases
    std::array<Scalar, ModelTraits::numFluidPhases()> diffCoeff_;

protected:
    const Impl &asImp_() const { return *static_cast<const Impl*>(this); }
    Impl &asImp_() { return *static_cast<Impl*>(this); }
};
/*!
 * \ingroup TwoPTwoCModel
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
 *        Specialization for chemical equilibrium
 */
template <class Traits, bool useConstraintSolver>
class TwoPTwoCVolumeVariablesImplementation<Traits, false, useConstraintSolver>
: public TwoPTwoCVolumeVariablesBase<Traits, TwoPTwoCVolumeVariablesImplementation<Traits, false, useConstraintSolver>>
{
    using ParentType = TwoPTwoCVolumeVariablesBase< Traits, TwoPTwoCVolumeVariablesImplementation<Traits, false, useConstraintSolver>>;
    using EnergyVolVars = EnergyVolumeVariables<Traits, TwoPTwoCVolumeVariables<Traits> >;

    using Scalar = typename Traits::PrimaryVariables::value_type;
    using ModelTraits = typename Traits::ModelTraits;

    static constexpr int numFluidComps = ParentType::numFluidComponents();
    // component indices
    enum
    {
        comp0Idx = Traits::FluidSystem::comp0Idx,
        comp1Idx = Traits::FluidSystem::comp1Idx,
        phase0Idx = Traits::FluidSystem::phase0Idx,
        phase1Idx = Traits::FluidSystem::phase1Idx
    };

    // phase presence indices
    enum
    {
        firstPhaseOnly = ModelTraits::Indices::firstPhaseOnly,
        secondPhaseOnly = ModelTraits::Indices::secondPhaseOnly,
        bothPhases = ModelTraits::Indices::bothPhases
    };

    // primary variable indices
    enum
    {
        switchIdx = ModelTraits::Indices::switchIdx,
        pressureIdx = ModelTraits::Indices::pressureIdx
    };

    // formulations
    static constexpr auto formulation = ModelTraits::priVarFormulation();

    using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, typename Traits::FluidSystem>;
    using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition< Scalar, typename Traits::FluidSystem >;
public:
    //! The type of the object returned by the fluidState() method
    using FluidState = typename Traits::FluidState;
    //! The fluid system used here
    using FluidSystem = typename Traits::FluidSystem;
    //! Export type of solid state
    using SolidState = typename Traits::SolidState;
    //! Export type of solid system
    using SolidSystem = typename Traits::SolidSystem;
    //! Export the primary variable switch
    using PrimaryVariableSwitch = TwoPNCPrimaryVariableSwitch;

    //! Return whether moles or masses are balanced
    static constexpr bool useMoles() { return ModelTraits::useMoles(); }
    //! Return the two-phase formulation used here
    static constexpr TwoPFormulation priVarFormulation() { return formulation; }

    // check for permissive combinations
    static_assert(useMoles() || (!useMoles() && useConstraintSolver), "if !UseMoles, UseConstraintSolver has to be set to true");

    // The computations in the explicit composition update most probably assume a liquid-gas interface with
    // liquid as first phase. TODO: is this really needed? The constraint solver does the job anyway, doesn't it?
    static_assert(useConstraintSolver || (!FluidSystem::isGas(phase0Idx) && FluidSystem::isGas(phase1Idx)),
                   "Explicit composition calculation has to be re-checked for NON-liquid-gas equilibria");

    /*!
     * \brief Sets complete fluid state.
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param problem The object specifying the problem which ought to
     *                be simulated
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     * \param fluidState A container with the current (physical) state of the fluid
     * \param solidState A container with the current (physical) state of the solid
     *
     * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies.
     */
    template<class ElemSol, class Problem, class Element, class Scv>
    void completeFluidState(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
                            FluidState& fluidState,
                            SolidState& solidState)
    {
        ParentType::completeFluidState(elemSol, problem, element, scv, fluidState, solidState);

        const auto& priVars = elemSol[scv.localDofIndex()];
        const auto phasePresence = priVars.state();
487
488
489
490

        // calculate the phase compositions
        typename FluidSystem::ParameterCache paramCache;

491
492
493
        // If constraint solver is not used, get the phase pressures and set the fugacity coefficients here
        if(!useConstraintSolver)
        {
494
            for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++ phaseIdx)
495
            {
496
                assert(FluidSystem::isIdealMixture(phaseIdx));
497
                for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++ compIdx) {
498
499
500
501
502
503
504
                    Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
                    fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
                }
            }
        }

        // now comes the tricky part: calculate phase compositions
505
506
        const Scalar p0 = fluidState.pressure(phase0Idx);
        const Scalar p1 = fluidState.pressure(phase1Idx);
507
508
        if (phasePresence == bothPhases)
        {
509
510
511
            // both phases are present, phase compositions are a result
            // of the equilibrium between the phases. This is the job
            // of the "MiscibleMultiPhaseComposition" constraint solver
512
            if(useConstraintSolver)
513
                MiscibleMultiPhaseComposition::solve(fluidState,
514
                                                     paramCache);
515
            // ... or calculated explicitly this way ...
516
517
            else
            {
518
519
520
521
522
523
524
                // get the partial pressure of the main component of the first phase within the
                // second phase == vapor pressure due to equilibrium. Note that in this case the
                // fugacityCoefficient * p is the vapor pressure (see implementation in respective fluidsystem)
                const Scalar partPressLiquid = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx)*p0;

                // get the partial pressure of the main component of the gas phase
                const Scalar partPressGas = p1 - partPressLiquid;
525

526
                // calculate the mole fractions of the components within the nonwetting phase
527
528
                const Scalar xnn = partPressGas / p1;
                const Scalar xnw = partPressLiquid / p1;
529
530

                // calculate the mole fractions of the components within the wetting phase
531
                // note that in this case the fugacityCoefficient * p is the Henry Coefficient
532
                // (see implementation in respective fluidsystem)
533
534
                const Scalar xwn = partPressGas / (FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx)*p0);
                const Scalar xww = 1.0 - xwn;
535

536
                // set all mole fractions
537
538
539
540
                fluidState.setMoleFraction(phase0Idx, comp0Idx, xww);
                fluidState.setMoleFraction(phase0Idx, comp1Idx, xwn);
                fluidState.setMoleFraction(phase1Idx, comp0Idx, xnw);
                fluidState.setMoleFraction(phase1Idx, comp1Idx, xnn);
541
542
            }
        }
543
        else if (phasePresence == secondPhaseOnly)
544
        {
545
            // only the second phase is present, composition is stored explicitly.
546
            if( useMoles() )
547
            {
548
549
                fluidState.setMoleFraction(phase1Idx, comp1Idx, 1 - priVars[switchIdx]);
                fluidState.setMoleFraction(phase1Idx, comp0Idx, priVars[switchIdx]);
550
            }
551
            // setMassFraction() has only to be called 1-numComponents times
552
            else
553
                fluidState.setMassFraction(phase1Idx, comp0Idx, priVars[switchIdx]);
554
555
556
557

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). This is the job
            // of the "ComputeFromReferencePhase" constraint solver
558
            if (useConstraintSolver)
559
560
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
561
                                                 phase1Idx);
562
            // ... or calculated explicitly this way ...
563
564
            else
            {
565
566
                // note that the water phase is actually not existing!
                // thus, this is used as phase switch criterion
567
568
                const Scalar xnw = priVars[switchIdx];
                const Scalar xnn = 1.0 - xnw;
569

570
                // first, xww:
571
572
573
574
575
576
                // xnw * pn = "actual" (hypothetical) vapor pressure
                // fugacityCoefficient * pw = vapor pressure given by thermodynamic conditions
                // Here, xww is not actually the mole fraction of water in the wetting phase
                // xww is only the ratio of "actual" vapor pressure / "thermodynamic" vapor pressure
                // If xww > 1 : gas is over-saturated with water vapor,
                // condensation takes place (see switch criterion in model)
577
                const Scalar xww = xnw*p1/( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx)*p0 );
578
579
580
581
582

                // second, xwn:
                // partialPressure / xwn = Henry
                // partialPressure = xnn * pn
                // xwn = xnn * pn / Henry
583
                // Henry = fugacityCoefficient * pw
584
                const Scalar xwn = xnn*p1/( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx)*p0 );
585

586
587
                fluidState.setMoleFraction(phase0Idx, comp0Idx, xww);
                fluidState.setMoleFraction(phase0Idx, comp1Idx, xwn);
588
589
            }
        }
590
        else if (phasePresence == firstPhaseOnly)
591
        {
592
593
            // only the wetting phase is present, i.e. wetting phase
            // composition is stored explicitly.
594
            if( useMoles() ) // mole-fraction formulation
595
            {
596
597
                fluidState.setMoleFraction(phase0Idx, comp0Idx, 1-priVars[switchIdx]);
                fluidState.setMoleFraction(phase0Idx, comp1Idx, priVars[switchIdx]);
598
            }
599
            // setMassFraction() has only to be called 1-numComponents times
600
            else // mass-fraction formulation
601
                fluidState.setMassFraction(phase0Idx, comp1Idx, priVars[switchIdx]);
602
603
604
605

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). This is the job
            // of the "ComputeFromReferencePhase" constraint solver
606
            if (useConstraintSolver)
607
608
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
609
                                                 phase0Idx);
610
            // ... or calculated explicitly this way ...
611
612
            else
            {
613
614
                // note that the gas phase is actually not existing!
                // thus, this is used as phase switch criterion
615
                const Scalar xwn = priVars[switchIdx];
616

617
618
619
                // first, xnw:
                // psteam = xnw * pn = partial pressure of water in gas phase
                // psteam = fugacityCoefficient * pw
620
                const Scalar xnw = ( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx)*p0 )/p1;
621

622
                // second, xnn:
623
624
625
626
627
                // xwn = partialPressure / Henry
                // partialPressure = pn * xnn
                // xwn = pn * xnn / Henry
                // xnn = xwn * Henry / pn
                // Henry = fugacityCoefficient * pw
628
                const Scalar xnn = xwn*( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx)*p0 )/p1;
629

630
631
                fluidState.setMoleFraction(phase1Idx, comp1Idx, xnn);
                fluidState.setMoleFraction(phase1Idx, comp0Idx, xnw);
632
633
634
            }
        }

635
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
636
        {
Timo Koch's avatar
Timo Koch committed
637
            // set the viscosity and desity here if constraintsolver is not used
638
639
            if(!useConstraintSolver)
            {
Timo Koch's avatar
Timo Koch committed
640
641
642
                paramCache.updateComposition(fluidState, phaseIdx);
                const Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
                fluidState.setDensity(phaseIdx, rho);
643
644
                Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
                fluidState.setMolarDensity(phaseIdx, rhoMolar);
645
            }
646

647
            // compute and set the enthalpy
648
649
            const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
            fluidState.setViscosity(phaseIdx,mu);
650
            Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx);
651
652
653
654
            fluidState.setEnthalpy(phaseIdx, h);
        }
   }

655
656
657
658
659
660
661
662
663
664
665
666
667
668
};
/*!
 * \ingroup TwoPTwoCModel
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
 *        Specialization for chemical non-equilibrium.
 *        The equilibrium mole fraction is calculated using Henry's and Raoult's law
 */
template <class Traits, bool useConstraintSolver>
class TwoPTwoCVolumeVariablesImplementation<Traits, true, useConstraintSolver>
: public TwoPTwoCVolumeVariablesBase<Traits, TwoPTwoCVolumeVariablesImplementation<Traits, true, useConstraintSolver>>
{
    using ParentType = TwoPTwoCVolumeVariablesBase< Traits, TwoPTwoCVolumeVariablesImplementation<Traits, true, useConstraintSolver>>;
    using EnergyVolVars = EnergyVolumeVariables<Traits, TwoPTwoCVolumeVariables<Traits> >;
669

670
671
    using Scalar = typename Traits::PrimaryVariables::value_type;
    using ModelTraits = typename Traits::ModelTraits;
672

673
674
675
676
677
678
679
680
681
    static constexpr int numFluidComps = ParentType::numFluidComponents();
    // component indices
    enum
    {
        comp0Idx = Traits::FluidSystem::comp0Idx,
        comp1Idx = Traits::FluidSystem::comp1Idx,
        phase0Idx = Traits::FluidSystem::phase0Idx,
        phase1Idx = Traits::FluidSystem::phase1Idx
    };
682

683
684
685
686
687
688
689
    // phase presence indices
    enum
    {
        firstPhaseOnly = ModelTraits::Indices::firstPhaseOnly,
        secondPhaseOnly = ModelTraits::Indices::secondPhaseOnly,
        bothPhases = ModelTraits::Indices::bothPhases
    };
690

691
692
693
694
695
696
    // primary variable indices
    enum
    {
        switchIdx = ModelTraits::Indices::switchIdx,
        pressureIdx = ModelTraits::Indices::pressureIdx
    };
697

698
699
    // formulations
    static constexpr auto formulation = ModelTraits::priVarFormulation();
700

701
702
703
    using PermeabilityType = typename Traits::PermeabilityType;
    using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, typename Traits::FluidSystem>;
    using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition< Scalar, typename Traits::FluidSystem >;
704

705
706
707
708
709
710
711
712
713
714
715
public:
    //! The type of the object returned by the fluidState() method
    using FluidState = typename Traits::FluidState;
    //! The fluid system used here
    using FluidSystem = typename Traits::FluidSystem;
    //! Export type of solid state
    using SolidState = typename Traits::SolidState;
    //! Export type of solid system
    using SolidSystem = typename Traits::SolidSystem;
    //! Export the primary variable switch
    using PrimaryVariableSwitch = TwoPNCPrimaryVariableSwitch;
716

717
718
719
720
    //! Return whether moles or masses are balanced
    static constexpr bool useMoles() { return ModelTraits::useMoles(); }
    //! Return the two-phase formulation used here
    static constexpr TwoPFormulation priVarFormulation() { return formulation; }
721

722
723
724
725
726
727
728
729
    // check for permissive combinations
    static_assert(useMoles() || (!useMoles() && useConstraintSolver), "if !UseMoles, UseConstraintSolver has to be set to true");
    static_assert((formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0), "Chosen TwoPFormulation not supported!");

    // The computations in the explicit composition update most probably assume a liquid-gas interface with
    // liquid as first phase. TODO: is this really needed? The constraint solver does the job anyway, doesn't it?
    static_assert(useConstraintSolver || (!FluidSystem::isGas(phase0Idx) && FluidSystem::isGas(phase1Idx)),
                   "Explicit composition calculation has to be re-checked for NON-liquid-gas equilibria");
730
731

    /*!
732
     * \brief Sets complete fluid state.
733
     *
734
735
736
737
738
739
740
741
742
     * \param elemSol A vector containing all primary variables connected to the element
     * \param problem The object specifying the problem which ought to
     *                be simulated
     * \param element An element which contains part of the control volume
     * \param scv The sub-control volume
     * \param fluidState A container with the current (physical) state of the fluid
     * \param solidState A container with the current (physical) state of the solid
     *
     * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies.
743
     */
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
    template<class ElemSol, class Problem, class Element, class Scv>
    void completeFluidState(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
                            FluidState& fluidState,
                            SolidState& solidState)
    {
        ParentType::completeFluidState(elemSol, problem, element, scv, fluidState, solidState);

        const auto& priVars = elemSol[scv.localDofIndex()];

        /////////////
        // set the fluid compositions
        /////////////
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState);

        updateMoleFraction(fluidState,
                           paramCache,
                           priVars);


        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
        {
            // compute and set the enthalpy
            const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
            fluidState.setViscosity(phaseIdx,mu);
            Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx);
            fluidState.setEnthalpy(phaseIdx, h);
        }
   }
776
777

    /*!
778
     * \brief Updates composition of all phases from the primary variables.
779
     *
780
781
782
     *        \param actualFluidState Container for all the secondary variables concerning the fluids
     *        \param paramCache Container for cache parameters
     *        \param priVars The primary Variables
783
     */
784
785
786
787
788
   void updateMoleFraction(FluidState &actualFluidState,
                           typename Traits::FluidSystem::ParameterCache & paramCache,
                           const typename Traits::PrimaryVariables& priVars)
   {
        const auto phasePresence = priVars.state();
789

790
791
792
793
        Scalar xwnNonEquil = 0.0;
        Scalar xwwNonEquil = 0.0;
        Scalar xnwNonEquil = 0.0;
        Scalar xnnNonEquil = 0.0;
794

795
796
797
798
799
        if (phasePresence == bothPhases)
        {
            xwnNonEquil = priVars[ModelTraits::numFluidPhases()];
            xwwNonEquil = 1-xwnNonEquil;
            xnwNonEquil = priVars[ModelTraits::numFluidPhases()+comp1Idx];
800

801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
            //if the capillary pressure is high, we need to use Kelvin's equation to reduce the vapor pressure
            if (actualFluidState.saturation(phase0Idx) < 0.01)
            {
                const Scalar p1 = actualFluidState.pressure(phase1Idx);
                const Scalar partPressLiquid = FluidSystem::vaporPressure(actualFluidState, comp0Idx);
                xnwNonEquil =std::min(partPressLiquid/p1, xnwNonEquil);
            }
            xnnNonEquil = 1- xnwNonEquil;

            actualFluidState.setMoleFraction(phase0Idx, comp0Idx, xwwNonEquil);
            actualFluidState.setMoleFraction(phase0Idx, comp1Idx, xwnNonEquil);
            actualFluidState.setMoleFraction(phase1Idx, comp0Idx, xnwNonEquil);
            actualFluidState.setMoleFraction(phase1Idx, comp1Idx, xnnNonEquil);

            //check if we need this
            for(int phaseIdx=0; phaseIdx<ModelTraits::numFluidPhases(); ++phaseIdx)
                for (int compIdx = 0; compIdx < numFluidComps; ++compIdx)
                {
                    const Scalar phi = FluidSystem::fugacityCoefficient(actualFluidState,
                                                                        paramCache,
                                                                        phaseIdx,
                                                                        compIdx);
                    actualFluidState.setFugacityCoefficient(phaseIdx,
                                                            compIdx,
                                                            phi);
                }

            FluidState equilFluidState; // the fluidState *on the interface* i.e. chemical equilibrium
            equilFluidState.assign(actualFluidState) ;
            // If constraint solver is not used, get the phase pressures and set the fugacity coefficients here
            if(!useConstraintSolver)
            {
                for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++ phaseIdx)
                {
                    assert(FluidSystem::isIdealMixture(phaseIdx));
                    for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++ compIdx) {
                        Scalar phi = FluidSystem::fugacityCoefficient(equilFluidState, paramCache, phaseIdx, compIdx);
                        equilFluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
                    }
                }
            }

            // now comes the tricky part: calculate phase compositions
            const Scalar p0 = equilFluidState.pressure(phase0Idx);
            const Scalar p1 = equilFluidState.pressure(phase1Idx);
            // both phases are present, phase compositions are a result
            // of the equilibrium between the phases. This is the job
            // of the "MiscibleMultiPhaseComposition" constraint solver
            if(useConstraintSolver)
                MiscibleMultiPhaseComposition::solve(equilFluidState,
                                                     paramCache);
            // ... or calculated explicitly this way ...
            else
            {
                // get the partial pressure of the main component of the first phase within the
                // second phase == vapor pressure due to equilibrium. Note that in this case the
                // fugacityCoefficient * p is the vapor pressure (see implementation in respective fluidsystem)
                const Scalar partPressLiquid = FluidSystem::fugacityCoefficient(equilFluidState, phase0Idx, comp0Idx)*p0;

                // get the partial pressure of the main component of the gas phase
                const Scalar partPressGas = p1 - partPressLiquid;

                // calculate the mole fractions of the components within the nonwetting phase
                const Scalar xnn = partPressGas / p1;
                const Scalar xnw = partPressLiquid / p1;

                // calculate the mole fractions of the components within the wetting phase
                // note that in this case the fugacityCoefficient * p is the Henry Coefficient
                // (see implementation in respective fluidsystem)
                const Scalar xwn = partPressGas / (FluidSystem::fugacityCoefficient(equilFluidState, phase0Idx, comp1Idx)*p0);
                const Scalar xww = 1.0 - xwn;

                // set all mole fractions
                equilFluidState.setMoleFraction(phase0Idx, comp0Idx, xww);
                equilFluidState.setMoleFraction(phase0Idx, comp1Idx, xwn);
                equilFluidState.setMoleFraction(phase1Idx, comp0Idx, xnw);
                equilFluidState.setMoleFraction(phase1Idx, comp1Idx, xnn);
            }

            // Setting the equilibrium composition (in a kinetic model not necessarily the same as the actual mole fraction)
            for(int phaseIdx=0; phaseIdx<ModelTraits::numFluidPhases(); ++phaseIdx){
                for (int compIdx=0; compIdx< numFluidComps; ++ compIdx){
                    xEquil_[phaseIdx][compIdx] = equilFluidState.moleFraction(phaseIdx, compIdx);
                }
            }
        }
        else
        {
            DUNE_THROW(Dune::InvalidStateException, "nonequilibrium is only possible for 2 phases present ");
        }
        // compute densities of all phases
        for(int phaseIdx=0; phaseIdx<ModelTraits::numFluidPhases(); ++phaseIdx)
        {
            const Scalar rho = FluidSystem::density(actualFluidState, paramCache, phaseIdx);
            actualFluidState.setDensity(phaseIdx, rho);
            const Scalar rhoMolar = FluidSystem::molarDensity(actualFluidState, paramCache, phaseIdx);
            actualFluidState.setMolarDensity(phaseIdx, rhoMolar);
        }
   }
900

901
    /*!
902
903
904
905
906
     * \brief The mole fraction we would have in the case of chemical equilibrium /
     *        on the interface.
     *
     * \param phaseIdx The index of the fluid phase
     * \param compIdx The local index of the component
907
     */
908
    const Scalar xEquil(const unsigned int phaseIdx, const unsigned int compIdx) const
909
    {
910
        return xEquil_[phaseIdx][compIdx] ;
911
    }
912

913
private:
914
    std::array<std::array<Scalar, numFluidComps>, ModelTraits::numFluidPhases()> xEquil_;
915
916
};

917
} // end namespace Dumux
918
919

#endif