Baudynamik 18

Baudynamik

Christian Bucher

Forschungsbereich für Baumechanik und Baudynamik

Institut für Hochbau und Technologie
Technische Universität Wien
Letzte Überarbeitung: 9. März 2018

Vorlesung am 20. März 2018
Modalanalyse
Bewegungsgleichung
Klassische Modalanalyse
Modale Dämpfung
Software slangTNG
http://info.tuwien.ac.at/bucher/Private/slangTNG.html
MacOS
WIN32
Linux
iOS
2/128 © Christian Bucher 2007-2018

Systeme mit mehreren Freiheitsgraden
Beispiel: System mit 2 Freiheitsgraden
x1(t) x2(t)
k1 f1(t) k2 f2(t)
m m
c1 c2
Dynamische Gleichgewichtsbedingungen
x1 f1(t) x2 f2(t)
k1x1 m1ẍ1 k2(x2 – x1) m2ẍ2
c1ẋ1 c2(ẋ2 – ẋ1)
Anordnung in Matrix-Vektor-Form
m1 0 ẍ1 c1 + c2 –c2 ẋ1

" #" # " #" #
+ +
0 m2 ẍ2 –c2 c2 ẋ2
k1 + k2 –k2 x1 f1 (t)
" #" # " #
+ =
–k2 k2 x2 f2 (t)

Allgemeine Form
Für n Freiheitsgrade
Mẍ + Cẋ + Kx = f(t); x(0) = x0 ; ẋ(0) = v0 (1)
x . . . Verschiebungsvektor, Größe n
f . . . Kraft-(Last-)vektor, Größe n
M . . . Massenmatrix, (n × n)
C . . . Dämpfungsmatrix (n × n)
K . . . Steifigkeitsmatrix(n × n)
Alle Matrizen sind symmetrisch and positiv semi-definit (M ist streng
positiv definit)

Interpretation der Systemmatrizen
Elemente sind Kräfte infolge kinematischer Einheitsgrößen

Kjk . . . Kraft am Freiheitsgrad j infolge einer Einheitsverschiebung am
Freiheitsgrad k
Cjk . . . Kraft am Freiheitsgrad j infolge einer Einheitsgeschwindigkeit am
Freiheitsgrad k
Mjk . . . Kraft am Freiheitsgrad j infolge einer Einheitsbeschleunigung am
Freiheitsgrad k
Zufolge Maxwell-Betti können die Indizes j and k vertauscht werden
Die Matrizen K and M können aus einem Finite-Elemente-Modell erstellt
werden
Die Bestimmung von C erfolgt meist auf der Grundlage experimenteller
Daten

Klassische Modalanalyse (1)
Ziel ist dieReduktion der Bewegungsgleichungen (System gekoppelter

Differentialgleichungen in eine einfache lösbare Form)
Einfachste Form: vollständig entkoppelte Gleichungen
Zu jedem Paar symmetrischer und positiv definiter Matrizen K und M der
Größe n × n existiert eine nicht-singuläre Matrix Φ der Größe n × n
(Modalmatrix) mit den Eigenschaften
ΦT MΦ = I; ΦT KΦ = diag(ω 2k ); k = 1...n (2)
Die Spalten φ k der Modalmatrix Φ sind Lösungsvektoren des

Eigenwertproblems
K – ω 2k M φ k = 0;

k = 1...n (3)
Eine nicht-triviale Lösung kann nur dann existieren, wenn
det K – ω 2k M = 0

(4)

Klassische Modalanalyse (2)
ωk . . . Eigenkreisfrequenzen, φ k . . . Eigenschwingungsformen (mode
shapes)
Nun wird eine lineare Transformation auf Grundlage der Modalmaatrix
eingeführt:
x = Φy; ẍ = Φÿ (5)
Setzt man dies in die Bewegungsgleichung (1) ein und multipliziert man
von links mit ΦT , so entsteht
ΦT MΦÿ + ΦT KΦy = ΦT f(t) = g(t)
Die Gleichungen für die neuen Variablen y sind entkoppelt
ÿk + ω 2k yk = gk (t); k = 1...n (6)
Die zugehörigen Anfangsbedingungen ergeben sich durch Inversion von (5)
y = Φ–1 x (7)
→ y0 = ΦT Mx0 ; ẏ0 = ΦT Mv0 (8)

Beispiel - ebener Rahmen (1)
Modell mit linear-elastischen, beidseitig eingespannten Stützen und
starren Deckenplatten
m x2
H EI EI
2m x1
H EI EI
Systemmatrizen
2 0 12EI 4 –2
" # " #
M=m ; K= 3
0 1 H –2 2
Die charakteristische Gleichung (4) ist
12EI 4 –2 2 0
" # " #!
2
det – mω =0
H3 –2 2 0 1
2 3
ω H
Mit der dimensionslosen Variablen µ = m12EI
4 –2 2 0 4 – 2µ –2
" # " #!
det –µ = =0
–2 2 0 1 –2 2–µ
Entwicklung der Determinante
µ2 – 4µ + 2 = 0
Lösungen √
µ 1,2 = 2 ∓ 2
Eigenkreisfrequenzen
r r
EI EI
ω 1 = 2.651 ; ω2 = 6.400
mH3 mH3
Eigenschwingungsformen (massennormiert)
1 0.500 1 0.500
" # " #
φ1 = √ ; φ2 = √
m 0.707 m –0.707

Skizze
φ1 φ2

slangTNG code
1 - -[[
2 Simple plane frame with two DOFs
3 ( c ) 2015 -2016 Christian Bucher , Vienna University of Technology
4 - -]]
5
6 -- System data
7 EI = 1
8 m = 1
9 H = 1
10
11 -- Stiffness matrix
12 K = tmath.Matrix ({
13 {6 , -2} ,
14 { -2 , 2}
15 }) *12* EI / H ^3
16
17 -- Mass matrix
18 M = tmath.Matrix ({
19 {2 , 0} ,
20 {0 , 1}
21 }) * m
22
23 -- Eigensolution
24 eig = t m a t h . M a t r i x E i g e n S y m (K , M )
25 eval = eig : Eigenvalues ()
26 evec = eig : Eigenvectors ()
27
28 -- Take square root of eigenvalues to get omegas
29 omegas = eval : CW () : Sqrt ()
30
31 -- Print results
32 print ( " omegas " , omegas )
33 print ( " evec " , evec )

Hochhaus (1)
Realitätsnahes Tragwerksmodell (Platten- und Balkenelemente, 1452
Knoten, 3840 Elemente, 8448 Freiheitsgrade)
Mode 1 (Biegung/Schub in y) Mode 2 (Biegung/Schub in x)

Hochhaus (2)
Dritte und vierte Schwingungsform
Mode 3 (Torsion um z) Mode 4 (Biegung/Schub in y)

Beispiel
Unsymmetrischer Rahmen
Modell
Eigenschwingungsformen 1–3
Livedemo - slangTNG

Modale Dämpfung (1)
Nach Entkopplung der Bewegungsgleichungen wird häufig “im Nachgang” zu
jeder einzelnen Gleichung in (6) ein modales Lehr’sches Dämpfungsmaß ζ k
festgelegt (z.B. aus Messungen, Erfahrung, etc.) Die entkoppelten Gleichungen
sind damit Form:
ÿk + 2 ζ k ω k ẏk + ω 2k yk = gk (t); k = 1...n (5)
Dies impliziert eine Dämpfungsmatrix, für die gilt (vgl. dazu (1) und (2)):
ΦT CΦ = diag(2ζk ω k ); k = 1...n (6)
Gl.(6) ist speziell gültig für alle C, die gemäß folgender Vorschrift aufgebaut
werden:
C = αK + βM (7)
Diese sogenannte Rayleighdämpfung wird auch als
“Bequemlichkeitshypothese” bezeichnet. Die Anwendung der Gl.(6) darauf
ergibt:
ΦT CΦ = α diag( ω2k ) + β I

Modale Dämpfung (2)
woraus für die Rayleighdämpfung unmittelbar folgt:
1 β
!
ζk = αωk +
2 ωk
Der Verlauf dieser Abhängigkeit ist in der folgenden Abbildung dargestellt. Aus
dieser Kurve wird erkennbar, dass die modalen Dämpfungsmaße mit
zunehmender Eigenkreisfrequenz asymptotisch linear ansteigen, und daher
Werte ζ k > 1 auftreten können. Daher wird in Anwendungen oft nur der
massenproportionale Anteil berücksichtigt, und der steifigkeitsproportionale
Anteil zu Null gesetzt ( α = 0).
ζk
ωk

Freie gedämpfte Schwingung
Für gegebene Anfangsbedingungen x(0) = x0 und ẋ(0) = v0
y0 = ΦT Mx0 ; ẏ0 = ΦT Mv0
Modale Reaktionen
q
ωk0 = ω k 1 – ζ2k
ẏ0,k + ω k y0,k (8)
yk (t) = e–ζk ωk t *y0,k cos ω k0 t + sin ω k0 t+
, ωk0 -
Rücktransformation
x(t) = Φy(t)

Vorlesung am 10. April 2018
Erzwungene Schwingung mittels Modalanalyse

Direkte Integrationsverfahren für die Lösung der Bewegungsgleichungen
Euler’sches Verfahren
Zentrale Differenzen
Newmark’sches Verfahren
Vergleiche

Modalanalyse - Erzwungene Schwingung (1)
Allgemeine geschlossene Lösung: Duhamel-Integral X

Numerische Lösung für einen zeit-diskreten Belastungsvektor: g(ti ) = g(i)
gk
gk(i)
gk(i+1)
∆t
t
ti ti+1
Annahme: Belastung ist in einem Zeitintervall ∆t konstant, Ansatz einer

konstanten Partikulärlösung für jede modale Koordinate
g(i)
yk,p (t) = C = k2
ωk

Modalanalyse - Erzwungene Schwingung (2)
Homogene Lösung
yk,h (t) = e–ζk ωk t A cos ω k0 t + B sin ω k0 t

Ermittlung von A und B aus Anfangs- bzw. Übergangsbedingungen
g(i) ζ k ωk A + ẏ(i)
A = y(i)
k
– k
2
; B = k
ωk ωk
0
Lösung am Ende des Zeitintervalls
g(i)
y(i+1) k
+ e–ζk ωk ∆t A cos ω k0 ∆t + B sin ω k0 ∆t

k
= 2ωk
ẏ(i+1) = e–ζk ωk ∆t
f g
k
– ζ k ω k A + ω k0 B cos ω k0 ∆t – ζ k ω k B + ω k0 A sin ω k0 ∆t
Numerisch effizienter als Duhamel-Integral

Direkte Lösungsverfahren
Vereinfachende Annahmen für x(t) innerhalb eines (kleinen) Zeitintervalls

∆t
Lösung der Differentialgleichungen und diesen Annahmen so, dass für
gegebenes x(t) die Reaktion am Ende des Zeitintervalls x(t + ∆t) berechnet
wird.
Die Gleichungen können als Gleichungen zweiter Ordnung (abhängig von
den Verschiebungen x) oder als Gleichungen erster Ordnung (abhängig von
Verschiebungen x und Geschwindigkeiten ẋ) geschrieben werden.
Die Formulierungen können explizit sein (die Werte am Ende des Intervalls
brauchen nicht bekannt zu sein) oder implizit (die Werte am Ende des
Intervalls müssen bekannt sein, daher muss normalerweise iteriert
werden)
Implizite Methoden sind i.a. numerisch stabiler, erfordern aber höheren
Berechnungsaufwand

Euler’sche Methode
Annahme: x(t) ist eine lineare Funktion von t innerhalb des Zeitintervalls ∆t
x
t
t t + ∆t
Verschiebung und Geschwindigkeit zum Zeitpunkt t + ∆t
x(t + ∆t) = x(t) + ẋ(t)∆t

ẋ(t + ∆t) = ẋ(t) + ẍ(t)∆t
Die Beschleunigung ẍ(t) wird aus der Bewegungsgleichung zum Zeitpunkt t

berechnet.

Beispiel - SDOF-System
Freie ungedämpfte Schwingung, Masse m = 1 kg, Steifigkeit k = 1 N/m2 ,

Anfangsverschiebung x0 = 1, Anfangsgeschwindigkeit v0 = 0
2.0
exact
Displacement [mm]
1.0
∆t = 0.05
0.0
∆t = 0.1
-1.0
∆t = 0.2
-2.0
0 3 6 9 12 15
Time [s]
Energiefehler wächst exponential mit der Zeit

Methode ist nicht langzeitstabil.

Zentrales Differenzenverfahren (1)
Annahme: x(t) ist eine quadratische Funktion der Zeit t im Zeitintervall ∆t.
x
t
t – ∆t t t + ∆t
Geschwindigkeit und Beschleunigung zum Zeitpunkt t
x(t + ∆t) – x(t – ∆t)

ẋ(t) =
2∆t
(9)
x(t + ∆t) – 2x(t) + x(t – ∆t)
ẍ(t) =
∆t2

Zentrales Differenzenverfahren (2)
Zusammen mit (1) ausgewertet zum Zeitpunkt t erhalten wir:
1 1
2

M + C x(t + ∆t) = f(t) – K – M x(t)
∆t2 2∆t ∆t2
1 1

– 2
M– C x(t – ∆t)
∆t 2∆t
Für gegebene Werte x(t) und x(t – ∆t) können die Werte x(t + ∆t) daraus
durch Lösen eines linearen Gleichungssystems berechnet werden.
Für den Start des Algorithmus müssen die Werte x(–∆t) bekannt sein.
Dafür werden die Anfangswerte x0 und v0 zum Zeitpunkt t = 0
herangezogen. Aus (9):
∆t2
x(–∆t) = x0 – v0 ∆t + ẍ(0)
2
Die Beschleunigung ẍ(0) wird aus (1) für t = 0 ermittelt:
ẍ(0) = –M–1 Cv0 + Kx0 – f(0)



1.0
exact
Displacement [mm]
0.5
∆t = 0.5
0.0
∆t = 1.0
-0.5
∆t = 2.0
-1.0
0 3 6 9 12 15
Time [s]
Gesamtenergie bleibt über die Zeit konstant

Frequenzfehler nimmt mit der Größe des Zeitschritts zu (die Frequenz
steigt an)
DEMO: slangTNG

Numerische Stabilität (1)
Das zentrale Differenzenverfahren ist nur bedingt stabil, dh. für ∆t > ∆tcrit
divergiert die numerische Lösung exponentiell.
Der kritische Zeitschritt ist gegeben durch
2
∆t ≤ ∆tcrit =
ω0
Zur Herleitung wird das rekursive zentrale Differenzenschema für die

ungedämpfte frei Schwingung eines SDOF-Systems (Masse m, Steifigkeit
k)als Matrix-Vektor-Operation geschrieben:
k 2
x(t + 2∆t) 2– m ∆t –1 x(t + ∆t) x(t + ∆t)
" # " #" # " #
= =Z
x(t + ∆t) 1 0 x(t) x(t)
Zur Berechnung von mehreren Zeitschritten wird dieses Schema mehrfach

ausgeführt, dh. der Zustandsvekrot wird immer wieder mit der Matrix Z
multipliziert. Dies führt zu einem exponentiellen Wachstum der
berechneten Lösung, wenn mindestens ein Eigenwert λmax der Matrix Z
einen Betrag |λmax | > 1 besitzt.
Numerische Stabilität (2)
Die Eigenwerte λ1,2 sind gegeben durch
r
a a2 k 2
λ1,2 = 1 – ± –a + ; a= ∆t
2 4 m
Für a < 4 sind die Eigenwerte komplex, und wir erhalten Eigenwerte vom
Betrag 1:
a 2 a2

|λ1,2 |2 = 1 – +a– =1
2 4
Für den Fall a ≥ 4 sind die Eigenwerte reell. Dann ist λ2 jener mit dem
größeren Betrag und aus
r r
a a2 a a2
|λ 2 | = 1 – – –a + = + –a + –1>1
2 4 2 4
erhalten wir die Bedingung für numerische Instabilität:

r
a a2 m 2
r
+ –a + > 2 → a > 4 also ∆t > 2 =
2 4 k ω0

Newmark’sche Methode
Implizites Verfahren, das unbedingt stabil ist
Kinematische Annahmen:
∆t
ẋ(t + ∆t) = ẋ(t) + [ẍ(t) + ẍ(t + ∆t)] (10)
2
∆t2
x(t + ∆t) = x(t) + ẋ(t)∆t + [ẍ(t) + ẍ(t + ∆t)]
4
Zusammen mit (1) zum Zeitpunkt t + ∆t:
Mẍ(t + ∆t) + Cẋ(t + ∆t) + Kx(t + ∆t) = f(t + ∆t)
ergibt sich für lineare Systeme
4 2

K+ M+ C x(t + ∆t) = f(t + ∆t)
∆t2 ∆t
(11)
4 4
2
+M x(t) + ẋ(t) + ẍ(t) + C x(t) + ẋ(t)
∆t2 ∆t ∆t


1.0
exact
Displacements [mm]
0.5
∆t = 0.5
0.0
∆t = 1.0
-0.5
∆t = 2.0
-1.0
0 3 6 9 12 15
Time [s]
Gesamtenergie bleibt über die Zeit konstant

Frequenzfehler nimmt mit dem Zeitschritt zu (abnehmende Frequenz)
DEMO: slangTNG

Vorlesung am 15. Mai 2018
Nichtklassische Modalanalyse
Nichtlinearitäten
Geometrische Nichtlinearität
Physikalische Nichtlinearität
Kontakt
Numerische Verfahren
Newmark’sches Verfahren
Newton-Raphson-Iteration

Nichtklassische Modalanalyse - Motivation
Zwei ungekoppelte ungedämpfte SDOF-Systeme

Kopplung durch viskosen Dämpfer
x1(t) x2(t)
f1(t) f2(t)
m1 m2
(a) k1 k2
(b) m1 m2
k1 c k2
Bewegungsgleichungen
m1 0 ẍ1 c –c ẋ1 k1 0 x1 f1
" #" # " #" # " #" # " #
+ + =
0 m2 ẍ2 –c c ẋ2 0 k2 x2 f2
Dämpfung bewirkt Kopplung

Sonderfall
Identische Teilsysteme (k1 = k2 , m1 = m2 )

Freie Schwingungen ohne Anfangsgeschwindigkeiten
Fall 1: Anfangsverschiebungen x1 (0) = x2 (0) = x0
führt zu vollständig synchroner Bewegung beider Systeme, Dämpfer ist
nicht aktiv.
Fall 2: Anfangsverschiebungen x1 (0) = –x2 (0) = x0
führt zu gegenläufiger Bewegung, Dämpfungseffekt ist maximal.
Unter der Annahme modaler Dämpfung würde der Koppelterm
vernachlässigt werden.
Damit wären Fall1 und Fall 2 gedämpft

Freie Schwingung (1)
Formulierung der Bewegungsgleichungen im Phasenraum

(Verschiebungen/Geschwindigkeiten)
System von 2n Differentialgleichungen erster Ordnung
" # " # " #" #
ẋ ẋ 0 I x
ẏ = = =
ẍ –M–1 Càx – M–1 Kx –M–1 K –M–1 C ẋ
Exponentialansatz
y = φ eλt ; ẏ = λφ eλt ; φ . . . 2n – dimensionalerVektor
Führt zum Eigenwertproblem
(G – λI) φ = 0
mit 2n Lösungen (λk , φ k ), k = 1 . . . 2n

Freie Schwingung (2)
Zusammenfassung der rechtsseitigen Eigenvektoren φ k spaltenweise in
der komplexen Modalmatrix Φ mit der Eigenschaft
Φ–1 GΦ = diag(λk ); k = 1 . . . 2n
Hier wird vorausgesetzt, dass Φ invertierbar ist.

Einführung neuer Variabler
y = Φz; z = Φ–1 y
ergibt entkoppelte Differentialgleichungen

Lösung des Anfangswertproblems mit y(0) = y0
y(t) = Φ diag(eλk t )Φ–1 y0 = Z(t)y0
Die reellwertige Matrix
Z(t) = Φ diag(eλk t )Φ–1
heisst Fundamentalmatrix. Sie genügt der Differentialgleichung Ż = GZ.

Erzwungene Schwingung
Im Phasenraum sind die Bewegungsgleichungen
" # " #
0 I 0
ẏ – y=
–M–1 K –M–1 C M–1 f(t)
In kompakter Form
ẏ – Gy = g(t)
Partikulärlösung (Variation der Konstanten)
t t
yp (t) = Z(t) –1
Z (τ)g(τ)dτ = Φ diag(eλk (t–τ) )Φ–1 g(τ)dτ
0 0
Sonderfall: Konstante Anregung im Zeitintervall ∆t, d.h.

g(τ) = gi ; 0 ≤ τ ≤ ∆t:
1
yp (t) = Φ diag(1 – eλk ∆t ) diag( )Φ–1 gi = T(∆t)gi
λk
Partikuläre und homogene Lösung
yi+1 = Z(∆t)yi + T(∆t)gi

Beispiel - 2DOF System unter harmonischer Anregung
Harmonische Anregung auf Masse m2 : f2 (t) = f0 sin ω t mit Zahlenwerten

f0 = 1 N, m1 = m2 = 1 kg, k1 = k2 = 1N/m, c = 0.1 Ns/m
Systemmatrix G
 0.00 0.00 1.00 0.00 
0.00 0.00 0.00 1.00
 
G =  
 –1.00 0.00 –0.02 0.02 
 0.00 –1.00 0.02 –0.02 
Eigenwerte von G
λ1,2 = –0.02 ± 0.9998 ı; λ3,4 = ± ı
Zwei Realteile sind Null → ungedämpfte Schwingung ist möglich

Lösung
Belastung in Resonanz mit dem ungedämpften Primärsystem

Dämpfung bewirkt Schwingung im Sekundärsystem
40
X1 X2
Response X1, X2 [m]
20
-20
-40
0 20 40 60 80 100
Time t [s]

Übungsbeispiel - freie Schwingung
Zweifreiheitsgradsystem wie zuvor

x1(t) x2(t)
m m
k1 c k2
Systemdaten: m1 = m2 = 1 kg, k1 = k2 = 1N/m

Anfangsbedingungen x1 (0) = –x2 (0) = x0
Gesucht: Dämpferkonstante c so, dass die Verschiebungen x1 (t), x2 (t)
kleiner sind als 0.01 x0 für t > 60 s.
Lösung: c > 0.0384 Ns/m

Nichtlineare Systeme
Geometrische Nichtlinearität (vor allem Stabilitätseffekte)

Werkstoffnichtlinearität (Plastizität, fortschreitende Schädigung)
Kontakt (wechselnde kinematische Bindungen)
Plastizität bringt
erhöhte Verformungen → Steifigkeit wird reduziert
Energiedissipation → Dämpfung wird vergrößert
Superpositionsprinzip ist nicht mehr gültig, daher kann die Modalanalyse
nicht angewandt werden
Nur direkte Integrationsverfahren für die Lösung der Bewegungsgleichung
verfügbar

Geometrische Nichtlinearität
Momentengleichgewicht am verformten
System
FH cos ϕ + mgH sin ϕ = kϕ ϕ

mg
Für kleine Winkel ϕ gilt cos ϕ ≈ 1, F
sin ϕ ≈ ϕ , und daher m
kϕ
!
F= – mg ϕ = keff ϕ
H
H ϕ
kϕ
Die effektive Steifigkeit wird durch das
Eigengewicht reduziert. Dies wird
allgemein als P – ∆-Effekt bezeichnet.
Die obige Linearisierung entspricht der
Betrachtung nach Theorie II. Ordnung.

Einfaches Plastizitätsmodell
linear elastisch - plastisch mit kinematischer Verfestigung
σ
σu
αE
βF
σ0
E
ε
σl
Werkstoff hat unter zyklischer Belastung eine stabile Hysterese

Spannung ist eine Funktion der Dehnung und der
Dehnungsgeschwindigkeit

Infinitesimale
Spannungs-Dehnungs-Beziehung
Fall A (elastisches Inkrement)
σl (ε ) < σ < σu (ε ) : dσ = Ed ε
Fall B (an der oberen plastischen Grenze)
 ε̇ > 0 : dσ = α Edε

σ = σu (ε ) : 
 ε̇ ≤ 0 : dσ = Ed ε

Fall C (an der unteren plastischen Grenze)
 ε̇ ≥ 0 : dσ = Ed ε

σ = σl ( ε) : 
 ε̇ < 0 : dσ = α Edε

Wichtig: Spannung darf plastische Grenzen niemals über/unterschreiten

Bewegungsgleichung
Die Gleichgewichtsbedingungen an den Knoten egeben
Mẍ + Cẋ + r(x, ẋ) = F(t)
Dabei hängt der Vektor der inneren Kräfte r oft sowohl von der
Verschiebung als auch der Geschwindigkeit ab
Darüber hinaus können auch noch interne Variable eine Rolle spielen (z. B.
plastische Dehnungen/Verformungen), die nicht explizit in der
Bewegungsgleichung auftauchen, aber die internen Kräfte beeinflussen.
Prinzipiell muss dieses Differentialgleichungssystem iterativ gelöst werden.
Außer in Sonderfällen stehen dafür nur rein numerische Verfahren zur
Verfügung.

Beispiel - SDOF-System mit
linear-elastisch/ideal-plastischer Feder
Plastizität wird über ein Reibungselement mit Grenzkraft r dargestellt,
kinematische Verfestigung durch Feder k1
Rutschen (dh. plastische Deformation) tritt ein, wenn |k2 (x – z)| ≥ r
x(t)
k1
m f(t)
r k2
z(t)
Akkumulierte plastische Deformationen werden durch eine

Zustandsvariable z beschrieben. Diese genügt der Gleichung ż = h(x, ẋ, z)ẋ.
Dabei ist die Funktion h() gegeben durch:


 0 |k2 (x – z)| < r

0 k2 (x – z) ≥ r ∧ ẋ < 0



h=



 0 k2 (x – z) ≤ –r ∧ ẋ > 0

1 sonst


slangTNG code (1)
1 - -[[
2 Elastic - plastic SDOF system
3 Harmonic response using explicit forward Euler method
4 ( c ) 2012 -2016 Christian Bucher , TU Wien
5 - -]]
6
7 -- Function to compte the derivatives of the state variables
8 function derivative (y , f )
9 local x = y [0];
10 local v = y [1];
11 local z = y [2];
12 local fr = k2 *( x - z ) ;
13 local ft = k1 * x + fr ;
14 local vd =f - ft ;
15 local xd = v ;
16 if ( math.abs ( fr ) <r or fr > r and v <0 or fr < - r and v >0) then zd =0;
17 else zd = v ; end
18 local ydot = tmath.Matrix (3)
19 ydot [0] = xd ;
20 ydot [1] = vd ;
21 ydot [2] = zd ;
22 return ydot
23 end
24
25 -- Main program starts here
26 k1 = .1 ;
27 k2 = .9 ;
28 r =1;
29 f0 = 1
30
31 N = 5000
32 dt = 0 .0125
33 x = tmath.ZeroMatrix (3 ,1)
34 resp = tmath.ZeroMatrix (N , 4)
35 for i =1 ,N -1 do
36 t = i * dt ;

slangTNG code (2)
37 force = f0 * math.sin ( t ) ;
38 xd = derivative (x , force ) ;
39 x = x + xd * dt
40 resp [{ i ,0}] = t
41 resp [{ i ,1}] = x [0]
42 resp [{ i ,2}] = x [1]
43 resp [{ i ,3}] = x [2]
44 end
45 v = graph.Graph ( " Response " , " Bright " )
46 v : Plot ( resp : GetCols (0) , resp : GetCols (1) , 2 , " x " )
47 v : Plot ( resp : GetCols (0) , resp : GetCols (3) , 2 , " z " )

Ergebnis
Anfänglich verhält sich das System wie ungedämpft (Resonanz)

Bei größeren Verschiebungen rutscht das Reibungselement durch
→ plastische Verformungen
Verschiebungsreaktion erreicht infolge der Energiedissipation durch
Plastizität einen stationären Zustand
20
x
Displacement [m]
10
z
0
-10
-20
0 10 20 30 40 50
Time [s]

Newmark-Verfahren
für nichtlineare Systeme
Innere Kräfte werden nun nicht mehr als Produkt von Steifigkeitsmatrix
und Verschiebungsvektor ausgedrückt, sondern direkt durch Integration
über die Finiten Elemente erhalten.
Kx wird ersetzt durch r(x, ẋ).
Es entsteht das Gleichungssystem
4 2

r x(t + ∆t), ẋ(t + ∆t) + + x(t + ∆t) = f(t + ∆t)

M C
∆t2 ∆t
(12)
4 4
2
+M x(t) + ẋ(t) + ẍ(t) + C x(t) + ẋ(t)
∆t2 ∆t ∆t
Dieses System kann nur iterativ gelöst werden (üblicherweise durch
Newton-Raphson-Iteration).

Newton-Raphson-Iteration
Laststeigerung in Schritten ∆f
Verschiebungsinkremente ∆x werden so berechnet, dass die
Gleichgewichtsbedingungen erfüllt sind (innere Kräfte r sind gleich groß
wie die äußeren Lasten f).
Erfordert die Berechung und Faktorisierung der tangentiellen
Steifigkeitsmatrix KT
r r
original modified
r1 k(1) k(2) r1 k(1)
r(2) r(2)
r(1) r(1)
r0 r0
x(1) x(2) x(1) x(2)

x x
x0 x1 x0 x1

Iterationsschema
1. Start von einer Gleichgewichtslage x0 für die gilt: r0 = f0
2. Steigerung der aufgebrachten Belastung f = f0 + ∆f, die
Ungleichgewichtskraft ist dann ∆r = ∆f,
3. Berechnung der tangentiellen Steifigkeitsmatrix abhängig von der
aktuellen Iteration der Verschiebungen
∂r
KT =
∂x x=x0
4. Vorausberechnung des Verschiebungsinkrements am Ende des Lastschritts

(Lösung eines linearen Gleichungssystems)
∆x(1) = K–1
T ∆r, x
(1)
= x0 + ∆x(1)
5. Berechnung der inneren Kräfte r(1) und Neuberechnung der

Ungleichgewichtskraft
∆r(1) = f – r(1)
6. Wiederholung ab Schritt 3 (oder ab Schritt 4 beim modifizierten
NR-Verfahren) bis die Ungleichgewichtskraft ausreichend klein wird.
Einfacher Rahmen - Modelldefinition
Geometrie
F(t) = F0 sin νt
1.50 x 0.25
2x8
0.30 x 0.30
10
Material Stützen: linear-elastisch/ideal-plastisch mit 1% kinematischer

Verfestigung (E = 30 GPa, ν = 0.15, β y = 20 MPa, ρ = 2500 kg/m3 .
Material Riegel: linear-elastisch, E = 30 GPa, ν = 0.16, ρ = 2500 kg/m3 .

Einfacher Rahmen - slangTNG code (1)
1 - -[[
2 This function defines a 2 - story frame structure
4 - -]]
5
6 -- define the frame in terms of a function
7 function frame1 ()
8 structure = fem.Structure ( " 2 story frame " )
9 blue = tmath.Matrix ({{0 , 100 , 170 , 255}}) -- define blue color
10 orange = tmath.Matrix ({{255 , 100 , 0 , 255}}) -- define orange color
11
12 s_column = structure : AddSection (101 , " RECT " , 0) ; -- define column cross section with ID 101
13 s_column : SetData ( tmath.Matrix ({{0 .3 , 0 .3 }}) )
14 s_column : SetColor ( blue )
15
16 s_beam = structure : AddSection (102 , " RECT " , 0) ; -- define beam cross section with ID 102
17 s_beam : SetData ( tmath.Matrix ({{0 .2 , 1 .5 }}) )
18 s_beam : SetColor ( orange )
19
20 m_column = structure : AddMaterial (201 , " E LA ST IC _PL AS TI C_ 1D " )
21 m_column : SetData ( tmath.Matrix ({{3 e10 , .16 , 2500 , 2 e7 , 0 .01 }}) )
22 m_beam = structure : AddMaterial (202 , " LINEAR_ELASTIC " )
23 m_beam : SetData ( tmath.Matrix ({{3 e10 , .16 , 200}}) )
24
25 -- nodes
26 B = 10; H = 8 -- width and height of one story
27 nodes = tmath.Matrix ( { -- define nodes
28 {1001 , 0 , 0 , 0} ,
29 {1002 , B , 0 , 0} ,
30 {1003 , 0 , 0 , H } ,
31 {1004 , B , 0 , H } ,
32 {1005 , 0 , 0 , 2* H } ,
33 {1006 , B , 0 , 2* H } ,
34 {1100 , B /2 , 0 , 0} -- This is the reference node for all beams
35 })
36 structure : AddNodes ( nodes ) -- Add all nodes

37
38 -- elements
39 columns = tmath.Matrix ({
40 {2001 , 1001 , 1003 , 1100} ,
41 {2002 , 1002 , 1004 , 1100} ,
42 {2003 , 1003 , 1005 , 1100} ,
43 {2004 , 1004 , 1006 , 1100}
44 })
45 structure : AddElements ( " RECT " , 201 , 101 , columns ) -- Add column elements
46
47 beams = tmath.Matrix ({
48 {3001 , 1003 , 1004 , 1100} ,
49 {3002 , 1005 , 1006 , 1100}
50 })
51 structure : AddElements ( " RECT " , 202 , 102 , beams ) -- Add beam elements
52
53 -- support conditions and global degrees of freedom
54 structure : SetAvailDof ( tmath.Matrix ({{1 , 0 , 1 , 0 , 1 , 0}}) ) -- plane problem (x - z )
55 n1001 = structure : GetNode (1001)
56 n1001 : SetAvailDof ( tmath.Matrix ({{0 , 0 , 0 , 0 , 0 , 0}}) ) -- Clamped at node 1001
58 n1002 : SetAvailDof ( tmath.Matrix ({{0 , 0 , 0 , 0 , 1 , 0}}) ) -- Hinge at node 1002
60 n1100 : SetAvailDof ( tmath.Matrix ({{0 , 0 , 0 , 0 , 0 , 0}}) ) -- No dof at reference node
61
62 -- Assemble all degrees of freedom
63 ndof = structure : GlobalDof ()
64 return structure
65 end

Einfacher Rahmen - dynamische Analyse
Kurve basiert auf ausgedehnten plastischen Zonen in den Stützen
Top displacement
-40.0 -30.0 -20.0 -10.0 0.0 10.0 20.0 30.0 40.0

Load [N] (*E3)
-30.0 -20.0 -10.0 0.0 10.0 20.0 30.0

Displacement [m] (*E-2)

slangTNG code (1)
1 - -[[
2 Dynamic analysis of simple linear structures
4 All rights reserved , free for educational use.
5 - -]]
6
7 -- load the definition of the frame structure
8 dofile ( " frame1.tng " )
9
10 -- Consider geometrical nonlinearity ?
11 geo = 0
12 F0 = 50000
13
14 -- create the structure
15 struct = frame1 ()
16 -- Assemble stiffness and mass matrices ( sparse )
17 KK = struct : SparseStiffness ()
18 MM = struct : SparseMass ()
19
20 -- Compute first 2 eigenvalues and mode shapes
21 val , vec = KK : Eigen ( MM , 2)
22 omega = math.sqrt ( val [0]) / math.pi
23 print ( " omega " , omega )
24 nu = 0 .9 * omega
25
26 -- Construct load vector for dynamic load
27 F1 = struct : Ge t A ll D i sp l a ce m e nt s ()
28 F1 : SetZero ()
29 load_node = 1006
30 load_index = struct : GetNodeIndex ( load_node )
31 F1 [{ load_index , 0}] = F0 -- load in x - direction
32
33 -- Construct load vector for gravity
35 F2 : SetZero ()
36 col = F2 : GetCols (2)

slangTNG code (2)
37 F2 : SetZero ()
38 col : SetConstant (1)
39 F2 : SetCols ( col , 2)
40
41 -- Bring load cases into vector form and multiply unit acceleration with mass matrix
42 F = struct : T oD ofD is pl ac em en ts ( F1 )
43 F_G = struct : To Do fD is pl ace me nt s ( F2 )
44 G = MM * F_G *( -9 .81 ) *12
45 UG = KK : Solve ( G )
46 print ( " UG " , UG )
47 struct : S et D o fD i s pl a c em e n ts ( UG )
48 if ( geo ==1) then
49 KG = struct : S pa rs eG eo St if fn es s ()
50 beig , bvec = KK : EigenBuckling ( KG , 4 , -1)
51 print ( " beig " , beig )
52 struct : Se t D of D i sp l a ce m e n ts ( UG *0)
53 end
54 if ( geo ==1) then KK = KK : Add ( KG , 1) end
55
56 T = 30
57 dt = 0 .05
58 N = T / dt
59
60 -- Prepare Newmark
61 a0 = 4/ dt ^2
62 a1 = 2/ dt
63 a2 = 4/ dt
64 a3 = 1
65 a4 = 1
66 a5 = 0
67 a6 = dt /2
68 a7 = dt /2
69
70 -- Effective " stiffness " for Newmark method
71 Keff = KK : Add ( MM , a0 )
72

slangTNG code (3)
73 U = tmath.Matrix ( F )
74 U : SetZero ()
75 V = tmath.Matrix ( U )
76 A = tmath.Matrix ( V )
77
78 -- Prepare graphics
79 g = graph.Graph3D ( " Frame " )
80 g : Rotate ( -60 , 1 , 0 , 0)
81 plot = struct : Draw (1)
82 g : Triangles ( plot )
83 g : Autoscale ()
84
85 disp_load = tmath.Matrix (N , 2)
86 disp_load [{0 ,0}] = 0
87 disp_load [{0 ,1}] = 0
88 supportlist = tmath.Matrix ({{0 , 1}})
89 topnode = 4
90
91 for i =0 ,N -1 do
92 factor = math.sin ( nu * i * dt )
93 R1 = F * factor + MM *( A + V * a2 + U * a0 )
94 U1 = Keff : Solve ( R1 )
95 for k =0 ,5 do
96 struct : Se t D of D i sp l a ce m e n ts ( U1 )
97 V1 = U1 * a1 - U * a1 - V
98 struct : SetDofVelocities ( V1 )
99 R = R1 - struct : GlobalResForce () - MM * U1 * a0
100 if ( geo ==1) then R = R - KG * U1 end
101 Rnorm = tmath.Norm ( R )
102 if ( Rnorm < 1) then break end
103 DU = Keff : Solve ( R )
104 U1 = U1 + DU
105 end
106 V1 = U1 * a1 - U * a1 - V
107 A1 = V1 * a1 - V * a1 - A
108 U = tmath.Matrix ( U1 )

slangTNG code (4)
109 V = tmath.Matrix ( V1 )
110 A = tmath.Matrix ( A1 )
111 struct : S et D o fD i s pl a c em e n ts ( U )
112 struct : SetDofVelocities ( V )
113
114 struct : GlobalUpdate ()
115
116 plot = struct : Draw (1) -- Deformation with scale 1
117 g : Clear ()
119 g : Render ()
120
121 disp = struct : G et A l lD i s pl a c em e n ts ()
122 disp_load [{ i ,0}] = disp [ topnode ]
123 disp_load [{ i ,1}] = F0 * factor
124
125 end -- for
126
127 v2 = graph.Graph ( " Response " , " Bright " )
128 v2 : AxisLabels ( " Displacement [ m ] " , " Load [ N ] " )
129 v2 : Plot ( disp_load : GetCols (0) , disp_load : GetCols (1) , 1 , " Top displacement " )
130 v2 : PDF ( " f ra m e 1_ r e sp o n se . p df " )

Vorlesung am 29. Mai 2018
Zufallsvariable
Autokovarianzfunktion
Leistungsspektraldichtefunktion
Spektrale Schätzung
Reaktion linearer Systeme auf zufällige Anregung

Leistungsspektralmethode und Näherung durch weißes Rauschen

Zufallsvariable (1)
Wir betrachten das Ereignis A dass eine zufällig beobachtete Größe X
kleiner ist als ein gegebener Wert x
A = {X|X < x}
Wahrscheinlichkeit von A hängt vom Wert von x ab
Prob[A] = F(x)
Wahrscheinlichskeitsverteilung von X
FX (x) = Prob[X < x]
Grenzwerte
lim FX (x) = 0; lim FX (x) = 1
x→–∞ x→+∞
Wahrscheinlichkeitsdichtefunktion von X
d
fX (x) = FX (x)
dx
Zufallsvariable (2)
Es folgt
∞
fX (x)dx = 1
–∞
Grafische Darstellung
Probability Distribution Function

0.6 1.00
Probability Density Function
fX(x) FX(x)
0.5 0.75
0.3 0.50
0.2 0.25
0.0 0.00
-1.0 0.0 1.0 2.0 3.0
Value of Random Variable

Momente einer Zufallsvariablen (bis zur 2. Ordnung
Mittelwert (Erwartungswert)
∞
µX = E[X] = xfX (x)dx
–∞
Varianz (Quadrat der Standardabweichung)
∞
σX2 = E[(X – µ X ) ] = (x – X̄)2 fX (x)dx
2
–∞
Der Erwartungswertoperator E[Z] bezeichnet einen sogenannten

Ensemblemittelwert, d.h. es wird über alle theoretisch möglichen
Realisationen der Größe Z gemittelt.
E ist ein linearer Operator
E[Y + Z] = E[Y] + E[Z], E[λZ] = λE[Z]; λ ∈ Ò und deterministisch

Zufallsprozess
Ein Zufallsprozess X(t)ist das Ensemble aller Realisationen X(t, γ) (t . . . Zeit,

γ . . . Zufall)
Zu jedem Zeitpunkt t können die statistischen Eigenschaften des Prozesses
X(t) verschieden sein
Ensemblemittelwerte zu den Zeitpunkten t und s
X(t, γ)
t
s t, s

Prozessstatistik
Mittelwertfunktion
µ X (t) = E[X(t)]
Autokovarianzfunktion
RXX (t, s) = E[(X(t) – µ X (t))(X(s) – µ X (s))]
Sonderfall t = s
RXX (t, t) = E[(X(t) – µ X (t))2 ] = σX2 (t)
Verteilungsfunktion für einen Zeitpunkt
FX(t) (x) = Prob[X(t) < x]
Verteilungsfunktion für mehrere Zeitpunkte
FX(t1 ),X(t2 ),...X(tn ) (x1 , x2 , . . . xn ) = Prob[(X(t1 ) < x1 ) ∧ (X(t2 ) < x2 ) ∧ . . .

∧ (X(tn ) < xn )]
Gauß’scher Prozess: alle Verteilungsfunktionen sind normal.

Schwach stationärer Prozess
Ein Zufallsprozess heißt schwach

stationär wenn
µX (t) = µX = const. RXX

RXX (t, s) = RXX (s – t) = RXX (τ)
σx2
Für einen derartigen Prozess gilt
τ
RXX (τ) = RXX (–τ)
max |RXX (τ)| = RXX (0) = σX2
τ∈Ò –σx2
Die Autokovarianzfunktion ist

symmetrisch in τ and hat ihren
größten Absolutwert bei τ = 0

Leistungsspektraldichtefunktion
Intiuitiv kann erwartet werden, dass lim RXX (τ) = 0
τ→±∞
Bei ausreichend schneller Konvergenz gegen 0 existiert die
Fouriertansformation (Leistungsspektraldichte, power spectral density
function, PSD)
∞
1
SXX ( ω ) = RXX (τ)eiωτ dτ
2π
–∞
Die inverse Beziehung ist
∞
RXX (τ) = SXX ( ω )e–iωτ d ω
–∞
Sonderfall τ = 0
∞
σX2 = RXX (0) = SXX ( ω )d ω
–∞
PSD kann als Verteilung der Varianz 2
σX über die Frequenzen interpretiert
werden.
Bandbreite eines Zufallsprozesses
Breitbandprozess: Varianz ist über einen großen Frequenzbereich verteilt,
die Autokovarianzfunktion fällt schnell gegen Null ab.
Schmalbandprozess: Varianz ist eng um eine Frequenz konzentriert, die
Auokovarianzfunktion geht langsam gegen Null
SXX RXX
Wide band
ω τ
SXX RXX
Narrow band
ω τ

Spektrale Schätzung
PSD kann mittels Fouriertransformation direkt aus Zeitreihen

(Messwerten) ermittelt werden
Diskrete finite Fouriertransformation (meist FFT)
T
X̂( ω , T) = x(t) exp(i ω t)dt
0
PSD S is estimated from the periodogram P
1 f
E |X̂( ω , T)|2 = lim PXX ( ω , T)
g
SXX ( ω ) = lim
T→∞ T T→∞
Erfordert mehrere lange Zeitreihen

Beispiel - Windaufzeichungen
Stündliche Mittelwerte der Windgeschwindigkeit am Logan Airport

(Boston) im Jahr 2008 (8784 data points)
16
Wind speed [m/s]
12
0
0 5 10 15 20 25 30 35
Time t [s·106]

Schätzung der PSD
Zerlegung der Zeitreihe in 12 Teile und Mittelung der 12 Periodogramme)
70
60
PSD [m2/s · 103]
50
40
30
20
10
0
0 2 4 6 8 10
Circular frequency ω [rad/s·10–4]

Unterer Frequenzbereich
Spitzen bei 1.67·10–5 and 7.12·10–5 rad/s
70
60
PSD [m2/s · 103]
50
40
30
20
10
0
0 2 4 6 8 10
Circular frequency ω [rad/s·10–5]

Dynamische Reaktion eines SDOF-Systems
System mit deterministischen Eigenschaften m, c, k unter einer schwach

stationären Zufallsanregung F(t)
Mittelwert F̄, Autokovarianzfunktion RFF (τ)
Reaktion X(t) ist ein Zufallsprozess
Wir wollen Erwartungswert und Autokovarianzfunktion von X(t) berechnen
X(t)
k
m F(t)
c

Erwartungswert (1)
Formale Lösung mittels Duhamel-Integral
t
X(t) = h(t – w)F(w)dw
0
Anwendung des Erwartungswertoperators

 t 
E[X(t)] = µ X (t) = E  h(t – w)F(w)dw
 
 
0 
t t
= h(t – w)E[F(w)]dw = µ F h(t – w)dw
0 0
Variablensubstitution u = t – w und anschließende Integration
µF
" q !#
µ X (t) = 1 – exp(– ζω t) ζ sin ω 0
t + 1 – ζ 2 cos ω 0 t
0
m ω 20
Erwartungswert (2)
Grenzwert
µF
lim µ X (t) = = µ X,∞
t→∞ k
2.0
1.5
X̄(t) [m]
1.0
0.5
0.0
0 20 40 60 80 100
Time [s]

Stationärer Zustand
Für ausreichend lange Anregungsdauer ist der stationäre Zustand

ausreichend um die Reaktion X(t) zu charakterisieren.
Um den stationären Zustand exakt zu erreichen, müsste die Anregung
schon seit unendliche langer Zeit einwirken.
Die Impulsreaktionsfunktion h(u) ist Null für u < 0 (Kausalitätsprinzip)
∞
→ X(t) = h(t – w)F(w)dw
–∞
Stationärer Zustand für den Mittelwert

∞ ∞
µF

E[X(t)] = µ X (t) = E  h(t – w)F(w)dw = µ F h(u)du =
 
–∞ k
0


Autokovarianz der Reaktion
Anwendung des Duhamel-Integrals für den stationären Zustand
∞ ∞ 
E[X(t)X(s)] = E  h(t – w)F(w)dw · h(s – z)F(z)dz
 
–∞ –∞

∞ ∞ 
= E  h(t – w)h(s – z)F(w)F(z)dwdz
 
–∞ –∞
 
∞ ∞
= h(t – w)h(s – z)E[F(w)F(z)]dwdz
–∞ –∞
Subtraktion von µ F und µ X an den entsprechenden Stellen führt zu
∞ ∞
RXX (t, s) = h(t – w)h(s – z)RFF (w, z)dwdz
–∞ –∞

PSD der Reaktion (1)
Nach Definition
∞ ∞ ∞
1
SXX ( ω ) = h(t – w)h(t + τ – z)RFF (z – w)eiωτ dwdzdτ
2π
–∞ –∞ –∞
Variablensubstitution u1 = z – w; u2 = t – w; u3 = t + τ – z mit der

Jacobi-Determinante
1 –1 0
∂ (u1 , u2 , u3 )
|J| = = 0 –1 0 = | – 1| = 1
∂ (z, w, τ)
–1 0 1

Umgeformt
∞
1
SXX ( ω ) = RFF (u1 )eiωu1 du1 ·
2π
–∞
∞ ∞
–i ω u2
· h(u2 )e du2 · h(u3 )eiωu3 du3
–∞ –∞
Die erste Zeile ist die PSD SFF ( ω ) der Anregung, die anderen Integrale sind
die komplexe Übertragungsfunktion H( ω ) und die dazu konjugiert
komplexe Funktion H∗ ( ω )
∞
H( ω ) = h(u)e–iωu du
–∞
Schließlich erhält man
SXX ( ω ) = SFF ( ω )H( ω )H∗ ( ω ) = SFF ( ω )|H( ω )|2

Für das SDOF-System ist dabei
1
H( ω ) =
k – m ω 2 + ıc ω
PSD der Reaktion
1
SXX ( ω ) = SFF ( ω ) 2
k + (c – 2km) ω 2 + m2 ω 4
2
Varianz
∞
σX2 = SFF ( ω )|H( ω )|2 d ω
–∞

Beispiel - Kragträger unter horizontaler Belastung
X(t)
Eigenschaften: H = 4 m, EI = 3600 kN/m2 F(t) m
und m = 1 t
Horizontalsteifigkeit k = 3EI
H3
= 400kN/m
Mittelwert der Last F̄
PSD der Last
EI H
σF2 a
SFF ( ω ) =
π(a2 + ω 2 )
mit σF = 0.2F̄ und a = 12 rad/s

Lösung
Mittelwert der Reaktion
µF
µX = = 0.0025 µ F
k
PSD von Last und Reaktion

1.5 2.0
1.5
SXX[·10–7]
1.0
SFF[·10–3]
2
`H( ω )` 1.0
0.5
0.5
0.0 0.0
0 10 20 30 40 50 0 10 20 30 40 50
ω ω
Varianz der Verschiebungsreaktion

13
σX2 = µ 2 = 2.338 · 10–6 µ 2F ; → σX = 1.529 · 10–3 µ F = 0.61 µX
5560000 F
Näherung durch weißes Rauschen
Wegen der ausgeprägten Spitze der Übertragungsfunktion hat die PSD der
Reaktion große Werte vor allem in der Nähe der Eigenfrequenz des
SDOF-Systems
Varianz kann daher approximiert wrden
∞ ∞
σx2 ≈ 2
SFF ( ω 0 )|H( ω )| d ω = SFF ( ω 0 ) |H( ω )|2 d ω
–∞ –∞
SFF( ω )
SXX( ω )
`H( ω )`2
SFF( ω 0)
ω ω
ω0 ω0

Analytische Lösung
Integration, z.B. mittels Residuensatz
∞ ∞
2 1 π
|H( ω )| d ω = dω =
k2 + (c2 – 2km) ω 2 + m2 ω 4 kc
–∞ –∞
Konstante PSD impliziert Dirac-Delta für die Autokovarianzfunktion
SFF = S0 = const. → RFF = 2πS0 δ (τ)
∞
δ(τ) = 0 : τ , 0; δ(τ)dτ = 1
–∞

Vorlesung am 5. Juni 2018
Pushover Analyse
Nichtlineare dyamische Analyse

Laststeigerung
Inkrementell gesteigerte statische Horizontallast repräsentiert
Erdbebenbelastung bis zur Kapazitätsgrenze (Bildung einer kinematischen
Kette bzw. Erreichen der statischen Stabilitätsgrenze zufolge Eigengewicht)

Last-Verformungs-Kurve
Kurve basierend auf Fließgelenktheorie

Signifikanter Einfluss des Eigengewichts (Stabilität)
Horizontal load
with dead load (P-∆)

w/o dead load
Horizontal deflection

Einfacher Rahmen - Modelldefinition
Geometrie
1.50 x 0.25
2x8
0.30 x 0.30
10
Material Stützen: linear-elastisch/ideal-plastisch mit 1% kinematischer

Verfestigung (E = 30 GPa, ν = 0.15, β y = 20 MPa, ρ = 2500 kg/m3 .
Material Riegel: linear-elastisch, E = 30 GPa, ν = 0.16, ρ = 2500 kg/m3 .

1 - -[[
2 This function defines a 2 - story frame structure
4 - -]]
5
6 -- define the frame in terms of a function
7 function frame1 ()
8 structure = fem.Structure ( " 2 story frame " )
9 blue = tmath.Matrix ({{0 , 100 , 170 , 255}}) -- define blue color
10 orange = tmath.Matrix ({{255 , 100 , 0 , 255}}) -- define orange color
11
12 s_column = structure : AddSection (101 , " RECT " , 0) ; -- define column cross section with ID 101
13 s_column : SetData ( tmath.Matrix ({{0 .3 , 0 .3 }}) )
14 s_column : SetColor ( blue )
15
16 s_beam = structure : AddSection (102 , " RECT " , 0) ; -- define beam cross section with ID 102
17 s_beam : SetData ( tmath.Matrix ({{0 .2 , 1 .5 }}) )
18 s_beam : SetColor ( orange )
19
20 m_column = structure : AddMaterial (201 , " E LA ST IC _PL AS TI C_ 1D " )
21 m_column : SetData ( tmath.Matrix ({{3 e10 , .16 , 2500 , 2 e7 , 0 .01 }}) )
22 m_beam = structure : AddMaterial (202 , " LINEAR_ELASTIC " )
23 m_beam : SetData ( tmath.Matrix ({{3 e10 , .16 , 200}}) )
24
25 -- nodes
26 B = 10; H = 8 -- width and height of one story
27 nodes = tmath.Matrix ( { -- define nodes
28 {1001 , 0 , 0 , 0} ,
29 {1002 , B , 0 , 0} ,
30 {1003 , 0 , 0 , H } ,
31 {1004 , B , 0 , H } ,
32 {1005 , 0 , 0 , 2* H } ,
33 {1006 , B , 0 , 2* H } ,
34 {1100 , B /2 , 0 , 0} -- This is the reference node for all beams
35 })
36 structure : AddNodes ( nodes ) -- Add all nodes

37
38 -- elements
39 columns = tmath.Matrix ({
40 {2001 , 1001 , 1003 , 1100} ,
41 {2002 , 1002 , 1004 , 1100} ,
42 {2003 , 1003 , 1005 , 1100} ,
43 {2004 , 1004 , 1006 , 1100}
44 })
45 structure : AddElements ( " RECT " , 201 , 101 , columns ) -- Add column elements
46
47 beams = tmath.Matrix ({
48 {3001 , 1003 , 1004 , 1100} ,
49 {3002 , 1005 , 1006 , 1100}
50 })
51 structure : AddElements ( " RECT " , 202 , 102 , beams ) -- Add beam elements
52
53 -- support conditions and global degrees of freedom
54 structure : SetAvailDof ( tmath.Matrix ({{1 , 0 , 1 , 0 , 1 , 0}}) ) -- plane problem (x - z )
56 n1001 : SetAvailDof ( tmath.Matrix ({{0 , 0 , 0 , 0 , 0 , 0}}) ) -- Clamped at node 1001
58 n1002 : SetAvailDof ( tmath.Matrix ({{0 , 0 , 0 , 0 , 1 , 0}}) ) -- Hinge at node 1002
60 n1100 : SetAvailDof ( tmath.Matrix ({{0 , 0 , 0 , 0 , 0 , 0}}) ) -- No dof at reference node
61
62 -- Assemble all degrees of freedom
63 ndof = structure : GlobalDof ()
64 return structure
65 end

Lineare Modalanalyse
Antwortspektrenmethode
1. Eigenschwingungsform
Bodenbeschleunigung: El Centro 1940 (NS)

4
Bodenbeschleunigung a(t) [m/s2]

3
2
ω 1 = 7.37 rad/s, T1 = 0.85 s 1

0
2. Eigenschwingungsform -1
-2
-3
-4
0 6 12 18 24 30
Zeit t [s]
ω 2 = 29.33 rad/s, T2 = 0.21 s

Modalanalyse - slangTNG code (1)
1 - -[[
2 Modal analysis of simple linear structures
3 ( c ) 2012 - 2017 Christian Bucher , TU Wien
5 - -]]
6
9
13 K = struct : SparseStiffness ()
14 M = struct : SparseMass ()
15
17 val , vec = K : Eigen (M , 8)
18 freq = val : CW () : Sqrt () /2/ math.pi
19 print ( " freq " , freq )
20 -- Show first 2 mode shapes
21 factor = 200
22 g ={}
23 for i =0 ,1 do
24 shape = vec : GetCols ( i )
25 struct : Se t D of D i sp l a ce m e nt s ( shape )
26 g [ i ] = graph.Graph3D ( " Mode " ..i )
27 g [ i ]: Rotate ( -70 , 1 , 0 , 0)
28 plot = structure : Draw ( factor )
29 g [ i ]: Triangles ( plot )
30 g [ i ]: Autoscale ()
31 g [ i ]: Render ()
32 end
33
34 -- Read ground acceleration data
35 Earthquake = tmath .Matri xInput ( " elcentro.txt " ) ;
36 dt = Earthquake [1] - Earthquake [0]

Modalanalyse - slangTNG code (2)
37 N = Earthquake : Rows ()
38 quake = Earthquake : GetCols (1)
39
40 -- Influence vector for ground acceleration in x - direction
41 F1 = struct : G et A l lD i s pl a c em e n ts ()
43 F1 : SetZero ()
46 -- Bring load case into vector form and multiply unit acceleration with mass matrix
47 Unit = struct : T oD of Di sp la ce me nt s ( F1 )
48 P = M * Unit
49 print ( " P " , P )
50
51 p_modal = vec : Transpose () * P
52
53 -- Evaluate response spectrum for first 2 modes with 5 % damping
54 for k =0 , 1 do
55 T = 1/ freq [ k ]
56 zeta = 0 .05
57 om2 = 4* math.pi ^2/ T ^2
58 phi = vec : GetCols ( k )
59 avd = s p e c t r a l . R e s p o n s e S p e c t r u m ( quake , dt , T , zeta )
60 Sa = avd [0]
61 Sd = avd [2]
62 f_modal = phi : Transpose () * P
63 print ( " Mode " , k , " omega " , math.sqrt ( om2 ) , " T " , T , " f_modal " , f_modal , " Sa " , Sa , " Sd " , Sd )
64 x_modal = phi * phi : Transpose () * P * Sa / om2
65 print ( " x_modal " , x_modal )
66 end

Nichtlineare statische Laststeigerung
Belastung proportional zur 1. Eigenschwingungsform
Kurve basiert auf ausgedehnten plastischen Zonen in den Finiten
Elementen
100
75
Base Shear [kN]
50 with P – ∆
w/o P – ∆
25
0
0 0.5 1 1.5
Top Displacement [m]
Plastische Kapazität wird durch Eigengewicht (Stabilität) signifikant

reduziert, Grenze wird durch Divergenz der Newton-Raphson-Iteration
angezeigt
slangTNG code (1)
1 - -[[
2 Static analysis of simple linear structures
3 ( c ) 2012 - 2017 Christian Bucher , TU Wien
5 - -]]
6
9
10 -- Prepare graphics
11 g = graph.Graph3D ( " Frame " )
12 g : Rotate ( -70 , 1 , 0 , 0)
13
14 v2 = graph.Graph ( " Response " , " Bright " )
15 v2 : AxisLabels ( " Displacement [ m ] " , " Horizontal Load [ N ] " )
16
17 for geo = 0 , 1 do
18 umax = 30000
19
24 M = struct : SparseMass ()
25
27 val , vec = K : Eigen (M , 4)
28 freq = val : CW () : Sqrt () /2/ math.pi
29 print ( " freq " , freq )
30
31 -- Influence vector gravity load
34 F1 : SetZero ()

slangTNG code (2)
37
38 -- Bring load case into vector form and multiply unit acceleration with mass matrix
39 F = struct : T oD ofD is pl ac em en ts ( F1 )
40 G = M * F *( -9 .81 ) *5
41 UG = K : Solve ( G )
42 print ( " UG " , UG )
43 struct : S et D o fD i s pl a c em e n ts ( UG )
45 KG = struct : S pa rs eG eo St if fn es s ()
46 beig , bvec = K : EigenBuckling ( KG , 4 , -1)
47 print ( " beig " , beig )
48 kg = KG : Expand ()
49 struct : Se t D of D i sp l a ce m e n ts ( UG *0)
50 -- K = K : Add ( KG , 1)
51 end
52 -- Load vector according to mode shape k
53 k = 1
54 shape = vec : GetCols (k -1)
55 shape_all = struct : To Al lD is pl ac em en ts ( shape )
56 -- remove everything except x - forces
57 xforce = shape_all : GetCols (0)
58 shape_all : SetZero ()
59 shape_all : SetCols ( xforce , 0)
60 print ( " shape_all " , shape_all )
61 F = struct : T oD ofD is pl ac em en ts ( shape_all )
62 column , row , fmax = F : CW () : Abs () : MaxCoeff ()
63 F0 = F / F [{ row , column }]
64 print ( " F0 " , F0 )
65
66 -- Number of load steps up / down
67 N = 400
68 Nm = N +1
69 du = umax / N
70
71 disp_load = tmath.Matrix ( N + Nm , 2)
72 supportlist = tmath.Matrix ({{0 , 1}})

slangTNG code (3)
73 topnode = 4
74
75 plot = struct : Draw (1)
77 g : Autoscale ()
78
79
80 disp_load [{0 ,0}] = 0
81 disp_load [{0 ,1}] = 0
82
83 F = F0 *0 -- Initial load is zero
84 U1 = struct : G et D o fD i s pl a c em e n ts ()
85
86 NN = N + Nm
87
88 for i =1 , N + Nm -1 do
89 u = i * du
90 DF = F0 * du
91 if (i > N ) then u = (2* N - i ) * du DF = DF *( -1) end
92 g : Clear ()
93 plot = struct : Draw (5) -- Deformation with scale 5
95 g : Autoscale ()
96 g : Render ()
97
98 F = F + DF
99 R1 = F *1
100 DU = K : Solve ( DF )
101 U1 = U1 + DU
102 success = true
103 for k =0 ,100 do -- Full Newton - Raphson iteration
104 struct : Se t D of D i sp l a ce m e nt s ( U1 )
106 R = R1 - struct : GlobalResForce ()
108 -- K = K : Add ( KG , -1)

slangTNG code (4)
109 R = R - KG * U1
110 end
111 Rnorm = tmath.Norm ( R )
112 if ( Rnorm < .01 * tmath.Norm ( F ) or Rnorm < 1e -3) then break end
113 DU = K : Solve ( R )
114 U1 = U1 + DU
115 if (k >98) then print ( " *** " ) success = false end
116 end
117 if ( success == false ) then NN = i ; break ; end
118 Rold = R1 *1
119 struct : S et D o fD i s pl a c em e n ts ( U1 )
120 U = struct : T oA llD is pl ac em en ts ( U1 )
121 disp_load [{ i ,0}] = U [{ topnode ,0}]
122 K = struct : SparseStiffness () ;
123 -- if ( geo ==1) then K = K : Add ( KG , -1) end
124 kk = K : Expand ()
125 struct : GlobalUpdate ()
126 -- print (" kk " , kk [0])
127
128
129 supports = struct : NodeForce ( supportlist )
130 -- print (" supports " , supports )
131 sum = 0
132 for k =0 , supports : Rows () -1 do
133 sum = sum + supports [{ k , 0}]
134 end
135 disp_load [{ i ,1}] = - sum
136 disp_load [{ i ,1}] = u
137 print ( " - sum " , - sum )
138
139 end -- for
140
141 disp_load = disp_load : GetRows (0 , NN )
142
143 legend = " w / o P - Delta "
144 if ( geo == 1) then legend = " with P - Delta " end

slangTNG code (5)
145 v2 : Plot ( disp_load : GetCols (0) , disp_load : GetCols (1) , 1 , legend )

146 tmath.Output ( disp_load , " disp_load_ " ..geo.. " .txt " )
147 -- v2 : SetRange (0 , 1 .3 , 0 , 1 .3e5 )
148 -- v2 : Render ()
149 end

Äquivalentes SDOF-System (1)
Elastisch-plastisches Verhalten des SDOF-Systems unter statischer

Belastung
Elastische Grenzlast Fy , elastische Verschiebungsgrenze xy ,
Anfangssteifigkeit k, Steifigkeit nach Fließbeginn α k
Duktilität µ = xmax
xy
F
Yield strength reduction factor R = Feqy
F
Feq
Fy
Emax
k
x
xy xeq xmax

Äquivalentes SDOF-System (2)
Maximale Verformungsenergie Emax

Energieäquivalentes lineares System
r
2Emax
Emax =
k
Maximale Kraft im äquivalenten linearen System
1 2Emax
Emax = Feq xeq → Feq =
2 xeq
Maximale Verschiebung des äquivalenten linearen Systems auf

Erdbebenanregung
m
r
xeq = mSa (T); T = 2π
k

Beispiel Rahmen
Idealisiertes SDOF-System (1)
Last-Verformungskurven sowie maximale Formänderungsenergie sollen

übereinstimmen
F
Fmax
Fy
idealized SDOF
nonlinear MDOF
(FEM)
x
xy xmax
Steifigkeit und Eigenkreisfrequenz sind daher

r
Fy k
k = ; ω0 =
xy m

Beispiel Rahmen
Idealisiertes SDOF-System (2)
ω0 folgt aus der Modalanalyse des linearen MDOF-Systems → äquivalente
Masse m kann damit berechnet werden
ω 20
m=
k
2π
Für die Schwingungsperiode T = ω 0
wird das (inelastische)
Beschleunigungsantwortspektrum Sa (T) ermittelt. Der Reduktionsfaktor R
ist dann gegeben durch
mSa (T)
R=
Fy
Zielverschiebung
xt = Rxy
Vergleich zur Maximalverschiebung:
xmax
Rxy < xmax → R <
xy

Rahmenbeispiel (1)
Grundfrequenz f0 = 1.17 Hz, T = 0.855 s
Fy = 83 kN, xy = 0.40 m → keq = 208 kN/m
meq = 3.85 t
Maximale Duktilität µ max = 1.10/0.40 = 2.75
Inelastisches Antwortspektrum für µ = 2:
Sa (T = 0.855, ζ = 0.05, µ = 2, α = 0) = 3m/s2
R = 3.85 ·3
83 = 0.14
xt = 0.14 · 0.40 = 0.06 m → Duktilität nicht erreicht.
Neuberechnung für µ = 1:
Sa (T = 0.855, ζ = 0.05, µ = 1.0, α = 0) = 5.8m/s2
R = 3.8583·5.8 = 0.27
xt = 0.27 · 0.40 = 0.11 m → bleibt elastisch.
Verschiebung an der Struktur
x = φ 1 φ T1 pxt = [0.084, 0.131]T
Zum Vergleich: Maximale Verschiebung aus dynamischer Berechnung nach

Newmark: xdyn = 0.13 m
Inelastische Antwortspektren (1)
Dargestellt als Funktion der vorgegebenen Duktilität µ

Beispiel El Centro 1940 (NS), ζ = 5%, α = 0
10
8
α=0
µ=1
µ=2
Sa [m/s2]
6
µ=3
µ=5
4
0
0 0.5 1 1.5 2 2.5 3
Period T [s]

Inelastische Antwortspektren (2)
Beispiel El Centro 1940 (NS), ζ = 5%, α = 0.5

10
8
α = 0.5
µ=1
µ=2
Sa [m/s2]
6
µ=3
µ=5
4
0
0 0.5 1 1.5 2 2.5 3
Period T [s]

Künstliche Erdbebenschriebe (1)
Für nichtlineare Zeitverlaufsberechnungen werden
Bodenbeschleunigungen entsprechend vorgegebener (linearer)
Antwortspektren benötigt.
Dies kann durch einen einfachen iterativen Prozess erfolgen.
Vorgabe einer Einhüllenden (Intensität, Bebendauer)
Simulation von Amplituden und Phasenwinkel zu Frequenzkomponenten
(Fourier-Synthese)
Modulation mit der Einhüllenden
Iterative Modifikation der Amplituden zu den einzelnen
Frequenzkomponenten so, dass das Antwortspektrum angenähert wird.
Abbruch nach wenigen (≈ 7) Iterationen

Künstliche Erdbebenschriebe (2)
Zwei ausgewählte Schriebe:
Modulating Function
Acceleration a [m/s2]
-3
-6
0 5 10 15 20 25 30
Time t [s]

Bewertung künstlicher Erdbebenschriebe
Vergleich des Mittelwerts aus 10 generierten Schrieben mit dem
vorgegebenen Antwortspektrum nach EC8 (PGA = 5 m/s2 ).
15
Response Spectrum Sa [m/s2]
Mean of 10 records
10 Target
0
0 1 2 3 4
Period T [s]

Unsymmetrischer Rahmen
Modell
Nichtlineare dynamische Berechnung

Livedemo slangTNG

Vorlesung am 19. Juni 2018
Modelle für natürlichen Wind

Reaktion in Windrichtung (böenerregte Schwingungen)
Reaktion quer zur Windrichtung (wirbelinduzierte Schwingungen)

Dynamische Reaktion auf Windlasten
Linearisierung und Anwendung der Leistungsspektralmethode

Berechung der Reaktion auf ständige Lasten und mittlere Windlast, ggf. unter
Einbeziehung geometrischer und physikalischer Nichtlinearitäten
Ermittlung der tangentiellen Steifigkeitsmatrix und Eigenschwingungsformen
mittels Linearisierung am Gleichgewichtspunkt
Annahme eines stochastischen Turbulenzfelds für die longitudinale
Komponente der Windgeschwindigkeit
Berücksichtigung der unvollständigen räumlichen Kohärenz des
Turbulenzfeldes (wirkt lastmindernd)
Anwendung der Leistungspektralmethode zur Berechnung der Varianz der
Reaktion
Zeitverlaufsberechnung
Annahme eines stochastischen Turbulenzfelds für die Windgeschwindigkeiten
Überlagerung von ständigen Lasten und Windlasten
Berechnung der Reaktion auf kombinierte Belastung, ggf. mit Hilfe
nichtlinearer dynamischer Analyse

Theoretische Grundlagen
Querschnitt in stationärer Strömung
qL
v S qD
m
Druckdifferenz im Staupunkt
ρ v2
pu =
2
Longitudinale und laterale Lasten pro Längeneinheit
ρv2 ρv2 ρv2 2
qD = BCD ; qL = BCL ; m= B CM
2 2 2
B ist eine charakteristische Längenabmessung; CD , CL und CM sind
dimensionlose Koeffizienten, die experimentell bestimmt werden
(Windkanal). Abhängig vom Form und Orientierung des Querschnitts.
Räumlich-zeitliche Schwankungen
Turbulenter Wind v(t) z
v(t) = v̄ + v 0(t)
v̄ . . . mittlere
Windgeschwindigkeit v(z)
v̄(z)
(Mittelungsintervall 10 min),
v 0(t) . . . mittelwertfreie zufällige
Schwankung
Veränderung der mittleren
Windgeschwindigkeit 10 m
!α
z
v̄(z) = v̄ref
zref

Davenport’sches Konzept (1)
Leistungsspektraldichte (PSD) von v
σv2 x2 Lω
Svv ( ω ) = ; x=
3 ω (1 + x2 )4/3 2πv̄10
σv . . . Standardabweichung der Windgeschwindigkeit

L = 1200m . . . Längenmaß der atmosphärischen Turbulenz
v̄10 . . . mittlere Windgeschwindigkeit in der Referenzhöhe von 10 m
100
PSD of wind velocity [m2/s]
10
0.1
0.01
0 10 20 30
Circular frequency ω [rad/s]

Turbulenzintensität I
σv
I=
v̄
Vertikale Kohärenzfunktion (Korrelation der Windgeschwindigkeiten in den
Höhen zm und zn ), frequenzabhängig
10|zm – zn | ω
!
R(zm , zn , ω ) = exp –
π(v̄m + v̄n )
Aerodynamischer Druck
ρv(t)2
p(t) = = p̄ + p 0(t)
2
Linearisierung (Vernachlässigung von v 02 )
1 2 0
p̄ = ρv̄ , p (t) = ρv̄v 0(t)
2

Mittlere aerodynamische Last
1
q̄D = ρBCD v̄2
2
Leistungsspektraldichte der aerodynamischen Last qd
Sqq ( ω ) = ( ρ Bv̄)2 χF2 ( ω )Svv ( ω )
Aerodynamische Admittanzfunktion
1 1
χF2 ( ω) = C2D 20 ω B
·
1+ 4πv̄ 1 + 84ωπv̄B

χF2( ω)
C2D
0.8
Normalized admittance
B = 5 m, v̄ = 25 m/s
0.6
0.4
0.2
0
0 10 20 30
Circular frequency ω [rad/s]
Kreuzleistungsspektraldichte in verschiedenen Höhen zm und zn

q
Spm pn ( ω ) = ρ 2 Bm Bn v̄m v̄n χFm χFn Svm vm ( ω )Svn vn ( ω )R(zm , zn , ω )

Modale, auf die Eigenformen φ j bezogene Belastung Fj (t)

Fj (t) = φ j (z)p(z, t)dz
L
Mittlere modale Belastung F̄j (t):

F̄j (t) = φ j (z)p̄(z)dz
L
Kreuzleistungsspektraldichte zweier modaler Lastgrößen Fj und Fk :

SFj Fk ( ω ) = φ j (x1 )φ k (x2 )Spm pn (x1 , x2 , ω )dx1 dx2
L L
2 durch Integration
Berechnung der Varianz der Verschiebungsreaktion σw
der Leistungsspektraldichte im Frequenzbereich
Erwartungswert der Reaktionsamplituden
E[wmax ] = w̄ + gσw
Spitzenfaktor(ν = ω 0 /2π)
p 0.6
g= 2 ln(νT) + √
2 ln(νT)
Böreaktionsfaktor (Gust response factor)
E[wmax ] gσw
GRF = =1+
w̄ w̄

Leistungsspektralmethode
Trennung der Varianz der Reaktion in Hintergrund- und Resonanzanteil

Varianz kann approximiert werden
∞ ∞
σ2 πSFF ( ω0 )
σx2 ≈ 2
SFF ( ω )|H(0)| d ω + SFF ( ω 0 ) |H( ω )|2 d ω = 2F +
k kc
–∞ –∞
SFF( ω )
`H( ω )`2 Resonance
SXX( ω )
SFF( ω 0)
Background
ω ω
ω0 ω0

Eurocode 1 Fragment

Übungsbeispiel - Wasserturm
Gesamtmasse aus Struktur und

Wasser (m) m
Linear-elastisches Verhalten des F(t) x(t)
Schafts, bestimmt durch Geometrie
(R, r) und Werkstoff (E)
Dämpfung muss angenommen
werden (z.B. nach Norm)
Belastung F(t) infolge Wind
H
Gesucht: Böreaktionsfaktor für T = 10 R

min
r
Annahme: Windlast wirkt nur auf den
ellipsoidförmigen Tank mit Breite B =
10 m und Höhe A = 5 m. Höhe über
dem Bode ist H = 25 m

Daten
Strukturdaten
Masse m = 400 t
Elastizitätsmodul E = 30 GPa
Querschnitt R = 1.2 m, r = 1.08 m.
Dämpfung ζ = 0.05
Wind- und aerodynamische Daten
Windgeschwindigkeit v̄10 = 25 m/s
Turbulenzintensität I = 0.15
Exponent des Höhenprofils α = 0.4
Widerstandsbeiwert CD = 1.2
Gesucht: Böreaktionsfaktor
Annahme: masseloser Schaft
Annahme: Konzentrierte Windlast in Tankmitte
Lösung: GRF = 2.92

slangTNG code (1)
1 R = 1 .2
2 r = 0 .9 * R
3 EI = 30 e9 * math.pi *( R ^4 - r ^4) /4
4 m = 400 e3
5 H = 25
6 zeta = 0 .05
7
8 k = 3* EI / H ^3
9
10 om0 = math.sqrt ( k / m )
11 print ( " om0 " , om0 )
12 c = 2* zeta * math.sqrt ( k * m )
13
14 I = 0 .15
15 alpha = 0 .4
16 v10 = 25
17 CD = 1 .2
18 B = 10
19 A = 5
20 Area = A * B * math.pi /4
21 vm = 25*( H /10) ^ alpha
22 rho = 1 .26
23 sigma_v = vm * I
24 chi2 = CD ^2/(1+20* om0 * B /4/ math.pi / vm ) /(1+8* om0 * B /4/ math.pi / vm )
25
26 print ( " vm " , vm )
27 print ( " chi2 " , chi2 )
28
29 x = 1200* om0 /2/ math.pi / v10
30 print ( " x " , x )
31 dave = sigma_v ^2/3/ om0 * x ^2/(1+ x ^2) ^(4/3)
32 print ( " dave " , dave )
33
34 Fm = 0 .5 * rho * CD * vm ^2* Area
35 sigma_F = rho * CD * vm * sigma_v * Area
36 print ( " Fm " , Fm , " sigma_F " , sigma_F )

slangTNG code (2)
37 wm = Fm / k
38 back = sigma_F ^2/ k ^2
39 reso = 4* rho ^2* vm ^2* CD ^2* chi2 * dave * Area ^2* math.pi / c / k
40 print ( " back " , back , " reso " , reso )
41 sigma_w = math.sqrt ( back + reso )
42 print ( " wm " , wm , " sigma_w " , sigma_w )
43
44 nuT = 10*60* om0 /2/ math.pi
45 g = math.sqrt (2* math.log ( nuT ) ) +0 .6 / math.sqrt (2* math.log ( nuT ) )
46 GRF = 1 + g * sigma_w / wm
47 print ( " g " , g , " GRF " , GRF )

Zeitverlaufsberechnung
Sehr hohes Gebäude, 800 m Höhe
Sehr schlank und sehr niedrige Eigenfrequenzen (ca. 0.11 Hz)
Wesentlicher Einfluss unvollständiger Kohärenz über die Höhe

Querschwingungen
Für einen Zylinder mit Durchmesser d

PSD der lateralen Windlast
% 2
C2L v̄10 d2
2 ( ω – ω S )2 
 
Sqq ( ω ) = √ exp –
2πIv ω S  2I2v ω2S 
Wirbelablösefrequenz (Kármán’sche Wirbel)

2π v̄10
ωS = Sr
d
CL . . . Auftriebsbeiwert
Sr . . . Strouhalzahl
Vertikale Kohärenz der lateralen Windlasten in Höhen zi und zk
2|zi – zk | |zi – zk |
! " #
γq q = cos exp –
i k 3d 3d

Baudynamik 18

Hochgeladen von

Copyright:

Verfügbare Formate

Baudynamik 18

Hochgeladen von

Dokumentinformationen

Originaltitel

Copyright

Verfügbare Formate

Dieses Dokument teilen

Dokument teilen oder einbetten

Freigabeoptionen

Stufen Sie dieses Dokument als nützlich ein?

Sind diese Inhalte unangemessen?

Copyright:

Verfügbare Formate

Baudynamik 18

Hochgeladen von

Copyright:

Verfügbare Formate

Baudynamik

Forschungsbereich für Baumechanik und Baudynamik

Letzte Überarbeitung: 9. März 2018

2/128 © Christian Bucher 2007-2018

m1 0 ẍ1 c1 + c2 –c2 ẋ1

3/128 © Christian Bucher 2007-2018

Mẍ + Cẋ + Kx = f(t); x(0) = x0 ; ẋ(0) = v0 (1)

4/128 © Christian Bucher 2007-2018

Elemente sind Kräfte infolge kinematischer Einheitsgrößen

5/128 © Christian Bucher 2007-2018

Ziel ist dieReduktion der Bewegungsgleichungen (System gekoppelter

ΦT MΦ = I; ΦT KΦ = diag(ω 2k ); k = 1...n (2)

Die Spalten φ k der Modalmatrix Φ sind Lösungsvektoren des

Eine nicht-triviale Lösung kann nur dann existieren, wenn

6/128 © Christian Bucher 2007-2018

ΦT MΦÿ + ΦT KΦy = ΦT f(t) = g(t)

Die Gleichungen für die neuen Variablen y sind entkoppelt

ÿk + ω 2k yk = gk (t); k = 1...n (6)

Die zugehörigen Anfangsbedingungen ergeben sich durch Inversion von (5)

→ y0 = ΦT Mx0 ; ẏ0 = ΦT Mv0 (8)

Die charakteristische Gleichung (4) ist

Entwicklung der Determinante

9/128 © Christian Bucher 2007-2018

10/128 © Christian Bucher 2007-2018

11/128 © Christian Bucher 2007-2018

Mode 1 (Biegung/Schub in y) Mode 2 (Biegung/Schub in x)

Mode 3 (Torsion um z) Mode 4 (Biegung/Schub in y)

13/128 © Christian Bucher 2007-2018

14/128 © Christian Bucher 2007-2018

ÿk + 2 ζ k ω k ẏk + ω 2k yk = gk (t); k = 1...n (5)

ΦT CΦ = diag(2ζk ω k ); k = 1...n (6)

15/128 © Christian Bucher 2007-2018

16/128 © Christian Bucher 2007-2018

Für gegebene Anfangsbedingungen x(0) = x0 und ẋ(0) = v0

y0 = ΦT Mx0 ; ẏ0 = ΦT Mv0

17/128 © Christian Bucher 2007-2018

Erzwungene Schwingung mittels Modalanalyse

18/128 © Christian Bucher 2007-2018

Allgemeine geschlossene Lösung: Duhamel-Integral X

Annahme: Belastung ist in einem Zeitintervall ∆t konstant, Ansatz einer

19/128 © Christian Bucher 2007-2018

yk,h (t) = e–ζk ωk t A cos ω k0 t + B sin ω k0 t

Ermittlung von A und B aus Anfangs- bzw. Übergangsbedingungen

Lösung am Ende des Zeitintervalls

Numerisch eﬃzienter als Duhamel-Integral

20/128 © Christian Bucher 2007-2018

Vereinfachende Annahmen für x(t) innerhalb eines (kleinen) Zeitintervalls

21/128 © Christian Bucher 2007-2018

Verschiebung und Geschwindigkeit zum Zeitpunkt t + ∆t

x(t + ∆t) = x(t) + ẋ(t)∆t

Die Beschleunigung ẍ(t) wird aus der Bewegungsgleichung zum Zeitpunkt t

22/128 © Christian Bucher 2007-2018

Freie ungedämpfte Schwingung, Masse m = 1 kg, Steiﬁgkeit k = 1 N/m2 ,

Energiefehler wächst exponential mit der Zeit

23/128 © Christian Bucher 2007-2018

Geschwindigkeit und Beschleunigung zum Zeitpunkt t

x(t + ∆t) – x(t – ∆t)

24/128 © Christian Bucher 2007-2018

Die Beschleunigung ẍ(0) wird aus (1) für t = 0 ermittelt: