mpncmodel.hh 8.41 KB
Newer Older
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
3
/*****************************************************************************
4
 *   Copyright (C) 2009-2012 by Andreas Lauser                               *
Andreas Lauser's avatar
Andreas Lauser committed
5
 *   Institute for Modelling Hydraulic and Environmental Systems             *
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
 *   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/>.   *
 *****************************************************************************/
#ifndef DUMUX_MPNC_MODEL_HH
#define DUMUX_MPNC_MODEL_HH

25
26
#include "mpncproperties.hh"
#include "mpncvtkwriter.hh"
27
28

#include <dumux/boxmodels/common/boxmodel.hh>
29
#include <tr1/array>
30
31
32
33
34

namespace Dumux
{
/*!
 * \ingroup MPNCModel
35
36
 * \brief A fully implicit model for M-phase, N-component flow using
 *        vertex centered finite volumes.
37
 *
38
39
40
41
 * This model implements a \f$M\f$-phase flow of a fluid mixture
 * composed of \f$N\f$ chemical species. The phases are denoted by
 * lower index \f$\alpha \in \{ 1, \dots, M \}\f$. All fluid phases
 * are mixtures of \f$N \geq M - 1\f$ chemical species which are
42
 * denoted by the upper index \f$\kappa \in \{ 1, \dots, N \} \f$.
43
 *
44
 * The standard multi-phase Darcy law is used as the equation for
45
46
47
 * the conservation of momentum:
 * \f[
     v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \boldsymbol{K}
Andreas Lauser's avatar
Andreas Lauser committed
48
     \left(
49
       \text{grad}\left(p_\alpha - \varrho_{\alpha} g\right)
Andreas Lauser's avatar
Andreas Lauser committed
50
     \right)
Andreas Lauser's avatar
Andreas Lauser committed
51
     \f]
52
53
 *
 * By inserting this into the equations for the conservation of the
54
55
 * mass of each component, one gets one mass-continuity equation for
 * each component \f$\kappa\f$
Andreas Lauser's avatar
Andreas Lauser committed
56
 * \f[
57
58
 \sum_{\kappa} \left(
    \phi \frac{\partial \varrho_\alpha x_\alpha^\kappa S_\alpha}{\partial t}
59
    -
60
    \mathrm{div}\;
61
    \left\{
62
       \frac{\varrho_\alpha}{\overline M_\alpha} x_\alpha^\kappa
63
       \frac{k_{r\alpha}}{\mu_\alpha} \boldsymbol{K}
64
       \mathbf{grad}\left( p_\alpha - \varrho_{\alpha} g\right)
65
    \right\}
Andreas Lauser's avatar
Andreas Lauser committed
66
    \right)
67
    = q^\kappa
Andreas Lauser's avatar
Andreas Lauser committed
68
    \f]
69
70
71
 * with \f$\overline M_\alpha\f$ being the average molar mass of the
 * phase \f$\alpha\f$: \f[ \overline M_\alpha = \sum_\kappa M^\kappa
 * \; x_\alpha^\kappa \f]
72
 *
73
74
 * For the missing \f$M\f$ model assumptions, the model assumes that
 * if a fluid phase is not present, the sum of the mole fractions of
75
76
77
78
 * this fluid phase is smaller than \f$1\f$, i.e.  
 * \f[
 * \forall \alpha: S_\alpha = 0 \implies \sum_\kappa x_\alpha^\kappa \leq 1
 * \f]
79
 *
80
81
82
 * Also, if a fluid phase may be present at a given spatial location
 * its saturation must be positive:
 * \f[ \forall \alpha: \sum_\kappa x_\alpha^\kappa = 1 \implies S_\alpha \geq 0 \f]
83
 *
84
 * Since at any given spatial location, a phase is always either
85
 * present or not present, the one of the strict equalities on the
86
87
88
89
90
 * right hand side is always true, i.e.
 * \f[ \forall \alpha: S_\alpha \left( \sum_\kappa x_\alpha^\kappa - 1 \right) = 0 \f]
 * always holds.
 *
 * These three equations constitute a non-linear complementarity
91
92
 * problem, which can be solved using so-called non-linear
 * complementarity functions \f$\Phi(a, b)\f$ which have the property
93
94
95
96
97
 * \f[\Phi(a,b) = 0 \iff a \geq0 \land b \geq0  \land a \cdot b = 0 \f]
 *
 * Several non-linear complementarity functions have been suggested,
 * e.g. the Fischer-Burmeister function
 * \f[ \Phi(a,b) = a + b - \sqrt{a^2 + b^2} \;. \f]
98
 * This model uses
Andreas Lauser's avatar
Andreas Lauser committed
99
 * \f[ \Phi(a,b) = \min \{a,  b \}\;, \f]
100
101
102
103
104
105
106
 * because of its piecewise linearity.
 *
 * These equations are then discretized using a fully-implicit vertex
 * centered finite volume scheme (often known as 'box'-scheme) for
 * spatial discretization and the implicit Euler method as temporal
 * discretization.
 *
107
108
109
 * The model assumes local thermodynamic equilibrium and uses the
 * following primary variables:
 * - The component fugacities \f$f^1, \dots, f^{N}\f$
110
111
112
 * - The pressure of the first phase \f$p_1\f$
 * - The saturations of the first \f$M-1\f$ phases \f$S_1, \dots, S_{M-1}\f$
 * - Temperature \f$T\f$ if the energy equation is enabled
113
114
115
116
117
118
 */
template<class TypeTag>
class MPNCModel : public BoxModel<TypeTag>
{
    typedef BoxModel<TypeTag> ParentType;

119
120
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
121
122
123
    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;

124
125
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
126
127
128
129

    typedef Dumux::MPNCVtkWriter<TypeTag> MPNCVtkWriter;

    enum {
130
131
132
133
134
135
136
137
138
139
        enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy),
        enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion),
        enableKinetic = GET_PROP_VALUE(TypeTag, EnableKinetic),
        enableKineticEnergy = GET_PROP_VALUE(TypeTag, EnableKineticEnergy),
        enableSmoothUpwinding = GET_PROP_VALUE(TypeTag, EnableSmoothUpwinding),
        enablePartialReassemble = GET_PROP_VALUE(TypeTag, EnablePartialReassemble),
        enableJacobianRecycling = GET_PROP_VALUE(TypeTag, EnableJacobianRecycling),
        numDiffMethod = GET_PROP_VALUE(TypeTag, NumericDifferenceMethod),
        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
        numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
140
        numEq = GET_PROP_VALUE(TypeTag, NumEq)
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
    };

public:
    ~MPNCModel()
    { delete vtkWriter_; };

    void init(Problem &problem)
    {
        ParentType::init(problem);
        vtkWriter_ = new MPNCVtkWriter(problem);

        if (this->gridView_().comm().rank() == 0)
            std::cout
                << "Initializing M-phase N-component model: \n"
                << "    phases: " << numPhases << "\n"
                << "    components: " << numComponents << "\n"
                << "    equations: " << numEq << "\n"
                << "    kinetic mass transfer: " << enableKinetic<< "\n"
                << "    kinetic energy transfer: " << enableKineticEnergy<< "\n"
                << "    diffusion: " << enableDiffusion << "\n"
                << "    energy equation: " << enableEnergy << "\n"
                << "    smooth upwinding: " << enableSmoothUpwinding << "\n"
                << "    partial jacobian reassembly: " << enablePartialReassemble << "\n"
                << "    numeric differentiation method: " << numDiffMethod << " (-1: backward, 0: central, +1 forward)\n"
                << "    jacobian recycling: " << enableJacobianRecycling << "\n";
    }

    /*!
     * \brief Compute the total storage inside one phase of all
     *        conservation quantities.
     */
    void globalPhaseStorage(PrimaryVariables &dest, int phaseIdx)
    {
        dest = 0;

        ElementIterator elemIt = this->gridView_().template begin<0>();
        const ElementIterator elemEndIt = this->gridView_().template end<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
            this->localResidual().addPhaseStorage(dest, *elemIt, phaseIdx);
        };

182
183
        if (this->gridView_().comm().size() > 1)
            dest = this->gridView_().comm().sum(dest);
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
    }

    /*!
     * \brief Add the result of the current timestep to the VTK output.
     */
    template <class MultiWriter>
    void addOutputVtkFields(const SolutionVector &sol,
                            MultiWriter &writer)
    {
        vtkWriter_->addCurrentSolution(writer);
    }

    MPNCVtkWriter *vtkWriter_;
};

}

201
#include "mpncpropertydefaults.hh"
202
203

#endif