1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221 | // *****************************************************************************
/*!
\file src/PDE/EoS/ThermallyPerfectGas.cpp
\copyright 2012-2015 J. Bakosi,
2016-2018 Los Alamos National Security, LLC.,
2019-2021 Triad National Security, LLC.
All rights reserved. See the LICENSE file for details.
\brief Thermally perfect gas equation of state
\details This file defines functions for the thermally perfect gas equation
of state for the compressible flow equations.
*/
// *****************************************************************************
#include <cmath>
#include <iostream>
#include "EoS/ThermallyPerfectGas.hpp"
using inciter::ThermallyPerfectGas;
ThermallyPerfectGas::ThermallyPerfectGas(
tk::real gamma,
tk::real R,
std::vector< tk::real > cp_coeff) :
m_gamma(gamma),
m_R(R),
m_cp_coeff(cp_coeff)
// *************************************************************************
// Constructor
//! \param[in] gamma Ratio of specific heats
//! \param[in] R gas constant
//! \param[in] cp_coeff NASA Glenn polynomials coefficients for cp fit
// *************************************************************************
{ }
tk::real
ThermallyPerfectGas::density(
tk::real pr,
tk::real temp ) const
// *************************************************************************
//! \brief Calculate density from the material pressure and temperature
//! using the stiffened-gas equation of state
//! \param[in] pr Material pressure
//! \param[in] temp Material temperature
//! \return Material density calculated using the stiffened-gas EoS
// *************************************************************************
{
tk::real R = m_R;
tk::real rho = pr / (R * temp);
return rho;
}
tk::real
ThermallyPerfectGas::pressure(
tk::real rho,
tk::real u,
tk::real v,
tk::real w,
tk::real rhoE,
tk::real,
std::size_t,
const std::array< std::array< tk::real, 3 >, 3 >& ) const
// *************************************************************************
//! \brief Calculate pressure from the material density, momentum and total
//! energy using the thermally perfect gas equation of state
//! \param[in] rho density
//! \param[in] u X-velocity
//! \param[in] v Y-velocity
//! \param[in] w Z-velocity
//! \param[in] rhoE total energy
//! \return Pressure calculated using the thermally perfect gas EOS
// *************************************************************************
{
tk::real R = m_R;
tk::real temp = temperature(rho, u, v, w, rhoE);
tk::real pres = rho * R * temp;
return pres;
}
std::array< std::array< tk::real, 3 >, 3 >
ThermallyPerfectGas::CauchyStress(
tk::real,
tk::real,
tk::real,
tk::real,
tk::real,
tk::real,
std::size_t,
const std::array< std::array< tk::real, 3 >, 3 >& ) const
// *************************************************************************
//! \brief Calculate the Cauchy stress tensor from the material density,
//! momentum, and total energy
//! \return Material Cauchy stress tensor (alpha_k * sigma_k)
// *************************************************************************
{
std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};
// No elastic contribution
return asig;
}
tk::real
ThermallyPerfectGas::soundspeed(
tk::real rho,
tk::real pr,
tk::real,
std::size_t,
const std::array< std::array< tk::real, 3 >, 3 >&,
const std::array< tk::real, 3 >& ) const
// *************************************************************************
//! Calculate speed of sound from the material density and material pressure
//! \param[in] rho density
//! \param[in] pr pressure
//! \return Material speed of sound using the ideal gas EoS
// *************************************************************************
{
auto g = m_gamma;
tk::real a = std::sqrt( g * pr / rho );
return a;
}
tk::real
ThermallyPerfectGas::totalenergy(
tk::real rho,
tk::real u,
tk::real v,
tk::real w,
tk::real pr,
const std::array< std::array< tk::real, 3 >, 3 >& ) const
// *************************************************************************
//! \brief Calculate material specific total energy from the material
//! density, momentum and material pressure
//! \param[in] rho density
//! \param[in] u X-velocity
//! \param[in] v Y-velocity
//! \param[in] w Z-velocity
//! \param[in] pr pressure
//! \return specific total energy using the thermally perfect gas EoS
// *************************************************************************
{
auto R = m_R;
tk::real temp = pr / (rho * R);
tk::real e = R * (-m_cp_coeff[0] * std::pow(temp, -1) +
m_cp_coeff[1] * std::log(temp) + (m_cp_coeff[2] - 1) * temp +
m_cp_coeff[3] * std::pow(temp, 2) / 2 +
m_cp_coeff[4] * std::pow(temp, 3) / 3 +
m_cp_coeff[5] * std::pow(temp, 4) / 4 +
m_cp_coeff[6] * std::pow(temp, 5) / 5 + m_cp_coeff[7]);
return (rho * e + 0.5 * rho * (u*u + v*v + w*w));
}
tk::real
ThermallyPerfectGas::temperature(
tk::real rho,
tk::real u,
tk::real v,
tk::real w,
tk::real rhoE,
tk::real,
const std::array< std::array< tk::real, 3 >, 3 >& ) const
// *************************************************************************
//! \brief Calculate material temperature from the material density
//! \param[in] rho density
//! \param[in] u X-velocity
//! \param[in] v Y-velocity
//! \param[in] w Z-velocity
//! \param[in] rhoE total energy
//! \return Material temperature using the thermally perfect gas EoS
// *************************************************************************
{
auto R = m_R;
// Solve for internal energy
tk::real e = rhoE / rho - 0.5 * (u*u + v*v + w*w);
// Solve for temperature - Newton's method
tk::real temp = 1000; // Starting guess
tk::real tol = 1e-8; // Stopping condition
tk::real err = 1e8;<--- Variable 'err' is assigned a value that is never used.
std::size_t maxiter = 1000;
std::size_t i(0);
while (i < maxiter) {
tk::real f_T = R * (-m_cp_coeff[0] * std::pow(temp, -1) +
m_cp_coeff[1] * std::log(temp) + (m_cp_coeff[2] - 1) * temp +
m_cp_coeff[3] * std::pow(temp, 2) / 2 +
m_cp_coeff[4] * std::pow(temp, 3) / 3 +
m_cp_coeff[5] * std::pow(temp, 4) / 4 +
m_cp_coeff[6] * std::pow(temp, 5) / 5 + m_cp_coeff[7]) - e;
err = abs(f_T);
// Get derivative - df/dT. For loop is working through polynomial.
tk::real fp_T = 0;
tk::real power = -2;
for (std::size_t k=0; k<m_cp_coeff.size()-1; ++k)
{
fp_T += m_cp_coeff[k] * std::pow(temp, power);
if (k == 2) fp_T += -1;
power += 1;
}
fp_T = fp_T * R;
temp = temp - f_T / fp_T;
if (err <= tol) break;
i++;
if ( i == maxiter ) {
Throw("ThermallyPerfectGas Newton's Method for temperature failed to converge after iterations "
+ std::to_string(i));
}
}
return temp;
}
|