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 : 0 : SmallShearSolid::SmallShearSolid(
39 : : tk::real gamma,
40 : : tk::real pstiff,
41 : : tk::real cv,
42 : 0 : tk::real mu ) :
43 : : m_gamma(gamma),
44 : : m_pstiff(pstiff),
45 : : m_cv(cv),
46 : 0 : 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 : 0 : { }
55 : :
56 : : void
57 : 0 : SmallShearSolid::setRho0( tk::real rho0 )
58 : : // *************************************************************************
59 : : // Set rho0 EOS parameter; i.e. the initial density
60 : : //! \param[in] rho0 Initial material density that needs to be stored
61 : : // *************************************************************************
62 : : {
63 : 0 : m_rho0 = rho0;
64 : 0 : }
65 : :
66 : : tk::real
67 : 0 : SmallShearSolid::density(
68 : : tk::real pr,
69 : : tk::real temp ) const
70 : : // *************************************************************************
71 : : //! \brief Calculate density from the material pressure and temperature
72 : : //! using the SmallShearSolid equation of state
73 : : //! \param[in] pr Material pressure
74 : : //! \param[in] temp Material temperature
75 : : //! \return Material density calculated using the SmallShearSolid EoS
76 : : // *************************************************************************
77 : : {
78 : 0 : tk::real g = m_gamma;
79 : 0 : tk::real p_c = m_pstiff;
80 : 0 : tk::real c_v = m_cv;
81 : :
82 : 0 : tk::real rho = (pr + p_c) / ((g-1.0) * c_v * temp);
83 : 0 : return rho;
84 : : }
85 : :
86 : : tk::real
87 : 0 : SmallShearSolid::pressure(
88 : : tk::real arho,
89 : : tk::real u,
90 : : tk::real v,
91 : : tk::real w,
92 : : tk::real arhoE,
93 : : tk::real alpha,
94 : : std::size_t imat,
95 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
96 : : // *************************************************************************
97 : : //! \brief Calculate pressure from the material density, momentum, total energy
98 : : //! and the inverse deformation gradient tensor using the SmallShearSolid
99 : : //! equation of state
100 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
101 : : //! \param[in] u X-velocity
102 : : //! \param[in] v Y-velocity
103 : : //! \param[in] w Z-velocity
104 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
105 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
106 : : //! the single-material system, this argument can be left unspecified by
107 : : //! the calling code
108 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
109 : : //! for the single-material system, this argument can be left unspecified
110 : : //! by the calling code
111 : : //! \param[in] defgrad Material inverse deformation gradient tensor
112 : : //! (g_k). Default is 0, so that for the single-material system,
113 : : //! this argument can be left unspecified by the calling code
114 : : //! \return Material partial pressure (alpha_k * p_k) calculated using the
115 : : //! SmallShearSolid EoS
116 : : // *************************************************************************
117 : : {
118 : : // obtain elastic contribution to energy
119 : : tk::real eps2;
120 [ - - ]: 0 : auto arhoEe = alpha*elasticEnergy(defgrad, eps2);
121 : : // obtain hydro contribution to energy
122 : 0 : auto arhoEh = arhoE - arhoEe;
123 : :
124 : : // use stiffened gas eos to get pressure
125 : 0 : tk::real partpressure = (arhoEh - 0.5 * arho * (u*u + v*v + w*w) -
126 : 0 : alpha*m_pstiff) * (m_gamma-1.0) - alpha*m_pstiff;
127 : :
128 : : // check partial pressure divergence
129 [ - - ]: 0 : if (!std::isfinite(partpressure)) {
130 [ - - ][ - - ]: 0 : std::cout << "Material-id: " << imat << std::endl;
[ - - ]
131 [ - - ][ - - ]: 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
[ - - ]
132 [ - - ][ - - ]: 0 : std::cout << "Partial density: " << arho << std::endl;
[ - - ]
133 [ - - ][ - - ]: 0 : std::cout << "Total energy: " << arhoE << std::endl;
[ - - ]
134 [ - - ][ - - ]: 0 : std::cout << "Hydro energy: " << arhoEh << std::endl;
[ - - ]
135 [ - - ][ - - ]: 0 : std::cout << "det(defgrad): " << tk::determinant(defgrad) << std::endl;
[ - - ]
136 [ - - ][ - - ]: 0 : std::cout << "Velocity: " << u << ", " << v << ", " << w
[ - - ][ - - ]
[ - - ][ - - ]
137 [ - - ]: 0 : << std::endl;
138 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) +
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
139 : : " has nan/inf partial pressure: " + std::to_string(partpressure) +
140 : : ", material volume fraction: " + std::to_string(alpha));
141 : : }
142 : :
143 : 0 : return partpressure;
144 : : }
145 : :
146 : : std::array< std::array< tk::real, 3 >, 3 >
147 : 0 : SmallShearSolid::CauchyStress(
148 : : tk::real,
149 : : tk::real,
150 : : tk::real,
151 : : tk::real,
152 : : tk::real,
153 : : tk::real alpha,
154 : : std::size_t /*imat*/,
155 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
156 : : // *************************************************************************
157 : : //! \brief Calculate the elastic Cauchy stress tensor from the material density,
158 : : //! momentum, total energy, and inverse deformation gradient tensor using the
159 : : //! SmallShearSolid equation of state
160 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
161 : : //! the single-material system, this argument can be left unspecified by
162 : : //! the calling code
163 : : // //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
164 : : // //! for the single-material system, this argument can be left unspecified
165 : : // //! by the calling code
166 : : //! \param[in] defgrad Material inverse deformation gradient tensor (g_k).
167 : : //! \return Material Cauchy stress tensor (alpha_k * sigma_k) calculated using
168 : : //! the SmallShearSolid EoS
169 : : // *************************************************************************
170 : : {
171 : 0 : std::array< std::array< tk::real, 3 >, 3 > asig{{{0,0,0}, {0,0,0}, {0,0,0}}};
172 : :
173 : : // obtain elastic contribution to energy
174 : : tk::real eps2;
175 [ - - ]: 0 : elasticEnergy(defgrad, eps2);
176 : :
177 : : // p_mean
178 : 0 : auto pmean = - alpha * m_mu * eps2;
179 : :
180 : : // Volumetric component of Cauchy stress tensor
181 : 0 : asig[0][0] = -pmean;
182 : 0 : asig[1][1] = -pmean;
183 : 0 : asig[2][2] = -pmean;
184 : :
185 : : // Deviatoric (trace-free) part of volume-preserving left Cauchy-Green tensor
186 [ - - ]: 0 : auto devbt = tk::getLeftCauchyGreen(defgrad);
187 : 0 : auto detb = std::pow(tk::determinant(devbt), 1.0/3.0);
188 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
189 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
190 : 0 : devbt[i][j] /= detb;
191 : : }
192 : 0 : auto trbt = (devbt[0][0]+devbt[1][1]+devbt[2][2])/3.0;
193 : 0 : devbt[0][0] -= trbt;
194 : 0 : devbt[1][1] -= trbt;
195 : 0 : devbt[2][2] -= trbt;
196 : :
197 : : // Add deviatoric component of Cauchy stress tensor
198 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
199 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
200 : 0 : asig[i][j] += m_mu*alpha*devbt[i][j];
201 : : }
202 : :
203 : 0 : return asig;
204 : : }
205 : :
206 : : tk::real
207 : 0 : SmallShearSolid::soundspeed(
208 : : tk::real arho,
209 : : tk::real apr,
210 : : tk::real alpha,
211 : : std::size_t imat,
212 : : const std::array< std::array< tk::real, 3 >, 3 >& /*defgrad*/,
213 : : const std::array< tk::real, 3 >& /*asigman*/ ) const
214 : : // *************************************************************************
215 : : //! Calculate speed of sound from the material density and material pressure
216 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
217 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
218 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
219 : : //! the single-material system, this argument can be left unspecified by
220 : : //! the calling code
221 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
222 : : //! for the single-material system, this argument can be left unspecified
223 : : //! by the calling code
224 : : //! (alpha * sigma_ij * n_j) projected onto the normal vector. Default is 0,
225 : : //! so that for the single-material system, this argument can be left
226 : : //! unspecified by the calling code
227 : : // //! \param[in] defgrad Material inverse deformation gradient tensor
228 : : // //! (g_k) with the first dimension aligned to direction in which
229 : : // //! wave speeds are required. Default is 0, so that for the single-material
230 : : // //! system, this argument can be left unspecified by the calling code
231 : : // //! \param[in] asigman Material traction vector in normal direction
232 : : // //! (alpha * sigma_ij * n_j ). Default is 0, so that for the single-material
233 : : // //! system, this argument can be left unspecified by the calling code
234 : : //! \return Material speed of sound using the SmallShearSolid EoS
235 : : // *************************************************************************
236 : : {
237 : : /* Rigorous approach: Eigenvalues of full elastic tensor
238 : :
239 : : // deformation gradient
240 : : tk::real g__11 = defgrad[0][0];
241 : : tk::real g__12 = defgrad[0][1];
242 : : tk::real g__13 = defgrad[0][2];
243 : : tk::real g__21 = defgrad[1][0];
244 : : tk::real g__22 = defgrad[1][1];
245 : : tk::real g__23 = defgrad[1][2];
246 : : tk::real g__31 = defgrad[2][0];
247 : : tk::real g__32 = defgrad[2][1];
248 : : tk::real g__33 = defgrad[2][2];
249 : :
250 : : std::array< std::array< tk::real, 3 >, 3> dsigdg;
251 : :
252 : : // d(sigma_11)/d(g_11)
253 : : dsigdg[0][0] = ( (((g__13 * g__13 * g__22 - 3 * g__12 * g__13 * g__23
254 : : - 2 * g__22 * (g__12 * g__12 + g__22 * g__22 + g__23 * g__23)) * g__33
255 : : - (-2 * g__13 * g__13 * g__23 - 3 * g__12 * g__13 * g__22
256 : : + g__23 * (g__12 * g__12 - 2 * g__22 * g__22 - 2 * g__23 * g__23)) * g__32)
257 : : * g__31 * g__31) + (( ((3 * g__12 * g__23 - 2 * g__13 * g__22) * g__33
258 : : * g__33) + (g__32 * (g__12 * g__22 - g__13 * g__23) * g__33) - (3 * g__22
259 : : * (g__22 * g__22 + g__23 * g__23 + g__32 * g__32) * g__13) + 0.3e1 * g__12
260 : : * g__23 * ( (g__22 * g__22) + (g__23 * g__23) + 0.2e1 / 0.3e1 * g__32
261 : : * g__32)) * g__11 + 0.3e1 * (( (g__12 * g__13) + 0.4e1 / 0.3e1 * g__22
262 : : * g__23) * (g__33 * g__33) + g__32 * ( (g__12 * g__12) - (g__13 * g__13)
263 : : + 0.4e1 / 0.3e1 * g__22 * g__22 - 0.4e1 / 0.3e1 * g__23 * g__23) * g__33
264 : : + (g__13 * g__13 * g__22 * g__23) + (g__12 * (g__22 * g__22 - g__23 * g__23
265 : : - g__32 * g__32) * g__13) - g__22 * g__23 * ( (g__12 * g__12)
266 : : + 0.4e1 / 0.3e1 * g__32 * g__32)) * g__21) * g__31 + (g__33 * g__22
267 : : - g__23 * g__32) * (g__22 * g__22 + g__23 * g__23 + g__32 * g__32 + g__33
268 : : * g__33) * g__11 * g__11 - 0.3e1 * ( ( std::pow( g__33, 3) * g__12)
269 : : - (g__33 * g__33 * g__13 * g__32) + (- (g__13 * g__22 * g__23) / 0.3e1
270 : : + 0.2e1 / 0.3e1 * g__12 * ( (g__22 * g__22) + 0.3e1 / 0.2e1 * g__23
271 : : * g__23 + 0.3e1 / 0.2e1 * g__32 * g__32)) * g__33 + (((-3 * g__22
272 : : * g__22 - 2 * g__23 * g__23 - 3 * g__32 * g__32) * g__13 + g__12 * g__22
273 : : * g__23) * g__32) / 0.3e1) * g__21 * g__11 + (- (14 * g__12 * g__12 * g__22)
274 : : - 0.2e1 * g__21 * g__21 * g__22 - (14 * std::pow( g__22, 3)))
275 : : * std::pow( g__33, 3) + 0.14e2 * ( (2 * g__12 * g__13 * g__22)
276 : : + g__23 * ( (g__12 * g__12) + g__21 * g__21 / 0.7e1
277 : : + (3 * g__22 * g__22))) * g__32 * (g__33 * g__33)
278 : : + (-0.2e1 * g__22 * (g__21 * g__21 + (7 * g__22 * g__22)
279 : : + (7 * g__32 * g__32)) * (g__13 * g__13) + 0.3e1 * g__12
280 : : * g__23 * (g__21 * g__21 + 0.28e2 / 0.3e1 * g__22 * g__22
281 : : - 0.28e2 / 0.3e1 * g__32 * g__32) * g__13 + g__22
282 : : * ((g__21 * g__21 - (14 * g__23 * g__23)) * (g__12 * g__12)
283 : : - 0.2e1 * (g__32 * g__32) * (g__21 * g__21 + (21 * g__23 * g__23))))
284 : : * g__33 + 0.2e1 * g__32 * (- g__23 * (g__21 * g__21 - (14 * g__22
285 : : * g__22) - (14 * g__32 * g__32)) * (g__13 * g__13) / 0.2e1
286 : : - 0.3e1 / 0.2e1 * g__22 * g__12 * (g__21 * g__21 + 0.28e2 / 0.3e1
287 : : * g__23 * g__23) * g__13 + g__23 * (g__21 * g__21
288 : : + (7 * g__23 * g__23)) * (g__12 * g__12 + g__32 * g__32)))
289 : : * m_mu * std::pow(std::pow( g__33 * (g__11 * g__22 - g__12 * g__21)
290 : : + (g__12 * g__23 * g__31) + g__13 * (g__21 * g__32 - (g__22 * g__31))
291 : : - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ( ((g__12 * g__23
292 : : - g__13 * g__22) * g__31) + (g__33 * g__22 - g__23 * g__32) * g__11
293 : : - g__21 * (g__33 * g__12 - g__13 * g__32)) / 0.9e1;
294 : :
295 : : // d(sigma_11)/d(g_21)
296 : : dsigdg[0][1] = 0.2e1 / 0.9e1 * m_mu * (((-g__33 * g__12/0.2e1 - g__13 * g__32)
297 : : * g__31 * g__31 + 0.3e1 / 0.2e1 * g__11 * (-g__12 * g__13 + g__32 * g__33)
298 : : * g__31 + g__12 * (g__11 * g__11 + 0.7e1 * g__12 * g__12 + 0.7e1 * g__32
299 : : * g__32) * g__33 + g__13 * g__32 * (g__11 * g__11 - 0.14e2 * g__12 * g__12
300 : : - 0.14e2 * g__32 * g__32) / 0.2e1) * g__23 * g__23 + (0.3e1 / 0.2e1 * g__22
301 : : * (-g__12 * g__32 + g__13 * g__33) * g__31 * g__31 + ((g__33 * g__33 * g__12
302 : : + g__33 * g__13 * g__32 / 0.2e1 + 0.3e1 / 0.2e1 * g__12 * (g__12 * g__12
303 : : + g__13 * g__13 + g__32 * g__32)) * g__21 - 0.3e1 / 0.2e1 * g__11 * g__22
304 : : * (g__12 * g__12 - g__13 * g__13 - g__32 * g__32 + g__33 * g__33)) * g__31
305 : : - 0.3e1 / 0.2e1 * g__11 * (g__33 * g__33 * g__32 + g__12 * g__13
306 : : * g__33 / 0.3e1 + g__32 * (g__12 * g__12 + 0.2e1 / 0.3e1 * g__13 * g__13
307 : : + g__32 * g__32)) * g__21 - 0.14e2 * (g__33 * g__33 * g__12 * g__32
308 : : + (0.3e1 / 0.28e2 * g__11 * g__11 * g__13 + g__12 * g__12 * g__13 - g__13
309 : : * g__32 * g__32) * g__33 - 0.3e1 / 0.28e2 * (g__11 * g__11 + 0.28e2 / 0.3e1
310 : : * g__13 * g__13) * g__12 * g__32) * g__22) * g__23
311 : : + (g__12 * (g__12 * g__12 + g__13 * g__13 + g__22 * g__22) * g__33
312 : : - g__13 * (g__12 * g__12 + g__13 * g__13 - g__22 * g__22 / 0.2e1) * g__32)
313 : : * g__31 * g__31 + (-0.3e1 / 0.2e1 * g__22 * (g__33 * g__33 * g__13 + g__33
314 : : * g__12 * g__32 / 0.3e1 + g__13 * (g__12 * g__12 + g__13 * g__13
315 : : + 0.2e1 / 0.3e1 * g__32 * g__32)) * g__21 - 0.2e1 * g__11 * (g__33 * g__33
316 : : * g__12 * g__13 + g__32 * (g__12 * g__12 - g__13 * g__13 + 0.3e1 / 0.4e1
317 : : * g__22 * g__22) * g__33 - 0.3e1 / 0.4e1 * g__13 * g__12 * (g__22 * g__22
318 : : + 0.4e1 / 0.3e1 * g__32 * g__32))) * g__31 - (g__33 * g__12 - g__13 * g__32)
319 : : * (g__12 * g__12 + g__13 * g__13 + g__32 * g__32 + g__33 * g__33) * g__21
320 : : * g__21 / 0.2e1 + 0.3e1 / 0.2e1 * g__22 * g__11 * (std::pow(g__33, 0.3e1)
321 : : + (0.2e1 / 0.3e1 * g__12 * g__12 + g__13 * g__13 + g__32 * g__32) * g__33
322 : : + g__12 * g__13 * g__32 / 0.3e1) * g__21 + g__12 * (g__11 * g__11 + 0.7e1
323 : : * g__12 * g__12 + 0.7e1 * g__22 * g__22) * std::pow(g__33, 0.3e1) - g__13
324 : : * g__32 * (g__11 * g__11 + 0.21e2 * g__12 * g__12 + 0.7e1 * g__22 * g__22)
325 : : * g__33 * g__33 - ((g__11 * g__11 - 0.14e2 * g__13 * g__13) * g__22 * g__22
326 : : - 0.2e1 * g__32 * g__32 * (g__11 * g__11 + 0.21e2 * g__13 * g__13)) * g__12
327 : : * g__33 / 0.2e1 - g__13 * g__32 * (g__22 * g__22 + g__32 * g__32) * (g__11
328 : : * g__11 + 0.7e1 * g__13 * g__13)) * std::pow(std::pow(g__33 * (g__11 * g__22
329 : : - g__12 * g__21) + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22
330 : : * g__31) - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ((-g__11
331 : : * g__32 + g__12 * g__31) * g__23 - g__13 * g__22 * g__31
332 : : + (-g__33 * g__12 + g__13 * g__32) * g__21 + g__22 * g__33 * g__11);
333 : :
334 : : // d(sigma_11)/d(g_31)
335 : : dsigdg[0][2] = -0.2e1 / 0.9e1 * m_mu * ((-g__12 * g__31 * g__31 / 0.2e1
336 : : + 0.3e1 / 0.2e1 * g__11 * g__31 * g__32 + g__12 * (g__11 * g__11
337 : : + 0.7e1 * g__12 * g__12 + 0.7e1 * g__32 * g__32)) * std::pow(g__23, 0.3e1)
338 : : + (((-0.3e1 / 0.2e1 * g__32 * g__11 + g__12 * g__31) * g__33
339 : : - 0.2e1 * (g__11 * g__12 + 0.3e1 / 0.4e1 * g__31 * g__32) * g__13) * g__21
340 : : - 0.3e1 / 0.2e1 * ((g__11 * g__31 + 0.28e2 / 0.3e1 * g__12 * g__32) * g__33
341 : : + 0.2e1 / 0.3e1 * g__13 * (g__11 * g__11 + 0.21e2 * g__12 * g__12
342 : : - g__31 * g__31 / 0.2e1 + 0.7e1 * g__32 * g__32)) * g__22) * g__23 * g__23
343 : : + ((-g__33 * g__33 * g__12 / 0.2e1 + 0.3e1 / 0.2e1 * g__33 * g__13 * g__32
344 : : + g__12 * (g__12 * g__12 + g__13 * g__13 + g__32 * g__32)) * g__21 * g__21
345 : : + 0.3e1 / 0.2e1 * (g__33 * g__33 * g__11 + g__33 * g__13 * g__31 / 0.3e1
346 : : - g__12 * g__31 * g__32 / 0.3e1 - 0.4e1 / 0.3e1 * g__11 * (g__12 * g__12
347 : : - g__13 * g__13 + 0.3e1 / 0.4e1 * g__32 * g__32)) * g__22 * g__21
348 : : + g__12 * (g__11 * g__11 + 0.7e1 * g__12 * g__12 + 0.7e1 * g__22 * g__22)
349 : : * g__33 * g__33 - 0.3e1 / 0.2e1 * g__13 * (g__11 * g__12 * g__31 / 0.3e1
350 : : + g__32 * (g__11 * g__11 + 0.28e2 / 0.3e1 * g__12 * g__12 - 0.28e2 / 0.3e1
351 : : * g__22 * g__22)) * g__33 - g__12 * (g__12 * g__12 + g__13 * g__13
352 : : + g__22 * g__22) * g__31 * g__31 / 0.2e1 + g__32 * g__11 * (g__12 * g__12
353 : : + 0.3e1 / 0.2e1 * g__13 * g__13 + 0.3e1 / 0.2e1 * g__22 * g__22) * g__31
354 : : + g__12 * ((g__22 * g__22 - g__32 * g__32 / 0.2e1) * g__11 * g__11
355 : : + (0.21e2 * g__22 * g__22 + 0.7e1 * g__32 * g__32) * g__13 * g__13))
356 : : * g__23 - g__22 * (g__33 * g__33 * g__13 + 0.3e1 / 0.2e1 * g__33
357 : : * g__12 * g__32 + g__13 * (g__12 * g__12 + g__13 * g__13 - g__32
358 : : * g__32 / 0.2e1)) * g__21 * g__21 + (-0.3e1 / 0.2e1 * g__33 * g__33
359 : : * g__11 * g__12 * g__13 + (0.3e1 / 0.2e1 * g__12 * (g__12 * g__12
360 : : + g__13 * g__13 + g__22 * g__22) * g__31 - 0.3e1 / 0.2e1 * g__32
361 : : * g__11 * (g__12 * g__12 - g__13 * g__13 - g__22 * g__22)) * g__33
362 : : + 0.2e1 * g__13 * (-0.3e1 / 0.4e1 * (g__12 * g__12 + g__13 * g__13
363 : : + 0.2e1 / 0.3e1 * g__22 * g__22) * g__32 * g__31 + g__11 * g__12
364 : : * (g__22 * g__22 + 0.3e1 / 0.4e1 * g__32 * g__32))) * g__21 + g__22
365 : : * (g__13 * (g__11 * g__11 - 0.14e2 * g__12 * g__12 - 0.14e2 * g__22 * g__22)
366 : : * g__33 * g__33 + (-0.3e1 * g__11 * (g__12 * g__12 + 0.2e1 / 0.3e1 * g__13
367 : : * g__13 + g__22 * g__22) * g__31 + 0.3e1 * (g__11 * g__11 + 0.28e2 / 0.3e1
368 : : * g__13 * g__13) * g__12 * g__32) * g__33 - 0.2e1 * g__13 * ((-g__12
369 : : * g__12 / 0.2e1 - g__13 * g__13 / 0.2e1 - g__22 * g__22 / 0.2e1)
370 : : * g__31 * g__31 - g__11 * g__12 * g__31 * g__32 / 0.2e1 + (g__22 * g__22
371 : : + g__32 * g__32) * (g__11 * g__11 + 0.7e1 * g__13 * g__13))) / 0.2e1)
372 : : * std::pow(std::pow(g__33 * (g__11 * g__22 - g__12 * g__21) + g__12
373 : : * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22 * g__31) - g__11
374 : : * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ((-g__32 * g__11 + g__12 * g__31)
375 : : * g__23 + (-g__33 * g__12 + g__13 * g__32) * g__21 + (g__11 * g__33 - g__13
376 : : * g__31) * g__22);
377 : :
378 : : // d(sigma_21)/d(g_11)
379 : : dsigdg[1][0] = ( ((g__33 * (3 * g__13 * g__13 + 4 * g__23 * g__23) * g__21
380 : : + ((-4 * g__13 * g__13 - 4 * g__23 * g__23) * g__31 + g__33 * g__11
381 : : * g__13) * g__23) * g__32 * g__32) + ((((g__13 * g__23 * g__23
382 : : - 6 * g__33 * g__33 * g__13) * g__21 - g__23 * (-7 * g__33 * g__13
383 : : * g__31 + g__11 * (g__23 * g__23 + g__33 * g__33))) * g__12 - ((g__13
384 : : * g__13 * g__23 + 8 * g__23 * g__33 * g__33) * g__21 - g__33 * (g__13
385 : : * g__13 + 8 * g__23 * g__23) * g__31 + g__13 * g__11 * (g__33 - g__23)
386 : : * (g__33 + g__23)) * g__22) * g__32) + (3 * (g__23 * g__23 + g__33 * g__33)
387 : : * (g__21 * g__33 - g__23 * g__31) * g__12 * g__12) + (g__22 * (-7 * g__33
388 : : * g__13 * g__23 * g__21 - g__13 * (-6 * g__23 * g__23 + g__33 * g__33)
389 : : * g__31 + g__33 * g__11 * (g__23 * g__23 + g__33 * g__33)) * g__12)
390 : : + 0.4e1 * (g__22 * g__22) * ( ((g__13 * g__13 * g__33
391 : : + std::pow( g__33, 3)) * g__21) - (( (g__33 * g__33) + 0.3e1 / 0.4e1
392 : : * g__13 * g__13) * g__31 + (g__33 * g__11 * g__13) / 0.4e1)
393 : : * g__23)) * m_mu * std::pow( std::pow( (g__33 * (g__11 * g__22
394 : : - g__12 * g__21) + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32
395 : : - g__22 * g__31) - g__11 * g__23 * g__32), 2), -0.2e1 / 0.3e1)
396 : : / ((-g__11 * g__23 + g__13 * g__21) * g__32 + (-g__21 * g__33
397 : : + g__23 * g__31) * g__12 + (g__11 * g__33 - g__13 * g__31)
398 : : * g__22) / 0.3e1;
399 : :
400 : : // d(sigma_21)/d(g_21)
401 : : dsigdg[1][1] = -0.4e1 / 0.3e1 * ((g__33 * (g__13 * g__13
402 : : + 0.3e1 / 0.4e1 * g__23 * g__23) * g__11
403 : : + ((-0.4e1 * g__13 * g__13 - 0.4e1 * g__23 * g__23)
404 : : * g__31 + g__33 * g__21 * g__23) * g__13 / 0.4e1) * g__32
405 : : * g__32 + ((-0.3e1 / 0.2e1 * (g__33 * g__33 - g__13 * g__13 / 0.6e1)
406 : : * g__23 * g__11 - g__13 * (-0.7e1 * g__33 * g__23 * g__31 + g__21
407 : : * (g__13 * g__13 + g__33 * g__33)) / 0.4e1) * g__22 - 0.2e1 * (g__13
408 : : * (g__33 * g__33 + g__23 * g__23 / 0.8e1) * g__11 - g__33 * (g__13
409 : : * g__13 + g__23 * g__23 / 0.8e1) * g__31 + g__21 * g__23 * (g__33 - g__13)
410 : : * (g__33 + g__13) / 0.8e1) * g__12) * g__32 + 0.3e1 / 0.4e1
411 : : * (g__13 * g__13 + g__33 * g__33) * (g__11 * g__33 - g__13 * g__31)
412 : : * g__22 * g__22 + g__12 * (-0.7e1 * g__33 * g__11 * g__13 * g__23 - g__23
413 : : * (-0.6e1 * g__13 * g__13 + g__33 * g__33) * g__31 + g__33 * g__21
414 : : * (g__13 * g__13 + g__33 * g__33)) * g__22 / 0.4e1 + ((g__23 * g__23
415 : : * g__33 + std::pow(g__33, 0.3e1)) * g__11 - g__13 * ((g__33 * g__33
416 : : + 0.3e1 / 0.4e1 * g__23 * g__23) * g__31 + g__33 * g__21 * g__23 / 0.4e1))
417 : : * g__12 * g__12) * m_mu * std::pow(std::pow(g__33 * (g__11 * g__22 - g__12
418 : : * g__21) + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22 * g__31)
419 : : - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ((-g__11 * g__23 + g__13
420 : : * g__21) * g__32 + (g__11 * g__33 - g__13 * g__31) * g__22 - g__12
421 : : * (g__21 * g__33 - g__23 * g__31));
422 : :
423 : : // d(sigma_21)/d(g_31)
424 : : dsigdg[1][2] = 0.4e1 / 0.3e1 * m_mu * (((-std::pow(g__13, 0.3e1)
425 : : - g__33 * g__33 * g__13) * g__21 + 0.3e1 / 0.4e1 * ((g__33 * g__33
426 : : + 0.4e1 / 0.3e1 * g__13 * g__13) * g__11 + g__33 * g__13 * g__31 / 0.3e1)
427 : : * g__23) * g__22 * g__22 + ((0.7e1 / 0.4e1 * g__33 * g__13 * g__23 * g__21
428 : : + g__33 * (g__13 * g__13 - 0.6e1 * g__23 * g__23) * g__11 / 0.4e1
429 : : - g__13 * g__31 * (g__13 * g__13 + g__23 * g__23) / 0.4e1) * g__32
430 : : - (-g__23 * (0.8e1 * g__13 * g__13 + g__33 * g__33) * g__21
431 : : + g__13 * (0.8e1 * g__23 * g__23 + g__33 * g__33) * g__11
432 : : - g__33 * g__31 * (g__13 - g__23) * (g__13 + g__23)) * g__12 / 0.4e1)
433 : : * g__22 + 0.3e1 / 0.4e1 * (g__13 * g__13 + g__23 * g__23)
434 : : * (g__11 * g__23 - g__13 * g__21) * g__32 * g__32 - 0.7e1 / 0.4e1
435 : : * (-0.6e1 / 0.7e1 * g__33 * (g__13 * g__13 - g__23 * g__23 / 0.6e1)
436 : : * g__21 + (g__33 * g__11 * g__13 - g__31 * (g__13 * g__13 + g__23
437 : : * g__23) / 0.7e1) * g__23) * g__12 * g__32 + g__12 * g__12
438 : : * ((-0.3e1 / 0.4e1 * g__33 * g__33 * g__13 - g__13 * g__23 * g__23)
439 : : * g__21 + g__23 * (g__11 * (g__23 * g__23 + g__33 * g__33)
440 : : - g__33 * g__13 * g__31 / 0.4e1))) * std::pow(std::pow(g__33
441 : : * (g__11 * g__22 - g__12 * g__21) + g__12 * g__23 * g__31 + g__13
442 : : * (g__21 * g__32 - g__22 * g__31) - g__11 * g__23 * g__32, 0.2e1),
443 : : -0.2e1 / 0.3e1) / ((-g__11 * g__23 + g__13 * g__21) * g__32
444 : : + (g__11 * g__33 - g__13 * g__31) * g__22 - g__12 * (g__21
445 : : * g__33 - g__23 * g__31));
446 : :
447 : : // d(sigma_31)/d(g_11)
448 : : dsigdg[2][0] = - ((g__32 * (3 * g__12 * g__12 + 4 * g__22 * g__22)
449 : : * g__21 + g__22 * ((-4 * g__12 * g__12 - 4 * g__22 * g__22) * g__31
450 : : + g__11 * g__12 * g__32)) * g__33 * g__33 + ((g__12 * (g__22 * g__22
451 : : - 6 * g__32 * g__32) * g__21 - g__22 * (-7 * g__12 * g__31 * g__32
452 : : + g__11 * (g__22 * g__22 + g__32 * g__32))) * g__13
453 : : + (-g__22 * (g__12 * g__12 + 8 * g__32 * g__32) * g__21 + g__32
454 : : * (g__12 * g__12 + 8 * g__22 * g__22) * g__31 + g__11 * g__12
455 : : * (g__22 - g__32) * (g__22 + g__32)) * g__23) * g__33
456 : : + 3 * (g__22 * g__22 + g__32 * g__32) * (g__21 * g__32 - g__22 * g__31)
457 : : * g__13 * g__13 + g__23 * (-7 * g__12 * g__22 * g__32 * g__21
458 : : + (6 * g__12 * g__22 * g__22 - g__12 * g__32 * g__32) * g__31
459 : : + g__32 * g__11 * (g__22 * g__22 + g__32 * g__32)) * g__13 - g__23
460 : : * g__23 * ((-4 * g__12 * g__12 * g__32 - 4 * std::pow( g__32, 3))
461 : : * g__21 + g__22 * (g__31 * (3 * g__12 * g__12 + 4 * g__32 * g__32)
462 : : + g__11 * g__12 * g__32))) * m_mu * std::pow( std::pow( (g__33
463 : : * (g__11 * g__22 - g__12 * g__21) + g__12 * g__23 * g__31 + g__13
464 : : * (g__21 * g__32 - g__22 * g__31) - g__11 * g__23 * g__32), 2),
465 : : -0.2e1 / 0.3e1) / (g__33 * (g__11 * g__22 - g__12 * g__21) + g__13
466 : : * (g__21 * g__32 - g__22 * g__31) - g__23 * (g__32 * g__11 - g__12
467 : : * g__31)) / 0.3e1;
468 : :
469 : : // d(sigma_31)/d(g_21)
470 : : dsigdg[2][1] = 0.4e1 / 0.3e1 * (((-std::pow(g__12, 0.3e1)
471 : : - g__12 * g__22 * g__22) * g__31 + ((g__12 * g__12 + 0.3e1 / 0.4e1
472 : : * g__22 * g__22) * g__11 + g__12 * g__21 * g__22 / 0.4e1) * g__32)
473 : : * g__33 * g__33 + ((0.7e1 / 0.4e1 * g__12 * g__22 * g__31 * g__32
474 : : + g__22 * (g__12 * g__12 - 0.6e1 * g__32 * g__32) * g__11 / 0.4e1
475 : : - g__12 * g__21 * (g__12 * g__12 + g__32 * g__32) / 0.4e1) * g__23
476 : : - g__13 * ((-0.8e1 * g__12 * g__12 - g__22 * g__22) * g__32 * g__31
477 : : + g__12 * (g__22 * g__22 + 0.8e1 * g__32 * g__32) * g__11 - g__21
478 : : * g__22 * (g__12 - g__32) * (g__12 + g__32)) / 0.4e1) * g__33
479 : : + 0.3e1 / 0.4e1 * (g__12 * g__12 + g__32 * g__32) * (g__32 * g__11
480 : : - g__12 * g__31) * g__23 * g__23 - 0.7e1 / 0.4e1 * g__13
481 : : * (-0.6e1 / 0.7e1 * g__22 * (g__12 * g__12 - g__32 * g__32 / 0.6e1)
482 : : * g__31 + g__32 * (g__11 * g__12 * g__22 - g__21 * (g__12 * g__12
483 : : + g__32 * g__32) / 0.7e1)) * g__23 + g__13 * g__13 * ((-0.3e1 / 0.4e1
484 : : * g__12 * g__22 * g__22 - g__12 * g__32 * g__32) * g__31
485 : : + (g__11 * (g__22 * g__22 + g__32 * g__32) - g__12 * g__21
486 : : * g__22 / 0.4e1) * g__32)) * m_mu * std::pow(std::pow(g__33
487 : : * (g__11 * g__22 - g__12 * g__21) + g__12 * g__23 * g__31 + g__13
488 : : * (g__21 * g__32 - g__22 * g__31) - g__11 * g__23 * g__32, 0.2e1),
489 : : -0.2e1 / 0.3e1) / (g__13 * (g__21 * g__32 - g__22 * g__31)
490 : : + (-g__32 * g__11 + g__12 * g__31) * g__23 + g__33
491 : : * (g__11 * g__22 - g__12 * g__21));
492 : :
493 : : // d(sigma_31)/d(g_31)
494 : : dsigdg[2][2] = -(((-0.4e1 / 0.3e1 * g__12 * g__32 * g__32
495 : : - 0.4e1 / 0.3e1 * std::pow(g__12, 0.3e1)) * g__21 + 0.4e1 / 0.3e1
496 : : * g__22 * ((g__12 * g__12 + 0.3e1 / 0.4e1 * g__32 * g__32) * g__11
497 : : + g__12 * g__31 * g__32 / 0.4e1)) * g__23 * g__23 + ((0.7e1 / 0.3e1
498 : : * g__12 * g__22 * g__32 * g__21 + g__32 * (g__12 * g__12 - 0.6e1
499 : : * g__22 * g__22) * g__11 / 0.3e1 - g__12 * g__31 * (g__12 * g__12
500 : : + g__22 * g__22) / 0.3e1) * g__33 - 0.8e1 / 0.3e1 * g__13 * (-g__22
501 : : * (g__12 * g__12 + g__32 * g__32 / 0.8e1) * g__21 + g__12 * (g__22
502 : : * g__22 + g__32 * g__32 / 0.8e1) * g__11 - g__31 * g__32 * (g__12
503 : : - g__22) * (g__12 + g__22) / 0.8e1)) * g__23 + (g__12 * g__12
504 : : + g__22 * g__22) * (g__11 * g__22 - g__12 * g__21) * g__33 * g__33
505 : : - 0.7e1 / 0.3e1 * (-0.6e1 / 0.7e1 * g__32 * (g__12 * g__12 - g__22
506 : : * g__22 / 0.6e1) * g__21 + (g__11 * g__12 * g__32 - g__31 * (g__12
507 : : * g__12 + g__22 * g__22) / 0.7e1) * g__22) * g__13 * g__33
508 : : + 0.4e1 / 0.3e1 * (-g__12 * (g__22 * g__22 + 0.3e1 / 0.4e1
509 : : * g__32 * g__32) * g__21 + g__22 * (g__11 * (g__22 * g__22
510 : : + g__32 * g__32) - g__12 * g__31 * g__32 / 0.4e1)) * g__13 * g__13)
511 : : * m_mu * std::pow(std::pow(g__33 * (g__11 * g__22 - g__12 * g__21)
512 : : + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22 * g__31)
513 : : - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / (g__13 * (g__21
514 : : * g__32 - g__22 * g__31) + (-g__32 * g__11 + g__12 * g__31) * g__23
515 : : + g__33 * (g__11 * g__22 - g__12 * g__21));
516 : :
517 : : // Define amat
518 : : double amat[9];
519 : : for (std::size_t i=0; i<3; ++i)
520 : : for (std::size_t j=0; j<3; ++j)
521 : : {
522 : : amat[i*3+j] = 0.0;
523 : : for (std::size_t k=0; k<3; ++k)
524 : : {
525 : : amat[i*3+j] -= defgrad[k][j]*dsigdg[i][k];
526 : : }
527 : : amat[i*3+j] /= (arho/alpha);
528 : : }
529 : :
530 : : double eig_real[3], eig_imag[3];
531 : : double vl[3], vr[3];
532 : : #ifndef NDEBUG
533 : : lapack_int ierr =
534 : : #endif
535 : : LAPACKE_dgeev (LAPACK_ROW_MAJOR, 'N', 'N', 3, amat, 3,
536 : : eig_real, eig_imag, vl, 3, vr, 3);
537 : : Assert(ierr==0, "Lapack failed to compute eigenvalues");
538 : :
539 : : // Manually find max
540 : : tk::real eig_max = eig_real[0];
541 : : for (std::size_t i=1; i<3; i++)
542 : : if (eig_real[i] > eig_max)
543 : : eig_max = eig_real[i];
544 : :
545 : : */
546 : :
547 : : // Approximated elastic contribution, from Barton, P. T. (2019).
548 : : // An interface-capturing Godunov method for the simulation of compressible
549 : : // solid-fluid problems. Journal of Computational Physics, 390, 25-50
550 : 0 : tk::real a = (4.0/3.0) * m_mu * alpha / arho;
551 : :
552 : : // hydrodynamic contribution
553 : 0 : auto p_eff = std::max( 1.0e-15, apr+(alpha*m_pstiff) );
554 : 0 : a += m_gamma * p_eff / arho;
555 : :
556 : : // Compute square root
557 : 0 : a = std::sqrt(a);
558 : :
559 : : // check sound speed divergence
560 [ - - ]: 0 : if (!std::isfinite(a)) {
561 : 0 : std::cout << "Material-id: " << imat << std::endl;
562 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
563 : 0 : std::cout << "Partial density: " << arho << std::endl;
564 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
565 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
566 : : + std::to_string(a) + ", material volume fraction: " +
567 : : std::to_string(alpha));
568 : : }
569 : :
570 : 0 : return a;
571 : : }
572 : :
573 : : tk::real
574 : 0 : SmallShearSolid::shearspeed(
575 : : tk::real arho,
576 : : tk::real alpha,
577 : : std::size_t imat ) const
578 : : // *************************************************************************
579 : : //! Calculate speed of sound from the material density and material pressure
580 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
581 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
582 : : //! the single-material system, this argument can be left unspecified by
583 : : //! the calling code
584 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
585 : : //! for the single-material system, this argument can be left unspecified
586 : : //! by the calling code
587 : : //! \return Material shear-wave speed speed using the SmallShearSolid EoS
588 : : // *************************************************************************
589 : : {
590 : : // Approximate shear-wave speed. Ref. Barton, P. T. (2019).
591 : : // An interface-capturing Godunov method for the simulation of compressible
592 : : // solid-fluid problems. Journal of Computational Physics, 390, 25-50.
593 : 0 : tk::real a = std::sqrt(alpha*m_mu/arho);
594 : :
595 : : // check shear-wave speed divergence
596 [ - - ]: 0 : if (!std::isfinite(a)) {
597 : 0 : std::cout << "Material-id: " << imat << std::endl;
598 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
599 : 0 : std::cout << "Partial density: " << arho << std::endl;
600 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf shear-wave speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
601 : : + std::to_string(a) + ", material volume fraction: " +
602 : : std::to_string(alpha));
603 : : }
604 : :
605 : 0 : return a;
606 : : }
607 : :
608 : : tk::real
609 : 0 : SmallShearSolid::totalenergy(
610 : : tk::real rho,
611 : : tk::real u,
612 : : tk::real v,
613 : : tk::real w,
614 : : tk::real pr,
615 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
616 : : // *************************************************************************
617 : : //! \brief Calculate material specific total energy from the material
618 : : //! density, momentum and material pressure
619 : : //! \param[in] rho Material density
620 : : //! \param[in] u X-velocity
621 : : //! \param[in] v Y-velocity
622 : : //! \param[in] w Z-velocity
623 : : //! \param[in] pr Material pressure
624 : : //! \param[in] defgrad Material inverse deformation gradient tensor
625 : : //! g_k. Default is 0, so that for the single-material system,
626 : : //! this argument can be left unspecified by the calling code
627 : : //! \return Material specific total energy using the SmallShearSolid EoS
628 : : // *************************************************************************
629 : : {
630 : : // obtain hydro contribution to energy
631 : 0 : tk::real rhoEh = (pr + m_pstiff) / (m_gamma-1.0) + 0.5 * rho *
632 : 0 : (u*u + v*v + w*w) + m_pstiff;
633 : : // obtain elastic contribution to energy
634 : : tk::real eps2;
635 [ - - ]: 0 : tk::real rhoEe = elasticEnergy(defgrad, eps2);
636 : :
637 : 0 : return (rhoEh + rhoEe);
638 : : }
639 : :
640 : : tk::real
641 : 0 : SmallShearSolid::temperature(
642 : : tk::real arho,
643 : : tk::real u,
644 : : tk::real v,
645 : : tk::real w,
646 : : tk::real arhoE,
647 : : tk::real alpha,
648 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
649 : : // *************************************************************************
650 : : //! \brief Calculate material temperature from the material density, and
651 : : //! material specific total energy
652 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
653 : : //! \param[in] u X-velocity
654 : : //! \param[in] v Y-velocity
655 : : //! \param[in] w Z-velocity
656 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
657 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
658 : : //! the single-material system, this argument can be left unspecified by
659 : : //! the calling code
660 : : //! \param[in] defgrad Material inverse deformation gradient tensor
661 : : //! (g_k). Default is 0, so that for the single-material system,
662 : : //! this argument can be left unspecified by the calling code
663 : : //! \return Material temperature using the SmallShearSolid EoS
664 : : // *************************************************************************
665 : : {
666 : : // obtain elastic contribution to energy
667 : : tk::real eps2;
668 [ - - ]: 0 : auto arhoEe = alpha*elasticEnergy(defgrad, eps2);
669 : : // obtain hydro contribution to energy
670 : 0 : auto arhoEh = arhoE - arhoEe;
671 : :
672 : 0 : tk::real t = (arhoEh - 0.5 * arho * (u*u + v*v + w*w) - alpha*m_pstiff)
673 : 0 : / (arho*m_cv);
674 : 0 : return t;
675 : : }
676 : :
677 : : tk::real
678 : 0 : SmallShearSolid::min_eff_pressure(
679 : : tk::real min,
680 : : tk::real,
681 : : tk::real ) const
682 : : // *************************************************************************
683 : : //! Compute the minimum allowed pressure
684 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
685 : : //! \return Minimum pressure allowed by physical constraints
686 : : // *************************************************************************
687 : : {
688 : : // minimum pressure is constrained by zero soundspeed.
689 : 0 : return (min - m_pstiff);
690 : : }
691 : :
692 : : tk::real
693 : 0 : SmallShearSolid::elasticEnergy(
694 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad,
695 : : tk::real& eps2 ) const
696 : : // *************************************************************************
697 : : //! \brief Calculate elastic contribution to material energy from the material
698 : : //! density, and deformation gradient tensor
699 : : //! \param[in] defgrad Material inverse deformation gradient tensor
700 : : //! \param[in/out] eps2 Elastic shear distortion
701 : : //! \return Material elastic energy using the SmallShearSolid EoS
702 : : //! \details This function returns the material elastic energy, and also stores
703 : : //! the elastic shear distortion for further use
704 : : // *************************************************************************
705 : : {
706 : : // compute Right Cauchy-Green strain tensor
707 [ - - ]: 0 : auto Ct = tk::getRightCauchyGreen(defgrad);
708 : 0 : auto detC = std::pow(tk::determinant(Ct), 1.0/3.0);
709 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
710 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
711 : 0 : Ct[i][j] /= detC;
712 : : }
713 : :
714 : : // compute elastic shear distortion
715 : 0 : eps2 = 0.5 * (Ct[0][0]+Ct[1][1]+Ct[2][2] - 3.0);
716 : :
717 : : // compute elastic energy
718 : 0 : auto rhoEe = m_mu * eps2;
719 : :
720 : 0 : return rhoEe;
721 : : }
|
