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