结构动力学作业姓名：学号：目录1.力插值法 (1)1.1分段常数插值法 (1)1.2分段线性插值法 (4)2.加速度插值法 (7)2.1常加速度法 (7)2.2线加速度法 (9)附录 (12)分段常数插值法源程序 (12)分段线性插值法源程序 (12)常加速度法源程序 (13)线加速度法源程序 (13)1.力插值法力插值法对结构的外荷载进行插值，分为分段常数插值法和分段线性插值法，这两种方法均适用于线性结构的动力反应计算。

1.1分段常数插值法图1-1为一个单自由度无阻尼系统，结构的刚度为k ，质量为m ，位移为y (t )，施加的外力为P (t )。

图1-2为矩形脉冲荷载的示意图，图中t d 表示作用的时间，P 0表示脉冲荷载的大小。

图1-1 单自由度无阻尼系统示意图图1-2 矩形脉冲荷载示意图对于一个满足静止初始条件的无阻尼单自由度体系来说，当施加一个t d 时间的矩形脉冲荷载，此时结构在t d 时间内的位移反应可以用杜哈梅积分得到：0()sin ()2 (1cos )(1cos ) (0)tst st d P y t t d m ty t y t t Tωττωπω=-=-=-≤≤⎰(1-1)如果结构本身有初始的位移和速度，那么叠加上结构自由振动的部分，结构的位移反应为：02()cos sin (1cos) (0)st d y ty t y t t y t t Tπωωω=++-≤≤ (1-2)图1-3 分段常数插值法微段示意图对于施加于结构任意大小的力，将其划分为Δt 的微段，每一段的荷载都为一个常数(每段相当于一个矩形的脉冲荷载)，如图1-3所示，则将每一段的位移和速度写成增量的形式为：1cos t sin t (1cos t)iii i y P y y kωωωω+=∆+∆+-∆ (1-3)i+1/sin t cos t sin t iii y P y y kωωωωω=-∆+∆+∆ (1-4)程序流程图如下i+1cos t sin t (1cos t)iii y P y y kωωωω=∆+∆+-∆i+1/sin t cos t sin ti i i y Py y kωωωωω=-∆+∆+∆图1-4 分段常数插值法流程图根据流程图可以编写相应的算法，利用MATLAB 进行编程，程序源代码见附录。

为了验证程序的正确性，本文选取的以下的例题进行验证。

对于一个单自由度的无阻尼结构，当其受到一个周期荷载时，其结构响应分为稳态解和瞬态解，由于没有阻尼的影响，其瞬态解并不会衰减，其理论表达式为：021()()(sin sin )1p x t t t k ωβωβ=-- (1-5)式中，()x t 为位移响应，0p 为激励，k 为刚度，β为荷载频率与固有振动频率之比，ω为荷载频率，ω为结构固有频率。

现令0p 为1，k 为1，则ω为1，ω取为2/3。

程序求得的解与解析解对比如图1-5所示（由于理论解与程序基本重合，所以将理论解乘以-1，方便比较）：位移y时间ta )位移速度v时间tb )速度图1-5 分段常数插值法结果验证由图1-5可知理论解与程序算得的解基本重合，可以验证程序的准确性。

1.2分段线性插值法与分段常数插值法不同，分段线性插值法将每一微段的力当成一个线性的直线，对于每一个微段，可看成一个矩形和一个三角形脉冲的叠加。

图1-6为分段线性插值微段示意图。

图1-6 分段线性插值法微段示意图对于无阻尼的体系，后一个时间步的位移和速度可由前一个时间步相应的值求得：11cos sin (1cos )(1sin )ii i i i y P P y y t t t t k k tωωωωωω+∆=∆+∆+-∆+-∆∆(1-6) 11/sin cos sin (1cos )i i i i i y P P y y t t t t k k tωωωωωωω+∆=-∆+∆+∆+-∆∆(1-7) 分段线性插值法的流程图如图1-7所示，与分段常数插值法仅仅是迭代的方式有所不一样。

11cos sin (1cos )(1sin )ii i i i y P P y y t t t t k k tωωωωωω+∆=∆+∆+-∆+-∆∆11/sin cos sin (1cos )i i i i i y P P y y t t t t k k tωωωωωωω+∆=-∆+∆+∆+-∆∆图1-7 分段线性插值法流程图程序源代码见附录，同样利用1.1节的算例进行验证，所得的结果如图1-8所示。

位移y时间ta )位移速度v时间tb )速度图1-8 分段线性插值法结果验证由上图可知程序的正确性。

2.加速度插值法加速度插值法也叫逐步积分法，其对加速度进行插值，可分为常加速度法和线加速度法。

2.1常加速度法图2-1常加速度法微段示意图对于一个单自由度结构，其运动方程为：my()y()()()t c t ky t p t ++=(2-1)将式（1-1）转变为增量方程：m y c y k y p ∆+∆+∆=∆(2-2)在通过逐步积分，将时间转化为一系列微小的时间段t ∆ ，如图2-1所示，现令111()()2i i i i y t y y t t t ++=+<≤，则t 时间的速度可表示为：111()(()())21()()(()())2i i t ti i t t i i i y t dt y t y t dty t y t y t y t t++=+-=+∆⎰⎰令t =t i +1,则i +1时刻速度可以表示为：111()2ii i i y y y y t ++=++∆(2-3)同理，位移可以表示为：2111()4i i i i i y y y t y y t ++=+∆++∆(2-4)将式(2-3)、(2-4)代入式(2-1)，即：1111my i i i i cy ky p ++++++=(2-5)此时，式(2-5)中只有1i y +为未知变量，可直接求出1i y +，之后再利用式(2-3)、(2-4)，可求出t i +1时刻的速度与位移。

算法的流程图如下所示：将荷载作用的时间划分为0243y t =∆t t∆∆*4(2)2i i i i mP P c y my t∆=∆+++∆22i i iy y y t∆=∆-∆24()2i i i i y y y t y t ∆=∆-∆-∆图2-2 常加速度法流程图算例验证的结果如下图所示，说明了该程序的正确性。

由于需要对加速度进行插值，此处增加了加速度验证。

位移y时间ta )位移速度v时间tb )速度加速度a时间tc )加速度图2-3 常加速法结果验证2.2线加速度法线加速度法与常加速度法原理类似，其速度与位移的增量方程与常加速度法对应的方程略有不同，图2-4为线加速度法微段示意图。

a )位移b )速度图2-4线加速度法微段示意图关于具体原理不在赘述，下图为线加速度法的流程图。

将荷载作用的时间划分为n 个分段，并计算263()()()()k t k t m c t t t=++∆∆6()()[()3()]()[3()()]2tP t P t m y t y t c t y t y t t ∆∆=∆++++∆/y P k∆=∆266()()()3()()y t y t y t y t tt∆=∆--∆∆3()()3()()2ty t y t y t y t t ∆∆=∆--∆0002332()2y y y t t=-∆∆图2-5 线加速度法流程图算例验证的结果如图2-6所示。

位移y时间ta )位移 速度v时间tb )速度加速度a时间tc )加速度图2-6 线加速法结果验证附录分段常数插值法源程序function x=inter_force_constant(p,w,dt,k,v0,y0)%分段常系数插值法%p代表输入的力，w为结构基本频率，dt为时间间隔,k为结构刚度%v0为初始的速度，y0为初始的位移%输出的矩阵第一列代表时间，第二列代表位移，第三列代表速度[r,~]=size(p);x=NaN(r,3);x(:,1)=p(:,1);x(1,2)=y0;x(1,3)=v0;for i=1:r-1x(i+1,2)=x(i,2)*cos(w*dt)+x(i,3)/w*sin(w*dt)+p(i,2)/k*(1-cos(w*dt));x(i+1,3)=(-x(i,2)*sin(w*dt)+x(i,3)/w*cos(w*dt)+p(i,2)/k*sin(w*dt))*w;end分段线性插值法源程序function x=inter_force_line(p,w,dt,k,v0,y0)%分段线性插值法%p代表输入的力，w为结构基本频率，dt为时间间隔,k为结构刚度%v0为初始的速度，y0为初始的位移%输出的矩阵第一列代表时间，第二列代表位移，第三列代表速度[r,~]=size(p);x=NaN(r,3);x(:,1)=p(:,1);x(1,2)=y0;x(1,3)=v0;for i=1:r-1x(i+1,2)=x(i,2)*cos(w*dt)+x(i,3)/w*sin(w*dt)+p(i,2)/k*(1-cos(w*dt))+(p(i+1,2)-p(i,2))/k/w/dt*(w*dt-sin(w*dt));x(i+1,3)=(-x(i,2)*sin(w*dt)+x(i,3)/w*cos(w*dt)+p(i,2)/k*sin(w*dt)+(p(i+1,2)-p(i,2))/k/w/dt*(1-cos(w*dt)))*w;end常加速度法源程序function x=inter_a_constant(p,w,m,keci,dt,v0,y0,k)%p代表输入的荷载，w为结构基本频率，keci为阻尼比，dt为时间间隔,m为结构质量,k 为结构刚度%v0为初始的速度，y0为初始的位移%输出的矩阵第一列代表时间，第二列代表位移,第三列代表速度，第四列代表加速度[r,~]=size(p);x=NaN(r,4);x(:,1)=p(:,1);x(1,2)=y0;x(1,3)=v0;x(1,4)=4/3/dt/dt*y0;c=2*keci*w;for i=1:r-1K=k+2*c/dt+4*m/dt/dt;dP=p(i+1,2)-p(i,2)+(4*m/dt+2*c)*x(i,3)+2*m*x(i,4);dY=dP/K;dV=2/dt*dY-2*x(i,3);dA=4/dt/dt*(dY-x(i,3)*dt)-2*x(i,4);x(i+1,2)=x(i,2)+dY;x(i+1,3)=x(i,3)+dV;x(i+1,4)=x(i,4)+dA;end线加速度法源程序function x=inter_a_line(p,w,m,keci,dt,v0,y0,k)%p代表输入的荷载，w为结构基本频率，keci为阻尼比，dt为时间间隔,m为结构质量,k 为结构刚度%v0为初始的速度，y0为初始的位移%输出的矩阵第一列代表时间，第二列代表位移,第三列代表速度，第四列代表加速度[r,~]=size(p);x=NaN(r,4);x(:,1)=p(:,1);x(1,2)=y0;x(1,3)=v0;x(1,4)=3/2/dt/dt*y0-3/2/dt*v0;c=2*keci*w;for i=1:r-1K=k+3*c/dt+6*m/dt/dt;dP=p(i+1,2)-p(i,2)+m*(6/dt*x(i,3)+3*x(i,4))+c*(3*x(i,3)+dt/2*x(i,4));dY=dP/K;dV=3/dt*dY-3*x(i,3)-dt/2*x(i,4);dA=6/dt/dt*dY-6/dt*x(i,3)-3*x(i,4);x(i+1,2)=x(i,2)+dY;x(i+1,3)=x(i,3)+dV;x(i+1,4)=x(i,4)+dA;end。