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