OGS
WaterVapourLatentHeatWithCriticalTemperature.cpp
Go to the documentation of this file.
1
13
14#include <array>
15#include <cmath>
16#include <numeric>
17
20
21namespace MaterialPropertyLib
22{
24constexpr double T_c =
26constexpr double alpha = 1. / 8.;
27constexpr double beta = 1. / 3.;
28constexpr double Delta = 0.79 - beta;
29// Coefficients a_1, a_2, a_4, b_1, ... b_5
30constexpr std::array c = {1989.41582, 11178.45586, 26923.68994,
31 -28989.28947, -19797.03646, 28403.32283,
32 -30382.306422, 15210.380};
33
35 const VariableArray& variable_array,
36 const ParameterLib::SpatialPosition& /*pos*/, const double /*t*/,
37 const double /*dt*/) const
38{
39 const double T = std::get<double>(
40 variable_array[static_cast<int>(Variable::temperature)]);
41
42 if (T >= T_c)
43 {
44 return 0.0;
45 }
46
47 double const tau = (T_c - T) / T_c;
48 double const tau_p2 = tau * tau;
49 double const tau_p3 = tau_p2 * tau;
50 double const tau_p4 = tau_p3 * tau;
51 double const tau_p5 = tau_p4 * tau;
52
53 std::array v = {std::pow(tau, beta),
54 std::pow(tau, beta + Delta),
55 std::pow(tau, 1 - alpha + beta),
56 tau,
57 tau_p2,
58 tau_p3,
59 tau_p4,
60 tau_p5};
61
62 // The formula gives the value in kJ/kg, and the return value is in the
63 // units of J/kg.
64
65 return 1000.0 *
66#if __GNUC__ < 9 || (__GNUC__ == 9 && (__GNUC_MINOR__ < 3))
67 std::inner_product
68#else
69 std::transform_reduce
70#endif
71 (begin(c), end(c), begin(v), 0.);
72}
73
75 VariableArray const& variable_array, Variable const primary_variable,
76 ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
77 double const /*dt*/) const
78{
79 if (primary_variable != Variable::temperature)
80 {
82 "WaterVapourLatentHeatWithCriticalTemperature::dValue is "
83 "implemented "
84 "for "
85 "the derivative with respect to temperature only.");
86 }
87 const double T = std::get<double>(
88 variable_array[static_cast<int>(Variable::temperature)]);
89
90 if (T >= T_c)
91 {
92 return 0.0;
93 }
94
95 constexpr std::array dc = {c[0] * beta,
96 c[1] * (beta + Delta),
97 c[2] * (1 - alpha + beta),
98 c[3],
99 c[4] * 2,
100 c[5] * 3,
101 c[6] * 4,
102 c[7] * 5};
103
104 double const tau = (T_c - T) / T_c;
105 double const tau_p2 = tau * tau;
106 double const tau_p3 = tau_p2 * tau;
107 double const tau_p4 = tau_p3 * tau;
108 std::array v = {std::pow(tau, beta - 1),
109 std::pow(tau, beta + Delta - 1),
110 std::pow(tau, -alpha + beta),
111 1.,
112 tau,
113 tau_p2,
114 tau_p3,
115 tau_p4};
116
117 // The formula gives the value in kJ/kg/K, and the value is return in
118 // the unit of J/kg/K.
119 return -1000.0 *
120#if __GNUC__ < 9 || (__GNUC__ == 9 && (__GNUC_MINOR__ < 3))
121 std::inner_product
122#else
123 std::transform_reduce
124#endif
125 (begin(dc), end(dc), begin(v), 0.) /
126 T_c;
127}
128
129} // namespace MaterialPropertyLib
#define OGS_FATAL(...)
Definition: Error.h:26
virtual PropertyDataType value() const
Definition: Property.cpp:76
PropertyDataType dValue(VariableArray const &variable_array, Variable const primary_variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override
constexpr double CelsiusZeroInKelvin
Zero degrees Celsius in Kelvin.
constexpr double T_c
Critical temperature.
std::variant< double, Eigen::Matrix< double, 2, 1 >, Eigen::Matrix< double, 3, 1 >, Eigen::Matrix< double, 2, 2 >, Eigen::Matrix< double, 3, 3 >, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 >, Eigen::MatrixXd > PropertyDataType
Definition: Property.h:35
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Definition: VariableType.h:110
static const double v