Alternative active nonlinear total Lagrangian truss finite element applied to the analysis of cable nets and long span suspension bridges

Abstract

An alternative geometrically nonlinear total Lagrangian finite element is presented and applied to solve cable, cable nets and a very long suspended bridge in both three and two-dimensional spaces from its setting-up through its response to earthquake. It includes dynamics, pseudo-dynamics regularization, elastic actuators and automatic stress calibration. Dynamics and pseudo-dynamics are used to perform transient structural analysis and the setting-up of very unstable structures. Elastic actuators allow pre-stressing structural members for the iterative structural design and cables natural length definition. Automatic stress calibration comprises continuous cables in complicated structures without sliding contact devices. The formulation is applied to model main cables of suspended bridges passing through saddle points. Inertial terms are introduced by an alternative mathematical way. Two simple examples are used to validate all aspects of the proposed formulation. Finally, a representative application is performed, i.e., the numerical design and analysis of a very long span suspension bridge by the proposed strategy.

Keywords:
FEM; Suspension bridges; Cable nets; Dynamics; Earthquake analysis

Graphical Abstract

1 INTRODUCTION

Wide span structures are desired when large areas without intermediate supports are needed. One economic solution to accomplish this task is the use of cable structures. This is because the final stress levels and the self-weight of cables are lower than those presented by bending structural elements. Given the importance of this structural system, the number of scientific papers - related to the development of numerical methods and strategies to solve cable structures - is very large and, therefore, it is not our intention to make an extensive review of the subject.

This study presents a position based finite element strategy capable of a complete solution of structures composed of cables (cable nets) and truss bars, culminating here in the analysis of a very long and slender suspended bridge, from its setting-up stage to its response to earthquakes.

Cable structures are known to exhibit highly nonlinear behavior and their equilibrium always depends on its current configuration (Greco and Cuomo 2012, Kim et al. 2002, Veenendaal and Block 2012). Thus, a first challenge is determining the natural initial configuration of the structure, usually called form-finding. According to Greco and Cuomo (2012) there are several form-finding strategies in literature, but the most commonly used types are dynamic relaxation (DR), minimal surface method (MSM) and force density method (FDM). In Veenendaal and Block (2012) an interesting review of the form-finding methods with their evolution in time can be seen. In Veenendaal and Block (2012) the form-finding methods are also divided into 3 categories stiffness matrix method (Siev and Eidelman 1964), geometric stiffness method (FDM) (Schek 1974) and dynamic equilibrium method (relaxation) (Barnes 1977), where the mentioned references are considered precursor works. Also in Greco and Cuomo (2012) it is mentioned that the equivalence between FDM and MSM is proven in Wüchner and Bletzinger (2005), therefore, there is no contradiction between nomenclatures.

Regardless the different denominations of form-finding strategies, it is important to note that their basic ideas are used in analytical, numerical, or any other solution scheme. In the majority of these strategies the determination of a possible equilibrium position for cables depends on the pre-definition of a ratio between the force and length of the network parts. The choice of this ratio is not a trivial task, as the initial configuration of the cable net depends on it.

Following Crusells-Girona et al. (2017), as far as the FEM is concerned, the simplest way to model cables is by the use of straight finite elements whose geometric nonlinearity is taken into account (in most cases) by corrotational formulations (Borst et al. 2012). Most works use the above mentioned form finding formulation in order to determine of the natural position of cables as the starting point for FEM analyses.

In the present study, non-corrotational juxtaposed truss finite elements are used to model cables. We adopt a total Lagrangian finite element formulation based on positions that does not use rotation matrices (Greco et al. 2006, Greco et al. 2013, Coda 2009). One of the criticisms of this type of element is the high number of degrees of freedom used for numerical simulation. However, due to the high degree of sparsity of the truss finite element matrices, the computational cost of cable assemble is very small when compared to the remaining structural degrees of freedom. One example of this statement is the analysis of suspension bridges, as a particular interest of this study. Moreover, in this study, we adopt special sparse matrix solvers to reduce numerical efforts (The Harwell Subroutine Library Mathematical Software Library 2014).

Many works make use of curved finite elements or semi-analytical elements for form-finding or even static analysis of cables and cable networks (Greco and Cuomo 2012, Impollonia et al. 2011, Andreu et al. 2006, Such et al. 2009, Yang and Tsay 2007, Kim et al. 2016, Ahmadizadeh 2013). These elements are intended to be computationally economic, but still need improvements to present well-established performance in dynamic analysis, one of the objectives of the present study.

From the known natural position of cables, an interesting discussion in which solution models for cable and cable nets are classified in two categories - finite element methods and classic elastic catenary expressions - can be found in Abad et al. (2013). In this reference various benchmark examples (and important references) for static analysis of cable and cable nets can also be found. Regarding optimization of cable forces in stayed bridges, valuable information can be found in Zhang and Au (2014) and Thai and Kim (2012). Dynamic analysis of isolated cables can be seen in Zhou et al. (2017) and Thai and Kim (2011) among others.

As already mentioned, in this work we adopt the position-based finite element method, which is a good alternative for structural geometrically nonlinear analysis. This technique presents general and simple numerical operations for the static and nonlinear dynamic analysis of various types of structural elements and structures (Coda 2018, Soares et al. 2019, Coda 2015). Although it is not the objective of this study, it is important to mention that the simplicity of considering large deformations is also applied to large strains, see Pascon and Coda (2015,2013) for instance. In the present study, we propose an alternative active element (elastic actuator) that controls the initial length of all elements involved, enabling (i) pre-tensioning or loosening cables in the process of determining natural initial configurations and (ii) optimal stress distribution in cable structures and contiguous cables, reducing the need to model saddles or pulleys in cable connections.

An alternative dynamic relaxation technique is also provided to enable the analysis of highly nonlinear problems including setting-up problems. Two simple examples are used to validate the formulation and the motivating example of the study, a very long span suspension bridge, is presented in detail. Conclusions on the validity of the formulation and the mechanical aspects of suspension bridge analysis are shown at the end of the article.

2 TOTAL LAGRANGIAN FINITE ELEMENT

As mentioned above the truss finite element - to be used here and, therefore, adapted for the search of initial (natural) cable position and for stress balancing (calibration) in mixed structures - is based on positions rather than displacements and its formulation is briefly presented in this section.

2.1 Mechanical energy

The adopted truss element is shown in Figure 1 in a suitable form to describe the FEM based on positions. The element undergoes a change of configuration, or position, so that the associated strain energy is related only to its change in length which is a function of the initial coordinates ($i$) $X_i^{\alpha }$ and current position $Y_i^{\alpha }$ of nodes $\alpha$ (Cartesian coordinates). The truss element does not resist to transverse loading, thus we consider that only nodal forces are possible and, in this work, this property is extended to inertial forces, therefore, lumped mass is adopted. Although Fig. 1 is a two-dimensional representation, the formulation is three-dimensional and, as nodal positions are the unknown, the bar element has six degrees of freedom (three by node).

Figure 1
Change of configuration of a truss finite element

The total mechanical energy of the structural system can be expressed as

$\Pi =P+U+K$ (1)

in which the potential of conservative external loads $F_i^{\alpha }$ is given by

$P=-F_i^{\alpha }{Y}_i^{\alpha }$ (2)

and the kinetic energy by

$K=\frac{1}{2}{M}_{(\alpha )}{\dot{Y}}_i^{\alpha }{\dot{Y}}_i^{\alpha }$ (3)

The strain energy $U_e^j\left(Y_i^{\alpha }\right)$ is also a function of current nodal positions and is described as follows. Using the one dimensional Green-Lagrange strain

$\mathbb{E}=\frac{1}{2}\left(\frac{{\mathcal{l}}^2}{{\mathcal{l}}_0^2}-1\right)$ (4)

and the Saint-Venant-Kirchhoff constitutive model

$u_e^j=\frac{1}{2}K{\mathbb{E}}^2$ (5)

in which $K$ is equivalent to the Young modulus for small strains. The strain energy of element $j$ (stored in its initial configuration - Lagrangian description) is written as

$U_e^j\left(Y_i^{\alpha }\right)=A_0^{(j)}{\mathcal{l}}_0^{(j)}{u}_e^j\left(\mathbb{E}\left(Y_i^{\alpha }\right)\right)=\frac{A_0^{(j)}{\mathcal{l}}_0^{(j)}}{8}K{\left(\frac{{\mathcal{l}}^2}{{\mathcal{l}}_0^2}-1\right)}^2$ (6)

in which $A_0^{(j)}$ is the cross section area of the element, ${\mathcal{l}}_0^{(j)}$ is its initial length, ${\mathcal{l}}_{}^{(j)}$ is the current length, $u_e^j$ is the specific strain energy and $\mathbb{E}\left(Y_i^{\alpha }\right)$ is the Green-Lagrange strain.

Writing the square of the current length as

${\mathcal{l}}_{}^2={\left(Y_1^2-Y_1^1\right)}^2+{\left(Y_2^2-Y_2^1\right)}^2+{\left(Y_3^2-Y_3^1\right)}^2$ (7)

and substituting it in Eq. (5), results

$U_e^j\left(Y_i^{\alpha }\right)=\frac{A_0^{(j)}{\mathcal{l}}_0^{(j)}}{8}K{\left(\frac{{\left(Y_1^2-Y_1^1\right)}^2+{\left(Y_2^2-Y_2^1\right)}^2+{\left(Y_3^2-Y_3^1\right)}^2}{{\mathcal{l}}_0^2}-1\right)}^2$ (8)

Recalling that the total mechanical energy is the sum of the portions depending on nodes $\alpha$ and elements $j$ ,i.e., introducing Eqs. (2), (3) and (8) into Eq. (1), one writes the general expression of $\Pi \left(Y_i^{\alpha }\right)$ depending only upon nodal positions.

2.2. Motion equation

The mechanical equilibrium equation, or the motion equation, is achieved here applying the mechanical energy stationary principle over $\Pi$. In dynamics, most references prefer to do this using the Hamilton principle or the virtual work principle (Clough and Penzien 1993), however a simple and elegant procedure is used here for, in an equivalent way, arrive to the desired equation when the unknown parameters are positions. Thus the equilibrium equation is

$\delta \Pi =\delta U+\delta K+\delta P=0$ (9)

Since all members of Eq. (1) are written as a function of nodal positions, the variation described in Eq. (9) can be rewritten as

$\frac{\partial \Pi }{\partial Y_k^{\beta }}\delta Y_k^{\beta }=\left(\frac{\partial U}{\partial Y_k^{\beta }}+\frac{\partial K}{\partial Y_k^{\beta }}+\frac{\partial P}{\partial Y_k^{\beta }}\right)\delta Y_k^{\beta }=0$ (10)

and, as the nodal position variation $\delta Y_k^{\beta }$ is arbitrary, one writes the set of nonlinear equilibrium equations as

$\frac{\partial U}{\partial Y_k^{\beta }}+\frac{\partial K}{\partial Y_k^{\beta }}+\frac{\partial P}{\partial Y_k^{\beta }}=0_k^{\beta }$ (11)

In what follows the terms of Eq. (11) are developed. For the potential of conservative external forces, one writes

$\frac{\partial P}{\partial Y_k^{\beta }}=\frac{\partial \left(-F_i^{\alpha }{Y}_i^{\alpha }\right)}{\partial Y_k^{\beta }}=-F_i^{\alpha }\frac{\partial Y_i^{\alpha }}{\partial Y_k^{\beta }}=-F_i^{\alpha }{\delta }_{\alpha \beta }{\delta }_{ik}=-F_k^{\beta }$ (12)

in which $F_k^{\beta }$ is the external force applied in $k$ direction at node $\beta$.

For strain energy one uses the energy conjugate concept among internal force and position (Coda 2018, Soares et al. 2019, Coda 2015, Pascon and Coda 2015, Pascon and Coda 2013) written in a compact form as

$F_k^{\beta (\mathrm{int})}=\frac{\partial U_e}{\partial Y_k^{\beta }}=\frac{\partial {\displaystyle \sum _{j=1}^{nel}{U}_e^j}\left(Y_i^{\alpha }\right)}{\partial Y_k^{\beta }}={\displaystyle \sum _{j=1}^{nel}\frac{\partial U_e^j\left(Y_i^{\alpha }\right)}{\partial Y_k^{\beta }}}$ (13)

From Eq. (8), the last term of Eq. (13) is developed for each finite element $j$ as

${\left(F_k^{\beta (\mathrm{int})}\right)}^j=\frac{\partial U_e^j}{\partial Y_k^{\beta }}=A_0^{j}S\frac{{(-1)}^{\beta }}{{\mathcal{l}}_0^{}}\left(Y_k^2-Y_k^1\right)$ (14)

where $\beta$ assumes values 1 or 2 of the local numbering of nodes. The internal force is assembled by a cumulative process respecting the relation among local and global numbering of nodes, an usual finite element procedure.

In Eq. (14) $S$ is the second Piola-Kirchhoff stress and is given by

$S(\mathbb{E})=\frac{\partial u_e}{\partial \mathbb{E}}=K\mathbb{E}$ (15)

see Eq. (5). For the one dimensional model the relation among the second Piola-Kirchhoff stress $S$ and the Cauchy stress $\sigma$ is simply $\sigma =S\mathcal{l}/{\mathcal{l}}_0$.

For kinetic energy the differential procedure is based on the independence of the space and time variables, inherent to Newtonian mechanics. Thus the variation of kinetic energy can be written as

$\delta K=\frac{\partial K}{\partial Y_k^{\beta }}\delta Y_k^{\beta }orequivalentlyas\delta K=\frac{dK}{dt}dt$ (16)

Developing the second form of Eq. (16) and matching, one finds

$\frac{dK}{dt}dt=\frac{1}{2}{M}_{(\beta )}\left({\ddot{Y}}_k^{\beta }{\dot{Y}}_k^{\beta }+{\dot{Y}}_k^{\beta }{\ddot{Y}}_k^{\beta }\right)dt=M_{(\beta )}{\ddot{Y}}_k^{\beta }\left({\dot{Y}}_k^{\beta }dt\right)=M_{(\beta )}{\ddot{Y}}_k^{\beta }\delta Y_k^{\beta }=\frac{\partial K}{\partial Y_k^{\beta }}\delta Y_k^{\beta }$ (17)

in which the property $\delta Y_i^{\alpha }={\dot{Y}}_i^{\alpha }dt$ has been used. Thus, from Eqs. (10), (11) and (17) results the inertial force (d'Alembert principle (Clough and Penzien 1993))

$F_k^{\beta (iner)}=M_{(\beta )}{\ddot{Y}}_k^{\beta }$ (18)

Substituting Eqs. (12), (13) and (18) into Eq. (11) yields the set of nonlinear dynamic equilibrium (or motion) equations in compact form

$F_k^{\beta (\mathrm{int})}+F_k^{\beta (iner)}-F_k^{\beta }=0_k^{\beta }$ (19)

or passing the notation from node $\beta$ and direction $k$ to degree of freedom $i=d\left(\beta -1\right)+k$ one writes Eq. (19) as

$F_i^{\mathrm{int}}+M_i{\ddot{Y}}_i-F_i^{ext}=0_i^{}or{\overrightarrow{F}}_{}^{\mathrm{int}}+M\cdot \overrightarrow{\ddot{Y}}-{\overrightarrow{F}}_{}^{ext}=\overrightarrow{0}$ (20)

in which $M$ is a diagonal mass matrix stored in a single vector and $d$ stands for dimension (3 or 2). It is interesting to recall that the internal force, Eq. (14) is calculated at the element level and globally assembled.

2.3. Solution process - transient problem

The equation of motion (20) is rewritten using compact notation as

$\overrightarrow{g}=\frac{\partial \Pi }{\partial \overrightarrow{Y}}=\frac{\partial U}{\partial \overrightarrow{Y}}(\overrightarrow{Y})+M\cdot \ddot{\overrightarrow{Y}}+C\cdot \dot{\overrightarrow{Y}}-{\overrightarrow{F}}_{}^{ext}(t)=\overrightarrow{0}$ (21)

where, for convenience, one recovers the meaning of the internal force as the derivative of the strain energy regarding nodal positions, Eq. (14). The damping matrix is included as simply mass proportional, as it is not the objective of this study to have further discussions on this subject. Equation (21) is valid at any instant of the structural analysis and time is a continuous variable. However, the numerical solution requires time to be treated discretely, i.e., the current instant is calculated as the previous instant plus the time interval, such as

$t_{s+1}=t_s+\Delta t$ (22)

in which $t_{s+1}$ is the current instant. Equation (21) is valid at the current instant and is rewritten as

$\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)={\frac{\partial \Pi }{\partial \overrightarrow{Y}}|}_{s+1}={\frac{\partial U}{\partial \overrightarrow{Y}}|}_{s+1}+M\cdot {\ddot{\overrightarrow{Y}}}_{s+1}+C\cdot {\dot{\overrightarrow{Y}}}_{s+1}-{\overrightarrow{F}}_{s+1}^{ext}=\overrightarrow{0}$ (23)

where, by simplicity, $t$ is omitted. The Newmark´s time integrator is used to solve Eq. (23) over time in the following form Warburton (1976)

${\overrightarrow{Y}}_{s+1}={\overrightarrow{Y}}_s+{\overrightarrow{\dot{Y}}}_s\Delta t+\left[\left(\frac{1}{2}-\beta \right){\overrightarrow{\ddot{Y}}}_s+\beta {\overrightarrow{\ddot{Y}}}_{s+1}\right]\Delta t^2$ (24)

${\overrightarrow{\dot{Y}}}_{s+1}={\overrightarrow{\dot{Y}}}_s+\left(1-\gamma \right)\Delta t{\overrightarrow{\ddot{Y}}}_s+\gamma \Delta t{\overrightarrow{\ddot{Y}}}_{s+1}$ (25)

in which $\beta$ and $\gamma$ are free parameters of the method, adopted here as $\beta =1/4$ and $\gamma =1/2$ (Warburton 1976).

It is of interest to write current velocity and acceleration as a function of current positions and known values of the past. In order to do that, it is enough to isolate these variables in Eqs. (24) and (25), as

${\overrightarrow{\ddot{Y}}}_{s+1}=\frac{{\overrightarrow{Y}}_{s+1}}{\beta \Delta t^2}-\left(\frac{{\overrightarrow{Y}}_s}{\beta \Delta t^2}+\frac{{\dot{\overrightarrow{Y}}}_s}{\beta \Delta t^2}+\left(\frac{1}{2\beta }-1\right){\ddot{Y}}_s\right)=\frac{{\overrightarrow{Y}}_{s+1}}{\beta \Delta t^2}-{\overrightarrow{Q}}_s$ (26)

${\dot{\overrightarrow{Y}}}_{s+1}=\frac{\gamma }{\beta \Delta t}{\overrightarrow{Y}}_{s+1}+\left[{\dot{\overrightarrow{Y}}}_s+\Delta t\left(1-\gamma \right){\ddot{\overrightarrow{Y}}}_s\right]-\left(\frac{{\overrightarrow{Y}}_s}{\beta \Delta t^2}+\frac{{\dot{\overrightarrow{Y}}}_s}{\beta \Delta t^2}+\left(\frac{1}{2\beta }-1\right){\ddot{Y}}_s\right)\gamma \Delta t=\frac{\gamma }{\beta \Delta t}{\overrightarrow{Y}}_{s+1}+{\overrightarrow{R}}_s-\gamma \Delta t{\overrightarrow{Q}}_s$ (27)

with

${\overrightarrow{Q}}_s=\left(\frac{{\overrightarrow{Y}}_s}{\beta \Delta t^2}+\frac{{\dot{\overrightarrow{Y}}}_s}{\beta \Delta t^{}}+\left(\frac{1}{2\beta }-1\right){\ddot{Y}}_s\right)and{\overrightarrow{R}}_s=\left[{\dot{\overrightarrow{Y}}}_s+\Delta t\left(1-\gamma \right){\ddot{\overrightarrow{Y}}}_s\right]$ (28)

Substituting Eqs. (26) and (27) into Eq. (23) results

$\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)={\frac{\partial U}{\partial \overrightarrow{Y}}|}_{s+1}+\frac{M}{\beta \Delta t^2}{\overrightarrow{Y}}_{s+1}-M\cdot {\overrightarrow{Q}}_s+\frac{\gamma C}{\beta \Delta t}\cdot {\overrightarrow{Y}}_{s+1}+C\cdot {\overrightarrow{R}}_s-\gamma \Delta tC\cdot {\overrightarrow{Q}}_s-{\overrightarrow{F}}_{s+1}^{ext}(t)=\overrightarrow{0}$ (29)

Equation (29) is reduced to $\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)=0$, resulting in a nonlinear set of equations with the current position vector $\left({\overrightarrow{Y}}_{s+1}\right)$ being the unknown. It will be solved by the Newton-Raphson procedure.

A first-order truncated Taylor series expansion is performed as follows

$\overrightarrow{0}=\overrightarrow{g}({\overrightarrow{Y}}_{s+1}^{})\cong \overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)+\nabla \overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)\Delta \overrightarrow{Y}$ (30)

in which

$\nabla \overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)=H={\frac{{\partial }^2\Pi }{\partial {\overrightarrow{Y}}^2}|}_{s+1}={\frac{{\partial }^{2}U}{\partial Y^2}|}_{s+1}+\frac{M}{\beta \Delta t^2}+\frac{\gamma C}{\beta \Delta t}=H^{estat}+\frac{M}{\beta \Delta t^2}+\frac{\gamma C}{\beta \Delta t}$ (31)

being $H^{estat}$the static Hessian matrix given for a finite element and written for local nodes and global directions, see Eq. (8), as

${\left(H_{ik}^{\alpha \beta }\right)}^j={(-1)}^{\beta }{(-1)}^{\alpha }\frac{A_0^{(j)}}{{\mathcal{l}}_0^{}}\left(K\frac{\left(Y_i^2-Y_i^1\right)}{{\mathcal{l}}_0}\frac{\left(Y_k^2-Y_k^1\right)}{{\mathcal{l}}_0}+S{\delta }_{ik}\right)$ (32)

where $\alpha$ and $\beta$ are local nodes numbers (1 or 2) and ${\delta }_{ik}$ is the Kronecker delta (or identity matrix). The global Hessian matrix is assembled in the same way of internal forces, but respecting advanced strategies for using sparse matrices solvers (The Harwell Subroutine Library Mathematical Software Library 2014).

From Eq. (31) one writes the linear system used to find a position correction at the analyzed instant, i.e.

$\nabla \overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)\cdot \Delta \overrightarrow{Y}=-\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)orH\cdot \Delta \overrightarrow{Y}=-\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)$ (33)

where ${\overrightarrow{Y}}_{s+1}^0$ is a trial position. At the beginning of a time step the trial position is the solution of the previous step, i.e., ${\overrightarrow{Y}}_s^{}$. Solving the position correction $\Delta \overrightarrow{Y}$ by Eq. (33) a new trial solution for ${\overrightarrow{Y}}_{s+1}$ is achieved by

${\overrightarrow{Y}}_{s+1}^0={\overrightarrow{Y}}_{s+1}^0+\Delta \overrightarrow{Y}$ (34)

The velocity and acceleration should be recalculated, at each iteration, using Eqs. (26) and (27) rewritten in a compact form as

${\overrightarrow{\ddot{Y}}}_{s+1}=\frac{{\overrightarrow{Y}}_{s+1}}{\beta \Delta t^2}-{\overrightarrow{Q}}_{S}and{\dot{\overrightarrow{Y}}}_{s+1}=\frac{\gamma }{\beta \Delta t}{\overrightarrow{Y}}_{s+1}+{\overrightarrow{R}}_s-\gamma \Delta t{\overrightarrow{Q}}_s$ (35)

It is interesting to remember that ${\overrightarrow{Q}}_s$ and $R_s$ remain constant during iterations as they are values of the past. They are updated at the end of a time step. The stop criterion of the iterative procedure is given by

$\frac{\Vert \overrightarrow{g}(Y^0)\Vert }{\Vert {\overrightarrow{F}}^{ext}\Vert }\le TOLor\frac{\Vert \Delta \overrightarrow{Y}\Vert }{\Vert \overrightarrow{X}\Vert }\le TOL$ (36)

in which $TOL$ is the tolerance defined by the designer.

When results converge ${\overrightarrow{Y}}_{s+1}^0$ becomes the proper solution ${\overrightarrow{Y}}_{s+1}$ that, for the next step, will again be the first trial solution

At the first time step the acceleration should be calculated by the equation of motion as

${\overrightarrow{\ddot{Y}}}_0=M^{-1}\left[{{\overrightarrow{F}}_0^{ext}-\frac{\partial U}{\partial Y}|}_0-C{\dot{Y}}_0\right]$ (37)

It is worth noting that the static process is the same, only all values containing mass and damping are neglected.

3. ACTIVE ELEMENT (ACTUADTOR) AND ALTERNATIVE DUNAMIC RELAXATION PROCESS

The search for the initial (natural) configuration of cables and cable networks will be made by the same element used to solve the mechanical problem. In order to do so a small modification of the presented formulation is necessary. The numerical procedure is inspired in mechanical actuators whose change in nominal (initial) length is made by controlling the volume of fluid inserted into their hydraulic chamber, regardless of whether or not external loading exists. In our finite element, we perform the control of the nominal length by imposing a length variation along the searching of cables natural initial position analysis.

3.1. Imposed length variation

In order to impose the length variation one applies (divided by load steps or time steps) an increment $\Delta \mathcal{l}$ upon the original numeric length ${\mathcal{l}}_0$ of the finite element, creating the nominal length of the element, called here the natural initial length of the cable element ${\mathcal{l}}_{0n}$, as

${\mathcal{l}}_{0n}={\mathcal{l}}_0+\Delta \mathcal{l}(t)$ (38)

in which $t$ represents time or load step and $\Delta \mathcal{l}(0)=0$ and ${\mathcal{l}}_{0n}(0)={\mathcal{l}}_0$. The increment is applied until ${\mathcal{l}}_{0n}$ reaches its final desired value. The natural initial length substitutes the original initial length in all presented formulae. For example, the Green strain, Eq. (4) becomes

$\mathbb{E}=\frac{1}{2}\left(\frac{{\mathcal{l}}^2}{{\mathcal{l}}_{0n}^2}-1\right)$ (39)

and the strain energy stored in a finite element, Eq. (6) becomes

$U_e^j\left(Y_i^{\alpha }\right)=A_0^{(j)}{\mathcal{l}}_{0n}^{(j)}{u}_e^j\left(\mathbb{E}\left(Y_i^{\alpha }\right)\right)=\frac{A_0^{(j)}{\mathcal{l}}_{0n}^{(j)}}{8}K{\left(\frac{{\mathcal{l}}^2}{{\mathcal{l}}_{0n}^2}-1\right)}^2$ (40)

that is, the referring space becomes the natural initial space (current undeformed length) of the structural element (cable). In practice, changing the length of the element changes the value of the internal force influencing the equilibrium equation, which changes the equilibrium configuration of the structure (cable) until achieving the natural position of the cable.

3.2. Stress calibration

For stress calibration on a cable element, we also use the concept of natural initial length, but the length variation to be imposed is calculated using the difference between the calculated stress $S_c$ and the design stress $S_d$ as

$\Delta \mathcal{l}=\frac{(S_c-S_d)}{\left|S_c-S_d\right|}\left(\sqrt{{\mathcal{l}}^2+(2{\mathcal{l}}_{0n}^2\left|(S_c-S_d)/K\right|)}-\mathcal{l}\right)$ (41)

The calculated $\Delta \mathcal{l}$ is applied in Eq. (38) throughout the Newton-Raphson iterative process until an acceptable $\Delta \mathcal{l}$ is achieved. In this case an additional stop criterion is introduced

$\frac{\left|S_c-S_d\right|}{\left|S_c+S_d\right|}<NTOL$ (42)

for which a new tolerance ($NTOL$) is also defined by the user.

3.3. Stress balance

The difference between the stress balance technique and the stress calibration is that in the latter the stress level intended for a particular cable is predefined by the user, whereas in the former the final stress level is unknown and should be calculated in order to cables have the same final stress.

Thus, for two finite elements, during each iteration of the Newton-Raphson method, the trial balanced stress ($S_d$) stress is defined as a weighted average of the calculated stresses as

$S_d=\left(\frac{S_{c1}}{A_1{K}_1}+\frac{S_{c2}}{A_2{K}_2}\right)\frac{A_1{K}_1{A}_2{K}_2}{A_1{K}_1+A_2{K}_2}$ (43)

in which indices 1 and 2 represent local numbering for the two involved finite elements. From the knowledge of the trial balanced stress level one applies Eq. (41) for each element calculating $\Delta {\mathcal{l}}_1$ and $\Delta {\mathcal{l}}_2$. These length changes are introduced in Eq. (38) throughout the Newton-Raphson iterative process until Eqs. (36) and (42) are satisfied.

When the cable has high stiffness the required length change to achieve the desired stress is small and vice versa. Thus, the weighting proposed in Eq. (43) is such that the stress level is closer to the one presented by the cable element of lower stiffness.

3.4. Alternative dynamic relaxation

As it is well known, see Eq. (32), when cables stress levels are near zero or negative the static Hessian matrix is near singular or even singular. The alternative dynamic relaxation used here (called pseudo-dynamic) is based on usual dynamic techniques (Barnes 1977), in which the presence of mass and damping matrices in the total Hessian matrix (Eq. (31)) eliminates singularity.

In classical dynamic relaxation it is necessary to calibrate the damping level so that it approaches the critical damping leading to a time solution that asymptotically leads the dynamic analysis of cables to the static natural initial position of the cable (Barnes 1977, Lewis et al. 1984).

In the proposed alternative pseudo-dynamic relaxation technique the dynamic procedure is adapted by removing the inertial and damping forces, maintaining only the mass matrix and, if desired, the damping matrix, see Eq. (31). This procedure eliminates negative stresses (compression) that can be generated by inertial forces in the static portion of total Hessian. This procedure also eliminates vibration.

In the proposed pseudo-dynamic relaxation, the total Hessian matrix results

$H=H^{estat}+\frac{1}{ite}\frac{M}{\Delta t^2}$ (44)

in which $ite$ is the iteration number of the Newton-Raphson procedure that gradually reduces the mass influence and $\Delta t$ is a pseudo-dynamic time step. The experience of the user (structural designer) in dynamic analyses facilitates the choice of this pseudo $\Delta t$.

The solution process is the same described in section 2.3, however for each iteration the matrix dynamic contribution is reduced and, after convergence, an additional iteration is performed without the mass matrix, i.e., with $H=H^{estat}$, that is, finding the static result for each pseudo time step.

3.5 - Basic algorithms:

The algorithms presented in this section are separated by their use. The first summarizes the solution of the dynamic problem after finding the initial position of the structure. The second describes the static solution including the pseudo dynamic regularization to find the natural position. The last algorithm introduces both stress calibration and stress balance over the static solution. The first and the second use $\Delta \mathcal{l}$ as input value and the last calculates it to find the design stress.

3.5.1 - Dynamic Solution

It is important to recall that a time step is called $s$.

• a) The first trial solutions for position, velocity and acceleration, are assumed as the solution of the previous time step, i.e., ${\overrightarrow{Y}}_{s+1}^{}={\overrightarrow{Y}}_{s+1}^0={\overrightarrow{Y}}_s$, ${\overrightarrow{\dot{Y}}}_{s+1}^{}={\overrightarrow{\dot{Y}}}_{s+1}^0={\overrightarrow{\dot{Y}}}_s$, ${\overrightarrow{\ddot{Y}}}_{s+1}^{}={\overrightarrow{\ddot{Y}}}_{s+1}^0={\overrightarrow{\ddot{Y}}}_s$. The updated value and the trial value are stored in the same vector as, at the end of iterations, they coincide with the step solution.

• b) Take the applied external force ${\overrightarrow{F}}_{s+1}^{ext}(t_{s+1})$, the prescribed position ${\overrightarrow{Y}}_{s+1}^{}(t_{s+1})$ and/or the prescribed $\Delta \mathcal{l}$ and calculate ${\overrightarrow{Q}}_s$ and ${\overrightarrow{R}}_s$ using expressions (28).

• c) Using the trial position and the imposed $\Delta \mathcal{l}$, compute the internal force $F_{s+1}^{\mathrm{int}}$ by Eq. (14) and the static Hessian matrix $H_{s+1}^{estat}$ by Eq. (32).

• d) Assemble the global internal force and static Hessian matrix, complete the Hessian matrix with dynamic parts using Eq. (31).

• e) Assemble the unbalanced vector $\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^{}\right)$ by Eq. (29), observe that at this point the inertial force is introduced.

• f) Solve $\Delta \overrightarrow{Y}$ by Eq.(33)

• g) Update the current position by Eq. (34)

• h) Update acceleration and velocity by Eq. (35)

• i) Evaluate $\left|\Delta \overrightarrow{Y}\right|/\left|\overrightarrow{X}\right|$.

  • i1) If $\left|\Delta \overrightarrow{Y}\right|/\left|\overrightarrow{X}\right|<tol$ the equilibrium position is found at this time step. Store the solution of this step as past values for the next step (first trial), i.e., ${\overrightarrow{Y}}_S={\overrightarrow{Y}}_{S+1}$, ${\overrightarrow{\dot{Y}}}_S={\overrightarrow{\dot{Y}}}_{S+1}$ and ${\overrightarrow{\ddot{Y}}}_S={\overrightarrow{\ddot{Y}}}_{S+1}$, go back to item (b) for the next time step.

  • i2)if$\left|\Delta \overrightarrow{Y}\right|/\left|\overrightarrow{X}\right|\ge tol$, go to item (c) and start a new iteration in order to improve the solution precision.

3.5.2 - Static solution with pseudo dynamics

In the static solution $s$ means a load or position step.

• a) The first trial solution for position is the solution of the previous step, i.e., ${\overrightarrow{Y}}_{s+1}^{}={\overrightarrow{Y}}_{s+1}^0={\overrightarrow{Y}}_s$. The updated value and the trial value are stored in the same vector as, at the end of iterations, they coincide with the step solution. Assume the first iteration as $ite=1$.

• b) Take the applied external force ${\overrightarrow{F}}_{s+1}^{ext}(t_{s+1})$, the prescribed position ${\overrightarrow{Y}}_{s+1}^{}(t_{s+1})$ and/or the prescribed $\Delta \mathcal{l}$ all input values.

• c) Using the trial position and the imposed $\Delta \mathcal{l}$, compute the internal force $F_{s+1}^{\mathrm{int}}$ by Eq. (14) and the static Hessian matrix $H_{s+1}^{estat}$ by Eq. (32).

• d) Assemble the global internal force and static Hessian matrix

• e) Complete the Hessian matrix with the pseudo dynamic part using Eq. (44).

• f) Assemble the unbalanced vector using this equation $\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)={\frac{\partial U}{\partial \overrightarrow{Y}}|}_{s+1}-{\overrightarrow{F}}_{s+1}^{ext}(t)$

• g) Solve $\Delta \overrightarrow{Y}$ at Eq.(33)

• h) Update the current position by Eq. (34)

• i) Evaluate $\left|\Delta \overrightarrow{Y}\right|/\left|\overrightarrow{X}\right|$.

  • i1) If $\left|\Delta \overrightarrow{Y}\right|/\left|\overrightarrow{X}\right|<tol$ go to item (c) and start the last iteration without executing item (e). The equilibrium position is found, store the solution of this step as previous values to the next step (first trial), i.e., ${\overrightarrow{Y}}_S={\overrightarrow{Y}}_{S+1}$ go back to item (b) for the next step.

  • i2) If $\left|\Delta \overrightarrow{Y}\right|/\left|\overrightarrow{X}\right|\ge tol$, update $ite=ite+1$and go to item (c) and start a new iteration in order to improve the solution precision.

3.5.3 - Static solution with stress calibration and stress balance

This analysis occurs after the first natural static position is found by algorithm 3.5.2. It means that the input data are constant and satisfies equilibrium with the internal force. Additional $\Delta \mathcal{l}$ are the new unknown for stress adjustments.

• a) The first trial solution for position is the solution of the previous step, i.e., ${\overrightarrow{Y}}_{s+1}^{}={\overrightarrow{Y}}_{s+1}^0={\overrightarrow{Y}}_s$. The updated value and the trial value are stored in the same vector as, at the end of iterations, they coincide with the step solution.

• b) Take the applied external force ${\overrightarrow{F}}_{s+1}^{ext}(t_{s+1})$, the prescribed position ${\overrightarrow{Y}}_{s+1}^{}(t_{s+1})$ and/or the prescribed $\Delta \mathcal{l}$ (that compose ${\mathcal{l}}_{0n}$) all input values that are now constant.

• c) Using the trial position and prescribed values, calculate stress $S_c$ at element level by Eq. (15) considering (39).

• c1) If there is stress calibration calculate the trial $\Delta \mathcal{l}$ by Eq. (41) and update ${\mathcal{l}}_{0n}$ using (38)

• c2) If there is stress balance calculate $S_d$ by Eq. (43), the trial $\Delta \mathcal{l}$ by Eq. (41) and update ${\mathcal{l}}_{0n}$ using (38)

• d) Compute the internal force $F_{s+1}^{\mathrm{int}}$ by Eq. (14) and the static Hessian matrix $H_{s+1}^{estat}$ by Eq. (32).

• e) Assemble the global internal force and static Hessian matrix

• f) Assemble the unbalanced vector by equation $\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)={\frac{\partial U}{\partial \overrightarrow{Y}}|}_{s+1}-{\overrightarrow{F}}_{s+1}^{ext}(t)$

• g) Solve $\Delta \overrightarrow{Y}$ at Eq. (33)

• h) Update the current position by Eq. (34)

• i) Evaluate $\left|\Delta \overrightarrow{Y}\right|/\left|\overrightarrow{X}\right|$.

  • i1) If $\left|\Delta \overrightarrow{Y}\right|/\left|\overrightarrow{X}\right|<tol$ and $\left|S_c-S_d\right|/\left|S_c+S_d\right|<NTOL$ the equilibrium position and the designed stress are achieved.

  • i2) If $\left|\Delta Y_k\right|/\left|X_k\right|\ge tol$ or $\left|S_c-S_d\right|/\left|S_c+S_d\right|\ge NTOL$, go to item (c) and start a new iteration in order to improve the solution precision.

4. EXAMPLES

The example section is divided in three cases. The first shows that the proposed active element and the pseudo-dynamic relaxation technique used together make possible to find the natural static position of cables and also give very good results in force. The second example shows that the proposed formulation is capable to solve loose cable nets, finding the natural static position and solving (in position and force) complicate loading situations. Finally, the third example is a particular interest of the proposed study by the importance of suspension bridges in the current engineering state of art. It is a study of a very long suspended bridge from its setting-up through an earthquake analysis. Readers interested in more benchmarks to check static formulations are invited to see, for example, Abad et al. (2013).

4.1. Simple cable modelling

This is a very simple example and is used to check the active element as a tool to find cable natural position (form-finding) and the positional truss finite element performance. The presented formulation is used to model a steel cable with the following properties: elastic modulus $K=E=210GPa$, density $\rho =7000kg/m^3$, cross section radius $r={10}^{-2}m$. Two values are chosen for the cable length, $L_1=10.5m$ and $L_2=11m$. These cables should overpass a horizontal distance of $\mathcal{l}=10m$ with a gravity acceleration of $g=10m/s^2$, see Fig. 2. The numerical simulation is divided in three phases.

Figure 2
Horizontal cable subjected to self-weight

In the first phase the pseudo-dynamic regularization is applied making $H_{tot}=H_{est}+M/\Delta t^2/ite$ in which $ite$ is the number of numerical iteration of the Newton-Raphson procedure and $\Delta t=1.0s$ is the fictitious time step. For this case only one load step is sufficient to find the final configuration. The initial configuration of the cable is straight with a length of $10m$. One observes that after the convergence ($tol=1.0x{10}^{-9}$ in positions) an additional iteration is performed without regularization. The transverse displacement of the center of the cable (sag) is, in this phase, $f=5.39cm$.

In the second phase all cable elements are considered elastic actuators (actives) and their length are changed in a smooth way in order to avoid compression. The total length change is $\Delta L=0.5m$ (for $L_1=\mathrm{10,5}m$), that is, $\Delta L_j=0.025m$ for each element (20 finite elements) divided in 100 pseudo time steps $\Delta t=0.01s$ used to regularize the Hessian matrix. For $L_2=11m$ we adopt $\Delta L_j=0.05m$ for each element (20 finite elements) divided in 200 pseudo time steps of $\Delta t=0.005s$.

The second phase result corresponds to the final cable position including the self-weight load and the flexibility of the cable. In the third phase we withdraw the load (one step) without regularization, in order to find the inextensible (natural) position of the cable.

For $L_1=10.5m$, the final displacement with the load is $f_{load}=1.39481m$ and without the load is $f_{unload}=1.39476m$, i.e., the difference is less than $5x{10}^{-5}m$ and, thus, the flexibility of this cable is negligible for the self-weight.

In Table 1 we present comparisons between numerical and analytical solutions (displacements and forces), in which $H$ represents horizontal forces and $V$ vertical forces at supports. Analytical solutions are well known and written in hyperbolic cosines and sines. One observes that the differences between solutions are very small.

Table 1
Values comparison for 20 finite elements

$L_i$	$f_{ana}$	$f_{unload}$	$f_{load}$	$H_{ana}$	$H_{load}$	$V_{ana}$	$V_{load}$
10.5m	1.393m	1.395m	1.395m	202.228N	202.020N	115.448N	115.453N
11.0m	2.003m	2.005m	2.006m	144.027N	143.893N	120.945N	120.945N

Only to be complete, keeping the same physical properties, we modeled a cable with length $L=11.0m$ that should overpass the slopping distance between points (0,1)m and (10,2)m, see Fig. 3. The distance among supports is $d_0=10.0499m$ and is used to define the initial configuration of the cable, see Fig. 3.

Figure 3
Cable overpassing a slopping distance

We use the same phases of the horizontal cable and$L_2=11.0m$ to run this case. For this length, we adopted $\Delta L_j=0.04750622m$ for 20 elements and 200 pseudo time steps of $\Delta t=0.005s$. As a result in force we have $H_{ana}=147.328N$, $H_{num}=147.187N$, $V_{ana}^1=103.573N$, $V_{num}^1=103.589N$, $V_{ana}^2=138.317N$, $V_{num}^2=138.301N$.

Both solutions are in equilibrium with the applied load and the maximum relative difference is less than 0.1%. Thus, we do not see any reason to use more than 20 finite elements for each cable to solve usual engineering problems. It is important to mention that the calculated vertical reaction is the automatic internal force plus the external force corresponding to half element connected to the support.

The number of variables stored to assemble the sparse Hessian matrix of 20 finite elements is 396 for a 3D analysis and 176 for a 2D analysis (The Harwell Subroutine Library Mathematical Software Library 2014). This number of variables is irrelevant for the current computational resources. The formulas $ST=18\cdot \left(nn-2\right)+54$ for 3D and $ST=8\cdot \left(nn-2\right)+24$ for 2D analyses can be used to calculate the number of stored variables ($ST$) in the system of equations to solve a single cable with ($nn$) nodes.

4.2. Cable net

This example was proposed by Greco et. al. (2014) and is a 5-cable network as shown in Fig. 4a. It is an important example for the validation of cable structure solution codes. The authors used C-FDM to find an initial configuration and then solved load cases, one of which was chosen for comparison.

As our formulation does not use a classical form-finding process, we will analyze the same problem using the active elements to find the natural initial position. The initial numerical problem is plane and has the following node initial coordinates $P1=(0.5;0.25;0)$, $P2=(0.5;0.75;0)$, $P3=(0;0;0)$, $P4=(0;1;0)$, $P5=(1;0;0)$ and $P6=(1;1;0)$, see Fig. 4a. We use 4 different meshes, for all of them the division of cables 1, 2, 4 and 5 has two times more elements than cable 3. The number of elements adopted for cable 3 that defines the discretization is 5, 10, 20 and 40, resulting in discretization with 46, 91, 181 and 371 nodes, respectively.

Figure 4
(a) Initial configuration, (b) Final configuration for 91 nodes

The cables properties are $K=E=175MPa$, $\rho =7000kg/m^3$ and $A=2.857x{10}^{-4}{m}^2$. The adopted gravity acceleration is $g=10m/s^2$, resulting in the same load used by reference, i.e., $q=20N/m$ accompanying the current position.

The search of the natural self-weight configuration is divided in three phases. The first uses the pseudo-dynamic regularization with fictitious time step $\Delta t=1.0s$ with a total of $10$ steps. Fig. 5 shows the equilibrium configuration and cable force distribution at the end of phase 1.

Figure 5
Displacements and cable normal force for phase 1 - 91 nodes

In phase 2 we applied 100 static steps without Hessian regularization pulling support 6 until position $P6=(1;1;1)$. In Fig. 6 final position and normal forces for all cables at the end of phase 2 are shown.

Figure 6
Displacements and cable normal force for phase 2 - 91 nodes

In phase 3, we apply the pseudo-dynamic Hessian regularization with $\Delta t=0.1s$ and apply a change of length for all finite elements in 100 steps. The initial lengths of the finite elements are calculated from the initial configuration (Fig. 4) and the final lengths are given by the length of cables (Greco et al. 2014).To facilitate the reproduction of this example, the initial and final lengths are given in Table 2. Figure 7 presents the final displacement results and normal forces in cables.

Table 2
Initial and final lengths of cables

$L(m)/Cable$	1	2	3	4	5
$l_0$	0.5590	0.5590	0.5000	0.5509	0.5590
$l_{n0}=l_f$	1.2887	1.2887	0.5912	1.1874	2.0978

Figure 7
Displacements and forces in cables - 381 nodes

Table 3 compares horizontal forces at points 3, 4, 5 and 6 (see Fig. 4) with values given by Greco et al. (2014). Table 4 compares the vertical forces.

Table 3
Horizontal forces (N)

Discr./point	Node 3	Node 4	Node 5	Node 6
46	5.87573	5.87964	5.87715	5.87748
91	5.86845	5.87357	5.87297	5.86665
181	5.86672	5.87203	5.87192	5.86404
371	5.86630	5.87168	5.87168	5.86342
Reference	5.87000	5.87000	5.87000	5.87000

Table 4
Vertical forces (N)

Discr./point	Node 3	Node 4	Node 5	Node 6
46	27.93301	25.33133	27.93467	47.87679
91	27.92903	25.33023	24.93427	47.88235
181	27.92804	25.32961	24.93406	47.88422
371	27.92777	25.32947	24.93399	47.88472
Reference	27.92800	25.32800	24.93400	47.88700

Using the values of Tables 3 and 4, it can be seen in Fig. 8 that, for any discretization, forces are practically the same as those presented by the reference for any discretization used, i.e., the convergence of results is very fast.

Figure 8
Relative difference among forces at extremities $rd=100*\left|F_{discr}-F_{ref}\right|/F_{ref}$

Table 5 compares the final positions found for points 1 and 2 of Fig. 4 (elastic case) with those given by Greco et al. (2014). As can be seen, there is no noticeable difference.

Table 5
Comparison of free node positions for various discretizations (m)

Coord./Node	X	Y	Z
$P_1^{ref}$	0.499300	0.249900	-1.11480
$P_1^{46}$	0.499984	0.249294	-1.11655
$P_1^{91}$	0.499947	0.249757	-1.11524
$P_1^{181}$	0.499938	0.249872	-1.11491
$P_1^{371}$	0.499936	0.249905	-1.11482
$P_2^{ref}$	0.499400	0.750000	-0.99630
$P_2^{46}$	0.499856	0.750818	-0.99789
$P_2^{91}$	0.499540	0.750211	-0.99673
$P_2^{181}$	0.499470	0.750060	-0.99643
$P_2^{371}$	0.499452	0.750019	-0.99636

The analysis is continued by applying an increasing horizontal force at point 2 until the value $\overrightarrow{F}=\left(\mathrm{0,}-\mathrm{100,0}\right)N$. We adopted the discretization called 371 and 100 for pseudo time steps $\Delta t=0.1s$. Figure 9 presents graphics that relate the forces at points 3, 4, 5 and 6 to the intensity of the applied force at point 2. It is important to mention that these results are in very good agreement with ones presented by Greco et al. (2014).

Figure 9
Horizontal and vertical forces at supports

In order to show that coarser discretizations also give excellent results, we show in Fig. 10 the relative difference among forces at the extremity of cable 4 ($dr=100\cdot \left|F_{371}-F_i\right|/F_{371}$) for discretization $i$ (number of nodes) regarding discretization 371 for the same increasing loading. As one can see the largest difference is of $0.35\%$ in the vertical reaction for discretization 46.

Figure 10
Relative difference at extremity of cable 4

Finally, Fig. 11 shows the vertical displacements of nodes 1 and 2 for the same increasing force applied at node 2 and Fig. 12 shows some snapshots for the load levels. $10N$, $30N$, $60N$ and $100N$.

Figure 11
Displacements for nodes 1 and 2 for horizontal increasing applied force.

Figure 12
Snapshots for different load levels

These first two examples show how the solution of cables and cable networks can be accurately made using positional geometric nonlinear truss elements associated with pseudo dynamic regularization and active element strategy. Moreover, it is observed that the solution convergence is very fast and that only 10 truss elements per cable length is sufficient to have very satisfactory answers to static problems. Next example is important to show the capacity of the proposed formulation, i.e., to solve a very flexible suspended structures (in this case a bridge) from this setting-up through an earthquake response.

It is important to note that the number of variables stored to constitute the sparse Hessian matrix of 10 finite elements is 216 for a 3D analysis and 96 for a 2D analysis (The Harwell Subroutine Library Mathematical Software Library 2014). This number of variables is irrelevant for the current computational resources.

4.3. Very long suspended bridge

After verifying that the proposed active truss element (associated with the pseudo-dynamic regularization technique adopted here) are able to solve cable and cable network problems, we solve a 4020m long suspension bridge with 2000m center span and 1000m of secondary spans, see Figs. 13, 14 and 15. Due to scale, some dimensions are not shown in figures, the pillars width is 10 m, the height of the girder is 15m and the distance of the lower points of main cables to the girder surface is 15m (at extremities of the bridge). At the center of the bridge the initial distance of the lower point of the main cable to the girder surface is 15m, but this distance will change in the process of the bridge setting-up. Moreover, pillars are square constituted by truss elements and the girder has 40m of width. This example is inspired in the Akashi-Kaikyo bridge Japan. Although not being the same problem, some mechanical information are adapted (not equal) to the ones given in Furuya et al. (1994), Yim (2007) and Miyata and Yamaguchi (1993) and other cited therein.

Figure 13
Basic dimensions of the modeled bridge - 2D model mesh.

Figure 14
General vision of the adopted 3D mesh

Figure 15
Detail of girder, pillar, main and suspension cables - 3D mesh and nodes

The structural elements have the following cross sections: main cables $A_{mc}=1.131m^2$, suspension cables $A_{sc}=0.0341m^2$, pillar truss elements $A_{pb}=0.784m^2$, main girder bars $A_{mg}=0.085m^2$ and secondary girder bars $A_{sg}=0.0304m^2$. The adopted material has elastic modulus of $K=E=200GPa$ and mass density $\rho =7000kg/m^3$. In order to calculate the self-weight we used $g=10m/s^2$ as gravity acceleration. We increased in 2.5% the girder (3D model) bars self-weight as additional loading (traffic and equipments).

For both discretizations (2D or 3D) we used 50 truss elements to model each main cable of secondary spans and 100 elements to model the central span main cables.

The solution of the problem involves (i) positioning the main cables using a 2D model and performing stress balance, (ii) the static solution of the 2D and 3D models without adjusting the suspension cables, (iii) adjusting the suspension cables to the proper positioning of the girder in its ideal static configuration (2D and 3D), (iv) shortening of loose cables (2D and 3D), (v) the transient dynamic analysis of the bridge under 3D real seismic shake and, (vi) determining the frequencies and modes of vibration of the positioned 3D structure susceptible to the shake.

i) Positioning the main cables

In order to find the natural initial position of main cables we begin with a two-dimensional analysis in which these cables are modeled with straight initial configuration, as shown in Fig. 16.

Figure 16
Initial straight and deformed loaded cable with final stress level (GPa)

In the first phase we use the pseudo dynamic regularization with time step $\Delta t={10}^{-2}s$, self-weight load of $q=79.17kN/m$. All cable elements are considered actives with time proportional increase from initial numerical straight length until reaching the full length. For each element of the central span we use $\Delta {\mathcal{l}}_c=0.63875m$ and for each element of the secondary spans we use $\Delta {\mathcal{l}}_{\mathcal{l}}=0.15m$.

With these elongation values, the normal stresses found (in the loaded phase) in the juxtaposed cable elements at the top of the towers are $S_{\mathcal{l}}=168.9MPa$ and $S_c=168.4MPa$, that is, almost equilibrated. In the second phase, keeping the loading and the initial stretching, the stress balancing $S_{\mathcal{l}}=S_c$ is imposed by Eq. (43) with a tolerance of $tol={10}^{-5}$. As the structure is stable, this phase is statically solved resulting in equal stresses $S_{\mathcal{l}}=S_c=168.48583MPa$ with a small length alteration for the connected cables at pillars, i.e., $\Delta {\mathcal{l}}_c=-0.0422m$ and $\Delta {\mathcal{l}}_{\mathcal{l}}=0.0423m$. Thus, the central cables become with total initial loaded length of ${\mathcal{l}}_c=2063.7906m$ and secondary cables with ${\mathcal{l}}_{\mathcal{l}}=1007.5423m$.

In the third phase, using the pseudo-dynamic regularization with the same number of time steps, we removed the self-weight loading to achieve the node coordinates to be used in the bridge discretization of Figs. 17 and 18 (initial natural position). As additional information, the sag at the central cable is 224.110m when loaded and e 221.311m unloaded.

ii) Static solution of the 2D and 3D models without adjusting the suspension cables

Using the natural coordinates of the isolated cable analysis, we generate the discretizations presented in Figs. 17 and 18 for the 2D and 3D models respectively. The 2D model has 1676 nodes and 2757 elements. The 3D model has 5986 nodes and 19405 elements. The 2D model girder loading is calculated from the loading of the 3D model to give the same overall weight for the complete structure. This load is equally distributed in all 2D girder elements as $q_t=20kN/m$ and in the pillars elements as $q_p=276.5kN/m$.

The 2D solution is performed in a single phase with pseudo dynamic regularization with only one time step of $\Delta t=1.0s$. The deformed configuration is presented in Fig. 17 with details for loose cables.

The 3D model is less stable than the 2D model in this phase of analysis. Thus, a restriction of displacement in the z direction is applied for all girder points until the equilibrium position is determined. A pseudo time interval of $\Delta t={10}^{-2}s$ is used along 100 steps. One may observe that the loading is applied proportionally to the time achieving 100% at the end of the analysis. The applied load is of self-weight plus 2.5% in the girder region, as previously described. The z displacements are released and a 3D static step is performed revealing the stability of the equilibrium configuration. The achieved results are shown in Fig. 18.

Figure 17
2D model vertical displacement profile, unit (m)

Figure 18
3D model vertical displacement profile, unit (m)

As it can be seen, the 2D and 3D vertical displacement profiles are similar. Contrary to what is expected, the 3D displacement is slightly higher than the 2D due to a redistribution of loading in the main cables resulting in less displacement near the pillars and extremities and greater at the central region. This redistribution is promoted by the stiffness of the girder that is secondary in the global stiffness of the structure.

iii) Adjusting the suspension cables to the proper positioning of the girder in its ideal static configuration (2D and 3D)

We impose a reduction of the suspension cables lengths equal to 1.2 times the girder displacement, immediately below the respective cable. This reduction is applied in just one pseudo dynamic step $\Delta t={10}^{-2}s$ with pseudo-dynamic regularization for both 2D and 3D models. Figures 19 and 20 show the new positions of the central girder with details for loose cables.

Figure 19
Vertical displacement profile for the 2D model after shortening of cables

Figure 20
Vertical displacement profile for the 3D model after shortening of cables

iv) Shortening of loose cables (2D and 3D)

In this item we apply stress control over loose cables. For both 2D and 3D models we take 7 cables from pillars through the main and secondary spans (for the 3D model it is a total of 28 cables) and impose that the stress at cables should be at least $S=25MPa$. For both models the relative tolerance in stress is $NTOL=0.02$ see Eq. (41). Figure 21 presents stress results for the two models in suspension cables near pillars. Figure 22 presents stresses at cables in the central region of the longest span. There is no noticeable change of overall displacements.

Figure 21
Stress values for cables near pillars in GPa.

Figure 22
Stress values at the center of the main span GPa

The maximum stress at main cables occurs near pillars and are $S_c=426MPa$ and $S_{\mathcal{l}}=425MPa$ for the 3D model. For the 2D model values are $S_c=407MPa$ and $S_{\mathcal{l}}=405MPa$ the difference is explained by the difference of girder stiffness.

Figure 23 shows an overview of the final stresses in the structure at its static configuration. Figure 24 shows a 3D view of the final static equilibrium configuration with all cables stretched.

Figure 23
Three-dimensional view of the stress distribution - GPa

Figure 24
Final configuration after stretching suspension cables (displacements (m)).

v) Transient dynamic analysis of the bridge under 3D real seismic shake.

This item shows the formulation's ability to simulate real structures subjected to actual seismic events. The chosen seismic event is that of Superstation Hill in 1987. The original unscaled time series record has been extracted from PEER Ground Motion Database (2014) with a time interval $\Delta t=5ms$ and is used in the FEM processing. Both vertical (UP) and horizontal (90) and (360) accelerations are transformed into displacements (m) and applied as base movements in the time domain. The horizontal base motion (90) is applied in the longitudinal direction of the bridge and the (360) in the transverse direction. As the bridge is very long we considered the propagation velocity of the earthquake approximately $V=2800m/s$ for all wave components. It introduces a difference of phase in the base motions, shown in Fig. 25. The wave arrives at one extremity of the bridge and follows its longitudinal direction (plane wave far from its origin).

Figure 25
Vertical base motion at the four supports.

Figure 26 shows the horizontal motion at the top of the first tower compared with the horizontal motion (earthquake input) pattern of the first support. This movement is the most significant for the tower.

Figure 26
First pillar horizontal displacement (bridge direction).

Figure 27 shows the vertical displacement at the center of the bridge (at one face) compared to the vertical movement pattern imposed on the first support. It is possible to estimate the natural frequency of the bridge from the presented motion $w_{vert}=1.518rad/s$ or $f_{vert}=0.242Hz$. The estimated period is $T=4.14s$

Figure 27
Vertical displacement at the bridge center

We observe that the displacement in the center of the bridge is different on each face, thus the difference between these displacements is shown in Fig. 28, which shows the excitation of the "torsion mode" of the bridge. Only as a reference value we show the vertical input in the first span support.

Figure 28
Difference of vertical displacements at the bridge center

In Fig. 29 we show the normal stress behavior for the two vertical cables in the center of the central span. As observed, the behavior presents a very high frequency not directly related to the central motion indicated in Fig. 27. However, a visual analysis of the mean tension reveals that this is following the main movement. Another detail is that there is no inversion of stress and, therefore, the cable is always tensioned.

Figure 29
Central suspension cables stress behavior.

Figure 30 shows the normal stress for the main cables (in the center of the longest span). In this case it is clear that the mean value of the stresses follows the bridge's vertical vibrating frequency. The largest difference between dynamic and static stress is approximately 7.5%.

Figure 30
Stress behavior of the main cables (center of main span)

Fig. 31 shows stresses for bars at the meeting between the pillar and the central span. The vertical pillar bar (with negative tension) has no stress inversion, while the horizontal connection bar (in the girder) has little stress inversion. However, it is a very robust bar and the stress level is very low. Both bars present high frequencies in the stress behavior, that is, the more significant displacement (vertical at the center of the bridge) does not induce high difference in stress at these points.

Figure 31
Pillar - girder stress level at truss bars.

After analyzing the vibration modes of the structure (see next sub-item) we chose to observe the vertical displacement of the central point of the first span, see Fig. 32. This motion shows that the critical point for vertical displacement of the analyzed structure is the central point of lateral (secondary) spans, not at central span as depicted in Fig. 27.

Figure 32
Vertical displacement at the center of the first span

vi) Determining the frequencies and modes of vibration of the positioned 3D structure

Taking the 3D structure in its final static (loaded) configuration, shown in Fig. 24, we analyze the natural frequencies and associated modes of vibration of the structure (case1) in the range $1.4171rad/s<w<1.5794rad/s$ in order to check the activated modes during the imposed earthquake of item iv, see text before Fig. 28. We achieved 73 vibration modes within the interval and, the most important ones (due to the amplitude of the girder motion when compared to cables) are depicted in Fig. 33. This result confirms the conclusion of both modes identified in the transient analysis, Figs. 27 and 28.

Figure 33
Two most important vibration modes with suspension cables

It is interesting to note the suspension cables movements that take place when the bridge is vibrating, Fig. 33. There are various intermediate frequencies dominated by suspension cables movements and there are also movements dominated by the main cables movements. In Fig. 34, omitting suspension cables, more detailed views of the two important modes are shown.

Figure 34
Vibration modes without suspension cables

In the analyzed interval there are modes presenting horizontal motion of main cables with small influence in the girder movement, see, for example, Figs. 35, 36 and 37.

Figure 35
Main cables superior view for $w_n=1.4841rad/s$

Figure 36
Main cables superior view for $w_n=1.5767rad/s$

Figure 37
Main cables superior view for $w_n=1.5768rad/s$

5. CONCLUSIONS

In this study we presented an alternative total Lagrangian finite element strategy to analyze cable structures, cable nets and mixed structures (cable/truss). Some functionalities as stress balance for contiguous elements and stress calibration for shortening loose cables are based on the presented active element idea. Other important feature of the finite element strategy is the presented pseudo-dynamic relaxation, which - eliminating inertial forces - reduces the presence of compression in finite element form finding analysis at the same time that improves the Hessian matrix conditioning.

The first two examples show that the formulation can analyze simple and complex cable and loose net problems. The adopted sparse matrix computation scheme shows that for nowadays computational resources the number of variables is small, even using truss elements to model cables. It also can be seen, in the third example, that the number of nodes related to cables is much smaller than the ones associated to the truss (girder) structure.

The third example (suspension bridge with a 2000m central span) demonstrates the applicability of the developed proposed active element. It shows the capacity and adequacy of the proposed formulation in modeling highly nonlinear structures, from its setting-up process through its dynamic response (earthquake excitation).

Using the proposed active finite element it is possible, in the setting-up analysis, the determination of the initial position of main cables, the stress balance at main cables through pillars, the stretching of loose suspension cables and the determination of the bridge form for utilization purposes.

The transient analysis is also successfully performed. A real seismic excitation is imposed on the structure and a preliminary value for the important vibration frequency could be determined. A study of the natural frequencies in an interval estimated by the transient analysis makes possible to find the important frequencies that was excited by the earthquake, improving the understanding and the future design of this kind of structure. Thus, the proposed FEM strategy is simple and general and can be used in the design process of outstanding engineering problems.

Acknowledgements

We would like to acknowledge Prof. Nicola Impollonia (University of Catania) for the valuable discussion on the second example and the Pacific Earthquake Engineering Research Center (PEER) for the seismic data. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

References

• Abad, M.S.A., Shooshtari, A., Esmaeili, V., Riabi, A.N., (2013). Nonlinear analysis of cable structures under general loadings, Finite Elements in Analysis and Design 7311-19.
• Ahmadizadeh, M., (2013). Three-dimensional geometrically nonlinear analysis of slack cable structures. Computers & Structures 128160-9.
• Andreu, A., Gil, L., Roca, P., (2006). A new deformable catenary element for the analysis of cable net structures, Computers & Structures 841882-90.
• Barnes, M.R., (1977). Form-finding and analysis of tension space structures by dynamic relaxation, Ph.D. thesis, City University London, United Kingdom.
• Borst, R., Crisfield, M., Remmers, J., Verhoosel, C., (2012). Nonlinear finite element analysis of solids and structures. 2nd ed. John Wiley & Sons, Ltd.
• Clough, R.W., Penzien, J., (1993). Dynamics of structures, 2nd edition, McGraw Hill, New York.
• Coda, H.B., (2009). Two dimensional analysis of inflatable structures by the positional FEM, Latin American Journal of Solids and Structures 6(3)187-212.
• Coda, H.B., (2015). Continuous inter-laminar stresses for regular and inverse geometrically nonlinear dynamic and static analyses of laminated plates and shells, Composite Structures 132406-422.
• Coda, H.B., (2018). Positional Finite Element Solids and structures - Geometrical nonlinearity and dynamics (in Portuguese), São Carlos EESC-USP.
• Crusells-Girona, M., Filippou, F.C., Taylor, R.L., (2017). A mixed formulation for nonlinear analysis of cable structures, Computers and Structures 18650-61.
• Furuya, N., Yamaoka, R., Paulson Jr, B.C., (1994). Construction of Akashi-Kaikyo Bridge West Anchorage, Journal of Construction Engineering and Management 120(2)337-356.
• Greco, L., Cuomo, M., (2012). On the force density method for slack cable nets, International Journal of Solids and Structures 49(13)1526-1540.
• Greco, L., Impollonia, N., Cuomo, M., (2014). A procedure for the static analysis of cable structures following elastic catenary theory, International Journal of Solids and Structures 511521-1533.
• Greco, M., Ferreira, I.P., Barros, F.B., (2013). A classical time integration method applied for solution of nonlinear equations of a double-layer tensegrity, Journal of the Brazilian Society of Mechanical Sciences and Engineering 35(1)41-50.
• 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.
• Impollonia, N., Ricciardi, G., Saitta, F., (2011). Statics of elastic cables under 3D point forces, International Journal of Solids and Structures 481268-1276.
• Kim, H.K., Lee, M.J., Chang, S.P., (2002). Nonlinear shape-finding analysis of a self-anchored suspension bridge, Engineering Structures 241547-1559.
• Kim, N.-I., Thai, S., Lee, J., (2016). Nonlinear elasto-plastic analysis of slack and taut cable structures. Engineering with Computers 32615-627.
• Lewis, W., Jones, M., Rushton, K., (1984). Dynamic relaxation analysis of the nonlinear static response of pretensioned cable roofs, Computers & Structures 18(6)989-997.
• Miyata, T., Yamaguchi, K., (1993). Aerodynamics of wind effects on the Akashi-Kaikyo Bridge, Journal of Wi Engineering and Industrial Aerodynamics 48287-315.
• Pascon, J.P., Coda, H.B., (2013). High-order tetrahedral finite elements applied to large deformation analysis of functionally graded rubber-like materials. Applied Mathematical Modelling 37(20-21)8757-8775.
• Pascon, J.P., Coda, H.B., (2015). Large deformation analysis of functionally graded elastoplastic materials via solid tetrahedral finite elements. Computers & Structures 14659-75.
• PEER Ground Motion Database, Pacific Earthquake Engineering Research Center. 2014 (accessed March 2014)
• Schek, H.-J., (1974). The force density method for form finding and computation of general networks. Computer Methods in Applied Mechanics and Engineering 3115-134.
• Siev, A., Eidelman, J., (1964). Stress analysis of prestressed suspended roofs, Journal of the Structural Division, Proceedings of the American Society of Civil Engineers 103-121.
• Soares, H.B., Paccola, R.R., Coda, H.B., (2019). Unconstrained vector positional shell FEM formulation applied to thin-walled members instability analysis. Thin-Walled Structures 136246-257.
• Such, M., Jimenez-Octavio, J., Carnicero, A., Lopez-Garcia, O., (2009). An approach based on the catenary equation to deal with static analysis of three dimensional cable structures, Engineering Structures 312162-70.
• Thai, H.-T., Kim, S.-E., (2011). Nonlinear static and dynamic analysis of cable structures, Finite Elements in Analysis and Design 47(3)237-246.
• Thai, H.-T., Kim, S.-E., (2012). Second-order inelastic analysis of cable-stayed bridges, Finite Elements in Analysis and Design 5348-55.
• The Harwell Subroutine Library Mathematical Software Library - HSL, (2014). A collection of Fortran codes for large-scale scientific computation. http//www.hsl.rl.ac.uk (accessed March 2014).
  » http//www.hsl.rl.ac.uk
• Veenendaal, D., Block, P., (2012), An overview and comparison of structural form finding methods for general networks, International Journal of Solids and Structures 493741-3753.
• Warburton, G.B., (1976). The Dynamical Behaviour of Structures, Pergamon Press.
• Wüchner, R., Bletzinger, K.U., (2005). Stress adapted numerical form finding of prestressed surfaces by the updated reference strategy. International Journal for Numerical Methods in Engineering 64143-166.
• Yang, Y.B., Tsay, J.Y., (2007). Geometric nonlinear analysis of cable structures with a twonode cable element by generalized displacement control method, International Journal of Structural Stability and Dynamics 7(4)571-588.
• Yim, W.T., (2007). The bridge engineering 2 Conference Akashi Bridge, Proceedings of Bridge Engineering 2 Conference 2007, University of Bath, Bath, UK.
• Zhang, J., Au, F.T.K., (2014). Calibration of initial cable forces in cable-stayed bridge based on Kriging approach, Finite Elements in Analysis and Design 9280-92.
• Zhou, X., Li, J., Liu, J., Chen, Y.F., (2017). Dynamic Performance Characteristics of Pre-Stressed Cable RC Truss Floor System under Human-Induced Loads, International Journal of Structural Stability and Dynamics 17(4)17500491-20.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  20 May 2020
• Date of issue
  2020

History

• Received
  03 Oct 2019
• Reviewed
  23 Mar 2020
• Accepted
  31 Mar 2020
• Published
  09 Apr 2020