Wave Propagation in Compressible Polymeric Solids: Shock Physics

Karan S. Surana, Elie Abboud

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

DOI: 10.4236/ajcm.2025.153017                

Abstract

In the work presented here, we investigate wave physics in compressible polymeric solids. The nonlinear system of partial differential equations is derived using classical continuum mechanics with appropriate strain and stress measures. Constitutive theories are derived using the representation theorem incorporating ordered rate damping and memory. The complete mathematical model is thermodynamically and mathematically consistent. The solutions of the initial value problems (IVPs) described by the mathematical model are obtained using a space-time coupled finite element method. Model problems and their solutions illustrate the nonlinear physics of wave propagation and shock formation. The research presented here on wave propagation in compressible polymeric solids for compressive loading is not available in the published works.

Keywords

Shock Physics, Finite Strain, Finite Deformation, Rheology, Dissipation, Stress and Density Shocks, Finite Element Method, Space-Time Coupled

1. Published Works, Shock Physics in Polymeric Solids

In our recent paper [1] on waves in compressible thermoviscoelastic, we pointed out that there are virtually no published works in this area that utilize mathematical models using classical continuum mechanics. The reported solutions of the IVPs are rarely shown to be accurate and converged. Published works on thermoviscoelastic solids with rheology are even rarer than those for thermoviscoelastic solids without memory. For the sake of completeness, a list of references [2]-[9] cited and discussed in our earlier paper [1] is included here. Interested readers can refer to reference [1] for more details regarding these references. The physics of wave propagation in a compressible solid medium has been demonstrated and explained by the authors in reference [1] for thermoelastic solids in which there is no entropy generation. In the case of polymeric solids, the viscosity and rheology influence the shape of the propagating wave, and some mechanical work is converted into entropy due to dissipation, resulting in a progressively diminishing peak of the wave and progressively elongated base. Figure 1 shows how an applied pulse in 1D wave propagation forms a stress and density shock front behind the peak of the wave and a rarefaction ahead of the peak of the wave.

(a)

(b)

Figure 1. 1D stress wave propagation. (a) Evolution of ${}_d{\sigma }_{11}$ ; (b) Evolution of density $\rho$ .

In the presence of memory, i.e., in TVE solids with rheology, the details discussed above still hold, with the only difference being that due to memory, the stresses do not become zero immediately at the cessation of disturbance. A TVE solid with a larger value of relaxation time would require more time to completely relax (achieve a stress-free state) upon cessation of external stimulus. Thus, this physics will influence stresses in the medium, especially in the vicinity of the sharp stress fronts, hence influencing the evolution. The study of shock fronts, their propagation, their speed of propagation, and reflections in TVE solids with rheology is the subject of investigation in this paper. The shock front physics and the solutions of the models problems calculated here are compared with TVE solids without memory to clearly demonstrate the impact of rheology on shock fronts.

In this paper, wave propagation in compressible polymeric solids is investigated. The mathematical model consists of the CBL of CCM and the constitutive theories derived using the theory of isotropic tensor [10]-[21]. The damping physics are incorporated using strain rates ${\epsilon }_{\left[i\right]},i=1,2,\cdots ,n$ and stress rates ${}_d{\sigma }^{\left[j\right]},j=1,2,\cdots ,m$ . This leads to a spectrum of dissipation coefficients based on the strain rates and a spectrum of relaxation times based on stress rates.

Numerical studies are conducted using a space-time coupled finite element method based on a space-time residual functional yielding space-time variationally consistent integral form. Complete evolution is obtained by using a space-time strip with time marching. For space-time finite elements, the local approximations are p-version hierarchical with higher-order global differentiability in space and time. Computed solutions are considered converged for each space-time strip when the space-time residual functional for each space-time strip is $\mathcal{O}\left({10}^{-8}\right)$ or lower, ensuring accurate evolution.

We remark that the shock physics presented in this paper for compressible solids is quite different from tensile shock physics considered in reference [22].

2. Mathematical Model

Using ${\sigma }^{\left[0\right]}$ and ${\epsilon }_{\left[0\right]}$ measures, the conservation and the balance laws in Lagrangian description can be written as [23] [24]:

${\rho }_0\left(x,t\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]}^{\text{T}}\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)

Consider additive decomposition of ${\sigma }^{\left[0\right]}$

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

Following references [23] [24]:

$\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)

The entropy inequality [24] becomes:

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

$p\left(\rho ,\theta \right)$ is thermodynamic pressure. We consider dissipation to be dependent on ${\epsilon }_{\left[i\right]},i=1,2,\cdots ,n$ and memory to be due to ${}_d{\sigma }^{\left[j\right]};j=0,1,\cdots ,m$ Hence, we can write the following for the constitutive tensor ${}_d{\sigma }^{\left[m\right]}$ :

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

Let ${}^{\sigma }{\underset{˜}{\mathcal{G}}}^i;i=1,2,\cdots ,N$ be the combined generators of the argument tensors of ${}_d{\sigma }^{\left[m\right]}$ in (9) that are symmetric tensors of rank two and let $I^j;j=1,2,\cdots ,M$ be the combined invariants of the same argument tensors in (9), all in the current configuration. Then ${}_d{\sigma }^{\left[m\right]}$ can be expressed as a linear combination of $I,{}^{\sigma }\underset{˜}{\mathcal{G}}{}^i;i=1,2,\cdots ,N$ in the current configuration

${}_d{\sigma }^{\left[m\right]}={}^{\sigma }\underset{˜}{\alpha }{}^{0}I+{\displaystyle \sum }_{i=1}^N{}^{\sigma }\underset{˜}{\alpha }{}^i\left({}^{\sigma }\underset{˜}{\mathcal{G}}{}^i\right)$ (10)

in which

${}^{\sigma }\underset{˜}{\alpha }{}^i={}^{\sigma }\underset{˜}{\alpha }{}^i\left({}^{\sigma }\underset{˜}{I}{}^j,\theta \right);j=1,2,\cdots ,M;i=0,1,\cdots ,N$ (11)

The material coefficients are obtained using the standard approach of Taylor series expansion of ${}^{\sigma }\underset{˜}{\alpha }{}^i;i=0,1,\cdots ,N$

${}^{\sigma }\underset{˜}{\alpha }{}^i={{}^{\sigma }\underset{˜}{\alpha }{}^i|}_{\underset{\_}{\Omega }}+{\displaystyle \sum }_{j=1}^M{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\left( {}^{\sigma }\underset{˜}{I}{}^j-{{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\right)+{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial \theta }|}_{\underset{\_}{\Omega }}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right);i=0,1,\cdots ,N$ (12)

Substitute ${}^{\sigma }\underset{˜}{\alpha }{}^i;i=0,1,\cdots ,N$ from (10) in (12)

$\begin{array}{c}{}_d{\sigma }^{\left[m\right]}=\left( {{}^{\sigma }\underset{˜}{\alpha }{}^0|}_{\underset{\_}{\Omega }}+{\displaystyle \sum }_{j=1}^M{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^0}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\left( {}^{\sigma }\underset{˜}{I}{}^j-{{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\right)+{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^0}{\partial \theta }|}_{\underset{\_}{\Omega }}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)\right)I\\ +{\displaystyle \sum }_{i=1}^M\left({{}^{\sigma }\underset{˜}{\alpha }{}^i|}_{\underset{\_}{\Omega }}+{\displaystyle \sum }_{j=1}^M{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\left( {}^{\sigma }\underset{˜}{I}{}^j-{{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\right)+{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial \theta }|}_{\underset{\_}{\Omega }}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)\right){}^{\sigma }\underset{˜}{\mathcal{G}}{}^i\end{array}$ (13)

Defining the following material coefficients and collecting coefficients of $I,{}^{\sigma }\underset{˜}{I}{}^{j}I,{}^{\sigma }\underset{˜}{\mathcal{G}}{}^i,\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right){}^{\sigma }\underset{˜}{\mathcal{G}}{}^i$ and $\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)I$ and defining new coefficients we can write

$\begin{array}{l}{\sigma }_0={{}^{\sigma }{\underset{˜}{\alpha }}^0|}_{\underset{\_}{\Omega }}-{\displaystyle \sum _{j=1}^M}{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^0}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\left({{}^{\sigma }{\underset{˜}{I}}^j|}_{\underset{\_}{\Omega }}\right);{}^{\sigma }{\underset{˜}{a}}_j={\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^0}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\\ {}^{\sigma }{\underset{˜}{b}}_i={{}^{\sigma }{\underset{˜}{\alpha }}^i|}_{\underset{\_}{\Omega }}+{\displaystyle \sum _{j=1}^{\text{}M}}{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\left({{}^{\sigma }{\underset{˜}{I}}^j|}_{\underset{\_}{\Omega }}\right);{}^{\sigma }{\underset{˜}{c}}_{ij}={\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial {}^{\sigma }\underset{˜}{I}{}^j}|}_{\underset{\_}{\Omega }}\\ {}^{\sigma }\underset{˜}{d}{}_i=-{\frac{\partial {}^{\sigma }\underset{˜}{\alpha }{}^i}{\partial \theta }|}_{\underset{\_}{\Omega }};{}^{\sigma }{\alpha }_{tm}=-{\frac{{\partial }^{\sigma }{\underset{˜}{\alpha }}^0}{\partial \theta }|}_{\underset{\_}{\Omega }}\\ i=1,2,\cdots ,N;j=1,2,\cdots ,M\end{array}$ (14)

$\begin{array}{c}{}_d{\sigma }^{\left[m\right]}={\sigma }_{0}I+{\displaystyle \sum }_{j=1}^M{}^{\sigma }\underset{˜}{a}{}_j\left({}^{\sigma }\underset{˜}{I}{}^j\right)I+{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^M{}^{\sigma }\underset{˜}{c}{}_{ij}\left({}^{\sigma }\underset{˜}{I}{}^j\right){}^{\sigma }\underset{˜}{\mathcal{G}}{}^i\\ +{\displaystyle \sum }_{i=1}^N{}^{\sigma }\underset{˜}{b}{}_i\left({}^{\sigma }\underset{˜}{\mathcal{G}}{}^i\right)-{\displaystyle \sum }_{i=1}^N{}^{\sigma }\underset{˜}{d}{}_i\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right){}^{\sigma }\underset{˜}{\mathcal{G}}{}^i-{}^{\sigma }\underset{˜}{\alpha }{}_{tm}\left(\theta -{\theta |}_{\underset{\_}{\Omega }}\right)I\end{array}$ (15)

The coefficients ${}^{\sigma }\underset{˜}{a}{}_j,{}^{\sigma }\underset{˜}{c}{}_{ij},{}^{\sigma }\underset{˜}{b}{}_i,{}^{\sigma }\underset{˜}{d}{}_i$ , and ${}^{\sigma }{\underset{˜}{\alpha }}_{tm}$ are the material coefficients. This constitutive theory requires $\left(M+N+MN\right)$ material coefficients.

A simplified linear constitutive theory

Neglecting the thermal term, a constitutive theory linear in the components of the argument tensors is given by:

${}_d{\sigma }^{\left[0\right]}+{\displaystyle \sum }_{i=1}^m{\lambda }_i\left({}_d{\sigma }^{\left[i\right]}\right)={\sigma }_{0}I+{\displaystyle \sum }_{i=0}^{n}2{\eta }_i\left({\epsilon }_{\left[i\right]}\right)+{\displaystyle \sum }_{i=0}^n{\kappa }_i\text{tr}\left({\epsilon }_{\left[i\right]}\right)I$ (16)

For $n=m=1$ and neglecting ${\sigma }_0$ , we can write (16) as follows

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

Fourier heat conduction law, a linear form of the constitutive theory for $q$ derived using integrity is given by:

$q=-\kappa \cdot g$ (18)

For 1D finite deformation, finite strain, we have:

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

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

${\dot{\epsilon }}_{\left[0\right]}={\epsilon }_{\left[1\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};{}_d{\sigma }_{11}^{\left[1\right]}=\frac{\partial }{\partial t}\left({}_d{\sigma }_{11}^{\left[0\right]}\right)$(21)

where

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

Using (19)-(22), the mathematical model for 1D wave propagation in a compressible polymeric solid can be written as:

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

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

${}_d{\sigma }_{11}^{\left[0\right]}+{\lambda }_1\left({}_d{\sigma }_{11}^{\left[1\right]}\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)$ (25)

and

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

in which $C_1=2{\eta }_0+{\kappa }_0$ and $C_2=2{\eta }_1+{\kappa }_1$ .

We consider the following equation of state:

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

Substituting (27) in (24) and then (24) in the third term of BLM (23), 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 p\left(\rho \right)}{\partial x_1}$ (28)

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]$ (29)

We substitute (29) in (28) 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)$ (30)

Using (30) in the balance of linear momentum we obtain the final mathematical model.

${\rho }_0\frac{\partial v_1}{\partial t}-{\rho }_0{F}_1^b+C\left(-\frac{1}{{\left(1+\frac{\partial u_1}{\partial x_1}\right)}^2}\frac{{\partial }^2{u}_1}{\partial x_1^2}\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$ (31)

${}_d{\sigma }_{11}^{\left[0\right]}+{\lambda }_1\frac{\partial \left({}_d{\sigma }_{11}^{\left[0\right]}\right)}{\partial t}=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)$ (32)

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

in which $C_1$ is elastic coefficient and $C_2$ is damping coefficient.

3. Solutions of the PDEs in the Mathematical Model

Equations (31)-(33) are nonlinear PDEs in dependent variables $u_1$ , ${}_d{\sigma }_{11}^{\left[0\right]}$ , and $v_1$ . We obtain solutions of these using a space-time coupled finite element method based on a space-time residual functional [22] [25] that yields a space-time variationally consistent space-time integral form.

We consider discretization for a space-time strip and time march upon obtaining a converged solution [22] [25] to obtain entire evolutions. 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 Figure 2(a), Figure 2(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. See references [22] [25] for the space-time finite element method based on space-time residual functional for a space-time strip with time marching.

(a)

(b)

(c)

(d)

Figure 2. 1D solid domain, idealization of 1D solid domain, discretization of first space-time strip with space-time finite elements and applied disturbance. (a) polymeric solid domain; (b) Mathematical idealization of solid domain; (c) Disretization of first space-time strip using 10 space-time p-version finite elements; (d) Velocity (v₁) pulse of duration 2Δt.

4. Numerical Studies

To non-dimensionalize PDEs (31)-(33), we write them in the following form (with a hat (^) on each quantity) in which all quantities have their usual dimensions.

${\widehat{\rho }}_0\frac{\partial {\widehat{v}}_1}{\partial \widehat{t}}-{\widehat{\rho }}_0{\widehat{F}}_1^b+\widehat{C}\left(-\frac{1}{{\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)-\frac{\partial }{\partial {\widehat{x}}_1}\left[{\left(1+\frac{\partial {\widehat{u}}_1}{\partial {\widehat{x}}_1}\right)}_d{\widehat{\sigma }}_{11}^{\left[0\right]}\right]=0$ (34)

${}_d{\widehat{\sigma }}_{11}^{\left[0\right]}+{\widehat{\lambda }}_1\frac{\partial \left({}_d{\widehat{\sigma }}_{11}^{\left[0\right]}\right)}{\partial \widehat{t}}={\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)$ (35)

${\widehat{v}}_1=\frac{\partial {\widehat{u}}_1}{\partial \widehat{t}}\forall \left(\widehat{x},\widehat{t}\right)\in {\Omega }_{\widehat{x}\widehat{t}}={\Omega }_{\widehat{x}}\times {\Omega }_{\widehat{t}}$(36)

We choose the following reference and dimensionless quantities:

$\begin{array}{llll}{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 & \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 & \hfill \\ v_0=\sqrt{\frac{E_0}{{\left({\rho }_0\right)}_{\text{ref}}}},\hfill & E=\frac{{\widehat{C}}_1}{E_0},\hfill & \eta =\frac{{\widehat{C}}_2}{{\eta }_0},\hfill & \lambda =\frac{{\widehat{\lambda }}_1}{t_0}=De,\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 & \hfill \end{array}\}$(37)

Using (37) in (34)-(36), we can obtain their following dimensionless forms:

${\rho }_0\frac{\partial v_1}{\partial t}-{\rho }_0{F}_1^b+C\left(-\frac{1}{{\left(1+\frac{\partial u_1}{\partial x_1}\right)}^2}\frac{{\partial }^2{u}_1}{\partial x_1^2}\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$ (38)

${}_d{\sigma }_{11}^{\left[0\right]}+De\frac{\partial \left({}_d{\sigma }_{11}^{\left[0\right]}\right)}{\partial t}=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)$ (39)

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

where

$Re=\frac{{\rho }_{\text{ref}}{v}_0{L}_0}{{\eta }_0}$ (41)

Figures 2(a)-(d) show a schematic of the rod, its mathematical idealization, space-time finite element discretization of the first space-time strip, and the applied disturbance (pulse) at $x=L$ . We note that $k_1=3$ and $k_2=3$ define a minimally conforming scalar product space for which integrals are Riemann.

Material properties and reference quantities used are listed below:

$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$ , $De=\{\begin{array}{l}\text{0}\text{.0004}\text{when}{\widehat{\lambda }}_1=3.44677233\times {10}^{-6}\\ \text{0}\text{.0006}\text{when}{\widehat{\lambda }}_1=5.17015854\times {10}^{-6}\end{array}$

We consider the following numerical studies:

4.1. Wave Physics with Dissipation but without Memory

4.2. Wave Physics with Dissipation and Rheology

4.3. Wave Physics with Dissipation and Memory for a Rectangular Pulse

In all numerical studies, we consider a 30-element uniform mesh for the space-time strip, with p-level of 7 in space and time and $\text{Δ}t=0.1$ , unless specified otherwise. Newton’s linear method is considered converged using a tolerance of $O\left({10}^{-6}\right)$ . The space-time residual functional $\le O\left({10}^{-8}\right)$ is considered sufficient for good accuracy of the computed solution.

4.1. Wave Physics with Dissipation but without Rheology

We note that dissipation (resistance to motion) and rheology (resistance to the motion of long chain molecules), both add stresses in addition to elastic stresses. Figure 3(a), Figure 3(b) show the plots of ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=6\text{Δ}t$ and $15\text{Δ}t$ . We clearly observe larger peak values of ${}_d{\sigma }_{11}^{\left[0\right]}$ for TVE solids with memory. We clearly note progressive peak increase and base reduction from the results for TE solids, TVE, and TVE solids with rheology in Figure 3(a) and

Figure 3(b). Even for linear case (nearly incompressible matter), $\frac{\partial u_1}{\partial x_1}\ne 0$

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

resulting in density change that progressively increases with increasing Deborah number, shown in Figure 4(a) and Figure 4(b) for $t=6\text{Δ}t$ and $15\text{Δ}t$ . Increasing peak values are accompanied by progressively reducing support of the wave for progressively increasing $De$ .

(a) (b)

Figure 3. (a) ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=6\text{Δ}t$ -LinearCase-Thermoviscoelastic with $C_2=0.0009$ ; (b) ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=15\text{Δ}t$ -LinearCase-Thermoviscoelastic with $C_2=0.0009$ .

(a) (b)

Figure 4. (a) $\rho$ versus $x_1$ at $t=6\text{Δ}t$ -LinearCase-Thermoviscoelastic with $C_2=0.0009$ ; (b) $\rho$ versus $x_1$ at $t=15\text{Δ}t$ -LinearCase-Thermoviscoelastic with $C_2=0.0009$ .

4.2. Wave Physics with Dissipation and Rheology

In this study, we consider a velocity pulse of amplitude $v^{\text{*}}=0.1$ , with a damping coefficient of 0.0009, and Deborah numbers $De=0.0004$ and $De=0.0006$ . We monitor ${}_d{\sigma }_{11}^{\left[0\right]}$ at $t=2\text{Δ}t,5\text{Δ}t,10\text{Δ}t$ , and $15\text{Δ}t$ along the length of the rod (Figures 5(a)-(d)).

In Figure 5(b), a mild shock front is observed, which steepens upon further propagation ($t=10\text{Δ}t$ , $t=15\text{Δ}t$ ). At $t=10\text{Δ}t$ , the waves are at the fixed boundary. The reflected waves with shock fronts are shown in Figure 5(d).

Increasing Deborah number results in steeper shock fronts due to increased stress, higher peaks, and reduced support. In reference [7], we showed that oscillations present behind and ahead of the front in purely elastic cases are completely eliminated by the presence of damping. In Figure 5(a) and Figure 5(d), we observe the reappearance of pre- and post-wave oscillations due to the presence of unrelaxed stress field. These oscillations increase in magnitude and spread over longer domains with increasing Deborah number.

Corresponding graphs of the density evolution are shown in Figures 6(a)-(d), with exactly the same behavior as observed in Figures 5(a)-(d) containing ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ graphs.

(a) (b)

(c) (d)

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

4.3. Wave Physics with Dissipation and Rheology for a Rectangular Velocity Pulse

We consider a rectangular velocity pulse (Figure 7) in this study. Plots of ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ are shown in Figures 8(a)-(d). Corresponding density $\rho$ versus $x_1$ graphs are shown in Figures 9(a)-(d). We observed similar behavior as discussed earlier for a velocity pulse of duration $2\text{Δ}t$ . We summarize key observations in the following:

(a) (b)

(c) (d)

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

Figure 7. Rectangular velocity pulse of duration $2\text{Δ}t$ .

(a) (b)

(c) (d)

Figure 8. (a) ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=4\text{Δ}t$ -$C_2=0.0009$ -$De=0.0004$ ; (b) ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=8\text{Δ}t$ -$C_2=0.0009$ -$De=0.0004$ ; (c) ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=15\text{Δ}t$ -$C_2=0.0009$ -$De=0.0004$ ; (d) ${}_d{\sigma }_{11}^{\left[0\right]}$ versus $x_1$ at $t=17\text{Δ}t$ -$C_2=0.0009$ -$De=0.0004$ .

1) Lower wave speed for compressible (nonlinear) case compared to linear case is observed due to increased density.

2) After reflection, we observe distinct shock front formation behind the wave.

3) Due to rheology, the unrelaxed stress fields result in oscillations behind the shock front.

4) The constant ${}_d{\sigma }_{11}^{\left[0\right]}$ region corresponding to constant $v^{\text{*}}$ (as seen in Figure 8(a) at $t=4\text{Δ}t$ ) slowly disappears upon further evolution, resulting in oscillations ahead of the shock front as well.

5) We note in Figure 8(c) and Figure 8(d) that the constant ${}_d{\sigma }_{11}^{\left[0\right]}$ region of the rectangular pulse completely disappeared.

6) The integrated sum of squares of the space-time residual function is $\mathcal{O}\left({10}^{-8}\right)$ or lower for minimally conforming spaces ensures that the PDEs are satisfied accurately across the discretized space-time domain.

(a) (b)

(c) (d)

Figure 9. (a) $\rho$ versus $x_1$ at $t=4\text{Δ}t$ - $C_2=0.0009$ -$De=0.0004$ ; (b) $\rho$ versus $x_1$ at $t=8\text{Δ}t$ -$C_2=0.0009$ -$De=0.0004$ ; (c): $\rho$ versus $x_1$ at $t=15\text{Δ}t$ -$C_2=0.0009$ -$De=0.0004$ ; (d): $\rho$ versus $x_1$ at $t=17\text{Δ}t$ -$C_2=0.0009$ -$De=0.0004$ .

5. Summary and Conclusions

We have presented wave propagation and shock formation in compressible polymeric solids with dissipation and rheology. The following is a brief summary of the work and some conclusions drawn from it.

1) The mathematical model is derived strictly using the CBL of CCM.

2) The constitutive theories are initiated using conjugate pairs in the entropy inequality and are derived using the representation theorem.

3) The complete mathematical model is thermodynamically and mathematically consistent.

4) Space-time coupled finite element method based on space-time residual functional yields unconditionally stable computations. Use of higher-order scalar product spaces in space and time is essential to maintain Riemann integrals over space-time discretization and when the integrated sum of squares of space-time residual $I\to 0$ , the calculated evolution approaches the theoretical solution.

5) A unique aspect of volumetric deformation physics in solids is explained. It is shown that volumetric deformation in solids is due to $\left|J\right|$ and not the equation of state. Volumetric deformation due to $\left[J\right]$ has associated with thermodynamic pressure that defines the equilibrium stress tensor necessary for correct force balance in BLM.

6) The oscillation free solution for stress and density behind and ahead of the wave for the TVE case regains oscillation due to the presence of unrelaxed stresses. With time, their magnitude reduces. Higher De leads to higher stresses, resulting in higher peaks and more pronounced oscillations in the vicinity of the front.

7) Due to compressive loading, higher density results in lower wave speed compared to TE solids or TVE solids without rheology.

8) Shock front formation behind the peak of the wave due to piling up of compression waves moving at higher speeds is clearly observed. Rarefaction ahead of the front is also observed.

9) Reflection of the waves with shock fronts from the clamped boundary and reversal of the shock front so that it is always behind the peak of the waves is clearly observed.

10) This work sheds light on intricate aspects of shock physics in the presence of dissipation and memory that are helpful in the design of polymer-based applications.

11) The applications of this work to actual practical problems of interest in design will require extending this work to 2D and 3D.

12) The novel aspect of this work is shock physics in polymeric solids due to compressive loading compared to reference [22] in which tensile shock physics is investigated.

13) In the case of a rectangular wave, reduction in the constant stress and density region during evolution, presence of oscillations on both sides of the shock, and finally almost complete disappearance of the constant stress and density region followed by pronounced oscillations (note-worthy aspect of the study) are clearly observed.

14) All reported solutions are time-accurate as the computed solution satisfies PDEs in the pointwise sense over the discretized space-time domain due to the use of minimally conforming spaces.

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.

List of Symbols

$C$ :	Bulk Modulus
$C_1$ :	Modulus of Elasticity
$C_2$ :	Dissipation Coefficient
$\overline{x}$ , ${\overline{x}}_i$ , $\left\{\overline{x}\right\}$ :	Coordinates in the current configuration
$x$ , $x_i$ , $\left\{x\right\}$ :	Coordinates in the reference configuration
$\left[J\right]$ :	Deformation gradient tensor
$\left|J\right|$ :	Determinant of $\left[J\right]$
$u$ , $u_i$ , $\left\{u\right\}$ :	Displacements (Lagrangian description)
$v$ , $v_i$ , $\left\{v\right\}$ :	Velocities (Lagrangian description)
${\rho }_0$ :	Density in reference configuration
${\left({\rho }_0\right)}_{\text{ref}}$ :	Reference density
$\rho$ :	Density (current configuration)
$\eta$ :	Entropy density
${\eta }_0$ :	Reference damping
$\lambda$ :	Relaxation time
$t_0$ :	Reference time
$e$ :	Internal energy density
$p$ :	Pressure
$\theta$ :	Temperature
$E$ :	Young’s modulus
$q$ , $q_i$ , $\left\{q\right\}$ :	Heat vector in Lagrangian description
$g$ , $g_i$ , $\left\{g\right\}$ :	Temperature gradients
${\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\}$ :	Cauchy stress (Equilibrium)
${}_d{\sigma }^{\left(0\right)}$ , $\left\{{}_d{\sigma }^{\left(0\right)}\right\}$ :	Cauchy stress (Deviatoric)
${\epsilon }_{\left[0\right]}$ , $\left\{{\epsilon }_{\left[0\right]}\right\}$ :	Green’s strain
${\epsilon }_{\left[i\right]}$ , $\left\{{\epsilon }_{\left[i\right]}\right\}$ :	Material derivative of ${\epsilon }_{\left[0\right]}$ of order $i$

Conflicts of Interest

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

References

[1]	Surana, K.S. and Abboud, E. (2025) Shock Physics in Compressible Thermoelastic and Thermoviscoelastic Solids. Meccanica, 60, 755-783.[CrossRef]
[2]	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.
[3]	Dobróka, M. and Molnár, J. (2014) An Introduction to Continuum Mechanics and Elastic Wave Propagation. University of Miskolc. https://geofizika.uni-miskolc.hu/Engineering%20physics%20part%20I..pdf
[4]	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.[CrossRef]
[5]	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 ID: 114971.[CrossRef]
[6]	Bechir, H. and Benslimane, A. (2017) On the Propagation of Weak Shock Waves in Compressible Thermohyperelastic Solids. Acta Mechanica, 229, 87-97.[CrossRef]
[7]	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.[CrossRef]
[8]	Shams, M. and Ejaz, K. (2022) Interfacial Wave Propagation in Initially Stressed Compressible Hyperelastic Materials. Mathematics and Mechanics of Solids, 27, 2510-2531.[CrossRef]
[9]	Wang, Y., Wang, Y. and Laude, V. (2015) Wave Propagation in Two-Dimensional Viscoelastic Metamaterials. Physical Review B, 92, Article ID: 104110.[CrossRef]
[10]	Smith, G.F. (1965) On Isotropic Integrity Bases. Archive for Rational Mechanics and Analysis, 18, 282-292.[CrossRef]
[11]	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.[CrossRef]
[12]	Smith, G.F. (1971) On Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 9, 899-916.[CrossRef]
[13]	Spencer, A.J.M. (1971) Theory of Invariants. In: Eringen, A.C., Ed., Mathematics, Elsevier, 239-353.[CrossRef]
[14]	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.[CrossRef]
[15]	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.[CrossRef]
[16]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part I. Archive for Rational Mechanics and Analysis, 33, 249-267.[CrossRef]
[17]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part II. Archive for Rational Mechanics and Analysis, 33, 268-287.[CrossRef]
[18]	Wang, 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.[CrossRef]
[19]	Wang, C. (1971) Corrigendum to My Recent Papers on “Representations for Isotropic Functions”. Archive for Rational Mechanics and Analysis, 43, 392-395.[CrossRef]
[20]	Zheng, Q. (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.[CrossRef]
[21]	Zheng, Q. (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.[CrossRef]
[22]	Surana, K.S. and Abboud, E. (2024) Tensile Shock Physics in Compressible Thermoviscoelastic Solids with Rheology. Applied Mathematics, 15, 856-886.[CrossRef]
[23]	Surana, K.S. (2015) Advanced Mechanics of Continua. CRC/Taylor and Francis.
[24]	Surana, K.S. (2022) Classical Continuum Mechanics. 2nd Edition, CRC/Taylor and Francis.
[25]	Surana, K.S. and Reddy, J.N. (2018) The Finite Element Method for Initial Value Problems. CRC/Taylor and Francis.