Bouc-wen模型模拟隔震支座

基础隔震结构简化成串联多质点系计算动力响应，采用Bouc-wen模型模拟隔震支座，用Newmark逐步积分法计算。

%%%%———结构动力学多自由度体系计算，包括隔震支座，采用Newmark法为基础，—————
%————考虑Bouc-Wen模型进行计算——————————————————
clear; %清除环境中存在的变量
clc; %清屏
tic; %计算程序运行时间
%—————————————初始量设定——————————————————————
%————质量矩阵——————————————
m=[400 520 520 520 520 520 520]*1e3; %各层的集中质量单位为kg
M=diag（m）; %质量矩阵
cn=length（m）;       %计算质点数
%————————刚度矩阵——————————
kkk=[218.7 800 800 800 800 800 800]*1e6; %各层层间刚度，单位为N/m
K=matrixju（kkkcn）; %初始刚度矩阵

%——————阻尼矩阵——————————————
c（1）=2*0.272*sqrt（m（1）*kkk（1））;
for i=2:cn
    c（i）=2*0.05*sqrt（m（i）*kkk（i））;%各层阻尼系数，单位是N*sec/m
end
 %为形成阻尼矩阵所用
C=matrixju（ccn）; %初始刚度矩阵

%————————————考虑隔震支座计算中涉及的参数取值————————————
mb=m（1）; %隔震支座的质量，单位为 kg
kb=kkk（1）; %隔震支座的初始刚度，单位为 N/m
cb=c（1）; %隔震支座的阻尼，单位为 N*sec/m
alphaB=0.1;
DyB=0.04; %单位换算成m
AB=1.0;
betaB=0.8;
gamaB=0.2;
nB=1;
%———————————Newmark法计算多自由度结构反应———————————————
%——————————————基本参数设定————————————————————
dt=0.02; %计算时间步长，同地震波数据时间步长一致
%for i=1:500
  %acceleration0（i）=4*sin（15.7*dt*i）; 
 %end
global EL;%导入加速度数据
acceleration0=EL*0.01*0.7/2.2; %换算成m/s^2 （1:100）
t=0:dt:（length（acceleration0）-1）*dt; %计算时长
gama=0.5; %定义gama值
beta=0.25; %定义beta值
L=ones（cn1）; %影响系数取值
p=-M*L*acceleration0‘; %地震力作用
%———————Newmark法———————————————————————
%——————参数计算——————
Gaa=M/（beta*dt）+gama*C/beta;
Gbb=M/（2*beta）+dt*C*（gama/（2*beta）-1）;
%————————初始量赋值——————————
displacement（:1）=zeros（cn1）; %初始相对位移向量
velocity（:1）=zeros（cn1）; %初始相对速度向量
acceleration（:1）=inv（M）*（p（:1）-K*displacement（:1）-C*velocity（:1））; %初始相对加速度向量

%——————隔震支座的初始量赋值——————————
RestoringForce（1）=0; % The stiffness restoring force of the lead-core rubber-bearing
xb（1）=0; %位移
dxb（1）=0; %速度，位移对时间的导数
vb（1）=0; %跟时间有关的参数
dvb（1）=0; %vb对时间的导数
%——————做循环计算—————————————————————————————
N=length（acceleration0）-1;
for i=1:N  
%——————————————————————————————————————
GKK=K+gama*C/（beta*dt）+M/（beta*dt^2）;
dp（:i）=（p（:i+1）-p（:i））+Gaa*velocity（:i）+Gbb*acceleration（:i）;
dq（:i）=GKK\dp（:i）;
dq1（:i）=gama*dq（:i）/（beta*dt）-gama*velocity（:i）/beta+dt*acceleration（:i）*（1-gama/（2*beta））;
dq2（:i）=dq（:i）/（beta*dt^2）-velocity（:i）/（beta*dt）-acceleration（:i）/（2*beta）;
displacement（:i+1）=displacement（:i）+dq（:i）;
velocity（:i+1）=velocity（:i）+dq1（:i）;
acceleration（:i+1）=acceleration（:i）+dq2（:i）;
%%————————隔震支座刚度矩阵的修正计算—————————————————
xb（i+1）=displacement（1i+1）; %隔震支座的相对位移
dxb（i+1）=velocity（1i+1）; %隔震支座的相对速度
vb（i+1）=vb（i）+dvb（i）*dt;
dvb（i+1）=（AB*dxb（i）-betaB*abs（dxb（i））*abs（vb（i））^（nB-1）*vb（i）-gamaB*dxb（i）*abs（vb（i）^nB））/DyB;
dvbdvx=（AB-betaB*abs（dxb（i+1））/dxb（i+1）*（abs（vb（i+1）））^（nB-1）*vb（i+1）-gamaB*（abs（vb（i+1））^nB））/DyB;
kbb=alphaB*kb+（1-alphaB）*kb*DyB*dvbdvx;
K（11）=kkk（2）+kbb; %修正后的刚度并进入下一次的循环计算中
%%——————————————————————————————————————
RestoringForce（i+1）=alphaB*kb*xb（i+1）+（1-alphaB）*kb*DyB*vb（i+1）;
end
%——————————————计算各层的绝对加速度————————————————
for i=1:length（m）
absacceleration（i:）=acceleration（i:）+acceleration0‘; %各层绝对加速度
q（i）=max（abs（absacceleration（i:）））/0.7;
end
qq=q‘   %加速度放大比

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       5238  2014-03-04 00:39  Is.m

----------- ---------  ---------- -----  ----

                 5238                    1