OGS
WaterVapourLatentHeatWithCriticalTemperature.cpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
2// SPDX-License-Identifier: BSD-3-Clause
3
5
6#include <array>
7#include <cmath>
8#include <numeric>
9
12
13namespace MaterialPropertyLib
14{
16constexpr double T_c =
18constexpr double alpha = 1. / 8.;
19constexpr double beta = 1. / 3.;
20constexpr double Delta = 0.79 - beta;
21// Coefficients a_1, a_2, a_4, b_1, ... b_5
22constexpr std::array c = {1989.41582, 11178.45586, 26923.68994,
23 -28989.28947, -19797.03646, 28403.32283,
24 -30382.306422, 15210.380};
25
27 const VariableArray& variable_array,
28 const ParameterLib::SpatialPosition& /*pos*/, const double /*t*/,
29 const double /*dt*/) const
30{
31 const double T = variable_array.temperature;
32
33 if (T >= T_c)
34 {
35 return 0.0;
36 }
37
38 double const tau = (T_c - T) / T_c;
39 double const tau_p2 = tau * tau;
40 double const tau_p3 = tau_p2 * tau;
41 double const tau_p4 = tau_p3 * tau;
42 double const tau_p5 = tau_p4 * tau;
43
44 std::array v = {std::pow(tau, beta),
45 std::pow(tau, beta + Delta),
46 std::pow(tau, 1 - alpha + beta),
47 tau,
48 tau_p2,
49 tau_p3,
50 tau_p4,
51 tau_p5};
52
53 // The formula gives the value in kJ/kg, and the return value is in the
54 // units of J/kg.
55
56 return 1000.0 *
57#if __GNUC__ < 9 || (__GNUC__ == 9 && (__GNUC_MINOR__ < 3))
58 std::inner_product
59#else
60 std::transform_reduce
61#endif
62 (begin(c), end(c), begin(v), 0.);
63}
64
66 VariableArray const& variable_array, Variable const variable,
67 ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
68 double const /*dt*/) const
69{
70 if (variable != Variable::temperature)
71 {
73 "WaterVapourLatentHeatWithCriticalTemperature::dValue is "
74 "implemented "
75 "for "
76 "the derivative with respect to temperature only.");
77 }
78 const double T = variable_array.temperature;
79
80 if (T >= T_c)
81 {
82 return 0.0;
83 }
84
85 constexpr std::array dc = {c[0] * beta,
86 c[1] * (beta + Delta),
87 c[2] * (1 - alpha + beta),
88 c[3],
89 c[4] * 2,
90 c[5] * 3,
91 c[6] * 4,
92 c[7] * 5};
93
94 double const tau = (T_c - T) / T_c;
95 double const tau_p2 = tau * tau;
96 double const tau_p3 = tau_p2 * tau;
97 double const tau_p4 = tau_p3 * tau;
98 std::array v = {std::pow(tau, beta - 1),
99 std::pow(tau, beta + Delta - 1),
100 std::pow(tau, -alpha + beta),
101 1.,
102 tau,
103 tau_p2,
104 tau_p3,
105 tau_p4};
106
107 // The formula gives the value in kJ/kg/K, and the value is return in
108 // the unit of J/kg/K.
109 return -1000.0 *
110#if __GNUC__ < 9 || (__GNUC__ == 9 && (__GNUC_MINOR__ < 3))
111 std::inner_product
112#else
113 std::transform_reduce
114#endif
115 (begin(dc), end(dc), begin(v), 0.) /
116 T_c;
117}
118
119} // namespace MaterialPropertyLib
#define OGS_FATAL(...)
Definition Error.h:19
virtual PropertyDataType value() const
PropertyDataType dValue(VariableArray const &variable_array, Variable const 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