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 : 36 : ThermallyPerfectGas::ThermallyPerfectGas(
21 : : tk::real gamma,
22 : : tk::real R,
23 : : std::vector< std::vector< tk::real > > cp_coeff,
24 : : std::vector< tk::real > t_range,
25 : 36 : tk::real dH_ref) :
26 : : m_gamma(gamma),
27 : : m_R(R),
28 : : m_cp_coeff(cp_coeff),
29 : : m_t_range(t_range),
30 [ + - ]: 36 : m_dH_ref(dH_ref)
31 : : // *************************************************************************
32 : : // Constructor
33 : : //! \param[in] gamma Ratio of specific heats
34 : : //! \param[in] R gas constant
35 : : //! \param[in] cp_coeff NASA Glenn polynomials coefficients for cp fit
36 : : //! \param[in] t_range temperature range where polynomial coeffs are valid
37 : : //! \param[in] dH_ref reference enthalpy, h(t = 298.15 K) - h(t = 0 K)
38 : : // *************************************************************************
39 : 36 : { }
40 : :
41 : : tk::real
42 : 363816 : ThermallyPerfectGas::density(
43 : : tk::real pr,
44 : : tk::real temp ) const
45 : : // *************************************************************************
46 : : //! \brief Calculate density from the material pressure and temperature
47 : : //! using the stiffened-gas equation of state
48 : : //! \param[in] pr Material pressure
49 : : //! \param[in] temp Material temperature
50 : : //! \return Material density calculated using the stiffened-gas EoS
51 : : // *************************************************************************
52 : : {
53 : 363816 : tk::real R = m_R;
54 : :
55 : 363816 : tk::real rho = pr / (R * temp);
56 : 363816 : return rho;
57 : : }
58 : :
59 : : tk::real
60 : 3504600 : ThermallyPerfectGas::pressure(
61 : : tk::real rho,
62 : : tk::real u,
63 : : tk::real v,
64 : : tk::real w,
65 : : tk::real rhoE,
66 : : tk::real,
67 : : std::size_t,
68 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
69 : : // *************************************************************************
70 : : //! \brief Calculate pressure from the material density, momentum and total
71 : : //! energy using the thermally perfect gas equation of state
72 : : //! \param[in] rho density
73 : : //! \param[in] u X-velocity
74 : : //! \param[in] v Y-velocity
75 : : //! \param[in] w Z-velocity
76 : : //! \param[in] rhoE total energy
77 : : //! \return Pressure calculated using the thermally perfect gas EOS
78 : : // *************************************************************************
79 : : {
80 : 3504600 : tk::real R = m_R;
81 : :
82 [ + - ]: 3504600 : tk::real temp = temperature(rho, u, v, w, rhoE);
83 : 3504600 : tk::real pres = rho * R * temp;
84 : :
85 : 3504600 : return pres;
86 : : }
87 : :
88 : : std::array< std::array< tk::real, 3 >, 3 >
89 : 0 : ThermallyPerfectGas::CauchyStress(
90 : : tk::real,
91 : : tk::real,
92 : : tk::real,
93 : : tk::real,
94 : : tk::real,
95 : : tk::real,
96 : : std::size_t,
97 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
98 : : // *************************************************************************
99 : : //! \brief Calculate the Cauchy stress tensor from the material density,
100 : : //! momentum, and total energy
101 : : //! \return Material Cauchy stress tensor (alpha_k * sigma_k)
102 : : // *************************************************************************
103 : : {
104 : 0 : std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};
105 : :
106 : : // No elastic contribution
107 : :
108 : 0 : return asig;
109 : : }
110 : :
111 : : tk::real
112 : 3504600 : ThermallyPerfectGas::soundspeed(
113 : : tk::real rho,
114 : : tk::real pr,
115 : : tk::real,
116 : : std::size_t,
117 : : const std::array< std::array< tk::real, 3 >, 3 >&,
118 : : const std::array< tk::real, 3 >& ) const
119 : : // *************************************************************************
120 : : //! Calculate speed of sound from the material density and material pressure
121 : : //! \param[in] rho density
122 : : //! \param[in] pr pressure
123 : : //! \return Material speed of sound using the ideal gas EoS
124 : : // *************************************************************************
125 : : {
126 : 3504600 : auto g = m_gamma;
127 : :
128 : 3504600 : tk::real a = std::sqrt( g * pr / rho );
129 : :
130 : 3504600 : return a;
131 : : }
132 : :
133 : : tk::real
134 : 243294 : ThermallyPerfectGas::totalenergy(
135 : : tk::real rho,
136 : : tk::real u,
137 : : tk::real v,
138 : : tk::real w,
139 : : tk::real pr,
140 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
141 : : // *************************************************************************
142 : : //! \brief Calculate material specific total energy from the material
143 : : //! density, momentum and material pressure
144 : : //! \param[in] rho density
145 : : //! \param[in] u X-velocity
146 : : //! \param[in] v Y-velocity
147 : : //! \param[in] w Z-velocity
148 : : //! \param[in] pr pressure
149 : : //! \return specific total energy using the thermally perfect gas EoS
150 : : // *************************************************************************
151 : : {
152 : 243294 : auto R = m_R;
153 : :
154 : 243294 : tk::real temp = pr / (rho * R);
155 : : // Identify what temperature range this falls in
156 : 243294 : std::size_t t_rng_idx = get_t_range(temp);
157 : :
158 : : // h = h_poly(T) + h_ref = e + R T (perfect gas)
159 : 243294 : tk::real e = R * (-m_cp_coeff[t_rng_idx][0] * std::pow(temp, -1) +
160 : 243294 : m_cp_coeff[t_rng_idx][1] * std::log(temp) + (m_cp_coeff[t_rng_idx][2] - 1) * temp +
161 : 243294 : m_cp_coeff[t_rng_idx][3] * std::pow(temp, 2) / 2 +
162 : 243294 : m_cp_coeff[t_rng_idx][4] * std::pow(temp, 3) / 3 +
163 : 243294 : m_cp_coeff[t_rng_idx][5] * std::pow(temp, 4) / 4 +
164 : 243294 : m_cp_coeff[t_rng_idx][6] * std::pow(temp, 5) / 5 + m_cp_coeff[t_rng_idx][7]) +
165 : 243294 : m_dH_ref;
166 : :
167 : 243294 : return (rho * e + 0.5 * rho * (u*u + v*v + w*w));
168 : : }
169 : :
170 : : tk::real
171 : 3504600 : ThermallyPerfectGas::temperature(
172 : : tk::real rho,
173 : : tk::real u,
174 : : tk::real v,
175 : : tk::real w,
176 : : tk::real rhoE,
177 : : tk::real,
178 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
179 : : // *************************************************************************
180 : : //! \brief Calculate material temperature from the material density
181 : : //! \param[in] rho density
182 : : //! \param[in] u X-velocity
183 : : //! \param[in] v Y-velocity
184 : : //! \param[in] w Z-velocity
185 : : //! \param[in] rhoE total energy
186 : : //! \return Material temperature using the thermally perfect gas EoS
187 : : // *************************************************************************
188 : : {
189 : 3504600 : auto R = m_R;
190 : :
191 : : // Solve for internal energy
192 : 3504600 : tk::real e = rhoE / rho - 0.5 * (u*u + v*v + w*w);
193 : :
194 : : // Solve for temperature - Newton's method
195 : 3504600 : tk::real temp = 1500; // Starting guess
196 : 3504600 : tk::real tol = 1e-8 * e; // Stopping condition
197 : : tk::real err;
198 : 3504600 : std::size_t maxiter = 10;
199 : 3504600 : std::size_t i(0);
200 [ + - ]: 11371363 : while (i < maxiter) {
201 : : // Identify what temperature range the current guess is in
202 : 11371363 : std::size_t t_rng_idx = get_t_range(temp);
203 : :
204 : : // With correct polynomial coefficients, construct e(temp) and de(temp)/dT
205 : 11371363 : tk::real f_T = R * (-m_cp_coeff[t_rng_idx][0] * std::pow(temp, -1) +
206 : 11371363 : m_cp_coeff[t_rng_idx][1] * std::log(temp) + (m_cp_coeff[t_rng_idx][2] - 1) * temp +
207 : 11371363 : m_cp_coeff[t_rng_idx][3] * std::pow(temp, 2) / 2 +
208 : 11371363 : m_cp_coeff[t_rng_idx][4] * std::pow(temp, 3) / 3 +
209 : 11371363 : m_cp_coeff[t_rng_idx][5] * std::pow(temp, 4) / 4 +
210 : 11371363 : m_cp_coeff[t_rng_idx][6] * std::pow(temp, 5) / 5 + m_cp_coeff[t_rng_idx][7]) +
211 : 11371363 : m_dH_ref - e;
212 : :
213 : 11371363 : err = abs(f_T);
214 : :
215 : : // Get derivative - df/dT. For loop is working through polynomial.
216 : 11371363 : tk::real fp_T = 0;
217 : 11371363 : tk::real power = -2;
218 [ + + ]: 90970904 : for (std::size_t k=0; k<m_cp_coeff[t_rng_idx].size()-1; ++k)
219 : : {
220 : 79599541 : fp_T += m_cp_coeff[t_rng_idx][k] * std::pow(temp, power);
221 [ + + ]: 79599541 : if (k == 2) fp_T += -1;
222 : 79599541 : power += 1;
223 : : }
224 : 11371363 : fp_T = fp_T * R;
225 : :
226 : : // Calculate next guess
227 : 11371363 : temp = temp - f_T / fp_T;
228 : :
229 [ + + ]: 11371363 : if (err <= tol) break;
230 : 7866763 : i++;
231 [ - + ]: 7866763 : if ( i == maxiter ) {
232 [ - - ][ - - ]: 0 : Throw("ThermallyPerfectGas Newton's Method for temperature failed to converge after iterations "
[ - - ][ - - ]
233 : : + std::to_string(i));
234 : : }
235 : : }
236 : :
237 : 3504600 : return temp;
238 : : }
|
