Tensile Shock Physics in Compressible Thermoviscoelastic Solid Medium

Karan S. Surana, Elie Abboud

Department of Mechanical Engineering, University of Kansas, Lawrence, USA.

DOI: 10.4236/am.2024.1510042                

Abstract

This paper addresses tensile shock physics in thermoviscoelastic (TVE) solids without memory. The mathematical model is derived using conservation and balance laws (CBL) of classical continuum mechanics (CCM), incorporating the contravariant second Piola-Kirchhoff stress tensor, the covariant Green’s strain tensor, and its rates up to order n. This mathematical model permits the study of finite deformation and finite strain compressible deformation physics with an ordered rate dissipation mechanism. Constitutive theories are derived using conjugate pairs in entropy inequality and the representation theorem. The resulting mathematical model is both thermodynamically and mathematically consistent and has closure. The solution of the initial value problems (IVPs) describing evolutions is obtained using a variationally consistent space-time coupled finite element method, derived using space-time residual functional in which the local approximations are in hpk higher-order scalar product spaces. This permits accurate description problem physics over the discretization and also permits precise a posteriori computation of the space-time residual functional, an accurate measure of the accuracy of the computed solution. Model problem studies are presented to demonstrate tensile shock formation, propagation, reflection, and interaction. A unique feature of this research is that tensile shocks can only exist in solid matter, as their existence requires a medium to be elastic (presence of strain), which is only possible in a solid medium. In tensile shock physics, a decrease in the density of the medium caused by tensile waves leads to shock formation ahead of the wave. In contrast, in compressive shocks, an increase in density and the corresponding compressive waves result in the formation of compression shocks behind of the wave. Although these are two similar phenomena, they are inherently different in nature. To our knowledge, this work has not been reported in the published literature.

Keywords

Tensile Shock Physics, Tensile Waves, Elastic, Viscoelastic Solids, Variationally Consistent, Space-Time Coupled, Space-Time Residual Functional, A Posteriori, Finite Element Method, Wave Speed, Conservation and Balance Laws

1. Introduction, Literature Review, and Scope of Work

1.1. Introduction and Literature Review

In recent papers [1] [2], Surana et al. presented shock physics in compressible thermoelastic (TE) and thermoviscoelastic (TVE) solid media without rheology, as well as shock physics in TVE with rheology. In these works, compressive stress pulse or velocity pulse of duration 2∆t is applied in one dimensional studies of pure 1D deformation of axial rod. It was shown that, along the base of the wave (i.e. along the length of the rod containing the wave), there is a continuous increase in density as the stress magnitude increases, reaching a peak increase in density corresponding to the peak value of the stress. This is followed by a continuous decrease in density upstream of the peak of the stress wave. For fixed elasticity, this change in density results in a continuous change in wave speed, keeping in mind that increasing density results in a reduction of wave speed, and vice versa. These studies focused on shock physics due to compressive loading always resulting in increase in density, hence reduced wave speed. Continuously varying density over compressive velocity pulse was used to study compression shock physics.

In compressible fluids, shock physics is also volumetric, but pressure waves cannot exist without an equation of state. The equation of state causes volumetric; hence density changes. The piling up of the faster moving compression waves behind the peak results in shock formation. The Riemann shock tube [3] [4] is an example of this physics. Since fluid medium lacks elasticity, thus lacks strain, it cannot support tensile shock physics, in which the external tensile stimulus must result in a decrease in density of the medium. This aspect of tensile shock physics is unique to solid matter and is quite different in terms of physics and its consequences on the resulting nature of the shocks.

For example, a propagating tensile stress wave will result in a continuous change in density over the support of the wave (as it does in the case of compressive stress waves [1]), resulting in a continuous change in wave speeds over the support of the wave when the elastic constants remain unchanged. Consider the evolution of deviatoric stress wave (_dσ₁₁) and density ρ in 1D wave propagation of a tensile stress pulse of duration 2∆t in a compressible TVE solid medium. The dimensionless wave speed is unity. At time 2∆t, the tensile stress pulse is completely in the solid medium, as shown in Figure 1(a) (wave ABC). This causes a change in density, as shown in Figure 1(b) (abc). Density decreases from its reference value ρ₀ to its lowest value at b, corresponding to the peak stress at B shown in Figure 1(b), and then increases from its lowest value at b back to the reference value ρ₀ at c (following the change in stress from B to C in Figure 1(a)). Thus, from A to B, the wave speed continuously increases from its reference value (lowest value) at A (or a) to B (or b), its highest value, and then decreases from B (or b) back to C (or c) to its reference value. Therefore, from A to B, the wave speed is continuously increasing, whereas from B to C, it is continuously decreasing. Thus, during further evolution in the portion AB of the wave, as we move from A to B, wave speed is increasing, i.e., points closer to A are moving slower compared to points closer to B, resulting in shallowing of the portion AB of the wave (reduced slope). The opposite happens in the region BC. From B to C, wave speed is decreasing, meaning that the part of the wave closer to B is moving faster than the part of the wave closer to C. This results in piling up of the waves (as in the Riemann shock tube [3] [4]), causing steepening or increased slope of the portion BC of the wave. We observe from Figure 1 that at t > 2∆t (red color), shock formation occurs ahead of the peak of the wave (portion BC) and rarefaction or shallowing of the wave occurs behind the peak of the wave. When comparing these results with compressive shock physics studies in ref. [1], this phenomenon is exactly opposite of what is observed in [1], where shock formation occurs behind the peak of the wave.

Figure 1. Schematic of 1D deviatoric tensile stress (_dσ₁₁) wave propagation: Evolution of (_dσ₁₁) and density ρ.

Unfortunately, there is not much published work on wave propagation in compressible TE or TVE solids. The recent works published by the author [1] [2] on compressive shock waves in compressible TVE solids with and without memory cite some references [5]-[12] that do not exactly address the same physics as in references [1] [2], but are cited here to provide a sample of the published works. Detailed discussion of these references is of little value to the subject matter considered in this paper, as these references only bear some resemblance to our work in the sense that they do consider wave propagation but mostly in incompressible solid medium.

1.2. Scope of Work

The work presented in this paper investigates tensile shock physics of stress and density waves in compressible thermoviscoelastic solid medium without memory under tensile loading. Tensile shocks of stress, hence density, are formed due to variations in stress and density along the base of the applied velocity or stress pulse. Investigation of formation, propagation, reflection, interaction and transmission of stress and density waves with shock fronts under tensile loading in TVE solids without memory is the main focus of this work. The mathematical model consists of conservation and balance laws of CCM for finite deformation, finite strain deformation physics derived using the contravariant second Piola-Kirchhoff stress tensor and Green’s strain tensor in Lagrangian description. The second Piola-Kirchhoff stress tensor is additively decomposed into equilibrium and deviatoric stress tensors ([13] [14]) to address volumetric and distortional de- formation physics (mutually exclusive) using constitutive theories for equilibrium and deviatoric stress tensors. The constitutive theory for the equilibrium contravariant second Piola-Kirchhoff stress tensor is derived using contravariant equilibrium Cauchy stress tensor for which constitutive theory is derived using Helmholtz free energy density in Eulerian description [13] [14]. The constitutive theory for the deviatoric contravariant second Piola-Kirchhoff stress tensor is derived using the representation theorem [15]-[26]. Green’s strain tensor and its rates up to order n are considered as argument tensors of the deviatoric contravariant second Piola-Kirchhoff stress tensor, the constitutive tensor. Each rate of Green’s strain tensor up to order n provides ordered rate dissipation mechanism.

This mathematical model is thermodynamically and mathematically consistent. In the model problem studies, we intentionally consider one-dimensional stress wave propagation, thus keeping the physics simple so that intricate aspects of the deviatoric stress and density wave propagation, shock formation, reflection of the waves containing shocks, and their interactions can be clearly illustrated and discussed.

Solutions of the mathematical models are obtained using a space-time coupled finite element method based on the space-time residual functional for a space-time strip with time marching [27]-[30]. The space-time integral form from this approach yields unconditionally stable computational processes. The local approximations for the space-time elements are p-version hierarchical in higher order scalar product $H_{{}_{\overline{\text{Ω}}}}^{\left(p,k\right)}$ spaces in which $p=\left(p_1,p_2\right)$ and $k=\left(k_1,k_2\right)$ , in space and time, ensuring that the local approximations in space and time yield the desired higher-order global differentiability with the desired degree of polynomials. In this computational methodology with minimally conforming approximation spaces [27], accurate a posteriori computation of the L₂ norm of the space-time residual is possible. The proximity of this norm to zero is a measure of the accuracy of the solution. Generally, computed solutions with this norm of the order of $O\left({10}^{-6}\right)$ or lower yield sufficiently converged solutions. Only upon obtaining a converged solution for the current space-time strip, the solution is time-marched to the next space-time strip. This approach ensures a converged solution for each space-time strip, hence converged solutions for the entire evolution.

We point out that compressibility in solids, in the Lagrangian description, is governed by the deformation gradient tensor $J$ . This is evident from the conservation of mass (CM), where ${\rho }_0\left(x\right)=\left|J\right|\rho \left(x,t\right)$ . However, the density change due to CM also induces a corresponding pressure field, which is defined by the equation of state and is transformed into the equilibrium second Piola-Kirchhoff stress tensor. Thus, in solids (in the Lagrangian description), the equation of state is needed as it describes the pressure field created due to the density change by CM. Its presence is necessary for the correct force balance in BLM. In the absence of the equation of state, the force balance in BLM is incorrect as it does not account for the force vector due to the pressure field corresponding to the density change due to CM, which would result in incorrect evolution. In conclusion, while it is possible to compute the evolution in compressible solids without using the equation of state (since the equation of state does not directly control compressibility), the computed evolution would be in error because the force balance in the balance of linear momentum (BLM) would be incorrect, ultimately affecting the entire evolution.

2. Mathematical Model

For finite deformation, finite strain physics in TVE compressible solid matter, we consider CBL of CCM derived using the contravariant second Piola-Kirchhoff stress tensor ${\sigma }^{\left[0\right]}$ and the covariant Green’s strain tensor ${\epsilon }_{\left[0\right]}$ in Lagrangian description. The dissipation mechanism is described by the time derivatives of Green’s strain tensor up to orders $n$ (${\epsilon }_{\left[i\right]};i=1,2,\cdots ,n$ ). ${\epsilon }_{\left[0\right]}$ and ${\epsilon }_{\left[i\right]};i=1,2,\cdots ,n$ being argument tensors of ${\sigma }^{\left[0\right]}$ in addition to temperature $\theta$ . Conservation of mass (CM), balance of linear momentum (BLM), balance of angular momentum (BAM), and the first and second laws of thermodynamics are given by [13] [14]:

${\rho }_0\left(x\right)=\left|J\right|\rho \left(x,t\right)$ (1)

${\rho }_0\frac{{\partial }^2\left\{u\right\}}{\partial t^2}-{\rho }_0\left\{F^b\right\}-\left(\left[J\right]\left[{\sigma }^{\left[0\right]}\right]\right)\{\nabla \}=0$ (2)

${ϵ}_{ijk}{\sigma }_{ij}^{\left(0\right)}=0$ (3)

${\rho }_0{c}_v\frac{De}{Dt}-\nabla \cdot q-{\sigma }^{\left[0\right]}:{\dot{\epsilon }}_{\left[0\right]}=0$ (4)

${\rho }_0\left(\frac{D\varphi }{Dt}-\eta \frac{D\theta }{Dt}\right)-{\sigma }^{\left[0\right]}:{\dot{\epsilon }}_{\left[0\right]}+\frac{q\cdot g}{\theta }\le 0$ (5)

where ${\rho }_0\left(x\right)$ is the density at a material point in the reference (or initial) configuration, $\rho \left(x,t\right)$ is the density of the material point at $x_i$ in the current configuration, $\left[J\right]$ is the deformation gradient tensor, $u$ are displacements in the fixed x-frame, $F^b$ is the body force vector per unit mass, ${\sigma }^{\left[0\right]}$ is the symmetric contravariant second Piola-Kirchhoff stress tensor, $\nabla$ is the gradient operator, ${ϵ}_{ijk}$ is the permutation tensor, $c_v$ is the specific heat, $e$ is specific internal energy density, $q$ is the heat vector, $\varphi$ is the Helmholtz free energy density, $\eta$ is the entropy density, $\theta$ is the absolute temperature, and $g$ is the temperature gradient vector. We consider additive decomposition of the stress tensor ${\sigma }^{\left[0\right]}$ into equilibrium (${}_e\sigma {}^{\left[0\right]}$ ) and deviatoric (${}_d\sigma {}^{\left[0\right]}$ ) stress tensors:

${\sigma }^{\left[0\right]}={}_e\sigma {}^{\left[0\right]}+{}_d\sigma {}^{\left[0\right]}$ (6)

The constitutive theory for ${}_e\sigma {}^{\left[0\right]}$ is related to volumetric deformation and the constitutive theory for ${}_d\sigma {}^{\left[0\right]}$ describes distortional deformation physics. Following references [13] [14], the constitutive theory for ${}_e\sigma {}^{\left[0\right]}$ can be derived using thermodynamic pressure $p\left(\rho ,\theta \right)$ , and we can write:

$\left[{}_e{\sigma }^{\left[0\right]}\right]=\left|J\right|p\left(\rho ,\theta \right){\left[{\left[J\right]}^{\text{T}}\left[J\right]\right]}^{-1}$ (7)

Thus, we note that ${}_e\sigma {}^{\left[0\right]}$ is not a pressure field. The constitutive theory for ${}_d\sigma {}^{\left[0\right]}$ can be derived using $\theta$ , $\left[{\epsilon }_{\left[0\right]}\right]$ , $\left[{\epsilon }_{\left[i\right]}\right]$ ; $i=1,2,\cdots ,n$ as argument tensors of ${}_d\sigma {}^{\left[0\right]}$ and the representation theorem [14]. A constitutive theory for ${}_d\sigma {}^{\left[0\right]}$ based on integrity (complete basis consisting of the combined generators of the argument tensors of ${}_d\sigma {}^{\left[0\right]}$ , including $\left[I\right]$ ) can be derived (see reference [13] [14] for details). In the following, we present a simple constitutive theory for ${}_d\sigma {}^{\left[0\right]}$ that is linear in the components of the argument tensors ${\epsilon }_{\left[0\right]}$ and ${\epsilon }_{\left[i\right]}$ ; $i=1,2,\cdots ,n$ .

$\begin{array}{c}\left[{}_d{\sigma }^{\left[0\right]}\right]={{\sigma }_0|}_{\Omega }\left[I\right]+{2\eta |}_{\underset{\_}{\Omega }}\left[{\epsilon }_{\left[0\right]}\right]+{\lambda |}_{\underset{\_}{\Omega }}\left(tr\left[{\epsilon }_{\left[0\right]}\right]\right)\left[I\right]\\ +{\displaystyle \sum }_{i=1}^n{2{\eta }_i|}_{\underset{\_}{\Omega }}\left[{\epsilon }_{\left[i\right]}\right]+{\displaystyle \sum }_{i=1}^n{{\kappa }_i|}_{\underset{\_}{\Omega }}\left(tr\left[{\epsilon }_{\left[i\right]}\right]\right)\left[I\right]-{\alpha }_{tm}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)\left[I\right]\end{array}$ (8)

In (8), product terms have not been considered as they are quadratic, hence of degree higher than one. $\underset{\_}{\Omega }$ refers to a known configuration in which the material coefficients $\eta$ , $\lambda$ , ${\eta }_i$ , ${\kappa }_i$ , and ${\alpha }_{tm}$ can be evaluated if they are functions of the invariants of the argument tensors of ${}_d\sigma {}^{\left[0\right]}$ and temperature $\theta$ . If we further limit the dependence of ${}_d\sigma {}^{\left[0\right]}$ in (8) to $\left[{\epsilon }_{\left[0\right]}\right]$ , $tr\left[{\epsilon }_{\left[0\right]}\right]$ , $\left[{\epsilon }_{\left[1\right]}\right]$ and $tr\left[{\epsilon }_{\left[1\right]}\right]$ , neglect initial stress ${{\sigma }_0|}_{\underset{\_}{\Omega }}$ and the thermal expansion term (last term in (8)), then we can obtain a further simplified but more practical and usable constitutive theory for ${}_d{\sigma }^{\left[0\right]}$ (we drop $\underset{\_}{\Omega }$ for convenience of notation):

$\left[{}_d{\sigma }^{\left[0\right]}\right]=2\eta \left[{\epsilon }_{\left[0\right]}\right]+\lambda \left(tr\left[{\epsilon }_{\left[0\right]}\right]\right)\left[I\right]+2{\eta }_1\left[{\epsilon }_{\left[1\right]}\right]+{\kappa }_1\left(tr\left[{\epsilon }_{\left[1\right]}\right]\right)\left[I\right]$ (9)

in which

$\left[{\epsilon }_{\left[0\right]}\right]=\frac{1}{2}\left({\left[J\right]}^{\text{T}}\left[J\right]-\left[I\right]\right),\left[J\right]=\left[{}^{d}J\right]+\left[I\right]\text{and}\left[{}^{d}J\right]=\left[\frac{\partial \left\{u\right\}}{\partial \left\{x\right\}}\right]$ (10)

with $\left[{}^{d}J\right]$ being the displacement gradient tensor and

$\left[{\epsilon }_{\left[1\right]}\right]=\frac{D}{Dt}\left[{\epsilon }_{\left[0\right]}\right]=\frac{\partial }{\partial t}\left[{\epsilon }_{\left[0\right]}\right]=\left({\left[\dot{J}\right]}^{\text{T}}\left[J\right]+{\left[J\right]}^{\text{T}}\left[\dot{J}\right]\right)$ (11)

$\left[\dot{J}\right]$ being the velocity gradient tensor in Lagrangian description. We can assume simple Fourier heat conduction law as the constitutive theory for $q$ (see reference [13] [14] for constitutive theory based on integrity):

$\left\{q\right\}=-k\left\{g\right\}$ (12)

This mathematical model ((2), (3), (4), (9), and (12)) consists of 13 equations: BLM (3), FLT (1), constitutive theories for ${}_d\sigma {}^{\left[0\right]}$ (6) and $q$ (3) in thirteen variables: $u$ (3), $\theta$ (1), ${}_d\sigma {}^{\left[0\right]}$ (6) and $q$ (3), hence the mathematical model has closure.

This mathematical model can be used to study 3D compressible deformation physics of TVE solids without memory. To illustrate the shock physics of deviatoric stress waves and density more clearly, it is necessary to consider physics in which dispersion and boundary effects causing transmission and reflection do not contaminate the solution. For this reason, we consider 1D deviatoric contravariant second Piola-Kirchhoff stress wave propagation in compressible thermoviscoelastic medium without memory. Furthermore, if we consider an insulated system, then the entropy generation is only due to viscous dissipation, which leads to extremely small temperature changes (for non-cyclic loads) that has virtually no effect on mechanical deformation. Thus, with this assumption, we can eliminate the energy equation (4). The entropy inequality has already been used in deriving constitutive theories for ${}_e\sigma {}^{\left[0\right]}$ , ${}_d\sigma {}^{\left[0\right]}$ , and $q$ . The reduced form of entropy inequality (5) is given by:

$-{}_d\sigma {}^{\left[0\right]}:{\dot{\epsilon }}_{\left[0\right]}+\frac{q\cdot g}{\theta }\le 0$ (13)

As mentioned earlier, in the Lagrangian description, CM is not part of the mathematical model as $\rho \left(x,t\right)$ in the current configuration is deterministic using CM (1). Thus, for the 1D case, we have displacement $u_1$ and deviatoric contravariant second Piola-Kirchhoff stress ${}_d{\sigma }_{11}^{\left[0\right]}$ as the only two dependent variables in the balance of linear momenta in the $x_1$ direction and the constitutive theory for ${}_d{\sigma }_{11}^{\left[0\right]}$ defined using (9). Thermodynamic pressure is only a function of $\rho$ and is known when $\left[J\right]$ is known (Equation (7)).

For 1D finite strain, finite deformation we have:

$\left[J\right]=\left|J\right|=\left(1+\frac{\partial u_1}{\partial x_1}\right)$ (14)

${\left[{\epsilon }_{\left[0\right]}\right]}_{11}=\frac{\partial u_1}{\partial x_1}+\frac{1}{2}{\left(\frac{\partial u_1}{\partial x_1}\right)}^2$ (15)

${\left[{\epsilon }_{\left[1\right]}\right]}_{11}=\frac{\partial {\left[{\epsilon }_{\left[0\right]}\right]}_{11}}{\partial t}=\frac{\partial v_1}{\partial x_1}+\frac{\partial u_1}{\partial x_1}\frac{\partial v_1}{\partial x_1}$ (16)

where

$\frac{\partial v_1}{\partial x_1}=\frac{\partial }{\partial t}\left(\frac{\partial u_1}{\partial x_1}\right)$ (17)

Using (14)-(17) and (6), the mathematical model for 1D wave propagation in a compressible TVE medium without memory can be written as follows. We have changed the sign of the second term in (18) by assuming that compressive pressure and the corresponding ${}_e{\sigma }_{11}^{\left[0\right]}$ is positive when compressive.

${\rho }_0\frac{{\partial }^2{u}_1}{\partial t^2}-{\rho }_0{F}_1^b+\frac{\partial }{\partial x_1}\left[\left(1+\frac{\partial u_1}{\partial x_1}\right){}_e\sigma {}_{11}^{\left[0\right]}\right]-\frac{\partial }{\partial x_1}\left[\left(1+\frac{\partial u_1}{\partial x_1}\right){}_d\sigma {}_{11}^{\left[0\right]}\right]=0$ (18)

${}_e{\sigma }_{11}^{\left[0\right]}=\frac{p\left(\rho \right)}{\left(1+\frac{\partial u_1}{\partial x_1}\right)}$ (19)

${}_d{\sigma }_{11}^{\left[0\right]}=C_1\left(\frac{\partial u_1}{\partial x_1}+\frac{1}{2}{\left(\frac{\partial u_1}{\partial x_1}\right)}^2\right)+C_2\left(\frac{\partial v_1}{\partial x_1}+\frac{\partial u_1}{\partial x_1}\frac{\partial v_1}{\partial x_1}\right)$ (20)

and

$v_1=\frac{\partial u_1}{\partial t}\forall \left(x,t\right)\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t$ (21)

in which $C_1=2\eta +\lambda$ and $C_2=2{\eta }_1+{\kappa }_1$ . $C_1$ and $C_2$ are material coefficients for elasticity and dissipation. Further details on the equilibrium contravariant second Piola-Kirchhoff stress tensor are given in the following.

In (19), we choose simple equation of state $p\left(\rho \right)$ given in the following:

$p\left(\rho \right)=C\left(\frac{\rho }{{\rho }_{{}_{{}_{\text{}\text{}0}}}}-1\right)$ (22)

in which $C$ is the bulk modulus. From (22), we note that increasing $\rho$ implies

$\frac{\rho }{{\rho }_0}>1$ , yielding compressive positive pressure $p$ , confirming that (22) has correct physics.

Substituting (22) in (19) and then (19) in the third term of BLM (18), we can obtain the following (assuming $C$ to be constant):

$\frac{\partial }{\partial x_1}\left[\left(1+\frac{\partial u_1}{\partial x_1}\right){}_e\sigma {}^{\left[0\right]}\right]=\frac{\partial }{\partial x_1}p\left(\rho \right)$ (23)

in which $p\left(\rho \right)$ can be written as follows (using CM):

$p\left(\rho \right)=C\left(\frac{\rho }{{\rho }_0}-1\right)=C\left[{\left(1+\frac{\partial u_1}{\partial x_1}\right)}^{-1}-1\right]$ (24)

We substitute (24) in (23) to obtain

$\frac{\partial }{\partial x_1}\left[\left(1+\frac{\partial u_1}{\partial x_1}\right){}_e\sigma {}_{11}^{\left[0\right]}\right]=C\left(-\frac{1}{{\left(1+\frac{\partial u_1}{\partial x_1}\right)}^2}\frac{{\partial }^2{u}_1}{\partial x^2}\right)$ (25)

Equation (25) can be used in the balance of linear momentum (BLM) equation (18) to obtain the final form of the BLM.

The final model consists of the following BLM, constitutive theory for ${}_d{\sigma }_{11}^{\left[0\right]}$ , and velocity $v_1$ and equation of state:

${\rho }_0\frac{\partial v_1}{\partial t}-{\rho }_0{F}_1^b+\left(-C\left[{\left(1+\frac{\partial u_1}{\partial x_1}\right)}^{-2}\frac{{\partial }^2{u}_1}{\partial x_1{}^2}\right]\right)-\frac{{\partial }^2{u}_1}{\partial x^2}{}_d\sigma {}_{11}^{\left[0\right]}-\left(1+\frac{\partial u_1}{\partial x_1}\right)\frac{\partial {}_d\sigma {}^{\left[0\right]}}{\partial x_1}=0$ (26)

${}_d{\sigma }_{11}^{\left[0\right]}=C_1\left(\frac{\partial u_1}{\partial x_1}+\frac{1}{2}{\left(\frac{\partial u_1}{\partial x_1}\right)}^2\right)+C_2\left(\frac{\partial v_1}{\partial x_1}+\frac{\partial u_1}{\partial x_1}\frac{\partial v_1}{\partial x_1}\right)$ (27)

and

$v_1=\frac{\partial u_1}{\partial t}\forall \left(x,t\right)\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t$ (28)

3. Solutions of the PDEs in the Mathematical Model

The system of three partial differential equations in (26)-(28) are systems of nonlinear partial differential equations (PDEs) that constitute initial value problem (IVP) in which dependent variables $u_1$ , ${}_d{\sigma }_{11}^{\left[0\right]}$ , and $v_1$ exhibit simultaneous dependence on space coordinates $x_1$ and time $t$ . We consider a space-time coupled finite element method based on a space-time residual functional [3] [27]-[30] for a space-time strip with time marching for obtaining solutions of the PDEs in the three models. It is shown in reference [27] that the space-time differential operators in the IVPs are either non-self-adjoint or nonlinear, for which all other methods of approximation (Space-Time Galerkin Method (STGM), Space-Time Petrov-Galerkin Method (STPGM), Space-Time Weighted Residual Method (STWRM), Space-Time Galerkin Method/Weak Form (STGM/WF)) yield space-time variationally inconsistent integral forms that do not ensure unconditionally stable computations during the evolution. Only the space-time integral form based on space-time residual functional is space-time variationally consistent, hence yields unconditionally stable computational processes during the entire evolution. Basic steps of the computational process are given in the following (see reference [27] for more details).

Let $t=0$ be the time at which the evolution commences. Let $\text{Δ}t$ be an increment of time. The space-time domain ${\overline{\Omega }}_{xt}^{\left(1\right)}={\overline{\Omega }}_x\times {\overline{\Omega }}_t^{\left(1\right)}=\left[0,L\right]\times \left[0,\Delta t\right]$ , of the first space-time strip shown in Figures 2(a)-(b), is discretized using nine-node p-version hierarchical space-time finite elements in which the local approximations are of higher degree and higher order, yielding higher order global differentiability in space and time, i.e., local approximations are in higher order scalar product space $H_{{}_{\overline{\text{Ω}}}}^{\left(p,k\right)}$, $p=\left(p_1,p_2\right)$ , $k=\left(k_1,k_2\right)$ , where one and two refer to space and time, respectively.

Let ${\left(u_1\right)}_h^e$ , ${\left(v_1\right)}_h^e$ , and ${\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e$ be the local approximations of $u_1$ , $v_1$ , and ${}_d{\sigma }_{11}^{\left[0\right]}$ over space-time element ${\overline{\Omega }}_{xt}^e$ such that

${\left(u_1\right)}_h={\displaystyle \underset{e}{\cup }{\left(u_1\right)}_h^e},{\left(v_1\right)}_h={\displaystyle \underset{e}{\cup }{\left(v_1\right)}_h^e},{\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h={\displaystyle \underset{e}{\cup }{\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e}$ (29)

in which ${\left(u_1\right)}_h$ , ${\left(v_1\right)}_h$ , and ${\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h$ are approximations of $u_1$ , $v_1$ , and ${}_d{\sigma }_{11}^{\left[0\right]}$

over the discretization ${\left({\overline{\Omega }}_{xt}^{\left(1\right)}\right)}^T$ of ${\overline{\Omega }}_{xt}^{\left(1\right)}$ .

Figure 2. 1D solid domain, idealization of 1D solid domain, discretization of first space-time strip with space-time finite elements and applied disturbance.

${\left({\overline{\Omega }}_{xt}^{\left(1\right)}\right)}^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_{xt}^e}$ (30)

The local approximations over a space-time element ${}_{\overline{\text{Ω}}}^{xt}$ can be written as

$\begin{array}{l}{\left(u_1\right)}_h^e={\displaystyle {\sum }_{i=1}^{n_u}{N}_i^u}\left(\xi ,\eta \right){\left\{{\delta }_i^e\right\}}_u\\ {\left(v_1\right)}_h^e={\displaystyle {\sum }_{i=1}^{n_v}{N}_i^v}\left(\xi ,\eta \right){\left\{{\delta }_i^e\right\}}_v\\ {\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e={\displaystyle {\sum }_{i=1}^{n_{\sigma }}{N}_i^{\sigma }}\left(\xi ,\eta \right){\left\{{\delta }_i^e\right\}}_{\sigma }\end{array}$ (31)

where $N_i^u$ , $N_i^v$ , and $N_i^{\sigma }$ are approximation functions for an element ${\overline{\Omega }}_{xt}^e$ mapped in the natural coordinate space $\left(\xi ,\eta \right)$ in a two-unit square, and ${\left\{{\delta }_i^e\right\}}_u$ ,

${\left\{{\delta }_i^e\right\}}_v$ , and ${\left\{{\delta }_i^e\right\}}_{\sigma }$ are the corresponding degrees of freedom for $u$ , $v$ and ${}_d{\sigma }_{11}^{\left[0\right]}$ . Substituting (31) in (26)-(28) yields three space-time residual equations

(functions) $E_1^e$ , $E_2^e$ , and $E_3^e$ $\forall$ $x,t\in {\overline{\Omega }}_{xt}$ or $\xi ,\eta \in {\overline{\Omega }}_{\xi ,\eta }$ . The space-time finite element method based on the residual functional utilizes residual equations $E_1^e$ , $E_2^e$ , and $E_3^e$ over an element ${\overline{\Omega }}_{xt}^e$ . Following references [27], we can write the residual functional for an element $e$ with domain ${\overline{\Omega }}_{xt}^e$ .

$I^e={\displaystyle \sum }_{i=1}^3{\left(E_i^e,E_i^e\right)}_{{\overline{\Omega }}_{xt}^e}$ (32)

and the residual functional for the discretization ${\left({\overline{\Omega }}_{xt}^{\left(1\right)}\right)}^T$ is given by

$I={\displaystyle \sum }_e{I}^e={\displaystyle \sum }_e{\displaystyle \sum }_{i=1}^3{\left(E_i^e,E_i^e\right)}_{{\overline{\Omega }}_{xt}^e}$ (33)

An extremum of $I$ requires that we set $\delta I=0$ (a necessary condition) provided $I$ is differentiable in its arguments.

$\begin{array}{l}\therefore \delta I={\displaystyle \sum }_e\delta I_e=2{\displaystyle \sum }_e{\displaystyle \sum }_{i=1}^3{\left(E_i^e,\delta E_i^e\right)}_{{\overline{\Omega }}_{xt}^e}\\ =2{\displaystyle \sum }_e{\displaystyle \sum }_{i=1}^3{\left\{g_e\right\}}_i=2{\displaystyle \sum }_e{\left\{g\right\}}_e=\left\{g\right\}=0\end{array}$ (34)

Let

$\left\{\delta \right\}={\left\{{\delta }^e\right\}}_u\cup {\left\{{\delta }^e\right\}}_{\sigma }\cup {\left\{{\delta }^e\right\}}_v$ (35)

be the total degrees of freedom for ${\left({\overline{\Omega }}_{xt}^{\left(1\right)}\right)}^T$ . Then $\left\{g\right\}$ in (34) is a nonlinear function of $\left\{\delta \right\}$ as the PDEs in the mathematical models are nonlinear. We obtain a solution $\left\{\delta \right\}$ that satisfies (34) using Newton’s linear method with line search [27].

4. Model Problem and Their Solutions

First, we non-dimensionalize the PDEs in the mathematical model ((26)-(28)). We write the BLM, the constitutive theory and other equations with a hat (^{^}) on all quantities, indicating they all have their usual dimensions or units.

${\widehat{\rho }}_0\frac{\partial {\widehat{v}}_1}{\partial \widehat{t}}-{\widehat{\rho }}_0{\widehat{F}}_1^b+\left(-\widehat{C}\left[{\left(1+\frac{\partial {\widehat{u}}_1}{\partial {\widehat{x}}_1}\right)}^{-2}\frac{{\partial }^2{\widehat{u}}_1}{\partial {\widehat{x}}_1{}^2}\right]\right)-\frac{{\partial }^2{\widehat{u}}_1}{\partial {\widehat{x}}^2}{}_d\widehat{\sigma }{}_{11}^{\left[0\right]}-\left(1+\frac{\partial {\widehat{u}}_1}{\partial {\widehat{x}}_1}\right)\frac{{\partial }_d{\widehat{\sigma }}^{\left[0\right]}}{\partial {\widehat{x}}_1}=0$ (36)

${}_d{\widehat{\sigma }}_{11}^{\left[0\right]}={\widehat{C}}_1\left(\frac{\partial {\widehat{u}}_1}{\partial {\widehat{x}}_1}+\frac{1}{2}{\left(\frac{\partial {\widehat{u}}_1}{\partial {\widehat{x}}_1}\right)}^2\right)+{\widehat{C}}_2\left(\frac{\partial {\widehat{v}}_1}{\partial {\widehat{x}}_1}+\frac{\partial {\widehat{u}}_1}{\partial {\widehat{x}}_1}\frac{\partial {\widehat{v}}_1}{\partial {\widehat{x}}_1}\right)$ (37)

${\widehat{v}}_1=\frac{\partial {\widehat{u}}_1}{\partial \widehat{t}}$ (38)

We choose the following reference and dimensionless quantities:

$\begin{array}{l}\begin{array}{lll}{x}_1=\frac{{\widehat{x}}_1}{L_0},\hfill & {\rho }_0=\frac{{\widehat{\rho }}_0}{{\left({\rho }_0\right)}_{\text{ref}}},\hfill & u_1=\frac{{\widehat{u}}_1}{L_0},\hfill \\ {\sigma }_{11}^{\left[0\right]}=\frac{{\widehat{\sigma }}_{11}}{{\tau }_0},\hfill & t_0=\frac{L_0}{v_0},\hfill & {\tau }_0=p_0={\left({\rho }_0\right)}_{\text{ref}}{v}_0^2,\hfill \\ v_0=\sqrt{\frac{E_0}{{\left({\rho }_0\right)}_{\text{ref}}}},\hfill & E=\frac{\widehat{C}}{E_0},\hfill & \eta =\frac{{\widehat{C}}_2}{{\eta }_0},\hfill \\ F_1^b=\frac{{\widehat{F}}_1^b}{F_0},\hfill & F_0=\frac{v_0^2}{L_0},\hfill & C=\frac{\widehat{C}}{{\tau }_0}=\frac{\widehat{C}}{E_0}\hfill \end{array}\hfill \end{array}\}$ (39)

Using (39) in (36)-(38), we can obtain their following dimensionless forms:

${\rho }_0\frac{\partial v_1}{\partial t}-{\rho }_0{F}_1^b+\left(-C\left[{\left(1+\frac{\partial u_1}{\partial x_1}\right)}^{-2}\frac{{\partial }^2{u}_1}{\partial x_1{}^2}\right]\right)-\frac{{\partial }^2{u}_1}{\partial x^2}{}_d\sigma {}_{11}^{\left[0\right]}-\left(1+\frac{\partial u_1}{\partial x_1}\right)\frac{{\partial }_d{\sigma }^{\left[0\right]}}{\partial x_1}=0$ (40)

${}_d{\sigma }_{11}^{\left[0\right]}=E\left(\frac{\partial u_1}{\partial x_1}+\frac{1}{2}{\left(\frac{\partial u_1}{\partial x_1}\right)}^2\right)+\frac{\eta }{Re}\left(\frac{\partial v_1}{\partial x_1}+\frac{\partial u_1}{\partial x_1}\frac{\partial v_1}{\partial x_1}\right)$ (41)

$v_1=\frac{\partial u_1}{\partial t}\text{where}Re=\frac{{\rho }_{\text{ref}}{v}_0{L}_0}{{\eta }_0}$ (42)

$Re$ is Reynolds number and $\frac{\eta }{Re}$ is the dimensionless damping coefficient. $E$

is the dimensionless modulus of elasticity.

We consider a one-dimensional rod of constant cross-section of compressible TVE solid material with a dimensionless length $L=1$ , fixed at the left end ($x_1=0$ ), and subjected to a tensile stress or a velocity pulse of duration $2\text{Δ}t$ at the right end $x_1=L$ . Figure 2(a) shows the schematic and the loading at $x_1=L$ . If we assume that all points at a cross-section of the rod displace in the $x_1$ direction by the same amount, then the problem becomes one-dimensional and we can idealize it by a line from $x_1=0$ to $x_1=L$ (Figure 2(b)). Figure 2(c) shows the first space-time strip ${\overline{\Omega }}_{xt}^{\left(1\right)}=\left[0,L\right]\times \left[0,\Delta t\right]$ .

Figure 3. Discretization of space-time domain into space-time strips, discretization of $n^{th}$ space-time strip into space-time elements.

This strip is discretized (uniformly) using nine-node p-version hierarchical space-time finite elements with higher global differentiability. Since the initial value problem contains up to second-order derivatives in space and time, $k_1=k_2=3$ (solution of class $C^2$ in space and time) is the minimally conforming space for which the space-time integrals over the discretization ${\left({\overline{\Omega }}_{\left(xt\right)}^{\left(1\right)}\right)}^T$ are Riemann integrable. When $k_1=k_2=2$ (solutions of class $C^1$ in space and time), the space-time integrals over ${\left({\overline{\Omega }}_{\left(xt\right)}^{\left(1\right)}\right)}^T$ are Lebesgue integrable. The applied disturbance is a stress or a velocity pulse of duration $2\text{Δ}t$ at $x=1$ . If we

consider a velocity pulse of $v_1$ , then $v_1=0$ , $\frac{\partial v_1}{\partial t}=0$ at $t=0$ ; $v_1=v_1^{\text{*}}$ and $\frac{\partial v_1}{\partial x_1}=0$ at $t=\text{Δ}t$ ; $v_1=0$ and $\frac{\partial v_1}{\partial x_1}=0$ at $t=2\text{Δ}t$ ; and $v_1=0$ for $t\ge 2\text{Δ}t$ .

Similar conditions hold for the stress pulse ${}_d{\sigma }_{11}^{\left[0\right]}$ . This mathematical model ((40)-(42)) contains first order derivatives of $u_1$ and $v_1$ in time, hence the solutions of class $C^1$ in time is the minimum requirement for the space-time integral to be Riemann over the space-time discretization.

We choose the following material properties and reference quantities:

$L_0=1$ ${\left({\rho }_0\right)}_{\text{ref}}=1240$ $v_0=116.05$ ${\eta }_0=372$

$\widehat{\rho }=1240$ $\widehat{E}=16.7\times {10}^6$ $\widehat{\eta }=129.51$

$E_0={\left({\rho }_0\right)}_{\text{ref}}\times v_0^2=16.7\times {10}^6$

5. Model Problem Studies

In this section we consider the following studies:

1. Wave propagation in linear TVE solid

2. Shock fronts in compressible TE and TVE solids

3. Shock fronts for increasing $v^{\ast }$ with fixed damping

4. Shock fronts for fixed $v^{\ast }$ with variable damping

5. Shock fronts formation, reflection and transmission.

Details of each of these investigations are presented in the following. In all studies, we use a uniform 30 nine-node p-version space-time element discretization, $\text{Δ}t=0.1$ , and p-levels of seven in space and time (unless specified otherwise) for a space-time strip. Convergence studies conducted confirm that increasing p-levels beyond seven does not result in measurable improvements in the calculated solution. Newton’s linear method for solving nonlinear algebraic equations is considered converged when $\text{max}\left|g_i\right|\le \text{Δ}$ , with $\text{Δ}=O\left({10}^{-6}\right)$ . Computed evolution for a space-time strip is considered converged to the true solution of the IVP when $I\le O\left({10}^{-8}\right)$ or lower for the discretization ${\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T$ for the i-th space-time strip. All computed solutions reported in the paper satisfy these criteria. In all studies we use $\text{Δ}t=0.1$ unless specified otherwise. In the graph presented for model problems, the pair of numbers in the upper left or right corner represent the abscissa (x-coordinate) and ordinate (value of the quantity) represents the peak values. Color legends correspond to the graphs. This information is included for the quantitative comparison of the peak values and their locations.

5.1. Wave Propagation in Linear TVE Solid

We apply a velocity pulse of peak value $v^{\ast }=0.1$ and its support of $2\text{Δ}t$ at $x_1=1.0$ . In linear TVE solids, due to small strain and small deformation physics, the material is nearly incompressible, hence density remains constant during the evolution. If the stiffness does not change during deformation, which is the case in linear viscoelasticity, then the wave speed remains the same for TE solids as well as for TVE solids, regardless of the extent of dissipation. However, the stresses will differ in TE solids and TVE solids. We expect TVE solids to have lower stress compared to TE solids due to viscous dissipation converting some of the mechanical work into entropy.

Figures 4(a)-(d) show the plots of ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=2\text{Δ}t$ , $t=5\text{Δ}t$ , $t=10\text{Δ}t$ and $t=16\text{Δ}t$ . In case of TE solids, the wave remains unchanged during evolution as there is no mechanism of energy conversion. We clearly observe lower peak values of ${}_d{\sigma }_{11}^{\left[0\right]}$ for TVE solids compared to TE solids due to dissipation. We also note amplitude decay and base elongation in case of TVE solids in Figure 4(a) and Figure 4(d). Small changes in the amplitude and base for TE solids can be further reduced by increasing $p-$ levels.

Even though for infinitesimal deformation, the material is assumed incompressible, however we know that since $\frac{\partial u_1}{\partial x_1}\ne 0$ , which must be the case to have strain and stress. Thus, a change in density due to $\frac{\partial u_1}{\partial x_1}$ is inevitable.

$\rho \left(x,t\right)=\frac{{\rho }_0\left(x,t\right)}{\left(1+\frac{\partial u_1}{\partial x_1}\right)}$

Figure 5(a) and Figure 5(d) show plots of $\rho$ versus $x_1$ at $t=2\text{Δ}t$ , $t=5\text{Δ}t$ , $t=10\text{Δ}t$ and $t=16\text{Δ}t$ . Density values follow stress ${}_d{\sigma }_{11}^{\left[0\right]}$ , i.e., an increase in ${}_d{\sigma }_{11}^{\left[0\right]}$ results in a decrease in $\rho$ . By reducing $v^{\text{*}}$ , we can obviously reduce the change in density, suggesting that we are approaching the incompressible case. We clearly observe the same wave speed for TE and TVE solids.

5.2. Shock Fronts in Compressible TE and TVE Solids

We apply a velocity pulse of maximum amplitude $v^{\text{*}}=0.1$ and of support $2\text{Δ}t$ at the right end of the rod, $x_1=1.0$ . Figures 6(a)-(d) plots of ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1=1.0$ for TE and for TVE solids with damping coefficient $C_2=0.0009$ at $t=2\text{Δ}t$ , $5\text{Δ}t$ , $10\text{Δ}t$ and $15\text{Δ}t$ . At $t=2\text{Δ}t$ , the stress pulse has just entered the rod (Figure 6(a)). At $t=5\text{Δ}t$ , we observe formation of shock front ahead of the wave, lower peak and increased base of the wave for damped case is seen clearly. Mild oscillations appear next to the shock front for the undamped case as in this case there is no mechanism to dampen the vibrations except lack of proximity to the shock front. For damped case, the wave is relatively oscillation-free as the vibration of the material points are damped. Figure 6(c) shows partial reflection process of the two fronts from the impenetrable boundary at $x_1=0.0$ . At $t=15\text{Δ}t$ , the reflected fronts are moving towards $x_1=1.0$ . First, we observe that the shock front in the reflected wave is ahead of the wave as it was before reflection. We observe oscillation in the undamped case, both behind and ahead of the shock, but more pronounced ahead of the shock as they should be. These are undamped vibrations of the material points. When damping is present (even though very small, $C_2=0.0009$ ) the wave is oscillation-free, peak is lower, base slightly increased and the shock front strength (Slope) is slightly reduced, as expected. Figures 7(a)-(d) show plots of $\rho$ versus $x_1$ for $t=2\text{Δ}t$ , $5\text{Δ}t$ , $10\text{Δ}t$ and $15\text{Δ}t$ . We observe behavior parallel to ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ . In the absence of damping, pre- and post-shock oscillations are present, post shock oscillation being more dominant as expected, lower density peak and larger support of the density wave are also observed. In the presence of damping, evolution of $\rho$ is smooth for all values of time, shock fronts ahead of the wave are observed in the incident as well as reflected waves. Amplitude decay and base elongation of the waves are observed in the presence of damping.

5.3. Shock Fronts for Increasing $v^{\ast }$ with Fixed Damping

In this study, we consider $v^{\text{*}}=0.1$ and 0.2 for $C_2=0.001$ , i.e., we investigate the influence of higher magnitude of the velocity pulse for fixed damping. Figures 8(a)-(d) show plots of ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ , and Figures 9(a)-(d) show plots of $\rho$ versus $x_1$ for $t=2\text{Δ}t$ , $5\text{Δ}t$ , $10\text{Δ}t$ , and $15\text{Δ}t$ . First, due to the presence of damping all evolutions are completely oscillation free. Formation of shock fronts for ${}_d{\sigma }_{11}^{\left[0\right]}$ and $\rho$ are observed for both values of $v^{\text{*}}$ . The pulse with $v^{\text{*}}=0.2$ has much higher energy content, hence results in larger peaks of stress and density and higher base elongation with much steeper shock front. We clearly observe the shock fronts of ${}_d{\sigma }_{11}^{\left[0\right]}$ and $\rho$ are always ahead of the incident and the reflected waves. There is no change in wave speed due to increase in $v^{\text{*}}$ , as expected.

5.4. Shock Fronts for Fixed $v^{\ast }$ with Variable Damping

In this study, we consider fixed $v^{\text{*}}$ , $v^{\text{*}}=0.2$ and choose two values of damping, $C_2=0.001$ and $C_2=0.0015$ . Fixed peak velocity pulse of duration $2\text{Δ}t$ imparts fixed energy to the rod for each value of damping. Figures 10(a)-(d) show plots of ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ for $t=2\text{Δ}t$ , $5\text{Δ}t$ , $10\text{Δ}t$ , and $15\text{Δ}t$ . Due to the presence of damping, all evolutions are oscillation free. The shock fronts of ${}_d{\sigma }_{11}^{\left[0\right]}$ and $\rho$ are observed for both values of damping and always appear ahead of the incident and reflected waves. Larger value of damping results in lower peaks, elongated base and evolution with less or no oscillations. There is no change in wave speed due to changes in damping, as expected. Figure 11(a)-(d) show graphs of density $\rho$ versus $x_1$ for the same values of time as for ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ . Density $\rho$ versus $x_1$ show similar behavior as ${}_d{\sigma }_{11}^{\left[0\right]}$ vs $x_1$ . Since the two values of the damping coefficient are relatively close, significant base elongation is not observed for $C_2=0.0015$ ; however, amplitude decay is quite noticeable compared to $C_2=0.001$ .

5.5. Formation of Shock Fronts and Their Reflection and Transmission

In this study we consider formation, propagation, interaction and transmission of waves with shock fronts. We choose TVE solids with $C_2=0.0009$ . A velocity pulse of duration $2\text{Δ}t$ with $v^{\text{*}}=0.1$ is applied at both ends of a TVE rod of dimensionless length of two units. The space-time strip is discretized using a 60 element uniform mesh of p-version hierarchical space-time elements with p = 7 in space and time and $k_1=k_2=2$ i.e., solutions of class $C_1$ in space and time. Figure 12(a) shows the evolution at $t=2\text{Δ}t$ ; both velocity pulses are in the TVE rod. At $t=8\text{Δ}t$ and $t=9\text{Δ}t$ (Figure 12(b) and Figure 12(c) respectively), the waves are at the center of the rod (as expected due to wave speed = 1) with shock fronts behind them. At $t=14\text{Δ}t$ (Figure 12(c)), the waves, after interaction, have resumed their original identities but with some base elongation and amplitude decay due to dissipation. Their interaction at $x=1$ is a spatial local change in the wave at $x=1$ for the particular value of time, but upon further evolution, the waves resume their original identities as expected and shown in Figure 12(d).

6. Summary and Conclusion

This paper presents investigations of various aspects of tensile shock physics in compressible TVE solid matter without rheology using mathematical models for finite deformation, finite strain deformation physics based on the CBL of CCM derived using the contravariant second Piola-Kirchhoff stress tensor, Green’s strain tensor and its rates up to order n. The solutions of the resulting IVPs are obtained using a space-time coupled variationally consistent, unconditionally stable, space-time finite element method based on space-time residual functional. To the authors’ knowledge, this work has not been reported in the published literature. We summarize the work presented in this paper and draw some conclusions:

1. Existence of tensile shock physics requires elastic medium; thus, this physics is possible in compressible TE and TVE solids with and without memory. Basic mechanism of tensile shock physics is due to the change in density along the base of the tensile wave. For an applied tensile velocity pulse, decrease in density behind the pulse and progressively increasing density ahead of the pulse result in piling up of the waves ahead of the pulse, hence formation of a tensile shock ahead of the pulse. Due to lack of elasticity, this physics cannot exist in compressible fluids.

2. Mathematical model is based on CBL of CCM and the constitutive theories are derived using conjugate pairs in entropy inequality and the representation theorem; hence, the mathematical model is thermodynamically and mathematically consistent.

3. Solutions of the equations in the mathematical model describing IVPs are obtained using space-time coupled finite element formulation in higher order spaces (hpk) based on space-time residual functional; hence, the resulting computations are unconditionally stable and permit accurate a posteriori error computation.

4. In tensile shock physics, shocks always form ahead of the wave as opposed to compression shock physics in which the shocks always form behind the wave.

5. We have shown that high-energy impact results in pre- and post-shock oscillations in the evolution of ${\sigma }_{11}^{\left[0\right]}$ and $\rho$ in the undamped compressible medium. Due to the lack of dissipation, the vibrations of material points persist. These oscillations are naturally more pronounced in the vicinity of the shock, where there is higher energy for vibrations. Even the slightest amount of damping effectively suppresses the vibrations of the material points, thereby eliminating the oscillations.

6. The presence of damping whether low or high and the strength of the applied disturbance (high $v^{\text{*}}$ ) do not affect the wave speed. Higher damping reduces the peak amplitudes, results in more elongation of the base of the wave, and results in a shallower wave front. Conversely, an increase in $v^{\text{*}}$ raises the peak amplitudes, while the base remains unaffected (as long as the damping remains unchanged). However, the shock fronts become steeper due to the higher energy content of the applied pulse.

In summary, the work presented in this paper addresses all aspects of tensile shock physics of stress and density waves in compressible TE and TVE solids without memory: shock front formation, propagation, reflection, interaction, influence of dissipation and influence of the strength of external stimulus are all addressed in the paper. The most significant aspects of the work are: (i) mathematical models are precisely based on CBL of CCM and the constitutive theories are derived using representation theorem; hence, the complete mathematical models are thermodynamically and mathematically consistent. (ii) The space-time coupled finite element method based on the space-time residual functional in higher-order scalar product spaces with minimally conforming spaces ensures that all space-time integrals over discretized space-time domain are Riemann. This permits rapid convergence of the computed solutions to the true solution and accurate computations of a posteriori errors based on space-time residual functional. The solutions reported in the paper are nearly as accurate as theoretical solutions, with the space-time residual functional for each space-time discretization consistently achieving $O\left({10}^{-8}\right)$ or lower. This level of accuracy is ensured by employing minimally conforming spaces and higher p-levels. Control over the p-level and the order of the approximation space allows for a more accurate description of the physics in the computational process. (iii) Complex process of tensile shock formation, propagation, reflection, interaction and transmission has been simulated with higher degree of accuracy using (2) and (3).

(a) (b)

(c) (d)

Figure 4. (a) ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=2\text{Δ}t$ ; (b): ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=5\text{Δ}t$ ; (c): ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=10\text{Δ}t$ ; (d): ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=16\text{Δ}t$ - LinearCase - Thermoelastic and Thermoviscoelastic with $C_2=0.0009$ .

(a) (b)

(c) (d)

Figure 5. (a): $\rho$ versus $x_1$ at $t=2\text{Δ}t$ ; (b): $\rho$ versus $x_1$ at $t=5\text{Δ}t$ ; (c): $\rho$ versus $x_1$ at $t=10\text{Δ}t$ ; (d): $\rho$ versus $x_1$ at $t=16\text{Δ}t$ - LinearCase - Thermoelastic and Thermoviscoelastic with $C_2=0.0009$ .

(a) (b)

(c) (d)

Figure 6. (a): ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=2\text{Δ}t$ ; (b): ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=5\text{Δ}t$ ; (c): ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=10\text{Δ}t$ ; (d): ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=15\text{Δ}t$ - NonLinearCase - Thermoelastic and Thermoviscoelastic with $C_2=0.0009$ .

(a) (b)

(c) (d)

Figure 7. (a): $\rho$ versus $x_1$ at $t=2\text{Δ}t$ - $NonLinearCase$ ; (b): $\rho$ versus $x_1$ at $t=5\text{Δ}t$ ; (c): $\rho$ versus $x_1$ at $t=10\text{Δ}t$ ; (d): $\rho$ versus $x_1$ at $t=15\text{Δ}t$ - NonLinearCase - Thermoelastic and Thermoviscoelastic with $C_2=0.0009$ .

(a) (b)

(c) (d)

Figure 8. (a): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus x1 at $t=1\text{Δ}t$ ; (b): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus x1 at $t=7\text{Δ}t$ ; (c): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus x1 at $t=10\text{Δ}t$ ; (d): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus x1 at $t=15\text{Δ}t$ (influence of increasing $v^{\text{*}}$ ).

(a) (b)

(c) (d)

Figure 9. (a): $\rho$ versus $x_1$ at $t=1\text{Δ}t$ ; (b): $\rho$ versus $x_1$ at $t=7\text{Δ}t$ ; (c): $\rho$ versus $x_1$ at $t=10\text{Δ}t$ ; (d): $\rho$ versus $x_1$ at $t=15\text{Δ}t$ (influence of increasing $v^{\text{*}}$ ).

(a) (b)

(c) (d)

Figure 10. (a): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus $x_1$ at $t=2\text{Δ}t$ ; (b): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus $x_1$ at $t=7\text{Δ}t$ ; (c): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus $x_1$ at $t=10\text{Δ}t$ ; (d): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus $x_1$ at $t=15\text{Δ}t$ (influence of damping).

(a) (b)

(c) (d)

Figure 11. (a): $\rho$ versus $x_1$ at $t=2\text{Δ}t$ ; (b): $\rho$ versus $x_1$ at $t=7\text{Δ}t$ ; (c): $\rho$ versus $x_1$ at $t=10\text{Δ}t$ ; (d): $\rho$ versus $x_1$ at $t=15\text{Δ}t$ (influence of damping).

(a) (b)

(c) (d)

Figure 12. (a): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus $x_1$ at $t=2\text{Δ}t$ ; (b): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus $x_1$ at $t=8\text{Δ}t$ ; (c): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus $x_1$ at $t=9\text{Δ}t$ ; (d): ${}^d{\sigma }_{11}{}^{\left[0\right]}$ versus $x_1$ at $t=14\text{Δ}t$ (propagation, interaction of shock waves).

List of Symbols

$C$ : Bulk Modulus

$C_1$ : Elastic Coefficient

$C_2$ : Damping Coefficient

$\overline{x}$ , ${\overline{x}}_i$ , $\left\{\overline{x}\right\}$ : deformed Coordinates

$x$ , $x_i$ , $\left\{x\right\}$ : undeformed Coordinates

$\left[J\right]$ : Deformation gradient tensor

$\left|J\right|$ : Determinant of $\left[J\right]$

$u$ , $u_i$ , $\left\{u\right\}$ : displacements in Lagrangian description

$v$ , $v_i$ , $\left\{v\right\}$ : velocities in Lagrangian description

${\rho }_0$ : density at time $t_0$

${\left({\rho }_0\right)}_{\text{ref}}$ : reference density

$\rho$ : density in Lagrangian description in the current configuration

$\eta$ : specific entropy in Lagrangian description

${\eta }_0$ : Reference Dissipation

$\lambda$ : Relaxation time

$t_0$ : Reference time

$e$ : specific internal energy in Lagrangian description

$p$ : thermodynamic or Mechanical Pressure in Lagrangian description

$\theta$ : temperature in Lagrangian description

$E$ : Elastic modulus

$q$ , $q_i$ , $\left\{q\right\}$ : heat vector in Lagrangian description

$g$ , $g_i$ , $\left\{g\right\}$ : temperature gradient tensor in Lagrangian description

${\sigma }^{\left[0\right]}$ , $\left\{{\sigma }^{\left[0\right]}\right\}$ , ${\sigma }_{11}^{\left[0\right]}$ : Second Piola-Kirchhoff stress

${}_e{\sigma }^{\left[0\right]}$ , $\left\{{}_e{\sigma }^{\left[0\right]}\right\}$ : Equilibrium part of the Second Piola-Kirchhoff stress

${}_d{\sigma }^{\left[0\right]}$ , $\left\{{}_d{\sigma }^{\left[0\right]}\right\}$ : Deviatoric part of the Second Piola-Kirchhoff stress

${}_d{\sigma }^{\left[j\right]}$ , $\left\{{}_d{\sigma }^{\left[j\right]}\right\}$ : Convected time derivative of order j of the Deviatoric Second Piola-Kirchhoff stress tensor

${\sigma }^{\left(0\right)}$ , $\left\{{\sigma }^{\left(0\right)}\right\}$ , ${\sigma }_{11}^{\left(0\right)}$ : Cauchy stress

${}_e{\sigma }^{\left(0\right)}$ , $\left\{{}_e{\sigma }^{\left(0\right)}\right\}$ : Equilibrium part of the Cauchy stress

${}_d{\sigma }^{\left(0\right)}$ , $\left\{{}_d{\sigma }^{\left(0\right)}\right\}$ : Deviatoric part of the Cauchy stress

${\epsilon }_{\left[0\right]}$ , $\left\{{\epsilon }_{\left[0\right]}\right\}$ : Green’s strain tensor (convected time derivative of order zero)

${\epsilon }_{\left[i\right]}$ , $\left\{{\epsilon }_{\left[i\right]}\right\}$ : Convected time derivative of order i of the covariant Green’s strain tensor

Abbreviations

CBL: Conservation and Balance Laws

CCM: Classical Continuum Mechanics

TE: Thermoelastic Solid

TVE: Thermoviscoelastic Solid

BVPs: Boundary Value Problems

IVPs: Initial Value Problems

CM: Conservation of Mass

BLM: Balance of Linear Momenta

BAM: Balance of Angular Momenta

Acknowledgements

The first author is grateful for his endowed professorships and the department of mechanical engineering of the University of Kansas for providing financial support to the second author. The computational facilities provided by the Computational Mechanics Laboratory of the mechanical engineering department are also acknowledged.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Surana, K.S. and Abboud, E. (2024) Shock Physics in Compressible Thermoelastic and Thermoviscoelastic Solids. Meccanica.
[2]	Surana, K.S. and Abboud, E. (2024) Shock Physics in Compressible Thermoviscoelastic Solid Medium with Rheology. Acta Meccanica.
[3]	Surana, K.S., Kendall, J.K. and Carranga, C.H. (2022) Formation, Propagation and Reflection of 1D Normal Shocks in Riemann Shock Tube. Applied Mathematics, 13, 295-323. https://doi.org/10.4236/am.2022.133022
[4]	Surana, K.S., Reddy, K.P.J., Joy, A.D. and Reddy, J.N. (2014) Riemann Shock Tube: 1D Normal Shocks in Air, Simulations and Experiments. International Journal of Computational Fluid Dynamics, 28, 251-271. https://doi.org/10.1080/10618562.2014.927056
[5]	Dobróka, M. and Molnár, J.S. (2014) An Introduction to Continuum Mechanics and Elastic Wave Propagation. https://geofizika.uni-miskolc.hu/Engineering%20physics%20part%20I..pdf
[6]	Bechir, H. and Benslimane, A. (2017) On the Propagation of Weak Shock Waves in Compressible Thermohyperelastic Solids. Acta Mechanica, 229, 87-97. https://doi.org/10.1007/s00707-017-1961-x
[7]	Shams, M. and Ejaz, K. (2022) Interfacial Wave Propagation in Initially Stressed Compressible Hyperelastic Materials. Mathematics and Mechanics of Solids, 27, 2510-2531. https://doi.org/10.1177/10812865221074304
[8]	Zaid, A.I.O. (2016) Stress Waves in Solids, Transmission, Reflection and Interaction and Fractures Caused by Them: State of the Art. International Journal of Theoretical and Applied Mechanics, 1, 155-164.
[9]	Wang, Y., Wang, Y. and Laude, V. (2015) Wave Propagation in Two-Dimensional Viscoelastic Metamaterials. Physical Review B, 92, Article 104110. https://doi.org/10.1103/physrevb.92.104110
[10]	Gulizzi, V. and Saye, R. (2022) Modeling Wave Propagation in Elastic Solids via High-Order Accurate Implicit-Mesh Discontinuous Galerkin Methods. Computer Methods in Applied Mechanics and Engineering, 395, Article 114971. https://doi.org/10.1016/j.cma.2022.114971
[11]	Zhang, G.M. and Batra, R.C. (2007) Wave Propagation in Functionally Graded Materials by Modified Smoothed Particle Hydrodynamics (MSPH) Method. Journal of Computational Physics, 222, 374-390. https://doi.org/10.1016/j.jcp.2006.07.028
[12]	Surana, K.S., Knight, J. and Reddy, J.N. (2015) Nonlinear Waves in Solid Continua with Finite Deformation. American Journal of Computational Mathematics, 5, 345-386. https://doi.org/10.4236/ajcm.2015.53032
[13]	Surana, K.S. (2015) Advanced Mechanics of Continua. CRC/Taylor and Francis.
[14]	Surana, K.S. (2022) Classical Continum Mechanics. 2nd Edition, CRC/Taylor and Francis.
[15]	Smith, G.F. (1965) On Isotropic Integrity Bases. Archive for Rational Mechanics and Analysis, 18, 282-292. https://doi.org/10.1007/bf00251667
[16]	Smith, G.F. (1970) On a Fundamental Error in Two Papers of C.-C. Wang “on Representations for Isotropic Functions, Parts I and II”. Archive for Rational Mechanics and Analysis, 36, 161-165. https://doi.org/10.1007/bf00272240
[17]	Smith, G.F. (1971) On Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 9, 899-916. https://doi.org/10.1016/0020-7225(71)90023-1
[18]	Spencer, A.J.M. (1971) Theory of Invariants. In: Mathematics, Elsevier, 239-353. https://doi.org/10.1016/b978-0-12-240801-4.50008-x
[19]	Spencer, A.J.M. and Rivlin, R.S. (1958) The Theory of Matrix Polynomials and Its Application to the Mechanics of Isotropic Continua. Archive for Rational Mechanics and Analysis, 2, 309-336. https://doi.org/10.1007/bf00277933
[20]	Spencer, A.J.M. and Rivlin, R.S. (1959) Further Results in the Theory of Matrix Polynomials. Archive for Rational Mechanics and Analysis, 4, 214-230. https://doi.org/10.1007/bf00281388
[21]	Wang, C.C. (1969) On Representations for Isotropic Functions. Archive for Rational Mechanics and Analysis, 33, 249-267. https://doi.org/10.1007/bf00281278
[22]	Wang, C.C. (1969) On Representations for Isotropic Functions. Archive for Rational Mechanics and Analysis, 33, 268-287. https://doi.org/10.1007/bf00281279
[23]	Wang, C.-C. (1970) A New Representation Theorem for Isotropic Functions: An Answer to Professor G. F. Smith's Criticism of My Papers on Representations for Isotropic Functions. Archive for Rational Mechanics and Analysis, 36, 166-197. https://doi.org/10.1007/bf00272241
[24]	Wang, C.-C. (1971) Corrigendum to My Recent Papers on “Representations for Isotropic Functions”. Archive for Rational Mechanics and Analysis, 43, 392-395. https://doi.org/10.1007/bf00252004
[25]	Zheng, Q.-S. (1993) On the Representations for Isotropic Vector-Valued, Symmetric Tensor-Valued and Skew-Symmetric Tensor-Valued Functions. International Journal of Engineering Science, 31, 1013-1024. https://doi.org/10.1016/0020-7225(93)90109-8
[26]	Zheng, Q.-S. (1993) On Transversely Isotropic, Orthotropic and Relative Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. Part I: Two Dimensional Orthotropic and Relative Isotropic Functions and Three Dimensional Relative Isotropic Functions. International Journal of Engineering Science, 31, 1399-1409. https://doi.org/10.1016/0020-7225(93)90005-f
[27]	Surana, K.S. and Reddy, J.N. (2018) The Finite Element Method for Initial Value Problems. CRC/Taylor and Francis.
[28]	Bell, B.C. and Surana, K.S. (1994) A Space-Time Coupled P-Version Least-Squares Finite Element Formulation for Unsteady Fluid Dynamics Problems. International Journal for Numerical Methods in Engineering, 37, 3545-3569. https://doi.org/10.1002/nme.1620372008
[29]	Surana, K.S., Allu, S., Reddy, J.N. and Tenpas, P.W. (2008) Least-Squares Finite Element Processes in h, P, K Mathematical and Computational Framework for a Non-Linear Conservation Law. International Journal for Numerical Methods in Fluids, 57, 1545-1568. https://doi.org/10.1002/fld.1695
[30]	Surana, K.S., Allu, S., Romkes, A. and Reddy, J.N. (2009) Evolution, Propagation, Reflection, and Interactions of 1d-Normal Shocks in Air and FC70 Using Hpkfinite Element Computational Framework. International Journal for Computational Methods in Engineering Science and Mechanics, 10, 370-392. https://doi.org/10.1080/15502280903106523