Elcentro地震波E-W方向
FLAC3D在地震边坡稳定性分析中的应用
第２９卷第５期江西理工大学学报ｖ。
１．２９，Ｎ。
．５２００８年１０月ＪＯＵＲＮＡＬＯＦＪＩＡＮＧＸＩＵＮＩＶＥＲＳＩＴＹＯＦＳＣＩＥＮＣＥＡＮＤＴＥＣＨＮＯＬＯＧＹＯｃｔ．２００８文章编号：１００７—１２２９（２００８）０５－００２３—０４ＦＬＡＣ３Ｄ在地震边坡稳定性分析中的应用张友锋，袁海平（江西理工大学应用科学学院，江西赣州３４１０００）摘要：介绍了ｎＡＣ如基本原理和动力分析的特点；通过对某土质边坡的地震荷载动力分析，讨论了利用ｎＡＣ∞进行边坡地震动力分析时如何设置边界条件、地质体阻尼的选取、地震波的时程转化和输入．在此基础上，依据边坡地震荷栽作用下的边坡位移变形和塑性区域的贯通情况对边坡稳定性进行了分析．关键词：ＦＬＡＣＷ；地震；边坡；稳定性分析；数值模拟中图分类号：ＴＵ４５２文献标识码：ＡＴｈｅＡｐｐｌｉｃａｔｉｏｎｏｆＦＬＡＣ如ｔｏｔｈｅＡｎａｌｙｓｉｓｏｆＳｌｏｐｅＳｔａｂｉｌｉｔｙｉｎＥａｒｔｈｑｕａｋｅＺＨＡＮＧＹｏｕ－ｆｅｎｇ，ＹＵＡＮＨａｉ－ｐｉｎｇ（ＦａｃｕｌｔｙｏｆＡｐｐｌｉｅｄＳｃｉｅｎｃｅ，Ｊｉａｎ弘ｉＵｎｉｖｅｒｓｉｔｙｄＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙ，Ｇａｎｚｈｏｕ３４１０００，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：ＴｈｅｓｔａｂｉｌｉｔｙａｎａｌｙｓｉｓｏｆｓｌｏｐｅｗｉｔｈＦＬＡＣ３０ｉｎｅａｒｔｈｑｕａｋｅｉｓｉｎｔｒｏｄｕｃｅｄｉｎｔｈｉｓｐａｐｅｒ．ＴｈｅｉｓｓｕｅｓｏｆｂａｓｉｃｔｈｅｒｏｙｏｆｄｙｎａｍｉｃａｎａｌｙｓｉｓｗｉｔｈＦＬＡＣ３Ｄ，ｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎｓ，ｄｙｎａｍｉｃｗａｖｅｔｒａｎｓｆｒｏｍａｎｄｌｏａｄ，ａｎｄｄａｍｐｉｎｇｆｏｒｇｅｏｌｏｇｉｃａｌｂｏｄｙａｒｅｄｉｓｃｕｓｓｅｄ．Ｔｈｅｓｔａｂｉｌｉｔｙｏｆｓｌｏｐｅｉｎｅａｒｔｈｑｕａｋｅｉｓｅｓｔｉｍａｔｅｄｂａｓｅｄｏｎｔｈｅｄｉｓｐｌａｃｅｍｅｎｔａｎｄｐｌａｓ－－ｔｉｃｄｉｓｔｏｒｔｉｏｎａｒｅａｔｒａｎｓｆｉｘｉｏｎｗｉｔｈＦＬＡＧ３Ｄ．Ｋｅｙｗｏｒｄｓ：ＦＬＡＣ如；ｅａｒｔｈｑｕａｋｅ；ａｎａｌｙｓｉｓｏｆｓｌｏｐｅｓｔａｂｉｌｉｔｙ；ｎｕｍｅｒｉｃａｌｖａｌｕｅｓｉｍｕｌａｔｉｏｎ地震荷载是触发边坡失稳的重要原因之一，是国内外工程界和学术界一直普遍关注的课题．而稳定性问题则是研究的核心，研究地震边坡稳定性的方法归纳起来主要有以下几种：拟静力法、滑块分析法、实验法、地震边坡的概率分析方法、数值方法【１】．目前，对边坡地震稳定性分析常采用的数值方法有：有限元法、离散元法和快速拉格朗日法．近年来，拉格朗日法由于能解决大变形问题和动力分析而倍受青睐．本文利用ＦＬＡＣ３Ｄ对地震边坡的稳定性进行数值模拟分析．１ＦＬＡＣ∞基本原理和特点ＦＬＡＣ３Ｄ软件［２１是由美国Ｉｔａｓｃａ公司开发的工程计算软件，其基本原理是拉格朗日差分法．拉格朗日法源于流体力学，它通过研究流体各个质点的运动参数（位置坐标、速度和加速度等）随时间变化的规律，综合所有流体质点运动参数的变化，得到整个流体的运动规律．把拉格朗日法移植到固体力学中，把所研究的区域划分成空间网格，其结点就相当于流体的质点，设施柏、口；及ｄｖ／ｄｒ（ｉ＝１，２，３）分别为结点的空问位置、位移、速度及加速度向量元，略为应力张量，则应变速率张量、旋转速率张量可表示为：祭（矽一珞），２，∞Ｆ（ｙ广珞），２当施加载荷时，单元结点的运动方程可表示为：ｏ＇＃ｆｈｏｂ；＝ｐｄｖ／ｄｔ，其中Ｐ为介质密度，ｂ为体力密度，本构方程可表示为：眵月＝凰（听，白，ｋ）（ｉＪ，ｋ＝ｌ，２，３），其中同是应力速率张量，嘲是给定函数表达式，ｋ是考虑加载历史的参数．联立上述方程即可求得应力速度张量、应变速度张量和速度矢量．收稿日期：２００８－０７—１４作者简介：张友锋（１９７９一），男，助教．２４江西理工大学学报２００８牟１０月ＦＬＡＣ３０方法在求解时使用以下３种计算方法：①离散模型方法：连续介质被离散为若干互相连接的六面体单元，作用力均被集中在节点上；②有限差分方法：变量关于空间和时问的一阶导数均用有限差分来近似，运动方程和动力方程均采用显式方法求解；③动态松弛方法：应用质点运动方程求解，通过阻尼使系统运动衰减至平衡状态．ＦＬＡＣ３Ｄ有如下特点：①采用混合离散化方法模拟塑性破裂与塑性流动，比采用有限元法更合理；②采用动态运动方程求解更适合解决物理上的不稳定过程问题；③采用显式解法，不需要存储刚度矩阵，与普通隐式解法相比，大大节约了内存和计算时间；④其运动总方程的显式时间逼近解法对于岩土体的渐进破坏与失稳，以及大变形分析较为适用．由于拉格郎日法基于动力学方程，采用了动态求解方法，因此能够更好模拟动态问题．ＦＬＡＣ３Ｄ动力分析原理中，考虑到结构材料的力学性质和大变形影响，采用非线性振动分析和等价线性振动分析等两种方法．等价线性方法根据试验和工程类比来给定材料的阻尼比和剪切模量进而计算其动力反应，由于建模简单而被广泛用于地震工程学中，模拟地震波在岩土体中传播以及岩土体与结构物间的动力相互作用．非线性动力分析则考虑材料物理力学性质空间和时间上的非线性，模型中各个单元不同的变形破坏阶段采用不同的阻尼比和剪切模量来计算动力反应．在ＦＬＡＣ３Ｄ的动力计算中即采取非线性震动分析法，能够真实地模拟地质体的应力一应变关系．动力分析过程一般分为以下两个步骤：一定地质条件下的静力平衡计算和施加动力荷载后的动力反应分析．在第一步中，确定模型范围、初始条件、材料类型、本构模型以及模型的填筑、开挖、衬砌等，也就是静力作用下的平衡计算．第二步，是在第一步计算的基础上，施加动荷载，应考虑以下３个方面的内容：①动力加载和动力边界；②力学阻尼；③地震波在介质中的传播．２工程地质概况某边坡为云南昆明某建筑土质边坡．该边坡走向大致为ＷＥ向，坡面倾向南，呈单面边坡，边坡高度为１０ｍ，坡角为５０。
12层钢筋混凝土标准框架振动台模型试验报告(PDF)-1
弹模均值(MPa) 7.751×103
名称 铁丝
型号
20＃ 18＃ 14＃
表 4 钢筋的材性试验结果
直径 (mm)
面积 (mm2)
0.90
0.63
1.20
1.13
2.11
3.50
屈服强度 （MPa）
327 347 391
极限强度 （MPa）
397 420 560
2.5 测点布置
试验中采用加速度计、应变传感器量测模型结构的动力响应。加速度计的方向有 X、Y、 Z 三个方向。
4.1 AutoCAD 文件 ............................................................................................................. 12 4.2 输入地震波数据文件............................................................................................... 12 4.3 测点记录数据文件 ................................................................................................... 12 4.4 传递函数数据文件 ................................................................................................... 12
2 试验设计
2.1 试验装置
地震模拟振动台主要性能参数：
台面尺寸
el-centro地震波积分
EL centro 地震波的积分及相关问题1.地震波的选取本作业的EL centro地震波来源于（NISEE—the national information service for earthquake engineering University of California, Berkeley）。
具体链接为：/data/strong_motion/caltech/volume2.d/el_centro_1940_ s00e这是一条南北向的地震波波，已经过一定的修正，网站中给出了地震波的初试速度和初始位移，原文具体说明如下：CORRECTED ACCELEROGRAM IIA001 40.001.0 COMP S00E FILE 1 CORRESPONDING TOFILE 1 OF UNCORRECTED ACCELEROGRAM DATA OF VOL. I-A, EERL 70-20 IMPERIAL VALLEY EARTHQUAKEMAY 18, 1940 - 2037 PSTIA001 40.001.0 S 18 STATION NO. 117 32 47 43N,115 32 55W 38EL CENTRO SITE IMPERIAL VALLEY IRRIGATION DISTRICT50COMP S00E 9 IMPERIAL VALLEY EARTHQUAKE MAY 18, 1940 - 2037 PST 52 EPICENTER 32 44 00N,115 27 00W 31 INSTR PERIOD = 0.0990 SEC DAMPING = 0.552 42 NO. OF POINTS = 985 DURATION = 53.73 SEC 42 UNITS ARE SEC AND G/10. 23 RMS ACCLN OF COMPLETE RECORD = 0.4876 G/10. 43 ACCELEROGRAM IS BAND-PASS FILTERED BETWEEN 0.070 AND 25.000CYC/SEC2688 INSTRUMENT AND BASELINE CORRECTED DATAAT EQUALLY-SPACED INTERVALS OF 0.02 SEC.PEAK ACCELERATION = 341.69531 CMS/SEC/SEC AT 2.1200 SECPEAK VELOCITY = 33.44914 CMS/SEC AT 2.1800 SECPEAK DISPLACEMENT = 10.86678 CMS AT 8.5800 SECINITIAL VELOCITY = -4.66421 CMS/SEC INITIAL DISP. = 2.15852 CMSIMPERIAL VALLEY EARTHQUAKE MAY 18, 1940 - 2037 PST读取excel 中的EL CENTRO 地震波南北向的数据，将时间数组命名为t ，加速度数组命名为a ； plot(t,a,'r -')加速度波形图如下：图1. EL centro 地震波（加速度时程）图（单位2/mm s ）2.积分方法本次作业数据为离散的数据，下面将分别采用梯形公式以及辛普森公式进行数值积分。
基于特征结构配置的结构主动控制及仿真
第18卷第6期系统仿真学报©V ol. 18 No. 6 2006年6月Journal of System SimulationJun., 2006基于特征结构配置的结构主动控制基于特征结构配置的结构主动控制及仿真及仿真王国胜1吕强1梁冰2段广仁2(1.装甲兵工程学院控制工程系, 北京 100072;2.哈尔滨工业大学控制理论与制导技术研究中心, 哈尔滨150001)摘要考虑了具有最优控制的结构主动控制问题目的是在系统满足闭环特性的前提下设计状态反馈控制器使得系统性能泛函极小化利用特征结构配置方法提供的自由度给出了性能泛函的显示参数化表示从而该问题转化为带有约束条件的优化问题参与优化的变量仅为一组参量并给出了求解该优化问题的算法该算法直接基于结构系统矩阵故其简单性为工程应用提供方便地震作用下对三层剪切结构建筑模型进行仿真分析结果表明所提结构主动控制方法的有效性关键词结构系统主动控制特征结构配置地震控制优化中图分类号TP13; TP271 文献标识码 A 文章编号1004-731X (2006) 06-1605-04 Structural Active Control Based on Eigenstructure Assignment and Its SimulationsWANG Guo-sheng1, LV Qiang1, LIANG Bing 2, DUAN Guang-ren 2(1.Department of Control Engineering, Academy of Armored Force Engineering, Beijing 100072, China;2.Center for Control Theory and Guidance Technology, Harbin Institute of Technology, Harbin 150001, China)Abstract: The design of structural active control with minimum control effort was investigated. The aim is to design a statefeedback controller, that the closed-loop system has desired eigenvalues, and a system performance function is optimized. By utilizing design degrees of freedom offered by parametric eigenstructure assignment, a parametric expression for the system performance index was proposed. Thus the optimization problem was changed into a minimization problem with someconstraints and the optimized variables are a group of parameters. An algorithm was proposed for this minimization and utilized the original system data, and thus it is simple to use in applications. A three-story shearing model under earthquake excitation was analyzed by using the proposed algorithm and the simulation results show the effect of this algorithm.Key words: structural systems; active control; eigenstructure assignment; earthquake control; optimization.引言自从1972年美国Yao[1]结合现代控制理论提出了土木工程结构振动控制的概念开创了结构振动的主动控制研究的历程结构振动控制从理论到应用都取得了很大进展结构振动控制方法按照控制系统有无能源输入分为主动控制被动控制半主动控制和混合控制等其中主动控制是一种积极的抗振手段具有效果好适用范围广等优点成为国内外相关领域研究的前沿课题[2-3]近30年来应用和发展起来的适用于土木结构的主动控制算法主要包括二次型最优控制独立模态最优控制极点配置和滑动模态控制等极点配置或特征结构配置作为土木工程结构的主动控制算法之一虽然很早被提出但在土木结构领域中的应用却很少可查到的文献很少见文献[4]本文则把特征结构配置参数化方法[5-6]和最优控制问题相结合引入土木结构领域中考虑了具有最优控制力的结构主动控制问题其目的是设计状态反馈控制器使得闭环系统具有希望的极点外还使得系统性能泛函极小化利用收稿日期2005-04-28修回日期2006-03-03基金项目国家杰出青年基金(69925308)作者简介王国胜(1975-), 男, 河北唐山人, 讲师, 博士, 研究方向为线性系统理论结构控制理论鲁棒控制理论及应用; 吕强(1962-), 男, 黑龙江人, 教授, 博导, 博士, 研究方向为神经网络控制火力控制及应用特征结构配置方法提供的自由度给出了性能指标参数化表达式把优化问题最终转化为含有约束条件的极小化问题参与优化的变量为特征结构配置方法提供的自由参量给出了解决该优化问题的方法该方法直接基于结构系统矩阵不涉及系统增广或变换其简单性为工程应用提供了方便最后应用该算法设计了地震作用下三层剪切结构建筑模型的状态反馈控制器仿真结果表明了本文所提方法的有效性1 结构系统状态空间模型考虑在水平地震地面运动加速度)(tx g 作用下n 自由度的层间剪切型结构模型其运动方程为)()()()()(1txMItHutKXtXCtXM gn×−=++(1)式中)(tX )(tX 和)(tX 分别为各楼层相对于地面的位移速度加速度向量)(tu为r维控制力向量g x 为地震地面运动加速度M, C和K分别为结构系统的nn×阶质量矩阵阻尼矩阵和刚度矩阵H为控制力作用位置矩阵1×n I为n行元素均为1的列向量要求M和H均为满秩且矩阵对),(11HMCM−−−可控即nHMsICM n n=−−−×−][rank11C∈∀s(2)系统(1)的等价状态方程为)()()()(tx EtButAZtZg−+=(3)2006年6月 系 统 仿 真 学 报 Jun., 2006式中−−=−−C M K M I A n 110,=−H M B 10, =n I E 0,=X X Z 选取状态反馈控制律)()()(10t FZ t XF t X F u =+=][10F F F = (4)式中F 为n r 2×的反馈增益矩阵反馈控制作用是状态变量(速度和位移)的线性组合此时闭环系统为)()()(t xE t Z A t Z g c −=BF A A c += (5)2 特征结构配置控制算法因非亏损矩阵较亏损矩阵对系统参数扰动具有良好的鲁棒性故本文仅考虑闭环系统矩阵c A 的特征值为互异且自共轭情形记特征值为C ∈is n i 2,,2,1"=其对应的特征向量分别为iv ni 2,,2,1"=则有i c i i v A v s =, n i 2,,2,1"= (6)记),,,(diag 221n s s s "=Λ,][221n v v v V "= (7)则方程(6)等价于Λ=V V A c 或Λ=+V BFV AV (8)因矩阵对),(11H M C M −−−可控故对矩阵][11H M sI C M −−−−进行初等变换可得n n ×阶单模阵)(s P 和)()(r n r n +×+阶单模阵)(s Q 满足下式]0[)(][)(11I s Q H M sIC M s P =−−−−C ∈∀s (9)对)(s Q进行如下分块=)()()()()(22211211s Q s Q s Q s Q s Q (10) 其中)(11s Q 为r n ×阶矩阵从而有下述定理它给出了方程(8)中增益阵F 的参数化表达式其证明过程详见文献[5-6]关于更多的特征结构配置参数化方法及其应用研究可参见文献[7-13]定理1 给定二阶动力学系统(1)那么1) 若矩阵对),(11H M C M −−−可控则矩阵对),(B A 仍可控的充要条件是存在n n ×阶单模阵)(s H 和)()(r n r n +×+阶单模阵=)()()()()(22211211s L s L s L s L s L (11) 其中)(11s L 为r n ×阶矩阵满足下式]0[)()]()()([)(11112n I s L s Q sI K M s P s Q s H =−+−− (12)2) 若上述条件成立那么满足(8)的状态反馈增益阵F 可如下给出1−=WV F (13)式中1112211()[],()i n i ii i L s V v v v v g s L s ==" (14)][221n w w w W "=i i i i i i i g s KL M s P s Q s L s Q w )]()()()()([111222121−+= (15)其中g i , n i 2,,2,1"=是一组1×r 阶参量需满足0)()()()(det 221121111122111111≠n n n n n g s L s g s L s g s L g s L "! (16)j i s s =⇔ ji g g=n j i 2,,2,1,"= (17)综合上述特征结构配置参数化结果其优越性可以归纳为如下几点1) 该方法给出了满足方程(8)的所有状态反馈增益阵和闭环特征向量矩阵的参数化表示其含有的参量可进一步用来满足系统设计中其它性能指标如鲁棒性等2) 该方法计算过程中只涉及层间剪切型结构模型(1)中矩阵MC K 和H 并不涉及增广系统(3)中矩阵A 和B故便于工程应用3 最优结构主动控制设计本文考虑的具有最优控制力的结构主动控制设计问题可如下描述: 给定层间剪切建筑结构系统(1)以及一组自共轭且互异的复数C ∈is n i 2,,2,1"=确定形如(4)的状态反馈控制律)(t FZ u=对于任意正定对称矩阵R满足下述条件:1) 闭环系统矩阵c A 的特征值为C ∈is n i 2,,2,1"=2) ))((min F P tr ;其中正定矩阵P 是下述Lyapunov方程的解RF F PA P A Tc T c −=+. (18)不难发现上述优化问题等价于极小化下述二次型性能指标函数dt t Ru t u I T )()(0∫∞=(19)为求解问题ESA 我们首先给出如下结论给定系统(1)以及一组共轭互异的复数i s ,n i 2,,2,1"=, 若矩阵对),(11H M C M −−−和),(B A 均可控那么对于任一正定对称矩阵r r R ×∈R , 方程(18)中矩阵P 的解为122)()(21−×−+−=V s s g s RM s M g V P nn j i j j i T T i T (20) 式中)()()()()()(111222121i i i i i i s KL M s P s Q s L s Q s M −+= (21)矩阵V 由(14)决定r i g C ∈, n i 2,,2,1"=,为一组满足(16)和(17)的自由参量若记)(i g V V =由(20)易知))2,,2,1,((tr ))(tr(n i g P F P i "== (22)式中)2,,2,1,(n i g P i "=由(20)给出从而问题ESA 转化为))2,,2,1,,((tr min n i s g P i i "=,s. t . (16)和(17) (23)2006年6月 王国胜, 等基于特征结构配置的结构主动控制及仿真 Jun., 2006综上分析问题ESA的求解过程可归纳为如下步骤我们称之为特征结构配置方法(以下简称算法ESA)1) 算满足(9)的单模阵)(s P 和)(s Q 如(10)对)(s Q 进行分块2) 计算满足(12)式的单模阵)(s H 和)(sL 如(11)对)(s L 进行分块3) 设定ig n i 2,,2,1"=的参量表示根据(14)和(15)分别计算矩阵V 和W的参量表达式4) 求解优化问题(23)确定满足(16)和(17)的一组参量ig n i 2,,2,1"=将其代回上步计算矩阵V 和W5) 根据(13)计算状态反馈增益阵F4 数值仿真分析考虑如图1所示三层剪切型结构模型[2-3]该模型的结构参数取自三层Benchmark 模型但与标准Benchmark 不同的是采用在底层和中间层设置两个主动拉索控制装置结构系统矩阵为)kg (981981981=M−=001011H s/m)(N 3.43763.26.6163.27.4563.576.613.577.382⋅−−−−=C62.741 1.6410.3691.641 3.021 1.62410(N/m)0.369 1.624 1.333K −=−−× − ,假设待配置的特征值为i s 1432,1±−=is 4364,3±−=is 7296,5±−=根据算法ESA有如下结果1) 由奇异值分解易算得满足(9)的单模阵)(i s P 和)(i sQ 6,,2,1"=i 并如(10)对)(i s Q 进行分块;2) 由奇异值分解易求得满足(12)的单模阵)(i s H 和)(i sL 6,,2,1"=i 并如(11)对)(i s L 进行分块;3) 设定=i i i b ag 6,,2,1"=i 由(14)和(15)算得+−+−−−+−−−−++++++−++=−11111111111141)4539.00352.0()0658.03754.0()3566.00186.0()3050.00180.0()0661.014960.0()1538.00481.0()0315.00042.0()0266.00010.0()0246.00036.0()0210.00032.0()0023.00112.0()0098.00055.0(10b i a i b i a i b i a i b i a i b i a i b i a i v 21v v =+−++++−−++−−−+−−++−+−+=−33333333333343)1471.02036.0()0149.03493.0()0843.00867.0()1189.02153.0()0789.03285.0()0400.05438.0()0051.00027.0()0008.00080.0()0022.00016.0()0053.00020.0()0072.00028.0()0125.00008.0(10b i a i b i a i b i a i b i a i b i a i b i a i v43v v =−+−−+−−+++−+−++++−+−−−=−55555555555545)0741.02663.0()0496.01058.0()5560.00227.0()0645.02279.0()0799.04484.0()1424.02206.0()0035.00014.0()0009.00014.0()0076.00006.0()0030.00013.0()0063.00003.0()0033.00016.0(10b i a i b i a i b i a i b i a i b i a i b i a i v65v v =×−+−−−×−=−−−−1919191101)107765.21()2324.31087.3(10)1088.32297.3(10)108036.91(b i a i b i a i w 21w w =×++−−−×+=−−−−393113113103)102100.11()529.523967.8(10)4043.8529.52(10)109636.11(b i a i b i a i w 43w w =×+++−+×+=−−−−595115115105)102100.11()529.523967.8(10)529.524043.8(10)109636.11(b i a i b i a i w 65w w =4) 由matlab 优化工具箱中函数fmincon 算得优化问题(23)的优化参数为++==i i g g 9283.99984.91414.88380.721,−−+==i i g g 6354.44319.29593.61988.343,+−−==i i g g 0498.45629.71648.28124.365 将其代入第三步易得矩阵V 和W5) 根据(13)算得状态反馈增益阵为50.86260.95530.54900.13530.04420.0198100.4901 1.43260.77330.08910.11590.0786F −−−−− =−−−为进一步验证算法有效性选取输入地震波为El Centro (S00E)波图23和4给出了无控和F 控制结构系统的各层位移速度和加速度反应曲线图5给出了相应的控制力时程曲线仿真结果表明El Centro 地震输入下本文所提算法对结构的位移速度和加速度响应均能起到良好的控制作用同样的El Centro(S00E)地震波输入下图2表明采用本文设计控制律的结构位移要远小于无控下的结构位移图3表明采用本文设计控制律的结构速度要远小于无控下的结构速度图4表明采用本文设计控制律的结构加速度要远小2006年6月 系 统 仿 真 学 报 Jun., 2006图2 Elcentro 波作用下无控和F 控制的结构位移比较 图3 Elcentro 波作用下无控和F 控制的各层速度比较图4 Elcentro 波作用下无控和F 控制的各层加速度比较 图5 Elcentro 波作用下F 控制的各层控制力时程于无控下的结构加速度图5表明得到上述很好的控制效果却所采用了较小控制力5 结论将现代控制理论中的特征结构配置方法引入土木结构中考虑了结构系统的最优控制问题基于特征结构配置参数化方法提供的自由度将该问题转化为含有约束条件的优化问题并给出了一种简单有效的算法最后把利用该算法设计的控制器应用于地震作用下的三层剪切结构建筑模型并进行了仿真分析仿真结果表明El centro 地震输入下本文所提算法对结构的位移速度和加速度响应均能起到良好的控制作用同时也表明该特征结构配置方法在实际应用中的简单且有效性其在土木工程中的进一步应用将是今后研究工作的重点参考文献[1] Yao J T P. Concept of structure control [J]. ASCE Journal of StructureDivision (S0733-9453). 1972, 98(ST7): 1567-1574. [2] 欧进萍. 结构振动控制:主动半主动和智能控制 [D]. 北京: 科学出版社, 2003.[3] 张春巍, 欧进萍. 结构振动控制Benchmark 研究发展综述 [C]//现代土木工程理论与实践, 南京: 河海大学出版社, 2003: 489-496. [4] Mohamed A R, Horst H L. Structural Control by Pole AssignmentMethod [J]. Journal of the Engineering Mechanics Division (S0044- 7951). 1978, 104(5): 1159-1175.[5] Duan G R, Liu G P . Complete parametric approach for eigenstructureassignment in a class of second-order linear systems [J]. Automatica (S0005-1098). 2002, 38(4): 725-729.[6] Duan G R, Wang G S, Liu G P. Eigenstructure assignment in a class ofsecond-order linear systems: a complete parametric approach [C]// Proceedings of CACSCUK, Manchester, UK, 2002, 89-96.[7] Wang G S, Duan G R. Robust pole assignment via P-D feedback in aclass of second-order dynamic systems [C]//International Conference of Automation, Robots and Computer Vision, Kunming, China, 2004, 1152-1156. [8] Wang G S, Liang B, Duan G R Reconfiguring second-order dynamicsystems via state feedback eigenstructure assignment [J]. International Journal of Control, Automation, and Systems (S1598-6446). 2005, 3(1): 109-116.[9] 王国胜, 柴庆宣, 段广仁. 二阶动力学系统分散状态反馈特征结构配置 [J]. 哈尔滨工业大学学报, 2004, 36(12): 1585-1588. [10] 王国胜, 段广仁. 二阶动力学系统输出反馈特征结构配置 [J]. 系统工程与电子技术, 2004, 26(8): 1080-1083.[11] 王国胜, 王子华, 段广仁. 二阶动力学系统的干扰解耦与抑制 [J].系统工程与电子技术, 2004, 26(11): 1640-1643.[12] 段广仁, 王国胜. 矩阵方程EVJ 2+A VJ+CV=BW 的两种解析通解[J]. 哈尔滨工业大学学报, 2005, 37(1): 1-4.[13] Duan G R. Parametric eigenstructure assignment in second-orderdescriptor linear systems [J]. IEEE Transactions on Automatic Control (S0018-9286). 2004, 49(10): 1789-1795.。
地震工程学-反应谱和地震时程波的相互转化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生成的人造地震波如图所示。
东南大学第十七届结构创新竞赛暨南京高校邀请赛
“东南大学第十七届结构创新竞赛暨第七届南京高校邀请赛〞加载组 A 竞赛题目及细那么一、根本竞赛规那么本次加载组竞赛所施加荷载为利用大型振动台模拟的实质水平川震作用。
参加加载组竞赛的每一个队伍需要设计并制作一个房屋模型，用以承受附加铁块载重，并抵抗振动台所产生的人工地震作用。
在正式加载时，参赛模型将被分批置放于振动台上，并施加假设干级不同样大小的地震作用。
地震作用将先由小地震开始，逐级加大，最大的地震加速度将到达1600gal (1600cm/s 2 )。
地震加载如右图所示。
全部参赛队伍的模型自己质量(M M)、承载铁块的加权重量 (W e)，与测试后之模型能承受的最大地震力大小(I)(详见模型损坏准那么)，都将被记录下来，用以计算每个参与队伍模型的效率比(efficiency ratio ，详尽计算方法见赛题第四、五局部)。
图 1加载表示图依照模型的结构造型与系统、率比结果三个方面进行综合评分模型制作工艺和加载效(详尽评分细那么见赛题第七局部)，依照综合评分排序决出加载组竞赛获奖队伍。
二、模型资料(1)木材：用于制作结构构件。
尺寸：长度 1000mm ，截面有 50mm× 1mm、50mm× 2mm 、42mm × 2mm、6mm× 6mm；性能参照值：顺纹弹性模量×10 MPa，顺纹抗拉强度30MPa。
(2)502 胶水：用于构件之间的连接、模型与底板的连接。
(3)热熔胶：用于铁块的固定。
(4) 模型固定底板：厚度约 6mm 的木板，长与宽分别为 29cm 与 29cm。
底板上除组委会预设的孔洞外，严禁参赛队伍另行钻孔。
三、模型要求3.1 几何尺寸要求(1)根本结构 : 模型必定最少包括一个房屋的根本骨架，可采用框架结构系统或外框架-内核心筒结构系统。
全部模型只能采用竞赛供应的资料制作。
(2)模型大小：模型总高度不得高出 75cm。
总高度为底板顶面至模型顶面的垂直距离。