mpncmodel.hh 8.41 KB
Newer Older
 Andreas Lauser committed Dec 16, 2011 1 2 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- // vi: set et ts=4 sw=4 sts=4:  Andreas Lauser committed Sep 26, 2011 3 /*****************************************************************************  Andreas Lauser committed Jan 10, 2012 4  * Copyright (C) 2009-2012 by Andreas Lauser *  Andreas Lauser committed Jan 31, 2012 5  * Institute for Modelling Hydraulic and Environmental Systems *  Andreas Lauser committed Sep 26, 2011 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  * University of Stuttgart, Germany * * email: .@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 . * *****************************************************************************/ #ifndef DUMUX_MPNC_MODEL_HH #define DUMUX_MPNC_MODEL_HH  Andreas Lauser committed Feb 13, 2012 25 26 #include "mpncproperties.hh" #include "mpncvtkwriter.hh"  Andreas Lauser committed Sep 26, 2011 27 28  #include  Andreas Lauser committed Dec 16, 2011 29 #include  Andreas Lauser committed Sep 26, 2011 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.  Andreas Lauser committed Sep 26, 2011 37  *  Andreas Lauser committed Jan 10, 2012 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  Andreas Lauser committed Feb 22, 2012 42  * denoted by the upper index \f$\kappa \in \{ 1, \dots, N \} \f$.  Andreas Lauser committed Sep 26, 2011 43  *  Andreas Lauser committed Jan 10, 2012 44  * The standard multi-phase Darcy law is used as the equation for  Andreas Lauser committed Sep 26, 2011 45 46 47  * the conservation of momentum: * \f[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \boldsymbol{K}  Andreas Lauser committed Jan 20, 2012 48  \left(  Andreas Lauser committed Feb 22, 2012 49  \text{grad}\left(p_\alpha - \varrho_{\alpha} g\right)  Andreas Lauser committed Feb 03, 2012 50  \right)  Andreas Lauser committed Jan 20, 2012 51  \f]  Andreas Lauser committed Sep 26, 2011 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 committed Jan 20, 2012 56  * \f[  Andreas Lauser committed Jan 10, 2012 57 58  \sum_{\kappa} \left( \phi \frac{\partial \varrho_\alpha x_\alpha^\kappa S_\alpha}{\partial t}  Andreas Lauser committed Sep 26, 2011 59  -  Andreas Lauser committed Feb 22, 2012 60  \mathrm{div}\;  Andreas Lauser committed Sep 26, 2011 61  \left\{  Andreas Lauser committed Jan 10, 2012 62  \frac{\varrho_\alpha}{\overline M_\alpha} x_\alpha^\kappa  Andreas Lauser committed Sep 26, 2011 63  \frac{k_{r\alpha}}{\mu_\alpha} \boldsymbol{K}  Andreas Lauser committed Feb 22, 2012 64  \mathbf{grad}\left( p_\alpha - \varrho_{\alpha} g\right)  Andreas Lauser committed Sep 26, 2011 65  \right\}  Andreas Lauser committed Jan 20, 2012 66  \right)  Andreas Lauser committed Jan 10, 2012 67  = q^\kappa  Andreas Lauser committed Jan 20, 2012 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]  Andreas Lauser committed Sep 26, 2011 72  *  Andreas Lauser committed Jan 10, 2012 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]  Andreas Lauser committed Sep 26, 2011 79  *  Andreas Lauser committed Jan 10, 2012 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]  Andreas Lauser committed Sep 26, 2011 83  *  Andreas Lauser committed Jan 10, 2012 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  Andreas Lauser committed Jan 10, 2012 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  Andreas Lauser committed Jan 10, 2012 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 committed Jan 20, 2012 99  * \f[ \Phi(a,b) = \min \{a, b \}\;, \f]  Andreas Lauser committed Jan 10, 2012 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$  Andreas Lauser committed Jan 10, 2012 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  Andreas Lauser committed Sep 26, 2011 113 114 115 116 117 118  */ template class MPNCModel : public BoxModel { typedef BoxModel ParentType;  Andreas Lauser committed Jan 05, 2012 119 120  typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;  Andreas Lauser committed Sep 26, 2011 121 122 123  typedef typename GridView::template Codim<0>::Entity Element; typedef typename GridView::template Codim<0>::Iterator ElementIterator;  Andreas Lauser committed Jan 05, 2012 124 125  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;  Andreas Lauser committed Sep 26, 2011 126 127 128 129  typedef Dumux::MPNCVtkWriter MPNCVtkWriter; enum {  Andreas Lauser committed Jan 05, 2012 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),  Andreas Lauser committed Jan 30, 2012 140  numEq = GET_PROP_VALUE(TypeTag, NumEq)  Andreas Lauser committed Sep 26, 2011 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); };  Bernd Flemisch committed Dec 21, 2011 182 183  if (this->gridView_().comm().size() > 1) dest = this->gridView_().comm().sum(dest);  Andreas Lauser committed Sep 26, 2011 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 void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer) { vtkWriter_->addCurrentSolution(writer); } MPNCVtkWriter *vtkWriter_; }; }  Andreas Lauser committed Feb 13, 2012 201 #include "mpncpropertydefaults.hh"  Andreas Lauser committed Sep 26, 2011 202 203  #endif