Finite Deformation, Finite Strain Nonlinear Dynamics and Dynamic Bifurcation in TVE Solids ()
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
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
with
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
and their time derivatives for the spatial discretization.
(1)
The mass matrix
is a constant coefficient nondiagonal matrix (often called consistent mass matrix as it is consistent with the local approximation of the displacement field).
is the damping matrix corresponding to strain rate
(material time derivative of order i of the Green’s strain tensor
) and
is the stiffness matrix. Coefficients of
and
are up to quadratic functions of the displacement gradients.
and
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
and
can be additively decomposed into (i.e., consist of):
(2)
(3)
in which coefficients of
are constant, while the coefficients of
and
are linear and quadratic functions of
. By substituting (2) and (3) in (1) and regrouping terms, we can obtain:
(4)
If we choose
to define the right hand side of (4) and choose
with
, then (4) is written as:
(5)
In (5), the mass, damping and stiffness matrices have constant coefficients. All nonlinearities are lumped into
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
, whereas in (4), both linear and nonlinear stiffness and damping terms are clearly identifiable. We remark that even though the
and
in (5) have constant coefficients, the
and
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
are clearly identified as we have done in the subsequent sections of this paper. Additionally, there are no reported solutions of (5) in which
and
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),
and
correspond to nonlinear ODEs in time. We consider linear undamped eigenvalue problem corresponding to linear ODEs in time describing linear dynamics with damping.
(6)
Let the columns of
represent the mass normalized eigenvectors of the undamped eigenvalue problem corresponding to (6), then (6) can be transformed to modal basis using
(7)
and
(8)
(9)
This procedure is standard [32] . However, transformation to modal basis using (7)~(9) cannot be applied to (5) as
and
in (5) are not the same as
and
in (6) and (7)~(9). We can only apply modal basis transformation to the following nonlinear PDEs:
(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
and
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
to obtain the following ODEs in time.
(11)
If
are the eigenvalues of the eigenvalue problem with linear
and
, then (11) reduces to
(12)
in which
identity and
is a diagonal matrix containing
on the diagonal.
We note that in the modal basis,
is nondiagonal as is the matrix
in (7). Rayleigh damping is used to diagonalize the matrix
. In the transformation (7)~(9) to modal basis,
and
are modal participation factors. This implies the actual response
and
of the nonlinear dynamics problem consists of the superposition of the mass normalized eigenvectors of the linear problem using modal participation factors
and
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
is a linear combination of
and constant coefficient
i.e.
(13)
where
and
are functions of
(natural frequencies of the linear system). Transforming (13) into modal basis, we have:
(14)
We note that right hand side of (14) is a diagonal matrix containing
on the diagonals. Thus,
has been diagonalized using (13) in modal basis (14).
If we assume that
is the dimensionless damping coefficient corresponding to frequency
, then based on 1D spring-mass-damper system, we can write
(15)
In (15), the diagonal entries of
are given by
in the modal basis. Thus, we have:
(16)
Equation (16) can be used to obtain
from
versus
experimental data (only two data points are needed to determine
and
). By substituting (15) into (11), we obtain
(17)
in which
is a diagonal damping matrix containing
on its diagonal.
Equations (17) represent a decoupled system of “n” second order nonlinear ODEs in time. We can also express (17) as follows:
(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
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
. This suggests that calculating modal participation factors
using nonlinear decoupled ODEs and then using them in (7)~(9) is sufficient for (7)~(9) to yield
and
of nonlinear problem instead of
of linear problem.
(3) Use of superposition for a nonlinear problem raises the question regarding the validity of the resulting solutions
and
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
in (1). Assuming that
or
is proportional to
and
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
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
) 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
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:
(19)
(20)
(21)
(22)
(23)
(24)
(25)
Constitutive theories for
and
based on integrity
(26)
(27)
Simplified constitutive theories for
and
(in Voigt’s notation) that are linear in tensors (
and
) are given by:
(28)
(29)
in which
are displacements,
are body forces per unit mass,
is contravariant second Piola-Kirchhoff stress tensor, e is internal energy density,
is heat vector,
and
are equilibrium and deviatoric contravariant second Piola-Kirchhoff stress tensors.
is equation of state,
is mechanical pressure.
and
are combined generators and invariants of the argument tensors
and
of the constitutive tensor
,
is a temperature gradient tensor (argument tensor of constitutive tensor
),
and
are elastic constants,
and
are dissipation coefficients for strain rate
. This mathematical model has closure. It consists of thirteen equations BLM(3), FLT(1), constitutive theories for
(6),
(3) in thirteen dependent variables
(3),
(6),
(3),
(1). We note that
.
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
and
. The coefficients of
are constant while the coefficients of
and
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
and
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
and
. We also note that
and
are symmetric but
is nonsymmetric.
(2) The damping matrices
(due to
) are also composed of addition of
and
. The coefficients of
are constant, while those of
and
are linear and quadratic functions of the components of displacement gradient tensor. Matrices
and
are symmetric but
is nonsymmetric.
(3) It is well understood that constant coefficient damping
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
and
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
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
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).
(30)
A finite element formulation based on integral form of (30), integrated over a discretization
of
using fundamental lemma of the calculus of variations followed by GM/WF yields (also see ref [6] ):
(31)
The following details can be found in ref [6] :
Displacement gradient tensor
is given by
(32)
in which
(33)
Let
(34)
Using local approximations for
, we can write:
(35)
We decompose Green’s strain tensor
into linear
and non-linear components
(36)
or
(37)
(38)
Using (35) in (36)
(39)
Thus,
(40)
or
(41)
where
(42)
(43)
In (31),
is a vector of secondary variables and
is the load vector (due to body forces in case of (30)). We assume that
is conservative and
are nodal degrees of freedom for an element e.
Recall
(44)
(45)
Using (43) and (45) in (31) we can write
(46)
Expanding the integrand in (46), we can write
(47)
where
(48)
(49)
(50)
are the element stiffness matrices.
is a constant coefficient matrix and the coefficients of
and
are linear and quadratic functions of the coefficients of
, hence linear and quadratic in dofs
.
and
are symmetric but
is non-symmetric. As expected, these are exactly same as in the case of corresponding IVP [6] .
We can write (46) or (47) as
(51)
Equations (51) represent a system of nonlinear algebraic equations. In order to find a solution
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
be the initial solution (assumed or guess), then
(52)
Let
be a change in
such that
(53)
We expand
in (53) in Taylor series about
and retain only upto linear terms in
.
(54)
From (54), we can determine
(55)
Improved solution
is obtained using
(56)
If
(57)
then the solution
is considered the converged solution. Otherwise, we set
and repeat the iteration process from (55)~(57).
is a preset tolerance for the computed zero.
Surana et al. [6] have shown that variation of
leads to
(58)
Substituting for
in (58) from (43)
(59)
Let
(60)
(61)
(62)
(63)
and
are all symmetric.
is called the tangent stiffness matrix for element “e”. Assembly of the element equations follows the usual procedure.
is called tangent stiffness matrix for the discretization.
Remarks
(1) In the neighborhood of static bifurcation (load at which amplitude of motion experiences jump),
is momentarily singular causing instability. At bifurcation load,
). That is, at static bifurcation point, compressive or negative
is exactly equal to the sum of
,
and
. Thus, significance of
is now rather obvious.
(2) We note that in calculating the increment solution of
using (55) and (59), there are four stiffness matrices involved.
and
have constant coefficients, linear and quadratic coefficients of the displacement gradient tensor (or
, hence
). While
is same as in (47) but
and
are not the same as
and
in (49) and (50).
(3) We also note that in
and
, the dependence of the coefficients on the gradients of displacement (i.e., the degree) is same as in
and
.
(4) The matrix
(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
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
, causing decrease in stiffness of the deforming volume of matter.
(5) In static and dynamic bifurcation physics,
plays a crucial role in the mechanism of sudden energy release that leads to bifurcation. Without
, bifurcation cannot be observed in the nonlinear dynamics. Additionally, from the derivation of
, we observe that
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
also controls dynamic bifurcation physics in nonlinear dynamics. Static bifurcation cannot exist without
, 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
and
, 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
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
, which enhances the total stiffness of the medium. On the other hand, a compressive stress field results in negative
which reduces the total stiffness of the medium. Thus
plays a major and significant role in the static as well as dynamic bifurcation physics.
(4) As discussed earlier, translational inertial physics due to
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
and the change in stiffness due to
all play a role in dynamic response. We note that negative
and
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
is always positive, dynamic bifurcation will require the minimum possible damping (to reduce stiffness) and enough decrease in stiffness due to
. On the other hand, when
is negative,
along with some decrease in stiffness due to
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,
plays a crucial role in bifurcation phenomenon. In the absence of
(corresponding to infinitesimal deformation physics), bifurcation is not possible. Presence of
requires finite deformation physics, which is essential for the occurrence of static or dynamic bifurcation. This highlights the significance of the role played by
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 (
).
(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
(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
. Presence of positive
and presence of negative
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).
(64)
In (64), m is mass,
are damping coefficients associated with constant, linear and quadratic damping in u. Likewise,
and
are similar stiffness coefficients.
Comparing (64) and (1), we note the following:
(a) Damping coefficients
and stiffness coefficients
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
which remain the same for constant, linear, and quadratic damping in displacement gradients. Same is true for stiffness. That is, in (64)
and
are three coefficients associated with constant, linear and quadratic stiffnesses in u, where as in model (1), only two material coefficients
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
and
, respectively (due to one dimensional nature of the model). But in the case of (1), two material coefficients,
and
, are needed for linear stiffness as well as two material coefficients,
and
, 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),
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
, 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
, 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
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
of spatial domain
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
and
. Using
, BLM can be written as:
(65)
in which
(66)
Here,
is the equation of state for compressible matter (thermodynamic pressure).
(67)
We can substitute (66) and (67) in (65). The resulting equation constitute a system of PDEs in displacements
. The mathematical model (65)-(67) consists of three balance of linear momenta equations and six equations from the constitutive theory for
), resulting in a total of nine equations in
and
, nine dependent variables. Equations (65) and (67) hold for
the space time domain of the IVP defined by (65) and (67). We divide space-time
domain into space-time strips
.
is a typical ith space-time strip
for
;
. Consider a discretization
of nth space-time strip
into p-version hierarchical space-time finite elements with higher order global differentiability local approximations.
(68)
is the space-time element “e”. We calculate converged solution for a space-time strip and time march to obtain entire evolution.
Let
be local approximation of
,
(equal order, equal degree) over
(see reference [23] for details). Then
is approximation
of
over
. By substituting
in the nine PDEs (equations (65) and (67)), we obtain nine residual equations
for
in which
are the residual functions from the PDEs. Let
be the degrees
of freedom in the local approximation
, then
are the degrees
of freedom for
approximation of
over
. The space-time residual functional I over
is constructed using
(69)
in which
are the residual functions obtained from (65) and (67) when
are substituted in them. Functional I(.) is a nonlinear function of
. Following reference [32] , we proceed as follow. If I(.) is differentiable in
, then
is unique and
is a necessary condition for an extremum (minimum in this case) of I(.) i.e., the following holds.
(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
yielding extremum of (70) is also a solution of Euler’s equation i.e., PDEs (65) and (67). Since PDEs are nonlinear,
in (70) is a nonlinear function of
. We use Newton’s linear method with line search to obtain solution
from (70).
Let
be an assumed (or guess) solution, then
(71)
Let
be incremental change in
such that
(72)
We expand
in (72) using a Taylor series expansion about
and retain only upto linear terms in
:
(73)
From (73), we can solve for
:
(74)
Following reference [6] [32]
is obtained using
(75)
Improved solution is given by (using line search [32] )
(76)
If
(77)
then
is the converged solution. Otherwise, we set
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
, hence we can write (1) as (using
for
)
(78)
following reference [6] ,
and
defined by:
(79)
(80)
is the usual consistent mass matrix (defined in ref [6] with constant coefficients).
The coefficients of
and
are up to quadratic functions of the gradients of displacements
.
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
(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
and not at
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
and
as displacements, velocities and accelerations (since
are degrees of freedom associated with displacements and
are first and second time derivatives of
). 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
,
and
be known solutions at time t and
(or
),
(or
) and
(or
) be unknown solution at time
. Let the acceleration be linear between t to
. If
is time measured from time t, then based on linear acceleration assumption in interval
, we can write
(81)
Integrating (81) with respect to (w.r.t)
and using
to determine constant of integration, we can obtain the following for velocity at
.
(82)
Integrating (82) w.r.t
and using
to determine constant of integration, we can obtain the following for displacement
i.e., we consider the following.
(83)
In Newmark linear acceleration method, we satisfy (78) at
i.e.,
(84)
In (84), the coefficients of
and
are up to quadratic functions of
. In (82) and (83), we substitute
to obtain
(85)
(86)
From (86), we determine
and then substitute it into (85) to determine
.
(87)
and
(88)
Substituting (87) and (88) in (84) and regrouping terms gives
(89)
Equation (89) is a system of nonlinear algebraic equations in
. We use Newton’s linear method to obtain solution
of (89). Details are given in the following.
We rewrite (89) as follows:
(90)
Let
be initial or known guess of
at
(commencement of evolution), then
(91)
Let
be incremental change in
such that
(92)
We expand
in (92) using Taylor series expansion about
and retain only up to linear terms in
(93)
From (93), we can solve for
(94)
We refer to the coefficient matrix
as tangent matrix
, equation (105) that controls determination of incremental solution. Improved solution
is obtained using
(95)
We check for convergence of the Newton’s linear method.
If
(96)
then
is the converged solution. Otherwise, we set
and repeat (94)-(96) till (96) is satisfied. Δ in (96) is a preset tolerance for the computed zero.
We still need
in (94) to compute
. This is obtained using (90). First, we note
(97)
(98)
(99)
in which
(100)
and
(101)
Following [24] , we can derive the following
(102)
(103)
and
are tangent damping and stiffness matrices for spatial discretization.
is given by
(104)
is obtained by replacing
in (104) by
(defined in (100)).
By substituting (102) and (103) into (98) and (99), and then using these in
, we obtain the following expression in which
:
(105)
The matrix
is called tangent matrix.
We still need
.
(106)
(107)
Let
Then
and
(108)
where
(109)
Substituting (109) in (105), we have
, hence
can be calculated using (94).
Remarks
Since, now we know the dependence of incremental solution on ‘tangent matrix’
, we can discuss various physics in
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
defined in (105).
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
.
(2) Now we can conclusively determine the influence of various physics through
on nonlinear dynamic response and the dynamic bifurcation phenomenon.
(3) Nonlinear dynamic response depends upon symmetric mass matrix
, tangent damping matrix
and tangent stiffness matrix
(as well as the last term in (105) containing variation of
). We recall that
(110)
Coefficients of
are constant and those of
and
are linear and quadratic functions of dofs
.
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
,
and
are all zero. Negative
(due to compressive field) and positive
(due to tensile stress field) are likely to result in different dynamic response, hence different dynamic bifurcations.
(4) Tangent damping matrix
has identical structure as
only the material coefficients are different.
(111)
Matrix
is similar to
. It uses
as opposed to
. Thus, in calculating
,
and
in
are replaced by
,
corresponding to strain rate
. Thus, constant, linear and quadratic damping as well as
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
, 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
(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:
(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
(left end) and subjected to harmonic excitation
at
(right end), as shown in Figure 1(a). If all points in each cross section of the rod deform in
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
and the constitutive theory for deviatoric second Piola-Kirchhoff stress
are the only two equations that constitute the mathematical model, assuming equilibrium stress
to be zero. We can write the following for BLM and the constitutive theory (in Lagrangian description using Equations (19)-(29)):
(113)
(114)
(115)
(116)
We substitute (115) and (116) into (114), then
from (114) into (113) to obtain a single PDE in
(letting
)
(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
be spatial discretization of
(Figure 1(c)).
Let
be local approximation of
over
. Then
is the approximation of
over
. The integral
form of (117) over
based on fundamental lemma of the calculus of variations using test function v is given by:
(118)
or
(119)
We consider scalar product
, where
(120)
is approximation functions matrix and
are time dependent nodal degrees of freedom for an element e. The total degrees of freedom
for
are given by:
(121)
Using (121) in the integral form for an element “e” with domain
, and following reference [23] , i.e., using GM/WF in space for
followed by assembly of element equations based on (119), we obtain the following (similar to (1))
(122)
where [24]
(123)
(124)
and
have constant coefficients where as coefficients of
and
are linear and quadratic functions of displacement gradients
(
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),
are purely hypothetical or phenomenological, whereas in (122),
used in
and
used in
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
of
that we obtain a system of nonlinear ODEs in that contain degrees of freedom
and
corresponding to the spatial discretization
. 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
defined by (105) i.e., inverse of the following matrix for calculating incremental change (or improvement) in the assumed solution:
(125)
Matrices
and
are derived considering
, which corresponds to
and
in (122).
and
consist of
(126)
(127)
and
are all symmetric.
and
have constant coefficients but coefficients of
and
are linear and quadratic in displacement gradients.
describes the influence of stress field on the stiffness.
likewise has similar structure as that of
(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.
and
are implicitly contained in (123) and (124) but can be explicitly observed in (126) and (127). Thus, in the following, we can view
in (122) is due to constant, linear and quadratic functions of displacement gradients as well as due to
, the influence of stress field on stiffness. Similar view holds for damping matrices
. Compressive stress field results in negative
, causing a reduction in stiffness. On the other hand, tensile stress field resulting in positive
will cause an increase in stiffness. Explicit appearance of
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
is not possible. However, the presence of
in stiffness and damping as well as extraction of
and
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
and
must be associated with positive or negative values of
and
. At this stage, it is intuitively plausible to assume that a negative
will result in reduction of stiffness (similar to negative
in model problem II) and positive
will result in increase of stiffness (similar to positive
) 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
. In this model problem, there is no mechanism of positive
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
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
. Considering zero body forces and
(for simplicity):
(128)
(129)
(130)
Equations (128)-(130) are three nonlinear PDEs in
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
, the discretization of nth
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
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
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:
(131)
in which F is the applied force.
Model problem II for BVP can be written (in the absence of body forces) as
(132)
(133)
(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
for all studies. In the first study, we choose
. For this choice, we have tension in the spring resulting in positive
(when compared with axial rod). In the second case, we choose
and
, simulating a negative
(in analogy to axial rod).
Force F versus u graphs are shown in Figure 2. When
is negative, the relationship between F and u shows progressively increasing displacement u with
increasing F. Additionally, the derivative
also increases as F increases for all
.
When
is positive, increasing F results in progressively increasing positive
, resulting in progressively decreasing displacement u. When
,
increases as F increases for all
. In the vicinity of
(between
and
), u increases more sharply with increasing F indicating progressive loss of stiffness due to negative
. Between
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.
(135)
in (135) is infinity (condition for bifurcation) when
(136)
(a)(b)
Figure 2. (a) Model problem I: Force F versus u1 plot for
and
; (b) Model problem II:
versus u1.
Using
and
we obtain solutions
and
. Clearly,
is not admissible, hence
becomes infinite at
, 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
and
.
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 u1 and
as dependent variables. A dimensionless rod of one unit is considered, fixed at left end (
) and subjected to
at
. Solutions are calculated for
. Graph of
versus u1 is shown in Figure 2(b). This graph is similar to a scaled graph of force F versus u1.
For
, the positive
contributes to progressive increase in the stiffness of the rod as
increases. This results in progressively reduced displacement u1 for the same increment of
during loading.
On the other hand, when
is negative, the negative
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
during loading.
In the vicinity of
, there is total loss of stiffness. As a result, computations fail beyond
. We note that when loading is compressive, there is total loss of stiffness in the rod for some value of
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
is clearly observed in model problem II. In case of model problem I, there is no
, but
is assumed to simulate
.
(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
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 u1) 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
at the free end of the rod, the rod experiences a steady compression. As a consequence,
is negative. Thus in model problem II, we can only study negative
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
(model problem I) or
(model problem II) and dissipation or dissipation coefficient c, appropriate combination of
(or
) 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
(or
) and higher values of c with low values of
(
) 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.
(137)
Choice of
and
are discussed in the following. We consider
in increments of 0.1. For a circular frequency
,
is
chosen as
i.e.,
is
of the time period of
. Computation of
(a)(b)(c)
Figure 3. (a) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
; (b) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
; (c) Model problem I: Frequency ω versus maximum amplitude u1 response for
.
evolution is commenced at
and continued till the response becomes cyclic (not necessarily periodic). At this point, maximum value of u is recorded for
. Frequency
is changed to
maintaining the final state of deformation at
as initial state at
. Computations are continued till the response is cyclic for
. At this point, maximum value of u is recorded for
. This process is continued for progressively increasing frequencies till we reach
. This frequency sweep from
is referred to as
(or
). To determine any path dependencies in the frequency response, the computations are started at
with progressively decreasing frequencies in the same manner as for
. We refer to these computations with progressively decreasing frequencies as right to left
(or
to
,
).
Choice of
and
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
and
. Thus, in most studies we adhere to the unit values for these parameters. In model problem II, the material coefficients in
and
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
and
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
,
. However, there is one exception: if we want to simulate negative
and capture the physics similar to compression in the axial rod, we need to choose a negative value for
.
8.2.1. Case (a): Simulations for Negative
i.e.,
.
We choose
but
to simulate negative
(as in case of axial rod). In this section, we present two studies. In the first study, we choose a fixed
and progressively reduce
till computations are no longer possible due to singular coefficient matrix. In the second study, we choose fixed values of
and vary
till computations cease due to singular coefficient matrix in the computations. These studies are aimed at determining if bifurcation is possible when
with the choice of coefficients
and
considered here.
In the first study, we choose
and perform computations for progressively decreasing
. The evolutions are computed for progressively increasing frequencies
(
) as well as for progressively decreasing frequencies
(or
) of excitation in increment
of 0.1. We choose integration time step
. 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
and
for
. Maximum positive peaks as well as absolute value of maximum negative peaks are shown in the graphs. In the second study we choose
and progressively increase
and compute frequency response for
(
) and
(
) using
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
and
and 0.035.
Discussion of results
(1) From Figures 3(a)-(c) we note that progressively reduced values of
for fixed
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
. Path dependency increases with reducing values of
. For
and
we observe no bifurcation in Figures 3(a)-(c). We remark that for
, choice of
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
. 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 u1 response for
,
; (b) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
; (c) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
.
(2) From Figures 4(a)-(c), we note that for fixed
, progressively increasing values of
results in progressively increasing maximum amplitude of motion, as expected. For
(Figure 4(a)) positive and negative peaks for each path show some differences. Positive peaks and the negative peaks for
and
are relatively close together. At
these differences increase (Figure 4(b)). We also note existence of mild bifurcation (more rapid changes) that is more pronounced in both peaks for
path.
path also shows mild bifurcation but at a different frequency than for
path. Path dependency can be clearly observed. For
(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.,
and
paths yield same frequency response. We also point out that in this case also for
computations fail if
is increased beyond 0.035.
(3) Studies in (1) and (2) show that when
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
(both needed for dynamic bifurcation) computations fail.
(4) We recall that static bifurcation for this model problem I, with
and
occurs when
. 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
and
, hence the reason for choosing the values given in the studies presented here.
(6) We note that for negative
, bifurcation occurs to the left of the peak in all studies.
8.2.2. Case (b): Simulations for Positive
i.e.,
: Model Problem I
We consider solutions of model problem I for positive
i.e.,
. In this case, the influence of the stress field is to increase total stiffness. We choose
in all studies presented in this section and consider different choices of
and
. 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
. Furthermore, in case of axial rod, model problem II, positive
cannot be simulated due to
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
and
. For each value of
, frequency response is obtained for the frequencies
(
) and for
(
) in increments of 0.1. For each
, at each frequency, positive and negative peaks are recorded
when the motion is steady and cyclic. Integration time step
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
and 0.01.
(1) Figure 5(a) shows frequency response for
. For
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
path, we also observe dynamic bifurcation (not as distinct as in case of
) but at a different frequency than for
. Path dependency can be clearly observed in all three figures.
(2) Frequency response for
and
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 f0 can be eliminated by lowering f0 values, thus reducing presence of nonlinearity and stiffness. Likewise, for a fixed dissipation physics progressively increasing f0 could likely result in dynamic bifurcation phenomenon.
(a)(b)(c)
Figure 5. (a) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
; (b) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
; (c) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
.
8.2.3. Case (c)
In this second study, we keep
and
fixed and choose
and 0.1 i.e., progressively reducing dissipation. Frequency response for all three values of damping for
and
frequency range are shown in Figures 6(a)-(c). From Figure 6(a) and Figure 6(b), we note that for
and 0.5, negative and positive peaks for
and
are almost identical but we clearly observe path dependency i.e.,
and
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
. First, we note that maximum amplitude in this case is roughly four times that of for
. We clearly observe bifurcation for both
and
loading paths. Negative peaks are always higher than positive (as explained earlier) and both peaks show presence of bifurcation in case of
bifurcation occurs at a higher frequency compared to
. Both positive and negative peaks at bifurcation in case of
are higher than the corresponding peaks for
. As expected, since the positive
implies positive
, the bifurcation occurs to the right of the peak amplitude.
In this study, we clearly see the role of
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
due to stress field.
(a)(b)(c)
Figure 6. (a) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
; (b) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
; (c) Model problem I: Frequency ω versus maximum amplitude u1 response for
,
.
8.3. Model Problem II: Numerical Studies
We consider dimensionless form of model II (Equations (128)-(130)). First, we write (128)-(130)) using hat (
) on all quantities indicating that they have their usual dimensions or units.
(138)
(139)
(140)
We introduce the following reference quantities (with subscript zero or “ref”) and dimensionless quantities:
(141)
Using (141) in (138)-(140), we obtain the dimensionless form of the equations (138)-(140).
(142)
(143)
(144)
in which
if we choose
and
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
;
being final value of time. A discretization
of the space-time domain
into space-time strips
is shown in Figure 7(c).
Figure 7(d) shows a discretization
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
,
at
, ICs:
,
at
. BCs at
consists of
. Using
values of the degrees of the freedom are calculated for the element nodes at
and for
. For the second space-time strip, these dofs are calculated for
and so on. This process enables application of
at
in time accurately. In this process, a change in ω is done at a value of time t when
, so that change in ω does not cause unnecessary discontinuity in the application of
at
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 (
, order of the approximation space) are sufficient to obtain space-time residual functional I for
of
or lower indicating accurate evolution for
.
Evolution is continued using time marching till steady cyclic response is reached
for each frequency. We use integration time step of
. 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
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
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
and
(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.,
and
frequency responses are the same.
Figure 8(b) shows frequency response for
and damping values of
and 0.0035. For
, we observe lower peak values. Frequency response for
and
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
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
In this study we choose a fixed value of
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
and
.
Figures 9(a)-(d) show frequency response for
and 0.002. For all four values of damping, (
), (
) frequency response are shown using positive and negative peaks. In Figures 9(a)-(c), we observe that for each value of
that
and
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
in Figure 9(d) shows distinct bifurcation to the left of the peak (due to negative
). 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
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 (
and
).
For this choice of
, we could observe bifurcation before the computations cease. In Figure 9(d) we observe that
and
frequency response show mild path dependency. This phenomenon of mild path dependency has also been reported in ref [26] .
8.3.3. Case (iii): Progressively Increasing
for a Fixed Damping c
In this study, we choose a fixed value of damping coefficient
and determine frequency response for progressively increasing
, 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
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
frequency response for
, 0.035, 0.045, 0.055. In Figure 10(a) and Figure 10(b) for
and 0.030,
and
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
, 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.,
and
frequency response are almost same. In Figure 10(d) for
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
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
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
Is Positive
In the model problem studies for axial rod (model problem II) subjected to harmonic excitation,
due to stress field is always negative resulting in reduction of total stiffness. In the type of applications, where
is negative, the inertial physics due to linear acceleration and
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
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
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
, just like in case of axial rod where
is always negative for excitation at the end of the rod, in such applications
is always positive due to midplane always being in tension during dynamic motion. Positive
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
studies and studies for positive
in case of model problem I, we expect that when
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
in BLM and only consider isothermal physics, then we obtain three PDEs in displacements
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
and
, 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
, stiffness
and damping
matrices. Only
is symmetric.
and
are nonsymmetric and are often referred to as secant stiffness and damping matrices. We have shown that
, the stiffness matrix for element e, can be decomposed into
and
(Equations (48)-(50)),
with constant coefficients and
with coefficients that are linear and quadratic functions of
, dofs for an element e. The same decomposition can be applied to each of
. At this stage, we have little more information, and nonlinear dynamics now depends upon stiffness
and damping matrices
whose coefficients are constant, linear and quadratic functions of dofs and of course on constant coefficient mass matrix
. 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
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
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
, the tangent matrix defined by (105). Thus, the dynamic solution depends on
that consists of symmetric
and
.
in terms of sum of
and
and likewise
is the sum of
and
.
have constant coefficients. Coefficients of
are linear functions of
and those of
and
are quadratic functions of
.
is the stiffness matrix due to presence of stress field. We note that
and
are not the same as
resulting from decoupling of space and time in which
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
and thus on the new
and
matrices. Concept of
in this approach is instrumental in understanding bifurcation physics. Going back to (64), linear and quadratic stiffness terms maybe viewed as synonymous to
and
in (1), but they fail to explain the physics due to
. It is worth noting that
is also realized in finite deformation, finite strain BVPs. Confirming that
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
, hence reduces total stiffness. On the other hand, tensile stress field will result in positive
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
,
(all symmetric) contribute to total stiffness. Influence of the degree of nonlinearity is reflected in
and
. For large finite strain and large finite deformation,
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
due to the stress field plays critical role in dynamic bifurcation, negative
enhances likelihood of the existence of dynamic bifurcation whereas positive
helps in inhibiting the likelihood of dynamic bifurcation due to increased stiffness. The last factor influencing dynamic bifurcation is inertial physics due to accelerations (
) 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
and on the other hand inertial physics due to acceleration and negative
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
is positive, it enhances stiffness, hence in this case damping, stiffnesses and
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
and negative
must be different. We have demonstrated this in the model problem studies that this indeed is the case. In case of negative
due to compressive stress field (model problem II, also model problem I with
) dynamic bifurcation is always to the left of the peak. When
is positive, the dynamic bifurcation is to the right of the peak (shown for model problem I with
).
(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
and
in (64) do not correspond to negative
and positive
. Derivation of
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
and their time derivatives
that correspond to the nonlinear problem. These cannot be transformed using modal basis determined from the linear form of ODEs (1) in which
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
becomes equal to the stiffness (
) 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
when
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
.
(ii) In static bifurcation, negative
is essential (compression). Positive
can never result in static bifurcation. Since inertial physics due to acceleration is absent in static bifurcation, it is a balance between
and the rest of the stiffnesses.
(iii) Dynamic bifurcation can occur with positive and negative
as in this case in addition to
the inertial physics due to acceleration (lowering stiffness) is present.
(iv) When
is negative, dynamic bifurcation occurs to the left of the peak and when
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
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
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.