diff --git a/MaterialLib/MPL/Properties/Enthalpy/GeneralWaterVapourLatentHeat.cpp b/MaterialLib/MPL/Properties/Enthalpy/GeneralWaterVapourLatentHeat.cpp new file mode 100644 index 0000000000000000000000000000000000000000..44d13cb15d9f27d99e05c93f91bff81424983fa5 --- /dev/null +++ b/MaterialLib/MPL/Properties/Enthalpy/GeneralWaterVapourLatentHeat.cpp @@ -0,0 +1,104 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + * Created on March 19, 2021, 11:49 AM + */ + +#include "GeneralWaterVapourLatentHeat.h" + +#include <algorithm> +#include <array> +#include <cmath> + +#include "MaterialLib/MPL/Medium.h" +#include "MaterialLib/MPL/VariableType.h" + +namespace MaterialPropertyLib +{ +static std::array<double, 3> constexpr a = {1989.41582, 11178.45586, + 26923.68994}; +static std::array<double, 3> constexpr a_exp = {0.3333333333333333, 0.79, + 1.2083333333333333}; + +static std::array<double, 5> constexpr b = { + -28989.28947, -19797.03646, 28403.32283, -30382.306422, 15210.380}; + +PropertyDataType GeneralWaterVapourLatentHeat::value( + const VariableArray& variable_array, + const ParameterLib::SpatialPosition& /*pos*/, const double /*t*/, + const double /*dt*/) const +{ + const double T = std::get<double>( + variable_array[static_cast<int>(Variable::temperature)]); + + if (T >= T_c_) + { + return 0.0; + } + + double const tau = (T_c_ - T) / T_c_; + double L_w = 0.0; + for (std::size_t i = 0; i < a.size(); i++) + { + L_w += a[i] * std::pow(tau, a_exp[i]); + } + + double tau_to_power_i = tau; + + std::for_each(b.begin(), b.end(), [&](double const b_i) { + L_w += b_i * tau_to_power_i; // b_i * tau^i + tau_to_power_i *= tau; + }); + + // The formula gives the value in kJ/kg, and the value is return in + // the unit of J/kg. + return 1000.0 * L_w; +} + +PropertyDataType GeneralWaterVapourLatentHeat::dValue( + VariableArray const& variable_array, Variable const primary_variable, + ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/, + double const /*dt*/) const +{ + if (primary_variable == Variable::temperature) + { + const double T = std::get<double>( + variable_array[static_cast<int>(Variable::temperature)]); + + if (T >= T_c_) + { + return 0.0; + } + + double const tau = (T_c_ - T) / T_c_; + double dL_w_dtau = 0.0; + for (std::size_t i = 0; i < a.size(); i++) + { + dL_w_dtau += a[i] * a_exp[i] * std::pow(tau, a_exp[i] - 1.0); + } + + double tau_to_power_i = 1.0; + int exponent_index = 0; + std::for_each(b.begin(), b.end(), [&](double const b_i) { + dL_w_dtau += + (exponent_index + 1) * b_i * tau_to_power_i; // b_i * dau^i/dT + tau_to_power_i *= tau; + exponent_index++; + }); + + // The formula gives the value in kJ/kg/K, and the value is return in + // the unit of J/kg/K. + return -1000.0 * dL_w_dtau / T_c_; + } + + OGS_FATAL( + "GeneralWaterVapourLatentHeat::dValue is implemented for " + "the derivative with respect to temperature only."); +} + +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/Enthalpy/GeneralWaterVapourLatentHeat.h b/MaterialLib/MPL/Properties/Enthalpy/GeneralWaterVapourLatentHeat.h new file mode 100644 index 0000000000000000000000000000000000000000..eb1623fd2f05a1b0b187e12417b71097702ae387 --- /dev/null +++ b/MaterialLib/MPL/Properties/Enthalpy/GeneralWaterVapourLatentHeat.h @@ -0,0 +1,56 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + * Created on March 19, 2021, 11:49 AM + */ + +#pragma once + +#include "MaterialLib/MPL/Property.h" +#include "MaterialLib/PhysicalConstant.h" + +namespace MaterialPropertyLib +{ +class Phase; + +class GeneralWaterVapourLatentHeat final : public Property +{ +public: + explicit GeneralWaterVapourLatentHeat(std::string name) + { + name_ = std::move(name); + } + + void checkScale() const override + { + if (!std::holds_alternative<Phase*>(scale_)) + { + OGS_FATAL( + "The property 'GeneralWaterVapourLatentHeat' is " + "implemented on the 'phase' scale only."); + } + } + + PropertyDataType value(VariableArray const& variable_array, + ParameterLib::SpatialPosition const& pos, + double const t, + double const dt) const override; + + PropertyDataType dValue(VariableArray const& variable_array, + Variable const primary_variable, + ParameterLib::SpatialPosition const& pos, + double const t, double const dt) const override; + +private: + /// Critical temperature. + static double constexpr T_c_ = + 373.92 + MaterialLib::PhysicalConstant::CelsiusZeroInKelvin; + ; +}; + +} // namespace MaterialPropertyLib