OGS
WaterViscosityIAPWS.cpp
Go to the documentation of this file.
1
13#include "WaterViscosityIAPWS.h"
14
15#include <array>
16#include <cmath>
17
18namespace MaterialLib
19{
20namespace Fluid
21{
22static const double Hi[4] = {1.67752, 2.20462, 0.6366564, -0.241605};
23static const double Hij[6][7] = {
24 {0.520094, 0.222531, -0.281378, 0.161913, -0.0325372, 0, 0},
25 {0.0850895, 0.999115, -0.906851, 0.257399, 0, 0, 0},
26 {-1.08374, 1.88797, -0.772479, 0, 0, 0, 0},
27 {-0.289555, 1.26613, -0.489837, 0, 0.0698452, 0, -0.00435673},
28 {0, 0, -0.25704, 0, 0, 0.00872102, 0},
29 {0, 0.120573, 0, 0, 0, 0, -0.000593264}};
30
31static double computeBarMu0Factor(const double barT);
32
33static std::array<double, 6> computeSeriesFactorTForMu1(const double barT);
34static std::array<double, 7> computeSeriesFactorRhoForMu1(const double bar_rho);
35static double computeBarMu1Factor(
36 const std::array<double, 6>& series_factorT,
37 const std::array<double, 7>& series_factorRho);
38
39static double computedBarMu_dbarT(const double barT, double bar_rho);
40static double computedBarMu_dbarRho(const double barT, double bar_rho);
41
42double WaterViscosityIAPWS::getValue(const ArrayType& var_vals) const
43{
44 const double bar_T =
45 var_vals[static_cast<unsigned>(PropertyVariableType::T)] / _ref_T;
46 const double bar_rho =
47 var_vals[static_cast<unsigned>(PropertyVariableType::rho)] / _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 const PropertyVariableType var_type) const
61{
62 const double bar_T =
63 var_vals[static_cast<unsigned>(PropertyVariableType::T)] / _ref_T;
64 const double bar_rho =
65 var_vals[static_cast<unsigned>(PropertyVariableType::rho)] / _ref_rho;
66
67 switch (var_type)
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:
74 return 0.;
75 }
76}
77
78double computeBarMu0Factor(const double barT)
79{
80 double sum_val = 0.;
81 double barT_i = 1.;
82 for (double value : Hi)
83 {
84 sum_val += (value / barT_i);
85 barT_i *= barT;
86 }
87 return sum_val;
88}
89
90std::array<double, 6> computeSeriesFactorTForMu1(const double barT)
91{
92 std::array<double, 6> series_factorT;
93 series_factorT[0] = 1.;
94 const double barT_fac = 1 / barT - 1.0;
95 for (int i = 1; i < 6; i++)
96 {
97 series_factorT[i] = series_factorT[i - 1] * barT_fac;
98 }
99
100 return series_factorT;
101}
102
103std::array<double, 7> computeSeriesFactorRhoForMu1(const double bar_rho)
104{
105 std::array<double, 7> series_factorRho;
106 series_factorRho[0] = 1.;
107 for (int i = 1; i < 7; i++)
108 {
109 series_factorRho[i] = series_factorRho[i - 1] * (bar_rho - 1.0);
110 }
111
112 return series_factorRho;
113}
114
115double computeBarMu1Factor(const std::array<double, 6>& series_factorT,
116 const std::array<double, 7>& series_factorRho)
117{
118 double sum_val = 0.;
119 for (int i = 0; i < 6; i++)
120 {
121 double sum_val_j = 0;
122 for (int j = 0; j < 7; j++)
123 {
124 sum_val_j += Hij[i][j] * series_factorRho[j];
125 }
126 sum_val += series_factorT[i] * sum_val_j;
127 }
128
129 return sum_val;
130}
131
132double computedBarMu_dbarT(const double barT, double bar_rho)
133{
134 const double mu0_factor = computeBarMu0Factor(barT);
135 const double sqrt_barT = std::sqrt(barT);
136
137 double dmu0_factor_dbarT = 0.0;
138 double barT_i = barT * barT;
139 for (int i = 1; i < 4; i++)
140 {
141 dmu0_factor_dbarT -= static_cast<double>(i) * (Hi[i] / barT_i);
142 barT_i *= barT;
143 }
144
145 const double dbar_mu0_dbarT =
146 50. / (mu0_factor * sqrt_barT) -
147 100. * sqrt_barT * dmu0_factor_dbarT / (mu0_factor * mu0_factor);
148
149 const auto& series_factorT = computeSeriesFactorTForMu1(barT);
150 const auto& series_factorRho = computeSeriesFactorRhoForMu1(bar_rho);
151
152 double dmu1_factor_dbarT = 0.0;
153 for (int i = 1; i < 6; i++)
154 {
155 double sum_val_j = 0;
156 for (int j = 0; j < 7; j++)
157 {
158 sum_val_j += Hij[i][j] * series_factorRho[j];
159 }
160 dmu1_factor_dbarT -= static_cast<double>(i) * series_factorT[i - 1] *
161 sum_val_j / (barT * barT);
162 }
163
164 const double mu1_factor =
165 computeBarMu1Factor(series_factorT, series_factorRho);
166 const double dbar_mu1_dbarT =
167 bar_rho * std::exp(bar_rho * mu1_factor) * dmu1_factor_dbarT;
168
169 return dbar_mu0_dbarT * std::exp(bar_rho * mu1_factor) +
170 dbar_mu1_dbarT * 100. * sqrt_barT / mu0_factor;
171}
172
173double computedBarMu_dbarRho(const double barT, double bar_rho)
174{
175 const auto& series_factorT = computeSeriesFactorTForMu1(barT);
176 const auto& series_factorRho = computeSeriesFactorRhoForMu1(bar_rho);
177
178 double dmu1_factor_dbar_rho = 0.0;
179 for (int i = 0; i < 6; i++)
180 {
181 double sum_val_j = 0;
182 for (int j = 1; j < 7; j++)
183 {
184 sum_val_j +=
185 static_cast<double>(j) * Hij[i][j] * series_factorRho[j - 1];
186 }
187 dmu1_factor_dbar_rho += series_factorT[i] * sum_val_j;
188 }
189
190 const double mu0 = 100. * std::sqrt(barT) / computeBarMu0Factor(barT);
191
192 const double mu1_factor =
193 computeBarMu1Factor(series_factorT, series_factorRho);
194 return mu0 * std::exp(bar_rho * mu1_factor) *
195 (mu1_factor + bar_rho * dmu1_factor_dbar_rho);
196}
197
198} // namespace Fluid
199} // namespace MaterialLib
Viscosity model according to Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Wate...
std::array< double, PropertyVariableNumber > ArrayType
Definition: FluidProperty.h:28
double getValue(const ArrayType &var_vals) const override
const double _ref_mu
reference viscosity in Pa.s
double getdValue(const ArrayType &var_vals, const PropertyVariableType var_type) const override
const double _ref_rho
reference density in kg/m^3
const double _ref_T
reference temperature in K
static const double Hij[6][7]
PropertyVariableType
Variable that determine the property.
@ rho
density. For some models, rho substitutes p
static double computeBarMu0Factor(const double barT)
static double computedBarMu_dbarRho(const double barT, double bar_rho)
static double computedBarMu_dbarT(const double barT, double bar_rho)
static std::array< double, 7 > computeSeriesFactorRhoForMu1(const double bar_rho)
static double computeBarMu1Factor(const std::array< double, 6 > &series_factorT, const std::array< double, 7 > &series_factorRho)
static std::array< double, 6 > computeSeriesFactorTForMu1(const double barT)
static const double Hi[4]