编制中心差分法程序报告(结构工程研究生作业)

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

北京工业大学结构工程

组员:胡建华 S201204111

马 恒 S201204112 陈相家S201204083 张力嘉S201204022 0前言

时域逐步积分法是数值分析方法，它只假设结构本构关系在一个微小的时间步距内是线性的。时域逐步积分法是结构动力问题中一个得到广泛研究的课题，并在结构动力反应计算中得到广泛应用。由于引进的假设条件不同，可以有各种不同的方法，比如中心差分法，线性加速度法，平均常加速度法，Wilson-θ法等，其中中心差分法精度最高。在本文中，通过编制中心差分法计算单自由度体系的动力反应的程序来了解其应用及稳定性。

1中心差分法原理

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

中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ∆=∆，则速度与加速度的中心差分近似为： t

u u u i i ∆-=

-+•

21

1 (a)

2

1

12t

u u u u i i i ∆+-=

-+•

• (b) 而离散时间点的运动为

)(),(),(i i i i i i t u u t u u t u u •

••

••

•

=== ( =i 0,1,2,3,……) 由体系运动方程为：0)()()(=++•

•

•t ku t u c t u m i (c)

将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t 时刻的运动方程：

0221

12

11=+∆-+∆+--+-+i i i i i i ku t u u c t u u u m

（d ）

在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：

12212)2()2()2(

-+∆-∆-∆--=∆+∆i i i u t c t

m u t m k u t c t m (e) 由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t 时刻的运动，如果需要可以用式（a ）和式（b ）求得体系的速度和加速度。

假设给定的初始条件为

),

0(),0(00•

•

==u u u u （g ）

由式（g ）确定1-u 。在零时刻速度和加速度的中心差分公式为：

t

u u u ∆-=

-•

21

10 （h ） `2

1

0102t

u u u u ∆+-=

-•

• （i ） 将式（i ）消去1u 得：020012

•

••

-∆+∆-=u t u t u u （j ） 而零时刻的加速度值0••u 可以用t ＝0时的运动方程 0000=++•

••ku u c u m 确定

即 )(1

000ku u c m

u --=•

•

• （k ）

这样就可以根据初始条件00,•

u u 和初始荷载0P ，就可以根据上式确定1-u 的值。

2中心差分法程序编制原理：

（1） 基本数据准备和初始条件计算

)(1

000ku u c m

u --=••

•

020012

•

••

-∆+∆-=u t u t u u （2） 计算等效刚度和中心差分计算公式中的相关系数

t c t m k ∆+∆=

22

22t m

k a ∆-=

t c

t

m b ∆-∆=22

（3） 根据i t 及i t 以前时刻的运动，计算1+i t 时刻的运动 1---=i i bu au P k P u i =+1 如果需要，可计算

t u u u i i ∆-=

-+•

21

1

2

1

12t

u u u u i i i ∆+-=

-+•

• （4）下一步计算用i+1代替i ，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。

以上为中心差分法逐步计算公式，其具有2阶精度，即误差)(02

t ∆∝ε；并且为有条件稳定，稳定条件为： π

n

T t ≤

∆

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

3算例

对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M=2000kg 、刚度K ＝50KN/m 、阻尼系数C ＝3KNs/m ，假设结构处于线弹性状态，用中心差分法计算结构的自由振动反应。

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

n

T t ≤

∆。

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

n

T t ωπ

2

=

≤

∆，而52000

50000

===

M K n ωrad/s,

则4.05

2

2

==

=

≤

∆n

n

T t ωπ

所以本次计算取t ∆＝0.1, 0.3, 0.4, 0.41, 0.42, 0.45分别进行计算

4 MATLAB 程序清单

[u,v,ac]=[M,C,K,u0,v0,time,dt]

% 本程序采用中心差分法计算结构的动力响应

% 本程序是既可以计算单自由度体系又可以计算多自由度体系，且均假设结构体系处于线