Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "PorousFlowBrooksCorey.h"
      11             : 
      12             : namespace PorousFlowBrooksCorey
      13             : {
      14             : Real
      15           6 : effectiveSaturation(Real pc, Real pe, Real lambda)
      16             : {
      17           6 :   if (pc < pe)
      18             :     return 1.0;
      19             :   else
      20           5 :     return std::pow(pc / pe, -lambda);
      21             : }
      22             : 
      23             : Real
      24           5 : dEffectiveSaturation(Real pc, Real pe, Real lambda)
      25             : {
      26           5 :   if (pc < pe)
      27             :     return 0.0;
      28             :   else
      29           4 :     return -lambda * std::pow(pc / pe, -lambda - 1.0) / pe;
      30             : }
      31             : 
      32             : Real
      33           3 : d2EffectiveSaturation(Real pc, Real pe, Real lambda)
      34             : {
      35           3 :   if (pc < pe)
      36             :     return 0.0;
      37             :   else
      38           2 :     return lambda * (lambda + 1.0) * std::pow(pc / pe, -lambda - 2.0) / pe / pe;
      39             : }
      40             : 
      41             : Real
      42     1342788 : capillaryPressure(Real seff, Real pe, Real lambda, Real pc_max)
      43             : {
      44     1342788 :   if (seff >= 1.0)
      45             :     return pe;
      46      584782 :   else if (seff <= 0.0)
      47        2745 :     return pc_max;
      48             :   else
      49      582037 :     return std::min(pe * std::pow(seff, -1.0 / lambda), pc_max);
      50             : }
      51             : 
      52             : Real
      53      480897 : dCapillaryPressure(Real seff, Real pe, Real lambda, Real pc_max)
      54             : {
      55      480897 :   if (seff <= 0.0 || seff > 1.0)
      56             :     return 0.0;
      57             :   else
      58             :   {
      59             :     // Return 0 if pc > pc_max
      60      478173 :     if (capillaryPressure(seff, pe, lambda, pc_max) >= pc_max)
      61             :       return 0.0;
      62             :     else
      63      478063 :       return -pe * std::pow(seff, -1.0 / lambda - 1.0) / lambda;
      64             :   }
      65             : }
      66             : 
      67             : Real
      68      138437 : d2CapillaryPressure(Real seff, Real pe, Real lambda, Real pc_max)
      69             : {
      70      138437 :   if (seff <= 0.0 || seff > 1.0)
      71             :     return 0.0;
      72             :   else
      73             :   {
      74             :     // Return 0 if pc > pc_max
      75      135735 :     if (capillaryPressure(seff, pe, lambda, pc_max) >= pc_max)
      76             :       return 0.0;
      77             :     else
      78      135735 :       return (lambda + 1.0) * pe * std::pow(seff, -1.0 / lambda - 2.0) / lambda / lambda;
      79             :   }
      80             : }
      81             : 
      82             : Real
      83      173602 : dRelativePermeabilityW(Real seff, Real lambda)
      84             : {
      85             :   // Guard against division by zero
      86      173602 :   if (seff <= 0.0 || seff >= 1.0)
      87             :     return 0.0;
      88             : 
      89      165536 :   return (2.0 + 3.0 * lambda) * std::pow(seff, (2.0 + 2.0 * lambda) / lambda) / lambda;
      90             : }
      91             : 
      92             : Real
      93      225898 : dRelativePermeabilityNW(Real seff, Real lambda)
      94             : {
      95             :   // Guard against division by zero
      96      225898 :   if (seff <= 0.0 || seff >= 1.0)
      97             :     return 0.0;
      98             : 
      99      216824 :   return seff * (2.0 + (seff * (2.0 + 3.0 * lambda) - 2.0 * lambda) *
     100      216824 :                            std::pow(1.0 - seff, 2.0 / lambda) / lambda);
     101             : }
     102             : } // namespace PorousFlowBrooksCorey