Nektar++
RedlichKwongEoS.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: RedlichKwongEoS.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Redlich-Kwong equation of state
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include "RedlichKwongEoS.h"
36
37namespace Nektar
38{
39
42 "RedlichKwong", RedlichKwongEoS::create,
43 "Redlich-Kwong equation of state.");
44
47 : EquationOfState(pSession)
48{
49 pSession->LoadParameter("Tcrit", m_Tc);
50 pSession->LoadParameter("Pcrit", m_Pc);
51
52 m_a = 0.42748 * m_gasConstant * m_gasConstant * m_Tc * m_Tc / m_Pc;
53 m_b = 0.08664 * m_gasConstant * m_Tc / m_Pc;
54}
55
57 const NekDouble &e)
58{
59 return GetTemperatureKernel(rho, e);
60}
61
63{
64 return GetTemperatureKernel(rho, e);
65}
66
68 const NekDouble &e)
69{
70 return GetPressureKernel(rho, e);
71}
72
74{
75 return GetPressureKernel(rho, e);
76}
77
79 const NekDouble &e)
80{
81 NekDouble T = GetTemperature(rho, e);
82 NekDouble logTerm = LogTerm(rho);
83 // Entropy for an ideal gas
84 NekDouble sIg =
85 m_gasConstant / (m_gamma - 1) * log(T) - m_gasConstant * log(rho);
86
87 NekDouble deltaS = m_gasConstant * log(1 - m_b * rho);
88 deltaS -= m_a * Alpha(T) * logTerm / (2 * m_b * T);
89
90 return sIg + deltaS;
91}
92
94 const NekDouble &e)
95{
96 NekDouble T = GetTemperature(rho, e);
97 NekDouble alpha = Alpha(T);
98 NekDouble dPde = GetDPDe_rho(rho, e);
99
100 // Calculate dPdrho_T
101 NekDouble dPdrho_T =
102 m_gasConstant * T / ((1.0 - m_b * rho) * (1.0 - m_b * rho)) -
103 m_a * alpha * rho * (m_b * rho + 2) /
104 ((1 + m_b * rho) * (1 + m_b * rho));
105
106 // Calculate dedrho_T
107 NekDouble dedrho_T = -3 * m_a * alpha / (2 * (1 + m_b * rho));
108
109 // The result is dPdrho_e = dPdrho_T - dPde_rho * dedrho_T
110 return dPdrho_T - dPde * dedrho_T;
111}
112
114 const NekDouble &e)
115{
116 NekDouble T = GetTemperature(rho, e);
117 NekDouble alpha = Alpha(T);
118 NekDouble logTerm = LogTerm(rho);
119
120 // First calculate the denominator 1/rho^2 + 2*b/rho - b^2
121 // and sqrt(Alpha) = 1+f_w*(1-sqrt(Tr))
122 NekDouble denom = 1.0 / (rho * rho) + m_b / rho;
123
124 // Compute cv = dedT_rho
125 NekDouble cv = m_gasConstant / (m_gamma - 1) +
126 3 * m_a * alpha * logTerm / (4 * m_b * T);
127
128 // Now we obtain dPdT_rho
129 NekDouble dPdT =
130 m_gasConstant / (1.0 / rho - m_b) + m_a * alpha / (denom * 2 * T);
131
132 // The result is dPde_rho = dPdT_rho / cv
133 return dPdT / cv;
134}
135
137 const NekDouble &p)
138{
139 NekDouble logTerm = LogTerm(rho);
140 // First calculate the temperature, which can be expressed as
141 // (T^1/2)^3 + A* T^1/2 + B = 0
142 NekDouble A, B;
143
144 A = -p * (1.0 / rho - m_b) / m_gasConstant;
145 B = -m_a * sqrt(m_Tc) * (1.0 / rho - m_b) /
146 (1.0 / (rho * rho) + m_b / rho) / m_gasConstant;
147
148 // Use ideal gas solution as starting guess for iteration
149 NekDouble sqrtT = sqrt(p / (rho * (m_gamma - 1)));
150 // Newton-Raphson iteration to find T^(1/2)
151 NekDouble tol = 1e-6;
152 NekDouble maxIter = 100;
153 NekDouble residual = 1;
154 NekDouble f, df;
155 unsigned int cnt = 0;
156 while (abs(residual) > tol && cnt < maxIter)
157 {
158 f = sqrtT * sqrtT * sqrtT + A * sqrtT + B;
159 df = 3 * sqrtT * sqrtT + A;
160 residual = f / df;
161 sqrtT -= residual;
162 ++cnt;
163 }
164 if (cnt == maxIter)
165 {
166 cout << "Newton-Raphson in RedlichKwongEoS::v_GetEFromRhoP did not "
167 "converge in "
168 << maxIter << " iterations (residual = " << residual << ")"
169 << endl;
170 }
171
172 // Calculate T
173 NekDouble T = sqrtT * sqrtT;
174
175 // Calculate internal energy
176 return m_gasConstant * T / (m_gamma - 1) -
177 3 * m_a * Alpha(T) / (2 * m_b) * logTerm;
178}
179
181 const NekDouble &T)
182{
183 // First solve for the compressibility factor Z using the cubic equation
184 // Z^3 + k1 * Z^2 + k2 * Z + k3 = 0
185 // for RedlichKwong:
186 // k1 = -1.0, k2 = A - B - B^2, k3 = -AB
187 // where A = a*alpha(T)*P/(RT)^2, B = bP/(RT)
188 NekDouble A = m_a * Alpha(T) * p / (m_gasConstant * m_gasConstant * T * T);
189 NekDouble B = m_b * p / (m_gasConstant * T);
190
191 NekDouble k1 = -1.0;
192 NekDouble k2 = A - B - B * B;
193 NekDouble k3 = -A * B;
194
195 // Use ideal gas (Z=1) as starting guess for iteration
196 NekDouble Z = 1.0;
197 // Newton-Raphson iteration to find Z
198 NekDouble tol = 1e-6;
199 NekDouble maxIter = 100;
200 NekDouble residual = 1;
201 NekDouble f, df;
202 unsigned int cnt = 0;
203 while (abs(residual) > tol && cnt < maxIter)
204 {
205 f = Z * Z * Z + k1 * Z * Z + k2 * Z + k3;
206 df = 3 * Z * Z + 2 * k1 * Z + k2;
207 residual = f / df;
208 Z -= residual;
209 ++cnt;
210 }
211 if (cnt == maxIter)
212 {
213 cout << "Newton-Raphson in RedlichKwongEoS::v_GetRhoFromPT did not "
214 "converge in "
215 << maxIter << " iterations (residual = " << residual << ")"
216 << endl;
217 }
218
219 // Now calculate rho = p/(ZRT)
220 return p / (Z * m_gasConstant * T);
221}
222
223} // namespace Nektar
Encapsulates equations of state allowing us to obtain thermodynamic properties: most relations are in...
NekDouble GetTemperature(const NekDouble &rho, const NekDouble &e)
Calculate the temperature.
NekDouble GetDPDe_rho(const NekDouble &rho, const NekDouble &e)
Calculate the partial derivative of P(rho,e) with respect to e.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
T GetTemperatureKernel(const T &rho, const T &e)
static std::string className
Name of the class.
NekDouble v_GetDPDrho_e(const NekDouble &rho, const NekDouble &e) final
RedlichKwongEoS(const LibUtilities::SessionReaderSharedPtr &pSession)
static EquationOfStateSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
T Alpha(const T &temp)
NekDouble v_GetRhoFromPT(const NekDouble &rho, const NekDouble &p) final
NekDouble v_GetDPDe_rho(const NekDouble &rho, const NekDouble &e) final
NekDouble v_GetTemperature(const NekDouble &rho, const NekDouble &e) final
NekDouble v_GetEntropy(const NekDouble &rho, const NekDouble &e) final
T GetPressureKernel(const T &rho, const T &e)
NekDouble v_GetPressure(const NekDouble &rho, const NekDouble &e) final
NekDouble v_GetEFromRhoP(const NekDouble &rho, const NekDouble &p) final
std::shared_ptr< SessionReader > SessionReaderSharedPtr
EquationOfStateFactory & GetEquationOfStateFactory()
Declaration of the equation of state factory singleton.
tinysimd::simd< NekDouble > vec_t
double NekDouble
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:294
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285