Transient analysis of trusses considering nonlinear elastic and viscoelastic material models

The use of simple bar elements in nonlinear structural finite element formulations has the academic advantage of uncoupling element technology issues from the structural phenomena to be observed. In this work, we present a finite element setting for the formulation of different nonlinear material models applied to the transient analysis of trusses. While nonlinear elasticity is considered by studying a Hooke-like linear relationship between different pairs of nonlinear measures of stress and strain, hyperelasticity is formulated using Ogden’s model. Viscoelasticity is introduced using a generalized Kelvin rheological model to account for strain rate effects. The finite kinematics is set in a corotational total Lagrangian description where the virtual work is described using the Second Piola-Kirchhoff and the Green Lagrange measures. Although the derivation is omitted, the consistent tangent moduli are given for all these cases. Numerical problems involving simultaneously different truss models are studied and made available as benchmarks since little comparative data is found in literature.


INTRODUCTION
Due to their relative simplicity, truss elements are frequently employed as a first step to study nonlinear formulations.Some examples are the works of Driemeier et al. (2005), Greco et al. (2006), Carniel et al. (2015), Pascon (2016), Rabelo et al. (2018), Guggenberger and Krenn (2019), Ananthapadmanabhan and Saravanan (2020), Endo et al. (2021) and Rezaiee-Pajand et al. (2023).In this work we describe the kinematics of a bar element using concepts of continuum mechanics and develop a finite element framework to study different material models in nonlinear elastic and viscoelastic transient analyses.Particularly, we extend the Saint-Venant-Kirchhoff material model, which relates linearly the Second Piola-Kirchhoff stress to the Green-Lagrange strain tensors, to consider linear relationships between other stress-strain measures.This behavior is possible provided the strains are finite, but relatively small.We also introduce the incompressible Ogden hyperelastic material law, which admits larger strains.Then, we use a generalized Kelvin rheological model to consider rate dependent effects in each of the nonlinear Hooke-like elastic material models.In the hyperelastic case, structural damping is introduced using a Rayleigh damping matrix.The formulation is set in a corotational total Lagrangian description where the virtual work is described using the Second Piola-Kirchhoff and the Green Lagrange measures.One of the contributions of this work is the presentation of the consistent tangent moduli for all these cases.Although the derivations are omitted, the expressions have been carefully validated by verification of quadratic convergence tests.
The bases for the formulation framework developed in this article are summarized as follows: the nonlinear analysis involving different constitutive models is based on Muñoz-Rojas (2023); the considerations involving Rayleigh damping are based on the discussions found in Charney (2008), Jehel et al. (2014) and Chen et al. (2015); the time discretization follows the principles of the Newmark family methods (Subbaraj and Dokainish, 1989); the hyperelastic material law is the Ogden model (Ogden, 1972) which is particularized to truss analyses as, for example, in Fonseca and Gonçalves (2022); as for the viscoelasticity algorithm, after studying several integration methods, such as those presented by Evangelista (2006), Araújo et al. (2010), Keramat and Ahmadi (2012), Kühl et al. (2017).Chazal and Moutou Pitti (2009), Chazal and Pitti (2011) and Ananthapadmanabhan and Saravanan (2020), the algorithm presented by Carniel et al. (2015) is adapted to an explicit approach.
Although there are many publications describing the limitations of linear relationships between nonlinear measures of stress and strain, particularly for bars (Crisfield, 1991), the corresponding discussion about a linear relationship between nonlinear measures of stress and strain rates is not so frequent.In this work we explore both issues with examples involving 2D and 3D trusses made of bars with different combinations of material models, all of them subjected to suddenly applied loads and undergoing snap-through instabilities.Under such conditions, the consideration of inertia forces is mandatory, implying the calculation of velocities and accelerations, which affect dissipation.The material parameters adopted in the numerical examples follow those used by Ogden (1972), Abdelrahman and El-Shafei (2021) and Endo et al. (2021).

BAR KINEMATICS
Consider a bar located in the three-dimensional space subjected to a uniform uniaxial stretch.The bar length, area and volume are denoted by  0 ,  0 and  0 in the undeformed configuration, and ,  and  in the current configuration.We define a global system of reference with origin  and base vectors { 1 ,  2 ,  3 }, and a local corotational system of reference with origin at the beginning of the bar and base vectors � 1  ,    ,  3  �.Following classical notation, we use capital letters to refer to coordinates in the original (undeformed) configuration and small letters to refer to coordinates in the current (deformed) configuration.With this convention, in the undeformed configuration, each material point of the bar has global and corotational local coordinate vectors  = (, , ) and   = (  ,   ,   ).Accordingly, in the current configuration, the coordinates of the same material point are given by the vectors  = (, , ) and   = (  ,   ,   ), as shown in Figure 1.
Under uniform uniaxial stretch, the location of a material point in the current configuration can be expressed with respect to the corotational local system of reference as   = (  ,   ,   ) = (  +   ,   +   ,   +   ) (1) where (  ,   ,   ) corresponds to the displacement suffered by the point.
Assuming that the bar behaves the same way in any direction of the cross-sectional plane, we have (Holzapfel, 2002;Dunne and Petrinic, 2005) Latin American Journal of Solids and Structures, 2024, 21(1), e521 3/31 where and Using the logarithmic strain measure to define the Poisson ratio, according to it turns out that Hence, the relation between the updated and initial cross section areas is and the strain gradient with respect to the corotational local system of reference is given by whose Jacobian is From the polar decomposition theorem, it turns out that the rotated right Cauchy strain tensor is Latin American Journal of Solids and Structures, 2024, 21(1), e521 4/31 where the elongation tensor   is used to define the Seth-Hill family of strain tensors Due to the hypothesis of uniaxial uniform stress state assumed for the bar, the only strain component that produces work and must be considered is the axial one.The axial component of the Seth-Hill family of strain tensors can be expressed as with  defining different measures.Some measures of strain are • Green-Lagrange strain ( = 2): • Rotated engineering (Biot) strain ( = 1): • Logarithmic strain ( = 0):

PRINCIPLE OF VIRTUAL WORK (PVW)
The virtual work principle establishes that, for every instant of time and every kinematically admissible virtual displacement field, a body is in equilibrium if the sum of the internal virtual work performed by internal forces and the inertial work equals the external virtual work performed by the external forces (Belytschko et al., 2014).Hence, for each instant of time, where   is the virtual work performed by the external load,   is the virtual work performed by internal forces and   is the virtual work performed by inertia forces.
In order to account for structural dissipation, normally associated to friction among assembled parts, one can introduce an additional term in Eq. ( 15), which then reads where   is the virtual work performed by structural damping forces.As depicted in Figure 2, for each bar element the internal work expression is given by with   being the axial internal force along the element length.Within each bar, the axial stress can be defined by alternative measures, for instance • the Cauchy stress: • the rotated engineering or First Piola-Kirchhoff stress: • the Second Piola-Kirchhoff stress: These stress measures can be confirmed by the application of the classical continuum mechanics expressions and where  and  are the First and Second Piola-Kirchhoff tensors respectively.Although Eq. ( 17) can be expressed using any pair of energetically conjugate stress-strain tensors, the derivations developed in this article employ the Green-Lagrange strain tensor and the Second Piola-Kirchhoff stress tensor.Hence, we have (Crisfield, 1991) Latin American Journal of Solids and Structures, 2024, 21(1), e521 6/31 The damping and inertia virtual work parcels mentioned in Eq. ( 16) are respectively and where  and  are the specific damping and mass coefficients.

SPACE DISCRETIZATION OF THE PRINCIPLE VIRTUAL WORK
For each bar we define the interpolation matrix based on a natural coordinate   [−1, 1] along the element axis.This matrix is used to approximate coordinates, displacements, velocities and accelerations by where the asterisk stands for continuum expressions being discretized and  �  ,  �  ,  �  ,  � ̇ and  � ̈ are the nodal values of the corresponding element vectors.For a linear two-noded bar, the interpolation functions are 2 () = 1 2 (1 + ). (33)

Discretization of the virtual internal work
Replacing Eq. ( 12) and ( 13) into (23), the virtual internal work can be discretized by implying with .
The transformation of the internal force vector to the global system of reference is performed by where  is the rotation matrix, composed by the direction cosines (Bathe, 2006).

Discretization of the inertia virtual work
Replacing Eq. ( 29) and (31) in Eq. ( 25), we obtain where   is the consistent element mass matrix and  () is the specific mass of the material of the element.Further, where The global mass matrix is assembled as usual, where  is the total number of elements and Λ is the assembly operator.

Discretization of the damping virtual work
Introducing Eq. ( 29) and (30) in Eq. ( 24), it results that where   is the consistent element damping matrix.Moreover, where As for the mass, the global damping matrix is obtained applying the assembly operator over all the elements, where   =      is the element damping matrix   rotated to the global system of reference.

Spatially discretized equilibrium equations
To trace the equilibrium path along the load application, the total load is divided into a finite number of time steps .As the PVW must be valid for an arbitrary virtual kinematically admissible displacement, for an external force in the corotational local system of reference  +1 we obtain that the discretized version of the equilibrium equations is which vanishes at equilibrium and takes non-zero values when an out-of-balance force occurs.

LINEARIZATION OF THE EQUILIBRIUM EQUATIONS
Owing to the term  +1 ( � +1 ), the residual is a nonlinear function of the global displacement vector, so that the solution of Eq. ( 47) requires the linearization where the derivative of internal force vector  +1 with respect to the global displacement vector defines the tangent matrix   and Δ � +1 ) .

Tangent matrix
Dropping the  + 1 index for neatness, we can express the tangent matrix as where    is the element tangent matrix with respect to the global system of reference.Employing the chain rule of differentiation on Eq. ( 49), this matrix can be written as the sum of two other matrices,   1 and   2 . where ( �  ) and  are given by (Muñoz-Rojas, 2023) is the identity matrix and  � is the vector of nodal current coordinates with respect to the global system of reference.
On the other hand, with where (Muñoz-Rojas, 2023) and   is the tangent modulus where ( 2 ,   ) is the conjugate pair adopted in this work.
Notice that the differentiation of  2 causes a dependence on the material model adopted.This important issue will be further discussed in Sec. 7.

Simplifications on the damping matrices
It is usual to replace the consistent damping matrix by the so-called Rayleigh matrix, which is defined as a global damping matrix obtained by the linear combination of global mass and stiffness matrices.For this, the literature displays three alternatives (Charney, 2008;Jehel et al., 2014), 1. constant coefficients  and  in the linear combination between mass and the initial tangent stiffness (equal to the linear stiffness matrix) 2. constant coefficients  and  in the linear combination between mass and the updated tangent stiffness 3. updated coefficients   and   in the linear combination between mass and the updated tangent stiffness In Eq. ( 58), ( 59) and ( 60),  and  are called Rayleigh constants and the lower index  indicates that the values are updated at each load step and iteration.
The application of Rayleigh matrix is mainly for numerical convenience but allows to consider, in a simplified way, both, material dissipation (viscous materials) and structural friction at joints and supports.However, it is frequent to disregard dependence on  or   adopting  = 0 (Chen et al., 2015).Moreover, viscous materials are treated in a more adequate way by employing appropriate material models, as those described in Section 7.3 .

TIME DISCRETIZATION
The Newmark family of time discretization methods is based on the expressions where   and   define the specific method adopted and Δ is the length of time steps.Isolating  � ̈+1 in Eq. ( 62), and substituting Eq. ( 62) and ( 63) in ( 61) we find that, for each iteration k, the effective stiffness matrix  �   and the effective external force vector  �  are defined as (Subbaraj and Dokainish, 1989) Latin American Journal of Solids and Structures, 2024, 21(1), e521 10/31 and the solution of the resulting global system via the Newton-Raphson method gives In Eq. ( 66),  �� �  � is the residual given by Eq. ( 48) with the effective external force in Eq. ( 65) in place of   .It is well known that, in the case of quasi-static problems, the use of the exact tangent matrix should result in quadratic convergence (close to the roots).On the other hand, Chang (2004) shows that when using the Newmark family algorithms to solve nonlinear problems involving inertia forces, quadratic convergence cannot be expected.
In this work we adopt the method of average accelerations, for which   = 1/4 and   = 1/2.

MATERIAL MODELS
In Section 5.1, the tangent modulus was defined by Eq. ( 57) as the derivative of the stress with respect to strain employed to describe the internal virtual work, which depends implicitly on the material behavior.It is important to remark that the constitutive behavior of a given material has nothing to do with the conjugate pair employed to describe the internal virtual work.Considering that we are adopting the conjugate pair ( 2 , ε  ) to express the internal virtual work, for any given stress-strain pair chosen to describe the material behavior, say ( * , ε * ), the tangent modulus will have the form (68)

Linear relationship between nonlinear stress and strain measures
Particularly, for a constitutive law that relates linearly the arbitrary stress and strain measures * σ and * ε , Table 1 shows the tangent moduli for constitutive laws that relate linearly different stress-strain pairs, where the well-known Saint-Venant-Kirchhoff model corresponds to σ 2 = ε  (Muñoz-Rojas, 2023).
Table 1.Tangent moduli for linear relationships between different stress-strain measures and internal virtual work evaluated using the conjugate pair ( 2 ,   ) .

Constitutive law
Noteworthy, the choice of a linear relationship between a given stress-strain pair implies in a nonlinear relationship when considering a different pair.Also, adoption of a linear relationship between stresses and strains normally produces an unphysical behavior for large compressive strains (finite stress for λ 1 = 0).The exception is σ = ε  , as shown in Figure 3. Notwithstanding, for moderate strains (represented by vertical lines in the range between  1 = 0.6 and  1 = 1.4 for illustrative purpose), such linear relationships can be used to model different material behaviors.Figure 4 shows the interesting fact that the linear relation σ  = ε  can model an auxetic Latin American Journal of Solids and Structures, 2024, 21(1), e521 11/31 material ( = −0.1)respecting the physical condition that the Cauchy stress σ → −∞ when λ 1 → 0. As expected, when small strain occurs, a linear relationship between any different stress-strain pairs leads practically to the same material behavior.

Hyperelasticity (Ogden's model)
When the work done by the stresses during a strain process depends only on the initial state at time  0 and the final configuration at time , the material is termed hyperelastic (Bonet and Wood, 2008).The hyperelastic model proposed by Ogden (Ogden, 1972) is defined in terms of the stretch tensor eigenvalues (principal stretches), λ 1 , λ 2 , λ 3 , where the elastic strain energy is given by and   and α  are material parameters.To obtain a physically consistent condition, the following constraints must be satisfied where  is the shear modulus in the reference (undeformed) configuration.It is usual to consider that hyperelastic materials are incompressible, in which case where   was defined in Eq. ( 8).
Introducing Eq. ( 6) into (74) we verify that for uniform uniaxial stretch, the incompressibility condition is met when  = 0.5, so that the strain energy is given by Moreover, the Second Piola-Kirchhoff stress is obtained differentiating the strain energy with respect to the Green-Lagrange strain, Figure 5 shows that this incompressible material law respects the physical constraint that the Cauchy stress  → −∞ when λ 1 → 0, hence allowing its use for large compressive strains.For comparison purpose, the relationships of σ  = ε  and σ 2 = ε  are also displayed for  = 0.5.The vertical lines at values  1 = 0.6 and  1 = 1.4 delimit a range of moderate stretches, where all material laws provide coherent, although different values.
Equation ( 76) is replaced in the expressions ( 35) and ( 51) of   and   1 .Differentiation of Eq. ( 76) provides the tangent modulus defined in Eq. ( 57) which must be introduced in the expression (54) for   2 .In linear viscoelasticity, for given constant constitutive parameters, the material response depends only on time, and it is possible to apply the Boltzmann superposition principle.This principle is given by the integrals and where () and () are the stress relaxation and the creep-compliance functions, respectively (Kühl et al, 2017).Equations ( 78) and ( 79) completely define the material response since the linear viscoelastic material properties are considered ideally constant.
Particularly, for creep-compliance response, a stress σ  is suddenly applied at time  0 and kept constant, while the stress is observed over time, that is and where ( −  0 ) and ( −  0 ) are the Heaviside and Dirac delta functions respectively, and σ  is the creep stress.Replacing Eq. ( 81) in ( 79) and applying the filter property, it turns out that or where frequently time  0 is taken as  0 = 0 (see Figure 6).Notice that the compliance is defined by Eq. ( 83) and in linear viscoelasticity this function does not depend on the creep stress σ  applied.A typical curve of linear creep-compliance (nonlinear in time) for any value of   is given in Figure 7.

Nonlinear Viscoelasticity
When we work with finite strains, the use of the pairs (  2 ,   ) and (   ,   ) in Eq. ( 82) implies a nonlinear relationship between the Cauchy stress and the logarithmic strain.For example, if the chosen constitutive measures are the second Piola-Kirchhoff stress and the Green-Lagrange strain, Eq. ( 82) becomes As it is possible to denote the Second Piola-Kirchhoff stress in terms of the Cauchy stress and the GL strain in terms of the logarithmic strain, it follows that which can be expressed as Similarly, if a linear viscoelastic behavior is defined between the engineering stress and strain, Eq. ( 82) in terms of the Cauchy stress and logarithmic strain is given by which can be rewritten as In both cases (Eq.( 87) and ( 88)) we can define the compliance as where   is the Cauchy creep stress.Figure 8 exemplifies the difference between the creep-compliance curves defined in terms of Cauchy stress and logarithmic strain when linearity is considered between other stress-strain pairs for given Cauchy stress values.

Generalized Kelvin-Voigt rheological model
The simple Kelvin-Voigt rheological model consists of a spring of elastic constant   arranged in parallel with a dashpot of viscosity constant   , as shown in the detail of Figure 9 (Ward and Sweeney, 2004).The generalized Kelvin-Voigt model consists in a series arrangement of a sole spring with a finite number of Kelvin-Voigt blocks, as shown in the same figure.For our purpose we consider that the springs obey a linear relationship between nonlinear stress and strain measures according to Eq.( 69).Damping in the dashpots is given by a linear relationship between the same measures for stress and strain rate.To avoid an excessive number of superscripts, in this Section we drop the asterisks that indicate the choice of a given stress-strain pair in Eq.( 69).
The stress-strain relationship for a single Kelvin block i is given by where   () and ̇  () are the viscoelastic strain and its time rate.Equation ( 90) can be rewritten as where is the relaxation time of the rheological block.and  is the number of rheological blocks involved in the material model.On the other hand, the stress of the sole spring is where the viscous strain corresponds to the sum of the viscous strains of each rheological block, so that If a stress is applied to the model, it acts on each rheological block and on the sole spring in the same way.Then Eq. ( 90) and ( 94) can be compared, resulting in where ω  =  0 /  .
In the case of viscous behavior, we define a linear relationship between stress and strain rate ( = ), as displayed in Figure 10.Notice that the strain rate depends on the stretch suffered by the material, yielding the expressions given in Eq. ( 97) The discussion made here on the implications of adopting different stress-strain rate pairs is relevant, since there are studies such as Douven et al. (1989) and Kaliske and Rothert (1997), that relate linearly pairs that are different than the Cauchy stress and the logarithmic strain rate.
Figure 11 shows the graphs of the auxiliary function ( 1 ) defined in Eq. ( 98) for the (,   ), (  ,   ) and ( 2 ,   ) material models.As well as when there is a linear relationship between stress and strain (Fig. 3), each of the constitutive laws represents a different material for moderate strains (for example,  1 between 0.6 and 1.4).In the case of large strains in compression, the relationships exhibit non-physical behavior.We would expect that, as the stretch approaches zero, the stresses would take on increasing values, as happens in the case of (,   ).However, for both constitutive laws (  ,   ) and ( 2 ,   ), when the stretch tends to 0, the values of ( 1 ) also tend to 0.
Latin American Journal of Solids and Structures, 2024, 21(1), e521  17/31 Since, according to Eq. ( 98) the stresses are obtained by multiplying ( 1 ) by the strain rates, it follows that the stresses tend to zero independent of  ̇1 values, contrary to the expected behavior.Figure 12 shows the behavior of ( 2 ,   ) and (  ,   ) for different values of .Although (  ,   ) presents physical sense for auxetic materials, the ( 2 ,   ) constitutive law does not give coherent results for any of the Poisson ratio values tested when large strains in compression are considered.

Numerical Integration of viscous strains
The numerical integration of viscous strains used in this work follows the explicit approach (Backward Euler Method), which was based on the adaptation of the three-dimensional formulation of Nedjar and Le Roy (2013) of implicit integration presented by Carniel et al. (2015).Other approaches are proposed by Zienkiewicz et al. (1968), Evangelista (2006) and Araújo et al (2010).
We assume that, for the numerical implementation, time is discretized into  steps of equal size Δ.For each step 0 <  < , the next step  + 1 will be equal to the time elapsed until the previous time step ,   , plus the time interval, i.e.,  +1 =   + Δ.Also, the total strain at step  + 1 is known and therefore will be constant throughout this time step.Then for time step  + 1, eq. ( 96) is such that Eq. ( 103) can be written as a linear system  +1  =  where Note that matrix  being diagonal implies that the cost of solving the system is zero.Despite this, there would be an operational cost due to the need to allocate space to store it, which is unnecessary.Therefore, although the equations are presented as a system, it is more efficient to solve them in a decoupled way, The consistent tangent modulus of a viscoelastic material (not to confuse with the tangent modulus defined in Eq.( 68)) in the time step  + 1 is, by definition, (remember that the asterisks have been dropped: in Eq. ( 68)) ( 107) where Using Eq. ( 105) we can write where   are the rows of the vector given by Latin American Journal of Solids and Structures, 2024, 21(1), e521 19/31 (110)

NUMERICAL EXAMPLES
Inspired in the work of Driemeier et al. (2005), Greco et al. (2006) and in the analyses of Carniel et al. (2015), we study two and three-dimensional trusses subjected to large displacements and strains under an unstable behavior (snap-through) along the strain path.The structures achieve high accelerations and require inclusion of inertia forces.Hence, within the framework of finite elements, we require to consider the kinematics of finite strains associated with the proper material model and transient geometric nonlinear formulation to obtain a response consistent with the physical phenomena involved.These structures are analyzed with different material laws and stress-strain pairs for material description.
For the viscoelastic analyses, we considered bars made of a hypothetical polymeric material introduced by Keramat and Ahmadi (2012) and used by Abdelrahman and El-Shafei (2021) with parameters shown in Table 2.As a numerical exercise this material was also analyzed neglecting the viscous effect in purely elastic simulations.In this case, in order to find the corresponding coefficient of elasticity, we considered the complete release of the dampers, leading to a regular series association of the springs.Thus, the value of the coefficient of elasticity considered for the elastic model corresponds to 10 9 MPa.The material parameters for hyperelastic bars correspond to natural rubber (Endo et al., 2021;Ogden, 1972).These parameters are given in Table 3. Damping can be introduced by inclusion of the Rayleigh damping matrix which, when considered, is taken as 10M.To simplify and standardize the nomenclature in the graphs that display results in this section, we have chosen to use the abbreviations Cauchy-Log, Eng -Eng and 2PK-GL, which correspond to Cauchy stress and logarithmic strain (and strain rate), engineering stress and strain (and strain rate) and Second Piola-Kirchhoff stress and Green-Lagrange strain (and strain rate) material laws, respectively.In addition, to make it easier to compare the behavior of the different stress measures over time, they were all converted to the Cauchy stress using Eq. ( 19) to (20).Viscous strain rate measures were all converted to logarithmic viscous strain rates, using Eq. ( 97).
Finally, all the algorithms used here have been validated with the results presented by Abdelrahman and El-Shafei (2021) for small strains.The results obtained for both small and large strains were consistent with those expected.

Two-bar symmetric structure
The first problem consists of the classical two-bar truss subjected to a vertical load (Figure 13).The simulation time is divided into two parts: for 0 ≤  ≤ 0.2, the timesteps length is Δ = 10 −4  and for 0.2 ≤  ≤ 1.0, Δ = 8 × 10 −5 .Figure 14 presents the vertical displacement history of node 2 when both bars are viscoelastic, using the properties of Table 2. Comparing the plots obtained for each of the linear relationships, it is noticeable that after the snap-through, the behavior for the three models is similar, although displaced with respect to time (Figs.14(a) and (b)).On the other hand, the oscillatory motion of the 2PK-GL material model appears first in relation to the others, followed by the Eng-Eng and the Cauchy-Log material models.As time evolves, the displacement values get closer, coinciding at 0.33s.From then on, the Cauchy-Log model presents a smaller displacement in relation to the others, tending to 0.33 m in the vertical direction.The 2PK-GL model is less rigid, tending to a displacement of 0.32 m (Fig. 14(c)).
Figure 15 shows that the Cauchy stress curves for the different constitutive laws present a behavior similar to the displacements in the sense that the curves are shifted with respect to time.The snap-through happens first for the 2PK-GL material model, which shows the lowest Cauchy stress value in compression, followed by the Eng-Eng material model and the Cauchy-Log material model.
Finally, Figure 16 presents the comparison between total and viscous logarithmic strain rates over time.Note that the viscous contribution is much smaller than the elastic one for all three material laws.Furthermore, while the rate of total deformation is nearly symmetric in traction and compression, this is not the case for the viscous strain rate, where the rate of traction is much higher than that of compression.Over time, the total strain rate and the viscous strain rate tend to 0.

Two elastic bars
As mentioned in the introduction of this section, we performed a numerical exercise in which the viscous part of the material was disregarded in the model described by Eq. ( 68).The displacements of node 2 are depicted in Fig. 17 and the Cauchy stress in Fig. 18.The main difference between the curves is the amplitude of the displacements, with the one corresponding to the Cauchy-Log model being larger, followed by the Eng -Eng model and the 2PK-GL model.The opposite trend is observed in the stress amplitude, as Fig. 18 shows increasing amplitude and nominal values in the same sequence of the material models.

Two hyperelastic bars
We now consider both members composed of the hyperelastic material defined by the properties in Table 3.This material model has no dissipation in the formulation, but the introduction of Rayleigh damping matrix, where  = 10 causes both the displacement of node 2 and the element stresses to decrease considerably over time.Note that in Figure 19 it is not possible to identify the snap-through.This is because, given the material properties, the snap-through occurs very quickly.For illustration, some intermediate configurations of the truss are shown for the undamped case.When damping is considered, the positions are different.

One bar viscoelastic and one bar hyperelastic
When the truss is composed of one viscoelastic and one hyperelastic members, node 2 moves in both directions and two vertical peak amplitudes can be identified, as shown in Fig. 20.
In this example, the choice of the material law has a negligible influence on the results found for both, the displacement, and the stress values over time.For this reason, we have chosen to display only the graphs that adopt the Cauchy-Log constitutive relation.Figure 21 shows the displacement of node 2 in the x-y plane as time elapses.Figure 21(a) displays a sequence of geometric configurations for the instants  = 0,  = 0.0068,  = 0.013,  = 0.7766 and  = 0.7786, where the path followed by node 2 is shown in blue (show the node numbers in the figure).The initial position of the truss is depicted in black.As time elapses, node 2 starts to move down and rightwards, producing traction in the viscoelastic bar and compression in the hyperelastic bar, until a snap-through occurs in the instant t = 0.068s (this configuration is depicted in purple).From then on, the strain energy stored in the bars is released, causing node 2 to move to the left of the initial position.As the stiffness of the hyperelastic bar is extremely low compared to the viscoelastic one, node 2 tends to stabilize oscillating with respect to a vertical line passing through node 1 (configurations depicted in red, orange and green in the plot).Notice that in this region, after a while, the length of the viscoelastic bar remains practically unchanged, and dissipation becomes negligible.Hence, the oscillatory motion does not vanish and the length of the hyperelastic bar continues to change along time.Figure 21(b) shows the same curve as above, now in 3 dimensions, but taking time into account.
Due to the aforementioned behavior, the dissipation effect fades out rapidly and the choice of the viscoelastic material law does not influence much the result.Figure 22 shows that the magnitude of the stress suffered by the viscoelastic element is much higher than that of the hyperelastic one.

Three-dimensional star truss
The second problem consists of a 3D truss whose central node is subjected to a vertically applied force.As in the previous example, the simulation time is divided into two parts: for 0 ≤  ≤ 0.2 we consider the timesteps length Δ = 10 −4  and for 0.2 ≤  ≤ 1.0, Δ = 8 × 10 −5  (see Figure 23).2).In the beginning, all the curves follow a similar pattern and present a lag according to the material law considered, as displayed in the detail given in Fig. 24(b).As time goes by, the models start to distance from each other.The material model 2PK-GL presents the highest stiffness, while the Cauchy-Log material model is the less stiff.The star-shaped truss has two snap-throughs that occur very close to each other, which are highlighted in Fig. 24 by a diamond and a square marker, respectively.The first snap-through occurs at  = 0.0044, followed by the second one at  = 0.0047.From then on, the structure starts to oscillate, with the amplitude of the movement decreasing over time.The viscous strain rate values (Figures 29 to 31) are much lower than those of the total strain rates, indicating the predominance of elastic behavior during the process.Over time, up to  = 0.2, they tend to a value close to 1.0/, while the total strain rate is already zero, i.e., they start to act strongly in damping the oscillations.From the instant  = 0.2 -which corresponds to the instant when the applied load becomes constant -the viscous rates decrease to zero.

All bars hyperelastic
Considering the entire 3D truss composed of hyperelastic bars (Table 3), a drastic difference is observed for the vertical displacement of node 7, as expected.Figure 32(a) shows such displacement for an undamped structure and the same displacement but with the introduction of a damping matrix  = 10.Notice that while elements 1 and 23 only suffer traction, element 10 suffers traction and compression, as it is not so demanded due to its position in the truss.

Mixed viscoelastic and hyperelastic bars
We now consider elements 1, 5, 8, 18, 19 e 22 being hyperelastic and the remaining ones viscoelastic.Once again, there is no significant difference between the displacement history of node 7 for different material models in the viscoelastic laws.However, in the presence of bars made of different materials distributed non symmetrically, horizontal displacements occur in directions X and Y, as shown in Figure 34.The initial position is marked with a red diamond while the final position is marked with a red square.It is noticeable that both frequency and amplitude of the stress oscillations are considerably higher for the viscoelastic bars compared to the hyperelastic ones.Furthermore, the difference in stiffness between the two materials is evident, since the stress values acting on the viscoelastic bar are much higher than on the hyperelastic bar. Figure 37 shows the geometric behavior of the 3D truss over time.

CONCLUDING REMARKS
In this paper, we present a finite element formulation for trusses made of different elastic and viscoelastic materials in transient nonlinear analyses.In the formulation, we emphasize the fact that the energy-conjugate pair considered in Latin American Journal of Solids and Structures, 2024, 21(1), e521 29/31 the VWP need not necessarily be the same as the stress-strain measure related by the constitutive equation and that the choice of the conjugate pair does not affect the results, regardless of the choice of the material law used.In this sense, we study different material behaviors by linearly relating different stress and strain measures, which we consider admissible if the compressive and tractive strains are moderate.We extend the derived geometric and material nonlinear formulation to account for rate effects using a generalized Kelvin rheological model.The numerical integration of the viscous strain rate is given by explicit uncoupled equations and is simpler than other methods found in the literature.We also present the Ogden hyperelastic model in the framework of the proposed formulation with dissipation introduced using a Rayleigh matrix proportional only to the mass.Transient problems are solved using the mean acceleration method.Two numerical examples involving different arrangements of trusses with members made of different materials are studied.The results obtained for the numerical problems studied are highly dependent on the properties of the materials adopted for the bars and the material law chosen to describe them.It is noteworthy that the usual practice of adopting the same stress-strain pair as a work-conjugate measure and for the material behavior is not always the most appropriate approach, since the former is arbitrary, and the latter should match experimental results.Furthermore, in the case of finite strains, the constitutive model does not always provide a physically realistic result.In the compressive case, the only material law that has a coherent behavior for both elastic and viscoelastic materials is the linear relationship between Cauchy stress and logarithmic strain (and strain rate).
The article is focused on bar elements as an attempt to observe the implications of different choices of nonlinear material models in a simple finite element setting.Interesting insights are obtained which can be useful for extending the analyses to continuum finite elements.For instance, Figs.14(a) and (c) of example 8.1.1 show that in the beginning of the phenomenon, the Cauchy-Log and the 2PK-GL material models yield the stiffest and the less stiff structures, respectively.However, this condition reverted completely along time, which is not directly intuitive.In example 8.2.1, another interesting issue can be remarked: according to the Hooke-like material model adopted, the viscous strain rate changes its qualitative pattern of evolution.Figures 29(b), 30(b) and 31(b) show increasing, constant and decreasing values for the Cauchy-Log, Eng-Eng and 2PK-GL material models, respectively.The relative simplicity of the bar element formulation, which does not demand complexities associated to element technology, allows focusing on structural effects of particular choices for the material model.

Figure 1 :
Figure 1: Global and corotational local systems of reference for the bar.

Figure 2 :
Figure 2: External and internal virtual work performed forces   and   due to a virtual displacement ().

Figure 3 :
Figure 3: Unphysical behavior for large compressive strains (for stretch approaching zero) in the cases of   =   and  2 =   .

Figure 6 :
Figure 6: A constant stress suddenly applied at instant  0 , the resulting strain profile in time and the corresponding creep-compliance function.

Figure 7 :
Figure 7: Compliance  = /  is independent of the creep stress level.

Figure 8 :
Figure 8: Compliance   =   /  for finite strains considering the linear relationship between (a) Green-Lagrange strain Second Piola-Kirchhoff stress and (b) engineering strain and stress, both with  = 0.

Figure 12 :
Figure 12 : Cauchy stress  versus stretch  1 for different values of the Poisson ratio  in the viscoelastic model: (a)   = ̇  and (b)  2 = ̇  .

Figure 14 :
Figure 14: 2D truss: displacement history of node 2 for both bars viscoelastic.(a) Whole time range evaluated (from 0 to 1 seconds), (b) detail in the neighborhood of the snap-throughs and (c) the moment when curves change position.

Figure 15 :
Figure 15: 2D truss: (a) Cauchy stress history for both bars viscoelastic; (b) detail of the plot between 0.15s and 0.3s.

Figure 16 :
Figure 16: 2D truss: (a) Total logarithmic strain rate and (b) viscous logarithmic strain rate histories of node 2 for both bars viscoelastic.

Figure 21 :
Figure 21: 2D truss structure with hyperelastic and viscoelastic materials: relationship between X and Y displacements

Figure 23 :
Figure 23: Three-dimensional star truss structure: initial geometry and corresponding (a) nodes and elements, (b) front and perspective views and (c) load history.

Figure 24 :
Figure 24: 3D truss: (a) vertical displacements of node 7 for all bars viscoelastic; (b) detail of the plot in the oscillatory region.

Figure 25
Figure 25 shows some configurations and the time at which they occur when the Cauchy-Log material model governs the behavior of the bars.

Figure 26 :
Figure 26: 3D truss: (a) Cauchy stress history for element 1.All bars viscoelastic; (b) detail of the plot in the oscillatory region.

Figure 27 :
Figure 27: 3D truss: (a) Cauchy stress history for element 10.All bars viscoelastic; (b) detail of the plot in the oscillatory region.

Figure 28 :
Figure 28: 3D truss: (a) Cauchy stress history for element 23.All bars viscoelastic; (b) detail of the plot in the oscillatory region.
Figure 32(b) shows a detail of the plot in the initial 0.3 seconds, where damping clearly reduces considerably the amplitude of the oscillations.Also, Figure 33 shows the evolution of the Cauchy stress in elements 1, 10 and 23 for the hyperelastic model with and without damping.Latin American Journal of Solids and Structures, 2024, 21(1), e521 27/31

Figure 32 :
Figure 32: 3D truss: (a) vertical displacements of node 7 for all bars hyperelastic; (b) detail of the plot in the first 0.3 seconds.

Figure 34 :
Figure 34: 3D truss: Relationship between the displacements of node 7 in the (a) X and Z, (b) Y and Z and (c) X and Y directions.Mixed viscoelastic and hyperelastic bars.

Figures
Figures 35 and 36 present the Cauchy stresses for the viscoelastic element 4 and for the hyperelastic element 1 respectively, which geometrically occupy similar positions in the truss.Once again, the choice of the constitutive model for the viscoelastic material does not influence the results.While Figures35(a) and 36(a) show the general behavior of the Cauchy stress curves along time, Figs.35(b) and 36(b) detail the corresponding plots from the beginning up to 0.1 seconds.It is noticeable that both frequency and amplitude of the stress oscillations are considerably higher for the viscoelastic bars compared to the hyperelastic ones.Furthermore, the difference in stiffness between the two materials is evident, since the stress values acting on the viscoelastic bar are much higher than on the hyperelastic bar.Figure37shows the geometric behavior of the 3D truss over time.

Figure 35 :
Figure 35: Mixed 3D truss: Cauchy stress for viscoelastic element 4; (b) detail of the plot in the first 0.1 seconds.

Figure 36 :
Figure 36: Mixed 3D truss: (a) Cauchy stress for hyperelastic element 1; (b) detail of the plot in the first 0.1 seconds.

Figure 37 :
Figure 37: 3D truss structure with hyperelastic and viscoelastic materials: behavior over time.