Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/EoS/WilkinsAluminum.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 Wilkins equation of state for aluminum
9 : : \details This file defines functions for the Wilkins equation of
10 : : state for solids and a hydro EoS for aluminum. These functions were
11 : : taken from Example 4 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/WilkinsAluminum.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::WilkinsAluminum;
42 : :
43 : 0 : WilkinsAluminum::WilkinsAluminum(
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 : WilkinsAluminum::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 : WilkinsAluminum::density(
73 : : tk::real pr,
74 : : tk::real ) const
75 : : // *************************************************************************
76 : : //! \brief Calculate density from the material pressure and temperature
77 : : //! using the WilkinsAluminum equation of state
78 : : //! \param[in] pr Material pressure
79 : : // //! \param[in] temp Material temperature
80 : : //! \return Material density calculated using the WilkinsAluminum 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 : WilkinsAluminum::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 WilkinsAluminum
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 : : //! WilkinsAluminum 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 : WilkinsAluminum::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 : : //! WilkinsAluminum 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 WilkinsAluminum 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 : WilkinsAluminum::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*/ ) const
209 : : // *************************************************************************
210 : : //! Calculate speed of sound from the material density and material pressure
211 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
212 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
213 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
214 : : //! the single-material system, this argument can be left unspecified by
215 : : //! the calling code
216 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
217 : : //! for the single-material system, this argument can be left unspecified
218 : : //! by the calling code
219 : : //! (alpha * sigma_ij * n_j) projected onto the normal vector. Default is 0,
220 : : //! so that for the single-material system, this argument can be left
221 : : //! unspecified by the calling code
222 : : // //! \param[in] defgrad Material inverse deformation gradient tensor
223 : : // //! (g_k) with the first dimension aligned to direction in which
224 : : // //! wave speeds are required. Default is 0, so that for the single-material
225 : : // //! system, this argument can be left unspecified by the calling code
226 : : //! \return Material speed of sound using the WilkinsAluminum EoS
227 : : // *************************************************************************
228 : : {
229 : 0 : tk::real a = 0.0;
230 : :
231 : : // Hydro contribution
232 : 0 : tk::real rho0 = m_rho0;
233 : 0 : tk::real rho = arho/alpha;
234 : 0 : a += std::max( 1.0e-15, 6*e2*std::pow(rho/rho0,2.0)/rho0
235 : 0 : + 2*e3*rho/(rho0*rho0) - e5/rho0 );
236 : :
237 : : // Shear contribution
238 : 0 : a += (4.0/3.0) * m_mu / (arho/alpha);
239 : :
240 : : // Compute square root
241 : 0 : a = std::sqrt(a);
242 : :
243 : : // check sound speed divergence
244 [ - - ]: 0 : if (!std::isfinite(a)) {
245 : 0 : std::cout << "Material-id: " << imat << std::endl;
246 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
247 : 0 : std::cout << "Partial density: " << arho << std::endl;
248 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
249 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
250 : : + std::to_string(a) + ", material volume fraction: " +
251 : : std::to_string(alpha));
252 : : }
253 : :
254 : 0 : return a;
255 : : }
256 : :
257 : : tk::real
258 : 0 : WilkinsAluminum::shearspeed(
259 : : tk::real arho,
260 : : tk::real alpha,
261 : : std::size_t imat ) const
262 : : // *************************************************************************
263 : : //! Calculate speed of sound from the material density and material pressure
264 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
265 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
266 : : //! the single-material system, this argument can be left unspecified by
267 : : //! the calling code
268 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
269 : : //! for the single-material system, this argument can be left unspecified
270 : : //! by the calling code
271 : : //! \return Material shear-wave speed speed using the SmallShearSolid EoS
272 : : // *************************************************************************
273 : : {
274 : : // Approximate shear-wave speed. Ref. Barton, P. T. (2019).
275 : : // An interface-capturing Godunov method for the simulation of compressible
276 : : // solid-fluid problems. Journal of Computational Physics, 390, 25-50.
277 : 0 : tk::real a = std::sqrt(alpha*m_mu/arho);
278 : :
279 : : // check shear-wave speed divergence
280 [ - - ]: 0 : if (!std::isfinite(a)) {
281 : 0 : std::cout << "Material-id: " << imat << std::endl;
282 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
283 : 0 : std::cout << "Partial density: " << arho << std::endl;
284 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf shear-wave speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
285 : : + std::to_string(a) + ", material volume fraction: " +
286 : : std::to_string(alpha));
287 : : }
288 : :
289 : 0 : return a;
290 : : }
291 : :
292 : : tk::real
293 : 0 : WilkinsAluminum::totalenergy(
294 : : tk::real rho,
295 : : tk::real u,
296 : : tk::real v,
297 : : tk::real w,
298 : : tk::real,
299 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
300 : : // *************************************************************************
301 : : //! \brief Calculate material specific total energy from the material
302 : : //! density, momentum and material pressure
303 : : //! \param[in] rho Material density
304 : : //! \param[in] u X-velocity
305 : : //! \param[in] v Y-velocity
306 : : //! \param[in] w Z-velocity
307 : : // //! \param[in] pr Material pressure
308 : : //! \param[in] defgrad Material inverse deformation gradient tensor
309 : : //! g_k. Default is 0, so that for the single-material system,
310 : : //! this argument can be left unspecified by the calling code
311 : : //! \return Material specific total energy using the WilkinsAluminum EoS
312 : : // *************************************************************************
313 : : {
314 : : // obtain hydro contribution to energy
315 : 0 : tk::real rho0 = m_rho0;
316 : 0 : tk::real rhoEh = (e1+e2*std::pow(rho/rho0,2.0)+e3*(rho/rho0)
317 : 0 : +e4*std::pow(rho/rho0,-1.0)-e5*std::log(rho/rho0))/rho0
318 : 0 : + 0.5*rho*(u*u + v*v + w*w);
319 : : // obtain elastic contribution to energy
320 : : std::array< std::array< tk::real, 3 >, 3 > devH;
321 [ - - ]: 0 : tk::real rhoEe = elasticEnergy(defgrad, devH);
322 : :
323 : 0 : return (rhoEh + rhoEe);
324 : : }
325 : :
326 : : tk::real
327 : 0 : WilkinsAluminum::temperature(
328 : : tk::real,
329 : : tk::real,
330 : : tk::real,
331 : : tk::real,
332 : : tk::real,
333 : : tk::real,
334 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
335 : : // *************************************************************************
336 : : //! \brief Calculate material temperature from the material density, and
337 : : //! material specific total energy
338 : : // //! \param[in] arho Material partial density (alpha_k * rho_k)
339 : : // //! \param[in] u X-velocity
340 : : // //! \param[in] v Y-velocity
341 : : // //! \param[in] w Z-velocity
342 : : // //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
343 : : // //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
344 : : // //! the single-material system, this argument can be left unspecified by
345 : : // //! the calling code
346 : : // //! \param[in] defgrad Material inverse deformation gradient tensor
347 : : // //! (g_k). Default is 0, so that for the single-material system,
348 : : // //! this argument can be left unspecified by the calling code
349 : : //! \return Material temperature using the WilkinsAluminum EoS
350 : : // *************************************************************************
351 : : {
352 : : // Temperature does not directly contribute to energy
353 : : // So we just set a value.
354 : 0 : tk::real t = 300.0;
355 : :
356 : 0 : return t;
357 : : }
358 : :
359 : : tk::real
360 : 0 : WilkinsAluminum::min_eff_pressure(
361 : : tk::real min,
362 : : tk::real,
363 : : tk::real ) const
364 : : // *************************************************************************
365 : : //! Compute the minimum allowed pressure
366 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
367 : : //! \return Minimum pressure allowed by physical constraints
368 : : // *************************************************************************
369 : : {
370 : : // minimum pressure is constrained by zero soundspeed.
371 : 0 : return min;
372 : : }
373 : :
374 : : tk::real
375 : 0 : WilkinsAluminum::elasticEnergy(
376 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad,
377 : : std::array< std::array< tk::real, 3 >, 3 >& devH ) const
378 : : // *************************************************************************
379 : : //! \brief Calculate elastic contribution to material energy from the material
380 : : //! density, and deformation gradient tensor
381 : : //! \param[in] defgrad Material inverse deformation gradient tensor
382 : : //! \param[in/out] devH Deviatoric part of the Hensky tensor
383 : : //! \return Material elastic energy using the WilkinsAluminum EoS
384 : : //! \details This function returns the material elastic energy, and also stores
385 : : //! the elastic shear distortion for further use
386 : : // *************************************************************************
387 : : {
388 : : // Compute deviatoric part of Hencky tensor
389 : 0 : devH = tk::getDevHencky(defgrad);
390 : :
391 : : // Compute elastic energy
392 : 0 : tk::real rhoEe = 0.0;
393 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
394 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
395 : 0 : rhoEe += m_mu*devH[i][j]*devH[i][j];
396 : :
397 : 0 : return rhoEe;
398 : : }
|
