弹塑性动力时程分析若干问题分析与探讨

《弹塑性动力时程分析若干问题分析与探讨》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

弹塑性动力时程分析若干问题的分析与探讨张剑（深圳大学建筑设计研究院，深圳518000）摘要本文分析与探讨了弹塑性动力时程的主流分析软件（ABAQUS）关于混凝土及钢管混凝土的本构关系、特殊情况下单元的选取、高振型阻尼比的影响、混凝土在重力作用下的本构关系和算法选择等方面的问题，并提出了相应解决方法，经理论分析与实际工程验证，其方法是可行的。另外，本文对弹塑性动力时程分析软件的使用作出了建议。关键词弹塑性动力时程分析,混凝土塑性损伤模型，阻尼，显式算法。ResearchonIssuesofelastic-plasticdynamictime-historyanalysisZHANGJian(ShenzhenUniversityInstituteofArchitecturalDesign&Research,Shenzhen518000,China)Abstract::Thispaperanalysesandprobesintoseveralproblemsofelastic-plasticdynamictime-historyanalysisaboutconstitutiverelationofconcreteandsteeltube-restrainedConcrete,selectionofelementinspecialcase,influenceofhighermodesdampingratio,selectionofalgorithmundertheloadofgravitity,andputforwardtheactualsolution.Besides,thispaperofferessomesuggestionaboutusingthesoftwareofelastic-plasticdynamictime-historyanalysis.Keywords:elastic-plasticdynamictime-historyanalysis,concretedamagedplasticity,damping,explicitalgorithm。1．引言随着人们对建筑功能与美观等方面的要求不断增加及建筑师的创作水平日益提高，建筑结构的形态日趋复杂。另一方面，人们对结构的安全性与经济性的要求越来越高，导致结构设计难度越来越大，传统的一些简化分析计算手段已难以满足复杂结构的设计要求。此外，结构受力变形过程本来就是一个非线性动态过程，要实现结构的设计目标，就要正确把握结构的真实行为，也就必须依据结构真实受力变形机制进行分析。基于上述原因，人们致力发展结构非线性动力分析方法，弹塑性动力时程分析属其中的一类分析方法。2．弹塑性动力时程分析的目的结构的抗震设计是结构设计的一项基本内容，抗震设计的基本目标是“小震不坏，中震可修，大震不倒”。另外，根据不同建筑的安全需求与经济方面的许可程度，按照性能化目标的思想，抗震设计目标可在基本目标下进一步细化与提高。一般来说，在安全与经济双重目标要求下，结构在小震状态下，它处于弹性状态，变形也较小，此时采用线弹性方法分析内力与变形是可行的。结构在中震状态下，少部分构件已进入塑性状态且变形加大。此时，若仍然采用线弹性方法来分析，则存在较大误差。结构在大震状态下，部分构件已进入塑性状态，并产生大变形，其P-△效应加剧，几何非线性程度加大，故计算分析不能采用线弹性方法，也不宜采用静力弹塑性方法，而应采用弹塑性动力时程分析方法。所以说弹塑性动力时程分析方法是实现结构大震与中震下结构性能目标的基本分析方法。 3．弹塑性动力时程分析的基本原理与基本过程结构分析的主要目标是获取结构的位移场、应变场及应力场，三者之间具有密切的关系，故我们仅需获得结构位移场即可。通过离散化的方法，按粘性阻尼理论，可将结构的弹塑性动力学方程表达如下：MuCuKuF(t)（1）式中：u为节点位移向量，结构连续体的位移场可通过节点位移向量求得。M为质量矩阵，C为阻尼矩阵，K为刚度矩阵，F为外力向量函数，t为时间变量，地震作用时，F(t)Mu，其中u为地面运动加速度,即地震波。gg由于在外力作用下，结构可能具有几何非线性与本构非线性（弹塑性本构是非线性本构中的一种），结构的形态、刚度矩阵及阻尼矩阵不断变化，使得上述方程的求解非常复杂。就刚度矩阵而言，它是由单元刚度矩阵Ke组装而成。根据弹塑性力学，单元刚度矩阵可表达为：TKe[B][D][B]dv（2）e式中：[B]为几何矩阵，通过几何矩阵，由位移可求得应变。[D]为本构矩阵，通过本构矩阵，由应变可求得应力。由于非线性效应，[B]、[D]矩阵是不断变化的。对弹塑性问题而言，一旦知道任何时刻的几何矩阵、本构矩阵，通过积分点的数值积分，即可得到单元刚度矩阵。也就是说，各积分点无论是处于弹性或塑性状态，我们可以都得到对应时刻的单元刚度矩阵。再通过边界条件，即可逐步求解得到节点位移向量，进一步可求得任意一处的位移、应变及应力，实现分析的目标。在上述离散化的过程，最一般的单元是三维实体单元，其位移模式可以是线性的或者是二次的，视精度与效率的要求而定。在具体问题中，由于受力与变形机制的特殊性，导致位移场与应力场具有一些特殊性，合理利用这些特殊性，并作出相应的假定，可大大提高计算效率和精度。如采用直法线等假定形成板单元，采用平截面等假定形成梁单元等。对弹塑性动力方程的求解，一般可分为两种求解算法：即隐式与显式。隐式算法常采用Newmark法，但它需要求解全区域的联立方程，因此，它不但求解过程复杂，而且容易导致结果不收敛的情况。显式算法采用中心差分法，对动力学方程进行时间积分，由一个时间增量步的动力学条件下求解下一时间增量步的动力学条件，当时间增量充分小时，不会产生结果不收敛的情况，可获得问题的解答，因此显式算法特别适合弹塑性动力时程分析，它的基本过程如下：先将弹塑性动力方程改写如下：MuF(t)CuKu（3）记PF(t)，ICuKu则MuPI（4）第一步：节点变量求解1、动力平衡方程1u|(t)(M)(P|I)(t)(t)上式中M为对角矩阵，易于求解。2、对时间显式积分(t|t|(tt)(t)u|u|u|tt(t)(t)(t)222 u|ut|u|(tt)(t)(tt)t(t)2第二步：单元应力计算1、由节点位移，可求单元位移场，由此可求出应变场及应变率场，进一步可求出各处应变增量d。2、根据本构关系计算应力|f(|,d)(tt)(t)3、由ICuKu可求节点内力I|(tt)第三步：将时间t变为(tt)，返回第一步。由上述内容可知，地震作用下弹塑性动力时程分析的主要问题是材料本构关系的建构、单元选用、大变形的描述、地震波选取、阻尼选取及算法选择。4．ABAQUS弹塑性动力时程分析ABAQUS是国际上最先进的大型通用有限元分析软件之一，它的非线性功能达到世界领先水平，从而使其具有力学系统仿真的功能，它广泛应用于工程和科研各个领域。近年来，我国工程界将它成功应用于弹塑性动力时程分析，其分析结果得到业内专家们广泛认可。4.1材料本构关系ABAQUS自带丰富的本构关系模型，可描述混凝土、钢材、岩土、高分子材料等物质的应力与应变关系。另外，ABAQUS还提供用户材料接口程序UMAT及UVMAT，因此可使用户自定义材料本构关系。对建筑结构来说，主要涉及混凝土与钢材两种材料，钢材本构关系可采用二折线或三折线弹塑性本构关系，由于钢材质地均匀、性能稳定，其动力滞回模型也较为简单,下面重点描述混凝土的本构关系。ABAQUS软件中，混凝土本构关系模型有混凝土弥散开裂模型、混凝土开裂模型及混凝土塑性损伤模型。其中混凝土塑性损伤模型（ConcreteDamagedPlasticity）可描述混凝土受动力往返作用下的力学行为，故而广泛用于地震作用下的弹塑性动力时程分析。下面通过考察单轴特征荷载作用下的应力与应变关系曲线来把握混凝土塑性损伤模型的主要概念。单轴特征荷载作用下的应力与应变关系曲线Stress-straincurveundertheactionofsingleaxialcharacteristicload 如图所示，曲线OAC代表单轴单向受拉的应力-应变曲线，曲线OLK代表单轴单向受压的应力-应变曲线，曲线OABDEFHIJ代表单轴特征荷载作用下的应力-应变曲线。直线OA段代是单轴受拉的弹性阶段，当应变增大并超过A点所对应的应变时，混凝土受拉屈服，并进入受拉塑性区；当应变增量为负时，混凝土进入受拉卸载阶段，其应力-应变曲线近似视为直线，其斜率与混凝土在B点的损伤程度相关，而受拉损伤程度（Dt）与B点的应变相关，混凝土受拉塑性应变越大，混凝土受拉损伤程度越大，其卸载时的刚度越小，这一点与我们的概念分析及实验结果是一致的。另外，在卸载直线范围内，若应变增量又为正时，则应力与应变值所构成的相点沿曲线DBC运动。当相点越过D点进入受压状态时，混凝土受压时的刚度可得到一定的恢复。当应变增量继续为负时，相点沿压缩曲线DE运动，当应变绝对值小于直线段最大应变绝对值时，混凝土处于受压弹性状态；否则，混凝土受压屈服，进入受压塑性状态。此时，若应变增量为正时，混凝土进入受压卸载阶段，相点沿FH直线段运动，其直线斜率与F点的受压损伤程度（Dc）相关。当应变增量继续为正时，相点越过H点进入受拉状态，并沿直线HI运动，此直线的斜率不断变小，不但与受拉损伤程度（Dt）相关，而且还与受压损伤程度（Dc）相关。当应变增量继续为正时，混凝土越过I点再次屈服进入受拉塑性状态。在上述过程中，值得指出的是，由于地震作用的随机性，相点随时有可能作反向运动；随着时间的推移，混凝土的受拉与受压损伤程度是不断增大的。上述内容是针对单轴（即一维）情况，对二维、三维情况，塑性损伤模型根据塑性势流动理论及单轴的应力-应变关系由塑性应变增量求出应力增量及应力。对钢管混凝土应考虑钢管的套箍效应对混凝土参数的影响,参考Mander约束混凝土模型，可认为具有套箍系数为SITA的约束混凝土的抗压强度标准值FC与混凝土峰值压应变EC将按下面关系得以提高。IF(SITA.LE.1.235)THENFC=FC*(1+SITA);EC=EC*(1+5*SITA)ELSEFC=FC*(1+SQRT(SITA));EC=EC*(1+5*SQRT(SITA))ENDIF梁单元受拉损伤时，吸收与抵抗剪力的能力降低，故在编制考虑剪切效应三维梁单元显式算法的VUMAT弹塑性材料用户程序时，应考虑梁和柱受拉损伤对梁和柱抗剪模量的不利影响。另外，对同一强度等级的混凝土来说，梁和柱单元通过用户程序定义的混凝土参数须与壳和块单元直接定义的混凝土参数一致。4.2单元选用从一般意义上来说，任何结构可采用三维实体单元进行分析，但也因此会带来收敛速度慢、计算效率低、后处理工作量大等问题，所以采用三维实体单元分析一般用于节点分析或结构的分析研究等方面。ABAQUS软件中的三维单元可用于隐式与显式分析，对钢筋混凝土结构来说，三维实体元可使用上述所有混凝土本构关系模型，还可采用rebar方式或embedded方式模拟混凝土中的钢筋。ABAQUS的壳单元可用来描述楼板、剪力墙及无扭转效应的梁(如连梁)，它亦可用于隐式与显式分析，并可采用各种混凝土的本构模型。可采用rebarlayer方式在层中定义钢筋。ABAQUS中考虑剪切效应的三维梁单元的使用受到一定限制，如它不能直接采用混凝土塑性损伤模型，但我们可运用VUMAT接口程序的方法来弥补这方面的不足。对钢筋混凝土结构来说，钢筋混凝土梁与柱由混凝土梁单元与钢筋梁单元叠合而成，钢筋梁单元截面为箱形截面，其几何尺寸由梁与柱的实际配筋量求得。在弹塑性动力时程分析中，ABAQUS无需定义塑性区，只需基于材料定义其塑性行为特性即可，这样，单元的刚度由单元积分求得，结构的塑性区及其塑性损伤程度由外部作用、本构关系及其力学规律来确定，因此它能较真实反映结构实际受力变形规律。上述的梁单元与壳单元也存在一定的缺陷，即梁单元的塑性效应不能考虑剪力的影响，壳单元的塑性效应不能考虑壳法向剪力的影响,所以只有当上述剪力较小时，其误差是可接受的。对一般具有剪力墙的高层建筑来说，剪力墙面内吸收大部分剪力，所以上述做法可满足精度的要求。但对一些特殊的情 况须作特殊处理，如结构中存在短柱，短柱吸收的剪力较大，易产生脆性的剪切破坏，故采用上述梁单元来模拟就不合适了,所以说，如单方向剪力较大的柱则采用壳单元模拟，双向剪力都较大的柱则应采用块单元模拟，块单元与梁和柱的相连接节点可用刚性单元连接。又如墙的厚度偏大时，采用壳单元模拟也不合适，而应采用块单元模拟。另外，连梁的竖向剪力较大，则应采用壳单元来模拟。4.3地震波选取地震波的选取首先应符合规范要求，除此之外，地震波的频谱特性应与场地的构造及其参数相协调。若结构扭转不规则时，计算分析须考虑双向水平地震。若结构存在跨度大于24m的结构或跨度大于8m的转换结构或悬挑长度大于5m的悬挑结构或高度超100m，则宜考虑竖向地震作用。双向水平地震和竖向地震同时作用，同一组地震波的两个水平分量和竖向分量的加速度比值可取1:0.85:0.65。对地震波来说，尽管峰值相同，但持续时间与频谱特性不一样，结构的效应会相差较大，另外，我们很难知道地震发生时的地震波是什么样，所以很难选择到恰当的地震波。一般来说，建议重点分析场地波，它毕竟与结构所在的场地情况和地层结构相关。4.4阻尼选取在结构动力学中，阻尼的选取对计算结果较为敏感，故应非常谨慎对待。一般情况可采用瑞利（Rayleigh）阻尼模式来定义阻尼，即（1）式中的阻尼矩阵C表达如下：CMK（5）式中：为质量阻尼系数，为刚度阻尼系数。与是难以直接确定的，但可根据它们与振型阻尼比的关系来间接确定。对多自由度力学系统，有如下关系：(/)/2（6）iii式中：i是系统园频率为i时的阻尼比,其值可根据特定的材料在自由振动的情况下振动幅值的衰减情况测得。对混凝土材料来说，可取0.05,对钢材料来说，可取0.02,而且可认为各种频率下的振型阻尼比是相同的。假定对混凝土结构来说，各处混凝土的与系数均相同，也即阻尼矩阵C仅由两个参数来决定，此时根据（6）式可知，阻尼矩阵C仅能保证两个振型的阻尼比为0.05,难以保证其它振型的阻尼比为0.05，此时也可通过(6)式求得系统园频率为j时的计算阻尼比j，如果j<0.05,则导致计算结果对振型j阻尼估计不足，计算效应偏大，可能导致设计偏保守；如果仅采用质量阻尼系数，并按结构基本频率来计算质量阻尼系数,即2,因此对高振型的阻尼比,有/(2),i1ii也即对高振型阻尼比估计偏小，导致结构效应偏大，设计偏保守；如果j>0.05,则导致计算结果对振型j阻尼估计过大，效应偏小，可能导致设计偏不安全，所以如何正确构造阻尼矩阵是一个至关重要的问题。对这个问题首先要建立这样的目标：合适的阻尼矩阵一是使方程（1）能解耦，二是使得结构各振型的阻尼比为指定值。令uZ(7)上式中φ为（1）式无阻尼自由振动方程的nxn阶正则化后的振型矩阵，Z为n阶广义自由度向量。T将（7）式代入（1）式并左乘,则有： TTTTMZCZKZF(t)(8)T由于的正交性，即有：M=I=diag([1,1,…,1])（9a）上式中,I为nxn阶单位矩阵,diag()为对角阵函数。TT且K为对角矩阵,因此要使方程(8)可解耦,C必为对角阵,亦即：TCdiag([c,c,...,c])(9)12n(8)式变为如下方程：ZcZkZTF(t)(10)iiiiTT式中：cC，kK，i1,2,,niiiiii令c/2,即第i振型的阻尼比。iii2T(10)式变为Z2ZZF(t)(11)iiiiii故有c2(13)iii由上可知,阻尼矩阵按如下方式构造必能实现预定的两个目标。TCM*diag([c,c,...,c])*M(14)12n上式中，符号*(非上标时)为矩阵运算中的乘号。T作一个验证，对(14)式分别左乘和右乘,并根据(9a)式，则有：TTTCM*diag([c,c,...,c])*Mdiag([c,c,...,c])(15)12n12n由(13)与（15）式，可知阻尼矩阵的构造实现了两个预定的目标。出于对计算效率的考虑，振型并不须取满，即取前m(m