结构动力学使用中心差分法计算单自由度体系动力反应的MATLAB程序

中心差分法计算单自由度体系动力反映的报告

前言

基于叠加原理的时域积分法与频域积分法一样，都假设结构在在全部反应过程中都是线性的。而时域逐步积分法只是假设结构本构关系在一个微小的时间步距内是线性的，相当于分段直线来逼近实际的曲线。时域逐步积分法是结构动力问题中研究并应用广泛的课题。中心差分法是一种目前发展的一系列结构动力反应分析的时域逐步积分法的一种，时域逐步积分法还包括分段解析法、平均常加速度法、线性加速度法、Newmarket−β和Wilson−θ法等。

中心差分法(central difference method)原理

中心差分法的基本思路

将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。

中心差分法是一种显示的积分法，它基于用有限差分代替位移对时间的求导（即速度和加速度）。如果采用等时间步长，∆t i=∆t（∆t为常数），则速度与加速度的中心差分近似为

u i=u i+1+u i−1

2∆t

(1)

üi=u i+1−2u i+u i−1

∆t2

(2)

用u表示位移，离散时间点的运动为：

u i=u(t i)，u i=u̇(t i)，u i=ü(t i)(i=0,1,2…)

体系的运动方程为

mü(t)+cu̇(t)+ku(t)=P(t)(3)

将速度和加速度的差分近似公式(1)和(2)代入(3)中得出在t i时刻的运动方程，将方程整理得到u i+1由u i 和u i−1表示的两步法的运动方程(4)：

(m ∆t2+c

2∆t

)u i+1=P i−(k−2m

∆t2

)u i−(m

∆t2

−c

2∆t

)u i−1(4)

这样就可以根据t i及以前的时刻的运动求得t i+1时刻的运动。

中心差分法属于两步法，用两步法计算时存在起步问题，必须要给出相邻的两个时刻的位移值，才能逐步计算。对于地震作用下结构的反应问题和一般的零初始条件下的动力问题，可以用(4)直接计算，因为总可以假设初始的两个时间点（一般取i=0,−1）的位移等于零。但是对应于非零初始条件或零时刻外荷载很大时，需要进行一定的分析，建立两个起步时刻（即i=0,−1）的位移值。

假设给定的初始条件为

u0=u(0)

u̇0=u̇(0)

}(5)

根据初始条件来确定u−1。根据中心差分公式

u̇0=u1+u−1

2∆t

ü0=u1−2u0+u−1

∆t2

}(6)

消去u1得到u−1的公式：

u−1=u0−∆tu̇0+∆t2

2

ü0(7)

其中零时刻加速度值ü0可以由t =0时的运动方程得到即

ü0=1m (P 0−cu̇0−ku 0) (8)

这样就可以根据初始条件得到u −1，然后再将初始条件应用于公式(4)中，逐步求出不同时刻的运动。 中心差分法分析时的具体计算步骤：

(1) 基本数据准备与初始条件计算

已知：初始位移u 0、u̇0和初始荷载值P 0来计算ü0和u −1

ü0=1m (P 0

−cu̇0−ku 0) u −1=u 0−∆tu̇0+∆t 22ü0

(2) 计算等效刚度和中心差分法计算公式中的系数

k ̂=m ∆t 2+c 2∆t

a =k −

2m ∆t 2

b =m ∆t 2−

c 2∆t 因此中心差分法计算公式可以表示为：k

̂u i+1=P i −au i −bu i−1 (3) 根据t i 及以前的时刻的运动求得t i+1时刻的运动

P

̂i =P i −au i −bu i−1 u i+1=P

̂i k ̂⁄ (4)下一步计算中用i +1代替i ，对于线弹性体系重复第3步计算步骤，对于非线性弹性体系，重复第2和第3计算步骤。

以上的中心差分法逐步计算公式具有2阶精度，即误差ϵ∝O(∆t 2)；并且是有条件稳定的，稳定条件为：

∆t ≤

T n π

式中，T n 为结构的自振周期，对于多自由度体系则为结构的最小自振周期。 算例

本算例根据结构动力学48页算例3.1数据编写，稳定条件为dt <=0.16s

对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M =9240kg 、刚度K ＝1460KN/m 、阻尼系数C ＝6.41KN ∙s/m ，对结构施加动力荷载P =73000∙sin 0.5πt 假设结构处于线弹性状态，用中心差分法计算结构的自由振动反应。

采用MATLAB 语言编程，并以单自由度体系为例进行计算,设初位移u0＝0.05m 和初速度v0=0，取不同的步长分别计算，以验证中心差分法的稳定条件πn

T t ≤∆。

先计算t ∆，由稳定条件n n

T t ωπ2

=≤∆，而12.57n ω=

==rad/s ，

则20.16n n

T t s πω∆≤==，所以本次计算取t ∆＝0.2，0.1，0.05分别进行计算。 计算结果与分析

1）当dt =0.05s 时，可以得到位移u ，速度v ，加速度ac 的时程曲线如下：

2）当dt =01s 时，可以得到位移u ，速度v ，加速度ac 的时程曲线如下：

3）当dt =0.2s 时，可以得到如下提示：

不满足稳定条件:dt<=Tn/pi,请重新输入符合稳定条件的时间步长dt 。