第48卷 第4期 2015年4月
天津大学学报(自然科学与工程技术版)
Journal of Tianjin University (Science and Technology )
V ol.48 No.4Apr. 2015
收稿日期：2013-09-02；修回日期：2014-03-13．
基金项目：国家自然科学基金资助项目(51178308，51278335)．
作者简介：姜忻良（1951— ），男，博士，教授，*********************． 通讯作者：张海顺，************************．
DOI:10.11784/tdxbz201309007
姜忻良1, 2，张海顺1, 2
(1. 天津大学建筑工程学院，天津 300072；
2. 滨海土木工程结构与安全教育部重点实验室(天津大学)，天津 300072) 摘 要：利用ANSYS 二次开发工具UPFs 编写成可执行程序文件，通过与MATLAB 进行联立计算，对线弹性和非线性弹塑性本构模型的结构分别实现了特性矩阵的近似提取．对线性本构关系模型的各特性矩阵提取方法进行阐述推导和验证，求得的各阶频率和动态响应时程曲线与ANSYS 直接计算法吻合．对于非线性弹塑性本构模型，将非线性本构关系在全局过程中的小段时间内进行等效线性化处理后提取等效刚度矩阵，进行静动力分析后所得各种响应时程曲线与ANSYS 直接计算法所得曲线基本吻合．结果表明了分段等效线性化模型的合理性和ANSYS 二次开发程序提取等效特性矩阵的可行性． 关键词：ANSYS 二次开发；弹塑性刚度矩阵；分段线性化；等效弹性模量；稀疏矩阵；模态分析；动力时程分析 中图分类号：TU317.1 文献标志码：A 文章编号：0493-2137(2015)04-0355-07
Extraction Method of Elastic -Plastic Stiffness Matrix in
Dynamic Time History Analysis
Jiang Xinliang 1, 2，Zhang Haishun 1, 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 )
Abstract ：Executable program for feature matrices extraction was written by using ANSYS secondary develop-
ment tools UPFs. With MATLAB simultaneous calculation, the feature matrices of the linear elastic and nonlinear elastic-p lastic constitutive model were resp ectively extracted. For the linear constitutive model, the extraction method of those feature matrices were derived and affirmed. All order frequencies and response curves were coin-cided with ANSYS direct calculation method on the linear constitutive model. Within short time of the entire pe-riod ，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 p iecewise equivalent linearization is rational and the ANSYS secondary development program is feasible and accurate.
Keywords ：ANSYS secondary develop ment ；
elastic-p lastic stiffness matrix ；p iecewise linearization ；equivalent elastic modulus ；sparse matrix ；modal analysis ；dynamic time history analysis
在采用商业软件进行各种结构动、静力分析时，研究人员根据不同目的而需要对刚度、质量、阻尼和荷载矩阵进行提取，以方便利用上述各类矩阵做后期分析与处理[1-3]．其中对于线弹性结构，研究人员可利用编程手段进行提取，方法较为简单[4-5]；
但对于非线性弹塑性结构，上述各类矩阵的提取较为困难，目前大型商业软件一般无直接简便的提取方法．这给采用变化的结构特性矩阵来研究结构某些阶段或整个阶段的受力特点带来困难．
本文利用ANSYS 二次开发工具手段来尝试进
·356·天津大学学报(自然科学与工程技术版)第48卷 第4期
行提取编程，而其中刚度矩阵(包括单元刚度矩阵、原始结构刚度矩阵和结构刚度矩阵等)的提取是最为研究人员所关心的，本文主要对提取结构刚度矩阵进行尝试分析．对线性本构关系模型的刚度矩阵基于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，其余Matrix-output.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二次开发方法所提取的刚
2015年4月姜忻良等：动力时程分析中弹塑性刚度矩阵的提取方法 ·357·度矩阵和质量矩阵的正确性，有必要做一般性验证
分析．对一简单弹性平面应变模型(见图1)进行模
态分析和天津人工地震波(见图2)
作用下的动力
分析．
图1网格模型
Fig.1
Unit grid model
图2天津人工地震波
Fig.2Tianjin artificial wave
模型采用平面应变PLANE42单元，弹性模量E＝80MPa，泊松比μ＝0.3，密度ρ＝1750kg/m3，尺寸为5m×12m，底部固定．上述分析均采用ANSYS直接计算法和ANSYS-MATLAB联立法计算并进行了对比分析．对比分析内容包括频率、动态响应时程曲线和计算时间．
用ANSYS直接进行模态分析计算并提取所有模态频率，再利用ANSYS二次开发程序所提取的刚度矩
阵K和质量矩阵M在MATLAB中进行特征值方程运算，求得各阶频率进行对比，见表1．从表1中可以看出两种方法计算所得的各阶频率相同，并且采用二次开发程序所运行的计算时间大为减少，表明此ANSYS二次开发提取矩阵运算的正确性和高效性．
动力时程分析中对整体结构施加天津人工地震波，提取其顶端中点的位移曲线待用，此过程由ANSYS直接计算；然后，将ANSYS和MATLAB 两个程序联立，同样利用二次开发方法提取特性矩阵(刚度、质量、阻尼和荷载矩阵)，求得其顶端中点位移时程曲线并与ANSYS的计算结果对比，如图3所示．由图3可以看出，提取的特性矩阵所计算的顶端中点位移时程曲线与ANSYS直接计算的结果吻合；速度、加速度曲线为位移曲线对时间求导所得，故也吻合．此外，二次开发程序所运行的动力
时程分析的计算时间也比ANSYS直接计算法的时
间大幅度缩减，这也同样说明了上述矩阵提取方法
的正确性和高效性．
表1两种方法计算所得模态频率
Tab.1Modal frequency of two methods
频率/Hz
模态阶数
ANSYS直接计算法 ANSYS-MATLAB联立法
1 1.119 1.119
2 4.59
3 4.593
3 4.707 4.707
660 566.051 566.051
计算时间/s
4.228 0.764
图3天津人工地震波矩阵验证
Fig.3Verification 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时刻
2
xy
2
1()()
()=
1
()
1()()
()=
1
()
()
()
()
王永庆的球童
=
x y
x
y
y
x
建材营销论坛y x
t t
t
E t
t t
t
E t
t
t
G
t
μ
μσσ
ε
μ
μ
μσσ
τ
ε
γ
μ
⎧−⎡⎤
−
⎪⎢⎥
−
⎣⎦
⎪
⎪−⎡⎤
⎪−
⎨⎢⎥
−
⎣⎦
⎪
⎪
⎪
⎪⎩
(1)
·358· 天津大学学报(自然科学与工程技术版) 第48卷 第4期
t ＋t Δ时刻
221(+)(+)(+)=1(+)1(+)(+)(+)=1(+)(+)(+)(+)=x y x y x y xy
xy t t t t t t E t t t t t t t t E t t t t t G t t t μμσσεμμμσσεγμτ⎧−⎡⎤
Δ−ΔΔ⎪⎢⎥−Δ⎣⎦⎪
⎪−⎡⎤⎪
Δ−ΔΔ⎨⎢⎥−Δ⎣⎦⎪
⎪Δ⎪ΔΔ⎪⎩
(2)
则t Δ时段的应变差为
{
2222
21(+)(+)(+)()=1(+)1()()1()1(+)()()(+)()11(+)(+)(+)()=1(+)1()()x y x x x y x x y y y x y y y t t t t t t t E t t t t E t t t t E t t t t t t t t t t t E t t t E t μμσσεεμμμσσμμσσμσσμμμσσεεμμσ−⎡⎤
Δ−ΔΔ−−⎢⎥−Δ⎣⎦−⎡⎤
−=⎢⎥−⎣⎦−Δ−−
Δ⎫Δ−⎡⎤⎬⎣⎦−⎭
−⎧Δ−Δ−
Δ−⎨−Δ⎩−−（）
{
[]2()11(+)()()(+)()1(+)()(+)()==(+)()(+)()()x y y x x xy xy xy xy xy xy t t t t E t t t t t t t t t t G t t G t t t t G t μσμμσσμ
σσμττγγττ⎧⎪⎪
⎪⎪⎪
⎪⎪⎪⎪
⎪⎪⎪⎪⎪⎪
⎨
⎫⎡⎤=⎬⎢⎥−⎣⎦⎭
−Δ−−Δ⎫Δ−⎪⎪⎬⎪−⎭ΔΔ−−
⎪⎪
⎪⎪⎪⎪
⎪⎪⎪
⎪⎪⎪⎩
ΔΔ−Δ（
）
(3)
式(3)为各向同性的物理方程，若为了提取刚度矩阵而修改E 、G 则变成各项异性，则物理增量方程应为 22(+)()(+)()=(1)()(+)()()1(+)()
(+)()=(1)()(+)()()1(+)()(+)()=()
x x x x
x y y y y y y y y x x x xy xy xy
xy t t t t t t E t t t t E t t t t t t t E t t t t E t t t t t t t G t σσεεμσσμμσσεεμσσμμττγγ⎧Δ−⎧−Δ−−⎨Δ⎩Δ−⎫⎡⎤⎪⎬
黄昏到寺蝙蝠飞⎢⎥Δ−⎪⎣⎦⎭
Δ−⎧⎪−Δ−−⎨⎨
Δ⎪⎩Δ−⎫⎡⎤⎬⎢⎥Δ−⎣⎦⎭
Δ−Δ−Δ⎪⎪
⎪
⎪⎪⎪⎪⎪
⎪
⎪
⎪⎪
⎪
⎪⎩
(4)式(4)可以简化为
2
2
()()()=(1)()1()()()()=(1)()1()()()=()y x x x y y x y y x xy
xy t t t E t E t t t t E t E t t t G t σσμεμμσσμεμμτγ⎧ΔΔ⎡⎤ΔΔ−ΔΔ−⎪⎢Δ−Δ⎪⎣⎦⎪
ΔΔ⎡⎤ΔΔ⎪−ΔΔ−⎨⎢⎥Δ−Δ⎣⎦⎪
⎪ΔΔ⎪ΔΔ⎪
Δ⎩
(5)
求解式(5)得到该单元节点的等效弹性模量和
等效剪切模量，即
2222
2222
()1(1)()=
()
()+111()1(1)()=
()()
+111()()=
()x x y x y y y x xy xy t E t t t t E t t t t G t t μσμεεμμμμμσμεεμμμμτγ⎧⎡⎤
ΔΔ−⎪⎢⎥
−⎣⎦⎪ΔΔΔ⎪ΔΔ⎪
−−−⎪⎪⎡⎤⎪
ΔΔ−⎨⎢⎥−⎣⎦⎪Δ⎪ΔΔΔΔ⎪
−−−⎪⎪ΔΔΔ⎪ΔΔ⎪⎩
厦门px(6)
若对平面应力问题进行等效弹性模量的推导，
则可以将2()/(1)E t μΔ−换成)E t Δ，(1)μμ−换成μ即可．依据上述推导，对于ANSYS 每一时刻分析
完成后，判断结构是否进入塑性阶段．若没有进入塑性阶段，则刚度矩阵仍保持不变；若已进入塑性阶段，则在t Δ时段内对各单元的每个节点提取其相应的应力和应变进行式(6)的计算，再取其平均值作为该单元的等效弹性模量()E t Δ和等效剪切模量
()G t Δ．当)E t Δ和()G t Δ小于初始弹性模量E 和G
时，表明结构进入非线性阶段．用每个单元的()E t Δ和()G t Δ代替E 和G 来进行分段等效线性化处理．在此分段等效线性化的本构关系下重新进行模态分析，然后如第2节所述从File.Full 文件中提取t Δ时段的等效刚度矩阵来调用MATLAB 程序进行后续的方程计算．计算完毕后，将材性变化的单元重新赋予初始E 和G ．随后需要利用ANSYS 的重启动命令继续进行下一时刻的时程分析．该命令功能将保留上一时刻的计算结果而继续往下计算，避免必须从初始零时刻重新进行计算．
注意在ANSYS 每个时刻Step 结束后调用MATLAB 时，必须设置只有在该Step 调用的MATLAB 程序完成后，才可以继续运行ANSYS ，即ANSYS 在调用MATLAB 程序未计算完毕时，必须等待而不能继续进行运算．这就需要两者同时运行并建立一个Flag 文件，通过在两者中读其内容来判断对方是否在运行．两者若运行完一个Step ，改
2015年4月姜忻良等：动力时程分析中弹塑性刚度矩阵的提取方法 ·359·变Flag，告诉对方自己当前运行结束，对方可以继
续运行，否则必须等待．
5 非线性模型特性矩阵提取方法验证
第4节推导出非线性塑性阶段分段等效刚度
矩阵的计算提取方法，下面分别应用Pushover静力
非线性分析和天津人工地震波的动力非线性分析
方法来进行该方法的验证．
模型同样采用图1所示的结构，其非线性本构
关系采用理想弹塑性的DP模型，黏聚力C为
20kPa，内摩擦角和膨胀角均为30°．在Pushover
静力非线性分析中，假定在某一荷载步时，某一单
元的最高等效塑性应变增量超过该单元前一步等
效塑性应变峰值的15%时，认为此结构进入强非线
性阶段，接近破坏而停止加载[9-10]．模型的von
Misis等效塑性应变云图如图4所示．
将ANSYS和MATLAB两个程序对接，在每一
个荷载步结束后进行后处理分析．依据式(6)计算
每个单元的等效弹性模量)
E tΔ和等效剪切模量
)
G tΔ，随后赋予该单元材料属性，继而进行模态分
dtfd
析并对整体结构提取其刚度矩阵与质量矩阵，再调
用MATLAB程序通过特征值方程依次求出前3阶
频率．全过程的前3阶频率变化趋势如图5
所示，
图4模型的von Misis等效塑性应变云图(Pushover) Fig.4von Misis equivalent plastic strain contour of the model(Pushover)
图5前3阶频率的变化曲线(Pushover)
Fig.5Variation curves of the 1st three order frequen-cies (Pushover)可以看出结构在进入塑性阶段后，各阶频率逐渐降低，说明其刚度在逐渐减小．
依据上述各个时刻的矩阵，利用MATLAB通过Hooke定律=
Ku F，求得其顶端中点位移时程曲线和底部最大塑性区节点位移曲线，与ANSYS 直接计算法的曲线进行对比，结果如图6所示．
（a）顶端中点位移
（b）底部最大塑性区节点位移
图6模型的位移曲线(Pushover)
Fig.6Displacement curves of the model(Pushover)
可以看出图6中两种方法的位移曲线比较吻合，说明在非线性弹塑性阶段所提取的等效弹性模量)
E tΔ和等效剪切模量)
G tΔ较为准确，继而说明式(6)推导的在tΔ时间内非线性弹塑性本构模型可以转变为分段等效线性化模型进行处理，也验证了本文ANSYS二次开发程序的可行性．
图7中为天津人工地震波动力分析的结果．由图7可见在地震波作用下，结构底部两端均进入塑性阶段，与实际情况较符合．此外，利用
ANSYS直
图7模型的von Misis等效塑性区应变云图(地震波) Fig.7von Misis equivalent plastic strain contour (earthquake wave)
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议