地震工程作业.doc

地震工程课程大作业

作业1
绘制1940 El Centro，N-S分量地震动的绝对加速度、相对速度和相对位移反应谱.地震动：在PEER Ground Motion Database自行下载经典的1940 El Centro，N—S分量。 要求：在此模板内完成，A4纸打印。自编程序与软件（Bispec或Seismosigna等）计算反应谱进行对比。提交自编写程序。
Matlab程序：
clear
fid = fopen（'E：\Earthquake\El centro。txt'）；
［Accelerate，count] = fscanf（fid,'%g’）; ％count 读入的记录的量
Accelerate=9.8＊Accelerate’; %单位统一为 m和s
time=0：0。02：(count—1）*0。02； ％单位 s
Displace=zeros(1，count）; ％相对位移
Velocity=zeros（1，count）； ％相对速度
AbsAcce=zeros（1，count）； %绝对加速度
DampA=［0.00，0.02,0。05］； %三个阻尼比
TA=0。0:0。02：4； ％TA=0.000001：0.02：4；结构周期
Dt=0。02； ％地震记录的步长
％记录计算得到的反应，MDis为某阻尼时最大相对位移
％MVel为某阻尼时最大相对速度，MAcc某阻尼时最大绝对加速度
MDis=zeros(3，length（TA)）;
MVel=zeros（3,length（TA)）;
MAcc=zeros（3，length(TA）)；
j=1；
for Damp=[0。00，0。05，0。1］
t=1；
for T=0。0:0.02:4
Frcy=2*pi/T ；
DamFrcy=Frcy*sqrt（1—Damp＊Damp）;
e_t=exp（—Damp＊Frcy＊Dt）；
s=sin（DamFrcy＊Dt)；
c=cos(DamFrcy*Dt）；
A=zeros（2，2）;
A（1，1)=e_t*（s＊Damp/sqrt（1—Damp*Damp)+c)；
A（1，2）=e_t＊s/DamFrcy;
A(2，1）=-Frcy＊e_t＊s/sqrt(1—Damp＊Damp)；
A(2，2)=e_t＊（—s＊Damp/sqrt（1-Damp＊Damp)+c）；
d_f=(2*Damp^2—1）/（Frcy^2*Dt）；
d_3t=Damp/（Frcy^3*Dt）;
B=zeros(2，2)；
B（1,1)=e_t*（(d_f+Damp/Frcy）＊s/DamFrcy+(2*d_3t+1/Frcy^2)＊c）-2＊d_3t；
B(1，2)=-e_t＊（d_f＊s/DamFrcy+2*d_3t＊c)-1/Frcy^2+2*d_3t； B（2，1）=e_t*（（d_f+Damp/Frcy)*(c—Damp/sqrt（1—Damp^2)＊s）—(2＊d_3t+1/Frcy^2)*
(DamFrcy＊s+Damp*Frcy＊c））+1/(Frcy^2＊Dt）；
B（2,2)=e_t＊（1/（Frcy^2＊Dt）＊c+s*Damp/（Frcy*DamFrcy*Dt）)-1/（Frcy^2*Dt）;
for i=1：（count-1）
Displace(i+1)=A（1,1）*Displace（i）+A(1，2）＊Velocity（i)+B(1，1)*Accelerate（i）+B（1，2）*
Accelerate（i+1）;
Velocity（i+1）=A(2，1）＊Displace(i）+A（2，2）*Velocity（i)+B(2，1）＊Accelerate(i）+B（2，2）*
Accelerate（i+1);
AbsAcce(i+1）=-2*Damp*Frcy＊Velocity（i+1）—Frcy^2＊Displace（i+1）；
end
MDis（j，t)=max（abs(Displace)）；
MVel(j,t）=max（abs(Velocity））；
if T==0。0
MAcc(j,t）=max（abs（Accelerate））;
else
MAcc(j,t)=max（abs(AbsAcce））；
end
Displace=zeros（1,count）;
Velocity=zeros(1，count)；
AbsAcce=zeros（1，count)；
t=t+1；
end
j=j+1；
end
close all
figure %绘制位移反应谱
plot（TA,MDis（1，:),’—b’,TA,MDis（2，：）,’-r’，TA，MDis（3，:），’:k’)
title（’Displacement’）
xlabel（'Tn(s）'）
ylabel(’Displacement（m)'）
legend（’ζ=0'，'ζ=0.02’，’ζ=0。05'）
grid
figure %绘制速度反应谱
plot(TA，MVel（1，：），’-b',TA，MVel（2,:），’—r',TA，MVel(3,：），':k'）
title('Velocity’）
xlabel('Tn（s)’)
ylabel(’velocity（m/s）’）
legend（'ζ=0'，'ζ=0。02’，’ζ=0。05’)
grid
figure ％绘制绝对加速度反应谱
plot（TA,MAcc（1，：)，'—b’，TA,MAcc（2，：)，’-r’，TA,MAcc（3，:），':k’)
title（’Absolute Acceleration')
xlabel（’Tn（s)’）
ylabel（'absolute acceleration(m/s^2）’)
legend(’ζ=0’,'ζ=0.02','ζ=0。05’）
grid
由此得到下图反应谱与由Seismosigna得到的反应谱作出比较
相对位移反应谱(上）与Seismosigna计算相对位移反应谱（下）的比较
相对速度反应谱（上）与Seismosigna计算相对速度反应谱（下)的比较
绝对加速度反应谱(左图）与Seismosigna计算绝对加反应谱（右图)的比较
作业2、
使用中心差分方法和Newmark方法分别计算下面单自由度结构的地震反应.题：考虑下面的具有2个自由度的一个简单系统,分析时间步长Δt=0。28sec的情况，求每时刻的加速度、速度、位移
求：
要 求:在此模板内完成，A4纸打印。
MATLAB程序：
clear
％求自振圆频率与振型
m=［2 0；0 1］;
k=[6 —2；—2 4］;
［v，d］=eig(inv(m）＊k）;
w=sqrt（d）;
fai1=v(:，1）。/v（1,1）；
fai2=v(：,2)./v（1,2）；
fai=［fai1 fai2]；
mm=fai’＊m＊fai;kk=fai’＊k＊fai;％等效质量矩阵及等效刚度矩阵
p0=[0;10］;
％中心差分法
s0=［0;0］；u0=［0;0］；%初速度及初位移均为0
％以下过程均采用q、q1、q2的形式分别表示位移、速度、加速度
q0=［fai1’*m*s0。/(fai1’＊m＊fai1);fai2’＊m＊s0。/(fai2’＊m*fai2）］；
q01=［fai1’＊m＊u0。/（fai1'＊m＊fai1）;fai2'＊m*u0。/(fai2'*m*fai2)］；
P0=fai’＊p0；
q02=mm^(-1)＊（P0—kk＊q0);
t=0。28；
q(:，1）=q0—t＊q01+t^2*q02/2;％q(:,1)表示q—1,即初始时刻前一时刻的位移
q（：，2）=[0；0］；％q（:，2)表示初始时刻的位移
Kh=mm/（t^2)；
a=mm/（t^2）；
b=kk—2＊mm/（t^2）;
P=fai'＊p0；
for i=2：7;
Ph（：，i）=P—a*q(:，i—1）-b＊q（：，i）；
q（：，i+1）=Kh^（—1）＊Ph（：,i);
q1（:，i）=(q（:,i+1）—q（:，i—1）)/（2＊t)；
q2(：，i）=(q(:，i+1）-2＊q（：,i)+q（:，i-1））/(t^2）；
end
u=fai＊q；
u1=fai＊q1;
u2=fai*q2；
for i=1：6
U（：，i）=u（：，i+1）；
U1（：，i）=u1（：,i+1)；
U2(：，i）=u2(:,i+1)；
end
U
U1
U2
％Newmark法（采用线加速度法）
gamma=1/2；
beta=1/6;
Kh1=kk+mm/（beta*t^2)；
a1=mm/(beta＊t);
b1=mm/（2*beta);
detaP=［0 0 0 0 0；0 0 0 0 0];
q1（：，1）=[0；0];
q11(：，1)=［0;0］；
q12（：,1）=q02;
for i=1:5
detaPh（：，i)=detaP(：，i）+a1＊q11(:，i）+b1＊q12（:，i）；
detaq（：，i)=Kh1^（-1）*detaPh（:,i）；
detaq1(:，i)=gamma＊detaq（：，i)/（beta*t）—gamma＊q11（：，i)/beta+t＊(1—gamma/（2＊beta））＊q12(：,i);
detaq2(：，i）=detaq(：,i)/(beta＊t^2)—q11（:，i)/（beta＊t）-q12（：，i）/(2＊beta）；
q1（：，i+1）=q1（：，i）+detaq（：,i）；
q11（：,i+1）=q11（：，i）+detaq1(：,i）;
q12（:，i+1）=q12(:，i)+detaq2(：，i）;
end
v=fai＊q1；
v1=fai*q11；
v2=fai＊q12;
for i=1：6
V（:，i)=v（:，i)；
V1（：，i）=v1(：，i）；
V2（：，i）=v2（:,i）；
end
V
V1
V2
表1中心差分法
Time
0。28s
0。56s
0。84s
1.12s
1。40s
位移
（m）
0
0。0307
0。1675
0。4871
1。0170
0.3920
1.4451
2。8338
4。1441
5.0152
速度
(m/s）
0。0549
0。2992
0。8149
1.5169
2.1675
2。5805
4。3603
4。8197
3。8954
1。9873
加速度
（m/s^2）
0.3920
1.3529
2.3312
2.6828
1。9642
8.4320
402812
-1。0001
—5.6022
-8.0268
表2 Newmark法
Time
0.28s
0.56s
0。84s
1.12s
1.40s
位移
(m)
0。0047
0.0444
0.1826
0。4850
0。9780
0.3726
1.3809
2。7317
4。0447
4。9744
速度
（m/s)
0。0502
0.2751
0。7555
1。4238
2.0720
2。5926
4。4244
4.9850
4.1771
2。3360
加速度
（m/s^2)
0。3586
1.2476
2。1839
2。5896
2。0403
8.5188
4.5654
—0.5615
—5。2088
—7。9416
作业3、如图1所示一3层钢筋混凝土框架结构FR3，各层层高均为4。2m.设计荷载:楼面恒荷5.0kN/m2，楼面活荷载2。0 kN/m2 ，屋面恒荷载5。0 kN/m2，屋面活荷载2。0kN/m2，结构构件尺寸和配筋见表1，该结构遭受一地震动作用。使用非线性静力分析方法或时程分析方法求取该地震动作用下此结构（图1右图FR3-A轴平面框架.不考虑楼板贡献。）的顶点最大侧向位移值（参考使用的软件：OpenSees、Abaqus、Ansys、Sap2000、Idarc等）。（20分左右）
地震动:1940 El Centro，N-S分量。
要 求：在此模板内完成，A4纸打印.

图1 3层结构的平立面布置图
表1a 框架结构柱的截面尺寸及配筋
结构编号

材料
尺寸(mm×mm）
主筋面积(mm2）/箍筋
混凝土
钢筋
（主筋）
边柱
中柱
边柱
中柱
FR3
1～3层
C30
HRB400
400×400
400×400
1608/φ8＠100
1608/φ8＠100
表1b 框架结构梁的截面尺寸及配筋
结构编号

材料
尺寸（宽×高）
（mm×mm）
主筋面积（mm2）/箍筋
混凝土
钢筋
（主筋）
边跨
中跨
跨中
支座
跨中
支座
FR3
1～3层
C30
HRB400
250×500
942/φ8＠200
1140/φ8＠100
942/φ8@200
1140/φ8＠100
利用SAP2000建立模型如下图所示，图中的grid point是所要求的关键点。
输入结构模型的物理数据，得到在El centro地震作用下框架侧移变形，如下图。
单位:m
表3 SAP2000输出的框架关键点位移表
TABLE： Joint Displacements

Joint
OutputCase
CaseType
StepType
U1
U2
U3
Text
Text
Text
Text
m
m
m
1
COMPOSITE
Combination
Max
0
0
0
1
COMPOSITE
Combination
Min
0
0
0
16
COMPOSITE
Combination
Max
0。193169
0。000101
0.000606
16
COMPOSITE
Combination
Min
—0.202617
—0。000168
—0。001481
表中U1、U2、U3分别表示所建模型的整体坐标系的x轴、y轴、z轴的正方向.由上述可得,所求FRA—3框架的顶点最大侧移量是0.2026m。