Static analysis of cable structures 1. Introduction Cable structures are lightweight structures of flexible type. The typical cables structure is usually composed of curved members. These members have only axial stiffness related to tension. They can carry neither bending nor compression. Some of them are load bearing ones, some are tightening. They allow spanning large distances. Therefore, they are used in the construction of large span roof girders, shell-like roofs and hanging bridges. 2. Types of cable structures Single cable q N
N
Double cable girders Initially tightened – they reach large stiffness. st iffness. Hence, they exhibit smaller displacements than single cables Load-bearing Load-bearing cable Tensile hanger
Tightening cable
Jahwert’s girder
There are also less frequently used double cable girders with compressive hangers.
Cable meshes
Load-bearing cables
Hyperbolic paraboloid Tightening cables
Hybrid structures – consisting of bars and cables or trusses and cables – suspension bridges single span Load-bearing cable Hanger
Pylon
Beam The main beam can be replaced by a truss multi-span (may have spans of several kilometres)
–
cable-stayed bridges
harp type
fan type
Pylons may have various cross-sections
–
masts with guys Mast
Guy (cable)
3. Basic equation for a single cable Let us consider a single cable subjected to the axial tightening force and transverse loading. q A
B f
N R A
=
ql 2
N C
R B
=
ql 2
l
Since the cable carries no bending we can write down the equation of equilibrium of moments with respect to the centre point C for a half of the cable
∑ M C left
Nf −
ql l ql l ⋅ + ⋅ 2 2 2 4
=
0
and get the basic relation between the tightening force and the uniform transverse loading N =
ql 2 8f
4. Analytical model The model for the cable analysis is based on the following assumptions – the material is linearly elastic – the strains are small but displacements are large – geometric non-linearity – loading is applied at nodes Let us consider a cable element e with the nodes i and k having the initial length L0. x y
e i
z
~ u i i
N ~ w i
~ v i
k
k e ~ y
~ w k
N
~ u k
~ x
~ v k
N – tightening force
L0
~ z The vectors of nodal displacements and nodal reactions in local co-ordinates are given as ~ U ~i V i ~ W ~ R e = ~ i U k ~ V ~k W k
~ u i ~ v i ~ w ~ qe = ~ i u k ~ v ~k w k
The stiffness matrix consists of elastic and geometric parts
~ K e =
1 0 0 0 0 0 EA L0
−1
0 0
0
0
0
0
1
0 0
0 0 0 0 1 0 0 0 0 N 1 0 + 0 L 0 0 0
0
0
−1
0 − 1
0
0
0 0
1
1
or ~ ~ ~ ~ K K ik 1,e K ii 2,e K ik 2,e ~ ~ ~ ii 1, e K e = K1e + K 2e = ~ + ~ ~ ~ K K K K ki 2,e kk 1,e kk 2, e ki 1,e
~ K ii 1,e = ~ − K ii 1,e
~
~ K ii e + ~ 2, − K ii 2,e
− K ii 1, e
~ K ii 1,e
~
− K ii 2, e
~ K ii 2,e
Note, that the geometric matrix includes the current element length L. The equilibrium equation for the element in local co-ordinates reads ~ ~ K e ~ qe = R e In order to analyse the entire system all quantities from local co-ordinates must be transformed to the global co-ordinates
T ~ K e = Te K e Te
where the transformation matrix is
C 0 Te = 0 C
T
with
c xx c xy c xz C = c yx c yy c yz c zx c zy c zz
The direction cosine c ab is defined as ~ c ab = cos ∠ a , b
( )
The first part of the element stiffness matrix can be transformed as follows ~ CT 0 K ii e K1e = ~ 1, T 0 C − K ii 1,e ~ CT K ii 1,e C = T ~ − C K ii 1,e C
~
~ CT K ii 1,e = T ~ 0 C − C K ii 1,e
− K ii 1,e C
~ K ii 1,e
0
T ~
K ii 1,e C 0 ~ = CT K ii 1,e 0 C
−C
T ~
K ii 1,e C ~ CT K ii 1,e C
−C
and similarly for the second part we get ~ CT K C K 2e = T ~ii 2,e − C K ii 2,e C
T ~
K ii 2,e C ~ CT K ii 2,e C
−C
We can substitute the explicit form of the matrices and get
c xx c xy c xz 1 0 0 c xx c yx c zx c xx 0 0 c xx c yx c zx EA EA ~ 0 0 0 c c c = c 0 0 c c c CT K ii 1,e C = c c c yx yy yz xy yy zy L yx xy yy zy L0 0 c zx c zy c zz 0 0 0 c xz c yz c zz c zx 0 0 c xz c yz c zz =
c xx 2 c xx c yx c xx c zx EA 2 c c c c c xx yx yx zx yx L0 2 c c c c c zx xx zx zx yx
~ CT K ii 2,e C =
=
=
c xx c xy c xz 0 0 0 c xx c yx c zx N c c c c c c 0 1 0 yx yy yz xy yy zy L c zx c zy c zz 0 0 1 c xz c yz c zz
(c xy 2 + c xz 2 ) N (c xy c yy + c xz c yz ) L (c xy c zy + c xz c zz )
=
0 c xy c xz c xx c yx c zx N 0 c c c c c yy yz xy yy zy L 0 c zy c zz c xz c yz c zz
(c xy c yy + c xz c yz ) (c xy c zy + c xz c zz ) (c yy 2 + c yz 2 ) (c yy c zy + c yz c zz ) (c yy c zy + c yz c zz ) (c zy 2 + c zz 2 )
The elements of the second submatrix can be rewritten if the orthogonality conditions for the direction cosines are used. We have for instance: 2
+ c xy + c xz = 1
2
+ c xz = 1 − c xx
c xx
2
2
so the first diagonal term can be given as c xy There is also
2
2
=
c xx c yx + c xy c yy + c xz c yz = 0 so the first off-diagonal term can be given as c xy c yy + c xz c yz = −c xx c yx Using the similar relations the entire second part of the stiffness matrix can be rewritten and finally we get
K ii 1,e
− K ii 1, e
− K ii 1,e
K ii 1,e
K e =
+
K ii 2,e − K ii 2,e
− K ii 2, e
K ii 2,e
where:
K ii 1,e =
K ii 2,e =
c xx 2 c xx c yx c xx c zx EA 2 c zx c yx c xx c yx c yx L0 2 c xx c zx c zx c yx c zx
(1 − c xx 2 ) − c xx c yx − c xx c zx N 2 − c xx c yx (1 − c yx ) − c yx c zx L 2 − c xx c zx − c yx c zx (1 − c zx )
~ Thus we have expressed the stiffness matrix in terms of direction cosines between the local axis x and the global axes x , y and z . Because of such a relatively simple form, usually the element matrices are directly given in the global co-ordinates using the co-ordinates of nodes. Let us consider the element
x x i
y i z i
x k y
y k
i L
e
z k k
z
~ x
The cosines required to express the element stiffness matrix can be given as
~) = x k − x i c xx = cos ∠(x , x L ~) = y k − y i c yx = cos ∠(y , x L ~) = z k − z i c zx = cos ∠(z , x L Thus, with these definitions there is no necessity for the transformation of co-ordinates. After the assembly of the global stiffness matrix K and inclusion of support conditions the global set of equations is formulated Kq = P The global stiffness matrix has the dimension n ×n , where n = 3s and s is the number of free (unsupported) nodes.
Solution of this equation gives the nodal displacements in the vector q. Then the equilibrium conditions in the deformed configuration must be checked. The new co-ordinates of nodes are obtained from x i ' = x i + u i y i ' = y i + v i z i ' = z i + w i x k ' = x k + u k y k ' = y k + v k z k ' = z k + w k
and the displacements u i, v i and w i come from q. The current length of the element is: L' =
(x k '− x i ')2 + (y k '−y i ')2 + (z k '−z i ')2
so the length increment is ∆L' =
(x k '− x i ')2 + (y k '−y i ')2 + (z k '−z i ')2 − L0
The current value of the axial force in the element is N ' = N + ∆N ' where the current force increment is obtained from the Hooke’s law ∆L ' Having found the current axial forces we can calculate the out-of-balance nodal ∆N ' = EA L0 forces ∆Q x', ∆Q y' and ∆Q z' assembled into a vector
Q x'
Q y'
∆Q x ' ∆Q' = ∆Q y ' ∆Q z '
x m
Q z'
N r'
y
N r ' = −∑ N e ' c e =1
z Let us denote the resultant from the axial forces in the m -elements coinciding at a node as N r'. Then, using the equilibrium conditions at the node we get m
~) = 0 Q x '+ ∑ N e ' cos ∠(x , x e = 1 m
~) = 0 Q y '+ ∑ N e ' cos ∠(y , x e = 1 m
~) = 0 Q z '+ ∑ N e ' cos ∠(z , x e =1
or in the matrix notation m
Q' = − ∑ N e ' c e =1
where
Q x ' Q' = Q y ' Q z '
~) c cos ∠(x , x xx ~ c = cos ∠(y , x ) = c yx ~) c cos ∠(z , x zx
The external forces at the nodes are P' = P0
+P
where P0 is the vector of initial (tightening) nodal forces and P is the vector of external imposed loads. The out-of-balance forces are m
∆Q' = P'−Q' = P0 + P +
∑ N e ' c e =1
In the next iteration the equilibrium equation reads
K' ∆q' = ∆Q' Solving of these equations yields the vector of displacement increments, which are used to calculate the new values of current displacements
q' = q + ∆q' From this point the calculations are repeated to get the new values of out-of-balance forces, etc. The iterations stop when the current value of the out-of-balance forces falls below a tolerance ∆Q' ≤
tol
The presented algorithm follows the line of the Newton method. Let us now briefly present the methods of iterative solving of non-linear equations. The Newton method
P - a new stiffness matrix at each iteration ∆Q '
P
Q ' q q
∆q '
∆q ''
q real The modified Newton method
P
- the initial stiffness matrix used at each iteration - slower convergence but low computation cost at each iteration
∆Q '
P Q ' q q
∆q '
q real
The Newton-Raphson method P - loading divided into increments - much better convergence in highly nonlinear cases
P 3 P P 2 P 1
q
5. Example We consider a plane cable system consisting of 3 elements. Initial configuration A
B
x
H
1 C y
H
3
R A
Node co-ordinates: x A = 0.0 y A = 0.0 x B = 10.0 y B = 0.0 x C = 3.0 y C = 3.0 x D = 7.0 y D = 3.0
2
3
R B
D
P 0
P 0
3
4
3
Data: EA = 8000 kN Initial tightening load – two nodal forces P 0 = 5 kN The reactions in the initial configuration are
∑ M B :
10R A
−5⋅7 − 5⋅3 =
∑ y : ∑ M C :
3R A
R B
− 3H =
=
0 ⇒ R A
=
5 kN
5 kN
0 ⇒ H = 5 kN
The axial forces in the inclined cables can be obtained from the equilibrium of the node A or B A H
N 1
∑ y :
N 1
R A
2 2
=
5 kN ⇒ N 1 = 7.07107 kN
=
N 2 ⇒ N 2
And from the symmetry N 3 = N 1. The force in the cable CD comes from the equilibrium of the node C N 1
C N 2 P 0
∑ x :
N 1
2 2
=
5 kN
The cable system in this configuration is now loaded by a set of two external loads P = 12 kN A
B
x
H
1
H
3
R A
C y
D
2
P 0
P 0
P
P
3
3
R B
4
3
Initial lengths for the elements are
(3 − 0 )2 + (3 − 0)2
L0,1 =
L0,2
=
=
4.24264m = L0,3
(7 − 3 )2 + (3 − 3)2
=
4.0m
The direction cosines are: c xx ,1 =
x C − x A L0,1
=
3−0 4.24264
=
0.707107
c yx ,1 =
y C − y A L0,1
=
3−0 4.24264
=
0.707107
c xx ,2
=
x D − x C L0,2
=
7−3 4.0
= 1.0
c yx ,2
=
y D − y C L0,2
=
3−3 4.0
=
0.0
c xx ,3
=
x B − x D L0,3
=
10 − 7 4.24264
=
c yx ,3
=
y B − y D L0,3
=
0−3 4.24264
= −0.707107
0.707107
The element stiffness matrices
K AA,1
− K AA,1
− K AA,1
K AA,1
K1 =
(
=
)
c xx ,1c yx ,1 N 1 1 − c xx ,1 EA c xx ,1 + 2 L0,1 c xx ,1c yx ,1 c yy ,1 L0,1 − c xx ,1c yx ,1 2
K AA,1 =
8000 4.24264
0.5 0.5 7.07107 0.5 0.5 0.5 + 4.24264 − 0.5 K2
=
K AA,2 − K AA,2
2
− 0.5
0.5
− K AA,2
K AA,2
=
− c xx ,1c yx ,1
(1 − c ) 2 yy ,1
=
943.643 941.976 941.976 943.643
c xx ,2 2 c xx ,2c yx ,2 N 2 K AA,2 = + 2 c c c xx ,2 yx ,2 L0,2 yy ,2 8000 1 0 5 0 0 2000 .0 = + = 4.0 0 0 4.0 0 1 0 EA L0,2
K3
=
K AA,3 − K AA,3
(1 − c xx ,2 2 ) − c xx ,2c yx ,2 0 1.250
− c xx ,2 c yx ,2
(1 − c ) 2
=
yy ,2
− K AA,3
K AA,3
2 c xx ,3 2 c xx ,3c yx ,3 N 3 (1 − c xx ,3 ) − c xx ,3c yx ,3 K AA,3 = + = 2 2 c yy ,3 L0,3 − c xx ,3 c yx ,3 (1 − c yy ,3 ) c xx ,3c yx ,3 8000 0.5 − 0.5 7.07107 0.5 0.5 943.643 − 941.976 = + = 4.24264 − 0.5 0.5 4.24264 0.5 0.5 − 941.976 943.643
EA L0,3
The pattern of global stiffness matrix assembly and support conditions inclusion is A A K=
C
C
D
B
K1 K2 K3
D B
and the remaining fragment of the matrix K corresponds to the free nodes C and D. This remaining fragment is
2943 .64 941.976 941.976 944.893 K= − 2000 .0 0. 0 − 1.250 0.0
0.0 − 1.250 2943 .64 − 941.976 944.893 − 941.976 − 2000.0
0.0
The load vector corresponding to the free nodes C and D has two forces P in the y direction and zero forces in the x direction
0.0 12.0 P= 0.0 12.0 The solution of the equilibrium equation Kq = P is
− 0.0029922 0.0157036 q= 0.0029922 0.0157036 Now we can find the new co-ordinates of free nodes
x C ' = 3.0 − 0.0029922 = 2.99701 y C ' = 3.0 + 0.0157036 = 3.01570 x D ' = 7.0 + 0.0029922 = 7.00299 y D ' = 3.0 + 0.0157036
=
3.01570
The current element lengths are L1 ' = L2 ' =
(2.99701 − 0.0)2 + (3.01570 − 0.0)2
=
4.25165m = L3 '
(7.00299 − 2.99701)2 + (3.01570 − 3.01570 )2
=
4.00598m
The length increments are ∆L0,1 ' =
4.25165 − 4.24264
∆L0,2 ' =
=
0.00901m = ∆L0,3'
4.00598 − 4.0 = 0.00598m
The increments of axial forces ∆N 1' =
EA
∆N 2 ' =
∆L0,1'
=
L0,1 EA
8000 ⋅
∆L0,2 '
=
L0,2
0.00901 4.24264
8000 ⋅
= 16.9894kN = ∆N 3 '
0.00598 4. 0
= 11.9680 kN
and the current values of the axial forces N 1' = N 1 + ∆N 1' = 7.07107 + 16.9894 = 24.0605kN = N 3 ' N 2 ' = N 2
+ ∆N 2 ' =
5.0 + 11.9680 = 16.9680kN
The current values of the direction cosines are
c xx ,1 = c yx ,1 =
c xx ,3 c yx ,3
2.99701 − 0 = 0.704905 4.25165 3.01570 − 0 = 0.709301 4.25165 c xx ,2
= 1.0
c yx ,2
=
0.0
10 − 7.00299 = 0.704905 4.25165 0 − 3.01570 = = −0.709301 4.25165 =
The out-of-balance forces can be found from the equilibrium of nodes
24.0605
C
x 16.9680
5 + 12 y
∆Q x ' = −24.0605 ⋅ 0.704905 + 16.9680 =
0.00763kN
∆Q y ' = −24.0605 ⋅ 0.709301 + 12 + 5 = −0.06614kN
And similarly, from the equilibrium of the node D we get ∆Q x ' = −0.00763 kN ∆Q y ' = −0.06614kN
to complete the vector
0.00763 − 0.06614 ∆Q' = − 0.00763 − 0.06614 The stiffness matrices in the current configuration are formed using the following submatrices
(
2
K AA,1' = =
)
c xx ,1c yx ,1 N 1' 1 − c xx ,1 EA c xx ,1 + 2 L0,1 c xx ,1c yx ,1 c yy ,1 L1' − c xx ,1c yx ,1 2
− c xx ,1c yx ,1
(1 − c ) 2 yy ,1
0.50 24.0605 0.503108 8000 0.496891 + 4.24264 0.50 0.503108 4.25165 − 0.50 K AA,2 ' = =
EA L0,2
=
939.794 939.980 = 0.496891 939.980 951.482 − 0.50
2 c xx ,2 2 c xx ,2c yx ,2 N 2 ' (1 − c xx ,2 ) + 2 ' L c c c − c xx ,2c yx ,2 2 yy ,2 xx ,2 yx ,2
8000 1 0 16.9680 0 0 + 4.0 0 0 4.00598 0 1
=
− c xx ,2 c yx ,2
(1 − c ) 2
=
yy ,2
0 2000.0 0 4.23567
2 c xx ,3 2 c xx ,3 c yx ,3 N 3 ' (1 − c xx ,3 ) − c xx ,3 c yx ,3 + = K AA,3 ' = 2 2 L ' ( ) − − c c c c c 1 c xx ,3 yx ,3 3 yy ,3 yy ,3 xx ,3 yx ,3 0.50 − 0.50 8000 0.496891 24.0605 0.503108 = + 4.24264 − 0.50 0.503108 4.25165 0.50 0.496891
EA L0,3
=
939.794 − 939.980
− 939.980
951.482
Following the same assembly and support conditions pattern we get the “active” part of the global stiffness matrix in the form
2939.79 939.980 K' = − 2000 .0 0.0
955.718 0. 0 − 4.23567 0. 0 2939 .79 − 939.980 − 4.23567 − 939.980 955.718 939.980
− 2000.0
and solving the equilibrium equations with out-of-balance forces K' ∆q' = ∆Q' we get the following increments of displacements
0.0000182 − 0.0000875 ∆q' = − 0.0000182 − 0.0000875 With these values the total current displacements can be found
0.0
u C ' − 0.0029922 + 0.0000182 v ' C = ∆q'+q = 0.0157037 − 0.0000852 u D ' 0.0029922 − 0.0000182 0.0157037 − 0.0000852 v D '
− 0.0029738 0.0156162 = 0.0029738 0.0156162
The co-ordinates of nodes in the current deformed configuration x C ' ' = 3.0 − 0.0029738
=
2.997026
y C ' ' = 3.0 + 0.0156162 = 3.016162 x D ' ' = 7.0 + 0.0029738
=
7.002974
y D ' ' = 3.0 + 0.0156162 = 3.016162 The current element lengths are L1' ' = 4.251600 m = L3 ' ' L2 ' ' = 4.005948m The length increments are ∆L0,1 ' ' =
0.008960 m = ∆L0,3 ' '
∆L0,2 ' ' =
0.005948 m
The increments of axial forces ∆N 1 ' ' = 16.8951kN = ∆N 3 ' ' ∆N 2 ' ' = 11.8960 kN
and the current values of the axial forces N 1' ' = N 1 + ∆N 1' ' = 23.9662kN = N 3 ' ' N 2 ' ' = N 2
+ ∆N 2 ' ' = 16.8960kN
The current values of the direction cosines are c xx ,1
=
0.704917
c yx ,1
=
0.709290
c xx ,2
= 1.0
c yx ,2
=
0.0
0.704917
c xx ,3
=
c yx ,3
= −0.709290
The current out-of-balance forces can be found from the equilibrium of nodes 23.9662
C
x 16.8960
5 + 12 y
∆Q x ' ' = −23.9662 ⋅ 0.704917 + 16.8960 = ∆Q y ' ' = −23.9662 ⋅ 0.709290 + 12 + 5 =
0.00182kN
0.00101kN
These values are several times smaller than the out-of-balance forces in the first iteration, what confirms the convergence behaviour of the solution process. They are also relatively small compared to the current values of forces in the elements (about 0.001%), so the iterations can be stopped at this point.
