WaterVapourDensity.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "WaterVapourDensity.h"
 
#include <cmath>
 
#include "BaseLib/Error.h"
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/PhysicalConstant.h"
 
namespace MaterialPropertyLib
{
static double saturatedVaporDensity(double const T)
{
    return 1.0e-3 * std::exp(19.819 - 4975.9 / T);
}

static double dsaturatedVaporDensitydT(double const T)
{
    return 4.9759 * std::exp(19.819 - 4975.9 / T) / (T * T);
}
 
static double humidity(double const T, double const p,
                       double const water_density)
{
    return std::exp(
        p / (MaterialLib::PhysicalConstant::SpecificGasConstant::WaterVapour *
             T * water_density));
}
 
PropertyDataType WaterVapourDensity::value(
    VariableArray const& variable_array,
    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
    double const /*dt*/) const
{
    double const p = variable_array.liquid_phase_pressure;
    double const T = variable_array.temperature;
    assert(T >= 273.);
    double const water_density = variable_array.density;
    assert(water_density > 0.);
 
    return humidity(T, p, water_density) * saturatedVaporDensity(T);
}
 
PropertyDataType WaterVapourDensity::dValue(
    VariableArray const& variable_array, Variable const variable,
    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
    double const /*dt*/) const
{
    double const p = variable_array.liquid_phase_pressure;
    double const T = variable_array.temperature;
    assert(T >= 273.);
    double const water_density = variable_array.density;
    assert(water_density > 0.);
 
    if (variable == Variable::temperature)
    {
        double const h = humidity(T, p, water_density);
        double const rho_v = h * saturatedVaporDensity(T);
        double const drho_vS_dT = dsaturatedVaporDensitydT(T);
 
        return h * drho_vS_dT - rho_v * p /
                                    (water_density * T * T *
                                     MaterialLib::PhysicalConstant::
                                         SpecificGasConstant::WaterVapour);
    }
 
    if (variable == Variable::liquid_phase_pressure)
    {
        double const factor =
            1.0 /
            (MaterialLib::PhysicalConstant::SpecificGasConstant::WaterVapour *
             T * water_density);
        double const dhumiditydp = factor * std::exp(factor * p);
 
        return dhumiditydp * saturatedVaporDensity(T);
    }
 
    OGS_FATAL(
        "WaterVapourDensity::dValue is implemented for derivatives with "
        "respect to temperature or liquid_phase_pressure only.");
}
 
}  // namespace MaterialPropertyLib