gmpMaterialVanGenuchten.h

/************************************************************************
**
** 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_VANGENUCHTEN_H_
#define _GEMA_PLUGIN_HYDRAULIC_MATERIAL_VANGENUCHTEN_H_
 
namespace GmpHydraulicMaterialVanGenuchten
{
//public:
    // -----------------------------------------------------------------------------------------------------------
    double liquidSaturation(double Pc, double Slr, double Sgr, double Pb, double lambda)
    {
        S_TRACE();
 
        // Sl can update according to the capillary pressure
        double Sl = 1.0;
        if (Pc >= Pb)
        {
            Sl = Slr + (1.0 - Slr - Sgr) * std::pow(1.0 + std::pow(Pc / Pb, lambda), (1.0 - lambda) / lambda);
        }
 
        // Apply limits
        Sl = std::min(Sl, 1.0 - Sgr - std::numeric_limits<double>::epsilon()); // Upper bound
        Sl = std::max(Sl, Slr + std::numeric_limits<double>::epsilon());       // Lower bound
 
        // return liquid saturation
        return Sl;
    }
 
    // -----------------------------------------------------------------------------------------------------------
    double derivativeLiquidSaturation(double Sl, double Slr, double Sgr, double Pb, double lambda)
    {
        S_TRACE();
 
        // Compute the effective liquid saturation degree
        double Se = (Sl - Slr) / (1.0 - Sgr - Slr);
 
        // Compute the derivative of the capillary pressure in relation to the liquid saturation
        double m = 1.0 - 1.0 / lambda;
        double v1 = std::pow(Se, (-1.0 / m));
        double v2 = std::pow(v1 - 1.0, -m);
        double dPcdSl = (Pb * (m - 1.0) * v1 * v2) / (m * (Sl - Slr));
 
        return 1.0 / dPcdSl;
    }
 
    // -----------------------------------------------------------------------------------------------------------
    double liquidRelativePermeability(double Se, double lambda)
    {
        S_TRACE();
 
        double m = 1.0 - 1.0 / lambda;
        double klr = std::sqrt(Se) * std::pow(1.0 - std::pow(1.0 - std::pow(Se, 1.0 / m), m), 2);
 
        return klr;
    }
 
    // -----------------------------------------------------------------------------------------------------------
    double gasRelativePermeability(double Se, double lambda)
    {
        S_TRACE();
 
        // Empirical curve fitting parameter
        double m = 1.0 - 1.0 / lambda;
        double kgr = std::pow(1.0 - Se, 1.0 / 3.0) * std::pow(1.0 - std::pow(Se, 1.0 / m), 2.0 * m);
 
        return kgr;
    }
};
#endif