gmpMaterialLiakopoulos.h

Go to the documentation of this file.

/************************************************************************
**
** Copyright (C) 2014 by Carlos Augusto Teixera Mendes
** All rights reserved.
**
** This file is part of the "GeMA" software. It's use should respect
** the terms in the license agreement that can be found together
** with this source code.
** It is provided AS IS, with NO WARRANTY OF ANY KIND,
** INCLUDING THE WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR
** A PARTICULAR PURPOSE.
**
************************************************************************/
 
#include <gmTrace.h>
#include <assert.h>
#include <cmath>
#include <limits>
#include <iostream>
 
#ifndef _GEMA_PLUGIN_HYDRAULIC_MATERIAL_LIAKOPOULOS_H_
#define _GEMA_PLUGIN_HYDRAULIC_MATERIAL_LIAKOPOULOS_H_
 
namespace GmpHydraulicMaterialLiakopoulos
{
  // -----------------------------------------------------------------------------------------------------------
  // ! Returns the liquid saturation
  // ! Considers that the capillary pressure Pc is given in kPa
  double liquidSaturation(double Pc)
    {
        S_TRACE();
 
        const double A = 1.9722e-11;
        const double B = 2.4279;
        const double SL_MIN = 0.2;
 
        // Sl can update according to the capillary pressure
        const double Sl = (Pc < 0.0) ? 1.0 : 1.0 - A * std::pow(1000 * Pc, B);
 
        // return liquid saturation
        return std::clamp(Sl, SL_MIN, 1.0);
    }
 
    // -----------------------------------------------------------------------------------------------------------
    double derivativeLiquidSaturation(double Pc)
    {
        S_TRACE();
 
        const double a = 1.9722e-11;
        const double b = 2.4279;
        const double Slmin = 0.2;
 
        // Sl can update according to the capillary pressure
        if (Pc < 0.0)
            return 0.0;
        else
        {
            double pcmax = std::pow((1.0 - Slmin) / a,(1.0 / b));
            return - 1000 * a * b * std::pow(std::min(1000 * Pc, pcmax), b-1);
        }
    }
 
    // -----------------------------------------------------------------------------------------------------------
    double liquidRelativePermeability(double Sl)
    {
        S_TRACE();
 
        const double a = 2.207;
        const double b = 1.0121;
        const double Slmin = 0.2;
 
        double klr;
 
        if (Sl < Slmin)
            return 0.0;
                else if (Sl > 1.0)
            return 1.0;
        else
            return 1.0 - a * std::pow(1.0 - Sl, b);
 
    }
 
    // -----------------------------------------------------------------------------------------------------------
    double derivativeLiquidRelativePermeability(double Sl)
    {
        S_TRACE();
 
        double a = 2.207;
        double b = 1.0121;
        double Slmin = 0.2;
 
        if (Sl < Slmin)
            return 0.0;
                else if(Sl > 1.0)
            return 0.0;
        else
            return a * b * std::pow(1.0 - Sl, b - 1.0);
 
    }
};
#endif