中心差分法在单由度中的应用.doc

中心差分法求解单自由度体系的自由振动问题 前言 时域逐步积分法是根据运动方程，引进某些假设，建立由t时刻状态向量、、到t＋时刻的状态向量、、的递推关系，从而从t＝0时刻的初始状态向量、、出发，逐步求出各时刻的状态向量，由于引进的假设条件不同，可以有各种不同的方法，下面主要介绍一种时域逐步积分方法－中心差分法。 中心差分法(central difference method)原理[1] 中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。 中心差分法只在相隔一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，，则速度与加速度的中心差分近似为： (a) (b) 而离散时间点的运动为 ( 0,1,2,3,……) 由体系运动方程为： (c) 将速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到时刻的运动方程： （d） 在（d）式中，假设和是已知的，即在及以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到： (e) 由式（e）就可以根据及以前时刻的运动，求得时刻的运动，如果需要可以用式（a）和式（b）求得体系的速度和加速度。 假设给定的初始条件为 （g） 由式（g）确定。在零时刻速度和加速度的中心差分公式为： （h） ` （i） 将式（i）消去得： （j） 而零时刻的加速度值可以用t＝0时的运动方程 确定 即 （k） 这样就可以根据初始条件和初始荷载，就可以根据上式确定的值。 下面给出采用中心差分法分析时的具体计算步骤： 基本数据准备和初始条件计算 计算等效刚度和中心差分计算公式中的相关系数 根据及以前时刻的运动，计算时刻的运动 如果需要，可计算 （4）下一步计算用i+1代替i，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。 以上为中心差分法逐步计算公式，其具有2阶精度，即误差；并且为有条件稳定，稳定条件为： 上式中，为结构的自振周期，对于多自由度结构体系则为结构的最小自振周期。 算例 对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M=2000kg、刚度K＝50KN/m、阻尼系数C＝3KNs/m，假设结构处于线弹性状态，用中心差分法计算结构的自由振动反应。 采用MATLAB语言编程，并以单自由度体系为例进行计算,设初位移u0＝0和初速度v0=0，取不同的步长分别计算，以验证中心差分法的稳定条件。 先计算，由稳定条件，而rad/s, 则 所以本次计算取＝0.1, 0.3, 0.4, 0.41, 0.42, 0.45分别进行计算 MATLAB程序清单 function [u,v,ac]=centraldifferent(M,C,K,u0,v0,time,dt) % 本程序采用中心差分法计算结构的动力响应 % 本程序是既可以计算单自由度体系又可以计算多自由度体系，且均假设结构体系处于线弹性状态； % ---------%%%%%输入参数%%%%%%%------------ % M------------质量矩阵 % C------------阻尼矩阵 % K------------刚度矩阵 % u0-----------初始位移 % v0-----------初始速度 % time---------模拟时间 % dt-----------时间步长 % ---%%%%%%输出值%%%%%%%%------ % u--------------位移 % v--------------速度 % ac-------------加速度 % -------%%%%%%%%中心差分法主要公式及原理%%%%%%%%%%----------- % MX+CX+KX=0 % M*(X(t+dt)-2*X(t)+X(t-