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

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.


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

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 (Π) can be determined by the Principle of Virtual Works as follows: in which  represents the total strain energy and  represents the potential energy of external applied forces, given by (2-3) The stress tensor is represented by , the strain tensor is represented by , the degrees of freedom are represented by   , the equivalent loads related to each degree of freedom are represented by   , and  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: (2-4) Considering the variations regarding the degrees of freedom (nodal parameters in positional formulation), one has: Considering that   represents an arbitrary variation in the degrees of freedom, a nonlinear system of ngl equations is obtained as follows: 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): Expanding Equation (2-7) through a first order Taylor series, one has: () ≅   () +   ,  ()∆  ≅ 0 (where q and r = 1, 2, ..., ndf) in which   () represents the components of the vector of residues,   ,  () represents the components of the Hessian matrix and ∆  represents the components of the nodal position correction vector: 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)  is taken as an undeformed configuration; (2) The components of the vector of residues   () are calculated using Eq.(2-9); Latin American Journal of Solids and Structures, 2022, 19(1), e418 5/21 (3) The components of the Hessian matrix   ,  () are calculated using Eq.(2-10); (4) The components of the nodal position correction vector ∆  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 ∆ 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 ( =  + ∆); (5.2) If the nodal position correction vector ∆ 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 ().
The convergence criterion is adopted based on a relation between Euclidian norms as follows: 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  and its derivatives ,  and ,  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.

STRESS-STRAIN RELATION
The differential form, presented in Equation (3-1), is generally used to describe viscoelastic behavior and stressstrain 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.
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 1 e E 2 .kinematics Juliano dos Santos Becho et al.Solids and Structures, 2022, 19(1), e418 6/21 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 given by the sum of the strain of the elastic element and the strain of the Kelvin-Voigt model, expressed by:

Latin American Journal of
(3-2) (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: in which   represents the Kronecker delta operator,  represents the first Lamé parameter,  represents the second Lamé parameter, and   and   represent the coefficients of viscosity.Since the first and second invariants of strain are given by: in which  represents the Poisson coefficient and E represents the longitudinal elastic modulus.
In a simplified way, for the isotropic case, the coefficients of viscosity   and   can be approximated as equals ( =   =   ), 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 .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: in which  �  represents the undimensionalized constitutive elastic tensor, given by: 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.

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  1 (varying from -1 to 1) and  2 (varying from -1 to 1), as illustrated in Figure 2. In the undeformed configuration a generic point (( 1 ,  2 ), ( 1 ,  2 )) 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 (( 1 ,  2 ), ( 1 ,  2 )), can be mapped after the configuration change of the element.Since the slopes of the crosssections can be determined by the directions of the unit vectors  �⃗ and  � �⃗ , respectively, in the undeformed and deformed configurations.
Thus, the  and  coordinates of a generic point , in the undeformed configuration, can be expressed, respectively, as follows: Latin American Journal of Solids and Structures, 2022, 19(1), e418 8/21 in which  � and  � 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 crosssection.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   ( 1 ).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 (  ) can be expressed, respectively, for the transformations from the dimensionless auxiliary configuration to the undeformed configuration (  0  ) and from the dimensionless auxiliary configuration to the deformed configuration (  1  ), as follows: 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 11 ) and for the transformation from auxiliary configuration to deformed configuration (  1 11 ), can be expressed by: The angles between two perpendicular directions at undeformed configuration and two directions in defomed configuration at the same considered point can be expressed by: Finally, the components of the strain tensor referring to the respective normal and shear strains can be expressed by: 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: The chain rule can be used to exchange integration variable and the total strain energy can be rewritten as: 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.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 ∆) as follows: in which indexes "s" e "s-1" represent current and previous instant of time, respectively, while ∆ represents an imposed time step.Therefore, the stress-strain relation for Boltzmann model can be expressed by: and total strain energy is represented by: 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.

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

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.
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 1 ) equal to 100 GPa, longitudinal elastic modulus (E 2 ) 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 1 is shown in Figure 5(a).It is possible to observe a non-linear relation between the displacements and the longitudinal elastic modulus E 1 .The influence of the longitudinal elastic modulus E 2 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 2 .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 1 or in the longitudinal elastic modulus E 2 .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 1 and the spring with longitudinal elastic modulus E 2 .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.

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.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 1 ) equal to 100 GPa, longitudinal elastic modulus (E 2 ) 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 crosssection 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.

𝑤𝑤(𝑡𝑡)
Latin American Journal of Solids and Structures, 2022, 19(1), e418 13/21 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).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.  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.

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 2 .
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 1 ) equal to 22,5757 kN/mm 2 , longitudinal elastic modulus (E 2 ) equal to 11,0 kN/mm 2 , modulus of viscosity (η) equal to 500,0 kN/mm 2 •day, and Poisson coefficient equal to 0. 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).

Cantilever beam analysis
This example deals with a cantilever beam loaded with a shear load applied on free end, as shown in Figure 12.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 2 , 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.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%.

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.The beam has length (L) equal to 10 m, rectangular cross-section with width equal to 0.01 m and height equal to h 0 at extremities and equal to h 1 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 1 ) equal to 98,00 MN/m 2 , longitudinal elastic modulus (E 2 ) equal to 24.50 MN/m 2 , modulus of viscosity (η) equal to 274,4 MN/m 2 •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 0 (at the ends) and height h 1 (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.
Latin American Journal of Solids and Structures, 2022, 19(1), e418 17/21 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: 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.

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

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

Figure 3
Figure 3 Schematic interpretation of uncoupled models.

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

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

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

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

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

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

Figure 14
Figure 14Supported beam with variable cross-section.
w(t) represents vertical displacement at mid-span and k represents the form factor (determined as 5/6 for rectangular cross-sections).

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

Table 1
Numerical results with both kinematics.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.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). *