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 >& ) const
175 : : // *************************************************************************
176 : : //! Calculate speed of sound from the material density and material pressure
177 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
178 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
179 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
180 : : //! the single-material system, this argument can be left unspecified by
181 : : //! the calling code
182 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
183 : : //! for the single-material system, this argument can be left unspecified
184 : : //! by the calling code
185 : : //! \return Material speed of sound using the stiffened-gas EoS
186 : : // *************************************************************************
187 : : {
188 : 0 : alpha = std::max(1e-14,alpha);
189 : : // limiting pressure to near-zero
190 : 0 : auto apr_eff = std::max(alpha*
191 : 0 : min_eff_pressure(1e-4*std::abs(apr/alpha), arho, alpha), apr);
192 : :
193 : 0 : auto co1 = m_rho0*alpha*alpha/(arho*arho);
194 : 0 : auto co2 = alpha*(1.0+m_w)/arho;
195 : :
196 : 0 : tk::real ss = m_a*(m_r1*co1 - co2) * exp(-m_r1*alpha*m_rho0/arho)
197 : 0 : + m_b*(m_r2*co1 - co2) * exp(-m_r2*alpha*m_rho0/arho)
198 : 0 : + (1.0+m_w)*apr_eff/arho;
199 : :
200 : 0 : auto ss2 = ss;
201 : 0 : ss = std::sqrt(ss);
202 : :
203 : : // check sound speed divergence
204 [ - - ]: 0 : if (!std::isfinite(ss)) {
205 : 0 : std::cout << "Material-id: " << imat << std::endl;
206 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
207 : 0 : std::cout << "Partial density: " << arho << std::endl;
208 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
209 : 0 : std::cout << "Min allowed pres: " << alpha*min_eff_pressure(0.0, arho,
210 : 0 : alpha) << std::endl;
211 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed. "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
212 : : "ss^2: " + std::to_string(ss2) + ", material volume fraction: " +
213 : : std::to_string(alpha));
214 : : }
215 : :
216 : 0 : return ss;
217 : : }
218 : :
219 : : tk::real
220 : 0 : JWL::totalenergy(
221 : : tk::real arho,
222 : : tk::real u,
223 : : tk::real v,
224 : : tk::real w,
225 : : tk::real apr,
226 : : tk::real alpha,
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] arho Material partial density
232 : : //! \param[in] u X-velocity
233 : : //! \param[in] v Y-velocity
234 : : //! \param[in] w Z-velocity
235 : : //! \param[in] apr Material partial pressure
236 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
237 : : //! the single-material system, this argument can be left unspecified by
238 : : //! the calling code
239 : : //! \return Material specific total energy using the stiffened-gas EoS
240 : : // *************************************************************************
241 : : {
242 : : //// reference energy (input quantity, might need for calculation)
243 : : //tk::real e0 = a/r1*exp(-r1*rho0/rho) + b/r2*exp(-r2*rho0/rho);
244 : :
245 : 0 : tk::real arhoE = arho*intEnergy( arho/alpha, apr/alpha )
246 : 0 : + 0.5*arho*(u*u + v*v + w*w);
247 : :
248 : 0 : return arhoE;
249 : : }
250 : :
251 : : tk::real
252 : 0 : JWL::temperature(
253 : : tk::real arho,
254 : : tk::real u,
255 : : tk::real v,
256 : : tk::real w,
257 : : tk::real arhoE,
258 : : tk::real alpha,
259 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
260 : : // *************************************************************************
261 : : //! \brief Calculate material temperature from the material density, and
262 : : //! material specific total energy
263 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
264 : : //! \param[in] u X-velocity
265 : : //! \param[in] v Y-velocity
266 : : //! \param[in] w Z-velocity
267 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
268 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
269 : : //! the single-material system, this argument can be left unspecified by
270 : : //! the calling code
271 : : //! \return Material temperature using the stiffened-gas EoS
272 : : // *************************************************************************
273 : : {
274 : 0 : alpha = std::max(1e-14,alpha);
275 : 0 : tk::real rho = arho/alpha;
276 : :
277 : : //// reference energy (input quantity, might need for calculation)
278 : : //tk::real e0 = a/r1*exp(-r1*rho0/rho) + b/r2*exp(-r2*rho0/rho);
279 : :
280 : 0 : tk::real t = ((arhoE - 0.5*arho*(u*u + v*v + w*w))/arho + m_de -
281 : 0 : 1.0/m_rho0*( m_a/m_r1*exp(-m_r1*m_rho0/rho)
282 : 0 : + m_b/m_r2*exp(-m_r2*m_rho0/rho) ))/m_cv;
283 : :
284 : 0 : return t;
285 : : }
286 : :
287 : : tk::real
288 : 0 : JWL::min_eff_pressure(
289 : : tk::real min,
290 : : tk::real arho,
291 : : tk::real alpha ) const
292 : : // *************************************************************************
293 : : //! Compute the minimum allowed pressure
294 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
295 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
296 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
297 : : //! the single-material system, this argument can be left unspecified by
298 : : //! the calling code
299 : : //! \return Minimum pressure allowed by physical constraints
300 : : // *************************************************************************
301 : : {
302 : 0 : alpha = std::max(1e-14,alpha);
303 : 0 : auto co1 = m_rho0*alpha*alpha/(arho*arho);
304 : 0 : auto co2 = alpha*(1.0+m_w)/arho;
305 : :
306 : : // minimum pressure is constrained by zero soundspeed.
307 : 0 : auto apr = -(m_a*(m_r1*co1 - co2) * exp(-m_r1*alpha*m_rho0/arho)
308 : 0 : + m_b*(m_r2*co1 - co2) * exp(-m_r2*alpha*m_rho0/arho))
309 : 0 : * arho/(1.0+m_w);
310 : :
311 : 0 : return apr/alpha + min;
312 : : }
313 : :
314 : : tk::real
315 : 0 : JWL::intEnergy(
316 : : tk::real rho,
317 : : tk::real pr ) const
318 : : // *************************************************************************
319 : : //! \brief Calculate specific internal energy using the JWL equation of
320 : : //! state
321 : : //! \param[in] rho Material density
322 : : //! \param[in] pr Material pressure
323 : : //! \return Material internal energy calculated using the JWL EoS
324 : : //! \details By inverting Eqn. 1 in 'JWL Equation of State', Menikoff,
325 : : //! LA-UR-15-29536
326 : : // *************************************************************************
327 : : {
328 : 0 : tk::real e = - m_de + 1.0/m_w/rho*( pr
329 : 0 : - m_a*(1.0 - m_w*rho/m_r1/m_rho0)*exp(-m_r1*m_rho0/rho)
330 : 0 : - m_b*(1.0 - m_w*rho/m_r2/m_rho0)*exp(-m_r2*m_rho0/rho) );
331 : :
332 : 0 : return e;
333 : : }
334 : :
335 : : tk::real
336 : 0 : JWL::bisection(
337 : : tk::real a,
338 : : tk::real b,
339 : : tk::real p_known,
340 : : tk::real t_known ) const
341 : : // *************************************************************************
342 : : //! \brief Calculate density from known pressure and temperature using
343 : : //! bisection root finding method for JWL equation of state
344 : : //! \param[in] a Left density bound for root finding
345 : : //! \param[in] b Right density bound for root finding
346 : : //! \param[in] p_known Known pressure
347 : : //! \param[in] t_known Known temperature
348 : : //! \return Material density calculated by inverting JWL pressure equation
349 : : // *************************************************************************
350 : : {
351 : 0 : tk::real tol = 1e-12;
352 : 0 : std::size_t maxiter = 1000;
353 : 0 : std::size_t i(0);
354 : : tk::real c;
355 : 0 : tk::real root(0);
356 : 0 : std::size_t idebug = 0;
357 : :
358 : : // Ensure that original bounds contain root
359 [ - - ]: 0 : while ( p_known < PfromRT( a, t_known)) {
360 : 0 : b = a;
361 : 0 : a *= 1e-6;
362 : : }
363 [ - - ]: 0 : while ( p_known > PfromRT( b, t_known)) {
364 : 0 : a = b;
365 : 0 : b *= 1e6;
366 : : }
367 : :
368 : : // function to minimize: fcn = p_known - PfromRT
369 : : // bounds b > a
370 : :
371 [ - - ]: 0 : while (i < maxiter)
372 : : {
373 : 0 : c = (a + b)/2.0;
374 : 0 : auto fcn = p_known - PfromRT( c, t_known);
375 [ - - ]: 0 : if ( idebug == 1)
376 : : {
377 : 0 : std::cout << "Bisection iter: " << i << std::endl;
378 : 0 : std::cout << "fcn: " << fcn << std::endl;
379 : 0 : std::cout << "(b - a)/2.0: " << (b - a)/2.0 << std::endl;
380 : : }
381 : :
382 [ - - ][ - - ]: 0 : if ( std::abs(fcn) <= 1e-16 or (b - a)/2.0 < tol )
[ - - ]
383 : : {
384 : 0 : root = c;
385 : 0 : break;
386 : : }
387 : :
388 : 0 : i++;
389 : 0 : if ( static_cast< int > (std::copysign( 1.0, p_known - PfromRT( c, t_known) )) ==
390 [ - - ]: 0 : static_cast< int > (std::copysign( 1.0, p_known - PfromRT( a, t_known) )) )
391 : : {
392 : 0 : a = c;
393 : : }
394 : : else
395 : : {
396 : 0 : b = c;
397 : : }
398 : :
399 [ - - ]: 0 : if ( i == maxiter )
400 : : {
401 [ - - ][ - - ]: 0 : Throw("JWL Bisection for density failed to converge after iterations "
[ - - ][ - - ]
402 : : + std::to_string(i));
403 : : }
404 : :
405 : : }
406 : 0 : return root;
407 : : }
408 : :
409 : :
410 : : tk::real
411 : 0 : JWL::PfromRT(
412 : : tk::real rho,
413 : : tk::real T ) const
414 : : // *************************************************************************
415 : : //! \brief Calculate pressure from density and temperature using JWL
416 : : //! equation of state
417 : : //! \param[in] rho Material density
418 : : //! \param[in] T Material temperature
419 : : //! \return Material pressure calculated using the JWL EoS
420 : : //! \details From Eqn. 14 in 'JWL Equation of State', Menikoff, LA-UR-15-29536
421 : : // *************************************************************************
422 : : {
423 : 0 : return ( m_a*exp(-m_r1*m_rho0/rho) + m_b*exp(-m_r2*m_rho0/rho) +
424 : 0 : m_w*(m_cv*T*rho) );
425 : : }
|
