第二十三届全国高层建筑结构学术会议论文2014年 高层建筑弹塑性动力时程分析的新型软件方法 李志山，和雪峰²，曹胜涛2 （1.广州建研数力建筑科技有限公司，广州510170：2.广州大学，广州510006） 提要：准确分析高层建筑在不同地震水准作用下的工作性能，是高层建筑抗震性能化设计的核心，各设计分 析软件所采用的分析模型和算法简繁不同.
本文介绍了一套基于全新GPUCPU硬件架构的干万自由度级高层 建筑结构动力弹塑性细粒度并行分析软件一PKPM-SAUSAGE.
软件采用健壮的有限单元模型、先进混凝土弹塑 性本构和高效的时程积分方法，实现了精细化大规模结构非线性地震动力时程响应数值仿真的目标.
软件与常 用设计软件无缝衔接，使复杂的大震分析成为善通结构工程师能快速胜任的桌面化工作任务.
关键词：高层建筑，弹塑性动力时程分析，PKPM-SAUSAGE，异构计算 1前言 性能化抗震设计要求对结构及构件在地震作用下的非线性性能作细致的量化评估".
我国《建筑抗震 设计规范》和《高层建筑混凝土结构技术规程》对建筑结构的地震位移、变形和内力等性能指标提出了量 化的控制要求.
要实现这些性能设计目标，理想的方法是实施精细的结构弹塑性动力时程分析，这已成为 研究和应用领域的共识.
然而，如何实施有效的结构弹塑性分析，以获取尽可能细致可靠的结构抗震性能 评价，长期以来却不得不面对来自分析模型、分析方法、软硬件架构等多方面的限制因素.
2高层建筑结构精细化非线性地震时程分析的思路 高层建筑结构非线性地震时程分析的模型、方法和软硬件架构，是相互影响的三个方面.
早期较为简 化的模型与算法，主要强调工程概念的合理性与正确性，来源是对试验现象和震害调查的总结.
随着计算 能力的增强，出现了更加灵活细致的算法和模型，设计分析能够掌握的性能指标也趋于多样面具体.
进行结构地震受力性能的非线性全过程分析，首先必须有往复荷载下材料或截面性能的本构关系，也 即地震恢复力模型2.
以往，钢筋混凝土结构地震恢复力模型多定义在楼层或构件的宏观层次.
楼层的恢 复力模型适于整体结构的抗震概念分析.
梁和柱的构件恢复力模型主要描述塑性铰的地震滞回受力特性， 目前在不少设计软件中均有应用，有双线性、三折线等常用形式.
剪力墙的构件恢复力模型主要应用于等 效梁、多垂直杆、多弹簧元件等剪力墙简化分析模型中.
这些宏观层次的地震恢复力模型具有概念简明的 特点，通常有着较为丰富的试验依据，但往往又容易受到试验情况的影响和限制，难以推广到复杂多样的 应用形式，并还存在多轴耦合分析等困难.
如果从整体结构的分析方法看，塑性较等宏观层次地震恢复力模型的产生和应用与传统结构矩阵分析 方法有密切关系.
在结构矩阵位移分析等传统结构矩阵分析中口，一个构件通常就是一个分析单元，刚度 矩阵也以构件为基础显式写出，弹塑性刚度往往通过修正构件弹性刚度而获得.
随着有限元结构分析技术 的发展和普及应用，一个构件通常划分多个单元，单元又采用多项式插值和积分处理，使描述更为细致的 构件弹塑性状态成为可能，同时也存在此必要.
为此，恢复力模型至少需要建立在截面层次.
截面层次恢复力模型有两种类型口.
其一，是依据试验建立构件截面位移、变形与内力的关系，如Coulgh 三线型等，但存在的间题与塑性较模型类似，而对板壳构件尤为不可行.
其次，是截面积分方法：对截面 基金项目：长江学者和创新团队发展计划资助，编号：IRT13057 作者黄合：李志山（1966一），男，博土，教授.
国家千人计划特聘专家.
第二十三届全国高层建筑结构学术会议论文2014年 上的点应用弹塑性本构，经截面积分（求和）获得刚度、内力等宏观物理量及其非线性关系.
采用截面积 分的梁通常被称为纤维梁，板壳称为分层壳”，它们早在上世纪五六十年代就被应用于结构弹塑性分析， 只是圈于当时的计算条件未能得到充分发展.
在理论上，梁、板、壳构件的截面刚度和内力，正是基于平 截面假定对材料点模量与应力的积分结果，因此这一方法是现代计算条件下对构件基本力学模型与方法的 回归.
截面积分方法还被用于传统塑性铰模型中塑性区域截面内力-变形关系的计算，称之为纤维铰模型.
当然，截面积分方法也存在局限，受平截面假定的限制，它不便体现钢筋与混凝土的相互作用，是一种刚 度组合式的钢筋混凝土材料非线性分析方法，适宜描述的弹塑性状态也以压弯耦合为主.
但另一方面，应 力-应变层次材料本构的研究显然更容易实施，而且成果丰富，便于应用：积分计算又不受多种材料组份及 具体截面形式的限制，适宜面豁然拓宽.
因此，截面积分方法是当前梁（柱）、板、壳宏观构件实现应力- 应变层次上精细化分析的有效途径.
早在2004年，本文第一作者利用现代CAE软件ABAQUS的梁、板、壳构件截面动态积分功能，研 究开发了基于ABAQUS的大规模结构动力弹塑性分析成套技术，实现了高层结构在应力-应变-材料损伤层 次上的精细化分析，先后应用于国内CCTV新大楼、上海环球金融中心等数十项重要高层、超高层结构的 罕遇地震抗震性能评估，获得良好效果”.
但这一成套技术，也有其局限性，主要问题是应用门槛相对过 高、软硬件成本较高，不利于广泛推广.
结构大规模精细化弹塑性分析必然导致计算规模的显著增长，更追论地震动力时程分析中数以万计时 步的任务重复.
这对计算机软硬件架构提出了新的发展要求.
基于对以往成功经验的反思，本文第一作者 根据计算机硬件技术的最新发展动态，研究开发了基于CPUGPU高性能异构架构的大规模结构精细化弹 塑性动力时程分析新型软件系统PKPM-SAUSAGE（以下简称SAUSAGE）.
测试和验证表明，SAUSAGE 的算法策略健壮高效，分析结果精细可靠，软件组织具有实用性和专业针对性，并在弹塑性动力时程的规 模和速度方面取得明显突破.
3大规模动力弹塑性分析软件SAUSAGE的技术特点 3.1CPUGPU高性能并行计算 并行计算是超大规模结构计算分析的必然技术途径.
传统的并行计算，主要是基于微机集群 （PC-Cluster）的并行模式，主要基于任务分块，通过消总传递编程接口（MPI）实现微处理器（CPU）进 程之间的组织和通信.
CPU处理任务的能力固然强大，但从硬件发展的进程看，决定CPU处理任务速度 的时钟频率进一步快速提升已面临瓶颈，曾经预言的摩尔定律不复重现.
在此情况下，要通过PC-Cluster 并行模式提高问题求解规模和速度，意味着要增加机群建设的成本.
以ABAQUS多CPU并行分析为例， 根据笔者经验，百万自由度规模的结构模型，采用ABAQUS-Explicit实施30秒地震波作用下的弹塑性动 力分析，若要在6小时以内完成，至少需配置48个CPU，单核频率不低于2.4G，内存不小于32G，再加 上维护成本和使用技术门槛，并非一般设计单位所能普遍承受.
CPUGPU高性能异构计算是近年兴起的新型并行计算技术.
CPU擅长处理不规则结构和不可预测 的数据存取模式，以及递归、分支密集型和顺序执行型的任务.
GPU，即显卡，则擅长处理规则数据结构 和可预测存储模式，任务处理具有天然的细粒度并行特征，通常是千万个线程同时并发，在光影处理、医 疗成像、基因分析等领域有广泛应用.
如果让CPU处理分析流程的管理调度和数据高速缓存，GPU则处 理规则简单的底层计算，则将能够实现优势互补，可使整体分析效率获得数量级的提升，并只需要市面常 见的单机配置-Nividia公司最先发展了基于GPU的高性能计算技术，推出了CUDA编程平台.之后，Nividia 和AMD等公司均推出支持新一代支持CPUGPU高性能异构计算的编程平台OpenCL.
SAUSAGE是基于CPUGPU异构架构开发的新型结构分析软件，并行模式具有显著的细粒度特点.
为在CPUGPU架构上实现高层建筑大规模精细化弹塑性动力时程分析，SAUSAGE对结构弹塑性动力时 程分析的有限元模型、算法和数据结构等进行了系统的重构与创新.
对此，SAUSAGE规划了以下六个方
第二十三届全国高层建筑结构学术会议论文2 2014年 面的GPU线程细粒度并行计算： ①材料点弹塑性状态并行分析：对梁、板、壳截面上材料点的弹塑性状态分析实施GPU线程细粒度 并行分析，包括单轴或多轴的应力-应变弹塑性分析： ②截面弹塑性状态并行分析：对单元积分截面的弹塑性状态实施GPU线程细粒读并行计算，主要是 对材料点分析结果的截面积分（求和）计算： ③单元弹塑性状态并行分析：对单元的弹塑性矩阵与向量实施GPU线程细粒度并行计算，主要是单 元内的积分计算： ④自由度组集的并行计算：根据结点和单元的关联性，对整体矩阵和向量实施GPU线程组集计算， 主要是求和计算： ③代数方程组的选代求解：对静力计算和特征值求解中的代数方程组问题，采用多波前预处理共轭梯 度和预处理共轭梯度法，实施自由度级的GPU线程细粒度并行选代求解： 动力方程的显式积分求解：对整体动力方程，采用基于中心差分的细小稳定步长显式递推格式，实 施自由度级的GPU线程细粒度并行求解.
上述由GPU执行的分析任务，根据GPU的架构特点对数据格式进行优化设计，避免出现错误和低效 的数据访存，计算整体流程的调度和大块数据的管理则由CPU完成.
总体上，SAUSAGE实现了千万数量 级线程任务的快速并行分析.
3.2SAUSAGE的截面类型和材料模型 SAUSAGE中杆件的截面按用途分为以下几类：钢筋混凝土、型钢、钢管混凝土、钢骨混凝土的梁截 面、柱截面、斜撑截面，以及方钢管形式的边缘构件、虚梁等辅助构件截面，同时也允许用户自定义截面 形式.
弹塑性动力时程分析中，杆件截面均采用截面精细积分处理，混凝土采用单轴损伤塑性本构， 钢筋、型钢及钢管等采用单轴双线性随动强化本构.
SAUSAGE中钢筋混凝土墙板的钢筋以单向配筋率作为输入参数，仍采用单轴的双线性随动强化本构.
内嵌钢板及纯钢板墙采用基于Von.Mises屈服准则并考虑随动硬化的金属材料平面本构，屈服后为理想弹 塑性，卸载不考虑刚度退化.
墙、板构件的混凝土材料采用平面的损伤塑性本构模型，但截面积分处理有 所差别.
板构件混凝土的面内受力（膜内力），按整体一层进行弹塑性分析，而面外（弯曲内力）则保持 弹性.
剪力墙的混凝土则细分多层（默认6层），每层均按平面应力模型进行弹塑性分析，以完整体现面 内面外受力的弹塑性耦合效应.
这样的处理，符合墙、板构件各自的受力特点，也符合工程概念.
关于混凝土材料损伤塑性本构，近年来在学术和工程应用领域获得了广泛关注，而ABAQUS软件也 应用了这一混凝土模型.
混凝土塑性损伤本构理论是由Lubliner等l和Lce等首先建立和完善的，它能 好地反应在不大于55（5倍单轴极限压应力）范围内的静水压力下混凝土材料的力学行为，适用于建筑结 构中通常使用的混凝土材料.
其基本特征如下：1）反映不同的压拉屈服应力（受压屈服应力约为受拉屈 服应力的10倍甚至更高）：2）反映受拉软化行为、受压硬化及软化行为：3）反映不同的拉压刚度退化特 征：4）反映往复加载时的刚度恢复现象：5）具有率敏感性，特别是应变速率较大时峰值强度有所提高.
值得一提，《混凝土结构设计规范GB50010-2010》依据我国学者的研究成果，在附录中加入了一种混 凝土损伤塑性模型.
附录中混凝土损伤塑性模型骨架曲线上一点的应力应变关系如下式： 6= (1-d)Ee (1) 其中，d为损伤指标，E为初始模量，6、e为应力和应变.
易知，上式对损伤指标的定义利用了该点的割 线模量E（=o/e）： d = 1 - E / E (2) 而Lubliner等应力应变关系采用下式： = (1 - d ) E (c - e) (3) 其中为塑性应变或等效塑性应变，含义是采用该点的弹性卸载刚度E（=o/（e-）定义损伤指标：
第二十三届全国高层建筑结构学术会议论文 2014年 d = 1 -E / E () 为与GB50010-2010保持一致，PKPM-SAUSAGE采用了附录的损伤指标定义式（1）形式，骨架曲线 及其参数也与附录一致.
而在弹塑性状态分析中，则采用Lubliner等所提出的屈服函数、非关联流动法则 和Drucker-Prager流动势函数，应力更新采用完全隐式的Euler向后积分算法.
此处恐繁不述.
单铂受控 单轴受压 双轴受拉 -agpi-a 双轴受压 ap)=a 图1混凝土平面损伤塑性模型屈服面 3.3SAUSAGE的单元类型 低阶次有限单元模型在弹塑性分析中容易获得较好的数值稳定性，精度可以通过型加密得到保证.
为此，SAUSAGE所采用的单元基本上以线性及双线性单元为主，通常选择的单元尺寸也在0.5m左右.
表 1列出了软件中各类构件的适用单元类型.
表1SAUSAGE构件与单元模型对应表 构件类型 单元类型 主梁、柱、斜撑 2结点12自由度厚薄通用梁单元 次梁 2结点12自由度欧拉梁单元 一端钦接梁 2结点较接梁单元 边缘构件、虚梁、连梁纵筋、两端铰接梁、两端铰接斜撑2结点桁杆单元 楼板、剪力墙 3结点18自由度厚薄通用平板壳元 4结点24自由度厚薄通用平板壳元 表1中梁和板壳单元的结点均具有6个自由度，能够全面地描述构件的空间受力行为.
其中，板壳单 元普遍考虑了面内旋转自由度.
关于面内旋转自由度，一般认为主要用于描述板壳的面内刚体转动行为.
旋转自由度的参与与否可使分析结果存在不小的差异，对基底剪力和楼层剪力的影响尤其不容忽视.
对此， SAUSAGE依据Belytschkol2和龙驭球3等的研究工作，参照ABAQUS等成熟软件，对旋转自由度作了 合乎工程概念的假设和处理，并进行了测试验证.
SAUSAGE中，除保持弹性假设的受力行为外，进行弹塑性计算的有限单元模型均采用了单点积 分.
单点积分不仅有助于够提高单元的稳定性和通用性，而且对大规模并行分析意义重大，一是可减少历 时量存储，二是能提高算法的并行度.
但单点积分必然导致附加零能模式（即沙漏模态），使计算结果失 真.
这一间题，在上世纪八九十年代就引起了Belytschko等有限元研究者的重视，提出了多种附加沙漏 刚度矩阵和沙漏力的方法，不仅使零能模态得到有效抑制，同时也使单元模型能够在单点积分的情况下获 得很高的计算精度.
SAUSAGE采用先进的方法对单点积分单元进行了有效的沙漏控制，并进行了严格验 证.
第二十三届全国高层建筑结构学术会议论文2014年 为保证剪力墙和楼板的有限元网格剖分质量，SAUSAGE除控制单元尺寸，还采用高效的网格自动生 成技术.
实施预处理后的设计模型，利用Delauney方法或循环内缩铺路法及结构化四边形网格剖分方法进 行多次细化选代，最后再由Laplace算法实施优化，效果示意如图2.
图2最终网格效果示意 3.4SAUSAGE的动静力求解方案 SAUSAGE具有高层结构模态分析和分层加载施工过程仿真的功能，这些任务涉及到大规模线性代数 方程组的求解问题.
基于矩阵分解的传统直接解法，算法并行度通常较低，为挖掘实现GPU线程细粒度 并行计算的潜在效率，软件采用选代法解法求解大规模线性代数方程组问题，研发了两种基于预处理共轭 梯度法的选代求解器，并对CPU调度和GPU访存做了专门优化.
测试表明，SAUSAGE的大规模线性方 程组选代求解器的计算效率高，并仍有进一步优化和提升的空间.
动力方程的直接积分求解是结构弹塑性动力时程分析的核心.
动力方程的求解向来有隐式和显式之 分，关键区别在于时域差分的形式".
虽然Newmark等传统隐式解法在工程中的应用由来已久，但在大 个大规模非线性代数方程组，存储量和计算量极大：二是弹塑性选代的收敛性难以控制，固然可以通过弧 长法等加以改善和控制，但至今仍缺乏具有通用性的解决途径.
弹塑性间题的隐式求解在形式上虽也较符 合工程师所谓“平衡方程求解"的概念，然而往往由于不易实现顺利求解，导致最终结果存在准确性隐患.
与隐式方法相比，显式直接积分方法表现为相邻时步间的递推计算，刚度矩阵不必显式存储，而是转 换为不平衡力.
如果采用集中质量矩阵，还可实现自由度之间的解耦，大大有利于细粒度并行计算.
然面， 显式直接积分方法在时域上是条件稳定的，因此必须要求差分步长小于临界稳定步长，亦即差分格式的谱 半径，它与结构刚度的特征参数有关.
以加速度二阶中心差分为代表的显式直接积分方法在波动问题求解 中的应用也由来已久，ABAQUS-Explicit、LS-DYNA等CAE软件还将显式动力求解方法应用于高速碰撞、 动力接触及流变成型等高度非线性的问题中.
上世纪八九十年代起，Belytschko等l就针对有限元动力方 程的显式求解作了系统研究，其中包括对求解结果可靠性的论证.
为充分发挥GPU细粒度并行的优势，SAUSAGE采用基于加速度二阶中心差分的细小稳定步长显式直 接积分方法，所依据的基本方程如下： （1M M -c Mu 了1 Cu-s (r²2△ r²）（△r²2r） (5) 式中，M为集中质量矩阵，C为阻尼矩阵，Q.为等效外力（包括基础加速度激励），K为刚度矩阵，Ku 即为：时刻等效内力.
其中，SAUSAGE中的阻尼C允许采用Rayleigh阻尼，也可采用进行了时频转换的 模态阻尼，后者可避免了通常的刚度阻尼问题困扰.
上述算法均经过了严格验证.
在动力求解过程中，SAUSAGE采用UL（UpdatedLagrange）格式考虑了大位移几何非线性效应，能
