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 : 1052 : StiffenedGas::StiffenedGas(
21 : : tk::real gamma,
22 : : tk::real pstiff,
23 : 1052 : tk::real cv ) :
24 : : m_gamma(gamma),
25 : : m_pstiff(pstiff),
26 : 1052 : 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 : 1052 : { }
34 : :
35 : : tk::real
36 : 2233602 : 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 : 2233602 : tk::real g = m_gamma;
48 : 2233602 : tk::real p_c = m_pstiff;
49 : 2233602 : tk::real c_v = m_cv;
50 : :
51 : 2233602 : tk::real rho = (pr + p_c) / ((g-1.0) * c_v * temp);
52 : 2233602 : return rho;
53 : : }
54 : :
55 : : tk::real
56 : 158862881 : 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 : 158862881 : tk::real g = m_gamma;
84 : 158862881 : tk::real p_c = m_pstiff;
85 : :
86 : 158862881 : tk::real partpressure = (arhoE - 0.5 * arho * (u*u + v*v + w*w)) * (g-1.0) -
87 : 158862881 : alpha*g*p_c;
88 : :
89 : : // check partial pressure divergence
90 [ - + ]: 158862881 : 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 : 158862881 : return partpressure;
103 : : }
104 : :
105 : : std::array< std::array< tk::real, 3 >, 3 >
106 : 0 : StiffenedGas::CauchyStress(
107 : : tk::real,
108 : : std::size_t,
109 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
110 : : // *************************************************************************
111 : : //! \brief Calculate the Cauchy stress tensor from the material
112 : : //! inverse deformation gradient tensor
113 : : //! \return Material Cauchy stress tensor (alpha_k * sigma_k)
114 : : // *************************************************************************
115 : : {
116 : 0 : std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};
117 : :
118 : : // No elastic contribution
119 : :
120 : 0 : return asig;
121 : : }
122 : :
123 : : tk::real
124 : 160174129 : StiffenedGas::soundspeed(
125 : : tk::real arho,
126 : : tk::real apr,
127 : : tk::real alpha,
128 : : std::size_t imat,
129 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
130 : : // *************************************************************************
131 : : //! Calculate speed of sound from the material density and material pressure
132 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
133 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
134 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
135 : : //! the single-material system, this argument can be left unspecified by
136 : : //! the calling code
137 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
138 : : //! for the single-material system, this argument can be left unspecified
139 : : //! by the calling code
140 : : //! \return Material speed of sound using the stiffened-gas EoS
141 : : // *************************************************************************
142 : : {
143 : 160174129 : auto g = m_gamma;
144 : 160174129 : auto p_c = m_pstiff;
145 : :
146 : 160174129 : auto al_eff = std::max( 1.0e-14, alpha );
147 : 160174129 : auto p_eff = std::max( 1.0e-15, apr+(al_eff*p_c) );
148 : :
149 : 160174129 : tk::real a = std::sqrt( g * p_eff / arho );
150 : :
151 : : // check sound speed divergence
152 [ - + ]: 160174129 : if (!std::isfinite(a)) {
153 : 0 : std::cout << "Material-id: " << imat << std::endl;
154 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
155 : 0 : std::cout << "Partial density: " << arho << std::endl;
156 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
157 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
158 : : + std::to_string(a) + ", material volume fraction: " +
159 : : std::to_string(alpha));
160 : : }
161 : :
162 : 160174129 : return a;
163 : : }
164 : :
165 : : tk::real
166 : 14773393 : StiffenedGas::totalenergy(
167 : : tk::real arho,
168 : : tk::real u,
169 : : tk::real v,
170 : : tk::real w,
171 : : tk::real apr,
172 : : tk::real alpha,
173 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
174 : : // *************************************************************************
175 : : //! \brief Calculate material specific total energy from the material
176 : : //! density, momentum and material pressure
177 : : //! \param[in] arho Material partial density
178 : : //! \param[in] u X-velocity
179 : : //! \param[in] v Y-velocity
180 : : //! \param[in] w Z-velocity
181 : : //! \param[in] apr Material partial pressure
182 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
183 : : //! the single-material system, this argument can be left unspecified by
184 : : //! the calling code
185 : : //! \return Material specific total energy using the stiffened-gas EoS
186 : : // *************************************************************************
187 : : {
188 : 14773393 : auto g = m_gamma;
189 : 14773393 : auto p_c = m_pstiff;
190 : :
191 : 14773393 : tk::real arhoE = (apr + alpha*g*p_c) / (g-1.0) + 0.5 * arho * (u*u + v*v + w*w);
192 : 14773393 : return arhoE;
193 : : }
194 : :
195 : : tk::real
196 : 2762480 : 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 : 2762480 : auto c_v = m_cv;
219 : 2762480 : auto p_c = m_pstiff;
220 : :
221 : 2762480 : tk::real t = (arhoE - 0.5 * arho * (u*u + v*v + w*w) - alpha*p_c)
222 : 2762480 : / (arho*c_v);
223 : 2762480 : return t;
224 : : }
225 : :
226 : : tk::real
227 : 98419266 : 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 : 98419266 : return (min - m_pstiff);
239 : : }
|
