Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/EoS/GodunovRomenski.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 Godunov-Romenski equation of state for solids
9 : : \details This file defines functions for the Godunov-Romenski equation of
10 : : state for solids and a hydro EoS for aluminum. These functions were
11 : : taken from Example 1 of Barton, Philip T. "An interface-capturing
12 : : Godunov method for the simulation of compressible solid-fluid
13 : : problems." Journal of Computational Physics 390 (2019): 25-50.
14 : : */
15 : : // *****************************************************************************
16 : :
17 : : #include <cmath>
18 : : #include <iostream>
19 : : #include "Vector.hpp"
20 : : #include "EoS/GodunovRomenski.hpp"
21 : : #include "EoS/GetMatProp.hpp"
22 : :
23 : : // // Lapacke forward declarations
24 : : // extern "C" {
25 : :
26 : : // using lapack_int = long;
27 : :
28 : : // #define LAPACK_ROW_MAJOR 101
29 : :
30 : : // lapack_int LAPACKE_dgeev(int, char, char, lapack_int, double*, lapack_int,
31 : : // double*, double*, double*, lapack_int, double*, lapack_int );
32 : :
33 : : // }
34 : :
35 : : using inciter::GodunovRomenski;
36 : :
37 : 0 : GodunovRomenski::GodunovRomenski(
38 : : tk::real gamma,
39 : : tk::real mu,
40 : : tk::real rho0,
41 : : tk::real alpha,
42 : 0 : tk::real K0 ) :
43 : : m_gamma(gamma),
44 : : m_mu(mu),
45 : : m_rho0(rho0),
46 : : m_alpha(alpha),
47 : 0 : m_K0(K0)
48 : : // *************************************************************************
49 : : // Constructor
50 : : //! \param[in] gamma Ratio of specific heats
51 : : //! \param[in] mu Constant shear modulus
52 : : //! \param[in] rho0 Unstressed density of material
53 : : //! \param[in] alpha Alpha parameter
54 : : //! \param[in] K0 K0 parameter
55 : : // *************************************************************************
56 : 0 : { }
57 : :
58 : : tk::real
59 : 0 : GodunovRomenski::density(
60 : : tk::real pr,
61 : : tk::real ) const
62 : : // *************************************************************************
63 : : //! \brief Calculate density from the material pressure and temperature
64 : : //! using the GodunovRomenski equation of state
65 : : //! \param[in] pr Material pressure
66 : : //! \return Material density calculated using the cold compression pressure
67 : : // *************************************************************************
68 : : {
69 : : // Quick Newton
70 : 0 : tk::real rho = m_rho0;
71 : 0 : std::size_t maxiter = 50;
72 : 0 : tk::real tol = 1.0e-04;
73 : 0 : tk::real err = tol + 1;
74 [ - - ]: 0 : for (std::size_t iter=0; iter<maxiter; ++iter)
75 : : {
76 : 0 : tk::real p = pressure_coldcompr(rho) - pr;
77 : 0 : auto dpdrho = DpccDrho(rho);
78 : 0 : auto delta = p/dpdrho;
79 : 0 : rho -= delta;
80 : 0 : err = std::sqrt(std::pow(p,2.0));
81 [ - - ]: 0 : if (err < tol) break;
82 : : }
83 : 0 : return rho;
84 : : }
85 : :
86 : : tk::real
87 : 0 : GodunovRomenski::pressure(
88 : : tk::real arho,
89 : : tk::real u,
90 : : tk::real v,
91 : : tk::real w,
92 : : tk::real arhoE,
93 : : tk::real alpha,
94 : : std::size_t imat,
95 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
96 : : // *************************************************************************
97 : : //! \brief Calculate pressure from the material density, momentum, total energy
98 : : //! and the inverse deformation gradient tensor using the GodunovRomenski
99 : : //! equation of state
100 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
101 : : //! \param[in] u X-velocity
102 : : //! \param[in] v Y-velocity
103 : : //! \param[in] w Z-velocity
104 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
105 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
106 : : //! the single-material system, this argument can be left unspecified by
107 : : //! the calling code
108 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
109 : : //! for the single-material system, this argument can be left unspecified
110 : : //! by the calling code
111 : : //! \param[in] defgrad Material inverse deformation gradient tensor
112 : : //! (g_k). Default is 0, so that for the single-material system,
113 : : //! this argument can be left unspecified by the calling code
114 : : //! \return Material partial pressure (alpha_k * p_k) calculated using the
115 : : //! GodunovRomenski EoS
116 : : //! \details The material pressure consists of the thermal component and the
117 : : //! cold-compression component. Pressure-effects due to shear are stored in
118 : : //! the elastic Cauchy stress tensor, not in the pressure.
119 : : // *************************************************************************
120 : : {
121 : : // obtain elastic contribution to energy
122 : : std::array< std::array< tk::real, 3 >, 3 > devH;
123 [ - - ]: 0 : auto arhoEe = alpha*elasticEnergy(defgrad, devH);
124 : : // obtain cold compression contribution to energy
125 : 0 : auto rho = arho/alpha;
126 [ - - ]: 0 : auto arhoEc = alpha*coldcomprEnergy(rho);
127 : : // obtain thermal contribution to energy
128 : 0 : auto arhoEt = arhoE - arhoEe - arhoEc - 0.5*arho*(u*u + v*v + w*w);
129 : :
130 : : // use Mie-Gruneisen form of Godunov-Romenski for pressure
131 [ - - ]: 0 : auto partpressure = pressure_coldcompr(arho,alpha) + m_gamma*arhoEt;
132 : :
133 : : // check partial pressure divergence
134 [ - - ]: 0 : if (!std::isfinite(partpressure)) {
135 [ - - ][ - - ]: 0 : std::cout << "Material-id: " << imat << std::endl;
[ - - ]
136 [ - - ][ - - ]: 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
[ - - ]
137 [ - - ][ - - ]: 0 : std::cout << "Partial density: " << arho << std::endl;
[ - - ]
138 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) +
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
139 : : " has nan/inf partial pressure: " + std::to_string(partpressure) +
140 : : ", material volume fraction: " + std::to_string(alpha));
141 : : }
142 : :
143 : 0 : return partpressure;
144 : : }
145 : :
146 : : std::array< std::array< tk::real, 3 >, 3 >
147 : 0 : GodunovRomenski::CauchyStress(
148 : : tk::real,
149 : : tk::real,
150 : : tk::real,
151 : : tk::real,
152 : : tk::real,
153 : : tk::real alpha,
154 : : std::size_t /*imat*/,
155 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
156 : : // *************************************************************************
157 : : //! \brief Calculate the elastic Cauchy stress tensor from the material density,
158 : : //! momentum, total energy, and inverse deformation gradient tensor using the
159 : : //! GodunovRomenski equation of state
160 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
161 : : //! the single-material system, this argument can be left unspecified by
162 : : //! the calling code
163 : : // //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
164 : : // //! for the single-material system, this argument can be left unspecified
165 : : // //! by the calling code
166 : : //! \param[in] defgrad Material inverse deformation gradient tensor (g_k).
167 : : //! \return Elastic part of material Cauchy stress tensor (alpha_k * sigma_k)
168 : : //! calculated using the GodunovRomenski EoS
169 : : // *************************************************************************
170 : : {
171 : 0 : std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};
172 : :
173 : : // obtain elastic contribution to energy and subtract it from pmean
174 : : std::array< std::array< tk::real, 3 >, 3 > devH;
175 : :
176 : : // p_mean
177 [ - - ]: 0 : auto p_se = -elasticEnergy(defgrad, devH);
178 : 0 : auto pmean = alpha * p_se;
179 : :
180 : : // Pressure due to shear
181 : 0 : asig[0][0] = -pmean;
182 : 0 : asig[1][1] = -pmean;
183 : 0 : asig[2][2] = -pmean;
184 : :
185 : : // Add deviatoric component of Cauchy stress tensor
186 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
187 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
188 : 0 : asig[i][j] += 2.0*m_mu*alpha*devH[i][j];
189 : : }
190 : :
191 : 0 : return asig;
192 : : }
193 : :
194 : : tk::real
195 : 0 : GodunovRomenski::soundspeed(
196 : : tk::real arho,
197 : : tk::real apr,
198 : : tk::real alpha,
199 : : std::size_t imat,
200 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
201 : : // *************************************************************************
202 : : //! Calculate speed of sound from the material density and material pressure
203 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
204 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
205 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
206 : : //! the single-material system, this argument can be left unspecified by
207 : : //! the calling code
208 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
209 : : //! for the single-material system, this argument can be left unspecified
210 : : //! by the calling code
211 : : //! (alpha * sigma_ij * n_j) projected onto the normal vector. Default is 0,
212 : : //! so that for the single-material system, this argument can be left
213 : : //! unspecified by the calling code
214 : : // //! \param[in] defgrad Material inverse deformation gradient tensor
215 : : // //! (g_k) with the first dimension aligned to direction in which
216 : : // //! wave speeds are required. Default is 0, so that for the single-material
217 : : // //! system, this argument can be left unspecified by the calling code
218 : : //! \return Material speed of sound using the GodunovRomenski EoS
219 : : // *************************************************************************
220 : : {
221 : 0 : tk::real a = 0.0;
222 : :
223 : : // Hydro contribution
224 : 0 : tk::real rho = arho/alpha;
225 : 0 : auto p_cc = pressure_coldcompr(arho, alpha);
226 : 0 : a += std::max( 1.0e-15, DpccDrho(rho) + (m_gamma+1.0) * (apr - p_cc)/arho );
227 : : // in the above expression, shear pressure is not included in apr in the first
228 : : // place, so should not subtract it
229 : :
230 : : // Shear contribution
231 : 0 : a += (4.0/3.0) * alpha * m_mu / arho;
232 : :
233 : : // Compute square root
234 : 0 : a = std::sqrt(a);
235 : :
236 : : // check sound speed divergence
237 [ - - ]: 0 : if (!std::isfinite(a)) {
238 : 0 : std::cout << "Material-id: " << imat << std::endl;
239 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
240 : 0 : std::cout << "Partial density: " << arho << std::endl;
241 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
242 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
243 : : + std::to_string(a) + ", material volume fraction: " +
244 : : std::to_string(alpha));
245 : : }
246 : :
247 : 0 : return a;
248 : : }
249 : :
250 : : tk::real
251 : 0 : GodunovRomenski::shearspeed(
252 : : tk::real arho,
253 : : tk::real alpha,
254 : : std::size_t imat ) const
255 : : // *************************************************************************
256 : : //! Calculate speed of sound from the material density and material pressure
257 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
258 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
259 : : //! the single-material system, this argument can be left unspecified by
260 : : //! the calling code
261 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
262 : : //! for the single-material system, this argument can be left unspecified
263 : : //! by the calling code
264 : : //! \return Material shear-wave speed speed using the SmallShearSolid EoS
265 : : // *************************************************************************
266 : : {
267 : : // Approximate shear-wave speed. Ref. Barton, P. T. (2019).
268 : : // An interface-capturing Godunov method for the simulation of compressible
269 : : // solid-fluid problems. Journal of Computational Physics, 390, 25-50.
270 : 0 : tk::real a = std::sqrt(alpha*m_mu/arho);
271 : :
272 : : // check shear-wave speed divergence
273 [ - - ]: 0 : if (!std::isfinite(a)) {
274 : 0 : std::cout << "Material-id: " << imat << std::endl;
275 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
276 : 0 : std::cout << "Partial density: " << arho << std::endl;
277 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf shear-wave speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
278 : : + std::to_string(a) + ", material volume fraction: " +
279 : : std::to_string(alpha));
280 : : }
281 : :
282 : 0 : return a;
283 : : }
284 : :
285 : : tk::real
286 : 0 : GodunovRomenski::totalenergy(
287 : : tk::real rho,
288 : : tk::real u,
289 : : tk::real v,
290 : : tk::real w,
291 : : tk::real pr,
292 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
293 : : // *************************************************************************
294 : : //! \brief Calculate material specific total energy from the material
295 : : //! density, momentum and material pressure
296 : : //! \param[in] rho Material density
297 : : //! \param[in] u X-velocity
298 : : //! \param[in] v Y-velocity
299 : : //! \param[in] w Z-velocity
300 : : //! \param[in] pr Material pressure
301 : : //! \param[in] defgrad Material inverse deformation gradient tensor
302 : : //! g_k. Default is 0, so that for the single-material system,
303 : : //! this argument can be left unspecified by the calling code
304 : : //! \return Material specific total energy using the GodunovRomenski EoS
305 : : // *************************************************************************
306 : : {
307 : : // obtain thermal and kinetic energy
308 [ - - ]: 0 : auto pt = pr - pressure_coldcompr(rho);
309 : : // in the above expression, shear pressure is not included in pr in the first
310 : : // place, so should not subtract it
311 : 0 : auto rhoEh = pt/m_gamma + 0.5*rho*(u*u + v*v + w*w);
312 : : // obtain elastic contribution to energy
313 : : std::array< std::array< tk::real, 3 >, 3 > devH;
314 [ - - ]: 0 : auto rhoEe = elasticEnergy(defgrad, devH);
315 [ - - ]: 0 : auto rhoEc = coldcomprEnergy(rho);
316 : :
317 : 0 : return (rhoEh + rhoEe + rhoEc);
318 : : }
319 : :
320 : : tk::real
321 : 0 : GodunovRomenski::temperature(
322 : : tk::real,
323 : : tk::real,
324 : : tk::real,
325 : : tk::real,
326 : : tk::real,
327 : : tk::real,
328 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
329 : : // *************************************************************************
330 : : //! \brief Calculate material temperature from the material density, and
331 : : //! material specific total energy
332 : : // //! \param[in] arho Material partial density (alpha_k * rho_k)
333 : : // //! \param[in] u X-velocity
334 : : // //! \param[in] v Y-velocity
335 : : // //! \param[in] w Z-velocity
336 : : // //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
337 : : // //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
338 : : // //! the single-material system, this argument can be left unspecified by
339 : : // //! the calling code
340 : : // //! \param[in] defgrad Material inverse deformation gradient tensor
341 : : // //! (g_k). Default is 0, so that for the single-material system,
342 : : // //! this argument can be left unspecified by the calling code
343 : : //! \return Material temperature using the GodunovRomenski EoS
344 : : // *************************************************************************
345 : : {
346 : : // Temperature as a function of energy is not known.
347 : : // So we just set a value.
348 : 0 : tk::real t = 300.0;
349 : :
350 : 0 : return t;
351 : : }
352 : :
353 : : tk::real
354 : 0 : GodunovRomenski::min_eff_pressure(
355 : : tk::real min,
356 : : tk::real arho,
357 : : tk::real alpha ) const
358 : : // *************************************************************************
359 : : //! Compute the minimum allowed pressure
360 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
361 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
362 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
363 : : //! \return Minimum pressure allowed by physical constraints
364 : : // *************************************************************************
365 : : {
366 : : // minimum pressure is constrained by zero soundspeed.
367 : 0 : auto rho = arho/alpha;
368 : : return min
369 : 0 : - rho/(m_gamma+1.0) * DpccDrho(rho)
370 : 0 : + pressure_coldcompr(arho, alpha)/alpha;
371 : : }
372 : :
373 : : tk::real
374 : 0 : GodunovRomenski::coldcomprEnergy( tk::real rho ) const
375 : : // *************************************************************************
376 : : //! \brief Calculate cold-compression contribution to material energy from the
377 : : //! material density
378 : : //! \param[in] rho Material density
379 : : //! \return Material cold compression energy using the GodunovRomenski EoS
380 : : // *************************************************************************
381 : : {
382 : 0 : auto rrho0a = std::pow(rho/m_rho0, m_alpha);
383 : 0 : return ( rho * m_K0/(2.0*m_rho0*m_alpha*m_alpha) * (rrho0a-1.0)*(rrho0a-1.0) );
384 : : }
385 : :
386 : : tk::real
387 : 0 : GodunovRomenski::pressure_coldcompr(
388 : : tk::real arho,
389 : : tk::real alpha ) const
390 : : // *************************************************************************
391 : : //! \brief Calculate cold-compression contribution to material pressure from the
392 : : //! material density
393 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
394 : : //! \param[in] alpha Material volume fraction. Default is 1.0.
395 : : //! \return Material cold compression pressure using the GodunovRomenski EoS
396 : : // *************************************************************************
397 : : {
398 : 0 : auto rho = arho/alpha;
399 : 0 : auto rrho0a = std::pow(rho/m_rho0, m_alpha);
400 : 0 : return alpha * (m_K0/m_alpha * (rrho0a*rho/m_rho0) * (rrho0a-1.0));
401 : : }
402 : :
403 : : tk::real
404 : 0 : GodunovRomenski::elasticEnergy(
405 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad,
406 : : std::array< std::array< tk::real, 3 >, 3 >& devH ) const
407 : : // *************************************************************************
408 : : //! \brief Calculate elastic contribution to material energy from the material
409 : : //! density, and deformation gradient tensor
410 : : //! \param[in] defgrad Material inverse deformation gradient tensor
411 : : //! \param[in/out] devH Deviatoric part of the Hensky tensor
412 : : //! \return Material elastic energy using the GodunovRomenski EoS
413 : : //! \details This function returns the material elastic energy, and also stores
414 : : //! the elastic shear distortion for further use
415 : : // *************************************************************************
416 : : {
417 : : // Compute deviatoric part of Hencky tensor
418 : 0 : devH = tk::getDevHencky(defgrad);
419 : :
420 : : // Compute elastic energy
421 : 0 : tk::real rhoEe = 0.0;
422 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
423 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
424 : 0 : rhoEe += m_mu*devH[i][j]*devH[i][j];
425 : :
426 : 0 : return rhoEe;
427 : : }
428 : :
429 : : tk::real
430 : 0 : GodunovRomenski::DpccDrho( tk::real rho ) const
431 : : // *************************************************************************
432 : : //! Calculate the derivative of the cold compression pressure wrt. density
433 : : //! \param[in] rho Material partial density (alpha_k * rho_k)
434 : : //! \return Derivative of the cold compression pressure wrt. density
435 : : // *************************************************************************
436 : : {
437 : 0 : auto rrho0a = std::pow(rho/m_rho0, m_alpha);
438 : 0 : return m_K0/(m_rho0*m_alpha) *
439 : 0 : ((2.0*m_alpha+1.0)*(rrho0a*rrho0a) - (m_alpha+1.0)*rrho0a);
440 : : }
|
