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