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