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 [ - - ]: 0 : partpressure = std::max(min_eff_pressure(1e-10, arho, alpha), partpressure);
143 : :
144 : : //// check partial pressure divergence
145 : : //if (!std::isfinite(partpressure)) {
146 : : // std::cout << "Material-id: " << imat << std::endl;
147 : : // std::cout << "Volume-fraction: " << alpha << std::endl;
148 : : // std::cout << "Partial density: " << arho << std::endl;
149 : : // Throw("Material-" + std::to_string(imat) +
150 : : // " has nan/inf partial pressure: " + std::to_string(partpressure) +
151 : : // ", material volume fraction: " + std::to_string(alpha));
152 : : //}
153 : :
154 : 0 : return partpressure;
155 : : }
156 : :
157 : : std::array< std::array< tk::real, 3 >, 3 >
158 : 0 : WilkinsAluminum::CauchyStress(
159 : : tk::real alpha,
160 : : std::size_t /*imat*/,
161 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
162 : : // *************************************************************************
163 : : //! \brief Calculate the elastic Cauchy stress tensor from the material
164 : : //! inverse deformation gradient tensor using the WilkinsAluminum EOS
165 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
166 : : //! the single-material system, this argument can be left unspecified by
167 : : //! the calling code
168 : : // //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
169 : : // //! for the single-material system, this argument can be left unspecified
170 : : // //! by the calling code
171 : : //! \param[in] defgrad Material inverse deformation gradient tensor (g_k).
172 : : //! \return Material Cauchy stress tensor (alpha_k * sigma_k) calculated using
173 : : //! the WilkinsAluminum EoS
174 : : // *************************************************************************
175 : : {
176 : 0 : std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};
177 : :
178 : : // obtain elastic contribution to energy and substract it from pmean
179 : : std::array< std::array< tk::real, 3 >, 3 > devH;
180 : :
181 : : // p_mean
182 [ - - ]: 0 : auto pmean = - alpha * elasticEnergy(defgrad, devH);
183 : :
184 : : // Pressure due to shear
185 : 0 : asig[0][0] = -pmean;
186 : 0 : asig[1][1] = -pmean;
187 : 0 : asig[2][2] = -pmean;
188 : :
189 : : // Add deviatoric component of Cauchy stress tensor
190 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
191 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
192 : 0 : asig[i][j] += 2.0*m_mu*alpha*devH[i][j];
193 : : }
194 : :
195 : 0 : return asig;
196 : : }
197 : :
198 : : tk::real
199 : 0 : WilkinsAluminum::soundspeed(
200 : : tk::real arho,
201 : : tk::real apr,
202 : : tk::real alpha,
203 : : std::size_t imat,
204 : : const std::array< std::array< tk::real, 3 >, 3 >& /*defgrad*/ ) const
205 : : // *************************************************************************
206 : : //! Calculate speed of sound from the material density and material pressure
207 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
208 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
209 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
210 : : //! the single-material system, this argument can be left unspecified by
211 : : //! the calling code
212 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
213 : : //! for the single-material system, this argument can be left unspecified
214 : : //! by the calling code
215 : : //! (alpha * sigma_ij * n_j) projected onto the normal vector. Default is 0,
216 : : //! so that for the single-material system, this argument can be left
217 : : //! unspecified by the calling code
218 : : // //! \param[in] defgrad Material inverse deformation gradient tensor
219 : : // //! (g_k) with the first dimension aligned to direction in which
220 : : // //! wave speeds are required. Default is 0, so that for the single-material
221 : : // //! system, this argument can be left unspecified by the calling code
222 : : //! \return Material speed of sound using the WilkinsAluminum EoS
223 : : // *************************************************************************
224 : : {
225 : 0 : tk::real a = 0.0;
226 : :
227 : : // Hydro contribution
228 : 0 : auto al_eff = std::max( 1.0e-14, alpha );
229 : 0 : tk::real rho0 = m_rho0;
230 : 0 : tk::real rho = arho/al_eff;
231 : 0 : a += std::max( 1.0e-15, 6*e2*std::pow(rho/rho0,2.0)/rho0
232 : 0 : + 2*e3*rho/(rho0*rho0) - e5/rho0 );
233 : :
234 : : // Shear contribution
235 : 0 : a += (4.0/3.0) * m_mu / (arho/al_eff);
236 : :
237 : : // Compute square root
238 : 0 : a = std::sqrt(a);
239 : :
240 : : // check sound speed divergence
241 [ - - ]: 0 : if (!std::isfinite(a)) {
242 : 0 : std::cout << "Material-id: " << imat << std::endl;
243 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
244 : 0 : std::cout << "Partial density: " << arho << std::endl;
245 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
246 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
247 : : + std::to_string(a) + ", material volume fraction: " +
248 : : std::to_string(alpha));
249 : : }
250 : :
251 : 0 : return a;
252 : : }
253 : :
254 : : tk::real
255 : 0 : WilkinsAluminum::shearspeed(
256 : : tk::real arho,
257 : : tk::real alpha,
258 : : std::size_t imat ) const
259 : : // *************************************************************************
260 : : //! Calculate speed of sound from the material density and material pressure
261 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
262 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
263 : : //! the single-material system, this argument can be left unspecified by
264 : : //! the calling code
265 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
266 : : //! for the single-material system, this argument can be left unspecified
267 : : //! by the calling code
268 : : //! \return Material shear-wave speed speed using the SmallShearSolid EoS
269 : : // *************************************************************************
270 : : {
271 : : // Approximate shear-wave speed. Ref. Barton, P. T. (2019).
272 : : // An interface-capturing Godunov method for the simulation of compressible
273 : : // solid-fluid problems. Journal of Computational Physics, 390, 25-50.
274 : 0 : auto al_eff = std::max( 1.0e-14, alpha );
275 : 0 : tk::real a = std::sqrt(al_eff*m_mu/arho);
276 : :
277 : : // check shear-wave speed divergence
278 [ - - ]: 0 : if (!std::isfinite(a)) {
279 : 0 : std::cout << "Material-id: " << imat << std::endl;
280 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
281 : 0 : std::cout << "Partial density: " << arho << std::endl;
282 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf shear-wave speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
283 : : + std::to_string(a) + ", material volume fraction: " +
284 : : std::to_string(alpha));
285 : : }
286 : :
287 : 0 : return a;
288 : : }
289 : :
290 : : tk::real
291 : 0 : WilkinsAluminum::totalenergy(
292 : : tk::real arho,
293 : : tk::real u,
294 : : tk::real v,
295 : : tk::real w,
296 : : tk::real,
297 : : tk::real alpha,
298 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
299 : : // *************************************************************************
300 : : //! \brief Calculate material specific total energy from the material
301 : : //! density, momentum and material pressure
302 : : //! \param[in] arho Material partial density
303 : : //! \param[in] u X-velocity
304 : : //! \param[in] v Y-velocity
305 : : //! \param[in] w Z-velocity
306 : : // //! \param[in] apr Material partial pressure
307 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
308 : : //! the single-material system, this argument can be left unspecified by
309 : : //! the calling code
310 : : //! \param[in] defgrad Material inverse deformation gradient tensor
311 : : //! g_k. Default is 0, so that for the single-material system,
312 : : //! this argument can be left unspecified by the calling code
313 : : //! \return Material specific total energy using the WilkinsAluminum EoS
314 : : // *************************************************************************
315 : : {
316 : : // obtain hydro contribution to energy
317 : 0 : tk::real rho0 = m_rho0;
318 : 0 : tk::real rho = arho/alpha;
319 : 0 : tk::real rhoEh = (e1+e2*std::pow(rho/rho0,2.0)+e3*(rho/rho0)
320 : 0 : +e4*std::pow(rho/rho0,-1.0)-e5*std::log(rho/rho0))/rho0
321 : 0 : + 0.5*rho*(u*u + v*v + w*w);
322 : : // obtain elastic contribution to energy
323 : : std::array< std::array< tk::real, 3 >, 3 > devH;
324 [ - - ]: 0 : tk::real rhoEe = elasticEnergy(defgrad, devH);
325 : :
326 : 0 : return alpha*(rhoEh + rhoEe);
327 : : }
328 : :
329 : : tk::real
330 : 0 : WilkinsAluminum::temperature(
331 : : tk::real,
332 : : tk::real,
333 : : tk::real,
334 : : tk::real,
335 : : tk::real,
336 : : tk::real,
337 : : const std::array< std::array< tk::real, 3 >, 3 >& ) const
338 : : // *************************************************************************
339 : : //! \brief Calculate material temperature from the material density, and
340 : : //! material specific total energy
341 : : // //! \param[in] arho Material partial density (alpha_k * rho_k)
342 : : // //! \param[in] u X-velocity
343 : : // //! \param[in] v Y-velocity
344 : : // //! \param[in] w Z-velocity
345 : : // //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
346 : : // //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
347 : : // //! the single-material system, this argument can be left unspecified by
348 : : // //! the calling code
349 : : // //! \param[in] defgrad Material inverse deformation gradient tensor
350 : : // //! (g_k). Default is 0, so that for the single-material system,
351 : : // //! this argument can be left unspecified by the calling code
352 : : //! \return Material temperature using the WilkinsAluminum EoS
353 : : // *************************************************************************
354 : : {
355 : : // Temperature does not directly contribute to energy
356 : : // So we just set a value.
357 : 0 : tk::real t = 300.0;
358 : :
359 : 0 : return t;
360 : : }
361 : :
362 : : tk::real
363 : 0 : WilkinsAluminum::min_eff_pressure(
364 : : tk::real min,
365 : : tk::real,
366 : : tk::real ) const
367 : : // *************************************************************************
368 : : //! Compute the minimum allowed pressure
369 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
370 : : //! \return Minimum pressure allowed by physical constraints
371 : : // *************************************************************************
372 : : {
373 : : // minimum pressure is constrained by zero soundspeed.
374 : 0 : return min;
375 : : }
376 : :
377 : : tk::real
378 : 0 : WilkinsAluminum::elasticEnergy(
379 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad,
380 : : std::array< std::array< tk::real, 3 >, 3 >& devH ) const
381 : : // *************************************************************************
382 : : //! \brief Calculate elastic contribution to material energy from the material
383 : : //! density, and deformation gradient tensor
384 : : //! \param[in] defgrad Material inverse deformation gradient tensor
385 : : //! \param[in/out] devH Deviatoric part of the Hensky tensor
386 : : //! \return Material elastic energy using the WilkinsAluminum EoS
387 : : //! \details This function returns the material elastic energy, and also stores
388 : : //! the elastic shear distortion for further use
389 : : // *************************************************************************
390 : : {
391 : : // Compute deviatoric part of Hencky tensor
392 : 0 : devH = tk::getDevHencky(defgrad);
393 : :
394 : : // Compute elastic energy
395 : 0 : tk::real rhoEe = 0.0;
396 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
397 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
398 : 0 : rhoEe += m_mu*devH[i][j]*devH[i][j];
399 : :
400 : 0 : return rhoEe;
401 : : }
|
