Nonlinear static and dynamic analysis of mixed cable elements

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.

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 [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 [2] that overcomes these limitations.

2 MIXED FORMULATION OF THE CATENARY PROBLEM

2.1 Kinematics

Fig. 1 shows the motion of a cable from a reference configuration ${\textstyle {\mathcal {P}}_{0}}$ to a current configuration ${\textstyle {\mathcal {P}}}$. Define an orthogonal frame ${\textstyle \{\mathbf {G} _{i}\}_{i=1}^{3}}$ with associated coordinates ${\textstyle \{\xi ^{i}\}_{i=1}^{3}}$ at any point ${\textstyle P\in {\mathcal {P}}_{0}}$, such that

$\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}\|}}$

(1)

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

Figure 1: Motion $\mathbf {x} ={\boldsymbol {\chi }}(\mathbf {X} (\xi ^{1}))$ of the cable ${\mathcal {C}}$.

The relevant stretch and Green-Lagrange strain of the problem, in the ${\textstyle \mathbf {g} _{1}}$ direction, are

$\lambda ={\frac {\|\mathbf {g} _{1}\|}{\|\mathbf {G} _{1}\|}}\;\;\;\;\;\;;\;\;\;\;\;\;\mathrm {E} ={\frac {1}{2}}(\lambda ^{2}-1)\|\mathbf {G} _{1}\|^{2}$

(2)

The displacement vector ${\textstyle \mathbf {u} }$ depends only on the curvilinear coordinate ${\textstyle \xi ^{1}}$,

$\mathbf {u} (\mathbf {X} (\xi ^{1}))=\mathbf {x} (\mathbf {X} (\xi ^{1}))-\mathbf {X} (\xi ^{1})=u_{A}(\xi ^{1})\mathbf {E} _{A}$

(3)

Therefore, the relationship between the displacement field ${\textstyle \mathbf {u} }$ and the relevant Green-Lagrange strain ${\textstyle \mathrm {E} }$ can be computed [2] as

$\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)$

(4)

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

$\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}$

(5)

Hence, the components

$\mathrm {\hat {E}} _{AB}=\mathrm {E} _{ij}(\mathbf {\hat {G}} _{A}\cdot \mathbf {G} ^{i})(\mathbf {\hat {G}} _{B}\cdot \mathbf {G} ^{j})$

(6)

are physical quantities.

2.2 Equilibrium and principle of virtual work

For expressing the equilibrium equation of the cable in finite deformations, let ${\textstyle \mathbf {n} }$ 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, ${\textstyle \mathbf {N} }$. It can be shown [2] that, using the orthonormal basis in Eq. 5, namely with components ${\textstyle \mathbf {n} =\mathrm {\hat {n}} \mathbf {\hat {g}} _{1}}$ and ${\textstyle \mathbf {N} =\mathrm {\hat {N}} \mathbf {\hat {G}} _{1}}$,

$\mathrm {\hat {n}} =\lambda \mathrm {\hat {N}}$

(7)

Denoting by ${\textstyle s}$ and ${\textstyle S}$ the arc-length coordinates in the current and reference configurations, respectively, the cable distributed load can be described as

$\mathbf {w} ds=\mathbf {\hat {W}} dS=\mathbf {W} d\xi ^{1}$

(8)

Then, global equilibrium for the cable in the current configuration ${\textstyle {\mathcal {P}}}$ states

$\mathbf {n} (s)-\mathbf {n} (0)+\int _{0}^{s}\mathbf {w} \,ds=\int _{0}^{s}\rho \mathbf {a} \,ds$

(9)

where ${\textstyle \rho }$ is the material density in the current configuration and ${\textstyle \mathbf {a} }$, the total acceleration. The corresponding local statement in the current configuration can be obtained with the fundamental theorem of calculus and the localization theorem

${\frac {d}{ds}}(\mathrm {n} \mathbf {g} _{1})+\mathbf {w} =\rho \mathbf {a}$

(10)

or, in the reference configuration,

${\frac {d}{dS}}(\mathrm {\hat {N}} {\sqrt {G^{11}}}\mathbf {g} _{1})+\mathbf {\hat {W}} =\rho _{0}\mathbf {a}$

(11)

where ${\textstyle G^{ij}=\mathbf {G} ^{i}\cdot \mathbf {G} ^{j}}$ represents the dual metric tensor and ${\textstyle \rho _{0}=\lambda \rho }$.

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

$\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.$

(12)

for ${\textstyle 0<S<L}$ equivalent to ${\textstyle \xi _{1}^{1}<\xi ^{1}<\xi _{2}^{1}}$.

The corresponding two-field weak statement of Eq. 12 can be obtained by considering any variation ${\textstyle \delta \mathbf {u} \in {\mathcal {V}}}$, the space of displacement test functions, any variation ${\textstyle \delta \mathrm {\hat {N}} \in {\mathcal {W}}}$, the space of axial force test functions, and integrating the equilibrium equation by parts,

$\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.$

(13)

where the constitutive relation is imposed strongly. The spaces for the trial solutions of the displacements and axial forces, ${\textstyle {\mathcal {S}}}$ and ${\textstyle {\mathcal {N}}}$, respectively, are

${\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}}$

(15)

Similarly, the spaces for the test functions of the displacements and the axial forces, ${\textstyle {\mathcal {V}}}$ and ${\textstyle {\mathcal {W}}}$, respectively, become

${\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}}$

(17)

where ${\textstyle H^{k}(\Omega )}$ is the Sobolev space for the ${\textstyle k-}$th weak derivative in the ${\textstyle L^{2}(\Omega )}$ norm.

As a result, there are no continuity requirements for the axial force field ${\textstyle \mathrm {\hat {N}} }$. 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 ${\textstyle \mathrm {\hat {N}} (\xi ^{1})}$, ${\textstyle \mathbf {u} (\xi ^{1})}$ and ${\textstyle \mathbf {a} (\xi ^{1})}$. Assume a ${\textstyle k}$-th order Galerkin interpolation for the axial forces

$\mathrm {\hat {N}} ={\boldsymbol {\varphi }}^{t}\mathbf {\hat {N}} =\mathbf {\hat {N}} ^{t}{\boldsymbol {\varphi }}$

(18)

and an ${\textstyle l}$-th order Galerkin interpolation for the displacement and acceleration fields

$\mathbf {u} ={\boldsymbol {\phi }}^{t}\mathbf {\hat {u}} \;\;\;\;\;\;;\;\;\;\;\;\;\mathbf {a} ={\boldsymbol {\phi }}^{t}\mathbf {\hat {a}}$

(19)

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

$\mathbf {x} =\mathbf {X} +\mathbf {u} ={\boldsymbol {\phi }}^{t}(\mathbf {\hat {X}} +\mathbf {\hat {u}} )={\boldsymbol {\phi }}^{t}\mathbf {\hat {x}}$

(20)

With these interpolation functions, one can discretize the weak statement in Eq. 13 for a finite element ${\textstyle \Omega _{e}}$ and a time step ${\textstyle n}$ as

$\left\{{\begin{array}{r}\displaystyle \int _{\Omega _{e}}\delta \mathbf {\hat {N}} ^{t}{\boldsymbol {\varphi }}\left\{G^{11}\mathbf {\hat {u}} _{n}^{t}{\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.$

(21)

where ${\textstyle (\cdot )'}$ represents the derivative with respect to the curvilinear coordinate ${\textstyle \xi ^{1}}$ and ${\textstyle \mathbf {\hat {g}} _{n}=\mathbf {G} _{1}+({\boldsymbol {\phi }}')^{t}\mathbf {\hat {u}} _{n}}$ is the numerical counterpart to ${\textstyle \mathbf {g} _{1}}$.

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 ${\textstyle e}$ elements of ${\textstyle k}$-th order in axial forces and ${\textstyle l}$-th order in displacements. By defining the expanded stress divergence term ${\textstyle \mathbf {R} =(\mathbf {R} _{1},\mathbf {R} _{2})}$ with components

${\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}} _{n}^{t}{\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}}$

(23)

and the mass matrix ${\textstyle \mathbf {M} }$ as

$\mathbf {M} =\int _{\Omega _{e}}\rho _{0}{\boldsymbol {\phi }}{\boldsymbol {\phi }}^{t}\,dS$

(24)

one can rewrite Eq. 21 in an implicit scheme as

$\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]$

(25)

where ${\textstyle \mathbf {F} _{ext,n+1}}$ refers to the external forces considered at the time step ${\textstyle n+1}$.

Introducing Newmark's time integrator [3], one obtains the system of equations

${\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}}$

(32)

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

${\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}$

(33)

where the Fréchet derivative ${\textstyle \partial {\boldsymbol {\Phi }}/\partial \mathbf {\hat {V}} _{n+1}|_{\mathbf {\hat {V}} _{n+1}^{(k)}}}$ corresponds to the dynamic stiffness ${\textstyle \mathbf {K} }$ of the problem, with components

${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$ $\mathbf {K} _{\mathbf {N} \mathbf {u} }=\int _{\Omega _{e}}G^{11}{\boldsymbol {\varphi }}\mathbf {\hat {g}} _{n+1}^{t}({\boldsymbol {\phi '}})^{t}dS=\mathbf {K} _{\mathbf {u} \mathbf {N} }^{t}$ $\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$ $\mathbf {K} _{\mathbf {u} \mathbf {u} }^{d}={\frac {1}{\beta \Delta t_{n}^{2}}}\mathbf {M} +\mathbf {K} _{\mathbf {u} \mathbf {u} }^{s}$

(34)

in the form

$\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]$

(35)

In order to satisfy stability of the solution scheme, it is necessary [2] that

$\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}$

(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 ${\textstyle e}$ elements of ${\textstyle k}$-th order in axial forces and ${\textstyle l}$-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 ${\textstyle \mathrm {\hat {N}} \in H^{0}(0,L)}$. The stress divergence term ${\textstyle \mathbf {R} (\mathbf {\hat {N}} _{n},\mathbf {\hat {u}} _{n})}$ is then understood as

$\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$

(37)

and one can rewrite Eq. 21 in an implicit scheme as

$\mathbf {R} (\mathbf {\hat {u}} _{n+1})+\mathbf {M} \mathbf {\hat {a}} _{n+1}=\mathbf {F} _{ext,n+1}$

(38)

Introducing Newmark's time integrator [3], one obtains the system of equations

${\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}$

(39)

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

${\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}$

(40)

where the Fréchet derivative ${\textstyle \partial {\boldsymbol {\Phi }}/\partial \mathbf {\hat {u}} _{n+1}|_{\mathbf {\bar {u}} _{n+1}^{(k)}}}$ corresponds to the condensed dynamic stiffness ${\textstyle \mathbf {K} }$ of the problem

$\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} }$

(41)

with the components defined in Eq. 34.

The stability condition of the solution scheme reads [2] in this case as

$\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}$

(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 [4] and Matlab toolbox FEDEASLab [5], use a linear approximation for the axial forces (${\textstyle k=1}$) and a quadratic approximation for the displacements ${\textstyle (l=2)}$. 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 ${\textstyle {\mathcal {U}}}$ in terms of the stretch ${\textstyle \lambda }$ and the generalized Young's modulus ${\textstyle E}$

${\mathcal {U}}={\frac {E}{8}}(\lambda ^{2}-1)^{2}$

(43)

Thus, if ${\textstyle A}$ is the area of the cross section and ${\textstyle \mathrm {\hat {N}} _{0}}$ the prestressing force,

$\mathrm {\hat {N}} -\mathrm {\hat {N}} _{0}=(EA)\mathrm {\hat {E}}$

(44)

with constant material stiffness ${\textstyle \partial \mathrm {\hat {N}} /\partial \mathrm {\hat {E}} =EA}$.

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 [6]. The structural model, whose geometric and material properties are shown in Table 1, consists of a cable anchored at both ends and supported by an intermediate roller. In this case, the inertia forces in Eq. 21 and the mass term of the stiffness are not considered as the problem is analyzed in a static manner.

Table. 1 Geometry and Material Properties for Example 1.

Property	Value
Cross-sectional area	805 mm${\textstyle ^{2}}$
Elastic modulus	16.0 kN/mm${\textstyle ^{2}}$
Cable self-weight	62.0679 N/m
Cable length	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 ${\textstyle \xi ^{1}}$ of the pulley as problem unknown. This curvilinear coordinate ${\textstyle \xi ^{1}}$ is used to construct the finite element mesh in each iteration.

Following the form finding procedure by Argyris et al [2,7], the analysis starts from a straight reference configuration, and imposes a displacement ${\textstyle \mathbf {u} =(-200,0,50)}$ m at the right support and a pair of displacements ${\textstyle u_{2}=50}$ m and ${\textstyle u_{3}=100}$ 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 ${\textstyle \xi ^{1}}$ so that the jump in the Cauchy axial force at the pulley becomes zero.

Table. 2 Results for Example 1 from different studies.

	Impollonia et al [6]	Present work (continuous)	Present work (discontinuous)
${\textstyle \xi _{1}^{1}}$ (m)	126.12	126.26	126.25
${\textstyle N_{1}}$ (kN)	14.12	8.31-13.99	8.31-13.99
${\textstyle \xi _{2}^{1}}$ (m)	219.98	219.46	219.47
${\textstyle N_{2}}$ (kN)	10.79	4.02-10.68	4.02-10.68
${\textstyle \xi _{3}^{1}}$ (m)	424.76	424.70	424.70
${\textstyle N_{3}}$ (kN)	17.42	10.25-17.28	10.25-17.28

Table 2 summarizes the results for the equilibrium configurations with ${\textstyle \xi _{i}^{1}}$ refering to the curvilinear coordinate of the pulley and ${\textstyle N_{i}}$, to the axial force. Because the study by Impollonia et al [6] does not consider finite deformations, ${\textstyle \xi _{i}^{1}}$ and ${\textstyle N_{i}}$ correspond to infinitesimal deformations. For the present study, ${\textstyle \xi _{i}^{1}}$ corresponds to the reference confi-guration and ${\textstyle N_{i}}$, to the Cauchy axial force. While the values of the present study agree well with those by Impollonia et al [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 2. In contrast, the model in [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. 2 shows: three stable confi-gurations denoted with solid lines (${\textstyle C_{1}}$ and ${\textstyle C_{3}}$), and one unstable configuration, denoted with a dashed line (${\textstyle C_{2}}$), as reflected in the change of direction for the horizontal component of the reaction at the pulley. The ${\textstyle x_{1}}$ positions of the pulley for these equilibrium states in Fig. 2 are ${\textstyle x_{1}^{1}=56.54/56.53}$ m, ${\textstyle x_{1}^{2}=134.00/134.01}$ m and ${\textstyle x_{1}^{3}=274.31/274.31}$ m for the continuous and the discontinuous formulations, respectively.

Figure 2: Deformed shape (30 elements) of equilibrium states for Example 1.

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 [8]. The structural model consists of a cable anchored at both ends and spanning 850 m in both cases. Table 3 summarizes the geometric and material properties of the cables denoted by C1 and C2.

Table. 3 Geometry and Material Properties for Example 1.

Property	C1	C2
Cross-sectional area	0.1159 m${\textstyle ^{2}}$	0.1159 m${\textstyle ^{2}}$
Elastic modulus	17.94 GPa	17.94 GPa
Density	8337.9 kg/m${\textstyle ^{3}}$	8337.9 kg/m${\textstyle ^{3}}$
Cable length	840.48 m	870.51 m
Prestressing	-	345 kN

First, following the shape finding procedure by Argyris et al [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

$\det[\mathbf {K} (\mathbf {u} _{eq})-w^{2}\mathbf {M} ]=0$

(45)

where ${\textstyle w}$ is the angular frequency, ${\textstyle \mathbf {u} _{eq}}$ refers to the displacement field at the equilibrium state, and ${\textstyle \mathbf {K} }$ and ${\textstyle \mathbf {M} }$ correspond to the static stiffness and mass matrices of the formulation in Sec. 3.2. Both cables are discretized with a mesh of 14 elements. Results are presented in Table 4, where the end tension is given in the Cauchy representation, and "S" and "A" refer to the symmetric and antisymmetric modes, respectively. Fig. 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.

Figure 3: Normalized vertical eigenvectors for Example 2.

Table. 4 Results for equilibrium configurations and natural vibration in Example 2.

	C1		C2
	Present	Srinil et al [8]	Present	Srinil et al [8]
Sag [m]	28.01	28.39	89.28	89.57
Sag/span [-]	1/30	1/30	1/9.5	1/9.5
End tension [kN]	30432	30000	10500	10500
Frequency (1st S) [Hz]	0.124	0.123	0.158	0.158
Frequency (1st A) [Hz]	0.208	0.206	0.112	0.112

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

Fig. 4(a) shows the normalized vertical displacements and Cauchy axial forces for cable C1 and ${\textstyle \alpha =15}$. The evolution of the energy for this case is presented in Fig. 5(a). Likewise, Fig. 4(b) and Fig. 5(b) present the normalized vertical displacements and Cauchy axial forces, and energy evolution, respectively, for cable C2 and ${\textstyle \alpha =15}$. 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 [6]. The total energy is conserved for cable C1, whereas it shows minor oscillations for cable C2.

Figure 4: Normalized vertical displacements and Cauchy axial forces for Example 2.

Figure 5: Energy evolution for Example 2.

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

[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–1890.

[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–61.

[3] Newmark, N. A method of computation for structural dynamics. J. Eng. Mech.-ASCE (1959) 85:67–94.

[4] Taylor, R.L. FEAP - Finite Element Analysis Program. Univ. of California, Berkeley (2014). http://www.ce.berkeley/feap.

[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.

[6] Impollonia, N., Ricciardi, G. and Saitta, F. Statics of elastic cables under 3D point forces. Int. J. Solids Struct. (2011) 48:1268–1276.

[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–149.

[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–852.