volumevariables.hh 20.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
21
 * \ingroup TwoPNCModel
22
23
24
25
26
27
28
29
30
31
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase, n-component model.
 */
#ifndef DUMUX_2PNC_VOLUME_VARIABLES_HH
#define DUMUX_2PNC_VOLUME_VARIABLES_HH

#include <iostream>
#include <vector>

#include <dumux/common/math.hh>
32
#include <dumux/common/properties.hh>
Timo Koch's avatar
Timo Koch committed
33
#include <dumux/discretization/methods.hh>
34

35
#include <dumux/material/fluidstates/compositional.hh>
36
#include <dumux/porousmediumflow/volumevariables.hh>
37
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
38
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
39

Timo Koch's avatar
Timo Koch committed
40
#include "indices.hh" // for formulation
41

42
namespace Dumux {
43
44
45
46
47
48
49

/*!
 * \ingroup TwoPNCModel
 * \brief Contains the quantities which are are constant within a
 *        finite volume in the two-phase, n-component model.
 */
template <class TypeTag>
50
class TwoPNCVolumeVariables : public PorousMediumFlowVolumeVariables<TypeTag>
51
{
52
    using ParentType = PorousMediumFlowVolumeVariables<TypeTag>;
53
54
55
56
    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
57
    using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
58
    using PermeabilityType = typename SpatialParams::PermeabilityType;
59
    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
60
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
61
62
63
64
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
    using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);

65
66
    enum
    {
67
68
69
        numPhases = GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases(),
        numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents(),
        numMajorComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numMajorComponents(),
70
71
72

        // formulations
        formulation = GET_PROP_VALUE(TypeTag, Formulation),
73
74
        pwsn = TwoPNCFormulation::pwsn,
        pnsw = TwoPNCFormulation::pnsw,
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

        // phase indices
        wPhaseIdx = FluidSystem::wPhaseIdx,
        nPhaseIdx = FluidSystem::nPhaseIdx,

        // component indices
        wCompIdx = FluidSystem::wCompIdx,
        nCompIdx = FluidSystem::nCompIdx,

        // phase presence enums
        nPhaseOnly = Indices::nPhaseOnly,
        wPhaseOnly = Indices::wPhaseOnly,
        bothPhases = Indices::bothPhases,

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

94
    using Element = typename GridView::template Codim<0>::Entity;
95
    using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
96
    using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, FluidSystem>;
97
98
    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
    static_assert(useMoles, "use moles has to be set true in the 2pnc model");
99

100
101
public:

102
    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
103
104

    /*!
105
106
107
108
109
110
111
112
     * \brief Update all quantities for a given control volume
     *
     * \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
    */
113
114
    template<class ElementSolution>
    void update(const ElementSolution &elemSol,
115
116
                const Problem &problem,
                const Element &element,
117
                const SubControlVolume& scv)
118
    {
119
        ParentType::update(elemSol, problem, element, scv);
120

121
        Implementation::completeFluidState(elemSol, problem, element, scv, fluidState_);
122
123
124
125

        /////////////
        // calculate the remaining quantities
        /////////////
126
        const auto &materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
127

128
        // Second instance of a parameter cache.
129
130
131
132
133
        // Could be avoided if diffusion coefficients also
        // became part of the fluid state.
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState_);
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
134
135
        {
            // relative permeabilities
136
137
138
            Scalar kr;
            if (phaseIdx == wPhaseIdx)
                kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
139
            else // ATTENTION: krn requires the wetting-phase saturation as parameter!
140
141
142
                kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));

            mobility_[phaseIdx] = kr / fluidState_.viscosity(phaseIdx);
143

144
145
146
            int compIIdx = phaseIdx;
            for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
            {
Thomas Fetzer's avatar
Thomas Fetzer committed
147
                // binary diffusion coefficients
148
149
                if(compIIdx!= compJIdx)
                {
150
151
152
153
154
155
                    setDiffusionCoefficient_(phaseIdx, compJIdx,
                                             FluidSystem::binaryDiffusionCoefficient(fluidState_,
                                                                                     paramCache,
                                                                                     phaseIdx,
                                                                                     compIIdx,
                                                                                     compJIdx));
156
157
                }
            }
158
        }
159

160
161
162
        // porosity & permeability
        porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
        permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
163
164
    }

165
    /*!
166
167
168
169
170
171
172
173
174
     * \brief Complete the fluid state
     *
     * \param elemSol A vector containing all primary variables connected to the element
     * \param problem The problem
     * \param element The element
     * \param scv The sub control volume
     * \param fluidState The fluid state
     *
     * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies.
175
     */
176
177
    template<class ElementSolution>
    static void completeFluidState(const ElementSolution& elemSol,
178
179
                                   const Problem& problem,
                                   const Element& element,
180
181
                                   const SubControlVolume& scv,
                                   FluidState& fluidState)
182
    {
183
        Scalar t = ParentType::temperature(elemSol, problem, element, scv);
184
185
        fluidState.setTemperature(t);

186
        auto&& priVars = ParentType::extractDofPriVars(elemSol, scv);
187
        const auto phasePresence = priVars.state();
188
189
190
191
192

        /////////////
        // set the saturations
        /////////////

193
        Scalar Sn;
194
        if (phasePresence == nPhaseOnly)
195
        {
196
            Sn = 1.0;
197
198
199
        }
        else if (phasePresence == wPhaseOnly)
        {
200
            Sn = 0.0;
201
        }
202
203
        else if (phasePresence == bothPhases)
        {
204
            if (formulation == pwsn)
205
                Sn = priVars[switchIdx];
206
            else if (formulation == pnsw)
207
                Sn = 1.0 - priVars[switchIdx];
Timo Koch's avatar
Timo Koch committed
208
209
            else
                DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
210
        }
211
212
213
214
215
        else
        {
            DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
        }

216
217
        fluidState.setSaturation(nPhaseIdx, Sn);
        fluidState.setSaturation(wPhaseIdx, 1.0 - Sn);
218
219
220
221
222
223

        /////////////
        // set the pressures of the fluid phases
        /////////////

        // calculate capillary pressure
224
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
225
        auto pc = MaterialLaw::pc(materialParams, 1 - Sn);
226
227

        // extract the pressures
228
        if (formulation == pwsn)
229
        {
230
231
            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
            if (priVars[pressureIdx] + pc < 0.0)
232
                 DUNE_THROW(NumericalProblem,"Capillary pressure is too low");
233
            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx] + pc);
234
        }
235
        else if (formulation == pnsw)
236
        {
237
            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx]);
238
            // Here we check for (p_g - pc) in order to ensure that (p_l > 0)
239
            if (priVars[pressureIdx] - pc < 0.0)
240
            {
241
                std::cout<< "p_n: "<< priVars[pressureIdx]<<" Cap_press: "<< pc << std::endl;
242
                DUNE_THROW(NumericalProblem,"Capillary pressure is too high");
243
            }
244
            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx] - pc);
245
        }
246
247
248
249
        else
        {
            DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
        }
250
251
252
253
254

        /////////////
        // calculate the phase compositions
        /////////////

255
        typename FluidSystem::ParameterCache paramCache;
256
257

        // now comes the tricky part: calculate phase composition
258
259
        if (phasePresence == bothPhases)
        {
260
            // both phases are present, phase composition results from
261
            // the nonwetting <-> wetting equilibrium. This is
262
263
264
265
            // the job of the "MiscibleMultiPhaseComposition"
            // constraint solver

            // set the known mole fractions in the fluidState so that they
266
            // can be used by the MiscibleMultiPhaseComposition constraint solver
267
268
269
270
271
272
            unsigned int knownPhaseIdx = nPhaseIdx;
            if (GET_PROP_VALUE(TypeTag, SetMoleFractionsForWettingPhase))
            {
                knownPhaseIdx = wPhaseIdx;
            }

273
274
            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
            {
275
                fluidState.setMoleFraction(knownPhaseIdx, compIdx, priVars[compIdx]);
276
277
            }

278
            MiscibleMultiPhaseComposition::solve(fluidState,
Timo Koch's avatar
Timo Koch committed
279
280
                                           paramCache,
                                           /*setViscosity=*/true,
281
282
                                           /*setEnthalpy=*/false,
                                           knownPhaseIdx);
283
        }
284
285
        else if (phasePresence == nPhaseOnly)
        {
286
287
            Dune::FieldVector<Scalar, numComponents> moleFrac;

Martin Schneider's avatar
Martin Schneider committed
288
            moleFrac[wCompIdx] = priVars[switchIdx];
289
290

            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
291
                moleFrac[compIdx] = priVars[compIdx];
292

293
            Scalar sumMoleFracOtherComponents = 0;
294
            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
Martin Schneider's avatar
Martin Schneider committed
295
                sumMoleFracOtherComponents += moleFrac[compIdx];
296

297
298
            sumMoleFracOtherComponents += moleFrac[wCompIdx];
            moleFrac[nCompIdx] = 1 - sumMoleFracOtherComponents;
299
300
301

            // Set fluid state mole fractions
            for (int compIdx=0; compIdx<numComponents; ++compIdx)
302
                fluidState.setMoleFraction(nPhaseIdx, compIdx, moleFrac[compIdx]);
303
304
305

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). this is the job
306
307
308
            // of the "ComputeFromReferencePhase" constraint solver
            ComputeFromReferencePhase::solve(fluidState,
                                             paramCache,
Martin Schneider's avatar
Martin Schneider committed
309
                                             nPhaseIdx,
310
311
                                             /*setViscosity=*/true,
                                             /*setEnthalpy=*/false);
312
        }
Martin Schneider's avatar
Martin Schneider committed
313
314
        else if (phasePresence == wPhaseOnly)
        {
315
        // only the wetting phase is present, i.e. wetting phase
316
        // composition is stored explicitly.
317
        // extract _mass_ fractions in the non-wetting phase
318
319
            Dune::FieldVector<Scalar, numComponents> moleFrac;

320
321
            moleFrac[nCompIdx] = priVars[switchIdx];

322
            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
323
                moleFrac[compIdx] = priVars[compIdx];
324

325
            Scalar sumMoleFracOtherComponents = 0;
326
            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
327
                sumMoleFracOtherComponents += moleFrac[compIdx];
328

329
330
            sumMoleFracOtherComponents += moleFrac[nCompIdx];
            moleFrac[wCompIdx] = 1 - sumMoleFracOtherComponents;
331
332
333
334
335
336
337
338
339

            // convert mass to mole fractions and set the fluid state
            for (int compIdx=0; compIdx<numComponents; ++compIdx)
            {
                fluidState.setMoleFraction(wPhaseIdx, compIdx, moleFrac[compIdx]);
            }

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). this is the job
340
341
342
343
344
345
            // of the "ComputeFromReferencePhase" constraint solver
            ComputeFromReferencePhase::solve(fluidState,
                                             paramCache,
                                             wPhaseIdx,
                                             /*setViscosity=*/true,
                                             /*setEnthalpy=*/false);
346
347
348
349
350
351
        }
        paramCache.updateAll(fluidState);
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
            Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
            Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
Timo Koch's avatar
Timo Koch committed
352
            Scalar h = ParentType::enthalpy(fluidState, paramCache, phaseIdx);
353
354
355

            fluidState.setDensity(phaseIdx, rho);
            fluidState.setViscosity(phaseIdx, mu);
Timo Koch's avatar
Timo Koch committed
356
            fluidState.setEnthalpy(phaseIdx, h);
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
        }
    }

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

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

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume.
     *
     * \param phaseIdx The phase index
     */
    Scalar density(int phaseIdx) const
    {
        if (phaseIdx < numPhases)
            return fluidState_.density(phaseIdx);
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
Timo Koch's avatar
Timo Koch committed
387
388
389
390
391
392
393
394
395
396
397
398
399
400
    }

    /*!
     * \brief Returns the kinematic viscosity of a given phase within the
     *        control volume.
     *
     * \param phaseIdx The phase index
     */
    Scalar viscosity(int phaseIdx) const
    {
        if (phaseIdx < numPhases)
            return fluidState_.viscosity(phaseIdx);
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
    }

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume.
     *
     * \param phaseIdx The phase index
     */
    Scalar molarDensity(int phaseIdx) const
    {
        if (phaseIdx < numPhases)
            return fluidState_.molarDensity(phaseIdx);
        else
            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
    }

    /*!
     * \brief Returns the effective pressure of a given phase within
     *        the control volume.
     *
     * \param phaseIdx The phase index
     */
    Scalar pressure(int phaseIdx) const
424
    { return fluidState_.pressure(phaseIdx); }
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442

    /*!
     * \brief Returns temperature inside the sub-control volume.
     *
     * 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 effective mobility of a given phase within
     *        the control volume.
     *
     * \param phaseIdx The phase index
     */
    Scalar mobility(int phaseIdx) const
443
    { return mobility_[phaseIdx]; }
444
445
446
447
448
449
450
451
452
453
454
455
456
457

    /*!
     * \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 fluidState_.pressure(FluidSystem::nPhaseIdx) - fluidState_.pressure(FluidSystem::wPhaseIdx); }

    /*!
     * \brief Returns the average porosity within the control volume.
     */
    Scalar porosity() const
    { return porosity_; }

458
459
460
    /*!
     * \brief Returns the permeability within the control volume.
     */
461
    const PermeabilityType& permeability() const
462
463
    { return permeability_; }

464
465

    /*!
Martin Schneider's avatar
Martin Schneider committed
466
     * \brief Returns the diffusion coefficient
467
     */
468
469
470
471
472
473
474
    Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
    {
        if (compIdx < phaseIdx)
            return diffCoefficient_[phaseIdx][compIdx];
        else if (compIdx > phaseIdx)
            return diffCoefficient_[phaseIdx][compIdx-1];
        else
Thomas Fetzer's avatar
Thomas Fetzer committed
475
            DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
476
    }
477
478
479
480
481
482
483
484

    /*!
     * \brief Returns the molarity of a component in the phase
     *
     * \param phaseIdx the index of the fluid phase
     * \param compIdx the index of the component
     */
     Scalar molarity(int phaseIdx, int compIdx) const // [moles/m^3]
485
    { return fluidState_.molarity(phaseIdx, compIdx);}
486
487
488
489
490
491
492
493

     /*!
      * \brief Returns the mass fraction of a component in the phase
      *
      * \param phaseIdx the index of the fluid phase
      * \param compIdx the index of the component
      */
     Scalar massFraction(int phaseIdx, int compIdx) const
494
     { return fluidState_.massFraction(phaseIdx, compIdx); }
495
496
497
498
499
500
501
502

     /*!
      * \brief Returns the mole fraction of a component in the phase
      *
      * \param phaseIdx the index of the fluid phase
      * \param compIdx the index of the component
      */
     Scalar moleFraction(int phaseIdx, int compIdx) const
503
     { return fluidState_.moleFraction(phaseIdx, compIdx); }
504
505
506
507
508

protected:
    FluidState fluidState_;

private:
509
510
511
512
513
514
515
    void setDiffusionCoefficient_(int phaseIdx, int compIdx, Scalar d)
    {
        if (compIdx < phaseIdx)
            diffCoefficient_[phaseIdx][compIdx] = std::move(d);
        else if (compIdx > phaseIdx)
            diffCoefficient_[phaseIdx][compIdx-1] = std::move(d);
        else
Thomas Fetzer's avatar
Thomas Fetzer committed
516
            DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient for phaseIdx = compIdx doesn't exist");
517
518
    }

519
520
521
    Scalar porosity_;               //!< Effective porosity within the control volume
    PermeabilityType permeability_; //!> Effective permeability within the control volume
    Scalar mobility_[numPhases];    //!< Effective mobility within the control volume
522
    std::array<std::array<Scalar, numComponents-1>, numPhases> diffCoefficient_;
523
524
};

525
} // end namespace Dumux
526
527

#endif