结构动力学振型分析

MATALAB 作业

某三层钢筋混凝土结构，结构的各层特性参数为：第一层到第三层质量m 分别为2400kg ，1200kg ，1200kg ，第一层到第三层刚度k 分别为3.3*10^4N/m,1.1*10^4N/m,0.66^4N/m.。地震采用acc_ElCentro_0.34g ，采样周期为0.02。

M3=1200kg

K3=0.66*10^4N/m.

M2=1200kg

K2=1.1*10^4N/m

M1=2400kg

K1=3.3*10^4N/m

用振型分解法求解结构地震反应的MATLAB 层序如下，编制该程序的程序框图以下所示

%振型分解法求解结构地震反应；主程序 clear

开始 输入地震参数和结构参数 计算结构振型与自振型频率 计算振型参与系数 计算单自由度体系的地震反应

求解结构的地震反应

输出结果

结束

clc

%地震波数据

xs=2*0.287;

dzhbo=load('acc_ElCentro_0.34g_0.02s.txt');

ag=dzhbo*0.01*xs;

dt=0.02;

ndzh=400;

cn=3; %cn为结构的层数，即质点数

m0=[2.4 1.2 1.2]*1e+3; %结构各层质量

k0=[3.3 1.1 0.66]*1e+5; %结构各层刚度

l=diag(ones(cn));

m=diag(m0); %计算质量矩阵

[ik]=matrixju(k0,cn); %计算刚度矩阵

[x,d]=eig(ik,m); %结构动力特性求解

d=diag(sqrt(d)); %求解结构圆频率

for i=1:cn;

[d1(i),j]=min(d);

xgd(:,i)=x(:,j);

d(j)=max(d)+1;

end %以此循环对所求频率和振型进行排序w=d1; %所求自振频率

x=xgd; %所求结构主振型

a1=2*w(1)*w(2)*(0.05*w(2)-0.07*w(1))/(w(2)^2-w(1)^2); a2=2*(0.07*w(2)-0.05*w(1))/(w(2)^2-w(1)^2);

for j=1:cn

x(:,j)=x(:,j)/x(cn,j);

znb0(j)=(a1+a2*w(j)^2)/2/w(j);

zhcan(j)=(x(:,j))'*m*l/((x(:,j))'*m*x(:,j));

%求解振型参数

[dlt(j,:),dltacceler(j,:)]=zxzj(znb0(j),w(j),ag);

end

%求解结构各层的地震反应

for i=1:cn;

disp1=0;

accel1=0;

for j=1:cn

disp0=zhcan(j)*dlt(j,:)*x(i,j);

accel0=zhcan(j)*dltacceler(j,:)*x(i,j);

disp1=disp1+disp0;

accel1=accel1+accel0;

end

disp(i,:)=disp1;

accel(i,:)=accel1;

end

t=0:dt:ndzh*dt;

%subplot(2,2,1)

%plot(t,disp(3,:)*1e+3,'k-')

%subplot(2,2,2)

plot(t,accel(3,:),'k-')

%子程序

%用于求解单自由度结构体系的地震反应

function [bx,acceler]=zxzj(znb,w,dag)

dt=0.02;

n=400;

x(1)=0;

dx(1)=0;

ddx(1)=0;

s=1+znb*dt*w+w^2*dt^2/6; %中间参数s

for i=1:n

a(i)=x(i)+dx(i)*dt+ddx(i)*dt^2/3;

b(i)=dx(i)+ddx(i)*dt/2;

ddx(i+1)=-1*(dag(i+1)+1*znb*w*b(i)+w^2*a(i))/s;%加速度dx(i+1)=b(i)+ddx(i+1)*dt/2;%速度

x(i+1)=a(i)+ddx(i+1)*dt^2/6;%位移

end

bx=x;

acceler=ddx;

%子程序

%刚度和阻尼矩阵的聚合

function [kcju]=matrixju(korc,cn)

kcju=zeros(cn);

for i=1:cn-1

kcju(i,i)=korc(i)+korc(i+1);

kcju(i,i+1)=-korc(i+1);

kcju(i+1,i)=-korc(i+1);

end

kcju(cn,cn)=korc(cn);

经程序求解，该结构的自振频率为

W= 4.8683 11.0666 15.4541

结构的振型矩阵为

X=

0.1634 -0.9238 2.7604

0.5691 -1.2267 -3.3423

1.0000 1.0000 1.0000

移位/ mm

时间/s

(a) 顶层位移反应

加速度/(m/s^2)

时间/s

(b)顶层位移反应