中心差分法计算程序编程.doc

中心差分法计算程序编程 姓名：张泽伟 学号： 电话： 一、中心差分法程序原理说明 1.1 中心差分法思路 中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。 1.2 中心差分法原理 中心差分法只在相隔一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，，则速度与加速度的中心差分近似为： (a) (b) 而离散时间点的运动为 ( 0,1,2,3,……) 由体系运动方程为： (c) 将速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到时刻的运动方程： （d） 在（d）式中，假设和是已知的，即在及以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到： (e) 由式（e）就可以根据及以前时刻的运动，求得时刻的运动，如果需要可以用式（a）和式（b）求得体系的速度和加速度。 1.3 初始条件转化 假设给定的初始条件为 （g） 由式（g）确定。在零时刻速度和加速度的中心差分公式为： （h） ` （i） 将式（i）消去得： （j） 而零时刻的加速度值可以用t＝0时的运动方程 确定 即 （k） 这样就可以根据初始条件和初始荷载，就可以根据上式确定的值。 1.4 中心差分法编程思路 ① 基本数据准备和初始条件计算： ② 计算等效刚度和中心差分计算公式中的相关系数： ③ 根据及以前时刻的运动，计算时刻的运动： ④ 下一步计算用i+1代替i，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。 1.5 中心差分法稳定条件 以上为中心差分法逐步计算公式，其具有2阶精度，即误差；并且为有条件稳定，稳定条件为： 程序框图 根据中心差分法的原理，可以得出本程序的主要程序思想，以下面框图的形式展示出来： 编辑回调函数 编辑回调函数 调用输入数据 调用输入数据 计算初始加速度和位移前一时刻值 计算初始加速度和位移前一时刻值 开始循环体（步数控制） 开始循环体（步数控制） 计算等效刚度和系数a、b 计算等效刚度和系数a、b 计算位移、速度、加速度并保存 计算位移、速度、加速度并保存 分析稳定条件 分析稳定条件 绘图 绘图 程序清单 %m,k,c分别为质量、刚度、阻尼 %p0,dt,t分别为外荷载幅值、时间步距、总时间 %u0,v0为初始条件初位移和初速度 %u,v,ac分别为位移、速度、加速度反应 ek=等效刚度；p=荷载；ep=等效荷载 %定义矩阵X0=input(请按格式和顺序输入初始矩阵，如X0=[m,k,c,u0,v0,t,P0,dt], m=X0(1,1);k=X0(1,2);c=X0(1,3);u0=X0(1,4); %分别取出其中的参数： v0=X0(1,5);t=X0(1,6);P0=X0(1,7);dt=X0(1,8) t=[0:dt:t]; %将时间分步，采用等时间步长； [mm,nn]=size(t); %计算t的向量长度，得出步数； u=zeros(size(t)); %设定存储u的矩阵； v=zeros(size(t)); %设定存储v的矩阵； ac=zeros(size(t)); %设定存储ac的矩阵； u(:,2)=u0; %赋值向量第2项为u0； v(:,2)=v0; %赋值向量第2项为v0； ac(:,2)=(P0-c*v(:,2)-k*u(:,2))/m; %求出初始加速度ac0； u(:,1)=u(:,2)-dt*v(:,2)+((dt)^2)*ac(:,2)/2;