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

virtual NekDouble v_GetTemperature (const NekDouble &rho, const NekDouble &e) override final
 
virtual vec_t v_GetTemperature (const vec_t &rho, const vec_t &e) override final
 
virtual NekDouble v_GetPressure (const NekDouble &rho, const NekDouble &e) override final
 
virtual vec_t v_GetPressure (const vec_t &rho, const vec_t &e) override final
 
virtual NekDouble v_GetEntropy (const NekDouble &rho, const NekDouble &e) override final
 
virtual NekDouble v_GetDPDrho_e (const NekDouble &rho, const NekDouble &e) override final
 
virtual NekDouble v_GetDPDe_rho (const NekDouble &rho, const NekDouble &e) override final
 
virtual NekDouble v_GetEFromRhoP (const NekDouble &rho, const NekDouble &p) override final
 
virtual NekDouble v_GetRhoFromPT (const NekDouble &rho, const NekDouble &p) override 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_GetSoundSpeed (const NekDouble &rho, const NekDouble &e)
 

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)
 
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  )
inlineprivate

Definition at line 106 of file RedlichKwongEoS.h.

106 {};

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 112 of file RedlichKwongEoS.h.

113  {
114  return 1.0 / sqrt(temp / m_Tc);
115  }
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References tinysimd::sqrt().

Referenced by 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 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 172 of file RedlichKwongEoS.h.

173  {
174  T temp = GetTemperatureKernel(rho, e);
175  T oneOrho = 1.0 / rho;
176  T p = m_gasConstant * temp / (oneOrho - m_b) -
177  m_a * Alpha(temp) / (oneOrho * (oneOrho + m_b));
178  return p;
179  }
T GetTemperatureKernel(const T &rho, const T &e)
T Alpha(const T &temp)

References 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 129 of file RedlichKwongEoS.h.

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

References tinysimd::abs(), and tinysimd::sqrt().

Referenced by 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 121 of file RedlichKwongEoS.h.

122  {
123  return log(1.0 + m_b * rho);
124  }
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:303

References tinysimd::log().

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

◆ v_GetDPDe_rho()

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

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 
)
finaloverrideprotectedvirtual

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 
)
finaloverrideprotectedvirtual

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 
)
finaloverrideprotectedvirtual

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 
)
finaloverrideprotectedvirtual

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 
)
finaloverrideprotectedvirtual

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 
)
finaloverrideprotectedvirtual

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 
)
finaloverrideprotectedvirtual

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 
)
finaloverrideprotectedvirtual

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.
Definition: NekFactory.hpp:198
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 RedlichKwongEoS(), and v_GetEFromRhoP().