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 : 788 : StiffenedGas::StiffenedGas(
21 : : tk::real gamma,
22 : : tk::real pstiff,
23 : 788 : tk::real cv ) :
24 : : m_gamma(gamma),
25 : : m_pstiff(pstiff),
26 : 788 : 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 : 788 : { }
34 : :
35 : : tk::real
36 : 2176098 : 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 : 2176098 : tk::real g = m_gamma;
48 : 2176098 : tk::real p_c = m_pstiff;
49 : 2176098 : tk::real c_v = m_cv;
50 : :
51 : 2176098 : tk::real rho = (pr + p_c) / ((g-1.0) * c_v * temp);
52 : 2176098 : return rho;
53 : : }
54 : :
55 : : tk::real
56 : 172473923 : 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 : 172473923 : tk::real g = m_gamma;
84 : 172473923 : tk::real p_c = m_pstiff;
85 : :
86 : 172473923 : tk::real partpressure = (arhoE - 0.5 * arho * (u*u + v*v + w*w) -
87 : 172473923 : alpha*p_c) * (g-1.0) - alpha*p_c;
88 : :
89 : : // check partial pressure divergence
90 [ - + ]: 172473923 : 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 : 172473923 : 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 : 172633893 : 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 >&,
135 : : const std::array< tk::real, 3 >&,
136 : : const std::array< tk::real, 3 >& ) const
137 : : // *************************************************************************
138 : : //! Calculate speed of sound from the material density and material pressure
139 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
140 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
141 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
142 : : //! the single-material system, this argument can be left unspecified by
143 : : //! the calling code
144 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
145 : : //! for the single-material system, this argument can be left unspecified
146 : : //! by the calling code
147 : : //! \return Material speed of sound using the stiffened-gas EoS
148 : : // *************************************************************************
149 : : {
150 : 172633893 : auto g = m_gamma;
151 : 172633893 : auto p_c = m_pstiff;
152 : :
153 : 172633893 : auto p_eff = std::max( 1.0e-15, apr+(alpha*p_c) );
154 : :
155 : 172633893 : tk::real a = std::sqrt( g * p_eff / arho );
156 : :
157 : : // check sound speed divergence
158 [ - + ]: 172633893 : if (!std::isfinite(a)) {
159 : 0 : std::cout << "Material-id: " << imat << std::endl;
160 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
161 : 0 : std::cout << "Partial density: " << arho << std::endl;
162 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
163 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
164 : : + std::to_string(a) + ", material volume fraction: " +
165 : : std::to_string(alpha));
166 : : }
167 : :
168 : 172633893 : return a;
169 : : }
170 : :
171 : : tk::real
172 : 14783292 : StiffenedGas::totalenergy(
173 : : tk::real rho,
174 : : tk::real u,
175 : : tk::real v,
176 : : tk::real w,
177 : : tk::real pr,
178 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
179 : : // *************************************************************************
180 : : //! \brief Calculate material specific total energy from the material
181 : : //! density, momentum and material pressure
182 : : //! \param[in] rho Material density
183 : : //! \param[in] u X-velocity
184 : : //! \param[in] v Y-velocity
185 : : //! \param[in] w Z-velocity
186 : : //! \param[in] pr Material pressure
187 : : //! \return Material specific total energy using the stiffened-gas EoS
188 : : // *************************************************************************
189 : : {
190 : 14783292 : auto g = m_gamma;
191 : 14783292 : auto p_c = m_pstiff;
192 : :
193 : 14783292 : tk::real rhoE = (pr + p_c) / (g-1.0) + 0.5 * rho * (u*u + v*v + w*w) + p_c;
194 : 14783292 : return rhoE;
195 : : }
196 : :
197 : : tk::real
198 : 2523784 : StiffenedGas::temperature(
199 : : tk::real arho,
200 : : tk::real u,
201 : : tk::real v,
202 : : tk::real w,
203 : : tk::real arhoE,
204 : : tk::real alpha,
205 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
206 : : // *************************************************************************
207 : : //! \brief Calculate material temperature from the material density, and
208 : : //! material specific total energy
209 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
210 : : //! \param[in] u X-velocity
211 : : //! \param[in] v Y-velocity
212 : : //! \param[in] w Z-velocity
213 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
214 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
215 : : //! the single-material system, this argument can be left unspecified by
216 : : //! the calling code
217 : : //! \return Material temperature using the stiffened-gas EoS
218 : : // *************************************************************************
219 : : {
220 : 2523784 : auto c_v = m_cv;
221 : 2523784 : auto p_c = m_pstiff;
222 : :
223 : 2523784 : tk::real t = (arhoE - 0.5 * arho * (u*u + v*v + w*w) - alpha*p_c)
224 : 2523784 : / (arho*c_v);
225 : 2523784 : return t;
226 : : }
227 : :
228 : : tk::real
229 : 103692166 : StiffenedGas::min_eff_pressure(
230 : : tk::real min,
231 : : tk::real,
232 : : tk::real ) const
233 : : // *************************************************************************
234 : : //! Compute the minimum allowed pressure
235 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
236 : : //! \return Minimum pressure allowed by physical constraints
237 : : // *************************************************************************
238 : : {
239 : : // minimum pressure is constrained by zero soundspeed.
240 : 103692166 : return (min - m_pstiff);
241 : : }
|
