Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/EoS/StiffenedGas.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 Stiffened-gas equation of state
9 : : \details This file defines functions for the stiffened gas equation of
10 : : state for the compressible flow equations.
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include <cmath>
15 : : #include <iostream>
16 : : #include "EoS/StiffenedGas.hpp"
17 : :
18 : : using inciter::StiffenedGas;
19 : :
20 : 750 : StiffenedGas::StiffenedGas(
21 : : tk::real gamma,
22 : : tk::real pstiff,
23 : 750 : tk::real cv ) :
24 : : m_gamma(gamma),
25 : : m_pstiff(pstiff),
26 : 750 : m_cv(cv)
27 : : // *************************************************************************
28 : : // Constructor
29 : : //! \param[in] gamma Ratio of specific heats
30 : : //! \param[in] pstiff Stiffened pressure term
31 : : //! \param[in] cv Specific heat at constant volume
32 : : // *************************************************************************
33 : 750 : { }
34 : :
35 : : tk::real
36 : 2147340 : StiffenedGas::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 : 2147340 : tk::real g = m_gamma;
48 : 2147340 : tk::real p_c = m_pstiff;
49 : 2147340 : tk::real c_v = m_cv;
50 : :
51 : 2147340 : tk::real rho = (pr + p_c) / ((g-1.0) * c_v * temp);
52 : 2147340 : return rho;
53 : : }
54 : :
55 : : tk::real
56 : 150798511 : StiffenedGas::pressure(
57 : : tk::real arho,
58 : : tk::real u,
59 : : tk::real v,
60 : : tk::real w,
61 : : tk::real arhoE,
62 : : tk::real alpha,
63 : : std::size_t imat,
64 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
65 : : // *************************************************************************
66 : : //! \brief Calculate pressure from the material density, momentum and total
67 : : //! energy using the stiffened-gas equation of state
68 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
69 : : //! \param[in] u X-velocity
70 : : //! \param[in] v Y-velocity
71 : : //! \param[in] w Z-velocity
72 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
73 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
74 : : //! the single-material system, this argument can be left unspecified by
75 : : //! the calling code
76 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
77 : : //! for the single-material system, this argument can be left unspecified
78 : : //! by the calling code
79 : : //! \return Material partial pressure (alpha_k * p_k) calculated using the
80 : : //! stiffened-gas EoS
81 : : // *************************************************************************
82 : : {
83 : 150798511 : tk::real g = m_gamma;
84 : 150798511 : tk::real p_c = m_pstiff;
85 : :
86 : 150798511 : tk::real partpressure = (arhoE - 0.5 * arho * (u*u + v*v + w*w) -
87 : 150798511 : alpha*p_c) * (g-1.0) - alpha*p_c;
88 : :
89 : : // check partial pressure divergence
90 [ - + ]: 150798511 : if (!std::isfinite(partpressure)) {
91 : 0 : std::cout << "Material-id: " << imat << std::endl;
92 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
93 : 0 : std::cout << "Partial density: " << arho << std::endl;
94 : 0 : std::cout << "Total energy: " << arhoE << std::endl;
95 : 0 : std::cout << "Velocity: " << u << ", " << v << ", " << w
96 : 0 : << std::endl;
97 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) +
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
98 : : " has nan/inf partial pressure: " + std::to_string(partpressure) +
99 : : ", material volume fraction: " + std::to_string(alpha));
100 : : }
101 : :
102 : 150798511 : return partpressure;
103 : : }
104 : :
105 : : std::array< std::array< tk::real, 3 >, 3 >
106 : 0 : StiffenedGas::CauchyStress(
107 : : tk::real,
108 : : tk::real,
109 : : tk::real,
110 : : tk::real,
111 : : tk::real,
112 : : tk::real,
113 : : std::size_t,
114 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
115 : : // *************************************************************************
116 : : //! \brief Calculate the Cauchy stress tensor from the material density,
117 : : //! momentum, and total energy
118 : : //! \return Material Cauchy stress tensor (alpha_k * sigma_k)
119 : : // *************************************************************************
120 : : {
121 : 0 : std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};
122 : :
123 : : // No elastic contribution
124 : :
125 : 0 : return asig;
126 : : }
127 : :
128 : : tk::real
129 : 149977855 : StiffenedGas::soundspeed(
130 : : tk::real arho,
131 : : tk::real apr,
132 : : tk::real alpha,
133 : : std::size_t imat,
134 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
135 : : // *************************************************************************
136 : : //! Calculate speed of sound from the material density and material pressure
137 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
138 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
139 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
140 : : //! the single-material system, this argument can be left unspecified by
141 : : //! the calling code
142 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
143 : : //! for the single-material system, this argument can be left unspecified
144 : : //! by the calling code
145 : : //! \return Material speed of sound using the stiffened-gas EoS
146 : : // *************************************************************************
147 : : {
148 : 149977855 : auto g = m_gamma;
149 : 149977855 : auto p_c = m_pstiff;
150 : :
151 : 149977855 : auto p_eff = std::max( 1.0e-15, apr+(alpha*p_c) );
152 : :
153 : 149977855 : tk::real a = std::sqrt( g * p_eff / arho );
154 : :
155 : : // check sound speed divergence
156 [ - + ]: 149977855 : if (!std::isfinite(a)) {
157 : 0 : std::cout << "Material-id: " << imat << std::endl;
158 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
159 : 0 : std::cout << "Partial density: " << arho << std::endl;
160 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
161 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
162 : : + std::to_string(a) + ", material volume fraction: " +
163 : : std::to_string(alpha));
164 : : }
165 : :
166 : 149977855 : return a;
167 : : }
168 : :
169 : : tk::real
170 : 13809909 : StiffenedGas::totalenergy(
171 : : tk::real rho,
172 : : tk::real u,
173 : : tk::real v,
174 : : tk::real w,
175 : : tk::real pr,
176 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
177 : : // *************************************************************************
178 : : //! \brief Calculate material specific total energy from the material
179 : : //! density, momentum and material pressure
180 : : //! \param[in] rho Material density
181 : : //! \param[in] u X-velocity
182 : : //! \param[in] v Y-velocity
183 : : //! \param[in] w Z-velocity
184 : : //! \param[in] pr Material pressure
185 : : //! \return Material specific total energy using the stiffened-gas EoS
186 : : // *************************************************************************
187 : : {
188 : 13809909 : auto g = m_gamma;
189 : 13809909 : auto p_c = m_pstiff;
190 : :
191 : 13809909 : tk::real rhoE = (pr + p_c) / (g-1.0) + 0.5 * rho * (u*u + v*v + w*w) + p_c;
192 : 13809909 : return rhoE;
193 : : }
194 : :
195 : : tk::real
196 : 2365840 : StiffenedGas::temperature(
197 : : tk::real arho,
198 : : tk::real u,
199 : : tk::real v,
200 : : tk::real w,
201 : : tk::real arhoE,
202 : : tk::real alpha,
203 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
204 : : // *************************************************************************
205 : : //! \brief Calculate material temperature from the material density, and
206 : : //! material specific total energy
207 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
208 : : //! \param[in] u X-velocity
209 : : //! \param[in] v Y-velocity
210 : : //! \param[in] w Z-velocity
211 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
212 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
213 : : //! the single-material system, this argument can be left unspecified by
214 : : //! the calling code
215 : : //! \return Material temperature using the stiffened-gas EoS
216 : : // *************************************************************************
217 : : {
218 : 2365840 : auto c_v = m_cv;
219 : 2365840 : auto p_c = m_pstiff;
220 : :
221 : 2365840 : tk::real t = (arhoE - 0.5 * arho * (u*u + v*v + w*w) - alpha*p_c)
222 : 2365840 : / (arho*c_v);
223 : 2365840 : return t;
224 : : }
225 : :
226 : : tk::real
227 : 92086930 : StiffenedGas::min_eff_pressure(
228 : : tk::real min,
229 : : tk::real,
230 : : tk::real ) const
231 : : // *************************************************************************
232 : : //! Compute the minimum allowed pressure
233 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
234 : : //! \return Minimum pressure allowed by physical constraints
235 : : // *************************************************************************
236 : : {
237 : : // minimum pressure is constrained by zero soundspeed.
238 : 92086930 : return (min - m_pstiff);
239 : : }
|
