八层框架的地震响应计算和人工波生成的matlab实现及所需曲线图的自动存储.docx
《八层框架的地震响应计算和人工波生成的matlab实现及所需曲线图的自动存储.docx》由会员分享，可在线阅读，更多相关《八层框架的地震响应计算和人工波生成的matlab实现及所需曲线图的自动存储.docx（38页珍藏版）》请在冰点文库上搜索。
一、作业概况
结构基本参数：
层间剪切型结构，采用Rayleigh阻尼，第一、第二阶阻尼比分别取3%、5%。
图1结构基本形状
表1各层集中质量（105kg）
层号
1
2
3
4
5
6
7
8
质量
3.40
3.40
3.20
3.20
2.80
2.80
2.70
2.60
表2各层层间刚度（×108N/m）
层号
1
2
3
4
5
6
7
8
层间刚度
2.00
2.00
1.80
1.80
1.80
1.80
1.60
1.60
二、频率及振型计算
根据层间模型的假定，可以建立结构的质量矩阵以及刚度矩阵如下。
根据上面求得的质量、刚度矩阵，即可求解特征方程：
（1）
求解自振频率以及阵型向量已经演变成为典型的求解矩阵特征值以及特征向量的问题，特征值即为圆频率，特征向量即为振型向量。
根据式
（1）利用matlab编程计算，求解矩阵的特征值以及特征向量，进而可以得到结构前八阶的自振圆频率、自振频率、自振周期如表2.1所示。
表2.1结构前八阶振型自振（圆）频率及周期
振型
1
4.7605
0.7577
1.3198
2
13.3620
2.1266
0.4702
3
21.5638
3.4320
0.2914
4
29.1249
4.6354
0.2157
5
36.1617
5.7553
0.1738
6
41.2494
6.5650
0.1523
7
44.8736
7.1418
0.1400
8
48.0355
7.6451
0.1308
而上述特征值问题求得的幅值向量实质即为振型向量，在下文中记为，振型矩阵记为。
振型向量经过归一化处理之后如表2.2所示。
表2.2结构前八阶振型
层号
1阶振型
2阶振型
3阶振型
4阶振型
5阶振型
6阶振型
7阶振型
8阶振型
1
0.175
-0.482
-0.784
0.830
-0.970
0.714
-0.703
0.144
2
0.344
-0.817
-0.948
0.463
0.216
-0.637
1.000
-0.277
3
0.516
-0.914
-0.298
-0.687
1.000
-0.091
-0.912
0.463
4
0.668
-0.721
0.599
-0.801
-0.541
0.730
0.440
-0.696
5
0.793
-0.299
1.000
0.293
-0.824
-0.658
0.216
1.000
6
0.890
0.206
0.678
1.000
0.569
-0.305
-0.685
-0.893
7
0.963
0.710
-0.236
0.311
0.834
1.000
0.715
0.584
8
1.000
1.000
-0.965
-0.823
-0.742
-0.567
-0.315
-0.212
根据表2.2，可以绘制结构的前八阶振型图，如图2.1所示。
图2.1结构前八阶振型图
三、采用振型分解法进行地震时程计算
多自由度体系结构的位移反应可以表示为：
（2）
其中为广义坐标，为广义坐标向量。
故线性结构的动力方程可以表示为式（3）所示：
（3）
式中，为振型参与系数，满足
而阻尼矩阵采用瑞利阻尼假定，满足：
（4）
其中
利用振型的正交性，即：
（5）
在（3）式两端左乘，并与（4）、（5）联立，可得：
（6）
又因为之间满足：
（7）
（10）、（11）联立，两边同除，并记，可得：
（8）
至此，完成了对多自由度耦合的动力方程的解耦，形成了若干单自由度体系动力方程，可以据此利用时程分析方法进行各单自由度体系的时程计算。
下面选用了纽马克-法进行单自由度体系的时程计算。
在该方法中假定：
（9）
其中、以及分别表示对应第i个广义坐标在第k个时间点的相对位移、速度和加速度；而、以及则分别表示对应增量。
和为参数，分别取1/2和1/6。
又根据式（8），其增量形式为：
（10）
所以将（9）带入（10）中有：
（11）
在利用matlab编制程序时，记，
所以有：
（12）
即当已知地震加速度时程以及上一时刻单自由度体系加速度、速度以及位移时，利用式
（1）以及（5）即可求得体系在下一时刻的加速度、速度以及位移。
不断递推计算最终可以求得体系在时域的地震响应。
计算前，假定在初始0时刻，有：
（13）
利用matlab编制程序（程序见附），求得了八个广义坐标的地震反应时程。
采用此方法计算得到的位移时程为相对位移时程。
因此为了得到结构顶部的绝对位移时程，还需要计算出场地位移时程后，与顶部相对位移时程相加。
由于地震仪在记录地震动时记录纸的蛇行运动和放大器的不稳定等，记录的零线会产生很小的摇摆错位，因此在时程计算之前需要对地震波进行基线校正。
采用线性修正的方法，地震波位移、速度、加速度的修正值如式（14）所示。
（14）
其中，系数由下面的式子计算：
（15）
经过上述过程编程计算后，结构顶部位移时程如图3.1所示。
图3.1EL-Centro波作用下结构顶层位移时程
各层层间位移时程如图3.2所示。
图3.2EL-Centro波作用下结构各层层间位移时程
底层层间剪力时程如图3.3所示。
图3.3EL-Centro波作用下结构底层层间剪力时程
各层的层间位移与层间剪力绝对值包络图如图3.4所示。
（a）层间位移绝对值包络图（m）
（b）层间剪力绝对值包络图（kN）
图3.4EL-Centro波作用下结构各层层间位移与层间剪力绝对值包络图
下面对上述计算程序及结果的正确性进行简要的验证。
验证采用的方法为振型分解反应谱法，即利用该方法得到结构在EL-Centro波作用下各层层间剪力的最大值，与上文中时程计算得到的剪力绝对值包络图进行对比。
根据相关规范，在罕遇地震作用下结构阻尼取5%计算得到EL-Centro波NS向的加速度反应谱，反应谱图如图3.5所示。
图3.5EL-Centro波在5%阻尼比下绝对加速度反应谱
将结构前八阶振型的周期域反应谱进行比对（利用线性内插法），可以得到相应的结构最大绝对加速度。
结构的前八阶振型以及振型参与系数、最大绝对加速度通过计算列于表3.1中。
表3.1振型分解反应谱法所需结构参数
振型
1
1.3198
1.3033
2.054
2
0.4702
-0.4600
7.926
3
0.2914
-0.2548
6.696
4
0.2157
0.1673
6.343
5
0.1738
0.1090
6.864
6
0.1523
0.0852
5.229
7
0.1400
-0.0610
6.561
8
0.1308
0.0145
7.334
所以利用表3.1中的数据，利用下式可以求得各振型各质点的最大地震作用（其中j表示振型编号，i表示由下自上的质点号）。
（16）
列表如3.2所示。
表3.2各振型各质点的最大地震作用（单位：
kN）
j
i
1
2
3
4
5
6
7
8
1
159.5
597.2
427.8
281.9
203.1
89.1
75.9
4.0
2
312.8
1013.1
517.4
157.3
-45.3
-79.5
-108.1
-7.7
3
469.8
1133.6
162.5
-233.2
-209.4
-11.3
98.5
12.8
4
607.9
894.2
-326.7
-272.0
113.3
91.1
-47.6
-19.2
5
721.4
371.1
-545.9
99.4
172.6
-82.1
-23.4
27.6
6
809.6
-255.2
-370.2
339.7
-119.2
-38.0
74.0
-24.7
7
876.6
-880.0
128.7
105.8
-174.7
124.7
-77.3
16.1
8
910.1
-1239.6
526.6
-279.5
155.3
-70.7
34.0
-5.9
因此，通过表3.2可以得到各个振型下各层间剪力的最大值，如表3.3所示。
表3.3各振型各层间剪力的最大值（单位：
kN）
j
i
1
2
3
4
5
6
7
8
1
4867.8
1634.4
520.3
199.3
95.7
23.4
26.2
3.1
2
4708.3
1037.2
92.5
-82.7
-107.4
-65.7
-49.7
-0.9
3
4395.5
24.1
-425.0
-240.0
-62.1
13.8
58.3
6.8
4
3925.7
-1109.5
-587.5
-6.8
147.3
25.1
-40.2
-6.0
5
3317.8
-2003.7
-260.8
265.3
34.0
-66.0
7.4
13.2
6
2596.4
-2374.8
285.1
165.9
-138.6
16.1
30.8
-14.4
7
1786.8
-2119.6
655.3
-173.7
-19.4
54.1
-43.3
10.3
8
910.1
-1239.6
526.6
-279.5
155.3
-70.7
34.0
-5.9
由于各振型下层间剪力的最大值并不一定同时出现，因此采用平方和开方的方式，求得地震波作用下层间剪力绝对值的最大值，如表3.4所示。
表3.4各层间剪力的最大值（单位：
kN）
编号
1
2
3
4
5
6
7
8
5166.0
4824.7
4423.4
4124.5
3894.5
3537.0
2854.8
1658.6
从表3.4可以得到利用振型分解反应谱法求得的在ElCentro波NS向作用下各层间剪力绝对值的包络图，如图3.6（a）所示，3.6（b）为上文中求得的时程计算得到的各层间剪力绝对值的包络图。
图3.6两种方法计算EL-Centro波作用下各层层间剪力绝对值包络图
通过对比可以发现时程计算的结果与振型分解反应谱法的计算结果在底层吻合的很好，在结构高程的差距相对较大，但整体上两者的结果是在同一数量级，且差值合理的，最大差值为26.9%。
考虑到振型分解反应谱法本身上是一种为设计服务求得最大地震作用的近似参考方法，所以通过与本作业中编程进行时程计算的结果相比，在一定程度上从侧面印证了前文时程计算结果的合理性与正确性。
四、人工地震波合成
本文采用经典谱表达方法来进行强震加速度模拟。
经典谱表达方法如式（17）所示。
（17）
式中：
为时间包络函数，表达式如下所示：
（18）
相关参数选择按照表4.1选取：
表4.1时间包络函数的参数值及持时
为主振平稳段段首时间，对设计的两条人工波假定分别对应二类及三类场地的第一组地震波，可分别取0.8s以及1.2s；
为主振平稳段段尾时间，对两条人工波分别取7s以及9s；
为衰减系数，对两条人工波分别取0.35s-1以及0.25s-1。
式（17）中的可由下式计算。
（19）
即为频域样本点，为对应样本点的相位，为所选用的功率谱。
满足：
（20）
其中为截断频率，为频率下界，为频域的样本点数目。
所以，在上式的基础上，需要功率谱密度函数模型以及相位谱，才能够生成人工地震波。
而根据作业要求，选用功率谱密度函数模型Clough-Penzien谱模型，其表达式为：
（21）
式中相关参数选择按照表4.2选取：
表4.2Clough-Penzien模型场地土参数设计值
为场地卓越周期，对设计的两条人工波假定分别对应二类及三类场地的第一组地震波，故分别取20.94与13.96；
为二次过滤参数，取为0.15（一般可取0.1-0.2）；
为场地地基阻尼，对两条地震波分别取为0.72以及0.8；
为二次过滤参数，取=；
为谱强度因子，在本文中，为使两条人工地震波的加速度峰值在0.4g左右，经过调整试算，分别取0.025、0.032。
地震动持时采用相对分数持时，即地震动加速度记录最后一次达到峰值加速度0.05倍时所对应的时间。
（上述参数的选择参考文献：
张猛,张哲,李天.与规范反应谱相对应Clough-Penzien模型参数研究[J].世界地震工程,2007,01:
56-60.）
两条地震波对应的Clough-Penzien模型如图4.1所示。
（a）人工波1
（b）人工波2
图4.1人工波Clough-Penzien模型功率谱图
相位谱则采用随机方式生成随机相位谱。
故综上利用matlab编制人工地震波合成程序（程序见附），可以得到两条人工地震波（地震波数据见附），两条波的加速度时程如图4.2所示。
（a）人工波1
（b）人工波2
图4.2人工波加速度时程
五、人工波时程计算结果
将四中人工生成的两条地震波用前文的方法进行时程计算，得到了响应的计算结果。
人工波1作用下，结构顶层位移时程如图5.1所示。
图5.1人工波1作用下结构顶层位移时程
人工波1作用下，结构各层层间位移时程如图5.2所示。
图5.2人工波1作用下结构各层层间位移时程
在人工波1作用下，结构底层剪力时程如图5.3所示。
图5.3人工波1作用下结构底层剪力时程
在人工波1作用下，各层的层间位移与层间剪力绝对值包络图如图5.4所示。
（a）层间位移绝对值包络图
（b）层间剪力绝对值包络图
图5.4人工波1作用下结构各层层间位移与层间剪力绝对值包络图
人工波2作用下，结构顶层位移时程如图5.5所示。
图5.5人工波2作用下结构顶层位移时程
人工波2作用下，结构各层层间位移时程如图5.6所示。
图5.6人工波2作用下结构各层层间位移时程
在人工波2作用下，结构底层剪力时程如图5.7所示。
图5.7人工波2作用下结构底层剪力时程
在人工波2作用下，各层的层间位移与层间剪力绝对值包络图如图5.8所示。
（a）层间位移绝对值包络图
（b）层间剪力绝对值包络图
图5.8人工波2作用下结构各层层间位移与层间剪力绝对值包络图
1、程序在2013Ra上可以正常运行
2、运行时需要把‘ThreeKindsWaves’复制到当前工作路径下
3、文件的存储地址可以根据自身情况进行修改，或者不需要可以直接注释
%利用解特征方程的方法得到结构的8阶自振周期和振型，并归一化，然后利用循环产生阵型图，保存在文件‘阵型图’内
tic%计时开始
clearall;%清除变量
digits（4）%设定精度
M=10^5*diag（[3.43.43.23.22.82.82.72.6]）;%输入质量矩阵
K=10^8*（diag（[43.83.63.63.63.43.21.6]）+diag（[-2-1.8-1.8-1.8-1.8-1.6-1.6],1）+...
diag（[-2-1.8-1.8-1.8-1.8-1.6-1.6],-1））;%输入刚度矩阵
[fai,omiga]=eig（M^-1*K）;%利用matlab内置函数解特征方程，周期平方存为omiga，振型存为fai
omiga=omiga.^0.5;%得到周期
transfer=fai（:
8）;
fai（:
8）=fai（:
7）;
fai（:
7）=transfer;
file_id='C:
\Users\lenovo\Desktop\第二次作业用图\振型图\';
mkdir（file_id）;
fori=1:
8%归一化振型并画图
[val,poi]=max（abs（fai（:
i）））;
fai（:
i）=fai（:
i）/fai（poi,i）;
f（i）=omiga（i,i）/（2*pi）;
figure（'color','white'）;%建立图形窗口，背景白色
axisequal
axisoff
holdon
plot（[00],[0,4],'k-'）
x=[0fai（:
i）'];
y=linspace（0,4,9）;
yi=0:
0.001:
4;
pp=interp1（y,x,yi,'cubic'）;%利用内置插值函数使连线光滑
plot（pp,yi,'k--'）
plot（fai（:
i）,0.5:
0.5:
4,'o','MarkerFaceColor','k','MarkerSize',20）%控制点的形状和大小代表集中质量
text（fai（:
i）+0.2,0.6:
0.5:
4.1,num2str（fai（:
i）））%标注振型数值
forj=1:
8
plot（[0fai（j,i）],[0.5*j0.5*j],'k-'）%画相关的辅助线
end
name=strcat（'第',num2str（9-i）,'阶振型图'）;
title（name,'FontSize',16,'position',[0-0.3]）%标注图名
F=getframe（figure（i））;
imwrite（F.cdata,[file_id,strcat（name,'.png'）]）;%保存振型图片
end
f=f';
T=1./f;%输出周期和频率
Toc
%通过解特征方程得到振型和周期，并进行归一化，为下面进行时程分析做铺垫。
读取Excel表格存储的地震动数据，
%利用纽马克方法进行时程分析得到地面速度和地面位移，利用相关公式对数据进行基线校正得到更为合理的数据。
%利用纽马克方法进行时程分析得到结构各阶振型的响应，然后通过公式X=A*u得到结构总的位移响应，然后利用
%结构力学相关知识得到剪力和剪力包络图。
tic;clearall;clf;%计时开始，并清楚其他无关变量
M=10^5*diag（[3.43.43.23.22.82.82.72.6]）;%输入质量矩阵
K=10^8*（diag（[43.83.63.63.63.43.21.6]）+diag（[-2-1.8-1.8-1.8-1.8-1.6-1.6],1）+...
diag（[-2-1.8-1.8-1.8-1.8-1.6-1.6],-1））;%输入刚度矩阵
[fai,w]=eig（inv（K）*M）;%特征值问题求解
omiga=1./diag（w.^0.5）;%得到周期
fori=1:
8%归一化振型
[val,poi]=max（abs（fai（:
i）））;
fai（:
i）=fai（:
i）/fai（poi,i）;
end
gama=fai'*M*ones（8,1）./diag（fai'*M*fai）;%计算振型参与系数
kesai=（0.1329+0.0067*omiga.^2）/2./omiga;%计算各阶阻尼
file_id='C:
\Users\lenovo\Desktop\第二次作业用图\地震动响应图\';
mkdir（file_id）;%创建文件夹保存生成的所有图
A=xlsread（'ThreeKindsWaves','ELcentro-NS'）;%读取地震波数据
N=length（A（:
1））;
mint=A（2,1）-A（1,1）;%求时间间隔
T=（N-1）*mint;
allmotion=zeros（N,8）;%用来存储所有时刻的位移
allacceleration=zeros（N,8）;%用来存储所有时刻的加速度值
allacceleration（1,:
）=A（1,2）;%对0时刻加速度矩阵初始化
allvelocity=zeros（N,8）;%用于存储所有时刻速度
vg
（1）=0;
xg
（1）=0;
form=1:
N-1%计算场地速度、位移用于对加速度进行基线修正
vg（m+1）=vg（m）+A（m,2）*mint+（A（m+1,2）-A（m,2））*mint/2;
xg（m+1）=xg（m）+vg（m）*mint+A（m,2）*mint^2/2+mint^2*（A（m+1,2）-A（m,2））/6;
end
c1=28/13/T^2*（2*vg（N）-15/T^5*xg*（3*T*A（:
升级会员 