Nektar++
Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Private Member Functions | Friends | List of all members
Nektar::RedlichKwongEoS Class Reference

Redlich-Kwong equation of state: p = RT/(1/rho - b) - a/( sqrt(T / Tc) * (1/rho^2 + b/rho) with a = 0.42748 * (R*Tc)^2 / Pc b = 0.08664 * (R*Tc) / Pc. More...

#include <RedlichKwongEoS.h>

Inheritance diagram for Nektar::RedlichKwongEoS:
[legend]

Static Public Member Functions

static EquationOfStateSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession)
 Creates an instance of this class. More...
 

Static Public Attributes

static std::string className
 Name of the class. More...
 

Protected Member Functions

NekDouble v_GetTemperature (const NekDouble &rho, const NekDouble &e) final
 
vec_t v_GetTemperature (const vec_t &rho, const vec_t &e) final
 
NekDouble v_GetPressure (const NekDouble &rho, const NekDouble &e) final
 
vec_t v_GetPressure (const vec_t &rho, const vec_t &e) final
 
NekDouble v_GetEntropy (const NekDouble &rho, const NekDouble &e) final
 
NekDouble v_GetDPDrho_e (const NekDouble &rho, const NekDouble &e) final
 
NekDouble v_GetDPDe_rho (const NekDouble &rho, const NekDouble &e) final
 
NekDouble v_GetEFromRhoP (const NekDouble &rho, const NekDouble &p) final
 
NekDouble v_GetRhoFromPT (const NekDouble &rho, const NekDouble &p) final
 
- Protected Member Functions inherited from Nektar::EquationOfState
 EquationOfState (const LibUtilities::SessionReaderSharedPtr &pSession)
 Constructor. More...
 
 EquationOfState (const NekDouble &gamma, const NekDouble &gasConstant)
 Programmatic Constructor. More...
 
virtual NekDouble v_GetTemperature (const NekDouble &rho, const NekDouble &e)=0
 
virtual vec_t v_GetTemperature (const vec_t &rho, const vec_t &e)=0
 
virtual NekDouble v_GetPressure (const NekDouble &rho, const NekDouble &e)=0
 
virtual vec_t v_GetPressure (const vec_t &rho, const vec_t &e)=0
 
virtual NekDouble v_GetSoundSpeed (const NekDouble &rho, const NekDouble &e)
 
virtual NekDouble v_GetEntropy (const NekDouble &rho, const NekDouble &e)=0
 
virtual NekDouble v_GetDPDrho_e (const NekDouble &rho, const NekDouble &e)=0
 
virtual NekDouble v_GetDPDe_rho (const NekDouble &rho, const NekDouble &e)=0
 
virtual NekDouble v_GetEFromRhoP (const NekDouble &rho, const NekDouble &p)=0
 
virtual NekDouble v_GetRhoFromPT (const NekDouble &rho, const NekDouble &p)=0
 

Protected Attributes

NekDouble m_a
 
NekDouble m_b
 
NekDouble m_Tc
 
NekDouble m_Pc
 
- Protected Attributes inherited from Nektar::EquationOfState
NekDouble m_gamma
 
NekDouble m_gasConstant
 
NekDouble m_gammaMone
 
NekDouble m_gammaMoneOgasConst
 

Private Member Functions

 RedlichKwongEoS (const LibUtilities::SessionReaderSharedPtr &pSession)
 
 ~RedlichKwongEoS (void) override
 
template<class T , typename = typename std::enable_if< std::is_floating_point<T>::value || tinysimd::is_vector_floating_point<T>::value>::type>
Alpha (const T &temp)
 
template<class T , typename = typename std::enable_if< std::is_floating_point<T>::value || tinysimd::is_vector_floating_point<T>::value>::type>
LogTerm (const T &rho)
 
template<class T , typename = typename std::enable_if< std::is_floating_point<T>::value || tinysimd::is_vector_floating_point<T>::value>::type>
GetTemperatureKernel (const T &rho, const T &e)
 
template<class T , typename = typename std::enable_if< std::is_floating_point<T>::value || tinysimd::is_vector_floating_point<T>::value>::type>
GetPressureKernel (const T &rho, const T &e)
 

Friends

class MemoryManager< RedlichKwongEoS >
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::EquationOfState
virtual ~EquationOfState ()
 
NekDouble GetTemperature (const NekDouble &rho, const NekDouble &e)
 Calculate the temperature. More...
 
vec_t GetTemperature (const vec_t &rho, const vec_t &e)
 
NekDouble GetPressure (const NekDouble &rho, const NekDouble &e)
 Calculate the pressure. More...
 
vec_t GetPressure (const vec_t &rho, const vec_t &e)
 
NekDouble GetSoundSpeed (const NekDouble &rho, const NekDouble &e)
 Calculate the sound speed. More...
 
NekDouble GetEntropy (const NekDouble &rho, const NekDouble &e)
 Calculate the entropy. More...
 
NekDouble GetDPDrho_e (const NekDouble &rho, const NekDouble &e)
 Calculate the partial derivative of P(rho,e) with respect to rho. More...
 
NekDouble GetDPDe_rho (const NekDouble &rho, const NekDouble &e)
 Calculate the partial derivative of P(rho,e) with respect to e. More...
 
NekDouble GetEFromRhoP (const NekDouble &rho, const NekDouble &p)
 Obtain the internal energy from rho and P. More...
 
NekDouble GetRhoFromPT (const NekDouble &p, const NekDouble &T)
 Obtain the density from P and T. More...
 

Detailed Description

Redlich-Kwong equation of state: p = RT/(1/rho - b) - a/( sqrt(T / Tc) * (1/rho^2 + b/rho) with a = 0.42748 * (R*Tc)^2 / Pc b = 0.08664 * (R*Tc) / Pc.

Definition at line 53 of file RedlichKwongEoS.h.

Constructor & Destructor Documentation

◆ RedlichKwongEoS()

Nektar::RedlichKwongEoS::RedlichKwongEoS ( const LibUtilities::SessionReaderSharedPtr pSession)
private

Definition at line 47 of file RedlichKwongEoS.cpp.

49 : EquationOfState(pSession)
50{
51 pSession->LoadParameter("Tcrit", m_Tc);
52 pSession->LoadParameter("Pcrit", m_Pc);
53
54 m_a = 0.42748 * m_gasConstant * m_gasConstant * m_Tc * m_Tc / m_Pc;
55 m_b = 0.08664 * m_gasConstant * m_Tc / m_Pc;
56}
EquationOfState(const LibUtilities::SessionReaderSharedPtr &pSession)
Constructor.

References m_a, m_b, Nektar::EquationOfState::m_gasConstant, m_Pc, and m_Tc.

◆ ~RedlichKwongEoS()

Nektar::RedlichKwongEoS::~RedlichKwongEoS ( void  )
inlineoverrideprivate

Definition at line 97 of file RedlichKwongEoS.h.

97{};

Member Function Documentation

◆ Alpha()

template<class T , typename = typename std::enable_if< std::is_floating_point<T>::value || tinysimd::is_vector_floating_point<T>::value>::type>
T Nektar::RedlichKwongEoS::Alpha ( const T &  temp)
inlineprivate

Definition at line 103 of file RedlichKwongEoS.h.

104 {
105 return 1.0 / sqrt(temp / m_Tc);
106 }
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References m_Tc, and tinysimd::sqrt().

Referenced by GetPressureKernel(), v_GetDPDe_rho(), v_GetDPDrho_e(), v_GetEFromRhoP(), v_GetEntropy(), and v_GetRhoFromPT().

◆ create()

static EquationOfStateSharedPtr Nektar::RedlichKwongEoS::create ( const LibUtilities::SessionReaderSharedPtr pSession)
inlinestatic

Creates an instance of this class.

Definition at line 59 of file RedlichKwongEoS.h.

61 {
64 return p;
65 }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::shared_ptr< EquationOfState > EquationOfStateSharedPtr
A shared pointer to an equation of state object.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ GetPressureKernel()

template<class T , typename = typename std::enable_if< std::is_floating_point<T>::value || tinysimd::is_vector_floating_point<T>::value>::type>
T Nektar::RedlichKwongEoS::GetPressureKernel ( const T &  rho,
const T &  e 
)
inlineprivate

Definition at line 163 of file RedlichKwongEoS.h.

164 {
165 T temp = GetTemperatureKernel(rho, e);
166 T oneOrho = 1.0 / rho;
167 T p = m_gasConstant * temp / (oneOrho - m_b) -
168 m_a * Alpha(temp) / (oneOrho * (oneOrho + m_b));
169 return p;
170 }
T GetTemperatureKernel(const T &rho, const T &e)
T Alpha(const T &temp)

References Alpha(), GetTemperatureKernel(), m_a, m_b, Nektar::EquationOfState::m_gasConstant, and CellMLToNektar.cellml_metadata::p.

Referenced by v_GetPressure().

◆ GetTemperatureKernel()

template<class T , typename = typename std::enable_if< std::is_floating_point<T>::value || tinysimd::is_vector_floating_point<T>::value>::type>
T Nektar::RedlichKwongEoS::GetTemperatureKernel ( const T &  rho,
const T &  e 
)
inlineprivate

Definition at line 120 of file RedlichKwongEoS.h.

121 {
122 // First we need to evaluate the log term
123 // ln[1 + b*rho]
124 T logTerm = LogTerm(rho);
125
126 // The temperature can be expressed as an equation in the form
127 // (T^1/2)^3 + A* T^1/2 + B = 0, which we solve iteratively
128 T A = e * (1.0 - m_gamma) / m_gasConstant;
129 T B = -3.0 * m_a / (2.0 * m_b * m_gasConstant) * (m_gamma - 1) *
130 sqrt(m_Tc) * logTerm;
131
132 // Use ideal gas solution as starting guess for iteration
133 T sqrtT = sqrt(e * (m_gamma - 1) / m_gasConstant);
134 // Newton-Raphson iteration to find T^(1/2)
135 T tol = 1e-6;
136 T residual = 1;
137 unsigned int maxIter = 100;
138 unsigned int cnt = 0;
139 while (abs(residual) > tol && cnt < maxIter)
140 {
141 T f = sqrtT * sqrtT * sqrtT + A * sqrtT + B;
142 T df = 3 * sqrtT * sqrtT + A;
143 residual = f / df;
144 sqrtT -= residual;
145 ++cnt;
146 }
147 if (cnt == maxIter)
148 {
149 std::cout << "Newton-Raphson in RedlichKwongEoS::v_GetTemperature "
150 "did not "
151 "converge in "
152 << maxIter << " iterations (residual = " << residual
153 << ")" << std::endl;
154 }
155
156 // Calculate the temperature
157 return sqrtT * sqrtT;
158 }
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298

References tinysimd::abs(), LogTerm(), m_a, m_b, Nektar::EquationOfState::m_gamma, Nektar::EquationOfState::m_gasConstant, m_Tc, and tinysimd::sqrt().

Referenced by GetPressureKernel(), and v_GetTemperature().

◆ LogTerm()

template<class T , typename = typename std::enable_if< std::is_floating_point<T>::value || tinysimd::is_vector_floating_point<T>::value>::type>
T Nektar::RedlichKwongEoS::LogTerm ( const T &  rho)
inlineprivate

Definition at line 112 of file RedlichKwongEoS.h.

113 {
114 return log(1.0 + m_b * rho);
115 }
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:303

References tinysimd::log(), and m_b.

Referenced by GetTemperatureKernel(), v_GetDPDe_rho(), v_GetEFromRhoP(), and v_GetEntropy().

◆ v_GetDPDe_rho()

NekDouble Nektar::RedlichKwongEoS::v_GetDPDe_rho ( const NekDouble rho,
const NekDouble e 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 115 of file RedlichKwongEoS.cpp.

117{
118 NekDouble T = GetTemperature(rho, e);
119 NekDouble alpha = Alpha(T);
120 NekDouble logTerm = LogTerm(rho);
121
122 // First calculate the denominator 1/rho^2 + 2*b/rho - b^2
123 // and sqrt(Alpha) = 1+f_w*(1-sqrt(Tr))
124 NekDouble denom = 1.0 / (rho * rho) + m_b / rho;
125
126 // Compute cv = dedT_rho
127 NekDouble cv = m_gasConstant / (m_gamma - 1) +
128 3 * m_a * alpha * logTerm / (4 * m_b * T);
129
130 // Now we obtain dPdT_rho
131 NekDouble dPdT =
132 m_gasConstant / (1.0 / rho - m_b) + m_a * alpha / (denom * 2 * T);
133
134 // The result is dPde_rho = dPdT_rho / cv
135 return dPdT / cv;
136}
NekDouble GetTemperature(const NekDouble &rho, const NekDouble &e)
Calculate the temperature.
double NekDouble

References Alpha(), Nektar::EquationOfState::GetTemperature(), LogTerm(), m_a, m_b, Nektar::EquationOfState::m_gamma, and Nektar::EquationOfState::m_gasConstant.

◆ v_GetDPDrho_e()

NekDouble Nektar::RedlichKwongEoS::v_GetDPDrho_e ( const NekDouble rho,
const NekDouble e 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 95 of file RedlichKwongEoS.cpp.

97{
98 NekDouble T = GetTemperature(rho, e);
99 NekDouble alpha = Alpha(T);
100 NekDouble dPde = GetDPDe_rho(rho, e);
101
102 // Calculate dPdrho_T
103 NekDouble dPdrho_T =
104 m_gasConstant * T / ((1.0 - m_b * rho) * (1.0 - m_b * rho)) -
105 m_a * alpha * rho * (m_b * rho + 2) /
106 ((1 + m_b * rho) * (1 + m_b * rho));
107
108 // Calculate dedrho_T
109 NekDouble dedrho_T = -3 * m_a * alpha / (2 * (1 + m_b * rho));
110
111 // The result is dPdrho_e = dPdrho_T - dPde_rho * dedrho_T
112 return dPdrho_T - dPde * dedrho_T;
113}
NekDouble GetDPDe_rho(const NekDouble &rho, const NekDouble &e)
Calculate the partial derivative of P(rho,e) with respect to e.

References Alpha(), Nektar::EquationOfState::GetDPDe_rho(), Nektar::EquationOfState::GetTemperature(), m_a, m_b, and Nektar::EquationOfState::m_gasConstant.

◆ v_GetEFromRhoP()

NekDouble Nektar::RedlichKwongEoS::v_GetEFromRhoP ( const NekDouble rho,
const NekDouble p 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 138 of file RedlichKwongEoS.cpp.

140{
141 NekDouble logTerm = LogTerm(rho);
142 // First calculate the temperature, which can be expressed as
143 // (T^1/2)^3 + A* T^1/2 + B = 0
144 NekDouble A, B;
145
146 A = -p * (1.0 / rho - m_b) / m_gasConstant;
147 B = -m_a * sqrt(m_Tc) * (1.0 / rho - m_b) /
148 (1.0 / (rho * rho) + m_b / rho) / m_gasConstant;
149
150 // Use ideal gas solution as starting guess for iteration
151 NekDouble sqrtT = sqrt(p / (rho * (m_gamma - 1)));
152 // Newton-Raphson iteration to find T^(1/2)
153 NekDouble tol = 1e-6;
154 NekDouble maxIter = 100;
155 NekDouble residual = 1;
156 NekDouble f, df;
157 unsigned int cnt = 0;
158 while (abs(residual) > tol && cnt < maxIter)
159 {
160 f = sqrtT * sqrtT * sqrtT + A * sqrtT + B;
161 df = 3 * sqrtT * sqrtT + A;
162 residual = f / df;
163 sqrtT -= residual;
164 ++cnt;
165 }
166 if (cnt == maxIter)
167 {
168 cout << "Newton-Raphson in RedlichKwongEoS::v_GetEFromRhoP did not "
169 "converge in "
170 << maxIter << " iterations (residual = " << residual << ")"
171 << endl;
172 }
173
174 // Calculate T
175 NekDouble T = sqrtT * sqrtT;
176
177 // Calculate internal energy
178 return m_gasConstant * T / (m_gamma - 1) -
179 3 * m_a * Alpha(T) / (2 * m_b) * logTerm;
180}

References tinysimd::abs(), Alpha(), LogTerm(), m_a, m_b, Nektar::EquationOfState::m_gamma, Nektar::EquationOfState::m_gasConstant, m_Tc, CellMLToNektar.cellml_metadata::p, and tinysimd::sqrt().

◆ v_GetEntropy()

NekDouble Nektar::RedlichKwongEoS::v_GetEntropy ( const NekDouble rho,
const NekDouble e 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 80 of file RedlichKwongEoS.cpp.

82{
83 NekDouble T = GetTemperature(rho, e);
84 NekDouble logTerm = LogTerm(rho);
85 // Entropy for an ideal gas
86 NekDouble sIg =
87 m_gasConstant / (m_gamma - 1) * log(T) - m_gasConstant * log(rho);
88
89 NekDouble deltaS = m_gasConstant * log(1 - m_b * rho);
90 deltaS -= m_a * Alpha(T) * logTerm / (2 * m_b * T);
91
92 return sIg + deltaS;
93}

References Alpha(), Nektar::EquationOfState::GetTemperature(), tinysimd::log(), LogTerm(), m_a, m_b, Nektar::EquationOfState::m_gamma, and Nektar::EquationOfState::m_gasConstant.

◆ v_GetPressure() [1/2]

NekDouble Nektar::RedlichKwongEoS::v_GetPressure ( const NekDouble rho,
const NekDouble e 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 69 of file RedlichKwongEoS.cpp.

71{
72 return GetPressureKernel(rho, e);
73}
T GetPressureKernel(const T &rho, const T &e)

References GetPressureKernel().

◆ v_GetPressure() [2/2]

vec_t Nektar::RedlichKwongEoS::v_GetPressure ( const vec_t rho,
const vec_t e 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 75 of file RedlichKwongEoS.cpp.

76{
77 return GetPressureKernel(rho, e);
78}

References GetPressureKernel().

◆ v_GetRhoFromPT()

NekDouble Nektar::RedlichKwongEoS::v_GetRhoFromPT ( const NekDouble rho,
const NekDouble p 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 182 of file RedlichKwongEoS.cpp.

184{
185 // First solve for the compressibility factor Z using the cubic equation
186 // Z^3 + k1 * Z^2 + k2 * Z + k3 = 0
187 // for RedlichKwong:
188 // k1 = -1.0, k2 = A - B - B^2, k3 = -AB
189 // where A = a*alpha(T)*P/(RT)^2, B = bP/(RT)
190 NekDouble A = m_a * Alpha(T) * p / (m_gasConstant * m_gasConstant * T * T);
191 NekDouble B = m_b * p / (m_gasConstant * T);
192
193 NekDouble k1 = -1.0;
194 NekDouble k2 = A - B - B * B;
195 NekDouble k3 = -A * B;
196
197 // Use ideal gas (Z=1) as starting guess for iteration
198 NekDouble Z = 1.0;
199 // Newton-Raphson iteration to find Z
200 NekDouble tol = 1e-6;
201 NekDouble maxIter = 100;
202 NekDouble residual = 1;
203 NekDouble f, df;
204 unsigned int cnt = 0;
205 while (abs(residual) > tol && cnt < maxIter)
206 {
207 f = Z * Z * Z + k1 * Z * Z + k2 * Z + k3;
208 df = 3 * Z * Z + 2 * k1 * Z + k2;
209 residual = f / df;
210 Z -= residual;
211 ++cnt;
212 }
213 if (cnt == maxIter)
214 {
215 cout << "Newton-Raphson in RedlichKwongEoS::v_GetRhoFromPT did not "
216 "converge in "
217 << maxIter << " iterations (residual = " << residual << ")"
218 << endl;
219 }
220
221 // Now calculate rho = p/(ZRT)
222 return p / (Z * m_gasConstant * T);
223}

References tinysimd::abs(), Alpha(), m_a, m_b, Nektar::EquationOfState::m_gasConstant, and CellMLToNektar.cellml_metadata::p.

◆ v_GetTemperature() [1/2]

NekDouble Nektar::RedlichKwongEoS::v_GetTemperature ( const NekDouble rho,
const NekDouble e 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 58 of file RedlichKwongEoS.cpp.

60{
61 return GetTemperatureKernel(rho, e);
62}

References GetTemperatureKernel().

◆ v_GetTemperature() [2/2]

vec_t Nektar::RedlichKwongEoS::v_GetTemperature ( const vec_t rho,
const vec_t e 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 64 of file RedlichKwongEoS.cpp.

65{
66 return GetTemperatureKernel(rho, e);
67}

References GetTemperatureKernel().

Friends And Related Function Documentation

◆ MemoryManager< RedlichKwongEoS >

friend class MemoryManager< RedlichKwongEoS >
friend

Definition at line 1 of file RedlichKwongEoS.h.

Member Data Documentation

◆ className

std::string Nektar::RedlichKwongEoS::className
static
Initial value:
=
"RedlichKwong", RedlichKwongEoS::create,
"Redlich-Kwong equation of state.")
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static EquationOfStateSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
EquationOfStateFactory & GetEquationOfStateFactory()
Declaration of the equation of state factory singleton.

Name of the class.

Definition at line 68 of file RedlichKwongEoS.h.

◆ m_a

NekDouble Nektar::RedlichKwongEoS::m_a
protected

◆ m_b

NekDouble Nektar::RedlichKwongEoS::m_b
protected

◆ m_Pc

NekDouble Nektar::RedlichKwongEoS::m_Pc
protected

Definition at line 74 of file RedlichKwongEoS.h.

Referenced by RedlichKwongEoS().

◆ m_Tc

NekDouble Nektar::RedlichKwongEoS::m_Tc
protected

Definition at line 73 of file RedlichKwongEoS.h.

Referenced by Alpha(), GetTemperatureKernel(), RedlichKwongEoS(), and v_GetEFromRhoP().