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