stokesncmodel.hh 8.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// -*- 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
 *
Thomas Fetzer's avatar
[docu]    
Thomas Fetzer committed
22
 * \brief Adaption of the box scheme to the n-component Stokes model.
23
 */
24
25
#ifndef DUMUX_STOKESNC_MODEL_HH
#define DUMUX_STOKESNC_MODEL_HH
26

27
28
#include <dune/common/version.hh>

29
30
#include <dumux/freeflow/stokes/stokesmodel.hh>

31
32
#include "stokesnclocalresidual.hh"
#include "stokesncproperties.hh"
33
34
35

namespace Dumux {
/*!
36
 * \ingroup BoxStokesncModel
Thomas Fetzer's avatar
[docu]    
Thomas Fetzer committed
37
 * \brief Adaption of the box scheme to the compositional Stokes model.
38
 *
39
 * This model implements an isothermal n-component Stokes flow of a fluid
Thomas Fetzer's avatar
Thomas Fetzer committed
40
41
42
 * solving a momentum balance, a mass balance and conservation equations for \f$n-1\f$
 * components. When using mole fractions naturally the densities represent molar
 * densities
43
 *
Thomas Fetzer's avatar
Thomas Fetzer committed
44
 * The momentum balance:
45
 * \f[
Thomas Fetzer's avatar
Thomas Fetzer committed
46
47
48
49
50
 *    \frac{\partial \left(\varrho_g {\boldsymbol{v}}_g\right)}{\partial t}
 *    + \text{div} \left( p_g {\bf {I}}
 *    - \mu_g \left( \textbf{grad}\, \boldsymbol{v}_g
 *                   + \textbf{grad}\, \boldsymbol{v}_g^T \right) \right)
 *    - \varrho_g {\bf g} = 0
51
 * \f]
Thomas Fetzer's avatar
Thomas Fetzer committed
52
53
54
55
56
57
 * By setting the property <code>EnableNavierStokes</code> to <code>true</code> the Navier-Stokes
 * equation can be solved. In this case an additional term
 * \f[
 *    + \text{div} \left( \varrho_g \boldsymbol{v}_g \boldsymbol{v}_g \right)
 * \f]
 * is added to the momentum balance equation.
58
 *
Thomas Fetzer's avatar
Thomas Fetzer committed
59
 * The mass balance equation:
60
 * \f[
Thomas Fetzer's avatar
Thomas Fetzer committed
61
62
 *    \frac{\partial \varrho_g}{\partial t}
 *    + \text{div} \left(\varrho_g {\boldsymbol{v}}_g\right) - q_g = 0
63
64
 * \f]
 *
Thomas Fetzer's avatar
Thomas Fetzer committed
65
 * The component mass balance equations:
66
 * \f[
Thomas Fetzer's avatar
Thomas Fetzer committed
67
68
69
70
 *    \frac{\partial \left(\varrho_g X_g^\kappa\right)}{\partial t}
 *    + \text{div} \left( \varrho_g {\boldsymbol{v}}_g X_g^\kappa
 *    - D^\kappa_g \varrho_g \frac{M^\kappa}{M_g} \textbf{grad}\, x_g^\kappa \right)
 *    - q_g^\kappa = 0
71
 * \f]
Thomas Fetzer's avatar
Thomas Fetzer committed
72
73
 * Please note that, even though it is n-component model, the diffusive
 * fluxes are still calculated with binary diffusion.
74
 *
Thomas Fetzer's avatar
Thomas Fetzer committed
75
76
 * This is discretized by a fully-coupled vertex-centered finite volume
 * (box) scheme in space and by the implicit Euler method in time.
77
78
 */
template<class TypeTag>
79
class StokesncModel : public StokesModel<TypeTag>
80
81
82
83
84
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;

85
86
87
88
    enum {  dim = GridView::dimension,
			transportCompIdx = Indices::transportCompIdx,
			phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx),
			useMoles = GET_PROP_VALUE(TypeTag, UseMoles),
89
			numComponents = Indices::numComponents
90
	};
91
92
93
94
95
96
97
98
99

    typedef typename GridView::template Codim<0>::Iterator ElementIterator;

    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;

    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
100

101
102
103
104
105
106
107

public:
    //! \copydoc ImplicitModel::addOutputVtkFields
    template <class MultiWriter>
    void addOutputVtkFields(const SolutionVector &sol,
                            MultiWriter &writer)
    {
108

109
		typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
110
111
112
113
114
115
        typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > VelocityField;

        const Scalar scale_ = GET_PROP_VALUE(TypeTag, Scaling);

        // create the required scalar fields
        unsigned numVertices = this->gridView_().size(dim);
116
        ScalarField &pN = *writer.allocateManagedBuffer(numVertices);
117
        ScalarField &delP = *writer.allocateManagedBuffer(numVertices);
118
119
120
121
122
		ScalarField &T = *writer.allocateManagedBuffer(numVertices);

		ScalarField *moleFraction[numComponents];
            for (int i = 0; i < numComponents; ++i)
                moleFraction[i] = writer.template allocateManagedBuffer<Scalar, 1>(numVertices);
123

124
125
126
		ScalarField *massFraction[numComponents];
		for (int i = 0; i < numComponents; ++i)
			massFraction[i] = writer.template allocateManagedBuffer<Scalar, 1>(numVertices);
127

128
		ScalarField &rho = *writer.allocateManagedBuffer(numVertices);
129
130
131
132
133
134
135
136
137
138
        ScalarField &mu = *writer.allocateManagedBuffer(numVertices);
        VelocityField &velocity = *writer.template allocateManagedBuffer<Scalar, dim> (numVertices);

        unsigned numElements = this->gridView_().size(0);
        ScalarField &rank = *writer.allocateManagedBuffer(numElements);

        FVElementGeometry fvGeometry;
        VolumeVariables volVars;
        ElementBoundaryTypes elemBcTypes;

Hao Wu's avatar
Hao Wu committed
139
140
141
        ElementIterator eIt = this->gridView_().template begin<0>();
        ElementIterator eEndIt = this->gridView_().template end<0>();
        for (; eIt != eEndIt; ++eIt)
142
        {
143
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
Hao Wu's avatar
Hao Wu committed
144
            int idx = this->elementMapper().index(*eIt);
145
#else
Hao Wu's avatar
Hao Wu committed
146
            int idx = this->elementMapper().map(*eIt);
147
#endif
148
            rank[idx] = this->gridView_().comm().rank();
149

Hao Wu's avatar
Hao Wu committed
150
151
            fvGeometry.update(this->gridView_(), *eIt);
            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
152

153
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
Hao Wu's avatar
Hao Wu committed
154
            int numLocalVerts = eIt->subEntities(dim);
155
#else
Hao Wu's avatar
Hao Wu committed
156
            int numLocalVerts = eIt->template count<dim>();
157
#endif
158
159
            for (int i = 0; i < numLocalVerts; ++i)
            {
160
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
Hao Wu's avatar
Hao Wu committed
161
                int vIdxGlobal = this->vertexMapper().subIndex(*eIt, i, dim);
162
#else
Hao Wu's avatar
Hao Wu committed
163
                int vIdxGlobal = this->vertexMapper().map(*eIt, i, dim);
164
#endif
165
                volVars.update(sol[vIdxGlobal],
166
                               this->problem_(),
Hao Wu's avatar
Hao Wu committed
167
                               *eIt,
168
169
170
171
                               fvGeometry,
                               i,
                               false);

172
173
                pN[vIdxGlobal] = volVars.pressure()*scale_;
                delP[vIdxGlobal] = volVars.pressure()*scale_ - 1e5;
174
175
				for (int compIdx = 0; compIdx < numComponents; ++compIdx)
                    {
176
177
178
179
                        (*moleFraction[compIdx])[vIdxGlobal]= volVars.moleFraction(compIdx);
						(*massFraction[compIdx])[vIdxGlobal]= volVars.massFraction(compIdx);
                        Valgrind::CheckDefined((*moleFraction[compIdx])[vIdxGlobal]);
						Valgrind::CheckDefined((*massFraction[compIdx])[vIdxGlobal]);
180
					}
181

182
				T   [vIdxGlobal] = volVars.temperature();
183

184
185
186
187
				rho[vIdxGlobal] = volVars.density()*scale_*scale_*scale_;
                mu[vIdxGlobal] = volVars.dynamicViscosity()*scale_;
                velocity[vIdxGlobal] = volVars.velocity();
                velocity[vIdxGlobal] *= 1/scale_;
188
189
            }
        }
190
191
		writer.attachVertexData(T, "temperature");
        writer.attachVertexData(pN, "P");
192
        writer.attachVertexData(delP, "delP");
Thomas Fetzer's avatar
Thomas Fetzer committed
193
194
195
196
197
198
199
200
201
202
203
204
205

        for (int j = 0; j < numComponents; ++j)
        {
            std::ostringstream moleFrac, massFrac;
            moleFrac << "x_" << FluidSystem::phaseName(phaseIdx)
                     << "^" << FluidSystem::componentName(j);
            writer.attachVertexData(*moleFraction[j], moleFrac.str().c_str());

            massFrac << "X_" << FluidSystem::phaseName(phaseIdx)
                     << "^" << FluidSystem::componentName(j);
            writer.attachVertexData(*massFraction[j], massFrac.str().c_str());
        }

206
207
208
209
210
211
212
213
        writer.attachVertexData(rho, "rho");
        writer.attachVertexData(mu, "mu");
        writer.attachVertexData(velocity, "v", dim);
    }
};

}

214
#include "stokesncpropertydefaults.hh"
215
216

#endif