1          引言

土木工程结构在其漫长的使用寿命中，会遇到各类极端灾害的作用。由于土木工程结构自身体量庞大、造价高昂、结构复杂，完全依赖物理试验手段研究其灾变过程难度很大。即使采用缩比模型，也依然存在尺寸效应、相似比设计等诸多困难。与此同时，土木工程也是城市功能的主要载体，当工程结构的灾变研究从单体发展到城市区域规模时，采用物理实验手段更是无能为力。因此，计算机数值模拟作为一种重要的科学研究手段，在工程结构抗震防灾领域得到日益广泛的应用。

工程结构计算机数值模拟的核心工作，是将工程结构的各种复杂行为（力学、热学等），建立相应的数学方程，而后通过计算机对这些方程进行求解，以预测相应的工程结构响应。工程结构数值模拟包括三个主要的构成部分：

(1) 工程结构的数值计算模型；

(2) 数学方程的求解算法；

(3) 完成工程结构数值模拟所需的计算机硬件平台。

其中高性能的硬件平台是基础，高效的求解算法是重要手段，而工程结构的数值计算模型是核心研究内容。

由于土木工程结构自身体量的庞大和行为的复杂，使得准确描述其复杂非线性受力行为的数值模型的计算量非常大。例如，虽然研究早已发现，实体单元是描述三维物体受力行为最为合适的单元类型，但是基于实体单元的建筑结构非线性计算，即便是利用当前最先进的计算平台，也只能完成一些简单的多、高层单体结构的非线性分析。例如，Yamashita等[1]利用实体单元建立了一座高129.7m的规则高层钢结构的计算模型，单元总数超过了1600万，这样的计算量，即便是对于E-simulator这样的超级计算机，也是一个非常有挑战性的工作。计算机有限的计算能力与工程结构数值模拟几乎无限的计算量需求，构成了工程结构数值模拟的一个主要矛盾，同时也成为工程结构数值模拟不断进步的一个重要源动力。

1.      基于精细实体单元的结构非线性分析[1]

实际上，科学研究/工程应用的需求与试验能力（包括物理试验和数值模拟试验）的限制之间的矛盾，无论是对于物理试验还是对于数值模拟试验都同样存在。然而，计算机技术日新月异的发展，为突破数值模拟计算能力限制不断提供新的手段。与此相对的，物理试验能力的发展却遇到了巨大的困难。例如，以振动台试验为例，目前世界上最大的振动台为1995年落成的日本E-Defense振动台。此后近20年，振动台试验能力都很难进一步得到提高。而世界上最快的超级计算机的头衔，几乎每年都在变化。甚至家家户户使用的桌面电脑的速度，已经可以和15年前世界上最快的超级计算机相媲美。日趋庞大而廉价的计算机数值计算能力，正沿着摩尔定律飞快发展，并不断为工程结构的数值模拟提供强有力的推动力。

世界各国的研究者，都从高性能计算技术的飞速发展中看到了工程结构数值模拟的美好前景。很多研究团队都做出了非常出色的研究成果。例如，日本东京大学在工程结构精细化数值模拟、城市区域震害数值模拟等方面做出了非常出色的研究成果，建立了IES等先进震害模拟系统。美国加州伯克利大学等领导开发的OpenSees计算程序结合NEEShub高性能计算平台，对促进工程结构数值模拟起到了重要的推动作用。其他例如Carnegie Mellon大学的J. Bielak教授等，在这一领域也都做出了非常出色的研究成果，极大的深化了人们对工程震害机理的理解及防灾减灾对策的研究。2009年，中国住建部、美国国家科学基金会（NSF）和日本文部省组织中美日三国50余位地震工程和结构工程知名专家在广州举行“中美日建筑结构抗震减灾研讨会”，探讨未来结构工程和地震工程的发展方向。与会专家一致认定：“基于超大规模计算的区域综合震害预测是未来地震工程领域具有重大价值的研究方向”。

本文将结合作者近年来在工程结构和城市区域震害数值模拟及防震减灾对策方面的科学研究和工程应用，介绍在计算模型、求解技术、硬件平台等方面的一些最新成果，为国内外相关领域的研究提供参考。

Equation Chapter (Next) Section 1Equation Chapter (Next) Section 1

2          重要工程地震灾变模拟

超高层建筑和超大跨桥梁等重大基础工程设施的建设，不仅反映了一个国家或地区的经济繁荣和社会进步，也是一个国家或地区建设科技发展水平的重要标志。随着我们城市化进程的不断推进，超高层建筑和超大跨桥梁等重大基础工程设施的建设势头强劲。然而，重大工程结构体系庞大，单元种类繁多，当面临强震作用时，其动力灾变效应无疑将更为复杂；同时，重大工程结构的损伤和破坏会给整个社会带来巨大的负面影响。因此，确保重大工程的地震安全是地震工程界的核心任务之一。而要保障重大工程的地震安全，首先要深入理解重大工程在强震作用下的损伤演化与倒塌灾变过程。故本节将通过建立高效的重大工程数值分析模型、提出倒塌灾变模拟方法，实现重大工程结构在强震下的损伤演化、灾变机理及倒塌机制的预测与评估，为重大工程抗震设计理论的进一步完善提供参考。

2.1   计算模型和算法

强大的软硬件计算平台是模拟工程结构地震灾变行为的基础，而合理的计算模型是高效准确的得到工程结构地震灾变行为的关键。目前工程结构灾变行为模拟的实现手段主要有两大类：

(1) 基于成熟的商用有限元软件平台，利用其二次开发功能，开发针对性的构件或材料数值模型，以满足工程结构地震灾变模拟的需求。其优势在于，成熟的商用平台具有稳定且强大的求解能力以及方便的前后处理功能，研究人员可以专心于开发专用的构件和材料模型，进而提高研究工作的效率，特别是在解决一些工程问题时，利用成熟的商用有限元软件平台无疑是一个很好的选择。本课题组以通用有限元软件MSC.Marc为基础，在材料、构件及单元生死准则上开展了一系列的研究工作，并在多个重要工程结构分析中得到了成功的应用。采用商用有限元软件平台的一个主要不足之处在于，商用有限元软件是一个封闭的“黑盒子”，一些内部机理及深入的开发受到很多限值。

(2) 利用具有针对性的开源有限元软件平台，并开发完善其特定功能，以满足工程结构地震灾变模拟的需求。典型的开源有限元软件平台即OpenSees，其优势在于，其绝大部分源程序都是公开的，研究人员可以深入研究其程序执行机理并进行一些深层次的开发。当然，其缺点在于其稳定性、易用性与商用有限元软件相比还有一定的差距。本课题组在OpenSees开源有限元软件平台上，开发了相应的分层壳单元模型及混凝土本构模型，使其可以实现重大工程结构地震灾变的模拟。并进一步利用其可深层开发的特点，加入了基于GPU的高性能矩阵求解算法，从而使得其计算效率得到了极大的提高。

2.1.1   纤维梁和分层壳模型

在建筑结构的地震响应分析中，学者们提出了很多数值模型来模拟梁、柱和墙等构件。如利用集中塑性铰模拟梁、柱构件[2]；三垂杆[3]和多垂杆模型[4]模拟剪力墙构件等。为了能够准确把握结构和构件在地震作用下发生灾变前后整个过程的非线性响应，特别是柱中复杂的轴力-弯矩的耦合，剪力墙中平面内外弯曲、平面内弯矩、剪力耦合的非线性行为，本研究将采用纤维梁单元模拟梁柱构件，分层壳单元模拟剪力墙构件，两个模型的简要介绍如下。

纤维梁单元已经被广泛地运用到了弯曲破坏的梁柱构件地震响应模拟[5, 6]。在纤维梁模型中，梁柱截面被离散成若干个纤维(如图2.1a所示)，每个纤维可以赋予不同的单轴应力-应变关系，同一截面上的所有纤维满足平截面假定。纤维梁模型能够很好地模拟梁柱构件中轴力和弯矩的耦合作用，同时也能适用不同的截面形状。此外，通过赋予核心区混凝土纤维单轴约束本构，可以方便的考虑柱截面中箍筋的约束作用。而梁柱构件中的剪切破坏将以弹脆性破坏模型来近似考虑，即：当构件的剪力超过抗剪承载力时，构件的发生剪切破坏，其抗剪强度和刚度设定为0。

1.      纤维梁截面及材料本构模型

纤维梁模型中，对于混凝土材料，本课题组在Légeron & Paultre模型[7]的基础上加以改进，考虑了约束效应、裂面效应、滞回效应等影响[8]，其单轴应力应变关系曲线如图2.1b所示，可以较好反映约束效应、软化行为，以及反复受力下的滞回和刚度退化的特性；对于钢筋材料，在Esmaeily & Xiao等提出的模型[9]基础上进行改进，可反映钢筋单调加载时的屈服、硬化和软化现象，并合理考虑了钢筋的Bauschinger效应[8]（如图2.1c）。

基于以上的材料模型，通过UBEAM子程序的二次开发将纤维梁模型内嵌到了通用有限元软件MSC.Marc中，通常每个混凝土梁柱构件被分成6~8个单元来保证分析结果的精度。大量分析结果表明，该纤维梁模型的模拟结果与试验吻合良好，能较好把握梁柱构件的弹塑性行为[10-13]。为了进一步验证钢筋混凝土柱临近倒塌状态，强度和刚度迅速退化时纤维梁模型的准确性，对唐代远等[14]人进行的某一具有明显退化现象的钢筋混凝土压弯柱的滞回性能进行了模拟。构件的主要尺寸和配筋如图2.2所示，构件的纵向配筋率约为1.29%，轴压比为0.348。滞回曲线的比较如图2.2所示，表明该纤维梁模型能较好的考虑构件临近倒塌状态时的强度和刚度退化现象。

2.      纤维梁模型模拟结果与试验结果比较

本研究中采用分层壳模型模拟剪力墙。分层壳模型基于复合材料的基本理论，可以较好地模拟剪力墙面内、面外的弯曲耦合以及面内的弯曲剪力耦合。与纤维梁的原理类似，分层壳模型将剪力墙沿厚度方向分成若干不同厚度的层，每一层可以具有不同的材料属性，同一截面的所有层依然满足平截面假定。通过将剪力墙中的纵横向钢筋网弥散成等效钢筋层来模型剪力墙中的纵横向配筋，如图2.3所示。对于剪力墙中的边缘约束构件，将纵向配筋离散成一系列杆单元进行模拟，通过共节点的方式保证两者的变形协调。大量的计算分析表明该模型可以较好模拟剪力墙构件的受力特点[10, 15]。同样，为了进一步验证分层壳模型在剪力墙临近倒塌阶段的适用性，对一具有明显承载力退化现象的矮墙滞回试验进行了模拟对比。该矮墙高1050mm，宽1100mm，厚度约为120mm，构件的具体尺寸和加载示意图如图2.4所示[16]。该构件滞回曲线的对比如图2.5所示，可见，分层壳模型仍能较好地反映剪力墙临界倒塌阶段的承载力迅速退化现象。

3.      分层壳示意图

4.      钢骨混凝土剪力墙体主要尺寸、配筋及材料属性

5.      典型剪力墙滞回曲线模拟与试验比较

OpenSees作为一款开源的有限元软件，在建筑结构抗震研究领域得到了广泛的运用。OpenSees程序中的单元以纤维梁柱模型为主，该模型能较好的考虑双向弯矩和轴力的耦合作用，但不能方便准确的考虑剪力墙的受剪特性。本课题组基于分层壳理论，在OpenSees中开发了二维混凝土本构模型和剪力墙分层壳模型。本研究提出采用基于损伤力学和弥散裂缝模型来模拟分层壳单元中混凝土的二维受力行为。该模型具有形式简单、计算稳定性好等优点，适合作为基本模型集成在开源程序中供研究人员进一步发展和完善[17]。

混凝土的二维本构模型的基本方程可以表示为：

                              \* MERGEFORMAT (2-1)

式中，、分别为主应力坐标系下混凝土的应力和应变，D1、D2为混凝土在主应力坐标系下的损伤标量，受拉与受压损伤分开考虑。其中受压损伤参考Løland建议的受压损伤演化曲线[18]，受拉损伤参考Mazars建议的受拉损伤演化曲线[19]。对于剪力墙内的分布钢筋，本研究基于OpenSees中现有的钢筋模型开发了相应的多维钢筋模型。

结合各层厚度和各层所使用的多维钢筋、混凝土材料可生成RC剪力墙分层壳截面，进而将该截面赋予4节点壳单元完成分层壳剪力墙的定义，其流程如图2.6所示。

6.      OpenSees中分层壳模型框架

2.1.2   生死单元技术

建筑结构在特大地震作用下发生倒塌破坏的过程中，构件会逐渐破坏失效而退出工作。在有限元分析中，这一过程可以采用生死单元技术来模拟。生死单元技术是通过修改单元的刚度矩阵来实现的，当某一单元的响应超过某一设定的阈值时将会被“杀死”，通过调整该单元的刚度矩阵，赋予被“杀死”单元一个极小的刚度，避免直接归零带来的刚度矩阵奇异的问题。当该单元被杀死，其相应的应力和应变将释放为0，在整个结构的质量和阻尼矩阵中该单元所对应的元素也将归为0。在MSC.Marc通用有限元分析中，可以方便的采用UACTIVE子程序接口，通过编写相应的代码来实现这一过程。

由于纤维梁和分层壳模型都是基于材料本构层次的数值模型，可以建立基于材料层次的单元失效准则来监控整个地震灾变过程各构件的应力-应变行为。对于纤维梁单元，每个单元至少包含36个混凝土纤维和4个钢筋纤维，而每个纤维沿长度方向有3个高斯积分点；类似，对于分层壳模型，每个壳单元至少分成10层（具体层数根据剪力墙的实际配筋情况而定），每层具有4个高斯积分点。当每一个纤维或者层上任意积分点的应变值超过了预设的失效准则，该纤维和层的应力和刚度将被杀死，当某一单元所有的纤维或者层均被杀死，则认为该构件完全失效，将其从整个模型中杀死并移除。

2.1.3   基于GPU的高性能矩阵并行求解技术

随着有限元计算在土木工程领域的不断应用，越来越多的工程采用有限元计算进行抗震分析，其中不乏超高层建筑、大跨建筑等，从而导致对于有限元软件计算能力的需求也日渐提高。现如今，并行计算技术已经成为了有限元计算加速的主流，而日益发展的GPU技术也将并行计算加速提升到了一个新的高度。OpenSees作为开源有限元软件，虽然拥有良好的可扩展性，但是计算效率一直以来都是它的瓶颈。通过分析OpenSees求解大型模型的时间分布，可以发现：影响OpenSees计算效率的主要因素是线性方程组求解模块。因此，本研究采用GPU加速技术，为OpenSees编写了线性方程组GPU加速求解器模块，将OpenSees的计算效率大幅提高，满足大型模型地震下动力时程分析的时效性要求。

在OpenSees中，求解线性方程组需要两个基本类：

(1)              线性方程组类（LinearSOE）。该类用于存储、集成和读取线性方程组。对于一个特定的线性方程组Ax=b，LinearSOE将A、x、b分别储存，并包含了相应的数据读取方法，使之可以被其他友元类调用。随着A的矩阵类型不同，LinearSOE派生了许多子类，其基本类型分为一般矩阵、带状矩阵、分块矩阵和稀疏矩阵等等。

(2)              线性方程组求解器类（LinearSOESolver类）。该类用于特定的LinearSOE求解。根据求解方式的不同，分为直接法和迭代法两类。一种LinearSOESolver往往仅能求解一种LinearSOE，然而一种LinearSOE可以采用多种LinearSOESolver进行求解。

在结构动力方程中，刚度阵、质量阵和阻尼阵往往都具有一定的稀疏矩阵的特征，因此，动力方程迭代求解所采用的线性方程组Ax=b中，矩阵A也具有稀疏矩阵的特征，适合采用稀疏矩阵求解。同时，稀疏矩阵可以有效减小数据所需存储空间，适合超大规模模型的计算。在原版OpenSees中，针对稀疏矩阵，有SparseSYM，SuperLU等许多LinearSOESolver类可供选择。

7.      GPU加速OpenSees方程组求解流程

图2.7为GPU加速OpenSees方程组求解流程。其求解模块与OpenSees原版主要有以下几个区别：

(1)              在CPU线程中集成矩阵，拷贝到显存中，再并行计算。在OpenSees原版的CPU计算模块中，LinearSOE集成方程组之后，LinearSOESolver直接调用LinearSOE存储于内存中的数据进行计算。而在GPU加速求解模块下，矩阵集成将在CPU线程中完成，之后LinearSOESolver会先将方程组矩阵和向量拷贝至显存中，再调用显卡并行求解方程组。

(2)              采用迭代法计算。原版OpenSees中，对于稀疏矩阵求解，多数采用直接求解法进行计算。直接求解法逻辑计算很多，并行度低，并不适合GPU计算。相对的，采用迭代法计算则可以很好地发挥GPU并行计算的能力，提升计算效率。因此，在GPU加速求解模块中采用迭代法进行稀疏矩阵线性方程组求解。

此外，本研究引入两个基于GPU加速的稀疏矩阵方程组求解库，用于OpenSees中稀疏矩阵方程组加速求解：

(1)              CulaSparse. CulaSparse是一个基于GPU加速的线性代数函数库[20]，用于迭代求解稀疏矩阵方程组。它由EM Photonics公司开发和维护，基于英伟达公司（NVIDIA）开发的CUDA（Compute Unified Device Architecture，统一计算设备架构）编写。它支持许多常用的稀疏矩阵方程组求解器、预处理器以及存储格式等等，其加速效率最高可达到10倍以上；此外，CulaSparse提供C、C++和Fortran语言接口，支持Linux、Windows、OSX等操作系统。CulaSparse为研究人员免费提供研究许可，目前最新的CulaSparse版本为CulaSparseS5，支持NVIDIA CUDA 5.0。

(2)              CuSP. CuSP是一个开源的C++稀疏矩阵函数模板库[21]，可以进行多种稀疏矩阵运算。CuSP同样基于CUDA编写，它有着方便、快捷的程序结构和函数调用接口，容易上手；同时，由于其以模板库的形式与CUDA合并，因此其编译较为方便，仅依赖CUDA库即可。目前最新的CuSP版本为0.4.0，支持NVIDIA CUDA 5.5。

8.      求解器CulaSparseSolver/CuSPSolver UML类图

通过调用以上两个GPU加速求解库，可以快速进行稀疏矩阵方程组的求解。本研究基于OpenSees架构，为OpenSees编写了CulaSparseSolver和CuSPSolver两个求解器（均为SparseGenRowLinSolver的派生类，对应于SparseGenRowLinSOE）。图2.8为这两个求解器的UML类图。由于它们继承自OpenSees自身的求解器基类，因此调用简单方便，仅需要在OpenSees命令流中，替换掉System部分的命令流即可。

采用上述配置对OpenSees的两个GPU加速求解器性能进行了测试。测试平台参见表2.1，二者价格基本相当。其中Intel Core i7-3970X是2013年市场上可以买到的最快的CPU型号。三种求解器的计算用时间比较见表2.1，顶点时程曲线比较如图2.9所示。可以看出，采用两种GPU加速求解器进行弹塑性时程分析，可以使计算效率提升10~15倍；同时，计算结果几乎完全重合。由此可以证明，采用GPU对OpenSees弹塑性时程分析进行加速是可靠且高效的。

1.      测试平台硬件配置及求解器

2.      三种求解器计算用时

9.      三种求解器顶点时程曲线结果

2.2   超高层建筑地震灾变模拟

2.2.1   上海中心大厦简介

上海中心大厦位于上海陆家嘴，是一栋以甲级写字楼为主的综合性超高层建筑(图2.10a)，主体塔顶建筑高度632m，结构屋顶高度约580m，共124层，采用“巨柱－核心筒－伸臂桁架”的混合抗侧力体系(如图2.10b)，该体系的主要组成如下：

(1) 核心筒主体为一个边长约30m的方形钢筋混凝土筒体，核心筒底部翼墙厚1.2m，随高度增加核心筒墙厚逐渐减小，顶部厚0.5m；核心筒内腹墙厚度由底部的0.9m逐减薄至顶部的0.5m。由于建筑功能的要求，核心筒的角部在第五区以上被逐步切去，最终形成一个十字形核心筒[22]。

(2) 巨柱系统由12根型钢－混凝土巨柱组成[22]，最大柱截面约为5300mm×3700mm。8根巨柱贯穿整个结构高度，柱截面尺寸随着高度的增加逐渐减小，最终减小为2400mm×1900mm；其余4根角柱仅延伸至结构第5节段。

10.    上海中心大厦模型图

(3) 桁架系统位于结构的加强层位置，由环形桁架和伸臂桁架共同组成，高度约为9.9m，所有桁架杆件均为工字型截面钢梁。

对于上海中心大厦中的剪力墙，采用分层壳单元进行模拟，外框架、环向桁架、伸臂桁架以及塔顶均采用工字型钢梁，均采用纤维梁单元进行模拟，对于巨型钢骨混凝土柱，采用了分层壳和杆单元的组合式模型进行模拟，并以巨柱构件的精细实体有限元模型为基准，通过纯压、纯弯、单向压弯、双向压弯等多种工况的数值试验对组合式模型进行了验证和标定。更详细的模型信息及建模方法见文献[23]，最终完成的上海中心大厦的有限元模型如图2.10c所示。

2.2.2   典型地震倒塌过程

上海中心大厦通过精心设计和反复论证，已经可以保证在设计罕遇地震下的安全性[22]，但是为了更好的研究此类结构的倒塌机理和倒塌过程，本研究中人为加大地震动强度，直至结构发生倒塌。虽然这样大的地震在上海地区基本不可能发生，但是由此得到的结构倒塌模式和破坏机理，对认识超高层结构体系受力特性有一定的科学意义。

上海中心大厦倒塌分析时，以科研中广泛采用的1940年的El-Centro地震动记录作为典型输入，采用PGA调幅方法将其PGA调至19.6m/s2，约为实际地震动的6.4倍(实际地震动的PGA=0.313g)，且沿x方向单向输入。倒塌分析采用经典的Rayleigh阻尼，阻尼比取5%，结构最终倒塌模式如图2.11所示。

11.    上海中心大厦典型倒塌模式

其倒塌过程如图2.12所示：当t=2.58s时(图2.12a)，核心筒的部分连梁发生破坏，且由于核心筒从第6节段到第7节段的转变中，洞口布局有改变，核心筒开洞从中间移至外缘，存在刚度的突变，极大削弱了核心筒的抗弯能力，因此，第7节段底部核心筒外缘的剪力墙被压溃；当t=3.90s(图2.12b)时，由于核心筒从第4节段到第5节段截面形状改变，四个角部被切除，因此第五区段底部剪力墙开始被压溃；当t=5.88s(图2.12c)时，由于第五区段部分剪力墙大量破坏，内力产生重分配，巨柱受到的侧向力和竖向力逐渐增大，巨柱开始压弯破坏；当t=6.18s时(图2.12d)，第五区段核心筒和巨柱完全破坏，结构倒塌开始发生。

12.    上海中心大厦在El-Centro波下的潜在倒塌过程(PGA=19.6m/s2)

13.    上海中心大厦在El-Centro波作用下灾变位移响应

El-Centro地震动输入下结构顶点的位移时程曲线和楼层位移分布如图2.13所示。由于该高层结构的第一、二阶平动周期非常长，相应的地震力也比较小，因此，在地震作用下的破坏主要由高阶振型(水平向三阶振型)控制。故而结构临近倒塌时，结构变形模式呈高阶振型形状。破坏部位以上结构的重心并未有显著偏移(图2.13b)，导致结构倒塌以竖向倒塌模式为主，而非侧向倾覆倒塌。

显然，该超高层结构在El-Centro波作用下，第五区段以上部位破坏比较严重，主要集中在第五、六、七区段，而最终在第五区段发生折断，整个结构断成两截，可见第五区段是引起结构倒塌的潜在薄弱部位，在设计过程中应该给与更多关注。在上述倒塌过程中，连梁最先发生破坏，是结构的第一道防线，耗散部分地震能量，随后核心筒内的剪力墙和巨柱开始发生破坏，只有在同一区段的剪力墙和巨柱都发生破坏后，结构才可能发生倒塌。而通过常规弹塑性分析得到的初始屈服部位可能并非是引起结构倒塌的关键部位，如果设计不当，对初始屈服部位进行加强，可能使得耗能构件不能充分耗散能量，反而使结构变得不安全，这进一步说明了倒塌分析的重要性。

2.2.3   关键构件应力-应变曲线

结构倒塌过程是结构在地震作用下的宏观响应，为了进一步研究结构倒塌过程中关键构件力学行为，本节对关键构件的应力-应变行为进行了研究。在El-Centro地震动单向输入下计算得到的结构典型破坏部位如图2.11所示。选取破坏部位的典型结构单元(图2.14)，得到相应的倒塌过程中混凝土和钢材应力-应变滞回关系如图2.15所示。

14.    典型单元分布

对于巨柱单元，其在地震输入下受力过程的混凝土应力-应变滞回曲线如图2.15a所示。显然，巨柱在整个破坏过程中，主要承受压应力，仅在较少时刻出现了拉应力，最终混凝土被压溃，巨柱发生压弯破坏；对于连梁，其在地震输入下受力过程的剪力与剪应变的关系如图2.15b所示，最终发生剪切破坏；核心筒和伸臂桁架在地震输入下受力过程的应力-应变关系如图2.15c~f所示。

15.    关键构件灾变全过程力-变形曲线

2.3   特大跨桥梁地震灾变模拟

2.3.1   琼州海峡大桥简介

琼州海峡位于广东雷州半岛和海南岛之间，东西长约80km，南北宽约30km，平均水深44m，最大深度114m。本研究的主要对象是琼州海峡大桥前期方案中的斜拉桥方案之一，该方案全桥总长4304m，主跨1500m，中塔标高46m，边塔标高386m，大桥两端各设置4个辅助墩，墩高62m，航道净空为1270m×81m，结构总体布置如图2.16。

16.    琼州海峡大桥斜拉桥设计方案立面布置图（单位：cm）

琼州海峡大桥主要包括桥塔（墩）、桥面板、拉索三个部分。调查表明，地震中斜拉桥的破坏主要集中在桥台、桥墩、支座、主梁、结构基础等部位，而主梁的破坏也主要由于墩、台的滑移或者倾斜造成[24, 25]，因此准确模拟桥塔（墩）在地震下的响应是大跨斜拉桥地震分析的关键。琼州海峡大桥设计方案采用空心桥塔，研究中可以采用分层壳单元模拟。为验证该单元的合理性，本研究对文献中的空心桥塔试验进行了分析[26]，构件沿高度方向划分20段，每一面沿长度方向划分4段，单元沿厚度方向划分16层，截面配筋及有限元模型如图2.17所示，计算结果如图2.18所示。结果表明，MSC.Marc及OpenSees的分层壳单元计算结果均与试验均吻合良好，有较高的计算精度，可以采用该单元模拟琼州海峡大桥的桥塔构件在地震下的受力行为。

17.    桥墩试验截面配筋及有限元模型示意图[26]

18.    桥墩试验模拟结果（实验数据[26]）

由于OpenSees缺少友好的前后处理界面，在处理大型复杂结构时建模工作量极大，因此本研究利用已经发展相当成熟的有限元软件MSC.Marc对该斜拉桥进行建模，通过转换程序获得对应的OpenSees模型，并利用MSC.Marc的计算结果与OpenSees进行对比，检验OpenSees在计算大型桥梁时的可行性及可靠性。在此基础上，建立MSC.Marc与OpenSees中单元及连接的对应关系，编写转换程序完成Marc模型到OpenSees模型的转化（如图2.19所示）。最终的有限元模型如图2.19所示。

19.    琼州海峡大桥有限元模型示意图

2.3.2   典型地震倒塌模拟

琼州海峡大桥的模态分析表明，结构顺桥向和横桥向的一阶周期分别为15.7s和8.5s，因此，进行强震倒塌分析时，选择了科研中常用的长周期分量较高的ChiChi 1999及Tōnankai 1944进行地震波输入，计算结果如图2.20所示。

20.    琼州海峡大桥地震倒塌模拟结果

结果表明，在上述两条地震波作用下，琼州海峡大桥均在主塔的下横梁与桥墩交界处发生了混凝土的压坏，进而触发了结构的整体倒塌，从主塔顶点位移时程曲线也可以看出，塔顶沿着一个方向发生倾斜。

2.4   灾变模拟的工程价值

倒塌模拟作为一种工程结构地震响应分析的重要手段，不仅可以对重大工程在强震下的损伤演化和倒塌灾变全过程进行预测，判断结构的薄弱部位，还能有助于完善重大工程的抗震设计理论与设计方法。

2.4.1   超高层建筑地震动强度指标

地震动强度指标是基于性能设计方法中的重要组成部分，是联系地震动危险性与结构响应的桥梁。合理的地震动强度指标能够有效减小结构响应预测结果的离散性，众多学者已经对适用于常规结构的地震动强度指标展开了大量的研究[27-30]。但是这些研究和提出的地震动强度指标基本都是针对0~4s的常规结构而展开的，而对于大量涌现的超高层建筑，这些地震动强度指标的适用性待考证，进一步提出适用于超高层建筑地震动强度指标也是十分必要的。

与常规结构相比，超高层建筑结构的显著特点就是基本周期很长，在地震动作用下结构的响应中高阶振型参与明显，结构的破坏模式甚至受高阶振型控制[23]。因此，适用于超高层建筑结构的地震动强度指标必须能够反映结构响应中高阶振型的影响；其次，作为一个合理的地震动强度指标，还要形式简单，便于工程人员所接受，现行规范的抗震设计中都是基于加速度反应谱的，且地震危险性的衰减关系也是基于谱加速度的，故本研究中提出的地震动强度指标将以谱加速度Sa作为基本参数；最后，由于超高层建筑结构都采用基于性能设计，性能水准比较高，完全满足规范中所规定的大震需求，且具有较高的安全储备[22]。在大震下，结构的主要抗侧力构件基本保持弹性或者只有少量的屈服，整个结构进入非线性的程度并不深，因此，本研究在提出新的地震动强度指标时，初步不考虑结构进入非线性周期变长带来的影响。鉴于上述原因，本研究提出的适用于超高层建筑的地震动强度指标的一般形式如式\* MERGEFORMAT (2-2)所示：

                                       \* MERGEFORMAT (2-2)

式中，Sa(Ti)是第i阶周期对应的谱加速度；n是所考虑的结构的平动周期数，与结构的基本周期相关的待定参数；a表示新的地震动强度指标，其反映了前n阶周期对应的结构的谱加速度的几何平均值，对于一个给定的结构，该指标可以直接通过加速度反应谱得到，且与常用的Sa(T1)指标有较好的延续性。显然，当结构周期非常短时，即n=1时，该指标退化成一阶周期对应的谱加速度Sa(T1)。卢啸等[31]通过大量参数研究，给出了关键参数n与结构基本周期的关系如式\* MERGEFORMAT (2-3)所示。

                              \* MERGEFORMAT (2-3)

为了验证所提出地震动强度指标的优越性，以上海中心大厦（即模型A）和另一580m的超高层建筑结构（即模型B）的三维有限元模型为基础，以FEMA p695中推荐的22组远场地震动记录作为基本输入，进行弹塑性增量时程分析，逐步增大地震动强度，直至结构发生倒塌，找到引起结构发生倒塌的临界地震动强度，比较该指标与其他较常用的地震动强度指标所对应的临界倒塌强度IMcritical的变异性系数COV，如表2.3所示。变异性系数越小，表示该指标对超高层建筑结构越有效，越适用于超高层建筑的抗震设计和分析。

3.      本研究中所比较的地震动强度指标

22.    倒塌临界状态各地震动强度指标的变异系数

图2.21为两栋超高层建筑发生倒塌的临界地震动强度指标的变异系数。可见，在利用不同地震动强度指标表征模型A发生倒塌的临界值时，指标a对应的变异性系数最小，约为0.207，仅为Sa(T1)对应变异性系数的37.5%；PGV指标所对应的变异性系数次之，约为0.249，为Sa(T1)对应变异性系数的45.1%；Sa(T1)的变异性系数最大，说明仅考虑结构一阶模态是不能反映超高层结构地震响应的；IM1E&2E考虑了结构前两阶振动模态，其变异性系数比Sa(T1)有所减小，但仍大于a的变异性系数，说明仅考虑前两阶模态依然不足以反映超高层建筑结构的地震响应。同理，图2.21b可以看出，在利用不同的地震动强度指标表征模型B发生倒塌的临界值时，指标PGV对应变异性系数最小，指标a的变异性系数次之，略大于PGV的变异性系数；PGD和Sa(T1)的变异性系数最大。总的说来，PGV和本研究提出的指标a对超高层结构的都具有较好的适用性。但是，根据叶列平等[32]的研究，对于短周期结构，PGV指标精度较差，而Sa(T1)的精度最好，本研究提出的指标a在短周期情况下自动退化为Sa(T1)指标，因此比PGV指标有着更好的通用性。

另外，本研究提出的地震动强度指标是基于谱加速度的，能较好地利用规范的反应谱，也便于地震风险性分析中直接利用现有的地震动衰减关系，因此，本研究提出的地震动强度指标a也是适用于超高层建筑结构的抗震设计的地震动强度指标之一。

2.4.2   最小地震剪力系数对超高层建筑抗倒塌能力影响

目前在建筑结构抗震设计时，基本都是基于加速度反应谱确定地震作用力。长周期地震动记录中的速度或者位移特性对结构可能具有更大的破坏性，基于规范加速度反应谱的底部剪力法或振型分解法所确定的地震力，并不能很好地反映长周期地震动对结构的影响。对于超高层长周期结构，基于规范加速度反应谱确定的底部剪力会很小。出于对结构安全的考虑，国内外的抗震相关规范中都提出了最小地震剪力系数（或最小基底剪力）的要求。而近年来，我国大量的超高层建筑设计实例表明[33]对于基本周期很长的超高层建筑，按振型分解法计算的地震剪力系数并不满足我国现行规范的最小地震剪力系数的要求。显然最小地震剪力系数已经成为了超高层建筑抗震设计中的关键难题。而通过倒塌分析方法，则可定量地研究最小地震剪力系数对超高层建筑抗倒塌能力的影响，进而完善超高层建筑的抗震设计理论。

为讨论抗震规范中最小地震剪力系数对超高层建筑结构抗震性能及抗特大震倒塌能力的影响，以我国8度区的某超高层为例，按照表2.4中两种不同设计地震剪力系数取值方法，设计了2个模型，其设计地震剪力的取值方法具体介绍如下：

4.      两个方案设计地震剪力取值差异比较

（1）模型A：振型分解法计算得到的楼层剪力完全满足我国现行《建筑抗震设计规范》(GB 50011-2010) 和《高层建筑混凝土结构技术规程》(JGJ 3-2010)中最小地震剪力系数要求；

（2）模型B：振型分解法计算得到的楼层剪力不满足抗震规范最小地震剪力系数要求，但基底总剪力不小于最小基底剪力限值的80%，对于不满足规范的楼层，按照最小地震剪力系数调整得到的楼层剪力进行内力组合和构件设计。

24.    结构三维示意图

该超高层建筑位于我国抗震设防8度区II类场地第一分组，场地特征周期为Tg=0.35 s，主体结构总高度为439 m，共97层，顶部为4层观光层，地下4层，约18 m。该结构平面尺寸约为36.1 m×53.7 m，外围为16根矩形钢管混凝土巨柱，内部为21 m×37 m的矩形核心筒，在结构的长边方向采用了“巨型钢管混凝土柱－核心筒－桁架”的混合抗侧力体系；为保证结构短边方向的刚度，在结构短边方向沿结构全高布置了“X”型巨型钢斜撑抵抗水平荷载。5道腰桁架均匀分布在结构的高度方向，每道桁架体系高约8.3 m，结构整体示意图如图2.22所示。

模型A与模型B的结构布置完全一致。为使得模型A振型分解法计算的楼层剪力满足我国规范最小地震剪力系数的要求，必须尽可能增大结构刚度，减小结构周期。虽然采用了加大框架柱边长、增加剪力墙厚度、剪力墙内插钢板等一系列手段，但是仍无法满足周期要求。因此，从研究的角度出发，通过反复调整，最终模型A的核心筒采用C100钢纤维高强混凝土，使得结构周期可以降低到6 s以下。虽然超高强混凝土在实际工程中还较少采用，但已经有了较多的研究[34]。从表2.5的对比可以看出，毫无疑问，模型A的材料用量和建设难度要远大于模型B。

5.      模型A和模型B典型截面差异比较

模型A和模型B的基本动力特性比较如表2.6所示。显然，对于采用刚度调整方案的模型A，其基本自振周期明显小于模型B，结构整体抗侧刚度明显大于模型B；由于模型A外围钢管混凝土柱和剪力墙的尺寸要明显大于模型B，因此，模型A在增大结构抗侧刚度的同时，其重力荷载代表值也有了明显增加，是模型B重力荷载代表值的1.24倍，也就意味着模型A的材料用量要明显多于模型B。显然，要使得振型分解法计算的地震剪力完全满足规范最小地震剪力系数的要求，必须要付出很多的代价，且这样的设计在实际工程中实现难度也非常大，经济性也较差。

6.      模型A和模型B基本动力特性比较

为了进一步比较两种不同剪重比控制方案模型的抗倒塌能力差异，利用倒塌分析方法对上述两个模型（模型A和B）进行了典型地震波下的倒塌模拟。从FEAM p695推荐的22组远场地震动记录中选择Kocaeli_DZC180为基本地震动输入，其5%阻尼比的弹性反应谱与规范反应谱的比较如图2.23所示。倒塌分析方法为，沿模型x轴单向输入逐步增大的地震力，直至结构发生倒塌，从而得到引起结构倒塌的临界地震动强度，并以此来分析结构的抗倒塌能力。

对于模型A，当Kocaeli_DZC180地震动记录的PGA增大到1.4g时，结构发生倒塌；对于模型B，当Kocaeli_DZC180地震动记录的PGA增大到2.0g时，结构发生倒塌。在建筑结构抗倒塌研究中，一般采用抗倒塌安全储备CMR来量化反映结构的抗倒塌能力[35]，可按式\* MERGEFORMAT (2-4)计算。

                                             \* MERGEFORMAT (2-4)

其中，PGACollapse为引起结构倒塌的临界地震动强度，PGAMCE为结构大震所对应的地震动强度。从而可得到模型A和模型B在Kocaeli_DZC180地震动记录输入下抗倒塌安全储备系数(CMR)的比较如图2.24所示。模型A在Kocaeli_DZC180地震动记录输入的抗倒塌安全储备系数CMR约为3.5（1.4g/0.4g），而模型B在Kocaeli_DZC180地震动记录输入的抗倒塌安全储备系数CMR约为5.0（2.0g/0.4g）。可见，虽然模型A振型分解法计算的地震剪力完全满足抗震规范中最小地震剪力系数的要求，结构抗侧刚度、构件和结构整体的绝对承载能力均高于模型B，但是，由于结构抗侧刚度大同时也使得结构在地震作用下的水平地震剪力大幅增大，地震力需求明显提高，导致模型A的抗震能力－需求比并没有明显提高。反而模型B在振型分解法计算的地震剪力不满足规范最小地震剪力系数的要求时，通过直接增大楼层地震剪力的方法进行结构设计，在不改变结构地震力需求的同时增大了结构的抗震能力，使其能力－需求比明显提高，因此，在Kocaeli_DZC180地震动记录输入下，模型B体现出了更高的抗倒塌安全储备。故而在相关超高层建筑结构的抗震设计中，当最小地震剪力系数不能满足规范要求时，建议采用强度调整的办法，直接增大设计地震剪力的方法进行设计。Equation Chapter (Next) Section 1

3          城市区域的地震灾变模拟

近年来，全球范围内地震频发，从2008年到2012年，全球上7级以上大地震共发生了111起[36]，其中有两起大地震发生在我国内陆地区，分别是2008年的汶川地震和2010年的玉树地震。不仅如此，中、小型地震也频频发生。从2008年到2012年，我国共发生了5级以上地震206起[36]，小震更是不计其数。城市人口密度和财产密度都非常高，尤其是特大型城市，一旦发生强烈地震，必然会造成严重的人员伤亡和财产损失。中小地震引起的非结构构件破坏和室内设备损坏，也可能导致非常严重的损失。因此，科学、合理地进行城市区域震害预测是目前我们国家乃至世界范围内的一项迫切需求。

3.1   城市区域建筑震害模拟的精细化建模

3.1.1   精细化建模需求

在过去的30年，基于易损性矩阵的区域建筑震害预测方法得到广泛应用[37]，该方法通过给出不同场地烈度下某类型建筑达到不同破坏状态的概率（此概率一般可通过历史震害经验调查得到）来对震害损失进行预测。但是易损性矩阵方法只能给出宏观和统计意义上的预测结果，无法反映具体建筑物的详细损伤情况，也难以反映某一具体地震事件的特点。因此，基于能力－需求分析的建筑震害预测方法也得到了大量研究[38-40] 。例如，Hazus 99 的Advanced Engineering Building Modules (AEBM)[39]通过计算建筑弹塑性能力谱和地震需求谱的交点得到性能点，进而根据性能点来判断建筑物的损失情况。该方法相对易损性矩阵方法有了很大的进步，它能够计算出每栋建筑的位移和最大绝对加速需求，因此可以更加准确地对结构构件和非结构构件进行损失预测，也可以考虑不同地震动之间的差异。然而，该方法仍存在以下两方面的问题：(1) 将实际建筑物简化为单自由度体系，无法考虑高阶振型影响；(2) 采用基于静力推覆的能力－需求分析，无法考虑地震动的一些特性（例如地震动时域特性和速度脉冲）的影响[41]。

近年来，随着计算机数值模拟技术的不断发展，精细化模型和非线性时程分析越来越多地被应用于城市区域结构震害分析，但带来的问题是计算量过大。为了解决计算精度和计算效率之间的矛盾，本文作者[42]提出了采用GPU作为计算平台，以剪切层模型（如图3.1和图3.2）作为计算模型，采用非线性时程计算作为分析手段，进行高效城市区域震害分析的方法。剪切层模型有着比较适中的简化度以及相对较低的计算量[43]，相对于单自由度体系，采用集中质量剪切模型可以得到结构每一层的破坏状态，结构高阶振型以及地震波的速度脉冲的影响也能够较好地考虑。国内外许多研究都证明了采用该模型进行多层建筑的地震响应模拟的可靠性[44-47]。然而，如何确定层模型中的层间滞回行为对保障预测结果的可靠性至关重要。对于城市区域，想要精确得到其中每座建筑每一层的侧移刚度和滞回参数是十分困难的，多数情况下仅能得到一些结构的宏观信息。国内外学者在此方面做了大量的研究[48-52]，但是该问题尚未得到圆满解决。

本研究以Hazus软件作为基础，提出了相应的城市区域建筑物震害预测剪切层模型，并基于Hazus软件中通过大量调查统计得到的建筑抗震性能参数，提出了剪切层模型参数的确定方法，为精细化的城市区域建筑震害预测提供参考。

1.      单体建筑集中质量剪切模型

2.      群体建筑集中质量剪切模型

3.1.2   结构类型和设防等级

Hazus将城市区域中一般建筑按照结构类型和高度分为36个建筑类型[53]。本研究考虑我国工程实际情况，选择了其中19种建筑类型（见表3.1）。对于同一建筑类型，Hazus考虑到不同建造年代设计规范的差异，其抗震性能参数也有所不同。本研究采用了Lin等[54]提出的Hazus抗震设计等级与我国建筑的对照关系（表3.2），用来确定我国建筑与Hazus所提供参数的对应关系。

1.      本研究使用的建筑类型[53]

2.      Hazus抗震设计等级与我国建筑对应关系[54]

3.1.3   层间恢复力模型

Hazus建议结构的骨架曲线可采用三线性模型（图3.3），本研究层间滞回模型的骨架线也采用此模型。该模型可以分别描述结构层间的弹性阶段，屈服阶段和塑性阶段。整个骨架曲线需要5个参数来标定：k0（初始侧向刚度），Vy （层间剪切屈服强度），η （强化系数），β （极限强度和屈服强度比）和 Δc （结构的层间极限位移）。

4.      层间骨架线

不同结构在地震下往复受力行为有较大差异。根据结构类型的不同，本研究采用以下三种不同的层间滞回模型来描述不同类型的结构。

1. 修正的Clough模型

修正的Clough模型[55]（图3.4 a）被广泛应用于模拟弯曲破坏的钢筋混凝土框架结构，因此本研究采用该模型模拟钢筋混凝土框架结构的层间滞回关系（也就是C1类型）。

2. 理想弹塑性模型

理想弹塑性模型（图3.4b）可以用于模拟钢框架结构。因此本研究采用该模型模拟钢框架结构的层间滞回关系（也就是S1和S3类型）。

3. 捏拢模型

对于剪切破坏起主导作用的结构，采用捏拢模型（图3.4c）进行模拟是一个可行的方法。本研究采用Steelman和Hajjar[56]提出的捏拢模型。在该模型中，卸载刚度等于初始加载刚度，加载捏拢的转折点将出现在固定的转折线上，如图3.5所示。因此，本模型仅需1个描述结构退化程度的参数即可确定其滞回行为，如式\* MERGEFORMAT (3-1)所示：

                                        \* MERGEFORMAT (3-1)

其中 Ap和 Ab 分别是捏拢滞回包络面积和理想弹塑性滞回包络面积。选用该捏拢模型是因为其参数τ可以通过Hazus技术手册[53]提供的退化系数κ比较方便的确定[56]。如果采用其他多参数的捏拢模型(例如Ibarra捏拢模型[57])，参数标定的工作将非常复杂，不适用于区域震害预测的实现。在本研究中，除C1、S1、S3之外的所有建筑类型均采用此模型。

6.      捏拢模型滞回行为

3.1.4   一般建筑剪切层模型参数确定

在Hazus中，在每一个抗震等级下，每一种结构类型都有一个典型建筑，可通过该典型建筑的已知的参数来标定该类型结构的抗震性能。这些参数包括：①结构初始弹性一阶周期；②结构性能曲线；③结构振型质量系数；④结构高度、楼层数；⑤结构材料弹性阻尼；⑥结构退化系数；⑦结构各个破坏状态所对应的层间位移角限值。

在本研究中，对于区域中的某一个目标建筑，仅需知道其结构类型、层高、层数、建筑年代、面积，就能通过Hazus已知的参数来确定其集中质量剪切模型的参数，具体方法如下：

（1） 结构周期

对于区域中的一般多层建筑，其地震响应主要受前两阶振型控制[58]。对于城市区域中的一般建筑，前两阶周期通过式\* MERGEFORMAT (3-2)进行计算[59]：

其中T1和T2分别是目标建筑的一阶和二阶振动周期，N是目标建筑的层数。N0和T0分别是Hazus给出的典型建筑的层数和一阶周期。

（2）结构骨架线参数

对于一般多层建筑，其层质量和层初始刚度可认为是各层基本相同的[39]。因此，剪切层模型的刚度矩阵和质量矩阵可写为：

                    \* MERGEFORMAT (3-3)

                          \* MERGEFORMAT (3-4)

其一阶圆频率可由式\* MERGEFORMAT (3-5)计算[58]

                 \* MERGEFORMAT (3-5)

其中[Φ1]是结构的一阶振型向量。如果结构的楼层数确定，则无论k0和m如何取值，[Φ1]的各自由度比例是确定的，且其值可以用广义特征值法很简便地加以确定[58]。

于是结构的层间初始刚度可由式\* MERGEFORMAT (3-6)计算

                               \* MERGEFORMAT (3-6)

                                     \* MERGEFORMAT (3-7)

取值与层数的关系见表3.3（为计算方便，表中给出的数值为的数值）。

3.       取值与层数之间的关系（1至12层）

于是，目标建筑第i层的层间骨架线参数（图3.3）可由式\* MERGEFORMAT (3-8)确定：

                                   \* MERGEFORMAT (3-8) (a)

                          (3-8) (b)

                             (3-8) (c)

                                        (3-8) (d)

其中m是结构每层的质量，可以根据结构层面积和建筑用途确定[60]；g 是重力加速度；δCo是Hazus建议的结构完全破坏所对应的层间位移角限值[53]；h 是结构层高；(SDy, SAy)，(SDu , SAu) 是Hazus中给出的典型结构性能曲线的屈服点和极限点[53]；α1是Hazus给出的典型建筑的振型质量系数[57]；Γi是结构第i层的设计抗剪强度Vy,i与基底设计抗剪强度Vy,1的比值：

                                           \* MERGEFORMAT (3-9)

在我国规范中，设计地震力随结构高度基本呈倒三角形[61]。在本研究中，Γi由式\* MERGEFORMAT (3-10)确定:

                         \* MERGEFORMAT (3-10)

其中Wj，Wk是结构第j，k层的重量，Hj， Hk则是结构第j，k层所在平面距离地面的高程。

（3）结构滞回模型参数

对于修正的Clough模型和理想弹塑性模型，其行为仅需要骨架线就可以确定。对于本研究给出的捏拢模型，则需要一个参数τ来确定其行为。参数τ的取值可以参照Steelman和Hajjar[56]提出的方法，取所对应结构类型在Hazus中的退化参数κm加以确定[53]。

（4）结构阻尼

本研究中，结构阻尼采用瑞雷阻尼进行计算[58]，其阻尼比确定参照表3.4。

4.      结构阻尼比取值[53]

（5）结构破坏状态限值

对于城市区域中一般建筑，其结构破坏状态采用Hazus定义的4种破坏状态：轻微破坏，中等破坏，严重破坏，完全破坏。而确定结构破坏状态的指标则采用位移角进行确定，其对应的取值可在Hazus手册[53]表5.9内查到。

3.1.5   特殊建筑剪切层模型的参数确定

对于一些特殊建筑（如新型结构类型，混合结构，不规则结构），其抗震性能与一般多层建筑差别很大，Hazus软件对此类建筑也没有相应的规定。对于这类建筑，可采用更加精细的模型（如纤维梁模型或分层壳模型[11, 62-65]）来确定其参数，如图3.6所示。

8.      特殊建筑剪切层模型的参数确定流程

3.1.6   震害分析模型的数据获取

本研究中仅需获取结构类型、层高、层数、建筑年代、面积等5个参数，就能通过Hazus已知的参数来确定其建筑震害分析模型。然而，城市建筑数量庞大，每一栋建筑都需要采集相应的参数，这就导致城市区域建筑震害分析工作需要依赖详尽的城市基础数据库。当城市基础数据库的信息详尽程度难以满足震害分析需求时，如何高效、快速、准确地获取城市建筑的基础资料成为了一个重要问题。

近年来Google Earth的城市区域3D建筑模型在全球各大城市依次上线，如图3.7。与此同时，国内的一些商业公司也逐渐完成了国内多个城市的区域3D数字模型，如图3.8。城市区域建筑的3D数字虚拟化正在成为一种趋势。城市区域3D建筑模型提供了精细的城市建筑尺寸信息。充分利用该数据，为城市区域建筑震害的模拟提供了一种新的数据获取思路。

为此，本研究提出了一种基于Google Earth城市区域3D建筑模型的建筑参数获取方法。该方法可以通过Google的Dae三维数字模型格式，获得分析区域内每一栋建筑的三维模型参数，如层高、楼面形状等，进而可以获得层高、层数、楼面面积等震害分析所需数据。

（1）方法介绍

在Google Earth城市区域3D建筑模型进行城市区域3D建筑几何参数的抓取，面临如下问题：

1. 城市三维数字模型本质上是由多边形或者三角形组成的。Dae模型文件中仅包含所有多边形的端点坐标信息以及连接关系信息。事先通常不知道某栋建筑是由哪些多边形构成的，如图3.9所示。
2. 城市三维数字模型中除了建筑模型，通常还包括地表，道路，植被等模型。因此需要剔除掉其他的模型，单独识别出构成建筑的多边形，如图3.9所示。
3. 为了实现建筑的地震响应模拟，在获得到了组成某栋建筑的所有多边形之后，需要根据这些多边形数据，抽取出结构外形的参数。包括结构高度，结构各层的平面布置情况。

11.    城市区域3D建筑模型示意

面对以上问题，本研究采用了基于城市切片的建筑参数抽取方法，具体实现步骤如下所示：

1. 对整个城市区域3D模型中的所有多边形进行切片。对于每个多边形，得到对应该高程的一条直线。
2. 提取出能围成闭合多边形的直线，将该多边形作为某栋建筑的层外形多边形，如图3.10所示。重复步骤2，得到各个高程的建筑层外形多边形，如图3.11所示。
3. 判断各个多边形的竖向重叠关系，对于竖向重叠的多边形认为这些多边形组成一栋建筑，这样就获取到了区域中建筑的各层外形数据。
4. 对于已经识别的建筑，计算其多边形的高程差，从而估算得到建筑的高度。

12.    对模型的所有多边形进行切片

13.    得到各个高程的建筑层外形多边形

（2）实例验证

通过上述的步骤，就能基本获得区域中建筑的各层外形参数以及各建筑的高度参数。为了校验本方法的可用性，首先可以对单栋建筑进行识别验证，如图3.12所示。图中半透明的模型为原建筑3D模型。橙色的部分为识别得到的各层建筑外形多边形。

此外对于整个区域，为了校验识别的效果，也可以与原模型进行对比，如图3.13所示。图中黄色部分为识别得到的建筑层外形多边形，可以看出虽然原模型中包含地形等无关的多边形数据。但是该方法依然能较好的识别出每一栋建筑。

通过识别出的建筑以及楼层多边形，可以快速计算楼层高度、楼层数量、楼层面积等建筑外形参数，为城市震害分析提供了重要的数据支持。

14.    单栋建筑层外形多边形识别效果

15.    区域建筑层外形多边形识别效果

3.2   城市区域震害模拟的高性能计算

3.2.1   基于GPU粗粒度并行的建筑震害计算

虽然基于精细化模型的非线性时程分析在单体建筑中已经得到了广泛的应用[66]，但是城市区域拥有海量的建筑，该方法所需的计算量巨大。2008年Hori和Ichimura[67]基于精细化结构模型和非线性时程分析开发了Integrated Earthquake Simulation (IES) 平台，实现了整个东京的震害预测。但是该平台必须依赖超级计算机，这意味着高昂的运营费用和维护成本。因此，为了更好地开展城市区域建筑震害的精细化模拟，迫切需要一种低成本、高性能的计算平台。

最近几年，图形处理单元(GPU)技术飞速发展并在通用计算领域得到广泛应用。虽然GPU单核的计算能力相对CPU较弱，但是GPU的处理核心数远多于相同价格的CPU，因此相同价格的GPU相对CPU具有更高的计算性能[68]。随着NVIDIA公司发布了CUDA (Compute Unified Device Architecture )，GPU通用计算的编程难度被大大降低。目前GPU计算在生物、电磁场、地理等领域得到了大量的应用[69-71]。

最能发挥GPU性能的是细粒度的并行计算[72]，即将每个子任务划分成许多更细小的操作步，然后GPU在操作步层面对任务进行并行计算。这种并行方式在神经网络以及有限元领域得到了很广泛的应用[73, 74]。但是这种并行计算方式需要对GPU程序进行细致的调试以实现对大量的操作步高效的运算，程序实现难度较大。而相对于细粒度的并行，粗粒度的并行方式则更容易实现。粗粒度的并行计算方式采用的是基于相对独立的子任务的并行而不是操作步的并行。这种并行方式在优化和控制领域都有应用[75]。当满足以下情况时，粗粒度和细粒度的并行效率相当：

(1) 子任务的数量远多于GPU计算的核心数；

(2) 每个子任务的计算量比较适中，能够独立在一个GPU核心上完成；

(3) 各个子任务之间不需要太多的数据交换。并且计算过程中不需要全局同步，每个GPU核心上的子任务可以按顺序逐个的进行计算。

幸运的是区域城市建筑震害预测分析可以满足以上三个特性。虽然整个城市有成千上万栋的建筑，但是如果通过采用适当的计算模型使每栋建筑结构的计算量不超过GPU单核的计算能力，则每栋单体建筑的计算就能被看作一个相对独立的计算子任务。此外不同建筑的地震响应相对独立，不同GPU线程之间几乎不需要进行数据交换。因此GPU数百计的计算核心可以实现几百栋建筑同时计算，这样计算效率将非常高。

本研究提出适合多自由度剪切层模型的GPU/CPU协同粗粒度并行计算的区域建筑震害预测方法，对并行计算的效率进行了详细的讨论，并对一个中等大小的城市进行了精细化的震害模拟。

3.2.2   CPU/GPU协同的程序构架

整个程序包括3个模块：前处理模块、结构分析模块和后处理模块。前处理模块的主要任务是获取计算模型的主要参数，并且选择地震场景进行非线性时程分析。后处理模块主要展示模拟结果并预测损失。

计算分析模块是程序的核心，由CPU和GPU协同完成计算任务。其中，CPU负责文件读写、任务分配等工作；GPU开启大量线程，每一条线程负责完成每一栋单体建筑的结构非线性时程计算，如图3.14所示。

本研究采用多自由度剪切层模型来模拟每个单体建筑物[42]，它可以很好地平衡计算量和计算精度之间的关系。为了避免隐式动力计算的收敛性问题，并提高GPU并行计算效率，本研究采用中心差分法求解动力方程[58]。动力方程中的阻尼采用经典瑞利阻尼。

3.2.3   GPU计算性能分析

为了对比本研究提出的基于GPU/CPU协同粗粒度并行计算的结构分析模块的性能，在传统CPU平台上也开发了一个功能相同计算程序。本研究采用的对比算例和对比平台如下所示：

(1) 根据Hazus中的建筑类型，随机生成不同结构类型的1024栋建筑，如表3.5所示。

(2) 采用广泛使用的El Centro 地震动时程记录[76]进行对比分析。地震动时程记录的峰值加速按照200 cm/s2进行调幅。

(3) 计算时长为40s，时间步为8000步。

(4) 为了避免读写速度的影响，下文中的计算时间不包括硬盘的读写时间。

计算平台如表3.6所示。所采用的两套计算平台均为2011年采购，在当时两者价格基本相同。因此采用以上两套平台可以用于比较相同价格情况下的性能差异。

5.      性能对比分析中使用的建筑结构类型分布情况

6.      对比采用的计算平台

首先采用单栋建筑对比CPU平台和CPU/GPU协同计算平台的计算效率。楼层数和单栋建筑平均计算时间的关系如图3.15所示。从图中可以看出对于一栋建筑GPU/CPU协同计算的时间要长于CPU计算的时间，即单颗GPU核心的计算能力要明显弱于CPU的计算能力。此外还可以注意到对于一栋10层楼建筑，采用单精度计算GPU/CPU协同计算平台需要5秒，如果采用双精度计算则需要8秒。然而CPU平台的单精度和双精度计算时间非常接近（双精度还稍快），这是因为使用的CPU是基于64位构架，其默认计算精度就是双精度。

对于1024栋建筑物，如果改变CUDA的block size（每个计算块的线程数量），得到的计算时间如图3.16所示，可见当程序的block size为32时程序能达到最高的计算能力。这是因为，当采用CUDA进行并行计算的时候每32个线程组成一个wrap，而GPU计算任务的最小单位是一个wrap[77]。因此，如果block size小于32时，GPU的性能不能充分发挥。而当单精度的block size大于520或者双精度的block size 大于260时，将超出GPU中速度最快的寄存器的容量（本平台GPU每个block只能使用32K的32位寄存器，可存储32,768个单精度数或16,384个双精度数），因此也会导致并行效率下降。

在1024栋建筑中，随机选取不同的建筑数量在CPU平台和GPU/CPU平台上进行性能对比（GPU/CPU平台的block size是32），结果如图3.17所示。CPU的计算时间与建筑的数量基本上是线性的关系。然而GPU/CPU协同计算平台的计算时间则主要取决于计算时间最长的那个线程（图3.17）。此外GPU/CPU协同计算平台双精度计算的曲线上有几个曲线的突跃，这是由于GTX460的双精度计算是在特殊函数单元中完成的 (Special Function Unit)，其数量是CUDA核心数的1/6。这个问题可望在Fermi架构的Tesla GPU上得到修正，因为其双精度计算是直接在CUDA核心上完成的[77]。

由GPU/CPU平台与CPU平台的计算时间对比可以看出GPU对繁重的非线性时程计算具有巨大的速度优势。如图3.18所示，当对1024栋建筑进行单精度计算时，采用GPU/CPU协同计算的时间仅是CPU计算的1/39，如果采用双精度则是1/21。此外，本文作者[42]对比了GPU单精度和双精度计算结果，发现差异不足0.1%，因此对于区域建筑震害模拟而言，单精度计算已经可以满足要求。

3.2.4   中等城市的建筑震害模拟

以某中等城市为例应用本研究中提出的方法对该城市进行了建筑震害模拟，并讨论本研究方法的优势。该中等城市总共有4255栋建筑。将该地区的建筑按照表3.7对不同年代的建筑进行了抗震等级分类[54]。对不同的抗震设计等级将根据Hazus中提出的方法确定建筑结构的层间滞回参数[53]。不同结构类型建筑的统计信息如表3.8所示。

7.      结构抗震标准的划分

8.      该区域建筑的统计情况

分析采用了3类地震动数据：远场地震动，近场无速度脉冲地震动，近场有速度脉冲地震动。对于每一类地震动，从FEMA p695推荐的地震动数据库中分别选取5条（表3.9）[78]。根据该城市的设防烈度，对于小震、中震、大震分别采用70, 200和 400 cm/s2的地面峰值加速度对选取的地震动记录进行调幅。由于本研究主要讨论建筑震害模拟方法，因此没有考虑场地对地震动的影响，每栋建筑都采用相同的地震动记录。

9.      从 PEER-NGA 数据库中选择的地震动时程记录

为了体现本研究采用的多自由度剪切层模型的优势，将分别使用本研究的多自由度模型方法和传统的单自由度模型方法对该城市进行模拟，单自由度模型的参数根据Steelman和Hajjor建议的方法[56]确定。对于这样一个中等城市，多自由度剪切层模型和单自由度模型完成一次地震模拟的计算时间完全可以满足应急震害预测的速度要求。

计算得到的最大楼层损伤状态的分析结果如图3.19所示。图中结果显示采用多自由度剪切层模型得到的楼层损伤结果普遍比单自由度模型的预测结果大一些。这是由于多自由度剪切层模型可以得到各个楼层的响应，从而能够考虑损伤集中现象。对于200gal远场地震波记录和近场无速度脉冲地震动记录，单自由度模型和多自由度剪切层模型的结果差异较小。但是对于近场有脉冲地震动记录，单自由度模型和多自由度剪切层模型的结果具有较大的差别。这是由于速度脉冲能够激励出高阶阵型从而导致结构更加严重的破坏，而单自由度模型无法考虑高阶阵型的影响。

以峰值加速度为200 cm/s2 的IMPVALL/H- BCR_233地震动为例，建筑破坏的预测结果如图3.20、图3.21所示。结果显示采用多自由度剪切层模型，能够准确获得建筑损伤的楼层位置和最大加速度位置，显著优于单自由度模型的计算结果。

本研究的结果可以证明GPU/CPU协同粗粒度并行计算具有很高的计算效率，其性能价格比能够达到传统CPU计算的39倍。并且，采用的多自由度剪切层模型可以考虑地震动中速度脉冲的影响，也能确定不同楼层的损伤情况，从而能够更加准确的进行损失预测。

3.3   城市区域震害的高真实度展示

3.3.1   震害真实感展示的需求

城市区域震害模拟的高真实度展示对防灾规划、虚拟应急训练、成果交流等都具有重要意义。在这方面，国内外很多研究者已经开展了基于三维城市模型的震害展示研究，如加拿大约克大学[79]基于分布式GIS建立了温哥华的三维地形和城市模型，并对城市整体震害进行评估；土耳其中东技术大学[80]、清华大学[81]利用CAD和纹理技术建立了真实感的3D城市，通过静态方式展现地震易损性或震害效果。但是这些研究所采用的震害预测模型往往过于简单，不能得到建筑物破坏的细节信息，因而其震害展示效果真实感不够。而近年来得到迅速发展的城市建筑群精细预测模型[60, 67]，可以得到更加清晰准确的震害计算结果（例如整个结构的震害时程以及楼层的局部损伤等），为提高震害模拟的真实感提供了关键的数据支持。

但是，精细化的区域建筑震害预测模型一般也难以包含建筑倒塌过程、残骸运动等倒塌细节。缺乏倒塌细节不仅影响震害模拟的真实感，也无法评价建筑倒塌对交通、救援等的影响。物理引擎[82]是近些年计算机图形学的新技术，专门用于计算场景中物体的复杂物理行为，可以模拟震害场景中建筑物的开裂、破碎、碰撞、堆积等一系列破坏过程，在科学仿真、电影制作上有相关应用[83-85]。应用物理引擎技术，可以弥补数值模拟缺失的细节，如清华大学[86]曾用物理引擎技术模拟桥梁垮塌全过程中残骸的运动，弥补有限元“生死单元”技术的不足。因此，物理引擎技术也可以应用在城市震害模拟中，为城市工程震害的动态破坏细节过程提供了技术支持。然而，国内外尚缺乏基于物理引擎等技术的城市工程震害视景仿真方面的深入研究，因此，在震灾真实感展示中，垮塌过程的真实感特效问题面临较大挑战。

本研究基于城市区域建筑震害精细化数值分析结果，提出真实感的城市区域建筑三维建模方法以及震害过程的动画实现方法，以准确而又真实感的表现震害分析结果；利用物理引擎，提出建筑震害过程中的垮塌效果的实现算法，以弥补数值模拟缺乏的细节，增加模拟完整性和真实感。并且，对我国一个中等规模城市的建筑震害过程进行综合模拟，为城市防灾减灾策略提供参考意见。

3.3.2   真实感建模与震害动态演示

（1）城市建筑群的真实感建模

城市建筑群的真实感建模主要包括三维建模和纹理映射两个部分，如图3.22所示。区域建筑震害预测模型将楼层简化为集中质量点，并不包含完整的三维信息。因此，需要利用城市区域的地理信息系统GIS数据，获取位置、建筑外形、楼层高度等扩充信息来建立建筑三维模型。本研究采用徐峰等人的基于GIS三维城市建筑建模方法[81]，利用GIS中建筑轮廓多边形在高度方向进行拉伸来建立不同的建筑三维模型。同时，为了使建筑三维模型能够反映层间的局部震害特征，建筑三维模型应与震害预测模型保持一致，同样以楼层为基本单元。

为增加建筑模型的真实感，需要对三维模型进行纹理映射。本研究基于Tsai和Lin的方法[87]，采用多重纹理技术，对建筑物的不同表面分别进行贴图，以更加细致地体现建筑物的特征。根据图3.22所示方法，可以建立城市区域的真实感模型，为震害动态演示奠定基础。

（2）基于精细化预测的震害动态演示

如何准确、高效地表现大规模的震害预测的时变数据是本研究城市区域建筑震害动态演示的主要问题。为了对图形表现过程进行充分控制，本研究选择开源图形引擎OSG（Open Scene Graph）作为视景模拟的主要平台[88]。在OSG中，一般会采用回调机制(Callback)来完成每一帧渲染前需要的工作，如空间变化、顶点更新、视点变化等。利用OSG的回调机制，本研究直接用震害预测数据在每一帧更新建筑模型的顶点坐标，以动态地表现建筑震害过程，具体算法如图3.23所示。

图3.23所示算法包含三个重要步骤：

1）更新器的创建与初始化

为了实现对图形顶点的操作，需要创建更新器updater。它需要继承OSG中图形的更新回调类UpdateCallback，以实现对回调过程的全面控制。在updater初始化阶段，更新器需要在动画渲染前加载对应建筑震害的节点位移数据。

2）更新器的动态执行

在回调过程中，需要将震害预测中的时间步与动态演示中的帧数相对应。在动态演示的每一帧中，更新器将利用对应时间步的震害数据来更新建筑的顶点坐标。重载更新器中的Update()函数可以实现对图形顶点的动态修改功能。通过回调，不断修改建筑模型，可形成区域建筑震害的动态演示过程。

3）清空更新器

当动态演示达到最大的时间步时，清空更新器，建筑模型不再更新，震害过程的动态演示结束。

3.3.3   基于物理引擎的倒塌模拟

区域震害预测模型只能给出建筑倒塌的判断，缺乏细节的倒塌过程。细节的倒塌过程可以增强震害模拟的真实感，也为地震疏散、救援调度等提供重要的参考，因此，倒塌模拟在地震模拟中是不可缺少的。而且，倒塌过程需要符合物理规律，并不是纯粹的图形表现。

为匹配上述需求，本研究选择物理引擎模拟建筑倒塌过程。由于在多刚体动力学和碰撞检测上的优势，物理引擎善于模拟倒塌过程中大量残骸的复杂运动和相互碰撞过程。本研究对比了三大物理引擎Havok，PhysX和Bullet[89-91]，发现由于PhysX可以采用GPU来加速计算过程，更适合大量物体的计算，而在区域震害模拟中，建筑的规模是庞大的，因此，PhysX更适合用于区域建筑的倒塌过程模拟。

为了保证物理计算与图形显示的准确对应，需要根据震害预测数据在物理空间中建立与图形一致的建筑倒塌模型。在本研究建筑模型场景层次中，楼层图形由OSG中的Geode节点储存和管理[88]。在PhysX中，物理角色Actor类是基本的计算单元，代表着物理世界中所有的物体[91]。因此，对应于Geode节点，在物理引擎中楼层由Actor建立。物理引擎的楼层模型，需要楼层的材料属性（如弹性系数、线阻尼等），也需要楼层倒塌时的状态参数（如位移、速度等）。材料属性需要根据不同建筑结构类型，查询相关实验结果确定，而状态参数通过数值分析数据获得，以保证模拟的准确性。在物理引擎中，楼层的外形参数通过对应的楼层的Geode节点获得。为减少震害模拟过程中的耗时，楼层的物理模型可以在渲染前建立，在需要物理计算时被激活。

物理引擎与图形引擎是两个独立的程序，需要构建两者的动态联系才能模拟建筑的倒塌过程。目前，关于OSG与PhysX的结合工作研究还非常有限。本研究中基于OSG回调机制，在每一帧中，利用PhysX实时计算数据不断更新楼层的顶点坐标，以形成楼层倒塌过程的动态演示[88]。为了建立OSG与PhysX的对应关系，创建了一个专门的类Stories。在Stories类中，同时包含楼层对应的Actor对象和Geode节点，以建立楼层物理模型与图形模型的联系。为了实现建筑倒塌的动态过程，需要控制Stories类的回调过程，因此，设计了一个新类Move-callback。由于继承了Callback机制的控制类，Move-callback类可以控制Stories类的渲染主循环。Move-callback类从PhysX的Actor中获取运动信息，然后应用运动信息来更新楼层图形。通过Stories类和Move-callback类实现的渲染循环过程形成了PhysX计算支持的建筑倒塌模拟，如图3.24所示。

建立OSG与PhysX的动态联系后，可以通过物理引擎模拟建筑倒塌过程。但是，为表现完整的建筑震害过程，震害动态演示与倒塌过程必须形成协同模拟的整体。根据震害预测给出的倒塌判断，将整个建筑震害过程分为两个阶段：倒塌前和倒塌后。首先，建立区域建筑群的真实感模型。在倒塌前，利用位移时程数据来动态地演示建筑的位移和变形。发生垮塌后，清空动态演示的更新器，根据前述方法建立楼层图形Geode节点与物理角色Actor的动态联系。然后，将Geode的运动状态信息传递给Actor，并激活Actor。物理引擎将不断计算Actor的倒塌过程，计算数据也会同步传递给Geode，以实时地楼表现建筑的倒塌过程。上述过程（如图3.25所示）实现了完整的城市区域建筑震害高真实度模拟。

3.3.4   中等城市的震害真实感展示

本研究以前文中的中型城市震害模拟作为算例，基于精细化震害预测结果和GIS数据建立了真实感的三维城市模型（图3.26a），并展示了震害过程中不同楼层间的位移和变形（图3.26b）。由于城市区域的建筑数量众多，所有建筑倒塌过程的非线性时程分析是不可能通过普通计算机模拟的。然而，通过使用PhysX，可以在桌面型计算机中快速模拟建筑倒塌过程（图3.26c 和 图3.26d）。倒塌过程模拟很好地弥补了非线性时程分析的不足，使整个震害模拟过程更加真实、完整。

为进一步增加震害场景的完整性和真实感，本研究将震害场景加入了地形、天空等环境模型，并且利用立体投影设备进行了区域建筑震害过程的立体感展示，进一步加强了震害过程模拟的沉浸感，如图3.27所示。图3.27的结果可以用于虚拟环境下地震疏散、应急调度等应用，综合提升防震减灾能力。

本研究基于精细化的结构模型和物理引擎技术，实现了城市区域建筑震害的真实感显示，对我国一个中型城市进行了区域建筑震害模拟。本研究不仅准确、真实感地表现了区域建筑震害过程，同时，也给出符合物理规律的倒塌过程，弥补了震害模拟的细节缺失。研究结果可以用于震害预测、灾后应急、救援训练等，为城市防震减灾提供了重要技术支持。

4          结论

Reitherman指出，地震工程学所面临的三个最大困难在于“风险、非弹性和动力学（Risk, inelasticity and dynamics）”[92]。虽然地震工程很多重要思想，在20世纪中期、甚至19世纪末期，就已经提出来了。但是离开了计算机的帮助，定量且精确的研究上述三个问题中的任何一个都是非常困难的。而自从计算机出现后，结构工程和地震工程的研究取得了翻天覆地的巨大变化。Roesset and Yao评价：“20世纪及未来若干年结构工程学最主要的改变是因为数字计算机的发展，包括将计算机作为一个强大的工具进行繁琐的计算以及作为一个新的通讯工具供设计团队、教授和学生、及其他人员开展交流”[93]。充分利用计算机科学和技术方面取得的最新成果，加深对工程结构和城市在强烈地震下的灾变演化规律和机理的理解，进而提出高效可靠的工程对策，对工程防灾减灾有着重大而深远的意义。

致谢

感谢国家自然科学基金（No. 51222804，51178249，51378299，51308321，90815025，51261120377），北京市自然科学基金(No. 8142024)，国家科技支撑计划课题（No. 2013BAJ08B02），交通运输部交通运输建设重大科技专项（No. 2011-318-223-170），霍英东教育基金(No. 131071)和教育部新世纪优秀人才支持计划（No. NCET-10-0528）等项目的资助。并感谢清华大学钱稼茹教授、任爱珠教授，澳大利亚Griffith大学H Guan博士，日本东京大学M Hori教授，美国斯坦福大学KH Law教授，美国加州大学伯克利分校F McKenna博士，ARUP工程咨询公司YL Huang博士等对本文研究的贡献。参与本文研究的还有清华大学博士研究生汪训流、解琳琳、曾翔，硕士研究生唐代远、张万开等。

参考文献：

[1]              Yamashita T, Kajiwara K, Hori M. Petascale Computation for Earthquake Engineering. Computing in Science & Engineering, 2011,13 (4):44-49.

[2]              Ibarra LF, Krawinkler H. Global Collapse of Frame Structures Under Seismic Excitations. Pacific Earthquake Engineering Research Center, 2005.

[3]              Kabeyasawa T, Shiohara H, Otani S, Aoyama H. Analysis of the Full-Scale Seven-Story Reinforced Concrete Test Structure: Proc. 3rd JTCC, US-Japan Cooperative Earthquake Research Program, BRI, Tsukuba, Japan, 1982.

[4]              李国强, 周向明, 丁翔. 钢筋混凝土剪力墙非线性动力分析模型. 世界地震工程, 2000(02):13-18.

[5]              Taucer FF, Spacone E, Filippou FC. A Fiber Beam-Column Element for Seismic Response Analysis of Reinforced Concrete Structures. Report No. UCB/EERC-91/17, UC Berkeley, America, 1991.

[6]              Spacone E, Filippou FC, Taucer FF. Fibre Beam-Column Model for Non-Linear Analysis of R/C Frames .1. Formulation. Earthquake Engineering & Structural Dynamics, 1996,25 (7):711-725.

[7]              Legeron F, Paultre P. Uniaxial Confinement Model for Normal- and High-Strength Concrete Columns. Journal of Structural Engineering-ASCE, 2003,129 (2):241-252.

[8]              汪训流, 陆新征, 叶列平. 往复荷载下钢筋混凝土柱受力性能的数值模拟. 工程力学, 2007(12):76-81.

[9]              Esmaeily A, Xiao Y. Behavior of Reinforced Concrete Columns Under Variable Axial Loads: Analysis. ACI Structural Journal, 2005,102 (5):736-744.

[10]              叶列平, 陆新征, 马千里, 汪训流, 缪志伟. 混凝土结构抗震非线性分析模型、方法及算例. 工程力学, 2006(S2):131-140.

[11]              Li Y, Lu X, Guan H, Ye L. An Improved Tie Force Method for Progressive Collapse Resistance Design of Reinforced Concrete Frame Structures. Engineering Structures, 2011,33 (10):2931-2942.

[12]              Lu X, Lu X, Guan H, Ye L. Collapse Simulation of Reinforced Concrete High-Rise Building Induced by Extreme Earthquakes. Earthquake Engineering & Structural Dynamics, 2013,42 (5):705-723.

[13]              Ren P, Li Y, Guan H, Lu X. Progressive Collapse Resistance of Two Typical High-Rise RC Frame Shear Wall Structures. Journal of Performance of Constructed Facilities-ASCE, 2014.

[14]              唐代远. 等跨 RC 框架结构抗地震倒塌性能的试验研究及理论分析. 北京: 清华大学硕士学位论文, 2011.

[15]              门俊, 陆新征, 宋二祥, 陈肇元. 分层壳模型在剪力墙结构计算中的应用. 防护工程, 2006,28(3):9-13.

[16]              魏勇, 钱稼茹, 赵作周, 蔡益燕, 郁银泉, 申林. 高轴压比钢骨混凝土矮墙水平加载试验. 工业建筑, 2007(06):76-79.

[17]              解琳琳, 黄羽立, 陆新征, 林楷奇, 叶列平. 基于 OpenSees 的 RC 框架-核心筒超高层建筑抗震弹塑性分析. 工程力学, 2014(01):64-71.

[18]              Loland KE. Continuous Damage Model for Load-Response Estimation of Concrete. Cement and Concrete Research, 1980,10 (3):395-402.

[19]              Mazars J. A Description of Microscale and Macroscale Damage of Concrete Structures. Engineering Fracture Mechanics, 1986,25 (5-6):729-737.

[20]              Humphrey JR, Price DK, Spagnoli KE, Paolini AL, Kelmelis EJ. CULA: Hybrid GPU Accelerated Linear Algebra Routines, 1000 20TH ST, PO BOX 10, BELLINGHAM, WA 98227-0010 USA, 2010.

[21]              CuSP Home Page. http://cusplibrary.github.io/.

[22]              蒋欢军, 和留生, 吕西林, 丁洁民, 赵昕. 上海中心大厦抗震性能分析和振动台试验研究. 建筑结构学报, 2011(11):55-63.

[23]              Lu X, Lu X, Zhang W, Ye L. Collapse Simulation of a Super High-Rise Building Subjected to Extremely Strong Earthquakes. Science China-Technological Sciences, 2011,54 (10):2549-2560.

[24]              范立础, 胡世德, 叶爱君. 大跨度桥梁抗震设计. 人民交通出版社, 2001.

[25]              陈幼平, 周宏业. 斜拉桥地震破坏的计算研究. 地震工程与工程振动, 1995(03):127-134.

[26]              Yeh Y, Mo YL, Yang CY. Seismic Performance of Rectangular Hollow Bridge Columns. Journal of Structural Engineering-ASCE, 2002,128 (1):60-68.

[27]              Shome N, Cornell CA, Bazzurro P, Carballo JE. Earthquakes, Records, and Nonlinear Responses. Earthquake Spectra, 1998,14 (3):469-500.

[28]              Tothong P, Luco N. Probabilistic Seismic Demand Analysis Using Advanced Ground Motion Intensity Measures. Earthquake Engineering & Structural Dynamics, 2007,36 (13):1837-1860.

[29]              Luco N, Cornell CA. Structure-Specific Scalar Intensity Measures for Near-Source and Ordinary Earthquake Ground Motions. Earthquake Spectra, 2007,23 (2):357-392.

[30]              Baker JW, Cornell CA. A Vector-Valued Ground Motion Intensity Measure Consisting of Spectral Acceleration and Epsilon. Earthquake Engineering & Structural Dynamics, 2005,34 (10):1193-1217.

[31]              Lu X, Ye L, Lu X, Li M, Ma X. An Improved Ground Motion Intensity Measure for Super High-Rise Buildings. Science China-Technological Sciences, 2013,56 (6):1525-1533.

[32]              Ye L, Ma Q, Miao Z, Guan H, Yan Z. Numerical and Comparative Study of Earthquake Intensity Indices in Seismic Analysis. Structural Design of Tall and Special Buildings, 2013,22 (4):362-381.

[33]              汪大绥, 周建龙, 姜文伟, 王建, 江晓峰. 超高层结构地震剪力系数限值研究. 建筑结构, 2012(05):24-27.

[34]              李俊, 王震宇. 超高强混凝土单轴受压应力-应变关系的试验研究. 混凝土, 2008(10):11-14.

[35]              陆新征, 叶列平. 基于IDA分析的结构抗地震倒塌能力研究. 工程抗震与加固改造, 2010(01):13-18.

[36]              中国地震局. 历史地震目录. http://www.cea.gov.cn/publish/dizhenj/468/496/index.html.

[37]              Onur T, Ventura CE, Finn WDL. A Comparison of Two Regional Seismic Damage Estimation Methodologies. Canadian Journal of Civil Engineering, 2006,33 (11):1401-1409.

[38]              Earthquake Loss Estimation Methodology - HAZUS97. Technical Manual. Washington, D.C.: Federal Emergency Management Agency – National Institute of Building Sciences, 1997.

[39]              Earthquake Loss Estimation Methodology - HAZUS99. Technical Manual. Washington, D.C.: Federal Emergency Management Agency – National Institute of Building Sciences, 1999.

[40]              Sextos AG, Kappos AJ, Stylianidis KC. Computer-Aided Pre- and Post-Earthquake Assessment of Buildings Involving Database Compilation, GIS Visualization, and Mobile Data Transmission. Computer-Aided Civil and Infrastructure Engineering, 2008,23 (1):59-73.

[41]              Iwan WD. Drift Spectrum: Measure of Demand for Earthquake Ground Motions. Journal of Structural Engineering-ASCE, 1997,123 (4):397-404.

[42]              Lu X, Han B, Hori M, Xiong C, Xu Z. A Coarse-Grained Parallel Approach for Seismic Damage Simulations of Urban Areas Based On Refined Models and GPU/CPU Cooperative Computing. Advances in Engineering Software, 2014,70 (0):90-103.

[43]              Moghaddam H, Mohammadi RK. Ductility Reduction Factor of MDOF Shear-Building Structures. Journal of Earthquake Engineering, 2001,5 (3):425-440.

[44]              Veletsos AS, Vann WP. Response of Ground–Excited Elasto-Plastic Systems. Journal of Structural Division-ASCE, 1971,97 (4):1257-1281.

[45]              Diaz O, Mendoza E, Esteva L. Seismic Ductility Demands Predicted by Alternate Models of Building Frames. Earthquake Spectra, 1994,10 (3):465-487.

[46]              Pampanin S, Priestley M, Sritharan S. Analytical Modelling of the Seismic Behaviour of Precast Concrete Frames Designed with Ductile Connections. Journal of Earthquake Engineering, 2001,5 (3):329-367.

[47]              经杰, 叶列平, 钱稼茹. 基于能量概念的剪切型多自由度体系弹塑性地震位移反应分析. 工程力学, 2003(03):31-37.

[48]              Aoyama H. A Method for the Evaluation of the Seismic Capacity of Existing Reinforced Concrete Buildings in Japan. Bulletin of the New Zealand National Society for Earthquake Engineering, 1981,14 (3):105-130.

[49]              Hajirasouliha I, Doostan A. A Simplified Model for Seismic Response Prediction of Concentrically Braced Frames. Advances in Engineering Software, 2010,41 (3):497-505.

[50]              Muto K. Dynamic Design for Structures. Tokyo: Maruzen Publishing, 1977.

[51]              何广乾, 魏琏, 戴国莹. 论地震作用下多层剪切型结构的弹塑性变形计算. 土木工程学报, 1982(03):10-19.

[52]              陈光华. 地震作用下多层剪切型结构弹塑性位移反应的简化计算. 建筑结构学报, 1984(02):45-57.

[53]              FEMA. Multi-Hazard Loss Estimation Methodology Hazus –MH 2.1 Advanced Engineering Building Module (AEBM) Technical and User’S Manual. Washington, D.C.: Federal Emergency Management Agency, 2012.

[54]              Lin S, Xie L, Gong M, Li M. Performance-Based Methodology for Assessing Seismic Vulnerability and Capacity of Buildings. Earthquake Engineering and Engineering Vibration, 2010,9 (2):157-165.

[55]              Mahin SA, Lin J. Construction of Inelastic Response Spectra for Single-Degree-Of-Freedom Systems, Computer Program and Applications, Report No. UCB/EERC-83/17. Berkeley: University of California, Earthquake Engineering Research Center, 1984.

[56]              Steelman JS, Hajjar JF. Influence of Inelastic Seismic Response Modeling On Regional Loss Estimation. Engineering Structures, 2009,31 (12):2976-2987.

[57]              Ibarra LF, Medina RA, Krawinkler H. Hysteretic Models that Incorporate Strength and Stiffness Deterioration. Earthquake Engineering & Structural Dynamics, 2005,34 (12):1489-1511.

[58]              Chopra AK. Dynamics of Structures. Prentice Hall New Jersey, 1995.

[59]              GB 50009-2001 建筑结构荷载规范. 北京: 中国建筑工业出版社, 2002.

[60]              Sobhaninejad G, Hori M, Kabeyasawa T. Enhancing Integrated Earthquake Simulation with High Performance Computing. Advances in Engineering Software, 2011,42 (5SI):286-292.

[61]              GB 50011-2010 建筑结构抗震设计规范. 北京: 中国建筑工业出版社, 2010.

[62]              马玉虎, 陆新征, 叶列平, 唐代远, 李易. 漩口中学典型框架结构震害模拟与分析. 工程力学, 2011(05):71-77.

[63]              Miao Z, Ye L, Guan H, Lu X. Evaluation of Modal and Traditional Pushover Analyses in Frame-Shear-Wall Structures. Advances in Structural Engineering, 2011,14 (5):815-836.

[64]              Tang B, Lu X, Ye L, Shi W. Evaluation of Collapse Resistance of RC Frame Structures for Chinese Schools in Seismic Design Categories B and C. Earthquake Engineering and Engineering Vibration, 2011,10 (3):369-377.

[65]              施炜, 叶列平, 陆新征, 唐代远. 不同抗震设防RC框架结构抗倒塌能力的研究. 工程力学, 2011(03):41-48.

[66]              Krawinkler H. Van Nuys Hotel Building Testbed Report: Exercising Seismic Performance Assessment.Pacific Earthquake Engineering Research Center, College of Engineering, University of California, Berkeley, 2005.

[67]              Hori M, Ichimura T. Current State of Integrated Earthquake Simulation for Earthquake Hazard and Disaster. Journal of Seismology, 2008,12 (2):307-321.

[68]              Owens JD, Luebke D, Govindaraju N, Harris M, Krüger J, Lefohn AE, Purcell TJ. A Survey of General‐Purpose Computation On Graphics Hardware, 2007.

[69]              Stone JE, Phillips JC, Freddolino PL, Hardy DJ, Trabuco LG, Schulten K. Accelerating Molecular Modeling Applications with Graphics Processors. Journal of Computational Chemistry, 2007,28 (16):2618-2640.

[70]              Weldon M, Maxwell L, Cyca D, Hughes M, Whelan C, Okoniewski M. A Practical Look at GPU-Accelerated FDTD Performance. Applied Computational Electromagnetics Society Journal, 2010,25 (4SI):315-322.

[71]              Lv MH, Wei X, Lei C. A GPU-Based Parallel Processing Method for Slope Analysis in Geographical Computation. Advanced Materials Research, 2012,538:625-631.

[72]              Che S, Boyer M, Meng J, Tarjan D, Sheaffer JW, Skadron K. A Performance Study of General-Purpose Applications On Graphics Processors Using CUDA. Journal of Parallel and Distributed Computing, 2008,68 (10):1370-1380.

[73]              Hung SL, ADELI H. A Parallel Genetic/Neural Network Learning Algorithm for Mimd Shared-Memory Machines. IEEE Transactions On Neural Networks, 1994,5 (6):900-909.

[74]              Mackerle J. FEM and BEM Parallel Processing: Theory and Applications - a Bibliography (1996-2002). Engineering Computations, 2003,20 (3-4):436-484.

[75]              Adeli H. High-Performance Computing for Large-Scale Analysis, Optimization, and Control. Journal of Aerospace Engineering, 2000,13 (1):1-10.

[76]              Cundumi O, Suarez LE. Numerical Investigation of a Variable Damping Semiactive Device for the Mitigation of the Seismic Response of Adjacent Structures. Computer-Aided Civil and Infrastructure Engineering, 2008,23 (4):291-308.

[77]              NVIDIA CUDA Programming Guide. http://docs.nvidia.com/cuda/pdf/CUDA_C_Programming_Guide.pdf.

[78]              FEMA P695 - Quantification of Building Seismic Performance Factors. Washington, D.C.: Federal Emergency Management Agency, 2009.

[79]              Abdalla R, Tao V. Integrated Distributed GIS Approach for Earthquake Disaster Modeling and Visualization//Berlin: Springer, 2005:1183-1192.

[80]              Duzgun HSB, Yucemen MS, Kalaycioglu HS, Celik K, Kemec S, Ertugay K, Deniz A. An Integrated Earthquake Vulnerability Assessment Framework for Urban Areas. Natural Hazards, 2011,59 (2):917-947.

[81]              Xu Z, Lu X, Guan H, Han B, Ren A. Seismic Damage Simulation in Urban Areas Based On a High-Fidelity Structural Model and a Physics Engine. Natural Hazards, 2014,71 (3):1679-1693.

[82]              Millington I. Game Physics Engine Development (The Morgan Kaufmann Series in Interactive 3D Technology). San Francisco: CRC Press, 2007.

[83]              Boeing A, Bräunl T. Evaluation of Real-Time Physics Simulation Systems, 2007.

[84]              Joselli M, Clua E, Montenegro A, Conci A, Pagliosa P. A New Physics Engine with Automatic Process Distribution Between CPU-GPU, 2008.

[85]              Maciel A, Halic T, Lu Z, Nedel LP, De S. Using the PhysX Engine for Physics-Based Virtual Surgery with Force Feedback. International Journal of Medical Robotics and Computer Assisted Surgery, 2009,5 (3):341-353.

[86]              Xu Z, Lu X, Guan H, Ren A. Physics Engine-Driven Visualization of Deactivated Elements and its Application in Bridge Collapse Simulation. Automation in Construction, 2013,35:471-481.

[87]              Tsai F, Lin HC. Polygon-Based Texture Mapping for Cyber City 3D Building Models. International Journal of Geographical Information Science, 2007,21 (9):965-981.

[88]              OSG Community, OpenSceneGraph. [2012.05.26]. http://www.openscenegraph.org/projects/osg/.

[89]              Havok. Havok Physics. [2012.07.24]. http://www.havok.com/products/physics.

[90]              Game Physics Simulation, Bullet Physics Library. [2012.07.13]. http://www.bulletphysics.com.

[91]              NVIDIA Corporation, PhysX. [2012.08.11]. http://www.nvidia.cn/object/physx_new_cn.html.

[92]              Reitherman R. Earthquakes and Engineers: An International History. American Society of Civil Engineers, 2012.

[93]              Roesset JM, Yao J. State of the Art of Structural Engineering. Journal of Structural Engineering-ASCE, 2002,128 (8):965-975.