Karan S. Surana^*, Stacia D. Zahn
Department of Mechanical Engineering, University of Kansas, Lawrence, Kansas, USA.
DOI: 10.4236/jamp.2026.148140   PDF    HTML   XML   3 Downloads   42 Views   Citations

Abstract

This paper considers shock physics, entropy generation and associated thermal physics in compressible elastic, and thermoviscoelastic solid medium with and without rheology. The mathematical model containing nonlinear partial differential equations consists of conservation of mass and balance of linear momenta and energy equation derived using contravariant second Piola-Kirchhoff stress tensor and covariant Green’s strain tensor. This mathematical model is augmented with constitutive theories for contravariant second Piola-Kirchhoff stress tensor, heat vector and equation of state, thermodynamic pressure. The dissipation and rheology mechanisms are incorporated using rates of Green’s strain tensor up to orders n and rates of second Piola-Kirchhoff stress tensor up to orders m. Hence, the mathematical model contains spectra of dissipation coefficients and relaxation times. The solution of the mathematical model is obtained using space-time coupled finite element method based on space-time residual functional in which space-time local approximations are p-version hierarchical in higher order scalar product spaces, thus permitting desired orders of global differentiability of the approximations in space and time. The evolution is computed using a space-time strip for an increment of time followed by time marching. The space-time integral form in this approach is space-time variationally consistent, hence unconditionally stable computations are ensured during the entire evolution. Model problem studies are presented for one-dimensional wave propagation. Besides unconditional stability, other meritorious features of this computational methodology used here are discussed in the paper. Monitoring complex entropy generation and associated thermal physics in the presence of shock waves is a significant feature of the work presented here.

Keywords

Elastic, Viscoelastic, Dissipation, Rheology, Entropy Generation, Shock Physics, Space-Time Finite Elements, Variationally Consistent, Wave Physics, Classical Continuum Mechanics, Conservation and Balance Laws

Share and Cite:

1. Introduction and Literature Review

The stress wave physics in incompressible, linear elastic, isotropic and homogeneous solids is perhaps the simplest possible physics. In these solids, the incident stress wave maintains its amplitude and base (support) during propagation, thus preserving its mechanical energy. In incompressible, linear elastic, isotropic and homogeneous solids with description, i.e. in linear viscoelastic solids, the fixed energy content (mechanical energy) of a wave such as a stress or velocity pulse continuously results in entropy generation during propagation causing amplitude decay and base elongation of the wave or the pulse, but maintaining the wave shape and wave speed due to incompressibility of the medium. In such physics upon continued evolution, the entire energy of the pulse is converted into entropy resulting in temperature rise of the medium, leading to a constant temperature equilibrium as stationary state of the evolution of wave propagation.

In compressible, isotropic, homogeneous, inviscid solid medium, the compressibility causes change in density, thus for fixed modulus of elasticity, the wave speed $v=\sqrt{E/\rho }$ changes in the medium as the density changes, compression resulting in higher density reduces wave speed while tension reducing density causes increase in the wave speed. We consider a simple illustrative example to demonstrate the consequences of changing wave speed. Consider an axial rod consisting of isotropic, homogeneous elastic matter fixed at the left and subjected to a compressive stress pulse of duration $2\text{Δ}t$ at the right end. Let the dimensionless wave speed be one. At the end of $t=2\text{Δ}t$ , the wave is completely in the solid rod (Figure 1(a)). As the compressive stress increases from A to B (Figure 1(a)), the density increases from ${\rho }_0$ at point A (initial density) to maximum value $\rho >{\rho }_0$ at point B resulting in progressively reducing wave speed from A to B. Along B to C, the progressively reducing stress causes progressively reduced value of density from its peak value $\rho$ at point B and eventually at point C, we have density ${\rho }_0$ (Figure 1(b)). Thus, from A to B (or a to b), we have progressively reducing wave speed, the minimum value being at point B (or b) and from B to C (or b to c), we have progressively increasing wave speed and at C, we have reference wave speed (same as at A). Thus in the support AB of the pulse, the stress wave at $A_2$ is moving faster than the wave at B and the wave at $A_1$ is moving faster than the wave at $A_2$ resulting in piling up of the stress waves from A to B upon continued evolution causing steepening of wavefront AB as shown in Figure 1(a) at a later value of time ($t_n$ ). This is often referred to as shock front of the stress wave. On the other hand, from BC, the progressively increasing wave speed results in the waves ahead of a wave moving at a faster speed, i.e. wave at $B_1$ , is moving faster than the wave at B and the wave at $B_2$ is moving faster than the

Figure 1. Evolution of ${}_d{\sigma }_{11}^{\left[0\right]}$ and $\rho$ for compressive ${}_d{\sigma }_{11}^{\left[0\right]}$ . (a) Evolution of compressive ${}_d{\sigma }_{11}^{\left[0\right]}$ ; (b) Evolution of $\rho$ for compressive ${}_d{\sigma }_{11}^{\left[0\right]}$ .

wave at $B_1$ resulting in stretching or shallowing of the wave front BC shown in Figure 1(a) at a time $t_n>t_1$ . Thus, we see that in case of a compressive stress wave propagating in compressive solid medium the shock formation in the wave occurs behind the peak of the wave and the swallowing or rarefaction occurs ahead of the wave. Figure 1(b) shows initial density wave at time $t_1=2\text{Δ}t$ and at time $t_n>t_1$ the density wave with shock front and rarefaction occurring in the same manner as in case of stress wave.

Wave propagation studies in linear elastic solids have been a subject of investigation for a long time. Wave physics in compressible solids with dissipation, rheology and isothermal as well as nonisothermal physics have not been investigated much, hence published works in this area are very limited. Recently, Surana et al. [1]-[3] investigated shock physics in compressible solids with dissipation and rheology using simple equation of state under isothermal assumption, i.e. in these studies conversion of mechanical energy into entropy and associated thermal effects existed but were not monitored by not considering energy equation in the mathematical model. In these studies, the wave propagation domain was considered insulated, the entropy generation and thermal effects are only due to dissipation that are small, thus had virtually no effect in the material coefficients, i.e. deformation field and thermal physics were decoupled in these studies. Authors in references [1]-[3] have reported pertinent literature review related to wave propagation, mostly linear wave propagation with some work related to weak shock waves. The mathematical models in almost all cases are compromised descriptions compared to the conservation and balance laws of classical continuum mechanics and constitutive theories are almost never based on representation theorem [4]-[15]. The literature review presented by Surana et al. [1]-[3] and Abboud [16] is quite comprehensive but not repeated here for the sake of brevity. Interested readers can see references [1]-[3] [16].

The mathematical model used in the present work includes energy equation as part of the group of PDEs constituting the complete mathematical model. This aspect of the mathematical model is not considered in references [1]-[3]. The energy equation is essential to determine the conversion of mechanical work (some) into heat energy through entropy generation. Tracking of entropy generation through various phases of shock physics is helpful and essential in understanding the extend to which mechanical work is converted into thermal energy especially during shock formation and shock reflections. Higher entropy values indicate more mechanical energy being converted into thermal energy. The solutions of the mathematical model consisting of nonlinear PDEs in space and time are obtained by using space-time residual functional based space-time finite element formulation established and used extensively by Surana et al. [17] in which the space-time local approximation in hpk scalar product spaces permits higher order global differentiability in space and time, a novel feature of the computational framework.

2. Present Study

The research presented here is extension of the works of Surana et al. [1]-[3] for nonisothermal shock physics. In the work presented here, the energy equation is integral part of the mathematical model, hence entropy generation due to mechanical work and due to thermal physics is monitored and the resulting thermal physics associated with the shock waves is considered in addition to mechanical deformation. The three distinct aspects of this work presented here are described in the following.

The first aspect is the derivation of the details of the mathematical model consisting of: conservation and balance laws in Lagrangian description based on classical continuum mechanics: conservation of mass, balance of linear momenta, balance of angular momenta, energy equation, entropy inequality, constitutive theories and equation of state. In the derivation of the conservation and balance laws compressibility physics requires consideration of finite deformation, finite strain physics, hence we consider contravariant second Piola-Kirchhoff stress tensor and covariant Green’s strain tensor as measures of stress and strain. Dissipation physics is described by rates of Green’s strain tensor up to order n, hence ordered rate dissipation mechanism, yielding a spectrum of dissipation coefficients corresponding to the strain rates. Rheology or memory mechanism is incorporated using rates of contravariant second Piola-Kirchhoff stress tensor up to orders m, hence yielding ordered rate rheology mechanism with a spectrum of relaxation times. Energy equation accounts for the conversion of mechanical energy into entropy as well as other aspects of entropy generation. Entropy inequality cast in Helmholtz free energy density provides mechanism for determining constitutive tensors and the initial determination of their argument tensors. These argument tensors are augmented due to dissipation physics and rheology and some constitutive tensors are redefined due to consideration of rheology. Constitutive theories are derived using representation theorem [4]-[15] [18] [19]. Simple equation of state is considered in which the thermodynamic pressure only depends upon density. Equation of state is used to define equilibrium contravariant second Piola-Kirchhoff stress tensor. This mathematical model consists of nonlinear partial differential equations in dependent variable, space coordinates and time, hence constitutes IVP. The mathematical model in $R^3$ has closure. This mathematical model in $R^3$ is reduced to $R^1$ for pure one dimension wave physics studies in compressible nonlinear elastic medium with dissipation and rheology.

The second aspect of the work presented here is the unconditionally stable computational infrastructure used to obtain accurate solutions of IVPs that has built in measure of error in the computed solution (without the knowledge of theoretical solution) and built in adaptivity mechanism to improve the accuracy of the computed solutions. Based on Surana et al. [17], the space-time coupled method in which the integral forms are constructed using space-time residual functional and calculus of variation and the space-time local approximations are in higher-order scalar product spaces permitting higher degree of polynomials as well as higher order global differentiability is highly meritorious. In this approach, when the space-time integrals over the space-time discretization are Riemann and when the integrated sum of squares of space-time residual functional (I) approaches zero i.e. of the order of $\mathcal{O}\left({10}^{-8}\right)$ or lower, the PDEs in the mathematical model are satisfied accurately in the point-wise sense. Thus, proximity of residual functional (I) to zero is an absolute measure of error. This approach obviously does not require theoretical solutions. Adaptively to improve the solution to achieve lower values of I based on h,p,k is inherent in this computational method. Converged solution is obtained for the first space-time strip (or slab). This is followed by computation of the solution for the second space-time strip (or slab) using ICs from the first time strip. This is continued till the desired time is reached. This approach is efficient compared to space-time mesh. In this approach only after obtaining converged solution for the current space-time strip, the solution is advanced to the next space-time strip, hence ensuring converged solution for the all space-time strips, hence the entire space-time domain.

The third aspect of this work is model problem studies using 1D wave physics in compressible and incompressible solid medium with dissipation and rheology. In each study wave propagation, shock formation, propagation of waves with shocks, reflection of waves with shocks, interaction of waves, entropy generation during evolution resulting in propagating temperature waves with and without shocks, their reflections and interactions one considered.

The mathematical model in Lagrangian description consisting of conservation and balance laws of classical continuum mechanics, constitutive theories, and their dimensionless forms for nonlinear and nonisothermal 3D deformation physics of thermoelastic and thermoviscoelastic solid mediums with and without rheology is presented. This mathematical model is specialized for 1D wave propagation. Details of the space-time coupled finite element method based om space-time residual functional and calculus of variations and the solution procedure for nonlinear algebraic equations for a space-time strip with time marching are also presented in the paper. Extensive model problem studies are presented to clearly illustrate accurate simulation of complex shock physics. Summary and conclusion are given in the last section of the paper.

3. Mathematical Model

3.1. Conservation and Balance Laws of Classical Continuum Mechanics, Equation of State

The mathematical model is based on classical continuum thermodynamics for compressible solid matter in Lagrangian description. First, we consider the mathematical model in $R^3$ in which density change during deformation is necessary for shock formation, hence we consider finite deformation, finite strain deformation physics. Model problems use 1D form of this model in $R^3$ . This is not to be viewed as reduction of 3D model to 1D which implies that some physics that exists in $R^3$ has been compromised in deriving the model in $R^1$ . This is not the case in the present work. In this work, we consider constant elastic material properties, hence they are neither dependent on the invariants of the strain tensor nor on the temperature. We remark that shock physics i.e. formation of shock waves can only exist in compressible solid media. Basic mechanism of shock formation is the continuous changing density for fixed modulus of elasticity resulting continuously changing wave speed during evolution that causes piling up of waves which eventually forms a shock (explained in introduction).

We consider conservation and balance laws of classical continuum mechanics in Lagrangian description using contravariant second Piola-Kirchhoff stress tensor ${\sigma }^{\left[0\right]}$ and covariant Green’s strain tensor ${\epsilon }_{\left[0\right]}$ that are valid measures for nonlinear deformation (finite deformation finite strain). We consider additive decomposition of ${\sigma }^{\left[0\right]}$ , ${\sigma }^{\left[0\right]}={}_e{\sigma }^{\left[0\right]}+{}_d{\sigma }^{\left[0\right]}$ in which ${}_e{\sigma }^{\left[0\right]}$ is contravariant second Piola-Kirchhoff equilibrium stress tensor and ${}_d{\sigma }^{\left[0\right]}$ is contravariant second Piola-Kirchhoff deviatoric stress tensor. The constitutive theory for ${}_e{\sigma }^{\left[0\right]}$ describes volumetric deformation physics and the constitutive theory for ${}_d{\sigma }^{\left[0\right]}$ addresses distortional deformation. The conservation and balance laws of classical continuum mechanics [18] [19]: conservation of mass, balance of linear momenta, balance of angular momenta, energy equations and entropy inequality are given in the following.

${\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\right\}}^b-\left[\left[J\right]\left(\left[{}_e{\sigma }^{\left[0\right]}\right]+\left[{}_d{\sigma }^{\left[0\right]}\right]\right)\right]\{\nabla \}=0$ (2)

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

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

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

We use the following equation of state

$p\left(\rho \right)=C\left(\frac{\rho }{{\rho }_0}-1\right)$ (6)

where $C$ is the bulk modulus, $\left\{u\right\}$ are displacements, ${\rho }_0$ is density in the reference or initial configuration, $\left\{F^b\right\}$ are body forces per unit mass, $\left[J\right]$ is deformation gradient tensor, ${}_e{\sigma }^{\left[0\right]}$ and ${}_d{\sigma }^{\left[0\right]}$ are contravariant equilibrium and deviatoric second Piola-Kirchhoff stress tensor, $ϵ$ is permutation tensor, $e$ is specific internal energy, $q$ is heat flux, $\nabla$ is gradient operator, ${\dot{\epsilon }}_{\left[0\right]}$ is rate of Green’s strain tensor, $\varphi$ is Helmholtz free energy density, $\eta$ is entropy density, $\theta$ is absolute temperature, $g$ is temperature gradient tensor, $p\left(\rho \right)$ is thermodynamic pressure and ${\sigma }^{\left(0\right)}$ is contravariant Cauchy stress tensor.

The choice of the equation of state depends upon the specific solid matter under consideration. The present equation of state describes compressibility physics correctly, hence is valid and is used here for illustration purposes to show the influence of compressibility on shock physics. A different equation of state will result in different thermodynamic pressure density relation altering the nature of shock physics but still preserving the shock physics.

3.2. Constitutive Theories

Following references [18] [19], for the nonlinear elastic solids, we have the following constitutive tensors and their argument tensors

${}_d{\sigma }^{\left[0\right]}={}_d{\sigma }^{\left[0\right]}\left({\epsilon }_{\left[0\right]},\theta \right)$ (7)

$q=q\left(g,\theta \right)$ (8)

${}_e{\sigma }^{\left[0\right]}={}_e{\sigma }^{\left[0\right]}\left(\rho ,\theta \right),\text{symbolic}$ (9)

Arguments of ${}_e{\sigma }^{\left[0\right]}$ are symbolic because $\rho$ cannot be an argument tensor in Lagrangian description. We assume that dissipation mechanism is described by ${\epsilon }_{\left[i\right]};i=1,2,\cdots ,n$ , the strain rates i.e. rates of Green’s strain tensor up to orders n and the rheology is due to rates of deviatoric contravariant second Piola-Kirchhoff stress tensor up to orders m i.e. due to ${}_d{\sigma }^{\left[j\right]};j=0,1,\cdots ,m$ . Constitutive tensor and the argument tensors in (7) can now be modified using

$d{\sigma }^{\left[m\right]}={}_d\sigma \left[m\right]\left(\epsilon \left[0\right],\epsilon \left[i\right],d\sigma \left[j\right],\theta \right);i=1,2,\cdots ,n;j=0,1,\cdots ,m-1$ (10)

We remark that $\varphi$ and $\eta$ are not constitutive tensor [18] [19], but are dependent on $\rho$ and $\theta$ .

3.2.1. Constitutive Theory for ${}_e{\sigma }^{\left[0\right]}$

The constitutive theory for ${}_e{\sigma }^{\left[0\right]}$ cannot be derived in Lagrangian description as $\rho$ is not a dependent variable. Following reference [18] [19] we can derive the constitutive theory for ${}_e{\overline{\sigma }}^{\left(0\right)}$ equilibrium Cauchy stress tensor using entropy inequality in Eulerian description.

${}_e{\overline{\sigma }}^{\left(0\right)}=-\overline{p}\left(\overline{\rho },\overline{\theta }\right)\delta$ (11)

$\overline{p}\left(\overline{\rho },\overline{\theta }\right)=-\left({\overline{p}}^2\right)\frac{\partial \overline{\varphi }}{\partial \overline{\rho }}$ (12)

in which $\overline{p}$ is thermodynamic pressure for incompressible matter we have

${}_e{\overline{\sigma }}^{\left(0\right)}=-\overline{p}\left(\overline{\theta }\right)\delta$ (13)

in which $\overline{p}\left(\theta \right)$ is mechanical pressure. In case of Lagrangian description, (11, 12, 13) can be written as (for compressible matter)

${}_e{\sigma }^{\left(0\right)}=-p\left(\rho ,\theta \right)\delta$ (14)

$p\left(\rho ,\theta \right)=-\left(p^2\right)\frac{\partial \varphi }{\partial \rho }$ (15)

For incompressible solid matter

${}_e{\sigma }^{\left(0\right)}=-p\left(\theta \right)\delta$ (16)

From (14) and (16), we can obtain equilibrium contravariant second Piola-Kirchhoff stress tensor ${}_e{\sigma }^{\left[0\right]}$ for compressible and incompressible solid matter.

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

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

Equations (17) and (18) are the constitutive theories for ${}_e{\sigma }^{\left[0\right]}$ for compressible and incompressible solid matter elastic and viscoelastic matter with and without rheology.

3.2.2. Constitutive Theories for ${}_e{\sigma }^{\left[0\right]}$ and $q$

Theories for ${}_d{\sigma }^{\left[m\right]}$ and $q$ have been presented in references [18] [19] using representation theorem and integrity, complete basis of the space of tensor ${}_d{\sigma }^{\left[m\right]}$ and $q$ . We consider simplified constitutive theories to ${}_d{\sigma }^{\left[m\right]}$ and $q$ . We consider a constitutive theory for deviatoric second Piola-Kirchhoff stress tensor based on $n=1$ and $m=1$ in which constitutive theory for ${}_d{\sigma }^{\left[1\right]}$ linear in its argument tensors. We can write this constitutive theory as follows (neglecting initial stress field and thermal i.e. temperature terms).

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

Following reference [18] [19], a linear constitutive theory for $q$ is given by

$q=-kg$ (20)

in which k is thermal conductivity and $g$ is temperature gradient. The reduced form of entropy inequality is given by

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

3.3. Complete Mathematical Model in $R^3$

The conservation and the balance laws (1)-(4), (21) and constitutive theories (17), (19) and (20) constitutive complete mathematical listed below:

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

${\rho }_0\frac{{\partial }^2\left\{u\right\}}{\partial t^2}-{\rho }_0\left\{F_b\right\}-\left[\left[J\right]\left(\left[{}_e{\sigma }^{\left[0\right]}\right]+\left[{}_d{\sigma }^{\left[0\right]}\right]\right)\right]\{\nabla \}=0$ (23)

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

$\begin{array}{l}{\rho }_0\frac{\partial e}{\partial t}+\nabla \cdot q-\left({}_e{\sigma }^{\left[0\right]}+{}_d{\sigma }^{\left[0\right]}\right):{\dot{\epsilon }}_{\left[0\right]}=0\\ \frac{q\cdot g}{\theta }-{}_d{\sigma }^{\left[0\right]}:{\dot{\epsilon }}_{\left[0\right]}\le 0\end{array}$(25)

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

${\sigma }^{\left[0\right]}+{\lambda }_1 \frac{\partial \left({}_d{\sigma }^{\left[0\right]}\right)}{\partial t}={\sigma }_{0}I+2\mu {\epsilon }_{\left[0\right]}+\lambda \left(tr{\epsilon }_{\left[0\right]}\right)I+2{\eta }_1{\epsilon }_{\left[1\right]}+k_1\left(\left(tr{\epsilon }_{\left[1\right]}\right)\right)I$ (27)

$q=-kg$ (28)

This mathematical model is a system of thirteen partial differential equations: balance of linear momenta (3), energy Equation (1), constitutive theories for ${}_d{\sigma }^{\left[0\right]}\left(6\right)$ , and $q\left(3\right)$ in thirteen variables: $u\left(3\right)$ , ${}_d{\sigma }^{\left[0\right]}\left(6\right)$ , $q\left(3\right)$ , $\theta \left(1\right)$ , thus this mathematical model has closure.

3.4. Dimensionless Form of Mathematical Model

When using methods of approximations such as space-time coupled finite element method based on space-time residual functional, it is necessary nondimensionalized the mathematical model to avoid round off errors in the computations. In doing so, we first rewrite (21)-(28) using hat ($ˆ$ ) on all quantitative indicating that they have their usual dimensions or units. We choose reference quantities with subscript zero and define dimensionless variables.

$\begin{array}{lll}x=\frac{\widehat{x}}{L_0},\hfill & u=\frac{\widehat{u}}{L_0},\hfill & {}_e{\sigma }^{\left[0\right]}=\frac{{}_e{\widehat{\sigma }}^{\left[0\right]}}{{\tau }_0}\hfill \\ p=\frac{\widehat{p}}{p_0},\hfill & t_0=\frac{L_0}{v_0},\hfill & v_0=\sqrt{\frac{E_0}{{\left({\rho }_0\right)}_{\text{ref}}}}\hfill \\ E=\frac{\widehat{E}}{E_0},\hfill & \tau =p_0={\left({\rho }_0\right)}_{\text{ref}}{v}_0^2,\hfill & De=\frac{\lambda }{t_0}\hfill \\ F^b=\frac{{\widehat{F}}^b}{F_0},\hfill & {\eta }_1=\frac{{\widehat{\eta }}_1}{{\eta }_0},\hfill & k_1=\frac{{\widehat{k}}_1}{{\eta }_0}\hfill \\ \rho =\frac{\widehat{\rho }}{{\left({\rho }_0\right)}_{\text{ref}}},\hfill & k=\frac{\widehat{k}}{k_0},\hfill & \theta =\frac{\widehat{\theta }}{{\theta }_0}\hfill \\ \hfill & c_v=\frac{{\widehat{c}}_v}{c_{v_0}}\hfill & \hfill \end{array}$(29)

We define

$Ec=\frac{v_0^2}{c_{v_0}{\theta }_0},Re=\frac{{\rho }_0{v}_0{L}_0}{{\eta }_0},Br=\frac{{\eta }_0{v}_0^2}{k_0{\tau }_0}$ (30)

In which Ec, Re, and Br are Ecket’s number, Reynold’s number and Brinckman number, then balance of linear momenta, energy equation, constitutive theories in $R^3$ and the equation of state can be written as (neglecting initial stress ${\sigma }_0$ ).

${\rho }_0\frac{{\partial }^2\left\{u\right\}}{\partial t^2}--{\rho }_0\left\{F_b\right\}-\left[\left[J\right]\left(\left[{}_e{\sigma }^{\left[0\right]}\right]+\left[{}_d{\sigma }^{\left[0\right]}\right]\right)\right]\{\nabla \}=0$ (31)

$\frac{{\rho }_0}{Ec}{c}_v\frac{\partial \theta }{\partial t}+\frac{1}{ReBr}\nabla \cdot q-\left({}_e\sigma +{}_d\sigma \right):{\dot{\epsilon }}_0=0$(32)

$q=-kg$ (33)

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

$\begin{array}{l}{}_d{\sigma }^{\left[0\right]}+\lambda \frac{\partial }{\partial t}\left({}_d{\sigma }^{\left[0\right]}\right)\\ =2\mu {\epsilon }_{\left[0\right]}+\lambda \left(tr{\epsilon }_{\left[0\right]}\right)I+\frac{2{\eta }_1}{Re}\frac{\partial {\epsilon }_{\left[0\right]}}{\partial t}+\frac{k_1}{Re}\left(tr\left(\frac{\partial {\epsilon }_{\left[0\right]}}{\partial t}\right)\right)I\end{array}$ (35)

$p\left(\rho \right)=C\left(\frac{\rho }{{\rho }_0}-1\right)$ (36)

3.4.1. Dimensionless Mathematical Model in $R^1$

For 1D wave physics, the mathematical model assuming compressive pressure to be positive (37)-(42) $R^3$ can be reduced to the following in $R^1$ (considering compressive ${}_d{\sigma }^{\left[0\right]}$ to be positive).

$\begin{array}{l}J=\left|J\right|=J=\left(1+\frac{\partial u_1}{\partial x_1}\right)\\ {\rho }_0\left(x_1\right)=\left|J\right|\rho \left(x_1,t\right)\end{array}$ (37)

${\rho }_0\frac{\partial v_1}{\partial t}-{\rho }_0{F}_1^b+\frac{\partial }{\partial x_1}\left(\left[J\right]\left({}_e{\sigma }_{11}^{\left[0\right]}\right)\right)-\frac{\partial }{\partial x_1}\left(\left[J\right]\left({}_d{\sigma }_{11}^{\left[0\right]}\right)\right)$ (38)

$\frac{p_0{c}_v}{E_c}\frac{\partial \theta }{\partial t}-\frac{1}{ReBr}\frac{\partial q_1}{\partial x_1}-\left({}_e{\sigma }_{11}^{\left[0\right]}+{}_d{\sigma }_{11}^{\left[0\right]}\right)\frac{\partial }{\partial t}\left({\left({\epsilon }_{\left[0\right]}\right)}_{11}\right)=0$ (39)

${}_e{\sigma }_{11}^{\left[0\right]}=\left|J\right|\left(J^{-1}\right)\left(p\left(\rho \right)\right)\left(J^{-1}\right)=\left(J^{-1}\right)p\left(\rho \right)$ (40)

${}_d{\sigma }_{11}^{\left[0\right]}+De\frac{\partial \left({}_d{\sigma }_{11}^{\left[0\right]}\right)}{\partial t}=E{\left({\epsilon }_{\left[0\right]}\right)}_{11}+C_2\frac{\partial }{\partial t}\left({\left({\epsilon }_{\left[0\right]}\right)}_{11}\right)$ (41)

$q_1=-k{g}_1$ (42)

$v_1=\frac{\partial u_1}{\partial t}$ (43)

$p\left(\rho \right)=C\left(\frac{\rho }{{\rho }_0}-1\right)$ (44)

in which $C_2=\frac{\eta }{Re}$ is the dimensionless damping coefficient. Various quantities appearing in (37)-(43) are explicitly defined in the following

$p\left(\rho \right)=C\left(J^{-1}-1\right)$ (45)

and we have the following

${\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$ (46)

Therefore

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

${}_e{\sigma }_{11}^{\left[0\right]}=C\left(J^{-1}\right)\left(J^{-1}-1\right)$ (48)

$J\left({}_e{\sigma }_{11}^{\left[0\right]}\right)=C\left(J^{-1}-1\right)$ (49)

$\frac{\partial }{\partial x_1}\left(J\left({}_e{\sigma }_{11}^{\left[0\right]}\right)\right)=C\frac{\partial }{\partial x_1}\left(J^{-1}\right)$ (50)

and

$\frac{\partial }{\partial x_1}\left(J\left({}_d{\sigma }_{11}^{\left[0\right]}\right)\right)=\frac{\partial J}{\partial x_1}\left({}_d{\sigma }_{11}^{\left[0\right]}\right)+J\left(\frac{\partial \left({}_d{\sigma }_{11}^{\left[0\right]}\right)}{\partial x_1}\right)$ (51)

The mathematical model in $R^1$ (37)-(44) in conjunction with (45)-(51) holds for compressible thermoviscoelastic solids with dissipation and rheology.

3.4.2. Entropy Generation

Consider entropy inequality in $R^1$ (obtained using (5))

${\rho }_0\left(\frac{D\varphi }{Dt}+\eta \frac{D\theta }{Dt}\right)+\frac{q_1\cdot g_1}{\theta }-\left({\sigma }_{11}^{\left[0\right]}\right)\left({\left({\dot{\epsilon }}^{\left[0\right]}\right)}_{11}\right)\le 0$ (52)

in which

$q_1=-k\frac{\partial \theta }{\partial x_1}$ (53)

$g_1=\frac{\partial \theta }{\partial x_1}$ (54)

In the derivation of (52) (See references [18] [19]), we had multiplied throughout by $\theta$ and had changed sign through and, thus the rate of entropy density generation (per unit volume) i.e. specific entropy, due to heat vector and mechanical work, is given by

${\eta }_q=-\frac{k\frac{\partial \theta }{\partial x_1}}{{\theta }^2},{\eta }_w=\frac{\left({\sigma }_{11}^{\left[0\right]}\right)\left({\left({\dot{\epsilon }}_{\left[0\right]}\right)}_{11}\right)}{\theta }$ (55)

Total rate of entropy density generation $\eta$ (per unit volume) i.e. specific entropy, is given

$\eta ={\eta }_q+{\eta }_w$ (56)

In which $q$ in (55) is straight forward. We explain details of ${\eta }_w$ due to mechanical work in the following.

${\eta }_w=\frac{1}{\theta }\left({}_e{\sigma }_{11}^{\left[0\right]}+{}_d{\sigma }_{11}^{\left[0\right]}\right){\left({\dot{\epsilon }}_{\left[0\right]}\right)}_{11}$(57)

substituting ${}_d{\sigma }_{11}^{\left[0\right]}$ from (41) in (52) and expanding, we can write the following.

$\begin{array}{c}{\eta }_w=\frac{1}{\theta }\left(\left({}_e{\sigma }_{11}^{\left[0\right]}\right){\left({\dot{\epsilon }}_{\left[0\right]}\right)}_{11}\right)+\frac{1}{\theta }\left(-De\frac{\partial }{\partial t}\left({}_d{\sigma }_{11}^{\left[0\right]}\right)\right)\\ +\frac{E}{\theta }\left({\epsilon }_{11}^{\left[0\right]}\right){\left({\dot{\epsilon }}_{\left[0\right]}\right)}_{11}+\frac{\eta }{\theta Re}\left({\left({\dot{\epsilon }}_{\left[0\right]}\right)}_{11}\right){\left({\dot{\epsilon }}_{\left[0\right]}\right)}_{11}\end{array}$(58)

and

${\eta }_w={\displaystyle \sum }_{i=1}^4{\eta }_w^i$ (59)

where ${\eta }_w^i,i=1,2,3,4$ in (58) correspond to the corresponding terms in (56). ${\eta }_w^1$ is due to compressibility i.e. due to volumetric change, ${\eta }_w^2$ is due to rheology, ${\eta }_w^3$ is due to elasticity and strain rate and ${\eta }_w^4$ is due to dissipation.

4. Remarks

Thus we observe that rate of entropy density generation is quite complex in polymeric solids. We note the following:

1) In the absence of rheology, ${\eta }_w^2=0$ .

2) In the absence of dissipation and rheology, ${\eta }_w^2=0$ and ${\eta }_w^4=0$ .

3) If the matter is incompressible, then ${\eta }_w^1=0$ . In this case, the stress, strain, and strain rate measures become ${}_d{\sigma }_{11}$ , ${\epsilon }_{11}$ , and ${\dot{\epsilon }}_{11}$ , in which ${}_d{\sigma }_{11}$ is Cauchy stress, ${\epsilon }_{11}$ , ${\dot{\epsilon }}_{11}$ are linear strain and strain rate.

4) We note that rate of entropy density generation ${\eta }_w^3$ is never zero in elastic solids as it is due to strain and strain rate. This aspect is absent in fluent media due to absence of strain.

5) When the medium is inviscid and incompressible only ${\eta }_w^3$ is nonzero.

6) Rate of entropy density $\eta$ and total entropy Ψ.

All conservation and balance laws are independent of the volume of matter, hence can be viewed as being valid for unit volume. $\eta$ , hence ${\eta }_q$ and ${\eta }_w$ represents rate of specific entropy i.e. entropy per unit volume. We integrate $\eta$ over each space-time finite element and then sum them over the elements to obtain Ψ, total entropy as referred and used in the paper.

In the model problem studies we integrate $\eta$ , hence ${\eta }_g$ and ${\eta }_w^i;i=1,2,3,4$ over the discretization of each space-time strip $\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right);i=1,2,\cdots ,n$ obtain entropy Ψ and entropies ${\Psi }_q$ , ${\Psi }_w^i;i=1,2,3,4$ . For ${\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T$ we can write

${\displaystyle {\int }_{{\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T}\eta \text{d}{\Omega }_{xt}}={\displaystyle {\int }_{{\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T}{\eta }_q\text{d}{\Omega }_{xt}}+{\displaystyle {\int }_{{\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T}{\eta }_w\text{d}{\Omega }_{xt}}$(60)

or

$\Psi ={\Psi }_q+{\Psi }_w$ (61)

Substitution of ${\eta }_w$ from (58) in (60) follows. In the model problem studies Ψ is reported.

5. Solutions of Partial Differential Equations in the Mathematical Model

The mathematical model considered for the model problem studies consists of dimensionless partial differential Equations (38)-(43) in dependent variable $u_1,{}_d{\sigma }_{11}^{\left[0\right]},q_1,\theta$ and $v_1$ that are functions of space coordinates $x_1$ and time $t$ . The solutions of the initial value problems defined by these equations is obtained by using space-time coupled finite element method based of space-time residual functional [17] in which the local approximation are ${\left(u_1\right)}_h^e,{\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e,{\left(q_1\right)}_h^e,{\left(\theta \right)}_h^e$ and ${\left(v_1\right)}_h^e$ for the dependent variables $u_1,{}_d{\sigma }_{11}^{\left[0\right]},q_1,\theta$ and $v_1$ are over a space-time element ${\overline{\Omega }}_{xt}^e$ . The local approximations are p-version hierarchical with higher order global differentiability in space as well as time, hence are in $H_{\left({\overline{\Omega }}_{xt}^e\right)}^{\left(p,k\right)}$ scalar product spaces.

Figure 2(a) shows space time domain ${\overline{\Omega }}_{xt}={\overline{\Omega }}_x\times {\overline{\Omega }}_t$ . Figure 2(b) shows discretization ${\overline{\Omega }}_{xt}^T$ of ${\overline{\Omega }}_{xt}$ in space-time strip ${\overline{\Omega }}_{xt}^{\left(i\right)};i=1,2,\cdots ,n$ given by the

Figure 2. Space-time domain and its discretizations using space-time elements and space-time strips. (a) Space-time domain ${\overline{\Omega }}_{xt}$ ; (b) Discretization ${\overline{\Omega }}_{xt}$ into a space-time strip; (c) Discretization of $i^{\text{th}}$ space-time strip ${\overline{\Omega }}_{xt}$ using space-time coupled finite elements.

following.

${\overline{\Omega }}_{xt}^T={\displaystyle \underset{i}{\cup }{\overline{\Omega }}_{xt}^{\left(i\right)}}$ (62)

Consider $i^{th}$ space-time strip discretized using space-time finite elements (Figure 2(c)). Let the local approximations for $u_1,{}_d{\sigma }_{11}^{\left[0\right]},q_1,\theta$ and $v_1$ be given by (equal order equal degree)

$\begin{array}{l}{\left(u_1\right)}_h^e={\displaystyle \sum }_{i=1}^n{N}_i\left(\xi ,\eta \right){\left({\delta }_i^e\right)}_{u_1},{\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e={\displaystyle \sum }_{i=1}^n{N}_i\left(\xi ,\eta \right){\left({\delta }_i^e\right)}_{{}_d{\overline{\sigma }}_{11}^{\left[0\right]}}\hfill \\ {\left(q_1\right)}_h^e={\displaystyle \sum }_{i=1}^n{N}_i\left(\xi ,\eta \right){\left({\delta }_i^e\right)}_{q_1},{\left(\theta \right)}_h^e={\displaystyle \sum }_{i=1}^n{N}_i\left(\xi ,\eta \right){\left({\delta }_i^e\right)}_{\theta }\hfill \\ {\left(v_1\right)}_h^e={\displaystyle \sum }_{i=1}^n{N}_i\left(\xi ,\eta \right){\left({\delta }_i^e\right)}_{v_1}\hfill \end{array}$(63)

Approximations ${\left(u_1\right)}_h^e,{\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e,{\left(q_1\right)}_h^e,{\left(\theta \right)}_h^e$ and ${\left(v_1\right)}_h^e$ over ${\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T$ are given by

$\begin{array}{l}{\left(u_1\right)}_h={\displaystyle \underset{e}{\cup }{\left(u_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};{\left(q_1\right)}_h={\displaystyle \underset{e}{\cup }{\left(q_1\right)}_h^e};\\ {\left(\theta \right)}_h={\displaystyle \underset{e}{\cup }{\left(\theta \right)}_h^e};{\left(v_1\right)}_h={\displaystyle \underset{e}{\cup }{\left(v_1\right)}_h^e};\end{array}$ (64)

in which ${\left\{\delta \right\}}_{u_1},{\left\{\delta \right\}}_{{}_d{\sigma }_{11}^{\left[0\right]}},{\left\{\delta \right\}}_{q_1},{\left\{\delta \right\}}_{\theta }$ and ${\left\{\delta \right\}}_{v_1}$ be degrees of freedom for an element ${\overline{\Omega }}_{xt}^e$ for $u_1,{}_d{\sigma }_{11}^{\left[0\right]},q_1,\theta$ and $v_1$ and the degrees of freedom for $u_1,{}_d{\sigma }_{11}^{\left[0\right]},q_1,\theta$ and $v_1$ for ${\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T$ are given by

$\begin{array}{l}{\left\{\delta \right\}}_{u_1}={\displaystyle \underset{e}{\cup }{\left\{{\delta }^e\right\}}_{u_1}},{\left\{\delta \right\}}_{{ }_d{\sigma }_{11}^{\left[0\right]}}={\displaystyle \underset{e}{\cup }{\left\{{\delta }^e\right\}}_{{}_d{\sigma }_{11}^{\left[0\right]}}},{\left\{\delta \right\}}_{q_1}={\displaystyle \underset{e}{\cup }{\left\{{\delta }^e\right\}}_{q_1}},\\ {\left\{\delta \right\}}_{\theta }={\displaystyle \underset{e}{\cup }{\left\{{\delta }^e\right\}}_{\theta }},{\left\{\delta \right\}}_{v_1}={\displaystyle \underset{e}{\cup }{\left\{{\delta }^e\right\}}_{v_1}}\end{array}$ (65)

The total degrees of freedom for ${\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T$ are given by $\left\{\delta \right\}$

${\left\{\delta \right\}}^T=\left[{\left\{\delta \right\}}_{u_1}^T,{\left\{\delta \right\}}_{{}_d{\sigma }_{11}^{\left[0\right]}}^T,{\left\{\delta \right\}}_{q_1}^T,{\left\{\delta \right\}}_{\theta }^T,{\left\{\delta \right\}}_{v_1}^T\right]$ (66)

substituting ${\left(u_1\right)}_h^e,{\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e,{\left(q_1\right)}_h^e,{\left(\theta \right)}_h^e$ and ${\left(v_1\right)}_h^e$ in Equations (38), (39) and (41)-(43), we obtain residual equations $E_i^e;i=1,2,\cdots ,5$ over ${\overline{\Omega }}_{xt}^e$ .

$E_1^e={\rho }_0\frac{\partial {\left(u_1\right)}_h^e}{\partial t}-{\rho }_0{F}_1^b-\frac{\partial \left(J{\left({}_e{\sigma }_{11}^{\left[0\right]}\right)}_h^e\right)}{\partial x_1}-\frac{\partial \left(J{\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e\right)}{\partial x_1}$ (67)

$E_2^e=\frac{{\rho }_0{c}_v}{Ec}\frac{\partial {\left(\theta \right)}_h^e}{\partial t}-\frac{1}{ReBr}\frac{\partial {\left(q_1\right)}_h^e}{\partial x_1}-\frac{\eta }{Re}\left({\left({}_e{\sigma }_{11}^{\left[0\right]}\right)}_h^e+{\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e\right)\frac{\partial {\left({\left({\epsilon }_{\left[0\right]}\right)}_{11}\right)}_h^e}{\partial t}$ (68)

$E_3^e={\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e+De\frac{\partial {\left({}_d{\sigma }_{11}^{\left[0\right]}\right)}_h^e}{\partial t}-E{\left({\left({\epsilon }_{\left[0\right]}\right)}_{11}\right)}_h^e-C_2\frac{\partial {\left({\left({\epsilon }_{\left[0\right]}\right)}_{11}\right)}_h^e}{\partial t}$ (69)

$E_4^e={\left(q_1\right)}_h^e+k\frac{\partial {\left(\theta \right)}_h^e}{\partial x_1}$ (70)

$E_5^e={\left(v_1\right)}_h^e-k\frac{\partial {\left(u_1\right)}_h^e}{\partial x_1}$ (71)

these hold $\forall x_1\in {\overline{\Omega }}_{xt}^e$ , the domain of space-time element $e$ . We construct space-time residual functional I over ${\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T$ given by

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

In which $E_i^e;i=1,2,\cdots ,5$ are residual equations for an element e with domain ${\overline{\Omega }}_{xt}^e$ given by (67)-(71).

A necessary condition for obtaining an extremum of I is that its first variation must vanish at the extremum provided functional I is differentiable in its arguments. Using (72), we can write

$\begin{array}{c}\delta I={\displaystyle \sum }_e\delta I^e=2{\displaystyle \sum }_{i=1}^5{\left(E_i,\delta E_i\right)}_{{\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T}\\ ={\displaystyle \sum }_e{\displaystyle \sum }_{i=1}^{5}2{\left(E_i^e,\delta E_i^e\right)}_{{\overline{\Omega }}_{xt}^e}={\displaystyle \sum }_e\left\{g^e\right\}=\left\{g\right\}=0\end{array}$(73)

Let

${\left\{{\delta }^e\right\}}^T=\left[{\left\{{\delta }^e\right\}}_{u_1}^T,{\left\{{\delta }^e\right\}}_{{}_d{\sigma }_{11}^{\left[0\right]}}^T,{\left\{{\delta }^e\right\}}_{q_1}^T,{\left\{{\delta }^e\right\}}_{\theta }^T,{\left\{{\delta }^e\right\}}_{v_1}^T\right]$ (74)

and

$\left\{\delta \right\}=\left[{\left\{\delta \right\}}_{u_1},{\left\{\delta \right\}}_{{}_d{\sigma }_{11}^{\left[0\right]}},{\left\{\delta \right\}}_{q_1},{\left\{\delta \right\}}_{\theta },{\left\{\delta \right\}}_{v_1}\right]$ (75)

In which $\left\{g^e\right\}$ is a nonlinear function of $\left\{{\delta }^e\right\}$ . Following Surana et al. [17] we find a solution $\left\{\delta \right\}$ that satisfies $\left\{g\right\}=0$ iteratively using Newton’s linear method with line search. Details are given in the following. Let ${\left\{\delta \right\}}_0$ be an assumed or starting solution, then

$\left\{g\left({\left\{\delta \right\}}_0\right)\right\}\ne 0$ (76)

Let $\left\{\text{Δ}\delta \right\}$ be corrections to ${\left\{\delta \right\}}_0$ such that

$\left\{g\left({\left\{\delta \right\}}_0+\left\{\text{Δ}\delta \right\}\right)\right\}=0$ (77)

We expand (77) in Taylor Series about ${\left\{\delta \right\}}_0$ and retain only up to linear term in $\left\{\text{Δ}\delta \right\}$ (Newton’s linear method or Newton-Raphson method).

$\left\{g\left({\left\{\delta \right\}}_0+\left\{\text{Δ}\delta \right\}\right)\right\}=\left\{g\left({\left\{\delta \right\}}_0\right)\right\}+{\left[\frac{\partial \left\{g\right\}}{\partial \left\{\delta \right\}}\right]}_{{\left\{\delta \right\}}_0}\left\{\text{Δ}\delta \right\}=0$ (78)

$\therefore \left\{\text{Δ}\delta \right\}=-{\left[\frac{\partial \left\{g\right\}}{\partial \left\{\delta \right\}}\right]}_{{\left\{\delta \right\}}_0}^{-1}\left\{g\left({\left\{\delta \right\}}_0\right)\right\}$ (79)

in which

$\left[\frac{\partial \left\{g\right\}}{\partial \left\{\delta \right\}}\right]=\delta \left\{g\right\}=\delta \left(\delta I\right)={\delta }^{2}I$ (80)

which is second variation of I.

Using (73)

$\begin{array}{c}{\delta }^{2}I=\delta \left\{g\right\}={\displaystyle \sum }_{i=1}^{5}2\left({\left(\delta E_i,\delta E_i\right)}_{{\overline{\Omega }}_{xt}}+{\left(E_i,{\delta }^2{E}_i\right)}_{{\overline{\Omega }}_{xt}}\right)\\ ={\displaystyle \sum }_e\left({\displaystyle \sum }_{i=1}^{5}2\left({\left(\delta E_i^e,\delta E_i^e\right)}_{{\overline{\Omega }}_{xt}}+{\left(E_i^e,{\delta }^2{E}_i\right)}_{{\overline{\Omega }}_{xt}^e}\right)\right)={\displaystyle \sum }_e\delta \left\{g^e\right\}\end{array}$(81)

Surana et al. [17] have shown that a sufficient condition or an extremum principle is possible if ${\delta }^2{E}_i,{\delta }^2{E}_i^e$ is neglected in (81). Justification and validity of this approximation are described in reference [17]. Thus, now we have

${\delta }^{2}I\equiv {\displaystyle \sum }_e\left({\displaystyle \sum }_{i=1}^{5}2{\left(\delta E_i^e,\delta E_i^e\right)}_{{}_{\overline{\text{Ω}}}^{xt}}\right)>0\text{for}\text{all}{E}_i^e\ne 0;i=1,2,\cdots ,5$(82)

or

${\delta }^{2}I={\displaystyle \sum }_e\left[K^e\left(\left\{{\delta }^e\right\}\right)\right]=\left[K\right]$ (83)

where

$\left[K^e\left(\left\{{\delta }^e\right\}\right)\right]={\displaystyle \sum }_{i=1}^{5}2{\left(\delta E_i^e,\delta E_i^e\right)}_{{\overline{\Omega }}_{xt}^e}$(84)

We calculate $\left\{\text{Δ}\delta \right\}$ using

$\left\{\text{Δ}\delta \right\}=-{\left[K\right]}_{{\left\{\delta \right\}}_0}^{-1}\left\{g\left(\left\{{\delta }_0\right\}\right)\right\}$ (85)

and the improved solution is given by (using line search [17])

$\left\{\delta \right\}={\left\{\delta \right\}}_0+{\alpha }^{\text{*}}\left\{\text{Δ}\delta \right\}$ (86)

In which ${\alpha }^{\text{*}}$ is determined using the smallest value of I.

$I\left(\left\{{\delta }_0\right\}+\alpha \left\{\Delta \delta \right\}\right)<I\left(\text{of}\text{solution}\right)\forall \alpha \in \left(0,2\right).$ (87)

This is referred to as line search. We check for convergence i.e. when ${\left|g_{{ }_i}\right|}_{\max }\le \Delta$ , $\Delta =\mathcal{O}\left({10}^{-6}\right)$ or lower the solution $\left\{\delta \right\}$ is the solution.

Main steps in computations:

1) Assume ${\left\{\delta \right\}}_0$ , starting or assumed solution

2) Calculate $\left\{g^e\right\}={\displaystyle {\sum }_{i=1}^5\left(E_i^e,\delta E_i^e\right)}=\left\{f^e\right\}$

3) Assemble $\left\{g^e\right\}$ to obtain $\left\{g\right\}={\displaystyle {\sum }_e\left\{g^e\right\}}$

4) Calculate ${\delta }^2{I}^e={\displaystyle {\sum }_{i=1}^5{\left(\delta E_i^e,\delta E_i^e\right)}_{{\overline{\Omega }}_{xt}^e}}$

5) Assemble ${\delta }^2{I}^e$ to obtain ${\delta }^{2}I={\displaystyle {\sum }_e{\delta }^2{I}^e}={\displaystyle {\sum }_e\left[K^e\right]}$

6) Calculate $\left\{\text{Δ}\delta \right\}$ using $\left\{\text{Δ}\delta \right\}=-{\left[{\delta }^{2}I\right]}_{{\left\{\delta \right\}}_0}^{-1}\left\{g\left({\left\{\delta \right\}}_0\right)\right\}$

7) Calculate ${\alpha }^{\text{*}}$

8) Calculate improved solution $\left\{\delta \right\}$ using

$\left\{\delta \right\}={\left\{\delta \right\}}_0+{\alpha }^{\text{*}}\left\{\text{Δ}\delta \right\}$ (88)

9) Calculate $\left\{g\left(\left\{\delta \right\}\right)\right\}$

10) Check for convergence of Newton’s linear method

${\left|g_i\right|}_{\max }\le \Delta$ , a present tolerance of computed zero(89)

$\text{Δ}=\mathcal{O}\left({10}^{-6}\right)$ or lower is generally satisfactory. If converged then $\left\{\delta \right\}$ in (88) is the converged solution. If not converged, then set ${\left\{\delta \right\}}_0=\left\{\delta \right\}$ and repeat steps 2 - 10 till converged.

6. Significant Aspects of the Computational Method Used in the Present Work

The space-time coupled finite element method using a space-time strip or a slab is most meritorious for obtaining solutions of initial value problems [17], hence is used in the present work.

1) Concurrent dependence of the solution on space and time as required by the physics in IVPs is preserved in this methodology.

2) The space-time integral form resulting from this approach is space-time variationally consistent [17], hence computations are unconditionally stable.

3) Use of higher order, higher degree p-version hierarchical local approximations in higher order space-time scalar product spaces ensure that the space-time integrals over the space-time discretization are always Riemann with the choice of minimally conforming (or orders higher than minimally conforming) spaces. Thus when the integrated sum of squares of the space-time residual (I) over the whole discretization approaches zero the PDEs in the mathematical model are satisfied in the point wise sense (i.e. everywhere) over the whole discretized space-time domain. Thus, proximity of I to zero is a measure error in the solution and when I is $\mathcal{O}\left({10}^{-8}\right)$ or lower, the computed solution is as good as the theoretical solution.

4) Use of space-time strip with time marching has two main advantages: 1) for large values of time, solutions are obtained by solving a very small problem (in terms of degrees of freedom) for each time strip, thus computationally much more efficient than the space-time mesh for ${\overline{\Omega }}_{xt}$ . 2) For each space-time strip only upon obtaining a converged solution, the solution is marched to the next space-time strip, thus ensuring accurate solution for each space-time strip and therefore for the entire space-time domain.

7. Model Problem Studies

We consider the following 1D wave studies in this paper for incompressible and compressible solid medium.

1) Incompressible isothermal: solid medium with and without dissipation, but no rheology.

2) Incompressible isothermal: elastic solid medium, with dissipation and rheology.

3) Compressible isothermal and nonisothermal: solid medium without dissipation and rheology.

4) Compressive nonisothermal: solid medium with and without dissipation, but no rheology.

5) Compressive nonisothermal: solid medium with dissipation and rheology.

6) Compressive nonisothermal: wave propagation, transmission, reflection, and interaction: studies with bimaterial interface of same and different materials.

Numerical solutions of the partial differential equations are obtained using space-time coupled finite element method for a space-time strip with time marching. Figure 3(a) shows a schematic of the rod of uniform cross section completely clamped at the left end ($x_1=0$ ) and is subjected to a compressive velocity pulse of duration $2\text{Δ}t$ and peak values $v^{\text{*}}$ (Figure 3(c)). Figure 3(b) shows mathematical idealization of the rod (described by line) assuming that all points in each cross section of the rod displace in the $x_1$ direction by the same amount. Figure 3(d) shows

Figure 3. Schematic, mathematical idealization, a space-time strip, boundary conditions and initial conditions. (a) Schematic of axial bar or rod; (b) Mathematical idealization of (a); (c) A 30 element uniform discretization of ${\overline{\Omega }}_{}^1{}_{xt}$ , first space-time strip; (d) A Compressive velocity pulse of duration 2Δt and peak −υ*.

uniform discretization of the first space-time strip ${\overline{\Omega }}_{xt}^{\left(1\right)}=\left[0,L\right]\times \left[0,\Delta t\right]$ using a thirty nine node p-version hierarchical space-time finite elements with higher order global differentiability. When $q_1$ is substituted in energy equation and $v_1$ is not used as auxiliary variable, the mathematical model consists of partial differential equation containing up to second order derivative of $u_1$ and T in space $x_1$ , but only first order time derivatives. Thus $k_1=3$ and $k_2=2$ are minimally confirming order of Hilbert Space $H_{{\overline{\Omega }}_{xt}^e}^{\left(k_1,k_2\right)}$ in space and time but use of $k_1=k_2=3$ is permissible without detrimental effects. For this choice, all space-time integrals over the discretization ${\left({\overline{\Omega }}_{xt}^{\left(1\right)}\right)}^T$ of space time strip ${\overline{\Omega }}_{xt}^{\left(1\right)}$ are Riemann. $k_1$ and $k_2$ higher than 3 are permissible as well. When $k_1=k_2=k=2$ then space-time integrals over ${\left({\overline{\Omega }}_{xt}^{\left(1\right)}\right)}^T$ are Lebesgue in space but Riemann in time. When the solutions are analytical this is permissible. The mathematical model used here in the calculations ((38)-(43)) is a system of first order equations for which integrals over ${\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T$ are Riemann when $k_1=k_2=2$ , hence used in the present work. For the compressive velocity pulse we have:

$\begin{array}{l}{v}_1=0,\frac{\partial v_1}{\partial t}=0\text{at}t=0\\ v_1=v^{\text{*}},\frac{\partial v_1}{\partial t}=0\text{at}t=\Delta t\\ v_1=0,\frac{\partial v_1}{\partial t}=0\text{at}t=2\Delta t\\ v_1=0,\frac{\partial v_1}{\partial t}=0\text{for}t\ge 2\Delta t\end{array}$ (90)

In all numerical studies that follow, we always use compressive velocity pulse of $v^{\text{*}}=0.1$ defined by (90) at $x_1=L$ . We choose $\text{Δ}t=0.1$ . The solution computations are initiated for the first space-time strip with the 30 space-time element uniform discretization. Due to smoothness of the solution $k=2$ i.e. solutions of class $C^1$ in space and time yield accurate results upon convergence. We have confirmed this using $k=3$ . Thus, in the studies presented here, we use $k_1=k_2=2$ solution of class $C^1$ in space and time. For the first space-time strip boundary conditions and the initial condition shown in Figure 3(c) are imposed and a convergence study is conducted starting with p-level of 3 and progressively increasing the p-level. At p-level of 7 in space and time the residual functional $\mathcal{O}\left({10}^{-8}\right)$ or lower is achieved indicating the solution is converged. Further increase in p-level does not result is measurable reduction in I as well as no appreciable improvement in the solution. With the tolerance Δ of $\mathcal{O}\left({10}^{-6}\right)$ , $\left|g_i\right|\le \text{Δ}$ is achieved in Newton’s linear method in 3 to 5 iterations. Upon obtaining converged solution for the first space-time strip, the solution is time marched to the second space time strip using initial conditions from the converged solution for the first space time strip, thus ensuring converged solution for each space-time strip, hence for the entire space-time domain. In the numerical studies rate of entropy density generation $\eta$ is integrated over ${\left({\overline{\Omega }}_{xt}^{\left(i\right)}\right)}^T$ to obtain entropy Ψ reported in the graphs.

In the numerical studies, we have used $1/Re=0.002585$ , $Br=0.2$ , $E_c=0.4$ , $c_v=2.5$ , $k=0.1$ and bulk modulus of 0.333.

7.1. Incompressible Solid Medium with and without Dissipation, in the Absence of Rheology (De = 0.0)

This is a linear wave propagation study in purely elastic medium and the elastic medium with dissipation. Due to incompressibility, density remains unaffected during wave propagation. Entropy generation due to mechanical work is not monitored, hence the studies presented here are isothermal as the energy equation is not part of the mathematical model. In this study equation of state is not needed either. We choose $C_2=0.0$ and 0.005 ($De=0.0$ ).

The main purpose of including this study is to compare this wave propagation physics with the wave studies in compressible solid medium to demonstrate sharp contrast between the two. Figure 4 shows the ${}_d{\sigma }_{11}$ stress waves in the rod for different values of time. Since the wave speed is one, the waves reaches the impermeable boundary at $x_1=0$ in ten time increments. The wave reflection at $t=11\text{Δ}t$ is shown in Figure 4(c), the peak magnitudes double. At $t=12\text{Δ}t$ , the waves recover their original shape and begin to propagate toward $x_1=1$ (Figure 4(e) and Figure 4(f)). In the absence of dissipation, the peak and the base of the wave do not change, as expected. In the presence of dissipation, continuous reduction in its peak and the base elongation are observed during the wave propagation. Very slight differences in wave speed are observed after reflection, too small to be measurable. In the presence of dissipation wave remains symmetric with respect to its peak. Important points to note are:

1) Due to fixed wave speed (as there is no change in density) wave shape is not distorted during propagation.

2) Dissipation results in progressive base elongation and amplitude decay.

3) In absence of dissipation, wave remains completely unaffected during propagation.

4) We do not observe formation of shock fronts in the waves during propagation due to constant density.

7.2. Incompressible Isothermal: Solid Medium with Dissipation and Rheology

This case is also linear wave propagation study under isothermal conditions. The purpose of this study is same as described in Section 7.1 i.e. to compare these results with compressible wave physics. In this study, we consider damping coefficient $C_2=0.005$ for $De=0,0.002,0.004$ . Figure 5 shows plots of Cauchy Stress ${}_d{\sigma }_{11}^{\left(0\right)}$ versus $x_1$ for different values of time, $t=2\Delta t,6\Delta t,11\Delta t,12\Delta t,14\Delta t$ and $18\Delta t$ . For $De=0.002$ increased peak value and reduced base of the wave compare to $De=0.0$ is observed during the entire evolution. For $De=0.004$ even higher amplitude and reduced base compared to $De=0.002$ is clearly observed during the entire evolution, Slightly reduced wave speed is observed after reflection when $De\ne 0$ clearly seen for $De=0.004$ (Figure 5(e) and Figure 5(f)). With increasing De, the increased peak of ${}_d{\sigma }_{11}$ are due to additional stresses present in the material that is not yet relaxed. Increasing peak with increasing De is due to increased rheology (higher relaxation time, hence higher De). Important observations are:

1) Increasing De i.e. increasing relaxation time results in increasing amplitudes of ${}_d{\sigma }_{11}^{\left(0\right)}$ waves but with progressively reducing base of the wave.

2) Symmetry of the wave about the peak is preserved in the presence of dissipation and rheology.

3) Wave speed is effected with increasing De, but the change is not measurable in the present study. Computations of evolution for high De and lower damping can be performed to study increasing wave speed with increasing De. Since

Figure 4. Evolution of ${}_d{\sigma }^{\left(0\right)}$ : Incompressible solid medium with damping ($C_2=0.005$ ) and without damping ($C_2=0.0$ ). (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 5. Evolution of ${}_d{\sigma }^{\left(0\right)}$ : Incompressible solid medium with damping ($C_2=0.005$ ) and rheology ($De=0,0.002,0.004$ ). (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

rheology adds additional elasticity due to the long chain molecules, increasing wave speed for increasing De is to be expected.

7.3. Compressible Isothermal and Nonisothermal, Elastic Solid Medium, without Dissipation and Rheology

In this case, we study wave physics in compressible solid matter. In the isothermal study, the entropy density generation is present but is not monitored, hence the energy equation was not part of the mathematical model. In the nonisothermal case the rate of entropy density generation due to the rate of mechanical work and other sources is monitored by including the energy equation as integral part of the mathematical model. In the nonisothermal case material coefficients are not function of temperature, hence the mechanical deformation and thermal field are decoupled implying that the deformation physics of the wave (stress ${}_d{\sigma }_{11}^{\left[0\right]}$ , density $\rho$ ) must be same in the present study as in isothermal studies. Entropy density generation $\eta$ results in thermal field causing changes in temperature of the medium.

Figure 6 shows evolution of ${}_d{\sigma }_{11}^{\left[0\right]}$ along the length of the rod, for both isothermal and nonisothermal studies. Figure 7 shows evolution of $\rho$ , for both isothermal and nonisothermal studies. Results from isothermal and nonisothermal studies compare quite well. This is expected because the deformation physics and the thermal physics are decoupled in the present study. Figure 8 shows evolution of temperature for the nonisothermal case.

From ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ we observe progressive steepening of the wave behind the peak (as explained in the introduction) during propagation. In Figure 6(b), at $t=6\text{Δ}t$ the steep front behind the peak is clearly seen. This steepening of the wave is referred to as shock front of the wave, hence the name “shock wave”. Upon reflection the wave direction reverses shown at $t=11\text{Δ}t$ in Figure 6(c), the shock front remains behind the peak of the wave. In Figure 6(d) and Figure 6(e) we clearly observe shock fronts behind the peak of the wave. In this study due to absence of dissipation, the oscillations in the vicinity of the shock front are physical due to the vibration of material particles as there is no mechanism of dampening them. With time the vibrational energy in the oscillations is transferred to the neighboring particles as the shock front advances. From Figure 6(e) the oscillation between $x_1=0.0\text{-}0.4$ , are not present in Figure 6(f) as the shock format is now located near $x_1=0.8$ , far removed from $x_1=0.0\text{-}0.4$ .

The density and the temperature evolutions follow evolution of ${}_d{\sigma }_{11}^{\left[0\right]}$ , progressive formation of shock front behind the peak, oscillations only in the vicinity of the shock front, the shock front behind the peak of the wave before and after reflection etc. (Figure 7 and Figure 8).

Figure 9(a) shows plots of entropy generation Ψ for each space-time strip during evolution. Figure 9(b) shows entropy generation Ψ as a function of time.

From Figure 9(a), we note Ψ increases from zero to a value at the end of the first time strip, remains constant from first to second time step i.e. the same

Figure 6. Evolution of ${}_d{\sigma }^{\left[0\right]}$ : Compressible solid medium without damping or memory ($C_2=0.0$ , $De=0.0$ ), isothermal and nonisothermal. (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 7. Evolution of $\rho$ : Compressible solid medium without damping or memory ($C_2=0.0$ , $De=0.0$ ), isothermal and nonisothermal. (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 8. Evolution of $\theta$ : Compressible solid medium without damping or memory ($C_2=0.0$ , $De=0.0$ ) isothermal and nonisothermal. (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 9. Entropy Ψ versus time step number and entropy Ψ versus time $t$ . (a) Entropy Ψ versus time step number; (b) Entropy Ψ versus time $t$ .

entropy generation in the second time step as the first time step. Ψ is zero for the third time step until the wave reaches the impermeable boundary at $x_1=0$ , between 9^th and 10^th time steps the rate of entropy generation Ψ increases sharply due to reflection followed by slight increase from time step 10 to 11. When the wave recovers its incident shape with shock front behind the peak. After reflection the entropy generation Ψ drops back to zero and remains zero there after until the wave reaches at $x_1=1.0$ . The cumulative entropy Ψ as a function of time in Figure 9(b) shows continuous increase in Ψ during first two time steps during which the wave enters the medium. Ψ remains unchanged for $t=2\text{Δ}t$ to $t=9\text{Δ}t$ . From $t=9\text{Δ}t$ to $t=11\text{Δ}t$ , the wave reflection and recovery results in continuous rise in Ψ, for $t>11\text{Δ}t$ , Ψ remains unchanged as there is no conversion of mechanical energy into rate of entropy.

Since the deformation field and the thermal physics are decoupled in the studies presented in the paper, the results from the two studies should match, but we observe slight differences in the results in Figure 6(c), peaks at $\left(0.08333,-0.087189\right)$ vs $\left(0.075000,-0.089051\right)$ and similarly in Figure 7(c) $\rho =1.100910$ vs $\rho =1.103080$ . These are attributed due to:

The two sets of results are obtained using two different mathematical models due to inclusion of energy equation in the nonisothermal case. Different number of equations, different round-off errors in computations and use of same tolerance in both cases to check convergence of Newton’s linear method can result in such small deviations. Computational studies with larger word size and very low tolerances for convergence of Newton linear method using higher p and k can be performed to confirm this. This has been confirmed in our studies presented in our earlier works.

7.4. Compressible Nonisothermal: Solid Medium with and without Dissipation, No Rheology (De = 0.0)

In this case, we consider wave study in compressible solid medium without dissipation ($C_2=0.0$ ) and also with dissipation ($C_2=0.005$ ), but no rheology. Figures 10-12 show evolutions of ${}_d{\sigma }_{11}^{\left[0\right]}$ , $\rho$ , $\theta$ along the length of the rod. Figure 13 shows entropy generation Ψ for each time step and also as a function of time. From Figure 10, we note that evolution of ${}_d{\sigma }_{11}^{\left[0\right]}$ for $C_2=0.0$ contains sharp shock fronts and oscillations in the vicinity of the shock front. The evolution of ${}_d{\sigma }_{11}^{\left[0\right]}$ with damping ($C_2=0.005$ ) dampens the oscillations due to the mechanical energy of oscillations being converted into entropy, as a result ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ is smooth (oscillation free) for all values of time. Due to dissipation, there is appreciable amplitude decay and non-symmetric base elongation of the ${}_d{\sigma }_{11}^{\left[0\right]}$ as a consequence of varying wave speed along the base of the wave.

We observe exactly identical behavior of the evolution of density $\rho$ and temperature $\theta$ along the length of the rod (Figure 11 and Figure 12). In Figure 12, lower peaks of $\theta$ for $C_2=0.005$ result in higher $\theta$ behind and ahead of the wave due to conduction. In compressible isothermal studies the thermal physics is present but not monitored due to absence of energy equation in the mathematical model.

Figure 13(a) shows plots of Ψ versus time step number i.e. entropy Ψ generation for each time step. Figure 13(b) presents Ψ versus time $t$ during the evolution. Entropy graphs follow ${}_d{\sigma }_{11}^{\left[0\right]}$ evolution. Higher peak of ${}_d{\sigma }_{11}^{\left[0\right]}$ results in

Figure 10. Evolution of ${}_d{\sigma }^{\left[0\right]}$ : Compressible solid medium with damping ($C_2=0.0$ , $C_2=0.005$ ) and without memory ($De=0.0$ ). (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 11. Evolution of $\rho$ : Compressible solid medium with damping ($C_2=0.0$ , $C_2=0.005$ ) and without memory ($De=0.0$ ). (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 12. Evolution of $\theta$ : Compressible solid medium with damping ($C_2=0.0$ , $C_2=0.005$ ) and without memory ($De=0.0$ ). (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 13. Entropy Ψ versus time step number and entropy Ψ versus time $t$ . (a) Entropy Ψ versus time step number; (b) Entropy Ψ versus time $t$ .

lower dissipation, hence lower entropy Ψ. From Figure 13(a) we note that even between second and third time step Ψ is higher for $C_2=0.005$ compared to $C_2=0.0$ . From time step 3 to 8 there is continuous production of Ψ when $C_2=0.005$ but Ψ is zero for the same time steps when $C_2=0.0$ . Since $C_2=0.005$ results in continuous amplitude decay, by the time the wave reaches the impermeable boundary its peak is substantially reduced (implying reduced mechanical work), hence we see Ψ substantially lower for $C_2=0.005$ compared to $C_2=0.0$ in Figure 13(a). In Figure 13(b), Ψ versus $t$ for $C_2=0.0$ is same as Ψ versus $t$ in Figure 9(b). We note that Ψ versus $t$ in Figure 13(b) is always higher for $C_2=0.005$ compared to $C_2=0.0$ for $0<t<11\text{Δ}t$ . Strong reflection for $C_2=0.0$ (due to higher amplitude of ${}_d{\sigma }_{11}^{\left[0\right]}$ for $C_2=0.0$ ) results in higher value of Ψ at $t=11\text{Δ}t$ but remains constant there after due to $C_2=0.0$ , but Ψ versus $t$ for $C_2=0.005$ continuously increases for $t>11\text{Δ}t$ till the wave reached $x_1=1.0$ .

7.5. Compressible Nonisothermal: Solid Medium with Dissipation and Rheology

We consider compressible thermoviscoelastic solid medium with damping coefficient $C_2=0.005$ and Deborah number $De=0.0,0.002,0.004$ to demonstrate the influence of progressively increasing rheology on shock physics for fixed dissipation. Figures 14-16 show evolutions of ${}_d{\sigma }_{11}^{\left[0\right]}$ , $\rho$ , $\theta$ along the length of the rod. Figure 14 presents details of entropy Ψ during the evolution.

From Figure 14(b), we note that the evolution of stress wave ${}_d{\sigma }_{11}^{\left[0\right]}$ is free of oscillations for all three Deborah numbers due to presence of viscosity but at the expense of diffusing the sharp shock fronts compared to $C_2=0.0$ . Shock formation is always behind the peak of the wave, before and after reflection (Figure 14). Increasing rheology for increasing De produces more pronounced resident stresses that need to be relaxed, hence resulting in progressively increasing amplitude, progressively reducing base and steeper shock front behind the wave.

Evolutions of density $\rho$ and temperature $\theta$ follow exactly the same pattern a evolution of ${}_d{\sigma }_{11}^{\left[0\right]}$ (Figure 15 and Figure 16). From Figure 16, we note that higher temperature peaks in the evolution have slower temperature behind and in front of the wave due to higher ${}_d{\sigma }_{11}^{\left[0\right]}$ caused by lower stress relaxation associated with higher De.

Entropy generation Ψ versus time step number during the evolution is presented in Figure 17(a). Figure 17(b) shows entropy production as a function of time. From Figure 17(a) we note that entropy generation Ψ for each time step is highest when $De=0.0$ and is progressively reduced for progressively increasing De due to lowest stress peaks at $De=0.0$ that increases with increasing De for time steps one to nine. At reflection highest amplitude of ${}_d{\sigma }_{11}^{\left[0\right]}$ for $De=0.004$ results in largest entropy generation (time steps 9 to 12). Graphs of Ψ versus $t$ in Figure 17(b) follow Figure 17(a) i.e. Ψ is largest for all time values for $De=0.0$ but after reflection Ψ for $De=0.004$ dominates i.e. is largest followed by Ψ for $De=0.002$ and $De=0.0$

7.6. Wave Physics with Bimaterial Interface

In the two studies presented here, we choose $C_2=0.0$ and $C_2=0.0009$ and $De=0.0$ . Lower values of $C_2$ maintain the sharp shock fronts in the waves. Choice of De in this study does not matter as we are studying wave propagation,

Figure 14. Evolution of ${}_d{\sigma }^{\left[0\right]}$ : Compressible solid medium with damping ($C_2=0.005$ ) and memory ($De=0.0$ , $De=0.0002$ , $De=0.004$ ). (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 15. Evolution of $\rho$ : Compressible solid medium with damping ($C_2=0.005$ ) and memory ($De=0.0$ , $De=0.0002$ , $De=0.004$ ). (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 16. Evolution of $\theta$ : Compressible solid medium with damping ($C_2=0.005$ ) and memory ($De=0.0$ , $De=0.0002$ , $De=0.004$ ). (a) $t=1\text{Δ}t$ ; (b) $t=6\text{Δ}t$ ; (c) $t=11\text{Δ}t$ ; (d) $t=12\text{Δ}t$ ; (e) $t=14\text{Δ}t$ ; (f) $t=18\text{Δ}t$ .

Figure 17. Entropy Ψ versus time step number and entropy Ψ versus time $t$ . (a) Entropy Ψ versus time step number; (b) Entropy Ψ versus time $t$ .

transmission, reflection, and interaction.

In both studies we choose a rod of 2 unit length, the bimaterial interface is at the center of the rod. Both ends of the rod are subjected to compressive velocity pulse of peak value $v^{\text{*}}=0.1$ and $2\text{Δ}t$ base or support. The rode is discretized using 60 space-time p-version hierarchical space-time finite elements with higher order global differentiability in space and time (Surana et al. [17]). As is all other studies, we choose $p_{\xi }=p_{\Psi }=1$ in space and time with $k_1=k_2=2$ i.e. solution of class $C^1$ in space and time. We labeled the material $M_1$ and $M_2$ to the left and the right of the interface only to identify modulus of elasticity, that is $M_1$ and $M_2$ denote elasticity modulus of the materials (only) of the rod between $0\le x_1\le 1$ and $1\le x_1\le 2$ and do not refer to damping coefficient or Deborah number. In the first study, we choose $M_1=M_2$ and we use the same material coefficients for $M_1$ and $M_2$ as used in earlier studies i.e. $E=1$ , and reference density $\rho =1$ etc. Damping coefficient $C_2$ and $De$ used are shown in the graphs of the results for the two halves of the rod. In this case wave in $M_1$ and $M_2$ are identical. Figures 18-21 show evolutions of ${}_d{\sigma }_{11}^{\left[0\right]}$ , $\rho$ , $\theta$ and entropy Ψ. Waves in $M_1$ and $M_2$ propagate toward the interface. Figures 18(b)-20(b) show the waves at the interface. Since $M_1=M_2$ , in this study there is only wave transmission without reflection. Figure 18(c), Figure 18(d)-Figure 20(c), Figure 20(d) show interaction of the waves at the interface. After transmission the waves resume their own identities and propagate toward $x_1=0$ and $x_1=2$ (Figure 18(e), Figure 18(f)-Figure 20(e), Figure 20(f)). Figure 21(a) and Figure 21(b) show entropy generation Ψ as a function of time step number and time $t$ . For the first two time steps i.e. $0\le t\le 2\text{Δ}t$ , there is constant rise in Ψ as the waves enter the medium followed by entropy production only due to the wave in $1\le x_1\le 2$ as it is propagates in the medium with dissipation resulting in continuous rise in Ψ in Figure 21(b). Between time steps 9 - 11, wave interaction results in substantial rise in Ψ followed by slow increase due to the medium on the right of the interface as it is viscous.

In the second study we choose $E=1$ in $M_1$ and $E=0.64$ in $M_2$ with reference density $\rho =1$ for both materials. This gives us wave speed of one in material $M_1$ and wave speed of $0.8$ in material $M_2$ . We choose $C_2=0.0009$ and $De=0.0$ in both materials. Figures 22-25 show evolutions of ${}_d{\sigma }_{11}^{\left[0\right]}$ , $\rho$ , $\theta$ and Ψ. Consider evolution of ${}_d{\sigma }_{11}^{\left[0\right]}$ in Figure 22. Figure 22(a) and Figure 22(b) show stress waves at $t=2\text{Δ}t$ and $t=8\text{Δ}t$ , wave in $M_1$ is propagating faster than the wave in $M_2$ . In Figure 22(b) the wave in $M_1$ is closer to the interface compared to the wave in $M_2$ . In Figures 22(c)-(e), we see the transmission, reflection, and interaction of the two waves. Figure 22(f) shows final waves in the two media propagating toward $x_1=0$ and $x_1=2$ . The density and temperature waves follow exact same pattern (Figure 23 and Figure 24) as the stress wave. Entropy generation Ψ graphs shown in Figure 25(a) and Figure 25(b) follow stress graphs in Figure 22 and are similar to entropy generation graphs shown in Figure 21.

8. Summary and Conclusions

The work presented in this paper focuses on investigating shock physics in compressible elastic and viscoelastic solid matter with and without rheology and monitoring of rate of entropy generation and associated thermal physics. Due to nonlinear deformation (finite strain finite deformation) in the shock physics

Figure 18. Evolution of ${}_d{\sigma }_{11}^{\left[0\right]}$ : Compressible solid medium. (a) $t=2\text{Δ}t$ ; (b) $t=9\text{Δ}t$ ; (c) $t=10\text{Δ}t$ ; (d) $t=11\text{Δ}t$ ; (e) $t=12\text{Δ}t$ ; (f) $t=14\text{Δ}t$ .

Figure 19. Evolution of $\rho$ : Compressible solid medium. (a) $t=2\text{Δ}t$ ; (b) $t=9\text{Δ}t$ ; (c) $t=10\text{Δ}t$ ; (d) $t=11\text{Δ}t$ ; (e) $t=12\text{Δ}t$ ; (f) $t=14\text{Δ}t$ .

Figure 20. Evolution of $\theta$ : Compressible solid medium. (a) $t=2\text{Δ}t$ ; (b) $t=9\text{Δ}t$ ; (c) $t=10\text{Δ}t$ ; (d) $t=11\text{Δ}t$ ; (e) $t=12\text{Δ}t$ ; (f) $t=14\text{Δ}t$ .

Figure 21. Entropy Ψ versus time step number and entropy Ψ versus time $t$ . (a) Entropy Ψ versus time step; (b) Entropy Ψ versus time $t$ .

Figure 22. Evolution of ${}_d{\sigma }_{11}^{\left[0\right]}$ : Compressible solid medium, $C_2=0.0009$ , $De=0$ ; $E_1=1.0$ , $E_2=0.64$ . (a) $t=2\text{Δ}t$ ; (b) $t=8\text{Δ}t$ ; (c) $t=10\text{Δ}t$ ; (d) $t=11\text{Δ}t$ ; (e) $t=12\text{Δ}t$ ; (f) $t=14\text{Δ}t$ .

Figure 23. Evolution of $\rho$ : Compressible solid medium, $C_2=0.0009$ , $De=0$ ; $E_1=1.0$ , $E_2=0.64$ . (a) $t=2\text{Δ}t$ ; (b) $t=8\text{Δ}t$ ; (c) $t=10\text{Δ}t$ ; (d) $t=11\text{Δ}t$ ; (e) $t=12\text{Δ}t$ ; (f) $t=14\text{Δ}t$ .

Figure 24. Evolution of $\theta$ : Compressible solid medium, $C_2=0.0009$ , $De=0$ ; $E_1=1.0$ , $E_2=0.64$ . (a) $t=2\text{Δ}t$ ; (b) $t=9\text{Δ}t$ ; (c) $t=10\text{Δ}t$ ; (d) $t=11\text{Δ}t$ ; (e) $t=12\text{Δ}t$ ; (f) $t=14\text{Δ}t$ .

Figure 25. Entropy Ψ versus time step number and entropy Ψ versus time $t$ . (a) Entropy Ψ versus time step; (b) Entropy Ψ versus time $t$ .

study, conservation and balance laws of classical mechanics and the constitutive theories are considered in contravariant second Piola-Kirchhoff stress tensor and covariant Green’s strain tensor as measures of stresses and strains. Conservation and balance laws of classical continuum mechanics constituting the mathematical model are first presented in $R^3$ : conservation of mass, balance of linear momenta, balance of angular momenta, energy equation, entropy inequality and the constitutive theories for equilibrium and deviatoric contravariant second Piola-Kirchhoff stress tensors and heat vector and equation of state, thermodynamic pressure. Dissipation mechanism is incorporated using rates of Green’s strain tensor up to order n i.e. using ${\epsilon }_{\left[i\right]};i=1,2,\cdots ,n$ . Rheology mechanism is incorporated using raters of deviatoric contravariant second Piola-Kirchhoff stress tensor up to order m i.e. using ${}_d{\sigma }^{\left[j\right]};j=1,2,\cdots ,m$ . This yields constitutive theory for deviatoric second Piola-Kirchhoff stress tensor in which dissipation is an ordered strain rate mechanism containing a spectrum of dissipation coefficients corresponding to strain rates ${\epsilon }_{\left[j\right]};j=1,2,\cdots ,n$ . The rheology is also an ordered stress rate mechanism yielding a spectrum of relaxation times corresponding to the stress rates ${}_d{\sigma }^{\left[j\right]};j=1,2,\cdots ,m$ . These constitutive theories are of orders n and m in strain and stress tensors. A simplified yet general constitutive theory for deviatoric contravariant second Piola-Kirchhoffstress tensor and heat vector that are linear in the argument of constitutive tensors is also presented. This is followed by constitutive theories of orders $n=1$ and $m=1$ that are commonly used in applications. The dimensionless form of the partial differential equations in $R^3$ is derived containing dimensionless parameters Ec, Re, Br. This mathematical model in $R^3$ is specialized for $R^1$ to study 1D shock physics in elastic and thermoviscoelastic media.

This mathematical model consists of five equations in five dependent variables $u_1,{}_d{\sigma }_{11}^{\left[0\right]},q_1,\theta$ and $v_1$ constitutes an initial value problem (IVP). The solution of the mathematical model is obtained using space-time coupled finite element method in which is constructed using space-time residual functional. This space-time integral form is space-time variationally consistent for linear as well as nonlinear space-time differential operators, hence the computations are unconditionally stable for all choices of discretization lengths in space and time and the dimensionless parameters in the mathematical model. Newton’s linear method with line search is used to obtain solutions of nonlinear algebraic equations resulting from the space-time finite element formulation.

Extensive numerical studies are presented in pure elastic medium, medium with elasticity and dissipation, and solid medium with elasticity dissipation and rheology. Shock formation, shock wave propagation and reflection are illustrated in elastic medium followed by studies to demonstrate influence of dissipation and rheology on shock fronts, shock waves, and entropy generation. In all studies, evolution of ${}_d{\sigma }_{11}^{\left[0\right]},\rho ,\theta$ and entropy production Ψ are presented and described for the evolutions. It is shown that waves of density and temperature also contain shock fronts behind the peaks of the waves similar to the stress waves of ${}_d{\sigma }_{11}^{\left[0\right]}$ .

In the present work, inclusion of energy equation in the mathematical model allows monitoring of entropy productions and associated thermal field. Various features of rate of entropy generation during wave propagation, shock formulations, shock wave propagation, reflection and propagation after reflection are clearly illustrated. This is the unique aspect of the present work. Various sources of entropy generation are presented and described. The entropy production is monitored in the model problem studies to demonstrate complex nature of rentropy generation physics when the medium is thermoviscoelastic with rheology.

Model problem studies for shock waves in the presence of bimaterial interface are also presented to illustrate propagation, transmission, reflection, and interaction of stress, density, and temperature shock waves. Details of the entropy generation are also presented for these studies.

Acknowledgements

First author is grateful to 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 departments are also acknowledged.