vangenuchten.hh 12.8 KB
 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:  Bernd Flemisch committed Jul 14, 2010 3 /*****************************************************************************  Bernd Flemisch committed Aug 27, 2012 4  * See the file COPYING for full copying permissions. *  Bernd Flemisch committed Jul 14, 2010 5  * *  Andreas Lauser committed Jan 04, 2011 6  * This program is free software: you can redistribute it and/or modify *  Bernd Flemisch committed Jul 14, 2010 7  * it under the terms of the GNU General Public License as published by *  Timo Koch committed Dec 17, 2018 8  * the Free Software Foundation, either version 3 of the License, or *  Andreas Lauser committed Jan 04, 2011 9  * (at your option) any later version. *  Bernd Flemisch committed Jul 14, 2010 10  * *  Andreas Lauser committed Jan 04, 2011 11 12  * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of *  Christoph Grueninger committed Sep 18, 2012 13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *  Andreas Lauser committed Jan 04, 2011 14 15 16 17  * 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 18 19  *****************************************************************************/ /*!  Bernd Flemisch committed Oct 21, 2010 20  * \file  21  * \ingroup Fluidmatrixinteractions  Philipp Nuske committed Nov 11, 2010 22 23  * \brief Implementation of the capillary pressure and * relative permeability <-> saturation relations according to van Genuchten.  Bernd Flemisch committed Jul 14, 2010 24  */  Timo Koch committed Jan 31, 2020 25 26 #ifndef DUMUX_VAN_GENUCHTEN_HH #define DUMUX_VAN_GENUCHTEN_HH  Bernd Flemisch committed Jul 14, 2010 27 28 29 30  #include "vangenuchtenparams.hh" #include  Christoph Grueninger committed Jul 31, 2012 31 32 #include #include  Bernd Flemisch committed Jul 14, 2010 33   Simon Scholz committed Dec 17, 2018 34 35 namespace Dumux {  Bernd Flemisch committed Jul 14, 2010 36 /*!  37  * \ingroup Fluidmatrixinteractions  Philipp Nuske committed Nov 11, 2010 38  * \brief Implementation of the van Genuchten capillary pressure <->  Bernd Flemisch committed Jul 14, 2010 39 40  * saturation relation. This class bundles the "raw" curves * as static members and doesn't concern itself converting  Philipp Nuske committed Nov 11, 2010 41  * absolute to effective saturations and vice versa.  Bernd Flemisch committed Jul 14, 2010 42  *  Philipp Nuske committed Nov 11, 2010 43 44  * For general info: EffToAbsLaw *  Timo Koch committed Jan 31, 2020 45 46 47  * \note Capillary pressure model from van Genuchten (1980), * relative permeability model from Mualem (1976) *  Philipp Nuske committed Nov 11, 2010 48  * \see VanGenuchtenParams  Bernd Flemisch committed Jul 14, 2010 49 50 51 52 53  */ template > class VanGenuchten { public:  Bernd Flemisch committed Dec 16, 2017 54 55  using Params = ParamsT; using Scalar = typename Params::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  * * Van Genuchten's empirical capillary pressure <-> saturation * function is given by  Alexander Kissinger committed Sep 08, 2015 62  * \f$\mathrm{  Bernd Flemisch committed Jul 14, 2010 63  p_C = (\overline{S}_w^{-1/m} - 1)^{1/n}/\alpha  Alexander Kissinger committed Sep 08, 2015 64 65 66  }\f$ * \param swe Effective saturation of the wetting phase \f$\mathrm{\overline{S}_w}\f$ * \param params A container object that is populated with the appropriate coefficients for the respective law.  Kai Wendel committed Nov 19, 2019 67  * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container  Philipp Nuske committed Nov 11, 2010 68  * is constructed accordingly. Afterwards the values are set there, too.  Kilian Weishaupt committed Feb 06, 2017 69 70  * \note Instead of undefined behaviour if swe is not in the valid range, we return a valid number, * by clamping the input. Note that for pc(swe = 0.0) = inf, have a look at RegularizedVanGenuchten if this is a problem.  Bernd Flemisch committed Jul 14, 2010 71  */  Bernd Flemisch committed May 29, 2013 72  static Scalar pc(const Params ¶ms, Scalar swe)  Bernd Flemisch committed Jul 14, 2010 73  {  Kilian Weishaupt committed Feb 06, 2017 74 75 76 77 78 79 80 81  using std::pow; using std::min; using std::max; swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 const Scalar pc = pow(pow(swe, -1.0/params.vgm()) - 1, 1.0/params.vgn())/params.vgAlpha(); return pc;  Bernd Flemisch committed Jul 14, 2010 82 83 84  } /*!  Philipp Nuske committed Nov 11, 2010 85  * \brief The saturation-capillary pressure curve according to van Genuchten.  Bernd Flemisch committed Jul 14, 2010 86 87  * * This is the inverse of the capillary pressure-saturation curve:  Alexander Kissinger committed Sep 08, 2015 88  * \f$\mathrm{  Bernd Flemisch committed Jul 14, 2010 89  \overline{S}_w = {p_C}^{-1} = ((\alpha p_C)^n + 1)^{-m}  Alexander Kissinger committed Sep 08, 2015 90  }\f$  Bernd Flemisch committed Jul 14, 2010 91  *  Alexander Kissinger committed Sep 08, 2015 92 93  * \param pc Capillary pressure \f$\mathrm{p_C}\f$ in \f$\mathrm{[Pa]}\f$ * \param params A container object that is populated with the appropriate coefficients for the respective law.  Kai Wendel committed Nov 19, 2019 94  * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container  Philipp Nuske committed Nov 11, 2010 95  * is constructed accordingly. Afterwards the values are set there, too.  Alexander Kissinger committed Sep 08, 2015 96  * \return The effective saturation of the wetting phase \f$\mathrm{\overline{S}_w}\f$  Kilian Weishaupt committed Feb 06, 2017 97 98  * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, * i.e. sw(pc < 0.0) = 0.0, by clamping the input to the physical bounds.  Bernd Flemisch committed Jul 14, 2010 99  */  Bernd Flemisch committed May 29, 2013 100  static Scalar sw(const Params ¶ms, Scalar pc)  Bernd Flemisch committed Jul 14, 2010 101  {  Kilian Weishaupt committed Feb 06, 2017 102 103 104 105  using std::pow; using std::max; pc = max(pc, 0.0); // the equation below is undefined for negative pcs  Bernd Flemisch committed Jul 14, 2010 106   Kilian Weishaupt committed Feb 06, 2017 107 108  const Scalar sw = pow(pow(params.vgAlpha()*pc, params.vgn()) + 1, -params.vgm()); return sw;  Bernd Flemisch committed Jul 14, 2010 109 110  }  Timo Koch committed Feb 17, 2017 111 112 113 114  /*! * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure * * \param params A container object that is populated with the appropriate coefficients for the respective law.  Kai Wendel committed Nov 19, 2019 115 116  * Therefore, in the (problem specific) spatialParameters first, * the material law is chosen, and then the params container  Timo Koch committed Feb 17, 2017 117 118 119 120 121  * is constructed accordingly. Afterwards the values are set there, too. */ static Scalar endPointPc(const Params ¶ms) { return 0.0; }  Bernd Flemisch committed Jul 14, 2010 122  /*!  Philipp Nuske committed Nov 11, 2010 123 124  * \brief The partial derivative of the capillary * pressure w.r.t. the effective saturation according to van Genuchten.  Bernd Flemisch committed Jul 14, 2010 125 126  * * This is equivalent to  Alexander Kissinger committed Sep 08, 2015 127  * \f$\mathrm{  Bernd Flemisch committed Jul 14, 2010 128 129 130  \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  Alexander Kissinger committed Sep 08, 2015 131  }\f$  Bernd Flemisch committed Nov 11, 2010 132  *  Alexander Kissinger committed Sep 08, 2015 133 134  * \param swe Effective saturation of the wetting phase \f$\mathrm{\overline{S}_w}\f$ * \param params A container object that is populated with the appropriate coefficients for the respective law.  Kai Wendel committed Nov 19, 2019 135 136  * Therefore, in the (problem specific) spatialParameters * first, the material law is chosen, and then the params container  Philipp Nuske committed Nov 11, 2010 137  * is constructed accordingly. Afterwards the values are set there, too.  Kilian Weishaupt committed Feb 06, 2017 138 139 140 141  * * \note Instead of undefined behaviour if swe is not in the valid range, we return a valid number, * by clamping the input. */  Kilian Weishaupt committed Feb 05, 2016 142  static Scalar dpc_dswe(const Params ¶ms, Scalar swe)  Bernd Flemisch committed Jul 14, 2010 143  {  Kilian Weishaupt committed Feb 06, 2017 144 145 146  using std::pow; using std::min; using std::max;  Bernd Flemisch committed Jul 14, 2010 147   Kilian Weishaupt committed Feb 06, 2017 148 149 150  swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 const Scalar powSwe = pow(swe, -1/params.vgm());  151  return - 1.0/params.vgAlpha() * pow(powSwe - 1, 1.0/params.vgn() - 1)/params.vgn()  Kilian Weishaupt committed Feb 06, 2017 152  * powSwe/swe/params.vgm();  Bernd Flemisch committed May 28, 2013 153 154  }  Bernd Flemisch committed Jul 14, 2010 155  /*!  Philipp Nuske committed Nov 11, 2010 156 157  * \brief The partial derivative of the effective * saturation to the capillary pressure according to van Genuchten.  Bernd Flemisch committed Nov 11, 2010 158  *  Alexander Kissinger committed Sep 08, 2015 159 160  * \param pc Capillary pressure \f$\mathrm{p_C}\f$ in \f$\mathrm{[Pa]}\f$ * \param params A container object that is populated with the appropriate coefficients for the respective law.  Kai Wendel committed Nov 19, 2019 161 162  * Therefore, in the (problem specific) spatialParameters * first, the material law is chosen, and then the params container  Philipp Nuske committed Nov 11, 2010 163  * is constructed accordingly. Afterwards the values are set there, too.  Kilian Weishaupt committed Feb 06, 2017 164 165 166  * * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, * by clamping the input.  Bernd Flemisch committed Jul 14, 2010 167  */  Kilian Weishaupt committed Feb 05, 2016 168  static Scalar dswe_dpc(const Params ¶ms, Scalar pc)  Bernd Flemisch committed Jul 14, 2010 169  {  Kilian Weishaupt committed Feb 06, 2017 170 171 172 173  using std::pow; using std::max; pc = max(pc, 0.0); // the equation below is undefined for negative pcs  Bernd Flemisch committed Jul 14, 2010 174   Kilian Weishaupt committed Feb 06, 2017 175 176  const Scalar powAlphaPc = pow(params.vgAlpha()*pc, params.vgn()); return -pow(powAlphaPc + 1, -params.vgm()-1)*params.vgm()*powAlphaPc/pc*params.vgn();  Bernd Flemisch committed May 28, 2013 177 178  }  Bernd Flemisch committed Jul 14, 2010 179 180  /*! * \brief The relative permeability for the wetting phase of  Timo Koch committed Jan 31, 2020 181  * the medium implied by van Genuchten / Mualem  Bernd Flemisch committed Jul 14, 2010 182 183  * parameterization. *  Alexander Kissinger committed Sep 08, 2015 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.  Kai Wendel committed Nov 19, 2019 186 187  * Therefore, in the (problem specific) spatialParameters * first, the material law is chosen, and then the params container  Kilian Weishaupt committed Feb 06, 2017 188 189 190 191 192  * is constructed accordingly. Afterwards the values are set there, too. * * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, * by clamping the input. */  Bernd Flemisch committed May 29, 2013 193  static Scalar krw(const Params ¶ms, Scalar swe)  Bernd Flemisch committed Jul 14, 2010 194  {  Kilian Weishaupt committed Feb 06, 2017 195 196 197 198 199  using std::pow; using std::min; using std::max; swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0  Bernd Flemisch committed Jul 14, 2010 200   Kilian Weishaupt committed Feb 06, 2017 201  const Scalar r = 1.0 - pow(1.0 - pow(swe, 1.0/params.vgm()), params.vgm());  Timo Koch committed Jan 31, 2020 202  return pow(swe, params.vgl())*r*r;  Christoph Grueninger committed Feb 10, 2014 203  }  Bernd Flemisch committed Jul 14, 2010 204 205 206 207  /*! * \brief The derivative of the relative permeability for the * wetting phase in regard to the wetting saturation of the  Timo Koch committed Jan 31, 2020 208  * medium implied by the van Genuchten / Mualem parameterization.  Bernd Flemisch committed Jul 14, 2010 209  *  Alexander Kissinger committed Sep 08, 2015 210 211  * \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.  Kai Wendel committed Nov 19, 2019 212 213  * Therefore, in the (problem specific) spatialParameters * first, the material law is chosen, and then the params container  Philipp Nuske committed Nov 11, 2010 214  * is constructed accordingly. Afterwards the values are set there, too.  Kilian Weishaupt committed Feb 06, 2017 215 216 217  * * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, * by clamping the input.  Bernd Flemisch committed Jul 14, 2010 218  */  Kilian Weishaupt committed Feb 05, 2016 219  static Scalar dkrw_dswe(const Params ¶ms, Scalar swe)  Bernd Flemisch committed Jul 14, 2010 220  {  Kilian Weishaupt committed Feb 06, 2017 221 222 223  using std::pow; using std::min; using std::max;  Bernd Flemisch committed Jul 14, 2010 224   Kilian Weishaupt committed Feb 06, 2017 225 226 227 228  swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0 const Scalar x = 1.0 - pow(swe, 1.0/params.vgm()); const Scalar xToM = pow(x, params.vgm());  Timo Koch committed Jan 31, 2020 229  return (1.0 - xToM)*pow(swe, params.vgl()-1) * ( (1.0 - xToM)*params.vgl() + 2*xToM*(1.0-x)/x );  Christoph Grueninger committed Feb 10, 2014 230  }  Bernd Flemisch committed Jul 14, 2010 231 232 233 234 235 236  /*! * \brief The relative permeability for the non-wetting phase * of the medium implied by van Genuchten's * parameterization. *  Alexander Kissinger committed Sep 08, 2015 237 238  * \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.  Kai Wendel committed Nov 19, 2019 239 240  * Therefore, in the (problem specific) spatialParameters * first, the material law is chosen, and then the params container  Philipp Nuske committed Nov 11, 2010 241  * is constructed accordingly. Afterwards the values are set there, too.  Kilian Weishaupt committed Feb 06, 2017 242 243 244  * * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, * by clamping the input.  Timo Koch committed Jan 31, 2020 245  * \note See e.g. Dury, Fischer, Schulin (1999) for application of Mualem model to non-wetting rel. perm.  Bernd Flemisch committed Jul 14, 2010 246  */  Bernd Flemisch committed May 29, 2013 247  static Scalar krn(const Params ¶ms, Scalar swe)  Bernd Flemisch committed Jul 14, 2010 248  {  Kilian Weishaupt committed Feb 06, 2017 249 250 251  using std::pow; using std::min; using std::max;  Bernd Flemisch committed Jul 14, 2010 252   Kilian Weishaupt committed Feb 06, 2017 253 254  swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0  Timo Koch committed Jan 31, 2020 255  return pow(1 - swe, params.vgl()) * pow(1 - pow(swe, 1.0/params.vgm()), 2*params.vgm());  Christoph Grueninger committed Feb 10, 2014 256  }  Bernd Flemisch committed Jul 14, 2010 257 258 259 260 261 262 263  /*! * \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. *  Alexander Kissinger committed Sep 08, 2015 264 265  * \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.  Kai Wendel committed Nov 19, 2019 266 267  * Therefore, in the (problem specific) spatialParameters * first, the material law is chosen, and then the params container  Philipp Nuske committed Nov 11, 2010 268  * is constructed accordingly. Afterwards the values are set there, too.  Kilian Weishaupt committed Feb 06, 2017 269 270 271  * * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number, * by clamping the input.  Bernd Flemisch committed Jul 14, 2010 272  */  Kilian Weishaupt committed Feb 05, 2016 273  static Scalar dkrn_dswe(const Params ¶ms, Scalar swe)  Bernd Flemisch committed Jul 14, 2010 274  {  Kilian Weishaupt committed Feb 06, 2017 275 276 277  using std::pow; using std::min; using std::max;  Bernd Flemisch committed Jul 14, 2010 278   Kilian Weishaupt committed Feb 06, 2017 279 280  swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0  Timo Koch committed Jan 31, 2020 281 282  const auto sne = 1.0 - swe; const auto x = 1.0 - pow(swe, 1.0/params.vgm());  Theresa Schollenberger committed Mar 31, 2020 283  return -pow(sne, params.vgl()-1.0) * pow(x, 2*params.vgm() - 1.0) * ( params.vgl()*x + 2.0*sne/swe*(1.0 - x) );  Bernd Flemisch committed Jul 14, 2010 284 285 286  } };  Kilian Weishaupt committed Feb 06, 2017 287 288  } // end namespace Dumux  Bernd Flemisch committed Jul 14, 2010 289 290  #endif