Tensile Shock Physics in Compressible Thermoviscoelastic Solids with Rheology

Abstract

The paper addresses tensile shock physics in compressible thermoviscoelastic solids with rheology i.e., in compressible polymeric solids. Polymeric solids have elasticity, dissipation mechanisms and relaxation phenomena due to the presence of long chain molecules in the viscous medium. Thus, the main focus of this investigation is to study how the presence of rheology influences tensile shock physics compared to the tensile shock physics in thermoelastic solids with same elasticity and dissipation but without rheology. A traveling stress wave in polymeric solids leaves nonzero stress signature behind it that naturally influences density. Complete relaxation of the nonzero signatures of stress and associated density change depends upon viscosity of the medium, Deborah number, strength of the stress or velocity wave etc. These aspects of the tensile shock physics in TVES with rheology are investigated in the paper. The mathematical model for finite deformation, finite strain is derived using CBL and CCM, and constitutive theories are derived using conjugate pairs in the entropy inequality and the representation theorem. This mathematical model is thermodynamically and mathematically consistent and has closure. The solution of the IVPs described by this mathematical model related to tensile shock physics in TVES with memory are obtained using space-time coupled finite element method based on space-time residual functional for a space-time strip with time marching. p-version hierarchical space-time local approximations with higher order global differentiability in hpk-scalar product spaces and use of minimally conforming spaces ensure that all space-time integrals over space-time discretization are Riemann. This facilitates more accurate description of the physics in the computational process. Model problem studies are presented to illustrate various aspects of tensile shock physics in compressible TVES with rheology.

Share and Cite:

Surana, K. and Abboud, E. (2024) Tensile Shock Physics in Compressible Thermoviscoelastic Solids with Rheology. Applied Mathematics, 15, 856-886. doi: 10.4236/am.2024.1512050.

1. Introduction, Literature Review, and Scope of Work

1.1. Introduction and Literature Review

In a recent paper, Surana et al. [1] presented shock physics in compressible thermoviscoelastic solid (TVES) without rheology due to compressive stress and velocity stimuli (pulse), using axial deformation of a TVE rod, fixed at one end and subjected to a velocity pulse of duration 2Δt at the other end. Authors demonstrated formation of compressive shockss of stress and density, their propagation, reflection, and interactions. The key element of this shock physics is the continuous change in stress and density over the support of the wave, which results in continuous change in the wave speed along the base of the pulse. Faster moving compression waves behind a wave result in piling up of the compression waves resulting in a compression shock of stress as well as density. On the other hand, faster moving waves ahead of a wave result in shallowing or rarefaction of the wave. Authors showed that persistent oscillations behind a stress and density wave in a TES are due to inability of the medium to dampen the vibrations of the material points. The slightest amount of damping completely eliminates these oscillations but at the expense of amplitude decay and base elongation of the stress and density waves. Damping only influences amplitude and the base of the wave, the wave speed remains unaffected. In another recent paper [2], authors investigated shock physics in compressible thermoviscoelastic solid (TVES) with memory due to compressive stress and velocity pulse. In a completely oscillation-free stress and density evolutions of the waves in TVE solids, the oscillation behind the shock wave reappears due to addition of rheology. Authors showed that these oscillations are due to unrelaxed stresses behind the wave, i.e., due to rheology of the solid medium. In order for the stresses behind a shock to completely relax (i.e., become zero), a finite amount of time is required that depends upon relaxation modulus, which in turn is a function of relaxation time (Deborah number, dimensionless relaxation time). Higher De results in more pronounced oscillations behind the stress or density wave that require longer time to relax (become zero). In references [1] [2], it was shown that under compressive pulse loading the compressive shocks of stress and density always appear behind the peak of the wave and shallowing or rarefaction of the waves always occur ahead of the peak of the wave. This is illustrated clearly and is explained to be due to density change along the base of the pulse that results in changes in wave speed.

In a third recent paper [3], authors investigated tensile shock physics in compressible TVES without memory. Axial deformation of a 1D TVE rod without memory was used to illustrate various aspects of shock formation, propagation, reflection, and interaction due to a tensile velocity pulse loading. In this case, density always remain below its initial value before the commencement of the evolution during evolution of the waves. The largest values of density are at the commencement of the evolution. This phenomenon changes the shock formation physics completely compared to [1] in which the pulse loading is compressive. The authors simulated tensile shock formation, its propagation, reflection and interaction. Stress and density variation along the base of the tensile pulse result in tensile shock ahead of the peak of the pulse. The shallowing of the wave or rarefaction occurs behind the peak of the pulse, this is precisely reverse of what was reported by the authors in reference [1] for compressive loading. Since the stress and density change along the base of the pulse are instrumental in the formation of tensile shocks. It is perhaps beneficial to discuss this using an illustrative tensile stress and density pulse and its propagation. This material that follows has been presented in reference [3], but it is vital to include here for understanding the tensile shock formation process, hence is included here. For example, a propagating tensile stress wave will result in a continuous change in density over the support of the wave (as it does in the case of compressive stress waves [1]), resulting in a continuous change in wave speeds over the support of the wave when the elastic constants remain unchanged. Consider the evolution of deviatoric stress wave ( d σ 11 ) and density ρ in 1D wave propagation of a tensile stress pulse of duration 2Δt in a compressible TVE solid medium. The dimensionless wave speed is unity. At time 2Δt , the tensile stress pulse is completely in the solid medium, as shown in Figure 1(a) (wave ABC). This causes a change in density, as shown in Figure 1(b) (abc). Density decreases from its reference value ρ 0 to its lowest value at b, corresponding to the peak stress at B shown in Figure 1(b), and then increases from its lowest value at b back to the reference value ρ 0 at c (following the change in stress from B to C in Figure 1(a)). Thus, from A to B, the wave speed continuously increases from its reference value (lowest value) at A (or a) to B (or b), its highest value, and then decreases from B (or b) back to C (or c) to its reference value. Therefore, from A to B, the wave speed is continuously increasing, whereas from B to C, it is continuously decreasing. Thus, during further evolution in the portion AB of the wave, as we move from A to B, wave speed is increasing, i.e., points closer to A are moving slower compared to points closer to B, resulting in shallowing of the portion AB of the wave (reduced slope). The opposite happens in the region BC. From B to C, wave speed is decreasing, meaning that the part of the wave closer to B is moving faster than the part of the wave closer to C. This results in piling up of the waves (as in the Riemann shock tube [4] [5]), causing steepening or increased slope of the portion BC of the wave. We observe from Figure 1 that at t>2Δt (red color), shock formation occurs ahead of the peak of the wave (portion BC) and rarefaction or shallowing of the wave occurs behind the peak of the wave. When comparing these results with compressive shock physics studies in ref [1], this phenomenon is exactly opposite of what is observed in [1], where shock formation occurs behind the peak of the wave.

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

Published works on wave propagation and shock physics in compressible solid continua are virtually nonexistent. Stress waves in incompressible solid continua without dissipation have been studied most abundantly, for example by Surana et al. [6] [7]. Achenbach [8] presented wave propagation in elastic solids; in [9], the authors provide an introduction to continuum mechanics and elastic wave propagation; in [10], the authors discuss the propagation of weak shock waves in compressible hyperelastic solids using the theory of singular surfaces. Interfacial wave propagation in initially stressed compressible hyperelastic materials is discussed in [11]. Stress waves in solids, their transmission, reflection, and interaction, as well as fracture caused by them, are presented in [12]. Wave propagation in 2D viscoelastic metamaterials is considered in [13]. Modeling of wave propagation in elastic solids via a higher-order accurate implicit mesh discontinuous Galerkin method is presented in [14]. Wave propagation in functionally graded materials by the modified smooth particle hydrodynamics (MSPH) method is discussed in [15]. While the equation of state is part of the mathematical model in reference [15], the influence of compressibility on wave propagation and density is not demonstrated in the paper. Linear and nonlinear elastic wave propagation is also reported in [7]. We point out that there are many other published works on wave propagation that are not cited here as they are not related to the work presented in this paper.

To our knowledge, the mathematical models based on CBL of CCM with finite strain, finite deformation for compressible TVES with memory and their accurate solution for wave propagation have not been reported in the published literature. This in fact is the motivation for the present work.

1.2. Scope of Work

This paper presents investigations related to the tensile shock physics of stress and density waves in compressible thermoviscoelastic solids with memory (polymeric solids) under tensile loading. Tensile shocks of stress and density are formed due to variations of stress and density along the base of the applied velocity or stress pulse. Investigation of the formation, propagation, reflection, and interaction of stress and density waves with tensile shocks due to applied tensile loading in compressible TVES with rheology is the main objective of the work presented in this paper. The mathematical model consists of the CBL of CCM for finite deformation, finite strain, compressible thermoviscoelastic solids with memory derived using contravariant second Piola-Kirchhoff stress tensor and Green’s strain tensor. The nonlinear dissipation mechanism is incorporated using the rate of Green’s strain tensor up to order n. Rheology or memory mechanism is due to contravariant second Piola-Kirchhoff stress tensor and its time derivatives of order m. In the present work, we consider n = 1 and m = 1 in the model problem (for simplicity). In the mathematical model, additive decomposition of the σ [ 0 ] is considered to address volumetric and distortion deformation physics that are mutually exclusive. Constitutive theory for equilibrium stress tensor e σ [ 0 ] , is derived using contravariant Cauchy stress tensor in Eulerian description and the constitutive theory for contravariant Cauchy stress tensor in tern is derived using entropy inequality in Eulerian description and Helmholtz free energy density. Constitutive theory for d σ [ 0 ] is derived using conjugate pairs in the entropy equality and representation theorem [16]-[27].

The mathematical model derived in this paper is both thermodynamically and mathematically consistent, has closure, and is ideally suited for investigating tensile shock physics in thermoviscoelastic solids (TVES) with rheology. Solutions of the mathematical models are obtained through a space-time coupled finite element method, leveraging the space-time residual functional for a space-time strip with time marching [28]. The integral form derived from this method guarantees unconditionally stable computational processes. The local approximations for the space-time elements are p-version hierarchical in higher order scalar product H Ω ¯ xt ( p,k ) spaces in which p=( p 1 , p 2 ) and k=( k 1 , k 2 ) , in space and time, ensure that the local approximations in space and time maintain the desired higher-order global differentiability with the desired degree of polynomials. In this computational methodology with minimally conforming approximation spaces [28] [29], accurate a posteriori computation of the L 2 norm of the space-time residual functional is possible. The proximity of this norm to zero is a measure of the accuracy of the solution. Generally, computed solutions with this norm of the order of O( 10 6 ) or lower yield sufficiently converged solutions. Only upon obtaining a converged solution for the current space-time strip, the solution is time-marched to the next space-time strip. This approach ensures a converged solution for each space-time strip, hence converged solutions for the entire evolution.

It was discussed in reference [1] that the compressibility in solids is due to the deformation gradient tensor. The density in the current configuration is determined by the conservation of mass. Consequently, thermodynamic pressure in compressible solids arises from density changes due to CM. The inclusion of e σ [ 0 ] , which depends on thermodynamic pressure, is essential in the mathematical model because its absence would result in incorrect force balance in the balance of linear momenta (BLM). While the study of tensile shock physics can be conducted without e σ [ 0 ] , but the results obtained would be in error due to incorrect force balance in the BLM. In the present work, we have employed a simple equation of state that is dependent solely on density. Model problem consisting of an axial rod of TVE material with rheology fixed at left end and subjected to a velocity pulse of time duration 2Δt at the right end is used to simulate and demonstrate various aspects of tensile shock physics.

2. Mathematical Model

For finite deformation, finite strain physics in TVE compressible solid matter with rheology, we consider CBL of CCM derived using the contravariant second Piola-Kirchhoff stress tensor σ [ 0 ] and the covariant Green’s strain tensor ε [ 0 ] in Lagrangian description. The dissipation mechanism is described by the convected time derivatives of the Green’s strain tensor up to order n ( ε [ i ] ;i=1,2,,n ). The relaxation or rheology physics is incorporated by considering convected time derivatives d σ [ i ] ;i=1,2,,m of the deviatoric contravariant second Piola-Kirchhoff stress tensor, each rate results in a characteristic time constant, relaxation time. Thus, in this constitutive theory, we have a spectrum of relaxation times, a more suitable model of the physics of rheology [30]. d σ [ m ] is used as the constitutive tensor with ε [ i ] ;i=0,1,,n , d σ [ j ] ;j=0,1,,m1 and θ as its argument tensors. Conservation of mass (CM), balance of linear momentum (BLM), balance of angular momentum (BAM), and the first and second laws of thermodynamics [31] [32] are given in the following:

ρ 0 ( x )=| J |ρ( x,t ) (1)

ρ 0 2 { u } t 2 ρ 0 { F b }( [ J ][ σ [ 0 ] ] ){ }=0 (2)

ϵ ijk σ ij ( 0 ) =0 (3)

ρ 0 c v De Dt q σ [ 0 ] : ε ˙ [ 0 ] =0 (4)

ρ 0 ( Dϕ Dt η Dθ Dt ) σ [ 0 ] : ε ˙ [ 0 ] + qg θ 0 (5)

where ρ 0 is the density at a material point in the reference (or initial) configuration, ρ( x,t ) is the density at a material point at x i in the current configuration, [ J ] is the deformation gradient tensor, u are displacements in the fixed x-frame, F b is the body force vector per unit mass, σ [ 0 ] is the symmetric contravariant second Piola-Kirchhoff stress tensor, is the gradient operator, ϵ ijk is the permutation tensor, c v is the specific heat, e is specific internal energy density, q is the heat vector, ϕ is the Helmholtz free energy density, η is the entropy density, θ is the absolute temperature, and g is the temperature gradient vector. We consider the additive decomposition of the stress tensor σ [ 0 ] into equilibrium ( e σ [ 0 ] ) and deviatoric ( d σ [ 0 ] ) stress tensors:

σ [ 0 ] = σ e [ 0 ] + σ d [ 0 ] (6)

The constitutive theory for e σ [ 0 ] describes volumetric deformation and the constitutive theory for d σ [ 0 ] describes distortional deformation physics. Following references [31] [32], the constitutive theory for e σ [ 0 ] can be derived using thermodynamic pressure p( ρ,θ ) and contravariant Cauchy stress tensor, and we can write:

[ e σ [ 0 ] ]=| J |p( ρ,θ ) [ [ J ] T [ J ] ] 1 (7)

Thus, we note that e σ [ 0 ] is not a pressure field. The reduced form of the entropy inequality [32] is given by:

σ [ 0 ] : ε ˙ [ 0 ] + qg θ 0 (8)

We note that (8) is satisfied if σ [ 0 ] : ε ˙ [ 0 ] >0 (i.e., the rate of work is positive)

and qg θ 0 . From the rate of work conjugate pair in conjunction with the

axioms of constitutive theory [31] [32], we conclude that the deviatoric contravariant second Piola-Kirchhoff stress tensor is the constitutive variable and the Green’s

strain tensor is its argument tensor. Likewise, the conjugate pair qg θ suggests

that q is constitutive tensor and g and θ its argument tensors. The presence of memory effects necessitates the inclusion of a memory modulus, which dictates that the constitutive theory for d σ [ 0 ] must be represented as a differential equation in time. Consequently, it is essential to consider both d σ [ 0 ] and d σ [ 1 ] , with d σ [ 1 ] serving as a constitutive tensor that incorporates d σ [ 0 ] as one of its argument tensors, among others. We generalize this to include d σ [ j ] ;j=0,1,,m convected time derivatives in place of d σ [ 0 ] as argument tensors. The constitutive theory for the deviatoric contravariant second Piola-Kirchhoff stress tensor is derived by considering d σ [ m ] , convected time derivatives of e σ [ 0 ] of order m, as the constitutive tensor with the following argument tensors:

σ d [ m ] = σ d [ m ] ( ε [ i ] ; σ d [ j ] ;θ );i=0,1,,n;j=0,1,,m1 (9)

In an insulated system, entropy production and the resulting temperature rise are solely attributed to dissipation. This dissipation is minimal for non-cyclic external stimuli, indicating that non-isothermal effects may not need to be considered in the context of deformation physics. However, we present a derivation of the constitutive theories that includes the variable θ . Using (9) in conjunction with the representation theorem, the constitutive theory for d σ [ m ] can be derived using the derivation presented in [31] [32]. The basic steps are outlined in the following.

Let σ G ˜ i ;i=1,2,,N be the combined generators of the argument tensors of d σ [ m ] in (9) that are symmetric tensors of rank two and let σ I j ;j=1,2,,M be the combined invariants of the same argument tensors in (9), all in the current configuration. Then I, G ˜ σ i ;i=1,2,,N constitute the complete basis (integrity) of the space containing constitutive tensor d σ [ m ] . Hence, d σ [ m ] can be expressed as a linear combination of I, G ˜ σ i ;i=1,2,,N in the current configuration

d σ [ m ] = α ˜ σ 0 I+ i=1 N α ˜ σ i ( G ˜ σ i ) (10)

in which

α ˜ σ i = α ˜ σ i ( I ˜ σ j ,θ );j=1,2,,M;i=0,1,,N (11)

We note that in the constitutive theory (10), α ˜ σ i ;i=0,1,,N are not material coefficients. These are simply coefficients in the linear combination in (10) that show dependence on invariants I ˜ σ j ;j=1,2,,M and θ (for non-isothermal case). To determine material coefficients, α ˜ σ i ;i=0,1,,N are expanded in Taylor series in I ˜ σ j ;j=1,2,,M and θ about a known configuration Ω _ . These are then substituted in (10) and terms are rearranged so that each term has quantities in the current configuration multiplied with quantities or coefficients in the known configuration Ω _ . The quantities in the known configuration are the material coefficients. We can obtain the following [31] [32].

Consider constitutive theory (10). We expand each α ˜ σ i ;i=0,1,,N in Taylor series in I ˜ σ j ;j=1,2,,M and θ about a known configuration Ω and retain only up to linear terms in I ˜ σ j ;j=1,2,,M and θ (for simplicity).

α ˜ σ i = α ˜ σ i | Ω _ + j=1 M α ˜ σ i I ˜ σ j | Ω _ ( I ˜ σ j I ˜ σ j | Ω _ ) + α ˜ σ i θ | Ω _ ( θ θ| Ω _ );i=0,1,,N (12)

Substitute α ˜ σ i ;i=0,1,,N from (12) in (10)

d σ [ m ] =( α ˜ σ 0 | Ω _ + j=1 M α ˜ σ 0 I ˜ σ j | Ω _ ( I ˜ σ j I ˜ σ j | Ω ) + α ˜ σ 0 θ | Ω _ ( θ θ| Ω _ ) )I + i=1 N ( α ˜ σ i | Ω _ + j=1 M α ˜ σ i I ˜ σ j | Ω _ ( I ˜ σ j I ˜ σ j | Ω _ ) + α ˜ σ i θ | Ω _ ( θ θ| Ω _ ) ) G ˜ σ i (13)

collecting coefficients (those terms that are defined in known configuration Ω _ ) of I, I ˜ σ j I, G ˜ σ i ,( θ θ| Ω _ ) G ˜ σ i and ( θ θ| Ω _ )I and defining new coefficients we can write the following.

σ ˜ i 0 | Ω _ = α ˜ σ 0 | Ω _ j=1 M α ˜ σ 0 σ I j | Ω _ ( I ˜ σ j | Ω _ ); a ˜ σ j = α ˜ σ 0 σ I j | Ω _ b ˜ σ i = α ˜ σ i | Ω _ j=1 M α ˜ σ i I ˜ σ j | Ω _ ( I ˜ σ j | Ω _ ); c ˜ σ ij = α ˜ σ i I ˜ σ j | Ω _ i=1,2,,N; j=1,2,,M (14)

d σ [ m ] = σ 0 | Ω _ I+ j=1 M a ˜ σ j ( I ˜ σ j )I+ i=1 N j=1 M c ˜ σ ij ( I ˜ σ j ) G ˜ σ i + i=1 N b ˜ σ i ( G ˜ σ i ) (15)

This constitutive theory (15) is based on integrity (complete basis), the representation theorem, and the conjugate pairs in the entropy inequality, hence it is thermodynamically and mathematically consistent constitutive theory. The coefficients a ˜ σ j , c ˜ σ ij , b ˜ σ i are the material coefficients that can be functions of

I ˜ σ j | Ω _ ;j=1,2,,M and θ| Ω _ . This constitutive theory requires ( M+N+MN )

material coefficients.

A simplified linear constitutive theory

We consider an ordered rate constitutive theory of up to orders m and n that is linear in the components of ε [ 0 ] , ε [ i ] ; i=1,2,,n and d σ [ j ] ; j=0,1,,m1 tensors. This constitutive theory will contain generators I , ε [ 0 ] , ε [ i ] ; i=1,2,,n , d σ [ j ] ; j=0,1,,m1 and invariants tr( ε [ 0 ] ) , tr( ε [ i ] ) ; i=1,2,,n , tr( d σ [ j ] ) ; j=0,1,,m1 . Using (15) and only retaining these generators and invariants, we can write (after introducing new notation for material coefficients):

d σ [ m ] = σ ˜ 0 | Ω _ I+ a 1 ε [ 0 ] + a 2 tr( ε [ 0 ] )I+ i=1 n b i 1 ε [i] + i=1 n b i 2 tr( ε [ i ] )I + j=0 m1 c j 1 ( d σ [ j ] )+ j=0 m1 c j 2 tr( d σ [ j ] )I (16)

We can rewrite (16) in the following form using more familiar notation for the material coefficients and neglecting the last term (as done usually in the published works).

d σ [ 0 ] + i=1 m λ i ( d σ [ i ] )= σ ˜ 0 | Ω _ I+ i=0 n 2 η i ( ε [ i ] )+ i=0 n κ i tr( ε [ i ] )I (17)

λ i ,i=1,2,,m , in the spectrum of relaxation times. η i and κ i are the first and second viscosities associated with strain rates ε [ i ] ( i=1,2,,n ), where η 0 and κ 0 are similar to Lamé’s constants in linear elasticity. For heat flux q the following simple heat conduction law is commonly used (see [31] [32] for a more comprehensive constitutive theory for q ).

q=kg (18)

The simplest possible constitutive theory can be obtained from (17) by considering n=1 and m=1 (neglecting σ 0 ):

d σ [ 0 ] + λ 1 ( d σ [ 1 ] )=2 η 0 ( ε [ 0 ] )+ κ 0 tr( ε [ 0 ] )I+2 η 1 ( ε [ 1 ] )+ κ 1 tr( ε [ 1 ] )I (19)

In which:

ε [ 0 ] = 1 2 ( [ J ] T [ J ][ I ] ),[ J ]=[ d J ]+[ I ],[ d J ]=[ { u } { x } ]

with [ d J ] being the displacement gradient tensor and

ε [ 1 ] = D Dt [ ε [ 0 ] ]= t [ ε [ 0 ] ]=( [ J ˙ ] T [ J ]+ [ J ] T [ J ˙ ] )

[ J ˙ ] being the velocity gradient tensor in Lagrangian description.

This mathematical model ((2), (4), (18), and (19)) consists of 13 equations: BLM (3), FLT (1), constitutive theories for d σ [ 0 ] (6) and q (3) in thirteen variables: u (3), θ (1), d σ [ 0 ] (6) and q (3), hence the mathematical model has closure.

This mathematical model can be used to study 3D compressible deformation physics of TVES with memory. To illustrate the shock physics of deviatoric stress and density waves more clearly, it is necessary to consider a model problem in which dispersion and boundary effects causing transmission and reflection do not contaminate the solution. For this reason, we consider 1D stress wave propagation in compressible thermoviscoelastic medium with memory. Furthermore, if we consider an insulated system, then the entropy generation is only due to viscous dissipation, which leads to extremely small temperature changes (for non-cyclic loads) that have virtually no effect on mechanical deformation. Thus, with this assumption, we can eliminate the energy Equation (4) and temperature θ as a dependent variable in the mathematical model.

As pointed earlier, in the Lagrangian description, CM is not part of the mathematical model as ρ( x,t ) in the current configuration is deterministic using CM (1). Thus, for the 1D case, we have displacement u 1 and deviatoric contravariant second Piola-Kirchhoff stress d σ 11 [ 0 ] as the only two dependent variables in the balance of linear momenta in the x 1 direction and the constitutive theory for d σ 11 [ 0 ] defined using (19). Thermodynamic pressure is only a function of ρ and is known when [ J ] is known (Equation (7)).

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

[ J ]=J=| J |=( 1+ u 1 x 1 ) (20)

ε [ 0 ] 11 = u 1 x 1 + 1 2 ( u 1 x 1 ) 2 (21)

ε ˙ [ 0 ] = ε [ 1 ] 11 = [ ε [ 0 ] ] 11 t = v 1 x 1 + u 1 x 1 v 1 x 1 ; d σ 11 [ 1 ] = t ( d σ 11 [ 0 ] ) (22)

where

v 1 x 1 = t ( u 1 x 1 ) (23)

Using (20)-(23) and (6), the mathematical model for 1D wave propagation in a compressible TVE medium with memory can be written as:

ρ 0 2 u 1 t 2 ρ 0 F 1 b x 1 [ ( 1+ u 1 x 1 ) σ e 11 [ 0 ] ] x 1 [ ( 1+ u 1 x 1 ) σ d 11 [ 0 ] ]=0 (24)

e σ 11 [ 0 ] = p( ρ ) 1+ u 1 x 1 (25)

d σ 11 [ 0 ] + λ 1 ( d σ 11 [ 1 ] )= C 1 ( u 1 x 1 + 1 2 ( u 1 x 1 ) 2 )+ C 2 ( v 1 x 1 + u 1 x 1 v 1 x 1 ) (26)

and

v 1 = u 1 t ( x,t ) Ω xt = Ω x × Ω t (27)

in which C 1 =2 η 0 + κ 0 and C 2 =2 η 1 + κ 1 . C 1 and C 2 are material coefficients for elasticity and dissipation. Further details on the equilibrium contravariant second Piola-Kirchhoff stress tensor are given in the following.

In (25), we choose a simple equation of state p( ρ ) given in the following:

p( ρ )=C( ρ ρ 0 1 ) (28)

in which C is the bulk modulus. From (28), we note that increasing ρ implies ρ ρ 0 >1 , yielding positive p, confirming that (28) has correct physics.

Substituting (28) in (25) and then (25) in the third term of BLM (24), we can obtain the following (assuming C to be constant):

x 1 [ ( 1+ u 1 x 1 ) σ e [ 0 ] ]= p( ρ ) x 1 (29)

in which p( ρ ) can be written as follows (using CM):

p( ρ )=C( ρ ρ 0 1 )=C[ ( 1+ u 1 x 1 ) 1 1 ] (30)

We substitute (30) in (29) to obtain

x 1 [ ( 1+ u 1 x 1 ) e σ 11 [ 0 ] ]=C( 1 ( 1+ u 1 x 1 ) 2 2 u 1 x 2 ) (31)

Equation (31) can be used in the balance of linear momentum (Equation (24)) to obtain the final form of the BLM. If we consider thermodynamic pressure to be positive when compressive (28) then we should change the sign of the third term in (24). We do so in (24) and then write (32). The final form of the CBL and the constitutive theory are given in the following for 1D wave propagation in a TVES medium with memory.

ρ 0 v 1 t ρ 0 F 1 b +C( 1 ( 1+ u 1 x 1 ) 2 2 u 1 x 1 2 ) x 1 [ ( 1+ u 1 x 1 ) σ 11 [ 0 ] ]=0 (32)

σ 11 [ 0 ] + λ 1 σ 11 [ 0 ] t = C 1 ( u 1 x 1 + 1 2 ( u 1 x 1 ) 2 )+ C 2 ( v 1 x 1 + u 1 x 1 v 1 x 1 ) (33)

v 1 = u 1 t ( x,t ) Ω xt = Ω x × Ω t (34)

in which C 1 is the elastic coefficient and C 2 is the damping coefficient.

3. Solutions of the PDEs in the Mathematical Model

The system of three partial differential equations in the mathematical Models (32)-(34) are systems of nonlinear partial differential equations (PDEs) that constitute initial value problem (IVP) in which the dependent variables u 1 , d σ 11 [ 0 ] , and v 1 exhibit simultaneous dependence on space coordinates x 1 and time t. We consider a space-time coupled finite element method based on a space-time residual functional for a space-time strip with time marching for obtaining solutions of the PDEs (32)-(34) constituting the mathematical model. It is shown in reference [28] that the space-time differential operators in the IVPs are either non-self-adjoint or nonlinear, for which all other space-time methods of approximation (Space-Time Galerkin Method (STGM), Space-Time Petrov-Galerkin Method (STPGM), Space-Time Weighted Residual Method (STWRM), Space-Time Galerkin Method/Weak Form (STGM/WF)) yield space-time variationally inconsistent integral forms that do not ensure unconditionally stable computations during the evolution. Only the space-time integral form based on space-time residual functional is space-time variationally consistent, hence yields unconditionally stable computational processes during the entire evolution. Basic steps of the computational process based on space-time residual funtional used in this work are given in the following (see reference [28] for more details).

Let t=0 be the time at which the evolution commences. Let Δt be an increment of time. The space-time domain Ω ¯ xt ( 1 ) = Ω ¯ x × Ω ¯ t ( 1 ) =[ 0,L ]×[ 0,Δt ] , of the first space-time strip shown in Figure 2(a) and 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. Hence, they are in higher-order scalar product space H Ω ¯ xt ( p,k ) , p=( p 1 , p 2 ) , k=( k 1 , k 2 ) , where one and two refer to space and time, respectively.

Let ( u 1 ) h e , ( v 1 ) h e , and ( d σ 11 [ 0 ] ) h e be the local approximations of u 1 , v 1 , and d σ 11 [ 0 ] over space-time element Ω ¯ xt e such that

( u 1 ) h = e ( u 1 ) h e , ( v 1 ) h = e ( v 1 ) h e , ( d σ 11 [ 0 ] ) h = e ( d σ 11 [ 0 ] ) h e (35)

Figure 2. Space-time domain; discretization of nth space-time strip.

in which ( u 1 ) h , ( v 1 ) h , and ( d σ 11 [ 0 ] ) h are approximations of u 1 , v 1 , and d σ 11 [ 0 ] over the discretization ( Ω ¯ xt ( 1 ) ) T of Ω ¯ xt ( 1 ) .

( Ω ¯ xt ( 1 ) ) T = e Ω ¯ xt e (36)

The local approximations over a space-time element Ω ¯ xt e can be written as

( u 1 ) h e = i=1 n u N i u ( ξ,η ) { δ i e } u ( v 1 ) h e = i=1 n v N i v ( ξ,η ) { δ i e } v ( d σ 11 [ 0 ] ) h e = i=1 n σ N i σ ( ξ,η ) { δ i e } σ (37)

where N i u , N i v , and N i σ are approximation functions for an element Ω ¯ xt e mapped in the natural coordinate space ( ξ,η ) in a two-unit square, and { δ i e } u , { δ i e } v , and { δ i e } σ are the corresponding degrees of freedom for u, v, and d σ 11 [ 0 ] . Substituting (37) in either (32)-(34) yields three space-time residual equations (functions) E 1 e , E 2 e , and E 3 e x,t Ω ¯ xt e or ξ,η Ω ¯ ξ,η . The space-time finite element method based on the residual functional utilizes residual equations E 1 e , E 2 e , and E 3 e over an element Ω ¯ xt e . Following references [28], we can write the space-time residual functional for an element e with domain Ω ¯ xt e .

I e = i=1 3 ( E i e , E i e ) Ω ¯ xt e (38)

and the space-time residual functional for the discretization ( Ω ¯ xt ( 1 ) ) T of the first space-time strip is given by

I= e I e = e i=1 3 ( E i e , E i e ) Ω ¯ xt e (39)

An extremum of I requires that we set δI=0 (a necessary condition) provided I is differentiable in its arguments, then δI is unique.

δI= e δ I e =2 e i=1 3 ( E i e ,δ E i e ) Ω ¯ xt e =2 e i=1 3 { g e } i =2 e { g } e ={ g }=0 (40)

Let

{ δ }= { δ e } u { δ e } v { δ e } σ (41)

be the total degrees of freedom for ( Ω ¯ xt ( 1 ) ) T . Then { g } in (40) is a nonlinear function of { δ } as the PDEs in the mathematical models are nonlinear. We obtain a solution { δ } that satisfies (40) using Newton’s linear method with line search (see reference [28] for details).

4. Model Problem and Their Solutions

First, we non-dimensionalize the PDEs (32)-(34) in the mathematical model. We write the PDEs (32)-(34) with a hat (^) on all quantities, indicating they all have their usual dimensions or units.

ρ ^ 0 v ^ 1 t ^ ρ ^ 0 F ^ 1 b + C ^ ( 1 ( 1+ u ^ 1 x ^ 1 ) 2 2 u ^ 1 x ^ 1 2 ) x ^ 1 [ ( 1+ u ^ 1 x ^ 1 ) σ ^ d 11 [ 0 ] ]=0 (42)

d σ ^ 11 [ 0 ] + λ ^ 1 ( d σ ^ 11 [ 0 ] ) t ^ = C ^ 1 ( u ^ 1 x ^ 1 + 1 2 ( u ^ 1 x ^ 1 ) 2 )+ C ^ 2 ( v ^ 1 x ^ 1 + u ^ 1 x ^ 1 v ^ 1 x ^ 1 ) (43)

v ^ 1 = u ^ 1 t ^ ( x ^ , t ^ ) Ω x ^ t ^ = Ω x ^ × Ω t ^ (44)

We choose the following reference and dimensionless quantities:

x 1 = x ^ 1 L 0 , ρ 0 = ρ ^ 0 ( ρ 0 ) ref , u 1 = u ^ 1 L 0 , σ 11 [ 0 ] = σ ^ 11 τ 0 , t 0 = L 0 v 0 , τ 0 = p 0 = ( ρ 0 ) ref v 0 2 , v 0 = E 0 ( ρ 0 ) ref , E= C ^ 1 E 0 , η= C ^ 2 η 0 , λ= λ ^ 1 t 0 =De, F 1 b = F ^ 1 b F 0 , F 0 = v 0 2 L 0 , C= C ^ τ 0 = C ^ E 0 } (45)

Using (45) in (42)-(44), we can obtain their following dimensionless forms:

ρ 0 v 1 t ρ 0 F 1 b +C( 1 ( 1+ u 1 x 1 ) 2 2 u 1 x 1 2 ) x 1 [ ( 1+ u 1 x 1 ) σ d 11 [ 0 ] ]=0 (46)

d σ 11 [ 0 ] +De ( d σ 11 [ 0 ] ) t =E( u 1 x 1 + 1 2 ( u 1 x 1 ) 2 )+ η Re ( v 1 x 1 + u 1 x 1 v 1 x 1 ) (47)

v 1 = u 1 t ( x,t ) Ω xt = Ω x × Ω t (48)

where

Re= ρ ref v 0 L 0 η 0 (49)

Re is the Reynolds number and η Re is the dimensionless damping coefficient (we note that η Re will be referred to as C 2 in all the later studies). E

is the dimensionless modulus of elasticity. Deborah number De is the dimensionless relaxation time.

We consider a one-dimensional rod with a constant cross-section along its length, made of compressible thermoviscoelastic (TVE) solid material with memory. The rod has a dimensionless length L=1 and is fixed at the left end ( x 1 =0 ). It is subjected to a compressive stress or a velocity pulse of duration 2Δt at the right end ( x 1 =L ). Figure 3(a) shows the schematic and the loading at x 1 =L . If we assume that all points at a cross-section of the rod displace in the x 1 direction by the same amount, then the problem becomes one-dimensional and we can idealize it by a line from x 1 =0 to x 1 =L (Figure 3(b)). Figure 3(c) shows the first space-time strip Ω ¯ xt ( 1 ) =[ 0,L ]×[ 0,Δt ] .

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

This strip is discretized (uniformly) using nine-node p-version hierarchical space-time finite elements with higher global differentiability. Since the initial value problem contains up to second-order derivatives in space and time, k 1 = k 2 =3 (solution of class C 2 in space and time) is the minimally conforming space for which the space-time integrals over the discretization ( Ω ¯ ( xt ) ( 1 ) ) T are Riemann integrable. When k 1 = k 2 =2 (solutions of class C 1 in space and time), the space-time integrals over ( Ω ¯ ( xt ) ( 1 ) ) T are Lebesgue integrable. The applied disturbance v 1 is a tensile velocity pulse of duration 2Δt at x=1 .

At t=0 , the conditions are as follows: v 1 =0 and v 1 t =0 . At t=Δt , we have v 1 = v 1 * and v 1 x 1 =0 . At t=2Δt , the conditions revert to v 1 =0 and v 1 x 1 =0 . For t2Δt , it remains that v 1 =0 . Similar conditions would hold if we

were to consider a tensile stress pulse of d σ 11 [ 0 ] . This description of velocity v 1 (or stress d σ 11 [ 0 ] ) is of class C 1 in time.

We choose the following material properties and reference quantities:

L 0 =1 ( ρ 0 ) ref =1240 v 0 =116.05 η 0 =372

ρ ^ =1240 E ^ =16.7× 10 6 η ^ =129.51

E 0 = ( ρ 0 ) ref × v 0 2 =16.7× 10 6 t 0 =8.6169× 10 3

De={ 0.0004when λ ^ 1 =3.44677233× 10 6 0.0006when λ ^ 1 =5.17015854× 10 6

5. Model Problem Studies

In this section we consider the following studies.

1) Tensile stress and density waves in linear TVES with rheology

2) Tensile shock front in compressible TVES with memory (non linear TVES with memory)

3) Tensile shock physics in compressible TVES with rheology due to rectangular velocity pulse

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

5.1. Tensile Stress and Density Waves in Linear TVES with Rheology

In linear TVES with rheology, due to small strain and small deformation physics, the material is nearly incompressible, hence density remains approximately constant during the evolution. If the elasticity, i.e., stiffness does not change during deformation, which is the case in linear viscoelasticity, then the wave speed remains the same for TES as well as for TVES with and without memory, regardless of the extent of dissipation and memory. However, the stresses will differ in TES, TVES without memory, and TVES with memory. We expect TVES without memory to have lower stress compared to TES due to viscous dissipation that causes amplitude decay and base elongation, and TVES with memory to have higher stresses compared to TVES without memory due to additional elasticity of the long-chain molecules.

Figure 4(a) and Figure 4(b) show the plots of d σ 11 [ 0 ] versus x 1 at t=7Δt and 18Δt (after reflection from the boundary at x 1 =0 ). We clearly observe larger peak values of d σ 11 [ 0 ] for TVES with memory compared to TVES without memory ( De=0 ). Since the energy imparted to the rod by velocity pulse is fixed, larger peak values of d σ 11 [ 0 ] must be accompanied by a reduction in the base of the wave so that energy content of the wave remains the same as for the imparted velocity pulse. We clearly note progressive peak increase and base reduction from the results in Figure 4(a) and Figure 4(b) as De is increased.

Even though for infinitesimal deformation, the material is assumed incompressible

we know that u 1 x 1 0 , which must be the case to have strain and stress. Thus,

a change in density due to

ρ( x,t )= ρ 0 ( x,t ) 1+ u 1 x 1

is inevitable. Figure 5(a) and Figure 5(b) show plots of ρ versus x 1 at t=7Δt and 18Δt . Density values follow stress d σ 11 [ 0 ] , i.e., an increase in tensile d σ 11 [ 0 ] results in an reduction in ρ . By reducing v * , we can obviously reduce the change in density, suggesting that we are approaching the incompressible case. We clearly observe the same wave speed for TES, TVES with and without memory. As expected, an increased Deborah number results in increased peak stress value but further reduction in the base of the wave. We shall see in the study presented in section 5.2 that finite deformation, finite strain, deformation physics for compressible TVES with rheology is drastically different from what we have observed here for infinitesimal deformation (linear viscoelasticity with memory).

Figure 4. (a): d σ 11 [ 0 ] versus x 1 at t=7Δt ; (b): d σ 11 [ 0 ] versus x 1 at t=15Δt —Linear Case—Thermoviscoelastic with C 2 =0.0009 .

5.2. Tensile Shock Front in Compressible TVES with Memory (Non Linear TVES with Memory)

We apply a tensile velocity pulse of maximum amplitude v * =0.1 and of base 2Δt at the right end of the rod ( x 1 =1.0 ). Figures 6(a)-(d) show plots of d σ 11 [ 0 ] versus x 1 for t=2Δt,5Δt,10Δt , and 15Δt . We have chosen a damping coefficient of 0.0009 and Deborah numbers of 0.0004 and 0.0006 for this study. At t=2Δt , the tensile stress pulse has just entered the rod (Figure 6(a)). Even at this early stage in the evolution, the stress wave in TVES without rheology is moving

Figure 5. (a): ρ versus x 1 at t=7Δt ; (b): ρ versus x 1 at t=15Δt —Linear Case—Thermoviscoelastic with C 2 =0.0009 .

slower than the stress wave in TVES with rheology. In Figure 6(b), we observe the formation of a stress front or shock front ahead the peak of the wave. The shock front is more pronounced in the case of TVES with rheology due to higher stresses (obvious from the peak values in Figure 6(a) and Figure 6(b)). Increased stiffness due to the tensile stress field [1] [33] and reduced density result in higher wave speed in the case of TVES with memory.

In our earlier paper [3] on tensile shock physics, we showed that in the case of compressible TES, there are oscillations in the evolution of d σ 11 [ 0 ] behind the

Figure 6. (a): d σ 11 [ 0 ] versus x 1 at t=2Δt ; (b): d σ 11 [ 0 ] versus x 1 at t=5Δt ; (c): d σ 11 [ 0 ] versus x 1 at t=10Δt ; (d): d σ 11 [ 0 ] versus x 1 at t=15Δt —Nonlinear Case—Thermoviscoelastic with C 2 =0.0009 .

shock front. We explained that their presence is due to vibrations of material points in the neighborhood of the shock front. In the absence of dissipation, there is no mechanism of dampening these except when the shock moves further from these locations, resulting in these material points receiving less energy for excitation. In reference [3], we also showed that the addition of even very mild dissipation completely eliminates these oscillations, as the vibrational energy of these material points is now converted into entropy. From Figures 6(a)-6(d), we note that for De=0 , the d σ 11 [ 0 ] versus x 1 plot has no oscillations behind or ahead of the peak due to dissipation converting vibrational energy into entropy. However, when De0 , we observe the appearance of oscillations ahead of the peak as well as behind it but more pronounced ahead of the peak. These are stresses that have not relaxed yet. We see more intense stress and density oscillations close to the shock front. These diminish as we move away from the shock front, as at these locations stresses have had time to relax some. Furthermore, for De=0.0006 , the magnitude of the oscillation is higher compared to De=0.0004 , as expected, because at De=0.0006 , longer time is required for relaxation compared to De=0.0004 .

Figure 6(c) at t=10Δt shows the shock front at the impermeable boundary ( x 1 =0 ). The shock front for De=0 is moving slower and hence has just reached the boundary x 1 =0 at t=10Δt , whereas for De0 , the shock front has already reflected and we see the reflected waves. Reflected shock fronts propagating towards x 1 =1.0 are shown in Figure 6(d) for t=15Δt . We observe the shock front is always ahead of the wave with oscillation in d σ 11 [ 0 ] in the vicinity of the shock front, with the highest peak immediately adjacent to the shock front. We note that a higher De results in a higher magnitude of oscillations. Lower wave speed of d σ 11 [ 0 ] versus x 1 for De=0 compared to De0 can be observed. Since De=0.0004 and De=0.0006 are fairly close, we don’t observe appreciable increase in wave speed for De=0.0006 compared to De=0.0004 . Figures 7(a)-(d) show the density ρ versus x 1 graphs for t=2Δt,5Δt,10Δt , and 15Δt . These follow exactly the same pattern as the d σ 11 [ 0 ] versus x 1 graphs in Figures 6(a)-(d).

5.3. Tensile Shock Physics in Compressible TVES with Rheology Due to Rectangular Velocity Pulse

In this study, we consider a rectangular tensile velocity pulse of duration 4Δt with a peak value of v * =0.1 , as shown in Figure 8. For 0tΔt and 3Δtt4Δt , v( t ) is a cubic distribution with the values 0, v * and v * , 0, and

zero slopes ( d v 1 dt =0 ) at t=0 , t=Δt , t=3Δt , and t4Δt . For Δtt3Δt ,

v= v * and v( t )=0 for t4Δt . This problem, in similar form, has been used as a model problem in reference [15] for compressive loading. The physics considered in ref [15] is not clear, hence the solutions presented in reference [15] cannot be compared with those presented here, also the loading in the present study is tensile. We elaborate more in the following.

Figure 7. (a): ρ versus x 1 at t=2Δt ; (b): ρ versus x 1 at t=5Δt ; (c): ρ versus x 1 at t=10Δt ; (d): ρ versus x 1 at t=15Δt —Nonlinear Case—Thermoviscoelastic with C 2 =0.0009 .

Figure 8. Rectangular tensile velocity pulse of duration 2Δt .

The discretization for a space-time strip, p-levels in space and time, as well as orders of approximation space in the space and time, are considered the same as used for the first model problem. Evolution is computed by time marching. We use designation “linear” on the graphs to refer to small deformation incompressible physics but with dissipation and memory, and the designation “nonlinear” for finite deformation, finite strain, compressible physics with dissipation and memory. Figures 9(a)-(d) show plots of d σ 11 [ 0 ] versus x 1 for t=4Δt,8Δt,15Δt , and 17Δt . At t=4Δt , the stress wave is completely in the medium. Figure 9(b) shows further evolution of d σ 11 [ 0 ] . Figure 9(c) shows the stress wave after reflection from the impermeable boundary at x 1 =0 . Figure 9(d) show further evolution of stress wave after reflection. Corresponding density ρ versus x 1 graphs are shown in Figures 10(a)-(d). We observe similar behavior as discussed earlier for a velocity pulse of duration 2Δt . We summarize these in the following.

1) Wave speed is higher in nonlinear TVES with memory compared to the linear case due to the tensile stress field resulting in increased stiffness and reduced density.

2) Formation of shock front ahead of the wave before and after reflection.

3) Persistent oscillations in the neighborhood of the shock front due to rheology (unrelaxed stress).

4) The constant v * (hence constant stress d σ 11 [ 0 ] ) region initially persists

Figure 9. (a): d σ 11 [ 0 ] versus x 1 at t=4Δt ; (b): d σ 11 [0] versus x 1 at t=8Δt ; (c): d σ 11 [0] versus x 1 at t=15Δt ; (d): d σ 11 [ 0 ] versus x 1 at t=17Δt C 2 =0.0009 De=0.0004 .

(Figure 9(a) at t=4Δt ) but slowly disappears during evolution resulting in oscillations ahead of the shock front as well. In Figures 9(b)-(d), we see oscillations in d σ 11 [ 0 ] behind the shock front as well as ahead of the shock front.

5) In Figure 9(c) and Figure 9(d), the flat portion of the wave (constant d σ 11 [ 0 ] region) completely disappears.

6) In all reported solutions, the space-time residual functional for the discretization of each space-time strip is O( 10 8 ) or lower, and the use of minimally conforming space ensures that the reported solutions indeed satisfy PDEs with higher degrees of accuracy.

6. Summary and Conclusion

This paper presents investigations of various aspects of the tensile shock physics in compressible TVES with rheology using mathematical model for finite deformation, finite strain deformation physics based on CBL and CCM. The mathematical model considers Green’s strain tensor ε [ 0 ] and its rates ε [ i ] ; i=1,2,,n up to order n for describing ordered rate mechanism of dissipation and deviatoric contravariant second Piola-Kirchhoff stress tensor d σ [ 0 ] and its rates d σ [ j ] ; j=0,1,,m up to order m describing spectrum of relaxation times. The constitutive theories are derived using entropy inequality and the representation theorem. The mathematical model describing IVP is thermodynamically and mathematically consistent and has closure. The solutions of the IVPs resulting from the mathematical model is obtained using space-time coupled finite element methods based on space-time residual functional for a space-time strip with time marching. To the author’s knowledge, the work presented in this paper has not been reported in the published literature. We provide a summary of the research and offer conclusions based on our findings.

1) The mathematical model used is thermodynamically and mathematically consistent.

2) Space-time coupled finite element method based 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 space-time residual functional I0 , the calculated evolution approaches the theoretical solution.

3) We have shown that in case of TVES without memory, there are no oscillations in the vicinity of the shock fronts of stress or density. This is due to dissipation dampening the vibrations of the material points. The oscillations in stress and density reappear for TVES with memory in the vicinity of the shock front due to the presence of unrelaxed stresses. When the shock front moves away from a location, stress and density relax at that location based on memory modulus which is the function of relaxation time. We have shown in the model problem studies that higher relaxation time resulting in higher Deborah numbers has more pronounced peaks of the oscillations in stress and density behind the shock front compared to lower De due to higher values of unrelaxed stress and density.

Figure 10. (a): ρ versus x 1 at t=4Δt ; (b): ρ versus x 1 at t=8Δt ; (c): ρ versus x 1 at t=15Δt ; (d): ρ versus x 1 at t=17Δt C 2 =0.0009 De=0.0004 .

4) Higher stiffness, hence higher elasticity, and lower density due to tensile loading will result in higher wave speed. In the present model problem studies, tensile stress increases stiffness and lowers density rsulting in higher wave speed compared to TVES without memory.

5) Shock front formation due to piling up of higher speed waves is clearly demonstrated in the model problem studies.

6) Propagation of shock fronts, their reflection from the impermeable boundary and reversal of the shock front so that is always ahead of the traveling wave is clearly illustrated in the model problem studies.

7) Study of square wave propagation illustrates interesting aspect of the disappearance of the constant stress and density (corresponding to constant velocity) portion of the wave upon evolution (Figures 9(a)-(d), Figures 10(a)-(d)). We also note that in this study oscillation exists on both sides of the shock (i.e., behind the shock as well as ahead of the shock), Figure 9(c) and Figure 9(d), due to constant velocity or stress region forcing a sharp change.

In conclusion, the work presented in this paper addresses all aspects of tensile shock physics of stress and density waves in compressible TVES with rheology: formation, propagation, reflection, interaction, influence of dissipation, and role of change of stiffness on shock speed. We note that: 1) mathematical models are precisely based on CBL of CCM and the constitutive theories are derived using conjugate pairs in entropy inequality and the representation theorem, hence the complete mathematical model(s) are thermodynamically and mathematically consistent. 2) The space-time coupled finite element method based on the space-time residual functional in higher-order scalar product spaces with minimally conforming spaces ensures that all space-time integrals over discretized space-time domain are Riemann. This improves convergence of the computed solutions to the true solution and permits accurate computations of a posteriori errors based on space-time residual functional. The solutions reported in the paper are nearly as accurate as theoretical solutions, with the space-time residual functional for each space-time discretization consistently achieving O( 10 8 ) or lower. The high accuracy of the computed solutions is achieved by using minimally conforming spaces, which ensure that the physics is correctly represented in the computational process.

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

Elastic Coefficient

C 2

Damping Coefficient

x ¯ , x ¯ i , { x ¯ }

deformed Coordinates

x , x i , { x }

undeformed Coordinates

[ J ]

Deformation gradient tensor

| J |

Determinant of [ J ]

u , u i , { u }

displacements in Lagrangian description

v , v i , { v }

velocities in Lagrangian description

ρ 0

density at time t 0

( ρ 0 ) ref

reference density

ρ

density in Lagrangian description in the current configuration

η

specific entropy in Lagrangian description

η 0

Reference Dissipation

λ

Relaxation time

t 0

Reference time

e

specific internal energy in Lagrangian description

p

thermodynamic or Mechanical Pressure in Lagrangian description

θ

temperature in Lagrangian description

E

Elastic modulus

q , q i , { q }

heat vector in Lagrangian description

g , g i , { g }

temperature gradient tensor in Lagrangian description

σ [ 0 ] , σ 11 [ 0 ]

Second Piola-Kirchhoff stress tensor

e σ [ 0 ] , e σ 11 [ 0 ]

Equilibrium part of the Second Piola-Kirchhoff stress tensor

d σ [ 0 ] , d σ 11 [ 0 ]

Deviatoric part of the Second Piola-Kirchhoff stress tensor

σ ( 0 ) , σ 11 ( 0 )

Cauchy stress tensor

e σ ( 0 ) , e σ 11 ( 0 )

Equilibrium part of the Cauchy stress tensor

d σ ( 0 ) , d σ 11 ( 0 )

Deviatoric part of the Cauchy stress tensor

ε [ 0 ]

Green’s strain tensor (convected time derivative of order zero)

ε [ 1 ]

first convected time derivative of the covariant Green’s strain tensor

ε [ i ]

convected time derivative of order i of the covariant Green’s strain tensor

Abbreviations

CBL

Conservation and Balance Laws

CCM

Classical Continuum Mechanics

TES

Thermoelastic Solid

TVES

Thermoviscoelastic Solid

BVPs

Boundary Value Problems

IVPs

Initial Value Problems

CM

Conservation of Mass

BLM

Balance of Linear Momenta

BAM

Balance of Angular Momenta

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

[1] Surana, K.S. and Abboud, E. (2024) Shock physics in Compressible Thermoelastic and Thermoviscoelastic Solids. Meccanica.
[2] Surana, K.S. and Abboud, E. (2024) Shock Physics in Compressible Thermoviscoelastic Solid Medium with Rheology. Acta Meccanica.
[3] Surana, K.S. and Abboud, E. (2024) Tensile Shock Physics in Compressible Thermoviscoelastic Solid Medium. Applied Mathematics, 15, 719-744.
https://doi.org/10.4236/am.2024.1510042
[4] Surana, K.S., Kendall, J.K. and Carranga, C.H. (2022) Formation, Propagation and Reflection of 1D Normal Shocks in Riemann Shock Tube. Applied Mathematics, 13, 295-323.
https://doi.org/10.4236/am.2022.133022
[5] Surana, K.S., Reddy, K.P.J., Joy, A.D. and Reddy, J.N. (2014) Riemann Shock Tube: 1D Normal Shocks in Air, Simulations and Experiments. International Journal of Computational Fluid Dynamics, 28, 251-271.
https://doi.org/10.1080/10618562.2014.927056
[6] Surana, K.S., Maduri, R.K., TenPas, P.W. and Reddy, J.N. (2006) Elastic Wave Propagation in Laminated Composites Using the Space-Time Least-Squares Formulation in h, p, k Framework. Mechanics of Advanced Materials and Structures, 13, 161-196.
https://doi.org/10.1080/15376490500451809
[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.
https://doi.org/10.4236/ajcm.2015.53032
[8] Achenbach, J.D. (1973) Wave Propagation in Elastic Solids. Elsevier.
[9] Dobróka, M. and Molnár, J.S. (2014) An Introduction to Continuum Mechanics and Elastic Wave Propagation.
https://geofizika.uni-miskolc.hu/Engineering%20physics%20part%20I..pdf
[10] Bechir, H. and Benslimane, A. (2017) On the Propagation of Weak Shock Waves in Compressible Thermohyperelastic Solids. Acta Mechanica, 229, 87-97.
https://doi.org/10.1007/s00707-017-1961-x
[11] Shams, M. and Ejaz, K. (2022) Interfacial Wave Propagation in Initially Stressed Compressible Hyperelastic Materials. Mathematics and Mechanics of Solids, 27, 2510-2531.
https://doi.org/10.1177/10812865221074304
[12] 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.
[13] Wang, Y., Wang, Y. and Laude, V. (2015) Wave Propagation in Two-Dimensional Viscoelastic Metamaterials. Physical Review B, 92, Article ID: 104110.
https://doi.org/10.1103/physrevb.92.104110
[14] 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.
https://doi.org/10.1016/j.cma.2022.114971
[15] Zhang, G.M. and Batra, R.C. (2007) Wave Propagation in Functionally Graded Materials by Modified Smoothed Particle Hydrodynamics (MSPH) Method. Journal of Computational Physics, 222, 374-390.
https://doi.org/10.1016/j.jcp.2006.07.028
[16] Smith, G.F. (1965) On Isotropic Integrity Bases. Archive for Rational Mechanics and Analysis, 18, 282-292.
https://doi.org/10.1007/bf00251667
[17] Smith, G.F. (1970) On a Fundamental Error in Two Papers of C.-C. Wang “on Representations for Isotropic Functions, Parts I and II”. Archive for Rational Mechanics and Analysis, 36, 161-165.
https://doi.org/10.1007/bf00272240
[18] Smith, G.F. (1971) On Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors and Vectors. International Journal of Engineering Science, 9, 899-916.
https://doi.org/10.1016/0020-7225(71)90023-1
[19] Spencer, A.J.M. (1971) Theory of Invariants. In: Eringen, A.C., Ed., Mathematics, Elsevier, 239-353.
https://doi.org/10.1016/b978-0-12-240801-4.50008-x
[20] Spencer, A.J.M. and Rivlin, R.S. (1958) The Theory of Matrix Polynomials and Its Application to the Mechanics of Isotropic Continua. Archive for Rational Mechanics and Analysis, 2, 309-336.
https://doi.org/10.1007/bf00277933
[21] Spencer, A.J.M. and Rivlin, R.S. (1959) Further Results in the Theory of Matrix Polynomials. Archive for Rational Mechanics and Analysis, 4, 214-230.
https://doi.org/10.1007/bf00281388
[22] Wang, C.C. (1969) On Representations for Isotropic Functions, Part I. Archive for Rational Mechanics and Analysis, 33, 249-267.
[23] Wang, C.C. (1969) On Representations for Isotropic Functions, Part II. Archive for Rational Mechanics and Analysis, 33, 268-287.
[24] Wang, C. (1970) A New Representation Theorem for Isotropic Functions, Part I and Part II. Archive for Rational Mechanics and Analysis, 36, 166-223.
[25] Wang, C.C. (1971) Corrigendum to ‘Representations for Isotropic Functions’. Archive for Rational Mechanics and Analysis, 43, 392-395.
[26] 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.
https://doi.org/10.1016/0020-7225(93)90109-8
[27] Zheng, Q.S. (1993) On Transversely Isotropic, Orthotropic and Relatively Isotropic Functions of Symmetric Tensors, Skew-Symmetric Tensors, and Vectors. International Journal of Engineering Science, 31, 1399-1453.
[28] Surana, K.S. and Reddy, J.N. (2018) The Finite Element Method for Initial Value Problems. CRC/Taylor and Francis.
[29] Surana, K.S. and Reddy, J.N. (2016) The Finite Element Method for Boundary Value Problems. CRC/Taylor and Francis.
[30] Bird, R.B., Armstrong, R.C. amd Hassager, O. (1987) Dynamics od Polymeric Liquids, Volumes 1 and 2, Fluid Mechanics.2nd Edition, John Wiley and Sons.
[31] Surana, K.S. (2015) Advanced Mechanics of Continua. CRC/Taylor and Francis.
[32] Surana, K.S. (2022) Classical Continum Mechanics. 2nd Edition, CRC/Taylor and Francis.
[33] Surana, K.S. and Mathi, S.S.C. (2023) Finite Deformation, Finite Strain Nonlinear Dynamics and Dynamic Bifurcation in TVE Solids. Applied Mathematics, 14, 773-838.
https://doi.org/10.4236/am.2023.1412047

Copyright © 2025 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.