CIVE50003 Computational Methods II

Coursework – Dynamic analysis of a cable-stayed bridge

This project is to be carried out individually using the Matlab programming
environment. Please make an electronic submission on Blackboard of a report (no
more than 12 pages in pdf format) and your complete Matlab code (all m files). Your
report should include a detailed critique of your findings in the context of finite
element and structural theory, including computational considerations where
appropriate. Please write carefully and professionally and use appropriate formatting
in your document. You may collaborate and use any code released to you, but you
may not share any files or results. The deadline is 10 am on Friday 25th March 2022.

You are asked to investigate the dynamic response of a cable-stayed bridge shown in
Fig. 1. The structure shown is an idealisation of the main portion of the Theodor
Heuss Bridge in Düsseldorf. It is famous for being the first cable-stayed bridge in
Germany, built between 1953 and 1957.

Assume that L = 108 m and H = 40 m. It is proposed to model the bridge deck
ABMCD with 2D general beam elements d1 to d14, for which you can assume A =
1.88 m2 and I = 1.8036 m4. The pylons BE and CF can be modelled with 2D general
beam elements p1 to p6, for which you can assume A = 0.49 m2 and I = 0.4904 m4.
The pin-ended cables can be modelled with 2D bar elements c1 to c6 for which you
can assume A = 0.0531 m2 (why do you not need I?). There should be 32 elements in
total. These properties assume only one plane through the bridge. Assume a modulus
of elasticity of steel of 200 GPa and a density of ρ = 7850 kg/m3. Assume that the
bridge is already prestressed to support its own self-weight and any other ‘permanent’
loads, which you do not need to consider. However, in Q1 & Q3 the bridge deck is
acted on by a ‘live’ vertical point load P = 100 kN (i.e. a truck) at the midpoint M.

You will learn more about bridges and dynamics in years 3 & 4, the aim here is to
write structural analysis software and to develop an intuition for dynamic behaviour.
2

Fig. 1 – Idealisation of the Theodor Heuss Bridge in Düsseldorf, Germany
3
The Matlab code that you submit must make use of object-oriented programming
through the creation of a single Class Bridge with appropriate properties and methods
(constructor, assembler, various FE solvers and post-processors / plotters). This Class
should then be instantiated inside other calling scripts, one per Q below. A premium is
placed on minimising code repetition and good code structure. You may find it
beneficial to simply get an FE model running correctly before turning it into a Class.

Static linear stress analysis
This analysis involves the solution of the static equilibrium equation below in the
same way as we have been doing through the module so far (with matrix partitioning).
[ ]{ } { }k d f=
Q1. Perform a static linear stress analysis of the structure under the action of the
midspan point load P. Verify that the reaction forces satisfy global equilibrium.
Report the vertical sag at M and the horizontal sway at E. Plot the deformed shape of
the structure and comment.
Note: the vertical sag at M should be 6.8018 mm and horizontal sway at E should be
1.1586 mm. Make sure you can compute these values exactly (for validation only,
almost nothing in civil engineering is ever known to within 10-4 mm or 0.1 μm!).

Linear frequency analysis
This analysis involves the solution of the eigenvalue problem
[ ]{ } [ ]{ }2FF F FF Fk mφ ω φ=
where ω is the angular frequency of vibration in radians / second, which you should
then transform into a cyclic frequency f = ω/2π in cycles / second (or Hertz) and a
period T = 1/f in seconds. [mFF] is a partitioned sub-matrix of [m], the global mass
matrix, and {ϕF} is a vector of free global dofs representing the vibration mode
associated with the ω vibration frequency. A vibration mode is a deformation state
where the internal elastic forces (caused by stiffness-resisting deformation) are in
balance with the inertia forces (caused by mass-resisting acceleration).

Q2. Perform a linear frequency analysis to compute all of the eigenvalues (how many
are there?) and plot the eigenmodes corresponding to the lowest four frequencies and
comment. Compute also the total weight of the structure.
4
Note: the lowest four vibration frequencies should be 0.47347, 0.63481, 0.97582 and
1.2728 Hz. The highest frequency should be 205.7568 Hz. The structure should weigh
7716581.2295 kg. Make sure you can compute these values exactly.

Explicit response-history analysis
This analysis involves the incremental solution of the dynamic equilibrium equation
[ ]{ } [ ]{ } [ ]{ } { }t t t tm d c d k d f′′ ′+ + =
where {d}t, {d′}t and {d″}t are the nodal displacement, velocity and acceleration
vectors respectively at time t, [k], [c] and [m] are the stiffness, damping and mass
matrices respectively and {f}t is the vector of external forces at time t. Specifically,
solution requires integrating the above to track the evolution of {d}t with time. An
explicit integration scheme uses, say, central differences to approximate the velocity
and acceleration vectors as
{ } { } { }( )12t t t t td d dt +∆ −∆′ ≈ −∆ and { } { } { } { }( )2
1 2
t t t t t t
d d d d
t +∆ −∆
′′ ≈ − +
∆

respectively to express {d}t+Δt (displacement vector at a time step Δt in the future)
explicitly in terms of its present {d}t and previous {d}t–Δt values. Convince yourself
that this gives:
{ } [ ] [ ] { } [ ] [ ] { } [ ] [ ] { }
1
2 2 2
1 1 2 1 1
2 2t t t t t t
d m c f k m d m c d
t t t t t
−
+∆ −∆
      = + − − − −      ∆ ∆ ∆ ∆ ∆      
In what follows, you may assume that:
• The damping matrix [c] is given by
[ ] [ ] [ ]c m kα β= + where α = 0.0583 and β = 3.6×10-5
which are the Rayleigh damping parameters to get 1% damping for relevant
vibration modes (more about this in your 3rd year dynamics module).
• The time increment Δt does not exceed
1
2 cr
t t∆ = ∆ where
max
2
crt ω
∆ =
where ωmax is the highest cyclic vibration frequency (in radians/second) from
Q2. If you exceed Δtcr (a ‘critical’ time step), your solution for {d}t may diverge
because explicit integration schemes are only conditionally stable.
• For t ≤ 0, initialise the displacement vector to zero:
{ } { } { }0 0t t td d= =−∆= =
5
Q3. In Matlab, create a vector of times ts from 0 to 50 seconds in steps of Δt. Create
a second vector fs that computes P·(1 + 0.1·sin(Ωt)) at each time point (i.e. a time-
history record for a ‘cyclic truck’ point load that is vibrating vertically at an angular
frequency of Ω in radians / second). This record is fed into the correct dof of {f}t, all
others zero. You must use [kFF], [cFF], [mFF], {dF} and {fF} during integration.

Perform an explicit response-history analysis of the bridge under a midspan cyclic
point load assuming Ω = 1, 2, 3, 4, 5 and 6 radians / second (separate analysis for
each). Compare the time-history of the structural response (the horizontal sway at E
and the vertical sway at M) under each vibrating frequency on a plot (also show the
static value from Q1) and comment. Hint: you are investigating ‘resonance’.

Q4. Download the ElCentro1940.m and LomaPrieta1989.m files and explore them.
Given an input Δt, they will output horizontal (x-direction) and vertical (y-direction)
ground acceleration records ax(t) and ay(t) in m/s2 for two strong earthquakes that hit
California in the 20th century. Add these into the external loading vector as
{ } [ ] { } ( ) { } ( )( )x x y ytf m a t a tι ι= +
where {ι}x is a vector of 1’s in the locations of only x-direction dofs and 0’s
elsewhere, while {ι}y is a vector of 1’s in the locations of only y-direction dofs and 0’s
elsewhere (these are known as influence vectors). No other loading should be present
(ignore the truck), and you must use the full unpartitioned [k], [c], [m], {d} and {f}
during integration.

Perform an explicit response-history analysis of the bridge for both earthquake
records (separately). Compare the time-history of the structural response (the
horizontal sway at E and the vertical sway at M) under both records on a plot and
comment. Which earthquake is worse for the structure and why do you think that is?

END AJS Jan’ 22 with very special thanks to Alfredo Camara
6