结构动力学

m ∆t2
+
c 2∆t
⎞⎠⎟ui+1
=
pi
−
⎛ ⎝⎜
k
−
2m ∆t2
⎞⎠⎟ui
−
⎛ ⎝⎜
m ∆t2
−
c 2∆t
⎟⎞⎠ui−1
中心差分法的数值稳定性证明
设体系为无阻尼，并设外荷载p=0 (算法的稳定性与外荷载无关),则 中心差逐步积分法的递推公式可以写成如下形式：
u&0
=
u1 − u−1 2∆t
u&&0
=
u1
−
2u0 + u−1 ∆t 2
u−1
=
u0
−
∆tu&0
+
∆t 2 2
u&&0
u&&0
=
1 m
(
p0
−
cu&0
−
ku0 )
⎛
5中.3心中差心分差法分计法算步骤：⎜⎝
m ∆t 2
+
c 2∆t
⎞ ⎟⎠
ui
+1
=
pi
−
⎛ ⎜⎝
k
−
2m ∆t 2
⎞ ⎟⎠
ui
−
⎛ ⎜⎝
u&(τ ) = A1 + (ωD A3 − ζωn A2 )e−ζωnτ cosωDτ − (ωD A2 + ζωn A3 )e−ζωnτ sin ωDτ
其中，
A0
=
pi k
− 2ζαi kωn
,
A1
=
αi
k
,
A2 = ui − A0,
A3
=
1
ωD
[u&i
+
ζωn
A2
−
αi
k
]
5.2 分段解析法
{u&}i = {u&(ti )}
{u&&}i = {u&&(ti )}
{ p} i
=
{p(ti )}
5.3 中心差分法
⎛ ⎜⎝
m ∆t 2
+
c 2∆t
⎞ ⎟⎠
ui
+1
=
pi
−
⎛ ⎜⎝
k
−
2m ∆t 2
⎞ ⎟⎠
ui
−
⎛ ⎜⎝
m ∆t 2
−
c 2∆t
⎞ ⎟⎠
ui
−1
单步法和多步法的概念
ui+1 = Aui + Bu&i + Cpi + Dpi+1 u&i+1 = A′ui + B′u&i + C′pi + D′pi+1
p
实际荷载 pi+1
pi 插值荷载：p(τ)
当τ=Δti时，得到
ui+1 = Aui + Bu&i + Cpi + Dpi+1 u&i+1 = A′ui + B′u&i + C′pi + D′pi+1
△ti
ti
ti+1
t
τ
其中系数A—D′是结构刚度k，自振频率ωn，阻尼比ζ
和时间步长Δt的函数。
上式给出了根据i时刻运动及外力计算i+1时刻运动的递推
5.2 分段解析法（Piecewise Exact Method）
分段解析算法假设 p 在ti≤t≤ti+1时段内
p(τ ) = pi + αiτ αi = ( pi+1 − pi ) / ∆ti
如果荷载p(t)采用 计算机采样，即 离散数值采样， 则以上定义可 认为是“精确”的。
实际荷载 pi
pi+1 插值荷载：p(τ)
△ti
ti
ti+1
t
τ
分段解析法对外荷载的离散
5.2 分段解析法
p
实际荷载 pi+1
pi 插值荷载：p(τ)
△ti
在ti≤t≤ti+1时段内体系的运动方程：
ti
ti+1
t
τ
mu&&(τ ) + cu&(τ ) + ku(τ ) = p(τ ) = pi + α iτ
初值条件：
u(τ ) τ =0 = ui , u&(τ ) τ =0 = u&i
pˆi = pi − aui − bui−1
u&i
=
ui+1 − ui−1 2∆t
ui+1 = pˆi / kˆ
u&&i
=
ui +1
− 2ui ∆t 2
+ ui−1
(4). 下一步计算用i+1代替i，重复(2)至(3)中的计算步骤。
5.3 中心差分法
中心差分法的精度和数值稳定性
以上给出的中心差分逐步积分公式， ① 是收敛的；
单步法：采用时域逐步积分法计算某一时刻的运动时，
仅需已知前一时刻的运动。
多步法：需要前两个或两个以上时刻的运动。
中心差分法在计算ti+1时刻的运动ui+1时，需要已知ti和ti-1两个 时刻的运动ui和ui-1，因此，中心差分法属于两步法；而分 段解析法为单步法。
用两步法进行计算时存在起步问题，因为仅根据已知的初始 位移和速度，并不能自动进行运算，而必需给出两个相邻 时刻的位移值，方可开始逐步计算。
而这种离散化正符合计算机存贮的特点。
与运动变量的离散化相对应，体系的运动微分方程也不 一定要求在全部时间上都满足，而仅要求在离散时间 点上满足，这相当于放松了对运动变量的约束。
5.1 数值算法中的基本问题
离散的定义？
采用等时间步长离散时，ti=iΔt，i=1, 2, 3,…。
体系的运动微分方程仅要求在离散时间点上满足。 Δt——离散时间步长
运动方程的特解:
u p (τ )
=
1 k
( pi
+ αiτ )
−
αi
k2
c
运动方程的通解:
uc (τ ) = e−ζωnτ ( AcosωDτ + B sin ωDτ )
5.2 分段解析法
将全解
u(τ ) = up (τ ) + uc (τ )
代入边界条件确定系数A、B，最后得：
u(τ ) = A0 + A1τ + A2e−ζωnτ cosωDτ + A3e−ζωnτ sin ωDτ
m ∆t 2
−
c 2∆t
⎞ ⎟⎠
ui
−1
(1). 基本数据准备和初始计算
u0 和 u&0 已知
u−1
=
u0
−
∆tu&0
+
∆t 2 2m
( p0
−
cu&0
−
ku0 )
(2). 计算等效刚度和中心差分计算公式中的系数
kˆ
=
m ∆t 2
+
c 2∆t
,
a
=
k
−
2m ∆t 2
,
b
=
m ∆t 2
−
c 2∆t
(3). 根据i及i以前时刻的运动，计算i+1时刻的运动
结构动力反应分析的时域直接数值计算方法：
（1）分段解析法； （2）中心差分法； （3）平均常加速度法； （4）线性加速度法； （5）Newmark－β法； （6）Wilson－θ法。
•••••••••
时域逐步积分法是结构动力分析问题中一个得到广泛研 究的课题，也是得到广泛应用的计算方法。
5.1 数值算法中的基本问题
⎜⎛ ⎜
⎝
ωn 1−ζ
2
sin
ω
D
∆t
⎟⎞ ⎟
⎠
B′
=
e
−ζωn∆t
⎜⎛ ⎜
cos
ω
D
∆t
−
⎝
ζ 1−ζ
2
sin
ω
D
∆t
⎟⎞ ⎟
⎠
C′
=
1 k
⎪⎨⎧− ⎪⎩
1 ∆t
+
e −ζωn∆t
⎢⎡⎜⎛ ⎢⎣⎜⎝
ωn + 1 − ζ 2 ∆t
ζ 1−ζ
2
⎟⎞ ⎟
sin
ω
D
∆t
⎠
+
1 ∆t
cos ω D ∆t ⎥⎥⎦⎤ ⎪⎭⎪⎬⎫
m ui+1
− 2ui ∆t 2
+ ui−1
+ c ui+1 − ui−1 2∆t
+ kui
=
pi
ui = u(ti ) u&i = u&(ti ) u&&i = u&&(ti ) pi = p(ti )
⎛ ⎜⎝
m ∆t 2
+
c 2∆t
⎞ ⎟⎠
ui
+1
=
pi
−
⎛ ⎜⎝
k
−
2m ∆t 2
⎞ ⎟⎠
ui
−
5.3 中心差分法（Central Difference Method）
中心差分方法用有限差分代替位移对时间的求导（即速 度和加速度）。如果采用等步长，Δti=Δt，则i时刻 速度和加速度的中心差分近似为：
u&i
=
ui+1 − ui−1 2∆t
u&&i
=
ui+1
− 2ui ∆t 2
+
ui−1
mu&&(ti ) + cu&(ti ) + ku(ti ) = p(ti )
计算公式。如果结构是线性的，并采用等时间步长，
则A—D′均为常数，其计算效率非常高，在p(t)为离
散采样的定义下是精确解，但如果是非线性问题，则
A—D′均为变量，计算效率会大为降低。
5.2 分段解析法
A
=
e −ζωn∆t
⎜⎛ ⎜
⎝
ζ 1−ζ
2
sin
ω
D
∆t
+
cosω
D
∆t
⎟⎞ ⎟
⎠
分段解析法 计算公式中
¾ 计算效率：所花费的计算时间的多少。
一个好的方法首先必须是收敛的、有足够的精度（例如2 阶，满足工程要求）、良好的稳定性、较高的计算效 率。在发展逐步积分法中，也的确发展了一些高精度 但很费时的方法，在实际中得不到应用和推广。
5.1 数值算法中的基本问题
根据是否需要联立求解耦联方程组，逐步积分法可分为 两大类：
1 k
⎡ ⎢1 − ⎢⎣
2ζ ω n ∆t
+
e−ζωn∆t ⎜⎜⎝⎛
2ζ ω
2 D
−1 ∆t
sin
ω
D
∆t
+
2ζ ω n ∆t
cosωD∆t ⎟⎟⎠⎞⎥⎥⎦⎤
ui+1 = Aui + Bu&i + Cpi + Dpi+1
u&i+1 = A′ui + B′u&i + C′pi + D′pi+1
A′
=
−e −ζωn∆t
采用叠加原理的时域和频域分析方法（Duhamel积分， Fourier变换），假设结构在全部反应过程中都是线性 的，而时域逐步积分法，只假设在一个时间步距内是 线性的，相当于用分段直线来逼近实际的曲线。
时域逐步积分法研究的是离散时间点上的值，例如位移 和速度为：
ui = u(ti ) , u&i = u&(ti ) , i = 1, 2, L
在初始时刻需要建立两个起步时刻（即i＝0, -1）的位移值， 这即是逐步积分的起步问题。
5.3 中心差分法
中心差分方法计算中的起步处理方法
⎛ ⎜⎝
m ∆t 2
+
c 2∆t
⎞ ⎟⎠
ui
+1
=
pi
−
⎛ ⎜⎝
k
−
2m ∆t 2
⎞ ⎟⎠
ui
−
⎛ ⎜⎝
m ∆t 2
−
c 2∆t
⎞ ⎟⎠
ui
−1
初始条件为 ： u0 = u(0), u&0 = u&(0)
隐式方法：逐步积分计算公式是耦联的方程组，需联立 求解，计算工作量大，通常增加的工作量与自由度的 平方成正比，例如Newmark—β法、Wilson —θ法。
显式方法：逐步积分计算公式是解耦的方程组，无需联 立求解，计算工作量小，增加的工作量与自由度成线 性关系，如中心差分方法（无阻尼时）。
下面先介绍分段解析算法，然后重点介绍两种常用的时 域逐步积分法—中心差分法和Newmark—β法，同时 也介绍Wilson —θ法，最后介绍非线性问题分析方法。
② 具有2阶精度，即误差 ε∝O(Δt2) ；
③ 是有条件稳定，稳定条件Δt≤Tn/π； ④ 具有较高的计算效率。
5.3 中心差分法 中心差分法的数值稳定性 (Δt≤Tn/π)
稳定性的含义，当满足稳定性条件时，计算值u为有限值； 当不满足稳定性条件时，随着t→∞，u→∞。
5.3 中心差分法
⎛ ⎜⎝
=
{ p} i
−
⎛ ⎜⎝
[
K
]
−
2 ∆t 2
[M
]⎟⎞⎠{u}i
−
⎛ ⎜⎝
1 ∆t 2
[M
]
−
1 2∆t
[C ] ⎟⎞⎠{u}i −1
( ) {u&}i
=1 2∆t
{u} i
+1
−
{u} i
−1
( ) {u&&}i
=
1 ∆t 2
{u} i
+1
−
2{u} i
+
{u} i
−1
{u} i
=
{u(ti )}
这两种分析方法的特点是均基于叠加原理，要求结构体 系是线弹性的，当外荷载较大时，结构反应可能进入 物理非线性（弹塑性），或结构位移较大时，结构可 能进入几何非线性，这时叠加原理将不再适用。此时 可以采用时域逐步积分法求解运动微分方程。
5.1 数值算法中的基本问题
时域逐步积分法——Step-by-step methods
D′
=
1 k∆t
⎡ ⎢1 ⎢
−
e
−ζωn
∆t
⎜⎛ ⎜
⎣
⎝
ζ 1−ζ
2
sin
ω
D
∆t
+
cos
ω
D
∆t
⎟⎞⎥⎤ ⎟⎥
⎠⎦
5.2 分段解析法
分段解析法的误差仅来自对外荷载的假设， 而 在连续时间轴上严格满足运动微分方程。
一般的时域逐步积分法将进一步放松要求，仅要 求在离散的时间点上满足运动方程，即放松了 对运动的约束。
结构动力学
清华大学土木工程系 刘晶波
2005年秋
结构动力学
第5章 结构动力反应 数值分析方法
5.1 数值算法中的基本问题