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

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

1. 算例

首先，本文给出一算例， 结构在外力谐振荷载 P (t ) = P 0 sin θt 作用下，分别利用理论解法，杜哈梅积分， Wilson-θ 法求出该结构的位移时程反应。其中：

m = 3.5×103 kg ， P 0 = 1.0×104 N ， k =1.3584515×107 ，ξ=0.05 ，θ=52.3s ?1 ，ω=62.3s ?1 ，

?

D ω= ω 1-ξ2

=62.222 ，初始位移、速度v (0) = 0 ，v (0) = 0 ；

2. 算法验证

2.1 理论解法

运动方程为： mv+cv+kv=0P sin t ?由线性代数解出其理论解为:

]cos 2)sin -[(]

4)-[(t]

sin )

-(2cos [2]

4)-[(t]

sin )

0()0(cos )0([2

2222222022222

2

2

2

22

t t m P t m P e

v v t v e v D D

D t

D D

D t θξωθθθωθωξθωωωθωωξωξωθωξθωθωωξωωξωξω-++-+

?++++=--

由于初始位移v(0) =0 ，v(0) =0 ；则：

]cos 2)sin -[(]

4)-[(t]

sin )

-(2cos [2]

4)-[(222

222220

22222222220

t t m P t m P e v D D

D t

θξωθθθωθωξθωωωθωωξωξωθωξθωθξω-++

-+

?+=-

v(t ) =-3.115t

e

? 1.05269898×410-?[6.230cos62.222t ?18.106sin 62.222t]

+2.012808757 6

10-??[1146sin 52.3t ?325.829cos52.3t]

可以用 MA TLAB 进行编程分析，画图位移时程图，详细程序见附录。

2.2Wilson-θ法

Wilson-θ法是Wilson 于1966年基于线性加速度法的基础上提出一种无条件收敛的计算方法。该方法假定在θt ?时程步长内，体系的加速度反应按线性变化。对于地震持续时间内

的每一个微小时 段t ? ，从第一时段开始到最后一个时段，逐一的重复以下计算步骤，即得到结构地震反应的全过程。下面以第i+1时段（时刻时刻到1+i i t t ）为例：

{}{}[]{}

[][][][]

{}

[]{}[]{}{}[]{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}

{}{}{}{}{}{}{}[][]{}[][]{}附录。

移时程图，详细程序见进行编程分析，画图位可以用；时刻的相对加速度反应，求得第代入结构运动微分方程和最后将和速度增量移反应时刻的各质点的相对位得到第，和速度增量，分别加上位移增量时刻的地震反应为基础再以各质点在第；和速度增量内位移增量正常时段再求解第；内的加速度增量正常时段计算出第先利用位移增量的加长时段内各质点的计算MATLAB x x

C x

x x x x x x x

x x x x

x x x t x

t x t x t x x

x

t x x x i x

x

x x

x x x

x C x x x

P C K P K x x i i g i i i i i i i i i i i i i i i i i i i i i g i

i

i )(-t )5(t )4(6

31)(2t 1)3()2

(6

t 1i )2()2

3()36

(136;

t )1(1i 1

1i 1

1i 1i 11i 1i 1i 11i 1i 112

21

2

2

2

1

1

,+-+-+++++++++++++*

*

*

-*K M +M +=?+=?+=???+?+?=?+?=????+--?=???+?+++M +?M -=?+M +K =?=???=+ ττ

τ

ττθτ

τ

τ

ττθττττ

ττ

τ

2.3 杜哈梅积分

杜哈梅积分在考虑阻尼的情况是：

ττωτωωωξωωτξωξωd t p e m v v t v e v D t

t D

D D

D t )(sin )(1t]sin )

0()0(cos )0([0

)(-+

++

=?

---

可以用 MA TLAB 进行编程分析，画图位移时程图，详细程序见附录。

3. 位移时程反应对比分析

利用MATLAB将理论解法，杜哈梅积分，Wilson-θ法求解出来的位移时程反应画在同一张图中，进行比较分析。

从图中可以看出，以上三种方法得出来的位移时程曲线基本吻合，误差基本保持在5%以内，所以以上几种方法在求解相关问题上都具有一定的作用效力。

4.结论

本文通过一个简单的单自由度系统动力分析算例（仅作位移分析，其它分析雷同），基于MATLAB，将理论解法，杜哈梅积分法，逐步积分法（本文采用Wilson-θ法）进行相互验证，从最后的位移分析图对比上，可以很好的看出三种方法均能很好的彼此验证，从而说明了三种方法在相关问题上的作用效力。

附录：MA TLAB 源程序

%理论解，杜哈梅积分，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*1 0^(-6)*(1146*sin(52.3*t)-325.829*cos(52.3*t)); plot(t,v,'k')

hold on

%杜哈梅积分法

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))/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 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=kxin\dpxin;

dtjsdu=6*dxtao/(ct*(ct^2*dt^2))-6*sdu1/(ct*ct*dt)

-(3/ct)*jsdu1;

jsdu=jsdu1+dtjsdu;

dtsdu=(dt/2)*(jsdu+jsdu1);

sdu=sdu1+dtsdu;

dtwyi=dt*sdu1+(1/3)*dt^2*jsdu1+(dt^2/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('\fontsize{9}\fontname{ 黑体} 理论解','\fontsize{9}\fontname{黑体} 杜哈梅积分法','\fontsize{9}\fontname{黑体} wilson-{\theta}法')

%基于双线性滞回模型的单自由度体系的

地震能量分析程序

%质量57041kg，阻尼36612 N·s/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);

%加速度速度位移

Ei=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表示反向弹塑性，-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))*dth^2/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;

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: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))*dth^2/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;

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)<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)*dth^2/2+(a(ii)-a(ii-1))*dth^3/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/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;

pxmax=xp;

end

end

if pd==2 %反向弹性

k=k1;

if x(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))*dth^2/h/2);

ap=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)<(pxmax-2*Xy)

pd=-1;

b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)

x(ii-1)-pxmax+2*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))*dth^2/h/2);

ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;

pp=m*ap+c*vp+(pxmax*k2-k1*Xy-k2*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-2*Xy+dtx;

v(ii)=vp+dtv;

dtf=k2*dtx;

hfl(ii)=pxmax*k2-k1*Xy-k2*Xy+dtf;

a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;

end

end

if pd==-1 %反向弹塑性

k=k2;

if v(ii)>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)*dth^2/2+(a(ii)-a(ii-1))*dth^3/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/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)

pd=-1;

b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)

x(ii-1)-nxmax];% 拐点处理

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))*dth^2/h/2);

ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;

pp=m*ap+c*vp-k1*Xy+k2*(nxmax+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)=nxmax+dtx;

v(ii)=vp+dtv;

dtf=k2*dtx;

hfl(ii)=-k1*Xy+k2*(nxmax+Xy)+dtf;

a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;

elseif 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];% 拐点处理

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))*dth^2/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)/2);

dtx=dtpd/kd;

dtv=3*dtx/h-3*v(ii)-h*a(ii)/2;

x(ii+1)=x(ii)+dtx;

v(ii+1)=v(ii)+dtv;

dtf=k*dtx;

hfl(ii+1)=hfl(ii)+dtf;

a(ii+1)=(P(ii+1)-hfl(ii+1)-c*v(ii+1))/m;

Ei(ii+1)=-m*ug(ii+1)*v(ii+1)*h;

Ed(ii+1)=c*v(ii+1)^2*h;

Ek(ii+1)=m*a(ii+1)*v(ii+1)*h;

Eh(ii+1)=hfl(ii+1)*v(ii+1)*h;

for jj=1:num-1

EI(jj+1)=EI(jj)+Ei(jj+1);

ED(jj+1)=ED(jj)+Ed(jj+1);

EK(jj+1)=EK(jj)+Ek(jj+1);

EH(jj+1)=EH(jj)+Eh(jj+1);

end

end

for j=1:num

time(j)=h*(j-1);

end

fid=fopen('EL波弹塑性各种能.txt','wt'); %输出结果到EL波弹塑性各种能.txt fprintf(fid,'%g %g %g %g%g\n',[time;EI;ED;EH;EK]);

fclose(fid);

plot(time,EI);

hold on

plot(time,ED);

hold on

plot(time,EK);

hold on

plot(time,EH);