八层框架的地震响应计算和人工波生成的matlab实现及所需曲线图的自动存储

一、 作业概况结构基本参数：层间剪切型结构，采用Rayleigh 阻尼，第一、第二阶阻尼比分别取3%、5%。

图1 结构基本形状表1 各层集中质量 ( 105kg)层号 1 2 3 4 5 6 7 8 质量3.403.403.203.202.802.802.702.60表2 各层层间刚度 (×108N/m)层号 1 2 3 4 5 6 7 8 层间刚度 2.002.001.801.801.801.801.601.60m m m m m m m m ()g x t二、 频率及振型计算根据层间模型的假定，可以建立结构的质量矩阵以及刚度矩阵如下。

12345678000000000000000000000000000000000000000000000000000000003.400000000 3.400000000 3.200000000 3.20000 =0000 2.800000000 2.800000000 2.700000000 2.6m m m m m m m m ⎛⎫ ⎪ ⎪ ⎪ ⎪⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎛⎫ ⎝M 510kg ⎪⎪⎪⎪⎪⨯⎪⎪⎪⎪⎪⎪⎭ 11121314151617182122232425262728313233343536373841424344454647485152535455565758616263646566676871727374757677788182838485868788k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k ⎛ =⎝K 8420000002 3.8 1.8000000 1.8 3.6 1.8000000 1.8 3.6 1.8000 =10/000 1.8 3.6 1.8000000 1.8 3.4 1.6000000 1.6 3.2 1.6000000 1.6 1.6N m ⎫⎪⎪⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭-⎛⎫ ⎪-- ⎪ ⎪-- ⎪-- ⎪⨯ ⎪-- ⎪-- ⎪ ⎪-- ⎪ ⎪-⎝⎭根据上面求得的质量、刚度矩阵，即可求解特征方程：20K M（1）求解自振频率以及阵型向量已经演变成为典型的求解矩阵特征值以及特征向量的问题，特征值即为圆频率2，特征向量即为振型向量。

根据式（1）利用matlab编程计算，求解矩M K的特征值以及特征向量，进而可以得到结构前八阶的自振圆频率、自振频率、自阵1振周期如表2.1所示。

表2.1 结构前八阶振型自振（圆）频率及周期而上述特征值问题求得的幅值向量实质即为振型向量，在下文中记为φj，振型矩阵记为Φ。

振型向量经过归一化处理之后如表2.2所示。

表2.2 结构前八阶振型根据表2.2，可以绘制结构的前八阶振型图，如图2.1所示。

图2.1 结构前八阶振型图三、 采用振型分解法进行地震时程计算多自由度体系结构的位移反应{}()x t 可以表示为：{}{}1()()()φ===Φ∑nj j j x t u t u t（2）其中()j u t 为广义坐标，{}()u t 为广义坐标向量。

故线性结构的动力方程可以表示为式（3）所示：1111()()()()γ====⎛⎫⎛⎫⎛⎫⎛⎫++=- ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭⎝⎭∑∑∑∑nnnnj j j j j j j j g j j j j u t u t u t x t M C K M （3） 式中，=γT j j T jjMIM为振型参与系数，满足{}1=1γ=∑nj jj而阻尼矩阵C 采用瑞利阻尼假定，满足： =αβ+C M K（4）其中22221122212222210.03 4.760513.36200.0513.3620 4.7605=2=2=0.132913.3620 4.7605a ζωωζωωαωω-⨯⨯-⨯⨯⨯--22112222210.0513.36200.03 4.7605=2=2=0.00673813.3620 4.7605ζωζωβωω-⨯-⨯⨯-- 利用振型的正交性，即：0()0()⎧=≠⎪⎨=≠⎪⎩Ti jT iji j i j M K（5）在（3）式两端左乘Ti，并与（4）、（5）联立，可得：()()()()()αβγ+++=-ii i ii ii i i i i i g u t u t u t x t M M K K M（6）又因为、M K 之间满足：2=ωi ii M K（7）（10）、（11）联立，两边同除i i M ，并记2=2αβωζω+i i i ，可得：22ζωωγ++=-i i i i i i i g u u u x（8）至此，完成了对多自由度耦合的动力方程的解耦，形成了若干单自由度体系动力方程，可以据此利用时程分析方法进行各单自由度体系的时程计算。

下面选用了纽马克-β法进行单自由度体系的时程计算。

在该方法中假定：k 22,12ηβ∆=∆+∆∆⎧⎪⎨∆=∆+∆+∆∆⎪⎩ik i ik ik ik ik ik u u t u t u u t u t u t （9）其中i k u 、ik u 以及k i u 分别表示对应第i 个广义坐标在第k 个时间点的相对位移、速度和加速度；而∆i k u 、∆ik u 以及∆ik u 则分别表示对应增量。

η和β为参数，分别取1/2和1/6。

又根据式（8），其增量形式为：22ζωωγ∆+∆+∆=-∆i i i i i i i g u u u x（10）所以将（9）带入（10）中有：()222211222i i i ik i gk i ik ik ik t t u x u t u t u t ηωζβωγωζω⎡⎤⎛⎫+∆+∆∆=-∆-∆+∆+∆ ⎪⎢⎥⎝⎭⎣⎦（11）在利用matlab 编制程序时，记22=12ηωβζω+∆+∆i i i S t t ，22=122ωζω⎛⎫∆+∆+∆ ⎪⎝⎭ik ik ik i Q u t u t u t所以有：i k i gk S u x Q γ∆=-∆-（12）即当已知地震加速度时程以及上一时刻单自由度体系加速度、速度以及位移时，利用式（1）以及（5）即可求得体系在下一时刻的加速度、速度以及位移。

不断递推计算最终可以求得体系在时域的地震响应。

计算前，假定在初始0时刻，有：000000=-∆==i g i i u x u u （13）利用matlab 编制程序（程序见附），求得了八个广义坐标的地震反应时程。

采用此方法计算得到的位移时程为相对位移时程。

因此为了得到结构顶部的绝对位移时程，还需要计算出场地位移时程后，与顶部相对位移时程相加。

由于地震仪在记录地震动时记录纸的蛇行运动和放大器的不稳定等，记录的零线会产生很小的摇摆错位，因此在时程计算之前需要对地震波进行基线校正。

采用线性修正的方法，地震波位移、速度、加速度的修正值如式（14）所示。

23012010111ˆ()()()261ˆ()()()2ˆ()()()g g g g g g xt x t a t a t x t x t a t a t x t x t a a t =-+=-+=-+ （14）其中，系数01a a 、由下面的式子计算：()()()210231250/2281152()3213T y T a T a Ta y T y t Tt t dt T T -=⎡⎤=--⎢⎥⎣⎦⎰（15）经过上述过程编程计算后，结构顶部位移时程如图3.1所示。

图3.1 EL-Centro波作用下结构顶层位移时程各层层间位移时程如图3.2所示。

图3.2 EL-Centro波作用下结构各层层间位移时程底层层间剪力时程如图3.3所示。

图3.3 EL-Centro波作用下结构底层层间剪力时程各层的层间位移与层间剪力绝对值包络图如图3.4所示。

(a)层间位移绝对值包络图（m）（b）层间剪力绝对值包络图（kN）图3.4 EL-Centro波作用下结构各层层间位移与层间剪力绝对值包络图下面对上述计算程序及结果的正确性进行简要的验证。

验证采用的方法为振型分解反应谱法，即利用该方法得到结构在EL-Centro波作用下各层层间剪力的最大值，与上文中时程计算得到的剪力绝对值包络图进行对比。

根据相关规范，在罕遇地震作用下结构阻尼取5%计算得到EL-Centro波NS向的加速度反应谱，反应谱图如图3.5所示。

图3.5 EL-Centro波在5%阻尼比下绝对加速度反应谱将结构前八阶振型的周期域反应谱进行比对（利用线性内插法），可以得到相应的结构最大绝对加速度S。

结构的前八阶振型以及振型参与系数、最大绝对加速度通过计算列于a表3.1中。

表3.1 振型分解反应谱法所需结构参数所以利用表3.1中的数据，利用下式可以求得各振型各质点的最大地震作用max ji F （其中j 表示振型编号，i 表示由下自上的质点号）。

max =ji j ji aj i F X S m γ（16）列表如3.2所示。

表3.2 各振型各质点的最大地震作用max ji F （单位：kN ）因此，通过表3.2可以得到各个振型下各层间剪力的最大值，如表3.3所示。

表3.3 各振型各层间剪力的最大值max ji V （单位：kN ）由于各振型下层间剪力的最大值并不一定同时出现，因此采用平方和开方的方式，求得地震波作用下层间剪力绝对值的最大值，如表3.4所示。

表3.4 各层间剪力的最大值max V （单位：kN ）从表3.4可以得到利用振型分解反应谱法求得的在El Centro 波NS 向作用下各层间剪力绝对值的包络图，如图3.6（a ）所示，3.6（b ）为上文中求得的时程计算得到的各层间剪力绝对值的包络图。

图3.6 两种方法计算EL-Centro 波作用下各层层间剪力绝对值包络图通过对比可以发现时程计算的结果与振型分解反应谱法的计算结果在底层吻合的很好，在结构高程的差距相对较大，但整体上两者的结果是在同一数量级，且差值合理的，最大差值为26.9%。

考虑到振型分解反应谱法本身上是一种为设计服务求得最大地震作用的近似参考方法，所以通过与本作业中编程进行时程计算的结果相比，在一定程度上从侧面印证了前文时程计算结果的合理性与正确性。

四、 人工地震波合成本文采用经典谱表达方法来进行强震加速度模拟。

经典谱表达方法如式（17）所示。

()1()cos()Ng i i i i X t A t g t ωφ==+⋅∑（17）式中：()g t 为时间包络函数，表达式如下所示：()()221112()2/ 0g = 1 c t t t t t t t t t t et t --⎧≤<⎪⎪≤<⎨⎪≥⎪⎩（18）相关参数选择按照表4.1选取：表4.1 时间包络函数的参数值及持时1t 为主振平稳段段首时间，对设计的两条人工波假定分别对应二类及三类场地的第一组地震波，1t 可分别取0.8s 以及1.2s ；2t 为主振平稳段段尾时间，对两条人工波分别取7s 以及9s ；c 为衰减系数，对两条人工波分别取0.35s -1以及0.25 s -1。