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.
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:where U represents total strain energy and W represents the potential energy of external applied forces, given by
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:
Considering the variations regarding the degrees of freedom (nodal parameters in positional formulation), one has
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:
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:
Expanding equation (6) through a first-order Taylor series, one haswhere 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:leading to
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) is taken as an undeformed configuration(2)The components of the vector of residues are calculated using equation (8)(3)The components of the Hessian matrix are calculated using equation (9)(4)The components of the nodal position correction vector 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 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 ()(ii)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 adopted is based on a relation between Euclidian norms as follows: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 energy and its derivatives and 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 Figure 2, 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.
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: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:where E represents the modulus of elasticity and represents the modulus of viscosity.
From equations (12) and (13) and considering stress-strain relations, one haswhich 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 , one has
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 () with a particular solution () as follows:where represents a constant to be evaluated and represents time variable.
Considering an undeformed configuration, the model initial condition is given by the following equation:
Thus, equation (17) can be rewritten as follows:
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:
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 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 and reaches zero over a sufficiently long time.
As shown in Figure 3, 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 (). Thus, the retardation time can be evaluated as follows:
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 equal to the retardation time , one haswhere 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 , one has
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 (varying from 0 to 1) and (varying from −1 to 1), as illustrated in Figure 4.
In the undeformed configuration, a generic point 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 , 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:where and 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 and represent the unit vectors that define the inclinations of the cross-sections, respectively, in undeformed and deformed configurations, as follows:
Thus, from equations (25)–(28) and Figure 4, the and coordinates of a generic point p, in the undeformed configuration, can be expressed, respectively, as follows:where and 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, the X and Y coordinates of a generic point P, in deformed configuration, can be expressed, respectively, as follows:where and 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 Figure 4.
As shown in Figure 4, 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, , , and represent the nodal parameters of node 1, and , , and represent the nodal parameters of node 2 in the undeformed configuration. On the other hand, , , and represent the nodal parameters of node 1, and , , and 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 and the axis and a cubic relation between and the axis. Thus, the coordinates of a point can be obtained, respectively, as follows:where parameters c, d, e, and f can be evaluated from the element boundary conditions in auxiliary configuration, as follows:
From equations (37) and (38), one has
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 , can be obtained as follows:
Nodal parameters presented in equations (33)–(41) are related to deformed configuration. The same equations are valid to undeformed configuration, just considering nodal parameters , , , , , and instead of , , , , , and , respectively. Furthermore, the same consideration is valid to obtain coordinates and and rotation . 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 length 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 and , respectively.
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:
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:
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:where represents the exact neutral axis curvature in deformation configuration, represented by the following equation:
From equations (43)–(45), it is possible to calculate the normal strain in a fiber, considering the nondimensional coordinates ( and ) and the nodal parameters (, , , , , and ).
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:
The chain rule can be used to exchange integration variables as follows:
Thus, the total strain energy can be rewritten as follows:
Performing the first derivative of total strain energy regarding nodal parameters, one has
Performing the second derivative of total strain energy regarding nodal parameters, one has
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:where indexes “” and “” represent current and previous instant of time, respectively, while represents an imposed time step.
Therefore, the stress-strain relation of the Kelvin–Voigt model can be expressed as follows:and total strain energy is represented by the following equation:
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:where represents position rate, evaluated using finite differences as follows:
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 (Figure 6) are presented. The Kelvin–Voigt viscoelastic model is considered, with the modulus of elasticity equal to and modulus of viscosity . A tension load is applied at the free end.
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, Figures 7–10 present the evolution of the residues vector and the evolution of the nodal position correction vector 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 is not presented because its value remains constant and equal to 10,000 MN/m.
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 Figures 7 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. Figure 11 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).
From Figures 11 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 can be expressed as follows:where represents the final viscoelastic strain of the creep analysis and represents the cross-sectional area of the bar.
Thus, equation (7) can be rewritten as follows:
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
Considering the retardation time, equation (58) can be rewritten as follows:
Considering the strain rate approximate by finite differences, one has
Furthermore, it is important to note that the time step can be expressed as a fraction of the retardation time as follows:where is a retardation factor.
Therefore, equation (60) can be rewritten as follows:
Considering the tensioned bar shown in Figure 6, the first derivative of strain is given by , and the second derivative of strain is given by . Thus, after evaluating the integrals, equation (62) can be rewritten as follows:
In the first iteration, . 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 () 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 would be a feasible strain, as shown in Figure 13. Thus, multiplying both sides of equation (63) by the retardation factor (), one has a nodal position correction vector penalized and compatible with a feasible 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 Figure 14. 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.
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.
From Table 2, 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 equal to 1 minus one small residue, such as .
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 , and it can assume a maximum value equal to 1 (in which case, ). 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 (), 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 Figure 16.
Additionally, from Figure 12, 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.
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 Figure 16.
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 Figure 18.
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 () 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 Figure 19. 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).
From Figure 19, 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.
Figure 20 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].
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 Figure 21.
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 Figure 22, 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.
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 Figure 23.
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 () 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 Figure 24. 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:
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 (Figure 23) is analyzed here. However, in this example, distributed loading is applied gradually, as shown in Figure 25.
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 Figure 26. 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.
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 (Figure 27) with thin thickness is submitted to uniform pressure.
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 ( = 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].
Figure 28 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.
Numerical results for radial displacements in the middle point of the arch are presented in Figure 29, during the inflation process. The results are in good agreement with numerical results available in the literature [41].
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 Section 6, 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.
Data Availability
The data used and analyzed in this paper are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
The authors would like to thank Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) for their financial support.