地震动反映谱的数值计算（含fortran程序）.doc

地震动反映谱的数值算冯宏光（云南大学资源环境与地球科学学院，云南 昆明 650091）摘要：本文从单质点系入手，通过解方程而得到直接积分法中的线性加速度法的递推公式，通过递推公式以及线性加速度发的主要思想，利用fortran编写程序得到地震动加速度反应谱的计算程序进而计算出地震动加速度度反应谱。关键词：单质点系、线性加速度法、地震动反应谱Compute the Date of Earthquake Response SpectrumFeng hongguang（Yunnan University, School of Resources Environment and Earth Scitnces,Kunming Yunnan China 650091）Abstract: Thesis is base on the single quality point department, and square distance in solution got the formula of linear acceleration method, and make use of the main thought of linear acceleration method, using the procedure of fortran got the earthquake response spectrum.Key words：single quality point department、linear acceleration method、earthquake response spectrum 序言地震是一项难以预报的自然灾害，随着人类文明的发展地震所带来的危害被更多人重视，其一次性所带来的灾害往往是难以用简单的经济数据来描述的，而地震中危害最大的就是房屋的倒塌所带来的损失。进入二十一世纪以来地震安全性评估就提出了小震不坏，中震可修，大震不倒的宗旨。因此伴随着计算机技术的急速发展利用计算机编程计算地震动反应谱成为了普遍方法。通过研究地震动反应谱我们可以的道地震动与地面建筑响应的关系，通过了解这种关系就可以指导我们修建地面建筑，因此反应谱的研究无疑是为我们地面建筑的修建指明了一条光明大道。因此我们可以看到研究地震动反应谱的必要性。第一章：单质点系的震动我们知道地震波是弹性波，弹性波最简单的模型便是单质点系的振动，同时反应谱本身是以单质点系振动理论为基础建立的，单质点系的振动对于研究反应谱有很重要的意义。我们采用Kelvin Voit模型来研究单质点振动，此模型为一弹簧与一个活塞并联。我们设质点的质量为m，弹簧的弹性系数为k，阻尼器的阻尼系数为c，并且设该模型只可以在一个方向上是自由的。如上图所示x+y为质点的绝对位移，x为质点的相对位移。那么根据牛顿第三定律我们可以得到： （1.1）1.1考虑了模型中同时存在阻力与粘滞力。当仅仅是考虑简单的一维无阻尼振动时我们可以得到：，我们可以解得： （1.2）根据初始条件：将A、B带入1.2式我们可以得到： （1.3）我们便可以得到一维无阻尼振动的特解当我们引入阻尼后我们便可以得到： （1.4）通过解1.4式我们可以解得： （1.5）1.5式是对于有阻尼的情况的通解，令；：我们便可以得到： （1.6）引入初始条件：我们可以得到： （1.7）第二章：地震动反应第一章的结论只考虑了单质点震动的问题，但是在实际情况中以单摆式地震仪为例，第一章的结果并没有考虑仪器的响应问题，也就是说仪器对于单质点产生的地震动会做何反应，大而言之我们可以认为地面建筑对于地震动的反应是如何的。在解决了单质点有阻尼的振动后我们便要考虑到波动响应的问题。对第一章中的（1.1）便是考虑了场地反应的问题。已知：，其中是由地震动引起的惯性力。该式的物理意义可以理解为，假定地面在空间中固定不动，而惯性力代替了地面的运动。惯性力与系统的弹性恢复力、阻尼力构成平衡力系，从而把运动状态变成了静止状态。在第一章中已知，那么1.1式可以成为： （2.1）惯性力是作用于质点上的随时间而变化的力，因此可令：于是（1.1）式为： （2.2）2.2的解事任意力F（t）作用于质点上，在任意时刻质点的相对位移x（t）。其解为：现在用替换F（t），得 （2.3）它是单质点阻尼系对地震动的位移反应。因为学x（t）是质点对地面的相对位移，所以一我们可以认为它是相对位移的反应。那么如果我们对2.3式对t求导我们是否可以得到相对速度反应。同理我们便可以得到地震动加速度反应。由于位移与速度都已经知道，那么我们便可以求出地震动的加速度反应：（2.4）加速度反应单质点阻尼系在地震加速度作用下的位移、速度、加速度反应可按2.4式求得。也可以利用脉冲反应函数把位移写成的形式方便计算： （2.5）其中， （2.6）它们分别是给定输入加速度时，单质点系关于位移、速度、加速度的脉冲反应函数。它们都是因果性的时间函数。具体的计算方法我们将在讨论完反应谱后集中探讨。第三章：地震动反应谱地震动反应谱可以定义为具有同一阻尼比的一系列单自由度体系（其自振周期为，i=1、2、3n）的最大反应绝对值S(,)与周期的关系，即S(,)它称为阻尼比为的反应谱。对于绝对加速度而言，它是绝对加速度反应谱，同理也存在速度反应谱，位移反应谱。三者存在下述近似关系： （2.7）地震动反应谱的具体计算公式如下： （2.8） （2.9）（2.10）虽然反应谱是通过单自由度体系的反应来定义的，正如同我们在前面所说到的那样。但是反应谱反应的是地震动过程a（t）的特性，单自由度体系仅仅具有滤波器的作用。我们知道地震波在高频（0.10.5s）部分主要反应加速度，中频（0.52s）反应速度，低频（大于2s）主要反应位移。那么我们通过单质点自由度体系根据我们所知道的共振原理通过控制摆长来控制自振周期达到过滤出加速度、位移、速度。首先反应谱的整个形状随a（t）而变，其次近震小震、坚硬场地上a（t）的反应谱峰值在高频部分。远震、大震、软厚场地上a（t）的反应谱峰值在低频部分。最后反应谱的形状并不随结构而变只是对于不同的阻尼比与不同自振周期T的结构，其反映了反应谱的该坐标值而已。因此来说：不同结构的反应取反应谱的不同坐标值，而不同地震动a（t）则有不同的反应；或者说反应谱的形状决定于地震动a（t）。结构特性决定其反应在反应谱中取特点坐标值。根据反应谱可以明了仅仅从记录上无法了解的地震波的种种特性，特别是地震波对结构的影响。不同的地震动有不同的选择性破坏，即使小震、近震在坚硬场地上的地震动容易使刚性结构产生震害，而大震远震软厚场地上的地震容易使高柔结构产生震害。这一规律只能从地震动的谱的特性去理解。第四章：计算方法的给出地震动反应谱的计算有多种方法，我们选择直接积分法作为计算方法，所谓直接积分法就是直接积分该方程式求反应的数值解得方法，这就是直接积分法。我们知道加速度地震计所记录到得数据不是一个连续的数据，而是具有一定时间间隔的离散值（这里我们设地震加速度为时间间隔的离散值）。因此线性加速度法的特点就是在加速度和之间按直线内插取值。那么我们之前所到处的单质点系的振动的问题所得到得方程可以写为： 通过求解该方程后我们可以得到如下递推公式：求得、后我们便可以求得绝对加速度为：这样我们便可取t=0时反应的初值为：得知初值后我们便可以计算出全部时刻的绝对加速度、相对速度反应和相对位移反应。根据这个公式我们便可以得到如下的计算方法：之前我们所给出的递推公式主要计算出的是绝对加速度的时间历程（如下图所给出的不同的固有周期下的阻尼为h=0.05的绝对加速度的时间历程）。通过不断比对每一个绝对加速度的大小就可以得到一个固有周期值所对应的绝对绝对加速度的最大值。这样我们便可以得到加速度反应谱。第五章：结论通过计算地震动反映谱我们可以得出，地震动当地面建筑的固有周期达到T=1.4s时地震加速度反应最为剧烈，地震动加速度达到1521.711GAL。同时我们可以通过该计算方法通过输入不同的地震记录数据便可以的到不同地震波的地震动反映谱，进而将这个计算方法运用到实际的工作中，计算实的地面建筑对于不同地震的地震动向应。致谢：论文结尾 我想特别感谢我的导师胡家富老师，该论文从理论研究到结论计算都受到了胡家富老师的鼎力帮助，特别是对于结果的判断上给予了极大的帮助，同时我还要感谢教授我知识的各位老师，没有老师们的细心教导我就无法有一定的理论基础去学习和研究新的知识内容。特此在论文结尾我特此感谢老师们的教导之恩。参考文献场地条件又J1设计反应谱最大值晌影响李广军 ，赵 艳 ，王文仲 ，张同伟(1佳木斯大学建筑工程学院，黑龙江佳木斯154007；2北京工业大学建筑工程学院，北京100022)非弹性反应谱计算的一个改进算法朱镜清 ，朱晓力 ，王东升(1中国地震局工程力学研究所，黑龙江哈尔滨150080；2，大连理工大学土木水利学院，辽宁大连116024)地震安全性评价技术教程胡聿贤 （地震出版社出版）Matlab程序计算反应谱的误差分析徐源 王毅(1_解放军理工大学工程兵工程学院，南京21000O；2东南大学信息科学与工程学院，南京210000)附录一：地震动反应谱的计算程序c 该程序中各个变量所代表的意义c t（96）是有外部输入的周期数据值c ddy（1024）是由外部输入的加速度地震记录器记录到的地震波数据c x（1024）是通过ddy所计算得到的地震动位移反映c c 定义段 dimension t(96),ddy(1024),x(1024),dx(1024) dimension accmax(96),ddxy(1024) real dt,w,wd,h,wy real Aa,Ab,Ac,Ad,Ba,Bb,Bc,Bdc 输入读入段c 地震波数据读入段 open(1,file=wave.dat,status=old) read(1,*)wy,wy,(ddy(i),i=1,1024)c 周期数据读入段 open(2,file=period.dat,status=old) read(2,*)wy,(t(j),j=1,96)c 输入输出测试段 write(*,*)(ddy(k),k=1,1024) write(*,*)(t(l),l=1,96)c 通过屏幕输入阻尼 write(*,*) please input the data of damp read(*,*)hc 实际计算段c 前期数据准备段 dt=0.02 do 60 m=1,96 w=6.283185/t(m) wd=w*SQRT(1-h*2) aa=exp(-h*w*dt)*(cos(wd)*dt+(h*w/wd)*sin(wd)*dt) ab=exp(-h*w*dt)*(1/wd)*sin(wd)*dt ac=-exp(-h*w*dt)*(w*2/wd)*sin(wd)*dt ad=exp(-h*w*dt)*(cos(wd)*dt-(h*w/wd)*sin(wd)*dt) ba=exp(-h*w*dt)*(1/w*2+2.*h/(w*3*dt)*dt*cos(wd)+ # (h/(w*wd)-(1.-2*h*2)/(w*2)*wd*dt)*dt*sin(wd)- # 2.*h/(w*3)*dt) bb=exp(-h*w*dt)*(-1*(2*h*dt*cos(wd)/(dt*w*3)+ # (1-2.*h*2)*(dt*sin(wd)/(wd*dt*w*2)*dt*sin(wd)- # 1/w*2+2.*h/(dt*w*3) bc=exp(-h*w*dt)*(-dt*cos(wd)/(dt*w*2)-(h/(w*wd*dt)+ # 1/wd)*dt*sin(wd)+1/(dt*w*2) bd=exp(-h*w*dt)*(cos(wd)*dt/(dt*w*2)+(h*dt*sin(wd)/ # (w*wd*dt)-1./(dt*w*2) write(*,*)Aa,Ab,Ac,Ad,Ba,Bb,Bc,Bdc read(*,*)wy x(1)=0 dx(1)=-ddy(1)*dt ddxy=2.*h*w*ddy(1)*dt time=0c 实际循环计算段 do 50 n=2,1024 x(n)=aa*x(n-1)+ab*dx(n-1)+ba*ddy(n-1)+bb*ddy(n) dx(n)=ac*x(n-1)+ad*dx(n-1)+bc*ddy(n-1)+bd*ddy(n) ddxy(n)=-1*(2*h*w*dx(n)+(w*2)*x(n)c if(m=20) thenc open(10,file=ddxy.dat,access=sequential)c time=time+0.02c write(10,*)time,ddxy(n)c elsec go to 80c end ifc 数据比较段 80 if(ac