基于 MATLAB 实现对结构动力响应的几种算法的验证

1、 基于 MATLAB 实现对结构动力响应的几种算法的验证1. 算例 首先，本文给出一算例， 结构在外力谐振荷载 P(t) = P0 sint 作用下，分别利用理论解法，杜哈梅积分， Wilson- 法求出该结构的位移时程反应。其中：m = 3.5103 kg ， P0 = 1.0104 N ， k =1.3584515107 ，=0.05 ，=52.3s1 ，=62.3s1 ， = 1-2 =62.222 ，初始位移、速度v(0) = 0 ，v(0) = 0 ；2. 算法验证 2.1 理论解法 运动方程为： mv+cv+kv=sin由线性代数解出其理论解为:由于初始位移v(0) =0 ，v(0

2、) =0 ；则： v(t ) =1.052698986.230cos62.222t 18.106sin 62.222t +2.012808757 1146sin 52.3t 325.829cos52.3t 可以用 MATLAB 进行编程分析，画图位移时程图，详细程序见附录。 2.2Wilson-法 Wilson-法是Wilson于1966年基于线性加速度法的基础上提出一种无条件收敛的计算方法。该方法假定在时程步长内，体系的加速度反应按线性变化。对于地震持续时间内的每一个微小时 段 ，从第一时段开始到最后一个时段，逐一的重复以下计算步骤，即得到结构地震反应的全过程。下面以第i+1时段（）为例：2

3、.3 杜哈梅积分 杜哈梅积分在考虑阻尼的情况是： 可以用 MATLAB 进行编程分析，画图位移时程图，详细程序见附录。 3. 位移时程反应对比分析利用 MATLAB 将理论解法，杜哈梅积分， Wilson- 法求解出来的位移时程反应画在同一张图 中，进行比较分析。 从图中可以看出，以上三种方法得出来的位移时程曲线基本吻合，误差基本保持在 5%以内，所以以上几种方法在求解相关问题上都具有一定的作用效力。4. 结论 本文通过一个简单的单自由度系统动力分析算例（仅作位移分析，其它分析雷同），基于 MATLAB，将理论解法，杜哈梅积分法，逐步积分法（本文采用 Wilson- 法）进行相互验证，从最后的

4、位移分析图对比上，可以很好的看出三种方法均能很好的彼此验证，从而说明了三种方法在相关问题上的作用效力。附录：MATLAB 源程序%理论解，杜哈梅积分，Wilson-法程序clc; clear h1=figure(8); set( h1, color,w) %理论解法t=0:0.01:1; v=1.05269898*10(-4)*exp(-3.115*t).*(6.230*cos(62.222*t)-18.106*sin(62.222*t)+2.012808757*10(-6)*(1146*sin(52.3*t)-325.829*cos(52.3*t); plot(t,v,k) hold on

5、%杜哈梅积分法aa=1;%输入时间长度bb=0.01;%输入精度t=bb:bb:aa; t1=t; theta=52.3;%输入荷载频率w=62.3;%输入自振频率m=3500;%输入质量p0=10000;%输入荷载幅值p0=p0*ones(1,aa/bb); p=p0.*sin(theta*t);%荷载函数for i=1:(aa/bb) for j=1:i canshu1(j)=p(j)/(m*w)*bb*sin(w*(t(i)-t1(j);%杜哈梅积分中的被积函数end y(i)=sum(canshu1);%位移值end for i=1:aa/bb-1 v1(i)=(y(i+1)-y(i)

6、/bb;%计算速度end for i=1:(aa/bb-2) a(i)=(v1(i+1)-v1(i)/bb;%计算加速度end hold on plot(t,y,r)%画位移图hold on %Wilson-法dt=0.01; ct=1.4;ndzh=100; k=13584515; c=21805; t=0:dt:ndzh*dt; ag=10000*sin(52.3*t); ag1=ag(1:ndzh); ag2=ag(2:ndzh+1); agtao=ct*(ag2-ag1); wyi1=0; sdu1=0; jsdu1=0; wyimt=0; sdumt=0; jsdumt=0; for

7、 i=1:ndzh kxin=k+(3/(ct*dt)*c+(6/(ct*ct*dt*dt)*m; %kxin为新的刚度dpxin=-m*agtao(i)+m*(6/(ct*dt)*sdu1+3*jsdu1)+c*(3*sdu1+ct*dt/2*jsdu1); %新的力增量dxtao=kxindpxin; dtjsdu=6*dxtao/(ct*(ct2*dt2)-6*sdu1/(ct*ct*dt)-(3/ct)*jsdu1; jsdu=jsdu1+dtjsdu; dtsdu=(dt/2)*(jsdu+jsdu1); sdu=sdu1+dtsdu; dtwyi=dt*sdu1+(1/3)*dt2

8、*jsdu1+(dt2/6)*jsdu; wyi=wyi1+dtwyi; jsdu=-m(m*ag2(i)+c*sdu+k*wyi); % 调整加速度wyi1=wyi; sdu1=sdu; jsdu1=jsdu; wyimt=wyimt wyi; sdumt=sdumt sdu; jsdumt=jsdumt jsdu; end plot(t,-wyimt/3000,b); hold off legend(fontsize9fontname 黑体 理论解,fontsize9fontname黑体 杜哈梅积分法,fontsize9fontname黑体 wilson-theta法)%基于双线性滞回模型

9、的单自由度体系的地震能量分析程序 %质量57041kg，阻尼36612 Ns/m,初始刚度2350000N/m,刚度折减 系数0.2，屈服位移0.01m，采用ELCENTRO波%参数替换直接在下面修改，然后运行clc format long; m=57041;%质量ug=importdata(ELCENTRO.txt);%地震波txt 文件ug=ug/100; P=-m*ug; num=size(P,1); c=36612;%阻尼k1=2350000;%初始刚度k2=k1*0.2; a=zeros(1,num);v=zeros(1,num);x=zeros(1,num);%加速度 速度 位移E

10、i=zeros(1,num);Ek=zeros(1,num);Ed=zeros(1,num);Eh=zeros(1,num);%输入能 动能 阻尼耗能 滞回耗能EI=zeros(1,num);EK=zeros(1,num);ED=zeros(1,num);EH=zeros(1,num);%累积的各种能量time=zeros(1,num); a(1)=P(1)/m; hfl=zeros(1,num); h=0.02;%地震波采样间隔Xy=0.01;%屈服位移pxmax=0; nxmax=0; pd=0;%双线型滞回模型折线的标识，0表示弹性，1表示正向弹塑性，2表示反向弹性，-1表示反向弹塑性，

11、-2表示正向弹性for ii=1:num if pd=0 %弹性阶段 k=k1; if x(ii)Xy pd=1; b=(a(ii)-a(ii-1)/6/h a(ii-1)/2 v(ii-1) x(ii-1)-Xy;% 拐点处理 d=roots(b); for j=1:length(d) if isreal(d(j)=1 dth=d(j); end end hp=h-dth; vp=v(ii-1)+a(ii-1)*dth+(a(ii)-a(ii-1)*dth2/h/2); ap=a(ii-1)+(a(ii)-a(ii-1)*dth/h; pp=m*ap+c*vp+k1*Xy; kd=k2+3*

12、c/hp+6*m/hp/hp; dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2); dtx=dtpd/kd; dtv=3*dtx/hp-3*vp-hp*ap/2; x(ii)=Xy+dtx; v(ii)=vp+dtv; dtf=k2*dtx; hfl(ii)=k1*Xy+dtf; a(ii)=(P(ii)-hfl(ii+1)-c*v(ii)/m; elseif x(ii)-Xy pd=-1; b=(a(ii)-a(ii-1)/6/h a(ii-1)/2 v(ii-1) x(ii-1)+Xy;% 拐点处理 d=roots(b); for j=1:

13、length(d) if isreal(d(j)=1 dth=d(j); end end hp=h-dth; vp=v(ii-1)+a(ii-1)*dth+(a(ii)-a(ii-1)*dth2/h/2); ap=a(ii-1)+(a(ii)-a(ii-1)*dth/h; pp=m*ap+c*vp-k1*Xy; kd=k2+3*c/hp+6*m/hp/hp; dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2); dtx=dtpd/kd; dtv=3*dtx/hp-3*vp-hp*ap/2; x(ii)=-Xy+dtx; v(ii)=vp+dtv;

14、 dtf=k2*dtx; hfl(ii)=-k1*Xy+dtf; a(ii)=(P(ii)-hfl(ii)-c*v(ii)/m; end end i f pd=1 %正向弹塑性 k=k2; if v(ii)pxmax pd=1; b=(a(ii)-a(ii-1)/6/h a(ii-1)/2 v(ii-1) x(ii-1)-pxmax;% 拐点处理 d=roots(b); for j=1:length(d) if isreal(d(j)=1 dth=d(j); end end hp=h-dth;vp=v(ii-1)+a(ii-1)*dth+(a(ii)-a(ii-1)*dth2/h/2); ap

15、=a(ii-1)+(a(ii)-a(ii-1)*dth/h; pp=m*ap+c*vp+k1*Xy+k2*(pxmax-Xy); kd=k2+3*c/hp+6*m/hp/hp; dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2); dtx=dtpd/kd; dtv=3*dtx/hp-3*vp-hp*ap/2; x(ii)=pxmax+dtx; v(ii)=vp+dtv; dtf=k2*dtx; hfl(ii)=k1*Xy+k2*(pxmax-Xy)+dtf; a(ii)=(P(ii)-hfl(ii)-c*v(ii)/m; elseif x(ii)

16、0 pd=-2; b=(a(ii)-a(ii-1)/2/h a(ii-1) v(ii-1);% 拐点处理 d=roots(b); for j=1:length(d) if isreal(d(j)=1 dth=d(j); end end hp=h-dth; xp=x(ii-1)+v(ii-1)*dth+a(ii-1)*dth2/2+(a(ii)-a(ii-1)*dth3/h/6; ap=a(ii-1)+(a(ii)-a(ii-1)*dth/h; pp=m*ap+c*0-k1*Xy+k2*(xp+Xy); kd=k1+3*c/hp+6*m/hp/hp; dtpd=P(ii+1)-pp+m*(6*0

17、/hp+3*ap)+c*(3*0+hp*ap/2); dtx=dtpd/kd; dtv=3*dtx/hp-3*0-hp*ap/2; x(ii)=xp+dtx; v(ii)=0+dtv; dtf=k1*dtx; hfl(ii)=-k1*Xy+k2*(xp+Xy)+dtf; a(ii)=(P(ii)-hfl(ii)-c*v(ii)/m; nxmax=x(ii); end end if pd=-2 %正向弹性 k=k1; if x(ii)(nxmax+2*Xy) pd=1; b=(a(ii)-a(ii-1)/6/h a(ii-1)/2 v(ii-1) x(ii-1)-nxmax-2*Xy;% 拐点处

18、理 d=roots(b); for j=1:length(d) if isreal(d(j)=1 dth=d(j); end end hp=h-dth; vp=v(ii-1)+a(ii-1)*dth+(a(ii)-a(ii-1)*dth2/h/2); ap=a(ii-1)+(a(ii)-a(ii-1)*dth/h; pp=m*ap+c*vp+(k1*Xy+nxmax*k2+Xy*k2); kd=k2+3*c/hp+6*m/hp/hp; dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2); dtx=dtpd/kd; dtv=3*dtx/hp-3*vp-hp*ap/2; x(ii)=nxmax+2*Xy+dtx; v(ii)=vp+dtv; dtf=k2*dtx; hfl(ii)=k1*Xy+nxmax*k2+Xy*k2+dtf; a(ii)=(P(ii)-hfl(ii)-c*v(ii)/m; end end i f ii=num break; end kd=k+3*c/h+6*m/h/h; dtpd=P(ii+1)-P(ii)+m*(6*v(ii)/h+3*a(ii)+c*(3*v(ii)+h*a(ii)/