动力弹塑性

动力弹塑性时程分析概述

孙雪帅 110940

摘要：本文对当前的弹塑性时程分析方法进行了概述，简要介绍了弹塑性动力时程分析的目的、原理和方法，对目前弹塑性动力时程分析软件的基本假定、模型、算法等方面作了扼要的描述和比较，。

关键词：动力弹塑性动力时程分析,ABAQUS，非线性算法。

结构抗震分析理论经历了从静力法到反应谱法再到弹塑性时程分析法的发展过程。反应谱理论考虑了结构动力特性和地震特性之间的动力关系，使结构动力特性对结构地震反应的影响得以体现，但是在进行结构抗震设计时，它仍然把地震惯性力作为静力来对待，无法准确反映地震对结构的实际影响。20世纪50年代末期，Houser实现了地震反应分析的动力计算方法，并将其成功应用与墨西哥城的拉丁美洲大厦设计，拉丁美洲大厦先后经过两次地震的考验，充分体现了动力时程分析法的优越性。动力弹塑性时程分析法应用于结构抗震设计主要是作为反应谱方法进行设计的补充手段。我国《建筑抗震设计规范》第5.2.2条规定，对于特别不规则的结构、板柱—抗震墙、底部框架砖房以及高度不大于150m的高层钢结构、7度三、四类场地和8度乙类建筑中的钢筋混凝土结构和钢结构宜进行弹塑性变形验算；对于高度大于150m的钢结构、甲类建筑等结构应进行弹塑性变形验算。

结构抗震设计的基本目标是“小震不坏，中震可修，大震不倒”。随着基于性态抗震设计方法的发展，根据不同建筑的安全需求与经济性等要求，按照性能化目标的思想，抗震设计目标在基本目标下被进一步细化和提高。一般来说，在安全与经济双重目标要求下，结构在小震状态下处于弹性状态，而且变形也较小，此时采用线弹性方法分析内力与变形误差较小，是可行的。在中震状态下，结构少部分构件已进入塑性状态且变形加大，此时若仍然采用线弹性的方法分析，则存在较大的误差。结构在打针状态下，大部分构件已进入塑性状态，并产生相当大的变形，其P-Δ效应加剧，几何非线性程度加大，所以计算分析不能采用线弹性方法，也不宜采用静力弹塑性方法，而应采用弹塑性动力时程分析方法。

结构分析的主要目标是获取结构的位移场、应变场及应力场，三者之间具有密切的关系，故我们 仅需获得结构位移场即可。通过离散化的方法，可将结构的弹塑性动力学方程表达如下：

...()MuCuKuFt

式中：u为节点位移向量，结构连续体的位移场可通过节点位移向量求得；M为质量矩阵；C为阻尼矩阵；K为刚度矩阵；F为外力向量函数；t为时间变量。地震作用时，..gFMu，其中..gu为地面运动加速度，即地震波。

由于在外力作用下，结构可能具有几何非线性与本构非线性（弹塑性本构是非线性本构中的一 种），结构的形态、刚度矩阵及阻尼矩阵不断变化， 使得上述方程的求解非常复杂。就刚度矩阵而言，它是由单元刚度矩阵eK组装而成。根据弹塑性力学，单元刚度矩阵可表达为：

[][][]TeeKBDBdv

式中：[]B为几何矩阵，通过几何矩阵[]B，由位移可求得应变；[]D为本构矩阵，通过本构矩阵[]D，由应变可求得应力。由于非线性效应，[]B、[]D矩 阵是不断变的。对弹塑性问题而言，一旦知道任 何时刻的几何矩阵、本构矩阵，通过积分点的数值 积分，即可得到单元刚度矩阵。也就是说，各积分点无论是处于弹性或塑性状态，我们都可以得到对 应时刻的单元刚度矩阵。再通过边界条件，即可逐 步求解得到节点位移向量，进一步可求得任意一处的位移、应变及应力，实现分析的目标。在上述离 散化的过程中，最一般的单元是三维实体单元，其 位移模式可以是线性的或者是二次的，视精度与效率的要求而定。在具体问题中，由于受力与变形机 制的特殊性，导致位移场具有一些特殊性，合理利 用这些特殊性，并作出相应的假定，可大大提高计算效率和精度。如采用直线法等假定形成板单元，采用平截面等假定形成梁单元等。 对弹塑性动力方程的求解，一般可分为两种求解算法：即隐式与显式。隐式算法常 采用 Newmark算法，它需要求解全区域的联立方程，求解过程复杂，并容易 导致结果不收敛的情况 。显式算法采用中心差分法，对动力学方程进行时间积分，由一个时间增量步的动力学条件下求解下一时间增量步的动力学条件，当时间增量充分小时，不会产生结果不收敛的情况，可获得问题的解答，因此显式算法特别适合弹塑性动力时程分析， 它的基本过程如下：

先将弹塑性动力方程改写如下 ：

...()MuFtCuKu

记P= ()Ft,I= .CuKu

则：

..Mu=P-I

第一步 ：节点变量求解

① 动力平衡方程

1()()tttuMPI

上式中ܯM为对角矩阵， 易于求解。

② 对时间显式积分

....()()()()()22()2ttttttttttuuu

.()()()()2tttttttuutu

第二步：单元应力计算

① 由节点位移，可求单元位移场，由此可求出应变场及应变率场，进一步求

出各处应变增量d。

② 根据本构关系计算应力 ()()(,)tttfd

③ 由I=.CuKu可求节点内力()ttI

第三步：将时间t变为()tt，返回第一步。

有上述内容可知，地震作用下弹塑性动力时程分析的主要关键点是材料本构关系的建构、单元选用、大变形的描述、地震波选取、阻尼选取及算法选择。

BAQUS弹塑性动力时程分析：

ABAQUS 是国际上ܯ先进的大型通用有限元分析软件之一，它的非线性功能达到世界领先水平，从而使其具有力学系统仿真的功能，它广泛应用于工程和科研各个领域。近年来，我国工程界将它成功 应用于弹塑性动力时程分析，其分析

结果得到业内专家们广泛认可。

1 材料本构关系

ABAQUS自带丰富的本构关系模型 ，可描述混凝土、钢材、岩土、高分子材料等物质的应力与应变关系。另外，ABAQUS还提供用户材料接口程序UMAT

及 UVMAT，因此可使用户自定义材料本构关系。对建筑结构来说，主要涉及混凝土 与钢材两种材料，钢材本构关系可采用二折线或三折线塑性本构关系，由于钢 材质地均匀 、性能稳定，其动力滞回模型也较为简单，下面重点介绍混凝土的本构关系。ABAQUS软件中，混凝土本构关系模型有混凝土弥散开裂模型、混凝土开裂模型及混凝土塑性损伤模型。其中混凝土塑性损伤模型（Concrete

Damage Plasticity）可描述混凝土受动力往返作用下的力学行为，故而故而广泛用于地震作用下的弹塑性动力时程分析。下面通过考察单轴特征荷载作用下的应力与应变关系曲线来把握混凝土塑性损伤模型的主要概念。

如图所示，曲线OAC代表单轴单向受拉的应力-应变曲线，曲线OLK代表单轴单向受压的应力-应变曲线，曲线OABDEFHIJ代表单轴往返荷载作用下的应力-应变曲线。直线OA段代表单轴受拉的弹性阶段，当应变增大并超过A点所对应的应变时，混凝土受拉屈服，并进入受拉塑性区；当应变增量为负时，混凝土浸入受拉卸载阶段，其应力-应变曲线近似视为直线，其斜率与混凝土在 B

点的损伤程 度相关，而受拉损伤程度（Dt）与 B 点的应变相关，混凝土受拉塑性应变越大，混凝土受拉损伤程度就越大，其卸载时的刚度也越小，这一点与我们的概念分析及实验结果是一致的。另外，在卸载直线范围内，若应变增量又为正时，则应力与应变值所构成的相点沿曲线DBC运动。

当相点越过 D 点进入受压状态时，混凝土受压时的刚度可得到一定程度的恢复。当应变增量继续为负时,相点沿弹性压缩曲线DE运动，当应变小于 E 的点对应值时，混凝土受压屈服，进入受压塑性状态。此后，若应变增量又转正时，混凝土进入受压卸载阶段，相点沿 FH 直线段运动，其直线斜率与 F 点的受压损伤程度(Dc)相关。当应变增量继续为正时，相点越过 H 点进入受拉状态，并沿直 线 HI 运动，此直线的斜率不断变小，不仅与受拉 损伤程度（Dt）相关，而且还与受压损伤程度(Dc)相 关。当应变增量继续为正时，混凝土越过 I

点再次 屈服进入受拉塑性状态。在上述过程中，值得指出的是，由于地震作用下的随机性，相点随时有可能作反向运动；随着时间的推移，混凝土的受拉与受压损伤程度是不断积累而增大的。上述内容是针对单轴应（即一维）情况的描述，对二维、三维情况，塑性损伤模型根据塑性势流动理论及单轴的应力‐

应变关系，由塑性应变增量求出应力增量及应力。

2 单元选用

从一般意义上来说，任何结构可采用三维实体单元进行分析，但也因此会带来收敛速度慢、计算效率低、后处理工作量大等问题，所以采用三维实体单元分析一般用于节点分析或结构的分析研究等方面。ABAQUS 软件中的三维单元可用于隐式与显式分析，对钢筋混凝土结构来说，三维实体元可使用前述所有混凝土本构关系模型 ，还可采用rebar 方式或 embedded 方式模拟混凝土中的钢筋。

ABAQUS 的壳单元可用来描述楼板、剪力墙及无扭转效应的梁，它可用于隐式与显式分析，并可采用各种混凝土的本构模型。可采用 rebar layer 方式在层中定义钢筋。ABAQUS 中考虑剪切效应的三维梁单元的使用受到一定限制，如它不能直接采用混凝土塑性损伤模型，这对 ABAQUS 直接应用于弹塑性动力时程分析的确算是一个遗憾，但我们可运用 VUMAT 接口程序的方法来弥补这方面的不足。在弹塑性动力时程分析中，ABAQUS 无需定义塑性区，只需在材料中定义其塑性行为特性即可，这样，单元的刚度由单元积分求得，结构的塑性区及其塑性损伤程度由外部作用、本构关系及其力学规律来确定，因此它能较真实地反映结构的实际受力与变形状态。

3 地震波选取

根据规范要求，地震波的选取首先应符合以下 要求：

①应按建筑场地类别和设计地震分组的相关 参数选用不少于两组天然波和一组人工波；

②地震波持续时间不宜小于12s，一般可取结 构基本自振周期的5-10倍，地震波数值化时距可取 0.01s或0.02s； ③地震波的ܯ大加速度应根据抗震设防烈度确定。

除此之外，地震波的频谱特性应与场地的构造 及其参数相协调。考虑双向水平地震和竖向地震同时作用，同一组地震波的两个水平分量和竖向分量的加速度比值可取1: 0.85: 0.65。

4 阻尼选取

在结构动力学中，阻尼的选取对计算结果较为敏感，故应谨慎对待。一般可采用瑞利（Rayleigh）阻尼模式，即对（1）式中的阻尼矩阵C表达如下：

CMK

式中：为质量阻尼系数；为刚度阻尼系数。

对多自由度力学系统，有如下关系：

(/)/2iii

min

min/

式中：i是系统园频率为i时的阻尼比。

根据上述关系可确定与。注意对用户定义的材料 ，只能选用质量阻尼系数 ，即取 = 0 ，2ii。

5 算法选择

前面已提到，隐式算法在求解非线性方程过程中存在不收敛的问题，故进行弹塑性动力时程分析宜采用显式算法。另外，显式算法不断更新变形后的结构形态，所以它可自动考虑几何非线性的问题。显式算法中的一个关键参数是稳定时间增量。若求解过程中，时间增量过大，并超过稳定时间增量，则可能导致计算结果不稳定。稳定时间增量可按下式进行估计：