地震工程学大作业代码4求人工波位移时程
减隔震支座及黏滞阻尼器减震效果分析研究
２０２１年８月第１２卷第４期高 速 铁 路 技 术ＨＩＧＨＳＰＥＥＤＲＡＩＬＷＡＹＴＥＣＨＮＯＬＯＧＹＮｏ．４，Ｖｏｌ．１２Ａｕｇ．２０２１ 收稿日期：２０２０ １０ ２７作者简介：李晓波（１９９２ ），男，工程师。
基金项目：中国铁路设计集团有限公司科技开发课题（２０２０ＹＹ３４０６４１）引文格式：李晓波．减隔震支座及黏滞阻尼器减震效果分析研究［Ｊ］．高速铁路技术，２０２１，１２（４）：３４－３８．ＬＩＸｉａｏｂｏ．ＡｎａｌｙｓｉｓｏｆｔｈｅＤａｍｐｉｎｇＥｆｆｅｃｔｏｆＳｅｉｓｍｉｃＭｉｔｉｇａｔｉｏｎａｎｄＩｓｏｌａｔｉｏｎＢｅａｒｉｎｇｓａｎｄＦｌｕｉｄＶｉｓｃｏｕｓＤａｍｐｅｒｓ［Ｊ］．ＨｉｇｈＳｐｅｅｄＲａｉｌｗａｙＴｅｃｈｎｏｌｏｇｙ，２０２１，１２（４）：３４－３８．文章编号：１６７４—８２４７（２０２１）０４—００３４—０５ＤＯＩ：１０．１２０９８／ｊ．ｉｓｓｎ．１６７４－８２４７．２０２１．０４．００７减隔震支座及黏滞阻尼器减震效果分析研究李晓波（中国铁路设计集团有限公司， 天津３００３０８）摘 要：本文基于高烈度震区某（７２＋１２８＋７２）ｍ大跨连续梁桥，采用非线性时程分析方法，研究了铁路连续梁常用双曲面球型减隔震支座的减震效果及减震机理，并采用黏滞阻尼器控制结构过大的地震位移响应。
结果表明：（１）双曲面球型减隔震支座均可大幅减小结构地震内力响应，纵桥向墩底弯矩减震率在９０％左右，横桥向墩底弯矩减震率在８５％左右，墩梁相对位移呈非线性增大趋势；（２）在不影响双曲面球型减隔震支座减震效果的情况下，黏滞阻尼器可有效控制墩梁相对位移，但当地震高烈度较高且场地类别较差时，建议适当增大位移限值。
研究结果为高速铁路大跨连续梁桥减隔震设计提供了依据。
关键词：大跨度连续梁桥；双曲面球型减隔震支座；黏滞阻尼器；减震效果中图分类号：Ｕ４４２．５＋５ 文献标志码：Ａ ＡｎａｌｙｓｉｓｏｆｔｈｅＤａｍｐｉｎｇＥｆｆｅｃｔｏｆＳｅｉｓｍｉｃＭｉｔｉｇａｔｉｏｎａｎｄＩｓｏｌａｔｉｏｎＢｅａｒｉｎｇｓａｎｄＦｌｕｉｄＶｉｓｃｏｕｓＤａｍｐｅｒｓＬＩＸｉａｏｂｏ（ＣｈｉｎａＲａｉｌｗａｙＤｅｓｉｇｎＣｏｒｐｏｒａｔｉｏｎ，Ｔｉａｎｊｉｎ ３００３０８，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：Ｂａｓｅｄｏｎａ（７２＋１２８＋７２）ｍｌｏｎｇ ｓｐａｎｃｏｎｔｉｎｕｏｕｓｇｉｒｄｅｒｂｒｉｄｇｅｉｎａｈｉｇｈ ｉｎｔｅｎｓｉｔｙｅａｒｔｈｑｕａｋｅｚｏｎｅ，ｔｈｉｓｐａｐｅｒｓｔｕｄｉｅｓｔｈｅｄａｍｐｉｎｇｅｆｆｅｃｔａｎｄｍｅｃｈａｎｉｓｍｏｆｈｙｐｅｒｂｏｌｏｉｄｓｐｈｅｒｉｃａｌｓｅｉｓｍｉｃｉｓｏｌａｔｉｏｎｂｅａｒｉｎｇｓｃｏｍｍｏｎｌｙｕｓｅｄｉｎｒａｉｌｗａｙｃｏｎｔｉｎｕｏｕｓｇｉｒｄｅｒｓｂｙｕｓｉｎｇｎｏｎｌｉｎｅａｒｔｉｍｅｈｉｓｔｏｒｙａｎａｌｙｓｉｓｍｅｔｈｏｄ，ａｎｄｆｌｕｉｄｖｉｓｃｏｕｓｄａｍｐｅｒｓａｒｅｕｓｅｄｆｏｒｃｏｎｔｒｏｌｌｉｎｇｔｈｅｅｘｃｅｓｓｉｖｅｓｅｉｓｍｉｃｄｉｓｐｌａｃｅｍｅｎｔｒｅｓｐｏｎｓｅｏｆｔｈｅｓｔｒｕｃｔｕｒｅ．Ｔｈｅｒｅｓｕｌｔｓｓｈｏｗｔｈａｔ：（１）Ｈｙｐｅｒｂｏｌｏｉｄｓｐｈｅｒｉｃａｌｓｅｉｓｍｉｃｉｓｏｌａｔｉｏｎｂｅａｒｉｎｇｓｃａｎｇｒｅａｔｌｙｒｅｄｕｃｅｔｈｅｓｔｒｕｃｔｕｒａｌｉｎｔｅｒｎａｌｆｏｒｃｅｒｅｓｐｏｎｓｅ，ａｎｄｔｈｅｂｅｎｄｉｎｇｍｏｍｅｎｔｄａｍｐｉｎｇｒａｔｅａｔｔｈｅｐｉｅｒｂｏｔｔｏｍｉｓａｂｏｕｔ９０％ｉｎｔｈｅｌｏｎｇｉｔｕｄｉｎａｌｄｉｒｅｃｔｉｏｎｏｆｔｈｅｂｒｉｄｇｅ，ａｎｄａｂｏｕｔ８５％ｉｎｔｈｅｔｒａｎｓｖｅｒｓｅｄｉｒｅｃｔｉｏｎｏｆｔｈｅｂｒｉｄｇｅ，ａｎｄｔｈｅｒｅｌａｔｉｖｅｄｉｓｐｌａｃｅｍｅｎｔｏｆｐｉｅｒｓａｎｄｇｉｒｄｅｒｓｉｎｃｒｅａｓｅｓｎｏｎｌｉｎｅａｒｌｙ．（２）Ｆｌｕｉｄｖｉｓｃｏｕｓｄａｍｐｅｒｓｃａｎｅｆｆｅｃｔｉｖｅｌｙｃｏｎｔｒｏｌｔｈｅｒｅｌａｔｉｖｅｄｉｓｐｌａｃｅｍｅｎｔｏｆｐｉｅｒｓａｎｄｇｉｒｄｅｒｓｗｉｔｈｏｕｔａｆｆｅｃｔｉｎｇｔｈｅｄａｍｐｉｎｇｅｆｆｅｃｔｏｆｈｙｐｅｒｂｏｌｏｉｄｓｐｈｅｒｉｃａｌｂｅａｒｉｎｇｓ．Ｈｏｗｅｖｅｒ，ｉｔｉｓｓｕｇｇｅｓｔｅｄｔｈａｔａｎａｐｐｒｏｐｒｉａｔｅｌｙｌａｒｇｅｒｄｉｓｐｌａｃｅｍｅｎｔｌｉｍｉｔｓｈｏｕｌｄｂｅｔａｋｅｎｆｏｒｔｈｅｓｉｔｅｓｏｆｌｏｗｅｒｃａｔｅｇｏｒｙｉｎｈｉｇｈ ｉｎｔｅｎｓｉｔｙｅａｒｔｈｑｕａｋｅａｒｅａｓ．Ｔｈｅｒｅｓｕｌｔｓｐｒｏｖｉｄｅａｂａｓｉｓｆｏｒｔｈｅｓｅｉｓｍｉｃｉｓｏｌａｔｉｏｎｄｅｓｉｇｎｏｆｌｏｎｇ ｓｐａｎｃｏｎｔｉｎｕｏｕｓｇｉｒｄｅｒｂｒｉｄｇｅｓｏｆｈｉｇｈ ｓｐｅｅｄｒａｉｌｗａｙｓ．Ｋｅｙｗｏｒｄｓ：ｌｏｎｇ ｓｐａｎｃｏｎｔｉｎｕｏｕｓｇｉｒｄｅｒｂｒｉｄｇｅ；ｈｙｐｅｒｂｏｌｉｃｓｐｈｅｒｉｃａｌｓｅｉｓｍｉｃｉｓｏｌａｔｉｏｎｂｅａｒｉｎｇ；ｆｌｕｉｄｖｉｓｃｏｕｓｄａｍｐｅｒ；ｄａｍｐｉｎｇｅｆｆｅｃｔ１ 研究背景高速铁路大跨连续梁桥抗震设计目前主要还是以强度设计为主，当地震烈度过大或场地条件等导致设计方案难以实施时，则转向减隔震设计。
高层建筑地震作用时程分析法计算要点
建筑技术开发Building Technology Development 建筑设计Architectural Design 第48卷第5期2021年3月(1 )髙度不超过40 m 且以剪切变形为主并且质点和刚度沿高度分布均匀的结构(2)近似于单质点体系的结构(1 >不满足底部剪力法适用条件(2)高层建筑(3)质M 和刚度不对称不均匀的结构、超过 100 m 的髙层应采用考虑扭转耦联振动影响的方法（CQC )(4) _度大于24m 的楼盖、跨度大于12 m 的转换与连抹结构、悬挑长度大于5 m 的悬挑结构，竖向地震作用效应标准值| (丨）特别不规则的结构(2)甲类建筑(3 ) 7-9度时，髙规所列高度的乙丙类建筑 | (4)不满足高规所列高度的竖向不规则结构)(8 )平面投影尺度很大的空间结构（跨度大于120m 或长度大于300m 或悬臂大于40m ),7度III 和IV 类场地和8、9度时，用此法计算i f f B 级高度高层、混合结构和复杂高层建筑竖向）[静力法1—取结构或构件重力的一定百分数作为竖向地震作用地震作用计算方法J 1反应谱法按阵型分解反应谱法计算竖向地震作用f 百分数法规定结构或构件所受到的竖向地震作用为水平地震作用的某一百分数图1地震作用计算方法2.4 反应谱不同振型分解法采用的是考虑了震动强度与平均频谱特性的 设计谱，时程分析法全面反映了地震动强度、谱特征与持续时间三要素。
|(5) B 级高度的高层、混合结构和复杂高层建筑||(6)结构顶层取消部分墙.柱形成的空旷房间时1(7 >跨度大于24m 的楼盖，跨度大于12tn 的转换与连体结构.悬桃长度大于5m 的悬挑结构，竖向地震作用效应标准值高层建筑地震作用计算方法包括底部剪力法、振型分解 反应谱法（以下简称反应谱法）、时程分析法（以下简称时程 法）、弹塑性静力或动力分析法、静力法及百分数法。
其中底部剪力法和反应谱法是基本方法，时程分析法则是高层建筑 地震作用计算中有效的补充计算方法。
地震工程学
back
地震动观测仪器
两种仪器——地震仪与强震仪
• 二者不同点 • 二者共同点
拾震器
–运动方程——用拾震器的位移表示地面运动
放大器
记录器
mx’’+cx’+kx=-mxg’’
x’’+2EWx’+W2x=-xg’’
仪器
地震类型
按成因（构造EQ、火山EQ、陷落EQ、诱发EQ） 按发震位置（板边EQ、板内EQ） 按震源深度（浅源EQ、中源EQ、深源EQ） 按地震序列（主震余震型、震群型、单发型）
back
2.4 几个有关名词
震源 • 震中 • 震源深度 • 震源距 • 震中距
back
2.5 地震分布
20
25
30
35
40
t (s)
back
频谱特性
三种谱表述方法 简要评价
back
三种谱表述方法
傅立叶谱 功率谱 反应谱
back
a(t) (m/s2)
傅立叶谱
2.5 1.5 0.5 -0.5 -1.5 -2.5
0
5
10
15
20
25
30
35
40
t (s)
y(t)=sigma{Aisin(Wit+Qi)}
傅立叶变换
Ai
Qi
Wi
Wi back
反应谱
单自由度弹性体系的地震反应 反应谱的定义 反应谱的性质 反应谱的种类 反应谱的影响因素及规律
back
单自由度弹性体系的地震反应
单自由度弹性体系
p-4 人造地震动和地震动数值模拟
合 成 方 法
3、基于相位差谱的方法
相位差谱：指将相位谱在频率轴上等间距离散化后 两个相邻相位角的差，相位差随频率变化的函数称 为相位差谱。 研究表明，相位差谱的频数分布与加速度包络函数 形状相似，近似呈现为正态分布。 基于相位差谱的方法的特点即为考虑了相位差因素 及分布特点。 7
合 成 方 法
用三角级数之和构造一个近似的平稳高斯 过程，然后乘以强度包线，以得到非平稳的地 面运动加速度时程。
10
三角级数法的常用模型
11
三角级数法的步骤
12
三角级数法的优缺点
优点 实现方便，适应性强 缺点 （1）在迭代的过程中，可能会出现个别控制点 的反应谱不收敛于目标谱的情况，尤其是当反应谱 的控制点较多时容易出现。 （2） 时常看到合成的地震波虽然满足了目标谱 的要求，但是波形极差。例如：出现“大尖”或是 “平台型峰值”，或者出现时间强度包线明显截头 的情况。 （3）人造地震波与真实地震波主要在基本特征 上有相当的差距：频率性征和强度性征。
人造地震动和数值模拟地震动
Synthetic and Numerical Simulation on ）时程分析法是当前在地震工程中获得最广泛应用 的动力分析方法。在应用时程分析法时，要选择适当 的地震动作为输入。由于受地震观测的客观条件的限 制，取得的强震加速度记录’缺少’，不能满足工程中的 多方面要求。 （2）在随机振动理论中，一个困难是实际可以得到的、 而且又满足同一集系的地震动次数常常过少，使严格 的统计分析难以进行；且地震是稀有的自然事件，真 实地震动非人为所能制造产生。 （3）。。。
2
ln(Te ) d eM w 2
A1
Main pulse
大震下位移计算书shishenme
大震下位移计算书shishenme 大震下位移计算书是地震工程学中一个重要的工具，用于估计建筑物、结构物以及土体中的位移变形情况。
在地震工程中，地震是一种导致地壳变形和地表震动的自然现象，而位移计算则是预测和评估结构物对地震响应的关键步骤之一。
位移计算是地震工程中的一个关键问题，因为地震引发的地壳运动将会导致建筑物、结构物和土体中的位移变形。
这些位移变形不仅会对结构物本身造成影响，还可能会对周围环境和人员安全带来潜在威胁。
因此，准确地预测和评估结构物的位移变形情况是地震工程中重要的任务之一。
在位移计算中，地震动作为输入参数，通过数值模拟方法或解析方法，计算出结构物在地震荷载作用下的位移响应。
位移计算通常包括两个主要方面：一是计算建筑物或结构物的整体位移变形，二是计算土体的变形和位移。
这两个方面的计算相互关联，共同决定了建筑物和土体的相对位移。
对于建筑物或结构物来说，位移计算主要涉及其整体刚度和弹性模量的关系，以及地震荷载的传递和分布，通过求解结构物的动态方程或使用分析方法，可以得到结构物的位移响应。
针对复杂结构物，如桥梁、高层建筑等，通常需要使用数值模拟方法，如有限元分析，来模拟结构物的非线性响应，从而获得更准确的位移计算结果。
对于土体来说，位移计算涉及土体的饱和度、密度、剪切强度等参数，以及地震荷载的传递和分布。
通常使用土体模型来描述土体的力学行为，如弹性模型、弹塑性模型等。
通过求解土体的力学方程，可以得到土体的位移变形情况。
土体的位移计算对于地震工程中的土体液化、土体侧移等问题尤为重要，可以帮助工程师评估基础的稳定性和土体的变形情况。
在位移计算中，还需要考虑到地震荷载的时程特性和频率特性。
通常使用地震波记录来表示地震荷载，通过对地震波进行处理和分析，可以获得地震荷载的频谱特性、峰值加速度等参数。
在位移计算中，需要将地震荷载的时程输入到结构物或土体模型中，然后计算出相应的位移响应。
综上所述，大震下位移计算书是地震工程学中一个重要的工具，用于估计建筑物、结构物以及土体中的位移变形情况。
地震工程作业哈工大 MATLAB程序
地震工程大作业哈工大李金平弹性反应谱原创性声明，主程序由本人独立编写完成，支持各种检验。
所选地震动RSN2345CHICHI.AT2 peer索取号：2345:010203040506070 -0.25-0.2-0.15-0.1-0.050.050.10.150.20.25时间 (s)加速度(m/s2)位移反应谱局部放大12345678910周期T (s)位移 (m m )11.021.041.06 1.081.11.129.29.259.39.359.49.459.59.559.69.659.7周期T (s)位移 (m m )求比值123456789100.940.950.960.970.980.9911.011.021.031.04周期T (s)比值绝对加速度反应谱局部放大0.20.40.60.81 1.2 1.41.61.8200.10.20.30.40.50.60.70.80.91周期T (s)S a (m /s 2)-0.25-0.2-0.15-0.1-0.0500.050.10.150.20.250.150.20.250.30.350.40.45周期T (s)S a (m /s 2)求比值进行比较速度反应谱123456789100.940.950.960.970.980.9911.011.021.031.04周期T (s)比值1234567891000.010.020.030.040.050.060.070.080.090.1周期T (s)速度 (m /s )局部放大观察差异求比值进行比较。
00.020.040.060.080.10.12-4周期T (s)速度 (m /s )123456789100.70.80.911.11.21.3周期T (s)比值结论：对比：拟合效果非常好，短周期0s<T<0.2s误差较大，其中相对速度谱的误差较大，接近30%，加速度和位移误差都在3%以内，周期0.2s<T<10s各反应谱的误差都保持在1%以内，精度相同，猜想，两种计算程序应该用的同一种计算方法（newmark法）。
地震工程学-反应谱和地震时程波的相互转化matlab编程#(精选.)
地震工程学作业课程名称:地震工程学______ 指导老师:_______翟永梅_________ 姓名:史先飞________ 学号:1232627________一、地震波生成反应谱1 所取的地震波为Elcentro地震波加速度曲线，如图1所示。
图1 Elcentro地震波加速度曲线2 所调用的Matlab程序为：% ***********读入地震记录***********ElCentro;Accelerate= ElCentro(:,1)*9.8067;%单位统一为m和sN=length(Accelerate);%N 读入的记录的量time=0:0.005:(N-1)*0.005; %单位 s%初始化各储存向量Displace=zeros(1,N); %相对位移Velocity=zeros(1,N); %相对速度AbsAcce=zeros(1,N); %绝对加速度% ***********A,B矩阵***********Damp=0.02; %阻尼比0.02TA=0.0:0.05:6; %TA=0.000001:0.02:6; %结构周期Dt=0.005; %地震记录的步长%记录计算得到的反应，MaxD为某阻尼时最大相对位移，MaxV为某阻尼最大相对速度，MaxA某阻尼时最大绝对加速度，用于画图MaxD=zeros(3,length(TA));MaxV=zeros(3,length(TA));MaxA=zeros(3,length(TA));t=1;for T=0.0:0.05:6NatualFrequency=2*pi/T ; %结构自振频率DampFrequency=NatualFrequency*sqrt(1-Damp*Damp); %计算公式化简e_t=exp(-Damp*NatualFrequency*Dt);s=sin(DampFrequency*Dt);c=cos(DampFrequency*Dt);A=zeros(2,2);A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);A(1,2)=e_t*s/DampFrequency;A(2,1)=-NatualFrequency*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)/(NatualFrequency^2*Dt);d_3t=Damp/(NatualFrequency^3*Dt);B=zeros(2,2);B(1,1)=e_t*((d_f+Damp/NatualFrequency)*s/DampFrequency+(2*d_3t+1/NatualFrequency^2)*c)-2*d_3 t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)-1/NatualFrequency^2+2*d_3t;B(2,1)=e_t*((d_f+Damp/NatualFrequency)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/NatualFrequency^2 )*(DampFrequency*s+Damp*NatualFrequency*c))+1/(NatualFrequency^2*Dt);B(2,2)=e_t*(1/(NatualFrequency^2*Dt)*c+s*Damp/(NatualFrequency*DampFrequency*Dt))-1/(NatualF requency^2*Dt);for i=1:(N-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*NatualFrequency*Velocity(i+1)-NatualFrequency^2*Displace(i+1);endMaxD(1,t)=max(abs(Displace));MaxV(1,t)=max(abs(Velocity));if T==0.0MaxA(1,t)=max(abs(Accelerate));elseMaxA(1,t)=max(abs(AbsAcce));endDisplace=zeros(1,N);%初始化各储存向量，避免下次不同周期计算时引用到前一个周期的结果Velocity=zeros(1,N);AbsAcce=zeros(1,N);t=t+1;End% ***********PLOT***********close allfigure %绘制地震记录图plot(time(:),Accelerate(:))title('PEER STRONG MOTION DATABASE RECORD')xlabel('time(s)')ylabel('acceleration(g)')gridfigure %绘制位移反应谱plot(TA,MaxD(1,:),'-.b',TA,MaxD(2,:),'-r',TA,MaxD(3,:),':k')title('Displacement')xlabel('Tn(s)')ylabel('Displacement(m)')legend('ζ=0.02')Gridfigure %绘制速度反应谱plot(TA,MaxV(1,:),'-.b',TA,MaxV(2,:),'-r',TA,MaxV(3,:),':k') title('Velocity')xlabel('Tn(s)')ylabel('velocity(m/s)')legend('ζ=0.02')Gridfigure %绘制绝对加速度反应谱plot(TA,MaxA(1,:),'-.b',TA,MaxA(2,:),'-r',TA,MaxA(3,:),':k') title('Absolute Acceleration')xlabel('Tn(s)')ylabel('absolute acceleration(m/s^2)')legend('ζ=0.02')Grid3 运行的结果得到的反应谱图2 位移反应谱图3 速度反应谱图4 加速度反应谱一、反应谱生成地震波1所取的反应谱为上海市设计反应谱图5 上海市设计反应谱2反应谱取值程序为：%%规范反应谱取值程序参照01年抗震规范function rs_z=r_s_1(pl,zn,ld,cd,fz) %%%pl 圆频率,zn阻尼比,ld烈度,cd场地类型,场地分组fz %%%%烈度选择if ld==6arfmax=0.11;endif ld==7arfmax=0.23;endif ld==8arfmax=0.45;endif ld==9arfmax=0.90;end%%%%场地类别，设计地震分组选择if cd==1if fz==1Tg=0.25;endif fz==2Tg=0.30;endif fz==3Tg=0.35;endendif cd==2if fz==1Tg=0.35;if fz==2Tg=0.40;endif fz==3Tg=0.45;endendif cd==3if fz==1Tg=0.45;endif fz==2Tg=0.55;endif fz==3Tg=0.65;endendif cd==4if fz==1Tg=0.65;endif fz==2Tg=0.75;endif fz==3Tg=0.90;endend%%%%%%%%%ceita=zn; %%%%%阻尼比lmt1=0.02+(0.05-ceita)/8;if lmt1<0lmt1=0;endlmt2=1+(0.05-ceita)/(0.06+1.7*ceita); if lmt2<0.55lmt2=0.55;endsjzs=0.9+(0.05-ceita)/(0.5+5*ceita); %%%%%分段位置 T1 T2 T3T1=0.1;T2=Tg;T_jg=2*pi./pl;%%%% 第一段 0～T1if T_jg<=T1arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg;end%%%% 第二段 T1～T2if T1<T_jg&T_jg<=T2arf_jg=lmt2*arfmax;end%%%% 第三段 T2~T3if T2<T_jg&T_jg<=T3arf_jg=((Tg/T_jg)^sjzs)*lmt2*arfmax;end%%%% 第四段 T3～6.0if T3<T_jg&T_jg<=6.0arf_jg=(lmt2*0.2^sjzs-lmt1*(T_jg-5*Tg))*arfmax;end%%%% 第五段 6.0～if 6.0<T_jgarf_jg=(lmt2*0.2^sjzs-lmt1*(6.0-5*Tg))*arfmax;end%%%%%%反应谱值拟加速度值rs_z=arf_jg*9.8;end3生成人造地震波主程序：%%%主程序%%%%%%%%确定需要控制的反应谱Sa（T）(T=T1,...,TM)的坐标点数M，反应谱控制容差rc Tyz=[0.04:0.016:0.1,0.15:0.05:3.0,3.2:0.05:5.0];rc=0.06;nTyz=length(Tyz);ceita=0.035;%%%阻尼比：0.035for i=1:nTyzSyz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1); %%%%8度，2类场地，第1地震分组end%%%%%% 变换的频率差：2*pi*0.005(可以保证长周期项5s附近有5项三角级数)；%%%%频率变化范围 N1=30, 30*0.005*2*pi ;N2=3000, 5000*0.005*2*piplc=2*pi*0.005;pl=30*0.005*2*pi:0.005*2*pi:10000*0.005*2*pi;npl=length(pl);P=0.9; %%%保证率%%%%%%人造地震动持续时间40s，时间间隔：0.02sTd=40;dt=0.02;t=0:0.02:40;nt=length(t);%%%%%%% 衰减包络函数t1=8; %%%%上升段t2=8+24; %%%%%平稳段; 下降段则为40－32＝8sc=0.6; %%%%衰减段参数for i=1:ntif t(i)<=t1f(i)=(t(i)/t1)^2;endif t(i)>t1 & t(i)<t2f(i)=1;endif t(i)>=t2f(i)=exp(-c*(t(i)-t2));endend%%%%%%% 反应谱转换功率谱for i=1:nplSw(i)=(2*ceita/(pi*pl(i)))*r_s_1(pl(i),ceita,8,2,1)^2/(-2*log(-1*pi*log(P)/(pl(i)*Td))); Aw(i)=sqrt(4*Sw(i)*plc);end%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfai(i)=rand(1)*2*pi;for j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求,需要时程动力分析%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0)); dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=1;while max(rcrsp(:))>rc%%%%%循环体函数blxs=Syz./rspa;for xsxs=1:nplif 2*pi/pl(xsxs)<Tyz(1)blxs1(xsxs)=blxs(1);endfor sxsx=1:nTyz-1if (2*pi/pl(xsxs)>=Tyz(sxsx)) & (2*pi/pl(xsxs)<=Tyz(sxsx+1))blxs1(xsxs)=blxs(sxsx)+(blxs(sxsx+1)-blxs(sxsx))*(2*pi/pl(xsxs)-Tyz(sxsx))/(Tyz(sxsx+1)-Tyz(sxsx));endendif 2*pi/pl(xsxs)>Tyz(nTyz)blxs1(xsxs)=blxs(nTyz);endendAw=Aw.*blxs1;%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfor j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0)); dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=jsnum+1max(rcrsp(:))end%%%%%%% 最终的反应谱与规范谱%%%%%%%%%%%% response spectra of callidar%%%%%%% parameter%% Tjs=0.05:0.01:6;%% nTjs=length(Tjs);g=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0));dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i)/g;rspa_S(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1)/g;endsubplot(2,1,1);plot(t,at);subplot(2,1,2);plot(Tyz,rspa);hold on;plot(Tyz,rspa_S);4生成的人造地震波如图所示。
人造地震动实用数值方法
106
四川建筑
第 31 卷 3 期 2011． 06
· 工
程
结
构 ·
图1
k = 200 时的 EL － Centro 南北地震动波模拟
图2
k = 500 时的 EL － Centro 南北地震动波模拟
图3
k = 1000 时的 EL － Centro 南北地震动波模拟
图4
EL － centro 南北地震动波原始记录曲线
连接的节点应该在构造上予以特别加强 。 （ 5 ） 局部有限元分析 FEQ 计算 从 PMCAD 中截取一榀框架， 采用高精度平面有限元分 析程序 FEQ 进行二次应力分析， 得出加强部位的应力和配 筋。实际配筋时取 FEQ 和 SATWE 计算结果的较大值。 （ 6 ） 室外天桥计算 室外天桥计算采用同济大学 3D3S 软件进行分析与设 计。21 m 跨的人行天桥桁架上下弦为梁单元 ， 腹杆为杆单 元。计算结果表明： 桁架在各工况组合作用下 ， 最大竖向位 、 。 移 竖向挠度均满足规范要求
若函数 f（ x） 在满足 Dirichlet 条件下（ 在一个周期内连续 或只有有限个第一类间断点 ， 没有无穷间断点及振荡间断 点， 并且至多只有有限个极值点 ） 可表示成如下的三角级数 形式： a0 + 2 式中： an = bn = x dx ∫ f（ x） cos nπ l 1 n πx f（ x） sin dx l ∫ l
110
四川建筑
第 31 卷 3 期 2011． 06
（ 3） （ 4）
2 l
x dx ∫ f（ x） sin nπ l
0
l
㊆ g （ t） 在其历程区间上满足狄 由于被模拟的地震动函数 y 氏条件， 所以上述的三角级数肯定是收敛的 ， 但具体求值时 不能取无穷项， 一般取至前 500 项左右精度可以满足要求 ， 设取其前 k 或 n 项， 则有：
上海地震波-三向输入选取(说明)m
上海地区抗震设计输入地震时程说明（共8页）同济大学房结构工程与防灾研究所二〇一二年六月目录1 天然地震时程选取原则 (3)2 峰值调整 (3)3 频谱特性 (3)4 地震动持时 (3)5 人造地震动生成的方法 (3)6 目标反应谱的确定 (4)7 所选地震时程的基本信息 (4)8 地震时程反应谱与规范反应谱对比 (5)上海地区抗震设计输入地震时程说明1 天然地震时程选取原则天然地震动具有很强的随机性，随着输入地震波的不同结构的地震响应也会有很大的差异，故要保证时程分析结果的合理性，在选择地震波时必须遵循一定的原则。
一般而言，选择输入地震波时应以地震波的三要素（峰值、频谱特性、地震动持时）为主要考虑因素。
2 峰值调整地震波的峰值一定程度上反应了地震波的强度，因此要求输入结构的地震波峰值应与设防烈度要求的多遇地震或罕遇地震的峰值相当（峰值相当并非峰值相等，而是在峰值相近的情况下所选用地震波的反应谱与规范反应谱基本相符）。
3 频谱特性频谱是地面运动的频率成分及各频率的影响程度。
它与地震传播距离、区域、介质及结构所在的场地土性质有密切关系。
一般来说，在震中附近或岩石等坚硬场地土中，地震波中的短周期成分较多，在震中距较远或软弱场地土中，地震波的长期成分较多。
输入地震波的卓越周期应尽可能与拟建场地的特征周期一致，且在一定的周期段内与规范反应谱尽量接近。
对于天然地震记录而言，3个方向地震波同时都与规范反应谱很接近的条件是很难满足的，但应保证至少一个水平向地震波反应谱与规范反应谱基本吻合。
4 地震动持时地震持时也是结构破坏和倒塌的重要因素，工程实践中确定地震动持续时间的原则是：1）地震记录最强烈部分应包含在所选持续时间内，2）若对结构进行弹塑性地震反应分析（考虑累计损伤效应），持续时间可取长些。
另外，在截取地震波时尚需注意尽量在速度/位移零点处截断以尽量避免加速度积分时速度或位移的``漂移''现象。
人工地震动生成程序
clearclcclose all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fni=input('生成人工地震波-输入数据文件名(20041012):','s');fid=fopen(fni,'r');fs=fscanf(fid,'%f',1);%采样频率tu=fscanf(fid,'%f',1);%上升时间长度%上升时间包络线线形(1-直线、2-抛物线、3-指数曲线)iu=fscanf(fid,'%f',1);%上升时间包络线线形参数（只有指数曲线需要具体参数，其均为1）cu=fscanf(fid,'%f',1);ta=fscanf(fid,'%f',1);%持时时间长度td=fscanf(fid,'%f',1);%下降时间长度%下降时间包络线线形(1-直线、2-抛物线、3-指数曲线)id= fscanf(fid,'%f',1);%下降时间包络线线形（只有抛物线，指数曲线需要具体参数，其余为1）cd=fscanf(fid,'%f',1);dp=fscanf(fid,'%f',1);%阴尼比值p=fscanf(fid,'%f',1);%概率系数（一般可取P=0.85）nn=fscanf(fid,'%f',1);%迭代次数fno=fscanf(fid,'%f',1);%输出数据文件名%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对目标反应谱取值x=fscanf(fid,'%f',[2,inf]);%反应谱频率和幅值数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tatus=fclose(fid);%计算生成地震波的数据长度tl=tu+ta+td;%计算生成地震波的数据长度nt=round(fs*tl+1);%大于并最接近nt的2的幂次方为FFT长度nfft=2^nestpow2(nt)%计算频率间隔（Hz）df=fs/nfft%定义反应谱的离散频率向量f=0:df:(nfft/2-1)*df%计算时间间隔（s）dt=1/fs;%定义的离散时间向量t=0:dt:(nt-1)*dt%生成0到2PI的随机数为随机相位g=rand(1,nfft/s)*2*pi;%建立时间包络线%建立与地震波长度相同元素为1的向量en=ones(1,nt);%上升时间阶段%确定上升时间段的长度l=round(tu*fs)+1%产生上升时间段的包络线数组元素switch iucase 1 %直线en(1:l)=linspace(0,1,1);% y = linspace(a,b,n) generates a row vector y of n points linearlyspaced between and including a and b.case 2 %抛物线a=0:l-1;en(1:l)=(a/(l-1)).^2;case 3 %指数曲线a=0:l-1;en(1:l)=1-exp(-cu*a/(l-1));end%持续时间阶段%确定0时刻到持续时间结束时刻时间段的长度m=round((tu+ta)*fs)+1;%下降时间阶段%产生下降时间段的包络线数组元素switch idcase 1 %直线en(m:nt)=linspace(1,0,nt-m+1);case 2 %抛物线a=0:nt-m;en(m:nt)=1-cd*(a/(nt-m)).^2;case 3 %指数曲线a=0:nt-m;en(m:nt)=exp(-cd*a/(nt-m));end%按线性插值建立目标反应谱离散数据%按目标反应谱的长度生成元素为0的向量a0=zeros(x(1,:));%取目标反应谱数据的长度n=length(x(1,:));%四舍五入取整求反应谱最大频率对应数组元素的下标nb=round(x(1,n)/df)+1;for k=1:n-1%四舍五入取整求反应谱前一个频率数据对应数组元素的下标 l=round(x(1,K)/df)+1;%四舍五入取整求反应谱后一个频率数据对应数组元素的下标m=round(x(1,K+1)/df)+1;%线性插值产生前后两个频率数据间的反应谱数组元素 a0(1:m)=linspace(x(2,k),x(2,k+1),m-l+1)end%根据目标反应谱计算对应的近似功率谱a1=a0;s=zeros(1,nfft/2);k=nb:ne;s(k)=2*dp/(pi.*(a1(k).^2)./f(k)./(-2*log(-log(p)*pi/tl)./f(k)));%将功率谱转换成傅里叶幅值谱b1=sqrt(4*df*s)*nfft/2;%定义元素为0的反谱传递函数矩阵hf=zeros(ne,nfft);%计算加速度反应谱传递函数矩阵for j-0:ne-1w=2*pi*df*jwd=w*sqrt(1-dp*dp);e=exp(-t.*W*dp);a=t.*wd;s=sin(a)>*((1-2*dp*dp)/(1-dp*dp));c=cos(a).*(2*dp/sqrt(1-dp*dp));%计算加速度反应谱的脉冲响应函数向量h=wd*e.*(s+c)/fs;%通过FFT变换求加速度反应谱传递函数向量hf(j+1,:)=fft(h,nfft);endmm=nn%进行生成人工地震波迭代计算%100为最大迭代次数for k=1:100%将幅值谱和相位谱转化为实部和虚部c=b1.*exp(i*g);%将正负圆频率傅里叶谱向量组合成一仙向量d=[c,c(nfft/2:-1:1)];%IFFT变换，并取变换结果实部为生成的地震波e=ifft(d,nfft);%给生成的地震波加上强度包络线y=en.*real(e(1:nt));%计算反应谱%对生成的地震波进行FFT变换yf=fft(y,nfft);for j=1:ne%用地震波FFT变换结果和反应谱传递函数的乘积的逆变换做卷积运算 d=ifft(yf.*hf(j,:),nfft);%求各频率对应地震的最大响应al(j)=max(real(d(1:nt)));end%如果达到指定的迭代次数显于图形if k==mmsubplot(2,1,1);%m 同时显示生成的地震波的强度包络线plot(t,y,t,en,t,-en);xlabel('时间（s）');ylabel('加速度（g）');grid onsublpot(2,1,2);%同时显示期望反应谱与反应谱计算谱、l=1:neplot(f(l),a0(l),':'f(l),a1(l));xlabel('频率（Hz）');ylabel('加速度（g）');legend('目标谱','计算谱');grid on;ig=input('继续迭代次数[取值1-9,否则退出]:');if ig>0&ig<10 %如果输入数字是1-9mm=mm+igelsebreak;endendc=bl%期望谱与计算谱的比值来修改傅里叶值谱j=nb:ne;bl(j)=c(j).*a0(j)./a1(j);end%打开文件输入人工地震动数据fid=fopen(fno,'W');for K=1:nt%每一行输出两个实型数据,t为时间,y为人工地震动信号值 fprintf(fid,'%f%f\n',t(k),y(k));endstatus=fclose(fid);。
八层框架的地震响应计算和人工波生成的matlab实现及所需曲线图的自动存储
一、 作业概况结构大体参数：层间剪切型结构，采纳Rayleigh 阻尼，第一、第二阶阻尼比别离取3%、5%。
图1 结构大体形状表1 各层集中质量 ( 105kg)层号12 3 4 5 6 7 8 质量表2 各层层间刚度 (×108N/m)层号 1 2 3 4 5 6 7 8 层间刚度m 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 ⎫⎪⎪⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭-⎛⎫ ⎪-- ⎪ ⎪-- ⎪-- ⎪⨯ ⎪-- ⎪-- ⎪ ⎪-- ⎪ ⎪-⎝⎭依照上面求得的质量、刚度矩阵，即可求解特点方程： 20KM （1）求解自振频率和阵型向量已经演变成为典型的求解矩阵特点值和特点向量的问题，特点值即为圆频率2，特点向量即为振型向量。
八层框架的地震响应计算和人工波生成的matlab实现及所需曲线图的自动存储
一、 作业概况结构基本参数：层间剪切型结构，采用Rayleigh 阻尼，第一、第二阶阻尼比分别取3%、5%。
图1 结构基本形状表1 各层集中质量 ( 105kg)层号 12345678质量表2 各层层间刚度 (×108N/m)层号 1 2 3 4 5 6 7 8 层间刚度m 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，特征向量即为振型向量。
地震工程学反应谱和地震时程波的相互转化编程
地震工程学作业课程名称:地震工程学______指导老师:_______翟永梅_________姓名:史先飞________学号: 1232627________一、地震波生成反应谱1 所取的地震波为Elcentro地震波加速度曲线，如图1所示。
图1 Elcentro地震波加速度曲线2 所调用的Matlab程序为：% ***********读入地震记录***********ElCentro;Accelerate= ElCentro(:,1)*;%单位统一为m和sN=length(Accelerate);%N 读入的记录的量time=0::(N-1)*; %单位 s%初始化各储存向量Displace=zeros(1,N); %相对位移Velocity=zeros(1,N); %相对速度AbsAcce=zeros(1,N); %绝对加速度% ***********A,B矩阵***********Damp=; %阻尼比TA=::6; %TA=::6; %结构周期Dt=; %地震记录的步长%记录计算得到的反应，MaxD为某阻尼时最大相对位移，MaxV为某阻尼最大相对速度，MaxA某阻尼时最大绝对加速度，用于画图MaxD=zeros(3,length(TA));MaxV=zeros(3,length(TA));MaxA=zeros(3,length(TA));t=1;for T=::6NatualFrequency=2*pi/T ; %结构自振频率DampFrequency=NatualFrequency*sqrt(1-Damp*Damp); %计算公式化简e_t=exp(-Damp*NatualFrequency*Dt);s=sin(DampFrequency*Dt);c=cos(DampFrequency*Dt);A=zeros(2,2);A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);A(1,2)=e_t*s/DampFrequency;A(2,1)=-NatualFrequency*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)/(NatualFrequency^2*Dt);d_3t=Damp/(NatualFrequency^3*Dt);B=zeros(2,2);B(1,1)=e_t*((d_f+Damp/NatualFrequency)*s/DampFrequency+(2*d_3t+1/NatualFrequency^2)*c)-2*d_3 t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)-1/NatualFrequency^2+2*d_3t;B(2,1)=e_t*((d_f+Damp/NatualFrequency)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/NatualFrequency^2 )*(DampFrequency*s+Damp*NatualFrequency*c))+1/(NatualFrequency^2*Dt);B(2,2)=e_t*(1/(NatualFrequency^2*Dt)*c+s*Damp/(NatualFrequency*DampFrequency*Dt))-1/(NatualF requency^2*Dt);for i=1:(N-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*NatualFrequency*Velocity(i+1)-NatualFrequency^2*Displace(i+1);endMaxD(1,t)=max(abs(Displace));MaxV(1,t)=max(abs(Velocity));if T==MaxA(1,t)=max(abs(Accelerate));elseMaxA(1,t)=max(abs(AbsAcce));endDisplace=zeros(1,N);%初始化各储存向量，避免下次不同周期计算时引用到前一个周期的结果Velocity=zeros(1,N);AbsAcce=zeros(1,N);t=t+1;End% ***********PLOT***********close allfigure %绘制地震记录图plot(time(:),Accelerate(:))title('PEER STRONG MOTION DATABASE RECORD')xlabel('time(s)')ylabel('acceleration(g)')gridfigure %绘制位移反应谱plot(TA,MaxD(1,:),'',TA,MaxD(2,:),'-r',TA,MaxD(3,:),':k')title('Displacement')xlabel('Tn(s)')ylabel('Displacement(m)')legend('ζ=')Gridfigure %绘制速度反应谱plot(TA,MaxV(1,:),'',TA,MaxV(2,:),'-r',TA,MaxV(3,:),':k')title('Velocity')xlabel('Tn(s)')ylabel('velocity(m/s)')legend('ζ=')Gridfigure %绘制绝对加速度反应谱plot(TA,MaxA(1,:),'',TA,MaxA(2,:),'-r',TA,MaxA(3,:),':k')title('Absolute Acceleration')xlabel('Tn(s)')ylabel('absolute acceleration(m/s^2)')legend('ζ=')Grid3 运行的结果得到的反应谱图2 位移反应谱图3 速度反应谱图4 加速度反应谱一、反应谱生成地震波1所取的反应谱为上海市设计反应谱图5 上海市设计反应谱2反应谱取值程序为：%%规范反应谱取值程序参照01年抗震规范function rs_z=r_s_1(pl,zn,ld,cd,fz) %%%pl 圆频率,zn阻尼比,ld烈度,cd场地类型,场地分组fz %%%%烈度选择if ld==6arfmax=;endif ld==7arfmax=;endif ld==8arfmax=;endif ld==9arfmax=;end%%%%场地类别，设计地震分组选择if cd==1if fz==1Tg=;endif fz==2Tg=;endif fz==3Tg=;endendif cd==2if fz==1Tg=;endif fz==2Tg=;endif fz==3Tg=;endendif cd==3if fz==1Tg=;endif fz==2Tg=;endif fz==3Tg=;endendif cd==4if fz==1Tg=;endif fz==2Tg=;endif fz==3Tg=;endend%%%%%%%%%ceita=zn; %%%%%阻尼比lmt1=+/8;if lmt1<0lmt1=0;endlmt2=1+/+*ceita);if lmt2<lmt2=;endsjzs=+/+5*ceita);%%%%%分段位置 T1 T2 T3T1=;T2=Tg;T3=5*Tg;T_jg=2*pi./pl;%%%% 第一段 0～T1if T_jg<=T1arf_jg=*arfmax+(lmt2**arfmax)/*T_jg;end%%%% 第二段 T1～T2if T1<T_jg&T_jg<=T2arf_jg=lmt2*arfmax;end%%%% 第三段 T2~T3if T2<T_jg&T_jg<=T3arf_jg=((Tg/T_jg)^sjzs)*lmt2*arfmax;end%%%% 第四段 T3～if T3<T_jg&T_jg<=arf_jg=(lmt2*^sjzs-lmt1*(T_jg-5*Tg))*arfmax; end%%%% 第五段～if <T_jgarf_jg=(lmt2*^sjzs-lmt1**Tg))*arfmax;end%%%%%%反应谱值拟加速度值rs_z=arf_jg*;end3生成人造地震波主程序：%%%主程序%%%%%%%%确定需要控制的反应谱Sa（T）(T=T1,...,TM)的坐标点数M，反应谱控制容差rc Tyz=[::,::,::];rc=;nTyz=length(Tyz);ceita=;%%%阻尼比：for i=1:nTyzSyz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1); %%%%8度，2类场地，第1地震分组end%%%%%% 变换的频率差：2*pi*(可以保证长周期项5s附近有5项三角级数)；%%%%频率变化范围 N1=30, 30**2*pi ;N2=3000, 5000**2*piplc=2*pi*;pl=30**2*pi:*2*pi:10000**2*pi;npl=length(pl);P=; %%%保证率%%%%%%人造地震动持续时间40s，时间间隔：Td=40;dt=;t=0::40;nt=length(t);%%%%%%% 衰减包络函数t1=8; %%%%上升段t2=8+24; %%%%%平稳段; 下降段则为40－32＝8sc=; %%%%衰减段参数for i=1:ntif t(i)<=t1f(i)=(t(i)/t1)^2;endif t(i)>t1 & t(i)<t2f(i)=1;endif t(i)>=t2f(i)=exp(-c*(t(i)-t2));endend%%%%%%% 反应谱转换功率谱for i=1:nplSw(i)=(2*ceita/(pi*pl(i)))*r_s_1(pl(i),ceita,8,2,1)^2/(-2*log(-1*pi*log(P)/(pl(i)*Td))); Aw(i)=sqrt(4*Sw(i)*plc);end%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfai(i)=rand(1)*2*pi;for j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求,需要时程动力分析%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0));dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=1;while max(rcrsp(:))>rc%%%%%循环体函数blxs=Syz./rspa;for xsxs=1:nplif 2*pi/pl(xsxs)<Tyz(1)blxs1(xsxs)=blxs(1);endfor sxsx=1:nTyz-1if (2*pi/pl(xsxs)>=Tyz(sxsx)) & (2*pi/pl(xsxs)<=Tyz(sxsx+1))blxs1(xsxs)=blxs(sxsx)+(blxs(sxsx+1)-blxs(sxsx))*(2*pi/pl(xsxs)-Tyz(sxsx))/(Tyz(sxsx+1)-Tyz(sxsx));endendif 2*pi/pl(xsxs)>Tyz(nTyz)blxs1(xsxs)=blxs(nTyz);endendAw=Aw.*blxs1;%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfor j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0)); dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=jsnum+1max(rcrsp(:))end%%%%%%% 最终的反应谱与规范谱%%%%%%%%%%%% response spectra of callidar%%%%%%% parameter%% Tjs=::6;%% nTjs=length(Tjs);g=;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0)); dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i)/g;rspa_S(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1)/g;endsubplot(2,1,1);plot(t,at);subplot(2,1,2);plot(Tyz,rspa);hold on;plot(Tyz,rspa_S);4生成的人造地震波如图所示。
人工地震波生成程序简介
姓名：郭 勇 学号：0220020128人工地震波生成程序简介一、 程序设计内容及方法1、程序内容本程序根据特征周期、水平地震波影响系数最大值和地震波幅值等初始条件生成人工地震波，为结构动力分析的时程分析法提供地震波来源。
2、程序设计方法(1) 理论依据本程序采用三角级数法生成人工地震波。
对于给定的功率谱密度函数()x S ω，按照下面的公式可以方便的生成以()x S ω为功率谱密度函数、均值为零的高斯平稳过程()a t 。
1()cos()Nk k k k a t C t ωϕ==+∑ (1)式中：12[4()]()/1()2k x k u l k l C S N k ωωωωωωωω⎫⎪=∆⎪∆=-⎬⎪⎪=+-∆⎭(2)k ϕ为内均匀分(0,2)π布的随机相角；u ω，l ω分别为正域ω内的上、下限值，即认为的有()x S ω效功率在范(,)u l ωω围内，而范围外的()x S ω值可视为零。
为了反映地面运动的非平稳性，采用包络函数乘以平稳()f t 过程()a t ，()()()x t f t a t = (3)(3)式即为人工地震波模型。
()f t 可根据下式确定：2221112()233/01()0c t t t t t t t t t f t e t t t t t T--⎧≤<⎪≤<⎪=⎨≤<⎪⎪<≤⎩ (4)式中：c 为衰减系数，通常取值范围为0.1~1.0，本程序取0.15；1t ，2t 和根据不同3t 实际情况取值，T 为地震波持时，本程序取1t ，2t 分别为4s ，15s ，3t 和均为40T s 。
本程序采用《建筑抗震设计规范》(GB50011-2001)中的反应谱作为目标谱，通过Kau l 提出的平稳过程反应谱与功率谱的近似关系22()[()]/[2ln(ln )]Tx k a k kk dS S p T ξπωωπωω=--(5) 式中：()Ta k S ω为规范反应谱；ξ为阻尼比；d T 为地震动持时；p 为反应不超过反应谱值的概率，本程序取0.85。
SV波斜入射时双线并行地铁隧道横截面地震响应分析
Vol. 41 No. 5Oct. 2021第41卷第5期2021年10月地震工程与工程振动EARTHQUAKE ENGINEERING AND ENGINEERING DYNAMICS文章编号：1000 -1301 (2021)05 -0001 -12 DOI ：10.13197/j. eeev. 2021.05.1. zhousx. 001SV 波斜入射时双线并行地铁隧道横截面地震响应分析周双喜，叶国涛，张季(华东交通大学土木建筑学院，江西南昌330013)摘要:SV 波以大角度斜入射时，场地伴随更大的竖向地震作用，这很可能使得地铁隧道地震响应 特点异于小角度斜入射情况。
基于粘弹性人工边界理论,采用频域刚度矩阵法计算任意角度斜入射SV 波的地震动输入等效节点力，通过ABAQUS 有限元软件建立自由场模型，验证了 0。
、30。
、40。
和 50。
斜入射SV 波地震动输入的准确性;在此基础上建立任意角度斜入射SV 波作用下的地下双线并行圆形地铁隧道地震响应分析数值模型，从场地类别、隧道埋深和双隧道间距等方面分析SV 波入射角度对隧道结构横断面地震响应的影响。
研究表明:入射角度0。
及略大于SV 波临界角的入射角度 是浅埋圆形地铁隧道结构抗震的最不利角度，尤其对于II 类场地，若仅在临界角以内研究衬砌结构 的地震响应，将显著低估衬砌结构的动力响应;此外,隧道埋深对衬砌动弯矩、动轴力峰值具有显著 影响，这种影响与场地类别密切相关;II 类和III 类场地中小隧道间距条件下隧道间的动力相互作用对衬砌结构动力响应具有一定的增大效应，而IV 类场地中可忽略隧道间的动力相互作用对衬砌结 构动力响应所产生的影响。
关键词:SV 波斜入射;频域刚度矩阵法;SV 波临界角；隧道埋深;隧道间距中图分类号:TU315.9文献标识码:ATransverse seismic response analysis of twin parallel metro tunnelssubjected to obliquely incident SV waveZHOU Shuangxi, YE Guolao, ZHANG Ji(School of Civil Engineering and Architecture , East China Jiaotong University , Nanchang 330013 , China )Abstract : When SV wave is incident obliquely at a large angle , a greater vertical seismic action will be exerted onsite , which may make seismic response characteristics of the tunnels different from the case of oblique incidence at small angle ・ In this paper , based on viscoelastic artificial boundary theory , equivalent nodal forces of SV waves in cident obliquely at any angle are calculated by adopting frequency domain stiffness matrix method , and accuracy ofground motion input of 0°, 30°, 40° and 50° oblique incident SV waves is verified through numerical free field models established by ABAQUS finite element software. On above basis , numerical models for seismic response a-nalysis of underground twin parallel circular metro tunnels under obliquely incident SV wave are established , and effects of incident angle of SV wave on transverse seismic response of the tunnels are investigated from aspects of收稿日期：2020-ll -25；修订日期：2021 -01 -05基金项目：国家自然科学基金项目(52008170, 51968022)；江西省自然科学基金项目(20171BAB216049, 20181BAB206007)Supported by ： National Natural Science Foundation of China (52008170, 51968022) ； Natural Science Foundation of Jiangxi Province(20171BAB216049, 20181BAB206007)作者简介:周双喜(1973 -)，男，教授，主要从事混凝土结构耐久性与工程结构抗震方面研究.E-mail ： green. 55@ 163. com通讯作者:张 季(1985 -)，男,副教授，主要从事地震工程方面研究.E-mail ： jizhang@ ecjtu. edu. cn2地震工程与工程振动第41卷soil types,tunnel buried depths and tunnels spacings.It is shown that incident angle of0°and the angle slighdy greater than critical angle are the most unfavorable angles for seismic resistance of shallow buried circular tunnel structure,especially for class II site,if seismic response of the lining is studied with cases that the incident angles are not greater than the critical angle,dynamic response of tunnel lining would be significantly underestimated.Besides,the buried depth has significant influences,closely related to the soil types,on peaks of dynamic bending moments and axial forces；dynamic interaction between tunnels with medium and small tunnel spacings in class II and III sites may amplify dynamic response of the tunnels lining,while the interaction between tunnels in class IV site may exert an ignored effect on the dynamic response of lining.Key words:oblique incidence of SV wave；frequency domain stiffness matrix method；critical angle of SV wave；tunnel depth;tunnel spacing引言近些年来，地下隧道结构发生了诸如衬砌开裂、剥落、错台和坍塌等不同程度的震害现象⑴，隧道结构的抗震安全性已受到国内外学者和工程界的广泛关注和重视S3］。
安评师第二课地震工程
场地地震动参数确定
地震地质灾害评价
地震动峰值加速度复核
不同安评工作之间的关系
主要内容
关键内容
地震动衰减关系的统计回归
地震烈度和强震动资料的收集原则与内容
基岩地震动衰减模型建立
缺乏强震动观测资料地区的地震动衰减关系确定的原则与方法
地震烈度衰减关系选择原则与方法。
基岩地震动衰减关系选取与适用性分析
土动力试验的测试要求和内容
场地地震工程地质条件勘测
地震地质灾害的场地勘查 地基土液化场地勘查内容和要求 崩塌、滑坡、地裂缝和泥石流的勘查内容和要求 地表断层影响分析所需的资料内容
关键内容
1
2
场地地震工程地质条件勘测
勘测范围与钻孔布设
场地地震工程地质条件勘测
勘测范围 I、II级工作，场地范围可取其建设工程所覆盖的范围；小区划工作，场地范围可取区划所覆盖的范围 钻孔布设 Ⅰ级工作 钻孔布置应能控制工程场地的工程地质条件，控制孔应不少于3个 Ⅱ级工作 钻孔布置应能控制工程场地的工程地质条件，控制孔应不少于2个 地震小区划 钻孔布置应能控制土层结构和工程场地内不同工程地质单元，每个工程地质单元内应至少有1个控制孔
现场原位波速测量，试验室波速测量（海洋工程）
波速测试方法
场地地震工程地质条件勘测
测定：G/G0——曲线，——曲线
共振柱试验，10-6-10-4应变范围
动三轴试验，10-4-10-3应变范围
土体动力性能测定要求与方法
场地地震工程地质条件勘测
场地地震工程地质条件勘测
土体动力性能测定结果
地震动的工程特性分析
地震动衰减关系确定
场地地震工程地质条件勘测
区域性地震区划
地震工程学大作业代码(求ELCENTRO波的位移
M=10^5*diag([3.4 3.4 3.2 3.2 3.0 2.8 2.7 2.6]); %输入质量矩阵K=10^8*(diag([4.4 4.2 4 3.8 3.6 3.4 3.2 1.6])+...diag([-2.2 -2.0 -2.0 -1.8 -1.8 -1.6 -1.6],1)+...diag([-2.2 -2.0 -2.0 -1.8 -1.8 -1.6 -1.6],-1)); %刚度矩阵[fai,w]=eig(inv(K)*M); %特征值问题求解omiga=1./diag(w.^0.5); %得到周期for i=1:8 %归一化振型[val,poi]=max(abs(fai(:,i)));fai(:,i)=fai(:,i)/fai(poi,i);endgama=fai'*M*ones(8,1)./diag(fai'*M*fai); %计算振型参与系数kesai=(0.1347+0.006306*omiga.^2)/2./omiga; %计算各阶阻尼file_id='D:\地震动响应图\';mkdir(file_id); %创建文件夹保存生成的所有图A= textread('NS.txt','%n'); %读取地震波数据N=length(A);mint=0.01; %求时间间隔T=(N-1)*mint;B=xlsread('T');allmotion=zeros(N,8); %用来存储所有时刻的位移allacceleration=zeros(N,8); %用来存储所有时刻的加速度值allacceleration(1,:)=A(1); %对0时刻加速度矩阵初始化allvelocity=zeros(N,8); %用于存储所有时刻速度vg(1)=0;xg(1)=0;for m=1:N-1 %计算场地速度、位移用于对加速度进行基线修正vg(m+1)=vg(m)+A(m)*mint+(A(m+1)-A(m))*mint/2;xg(m+1)=xg(m)+vg(m)*mint+A(m)*mint^2/2+mint^2*(A(m+1)-A(m))/6;endc1=28/13/T^2*(2*vg(N)-15/T^5*xg*(3*T*B(:,1).^2-2*B(:,1).^3)*mint); %计算系数c1c0=(vg(N)-c1*T^2/2)/T; %计算系数c0A=max(A)/max(abs(A-(c0+c1*B(:,1))))*(A-(c0+c1*B(:,1))); %修正后的地震波加速度xg=xg'-(0.5*c0*B(:,1).^2+1/6*c1*B(:,1).^3); %修正后的场地位移for i=1:8 %第三层循环用于对不同周期值计算反应S=1+kesai(i)*omiga(i)*mint+((omiga(i)*mint)^2)/6;for j=2:N %对每个周期利用NewMark方法进行数值积分得到对应值（具体请查看word文档）Q=2*omiga(i)*kesai(i)*allacceleration(j-1,i)*mint+...(omiga(i)^2)*(allvelocity(j-1,i)*mint+0.5*allacceleration(j-1,i)*mint^2);deltaa=-(gama(i)*(A(j)-A(j-1))+Q)/S;allacceleration(j,i)=allacceleration(j-1,i)+deltaa;allvelocity(j,i)= allvelocity(j-1,i)+allacceleration(j-1,i)*mint+deltaa*mint/2;allmotion(j,i)= allmotion(j-1,i)+allvelocity(j-1,i)*mint+0.5*allacceleration(j-1,i)*mint^2+...deltaa*(mint^2)/6;end %积分计算反应并存放在前面的空矩阵里endmotion1=allmotion*fai'; %计算位移acceleration=allacceleration*fai'; %计算加速度velocity=allvelocity*fai'; %计算速度for n=1:9figure('color','white') %利用figure函数准确的控制画图窗口if n~=9axis([0,50,-0.03,0.03]) %控制坐标轴范围，便于比较endgrid onhold onbox off%annotation('arrow',[0.132 0.132],[0.8 1]);%annotation('arrow',[0.8 1],[0.108 0.108]); %产生坐标轴箭头（可能不够美观）if n==1 %计算相关数据和得到图名name='底层层间位移时程';motion(:,n)=motion1(:,n);elseif n==9name='顶层位移时程';motion(:,9)=motion1(:,8)+xg;elsemotion(:,n)=motion1(:,n)-motion1(:,n-1);name=strcat('第',num2str(n-1),'-',num2str(n),'层间位移时程');endplot(B(:,1),motion(:,n),'linewidth',2) %画相对位移时程xlabel('时间（s）','FontName','宋体','FontSize',16); %x轴标注ylabel('位移(m)','FontName','宋体','FontSize',16); %y轴标注title(name,'FontName','宋体','FontSize',20) %标注图名set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看shgF=getframe(gcf);imwrite(F.cdata,[file_id,strcat(name,'.png')]); %存储得到的时程图endGK=repmat(-[-2 -2 -1.8 -1.8 -1.8 -1.8 -1.6 -1.6]*10^8,N,1);%输入刚度矩阵motion(:,8)=motion1(:,8)-motion1(:,7); %计算顶层的层间相对位移Fq=motion(:,1:8).*GK; %计算剪力figure('color','white') %生成图形窗口，背景白色plot(B(:,1),Fq(:,1)) %画底层剪力时程图xlabel('时间（s）','FontName','宋体','FontSize',16); %x轴标注ylabel('剪力（N）','FontName','宋体','FontSize',16); %y轴标注title('底层剪力时程','FontName','宋体','FontSize',20); %标注图名set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看F=getframe(gcf); %存储生成的图形imwrite(F.cdata,[file_id,strcat('底层剪力时程图','.png')]);maxFq=zeros(1,16);for f=1:8 %计算所需剪力和位移包络图数据，并得到便于画图的矩阵maxmotion(2*f-1)=max(abs(motion(:,f)));maxmotion(2*f)=maxmotion(:,2*f-1);maxFq(2*f-1)=max(abs(Fq(:,f)));maxFq(2*f)=max(abs(Fq(:,f)));endfigure('color','white')plot(maxmotion,[0 0.5 0.5 1 1 1.5 1.5 2 2 2.5 2.5 3 3 3.5 3.5 4],'k-');%生成层间位移包络图title('层间位移绝对值最大值包络图','FontName','宋体','FontSize',20); %标注图名text(maxmotion(1:2:16),0.25:0.5:4,num2str(maxmotion(1:2:16)')) %标注位移数据set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看F=getframe(gcf);imwrite(F.cdata,[file_id,strcat('层间位移绝对值最大值包络图','.png')]); %存储图形figure('color','white')plot(maxFq/1000,[0 0.5 0.5 1 1 1.5 1.5 2 2 2.5 2.5 3 3 3.5 3.5 4],'k-');title('层间剪力绝对值最大值包络图','FontName','宋体','FontSize',20); %标注图名text(maxFq(1:2:16)/1000,0.25:0.5:4,num2str(0.001*maxFq(1,1:2:16)'),'FontSize',10);set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看F=getframe(gcf);imwrite(F.cdata,[file_id,strcat('层间剪力绝对值最大值包络图','.png')]);。
地震工程学-反应谱和地震时程波的相互转化matlab编程
地震工程学作业课程名称:地震工程学______ 指导老师:_______翟永梅_________姓名:史先飞________学号:1232627________一、地震波生成反应谱1 所取的地震波为Elcentro地震波加速度曲线，如图1所示。
图1 Elcentro地震波加速度曲线2 所调用的Matlab程序为：% ***********读入地震记录***********ElCentro;Accelerate= ElCentro(:,1)*9.8067;%单位统一为m和sN=length(Accelerate);%N 读入的记录的量time=0:0.005:(N-1)*0.005; %单位s%初始化各储存向量Displace=zeros(1,N); %相对位移Velocity=zeros(1,N); %相对速度AbsAcce=zeros(1,N); %绝对加速度% ***********A,B矩阵***********Damp=0.02; %阻尼比0.02TA=0.0:0.05:6; %TA=0.000001:0.02:6; %结构周期Dt=0.005; %地震记录的步长%记录计算得到的反应，MaxD为某阻尼时最大相对位移，MaxV为某阻尼最大相对速度，MaxA某阻尼时最大绝对加速度，用于画图MaxD=zeros(3,length(TA));MaxV=zeros(3,length(TA));MaxA=zeros(3,length(TA));t=1;for T=0.0:0.05:6NatualFrequency=2*pi/T ; %结构自振频率DampFrequency=NatualFrequency*sqrt(1-Damp*Damp); %计算公式化简e_t=exp(-Damp*NatualFrequency*Dt);s=sin(DampFrequency*Dt);c=cos(DampFrequency*Dt);A=zeros(2,2);A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);A(1,2)=e_t*s/DampFrequency;A(2,1)=-NatualFrequency*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)/(NatualFrequency^2*Dt);d_3t=Damp/(NatualFrequency^3*Dt);B=zeros(2,2);B(1,1)=e_t*((d_f+Damp/NatualFrequency)*s/DampFrequency+(2*d_3t+1/NatualFrequency^2)*c) -2*d_3t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)-1/NatualFrequency^2+2*d_3t;B(2,1)=e_t*((d_f+Damp/NatualFrequency)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/NatualFreque ncy^2)*(DampFrequency*s+Damp*NatualFrequency*c))+1/(NatualFrequency^2*Dt);B(2,2)=e_t*(1/(NatualFrequency^2*Dt)*c+s*Damp/(NatualFrequency*DampFrequency*Dt))-1/(N atualFrequency^2*Dt);for i=1:(N-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*NatualFrequency*Velocity(i+1)-NatualFrequency^2*Displace(i+1); endMaxD(1,t)=max(abs(Displace));MaxV(1,t)=max(abs(Velocity));if T==0.0MaxA(1,t)=max(abs(Accelerate));elseMaxA(1,t)=max(abs(AbsAcce));endDisplace=zeros(1,N);%初始化各储存向量，避免下次不同周期计算时引用到前一个周期的结果Velocity=zeros(1,N);AbsAcce=zeros(1,N);t=t+1;End% ***********PLOT***********close allfigure %绘制地震记录图plot(time(:),Accelerate(:))title('PEER STRONG MOTION DATABASE RECORD')xlabel('time(s)')ylabel('acceleration(g)')gridfigure %绘制位移反应谱plot(TA,MaxD(1,:),'-.b',TA,MaxD(2,:),'-r',TA,MaxD(3,:),':k')title('Displacement')xlabel('Tn(s)')ylabel('Displacement(m)')legend('ζ=0.02')Gridfigure %绘制速度反应谱plot(TA,MaxV(1,:),'-.b',TA,MaxV(2,:),'-r',TA,MaxV(3,:),':k')title('Velocity')xlabel('Tn(s)')ylabel('velocity(m/s)')legend('ζ=0.02')Gridfigure %绘制绝对加速度反应谱plot(TA,MaxA(1,:),'-.b',TA,MaxA(2,:),'-r',TA,MaxA(3,:),':k') title('Absolute Acceleration')xlabel('Tn(s)')ylabel('absolute acceleration(m/s^2)')legend('ζ=0.02')Grid3 运行的结果得到的反应谱图2 位移反应谱图3 速度反应谱图4 加速度反应谱一、反应谱生成地震波1所取的反应谱为上海市设计反应谱图5 上海市设计反应谱2反应谱取值程序为：%%规范反应谱取值程序参照01年抗震规范function rs_z=r_s_1(pl,zn,ld,cd,fz) %%%pl 圆频率,zn阻尼比,ld烈度,cd场地类型,场地分组fz %%%%烈度选择if ld==6arfmax=0.11;endif ld==7endif ld==8arfmax=0.45;endif ld==9arfmax=0.90;end%%%%场地类别，设计地震分组选择if cd==1if fz==1Tg=0.25;endif fz==2Tg=0.30;endif fz==3Tg=0.35;endendif cd==2if fz==1Tg=0.35;endif fz==2Tg=0.40;endif fz==3Tg=0.45;endendif cd==3if fz==1Tg=0.45;endif fz==2Tg=0.55;endif fz==3Tg=0.65;endendif cd==4if fz==1endif fz==2Tg=0.75;endif fz==3Tg=0.90;endend%%%%%%%%%ceita=zn; %%%%%阻尼比lmt1=0.02+(0.05-ceita)/8;if lmt1<0lmt1=0;endlmt2=1+(0.05-ceita)/(0.06+1.7*ceita);if lmt2<0.55lmt2=0.55;endsjzs=0.9+(0.05-ceita)/(0.5+5*ceita);%%%%%分段位置T1 T2 T3T1=0.1;T2=Tg;T3=5*Tg;T_jg=2*pi./pl;%%%% 第一段0～T1if T_jg<=T1arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg; end%%%% 第二段T1～T2if T1<T_jg&T_jg<=T2arf_jg=lmt2*arfmax;end%%%% 第三段T2~T3if T2<T_jg&T_jg<=T3arf_jg=((Tg/T_jg)^sjzs)*lmt2*arfmax;end%%%% 第四段T3～6.0if T3<T_jg&T_jg<=6.0arf_jg=(lmt2*0.2^sjzs-lmt1*(T_jg-5*Tg))*arfmax;end%%%% 第五段6.0～if 6.0<T_jgarf_jg=(lmt2*0.2^sjzs-lmt1*(6.0-5*Tg))*arfmax;end%%%%%%反应谱值拟加速度值rs_z=arf_jg*9.8;end3生成人造地震波主程序：%%%主程序%%%%%%%%确定需要控制的反应谱Sa（T）(T=T1,...,TM)的坐标点数M，反应谱控制容差rc Tyz=[0.04:0.016:0.1,0.15:0.05:3.0,3.2:0.05:5.0];rc=0.06;nTyz=length(Tyz);ceita=0.035;%%%阻尼比：0.035for i=1:nTyzSyz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1); %%%%8度，2类场地，第1地震分组end%%%%%% 变换的频率差：2*pi*0.005(可以保证长周期项5s附近有5项三角级数)；%%%%频率变化范围N1=30, 30*0.005*2*pi ;N2=3000, 5000*0.005*2*piplc=2*pi*0.005;pl=30*0.005*2*pi:0.005*2*pi:10000*0.005*2*pi;npl=length(pl);P=0.9; %%%保证率%%%%%%人造地震动持续时间40s，时间间隔：0.02sTd=40;dt=0.02;t=0:0.02:40;nt=length(t);%%%%%%% 衰减包络函数t1=8; %%%%上升段t2=8+24; %%%%%平稳段; 下降段则为40－32＝8sc=0.6; %%%%衰减段参数for i=1:ntif t(i)<=t1f(i)=(t(i)/t1)^2;endif t(i)>t1 & t(i)<t2f(i)=1;endif t(i)>=t2f(i)=exp(-c*(t(i)-t2));endend%%%%%%% 反应谱转换功率谱for i=1:nplSw(i)=(2*ceita/(pi*pl(i)))*r_s_1(pl(i),ceita,8,2,1)^2/(-2*log(-1*pi*log(P)/(pl(i)*Td)));Aw(i)=sqrt(4*Sw(i)*plc);end%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfai(i)=rand(1)*2*pi;for j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求,需要时程动力分析%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0));dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=1;while max(rcrsp(:))>rc%%%%%循环体函数blxs=Syz./rspa;for xsxs=1:nplif 2*pi/pl(xsxs)<Tyz(1)blxs1(xsxs)=blxs(1);endfor sxsx=1:nTyz-1if (2*pi/pl(xsxs)>=Tyz(sxsx)) & (2*pi/pl(xsxs)<=Tyz(sxsx+1))blxs1(xsxs)=blxs(sxsx)+(blxs(sxsx+1)-blxs(sxsx))*(2*pi/pl(xsxs)-Tyz(sxsx))/(Tyz(sxsx+1)-Tyz(sxs x));endendif 2*pi/pl(xsxs)>Tyz(nTyz)blxs1(xsxs)=blxs(nTyz);endendAw=Aw.*blxs1;%%%%%%%%%%%%%% 合成地震动at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfor j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%% 计算反应谱验证是否满足rc在5%的要求%%%%%%%%%%%% response spectra of callidar%%%%%%% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0));dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%%%%%%% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=jsnum+1max(rcrsp(:))end%%%%%%% 最终的反应谱与规范谱%%%%%%%%%%%% response spectra of callidar%%%%%%% parameter%% Tjs=0.05:0.01:6;%% nTjs=length(Tjs);g=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%% loadag=at; %%%%%%%修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0));dx=kk\pp;dv=dt\(3*dx)-3*v0-2\(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i)/g;rspa_S(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1)/g;endsubplot(2,1,1);plot(t,at);subplot(2,1,2);plot(Tyz,rspa);hold on;plot(Tyz,rspa_S);4生成的人造地震波如图所示。
时程分析时地震波的选取及地震波的反应谱化
时程分析时地震波的选取及地震波的反应谱化摘要：目前我国规范要求结构计算中地震作用的计算方法一般为振型分解反应谱法。
时程分析法作为补充计算方法，在不规则、重要或较高建筑中采用。
进行时程分析时，首先面临正确选择输入的地震加速度时程曲线的问题。
时程曲线的选择是否满足规范的要求，则需要首先将时程曲线进行单自由度反应计算，得到其反应谱曲线，并按规范要求和规范反应谱进行对比和取舍。
本文通过介绍常用的数值计算方法及计算步骤，实现将地震加速度时程曲线计算转化成反应谱曲线，从而为特定工程在时程分析时地震波的选取提供帮助。
关键词：时程分析，地震波，反应谱，动力计算1 地震反应分析方法的发展过程结构的地震反应取决于地震动和结构特性。
因此，地震反应分析的水平也是随着人们对这两个方面认识的深入而提高的。
结构地震反应分析的发展可以分为静力法、反应谱法、动力分析法这三个阶段。
在动力分析法阶段中又可分为弹性和非弹性（或非线性）两个阶段。
[1]目前，在我国和其他许多国家的抗震设计规范中，广泛采用反应谱法确定地震作用，其中以加速度反应谱应用得最多。
反应谱是指：单自由度弹性体系在给定的地震作用下，某个最大反应量（如加速度、速度、位移等）与体系自振周期的关系曲线。
反应谱理论是指：结构物可以简化为多自由度体系，多自由度体系的地震反应可以按振型分解为多个单自由度体系反应的组合，每个单自由度体系的最大反应可以从反应谱求得。
其优点是物理概念清晰，计算方法较为简单，参数易于确定。
反应谱理论包括如下三个基本假定：1、结构物的地震反应是弹性的，可以采用叠加原理来进行振型组合；2、现有反应谱假定结构的所有支座处地震动完全相同；3、结构物最不利的地震反应为其最大地震反应，而与其他动力反应参数，如最大值附近的次数、概率、持时等无关。
[1]时程分析法是对结构物的运动微分方程直接进行逐步积分求解的一种动力分析方法。
由于此法是对运动方程直接求解，又称直接动力分析法。
矿产
矿产资源开发利用方案编写内容要求及审查大纲
矿产资源开发利用方案编写内容要求及《矿产资源开发利用方案》审查大纲一、概述
㈠矿区位置、隶属关系和企业性质。
如为改扩建矿山, 应说明矿山现状、
特点及存在的主要问题。
㈡编制依据
(1简述项目前期工作进展情况及与有关方面对项目的意向性协议情况。
(2 列出开发利用方案编制所依据的主要基础性资料的名称。
如经储量管理部门认定的矿区地质勘探报告、选矿试验报告、加工利用试验报告、工程地质初评资料、矿区水文资料和供水资料等。
对改、扩建矿山应有生产实际资料, 如矿山总平面现状图、矿床开拓系统图、采场现状图和主要采选设备清单等。
二、矿产品需求现状和预测
㈠该矿产在国内需求情况和市场供应情况
1、矿产品现状及加工利用趋向。
2、国内近、远期的需求量及主要销向预测。
㈡产品价格分析
1、国内矿产品价格现状。
2、矿产品价格稳定性及变化趋势。
三、矿产资源概况
㈠矿区总体概况
1、矿区总体规划情况。
2、矿区矿产资源概况。
3、该设计与矿区总体开发的关系。
㈡该设计项目的资源概况
1、矿床地质及构造特征。
2、矿床开采技术条件及水文地质条件。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
矿产资源开发利用方案编写内容要求及审查大纲
矿产资源开发利用方案编写内容要求及《矿产资源开发利用方案》审查大纲一、概述
㈠矿区位置、隶属关系和企业性质。
如为改扩建矿山, 应说明矿山现状、
特点及存在的主要问题。
㈡编制依据
(1简述项目前期工作进展情况及与有关方面对项目的意向性协议情况。
(2 列出开发利用方案编制所依据的主要基础性资料的名称。
如经储量管理部门认定的矿区地质勘探报告、选矿试验报告、加工利用试验报告、工程地质初评资料、矿区水文资料和供水资料等。
对改、扩建矿山应有生产实际资料, 如矿山总平面现状图、矿床开拓系统图、采场现状图和主要采选设备清单等。
二、矿产品需求现状和预测
㈠该矿产在国内需求情况和市场供应情况
1、矿产品现状及加工利用趋向。
2、国内近、远期的需求量及主要销向预测。
㈡产品价格分析
1、国内矿产品价格现状。
2、矿产品价格稳定性及变化趋势。
三、矿产资源概况
㈠矿区总体概况
1、矿区总体规划情况。
2、矿区矿产资源概况。
3、该设计与矿区总体开发的关系。
㈡该设计项目的资源概况
1、矿床地质及构造特征。
2、矿床开采技术条件及水文地质条件。