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

Karan S. Surana, Sri Sai Charan Mathi

Department of Mechanical Engineering, Univeristy of Kansas, Lawrence, KS, USA.

DOI: 10.4236/am.2023.1412047                

Abstract

This paper presents the mathematical model consisting of conservation and balance laws (CBL) of classical continuum mechanics (CCM) and the constitutive theories derived using entropy inequality and representation theorem for thermoviscoelastic solids (TVES) matter without memory. The CBL and the constitutive theories take into account finite deformation and finite strain deformation physics. This mathematical model is thermodynamically and mathematically consistent and is ideally suited to study nonlinear dynamics of TVES and dynamic bifurcation and is used in the work presented in this paper. The finite element formulations are constructed for obtaining the solution of the initial value problems (IVPs) described by the mathematical models. Both space-time coupled as well as space-time decoupled finite element methods are considered for obtaining solutions of the IVPs. Space-time coupled finite element formulations based on space-time residual functional (STRF) that yield space-time variationally consistent space-time integral forms are considered. This approach ensures unconditional stability of the computations during the entire evolution. In the space-time decoupled finite element method based on Galerkin method with weak form for spatial discretization, the solutions of nonlinear ODEs in time resulting from the decoupling of space and time are obtained using Newmark linear acceleration method. Newton’s linear method is used to obtain converged solution for the nonlinear system of algebraic equations at each time step in the Newmark method. The different aspects of the deformation physics leading to the factors that influence nonlinear dynamic response and dynamic bifurcation are established using the proposed mathematical model, the solution method and their validity is demonstrated through model problem studies presented in this paper. Energy methods and superposition techniques in any form including those used in obtaining solutions are neither advocated nor used in the present work as these are not supported by calculus of variations and mathematical classification of differential operators appearing in nonlinear dynamics. The primary focus of the paper is to address various aspects of the deformation physics in nonlinear dynamics and their influence on dynamic bifurcation phenomenon using mathematical models strictly based on CBL of CCM using reliable unconditionally stable space-time coupled solution methods, which ensure solution accuracy or errors in the calculated solution are always identified. Many model problem studies are presented to further substantiate the concepts presented and discussed in the paper. Investigations presented in this paper are also compared with published works when appropriate.

Keywords

Thermodynamic Consistency, Dynamic Bifurcation, Static Bifurcation, Nonlinear Formulation, Finite Strain, Finite Deformation, Thermoviscoelastic, Classical Continuum Mechanics, Conservation and Balance Laws, Nonlinear Damping

1. Literature Review, Motivation and Scope of Work

In this section, we present some basic definitions that are used throughout the paper. This is followed by literature review, motivation and scope of work.

1.1. Bifurcation Phenomenon

In the static (BVPs) and dynamic (IVPs) finite deformation, finite strain studies in structural and solid mechanics, we utilize mathematical models (based on CCM) in which contravariant deviatoric second Piola-Kirchhoff stress tensor and covariant Green strain tensor are work conjugate pair. These mathematical models are a system of nonlinear PDEs in spatial coordinates for BVPs, and in spatial coordinates and time for IVPs. The solutions of these PDEs may exhibit path dependency (if irreversibility is present in the physics) and may also exhibit sudden change in the amplitude of motion or deformation in the loading path. This sudden change in amplitude of motion is called “bifurcation”. Bifurcation can happen during loading in case of BVPs in which case we refer to it as “static bifurcation”. If the bifurcation occurs in the IVPs, then we refer to it as “dynamic bifurcation”. In this paper, we study nonlinear dynamic response and dynamic bifurcation phenomenon in TVE solids without memory with finite deformation and finite strain.

1.2. Static Bifurcation

In case of BVPs, we are seeking stationary state of the corresponding IVPs i.e., the solutions of corresponding IVP when it ceases to change as time elapses. Thus, time and time dependent physics is absent in the mathematical description of BVPs. If one constructs finite element formulation using this mathematical model based on Galerkin method with weak form (GM/WF), then we realize that the resulting finite element formulation only contains stiffness (matrices), load vector and vector of secondary variables. This suggests that in case of BVPs, if static bifurcation phenomenon (also called snap-through [1] [2] [3] [4] [5] ) exists (not all BVPs have this phenomenon) it must be purely due to stiffness. In other words, in this case the total stiffness must consist of competing aspects in it that create instability for certain combination of loading and parameters, so a sudden change in deformation state must occur to restore stable configuration, which is called static bifurcation. Works published in references [1] [2] [3] [4] [5] are examples of “static bifurcation” (also called snap through phenomenon) purely related to stiffness. Presence of finite deformation, finite strain, hence nonlinearity in the mathematical model is essential for the existence of static bifurcation.

1.3. Dynamic Bifurcation

In case of IVPs, the PDEs in the mathematical model consist of dependent variables, space coordinates and time. That is, IVPs contain time dependent physics that evolves as time elapses. The study of the time dependent solution of the mathematical model with finite strain, finite deformation with consistent nonlinear dissipation physics with or without memory (or rhegology) in structural and solid mechanics is of course nonlinear dynamics. Solutions of such IVPs naturally may also exhibit bifurcation phenomenon. We refer to this as dynamic bifurcation. In case of IVPs, compared to BVPs, there are other aspects of the physics (due to dynamics) involved that either promote and allow dynamic bifurcation to exist or discourage and inhibit existence of dynamic bifurcation. One of the main focuses of this work is determination of the factors influencing dynamic bifurcation.

Authors in references [6] presented complete details of the mathematical model in ${ℝ}^3$ for TVES without memory for finite deformation and finite strain physics. This model is derived using CBL of CCM and constitutive theories are based on entropy inequality and representation theorem [7] - [15] . The resulting mathematical model is a system of nonlinear PDEs in space and time in which dependent variables naturally exhibit simultaneous dependence on both space and time. The finite element formulations for this model were also presented in ref [6] . In contrast, the mathematical models in the published works used to study nonlinear dynamics are either primarily phenomenological or if based on CBL of CCM, contain many phenomenological modifications or adjustments that are generally not supported by CCM. We make following observations:

(1) Energy functional or principal of virtual work are employed almost exclusively as a starting point in the derivation of the mathematical descriptions of the deformation physics and for obtaining the solutions of the IVPs using methods such as finite element method.

(2) The Euler’s equation (differential form of the mathematical model) is almost never derived from the integral forms. Instead, energy methods or the principle of virtual work are used directly to obtain solutions (based on extremum of the energy functional) either using finite element method or expressing solutions in terms of superposition of periodic functions. In these approaches we never know the differential form of the actual mathematical model of the physics under consideration.

(3) Constitutive theories are never derived using entropy inequality and representation theorem in the energy method and in the principle of virtual work. This is because these methods have no concept of entropy inequality, hence have no basis or mechanism of deriving constitutive theories.

(4) Phenomenological constitutive theories are often constructed using 1D spring and 1D dashpot in suitable series-parallel arrangements. These constitutive theories are purely in time, hence lumped in space. These constitutive theories have the following drawbacks:

(a) Lack of dependence on space precludes their use for continuous matter in which both space and time are intrinsically coupled during the evolution.

(b) Phenomenological mathematical models and the constitutive theories in time based on 1D springs and dashpots cannot be extended to the deformation of continuous matter.

(c) Since there are no laws or principles that support phenomenological constitutive theories, their extensions to include more comprehensive physics can only be done phenomenologically.

(5) Published work do not adhere to valid and mathematically sound specific concepts of space-time coupled or space-time decoupled approaches for obtaining solutions of the PDEs in the mathematical models.

Keeping these points in mind, we discuss some published works that are relevant in the context of this paper. In applied mathematics, empirical nonlinear ODEs in time are used to study their nonlinear dynamic response [16] . 1D spring, dashpot and mass based phenomenological mathematical models have been used frequently in engineering to study nonlinear dynamics [17] . In reference [18] , authors use a two-mass oscillator system with 1D springs and dashpots involving two degrees of freedom to study bifurcation and backbone curves. A nonlinear second order ODE $\ddot{x}+x+x{\stackrel{\dot{}}{x}}^2+f\left(x\right)=0$ with $f\left(x\right)=a_3{x}^3+a_5{x}^3+a_7{x}^7+\cdots$ is used to investigate backbone curve in ref. [19] . Physical significance of this mathematical model in the context of nonlinear structural dynamics is neither discussed nor explained in this paper. In reference [20] , a mathematical model consisting of ODEs in time with constant coefficient mass, stiffness and damping matrices is used for identification of nonlinear stiffness in backbone curves. This mathematical model is derived based on total potential energy of the system (Appendix B), and the equations are transformed using modal basis of the linear problem. All nonlinearities are lumped into the force vector. Reference [21] uses the same mathematical model as reference [20] to investigate the tracking of backbone curves. Reference [22] also uses same system of ODEs in time as in reference [20] [21] , with nonlinearities lumped into a vector for identification of backbone curves using control-based continuation. Reference [23] presents nonlinear vibrations of water filled circular cylindrical shells. The mathematical models are based on energy functional, and fluid-solid interactions intrinsic to the problem are only accounted for by using kinetic energy considerations related to the fluid. In reference [24] a phenomenological model based on 1D nonlinear springs, dashpot and mass (single nonlinear ODE in time) is used as a basis to study large amplitude vibrations of plates and shells. In reference [25] , 1D nonlinear spring, dashpot and mass model (a single nonlinear ODE) is used to study large amplitude vibrations. In reference [25] a system of second order ODEs in time (not derived using CBL of CCM or NCCM) is used to study backbone curves in nonlinear mechanical systems. System of ODEs in time similar to references [20] [21] [22] with nonlinearities lumped into a vector with periodic function expansion of the solution is considered in reference [26] for investigating nonlinear dynamic response. Other works based on similar mathematical models and techniques described above can be also be found in references [27] [28] [29] [30] .

In summary, the mathematical models used to investigate large amplitude nonlinear dynamics in the presence of damping fall into two categories. The first category comprises phenomenological models based on 1D nonlinear springs, dampers and masses. The second category consists of linear system of second order ODEs with constant coefficient mass, damping and stiffness matrices in which all nonlinear effects (which are not explained or defined in many works) are accounted for in the non homogeneous part i.e., force vector. To the best of our knowledge, there are no published studies on large amplitude nonlinear dynamics in the presence of damping that utilizes a mathematical model based on CBL of CCM with consistent constitutive theories derived from entropy inequality and representation theorem without any phenomenological or any other modifications.

In reference [6] authors presented a complete mathematical model for TVE solid undergoing finite deformation, finite strain using contravariant second Piola Kirchhoff stress tensor, Green’s strain and rates of Green’s strain tensor upto order n. The constitutive theory for deviatoric second Piola-Kirchhoff stress tensor was derived using entropy inequality and representation theorem [7] - [15] by considering Green’s strain tensor, rates of Green’s strain tensor and temperature as its argument tensors. The resulting mathematical model consists of a system of nonlinear PDEs in both space and time. By using a space-time decoupled finite element formulation and Galerkin method with weak form (GM/WF) in space, the authors showed that the nonlinear PDEs in the mathematical model can be converted to a system of second order nonlinear ODEs in time in degrees of freedom $\left\{\delta \right\}$ and their time derivatives for the spatial discretization.

$\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\}$ (1)

The mass matrix $\left[M\right]$ is a constant coefficient nondiagonal matrix (often called consistent mass matrix as it is consistent with the local approximation of the displacement field). $\left[C_{\left[i\right]}\right]$ is the damping matrix corresponding to strain rate ${\epsilon }_{\left[i\right]}$ (material time derivative of order i of the Green’s strain tensor $\left[{\epsilon }_{\left[0\right]}\right]$ ) and $\left[K\right]$ is the stiffness matrix. Coefficients of $\left[C_{\left[i\right]}\right]$ and $\left[K\right]$ are up to quadratic functions of the displacement gradients. $\left\{F\right\}$ and $\left\{P\right\}$ are load vector and vector of secondary variables resulting due to integration by parts due to GM/WF used for spatial discretization. These nonlinear ODEs (1) correspond precisely to the mathematical model based on CBL of CCM. Evolution i.e., nonlinear dynamics resulting from the solutions of (1) without modifications or approximation are not available in the published literature. Modified forms of (1) based on many assumptions and approximations are obtained and are used to study nonlinear dynamics. We discuss one such approach that is common in many published works [30] [31] in the following. It can be shown using Equation (115) in reference [6] that $\left[C_{\left[i\right]}\right]$ and $\left[K\right]$ can be additively decomposed into (i.e., consist of):

$\left[C_{\left[i\right]}\right]=\left[C_{\mathrm{<mo\; stretchy="false">\; [\; </mo>}i\mathrm{<mo\; stretchy="false">\; ]\; </mo>}}^1\right]+\left[C_{\left[i\right]}^2\right]+\left[C_{\left[i\right]}^3\right]$ (2)

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

in which coefficients of $\left[C_{\left[i\right]}^1\right]\mathrm{,}\left[K^1\right]$ are constant, while the 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[\frac{\partial \left\{u\right\}}{\partial \left\{x\right\}}\right]$ . By substituting (2) and (3) in (1) and regrouping terms, we can obtain:

$\begin{array}{l}\left[M\right]\left\{\ddot{\delta }\right\}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[C_{\mathrm{<mo\; stretchy="false">\; [\; </mo>}i\mathrm{<mo\; stretchy="false">\; ]\; </mo>}}^1\right]\left\{{\delta }_{\left[i\right]}\right\}+\left[K^1\right]\left\{\delta \right\}\\ =\left\{F\right\}+\left\{P\right\}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[\left[C_{\left[i\right]}^2\right]+\left[C_{\left[i\right]}^3\right]\right]\left\{{\delta }_{\left[i\right]}\right\}-\left[\left[K^2\right]+\left[K^3\right]\right]\left\{\delta \right\}.\end{array}$ (4)

If we choose $\left\{\underset{˜}{F}\left(\left\{\delta \right\}\right)\right\}$ to define the right hand side of (4) and choose $n=1$ with $\left[C_{\left[1\right]}^1\right]=\left[C\right],\left[K^1\right]=\left[K\right],\left\{{\delta }_{\left[1\right]}\right\}=\left\{\stackrel{\dot{}}{\delta }\right\}$ , then (4) is written as:

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

In (5), the mass, damping and stiffness matrices have constant coefficients. All nonlinearities are lumped into $\left\{\underset{˜}{F}\left(\left\{\delta \right\}\right)\right\}$ in (5). While (5) is equivalent to (4), it is not possible to clearly ascertain the factors that influence nonlinear dynamics solely from (5) without the knowledge and further assessment of (4). In (5), the nonlinear damping and stiffness physics are lumped together in $\left\{\underset{˜}{F}\right\}$ , whereas in (4), both linear and nonlinear stiffness and damping terms are clearly identifiable. We remark that even though the $\left[M\right]\mathrm{,}\left[C\right]$ and $\left[K\right]$ in (5) have constant coefficients, the $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ in (5) are same as those in (4), hence they correspond to the solution of (4) i.e. are due to nonlinear system of ODEs in time.

We note that the mass, stiffness and damping matrices in (5) are non diagonal matrices. To the best of our knowledge, we have not found published works in which the nonlinearities in $\left\{F\right\}$ are clearly identified as we have done in the subsequent sections of this paper. Additionally, there are no reported solutions of (5) in which $\left[M\right]\mathrm{,}\left[C\right]$ and $\left[K\right]$ are all nondiagonal and are not altered using further assumptions. Instead, we find that the ODEs in (5) are reduced to a decoupled system of ODEs in time using the following approach.

In (5), $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ correspond to nonlinear ODEs in time. We consider linear undamped eigenvalue problem corresponding to linear ODEs in time describing linear dynamics with damping.

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

Let the columns of $\left[\Phi \right]$ represent the mass normalized eigenvectors of the undamped eigenvalue problem corresponding to (6), then (6) can be transformed to modal basis using

${\left\{\delta \right\}}_l=\left[\Phi \right]\left\{x\right\}$ (7)

and ${\left\{\stackrel{\dot{}}{\delta }\right\}}_l=\left[\Phi \right]\left\{\stackrel{\dot{}}{x}\right\}$ (8)

${\left\{\ddot{\delta }\right\}}_l=\left[\Phi \right]\left\{\ddot{x}\right\}\mathrm{.}$ (9)

This procedure is standard [32] . However, transformation to modal basis using (7)~(9) cannot be applied to (5) as $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ in (5) are not the same as ${\left\{\delta \right\}}_l\mathrm{,}{\left\{\stackrel{\dot{}}{\delta }\right\}}_l$ and ${\left\{\ddot{\delta }\right\}}_l$ in (6) and (7)~(9). We can only apply modal basis transformation to the following nonlinear PDEs:

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

Clearly, (10) are not the same as (5), hence cannot describe nonlinear dynamic behavior in the same manner as (5). Another important and inherent assumption in (10) is that the response ${\left\{\delta \right\}}_l\mathrm{,}{\left\{\stackrel{\dot{}}{\delta }\right\}}_l$ and ${\left\{\ddot{\delta }\right\}}_l$ is due to superposition of linear modes of vibrations (as evident from (7)~(9)). Superposition only holds for linear system, hence certainly cannot be used for the nonlinear dynamics under consideration. In the published works (7)~(9) are substituted in (10) and (10) is premultiplied by ${\left[\Phi \right]}^T$ to obtain the following ODEs in time.

${\left[\Phi \right]}^T\left[M\right]\left[\Phi \right]\left\{\ddot{x}\right\}+{\left[\Phi \right]}^T\left[C\right]\left[\Phi \right]\left\{\stackrel{\dot{}}{x}\right\}+{\left[\Phi \right]}^T\left[K\right]\left[\Phi \right]\left\{x\right\}={\left[\Phi \right]}^T\left\{\underset{˜}{F}\left({\left\{\delta \right\}}_l\right)\right\}$ (11)

If ${\omega }_i^2$ are the eigenvalues of the eigenvalue problem with linear $\left[K\right]$ and $\left[M\right]$ , then (11) reduces to

$\left[I\right]\left\{\ddot{x}\right\}+{\left[\Phi \right]}^T\left[C\right]\left[\Phi \right]\left\{\stackrel{\dot{}}{x}\right\}+\left[\Lambda \right]\left\{x\right\}={\left[\Phi \right]}^T\left\{F\left({\left\{\delta \right\}}_l\right)\right\}$ (12)

in which $\left[I\right]$ identity and $\left[\Lambda \right]$ is a diagonal matrix containing ${\omega }_1^2\mathrm{,}{\omega }_2^2\mathrm{,}\cdots \mathrm{,}{\omega }_n^2$ on the diagonal.

We note that in the modal basis, ${\left[\Phi \right]}^T\left[C\right]\left[\Phi \right]$ is nondiagonal as is the matrix $\left[C\right]$ in (7). Rayleigh damping is used to diagonalize the matrix ${\left[\Phi \right]}^T\left[C\right]\left[\Phi \right]$ . In the transformation (7)~(9) to modal basis, $\left\{x\right\}\mathrm{,}\left\{\stackrel{\dot{}}{x}\right\}$ and $\left\{\ddot{x}\right\}$ are modal participation factors. This implies the actual response $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ of the nonlinear dynamics problem consists of the superposition of the mass normalized eigenvectors of the linear problem using modal participation factors $\left\{x\right\}\mathrm{,}\left\{\stackrel{\dot{}}{x}\right\}$ and $\left\{\ddot{x}\right\}$ as in (7)~(9). We discuss this in more details in the following. It is assumed that since equation (12) is in modal basis, the mechanism of dissipation must be related to the natural modes of vibrations (this assumption only holds if ODEs are linear). This leads to consideration of Rayleigh damping. We assume that $\left[C\right]$ is a linear combination of $\left[M\right]$ and constant coefficient $\left[K\right]$ i.e.

$\left[C\right]=\alpha \left[M\right]+\beta \left[K\right]$ (13)

where $\alpha$ and $\beta$ are functions of ${\omega }_i$ (natural frequencies of the linear system). Transforming (13) into modal basis, we have:

${\left[\Phi \right]}^T\left[C\right]\left[\Phi \right]=\alpha {\left[\Phi \right]}^T\left[M\right]\left[\Phi \right]+\beta {\left[\Phi \right]}^T\left[K\right]\left[\Phi \right]=\alpha \left[I\right]+\beta \left[\Lambda \right].$ (14)

We note that right hand side of (14) is a diagonal matrix containing $\alpha +\beta {\omega }_i^2;i=1,2,\cdots ,n$ on the diagonals. Thus, ${\left[\Phi \right]}^T\left[C\right]\left[\Phi \right]$ has been diagonalized using (13) in modal basis (14).

If we assume that ${\zeta }_i$ is the dimensionless damping coefficient corresponding to frequency ${\omega }_i$ , then based on 1D spring-mass-damper system, we can write

${\left[\Phi \right]}^T\left[C\right]\left[\Phi \right]=\alpha \left[I\right]+\beta \left[\Lambda \right]=\left[2{\zeta }_i{\omega }_i\left[I\right]\right].$ (15)

In (15), the diagonal entries of $\left[C\right]$ are given by $2{\zeta }_i{\omega }_i;i=1,2,\cdots ,n$ in the modal basis. Thus, we have:

$\alpha +\beta {\omega }_i^2=2{\zeta }_i{\omega }_i$ (16)

Equation (16) can be used to obtain $\alpha \mathrm{,}\beta$ from ${\omega }_i$ versus ${\zeta }_i$ experimental data (only two data points are needed to determine $\alpha$ and $\beta$ ). By substituting (15) into (11), we obtain

$\left[I\right]\left\{\ddot{x}\right\}+\left[2{\zeta }_i{\omega }_i\left[I\right]\right]\left\{\stackrel{\dot{}}{x}\right\}+\left[\Lambda \right]\left\{x\right\}={\left[\Phi \right]}^{\text{T}}\left\{\underset{˜}{F}\left(\left\{\delta \right\}\right)\right\}=\left\{Q\left(\left\{x\right\}\right)\right\}$ (17)

in which $\left[2{\zeta }_i{\omega }_i\left[I\right]\right]$ is a diagonal damping matrix containing $2{\zeta }_i{\omega }_i;i=1,2,\cdots ,n$ on its diagonal.

Equations (17) represent a decoupled system of “n” second order nonlinear ODEs in time. We can also express (17) as follows:

${\ddot{x}}_i+2{\zeta }_i{\omega }_i+{\omega }_i^2{x}_i=Q_i\left(\left\{x\right\}\right);i=1,2,\cdots ,n$ (18)

1.4. Remarks

In arriving at (18) from (5) there are many assumptions involved.

(1) The nonlinear ODEs cannot be transformed to modal basis as there is no modal basis for nonlinear ODEs. We can only transform (6) to modal basis as in this case ${\left\{\delta \right\}}_l\mathrm{,}{\left\{\stackrel{\dot{}}{\delta }\right\}}_l\mathrm{,}{\left\{\ddot{\delta }\right\}}_l$ correspond to linear case. However, (6) does not represent the nonlinear dynamics physics described by (5).

(2) Decoupled ODEs in time (18) are exactly same as they appear in linear dynamics if the nonhomogeneous part is not a function of $\left\{x\right\}$ . This suggests that calculating modal participation factors $x_i$ using nonlinear decoupled ODEs and then using them in (7)~(9) is sufficient for (7)~(9) to yield $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ of nonlinear problem instead of ${\left\{\delta \right\}}_l\mathrm{,}{\left\{\stackrel{\dot{}}{\delta }\right\}}_l\mathrm{,}{\left\{\ddot{\delta }\right\}}_l$ of linear problem.

(3) Use of superposition for a nonlinear problem raises the question regarding the validity of the resulting solutions $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ obtained using (18) and (7)~(9).

(4) Use of Rayleigh damping itself is highly questionable. We recall that in deriving nonlinear ODEs (1), damping or dissipation is considered to be dependent on strain rates in the constitutive theory for contravariant deviatoric second Piola-Kirchhoff stress tensor that leads to damping matrices $\left[C_{\left[i\right]}\right]$ in (1). Assuming that $\left[C_{\left[1\right]}\right]$ or $\left[C\right]$ is proportional to $\left[M\right]$ and $\left[K\right]$ in (13) is heuristic and has no physical basis and often leads to wrong evolution and the stationary state. Therefore, the validity of the resulting ODEs in time (18) is highly questionable. Surana et al. [33] [34] have shown that solutions obtained using (18) lead to incorrect stationary states in cases where the evolution possess a stationary state, solution of the corresponding boundary value problem.

(5) In decoupled system (18), each ODE is precisely a single degree of freedom spring, mass, damper system. Thus, it is perhaps not surprising that superposition of the solutions obtained from (18) using (7)~(9) exhibit similar characteristics as a single mass, spring, damper system.

(6) We point out that exact form of $\left\{Q\left(\left\{\delta \right\}\right)\right\}$ can be quite complex. What form is exactly used in the reported solution of (18) is not available in the published works.

(7) Regardless of many other details, it is conclusive that this approach described above (used dominantly in published works) cannot be used to study nonlinear dynamic response of continuous system with finite deformation, finite strain and consistent mechanism of dissipation.

2. Scope of Work

Based on the remarks in Section 1, our view is that the true mathematical description of the nonlinear dynamics physics for continuous TVE solids without memory is given by CBL of CCM derived using contravariant second Piola-Kirchhoff stress tensor, Green strain tensor, its rates and consistent constitutive theories derived using entropy inequality and representation theorem. The resulting system of nonlinear PDEs can be converted into a nonlinear system of second order ODEs (1) by using space-time decoupled finite element processes with GM/WF in space for the spatial discretization [6] . Therefore, the solutions of original nonlinear PDEs or the solution of the nonlinear ODEs (1) is the solution that describes the deformation physics of nonlinear dynamics in TVE solids with finite deformation physics and nonlinear damping. In this paper, we present such solutions.

When TVE solids are subjected to harmonic excitation resulting in finite deformation, finite strain deformation physics, we observe (also reported in many published works) the following in the frequency response:

(1) Progressively increasing frequency of excitation yields a frequency response curve in which there could be a sudden change in the maximum amplitude of motion at or in the neighborhood of a specific frequency.

(2) Progressively decreasing frequency of excitation also yields a frequency response curve in which there could also be a sudden change in the maximum amplitude of motion at or in the neighborhood of a specific frequency.

(3) First, we note that (1) and (2) do not always exist for all choice of parameters that influence the dynamic response. However, if (1) and (2) do exist, then the two frequency response curves generally differ except in the range of very low and very high frequencies. Secondly, the occurrence of jump (decrease or increase) in the maximum amplitude in (1) and (2) generally does not happen at the same frequency.

(4) These frequency response curves and associated backbone curves constructed using the information in these, have been reported for both single mass and two mass spring, dashpot systems with nonlinear stiffness and damping as well from the solutions of decoupled ODEs as in (10) but transformed to (18) using linear modal basis.

(5) The work in this paper focused on many basic aspects of nonlinear dynamics including identification of the physics in nonlinear dynamics that can result in dynamic bifurcation phenomenon in the frequency response.

(6) It is shown that using thermodynamically and mathematically consistent system of PDEs in space and time resulting from the CBL of CCM and the consistent constitutive theories, the factors influencing nonlinear dynamics cannot be illustrated explicitly, though these are implicitly present in the nonlinear PDEs.

(7) The nonlinear system of ODEs resulting from space-time decoupled finite element method is the first step in identifying some of the factors influencing nonlinear dynamic response.

(8) It is shown that consideration of finite deformation, finite strain physics in corresponding BVPs and the solution of resulting nonlinear algebraic equation using Newton’s linear method is essential in explicitly determining some aspects of the physics that are implicitly present in item (7).

(9) The solution of the nonlinear ODEs in item (7) (i.e. Equations (1)) are calculated using Newmark linear acceleration method with Newton’s linear method for each increment of time to obtain a converged solution. Use of Newton’s linear method reveals that incremental solution calculation procedure in the iterative process contains some of the same physics as in BVP, hence the same mathematical form associated with the discretization as in case of BVP discussed in (8). This facilitates identification of some aspects of the physics and their influence on dynamic bifurcation phenomenon in nonlinear dynamics.

(10) Solutions of model problems are presented using space-time coupled finite element method. Space time domain is discretized into space-time strips, each corresponding to increments in time. The first space-time strip (for $0\le t\le \Delta t$ ) is discretized using p-version hierarchical space-time finite elements with higher order global differentiability in space and time. PDEs in the mathematical model are used to construct space-time variationally consistent space-time integral form based on space-time residual functional for the discretization of a space-time strip. Space-time variational consistency of the space-time residual functional (IVP) ensures unconditionally stable computations [32] . A converged solution is obtained for appropriate choice of h and p with minimally conforming approximation space. Residual functional values for the discretization of the order of $O\left({10}^{-8}\right)$ or lower, ensure accurate time evolution. Upon convergence of the solution for the first space-time strip, the remaining evolution is computed by using subsequent space-time strips and time marching. In this approach, the use of space-time variationally consistent space time coupled finite element method eliminates some approximations and other issues of convergence that exist in space-time decoupled methods. All model problem studies in the paper are presented using this approach. Model problems are designed in such a way that their solutions presented in this paper can be compared with some published results as well as with their 1D equivalent models using mass, 1D nonlinear spring and 1D nonlinear dashpot when possible.

3. Determination of Factors Influencing Nonlinear Dynamic Response

Complete mathematical model for finite deformation, finite strain nonlinear dynamics of thermoviscoelastic solids (TVES) without memory using contravariant second Piola-Kirchhoff stress tensor and Green’s strain tensor and using ordered rate constitutive theory for dissipation based on material derivatives of the Green’s strain tensor up to order n has been presented by Surana et al. in reference [6] . The derivation of this mathematical model is not presented in this paper for the sake of brevity, but the equations constituting the mathematical model are given in the following:

${\rho }_0=\left|J\right|\rho \left(x\mathrm{,}t\right)\left(\text{CM}\right)$ (19)

${\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]}^T\right]\{\nabla \}=0\left(\text{BLM}\right)$ (20)

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

${\rho }_0\frac{De}{Dt}+{\{\nabla \}}^T\left\{q\right\}-\text{tr}\left(\left[{\sigma }^{\left[0\right]}\right]\left[\stackrel{\dot{}}{\epsilon }\right]\right)=0\left(\text{FLT}\right)$ (22)

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

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

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

Constitutive theories for $\left[{}_d{\sigma }^{\left[0\right]}\right]$ and $\left\{q\right\}$ based on integrity

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

$\left\{q\right\}=-{{\kappa }_1|}_{{}_{\underset{\_}{\Omega }}}\left\{g\right\}-{{\kappa }_2|}_{{}_{\underset{\_}{\Omega }}}\left({\left\{g\right\}}^{\text{T}}\left\{g\right\}\right)\left\{g\right\}$ (27)

Simplified constitutive theories for $\left[{}_d{\sigma }^{\left[0\right]}\right]$ and $\left\{q\right\}$ (in Voigt’s notation) that are linear in tensors ( $\left[\epsilon \right]\mathrm{,}\left[{\epsilon }_{\left[i\right]}\right]\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2},\cdots \mathrm{,}n$ and $\left\{q\right\}$ ) are given by:

$\left[{}_d{\sigma }^{\left[0\right]}\right]={\underset{˜}{\sigma }}^0\left[I\right]+2\mu \left[{\epsilon }_{\left[0\right]}\right]+\lambda \text{tr}\left(\left[{\epsilon }_{\left[0\right]}\right]\left[I\right]\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}2{\mu }_i\left[{\epsilon }_{\left[i\right]}\right]+{\displaystyle \underset{i\mathrm{<mo>\; =\; </mo>1}}{\overset{n}{\sum }}}{\lambda }_i\text{tr}\left(\left[{\epsilon }_{\left[i\right]}\right]\left[I\right]\right)$ (28)

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

in which $\left\{u\right\}$ are displacements, $\left\{F^b\right\}$ are body forces per unit mass, ${\sigma }^{\left[0\right]}$ is contravariant second Piola-Kirchhoff stress tensor, e is internal energy density, $\left\{q\right\}$ is heat vector, ${}_e\sigma {}^{\left[0\right]}$ and ${}_d\sigma {}^{\left[0\right]}$ are equilibrium and deviatoric contravariant second Piola-Kirchhoff stress tensors. $p\left(\rho \mathrm{,}\theta \right)$ is equation of state, $p\left(\theta \right)$ is mechanical pressure. $\left[{}^{\sigma }\underset{˜}{G}{}^i\right]\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}N$ and ${}^{\sigma }\underset{˜}{I}{}^J\mathrm{<mo>\; ;\; </mo>}j=\mathrm{1,2,}\cdots \mathrm{,}M$ are combined generators and invariants of the argument tensors $\left[{\epsilon }_{\left[0\right]}\right]$ and $\left[{\epsilon }_{\left[i\right]}\right]\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,}n$ of the constitutive tensor $\left[{}_d{\sigma }^{\left[0\right]}\right]$ , $\left\{g\right\}$ is a temperature gradient tensor (argument tensor of constitutive tensor $\left\{q\right\}$ ), $\mu$ and $\lambda$ are elastic constants, ${\mu }_i$ and ${\lambda }_i$ are dissipation coefficients for strain rate $\left[{\epsilon }_{\left[i\right]}\right]$ . This mathematical model has closure. It consists of thirteen equations BLM(3), FLT(1), constitutive theories for $\left[{}_d{\sigma }^{\left[0\right]}\right]$ (6), $\left\{q\right\}$ (3) in thirteen dependent variables $u_i$ (3), $\left[{}_d{\sigma }^{\left[0\right]}\right]$ (6), $\left\{q\right\}$ (3), $\theta$ (1). We note that $e=e\left(\rho \mathrm{,}\theta \right)$ . $\underset{\_}{\Omega }$ is a known configuration in which all material coefficients are described. For details of the derivation of the PDEs in the mathematical model, the reader can refer to references [6] [15] .

If we consider the mathematical model consisting of nonlinear PDEs in the CBL of CCM and consistent constitutive theories presented above, then all we can infer is that: (1) Green’s strain permits finite deformation and finite strain, (2) the nonlinear dissipation mechanism (due to Green’s strain rates) converts some mechanical energy into entropy, (3) the inertial terms related to mass in general result in a reduction in the stiffness of the system, thereby resulting in larger displacements dynamically compared to when the same loads applied statically and (4) the presence of dissipation also provides resistance to motion.

On the other hand, if we consider the nonlinear ODEs (1) resulting from space-time decoupled finite element method (with GM/WF in space), then we can infer the following:

(1) Stiffness matrix is the addition of three matrices, namely $\left[K^1\right]\mathrm{,}\left[K^2\right]$ and $\left[K^3\right]$ . The coefficients of $\left[K^1\right]$ are constant while the coefficients of $\left[K^2\right]$ and $\left[K^3\right]$ are linear and quadratic in the gradients of displacements (and thereby degrees of freedom), respectively. However, the specific influence of each matrix on nonlinear dynamics is not readily apparent. Therefore, at this stage we can only say that nonlinear dynamic response is dependent on a constant stiffness matrix and stiffness matrices with coefficients that are linear and quadratic functions of displacement gradients. Thus, the total stiffness is nonlinear. Significant contributions of $\left[K^2\right]$ and $\left[K^3\right]$ are essential in nonlinear dynamic response. In order for these contributions to be substantial, the applied loads must be sufficiently large so that finite deformation exists, resulting in significant contribution of $\left[K^2\right]$ and $\left[K^3\right]$ . We also note that $\left[K^1\right]$ and $\left[K^3\right]$ are symmetric but $\left[K^2\right]$ is nonsymmetric.

(2) The damping matrices $\left[C_{\left[i\right]}\right]$ (due to ${\epsilon }_{\left[i\right]}$ ) are also composed of addition of $\left[C_{\left[i\right]}^1\right]\mathrm{,}\left[C_{\left[i\right]}^2\right]$ and $\left[C_{\left[i\right]}^3\right]$ . The coefficients of $\left[C_{\left[i\right]}^1\right]$ are constant, while those of $\left[C_{\left[i\right]}^2\right]$ and $\left[C_{\left[i\right]}^3\right]$ are linear and quadratic functions of the components of displacement gradient tensor. Matrices $\left[C_{\left[i\right]}^1\right]$ and $\left[C_{\left[i\right]}^3\right]$ are symmetric but $\left[C_{\left[i\right]}^2\right]$ is nonsymmetric.

(3) It is well understood that constant coefficient damping $\left[C_{\left[i\right]}^1\right]$ results in base elongation and amplitude decay of an applied disturbance, conversion of mechanical energy into entropy and offers resistance to motion. We cannot conclusively make similar statements for $\left[C_{\left[i\right]}^2\right]$ and $\left[C_{\left[i\right]}^3\right]$ but intuitively this seems plausible. Regardless of their specific form, damping matrix must result in some level of dissipation and resistance to motion. The damping material coefficients in $\left[C_{\left[i\right]}\right]$ are actual material coefficients from the constitutive theory. These play critical role, as discussed in a later section, when Equations (1) are compared with phenomenological mathematical models in which there is no physics, hence material coefficients are arbitrary.

(4) Third important aspect is due to presence of $\left[M\right]\left\{\ddot{\delta }\right\}$ term in (1). In a vibrating system with a low but progressively increasing frequency below the first fundamental frequency (as observed in a linear system), the maximum amplitude of motion also progressively increases. This indicates that the dynamic stiffness of the system is lower than the static stiffness, resulting in an increased amplitude of motion compared to the static response for the same applied force. In a nonlinear system, this physics may not be as straight forward as in case of linear system. However, it is reasonable to assume that there might be some similarities to the observations stated above.

(5) Is the physics discussed in (1)~(3) and present in (1) sufficient to understand if and when the dynamic bifurcation in nonlinear dynamics can exist? Our view is that there is some additional physics and some crucial findings in the solution procedure over and beyond (1)~(3) that play an extremely important and crucial role in the existence or lack of dynamic bifurcation in nonlinear dynamics. We note that sudden change (increase or decrease) in the amplitude of motion at or in the neighborhood of a frequency is akin to the phenomenon of static bifurcation. This is accompanied by a sudden release of energy leading to a sudden change in the stiffness of the system, hence a sudden jump in the deformed state of the matter. If the bifurcation occurs in frequency response, then the aspects discussed so far (items (1)~(3)) must be inherently present in the mathematical model (1). However, they cannot be explicitly observed in the mathematical model itself. In the following, we demonstrate additional physics present in (1) that plays significant role in dynamic bifurcation.

3.1. Static Bifurcation

Consider balance of linear momenta (equation (10) in ref [6] ) for BVPs (finite deformation, finite strain).

${\rho }_0\left\{F^b\right\}+\left[\left[J\right]\left[{\sigma }^{\left[0\right]}\right]\right]\{\nabla \}=0\forall x\in {\Omega }_x$ (30)

A finite element formulation based on integral form of (30), integrated over a discretization ${\bar{\Omega }}_x^T={\displaystyle \underset{e}{\cup }{\bar{\Omega }}_x^e}$ of ${\bar{\Omega }}_x$ using fundamental lemma of the calculus of variations followed by GM/WF yields (also see ref [6] ):

${\displaystyle \underset{e}{\sum }}{\displaystyle {\int }_{{\bar{\Omega }}^e}}{\left[B\right]}^T\left\{{\sigma }^{\left[0\right]}\right\}\text{d}\Omega =\left\{P\right\}+\left\{F\right\}\mathrm{.}$ (31)

The following details can be found in ref [6] :

Displacement gradient tensor $\left[{}^{d}J\right]$ is given by

$\left[{}^{d}J\right]=\left[\frac{\partial \left\{u\right\}}{\partial \left\{x\right\}}\right]=\left[\left\{g_1\right\}\left\{g_2\right\}\left\{g_3\right\}\right]\mathrm{<mo>\; ;\; </mo>}\left[J\right]=\left[I\right]+\left[{}^{d}J\right]$ (32)

in which

$\begin{array}{l}{\left\{g_1\right\}}^T=\left[u_{1,1},u_{2,1},u_{3,1}\right]\\ {\left\{g_2\right\}}^T=\left[u_{1,1},u_{2,1},u_{3,1}\right]\\ {\left\{g_3\right\}}^T=\left[u_{1,1},u_{2,1},u_{3,1}\right]\end{array}$ (33)

Let

${\left\{g\right\}}^T=\left[{\left\{g_1\right\}}^T\mathrm{,}{\left\{g_2\right\}}^T\mathrm{,}{\left\{g_3\right\}}^T\right]\mathrm{.}$ (34)

Using local approximations for $u$ , we can write:

$\left\{g\right\}=\left[G\right]\left\{{\delta }^e\left(t\right)\right\}\mathrm{.}$ (35)

We decompose Green’s strain tensor ${\epsilon }_{\left[0\right]}$ into linear ${\epsilon }_{\left[0\right]}^l$ and non-linear components ${\epsilon }_{\left[0\right]}^{nl}$

$\left\{{\epsilon }_{\left[0\right]}^l\right\}=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0.5& 0& 0.5& 0\\ 0& 0& 0.5& 0& 0& 0& 0.5& 0& 0\\ 0& 0.5& 0& 0.5& 0& 0& 0& 0& 0\end{array}\right]\left\{g\right\}$ (36)

or

$\left\{{\epsilon }_{\left[0\right]}^l\right\}=\left[H\right]\left\{g\right\}=\left[H\right]\left[G\right]\left\{{\delta }^e\left(t\right)\right\}$ (37)

$\left\{{\epsilon }_{\left[0\right]}^{nl}\right\}=\frac{1}{2}\left[\begin{array}{ccc}{\left\{g_1\right\}}^T& 0& 0\\ 0& {\left\{g_2\right\}}^T& 0\\ 0& 0& {\left\{g_3\right\}}^T\\ 0& {\left\{g_3\right\}}^T& {\left\{g_2\right\}}^T\\ {\left\{g_3\right\}}^T& 0& {\left\{g_1\right\}}^T\\ {\left\{g_2\right\}}^T& {\left\{g_1\right\}}^T& 0\end{array}\right]\left\{g\right\}=\frac{1}{2}\left[A_g\right]\left\{g\right\}$ (38)

Using (35) in (36)

$\left\{{\epsilon }_{\left[0\right]}^{nl}\right\}=\frac{1}{2}\left[A_g\right]\left[G\right]\left\{{\delta }^e\left(t\right)\right\}$ (39)

Thus,

$\left\{{\epsilon }_{\left[0\right]}\right\}=\left\{{\epsilon }_{\left[0\right]}^l\right\}+\left\{{\epsilon }_{\left[0\right]}^{nl}\right\}=\left[H\right]\left[G\right]\left\{{\delta }^e\left(t\right)\right\}+\frac{1}{2}\left[A_g\right]\left[G\right]\left\{{\delta }^e\left(t\right)\right\}$ (40)

or

$\begin{array}{c}\left\{{\epsilon }_{\left[0\right]}\right\}=\left[B^l\right]\left\{{\delta }^e\left(t\right)\right\}+\frac{1}{2}\left[B^{nl}\right]\left\{{\delta }^e\left(t\right)\right\}\\ =\left[\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right]\left\{{\delta }^e\left(t\right)\right\}\end{array}$ (41)

where

$\begin{array}{l}\left[B^l\right]=\left[H\right]\left[G\right]\hfill \\ \left[B^{nl}\right]=\left[A_g\right]\left[G\right]\hfill \end{array}\}$ (42)

$\left[B\right]=\left[B^l\right]+\left[B^{nl}\right]$ (43)

In (31), $\left\{P\right\}$ is a vector of secondary variables and $\left\{F\right\}$ is the load vector (due to body forces in case of (30)). We assume that $\left\{F\right\}$ is conservative and $\left\{{\delta }^e\right\}$ are nodal degrees of freedom for an element e.

Recall

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

$\left\{{}_d{\sigma }^{\left[0\right]}\right\}=\left[D\right]\left\{{\epsilon }_{\left[0\right]}\right\}+\left[D\right]\left[\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right]\left\{{\delta }^e\right\}$ (45)

Using (43) and (45) in (31) we can write

$\left\{\Psi \right\}=\left\{P\right\}+\left\{F\right\}={\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}^e}{\int }}{\left[\left[B^l\right]+\left[B^{nl}\right]\right]}^T\left[D\right]\left[\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right]\left\{{\delta }^e\right\}\text{d}\Omega =0.$ (46)

Expanding the integrand in (46), we can write

$\left\{\Psi \right\}=\left\{P\right\}+\left\{F\right\}-{\displaystyle \underset{e}{\sum }}\left[\left[K_e^1\right]+\left[K_e^2\right]+\left[K_e^3\right]\right]\left\{{\delta }^e\right\}$ (47)

where

$\left[K_e^1\right]={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{\left[B^l\right]}^T\left[D\right]\left[B^l\right]\text{d}\Omega$ (48)

$\left[K_e^2\right]={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left(\frac{1}{2}{\left[B^l\right]}^T\left[D\right]\left[B^{nl}\right]+{\left[B^{nl}\right]}^T\left[D\right]\left[B^l\right]\right)\text{d}\Omega$ (49)

$\left[K_e^3\right]={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\frac{1}{2}{\left[B^{nl}\right]}^T\left[D\right]\left[B^{nl}\right]\text{d}\Omega$ (50)

are the element stiffness matrices. $\left[K_e^1\right]$ is a constant coefficient matrix and the coefficients of $\left[K_e^2\right]$ and $\left[K_e^3\right]$ are linear and quadratic functions of the coefficients of $\left[{}^{d}J\right]$ , hence linear and quadratic in dofs $\left\{{\delta }^e\right\}$ . $\left[K_e^1\right]$ and $\left[K_e^3\right]$ are symmetric but $\left[K_e^2\right]$ is non-symmetric. As expected, these are exactly same as in the case of corresponding IVP [6] .

We can write (46) or (47) as

$\left\{\Psi \left(\left\{\delta \right\}\right)\right\}-\left\{P\right\}+\left\{F\right\}-\left[K\left(\left\{\delta \right\}\right)\right]\left\{\delta \right\}=0$ (51)

Equations (51) represent a system of nonlinear algebraic equations. In order to find a solution $\left\{\delta \right\}$ that satisfies (51), an iterative method must be employed. Specifically, Newton’s linear method can be used to obtain a solution of the nonlinear algebraic Equations (51).

Let ${\left\{\delta \right\}}_0$ be the initial solution (assumed or guess), then

$\left\{\Psi \left({\left\{\delta \right\}}_0\right)\right\}\ne 0.$ (52)

Let $\left\{\Delta \delta \right\}$ be a change in ${\left\{\delta \right\}}_0$ such that

$\left\{\Psi \left({\left\{\delta \right\}}_0+\left\{\Delta \delta \right\}\right)\right\}=0.$ (53)

We expand $\left\{\Psi \left(\mathrm{.}\right)\right\}$ in (53) in Taylor series about ${\left\{\delta \right\}}_0$ and retain only upto linear terms in $\left\{\Delta \delta \right\}$ .

$\left\{\Psi \left({\left\{\delta \right\}}_0\right)+\left\{\Delta \delta \right\}\right\}=\left\{\Psi \left({\left\{\delta \right\}}_0\right)\right\}+{\left[\frac{\partial \left\{\Psi \right\}}{\partial \left\{\delta \right\}}\right]}_{{\left\{\delta \right\}}_0}\left\{\Delta \delta \right\}=0$ (54)

From (54), we can determine $\left\{\Delta \delta \right\}$

$\left\{\Delta \delta \right\}=-{\left[\frac{\partial \left\{\Psi \right\}}{\partial \left\{\delta \right\}}\right]}_{{\left\{\delta \right\}}_0}^{-1}\left\{\Psi \left({\left\{\delta \right\}}_0\right)\right\}=-{\left[\delta \left\{\Psi \right\}\right]}_{{\left\{\delta \right\}}_0}^{-1}\left\{\Psi \left({\left\{\delta \right\}}_0\right)\right\}.$ (55)

Improved solution $\left\{\delta \right\}$ is obtained using

$\left\{\delta \right\}={\left\{\delta \right\}}_0+\left\{\Delta \delta \right\}\mathrm{.}$ (56)

If

$\left|{\Psi }_i\left(\left\{\delta \right\}\right)\right|\le \Delta$ (57)

then the solution $\left\{\delta \right\}$ is considered the converged solution. Otherwise, we set ${\left\{\delta \right\}}_0=\left\{\delta \right\}$ and repeat the iteration process from (55)~(57). $\Delta$ is a preset tolerance for the computed zero.

Surana et al. [6] have shown that variation of $\left\{\Psi \left(\left\{\delta \right\}\right)\right\}$ leads to

$\delta \left\{\Psi \left(\left\{\delta \right\}\right)\right\}=-{\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}^e}{\int }}\left({\left[B\right]}^T\left[D\right]\left[B\right]+{\left[G\right]}^T\left[S\right]\left[G\right]\right)\text{d}\Omega =-{\displaystyle \underset{e}{\sum }}\left[K_T^e\right]=-\left[K_T\right].$ (58)

Substituting for $\left[B\right]$ in (58) from (43)

$\begin{array}{c}\delta \left\{\Psi (.)\right\}=-{\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}^e}{\int }}(\left[B^l\right]\left[D\right]\left[B^l\right]+\left({\left[B^l\right]}^T\left[D\right]\left[B^{nl}\right]+{\left[B^{nl}\right]}^T\left[D\right]\left[B^l\right]\right)\\ +{\left[B^{nl}\right]}^T\left[D\right]\left[B^{nl}\right]+{\left[G\right]}^T\left[S\right]\left[G\right])\text{d}\Omega \\ =\left[K_T\right]={\displaystyle \underset{e}{\sum }}\left[K_e^T\right].\end{array}$ (59)

Let

$\left[K_e^1\right]={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{\left[B^l\right]}^T\left[D\right]\left[B^l\right]\text{d}\Omega ;\left[K^1\right]={\displaystyle \underset{e}{\sum }}\left[K_e^1\right]$ (60)

$\left[{\underset{˜}{K}}_e^2\right]={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left({\left[B^l\right]}^T\left[D\right]\left[B^{nl}\right]+{\left[B^{nl}\right]}^T\left[D\right]\left[B^l\right]\right)\text{d}\Omega ;\left[{\underset{˜}{K}}^2\right]={\displaystyle \underset{e}{\sum }}\left[{\underset{˜}{K}}_e^2\right]$ (61)

$\left[{\underset{˜}{K}}_e^3\right]={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{\left[B^{nl}\right]}^T\left[D\right]\left[B^{nl}\right]\text{d}\Omega ;\left[{\underset{˜}{K}}^3\right]={\displaystyle \underset{e}{\sum }}\left[{\underset{˜}{K}}_e^3\right]$ (62)

$\left[K_{\sigma }^e\right]={\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{\left[G\right]}^T\left[\sigma \right]\left[G\right]\text{d}\Omega ;\left[K_{\sigma }\right]={\displaystyle \underset{e}{\sum }}\left[K_{\sigma }^e\right]$ (63)

$\left[K_T\right]\mathrm{,}\left[K_e^1\right]\mathrm{,}\left[{\underset{˜}{K}}_e^2\right]\mathrm{,}\left[{\underset{˜}{K}}_e^3\right]$ and $\left[K_{\sigma }^e\right]$ are all symmetric. $\left[K_e^T\right]$ is called the tangent stiffness matrix for element “e”. Assembly of the element equations follows the usual procedure. $\left[K_T\right]$ is called tangent stiffness matrix for the discretization.

Remarks

(1) In the neighborhood of static bifurcation (load at which amplitude of motion experiences jump), $\left[K_T\right]$ is momentarily singular causing instability. At bifurcation load, $\left[K_1\right]+\left[{\underset{˜}{K}}^2\right]+\left[{\underset{˜}{K}}^3\right]+\left[K_{\sigma }\right]=\left[0\right]$ ). That is, at static bifurcation point, compressive or negative $\left[K_{\sigma }\right]$ is exactly equal to the sum of $\left[K^1\right]$ , ${\left[\underset{˜}{K}\right]}^2$ and $\left[{\underset{˜}{K}}^3\right]$ . Thus, significance of $\left[K_{\sigma }\right]$ is now rather obvious.

(2) We note that in calculating the increment solution of $\left\{\Delta \delta \right\}$ using (55) and (59), there are four stiffness matrices involved. $\left[K_e^1\right]\mathrm{,}\left[{\underset{˜}{K}}_e^2\right]$ and $\left[{\underset{˜}{K}}_e^3\right]$ have constant coefficients, linear and quadratic coefficients of the displacement gradient tensor (or $\left\{{\delta }^e\right\}$ , hence $\left\{\delta \right\}={\displaystyle \underset{e}{\cup }{\delta }^e}$ ). While $\left[K_e^1\right]$ is same as in (47) but $\left[{\underset{˜}{K}}_e^2\right]$ and $\left[{\underset{˜}{K}}_e^3\right]$ are not the same as $\left[K_e^2\right]$ and $\left[K_e^3\right]$ in (49) and (50).

(3) We also note that in $\left[{\underset{˜}{K}}_e^2\right]$ and $\left[{\underset{˜}{K}}_e^3\right]$ , the dependence of the coefficients on the gradients of displacement (i.e., the degree) is same as in $\left[K_e^2\right]$ and $\left[K_e^3\right]$ .

(4) The matrix $\left[K_{\sigma }^e\right]$ (explicitly not present in (47)) is called stress stiffness matrix. Its presence in (59) indicates that the presence of stress field influences the stiffness of the deforming matter. When the stress field is tensile (positive), the resulting $\left[K_{\sigma }^e\right]$ likewise is positive, hence results in increased stiffness of the deforming volume of matter. On the other hand, a compressive stress field (negative) will result in a negative $\left[K_{\sigma }^e\right]$ , causing decrease in stiffness of the deforming volume of matter.

(5) In static and dynamic bifurcation physics, $\left[K_{\sigma }\right]={\displaystyle {\sum }_e}\left[K_{\sigma }^e\right]$ plays a crucial role in the mechanism of sudden energy release that leads to bifurcation. Without $\left[K_{\sigma }\right]$ , bifurcation cannot be observed in the nonlinear dynamics. Additionally, from the derivation of $\left[K_{\sigma }\right]$ , we observe that $\left[K^{\sigma }\right]$ only exists when finite deformation, finite strain deformation physics is present in the deforming matter.

(6) We shall see that the matrices (60)~(63) appear explicitly when integrating nonlinear ODEs (1) in time using Newmark linear acceleration method with Newton’s linear method for a time interval to obtain converged solution, confirming that stress stiffness matrix $\left[K_{\sigma }\right]$ also controls dynamic bifurcation physics in nonlinear dynamics. Static bifurcation cannot exist without $\left[K_{\sigma }\right]$ , hence deviatoric second Piola-Kirchhoff stress and Green’s strain as a conjugate pair are essential in the mathematical model.

3.2. Factors Influencing Dynamic Bifurcation Phenomenon in Nonlinear Dynamics

First, we note that factual determination of the factors influencing dynamic bifurcation can only be done when we determine how the incremental solution for each increment of time is computed. The coefficient matrix in the computations must contain all of the details related to the evolution of the dynamic problem.

In this section we list various aspects of the physics in the mathematical model that we have observed so far that are likely to influence nonlinear dynamics in general, and in particular, the existence of bifurcation phenomenon in the frequency response.

(1) Finite deformation, finite strain with dissipation is irreversible physics that leads to path (loading path) dependent response. This, of course, is present due to the stiffness matrices $\left[K^1\right]\mathrm{,}\left[K^2\right]$ and $\left[K^3\right]$ , as well as due to similar matrices for the dissipation mechanism. The mechanism of increase or decrease in stiffness resulting due to increasing or decreasing deformation and strain is present in the stiffness matrices $\left[K^1\right]\mathrm{,}\left[K^2\right]\mathrm{,}\left[K^3\right]$ although quantifying it precisely is not trivial. At this stage, all we can conclude is that these factors are essential for nonlinear dynamics and the existence of bifurcation physics.

(2) Dissipation provides resistance to motion, thereby contributes additional stiffness to the system. Increasing damping material coefficients increases stiffness of the system and vice versa. Thus, damping physics influences dynamic bifurcation physics. Reduced damping physics enhances the likelihood of the presence of dynamic bifurcation and vice versa. The presence of dissipation results in conversion of some mechanical work into entropy, thus heat, resulting in change of thermal field. In case of repeated cyclic loading over long periods of time, the thermal effects may need to be considered. In the present work, we only consider isothermal physics. Non-isothermal nonlinear dynamics investigation will be considered in a follow on paper.

(3) A tensile stress field results in positive $\left[K_{\sigma }\right]$ , which enhances the total stiffness of the medium. On the other hand, a compressive stress field results in negative $\left[K_{\sigma }\right]$ which reduces the total stiffness of the medium. Thus $\left[K_{\sigma }\right]$ plays a major and significant role in the static as well as dynamic bifurcation physics.

(4) As discussed earlier, translational inertial physics due to $\left[M\right]\left\{\ddot{\delta }\right\}$ results in reduction of the stiffness of the medium, thus enhancing the likelihood of the existence of dynamic bifurcation physics.

(5) Thus, nonlinear stiffness, linear and nonlinear damping, change in the stiffness due to $\left[K_{\sigma }\right]$ and the change in stiffness due to $\left[M\right]\left\{\ddot{\delta }\right\}$ all play a role in dynamic response. We note that negative $\left[K_{\sigma }\right]$ and $\left[M\right]\left\{\ddot{\delta }\right\}$ lead to a reduction in stiffness, while increasing dissipation and increasing nonlinearity contribute to an increase in stiffness. Therefore, when generating frequency response curves in the nonlinear deformation range, a balance between these two mechanisms of stiffness increase and decrease exist at every frequency when the motion becomes cyclic. A bifurcation point indicates an imbalance between these mechanisms, resulting in a sudden change in the energy and hence a sudden change in the amplitude of the motion to restore the balance. When $\left[K_{\sigma }\right]$ is always positive, dynamic bifurcation will require the minimum possible damping (to reduce stiffness) and enough decrease in stiffness due to $\left[M\right]\left\{\ddot{\delta }\right\}$ . On the other hand, when $\left[K_{\sigma }\right]$ is negative, $\left[K_{\sigma }\right]$ along with some decrease in stiffness due to $\left[M\right]\left\{\ddot{\delta }\right\}$ may result in enough decrease in total stiffness to allow for the dynamic bifurcation to exist. We point out that existence of dynamic bifurcation physics for some damping can possibly be eliminated by increasing damping, which leads to increase in stiffness, and vice versa. A sudden change in the state of energy of the deforming volume of matter at a frequency is necessary for the static as well as dynamic bifurcation to occur. Thus, $\left[K_{\sigma }\right]$ plays a crucial role in bifurcation phenomenon. In the absence of $\left[K_{\sigma }\right]$ (corresponding to infinitesimal deformation physics), bifurcation is not possible. Presence of $\left[K_{\sigma }\right]$ requires finite deformation physics, which is essential for the occurrence of static or dynamic bifurcation. This highlights the significance of the role played by $\left[K_{\sigma }\right]$ in enabling or contributing to static or dynamic bifurcation phenomena. When finite deformation nonlinearities are insignificant, either bifurcation cannot occur or its existence may require extremely low values of dissipation and significant weakening of stiffness due to translational inertial physics ( $\left[M\right]\left\{\ddot{\delta }\right\}$ ).

(6) In a continuous system where the response at a material point is influenced by the neighboring material points, the physics of evolution is dependent on position and time. This physics cannot be described by a lumped system that is independent of position. We discuss this aspect in more detail in a later section.

(7) We summarize the main points in the following:

(a) For bifurcation to exist, the deformation must be finite in order for the nonlinearities in stiffness and damping to be active in the dynamic response.

(b) Magnitude of the applied force (or disturbance) and dissipation are two key elements for the dynamic bifurcation physics in the nonlinear deformation range.

(c) For a given magnitude of the force (when deformation is nonlinear), there is a threshold value of dissipation for the existence of dynamic bifurcation. Below this threshold value of dissipation, dynamic bifurcation is always likely, while a damping value above the threshold value can eliminate the existence of dynamic bifurcation.

(d) Similarly, for a given damping value, there is a threshold value of the magnitude of applied force (disturbance) for the existence of dynamic bifurcation phenomenon. A force value below this threshold value, the dynamic bifurcation can be eliminated. But for forces higher than threshold value, dynamic bifurcation is likely to exist.

(e) The presence of $\left[K_{\sigma }\right]$ (which only exists for finite deformation, finite strain physics) plays crucial role in the static as well as dynamic bifurcation physics, which is always associated with sudden change in energy state and as a result, a sudden change in $\left[K_{\sigma }\right]$ . Presence of positive $\left[K_{\sigma }\right]$ and presence of negative $\left[K_{\sigma }\right]$ leads to dynamic bifurcation curves in the frequency response with different characteristics.

(f) When dynamic bifurcation exists, frequency response is generally path dependent i.e., progressively increasing frequencies and progressively decreasing frequencies will result in different frequency response curves, with the presence of bifurcation at different frequencies.

(g) In a following section we present derivation of the time integration of the nonlinear ODEs resulting from the space-time decoupling in which computations of the incremental solution for a time step contains tangent matrix with all necessary details of the factors influencing dynamic bifurcation.

4. Mathematical Models

In the following two sections, we consider a 1D phenomenological mathematical model used quite frequently in published works and a mathematical model based on CBL of CCM.

4.1. 1D Phenomenological Mathematical Models versus Mathematical Models Based on CBL of CCM

There is a fundamental question that we must address: can 1D phenomenological mathematical models with nonlinear damping and nonlinear stiffness ever replicate the nonlinear dynamic response of a TVE continuous solid media? This is the simplest possible exercise to determine if the phenomenological mathematical models reported in the literature are of any relevance in regard to the nonlinear dynamic response of TVE continuous matter.

Consider the one dimensional phenomenological mathematical model (used in many published works) with single degree of freedom u, consisting of single mass, nonlinear spring and a nonlinear dashpot in which stiffness and damping coefficients are up to quadratic functions of u [16] [25] [26] [30] (later referred to as model problem I).

$m\ddot{u}+\left(c_1+c_{2}u+c_3{u}^2\right)\stackrel{\dot{}}{u}+\left(k_1+k_{2}u+k_3{u}^2\right)u=f_0\sin \left(\omega t\right)$ (64)

In (64), m is mass, $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3$ are damping coefficients associated with constant, linear and quadratic damping in u. Likewise, $k_1\mathrm{,}{k}_2$ and $k_3$ are similar stiffness coefficients.

Comparing (64) and (1), we note the following:

(a) Damping coefficients $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3$ and stiffness coefficients $k_1\mathrm{,}{k}_2\mathrm{,}{k}_3$ are not material coefficients in (64). These are a result of the products of the material coefficients for damping and elasticity with the coefficients of damping and stiffness matrices. Thus, in this model, we have no idea about the definition of the material coefficients and their values as well as the coefficients in the damping and stiffness matrices that arise due to nonlinear stiffness and dissipation.

(b) In the mathematical model (1), each strain rate damping requires two material coefficients in $\left[D_i\right]$ which remain the same for constant, linear, and quadratic damping in displacement gradients. Same is true for stiffness. That is, in (64) $k_1\mathrm{,}{k}_2$ and $k_3$ are three coefficients associated with constant, linear and quadratic stiffnesses in u, where as in model (1), only two material coefficients $2\underset{˜}{\mu }\mathrm{,}\underset{˜}{\lambda }$ are needed, which remain the same for constant, linear and quadratic stiffnesses in displacement gradients. It is noteworthy that in (64), linear stiffness and linear damping only require a single material coefficient each, represented by $k_1$ and $c_1$ , respectively (due to one dimensional nature of the model). But in the case of (1), two material coefficients, $2\eta$ and $\lambda$ , are needed for linear stiffness as well as two material coefficients, $2{\eta }_i$ and ${\kappa }_i$ , are needed for linear damping.

(c) A mathematical model that is based on CBL of CCM that can possibly be viewed closed to (64) will be mathematical model based on CBl of CCM for 1D axial rod fixed at left end and subjected to harmonic excitation at the right end. In this case, the mathematical model is a reduction of 3D physics in (1) to 1D physics, hence contains material coefficients that are different than in model (1). In this model based on CBL of CCM, deformation intrinsically depends upon space and time. On the other hand, (64) represents a 1D lumped phenomenological model that is insensitive to spatial position. Does the phenomenological model (64) describe physics of nonlinear dynamics of the axial rod subjected to harmonic excitation? Other than this fundamental difference between (64) and (1), we have also pointed out enough other differences between the two mathematical models to conclude that (64) is not adequate to describe dynamic response of any continuous solid media. We present numerical studies in the following section to illustrate many points discussed here.

(d) In the mathematical model given by (64), stiffness and damping are up to quadratic functions of degrees of freedom (u), $\left[K_{\sigma }\right]$ physics must be intrinsic in the quadratic terms (if (64) describes bifurcation physics) but cannot be explicitly extracted from it as we have shown in case of model (I) based on CBL of CCM.

(e) Solution of (64) is obtained using Newmark linear acceleration method with Newton’s linear method for each time increment.

4.2. Mathematical Models Based on CBL of CCM (Model Problem II) and Solution Methods

The solutions of the PDEs in the mathematical model describing evolution derived using CBL of CCM can be obtained using different methods of approximation amongst which finite element method is the most suitable method due to its sound mathematical foundation. We can possibly consider space-time coupled finite element method based on space-time residual functional or space-time decoupled finite element method in which we generally use Galerkin method with weak form for spatial discretization to convert PDEs in the mathematical model to ODEs in time (nonlinear in this case) which are then integrated in time to obtain the solution. Merits and short comings of both finite element methods can be found in reference [32] and are briefly described in the following.

Space-time coupled finite element method maintains simultaneous dependence of solution on space and time (physics), hence uses space-time finite elements in which time can be viewed as another independent variable. Space-time coupled finite element method in which the space-time integral form is based on space-time residual functional is space-time variationally consistent, hence yields unconditionally stable computational processes during the entire evolution. In the most effective use of this method, a space-time strip or slab for an increment of time is discretized into space-time finite elements. Upon obtaining a converged solution for the first space-time strip, the solution is computed for second space-time strip for the next increment of time using initial conditions from the first space-time strip. This is continued till the final desired time is reached. In this method, we time march only after obtaining converged solution for the current space-time strip or slab. This ensures accuracy of the entire evolution. When the approximation spaces are minimally conforming [32] the space-time residual functional provides accurate measurement of error in the solution. The method potentially has the capability to provide the time accurate computed solutions and is used for obtaining solutions of PDEs in this paper.

In the space-time decoupled finite element method, we consider a discretization in space using Galerkin method with weak form. In this method, local approximation functions are functions of spatial coordinates and the nodal dofs are functions of time. Use of this local approximation for the spatial discretization and the integral form in space based on Galerkin method with weak form yields ODEs in time in the degrees of freedom. These ODEs are then integrated in time to obtain the solutions for the dofs. In this method, error due to decoupling of space and time are difficult to measure. For discretization in space, we can obtain correct or converged solution of ODEs in time, but this may not be the solution of IVP if the discretization and use of p-level in space is not adequate. In general for IVPs in ${ℝ}^3$ , space-time coupled method becomes cumbersome while space-time decoupled methods though more approximate than space-time coupled methods, but remain simple to use. In the present work, we only consider an IVP in ${ℝ}^1$ , hence space-time coupled method is practical, simple to use and provides full control over the accuracy of the computed evolution.

The mathematical model for nonlinear dynamics of TVE solids without memory, derived using CBL of CCM, consists of a system of nonlinear partial differential equations in independent variables $x_i$ and time t. In this paper, this mathematical model is referred to as ‘‘Model A”. By using space-time decoupled

finite element formulation over a discretization ${\bar{\Omega }}_x^T={\displaystyle \underset{e}{\cup }{\bar{\Omega }}_x^e}$ of spatial domain

${\bar{\Omega }}_x$ and by using GM/WF (with some integration by parts in space), the nonlinear PDEs can be reduced to a system of nonlinear ODEs in time. This reduced system of ODEs in time is referred to as ‘‘Model B”. Thus, we can study solutions of nonlinear dynamics problem for TVE solids either by considering the mathematical model A or the mathematical model B. The methods of obtaining solutions of the mathematical model A and B differ significantly and are discussed in the following:

4.2.1. Solutions of PDEs in Space and Time (Model A): Space-Time Coupled Finite Element Method

For simplicity of presenting details of the solution method, let us consider isothermal physics. Then the mathematical model from the CBL of CCM only consists of balance of linear momenta and constitutive theories for ${}_e\sigma {}^{\left[0\right]}$ and ${}_d\sigma {}^{\left[0\right]}$ . Using $\left[{\sigma }^{\left[0\right]}\right]=\left[{}_e{\sigma }^{\left[0\right]}\right]+\left[{}_d{\sigma }^{\left[0\right]}\right]$ , BLM can be written 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]}^T+{\left[{}_d{\sigma }^{\left[0\right]}\right]}^T\right)\right]\{\nabla \}=0.$ (65)

in which

$\left[{}_e{\sigma }^{\left[0\right]}\right]=p\left(\rho \mathrm{,}\theta \right)\left[I\right]\mathrm{.}$ (66)

Here, $p\left(\rho \mathrm{,}\theta \right)$ is the equation of state for compressible matter (thermodynamic pressure).

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

We can substitute (66) and (67) in (65). The resulting equation constitute a system of PDEs in displacements $u_i$ . The mathematical model (65)-(67) consists of three balance of linear momenta equations and six equations from the constitutive theory for $\left[{}_d{\sigma }^{\left[0\right]}\right]$ ), resulting in a total of nine equations in $u_i$ and ${}_d{\sigma }_{ij}^{\left[0\right]}$ , nine dependent variables. Equations (65) and (67) hold for $\forall \left(x_i\mathrm{,}t\right)\in {\bar{\Omega }}_{xt}$ the space time domain of the IVP defined by (65) and (67). We divide space-time

domain into space-time strips ${\bar{\Omega }}_{xt}={\displaystyle \underset{i}{\cup }{\bar{\Omega }}_{xt}^{\left(i\right)}}$ . ${\bar{\Omega }}_{xt}^{\left(i\right)}$ is a typical i^th space-time strip

for $\forall t\in \left[t_i,t_{i+1}\right]$ ; $t_{i+1}-t_i=\Delta t$ . Consider a discretization ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ of n^th space-time strip ${\bar{\Omega }}_{xt}^{\left(n\right)}$ into p-version hierarchical space-time finite elements with higher order global differentiability local approximations.

${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T={\displaystyle \underset{e}{\cup }}{\bar{\Omega }}_{xt}^e\mathrm{.}$ (68)

${\bar{\Omega }}_{xt}^e$ is the space-time element “e”. We calculate converged solution for a space-time strip and time march to obtain entire evolution.

Let $\left\{{\varphi }_h^e\right\}$ be local approximation of $u_i$ , ${}_d\sigma {}_{ij}^{\left[0\right]}$ (equal order, equal degree) over ${\bar{\Omega }}_{xt}^e$ (see reference [23] for details). Then $\left\{{\varphi }_h\right\}={\displaystyle \underset{e}{\cup }}\left\{{\varphi }_h^e\right\}$ is approximation

of $\left\{\varphi \right\}$ over ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ . By substituting $\left\{{\varphi }_h\right\}$ in the nine PDEs (equations (65) and (67)), we obtain nine residual equations $E_i\mathrm{<mo>\; ;\; </mo>}i=\mathrm{1,2,}\cdots \mathrm{,9}$ for ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ in which $E_i$ are the residual functions from the PDEs. Let $\left\{{\delta }^e\right\}$ be the degrees

of freedom in the local approximation $\left\{{\varphi }_h^e\right\}$ , then $\left\{\delta \right\}={\displaystyle \underset{e}{\cup }}\left\{{\delta }^e\right\}$ are the degrees

of freedom for $\left\{{\varphi }_h\right\}$ approximation of $\left\{\varphi \right\}$ over ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ . The space-time residual functional I over ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ is constructed using

$I\left(\left\{\delta \right\}\right)={\displaystyle \underset{1}{\overset{9}{\sum }}}{\left(E_i\mathrm{,}{E}_i\right)}_{{\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T}={\displaystyle \underset{e}{\sum }}\left({\displaystyle \underset{1}{\overset{9}{\sum }}}{\left(E_i^e\mathrm{,}{E}_i^e\right)}_{{\bar{\Omega }}_{xt}^e}={\displaystyle \underset{e}{\sum }}{I}^e\left(\left\{{\delta }^e\right\}\right)\right)$ (69)

in which $E_i^e$ are the residual functions obtained from (65) and (67) when $\left\{{\varphi }_h^e\right\}$ are substituted in them. Functional I(.) is a nonlinear function of $\left\{\delta \right\}$ . Following reference [32] , we proceed as follow. If I(.) is differentiable in $\left\{\delta \right\}$ , then $\delta I\left(\left\{\delta \right\}\right)$ is unique and $\delta I\left(\left\{\delta \right\}\right)=0$ is a necessary condition for an extremum (minimum in this case) of I(.) i.e., the following holds.

$\begin{array}{c}\delta I={\displaystyle \underset{1}{\overset{9}{\sum }}}2{\left(E_i,\delta E_i\right)}_{{\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T}={\displaystyle \underset{e}{\sum }}\left({\displaystyle \underset{1}{\overset{9}{\sum }}}2{\left(E_i^e,\delta E_i^e\right)}_{{\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T}\right)\\ ={\displaystyle \underset{e}{\sum }}\left\{g^e\left(\left\{\delta \right\}\right)\right\}=\left\{g\left(\left\{\delta \right\}\right)\right\}=0\end{array}$ (70)

We can show that the Euler’s equations resulting from (70) are in fact the same as PDEs in (65) and (67). Therefore, a $\left\{\delta \right\}$ yielding extremum of (70) is also a solution of Euler’s equation i.e., PDEs (65) and (67). Since PDEs are nonlinear, $\left\{g\left(\left\{\delta \right\}\right)\right\}$ in (70) is a nonlinear function of $\left\{\delta \right\}$ . We use Newton’s linear method with line search to obtain solution $\left\{\delta \right\}$ from (70).

Let ${\left\{\delta \right\}}_0$ be an assumed (or guess) solution, then

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

Let $\left\{\Delta \delta \right\}$ be incremental change in ${\left\{\delta \right\}}_0$ such that

$\left\{g\left({\left\{\delta \right\}}_0+\left\{\Delta \delta \right\}\right)\right\}=0$ (72)

We expand $\left\{g\left(\mathrm{.}\right)\right\}$ in (72) using a Taylor series expansion about ${\left\{\delta \right\}}_0$ and retain only upto linear terms in $\left\{\Delta \delta \right\}$ :

$\left\{g\left({\left\{\delta \right\}}_0+\left\{\Delta \delta \right\}\right)\right\}=\left\{g\left({\left\{\delta \right\}}_0\right)\right\}+{\left[\frac{\partial \left\{g\right\}}{\partial \left\{\delta \right\}}\right]}_{{\left\{\delta \right\}}_0}\left\{\Delta \delta \right\}=0.$ (73)

From (73), we can solve for $\left\{\Delta \delta \right\}$ :

$\left\{\Delta \delta \right\}=-{\left[\frac{\partial \left\{g\right\}}{\partial \left\{\delta \right\}}\right]}_{{\left\{\delta \right\}}_0}^{-1}\left\{g\left({\left\{\delta \right\}}_0\right)\right\}=-{\left[{\delta }^{2}I\right]}_{{\left\{\delta \right\}}_0}^{-1}\left\{g\left({\left\{\delta \right\}}_0\right)\right\}\mathrm{.}$ (74)

Following reference [6] [32] ${\delta }^{2}I$ is obtained using

${\delta }^{2}I\simeq {\displaystyle \underset{e}{\sum }}\left({\displaystyle \underset{1}{\overset{9}{\sum }}}{\left(\delta E_i^e\mathrm{,}\delta E_i^e\right)}_{{\bar{\Omega }}_{xt}^e}\right)$ (75)

Improved solution is given by (using line search [32] )

$\left\{\delta \right\}={\left\{\delta \right\}}_0+{\alpha }^{\mathrm{<mo>\; *\; </mo>}}\left\{\Delta \delta \right\}\mathrm{<mo>\; ;\; </mo>}0<{\alpha }^{*}<2\text{such}\text{that}I\left(\left\{\delta \right\}\right)\le I\left({\left\{\delta \right\}}_0\right)$ (76)

If

$\left|g_i\left(\left\{\delta \right\}\right)\right|\le \Delta$ (77)

then $\left\{\delta \right\}$ is the converged solution. Otherwise, we set ${\left\{\delta \right\}}_0=\left\{\delta \right\}$ and repeat (74)-(76).

Benefits and highly meritorious feature of STRF finite element method are well documented in [32] . In (77), Δ is a preset tolerance for the computed zero. The entire solution is computed by using converged solution for the current space-time strip and time marching using subsequent space-time strips [32] .

4.2.2. Solutions of Nonlinear ODEs in Time (Model B): Space-Time Decoupled Method

In this section, we consider solutions of system of nonlinear ODEs (1) resulting from space-time decoupled finite element method with GM/WF in space for a spatial discretization. For simplicity, we consider $n=1$ , hence we can write (1) as (using $\left[C\right]$ for $\left[C_{\left[1\right]}\right]$ )

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

following reference [6] , $\left[C\right]$ and $\left[K\right]$ defined by:

$\left[K\right]={\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left(\left[{\left[B^l\right]}^T+{\left[B^{nl}\right]}^T\right]\left[D\right]\left[\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right]\text{d}\Omega \right)$ (79)

$\left[C\right]={\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left(\left[{\left[B^l\right]}^T+{\left[B^{nl}\right]}^T\right]\left[D_1\right]\left[\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right]\text{d}\Omega \right)$ (80)

$\left[M\right]$ is the usual consistent mass matrix (defined in ref [6] with constant coefficients).

The coefficients of $\left[K\right]$ and $\left[C\right]$ are up to quadratic functions of the gradients of displacements $u_i$ .

We present details of Newmark linear acceleration method for obtaining solution of (78) in the following. The nonlinear ODEs in time (78) can be integrated numerically using a variety of methods [32] . In solid and structural mechanics, the decoupling of space and time is done using Galerkin method with weak form in space. The ODEs in time generally contain mass, stiffness and damping matrices as in (78). Wilson’s θ-method and Newmark’s constant average acceleration method and linear acceleration methods are well established time integration methods for (78). For good accuracy of computed solutions (no amplitude decay and base elongation), the integration time step generally needs to be much

smaller than $\frac{1}{20}\left(\frac{2\pi }{\omega }\right)$ (the time period is divided in twenty increments). For

such time steps, accuracy of Wilson’s θ-method and Newmark methods are comparable but Newmark method is more straight forward, as in this method, we consider equilibrium at $t+\Delta t$ and not at $t+\theta \Delta t$ as in Wilson’s θ-method, hence used in the present work. Secondly, since the ODEs ((78)) are nonlinear, for every increment of time we need to perform iterations to converge to the correct solution. Newton’s linear method has quadratic convergence, hence is meritorious to use. Thus, the choice of Newmark linear acceleration method (more accurate than constant average acceleration [32] ) with Newton’s linear method.

In the following, we refer to $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ as displacements, velocities and accelerations (since $\left\{\delta \right\}$ are degrees of freedom associated with displacements and $\left\{\stackrel{\dot{}}{\delta }\right\},\left\{\ddot{\delta }\right\}$ are first and second time derivatives of $\left\{\delta \right\}$ ). Details of Newmark linear acceleration methods are given in ref [32] for a system of linear second order ODEs. We present the details of Newmark linear acceleration method for a system of nonlinear second order ODEs (78) in the following.

Let ${\left\{\delta \right\}}_t$ , ${\left\{\stackrel{\dot{}}{\delta }\right\}}_t$ and ${\left\{\ddot{\delta }\right\}}_t$ be known solutions at time t and ${\left\{\delta \right\}}_{t+\Delta t}$ (or $\left\{\delta \left(t+\Delta t\right)\right\}$ ), ${\left\{\stackrel{\dot{}}{\delta }\right\}}_{t+\Delta t}$ (or ${\left\{\stackrel{\dot{}}{\delta }\right\}}_{t+\Delta t}$ ) and ${\left\{\ddot{\delta }\right\}}_{t+\Delta t}$ (or $\left\{\ddot{\delta }\left(t+\Delta t\right)\right\}$ ) be unknown solution at time $t+\Delta t$ . Let the acceleration be linear between t to $t+\Delta t$ . If $\tau$ is time measured from time t, then based on linear acceleration assumption in interval $\left[t\mathrm{,}t+\Delta t\right]$ , we can write

${\left\{\ddot{\delta }\right\}}_{t+\tau }={\left\{\ddot{\delta }\right\}}_t+\frac{\tau }{\Delta t}\left({\left\{\ddot{\delta }\right\}}_{t+\Delta t}-{\left\{\ddot{\delta }\right\}}_t\right)\mathrm{.}$ (81)

Integrating (81) with respect to (w.r.t) $\tau$ and using $\tau =0$ to determine constant of integration, we can obtain the following for velocity at $t+\tau$ .

${\left\{\stackrel{\dot{}}{\delta }\right\}}_{t+\tau }={\left\{\stackrel{\dot{}}{\delta }\right\}}_t+\tau {\left\{\ddot{\delta }\right\}}_t+\frac{{\tau }^2}{2\Delta t}\left({\left\{\ddot{\delta }\right\}}_{t+\Delta t}-{\left\{\ddot{\delta }\right\}}_t\right)\mathrm{.}$ (82)

Integrating (82) w.r.t $\tau$ and using $\tau =0$ to determine constant of integration, we can obtain the following for displacement ${\left\{\delta \right\}}_{t+\tau }$ i.e., we consider the following.

${\left\{\delta \right\}}_{t+\tau }={\left\{\delta \right\}}_t+\tau {\left\{\stackrel{\dot{}}{\delta }\right\}}_t+\frac{{\tau }^2}{2}{\left\{\ddot{\delta }\right\}}_t+\frac{{\tau }^3}{6\Delta t}\left({\left\{\ddot{\delta }\right\}}_{t+\Delta t}-{\left\{\ddot{\delta }\right\}}_t\right)$ (83)

In Newmark linear acceleration method, we satisfy (78) at $t+\Delta t$ i.e.,

$\left[M\right]{\left\{\ddot{\delta }\right\}}_{t+\Delta t}+\left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]{\left\{\stackrel{\dot{}}{\delta }\right\}}_{t+\Delta t}+\left[K\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]{\left\{\delta \right\}}_{t+\Delta t}=\left\{P\right\}+\left\{F\right\}\mathrm{.}$ (84)

In (84), the coefficients of $\left[K\right]$ and $\left[C\right]$ are up to quadratic functions of ${\left\{\delta \right\}}_{t+\Delta t}$ . In (82) and (83), we substitute $\tau =\Delta t$ to obtain

${\left\{\stackrel{\dot{}}{\delta }\right\}}_{t+\Delta t}={\left\{\stackrel{\dot{}}{\delta }\right\}}_t+\frac{\Delta t}{2}{\left\{\ddot{\delta }\right\}}_t+\frac{\Delta t}{2}{\left\{\ddot{\delta }\right\}}_{t+\Delta t}$ (85)

${\left\{\delta \right\}}_{t+\Delta t}={\left\{\delta \right\}}_t+\Delta t{\left\{\stackrel{\dot{}}{\delta }\right\}}_t+\frac{\Delta t^2}{3}{\left\{\ddot{\delta }\right\}}_t+\frac{\Delta t^2}{6}{\left\{\ddot{\delta }\right\}}_{t+\Delta t}\mathrm{.}$ (86)

From (86), we determine ${\left\{\ddot{\delta }\right\}}_{t+\Delta t}$ and then substitute it into (85) to determine ${\left\{\stackrel{\dot{}}{\delta }\right\}}_{t+\Delta t}$ .

${\left\{\ddot{\delta }\right\}}_{t+\Delta t}=\frac{6}{{\left(\Delta t\right)}^2}{\left\{\stackrel{\dot{}}{\delta }\right\}}_{t+\Delta t}-{\left\{\delta \right\}}_t-\frac{6}{\Delta t}{\left\{\stackrel{\dot{}}{\delta }\right\}}_t-2{\left\{\ddot{\delta }\right\}}_t$ (87)

and

${\left\{\stackrel{\dot{}}{\delta }\right\}}_{t+\Delta t}=\frac{3}{\Delta t}\left({\left\{\delta \right\}}_{t+\Delta t}-{\left\{\delta \right\}}_t\right)-2{\left\{\stackrel{\dot{}}{\delta }\right\}}_t-\frac{\Delta t}{2}{\left\{\ddot{\delta }\right\}}_t$ (88)

Substituting (87) and (88) in (84) and regrouping terms gives

$\begin{array}{l}\left[\frac{6}{{\left(\Delta t\right)}^2}\left[M\right]+\frac{3}{\Delta t}\left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]+\left[K\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]\right]{\left\{\delta \right\}}_{t+\Delta t}\\ ={\left\{F\right\}}_{t+\Delta t}+{\left\{P\right\}}_{t+\Delta t}+\left[\frac{6}{{\left(\Delta t\right)}^2}\left[M\right]+\frac{3}{\Delta t}\left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]\right]{\left\{\delta \right\}}_t\\ +\left[\frac{6}{\Delta t}\left[M\right]+2\left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]\right]{\left\{\stackrel{\dot{}}{\delta }\right\}}_t+\left[2\left[M\right]+\frac{\Delta t}{2}\left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]\right]{\left\{\ddot{\delta }\right\}}_t\mathrm{.}\end{array}$ (89)

Equation (89) is a system of nonlinear algebraic equations in ${\left\{\delta \right\}}_{t+\Delta t}$ . We use Newton’s linear method to obtain solution ${\left\{\delta \right\}}_{t+\Delta t}$ of (89). Details are given in the following.

We rewrite (89) as follows:

$\begin{array}{c}\left\{g\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right\}=\left[\frac{6}{\Delta t^2}\left[M\right]+\frac{3}{\Delta t}\left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]+\left[K\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]\right]{\left\{\delta \right\}}_{t+\Delta t}\\ -\left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]\left(\frac{3}{\Delta t}{\left\{\delta \right\}}_t+2{\left\{\stackrel{\dot{}}{\delta }\right\}}_t+\frac{\Delta t}{2}{\left\{\ddot{\delta }\right\}}_t\right)\\ -\left(\frac{6}{{\left(\Delta t\right)}^2}\left[M\right]{\left\{\delta \right\}}_t+\frac{6}{\Delta t}\left[M\right]{\left\{\stackrel{\dot{}}{\delta }\right\}}_t+2\left[M\right]{\left\{\ddot{\delta }\right\}}_t\right)\\ -{\left\{F\right\}}_{t+\Delta t}-{\left\{P\right\}}_{t+\Delta t}\\ =0.\end{array}$ (90)

Let ${\left\{\delta \right\}}_{t+\Delta t}^0$ be initial or known guess of ${\left\{\delta \right\}}_{t+\Delta t}$ at $t=0$ (commencement of evolution), then

$\left\{g\left({\left\{\delta \right\}}_{t+\Delta t}^0\right)\right\}\ne 0.$ (91)

Let ${\left\{\Delta \delta \right\}}_{t+\Delta t}$ be incremental change in ${\left\{\delta \right\}}_{t+\Delta t}^0$ such that

$\left\{g\left({\left\{\delta \right\}}_{t+\Delta t}^0+{\left\{\Delta \delta \right\}}_{t+\Delta t}\right)\right\}=0.$ (92)

We expand $\left\{g\left(\mathrm{.}\right)\right\}$ in (92) using Taylor series expansion about ${\left\{\delta \right\}}_{t+\Delta t}^0$ and retain only up to linear terms in ${\left\{\Delta \delta \right\}}_{t+\Delta t}$

$\left\{g\left({\left\{\delta \right\}}_{t+\Delta t}^0+{\left\{\Delta \delta \right\}}_{t+\Delta t}\right)\right\}=\left\{g\left({\left\{\delta \right\}}_{t+\Delta t}^0\right)\right\}+{\left[\frac{\partial \left\{g\right\}}{\partial {\left\{\delta \right\}}_{t+\Delta t}}\right]}_{{\left\{\delta \right\}}_{t+\Delta t}^0}{\left\{\Delta \delta \right\}}_{t+\Delta t}=0.$ (93)

From (93), we can solve for ${\left\{\Delta \delta \right\}}_{t+\Delta t}$

${\left\{\Delta \delta \right\}}_{t+\Delta t}=-{\left[\frac{\partial \left\{g\right\}}{\partial {\left\{\delta \right\}}_{t+\Delta t}}\right]}_{{\left\{\delta \right\}}_{t+\Delta t}^0}^{-1}\left\{g\left({\left\{\delta \right\}}_{t+\Delta t}^0\right)\right\}=-{\left[\delta \left\{g\right\}\right]}_{{\left\{\delta \right\}}_{t+\Delta t}^0}^{-1}\left\{g\left({\left\{\delta \right\}}_{t+\Delta t}^0\right)\right\}.$ (94)

We refer to the coefficient matrix $\delta \left\{g\right\}$ as tangent matrix ${\left[H\right]}_T$ , equation (105) that controls determination of incremental solution. Improved solution ${\left\{\delta \right\}}_{t+\Delta t}$ is obtained using

${\left\{\delta \right\}}_{t+\Delta t}={\left\{\delta \right\}}_{t+\Delta t}^0+{\left\{\Delta \delta \right\}}_{t+\Delta t}\mathrm{.}$ (95)

We check for convergence of the Newton’s linear method.

If

$\left|g_i\left(\left\{\delta \right\}\right)\right|\le \Delta$ (96)

then ${\left\{\delta \right\}}_{t+\Delta t}$ is the converged solution. Otherwise, we set ${\left\{\delta \right\}}_{t+\Delta t}^0={\left\{\delta \right\}}_{t+\Delta t}$ and repeat (94)-(96) till (96) is satisfied. Δ in (96) is a preset tolerance for the computed zero.

We still need $\delta \left\{g\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right\}$ in (94) to compute ${\left\{\Delta \delta \right\}}_{t+\Delta t}$ . This is obtained using (90). First, we note

$\delta \left(\left[M\right]{\left\{\delta \right\}}_{t+\Delta t}\right)=\left[M\right]$ (97)

$\delta \left(\left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]{\left\{\delta \right\}}_{t+\Delta t}\right)=\delta \left({\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{\left[B\right]}^T\left\{{\sigma }^1\right\}\text{d}\Omega \right)$ (98)

$\delta \left(\left[K\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]{\left\{\delta \right\}}_{t+\Delta t}\right)=\delta \left({\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{\left[B\right]}^T\left\{{}_d{\sigma }^{\left[0\right]}\right\}\text{d}\Omega \right)$ (99)

in which

$\left\{{}_d{\sigma }^1\right\}=\left[D_1\right]\left(\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right){\left\{\delta \right\}}_{t+\Delta t}$ (100)

and $\left\{{}_d{\sigma }^{\left[0\right]}\right\}=\left[D\right]\left(\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right){\left\{\delta \right\}}_{t+\Delta t}\mathrm{.}$ (101)

Following [24] , we can derive the following

$\begin{array}{c}\delta \left({\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{\left[B\right]}^T\left\{{\sigma }^1\right\}\text{d}\Omega \right)={\displaystyle \underset{e}{\sum }}\left({\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left({\left[B\right]}^T\left[D_1\right]\left[B\right]+{\left[G\right]}^T\left[S^1\right]\left[G\right]\right)\text{d}\Omega \right)\\ ={\displaystyle \underset{e}{\sum }}\left[C_T^e\right]=\left[C_T\right]\end{array}$ (102)

$\begin{array}{c}\delta \left({\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}{\left[B\right]}^T\left\{{\sigma }^{\left[0\right]}\right\}\text{d}\Omega \right)={\displaystyle \underset{e}{\sum }}\left({\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left({\left[B\right]}^T\left[D\right]\left[B\right]+{\left[G\right]}^T\left[S\right]\left[G\right]\text{d}\Omega \right)\right)\\ ={\displaystyle \underset{e}{\sum }}\left[K_T^e\right]=\left[K_T\right]\end{array}$ (103)

$\left[C_T\right]$ and $\left[K_T\right]$ are tangent damping and stiffness matrices for spatial discretization.

$\left[S\right]$ is given by

$\left[S\right]=\left[\begin{array}{ccc}{\sigma }_{11}I& {\sigma }_{12}I& {\sigma }_{13}I\\ & {\sigma }_{22}I& {\sigma }_{23}I\\ \begin{array}{c}\text{Symm}\end{array}& & {\sigma }_{33}I\end{array}\right]$ (104)

$\left[S^1\right]$ is obtained by replacing ${}_d{\sigma }_{ij}^{\left[0\right]}$ in (104) by ${\sigma }_{ij}^1$ (defined in (100)).

By substituting (102) and (103) into (98) and (99), and then using these in $\delta \left\{g\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right\}$ , we obtain the following expression in which $\left\{\tilde{\delta }\right\}=\frac{3}{\Delta t}{\left\{\delta \right\}}_t+2{\left\{\stackrel{\dot{}}{\delta }\right\}}_t+\frac{\Delta t}{2}{\left\{\ddot{\delta }\right\}}_t$ :

$\delta \left\{g{\left\{\delta \right\}}_{t+\Delta t}\right\}=\frac{6}{\Delta {\tau }^2}\left[M\right]+\frac{3}{\Delta t}\left[C_T\right]+\left[K_T\right]-\delta \left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\left\{\tilde{\delta }\right\}\right]=\left[H_T\right]\mathrm{.}$ (105)

The matrix $\left[H_T\right]$ is called tangent matrix.

We still need $\delta \left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]\left\{\tilde{\delta }\right\}$ .

$\delta \left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\left\{\tilde{\delta }\right\}\right]={\displaystyle \underset{e}{\sum }}\delta {\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}\left(\left[{\left[B^l\right]}^T+{\left[B^{nl}\right]}^T\right]\left[D_1\right]\left[\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right]\left\{{\tilde{\delta }}^e\right\}\right)\text{d}\Omega$ (106)

$\begin{array}{l}={\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}_x^e}{\int }}(\left[\delta {\left[B^{nl}\right]}^T\right]\left[D_1\right]\left[\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right]\left\{{\tilde{\delta }}^e\right\}\\ +\left[{\left[B^l\right]}^T+{\left[B^{nl}\right]}^T\right]\left[D_1\right]\delta \left(\frac{1}{2}\left[B^{nl}\right]\left\{{\tilde{\delta }}^e\right\}\right))\text{d}\Omega .\end{array}$ (107)

Let $\left[D_1\right]\left[\left[B^l\right]+\frac{1}{2}\left[B^{nl}\right]\right]\left\{\tilde{\delta }\right\}=\left\{\sigma \right\}$

Then $\delta {\left[B^{nl}\right]}^T\left\{\sigma \right\}=\delta {\left(\left[A_g\right]\left[G\right]\right)}^T\left\{\sigma \right\}={\left[G\right]}^T\delta {\left[A_g\right]}^T\left\{\sigma \right\}={\left[G\right]}^T\left[S\right]\left[G\right]$ and

$\begin{array}{l}\left[{\left[B^l\right]}^T+{\left[B^{nl}\right]}^T\right]\left[D_1\right]\delta \left(\frac{1}{2}\left[B^{nl}\right]\left\{{\tilde{\delta }}^e\right\}\right)\\ =\left[{\left[B^l\right]}^T+{\left[B^{nl}\right]}^T\right]\left[D_1\right]\delta \left(\left[A_g\right]\left[G\right]\left\{{\tilde{\delta }}^e\right\}\right)\\ =\left[{\left[B^l\right]}^T+{\left[B^{nl}\right]}^T\right]\left[D_1\right]\delta \left(\left[A_g\right]\left\{{\tilde{t}}^e\right\}\right)\end{array}$ (108)

where $\left\{{\tilde{t}}^e\right\}=\left[G\right]\left\{{\tilde{\delta }}^e\right\}$

$\begin{array}{c}\delta \left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]\left\{\tilde{\delta }\right\}={\displaystyle \underset{e}{\sum }}{\displaystyle \underset{{\bar{\Omega }}_x}{\int }}\left({\left[G\right]}^{T}S\left[G\right]+\left({\left[B^l\right]}^T+{\left[B^{nl}\right]}^T\right)\left[D_1\right]\delta \left(\left[A_g\right]\left\{\tilde{t}\right\}\right)\right)\text{d}\Omega \\ ={\displaystyle \underset{e}{\sum }}\left[{\widetilde{\stackrel{˜}{C}}}^e\right]=\left[\widetilde{\stackrel{˜}{C}}\right].\end{array}$ (109)

Substituting (109) in (105), we have $\delta \left\{g\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right\}$ , hence ${\left\{\Delta \delta \right\}}_{t+\Delta t}$ can be calculated using (94).

Remarks

Since, now we know the dependence of incremental solution on ‘tangent matrix’ $\left[H_T\right]$ , we can discuss various physics in $\left[H_T\right]$ that influences nonlinear dynamic response and dynamic bifurcation.

(1) We note that for an increment of time, incremental solution calculation (in Newton’s linear method) using (94) requires $\delta \left\{g\right\}$ defined in (105). $\delta \left\{g\right\}$ in (105) is the coefficient matrix that must be inverted to obtain incremental solution. Thus, the time response of a nonlinear dynamic system and dynamic bifurcation depends upon $\left[H_T\right]$ .

(2) Now we can conclusively determine the influence of various physics through $\left[H_T\right]$ on nonlinear dynamic response and the dynamic bifurcation phenomenon.

(3) Nonlinear dynamic response depends upon symmetric mass matrix $\left[M\right]$ , tangent damping matrix $\left[C_T\right]$ and tangent stiffness matrix $\left[K_T\right]$ (as well as the last term in (105) containing variation of $\left[C\right]$ ). We recall that

$\left[K_T\right]=\left[K_1\right]+\left[{\underset{˜}{K}}_2\right]+\left[{\underset{˜}{K}}_3\right]+\left[K_{\sigma }\right]$ (110)

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 dofs $\left\{\delta \right\}$ . $\left[K_{\sigma }\right]$ is the stiffness due to stress field. We note that the presence of finite deformation, finite strain physics is essential for dynamic bifurcation to exist, as in the absence of these $\left[{\underset{˜}{K}}_2\right]$ , $\left[{\underset{˜}{K}}_3\right]$ and $\left[K_{\sigma }\right]$ are all zero. Negative $\left[K_{\sigma }\right]$ (due to compressive field) and positive $\left[K_{\sigma }\right]$ (due to tensile stress field) are likely to result in different dynamic response, hence different dynamic bifurcations.

(4) Tangent damping matrix $\left[C_T\right]$ has identical structure as $\left[K_T\right]$ only the material coefficients are different.

$\left[C_T\right]=\left[C_1\right]+\left[{\underset{˜}{C}}_2\right]+\left[{\underset{˜}{C}}_3\right]+\left[C_{\sigma }\right]$ (111)

Matrix $\left[C_{\sigma }\right]$ is similar to $\left[K_{\sigma }\right]$ . It uses $\left[D_i\right]$ as opposed to $\left[D\right]$ . Thus, in calculating $\left[D_i\right]$ , $2\mu$ and $\lambda$ in $\left[D\right]$ are replaced by $2{\mu }_i$ , ${\lambda }_i$ corresponding to strain rate $\left[{\epsilon }_{\left[i\right]}\right]$ . Thus, constant, linear and quadratic damping as well as $\left[C_{\sigma }\right]$ influence nonlinear dynamics. In general, dissipation provides resistance to motion, reduced amplitude of motion, thus more stability. Hence, increasing dissipation presence is likely to inhibit existence of dynamic bifurcation.

(5) Presence of $\left[M\right]$ , hence translation inertial physics results in lowering of the stiffness of the system. Thus, the increasing influence of this physics increases likelihood of the existence of dynamic bifurcation.

(6) Static bifurcation in BVPs only depends upon $\left[K_{\sigma }\right]$ (precisely in the same form as defined here). Thus, we see that static bifurcation and dynamic bifurcation are totally two different physical phenomena. Static bifurcation is not a necessary condition for dynamic bifurcation (also shown later in the model problem studies) as there is obviously no connection between the two.

5. Model Problems and Their Solutions

We consider two model problems and present their solutions obtained using space-time coupled method based on space-time residual functional [32] . The first model problem is a 1D phenomenologically constructed mathematical model for finite deformation, finite strain nonlinear dynamics using a mass, a nonlinear spring and a nonlinear dashpot (64), which is reproduced here for convenience.

Model problem I:

$m\ddot{u}+\left(c_1+c_{2}u+c_3{u}^2\right)\stackrel{\dot{}}{u}+\left(k_1+k_{2}u+k_3{u}^2\right)u=f_0\sin \left(\omega t\right)$ (112)

In this mathematical model, both stiffness and damping are up to quadratic functions of the degrees of freedom u. This mathematical model is lumped in space; hence does not contain space coordinates. This is obviously nonphysical for a continuous system.

Model problem II:

A mathematical model that appears somewhat parallel to (112), but based on CBL of CCM, is considered in this model problem.

We consider a TVE axial rod 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). If all points in each cross section of the rod deform in $x_1$ direction by the same amount, we can idealize the rod as a line, as shown in Figure 1(b). For simplicity, we consider isothermal physics. Thus, balance of linear momenta in $x_1$ and the constitutive theory for deviatoric second Piola-Kirchhoff stress ${}_d{\sigma }_{11}^{\left[0\right]}$ are the only two equations that constitute the mathematical model, assuming equilibrium stress ${}_e{\sigma }^{\left[0\right]}$ to be zero. We can write the following for BLM and the constitutive theory (in Lagrangian description using Equations (19)-(29)):

${\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_i}{\partial x_1}\right]{}_d\sigma {}_{11}^{\left[0\right]}\right)=0\text{BLM}\forall x\mathrm{,}t\in {\Omega }_{xt}$ (113)

${}_d{\sigma }_{11}^{\left[0\right]}=\left(2\mu +\lambda \right){\left({\epsilon }_{\left[0\right]}\right)}_{11}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{c}_i{\left({\epsilon }_{\left[i\right]}\right)}_{11}$ (114)

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

${\left({\epsilon }_{\left[i\right]}\right)}_{11}=\frac{\partial }{\partial t}{\left({\epsilon }_{\left[i-1\right]}\right)}_{11};i=1,2,\cdots ,n.$ (116)

We substitute (115) and (116) into (114), then ${}_d{\sigma }_{11}^{\left[0\right]}$ from (114) into (113) to obtain a single PDE in $u_1$ (letting $2\mu +\lambda =E$ )

$\begin{array}{c}A{u}_1-f={\rho }_0\frac{{\partial }^2\text{}{u}_1}{\partial t{\text{}}^2}-{\rho }_0{F}_1^b-\frac{\partial }{\partial x_1}(\left(1+\frac{\partial u_1}{\partial x_1}\right)\left(E\left(\frac{\partial u_1}{\partial x_1}+\frac{1}{2}{\left(\frac{\partial u_1}{\partial x_1}\right)}^2\right)\right)\\ -{\displaystyle \underset{i\mathrm{<mo>\; =\; </mo>1}}{\overset{n}{\sum }}}{c}_i\frac{{\partial }^i}{\partial t^i}\left(\left(\frac{\partial u_1}{\partial x_1}\right)+\frac{1}{2}{\left(\frac{\partial u_1}{\partial x_1}\right)}^2\right))\\ =0.\end{array}$ (117)

Figure 1. 1D axial rod, spatial domain and discretization (Model Problem II).

5.1. Space-Time Decoupled Finite Element Method (for Model Problem II)

We construct a space-time decoupled finite element formulation (117). Let ${\bar{\Omega }}_x^T={\displaystyle \underset{e}{\cup }}{\bar{\Omega }}_x^e$ be spatial discretization of ${\bar{\Omega }}_x$ (Figure 1(c)).

Let ${\left(u_1\left(x\mathrm{,}t\right)\right)}_h^e$ be local approximation of $u_1$ over ${\bar{\Omega }}_x^e$ . Then

${\left(u_1\left(x\mathrm{,}t\right)\right)}_h={\displaystyle \underset{e}{\cup }}{\left(u_1\left(x\mathrm{,}t\right)\right)}_h^e$ is the approximation of $u_1$ over ${\bar{\Omega }}_x^T$ . The integral

form of (117) over ${\bar{\Omega }}_x^T$ based on fundamental lemma of the calculus of variations using test function v is given by:

${\left(A{\left(u_1\right)}_h-f,v\right)}_{{\bar{\Omega }}_x^T}=0$ (118)

or

${\left(A{\left(u_1\right)}_h-f,v\right)}_{{\bar{\Omega }}_x^T}={\displaystyle \underset{e}{\sum }}{\left(A{\left(u_1\right)}_h^e-f,v\right)}_{{\bar{\Omega }}_x^e}=0.$ (119)

We consider scalar product ${\left(A{\left(u_1\right)}_h^e-f\mathrm{,}v\right)}_{{\bar{\Omega }}_x^e}$ , where

${\left(u_1\right)}_h^e=\left[N\left(x\right)\right]\left\{{\delta }^e\left(t\right)\right\}\mathrm{.}$ (120)

$\left[N\left(x\right)\right]$ is approximation functions matrix and $\left\{{\delta }^e\left(t\right)\right\}$ are time dependent nodal degrees of freedom for an element e. The total degrees of freedom $\left\{\delta \right\}$ for ${\bar{\Omega }}_x^T$ are given by:

$\left\{\delta \right\}={\displaystyle \underset{e}{\cup }}\left\{{\delta }^e\right\}\mathrm{.}$ (121)

Using (121) in the integral form for an element “e” with domain ${\bar{\Omega }}_x^e$ , and following reference [23] , i.e., using GM/WF in space for ${\bar{\Omega }}_x^e$ followed by assembly of element equations based on (119), we obtain the following (similar to (1))

$\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\{P\right\}+\left\{\sigma \right\}$ (122)

where [24]

$\left[C_i\right]=\left[C_i^1\right]+\left[C_i^2\right]+\left[C_i^3\right]$ (123)

$\left[K\right]=\left[K^1\right]+\left[K^2\right]+\left[K^3\right]\mathrm{.}$ (124)

$\left[C_i^1\right]$ and $\left[K^1\right]$ have constant coefficients where as coefficients of $\left[C_i^2\right]\mathrm{,}\left[K^2\right]$ and $\left[C_i^3\right]\mathrm{,}\left[K^3\right]$ are linear and quadratic functions of displacement gradients

( $\frac{\partial u_1}{\partial x_1}$ in this case). System of second order ODEs in (122) based on CBL of CCM

(with GM/WF for spatial discretization) are a set of coupled nonlinear ODEs in time. In (112), $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3\mathrm{,}{k}_1\mathrm{,}{k}_2\mathrm{,}{k}_3$ are purely hypothetical or phenomenological, whereas in (122), $2{\eta }_i\mathrm{,}{\kappa }_i$ used in $\left[C_i\right]$ and $2\eta ,\lambda$ used in $\left[K\right]$ are actual material coefficients of TVE solid matter. At this stage, it might appear that we can study the trends in the physics of (122) by using phenomenological model (112).

We make some remarks:

(1) For a continuous system, mathematical models are always PDEs in space and time. It is only after decoupling space and time using GM/WF in space for a spatial discretization ${\bar{\Omega }}_x^T$ of ${\bar{\Omega }}_x$ that we obtain a system of nonlinear ODEs in that contain degrees of freedom $\left\{\delta \right\}$ and $\left\{\stackrel{\dot{}}{\delta }\right\}\mathrm{,}\left\{\ddot{\delta }\right\}$ corresponding to the spatial discretization ${\bar{\Omega }}_x^T$ . Thus, in the system of ODEs in time (122) information at spatial locations are preserved. In other words, system of ODEs in time (122) is a reflection of preservation of spatial information contained in original PDEs in (122) after decoupling space and time. This is always the case for continuous matter. For continuous system with a spatial discretization, space-time decoupling will never result in a single ODE such as (112).

(2) Based on the discussion in (1), it is conclusive that (112) cannot describe the same physics as (122).

(3) We note that, when integrating (112) using Newmark’s linear acceleration method we take the inverse of the tangent matrix $\left[H_T\right]$ defined by (105) i.e., inverse of the following matrix for calculating incremental change (or improvement) in the assumed solution:

$\left[H_T\right]=\frac{6}{\Delta {\tau }^2}\left[M\right]+\frac{3}{\Delta t}\left[C_T\right]+\left[K_T\right]-\delta \left[C\left({\left\{\delta \right\}}_{t+\Delta t}\right)\right]\left\{\tilde{\delta }\right\}$ (125)

Matrices $\left[C_T\right]$ and $\left[K_T\right]$ are derived considering $n=1$ , which corresponds to $\left[C_1\right]$ and $\left[K\right]$ in (122). $\left[C_T\right]$ and $\left[K_T\right]$ consist of

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

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

$\left[C^1\right]\mathrm{,}\left[K^1\right]\mathrm{,}\left[{\underset{˜}{C}}^i\right]\mathrm{,}\left[{\underset{˜}{K}}^i\right]\mathrm{<mo>\; ;\; </mo>}i=\mathrm{2,3}$ and $\left[C_{\sigma }\right]\mathrm{,}\left[K_{\sigma }\right]$ are all symmetric. $\left[C^1\right]$ and $\left[K^1\right]$ have constant coefficients but coefficients of $\left[{\underset{˜}{C}}^1\right]\mathrm{,}\left[{\underset{˜}{K}}^2\right]$ and $\left[{\underset{˜}{C}}^3\right]\mathrm{,}\left[{\underset{˜}{K}}^3\right]$ are linear and quadratic in displacement gradients. $\left[K_{\sigma }\right]$ describes the influence of stress field on the stiffness. $\left[C_{\sigma }\right]$ likewise has similar structure as that of $\left[K_{\sigma }\right]$ (as shown earlier).

(4) Thus we note that (123), (124) and (125), (126) are equivalent in terms of influencing stiffness of the volume of matter. $\left[C_{\sigma }\right]$ and $\left[K_{\sigma }\right]$ are implicitly contained in (123) and (124) but can be explicitly observed in (126) and (127). Thus, in the following, we can view $\left[K\right]$ in (122) is due to constant, linear and quadratic functions of displacement gradients as well as due to $\left[K_{\sigma }\right]$ , the influence of stress field on stiffness. Similar view holds for damping matrices $\left[C_i\right]$ . Compressive stress field results in negative $\left[K_{\sigma }\right]$ , causing a reduction in stiffness. On the other hand, tensile stress field resulting in positive $\left[K_{\sigma }\right]$ will cause an increase in stiffness. Explicit appearance of $\left[K_{\sigma }\right]$ provides additional information regarding influence of stress on total stiffness. This can only be observed and extracted explicitly in this mathematical model based on CCM and the approach of obtaining the solution presented here.

(5) In model problem I, explicit realization of $\left[K_{\sigma }\right]$ is not possible. However, the presence of $u^2$ in stiffness and damping as well as extraction of $\left[K_{\sigma }\right]$ and $\left[C_{\sigma }\right]$ in model problem II from the linear and quadratic terms in displacement gradients in the stiffness and damping matrices, suggest that positive and negative values of $\left[K_{\sigma }\right]$ and $\left[C_{\sigma }\right]$ must be associated with positive or negative values of $k_3$ and $c_3$ . At this stage, it is intuitively plausible to assume that a negative $k_3$ will result in reduction of stiffness (similar to negative $\left[K_{\sigma }\right]$ in model problem II) and positive $k_3$ will result in increase of stiffness (similar to positive $\left[K_{\sigma }\right]$ ) in model problem II.

(6) It is worth pointing out that in case of model problem II, the applied excitation causes compression, hence will only result in negative $\left[K_{\sigma }\right]$ . In this model problem, there is no mechanism of positive $\left[K_{\sigma }\right]$ as it would require a tensile stress field in this rod. However in case of model problem I, we can study dynamic bifurcation physics for both tension and compression simply by choosing positive or negative $k_3$ as this model is not sensitive to the physics as model problem II is.

5.2. Space-Time Coupled Finite Element Method (Model Problem II)

The mathematical model (113)-(116) is a system of first order PDEs that facilitate specifications of BCs and ICs. This form is also more convenient for using space-time coupled finite element method on a space-time strip in conjunction with time marching to obtain entire evolution for $t\in \left[\mathrm{0,}\tau \right]$ . Considering zero body forces and $n=1$ (for simplicity):

${\rho }_0\frac{\partial v_1}{\partial t}-\left(\frac{{\partial }^2\text{}{u}_1}{\partial x_1^2}\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)\right)-\frac{\partial }{\partial x_1}\left({}_d{\sigma }_{11}^{\left[0\right]}\right)=0$ (128)

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

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

Equations (128)-(130) are three nonlinear PDEs in $u_1\mathrm{,}{}_d\sigma {}_{11}^{\left[0\right]}\mathrm{,}{v}_1,x_1$ and t. Space-time coupled finite element method is used to obtain their solution using space-time integral form based on space-time residual functions (details are presented in Section 4.2.1).

Remarks

(1) Space-time residual functional I for ${\left({\bar{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ , the discretization of n^{th $n^{th}$} space-time strip provides a measure of accuracy of the solution (or the error in the solution). When I tends to zero, the computed solution approaches theoretical solution. For each space time step, I of the order of $O\left({10}^{-8}\right)$ or lower is achieved to ensure accurate evolution for each space-time strip.

(2) This feature (described in (1)) is missing in space-time decoupled methods. Accuracy of the solutions of ODEs in time can be measured or quantified but the error due to discretization in space and those due to decoupling of space and time cannot be quantified. As a result, it is difficult to assess the accuracy of the evolution of the PDEs in the space-time decoupled method.

(3) In the present work, we only consider space-time coupled finite element method based on space-time residual functional (STRF) due to high accuracy and unconditional stability [32] of the method.

6. Static Bifurcation

The purpose of these studies is to demonstrate the influence of matrix $\left[K_{\sigma }\right]$ due to stress field on the total stiffness and the resulting deformation. We study conditions under which bifurcation physics can exist. We consider model problem I as well as model problem II. BVP associated with model problem I (using (64)) is given by:

$\left(k_1+k_{2}u+k_3{u}^2\right)u=F$ (131)

in which F is the applied force.

Model problem II for BVP can be written (in the absence of body forces) as

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

${}_d{\sigma }_{11}^{\left[0\right]}=\left(2\mu +\lambda \right){\left({\epsilon }_{\left[0\right]}\right)}_{11}=E{\left({\epsilon }_{\left[0\right]}\right)}_{11}$ (133)

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

Numerical Solutions of Model Problem I and Model Problem II

Model problem I:

First, we consider model problem I, a cubic nonlinear algebraic equation. We vary the force F and calculate corresponding u for each value of F using Newton’s linear method. We choose $0\le F\le 2.5$ for all studies. In the first study, we choose $k_1=k_2=k_3=1$ . For this choice, we have tension in the spring resulting in positive $\left[K_{\sigma }\right]$ (when compared with axial rod). In the second case, we choose $k_1=k_2=1$ and $k_3=-1$ , simulating a negative $\left[K_{\sigma }\right]$ (in analogy to axial rod).

Force F versus u graphs are shown in Figure 2. When $k_3$ is negative, the relationship between F and u shows progressively increasing displacement u with

increasing F. Additionally, the derivative $\left(\frac{\text{d}u}{\text{d}F}\right)$ also increases as F increases for all $F\le 1.3$ .

When $k_3$ is positive, increasing F results in progressively increasing positive $\left[K_{\sigma }\right]$ , resulting in progressively decreasing displacement u. When $k_3=-1$ ,

$\left(\frac{\text{d}u}{\text{d}F}\right)$ increases as F increases for all $F\le 1.3$ . In the vicinity of $F=1.3$ (between $F=1.3$ and $F=1.4$ ), u increases more sharply with increasing F indicating progressive loss of stiffness due to negative $\left[K_{\sigma }\right]$ . Between $1.3\le F\le 1.4$ we observe static bifurcation (a sudden drop in the value of u from largest positive to negative). The precise value of u at bifurcation can be determined if we differentiate (131) with respect to F.

$\frac{\text{d}u}{\text{d}F}=\frac{1}{k_1+2k_{2}u+3k_3{u}^2}$ (135)

$\frac{\text{d}u}{\text{d}F}$ in (135) is infinity (condition for bifurcation) when

$k_1+2k_{2}u+3k_3{u}^2=0.$ (136)

(a)(b)

Figure 2. (a) Model problem I: Force F versus u₁ plot for $-k_3$ and $+k_3$ ; (b) Model problem II: ${}_d{\sigma }_{11}^{\left[0\right]}$ versus u₁.

Using $k_1=1,k_2=1$ and $k_3=-1$ we obtain solutions $u=0$ and $u=\frac{4}{3}$ . Clearly, $u=0$ is not admissible, hence $\frac{\text{d}u}{\text{d}F}$ becomes infinite at $u=\frac{4}{3}=1.3333$ , approximately the same value of u at which bifurcation happens in Figure 2(a).

Static bifurcation in tension is not possible, and we do not observe one either in Figure 2(a). Values of u at which bifurcation occurs in model problem I is not possible to achieve in the axial rod (shown below) model problem, pointing unrealistic stiffness physics at bifurcation in model problem I when $k_1=k_2=1$ and $k_3=-1$ .

Model problem II

Solution of model problem II consisting of (131)-(134) can be obtained using finite element method based on GM/WF or using finite element method based on residual functional (see [6] [32] and Section 4.2.2). Finite element based on residual functional (unconditionally stable for nonlinear differential operators) is used to calculate solution of (131)-(134) using u₁ and ${}_d{\sigma }_{11}^{\left[0\right]}$ as dependent variables. A dimensionless rod of one unit is considered, fixed at left end ( $x_1=0$ ) and subjected to ${}_d{\sigma }_{11}^{\left[0\right]}$ at $x_1=1$ . Solutions are calculated for $-0.5<{}_d\sigma {}_{11}^{\left[0\right]}\le 2$ . Graph of ${}_d{\sigma }_{11}^{\left[0\right]}$ versus u₁ is shown in Figure 2(b). This graph is similar to a scaled graph of force F versus u₁.

For ${}_d{\sigma }_{11}^{\left[0\right]}\ge 0$ , the positive $\left[K_{\sigma }\right]$ contributes to progressive increase in the stiffness of the rod as ${}_d{\sigma }_{11}^{\left[0\right]}$ increases. This results in progressively reduced displacement u₁ for the same increment of ${}_d{\sigma }_{11}^{\left[0\right]}$ during loading.

On the other hand, when ${}_d{\sigma }_{11}^{\left[0\right]}$ is negative, the negative $\left[K_{\sigma }\right]$ contributes to a reduction in the total stiffness of the rod. This reduction in stiffness leads to progressively larger displacement for the same increment of ${}_d{\sigma }_{11}^{\left[0\right]}$ during loading.

In the vicinity of ${}_d{\sigma }_{11}^{\left[0\right]}=-0.5$ , there is total loss of stiffness. As a result, computations fail beyond ${}_d{\sigma }_{11}^{\left[0\right]}<-0.5$ . We note that when loading is compressive, there is total loss of stiffness in the rod for some value of ${}_d{\sigma }_{11}^{\left[0\right]}$ but unlike in model problem I, there is no static bifurcation. On the other hand in model I, there is no complete loss of stiffness but bifurcation exists. Existence of different physics in static loading in model problems I and II does influence the dynamic response and dynamic bifurcation as static stiffness remains the same in dynamic case.

7. Remarks

(1) Increase and decrease in total stiffness due to positive and negative $\left[K_{\sigma }\right]$ is clearly observed in model problem II. In case of model problem I, there is no $\left[K_{\sigma }\right]$ , but $k_3$ is assumed to simulate $\left[K_{\sigma }\right]$ .

(2) In model problem II, the occurrence of static bifurcation, as observed in model problem I, cannot be simulated under compressive loading. This is because the static bifurcation physics does not exist in model problem II.

(3) Main points of these studies are

(a) to demonstrate influence of $\left[K_{\sigma }\right]$ on stiffness, hence deformation and possible static bifurcation.

(b) Differences and similarities between model problems I and II: Since in model problem I, there is no space coordinate, u can be arbitrarily small or large. Thus, in model problem I, finite deformation physics can be simulated easily as it is purely phenomenological. Thus, existence or lack of static bifurcation physics is easier to simulate. We show from the studies using model problem II (in a later section) that u (or u₁) beyond a certain value (about 0.2 units in compression for a rod of one unit length) becomes non-physical. Is the total loss of stiffness in axial rod in compression, onset of static bifurcation? Perhaps we can argue either way. But, we observe no distinct static bifurcation in model problem II.

(c) Due to application of ${\sigma }_0\sin \left(\omega t\right)$ at the free end of the rod, the rod experiences a steady compression. As a consequence, $\left[K_{\sigma }\right]$ is negative. Thus in model problem II, we can only study negative $k_3$ physics of model problem I.

8. Numerical Studies for Model Problems I and II: Dynamic Response and Dynamic Bifurcation

8.1. Path Dependency in Nonlinear Dynamics

We discuss some points here that are helpful when analyzing the results of the model problems. In the absence of dissipation, the nonlinear dynamics physics is a reversible process. A reversible process cannot exhibit path dependency. For example, in static bifurcation, solutions are independent of the size of increment of load [1] [2] [3] [4] [5] . The same is true in nonlinear dynamics of inviscid media. Introduction of dissipation causes conversion of some mechanical work into entropy that is not reversible upon unloading. The entropy production is obviously path dependent. Thus, in nonlinear dynamics of TVES (without memory), path dependency is inherent. The severity of path dependency is of course dependent on entropy production. Higher entropy production leads to stronger path dependency. Since the mechanical work depends upon magnitude of $F_0$ (model problem I) or ${\sigma }_0$ (model problem II) and dissipation or dissipation coefficient c, appropriate combination of $F_0$ (or ${\sigma }_0$ ) and c is necessary to cause sufficient entropy generation that will result in visible path dependency. As a guide, low values of c with higher values of $F_0$ (or ${\sigma }_0$ ) and higher values of c with low values of $F_0$ ( ${\sigma }_0$ ) is likely to result in milder path dependency. We shall observe this in the model problem studies. In case of linear PDEs, presence of irreversibility does not introduce path dependency (Figure 8(a), Figure 8(b)). For path dependency, the mathematical model has to be nonlinear with presence of irreversibility.

In the following sections we present numerical studies for model problems I and II.

8.2. Model Problem I: Numerical Studies

Evolution for model problem I is calculated using Newmark linear acceleration method with Newton’s linear method for each increment of time to obtain converged solution. Since (64) is a second order nonlinear ODE in time; we use following ICs.

$u\left(x_1,0\right)=0;\stackrel{\dot{}}{u}\left(x_1,0\right)=0\forall x_1\in \left[0,L\right]=\left[0,1\right]$ (137)

Choice of $f_0\mathrm{,}{c}_1\mathrm{,}{c}_2\mathrm{,}{c}_3\mathrm{,}{k}_1\mathrm{,}{k}_2\mathrm{,}{k}_3$ and $\omega$ are discussed in the following. We consider $\omega \in \left[{\omega }_1\mathrm{,}{\omega }_n\right]$ in increments of 0.1. For a circular frequency $\omega$ , $\Delta t$ is

chosen as $\frac{1}{20}\left(\frac{2\pi }{\omega }\right)$ i.e., $\Delta t$ is $\frac{1}{20}$ of the time period of $\omega$ . Computation of

(a)(b)(c)

Figure 3. (a) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.2$ , $f_0=0.014$ ; (b) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.1$ , $f_0=0.014$ ; (c) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.08$ $f_0=0.014$ .

evolution is commenced at ${\omega }_1$ and continued till the response becomes cyclic (not necessarily periodic). At this point, maximum value of u is recorded for ${\omega }_1$ . Frequency ${\omega }_1$ is changed to ${\omega }_2$ maintaining the final state of deformation at ${\omega }_1$ as initial state at ${\omega }_2$ . Computations are continued till the response is cyclic for ${\omega }_2$ . At this point, maximum value of u is recorded for ${\omega }_2$ . This process is continued for progressively increasing frequencies till we reach ${\omega }_n$ . This frequency sweep from ${\omega }_1\to {\omega }_n$ is referred to as $L\to R$ (or ${\omega }_1\to {\omega }_n$ ). To determine any path dependencies in the frequency response, the computations are started at ${\omega }_n$ with progressively decreasing frequencies in the same manner as for ${\omega }_1\to {\omega }_n$ . We refer to these computations with progressively decreasing frequencies as right to left $L\leftarrow R$ (or ${\omega }_n$ to ${\omega }_1$ , ${\omega }_1\leftarrow {\omega }_n$ ).

Choice of $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3$ and $k_1\mathrm{,}{k}_2\mathrm{,}{k}_3$ is arbitrary as these are not related to any particular physical model. But, their choices undoubtedly influence the frequency response curve. We do not have any particular experimental frequency response in mind that we want to match using appropriate choices of $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3$ and $k_1\mathrm{,}{k}_2\mathrm{,}{k}_3$ . Thus, in most studies we adhere to the unit values for these parameters. In model problem II, the material coefficients in $\left[C_i\right]$ and $\left[K\right]$ are actual material coefficients, but after space-time decoupling, they are multiplied by the approximation functions, their derivatives, etc. As a result, if one chooses only one ODE from the space-time decoupled mathematical model, then surely we expect $k_1\mathrm{,}{k}_2\mathrm{,}{k}_3$ and $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3$ to have different values than the material coefficients. Unfortunately, in model problem I, we do not have any basis to decide their numerical values other than $k_1=k_2=k_3$ , $c_1=c_2=c_3$ . However, there is one exception: if we want to simulate negative $\left[K_{\sigma }\right]$ and capture the physics similar to compression in the axial rod, we need to choose a negative value for $k_3$ .

8.2.1. Case (a): Simulations for Negative $\left[K_{\sigma }\right]$ i.e., $k_3=-1$ .

We choose $k_1=k_2=1$ but $k_3=-1$ to simulate negative $\left[K_{\sigma }\right]$ (as in case of axial rod). In this section, we present two studies. In the first study, we choose a fixed $f_0$ and progressively reduce $c_1=c_2=c_3$ till computations are no longer possible due to singular coefficient matrix. In the second study, we choose fixed values of $c_1=c_2=c_3$ and vary $f_0$ till computations cease due to singular coefficient matrix in the computations. These studies are aimed at determining if bifurcation is possible when $k_3=-1$ with the choice of coefficients $c_i$ and $k_i$ considered here.

In the first study, we choose $F_0=0.014$ and perform computations for progressively decreasing $c_1=c_2=c_3=0.2,0.1,0.08$ . The evolutions are computed for progressively increasing frequencies ${\omega }_1\to {\omega }_n$ ( $L\to R$ ) as well as for progressively decreasing frequencies $L\leftarrow R$ (or ${\omega }_1\leftarrow {\omega }_n$ ) of excitation in increment

of 0.1. We choose integration time step $\Delta t=\frac{1}{20}\left(\frac{2\pi }{\omega }\right)$ . Evolution is computed

for 25 cycles of each frequency. We remark that due to negative shift caused by compression, the positive and the negative peaks may differ, especially in the nonlinear range when deformation is finite. For each frequency, we record maximum negative as well as maximum positive peak values after the response becomes cyclic. In the frequency response graphs, absolute values of peaks are reported. Figures 3(a)-(c) show frequency response for progressively changing frequencies from $L\to R$ and $L\leftarrow R$ for $c_1=c_2=c_3=0.2,0.1,0.8$ . Maximum positive peaks as well as absolute value of maximum negative peaks are shown in the graphs. In the second study we choose $c_1=c_2=c_3=0.1$ and progressively increase $f_0=0.025,0.03,0.035$ and compute frequency response for

$L\to R$ ( ${\omega }_1\to {\omega }_n$ ) and $L\leftarrow R$ ( ${\omega }_1\leftarrow {\omega }_n$ ) using $\Delta t=\frac{1}{20}\left(\frac{2\pi }{\omega }\right)$ integration

time steps. For each frequency, we compute 25 cycles, found sufficient for the motion to be cyclic. Figures 4(a)-(c) show frequency response for $c_1=c_2=c_3=0.1$ and $f_0=0.025,0.03$ and 0.035.

Discussion of results

(1) From Figures 3(a)-(c) we note that progressively reduced values of $c_1=c_2=c_3$ for fixed $f_0$ yield progressively increasing maximum amplitude due to progressively reduced resistance to motion. The positive and the negative peak differ almost insignificantly. We do observe path dependency for all three values of $c_1=c_2=c_3$ . Path dependency increases with reducing values of $c_1=c_2=c_3$ . For $f_0=0.014$ and $c_1=c_2=c_3=0.2,0.1,0.08$ we observe no bifurcation in Figures 3(a)-(c). We remark that for $f_0=0.014$ , choice of $c_1=c_2=c_3=0.08$ is limiting value of damping. Below this value, computations fail. Computed solutions show no existence of dynamic bifurcation for these choices of damping coefficients and $f_0$ . This study confirms that presence of dynamic bifurcation requires appropriate combination of various aspects of the physics such that an unstable condition exists which forces a sudden change to restore stability. Evidently, in this case the choices of parameters only lead to stable deformation configurations.

(a)(b)(c)

Figure 4. (a) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.1$ , $f_0=0.025$ ; (b) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.1$ , $f_0=0.030$ ; (c) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.1$ , $f_0=0.035$ .

(2) From Figures 4(a)-(c), we note that for fixed $c_1=c_2=c_3=0.1$ , progressively increasing values of $f_0=0.025,0.03,0.035$ results in progressively increasing maximum amplitude of motion, as expected. For $f_0=0.025$ (Figure 4(a)) positive and negative peaks for each path show some differences. Positive peaks and the negative peaks for $L\to R$ and $L\leftarrow R$ are relatively close together. At $f_0=0.030$ these differences increase (Figure 4(b)). We also note existence of mild bifurcation (more rapid changes) that is more pronounced in both peaks for $L\leftarrow R$ path. $L\to R$ path also shows mild bifurcation but at a different frequency than for $L\leftarrow R$ path. Path dependency can be clearly observed. For $f_0=0.035$ (Figure 4(c)) distinct dynamic bifurcation is observed to the left of the peaks. In this case, there is virtually no path dependency i.e., $L\to R$ and $L\leftarrow R$ paths yield same frequency response. We also point out that in this case also for $c_1=c_2=c_3=0.1$ computations fail if $f_0$ is increased beyond 0.035.

(3) Studies in (1) and (2) show that when $k_3=-1$ mild dynamic bifurcation can be simulated but more distinct and sharper change in peaks is difficult to simulate due to the fact that for very low values of damping and for very high value of $f_0$ (both needed for dynamic bifurcation) computations fail.

(4) We recall that static bifurcation for this model problem I, with $k_1=k_2=1$ and $k_3=-1$ occurs when $1.3\le u\mathrm{,}f\le 1.4$ . Displacement value of u in this range is not possible for the choices of coefficients considered here for model problem I.

(5) In this phenomenological model, we do not have any basis for choosing $k_1\mathrm{,}{k}_2\mathrm{,}{k}_3$ and $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3$ , hence the reason for choosing the values given in the studies presented here.

(6) We note that for negative $\left[K_{\sigma }\right]$ , bifurcation occurs to the left of the peak in all studies.

8.2.2. Case (b): Simulations for Positive $\left[K_{\sigma }\right]$ i.e., $k_3=1$ : Model Problem I

We consider solutions of model problem I for positive $\left[K_{\sigma }\right]$ i.e., $k_3=1$ . In this case, the influence of the stress field is to increase total stiffness. We choose $k_1=k_2=k_3=1$ in all studies presented in this section and consider different choices of $c_1=c_2=c_3$ and $f_0$ . In this case, increasing stress field, increasing non linearity and increasing dissipation all contribute to increase in stiffness. The only mechanisms of reduction in stiffness is dynamic motion i.e., inertial physics and reduced dissipation. In the vicinity of bifurcation point, an unbalance between the two is followed by sudden change in the energy state causing dynamic bifurcation. As shown earlier in the static bifurcation studies, an axial rod in tension can never experience total loss of stiffness at any loading that produces positive $\left[K_{\sigma }\right]$ . Furthermore, in case of axial rod, model problem II, positive $\left[K_{\sigma }\right]$ cannot be simulated due to $f_0\sin \left(\omega t\right)$ acting at the free end (right) of the rod. This is another difference in the two model problems that confirms the fact that model problem I is not representative of model problem II.

In the first study, we choose damping coefficients $c_1=c_2=c_3=0.1$ and $f_0=0.1,0.05,0.01$ . For each value of $f_0$ , frequency response is obtained for the frequencies ${\omega }_1\to {\omega }_n$ ( $L\to R$ ) and for ${\omega }_1\leftarrow {\omega }_n$ ( $L\leftarrow R$ ) in increments of 0.1. For each $f_0$ , at each frequency, positive and negative peaks are recorded

when the motion is steady and cyclic. Integration time step $\Delta t=\frac{1}{20}\left(\frac{2\pi }{\omega }\right)$ is

used in all calculations. Evolution is computed till completion of 25 cycles for each frequency for motion to be repetitive cyclic motion. Frequency responses for this study are shown in Figures 5(a)-(c) for $f_0=0.1,0.05$ and 0.01.

(1) Figure 5(a) shows frequency response for $f_0=0.1$ . For $L\to R$ path, positive and negative peaks differ. Negative peaks being always higher, but both naturally show bifurcation at the same frequency. Negative peak always being higher, hence results in more distinct appearances of dynamic bifurcation. For $L\leftarrow R$ path, we also observe dynamic bifurcation (not as distinct as in case of $L\to R$ ) but at a different frequency than for $L\to R$ . Path dependency can be clearly observed in all three figures.

(2) Frequency response for $f_0=0.05$ and $f=0.01$ are shown in Figure 5(b) and Figure 5(c). We observe no existence of dynamic bifurcation for both force values. Low force values result in very mild nonlinearity, thus eliminating possibility of the existence of dynamic bifurcation. In Figure 5(b) and Figure 5(c), we observe no path dependency and the differences between the negative and positive peaks reduce with progressively reducing force. The results indicate only mild nonlinearity, not sufficient for dynamic bifurcation and path dependency.

(3) This study shows that for a fixed dissipation physics, existence of dynamic bifurcation at higher value(s) of f₀ can be eliminated by lowering f₀ values, thus reducing presence of nonlinearity and stiffness. Likewise, for a fixed dissipation physics progressively increasing f₀ could likely result in dynamic bifurcation phenomenon.

(a)(b)(c)

Figure 5. (a) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.1$ , $f_0=0.1$ ; (b) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.1$ , $f_0=0.05$ ; (c) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.1$ , $f_0=0.01$ .

8.2.3. Case (c)

In this second study, we keep $k_1=k_2=k_3=1$ and $f_0=0.3$ fixed and choose $c_1=c_2=c_3=1.0,0.5$ and 0.1 i.e., progressively reducing dissipation. Frequency response for all three values of damping for $L\to R$ and $L\leftarrow R$ frequency range are shown in Figures 6(a)-(c). From Figure 6(a) and Figure 6(b), we note that for $c_1=c_2=c_3=1$ and 0.5, negative and positive peaks for $L\to R$ and $L\leftarrow R$ are almost identical but we clearly observe path dependency i.e., $L\to R$ and $L\leftarrow R$ frequency response are different for both damping coefficient values. We also note absence of bifurcation for both damping coefficient values. Clearly, much higher peak amplitudes are observed for 0.5 damping coefficient compared to 1.0. Figure 6(c) shows frequency response for $c_1=c_2=c_3=0.1$ . First, we note that maximum amplitude in this case is roughly four times that of for $c_1=c_2=c_3=1.0$ . We clearly observe bifurcation for both $L\to R$ and $L\leftarrow R$ loading paths. Negative peaks are always higher than positive (as explained earlier) and both peaks show presence of bifurcation in case of $L\to R$ bifurcation occurs at a higher frequency compared to $L\leftarrow R$ . Both positive and negative peaks at bifurcation in case of $L\to R$ are higher than the corresponding peaks for $L\leftarrow R$ . As expected, since the positive $k_3$ implies positive $\left[K_{\sigma }\right]$ , the bifurcation occurs to the right of the peak amplitude.

In this study, we clearly see the role of $c_1=c_2=c_3$ in the bifurcation physics. For fixed magnitude of excitation force, progressively reduced values of damping enables existence of bifurcation. Progressively reducing damping results in enhanced motion, hence increased nonlinearity that influences $\left[K_{\sigma }\right]$ due to stress field.

(a)(b)(c)

Figure 6. (a) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=1.0$ , $f_0=0.3$ ; (b) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.5$ , $f_0=0.3$ ; (c) Model problem I: Frequency ω versus maximum amplitude u₁ response for $c_1=c_2=c_3=0.1$ , $f_0=0.3$ .

8.3. Model Problem II: Numerical Studies

We consider dimensionless form of model II (Equations (128)-(130)). First, we write (128)-(130)) using hat ( ${}^{\wedge }$ ) on all quantities indicating that they have their usual dimensions or units.

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

${}_d{\hat{\sigma }}_{11}^{\left[0\right]}-\hat{E}\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)-{\hat{c}}_1\left(\frac{\partial {\hat{v}}_1}{\partial {\hat{x}}_1}+\frac{\partial {\hat{u}}_1}{\partial {\hat{x}}_1}\frac{\partial {\hat{v}}_1}{\partial {\hat{x}}_1}\right)=0$ (139)

${\hat{v}}_1-\frac{\partial {\hat{u}}_1}{\partial \hat{t}}=0$ (140)

We introduce the following reference quantities (with subscript zero or “ref”) and dimensionless quantities:

$\begin{array}{l}\begin{array}{l}{\rho }_0=\frac{{\hat{\rho }}_0}{{\left({\rho }_0\right)}_{\text{ref}}}\mathrm{<mo>\; ;\; </mo>}{u}_1=\frac{{\hat{u}}_1}{L_0}\mathrm{<mo>\; ;\; </mo>}{v}_1=\frac{{\hat{v}}_1}{v_0}\mathrm{<mo>\; ;\; </mo>}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}^{\left[0\right]}}{{\sigma }_0};E=\frac{\hat{E}}{E_0}\\ {\sigma }_0=E_0={\rho }_0{v}_0^2;\underset{˜}{c}=\frac{\hat{c}}{c_0};\text{wechoose}{\rho }_0=0.001\end{array}\}\hfill \end{array}$ (141)

Using (141) in (138)-(140), we obtain the dimensionless form of the equations (138)-(140).

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

${}_d{\sigma }_{11}^{\left[0\right]}-E\left(\frac{\partial u_1}{\partial x_1}+\frac{1}{2}{\left(\frac{\partial u_1}{\partial x_1}\right)}^2\right)-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$ (143)

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

$\forall \left(x\mathrm{,}t\right)\in {\bar{\Omega }}_{xt}={\Omega }_x\times {\Omega }_t=\left(\mathrm{0,}L\right)\times \left(\mathrm{0,}\tau \right)$ in which $E=1$ if we choose $E_0=\hat{E}$ and $c=\left(\frac{\hat{c}}{{\sigma }_0}\frac{v_0}{L_0}\right)=\frac{\hat{c}}{{\rho }_0{v}_0{L}_0}=c\left(\frac{c_0}{{\rho }_0{v}_0{L}_0}\right)=\frac{c}{Re}$

c is dimensionless damping coefficient. Re is Reynolds number. We calculate solutions of (142)-(144) using space-time finite element based on space-time residual functional for a space-time strip with time marching. Figure 7(a) shows schematic of the rod of dimensionless length L (=1). Figure 7(b) shows space-time domain ${\bar{\Omega }}_{xt}={\bar{\Omega }}_x\times {\bar{\Omega }}_t=\left[0,L\right]\times \left[0,\tau \right]$ ; $\tau$ being final value of time. A discretization ${\bar{\Omega }}_{xt}={\displaystyle \underset{i}{\cup }{\bar{\Omega }}_{xt}^{\left(i\right)}}$ of the space-time domain ${\bar{\Omega }}_{xt}$ into space-time strips ${\bar{\Omega }}_{xt}^i;i=1,2,\cdots$ is shown in Figure 7(c).

Figure 7(d) shows a discretization ${\left({\bar{\Omega }}_{xt}^{\left(1\right)}\right)}^T={\displaystyle \underset{e}{\cup }{\bar{\Omega }}_{xt}^e}$ of the first space-time

step 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 7(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 for $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 ω is done at a value of time t when ${}_d{\sigma }_{11}^{\left[0\right]}=0$ , so that change in ω 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 observed confirming accurate evolution.

Figure 7. Schematic, discretization, BCs, and ICs for an axial rod (model problem II).

As in case of model problem I, 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.

8.3.1. Case (i)

First, to demonstrate validity and accuracy of computations we present frequency response for linear dynamic study corresponding to linear stiffness and linear damping. All nonlinear terms in (142)-(144) are set to zero. Figure 8(a) shows frequency response for $c=0.003$ and ${\sigma }_0=0.03$ (low damping, low force) in which the peak is at first fundamental frequency of the linear system, positive and negative peaks responses are identical confining small deformation, small strain physics. We observe no path dependency i.e., $L\to R$ and $L\leftarrow R$ frequency responses are the same.

Figure 8(b) shows frequency response for ${\sigma }_0=0.003$ and damping values of $c=0.003$ and 0.0035. For $c=0.0035$ , we observe lower peak values. Frequency response for $c=0.003$ and $c=0.0035$ only differs in the vicinity of the peak (but stay symmetric with respect to peak values in both cases) as expected. These studies are included to show that the computed evolution in the entire frequency range is accurate (residual functional I values $O\left({10}^{-8}\right)$ or lower). These studies confirm that when the mathematical model is linear, presence of irreversibility does not introduce path dependency in the solution.

8.3.2. Case (ii): Progressively Reduced Damping for a Fixed ${\sigma }_0$

In this study we choose a fixed value of ${\sigma }_0$ and compute frequency response for progressively reduced damping coefficient c until further reduction in the damping coefficient leads to singular matrices or lack of convergence of the Newton’s linear method. We choose ${\sigma }_0=0.03$ and $c=0.01,0.005,0.003,0.002$ .

Figures 9(a)-(d) show frequency response for $c=0.01,0.005,0.003$ and 0.002. For all four values of damping, ( $L\to R$ ), ( $L\leftarrow R$ ) frequency response are shown using positive and negative peaks. In Figures 9(a)-(c), we observe that for each value of $c=0.01,0.005,0.003$ that $L\to R$ and $L\leftarrow R$ yield same frequency response for negative as well as positive peaks, negative peak always being higher than positive peaks. With progressive reduction in c, we also observe increase in peak values of the response. Frequency response for $c=0.002$ in Figure 9(d) shows distinct bifurcation to the left of the peak (due to negative $\left[K_{\sigma }\right]$ ). Here also higher negative peak amplitude show more distinct drop in magnitude than positive peaks, but we clearly observe bifurcation in both positive and negative peaks. This study confirms that for a fixed ${\sigma }_0$ value, total lack of bifurcation at higher damping can be changed to distinct bifurcation by lowering the damping value. We remark that in some cases this process may result into failure of computations or lack of convergence of Newton’s linear method for lower value of damping before existence of bifurcation can be established.

(a)(b)

Figure 8. (a) Model problem II: Frequency ω versus maximum amplitude u response for linear stiffness and linear damping (c = 0.003); (b) Model problem II: Frequency ω versus maximum amplitude u response for linear stiffness and linear damping ( $c=0.0035$ and $c=0.003$ ).

For this choice of ${\sigma }_0$ , we could observe bifurcation before the computations cease. In Figure 9(d) we observe that $L\to R$ and $L\leftarrow R$ frequency response show mild path dependency. This phenomenon of mild path dependency has also been reported in ref [26] .

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

Figure 9. (a) Model problem II: Frequency ω vs Amplitude u response for $c=0.01$ and $c=0.005$ ; (b) Model problem II: Frequency ω vs Amplitude u response for $c=0.01$ and $c=0.005$ ; (c) Model problem II: Frequency ω vs Amplitude u response for $c=0.003$ ; (d) Model problem II: Frequency ω vs Amplitude u response for $c=0.002$ .

8.3.3. Case (iii): Progressively Increasing ${\sigma }_0$ for a Fixed Damping c

In this study, we choose a fixed value of damping coefficient $c=0.003$ and determine frequency response for progressively increasing ${\sigma }_0=0.025$ , 0.030, 0.045 and 0.055. The objective of this study is to demonstrate that for a fixed value of damping, presence of dynamic bifurcation is possible if ${\sigma }_0$ is increased beyond a threshold value. Details of computations and of generating frequency response curve remain the same as discussed in earlier studies. Figures 10(a)-(d) show

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

Figure 10. (a) Model problem II: Frequency ω vs Amplitude u response for ${\sigma }_0=0.030$ and 0.025; (b) Model problem II: Frequency ω vs Amplitude u response for ${\sigma }_0=0.030$ and 0.025; (c) Model problem II: Frequency ω vs Amplitude u response for ${\sigma }_0=0.045$ and 0.055; (d) Model problem II: Frequency ω vs Amplitude u response for ${\sigma }_0=0.045$ and 0.055.

frequency response for ${\sigma }_0=0.025$ , 0.035, 0.045, 0.055. In Figure 10(a) and Figure 10(b) for $\sigma =0.025$ and 0.030, $L\to R$ and $L\leftarrow R$ frequency response are almost the same for positive as well as negative peak and there is no evidence of dynamic bifurcation. As usual, negative peaks are higher than the positive peaks. In Figure 10(c), for ${\sigma }_0=0.045$ , we observe onset of dynamic bifurcation, more evident in the negative peak than the positive peak, but we observe lack of path dependency i.e., $L\to R$ and $L\leftarrow R$ frequency response are almost same. In Figure 10(d) for ${\sigma }_0=0.055$ we clearly observe bifurcation to the left of the positive and negative peaks, negative peak showing more distinct bifurcation due to larger maximum amplitude of motion. Since path dependency is due to irreversibility caused by damping, the path dependency is mild due to low damping value. It is combination of damping and force that determine irreversibility. High damping but low force and low damping but high force, both are likely to result in low irreversibility, hence milder path dependency.

In this study, we have seen that for a fixed c (0.003 in this study) there is a threshold value of ${\sigma }_0$ at which distinct bifurcation is observed. We clearly see that path dependency is not a necessary condition for the existence of bifurcation. Larger maximum amplitudes in Figure 10(c) and Figure 10(d) indicate that stronger nonlinearity, hence stronger influence of $\left[K_{\sigma }\right]$ is essential for bifurcation to exist. We also observe that positive peaks are not as effective as negative peaks in illustrating existence of bifurcation.

9. Nonlinear Dynamics of Continuous Solid Media When $\left[K_{\sigma }\right]$ Is Positive

In the model problem studies for axial rod (model problem II) subjected to harmonic excitation, $\left[K_{\sigma }\right]$ due to stress field is always negative resulting in reduction of total stiffness. In the type of applications, where $\left[K_{\sigma }\right]$ is negative, the inertial physics due to linear acceleration and $\left[K_{\sigma }\right]$ both result in reduction in stiffness. On the other hand, the stiffness matrices (all three of them) and dissipation mechanism offer resistance to the motion. Thus in such applications (in which $\left[K_{\sigma }\right]$ is negative), dynamic bifurcation is due to lack of balance between the two competing mechanism of stiffnesses. When there is unbalance between these two stiffness mechanisms, a state of instability, balance is restored by sudden change in energy state resulting in sudden change in the magnitude of motion. This of course is dynamic bifurcation. We have seen that when $\left[K_{\sigma }\right]$ is negative, the dynamic bifurcation is always to the left of the peak amplitude in the frequency response in case of model problem I as well as model problem II.

In many applications such as nonlinear dynamics of bending of fixed-fixed beams, bending of clamped plates, the midplane is always in tension, resulting in positive $\left[K_{\sigma }\right]$ , just like in case of axial rod where $\left[K_{\sigma }\right]$ is always negative for excitation at the end of the rod, in such applications $\left[K_{\sigma }\right]$ is always positive due to midplane always being in tension during dynamic motion. Positive $\left[K_{\sigma }\right]$ of course increases total stiffness. Thus, in this case reduced damping and increased inertial physics due to linear acceleration are the mechanisms that result in reduction of the total stiffness. When there is a lack of balance between the stiffness reduction and other mechanism causing increase in stiffness, dynamic bifurcation occurs to restore the balance. Based on negative $\left[K_{\sigma }\right]$ studies and studies for positive $k_3$ in case of model problem I, we expect that when $\left[K_{\sigma }\right]$ is positive, dynamic bifurcation will occur to the right of the peak values (negative as well as positive). Dynamic bifurcation studies for beams, plates and shells using mathematical models based on CBL of CCM and solution methods described in this paper are in progress and will be reported in a forthcoming paper.

10. Summary and Conclusions

(1) The CBL of CCM derived using contravariant second Piola-Kirchhoff stress and Green’s strain tensor augmented with consistent constitutive rate theories for deviatoric contravariant second Piola-Kirchoff stress tensor in terms of Green’s strain tensor and its rates up to order n completely define finite deformation, finite strain nonlinear dynamics of TVE solids. Unfortunately, the factors influencing nonlinear dynamics and in particular dynamic bifurcation are difficult to establish from this mathematical model.

If we substitute constitutive theory for ${}_d\sigma {}^{\left[0\right]}$ in BLM and only consider isothermal physics, then we obtain three PDEs in displacements $u_i$ from the balance of linear momenta. If we decouple space and time and consider integral form for a spatial discretization using fundamental lemma of the calculus of variations followed by Galerkin method with weak form in space, then we obtain a system of nonlinear ODEs (1) in time in $\left\{\delta \right\}\mathrm{,}\left\{\stackrel{\dot{}}{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ , the degrees of freedom for the spatial discretization and their first and second time derivatives. This is for the first time (due to GM/WF in space) we realize the concept of mass $\left[M\right]$ , stiffness $\left[K\right]$ and damping $\left[C_i\right]$ matrices. Only $\left[M\right]$ is symmetric. $\left[K\right]$ and $\left[C_i\right]$ are nonsymmetric and are often referred to as secant stiffness and damping matrices. We have shown that $\left[K_e\right]$ , the stiffness matrix for element e, can be decomposed into $\left[K_e^1\right]\mathrm{,}\left[K_e^2\right]$ and $\left[K_e^3\right]$ (Equations (48)-(50)), $\left[K_e^1\right]$ with constant coefficients and $\left[K_e^2\right]\mathrm{,}\left[K_e^3\right]$ with coefficients that are linear and quadratic functions of $\left\{{\delta }_e\right\}$ , dofs for an element e. The same decomposition can be applied to each of $\left[C_{\left[i\right]}\right]$ . At this stage, we have little more information, and nonlinear dynamics now depends upon stiffness $\left[K\right]$ and damping matrices $\left[C_{\left[i\right]}\right]$ whose coefficients are constant, linear and quadratic functions of dofs and of course on constant coefficient mass matrix $\left[M\right]$ . But, still this information is not sufficient to completely determine precisely the factors influencing dynamic bifurcation in nonlinear dynamics.

(2) The 1D phenomenological mass, spring, damper model (model problem I), Equation (64) is supposedly 1D equivalent to the model discussed in (1) (for $C_{\left[1\right]}$ only). This of course is not a valid representation of mathematical model (1). We remark that mathematical model (64) is a lumped model in space, thus cannot possibly describe deformation physics of a continuous solid media in which space and time are intrinsically coupled and always present. We point out that in mathematical model (1), the nonlinear coupled ODEs contain dofs and their derivatives that correspond to spatial positions. This is an essential feature of space-time decoupled description for a continuous media in which spatial details are preserved in the ODEs in time. This establishes without doubt that study of nonlinear dynamics using (64) may be interesting, but cannot possibly describe nonlinear dynamics of a continuous solid media.

(3) We have shown that when integrating (1) in time using a time integration method such a Newmark linear acceleration method and employing Newton’s linear method to obtain a converged solution for each time increment, we observe the appearance of tangent matrix $\left[H_T\right]$ whose inverse controls the incremental changes in the solution in Newton’s linear method. That is, computation of incremental solution in the computation of evolution requires inverse of $\left[H_T\right]$ , the tangent matrix defined by (105). Thus, the dynamic solution depends on $\left[H_T\right]$ that consists of symmetric $\left[M\right]\mathrm{,}\left[C_T\right]$ and $\left[K_T\right]$ . $\left[K_T\right]$ in terms of sum of $\left[K^1\right]\mathrm{,}\left[{\underset{˜}{K}}^2\right]\mathrm{,}\left[{\underset{˜}{K}}^3\right]$ and $\left[K_{\sigma }\right]$ and likewise $\left[C_T\right]$ is the sum of $\left[C^1\right]\mathrm{,}\left[{\underset{˜}{C}}^2\right]\mathrm{,}\left[{\underset{˜}{C}}^3\right]$ and $\left[C_{\sigma }\right]$ . $\left[M\right]\mathrm{,}\left[K^1\right]\mathrm{,}\left[C^1\right]$ have constant coefficients. Coefficients of $\left[{\underset{˜}{K}}^2\right]\mathrm{,}\left[{\underset{˜}{C}}^2\right]$ are linear functions of $\left\{\delta \right\}$ and those of $\left[{\underset{˜}{K}}^3\right]$ and $\left[{\underset{˜}{C}}^3\right]$ are quadratic functions of $\left\{\delta \right\}$ . $\left[K_{\sigma }\right]$ is the stiffness matrix due to presence of stress field. We note that $\left[{\underset{˜}{K}}^2\right]\mathrm{,}\left[{\underset{˜}{K}}^3\right]\mathrm{,}\left[{\underset{˜}{C}}^2\right]$ and $\left[{\underset{˜}{C}}^3\right]$ are not the same as $\left[K^2\right]\mathrm{,}\left[K^3\right]\mathrm{,}\left[C^2\right]\mathrm{,}\left[C^3\right]$ resulting from decoupling of space and time in which $\left[K^2\right]\mathrm{,}\left[C^2\right]$ are not symmetric. It is only after applying Newmark linear acceleration method with Newton’s linear method that we realize dependence of the solution on $\left[H_T\right]$ and thus on the new $\left[{\underset{˜}{K}}^2\right]\mathrm{,}\left[{\underset{˜}{K}}^3\right]\mathrm{,}\left[{\underset{˜}{C}}^2\right]\mathrm{,}\left[{\underset{˜}{C}}^3\right]\mathrm{,}\left[K_{\sigma }\right]$ and $\left[C_{\sigma }\right]$ matrices. Concept of $\left[K_{\sigma }\right]$ in this approach is instrumental in understanding bifurcation physics. Going back to (64), linear and quadratic stiffness terms maybe viewed as synonymous to $\left[K^2\right]$ and $\left[K^3\right]$ in (1), but they fail to explain the physics due to $\left[K_{\sigma }\right]$ . It is worth noting that $\left[K_T\right]$ is also realized in finite deformation, finite strain BVPs. Confirming that $\left[K_T\right]$ does not depend upon BVP or IVPs, it describes the physics resulting in different forms of stiffness contributions due to linear and nonlinear strain measure in Green’s strain tensor that contribute to the total stiffness and plays crucial role in static as well as dynamic bifurcations. Compressive stress field resulting in negative stresses causes negative $\left[K_{\sigma }\right]$ , hence reduces total stiffness. On the other hand, tensile stress field will result in positive $\left[K_{\sigma }\right]$ which would increase the total stiffness of the medium.

(4) Now we can identify the factors influencing nonlinear dynamics and dynamic bifurcation physics. Stiffness matrices $\left[K^1\right]$ $\left[{\underset{˜}{K}}^2\right]$ , $\left[{\underset{˜}{K}}^3\right]$ (all symmetric) contribute to total stiffness. Influence of the degree of nonlinearity is reflected in $\left[{\underset{˜}{K}}^2\right]$ and $\left[{\underset{˜}{K}}^3\right]$ . For large finite strain and large finite deformation, $\left[{\underset{˜}{K}}^3\right]$ plays significant role in influencing the stiffness. Dissipation mechanism (due to all damping matrices) provides resistance to motion, thus opposing the likelihood of dynamic instability (dynamic bifurcation). In order words, low damping assists existence of dynamic bifurcation and vice versa. The stiffness matrix $\left[K_{\sigma }\right]$ due to the stress field plays critical role in dynamic bifurcation, negative $\left[K_{\sigma }\right]$ enhances likelihood of the existence of dynamic bifurcation whereas positive $\left[K_{\sigma }\right]$ helps in inhibiting the likelihood of dynamic bifurcation due to increased stiffness. The last factor influencing dynamic bifurcation is inertial physics due to accelerations ( $\left[M\right]\left\{\ddot{\delta }\right\}$ ) which results in reducing stiffness, hence aids in the likelihood of the existence of dynamic bifurcation. Thus, in summary, on one hand, damping and stiffnesses $\left[K^1\right]\mathrm{,}\left[{\underset{˜}{K}}^2\right]\mathrm{,}\left[{\underset{˜}{K}}^3\right]$ and on the other hand inertial physics due to acceleration and negative $\left[K_{\sigma }\right]$ are two competing mechanisms in nonlinear dynamics and dynamic bifurcation physics. In stable dynamic response there is a balance between these two mechanisms. Any unbalance between these two mechanisms must result in sudden change in the deformation physics to restore stable balance, this of course is dynamic bifurcation. When $\left[K_{\sigma }\right]$ is positive, it enhances stiffness, hence in this case damping, stiffnesses and $\left[K_{\sigma }\right]$ equilibrate with inertial physics due to acceleration and the imbalance between these two mechanisms results in dynamic bifurcation. Thus, it is obvious that the nature of dynamic bifurcation for positive $\left[K_{\sigma }\right]$ and negative $\left[K_{\sigma }\right]$ must be different. We have demonstrated this in the model problem studies that this indeed is the case. In case of negative $\left[K_{\sigma }\right]$ due to compressive stress field (model problem II, also model problem I with $k_3=-1$ ) dynamic bifurcation is always to the left of the peak. When $\left[K_{\sigma }\right]$ is positive, the dynamic bifurcation is to the right of the peak (shown for model problem I with $k_3=1$ ).

(5) 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. Some of these are summarized in the following.

(a) We have established that 1D phenomenological mathematical model like (64) 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 (64) do not correspond to negative $\left[K_{\sigma }\right]$ and positive $\left[K_{\sigma }\right]$ . 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 modal basis determined from the linear form of ODEs (1) in which $\left[K\right]$ is constant. We have presented details in this paper to show that forcing this transformation on (1) followed by use of Rayleigh damping (see reference [33] [34] ) for issues with Rayleigh damping) leads to a system of decoupled ODEs that neither describe 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 straight forward space-time coupled or space-time decoupled methods of obtaining their solutions, like we have presented in this paper.

(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 superposition of linear modes of vibrations using modal participation factors that are determined using nonlinear ODEs in the participation factors, equation (31). 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 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 combination with undetermined coefficients satisfies 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 homogenous 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).

(f) We have shown in model problem II that static bifurcation does not exist for axial rod. Upon increasing compressive force at the right end, eventually negative $\left[K_{\sigma }\right]$ becomes equal to the stiffness ( $\left[K^1\right]+\left[K^2\right]+\left[K^3\right]$ ) resulting in total loss of stiffness, hence failure of computations. But, the static bifurcation does exist as we have shown in the paper. Nonetheless, dynamic bifurcation does exist in the axial rod problem. Thus, static bifurcation is not a necessary condition for dynamic bifurcation to exist. On the other hand in model problem I, we have static bifurcation at $u=1.3$ when $k_3=-1$ and we also have dynamic bifurcation as well, confirming that in some applications both static and dynamic bifurcation can exist.

(g) We draw some conclusions and make some final remarks based on the theory and model problem studies for model problem II based on CBL of CCM and consistent constitutive theories.

(i) In this paper, we have identified factors influencing dynamic bifurcation, (b) we have shown that both in the static bifurcation and the dynamic bifurcation are controlled by tangent stiffness matrix $\left[K_T\right]$ .

(ii) In static bifurcation, negative $\left[K_{\sigma }\right]$ is essential (compression). Positive $\left[K_{\sigma }\right]$ can never result in static bifurcation. Since inertial physics due to acceleration is absent in static bifurcation, it is a balance between $\left[K_{\sigma }\right]$ and the rest of the stiffnesses.

(iii) Dynamic bifurcation can occur with positive and negative $\left[K_{\sigma }\right]$ as in this case in addition to $\left[K_{\sigma }\right]$ the inertial physics due to acceleration (lowering stiffness) is present.

(iv) When $\left[K_{\sigma }\right]$ is negative, dynamic bifurcation occurs to the left of the peak and when $\left[K_{\sigma }\right]$ is positive, the dynamic bifurcation occurs to the right of the peak.

(v) Static bifurcation is not a necessary condition for dynamic bifurcation to exist.

(vi) Phenomenological 1D models in time such as (64) can never describe the nonlinear dynamics of continuous solid media.

(vii) We have shown that in the study of nonlinear dynamics and dynamic bifurcation physics of TVE solids, the mathematical models must consist of CBL of CCM and the constitutive theories derived using entropy inequality and the representation theorem. This is the only mathematical model that is thermodynamically and mathematically consistent in the study of nonlinear dynamics.

(viii) The space-time coupled and the space-time decoupled finite element methods presented in the paper are two viable and mathematically sound approaches of obtaining solutions of the mathematical models is the best approach for obtaining solutions of IVPs in nonlinear dynamics without introducing any approximation in them.

(ix) A space-time coupled finite element method based on space-time residual functional, the space-time integrals are space-time variationally consistent [32] . This methodology leads to unconditionally stabled computations during the entire evolution. This method when used for a space-time strip with time marching (presented in this paper) is highly efficient. 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 equations in the mathematical model are satisfied accurately at each point in the space-time domain. This ensures that calculated evolution is almost time accurate. This approach is used in calculating harmonic response for model problem II, thus the reported solutions, through numerical, undoubtedly ensure time accuracy and satisfy equations in the point wise sense.

(x) In space-time decoupled methods, approximations due to space-time decoupling, stability and accuracy of time integration schemes and over all inability of this approach to quantitatively judge if the equations (PDEs) in the mathematical model are satisfied in the point wise sense, makes this approach less attractive and in fact inferior to space-time coupled finite element method based on space-time residual functional discussed and used here.

(xi) The solutions of nonlinear PDEs with irreversible physics exhibit path dependency. In nonlinear dynamics of TVES (without memory) some mechanical work is used or converted into entropy production (due to dissipation physics). Thus, solutions of IVPs in TVES with finite deformation, finite strain will always exhibit path dependency. Severity of path dependency depends upon entropy production, hence for some choice of parameters, it is possible that significant path dependency may not be observed in the solution.

The paper presents thermodynamically and mathematically consistent mathematical models for nonlinear dynamics and dynamic bifurcation physics and method of obtaining their solution that are almost time accurate and have no issues of stability. Model problem studies using model problem II (CBL of CCM) illustrates all of the features discussed in item (viii) and (ix). Whether the solutions obtained using this approach match experimentally measured response is another issue. If they do, we have a computational tool that avoids experiments. If they do not, then obviously either the physics considered in the two differs or there are other sources of errors that are contaminating the results. We remark that in the computational approach based on (viii) and (ix), the only source of confusion can be, lack of consideration of some physics in the mathematical model that is present in the experiment. Computational technique in (ix) has no flaws when the space-time residual functional is $O\left({10}^{-8}\right)$ or lower for each space-time strip. In view of the consistent mathematical model based on CBL of CCM and rigorous computational infrastructure, we have discussed and presented in the paper for studying nonlinear dynamics and dynamic bifurcation physics, 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.

Conflicts of Interest

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

References

[1]	Surana, K.S. (1983) Geometrically Non-Linear Formulation for Two Dimensional Curved Beam Elements. Computers & Structures, 17, 105-114. https://doi.org/10.1016/0045-7949(83)90035-4
[2]	Mescall, J.F. (1965) Large Deflections of Spherical Shells under Concentrated Loads. Journal of Applied Mechanics, 32, 936-938. https://doi.org/10.1115/1.3627340
[3]	Surana, K.S. (1982) Geometrically Nonlinear Formulation for the Axisymmetric Shell Elements. International Journal for Numerical Methods in Engineering, 18, 477-502. https://doi.org/10.1002/nme.1620180402
[4]	Surana, K.S. (1983) Geometrically Nonlinear Formulation for the Curved Shell Elements. International Journal for Numerical Methods in Engineering, 19, 581-615. https://doi.org/10.1002/nme.1620190409
[5]	Ramm, E. (1977) A Plate/Shell Element for Large Deflections and Rotations. MIT Press, Cambridge, Chapter 10.
[6]	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
[7]	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
[8]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part I. Archive for Rational Mechanics and Analysis, 33, 249-267. https://doi.org/10.1007/BF00281278
[9]	Wang, C.C. (1969) On Representations for Isotropic Functions, Part II. Archive for Rational Mechanics and Analysis, 33, 268-287. https://doi.org/10.1007/BF00281279
[10]	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. https://doi.org/10.1007/BF00272241
[11]	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
[12]	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
[13]	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
[14]	Spencer, A.J.M. (1971) Chapter 3. 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
[15]	Surana, K.S. (2022) Classical Continuum Mechanics. 2nd Edition, CRC/Taylor and Francis, Boca Raton. https://doi.org/10.1201/9781003105336
[16]	Stoker, J.J. (1950) Nonlinear Vibrations in Mechanical and Electrical Systems. Interscience, New York.
[17]	Chattarjee, A. (2009) A Brief Introduction to Nonlinear Vibrations.
[18]	Cammarano, A., Hill, T.L., Neild, S.A. and Wagg, D.J. (2014) Bifurcations of Backbone Curves for Systems of Coupled Nonlinear Two-Mass Oscillators. Nonlinear Dynamics, 77, 311-320. https://doi.org/10.1007/s11071-014-1295-3
[19]	Kovacic, I. and Rand, R. (2013) Straight-Line Backbone Curve. Communications in Nonlinear Science and Numerical Simulation, 18, 2281-2288. https://doi.org/10.1016/j.cnsns.2012.11.031
[20]	Londoo, J.M., Cooper, J.E. and Neild, S.A. (2017) Identification of Systems Containing Nonlinear Stiffnesses Using Backbone Curves. Mechanical Systems and Signal Processing, 84, 116-135. https://doi.org/10.1016/j.ymssp.2016.02.008
[21]	Peter, S., Riethmüller, R. and Leine, R. (2016) Tracking of Backbone Curves of Nonlinear Systems Using Phaselocked-Loops. In: Kerschen, G., Ed., Nonlinear Dynamics, Volume 1, Springer International Publishing, Cham, 107-120. https://doi.org/10.1007/978-3-319-29739-2_11
[22]	Renson, L., Gonzalez-Buelga, A., Barton, D.A.W. and Neild, S.A. (2016) Robust Identification of Backbone Curves Using Control-Based Continuation. Journal of Sound and Vibration, 367, 145-158. https://doi.org/10.1016/j.jsv.2015.12.035
[23]	Amabili, M., Balasubramanian, P. and Ferrari, G. (2016) Travelling Wave and Non-Stationary Response in Nonlinear Vibrations of Water-Filled Circular Cylindrical Shells: Experiments and Simulations. Journal of Sound and Vibration, 381, 220-245. https://doi.org/10.1016/j.jsv.2016.06.026
[24]	Alijani, F., Amabili, M., Balasubramanian, P., Carra, S., Ferrari, G. and Garziera, R. (2016) Damping for Large Amplitude Vibrations of Plates and Curved Panels, Part 1: Modeling and Experiments. International Journal of Non-Linear Mechanics, 85, 23-40. https://doi.org/10.1016/j.ijnonlinmec.2016.05.003
[25]	Amabili, M., Alijani, F. and Delannoy, J. (2016) Damping for Large-Amplitude Vibrations of Plates and Curved Panels, Part 2: Identification and Comparisons. International Journal of Non-Linear Mechanics, 85, 226-240. https://doi.org/10.1016/j.ijnonlinmec.2016.05.004
[26]	Amabili, M. (2018) Nonlinear Damping in Large-Amplitude Vibrations: Modelling and Experiments. Nonlinear Dynamics, 93, 5-18. https://doi.org/10.1007/s11071-017-3889-z
[27]	Breunung, T. and Haller, G. (2018) Explicit Backbone Curves from Spectral Submanifolds of Forced-Damped Nonlinear Mechanical Systems. Proceedings of the Royal Society A, 474, Article ID: 20180083. https://doi.org/10.1098/rspa.2018.0083
[28]	Ghorbel, B., Zghal, A., Abdennadher, M., Walha, L. and Haddar, M. (2020) Analysis of Strongly Nonlinear Systems by Using HBM-AFT Method and Its Comparison with the Five-Order Runge-Kutta Method: Application to Duffing Oscillator and Disc Brake Model. International Journal of Applied and Computational Mathematics, 6, Article No. 50. https://doi.org/10.1007/s40819-020-0803-z
[29]	Balasubramanian, P., Ferrari, G. and Amabili, M. (2017) A Tool to Identify Damping during Large Amplitude Vibrations of Viscoelastic Structures. Proceedings of the ASME 2017 International Mechanical Engineering Congress and Exposition. Volume 4B: Dynamics, Vibration, and Control, Tampa, 3-9 November 2017, V04BT05A006.1. https://doi.org/10.1115/IMECE2017-70773
[30]	Amabili, M. (2008) Nonlinear Vibrations and Stability of Shells and Plates. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511619694
[31]	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
[32]	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
[33]	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
[34]	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