vangenuchten.hh 9.54 KB
 Bernd Flemisch committed Jul 14, 2010 1 2 3 4 5 6 /***************************************************************************** * Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch * * Institute of Hydraulic Engineering * * University of Stuttgart, Germany * * email: .@iws.uni-stuttgart.de * * *  Andreas Lauser committed Jan 04, 2011 7  * This program is free software: you can redistribute it and/or modify *  Bernd Flemisch committed Jul 14, 2010 8  * it under the terms of the GNU General Public License as published by *  Andreas Lauser committed Jan 04, 2011 9 10  * the Free Software Foundation, either version 2 of the License, or * * (at your option) any later version. *  Bernd Flemisch committed Jul 14, 2010 11  * *  Andreas Lauser committed Jan 04, 2011 12 13 14 15 16 17 18  * 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 . *  Bernd Flemisch committed Jul 14, 2010 19 20  *****************************************************************************/ /*!  Bernd Flemisch committed Oct 21, 2010 21 22  * \file *  Philipp Nuske committed Nov 11, 2010 23 24  * \brief Implementation of the capillary pressure and * relative permeability <-> saturation relations according to van Genuchten.  Bernd Flemisch committed Jul 14, 2010 25 26 27 28 29 30 31 32  */ #ifndef VAN_GENUCHTEN_HH #define VAN_GENUCHTEN_HH #include "vangenuchtenparams.hh" #include  Philipp Nuske committed Nov 11, 2010 33 34 #include #include  Bernd Flemisch committed Jul 14, 2010 35 36 37 38  namespace Dumux { /*!  Philipp Nuske committed Nov 11, 2010 39  * \ingroup fluidmatrixinteractionslaws  Bernd Flemisch committed Jul 14, 2010 40  *  Philipp Nuske committed Nov 11, 2010 41  * \brief Implementation of the van Genuchten capillary pressure <->  Bernd Flemisch committed Jul 14, 2010 42 43  * saturation relation. This class bundles the "raw" curves * as static members and doesn't concern itself converting  Philipp Nuske committed Nov 11, 2010 44  * absolute to effective saturations and vice versa.  Bernd Flemisch committed Jul 14, 2010 45  *  Philipp Nuske committed Nov 11, 2010 46 47 48  * For general info: EffToAbsLaw * * \see VanGenuchtenParams  Bernd Flemisch committed Jul 14, 2010 49 50 51 52 53  */ template > class VanGenuchten { public:  Philipp Nuske committed Nov 11, 2010 54 55  typedef ParamsT Params; typedef typename Params::Scalar Scalar;  Bernd Flemisch committed Jul 14, 2010 56 57  /*!  Philipp Nuske committed Nov 11, 2010 58  * \brief The capillary pressure-saturation curve according to van Genuchten.  Bernd Flemisch committed Jul 14, 2010 59 60 61 62 63 64  * * Van Genuchten's empirical capillary pressure <-> saturation * function is given by * \f[ p_C = (\overline{S}_w^{-1/m} - 1)^{1/n}/\alpha \f]  Philipp Nuske committed Nov 11, 2010 65 66 67 68  * \param Swe Effective saturation of the wetting phase \f$\overline{S}_w\f$ * \param params A container object that is populated with the appropriate coefficients for the respective law. * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too.  Bernd Flemisch committed Jul 14, 2010 69  */  Bernd Flemisch committed Nov 11, 2010 70  static Scalar pC(const Params ¶ms, Scalar Swe)  Bernd Flemisch committed Jul 14, 2010 71  {  Bernd Flemisch committed Nov 11, 2010 72 73  assert(0 <= Swe && Swe <= 1); return pow(pow(Swe, -1.0/params.vgM()) - 1, 1.0/params.vgN())/params.vgAlpha();  Bernd Flemisch committed Jul 14, 2010 74 75 76  } /*!  Philipp Nuske committed Nov 11, 2010 77  * \brief The saturation-capillary pressure curve according to van Genuchten.  Bernd Flemisch committed Jul 14, 2010 78 79 80 81 82 83  * * This is the inverse of the capillary pressure-saturation curve: * \f[ \overline{S}_w = {p_C}^{-1} = ((\alpha p_C)^n + 1)^{-m} \f] *  Philipp Nuske committed Nov 11, 2010 84 85 86 87 88  * \param pC Capillary pressure \f$p_C\f$ * \param params A container object that is populated with the appropriate coefficients for the respective law. * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too. * \return The effective saturation of the wetting phase \f$\overline{S}_w\f$  Bernd Flemisch committed Jul 14, 2010 89 90 91 92 93 94 95 96 97  */ static Scalar Sw(const Params ¶ms, Scalar pC) { assert(pC >= 0); return pow(pow(params.vgAlpha()*pC, params.vgN()) + 1, -params.vgM()); } /*!  Philipp Nuske committed Nov 11, 2010 98 99  * \brief The partial derivative of the capillary * pressure w.r.t. the effective saturation according to van Genuchten.  Bernd Flemisch committed Jul 14, 2010 100 101 102 103 104 105 106  * * This is equivalent to * \f[ \frac{\partial p_C}{\partial \overline{S}_w} = -\frac{1}{\alpha} (\overline{S}_w^{-1/m} - 1)^{1/n - } \overline{S}_w^{-1/m} / \overline{S}_w / m \f]  Bernd Flemisch committed Nov 11, 2010 107  *  Philipp Nuske committed Nov 11, 2010 108 109 110 111  * \param Swe Effective saturation of the wetting phase \f$\overline{S}_w\f$ * \param params A container object that is populated with the appropriate coefficients for the respective law. * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too.  Bernd Flemisch committed Jul 14, 2010 112  */  Bernd Flemisch committed Nov 11, 2010 113  static Scalar dpC_dSw(const Params ¶ms, Scalar Swe)  Bernd Flemisch committed Jul 14, 2010 114  {  Bernd Flemisch committed Nov 11, 2010 115  assert(0 <= Swe && Swe <= 1);  Bernd Flemisch committed Jul 14, 2010 116   Philipp Nuske committed Nov 11, 2010 117 118 119  Scalar powSwe = pow(Swe, -1/params.vgM()); return - 1/params.vgAlpha() * pow(powSwe - 1, 1/params.vgN() - 1)/params.vgN() * powSwe/Swe/params.vgM();  Bernd Flemisch committed Jul 14, 2010 120 121 122  } /*!  Philipp Nuske committed Nov 11, 2010 123 124  * \brief The partial derivative of the effective * saturation to the capillary pressure according to van Genuchten.  Bernd Flemisch committed Nov 11, 2010 125  *  Philipp Nuske committed Nov 11, 2010 126 127 128 129  * \param pC Capillary pressure \f$p_C\f$ * \param params A container object that is populated with the appropriate coefficients for the respective law. * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too.  Bernd Flemisch committed Jul 14, 2010 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144  */ static Scalar dSw_dpC(const Params ¶ms, Scalar pC) { assert(pC >= 0); Scalar powAlphaPc = pow(params.vgAlpha()*pC, params.vgN()); return -pow(powAlphaPc + 1, -params.vgM()-1)* params.vgM()*powAlphaPc/pC*params.vgN(); } /*! * \brief The relative permeability for the wetting phase of * the medium implied by van Genuchten's * parameterization. *  Philipp Nuske committed Nov 11, 2010 145 146 147 148  * \param Swe The mobile saturation of the wetting phase. * \param params A container object that is populated with the appropriate coefficients for the respective law. * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too. */  Bernd Flemisch committed Nov 11, 2010 149  static Scalar krw(const Params ¶ms, Scalar Swe)  Bernd Flemisch committed Jul 14, 2010 150  {  Bernd Flemisch committed Nov 11, 2010 151  assert(0 <= Swe && Swe <= 1);  Bernd Flemisch committed Jul 14, 2010 152   Bernd Flemisch committed Nov 11, 2010 153 154  Scalar r = 1. - pow(1 - pow(Swe, 1/params.vgM()), params.vgM()); return sqrt(Swe)*r*r;  Bernd Flemisch committed Jul 14, 2010 155 156 157 158 159 160 161  }; /*! * \brief The derivative of the relative permeability for the * wetting phase in regard to the wetting saturation of the * medium implied by the van Genuchten parameterization. *  Philipp Nuske committed Nov 11, 2010 162 163 164 165  * \param Swe The mobile saturation of the wetting phase. * \param params A container object that is populated with the appropriate coefficients for the respective law. * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too.  Bernd Flemisch committed Jul 14, 2010 166  */  Bernd Flemisch committed Nov 11, 2010 167  static Scalar dkrw_dSw(const Params ¶ms, Scalar Swe)  Bernd Flemisch committed Jul 14, 2010 168  {  Bernd Flemisch committed Nov 11, 2010 169  assert(0 <= Swe && Swe <= 1);  Bernd Flemisch committed Jul 14, 2010 170   Bernd Flemisch committed Nov 11, 2010 171  const Scalar x = 1 - std::pow(Swe, 1.0/params.vgM());  Bernd Flemisch committed Jul 14, 2010 172  const Scalar xToM = std::pow(x, params.vgM());  173  return (1 - xToM)/std::sqrt(Swe) * ( (1 - xToM)/2 + 2*xToM*(1-x)/x );  Bernd Flemisch committed Jul 14, 2010 174 175 176 177 178 179 180 181  }; /*! * \brief The relative permeability for the non-wetting phase * of the medium implied by van Genuchten's * parameterization. *  Philipp Nuske committed Nov 11, 2010 182 183 184 185  * \param Swe The mobile saturation of the wetting phase. * \param params A container object that is populated with the appropriate coefficients for the respective law. * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too.  Bernd Flemisch committed Jul 14, 2010 186  */  Bernd Flemisch committed Nov 11, 2010 187  static Scalar krn(const Params ¶ms, Scalar Swe)  Bernd Flemisch committed Jul 14, 2010 188  {  Bernd Flemisch committed Nov 11, 2010 189  assert(0 <= Swe && Swe <= 1);  Bernd Flemisch committed Jul 14, 2010 190 191  return  Bernd Flemisch committed Nov 11, 2010 192 193  pow(1 - Swe, 1.0/3) * pow(1 - pow(Swe, 1/params.vgM()), 2*params.vgM());  Bernd Flemisch committed Jul 14, 2010 194 195 196 197 198 199 200 201  }; /*! * \brief The derivative of the relative permeability for the * non-wetting phase in regard to the wetting saturation of * the medium as implied by the van Genuchten * parameterization. *  Philipp Nuske committed Nov 11, 2010 202 203 204 205  * \param Swe The mobile saturation of the wetting phase. * \param params A container object that is populated with the appropriate coefficients for the respective law. * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container * is constructed accordingly. Afterwards the values are set there, too.  Bernd Flemisch committed Jul 14, 2010 206  */  Bernd Flemisch committed Nov 11, 2010 207  static Scalar dkrn_dSw(const Params ¶ms, Scalar Swe)  Bernd Flemisch committed Jul 14, 2010 208  {  Bernd Flemisch committed Nov 11, 2010 209  assert(0 <= Swe && Swe <= 1);  Bernd Flemisch committed Jul 14, 2010 210   Bernd Flemisch committed Nov 11, 2010 211  const Scalar x = std::pow(Swe, 1.0/params.vgM());  Bernd Flemisch committed Jul 14, 2010 212 213  return -std::pow(1 - x, 2*params.vgM())  214  *std::pow(1 - Swe, -2/3)  Bernd Flemisch committed Nov 11, 2010 215  *(1.0/3 + 2*x/Swe);  Bernd Flemisch committed Jul 14, 2010 216 217 218 219 220 221  } }; } #endif