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 103 of file PengRobinsonEoS.h.

104  {
105  T sqrtAlpha = 1.0 + m_fw * (1.0 - sqrt(temp / m_Tc));
106  return sqrtAlpha * sqrtAlpha;
107  }
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291

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 71 of file PengRobinsonEoS.cpp.

72 {
73  return GetPressureKernel(rho, e);
74 }
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 76 of file PengRobinsonEoS.cpp.

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

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 147 of file PengRobinsonEoS.h.

148  {
149  T temp = GetTemperatureKernel(rho, e);
150  T oneOrho = 1.0 / rho;
151  T p = m_gasConstant * temp / (oneOrho - m_b) -
152  m_a * Alpha(temp) /
153  (oneOrho * oneOrho + 2.0 * m_b * oneOrho - m_b * m_b);
154 
155  return p;
156  }
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.

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

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 122 of file PengRobinsonEoS.h.

123  {
124  // First we need to evaluate the log term
125  // ln[(1/rho + b - b*sqrt(2)) / (1/rho + b + b*sqrt(2))]
126  NekDouble sqrt2 = sqrt(2.0);
127  T logTerm = LogTerm(rho);
128 
129  // The temperature can be expressed as an equation in the form
130  // A * (T^1/2)^2 + B * T^1/2 + C = 0
131 
132  NekDouble A = m_gasConstant / (m_gamma - 1.0);
133  NekDouble f1 = m_a / (m_b * 2.0 * sqrt2);
134  NekDouble f2 = (1.0 + m_fw);
135  T B = -f1 * logTerm / sqrt(m_Tc) * m_fw * f2;
136  T C = f1 * logTerm * f2 * f2 - e;
137 
138  // Solve for T^1/2 (positive root)
139  T sqrtT = (sqrt(B * B - 4 * A * C) - B) / (2 * A);
140  // Calculate the temperature
141  return sqrtT * sqrtT;
142  }
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 113 of file PengRobinsonEoS.h.

114  {
115  return log((1.0 / rho + m_b - m_b * sqrt(2)) /
116  (1.0 / rho + m_b + m_b * sqrt(2)));
117  }
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:300

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 125 of file PengRobinsonEoS.cpp.

127 {
128  NekDouble T = GetTemperature(rho, e);
129  NekDouble logTerm = LogTerm(rho);
130 
131  // First calculate the denominator 1/rho^2 + 2*b/rho - b^2
132  // and sqrt(Alpha) = 1+f_w*(1-sqrt(Tr))
133  NekDouble denom = 1.0 / (rho * rho) + 2.0 * m_b / rho - m_b * m_b;
134  NekDouble sqrtA = sqrt(Alpha(T));
135 
136  // Compute cv = dedT_rho
137  NekDouble cv = m_gasConstant / (m_gamma - 1) -
138  m_a / (2 * m_b * sqrt(8)) * logTerm * (m_fw * (1 + m_fw)) /
139  sqrt(T * m_Tc);
140 
141  // Now we obtain dPdT_rho
142  NekDouble dPdT = m_gasConstant / (1.0 / rho - m_b) +
143  m_a / sqrt(T * m_Tc) * m_fw * sqrtA / denom;
144 
145  // The result is dPde_rho = dPdT_rho / cv
146  return dPdT / cv;
147 }
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 100 of file PengRobinsonEoS.cpp.

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

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

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 81 of file PengRobinsonEoS.cpp.

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

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 177 of file PengRobinsonEoS.cpp.

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

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: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 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