OGS
WaterViscosityIAPWS.cpp
Go to the documentation of this file.
1
12#include "WaterViscosityIAPWS.h"
13
14#include <cmath>
15
16#include "BaseLib/Error.h"
18
19namespace MaterialPropertyLib
20{
21static const double Hi[4] = {1.67752, 2.20462, 0.6366564, -0.241605};
22static const double Hij[6][7] = {
23 {0.520094, 0.222531, -0.281378, 0.161913, -0.0325372, 0, 0},
24 {0.0850895, 0.999115, -0.906851, 0.257399, 0, 0, 0},
25 {-1.08374, 1.88797, -0.772479, 0, 0, 0, 0},
26 {-0.289555, 1.26613, -0.489837, 0, 0.0698452, 0, -0.00435673},
27 {0, 0, -0.25704, 0, 0, 0.00872102, 0},
28 {0, 0.120573, 0, 0, 0, 0, -0.000593264}};
29
30static double computeBarMu0Factor(const double barT);
31
32static std::array<double, 6> computeSeriesFactorTForMu1(const double barT);
33static std::array<double, 7> computeSeriesFactorRhoForMu1(const double bar_rho);
34static double computeBarMu1Factor(
35 const std::array<double, 6>& series_factorT,
36 const std::array<double, 7>& series_factorRho);
37
38static double computedBarMu_dbarT(const double barT, double bar_rho);
39static double computedBarMu_dbarRho(const double barT, double bar_rho);
40
42 VariableArray const& variable_array,
43 ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
44 double const /*dt*/) const
45{
46 const double bar_T = variable_array.temperature / ref_T_;
47 const double bar_rho = variable_array.density / ref_rho_;
48
49 const double mu0 = 100. * std::sqrt(bar_T) / computeBarMu0Factor(bar_T);
50
51 const auto& series_factorT = computeSeriesFactorTForMu1(bar_T);
52 const auto& series_factorRho = computeSeriesFactorRhoForMu1(bar_rho);
53 const double mu1 = std::exp(
54 bar_rho * computeBarMu1Factor(series_factorT, series_factorRho));
55
56 return mu0 * mu1 * ref_mu_;
57}
58
60 VariableArray const& variable_array, Variable const variable,
61 ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
62 double const /*dt*/) const
63{
64 const double bar_T = variable_array.temperature / ref_T_;
65 const double bar_rho = variable_array.density / ref_rho_;
66
67 switch (variable)
68 {
70 return ref_mu_ * computedBarMu_dbarT(bar_T, bar_rho) / ref_T_;
72 return ref_mu_ * computedBarMu_dbarRho(bar_T, bar_rho) / ref_rho_;
73 default:
75 "WaterViscosityIAPWS::dValue is implemented for "
76 "derivatives with respect to temperature and liquid density "
77 "only.");
78 }
79}
80
81double computeBarMu0Factor(const double barT)
82{
83 double sum_val = 0.;
84 double barT_i = 1.;
85 for (double value : Hi)
86 {
87 sum_val += (value / barT_i);
88 barT_i *= barT;
89 }
90 return sum_val;
91}
92
93std::array<double, 6> computeSeriesFactorTForMu1(const double barT)
94{
95 std::array<double, 6> series_factorT;
96 series_factorT[0] = 1.;
97 const double barT_fac = 1 / barT - 1.0;
98 for (int i = 1; i < 6; i++)
99 {
100 series_factorT[i] = series_factorT[i - 1] * barT_fac;
101 }
102
103 return series_factorT;
104}
105
106std::array<double, 7> computeSeriesFactorRhoForMu1(const double bar_rho)
107{
108 std::array<double, 7> series_factorRho;
109 series_factorRho[0] = 1.;
110 for (int i = 1; i < 7; i++)
111 {
112 series_factorRho[i] = series_factorRho[i - 1] * (bar_rho - 1.0);
113 }
114
115 return series_factorRho;
116}
117
118double computeBarMu1Factor(const std::array<double, 6>& series_factorT,
119 const std::array<double, 7>& series_factorRho)
120{
121 double sum_val = 0.;
122 for (int i = 0; i < 6; i++)
123 {
124 double sum_val_j = 0;
125 for (int j = 0; j < 7; j++)
126 {
127 sum_val_j += Hij[i][j] * series_factorRho[j];
128 }
129 sum_val += series_factorT[i] * sum_val_j;
130 }
131
132 return sum_val;
133}
134
135double computedBarMu_dbarT(const double barT, double bar_rho)
136{
137 const double mu0_factor = computeBarMu0Factor(barT);
138 const double sqrt_barT = std::sqrt(barT);
139
140 double dmu0_factor_dbarT = 0.0;
141 double barT_i = barT * barT;
142 for (int i = 1; i < 4; i++)
143 {
144 dmu0_factor_dbarT -= static_cast<double>(i) * (Hi[i] / barT_i);
145 barT_i *= barT;
146 }
147
148 const double dbar_mu0_dbarT =
149 50. / (mu0_factor * sqrt_barT) -
150 100. * sqrt_barT * dmu0_factor_dbarT / (mu0_factor * mu0_factor);
151
152 const auto& series_factorT = computeSeriesFactorTForMu1(barT);
153 const auto& series_factorRho = computeSeriesFactorRhoForMu1(bar_rho);
154
155 double dmu1_factor_dbarT = 0.0;
156 for (int i = 1; i < 6; i++)
157 {
158 double sum_val_j = 0;
159 for (int j = 0; j < 7; j++)
160 {
161 sum_val_j += Hij[i][j] * series_factorRho[j];
162 }
163 dmu1_factor_dbarT -= static_cast<double>(i) * series_factorT[i - 1] *
164 sum_val_j / (barT * barT);
165 }
166
167 const double mu1_factor =
168 computeBarMu1Factor(series_factorT, series_factorRho);
169 const double dbar_mu1_dbarT =
170 bar_rho * std::exp(bar_rho * mu1_factor) * dmu1_factor_dbarT;
171
172 return dbar_mu0_dbarT * std::exp(bar_rho * mu1_factor) +
173 dbar_mu1_dbarT * 100. * sqrt_barT / mu0_factor;
174}
175
176double computedBarMu_dbarRho(const double barT, double bar_rho)
177{
178 const auto& series_factorT = computeSeriesFactorTForMu1(barT);
179 const auto& series_factorRho = computeSeriesFactorRhoForMu1(bar_rho);
180
181 double dmu1_factor_dbar_rho = 0.0;
182 for (int i = 0; i < 6; i++)
183 {
184 double sum_val_j = 0;
185 for (int j = 1; j < 7; j++)
186 {
187 sum_val_j +=
188 static_cast<double>(j) * Hij[i][j] * series_factorRho[j - 1];
189 }
190 dmu1_factor_dbar_rho += series_factorT[i] * sum_val_j;
191 }
192
193 const double mu0 = 100. * std::sqrt(barT) / computeBarMu0Factor(barT);
194
195 const double mu1_factor =
196 computeBarMu1Factor(series_factorT, series_factorRho);
197 return mu0 * std::exp(bar_rho * mu1_factor) *
198 (mu1_factor + bar_rho * dmu1_factor_dbar_rho);
199}
200
201} // namespace MaterialPropertyLib
#define OGS_FATAL(...)
Definition Error.h:26
virtual PropertyDataType value() const
Definition Property.cpp:76
static constexpr double ref_mu_
reference viscosity in Pa.s
PropertyDataType dValue(VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override
static constexpr double ref_T_
reference temperature in K
static constexpr double ref_rho_
reference density in kg/m^3
static double computeBarMu0Factor(const double barT)
static const double Hi[4]
static double computedBarMu_dbarT(const double barT, double bar_rho)
static std::array< double, 7 > computeSeriesFactorRhoForMu1(const double bar_rho)
static const double Hij[6][7]
static double computeBarMu1Factor(const std::array< double, 6 > &series_factorT, const std::array< double, 7 > &series_factorRho)
static double computedBarMu_dbarRho(const double barT, double bar_rho)
static std::array< double, 6 > computeSeriesFactorTForMu1(const double barT)
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:31