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