Cantera  3.0.0
Loading...
Searching...
No Matches
WaterProps.cpp
Go to the documentation of this file.
1/**
2 * @file WaterProps.cpp
3 */
4
5// This file is part of Cantera. See License.txt in the top-level directory or
6// at https://cantera.org/license.txt for license and copyright information.
7
11
12namespace Cantera
13{
15 m_waterIAPWS(new WaterPropsIAPWS()),
16 m_own_sub(true)
17{
18}
19
21{
22 if (wptr) {
23 // object in slave mode; it doesn't own its own water evaluator.
24 m_waterIAPWS = wptr->getWater();
25 } else {
27 m_own_sub = true;
28 }
29}
30
32{
33 if (waterIAPWS) {
34 m_waterIAPWS = waterIAPWS;
35 } else {
37 m_own_sub = true;
38 }
39}
40
41WaterProps::~WaterProps()
42{
43 if (m_own_sub) {
44 delete m_waterIAPWS;
45 }
46}
47
48double WaterProps::density_T(double T, double P, int ifunc)
49{
50 static const double Tc = T - 273.15;
51 static const double U1 = 288.9414;
52 static const double U2 = 508929.2;
53 static const double U3 = 68.12963;
54 static const double U4 = -3.9863;
55
56 double tmp1 = Tc + U1;
57 double tmp4 = Tc + U4;
58 double t4t4 = tmp4 * tmp4;
59 double tmp3 = Tc + U3;
60 double rho = 1000. * (1.0 - tmp1*t4t4/(U2 * tmp3));
61
62 // Impose an ideal gas lower bound on rho. We need this to ensure positivity
63 // of rho, even though it is grossly unrepresentative.
64 double rhomin = P / (GasConstant * T);
65 if (rho < rhomin) {
66 rho = rhomin;
67 if (ifunc == 1) {
68 return - rhomin / T;
69 } else if (ifunc == 3) {
70 return rhomin / P;
71 } else if (ifunc == 2) {
72 return 2.0 * rhomin / (T * T);
73 }
74 }
75
76 if (ifunc == 1) {
77 double drhodT = 1000./U2 * (
78 - tmp4 * tmp4 / (tmp3)
79 - tmp1 * 2 * tmp4 / (tmp3)
80 + tmp1 * t4t4 / (tmp3*tmp3)
81 );
82 return drhodT;
83 } else if (ifunc == 3) {
84 return 0.0;
85 } else if (ifunc == 2) {
86 double t3t3 = tmp3 * tmp3;
87 double d2rhodT2 = 1000./U2 *
88 ((-4.0*tmp4-2.0*tmp1)/tmp3 +
89 (2.0*t4t4 + 4.0*tmp1*tmp4)/t3t3
90 - 2.0*tmp1 * t4t4/(t3t3*tmp3));
91 return d2rhodT2;
92 }
93 return rho;
94}
95
96double WaterProps::relEpsilon(double T, double P_pascal, int ifunc)
97{
98 static const double U1 = 3.4279E2;
99 static const double U2 = -5.0866E-3;
100 static const double U3 = 9.4690E-7;
101 static const double U4 = -2.0525;
102 static const double U5 = 3.1159E3;
103 static const double U6 = -1.8289E2;
104 static const double U7 = -8.0325E3;
105 static const double U8 = 4.2142E6;
106 static const double U9 = 2.1417;
107 double T2 = T * T;
108
109 double eps1000 = U1 * exp(U2 * T + U3 * T2);
110 double C = U4 + U5/(U6 + T);
111 double B = U7 + U8/T + U9 * T;
112 double Pbar = P_pascal * 1.0E-5;
113 double tmpBpar = B + Pbar;
114 double tmpB1000 = B + 1000.0;
115 double ltmp = log(tmpBpar/tmpB1000);
116 double epsRel = eps1000 + C * ltmp;
117
118 if (ifunc == 1 || ifunc == 2) {
119 double tmpC = U6 + T;
120 double dCdT = - U5/(tmpC * tmpC);
121 double dBdT = - U8/(T * T) + U9;
122 double deps1000dT = eps1000 * (U2 + 2.0 * U3 * T);
123 double dltmpdT = (dBdT/tmpBpar - dBdT/tmpB1000);
124 if (ifunc == 1) {
125 return deps1000dT + dCdT * ltmp + C * dltmpdT;
126 }
127 double T3 = T2 * T;
128 double d2CdT2 = - 2.0 * dCdT / tmpC;
129 double d2BdT2 = 2.0 * U8 / (T3);
130 double d2ltmpdT2 = (d2BdT2*(1.0/tmpBpar - 1.0/tmpB1000) +
131 dBdT*dBdT*(1.0/(tmpB1000*tmpB1000) - 1.0/(tmpBpar*tmpBpar)));
132 double d2eps1000dT2 = (deps1000dT * (U2 + 2.0 * U3 * T) + eps1000 * (2.0 * U3));
133
134 if (ifunc == 2) {
135 double d2epsReldT2 = (d2eps1000dT2 + d2CdT2 * ltmp + 2.0 * dCdT * dltmpdT
136 + C * d2ltmpdT2);
137 return d2epsReldT2;
138 }
139 }
140 if (ifunc == 3) {
141 double dltmpdP = 1.0E-5 / tmpBpar;
142 return C * dltmpdP;
143 }
144 return epsRel;
145}
146
147double WaterProps::ADebye(double T, double P_input, int ifunc)
148{
149 double psat = satPressure(T);
150 double P;
151 if (psat > P_input) {
152 P = psat;
153 } else {
154 P = P_input;
155 }
156 double epsRelWater = relEpsilon(T, P, 0);
157 double epsilon = epsilon_0 * epsRelWater;
158 double dw = density_IAPWS(T, P);
159 double tmp = sqrt(2.0 * Avogadro * dw / 1000.);
160 double tmp2 = ElectronCharge * ElectronCharge * Avogadro /
161 (epsilon * GasConstant * T);
162 double tmp3 = tmp2 * sqrt(tmp2);
163 double A_Debye = tmp * tmp3 / (8.0 * Pi);
164
165 // dAdT = - 3/2 Ad/T + 1/2 Ad/dw d(dw)/dT - 3/2 Ad/eps d(eps)/dT
166 // dAdT = - 3/2 Ad/T - 1/2 Ad/Vw d(Vw)/dT - 3/2 Ad/eps d(eps)/dT
167 if (ifunc == 1 || ifunc == 2) {
168 double dAdT = - 1.5 * A_Debye / T;
169
170 double depsRelWaterdT = relEpsilon(T, P, 1);
171 dAdT -= A_Debye * (1.5 * depsRelWaterdT / epsRelWater);
172
173 // calculate d(lnV)/dT _constantP, that is, the cte
174 double cte = coeffThermalExp_IAPWS(T, P);
175 double contrib2 = - A_Debye * (0.5 * cte);
176 dAdT += contrib2;
177
178 if (ifunc == 1) {
179 return dAdT;
180 }
181
182 if (ifunc == 2) {
183 // Get the second derivative of the dielectric constant wrt T
184 // -> we will take each of the terms in dAdT and differentiate
185 // it again.
186 double d2AdT2 = 1.5 / T * (A_Debye/T - dAdT);
187 double d2epsRelWaterdT2 = relEpsilon(T, P, 2);
188 d2AdT2 += 1.5 * (- dAdT * depsRelWaterdT / epsRelWater
189 - A_Debye / epsRelWater *
190 (d2epsRelWaterdT2 - depsRelWaterdT * depsRelWaterdT / epsRelWater));
191 double deltaT = -0.1;
192 double Tdel = T + deltaT;
193 double cte_del = coeffThermalExp_IAPWS(Tdel, P);
194 double dctedT = (cte_del - cte) / Tdel;
195 double contrib3 = 0.5 * (-(dAdT * cte) -(A_Debye * dctedT));
196 d2AdT2 += contrib3;
197 return d2AdT2;
198 }
199 }
200
201 // A_Debye = (1/(8 Pi)) sqrt(2 Na dw / 1000)
202 // (e e/(epsilon R T))^3/2
203 //
204 // dAdP = + 1/2 Ad/dw d(dw)/dP - 3/2 Ad/eps d(eps)/dP
205 // dAdP = - 1/2 Ad/Vw d(Vw)/dP - 3/2 Ad/eps d(eps)/dP
206 // dAdP = + 1/2 Ad * kappa - 3/2 Ad/eps d(eps)/dP
207 //
208 // where kappa = - 1/Vw d(Vw)/dP_T (isothermal compressibility)
209 if (ifunc == 3) {
210 double dAdP = 0.0;
211 double depsRelWaterdP = relEpsilon(T, P, 3);
212 dAdP -= A_Debye * (1.5 * depsRelWaterdP / epsRelWater);
213 double kappa = isothermalCompressibility_IAPWS(T,P);
214 dAdP += A_Debye * (0.5 * kappa);
215 return dAdP;
216 }
217 return A_Debye;
218}
219
221{
222 return m_waterIAPWS->psat(T);
223}
224
225double WaterProps::density_IAPWS(double temp, double press)
226{
227 return m_waterIAPWS->density(temp, press, WATER_LIQUID);
228}
229
231{
232 return m_waterIAPWS->density();
233}
234
235double WaterProps::coeffThermalExp_IAPWS(double temp, double press)
236{
237 double dens = m_waterIAPWS->density(temp, press, WATER_LIQUID);
238 if (dens < 0.0) {
239 throw CanteraError("WaterProps::coeffThermalExp_IAPWS",
240 "Unable to solve for density at T = {} and P = {}", temp, press);
241 }
242 return m_waterIAPWS->coeffThermExp();
243}
244
245double WaterProps::isothermalCompressibility_IAPWS(double temp, double press)
246{
247 double dens = m_waterIAPWS->density(temp, press, WATER_LIQUID);
248 if (dens < 0.0) {
249 throw CanteraError("WaterProps::isothermalCompressibility_IAPWS",
250 "Unable to solve for density at T = {} and P = {}", temp, press);
251 }
253}
254
255}
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
Header for a class used to house several approximation routines for properties of water.
Base class for exceptions thrown by Cantera classes.
Class for the liquid water pressure dependent standard state.
Definition PDSS_Water.h:50
WaterPropsIAPWS * getWater()
Get a pointer to a changeable WaterPropsIAPWS object.
Definition PDSS_Water.h:142
Class for calculating the equation of state of water.
double coeffThermExp() const
Returns the coefficient of thermal expansion.
double density(double temperature, double pressure, int phase=-1, double rhoguess=-1.0)
Calculates the density given the temperature and the pressure, and a guess at the density.
double isothermalCompressibility() const
Returns the coefficient of isothermal compressibility for the state of the object.
double psat(double temperature, int waterState=WATER_LIQUID)
This function returns the saturation pressure given the temperature as an input parameter,...
double satPressure(double T)
Returns the saturation pressure given the temperature.
double isothermalCompressibility_IAPWS(double T, double P)
Returns the isothermal compressibility of water.
WaterPropsIAPWS * m_waterIAPWS
Pointer to the WaterPropsIAPWS object.
Definition WaterProps.h:184
double relEpsilon(double T, double P_pascal, int ifunc=0)
Bradley-Pitzer equation for the dielectric constant of water as a function of temperature and pressur...
double density_IAPWS() const
Returns the density of water.
static double density_T(double T, double P, int ifunc)
Simple calculation of water density at atmospheric pressure.
WaterProps()
Default constructor.
double coeffThermalExp_IAPWS(double T, double P)
returns the coefficient of thermal expansion
double ADebye(double T, double P, int ifunc)
ADebye calculates the value of A_Debye as a function of temperature and pressure according to relatio...
bool m_own_sub
true if we own the WaterPropsIAPWS object
Definition WaterProps.h:187
const double Avogadro
Avogadro's Number [number/kmol].
Definition ct_defs.h:81
const double epsilon_0
Permittivity of free space [F/m].
Definition ct_defs.h:137
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition ct_defs.h:120
const double ElectronCharge
Elementary charge [C].
Definition ct_defs.h:90
const double Pi
Pi.
Definition ct_defs.h:68
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
Contains declarations for string manipulation functions within Cantera.