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*/ ) const
213 : : // *************************************************************************
214 : : //! Calculate speed of sound from the material density and material pressure
215 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
216 : : //! \param[in] apr Material partial pressure (alpha_k * p_k)
217 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
218 : : //! the single-material system, this argument can be left unspecified by
219 : : //! the calling code
220 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
221 : : //! for the single-material system, this argument can be left unspecified
222 : : //! by the calling code
223 : : //! (alpha * sigma_ij * n_j) projected onto the normal vector. Default is 0,
224 : : //! so that for the single-material system, this argument can be left
225 : : //! unspecified by the calling code
226 : : // //! \param[in] defgrad Material inverse deformation gradient tensor
227 : : // //! (g_k) with the first dimension aligned to direction in which
228 : : // //! wave speeds are required. Default is 0, so that for the single-material
229 : : // //! system, this argument can be left unspecified by the calling code
230 : : //! \return Material speed of sound using the SmallShearSolid EoS
231 : : // *************************************************************************
232 : : {
233 : : /* Rigorous approach: Eigenvalues of full elastic tensor
234 : :
235 : : // deformation gradient
236 : : tk::real g__11 = defgrad[0][0];
237 : : tk::real g__12 = defgrad[0][1];
238 : : tk::real g__13 = defgrad[0][2];
239 : : tk::real g__21 = defgrad[1][0];
240 : : tk::real g__22 = defgrad[1][1];
241 : : tk::real g__23 = defgrad[1][2];
242 : : tk::real g__31 = defgrad[2][0];
243 : : tk::real g__32 = defgrad[2][1];
244 : : tk::real g__33 = defgrad[2][2];
245 : :
246 : : std::array< std::array< tk::real, 3 >, 3> dsigdg;
247 : :
248 : : // d(sigma_11)/d(g_11)
249 : : dsigdg[0][0] = ( (((g__13 * g__13 * g__22 - 3 * g__12 * g__13 * g__23
250 : : - 2 * g__22 * (g__12 * g__12 + g__22 * g__22 + g__23 * g__23)) * g__33
251 : : - (-2 * g__13 * g__13 * g__23 - 3 * g__12 * g__13 * g__22
252 : : + g__23 * (g__12 * g__12 - 2 * g__22 * g__22 - 2 * g__23 * g__23)) * g__32)
253 : : * g__31 * g__31) + (( ((3 * g__12 * g__23 - 2 * g__13 * g__22) * g__33
254 : : * g__33) + (g__32 * (g__12 * g__22 - g__13 * g__23) * g__33) - (3 * g__22
255 : : * (g__22 * g__22 + g__23 * g__23 + g__32 * g__32) * g__13) + 0.3e1 * g__12
256 : : * g__23 * ( (g__22 * g__22) + (g__23 * g__23) + 0.2e1 / 0.3e1 * g__32
257 : : * g__32)) * g__11 + 0.3e1 * (( (g__12 * g__13) + 0.4e1 / 0.3e1 * g__22
258 : : * g__23) * (g__33 * g__33) + g__32 * ( (g__12 * g__12) - (g__13 * g__13)
259 : : + 0.4e1 / 0.3e1 * g__22 * g__22 - 0.4e1 / 0.3e1 * g__23 * g__23) * g__33
260 : : + (g__13 * g__13 * g__22 * g__23) + (g__12 * (g__22 * g__22 - g__23 * g__23
261 : : - g__32 * g__32) * g__13) - g__22 * g__23 * ( (g__12 * g__12)
262 : : + 0.4e1 / 0.3e1 * g__32 * g__32)) * g__21) * g__31 + (g__33 * g__22
263 : : - g__23 * g__32) * (g__22 * g__22 + g__23 * g__23 + g__32 * g__32 + g__33
264 : : * g__33) * g__11 * g__11 - 0.3e1 * ( ( std::pow( g__33, 3) * g__12)
265 : : - (g__33 * g__33 * g__13 * g__32) + (- (g__13 * g__22 * g__23) / 0.3e1
266 : : + 0.2e1 / 0.3e1 * g__12 * ( (g__22 * g__22) + 0.3e1 / 0.2e1 * g__23
267 : : * g__23 + 0.3e1 / 0.2e1 * g__32 * g__32)) * g__33 + (((-3 * g__22
268 : : * g__22 - 2 * g__23 * g__23 - 3 * g__32 * g__32) * g__13 + g__12 * g__22
269 : : * g__23) * g__32) / 0.3e1) * g__21 * g__11 + (- (14 * g__12 * g__12 * g__22)
270 : : - 0.2e1 * g__21 * g__21 * g__22 - (14 * std::pow( g__22, 3)))
271 : : * std::pow( g__33, 3) + 0.14e2 * ( (2 * g__12 * g__13 * g__22)
272 : : + g__23 * ( (g__12 * g__12) + g__21 * g__21 / 0.7e1
273 : : + (3 * g__22 * g__22))) * g__32 * (g__33 * g__33)
274 : : + (-0.2e1 * g__22 * (g__21 * g__21 + (7 * g__22 * g__22)
275 : : + (7 * g__32 * g__32)) * (g__13 * g__13) + 0.3e1 * g__12
276 : : * g__23 * (g__21 * g__21 + 0.28e2 / 0.3e1 * g__22 * g__22
277 : : - 0.28e2 / 0.3e1 * g__32 * g__32) * g__13 + g__22
278 : : * ((g__21 * g__21 - (14 * g__23 * g__23)) * (g__12 * g__12)
279 : : - 0.2e1 * (g__32 * g__32) * (g__21 * g__21 + (21 * g__23 * g__23))))
280 : : * g__33 + 0.2e1 * g__32 * (- g__23 * (g__21 * g__21 - (14 * g__22
281 : : * g__22) - (14 * g__32 * g__32)) * (g__13 * g__13) / 0.2e1
282 : : - 0.3e1 / 0.2e1 * g__22 * g__12 * (g__21 * g__21 + 0.28e2 / 0.3e1
283 : : * g__23 * g__23) * g__13 + g__23 * (g__21 * g__21
284 : : + (7 * g__23 * g__23)) * (g__12 * g__12 + g__32 * g__32)))
285 : : * m_mu * std::pow(std::pow( g__33 * (g__11 * g__22 - g__12 * g__21)
286 : : + (g__12 * g__23 * g__31) + g__13 * (g__21 * g__32 - (g__22 * g__31))
287 : : - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ( ((g__12 * g__23
288 : : - g__13 * g__22) * g__31) + (g__33 * g__22 - g__23 * g__32) * g__11
289 : : - g__21 * (g__33 * g__12 - g__13 * g__32)) / 0.9e1;
290 : :
291 : : // d(sigma_11)/d(g_21)
292 : : dsigdg[0][1] = 0.2e1 / 0.9e1 * m_mu * (((-g__33 * g__12/0.2e1 - g__13 * g__32)
293 : : * g__31 * g__31 + 0.3e1 / 0.2e1 * g__11 * (-g__12 * g__13 + g__32 * g__33)
294 : : * g__31 + g__12 * (g__11 * g__11 + 0.7e1 * g__12 * g__12 + 0.7e1 * g__32
295 : : * g__32) * g__33 + g__13 * g__32 * (g__11 * g__11 - 0.14e2 * g__12 * g__12
296 : : - 0.14e2 * g__32 * g__32) / 0.2e1) * g__23 * g__23 + (0.3e1 / 0.2e1 * g__22
297 : : * (-g__12 * g__32 + g__13 * g__33) * g__31 * g__31 + ((g__33 * g__33 * g__12
298 : : + g__33 * g__13 * g__32 / 0.2e1 + 0.3e1 / 0.2e1 * g__12 * (g__12 * g__12
299 : : + g__13 * g__13 + g__32 * g__32)) * g__21 - 0.3e1 / 0.2e1 * g__11 * g__22
300 : : * (g__12 * g__12 - g__13 * g__13 - g__32 * g__32 + g__33 * g__33)) * g__31
301 : : - 0.3e1 / 0.2e1 * g__11 * (g__33 * g__33 * g__32 + g__12 * g__13
302 : : * g__33 / 0.3e1 + g__32 * (g__12 * g__12 + 0.2e1 / 0.3e1 * g__13 * g__13
303 : : + g__32 * g__32)) * g__21 - 0.14e2 * (g__33 * g__33 * g__12 * g__32
304 : : + (0.3e1 / 0.28e2 * g__11 * g__11 * g__13 + g__12 * g__12 * g__13 - g__13
305 : : * g__32 * g__32) * g__33 - 0.3e1 / 0.28e2 * (g__11 * g__11 + 0.28e2 / 0.3e1
306 : : * g__13 * g__13) * g__12 * g__32) * g__22) * g__23
307 : : + (g__12 * (g__12 * g__12 + g__13 * g__13 + g__22 * g__22) * g__33
308 : : - g__13 * (g__12 * g__12 + g__13 * g__13 - g__22 * g__22 / 0.2e1) * g__32)
309 : : * g__31 * g__31 + (-0.3e1 / 0.2e1 * g__22 * (g__33 * g__33 * g__13 + g__33
310 : : * g__12 * g__32 / 0.3e1 + g__13 * (g__12 * g__12 + g__13 * g__13
311 : : + 0.2e1 / 0.3e1 * g__32 * g__32)) * g__21 - 0.2e1 * g__11 * (g__33 * g__33
312 : : * g__12 * g__13 + g__32 * (g__12 * g__12 - g__13 * g__13 + 0.3e1 / 0.4e1
313 : : * g__22 * g__22) * g__33 - 0.3e1 / 0.4e1 * g__13 * g__12 * (g__22 * g__22
314 : : + 0.4e1 / 0.3e1 * g__32 * g__32))) * g__31 - (g__33 * g__12 - g__13 * g__32)
315 : : * (g__12 * g__12 + g__13 * g__13 + g__32 * g__32 + g__33 * g__33) * g__21
316 : : * g__21 / 0.2e1 + 0.3e1 / 0.2e1 * g__22 * g__11 * (std::pow(g__33, 0.3e1)
317 : : + (0.2e1 / 0.3e1 * g__12 * g__12 + g__13 * g__13 + g__32 * g__32) * g__33
318 : : + g__12 * g__13 * g__32 / 0.3e1) * g__21 + g__12 * (g__11 * g__11 + 0.7e1
319 : : * g__12 * g__12 + 0.7e1 * g__22 * g__22) * std::pow(g__33, 0.3e1) - g__13
320 : : * g__32 * (g__11 * g__11 + 0.21e2 * g__12 * g__12 + 0.7e1 * g__22 * g__22)
321 : : * g__33 * g__33 - ((g__11 * g__11 - 0.14e2 * g__13 * g__13) * g__22 * g__22
322 : : - 0.2e1 * g__32 * g__32 * (g__11 * g__11 + 0.21e2 * g__13 * g__13)) * g__12
323 : : * g__33 / 0.2e1 - g__13 * g__32 * (g__22 * g__22 + g__32 * g__32) * (g__11
324 : : * g__11 + 0.7e1 * g__13 * g__13)) * std::pow(std::pow(g__33 * (g__11 * g__22
325 : : - g__12 * g__21) + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22
326 : : * g__31) - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ((-g__11
327 : : * g__32 + g__12 * g__31) * g__23 - g__13 * g__22 * g__31
328 : : + (-g__33 * g__12 + g__13 * g__32) * g__21 + g__22 * g__33 * g__11);
329 : :
330 : : // d(sigma_11)/d(g_31)
331 : : dsigdg[0][2] = -0.2e1 / 0.9e1 * m_mu * ((-g__12 * g__31 * g__31 / 0.2e1
332 : : + 0.3e1 / 0.2e1 * g__11 * g__31 * g__32 + g__12 * (g__11 * g__11
333 : : + 0.7e1 * g__12 * g__12 + 0.7e1 * g__32 * g__32)) * std::pow(g__23, 0.3e1)
334 : : + (((-0.3e1 / 0.2e1 * g__32 * g__11 + g__12 * g__31) * g__33
335 : : - 0.2e1 * (g__11 * g__12 + 0.3e1 / 0.4e1 * g__31 * g__32) * g__13) * g__21
336 : : - 0.3e1 / 0.2e1 * ((g__11 * g__31 + 0.28e2 / 0.3e1 * g__12 * g__32) * g__33
337 : : + 0.2e1 / 0.3e1 * g__13 * (g__11 * g__11 + 0.21e2 * g__12 * g__12
338 : : - g__31 * g__31 / 0.2e1 + 0.7e1 * g__32 * g__32)) * g__22) * g__23 * g__23
339 : : + ((-g__33 * g__33 * g__12 / 0.2e1 + 0.3e1 / 0.2e1 * g__33 * g__13 * g__32
340 : : + g__12 * (g__12 * g__12 + g__13 * g__13 + g__32 * g__32)) * g__21 * g__21
341 : : + 0.3e1 / 0.2e1 * (g__33 * g__33 * g__11 + g__33 * g__13 * g__31 / 0.3e1
342 : : - g__12 * g__31 * g__32 / 0.3e1 - 0.4e1 / 0.3e1 * g__11 * (g__12 * g__12
343 : : - g__13 * g__13 + 0.3e1 / 0.4e1 * g__32 * g__32)) * g__22 * g__21
344 : : + g__12 * (g__11 * g__11 + 0.7e1 * g__12 * g__12 + 0.7e1 * g__22 * g__22)
345 : : * g__33 * g__33 - 0.3e1 / 0.2e1 * g__13 * (g__11 * g__12 * g__31 / 0.3e1
346 : : + g__32 * (g__11 * g__11 + 0.28e2 / 0.3e1 * g__12 * g__12 - 0.28e2 / 0.3e1
347 : : * g__22 * g__22)) * g__33 - g__12 * (g__12 * g__12 + g__13 * g__13
348 : : + g__22 * g__22) * g__31 * g__31 / 0.2e1 + g__32 * g__11 * (g__12 * g__12
349 : : + 0.3e1 / 0.2e1 * g__13 * g__13 + 0.3e1 / 0.2e1 * g__22 * g__22) * g__31
350 : : + g__12 * ((g__22 * g__22 - g__32 * g__32 / 0.2e1) * g__11 * g__11
351 : : + (0.21e2 * g__22 * g__22 + 0.7e1 * g__32 * g__32) * g__13 * g__13))
352 : : * g__23 - g__22 * (g__33 * g__33 * g__13 + 0.3e1 / 0.2e1 * g__33
353 : : * g__12 * g__32 + g__13 * (g__12 * g__12 + g__13 * g__13 - g__32
354 : : * g__32 / 0.2e1)) * g__21 * g__21 + (-0.3e1 / 0.2e1 * g__33 * g__33
355 : : * g__11 * g__12 * g__13 + (0.3e1 / 0.2e1 * g__12 * (g__12 * g__12
356 : : + g__13 * g__13 + g__22 * g__22) * g__31 - 0.3e1 / 0.2e1 * g__32
357 : : * g__11 * (g__12 * g__12 - g__13 * g__13 - g__22 * g__22)) * g__33
358 : : + 0.2e1 * g__13 * (-0.3e1 / 0.4e1 * (g__12 * g__12 + g__13 * g__13
359 : : + 0.2e1 / 0.3e1 * g__22 * g__22) * g__32 * g__31 + g__11 * g__12
360 : : * (g__22 * g__22 + 0.3e1 / 0.4e1 * g__32 * g__32))) * g__21 + g__22
361 : : * (g__13 * (g__11 * g__11 - 0.14e2 * g__12 * g__12 - 0.14e2 * g__22 * g__22)
362 : : * g__33 * g__33 + (-0.3e1 * g__11 * (g__12 * g__12 + 0.2e1 / 0.3e1 * g__13
363 : : * g__13 + g__22 * g__22) * g__31 + 0.3e1 * (g__11 * g__11 + 0.28e2 / 0.3e1
364 : : * g__13 * g__13) * g__12 * g__32) * g__33 - 0.2e1 * g__13 * ((-g__12
365 : : * g__12 / 0.2e1 - g__13 * g__13 / 0.2e1 - g__22 * g__22 / 0.2e1)
366 : : * g__31 * g__31 - g__11 * g__12 * g__31 * g__32 / 0.2e1 + (g__22 * g__22
367 : : + g__32 * g__32) * (g__11 * g__11 + 0.7e1 * g__13 * g__13))) / 0.2e1)
368 : : * std::pow(std::pow(g__33 * (g__11 * g__22 - g__12 * g__21) + g__12
369 : : * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22 * g__31) - g__11
370 : : * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ((-g__32 * g__11 + g__12 * g__31)
371 : : * g__23 + (-g__33 * g__12 + g__13 * g__32) * g__21 + (g__11 * g__33 - g__13
372 : : * g__31) * g__22);
373 : :
374 : : // d(sigma_21)/d(g_11)
375 : : dsigdg[1][0] = ( ((g__33 * (3 * g__13 * g__13 + 4 * g__23 * g__23) * g__21
376 : : + ((-4 * g__13 * g__13 - 4 * g__23 * g__23) * g__31 + g__33 * g__11
377 : : * g__13) * g__23) * g__32 * g__32) + ((((g__13 * g__23 * g__23
378 : : - 6 * g__33 * g__33 * g__13) * g__21 - g__23 * (-7 * g__33 * g__13
379 : : * g__31 + g__11 * (g__23 * g__23 + g__33 * g__33))) * g__12 - ((g__13
380 : : * g__13 * g__23 + 8 * g__23 * g__33 * g__33) * g__21 - g__33 * (g__13
381 : : * g__13 + 8 * g__23 * g__23) * g__31 + g__13 * g__11 * (g__33 - g__23)
382 : : * (g__33 + g__23)) * g__22) * g__32) + (3 * (g__23 * g__23 + g__33 * g__33)
383 : : * (g__21 * g__33 - g__23 * g__31) * g__12 * g__12) + (g__22 * (-7 * g__33
384 : : * g__13 * g__23 * g__21 - g__13 * (-6 * g__23 * g__23 + g__33 * g__33)
385 : : * g__31 + g__33 * g__11 * (g__23 * g__23 + g__33 * g__33)) * g__12)
386 : : + 0.4e1 * (g__22 * g__22) * ( ((g__13 * g__13 * g__33
387 : : + std::pow( g__33, 3)) * g__21) - (( (g__33 * g__33) + 0.3e1 / 0.4e1
388 : : * g__13 * g__13) * g__31 + (g__33 * g__11 * g__13) / 0.4e1)
389 : : * g__23)) * m_mu * std::pow( std::pow( (g__33 * (g__11 * g__22
390 : : - g__12 * g__21) + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32
391 : : - g__22 * g__31) - g__11 * g__23 * g__32), 2), -0.2e1 / 0.3e1)
392 : : / ((-g__11 * g__23 + g__13 * g__21) * g__32 + (-g__21 * g__33
393 : : + g__23 * g__31) * g__12 + (g__11 * g__33 - g__13 * g__31)
394 : : * g__22) / 0.3e1;
395 : :
396 : : // d(sigma_21)/d(g_21)
397 : : dsigdg[1][1] = -0.4e1 / 0.3e1 * ((g__33 * (g__13 * g__13
398 : : + 0.3e1 / 0.4e1 * g__23 * g__23) * g__11
399 : : + ((-0.4e1 * g__13 * g__13 - 0.4e1 * g__23 * g__23)
400 : : * g__31 + g__33 * g__21 * g__23) * g__13 / 0.4e1) * g__32
401 : : * g__32 + ((-0.3e1 / 0.2e1 * (g__33 * g__33 - g__13 * g__13 / 0.6e1)
402 : : * g__23 * g__11 - g__13 * (-0.7e1 * g__33 * g__23 * g__31 + g__21
403 : : * (g__13 * g__13 + g__33 * g__33)) / 0.4e1) * g__22 - 0.2e1 * (g__13
404 : : * (g__33 * g__33 + g__23 * g__23 / 0.8e1) * g__11 - g__33 * (g__13
405 : : * g__13 + g__23 * g__23 / 0.8e1) * g__31 + g__21 * g__23 * (g__33 - g__13)
406 : : * (g__33 + g__13) / 0.8e1) * g__12) * g__32 + 0.3e1 / 0.4e1
407 : : * (g__13 * g__13 + g__33 * g__33) * (g__11 * g__33 - g__13 * g__31)
408 : : * g__22 * g__22 + g__12 * (-0.7e1 * g__33 * g__11 * g__13 * g__23 - g__23
409 : : * (-0.6e1 * g__13 * g__13 + g__33 * g__33) * g__31 + g__33 * g__21
410 : : * (g__13 * g__13 + g__33 * g__33)) * g__22 / 0.4e1 + ((g__23 * g__23
411 : : * g__33 + std::pow(g__33, 0.3e1)) * g__11 - g__13 * ((g__33 * g__33
412 : : + 0.3e1 / 0.4e1 * g__23 * g__23) * g__31 + g__33 * g__21 * g__23 / 0.4e1))
413 : : * g__12 * g__12) * m_mu * std::pow(std::pow(g__33 * (g__11 * g__22 - g__12
414 : : * g__21) + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22 * g__31)
415 : : - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / ((-g__11 * g__23 + g__13
416 : : * g__21) * g__32 + (g__11 * g__33 - g__13 * g__31) * g__22 - g__12
417 : : * (g__21 * g__33 - g__23 * g__31));
418 : :
419 : : // d(sigma_21)/d(g_31)
420 : : dsigdg[1][2] = 0.4e1 / 0.3e1 * m_mu * (((-std::pow(g__13, 0.3e1)
421 : : - g__33 * g__33 * g__13) * g__21 + 0.3e1 / 0.4e1 * ((g__33 * g__33
422 : : + 0.4e1 / 0.3e1 * g__13 * g__13) * g__11 + g__33 * g__13 * g__31 / 0.3e1)
423 : : * g__23) * g__22 * g__22 + ((0.7e1 / 0.4e1 * g__33 * g__13 * g__23 * g__21
424 : : + g__33 * (g__13 * g__13 - 0.6e1 * g__23 * g__23) * g__11 / 0.4e1
425 : : - g__13 * g__31 * (g__13 * g__13 + g__23 * g__23) / 0.4e1) * g__32
426 : : - (-g__23 * (0.8e1 * g__13 * g__13 + g__33 * g__33) * g__21
427 : : + g__13 * (0.8e1 * g__23 * g__23 + g__33 * g__33) * g__11
428 : : - g__33 * g__31 * (g__13 - g__23) * (g__13 + g__23)) * g__12 / 0.4e1)
429 : : * g__22 + 0.3e1 / 0.4e1 * (g__13 * g__13 + g__23 * g__23)
430 : : * (g__11 * g__23 - g__13 * g__21) * g__32 * g__32 - 0.7e1 / 0.4e1
431 : : * (-0.6e1 / 0.7e1 * g__33 * (g__13 * g__13 - g__23 * g__23 / 0.6e1)
432 : : * g__21 + (g__33 * g__11 * g__13 - g__31 * (g__13 * g__13 + g__23
433 : : * g__23) / 0.7e1) * g__23) * g__12 * g__32 + g__12 * g__12
434 : : * ((-0.3e1 / 0.4e1 * g__33 * g__33 * g__13 - g__13 * g__23 * g__23)
435 : : * g__21 + g__23 * (g__11 * (g__23 * g__23 + g__33 * g__33)
436 : : - g__33 * g__13 * g__31 / 0.4e1))) * std::pow(std::pow(g__33
437 : : * (g__11 * g__22 - g__12 * g__21) + g__12 * g__23 * g__31 + g__13
438 : : * (g__21 * g__32 - g__22 * g__31) - g__11 * g__23 * g__32, 0.2e1),
439 : : -0.2e1 / 0.3e1) / ((-g__11 * g__23 + g__13 * g__21) * g__32
440 : : + (g__11 * g__33 - g__13 * g__31) * g__22 - g__12 * (g__21
441 : : * g__33 - g__23 * g__31));
442 : :
443 : : // d(sigma_31)/d(g_11)
444 : : dsigdg[2][0] = - ((g__32 * (3 * g__12 * g__12 + 4 * g__22 * g__22)
445 : : * g__21 + g__22 * ((-4 * g__12 * g__12 - 4 * g__22 * g__22) * g__31
446 : : + g__11 * g__12 * g__32)) * g__33 * g__33 + ((g__12 * (g__22 * g__22
447 : : - 6 * g__32 * g__32) * g__21 - g__22 * (-7 * g__12 * g__31 * g__32
448 : : + g__11 * (g__22 * g__22 + g__32 * g__32))) * g__13
449 : : + (-g__22 * (g__12 * g__12 + 8 * g__32 * g__32) * g__21 + g__32
450 : : * (g__12 * g__12 + 8 * g__22 * g__22) * g__31 + g__11 * g__12
451 : : * (g__22 - g__32) * (g__22 + g__32)) * g__23) * g__33
452 : : + 3 * (g__22 * g__22 + g__32 * g__32) * (g__21 * g__32 - g__22 * g__31)
453 : : * g__13 * g__13 + g__23 * (-7 * g__12 * g__22 * g__32 * g__21
454 : : + (6 * g__12 * g__22 * g__22 - g__12 * g__32 * g__32) * g__31
455 : : + g__32 * g__11 * (g__22 * g__22 + g__32 * g__32)) * g__13 - g__23
456 : : * g__23 * ((-4 * g__12 * g__12 * g__32 - 4 * std::pow( g__32, 3))
457 : : * g__21 + g__22 * (g__31 * (3 * g__12 * g__12 + 4 * g__32 * g__32)
458 : : + g__11 * g__12 * g__32))) * m_mu * std::pow( std::pow( (g__33
459 : : * (g__11 * g__22 - g__12 * g__21) + g__12 * g__23 * g__31 + g__13
460 : : * (g__21 * g__32 - g__22 * g__31) - g__11 * g__23 * g__32), 2),
461 : : -0.2e1 / 0.3e1) / (g__33 * (g__11 * g__22 - g__12 * g__21) + g__13
462 : : * (g__21 * g__32 - g__22 * g__31) - g__23 * (g__32 * g__11 - g__12
463 : : * g__31)) / 0.3e1;
464 : :
465 : : // d(sigma_31)/d(g_21)
466 : : dsigdg[2][1] = 0.4e1 / 0.3e1 * (((-std::pow(g__12, 0.3e1)
467 : : - g__12 * g__22 * g__22) * g__31 + ((g__12 * g__12 + 0.3e1 / 0.4e1
468 : : * g__22 * g__22) * g__11 + g__12 * g__21 * g__22 / 0.4e1) * g__32)
469 : : * g__33 * g__33 + ((0.7e1 / 0.4e1 * g__12 * g__22 * g__31 * g__32
470 : : + g__22 * (g__12 * g__12 - 0.6e1 * g__32 * g__32) * g__11 / 0.4e1
471 : : - g__12 * g__21 * (g__12 * g__12 + g__32 * g__32) / 0.4e1) * g__23
472 : : - g__13 * ((-0.8e1 * g__12 * g__12 - g__22 * g__22) * g__32 * g__31
473 : : + g__12 * (g__22 * g__22 + 0.8e1 * g__32 * g__32) * g__11 - g__21
474 : : * g__22 * (g__12 - g__32) * (g__12 + g__32)) / 0.4e1) * g__33
475 : : + 0.3e1 / 0.4e1 * (g__12 * g__12 + g__32 * g__32) * (g__32 * g__11
476 : : - g__12 * g__31) * g__23 * g__23 - 0.7e1 / 0.4e1 * g__13
477 : : * (-0.6e1 / 0.7e1 * g__22 * (g__12 * g__12 - g__32 * g__32 / 0.6e1)
478 : : * g__31 + g__32 * (g__11 * g__12 * g__22 - g__21 * (g__12 * g__12
479 : : + g__32 * g__32) / 0.7e1)) * g__23 + g__13 * g__13 * ((-0.3e1 / 0.4e1
480 : : * g__12 * g__22 * g__22 - g__12 * g__32 * g__32) * g__31
481 : : + (g__11 * (g__22 * g__22 + g__32 * g__32) - g__12 * g__21
482 : : * g__22 / 0.4e1) * g__32)) * m_mu * std::pow(std::pow(g__33
483 : : * (g__11 * g__22 - g__12 * g__21) + g__12 * g__23 * g__31 + g__13
484 : : * (g__21 * g__32 - g__22 * g__31) - g__11 * g__23 * g__32, 0.2e1),
485 : : -0.2e1 / 0.3e1) / (g__13 * (g__21 * g__32 - g__22 * g__31)
486 : : + (-g__32 * g__11 + g__12 * g__31) * g__23 + g__33
487 : : * (g__11 * g__22 - g__12 * g__21));
488 : :
489 : : // d(sigma_31)/d(g_31)
490 : : dsigdg[2][2] = -(((-0.4e1 / 0.3e1 * g__12 * g__32 * g__32
491 : : - 0.4e1 / 0.3e1 * std::pow(g__12, 0.3e1)) * g__21 + 0.4e1 / 0.3e1
492 : : * g__22 * ((g__12 * g__12 + 0.3e1 / 0.4e1 * g__32 * g__32) * g__11
493 : : + g__12 * g__31 * g__32 / 0.4e1)) * g__23 * g__23 + ((0.7e1 / 0.3e1
494 : : * g__12 * g__22 * g__32 * g__21 + g__32 * (g__12 * g__12 - 0.6e1
495 : : * g__22 * g__22) * g__11 / 0.3e1 - g__12 * g__31 * (g__12 * g__12
496 : : + g__22 * g__22) / 0.3e1) * g__33 - 0.8e1 / 0.3e1 * g__13 * (-g__22
497 : : * (g__12 * g__12 + g__32 * g__32 / 0.8e1) * g__21 + g__12 * (g__22
498 : : * g__22 + g__32 * g__32 / 0.8e1) * g__11 - g__31 * g__32 * (g__12
499 : : - g__22) * (g__12 + g__22) / 0.8e1)) * g__23 + (g__12 * g__12
500 : : + g__22 * g__22) * (g__11 * g__22 - g__12 * g__21) * g__33 * g__33
501 : : - 0.7e1 / 0.3e1 * (-0.6e1 / 0.7e1 * g__32 * (g__12 * g__12 - g__22
502 : : * g__22 / 0.6e1) * g__21 + (g__11 * g__12 * g__32 - g__31 * (g__12
503 : : * g__12 + g__22 * g__22) / 0.7e1) * g__22) * g__13 * g__33
504 : : + 0.4e1 / 0.3e1 * (-g__12 * (g__22 * g__22 + 0.3e1 / 0.4e1
505 : : * g__32 * g__32) * g__21 + g__22 * (g__11 * (g__22 * g__22
506 : : + g__32 * g__32) - g__12 * g__31 * g__32 / 0.4e1)) * g__13 * g__13)
507 : : * m_mu * std::pow(std::pow(g__33 * (g__11 * g__22 - g__12 * g__21)
508 : : + g__12 * g__23 * g__31 + g__13 * (g__21 * g__32 - g__22 * g__31)
509 : : - g__11 * g__23 * g__32, 0.2e1), -0.2e1 / 0.3e1) / (g__13 * (g__21
510 : : * g__32 - g__22 * g__31) + (-g__32 * g__11 + g__12 * g__31) * g__23
511 : : + g__33 * (g__11 * g__22 - g__12 * g__21));
512 : :
513 : : // Define amat
514 : : double amat[9];
515 : : for (std::size_t i=0; i<3; ++i)
516 : : for (std::size_t j=0; j<3; ++j)
517 : : {
518 : : amat[i*3+j] = 0.0;
519 : : for (std::size_t k=0; k<3; ++k)
520 : : {
521 : : amat[i*3+j] -= defgrad[k][j]*dsigdg[i][k];
522 : : }
523 : : amat[i*3+j] /= (arho/alpha);
524 : : }
525 : :
526 : : double eig_real[3], eig_imag[3];
527 : : double vl[3], vr[3];
528 : : #ifndef NDEBUG
529 : : lapack_int ierr =
530 : : #endif
531 : : LAPACKE_dgeev (LAPACK_ROW_MAJOR, 'N', 'N', 3, amat, 3,
532 : : eig_real, eig_imag, vl, 3, vr, 3);
533 : : Assert(ierr==0, "Lapack failed to compute eigenvalues");
534 : :
535 : : // Manually find max
536 : : tk::real eig_max = eig_real[0];
537 : : for (std::size_t i=1; i<3; i++)
538 : : if (eig_real[i] > eig_max)
539 : : eig_max = eig_real[i];
540 : :
541 : : */
542 : :
543 : : // Approximated elastic contribution, from Barton, P. T. (2019).
544 : : // An interface-capturing Godunov method for the simulation of compressible
545 : : // solid-fluid problems. Journal of Computational Physics, 390, 25-50
546 : 0 : tk::real a = (4.0/3.0) * m_mu * alpha / arho;
547 : :
548 : : // hydrodynamic contribution
549 : 0 : auto p_eff = std::max( 1.0e-15, apr+(alpha*m_pstiff) );
550 : 0 : a += m_gamma * p_eff / arho;
551 : :
552 : : // Compute square root
553 : 0 : a = std::sqrt(a);
554 : :
555 : : // check sound speed divergence
556 [ - - ]: 0 : if (!std::isfinite(a)) {
557 : 0 : std::cout << "Material-id: " << imat << std::endl;
558 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
559 : 0 : std::cout << "Partial density: " << arho << std::endl;
560 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
561 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf sound speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
562 : : + std::to_string(a) + ", material volume fraction: " +
563 : : std::to_string(alpha));
564 : : }
565 : :
566 : 0 : return a;
567 : : }
568 : :
569 : : tk::real
570 : 0 : SmallShearSolid::shearspeed(
571 : : tk::real arho,
572 : : tk::real alpha,
573 : : std::size_t imat ) const
574 : : // *************************************************************************
575 : : //! Calculate speed of sound from the material density and material pressure
576 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
577 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
578 : : //! the single-material system, this argument can be left unspecified by
579 : : //! the calling code
580 : : //! \param[in] imat Material-id who's EoS is required. Default is 0, so that
581 : : //! for the single-material system, this argument can be left unspecified
582 : : //! by the calling code
583 : : //! \return Material shear-wave speed speed using the SmallShearSolid EoS
584 : : // *************************************************************************
585 : : {
586 : : // Approximate shear-wave speed. Ref. Barton, P. T. (2019).
587 : : // An interface-capturing Godunov method for the simulation of compressible
588 : : // solid-fluid problems. Journal of Computational Physics, 390, 25-50.
589 : 0 : tk::real a = std::sqrt(alpha*m_mu/arho);
590 : :
591 : : // check shear-wave speed divergence
592 [ - - ]: 0 : if (!std::isfinite(a)) {
593 : 0 : std::cout << "Material-id: " << imat << std::endl;
594 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
595 : 0 : std::cout << "Partial density: " << arho << std::endl;
596 [ - - ][ - - ]: 0 : Throw("Material-" + std::to_string(imat) + " has nan/inf shear-wave speed: "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
597 : : + std::to_string(a) + ", material volume fraction: " +
598 : : std::to_string(alpha));
599 : : }
600 : :
601 : 0 : return a;
602 : : }
603 : :
604 : : tk::real
605 : 0 : SmallShearSolid::totalenergy(
606 : : tk::real rho,
607 : : tk::real u,
608 : : tk::real v,
609 : : tk::real w,
610 : : tk::real pr,
611 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
612 : : // *************************************************************************
613 : : //! \brief Calculate material specific total energy from the material
614 : : //! density, momentum and material pressure
615 : : //! \param[in] rho Material density
616 : : //! \param[in] u X-velocity
617 : : //! \param[in] v Y-velocity
618 : : //! \param[in] w Z-velocity
619 : : //! \param[in] pr Material pressure
620 : : //! \param[in] defgrad Material inverse deformation gradient tensor
621 : : //! g_k. Default is 0, so that for the single-material system,
622 : : //! this argument can be left unspecified by the calling code
623 : : //! \return Material specific total energy using the SmallShearSolid EoS
624 : : // *************************************************************************
625 : : {
626 : : // obtain hydro contribution to energy
627 : 0 : tk::real rhoEh = (pr + m_pstiff) / (m_gamma-1.0) + 0.5 * rho *
628 : 0 : (u*u + v*v + w*w) + m_pstiff;
629 : : // obtain elastic contribution to energy
630 : : tk::real eps2;
631 [ - - ]: 0 : tk::real rhoEe = elasticEnergy(defgrad, eps2);
632 : :
633 : 0 : return (rhoEh + rhoEe);
634 : : }
635 : :
636 : : tk::real
637 : 0 : SmallShearSolid::temperature(
638 : : tk::real arho,
639 : : tk::real u,
640 : : tk::real v,
641 : : tk::real w,
642 : : tk::real arhoE,
643 : : tk::real alpha,
644 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad ) const
645 : : // *************************************************************************
646 : : //! \brief Calculate material temperature from the material density, and
647 : : //! material specific total energy
648 : : //! \param[in] arho Material partial density (alpha_k * rho_k)
649 : : //! \param[in] u X-velocity
650 : : //! \param[in] v Y-velocity
651 : : //! \param[in] w Z-velocity
652 : : //! \param[in] arhoE Material total energy (alpha_k * rho_k * E_k)
653 : : //! \param[in] alpha Material volume fraction. Default is 1.0, so that for
654 : : //! the single-material system, this argument can be left unspecified by
655 : : //! the calling code
656 : : //! \param[in] defgrad Material inverse deformation gradient tensor
657 : : //! (g_k). Default is 0, so that for the single-material system,
658 : : //! this argument can be left unspecified by the calling code
659 : : //! \return Material temperature using the SmallShearSolid EoS
660 : : // *************************************************************************
661 : : {
662 : : // obtain elastic contribution to energy
663 : : tk::real eps2;
664 [ - - ]: 0 : auto arhoEe = alpha*elasticEnergy(defgrad, eps2);
665 : : // obtain hydro contribution to energy
666 : 0 : auto arhoEh = arhoE - arhoEe;
667 : :
668 : 0 : tk::real t = (arhoEh - 0.5 * arho * (u*u + v*v + w*w) - alpha*m_pstiff)
669 : 0 : / (arho*m_cv);
670 : 0 : return t;
671 : : }
672 : :
673 : : tk::real
674 : 0 : SmallShearSolid::min_eff_pressure(
675 : : tk::real min,
676 : : tk::real,
677 : : tk::real ) const
678 : : // *************************************************************************
679 : : //! Compute the minimum allowed pressure
680 : : //! \param[in] min Numerical threshold above which pressure needs to be limited
681 : : //! \return Minimum pressure allowed by physical constraints
682 : : // *************************************************************************
683 : : {
684 : : // minimum pressure is constrained by zero soundspeed.
685 : 0 : return (min - m_pstiff);
686 : : }
687 : :
688 : : tk::real
689 : 0 : SmallShearSolid::elasticEnergy(
690 : : const std::array< std::array< tk::real, 3 >, 3 >& defgrad,
691 : : tk::real& eps2 ) const
692 : : // *************************************************************************
693 : : //! \brief Calculate elastic contribution to material energy from the material
694 : : //! density, and deformation gradient tensor
695 : : //! \param[in] defgrad Material inverse deformation gradient tensor
696 : : //! \param[in/out] eps2 Elastic shear distortion
697 : : //! \return Material elastic energy using the SmallShearSolid EoS
698 : : //! \details This function returns the material elastic energy, and also stores
699 : : //! the elastic shear distortion for further use
700 : : // *************************************************************************
701 : : {
702 : : // compute volume-preserving part of Right Cauchy-Green strain tensor
703 [ - - ]: 0 : auto Ct = tk::getIsochorRightCauchyGreen(defgrad);
704 : :
705 : : // compute elastic shear distortion
706 : 0 : eps2 = 0.5 * (Ct[0][0]+Ct[1][1]+Ct[2][2] - 3.0);
707 : :
708 : : // compute elastic energy
709 : 0 : auto rhoEe = m_mu * eps2;
710 : :
711 : 0 : return rhoEe;
712 : : }
|
