Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/EoS/SmallShearSolid.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 Small shear strain equation of state for solids
9 : : \details This file defines functions for the SmallShearSolid equation of
10 : : state for the compressible flow equations. These functions are
11 : : taken from Plohr, J. N., & Plohr, B. J. (2005). Linearized analysis
12 : : of Richtmyer–Meshkov flow for elastic materials. Journal of Fluid
13 : : Mechanics, 537, 55-89. The SmallShearSolid EOS uses a small-shear
14 : : approximation for the elastic contribution, and a stiffened gas EOS
15 : : for the hydrodynamic contribution of the internal energy.
16 : : */
17 : : // *****************************************************************************
18 : :
19 : : #include <cmath>
20 : : #include <iostream>
21 : : #include "Vector.hpp"
22 : : #include "EoS/SmallShearSolid.hpp"
23 : :
24 : : // Lapacke forward declarations
25 : : extern "C" {
26 : :
27 : : using lapack_int = long;
28 : :
29 : : #define LAPACK_ROW_MAJOR 101
30 : :
31 : : lapack_int LAPACKE_dgeev(int, char, char, lapack_int, double*, lapack_int,
32 : : double*, double*, double*, lapack_int, double*, lapack_int );
33 : :
34 : : }
35 : :
36 : : using inciter::SmallShearSolid;
37 : :
38 : 9 : SmallShearSolid::SmallShearSolid(
39 : : tk::real gamma,
40 : : tk::real pstiff,
41 : : tk::real cv,
42 : 9 : tk::real mu ) :
43 : : m_gamma(gamma),
44 : : m_pstiff(pstiff),
45 : : m_cv(cv),
46 : 9 : m_mu(mu)
47 : : // *************************************************************************
48 : : // Constructor
49 : : //! \param[in] gamma Ratio of specific heats
50 : : //! \param[in] pstiff Stiffness pressure term
51 : : //! \param[in] cv Specific heat at constant volume
52 : : //! \param[in] mu Constant shear modulus
53 : : // *************************************************************************
54 : 9 : { }
55 : :
56 : : tk::real
57 : 20102 : SmallShearSolid::density(
58 : : tk::real pr,
59 : : tk::real temp ) const
60 : : // *************************************************************************
61 : : //! \brief Calculate density from the material pressure and temperature
62 : : //! using the SmallShearSolid equation of state
63 : : //! \param[in] pr Material pressure
64 : : //! \param[in] temp Material temperature
65 : : //! \return Material density calculated using the SmallShearSolid EoS
66 : : // *************************************************************************
67 : : {
68 : 20102 : tk::real g = m_gamma;
69 : 20102 : tk::real p_c = m_pstiff;
70 : 20102 : tk::real c_v = m_cv;
71 : :
72 : 20102 : tk::real rho = (pr + p_c) / ((g-1.0) * c_v * temp);
73 : 20102 : return rho;
74 : : }
75 : :
76 : : tk::real
77 : 147368 : SmallShearSolid::pressure(
78 : : tk::real arho,
79 : : tk::real u,
80 : : tk::real v,
81 : : tk::real w,
82 : : tk::real arhoE,
83 : : tk::real alpha,
84 : : std::size_t imat,
85 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
86 : : // *************************************************************************
87 : : //! \brief Calculate pressure from the material density, momentum, total energy
88 : : //! and the inverse deformation gradient tensor using the SmallShearSolid
89 : : //! equation of state
90 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
91 : : //! \param[in] u X-velocity
92 : : //! \param[in] v Y-velocity
93 : : //! \param[in] w Z-velocity
94 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
95 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
96 : : //! the single-material system, this argument can be left unspecified by
97 : : //! the calling code
98 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
99 : : //! for the single-material system, this argument can be left unspecified
100 : : //! by the calling code
101 : : //! \param[in] defgrad Material inverse deformation gradient tensor
102 : : //! (g_k). Default is 0, so that for the single-material system,
103 : : //! this argument can be left unspecified by the calling code
104 : : //! \return Material partial pressure (alpha_k * p_k) calculated using the
105 : : //! SmallShearSolid EoS
106 : : // *************************************************************************
107 : : {
108 : : // obtain elastic contribution to energy
109 : : tk::real eps2;
110 [ + - ]: 147368 : auto arhoEe = alpha*elasticEnergy(defgrad, eps2);
111 : : // obtain hydro contribution to energy
112 : 147368 : auto arhoEh = arhoE - arhoEe;
113 : :
114 : : // use stiffened gas eos to get pressure
115 : 147368 : tk::real partpressure = (arhoEh - 0.5 * arho * (u*u + v*v + w*w) -
116 : 147368 : alpha*m_pstiff) * (m_gamma-1.0) - alpha*m_pstiff;
117 : :
118 : : // check partial pressure divergence
119 [ - + ]: 147368 : if (!std::isfinite(partpressure)) {
120 [ - - ][ - - ]: 0 : std::cout << "Material-id: " << imat << std::endl;
[ - - ]
121 [ - - ][ - - ]: 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
[ - - ]
122 [ - - ][ - - ]: 0 : std::cout << "Partial density: " << arho << std::endl;
[ - - ]
123 [ - - ][ - - ]: 0 : std::cout << "Total energy: " << arhoE << std::endl;
[ - - ]
124 [ - - ][ - - ]: 0 : std::cout << "Hydro energy: " << arhoEh << std::endl;
[ - - ]
125 [ - - ][ - - ]: 0 : std::cout << "det(defgrad): " << tk::determinant(defgrad) << std::endl;
[ - - ]
126 [ - - ][ - - ]: 0 : std::cout << "Velocity: " << u << ", " << v << ", " << w
[ - - ][ - - ]
[ - - ][ - - ]
127 [ - - ]: 0 : << std::endl;
128 [ - - ][ - - ]: 0 : 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 : 147368 : return partpressure;
134 : : }
135 : :
136 : : std::array< std::array< tk::real, 3 >, 3 >
137 : 102548 : SmallShearSolid::CauchyStress(
138 : : tk::real,
139 : : tk::real,
140 : : tk::real,
141 : : tk::real,
142 : : tk::real,
143 : : tk::real alpha,
144 : : std::size_t /*imat*/,
145 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
146 : : // *************************************************************************
147 : : //! \brief Calculate the elastic Cauchy stress tensor from the material density,
148 : : //! momentum, total energy, and inverse deformation gradient tensor using the
149 : : //! SmallShearSolid equation of state
150 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
151 : : //! the single-material system, this argument can be left unspecified by
152 : : //! the calling code
153 : : // //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
154 : : // //! for the single-material system, this argument can be left unspecified
155 : : // //! by the calling code
156 : : //! \param[in] defgrad Material inverse deformation gradient tensor (g_k).
157 : : //! \return Material Cauchy stress tensor (alpha_k * sigma_k) calculated using
158 : : //! the SmallShearSolid EoS
159 : : // *************************************************************************
160 : : {
161 : 102548 : std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};
162 : :
163 : : // obtain elastic contribution to energy
164 : : tk::real eps2;
165 [ + - ]: 102548 : elasticEnergy(defgrad, eps2);
166 : :
167 : : // p_mean
168 : 102548 : auto pmean = - alpha * m_mu * eps2;
169 : :
170 : : // Volumetric component of Cauchy stress tensor
171 : 102548 : asig[0][0] = -pmean;
172 : 102548 : asig[1][1] = -pmean;
173 : 102548 : asig[2][2] = -pmean;
174 : :
175 : : // Deviatoric (trace-free) part of volume-preserving left Cauchy-Green tensor
176 [ + - ]: 102548 : auto devbt = tk::getLeftCauchyGreen(defgrad);
177 : 102548 : auto detb = std::pow(tk::determinant(devbt), 1.0/3.0);
178 [ + + ]: 410192 : for (std::size_t i=0; i<3; ++i) {
179 [ + + ]: 1230576 : for (std::size_t j=0; j<3; ++j)
180 : 922932 : devbt[i][j] /= detb;
181 : : }
182 : 102548 : auto trbt = (devbt[0][0]+devbt[1][1]+devbt[2][2])/3.0;
183 : 102548 : devbt[0][0] -= trbt;
184 : 102548 : devbt[1][1] -= trbt;
185 : 102548 : devbt[2][2] -= trbt;
186 : :
187 : : // Add deviatoric component of Cauchy stress tensor
188 [ + + ]: 410192 : for (std::size_t i=0; i<3; ++i) {
189 [ + + ]: 1230576 : for (std::size_t j=0; j<3; ++j)
190 : 922932 : asig[i][j] += m_mu*alpha*devbt[i][j];
191 : : }
192 : :
193 : 205096 : return asig;
194 : : }
195 : :
196 : : tk::real
197 : 524832 : SmallShearSolid::soundspeed(
198 : : tk::real arho,
199 : : tk::real apr,
200 : : tk::real alpha,
201 : : std::size_t imat,
202 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad,
203 : : const std::array< tk::real, 3 >& /*adefgradn*/,
204 : : const std::array< tk::real, 3 >& /*asigman*/ ) 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 : : // //! \param[in] adefgradn Material inverse deformation gradient tensor in
223 : : // //! direction of vector n (alpha_k * g_ij * n_j). Default is 0, so that for
224 : : // //! the single-material system, this argument can be left unspecified by the
225 : : // //! calling code
226 : : // //! \param[in] asigman Material traction vector in normal direction
227 : : // //! (alpha * sigma_ij * n_j ). Default is 0, so that for the single-material
228 : : // //! system, this argument can be left unspecified by the calling code
229 : : //! \return Material speed of sound using the SmallShearSolid EoS
230 : : // *************************************************************************
231 : : {
232 : : // deformation gradient
233 : 524832 : tk::real g__11 = defgrad[0][0];
234 : 524832 : tk::real g__12 = defgrad[0][1];
235 : 524832 : tk::real g__13 = defgrad[0][2];
236 : 524832 : tk::real g__21 = defgrad[1][0];
237 : 524832 : tk::real g__22 = defgrad[1][1];
238 : 524832 : tk::real g__23 = defgrad[1][2];
239 : 524832 : tk::real g__31 = defgrad[2][0];
240 : 524832 : tk::real g__32 = defgrad[2][1];
241 : 524832 : tk::real g__33 = defgrad[2][2];
242 : :
243 : : std::array< std::array< tk::real, 3 >, 3> dsigdg;
244 : :
245 : : // d(sigma_11)/d(g_11)
246 : 524832 : dsigdg[0][0] = ( (((g__13 * g__13 * g__22 - 3 * g__12 * g__13 * g__23
247 : 524832 : - 2 * g__22 * (g__12 * g__12 + g__22 * g__22 + g__23 * g__23)) * g__33
248 : 524832 : - (-2 * g__13 * g__13 * g__23 - 3 * g__12 * g__13 * g__22
249 : 524832 : + g__23 * (g__12 * g__12 - 2 * g__22 * g__22 - 2 * g__23 * g__23)) * g__32)
250 : 524832 : * g__31 * g__31) + (( ((3 * g__12 * g__23 - 2 * g__13 * g__22) * g__33
251 : 524832 : * g__33) + (g__32 * (g__12 * g__22 - g__13 * g__23) * g__33) - (3 * g__22
252 : 524832 : * (g__22 * g__22 + g__23 * g__23 + g__32 * g__32) * g__13) + 0.3e1 * g__12
253 : 524832 : * g__23 * ( (g__22 * g__22) + (g__23 * g__23) + 0.2e1 / 0.3e1 * g__32
254 : 524832 : * g__32)) * g__11 + 0.3e1 * (( (g__12 * g__13) + 0.4e1 / 0.3e1 * g__22
255 : 524832 : * g__23) * (g__33 * g__33) + g__32 * ( (g__12 * g__12) - (g__13 * g__13)
256 : 524832 : + 0.4e1 / 0.3e1 * g__22 * g__22 - 0.4e1 / 0.3e1 * g__23 * g__23) * g__33
257 : 524832 : + (g__13 * g__13 * g__22 * g__23) + (g__12 * (g__22 * g__22 - g__23 * g__23
258 : 524832 : - g__32 * g__32) * g__13) - g__22 * g__23 * ( (g__12 * g__12)
259 : 524832 : + 0.4e1 / 0.3e1 * g__32 * g__32)) * g__21) * g__31 + (g__33 * g__22
260 : 524832 : - g__23 * g__32) * (g__22 * g__22 + g__23 * g__23 + g__32 * g__32 + g__33
261 : 524832 : * g__33) * g__11 * g__11 - 0.3e1 * ( ( std::pow( g__33, 3) * g__12)
262 : 524832 : - (g__33 * g__33 * g__13 * g__32) + (- (g__13 * g__22 * g__23) / 0.3e1
263 : 524832 : + 0.2e1 / 0.3e1 * g__12 * ( (g__22 * g__22) + 0.3e1 / 0.2e1 * g__23
264 : 524832 : * g__23 + 0.3e1 / 0.2e1 * g__32 * g__32)) * g__33 + (((-3 * g__22
265 : 524832 : * g__22 - 2 * g__23 * g__23 - 3 * g__32 * g__32) * g__13 + g__12 * g__22
266 : 1049664 : * g__23) * g__32) / 0.3e1) * g__21 * g__11 + (- (14 * g__12 * g__12 * g__22)
267 : 524832 : - 0.2e1 * g__21 * g__21 * g__22 - (14 * std::pow( g__22, 3)))
268 : 524832 : * std::pow( g__33, 3) + 0.14e2 * ( (2 * g__12 * g__13 * g__22)
269 : 524832 : + g__23 * ( (g__12 * g__12) + g__21 * g__21 / 0.7e1
270 : 524832 : + (3 * g__22 * g__22))) * g__32 * (g__33 * g__33)
271 : 524832 : + (-0.2e1 * g__22 * (g__21 * g__21 + (7 * g__22 * g__22)
272 : 524832 : + (7 * g__32 * g__32)) * (g__13 * g__13) + 0.3e1 * g__12
273 : 524832 : * g__23 * (g__21 * g__21 + 0.28e2 / 0.3e1 * g__22 * g__22
274 : 524832 : - 0.28e2 / 0.3e1 * g__32 * g__32) * g__13 + g__22
275 : 524832 : * ((g__21 * g__21 - (14 * g__23 * g__23)) * (g__12 * g__12)
276 : 524832 : - 0.2e1 * (g__32 * g__32) * (g__21 * g__21 + (21 * g__23 * g__23))))
277 : 524832 : * g__33 + 0.2e1 * g__32 * (- g__23 * (g__21 * g__21 - (14 * g__22
278 : 524832 : * g__22) - (14 * g__32 * g__32)) * (g__13 * g__13) / 0.2e1
279 : 524832 : - 0.3e1 / 0.2e1 * g__22 * g__12 * (g__21 * g__21 + 0.28e2 / 0.3e1
280 : 524832 : * g__23 * g__23) * g__13 + g__23 * (g__21 * g__21
281 : 524832 : + (7 * g__23 * g__23)) * (g__12 * g__12 + g__32 * g__32)))
282 : 524832 : * m_mu * std::pow(std::pow( g__33 * (g__11 * g__22 - g__12 * g__21)
283 : 524832 : + (g__12 * g__23 * g__31) + g__13 * (g__21 * g__32 - (g__22 * g__31))
284 : 524832 : - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ( ((g__12 * g__23
285 : 524832 : - g__13 * g__22) * g__31) + (g__33 * g__22 - g__23 * g__32) * g__11
286 : 524832 : - g__21 * (g__33 * g__12 - g__13 * g__32)) / 0.9e1;
287 : :
288 : : // d(sigma_11)/d(g_21)
289 : 524832 : dsigdg[0][1] = 0.2e1 / 0.9e1 * m_mu * (((-g__33 * g__12/0.2e1 - g__13 * g__32)
290 : 524832 : * g__31 * g__31 + 0.3e1 / 0.2e1 * g__11 * (-g__12 * g__13 + g__32 * g__33)
291 : 524832 : * g__31 + g__12 * (g__11 * g__11 + 0.7e1 * g__12 * g__12 + 0.7e1 * g__32
292 : 524832 : * g__32) * g__33 + g__13 * g__32 * (g__11 * g__11 - 0.14e2 * g__12 * g__12
293 : 524832 : - 0.14e2 * g__32 * g__32) / 0.2e1) * g__23 * g__23 + (0.3e1 / 0.2e1 * g__22
294 : 524832 : * (-g__12 * g__32 + g__13 * g__33) * g__31 * g__31 + ((g__33 * g__33 * g__12
295 : 524832 : + g__33 * g__13 * g__32 / 0.2e1 + 0.3e1 / 0.2e1 * g__12 * (g__12 * g__12
296 : 524832 : + g__13 * g__13 + g__32 * g__32)) * g__21 - 0.3e1 / 0.2e1 * g__11 * g__22
297 : 524832 : * (g__12 * g__12 - g__13 * g__13 - g__32 * g__32 + g__33 * g__33)) * g__31
298 : 524832 : - 0.3e1 / 0.2e1 * g__11 * (g__33 * g__33 * g__32 + g__12 * g__13
299 : 524832 : * g__33 / 0.3e1 + g__32 * (g__12 * g__12 + 0.2e1 / 0.3e1 * g__13 * g__13
300 : 524832 : + g__32 * g__32)) * g__21 - 0.14e2 * (g__33 * g__33 * g__12 * g__32
301 : 524832 : + (0.3e1 / 0.28e2 * g__11 * g__11 * g__13 + g__12 * g__12 * g__13 - g__13
302 : 524832 : * g__32 * g__32) * g__33 - 0.3e1 / 0.28e2 * (g__11 * g__11 + 0.28e2 / 0.3e1
303 : 524832 : * g__13 * g__13) * g__12 * g__32) * g__22) * g__23
304 : 524832 : + (g__12 * (g__12 * g__12 + g__13 * g__13 + g__22 * g__22) * g__33
305 : 524832 : - g__13 * (g__12 * g__12 + g__13 * g__13 - g__22 * g__22 / 0.2e1) * g__32)
306 : 524832 : * g__31 * g__31 + (-0.3e1 / 0.2e1 * g__22 * (g__33 * g__33 * g__13 + g__33
307 : 524832 : * g__12 * g__32 / 0.3e1 + g__13 * (g__12 * g__12 + g__13 * g__13
308 : 524832 : + 0.2e1 / 0.3e1 * g__32 * g__32)) * g__21 - 0.2e1 * g__11 * (g__33 * g__33
309 : 524832 : * g__12 * g__13 + g__32 * (g__12 * g__12 - g__13 * g__13 + 0.3e1 / 0.4e1
310 : 524832 : * g__22 * g__22) * g__33 - 0.3e1 / 0.4e1 * g__13 * g__12 * (g__22 * g__22
311 : 524832 : + 0.4e1 / 0.3e1 * g__32 * g__32))) * g__31 - (g__33 * g__12 - g__13 * g__32)
312 : 524832 : * (g__12 * g__12 + g__13 * g__13 + g__32 * g__32 + g__33 * g__33) * g__21
313 : 524832 : * g__21 / 0.2e1 + 0.3e1 / 0.2e1 * g__22 * g__11 * (std::pow(g__33, 0.3e1)
314 : 524832 : + (0.2e1 / 0.3e1 * g__12 * g__12 + g__13 * g__13 + g__32 * g__32) * g__33
315 : 524832 : + g__12 * g__13 * g__32 / 0.3e1) * g__21 + g__12 * (g__11 * g__11 + 0.7e1
316 : 524832 : * g__12 * g__12 + 0.7e1 * g__22 * g__22) * std::pow(g__33, 0.3e1) - g__13
317 : 524832 : * g__32 * (g__11 * g__11 + 0.21e2 * g__12 * g__12 + 0.7e1 * g__22 * g__22)
318 : 524832 : * g__33 * g__33 - ((g__11 * g__11 - 0.14e2 * g__13 * g__13) * g__22 * g__22
319 : 524832 : - 0.2e1 * g__32 * g__32 * (g__11 * g__11 + 0.21e2 * g__13 * g__13)) * g__12
320 : 524832 : * g__33 / 0.2e1 - g__13 * g__32 * (g__22 * g__22 + g__32 * g__32) * (g__11
321 : 524832 : * g__11 + 0.7e1 * g__13 * g__13)) * std::pow(std::pow(g__33 * (g__11 * g__22
322 : 524832 : - g__12 * g__21) + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22
323 : 524832 : * g__31) - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ((-g__11
324 : 524832 : * g__32 + g__12 * g__31) * g__23 - g__13 * g__22 * g__31
325 : 524832 : + (-g__33 * g__12 + g__13 * g__32) * g__21 + g__22 * g__33 * g__11);
326 : :
327 : : // d(sigma_11)/d(g_31)
328 : 524832 : dsigdg[0][2] = -0.2e1 / 0.9e1 * m_mu * ((-g__12 * g__31 * g__31 / 0.2e1
329 : 524832 : + 0.3e1 / 0.2e1 * g__11 * g__31 * g__32 + g__12 * (g__11 * g__11
330 : 524832 : + 0.7e1 * g__12 * g__12 + 0.7e1 * g__32 * g__32)) * std::pow(g__23, 0.3e1)
331 : 524832 : + (((-0.3e1 / 0.2e1 * g__32 * g__11 + g__12 * g__31) * g__33
332 : 524832 : - 0.2e1 * (g__11 * g__12 + 0.3e1 / 0.4e1 * g__31 * g__32) * g__13) * g__21
333 : 524832 : - 0.3e1 / 0.2e1 * ((g__11 * g__31 + 0.28e2 / 0.3e1 * g__12 * g__32) * g__33
334 : 524832 : + 0.2e1 / 0.3e1 * g__13 * (g__11 * g__11 + 0.21e2 * g__12 * g__12
335 : 524832 : - g__31 * g__31 / 0.2e1 + 0.7e1 * g__32 * g__32)) * g__22) * g__23 * g__23
336 : 524832 : + ((-g__33 * g__33 * g__12 / 0.2e1 + 0.3e1 / 0.2e1 * g__33 * g__13 * g__32
337 : 524832 : + g__12 * (g__12 * g__12 + g__13 * g__13 + g__32 * g__32)) * g__21 * g__21
338 : 524832 : + 0.3e1 / 0.2e1 * (g__33 * g__33 * g__11 + g__33 * g__13 * g__31 / 0.3e1
339 : 524832 : - g__12 * g__31 * g__32 / 0.3e1 - 0.4e1 / 0.3e1 * g__11 * (g__12 * g__12
340 : 524832 : - g__13 * g__13 + 0.3e1 / 0.4e1 * g__32 * g__32)) * g__22 * g__21
341 : 524832 : + g__12 * (g__11 * g__11 + 0.7e1 * g__12 * g__12 + 0.7e1 * g__22 * g__22)
342 : 524832 : * g__33 * g__33 - 0.3e1 / 0.2e1 * g__13 * (g__11 * g__12 * g__31 / 0.3e1
343 : 524832 : + g__32 * (g__11 * g__11 + 0.28e2 / 0.3e1 * g__12 * g__12 - 0.28e2 / 0.3e1
344 : 524832 : * g__22 * g__22)) * g__33 - g__12 * (g__12 * g__12 + g__13 * g__13
345 : 524832 : + g__22 * g__22) * g__31 * g__31 / 0.2e1 + g__32 * g__11 * (g__12 * g__12
346 : 524832 : + 0.3e1 / 0.2e1 * g__13 * g__13 + 0.3e1 / 0.2e1 * g__22 * g__22) * g__31
347 : 524832 : + g__12 * ((g__22 * g__22 - g__32 * g__32 / 0.2e1) * g__11 * g__11
348 : 524832 : + (0.21e2 * g__22 * g__22 + 0.7e1 * g__32 * g__32) * g__13 * g__13))
349 : 524832 : * g__23 - g__22 * (g__33 * g__33 * g__13 + 0.3e1 / 0.2e1 * g__33
350 : 524832 : * g__12 * g__32 + g__13 * (g__12 * g__12 + g__13 * g__13 - g__32
351 : 524832 : * g__32 / 0.2e1)) * g__21 * g__21 + (-0.3e1 / 0.2e1 * g__33 * g__33
352 : 524832 : * g__11 * g__12 * g__13 + (0.3e1 / 0.2e1 * g__12 * (g__12 * g__12
353 : 524832 : + g__13 * g__13 + g__22 * g__22) * g__31 - 0.3e1 / 0.2e1 * g__32
354 : 524832 : * g__11 * (g__12 * g__12 - g__13 * g__13 - g__22 * g__22)) * g__33
355 : 524832 : + 0.2e1 * g__13 * (-0.3e1 / 0.4e1 * (g__12 * g__12 + g__13 * g__13
356 : 524832 : + 0.2e1 / 0.3e1 * g__22 * g__22) * g__32 * g__31 + g__11 * g__12
357 : 524832 : * (g__22 * g__22 + 0.3e1 / 0.4e1 * g__32 * g__32))) * g__21 + g__22
358 : 524832 : * (g__13 * (g__11 * g__11 - 0.14e2 * g__12 * g__12 - 0.14e2 * g__22 * g__22)
359 : 524832 : * g__33 * g__33 + (-0.3e1 * g__11 * (g__12 * g__12 + 0.2e1 / 0.3e1 * g__13
360 : 524832 : * g__13 + g__22 * g__22) * g__31 + 0.3e1 * (g__11 * g__11 + 0.28e2 / 0.3e1
361 : 524832 : * g__13 * g__13) * g__12 * g__32) * g__33 - 0.2e1 * g__13 * ((-g__12
362 : 524832 : * g__12 / 0.2e1 - g__13 * g__13 / 0.2e1 - g__22 * g__22 / 0.2e1)
363 : 524832 : * g__31 * g__31 - g__11 * g__12 * g__31 * g__32 / 0.2e1 + (g__22 * g__22
364 : 524832 : + g__32 * g__32) * (g__11 * g__11 + 0.7e1 * g__13 * g__13))) / 0.2e1)
365 : 524832 : * std::pow(std::pow(g__33 * (g__11 * g__22 - g__12 * g__21) + g__12
366 : 524832 : * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22 * g__31) - g__11
367 : 524832 : * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ((-g__32 * g__11 + g__12 * g__31)
368 : 524832 : * g__23 + (-g__33 * g__12 + g__13 * g__32) * g__21 + (g__11 * g__33 - g__13
369 : 524832 : * g__31) * g__22);
370 : :
371 : : // d(sigma_21)/d(g_11)
372 : 524832 : dsigdg[1][0] = ( ((g__33 * (3 * g__13 * g__13 + 4 * g__23 * g__23) * g__21
373 : 524832 : + ((-4 * g__13 * g__13 - 4 * g__23 * g__23) * g__31 + g__33 * g__11
374 : 524832 : * g__13) * g__23) * g__32 * g__32) + ((((g__13 * g__23 * g__23
375 : 524832 : - 6 * g__33 * g__33 * g__13) * g__21 - g__23 * (-7 * g__33 * g__13
376 : 524832 : * g__31 + g__11 * (g__23 * g__23 + g__33 * g__33))) * g__12 - ((g__13
377 : 524832 : * g__13 * g__23 + 8 * g__23 * g__33 * g__33) * g__21 - g__33 * (g__13
378 : 524832 : * g__13 + 8 * g__23 * g__23) * g__31 + g__13 * g__11 * (g__33 - g__23)
379 : 524832 : * (g__33 + g__23)) * g__22) * g__32) + (3 * (g__23 * g__23 + g__33 * g__33)
380 : 524832 : * (g__21 * g__33 - g__23 * g__31) * g__12 * g__12) + (g__22 * (-7 * g__33
381 : 524832 : * g__13 * g__23 * g__21 - g__13 * (-6 * g__23 * g__23 + g__33 * g__33)
382 : 524832 : * g__31 + g__33 * g__11 * (g__23 * g__23 + g__33 * g__33)) * g__12)
383 : 1049664 : + 0.4e1 * (g__22 * g__22) * ( ((g__13 * g__13 * g__33
384 : 524832 : + std::pow( g__33, 3)) * g__21) - (( (g__33 * g__33) + 0.3e1 / 0.4e1
385 : 524832 : * g__13 * g__13) * g__31 + (g__33 * g__11 * g__13) / 0.4e1)
386 : 1049664 : * g__23)) * m_mu * std::pow( std::pow( (g__33 * (g__11 * g__22
387 : 524832 : - g__12 * g__21) + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32
388 : 524832 : - g__22 * g__31) - g__11 * g__23 * g__32), 2), -0.2e1 / 0.3e1)
389 : 524832 : / ((-g__11 * g__23 + g__13 * g__21) * g__32 + (-g__21 * g__33
390 : 524832 : + g__23 * g__31) * g__12 + (g__11 * g__33 - g__13 * g__31)
391 : 524832 : * g__22) / 0.3e1;
392 : :
393 : : // d(sigma_21)/d(g_21)
394 : 524832 : dsigdg[1][1] = -0.4e1 / 0.3e1 * ((g__33 * (g__13 * g__13
395 : 524832 : + 0.3e1 / 0.4e1 * g__23 * g__23) * g__11
396 : 524832 : + ((-0.4e1 * g__13 * g__13 - 0.4e1 * g__23 * g__23)
397 : 524832 : * g__31 + g__33 * g__21 * g__23) * g__13 / 0.4e1) * g__32
398 : 524832 : * g__32 + ((-0.3e1 / 0.2e1 * (g__33 * g__33 - g__13 * g__13 / 0.6e1)
399 : 524832 : * g__23 * g__11 - g__13 * (-0.7e1 * g__33 * g__23 * g__31 + g__21
400 : 524832 : * (g__13 * g__13 + g__33 * g__33)) / 0.4e1) * g__22 - 0.2e1 * (g__13
401 : 524832 : * (g__33 * g__33 + g__23 * g__23 / 0.8e1) * g__11 - g__33 * (g__13
402 : 524832 : * g__13 + g__23 * g__23 / 0.8e1) * g__31 + g__21 * g__23 * (g__33 - g__13)
403 : 524832 : * (g__33 + g__13) / 0.8e1) * g__12) * g__32 + 0.3e1 / 0.4e1
404 : 524832 : * (g__13 * g__13 + g__33 * g__33) * (g__11 * g__33 - g__13 * g__31)
405 : 524832 : * g__22 * g__22 + g__12 * (-0.7e1 * g__33 * g__11 * g__13 * g__23 - g__23
406 : 524832 : * (-0.6e1 * g__13 * g__13 + g__33 * g__33) * g__31 + g__33 * g__21
407 : 524832 : * (g__13 * g__13 + g__33 * g__33)) * g__22 / 0.4e1 + ((g__23 * g__23
408 : 524832 : * g__33 + std::pow(g__33, 0.3e1)) * g__11 - g__13 * ((g__33 * g__33
409 : 524832 : + 0.3e1 / 0.4e1 * g__23 * g__23) * g__31 + g__33 * g__21 * g__23 / 0.4e1))
410 : 524832 : * g__12 * g__12) * m_mu * std::pow(std::pow(g__33 * (g__11 * g__22 - g__12
411 : 524832 : * g__21) + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22 * g__31)
412 : 524832 : - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ((-g__11 * g__23 + g__13
413 : 524832 : * g__21) * g__32 + (g__11 * g__33 - g__13 * g__31) * g__22 - g__12
414 : 524832 : * (g__21 * g__33 - g__23 * g__31));
415 : :
416 : : // d(sigma_21)/d(g_31)
417 : 524832 : dsigdg[1][2] = 0.4e1 / 0.3e1 * m_mu * (((-std::pow(g__13, 0.3e1)
418 : 524832 : - g__33 * g__33 * g__13) * g__21 + 0.3e1 / 0.4e1 * ((g__33 * g__33
419 : 524832 : + 0.4e1 / 0.3e1 * g__13 * g__13) * g__11 + g__33 * g__13 * g__31 / 0.3e1)
420 : 524832 : * g__23) * g__22 * g__22 + ((0.7e1 / 0.4e1 * g__33 * g__13 * g__23 * g__21
421 : 524832 : + g__33 * (g__13 * g__13 - 0.6e1 * g__23 * g__23) * g__11 / 0.4e1
422 : 524832 : - g__13 * g__31 * (g__13 * g__13 + g__23 * g__23) / 0.4e1) * g__32
423 : 524832 : - (-g__23 * (0.8e1 * g__13 * g__13 + g__33 * g__33) * g__21
424 : 524832 : + g__13 * (0.8e1 * g__23 * g__23 + g__33 * g__33) * g__11
425 : 524832 : - g__33 * g__31 * (g__13 - g__23) * (g__13 + g__23)) * g__12 / 0.4e1)
426 : 524832 : * g__22 + 0.3e1 / 0.4e1 * (g__13 * g__13 + g__23 * g__23)
427 : 524832 : * (g__11 * g__23 - g__13 * g__21) * g__32 * g__32 - 0.7e1 / 0.4e1
428 : 524832 : * (-0.6e1 / 0.7e1 * g__33 * (g__13 * g__13 - g__23 * g__23 / 0.6e1)
429 : 524832 : * g__21 + (g__33 * g__11 * g__13 - g__31 * (g__13 * g__13 + g__23
430 : 524832 : * g__23) / 0.7e1) * g__23) * g__12 * g__32 + g__12 * g__12
431 : 524832 : * ((-0.3e1 / 0.4e1 * g__33 * g__33 * g__13 - g__13 * g__23 * g__23)
432 : 524832 : * g__21 + g__23 * (g__11 * (g__23 * g__23 + g__33 * g__33)
433 : 524832 : - g__33 * g__13 * g__31 / 0.4e1))) * std::pow(std::pow(g__33
434 : 524832 : * (g__11 * g__22 - g__12 * g__21) + g__12 * g__23 * g__31 + g__13
435 : 524832 : * (g__21 * g__32 - g__22 * g__31) - g__11 * g__23 * g__32, 0.2e1),
436 : 524832 : -0.2e1 / 0.3e1) / ((-g__11 * g__23 + g__13 * g__21) * g__32
437 : 524832 : + (g__11 * g__33 - g__13 * g__31) * g__22 - g__12 * (g__21
438 : 524832 : * g__33 - g__23 * g__31));
439 : :
440 : : // d(sigma_31)/d(g_11)
441 : 524832 : dsigdg[2][0] = - ((g__32 * (3 * g__12 * g__12 + 4 * g__22 * g__22)
442 : 524832 : * g__21 + g__22 * ((-4 * g__12 * g__12 - 4 * g__22 * g__22) * g__31
443 : 524832 : + g__11 * g__12 * g__32)) * g__33 * g__33 + ((g__12 * (g__22 * g__22
444 : 524832 : - 6 * g__32 * g__32) * g__21 - g__22 * (-7 * g__12 * g__31 * g__32
445 : 524832 : + g__11 * (g__22 * g__22 + g__32 * g__32))) * g__13
446 : 524832 : + (-g__22 * (g__12 * g__12 + 8 * g__32 * g__32) * g__21 + g__32
447 : 524832 : * (g__12 * g__12 + 8 * g__22 * g__22) * g__31 + g__11 * g__12
448 : 524832 : * (g__22 - g__32) * (g__22 + g__32)) * g__23) * g__33
449 : 524832 : + 3 * (g__22 * g__22 + g__32 * g__32) * (g__21 * g__32 - g__22 * g__31)
450 : 524832 : * g__13 * g__13 + g__23 * (-7 * g__12 * g__22 * g__32 * g__21
451 : 524832 : + (6 * g__12 * g__22 * g__22 - g__12 * g__32 * g__32) * g__31
452 : 1049664 : + g__32 * g__11 * (g__22 * g__22 + g__32 * g__32)) * g__13 - g__23
453 : 524832 : * g__23 * ((-4 * g__12 * g__12 * g__32 - 4 * std::pow( g__32, 3))
454 : 524832 : * g__21 + g__22 * (g__31 * (3 * g__12 * g__12 + 4 * g__32 * g__32)
455 : 1049664 : + g__11 * g__12 * g__32))) * m_mu * std::pow( std::pow( (g__33
456 : 524832 : * (g__11 * g__22 - g__12 * g__21) + g__12 * g__23 * g__31 + g__13
457 : 524832 : * (g__21 * g__32 - g__22 * g__31) - g__11 * g__23 * g__32), 2),
458 : 524832 : -0.2e1 / 0.3e1) / (g__33 * (g__11 * g__22 - g__12 * g__21) + g__13
459 : 524832 : * (g__21 * g__32 - g__22 * g__31) - g__23 * (g__32 * g__11 - g__12
460 : 524832 : * g__31)) / 0.3e1;
461 : :
462 : : // d(sigma_31)/d(g_21)
463 : 524832 : dsigdg[2][1] = 0.4e1 / 0.3e1 * (((-std::pow(g__12, 0.3e1)
464 : 524832 : - g__12 * g__22 * g__22) * g__31 + ((g__12 * g__12 + 0.3e1 / 0.4e1
465 : 524832 : * g__22 * g__22) * g__11 + g__12 * g__21 * g__22 / 0.4e1) * g__32)
466 : 524832 : * g__33 * g__33 + ((0.7e1 / 0.4e1 * g__12 * g__22 * g__31 * g__32
467 : 524832 : + g__22 * (g__12 * g__12 - 0.6e1 * g__32 * g__32) * g__11 / 0.4e1
468 : 524832 : - g__12 * g__21 * (g__12 * g__12 + g__32 * g__32) / 0.4e1) * g__23
469 : 524832 : - g__13 * ((-0.8e1 * g__12 * g__12 - g__22 * g__22) * g__32 * g__31
470 : 524832 : + g__12 * (g__22 * g__22 + 0.8e1 * g__32 * g__32) * g__11 - g__21
471 : 524832 : * g__22 * (g__12 - g__32) * (g__12 + g__32)) / 0.4e1) * g__33
472 : 524832 : + 0.3e1 / 0.4e1 * (g__12 * g__12 + g__32 * g__32) * (g__32 * g__11
473 : 524832 : - g__12 * g__31) * g__23 * g__23 - 0.7e1 / 0.4e1 * g__13
474 : 524832 : * (-0.6e1 / 0.7e1 * g__22 * (g__12 * g__12 - g__32 * g__32 / 0.6e1)
475 : 524832 : * g__31 + g__32 * (g__11 * g__12 * g__22 - g__21 * (g__12 * g__12
476 : 524832 : + g__32 * g__32) / 0.7e1)) * g__23 + g__13 * g__13 * ((-0.3e1 / 0.4e1
477 : 524832 : * g__12 * g__22 * g__22 - g__12 * g__32 * g__32) * g__31
478 : 524832 : + (g__11 * (g__22 * g__22 + g__32 * g__32) - g__12 * g__21
479 : 524832 : * g__22 / 0.4e1) * g__32)) * m_mu * std::pow(std::pow(g__33
480 : 524832 : * (g__11 * g__22 - g__12 * g__21) + g__12 * g__23 * g__31 + g__13
481 : 524832 : * (g__21 * g__32 - g__22 * g__31) - g__11 * g__23 * g__32, 0.2e1),
482 : 524832 : -0.2e1 / 0.3e1) / (g__13 * (g__21 * g__32 - g__22 * g__31)
483 : 524832 : + (-g__32 * g__11 + g__12 * g__31) * g__23 + g__33
484 : 524832 : * (g__11 * g__22 - g__12 * g__21));
485 : :
486 : : // d(sigma_31)/d(g_31)
487 : 524832 : dsigdg[2][2] = -(((-0.4e1 / 0.3e1 * g__12 * g__32 * g__32
488 : 524832 : - 0.4e1 / 0.3e1 * std::pow(g__12, 0.3e1)) * g__21 + 0.4e1 / 0.3e1
489 : 524832 : * g__22 * ((g__12 * g__12 + 0.3e1 / 0.4e1 * g__32 * g__32) * g__11
490 : 524832 : + g__12 * g__31 * g__32 / 0.4e1)) * g__23 * g__23 + ((0.7e1 / 0.3e1
491 : 524832 : * g__12 * g__22 * g__32 * g__21 + g__32 * (g__12 * g__12 - 0.6e1
492 : 524832 : * g__22 * g__22) * g__11 / 0.3e1 - g__12 * g__31 * (g__12 * g__12
493 : 524832 : + g__22 * g__22) / 0.3e1) * g__33 - 0.8e1 / 0.3e1 * g__13 * (-g__22
494 : 524832 : * (g__12 * g__12 + g__32 * g__32 / 0.8e1) * g__21 + g__12 * (g__22
495 : 524832 : * g__22 + g__32 * g__32 / 0.8e1) * g__11 - g__31 * g__32 * (g__12
496 : 524832 : - g__22) * (g__12 + g__22) / 0.8e1)) * g__23 + (g__12 * g__12
497 : 524832 : + g__22 * g__22) * (g__11 * g__22 - g__12 * g__21) * g__33 * g__33
498 : 524832 : - 0.7e1 / 0.3e1 * (-0.6e1 / 0.7e1 * g__32 * (g__12 * g__12 - g__22
499 : 524832 : * g__22 / 0.6e1) * g__21 + (g__11 * g__12 * g__32 - g__31 * (g__12
500 : 524832 : * g__12 + g__22 * g__22) / 0.7e1) * g__22) * g__13 * g__33
501 : 524832 : + 0.4e1 / 0.3e1 * (-g__12 * (g__22 * g__22 + 0.3e1 / 0.4e1
502 : 524832 : * g__32 * g__32) * g__21 + g__22 * (g__11 * (g__22 * g__22
503 : 524832 : + g__32 * g__32) - g__12 * g__31 * g__32 / 0.4e1)) * g__13 * g__13)
504 : 524832 : * m_mu * std::pow(std::pow(g__33 * (g__11 * g__22 - g__12 * g__21)
505 : 524832 : + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22 * g__31)
506 : 524832 : - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / (g__13 * (g__21
507 : 524832 : * g__32 - g__22 * g__31) + (-g__32 * g__11 + g__12 * g__31) * g__23
508 : 524832 : + g__33 * (g__11 * g__22 - g__12 * g__21));
509 : :
510 : : // Define amat
511 : : double amat[9];
512 [ + + ]: 2099328 : for (std::size_t i=0; i<3; ++i)
513 [ + + ]: 6297984 : for (std::size_t j=0; j<3; ++j)
514 : : {
515 : 4723488 : amat[i*3+j] = 0.0;
516 [ + + ]: 18893952 : for (std::size_t k=0; k<3; ++k)
517 : : {
518 : 14170464 : amat[i*3+j] -= defgrad[k][j]*dsigdg[i][k];
519 : : }
520 : 4723488 : amat[i*3+j] /= (arho/alpha);
521 : : }
522 : :
523 : : double eig_real[3], eig_imag[3];
524 : : double vl[3], vr[3];
525 : : #ifndef NDEBUG
526 : : lapack_int ierr =
527 : : #endif
528 [ + - ]: 524832 : LAPACKE_dgeev (LAPACK_ROW_MAJOR, 'N', 'N', 3, amat, 3,
529 : : eig_real, eig_imag, vl, 3, vr, 3);
530 [ - + ][ - - ]: 524832 : Assert(ierr==0, "Lapack failed to compute eigenvalues");
[ - - ][ - - ]
531 : :
532 : : // Manually find max
533 : 524832 : tk::real eig_max = eig_real[0];
534 [ + + ]: 1574496 : for (std::size_t i=1; i<3; i++)
535 [ + + ]: 1049664 : if (eig_real[i] > eig_max)
536 : 1112 : eig_max = eig_real[i];
537 : 524832 : tk::real a = std::sqrt(eig_max);
538 : :
539 : : // hydrodynamic contribution
540 : 524832 : auto p_eff = std::max( 1.0e-15, apr+(alpha*m_pstiff) );
541 : 524832 : a += std::sqrt( m_gamma * p_eff / arho );
542 : :
543 : : // check sound speed divergence
544 [ - + ]: 524832 : if (!std::isfinite(a)) {
545 [ - - ][ - - ]: 0 : std::cout << "Material-id: " << imat << std::endl;
[ - - ]
546 [ - - ][ - - ]: 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
[ - - ]
547 [ - - ][ - - ]: 0 : std::cout << "Partial density: " << arho << std::endl;
[ - - ]
548 [ - - ][ - - ]: 0 : std::cout << "Partial pressure: " << apr << std::endl;
[ - - ]
549 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
550 : : + std::to_string(a) + ", material volume fraction: " +
551 : : std::to_string(alpha));
552 : : }
553 : :
554 : 524832 : return a;
555 : : }
556 : :
557 : : tk::real
558 : 569366 : SmallShearSolid::totalenergy(
559 : : tk::real rho,
560 : : tk::real u,
561 : : tk::real v,
562 : : tk::real w,
563 : : tk::real pr,
564 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
565 : : // *************************************************************************
566 : : //! \brief Calculate material specific total energy from the material
567 : : //! density, momentum and material pressure
568 : : //! \param[in] rho Material density
569 : : //! \param[in] u X-velocity
570 : : //! \param[in] v Y-velocity
571 : : //! \param[in] w Z-velocity
572 : : //! \param[in] pr Material pressure
573 : : //! \param[in] defgrad Material inverse deformation gradient tensor
574 : : //! g_k. Default is 0, so that for the single-material system,
575 : : //! this argument can be left unspecified by the calling code
576 : : //! \return Material specific total energy using the SmallShearSolid EoS
577 : : // *************************************************************************
578 : : {
579 : : // obtain hydro contribution to energy
580 : 569366 : tk::real rhoEh = (pr + m_pstiff) / (m_gamma-1.0) + 0.5 * rho *
581 : 569366 : (u*u + v*v + w*w) + m_pstiff;
582 : : // obtain elastic contribution to energy
583 : : tk::real eps2;
584 [ + - ]: 569366 : tk::real rhoEe = elasticEnergy(defgrad, eps2);
585 : :
586 : 569366 : return (rhoEh + rhoEe);
587 : : }
588 : :
589 : : tk::real
590 : 44820 : SmallShearSolid::temperature(
591 : : tk::real arho,
592 : : tk::real u,
593 : : tk::real v,
594 : : tk::real w,
595 : : tk::real arhoE,
596 : : tk::real alpha,
597 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
598 : : // *************************************************************************
599 : : //! \brief Calculate material temperature from the material density, and
600 : : //! material specific total energy
601 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
602 : : //! \param[in] u X-velocity
603 : : //! \param[in] v Y-velocity
604 : : //! \param[in] w Z-velocity
605 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
606 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
607 : : //! the single-material system, this argument can be left unspecified by
608 : : //! the calling code
609 : : //! \param[in] defgrad Material inverse deformation gradient tensor
610 : : //! (g_k). Default is 0, so that for the single-material system,
611 : : //! this argument can be left unspecified by the calling code
612 : : //! \return Material temperature using the SmallShearSolid EoS
613 : : // *************************************************************************
614 : : {
615 : : // obtain elastic contribution to energy
616 : : tk::real eps2;
617 [ + - ]: 44820 : auto arhoEe = alpha*elasticEnergy(defgrad, eps2);
618 : : // obtain hydro contribution to energy
619 : 44820 : auto arhoEh = arhoE - arhoEe;
620 : :
621 : 44820 : tk::real t = (arhoEh - 0.5 * arho * (u*u + v*v + w*w) - alpha*m_pstiff)
622 : 44820 : / (arho*m_cv);
623 : 44820 : return t;
624 : : }
625 : :
626 : : tk::real
627 : 2393828 : SmallShearSolid::min_eff_pressure(
628 : : tk::real min,
629 : : tk::real,
630 : : tk::real ) const
631 : : // *************************************************************************
632 : : //! Compute the minimum allowed pressure
633 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
634 : : //! \return Minimum pressure allowed by physical constraints
635 : : // *************************************************************************
636 : : {
637 : : // minimum pressure is constrained by zero soundspeed.
638 : 2393828 : return (min - m_pstiff);
639 : : }
640 : :
641 : : tk::real
642 : 864102 : SmallShearSolid::elasticEnergy(
643 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad,
644 : : tk::real& eps2 ) const
645 : : // *************************************************************************
646 : : //! \brief Calculate elastic contribution to material energy from the material
647 : : //! density, and deformation gradient tensor
648 : : //! \param[in] defgrad Material inverse deformation gradient tensor
649 : : //! \param[in/out] eps2 Elastic shear distortion
650 : : //! \return Material elastic energy using the SmallShearSolid EoS
651 : : //! \details This function returns the material elastic energy, and also stores
652 : : //! the elastic shear distortion for further use
653 : : // *************************************************************************
654 : : {
655 : : // compute Right Cauchy-Green strain tensor
656 [ + - ]: 864102 : auto Ct = tk::getRightCauchyGreen(defgrad);
657 : 864102 : auto detC = std::pow(tk::determinant(Ct), 1.0/3.0);
658 [ + + ]: 3456408 : for (std::size_t i=0; i<3; ++i) {
659 [ + + ]: 10369224 : for (std::size_t j=0; j<3; ++j)
660 : 7776918 : Ct[i][j] /= detC;
661 : : }
662 : :
663 : : // compute elastic shear distortion
664 : 864102 : eps2 = 0.5 * (Ct[0][0]+Ct[1][1]+Ct[2][2] - 3.0);
665 : :
666 : : // compute elastic energy
667 : 864102 : auto rhoEe = m_mu * eps2;
668 : :
669 : 864102 : return rhoEe;
670 : : }
|
