1.   引言

建筑结构连续倒塌是指由于偶然作用造成结构局部构件破坏，继而引起其余构件相继破坏，最终形成与初始局部破坏不成比例的结构大范围倒塌或者整体建筑倒塌的现象，属于一种低概率高损失的结构安全事故。建筑结构连续倒塌过程从力学上讲是一个从连续体到非连续体动态转化的非线性动力接触碰撞问题^[1]，涵盖物理、几何、接触多重非线性，其计算难点在于建筑材料由弹性到弹塑性进而断裂发展的材料非线性计算，梁板构件在材料非线性下的大位移大转动计算，以及梁板断裂破坏引起的碰撞接触非线性计算，涉及材料弹塑性及应变率效应、构件大位移大转动、接触碰撞等问题，带有明显的撞击—接触特点。鉴于建筑结构抗连续倒塌对建筑结构安全的重要意义，有必要针对性的研究相应的非线性计算理论和数值模拟方法。

2.   数值计算工具

现阶段常用的有限元软件主要分为两类，一类是隐式有限元，另一类是显式有限元，前者基于刚度矩阵建立平衡方程，后者则基于内外力平衡进行时间积分，常用的计算软件汇总见表 1。

表  1  商用有限元软件列表

软件分类	商用有限元软件
分析工具	隐式	显式
	ANSYS ADINA MARC SAP2000 ABAQUS/STANDARD	LSDYNA AUTODYN DYTRAN ABAQUS/EXPLICIT

下载: 导出CSV 

| 显示表格

对于力学上具有接触碰撞非连续体受力特点的连续倒塌问题，因为这些软件的求解能力各异，间接导致连续倒塌计算被分成两类，一类是连续倒塌的风险评定，另一类是连续倒塌过程的模拟。前者通过假定倒塌准则，一旦超过倒塌准则，即认为结构发生连续倒塌，典型的为弹性静力分析方法，其突出优点在于计算过程简单，计算效率高，其缺点就是将连续倒塌分析更多的引向了经验判断，对倒塌准则的依赖程度较高，结果的可信度和说服力有限；后者则通过全过程数值模拟，实现倒塌的范围和破坏程度鉴定，典型的为非线性动力分析方法，通过数值模拟结果来判断结构的损坏程度，其突出优点在于结果是对破坏过程的仿真再现，能够研究倒塌过程和倒塌模式，缺点就是对分析软件的计算能力要求很高，计算效率也相对较低，结果分析工作较为复杂，结果更依赖于数值计算结果的准确性。

3.   常用数值模拟方法

3.1   离散元方法

离散元法(Discrete Element Method，DEM)又称刚体有限元，是20世纪70年代初由Cundall首先提出的^[2]。最初，它的研究对象主要是岩石等非连续体的力学行为，其基本思想是把不连续体分离为刚性元素的集合，各个单元满足牛顿第二定律，用时步更新的方法求解各刚性元素的运动方程，继而求得不连续体的整体运动形态。离散单元法允许单元间的相对运动，不一定要满足位移连续和变形协调条件，尤其适合于大位移和大变形的非连续介质问题的求解。

离散元算法的一般求解过程为^[3]：首先将求解空间离散为离散元单元阵，然后根据实际问题用合理的连接件将相邻的两单元连接起来，单元间相对位移是基本变量，有利于相对位移的关系可以得到两个单元间法向方向和切向方向的作用力；对单元在各个方向上与其它单元的作用力以及其它物理场对单元的作用所引起的外力求合力和合力矩，根据牛顿第二定律，可以求得单元的加速度，对加速度进行时间积分，即可依次得到单元的速度、位移，进而得到所有单元在任意时刻的速度、加速度、角速度、角加速度，位移和转角等物理量。

采用离散元方法模拟建筑倒塌，可以充分利用该方法较为成熟的非连续体计算能力，因离散元已经建立较多的接触模型，能够较好地处理建筑结构倒塌非连续体阶段的接触碰撞计算；但是由于离散元无法采用连续介质力学的应力应变本构理论，只能采用在构件层次上的恢复力模型，虽能处理连续体的力学分析，但却无法应用有限元方法成熟的材料本构理论，所以离散元也需要进一步改造，以适应计算的需要。

3.2   有限元方法

有限元方法的基础是变分原理，在每个单元内，将微分方程中的变量改写成由各变量或其导数的结点值与所选用的插值函数组成的表达式，借助于变分原理，将微分方程离散求解。采用不同的插值函数形式，便构成不同的有限元方法。

3.2.1   隐式有限元

隐式有限元方法可以说是目前最为广泛流行的有限元方法，其基本思路和解题步骤可归纳如下^[4]：

(1) 按照最小势能原理建立积分方程，根据变分原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。

(2) 区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接的单元。

(3) 确定单元基函数，根据单元中结点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元基函数。

(4) 单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数(即单元中各结点的参数值)的代数方程组，称为单元有限元方程。

(5) 总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进行累加，形成总体有限元方程。

(6) 边界条件的处理：一般边界条件有三种形式，分为本质边界条件(狄里克雷边界条件)、自然边界条件(黎曼边界条件)、混合边界条件(柯西边界条件)。对于自然边界条件，一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法则对总体有限元方程进行修正满足。

(7) 解有限元方程：根据边界条件修正的总体有限元方程组，是含所有待定未知量的封闭方程组，采用适当的数值计算方法求解，可求得各结点的函数值。

对于几何非线性分析，隐式有限元首先通过虚位移或虚功原理来推导平衡方程，然后进行有限元离散，针对选用的单元列式，计算相应的单元刚度矩阵(弹塑性矩阵、几何刚度矩阵和初应力矩阵)，进而迭代求解相应的非线性方程。在处理连续倒塌过程的几何材料双重非线性并带有构件断裂破坏问题时，刚度矩阵将发生奇异，必然引起收敛困难，另一方面一旦处理接触碰撞问题，势必要求较小的时间步长，甚至接近显式有限元的积分步长，将无法发挥隐式有限元时间积分大步长无条件稳定的优点，因而隐式有限元应用于连续倒塌评定工作具有突出优势，即进行弹性静力、弹性动力及非线性静力(位移加载)计算，而应用于连续倒塌的模拟分析工作，将较难发挥隐式有限元的优势，需要优化算法，采用高性能计算。

3.2.2   显式有限元

显式有限元的理论研究开始于20世纪60年代中期，为了尽可能减少计算量，提高计算效率，显式有限元逐渐建立了一套与显式时间积分算法相配套的单元构造理论，这些单元支持所有非线性选项，并且由于动力平衡方程建立于当前时刻，利用当前时刻的平衡条件获得下一时刻的位移、速度和加速度反应值，不需要迭代计算即能不断推进计算，因而显式有限元拥有极为强大的非线性计算能力，它具有如下突出特点：对于材料非线性计算问题，材料本构能够考虑弹塑性、损伤、断裂破坏，也能够考虑应变率效应，还能够考虑状态方程；对于几何非线性计算问题，能够处理大位移大转动且涉及材料非线性构件失稳的双重非线性问题；对于接触非线性分析，具有高效的接触判断和计算效率，不需要在接触与不接触的反复迭代试算；由于控制方程采用显式时间积分决定了不存在隐式迭代计算的收敛性问题，而且计算量与单元数目成线性关系；显式有限元采用应力积分的方式求解内力，因而避免了计算刚度矩阵，在计算流程上拥有更强健更稳定的能力^[5]。

但是，显式有限法对于连续倒塌模拟存在的问题有：1、目前商用显式有限元的梁元，普遍缺少合适的一维混凝土弹塑性本构关系，即钢筋混凝土梁柱很难直接采用梁元模拟，需要二次开发约束混凝土的一维本构；2、显式有限元作为一种连续介质分析方法，采用“生死”单元造成随着结构计算持续，结构的单元数目越算越少，表现为总体质量偏小堆积高度偏低，基于连续介质力学的有限单元方法在一定程度上计算结构倒塌阶段的单元之间的碰撞和堆积情况失真；3、建筑结构由于自重和活荷载的原因存在初始应力，而显式有限元主要采取两种方式解决初应力的施加问题，一是隐式-显式结合，另一种是采用动力松弛法。显隐式结合的方法，有限元在转换过程中存在诸多限制，而动力松弛法的计算耗时和结果精度却无法保证。

针对隐式有限元和显式有限元的计算特点表 2进行了对比分析，由于显式有限元具有强健的计算能力，其计算原理适合于结构连续倒塌过程的这类强非线性分析，只是由于材料本构和单元的限制，对于普通用户使用还比较有难度，对于高级用户需要本构二次开发才能实现模拟；另一方面，由于显式有限元法的单元消失问题，对于连续倒塌过程，会引起接触撞击力计算不准或者堆积高度不足等问题，需要引进非连续介质力学等方法来改进显式有限元方法。

表  2  隐式有限元与显式有限元计算特点对比表

	隐式有限元	显式有限元
时间积分稳定性	无条件稳定算法	条件稳定算法
积分步长	大步长	小步长，受Courant稳定条件限制
精度	二阶精度	二阶精度
内存要求	高	低
CPU要求	高	低
健壮性	低，刚度矩阵容易奇异、非线性方程求解容易发散	高，能够处理高度非线性及多场耦合问题
求解域	全局求解，步步收敛	区域求解，逐个单元计算
适用性	静力问题、动力问题、弱非线性问题	拟静力问题、动力问题、强非线性问题

下载: 导出CSV 

| 显示表格

3.3   适用性分析

从上面的隐式有限元和显式有限元的原理分析可以看出，采用有限元方法模拟建筑倒塌，可以归纳出如下优点：

(1) 成熟的单元技术，直接能够考虑大位移大转动力学计算问题，可以较好解决几何非线性问题；

(2) 成熟的材料本构技术，能够继承既有的研究成果，可以较好解决材料非线性问题；

(3) 较为成熟的有限元计算软件，具有较好的计算平台，在计算连续体受力时精度较高。

但同时也面临如下问题：

(1) 单元破坏后的处理问题，由连续体进入非连续体后的力学计算，是有限元方法不具备的能力，“杀死”单元会降低计算精度，对于单元破坏较多的结构，常导致结构最终质量偏小，同时可能由于接触力计算精度不够而引起较大的计算误差，可以说对于非连续体的力学计算问题，是采用有限元方法模拟结构倒塌的较为明显的弱点；

(2) 接触计算问题，有限元所采用的梁、柱及实体单元的接触算法，最终均转化为点面算法，或者面面算法，在计算效率上相对较低，对于构件众多的建筑结构，处理大规模的接触计算，仍然具有较大难度。

综上，由于有限元针对的是可变形的连续体分析，而离散元针对的是非连续体分析，因而二者对于处理连续体与非连续体共同受力时均存在一定的局限。采用有限元方法虽然能够在一定程度上模拟建筑结构倒塌问题，但是仍然存在一些不足，需要进一步改进方法；采用离散元方法不能够继承成熟的大变形单元技术和本构技术，在一定程度上增加了研究的难度；采用显式有限元与离散元的耦合方法虽能进一步提高计算精度，但是由于是两套系统单独运行，在计算效率方面存在较大考验，因此有必要进一步研究数值计算方法。

4.   质点元方法

为了解决有限元法和离散元方法的不足，针对建筑结构连续性倒塌的力学特点，作者提出了能够考虑多重非线性用于连续介质与非连续介质共同作用且动态转化计算的质点元方法，图 1和图 2分别给出一个通过质点元方法进行的结构倒塌计算^{[6, 7]}。

图  1  某建筑结构倒塌过程

下载: 全尺寸图片 幻灯片

图  2  考虑接触和结构自重某框剪结构倒塌过程示意图

下载: 全尺寸图片 幻灯片

质点元方法将显式有限元和离散元统一为以牛顿第二定律的质点运动学为基础的同一方法，通过结点形状描述结构形状，通过定义广义连接模型、构造连接模型转化法则和建立接触碰撞算法，使得显式有限元与离散元具有相同的计算框架，实现了连续体到非连续体的动态计算，能够进行建筑结构倒塌等强非线性的动力计算。该方法从实现过程上可以看作放弃了有限单元的单元形状而增加了结点形状以描述结构外形，并且将离散元的接触模型扩展到有限单元，从而消除了显式有限元方法模拟结构倒塌因单元消失导致总体质量减小的问题，并要求接触算法均在质点元间进行，避免了显式有限元向离散元转换的耦合算法中的有限单元与有限单元、有限单元与离散单元、离散单元与离散单元三类接触算法共存现象，具有效率优先、精度可调的特点，适合大型建筑结构高效计算的需要。

5.   结语

虽然建筑结构连续倒塌的研究越来越多，但现阶段对于连续倒塌问题研究仍然有许多工作要做，针对建筑结构的连续倒塌非线性动力计算，研究开发计算软件系统的很少，目前主要依赖于国外商用软件，在计算力学层面滞后于国外；同时由于通用分析工具与规范结合度较低，对设计指导性较差，且操作复杂，门槛较高，针对性的数值计算研究较少；而且非线性动力计算缺少合适的分析工具，隐式有限元收敛困难，显式有限元单元消失严重，作者新提出的质点元方法可以改进显式有限元的计算精度。