Nektar++
VanDerWaalsEoS.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: VanDerWaalsEoS.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Van der Waals equation of state
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include "VanDerWaalsEoS.h"
36
37namespace Nektar
38{
39
40std::string VanDerWaalsEoS::className =
42 "VanDerWaals", VanDerWaalsEoS::create,
43 "Van der Waals equation of state.");
44
47 : EquationOfState(pSession)
48{
49 NekDouble Tcrit, Pcrit;
50 pSession->LoadParameter("Tcrit", Tcrit);
51 pSession->LoadParameter("Pcrit", Pcrit);
52
53 m_a = 27.0 / 64.0 * m_gasConstant * m_gasConstant * Tcrit * Tcrit / Pcrit;
54 m_b = 1.0 / 8.0 * m_gasConstant * Tcrit / Pcrit;
55}
56
58 const NekDouble &e)
59{
60 return GetTemperatureKernel(rho, e);
61}
62
64{
65 return GetTemperatureKernel(rho, e);
66}
67
69 const NekDouble &e)
70{
71 return GetPressureKernel(rho, e);
72}
73
75{
76 return GetPressureKernel(rho, e);
77}
78
80{
81 NekDouble T = GetTemperature(rho, e);
82 NekDouble sIg =
83 m_gasConstant / (m_gamma - 1) * log(T) - m_gasConstant * log(rho);
84
85 return sIg + m_gasConstant * log(1 - m_b * rho);
86}
87
89 const NekDouble &e)
90{
91 NekDouble result;
92
93 result = (m_gamma - 1) * (e + 2 * m_a * rho - m_a * m_b * rho * rho);
94 result = result / ((1 - m_b * rho) * (1 - m_b * rho));
95 result = result - 2 * m_a * rho;
96
97 return result;
98}
99
101 [[maybe_unused]] const NekDouble &e)
102{
103 return (m_gamma - 1) / (1.0 / rho - m_b);
104}
105
107 const NekDouble &p)
108{
109 return (p + m_a * rho * rho) * (1.0 / rho - m_b) / (m_gamma - 1) -
110 m_a * rho;
111}
112
114{
115 // First solve for the compressibility factor Z using the cubic equation
116 // Z^3 + k1 * Z^2 + k2 * Z + k3 = 0
117 // for van der Waals:
118 // k1 = -(B+1), k2 = A, k3 = -AB
119 // where A = aP/(RT)^2, B = bP/(RT)
120 NekDouble A = m_a * p / (m_gasConstant * m_gasConstant * T * T);
121 NekDouble B = m_b * p / (m_gasConstant * T);
122
123 NekDouble k1 = -(B + 1);
124 NekDouble k2 = A;
125 NekDouble k3 = -A * B;
126
127 // Use ideal gas (Z=1) as starting guess for iteration
128 NekDouble Z = 1.0;
129 // Newton-Raphson iteration to find Z
130 NekDouble tol = 1e-6;
131 NekDouble maxIter = 100;
132 NekDouble residual = 1;
133 NekDouble f, df;
134 unsigned int cnt = 0;
135 while (abs(residual) > tol && cnt < maxIter)
136 {
137 f = Z * Z * Z + k1 * Z * Z + k2 * Z + k3;
138 df = 3 * Z * Z + 2 * k1 * Z + k2;
139 residual = f / df;
140 Z -= residual;
141 ++cnt;
142 }
143 if (cnt == maxIter)
144 {
145 std::cout << "Newton-Raphson in VanDerWaalsEoS::v_GetRhoFromPT did not "
146 "converge in "
147 << maxIter << " iterations (residual = " << residual << ")"
148 << std::endl;
149 }
150
151 // Now calculate rho = p/(ZRT)
152 return p / (Z * m_gasConstant * T);
153}
154} // namespace Nektar
Encapsulates equations of state allowing us to obtain thermodynamic properties: most relations are in...
NekDouble GetTemperature(const NekDouble &rho, const NekDouble &e)
Calculate the temperature.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
NekDouble v_GetPressure(const NekDouble &rho, const NekDouble &e) final
T GetTemperatureKernel(const T &rho, const T &e)
static EquationOfStateSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
NekDouble v_GetEntropy(const NekDouble &rho, const NekDouble &e) final
NekDouble v_GetTemperature(const NekDouble &rho, const NekDouble &e) final
T GetPressureKernel(const T &rho, const T &e)
NekDouble v_GetEFromRhoP(const NekDouble &rho, const NekDouble &p) final
VanDerWaalsEoS(const LibUtilities::SessionReaderSharedPtr &pSession)
NekDouble v_GetDPDe_rho(const NekDouble &rho, const NekDouble &e) final
NekDouble v_GetRhoFromPT(const NekDouble &rho, const NekDouble &p) final
NekDouble v_GetDPDrho_e(const NekDouble &rho, const NekDouble &e) final
static std::string className
Name of the class.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
EquationOfStateFactory & GetEquationOfStateFactory()
Declaration of the equation of state factory singleton.
tinysimd::simd< NekDouble > vec_t
double NekDouble
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:294