%%%%———結(jié)構(gòu)動(dòng)力學(xué)多自由度體系計(jì)算，包括隔震支座，采用Newmark法為基礎(chǔ)，—————
%————考慮Bouc-Wen模型進(jìn)行計(jì)算——————————————————
clear;?%清除環(huán)境中存在的變量
clc;?%清屏
tic;?%計(jì)算程序運(yùn)行時(shí)間
%—————————————初始量設(shè)定——————————————————————
%————質(zhì)量矩陣——————————————
m=[400?520?520?520?520?520?520]*1e3;?%各層的集中質(zhì)量單位為kg
M=diag（m）;?%質(zhì)量矩陣
cn=length（m）;???????%計(jì)算質(zhì)點(diǎn)數(shù)
%————————?jiǎng)偠染仃嚒?br/>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））;%各層阻尼系數(shù)，單位是N*sec/m
end
?%為形成阻尼矩陣所用
C=matrixju（ccn）;?%初始剛度矩陣

%————————————考慮隔震支座計(jì)算中涉及的參數(shù)取值————————————
mb=m（1）;?%隔震支座的質(zhì)量，單位為?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法計(jì)算多自由度結(jié)構(gòu)反應(yīng)———————————————
%——————————————基本參數(shù)設(shè)定————————————————————
dt=0.02;?%計(jì)算時(shí)間步長(zhǎng)，同地震波數(shù)據(jù)時(shí)間步長(zhǎng)一致
%for?i=1:500
??%acceleration0（i）=4*sin（15.7*dt*i）;?
?%end
global?EL;%導(dǎo)入加速度數(shù)據(jù)
acceleration0=EL*0.01*0.7/2.2;?%換算成m/s^2?（1:100）
t=0:dt:（length（acceleration0）-1）*dt;?%計(jì)算時(shí)長(zhǎng)
gama=0.5;?%定義gama值
beta=0.25;?%定義beta值
L=ones（cn1）;?%影響系數(shù)取值
p=-M*L*acceleration0‘;?%地震力作用
%———————Newmark法———————————————————————
%——————參數(shù)計(jì)算——————
Gaa=M/（beta*dt）+gama*C/beta;
Gbb=M/（2*beta）+dt*C*（gama/（2*beta）-1）;
%————————初始量賦值——————————
displacement（:1）=zeros（cn1）;?%初始相對(duì)位移向量
velocity（:1）=zeros（cn1）;?%初始相對(duì)速度向量
acceleration（:1）=inv（M）*（p（:1）-K*displacement（:1）-C*velocity（:1））;?%初始相對(duì)加速度向量

%——————隔震支座的初始量賦值——————————
RestoringForce（1）=0;?%?The?stiffness?restoring?force?of?the?lead-core?rubber-bearing
xb（1）=0;?%位移
dxb（1）=0;?%速度，位移對(duì)時(shí)間的導(dǎo)數(shù)
vb（1）=0;?%跟時(shí)間有關(guān)的參數(shù)
dvb（1）=0;?%vb對(duì)時(shí)間的導(dǎo)數(shù)
%——————做循環(huán)計(jì)算—————————————————————————————
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）;
%%————————隔震支座剛度矩陣的修正計(jì)算—————————————————
xb（i+1）=displacement（1i+1）;?%隔震支座的相對(duì)位移
dxb（i+1）=velocity（1i+1）;?%隔震支座的相對(duì)速度
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;?%修正后的剛度并進(jìn)入下一次的循環(huán)計(jì)算中
%%——————————————————————————————————————
RestoringForce（i+1）=alphaB*kb*xb（i+1）+（1-alphaB）*kb*DyB*vb（i+1）;
end
%——————————————計(jì)算各層的絕對(duì)加速度————————————————
for?i=1:length（m）
absacceleration（i:）=acceleration（i:）+acceleration0‘;?%各層絕對(duì)加速度
q（i）=max（abs（absacceleration（i:）））/0.7;
end
qq=q‘???%加速度放大比

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----

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

-----------?---------??----------?-----??----

?????????????????5238????????????????????1