Modeling the creep behavior of GRFP truss structures with Positional Finite Element Method

Abstract

This paper presents the development of a formulation, based on Positional Finite Element Method, to describe the viscoelastic mechanical behavior of space trusses. The numerical method used was chosen due to its efficiency in the applications concerning nonlinear numerical analyses. The formulation describes the positional variation over time under constant stress state (creep). The objective is to provide a way to quantify the creep behavior for space truss structures and thus contribute to the encouragement of GFRP usage in such structural components. Time-dependent behavior of such materials is one the most important factors for their use in design of structures, demanding studies about the deformations expected within the operational life of the structural systems. To perform this study, the proposed methodology considers a standard solid rheological model to describe stress-strain time-dependent law. This model is implemented in the formulation for quantify the total strain energy. The effects of the model parameters in the mechanical response of the structure with accentuated geometric nonlinearity were presented. In this analysis, it was possible to identify the influence of the elastic and the viscous moduli on the creep response. Model calibration was performed using test data obtained from literature and a GFRP transmission line tower cross-arm was simulated to predict the evolution of displacements under real operational loads. From the results, it was possible to observe a fast evolution of displacements due to the creep effect in the first 7,500 h. This increase was close to 0.6% in relation to the displacement obtained in the elastic behavior. The presented methodology provided a simple and efficient way to quantify the creep phenomenon in viscoelastic GFRP composites truss structures, as can be seen in the developed analyses.

Keywords
Creep; Positional Finite Element Method; Rheological Model; Nonlinear Analysis; GFRP

1 INTRODUCTION

Composites have great potential to be employed in the field of Civil, Mechanical, Maritime and Aeronautical Engineering. Specific characteristics, such as high stiffness and strength, associated with low specific weight make their use very attractive, compensating the higher costs of production.

Glass Fiber Reinforced Plastics (GFRP) have been widely used in construction of structures, replacing the usual steel elements, particularly in truss structures such as electrical transmission lattice towers. These materials are essentially composed of glass-fibers embedded in a resin matrix polymer. GFRP prismatic components can be manufactured with different cross sections via the process of pultrusion, in which the fibers are wetted in a viscoelastic matrix (resin) and subsequently pulled through a die for compacting and curing.

Specific properties of these materials make them advantageous to be considered as an alternative to the usual metallic materials. For instance, proprieties such as electrical and magnetic insulation, controllable thermal expansion, fatigue strength, damping characteristics, high strength-to-weight ratio and adequate tensile and compression strengths make GFRP a competitive material to be used in some applications to replace the usual steel components ( Benmokrane et al., 1995 ). Furthermore, GFRP can be easily subjected to recycling: the waste could be incorporated, for example, into based mortars, as sand aggregates and filler replacements, which is a benefit to the environment and also improves the mechanical properties of the host material ( Meira Castro et al. 2013 ).

These characteristics make GFRP an excellent candidate for Transmission Line Towers (TLT). Selvaraj et al. (2012) and Selvaraj et al. (2013) developed experimental studies on X-braced panels made from GFRP pultruded sections and also on a composite cross arm for TLTs, encouraging the use of GFRP structural profiles as an alternative material for steel members. The use of GFRP demonstrated to be efficient in reducing the hallway practiced in India for power utilities. Godat et al. (2013) also investigated the replacement of traditional materials (steel, wood and concrete) in TLTs by fiberglass pultruded members, leading to a better understanding of the behavior of GFRP in such structures. They analyzed the behavior of GFRP angle, square, rectangular, W and I sections under axial load. According Godat et al. (2013) , the generalized buckling appears at low stresses and prevents FRP (Fiber Reinforced Plastic) profiles to reach their full-strength capacity. Izumi et al. (2000) developed a type FRP insulation arm for a 154 kV line post and tested it for mechanical and electrical performance. Yeh and Yang (1997) and Yeh and Yeh (2001) studied the feasibility of building a transmission tower from a composite material and indicated the need to look into the important aspects of creep.

One of the challenges faced in utilizing the GFRP designs is to evaluate the strains of pultruded elements under constant loading during the life of the structure. This time-dependent phenomenon, known as creep, is not simple to characterize and to quantify. The viscoelastic properties are of extreme importance, especially in the study of vibration damping ( Melo and Radford, 2003 ), displacements and failure ( Sá et al., 2011a ; Sá et al., 2011b ) of structures. Numerical techniques and analytical models have been employed, since experimental tests require long periods and are therefore costly.

Regarding the creep phenomenon, according to Benmokrane et al. (1995) and Ascione et al. (2012) , in the creep behavior of GFRP strains due to viscoelastic effects are mainly caused by the resin, i.e., strains related to glass fibers are negligible, and further, the volume and the orientation of fibers have a great influence.

Some examples of techniques applying rheological formulations implemented in numerical models can be found in the specialized literature to predict the time-dependent behavior of these materials. Sá et al. (2011a) and Sá et al. (2011b) proposed a formulation considering the Maxwell rheological model and two rheological Kelvin models connected in series (Bruger-Kelvin model) to express the viscoelastic behavior of an I-section beam. This model, when compared with the experimental data, presented a similar behavior up to the first 1,000 h. Argyris et al. (1991) and Argyris et al. (1992) developed an appropriate numerical model using rheological models of Maxwell and Kelvin-Voigt for the analysis of membrane structures in PVC-coated fabric. Fritsch et al. (2009) developed a novel rheological material model that features a decomposition of the stress into a time independent quasi-static component and a time and strain dependent viscous component, correctly reproducing the quasi-static and dynamic stress-strain behavior of the fiber-reinforced polypropylene. Kaliske (2000) introduced an anisotropic constitutive formulation discretized into finite elements, using the generalized Maxwell model for application to transversely isotropic materials, performing static and dynamic analyses of U and laminates profiles, and extending the application for mechanical biology. Ascione et al. (2012) presented a program of creep tests to validate a mechanical model (based on Maxwell and Kelvin-Voigt rheological models) formulated for the analysis of viscous properties of FRP laminates.

Other recent works related to viscoelastic behavior, focused on the development of numerical formulations, adopt different rheological models and confirm the relevance and actuality of the theme.

Panagiotopoulos et al. (2014) use a formulation based on the Boundary Element Method for the analysis of three-dimensional solids to compare the responses of the Hooke, Kelvin-Voigt and Boltzmann solid models, which respectively represent the instantaneous elastic, the damped elastic and the creep behaviors. Besides these, the responses of the fluid models of Maxwell, Jeffreys and Burgers were also implemented and compared.

Kühl et al. (2016) propose a formulation based on integral equations to describe the nonlinear viscoelastic behavior of High Density Polyethylene (HDP). For the viscoelastic part of the strain, a generalized model of Kelvin-Voigt is adopted, and for the viscoelastoplastic part the equation of Zapas-Crissman was considered. The viscoelastic parameters are obtained through a curve fitting based on the method of Particle Swarm Optimization, whereas the viscoplastic parameters are obtained through a linear regression of the Least Squares Method.

Pascon and Coda (2017) propose an alternative solid finite formulation for large deformations analysis of viscoelastic materials. In this alternative methodology the neo-Hookean hyperelastic law is taken into account together with the standard solid rheological model. This new approach is based on positions and can be used to reproduce creep, stress relaxation and viscoelastic rate dependent stiffening at large strains, which are usually observed in polymeric materials.

Furthermore, some recent studies about viscoelasticity present developments related to the use of fractional derivatives to obtain simpler rheological models and with a more precise representation of the complex behavior of real materials, such as composite materials, polymeric materials, biological tissues, among others ( Bahraini et al., 2013 ; Pérez Zerpa et al., 2015 ; Ciniello et al., 2017 ; Borges et al., 2017 ).

The present paper proposes a model for the analyses of creep phenomenon in viscoelastic composites truss structures. For this purpose, a positional formulation of geometrical nonlinear analysis is used to introduce the rheological model and describe the behavior of viscoelastic materials in structures with time. The objective is to provide a simple and efficient way to quantify the creep behavior for space truss structures and thus contribute to the encouragement of the use of GFRP in such structural components. As an application example, it is presented the calibration of the model based on experimental tests performed by Youssef (2010) at different stress levels and the simulation of the time-dependent deformation evolution resulting from the viscoelastic effect for a typical structure of a power TLT cross-arm.

2 VISCOELASTIC COMPOSITE MATERIALS AND CREEP PHENOMENON

In general, polymeric matrix composites present viscoelastic behavior. This time-dependent behavior of such materials is one the most important factors for their use in design of structures, demanding studies about the deformations expected within the operational life of the structural systems. It is important to notice that fiber-matrix slippage may occur due the viscosity of the resins employed, causing additional plastic deformation over time ( Scott et al., 1995 ). This viscoelastic behavior is a consequence of the Second Thermodynamics’ Law according to which a portion of the imparted energy of deformation is always dissipated as heat by viscous forces even while the rest is resiliently stored. Thus, once the material is loaded for a defined finite time, there is no return to its original position and permanent deformations occur. The process is neither instantaneous nor infinitely slow: it is a gradual process ( Christensen, 1982 ). This phenomenon depends on the duration of loading and the environment conditions that structures are subjected, affecting their functionality and durability. Manufacturing conditions and exposition to temperature changes, pressure and humidity are examples that can affect the viscoelastic behavior of materials. The loading conditions can also interfere in the proprieties, and high stresses levels can lead to rupture by creep over time ( Figure 1 ).

Figure 1
Effect of loading on the creep test.

The viscoelastic behavior sits between two extremes: purely elastic materials and purely viscous materials. However, for the developed model in this work, the contribution of the material plastic deformation due to the viscosity of it is not taken into account. All the strain attained in the loading procedure can be slowly recovered after the structural unloading.

Macro-mechanical theories concerning the viscoelastic behavior of FRP composites have attempted to predict the behavior of these materials. For example, two theories are already consolidated: Findley power law, that describe the time-dependent behavior of FRP composites under a constant stress ( Findley et al., 1976 ), and Schapery single integral equation ( Schapery, 1969 ), that uses the principles of thermodynamics to describe the viscoelastic response of many polymers by four stress-dependent parameters. A thorough review on the creep phenomenon is performed in Scott et al. (1995) .

Sá et al. (2011a) and Sá et al. (2011b) reviewed the creep phenomenon and showed that most of the expressions obtained in studies by other authors have indicated reductions of viscoelastic stiffness after 30 years up to 68% for tensile loading and of 36% for shear loading, when compared to the corresponding elastic values. Creep deflection increased 50% over the same period. This study also showed rapid evolutions of the deformations especially in the first 24 h, which was noted by the premature failure of specimens loaded at 50% of their ultimate strength. Therefore, it is possible to observe that there is a need to study this phenomenon in more detail and to consider the effect of creep in the design of composites structural components.

3 NONLINEAR POSITIONAL FORMULATION FOR VISCOELASTIC SPACE TRUSS ANALYSIS

The positional formulation of the Finite Element Method used for this paper is based on the formulations found in Coda and Greco (2004) and Greco et al. (2006) . The method was selected due to its efficiency in the applications and developments shown in recent published work concerning nonlinear numerical analysis ( Carrazedo and Coda, 2010 ; Jian et al., 2012 ; Ma et al., 2012 ). It is supported by physics principles based on energy equilibrium that are both consistent and easy to understand. This formulation considers the nodal positions of a structure as variables, rather than nodal displacements, regarding a system of reference fixed in the space in order to describe the kinematics of the finite elements. The position description makes use of an intermediate non-dimensional space that enables the determination of a non-linear “engineering” strain calculated from a relative unitary length for different positions in the structure ( Coda and Greco, 2004 ).

A geometric nonlinear formulation for the analysis of truss structures using spatial bar finite elements is further described. The kinematics of this adopted element can be seen in Figure 2 .

Figure 2
Parametrization of a bar element for trusses geometric configuration (Ω ₀ is the initial configuration and Ω is the deformed configuration).

This kinematics can be parameterized as a function of a non-dimensional variable ξ , that assumes values between 0 and 1. Note that X₁ , Y₁ and Z₁ are the coordinates of the initial node and X₂, Y₂ and Z₂ the coordinates of final node of the bar finite element ( Figure 2 ) and represent its six degrees of freedom. It is possible to write the parametric equations of the element as:

$x=X_1+\left(X_2-X_1\right)\xi$ (1)

$y=Y_1+\left(Y_2-Y_1\right)\xi$ (2)

$z=Z_1+\left(Z_2-Z_1\right)\xi$ (3)

Assuming ds as an infinitesimal linear dimension calculated by:

$ds=\sqrt{d{x}^2+d{y}^2+d{z}^2}$ (4)

the differential of ds in terms of ξ can be written as:

$\frac{ds}{d\xi }=\sqrt{{\left(\frac{dx}{d\xi }\right)}^2+{\left(\frac{dy}{d\xi }\right)}^2+{\left(\frac{dz}{d\xi }\right)}^2}=\sqrt{{\left(X_2-X_1\right)}^2+{\left(Y_2-Y_1\right)}^2+{\left(Z_2-Z_1\right)}^2}=l$ (5)

that determines the length of the element (l).

By varying the length of the element caused by the deformation of the structure it is possible to determine the Total Strain Energy U accumulated for each structural element. The term U is described by the integral of the deformation energy in a volume element, the Specific Strain Energy (u), which considers viscoelastic effects:

$U={\displaystyle \underset{V}{\int }udV=}{\displaystyle \underset{V}{\int }{\displaystyle \underset{\epsilon }{\int }\sigma d\epsilon }}dV$ (6)

where σ is defined as the Cauchy stress and ε is the engineering strain measure. The strain energy is assumed to be zero in a reference position or non-deformed position of the structure ( Coda and Greco, 2004 ).

3.1 The standard solid rheological model

It was considered the standard solid rheological model to take into account the viscoelastic behavior of the material. This model is formed by the association of two elastic elements and one viscous element ( Marques and Creus, 2012 ), as shown in Figure 3 . The model describes the rheological law which relates stress, strain and time.

Figure 3
Standard solid rheological model (uniaxial representation).

From Figure 3 it is possible to assume the following relations ( Equations (7) and (8) ) for the total stress and strain, respectively:

$\sigma ={\sigma }_1+{\sigma }_2$ (7)

$\epsilon ={\epsilon }_1={\epsilon }_2={\epsilon }_2^e+{\epsilon }_2^v$ (8)

The subscripts ()₁ and ()₂ refer to the path 1 (elastic element) and path 2 (elastic and viscous elements) of the model ( Figure 3 ) and the superscripts ()^e and ()^v refer to the elastic and viscous components, respectively.

The required terms of Equations (7) and (8) are given by:

${\sigma }_1=E_1\epsilon$ (9)

${\sigma }_2={\sigma }_2^e={\sigma }_2^v$ (10)

${\epsilon }_2^e=\frac{{\sigma }_2^e}{E_2}$ (11)

${\dot{\epsilon }}_2^v=\frac{{\sigma }_2^v}{{\eta }_2}$ (12)

where E₁, E₂ and η are physical parameters of elements presented in Figure 3 and represent the mechanical properties of the material.

Deriving the Equation (8) as a function of time, combining the terms obtained by Equations (11) and (12) , and considering the time-independent elastic modulus, results in:

$\dot{\epsilon }=\frac{{\dot{\sigma }}_2}{E_2}+\frac{{\sigma }_2}{{\eta }_2}$ (13)

which $\dot{\epsilon }$ represents the deformation rate of the model. Replacing the Equations (7) and (9) in Equation (13):

$\dot{\epsilon }=\frac{\dot{\sigma }-E_1\dot{\epsilon }}{E_2}+\frac{\sigma -E_1\epsilon }{{\eta }_2}$ (14)

and rearranging the terms of Equation (14) for σ, follows:

$\sigma =E_1\epsilon +\frac{{\eta }_2(E_1+E_2)\dot{\epsilon }}{E_2}-\frac{{\eta }_2\dot{\sigma }}{E_2}$ (15)

Equation (15) represents the rheological relation, which describes viscoelastic mechanical behavior of the model.

Considering the used nonlinear formulation, the total strain ε can be written as follows:

$\epsilon =\left(\frac{\sqrt{B}}{l_0}-1\right)$ (16)

where l₀ is the initial length of the element, and B is defined as:

$B={\left(X_2-X_1\right)}^2+{\left(Y_2-Y_1\right)}^2+{\left(Z_2-Z_1\right)}^2$ (17)

It is important to note that the positional formulation considers the nodal positions as variables of the problem while the rheological relation is a function of the deformation rate. Thus, by using the chain rule, the strain rate $\dot{\epsilon }$ of the structure is replaced by the product between the derivative of the strain with respect to the nodal parameters of the positional formulation and the rate of variation of these parameters, as follows:

$\dot{\epsilon }=\frac{\partial \epsilon }{\partial x}\frac{\partial x_i}{\partial t}=\epsilon {,}_i{\dot{x}}_i$ (18)

where ${\dot{x}}_i$ represents the velocity of positional variation of the finite element in each degree of freedom. An auxiliary variable, ${\dot{x}}_{axial}$ , can be used to describe the length variation of the finite element at a time variation ∆ t, as described by the following expression:

$\begin{array}{l}{\dot{x}}_{axial}=\frac{1}{\Delta t}(\sqrt{{\left(X_2^{t=1}-X_1^{t=1}\right)}^2+{\left(Y_2^{t=1}-Y_1^{t=1}\right)}^2+{\left(Z_2^{t=1}-Z_1^{t=1}\right)}^2}+\\ -\sqrt{{\left(X_2^{t=0}-X_1^{t=0}\right)}^2+{\left(Y_2^{t=0}-Y_1^{t=0}\right)}^2+{\left(Z_2^{t=0}-Z_1^{t=0}\right)}^2}\end{array}$ (19)

where $X_1^{t=1}$ , $Y_1^{t=1}$ , $Z_1^{t=1}$ , $X_2^{t=1}$ , $Y_2^{t=1}$ , and $Z_2^{t=1}$ are the positions in the current instant moment and $X_1^{t=0}$ , $Y_1^{t=0}$ , $Z_1^{t=0}$ , $X_2^{t=0}$ , $Y_2^{t=0}$ , and $Z_2^{t=0}$ are the positions in the previous instant. Thus, the terms of the velocity of positional variation in each degree of freedom can be obtained by means of an appropriate matrix of transformation from the local (axial) system to the global system.

3.2 The deformation described in terms of time variable

The rheological model for the viscosity analysis of structural systems and the development of the positional formulation for the deformation described in terms of time variable is following presented. Rewriting Equation (6) using Equation (15):

$U={\displaystyle \underset{V}{\int }\left({\displaystyle \underset{\epsilon }{\int }{E}_1\epsilon d\epsilon +{\displaystyle \underset{\epsilon }{\int }\frac{\eta (E_1+E_2)\dot{\epsilon }}{E_2}d\epsilon -{\displaystyle \underset{\epsilon }{\int }\frac{\eta \dot{\sigma }}{E_2}d\epsilon }}}\right)}dV$ (20)

In this equation, the strain terms must be integrated in the volume element. Substituting Equation (16) in Equation (20) and writing dε in terms of dX_i results in:

$U={\displaystyle \underset{V}{\int }(\frac{E_1}{2}{\left(\frac{\sqrt{B}}{l_0}-1\right)}^2}+{\displaystyle \underset{X_i}{\int }\frac{\eta \left(E_1+E_2\right)}{4l_0^2{E}_2}\frac{B{,}_i^2}{B}{\dot{x}}_{i}d{X}_i}-\eta {\dot{\epsilon }}_2^e\left(\frac{\sqrt{B}}{l_0}-1\right)-{\displaystyle \underset{X_i}{\int }\frac{\eta E_1}{4l_0^2{E}_2}\frac{B{,}_i^2}{B}{\dot{x}}_{i}d{X}_i})dV$ (21)

Considering the cross-section area (A) invariant, and a simple mapping in the dimensionless variable ξ domain, it is possible to achieve the integration of Equation (21) along the length of the finite element, as follows:

$U=l_0{\displaystyle \underset{\xi }{\int }A}(\frac{E_1}{2}{\left(\frac{\sqrt{B}}{l_0}-1\right)}^2+{\displaystyle \underset{X_i}{\int }\frac{\eta \left(E_1+E_2\right)}{4l_0^2{E}_2}\frac{B{,}_i^2}{B}{\dot{x}}_{i}d{X}_i}-\eta {\dot{\epsilon }}_2^e\left(\frac{\sqrt{B}}{l_0}-1\right)-{\displaystyle \underset{X_i}{\int }\frac{\eta E_1}{4l_0^2{E}_2}\frac{B{,}_i^2}{B}{\dot{x}}_{i}d{X}_i})d\xi$ (22)

Equation (22) can be rewritten in a compact form;

$U={\displaystyle {\int }_0^1{l}_0{u}_l}d\xi$ (23)

where the variable u_l is a term associated with the specific energy of deformation per unit of length.

For a conservative structural problem associated with a reference system fixed in space, it is possible to write the stationary Total Potential Energy (П ) as a function of the Total Strain Energy (U) and the Potential Energy of Applied Forces (P) acting on each node of a finite element, as follows:

$\Pi =U-P$ (24)

The principle of Minimum Potential Energy is applied, in which the equilibrium of the structure will occur when the differential of the Total Potential Energy (П ) for each degree of freedom of the structure is equal to 0, i.e., when the change rate of the Total Potential Energy is equal to zero. Writing in indicial notation (i = 1, 2, ..., 6):

$\frac{\partial \Pi }{\partial X_i}=l_o{\displaystyle \underset{0}{\overset{1}{\int }}\frac{\partial u_l}{\partial X_i}}d\xi -F_i=0$ (25)

The Equation (25) represents a system with six equations whose adopted solution consists in a numerical strategy in which the derivatives are developed within the integral and analytically integrated for the dimensionless variable ξ. It can be noticed that the resulting numerical integral is nonlinear for the nodal positions and corresponds to a system of six equations by finite element due to its the six degrees of freedom ( Greco et al., 2006 ). Therefore, the system of Equation (25) can be written using indicial notation (free index i = 1 to 6 and dummy index j = 1 to 6):

$\frac{\partial \Pi }{\partial X_i}=g_i(X_j,F_i)=f_i(X_j)-F_i=0$ (26)

or in a tensor representation:

$g_i(X)=f_i-F_i=0$ (27)

Note that, in this study, the applied forces are independent of coordinates (i.e. conservative forces). The vector function g_i(X) is nonlinear regarding the nodal parameters. To solve Equation (27) , the Newton-Raphson procedure can be used ( Kleiber, 1989 ):

$g_i\left(X\right)\widetilde{=}{g}_i\left(X_0\right)+g_i{,}_k\left(X_0\right)\Delta X\widetilde{=}0$ (28)

where X is any nodal position and X₀ is the initial nodal position.

The iterative Newton-Raphson process is summarized as follows ( Greco et al., 2006 ):

1. 1

   Assume X₀ as the initial configuration. Calculate g_i(X₀) by Equation (29):

   $g_i(X_0)=l_0{\displaystyle \underset{0}{\overset{1}{\int }}{u}_l{,}_i}d\xi |{}_{X_0}-F_i$(29)

2. 2

   For this X₀, calculate the Hessian matrix. Integrate this value, as indicated in Equation (30) , and the result is the gradient of g_i at X ₀:

   $g_i{,}_k(X_0)=l_0{\displaystyle \underset{0}{\overset{1}{\int }}{u}_l{,}_{ik}}d\xi |{}_{X_0}$(30)

3. 3

   Solve the system of Equation (28) and determine ∆X.

4. 4

   Update position X₀ = X₀ + ΔX. Return to step (1) until ΔX to be less than the tolerance.

3.3 Necessary algebraic development

To implement this formulation, the variables involved must be numerically determined, thus, the term l₀u_l of Equation (23) is rewritten as:

$l_0{u}_l=\frac{E_{1}A{l}_0}{2}{\left(\frac{\sqrt{B}}{l_0}-1\right)}^2+{\displaystyle \underset{X_i}{\int }\frac{\eta A\left(E_1+E_2\right)}{4l_0{E}_2}\frac{B{,}_i^2}{B}{\dot{x}}_{i}d{X}_i}-\eta A{l}_0{\dot{\epsilon }}_2^e\left(\frac{\sqrt{B}}{l_0}-1\right)-{\displaystyle \underset{X_i}{\int }\frac{\eta A{E}_1}{4l_0{E}_2}\frac{B{,}_i^2}{B}}{\dot{x}}_{i}d{X}_i$ (31)

The first derivative of Equation (31) , related to the nodal parameter i, is expressed by:

$l_0{u}_l{,}_i=\frac{E_{1}A}{2}\frac{B{,}_i}{\sqrt{B}}\left(\frac{\sqrt{B}}{l_0}-1\right)+\frac{\eta \left(E_1+E_2\right)}{E_2}\frac{A}{4l_0}\frac{B{,}_i^2{\dot{x}}_i}{B}-\frac{\eta A{\dot{\epsilon }}_2^{e}B{,}_i}{2\sqrt{B}}-\frac{\eta A{E}_{1}B{,}_i^2{\dot{x}}_i}{4l_0{E}_{2}B}$ (32)

Analogously the derivative of Equation (32) related to the nodal parameter k is given by:

$$\begin{array}{l}{l}_0{u}_l{,}_{ik}=\frac{\text{A}}{4B}[\frac{E_1}{\sqrt{B}}\left(2BB{,}_{ik}\left(\frac{\sqrt{B}}{l_o}-1\right)+B{,}_{i}B{,}_k\right)+\frac{\eta (E_1+E_2)B{,}_i}{B{l}_0{E}_2}\left(2BB{,}_{ik}{\dot{x}}_i+BB{,}_i{\dot{x}}_i{,}_k+B{,}_{i}B{,}_k{\dot{x}}_i\right)+\\ +\frac{2\eta }{\sqrt{B}}\left(BB{,}_i{\dot{\epsilon }}_2^e{,}_k+{\dot{\epsilon }}_2^e\left(B{,}_{i}B{,}_k-BB{,}_{ik}\right)\right)]\end{array}$$ (33)

With these equations and the iterative Newton-Raphson process, all of the terms required for obtaining the nodal parameters X₁, Y ₁, Z₁, X₂, Y₂ and Z₂ can be calculated, taking into account the effects of viscoelastic behavior. Thus, it is possible to define the equilibrium positions of a truss structure of GFRP over time by adopting suitable material parameters (E₁, E₂ and η ). These parameters for describing the viscoelastic behavior can be obtained by calibrating the proposed formulation by means of experimental results of traction creep tests of GFRP bars.

It should be noted that no coordinate transformation was performed using this formulation, since no local coordinate systems are used. Therefore, the derivatives are calculated in a single coordinate system. Normal loads acting on the elements can be calculated by the Cauchy equation ( Greco et al., 2006 ).

4 NUMERICAL EXAMPLES AND ANALYSIS OF RESULTS

For the analyses presented, two assumptions were adopted. The first is that the rate of stress $\dot{\sigma }$ in Equation (15) is set equal to zero. This approximation is valid for the study of creep behavior, in which the stress is kept constant over time. For relaxation, this assumption is not valid because, in this case, the deformation is taken as constant and the stress is varying in time. The second refers to the elastic parameters shown in Equation (34). It is assumed that the sum of the modules E₁ and E₂ is always constant and represents the material longitudinal elasticity modulus E.

$E=E_1+E_2$ (34)

4.1 Analysis of viscoelastic parameters

We would like to demonstrate the effects of the viscoelastic rheological model parameters in the mechanical response of the structure with accentuated geometric nonlinearity. The application considered as an example is a typical truss structure ( Figure 4 ) subjected to the instability phenomenon known as snap-through. In this analysis, prescribed downward displacements at the load point were applied varying from 0 to 0.225 m in steps of 0.001 m.

Figure 4
Truss geometry (dimensions in m).

A variation in the elastic modulus E₁ was considered in the first case. Figure 5 presents the results obtained for the labeled node where the prescribed position was applied. It can be noted that as the elastic modulus E₁ decreases, the force required to achieve the same position is smaller. This follows from the simple fact that the effect of this viscoelastic parameter causes a lower structural stiffness and therefore the structure comes into instability to a lower stress level. The elastic curve is represented by E₁ = E, i.e., without considering the viscosity of the material. The other fixed parameters used were A = 0.0009 m², E = 46.9 GPa, η = 10.0 GPa·s and Δt = 100 s.

Figure 5
Influence of the elastic modulus E₁ on the equilibrium path (prescribed position).

In the second case, a variation in the viscous modulus η was considered, while the elastic moduli, E₁ and E₂, were kept respectively fixed as 0.999 and 0.001 of the material longitudinal elasticity modulus E. Furthermore, the other parameters were considered the same as in the previous case. Figure 6 presents the results obtained for the node where the prescribed position is applied. The elastic curve was represented with η tending to infinity, i.e., the dashpot was considered as a rigid bar in the rheological model. Figure 7 is a detail of Figure 6 , magnifying the region of the snap-through first peak. It is possible to note that the variation in the viscous modulus causes slight changes on the obtained results, mainly due to the influence of the viscous modulus on the velocity to obtain the equilibrium position for a given applied load, while the value of the equilibrium position is defined by the elastic modulus E₁.

Figure 6
Influence of the viscous modulus η on the equilibrium path (prescribed position).

Figure 7
Influence of the viscous modulus η on the equilibrium path in detail (prescribed position).

Figure 8 shows the contribution of the viscous modulus to the viscoelastic behavior, in the analysis of the same structure. In this analysis, a downward static load of P = 14.7 kN was applied and maintained constant during 30 time steps of Δt = 100 s. The elastic modules, E₁ and E ₂, were kept respectively 0.9 and 0.1 of the material longitudinal elasticity modulus E. The influence of the viscous modulus on the velocity to obtain the equilibrium position can be clearly seen. The equilibrium position is reached faster for lower viscous moduli materials. Furthermore, it is possible to notice the phenomenon of the snap-through for load levels lower than expected in the elastic snap-through behavior, considering the creep behavior.

Figure 8
Influence of the viscous modulus η on displacement over time.

4.2 Model calibration

In order to access the viscoelastic behavior of GFRP, a model calibration was performed based on the test data obtained by Youssef (2010) . In his work, he performed an experimental study on the creep behavior of GFRP rods for concrete reinforcement over the first 10,000 hours at different stress levels: 15%, 30%, 45% and 60% of ultimate strength. The tested bar is made of high-strength E-glass fibres (77.3% fibres by volume; 59.9% fibres by weight) impregnated in vinylester resin. The bar's circular cross section has 9.5 mm diameter; manufactured by a process that involves pultrusion with sand-coating along the external surface of the bar. The surrounding environment whilst conducting the long-term creep tests is standard laboratory atmosphere (23 ± 3 °C and 50 ± 10% relative humidity).

As presented in Youssef (2010) , all samples were prepared according to CAN/CSA S806-02 ( Canadian Standards Association, 2002 ) and ACI 440.3R-04 ( American Concrete Institute, Committee 440, 2004 ). The GFRP bars were cut into a variety of lengths (1170 mm, 1270 mm and 1470 mm) to fit into three different frame sizes of heights 1550 mm, 1750 mm and 1880 mm, respectively. Each one of the bar sample extremities was fitted into a 410 mm-long steel tube (grip) using Bristar 100 expansive grout. The steel pipe that grips on both ends have a 50 mm hollow portion, threaded on the inside, to screw/fit onto spherical nuts that keep the sample intact with the frame. Two Kyowa 10-mm strain gauges (120 ohm resistance and a gauge factor of 2.10) were attached on opposite sides at the middle of the free portion of each sample. For gauge installation, M-bond AE adhesive was used and gauges were properly aligned in the longitudinal direction of the bar.

The schematic of a sustained-load frame, presented in Youssef (2010) and shown in Figure 9 , illustrates the comprising frame elements and the location of the sample within. The objective of such frames is to maintain a constant tensile load sustained along the bar's length for extended time durations. The load should be maintained perfectly axial; assurance that no eccentricity or bending occurs to the bar is imperative. The associated load magnifying system (the two lever arms and the sustained weight-pan) multiplies the kept-on-pan weight to reach a significant percentage of the sample's ultimate tensile capacity.

Figure 9
Schematic of bar sustained load frame ( Youssef, 2010 ).

Youssef (2010) noticed that the failure of the material by creep in this time range occurred for stresses close to 60% of the tensile strength of the material. The properties of the GFRP samples extracted from Youssef (2010) are shown in Table 1 .

Table 1
Mechanical and geometrical properties ( Youssef, 2010 ).

Ultimate tensile stress	MPa	854
Modulus of elasticity	GPa	46.9
Cross-section área	m2	7.09 · ${10}^{-5}$
Length	m	1

The calibration of the formulation was performed by obtaining the relations between the viscoelastic parameters (E₁, E₂ and η) and the stress level (σ) to which the material was exposed. Thus, suitable parameters were numerically adjusted to approximately match the experimental results obtained by Youssef (2010) . For this, from the suitable viscoelastic parameters for each stress level, a quadratic regression was performed, using the least squares method, obtaining the following relations:

$\eta (\sigma )=\text{1}\text{.4375}\cdot {\text{10}}^{\text{13}}-\text{2}\text{.7127}\cdot {\text{10}}^{\text{4}}\sigma +7.6175\cdot {\text{10}}^{\text{-6}}{\sigma }^2$ (35)

$E_1(\sigma )=(\text{0}\text{.9945}-\text{4}\text{.6972}\cdot {\text{10}}^{\text{-10}}\sigma +5.5717\cdot {\text{10}}^{\text{-19}}{\sigma }^2)E$ (36)

where σ is the acting stress in Pascal (Pa). Figures 10 10(b) show the curves describing the behavior of the parameters η and E₁ as function of the applied stress level σ , described respectively by Equations (35) and (36) .

Figure 10
(a) Relation between the viscous modulus η and the stress level σ for GFRP material; (b) Relation between the relative ratio E1/E and the stress level σ for GFRP material.

The results of the model calibration using Equations (34, 35 and 36)))) are shown in Figure 11 , comparing the different experimental Ultimate Stress levels from Youssef (2010) with the model predictions of the strain time-depending behavior.

Figure 11
Calibration of the evolution of the displacements caused by the creep phenomenon based on Youssef (2010) experimental tests.

Using the calibrated model, it is possible to predict the creep behavior of a GFRP structure for active different stress levels. An example of this application is further presented, showing the evolution of the deformations for a TLT cross-arm.

4.3 Analyzing the viscoelastic behavior of TLT cross-arm structure

The cross-arm analyzed is a typical suspension TLT subjected to the loading hypothesis that considered the line in service without any exceptional loads, only the permanent active operational loads. The goal of this analysis is to verify the increase in displacement due to viscoelastic effects. The geometry of the analyzed structure, presented in Figure 12 , is divided into three groups of bars ( Table 2 ).

Figure 12
Cross-arm transmission line suspension tower (dimensions in m).

Table 2
Structure elements.

Bars	Type of stress	Cross-section area (10-4 m)	Active stress (MPa)
1, 2, 3, 4	tension	4.00	15.9
5, 6, 7, 8	compression	9.48	6.1
9, 10, 11, 12, 13, 14, 15, 16	tension or compression	2.77	< 0.1

The structure was modeled with 16 bar finite elements and 9 nodes, from which 4 are restricted, resulting in 15 degrees of freedom. The load, with magnitude of 5,394 N, is applied instantaneously in the vertical direction (negative Z axis) at the tip of the structure. It corresponds to the dead weight of the conductor cables, insulation chain, and accessories for the class of the TLT considered (138 kV lines). The cross-sections of the bars were statically dimensioned for tensile loading according to the Tsai-Wu criteria ( Tsai and Wu, 1971 ), for local buckling according to the methodology defined by Pecce and Cosenza (2000) and for global buckling according to Engesser’s equation ( Zureick and Scott, 1997 ).

The vertical displacement at the loading node is shown in Figure 13 . To generate this response, it was used Δt = 100 h and the viscoelastic parameters η, E₁ and E₂ obtained by Equations (34, 35 and 36)))). Table 3 shows the value of the displacement evolution for different increasing times.

Figure 13
Effect of loading on the creep test.

Table 3
Displacements after varied times.

	Displacement (mm)	Relative increase in displacement (%)
No creep (t = 0 h)	3.18751	–
Creep (t = 2,400 h)	3.19791	0.3
Creep (t = 7,500 h)	3.20751	0.6
Creep (t = 94,800 h)	3.21620	0.9
Creep (t = 222,100 h)	3.21939	1.0

From Figure 13 and Table 3 , it is possible to observe a fast evolution of vertical displacement due to the creep effect in the first 7,500 h. This increase is close to 0.6% in relation to the displacement obtained in the elastic behavior. After 7,500 h (~ 313 days) the evolution of displacement decreases considerably, reaching 1% only after 25 years, after which the structure is considered stable. More pronounced creep contributions can be obtained for higher stress levels, i.e., for a load of 100 kN (~ 20 times higher) an increase of approximately 9% of the displacement is predicted over 2,000 h (83 days). In this case, the maximum stress is 297 MPa (~ 35% of the material ultimate stress).

The analysis was performed for the first 300,000 h and the results obtained are of considerable relevance for the designers to realistically access the behavior of viscoelastic structural components. However, for considerably long-time interval other aspects such as environmental conditions (moisture, temperature, radiation, etc.) must be introduced into the model to represent the behavior of the material with more reliability. Furthermore, it is important to note that the creep behavior was considered symmetric with respect to the stress level for simplification, i.e., the creep behavior for the tensile rods is equal to the compressed rods.

5 CONCLUSIONS

A formulation of geometrical nonlinear analysis was successfully implemented to analyze the creep phenomenon in viscoelastic GFRP composites truss structures. A standard rheological model was used together to the Positional Finite Element Method to predict the evolution of deformation of viscoelastic structural components over time. The formulation was applied in a relatively simple way, since it is based on physical concepts of Total Potential Energy balance, which make it easy to understand. The viscoelastic behavior was introduced in a natural way through the implementation of a rheological relation for obtaining the Total Deformation Energy, since this relation is based on the physical behavior of springs and dashpots associations that provides equations known as phenomenological.

A parametric study was conducted to identify the influence of the elastic and the viscous moduli on the creep response of a structure with snap-through behavior. From the obtained results, it was concluded that elastic modulus is clear responsible to represent the structural stiffness. Thus, the elastic modulus decrease leads to a reduction in the force required to achieve the same position and therefore the structure presents a snap-through instability with a lower stress level. On the other hand, the variation in the viscous modulus causes slight changes on the obtained equilibrium position. This is due to the influence of the viscous modulus is mainly related to the velocity necessary to obtain the equilibrium position for a given applied load. Thus, the viscous modulus decrease leads to a reduction in the time required to achieve the same equilibrium position and, therefore, the snap-through instability phenomena can be anticipated.

In order to model the viscoelastic behavior of GFRP, a calibration of the model was performed in a simple way based on tensile creep test results obtained from the literature. The results of the calibration demonstrated fairly good approximation between numerical and experimental results.

A GFRP TLT cross-arm structure, subjected to the real operational loading was selected to be analyzed and numerical simulation are performed to evaluate increase vertical tip displacement over time. The model was able to simulate the evolution of displacements throughout the structure and on critical nodes, such as load-application point, which are essential for safe distance restrictions and structural tolerances determinations. From the results, it is observed a fast evolution of the vertical displacements due to the creep effect in the first hours as expected from the theory of viscoelasticity. This increase is close to 0.6% in relation to the displacement obtained in the elastic behavior in the first 7,500h. After 7,500 h the evolution of displacement decreases considerably, reaching 1% only after 25 years, after which the structure is considered stable. It is concluded that the analyzed GFRP TLT cross-arm is stable and present slight creep displacements under real operational loads. However, special caution must be taken, because a pronounced creep contribution was obtained for higher stress levels, even though relatively less than the material ultimate stress.

Acknowledgements

The authors acknowledge the financial supports granted by CNPq (National Council for Scientific and Technological Development), CAPES (Coordination of Improvement of Higher Education Personnel), FAPEMIG (Foundation for Research Support of Minas Gerais State) and Federal University of Minas Gerais (UFMG), under Grant Numbers 302376/2016-0, TEC-PPM-00354-14 and TEC-PPM-00409-16.

• http://dx.doi.org/10.1590/1679-78254432

References

• American Concrete Institute, Committee 440, (2004). Guide Test Methods for Fibre Reinforced Polymers (FRPs) for Reinforcing or Strengthening Concrete Structures. ACI 440.3R-04, American Concrete Institute, Farmington Hills, Mich., 40p.
• Argyris, J., Doltsinis, I.S., da Silva, V.D. (1992). Constitutive modelling and computation of non-linear viscoelastic solids. Part II: Application to orthotropic PVC-coated fabrics. Computer Methods in Applied Mechanics and Engineering, 98(2): 159–226.
• Argyris, J., Doltsinis, I.S., da Silva, V.D. (1991). Constitutive modelling and computation of non-linear viscoelastic solids. Part I: Rheological models and numerical integration techniques. Computer Methods in Applied Mechanics and Engineering, 88(2): 135–163.
• Ascione, L., Berardi, V.P., D’Aponte, A. (2012). Creep phenomena in FRP materials. Mechanics Research Communications, 43: 15–21.
• Bahraini, S.M.S., Eghtesad, M., Farid, M., Ghavanloo, E. (2013). Large deflection of viscoelastic beams using fractional derivative model. Journal of Mechanical Science and Technology, 27(4): 1063-1070.
• Benmokrane, B., Chaallal, O., Masmoudi, R. (1995). Glass fiber reinforced plastic (GFRP) rebars for concrete structures. Construction and Building Materials, 9(6): 353–364.
• Borges, A.S., Faria, A., Rade, D., Sales, T. (2017). Time Domain Modeling and Simulation of Nonlinear Slender Viscoelastic Beams Associating Cosserat Theory and a Fractional Derivative Model. Latin American Journal of Solids and Structures, 14(1): 153-173.
• Canadian Standards Association, (2002). Design and Construction of Building Components with Fibre Rei forced Polymers. CAN/CSA S806-02, Rexdale, Ontario, Canada, 177p.
• Carrazedo, R., Coda, H.B. (2010) Alternative positional FEM applied to thermomechanical impact of truss structures. Finite Elements in Analysis and Design, 46(11) 1008-1016.
• Ciniello, A.P.D., Bavastri, C.A., Pereira, J.T. (2017). Identifying Mechanical Properties of Viscoelastic Materials in Time Domain Using the Fractional Zener Model. Latin American Journal of Solids and Structures, 14(1): 131-152.
• Coda, H.B., Greco, M. (2004). A simple FEM formulation for large deflection 2D frame analysis based on position description. Computer Methods in Applied Mechanics and Engineering, 193(33): 3541–3557.
• Christensen, R.M. (1982) Theory of Viscoelasticity: An Introduction. 2^nd ed., Academic Press (New York).
• Findley, W.N., Lai, J.S., Onaran, K. (1976). Creep and Relaxation of Nonlinear Viscoelastic Materials. Dover Publications (New York).
• Fritsch, J., Hiermaier, S., Strobl, G. (2009). Characterizing and modeling the non-linear viscoelastic tensile deformation of a glass fiber reinforced polypropylene. Composites Science and Technology, 69(14): 2460–2466.
• Godat, A., Légeron, F., Gagné, V., Marmion, B. (2013). Use of FRP pultruded members for electricity transmission towers. Composite Structures, 105: 408–421.
• Greco, M., Gesualdo, F.A.R., Venturini, W.S., Coda, H.B. (2006). Nonlinear positional formulation for space truss analysis. Finite Elements in Analysis and Design, 42(12): 1079–1086.
• Izumi, K., Takahashi, T., Homma, H., Kuroyagi, T. (2000). Development of line post type polymer insulation arm for 154 kV. IEEE Transactions on Power Delivery, 15(4): 1304–1310.
• Jian, L., Shenjie, Z., Meiling, D., Yugin, Y. (2012). Three-node Euler-Bernoulli beam element based on positional FEM. Procedia Engineering, 29: 3703-3707.
• Kaliske, M. (2000). A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains. Computer Methods in Applied Mechanics and Engineering, 185(2): 225–243.
• Kleiber, M. (1989). Incremental Finite Element Modeling in Non-linear Solid Mechanics, Ellis Horwood (England).
• Kühl, A., Muñoz-Rojas, P.A., Barbieri, R., Benvenutti, I.J. (2016). A procedure for modeling the nonlinear viscoelastoplastic creep of HDPE at small strains. Polymer Engineering and Science, 57(2): 144-152.
• Ma, L.Z., Liu, J., Diao, X.L., Yan, Y.Q. (2012). The Nonlinear Static Formulation of Variable Cross-Section Beam Element Based on Positional Description. Advanced Materials Research, 557: 2367-2370.
• Marques, S.P.C., Creus, G.J. (2012). Computational viscoelasticity. Springer (Heidelberg).
• Meira Castro, A.C., Ribeiro, M.C.S., Santos, J., Meixedo, J.P., Silva, F.J., Fiúza, A., Alvim, M.R. (2013). Sustainable waste recycling solution for the glass fiber reinforced polymer composite materials industry. Construction and Building Materials, 45: 87–94.
• Melo, J.D.D., Radford, D.W. (2003). Viscoelastic Characterization of Transversely Isotropic Composite Laminae. Journal of Composite Materials, 37(2): 129–145.
• Panagiotopoulos, C.G., Mantič, V., Roubíček, T. (2014). A simple and efficient BEM implementation of quasistatic linear visco-elasticity. International Journal of Solids and Structures, 51(13): 2261-2271.
• Pascon, J.P., Coda, H.B. (2017). Finite deformation analysis of visco-hyperelastic materials via solid tetrahedral finite elements. Finite Elements in Analysis and Design, 133: 25-41.
• Pecce, M., Cosenza, E. (2000). Local buckling curves for the design of FRP profiles. Thin-Walled Structures, 37(3): 207–222.
• Pérez Zerpa, J., Canelas, A., Sensale, B., Santana, D.B., Armentano, R.L. (2015). Modeling the arterial wall mechanics using a novel high-order viscoelastic fractional element. Applied Mathematical Modelling, 39(16): 4767-4780.
• Sá, M.F., Gomes, A.M., Correia, J.R., Silvestre, N. (2011a). Creep behavior of pultruded GFRP elements - Part 1: Literature review and experimental study. Composite Structures, 93(10): 2450–2459.
• Sá, M.F., Gomes, A.M., Correia, J.R., Silvestre, N. (2011b). Creep behavior of pultruded GFRP elements - Part 2: Analytical study. Composite Structures, 93(9): 2409–2418.
• Schapery, R.A. (1969). On the Characterization of Nonlinear Viscoelastic Materials. Polymer Engineering and Science, 9(4): 295-310.
• Scott, D.W., Lai, J.S., Zureick, A. (1995). Creep Behavior of Fiber-Reinforced Polymeric Composites: A Review of the Technical Literature. Journal of Reinforced Plastics and Composites, 14(6): 588-617.
• Selvaraj, M., Kulkarni, S., Babu, R.R. (2013). Analysis and experimental testing of a built-up composite cross arm in a transmission line tower for mechanical performance. Composite Structures, 96: 1–7.
• Selvaraj, M., Kulkarni, S., Babu, R.R. (2012). Structural evaluation of FRP Pultruded Sections in overhead transmission line Towers. International Journal of Civil Structural Engineering, 2(3): 943–949.
• Tsai, S.W., Wu, E.M. (1971). A General Theory of Strength for Anisotropic Materials. Journal of Composite Materials, 5(1): 58–80.
• Yeh, H.Y., Yang, S.C. (1997). Building of a Composite Transmission Tower. Journal of Reinforced Plastics and Composites, 16(5): 414-424.
• Yeh, H.Y., Yeh, H.L. (2001). A Simple Failure Analysis of the Composite Transmission Tower. Journal of Reinforced Plastics and Composites, 20(12): 1054–1065.
• Youssef, T.A. (2010). Time-Dependent Behavior of Fibre Reinforced Polymer (FRP) Bars and FRP Reinforced Concrete Beams Under Sustained Load. Doctor’s Thesis. University of Sherbrooke, Canada.
• Zureick, A., Scott, D. (1997). Short-Term Behavior and Design of Fiber-Reinforced Polymeric Slender Members under Axial Compression. Journal of Composites Construction, 1(4): 140–149.

Publication Dates

• Publication in this collection
  2018

History

• Received
  26 Aug 2017
• Reviewed
  22 Dec 2017
• Accepted
  22 Feb 2018