Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "HelmholtzFluidProperties.h"
11 : #include "BrentsMethod.h"
12 : #include "libmesh/utility.h"
13 :
14 : InputParameters
15 371 : HelmholtzFluidProperties::validParams()
16 : {
17 371 : InputParameters params = SinglePhaseFluidProperties::validParams();
18 371 : params.addClassDescription("Base class for Helmholtz free energy fluid EOS");
19 371 : return params;
20 0 : }
21 :
22 199 : HelmholtzFluidProperties::HelmholtzFluidProperties(const InputParameters & parameters)
23 199 : : SinglePhaseFluidProperties(parameters)
24 : {
25 199 : }
26 :
27 : Real
28 694 : HelmholtzFluidProperties::rho_from_p_T(Real pressure, Real temperature) const
29 : {
30 : Real density;
31 : // Initial estimate of a bracketing interval for the density
32 694 : Real lower_density = 1.0e-2;
33 694 : Real upper_density = 100.0;
34 :
35 : // The density is found by finding the zero of the pressure
36 6306 : auto pressure_diff = [&pressure, &temperature, this](Real x)
37 6306 : { return this->p_from_rho_T(x, temperature) - pressure; };
38 :
39 694 : BrentsMethod::bracket(pressure_diff, lower_density, upper_density);
40 694 : density = BrentsMethod::root(pressure_diff, lower_density, upper_density);
41 :
42 694 : return density;
43 : }
44 :
45 : void
46 41 : HelmholtzFluidProperties::rho_from_p_T(
47 : Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const
48 : {
49 41 : rho = this->rho_from_p_T(pressure, temperature);
50 :
51 : // Scale the density and temperature
52 41 : const Real delta = rho / criticalDensity();
53 41 : const Real tau = criticalTemperature() / temperature;
54 41 : const Real da_dd = dalpha_ddelta(delta, tau);
55 41 : const Real d2a_dd2 = d2alpha_ddelta2(delta, tau);
56 :
57 41 : drho_dp = molarMass() / (_R * temperature * delta * (2.0 * da_dd + delta * d2a_dd2));
58 41 : drho_dT = rho * (tau * d2alpha_ddeltatau(delta, tau) - da_dd) / temperature /
59 : (2.0 * da_dd + delta * d2a_dd2);
60 41 : }
61 :
62 : Real
63 910377 : HelmholtzFluidProperties::e_from_p_T(Real pressure, Real temperature) const
64 : {
65 : // Require density first
66 910377 : const Real density = rho_from_p_T(pressure, temperature);
67 : // Scale the input density and temperature
68 910377 : const Real delta = density / criticalDensity();
69 910377 : const Real tau = criticalTemperature() / temperature;
70 :
71 910377 : return _R * temperature * tau * dalpha_dtau(delta, tau) / molarMass();
72 : }
73 :
74 : void
75 15 : HelmholtzFluidProperties::e_from_p_T(
76 : Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const
77 : {
78 15 : e = this->e_from_p_T(pressure, temperature);
79 :
80 : // Require density first
81 15 : const Real density = rho_from_p_T(pressure, temperature);
82 : // Scale the input density and temperature
83 15 : const Real delta = density / criticalDensity();
84 15 : const Real tau = criticalTemperature() / temperature;
85 :
86 15 : const Real da_dd = dalpha_ddelta(delta, tau);
87 15 : const Real d2a_dd2 = d2alpha_ddelta2(delta, tau);
88 15 : const Real d2a_ddt = d2alpha_ddeltatau(delta, tau);
89 :
90 15 : de_dp = tau * d2a_ddt / (density * (2.0 * da_dd + delta * d2a_dd2));
91 30 : de_dT = -_R *
92 30 : (delta * tau * d2a_ddt * (da_dd - tau * d2a_ddt) / (2.0 * da_dd + delta * d2a_dd2) +
93 15 : tau * tau * d2alpha_dtau2(delta, tau)) /
94 15 : molarMass();
95 15 : }
96 :
97 : Real
98 900298 : HelmholtzFluidProperties::c_from_p_T(Real pressure, Real temperature) const
99 : {
100 : // Require density first
101 900298 : const Real density = rho_from_p_T(pressure, temperature);
102 : // Scale the input density and temperature
103 900298 : const Real delta = density / criticalDensity();
104 900298 : const Real tau = criticalTemperature() / temperature;
105 :
106 900298 : const Real da_dd = dalpha_ddelta(delta, tau);
107 :
108 900298 : Real w = 2.0 * delta * da_dd + delta * delta * d2alpha_ddelta2(delta, tau);
109 900298 : w -= Utility::pow<2>(delta * da_dd - delta * tau * d2alpha_ddeltatau(delta, tau)) /
110 900298 : (tau * tau * d2alpha_dtau2(delta, tau));
111 :
112 900298 : return std::sqrt(_R * temperature * w / molarMass());
113 : }
114 :
115 : Real
116 900343 : HelmholtzFluidProperties::cp_from_p_T(Real pressure, Real temperature) const
117 : {
118 : // Require density first
119 900343 : const Real density = rho_from_p_T(pressure, temperature);
120 : // Scale the input density and temperature
121 900343 : const Real delta = density / criticalDensity();
122 900343 : const Real tau = criticalTemperature() / temperature;
123 :
124 900343 : const Real da_dd = dalpha_ddelta(delta, tau);
125 :
126 900343 : const Real cp = _R *
127 900343 : (-tau * tau * d2alpha_dtau2(delta, tau) +
128 900343 : Utility::pow<2>(delta * da_dd - delta * tau * d2alpha_ddeltatau(delta, tau)) /
129 900343 : (2.0 * delta * da_dd + delta * delta * d2alpha_ddelta2(delta, tau))) /
130 900343 : molarMass();
131 :
132 900343 : return cp;
133 : }
134 :
135 : Real
136 900343 : HelmholtzFluidProperties::cv_from_p_T(Real pressure, Real temperature) const
137 : {
138 : // Require density first
139 900343 : const Real density = rho_from_p_T(pressure, temperature);
140 : // Scale the input density and temperature
141 900343 : const Real delta = density / criticalDensity();
142 900343 : const Real tau = criticalTemperature() / temperature;
143 :
144 900343 : return -_R * tau * tau * d2alpha_dtau2(delta, tau) / molarMass();
145 : }
146 :
147 : Real
148 900343 : HelmholtzFluidProperties::s_from_p_T(Real pressure, Real temperature) const
149 : {
150 : // Require density first
151 900343 : const Real density = rho_from_p_T(pressure, temperature);
152 : // Scale the input density and temperature
153 900343 : const Real delta = density / criticalDensity();
154 900343 : const Real tau = criticalTemperature() / temperature;
155 :
156 900343 : return _R * (tau * dalpha_dtau(delta, tau) - alpha(delta, tau)) / molarMass();
157 : }
158 :
159 : void
160 0 : HelmholtzFluidProperties::s_from_p_T(
161 : Real pressure, Real temperature, Real & s, Real & ds_dp, Real & ds_dT) const
162 : {
163 0 : s = this->s_from_p_T(pressure, temperature);
164 :
165 : // Require density first
166 0 : const Real density = rho_from_p_T(pressure, temperature);
167 : // Scale the input density and temperature
168 0 : const Real delta = density / criticalDensity();
169 0 : const Real tau = criticalTemperature() / temperature;
170 :
171 0 : const Real da_dd = dalpha_ddelta(delta, tau);
172 0 : const Real da_dt = dalpha_dtau(delta, tau);
173 0 : const Real d2a_dd2 = d2alpha_ddelta2(delta, tau);
174 0 : const Real d2a_dt2 = d2alpha_dtau2(delta, tau);
175 0 : const Real d2a_ddt = d2alpha_ddeltatau(delta, tau);
176 :
177 0 : ds_dp = tau * (d2a_ddt - da_dd) / (density * temperature * (2.0 * da_dd + delta * d2a_dd2));
178 0 : ds_dT = -_R * tau * (da_dt - alpha(delta, tau) + tau * (d2a_dt2 - da_dt)) /
179 0 : (molarMass() * temperature);
180 0 : }
181 :
182 : Real
183 910372 : HelmholtzFluidProperties::h_from_p_T(Real pressure, Real temperature) const
184 : {
185 : // Require density first
186 910372 : const Real density = rho_from_p_T(pressure, temperature);
187 : // Scale the input density and temperature
188 910372 : const Real delta = density / criticalDensity();
189 910372 : const Real tau = criticalTemperature() / temperature;
190 :
191 910372 : return _R * temperature * (tau * dalpha_dtau(delta, tau) + delta * dalpha_ddelta(delta, tau)) /
192 910372 : molarMass();
193 : }
194 :
195 : void
196 9 : HelmholtzFluidProperties::h_from_p_T(
197 : Real pressure, Real temperature, Real & h, Real & dh_dp, Real & dh_dT) const
198 : {
199 9 : h = this->h_from_p_T(pressure, temperature);
200 :
201 : // Require density first
202 9 : const Real density = rho_from_p_T(pressure, temperature);
203 : // Scale the input density and temperature
204 9 : const Real delta = density / criticalDensity();
205 9 : const Real tau = criticalTemperature() / temperature;
206 :
207 9 : const Real da_dd = dalpha_ddelta(delta, tau);
208 9 : const Real d2a_dd2 = d2alpha_ddelta2(delta, tau);
209 9 : const Real d2a_ddt = d2alpha_ddeltatau(delta, tau);
210 :
211 9 : dh_dp = (da_dd + delta * d2a_dd2 + tau * d2a_ddt) / (density * (2.0 * da_dd + delta * d2a_dd2));
212 18 : dh_dT = _R *
213 9 : (delta * da_dd * (1.0 - tau * d2a_ddt / da_dd) * (1.0 - tau * d2a_ddt / da_dd) /
214 18 : (2.0 + delta * d2a_dd2 / da_dd) -
215 9 : tau * tau * d2alpha_dtau2(delta, tau)) /
216 9 : molarMass();
217 9 : }
218 :
219 : Real
220 1 : HelmholtzFluidProperties::T_from_p_h(Real pressure, Real enthalpy) const
221 : {
222 : auto lambda = [&](Real pressure, Real current_T, Real & new_h, Real & dh_dp, Real & dh_dT)
223 5 : { h_from_p_T(pressure, current_T, new_h, dh_dp, dh_dT); };
224 1 : Real T = FluidPropertiesUtils::NewtonSolve(
225 1 : pressure, enthalpy, _T_initial_guess, _tolerance, lambda, name() + "::T_from_p_h")
226 1 : .first;
227 : // check for nans
228 1 : if (std::isnan(T))
229 0 : mooseError("Conversion from enthalpy (h = ",
230 : enthalpy,
231 : ") and pressure (p = ",
232 : pressure,
233 : ") to temperature failed to converge.");
234 1 : return T;
235 : }
236 :
237 : Real
238 89860218 : HelmholtzFluidProperties::p_from_rho_T(Real density, Real temperature) const
239 : {
240 : // Scale the input density and temperature
241 89860218 : const Real delta = density / criticalDensity();
242 89860218 : const Real tau = criticalTemperature() / temperature;
243 :
244 89860218 : return _R * density * temperature * delta * dalpha_ddelta(delta, tau) / molarMass();
245 : }
|
