Shear contribution on viscoelastic behavior of beams modeled using plane frame elements with Reissner kinematics

Abstract

The development of a numerical formulation is presented to describe viscoelastic behavior considering shear effects. The nonlinear positional formulation of the Finite Element Method is used considering plane frame elements with Reissner kinematics. The description of the viscoelastic behavior is considered through the adoption of a stress-strain relation based on Boltzmann rheological model. The used kinematics allows to describe the decoupling between the rotations and the displacements in element cross-sections. This approach allows to evaluate the contribution of shear effects in viscoelastic behavior in an original way. Based on the developments and the results obtained, it is possible to observe that strains and displacements due to viscoelastic behavior are significantly superior to the results obtained considering Bernoulli-Euler kinematics hypothesis. It is possible to notice a better agreement between the obtained numerical results and the results in the literature. Results obtained from the developed formulation allow the assessment of the shear effects on the viscoelastic behavior of plane frames.

Keywords
Viscoelasticity; Reissner kinematics; Finite element analysis; Positional formulation; Rheological models; Creep

Graphical Abstract

1 INTRODUCTION

The continuous search for engineering solutions with efficient performance reiterates the importance of researches related to the knowledge of complex mechanical behaviors and their more accurate computational simulations. Particularly in the structural engineering field, good strength/weight ratios are aimed, combined with safety, predictable and low cost structures. In this context, several researches are related to nonconventional materials behavior modeling and more realistic and precise analyses of structures and structural components, aiming to address the diverse demands of areas such as infrastructure, civil construction, mechanical industry, aerospace industry, among others. Nonetheless, complex analyses involve nonlinear behaviors, with few analytical solutions available and making the use of numerical methods necessary, most of them based on the Finite Element Method (FEM).

Among the nonlinear behaviors to which a quasi-static structure is subjected, it is possible to highlight the physical nonlinearities described by the constitutive or rheological relations. These relations are expressed by stress-strain equations that may include the dependence on specific variables, such as time, temperature, humidity, among others. These equations relate the components of the stress and strain tensors, enabling the description of various mechanical behaviors such as elastic, plastic, and viscous. In general, real materials can exhibit these three behaviors separately, depending on material properties and service conditions, or simultaneously, through intermediate behaviors between two or three of these, such as viscoelastic, elastoplastic, and viscoelastoplastic (Marques and Creus, 2012; Christensen, 2003; Findley et al., 1989).

In the present paper, interest is restricted to the viscoelastic mechanical behavior of solid materials. In this case, the relations between the components of stress and strain tensors must be represented by rheological relations which take into account the time variable and can be deduced from rheological models (Marques and Creus, 2012; Christensen, 2003; Findley et al., 1989).

The rheological models are schematic representations that combine elastic and viscous elements to obtain a physical interpretation of rheological relations. Thus, from the interpretation of these rheological models, it is possible to deduce relations between the components of the stress and strain tensors suitable for the considered mechanical behavior. Several researches are related to the modeling of viscoelastic structural behavior, considering various numerical formulations. Some of these researches adopt stress-strain relations deduced from rheological models, as demonstrated in the works of Aköz and Kadioǧlu (1999), Mesquita and Coda (2003), Liu et al. (2008), Panagiotopoulos et al. (2014), Carniel et al. (2015), Oliveira and Leonel (2017), Pascon and Coda (2017), Rabelo et al. (2018), and Kühl and Muñoz-Rojas (2019).

The interest of part of the recent researches in nonlinear behaviors and analysis methods is due to the developments and applications of micro-scale mechanical systems and nanostructures, as exposed in the works of Demir and Civalek (2017) and Ebrahimi et al. (2020). Some of these structures are nanobeams and compositions with nanobeams. These nanostructures are often made up of materials with complex behavior such as polymeric and composite materials due to their characteristics and mechanical requirements. Thus, some of these studies are dedicated to analyzing nanostructures with viscoelastic behavior, as can be seen in the works of Li et al. (2016), Xiao et al. (2017), Jalaei and Civalek (2019), and Dastjerdi et al. (2020).

The viscoelastic behavior is of special importance in polymeric materials and composite materials with polymeric matrix. Such materials have gained prominence in civil construction, mechanical industry, aerospace industry and others, mainly due to their good strength/weight ratio. For this reason, some works have presented studies referring to the analysis of the viscoelastic behavior of beams of Fiber Reinforced Polymer (FRP) (Sá et al., 2011a; Sá et al., 2011b; Alachek et al., 2019; Gand et al., 2013) and sandwich panels with polymeric core (Shenoi, et al., 1997; Galucio, et al., 2004; Garrido and Correia, 2018).

In addition, it is important to note that, due to the physical characteristics of materials that exhibit viscoelastic behavior, particularly for polymeric and composite materials with polymeric matrix, several experimental studies demonstrate that the shear effects on this behavior cannot be neglected, as demonstrated in the works of Bank and Mosallam (1992), Mottram (1993), Abdel-Magid et al. (2003), Shao and Shanmugam (2004), Sá et al. (2011a), and Garrido and Correia (2018). The shear effects can be preponderant if the material is considered almost incompressible (such as rubber-like materials) or even the only effect when volumetric deformations are completely suppressed (for instance, in geophysical applications). Riobom Neto et al. (2021) uses Boltzmann model in a Boundary Elements Formulation applied to plane and half-plane quasi-static viscoelastic problems considering only shear effects. Authors model the elastic behavior based on spherical tensors and the deviatoric constitutive relation based on Boltzmann model. Oliveira et al. (2018) investigate future structural behavior of two-dimensional structures considering Boltzmann model applied to creep analyses. Shear effects in a concrete viaduct segment is modeled.

In this context, this paper aims to present the development of a numerical formulation based on the positional formulation of the Finite Element Method and capable of describing the viscoelastic mechanical behavior in structures discretized by plane frame elements, considering the shear effects. To do this, the viscoelastic behavior is evaluated through the stress-strain relation deduced from the Boltzmann rheological model (Standard solid model), which takes into account the time variable. Thus, the choice of Boltzmann model is mainly due to its simplicity. As the model has only three parameters it is easy to do parametric analyses and isolate each parameter influence in the formulation. Rheological models obtained from associations of Kelvin and Maxwell models are interesting for solid with viscoelastic due to the ease of calibrating parameters based on experimental results. The calibration process for viscoelastic behavior may be a very complex task.

Meanwhile, the relevant effects of shear on viscoelastic behavior are taken into account in an original way by the use of plane frame elements with Reissner kinematics. The first reason to use the Reissner kinematics instead of a conventional 2D elastic model is related to its simplicity. Both regarding the number of degrees of freedom and the mathematical operations are simpler. For instance, in Becho and Greco (2021) a cylindrical pressure vessel made by viscoelastic material is simulated. The pressure is applied gradually and the radial displacement evolution along time is analyzed. Numerical results for radial displacements are compared with numerical results from Mesquita and Coda (2002) during inflation process. In Mesquita and Coda (2002) a bidimensional mesh with 200 free formulation discrete Kirchhoff triangle finite elements and the viscoelastic Boltzmann model were considered. Results obtained from Becho and Greco (2021), with only 10 plane frame finite elements, are in good agreement with results available in the literature. The second reason, not treated in current the manuscript, is regarding the model calibration and it is also a direct consequence from the first reason. Rabelo et al. (2018) presents a procedure to calibrate viscoelastic parameters for a Glass Fiber Reinforced Plastics (GFRP) structure using the Least Squares Method. It is known that experimental tests for different ultimate strength stress levels must be done for creep behavior modelling. The parameters calibration for numerical formulation can be done by relations between viscoelastic parameters and the stress level to which the material is exposed. Thus, suitable parameters were numerically adjusted to approximately match the experimental results obtained by experimental results. In Rabelo et al. (2018) the suitable viscoelastic parameters for each stress level were found considering quadratic functions of regression. The calibration process is exhausting and computationally expensive, even for a truss-like structure such as the one analyzed in Rabelo et al. (2018). Besides the authors' familiarity with the aforementioned positional formulation of the Finite Element Method, one of the motivations for using it in the present paper is due to its simplicity and applicability in analyses of nonlinear problems as highlight by the precursory authors of the formulation (Coda and Greco, 2004; Maciel and Coda, 2008). For this reason, although relatively recent, the positional formulation has been the object of study by several researchers, as demonstrated in the works of Coda and Paccola (2007), Greco and Ferreira (2009), Oliveira and Greco (2014), Sampaio et al. (2015), Cavalcante et al. (2017), Pascon and Coda (2017), Carrazedo and Coda (2017), Rabelo et al. (2018), Fernandes et al. (2018), and Ramos and Carrazedo (2020).

2 POSITIONAL FORMULATION OF THE FINITE ELEMENT METHOD

Originally presented in Coda and Greco (2004), the positional formulation of the Finite Element Method is based on variational concepts of the Minimal Total Potential Energy Principle. In this formulation, the main variables adopted to describe the finite element kinematics are the nodal positions of the structure in relation to a space-fixed coordinate system (Total Lagrangian description).

For quasi-static and conservative structural systems, as presented in Coda and Greco (2004), the functional of the total potential energy ($\mathrm{\Pi }$) can be determined by the Principle of Virtual Works as follows:

$\Pi =U+W$(2-1)

in which $U$ represents the total strain energy and $W$ represents the potential energy of external applied forces, given by

$U={\int }_V^{}udV={\int }_V^{}{\int }_{\epsilon }^{}\sigma \mathrm{}d\epsilon dV$(2-2)

$W=-\sum _{q=1}^{ngl}{P}_q{X}_q$(2-3)

The stress tensor is represented by $\sigma$, the strain tensor is represented by $\epsilon$, the degrees of freedom are represented by $X_q$, the equivalent loads related to each degree of freedom are represented by $P_q$, and $ndf$ represents the total number of degrees of freedom.

By the application of variational principles, it is possible to obtain the equilibrium configuration corresponding to the minimum of total potential energy functional. To do this, according to Dym and Shames (2013) the first variation of total potential energy functional must be zero as follows:

${\delta }^{\left(1\right)}\Pi =0$(2-4)

Considering the variations regarding the degrees of freedom (nodal parameters in positional formulation), one has:

${\delta }^{\left(1\right)}\Pi =\frac{\partial \Pi }{\partial X_q}\delta X_q=0$(2-5)

Considering that $\delta X_q$ represents an arbitrary variation in the degrees of freedom, a nonlinear system of ngl equations is obtained as follows:

$\frac{\partial \Pi }{\partial X_q}=\frac{\partial U}{\partial X_q}+\frac{\partial W}{\partial X_q}=0$ (where q = 1, 2, ..., ndf) (2-6)

It represents the Principle of Minimal Total Potential Energy (Crisfield, 1991). Since the guarantee of the application of this principle to provide the minimum value of the total potential energy can be based on the evaluation of the second derivative of the functional (or from the second variation of the functional).

Conceptually, before discretization, the structure has infinite degrees of freedom and, consequently, the system expressed by Equation (2-6) has infinite equations. After discretization, the number of degrees of freedom is reduced to the number of nodal parameters, making the numerical resolution of the system of equations possible.

Observing that the relation between the total strain energy and the nodal positions (generalized nodal parameters) is nonlinear, it is possible to infer that Equation (2-6) represents a system of nonlinear equations. In this case, this system needs to be solved by using an appropriate equation system solving method. From this, it is possible to determine the equilibrium position of a structure subjected to a specific loading.

In this paper, the system of nonlinear equations is solved using the Newton-Rapshon Iterative Method, briefly described below, considering the structure discretized by the nodal positions (generalized nodal parameters).

To describe the Newton-Rapshon Iterative Method, Equation (2-6) can be rewritten in a compact form as follows (Greco and Coda, 2006):

$\frac{\partial \Pi }{\partial X_q}=g_q\left(X\right)=0$(2-7)

Expanding Equation (2-7) through a first order Taylor series, one has:

$g_q\left(X\right)\cong g_q\left(X\right)+g_q{,}_r\left(X\right){∆X}_r\cong 0\mathrm{}$ (where q and r = 1, 2, ..., ndf) (2-8)

in which $g_q\left(X\right)$ represents the components of the vector of residues, $g_q{,}_r\left(X\right)$ represents the components of the Hessian matrix and ${∆X}_r$ represents the components of the nodal position correction vector:

$g_q\left(X\right)={U,}_q-P_q$(2-9)

${g_q,}_r\left(X\right)={U,}_{qr}$(2-10)

The solution for the nonlinear system of equation given by Equation (2-8) can be obtained by Newton-Raphson Iterative Method, summarized as follows:

(1)$X$ is taken as an undeformed configuration;

(2) The components of the vector of residues $g_q\left(X\right)$ are calculated using Eq. (2-9);

(3) The components of the Hessian matrix $g_q{,}_r\left(X\right)$ are calculated using Eq. (2-10);

(4) The components of the nodal position correction vector ${∆X}_r$ are calculated using Eq. (2-8);

(5) The adopted criterion of convergence is verified considering two situations:

(5.1) If the nodal position correction vector $∆X$ does not satisfy the adopted criterion of convergence, a new iterative step is done returning to item (2) with the nodal positions vector updated considering the correction vector ($X=X+\mathrm{}∆X$);

(5.2) If the nodal position correction vector $∆X$ satisfies the adopted criterion of convergence, the iterative process is finished and the equilibrium at deformed configuration is given by the nodal positions that were reached ($X$).

The convergence criterion is adopted based on a relation between Euclidian norms as follows:

$\frac{‖∆X‖}{‖X‖}\le tol$(2-11)

in which tol represents an adopted numerical tolerance.

From the basic formulation developed for the Principle of Minimal Total Potential Energy and for the Newton-Raphson Iterative Method, it is possible to obtain the equilibrium position of a quasi-static and conservative structural system subjected to a determined loading state. Therefore, it is necessary to particularize the formulation to take into account the implemented finite element kinematics and the mechanical behavior of interest.

The particularization of the positional formulation of the Finite Element Method is summarized in determining the total strain energy $U$ and its derivatives ${U,}_q$ and ${U,}_{qr}$ in relation to the nodal positions (generalized nodal parameters), present in Equations (2-9) and (2-10). The total strain energy must be described as a function of the stress-strain relation (characteristic of mechanical behavior) and the strain measure, which is described in terms of the kinematics of the considered finite element, expressed in terms of the shape functions and nodal positions.

3 STRESS-STRAIN RELATION

The differential form, presented in Equation (3-1), is generally used to describe viscoelastic behavior and stress-strain relations. However, the integration of Equation (3-1) is complex and, therefore, it is necessary to divide it into simpler models and proceed with the sum of strains in each model to obtain the total deformation (Christensen, 2003; Kassner et al., 2004) or use adequate numerical methods.

$\sigma +R_1\frac{\partial \sigma }{\partial t}+\dots +R_n\frac{{\partial }^n\sigma }{\partial t^n}=Q_0\epsilon +Q_1\frac{\partial \epsilon }{\partial t}+\dots +Q_n\frac{{\partial }^n\epsilon }{\partial t^n}$(3-1)

where the parameters 𝑅𝑛 and 𝑄𝑛 represent the material properties dependent on time, temperature, humidity and stress/strain levels. In the presented manuscript, temperature and humidity are not considered in the fomrulation.

In general, the rheological models used to represent the mechanical behavior of viscoelastic materials can be derived from the generalized Maxwell and Kelvin-Voigt models, depending on the definition of the parameters involved (Argyris et al., 1991). In order to determine the total strain energy, it is necessary to obtain the stress-strain relation for the model under consideration. In the case of viscoelastic behavior, a part of the stress components is associated with elastic behavior, while another portion is associated with viscous behavior. Thus, the stress-strain relation can be obtained considering the appropriate association between these parts obtained from the physical interpretation of the adopted rheological model.

The rheological model considered in the present paper is the Boltzmann model (Standard solid model), shown in Figure 1, as adopted in Mesquita and Coda (2001) and Oliveira and Leonel (2017). This model is capable of describing an instantaneous elastic behavior followed by a damped elastic behavior (viscoelastic behavior) with decreasing strain rate, typical of solid materials. For a sufficiently large period of time, the total viscoelastic strains converges to the instantaneous elastic strain response predicted by a linear elastic model with an elastic modulus equivalent to the series association of two elastic elements with longitudinal elastic modulus respectively equal to E₁ e E₂.

Figure 1
Boltzmann rheological model (Standard solid model).

The rheological relation for the Boltzmann model can be obtained considering that the total stress in the model is equal to the stress in the elastic element or the stress in the Kelvin-Voigt model (parallel association between an elastic element and a viscous element) and that the total strain of the model can be g iven by the sum of the strain of the elastic element and the strain of the Kelvin-Voigt model, expressed by:

${\sigma }_{ij}={\sigma }_{ij}^{\left(1\right)}={\sigma }_{ij}^{\left(2\right)}$(3-2)

${\epsilon }_{ij}={\epsilon }_{ij}^{\left(1\right)}+{\epsilon }_{ij}^{\left(2\right)}$(3-3)

in which indexes “(1)” and “(2)” refer to the single elastic element and to the Kelvin-Voigt model, respectively.

It is also important to note that the relations between stresses and strains in an elastic element and a viscous element in the isotropic case can be defined, respectively, by:

${\sigma }_{ij}=\left[\lambda {\delta }_{ij}{\delta }_{kl}+\mathrm{}\mu ({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})\right]{\epsilon }_{kl}$(3-4)

${\sigma }_{ij}=\left[\lambda {\alpha }_{\lambda }{\delta }_{ij}{\delta }_{kl}+\mathrm{}\mu {\alpha }_{\mu }({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})\right]{\dot{\epsilon }}_{kl}$(3-5)

in which ${\delta }_{ij}$ represents the Kronecker delta operator, $\lambda$ represents the first Lamé parameter, $\mu$ represents the second Lamé parameter, and ${\alpha }_{\lambda }$ and ${\alpha }_{\mu }$ represent the coefficients of viscosity. Since the first and second invariants of strain are given by:

$\lambda =\frac{E\nu }{\left(1+\nu \right)\mathrm{}\left(1-2\nu \right)}$(3-6)

$\mu =\frac{E}{2\left(1+\nu \right)}$(3-7)

in which $\nu$ represents the Poisson coefficient and E represents the longitudinal elastic modulus.

In a simplified way, for the isotropic case, the coefficients of viscosity ${\alpha }_{\lambda }$ and ${\alpha }_{\mu }$ can be approximated as equals ($\alpha ={\alpha }_{\lambda }={\alpha }_{\mu }$), as observed in Mesquita and Coda (2001). In addition, as observed in Oliveira and Leonel (2017), the product between the longitudinal elastic modulus and the coefficient of viscosity can be represent a viscosity parameter represented by $\eta$. This parameter represents the modulus of viscosity.

Thus, as presented in Mesquita and Coda (2001), considering Equations (3-2) to (3-7), the rheological relation of the Boltzmann model can be given by:

${\sigma }_{ij}=\frac{E_1{E}_2}{E_1+E_2}{{\stackrel{-}{C}}_{ijkl}\epsilon }_{kl}+\frac{\eta E_1}{E_1+E_2}{{\stackrel{-}{C}}_{ijkl}\dot{\epsilon }}_{kl}-\frac{\eta }{E_1+E_2}{\dot{\sigma }}_{ij}$(3-8)

in which ${\stackrel{-}{C}}_{ijkl}$ represents the undimensionalized constitutive elastic tensor, given by:

${\stackrel{-}{C}}_{ij}=\stackrel{-}{\lambda }{\delta }_{ij}{\delta }_{kl}+\mathrm{}\stackrel{-}{\mu }({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})$(3-9)

in which $\stackrel{-}{\lambda }$ represents the undimensionalized first Lamé parameter and $\stackrel{-}{\mu }$ represents the undimensionalized second Lamé parameter, expressed, respectively, by:

$\stackrel{-}{\lambda }=\frac{\lambda }{E}=\frac{\nu }{\left(1+\nu \right)\mathrm{}\left(1-2\nu \right)}$(3-10)

$\stackrel{-}{\mu }=\frac{\mu }{E}=\frac{1}{2\left(1+\nu \right)}$(3-11)

Similarly, it is possible to develop a stress-strain relation for other rheological models.

Considering the stress-strain relation, expressed by Equation (3-8), it is possible to determine the total strain energy, expressed by Equation (2-2), for that it is necessary to particularize the kinematics and the strain measure for the adopted finite element in the present paper. Thus, the following section presents the particularization of the kinematics and the strain measure for the plane frame element with Reissner kinematics.

4 PARTICULARIZATION OF THE POSITIONAL FORMULATION

In order to particularize the total strain energy, it is necessary to understand the kinematics (Reissner) of the considered finite element and its relationship with the adopted strain measure (Biot). To do this, it is necessary to map the configuration change of a finite element. Thus, in positional formulation, each finite element has its geometry mapped by the parameterization along length and height depending on, respectively, the dimensionless variables ${\xi }_1$ (varying from -1 to 1) and ${\xi }_2$ (varying from -1 to 1), as illustrated in Figure 2.

Figure 2
Geometry parameterization of a plane frame element with Reissner kinematics.

In the undeformed configuration a generic point $p(x\left({\xi }_1,{\xi }_2\right),y\left({\xi }_1,{\xi }_2\right))$ can be mapped by centroidal position and slope of the cross-section, as shown in Figure 2. Similarly, the same point in the deformed configuration, represented by $P(X\left({\xi }_1,{\xi }_2\right),Y\left({\xi }_1,{\xi }_2\right))$, can be mapped after the configuration change of the element. Since the slopes of the cross-sections can be determined by the directions of the unit vectors $\overrightarrow{n}$ and $\overrightarrow{N}$, respectively, in the undeformed and deformed configurations.

Thus, the $x$ and $y$ coordinates of a generic point $p$, in the undeformed configuration, can be expressed, respectively, as follows:

$x\left({\xi }_1,{\xi }_2\right)=\tilde{x}\left({\xi }_1\right)-\frac{h}{2}{\xi }_2\mathrm{}sen\mathrm{}\left(\theta \left({\xi }_1\right)\right)={\varphi }_n\left({\xi }_1\right)x_n-\frac{h}{2}{\xi }_2\mathrm{}sen\mathrm{}\left({\varphi }_n\left({\xi }_1\right){\theta }_n\right)$(4-1)

$y\left({\xi }_1,{\xi }_2\right)=\tilde{y}\left({\xi }_1\right)+\frac{h}{2}{\xi }_2\mathrm{}cos\mathrm{}\left(\theta \left({\xi }_1\right)\right)={\varphi }_n\left({\xi }_1\right)y_n+\frac{h}{2}{\xi }_2\mathrm{}cos\mathrm{}\left({\varphi }_n\left({\xi }_1\right){\theta }_n\right)$(4-2)

in which $\tilde{x}$ and $\tilde{y}$ represent coordinates of the common point between the cross-section and the neutral plane, θ represents the angle between the cross-section and the horizontal axis, and h represents the height of the element cross-section. Similar equations are obtained for the deformed configuration, in terms of the X and Y coordinates of a generic point P. The coordinates mapping of the points, both in the undeformed configuration and in the deformed configuration, can be rewritten in terms of nodal parameters and shape functions ${\varphi }_n\left({\xi }_1\right)$. In this paper, a plane frame finite element with four nodes is adopted and, consequently, a polynomial approximation of order 3 (cubic) is considered.

The rotation of a cross-section is approximated independently in the Reissner kinematics, as it is done for the other nodal parameters. Thus, the cross-sections do not necessarily remain perpendicular to the centroidal axis. Considering the geometry mapping of a plane frame element, it is possible to describe its deformation by the Biot strain measure. To do this, it is necessary to obtain the deformation gradient, which describes the relationships between transformations of the point coordinates of a domain when it passes from the undeformed configuration to the deformed configuration.

Thus, particularizing for the case of a finite element of plane frame with Reissner kinematics, the deformation gradient ($F_{ij}$) can be expressed, respectively, for the transformations from the dimensionless auxiliary configuration to the undeformed configuration (${{}^{0}F}_{ij}$) and from the dimensionless auxiliary configuration to the deformed configuration (${{}^{1}F}_{ij}$), as follows:

${{}^{0}F}_{ij}=\left[\begin{array}{cc}\partial x/\partial {\xi }_1& \partial x/\partial {\xi }_2\\ \partial y/\partial {\xi }_1& \partial y/\partial {\xi }_2\end{array}\right]=\left[\begin{array}{cc}{{}^{0}F}_{11}& {{}^{0}F}_{12}\\ {{}^{0}F}_{21}& {{}^{0}F}_{22}\end{array}\right]$ and ${{}^{1}F}_{ij}=\left[\begin{array}{cc}\partial X/\partial {\xi }_1& \partial X/\partial {\xi }_2\\ \partial Y/\partial {\xi }_1& \partial Y/\partial {\xi }_2\end{array}\right]=\left[\begin{array}{cc}{{}^{1}F}_{11}& {{}^{1}F}_{12}\\ {{}^{1}F}_{21}& {{}^{1}F}_{22}\end{array}\right]$(4-3)

From the mapping of the presented geometry and the deformation gradient it is possible to calculate the strains for the adopted finite element. To do this, in the present paper, the non-linear engineering strain measure is adopted, as presented in the works of Coda and Greco (2004) and Maciel et al. (2004), i.e. Biot strain measure for non-linear analysis is considered.

The stretches in the normal (longitudinal) direction for the transformation from auxiliary configuration to undeformed configuration (${{}^0\lambda }_{11}$) and for the transformation from auxiliary configuration to deformed configuration (${{}^1\lambda }_{11}$), can be expressed by:

${}^0{\lambda }_{11}=\left|{}^{0}F\mathrm{}\left[\begin{array}{c}1\\ 0\end{array}\right]\right|=\sqrt{{\left({}^0{F}_{11}\right)}^2+{\left({}^0{F}_{21}\right)}^2}$ and ${}^1{\lambda }_{11}=\left|{}^{1}F\mathrm{}\left[\begin{array}{c}1\\ 0\end{array}\right]\right|=\sqrt{{\left({}^1{F}_{11}\right)}^2+{\left({}^1{F}_{21}\right)}^2}$(4-4)

The angles between two perpendicular directions at undeformed configuration and two directions in defomed configuration at the same considered point can be expressed by:

$\theta =\arccos \left(\frac{\left[\begin{array}{cc}1& 0\end{array}\right]\left[{{}^{0}F}^T{}^{0}F\right]\mathrm{}\left[\begin{array}{c}0\\ 1\end{array}\right]}{{}^0{\lambda }_{11}\mathrm{}{}^0{\lambda }_{22}}\right)=\arccos \left(\frac{{}^0{F}_{11}{}^0{F}_{12}+{}^0{F}_{21}{}^0{F}_{22}}{\sqrt{{\left({}^0{F}_{11}\right)}^2+{\left({}^0{F}_{21}\right)}^2}}\right)$(4-5)

$\Theta =\arccos \left(\frac{\left[\begin{array}{cc}1& 0\end{array}\right]\left[{{}^{1}F}^T{}^{1}F\right]\mathrm{}\left[\begin{array}{c}0\\ 1\end{array}\right]}{{}^1{\lambda }_{11}\mathrm{}{}^1{\lambda }_{22}}\right)=\arccos \left(\frac{{}^1{F}_{11}{}^1{F}_{12}+{}^1{F}_{21}{}^1{F}_{22}}{\sqrt{{\left({}^1{F}_{11}\right)}^2+{\left({}^1{F}_{21}\right)}^2}}\right)$(4-6)

in which ${}^0{\lambda }_{22}={}^1{\lambda }_{22}=1$.

Finally, the components of the strain tensor referring to the respective normal and shear strains can be expressed by:

${\epsilon }_{11}={\lambda }_{11}-1=\frac{{}^1{\lambda }_{11}}{{}^0{\lambda }_{11}}-1=\frac{\sqrt{{\left({}^1{F}_{11}\right)}^2+{\left({}^1{F}_{21}\right)}^2}}{\sqrt{{\left({}^0{F}_{11}\right)}^2+{\left({}^0{F}_{21}\right)}^2}}-1$(4-7)

${\epsilon }_{22}={\lambda }_{22}-1=\frac{{}^1{\lambda }_{22}}{{}^0{\lambda }_{22}}-1=0$(4-8)

${\epsilon }_{12}=\frac{1}{2}\left[\arccos \left(\frac{{}^0{F}_{11}{}^0{F}_{12}+{}^0{F}_{21}{}^0{F}_{22}}{\sqrt{{\left({}^0{F}_{11}\right)}^2+{\left({}^0{F}_{21}\right)}^2}}\right)-\arccos \left(\frac{{}^1{F}_{11}{}^1{F}_{12}+{}^1{F}_{21}{}^1{F}_{22}}{\sqrt{{\left({}^1{F}_{11}\right)}^2+{\left({}^1{F}_{21}\right)}^2}}\right)\right]$(4-9)

Considering the stress-strain relation deduced for Boltzmann rheological model, expressed by Equation (3-8), and the appropriate strain components for plane frame finite element with Reissner kinematics, the total strain energy for the viscoelastic behavior can be expressed by:

$U=\mathrm{}{\int }_V^{}\left[{\int }_{\epsilon }^{}\left(\frac{E_1{E}_2}{E_1+E_2}{\epsilon }_{11}+\frac{\eta E_1}{E_1+E_2}{\dot{\epsilon }}_{11}-\frac{\eta }{E_1+E_2}{\dot{\sigma }}_{11}\right)\mathrm{}d\epsilon +\right.$

$+\left.{\int }_{\epsilon }^{}\left(2\frac{E_1{E}_2}{E_1+E_2}\stackrel{-}{\mu }{\epsilon }_{12}+2\frac{\eta E_1}{E_1+E_2}\stackrel{-}{\mu }{\dot{\epsilon }}_{12}-\frac{\eta }{E_1+E_2}{\dot{\sigma }}_{12}\right)\mathrm{}d\epsilon \right]dV$(4-10)

The chain rule can be used to exchange integration variable and the total strain energy can be rewritten as:

$U=\mathrm{}{\int }_V^{}\left[{\int }_X^{}\left(\frac{E_1{E}_2}{E_1+E_2}{\epsilon }_{11}+\frac{\eta E_1}{E_1+E_2}{\dot{\epsilon }}_{11}-\frac{\eta }{E_1+E_2}{\dot{\sigma }}_{11}\right){{\epsilon }_{11},}_{q}dX+\right.$

$+\left.{\int }_X^{}\left(2\frac{E_1{E}_2}{E_1+E_2}\stackrel{-}{\mu }{\epsilon }_{12}+2\frac{\eta E_1}{E_1+E_2}\stackrel{-}{\mu }{\dot{\epsilon }}_{12}-\frac{\eta }{E_1+E_2}{\dot{\sigma }}_{12}\right){{\epsilon }_{12},}_{q}dX\right]dV\mathrm{}\mathrm{}$(4-11)

From Equation (4-11), it is possible to note that the total strain energy can be expressed by two uncoupled parts. One of these parts is associated with the normal stress effects and the other is associated with the shear stress effects. In 2D elasticity problems the uncoupling between volumetric and deviatoric parts is usual (Oliveira et al., 2018). Thus, there is a schematic interpretation for the adoption of uncoupled models for each of these effects. In Figure 3, two uncoupled Boltzmann models are considered, one for the normal stress effects and the other for the shear stress effects. Similarly, a schematic interpretation can be obtained by adopting other rheological models, and it is even possible to consider the combination of different rheological models, one for the normal stress effects and the other for the shear stress effects.

Figure 3
Schematic interpretation of uncoupled models.

Due to the characteristic viscoelastic behavior of the Boltzmann model, it is important to note that the stress-strain relation, expressed by Equation (3-8), is time-dependent, as well as the total strain energy, expressed by Equation (4-11), and its derivatives. This time dependence is represented by the strain rates and stress rates. Thus, an adequate procedure is necessary to evaluate them over time.

In the present paper, these rates are evaluated using the backward Euler method with an equidistant partition of a time interval (a fixed time step $∆t$) as follows:

$\dot{\epsilon }=\mathrm{}\frac{{\epsilon }^s-{\epsilon }^{s-1}}{∆t}$ and $\dot{\sigma }=\mathrm{}\frac{{\sigma }^s-{\sigma }^{s-1}}{∆t}$(4-12)

in which indexes “s” e “s-1” represent current and previous instant of time, respectively, while $∆t$ represents an imposed time step.

Therefore, the stress-strain relation for Boltzmann model can be expressed by:

${\sigma }_{ij}^s={\frac{E_1{E}_2}{E_1+E_2}\stackrel{-}{C}}_{ijkl}{\epsilon }_{kl}^s+{\frac{\eta E_1}{E_1+E_2}\stackrel{-}{C}}_{ijkl}\frac{{\epsilon }_{kl}^s-{\epsilon }_{kl}^{s-1}}{∆t}-\frac{\eta }{E_1+E_2}\frac{{\sigma }_{kl}^s-{\sigma }_{kl}^{s-1}}{∆t}$(4-13)

and total strain energy is represented by:

$U^s=\mathrm{}{\int }_V^{}\left[{\int }_X^{}\left(\frac{E_1{E}_2}{E_1+E_2}{\epsilon }_{11}^s+\frac{\eta E_1}{E_1+E_2}\frac{{\epsilon }_{11}^s-{\epsilon }_{11}^{s-1}}{∆t}-\frac{\eta }{E_1+E_2}\frac{{\sigma }_{11}^s-{\sigma }_{11}^{s-1}}{∆t}\right){{\epsilon }_{11}^s,}_{q}dX+\right.$

$+\left.{\int }_X^{}\left(2\frac{E_1{E}_2}{E_1+E_2}\stackrel{-}{\mu }{\epsilon }_{12}^s+2\frac{\eta E_1}{E_1+E_2}\stackrel{-}{\mu }\frac{{\epsilon }_{12}^s-{\epsilon }_{12}^{s-1}}{∆t}-\frac{\eta }{E_1+E_2}\frac{{\sigma }_{12}^s-{\sigma }_{12}^{s-1}}{∆t}\right){{\epsilon }_{12}^s,}_{q}dX\right]dV$(4-14)

Finally, considering the first and second derivatives of the total strain energy, it is possible to proceed to the Newton-Raphson Iterative Method and apply the Principle of Minimal Total Potential Energy, as described in section 2. Thus, the presented formulation can be used to analyze structures that are discretized by plane frame elements with Reissner kinematics and that exhibit viscoelastic behavior characteristic of the Boltzmann model.

5 ANALYSES AND RESULTS

In this section a parametric analysis and an example of a beam with several height/span ratios are presented in order to demonstrate the consistency of the positional formulation using Reissner kinematics and its capacity to describe the viscoelastic behavior considering the shear effects. In addition, three examples are presented to assess the agreement of the numerical results obtained in relation to the numerical results available in the literature. The three parameter Boltzmann rheological model was used to model the viscoelastic solids. As the model is one of the simplest suitable for solids, the main mechanical parameters (except temperature and humidity) were chosen to be assess in the parametric analyses. Moreover, spatial discretization and time steps influence were included in the assessment. No inertial effects were considered and a numerical tolerance of 10^-8 was adopted for calculations (absolute in terms of nodal positions).

5.1 Parametric analysis

This analysis aims to evaluate the influence of physical parameters of the Boltzmann model and numerical parameters of the simulation in the description of the viscoelastic behavior. For this purpose, the numerical results obtained in the analysis of a supported beam with centered vertical load are evaluated, as shown in Figure 4.

Figure 4
Supported beam with vertical load applied on mid-span.

The beam has length (L) equal to 2 m and rectangular cross-section with height equal to 0.40 m and width equal to 0.01 m. The vertical load (P) is equal to 1000 kN. The rheological model considered to assess the viscoelastic behavior is the Boltzmann model with longitudinal elastic modulus (E₁) equal to 100 GPa, longitudinal elastic modulus (E₂) equal to 400 GPa, modulus of viscosity (η) equal to 5000 GPa·s, and Poisson coefficient equal to 0.3. Regarding spatial discretization, 10 finite elements are adopted, with 10 Gauss points along height and 10 Gauss points along length. Regarding temporal discretization, 10 time steps equal to 10 s are considered. In each of the following subsections, the variation of each of these parameters is evaluated.

The analyses presented in the following are based on the results of vertical displacement at mid-span over time when the physical parameters of the Boltzmann model and the numerical parameters of the simulation are varied.

The influence of the longitudinal elastic modulus E₁ is shown in Figure 5(a). It is possible to observe a non-linear relation between the displacements and the longitudinal elastic modulus E₁. The influence of the longitudinal elastic modulus E₂ is shown in Figure 5(b). Again, it is possible to observe a non-linear relation between the displacements and the longitudinal elastic modulus E₂. From the results presented in Figures 5(a) and 5(b), it can also be observed that the same final viscoelastic displacements are obtained with the same variations performed in the longitudinal elastic modulus E₁ or in the longitudinal elastic modulus E₂. This result is consistent since the viscoelastic response for a sufficiently long time (final viscoelastic response) of the Boltzmann model is determined by the equivalent longitudinal elastic modulus given by the series association between the spring with longitudinal elastic modulus E₁ and the spring with longitudinal elastic modulus E₂.

Figure 5
Vertical displacement at mid-span over time, considering several: (a) moduli of elasticity E₁; (b) moduli of elasticity E₂; (c) moduli of viscosity η; (d) Poisson coefficients ν; (e) time steps ∆t.

The influence of the modulus of viscosity η in the Boltzmann model response is shown in Figure 5(c). It is possible to observe a non-linear relation between the displacements and the modulus of viscosity η. The influence of the Poisson coefficient ν in the Boltzmann model response is shown in Figure 5(d). Now, it is possible to observe a linear relation between the displacements and the Poisson coefficient ν. It is important to note that, despite the physical inconsistency in the adoption of a Poisson coefficient equal to 0.5, a numerical inconsistency is not observed. This fact is due to the consideration of unidimensional frame elements with uncoupled between normal stress effects and shear stress effects. Therefore, the Poisson coefficient is considered only for assessing the parameters associated with shear through the second Lamé parameters, while the first Lamé parameter in the isotropy relation is not considered in the formulation.

Figure 5(e) presentes the evaluation of the adopted time step influence. It is possible to observe a dependence on the displacement results regarding temporal discretization. It is also possible to observe convergence of the results with the temporal discretization refinement. In addition, the instantaneous elastic response and the final viscoelastic response do not change with the size variation of the adopted time interval. These results are consistent with the approach adopted for assessing strain rates, based on the backward Euler method, and are in accordance with the results obtained in the works of several authors such as Mesquita and Coda (2007) and Oliveira and Leonel (2017).

To assess the influence of spatial discretization, the results obtained by adopting 2, 10, 16 and 64 finite elements are considered. In this case, by refining the mesh, it is possible to notice a small variation in the numerical results, observing a trend of convergence as the number of finite elements increases. Between the numerical results with two and ten elements the difference is of the order of 10^-5 m, between the numerical results with ten and sixteen elements the difference is of the order of 10^-9 m, and between the numerical results with sixteen and sixty-four elements the difference is on the order of 10^-10 m.

5.2 Short beam bending analysis

In this example, a supported beam with centered vertical load is analyzed, as shown in Figure 6, considering several height measurements of the cross-section. Thus, from the comparison between the results obtained using the formulation developed based on Reissner kinematics and the results obtained with the formulation presented in Becho et al. (2015) based on the Bernoulli-Euler kinematics, it is possible to evaluate the shear effects on viscoelastic behavior.

Figure 6
Supported beam with vertical load applied on mid-span.

The beam shown in Figure 6 has length (L) equal to 2 m and a rectangular cross-section of width (b) equal to 0.10 m. As for the height (h) of the cross-section, five cases with dimensions equal to 0.10 m, 0.20 m, 0.30 m, 0.40 m, and 0.50 m are considered. In addition, in order to maintain the same maximum stress level (150 MPa) in each case, applied load intensities (P) equal to 50 kN, 200 kN, 450 kN, 800 kN, and 1250 kN are adopted, respectively. This approach is adopted only to keep the displacement results in the same level of magnitude.

Regarding spatial discretization, 10 finite elements are adopted, with 10 Gauss points along the height and along the length. Regarding temporal discretization, 20 time steps equal to 5 s are adopted. In addition, the rheological model considered to assess the viscoelastic behavior is the Boltzmann model with longitudinal elastic modulus (E₁) equal to 100 GPa, longitudinal elastic modulus (E₂) equal to 400 GPa, modulus of viscosity (η) equal to 5000 GPa·s, and Poisson coefficient equal to 0.3.

The numerical results of vertical displacement at mid-span over time, for each of the cases of height of the cross-section and for each of the adopted kinematics, are presented in Figure 7 and Table 1. Additionally, Figure 7 and Table 2 show the analytical results for each of the cases of height of the cross-section. Such analytical results can be obtained from the Equations (5-1) and (5-2) without shear effects and with shear effects.

$w\left(t\right)=\frac{P\mathrm{}{L}^3}{4\mathrm{}b{h}^3}\left(\frac{1}{E_1}+\frac{1}{E_2}\left(1-e^{\frac{-E_2}{\eta }\mathrm{}t}\right)\right)$(5-1)

$w\left(t\right)=\frac{P\mathrm{}{L}^3}{4\mathrm{}b{h}^3}\left(1+2\frac{\left(1+\nu \right)}{k}{\left(\frac{h}{L}\right)}^2\right)\left(\frac{1}{E_1}+\frac{1}{E_2}\left(1-e^{\frac{-E_2}{\eta }\mathrm{}t}\right)\right)$(5-2)

in which w(t) represents vertical displacement at mid-span and k represents the form factor (determined as 5/6 for rectangular cross-sections). These analytical expressions are obtained based on fluency functions expressed in Prony series as shown in Chen (1995).

Figure 7
Vertical displacement at mid-span over time, considering several height/span ratios.

Table 1
Numerical results with both kinematics.

Height ($h$) [m]	Ratio ($h/L$)	Applied load [kN]	Instantaneous elastic displacement [mm]		Final viscoelastic displacement [mm]		Relative increase *
			Bernoulli-Euler	Reissner	Bernoulli-Euler	Reissner
0.10	5 %	50	9.9990	10.0640	12.4980	12.5792	25.81 %
0.20	10 %	200	4.9999	5.1299	6.2498	6.4122	28.25 %
0.30	15 %	450	3.3333	3.5283	4.1666	4.4103	32.31 %
0.40	20 %	800	2.5000	2.7600	3.1250	3.4500	38.00 %
0.50	25 %	1250	2.0000	2.3250	2.5000	2.9062	45.31 %

• ^*
  Relative increase of the final viscoelastic displacements (t = 100 s) using the Reissner kinematics compared with the instantaneous elastic displacements (t = 0) using the Bernoulli-Euler kinematics (simultaneous effects of the shear and the viscoelastic behavior).

Table 2
Analytical results considering and ignoring the shear effects.

Height ($h$) [m]	Ratio ($h/L$)	Applied load [kN]	Instantaneous elastic displacement [mm]		Final viscoelastic displacement [mm]		Relative increase *
			With shear	Without shear	With shear	Without Shear
0.10	5 %	50	10.0000	10.0780	12.4992	12.5967	25.97 %
0.20	10 %	200	5.0000	5.1560	6.2496	6.4446	28.89 %
0.30	15 %	450	3.3333	3.5673	4.1664	4.4589	33.77 %
0.40	20 %	800	2.5000	2.8120	3.1248	3.5148	40.59 %
0.50	25 %	1250	2.0000	2.3900	2.4998	2.9873	49.37 %

• ^*
  Relative increase of the final viscoelastic displacements (t = 100 s) using the Reissner kinematics compared with the instantaneous elastic displacements (t = 0) using the Bernoulli-Euler kinematics (simultaneous effects of the shear and the viscoelastic behavior).

From Figure 7, comparing the respective results obtained with the Reissner and Bernoulli-Euler kinematics, it is possible to observe that with the increase in the height/span ratio, the shear effects are more pronounced. In addition, by comparing the numerical results and the analytical results, it is possible to verify the agreement in each of the kinematics. For h/L ratio around 5% there is practically no influence of shear and curves for both kinematics begins at the same point. As the analysis evolves, slight differences between displacements are noted. Thus, is possible to state that there is a very small influence of viscoelastic behavior too between kinematics. As h/L ratio increases shear effects are increased and viscoelastic curves are almost parallel.

Finally, by evaluating the relative percentage increase between the final viscoelastic displacements (t = 100 s) obtained with the Reissner kinematics and the instantaneous elastic displacements (t = 0) obtained with the Bernoulli-Euler kinematics, it is possible to obtain the results shown in the last column of Table 1. These results represent the simultaneous effects of the shear and the viscoelastic behavior. Similar results can be obtained considering the relative percentage increase between the viscoelastic displacements at any time instant, obtained with the Reissner kinematics, and the instantaneous elastic displacements, obtained by the Bernoulli-Euler kinematics, for each height/span ratio. The results obtained in this way allow to evaluate the simultaneous effects of the shear and the viscoelastic behavior at any time instant and can be interpreted graphically as shown in Figure 8.

Figure 8
Simultaneous effects of the shear and the viscoelastic behavior at several time instants as a function of the height/span ratio.

From the results obtained in this example, it is possible to note the more pronounced influence of the shear effects in short beams, demonstrating the importance of adopting a kinematics that takes into account such effects. In addition, it is possible to observe that the simultaneous effects of the shear and the viscoelastic behavior can be relevant, depending on the geometric and physical properties, and should not be neglected.

5.3 Tensioned bar analysis

A classical example is analyzed here. This example deals with a tensioned bar, as shown in Figure 9, used by several authors to verify the viscoelastic models because it is a typical case of a plane stress state. The beam has length (L) equal to 800 mm and rectangular cross-section with height (h) equal to 100 mm and unitary width (b). The traction (P) is equal to 0.005 kN/mm².

Figure 9
Tensioned bar.

Regarding spatial discretization, 10 finite elements are adopted, with 10 Gauss points along the height and along the length. Regarding temporal discretization, several time steps are adopted (1 day, 5 days, 10 days, 25 days and 50 days). These distinct time step values are adopted in order to assess the influence of the refinement of temporal discretization. In addition, the rheological model considered to assess the viscoelastic behavior is the Boltzmann model with longitudinal elastic modulus (E₁) equal to 22,5757 kN/mm², longitudinal elastic modulus (E₂) equal to 11,0 kN/mm², modulus of viscosity (η) equal to 500,0 kN/mm²·day, and Poisson coefficient equal to 0.

Numerical results for axial displacements at free end are presented in Figure 10. The presented analytical results and numerical results are also available in the literature. The numerical results available in Mesquita and Coda (2007) are achieved with a mesh of 8x4 cubic boundary elements. The analytical results are obtained from the solution of the one-dimensional differential equation, which governs the behavior of the Boltzmann model, expressed by:

Figure 10
Axial displacements at free end over time, considering several time steps ∆t.

$\epsilon \left(t\right)=\left(\frac{P}{E_1}-\frac{E_1+E_2}{E_2{E}_1}P\right)\mathrm{}{e}^{-\frac{E_2}{\eta }\mathrm{}t}+\frac{E_1+E_2}{E_2{E}_1}P$(5-3)

From Figure 10, it is possible to observe that the reduction in the time step results in the approximation of the numerical results in relation to the analytical results, presenting a tendency of convergence of the numerical results with the temporal discretization refinement. In addition, this observed behavior is consistent and in accordance with the results available in Mesquita and Coda (2007), obtained using a different formulation than the one presented in this paper.

Finally, Figure 11 shows the numerical results for the viscoelastic deformation-recovery process. In this case, at the instant of time equal to 200 days, the traction (P) is removed, checking the process of recovery of viscoelastic displacements, considering the time step equal to 1 day. From Figure 11, it is possible to note that the results obtained with the developed formulation are consistent, showing agreement with the numerical results available in Mesquita and Coda (2007).

Figure 11
Viscoelastic deformation-recovery process of a tensioned bar.

5.4 Cantilever beam analysis

This example deals with a cantilever beam loaded with a shear load applied on free end, as shown in Figure 12.

Figure 12
Cantilever beam.

Geometry, rheological parameters, and spatial discretization are the same as the ones presented in section 5.3. The shear load is equal to 0.005 kN/mm², applied at initial time, and removed after 453 days. Further analysis continues up to 267 days without applied load in order to observe the viscoelastic recovery process. Regarding temporal discretization, 720 time steps equal to 1 day are considered.

Numerical results for vertical displacements in the middle point of the right face are presented in Figure 13, for viscoelastic deformation-recovery process. The numerical results available in Panagiotopoulos et al. (2014) made with a BEM are also presented here. In addition, the obtained numerical results using the formulation developed in Becho et al. (2015) are presented. It is important to note that in Becho et al. (2015) the formulation adopts the Bernoulli-Euler kinematics, therefore, it does not take into account the shear effects.

Figure 13
Viscoelastic deformation-recovery process of a cantilever beam.

From Figure 13, it is possible to observe the agreement with the numerical results available in the literature. In addition, it is possible to observe that adopting Reissner kinematics the displacements are slightly higher than those obtained by adopting Bernoulli-Euler kinematics. However, this difference is not very pronounced due to the height/span ratio of the beam being only 12.5%.

5.5 Supported beam with variable cross-section under bending

In this example, a supported beam with uniformly distributed load is analyzed, as shown in Figure 14, considering variable cross-section.

Figure 14
Supported beam with variable cross-section.

The beam has length (L) equal to 10 m, rectangular cross-section with width equal to 0.01 m and height equal to h₀ at extremities and equal to h₁ at mid-span (with linear variation along the length). The uniformly distributed load (P) is equal to 10 N/m.

Regarding spatial discretization, 40 finite elements are adopted, with 10 Gauss points both along the height and the length. Regarding temporal discretization, 500 time steps equal to 0.04 s are adopted. In addition, the rheological model considered to assess the viscoelastic behavior is the Boltzmann model with longitudinal elastic modulus (E₁) equal to 98,00 MN/m², longitudinal elastic modulus (E₂) equal to 24.50 MN/m², modulus of viscosity (η) equal to 274,4 MN/m²·s, and Poisson coefficient equal to 0.3.

Numerical results for vertical displacements at mid-span are presented in Figure 15. These results refer to the adoption of several variable cross-sections through relations between height h₀ (at the ends) and height h₁ (in the middle of the span). The numerical results available in Aköz and Kadioǧlu (1999) made with a classical FEM considering Timoshenko kinematics are also presented here. In addition, the obtained numerical results using the formulation developed in Becho et al. (2015) considering Bernoulli-Euler kinematics are presented.

Figure 15
Vertical displacement at mid-span over time, considering several variable cross-sections.

Finally, Figure 16 presents the results obtained considering the constant cross-sections with height equal to 0.7 m. In addition, the analytical results available in the literature are presented. The analytical results are obtained by the expressions presented in Chen (1995) and rewritten, without considering the shear effects and considering the shear effects, respectively, as follows:

$w\left(t\right)=\frac{5P{L}^4}{32\mathrm{}b{h}^3}\left(\frac{1}{E_1}+\frac{1}{E_2}\left(1-e^{\frac{-E_2}{\eta }\mathrm{}t}\right)\right)$(5-4)

$w\left(t\right)=\frac{5P{L}^4}{32\mathrm{}b{h}^3}\left(1+1.6\frac{\left(1+\nu \right)}{k}{\left(\frac{h}{L}\right)}^2\right)\left(\frac{1}{E_1}+\frac{1}{E_2}\left(1-e^{\frac{-E_2}{\eta }\mathrm{}t}\right)\right)$(5-5)

in which w(t) represents vertical displacement at mid-span and k represents the form factor (determined as 5/6 for rectangular cross-sections).

Figure 16
Vertical displacement at mid-span over time, considering constant cross-sections (h₀=h₁=0.7m).

From Figure 15, it is possible to observe the agreement between the results obtained and the results available in Aköz and Kadioǧlu (1999). In addition, it is possible to observe the difference between the results obtained with the Reissner kinematics and the results obtained with the Bernoulli-Euler kinematics. This difference, due to the contribution of the shear effects, is seen more clearly in Figure 16. From Figure 16, it is possible to observe the agreement between the numerical results obtained and the analytical results.

6 CONCLUSION

In the present paper, a numerical formulation based on the positional formulation of the Finite Element Method is developed to describe the viscoelastic mechanical behavior in structures discretized by plane frame elements, considering the shear effects. To assess the effects of the viscoelastic behavior, the stress-strain relation deduced from the Boltzmann rheological model is used. While for numerical evaluation of the contributions of the shear effects, the cross-section rotations are parameterized independently using the Reissner kinematics. Such an approach is considered original for the description of viscoelastic behavior.

Similar developments can be made considering other rheological models, such as the Kelvin-Voigt model and the Zener model. To do this, it is enough to develop the appropriate stress-strain relations and, based on these, to particularize the expressions to evaluate the strain energy and its derivatives.

Based on the developments presented and the results obtained, the formulation developed is considered relatively simple and capable of considering the relevant shear effects on viscoelastic behavior. This simplicity is attributed to the compact and intuitively obtained expressions based on energy concepts and the physical behavior of rheological models. The capacity to consider the shear effects is verified by comparing the results obtained by adopting the Reissner kinematics with the results obtained by adopting the Bernoulli-Euler kinematics. In addition, it is possible to observe the agreement of the numerical results obtained in relation to the analytical and numerical results available in the literature, when the shear effects are taken into account.

Finally, from the results obtained, it is possible to verify the relevance of the shear effects on the viscoelastic behavior, justifying the importance of using formulations that take these effects into account. These results are in accordance with several experimental studies available in the literature, reinforcing that the shear effects on viscoelastic behavior cannot be neglected, depending on the physical characteristics of the material, the geometric characteristics of the structural element and the imposed service conditions, as demonstrated in the experimental works of Bank and Mosallam (1992), Mottram (1993), Abdel-Magid et al. (2003), Shao and Shanmugam (2004), Sá et al. (2011a), and Garrido and Correia (2018).

Acknowledgments

The authors would like to thank Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) for their research financial support.

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

References

• Abdel-Magid, B., Lopez-Anido, R., Smith, G., Trofka, S., (2003). Flexure creep properties of E-glass reinforced polymers. Composite Structures 62:247-253. doi:10.1016/j.compstruct.2003.09.022
  » https://doi.org/10.1016/j.compstruct.2003.09.022
• Aköz, Y., Kadioǧlu, F., (1999). The mixed finite element method for the quasi‐static and dynamic analysis of viscoelastic Timoshenko beams. International Journal for Numerical Methods in Engineering 44:1909-1932. doi:10.1002/(SICI)1097-0207(19990430)44:12<1909::AID-NME573>3.0.CO;2-P
  » https://doi.org/10.1002/(SICI)1097-0207(19990430)44:12<1909::AID-NME573>3.0.CO;2-P
• Alachek, I., Reboul, N., Jurkiewiez, B., (2019). Experimental and finite element analysis of the long-term behaviour of GFRP-concrete hybrid beams fabricated using adhesive bonding. Composite Structures 207, 148-165. doi:10.1016/j.compstruct.2018.09.013
  » https://doi.org/10.1016/j.compstruct.2018.09.013
• Argyris, J., Doltsinis, I. St.; da Silva, V.D., (1991). Constitutive modelling and computation of non-linear viscoelastic solids. Part I: Rheological models and numerical integration techniques. Computer Methods in Applied Mechanics and Engineering, 88: 135-163. doi: 10.1016/0045-7825(91)90252-2.
  » https://doi.org/10.1016/0045-7825(91)90252-2
• Bank, L.C., Mosallam, A.S., (1992). Creep and failure of a full-size fiber-reinforced plastic pultruded frame. Composites Engineering, 2:213-227. doi:10.1016/0961-9526(92)90005-Q
  » https://doi.org/10.1016/0961-9526(92)90005-Q
• Becho, J.S., Greco, M., (2021). Influence of retardation time on viscoelastic behavior of plane beams modeled by the nonlinear positional formulation of the finite element method. Mathematical Problems in Engineering, 2021: Article ID-2684127. doi: 10.1155/2021/2684127
  » https://doi.org/10.1155/2021/2684127
• Becho, J.S., Barros, F.B., Greco, M., (2015). Positional formulation for creep mechanical behavior of beams and framed structures. Science & Engineering Journal 24:31-41. doi:10.14393/19834071.2015.30300
  » https://doi.org/10.14393/19834071.2015.30300
• Carniel, T.A., Muñoz-Rojas, P.A., Vaz, 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 Modeling 39:1725-1739. doi:10.1016/j.apm.2014.09.036
  » https://doi.org/10.1016/j.apm.2014.09.036
• Carrazedo, R., Coda, H.B., (2017). Triangular based prismatic finite element for the analysis of orthotropic laminated beams, plates and shells. Composite Structures 168:234-246. doi: 10.1016/j.compstruct.2017.02.027
  » https://doi.org/10.1016/j.compstruct.2017.02.027
• Cavalcante, J.P.B., Maciel, D.N., Greco, M., (2017). Impact response of flying objects modeled by Positional Finite Element Method. International Journal of Structural Stability and Dynamics 18:1850076-1. doi:10.1142/S0219455418500761
  » https://doi.org/10.1142/S0219455418500761
• Chen, T.M., (1995). The hybrid Laplace transform/finite element method applied to the quasi‐static and dynamic analysis of viscoelastic Timoshenko beams. International Journal for Numerical Methods in Engineering 38:509-522. doi:10.1002/nme.1620380310
  » https://doi.org/10.1002/nme.1620380310
• Christensen, R.M., (2003). Theory of viscoelasticity. Dover Publications, New York.
• Coda, H.B., Greco, M., (2004). A simple FEM formulation for large deflection 2D frame analysis based on position description. Computer Methods in Applied Mechanics and Engineering 193:3541-3557. doi:10.1016/j.cma.2004.01.005
  » https://doi.org/10.1016/j.cma.2004.01.005
• Coda, H.B., Paccola, R.R., (2007). An alternative positional FEM formulation for geometrically nonlinear analysis of shells-curved triangular isoparametric elements. Computational Mechanics 40:185-200. doi:10.1007/s00466-006-0094-1
  » https://doi.org/10.1007/s00466-006-0094-1
• Crisfield, M.A., (1991). Non-linear finite element analysis of solids and structures. v.1. John Wiley & Sons, England.
• Dastjerdi, S., Akgöz, B., Civalek, Ö., (2020). On the effect of viscoelasticity on behavior of gyroscopes. International Journal of Engineering Science 149:103236. doi:10.1016/j.ijengsci.2020.103236
  » https://doi.org/10.1016/j.ijengsci.2020.103236
• Demir, Ç., Civalek, Ö., (2017). On the analysis of microbeams. International Journal of Engineering Science 121(1):14-33. doi:10.1016/j.ijengsci.2017.08.016
  » https://doi.org/10.1016/j.ijengsci.2017.08.016
• Dym, C.L., Shames, I.H., (2013). Solid Mechanics: A variational approach. Springer-Verlag, New York.
• Ebrahimi, F., Barati, M.R., Civalek, Ö. (2020). Application of Chebyshev–Ritz method for static stability and vibration analysis of nonlocal microstructure-dependent nanostructures. Engineering with Computers 36(3):953-964. doi:10.1007/s00366-019-00742-z
  » https://doi.org/10.1007/s00366-019-00742-z
• Fernandes, W.L., Vasconcellos, D.B., Greco, M., (2018). Dynamic instability in shallow arches under transversal forces and plane frames with semirigid connections. Mathematical Problems in Engineering, Article ID 1985907. doi:10.1155/2018/1985907
  » https://doi.org/10.1155/2018/1985907
• Findley, W.N., Lai, J.S., Onaran, K., (1989). Creep and relaxation of nonlinear viscoelastic materials. Dover Publications, New York.
• Galucio, A.C., Deü, J.F., Ohayon, R., (2004). Finite element formulation of viscoelastic sandwich beams using fractional derivative operators. Computational mechanics 33(4):282-291. doi:10.1007/s00466-003-0529-x
  » https://doi.org/10.1007/s00466-003-0529-x
• Gand, A.K., Chan, T.M., Mottram, J.T. (2013). Civil and structural engineering applications, recent trends, research and developments on pultruded fiber reinforced polymer closed sections: a review. Frontiers of Structural and Civil Engineering 7(3):227-244. doi:10.1007/s11709-013-0216-8
  » https://doi.org/10.1007/s11709-013-0216-8
• Garrido, M., Correia, J.R., (2018). Elastic and viscoelastic behaviour of sandwich panels with glass-fibre reinforced polymer faces and polyethylene terephthalate foam core. Journal of Sandwich Structures & Materials 20:399-424. doi:10.1177/1099636216657388
  » https://doi.org/10.1177/1099636216657388
• Greco, M., Coda, H.B., (2006). Positional FEM formulation for flexible multi-body dynamic analysis. Journal of Sound and Vibration 290:1141-1174. doi:10.1016/j.jsv.2005.05.018
  » https://doi.org/10.1016/j.jsv.2005.05.018
• Greco, M., Ferreira, I.P., (2009). Logarithmic strain measure applied do the nonlinear positional formulation for space truss analysis. Finite Elements in Analysis and Design 45:632-639. doi:10.1016/j.finel.2009.05.005
  » https://doi.org/10.1016/j.finel.2009.05.005
• Jalaei, M.H., Civalek, Ӧ., (2019). On dynamic instability of magnetically embedded viscoelastic porous FG nanobeam. International Journal of Engineering Science 143:14-32. doi:10.1016/j.ijengsci.2019.06.013
  » https://doi.org/10.1016/j.ijengsci.2019.06.013
• Kassner, M.E., Pérez-Prado, M., (2004). Fundamentals of creep in metals and alloys. Elsevier, New York.
• Kühl, A., Muñoz-Rojas, P.A., (2019). Application of a master curve and the modified superposition principle for modeling creep and loading rate effects at small strains in high-density polyethylene. Mechanics of Advanced Materials and Structures 1:1-12. doi:10.1080/15376494.2019.1605008
  » https://doi.org/10.1080/15376494.2019.1605008
• Li L., Hu, Y., Ling, L., (2016). Wave propagation in viscoelastic single-walled carbon nanotubes with surface effect under magnetic field based on nonlocal strain gradient theory. Physica E: Low-dimensional Systems and Nanostructures 75:118-124. doi:10.1016/j.physe.2015.09.028
  » https://doi.org/10.1016/j.physe.2015.09.028
• Liu, H., Polak, M., Penlidis, A., (2008). A practical approach to modeling time-dependent nonlinear creep behavior of polyethylene for Structural applications. Polymer Engineering and Science 48:159-167. doi:10.1002/pen.20942
  » https://doi.org/10.1002/pen.20942
• Maciel, D.N., Coda, H.B., (2008). Positional Finite Element Methodology for geometrically nonlinear analysis of 2D frames. Minerva Pesquisa e Tecnologia 5:73-83.
• Maciel, D.N., Greco, M., Coda, H.B., (2004). Positional description for nonlinear 2D frames analysis by FEM with Reissner kinematics. In: XXV Iberian Latin American Congress on Computational Methods in Engineering, Recife, Brazil.
• Marques, S.P.C., Creus, G.J., (2012). Computational viscoelasticity. Springer, Heidelberg.
• Mesquita, A.D., Coda, H.B., (2001). An alternative time integration procedure for Boltzmann viscoelasticity: a BEM approach. Computers & Structures 79:1487-1496. doi:10.1016/S0045-7949(01)00042-6
  » https://doi.org/10.1016/S0045-7949(01)00042-6
• Mesquita, A.D., Coda, H.B., (2002). Alternative Kelvin viscoelastic procedure for finite elements. Applied Mathematical Modelling, 26: 501–516. doi: 10.1016/S0307-904X(01)00048-8.
  » https://doi.org/10.1016/S0307-904X(01)00048-8
• Mesquita, A.D., Coda, H.B., (2003). A simple Kelvin and Boltzmann viscoelastic analysis of three-dimensional solids by the boundary element method. Engineering Analysis with Boundary Elements 27:885-895. doi:10.1016/S0955-7997(03)00060-2
  » https://doi.org/10.1016/S0955-7997(03)00060-2
• Mesquita, A.D., Coda, H.B., (2007). A boundary element methodology for viscoelastic analysis: Part I with cells. Applied Mathematical Modeling 31:1149-1170. doi:10.1016/j.apm.2006.04.006
  » https://doi.org/10.1016/j.apm.2006.04.006
• Mottram, J.T., (1993). Short- and long-term structural properties of pultruded beam assemblies fabricated using adhesive bonding. Composite Structures 25:387-395. doi:10.1016/0263-8223(93)90186-T
  » https://doi.org/10.1016/0263-8223(93)90186-T
• Oliveira, F.M., Greco, M., (2014). Nonlinear dynamic analysis of beams with layered cross-sections under moving masses. Journal of the Brazilian Society of Mechanical Sciences and Engineering 37:451-462. doi:10.1007/s40430-014-0184-5
  » https://doi.org/10.1007/s40430-014-0184-5
• Oliveira, H.L., Leonel, E.D., (2017). A BEM formulation applied in the mechanical material modeling of viscoelastic cracked structures. International Journal of Advanced Structural Engineering 9:1-12. doi:10.1007/s40091-016-0144-7
  » https://doi.org/10.1007/s40091-016-0144-7
• Oliveira, H.L., Chateauneuf, A., Leonel, E.D., (2018). Probabilistic mechanical modelling of concrete creep based on the boundary element method Advances in Structural Engineering, 2018: 1–12. doi: 10.1177/1369433218788088.
  » https://doi.org/10.1177/1369433218788088
• Panagiotopoulos, C.G., Mantič, V., Roubíček, T., (2014). A simple and efficient BEM implementation of quasistatic linear visco-elasticity. International Journal of Solids and Structures 51:2261-2271. doi:10.1016/j.ijsolstr.2014.02.028
  » https://doi.org/10.1016/j.ijsolstr.2014.02.028
• Pascon, J.P., Coda, H.B., (2017). Finite deformation analysis of visco-hyperelastic materials via solid tetrahedral finite elements. Finite Elements in Analysis and Design 133:25-41. doi:10.1016/j.finel.2017.05.007
  » https://doi.org/10.1016/j.finel.2017.05.007
• Rabelo, J.M.G., Becho, J.S., Greco, M., Cimini, C.A.J., (2018). Modeling the creep behavior of GRFP truss structures with Positional Finite Element Method. Latin American Journal of Solids and Structures 15:e17. doi:10.1590/1679-78254432
  » https://doi.org/10.1590/1679-78254432
• Ramos, É.S., Carrazedo, R., (2020). Cross-section modeling of the non-uniform corrosion due to chloride ingress using the positional finite element method. Journal of the Brazilian Society of Mechanical Sciences and Engineering 42:1-18. doi:10.1007/s40430-020-02627-5
  » https://doi.org/10.1007/s40430-020-02627-5
• Riobom Neto, C.G., Santiago, J.A.F., Telles, J.C.F., Costa, E.G.A., (2021). An accurate Galerkin-BEM approach for the modeling of quasi-static iscoelastic problems. Engineering Analysis with Boundary Elements, 130: 94-108. doi: 10.1016/j.enganabound.2021.05.013
  » https://doi.org/10.1016/j.enganabound.2021.05.013
• Sá, M.F., Gomes, A.M., Correia, J.R., Silvestre, N., (2011a). Creep behavior of pultruded GFRP elements–Part 1: Literature review and experimental study. Composite Structures 93(9):2450-2459. doi:10.1016/j.compstruct.2011.04.013
  » https://doi.org/10.1016/j.compstruct.2011.04.013
• Sá, M.F., Gomes, A.M., Correia, J.R., Silvestre, N., (2011b). Creep behavior of pultruded GFRP elements–Part 2: Analytical study. Composite structures 93(9), 2409-2418. doi:10.1016/j.compstruct.2011.04.001
  » https://doi.org/10.1016/j.compstruct.2011.04.001
• Sampaio, M.S.M., Paccola, R.R., Coda, H.B., (2015). A geometrically nonlinear FEM formulation for the analysis of fiber reinforced laminated plates and shells. Composite Structures 119:799-814. doi:10.1016/j.compstruct.2014.09.009
  » https://doi.org/10.1016/j.compstruct.2014.09.009
• Shao, Y., Shanmugam, J., (2004). Deflection creep of pultruded composite sheet piling. Journal of composites for construction 8:471-479. doi:10.1061/(ASCE)1090-0268(2004)8:5(471)
  » https://doi.org/10.1061/(ASCE)1090-0268(2004)8:5(471)
• Shenoi, R.A., Allen, H.G., Clark, S.D., (1997). Cyclic creep and creep-fatigue interaction in sandwich beams. The Journal of Strain Analysis for Engineering Design 32(1):1-18. doi:10.1243/0309324971513175
  » https://doi.org/10.1243/0309324971513175
• Xiao, W., Li, L., Wang, M., (2017). Propagation of in-plane wave in viscoelastic monolayer graphene via nonlocal strain gradient theory. Applied Physics A 123(6):388. doi:10.1007/s00339-017-1007-1
  » https://doi.org/10.1007/s00339-017-1007-1

Edited by

Editor: Rogério José Marczak

Publication Dates

• Publication in this collection
  07 Jan 2022
• Date of issue
  2022

History

• Received
  06 Aug 2021
• Reviewed
  10 Oct 2021
• Accepted
  16 Nov 2021
• Published
  22 Nov 2021