volumevariables.hh 41.5 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
    using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
104
    using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
105
public:
Timo Koch's avatar
Timo Koch committed
106
    //! The type of the object returned by the fluidState() method
107
108
109
    using FluidState = typename Traits::FluidState;
    //! The fluid system used here
    using FluidSystem = typename Traits::FluidSystem;
110
111
    //! Export the indices
    using Indices = typename ModelTraits::Indices;
112
    //! Export type of solid state
113
    using SolidState = typename Traits::SolidState;
114
    //! Export type of solid system
115
    using SolidSystem = typename Traits::SolidSystem;
116
    //! Export the primary variable switch
117
    using PrimaryVariableSwitch = TwoPNCPrimaryVariableSwitch;
118

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

124
    // check for permissive combinations
125
126
    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!");
127
128
    static_assert((formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0), "Chosen TwoPFormulation not supported!");

129
    /*!
130
     * \brief Updates all quantities for a given control volume.
131
132
133
134
135
136
137
     *
     * \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
    */
138
139
    template<class ElemSol, class Problem, class Element, class Scv>
    void update(const ElemSol& elemSol, const Problem& problem, const Element& element, const Scv& scv)
140
    {
141
        ParentType::update(elemSol, problem, element, scv);
142
        asImp_().completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
143

144
145
        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
        const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
146

147
148
        const int wPhaseIdx = fluidState_.wettingPhase();
        const int nPhaseIdx = 1 - wPhaseIdx;
149
150

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

154
        // porosity & permeabilty
Katharina Heck's avatar
Katharina Heck committed
155
        updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
156
        EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
157
        permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
158

159
160
161
162
163
164
        auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
        {
            return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
        };

        effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
165
166

        EnergyVolVars::updateEffectiveThermalConductivity();
167
168
169
    }

    /*!
170
     * \brief Sets complete fluid state.
171
172
     *
     * \param elemSol A vector containing all primary variables connected to the element
173
174
175
176
177
178
     * \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
179
180
     *
     * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies.
181
     */
182
    template<class ElemSol, class Problem, class Element, class Scv>
183
184
185
186
    void completeFluidState(const ElemSol& elemSol,
                            const Problem& problem,
                            const Element& element,
                            const Scv& scv,
187
188
                            FluidState& fluidState,
                            SolidState& solidState)
189
    {
190
        EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState);
191

192
        const auto& priVars = elemSol[scv.localDofIndex()];
193
        const auto phasePresence = priVars.state();
194

195
196
        using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
197
198
        const auto wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
        fluidState.setWettingPhase(wPhaseIdx);
199

200
        // set the saturations
201
202
203
204
205
206
207
208
209
210
        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);
        }
211
212
        else if (phasePresence == bothPhases)
        {
213
214
215
216
217
            if (formulation == TwoPFormulation::p0s1)
            {
                fluidState.setSaturation(phase1Idx, priVars[switchIdx]);
                fluidState.setSaturation(phase0Idx, 1 - priVars[switchIdx]);
            }
218
            else
219
220
221
222
            {
                fluidState.setSaturation(phase0Idx, priVars[switchIdx]);
                fluidState.setSaturation(phase1Idx, 1 - priVars[switchIdx]);
            }
223
        }
224
225
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase presence.");
226

227
        // set pressures of the fluid phases
228
        pc_ = MaterialLaw::pc(materialParams, fluidState.saturation(wPhaseIdx));
229
230
231
        if (formulation == TwoPFormulation::p0s1)
        {
            fluidState.setPressure(phase0Idx, priVars[pressureIdx]);
232
            fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_
233
                                                                       : priVars[pressureIdx] - pc_);
234
        }
235
236
237
        else
        {
            fluidState.setPressure(phase1Idx, priVars[pressureIdx]);
238
            fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_
239
                                                                       : priVars[pressureIdx] + pc_);
240
        }
241
242
243
244
245
246
247
248
249
250
251
252
253
254
   }

    /*!
     * \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_; }

255
    /*!
Kilian Weishaupt's avatar
Kilian Weishaupt committed
256
     * \brief Returns the average molar mass \f$\mathrm{[kg/mol]}\f$ of the fluid phase.
257
258
259
260
261
262
     *
     * \param phaseIdx The phase index
     */
    Scalar averageMolarMass(int phaseIdx) const
    { return fluidState_.averageMolarMass(phaseIdx); }

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
    /*!
     * \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_; }

375
376
377
378
    /*!
     * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
     */
    Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
379
380
381
382
383
    {
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updatePhase(fluidState_, phaseIdx);
        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
    }
384
385
386
387

    /*!
     * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
     */
388
    Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
389
390
391
392
393
    { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }

    /*!
     * \brief Returns the wetting phase index
     */
394
    int wettingPhase() const
395
    { return fluidState_.wettingPhase(); }
396
397
398
399
400
401
402
403
404
405
406

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_;

407
408
    // Effective diffusion coefficients for the phases
    DiffusionCoefficients effectiveDiffCoeff_;
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
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509

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();
510
511
512
513

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

514
515
516
        // If constraint solver is not used, get the phase pressures and set the fugacity coefficients here
        if(!useConstraintSolver)
        {
517
            for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++ phaseIdx)
518
            {
519
                assert(FluidSystem::isIdealMixture(phaseIdx));
520
                for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++ compIdx) {
521
522
523
524
525
526
527
                    Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
                    fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
                }
            }
        }

        // now comes the tricky part: calculate phase compositions
528
529
        const Scalar p0 = fluidState.pressure(phase0Idx);
        const Scalar p1 = fluidState.pressure(phase1Idx);
530
531
        if (phasePresence == bothPhases)
        {
532
533
534
            // both phases are present, phase compositions are a result
            // of the equilibrium between the phases. This is the job
            // of the "MiscibleMultiPhaseComposition" constraint solver
535
            if(useConstraintSolver)
536
                MiscibleMultiPhaseComposition::solve(fluidState,
537
                                                     paramCache);
538
            // ... or calculated explicitly this way ...
539
540
            else
            {
541
542
543
544
545
546
547
                // 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;
548

549
                // calculate the mole fractions of the components within the nonwetting phase
550
551
                const Scalar xnn = partPressGas / p1;
                const Scalar xnw = partPressLiquid / p1;
552
553

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

559
                // set all mole fractions
560
561
562
563
                fluidState.setMoleFraction(phase0Idx, comp0Idx, xww);
                fluidState.setMoleFraction(phase0Idx, comp1Idx, xwn);
                fluidState.setMoleFraction(phase1Idx, comp0Idx, xnw);
                fluidState.setMoleFraction(phase1Idx, comp1Idx, xnn);
564
565
            }
        }
566
        else if (phasePresence == secondPhaseOnly)
567
        {
568
            // only the second phase is present, composition is stored explicitly.
569
            if( useMoles() )
570
            {
571
572
                fluidState.setMoleFraction(phase1Idx, comp1Idx, 1 - priVars[switchIdx]);
                fluidState.setMoleFraction(phase1Idx, comp0Idx, priVars[switchIdx]);
573
            }
574
            // setMassFraction() has only to be called 1-numComponents times
575
            else
576
                fluidState.setMassFraction(phase1Idx, comp0Idx, priVars[switchIdx]);
577
578
579
580

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). This is the job
            // of the "ComputeFromReferencePhase" constraint solver
581
            if (useConstraintSolver)
582
583
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
584
                                                 phase1Idx);
585
            // ... or calculated explicitly this way ...
586
587
            else
            {
588
589
                // note that the water phase is actually not existing!
                // thus, this is used as phase switch criterion
590
591
                const Scalar xnw = priVars[switchIdx];
                const Scalar xnn = 1.0 - xnw;
592

593
                // first, xww:
594
595
596
597
598
599
                // 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)
600
                const Scalar xww = xnw*p1/( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx)*p0 );
601
602
603
604
605

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

609
610
                fluidState.setMoleFraction(phase0Idx, comp0Idx, xww);
                fluidState.setMoleFraction(phase0Idx, comp1Idx, xwn);
611
612
            }
        }
613
        else if (phasePresence == firstPhaseOnly)
614
        {
615
616
            // only the wetting phase is present, i.e. wetting phase
            // composition is stored explicitly.
617
            if( useMoles() ) // mole-fraction formulation
618
            {
619
620
                fluidState.setMoleFraction(phase0Idx, comp0Idx, 1-priVars[switchIdx]);
                fluidState.setMoleFraction(phase0Idx, comp1Idx, priVars[switchIdx]);
621
            }
622
            // setMassFraction() has only to be called 1-numComponents times
623
            else // mass-fraction formulation
624
                fluidState.setMassFraction(phase0Idx, comp1Idx, priVars[switchIdx]);
625
626
627
628

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). This is the job
            // of the "ComputeFromReferencePhase" constraint solver
629
            if (useConstraintSolver)
630
631
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
632
                                                 phase0Idx);
633
            // ... or calculated explicitly this way ...
634
635
            else
            {
636
637
                // note that the gas phase is actually not existing!
                // thus, this is used as phase switch criterion
638
                const Scalar xwn = priVars[switchIdx];
639

640
641
642
                // first, xnw:
                // psteam = xnw * pn = partial pressure of water in gas phase
                // psteam = fugacityCoefficient * pw
643
                const Scalar xnw = ( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx)*p0 )/p1;
644

645
                // second, xnn:
646
647
648
649
650
                // xwn = partialPressure / Henry
                // partialPressure = pn * xnn
                // xwn = pn * xnn / Henry
                // xnn = xwn * Henry / pn
                // Henry = fugacityCoefficient * pw
651
                const Scalar xnn = xwn*( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx)*p0 )/p1;
652

653
654
                fluidState.setMoleFraction(phase1Idx, comp1Idx, xnn);
                fluidState.setMoleFraction(phase1Idx, comp0Idx, xnw);
655
656
657
            }
        }

658
        for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
659
        {
Timo Koch's avatar
Timo Koch committed
660
            // set the viscosity and desity here if constraintsolver is not used
661
662
            if(!useConstraintSolver)
            {
Timo Koch's avatar
Timo Koch committed
663
664
665
                paramCache.updateComposition(fluidState, phaseIdx);
                const Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
                fluidState.setDensity(phaseIdx, rho);
666
667
                Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
                fluidState.setMolarDensity(phaseIdx, rhoMolar);
668
            }
669

670
            // compute and set the enthalpy
671
672
            const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
            fluidState.setViscosity(phaseIdx,mu);
673
            Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx);
674
675
676
677
            fluidState.setEnthalpy(phaseIdx, h);
        }
   }

678
679
680
681
682
683
684
685
686
687
688
689
690
691
};
/*!
 * \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> >;
692

693
694
    using Scalar = typename Traits::PrimaryVariables::value_type;
    using ModelTraits = typename Traits::ModelTraits;
695

696
697
698
699
700
701
702
703
704
    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
    };
705

706
707
708
709
710
711
712
    // phase presence indices
    enum
    {
        firstPhaseOnly = ModelTraits::Indices::firstPhaseOnly,
        secondPhaseOnly = ModelTraits::Indices::secondPhaseOnly,
        bothPhases = ModelTraits::Indices::bothPhases
    };
713

714
715
716
717
718
719
    // primary variable indices
    enum
    {
        switchIdx = ModelTraits::Indices::switchIdx,
        pressureIdx = ModelTraits::Indices::pressureIdx
    };
720

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

724
725
726
    using PermeabilityType = typename Traits::PermeabilityType;
    using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, typename Traits::FluidSystem>;
    using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition< Scalar, typename Traits::FluidSystem >;
727

728
729
730
731
732
733
734
735
736
737
738
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;
739

740
741
742
743
    //! 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; }
744

745
746
747
748
749
750
751
752
    // 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");
753
754

    /*!
755
     * \brief Sets complete fluid state.
756
     *
757
758
759
760
761
762
763
764
765
     * \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.
766
     */
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
    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);
        }
   }
799
800

    /*!
801
     * \brief Updates composition of all phases from the primary variables.
802
     *
803
804
805
     *        \param actualFluidState Container for all the secondary variables concerning the fluids
     *        \param paramCache Container for cache parameters
     *        \param priVars The primary Variables
806
     */
807
808
809
810
811
   void updateMoleFraction(FluidState &actualFluidState,
                           typename Traits::FluidSystem::ParameterCache & paramCache,
                           const typename Traits::PrimaryVariables& priVars)
   {
        const auto phasePresence = priVars.state();
812

813
814
815
816
        Scalar xwnNonEquil = 0.0;
        Scalar xwwNonEquil = 0.0;
        Scalar xnwNonEquil = 0.0;
        Scalar xnnNonEquil = 0.0;
817

818
819
820
821
822
        if (phasePresence == bothPhases)
        {
            xwnNonEquil = priVars[ModelTraits::numFluidPhases()];
            xwwNonEquil = 1-xwnNonEquil;
            xnwNonEquil = priVars[ModelTraits::numFluidPhases()+comp1Idx];
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
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
            //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);
        }
   }
923

924
    /*!
925
926
927
928
929
     * \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
930
     */
931
    const Scalar xEquil(const unsigned int phaseIdx, const unsigned int compIdx) const
932
    {
933
        return xEquil_[phaseIdx][compIdx] ;
934
    }
935

936
private:
937
    std::array<std::array<Scalar, numFluidComps>, ModelTraits::numFluidPhases()> xEquil_;
938
939
};

940
} // end namespace Dumux
941
942

#endif