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

一、作业概况结构基本参数：层间剪切型结构，采用Rayleigh阻尼，第一、第二阶阻尼比分别取3%、5%。图1结构基本形状表1各层集中质量(105kg)层号12345678质量3.403.403.203.202.802.802.702.60表2各层层间刚度(×108N/m)层号12345678层间刚度2.002.001.801.801.801.801.601.601m2m3m4m5m6m7m8m1k2k3k4k5k6k7k8k()&&gxtto二、频率及振型计算根据层间模型的假定，可以建立结构的质量矩阵以及刚度矩阵如下。12345678000000000000000000000000000000000000000000000000000000003.4000000003.4000000003.2000000003.20000=00002.8000000002.8000000002.7000000002.6mmmmmmmmM510kg11121314151617182122232425262728313233343536373841424344454647485152535455565758616263646566676871727374757677788182838485868788kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkK84200000023.81.80000001.83.61.80000001.83.61.8000=10/0001.83.61.80000001.83.41.60000001.63.21.60000001.61.6Nm根据上面求得的质量、刚度矩阵，即可求解特征方程：20w-=KM（1）求解自振频率以及阵型向量已经演变成为典型的求解矩阵特征值以及特征向量的问题，特征值即为圆频率2w，特征向量即为振型向量。根据式（1）利用matlab编程计算，求解矩阵1MK-的特征值以及特征向量，进而可以得到结构前八阶的自振圆频率、自振频率、自振周期如表2.1所示。表2.1结构前八阶振型自振（圆）频率及周期振型/)（rads1（）fs（）Ts14.76050.75771.3198213.36202.12660.4702321.56383.43200.2914429.12494.63540.2157536.16175.75530.1738641.24946.56500.1523744.87367.14180.1400848.03557.64510.1308而上述特征值问题求得的幅值向量实质即为振型向量，在下文中记为j，振型矩阵记为Φ。振型向量经过归一化处理之后如表2.2所示。表2.2结构前八阶振型层号1阶振型2阶振型3阶振型4阶振型5阶振型6阶振型7阶振型8阶振型10.175-0.482-0.7840.830-0.9700.714-0.7030.14420.344-0.817-0.9480.4630.216-0.6371.000-0.27730.516-0.914-0.298-0.6871.000-0.091-0.9120.46340.668-0.7210.599-0.801-0.5410.7300.440-0.69650.793-0.2991.0000.293-0.824-0.6580.2161.00060.8900.2060.6781.0000.569-0.305-0.685-0.89370.9630.710-0.2360.3110.8341.0000.7150.58481.0001.000-0.965-0.823-0.742-0.567-0.315-0.212根据表2.2，可以绘制结构的前八阶振型图，如图2.1所示。图2.1结构前八阶振型图三、采用振型分解法进行地震时程计算多自由度体系结构的位移反应()xt可以表示为：1()()()njjjxtutut（2）其中()jut为广义坐标，()ut为广义坐标向量。故线性结构的动力方程可以表示为式（3）所示：1111()()()()nnnnjjjjjjjjgjjjjutututxtMCKMffff（3）式中，=TjjTjjMIMfff为振型参与系数，满足1=1njjjf而阻尼矩阵C采用瑞利阻尼假定，满足：=CMK（4）其中22221122212222210.034.760513.36200.0513.36204.7605=2=2=0.132913.36204.7605a22112222210.0513.36200.034.7605=2=2=0.00673813.36204.7605利用振型的正交性，即：0()0()TijTijijijMKffff（5）在（3）式两端左乘Tif，并与（4）、（5）联立，可得：()()()()iiiiiiiiiiiigutututxtMMKKMffffffff（6）又因为、MK之间满足：2=iiiMKff（7）（10）、（11）联立，两边同除iiMff，并记2=2iii，可得：22iiiiiiiguuux（8）至此，完成了对多自由度耦合的动力方程的解耦，形成了若干单自由度体系动力方程，可以据此利用时程分析方法进行各单自由度体系的时程计算。下面选用了纽马克-法进行单自由度体系的时程计算。在该方法中假定：k22,12ikiikikikikikuututuututut（9）其中iku、iku以及kiu分别表示对应第i个广义坐标在第k个时间点的相对位移、速度和加速度；而iku、iku以及iku则分别表示对应增量。和为参数，分别取1/2和1/6。又根据式（8），其增量形式为：22iiiiiiiguuux（10）所以将（9）带入（10）中有：222211222iiiikigkiikikikttuxututut（11）在利用matlab编制程序时，记22=12iiiStt，22=122ikikikiQututut所以有：ikigkSuxQ（12）即当已知地震加速度时程以及上一时刻单自由度体系加速度、速度以及位移时，利用式（1）以及（5）即可求得体系在下一时刻的加速度、速度以及位移。不断递推计算最终可以求得体系在时域的地震响应。计算前，假定在初始0时刻，有：000000igiiuxuu（13）利用matlab编制程序（程序见附），求得了八个广义坐标的地震反应时程。采用此方法计算得到的位移时程为相对位移时程。因此为了得到结构顶部的绝对位移时程，还需要计算出场地位移时程后，与顶部相对位移时程相加。由于地震仪在记录地震动时记录纸的蛇行运动和放大器的不稳定等，记录的零线会产生很小的摇摆错位，因此在时程计算之前需要对地震波进行基线校正。采用线性修正的方法，地震波位移、速度、加速度的修正值如式（14）所示。23012010111ˆ()()()261ˆ()()()2ˆ()()()ggggggxtxtatatxtxtatatxtxtaat（14）其中，系数01aa、由下面的式子计算：210231250/2281152()3213TyTaTaTayTytTttdtTT（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%阻尼比下绝对加速度反应谱将结构前八阶振型的周期域反应谱进行比对（利用线性内插法），可以得到相应的结构最大绝对加速度aS。结构的前八阶振型以及振型参与系数、最大绝对加速度通过计算列于表3.1中。表3.1振型分解反应谱法所需结构参数振型TsaS11.31981.30332.05420.4702-0.46007.92630.2914-0.25486.69640.21570.16736.34350.17380.10906.86460.15230.08525.22970.1400-0.06106.56180.13080.01457.334所以利用表3.1中的数据，利用下式可以求得各振型各质点的最大地震作用maxjiF（其中j表示振型编号，i表示由下自上的质点号）。max=jijjiajiFXSm（16）列表如3.2所示。表3.2各振型各质点的最大地震作用maxjiF（单位：kN）ji123456781159.5597.2427.8281.9203.189.175.94.02312.81013.1517.4157.3-45.3-79.5-108.1-7.73469.81133.6162.5-233.2-209.4-11.398.512.84607.9894.2-326.7-272.0113.391.1-47.6-19.25721.4371.1-545.999.4172.6-82.1-23.427.66809.6-255.2-370.2339.7-119.2-38.074.0-24.77876.6-880.0128.7105.8-174.7124.7-77.316.18910.1-1239.6526.6-279.5155.3-70.734.0-5.9因此，通过表3.2可以得到各个振型下各层间剪力的最大值，如表3.3所示。表3.3各振型各层间剪力的最大值maxjiV（单位：kN）ji1234567814867.81634.4520.3199.395.723.426.23.124708.31037.292.5-82.7-107.4-65.7-49.7-0.934395.524.1-425.0-240.0-62.113.858.36.843925.7-1109.5-587.5-6.8147.325.1-40.2-6.053317.8-2003.7-260.8265.334.0-66.07.413.262596.4-2374.8285.1165.9-138.616.130.8-14.471786.8-2119.6655.3-173.7-19.454.1-43.310.38910.1-1239.6526.6-279.5155.3-70.734.0-5.9由于各振型下层间剪力的最大值并不一定同时出现，因此采用平方和开方的方式，求得地震波作用下层间剪力绝对值的最大值，如表3.4所示。表3.4各层间剪力的最大值maxV（单位：kN）编号12345678maxV5166