Nonlinear static and dynamic analysis of mixed cable elements - View source for Crusells-Girona 2017a

==NONLINEAR STATIC AND DYNAMIC ANALYSIS OF MIXED CABLE ELEMENTS==

'''MIQUEL CRUSELLS-GIRONA<math>^{*}</math>, 

FILIP C. FILIPPOU<math>^{*}</math> AND ROBERT L. TAYLOR<math>^{*}</math>'''

{|
|-
|<math>^{*}</math> Department of Civil and Environmental Engineering
|-
| University of California, Berkeley
|-
| 94720 Berkeley, CA, USA
|-
| e-mail: [mailto:miquel.crusells@berkeley.edu miquel.crusells@berkeley.edu]
|}

==Abstract==

This paper presents a family of finite elements for the nonlinear static and dynamic analysis of cables based on a mixed variational formulation in curvilinear coordinates and finite deformations. This formulation identifies stress measures, in the form of axial forces, and conjugate deformation measures for the nonlinear catenary problem. The continuity requirements lead to two distinct implementations: one with a continuous axial force distribution and one with a discontinuous. Two examples from the literature on nonlinear cable analysis are used to validate the proposed formulation for St Venant-Kirchhoff elastic materials. These studies show that displacements and axial forces are captured with high accuracy for both the static and the dynamic case.

'''keywords''' Cables, Mixed FEM, Weak Compatibility, Nonlinear Dynamics

==1 INTRODUCTION==

Cable structures are of great interest in many engineering applications because they offer numerous advantages, such as high ultimate strength, light weight or prestressing capabilities, among others. Nonetheless, a highly nonlinear behavior arises in this type of structures because of their high flexibility. For analyzing cable structures, two families of elements have traditionally been considered: truss elements and catenary elements.

For truss elements, the cable is discretized in a series of straight 2-node elements. In this case, the geometric nonlinearity is often accounted for by a corotational formulation, involving the transformation of the node kinematic variables under large displacements. Truss elements suffer from excessive mesh refinement to obtain accurate results, especially when assuming a constant axial force distribution in the element. Moreover, they may exhibit snap-through instabilities at states of nearly singular stiffness.

Catenary elements use linear kinematics to discretize the cable into a series of curved elements that satisfy the catenary equation. These elements solve the global balance of linear momentum by explicit integration and assuming linear elasticity <span id='citeF-1'></span>[[#cite-1|[1]]]. As a result, loads are not adjusted with the cable elongation, so that these elements cannot be extended to nonlinear elasticity or inelasticity. Recently, the authors have proposed a general formulation for a catenary element in finite deformations and curvilinear coordinates <span id='citeF-2'></span>[[#cite-2|[2]]] that overcomes these limitations.

==2 MIXED FORMULATION OF THE CATENARY PROBLEM==

===2.1 Kinematics===

Fig.&nbsp;[[#img-1|1]] shows the motion of a cable from a reference configuration <math display="inline">\mathcal{P}_0</math> to a current configuration <math display="inline">\mathcal{P}</math>. Define an orthogonal frame <math display="inline">\{ \mathbf{G}_i\} _{i=1}^3</math> with associated coordinates <math display="inline">\{ \xi ^i\} _{i=1}^3</math> at any point <math display="inline">P\in \mathcal{P}_0</math>, such that

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{G}_1 = \frac{d\mathbf{X}}{d\xi ^1} \; \; ; \; \; \mathbf{G}_1\cdot \mathbf{G}_2 = 0 \; \; ; \; \; \| \mathbf{G}_2\|=1 \; \; ; \; \; \mathbf{G}_3=\frac{\mathbf{G}_1 \times \mathbf{G}_2}{\| \mathbf{G}_1 \times \mathbf{G}_2\| } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">\xi ^1</math> is the parameter describing the curve. Under the motion <math display="inline">\mathbf{x}=\boldsymbol{\chi }(\mathbf{X})</math>, this frame is convected to the orthogonal frame <math display="inline">\{ \mathbf{g}_i\} _{i=1}^3</math>. Let upper case letters denote variables in the reference configuration and lower case letters, variables in the current configuration.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Crusells-Girona_842065898-motion.png|792px|Motion \mathbfx=χ(\mathbfX(ξ¹)) of the cable \mathcalC.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Motion <math>\mathbf{x}=\boldsymbol{\chi }(\mathbf{X}(\xi ^1))</math> of the cable <math>\mathcal{C}</math>.
|}

The relevant stretch and Green-Lagrange strain of the problem, in the <math display="inline">\mathbf{g}_1</math> direction, are

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\lambda = \frac{\| \mathbf{g}_1\| }{\| \mathbf{G}_1\| } \;\;\;\;\;\; ; \;\;\;\;\;\; \mathrm{E}=\frac{1}{2}(\lambda ^2-1)\| \mathbf{G}_1\| ^2 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

The displacement vector <math display="inline">\mathbf{u}</math> depends only on the curvilinear coordinate <math display="inline">\xi ^1</math>,

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{u}(\mathbf{X}(\xi ^1))=\mathbf{x}(\mathbf{X}(\xi ^1))-\mathbf{X}(\xi ^1) = u_A(\xi ^1) \mathbf{E}_A </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

Therefore, the relationship between the displacement field <math display="inline">\mathbf{u}</math> and the relevant Green-Lagrange strain <math display="inline">\mathrm{E}</math> can be computed <span id='citeF-2'></span>[[#cite-2|[2]]] as

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{E}=\frac{d \mathbf{u}}{d \xi ^1}\cdot \mathbf{G}_1+\frac{1}{2}\left|\frac{d \mathbf{u}}{d \xi ^1}\right|^2=\frac{d \mathbf{u}}{d \xi ^1} \cdot \left(\mathbf{G}_1 +\frac{1}{2}\frac{d \mathbf{u}}{d \xi ^1} \right)=\frac{1}{2}\frac{d \mathbf{u}}{d \xi ^1}\cdot \left(\mathbf{G}_1+\mathbf{g}_1\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

It is relevant to observe that the frame <math display="inline">\{ \mathbf{G}_i\} _{i=1}^3</math> is orthogonal, but not orthornomal in general. Indeed, the metric tensor <math display="inline">G_{ij}=\mathbf{G}_i \cdot \mathbf{G}_j</math> is not necessarily the identity operator and the Green-Lagrange strain <math display="inline">\mathrm{E}</math> may not be physical. Nevertheless, one can construct an orthonormal basis <math display="inline">\{ \mathbf{\hat{G}}_A\} _{A=1}^3=\{ \mathbf{G}_i/\| \mathbf{G}_i\| \} _{i=1}^3</math> such that

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{E}=\mathrm{E}_{ij} \mathbf{G}^i \otimes \mathbf{G}^j = \mathrm{\hat{E}}_{AB}\mathbf{\hat{G}}_A \otimes \mathbf{\hat{G}}_B </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

Hence, the components

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{\hat{E}}_{AB}=\mathrm{E}_{ij}(\mathbf{\hat{G}}_A \cdot \mathbf{G}^i)(\mathbf{\hat{G}}_B \cdot \mathbf{G}^j) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

are physical quantities.

===2.2 Equilibrium and principle of virtual work===

For expressing the equilibrium equation of the cable in finite deformations, let <math display="inline">\mathbf{n}</math> denote the axial force in the current configuration, thus a Cauchy representation. Observe that the first Piola-Kirchhoff  and the Cauchy representations of the axial force coincide for the problem in hand, which does not account for changes in the cross section dimensions.

The Cauchy axial force can be pulled back to the reference configuration to obtain a second Piola-Kirchhoff representation of the axial force, <math display="inline">\mathbf{N}</math>. It can be shown <span id='citeF-2'></span>[[#cite-2|[2]]] that, using the orthonormal basis in Eq.&nbsp;[[#eq-5|5]], namely with components <math display="inline">\mathbf{n}=\mathrm{\hat{n}}\mathbf{\hat{g}}_1</math> and <math display="inline">\mathbf{N}=\mathrm{\hat{N}}\mathbf{\hat{G}}_1</math>,

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{\hat{n}}=\lambda \mathrm{\hat{N}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

Denoting by <math display="inline">s</math> and <math display="inline">S</math> the arc-length coordinates in the current and reference configurations, respectively, the cable distributed load can be described as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{w}ds = \mathbf{\hat{W}}dS = \mathbf{W}d\xi ^1 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

Then, global equilibrium for the cable in the current configuration <math display="inline">\mathcal{P}</math> states

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{n}(s)-\mathbf{n}(0)+\int _0^s\mathbf{w}\, ds = \int _0^s \rho \mathbf{a}\, ds </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

where <math display="inline">\rho </math> is the material density in the current configuration and <math display="inline">\mathbf{a}</math>, the total acceleration. The corresponding local statement in the current configuration can be obtained with the fundamental theorem of calculus and the localization theorem

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{d}{ds}(\mathrm{n}\mathbf{g}_1)+\mathbf{w}=\rho \mathbf{a} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

or, in the reference configuration,

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{d}{dS}(\mathrm{\hat{N}}\sqrt{G^{11}}\mathbf{g}_1)+\mathbf{\hat{W}}=\rho _0\mathbf{a} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

where <math display="inline">G^{ij}=\mathbf{G}^i \cdot \mathbf{G}^j</math> represents the dual metric tensor and <math display="inline">\rho _0 = \lambda \rho </math>.

In summary, if <math display="inline">\mathrm{\hat{N}}=\Psi (\mathrm{\hat{E}})</math> is a frame-indifferent constitutive relation between the phy-sical Green-Lagrange strain <math display="inline">\mathrm{\hat{E}}</math> and the physical 2nd Piola-Kirchhoff axial force <math display="inline">\mathrm{\hat{N}}</math>, the pair of fields (<math display="inline">\mathbf{u}</math>, <math display="inline">\mathrm{\hat{N}}</math>) will be the solution of the cable problem, if and only if, they satisfy

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\{ \begin{array}{rl}\displaystyle G^{11}\frac{d\mathbf{u}}{d\xi ^1} \cdot \left(\mathbf{G}_1+\frac{1}{2}\frac{d\mathbf{u}}{d\xi ^1}\right)-\mathrm{\hat{E}}=0 & \mathrm{in} \; \Omega=(0,L) \\[3.5mm] \displaystyle \frac{d}{dS}\left(\sqrt{G^{11}}\mathrm{\hat{N}}\mathbf{g}_1 \right)+\mathbf{\hat{W}}=\rho _0 \mathbf{a}  & \mathrm{in} \; \Omega=(0,L)\\[] \mathrm{\hat{N}}-\Psi (\mathrm{\hat{E}})=0 & \mathrm{in} \; \Omega=(0,L) \\[] \mathbf{u} = \mathbf{\bar{u}} & \mathrm{on} \; \Gamma _u \\[1mm] \displaystyle \sqrt{G^{11}}\mathrm{\hat{N}}\mathbf{g}_1= \mathbf{\bar{T}} & \mathrm{on} \; \Gamma _q \end{array} \right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

for <math display="inline">0<S<L</math> equivalent to <math display="inline">\xi ^1_1<\xi ^1<\xi ^1_2</math>.

The corresponding two-field weak statement of Eq.&nbsp;[[#eq-12|12]] can be obtained by considering any variation <math display="inline">\delta \mathbf{u}\in \mathcal{V}</math>, the space of displacement test functions, any variation <math display="inline">\delta \mathrm{\hat{N}}\in \mathcal{W}</math>, the space of axial force test functions, and integrating the equilibrium equation by parts,

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\{ \begin{array}{r}\displaystyle \int _0^L \delta \mathrm{\hat{N}}\left\{G^{11}\frac{d\mathbf{u}}{d\xi ^1} \cdot \left(\mathbf{G}_1+\frac{1}{2}\frac{d\mathbf{u}}{d\xi ^1}\right)-\mathrm{\hat{E}}\right\}dS=0 \\[3.6mm] \displaystyle \int _0^L \frac{d(\delta \mathbf{u})}{dS} \cdot \mathrm{\hat{N}}\sqrt{G^{11}}\mathbf{g}_1 \, dS + \int _0^L \delta \mathbf{u}\cdot \rho _0 \mathbf{a}\, dS =\left[\delta \mathbf{u}\cdot \mathbf{\bar{T}} \right]_{\Gamma _q}+\int _0^L \delta \mathbf{u}\cdot \mathbf{\hat{W}}\, dS \end{array} \right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where the constitutive relation is imposed strongly. The spaces for the trial solutions of the displacements and axial forces, <math display="inline">\mathcal{S}</math> and <math display="inline">\mathcal{N}</math>, respectively, are

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\mathcal{S}&=\{ \mathbf{u}\in H^1(0,L) \, | \,  \mathbf{u}=\mathbf{\bar{u}} \; \mathrm{on} \; \Gamma _u\} \\ \mathcal{N}&= \left\{\mathrm{\hat{N}}\in H^0(0,L) \, | \,  \mathrm{\hat{N}}>0, \; \mathrm{and} \; \mathrm{\hat{N}}=g^{11}\sqrt{G_{11}}\mathbf{\bar{T}}\cdot \mathbf{g}_1 \; \mathrm{on} \; \Gamma _q \right\} \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

Similarly, the spaces for the test functions of the displacements and the axial forces, <math display="inline">\mathcal{V}</math> and <math display="inline">\mathcal{W}</math>, respectively, become

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\mathcal{V}&=\{ \delta \mathbf{u}\in H^1(0,L) \, | \,  \delta \mathbf{u}=0 \; \mathrm{on} \; \Gamma _u\} \\ \mathcal{W}&=\{ \delta \mathrm{\hat{N}} \in H^0(0,L) \, | \, \delta \mathrm{\hat{N}}=0 \; \mathrm{on} \; \Gamma _q\}  \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

where <math display="inline">H^k(\Omega )</math> is the Sobolev space for the <math display="inline">k-</math>th weak derivative in the <math display="inline">L^2(\Omega )</math> norm.

As a result, there are no continuity requirements for the axial force field <math display="inline">\mathrm{\hat{N}}</math>. This implies the possibility of exploring cable finite elements with continuous or discontinuous axial force distribution.

==3 FINITE-ELEMENT IMPLEMENTATION==

===3.1 Discretization===

The discretization of the governing equations requires interpolations for <math display="inline">\mathrm{\hat{N}}(\xi ^1)</math>, <math display="inline">\mathbf{u}(\xi ^1)</math> and <math display="inline">\mathbf{a}(\xi ^1)</math>. Assume a <math display="inline">k</math>-th order Galerkin interpolation for the axial forces

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{\hat{N}}=\boldsymbol{\varphi }^t\mathbf{\hat{N}}=\mathbf{\hat{N}}^t \boldsymbol{\varphi } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

and an <math display="inline">l</math>-th order Galerkin interpolation for the displacement and acceleration fields

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{u}=\boldsymbol{\phi }^t \mathbf{\hat{u}} \;\;\;\;\;\; ; \;\;\;\;\;\; \mathbf{a}=\boldsymbol{\phi }^t \mathbf{\hat{a}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

Then, using the same shape functions for the reference configuration, the current configu-ration is obtained as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}=\mathbf{X}+\mathbf{u}=\boldsymbol{\phi }^t (\mathbf{\hat{X}}+\mathbf{\hat{u}})=\boldsymbol{\phi }^t \mathbf{\hat{x}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

With these interpolation functions, one can discretize the weak statement in Eq.&nbsp;[[#eq-13|13]] for a finite element <math display="inline">\Omega _e</math> and a time step <math display="inline">n</math> as

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\{ \begin{array}{r}\displaystyle \int _{\Omega _e} \delta \mathbf{\hat{N}}^t \boldsymbol{\varphi } \left\{G^{11}\mathbf{\hat{u}}^t_n \boldsymbol{\phi }' \left(\mathbf{G}_1+\frac{1}{2}(\boldsymbol{\phi }')^t\mathbf{\hat{u}}_n\right)-\mathrm{\hat{E}}(\mathbf{\hat{N}}_n)\right\}dS=0 \\[3.6mm] \displaystyle \int _{\Omega _e} \delta \mathbf{\hat{u}}^t \boldsymbol{\phi }' G^{11} \boldsymbol{\varphi }^t \mathbf{\hat{N}}_n \mathbf{\hat{g}}_n \, dS + \int _{\Omega _e} \delta \mathbf{\hat{u}}^t \rho _0 \boldsymbol{\phi } \boldsymbol{\phi }^t\mathbf{\hat{a}}_n\, dS =\left[\delta \mathbf{\hat{u}}^t\boldsymbol{\phi } \mathbf{\bar{T}}_n \right]_{\partial \Omega _e}+\int _{\Omega _e} \delta \mathbf{\hat{u}}^t \boldsymbol{\phi } \mathbf{\hat{W}}_n\, dS \end{array} \right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

where <math display="inline">(\cdot )'</math> represents the derivative with respect to the curvilinear coordinate <math display="inline">\xi ^1</math> and <math display="inline">\mathbf{\hat{g}}_n=\mathbf{G}_1+(\boldsymbol{\phi }')^t\mathbf{\hat{u}}_n</math> is the numerical counterpart to <math display="inline">\mathbf{g}_1</math>.

===3.2 Time integration and consistent linearization===

Once the discretization of the problem has been performed, the corresponding time-dependent equations need to be solved. As stated before, one can consider cable finite elements with a continuous or a discontinuous axial force field.

====3.2.1 Mixed cable element with continuous axial force====

For the element with continuous axial force distribution, the cable is subdivided into <math display="inline">e</math> elements of <math display="inline">k</math>-th order in axial forces and <math display="inline">l</math>-th order in displacements. By defining the expanded stress divergence term <math display="inline">\mathbf{R}=(\mathbf{R}_1,\mathbf{R}_2)</math> with components

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\mathbf{R}_1 (\mathbf{\hat{N}}_n, \mathbf{\hat{u}}_n) &= \int _{\Omega _e} \boldsymbol{\varphi } \left(G^{11} \mathbf{\hat{u}}^t_n\boldsymbol{\phi }' \left(\mathbf{G}_1+\frac{1}{2}(\boldsymbol{\phi }')^t\mathbf{\hat{u}}_n\right)-\mathrm{\hat{E}}(\mathbf{\hat{N}}_n)\right)dS \\[0.3em] \mathbf{R}_2 (\mathbf{\hat{N}}_n, \mathbf{\hat{u}}_n) &= \int _{\Omega _e} G^{11} \boldsymbol{\varphi }^t\mathbf{\hat{N}}_n \, \boldsymbol{\phi' } \mathbf{\hat{g}}_n \, dS \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

and the mass matrix <math display="inline">\mathbf{M}</math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M}=\int _{\Omega _e} \rho _0 \boldsymbol{\phi }\boldsymbol{\phi }^t\, dS </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

one can rewrite Eq.&nbsp;[[#eq-21|21]] in an implicit scheme as

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[ \begin{array}{c}\mathbf{R}_1(\mathbf{\hat{N}}_{n+1},\mathbf{\hat{u}}_{n+1}) \\ \mathbf{R}_2(\mathbf{\hat{N}}_{n+1},\mathbf{\hat{u}}_{n+1})  \end{array} \right]+ \left[ \begin{array}{c}\mathbf{0} \\ \mathbf{M}\mathbf{\hat{a}}_{n+1} \end{array} \right]= \left[ \begin{array}{c}\mathbf{0} \\ \mathbf{F}_{ext,n+1} \end{array} \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

where <math display="inline">\mathbf{F}_{ext,n+1}</math> refers to the external forces considered at the time step <math display="inline">n+1</math>.

Introducing Newmark's time integrator <span id='citeF-3'></span>[[#cite-3|[3]]], one obtains the system of equations

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\frac{1}{\beta \Delta t_n^2}\left[ \begin{array}{c}\mathbf{0} \\ \mathbf{M}\mathbf{\hat{u}}_{n+1} \end{array} \right]&+ \left[ \begin{array}{c}\mathbf{R}_1(\mathbf{\hat{N}}_{n+1},\mathbf{\hat{u}}_{n+1}) \\ \mathbf{R}_2(\mathbf{\hat{N}}_{n+1},\mathbf{\hat{u}}_{n+1})  \end{array} \right] \\ &= \left[ \begin{array}{c}\mathbf{0} \\ \mathbf{F}_{ext,n+1} \end{array} \right]+ \frac{1}{\beta \Delta t_n^2}\left[ \begin{array}{c}\mathbf{0} \\ \mathbf{M}(\mathbf{\hat{u}}_n + \Delta t_n \mathbf{\hat{v}}_n) \end{array} \right]+ \frac{1-2\beta }{2\beta }\left[ \begin{array}{c}\mathbf{0} \\ \mathbf{M}\mathbf{\hat{a}}_n \end{array} \right] \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

Hence the consistent linearization of the former equation, namely <math display="inline">\boldsymbol{\Phi }(\mathbf{\hat{u}}_{n+1},\mathbf{\hat{N}}_{n+1})=\mathbf{0}</math>, around a point <math display="inline">\mathbf{\bar{V}}_{n+1}=(\mathbf{\hat{u}}_{n+1},\mathbf{\hat{N}}_{n+1})</math> and for the <math display="inline">k</math>-th iterate establishes

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L}\boldsymbol{\Phi }=\boldsymbol{\Phi }|_{\mathbf{\bar{V}}_{n+1}}^{(k)}+\underbrace{\left. \frac{\partial \mathbf{\Phi }}{\partial \mathbf{\hat{V}}_{n+1}} \right|_{\mathbf{\bar{V}}_{n+1}^{(k)}} ( \mathbf{\hat{V}}_{n+1}^{(k+1)}-\mathbf{\bar{V}}_{n+1}^{(k)})}_{\mathrm{D}\boldsymbol{\Phi }(\mathbf{\bar{V}}_{n+1}^{(k)} ,\Delta \mathbf{V}_{n+1})}=\mathbf{0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

where the Fréchet derivative <math display="inline">\partial \boldsymbol{\Phi }/\partial \mathbf{\hat{V}}_{n+1}|_{\mathbf{\hat{V}}_{n+1}^{(k)}}</math> corresponds to the dynamic stiffness <math display="inline">\mathbf{K}</math> of the problem, with components

<span id="eq-34"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{2}   \mathbf{K}_{\mathbf{N}\mathbf{N}} = -\int _{\Omega _e} \boldsymbol{\varphi } \frac{\partial \mathrm{\hat{E}}}{\partial \mathbf{\hat{N}}_{n+1}} \, dS = -\int _{\Omega _e} \boldsymbol{\varphi } \frac{\partial \mathrm{\hat{E}}}{\partial \mathrm{\hat{N}}} \boldsymbol{\varphi }^t \, dS </math>
|-
| style="text-align: center;" | <math>   \mathbf{K}_{\mathbf{N}\mathbf{u}} =  \int _{\Omega _e} G^{11}\boldsymbol{\varphi }\mathbf{\hat{g}}^t_{n+1} (\boldsymbol{\phi' })^t  dS =  \mathbf{K}_{\mathbf{u}\mathbf{N}}^t </math>
|-
| style="text-align: center;" | <math>   \mathbf{K}_{\mathbf{u}\mathbf{u}}^s=  \int _{\Omega _e} G^{11} \boldsymbol{\varphi }^t \mathbf{\hat{N}}_{n+1}\boldsymbol{\phi' }(\boldsymbol{\phi' })^t dS </math>
|-
| style="text-align: center;" | <math>   \mathbf{K}_{\mathbf{u}\mathbf{u}}^d=  \frac{1}{\beta \Delta t_n^2}\mathbf{M}+ \mathbf{K}_{\mathbf{u}\mathbf{u}}^s </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

in the form

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{K}=\left.\frac{\partial \boldsymbol{\Phi }}{\partial \mathbf{\hat{V}}_{n+1}}\right|_{\mathbf{\bar V}_{n+1}^{(k)}}=\left[\begin{array}{cc}\mathbf{K}_{\mathbf{N}\mathbf{N}} & \mathbf{K}_{\mathbf{N}\mathbf{u}} \\ \mathbf{K}_{\mathbf{u}\mathbf{N}} & \mathbf{K}_{\mathbf{u}\mathbf{u}}^d \end{array} \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

In order to satisfy stability of the solution scheme, it is necessary <span id='citeF-2'></span>[[#cite-2|[2]]] that

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\ker \, (\mathbf{K}_{\mathbf{N}\mathbf{N}}-\mathbf{K}_{\mathbf{N}\mathbf{u}}(\mathbf{K}_{\mathbf{u}\mathbf{u}}^s)^{-1}\mathbf{K}_{\mathbf{N}\mathbf{u}}^t)=\mathbf{0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

====3.2.2 Mixed cable element with discontinuous axial force====

For the element with discontinuous axial force distribution, the cable is also subdivided into <math display="inline">e</math> elements of <math display="inline">k</math>-th order in axial forces and <math display="inline">l</math>-th order in displacements. In this case, however, the axial forces are treated as internal degrees of freedom, and are consequently condensed out at the element level before assembly of the element response. This generates a discontinuity in the axial forces, which is allowed by the condition <math display="inline">\mathrm{\hat{N}}\in H^0(0,L)</math>. The stress divergence term <math display="inline">\mathbf{R}(\mathbf{\hat{N}}_n,\mathbf{\hat{u}}_n)</math> is then understood as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{R}(\mathbf{\hat{N}}_n(\mathbf{\hat{u}}_n),\mathbf{\hat{u}}_n)=\int _{\Omega _e} G^{11}\boldsymbol{\varphi }^t\mathbf{\hat{N}}_n \, \boldsymbol{\phi' } \mathbf{\hat{g}}_n \, dS </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

and one can rewrite Eq.&nbsp;[[#eq-21|21]] in an implicit scheme as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{R}(\mathbf{\hat{u}}_{n+1})+\mathbf{M}\mathbf{\hat{a}}_{n+1}=\mathbf{F}_{ext,n+1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

Introducing Newmark's time integrator <span id='citeF-3'></span>[[#cite-3|[3]]], one obtains the system of equations

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{1}{\beta \Delta t_n^2}\mathbf{M}\mathbf{\hat{u}}_{n+1}+\mathbf{R}(\mathbf{\hat{u}}_{n+1})=\mathbf{F}_{ext,n+1}+\frac{1}{\beta \Delta t_n^2}\mathbf{M}(\mathbf{\hat{u}}_n+\Delta t_n \mathbf{\hat{v}}_n)+\frac{1-2\beta }{2\beta }\mathbf{M}\mathbf{\hat{a}}_n </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

Hence the consistent linearization of the former equation, namely <math display="inline">\boldsymbol{\Phi }(\mathbf{\hat{u}}_{n+1})=\mathbf{0}</math>, around a point <math display="inline">\mathbf{\bar{u}}_{n+1}</math> and for the <math display="inline">k</math>-th iterate establishes

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L}\boldsymbol{\Phi }=\boldsymbol{\Phi }|_{\mathbf{\bar{u}}_{n+1}^{(k)}}+\underbrace{\left.\frac{\partial \boldsymbol{\Phi }}{\partial \mathbf{\hat{u}}_{n+1}}\right|_{\mathbf{\bar{u}}_{n+1}^{(k)}}(\mathbf{\hat{u}}_{n+1}^{(k+1)}-\mathbf{\bar{u}}_{n+1}^{(k)})}_{\mathrm{D}\boldsymbol{\Phi }(\mathbf{\bar{u}}_{n+1}^{(k)},\Delta \mathbf{u}_{n+1})}=\mathbf{0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

where the Fréchet derivative <math display="inline">\partial \boldsymbol{\Phi }/\partial \mathbf{\hat{u}}_{n+1}|_{\mathbf{\bar{u}}_{n+1}^{(k)}}</math> corresponds to the condensed dynamic stiffness <math display="inline">\mathbf{K}</math> of the problem

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{K}=\mathbf{K}_{\mathbf{u}\mathbf{u}}^d-\mathbf{K}_{\mathbf{u}\mathbf{N}}\mathbf{K}_{\mathbf{N}\mathbf{N}}^{-1}\mathbf{K}_{\mathbf{N}\mathbf{u}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

with the components defined in Eq.&nbsp;[[#eq-34|34]].

The stability condition of the solution scheme reads <span id='citeF-2'></span>[[#cite-2|[2]]] in this case as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\ker \, (\mathbf{K}_{\mathbf{u}\mathbf{u}}^s-\mathbf{K}_{\mathbf{u}\mathbf{N}}\mathbf{K}_{\mathbf{N}\mathbf{N}}^{-1}\mathbf{K}_{\mathbf{N}\mathbf{u}})=\mathbf{0} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

==4 NUMERICAL EXAMPLES==

The proposed formulation is implemented in two cable elements with continuous and discontinuous axial force distributions. The elements, deployed in the general purpose finite element program FEAP <span id='citeF-4'></span>[[#cite-4|[4]]] and Matlab toolbox FEDEASLab <span id='citeF-5'></span>[[#cite-5|[5]]], use a linear approximation for the axial forces (<math display="inline">k=1</math>) and a quadratic approximation for the displacements <math display="inline">(l=2)</math>. The two-dimensional element results in eight degrees of freedom (DOFs), six displacement DOFs and two axial force DOFs, while the three-dimensional element results in eleven DOFs, nine displacement DOFs and two axial force DOFs.

Both elements are implemented with a St Venant - Kirchhoff elastic material model with stored energy <math display="inline">\mathcal{U}</math> in terms of the stretch <math display="inline">\lambda </math> and the generalized Young's modulus <math display="inline">E</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{U}=\frac{E}{8}(\lambda ^2-1)^2 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

Thus, if <math display="inline">A</math> is the area of the cross section and <math display="inline">\mathrm{\hat{N}}_0</math> the prestressing force,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathrm{\hat{N}}-\mathrm{\hat{N}}_0=(EA)\mathrm{\hat{E}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

with constant material stiffness <math display="inline">\partial \mathrm{\hat{N}}/\partial \mathrm{\hat{E}}=EA</math>.

===4.1 Example 1: Stability of a 3d pulley system===

The first example investigates the stability of a 3d cable supported by a pulley that was previously studied by Impollonia ''et al'' <span id='citeF-6'></span>[[#cite-6|[6]]]. The structural model, whose geometric and material properties are shown in Table&nbsp;[[#table-1|1]], consists of a cable anchored at both ends and supported by an intermediate roller. In this case, the inertia forces in Eq.&nbsp;[[#eq-21|21]] and the mass term of the stiffness are not considered as the problem is analyzed in a static manner.
{|  class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;"
|+ <span id='table-1'></span>Table. 1 <span style="text-align: center; font-size: 75%;">Geometry and Material Properties for Example 1.</span>
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  '''Property''' 
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | '''Value''' 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Cross-sectional area 
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 805 mm<math display="inline">^2</math> 
|-
| style="border-left: 2px solid;border-right: 2px solid;" | Elastic modulus 
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 16.0 kN/mm<math display="inline">^2</math> 
|-
| style="border-left: 2px solid;border-right: 2px solid;" | Cable self-weight 
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 62.0679 N/m 
|- style="border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | Cable length 
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 500 m 

|}

The objective of this example is to determine the equilibrium configurations of the cable under the assumptions that the pulley is free to move horizontally and that the pulley radius and friction are negligible. For the nonlinear analysis, the cable is subdivided into two segments, one for each span, with the reference curvilinear coordinate <math display="inline">\xi ^1</math> of the pulley as problem unknown. This curvilinear coordinate <math display="inline">\xi ^1</math> is used to construct the finite element mesh in each iteration.

Following the form finding procedure by Argyris ''et al'' <span id='citeF-2'></span><span id='citeF-7'></span>[[#cite-2|[2,7]]], the analysis starts from a straight reference configuration, and imposes a displacement <math display="inline">\mathbf{u}=(-200,0,50)</math> m at the right support and a pair of displacements <math display="inline">u_2 = 50</math> m and <math display="inline">u_3= 100</math> m at the intermediate roller. Because friction is not considered, the jump in the Cauchy axial force at the roller support must be zero. As a result, the problem is solved by iterating over the curvilinear coordinate <math display="inline">\xi ^1</math> so that the jump in the Cauchy axial force at the pulley becomes zero.
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;"
|+ <span id='table-2'></span>Table. 2 <span style="text-align: center; font-size: 75%;">Results for Example 1 from different studies.</span>
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | '''Impollonia et al''' <span id='citeF-6'></span>[[#cite-6|[6]]] 
| style="border-left: 2px solid;border-right: 2px solid;" | '''Present work (continuous)''' 
| style="border-left: 2px solid;border-right: 2px solid;" | '''Present work (discontinuous)''' 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">\xi ^1_1</math> (m) 
| style="border-left: 2px solid;border-right: 2px solid;" | 126.12 
| style="border-left: 2px solid;border-right: 2px solid;" | 126.26 
| style="border-left: 2px solid;border-right: 2px solid;" | 126.25 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">N_1</math> (kN) 
| style="border-left: 2px solid;border-right: 2px solid;" | 14.12 
| style="border-left: 2px solid;border-right: 2px solid;" | 8.31-13.99 
| style="border-left: 2px solid;border-right: 2px solid;" | 8.31-13.99 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">\xi ^1_2</math> (m) 
| style="border-left: 2px solid;border-right: 2px solid;" | 219.98 
| style="border-left: 2px solid;border-right: 2px solid;" | 219.46 
| style="border-left: 2px solid;border-right: 2px solid;" | 219.47 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">N_2</math> (kN) 
| style="border-left: 2px solid;border-right: 2px solid;" | 10.79 
| style="border-left: 2px solid;border-right: 2px solid;" | 4.02-10.68 
| style="border-left: 2px solid;border-right: 2px solid;" | 4.02-10.68 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">\xi ^1_3</math> (m) 
| style="border-left: 2px solid;border-right: 2px solid;" | 424.76 
| style="border-left: 2px solid;border-right: 2px solid;" | 424.70 
| style="border-left: 2px solid;border-right: 2px solid;" | 424.70  
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">N_3</math> (kN) 
| style="border-left: 2px solid;border-right: 2px solid;" | 17.42  
| style="border-left: 2px solid;border-right: 2px solid;" | 10.25-17.28 
| style="border-left: 2px solid;border-right: 2px solid;" | 10.25-17.28 

|}

Table&nbsp;[[#table-2|2]] summarizes the results for the equilibrium configurations with <math display="inline">\xi ^1_i</math> refering to the curvilinear coordinate of the pulley and <math display="inline">N_i</math>, to the axial force. Because the study by Impollonia ''et al'' <span id='citeF-6'></span>[[#cite-6|[6]]] does not consider finite deformations, <math display="inline">\xi ^1_i</math> and <math display="inline">N_i</math> correspond to infinitesimal deformations. For the present study, <math display="inline">\xi ^1_i</math> corresponds to the reference confi-guration and <math display="inline">N_i</math>, to the Cauchy axial force. While the values of the present study agree well with those by Impollonia ''et al'' <span id='citeF-6'></span>[[#cite-6|[6]]], it is worth noting the variation of the Cauchy axial force that the current formulation captures, as indicated by the range of axial force values in Table&nbsp;[[#table-2|2]]. In contrast, the model in <span id='citeF-6'></span>[[#cite-6|[6]]] overestimates the axial force by reporting a value corresponding to the maximum of the current formulation.

Three equilibrium states result from the analysis, as Fig.&nbsp;[[#img-2|2]] shows: three stable confi-gurations denoted with solid lines (<math display="inline">C_1</math> and <math display="inline">C_3</math>), and one unstable configuration, denoted with a dashed line (<math display="inline">C_2</math>), as reflected in the change of direction for the horizontal component of the reaction at the pulley. The <math display="inline">x_1</math> positions of the pulley for these equilibrium states in Fig.&nbsp;[[#img-2|2]] are <math display="inline">x_1^1 = 56.54 / 56.53</math> m, <math display="inline">x_1^2=134.00/134.01</math> m and <math display="inline">x_1^3=274.31/274.31</math> m for the continuous and the discontinuous formulations, respectively.

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Crusells-Girona_842065898-.png|600px|<span style="text-align: center; font-size: 75%;">Deformed shape (30 elements) of equilibrium states for Example 1.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' <span style="text-align: center; font-size: 75%;">Deformed shape (30 elements) of equilibrium states for Example 1.</span>
|}

===4.2 Example 2: Free vibration in finite deformations===

The second example studies the large-amplitude free vibration of two cables with diffe-rent sag/span ratio that were investigated by Srinil ''et al'' <span id='citeF-8'></span>[[#cite-8|[8]]]. The structural model consists of a cable anchored at both ends and spanning  850 m in both cases. Table&nbsp;[[#table-3|3]] summarizes the geometric and material properties of the cables denoted by C1 and C2.
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;"
|+ <span id='table-3'></span>Table. 3 <span style="text-align: center; font-size: 75%;">Geometry and Material Properties for Example 1.</span>
|- style="border-top: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  '''Property''' 
| style="border-left: 2px solid;border-right: 2px solid;" | '''C1''' 
| style="border-left: 2px solid;border-right: 2px solid;" | '''C2''' 
|- style="border-top: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  Cross-sectional area 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.1159 m<math display="inline">^2</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.1159 m<math display="inline">^2</math> 
|-
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Elastic modulus 
| style="border-left: 2px solid;border-right: 2px solid;" | 17.94 GPa 
| style="border-left: 2px solid;border-right: 2px solid;" | 17.94 GPa 
|-
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Density 
| style="border-left: 2px solid;border-right: 2px solid;" | 8337.9 kg/m<math display="inline">^3</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 8337.9 kg/m<math display="inline">^3</math> 
|-
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Cable length 
| style="border-left: 2px solid;border-right: 2px solid;" | 840.48 m 
| style="border-left: 2px solid;border-right: 2px solid;" | 870.51 m 
|- style="border-bottom: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Prestressing 
| style="border-left: 2px solid;border-right: 2px solid;" | - 
| style="border-left: 2px solid;border-right: 2px solid;" | 345 kN 

|}

First, following the shape finding procedure by Argyris ''et al'' <span id='citeF-7'></span>[[#cite-7|[7]]], the equilibrium confi-guration and the first two natural modes of vibration around this configuration are obtained for both cables by solving the standard eigenvalue problem

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\det [\mathbf{K}(\mathbf{u}_{eq})-w^2\mathbf{M}]=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

where <math display="inline">w</math> is the angular frequency, <math display="inline">\mathbf{u}_{eq}</math> refers to the displacement field at the equilibrium state, and <math display="inline">\mathbf{K}</math> and <math display="inline">\mathbf{M}</math> correspond to the static stiffness and mass matrices of the formulation in Sec.&nbsp;[[#3.2 Time integration and consistent linearization|3.2]]. Both cables are discretized with a mesh of 14 elements. Results are presented in Table&nbsp;[[#table-4|4]], where the end tension is given in the Cauchy representation, and "S" and "A" refer to the symmetric and antisymmetric modes, respectively. Fig.&nbsp;[[#img-3|3]] shows the normalized first symmetric and antisymmetric modes of both cables. From this figure, it is interesting to note that, when the sag/span ratio increases, the single extremum for the symmetric mode divides into three because of increasing horizontal displacements.

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Crusells-Girona_842065898-m1.png|546px|]]
|[[Image:draft_Crusells-Girona_842065898-m2.png|400px|<span style="text-align: center; font-size: 75%;">Normalized vertical eigenvectors for Example 2.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 3:''' <span style="text-align: center; font-size: 75%;">Normalized vertical eigenvectors for Example 2.</span>
|}
  
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;"
|+ <span id='table-4'></span>Table. 4 <span style="text-align: center; font-size: 75%;">Results for equilibrium configurations and natural vibration in Example 2.</span>
|- style="border-top: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |   
| colspan='2' style="border-left: 2px solid;border-right: 2px solid;" | '''C1'''
| colspan='2' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | '''C2'''
|-
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | Present 
| style="border-left: 2px solid;border-right: 2px solid;" | Srinil ''et al'' <span id='citeF-8'></span>[[#cite-8|[8]]] 
| style="border-left: 2px solid;border-right: 2px solid;" | Present 
| style="border-left: 2px solid;border-right: 2px solid;" | Srinil ''et al'' <span id='citeF-8'></span>[[#cite-8|[8]]]  
|- style="border-top: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" |  Sag [m] 
| style="border-left: 2px solid;border-right: 2px solid;" | 28.01 
| style="border-left: 2px solid;border-right: 2px solid;" | 28.39 
| style="border-left: 2px solid;border-right: 2px solid;" | 89.28 
| style="border-left: 2px solid;border-right: 2px solid;" | 89.57 
|-
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Sag/span [-] 
| style="border-left: 2px solid;border-right: 2px solid;" | 1/30 
| style="border-left: 2px solid;border-right: 2px solid;" | 1/30 
| style="border-left: 2px solid;border-right: 2px solid;" | 1/9.5 
| style="border-left: 2px solid;border-right: 2px solid;" | 1/9.5 
|-
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | End tension [kN] 
| style="border-left: 2px solid;border-right: 2px solid;" | 30432 
| style="border-left: 2px solid;border-right: 2px solid;" | 30000 
| style="border-left: 2px solid;border-right: 2px solid;" | 10500 
| style="border-left: 2px solid;border-right: 2px solid;" | 10500 
|-
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Frequency (1st S) [Hz] 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.124 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.123 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.158 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.158 
|- style="border-bottom: 2px solid;"
| style="text-align: left;border-left: 2px solid;border-right: 2px solid;" | Frequency (1st A) [Hz] 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.208 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.206 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.112 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.112 

|}

To evaluate the large-amplitude free vibration, an initial displacement field is imposed corresponding to an amplified first symmetric mode, <math display="inline">\mathbf{u}_0=\alpha \mathbf{u}_{m1}</math>, where <math display="inline">\mathbf{u}_{m1}</math> is the normalized first symmetric mode. The parameters for Newmark's method are <math display="inline">\beta=0.25</math>, <math display="inline">\gamma=0.5</math> and <math display="inline">\Delta t=0.05</math> s. No differences are observed between the continuous and the discontinuous formulations as the problem in hand is symmetric in geometry and loads.

Fig.&nbsp;[[#img-4|4]](a) shows the normalized vertical displacements and Cauchy axial forces for cable C1 and <math display="inline">\alpha=15</math>. The evolution of the energy for this case is presented in Fig.&nbsp;[[#img-5|5]](a). Likewise, Fig.&nbsp;[[#img-4|4]](b) and Fig.&nbsp;[[#img-5|5]](b) present the normalized vertical displacements and Cauchy axial forces, and energy evolution, respectively, for cable C2 and <math display="inline">\alpha = 15</math>. While cable C1 behaves linearly in displacements, cable C2 shows a high dependence on high-frequency modes. Also, high-frequency contributions are observed in both cases for the axial force, becoming more relevant for the large sag/span ratio, as observed by Srinil ''et al'' <span id='citeF-6'></span>[[#cite-6|[6]]]. The total energy is conserved for cable C1, whereas it shows minor oscillations for cable C2.

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Crusells-Girona_842065898-plC1.png|540px|]]
|[[Image:draft_Crusells-Girona_842065898-plC2.png|400px|<span style="text-align: center; font-size: 75%;">Normalized vertical displacements and Cauchy axial forces for Example 2.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 4:''' <span style="text-align: center; font-size: 75%;">Normalized vertical displacements and Cauchy axial forces for Example 2.</span>
|}

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Crusells-Girona_842065898-enC1.png|540px|]]
|[[Image:draft_Crusells-Girona_842065898-enC2.png|400px|<span style="text-align: center; font-size: 75%;">Energy evolution for Example 2.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 5:''' <span style="text-align: center; font-size: 75%;">Energy evolution for Example 2.</span>
|}

==5 CONCLUSIONS==

The paper presents a general formulation of catenary elements based on finite deformations and curvilinear coordinates for the nonlinear static and dynamic analysis of cables. From the weak statement of the problem, two implementations are derived: one with a continuous axial force distribution and one with a discontinuous.

As demonstrated by the first example, the formulation is capable of determining equilibrium configurations of three-dimensional cable arrangements with high accuracy, especially in axial forces, compared to other elements in the literature which do not distinguish between Cauchy and 2nd PK axial forces. Furthermore, the second example shows that the natural modes of vibration around equilibrium configurations can also be obtained by the proposed formulation. Because the energy is conserved in the analyzed range of sag/span ratios, Newmark's implicit method can be used to solve the nonlinear dynamic problem. Nevertheless, as observed in the literature, high-frequency contributions in the axial force appear in the analysis, with their amplitude increasing with the sag/span ratio.

In conclusion, because of their consistency and versatility, the proposed catenary elements seem well suited for the nonlinear static and dynamic analysis of nonlinear elastic cables under general loading.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  Andreu, A., Gil, L. and Roca, P. A new deformable catenary element for the analysis of cable net structures. ''Comp. Struct.'' (2006) '''84''':882&#8211;1890. <div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''  Crusells-Girona, M., Filippou, F.C. and Taylor, R.L. A mixed formuation for nonlinear analysis of cable structures. ''Comp. Struct.'' (2017) '''186''':50&#8211;61. <div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''  Newmark, N. A method of computation for structural dynamics. ''J. Eng. Mech.-ASCE'' (1959) '''85''':67&#8211;94. <div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  Taylor, R.L. ''FEAP - Finite Element Analysis Program''. Univ. of California, Berkeley (2014). http://www.ce.berkeley/feap. <div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  Filippou, F. C. ''FEDEASLab - Finite Elements in Design, Evaluation and Analysis of Structures''. Univ. of California, Berkeley (2007). http://www.ce.berkeley.edu/filippou/Courses/FEDEASLab.htm. <div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  Impollonia, N., Ricciardi, G. and Saitta, F. Statics of elastic cables under 3D point forces. ''Int. J. Solids Struct.'' (2011) '''48''':1268&#8211;1276. <div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  Argyris, J. H., Angelopoulos, T. and Bichat, B. A general method for the shape finding of lightweight tension structures. ''Comp. Meth. Appl. Mech. Eng.'' (1974) '''3''':135&#8211;149. <div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  Srinil, N., Rega, G. and Chucheepsakul, S. Three-dimensional non-linear coupling and dynamic tension in large-amplitude free vibrations of arbitrarily sagged cables. ''J. Sound Vib.'' (2004) '''269''':823&#8211;852.