Advanced Structural Analysis and Earthquake Engineering ENG5274
Dr Andrew McBride
University of Glasgow Room 733 Rankine Building andrew.mcbride@glasgow.ac.uk
2021 – Second Semester
Dr Andrew McBride
2021 – Second Semester
1 / 131
Course information
Dr Andrew McBride Course information 2021 – Second Semester 2 / 131
Course information
The course is composed of two main parts:
Part I is on structural elements; Part II on time-dependent problems
The aims of the course are to develop:
Structural elements:
an understanding of the mathematical models for structural elements including beams, plates and shells
an appreciation of boundary conditions and locking-related phenomena
an understanding of how to solve these models using the finite element method
Time-dependent problems: dynamics and earthquake engineering
a knowledge of mathematical models for waves
an understanding of elastodynamics (stress wave propagation)
an understanding and appreciation of the solution of time-dependent problems using
the finite difference and finite element methods
Intended learning outcomes
By the end of the course you should be able to:
apply the finite element method to systematically derive the solution procedure for a range of structural elements with an emphasis on beams
develop a basic finite element code to compute the deflection of a beam and the dynamic response of a rod
appreciate key aspects of earthquake engineering
Dr Andrew McBride Course information 2021 – Second Semester 3 / 131
Contents
The course is broken into two main sections (3–4)
1 Course information
2 The course in pictures
3 The finite element method: a refresher
4 Structural elements
Introduction and overview
Beam theory: notation and geometrically exact theory Euler–Bernoulli beam theory
Timoshenko beams
Summary of plate and shell theory
5 Time-dependent problems: Dynamics Overview
Introduction to wave propagation in elastic media Finite element method for wave propagation
Dr Andrew McBride Course information 2021 – Second Semester 4 / 131
Course management
We will use Moodle to manage the course
All announcements will be on Moodle
All submissions will be electronic and via Moodle
The majority of the lectures are pre-recorded and the links available on Moodle Weekly problem sets (not for submission) will be uploaded to Moodle
Some of the scheduled lecture slots will be used to answer questions and for worked examples
The primary mechanism to ask questions is via the Moodle forum
Due to the size of the class, please ask questions via the forum and not email The complete set of course notes is available on Moodle
These will be updated regularly
Dr Andrew McBride Course information 2021 – Second Semester 5 / 131
Lectures and laboratory sessions
Preliminary material
Lecture 1a Lecture 2a
Course in pictures
FEM for elasticity in 1D – Python code
Structural elements
Lecture 3a Lecture 3b
Overview
Notation and full theory
Time-dependent problems
Dr Andrew McBride
Course information
2021 – Second Semester
6 / 131
Course requirements and assessment
Refer to the course specification document on Moodle for more information
Minimum requirements for award of credits
Must attend the degree examination
Should attend the computational laboratory classes
Must submit the two reports (additional to requirements in course specification)
The programming language Python will be used to illustrate key aspects of the theory and for the reports. A basic understanding of Python and of the finite element method is assumed.
Assessment
Examination 70%
Two reports (assignments) 30% of overall mark (equally weighted)
Dr Andrew McBride Course information 2021 – Second Semester 7 / 131
Reading material and resources
Comprehensive lecture notes will be provided on the Moodle site. These will be updated regularly.
In addition, the following reference books are recommended:
J Fish and T Belytschko, “A First Course in Finite Elements” – this book is available online from the University of Glasgow library
TJR Hughes, “The Finite Element Method: Linear Static and Dynamic Finite Element Analysis”
Software:
We will be using Python version 3.8 or similar
Recommend that you download Anaconda Individual Edition We will also make use of Google Colaboratory
Dr Andrew McBride Course information 2021 – Second Semester 8 / 131
The course in pictures
Dr Andrew McBride The course in pictures 2021 – Second Semester 9 / 131
Course in pictures: structural elements
Beam structures
Shell structures
Approximate geometry, loading and boundary conditions to simplify analysis
3D continuum
2D approximation
1D approximation
Plate structures
Shell structures
Beam structures
Dr Andrew McBride The course in pictures 2021 – Second Semester 10 / 131
Plate structures
direct matrix solving methods are usually used with Newton-Rhapson for non-linearities. This section gives an introductory description of FEM and the relevant equations as used in this project. The underlying mechanics and procedure of formulating the FEM are explained and derived in more detail
C
n
c
t
u
u
e
nts
s
Fig. 4.2. FEM simulation superimposed on an deformation experiment (DHBW Stuttgart 2012) Numerical modelling (FEA)
o
u
llow
r
e
i
p
i
ng
boo
k
s: N
on
li
i
t
ne
ar
so
u
r
li
dm
e
s
:
P
e
cha
nics
:
a
con
a
in
Nonlinear continuum mechanics for finite element analysis[13]. For a description of the approaches used by ANSYS MECHANICAL please refer to the Mechanical APDL documentation [14][15]
the
fo
r
I
–
st
r
c
tin
uum
ap
p
roa
ch
f
or
en
gi
nee
rin
g
t
r
a
l
e
l
m
[1
2] an
d
e
20
Dr Andrew McBride The course in pictures
2021 – Second Semester
11 / 131
Course in pictures: Part I – structural elements
2D approximation
1D approximation
Plate structures Shell structures
Kirchhoff-Love plate theory Reissner-Mindlin plate theory
Euler-Bernoulli* beam theory (slender beams)
nonlinearities & buckling
Beam structures
Timoshenko beam theory (account for shear – stocky beams)
locking (numerical issue)
3D continuum
*Implement FE model in Matlab (project 1)
Dr Andrew McBride The course in pictures 2021 – Second Semester 12 / 131
Course in pictures: Part II – dynamics
liquids
gases
solids
Dr Andrew McBride The course in pictures 2021 – Second Semester 13 / 131
image correlation (DIC) to visualize the transverse wave propagation. In this study, we applied three-point- bending (TPB) technique to Kolsky compression bar to facilitate dynamic transverse loading on a glass
Course i
3-point bendinTghe dinynaamicKTPoBlsexkpyerimceontsmweprerecosndsuciotednwibthathre Kolsky compression bar at Sandia National Laboratories, California. Figure 1 shows the schematic of the testing section. A random pattern was painted on the composite surface for the purpose of applying DIC technique. The high-rate deformation of the composite beam was photographed with a Cordin 550 high speed digital camera. The stain gage signals on the pressure bars provide a velocity/displacement boundary transversally loaded on the composite beam. In this case, the transmitted signal is nearly negligible in comparison to the incident pulse. The displacement boundary on the span side can be considered stationary. The displacement history at the wedge end are calculated with
fi
ber/
ep
o
technique was employed to study the transverse wave
n
p
propagation.
ed
D
IC
i
c
xy
co
m
p
os
ite
be
a
m
.
T
he
hi
tu
r
e
s
:
P
a
r
t
I
I
–
d
gh
-s
pe
y
n
a
m
i
cs
Fig. 1
X =C03t(εI(τ)−εR(τ))dτ ≈2C03tεI(τ)dτ (1) 00
and is shown in Fig. 2. The impact speed is measured as 9.0 m/s.
When the stress wave in the incident bar travels to the composite beam, a longitudinal stress wave is generated and then propagates through the thickness direction. When the longitudinal stress wave arrives at the free back side of the composite beam, it reflects back and doubles the particle velocity. Due to the fast longitudinal stress wave speed and relatively small thickness in the composite beam, it takes only a few microseconds to achieve a nearly equilibrated state in displacement. In other words, the displacement gradient through the thickness direction can be neglected after a few microseconds. However, along the beam direction, there exists a significant gradient in displacement due to the propagation of transverse wave. The DIC results shown in Fig. 3 confirm the transverse wave propagation.
Schematictof the testing section. Fig. 2 Displacement history at the impact wedge end. 0
Y
Fig. 3 Propagation of transverse wave (time interval: 5 microseconds)
Dr Andrew McBride
The course in pictures 2021 – Second Semester 14 / 131
Course in pictures: Part II – dynamics
Dr Andrew McBride The course in pictures 2021 – Second Semester 15 / 131
Course in pictures: Part II – dynamics
3D wave equation
1D wave equation
P-waves S-waves
Computational model for waves in 1D
Finite element method in space Finite difference method in time
Aspects of earthquake engineering
Eigenvalue analysis Stability and accuracy Time integration schemes Numerical implementation*
Dr Andrew McBride The course in pictures 2021 – Second Semester 16 / 131
The finite element method: a refresher
Dr Andrew McBride The finite element method: a refresher 2021 – Second Semester 17 / 131
Key steps in the finite element method
1 Strong form (governing equations and boundary conditions)
2 Weak form
3 Approximation of test function and trial solution
4 Discrete FE system of equations
5 Solution and postprocessing
Dr Andrew McBride The finite element method: a refresher 2021 – Second Semester 18 / 131
Recap: FE formulation an axially-loaded elastic bar
b : external body force per unit length [N/m] p : internal force [N]
equilibrium (balance of forces)
p(x)+b(x+ x/2) x+p(x+ x) = 0 p(x+ x) p(x) +b(x+ x/2) = 0
x
=) dp + b(x) = 0 dx
kinematic (strain – displacement relation)
p(x)
b(x+ x/2)
x
p(x + x)
u
A(x) tx
“(x) = =)
u(x + x) u(x) x
t
⌦
b(x)
du(x) “(x)= dx
[ ] (x) = p(x) [N/m2]
A(x) constitutive relation
[N/m2]
kinetics (stress)
l
governing equation and boundary conditions
d ⇣ du⌘
dx AEdx +b=0
du
n = Edxn = t
u = u
⌦:=0
>0
⌦
<0
tension
n = +1
x
n = t
[N/m2]
n = 1
compression
t>0 t<0
⌦
x = t[ u t \ u = ;
Dr Andrew McBride The finite element method: a refresher 2021 - Second Semester 20 / 131
Weak form in 1D
Strong form for linear elasticity
d AEdu+b=0 dx dx
du σn= Edx n=t
0
■ ForstockybeamDsr(ALtn
cp =
Shear or S (secondary) waves
P waves travel faster than S waves!
Dr Andrew McBride Time-dependent problems: Dynamics 2021 – Second Semester 108 / 131
Where are we?
1 Course information
2 The course in pictures
3 The finite element method: a refresher
4 Structural elements
Introduction and overview
Beam theory: notation and geometrically exact theory Euler–Bernoulli beam theory
Timoshenko beams
Summary of plate and shell theory
5 Time-dependent problems: Dynamics Overview
Introduction to wave propagation in elastic media
Finite element method for wave propagation
Dr Andrew McBride Time-dependent problems: Dynamics 2021 – Second Semester 109 / 131
Recap: FE formulation an axially-loaded elastic bar
b : external body force per unit length [N/m] p : internal force [N]
equilibrium (balance of forces)
p(x)+b(x+ x/2) x+p(x+ x) = 0 p(x+ x) p(x) +b(x+ x/2) = 0
x
=) dp + b(x) = 0 dx
kinematic (strain – displacement relation)
p(x)
b(x+ x/2)
x
p(x + x)
u
A(x) tx
“(x) = =)
u(x + x) u(x) x
t
⌦
b(x)
du(x) “(x)= dx
[ ] (x) = p(x) [N/m2]
A(x) constitutive relation
[N/m2]
kinetics (stress)
l
governing equation and boundary conditions
d ⇣ du⌘
dx AEdx +b=0
du
n = Edxn = t
u = u
⌦:=0
We can estimate ∆tcrit using the maximum eigenvalue λneq from the global eigenvalue problem
this is however expensive! It can be shown that:
λneq ≤ max(λemax) =: λhmax e
where λemax is the maximum eigenvalue of element e, and λhmax is the maximum eigenvalue of the whole system.
Element eigenvalue problem
[Ke −λeMe]Ψe =0 To solve: convert to a standard eigenvalue problem:
[Me]−1Ke −λeIΨe = 0 where Me−1M = I
Ae
=⇒ det [Ae − λeI] = 0 (characteristic polynomial)
Dr Andrew McBride Time-dependent problems: Dynamics 2021 – Second Semester 119 / 131
Example: Eigenvalue problem for linear element and consistent mass matrix
e AeEe 1 −1 e ρle 2 1 (ρAe)le 2 1 K=le−11 andM=612=612.
6E1−1 21e
ρle2 −11−λ12Ψ=0
̃22 d:=6c /le
convert to a standard eigenvalue problem
̃12 −11 −1 1 =⇒ d3 −1 2 −1 1 −λ 0
̃1 −1 1 =⇒ d −1 1 −λ 0
0 e
1 Ψ =0
0 e
1 Ψ =0
− d ̃ ̃2 ̃
d ̃ − λ
=⇒ det −d ̃ d ̃−λ =0 =⇒ (d−λ)(d−λ)−d2 =0
−2dλ+λ = λ(λ−2d) = 0 2 e ̃ 12c2 c2
=⇒ ωmax =λmax =2d= le2 =12 le
̃ ̃ ̃
Note: the critical frequency depends on the ratio of the wave speed to the element length c/le. Its inverse relates to how long it takes a stress wave to pass through an element.
Dr Andrew McBride Time-dependent problems: Dynamics 2021 – Second Semester 120 / 131
Mass lumping
Recall: for problems in 1D we required ngp > (p + 1)/2 Gauss points to fully integrate a polynomial of degree p using Gaussian quadrature.
Consistent mass matrix
If the mass matrix is integrated using a sufficient quadrature rule and computed from an element contribution:
Me= ρNTeNedx Ωe
then it is consistent. For a linear element in 1D, p = 2 in the mass matrix and ngp ≡ 2. Sometimes it’s good to break the rules: mass lumping
a diagonal mass matrix is attractive as it can lead to very efficient algorithms when used in conjunction with explicit time-integration schemes
Dr Andrew McBride Time-dependent problems: Dynamics 2021 – Second Semester 121 / 131
Mass lumping by nodal quadrature Recall: consistent mass matrices
e ρle2 1 e ρle4 2 −1 Mcon=6 1 2 and Mcon=302 16 2
−1 2 4
To construct a lumped element mass matrix, the general idea is to choose a quadrature scheme where the quadrature points are at the nodes of an element.
For example, using the trapezoidal rule with linear shape functions:
Melumped = ρle 1 0 201
Using Simpson’s rule with quadratic shape functions:
e ρle1 0 0 Mlumped = 6 0 4 0
001
More generally we can use a Lobatto quadrature rule.
In both cases were are concentrating the nodal contributions into a single diagonal term!
Dr Andrew McBride Time-dependent problems: Dynamics 2021 – Second Semester 122 / 131
Algorithms for the problem of elastodynamics
24d035 v0 a0
t0
t1
t2
24dn35 ? 24dn+135 vn vn+1 an an+1
tn t tn+1 T
[t0 = 0, …,tn,tn+1,…,T] and ∆t := tn+1 − tn
Initial conditions
d(0) = d0 and v(0) = v0 We can compute a0 from
Ma0 =F(t=0)−Cv0 −Kd0
Discrete problem
Given the state [dn, vn, an] and Fn+1 solve:
Man+1 +Cvn+1 +Kdn+1 =Fn+1
for [dn+1, vn+1, an+1]
Dr Andrew McBride Time-dependent problems: Dynamics 2021 – Second Semester 123 / 131
One-step algorithm (a-form) with a Newmark time-integration scheme Discrete problem
Giventhestate[dn,vn,an]andFn+1 solve: Man+1 +Cvn+1 +Kdn+1 =Fn+1 (⋆)
Define predictors (β and γ determine the stability and accuracy): ∆t2
d = dn +∆tvn + 2 (1−2β)an and v = vn +(1−γ)∆tan . Correctors for Newmark scheme
dn+1 = d + β∆t2an+1 and vn+1 = v + γ∆tan+1 . Finite-difference formulae for evolution of approximate solution:
∆t2 2 dn+1 =dn +∆tvn + 2 (1−2β)an +β∆t an+1
vn+1 = vn + (1 − γ)∆tan + γ∆tan+1
Substitute into (⋆) – solve for an+1 and update state from correctors
M+γ∆tC+β∆t2Kan+1 =Fn+1−Cv−Kd.
A
R
Dr Andrew McBride
Time-dependent problems: Dynamics
2021 – Second Semester
124 / 131
The Newmark family of schemes
Method Average acceleration
Linear acceleration Fox Goodwin Central difference
Type β γ
Stability∗ Unconditional
Ωcrit = 2√3 Ωcrit = √6 Ωcrit = 2
Order accuracy 2
2 2 2
Implicit Implicit Implicit Explicit
11 42 11 62 11
12 2
01 2
∗Stability is based on analysis of
Damping increases the ∆tcrit. Ωcrit = (γ/2 − β)−1/2 (critical sampling frequency).
Stability
a single degree of freedom undamped problem.
Unconditional stability: 2β ≥ γ ≥ 1 2
Conditional stability: γ ≥ 1 , 22
ωh ∆t ≤ Ωcrit ↔ max
β < γ
∆t ≤ Ωcrit T 2π
Dr Andrew McBride Time-dependent problems: Dynamics
2021 - Second Semester
125 / 131
Accuracy
A time-integration scheme is n-th order accurate if the error e is proportional to the time-step size ∆t to the power of n
e(∆t) = C∆tn
All the schemes we considered were second order accurate: the errors decreases by a
where C is a constant.
factor of 4 if we halve the time step size.
Dr Andrew McBride Time-dependent problems: Dynamics 2021 - Second Semester 126 / 131
Explicit versus implicit schemes
Consider an undamped problem C = 0
Explicit scheme: e.g. Central difference (β = 0)
dn+1 = d
Man+1 =Fn+1 −Kd
The displacement dn+1 is computed using only information at tn If M is lumped (diagonal) then no matrix inversion required
leads to a very efficient scheme but requires a small time step
well suited for high strain rate problems (e.g. impact) that require a small time step
Dr Andrew McBride Time-dependent problems: Dynamics 2021 - Second Semester 127 / 131
Stability: some examples
Solved the element eigenvalue problem:
and have λhmax
[Ke −λeMe]Ψe =0 λneq ≤ max(λemax) =: λhmax
e
Assume we are using a central difference scheme: Ωcrit = 2 Two-node linear element with a lumped mass matrix:
ωh =2c thus ∆t≤ 2 =le max le ωmhax c
Two-node linear element with a consistent mass matrix:
h 2√3c 2 le ωmax = le thus ∆t ≤ ωh = √3c
max
Consistent mass matrices tend to yield smaller critical time steps than lumped mass matrices. Another advantage of mass lumping.
Exercise: Compute the critical time step for a three-node quadratic element.
Dr Andrew McBride Time-dependent problems: Dynamics 2021 - Second Semester 128 / 131
Algorithm for elastodynamics using Newmark scheme
Data: Material and load, Domain, Mesh, Time Assemble M and K;
foretone do
M ← Me; K ← Ke;
end
Result: A := M + β∆t2K
Data: Initial conditions d0 and v0 and loading F 0 Solveforinitialaccelerationa0: Ma0 =F0 −Kd0; Assemble R and solve Aan+1 = R for accelerations for all n; forn=1tont do
Data: an ←a,vn ←v,dn ←d d = dn + ∆tvn + ∆t2 (1 − 2β)an;
2
v = vn + (1 − γ)∆tan; R=Fn+1 −Kd;
Solve Aan+1 = R subject to boundary conditions;
dn+1 ← d + β∆t2an+1;
vn+1 ← v + γ∆tan+1; end
Result: Complete time history of problem at all nodal points
Dr Andrew McBride Time-dependent problems: Dynamics 2021 - Second Semester 129 / 131
Some comments on damping
The higher modes of the discrete equations are artefacts of the spatial discretisation (the mesh)
It is generally desirable and often necessary to remove the high-frequency modes The Newmark algorithm itself can introduce algorithmic damping:
γ > 1/2 introduces high-frequency dissipation
for fixed γ > 1/2 , select β to maximise high-frequency dissipation
an ideal value of β = (γ + 1/2)2/4
one drawback: only γ = 1/2 is second-order accurate!
Viscous damping (i.e. C ̸= 0):
Viscous damping damps an intermediate band of frequencies without having a
significant effect on the all-important high modes
Dr Andrew McBride Time-dependent problems: Dynamics 2021 – Second Semester 130 / 131
Stress waves and boundary conditions
A longitudinal stress wave can be compressive (- sign) or tensile (+ sign) A wave will reflect without inverting off a fixed (essential) boundary
A wave will reflect and invert of a free boundary
L
t=0
–
t=L
2c –
t = 3L 2c
t = 5L 2c
–
+
Dr Andrew McBride Time-dependent problems: Dynamics 2021 – Second Semester 131 / 131