# Influence of Retardation Time on Viscoelastic Behavior of Plane Beams Modeled by the Nonlinear Positional Formulation of the Finite Element Method

Authors: Juliano dos Santos Becho; Marcelo Greco
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

A numerical procedure is presented to avoid the divergence problem during the iterative process in viscoelastic analyses. This problem is observed when the positional formulation of the finite element method is adopted in association with the finite difference method. To do this, the nonlinear positional formulation is presented considering plane frame elements with Bernoulli–Euler kinematics and viscoelastic behavior. The considered geometrical nonlinearity refers to the structural equilibrium analysis in the deformed position using the Newton–Raphson iterative method. However, the considered physical nonlinearity refers to the description of the viscoelastic behavior through the adoption of the stress-strain relation based on the Kelvin–Voigt rheological model. After the presentation of the formulation, a detailed analysis of the divergence problem in the iterative process is performed. Then, an original numerical procedure is presented to avoid the divergence problem based on the retardation time of the adopted rheological model and the penalization of the nodal position correction vector. Based on the developments and the obtained results, it is possible to conclude that the presented formulation is consistent and that the proposed procedure allows for obtaining the equilibrium positions for any time step value adopted without presenting divergence problems during the iterative process and without changing the analysis of the final results.


---

## Body

## 1. Introduction

Different numerical formulations have been developed to analyze the viscoelastic behavior of solid materials, as demonstrated in the works of Chen [1], Aköz and Kadioǧlu [2], Beijer and Spoormaker [3], Zheng et al. [4], Mesquita and Coda [5], Galucio et al. [6], Bottoni et al. [7], Payette and Reddy [8], Panagiotopoulos et al. [9], Latorre and Montáns [10], dos Santos Becho et al. [11], Carniel et al. [12], Oliveira and Leonel [13], Pascon and Coda [14], Rabelo et al. [15], and Fernandes et al. [16]. This interest is due to the need for modeling behavior of unconventional materials and the more complex and realistic analyzes of structures and structural components. Thus, structural engineering aims to address the demands for technological advances in the most diverse areas, such as infrastructure, civil construction, mechanical industry, aerospace industry, and others. Part of these studies are related to the analyses of structures with viscoelastic behavior adopting stress-strain relation deduced from rheological models, as demonstrated in the works of Aköz and Kadioǧlu [2], Mesquita and Coda [5], Panagiotopoulos et al. [9], dos Santos Becho et al. [11], Carniel et al. [12], Oliveira and Leonel [13], Pascon and Coda [14], Rabelo et al. [15], and Fernandes et al. [16]. These rheological models are schematic representations that combine elastic and viscous elements to obtain a physical interpretation of stress-strain relationships. Therefore, from the interpretation of these rheological models, it is possible to deduce relations between the components of the stress and strain tensors that are appropriate to the mechanical behavior of interest and necessary for the assessment of the strain energy.It is worth mentioning that part of these studies uses the positional formulation of the finite element method to describe the viscoelastic behavior, as demonstrated in the works of dos Santos Becho et al. [11], Pascon and Coda [14], Rabelo et al. [15], and Fernandes et al. [16]. The use of such a formulation is justified by the capacity and simplicity of its application in the analysis of nonlinear problems as highlighted in Coda and Greco [17] and Greco and Coda [18]. For this reason, although relatively recent, the positional formulation has been the object of study of different researchers, as demonstrated in the works of Coda and Paccola [19], Maciel and Coda [20], Greco and Ferreira [21], de Oliveira and Greco [22], Sampaio et al. [23], de Barros Cavalcante et al. [24], Pascon and Coda [14], Carrazedo and Coda [25], Rabelo et al. [15], and Ramos and Carrazedo [26].However, as shown in the present study, when the positional formulation is used in association with the finite difference method to describe viscoelastic behavior based on rheological models, a divergence problem is observed in the iterative process. Such divergence problem refers to the departure of nodal positions in relation to equilibrium positions throughout the iterative process when time steps smaller than the characteristic delay time of the adopted rheological model are used. This divergence problem is not restricted to positional formulation and is also found in conventional finite element formulations when the nonlinear problem is solved by the implicit method of Newton–Raphson.The first option for complex structural analysis, especially for coupled multiphysics phenomena, is explicit integration schemes. In this case, the smaller the time step, the better the convergence. But for some analyses, it is suitable to use implicit integration schemes. Khan et al. [27] recently published an interesting paper about the influence of creep rate on the viscoelastic response of a polymeric beam. The paper relates the heat generation due to energy dissipation in viscoelastic materials under cyclic loading with its consequences on the overall time-dependent deformations. The authors used the Abaqus software to implement a viscoelastic constitutive model. Implicit time integration is considered at each time increment during nonlinear analysis. The obtained results showed that higher frequencies and higher load amplitudes provide the occurrence of instabilities earlier. It can be noted the model tendency to become unstable for higher frequencies due to the use of the implicit integration scheme. The instability problem cannot be avoided using standard finite element formulation, but it can be avoided using an explicit integration scheme when this type of strategy can be applied.The interest of part of the recent researches in nonlinear behaviors and analysis methods is due to the developments and applications of microscale mechanical systems and nanostructures. 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. [28] and Xiao et al. [29].The viscoelastic behavior is of special importance in polymeric materials and composite materials with the polymeric matrix. Such materials have gained prominence in civil construction, the mechanical industry, the 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) [30–33] and sandwich panels with polymeric core [6, 34–36].Thus, this paper presents a positional formulation of the finite element method that is capable of describing the viscoelastic mechanical behavior in structures discretized by plane frame elements with Bernoulli–Euler kinematics. To do this, the viscoelastic behavior is evaluated through the stress-strain relationship obtained from the Kelvin–Voigt rheological model, which takes into account the time variable. Then, using the presented positional formulation, an analysis is carried out to expose the divergence problem during the iterative process according to the adopted time step. Finally, a numerical procedure based on the retardation time and the penalization of the nodal position correction vector is presented to avoid the divergence problem.


## 2. Positional Formulation of the Finite Element Method

Originally presented in Coda and Greco [17], the positional formulation of the finite element method is based on variational concepts of the minimum 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 fixed coordinate system in space (total Lagrangian description).For quasi-static and conservative structural systems, as presented in Coda and Greco [17], the functional of the total potential energy (Π) can be determined by the principle of virtual works as follows:(1)Π=U+W,where U represents total strain energy and W represents the potential energy of external applied forces, given by(2)U=∫VudV=∫V∫εσdεdV,W=−∑q=1ndfPqXq,where the stress tensor is represented byσ, the strain tensor is represented by ε, the degrees of freedom are represented by Xq, the equivalent loads related to each degree of freedom are represented by Pq, 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 [36], the first variation of total potential energy functional must be zero as follows:(3)δ1Π=0.Considering the variations regarding the degrees of freedom (nodal parameters in positional formulation), one has(4)δ1Π=∂Π∂XqδXq=0.Considering thatδXq represents an arbitrary variation in the degrees of freedom, a linear system with the number of equations equal to the number of degrees of freedom (ndf) is obtained as follows:(5)∂Π∂Xq=∂U∂Xq+∂W∂Xq=0,whereq=1,2,…,ndf.Equation (5) represents the principle of minimum total potential energy. Conceptually, before discretization, the structure has infinite degrees of freedom, and consequently, the system 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.The principle of minimum total potential energy determines that among all possible configurations for a system constituted by a deformable body loaded by external forces, the stationary value of the total potential energy (Π) is the equilibrium configuration [37]. Thus, the equilibrium of the structure will occur when the partial derivative of the total potential energy in relation to the degrees of freedom (nodal parameters) is zero. 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).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 (5) 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 study, the system of nonlinear equations is solved using the Newton–Raphson iterative method, briefly described below, considering the structure discretized by the nodal positions (generalized nodal parameters).To describe the Newton–Raphson iterative method, equation (5) can be rewritten in a compact form as follows:(6)∂Π∂Xq=−gqX=0.Expanding equation (6) through a first-order Taylor series, one has(7)gqX≅gqX+gq,rXΔXr≅0whereq∧r=1,2,…,ndf,where gqX represents the components of the vector of residues, gq,rX represents the components of the Hessian matrix, and ΔXr represents the components of the nodal position correction vector:(8)gqX=Pq−U,q,(9)gq′rX=−U,qr,leading to(10)ΔXr=Pq−U,qU,qr−1.The solution of the nonlinear system of the equation given by equation (7) can be obtained by Newton–Raphson iterative method, graphically represented in Figure 1 and summarized as follows:(1)
X is taken as an undeformed configuration(2)
The components of the vector of residuesgqX are calculated using equation (8)(3)
The components of the Hessian matrixgq,rX are calculated using equation (9)(4)
The components of the nodal position correction vectorΔXr are calculated using equation (10)(5)
The adopted criterion of convergence is verified, considering the following two situations:(i)
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+ΔX)(ii)
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)Figure 1
Representational scheme iterative process of the Newton–Raphson method (force control). The superscripts inX and ΔX refer to the current iterative step.The convergence criterion adopted is based on a relation between Euclidian norms as follows:(11)ΔXX≤tol,where tol represents an adopted numerical tolerance.From the basic formulation developed for the principle of minimum total potential energy and 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 energyU and its derivatives U,q and U,qr in relation to the generalized nodal positions (degrees of freedom), presented in equations (8) and (9). 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. Kelvin–Voigt Rheological Model

In order to determine the strain energy, it is necessary to obtain the stress-strain relation for the model under consideration. In addition, to understand the divergence problem in the iterative process that is characteristic of positional formulation, it is necessary to characterize the mechanical response of the model.The Kelvin–Voigt model is the rheological model considered in the present study, shown in Figure2, as well as adopted in Aköz and Kadioǧlu [2], Mesquita and Coda [38], and Oliveira and Leonel [13]. This model is not capable of describing an instantaneous elastic behavior, describing only a damped elastic behavior (viscoelastic behavior) with a decreasing strain rate over time, typical of solid materials. Thus, for a sufficiently long period, the damped elastic deformation converges to the instantaneous elastic deformation response predicted by a linear elastic model with the same modulus of elasticity.Figure 2
Kelvin–Voigt rheological model.


### 3.1. Stress-Strain Relation

The rheological relation of the Kelvin–Voigt model can be obtained considering that the total strain is equal to the strain in each element; the total stress in the model can be given by the sum of stress in the elastic element and stress in the viscous element, expressed, respectively, by the following equations:(12)ε=εe=εv,(13)σ=σe+σv,where indexes “(e)” and “(v)” refer to elastic and viscous elements, respectively, and the stress-strain relation for each element is given by the following equations:(14)σe=Eεe,(15)σv=ηε˙v,where E represents the modulus of elasticity and η represents the modulus of viscosity.From equations (12) and (13) and considering stress-strain relations, one has(16)σ=Eε+ηε˙,which represents the rheological relation of the Kelvin–Voigt model. Similarly, it is possible to develop stress-strain relation for other rheological models.From the stress-strain relation, the mechanical response of a rheological model can be characterized by creep and relaxation behaviors. Thus, in the next two subsections, these behaviors are detailed for the Kelvin–Voigt model.


### 3.2. Creep Behavior

The creep behavior is related to the gradual evolution of strains over time when the solid is submitted to a constant stress state. Thus, from equation (16) and considering the constant stress state σ=σ0, one has(17)dεdt+Eηε=σ0η,which represents the nonhomogeneous differential equation that rules the creep behavior of the Kelvin–Voigt model.The solution of equation (17) can be obtained from the linear combination of a homogeneous solution (εh) with a particular solution (εp) as follows:(18)εh=Ae−E/ηt,εp=σ0E,εt=Ae−E/ηt+σ0E,where A represents a constant to be evaluated and t represents time variable.Considering an undeformed configuration, the model initial condition is given by the following equation:(19)ε0=A+σ0E=0.Thus, equation (17) can be rewritten as follows:(20)εt=σ0E1−e−E/ηt.Performing the first derivative of equation (20) in relation to time, the strain rate (ε˙) of the Kelvin–Voigt model is given by the following equation:(21)ε˙t=σ0ηe−E/ηt.From equations (20) and (21), it is possible to describe the strain evolution over time of the Kelvin–Voigt model, under a constant stress state as shown in Figure 3. From equation (20), it is possible to evaluate that strain varies over time, presenting zero initial strain and converging to σ0/E over a sufficiently long time. However, from equation (21), it is possible to observe that the strain rate varies over time in a decreasing way. Wherein, the initial deformation rate is equal to σ0/η and reaches zero over a sufficiently long time.Figure 3
Representation of creep by the Kelvin–Voigt rheological model.As shown in Figure3, for a constant strain rate equal to the initial strain rate, the strain would be equal to the maximum strain of the model (ε∞) at the retardation time (tε). Thus, the retardation time can be evaluated as follows:(22)tε=ηE.According to Findley et al. [39], the retardation time is a physical property of materials and is related to the initial strain rate obtained in a creep test. The retardation time is then defined as the time required for the strain, under a state of constant stress, to reach its maximum value if the strain rate remains constant and equal to the initial value. In other words, the retardation time is the necessary time to stop the creep process if the strain rate remains constant and equal to the initial value. According to Marques and Creus [40], the retardation time represents a characteristic of the creep rate. The lower its value, the faster the creep process occurs, and the material is classified as less viscous. The retardation time provides the time estimate required for the creep process to approach its end.Most of the creep process occurs within the retardation time. However, contrary to what is considered in its definition, for an elapsed time equal to the retardation time, the strain will still not be equal to the maximum strain. This is because the strain rate is decreasing, according to equation (21). In addition, from equations (20) and (22), it is possible to observe that for a time t equal to the retardation time tε, one has(23)εtε=σ0E1−1e≅0.632σ0E≅0.632ε∞,where ε∞ represents the maximum strain obtained over a sufficiently long time. Thus, after the retardation time is reached, the strain level is approximately equal to 63.2% of the maximum creep strain of the Kelvin–Voigt model. The exact same value is found for other rheological models, such as Boltzmann and Zener.The retardation time is important in the present study because it is used to control how much the model is able to deform physically in a certain time interval. Thus, avoiding problems of divergence during the iterative method adopted in the positional formulation.


### 3.3. Relaxation Behavior

The relaxation behavior for a specific model is related to the gradual reduction of stress state over time when the solid is submitted to a constant strain state. Thus, from equation (16) and considering the constant strain ε=ε0, one has(24)σ=Eε0.It is possible to note that as strain remains constant, the stress also remains constant over time. Thus, the Kelvin–Voigt model is not capable of reproducing the relaxation behavior.


## 3.1. Stress-Strain Relation

The rheological relation of the Kelvin–Voigt model can be obtained considering that the total strain is equal to the strain in each element; the total stress in the model can be given by the sum of stress in the elastic element and stress in the viscous element, expressed, respectively, by the following equations:(12)ε=εe=εv,(13)σ=σe+σv,where indexes “(e)” and “(v)” refer to elastic and viscous elements, respectively, and the stress-strain relation for each element is given by the following equations:(14)σe=Eεe,(15)σv=ηε˙v,where E represents the modulus of elasticity and η represents the modulus of viscosity.From equations (12) and (13) and considering stress-strain relations, one has(16)σ=Eε+ηε˙,which represents the rheological relation of the Kelvin–Voigt model. Similarly, it is possible to develop stress-strain relation for other rheological models.From the stress-strain relation, the mechanical response of a rheological model can be characterized by creep and relaxation behaviors. Thus, in the next two subsections, these behaviors are detailed for the Kelvin–Voigt model.


## 3.2. Creep Behavior

The creep behavior is related to the gradual evolution of strains over time when the solid is submitted to a constant stress state. Thus, from equation (16) and considering the constant stress state σ=σ0, one has(17)dεdt+Eηε=σ0η,which represents the nonhomogeneous differential equation that rules the creep behavior of the Kelvin–Voigt model.The solution of equation (17) can be obtained from the linear combination of a homogeneous solution (εh) with a particular solution (εp) as follows:(18)εh=Ae−E/ηt,εp=σ0E,εt=Ae−E/ηt+σ0E,where A represents a constant to be evaluated and t represents time variable.Considering an undeformed configuration, the model initial condition is given by the following equation:(19)ε0=A+σ0E=0.Thus, equation (17) can be rewritten as follows:(20)εt=σ0E1−e−E/ηt.Performing the first derivative of equation (20) in relation to time, the strain rate (ε˙) of the Kelvin–Voigt model is given by the following equation:(21)ε˙t=σ0ηe−E/ηt.From equations (20) and (21), it is possible to describe the strain evolution over time of the Kelvin–Voigt model, under a constant stress state as shown in Figure 3. From equation (20), it is possible to evaluate that strain varies over time, presenting zero initial strain and converging to σ0/E over a sufficiently long time. However, from equation (21), it is possible to observe that the strain rate varies over time in a decreasing way. Wherein, the initial deformation rate is equal to σ0/η and reaches zero over a sufficiently long time.Figure 3
Representation of creep by the Kelvin–Voigt rheological model.As shown in Figure3, for a constant strain rate equal to the initial strain rate, the strain would be equal to the maximum strain of the model (ε∞) at the retardation time (tε). Thus, the retardation time can be evaluated as follows:(22)tε=ηE.According to Findley et al. [39], the retardation time is a physical property of materials and is related to the initial strain rate obtained in a creep test. The retardation time is then defined as the time required for the strain, under a state of constant stress, to reach its maximum value if the strain rate remains constant and equal to the initial value. In other words, the retardation time is the necessary time to stop the creep process if the strain rate remains constant and equal to the initial value. According to Marques and Creus [40], the retardation time represents a characteristic of the creep rate. The lower its value, the faster the creep process occurs, and the material is classified as less viscous. The retardation time provides the time estimate required for the creep process to approach its end.Most of the creep process occurs within the retardation time. However, contrary to what is considered in its definition, for an elapsed time equal to the retardation time, the strain will still not be equal to the maximum strain. This is because the strain rate is decreasing, according to equation (21). In addition, from equations (20) and (22), it is possible to observe that for a time t equal to the retardation time tε, one has(23)εtε=σ0E1−1e≅0.632σ0E≅0.632ε∞,where ε∞ represents the maximum strain obtained over a sufficiently long time. Thus, after the retardation time is reached, the strain level is approximately equal to 63.2% of the maximum creep strain of the Kelvin–Voigt model. The exact same value is found for other rheological models, such as Boltzmann and Zener.The retardation time is important in the present study because it is used to control how much the model is able to deform physically in a certain time interval. Thus, avoiding problems of divergence during the iterative method adopted in the positional formulation.


## 3.3. Relaxation Behavior

The relaxation behavior for a specific model is related to the gradual reduction of stress state over time when the solid is submitted to a constant strain state. Thus, from equation (16) and considering the constant strain ε=ε0, one has(24)σ=Eε0.It is possible to note that as strain remains constant, the stress also remains constant over time. Thus, the Kelvin–Voigt model is not capable of reproducing the relaxation behavior.


## 4. Particularization of the Positional Formulation

Considering the stress-strain relation expressed in equation (16), it is possible to determine the strain energy. To do this, it is necessary to particularize the kinematics and the strain measure for the element finite model adopted in the present study. Thus, this section presents the particularization of the kinematics and the strain measure for the plane frame elements with Bernoulli–Euler kinematics.In order to determine the strain measure, 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 0 to 1) and ξ2 (varying from −1 to 1), as illustrated in Figure 4.Figure 4
Geometry parameterization of a plane frame element with Bernoulli–Euler kinematics.In the undeformed configuration, a generic pointpxξ1,ξ2,yξ1,ξ2 can be mapped by centroidal position and slope of the cross-section, as shown in Figure 4. Similarly, the same point in the deformed configuration, represented by PXξ1,ξ2,Yξ1,ξ2, can be mapped after the configuration change of the element. Thus, for undeformed and deformed configurations, respectively, we have the following expressions for mapping a generic point as a function of dimensionless variables and the slope of the cross-section:(25)pxξ1,ξ2,yξ1,ξ2=p˜xξ1,yξ1+h2ξ2n⟶ξ1,(26)PXξ1,ξ2,Yξ1,ξ2=P˜Xξ1,Yξ1+h2ξ2N⟶ξ1,where p˜ and P˜ represent the points that locate the centroid of cross-sections of the finite element, h represents the height of the cross-section of the element, and n⟶ and N⟶ represent the unit vectors that define the inclinations of the cross-sections, respectively, in undeformed and deformed configurations, as follows:(27)n⟶ξ1=sinθξ1,cosθξ1,(28)N⟶ξ1=sinΘξ1,cosΘξ1.Thus, from equations (25)–(28) and Figure 4, the x and y coordinates of a generic point p, in the undeformed configuration, can be expressed, respectively, as follows:(29)xξ1,ξ2=x˜ξ1−h2ξ2sinθξ1,(30)yξ1,ξ2=y˜ξ1+h2ξ2cosθξ1,where x˜ and y˜ represent coordinates of the common point between the cross-section and the neutral axis and θ represents the angle between the cross-section and the horizontal axis.Similarly, theX and Y coordinates of a generic point P, in deformed configuration, can be expressed, respectively, as follows:(31)Xξ1,ξ2=X˜ξ1−h2ξ2sinΘξ1,(32)Yξ1,ξ2=Y˜ξ1+h2ξ2cosΘξ1,where X˜ and Y˜ represent coordinates of the common point between the cross-section and the neutral axis and Θ represents the angle between the cross-section and the horizontal axis.In order to enable the analysis by the positional formulation of the finite element method, it is necessary to discretize the domain, that is, to leave it in the function of discrete parameters. In the present case, for discrete representation, the considered parameters are the nodal positions, and the formulation is said to be positional because of this. It is important to observe that the slope of the cross-section is a nodal parameter that is dependent on the nodal positions (nodal coordinates).In this formulation, the finite elements of two nodes are considered. Thus, the mapping of nodal coordinates, both in the undeformed and in the deformed configurations, can be rewritten in terms of the positions of these two nodes and of the functions that relate these nodes to dimensionless variables. To do this, the geometry parameterization is considered based on the dimensionless auxiliary configuration, as shown in Figure4.As shown in Figure4, ω represents the domain of an element with two nodes, in the undeformed configuration, and Ω represents the domain of the same element in the deformed configuration. Thus, x1, y1, and θ1 represent the nodal parameters of node 1, and x2, y2, and θ2 represent the nodal parameters of node 2 in the undeformed configuration. On the other hand, X1, Y1, and Θ1 represent the nodal parameters of node 1, and X2, Y2, and Θ2 represent the nodal parameters of node 2 in the deformed configuration.It is possible to describe the geometry that represents the centroidal line of the plane frame element with two nodes using a linear relation betweenξ1 and the X axis and a cubic relation between ξ1 and the Y axis. Thus, the coordinates of a point P˜ can be obtained, respectively, as follows:(33)X˜ξ1=X1+X2−X1ξ1,(34)Y˜ξ1=cξ13+dξ12+eξ1+f,where parameters c, d, e, and f can be evaluated from the element boundary conditions in auxiliary configuration, as follows:(35)Y˜ξ1=0=f=Y1,(36)dY˜dξ1ξ1=0=e=dY˜dX˜dX˜dξ1ξ1=0=tanΘ1X2−X1,(37)dY˜dξ1ξ1=1=3c+2d+e=dY˜dX˜dX˜dξ1ξ1=1=tanΘ2X2−X1,(38)Y˜ξ1=1=c+d+tanΘ1X2−X1+Y1=Y2.From equations (37) and (38), one has(39)c=tanΘ2+tanΘ1X2−X1−2Y2−Y1,(40)d=3Y2−Y1−tanΘ2+2tanΘ1X2−X1.It is important to note that considering the Bernoulli–Euler kinematics, in which the plane cross-section remains plane and orthogonal to the element’s neutral axis after deformation occurs, the cross-section rotation, represented byΘξ1, can be obtained as follows:(41)Θξ1=arctandY˜dX˜=arctandY˜dξ1dξ1dX˜=arctan3cξ12+2dξ1+eX2−X1.Nodal parameters presented in equations (33)–(41) are related to deformed configuration. The same equations are valid to undeformed configuration, just considering nodal parameters x1, x2, y1, y2, θ1, and θ2 instead of X1, X2, Y1, Y2, Θ1, and Θ2, respectively. Furthermore, the same consideration is valid to obtain coordinates x˜ξ1 and y˜ξ1 and rotation θξ1. Thus, mapping coordinates in undeformed and deformed configurations can be fully described in terms of dimensionless variables and nodal parameters considering equations (29)–(32).Through geometry mapping and considering Bernoulli–Euler kinematics, it is possible to write a strain measure for each finite element. To do this, an infinitesimal fiber with lengthdξ1 and parallel to the neutral plane in the auxiliary nondimensional configuration is initially considered, as shown in Figure 5. When the element passes from the auxiliary nondimensional configuration to undeformed configuration and deformed configuration, the fiber element’s length changes to ds and dS, respectively.Figure 5
Strain measure parametrization.Considering a straight finite element at the undeformed configuration, the stretch from auxiliary nondimensional configuration to undeformed configuration (λ˜) and the stretch from auxiliary nondimensional configuration to deformed configuration (Λ˜) at neutral axis are defined, respectively, by the following equations:(42)λ˜=dsdξ1=dx˜dξ12+dy˜dξ12=x2−x12+y2−y12,Λ˜=dSdξ1=dX˜dξ12+dY˜dξ12=X2−X12+3cξ12+2dξ1+e2.Thus, the stretch rate can be calculated from the relation between the stretch of deformed configuration and the stretch of undeformed configuration. Therefore, an objective strain measure, as Biot (ε˜), can be used for small strains in neutral axis as follows:(43)ε˜=Λ˜λ˜−1=X2−X12+3cξ12+2dξ1+e2x2−x12+y2−y12−1.It is worth mentioning that this formulation can be successfully applied to nonlinear problems with large displacements and rotations but with normal strain limited up to 5%.For another fiber not on neutral axis, longitudinal strain can be calculated according to the Bernoulli–Euler kinematics, as follows:(44)ε=ε˜−h2ξ21r,where 1/r represents the exact neutral axis curvature in deformation configuration, represented by the following equation:(45)1r=dX/dξ1d2Y/dξ12−d2X/dξ12dY/dξ1dX/dξ12+dY/dξ123=X2−X126cξ1+2dX2−X12+3cξ12+2dξ1+e23.From equations (43)–(45), it is possible to calculate the normal strain in a fiber, considering the nondimensional coordinates (ξ1 and ξ2) and the nodal parameters (X1, Y1, Θ1, X2, Y2, and Θ2).Considering Kelvin–Voigt stress-strain relation presented in equation (16) and strain measure presented in equation (44), the total strain energy of viscoelastic analysis can be expressed as follows:(46)U=∫V∫εEε+ηε˙dεdV.The chain rule can be used to exchange integration variables as follows:(47)dε=dεdXqdXq=ε,qdXq.Thus, the total strain energy can be rewritten as follows:(48)U=∫V∫XEε+ηε˙ε,qdXdV.Performing the first derivative of total strain energy regarding nodal parameters, one has(49)U,q=∫VEε+ηε˙ε,qdV.Performing the second derivative of total strain energy regarding nodal parameters, one has(50)U,qr=∫VEε,r+ηε˙,rε,q+Eε+ηε˙ε,qrdV.Finally, considering equations (49) and (50), it is possible to apply the Newton–Raphson iterative method as described in Section 2. Thus, the formulation can be used to analyze structures that are discretized by plane frame elements and with a characteristic viscoelastic behavior of the Kelvin–Voigt model.


## 5. Numerical Procedure to Assess Strain Rate

Due to the characteristic viscoelastic behavior of the model, the stress-strain relation is dependent on time, as well as the strain energy and its derivatives, expressed by equations (48)–(50). This time dependence is represented by the strain rate. Thus, an adequate procedure is necessary to evaluate the strain rate over time.In the present study, this assessment is performed using the regressive finite differences method, that is, the rates are assessed by the difference between the variable value at the current moment and the variable value at the previous instant divided by the time step itself. Thus, the strain rate can be expressed as follows:(51)ε˙=εs−εs−1Δt,where indexes “s” and “s−1” represent current and previous instant of time, respectively, while Δt represents an imposed time step.Therefore, the stress-strain relation of the Kelvin–Voigt model can be expressed as follows:(52)σs=Eεs+ηεs−εs−1Δt,and total strain energy is represented by the following equation:(53)Us=∫V∫XEεs+ηεs−εs−1Δtεs,qdXdV.It is important to note that the approach used here to evaluate strain rate is different from the approach used in the works of dos Santos Becho et al. [11] and Rabelo et al. [15]. In such works, the strain rate is evaluated by direct use of chain rule in terms of nodal positions as follows:(54)ε˙=dεdt=dεdXqdXqdt=ε,qX˙q,where X˙ represents position rate, evaluated using finite differences as follows:(55)X˙q=Xqs−Xqs−1Δt.From equation (55), it is possible to note that in the evaluation of position rate, the effects of rigid body movement are computed. Thus, such an approach can be used to analyze the mechanical behavior of finite truss elements, as long as the relative axial position rate between the element nodes is evaluated, as shown in Rabelo et al. [15]. However, for frame elements, as in the present study, the segregation of rigid body movement is not so simple, and the adoption of equation (51) to evaluate the strain rate is more appropriate.


## 6. Influence of the Adopted Time Step on the Creep Response

An original proposal to avoid the divergence problem during the iterative process in nonlinear viscoelastic analyzes is presented here. This problem is observed in nonlinear analyzes when the adopted time steps are smaller than the retardation time. In this case, there is a divergence between the subsequent nodal positions along the iterative process, making it impossible to obtain the equilibrium positions.To better explain the divergence problem, some results obtained from the numerical analysis of a tensioned viscoelastic fixed bar (Figure6) are presented. The Kelvin–Voigt viscoelastic model is considered, with the modulus of elasticity equal to E=100GPa and modulus of viscosity η=1000GPa·s. A tension load P=1000MN is applied at the free end.Figure 6
Tensioned viscoelastic fixed bar.To evaluate the evolution of the iterative process in viscoelastic behavior, there are adopted four distinct time steps (9 s, 10 s, 11 s, and 12 s). The obtained results are based on the equilibrium position of the free end node.To illustrate the influence of time step in the iterative process, Figures7–10 present the evolution of the residues vector gqX and the evolution of the nodal position correction vector ΔX of nodal positions, considering just the components related to the node at the free end. The first 20 iterations of first time step are presented for each adopted time step. The evolution of the Hessian matrix gq,rX is not presented because its value remains constant and equal to 10,000 MN/m.Figure 7
Evolution of the nodal position correction vector and the residue vector forΔt = 12 s.Figure 8
Evolution of the nodal position correction vector and the residue vector forΔt = 11 s.Figure 9
Evolution of the nodal position correction vector and the residue vector forΔt = 10 s.Figure 10
Evolution of the nodal position correction vector and the residue vector forΔt = 9 s.The retardation time is equal to 10 s in this specific analysis. Thus, it is possible to assess the convergence of the iterative process around this value. Besides, it is possible to observe that for time steps larger than retardation time, the iterative process converges, showing an increasing reduction of residues vector and nodal position correction vector, as presented in Figures7 and 8. On the other hand, for a time step equal to the retardation time, the iterative process does not converge. In this case, the same values of residues vector and nodal position correction vector are obtained over the iterative process, alternating only the signs. That is, the same correction obtained in the first iteration is repeated in the second iteration but with the changed signal. The position of the free end node of the bar oscillates between a corrected position and the initial position indefinitely, as shown in Figure 9. Finally, for time steps smaller than the retardation time, the iterative process diverges, gradually increasing the residues vector and nodal position correction vector, as shown in Figure 10. These same behaviors were observed for time steps higher and lower than those presented in this analysis and for other rheological models such as the Boltzmann model and the Zener model.The relation between time step and evolution of the iterative process is also observed in plots related to position versus time and force versus position. Figure11 presents the positions of the free end node during the first four iterations in the equilibrium position search process for the first time step. Figure 12 presents the same positions plotted in load versus position graphics. Additionally, the equilibrium positions at the end of each time step are shown in Figure 11. It is important to observe that for time steps smaller than retardation time, these equilibrium positions were only possible to obtain when the divergence problem was solved. Besides, there were plotted the analytical solutions for positions over time, obtained from equation (20).Figure 11
Iterative process in position versus time graphics.Figure 12
Iterative process in load versus position graphics.From Figures11 and 12, it is possible to observe that for the adoption of time step larger than the retardation time, positions of the free end node tend to convergence along the iterative process to the equilibrium value. While adopting a time step equal to the retardation time, the positions remain oscillating between the same values around the equilibrium position. Finally, adopting time steps smaller than the retardation time, the positions show a tendency of divergence along the iterative process, moving away from the equilibrium position. In addition, it is possible to observe that the smaller the time step, the closer the numerical equilibrium position is to the analytical response. This is due to the fact that for smaller time steps, the temporal discretization becomes more refined.To clarify how the divergence problem can be solved, it is necessary to return to the system of equations obtained from the Newton–Raphson procedure and expressed by equation (7). To solve this system, the vector of residues and the Hessian matrix need to be evaluated from equations (8) and (9). Therefore, as demonstrated, it is necessary to obtain the strain energy and its first and second derivatives, which, considering the Kelvin–Voigt model, are expressed by equations (48)–(50). However, the external force term Pq can be expressed as follows:(56)Pq=∫Aσ0dA=∫AEε∞dA,where ε∞ represents the final viscoelastic strain of the creep analysis and A represents the cross-sectional area of the bar.Thus, equation (7) can be rewritten as follows:(57)ΔXr=∫AEε∞dA−∫VEε+ηε˙ε,qdV∫VEε,r+ηε˙,rε,q+Eε+ηε˙ε,qrdV.Equation (57) defines the nodal position correction vector to update the deformed configuration at each iteration.Dividing both numerator and denominator by the modulus of elasticity, one has(58)ΔXr=∫Aε∞dA−∫Vε+η/Eε˙ε,qdV∫Vε,r+η/Eε˙,rε,q+ε+η/Eε˙ε,qrdV.Considering the retardation time, equation (58) can be rewritten as follows:(59)ΔXr=∫Aε∞dA−∫Vε+tεε˙ε,qdV∫Vε,r+tεε˙,rε,q+ε+tεε˙ε,qrdV.Considering the strain rate approximate by finite differences, one has(60)ΔXr=∫Aε∞dA−∫Vεs+tεεs−εs−1/Δtε,qdV∫Vε,rs+tεε,rs−ε,rs−1/Δtε,q+εs+tεεs−εs−1/Δtε,qrdV.Furthermore, it is important to note that the time step can be expressed as a fraction of the retardation time as follows:(61)Δt=fεtε,where fε is a retardation factor.Therefore, equation (60) can be rewritten as follows:(62)ΔXr=∫Aε∞dA−1/fε∫Vfεεs+εs−εs−1ε,qdV∫Vε,rs+1/fεε,rs−ε,rs−1ε,q+εs+1/fεεs−εs−1ε,qrdV.Considering the tensioned bar shown in Figure6, the first derivative of strain is given by ε,q=1/L, and the second derivative of strain is given by ε,qr=0. Thus, after evaluating the integrals, equation (62) can be rewritten as follows:(63)ΔXr=ε∞−1fεfεεs+εs−εs−1.In the first iteration,εs=εs−1=0. Therefore, considering equation (63) and Figure 13, it is possible to note that for a time step greater than retardation time, a strain ε∞ is a possible result. For a time step greater than retardation time, equation (62) gives a nodal position correction vector (ΔX) compliant with a feasible strain (ε∞). On the other hand, for a time step smaller than retardation time, a strain equal to ε∞ is a not possible result. However, a strain equal to fεε∞ would be a feasible strain, as shown in Figure 13. Thus, multiplying both sides of equation (63) by the retardation factor (fε), one has a nodal position correction vector penalized fεΔX and compatible with a feasible strain fεε∞.Figure 13
Relation between time step, retardation time, and possible strain.From this analysis, it is possible to observe that the iterative process must be carried out in two ways. The first one refers to cases in which the time step is larger than the retardation time. For this case, the nodal position correction vector does not need to be penalized by the retardation factor, as shown in Figure14. The second one refers to cases in which the time step is smaller than the retardation time. For this case, the nodal position correction vector needs to be penalized by the retardation factor, as shown in Figure 15, to avoid the divergence problem.Figure 14
Iterative process in case ofΔt>tε. The superscripts in X and ΔX refer to the current iterative step.Figure 15
Iterative process in case ofΔt<tε. The superscripts in X and ΔX refer to the current iterative step.Numerically, by adopting the time step equal to 12 s and using equation (63), it is possible to follow the results of the variables along the iterative process, as shown in Table 1. Similarly, results over the iterative process for the other adopted time steps can be obtained. Thus, Tables 2 and 3 present the values of the nodal position correction vector over the first 10 iterations for the time steps of 11 s, 10 s, and 9 s. To obtain the results presented in Table 2, the procedure for penalizing the nodal position correction vector is not considered, while to obtain the results presented in Table 3, this procedure is considered.Table 1
Results of the iterative process with a time step of 12 s and based on equation (63).
StepIterationε∞εsεs−1ΔXm110.5000,5120.50.5000000−0.416667130.50.08333300.347222140.50.4305550−0.289350150.50.14120500.241125160.50.3823300−0.200940⋮⋮⋮⋮⋮⋮11050.50.2727280−1.00E−0811060.50.27272808.40E−09210.50.2727280.2727280.227272220.50.5000000.272728−0.189395230.50.3106050.2727280.157828240.50.4684340.272728−0.131525⋮⋮⋮⋮⋮⋮21020.50.3966940.272728−4.00E−0821030.50.3966940.2727285.10E−09310.50.3966940.3966940.1033063⋮⋮⋮⋮⋮Table 2
Evolution of the nodal position correction vector in the iterative process without using the retardation factor.
IterationΔXmΔt=11sΔt=10sΔt=9s10.500000.500000.500002−0.45454−0.50000−0.5555530.413220.500000.617284−0.37566−0.50000−0.6858750.341500.500000.762076−0.31056−0.50000−0.8467670.282230.500000.940838−0.25657−0.50000−1.0453890.233250.500001.1615210−0.21204−0.50000−1.29059⋮⋮⋮⋮Table 3
Evolution of the nodal position correction vector in the iterative process using the retardation factor.
IterationfεΔXmΔt=11sΔt=10sΔt=9s10.550000.500000.450002−0.60500−0.50000−0.4050030.665500.500000.364504−0.73205−0.50000−0.3280550.805260.500000.295256−0.88578−0.50000−0.2657270.974360.500000.239158−1.07180−0.50000−0.2152491.178970.500000.1937110−1.29687−0.50000−0.17434⋮⋮⋮⋮From Table2, it is possible to observe that the iterative processes present convergence only for the time step equal to 11 s. On the other hand, from Table 3, it is possible to observe that the iterative processes present convergence only for a time step equal to 9 s. Thus, it is possible to note that by adopting the time step greater than the retardation time, the iterative process presents convergence using equation (63) without penalizing the nodal position correction vector. While adopting the time step smaller than the retardation time, the iterative process presents convergence using the penalty procedure of the nodal position correction vector. These same behaviors are obtained with time steps higher or lower than those adopted. Furthermore, it is possible to observe that by adopting a time step equal to the retardation time, the retardation factor is equal to 1. In this case, convergence is not observed either by adopting the penalty or without adopting the penalty procedure of the correction vector of the nodal positions. This behavior is attributed to the fact that the time step equal to the retardation time is at the limit between the convergence behavior and the divergence behavior. The convergence, in this case, is guaranteed by considering the retardation factor fε equal to 1 minus one small residue, such as fε=1−1⋅10−8.Finally, it is important to note that the retardation factor is dimensionless, obtained by a relation between a parameter related to the temporal discretization and the physical properties of the material. In addition, it is important to note that the retardation factor is a limit value. Thus, the retardation factor can be lower than or equal to the ratio between the time step and the retardation timeΔt/tε, and it can assume a maximum value equal to 1 (in which case, Δt>tε). The use of the retardation factor reduces the nodal position correction vector as a whole but does not change the nature of its components or the proportion between its components. Thus, the retardation factor does not interfere in the results of the equilibrium position itself but in the iterative process, ensuring its convergence. This can be verified by the equality of the results obtained when using a retardation factor smaller than the one calculated by the ratio between the time step and the retardation time. In this case, the equilibrium positions at the end of each step do not change. However, a greater number of iterations are required because the lower the retardation factor, the smaller the portion of the nodal position correction vector that is being considered in each iteration. Therefore, the retardation factor does not change the results; it only guarantees the convergence, reducing the march of the iterative process and, consequently, increasing the number of iterations. Similarly, in the case of a time step greater than the retardation time (Δt>tε), in which it is not necessary to use the retardation factor, it is possible to use a retardation factor smaller than 1, obtaining the same results, but with a greater number of iterations.The proposed procedure to avoid the divergence problem during the iterative process can be summarized as shown in Figure16.Figure 16
Simplified procedure to avoid the problem of divergence in the iterative process.Additionally, from Figure12, the geometric interpretation of the iterative process of the Newton–Raphson method for the viscoelastic behavior of a model with instantaneous elastic strain can be illustrated as shown in Figure 17.Figure 17
Geometric interpretation of the iterative process of the Newton–Raphson method. The superscripts inX and ΔX refer to the current iterative step, and the subscripts “(E)” and “(V)” refer to the instantaneous elastic response and the viscoelastic response, respectively.The developments and results presented in this section are based on the analysis of a tensioned bar with viscoelastic behavior described by the Kelvin–Voigt model. This is a didactic option, but likewise, it can also be developed for other rheological models and for analyzes involving bending and shear effects. However, a considerable algebraic effort is required due to the additional terms and effects, making it difficult to expose the problem of divergence.


## 7. Application Problems

In order to demonstrate the validity of the proposed procedure, five application problems are presented using the developed formulation and adopting the computational procedure described in Figure16.


### 7.1. Analysis of a Tensioned Viscoelastic Bar

A classical example is presented here. This example is used by different authors to calibrate viscoelastic models because it is a typical case of a plane state of stress. The example deals with a tensioned bar supported on the left and bottom faces. Other faces are free to move, and a traction force is applied on the free end, on the right, as shown in Figure18.Figure 18
Tensioned bar.The bar has a length (L) of 800 mm, a cross-section height (h) of 100 mm, and unitary width (b). The traction force is equal to 0.005 kN/mm2.The Kelvin–Voigt viscoelastic model is considered, with the modulus of elasticity (E) equal to 11 kN/mm2 and modulus of viscosity (η) equal to 500 kN/mm2·day. Thus, the retardation time (tε) is equal to 45.4545 days. Regarding spatial discretization, 10 finite elements are adopted, with 10 Gauss points along the height and 10 Gauss points along the length. Regarding temporal discretization, 5 different time steps were used (1 day, 5 days, 10 days, 25 days, and 50 days). No inertial effects were considered, and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions).Numerical results for axial displacements at the free end are presented in Figure19. The presented analytical results and numerical results are also available in Mesquita and Coda [38]. The analytical results are obtained from equation (20), and the results available in Mesquita and Coda [38] are achieved with a mesh of 8 × 4 cubic boundary elements (boundary element method).Figure 19
Axial displacements over time.From Figure19, it is possible to note that a reduction of time step makes numerical results converge to the analytical solution. This is due to the fact that the smaller the time steps, the more refined the temporal discretization becomes. Therefore, the mechanical behavior is consistent with the analytical solution and in agreement with the results available in Mesquita and Coda [38]. Besides, it is important to note that the convergence for time steps smaller than retardation time is only possible using the proposed procedure. Therefore, these results validate the proposed procedure to avoid numerical divergence for time steps smaller than the retardation time.Figure20 presents numerical results related to the creep-recovery process, considering a time step equal to 1 day. For the time greater than 200 days, the external applied force (P) is reduced to zero. It is possible to note that the obtained results for recovery analysis are consistent and in agreement with the results available in Mesquita and Coda [38].Figure 20
Creep-recovery process.


### 7.2. Cantilever Beam Loaded by Shear Force

This example deals with an analysis of a cantilever beam loaded by a shear force applied on the free end, as shown in Figure21.Figure 21
Cantilever beam.Geometry, rheological parameters, and spatial discretization are the same as the ones presented in the previous example. The shear force is equal to 0.005 kN/mm2, applied at the initial time, and removed after 453 days. Thereafter, the analysis continues for 267 days without applied forces in order to observe the viscoelastic recovery process. Regarding temporal discretization, 720 time steps equal to 1 day are considered. No inertial effects were considered, and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions).Numerical results for vertical displacements in the middle point of the right face are presented in Figure22, for the creep-recovery process. The numerical results available in Panagiotopoulos et al. [9] made with a boundary element method are also presented here. The results are in good agreement with the numerical results available in the literature. A portion of the difference between the results is attributed to the effects of shear, not computed by the adopted Bernoulli–Euler kinematics.Figure 22
Vertical displacements of the free end over time in the creep-recovery process.


### 7.3. Supported Beam with Uniformly Distributed Loading

This example deals with an analysis of a double supported beam loaded by a uniformly distributed force applied on the upper face, as shown in Figure23.Figure 23
Supported beam.The beam has the length (L) equal to 10 m, cross-section height (h) equal to 0.5 m, and width (b) equal to 2 m. The uniformly distributed force (P) is equal to 10 N/m, applied at the initial time, and removed after 3 s.The Kelvin–Voigt viscoelastic model is considered, with the modulus of elasticity (E) equal to 98 MN/m2 and modulus of viscosity (η) equal to 27.44 MN/m2·s. Thus, the retardation (tε) is equal to 0.28 s. Regarding spatial discretization, 10 finite elements are adopted, with 10 Gauss points along the height and 10 Gauss points along the length. Regarding temporal discretization, three different time steps were used (0.5 s, 0.1 s, and 0.01 s). No inertial effects were considered, and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions).The numerical results for vertical displacements at the midspan are presented in Figure24. The analytical results and numerical results available in Aköz and Kadioǧlu [2] are also presented here. Results are in good agreement with analytical solutions and numerical results available in the literature. It should be noted that the analytical solution is presented in Chen [1] and given by the following equation:(64)wt=5PL432bh31E1−e−Et/η.Figure 24
Vertical displacements at the midspan over time.It is important to note that the convergence for time steps smaller than retardation time is only possible using the proposed procedure to avoid numerical divergence. Besides, it is possible to note that the reduction of time step makes numerical results converge to analytical solutions. This is due to the fact that the smaller the time steps, the more refined the temporal discretization becomes.


### 7.4. Supported Beam with Uniformly Distributed Loading Applied Gradually

The same double supported beam presented in the previous example (Figure23) is analyzed here. However, in this example, distributed loading is applied gradually, as shown in Figure 25.Figure 25
Distributed loading applied gradually over time.Geometry, rheological parameters, and spatial discretization are the same as the ones presented in the previous example. Regarding temporal discretization, 500 time steps equal to 0.01 s were considered in the analysis. No inertial effects were considered, and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions).The numerical results for displacements at the midspan are presented in Figure26. The numerical results available in Aköz and Kadioǧlu [2] are also presented. Obtained results are in good agreement with the numerical results available in the literature.Figure 26
Vertical displacements over time at the midspan with distributed loading applied gradually.


### 7.5. Cylindrical Pressure Vessel

The last numerical example is an engineering application where the pressure is applied gradually, and the Kelvin–Voigt viscoelastic model is used to show radial displacement evolution along time. The cylindrical vessel (Figure27) with thin thickness is submitted to uniform pressure.Figure 27
Cylindrical pressure vessel with thin thickness.The geometrical properties of the vessel are: 2L = 600 mm, R = 300 mm, and h = 30 mm. The total pressure is equal to P = 260 kPa, totally applied from the initial time step. Regarding temporal discretization, four different time steps were used (0.5 s, 2.5 s, 5.0 s, and 7.5 s), and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions). The Kelvin–Voigt viscoelastic properties are: E = 200 MPa and η = 2,000 MPa s (tε = 10 s). The viscoelastic analysis allows to describe inflation process of the vessel. The parameters were chosen to compare results with the reference Mesquita and Coda [41].Figure28 presents the considered mesh (10 plane beams) used to model one-eighth of the cylindrical mesh, applied pressure, and boundary conditions considering problem symmetry. In original reference [41], a bidimensional mesh with 200 free formulation discrete Kirchhoff triangle finite elements and the viscoelastic Kelvin–Voigt model were considered with time steps equal to 0.5 s.Figure 28
Finite element mesh and initial and boundary conditions used to model one-eighth of the vessel.Numerical results for radial displacements in the middle point of the arch are presented in Figure29, during the inflation process. The results are in good agreement with numerical results available in the literature [41].Figure 29
Radial displacements over time during the inflation process.


## 7.1. Analysis of a Tensioned Viscoelastic Bar

A classical example is presented here. This example is used by different authors to calibrate viscoelastic models because it is a typical case of a plane state of stress. The example deals with a tensioned bar supported on the left and bottom faces. Other faces are free to move, and a traction force is applied on the free end, on the right, as shown in Figure18.Figure 18
Tensioned bar.The bar has a length (L) of 800 mm, a cross-section height (h) of 100 mm, and unitary width (b). The traction force is equal to 0.005 kN/mm2.The Kelvin–Voigt viscoelastic model is considered, with the modulus of elasticity (E) equal to 11 kN/mm2 and modulus of viscosity (η) equal to 500 kN/mm2·day. Thus, the retardation time (tε) is equal to 45.4545 days. Regarding spatial discretization, 10 finite elements are adopted, with 10 Gauss points along the height and 10 Gauss points along the length. Regarding temporal discretization, 5 different time steps were used (1 day, 5 days, 10 days, 25 days, and 50 days). No inertial effects were considered, and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions).Numerical results for axial displacements at the free end are presented in Figure19. The presented analytical results and numerical results are also available in Mesquita and Coda [38]. The analytical results are obtained from equation (20), and the results available in Mesquita and Coda [38] are achieved with a mesh of 8 × 4 cubic boundary elements (boundary element method).Figure 19
Axial displacements over time.From Figure19, it is possible to note that a reduction of time step makes numerical results converge to the analytical solution. This is due to the fact that the smaller the time steps, the more refined the temporal discretization becomes. Therefore, the mechanical behavior is consistent with the analytical solution and in agreement with the results available in Mesquita and Coda [38]. Besides, it is important to note that the convergence for time steps smaller than retardation time is only possible using the proposed procedure. Therefore, these results validate the proposed procedure to avoid numerical divergence for time steps smaller than the retardation time.Figure20 presents numerical results related to the creep-recovery process, considering a time step equal to 1 day. For the time greater than 200 days, the external applied force (P) is reduced to zero. It is possible to note that the obtained results for recovery analysis are consistent and in agreement with the results available in Mesquita and Coda [38].Figure 20
Creep-recovery process.


## 7.2. Cantilever Beam Loaded by Shear Force

This example deals with an analysis of a cantilever beam loaded by a shear force applied on the free end, as shown in Figure21.Figure 21
Cantilever beam.Geometry, rheological parameters, and spatial discretization are the same as the ones presented in the previous example. The shear force is equal to 0.005 kN/mm2, applied at the initial time, and removed after 453 days. Thereafter, the analysis continues for 267 days without applied forces in order to observe the viscoelastic recovery process. Regarding temporal discretization, 720 time steps equal to 1 day are considered. No inertial effects were considered, and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions).Numerical results for vertical displacements in the middle point of the right face are presented in Figure22, for the creep-recovery process. The numerical results available in Panagiotopoulos et al. [9] made with a boundary element method are also presented here. The results are in good agreement with the numerical results available in the literature. A portion of the difference between the results is attributed to the effects of shear, not computed by the adopted Bernoulli–Euler kinematics.Figure 22
Vertical displacements of the free end over time in the creep-recovery process.


## 7.3. Supported Beam with Uniformly Distributed Loading

This example deals with an analysis of a double supported beam loaded by a uniformly distributed force applied on the upper face, as shown in Figure23.Figure 23
Supported beam.The beam has the length (L) equal to 10 m, cross-section height (h) equal to 0.5 m, and width (b) equal to 2 m. The uniformly distributed force (P) is equal to 10 N/m, applied at the initial time, and removed after 3 s.The Kelvin–Voigt viscoelastic model is considered, with the modulus of elasticity (E) equal to 98 MN/m2 and modulus of viscosity (η) equal to 27.44 MN/m2·s. Thus, the retardation (tε) is equal to 0.28 s. Regarding spatial discretization, 10 finite elements are adopted, with 10 Gauss points along the height and 10 Gauss points along the length. Regarding temporal discretization, three different time steps were used (0.5 s, 0.1 s, and 0.01 s). No inertial effects were considered, and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions).The numerical results for vertical displacements at the midspan are presented in Figure24. The analytical results and numerical results available in Aköz and Kadioǧlu [2] are also presented here. Results are in good agreement with analytical solutions and numerical results available in the literature. It should be noted that the analytical solution is presented in Chen [1] and given by the following equation:(64)wt=5PL432bh31E1−e−Et/η.Figure 24
Vertical displacements at the midspan over time.It is important to note that the convergence for time steps smaller than retardation time is only possible using the proposed procedure to avoid numerical divergence. Besides, it is possible to note that the reduction of time step makes numerical results converge to analytical solutions. This is due to the fact that the smaller the time steps, the more refined the temporal discretization becomes.


## 7.4. Supported Beam with Uniformly Distributed Loading Applied Gradually

The same double supported beam presented in the previous example (Figure23) is analyzed here. However, in this example, distributed loading is applied gradually, as shown in Figure 25.Figure 25
Distributed loading applied gradually over time.Geometry, rheological parameters, and spatial discretization are the same as the ones presented in the previous example. Regarding temporal discretization, 500 time steps equal to 0.01 s were considered in the analysis. No inertial effects were considered, and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions).The numerical results for displacements at the midspan are presented in Figure26. The numerical results available in Aköz and Kadioǧlu [2] are also presented. Obtained results are in good agreement with the numerical results available in the literature.Figure 26
Vertical displacements over time at the midspan with distributed loading applied gradually.


## 7.5. Cylindrical Pressure Vessel

The last numerical example is an engineering application where the pressure is applied gradually, and the Kelvin–Voigt viscoelastic model is used to show radial displacement evolution along time. The cylindrical vessel (Figure27) with thin thickness is submitted to uniform pressure.Figure 27
Cylindrical pressure vessel with thin thickness.The geometrical properties of the vessel are: 2L = 600 mm, R = 300 mm, and h = 30 mm. The total pressure is equal to P = 260 kPa, totally applied from the initial time step. Regarding temporal discretization, four different time steps were used (0.5 s, 2.5 s, 5.0 s, and 7.5 s), and a numerical tolerance of 10−8 was adopted for calculations (absolute in terms of nodal positions). The Kelvin–Voigt viscoelastic properties are: E = 200 MPa and η = 2,000 MPa s (tε = 10 s). The viscoelastic analysis allows to describe inflation process of the vessel. The parameters were chosen to compare results with the reference Mesquita and Coda [41].Figure28 presents the considered mesh (10 plane beams) used to model one-eighth of the cylindrical mesh, applied pressure, and boundary conditions considering problem symmetry. In original reference [41], a bidimensional mesh with 200 free formulation discrete Kirchhoff triangle finite elements and the viscoelastic Kelvin–Voigt model were considered with time steps equal to 0.5 s.Figure 28
Finite element mesh and initial and boundary conditions used to model one-eighth of the vessel.Numerical results for radial displacements in the middle point of the arch are presented in Figure29, during the inflation process. The results are in good agreement with numerical results available in the literature [41].Figure 29
Radial displacements over time during the inflation process.


## 8. Conclusion

This paper presents a positional formulation of the finite element method to describe the viscoelastic mechanical behavior in structures discretized by plane frame elements with Bernoulli–Euler kinematics.To assess the effects of viscoelastic behavior, the stress-strain relation deduced from the Kelvin–Voigt rheological model is used. Similar developments can be made considering other rheological models, such as the Boltzmann model and the Zener model. To do so, it is enough to develop the appropriate stress-strain relations and, based on these, particularize the expressions to evaluate the strain energy and its derivatives.From the analysis and results presented in Section6, it is possible to observe the influence of the adopted time step in the iterative process of the used positional formulation. In addition, it is possible to observe the problem of divergence in the iterative process when adopting time steps smaller than the retardation time (characteristic of rheological models). Such a problem makes it impossible to determine the equilibrium positions. To avoid this problem, a penalization procedure of the nodal position correction vector based on the retardation time is proposed.From the results obtained in the analyzes and the presented examples, it is possible to conclude that the developed formulation is consistent and that the proposed procedure allows for obtaining the equilibrium positions for any time step value adopted without presenting divergence problems in the iterative process and without changing the final results of the analyzes.



---

*Source: 2684127-2021-10-04.xml*