Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/EoS/JWL.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 Jones, Wilkins, and Lee (JWL) equation of state
9 : : \details This file defines functions for the JWL equation of
10 : : state for the compressible flow equations. These functions are
11 : : taken from 'JWL Equation of State', Menikoff, LA-UR-15-29536.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include <cmath>
16 : : #include <iostream>
17 : : #include "EoS/JWL.hpp"
18 : :
19 : : using inciter::JWL;
20 : :
21 : 0 : JWL::JWL( tk::real w, tk::real cv, tk::real rho0, tk::real de, tk::real rhor,
22 : 0 : tk::real tr, tk::real pr, tk::real A, tk::real B, tk::real R1, tk::real R2 ) :
23 : : m_w(w),
24 : : m_cv(cv),
25 : : m_rho0(rho0),
26 : : m_de(de),
27 : : m_rhor(rhor),
28 : : m_tr(tr),
29 : : m_pr(pr),
30 : : m_a(A),
31 : : m_b(B),
32 : : m_r1(R1),
33 : 0 : m_r2(R2)
34 : : // *************************************************************************
35 : : // Constructor
36 : : //! \param[in] w Grueneisen coefficient
37 : : //! \param[in] cv Specific heat at constant volume
38 : : //! \param[in] rho0 Density of initial state
39 : : //! \param[in] de Heat of detonation for products. For reactants, it is
40 : : //! chosen such that the ambient internal energy (e0) is 0.
41 : : //! \param[in] rhor Density of reference state
42 : : //! \param[in] tr Temperature of reference state
43 : : //! \param[in] pr Pressure of reference state
44 : : //! \param[in] A Parameter A
45 : : //! \param[in] B Parameter B
46 : : //! \param[in] R1 Parameter R1
47 : : //! \param[in] R2 Parameter R2
48 : : // *************************************************************************
49 : : {
50 : : // reference density provided
51 [ - - ]: 0 : if (m_tr < 1e-8) {
52 : : // reference internal energy
53 : 0 : auto er = intEnergy(rhor, pr);
54 : : // reference temperature from Eqn (15)
55 : 0 : m_tr = 1.0/m_cv * (er + de -
56 : 0 : (m_a/m_r1*exp(-m_r1*m_rho0/m_rhor) +
57 : 0 : m_b/m_r2*exp(-m_r2*m_rho0/m_rhor)) / m_rho0);
58 : : }
59 : : // reference temperature provided
60 : : else
61 : : {
62 : 0 : m_rhor = density(m_pr, m_tr);
63 : : }
64 : 0 : }
65 : :
66 : : tk::real
67 : 0 : JWL::density(
68 : : tk::real pr,
69 : : tk::real temp ) const
70 : : // *************************************************************************
71 : : //! \brief Calculate density from the material pressure and temperature
72 : : //! using the stiffened-gas equation of state
73 : : //! \param[in] pr Material pressure
74 : : //! \param[in] temp Material temperature
75 : : //! \return Material density calculated using the stiffened-gas EoS
76 : : // *************************************************************************
77 : : {
78 : 0 : tk::real r_guessL = 1e-4*m_rho0; // left density bound
79 : 0 : tk::real r_guessR = 1e2*m_rho0; // right density bound
80 : : tk::real rho;
81 : :
82 : 0 : rho = bisection( r_guessL, r_guessR, pr, temp );
83 : :
84 : 0 : return rho;
85 : : }
86 : :
87 : :
88 : : tk::real
89 : 0 : JWL::pressure(
90 : : tk::real arho,
91 : : tk::real u,
92 : : tk::real v,
93 : : tk::real w,
94 : : tk::real arhoE,
95 : : tk::real alpha,
96 : : std::size_t imat,
97 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
98 : : // *************************************************************************
99 : : //! \brief Calculate pressure from the material density, momentum and total
100 : : //! energy using the stiffened-gas equation of state
101 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
102 : : //! \param[in] u X-velocity
103 : : //! \param[in] v Y-velocity
104 : : //! \param[in] w Z-velocity
105 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
106 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
107 : : //! the single-material system, this argument can be left unspecified by
108 : : //! the calling code
109 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
110 : : //! for the single-material system, this argument can be left unspecified
111 : : //! by the calling code
112 : : //! \return Material partial pressure (alpha_k * p_k) calculated using the
113 : : //! stiffened-gas EoS
114 : : //! \details From Eqn. 1 in 'JWL Equation of State', Menikoff, LA-UR-15-29536
115 : : // *************************************************************************
116 : : {
117 : : // specific internal energy
118 : 0 : tk::real e = (arhoE - 0.5*arho*(u*u + v*v + w*w))/arho;
119 : :
120 : : //// reference energy (input quantity, might need for calculation)
121 : : //tk::real e0 = a/r1*exp(-r1*rho0/rho) + b/r2*exp(-r2*rho0/rho);
122 : :
123 : 0 : alpha = std::max(1e-14,alpha);
124 : : tk::real partpressure =
125 : 0 : m_a*(alpha - m_w*arho/(m_rho0*m_r1))*exp(-m_r1*alpha*m_rho0/arho) +
126 : 0 : m_b*(alpha - m_w*arho/(m_rho0*m_r2))*exp(-m_r2*alpha*m_rho0/arho) +
127 : 0 : m_w*arho*(e + m_de);
128 : :
129 : : // check partial pressure divergence
130 [ - - ]: 0 : if (!std::isfinite(partpressure)) {
131 : 0 : std::cout << "Material-id: " << imat << std::endl;
132 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
133 : 0 : std::cout << "Partial density: " << arho << std::endl;
134 : 0 : std::cout << "Total energy: " << arhoE << std::endl;
135 : 0 : std::cout << "Velocity: " << u << ", " << v << ", " << w
136 : 0 : << std::endl;
137 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) +
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
138 : : " has nan/inf partial pressure: " + std::to_string(partpressure) +
139 : : ", material volume fraction: " + std::to_string(alpha));
140 : : }
141 : :
142 : 0 : return partpressure;
143 : : }
144 : :
145 : : std::array< std::array< tk::real, 3 >, 3 >
146 : 0 : JWL::CauchyStress(
147 : : tk::real,
148 : : tk::real,
149 : : tk::real,
150 : : tk::real,
151 : : tk::real,
152 : : tk::real,
153 : : std::size_t,
154 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
155 : : // *************************************************************************
156 : : //! \brief Calculate the Cauchy stress tensor from the material density,
157 : : //! momentum, and total energy
158 : : //! \return Material Cauchy stress tensor (alpha_k * sigma_k)
159 : : // *************************************************************************
160 : : {
161 : 0 : std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};
162 : :
163 : : // No elastic contribution
164 : :
165 : 0 : return asig;
166 : : }
167 : :
168 : : tk::real
169 : 0 : JWL::soundspeed(
170 : : tk::real arho,
171 : : tk::real apr,
172 : : tk::real alpha,
173 : : std::size_t imat,
174 : : const std::array< std::array< tk::real, 3 >, 3 >&,
175 : : const std::array< tk::real, 3 >& ) const
176 : : // *************************************************************************
177 : : //! Calculate speed of sound from the material density and material pressure
178 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
179 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
180 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
181 : : //! the single-material system, this argument can be left unspecified by
182 : : //! the calling code
183 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
184 : : //! for the single-material system, this argument can be left unspecified
185 : : //! by the calling code
186 : : //! \return Material speed of sound using the stiffened-gas EoS
187 : : // *************************************************************************
188 : : {
189 : 0 : alpha = std::max(1e-14,alpha);
190 : : // limiting pressure to near-zero
191 : 0 : auto apr_eff = std::max(alpha*
192 : 0 : min_eff_pressure(1e-4*std::abs(apr/alpha), arho, alpha), apr);
193 : :
194 : 0 : auto co1 = m_rho0*alpha*alpha/(arho*arho);
195 : 0 : auto co2 = alpha*(1.0+m_w)/arho;
196 : :
197 : 0 : tk::real ss = m_a*(m_r1*co1 - co2) * exp(-m_r1*alpha*m_rho0/arho)
198 : 0 : + m_b*(m_r2*co1 - co2) * exp(-m_r2*alpha*m_rho0/arho)
199 : 0 : + (1.0+m_w)*apr_eff/arho;
200 : :
201 : 0 : auto ss2 = ss;
202 : 0 : ss = std::sqrt(ss);
203 : :
204 : : // check sound speed divergence
205 [ - - ]: 0 : if (!std::isfinite(ss)) {
206 : 0 : std::cout << "Material-id: " << imat << std::endl;
207 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
208 : 0 : std::cout << "Partial density: " << arho << std::endl;
209 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
210 : 0 : std::cout << "Min allowed pres: " << alpha*min_eff_pressure(0.0, arho,
211 : 0 : alpha) << std::endl;
212 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed. "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
213 : : "ss^2: " + std::to_string(ss2) + ", material volume fraction: " +
214 : : std::to_string(alpha));
215 : : }
216 : :
217 : 0 : return ss;
218 : : }
219 : :
220 : : tk::real
221 : 0 : JWL::totalenergy(
222 : : tk::real rho,
223 : : tk::real u,
224 : : tk::real v,
225 : : tk::real w,
226 : : tk::real pr,
227 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
228 : : // *************************************************************************
229 : : //! \brief Calculate material specific total energy from the material
230 : : //! density, momentum and material pressure
231 : : //! \param[in] rho Material density
232 : : //! \param[in] u X-velocity
233 : : //! \param[in] v Y-velocity
234 : : //! \param[in] w Z-velocity
235 : : //! \param[in] pr Material pressure
236 : : //! \return Material specific total energy using the stiffened-gas EoS
237 : : // *************************************************************************
238 : : {
239 : : //// reference energy (input quantity, might need for calculation)
240 : : //tk::real e0 = a/r1*exp(-r1*rho0/rho) + b/r2*exp(-r2*rho0/rho);
241 : :
242 : 0 : tk::real rhoE = rho*intEnergy( rho, pr )
243 : 0 : + 0.5*rho*(u*u + v*v + w*w);
244 : :
245 : 0 : return rhoE;
246 : : }
247 : :
248 : : tk::real
249 : 0 : JWL::temperature(
250 : : tk::real arho,
251 : : tk::real u,
252 : : tk::real v,
253 : : tk::real w,
254 : : tk::real arhoE,
255 : : tk::real alpha,
256 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
257 : : // *************************************************************************
258 : : //! \brief Calculate material temperature from the material density, and
259 : : //! material specific total energy
260 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
261 : : //! \param[in] u X-velocity
262 : : //! \param[in] v Y-velocity
263 : : //! \param[in] w Z-velocity
264 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
265 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
266 : : //! the single-material system, this argument can be left unspecified by
267 : : //! the calling code
268 : : //! \return Material temperature using the stiffened-gas EoS
269 : : // *************************************************************************
270 : : {
271 : 0 : alpha = std::max(1e-14,alpha);
272 : 0 : tk::real rho = arho/alpha;
273 : :
274 : : //// reference energy (input quantity, might need for calculation)
275 : : //tk::real e0 = a/r1*exp(-r1*rho0/rho) + b/r2*exp(-r2*rho0/rho);
276 : :
277 : 0 : tk::real t = ((arhoE - 0.5*arho*(u*u + v*v + w*w))/arho + m_de -
278 : 0 : 1.0/m_rho0*( m_a/m_r1*exp(-m_r1*m_rho0/rho)
279 : 0 : + m_b/m_r2*exp(-m_r2*m_rho0/rho) ))/m_cv;
280 : :
281 : 0 : return t;
282 : : }
283 : :
284 : : tk::real
285 : 0 : JWL::min_eff_pressure(
286 : : tk::real min,
287 : : tk::real arho,
288 : : tk::real alpha ) const
289 : : // *************************************************************************
290 : : //! Compute the minimum allowed pressure
291 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
292 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
293 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
294 : : //! the single-material system, this argument can be left unspecified by
295 : : //! the calling code
296 : : //! \return Minimum pressure allowed by physical constraints
297 : : // *************************************************************************
298 : : {
299 : 0 : alpha = std::max(1e-14,alpha);
300 : 0 : auto co1 = m_rho0*alpha*alpha/(arho*arho);
301 : 0 : auto co2 = alpha*(1.0+m_w)/arho;
302 : :
303 : : // minimum pressure is constrained by zero soundspeed.
304 : 0 : auto apr = -(m_a*(m_r1*co1 - co2) * exp(-m_r1*alpha*m_rho0/arho)
305 : 0 : + m_b*(m_r2*co1 - co2) * exp(-m_r2*alpha*m_rho0/arho))
306 : 0 : * arho/(1.0+m_w);
307 : :
308 : 0 : return apr/alpha + min;
309 : : }
310 : :
311 : : tk::real
312 : 0 : JWL::intEnergy(
313 : : tk::real rho,
314 : : tk::real pr ) const
315 : : // *************************************************************************
316 : : //! \brief Calculate specific internal energy using the JWL equation of
317 : : //! state
318 : : //! \param[in] rho Material density
319 : : //! \param[in] pr Material pressure
320 : : //! \return Material internal energy calculated using the JWL EoS
321 : : //! \details By inverting Eqn. 1 in 'JWL Equation of State', Menikoff,
322 : : //! LA-UR-15-29536
323 : : // *************************************************************************
324 : : {
325 : 0 : tk::real e = - m_de + 1.0/m_w/rho*( pr
326 : 0 : - m_a*(1.0 - m_w*rho/m_r1/m_rho0)*exp(-m_r1*m_rho0/rho)
327 : 0 : - m_b*(1.0 - m_w*rho/m_r2/m_rho0)*exp(-m_r2*m_rho0/rho) );
328 : :
329 : 0 : return e;
330 : : }
331 : :
332 : : tk::real
333 : 0 : JWL::bisection(
334 : : tk::real a,
335 : : tk::real b,
336 : : tk::real p_known,
337 : : tk::real t_known ) const
338 : : // *************************************************************************
339 : : //! \brief Calculate density from known pressure and temperature using
340 : : //! bisection root finding method for JWL equation of state
341 : : //! \param[in] a Left density bound for root finding
342 : : //! \param[in] b Right density bound for root finding
343 : : //! \param[in] p_known Known pressure
344 : : //! \param[in] t_known Known temperature
345 : : //! \return Material density calculated by inverting JWL pressure equation
346 : : // *************************************************************************
347 : : {
348 : 0 : tk::real tol = 1e-12;
349 : 0 : std::size_t maxiter = 1000;
350 : 0 : std::size_t i(0);
351 : : tk::real c;
352 : 0 : tk::real root(0);
353 : 0 : std::size_t idebug = 0;
354 : :
355 : : // Ensure that original bounds contain root
356 [ - - ]: 0 : while ( p_known < PfromRT( a, t_known)) {
357 : 0 : b = a;
358 : 0 : a *= 1e-6;
359 : : }
360 [ - - ]: 0 : while ( p_known > PfromRT( b, t_known)) {
361 : 0 : a = b;
362 : 0 : b *= 1e6;
363 : : }
364 : :
365 : : // function to minimize: fcn = p_known - PfromRT
366 : : // bounds b > a
367 : :
368 [ - - ]: 0 : while (i < maxiter)
369 : : {
370 : 0 : c = (a + b)/2.0;
371 : 0 : auto fcn = p_known - PfromRT( c, t_known);
372 [ - - ]: 0 : if ( idebug == 1)
373 : : {
374 : 0 : std::cout << "Bisection iter: " << i << std::endl;
375 : 0 : std::cout << "fcn: " << fcn << std::endl;
376 : 0 : std::cout << "(b - a)/2.0: " << (b - a)/2.0 << std::endl;
377 : : }
378 : :
379 [ - - ][ - - ]: 0 : if ( std::abs(fcn) <= 1e-16 or (b - a)/2.0 < tol )
[ - - ]
380 : : {
381 : 0 : root = c;
382 : 0 : break;
383 : : }
384 : :
385 : 0 : i++;
386 : 0 : if ( static_cast< int > (std::copysign( 1.0, p_known - PfromRT( c, t_known) )) ==
387 [ - - ]: 0 : static_cast< int > (std::copysign( 1.0, p_known - PfromRT( a, t_known) )) )
388 : : {
389 : 0 : a = c;
390 : : }
391 : : else
392 : : {
393 : 0 : b = c;
394 : : }
395 : :
396 [ - - ]: 0 : if ( i == maxiter )
397 : : {
398 [ - - ][ - - ]: 0 : Throw("JWL Bisection for density failed to converge after iterations "
[ - - ][ - - ]
399 : : + std::to_string(i));
400 : : }
401 : :
402 : : }
403 : 0 : return root;
404 : : }
405 : :
406 : :
407 : : tk::real
408 : 0 : JWL::PfromRT(
409 : : tk::real rho,
410 : : tk::real T ) const
411 : : // *************************************************************************
412 : : //! \brief Calculate pressure from density and temperature using JWL
413 : : //! equation of state
414 : : //! \param[in] rho Material density
415 : : //! \param[in] T Material temperature
416 : : //! \return Material pressure calculated using the JWL EoS
417 : : //! \details From Eqn. 14 in 'JWL Equation of State', Menikoff, LA-UR-15-29536
418 : : // *************************************************************************
419 : : {
420 : 0 : return ( m_a*exp(-m_r1*m_rho0/rho) + m_b*exp(-m_r2*m_rho0/rho) +
421 : 0 : m_w*(m_cv*T*rho) );
422 : : }
|
