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=default
 
template<class T , typename = typename std::enable_if< std::is_floating_point_v<T> || tinysimd::is_vector_floating_point_v<T>>::type>
Alpha (const T &temp)
 
template<class T , typename = typename std::enable_if< std::is_floating_point_v<T> || tinysimd::is_vector_floating_point_v<T>>::type>
LogTerm (const T &rho)
 
template<class T , typename = typename std::enable_if< std::is_floating_point_v<T> || tinysimd::is_vector_floating_point_v<T>>::type>
GetTemperatureKernel (const T &rho, const T &e)
 
template<class T , typename = typename std::enable_if< std::is_floating_point_v<T> || tinysimd::is_vector_floating_point_v<T>>::type>
GetPressureKernel (const T &rho, const T &e)
 

Friends

class MemoryManager< RedlichKwongEoS >
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::EquationOfState
virtual ~EquationOfState ()=default
 
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 45 of file RedlichKwongEoS.cpp.

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}
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  )
overrideprivatedefault

Member Function Documentation

◆ Alpha()

template<class T , typename = typename std::enable_if< std::is_floating_point_v<T> || tinysimd::is_vector_floating_point_v<T>>::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:285

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_v<T> || tinysimd::is_vector_floating_point_v<T>>::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_v<T> || tinysimd::is_vector_floating_point_v<T>>::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:289

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_v<T> || tinysimd::is_vector_floating_point_v<T>>::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:294

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 113 of file RedlichKwongEoS.cpp.

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}
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 93 of file RedlichKwongEoS.cpp.

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}
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 136 of file RedlichKwongEoS.cpp.

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}

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 78 of file RedlichKwongEoS.cpp.

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}

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 67 of file RedlichKwongEoS.cpp.

69{
70 return GetPressureKernel(rho, e);
71}
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 73 of file RedlichKwongEoS.cpp.

74{
75 return GetPressureKernel(rho, e);
76}

References GetPressureKernel().

◆ v_GetRhoFromPT()

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

Implements Nektar::EquationOfState.

Definition at line 180 of file RedlichKwongEoS.cpp.

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}

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 56 of file RedlichKwongEoS.cpp.

58{
59 return GetTemperatureKernel(rho, e);
60}

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 62 of file RedlichKwongEoS.cpp.

63{
64 return GetTemperatureKernel(rho, e);
65}

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().