Nektar++
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 <boost/core/ignore_unused.hpp>
36
37#include "VanDerWaalsEoS.h"
38
39using namespace std;
40
41namespace Nektar
42{
43
44std::string VanDerWaalsEoS::className =
46 "VanDerWaals", VanDerWaalsEoS::create,
47 "Van der Waals equation of state.");
48
51 : EquationOfState(pSession)
52{
53 NekDouble Tcrit, Pcrit;
54 pSession->LoadParameter("Tcrit", Tcrit);
55 pSession->LoadParameter("Pcrit", Pcrit);
56
57 m_a = 27.0 / 64.0 * m_gasConstant * m_gasConstant * Tcrit * Tcrit / Pcrit;
58 m_b = 1.0 / 8.0 * m_gasConstant * Tcrit / Pcrit;
59}
60
62 const NekDouble &e)
63{
64 return GetTemperatureKernel(rho, e);
65}
66
68{
69 return GetTemperatureKernel(rho, e);
70}
71
73 const NekDouble &e)
74{
75 return GetPressureKernel(rho, e);
76}
77
79{
80 return GetPressureKernel(rho, e);
81}
82
84{
85 NekDouble T = GetTemperature(rho, e);
86 NekDouble sIg =
87 m_gasConstant / (m_gamma - 1) * log(T) - m_gasConstant * log(rho);
88
89 return sIg + m_gasConstant * log(1 - m_b * rho);
90}
91
93 const NekDouble &e)
94{
95 NekDouble result;
96
97 result = (m_gamma - 1) * (e + 2 * m_a * rho - m_a * m_b * rho * rho);
98 result = result / ((1 - m_b * rho) * (1 - m_b * rho));
99 result = result - 2 * m_a * rho;
100
101 return result;
102}
103
105 const NekDouble &e)
106{
107 boost::ignore_unused(e);
108 return (m_gamma - 1) / (1.0 / rho - m_b);
109}
110
112 const NekDouble &p)
113{
114 return (p + m_a * rho * rho) * (1.0 / rho - m_b) / (m_gamma - 1) -
115 m_a * rho;
116}
117
119{
120 // First solve for the compressibility factor Z using the cubic equation
121 // Z^3 + k1 * Z^2 + k2 * Z + k3 = 0
122 // for van der Waals:
123 // k1 = -(B+1), k2 = A, k3 = -AB
124 // where A = aP/(RT)^2, B = bP/(RT)
125 NekDouble A = m_a * p / (m_gasConstant * m_gasConstant * T * T);
126 NekDouble B = m_b * p / (m_gasConstant * T);
127
128 NekDouble k1 = -(B + 1);
129 NekDouble k2 = A;
130 NekDouble k3 = -A * B;
131
132 // Use ideal gas (Z=1) as starting guess for iteration
133 NekDouble Z = 1.0;
134 // Newton-Raphson iteration to find Z
135 NekDouble tol = 1e-6;
136 NekDouble maxIter = 100;
137 NekDouble residual = 1;
138 NekDouble f, df;
139 unsigned int cnt = 0;
140 while (abs(residual) > tol && cnt < maxIter)
141 {
142 f = Z * Z * Z + k1 * Z * Z + k2 * Z + k3;
143 df = 3 * Z * Z + 2 * k1 * Z + k2;
144 residual = f / df;
145 Z -= residual;
146 ++cnt;
147 }
148 if (cnt == maxIter)
149 {
150 cout << "Newton-Raphson in VanDerWaalsEoS::v_GetRhoFromPT did not "
151 "converge in "
152 << maxIter << " iterations (residual = " << residual << ")"
153 << endl;
154 }
155
156 // Now calculate rho = p/(ZRT)
157 return p / (Z * m_gasConstant * T);
158}
159} // 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.
Definition: NekFactory.hpp:198
virtual NekDouble v_GetDPDrho_e(const NekDouble &rho, const NekDouble &e) override final
T GetTemperatureKernel(const T &rho, const T &e)
virtual NekDouble v_GetTemperature(const NekDouble &rho, const NekDouble &e) override final
static EquationOfStateSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
virtual NekDouble v_GetEntropy(const NekDouble &rho, const NekDouble &e) override final
T GetPressureKernel(const T &rho, const T &e)
virtual NekDouble v_GetRhoFromPT(const NekDouble &rho, const NekDouble &p) override final
virtual NekDouble v_GetEFromRhoP(const NekDouble &rho, const NekDouble &p) override final
virtual NekDouble v_GetDPDe_rho(const NekDouble &rho, const NekDouble &e) override final
VanDerWaalsEoS(const LibUtilities::SessionReaderSharedPtr &pSession)
virtual NekDouble v_GetPressure(const NekDouble &rho, const NekDouble &e) override final
static std::string className
Name of the class.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
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:298
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:303