Finite Deformation, Finite Strain Nonlinear Dynamics and Dynamic Bifurcation in TVE Solids with Rheology

Karan S. Surana, Sri Sai Charan Mathi

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

DOI: 10.4236/am.2024.151009                

Abstract

This paper presents a mathematical model consisting of conservation and balance laws (CBL) of classical continuum mechanics (CCM) and ordered rate constitutive theories in Lagrangian description derived using entropy inequality and the representation theorem for thermoviscoelastic solids (TVES) with rheology. The CBL and the constitutive theories take into account finite deformation and finite strain deformation physics and are based on contravariant deviatoric second Piola-Kirchhoff stress tensor and its work conjugate covariant Green’s strain tensor and their material derivatives of up to order m and n respectively. All published works on nonlinear dynamics of TVES with rheology are mostly based on phenomenological mathematical models. In rare instances, some aspects of CBL are used but are incorrectly altered to obtain mass, stiffness and damping matrices using space-time decoupled approaches. In the work presented in this paper, we show that this is not possible using CBL of CCM for TVES with rheology. Thus, the mathematical models used currently in the published works are not the correct description of the physics of nonlinear dynamics of TVES with rheology. The mathematical model used in the present work is strictly based on the CBL of CCM and is thermodynamically and mathematically consistent and the space-time coupled finite element methodology used in this work is unconditionally stable and provides solutions with desired accuracy and is ideally suited for nonlinear dynamics of TVES with memory. The work in this paper is the first presentation of a mathematical model strictly based on CBL of CCM and the solution of the mathematical model is obtained using unconditionally stable space-time coupled computational methodology that provides control over the errors in the evolution. Both space-time coupled and space-time decoupled finite element formulations are considered for obtaining solutions of the IVPs described by the mathematical model and are presented in the paper. Factors or the physics influencing dynamic response and dynamic bifurcation for TVES with rheology are identified and are also demonstrated through model problem studies. A simple model problem consisting of a rod (1D) of TVES material with memory fixed at one end and subjected to harmonic excitation at the other end is considered to study nonlinear dynamics of TVES with rheology, frequency response as well as dynamic bifurcation phenomenon.

Keywords

Thermoviscoelasticity, Rheology, Memory, Finite Strain, Finite Deformation, Nonlinear Dynamics, Dynamic Bifurcation, Ordered Rate Theories

1. Introduction, Literature Review and Scope of Work

A literature review of published works on finite deformation, and finite strain nonlinear dynamics of TVES with or without memory reveals that there is no unified approach for deriving mathematical models for TVES with memory. Energy methods, purely phenomenological 1D approaches, and principles of virtual work are commonly employed. In most cases, the differential forms of the mathematical models are never derived from the integral forms (generally Euler’s equations corresponding to the functional). Instead, the integral forms are directly used in obtaining finite element or series solutions. Mathematical models based on CBL of CCM and consistent constitutive theory have not been used to our knowledge for nonlinear dynamics studies of TVES with or without memory. The purpose of including the following two sections is to present concepts and methodologies for IVPs and BVPs that are supported by the theory of differential operators and the calculus of variations. This material is critical in evaluating which published formulations, methodologies, etc are supported by the calculus of variations and the theory of differential operators.

1.1. Initial Value Problems (IVPs)

If we employ the CBL of CCM, we can derive a mathematical model [1] [2] that consists of conservation of mass, balance of linear momenta, balance of angular momenta, first and second law of thermodynamics and associated constitutive theories. These can be expressed as:

$A\varphi \left(x\mathrm{,}t\right)-f\left(x\mathrm{,}t\right)=0\forall \left(x\mathrm{,}t\right)\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t$ (1)

in which $A$ is a space-time differential operator, $\varphi \left(x\mathrm{,}t\right)$ is a vector of dependent variables exhibiting simultaneous dependence on space $x$ and time t. It can be shown that $A$ can be linear or nonlinear but cannot be symmetric over ${\bar{\Omega }}_{xt}$ . Thus, $A$ is a non-self adjoint or nonlinear but not self adjoint (linear and symmetric). Thus, space-time energy functional cannot be constructed for (1) using the fundamental lemma of the calculus of variations (space-time Galerkin method with weak form). Hence, for IVPs in solid mechanics, the use of energy methods in which one constructs an integral statement of energy functional directly without knowing (1) is invalid for constructing finite element processes. Also, if one does construct some space-time integral form, the Euler’s equation from the integral will not yield operator $A$ .

Surana et al. [3] have shown that a space-time integral statement of the principle of virtual work is indeed the same as the space-time integral statement resulting from the use of fundamental lemma for (1) i.e.,

${\left(A\varphi \left(x\mathrm{,}t\right)-f\left(x\mathrm{,}t\right)\mathrm{,}v\left(x\mathrm{,}t\right)\right)}_{{\bar{\Omega }}_{xt}}=0$ (2)

This of course requires knowledge of (1) which we do not have if we are constructing the space-time integral statement of virtual work directly (i.e., without the knowledge of (1)). Thus, we see that energy methods and the principle of virtual work have no role in the derivation of the mathematical models for IVPs or in the methods of approximation used for obtaining the solutions of IVPs.

Surana et al. [4] have shown that when considering space-time coupled methods of obtaining the solution of (1), only the space-time finite element method based on space-time residual functional is space-time variationally consistent integral form, hence unconditionally stable computational processes for the entire evolution.

We remark that space-time decoupled methods of approximation [4] can also be used for obtaining the solution of IVPs related to beams, plates and shells or any other physical system, but the methods require a differential form of the mathematical model (1), thus precluding use of energy methods and the principle of virtual work, only the mathematical model based on CBL of CCM remain a viable alternative.

1.2. Boundary Value Problems (BVPs)

The mathematical models for TVE solids with rheology naturally lead to initial value problems. Nonetheless, the literature review reveals that in many instances, the mathematical description of BVPs is established first generally using energy methods or the principle of virtual work and then terms added to these models phenomenologically (generally force balance) to arrive at the mathematical descriptions of the corresponding initial value problems. It is for this reason that the discussion of the differential operators in BVPs, possible integral form strategies and when energy methods are supported by calculus of variations and the mathematical classification of differential operators in BVPs is included in this section.

It is well known [4] [5] that valid BVP must correspond to the stationary state of the corresponding IVP, nonetheless study of BVP is helpful in the study of stationary processes. We consider BVPs in solid and structural mechanics. It can be shown [5] that when the mathematical models of BVPs are derived using CBL of CCM (or NCCM), the differential operators appearing in all BVPs regardless of their origin or complexity can be mathematically classified into three groups: self adjoint (linear and symmetric), non-self adjoint (linear but not symmetric) and nonlinear (neither linear nor symmetric). The mathematical description of infinitesimal, linear elastic and reversible deformation of solid media results in self adjoint operators. The inclusion of dissipation and memory mechanisms with infinitesimal deformation results in a mathematical model that contains non-self adjoint differential operators. Purely elastic deformation or deformation with dissipation and/or memory mechanisms even in the presence of finite deformation and finite strain result in mathematical models in which the differential operators are nonlinear. This classification is essential in demonstrating which methods of approximation are suitable choices for which class of operators.

Using calculus of variations [5] , we can easily show that if (consider a scalar equation for simplicity)

$A\varphi \left(x\mathrm{,}t\right)-f\left(x\right)=0\forall x\in {\Omega }_x$ (3)

is the BVP, then, when A is self adjoint, the intergal form (using test function v),

${\left(A\varphi -f\mathrm{,}v\right)}_{{\bar{\Omega }}_x}=0$ (4)

resulting from fundamental lemma can be shown to yield (using GM/WF) the weak form of (4),

${\left(A\varphi -f\mathrm{,}v\right)}_{{\bar{\Omega }}_e}=B\left(\varphi \mathrm{,}v\right)-l\left(v\right)=0$ (5)

in which, $B\left(\varphi \mathrm{,}v\right)$ is bilinear and symmetric. In this case, the quadratic functional $I\left(\varphi \right)$ can be constructed as

$I\left(\varphi \right)=\frac{1}{2}B\left(\varphi \mathrm{,}\varphi \right)-l\left(\varphi \right)$ (6)

The quadratic functional $I\left(\varphi \right)$ represents the total potential energy, $\frac{1}{2}B\left(\varphi \mathrm{,}\varphi \right)$ is the stored strain energy and $l\left(\varphi \right)$ is the potential energy of loads. In this case, one can show that the $\delta I\left(\varphi \right)=0$ yields (5), thus the Euler’s equations resulting from $\delta I\left(\varphi \right)=0$ is in fact (3), the BVP. Thus, for self adjoint differential operators, energy methods are valid. In this case, if the energy functional ( $I\left(\varphi \right)$ ) is constructed without the knowledge of the differential form of the mathematical model, we must show that the Euler’s equation resulting from the integral statement is indeed differential form of the BVP obtained using CBL of the CCM.

It is needless to discuss the principle of virtual work as this in fact is the integral form resulting from the use of fundamental lemma when the differential form of the mathematical model is known. That is, the virtual work statement cannot be constructed correctly without the knowledge of the differential operator A. When the differential operator is not self adjoint, the quadratic or energy functional $I\left(\varphi \right)$ does not exist, hence the use of energy methods for the BVPs described by non-self adjoint and nonlinear differential operator is not possible and will undoubtedly result in spurious integral forms and spurious finite element solutions based on these integral forms. As in the case of IVPs, here also, when the differential form of the mathematical model is derived using CBL of CCM (or NCCM), the integral form for this mathematical model based on residual functional yields unconditionally stable computational processes for all three types of differential operators. Clearly, energy methods and the principle of virtual work play absolutely no role in this.

1.3. Introduction and Literature Review

The study of the dynamic response (evolution, frequency response, etc) of TES matter experiencing finite deformation and finite strain is nonlinear dynamics. In this case, the stiffness is nonlinear (up to the quadratic function of the dofs) when space and time are decoupled using GM/WF in space. In the case of TVES with finite deformation, finite strain including dissipation physics, a consistent derivation of the constitutive theory for deviatoric second Piola-Kirchhoff stress tensor results in a dissipation mechanism dependent on Green’s strain rates. When decoupling space and time in this case and using GM/WF in space, we obtain stiffness as well as damping matrices, both up to quadratic functions of the dofs. Thus, in this case, we have nonlinear stiffness (as in TES) but we also have nonlinear damping. In TVES with rheology, the composition of the solid matter (rubber like materials) contains long chain molecules which introduce rheology physics in the solid matter. In such solid continua, upon cessation of applied disturbance, the stress-free state is not reached instantaneously, but requires additional time which depends upon relaxation time, a time constant (a material property) of the TVES. This mechanism is called relaxation or stress relaxation or rheology or short-term memory mechanism. In TVES with rheology, the constitutive theory for the deviatoric second Piola-Kirchhoff stress tensor is not a simple algebraic expression involving strains, strain rates and their traces. Instead, it is a differential equation in time containing time derivative(s) of the deviatoric second Piola-Kirchhoff stress tensor as well as the second Piola-Kirchhoff stress tensor itself. In this case, it is not possible to substitute deviatoric second Piola-Kirchhoff stress in BLM to obtain BLM purely in terms of displacements, as is the case for TES and TVES without rheology. Thus, the method of obtaining solutions of IVPs in TVES with rheology generally differs compared to those for TES and TVES without rheology. In this paper, we address details of the mathematical model as well as the method of obtaining the solution of the associated IVPs describing finite deformation, and finite strain nonlinear dynamics of TVES with rheology.

In the following, we discuss published works on conservation and balance laws and constitutive theories for finite deformation, and finite strain physics of thermoviscoelastic solids with memory. The initial value problems described by these mathematical models naturally contain nonlinear elasticity, nonlinear dissipation mechanism and nonlinear mechanism of rheology. The phenomenological constitutive models for polymeric liquids have long been in existence (see reference [6] ). These mathematical models provide useful mathematical descriptions of physics observed in experiments. To our knowledge, the first formal presentation of complete mathematical models and ordered rate constitutive theories for thermoviscoelastic solids with finite deformation physics and memory seem to have appeared in reference [7] . The authors presented the derivation of ordered rate constitutive theory based on entropy inequality and representation theorem. Due to the deviatoric contravariant second Piola-Kirchhoff stress and rate of Green strain tensor as conjugate rate of work pairs, nonlinear elasticity, nonlinear dissipation mechanism and rheology are inherent but are not highlighted in the paper. In a companion paper, Surana et al. [8] presented mathematical models and ordered rate constitutive theories for finite deformation thermoviscoelastic solids with memory using Gibbs potential. Details of the mathematical models and ordered rate constitutive theories for the thermoviscoelastic solid for finite deformation and finite strain can also be found in the textbook by Surana [1] [2] . The IVPs described by the work in references [9] [10] [11] [12] naturally have nonlinear elasticity, nonlinear dissipation mechanism and rheology. We discuss recently published works on TVEs with memory undergoing finite deformation and finite strain. First, we make a few remarks that are helpful when we discuss published works.

1) Energy methods can only be used for BVPs in which the differential operator is self adjoint [5] . This limits the use of energy methods to linear elasticity with reversible deformation physics (i.e., without dissipation and memory).

2) The phenomenological constitutive models are lumped models that have no concept of spatial coordinates, hence cannot possibly describe the constitutive behavior of a continuum in which space and time are intrinsically present.

3) Phenomenological models cannot be extended for continuous matter in ${ℝ}^1\mathrm{,}{ℝ}^2$ or ${ℝ}^3$ .

4) Material coefficients in phenomenological constitutive models are generally not physical, hence almost always show a lack of agreement with those in the constitutive theories derived using entropy inequality and representation theorem.

5) Energy methods cannot be used for IVPs [4] as the space-time differential operators are not self adjoint.

Before we present a review of the published works on the mathematical models for nonlinear dynamics of TVES with memory, we consider the following mathematical models and details that are helpful.

The mathematical models reported and used in the finite deformation, finite strain nonlinear dynamic studies of the TVES without rheology in the published work are primarily phenomenological. These mathematical models are not based on the CBL of CCM and the constitutive theories are not derived using entropy inequality and representation theorem. In many cases, the mathematical models are 1D lumped models in space, hence are purely in terms of a single variable and time. In general, the reported phenomenological models contain nonlinear stiffness and nonlinear damping mechanisms. The situation is much more complex in the case of TVES with memory. In the following, we employ the CBL of CCM and consistent constitutive theories to derive a phenomenological model, aiming to investigate its feasibility.

For simplicity, we consider isothermal physics. For finite deformation, finite strain physics, the BLM and the constitutive theory for deviatoric second Piola-Kirchhoff stress tensor ${}_d\sigma {}^{\left[0\right]}$ in terms of Green’s strain tensor and its rates $\left({\epsilon }_{\left[0\right]}\mathrm{,}{\epsilon }_{\left[1\right]}\right)$ constitute the complete mathematical models based on CBL of CCM.

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

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

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

$\forall \left(x\mathrm{,}t\right)\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t$

${}_e\sigma {}^{\left[0\right]}$ is equilibrium stress, generally referred to as equation of state. By substituting (8) and (9) in (7) and considering space-time decoupled finite element formulation for a spatial discretization ${\bar{\Omega }}_x^T$ of ${\bar{\Omega }}_x$ and using GM/WF in space, we can obtain the following system of nonlinear ODEs in time [13] [14] :

$\left[M\right]\left\{\ddot{\delta }\right\}+\left[C\right]\left\{\stackrel{\dot{}}{\delta }\right\}+\left[K\right]\left\{\delta \right\}=\left\{F\right\}+\left\{P\right\}$ (10)

in which $\left[M\right]$ is a symmetric, constant coefficient mass matrix. $\left[C\right]$ and $\left[K\right]$ are damping and stiffness matrices. Following references [13] [14]

$\left[C\right]=\left[C^1\right]+\left[C^2\right]+\left[C^3\right]$ (11)

$\left[K\right]=\left[K^1\right]+\left[K^2\right]+\left[K^3\right]$ (12)

Matrices $\left[C^1\right]$ and $\left[K^1\right]$ are symmetric and their coefficients are constant, $\left[C^2\right]$ and $\left[K^2\right]$ are nonsymmetric and their coefficients are linear functions of $\left\{\delta \right\}$ , $\left[C^3\right]$ and $\left[K^3\right]$ are also symmetric but their coefficients are quadratic functions of $\left\{\delta \right\}$ . Equations (10) are a system of coupled nonlinear ODEs in time in $\left\{\stackrel{\dot{}}{\delta }\right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ . The mathematical model (7) - (9) and (10) are both based on CCM and are equivalent, except that in (10), space and time are decoupled and integration in space for ${\bar{\Omega }}_x^{\text{T}}$ has eliminated presence of spatial coordinates $x_i$ . But the spatial details present in (7) - (9) are preserved in (10) due to $\left\{\stackrel{\dot{}}{\delta }\right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ corresponding to spatial locations.

A one dimensional phenomenological model in a single dependent variable cannot describe physics in (7) - (9) or in (10), this is rather obvious. However, if we insist in doing so, then in (10) we can retain only a single dependent variable $\delta$ (making others zero) to obtain:

$m\ddot{\delta }+\left(c_1+c_2\delta +c_3{\delta }^2\right)\stackrel{\dot{}}{\delta }+\left(k_1+k_2\delta +k_3{\delta }^2\right)\delta =\underset{˜}{f}+\underset{˜}{P}$ (13)

It is obvious a single ODE in time as in (13) cannot possibly describe the same physics as the system of ODEs (9), hence also can describe the nonlinear dynamics physics in (7) - (9). This mathematical model (10) has been discussed in ref [14] and is in widespread use in the nonlinear dynamics studies in many published works. In this illustration we note that in order to realize the concepts of mass $\left[M\right]$ , damping $\left[C\right]$ and stiffness $\left[K\right]$ , we must be able to substitute constitutive theory (9) in the balance of linear moment (7) and then use GM/WF for the spatial discretization.

Next, we consider CBL of CCM and constitutive theory for TVES with rheology for isothermal physics. BLM and the constitutive theory (CCM) are given by

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

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

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

In this case $\left[{}_d{\sigma }^{\left[0\right]}\right]$ cannot be substituted from (16) into (14), hence the nonlinear ODEs like (10) cannot be derived for TVES with memory. Thus, 1D phenomenological model like (13) is not possible for TVES with memory. A major consequence of this discussion is that in TVES with memory, explicit form of the equations containing mass, damping and stiffness matrices (in the form (10)) is not possible, primarily due to the fact that the constitutive theory for ${}_d\sigma {}^{\left[0\right]}$ is a differential equation in ${}_d\sigma {}^{\left[0\right]}$ and time t, hence ${}_d\sigma {}^{\left[0\right]}$ cannot be substituted in BLM from the constitutive equations.

Thus, for TVES with memory, we obviously need to consider different considerations and approaches in establishing 1D phenomenological models if we wish to do so. Unfortunately, in the published works on mathematical models for TVES with memory, finite deformation and finite strain physics, we only find use of phenomenological models that are not based on CBL of CCM (or NCCM). We discuss some of these in the following. In reference [9] , Amabili considers 1D phenomenological constitutive theories that have been borrowed from the references cited as [15] [16] [17] (1968, 1982, 2009) in reference [9] . From reference [9] , Equations (1), use of two strain rate terms has no basis (based on CCM). This constitutive theory (1D), if considered, should read (only in time):

$\tau \left(t\right)+{\lambda }_1\frac{\text{d}\tau \left(t\right)}{\text{d}t}=E\epsilon \left(t\right)+\eta \stackrel{\dot{}}{\epsilon }\left(t\right)$ (17)

${\lambda }_1$ is relaxation time, E is modulus of elasticity and $\eta$ is viscosity. $\tau$ is stress, $\epsilon$ and $\stackrel{\dot{}}{\epsilon }$ are strain and strain rate.

First, $\tau$ is deviatoric second Piola Kirchhoff stress and not simple Cauchy stress. Secondly, we cannot premultiply (17) by say “Area” (unit area) and view (17) as force balance equation as done in reference [9] :

$F\left(t\right)+\lambda \frac{\text{d}F\left(t\right)}{\text{d}t}=F_s\left(t\right)+F_d\left(t\right)$ (18)

where, $F_s\left(t\right)$ is the force due to strain or nonlinear springs and $F_d\left(t\right)$ is similarly force due to nonlinear dissipation.

If we assume that the force in the spring is up to a cubic function of $\delta$ , then based on reference [9] we have:

$E\epsilon =F_s=k_1\delta \left(t\right)+k_2{\delta }^2\left(t\right)+k_3{\delta }^3\left(t\right)$ (19)

$\stackrel{\dot{}}{\epsilon }=\frac{1}{E}\left(k_1+2k_2\delta +3k_3{\delta }^2\right)\stackrel{\dot{}}{\delta }$ (20)

$F_d=\eta \stackrel{\dot{}}{\epsilon }=\left(\frac{\eta }{E}{k}_1+\frac{2k_2\eta }{E}\delta +\frac{3k_3\eta }{E}{\delta }^2\right)\stackrel{\dot{}}{\delta }=\left(c_1+c_2\delta \left(t\right)+c_3{\delta }^2\left(t\right)\right)\stackrel{\dot{}}{\delta }$ (21)

then, (18) can be written as:

$F\left(t\right)+{\lambda }_1\frac{\text{d}F\left(t\right)}{\text{d}t}=\left(k_1+k_2\delta \left(t\right)+k_3{\delta }^2\left(t\right)\right)\delta \left(t\right)+\left(c_1+c_2\delta \left(t\right)+c_3{\delta }^2\left(t\right)\right)\stackrel{\dot{}}{\delta }$ (22)

In this derivation: 1) there is no basis for (18) from (17). 2) no basis for (19) as strains are spatial derivative of displacements. Nonetheless, if we assume that (18) and (19) hold, then (22) is the result of correct derivation using them. This is similar to Equation (18) in ref [9] (barring third term on the right hand side of (2)). 3) We recall that force balance in continuum mechanics (CM) is due to balance of linear momenta, in which gradients of stresses appear (per unit volume). Thus, the statement of force F(t) in (22) is erroneous. Authors in reference [9] state that if m is the mass and $\ddot{\delta }$ is the acceleration, then $F\left(t\right)$ acting on this mass is in equilibrium with external excitation and they write (Equation (3) in ref [9] ) the following as equation of motion of the system.

$m\ddot{\delta }+F\left(t\right)=\hat{F}\sin \left(\omega t+\varphi \right)$ (23)

Equation (23) is erroneous. We present details in the following using 1D case. We have seen in case of TVES without memory that substitution of ${}_d{\sigma }_{11}^{\left[0\right]}$ in BLM yield BLM purely in displacement $u_1$ . Upon decoupling of space and time and using GM/WF in space, we obtain a system of nonlinear ODEs (10) in which mass, stiffness and damping matrices are identified. Thus, in order to obtain equation with $m\ddot{\delta }$ as in (17), we must consider BLM (1D for simplicity), using $\tau$ for deviatoric second Piola-Kirchhoff stress, we can write:

${\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)\tau \right)=0\forall x\mathrm{,}t\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t$ (24)

If we consider space-time decoupled finite element formulation for discretization ${\bar{\Omega }}_x^T$ of ${\bar{\Omega }}_x$ and use GM/WF in space, then we can obtain the following system of ODEs:

$\left[M\right]\left\{\ddot{\delta }\right\}+\left\{F^{\mathrm{<mo>\; *\; </mo>}}\right\}=\left\{P\right\}+\left\{F\right\}$ (25)

For the sake of completeness and clarity, we present derivation of (25) using (24). Consider a discretization ${\bar{\Omega }}_x^T={\displaystyle \underset{e}{\cup }{\bar{\Omega }}_x^e}$ of spatial domain ${\bar{\Omega }}_x$ in which ${\bar{\Omega }}_x^e$ could be a three node p-version hierarchical element with higher order global differentiability ( ${\bar{\Omega }}_x^e=\left[x_e\mathrm{,}{x}_{e+1}\right]$ ). Let $u_h^e$ and ${\tau }_h^e$ be local approximations of $u_1$ and $\tau$ given by (unequal degree):

$\begin{array}{l}{u}_h^e={\displaystyle \underset{i=1}{\overset{n_u}{\sum }}}{N}_i^u\left(x_1\right){\delta }_i^u\left(t\right)=\left[N^u\right]\left\{{\delta }_e^u\right\}\\ {\tau }_h^e={\displaystyle \underset{i=1}{\overset{n_{\tau }}{\sum }}}{N}_i^{\tau }\left(x_1\right){\delta }_i^{\tau }\left(t\right)=\left[N^{\tau }\right]\left\{{\delta }_e^{\tau }\right\}\end{array}$ (26)

in which $N_i^u\left(x\right)$ and $N_i^{\tau }\left(x_1\right)$ are approximation functions for $u_h^e$ and ${\tau }_h^e$ . ${\delta }_i^u$ and ${\delta }_i^{\tau }$ are corresponding degrees of freedom.

Let $v^u\left(x\right)$ be test functions such that $v^u\left(x_1\right)=\delta u_h^e$ (GM/WF), thus $v^u\left(x_1\right)=N_j^u\left(x_1\right);j=1,2,\cdots ,n_u$ . Using fundamental lemma [4] [5] , we can write the following integral form of (24) over ${\bar{\Omega }}_x^T$ .

${\displaystyle \underset{e}{\sum }}{\left({\rho }_0\frac{{\partial }^2{u}_h^e}{\partial t^2}-{\rho }_0{F}_1^b-\frac{\partial }{\partial x_1}\left(\left(1+\frac{\partial u_h^e}{\partial x_1}\right){\tau }_h^e,v^u\left(x_1\right)\right)\right)}_{{\bar{\Omega }}_x^e}=0$ (27)

or

$\begin{array}{l}{\displaystyle \underset{e}{\sum }}{\left({\rho }_0\frac{{\partial }^2{u}_h^e}{\partial t^2}\mathrm{,}{v}^u\left(x_1\right)\right)}_{{\bar{\Omega }}_x^e}-{\displaystyle \underset{e}{\sum }}{\left({\rho }_0{F}_1^b\mathrm{,}{v}^u\left(x_1\right)\right)}_{{\bar{\Omega }}_x^e}\\ -{\displaystyle \underset{e}{\sum }}{\left(\frac{\partial }{\partial x_1}\left(\left(1+\frac{\partial u_h^e}{\partial x_1}\right){\tau }_h^e\mathrm{,}{v}^u\left(x_1\right)\right)\right)}_{{\bar{\Omega }}_x^e}=0\end{array}$ (28)

In the third term, we perform integration by parts once to transfer one order of differentiation with respect to $x_1$ to $v^u\left(x_1\right)$ .

$\begin{array}{l}{\displaystyle \underset{e}{\sum }}{\left({\rho }_0\frac{{\partial }^2\text{}{u}_h^e}{\partial t{\text{}}^2},v^u\left(x_1\right)\right)}_{{\bar{\Omega }}_x^e}-{\displaystyle \underset{e}{\sum }}{\left({\rho }_0{F}_1^b,v^u\left(x_1\right)\right)}_{{\bar{\Omega }}_x^e}\\ +{\displaystyle \underset{e}{\sum }}{\left(\left(1+\frac{\partial u_h^e}{\partial x_1}\right){\tau }_h^e,\frac{\partial v^u\left(x_1\right)}{\partial x_1}\right)}_{{\bar{\Omega }}_x^e}-{\displaystyle \underset{e}{\sum }}{\left[\left(1+\frac{\partial u_h^e}{\partial x_1}\right){\tau }_h^e{v}^u\left(x_1\right)\right]}_{x_e}^{x_{e+1}}=0\end{array}$ (29)

substituting from (26) into (29) and using $v^u\left(x_1\right)=N_j^u\left(x_1\right);j=1,2,\cdots ,n_u$ , we obtain:

$\begin{array}{l}{\displaystyle \underset{e}{\sum }}{\left({\rho }_0{\displaystyle \underset{i=1}{\overset{n_u}{\sum }}}{N}_i^u\left(x_1\right){\ddot{\delta }}_i^u,N_j^u\left(x_1\right)\right)}_{{\bar{\Omega }}_x^e}-{\displaystyle \underset{e}{\sum }}{\left({\rho }_0{F}_1^b,N_j^u\left(x_1\right)\right)}_{{\bar{\Omega }}_x^e}\\ +{\displaystyle \underset{e}{\sum }}{\left(\left(1+\frac{\partial u_h^e}{\partial x_1}\right){\displaystyle \underset{i=1}{\overset{n_{\tau }}{\sum }}}{N}_i^{\tau }\left(x_1\right){\delta }_i^{\tau },\frac{\partial N_j^u\left(x_1\right)}{\partial x_1}\right)}_{{\bar{\Omega }}_x^e}-{\displaystyle \underset{e}{\sum }}{\left[\left(1+\frac{\partial u_h^e}{\partial x_1}{\tau }_h^e\right)N_j^u\left(x_1\right)\right]}_{x^e}^{x^{e+1}}=0\end{array}$ (30)

or

${\displaystyle \underset{e}{\sum }}\left(\left[M^e\right]\left\{{\ddot{\delta }}_e^u\right\}-\left\{F^e\right\}+\left[B^e\right]\left\{{\delta }_e^{\tau }\right\}+\left\{P^e\right\}\right)=0$ (31)

$\left\{{\ddot{\delta }}^u\right\}-\left\{F\right\}+\left[B\right]\left\{{\delta }^{\tau }\right\}-\left\{P\right\}=0$ (32)

$\begin{array}{l}\left[M\right]={\displaystyle \underset{e}{\sum }}{m}^e\mathrm{,}\left\{F\right\}={\displaystyle \underset{e}{\sum }}\left\{F^e\right\}\mathrm{,}\left[B\right]={\displaystyle \underset{e}{\sum }}\left[B^e\right]\mathrm{,}\left\{P\right\}={\displaystyle \underset{e}{\sum }}\left\{P^e\right\}\\ \left\{{\ddot{\delta }}^u\right\}={\displaystyle \underset{e}{\cup }\left\{{\ddot{\delta }}_e^u\right\}}\text{and}\left\{{\delta }^{\tau }\right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }_e^{\tau }\right\}}\end{array}$ (33)

We can also write (32) as:

$\left[M\right]\left\{{\ddot{\delta }}^u\right\}+\left\{F^{\mathrm{<mo>\; *\; </mo>}}\right\}=\left\{F\right\}+\left\{P\right\}$ (34)

$\left\{P\right\}$ is a vector of secondary variables and $\left\{F^{\mathrm{<mo>\; *\; </mo>}}\right\}=\left[B\right]\left\{{\delta }^{\tau }\right\}$ .

This completes the derivation of (25).

Equation (34) is a system of $n_u$ coupled nonlinear ODEs in degrees of freedom $\left\{{\ddot{\delta }}^u\right\}$ and $\left\{{\delta }^{\tau }\right\}$ . $\left[M\right]$ is not a diagonal matrix. However, if we approximate M by a diagonal matrix $\left[m\right]$ , then we can write (34) as (using, $m_i$ as diagonal element of $\left[M\right]$ ):

$m_i{\ddot{\delta }}_i^u+F_i^{\mathrm{<mo>\; *\; </mo>}}=F_i+P_i\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}{n}_u$ (35)

From (34) a typical ODE can be viewed as (with harmonic excitation):

$m\ddot{\delta }+F^{\mathrm{<mo>\; *\; </mo>}}=f_0\sin \left(\omega t\right)$ (36)

We note that $F^{\mathrm{<mo>\; *\; </mo>}}$ in (36) has no connection with $F\left(t\right)$ in (22), thus we cannot write (23) by simply replacing $F^{\mathrm{<mo>\; *\; </mo>}}$ in (36) by $F\left(t\right)$ in (22). It is conclusive that (23) has no real physics and its derivation has questionable and serious issues. Rest of the details that follow Equation (23) in reference [9] are of little consequence.

Remarks

1) Constitutive Equation (17) is valid form for 1D case (based on (16)), except that $\tau$ must be function of space $x_1$ and time t instead of just time (otherwise strain has no meaning) i.e.,

$\tau \left(x_1\mathrm{,}t\right)+{\lambda }_1\frac{\partial \tau \left(x_1\mathrm{,}t\right)}{\partial t}=E\epsilon \left(x_1\mathrm{,}t\right)+\eta \stackrel{\dot{}}{\epsilon }\left(x_1\mathrm{,}t\right)$ (37)

Now, there is no need to go through (19) - (22) (phenomenological). We can actually define Green’s strain and its time derivative.

${\epsilon }_{11}\left(x\mathrm{,}t\right)=\frac{\partial u_1}{\partial x_1}+\frac{1}{2}{\left(\frac{\partial u_1}{\partial x_1}\right)}^2$ (38)

${\stackrel{\dot{}}{\epsilon }}_{11}\left(x\mathrm{,}t\right)=\frac{\partial }{\partial t}\left({\epsilon }_{11}\left(x\mathrm{,}t\right)\right)=\frac{\partial v_1}{\partial x_1}+\frac{\partial u_1}{\partial x_1}\frac{\partial v_1}{\partial x_1}\mathrm{<mo>\; ;\; </mo>}{v}_1=\frac{\partial u_1}{\partial t}$ (39)

2) Force balance is not defined by (23), but instead by balance of linear momenta (as we have shown):

${\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)\tau \left(x_1\mathrm{,}t\right)\right)=0$ (40)

In (37) and (40), $\tau \left(x\mathrm{,}t\right)$ is deviatoric contravariant second Piola-Kirchhoff stress. Now, instead of (22) and (23) as in reference [9] , we have (37) and (40).

We clearly see that $\tau \left(x\mathrm{,}t\right)$ from (37) cannot be substituted in (27). Hence, it is not possible to obtain balance of linear momenta purely in terms of displacement and its spatial and time derivatives. As a consequence, we cannot derive mass, stiffness and damping matrices from the BLM by using GM/WF for a spatial discretization (space-time decoupled) as in case of TVES with rheology.

3) Definition of force F using (18) and its use in (23) both are erroneous. We remark that in BLM (force balance) gradients of stresses appear, not forces defined using stresses based on some area (or unit area).

4) We had to present these details to illustrate that the mathematical model used in ref [9] has many inconsistencies and errors, hence it is not a valid phenomenological mathematical for 1D physics of TVES with memory. This model has been used subsequently in published works.

5) We keep in mind that 1D models of any kind that are purely in time cannot be used to study nonlinear dynamics of continuous TVES with memory in ${ℝ}^1$ or ${ℝ}^2$ and ${ℝ}^3$ . As we have emphasized, 1D phenomenological models in time are not even valid for 1D continua media [14] .

In references [10] [11] , the authors utilize exactly the same approach as in reference [9] for the mathematical model but the model uses fractional derivatives. The issues discussed regarding reference [9] apply to references [10] [11] as well. In reference [12] , authors utilize Hamilton’s principle and Rayleigh dissipation functional to establish mathematical model for TVES without memory. The resulting mathematical model cannot be supported by CBL of CCM, hence is thermodynamically inconsistent. In reference [13] authors consider the same 1D constitutive theory purely in time as used in reference [9] . The authors express $\epsilon \left(t\right)$ and $\tau \left(t\right)$ in terms of trigonometric series in $\omega t$ using different coefficients for $\epsilon \left(t\right)$ and $\sigma \left(t\right)$ . Why are nonlinear strains and stresses a result of superposition of various terms in Equations (2) and (3) in reference [10] ? In any case, if we proceed further we find that damping matrix is diagonalized (Equation (48) in ref [10] ). None of these steps are supported by CCM and the modal basis transformations are only applicable in linear dynamics. In our view, works reported in reference [10] are of little consequence due to the fact that for continuous TVES with memory, a constitutive model will have space as well as time, hence the treatment in reference [10] (Equations (1) - (16) on page 276) is of little consequence for TVES matter with memory. We remark that derivation of mathematical model in reference [10] (in Section 3, beginning with page 278) based on potential energy is not supported by calculus of variations and mathematical classification of differential operators [4] [5] . Authors in reference [10] eventually derive matrix and vector forms of equations of motion (Equation (44)) in which mass, stiffness and damping matrices appear explicitly. We have shown that this is not possible for TVES with rheology due to the fact that the constitutive theory for deviatoric stress is a differential equation in time, hence deviatoric stresses from the constitutive theory cannot be substituted in BLM, an absolutely essential step to obtain discretized BLM in which mass, stiffness and damping matrices appear explicitly when using space-time decoupled method with GM/WF for the spatial discretization.

In another reference [11] , authors present “a review of nonlinear dynamics of hyperelastic structures”. Our view is that if the physics has nonlinearity and irreversibility due to dissipation, use of energy methods cannot be supported by calculus of variations and differential operator classification for IVPs as well as BVPs. In reference [12] dynamics of clamped rubber plates are investigated using energy methods. Few other publications also follow similar approaches discussed above. To our knowledge, complete mathematical models consisting of conservation and balance laws of CCM and constitutive theories derived using entropy inequality and representation theorem for TVE solids undergoing finite deformation and finite strain (hence nonlinear mechanisms of elasticity, dissipation and rheology) and the methods of obtaining their solutions have not been investigated in the published works.

In a recent paper Surana et al. [13] presented thermodynamically and mathematically consistent mathematical model in Lagrangian description using CCM for thermoviscoelastic solids without memory using contravariant second Piola-Kirchhoff stress tensor, Green’s strain tensor (a covariant measure) and rates of Green’s strain tensor up to order n. This mathematical model was used to derive a space-time decoupled finite element formulation for plates and shells, thus incorporating finite deformation, finite strain and rates of finite strain up to order n with formulations. In this formulation, stiffness as well as damping matrices are up to quadratic functions of the dofs, hence are nonlinear. Rayleigh damping is also discussed and incorporated in the mathematical model. An iterative solution procedure is also presented for obtaining solutions of the BVP associated with the corresponding IVP.

In another recent paper by Surana et al. [14] various aspects of finite deformation, finite strain nonlinear dynamics and dynamic bifurcation in TVES without memory are considered in detail using currently used 1D phenomenological models as well as using the thermodynamically and mathematically consistent mathematical models based on CCM presented in reference [13] . The work presented in this paper is an extension, modification and application of many concepts and methodologies presented by Surana et al. in reference [13] , hence we summarize significant aspects of the findings and the conclusion of reference [13] in the following.

1) Considering the mathematical model for TVEs without memory based on CCM, if we use space-time decoupled finite element method in conjunction with GM/WF in space, and then we obtain a system of nonlinear ODEs in time.

$\left[M\right]\left\{\ddot{\delta }\right\}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[C_i\right]\left\{{\delta }_{\left[i\right]}\right\}+\left[K\right]\left\{\delta \right\}=\left\{F\right\}+\left\{P\right\}$ (41)

In (41) material derivative of Green’s strain rate up to order n have been used in describing dissipation physics. Matrices $\left[K\right]$ and $\left[C_i\right]$ consists of the addition of following three matrices.

$\left[C_i\right]=\left[C_{\left[i\right]}^1\right]+\left[C_{\left[i\right]}^2\right]+\left[C_{\left[i\right]}^3\right]$ (42)

$\left[K\right]=\left[K^1\right]+\left[K^2\right]+\left[K^3\right]$ (43)

$\left\{\delta \right\}$ in (41) are the degrees of freedom for the spatial discretization and $\left\{{\delta }_{\left[i\right]}\right\}$ is the time derivative of $\left\{\delta \right\}$ of order i.

Matrices $\left[K^1\right]\mathrm{,}\left[C_{\left[i\right]}^1\right]$ have constant coefficients, coefficients of $\left[C_{\left[i\right]}^2\right]\mathrm{,}\left[K^2\right]$ and $\left[C_{\left[i\right]}^3\right]\mathrm{,}\left[K^3\right]$ are linear and quadratic functions of $\left\{\delta \right\}$ . $\left[K^1\right]\mathrm{,}\left[K^3\right]\mathrm{,}\left[C_{\left[i\right]}^1\right]\mathrm{,}\left[C_i^3\right]$ are symmetric, but $\left[K^2\right],\left[C_{\left[i\right]}^2\right]$ are non symmetric.

2) When integrating (41) in time using methods such as Newmark linear acceleration with Newton’s linear method for each time step, we find that computation of incremental solution of $\left\{\Delta \delta \right\}$ requires inverse of a symmetric “tangent matrix $\left[H_T\right]$ ” that consists of the sum of $\left[M\right]$ , the tangent stiffness matrix $\left[K_T\right]$ , the tangent damping matrix $\left[C_T\right]$ and another term containing variation of damping matrix with appropriate multipliers (Equation (94) in ref [14] ). $\left[M\right]$ is a constant coefficient consistent mass matrix. Tangent stiffness matrix $\left[K_T\right]$ consists of addition of four symmetric matrices.

$\left[K_T\right]=\left[{\underset{˜}{K}}^1\right]+\left[{\underset{˜}{K}}^2\right]+\left[{\underset{˜}{K}}^3\right]+\left[K_{\sigma }\right]$ (44)

Coefficients of $\left[K^1\right]$ are constant and those of $\left[{\underset{˜}{K}}^2\right]$ and $\left[{\underset{˜}{K}}^3\right]$ are linear and quadratic functions of $\left\{\delta \right\}$ . Matrix $\left[K_{\sigma }\right]$ is the stiffness matrix due to presence of stress field. Symmetric tangent damping matrix $\left[C_T\right]$ has the same structure as $\left[K_T\right]$ (44). $\left[M\right]\mathrm{,}\left[C_T\right]$ and $\left[K_T\right]$ clearly identify the physics due to stiffness, damping and mass that influence nonlinear dynamics and dynamic bifurcation phenomenon in TVES without memory. We remark that this physics is implicitly present in the mathematical model based on CBL of CCM and the constitutive theories. However, the details presented above show that explicit dependence of nonlinear dynamic response on various aspects of the physics is only possible due to the use of Newmark time integration in conjunction with Newton’s linear method for each time step.

3) Stiffness matrix due to stress field $\left[K_{\sigma }\right]$ is shown to play a key role in static as well as dynamic bifurcation. The static bifurcation is purely controlled by $\left[K_T\right]$ and can only exist if $\left[K_{\sigma }\right]$ is negative (due to compressive stress field). When $\left[K_{\sigma }\right]$ is positive (due to tensile stress field), static bifurcation is not possible.

4) In the case of nonlinear dynamics, $\left[M\right]\mathrm{,}\left[C_T\right]$ and $\left[K_T\right]$ , all three control the evolution (solutions of IVP), hence dynamic bifurcation depends on:

a) Tangent stiffness matrix $\left[K_T\right]$ which consists of $\left[K_{\sigma }\right]$ . $\left[K^1\right],\left[{\underset{˜}{K}}^2\right]$ and $\left[{\underset{˜}{K}}^3\right]$ .

b) $\left[C_T\right]$ i.e., tangent damping matrix similar to $\left[K_T\right]$ defining dissipation mechanism. Dissipation provides resistance to motion, hence higher dissipation contributes to inhibiting dynamic bifurcation and vice versa.

c) Translational inertial physics due to $\left[M\right]$ causes reduction in stiffness, hence enhances dynamic response and therefore enhances likelihood of the existence of dynamic bifurcation phenomenon.

d) It was shown that by an appropriate combination of the magnitude of the applied force and damping, dynamic bifurcation presence or absence can be influenced.

5) The existence of static bifurcation is not a necessary condition for dynamic bifurcation.

6) 1D nonlinear phenomenological models in time can never describe nonlinear dynamics and dynamic bifurcation of continuous systems in which space and time are always intrinsically coupled and present. In space-time decoupled methods, nonlinear ODEs in time preserve spatial details through $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ , whereas a single ODE in time, a lumped description in space, cannot describe the nonlinear dynamic response of continuous media as all the physics related to space is lost in this description.

7) In TVES with finite deformation, finite strain, dissipation mechanism converts some mechanical work into entropy which is non recoverable, hence the deformation process in nonlinear dynamics of TVES is irreversible, making the solutions of the corresponding IVPs path dependent. The severity of path dependency upon the amount of entropy production is controlled by the dissipation coefficient and magnitude of applied disturbance. In the linear dynamics of TVES, there is entropy production, but since the IVPs are linear, their solution cannot exhibit path dependency. Frequency response for damped linear systems is well known example of path independent frequency response.

8) The presence of nonlinearity in the PDEs and some irreversible physics is essential for the solutions of PDEs to be path dependent. In static bifurcation (TES), the strain and deformation are finite but there is the absence of irreversibility (conversion of some mechanical energy into entropy), hence static bifurcation processes defined by nonlinear PDEs are path independent, that is the solutions are independent of the increment of the load.

9) Using nonlinear dynamics of a TVE axial rod fixed at one end and subjected to harmonic excitation at the other end, authors in reference [14] illustrated significance of $\left[K_T\right]$ and in particular $\left[K_{\sigma }\right]$ , damping mechanism due to $\left[C_T\right]$ and the translational inertial physics due to $\left[M\right]$ and linear acceleration.

a) When the excitation force magnitude is sufficient to cause nonlinearities in the deformation physics in TVES, the frequency response is always path dependent due to nonlinearities in the IVP and the presence of irreversibility due to damping.

b) In the case of TVE axial rod, $\left[K_{\sigma }\right]$ is always negative, hence dynamic bifurcation always occurs to the left of the peaks.

c) When $\left[K_{\sigma }\right]$ is positive (bending of plates, shells, beams etc.) the dynamic bifurcation will occur to the right of the peak.

d) Path dependency results in bifurcation at different frequencies in $L\to R$ ( ${\omega }_1\to {\omega }_n$ ) and $L\leftarrow R$ ( ${\omega }_1\leftarrow {\omega }_n$ ) paths.

e) Due to the nature of external load (harmonic excitation) the rod is always in steady compression due to which the positive and the negative peak amplitudes differ when the response has become cyclic and repetitive. Negative peaks are always larger in magnitude than the positive peak, thus negative peaks show bifurcation more distinctly. Both peaks naturally show bifurcation at the same frequency.

f) The magnitude of the applied force and damping coefficient are critical in dynamic bifurcation. Lower damping and higher force are the most preferred choice for the likelihood of the existence of dynamic bifurcation. For a fixed force, the likelihood of the existence of bifurcation is enhanced by progressively reduced damping. Likewise for a fixed damping, bifurcation can be facilitated by increasing force.

10) There are many concerns regarding the validity and rigor of some of the mathematical models and solution techniques used to study nonlinear dynamics and dynamic bifurcation physics of continuous solid media in the published works. Some of these are summarized in the following.

a) We have established that 1D phenomenological mathematical model like (53) of ref [14] cannot possibly describe the nonlinear dynamics and dynamic bifurcation physics of a continuous solid media.

b) We have clearly demonstrated that $k_3=-1$ and $k_3=1$ (in (53) of ref [14] ) do not correspond to negative $\left[K_{\sigma }\right]$ and positive $\left[K_{\sigma }\right]$ . The derivation of $\left[K_{\sigma }\right]$ in time integration of ODEs (1) with Newton’s linear method clearly demonstrates this.

c) The nonlinear ODEs (1) resulting from the CBL of CCM after using space-time decoupled finite element method with GM/WF for the spatial discretization contains degrees of freedom $\left\{\delta \right\}$ and their time derivatives $\left\{\stackrel{\dot{}}{\delta }\right\},\left\{\ddot{\delta }\right\}$ that correspond to the nonlinear problem. These cannot be transformed using a modal basis determined from the linear form of ODEs (1) in which $\left[K\right]$ is constant. Authors in ref [14] presented details in the paper to show that forcing this transformation on (1) followed by the use of Rayleigh damping (see reference [18] [19] for issues with Rayleigh damping) leads to a system of decoupled ODEs that neither describes linear dynamics nor nonlinear dynamics. Hence, our view is that all nonlinear dynamic studies based on this approach are not supported by the mathematical models based on CBL of CCM followed by rather straightforward space-time coupled or space-time decoupled methods of obtaining their solutions, like we have presented in this paper and ref [13] [14] .

d) Transformation of (1) to modal basis using linear modes of vibrations derived using linear form of (1) is obviously flawed mathematically, but there is another significant issue in this approach. In this approach, the solution of a nonlinear problem reduces to the superposition of linear modes of vibrations using modal participation factors that are determined using nonlinear ODEs in the participation factors, Equation (18) of ref [14] . Superposition of any sort cannot be entertained for a nonlinear problem. This is another serious and fundamental problem in this approach.

e) In many published works, the trigonometric series with undetermined coefficients are used as solutions of the mathematical models in nonlinear dynamics. These solutions are obviously based on superposition which is not valid for a nonlinear system. Furthermore, the trigonometric functions are chosen such that their linear combinations with undetermined coefficients satisfy the necessary BCs and ICs. The values of the undetermined coefficients are established based on the fact that this solution must be in agreement with some measured experimental response of interest. This approach is more like establishing an analytical expression using the mathematical model that matches experimental results, in other words, similar to curve fit to the experimental data. From the solutions of ODEs and PDEs, we know that the solution consists of the sum of complementary and particular solutions. Complementary solution has undetermined coefficients (as many as supported by BCs and ICs depending upon BVP or IVP) and satisfies homogeneous form of the equations. Undetermined coefficients in the complementary solution are determined by ensuring that the total solution satisfies all BCs and ICs. In this approach, the solution of the mathematical model has nothing to do with experimental measurements i.e., it is completely independent of the experiment. A comparison of the solution calculated using this approach with the experimental measurements is now helpful in determining if the mathematical model and the experiment are addressing the same physics (of course barring errors in either approach).

2. Scope of Work

In this paper, we derive mathematical model for finite deformation, finite strain nonlinear dynamics of TVES with memory in Lagrangian description that consists of 1) conservation and balance laws of CCM in which deviatoric second Piola-Kirchhoff stress tensor and rate of covariant Green’s strain tensor are rate of work conjugate pair. 2) The constitutive theories for deviatoric second Piola-Kirchhoff stress tensor and heat vector are derived using conjugate pairs in the entropy inequality and the representation theorem. The dissipation mechanism is based on material derivatives of the Green’s strain tensor up to order n. The rheology mechanism requires the constitutive theory for the deviatoric second Piola-Kirchhoff stress tensor to be a differential equation in time [6] . Generalization of this leads to consideration of the material derivatives of the deviatoric second Piola-Kirchhoff stress up to order m in which the m^th material derivative is considered as constitutive tensor, remaining material derivatives up to order (m − 1) become the argument tensors of m^th material derivative.

Two approaches of obtaining the solutions of the mathematical model (IVPs) are discussed in the paper: 1) space-time coupled finite element approach based on a space-time strip or a space-time slab with time marching. The approach is based on space-time residual functional (STRF) and is unconditionally stable for all classes of space-time differential operators [4] , 2) a space-time decoupled finite element method with GM/WF for the spatial discretization. This approach results in a system of nonlinear ODEs in time that can integrated in time using Newmark linear acceleration method with Newton’s linear method for each increment of time (described in detail in reference [14] ). Merits and shortcomings of these two methodologies of obtaining solution of IVPs are well documented in ref [4] .

Various factors influencing nonlinear dynamics and dynamic bifurcation in TVES with memory are identified and their influence on the nonlinear dynamic response and dynamic bifurcation is discussed and illustrated using the mathematical model as well as the model problem studies. The results of the model problem studies considered in this paper are also compared with those obtained for TVES without rheology to illustrate the influence of rheology on nonlinear dynamic response and dynamic bifurcation phenomenon.

3. Preliminary Considerations

Thermoviscoelastic solids with dissipation and memory mechanisms are polymeric solids (rubbers or rubber like materials). Such materials contain short chain molecules of the solvent and the long chain molecules of the polymer. Both provide elasticity and dissipation mechanisms, but the rheology (memory) is only due to the presence of long chain molecules [20] [21] .

If the composition of the polymeric solid is dominated by the short chain molecules, then the memory mechanism is weak and we refer to these as non-dense polymeric solids. If the polymeric solid composition is heavily dominated by the long chain molecules then we refer to them as dense polymeric solids with pronounced rheology.

When an external stimulus is applied to a polymeric solid, the long chain molecules in the coiled relax state begin to unwind. In doing so these must overcome the viscous drag due to solvent as well as polymer. Excessively large deformation of rubber like materials is due to the unwinding of the long chain molecules. Upon cessation of external stimulus, the stretched long chain molecules begin to retreat to their relaxed (stress free) coiled state. In doing so they must also overcome the viscous, resistive forces due to the solvent and the polymer. The result is that, after cessation of external stimulus, it takes finite amount of time for the long chain molecules to resume their stress free state. The physical time to achieve relaxed state is a function of a time constant of material referred to as relaxation time. It can be shown that relaxation modulus that is a function of relaxation time can be used to determine clock time for a polymeric solid to achieve complete stress relaxation. Such materials are referred to as having memory or rheology. A material with larger relaxation time takes longer to achieve stress free state as opposed to a polymeric solid with smaller relaxation time. This physics of rheology is often termed as short term memory physics.

In this paper we consider CBL of CCM with finite deformation, finite strain, compressible, non-isothermal physics. Constitutive theories for volumetric change, distortional change, dissipation and memory mechanisms are derived using entropy inequality in conjunction with representation theorem. This mathematical model (in ${ℝ}^3$ ) has strict adherence to CBL of CCM and the constitutive theories are strictly based on entropy inequality and representation theorem. Thus, this mathematical model is thermodynamically and mathematically consistent.

The notations and symbols used in this paper are same as in reference [1] [2] and the previous recent papers by Surana et al. [13] [14] , hence their explanation is not repeated in this paper, but list of nomenclature is provided. We consider CBL of CCM and constitutive theories in Lagrangian description. Consideration of finite deformation and finite strain necessitates use of contravariant second Piola-Kirchhoff stress tensor ${\sigma }^{\left[0\right]}$ and covariant Green’s strain tensor ${\epsilon }_{\left[0\right]}$ as well as their material derivatives (same as ordinary time derivative or convected time derivatives) of up to order m and n respectively, ${\sigma }^{\left[j\right]}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}m$ and ${\epsilon }_{\left[i\right]}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n$ . A quantity Q with over bar $\bar{Q}$ is either defined in the current configuration or is an Eulerian description of Q, hence depends upon deformed coordinates $\bar{x}$ and time t i.e. $\bar{Q}=\bar{Q}\left(\bar{x}\mathrm{,}t\right)$ . Likewise $Q=Q\left(x\mathrm{,}t\right)$ is Lagrangian description of Q in the current configuration that depends on undeformed coordinates x and time t. We follow this throughout the paper.

We remark that in Lagrangian description density $\rho \left(x\mathrm{,}t\right)$ in the current configuration is deterministic using conservation of mass, $\rho \left(x\mathrm{,}t\right)=\frac{{\rho }_0}{\left|J\right|}$ . Hence, density $\rho \left(x\mathrm{,}t\right)$ is not a dependent variable in the mathematical model. Thus, conservation of mass (45) cannot be considered as part of the mathematical model, hence (45) adds no additional equations to the mathematical model. Secondly, balance of angular momenta (47) only establishes symmetry of the Cauchy stress tensor ${\sigma }^{\left(0\right)}$ , hence symmetry of ${\sigma }^{\left[0\right]}$ but adds no additional equations to the mathematical model. Entropy inequality (11) is a condition that needs to be satisfied for thermodynamic equilibrium of the deforming matter, but it also adds no additional equations to the mathematical model. Thus, the mathematical model for non isothermal physics consists of BLM (3 equations) and energy equation (1 equation), a total of four PDEs in thirteen dependent variables: $u$ (3), ${\sigma }^{\left[0\right]}$ (6), $q\left(3\right)$ , $\theta \left(1\right)$ , thus nine additional equations are needed to provide closure to this mathematical model. These are obtained from the constitutive theories for ${\sigma }^{\left[0\right]}$ (6) and $q\left(3\right)$ .

4. Complete Mathematical Model

Complete mathematical model in Lagrangian description for finite deformation, finite strain nonlinear dynamics of TVES with rheology consists of: 1) CBL derived using CCM in which contravariant second Piola-Kirchhoff stress tensor and covariant Green strain tensor and their material derivatives of up to order m and n (respectively) and 2) and the desired constitutive theories derived using entropy inequality in conjunction with representation theorem. We present details of the CBL of CCM for finite deformation, finite strain and details of the derivation of the constitutive theories in the following sections.

4.1. Conservation and Balance Laws of CCM

Following references [1] [2] , we can derive the following in the current configuration for conservation of mass, balance of linear momenta, balance of angular momenta, first and second laws of thermodynamics based on CCM in Lagrangian description.

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

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

$\left[{\sigma }^{\left[0\right]}\right]={\left[{\sigma }^{\left[0\right]}\right]}^{\text{T}}$ (47)

${\rho }_0\frac{De}{Dt}+\nabla \cdot g-{\sigma }^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}=0$ (48)

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

The symbols used in (45) - (49), conservation of mass, balance of linear momenta, balance of angular momenta, first and second law of thermodynamics are defined in nomenclatures. $\left[J\right]$ , the deformation gradient tensor, $g$ the temperature gradients and $\left\{u\right\}$ the displacement vector are defined as:

$\left[J\right]=\left[\frac{\partial \left\{\bar{x}\right\}}{\partial \left\{x\right\}}\right]=\left[\left[I\right]+\left[\frac{\partial \left\{u\right\}}{\partial \left\{x\right\}}\right]\right]$ (50)

$g=\nabla \cdot \theta$ (51)

$\left\{u\right\}={\left[u_1\mathrm{,}{u}_2\mathrm{,}{u}_3\right]}^{\text{T}}$ (52)

In Lagrangian description, $\frac{D}{Dt}=\frac{\partial }{\partial t}$ holds. We remark that specific internal energy e and entropy density $\eta$ are not dependent variables in the mathematical model as:

$e\left(x\mathrm{,}t\right)=e\left(\rho \left(x\mathrm{,}t\right)\mathrm{,}\theta \left(x\mathrm{,}t\right)\right)$ (53)

Helmholtz free energy density $\Phi$ and $\eta$ are not dependent variables in the mathematical model as shown in Section 4.2.

4.2. Constitutive Theories

In this section, we present derivations of constitutive theories for contravariant second Piola-Kirchhoff stress tensor and heat vector. In a compressible polymeric solid (or any other deforming solid), the total deformation can be decomposed into volumetric and deviatoric. The volumetric deformation results in change of volume for a fixed mass without change in shape whereas distortion deformation results in distortion of the shape of the volume of matter without change in volume. Since both of these deformations are mutually exclusive, a single constitutive theory for ${\sigma }^{\left[0\right]}$ stress tensor cannot possibly be used to describe both deformation physics. This requires additive decomposition of ${\sigma }^{\left[0\right]}$ into equilibrium ( ${}_e\sigma {}^{\left[0\right]}$ ) and deviatoric ( ${}_d\sigma {}^{\left[0\right]}$ ) stress tensors.

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

The constitutive theory for ${}_e\sigma {}^{\left[0\right]}$ addresses volumetric deformation physics and the distortional physics is addressed by the constitutive theory for ${}_d\sigma {}^{\left[0\right]}$ . For a volume with fixed mass, a change in volume without change in mass results in the change of density of the matter in the volume. Also, a uniform temperature increase or decrease for a volume of matter with fixed mass will result in pure change in volume, thus influencing density. Thus, volumetric deformation is a function of density ( $\rho \left(x\mathrm{,}t\right)$ ) and temperature $\theta$ .

${}_e\sigma {}^{\left[0\right]}={}_e\sigma {}^{\left[0\right]}\left(\rho \left(x\mathrm{,}t\right)\mathrm{,}\theta \left(x\mathrm{,}t\right)\mathrm{,}t\right)$ (55)

Since $\rho \left(x\mathrm{,}t\right)$ is not a dependent variable in CBL in the Lagrangian description as it is deterministic from conservation of mass (45), the constitutive tensor ${}_e\sigma {}^{\left(0\right)}$ in (55) cannot include $\rho \left(x\mathrm{,}t\right)$ as its argument tensor. To correct this situation, we need to consider Eulerian description, in which density $\bar{\rho }\left(\bar{x}\mathrm{,}t\right)$ is a dependent variable in the CBL, hence can be included as argument tensor of ${}_e\sigma {}^{\left(0\right)}$ , equilibrium contravariant Cauchy stress tensor and we would have:

${}_e\sigma {}^{\left(0\right)}={}_e\sigma {}^{\left(0\right)}\left(\bar{\rho }\left(\bar{x}\mathrm{,}t\right)\mathrm{,}\bar{\theta }\left(\bar{x}\mathrm{,}t\right)\mathrm{,}t\right)\mathrm{.}$ (56)

Substituting (56) in (49), we obtain:

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

Conjugate pairs $\frac{q\cdot g}{\theta }\mathrm{,}{}_e\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}$ and ${}_d\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}$ in conjunction with axioms of constitutive theories [1] [2] suggest $q\mathrm{,}{}_e\sigma {}^{\left[0\right]}\mathrm{,}{}_d\sigma {}^{\left[0\right]}$ to be the possible choice of constitutive tensors with $g$ and ${\epsilon }_{\left[0\right]}$ to be the possible choice of their argument tensors. Additionally, $\theta$ is also included as an argument tensor of these constitutive tensors due to non isothermal physics. Thus, at this stage, we have (time t is assumed to be an argument for al constitutive tensors):

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

${}_e\sigma {}^{\left[0\right]}={}_e\sigma {}^{\left[0\right]}\left({\epsilon }_{\left[0\right]}\mathrm{,}\theta \right)$ (59)

$q=q\left(g\mathrm{,}\theta \right)$ (60)

Comparing (59) with (55) and (56), we clearly note that choice of ${\epsilon }_{\left[0\right]}$ as an argument tensor of ${}_e\sigma {}^{\left[0\right]}$ is not valid. Since the Lagrangian description does not have mechanism of deriving constitutive theory for ${}_e\sigma {}^{\left[0\right]}$ , insistence on doing so i.e., substitution of (54) in (49), has caused this situation. We can continue further using (49) if we keep in mind that ${}_e\sigma {}^{\left[0\right]},{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}$ is not a valid rate of work conjugate pair and that the constitutive theory for volumetric deformation physics must be derived using Eulerian description of CBL.

The mechanism of dissipation requires at least ${\epsilon }_{\left[1\right]}$ to be argument tensor of ${}_d\sigma {}^{\left[0\right]}$ . If we assume that strain rates ${\epsilon }_{\left[i\right]}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n$ contribute to dissipation, then they all must be included as argument tensors of ${}_d\sigma {}^{\left[0\right]}$ . Thus, (58) can be modified:

${}_d\sigma {}^{\left[0\right]}={}_d\sigma {}^{\left[0\right]}\left({\epsilon }_{\left[0\right]}\mathrm{,}{\epsilon }_{\left[i\right]}\mathrm{,}\theta \right)\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n\mathrm{.}$ (61)

From the constitutive theories for polymeric fluid [1] [2] [6] , we know that presence of memory necessitates that the constitutive theory for the deviatoric stress tensor must at least be a first order differential equation in ${}_d\sigma {}^{\left[0\right]}$ and time i.e., in the constitutive theory for deviatoric second Piola-Kirchhoff stress tensor, we must at least consider ${}_d\sigma {}^{\left[0\right]}$ as well as ${}_d\sigma {}^{\left[1\right]}$ , with ${}_d\sigma {}^{\left[1\right]}$ as constitutive tensor having ${}_d\sigma {}^{\left[0\right]}$ as its argument tensor. We can generalize this choice by considering stress rate tensors ${}_d{\sigma }^{\left[k\right]}\mathrm{<mo>\; ;\; </mo>}k=\mathrm{1,2,}\cdots \mathrm{,}m$ and stress tensor ${}_d\sigma {}^{\left[0\right]}$ , with ${}_d\sigma {}^{\left[m\right]}$ as constitutive tensor and ${}_d{\sigma }^{\left[k\right]}\mathrm{<mo>\; ;\; </mo>}k=\mathrm{0,2,}\cdots \mathrm{,}m-1$ as its argument tensors. Thus, (61) can be written as:

${}_d\sigma {}^{\left[m\right]}={}_d\sigma {}^{\left[m\right]}\left({\epsilon }_{\left[0\right]}\mathrm{,}{\epsilon }_{\left[i\right]}\mathrm{,}{}_d\sigma {}^{\left[j\right]}\mathrm{,}\theta \right)\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ (62)

Choice of $\varphi$ and $\eta$ as additional constitutive tensors is valid at this stage based on entropy inequality (57).

The choice of the argument tensors of $\varphi$ and $\eta$ is not so straight forward at this stage, but can use principle of equipresence to begin with, thus we can write the following:

$\Phi =\Phi \left({\epsilon }_{\left[0\right]}\mathrm{,}{\epsilon }_{\left[i\right]}\mathrm{,}{}_d\sigma {}^{\left[j\right]},g\mathrm{,}\theta \right)$ (63)

$\eta =\eta \left({\epsilon }_{\left[0\right]}\mathrm{,}{\epsilon }_{\left[i\right]}\mathrm{,}{}_d\sigma {}^{\left[j\right]},g\mathrm{,}\theta \right)$ (64)

${}_d\sigma {}^{\left[0\right]}={}_d\sigma {}^{\left[0\right]}\left({\epsilon }_{\left[0\right]}\mathrm{,}{\epsilon }_{\left[i\right]}\mathrm{,}{}_d\sigma {}^{\left[j\right]}\mathrm{,}\theta \right)$ (65)

$q=q\left(g\mathrm{,}\theta \right)$ (66)

where $i=\mathrm{1,2,}\cdots \mathrm{,}n\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ . ${}_e\sigma {}^{\left[0\right]}$ is omitted in the above list as constitutive theory for ${}_e\sigma {}^{\left[0\right]}$ cannot be derived in Lagrangian description.

Using arguments of $\Phi$ , we can obtain material derivative of $\Phi$ ,

$\frac{D\Phi }{Dt}=\stackrel{\dot{}}{\Phi }=\frac{\partial \Phi }{\partial \left({\epsilon }_{\left[0\right]}\right)}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{\partial \Phi }{\partial \left(e_{\left[i\right]}\right)}:{\stackrel{\dot{}}{\epsilon }}_{\left[i\right]}+{\displaystyle \underset{j=0}{\overset{m-1}{\sum }}}\frac{\partial \Phi }{\partial \left({}_d\sigma {}^{\left[j\right]}\right)}:{}_d\stackrel{\dot{}}{\sigma }{}^{\left[j\right]}+\frac{\partial \varphi }{\partial g}\cdot \stackrel{\dot{}}{g}+\frac{\partial \Phi }{\partial \theta }\stackrel{\dot{}}{\theta }$ (67)

substituting $\frac{D\Phi }{Dt}$ from (67) into (57):

$\begin{array}{l}{\rho }_0(\frac{\partial \Phi }{\partial \left({\epsilon }_{\left[0\right]}\right)}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{\partial \Phi }{\partial \left(e_{\left[i\right]}\right)}:{\stackrel{\dot{}}{\epsilon }}_{\left[i\right]}+{\displaystyle \underset{j=0}{\overset{m-1}{\sum }}}\frac{\partial \Phi }{\partial {}_d\sigma {}^{\left[j\right]}}:{}_d\stackrel{\dot{}}{\sigma }{}^{\left[j\right]}\\ +\frac{\partial \varphi }{\partial g}\cdot \stackrel{\dot{}}{g}+\frac{\partial \Phi }{\partial \theta }\stackrel{\dot{}}{\theta })+\frac{q\cdot g}{\theta }-{}_e\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}-{}_d\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\le 0\end{array}$ (68)

Regrouping terms in (68):

$\begin{array}{l}\left({\rho }_0\frac{\partial \Phi }{\partial \left({\epsilon }_{\left[0\right]}\right)}-{}_e\sigma {}^{\left[0\right]}\right):{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\rho }_0\frac{\partial \Phi }{\partial \left(e_{\left[i\right]}\right)}:{\stackrel{\dot{}}{\epsilon }}_{\left[i\right]}+{\displaystyle \underset{j=0}{\overset{m-1}{\sum }}}{\rho }_0\frac{\partial \Phi }{\partial {}_d\sigma {}^{\left[j\right]}}:{}_d\stackrel{\dot{}}{\sigma }{}^{\left[j\right]}\\ +{\rho }_0\left(\eta +\frac{\partial \Phi }{\partial \theta }\right)\stackrel{\dot{}}{\theta }+{\rho }_0\frac{\partial \varphi }{\partial g}\cdot \stackrel{\dot{}}{g}+\frac{q\cdot g}{\theta }\le 0\end{array}$ (69)

The entropy inequality (69) holds for arbitrary but admissible ${\stackrel{\dot{}}{\epsilon }}_{\left[i\right]}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n$ ; ${}_d\sigma {}^{\left[j\right]}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ , $\stackrel{\dot{}}{\theta }$ , and $\stackrel{\dot{}}{g}$ if their coefficients are set to zero i.e., if the following holds:

${\rho }_0\frac{\partial \Phi }{\partial {\stackrel{\dot{}}{\epsilon }}_{\left[i\right]}}=0\Rightarrow \Phi \ne \Phi \left(e_{\left[i\right]}\right)\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n$ (70)

${\rho }_0\frac{\partial \Phi }{\partial {}_d\sigma {}^{\left[j\right]}}=0\Rightarrow \Phi \ne \Phi \left({}_d\sigma {}^{\left[j\right]}\right)\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1},\cdots \mathrm{,}m-1$ (71)

${\rho }_0\left(\eta +\frac{\partial \Phi }{\partial \theta }\right)=0\Rightarrow \eta =-\frac{\partial \Phi }{\partial \theta }$ (72)

$\frac{\partial \Phi }{\partial g}=0\Rightarrow \Phi \ne \Phi \left(g\right)$ (73)

Remarks

From (69) - (73), we can conclude the following:

1) Equations (70) and (71) imply that $\Phi$ is not a function of ${\epsilon }_{\left[i\right]}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n$ and ${}_d\sigma {}^{\left[j\right]}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ .

2) Based on (73), $\Phi$ is not a function of $g$ either.

3) Equation (72) implies that $\eta$ is not a constitutive variable as $\eta$ is defined by $-\frac{\partial \Phi }{\partial \theta }$ .

4) Based on (73), $\Phi$ does not depend upon $g$ .

5) Thus, now we have the following constitutive tensors and their argument tensors:

$\Phi =\Phi \left({\epsilon }_{\left[i\right]}\mathrm{,}\theta \right)$ (74)

${}_d\sigma {}^{\left[m\right]}={}_d\sigma {}^{\left[m\right]}\left({\epsilon }_{\left[0\right]}\mathrm{,}{\epsilon }_{\left[i\right]}\mathrm{,}{}_d\sigma {}^{\left[j\right]}\mathrm{,}\theta \right)\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ (75)

$q=q\left(g\mathrm{,}\theta \right)$ (76)

and the entropy inequality (69) reduces to

$\left({\rho }_0\frac{\partial \Phi }{\partial {\epsilon }_{\left[0\right]}}-{}_e\sigma {}^{\left[0\right]}\right):{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}-{}_d\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}+\frac{q\cdot g}{\theta }\le 0$ (77)

We cannot set the coefficient of ${\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}$ in the first term of (77) to zero as it would imply:

${}_e\sigma {}^{\left[0\right]}={\rho }_0\frac{\partial \Phi \left({\epsilon }_{\left[0\right]}\mathrm{,}\theta \right)}{\partial {\epsilon }_{\left[0\right]}}$ (78)

which implies that ${}_e\sigma {}^{\left[0\right]}$ depends upon ${\epsilon }_{\left[0\right]}$ , but ${}_e\sigma {}^{\left[0\right]}$ depends on density and temperature only. Thus, (78) is obviously erroneous. Hence, constitutive theory for ${}_e\sigma {}^{\left[0\right]}$ cannot be derived using (77). This of course is due to the fact that conservation of mass in Lagrangian description permits calculation of density $\rho \left(x\mathrm{,}t\right)$ , hence $\rho \left(x\mathrm{,}t\right)$ is not a dependent variable in the mathematical model. We must consider equivalent of the first term in (77) in Eulerian description to derive constitutive theory for ${}_e\sigma {}^{\left[0\right]}$ .

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

In the following, we derive constitutive theories for ${}_e\sigma {}^{\left[0\right]}$ for compressible as well as incompressible thermoviscoelastic solid with memory.

Compressible matter

We consider entropy inequality in Eulerian description [1] [2] (using contravariant Cauchy stress tensor):

$\bar{\rho }\left(\frac{D\bar{\Phi }}{Dt}+\bar{\eta }\frac{D\bar{\theta }}{Dt}\right)+\frac{\bar{q}\cdot \bar{g}}{\theta }-{\bar{\sigma }}^{\left(0\right)}:\bar{D}\le 0$ (79)

${\bar{\sigma }}^{\left(0\right)}={}_e\bar{\sigma }{}^{\left(0\right)}+{}_d\bar{\sigma }{}^{\left(0\right)}$ (80)

We can write the following constitutive tensors and their argument tensors for thermoviscous case in which dissipation is dependent on $\bar{D}$ (symmetric part of velocity gradient tensor):

${}_e\bar{\sigma }{}^{\left(0\right)}={}_e\bar{\sigma }{}^{\left(0\right)}\left(\bar{\rho }\left(\bar{x}\mathrm{,}t\right)\mathrm{,}\bar{\theta }\left(\bar{x}\mathrm{,}t\right)\right)\mathrm{,}$ (81)

${}_d\bar{\sigma }{}^{\left(0\right)}={}_d\bar{\sigma }{}^{\left(0\right)}\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{\theta }\right)\mathrm{,}$ (82)

$\bar{q}=\bar{q}\left(\bar{g}\mathrm{,}\bar{\theta }\right)$ (83)

$\begin{array}{l}\bar{\Phi }=\bar{\Phi }\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{g}\mathrm{,}\bar{\theta }\right)\\ \bar{\eta }=\bar{\eta }\left(\bar{\rho }\mathrm{,}\bar{D}\mathrm{,}\bar{g}\mathrm{,}\bar{\theta }\right)\end{array}$ (84)

Using $\bar{\Phi }$ in (84) we can write the following for $\frac{D\bar{\Phi }}{Dt}$ ,

$\frac{D\bar{\Phi }}{Dt}=\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}\stackrel{\dot{}}{\bar{\rho }}+\frac{\partial \bar{\Phi }}{\partial \bar{D}}:\stackrel{\dot{}}{\bar{D}}+\frac{\partial \bar{\Phi }}{\partial \bar{g}}\cdot \stackrel{\dot{}}{\bar{g}}+\frac{\partial \bar{\Phi }}{\partial \bar{\theta }}\stackrel{\dot{}}{\bar{\theta }}\mathrm{.}$ (85)

From continuity equation in Eulerian description:

$\frac{D\bar{\rho }}{Dt}=\stackrel{\dot{}}{\bar{\rho }}=-\bar{\rho }\left(\bar{\nabla }\cdot \bar{v}\right)=-\bar{\rho }{\bar{D}}_{kk}=-\bar{\rho }{D}_{kl}{\delta }_{lk}=-\bar{\rho }\bar{D}:d\mathrm{.}$ (86)

substituting (86) in (85) and regrouping terms:

$\begin{array}{l}\left(-{\bar{\rho }}^2\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}d-{}_e\bar{\sigma }{}^{\left(0\right)}\right):\bar{D}+\bar{\rho }\frac{\partial \bar{\Phi }}{\partial \bar{D}}:\stackrel{\dot{}}{\bar{D}}+\bar{\rho }\frac{\partial \bar{\Phi }}{\partial \bar{g}}\cdot \stackrel{\dot{}}{\bar{g}}\\ +\bar{\rho }\left(\bar{\eta }+\frac{\partial \bar{\Phi }}{\partial \bar{\theta }}\right)\stackrel{\dot{}}{\bar{\theta }}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}-{}_d\sigma {}^{\left(0\right)}:\bar{D}\le 0.\end{array}$ (87)

The entropy inequality (87) holds for arbitrary but admissible choices of $\stackrel{\dot{}}{\bar{D}}\mathrm{,}\stackrel{\dot{}}{\bar{g}}$ and $\stackrel{\dot{}}{\bar{\theta }}$ , if the following conditions are satisfied:

$\bar{\rho }\frac{\partial \bar{\Phi }}{\partial \bar{D}}=0\Rightarrow \bar{\Phi }\ne \bar{\Phi }\left(\bar{D}\right)\mathrm{,}$ (88)

$\bar{\rho }\frac{\partial \bar{\Phi }}{\partial \bar{g}}=0\Rightarrow \bar{\Phi }\ne \bar{\Phi }\left(\bar{g}\right)\mathrm{,}$ (89)

$\bar{\rho }\left(\bar{\eta }+\frac{\partial \bar{\Phi }}{\partial \bar{\theta }}\right)=0\Rightarrow \bar{\eta }=-\frac{\partial \bar{\Phi }}{\partial \bar{\theta }}\mathrm{.}$ (90)

From (90), we can conclude that $\bar{\eta }$ is not a constitutive variable as it is defined by $-\frac{\partial \bar{\varphi }}{\partial \bar{\eta }}$ . The arguments of $\bar{\Phi }$ reduces to

$\bar{\varphi }=\bar{\varphi }\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)$ (91)

and the entropy inequality (87) reduces to

$\left(-{\bar{\rho }}^2\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}\delta -{}_e\bar{\sigma }{}^{\left(0\right)}\right):\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}-{}_e\bar{\sigma }{}^{\left(0\right)}:\bar{D}\le 0.$ (92)

The first term in (92) can be thought of as equivalent of first term in (77) in Lagrangian description. Since ${}_e\bar{\sigma }{}^{\left(0\right)}={}_e\bar{\sigma }{}^{\left(0\right)}\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)$ , we can set the coefficient of $\bar{D}$ in the first term of (92) to zero to obtain:

${}_e\bar{\sigma }{}^{\left(0\right)}=-{\bar{\rho }}^2\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}\delta =\bar{p}\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)\delta$ (93)

$\bar{p}\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)=-{\bar{\rho }}^2\frac{\partial \bar{\Phi }}{\partial \bar{\rho }}$ (94)

in which $\bar{p}\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)$ is thermodynamic pressure for compressible matter. From (93) we can write (Lagrangian description of (93)):

${}_e\bar{\sigma }{}^{\left(0\right)}=p\left(\rho \left(x\mathrm{,}t\right)\mathrm{,}\theta \left(x\mathrm{,}t\right)\right)\delta$ (95)

recall the relationship between Cauchy stress ${\bar{\sigma }}^{\left(0\right)}$ and second Piola-Kirchhoff stress ${\sigma }^{\left[0\right]}$ [1] [2] :

$\left[{}_e\sigma {}^{\left[0\right]}\right]=\left|J\right|{\left[J\right]}^{-1}\left[{}_e\sigma {}^{\left[0\right]}\right]{\left[{\left[J\right]}^{-1}\right]}^{\text{T}}$ (96)

hold also. Substituting $\left[{}_e\sigma {}^{\left[0\right]}\right]$ from (95) in (96), we finally have constitutive theory for ${}_e\sigma {}^{\left[0\right]}$ for compressible matter.

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

We note that ${}_e\sigma {}^{\left[0\right]}$ is not a pressure field as is the case in (95).

Incompressible matter

For incompressible matter, density remains constant i.e., $\bar{\rho }\left(\bar{x}\mathrm{,}t\right)=\rho \left(x\mathrm{,0}\right)={\rho }_0$ . Thus, from conservation of mass (in Eulerian description), we have:

$\stackrel{\dot{}}{\bar{\rho }}=-\bar{\rho }\left(\bar{\nabla }\cdot \bar{v}\right)=0$ (CM) (98)

and

$\frac{\partial \bar{\Phi }\left(\bar{\rho }\mathrm{,}\bar{\theta }\right)}{\partial \bar{\rho }}=0.$ (99)

Thus, in this case, entropy inequality (92) reduces to

$-{}_e\bar{\sigma }{}^{\left(0\right)}:\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}-{}_d\sigma {}^{\left(0\right)}:\bar{D}\le 0.$ (100)

In order to derive the constitutive theory for ${}_e\bar{\sigma }{}^{\left(0\right)}$ for an incompressible matter, we introduce incompressibility condition (98) in (100).

From conservation of mass (98), we have have the following for incompressible case:

$\bar{\nabla }\cdot \bar{v}={\bar{D}}_{kk}={\bar{D}}_{kl}{\delta }_{lk}=\delta :\bar{D}=0.$ (101)

If (101) holds, then the following holds too:

$\bar{p}\left(\bar{\theta }\right)\delta :\bar{D}=0.$ (102)

$\bar{p}\left(\bar{\theta }\right)$ is a Lagrange multiplier. Adding (102) to (100) and regrouping terms:

$\left(\bar{p}\left(\bar{\theta }\right)\delta -{}_e\bar{\sigma }{}^{\left(0\right)}\right):\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}-{}_e\bar{\sigma }{}^{\left(0\right)}:\bar{D}\le 0.$ (103)

Entropy inequality (103) is satisfied for arbitrary but admissible $\bar{D}$ if the coefficient of $\bar{D}$ in the first term in (103) is set to zero.

${}_e\bar{\sigma }{}^{\left(0\right)}=\bar{p}\left(\bar{\theta }\right)\delta \mathrm{.}$ (104)

Equation (104) is the constitutive theory for ${}_e\bar{\sigma }{}^{\left(0\right)}$ for incompressible non isothermal physics in Eulerian description. In Lagrangian description, (104) reduces to

${}_e\sigma {}^{\left(0\right)}=p\left(\theta \left(x\mathrm{,}t\right)\right)\delta \mathrm{.}$ (105)

Using (97) and (105) we obtain the desired constitutive theory for ${}_e\sigma {}^{\left[0\right]}$ in

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

We note that ${}_e\sigma {}^{\left[0\right]}$ in this case is also not a pressure field. In case of isothermal physics, (106) holds but $p\left(\theta \right)$ is replaced by p, mechanical pressure.

Thus, the entropy inequality (103) in Eulerian description reduces (reduced form) to the following:

$-{}_d\bar{\sigma }{}^{\left(0\right)}:\bar{D}+\frac{\bar{q}\cdot \bar{g}}{\bar{\theta }}\le 0.$ (107)

Likewise in (77), the entropy inequality in Lagrangian description, the first term is eliminated and we have the following reduced form of entropy inequality in Lagrangian description:

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

4.2.2. Constitutive Theory for Deviatoric Contravariant Second Piola-Kirchhoff Stress Tensor

We have already established ${}_d\sigma {}^{\left[m\right]}$ as constitutive tensor and its argument tensors in (62):

${}_d\sigma {}^{\left[m\right]}={}_d\sigma {}^{\left[m\right]}\left({\epsilon }_{\left[0\right]}\mathrm{,}{\epsilon }_{\left[i\right]}\mathrm{,}\theta \right)i=\mathrm{1,2,}\cdots \mathrm{,}n\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ (109)

We derive constitutive theory for ${}_d\sigma {}^{\left[m\right]}$ using representation theorem [2] [22] - [29] . If S be the space of tensor ${}_d\sigma {}^{\left[m\right]}$ , then the basis of this space must be determined using argument tensors of ${}_d\sigma {}^{\left[m\right]}$ in (109). Since ${}_d\sigma {}^{\left[m\right]}$ is a symmetric tensor of rank two, we determine combined generators ${}^{\sigma }\underset{˜}{G}{}^i\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}N$ of the argument tensors of ${}_d\sigma {}^{\left[m\right]}$ in (109) that are symmetric tensors of rank two. The tensor $I$ and the generators ${}^{\sigma }\underset{˜}{G}{}^i\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}N$ constitute integrity (irreducible set) that forms the basis of the space S of constitutive tensor ${}_d\sigma {}^{\left[m\right]}$ . Thus, ${}_d\sigma {}^{\left[m\right]}$ can be expanded using a linear combination of the tensors in the integrity i.e., the basis of the space “S”.

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

If ${}^{\sigma }\underset{˜}{I}{}^j\mathrm{<mo>\; :\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}M$ are combined invariants of the argument tensors of ${}_d\sigma {}^{\left[m\right]}$ in (109), then the coefficients ${\underset{˜}{\alpha }}^i\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,}\cdots \mathrm{,}N$ in the linear combination can be functions of these invariants and temperature $\theta$ i.e.,

${\underset{˜}{\alpha }}^i={\underset{˜}{\alpha }}^i\left({}^{\sigma }\underset{˜}{I}{}^J\mathrm{,}\theta \right)\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}M\mathrm{.}$ (111)

Using (110) and (111) we can derive material coefficients [1] [2] . We expand ${\underset{˜}{\alpha }}^i\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,}\cdots \mathrm{,}N$ in Taylor series in ${}^{\sigma }\underset{˜}{I}{}^j;j=\mathrm{1,2,}\cdots \mathrm{,}M$ about a known configuration $\underset{\_}{\Omega }$ and retain only up to linear terms in ${}^{\sigma }\underset{˜}{I}{}^j;j=\mathrm{1,2,}\cdots \mathrm{,}M$ (for simplicity). Taylor series expansion is not considered in $\theta$ , as the stresses due to thermal field have already been accounted for in the equilibrium stress.

${}^{\sigma }\alpha {}^i={{}^{\sigma }\alpha {}^i|}_{\underset{\_}{\Omega }}+{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{\frac{\partial \left({}^{\sigma }\alpha {}^i\right)}{\partial \left({}^{\sigma }\underset{˜}{I}{}^j\right)}|}_{\underset{\_}{\Omega }}\left({}^{\sigma }\underset{˜}{I}{}^j-{{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\right)\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,}\cdots \mathrm{,}N\mathrm{.}$ (112)

Substituting (112) in (110) and collecting coefficients of $I\mathrm{,}{}^{\sigma }\underset{˜}{I}{}^{j}I\mathrm{,}{}^{\sigma }\underset{˜}{G}{}^i$ and ${}^{\sigma }\underset{˜}{I}{}^j\left({}^{\sigma }\underset{˜}{G}{}^i\right)\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}M$ and $i=\mathrm{0,1,}\cdots \mathrm{,}N$ and defining new coefficients using:

$\begin{array}{l}{{\sigma }^0|}_{\underset{\_}{\Omega }}=\left({{}^{\sigma }\alpha {}^0|}_{\underset{\_}{\Omega }}+{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{\frac{\partial \left({}^{\sigma }\alpha {}^0\right)}{\partial \left({}^{\sigma }\underset{˜}{I}{}^j\right)}|}_{\underset{\_}{\Omega }}\right)\left(-{{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\right)\mathrm{<mo>\; ;\; </mo>}{\underset{˜}{a}}_j={\frac{\partial \left({}^{\sigma }\alpha {}^0\right)}{\partial \left({}^{\sigma }\underset{˜}{I}{}^j\right)}|}_{\underset{\_}{\Omega }}\\ {\underset{˜}{b}}_i={{}^{\sigma }\alpha {}^i|}_{\underset{\_}{\Omega }}+{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{\frac{\partial \left({}^{\sigma }\alpha {}^i\right)}{\partial \left({}^{\sigma }\underset{˜}{I}{}^j\right)}|}_{\underset{\_}{\Omega }}\left(-{{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\right)\mathrm{<mo>\; ;\; </mo>}{\underset{˜}{c}}_{ij}={\frac{\partial \left({}^{\sigma }\alpha {}^i\right)}{\partial \left({}^{\sigma }\underset{˜}{I}{}^j\right)}|}_{\underset{\_}{\Omega }}\mathrm{.}\end{array}$ (113)

for $i=\mathrm{1,2,}\cdots \mathrm{,}N$ and $j=\mathrm{1,2,}\cdots \mathrm{,}M$ .

We can write (110) as follows:

${}_d\sigma {}^{\left[m\right]}={\underset{˜}{\sigma }}^{0}I+{\displaystyle \underset{j=0}{\overset{M}{\sum }}}{\underset{˜}{a}}_j\left({}^{\sigma }\underset{˜}{I}{}^j\right)I+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{\underset{˜}{b}}_i\left({}^{\sigma }\underset{˜}{G}{}^i\right)+{\displaystyle \underset{i=1}{\overset{N}{\sum }}}{\displaystyle \underset{j=1}{\overset{M}{\sum }}}{\underset{˜}{c}}_{ij}\left({}^{\sigma }\underset{˜}{I}{}^j\right)\left({}^{\sigma }\underset{˜}{G}{}^i\right)\mathrm{.}$ (114)

${\underset{˜}{a}}_j\mathrm{,}{\underset{˜}{b}}_i$ and ${\underset{˜}{c}}_{ij}$ are material coefficients. These are functions of ${{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}M$ and ${\theta |}_{\underset{\_}{\Omega }}$ . This constitutive theory requires M + NM + N material coefficients. ${\underset{˜}{\sigma }}^0$ is known initial stress in the known configuration $\underset{\_}{\Omega }$ . The constitutive theory (114) is based on complete basis of the space S of the constitutive tensor ${}_d\sigma {}^{\left[m\right]}$ .

4.2.3. Simplified Constitutive Theories for Contravariant Deviatoric Second Piola-Kirchhoff Stress Tensor

First, we consider an ordered rate theory of orders m and n that is linear in the components of tensor ${\epsilon }_{\left[i\right]}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,}\cdots \mathrm{,}n$ and ${}_d\sigma {}^{\left[j\right]}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ and the invariants ${}^{\sigma }\underset{˜}{I}{}^j\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ . Thus, in this constitutive theory, we use generators ${\epsilon }_{\left[i\right]}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n$ ; ${}_d\sigma {}^{\left[j\right]}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ and the invariants ${}^{\sigma }\underset{˜}{I}{}^j$ that are given by $tr\left({\epsilon }_{\left[i\right]}\right)$ , $tr\left({}_d\sigma {}^{\left[j\right]}\right)$ ; $i=\mathrm{0,1,}\cdots \mathrm{,}n$ and $j=\mathrm{0,1,}\cdots \mathrm{,}m-1$ . Thus, in this case (114) reduces to (redefining material coefficients):

$\begin{array}{c}{}_d\sigma {}^{\left[m\right]}={\underset{˜}{\sigma }}^{0}I+a_1{\epsilon }_{\left[0\right]}+a_{2}tr\left({\epsilon }_{\left[i\right]}\right)I+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{b}_i^1\left({\epsilon }_{\left[i\right]}\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{b}_i^{2}tr\left({\epsilon }_{\left[i\right]}\right)I\\ +{\displaystyle \underset{j=0}{\overset{m-1}{\sum }}}{c}_j^1\left({}_d\sigma {}^{\left[j\right]}\right)+{\displaystyle \underset{j=0}{\overset{m-1}{\sum }}}{c}_j^{2}tr\left({}_d\sigma {}^{\left[j\right]}\right)I\end{array}$ (115)

The material coefficients in (115) can still be functions of invariants ${{}^{\sigma }\underset{˜}{I}{}^j|}_{\underset{\_}{\Omega }}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}m$ and ${\theta |}_{\underset{\_}{\Omega }}$ . We rewrite (115) as follows (neglects ${\underset{˜}{\sigma }}^0$ without loss of generality):

$\begin{array}{l}{}_d\sigma {}^{\left[m\right]}-{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left(c_j^1\left({}_d\sigma {}^{\left[j\right]}\right)+c_j^{2}tr\left({}_d\sigma {}^{\left[j\right]}\right)I\right)\\ =a_1{\epsilon }_{\left[0\right]}+a_{2}tr\left({\epsilon }_{\left[i\right]}\right)I+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(b_i^1\left({\epsilon }_{\left[i\right]}\right)+b_i^{2}tr\left({\epsilon }_{\left[i\right]}\right)I\right)\end{array}$ (116)

Divide throughout by $c_0^1$ (assuming $c_0^1\ne 0$ ) and introduce new notation for the material coefficients.

$\begin{array}{l}-\frac{1}{c_0^1}=\lambda ;\frac{c_0^2}{c_0^1}={\lambda }^0\\ -\frac{a_1}{c_0^1}=2\underset{˜}{\mu };-\frac{a_2}{c_0^1}=-\underset{˜}{\lambda }\\ -\frac{b_i^1}{c_0^1}=2\eta ;-\frac{b_i^2}{c_0^1}={\kappa }_i;i=1,2,\cdots ,n\\ -\frac{c_j^1}{c_0^1}={\lambda }_j^1;-\frac{c_2^j}{c_0^1}={\lambda }_j^2;j=1,2,\cdots ,m-1\end{array}$ (117)

$\begin{array}{l}{}_d\sigma {}^{\left[0\right]}+\lambda {}_d\sigma {}^{\left[m\right]}+{\lambda }^{0}tr\left({}_d\sigma {}^{\left[0\right]}\right)I+{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}\left({\lambda }_j^1\left({}_d\sigma {}^{\left[j\right]}\right)+{\lambda }_j^{2}tr\left({}_d\sigma {}^{\left[j\right]}\right)I\right)\\ =2\underset{˜}{\mu }{\epsilon }_{\left[0\right]}+\underset{˜}{\lambda }\left(tr{\epsilon }_{\left[0\right]}\right)I+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(2{\eta }_i{\epsilon }_{\left[i\right]}+{\kappa }_{i}tr\left({\epsilon }_{\left[i\right]}\right)I\right)\end{array}$ (118)

The term ${\lambda }^{0}tr\left({}_d\sigma {}^{\left[0\right]}\right)I$ is neglected in the currently used constitutive theories.

Remarks

Based on this constitutive theory (118) we make following remarks:

1) Elasticity requires two material coefficients in finite deformation, finite strain case also (similar to Lames constants).

2) Each strain rate ${\epsilon }_{\left[i\right]}$ contributing to nonlinear dissipation mechanism also requires two material coefficients.

3) $\lambda$ is primary relaxation time. Each stress rate ${}_d\sigma {}^{\left[k\right]}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}m-1$ is associated with their own relaxation time defined by ${\lambda }_j^1$ and ${\lambda }_j^2$ . This is a rather new and important phenomenon, the existence of which is only justified and revealed due to derivations of ordered rate constitutive theory presented here. We explain this further. The polymer molecules generally connect to neighbors, forming colonies and a colony may also be interconnected to neighboring colonies. In addition to this, some polymer molecules may remain isolated in the solvent by themselves. In this complicated network of polymer molecules, upon cessation of external stimulus, it is highly unlikely that relaxation of all polymer molecules can be described by a single relaxation time. The constitutive theory suggests that ( ${\lambda }_j^1\mathrm{,}{\lambda }_j^2$ ) are progressively decreasing relaxation times that are associated with relaxation at multiple time scales. Relaxation time $\lambda$ is a bulk relaxation time (aggregate) that accounts for majority of this relaxation phenomenon. Generally, the terms ${\lambda }_j^{2}tr\left({}_d\sigma {}^{\left[j\right]}\right)I$ are also neglected in (118). Thus, we now have (by redefining ${\lambda }_j^1$ as ${\lambda }_j$ ).

$\begin{array}{l}{}_d\sigma {}^{\left[0\right]}+\lambda \left({}_d\sigma {}^{\left[m\right]}\right)+{\displaystyle \underset{j=1}{\overset{m-1}{\sum }}}{\lambda }_j\left({}_d\sigma {}^{\left[j\right]}\right)\\ =2\underset{˜}{\mu }{\epsilon }_{\left[0\right]}+\underset{˜}{\lambda }\left(tr{\epsilon }_{\left[0\right]}\right)I+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(2{\eta }_i{\epsilon }_{\left[i\right]}+{\kappa }_i\left(tr{\epsilon }_{\left[i\right]}\right)I\right)\end{array}$ (119)

In (119), we clearly observe elasticity material coefficients $2\underset{˜}{\mu }\mathrm{,}\underset{˜}{\lambda }$ and dissipation material coefficients $2{\eta }_i\mathrm{,}{\kappa }_i$ corresponding to strain rate ${\epsilon }_{\left[i\right]}$ . Bulk relaxation time $\lambda$ accompanied by hierarchy of relaxation times ${\lambda }_j$ associated with stress rates ${}_d\sigma {}^{\left[j\right]}$ . In the absence of ${\lambda }_j\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}m-1$ , the relaxation is homogeneous i.e., it occurs in the entire deformed matter uniformly and is uniformly controlled by bulk relaxation time $\lambda$ . Presence of relaxation times ${\lambda }_j\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}m-1$ suggest inhomogeneous relaxation process, perhaps more realistic due to complexity of the network of polymer molecules in the solvent constituting the composition of the matter.

4) Commonly used constitutive theories are a subset of further simplifications of (119) for $m=1,n=1$ . In (119), if we choose $m=1$ and $n=1$ , then we obtain:

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

5) Constitutive theories (118) or (119) or (120) can be easily written in matrix and vector form [1] [2] using Voigt’s notation.

4.2.4. Constitutive Theory for q

We consider $q=q\left(g\mathrm{,}\theta \right)$ and use representation theorem [2] [22] - [29] . The combined generators of the tensors $g$ and $\theta$ that are tensors of rank one is just $g$ and the combined invariants of $g$ and $\theta$ is $g\cdot g$ (or ${}^q\underset{˜}{I}$ ). Thus, we can write the following in the current configuration:

$q=-{}^q\underset{˜}{\alpha }g$ (121)

in which

${}^q\underset{˜}{\alpha }={}^q\underset{˜}{\alpha }\left({}^q\underset{˜}{I}\mathrm{,}\theta \right)$ (122)

Material coefficients in the constitutive theory for $q$ given by (122) are obtained by considering Taylor series expansion of ${}^q\underset{˜}{\alpha }$ in ${}^q\underset{˜}{I}$ about a known configuration $\underset{\_}{\Omega }$ and only retaining up to linear terms in ${}^q\underset{˜}{I}$ (for simplicity).

${}^q\underset{˜}{\alpha }={{}^q\underset{˜}{\alpha }|}_{\underset{\_}{\Omega }}+{\frac{\partial {}^q\underset{˜}{\alpha }}{\partial {}^q\underset{˜}{I}}|}_{\underset{\_}{\Omega }}\left({}^q\underset{˜}{I}-{{}^q\underset{˜}{I}|}_{\underset{\_}{\Omega }}\right)$ (123)

Substituting (123) in (121) and collecting coefficients of the terms defined in the current configuration yield (after introducing new notation for the terms in $\underset{\_}{\Omega }$ ):

$q=-{{\kappa }_1|}_{\underset{\_}{\Omega }}g-{{\kappa }_2|}_{\underset{\_}{\Omega }}\left(g\cdot g\right)g$ (124)

This constitutive theory is based on integrity, hence uses complete basis of the space of tensor $q$ . Constitutive theory (124) is cubic in $g$ . ${\kappa }_1$ and ${\kappa }_2$ are material coefficients (thermal conductivities) defined in the known configuration $\underset{\_}{\Omega }$ . The material coefficients are functions of ${{}^q\underset{˜}{I}|}_{\underset{\_}{\Omega }}$ i.e. ${\left(g\cdot g\right)|}_{\underset{\_}{\Omega }}$ and ${\theta |}_{\underset{\_}{\Omega }}$ . From (124) we can obtain a constitutive theory for $q$ that is linear in $g$ by neglecting product terms of $g$ .

$q=-{{\kappa }_1|}_{\underset{\_}{\Omega }}g$ (125)

5. Complete Mathematical Model

The partial differential equations constituting complete mathematical model consists of balance of linear momenta, energy equation and constitutive theories (complete linear theory) for contravariant second Piola-Kirchhoff stress tensor (119) and the heat vector. Conservation of mass and entropy inequality are also part of the mathematical model, but do not provide additional equations, hence are not included in the following:

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

${\rho }_0\frac{\partial e}{\partial t}+\nabla \cdot q-{\sigma }^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}=0$ (127)

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

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

$\begin{array}{l}{}_d\sigma {}^{\left[m\right]}+\lambda {}_d\sigma {}^{\left[m\right]}+{\lambda }^{0}tr\left({}_d\sigma {}^{\left[m\right]}\right)I+{\displaystyle \underset{J=1}{\overset{m-1}{\sum }}}\left({\lambda }_j^1\left({}_d\sigma {}^{\left[j\right]}\right)+{\lambda }_j^{2}tr\left({}_d\sigma {}^{\left[j\right]}\right)I\right)\\ =2\mu {\epsilon }_{\left[0\right]}+\lambda \left(tr{\epsilon }_{\left[0\right]}\right)I+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(2{\eta }_i{\epsilon }_{\left[i\right]}+{\kappa }_{i}tr\left({\epsilon }_{\left[i\right]}\right)I\right)\end{array}$ (130)

$q=-{{\kappa }_1|}_{\underset{\_}{\Omega }}g-{{\kappa }_2|}_{\underset{\_}{\Omega }}\left(g\cdot g\right)g$ (131)

Remarks

1) A reduced form of (126) - (131) in ${ℝ}^2$ constitutes the mathematical model for 2D beams and plane stress, plane strain deformation physics.

2) Equations (126) - (131) is the mathematical model of nonlinear dynamic deformation physics of viscoelastic plates, shells and 3D beams with memory i.e., in ${ℝ}^3$ .

6. Methods of Solutions for IVPs Described by (126) - (131)

First we note that (126) - (131) are a system of nonlinear PDEs in $u\mathrm{,}\theta \mathrm{,}{}_d\sigma {}^{\left[0\right]}$ and $q$ in which all dependent variables are simultaneously dependent on space and time. Thus space-time coupled finite element methods of solution are ideally suited for obtaining their solution. Space-time decoupled finite element methods can also be used for obtaining solutions of the IVPs defined by (126) - (131). We make some remarks:

1) In case of thermoelastic and thermoviscoelastic (no memory) solid matter with dissipation and elasticity, the constitutive theories for ${}_d\sigma {}^{\left[0\right]}$ are not differential equations in time, hence ${}_d\sigma {}^{\left[0\right]}$ from the constitutive theory can be substituted in the balance of linear momenta as done by Surana et al. in references [13] [14] . The resulting balance of linear momenta is nonlinear PDEs in displacement (for isothermal physics). Application of space-time decoupled finite element method of approximation using GM/WF for a spatial discretization [4] yields a system of nonlinear familiar ODEs in degrees of freedom $\left\{\delta \right\}$ and their time derivatives $\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ . This form permits use of many standard solution methods for ODEs in time. Details of the mathematical model and the space-time decoupled solution method for IVPs and BVPs have been presented by Surana et al. [13] [14] .

2) In the present mathematical model it is not possible to substitute ${}_d\sigma {}^{\left[0\right]}$ from the constitutive theory into balance of linear momenta to obtain three PDEs purely in terms of spatial and time derivatives of displacement. This is rather a serious restriction as it precludes use of standard methods of solutions applicable to BLM in displacements.

3) We note that ${}_d\sigma {}^{\left[0\right]}$ in the constitutive theories are functions of space and time. In phenomenological constitutive models in 1D [9] ${}_d\sigma {}^{\left[0\right]}$ reduces to ${}_d\sigma {}_{11}^{\left[0\right]}$ and we have a single ODE in time for ${}_d\sigma {}_{11}^{\left[0\right]}$ . In this case time integration is possible as done in ref [1] [2] to obtain memory modulus, but the result of time integration is an integral form and not differential or algebraic form. Serious shortcomings of this approach are: 1) not applicable to continuous matter 2) and secondly, it cannot be extended to obtain a constitutive theory for continuous matter.

4) Thus, we conclude that we must seek solution methods that can be directly applied to PDEs (126) - (131). Of course substitutions of (128) in (126), (127) and (131) in (127) are made to obtain PDEs in desired dependent variables and to reduce the number of dependent variables.

6.1. Finite Element Formulations

As mentioned earlier, (126) - (131) is a system of nonlinear PDEs in dependent variables $u\mathrm{,}\theta \mathrm{,}{}_d\sigma {}^{\left[0\right]}$ and $q$ in which all dependent variables exhibit simultaneous dependence on space coordinates $x$ and time t. Thus, space-time coupled finite element processes in which simultaneous dependence of evolution on space coordinates $x$ and time t is maintained are ideally suited [4] . These methods are highly meritorious for IVPs in ${ℝ}^1$ and ${ℝ}^2$ (in space) but become rather demanding for IVPs in ${ℝ}^3$ . Thus, when addressing methods of solutions of the IVPs described in (126) - (131) i.e., methods of approximations such as finite element method, we must consider space time decoupled methods also in which simultaneous dependence of evolution in $x$ and t is not preserved. This is necessitated to avoid computational complexities of space-time coupled methods in ${ℝ}^3$ . We present details of both methods in the following sections.

Space-Time Coupled Finite Element Methods of Approximation

In deriving finite element formulations of (126) - (131), it is advantageous to substitute constitutive theory for heat flux $q$ in the energy Equation (127). Also, we note the following:

We substitute ${\sigma }^{\left[0\right]}={}_e\sigma {}^{\left[0\right]}+{}_d\sigma {}^{\left[0\right]}$ in (127). Thus, we can write BLM as:

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

or

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

Let

$\left[J\right]\left[{}_d\sigma {}^{\left[0\right]}\right]=\left[J^{\mathrm{<mo>\; *\; </mo>}}\right]$ (134)

Substituting (134) in (133) and defining $\left\{F^V\right\}$ (due to volumetric change), we obtain:

$\left\{F^V\right\}=\left[\left[J\right]\left[{}_e\sigma {}^{\left[0\right]}\right]\right]\{\nabla \}$ (135)

${\rho }_0\frac{{\partial }^2{u}_i}{\partial t^2}-{\rho }_0{F}_b^i-\frac{\partial J_{ij}^{\mathrm{<mo>\; *\; </mo>}}}{\partial x_j}=F_i^V$ (136)

for simplicity of illustration we choose $m=1$ in (130) and assume $e=c_v\theta$ ; $c_v$ is constant specific heat.

In finite element processes, matrix and vector forms of equations are readily usable, hence we rewrite constitutive theory (130) in Voigt’s notation.

$\left\{{}_d\sigma {}^{\left[0\right]}\right\}+\lambda \frac{\partial }{\partial t}\left\{{}_d\sigma {}^{\left[0\right]}\right\}=\left[\underset{˜}{D}\right]\left\{{\epsilon }_{\left[0\right]}\right\}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[{\underset{˜}{C}}_i\right]\left\{{\epsilon }_{\left[i\right]}\right\}$ (137)

Coefficients of $\left[D\right]$ are given by:

$\left[\underset{˜}{D}\right]=\left[\begin{array}{cccccc}2\underset{˜}{\mu }+\underset{˜}{\lambda }& 2\underset{˜}{\mu }& 2\underset{˜}{\mu }& 0& 0& 0\\ \underset{˜}{\lambda }& 2\underset{˜}{\mu }+\underset{˜}{\lambda }& 2\underset{˜}{\mu }& 0& 0& 0\\ \underset{˜}{\lambda }& \lambda & 2\underset{˜}{\mu }+\underset{˜}{\lambda }& 0& 0& 0\\ 0& 0& 0& 2\underset{˜}{\mu }& 0& 0\\ 0& 0& 0& 0& 2\underset{˜}{\mu }& 0\\ 0& 0& 0& 0& 0& 2\underset{˜}{\mu }\end{array}\right]$ (138)

Coefficients of $\left[{\underset{˜}{C}}_i\right]$ can be obtained by replacing $\underset{˜}{\eta }\mathrm{,}\underset{˜}{\lambda }$ in (138) by $2{\eta }_i\mathrm{,}{\kappa }_i$ .

${\left\{{}_d\sigma {}^{\left[0\right]}\right\}}^{\text{T}}=\left[{}_d\sigma {}_{11}^{\left[0\right]}\mathrm{,}{}_d\sigma {}_{22}^{\left[0\right]}\mathrm{,}{}_d\sigma {}_{33}^{\left[0\right]}\mathrm{,}{}_d\sigma {}_{23}^{\left[0\right]}\mathrm{,}{}_d\sigma {}_{31}^{\left[0\right]}\mathrm{,}{}_d\sigma {}_{12}^{\left[0\right]}\right]$ (139)

Elements of $\left[{\epsilon }_{\left[0\right]}\right]$ and $\left[{\epsilon }_{\left[i\right]}\right]$ in $\left\{{\epsilon }_{\left[0\right]}\right\}$ and $\left\{{\epsilon }_{\left[i\right]}\right\}$ are arranged in the same order as those of $\left\{{}_d\sigma {}^{\left[0\right]}\right\}$ in $\left\{{}_d\sigma {}^{\left[0\right]}\right\}$ .

It is convenient to use new notation for the components of $\left\{{}_d\sigma {}^{\left[0\right]}\right\},\left\{{\epsilon }_{\left[0\right]}\right\}$ and $\left\{{\epsilon }_{\left[i\right]}\right\}$ .

$\left[{}_d\sigma {}^{\left[0\right]}\right]=\left[{\sigma }^1{\sigma }^2{\sigma }^3{\sigma }^4{\sigma }^5{\sigma }^6\right]$ (140)

$\left[{\epsilon }_{\left[0\right]}\right]=\left[{\epsilon }^1{\epsilon }^2{\epsilon }^3{\epsilon }^4{\epsilon }^5{\epsilon }^6\right]$ (141)

$\left[{\epsilon }_{\left[i\right]}\right]=\left[{\epsilon }_i^1{\epsilon }_i^2{\epsilon }_i^3{\epsilon }_i^4{\epsilon }_i^5{\epsilon }_i^6\right]$ (142)

where i in (142) refers to time derivative of order i. Using (140) - (142) in (137), substituting (131) in (127) and using (128) we have the following mathematical model.

${\rho }_0\frac{{\partial }^2{u}_i}{\partial t^2}-{\rho }_0{F}_i^b-\frac{\partial J_{ij}^{\mathrm{<mo>\; *\; </mo>}}}{\partial x_j}=F_i^V$ (143)

${\rho }_0{c}_v\frac{\partial \theta }{\partial t}-{\kappa }_1{\nabla }^2\theta -{}_e\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}-{}_d\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}=0$ (144)

$\left\{\begin{array}{c}{\sigma }^1\\ {\sigma }^2\\ {\sigma }^3\\ {\sigma }^4\\ {\sigma }^5\\ {\sigma }^6\end{array}\right\}+\lambda \frac{\partial }{\partial t}\left\{\begin{array}{c}{\sigma }^1\\ {\sigma }^2\\ {\sigma }^3\\ {\sigma }^4\\ {\sigma }^5\\ {\sigma }^6\end{array}\right\}-\left[\underset{˜}{D}\right]\left\{\begin{array}{c}{\epsilon }^1\\ {\epsilon }^2\\ {\epsilon }^3\\ {\epsilon }^4\\ {\epsilon }^5\\ {\epsilon }^6\end{array}\right\}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[{\underset{˜}{C}}_i\right]\left\{\begin{array}{c}{\epsilon }_i^1\\ {\epsilon }_i^2\\ {\epsilon }_i^3\\ {\epsilon }_i^4\\ {\epsilon }_i^5\\ {\epsilon }_i^6\end{array}\right\}=0$ (145)

These are the final PDEs describing the IVPs related to beams, plates and shells in ten variables ( $u_i\mathrm{,}{\sigma }_i\mathrm{,}\theta$ ). These equations must be nondimensionalized before using them in computations. We do this when we present numerical studies. For deriving finite element formulations, (143) - (145) suffice as the nondimensionalization only introduces dimensionless parameters.

6.2. Space-Time Coupled Finite Element Method of Approximation

If we write (143) - (145) as:

$A\varphi -f=0\forall \left(x\mathrm{,}t\right)\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t$ (146)

in which $A$ is the space time differential operator and $\varphi$ is a vector of dependent variables then it is straight forward to see that $A$ is not linear, and since $A$ is space time differential operator, it cannot be symmetric. Thus, out of all space-time coupled finite element formulations [4] (space-time Galerkin method (STGM), space-time Galerkin method with weak form (STGM/WF), space-time Petrov Galerkin method (STPGM), space-time weighted residual method (STWRM) and space-time finite element method based on space-time residual functional (STRF)) only the finite element method based on STRF yields space-time integral form that is space-time variationally consistent [4] , hence the resulting computational processes are unconditionally stable during the entire evolution [4] . We present details of the STRF method in ${ℝ}^3$ (can be easily reduced to ${ℝ}^1$ and ${ℝ}^2$ ).

Let ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ be the discretization of the space-time domain ${\Omega }_{xt}^{\left(n\right)}$ of the nth space-time strip ${\bar{\Omega }}_x\times \left[t_n\mathrm{,}{t}_{n+1}\right]$ (see ref [4] ), then:

${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T={\displaystyle \underset{e}{\cup }{\bar{\Omega }}_{xt}^e}\text{and}{\bar{\Omega }}_{xt}^T={\displaystyle \underset{i}{\cup }{\left({\bar{\Omega }}_{xt}^{\left(i\right)}\right)}^T}$ (147)

in which ${\bar{\Omega }}_{xt}^e$ is a p-version hierarchical space-time finite element with higher order global differentiability in space and time.

Let ${\left\{\varphi \right\}}^{\text{T}}=\left[u_1{u}_2{u}_3{\sigma }^1{\sigma }^2{\sigma }^3{\sigma }^4{\sigma }^5{\sigma }^6\theta \right]$ be the order of dependent variables.

Let ${\left(u_i\left(x\mathrm{,}t\right)\right)}_h^e\mathrm{,}{\left({\sigma }^j\left(x\mathrm{,}t\right)\right)}_h^e$ and ${\theta }_h^e\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,3\; <mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,6}$ be the local approximations of $u_i\mathrm{,}{\sigma }^j$ and $\theta$ over a space time element ${\bar{\Omega }}_{xt}^e$ , then we have (using equal order equal degree approximation) the following:

$\begin{array}{l}{\left(u_i\left(x\mathrm{,}t\right)\right)}_h^e=\left[N\left(x\mathrm{,}t\right)\right]\left\{{\delta }_{u_i}^e\right\}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,3}\\ {\left({\sigma }^j\left(x\mathrm{,}t\right)\right)}_h^e=\left[N\left(x\mathrm{,}t\right)\right]\left\{{\delta }_{{\sigma }^j}^e\right\}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,6}\\ {\theta }_h^e=\left[N\left(x\mathrm{,}t\right)\right]\left\{{\delta }_{\theta }^e\right\}\\ {\left\{{\varphi }_h^e\right\}}^{\text{T}}=\left[{\left(u_i\right)}_h^e\mathrm{,}{\left({\sigma }^j\right)}_h^e\mathrm{,}{\theta }_h^e\right]\end{array}$ (148)

Substituting (148) in (143) - (145) we obtain residual equation: $E_k^e\mathrm{<mo>\; ;\; </mo>}k=\mathrm{1,2,}\cdots \mathrm{,10}\forall \left(x\mathrm{,}t\right)\in {\bar{\Omega }}_{xt}^e\mathrm{.}$

The approximation $\left\{{\varphi }_h\right\}$ over ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ is given by:

$\left\{{\varphi }_h\right\}={\displaystyle \underset{e}{\cup }\left\{{\varphi }_h^e\right\}}$ (149)

Substituting $\left\{{\varphi }_h\right\}$ in (143) - (145) yield residual equations $E_k\mathrm{<mo>\; ;\; </mo>}k=\mathrm{1,2,}\cdots \mathrm{,10}$ for the discretization ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ and we can construct STRF $I\left({\varphi }_h\right)$ over ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ .

$I\left({\varphi }_h\right)={\displaystyle \underset{k=1}{\overset{10}{\sum }}}{\left(E_k,E_k\right)}_{{\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T}={\displaystyle \underset{e}{\sum }}\left({\displaystyle \underset{k=1}{\overset{10}{\sum }}}{\left(E_k^e,E_k^e\right)}_{{\bar{\Omega }}_{xt}^e}\right)$ (150)

If the space-time functional $I\left(\mathrm{.}\right)$ is differentiable in its arguments, then $\delta I\left({\varphi }_h\right)$ is unique and $\delta \left({\varphi }_h\right)=0$ is a necessary condition for an extremum of $I\left({\varphi }_h\right)$ defined by (150).

$\begin{array}{c}\delta I\left({\varphi }_h\right)={\displaystyle \underset{k=1}{\overset{10}{\sum }}}2{\left(E_k,\delta E_k\right)}_{{\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T}={\displaystyle \underset{e}{\sum }}2\left({\displaystyle \underset{k=1}{\overset{10}{\sum }}}{\left(E_k^e,\delta E_k^e\right)}_{{\bar{\Omega }}_{xt}^e}\right)\\ ={\displaystyle \underset{e}{\sum }}\left\{g^e\right\}=\left\{g\right\}=0\end{array}$ (151)

If

$\left\{\delta \right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }^e\right\}};{\left\{{\delta }^e\right\}}^{\text{T}}=\left[{\left\{{\delta }_{u_i}^e\right\}}^{\text{T}},{\left\{{\delta }_{{\sigma }^j}^e\right\}}^{\text{T}},{\left\{{\delta }_{\theta }^e\right\}}^{\text{T}}\right]$ (152)

then $\left\{g\right\}$ in (151) is a nonlinear function of $\left\{\delta \right\}$ . We satisfy $\left\{g\left(\left\{\delta \right\}\right)\right\}=0$ iteratively by using Newton’s linear method with line search [4] . Thus, if ${\left\{\delta \right\}}_0$ is assumed starting solution, then improved solution $\left\{\delta \right\}$ is given by:

$\left\{\delta \right\}={\left\{\delta \right\}}_0+{\alpha }^{\mathrm{<mo>\; *\; </mo>}}\left\{\Delta \delta \right\}$ (153)

in which $\left\{\Delta \delta \right\}$ is calculated using:

$\left\{\Delta \delta \right\}=-{\left[{\delta }^{2}I\left({\varphi }_h\right)\right]}_{{\left\{\delta \right\}}_0}^{-1}\left\{g\left({\left\{\delta \right\}}_0\right)\right\}$ (154)

${\delta }^{2}I\left({\varphi }_h\right)$ is approximated by [4] .

${\delta }^{2}I\left({\varphi }_h\right)\simeq -{\displaystyle \underset{k=1}{\overset{10}{\sum }}}2{\left(\delta E_k,\delta E_K\right)}_{{\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T}={\displaystyle \underset{e}{\sum }}2\left({\displaystyle \underset{k=1}{\overset{10}{\sum }}}{\left(\delta E_k^e,\delta E_k^e\right)}_{{\bar{\Omega }}_{xt}^e}\right)$ .

${\alpha }^{\mathrm{<mo>\; *\; </mo>}}$ is determined using line search [4] .

Using new $\left\{\delta \right\}$ , we check if $\left\{g\left(\left\{\delta \right\}\right)\right\}\le \Delta$ holds. If it does, then $\left\{\delta \right\}$ is the converged solution, if not then we set ${\left\{\delta \right\}}_0=\left\{\delta \right\}$ and repeat the iteration process till convergence. $\Delta$ is a preset tolerance for computed zero. Computations begin with first space-time strip and then time marching using subsequent space-time strip till the final time is reached [4] .

Remarks

1) This method can be easily reduced to IVPs in ${ℝ}^1$ and ${ℝ}^2$ .

2) Newton’s linear method requires that starting or initial solution ${\left\{\delta \right\}}_0$ must be in close proximity of the true solution $\left\{\delta \right\}$ . Generally, physics of the IVPs provide means for appropriate choice of ${\left\{\delta \right\}}_0$ .

Space-Time Decoupled Finite Element Processes

As mentioned earlier, in this mathematical model, ${}_d\sigma {}^{\left[0\right]}$ from the constitutive theory cannot be substituted in the balance of linear momenta, thus standard form of ODEs resulting from space-time decoupled methods using GM/WF for a discretization giving mass, stiffness and damping matrices is not possible. Thus, we are left with no choice but construct a space-time decoupled finite element formulation of (143) - (145) as they are.

Let ${\bar{\Omega }}_{xt}^T={\displaystyle \underset{e}{\cup }{\bar{\Omega }}_x^e}$ be discretization of spatial domain ${\bar{\Omega }}_x$ in which ${\bar{\Omega }}_x^e$ is a typical p-version hierarchical finite element “e” in space with higher order global differentiability. We consider following space-time decoupled local approximation over an element e over ( ${\bar{\Omega }}_x^e$ ).

Local approximation over an element e over ( ${\bar{\Omega }}_x^e$ ):

$\begin{array}{l}{\left(u_i\left(x,t\right)\right)}_h^e=\left[N\left(x\right)\right]\left\{{\delta }_{u_i}^e\left(t\right)\right\};i=1,2,3\\ {\left({\sigma }^j\left(x,t\right)\right)}_h^e=\left[N\left(x\right)\right]\left\{{\delta }_{{\sigma }^j}^e\left(t\right)\right\};j=1,2,\cdots ,6\\ {\left(\theta \left(x,t\right)\right)}_h^e=\left[N\left(x\right)\right]\left\{{\delta }_{\theta }^e\left(t\right)\right\}\forall x,t\in {\bar{\Omega }}_{xt}^e\end{array}$ (155)

Let ${\beta }_1\mathrm{,}{\beta }_2\mathrm{,}{\beta }_3$ be test functions such that ${\beta }_i=\delta {\left(u_i\right)}_h^e$ and ${\beta }_4=\delta {\left({\theta }_h\right)}^e$ , ${\beta }_{j+4}=\delta {\left({\sigma }^j\right)}_h^e\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,6}$ . If we denote (143) - (145) as $A_i\left(\mathrm{.}\right)-f_i=0\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,3}$ , $A_4(.)-f_4=0$ and $A_k\left(\mathrm{.}\right)-f_k=0\mathrm{<mo>\; ;\; </mo>}k=\mathrm{5,}\cdots \mathrm{,10}$ , then for the discretization ${\bar{\Omega }}_x^T$ of ${\bar{\Omega }}_x$ , we can write the following using fundamental lemma of the calculus of variations.

${\left(A_i\left(\mathrm{.}\right)-f_i\mathrm{,}{\beta }_i\right)}_{{\bar{\Omega }}_x^T}=0\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,10}$ (156)

Since (156) are functionals, we can write:

${\left(A_i\left(\mathrm{.}\right)-f_i\mathrm{,}{\beta }_i\right)}_{{\bar{\Omega }}_x^T}={\displaystyle \underset{e}{\sum }}{\left(A_i^e\left(\mathrm{.}\right)-f_i^e\mathrm{,}{\beta }_i\right)}_{{\bar{\Omega }}_x^e}=0\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,10}$ (157)

$A_i^e$ and $f_i^e$ signify that they are over ${\bar{\Omega }}_x^e$ . Each ${\beta }_i=N_j\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,10}\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}n$ .

Consider ${\left(A_i^e\left(\mathrm{.}\right)-f_i^e\mathrm{,}{\beta }_i\right)}_{{\bar{\Omega }}_x^e}$ and substitute $A^e\left(\mathrm{.}\right)$ and $f_i^e$ from (143) - (145). First consider:

$\begin{array}{l}{\left(A_i^e\left(\mathrm{.}\right)-f_i^e\mathrm{,}{\beta }_i\right)}_{{\bar{\Omega }}_x^e}={\left({\rho }_0\frac{{\partial }^2{\left(u_i\right)}_h^e}{\partial t^2}-{\rho }_0{F}_i^b-\frac{\partial J_{ij}^{\mathrm{<mo>\; *\; </mo>}}}{\partial x_j}\mathrm{,}{\beta }_i\right)}_{{\bar{\Omega }}_x^e}\\ ={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left({\rho }_0\frac{{\partial }^2{\left(u_i\right)}_h^e}{\partial t^2}{\beta }_i-{\rho }_0{F}_i^b{\beta }_i-\frac{\partial J_{ij}^{\mathrm{<mo>\; *\; </mo>}}}{\partial x_j}\right)\text{d}x\\ \text{nosumon}i\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,3}\end{array}$ (158)

In the last term we perform integration by parts once to transfer differentiation to ${\beta }_i$ . We also note that $J^{*}$ is a function of gradients of $u_i$ and ${\sigma }^j$ .

$\begin{array}{c}{\left(A_i^e\left(\mathrm{.}\right)-f_i^e\mathrm{,}{\beta }_i\right)}_{{\bar{\Omega }}_x^e}={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left({\rho }_0\frac{{\partial }^2{\left(u_i\right)}_h^e}{\partial t^2}{\beta }_i-{\rho }_0{F}_i^b{\beta }_i-\frac{\partial {\beta }_i}{\partial x_j}{J}_{ij}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }_i\right)\text{d}x\\ +{\displaystyle {\oint }_{{\Gamma }^e}}{J}_{ij}^{\mathrm{<mo>\; *\; </mo>}}{\beta }_i{n}_j\text{d}\Gamma \text{nosumover}i\text{and}j\end{array}$ (159)

using local approximation for ${\left(u_i\right)}_h^e$ from (155), (159) can be written as:

${\left(A_i^e\left(\mathrm{.}\right)-f_i^e\mathrm{,}{\beta }_i\right)}_{{\bar{\Omega }}_x^e}=\left[M_u^e\right]\left\{{\ddot{\delta }}_u^e\right\}+\left[K_1^e\right]\left\{{\delta }_u^e\right\}+\left[K_2^e\right]\left\{{\delta }_{\sigma }^e\right\}-\left\{f_{u_i}^e\right\}-\left\{P_{u_i}^e\right\}$ (160)

$\left\{f_{u_i}^e\right\}$ is due to $F_i^b$ and $\left\{P_{u_i}^e\right\}$ is due to concomittant in (159).

Consider integral form (144) based on fundamental lemma, we can write the following:

$\begin{array}{l}\left(A_4^e\left(\mathrm{.}\right)-f_4^e\mathrm{,}{\beta }_4\right)\\ ={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left({\rho }_0{c}_v\frac{\partial \theta }{\partial t}{\beta }_4-{\kappa }_1{\nabla }^2\theta {\beta }_4-\left({}_e\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\right){\beta }_4-\left({}_d\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\right){\beta }_4\right)\text{d}x\end{array}$ (161)

In the second term we perform integration by parts once to transfer one order differentiation to ${\beta }_4$ .

$\begin{array}{c}\left(A_4^e\left(\mathrm{.}\right)-f_4^4\mathrm{,}{\beta }_4\right)={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}({\rho }_0{c}_v\frac{\partial \theta }{\partial t}{\beta }_4+{\kappa }_1\nabla {\beta }_4\cdot \nabla {\theta }_h^e-\left({}_e\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\right){\beta }_4\\ -\left({}_d\sigma {}^{\left[0\right]}:{\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\right){\beta }_4)\text{d}x-{\displaystyle \underset{{\Gamma }^e}{\oint }}\kappa \frac{\partial {\theta }_h^e}{\partial n}\text{d}\Gamma \end{array}$ (162)

or

$\left(A_4^e\left(\mathrm{.}\right)-f_4^4\mathrm{,}{\beta }_4\right)=\left[C^e\right]\left\{{\stackrel{\dot{}}{\delta }}_{\theta }^e\right\}+\left[H^e\right]\left\{{\delta }_{\theta }^e\right\}-\left[K_3^e\right]\left\{{\stackrel{\dot{}}{\delta }}_u^e\right\}-\left[K_4^e\right]\left\{{\delta }_{{\sigma }_j}^e\right\}-\left\{P_{\theta }^e\right\}$ (163)

in which $\left[C^e\right]$ and $\left[H^e\right]$ are symmetric matrices. The third and fourth term in (163) result due to trace terms in (162) and $\left\{P_{\theta }^e\right\}$ are secondary variables due to concomitant in (162).

Consider integral form of (145) based on fundamental lemma.

$\begin{array}{l}{\left(A_k^e\left(\mathrm{.}\right)-f_k^e\mathrm{,}{\beta }_{k+4}\right)}_{{\bar{\Omega }}_x^e}\\ ={\left({\left({\sigma }_k\right)}_h^e\mathrm{,}{\beta }_{k+4}\right)}_{{\bar{\Omega }}_x^e}+\lambda {\left(\frac{\partial {\left({\sigma }_k\right)}_h^e}{\partial t}\mathrm{,}{\beta }_{k+4}\right)}_{{\bar{\Omega }}_x^e}-{\left({\underset{˜}{D}}_{kj}{\left\{{\epsilon }^j\right\}}_h^e\mathrm{,}{\beta }_{k+4}\right)}_{{\bar{\Omega }}_x^e}\\ -{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\left({\left({\underset{˜}{c}}_i\right)}_{kp}{\left\{{\epsilon }_{\mathrm{,}i}^p\right\}}_h^e\mathrm{,}{\beta }_{k+4}\right)}_{{\bar{\Omega }}_x^e}\mathrm{<mo>\; ;\; </mo>}k=\mathrm{1,2,}\cdots \mathrm{,6}\end{array}$ (164)

Substituting local approximations in (164).

${\left(A_k^e\left(\mathrm{.}\right)-f_k^e\mathrm{,}{\beta }_{k+4}\right)}_{{\bar{\Omega }}_x^e}=\left[K_5^e\right]\left\{{\delta }_{\sigma }^e\right\}+\left[K_6^e\right]\left\{{\stackrel{\dot{}}{\delta }}_{\sigma }^e\right\}-\left[K^e\right]\left\{{\delta }_u^e\right\}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[C_i^e\right]\left\{{\delta }_{\mathrm{,}i}^u\right\}-\left\{P_{\sigma }^e\right\}$ (165)

where $\left\{{\delta }_{\sigma }\right\}$ are all degrees of freedom due to ${\left({\sigma }_j\right)}_h^e\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,6}$ .

Assembly of the element expression for this discretization (based on (160), (165) and (164)):

$\left[M_u^e\right]\left\{{\ddot{\delta }}_u^e\right\}+\left[K_1^e\right]\left\{{\delta }_u^e\right\}+\left[K_2^e\right]\left\{{\delta }_{\sigma }^e\right\}-\left\{f_{u_i}^e\right\}-\left\{P_{u_i}^e\right\}$ (166)

$\left\{{\delta }_{\sigma }^e\right\}+\left[K_6^e\right]\left\{{\stackrel{\dot{}}{\delta }}_{\sigma }^e\right\}-\left[K^e\right]\left\{{\delta }_u^e\right\}-{\displaystyle \underset{i\mathrm{<mo>\; =\; </mo>1}}{\overset{n}{\sum }}}\left[C_i^e\right]\left\{{\delta }_{\mathrm{,}i}^u\right\}-\left\{P_{\sigma }^e\right\}$ (167)

$\left\{{\delta }_{\theta }^e\right\}+\left[H^e\right]\left\{{\ddot{\delta }}_{\theta }^e\right\}-\left[K_3^e\right]\left\{{\stackrel{\dot{}}{\delta }}_u^e\right\}-\left[K_4^e\right]\left\{{\delta }_{{\sigma }_j}^e\right\}-\left\{P_{\theta }^e\right\}$ (168)

If we group ${\left\{\delta \right\}}^{\text{T}}=\left[{\left\{{\delta }_u\right\}}^{\text{T}}{\left\{{\delta }_{\sigma }\right\}}^{\text{T}}{\left\{{\delta }_{\theta }\right\}}^{\text{T}}\right]$ and consider only $n=1$ and define $\left\{{\delta }_{\mathrm{,1}}^u\right\}=\left\{{\stackrel{\dot{}}{\delta }}_u\right\}$ , then we can write (166) - (168) as:

$\begin{array}{l}\left[\begin{array}{cc}\left[M_u\right]& \left[0\right]\\ \left[0\right]& \left[0\right]\\ \left[0\right]& \left[0\right]\end{array}\right]\left\{\begin{array}{c}\left\{{\ddot{\delta }}_u\right\}\\ \left\{{\ddot{\delta }}_{\sigma }\right\}\\ \left\{{\ddot{\delta }}_{\theta }\right\}\end{array}\right\}+\left[\begin{array}{ccc}\left[0\right]& \left[0\right]& \left[0\right]\\ \left[C_1\right]& \left[K_6\right]& \left[0\right]\\ \left[K_3\right]& \left[0\right]& \left[C\right]\end{array}\right]\left\{\begin{array}{c}\left\{{\stackrel{\dot{}}{\delta }}_u\right\}\\ \left\{{\stackrel{\dot{}}{\delta }}_{\sigma }\right\}\\ \left\{{\stackrel{\dot{}}{\delta }}_{\theta }\right\}\end{array}\right\}\\ +\left[\begin{array}{ccc}\left[K_1\right]& \left[K_2\right]& \left[0\right]\\ \left[K\right]& \left[K_5\right]& \left[0\right]\\ \left[0\right]& \left[K_4\right]& \left[H\right]\end{array}\right]\left\{\begin{array}{c}\left\{{\delta }_u\right\}\\ \left\{{\delta }_{\sigma }\right\}\\ \left\{{\delta }_{\theta }\right\}\end{array}\right\}=\left\{\begin{array}{c}\left\{f_u\right\}+\left\{P_u\right\}\\ \left\{P_{\sigma }\right\}\\ \left\{P_{\theta }\right\}\end{array}\right\}\end{array}$ (169)

This is a system of nonlinear second order system of ODEs in time.

These can be integrated using any of the time integration methods. Based on our previous work [14] , Newmark linear acceleration method with Newton’s linear method is well suited. In these equations, initial conditions and initial solutions are always critical and important to determine to ensure that time integration methods yield physical evolution.

Initial conditions

For a system at rest and at NTP, we have the following initial conditions

$\left\{{\delta }_u\right\}=0,\left\{{\stackrel{\dot{}}{\delta }}_u\right\}=0$

If ${\theta |}_{t_0}={\theta }_0$ is known temperature, then $\left\{{\delta }_{\theta }\right\}$ are known.

For a system at rest, ${\sigma }^i$ are zero i.e., $\left\{{\delta }_{\sigma }\right\}$ are zero. Secondary variables for unknown dofs are known in (169)

In any time integration method for nonlinear ODEs in time we begin with a known solution. Assume $\left\{\delta \right\}={\left\{\delta \right\}}_0=0$ , then all matrices in (169) can be computed.

1) Since ${\left\{{\delta }_u\right\}}_{t_0=0}$ and ${\left\{{\delta }_{\sigma }\right\}}_{t_0=0}$ are known ${\left\{\ddot{\delta }\right\}}_{t_0=0}$ can be calculated using assembled form of (166).

$\left\{{\ddot{\delta }}_u\right\}={\left[M_u\right]}^{-1}\left(-\left[K_1\right]\left\{{\delta }_u\right\}-\left[K_2\right]\left\{{\delta }_{\sigma }\right\}+\left\{f_{u_i}\right\}+\left\{P_{u_i}\right\}\right)$ (170)

2) From assembled Equations (168), since ${\left\{{\delta }_u\right\}}_{t_0=0},{\left\{{\delta }_{\theta }\right\}}_{t_0=0}$ and ${\left\{{\stackrel{\dot{}}{\delta }}_u\right\}}_{t_0=0}$ are known, ${\left\{{\stackrel{\dot{}}{\delta }}_{\theta }\right\}}_{t_0=0}$ can be calculated using assembled of (168):

$\left\{{\delta }_{\theta }\right\}={\left[C\right]}^{-1}\left(-\left[H\right]\left\{{\delta }_{\theta }\right\}+\left[K_3\right]\left\{{\stackrel{\dot{}}{\delta }}_u\right\}+\left[K_4\right]\left\{{\delta }_{{\sigma }_j}\right\}+\left\{P_{\theta }\right\}\right)$ (171)

${\left\{{\ddot{\delta }}_{\theta }\right\}}_{t_0=0}$ is zero at $t_0=0$ as it is not part of the ODEs.

3) Knowing ${\left\{{\delta }_u\right\}}_{t_0=0},{\left\{{\stackrel{\dot{}}{\delta }}_u\right\}}_{t_0=0},{\left\{{\delta }_{\sigma }\right\}}_{t_0=0}$ (with $n=1$ ), ${\left\{{\stackrel{\dot{}}{\delta }}_{\sigma }\right\}}_{t_0=0}$ can be determined using assembled form of Equations (167).

${\left\{{\stackrel{\dot{}}{\delta }}_{\sigma }\right\}}_{t_0=0}={\left[K_6\right]}^{-1}\left(\left[K\right]{\left\{{\delta }_u\right\}}_{t_0=0}+\left[C_1\right]{\left\{{\stackrel{\dot{}}{\delta }}_u\right\}}_{t_0=0}-\left[K_5\right]{\left\{{\delta }_{\sigma }\right\}}_{t_0=0}-{\left\{P_{\sigma }\right\}}_{t_0=0}\right)$ (172)

${\left\{{\ddot{\delta }}_{\sigma }\right\}}_{t_0=0}$ is zero as $t_0=0$ as it is not part of the mathematical model. Then at $t_0=0$ , commencement of the evolution of all degrees of freedom and their first and second time derivatives are known.

Remarks

1) Details of the time integration using Newmark linear acceleration method in conjunction with Newton’s linear method to each integration time step are similar to those presented in ref [14] .

2) The mathematical model needs to be nondimensionalized, these details are presented in the section containing model problem studies.

7. Factors Influencing Non-Linear Dynamics and Dynamic Bifurcation in TVES with Rheology

From the complete mathematical model in Section 5, Equations (126) - (131), we note that the constitutive theory for deviatoric second Piola-Kirchhoff stress tensor is a partial differential equation in deviatoric second Piola-Kirchhoff stress tensor ( ${}_d\sigma {}^{\left[0\right]}$ ) and its time derivatives up to orders m. Thus, in this case ${}_d\sigma {}^{\left[0\right]}$ cannot be substituted in the balance of linear momenta (126) to obtain a system of three PDEs in displacements $u_i$ as done in case of TVES without memory [14] . As a consequence, decoupling of space and time in (126) using fundamental lemma followed by GM/WF for the spatial discretization obviously will not yield the same information regarding stiffness and damping matrices as in TVES without memory [14] as we have seen in the previous section. In the present mathematical model, space-time decoupling must be considered for (126) - (131) preferably without any substitutions (as done here), as it cannot be done with substitutions, hence the main benefit of substitution of deviatoric stress tensor in the BLM cannot be realized in this case. Hence, in the present mathematical model, integrating nonlinear ODEs in time resulting from the space-time decoupling using Newmark linear acceleration method with Newton’s linear method for each time increment will not yield the same tangent matrix $\left[H_T\right]$ as in case of TVE solids without memory [14] . In reference [14] , the factors influencing nonlinear dynamics and dynamic bifurcation in TVES without memory were clearly identified using tangent matrix $\left[H_T\right]$ , controlling the calculation of incremental solution for each time step. In the absence of the $\left[H_T\right]$ precisely in the same form as in ref [14] , can we still identify the factors influencing the nonlinear dynamics in TVES with memory using the mathematical model (126) - (131) and the methods of solutions PDEs in Section 5. We proceed as follows.

We note that in the absence of rheology ( ${}_d\sigma {}^{\left[i\right]}=\mathrm{0,}i=\mathrm{1,2,}\cdots \mathrm{,}m$ ), the constitutive theory for ${}_d\sigma {}^{\left[0\right]}$ for TVES with rheology given by (130) reduces to exactly same constitutive theory as in reference [14] for TVES without memory for which tangent matrix $\left[H_T\right]$ is valid. This suggests that all of the physics of TVES solid without memory is implicitly present in its entirety in the mathematical model (126) - (131). Presence of the physics of rheology may require consideration and/or additional factors over and beyond those due to $\left[H_T\right]$ in ref [14] that influence nonlinear dynamics and dynamic bifurcation in TVES with rheology (polymeric solids). Polymeric solids (rubber like solid materials) consist of solvent (short chain molecules) and a polymer (long chain molecules). In such solids, in the absence of external stimuli, the polymeric solid is in stress free state. The long chain polymer molecules in this stress free state are in relaxed configuration (coilded, observed in electron microscopy). The solvent has its own viscosity and so does the polymer. The viscosity of the polymeric solid is not the sum of the two, but must be determined experimentally.

The mechanical behavior of polymeric solids of course depends upon whether the polymeric solid constitution is dominated by solvent or by the polymer. When dominated by solvent, the polymeric solid roughly behave more like isotropic solid with mild rheology. In polymeric solids dominated by the long chain molecules, the rheological behavior is more pronounced.

Micro-photographs of polymeric solids show that in the relaxed state, the polymer molecules are connected to the neighboring molecules forming colonies of polymer molecules that are interim interconnected wit the neighboring colonies. Of course individual stray molecules are also observed in the solvent. Upon application of external stimuli to polymeric solids, the long chain polymer molecules exhibit complex motion (Brownian motion). 1) The polymer molecules begin to unwind from their relaxes coiled state, 2) in doing so they have to overcome the viscous drag (forces) due to viscosity of the solvent as well as the polymer. This resistance offered by the viscous forces creates dynamic stiffness. Polymer molecules collectively act as one dimensional springs in the direction of the disturbance. This dynamic stiffness obviously requires time dependent deformation of the polymeric solid. This dynamic stiffness is over and beyond the tangent stiffness $\left[K_T\right]$ (ref [14] [20] ). Thus, TVES with rheology have stiffness $\left[K_T\right]$ plus dynamic stiffness compared to TVES without rheology (only $\left[K_T\right]$ ). When compared to $\left[K_T\right]$ , the dynamic stiffness is small, hence only plays a minor role in influencing the total effective stiffness during nonlinear dynamic motion. Overall stiffness ( $\left[K_T\right]$ + dynamic stiffness) of a polymeric solid is lower than the corresponding TVES without rheology. Due to presence of long chain molecules TVES with rheology are easier to deform compared to TVES without memory. With progressively increasing rheology TVES will progressively loose stiffness.

Upon cessation of external stimuli, the polymer molecules try to return to their unstressed coiled state. In doing so, they have to overcome the viscous drag (resistance to motion) due to viscosity of the medium. Thus, upon cessation of external stimuli, the polymer molecules take finite amount of time to achieve their uncoiled stress free state. During this time, the stress in the polymeric solid is not zero but is continuously reducing with elapsed time. This phenomenon is known as stress relaxation, rheology or memory mechanism. Every polymeric solid has a characteristic time constant called “relaxation time”. The physical time it takes for a polymeric solid to achieve stress free state depends upon relaxation time but is not the relaxation time. It has been shown [1] [2] [6] that physical time for complete stress relaxation depends upon relaxation modulus which is a function of relaxation time. Higher relaxation time implies longer clock time for stress relaxation and vice versa. Thus, in TVES with rheology, the stress field at any time during evolution is affected by relaxation phenomenon. For example, a memory deviatoric stress pulse in TVES without memory leaves zero stress behind it during propagation, whereas in case of TVES with memory, the stress behind a propagating stress pulse is not zero. Its magnitude and base depends upon relaxation time and the speed of sound in the solid matter. In summary, we can conclude from this discussion that in TVES with memory the stress field is enhanced (higher) compared to TVES without memory. This, of course holds for compressive stress field as well as tensile stress field and that TVES with rheology have lower stiffness compared to TVES without rheology.

Remarks

A) In nonlinear dynamics, there are two opposing mechanisms: those that result in increase in the stiffness (positive $\left[K_{+}\right]$ ) and those that result in decrease of stiffness (negative $\left[K_{-}\right]$ ) during evolution or frequency response. When the magnitude of $\left[K_{+}\right]$ is greater than the magnitude of $\left[K_{-}\right]$ , we have stable dynamic response of the system. $\left[K_{+}\right]=\left[K_{-}\right]$ is the condition for instability. At this point, a sudden change must occur in the response (amplitude) of the system to restore stable configuration. This of course is dynamic bifurcation.

B) Thus, in nonlinear dynamics of TVES with rheology, we must identify the mechanism that contributes to $\left[K_{+}\right]$ and $\left[K_{-}\right]$ .

C) As discussed earlier, TVES with rheology have all of the mechanisms of stiffness, damping and mass as in TVES without memory [14] , but also contain some additional mechanisms due to presence of rheology due to long chain molecules.

D) In the absence of rheology ( ${}_d\sigma {}^{\left[i\right]}\mathrm{,}i=\mathrm{1,2,}\cdots \mathrm{,}m$ ), the constitutive theory for TVES with rheology (Equation (130)) reduces to the constitutive theory for TVES without memory (ref [14] ), confirming that the physics of TVES without memory is intact in the mathematical model for TVES with rheology. Thus, nonlinear dynamics in TVES with rheology also depends upon tangent matrix $\left[H_T\right]$ (Equation (94) in ref [14] ). Symmetric tangent matrix $\left[H_T\right]$ consists of mass matrix $\left[M\right]$ , tangent stiffness matrix $\left[K_T\right]$ and tangent damping matrix which in turn consists of [14] :

$\left[K_T\right]=\left[K^1\right]+\left[{\underset{˜}{K}}^2\right]+\left[{\underset{˜}{K}}^3\right]+\left[K_{\sigma }\right]$ (173)

$\left[C_T\right]=\left[C^1\right]+\left[{\underset{˜}{C}}^2\right]+\left[{\underset{˜}{C}}^3\right]+\left[C_{\sigma }\right]$ (174)

Thus, the nonlinear dynamic response due to $\left[H_T\right]$ depends upon all of the above.

1) Translational inertial physics ( $\left[M\right]\left\{\ddot{\delta }\right\}$ ) with increasing frequency of excitation results in progressively reducing stiffness.

2) The influence of tangent stiffness $\left[K_T\right]$ on the nonlinear dynamics is through $\left[K^1\right]$ , $\left[{\underset{˜}{K}}^2\right]$ , $\left[{\underset{˜}{K}}^3\right]$ and $\left[K_{\sigma }\right]$ . Coefficients of $\left[K^1\right]$ are constant, coefficients of $\left[{\underset{˜}{K}}^2\right]$ and $\left[{\underset{˜}{K}}^3\right]$ are linear and quadratic functions of $\left\{\delta \right\}$ (total dofs for the discretization). Matrix $\left[K_{\sigma }\right]$ is due to stress field. Increasing stiffness inhibits bifurcation. When the deformation is finite, contributions of $\left[{\underset{˜}{K}}^3\right]$ can be quite significant. Tensile stress field results in positive $\left[K_{\sigma }\right]$ , thus increasing stiffness, hence inhibiting dynamic bifurcation. On the other hand, compressive stress field results in negative $\left[K_{\sigma }\right]$ , thus reducing stiffness, hence promoting existence of dynamic bifurcation.

3) The influence of tangent damping matrix $\left[C_T\right]$ on nonlinear dynamics is through $\left[C^1\right]$ , $\left[{\underset{˜}{C}}^2\right]$ , $\left[{\underset{˜}{C}}^3\right]$ and $\left[C_{\sigma }\right]$ . $\left[C^1\right]$ is a constant coefficient matrix and the coefficients of $\left[{\underset{˜}{C}}^2\right]$ , $\left[{\underset{˜}{C}}^3\right]$ are linear and quadratic functions of $\left\{\delta \right\}$ . Matrix $\left[C_{\sigma }\right]$ in similar to $\left[K_{\sigma }\right]$ but contain stresses due to viscous forces, hence its coefficients are relatively small compared to $\left[K_{\sigma }\right]$ . It influences total dissipation in the dynamic response but its individual contribution is hard to quantify. Dissipation provides resistance to motion, hence increasing dissipation progressively inhibits likelihood of the existence of dynamic bifurcation and vice versa.

4) Thus, due to $\left[H_T\right]$ , when:

a) $\left[K_{\sigma }\right]$ is positive due to tensile stress field, stiffness $\left[K^1\right]$ , $\left[{\underset{˜}{K}}^2\right]$ , $\left[{\underset{˜}{K}}^3\right]$ , dissipation and $\left[K_{\sigma }\right]$ are contributing to $\left[K_{+}\right]$ and only translational inertial physics ( $\left[M\right]\left\{\ddot{\delta }\right\}$ ) contributes to $\left[K_{-}\right]$ .

b) When $\left[K_{\sigma }\right]$ is negative due to compressive stress field; $\left[K^1\right]$ , $\left[{\underset{˜}{K}}^2\right]$ , $\left[{\underset{˜}{K}}^3\right]$ and dissipation are contributing to $\left[K_{+}\right]$ whereas translation inertial physics ( $\left[m\right]\left\{\ddot{\delta }\right\}$ ) and $\left[K_{\sigma }\right]$ contribute to $\left[K_{-}\right]$ .

5) We keep in mind that reduced dissipation, negative $\left[K_{\sigma }\right]$ and progressively increasing influence of translational inertial physics, all contribute to $\left[K_{-}\right]$ . We have shown in ref [14] in the numerical studies for model problem II that lack of dynamic bifurcation for some ${\sigma }_0$ (model problem II) and damping can be overcome by reducing damping. Likewise, lack of dynamic bifurcation for some ${\sigma }_0$ and damping can also be overcome by increasing ${\sigma }_0$ .

6) In TVES with rheology, there is additional physics over and beyond $\left[H_T\right]$ (discussed in item (D)) that must be considered to determine how it influences nonlinear dynamics and dynamic bifurcation. Due to rheology, the additional physics are: a) additional stiffness due to stretching of long chain molecules but overall total lower stiffness compared to TVES without memory, b) the stress relaxation phenomenon.

Surana et al. [20] [21] have shown that additional stiffness in TVES with memory due to stretching of long chain molecules called dynamic stiffness is very small compared to the stiffness due to strain. Authors (in ref [20] ) isolated this physics in wave propagation studies to confirm that this is indeed true. Thus, the total stiffness of the deforming TVES with rheology is not affected appreciably by this dynamic stiffness. The stress relaxation phenomenon on the other hand needs more careful consideration. As discussed earlier, in TVES with rheology, upon cessation of external stimuli, finite amount of time is needed for the material to achieve stress free state. This physical clock time is determined using relaxation modulus related to relaxation time, a material property. TVES with longer relaxation time requires more clock time to relax and vice versa. Thus, in TVES with rheology, there is speed of sound and there is relaxation time or modulus. Thus, depending upon the speed of sound and relaxation time, we could have a situation in which a propagating disturbance leaves behind nonzero stress field, thus altering the stress field behind the propagating disturbance. A consequence of this is that due to change in the stress field, $\left[K_{\sigma }\right]$ , the stiffness due to stress field will change. We consider two possibilities that can exist.

E) If the external stimuli creates a compressive stress field, then the stress relaxation may leave additional compressive (negative) stress field behind the propagating disturbance (depends upon relaxation time) thereby increasing the magnitude of the total compressive stress field which enhances negative $\left[K_{\sigma }\right]$ , hence adds in $\left[K_{-}\right]$ . This additional reduction in stiffness promotes existence of dynamic bifurcation at an earlier (lower) frequency compared to TVES without memory as now we do not require as much reduction in stiffness due to translational inertial physics. In summary, for the same material coefficients and loading, dynamic bifurcation in TVES with memory can occur at a lower frequency compared to TVES without rheology. There could also be instances where TVES without rheology show no dynamic bifurcation for some choices of parameter, but for the same parameters dynamic bifurcation could exist in TVES with memory due to additional loss of stiffness by increased contribution of negative $\left[K_{\sigma }\right]$ . In the axial rod model problem considered in a later section, $\left[K_{\sigma }\right]$ is negative, hence the details presented here are applicable. In this case, we expect dynamic bifurcation for TVES with memory to occur at a frequency smaller than for TVES without memory.

F) When the external stimuli creates a tensile stress field, then the stress relaxation may leave additional tensile (positive) stress field behind the propagating disturbance thereby increasing the magnitude of the total tensile stress field which enhances positive $\left[K_{\sigma }\right]$ , hence aids to $\left[K_{+}\right]$ . This additional increase in stiffness need to be overcome by additional reduction in stiffness due to translational inertial physics, hence will require a larger value of frequency compared to TVES without memory at dynamic bifurcation. There could also be instances where TVES without rheology show bifurcation for some choices of parameters, but for the same choice of parameters dynamic bifurcation could be eliminated for TVES with rheology due to increase in stiffness because of higher positive $\left[K_{\sigma }\right]$ . In applications such as nonlinear dynamics of plates, beams etc. $\left[K_{\sigma }\right]$ is always positive (for simple BCs), hence we expect the dynamic bifurcation to occur at higher frequency for TVES with memory compared to TVES without memory.

G) We present a comparison of bifurcation phenomenon in TVES without memory and TVES with memory. In doing so, we assume that the structure of matrices in TVES without memory is implicitly present in case of TVES with memory with additional physics of rheology. This can be verified by discarding relaxation term in the constitutive theory. In case of TVES without memory BLM expressed in displacement followed by GM/WF, Newmark linear acceleration with Newton’s linear method reveals the elegant structure of various matrices associated with specific physics that facilitate understanding of dynamic bifurcation phenomenon. In TVES with memory, BLM cannot be expressed purely in terms of displacements, hence structure of matrices in TVES without memory cannot be explicitly obtained for TVES with memory. Nonetheless, consideration of common physics between TVES with and without memory is compelling and strong argument (in addition to obtaining same PDEs as for TVES without memory when ${\lambda }_1=0$ in the PDEs for TVE with memory) to use the $\left[H_T\right]$ ; $\left[K_T\right]$ and structure of other matrices in TVES without memory to explain the contribution of additional physics of rheology on dynamic bifurcation in case of TVES with memory.

Thus, when we compare the physics of TVES with rheology and TVES without rheology we find:

1) TVES with rheology have lower (linear and nonlinear) stiffness compared TVES without rheology due to presence of long chain molecules. As the concentration of long chain molecules increases (dense polymeric solid), progressively increasing rheology, the overall stiffness decreases. Thus, for the same external stimulus a TVES with progressively increasing rheology will yield progressively increasing deformation (displacements). Thus, in frequency response, we expect higher amplitudes for TVES with rheology compared to TVES without rheology.

2) Rheology increases the residual stress field in compression as well as tension, thus increases negative $\left[K_{\sigma }\right]$ as well as positive $\left[K_{\sigma }\right]$ , but only in nonlinear case as there is no $\left[K_{\sigma }\right]$ for linear case.

3) In TVES, the dynamic bifurcation depends upon the following:

a) We except same contributions to $\left[H_T\right]$ (though $\left[H_T\right]$ in the same form as for TVES without memory is not possible in this case) as in case of TVES without memory but the stiffness is reduced due to long chain molecules, thus effective stiffness is lower in TVES with rheology compared to TVES without memory.

b) $\left[K_{\sigma }\right]$ as in case of TVES without memory. When the stress field is compressive, presence of rheology enhances negative stress field which in turn increases negative $\left[K_{\sigma }\right]$ , thus reducing effective stiffness.

c) Translational inertial physics. Since in (a) effective stiffness is reduced and in (b) effective stiffness is also reduced, thus requires reduction in stiffness due to translational inertial physics, implying that dynamic bifurcation should happen at a lower frequency i.e., to the left of the TVES without memory.

d) When $\left[K_{\sigma }\right]$ is positive, the situation can be quite different. When compared with TVES without memory, we observe and note the following.

i) Lower stiffness (linear and nonlinear without $\left[K_{\sigma }\right]$ ) in case of TVES with memory. This reduction compared to TVES without memory depends upon rheology. Let $\left[\Delta K\right]$ be the reduction in stiffness due to rheology.

ii) Increase in tensile stress field due to rheology. Hence, increase in positive $\left[K_{\sigma }\right]$ compared to TVES without memory. Let $\left[\Delta K_{\sigma }\right]$ be the increase in stiffness due to additional stress field because of rheology.

iii) If $\left[\Delta K\right]$ is greater than $\left[\Delta K_{\sigma }\right]$ then the dynamic bifurcation will occur at a frequency lower than that for corresponding TVES without memory.

iv) If $\left[\Delta K_{\sigma }\right]$ is greater than $\left[\Delta K\right]$ , then the dynamic bifurcation will occur at a frequency greater than that for corresponding TVES without memory.

We remark that stress relaxation phenomenon is primarily responsible for change in the stress field that depends upon relaxation time. For smaller relaxation time, the stress free state is achieved quicker (in terms of clock time) and vice versa. Thus, in general, higher the relaxation time, more is the deviation between the bifurcation frequencies for TVES without memory and TVES with memory.

8. Model Problems

In references [13] [14] , we had considered phenomenological mathematical models as well as those derived using CBL of CCM. In case of TVES with memory, the deviatoric second Piola-Kirchhoff stress tensor constitutive theory cannot be substituted in the BLM as the constitutive theory is PDE in time. Nonetheless, we find that in almost all published works, this is done in some manner or other. Thus, our view is that these mathematical models have inconsistencies or serious assumptions. For this reason, in the present work we only consider mathematical models for the model problems based on CBL of CCM and consistent constitutive theories.

We consider the same model prolem of 1D axial TVE rod without memory considered in [14] but with rheology so that we can illustrate the influence of rheology on the nonlinear dynamics and dynamic bifurcation. We consider a TVE axial rod with rheology of length L, fixed at $x_1=0$ (left end) and subjected to harmonic excitation $f_0\sin \left(\omega t\right)$ at $x_1=L$ (right end), as shown in Figure 1(a). Assuming that all points in each cross section of the rod deform in the $x_1$ direction by the same amount, we can idealize the rod by a line shown in Figure 1(b). For simplicity, we consider isothermal physics. Thus, balance of linear momenta and the constitutive theory for deviatoric second Piola-Kirchhoff

Figure 1. Schematic, discretization, BCs and ICs for an axial rod (model problem II). (a) Schematic; (b) Space-time domain; (c) Discretization of ${\bar{\Omega }}_{xt}$ in space-time strip: ${\bar{\Omega }}_{xt}^T={\displaystyle \underset{i}{\cup }{\bar{\Omega }}_{xt}^{\left(i\right)}}$ ; (d) Discretization ${\left({\bar{\Omega }}_{xt}^{\left(1\right)}\right)}^T$ of the first space-time strip ${\bar{\Omega }}_{xt}^{\left(1\right)}$ : ${\left({\bar{\Omega }}_{xt}^{\left(1\right)}\right)}^T={\displaystyle \underset{e}{\cup }{\bar{\Omega }}_{xt}^e}$ ; (e) BCs and ICs.

stress are the only equations that constitute the mathematical model. The constitutive theory for the equilibrium second Piola-Kirchhoff stress is defined using equation of state [1] [2] , hence is known. In the derivation of the constitutive theories, we consider the dissipation mechanism to be dependent on Green’s strain rates up to order “n” and the rheology to be dependent upon the rates of the deviatoric second Piola-Kirchhoff stress tensor up to order m (see reference [1] [2] ). The deviatoric second Piola-Kirchhoff stress rate of order m is used a constitutive variable and ${}_d\sigma {}^{\left[i\right]}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{0,1,}\cdots \mathrm{,}m-1$ are considered as argument tensors of ${}_d\sigma {}^{\left[m\right]}$ . In the present numerical studies, we choose $m=1$ and $n=1$ for which we can obtain the following:

${\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){\sigma }_{11}^{\left[0\right]}\right)=0\left(\text{BLM}\right)$ (175)

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

${}_d\sigma {}_{11}^{\left[0\right]}+{\lambda }_1\frac{\partial \left({}_d\sigma {}_{11}^{\left[0\right]}\right)}{\partial t}=2\mu {\left({\epsilon }_{\left[0\right]}\right)}_{11}+\lambda {\left({\epsilon }_{\left[0\right]}\right)}_{11}+2{\mu }_1{\left({\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\right)}_{11}+{\lambda }_1{\left({\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\right)}_{11}$ (177)

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

${\left({\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\right)}_{11}=\frac{\partial }{\partial t}{\left({\epsilon }_{\left[0\right]}\right)}_{11}$ (179)

We note that for 1D case $2\mu +\lambda =E$ , modulus of elasticity and $\left(2{\mu }_1+{\lambda }_1\right)=\eta$ , viscosity. ${\lambda }_1$ in (177) is relaxation time. We nondimensionalize (175) - (179) with hat ( $\stackrel{^}{}$ ) on all quantities indicating they have their usual dimensions or units (neglecting equilibrium stress)

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

${}_d\hat{\sigma }{}_{11}^{\left[0\right]}+{\hat{\lambda }}_1\frac{\partial }{\partial {\hat{x}}_1}\left({}_d\hat{\sigma }{}_{11}^{\left[0\right]}\right)=\hat{E}{\left({\hat{\epsilon }}_{\left[0\right]}\right)}_{11}+\hat{\eta }{\left({\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\right)}_{11}$ (181)

${\left({\hat{\epsilon }}_{\left[0\right]}\right)}_{11}=\frac{\partial {\hat{u}}_1}{\partial {\hat{x}}_1}+\frac{1}{2}{\left(\frac{\partial {\hat{u}}_1}{\partial {\hat{x}}_1}\right)}^2$ (182)

${\left({\stackrel{\dot{}}{\epsilon }}_{\left[0\right]}\right)}_{11}=\frac{\partial }{\partial t}\left(\frac{\partial {\hat{u}}_1}{\partial {\hat{x}}_1}+\frac{1}{2}{\left(\frac{\partial {\hat{u}}_1}{\partial {\hat{x}}_1}\right)}^2\right)$ (183)

$\forall \left(x\mathrm{,}t\right)\in {\Omega }_{xt}$

$\begin{array}{l}\begin{array}{l}{\rho }_0=\frac{{\hat{\rho }}_0}{{\left({\rho }_0\right)}_{ref}};u_1=\frac{{\hat{u}}_1}{L_0};v_1=\frac{{\hat{v}}_1}{v_0};x=\frac{{\hat{x}}_1}{L_0};t_0=\frac{L_0}{v_0};\\ {}_d\sigma {}_{11}^{\left[0\right]}=\frac{{}_d\hat{\sigma }{}_{11}^{[0]}}{{\sigma }_0};E=\frac{\hat{E}}{E_0};\eta =\frac{\hat{\eta }}{{\eta }_0};{\lambda }_1=\frac{{\hat{\lambda }}_1}{t_0}=De\\ {\sigma }_0=E_0={\rho }_0{v}_0^2;v_0=\sqrt{\frac{E_0}{{\left({\rho }_0\right)}_{ref}}};F_1^b=\frac{{\hat{F}}_1^b}{F_0};F_0=\frac{L_0}{v_0^2}\end{array}\}\hfill \end{array}$ (184)

$v_0$ is reference speed of sound. ${\tau }_0\mathrm{,}{E}_0\mathrm{,}{\left({\rho }_0\right)}_{ref}$ and $t_0\mathrm{,}{\eta }_0$ are reference stress, modulus of elasticity, density, reference time and reference viscosity.

We substitute ${\left({\epsilon }_{\left[0\right]}\right)}_{11}$ and ${\left({\epsilon }_{\left[0\right]}^2\right)}_{11}$ form (182) and (183) into (181), expand the third term in (180), introduce $v_1=\frac{\partial u_1}{\partial t}$ and then nondimensionalize using (184) to obtain the following:

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

${}_d\sigma {}_{11}^{\left[0\right]}+De\frac{\partial }{\partial t}\left({}_d\sigma {}_{11}^{\left[0\right]}\right)-E\left(\frac{\partial u_1}{\partial x_1}+\frac{1}{2}{\left(\frac{\partial u_1}{\partial x_1}\right)}^2\right)-c\left(\frac{\partial v_1}{\partial x_1}+\frac{\partial u_1}{\partial x_1}\frac{\partial v_1}{\partial x_1}\right)=0$ (186)

$v_1-\frac{\partial u_1}{\partial t}=0$ (187)

$\forall \left(x\mathrm{,}t\right)\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t=\left(\mathrm{0,}L\right)\times \left(\mathrm{0,}\tau \right)$ in which $E=1$ if $E_0=\hat{E}$ and $c=\frac{\hat{\eta }{v}_0}{{\sigma }_0{L}_0}=\frac{\hat{\eta }{v}_0}{{\left({\rho }_0\right)}_{ref}{\sigma }_0{L}_0}=\eta \left(\frac{{\eta }_0}{{\left({\rho }_0\right)}_{ref}{\sigma }_0{L}_0}\right)=\frac{\eta }{Re}$ . where $Re=\frac{\left({\rho }_0\right)v_0{L}_0}{{\eta }_0}$ . c is dimensionless damping coefficient and $De=\frac{{\hat{\lambda }}_1}{t_0}$ is dimensionless relaxation time (Deborah number). We calculate solution of (185) - (187) using space-time finite element method based on space-time residual functional for a space-time strip with time marching. Figure 1(c) shows space-time domain ${\bar{\Omega }}_{xt}={\bar{\Omega }}_x\times {\bar{\Omega }}_t=\left[\mathrm{0,}L\right]\times \left[\mathrm{0,}\tau \right]$ , $\tau$ being value of time. A discretization ${\bar{\Omega }}_{xt}^T={\displaystyle \underset{i}{\cup }{\bar{\Omega }}_{xt}^{\left(i\right)}}$ of space time domain ${\bar{\Omega }}_{xt}$ into space-time strips ${\bar{\Omega }}_{xt}^{\left(i\right)}\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots$ is shown in Figure 1(c).

In the computational studies choose $L=1$ and dimensionless density ${\rho }_0=0.001$ . Thus, the dimensionless wave speed in the linear elastic medium is $v=\sqrt{\frac{E}{{\rho }_0}}=\sqrt{\frac{1}{0.001}}=31.6218$ and the natural frequencies of axial vibration (small strain, linear elasticity without dissipation) are given by ${\omega }_n=\frac{\left(2n-1\right)\pi v}{2L}=\left(2n-1\right)49.673\mathrm{<mo>\; ;\; </mo>}n=\mathrm{1,2,}\cdots$

Figure 1(d) shows a discretization ${\left({\bar{\Omega }}_{xt}^{\left(1\right)}\right)}^T={\displaystyle \underset{r}{\cup }{\bar{\Omega }}_{xt}^e}$ of the first space-time strip using p-version nine node space-time finite elements with hierarchical local approximations of higher degree and higher order global differentiability in space and time. Figure 1(e) shows BCs: $u_1=0,v_1=0$ at $x_1=0$ $\forall t\in \left[0,\Delta t\right]$ , ICs: $u_1=0,v_1=0$ at $t=0$ $\forall x_1\in \left[0,L\right]$ . BCs at $x_1=L$ consists of ${}_d\sigma {}_{11}^{\left[0\right]}={\sigma }_0\sin \left(\omega t\right)$ . Using ${}_d\sigma {}_{11}^{\left[0\right]}$ values of the degrees of the freedom are calculated for the element nodes at $x_1=L$ and $t\in \left[\mathrm{0,}\Delta t\right]$ . For the second space-time strip, these dofs are calculated for $t\in \left[\Delta t\mathrm{,2}\Delta t\right]$ and so on. This process enables application of ${}_d\sigma {}_{11}^{\left[0\right]}={\sigma }_0\sin \left(\omega t\right)$ at $x_1=L$ in time accurately. In this process, a change in $\omega$ is done at a value of time t when ${}_d\sigma {}_{11}^{\left[0\right]}=0$ , so that change in $\omega$ does not cause unnecessary discontinuity in the application of ${}_d\sigma {}_{11}^{\left[0\right]}$ at $x_1=L$ in time. Convergence studies using mesh refinement and p-level increase confirm that a 10 element discretization with p-level of 9 in space and time with global differentiability of order one ( $k=2$ , order of the approximation space) are sufficient to obtain space-time residual functional I for ${\left({\bar{\Omega }}_{xt}^{\left(1\right)}\right)}^T$ of $O\left({10}^{-8}\right)$ or lower indicating accurate evolution for ${\left({\bar{\Omega }}_{xt}^{\left(1\right)}\right)}^T$ .

Evolution is continued using time marching till steady cyclic response is reached for each frequency. We use integration time step of $\Delta t=\frac{1}{20}\left(\frac{2\pi }{w}\right)$ . Evolution is computed for 20 cycles of each frequency. This is found to be adequate enough to reach cyclic response. During the entire evolution, I values of $O\left({10}^{-8}\right)$ are achieved, ensuring accurate evolution.

As in case of model problem II in ref [14] , here also there is a negative shift (due to compression of the rod) in the peak positive and negative values observed in the cylic response, we record and report both. In this model problem $\left[K_{\sigma }\right]$ is always negative.

We present a number of numerical studies for various combinations of ${\sigma }_0$ , damping coefficient c and Deborah number De. The objective of these model problem studies is to demonstrate the influence of rheology on the nonlinear dynamic response and dynamic bifurcation. We intentionally choose same combination of ${\sigma }_0$ and c used for nonlinear dynamic response of TVES without memory in reference [14] , but incorporate rheology through De. In the numerical studies we choose two values of De to illustrate the influence of increasing rheology on nonlinear dynamic response and dynamic bifurcation.

8.1. Linear TVES with Memory

In this study, we consider small deformation, small strain physics with linear damping and memory. We choose ${\sigma }_0=0.025$ , $c=0.003$ and $De=0.0005$ . We observe no path dependency in this case i.e., $L\to R$ and $L\leftarrow R$ frequency responses are identical as shown in Figure 2.

8.2. Linear TVES with and without Rheology

In this study we consider small deformation, small strain physics with linear damping. We choose ${\sigma }_0=0.03$ , $c=0.003$ , same as used in ref. [14] for model problem II, Section 8.3.1. case (2), figure 8(b). When $De=0$ , the frequency response shown in Figure 3 is same as in figure 8(b) of ref [14] . We observe no path dependence i.e., $L\to R$ and $L\leftarrow R$ frequency responses are identical. We also perform numerical studies for ${\sigma }_0=0.03$ , $c=0.0035$ and $De=0.0005$ . Figure 3 also shows frequency response for $De=0.0005$ for $L\to R$ and $L\leftarrow R$ . In both studies (for $De=0$ , $De=0.0005$ ) we have kept some values of ${\sigma }_0$ and damping coefficient c (0.003) so that we can demonstrate the influence of rheology. Presence of long chain molecules decreases the stiffness of the solid matter resulting in high displacement values. This is quite visible in Figure 3. Peak amplitude of 0.12 for $De=0.0$ increases to 0.15 at $De=0.0005$ . We remark that presence of rheology also increases damping coefficient c, which increases resistance to motion, hence causes reduction in amplitude of motion. Thus, with actual damping of the polymeric solid, the peak values would be lower than 0.15 but we expect it to be higher than 0.12 (for $De=0.0$ ).

Figure 2. TVES with memory: Linear frequency $\omega$ vs amplitude u response for $c=0.003$ .

Figure 3. TVES without memory and with memory: Linear frequency $\omega$ vs amplitude u response for $c=0.003$ , $De=0.0,0.0005$ .

8.3. Fixed Values of σ₀, De with Decreasing c

In this study, we consider fixed values of ${\sigma }_0=0.030$ and two values of $c=0.003,0.002$ . We determine frequency response for $De=0.0$ and 0.0005 for $L\to R$ and $L\leftarrow R$ paths of excitation frequency for both values of c.

First, we consider ${\sigma }_0=0.030,c=0.003$ and calculate frequency response for $De=0.0$ and 0.0005. Figure 4(a) shows frequency response for $c=0.003$ and $De=0.0$ i.e., TVES without memory (Same as figure 9(c) in ref [14] ). We observe virtually no dynamic bifurcation and very little path dependencies. Figure 4(b) shows frequency response for the same values of ${\sigma }_0$ and c but for $De=0.0005$ . We observe distinct dynamic bifurcation in the neighborhood of $\omega =46$ Hz. Peak amplitude for both positive and negative peaks as higher in Figure 4(b) compared to Figure 4(a). This of course is due to lower stiffness of TVES with memory compared to TVES without memory and increased influence of negative $\left[K_{\sigma }\right]$ which results in presence of dynamic bifurcation, higher peak amplitudes for both negative and positive peaks. Wee also observe mild path dependency in the response.

Next, we consider ${\sigma }_0=0.030$ (same as before) but $c=0.002$ (reduced damping) and determine frequency response for $De=0.0005$ . Figure 4(c) shows frequency response for $De=0.0$ (same as figure 9(d) of ref [14] ). We observe clear dynamic bifurcation and some path dependency. Figure 4(d) shows frequency response for $De=0.0005$ . We see distinct bifurcation at $\omega =45$ Hz with both positive and negative peaks higher than in Figure 4(c) for the same reason as discussed for $c=0.003$ .

From this study, we clearly note that when $\left[K_{\sigma }\right]$ is compressive (as in this case), presence of rheology (hence lower stiffness) for some ${\sigma }_0$ and $c_0$ enhances dynamic bifurcation i.e., 1) more distinct in appearance; 2) results in higher peak amplitudes compare to the case without rheology due to reduced stiffness in TVES with memory and higher negative $\left[K_{\sigma }\right]$ .

8.4. Fixed Values of σ₀, c with Increasing De

In this study, we consider fixed values of ${\sigma }_0=0.025,c=0.003$ for $De=0.0,0.0005,0.001$ . First, we consider ${\sigma }_0=0.025,c=0.003$ and $De=0.0,0.0005$ . Figure 5(a) shows frequency response for $De=0.0$ . We do not observe dynamic bifurcation as well as significant path dependencies. Figure 5(b) shows frequency response for $De=0.0005$ . We observe mild dynamic bifurcation at $\omega =46.5$ Hz and very little path dependency. Both positive and negative peaks in Figure 5(b) for $De=0.0005$ are higher than those in Figure 5(a) for $De=0.0$ . Figure 5(c) shows dynamic response for $De=0.001$ . We observe distinct dynamic bifurcation at $\omega =46$ Hz. There is very little path dependency. Both positive and negative peaks are higher than for $De=0.0$ and $De=0.0005$ . Higher peaks are due to reduced stiffness for $De=0.001$ (compared to $De=0.0$ ) and enhance negative $\left[K_{\sigma }\right]$ .

(a)(b)(c)(d)

Figure 4. (a) TVES without memory: Frequency $\omega$ vs amplitude u response for $c=0.003$ , $De=0.0$ ; (b) TVES with memory: Frequency $\omega$ vs amplitude u response for $c=0.003$ ; (c) TVES without memory: Frequency $\omega$ vs amplitude u response for $c=0.002$ , $De=0.0$ ; (d) TVES with memory: Frequency $\omega$ vs amplitude u response for $c=0.003$ , $De=0.0005$ .

8.5. Fixed Values of σ₀, c with Increasing De

This study is similar to Section 8.3 except that in this case ${\sigma }_0=0.03$ , higher than 0.025 used in Section 8.3. We use $c=0.003$ and $De=0.0,0.0005,0.001$ . Figure 6(a) shows frequency response for $De=0.0$ . We observe mild bifurcation, but insignificant path dependency. Both positive and negative peaks are higher for $De\ne 0.0$ compared to $De=0.0$ as shown in Figure 6(b) and Figure 6(c). Furthermore, the positive and negative peaks for $De=0.001$ are even higher than those for $De=0.0005$ . This is due to the same reason stated in earlier two studies.

1) From the studies presented here for the linear dynamic with linear damping, we confirm that peak value of the amplitude in the frequency response is higher for $De\ne 0$ compared to $De=0$ . That is for a TVES without rheology the presence of long chain molecules result in reduced stiffness, hence increased peak amplitudes. We keep in mind that when $De\ne 0$ , the compressive stress field is enhanced due to rheology, but it cannot cause reduction in stiffness due to absence of $\left[K_{\sigma }\right]$ in linear dynamics. In this study, same damping coefficient c is used for $De=0.0$ and $De\ne 0$ . In actual case c for $De\ne 0$ may be slightly higher (or higher) than for $De=0$ . This will result in slight reduction in the peak amplitude for $De\ne 0$ .

2) In nonlinear dynamic studies, the presence of rheology for some ${\sigma }_0$ and c, enhances compressive stress field due to rheology, hence results in increased negative $\left[K_{\sigma }\right]$ causing reduction in stiffness (over and beyond reduced stiffness

(a)(b)(c)

Figure 5. (a) TVES without memory: Frequency $\omega$ vs amplitude u response for $De=0.0$ ; (b) TVES with memory: Frequency $\omega$ vs amplitude u response for $De=0.0005$ ; (c) TVES with memory: Frequency $\omega$ vs amplitude u response for $De=0.001$ .

(a)(b)(c)

Figure 6. (a) TVES without memory: Frequency $\omega$ vs amplitude u response for $c=0.003$ , $De=0.0$ ; (b) TVES with memory: Frequency $\omega$ vs amplitude u response for $c=0.003$ , $De=0.0005$ ; (c) TVES with memory: Frequency $\omega$ vs amplitude u response for $c=0.003$ , $De=0.001$ .

of the TVES with rheology due to long chain molecules) which in turn: a) results in increased positive and negative peak values. b) If dynamic bifurcation is already present for $De=0$ , then progressively increasing De, increases likelihood of dynamic bifurcation at a lower frequency than for $De=0.0$ progressively enhance dynamic bifurcation. c) If dynamic bifurcation is not present for $De=0$ , then progressively increasing De enhances the likelihood of the presence of dynamic bifurcation.

3) In all the studies presented here we have kept damping coefficient “c” same for TVES without and with memory, knowing fully well that in the presence of rheology the higher polymer viscosity ( ${\eta }_p$ ) (compared to the solvent viscosity ( ${\eta }_s$ )) will undoubtedly result in a polymeric solid with viscosity $\eta$ greater than ${\eta }_s$ , hence damping coefficient higher than that of the solvent ( $De=0$ ). This is intentional so that we could clearly demonstrate the influence of rheology without bringing in another aspect of physics due to change in damping coefficient c.

4) In reality for TVES without memory c is much lower than those of TVES with memory. Thus, if we use actual damping coefficient, a higher value of ${\sigma }_0$ may be required to reach unstable region of dynamic bifurcation, but we see no issues here. This has been confirmed in our studies.

5) Unfortunately, for these model problems, studies with positive $\left[K_{\sigma }\right]$ are not possible to conduct. These will be presented in a follow up paper on nonlinear dynamics of plates and shells.

9. Summary and Conclusions

In the following, we summarize the work presented in the paper and draw some conclusions.

1) The CBL of CCM derived using contravariant second Piola-Kirchhoff stress tensor and the Green’s strain tensor augmented with consistent ordered rate constitutive theories of orders (m, n) for deviatoric contravariant second Piola-Kirchhoff stress tensor based on Green’s strain rates up to orders n and the deviatoric second Piola-Kirchhoff stress tensor rates of up to order m completely define the finite deformation, finite strain nonlinear dynamics of TVES with rheology.

2) As in case of TVES without memory [14] , here also the factors influencing nonlinear dynamics and dynamic bifurcation are difficult to establish using the mathematical model described in (1).

3) In case of TVES with memory, the constitutive theory for deviatoric second Piola-Kirchhoff stress tensor is a differential equation in deviatoric second Piola-Kirchhoff stress tensor and time. Thus, in this case the deviatoric second Piola-Kirchhoff stress tensor cannot be substituted in BLM. The consequence of this is that we cannot obtain BLM purely in terms of spatial and time derivatives of the displacement. This eliminates the possibility of using GM/WF for BLM to obtain nonlinear ODEs in time containing mass, damping and stiffness matrices and dofs related to displacements, velocities and accelerations: ( $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ ). That is, the following system of ODEs (valid for TVES without rheology) are not possible for TVES with memory (shown for $n=1$ and $m=1$ ):

$\left[M\right]\left\{\ddot{\delta }\right\}+\left[C\right]\left\{\stackrel{\dot{}}{\delta }\right\}+\left[K\right]\left\{\delta \right\}=\left\{F\right\}+\left\{P\right\}$ (188)

4) Based on item (3), it is straight forward to conclude that all mathematical models in published works (phenomenological or otherwise) based on similar nonlinear ODEs as in item (3) for TVES with memory are in violation of CBL of CCM and the constitutive theories, hence cannot possibly describe the nonlinear dynamics of TVES with rheology correctly.

5) Item (3) clearly suggests that the approach used for identifying factors influencing nonlinear dynamics and dynamic bifurcation for TVES without memory (ref [14] ) based on (188) and the time integration scheme cannot be used for TVES with memory as in this case (188) cannot be derived. Thus, we must consider an alternative approach to determine the various factors influencing nonlinear dynamics and dynamic bifurcation in TVES with rheology.

6) We note only the presence of long chain molecules in TVES with rheology distinguishes them from the TVES without rheology. This suggests that perhaps all of the physics of TVES without rheology is also present in TVES with rheology. The presence of long chain molecules brings in new physics of rheology that is not present in TVES without memory and may result in some modifications regarding the intensities of the various physics present in TVES without memory. The validity of the above suggestion becomes obvious when we set relaxation time ${\lambda }_1=0$ in the mathematical model to recover the precise mathematical model for TVES without memory for which (188) and the factors influencing nonlinear dynamics and dynamic bifurcation presented in reference [14] remain valid. Thus, we can conclude that even though (188) is not possible in TVES with memory, but based on the basis of the arguments given above the factors influencing nonlinear dynamics and dynamic bifurcation in TVES without memory remain valid for TVES with memory. Additionally, we need to consider the influence of rheology which is not present in TVES without memory.

7) Presence of long chain molecules reduces stiffness, but increases resistance due to dissipation (higher damping coefficient). The overall stiffness of the polymeric solid is lower than its counterpart without rheology. The frequency response presented for linear TVES with and without memory illustrates this quite well when the damping coefficient is the same for TVES with and without memory. Increased damping coefficient for TVES with memory will result in some reduction in the amplitude but the response will still have amplitudes higher than those for TVES without memory. This is the first significant difference between TVES with memory compared to TVES without memory.

8) The second major difference between TVES with memory and those without memory is the presence of rheology or relaxation phenomenon in TVES with memory. In an externally applied moving disturbance, the stress behind the disturbance does not attain zero value immediately after the movement of the disturbance. The relaxation modulus which is a function of relaxation time controls the physical time for the stress behind the disturbance to reach zero stress state. The consequence of this is that both tensile and compressive stress fields are enhanced in the case of TVES with rheology compared to TVES without rheology due to the presence of additional stresses because of rheology. The precise increase in the stress field depends upon the relaxation time. In any case, $\left[K_{\sigma }\right]$ , both positive and negative are enhanced in the case of TVES with rheology compared to TVES without memory (of course for the same material coefficients).

9) The dynamic bifurcation in TVES with memory depends upon the same physics and factors as in the case of TVES without memory [14] but additionally we must take into account the reduced stiffness due to long chain molecules and enhanced $\left[K_{\sigma }\right]$ due to increased stress field. We consider the following possibilities.

a) When $\left[K_{\sigma }\right]$ is negative, but enhanced compared to $\left[K_{\sigma }\right]$ for TVES without memory, then reduced stiffness due to long chain molcules $\left[\Delta K\right]$ and increased negative $\left[K_{\sigma }\right]$ ( $\left[\Delta K_{\sigma }\right]$ ) will result in total stiffness that is lower than that for TVES without memory, hence reduced influence of the translational inertial physics is needed for dynamic bifurcation, implying that dynamic bifurcation will occur at a lower frequency compared to TVES without memory i.e., to the left of the dynamic bifurcation for corresponding TVES without rheology.

b) When $\left[K_{\sigma }\right]$ is positive, but enhanced compared to $\left[K_{\sigma }\right]$ for TVES without memory, then the reduced stiffness due to rheology ( $\left[\Delta K\right]$ ) and increased stiffness due to stress field $\left[\Delta K_{\sigma }\right]$ control bifurcation frequency. If $\left[\Delta K\right]$ is lower than $\left[\Delta K_{\sigma }\right]$ , then the dynamic bifurcation will occur at a higher frequency compared to corresponding TVES without memory and vice versa.

10) Model problem studies presented in the paper demonstrate all aspects discussed above here except positive $\left[K_{\sigma }\right]$ , as this is not possible for the model problem considered in the present study. Model problems clearly illustrate the influence of rheology on the dynamic bifurcation when compared with TVES without memory.

11) The work presented in this paper establishes the factors influencing dynamic bifurcation in TVES with memory compared to TVES without memory and illustrates their validity in the model problem studies.

12) The space-time coupled finite element method based on residual functional for a space-time strip or a space-time slab with time marching used in the present work is the best methodology for obtaining solutions of the nonlinear PDEs in space and time. In this methodology, space-time integral forms are space-time variationally consistent, hence ensuring unconditional stability of computations during the entire evolution. In this approach, use of minimally conforming approximation spaces in space and time and with appropriate choices of h and p, the space time residual functional can be ensured to be of the order of $O\left({10}^{-8}\right)$ or lower, ensuring that the PDEs in the mathematical model are satisfied accurately at each point in the space-time domain. This ensures that the calculated solutions are almost time accurate.

13) In the space-time decoupled methods of approximation, the approximations due to space-time decoupling, accuracy and stability of the time integration method and the inability of these techniques to quantitatively judge how well the PDEs are satisfied are of concern. This approach is undoubtedly inferior to the space-time coupled methods of approximation.

14) As in the case of TVES without memory, here also we expect path dependency due to nonlinear PDEs with irreversibility. The severity of path dependency depends upon the severity of entropy production.

This paper presents thermodynamically and mathematically consistent mathematical model based on CBL of CCM for nonlinear dynamics and dynamic bifurcation physics and methods of obtaining their solutions. Space-time coupled finite element method can be almost time accurate and has no issues of stability. Thus, reported solutions are true solutions of the PDEs in the mathematical model. As we mentioned in our previous paper on TVES without memory, whether the solutions reported here match experimental data is another issue. If they do, we have a computational tool that eliminates experiments. If they do not, then obviously either the physics considered in the two differ or there are other sources of errors that are contaminating the results. We remark that in the solutions presented in this paper, the only source of confusion can be a lack of consideration of some physics that is present in the experiment. In view of the fact that the material presented in this paper is strictly based on CBL of CCM and the solutions of the IVPs are obtained using space-time coupled finite element method that is unconditionally stable and provides control over the accuracy of the solution for nonlinear dynamics and dynamic bifurcation, phenomenological mathematical models and questionable methods of obtaining their solutions can be completely avoided.

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.

Appendix

Nomenclature

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

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

${\rho }_0$ : reference density

$\rho$ : density in Lagrangian description

$\bar{\rho }$ : density in Eulerian description

$\eta$ : specific entropy in Lagrangian description

$\bar{\eta }$ : specific entropy in Eulerian description

$e$ : specific internal energy in Lagrangian description

$\bar{e}$ : specific internal energy in Eulerian description

$J$ : deformation gradient tensor in Lagrangian description

${}_{s}J$ : symmetric part of deformation gradient tensor in Lagrangian description

${}_{a}J$ : skew-symmetric part of deformation gradient tensor in Lagrangian description

${}^{d}J$ : displacement gradient tensor in Lagrangian description

${}_a{}^{{}_i\Theta }J$ : skew-symmetric part of internal rotation gradient tensor in Lagrangian description

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

$\bar{q}$ , ${\bar{q}}_i$ , $\left\{\bar{q}\right\}$ : heat vector in Eulerian description

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

$\bar{v}$ , ${\bar{v}}_i$ , $\left\{\bar{v}\right\}$ : velocities in Eulerian description

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

$\bar{u}$ , ${\bar{u}}_i$ , $\left\{\bar{u}\right\}$ : displacements in Eulerian description

$P$ : average stress in Lagrangian description

$\bar{P}$ : average stress in Eulerian description

$\sigma$ , ${\sigma }_{ij}$ , $\left[\sigma \right]$ : Cauchy stress tensor in Lagrangian description

${\sigma }^{\left(0\right)}$ , ${\bar{\sigma }}_{ij}$ , $\left[\bar{\sigma }\right]$ : Cauchy stress tensor in Eulerian description

${}_s\sigma$ : symmetric part of Cauchy stress tensor tensor

${}_d\left({}_s\sigma \right)$ : deviatoric part of the symmetric Cauchy stress tensor tensor

${}_e\left({}_s\sigma \right)$ : equilibrium part of the symmetric Cauchy stress tensor tensor

$\theta$ : temperature in Lagrangian description

$\bar{\theta }$ : temperature in Eulerian description

$k$ : thermal conductivity in Lagrangian

$p$ : thermodynamic or Mechanical Pressure in Lagrangian description

$\bar{p}$ : thermodynamic or Mechanical Pressure in Eulerian description

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

$\bar{g}$ , ${\bar{g}}_i$ , $\left\{\bar{g}\right\}$ : temperature gradient tensor in Eulerian description

$\bar{L}$ : velocity gradient tensor in Eulerian description

$\bar{D}$ : symmetric part of velocity gradient tensor in Eulerian description

TES: thermoelastic solid

TVES: thermoviscoelastic solid

Conflicts of Interest

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

References

[1]	Surana, K.S. (2015) Advanced Mechanics of Continua. CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/b17959
[2]	Surana, K.S. (2022) Classical Continuum Mechanics. 2nd Edition, CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/9781003105336
[3]	Surana, K. and Alverio, E. (2020) Consistency and Validity of the Mathematical Models and the Solution Methods for BVPS and IVPS Based on Energy Methods and Principle of Virtual Work for Homogeneous Isotropic and Non-Homogeneous Non-Isotropic Solid Continua. Applied Mathematics, 11, 546-578. https://doi.org/10.4236/am.2020.117039
[4]	Surana, K.S. and Reddy, J.N. (2017) The Finite Element Method for Initial Value Problems. CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/b22512
[5]	Surana, K.S. and Reddy, J.N. (2017) The Finite Element Method for Boundary Value Problems: Mathematics and Computations. CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/b22512
[6]	Byron Bird, R., Armstrong, R.C. and Hassager, O. (1987) Dynamics of Polymeric Liquids: Fluid Mechanics. 2nd Edition, Wiley, New York.
[7]	Surana, K.S., Nunez, D., Reddy, J.N. and Romkes, A. (2014) Rate Constitutive Theories for Ordered Thermoviscoelastic Fluids: Polymers. Continuum Mechanics and Thermodynamics, 26, 143-181. https://doi.org/10.1007/s00161-013-0295-8
[8]	Surana, K.S., Reddy, J.N. and Nunez, D. (2015) Ordered Rate Constitutive Theories for Thermoviscoelastic Solids without Memory in Lagrangian Description Using Gibbs Potential. Continuum Mechanics and Thermodynamics, 27, 409-431. https://doi.org/10.1007/s00161-014-0366-5
[9]	Amabili, M. (2019) Derivation of Nonlinear Damping from Viscoelasticity in Case of Nonlinear Vibrations. Nonlinear Dynamics, 97, 1785-1797. https://doi.org/10.1007/s11071-018-4312-0
[10]	Amabili, M. (2018) Nonlinear Damping in Nonlinear Vibrations of Rectangular Plates: Derivation from Viscoelasticity and Experimental Validation. Journal of the Mechanics and Physics of Solids, 118, 275-292. https://doi.org/10.1016/j.jmps.2018.06.004
[11]	Khaniki, H.B., Ghayesh, M.H., Chin, R., et al. (2022) A Review on the Nonlinear Dynamics of Hyperelastic Structures. Nonlinear Dynamics, 110, 963-994. https://doi.org/10.1007/s11071-022-07700-3
[12]	Balasubramanian, P., Ferrari, G., Amabili, M. and Del Prado, Z.G.J.N. (2017) Experimental and Theoretical Study on Large Amplitude Vibrations of Clamped Rubber Plates. International Journal of Non-Linear Mechanics, 94, 36-45. https://doi.org/10.1016/j.ijnonlinmec.2016.12.006
[13]	Surana, K.S. and Mathi, S.S.C. (2020) A Thermodynamically Consistent Non-Linear Mathematical Model for Thermoviscoelastic Plates/Shells with Finite Deformation and Finite Strain Based on Classical Continuum Mechanics. International Journal of Non-Linear Mechanics, 126, Article ID: 103565. https://doi.org/10.1016/j.ijnonlinmec.2020.103565
[14]	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
[15]	Fung, Y.C. (1965) Foundations of Solid Mechanics. Prentice-Hall, Englewood Cliffs.
[16]	Christensen, R.M. (1982) Theory of Viscoelasticity. 2nd Edition, Dover, Mineola. https://doi.org/10.1016/B978-0-12-174252-2.50012-0
[17]	Lakes, R.S. (2009) Viscoelastic Materials. Cambridge University Press, New York. https://doi.org/10.1017/CBO9780511626722
[18]	Surana, K.S. and Mathi, S.S.C. (2020) A Thermodynamically Consistent Formulation for Dynamic Response of Thermoviscoelastic Plate/Shell Based on Classical Continuum Mechanics (CCM). International Journal of Structural Stability and Dynamics, 20, Article ID: 2050043. https://doi.org/10.1142/S0219455420430129
[19]	Surana, K.S., Mysore, D. and Carranza, C.H. (2020) A Thermodynamically Consistent Formulation for Bending of Thermoviscoelastic Beams for Small Deformation, Small Strain Based on Classical Continuum Mechanics. Mechanics of Advanced Materials and Structures, 27, 1120-1140. https://doi.org/10.1080/15376494.2020.1725987
[20]	Surana, K.S. and Abboud, E. (2023) Deviatoric Stress Waves Due to Rheology in Incompressible Thermoviscoelastic Solid Medium with Small Strain, Small Deformation Physics. Meccanica, 58, 755-779. https://doi.org/10.1007/s11012-023-01646-5
[21]	Surana, K. and Kitchen, M. (2022) Stress Waves in Polymeric Fluids. American Journal of Computational Mathematics, 12, 87-118. https://doi.org/10.4236/ajcm.2022.121007
[22]	Zheng, Q.S. (1993) On the Representations for Isotropic Vector-Valued, Symmetric Tensor-Valued and Skew-Symmetric Tensor-Valued Functions. International Journal of Engineering Science, 31, 1013-1024. https://doi.org/10.1016/0020-7225(93)90109-8
[23]	Wang, C.C. (1969) On Representations for Isotropic Functions: Part I. Isotropic Functions of Symmetric Tensors and Vectors. Archive for Rational Mechanics and Analysis, 33, 249-267. https://doi.org/10.1007/BF00281278
[24]	Wang, C.C. (1969) On Representations for Isotropic Functions: Part II. Isotropic Functions of Skew-Symmetric Tensors, Symmetric Tensors, and Vectors. Archive for Rational Mechanics and Analysis, 33, 268-287. https://doi.org/10.1007/BF00281279
[25]	Wang, C.C. (1970) A New Representation Theorem for Isotropic Functions, Part I and Part II. Archive for Rational Mechanics and Analysis, 36, 166-223.
[26]	Wang, C.C. (1971) Corrigendum to “Representations for Isotropic Functions”. Archive for Rational Mechanics and Analysis, 43, 392-395. https://doi.org/10.1007/BF00252004
[27]	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
[28]	Spencer, A.J.M. and Rivlin, R.S. (1959) 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
[29]	Spencer, A.J.M. (1971) Theory of Invariants. In: Eringen, A.C., Ed., Mathematics, Academic Press, New York, 239-353. https://doi.org/10.1016/B978-0-12-240801-4.50008-X