动力时程分析中弹塑性刚度矩阵的提取方法

姜忻良，张海顺

(1. 天津大学建筑工程学院，天津 300072；2. 滨海土木工程结构与安全教育部重点实验室(天津大学)，天津 300072)

动力时程分析中弹塑性刚度矩阵的提取方法

姜忻良1,2，张海顺1,2

(1. 天津大学建筑工程学院，天津 300072；2. 滨海土木工程结构与安全教育部重点实验室(天津大学)，天津 300072)

利用ANSYS 二次开发工具UPFs编写成可执行程序文件，通过与MATLAB进行联立计算，对线弹性和非线性弹塑性本构模型的结构分别实现了特性矩阵的近似提取．对线性本构关系模型的各特性矩阵提取方法进行阐述推导和验证，求得的各阶频率和动态响应时程曲线与 ANSYS直接计算法吻合．对于非线性弹塑性本构模型，将非线性本构关系在全局过程中的小段时间内进行等效线性化处理后提取等效刚度矩阵，进行静动力分析后所得各种响应时程曲线与 ANSYS直接计算法所得曲线基本吻合．结果表明了分段等效线性化模型的合理性和ANSYS二次开发程序提取等效特性矩阵的可行性．

ANSYS二次开发；弹塑性刚度矩阵；分段线性化；等效弹性模量；稀疏矩阵；模态分析；动力时程分析

在采用商业软件进行各种结构动、静力分析时，研究人员根据不同目的而需要对刚度、质量、阻尼和荷载矩阵进行提取，以方便利用上述各类矩阵做后期分析与处理[1-3]．其中对于线弹性结构，研究人员可利用编程手段进行提取，方法较为简单[4-5]；但对于非线性弹塑性结构，上述各类矩阵的提取较为困难，目前大型商业软件一般无直接简便的提取方法．这给采用变化的结构特性矩阵来研究结构某些阶段或整个阶段的受力特点带来困难．

本文利用ANSYS二次开发工具手段来尝试进行提取编程，而其中刚度矩阵(包括单元刚度矩阵、原始结构刚度矩阵和结构刚度矩阵等)的提取是最为研究人员所关心的，本文主要对提取结构刚度矩阵进行尝试分析．对线性本构关系模型的刚度矩阵基于ANSYS二次开发提取方法进行阐述推导和验证，分别由特征值方程计算求得频率和由地震波动力时程分析求得响应曲线，并与 ANSYS直接计算法所得到的频率和响应曲线进行对比分析，以验证所二次开发的程序文件提取矩阵的可行性．

1 ANSYS二次开发工具

用于ANSYS二次开发的工具主要有4个，即APDL、UPFs(user programmable features)、UIDL和Tcl/Tk，使用以上工具可以建立新的材料模型，构建新的单元类型，参数化建模，优化分析，构建流程化的 ANSYS分析平台，建立符合用户专业需求的ANSYS用户界面等[6-8]．

本文主要依据APDL和UPFs两种工具，其中APDL应用比较普遍，不再赘述介绍，而UPFs是在ANSYS提供的 Fortran源代码的基础上，修改其用户可编程子程序和函数，从源代码层次上对ANSYS进行二次开发的工具．用户需要在响应的Fortran语言编译器的支持下，将编译修改后的源代码与ANSYS库相连而形成用户所需的ANSYS可执行文件来进行后续操作．

2 结构刚度矩阵提取方法

在ANSYS有限元程序中提取整体刚度矩阵方法主要是 HBMAT命令法和超单元(super-element)法．HBMAT命令是ANSYS提取整体刚度矩阵的直接内部提取方法，但该命令仅适用于线弹性分析，对非线性分析不适用，且该命令采用索引存储的稀疏矩阵并以Harwell-Boeing格式来存储刚度矩阵下三角的非零数据，并需要后期借用 MATLAB等程序编写命令处理来还原其矩阵形式．但是同时HBMAT命令法也有严重不足，首先形成的刚度矩阵是无序的，对于不同结构的单元网格，研究人员很难推算出矩阵的某一行某一列的值是与哪个节点自由度所相关的，且不能说明整体刚度是如何从单元刚度矩阵“对号入座”组装来的，而此过程却恰恰是研究人员所关心的内容．

超单元法同样仅可提取线弹性整体刚度矩阵，基本步骤是：创建有限元模型并施加约束，定义分析类型为子结构，定义输入何种矩阵，选择并定义所有节点为主自由度，求解并列出矩阵．所列出的整体刚度矩阵为全部元素(全矩阵)，按行的顺序分别列出各列元素数值，但其问题是当结构节点较多时，数据量非常庞大，且后期提取数据需要大量人工手动处理，非常繁琐导致使用不方便．

上述HBMAT命令法和超单元法在提取整体刚度矩阵的应用上均有不便之处，所以本文尝试基于Fortran语言的UPFs工具在ANSYS有限元程序中进行二次开发来提取矩阵．其核心思想就是通过编写外部用户程序直接从ANSYS子空间计算方法的模态分析结果的File.Full二进制文件中提取矩阵的各行各列非零值，然后编写 MATLAB程序使其有序地按照节点编号由小到大的顺序“对号入座”组装成为完整的矩阵形式．

ANSYS二次开发环境为Compaq Fortran 6.5，其中主文件为 Matrix-extraction.For，其余 Matrixoutput.F90用于矩阵输出，Binlib.Lin为 ANSYS提供库文件，Binlib.Dll为动态链接库文件．运行编译后而形成Matrix-output.Exe文件即可得到质量矩阵(mass matrix)和刚度矩阵(stiffness matrix)文件．整个分析的具体的流程如下．

(1) 线弹性本构模型ANSYS-MATLAB联立计算．对于线弹性本构模型的结构，首先提取模型的各个节点编号，进行ANSYS模态分析后按节点编号顺序在MATLAB中对号入座生成初始刚度矩阵和质量矩阵，并由此计算Rayleigh阻尼矩阵和荷载矩阵．

(2) 弹塑性本构模型ANSYS-MATLAB联立计算．对于非线性弹塑性本构模型的结构，每个时刻计算完毕后判断结构是否进入塑性阶段，若未进入塑性阶段则仍利用方法(1)的线弹性计算步骤；若已经进入塑性阶段，则首先提取此时的各个塑性单元的应力和应变值，通过由公式推导求出的等效弹性模量和等效剪切模量，并赋予该单元新的材料属性；重新进行ANSYS模态分析后，通过MATLAB程序提取其等效刚度矩阵、质量矩阵、阻尼矩阵和荷载矩阵；上述步骤均完成之后再进行动力时程分析从而结束一个时刻的计算．随后将其本构模型恢复至初始状态，利用 ANSYS的重启命令继续下一时刻的计算，如此往复进行．

3 线性模型的特性矩阵提取方法的验证

为检验上述ANSYS二次开发方法所提取的刚度矩阵和质量矩阵的正确性，有必要做一般性验证分析．对一简单弹性平面应变模型(见图 1)进行模态分析和天津人工地震波(见图 2)作用下的动力分析．

图1 网格模型Fig.1 Unit grid model

图2 天津人工地震波Fig.2 Tianjin artificial wave

模型采用平面应变 PLANE42单元，弹性模量E＝80,MPa，泊松比μ＝0.3，密度ρ＝1,750,kg/m3，尺寸为 5,m×12,m，底部固定．上述分析均采用ANSYS直接计算法和ANSYS-MATLAB联立法计算并进行了对比分析．对比分析内容包括频率、动态响应时程曲线和计算时间．

用ANSYS直接进行模态分析计算并提取所有模态频率，再利用 ANSYS二次开发程序所提取的刚度矩阵K和质量矩阵M在MATLAB中进行特征值方程运算，求得各阶频率进行对比，见表 1．从表 1中可以看出两种方法计算所得的各阶频率相同，并且采用二次开发程序所运行的计算时间大为减少，表明此 ANSYS二次开发提取矩阵运算的正确性和高效性．

动力时程分析中对整体结构施加天津人工地震波，提取其顶端中点的位移曲线待用，此过程由ANSYS直接计算；然后，将 ANSYS和 MATLAB两个程序联立，同样利用二次开发方法提取特性矩阵(刚度、质量、阻尼和荷载矩阵)，求得其顶端中点位移时程曲线并与 ANSYS的计算结果对比，如图3所示．由图 3可以看出，提取的特性矩阵所计算的顶端中点位移时程曲线与ANSYS直接计算的结果吻合；速度、加速度曲线为位移曲线对时间求导所得，故也吻合．此外，二次开发程序所运行的动力时程分析的计算时间也比ANSYS直接计算法的时间大幅度缩减，这也同样说明了上述矩阵提取方法的正确性和高效性．

表1 两种方法计算所得模态频率Tab.1 Modal frequency of two methods

图3 天津人工地震波矩阵验证Fig.3 Verification of matrix of Tianjin artificial wave

4 非线性等效刚度矩阵处理方法

第3节中的分析验证了ANSYS二次开发在线弹性模型下所提取的矩阵真实有效．若考虑结构在非线性塑性阶段提取其各个时刻的刚度矩阵，可以近似地认为非线性塑性本构关系是分段逐步递减的在小时间段 Δt内的线性本构关系的组合．因此，若要在动力时程分析的塑性阶段生成刚度矩阵，就需要将 Δt时间内视为线性本构关系．

对于进入非线性阶段的结构，在小时间段 Δt内，提取t和t+ Δt时刻的各个单元节点的x、y方向的正应力 σx、σy和剪应力τxy，以及与之对应的 x、y方向的正应变 εx、εy和剪应变 γxy．对于平面应变问题，依据如下公式进行推导．各单元的每一个节点的平面应变的物理方程为

t时刻

t＋ Δt时刻

则 Δt时段的应变差为

求解式(5)得到该单元节点的等效弹性模量和等效剪切模量，即

式(3)为各向同性的物理方程，若为了提取刚度矩阵而修改E、G则变成各项异性，则物理增量方程应为

式(4)可以简化为

注意在 ANSYS每个时刻 Step结束后调用MATLAB时，必须设置只有在该 Step调用的MATLAB程序完成后，才可以继续运行 ANSYS，即ANSYS在调用MATLAB程序未计算完毕时，必须等待而不能继续进行运算．这就需要两者同时运行并建立一个Flag文件，通过在两者中读其内容来判断对方是否在运行．两者若运行完一个 Step，改变 Flag，告诉对方自己当前运行结束，对方可以继续运行，否则必须等待．

5 非线性模型特性矩阵提取方法验证

第 4节推导出非线性塑性阶段分段等效刚度矩阵的计算提取方法，下面分别应用Pushover静力非线性分析和天津人工地震波的动力非线性分析方法来进行该方法的验证．

模型同样采用图 1所示的结构，其非线性本构关系采用理想弹塑性的 DP模型，黏聚力 C为20,kPa，内摩擦角和膨胀角均为 30°．在 Pushover静力非线性分析中，假定在某一荷载步时，某一单元的最高等效塑性应变增量超过该单元前一步等效塑性应变峰值的 15%时，认为此结构进入强非线性阶段，接近破坏而停止加载[9-10]．模型的 von Misis等效塑性应变云图如图4所示．

将ANSYS和MATLAB两个程序对接，在每一个荷载步结束后进行后处理分析．依据式(6)计算每个单元的等效弹性模量和等效剪切模量随后赋予该单元材料属性，继而进行模态分析并对整体结构提取其刚度矩阵与质量矩阵，再调用MATLAB程序通过特征值方程依次求出前3阶频率．全过程的前 3阶频率变化趋势如图 5所示，可以看出结构在进入塑性阶段后，各阶频率逐渐降低，说明其刚度在逐渐减小．

图4 模型的von Misis等效塑性应变云图(Pushover)Fig.4 von Misis equivalent plastic strain contour of the model(Pushover)

图5 前3阶频率的变化曲线(Pushover)Fig.5 Variation curves of the 1st three order frequencies (Pushover)

依据上述各个时刻的矩阵，利用 MATLAB通过 Hooke定律 Ku =F，求得其顶端中点位移时程曲线和底部最大塑性区节点位移曲线，与 ANSYS直接计算法的曲线进行对比，结果如图6所示．

图6 模型的位移曲线(Pushover)Fig.6 Displacement curves of the model(Pushover)

可以看出图 6中两种方法的位移曲线比较吻合，说明在非线性弹塑性阶段所提取的等效弹性模量和等效剪切模量较为准确，继而说明式(6)推导的在 Δt时间内非线性弹塑性本构模型可以转变为分段等效线性化模型进行处理，也验证了本文ANSYS二次开发程序的可行性．

图 7中为天津人工地震波动力分析的结果．由图 7可见在地震波作用下，结构底部两端均进入塑性阶段，与实际情况较符合．此外，利用 ANSYS直接计算法和ANSYS-MATLAB联立法得到的结构顶部中点位移、速度、加速度时程曲线见图 8．这些动力响应(位移、速度、加速度)时程曲线基本上吻合，说明采用等效弹性模量和等效剪切模量所提取的结构逐渐变化的刚度矩阵较为准确.

图7 模型的von Misis等效塑性区应变云图(地震波)Fig.7 von Misis equivalent plastic strain contour (earthquake wave)

图8 两种方法的顶部中点动力响应时程曲线对比Fig.8 Contrast of dynamic response time history curves of top midpoint between two methods

图 9为地震波作用下结构前 3阶频率的变化过程，可以看出结构在弹性阶段内频率保持不变，在进入塑性阶段时，结构的频率呈阶梯状降低．

图9 地震波作用下前3阶频率的曲线Fig.9 Curves of the 1st three order frequencies under earthquake action

表 2和表 3汇集了上述 Pushover静力非线性分析和天津人工地震波动力非线性分析的前 3阶频率在整个过程的某些时间点的变化情况．从Pushover静力非线性分析中可以看出，随着荷载的逐步增加，结构的前 3阶频率均逐步降低，而第 1阶频率降低最为突出．在天津人工波动力非线性分析中，在地震波达到峰值的时刻前 3阶频率的降低速度较大，当地震波趋于平稳的时候其前 3阶频率保持恒定，同样也是第1阶频率降低得最多．

表2 前3阶频率变化情况(Pushover)Tab.2 Variation of the 1st three order frequencies(Pushover)

表3 前3阶频率变化情况(地震波)Tab.3 Variation of the 1st three order frequencies(Earthquake wave)

6 结 语

采用分段等效线性化的处理方法，对非线性弹塑性本构模型在 Pushover静力方法与时程分析中进行各种矩阵的提取．结果表明在非线性弹塑性阶段所提取并采用的等效弹性模量和等效剪切模量较为准确，说明在 Δt时间段内将非线性弹塑性本构模型转变为分段等效线性化模型的合理性，同时也验证了本文 ANSYS二次开发程序的可行性．

本文方法为采用变化的特性矩阵对结构进入塑性阶段后深层次的分析研究提供了有效的方法．

［1］ 张令心，石 磊. 土-结构相互作用地震反应分析软件及其二次开发[J]. 地震工程与工程振动，2006，26(3)：225-227.

Zhang Lingxin，Shi Lei. Software of seismic soil structure interaction analysis and its re-developing[J]. Earthquake Engineering and Engineering Vibration，2006，26(3)：225-227(in Chinese).

［2］ 刘艳萍，杨新华，杨文兵. 预应力钢筋混凝土局部有限元分析的 ANSYS二次开发[J]. 华中科技大学学报(城市科学版)，2005，22(增)：87-90.

Liu Yanping，Yang Xinhua，Yang Wenbing. ANSYS secondary development for the finite element analysis of prestressed reinforced concrete structures[J]. Journal of Huazhong University of Science and Technology(Urban Science Edition)，2005，22(Suppl)：87-90(in Chinese).

［3］ 王一功，杨佑发. 针对场地地震反应分析的 ANSYS二次开发[J]. 地震工程与工程振动，2004，24(2)：42-45.

Wang Yigong，Yang Youfa. A secondary development of ANSYS for site earthquake response[J]. Earthquake Engineering and Engineering Vibration，2004，24(2)：42-45(in Chinese).

［4］ 刘光栋，王解君，何放龙. 空间梁单元的几何非线性刚度矩阵的分解形式[J]. 湖南大学学报，1992，19(1)：60-71.

Liu Guangdong，Wang Jiejun，He Fanglong. Resolvedformulation of geometrical nonlinear stiffness matrix for three-dimensional beam element[J]. Journal of Hunan University，1992，19(1)：60-71(in Chinese).

［5］ 刘齐茂，燕柳斌. 基于 Newmark 法敏度计算的刚架结构动力优化[J]. 工程力学，2010，27(3)：145-154.

Liu Qimao，Yan Liubin. Dynamic optimization of frame structures using sensitivity calculation based on newmark method[J]. Engineering Mechanics，2010，27(3)：145-154(in Chinese).

［6］ 李 妍，吴 斌，欧进萍. 弹塑性结构等效线性化方法的对比研究[J]. 工程抗震与加固改造，2005，27(1)：1-6.

Li Yan，Wu Bin，Ou Jinping. Comparison of equivalent linearization methods for inelastic structures[J]. Earthquake Resistant Engineering and Retrofitting，2005，27(1)：1-6(in Chinese).

［7］ 曲 哲，叶列平. 建筑结构弹塑性地震响应计算的等价线性化法研究[J]. 建筑结构学报，2010，31(9)：95-102.

Qu Zhe，Ye Lieping. Equivalent linear analysis in estimating nonlinear seismic responses of building structures[J]. Journal of Building Structures，2010，31(9)：95-102(in Chinese).

［8］ 程光煜，叶列平. 基于等效线性化方法弹塑性

SDOF系统能量谱的研究[J]. 建筑结构，2007，37(8)：74-77. Cheng Guangyu，Ye Lieping. Estimation of input energy spectra of inelastic SDOF systems with equivalent linear method[J]. Building Structures，2007，37(8)：74-77(in Chinese).

［9］ 白建方. 复杂场地土层地震反应分析的并行有限元方法[D]. 上海：同济大学土木工程学院，2007.

Bai Jianfang. Parallel Finite Element Method for Seismic Response Analysis of Irregular Site[D]. Shanghai：School of Civil Engineering，Tongji University，2007(in Chinese).

［10］ 张国栋. 土与结构相互作用体系随机地震反应分析[D]. 武汉：武汉大学土木建筑工程学院，2004.

Zhang Guodong. Stochastic Seismic Response Analysis for Soil-Structure Interaction System[D]. Wuhan：School of Civil Engineering，Wuhan University，2004(in Chinese).

(责任编辑：樊素英)

Extraction Method of Elastic-Plastic Stiffness Matrix in Dynamic Time History Analysis

Jiang Xinliang1,2，Zhang Haishun1,2
(1. School of Civil Engineering，Tianjin University，Tianjin 300072，China；2. Key Laboratory of Coastal Civil Engineering Structure and Safety of Ministry of Education (Tianjin University)，Tianjin 300072，China)

Executable program for feature matrices extraction was written by using ANSYS secondary development tools UPFs. With MATLAB simultaneous calculation, the feature matrices of the linear elastic and nonlinear elastic-plastic constitutive model were respectively extracted. For the linear constitutive model, the extraction method of those feature matrices were derived and affirmed. All order frequencies and response curves were coincided with ANSYS direct calculation method on the linear constitutive model. Within short time of the entire period，the nonlinear constitutive model was transformed through equivalent linearization for extracting the equivalent stiffness matrix. The static and dynamic response time history curves were basically coincided with the ANSYS direct calculation method. Results show that the piecewise equivalent linearization is rational and the ANSYS secondary development program is feasible and accurate.

ANSYS secondary development；elastic-plastic stiffness matrix；piecewise linearization；equivalent elastic modulus；sparse matrix；modal analysis；dynamic time history analysis

TU317.1

A

0493-2137(2015)04-0355-07

10.11784/tdxbz201309007

2013-09-02；

2014-03-13．

国家自然科学基金资助项目(51178308，51278335)．

姜忻良（1951— ），男，博士，教授，jiangxinliang@126.com．

张海顺，zhangyibiao_0216@163.com．