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 : auto p_eff = std::max( 1.0e-15, pr );
129 : 3504600 : tk::real a = std::sqrt( g * p_eff / rho );
130 : :
131 : 3504600 : return a;
132 : : }
133 : :
134 : : tk::real
135 : 243294 : ThermallyPerfectGas::totalenergy(
136 : : tk::real rho,
137 : : tk::real u,
138 : : tk::real v,
139 : : tk::real w,
140 : : tk::real pr,
141 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
142 : : // *************************************************************************
143 : : //! \brief Calculate material specific total energy from the material
144 : : //! density, momentum and material pressure
145 : : //! \param[in] rho density
146 : : //! \param[in] u X-velocity
147 : : //! \param[in] v Y-velocity
148 : : //! \param[in] w Z-velocity
149 : : //! \param[in] pr pressure
150 : : //! \return specific total energy using the thermally perfect gas EoS
151 : : // *************************************************************************
152 : : {
153 : 243294 : auto R = m_R;
154 : :
155 : 243294 : tk::real temp = pr / (rho * R);
156 : :
157 : : // h = h_poly(T) + h_ref = e + R T (perfect gas)
158 : 243294 : auto h = calc_h(temp) * R * temp + m_dH_ref; // dimensionalize the results out of the calc_h func
159 : 243294 : tk::real e = h - R * temp;
160 : :
161 : 243294 : return (rho * e + 0.5 * rho * (u*u + v*v + w*w));
162 : : }
163 : :
164 : : tk::real
165 : 3504600 : ThermallyPerfectGas::temperature(
166 : : tk::real rho,
167 : : tk::real u,
168 : : tk::real v,
169 : : tk::real w,
170 : : tk::real rhoE,
171 : : tk::real,
172 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
173 : : // *************************************************************************
174 : : //! \brief Calculate material temperature from the material density
175 : : //! \param[in] rho density
176 : : //! \param[in] u X-velocity
177 : : //! \param[in] v Y-velocity
178 : : //! \param[in] w Z-velocity
179 : : //! \param[in] rhoE total energy
180 : : //! \return Material temperature using the thermally perfect gas EoS
181 : : // *************************************************************************
182 : : {
183 : 3504600 : auto R = m_R;
184 : :
185 : : // Solve for internal energy
186 : 3504600 : tk::real e = rhoE / rho - 0.5 * (u*u + v*v + w*w);
187 [ - + ]: 3504600 : if (e < 1e-8) e = 1e-8; // TODO: standin until positivity is implemented
188 : :
189 : : // Solve for temperature - Newton's method
190 : 3504600 : tk::real temp = 1500; // Starting guess
191 : 3504600 : tk::real tol = std::max(1e-8, 1e-8 * e); // Stopping condition
192 : : tk::real err;
193 : 3504600 : std::size_t maxiter = 10;
194 : 3504600 : std::size_t i(0);
195 [ + - ]: 11371363 : while (i < maxiter) {
196 : : // Construct e(temp) and de(temp)/dT
197 : 11371363 : tk::real h = calc_h(temp) * R * temp + m_dH_ref;
198 : 11371363 : tk::real f_T = h - R * temp - e;
199 : :
200 : 11371363 : err = abs(f_T);
201 : :
202 : : // Get derivative - df/dT.
203 : 11371363 : tk::real cp = calc_cp(temp) * R;
204 : 11371363 : tk::real fp_T = cp - R;
205 : :
206 : : // Calculate next guess
207 : 11371363 : temp = temp - f_T / fp_T;
208 : :
209 [ + + ]: 11371363 : if (err <= tol) break;
210 : 7866763 : i++;
211 [ - + ]: 7866763 : if ( i == maxiter ) {
212 [ - - ][ - - ]: 0 : Throw("ThermallyPerfectGas Newton's Method for temperature failed to converge after iterations "
[ - - ][ - - ]
213 : : + std::to_string(i));
214 : : }
215 : : }
216 : :
217 : 3504600 : return temp;
218 : : }
|
