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

Peng-Robinson equation of state: p = RT/(1/rho - b) - a*Alpha(T/Tc) / (1/rho^2 + 2*b/rho - b^2) with a = 0.45724 * (R*Tc)^2 / Pc b = 0.0778 * (R*Tc) / Pc Alpha(T/Tc) = [1 + fw * (1 - sqrt(T/ Tc))]^2 fw = 0.37464 + 1.54226*omega - 0.2699*omega*omega. More...

#include <PengRobinsonEoS.h>

Inheritance diagram for Nektar::PengRobinsonEoS:
[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 GetTemperature (const NekDouble &rho, const NekDouble &e) final
 Calculate the temperature. More...
 
vec_t GetTemperature (const vec_t &rho, const vec_t &e) final
 
NekDouble GetPressure (const NekDouble &rho, const NekDouble &e) final
 Calculate the pressure. More...
 
vec_t 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_GetSoundSpeed (const NekDouble &rho, const NekDouble &e)
 

Protected Attributes

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

Private Member Functions

 PengRobinsonEoS (const LibUtilities::SessionReaderSharedPtr &pSession)
 
 ~PengRobinsonEoS (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< PengRobinsonEoS >
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::EquationOfState
virtual ~EquationOfState ()
 
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

Peng-Robinson equation of state: p = RT/(1/rho - b) - a*Alpha(T/Tc) / (1/rho^2 + 2*b/rho - b^2) with a = 0.45724 * (R*Tc)^2 / Pc b = 0.0778 * (R*Tc) / Pc Alpha(T/Tc) = [1 + fw * (1 - sqrt(T/ Tc))]^2 fw = 0.37464 + 1.54226*omega - 0.2699*omega*omega.

Definition at line 51 of file PengRobinsonEoS.h.

Constructor & Destructor Documentation

◆ PengRobinsonEoS()

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

Definition at line 47 of file PengRobinsonEoS.cpp.

49  : EquationOfState(pSession)
50 {
51  pSession->LoadParameter("Tcrit", m_Tc);
52  pSession->LoadParameter("Pcrit", m_Pc);
53  pSession->LoadParameter("AcentricFactor", m_omega);
54 
55  m_a = 0.45724 * m_gasConstant * m_gasConstant * m_Tc * m_Tc / m_Pc;
56  m_b = 0.0778 * m_gasConstant * m_Tc / m_Pc;
57  m_fw = 0.37464 + 1.54226 * m_omega - 0.2699 * m_omega * m_omega;
58 }
EquationOfState(const LibUtilities::SessionReaderSharedPtr &pSession)
Constructor.

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

◆ ~PengRobinsonEoS()

Nektar::PengRobinsonEoS::~PengRobinsonEoS ( void  )
inlineprivate

Definition at line 97 of file PengRobinsonEoS.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::PengRobinsonEoS::Alpha ( const T &  temp)
inlineprivate

Definition at line 106 of file PengRobinsonEoS.h.

107  {
108  T sqrtAlpha = 1.0 + m_fw * (1.0 - sqrt(temp / m_Tc));
109  return sqrtAlpha * sqrtAlpha;
110  }
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267

References m_fw, 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::PengRobinsonEoS::create ( const LibUtilities::SessionReaderSharedPtr pSession)
inlinestatic

Creates an instance of this class.

Definition at line 57 of file PengRobinsonEoS.h.

59  {
62  return p;
63  }
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.

◆ GetPressure() [1/2]

NekDouble Nektar::PengRobinsonEoS::GetPressure ( const NekDouble rho,
const NekDouble e 
)
finalprotectedvirtual

Calculate the pressure.

Implements Nektar::EquationOfState.

Definition at line 72 of file PengRobinsonEoS.cpp.

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

References GetPressureKernel().

◆ GetPressure() [2/2]

vec_t Nektar::PengRobinsonEoS::GetPressure ( const vec_t rho,
const vec_t e 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 78 of file PengRobinsonEoS.cpp.

80 {
81  return GetPressureKernel(rho, e);
82 }

References GetPressureKernel().

◆ 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::PengRobinsonEoS::GetPressureKernel ( const T &  rho,
const T &  e 
)
inlineprivate

Definition at line 159 of file PengRobinsonEoS.h.

160  {
161  T temp = GetTemperatureKernel(rho, e);
162  T oneOrho = 1.0 / rho;
163  T p = m_gasConstant * temp / (oneOrho - m_b) - m_a * Alpha(temp) /
164  (oneOrho * oneOrho + 2.0 * m_b * oneOrho - m_b * m_b);
165 
166  return p;
167  }
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 GetPressure().

◆ GetTemperature() [1/2]

NekDouble Nektar::PengRobinsonEoS::GetTemperature ( const NekDouble rho,
const NekDouble e 
)
finalprotectedvirtual

Calculate the temperature.

Implements Nektar::EquationOfState.

Definition at line 60 of file PengRobinsonEoS.cpp.

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

References GetTemperatureKernel().

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

◆ GetTemperature() [2/2]

vec_t Nektar::PengRobinsonEoS::GetTemperature ( const vec_t rho,
const vec_t e 
)
finalprotectedvirtual

Implements Nektar::EquationOfState.

Definition at line 66 of file PengRobinsonEoS.cpp.

68 {
69  return GetTemperatureKernel(rho, e);
70 }

References GetTemperatureKernel().

◆ 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::PengRobinsonEoS::GetTemperatureKernel ( const T &  rho,
const T &  e 
)
inlineprivate

Definition at line 131 of file PengRobinsonEoS.h.

132  {
133  // First we need to evaluate the log term
134  // ln[(1/rho + b - b*sqrt(2)) / (1/rho + b + b*sqrt(2))]
135  NekDouble sqrt2 = sqrt(2.0);
136  T logTerm = LogTerm(rho);
137 
138  // The temperature can be expressed as an equation in the form
139  // A * (T^1/2)^2 + B * T^1/2 + C = 0
140 
141  NekDouble A = m_gasConstant / (m_gamma - 1.0);
142  NekDouble f1 = m_a / (m_b * 2.0 * sqrt2);
143  NekDouble f2 = (1.0 + m_fw);
144  T B = -f1 * logTerm / sqrt(m_Tc) * m_fw * f2;
145  T C = f1 * logTerm * f2 * f2 - e;
146 
147  // Solve for T^1/2 (positive root)
148  T sqrtT = (sqrt(B * B - 4 * A * C) - B) / (2 * A);
149  // Calculate the temperature
150  return sqrtT * sqrtT;
151  }
double NekDouble

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

Referenced by GetPressureKernel(), and 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::PengRobinsonEoS::LogTerm ( const T &  rho)
inlineprivate

Definition at line 119 of file PengRobinsonEoS.h.

120  {
121  return log((1.0 / rho + m_b - m_b * sqrt(2)) /
122  (1.0 / rho + m_b + m_b * sqrt(2)));
123  }
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:278

References tinysimd::log(), m_b, and tinysimd::sqrt().

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

◆ v_GetDPDe_rho()

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

Implements Nektar::EquationOfState.

Definition at line 128 of file PengRobinsonEoS.cpp.

130 {
131  NekDouble T = GetTemperature(rho, e);
132  NekDouble logTerm = LogTerm(rho);
133 
134  // First calculate the denominator 1/rho^2 + 2*b/rho - b^2
135  // and sqrt(Alpha) = 1+f_w*(1-sqrt(Tr))
136  NekDouble denom = 1.0 / (rho * rho) + 2.0 * m_b / rho - m_b * m_b;
137  NekDouble sqrtA = sqrt(Alpha(T));
138 
139  // Compute cv = dedT_rho
140  NekDouble cv = m_gasConstant / (m_gamma - 1) -
141  m_a / (2 * m_b * sqrt(8)) * logTerm * (m_fw * (1 + m_fw)) /
142  sqrt(T * m_Tc);
143 
144  // Now we obtain dPdT_rho
145  NekDouble dPdT = m_gasConstant / (1.0 / rho - m_b) +
146  m_a / sqrt(T * m_Tc) * m_fw * sqrtA / denom;
147 
148  // The result is dPde_rho = dPdT_rho / cv
149  return dPdT / cv;
150 }
NekDouble GetTemperature(const NekDouble &rho, const NekDouble &e) final
Calculate the temperature.

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

◆ v_GetDPDrho_e()

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

Implements Nektar::EquationOfState.

Definition at line 103 of file PengRobinsonEoS.cpp.

105 {
106  NekDouble T = GetTemperature(rho, e);
107  NekDouble dPde = GetDPDe_rho(rho, e);
108 
109  // First calculate the denominator 1/rho^2 + 2*b/rho - b^2
110  // and alpha = [1+f_w*(1-sqrt(Tr))]^2
111  NekDouble denom = 1.0 / (rho * rho) + 2.0 * m_b / rho - m_b * m_b;
112  NekDouble alpha = Alpha(T);
113 
114  // Calculate dPdrho_T
115  NekDouble dPdrho_T =
116  m_gasConstant * T / ((1.0 - m_b * rho) * (1.0 - m_b * rho)) -
117  2 * m_a * alpha * rho * (1 + m_b * rho) /
118  ((denom * rho * rho) * (denom * rho * rho));
119 
120  // Calculate dedrho_T
121  NekDouble dedrho_T =
122  -m_a * sqrt(alpha) * (1.0 + m_fw) / (denom * rho * rho);
123 
124  // The result is dPdrho_e = dPdrho_T - dPde_rho * dedrho_T
125  return dPdrho_T - dPde * dedrho_T;
126 }
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(), GetTemperature(), m_a, m_b, m_fw, Nektar::EquationOfState::m_gasConstant, and tinysimd::sqrt().

◆ v_GetEFromRhoP()

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

Implements Nektar::EquationOfState.

Definition at line 152 of file PengRobinsonEoS.cpp.

154 {
155  NekDouble denom = 1.0 / (rho * rho) + 2.0 * m_b / rho - m_b * m_b;
156  NekDouble logTerm = LogTerm(rho);
157  // First we solve for the temperature, which can be expressed as
158  // A * (T^1/2)^2 + B * T^1/2 + C = 0
159  NekDouble A, B, C;
160 
161  A = m_gasConstant / (1.0 / rho - m_b) -
162  (m_a * m_fw * m_fw) / (denom * m_Tc);
163  B = 2 * m_a / denom * m_fw * (1.0 + m_fw) / sqrt(m_Tc);
164  C = -m_a * (1.0 + m_fw) * (1 + m_fw) / denom - p;
165 
166  // Solve for T^1/2 (positive root)
167  NekDouble T = (-B + sqrt(B * B - 4 * A * C)) / (2 * A);
168  T = T * T;
169 
170  // Calculate alpha(T))
171  NekDouble alpha = Alpha(T);
172  // sqrt(Tr)
173  NekDouble sqrtTr = sqrt(T / m_Tc);
174  // Calculate internal energy
175  return m_gasConstant * T / (m_gamma - 1) +
176  m_a / (m_b * sqrt(8)) * logTerm *
177  (alpha + sqrt(alpha) * m_fw * sqrtTr);
178 }

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

◆ v_GetEntropy()

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

Implements Nektar::EquationOfState.

Definition at line 84 of file PengRobinsonEoS.cpp.

86 {
87  NekDouble T = GetTemperature(rho, e);
88  NekDouble logTerm = LogTerm(rho);
89  // Entropy for an ideal gas
90  NekDouble sIg =
91  m_gasConstant / (m_gamma - 1) * log(T) - m_gasConstant * log(rho);
92 
93  // First sqrt(Alpha) = 1+_w*(1-sqrt(Tr)) and sqrt(Tr)
94  NekDouble sqrtA = sqrt(Alpha(T));
95  NekDouble sqrtTr = sqrt(T / m_Tc);
96 
97  NekDouble deltaS = m_gasConstant * log(1 - m_b * rho);
98  deltaS += m_a * sqrtA * m_fw * logTerm * (sqrtTr / T) / (m_b * sqrt(8));
99 
100  return sIg + deltaS;
101 }

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

◆ v_GetRhoFromPT()

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

Implements Nektar::EquationOfState.

Definition at line 180 of file PengRobinsonEoS.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 PengRobinson:
186  // k1 = B-1, k2 = A - 2*B - 3*B^2, k3 = - AB + B^2 + B^3
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 = B - 1.0;
192  NekDouble k2 = A - 2 * B - 3 * B * B;
193  NekDouble k3 = -A * B + B * B + B * B * 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 ((fabs(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 PengRobinsonEoS::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 Alpha(), m_a, m_b, Nektar::EquationOfState::m_gasConstant, and CellMLToNektar.cellml_metadata::p.

Friends And Related Function Documentation

◆ MemoryManager< PengRobinsonEoS >

friend class MemoryManager< PengRobinsonEoS >
friend

Definition at line 1 of file PengRobinsonEoS.h.

Member Data Documentation

◆ className

std::string Nektar::PengRobinsonEoS::className
static
Initial value:
=
"PengRobinson", PengRobinsonEoS::create,
"Peng-Robinson equation of state.")
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:200
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 66 of file PengRobinsonEoS.h.

◆ m_a

NekDouble Nektar::PengRobinsonEoS::m_a
protected

◆ m_b

NekDouble Nektar::PengRobinsonEoS::m_b
protected

◆ m_fw

NekDouble Nektar::PengRobinsonEoS::m_fw
protected

◆ m_omega

NekDouble Nektar::PengRobinsonEoS::m_omega
protected

Definition at line 73 of file PengRobinsonEoS.h.

Referenced by PengRobinsonEoS().

◆ m_Pc

NekDouble Nektar::PengRobinsonEoS::m_Pc
protected

Definition at line 72 of file PengRobinsonEoS.h.

Referenced by PengRobinsonEoS().

◆ m_Tc

NekDouble Nektar::PengRobinsonEoS::m_Tc
protected