Transient analysis of trusses considering nonlinear elastic and viscoelastic material models

Abstract

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.

Keywords Hyperelasticity; Viscoelasticity; Transient analysis

Graphical Abstract

1 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).

2 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 $L_0,A_0$ and $V_0$ in the undeformed configuration, and $L,A$ and $V$ in the current configuration. We define a global system of reference with origin $O$ and base vectors ${\{\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3\}$, and a local corotational system of reference with origin at the beginning of the bar and base vectors $\left\{{\mathit{e}}_1^C,{\mathit{e}}_2^{\mathit{C}},{\mathit{e}}_3^C\right\}$. 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 $\mathit{X}=\left(X,Y,Z\right)$ and ${\mathit{X}}_l=\left(X_l,Y_l,Z_l\right)$. Accordingly, in the current configuration, the coordinates of the same material point are given by the vectors $\mathit{x}=\left(x,y,z\right)$ and ${\mathit{x}}_l=\left(x_l,y_l,z_l\right)$, as shown in Figure 1.

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

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

${\mathit{x}}_l=\left(x_l,y_l,z_l\right)=\left(X_l+u_l,Y_l+v_l,Z_l+w_l\right)$ (1)

where $\left(u_l,v_l,w_l\right)$ 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)

${\mathit{x}}_l=\left({\lambda }_1{X}_l,{\lambda }_2{Y}_l,{\lambda }_3{Z}_l\right)$ (2)

where

${\lambda }_1=\frac{d{x}_l}{d{X}_l}\mathrm{}$ (3)

and

${\lambda }_2={\lambda }_3.$ (4)

Using the logarithmic strain measure to define the Poisson ratio, according to

$\nu =-\mathrm{}\frac{\ln \left({\lambda }_2\right)}{\ln \left({\lambda }_1\right)}=-\mathrm{}\frac{\ln \left({\lambda }_3\right)}{\ln \left({\lambda }_1\right)}$,(5)

it turns out that

${\lambda }_2={\lambda }_3={\lambda }_1^{-\nu }$.(6)

Hence, the relation between the updated and initial cross section areas is

$dA=d{y}_{l}d{z}_l=\left({\lambda }_{2}d{Y}_l\right)\left({\lambda }_{3}d{Z}_l\right)={\lambda }_2{\lambda }_{3}d{Y}_{l}d{Z}_l={\lambda }_1^{-2\nu }d{A}_0$ (7)

and the strain gradient with respect to the corotational local system of reference is given by

${\mathit{F}}_C=\frac{\partial {\mathit{x}}_l}{\partial {\mathit{X}}_l}=\left[\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right]=\left[\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_1^{-\nu }& 0\\ 0& 0& {\lambda }_1^{-\nu }\end{array}\right]$ (8)

whose Jacobian is $J_{{\mathit{F}}_c}=det\left({\mathit{F}}_c\right)={\lambda \mathrm{}}_1^{-(2\nu +1)}$.

From the polar decomposition theorem, it turns out that the rotated right Cauchy strain tensor is

$\mathit{C}={\mathit{F}}_C^T{\mathit{F}}_{\mathit{C}}={\mathit{U}}_C^2,$ (9)

where the elongation tensor ${\mathit{U}}_c$ is used to define the Seth-Hill family of strain tensors

${\mathit{E}}_C^{\left(m\right)}=\frac{1}{m}\left({\mathit{U}}_C^{\left(m\right)}-\mathit{I}\right)=\left[\begin{array}{ccc}\frac{1}{m}\mathrm{}\left({\lambda }_1^{\left(m\right)}-1\right)& 0& 0\\ 0& \frac{1}{m}\mathrm{}\left[{\left({\lambda }_1^{-\nu }\right)}^{\left(m\right)}-1\right]& 0\\ 0& 0& \frac{1}{m}\mathrm{}\left[{\left({\lambda }_1^{-\nu }\right)}^{\left(m\right)}-1\right]\end{array}\right]$.(10)

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

${\epsilon }^{\left(m\right)}=\left\{\begin{array}{cc}\frac{1}{m}\left({\lambda }_1^m-1\right),& m\ne 0\\ \ln \left({\lambda }_1\right),& m=0\end{array}\right.$ (11)

with $m$ defining different measures. Some measures of strain are

• Green-Lagrange strain ($m=2$):

  ${\epsilon }^{GL}={\epsilon }^{\left(2\right)}=\frac{1}{2}\left({\lambda }_1^2-1\right)$(12)

• Rotated engineering (Biot) strain ($m=1$):

  ${\epsilon }^{RE}={\epsilon }^{\left(1\right)}={\lambda }_1-1$(13)

• Logarithmic strain ($m=0$):

  ${\epsilon }^L={\epsilon }^{\left(0\right)}=\ln \left({\lambda }_1\right)$(14)

3 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,

$\delta W_{ext}=\delta W_{int}+\delta W_m,$ (15)

where $\delta W_{ext}$ is the virtual work performed by the external load, $\delta W_{int}$ is the virtual work performed by internal forces and $\delta W_m$ 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

$\delta W_{ext}=\delta W_{int}+\delta W_m+\delta W_c$,(16)

where $\delta W_c$ 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

$\delta W_{int}={\int }_0^{L_0}{F}^{int}\left({\mathit{X}}_l\right)d\left(\delta u_l\right)={\int }_0^{L\mathrm{}}{F}^{int}\left({\mathit{x}}_l\right)d\left(\delta u_l\right)$ (17)

with $F^{int}$ being the axial internal force along the element length.

Figure 2
External and internal virtual work performed forces $F^{ext}$ and $F^{int}$ due to a virtual displacement $\delta u\left(L\right)$.

Within each bar, the axial stress can be defined by alternative measures, for instance

• the Cauchy stress:

  $\sigma =F^{int}/A$;(18)

• the rotated engineering or First Piola-Kirchhoff stress:

  ${\sigma }^{RE}=F^{int}/A_0={\lambda }_1^{-2\nu }\sigma$;(19)

• the Second Piola-Kirchhoff stress:

  ${\sigma }^{2PK}={\lambda }_1^{-\left(2\nu +1\right)}\sigma$.(20)

These stress measures can be confirmed by the application of the classical continuum mechanics expressions

$\mathit{P}=J_{{\mathit{F}}_c}{\mathit{F}}_c^{-1}\mathit{\sigma }={\lambda }_1^{-\left(2\nu -1\right)}\left[\begin{array}{ccc}{\lambda }_1^{-1}& 0& 0\\ 0& {\lambda }_1^{\nu }& 0\\ 0& 0& {\lambda }_1^{\nu }\end{array}\right]\left[\begin{array}{ccc}\sigma & 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]=\left[\begin{array}{ccc}{\sigma }^{RE}={\lambda }_1^{-2\nu }\sigma & 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$ (21)

and

$\mathit{S}=\mathit{P}{\mathit{F}}_c^{-T}=\left[\begin{array}{ccc}{\lambda }_1^{2\nu }\sigma & 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{ccc}{\lambda }_1^{-1}& 0& 0\\ 0& {\lambda }_1^{\nu }& 0\\ 0& 0& {\lambda }_1^{\nu }\end{array}\right]=\left[\begin{array}{ccc}{{\sigma }^{2PK}=\lambda }_1^{-(2\nu +1)}\sigma & 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$ (22)

where $\mathit{P}$ and $\mathit{S}$ 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)

$\delta W_{int}={\int }_{V_0}{\sigma }^{2\mathrm{P}\mathrm{K}}\delta {\epsilon }^{GL}d{V}_0$.(23)

The damping and inertia virtual work parcels mentioned in Eq. (16) are respectively

$\delta W_c={\int }_{V_0}c\dot{u}\delta u_{l}d{V}_0$ (24)

and

$\delta W_m={\int }_{V_0}\rho \mathrm{}\ddot{u}\delta u_{l}d{V}_0$ (25)

where $\mathit{c}\mathit{}$and $\mathit{\rho }$ are the specific damping and mass coefficients.

4 SPACE DISCRETIZATION OF THE PRINCIPLE VIRTUAL WORK

For each bar we define the interpolation matrix

$\mathit{N}\left(\xi \right)=\left[\begin{array}{ccc}\begin{array}{cc}{N}_1\left(\xi \right)& 0\end{array}& \begin{array}{cc}0& N_2\left(\xi \right)\end{array}& \begin{array}{cc}0& 0\end{array}\end{array}\right]$ (26)

based on a natural coordinate $\xi ϵ[-1,1]$ along the element axis. This matrix is used to approximate coordinates, displacements, velocities and accelerations by

${\mathit{X}}_l^{\mathrm{*}}=\mathit{N}{\stackrel{-}{\mathit{X}}}_l,$ (27)

${\mathit{x}}_l^{\mathrm{*}}=\mathit{N}{\stackrel{-}{\mathit{x}}}_l,$ (28)

${\mathit{u}}_l^{\mathrm{*}}=\mathit{N}{\stackrel{-}{\mathit{u}}}_l$,(29)

${\dot{\mathit{u}}}_l^{\mathrm{*}}=\mathit{N}\dot{{\stackrel{-}{\mathit{u}}}_l}$,(30)

${\ddot{\mathit{u}}}_l^{\mathrm{*}}=\mathit{N}{\ddot{\stackrel{-}{\mathit{u}}}}_l$,(31)

where the asterisk stands for continuum expressions being discretized and ${\stackrel{-}{\mathit{X}}}_l,{\stackrel{-}{\mathit{x}}}_l,{\stackrel{-}{\mathit{u}}}_l,{\dot{\stackrel{-}{\mathit{u}}}}_l$ and ${\ddot{\stackrel{-}{\mathit{u}}}}_l$ are the nodal values of the corresponding element vectors. For a linear two-noded bar, the interpolation functions are

$N_1\left(\xi \right)=\frac{1}{2}\left(1-\xi \right)$,(32)

$N_2\left(\xi \right)=\frac{1}{2}\left(1+\xi \right)$.(33)

4.1 Discretization of the virtual internal work

Replacing Eq. (12) and (13) into (23), the virtual internal work can be discretized by

$\delta W_{int}={\int }_{-1}^1{\left(\delta {\stackrel{-}{\mathit{u}}}_l\right)}^T{{\mathit{B}}^{\mathrm{*}}}^T\left({\stackrel{-}{\mathit{X}}}_l\right){\lambda }_1{\sigma }^{2PK}J\left({\stackrel{-}{\mathit{X}}}_l\right)A_{0}d\xi$,(34)

implying

${\mathit{q}}_l={\int }_{-1}^1{{\mathit{B}}^{\mathbf{*}}}^{\mathit{T}}\left({\stackrel{-}{\mathit{X}}}_l\right){\lambda }_1{\sigma }^{2PK}J\left({\stackrel{-}{\mathit{X}}}_l\right)A_{0}d\xi$,(35)

with

${\mathit{B}}^{\mathrm{*}}\left({\stackrel{-}{\mathit{X}}}_l\right)=\frac{1}{J\left({\stackrel{-}{\mathit{X}}}_l\right)}\frac{d\left(\mathit{N}\left(\xi \right)\right)}{d\xi }=\frac{1}{J\left({\stackrel{-}{\mathit{X}}}_l\right)}\left[\begin{array}{ccc}\begin{array}{cc}-1& 0\end{array}& \begin{array}{cc}0& 1\end{array}& \begin{array}{cc}0& 0\end{array}\end{array}\right]$ (36)

and $J\left({\stackrel{-}{\mathit{X}}}_l\right)=\frac{d{\stackrel{-}{\mathit{X}}}_l}{d\xi }$.

The transformation of the internal force vector to the global system of reference is performed by

$\mathit{q}={\mathit{T}}^T{\mathit{q}}_l$,(37)

where $\mathit{T}$ is the rotation matrix, composed by the direction cosines (Bathe, 2006).

4.2 Discretization of the inertia virtual work

Replacing Eq. (29) and (31) in Eq. (25), we obtain

$\delta W_m={\int }_{V_0}\rho \ddot{u}\delta u_{l}d{V}_0=\delta {\stackrel{-}{\mathit{u}}}_l^T{\int }_{-1}^1{\rho }^{\left(e\right)}\mathit{N}{\left({\stackrel{-}{\mathit{X}}}_l\right)}^T\mathit{N}\left({\stackrel{-}{\mathit{X}}}_l\right)A_{0}J\left({\stackrel{-}{\mathit{X}}}_l\right)d\xi {\ddot{\stackrel{-}{\mathit{u}}}}_l=\delta {\stackrel{-}{\mathit{u}}}_l^T{\mathit{m}}_e{\ddot{\stackrel{-}{\mathit{u}}}}_l$,(38)

where ${\mathit{m}}_e$ is the consistent element mass matrix and ${\rho }^{\left(e\right)}$ is the specific mass of the material of the element. Further,

$\delta W_m=\delta {\stackrel{-}{\mathit{u}}}_l^T{\mathit{m}}_e{\ddot{\stackrel{-}{\mathit{u}}}}_l=\delta {\stackrel{-}{\mathit{u}}}_l^T{\mathit{f}}_l^m,$ (39)

where

${\mathit{f}}_l^m={\mathit{m}}_e{\ddot{\stackrel{-}{\mathit{u}}}}_l$.(40)

The global mass matrix is assembled as usual,

$\mathit{M}=\underset{e=1}{\overset{nel}{\bigwedge }}{\mathit{m}}_e$,(41)

where $nel$ is the total number of elements and $\mathrm{\Lambda }$ is the assembly operator.

4.3 Discretization of the damping virtual work

Introducing Eq. (29) and (30) in Eq. (24), it results that

$\delta W_c={\int }_{V_0}c\dot{u}\delta u_{l}d{V}_0=\delta {\stackrel{-}{\mathit{u}}}_l^T{\int }_{-1}^1{c}^{\left(e\right)}\mathit{N}{\left({\stackrel{-}{\mathit{X}}}_l\right)}^T\mathit{N}\left({\stackrel{-}{\mathit{X}}}_l\right)A_{0}J\left({\stackrel{-}{\mathit{X}}}_l\right)d\xi {\dot{\stackrel{-}{\mathit{u}}}}_l=\delta {\stackrel{-}{\mathit{u}}}_l^T{\mathit{c}}_l{\dot{\stackrel{-}{\mathit{u}}}}_l$,(42)

where ${\mathit{c}}_l$ is the consistent element damping matrix. Moreover,

$\delta W_c=\delta {\stackrel{-}{\mathit{u}}}_l^T{\mathit{f}}_l^c$,(43)

where

${\mathit{f}}_l^c={\mathit{c}}_l{\dot{\stackrel{-}{\mathit{u}}}}_l$.(44)

As for the mass, the global damping matrix is obtained applying the assembly operator over all the elements,

$\mathit{C}=\underset{e=1}{\overset{nel}{\bigwedge }}{\mathit{c}}_e$,(45)

where ${\mathit{c}}_e={\mathit{T}}^T{\mathit{c}}_{\mathit{l}}\mathit{T}$ is the element damping matrix ${\mathit{c}}_l$ rotated to the global system of reference.

4.4 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 $\Delta t$. As the PVW must be valid for an arbitrary virtual kinematically admissible displacement, for an external force in the corotational local system of reference ${\mathit{f}}_{n+1}^{ext}$ we obtain that the discretized version of the equilibrium equations is

${\mathit{q}}_{n+1}\left({\stackrel{-}{\mathit{u}}}_{n+1}\right)+{\mathit{f}}_{n+1}^c\left({\dot{\stackrel{-}{\mathit{u}}}}_{n+1}\right)+{\mathit{f}}_{n+1}^m\left({\ddot{\stackrel{-}{\mathit{u}}}}_{n+1}\right)-{\mathit{f}}_{n+1}^{ext}=0$ (46)

valid for an arbitrary time or load step n+1. Summing up element contributions we can assemble the global form of Eq. (46)(45), obtaining a residual equation for force balance as

${\mathit{r}}_{n+1}\left({\stackrel{-}{\mathit{u}}}_{n+1}\right)={\mathit{Q}}_{n+1}\left({\stackrel{-}{\mathit{u}}}_{n+1}\right)+\mathit{C}{\dot{\stackrel{-}{\mathit{u}}}}_{n+1}+\mathit{M}{\ddot{\stackrel{-}{\mathit{u}}}}_{n+1}-{\mathit{F}}_{n+1}^{ext}$,(47)

which vanishes at equilibrium and takes non-zero values when an out-of-balance force occurs.

5 LINEARIZATION OF THE EQUILIBRIUM EQUATIONS

Owing to the term ${\mathit{Q}}_{n+1}\left({\stackrel{-}{\mathit{u}}}_{n+1}\right)$, the residual is a nonlinear function of the global displacement vector, so that the solution of Eq. (47) requires the linearization

${\mathit{r}}_{n+1}^{\left(k+1\right)}\left({\stackrel{-}{\mathit{u}}}_{n+1}^{\left(k+1\right)}\right)={\mathit{Q}}_{n+1}^{\left(k\right)}\left({\stackrel{-}{\mathit{u}}}_{n+1}^{\left(k\right)}\right)+\frac{d{Q}_{n+1}^{\left(k\right)}}{{\stackrel{-}{\mathit{u}}}_{n+1}^{\left(k\right)}}\mathrm{\Delta }{\stackrel{-}{\mathit{u}}}_{n+1}^{\left(k+1\right)}+\mathit{C}{\dot{\stackrel{-}{\mathit{u}}}}_{n+1}+\mathit{M}{\ddot{\stackrel{-}{\mathit{u}}}}_{n+1}-{\mathit{F}}_{n+1}^{ext}=0$,(48)

where the derivative of internal force vector ${\mathit{Q}}_{n+1}$ with respect to the global displacement vector defines the tangent matrix ${\mathit{K}}_T$ and $\mathrm{\Delta }{\stackrel{-}{\mathit{u}}}_{n+1}^{\left(k+1\right)}={\stackrel{-}{\mathit{u}}}_{n+1}^{\left(k+1\right)}-{\stackrel{-}{\mathit{u}}}_{n+1}^{\left(k\right)}$.

5.1 Tangent matrix

Dropping the $n+1$ index for neatness, we can express the tangent matrix as

${\mathit{K}}_T^k=\frac{d\mathit{Q}\left({\stackrel{-}{\mathit{u}}}^k\right)}{d{\stackrel{-}{\mathit{u}}}^k}=\underset{e=1}{\overset{nel}{\bigwedge }}{\mathit{k}}_{T_e}$,(49)

where ${\mathit{k}}_{T_e}$ 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, ${\mathit{k}}_{T_1}$ and ${\mathit{k}}_{T_2}$ .

${\mathit{k}}_{T_e}={\mathit{k}}_{T_1}^{\left(e\right)}+{\mathit{k}}_{T_1}^{\left(e\right)}$,(50)

where

${\mathit{k}}_{T_1}^{\left(e\right)}=\frac{d{\mathit{T}}^T{\mathit{q}}_{\mathit{l}}}{d\stackrel{-}{\mathit{u}}}={\int }_{-1}^1{\mathit{B}}^T\left({\stackrel{-}{\mathit{X}}}_l\right){\sigma }^{2PK}\mathit{H}\mathit{B}\left({\stackrel{-}{\mathit{X}}}_l\right)A_{0}J\left({\stackrel{-}{\mathit{X}}}_l\right)d\xi$,(51)

$\mathit{B}\left({\stackrel{-}{\mathit{X}}}_l\right)$ and $\mathit{H}$ are given by (Muñoz-Rojas, 2023)

$\mathit{B}\left({\stackrel{-}{\mathit{X}}}_l\right)=\frac{1}{2J\left({\stackrel{-}{\mathit{X}}}_l\right)}\left[\begin{array}{c}-1\\ \mathrm{}0\mathrm{}\\ 0\end{array}\begin{array}{c}\mathrm{}0\\ \mathrm{}-1\mathrm{}\\ 0\end{array}\begin{array}{c}0\mathrm{}\\ 0\mathrm{}\\ -1\mathrm{}\end{array}\begin{array}{c}\mathrm{}1\mathrm{}\\ 0\mathrm{}\\ 0\mathrm{}\end{array}\begin{array}{c}0\\ 1\\ 0\end{array}\begin{array}{c}\mathrm{}0\mathrm{}\\ \mathrm{}0\\ \mathrm{}1\end{array}\right]$,(52)

$\mathit{H}=\left[\mathit{I}-{\left(\frac{1}{{\lambda }_1}\right)}^2\mathit{B}\left({\stackrel{-}{\mathit{X}}}_l\right)\stackrel{-}{\mathit{x}}{\stackrel{-}{\mathit{x}}}^T{\mathit{B}}^T\left({\stackrel{-}{\mathit{X}}}_l\right)\right]$,(53)

$\mathit{I}$ is the identity matrix and $\stackrel{-}{\mathit{x}}$ is the vector of nodal current coordinates with respect to the global system of reference. On the other hand,

${\mathit{k}}_{T_2}^{\left(e\right)}={\mathit{T}}^T\frac{d{\mathit{q}}_{\mathit{l}}}{d\stackrel{-}{\mathit{u}}}={\mathit{T}}^T\frac{d{\mathit{q}}_{\mathit{l}}}{d{\stackrel{-}{\mathit{u}}}_{\mathit{l}}}\frac{d{\stackrel{-}{\mathit{u}}}_{\mathit{l}}}{d\stackrel{-}{\mathit{u}}}={\mathit{T}}^T\frac{d{\mathit{q}}_{\mathit{l}}}{d{\stackrel{-}{\mathit{u}}}_{\mathit{l}}}\mathit{T}={\mathit{T}}^T{\mathit{k}}_{T_{2}l}^{\left(e\right)}\mathit{T}$ , (54)

with

${\mathit{k}}_{T_{2}l}^e={\int }_{-1}^1{{\mathit{B}}^{\mathrm{*}}}^T\left({\stackrel{-}{\mathit{X}}}_l\right)\tilde{E}{\mathit{B}}^{\mathrm{*}}\left({\stackrel{-}{\mathit{X}}}_l\right)A_{0}J\left({\stackrel{-}{\mathit{X}}}_l\right)d\xi$,(55)

where (Muñoz-Rojas, 2023)

$\tilde{E}={\lambda }_1^2{E}_T+{\sigma }^{2PK}$,(56)

and $E_T$ is the tangent modulus

$E_T=\frac{\partial {\sigma }^{2PK}}{\partial {\epsilon }^{GL}}$ ,(57)

where $\left({\sigma }^{2PK},{\epsilon }^{GL}\right)$ is the conjugate pair adopted in this work.

Notice that the differentiation of ${\sigma }^{2PK}$ causes a dependence on the material model adopted. This important issue will be further discussed in Sec. 7.

5.2 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 $a$ and $b$ in the linear combination between mass and the initial tangent stiffness (equal to the linear stiffness matrix)

   $\mathit{C}=a\mathit{M}+b\mathit{K}$;(58)

2. constant coefficients $a$ and $b$ in the linear combination between mass and the updated tangent stiffness

   $\mathit{C}=a\mathit{M}+b{\mathit{K}}_T$;(59)

3. updated coefficients $a_t$ and $b_t$ in the linear combination between mass and the updated tangent stiffness

   $\mathit{C}=a_t\mathit{M}+b_t{\mathit{K}}_T$.(60)

In Eq. (58), (59) and (60), $a$ and $b$are called Rayleigh constants and the lower index $t$ 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 $\mathit{K}$ or ${\mathit{K}}_{\mathit{T}}$ adopting $b=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 .

6 TIME DISCRETIZATION

The Newmark family of time discretization methods is based on the expressions

$\mathit{M}{\ddot{\stackrel{-}{\mathit{u}}}}_{n+1}+\mathit{C}{\dot{\stackrel{-}{\mathit{u}}}}_{n+1}+{\mathit{K}}_T{\stackrel{-}{\mathit{u}}}_{n+1}={\mathit{F}}_{n+1}^{ext}$,(61)

${\stackrel{-}{\mathit{u}}}_{n+1}={\stackrel{-}{\mathit{u}}}_n+\mathrm{\Delta }t{\dot{\stackrel{-}{\mathit{u}}}}_n+{\left(\mathrm{\Delta }\mathrm{}t\right)}^2\left[\left(\frac{1}{2}-{\beta }_N\right){\ddot{\stackrel{-}{\mathit{u}}}}_n+{\beta }_N{\ddot{\stackrel{-}{\mathit{u}}}}_{n+1}\right],$ (62)

${\dot{\stackrel{-}{\mathit{u}}}}_{n+1}={\dot{\stackrel{-}{\mathit{u}}}}_n+\mathrm{\Delta }t\left[\left(1-{\gamma }_N\right){\ddot{\stackrel{-}{\mathit{u}}}}_n+{\gamma }_N{\ddot{\mathrm{}\stackrel{-}{\mathit{u}}}}_{n+1}\right],$ (63)

where ${\beta }_N$ and ${\gamma }_N$ define the specific method adopted and $\mathrm{\Delta }t$ is the length of time steps. Isolating ${\ddot{\stackrel{-}{\mathit{u}}}}_{n+1}$ in Eq. (62), and substituting Eq. (62) and (63) in (61) we find that, for each iteration k, the effective stiffness matrix ${\widehat{\mathit{K}}}_T^k$ and the effective external force vector ${\widehat{\mathit{F}}}^k$ are defined as (Subbaraj and Dokainish, 1989)

${\widehat{\mathit{K}}}_T=\frac{1}{{\beta }_N{\left(\mathrm{\Delta }t\right)}^2}\mathit{M}+\frac{{\gamma }_N}{{\beta }_N\mathrm{\Delta }t}\mathit{C}+{\mathit{K}}_{\mathit{T}},$ (64)

${\widehat{\mathit{F}}}^k={\left({\mathit{F}}^{ext}\right)}^k+\mathit{M}\left[\frac{1}{{\beta }_{\mathit{N}}{\left(\mathbf{\Delta }\mathit{t}\right)}^2}{\stackrel{-}{\mathit{u}}}^{\mathit{k}}+\frac{1}{{\beta }_{\mathit{N}}\mathbf{\Delta }\mathit{t}}{\dot{\stackrel{-}{\mathit{u}}}}^{\mathit{k}}+\left(\frac{1}{2{\beta }_{\mathit{N}}}-1\right){\ddot{\stackrel{-}{\mathit{u}}}}^{\mathit{k}}\right]++\mathit{C}\left[\frac{1}{{\beta }_{\mathit{N}}\mathrm{\Delta }\mathit{t}}{\stackrel{-}{\mathit{u}}}^{\mathit{k}}+\left(\frac{{\gamma }_{\mathit{N}}}{{\beta }_{\mathit{N}}}-1\right){\dot{\stackrel{-}{\mathit{u}}}}^{\mathit{k}}+\frac{\mathrm{\Delta }\mathit{t}}{2}\left(\frac{{\gamma }_{\mathit{N}}}{{\beta }_{\mathit{N}}}-2\right){\ddot{\stackrel{-}{\mathit{u}}}}^{\mathit{k}}\right],$ (65)

and the solution of the resulting global system via the Newton-Raphson method gives

${\widehat{\mathit{K}}}_T^k\mathrm{\Delta }{\stackrel{-}{\mathit{u}}}^{k+1}=-\widehat{\mathit{r}}\left({\stackrel{-}{\mathit{u}}}^{k+1}\right)$ (66)

${\stackrel{-}{\mathit{u}}}^{k+1}={\stackrel{-}{\mathit{u}}}^k+\mathrm{\Delta }{\stackrel{-}{\mathit{u}}}^{k+1}$ (67)

In Eq. (66), $\widehat{\mathit{r}}\left({\stackrel{-}{\mathit{u}}}^{\mathit{k}}\right)$ is the residual given by Eq. (48) with the effective external force in Eq. (65) in place of ${\mathit{F}}_{}^{ext}$. 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 ${\beta }_N=1/4$ and ${\gamma }_N=1/2$.

7 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 $\left({\sigma }^{2PK},{\mathrm{\epsilon }}^{GL}\right)$ to express the internal virtual work, for any given stress-strain pair chosen to describe the material behavior, say $\left({\sigma }^{*},{\mathrm{\epsilon }}^{*}\right)$, the tangent modulus will have the form

$E_T=\frac{\partial {\sigma }^{2PK}}{\partial {\epsilon }^{GL}}=\frac{\partial {\sigma }^{2PK}}{\partial {\sigma }^{\mathrm{*}}}\frac{\partial {\sigma }^{\mathrm{*}}}{\partial {\epsilon }^{\mathrm{*}}}\frac{\partial {\epsilon }^{\mathrm{*}}}{\partial {\epsilon }^{GL}}.$ (68)

7.1 Linear relationship between nonlinear stress and strain measures

Particularly, for a constitutive law that relates linearly the arbitrary stress and strain measures ${\sigma }^{*}$ and ${\epsilon }^{*}$,

${\sigma }^{\mathrm{*}}=E{\epsilon }^{\mathrm{*}}$ (69)

$E_T=\frac{\partial {\sigma }^{2PK}}{\partial {\epsilon }^{GL}}=\frac{\partial {\sigma }^{2PK}}{\partial {\sigma }^{\mathrm{*}}}E\frac{\partial {\epsilon }^{\mathrm{*}}}{\partial {\epsilon }^{GL}}.$ (70)

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 ${\mathrm{\sigma }}^{2PK}=E{\mathrm{\epsilon }}^{GL}$ (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 $({\sigma }^{2PK},{\epsilon }^{GL})$.

Constitutive law	$E_T$
${\sigma }^{RE}=E{\epsilon }^{RE}$	$E_T={\mathrm{\lambda }}_1^{-2}\left(E-{\mathrm{\sigma }}^{2PK}\right)$
${\sigma }^{2PK}=E{\epsilon }^{GL}$	$E_T=E$
$\mathrm{\sigma }=E{\epsilon }^L$	$E_T={\lambda }_1^{-2}\left(E{\lambda }_1^{-(2\nu +1)}-\mathrm{}\left(2\nu +1\right){\sigma }^{2PK}\right)$

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 ${\mathrm{\lambda }}_1=0$). The exception is $\mathrm{\sigma }=E{\mathrm{\epsilon }}^L$, as shown in Figure 3. Notwithstanding, for moderate strains (represented by vertical lines in the range between ${\lambda }_1=0.6$and ${\lambda }_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 ${\mathrm{\sigma }}^{RE}=E{\mathrm{\epsilon }}^{RE}$ can model an auxetic material $(\nu =-0.1)$ respecting the physical condition that the Cauchy stress $\mathrm{\sigma }\to -\mathrm{\infty }$ when ${\mathrm{\lambda }}_1\to 0.$ As expected, when small strain occurs, a linear relationship between any different stress-strain pairs leads practically to the same material behavior.

Figure 3
Unphysical behavior for large compressive strains (for stretch approaching zero) in the cases of ${\sigma }^{RE}=E{\epsilon }^{RE}$ and ${\sigma }^{2PK}=E{\epsilon }^{GL}$.

Figure 4
Cauchy stress $\sigma$ versus stretch ${\lambda }_1$ for different values of the Poisson ratio $\nu$: (a) ${\sigma }^{RE}=E{\epsilon }^{RE}$ and (b) ${\sigma }^{2PK}=E{\epsilon }^{GL}$.

7.2 Hyperelasticity (Ogden's model)

When the work done by the stresses during a strain process depends only on the initial state at time $t_0$ and the final configuration at time $t$, 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), ${\mathrm{\lambda }}_1,{\mathrm{\lambda }}_2,{\mathrm{\lambda }}_3$, where the elastic strain energy is given by

$W=\sum _{i=1}^N\frac{{\mu }_i}{{\alpha }_i}\left({\lambda }_1^{{\alpha }_i}+{\lambda }_2^{{\alpha }_i}+{\lambda }_3^{{\alpha }_i}-3\right),$ (71)

and ${\mu }_i$ and ${\mathrm{\alpha }}_i$ are material parameters. To obtain a physically consistent condition, the following constraints must be satisfied

${\mu }_i{\mathrm{\alpha }}_i>0,$ (72)

$\sum _{i=1}^N{\mu }_i{\alpha }_i=2\mu ,$ (73)

where $\mu$ is the shear modulus in the reference (undeformed) configuration. It is usual to consider that hyperelastic materials are incompressible, in which case

$det\left({\mathit{F}}_c\right)=J_{{\mathit{F}}_c}={\lambda }_1{\lambda }_2{\lambda }_3=1,$ (74)

where ${\mathit{F}}_c$ was defined in Eq. (8).

Introducing Eq. (6) into (74) we verify that for uniform uniaxial stretch, the incompressibility condition is met when $\nu =0.5$, so that the strain energy is given by

$W\left({\lambda }_1\right)=\sum _{i=1}^N\frac{{\mu }_i}{{\alpha }_i}\left({\lambda }_1^{{\alpha }_i}+2{\lambda }_1^{-{\alpha }_i/2}-3\right).$ (75)

Moreover, the Second Piola-Kirchhoff stress is obtained differentiating the strain energy with respect to the Green-Lagrange strain,

${\sigma }^{2PK}=\frac{dW}{d{\epsilon }^{GL}}=\frac{dW}{d{\lambda }_1}\frac{d{\lambda }_1}{d{\epsilon }^{GL}}={\sum }_{i=1}^N{\mu }_i\left({\lambda }_1^{{\alpha }_i-2}+2{\lambda }_1^{-\frac{{\alpha }_i}{2}-2}\right).$ (76)

Figure 5 shows that this incompressible material law respects the physical constraint that the Cauchy stress $\sigma \to -\mathrm{\infty }$ when ${\mathrm{\lambda }}_1\to 0$, hence allowing its use for large compressive strains. For comparison purpose, the relationships of ${\mathrm{\sigma }}^{RE}=E{\mathrm{\epsilon }}^{RE}$ and ${\mathrm{\sigma }}^{2PK}=E{\mathrm{\epsilon }}^{GL}$ are also displayed for $\nu =0.5$. The vertical lines at values ${\lambda }_1=0.6$and ${\lambda }_1=1.4$ delimit a range of moderate stretches, where all material laws provide coherent, although different values.

Figure 5
Cauchy stress $\sigma$ versus stretch ${\lambda }_1$ for incompressible hyperelastic and linear stress-strain relationships ($\nu =0.5)$.

Equation (76) is replaced in the expressions (35) and (51) of ${\mathit{q}}_l$ and ${\mathit{k}}_{T_1}$. Differentiation of Eq. (76) provides the tangent modulus defined in Eq. (57)

$E_T=\frac{d{\mathrm{\sigma }}^{2PK}}{d{\mathrm{\epsilon }}^{GL}}={\sum }_{i=1}^N{\mu }_i\left[\left({\mathrm{\alpha }}_i-2\right){\mathrm{\lambda }}_1^{{\mathrm{\alpha }}_i-4}+\left(\frac{{\mathrm{\alpha }}_i}{2}+2\right){\mathrm{\lambda }}_1^{-\frac{{\mathrm{\alpha }}_i}{2}-4}\right]$,(77)

which must be introduced in the expression (54) for ${\mathit{k}}_{T_2}$.

7.3 Viscoelasticity

7.3.1 Linear Viscoelasticity

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

$\sigma \left(t\right)={\int }_0^{t}G\left(t-\tau \right)\dot{\epsilon }\left(\tau \right)d\tau$ (78)

and

$\epsilon \left(t\right)={\int }_0^{t}J\left(t-\mathrm{\tau }\right)\dot{\mathrm{\sigma }}\left(\mathrm{\tau }\right)d\tau ,$ (79)

where $G\left(t\right)$ and $J\left(t\right)$ 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 ${\mathrm{\sigma }}_c$ is suddenly applied at time $t_0$ and kept constant, while the stress is observed over time, that is

$\mathrm{\sigma }\left(t\right)=H\left(t-t_0\right){\mathrm{\sigma }}_c$ (80)

and

$\dot{\mathrm{\sigma }}\left(t\right)=\frac{dH}{dt}\left(t-t_0\right){\mathrm{\sigma }}_c=\delta \left(t-t_0\right){\mathrm{\sigma }}_c$,(81)

where $H\left(t-t_0\right)$ and $\delta \left(t-t_0\right)$ are the Heaviside and Dirac delta functions respectively, and ${\mathrm{\sigma }}_c$ is the creep stress.

Replacing Eq. (81) in (79) and applying the filter property, it turns out that

$\epsilon \left(t\right)=J\left(t-t_0\right){\sigma }_c,$ (82)

or

$J\left(t-t_0\right)=\frac{\epsilon \left(t\right)}{{\sigma }_c},$ (83)

where frequently time $t_0$ is taken as $t_0=0$ (see Figure 6).

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

Notice that the compliance is defined by Eq. (83) and in linear viscoelasticity this function does not depend on the creep stress ${\mathrm{\sigma }}_c$ applied. A typical curve of linear creep-compliance (nonlinear in time) for any value of ${\sigma }_C$ is given in Figure 7.

Figure 7
Compliance $J=\epsilon /{\sigma }_c$ is independent of the creep stress level.

7.3.2 Nonlinear Viscoelasticity

When we work with finite strains, the use of the pairs $\left({\sigma }_c^{2PK},{\epsilon }^{GL}\right)$ and $\left({\sigma }_c^{RE},{\epsilon }^{RE}\right)$ 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

${\epsilon }^{GL}\left(t\right)=J\left(t-t_0\right){\sigma }_c^{2PK}.$ (84)

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

$e^{2\nu {\epsilon }^L+1}\left(e^{2\nu {\epsilon }^L}-1\right)=2J\left(t-t_0\right){\sigma }_c,$ (85)

which can be expressed as

$g\left({\epsilon }^L\right)=e^{2\nu {\epsilon }^L+1}\left(e^{2{\epsilon }^L}-1\right)-2J\left(t-t_0\right){\mathrm{\sigma }}_c=0.$ (86)

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

$e^{2\nu {\epsilon }^L}\left(e^{{\epsilon }^L}-1\right)=J\left(t-t_0\right){\mathrm{\sigma }}_c,$ (87)

which can be rewritten as

$g\left({\epsilon }^L\right)=e^{2\nu {\epsilon }^L}\left(e^{2{\epsilon }^L}-1\right)-J\left(t-t_0\right){\mathrm{\sigma }}_c=0.$ (88)

In both cases (Eq.(87) and (88)) we can define the compliance as

$J_{NL}\left(t-t_0\right)=\frac{Root\mathrm{}of\mathrm{}g\left({\epsilon }^L\right)}{{\mathrm{\sigma }}_c}$.(89)

where ${\sigma }_C$ 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.

Figure 8
Compliance $J_{NL}={\epsilon }^L/{\sigma }_c$ for finite strains considering the linear relationship between (a) Green-Lagrange strain Second Piola-Kirchhoff stress and (b) engineering strain and stress, both with $\nu =0$.

7.3.3 Generalized Kelvin-Voigt rheological model

The simple Kelvin-Voigt rheological model consists of a spring of elastic constant $E_i$ arranged in parallel with a dashpot of viscosity constant ${\eta }_i$, 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).

Figure 9
Simple and generalized Kelvin-Voigt rheological models.

The stress-strain relationship for a single Kelvin block $i$ is given by

$\sigma \left(t\right)=E_i{\epsilon }_i^{ve}\left(t\right)+{\eta }_i{\dot{\epsilon }}_i^{ve}\left(t\right)$,(90)

where ${\epsilon }^{ve}\left(t\right)$ and ${\dot{\epsilon }}^{ve}\left(t\right)$ are the viscoelastic strain and its time rate. Equation (90) can be rewritten as

${\dot{\epsilon }}_i^{ve}\left(t\right)=\left(\frac{1}{{\tau }_i{E}_i}\right)\sigma \left(t\right)-\left(\frac{1}{{\tau }_i}\right){\epsilon }^{ve}\left(t\right)$,(91)

where ${\mathrm{\tau }}_i=\frac{{\mathrm{\eta }}_i}{E_i}$ is the relaxation time of the rheological block.

Because we have a series arrangement, in each rheological block it holds that

$\sigma \left(t\right)={\sigma }_{S,i}\left(t\right)+{\sigma }_{D,i}\left(t\right)$,(92)

where

$\left\{\begin{array}{ccc}{\sigma }_{S,i}\left(t\right)=E_i{\epsilon }_i^{ve}\left(t\right);& i=\mathrm{1,2},...,nrb& \left(\mathrm{s}\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\right)\\ {\sigma }_{D,i}\left(t\right)={\eta }_i{\dot{\epsilon }}_i^{ve}\left(t\right);& i=\mathrm{1,2},...,nrb& \left(\mathrm{d}\mathrm{a}\mathrm{s}\mathrm{h}\mathrm{p}\mathrm{o}\mathrm{t}\mathrm{}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\right)\end{array}\right.$ (93)

and $nrb$ is the number of rheological blocks involved in the material model. On the other hand, the stress of the sole spring is

$\sigma \left(t\right)=E_0\left[\epsilon \left(t\right)-{\epsilon }^{ve}\left(t\right)\right],$ (94)

where the viscous strain corresponds to the sum of the viscous strains of each rheological block, so that

$\sigma \left(t\right)=E_0\left[\epsilon \left(t\right)-\sum _{j=1}^{nrb}{\epsilon }_j^{ve}\left(t\right)\right].$ (95)

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

${\dot{\epsilon }}_i^{ve}\left(t\right)+\frac{1}{{\tau }_i}{\epsilon }_i^{ve}\left(t\right)+\frac{{\omega }_i}{{\tau }_i}{\epsilon }_i^{ve}\left(t\right)=\frac{{\omega }_i}{{\tau }_i}{\epsilon }_i\left(t\right),$ (96)

where ${\mathrm{\omega }}_i=E_0/E_i$.

In the case of viscous behavior, we define a linear relationship between stress and strain rate ($\sigma =\eta \dot{\epsilon }$), 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)

Figure 10
Purely viscous system.

$\left\{\begin{array}{cc}{\epsilon }^L=\ln \left({\lambda }_1\right)& \Rightarrow \mathrm{}{\dot{\epsilon }}^L=\frac{1}{{\lambda }_1}{\dot{\lambda }}_1\\ {\epsilon }^{RE}={\lambda }_1-1& \Rightarrow \mathrm{}{\dot{\epsilon }}^{RE}={\dot{\lambda }}_1\\ {\epsilon }^{GL}=\frac{1}{2}\left({\lambda }_i^2-1\right)& \Rightarrow {\mathrm{}{\dot{\epsilon }}^{GL}=\lambda }_1{\dot{\lambda }}_1\end{array}\right.$,(97)

which replaced in the Hooke-like material relationships lead to

$\left\{\begin{array}{ccc}\sigma =\eta {\dot{\epsilon }}^L& \Rightarrow \mathrm{}\sigma =\phi \left({\lambda }_1\right){\dot{\lambda }}_1,& \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\phi \left({\lambda }_1\right)=\frac{\eta }{{\lambda }_1}\mathrm{}\\ {\sigma }^{RE}=\eta {\dot{\epsilon }}^{RE}& \mathrm{}\Rightarrow \mathrm{}\sigma =\phi {\left({\lambda }_1\right)\dot{\lambda }}_1,& \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\phi \left({\lambda }_1\right)={\eta \lambda }_1^{2\nu }\\ {\sigma }^{2PK}=\eta {\dot{\epsilon }}^{GL}& \Rightarrow \mathrm{}\sigma =\phi \left({\lambda }_1\right)\mathrm{}{\dot{\lambda }}_1,& \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\phi \left({\lambda }_1\right)={\eta \lambda }_1^{2(\nu +1)}\end{array}\mathrm{}\right.$.(98)

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 $\phi \left({\lambda }_1\right)$ defined in Eq. (98) for the $(\sigma ,{\epsilon }^L)$, $\left({\sigma }^{RE},{\epsilon }^{RE}\right)$and $({\sigma }^{2PK},{\epsilon }^{GL})$ 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, ${\lambda }_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 $(\sigma ,{\epsilon }^L)$. However, for both constitutive laws $\left({\sigma }^{RE},{\epsilon }^{RE}\right)$ and $({\sigma }^{2PK},{\epsilon }^{GL})$, when the stretch tends to 0, the values of $\phi \left({\lambda }_1\right)$also tend to 0. Since, according to Eq. (98) the stresses are obtained by multiplying $\phi \left({\lambda }_1\right)$ by the strain rates, it follows that the stresses tend to zero independent of ${\dot{\lambda }}_1$ values, contrary to the expected behavior.

Figure 11
Unphysical behavior for large compressive strains (for stretch approaching zero) in the cases of ${\sigma }^{RE}=\eta {\dot{\epsilon }}^{RE}$ and ${\sigma }^{2PK}=\eta {\dot{\epsilon }}^{GL}$ considering $\eta =1$ and $\nu =0.5.$

Figure 12 shows the behavior of $\left({\sigma }^{2PK},{\epsilon }^{GL}\right)$ and $({\sigma }^{RE},{\epsilon }^{RE})$ for different values of $\nu$. Although $\left({\sigma }^{RE},{\epsilon }^{RE}\right)$ presents physical sense for auxetic materials, the $\left({\sigma }^{2PK},{\epsilon }^{GL}\right)$ constitutive law does not give coherent results for any of the Poisson ratio values tested when large strains in compression are considered.

Figure 12
: Cauchy stress $\sigma$ versus stretch ${\lambda }_1$ for different values of the Poisson ratio $\nu$ in the viscoelastic model: (a) ${\sigma }^{RE}=\eta {\dot{\epsilon }}^{RE}$ and (b) ${\sigma }^{2PK}=\eta {\dot{\epsilon }}^{GL}.$

7.3.4 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 $N$ steps of equal size $\mathrm{\Delta }t$. For each step $0<n<N$, the next step $n+1$ will be equal to the time elapsed until the previous time step $n$, $t_n$, plus the time interval, i.e., $t_{n+1}=t_n+\mathrm{\Delta }t$. Also, the total strain at step $n+1$ is known and therefore will be constant throughout this time step. Then for time step $n+1$, eq. (96) is such that

${\dot{\epsilon }}_{i,n+1}^{ve}+\frac{1}{{\tau }_i}{\epsilon }_{i,n}^{ve}+\frac{{\omega }_i}{{\tau }_i}{\epsilon }_n^{ve}=\frac{{\omega }_i}{{\tau }_i}{\epsilon }_{n+1},\mathrm{}i=\mathrm{1,2},\dots ,nrb$,(99)

where

${\dot{\epsilon }}_{i,n+1}^{ve}=\frac{{\epsilon }_{i,n+1}^{ve}-{\epsilon }_{i,n}^{ve}}{\mathrm{\Delta }t}$.(100)

Replacing Eq. (100) into (99), for each rheological block $i$,

${\epsilon }_{i,n+1}^{ve}={\epsilon }_{i,n}^{ve}+\mathrm{\Delta }t\frac{{\omega }_i}{{\tau }_i}{\left(\epsilon \right)}_{n+1}-\frac{\mathrm{\Delta }t}{{\tau }_i}{\epsilon }_{i,n}^{ve}-\mathrm{\Delta }t\frac{{\omega }_i}{{\tau }_i}{\epsilon }_n^{ve}$ (101)

with

${\epsilon }_n^{ve}=\sum _{j=1}^{nrb}{\epsilon }_{j,n}^{ve}$ (102)

or

${\mathrm{\tau }}_i{\mathrm{\epsilon }}_{i,n+1}^{ve}=\left[{\mathrm{\tau }}_i-\mathrm{\Delta }t\right]{\mathrm{\epsilon }}_{i,n}^{ve}+\mathrm{\Delta }t{\mathrm{\omega }}_i{\mathrm{\epsilon }}_{n+1}-\mathrm{\Delta }t{\mathrm{\omega }}_i{\sum }_{j=1}^{nrb}{\mathrm{\epsilon }}_{j,n}^{ve}=\left[{\tau }_i-\mathrm{\Delta }t\left(1-{\omega }_i\right)\right]{\epsilon }_{i,n}^{ve}+\mathrm{\Delta }t{\omega }_i{\epsilon }_{n+1}-\mathrm{\Delta }t{\omega }_i{\sum }_{j=1,j\ne i}^{nrb}{\epsilon }_{j,n}^{ve}$.(103)

Eq. (103) can be written as a linear system $\mathit{A}{\mathbf{\epsilon }}_{n+1}^{ve}=\mathit{b}$ where

$\mathit{A}=\left[\begin{array}{c}\begin{array}{c}{\tau }_1\\ 0\end{array}\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\mathrm{}\begin{array}{c}\begin{array}{c}0\\ {\tau }_2\end{array}\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\mathrm{}\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ {\tau }_3\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\mathrm{}\begin{array}{cc}\begin{array}{c}...\\ ...\end{array}& \begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}...\\ ...\\ ...\end{array}& \begin{array}{c}0\\ ⋮\\ {\tau }_{nrb}\end{array}\end{array}\right]$, ${\mathbf{}\mathbf{\epsilon }}_{n+1}^{ve}=\left[\begin{array}{c}\begin{array}{c}{\epsilon }_{1,n+1}^{ve}\\ {\epsilon }_{2,n+1}^{ve}\end{array}\\ \begin{array}{c}{\epsilon }_{3,n+1}^{ve}\\ ⋮\\ {\epsilon }_{nrb,n+1}^{ve}\end{array}\end{array}\right],$(104)

$\mathit{b}=\left[\begin{array}{cc}\begin{array}{cc}{\mathrm{\tau }}_1-\mathrm{\Delta }t\left(1+{\mathrm{\omega }}_1\right)& -{\mathrm{\omega }}_1\mathrm{\Delta }t\\ -{\mathrm{\omega }}_2\mathrm{\Delta }t& {\mathrm{\tau }}_2-\mathrm{\Delta }t\left(1+{\mathrm{\omega }}_2\right)\end{array}& \begin{array}{cc}...\mathrm{}& -{\mathrm{\omega }}_1\mathrm{\Delta }t\\ ...\mathrm{}& -{\mathrm{\omega }}_2\mathrm{\Delta }t\end{array}\\ \begin{array}{cc}⋮\mathrm{}& ⋮\\ -{\mathrm{\omega }}_{nrb}\mathrm{\Delta }t\mathrm{}& -{\mathrm{\omega }}_{nrb}\mathrm{\Delta }t\end{array}& \begin{array}{cc}\mathrm{}⋮& \mathrm{}⋮\\ \mathrm{}...& {\mathrm{}\mathrm{\tau }}_{nbr}-\mathrm{\Delta }t\left(1+{\mathrm{\omega }}_{nrb}\right)\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\epsilon }_{1,n}^{ve}\\ {\epsilon }_{2,n}^{ve}\end{array}\\ \begin{array}{c}⋮\\ {\epsilon }_{nrb,n}^{ve}\end{array}\end{array}\right]+\mathrm{\Delta }\mathrm{}t\left[\begin{array}{c}\begin{array}{c}{\omega }_1\\ {\omega }_2\end{array}\\ \begin{array}{c}⋮\\ {\omega }_{nrb}\end{array}\end{array}\right]{\epsilon }_{n+1}$ (105)

Note that matrix $\mathit{A}$ 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,

${\mathrm{\epsilon }}_{i,n+1}^{ve}=\frac{1}{{\mathrm{\tau }}_i}\left[\left({\tau }_i-\mathrm{\Delta }t\right){\epsilon }_{i,n}^{ve}+\mathrm{\Delta }t{\omega }_i\left({\epsilon }_{n+1}-{\sum }_{j=1}^{nrb}{\epsilon }_{j,n}^{ve}\right)\right],\mathrm{}i=\mathrm{1,2},\dots ,nrb.$ (106)

The consistent tangent modulus of a viscoelastic material (not to confuse with the tangent modulus defined in Eq.(68)) in the time step $n+1$ is, by definition,

$E_{n+1}^{ve}=\frac{\partial {\sigma }_{n+1}}{\partial {\epsilon }_{n+1}}$ (remember that the asterisks have been dropped: $\frac{\partial {\sigma }_{n+1}^{*}}{\partial {\epsilon }_{n+1}^{*}}$ in Eq. (68))(107)

where ${\mathrm{\sigma }}_{n+1}=E_0\left({\mathrm{\epsilon }}_{n+1}-{\mathrm{\epsilon }}_{n+1}^{ve}\right)$. Then

$E_{n+1}^{ve}=E_0\left(1-{\sum }_{j=1}^{nrb}\frac{\partial {\epsilon }_{j,n+1}^{ve}}{\partial {\epsilon }_{n+1}}\right).$ (108)

Using Eq. (105) we can write

$E_{n+1}^{ve}=E_0\left(1-{\sum }_{j=1}^{nrb}{a}_j\right),$ (109)

where $a_j$ are the rows of the vector given by

$\mathit{a}=\left[\begin{array}{c}\begin{array}{c}{\omega }_1/{\tau }_1\\ {\omega }_2/{\tau }_2\end{array}\\ \begin{array}{c}{\omega }_3/{\tau }_3\\ ⋮\\ {\omega }_{nrb}/{\tau }_{nrb}\end{array}\end{array}\right]\mathrm{\Delta }t$.(110)

8 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.

Table 2
Viscoelastic material properties (Abdelrahman and El-Shafei, 2021)

$\nu$	$\rho \left[kg/m^3\right]$	$E_0\left[MPa\right]$	$E_i\left[MPa\right]$	${\mathrm{\tau }}_i\left[s\right]$
0.50	12000.00	1.0035 E +13	1.086957 E +11	1.101 E -03
			9.049774 E +09	3.0115 E -02
			1.281558 E +09	1.50784 E -01

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.

Table 3
Hyperelastic material properties of natural rubber (Endo et al., 2021)

${\mu }_i\left[\mathrm{M}\mathrm{P}\mathrm{a}\right]$	${\mathrm{\alpha }}_i$
0.7.7817 E +05	2.7971
-1.1229 E +04	-2.7188
1.269 E -01	10.505
1.6169 E +07	0.33382

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.

8.1 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\le t\le 0.2s$, the timesteps length is $\mathrm{\Delta }t={10}^{-4}s$ and for $0.2s\le t\le 1.0s$, $\mathrm{\Delta }t={8\times 10}^{-5}s$.

Figure 13
(a) Two-dimensional truss structure: initial geometry and (b) load history.

8.1.1 Two viscoelastic bars

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 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 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.

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

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.

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

8.1.2 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.

Figure 17
2D truss: Displacement history of node 2 for both elastic bars with (a) Cauchy-Log, (b) Eng-Eng and (c) 2PK-GL material models.

Figure 18
2D truss: Cauchy stress history for both elements with (a) Cauchy-Log, (b) Eng-Eng and (c) 2PK-GL material models.

8.1.3 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 $\mathit{C}=10\mathit{M}$ 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.

Figure 19
2D truss:(a) Displacement and (b) stress history of node 2 for both bars hyperelastic (undamped and damped).

8.1.4 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.

Figure 20
2D truss structure with hyperelastic and viscoelastic materials: (a) horizontal and (b) vertical displacements of node 2.

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 $t=0$, $t=0.0068$, $t=0.013$, $t=0.7766$ and $t=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.

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

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.

Figure 22
2D truss: Cauchy stress history (a) Viscoelastic element (b) Hyperelastic element.

8.2 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\le t\le 0.2$ we consider the timesteps length $\mathrm{\Delta }t={10}^{-4}s$ and for $0.2\le t\le 1.0$, $\mathrm{\Delta }t=8\times {10}^{-5}s$ (see 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.

8.2.1 All bars viscoelastic

Figure 24(a) shows the vertical displacement of node 7 for the three Hooke-like material laws when all the bars are viscoelastic (Table 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 $t=0.0044s,$followed by the second one at $t=0.0047s$. From then on, the structure starts to oscillate, with the amplitude of the movement decreasing over time.

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 shows some configurations and the time at which they occur when the Cauchy-Log material model governs the behavior of the bars.

Figure 25
3D truss: some configurations for the viscoelastic case (Cauchy-Log material model).

Owing to their different positions in the truss layout, elements 1, 10 and 23 are representative of the whole structure. Figures 26 to 28 display the Cauchy stress for such elements, where we notice the difference between the compression and tension peaks of the different stress measures. Again, all the curves show a similar behavior, although shifted, with a small difference in amplitude according to the material law considered, which is evinced by Figures 26(b) to 28(b). Cauchy stress in elements 1 and 23 increase until $t=0.2s$, when it becomes constant. This instant coincides with the moment when the load applied to the truss becomes constant as well.

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
3D truss: (a) Cauchy stress history for element 10. All bars viscoelastic; (b) detail of the plot in the oscillatory region.

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

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 $t=0.2s$, they tend to a value close to $1.0MPa/s$, while the total strain rate is already zero, i.e., they start to act strongly in damping the oscillations. From the instant $t=0.2s$ - which corresponds to the instant when the applied load becomes constant - the viscous rates decrease to zero.

8.2.2 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 $\mathit{C}=10\mathit{M}$. 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.

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 33
3D truss: Stress history for element (a) 1, (b) 10 and (c) 23 with and without structural damping. All bars hyperelastic.

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.

8.2.3 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.

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 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 Figures 35(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. Figure 37 shows the geometric behavior of the 3D truss over time.

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

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
3D truss structure with hyperelastic and viscoelastic materials: behavior over time.

9 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 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 29
3D truss: (a) total and (b) viscous logarithmic strain rates for Cauchy-Log material model.

Figure 30
3D truss: (a) total and (b) viscous logarithmic strain rates for Eng-Eng material model.

Figure 31
3D truss: (a) total and (b) viscous logarithmic strain rate for 2PK-GL material model.

Acknowledgements

The authors acknowledge the financial support given by FAPESC grants number 2023TR000563.

• https://doi.org/10.1590/1679-78257918

References

• Abdelrahman, A.A., El-Shafei, A.G. (2021). Modeling and analysis of the transient response of viscoelastic solids, Waves in Random and Complex Media 31: 1990—2020.
• Ananthapadmanabhan, S., and Saravanan, U. (2020). Numerical techniques for solving truss problems involving viscoelastic materials, International Journal of Non-Linear Mechanics, 122: 103479.
• Araújo Junior, P.C., de Holanda, A.S., Parente, E., Evangelista, F. (2010). Dynamic Viscoelastic Analysis of Asphalt Pavements using a Finite Element Formulation, Road Materials and Pavement Design, 11:2, 409—433.
• Bathe, K. J. (2006). Finite element procedures. Prentice Hall.
• Belytschko, T., Liu, W., Moran, B., Elkhodary, K.I. (2014). Nonlinear finite elements for continua and structures, John Wiley & Sons, Ltd.
• Bonet, J., and Wood, R. D. (2008). Nonlinear continuum mechanics for finite element analysis, Cambridge university press.
• Chang, S. Y. (2004). Studies of Newmark method for solving nonlinear systems:(I) basic analysis, Journal of the Chinese Institute of Engineers, 27(5): 651-662.
• Carniel, T. A., Muñoz-Rojas, P. A., and Vaz Jr, M. (2015). A viscoelastic viscoplastic constitutive model including mechanical degradation: Uniaxial transient finite element formulation at finite strains and application to space truss structures., Applied Mathematical Modelling, 39(5-6: 1725-1739.
• Charney, F. A. (2008). Unintended consequences of modeling damping in structures, Journal of structural engineering, 134(4):581-592.
• Chazal, C. and Moutou Pitti, R. (2009). A new incremental formulation for linear viscoelastic analysis: creep differential approach, Journal of Theoretical and Applied Mechanics, 47(2): 397-409.
• Chazal, C. and Pitti, R. M. (2011). Integral approach for time dependent materials using finite element method, Journal of theoretical and applied mechanics, 49(4): 1029-1048.
• Chen, X., Duan, J., Li, Y. (2015). Mass proportional damping in nonlinear time-history analysis. In 3rd International Conference on Material, Mechanical and Manufacturing Engineering (IC3ME 2015) (pp. 567-571). Atlantis Press.
• Crisfield, M. A. (1991). Non-Linear Finite Element Analysis of Solids and Structures, Essentials. John Wiley and Sons.
• Douven, L. F. A., Schreurs, P. J. G. and Janssen, J. D. (1989). Analysis of viscoelastic behaviour of transversely isotropic materials, International journal for numerical methods in engineering, 28(4): 845-860.
• Driemeier, L., Proença, S. P. B., and Alves, M. (2005). A contribution to the numerical nonlinear analysis of three-dimensional truss systems considering large strains, damage and plasticity, Communications in nonlinear science and numerical simulation, 10(5): 515-535.
• Dunne, F. and Petrinic, N. (2005). Introduction to computational plasticity. OUP Oxford.
• Endo, V. T., Fancello, E. A., Muñoz-Rojas, P. A. (2021). A study on the computational effort of hyper-dual numbers to evaluate derivatives in geometrically nonlinear hyperelastic trusses, Journal of the Brazilian Society of Mechanical Sciences and Engineering, 43: 1-15.
• Evangelista Junior, F. (2006). Análise quasi-estática e dinâmica de pavimentos asfálticos. Universidade Federal do Ceará.
• Fonseca, F. M. and Gonçalves, P. B. (2022). Nonlinear behavior and instabilities of a hyperelastic von Mises truss, International Journal of Non-Linear Mechanics, 142: 103964.
• Greco, M., Gesualdo, F. A. R., Venturini, W. S., Coda, H. B. (2006). Nonlinear positional formulation for space truss analysis. Finite elements in analysis and design, 42(12): 1079-1086.
• Guggenberger, W. and Krenn, B. (2019). 3D truss model for large deformations and large strains: Nonlinear formulation and continuum mechanical aspects. In Advances in Engineering Materials, Structures and Systems: Innovations, Mechanics and Applications: Proceedings of the 7th International Conference on Structural Engineering, Mechanics and Computation (SEMC 2019), September 2-4, 2019, Cape Town, South Africa (p. 320). CRC Press.
• Holzapfel, G. A. (2002). Nonlinear solid mechanics: a continuum approach for engineering science, Kluwer Academic Publishers Dordrecht.
• Jehel, P., Léger, P. Ibrahimbegovic, A. (2014). Initial versus tangent stiffness‐based Rayleigh damping in inelastic time history seismic analyses, Earthquake Engineering and Structural Dynamics, 43(3): 467-484.
• Kaliske, M. and Rothert, H. (1997). Formulation and implementation of three-dimensional viscoelasticity at small and finite strains, Computational Mechanics, 19(3): 228-239.
• Keramat, A. and Ahmadi, A. (2012). Axial wave propagation in viscoelastic bars using a new finite-element-based method, Journal of Engineering Mathematics, 77(1): 105-117.
• Kühl, A., Muñoz‐Rojas, P. A., Barbieri, R., Benvenutti, I. J. (2017). A procedure for modeling the nonlinear viscoelastoplastic creep of HDPE at small strains. Polymer Engineering and Science, 57(2): 144-152.
• Muñoz-Rojas, P.A. (2023). Lecture notes.
• Nedjar, B. and Le Roy, R. (2013). An approach to the modeling of viscoelastic damage, Application to the long‐term creep of gypsum rock materials, International Journal for Numerical and Analytical Methods in Geomechanics, 37(9): 1066-1078.
• Ogden, R. W. (1972). Large deformation isotropic elasticity–on the correlation of theory and experiment for incompressible rubberlike solids, Proceedings of the Royal Society of London, A. Mathematical and Physical Sciences, 326(1567): 565-584.
• Pascon, J. P. (2016). Nonlinear analysis of hyperelastoplastic truss-like structures, Archive of Applied Mechanics, 86: 831-851.
• Rabelo, J. M., Becho, J. S., Greco, M. Cimini Jr, C. A. (2018). Modeling the creep behavior of GRFP truss structures with Positional Finite Element Method, Latin American Journal of Solids and Structures, 15.
• Rezaiee-Pajand, M., Masoodi, A. R., and Arabi, E. (2023). Geometric and Material Nonlinear Analyses of Trusses Subjected to Thermomechanical Loads, Structural Engineering International, 33(2): 302-313.
• Subbaraj, K. and Dokainish, M. (1989). A survey of direct time-integration methods in computational structural dynamics—II. Implicit methods, Computers and structures, 32(6): 1387-1401.
• Ward, I. M. and Sweeney, J. (2004). An introduction to the mechanical properties of solid polymers, John Wiley and Sons, Ltd.
• Zienkiewicz, O. C., Watson, M., King, I. P. (1968). A numerical method of visco-elastic stress analysis, International Journal of Mechanical Sciences, 10(10): 807-827.

Edited by

Editor: Marcílio Alves

Publication Dates

• Publication in this collection
  22 Jan 2024
• Date of issue
  2024

History

• Received
  08 Nov 2023
• Reviewed
  09 Dec 2023
• Accepted
  09 Dec 2023
• Published
  14 Dec 2023