%%%%———結構動力學多自由度體系計算，包括隔震支座，采用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