1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
// *****************************************************************************
/*!
  \file      src/PDE/EoS/JWL.cpp
  \copyright 2012-2015 J. Bakosi,
             2016-2018 Los Alamos National Security, LLC.,
             2019-2021 Triad National Security, LLC.
             All rights reserved. See the LICENSE file for details.
  \brief     Jones, Wilkins, and Lee (JWL) equation of state
  \details   This file defines functions for the JWL equation of
             state for the compressible flow equations. These functions are
             taken from 'JWL Equation of State', Menikoff, LA-UR-15-29536.
*/
// *****************************************************************************

#include <cmath>
#include <iostream>
#include "EoS/JWL.hpp"

using inciter::JWL;

JWL::JWL( tk::real w, tk::real cv, tk::real rho0, tk::real de, tk::real rhor,
  tk::real tr, tk::real pr, tk::real A, tk::real B, tk::real R1, tk::real R2 ) :
  m_w(w),
  m_cv(cv),
  m_rho0(rho0),
  m_de(de),
  m_rhor(rhor),
  m_tr(tr),
  m_pr(pr),
  m_a(A),
  m_b(B),
  m_r1(R1),
  m_r2(R2)
// *************************************************************************
//  Constructor
//! \param[in] w Grueneisen coefficient
//! \param[in] cv Specific heat at constant volume
//! \param[in] rho0 Density of initial state
//! \param[in] de Heat of detonation for products. For reactants, it is
//!   chosen such that the ambient internal energy (e0) is 0.
//! \param[in] rhor Density of reference state
//! \param[in] tr Temperature of reference state
//! \param[in] pr Pressure of reference state
//! \param[in] A Parameter A
//! \param[in] B Parameter B
//! \param[in] R1 Parameter R1
//! \param[in] R2 Parameter R2
// *************************************************************************
{
  // reference density provided
  if (m_tr < 1e-8) {
    // reference internal energy
    auto er = intEnergy(rhor, pr);
    // reference temperature from Eqn (15)
    m_tr = 1.0/m_cv * (er + de -
      (m_a/m_r1*exp(-m_r1*m_rho0/m_rhor) +
       m_b/m_r2*exp(-m_r2*m_rho0/m_rhor)) / m_rho0);
  }
  // reference temperature provided
  else
  {
    m_rhor = density(m_pr, m_tr);
  }
}

tk::real
JWL::density(
  tk::real pr,
  tk::real temp ) const
// *************************************************************************
//! \brief Calculate density from the material pressure and temperature
//!   using the stiffened-gas equation of state
//! \param[in] pr Material pressure
//! \param[in] temp Material temperature
//! \return Material density calculated using the stiffened-gas EoS
// *************************************************************************
{
  tk::real r_guessL = 1e-4*m_rho0;  // left density bound
  tk::real r_guessR = 1e2*m_rho0;   // right density bound
  tk::real rho;

  rho = bisection( r_guessL, r_guessR, pr, temp );

  return rho;
}


tk::real
JWL::pressure(
  tk::real arho,
  tk::real u,
  tk::real v,
  tk::real w,
  tk::real arhoE,
  tk::real alpha,
  std::size_t imat,
  const std::array< std::array< tk::real, 3 >, 3 >& ) const
// *************************************************************************
//! \brief Calculate pressure from the material density, momentum and total
//!   energy using the stiffened-gas equation of state
//! \param[in] arho Material partial density (alpha_k * rho_k)
//! \param[in] u X-velocity
//! \param[in] v Y-velocity
//! \param[in] w Z-velocity
//! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
//! \param[in] alpha Material volume fraction. Default is 1.0, so that for
//!   the single-material system, this argument can be left unspecified by
//!   the calling code
//! \param[in] imat Material-id who's EoS is required. Default is 0, so that
//!   for the single-material system, this argument can be left unspecified
//!   by the calling code
//! \return Material partial pressure (alpha_k * p_k) calculated using the
//!   stiffened-gas EoS
//! \details From Eqn. 1 in 'JWL Equation of State', Menikoff, LA-UR-15-29536
// *************************************************************************
{
  // specific internal energy
  tk::real e = (arhoE - 0.5*arho*(u*u + v*v + w*w))/arho;

  //// reference energy (input quantity, might need for calculation)
  //tk::real e0 = a/r1*exp(-r1*rho0/rho) + b/r2*exp(-r2*rho0/rho);

  tk::real partpressure =
    m_a*(alpha - m_w*arho/(m_rho0*m_r1))*exp(-m_r1*alpha*m_rho0/arho) +
    m_b*(alpha - m_w*arho/(m_rho0*m_r2))*exp(-m_r2*alpha*m_rho0/arho) +
    m_w*arho*(e + m_de);

  // check partial pressure divergence
  if (!std::isfinite(partpressure)) {
    std::cout << "Material-id:      " << imat << std::endl;
    std::cout << "Volume-fraction:  " << alpha << std::endl;
    std::cout << "Partial density:  " << arho << std::endl;
    std::cout << "Total energy:     " << arhoE << std::endl;
    std::cout << "Velocity:         " << u << ", " << v << ", " << w
      << std::endl;
    Throw("Material-" + std::to_string(imat) +
      " has nan/inf partial pressure: " + std::to_string(partpressure) +
      ", material volume fraction: " + std::to_string(alpha));
  }

  return partpressure;
}

std::array< std::array< tk::real, 3 >, 3 >
JWL::CauchyStress(
  tk::real,
  tk::real,
  tk::real,
  tk::real,
  tk::real,
  tk::real,
  std::size_t,
  const std::array< std::array< tk::real, 3 >, 3 >& ) const
// *************************************************************************
//! \brief Calculate the Cauchy stress tensor from the material density,
//!   momentum, and total energy
//! \return Material Cauchy stress tensor (alpha_k * sigma_k)
// *************************************************************************
{
  std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};

  // No elastic contribution

  return asig;
}

tk::real
JWL::soundspeed(
  tk::real arho,
  tk::real apr,
  tk::real alpha,
  std::size_t imat,
  const std::array< std::array< tk::real, 3 >, 3 >&,
  const std::array< tk::real, 3 >&,
  const std::array< tk::real, 3 >& ) const
// *************************************************************************
//! Calculate speed of sound from the material density and material pressure
//! \param[in] arho Material partial density (alpha_k * rho_k)
//! \param[in] apr Material partial pressure (alpha_k * p_k)
//! \param[in] alpha Material volume fraction. Default is 1.0, so that for
//!   the single-material system, this argument can be left unspecified by
//!   the calling code
//! \param[in] imat Material-id who's EoS is required. Default is 0, so that
//!   for the single-material system, this argument can be left unspecified
//!   by the calling code
//! \return Material speed of sound using the stiffened-gas EoS
// *************************************************************************
{
  // limiting pressure to near-zero
  auto apr_eff = std::max(alpha*
    min_eff_pressure(1e-4*std::abs(apr/alpha), arho, alpha), apr);

  auto co1 = m_rho0*alpha*alpha/(arho*arho);
  auto co2 = alpha*(1.0+m_w)/arho;

  tk::real ss = m_a*(m_r1*co1 - co2) * exp(-m_r1*alpha*m_rho0/arho)
              + m_b*(m_r2*co1 - co2) * exp(-m_r2*alpha*m_rho0/arho)
              + (1.0+m_w)*apr_eff/arho;

  auto ss2 = ss;
  ss = std::sqrt(ss);

  // check sound speed divergence
  if (!std::isfinite(ss)) {
    std::cout << "Material-id:      " << imat << std::endl;
    std::cout << "Volume-fraction:  " << alpha << std::endl;
    std::cout << "Partial density:  " << arho << std::endl;
    std::cout << "Partial pressure: " << apr << std::endl;
    std::cout << "Min allowed pres: " << alpha*min_eff_pressure(0.0, arho,
      alpha) << std::endl;
    Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed. "
      "ss^2: " + std::to_string(ss2) + ", material volume fraction: " +
      std::to_string(alpha));
  }

  return ss;
}

tk::real
JWL::totalenergy(
  tk::real rho,
  tk::real u,
  tk::real v,
  tk::real w,
  tk::real pr,
  const std::array< std::array< tk::real, 3 >, 3 >& ) const
// *************************************************************************
//! \brief Calculate material specific total energy from the material
//!   density, momentum and material pressure
//! \param[in] rho Material density
//! \param[in] u X-velocity
//! \param[in] v Y-velocity
//! \param[in] w Z-velocity
//! \param[in] pr Material pressure
//! \return Material specific total energy using the stiffened-gas EoS
// *************************************************************************
{
  //// reference energy (input quantity, might need for calculation)
  //tk::real e0 = a/r1*exp(-r1*rho0/rho) + b/r2*exp(-r2*rho0/rho);

  tk::real rhoE = rho*intEnergy( rho, pr )
                + 0.5*rho*(u*u + v*v + w*w);

  return rhoE;
}

tk::real
JWL::temperature(
  tk::real arho,
  tk::real u,
  tk::real v,
  tk::real w,
  tk::real arhoE,
  tk::real alpha,
  const std::array< std::array< tk::real, 3 >, 3 >& ) const
// *************************************************************************
//! \brief Calculate material temperature from the material density, and
//!   material specific total energy
//! \param[in] arho Material partial density (alpha_k * rho_k)
//! \param[in] u X-velocity
//! \param[in] v Y-velocity
//! \param[in] w Z-velocity
//! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
//! \param[in] alpha Material volume fraction. Default is 1.0, so that for
//!   the single-material system, this argument can be left unspecified by
//!   the calling code
//! \return Material temperature using the stiffened-gas EoS
// *************************************************************************
{
  tk::real rho = arho/alpha;

  //// reference energy (input quantity, might need for calculation)
  //tk::real e0 = a/r1*exp(-r1*rho0/rho) + b/r2*exp(-r2*rho0/rho);

  tk::real t = ((arhoE - 0.5*arho*(u*u + v*v + w*w))/arho + m_de -
    1.0/m_rho0*( m_a/m_r1*exp(-m_r1*m_rho0/rho)
               + m_b/m_r2*exp(-m_r2*m_rho0/rho) ))/m_cv;

  return t;
}

tk::real
JWL::min_eff_pressure(
  tk::real min,
  tk::real arho,
  tk::real alpha ) const
// *************************************************************************
//! Compute the minimum allowed pressure
//! \param[in] min Numerical threshold above which pressure needs to be limited
//! \param[in] arho Material partial density (alpha_k * rho_k)
//! \param[in] alpha Material volume fraction. Default is 1.0, so that for
//!   the single-material system, this argument can be left unspecified by
//!   the calling code
//! \return Minimum pressure allowed by physical constraints
// *************************************************************************
{
  auto co1 = m_rho0*alpha*alpha/(arho*arho);
  auto co2 = alpha*(1.0+m_w)/arho;

  // minimum pressure is constrained by zero soundspeed.
  auto apr = -(m_a*(m_r1*co1 - co2) * exp(-m_r1*alpha*m_rho0/arho)
             + m_b*(m_r2*co1 - co2) * exp(-m_r2*alpha*m_rho0/arho))
    * arho/(1.0+m_w);

  return apr/alpha + min;
}

tk::real
JWL::intEnergy(
  tk::real rho,
  tk::real pr ) const
// *************************************************************************
//! \brief Calculate specific internal energy using the JWL equation of
//!   state
//! \param[in] rho Material density
//! \param[in] pr Material pressure
//! \return Material internal energy calculated using the JWL EoS
//! \details By inverting Eqn. 1 in 'JWL Equation of State', Menikoff,
//!   LA-UR-15-29536
// *************************************************************************
{
  tk::real e = - m_de + 1.0/m_w/rho*( pr
                - m_a*(1.0 - m_w*rho/m_r1/m_rho0)*exp(-m_r1*m_rho0/rho)
                - m_b*(1.0 - m_w*rho/m_r2/m_rho0)*exp(-m_r2*m_rho0/rho) );

  return e;
}

tk::real
JWL::bisection(
  tk::real a,
  tk::real b,
  tk::real p_known,
  tk::real t_known ) const
// *************************************************************************
//! \brief Calculate density from known pressure and temperature using
//!   bisection root finding method for JWL equation of state
//! \param[in] a Left density bound for root finding
//! \param[in] b Right density bound for root finding
//! \param[in] p_known Known pressure
//! \param[in] t_known Known temperature
//! \return Material density calculated by inverting JWL pressure equation
// *************************************************************************
{
  tk::real tol = 1e-12;
  std::size_t maxiter = 1000;
  std::size_t i(0);
  tk::real c;
  tk::real root(0);
  std::size_t idebug = 0;<--- Assignment 'idebug=0', assigned value is 0
  auto a_o = a;
  auto b_o = b;

  // function to minimize: fcn = p_known - PfromRT
  // bounds b > a

  while (i < maxiter)
  {
    c = (a + b)/2.0;
    auto fcn = p_known - PfromRT( c, t_known);
    if ( idebug == 1)<--- Condition 'idebug==1' is always false
    {
      std::cout << "Bisection iter:      " << i << std::endl;
      std::cout << "fcn:  " << fcn << std::endl;
      std::cout << "(b - a)/2.0: " << (b - a)/2.0 << std::endl;
    }

    if ( std::abs(fcn) <= 1e-16 or (b - a)/2.0 < tol )
    {
      root = c;
      break;
    }

    i++;
    if ( static_cast< int > (std::copysign( 1.0, p_known - PfromRT( c, t_known) )) ==
         static_cast< int > (std::copysign( 1.0, p_known - PfromRT( a, t_known) )) )
    {
      a = c;
    }
    else
    {
      b = c;
    }

    if ( i == maxiter )
    {
      Throw("JWL Bisection for density failed to converge after iterations "
      + std::to_string(i));
    }
    if (std::abs(root-a_o) < 1e-16 || std::abs(root-b_o) < 1e-16)
    {
      Throw("JWL bisection for density resulted in left/right bound as "
      "solution. Extend bounds for correctness");
    }

  }
  return root;
}


tk::real
JWL::PfromRT(
  tk::real rho,
  tk::real T ) const
// *************************************************************************
//! \brief Calculate pressure from density and temperature using JWL
//!   equation of state
//! \param[in] rho Material density
//! \param[in] T Material temperature
//! \return Material pressure calculated using the JWL EoS
//! \details From Eqn. 14 in 'JWL Equation of State', Menikoff, LA-UR-15-29536
// *************************************************************************
{
  return ( m_a*exp(-m_r1*m_rho0/rho) + m_b*exp(-m_r2*m_rho0/rho) +
    m_w*(m_cv*T*rho) );
}