弹塑性矩阵推导
弹塑性矩阵推导(全面版)资料
弹塑性矩阵推导（全面版）资料考虑材料的塑性，其增量形式的本构关系可表达为p d σ=(D-D )d ε (1)式(1)中，D 为弹性矩阵，p D 为塑性矩阵。
弹性矩阵D 的形式为422000333242000333224000333000000000000K G K GK G K G K GK G K GK GK G G G G ⎡⎤+--⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥=⎢⎥--+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦D (2) 体积模量3(12)E K μ=-，剪切模量2(1)EG μ=+。
在应变空间内，塑性矩阵可表达为1()T f fA ∂∂=∂∂p D D D σσ(3) 式中，()()T T p f f f fA B ∂∂∂∂=--∂∂∂∂D D σσσσ(4) f为屈服函数；p σ为塑性应力，p p =σD ε；1/2(()())T p T p T p f f f f B f f f ωθε∂∂⎧⎪∂∂⎪∂∂⎪'=⎨∂∂⎪∂∂∂⎪⎪∂∂∂⎩σσI σσσ(5) [111000]T '=I ；p ω为塑性功；p θ为塑性体应变；p ε为等效塑性应变；κ为反映加载历史的参数。
对于Drucker-Prager 模型，其屈服条件为10f I α== (6)1x y z I σσσ=++，22222221()2x yz xy yz zx J S S S S S S =+++++，α为材料常数。
f α∂''=∂I σ (7) 222Txyzxyyzzx S S S S S S '⎡⎤=⎣⎦S(8)()3f K αα∂'''==+∂DD I I σ(9) 2()()(3)9T T T f f A K K G ααα∂∂'''===+∂∂D I I σσ (10) 1()T f fA ∂∂=∂∂p D D D σσ22222222112123113211212311321121231131121121121(3)(3)999(9)(T TT T T TK K K GK G K G K G J m mn ml S m S m S m mn n nl S n S n S n ml nl l S l S l S l S m S n S l ααααααβββββββββββββ''=+++''''=++++=p D I I I I S SS211211222311221321231231231122231231132232113113113113212113223113)()()S S S S S S m S n S lS S S S S S m S n S l S S S S S ββββββββββββββββββββ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⋅⋅⎢⎥⎢⎥⋅⋅⎢⎥⋅⋅⎢⎥⎣⎦令43p K G =+，23q K G =- 弹塑性矩阵可表达为2112123113211212311321121231132112112112112112223112213123123123112223()(p m q mn q mlS mS m S m q mnp n q nl S n S n S nq ml q nl p l S l S l S l S m S n S l G S S S S S S m S n S lS S G βββββββββββββββββββββββ------------------=-=-----⋅-⋅----⋅-ep p D D D 21231132232113113113113212113223113)()S S S S m S n S lS S S S G S ββββββββββ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⋅⎢⎥----⋅-⋅-⎢⎥⎣⎦令1β=，2β=1112m S ββ=+，1222n S ββ=+，大型矩阵系统◆ELL系列ELL**V##(视频矩阵)ELL**A##(音频矩阵)ELL**RV##(环通视频矩阵)通道选择输入：8XN+64 最大256路输出：4XN+8 最大32路产品介绍：ELL系列视、音频主机采用大规模音视频切换专用芯片作为音频视频切换矩阵电路的多路多通道切换设备，溶合了先进的矩阵切换技术和计算机技术，可以给用户提供卓越的整体性能。
建筑弹塑性分析PUSHOVER
2.需求谱法
结构抗震性能需求谱是在给定地震作用下， 不同周期结构的承载力和位移响应的需求 值。
先将能力曲线转化为A-D格式，能力谱曲线
将不同的周期结构的加速度响应需求Sa和位
移响应需求Sd也在A-D坐标系下给出，由此得
到的Sa-Sd关系曲线即为需求谱。对于弹性结
构，弹性谱加速度需求Sa可以采用地震弹性
其中 Dntqnt/，n D表n 示t 一个对应原结构
第n阶振型的单自由度体系在地震作用 下u g ( t ) 的位移响应，圆频率和阻尼比分别为 和 n 。
从而可n 求得结构第n阶振型的位移，内力，层
间位移等。
对前N阶振型都采用上述方法求算其最大响应 量，并采用某种方法进行组合（SASS法或 CQC法）—振型分解反应谱法。
Fass
T
ass
fs(D,signD)
aTssm ；对于地震响应由结构振型
向 量量成正控a s 比s制a s的s的荷弹载塑进性行结推构覆，，仍即采：用振型sa向ss mass
得到
Fass
Vb Mass
uroof
,DБайду номын сангаасass
roof ass
u u V
V
b
基底剪力， r o o顶f 点位移。 — r o 的o f 关系曲线称为
b
“结构的能力曲线”。或“推覆曲线”
为便于评价结构抗震性能是否达到要求，还
可以按照单阶振型反应谱法将推覆曲线上
各店的承载力和位移转化为谱加速度与谱 位移的关系曲线，得到结构的能力谱曲线，
即 S a S格d 式能力谱曲线。
Sa
Vb M
,
Sd
uroof
roof
静力弹塑性分析(Pushover分析)两种方法剖析
静力弹塑性分析(Pushover 分析)■ 简介Pushover 分析是考虑构件的材料非线性特点，分析构件进入弹塑性状态直至到达极限状态时结构响应的方法。
Pushover 分析是最近在地震研究及耐震设计中经常采用的基于性能的耐震设计(Performance-Based Seismic Design, PBSD)方法中最具代表性的分析方法。
所谓基于性能的耐震设计就是由用户及设计人员设定结构的目标性能(target performance)，并使结构设计能满足该目标性能的方法。
Pushover 分析前要经过一般设计方法先进行耐震设计使结构满足小震不坏、中震可修的规X 要求，然后再通过pushover 分析评价结构在大震作用下是否能满足预先设定的目标性能。
计算等效地震静力荷载一般采用如图2.24所示的方法。
该方法是通过反应修正系数(R)将设计荷载降低并使结构能承受该荷载的方法。
在这里使用反应修正系数的原因是为了考虑结构进入弹塑性阶段时吸收地震能量的能力，即考虑结构具有的延性使结构超过弹性极限后还可以承受较大的塑性变形，所以设计时的地震作用就可以比对应的弹性结构折减很多，设计将会更经济。
目前我国的抗震规X 中的反应谱分析方法中的小震影响系数曲线就是反应了这种设计思想。
这样的设计方法可以说是基于荷载的设计(force-based design)方法。
一般来说结构刚度越大采用的修正系数R 越大，一般在1~10之间。
但是这种基于荷载与抗力的比较进行的设计无法预测结构实际的地震响应，也无法从各构件的抗力推测出整体结构的耐震能力，设计人员在设计完成后对结构的耐震性能的把握也是模糊的。
基于性能的耐震设计中可由开发商或设计人员预先设定目标性能，即在预想的地震作用下事先设定结构的破坏程度或者耗能能力，并使结构设计满足该性能目标。
结构的耗能能力与结构的变形能力相关，所以要预测到结构的变形发展情况。
所以基于性能的耐震设计经常通过评价结构的变形来实现，所以也可称为基于位移的设计(displacement-based design)。
弹塑性_塑性力学基本方程和解法
在加载过程中物体各点处的偏应力分量 sij 保持比例不变。在工程允许精度下，也可推
广应用于稍为偏离简单加载的情况。
以上各种理论中涉及的一些假设，例如：塑性应变偏量的增在单一的函数关系等假设，都得到了常用金属材
料大量试验的验证。
z 强化规律 对于理想弹塑性材料，材料一旦屈服，其应力状态点在主应力空间中就落在屈服
变形， Hα 也不变，于是
∂f ∂σ ij
除等向强化外，有些强化材料表现为随动强化（图 7.7b），即，在强化过程中，屈
服面的大小和形状保持不变，只随塑性变形的发展而在应力空间中平移。还有些材料
在强化过程中随动强化与等向强化同时发生，称为混合强化。
由于在应力和强化参数空间中，表示应力状态的应力点只可能位于后继屈服面
（或加载面）上或其内，不可能位于曲面之外，若加载面是一个正则曲面，则有
⎯2⎯
研究生学位课弹塑性力学电子讲义
姚振汉
⎧ε = 0 ⎨⎩σ = σ s
当 σ <σs 当 ε >0
(2)
图 7.5 理想弹塑性和刚塑性
当考虑材料强化性质时，可在理想弹塑性模型的基础上加以改进，采用线性强化 弹塑性模型来近似：
⎧σ = Eε
⎨⎩σ = σ s +E1 (ε − εs )
当 ε ≤εs 当 ε >εs
(5)
⎯3⎯
第七章 塑性力学的基本方程与解法
其中 k 可由单向拉伸或其它材料试验测得的σ s 确定， k = σ s 2 。当不能确定主应力的 排序时，在以三个主应力为坐标轴的应力空间中，由特雷斯卡条件所包围的弹性状态 的应力空间为
σ1 −σ 2 ≤ 2k, σ 2 −σ 3 ≤ 2k, σ 3 −σ1 ≤ 2k
第二讲 弹塑性理论知识
J1 0 1 2 2 2 J 2 ( 1 2 3 1 2 2 3 3 1 ) 3 1 2 2 2 J 3 S1S 2 S 3 ( S1 S 2 S 3 ) 3
式中: J —偏应力第一不变量； 量； J —偏应力第二不变量
~
3
1 ~
2 ~
3 ~
我们 55 3
的直线为静水压力轴，在静水压力轴上每一点的
，球张量的几何解释是静水压力轴上的 一个分矢量。
1 2 3 P0
(3) 平面 过原点且垂直于静水压力轴的平面称为 平面。 它的方程是： 0 平面的法线矢量： n 1 (e e e )
因为主应力和坐标系的选择无关（即用主平面上的主 应力描述一点的应力状态不随坐标系而变化），因此 在坐标变换时也保持不变，故称 I , I , I 分别为应力张 量的第一、第二、第三不变量。 4. 偏应力张量及其不变量 由于任何张量 T 总可以分解为球张量和偏张量两部分， 即 T S Q S P
已知： 则：
n1
s1 s
n2
s2 s
n3
s3 s
P n1 11 n2 21 n3 31 1 P2 n1 12 n2 22 n3 32 P n n n 1 13 2 23 3 33 3
即: 式中 P 为 截面上的应力矢量的各分量，将Pj投影 到任意新的坐标轴上去，可得：
11 21 31 12 22 32 0 3 I1 2 I 2 I 3 0 13 23 33
I1 11 22 33 kk
弹塑性本构关系简介
2） 势能原理的数学表达
应变能
总势能
Ve=Vε+VP =1/2∫VσijεijdV 外力势能
-∫VFbiuidV- ∫SσFsiuidS = min
2 虚力原理
1）虚力原理的表述
给定位移状态协调的充分必要条件为：对 一切自平衡的虚应力，恒有如下虚功方程成 立（矩阵）
∫V[ε]Tδ[σ]dV=∫Su([L]δ[σ])T [u ]0dS
收敛准则
1、位移模式必须包含单元的刚体位移
2、位移模式必须能包含单元的常应变
3、位移模式在单元内要连续、并使相邻单元间的位移必须协调
满足条件1、2的单元为完备单元
满足条件3的单元为协调单元 多项式位移模式阶次的选择——按照帕斯卡三角形选
几何各向同性：位移模式应与局部坐标系的方位无关
多项式应有偏惠的坐标方向，多项式项数等于单元边界结点的自由度总
变间关系为 octσoct
GKtt
oct 3K s oct oct Gs oct
并有
Gs G
1
a
oct
B c
m
KGss
εoct
oct
K G e s
s (c oct ) p
KG
其中G、K分别为初始切线剪切和体积模量，
B c
为混凝土单轴抗压强度，a、m、c和p为由试验
确定的常数。
POCT
弹性张量Dijkl
ij
Dijkl kl
( 2G 1 2
ij kl
2Giklj ) kl
i 1, j 2, k 1,l 2
12
D1212 12
( 2G 1 2
1212
2G1122 )12
11 1 12 0 22 1
迈达斯之——静力弹塑性分析基本原理及方法
m i d a s C i v i lm i d a s C i v i lm i d a s C i v i l图2.8.38 基于位移设计法的结构抗震性能评价m i d a s C i v i l示。
m i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i l1n λ- : 前一步骤(n-1)的荷载因子1λ : 第1荷载步的荷载因子nstep : 总步骤数i : 等差增量步骤号当前步骤的外力向量如下。
0n n λ=⋅P P(10)(3) 第3阶段: 最终步骤的荷载增量(n nstep =) 最终荷载步骤(nstep )的外力向量如下、0nstep nstep λ=⋅P P ; 1.0nstep λ= (11)图2.8.43 自动调整荷载步长的例题(荷载因子结果)m i d a s C i v i l2. 点击步长控制选项 > 增量控制函数定义步长控制函数m i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lm i d a s C i v i lATC-40中对不同结构响应类型规定了谱折减系数的下限值(参见表2.8.7)。
第四章__弹塑性有限元法基本理论与模拟方法讲解
Pi( k ) P(k 1) Pi(k )
q
(k ) 1
k) (k ) qi( k ) qi( q 1 i
q(k )
q( k 1)
第四章 弹塑性有限元法基本理论与模拟方法
(3) 所有载荷段循环，并将结果进行累加
第四章 弹塑性有限元法基本理论与模拟方法
4.2 材料非线性问题及分类
为了与初始屈服应力相区别，我们称之为后继屈服应力。 与初始屈服应力不同，它不是一个材料常数，而是依赖 于塑性变形的大小和历史。 后继屈服应力是在简单拉伸下，材料在经历一定塑性变形 后再次加载时，变形是按弹性还是塑性规律变化的界限。
第四章 弹塑性有限元法基本理论与模拟方法 第四章 弹塑性有限元法基本理论与模拟方法
第四章 弹塑性有限元法基本理论与模拟方法
F
ห้องสมุดไป่ตู้ nom
s0
F nom A0 L nom L0
L
nom
s0
nom (1 nom ) p e ln(1 nom ) E
x1 x x 2 , xn
F(x)=0
f1 (x) f ( x) F ( x) 2 , f n ( x)
0 0 0 0
第四章 弹塑性有限元法基本理论与模拟方法 第四章 弹塑性有限元法基本理论与模拟方法
和简单应力状态相似，材料在复杂应力状态下同样 存在初始屈服和后继屈服的问题。
材料在复杂应力状态下，在经历初始屈服和发生塑性 变形后，此时卸载，将再次进入弹性状态（称为后继弹 性状态）。
第四章 弹塑性有限元法基本理论与模拟方法 第四章 弹塑性有限元法基本理论与模拟方法
弹塑性力学基本知识
非相关联的流动法则。在非相关联流动法则下，塑性应变增量与屈服面不正交，岩土材料的
塑性本构关系一般认为服从非关联的流动法则。 z Mises 屈服条件相关联的流动法则（理想弹塑性）
根据 Mises 屈服准则，结合式（43）
dε
P ij
=
dλ i sij
（46）
式（46）被称为 Prandtl-Reuss 本构关系。结合式（25），总的增量应变与增量应力的关
面在 π 平面上的投影为圆形。根据式（18）可知，Mises 屈服条件的物理意义为：当材料的
八面体剪应力达到一定值时，材料屈服；根据式（26）可知，Mises 屈服条件的物理意义也 为：当材料的剪切应变能达到一定值时，材料屈服。注意，Mises 屈服条件考虑了中间主应 力的影响，但也忽略了静水压力的影响。
（60）
对于 Mises 材料，设材料等向硬化，且内变量为累积塑性应变，结合式（51），有：
2 ∂f ∂f =1
3 ∂σ ij ∂σ ij
结合式（61），（59），（60），可得：
（61）
( ) dλ = dε p ； h = dσ
∫d dε p
（62）
根据式（58），式（32）：
dε ij
=
dε
e ij
可得：
dλ = dW p
∂f ∂σ
sij
ij
（60）
注意：累积塑性应变 ε p 和塑性功W p 都可以作为内变量。
z 基于 Mises 屈服准则的等向硬化模型
使用累积塑性应变作为内变量，即 ε p = dξβ ，结合式（57）和一致性条件式（35），可
得：
(∫ ) h = − ∂f ∂ dε p
2 ∂f ∂f 3 ∂σ ij ∂σ ij
第二章：(2)弹塑性一般知识讲解
T
D
d
A＋ f
T
Dg
＝ D
Dg
f
T
D
d
A＋ f
T
D g
＝Dep d
相适应f=g
2.6 土的剑桥模型（Cam-clay)
2.6 土的剑桥模型
2.6.1 正常固结粘土的物态边界面（state boundary surface)
2.6.2 超固结土及完全的物态边界面 2.6.3 弹性墙与剑桥模型的屈服函数 2.6.4 修正的剑桥模型
2.6.1 正常固结粘土的物态边界面
完全的物态边界面：
CS：v=常数的Roscoe 面 TS：超固结土的强度线-Hvorslev面 0T：零应力线 包括了正常固结土、重超固结土的 可能的（极限）应力状态
包括超固 结土的完 全的物态 边界面
vi-Ti-Si-Ni
HS
超固结
CS
正常 固结
2.6.3 弹性墙与屈服轨迹
1. 弹性墙 正常固结粘土与轻超固结粘土 （wet clay) 各向等压固结: 加载：NCL
NCL
CSL p
NCL
CSL
lnp
正常固结粘土的排水与不排水应力路径
物态边界面与临界状态线
p＝exp((-v)/ ) q=Mp=M exp((-v)/ ) 强度线，物态面与 应力路径的唯一性
v
v=N- lnp：初始加载 v=v- lnp：回弹曲线
lnp
2.6.2 超固结土及完全的物态边界面
土的弹塑性模型25土的弹塑性模型的一般原理254弹塑性本构模型的模量矩阵的一般表达式251塑性理论在土力学中的应用早在1776年库仑公式与土压力理论刚塑性借鉴金属塑性理论弹性理想完全塑性1960s弹塑性理论应用刚塑性perfectlyplastic弹性完全塑性elastoplastic增量弹塑性incrementalelastoplastic不同塑性模型的应用刚塑性理论极限平衡法
弹塑性矩阵推导
对于变形类型，需要判断是纯弹性变形还是弹塑性变形。
若为后者，需要进一步判断塑性应变的正负号及塑性应变的大小。
这是一个比较困难的问题。
MISES 提出的塑性势理论认为，经过应力空间（σ1，σ2，σ3）任何一点，必有一塑性位势等势性存在，它可以表示为：g(σij,H)=0,而塑性应变增量，其变形方向与塑性位势正交，即d εijp=d λ.(g'/σ')。
这个法则与理想流体的流动问题类似，因而称为流动法则（由于是一个梯度的表示，也称为正交法则）（江见鲸，有限元法及其应用P110；陈惠发著，余天庆译.弹性与塑性力学，P182）。
至于准确定出塑性状态下的本构关系矩阵，需要联合用到强化法则、流动法则、应变的弹性分解和塑性分解及弹性区边界上一致性条件（即边界点既在弹性区，又在塑性区）（江见鲸，有限元法及其应用P111-114；陈惠发著，余天庆译.弹性与塑性力学，P182-183）。
总结就是：非线性阶段材料的塑性增量如何确定？（简单说就是其正负号如何定？大小如何定？与应力的关系怎样）流动法则就是来解决这一问题的！考虑材料的塑性，其增量形式的本构关系可表达为p d σ=(D-D )d ε (1)式(1)中，D 为弹性矩阵，p D 为塑性矩阵。
弹性矩阵D 的形式为422000333242000333224000333000000000000K G K GK G K G K GK G K GK GK G G G G ⎡⎤+--⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥=⎢⎥--+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦D (2)体积模量3(12)E K μ=-，剪切模量2(1)EG μ=+。
在应变空间内，塑性矩阵可表达为1()T f fA ∂∂=∂∂p D D D σσ(3) 式中，()()T T p f f f fA B ∂∂∂∂=--∂∂∂∂D D σσσσ(4) f为屈服函数；p σ为塑性应力，p p =σD ε；1/2(()())T p T p T p f f f ff f f ωθε∂∂∂∂∂∂'∂∂∂∂∂∂∂∂σσI σσσ(5) [111000]T '=I ；p ω为塑性功；p θ为塑性体应变；p ε为等效塑性应变；κ为反映加载历史的参数。
T形短肢剪力墙空间弹塑性刚度矩阵
４ 室内不允许设柱， ｍ， 四周为砌体围护， 门窗根据需要开设； 抗震 设计 能迅速 、 有效地对结构 总体方 案进 行构思 、 比较与 选择 ， 为完 设防烈度为 ６度 ； 场地地形基本平坦 ， ． 以下的可塑状红粘 成后面的设计文件提供了正确的概念和思路， 一１５ｍ 也为保证施工的可
土层可作为基础直接持力层 ， 地基承载力特征值 ＝２ ０ｋ ａ ２ Ｐ 。
行性和合理性以确保施工质量、 安全经济地使用创造了条件。而
根据工程实践经验 ， 对于较大 面积 和跨 度的混凝土梁 板屋面 对 于较 复杂的结构设计 ， 往是 由多种 基本 体 系组 合 而成 ， 往 既要
结构 ， 在温度作用下屋 面通常 由于温度变 形会 产生裂 缝 ， 在设 计 考虑充 分发挥组成总结构体系后所具 有的力学特性 ， 利用各 又要 中有针对性地 采用 了以下构造措施防止裂缝 产生 ： ） １ 在屋 面板与 基本体系所 固有 的力学 特性 。因此运 用结 构概念设 计就 愈发显 圈梁 的支 承处 沿墙 体 两端 １４ 围 内设置 了两层 油毡 滑 石粉 。 得重要 。 ／范 这是“ 的措 施。２ 适 当加 大屋面 圈梁 的高度 。３ 在屋 面板 中 抗” ） ） 从上例中还可发现通过概念设计可以采取有效措施解决计 配足温度钢筋 ， 这是“ ” 抗 的措施 。４ 采 用水化热 小和收缩小 的水 算中没有涉及的问题。现今进行结构设计已普遍采用计算机辅 ） 泥（ 如矿渣水 泥 、 粉煤灰水 泥） 选用级配 蘸好 的骨料 ， ， 并严 格控制 助设计 ， 但是建筑结 构计 算机 辅助设计 也 是 以概 念设 计为基 础 ，
非线性有限元-9-弹塑性本构关系
屈服面：
对于单向应力状态，其屈服条件可以写成 s
可以看出，描述一维问题的屈服条件需要应力－应变曲线上的一个临界点
（屈服点），描写多维问题的屈服条件就需要应力或应变空间的一个临界曲面，该
曲面称为屈服面。
考虑到塑性变形与静
水压力无关的特点
f 1,2,3 C
F J2, J3 C
至今已出现许多屈服理论。俞茂宏教授在这方面做出了重要贡献。 屈服函数：
最大剪应力屈服条件。 1870年：圣维南（Saint-Venant）提出在平面情
况下理想刚塑性的应力-应变关系。假设最大剪应 力方向和最大剪应变率方向一致，求解了柱体中发 生部分塑性变形的扭转和弯曲问题、以及厚壁筒受 内压问题。 1871年：莱维（Levy）将塑性应力-应变关系推广 到三维情况。
3) 塑性阶段：继续加载，材料可承受 更大应力，称为材料强化，并伴随 出现塑性应变。至A点以前卸载， 路径接近直线，即坏点：继续加载至可承受的最大 极限应力，试件出现颈缩而破坏，
称为强度极限。
单轴试验下材料的弹塑性性态 （3/3）
强度限 b
A
弹性限 s
其它：1）在强化规律方面，除等向强化模型外， 普拉格（Prager）提出随动强化等模型；2）在实 验分析方面，运用光塑性法、云纹法、散斑干涉法 等能测量大变形的手段。等等。
单轴试验下材料的弹塑性性态 （1/3）
对塑性变形基本规律的认识来自于实验： 1) 从实验中找出在应力超出弹性极限后材料的特性; 2) 将这些特性进行归纳并提出合理的假设和简化模型，
25
二、塑性力学的基本法则
将上述单轴应力状态的基本概念推广到一般的应力 状态，需要利用塑性力学的增量理论。
初始屈服条件
弹塑性力学基本知识
面在 π 平面上的投影为圆形。根据式（18）可知，Mises 屈服条件的物理意义为：当材料的 八面体剪应力达到一定值时，材料屈服；根据式（26）可知，Mises 屈服条件的物理意义也 为：当材料的剪切应变能达到一定值时，材料屈服。注意，Mises 屈服条件考虑了中间主应 力的影响，但也忽略了静水压力的影响。
0
则材料稳定， （2） 加载面 f σ ij , ξ β = 0 外凸。这也可以由式（42）推出。 （3） 正交流动法则（ dλ 的物理意义：反映塑性应变增量的大小，称作比例因子。 ） ：
P dε ij = dλ
(
)
∂f ∂σ ij ∂f s ∂σ ij
（43）
或： dε ij = dλs
P
（44）
p
得：
h=−
∂
( ∫ dε )
p
∂f
2 ∂f
∂f
3 ∂σ ij ∂σ ij
（60）
对于 Mises 材料，设材料等向硬化，且内变量为累积塑性应变，结合式（51） ，有：
2 ∂f
∂f
3 ∂σ ij ∂σ ij
=1
（61）
结合式（61） ， （59） ， （60） ，可得：
dλ = d ε p ； h =
( 2σ 2 − σ 1 − σ 3 )
当采用极坐标表示时，则有：
⎧ rσ = x 2 + y 2 = 2 J 2 ⎪ ⎨ y 1 ⎛ 2σ 2 − σ 1 − σ 3 ⎞ 1 μσ ⎪ tan θσ = = ⎜ ⎟= x 3 ⎝ σ1 − σ 3 3 ⎩ ⎠
z Tresca 屈服条件 当 τ max =
（59）
结合式（43）和式（14） ， （注意：当屈服与静水压力无关，体积应力不产生塑性应变） ， 可得：
5_1弹塑性数值方法
(1.6)
在主应力空间米赛斯屈服面为一无限延伸的圆柱面(图 1)。
图1
主应力空间的 Mises 和 Tresca 屈服面
5
(2) 屈斯加(Tresca)屈服准则
max k
其中最大剪应力 max max 1 2 , 2 3 , 3 1 / 2 。 该准则又称最大剪应力准则， 即当最 材料开始屈 大剪应力达到极限值 k 时， 服。 在主应力空间中， 屈斯加屈服面为 一无限延伸的正六边形柱面(图 1)。
则流动矢量 a d 可写为
a d C1a1 C2a 2 C3a 3
(1.28)
22
对于不同的屈服准则，可将相应的屈服函数 F 代入 (1.26) 式，可以计算得到 C1 , C 2 , C 3 值(表 1)。
表1 计算流动矢量的有关常数
屈服准则 Tresca Mises Mohr-Coulomb Drucker-Prager
1 F m sin J 2 cos sin sin c cos 0 3
0
7
(1.9)
3 1 3J 3 J 2 其中 arcsin 3 2
3 2
。
π 平面上的摩尔－库仑屈服条件如图 3 所示。 试验表明，摩尔－库仑屈服准则较为符合岩土和混凝土材 料的屈服和破坏特征。
11
1.2 硬化法则 硬化法则规定材料进入塑性变形后的后继屈服函数(又称加 载函数或加载曲面)。一般来说加载函数可以采用以下形式
F ij , 0
(1.11)
其中 是硬化参数，它依赖于变形的历史。 对于理想弹塑性材料，因无硬化效应，显然后继屈服函数 和初始屈服函数一致，即
弹塑性有限元分析
角Ｂ 处的偏应力为
２（ 一２ｏ ＋ ，ｒ ４，ｒｃ２＋ｒ ４ ｒ， ｒｃ０． ｒ ， （ 一２ｏＢ ５ 一ｒ ． ， ｓ （ 一） ｒ ， ｓ ， ， ｒ ） ， ２ ４． ） ， ｒ ） （ ＝ Ｂ
的结果与试验资料能很好的吻合。特别是在低静水压力区与试验结果相当吻合
［。 ９ 矩形钢管混凝土析架节点中，核心混凝土处于三轴应力状态， １ ］ 管壁对核心
混凝土产生的约束作用不大，比较适合采用 Ｗｉ ｍ Ｗａ ｋ 的五参数破坏准则。 ｌ － ｍｅ ｌ ａ
该准则在主 应力 （； Ｑ， 空间 破坏面通常 Ｑ， ３ ２ Ｑ） 上的 表示为 静水压力咨 和与 静水 轴相垂直的 压力 偏平面 应力， 相 上偏 与 似角。 数， 的函 如图６ ６ － 所示。 、 咨
ａ ｂ
肪 湘 川
、 ！ ． … １ ． ． 一
应变Ｅ一直延续到ｅ １ １ １ ｚ ０ 的ｂ ＝６
点： ③强化阶段 （‘ ｂ 段） 用一条
直线来简化，该直线一直延续到 抗拉强度ｆ ｕ ，对应的应变为
ｏ ｅ ，
匀
图６ －５钢材的 应力￣应变关系曲 线
３ ０１ ｆ１ ６ １ ６ 取 ｝ ． ④二 塑 ＝ ０ ， ＝ 所。 次
Ｆ ｆａ）ｃ ＝ （； ＝ ， 一 ０ （３ ６０ －）
式中ｃ 为与材料有关的参数。上式退化为单轴应力状态时， 可以 ｃ 根据单轴应力
应变曲 线上定义的屈服强度或破坏强度确定。 此可以引 因 用等效应力和等效应变 的 概念来利用单应力应变曲 线的试验结果。用上式可以 绘制应力空间中的屈服 面， 对于脆性材料即为破坏包络曲 若Ｆ＜ ， 面。 ０ 则处于弹性工作状态； ＿ ， 若Ｆ＞ ０
由 ６ ２） 以 得 元 变 单 结 位 的 系 阵 】单 结 位 式（ ５可 求 单 应 与 元 点 移 关 矩 巨 ， 元 点 移 －
基于M_C准则下非线性材料弹塑性矩阵的推导_方从严
收稿日期： #%%# ’ %$ ’ %Q 作者简介： 方从严 ) $QJJ ’ * ， 男， 安徽桐城人， 讲师， 河海大学在读硕士研究生 -
・!"・ !% " !" " !! " !- .
安
徽
机
电
学
院
学
报
#$$# 年
&（ # ） & !% ! # $ （ !% $ ! # ） ’() " # % & *+’ " , # # &（ # ） & !% ! & $ （ !’ $ ! & ） ’() " # % & *+’ " ， # # &（ # ） & !# ! & $ （ !# $ ! & ） ’() " # % & *+’ " ， # #
$
+ ’ , 屈服准则
参见图 $， 基于 + ’ , 准则的屈服条件为： !$ " !# " $（ # ） $ ! $ !" $ （ !$ $ ! " ） 47F " # % & 584 " I # # $（ # ） $ ! $ !# $ （ !$ $ ! # ） 47F " # % & 584 " I # #
&（ # ） & !# ! % $ （ !# $ ! % ） ’() " # % & *+’ " / # # 式中： 为主应力， 令 !% . $，设 !& 、 !& " !# " !% ， " 为材料的内摩擦角 / 对于平面应力状态， !# 为分别大小主应力， 逆时针） 转到 !& 角度， 得： # 为由 ( 正向（ $ # " !( !) $ !( !) *+’# # $ $()’()# # # # ， ! ( $ ! ) # ! ( # !) " # *+’# # $()’()# # # # # 式中： !(、 !)、 $() 分别为单元体的正应力和剪应力 / 代入得： ’() " 0 & # 1 $ ’() " 0 & # 1 0 ’() " 0 & !& . *+’ # !( ’() # !) ’()# # 1 $() # % 1 *+’ # # # ’() " 0 *+’# # 1 0 ’() " 2 *+’# # 0 !# . !( !) ’()# # 1 $() # % 1 *+’ "， # # ’() " 2 & # 1 $ ’() " 2 & # 1 2 ’() " 2 & !% . ’() # !( *+’ # !) ’()# # 1 $() # % 1 *+’ # # # ’() " 2 & # 1 0 ’() " 2 & # 1 0 ’() " 2 & !" . *+’ # !( ’() # !) ’()# # 1 $() # % 1 *+’ # # # ’() " 2 *+’# # 1 0 ’() " 0 *+’# # 2 !! . !3 !) ’()# # & $() 2 % 1 *+’ "， # # ’() " 0 & # 1 0 ’() " 0 & # 1 # ’() " 0 & !- . ’() # !( *+’ # !) ’()# # & $() # % & *+’ # # #
动力时程分析中弹塑性刚度矩阵的提取方法_姜忻良
第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-07Extraction Method of Elastic -Plastic Stiffness Matrix inDynamic Time History AnalysisJiang 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.1Unit 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.1192 4.593 4.5933 4.707 4.707660 566.051 566.051计算时间/s4.228 0.764图3天津人工地震波矩阵验证Fig.3Verification of matrix of Tianjin artificial wave4 非线性等效刚度矩阵处理方法第3节中的分析验证了ANSYS二次开发在线弹性模型下所提取的矩阵真实有效．若考虑结构在非线性塑性阶段提取其各个时刻的刚度矩阵，可以近似地认为非线性塑性本构关系是分段逐步递减的在小时间段tΔ内的线性本构关系的组合．因此，若要在动力时程分析的塑性阶段生成刚度矩阵，就需要将tΔ时间内视为线性本构关系．对于进入非线性阶段的结构，在小时间段tΔ内，提取t和t+tΔ时刻的各个单元节点的x、y方向的正应力xσ、yσ和剪应力xyτ，以及与之对应的x、y方向的正应变xε、yε和剪应变xyγ．对于平面应变问题，依据如下公式进行推导．各单元的每一个节点的平面应变的物理方程为t时刻2xy21()()()=1()1()()()=1()()()()=x yxyyxy xt ttE tt ttE tttGtμμσσεμμμσστεγμ⎧−⎡⎤−⎪⎢⎥−⎣⎦⎪⎪−⎡⎤⎪−⎨⎢⎥−⎣⎦⎪⎪⎪⎪⎩(1)·358· 天津大学学报(自然科学与工程技术版) 第48卷 第4期t ＋t Δ时刻221(+)(+)(+)=1(+)1(+)(+)(+)=1(+)(+)(+)(+)=x y x y x y xyxy 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 Δ时段的应变差为{222221(+)(+)(+)()=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 xx y y y y y y y y x x x xy xy xyxy 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)可以简化为22()()()=(1)()1()()()()=(1)()1()()()=()y x x x y y x y y x xyxy t t t E t E t t t t E t E t t t G t σσμεμμσσμεμμτγ⎧ΔΔ⎡⎤ΔΔ−ΔΔ−⎪⎢Δ−Δ⎪⎣⎦⎪ΔΔ⎡⎤ΔΔ⎪−ΔΔ−⎨⎢⎥Δ−Δ⎣⎦⎪⎪ΔΔ⎪ΔΔ⎪Δ⎩(5)求解式(5)得到该单元节点的等效弹性模量和等效剪切模量，即22222222()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 μσμεεμμμμμσμεεμμμμτγ⎧⎡⎤ΔΔ−⎪⎢⎥−⎣⎦⎪ΔΔΔ⎪ΔΔ⎪−−−⎪⎪⎡⎤⎪ΔΔ−⎨⎢⎥−⎣⎦⎪Δ⎪ΔΔΔΔ⎪−−−⎪⎪ΔΔΔ⎪ΔΔ⎪⎩(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]．模型的vonMisis等效塑性应变云图如图4所示．将ANSYS和MATLAB两个程序对接，在每一个荷载步结束后进行后处理分析．依据式(6)计算每个单元的等效弹性模量)E tΔ和等效剪切模量)G tΔ，随后赋予该单元材料属性，继而进行模态分析并对整体结构提取其刚度矩阵与质量矩阵，再调用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)·360·天津大学学报(自然科学与工程技术版)第48卷 第4期接计算法和ANSYS-MATLAB联立法得到的结构顶部中点位移、速度、加速度时程曲线见图8．这（a）顶点位移（b）顶点速度（c）顶点加速度图8两种方法的顶部中点动力响应时程曲线对比Fig.8Contrast of dynamic response time history curves of top midpoint between two methods 些动力响应(位移、速度、加速度)时程曲线基本上吻合，说明采用等效弹性模量)E tΔ和等效剪切模量()G tΔ所提取的结构逐渐变化的刚度矩阵较为准确.图9为地震波作用下结构前3阶频率的变化过程，可以看出结构在弹性阶段内频率保持不变，在进入塑性阶段时，结构的频率呈阶梯状降低．图9地震波作用下前3阶频率的曲线Fig.9Curves of the 1st three order frequencies under earthquake action表2和表3汇集了上述Pushover静力非线性分析和天津人工地震波动力非线性分析的前3阶频率在整个过程的某些时间点的变化情况．从Pushover静力非线性分析中可以看出，随着荷载的逐步增加，结构的前3阶频率均逐步降低，而第1阶频率降低最为突出．在天津人工波动力非线性分析中，在地震波达到峰值的时刻前3阶频率的降低速度较大，当地震波趋于平稳的时候其前3阶频率保持恒定，同样也是第1阶频率降低得最多．表2前3阶频率变化情况(Pushover)Tab.2Variation of the 1st three order frequencies(Pushover)时间/s 第1阶频率/Hz 第1阶频率降低率/%第2阶频率/Hz第2阶频率降低率/%第3阶频率/Hz 第3阶频率降低率/%1.119 0 4.593 0 4.706 0 3 1.094 02.23 4.546 01.02 4.684 00.4760.834 25.47 3.955 13.89 4.319 08.228.20.340 69.62 1.965 57.22 2.456 47.81表3前3阶频率变化情况(地震波)Tab.3Variation of the 1st three order frequencies(Earthquake wave)时间/s 第1阶频率/Hz 第1阶频率降低率/%第2阶频率/Hz第2阶频率降低率/%第3阶频率/Hz 第3阶频率降低率/%1.119 0 4.593 0 4.706 0 3 1.114 00.45 4.584 00.20 4.701 00.1160.788 29.58 3.809 17.07 3.957 15.92 90.740 33.87 3.507 23.64 3.745 20.42 120.740 33.87 3.507 23.64 3.745 20.42 150.740 33.87 3.507 23.64 3.745 20.422015年4月姜忻良等：动力时程分析中弹塑性刚度矩阵的提取方法 ·361·6 结 语采用分段等效线性化的处理方法，对非线性弹塑性本构模型在Pushover静力方法与时程分析中进行各种矩阵的提取．结果表明在非线性弹塑性阶段所提取并采用的等效弹性模量()E tΔ和等效剪切模量()G tΔ较为准确，说明在tΔ时间段内将非线性弹塑性本构模型转变为分段等效线性化模型的合理性，同时也验证了本文ANSYS二次开发程序的可行性．本文方法为采用变化的特性矩阵对结构进入塑性阶段后深层次的分析研究提供了有效的方法．参考文献：［1］张令心，石 磊. 土-结构相互作用地震反应分析软件及其二次开发[J]. 地震工程与工程振动，2006，26(3)：225-227.Zhang Lingxin，Shi Lei. Software of seismic soilstructure 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. ANSYSsecondary development for the finite element analysisof prestressed reinforced concrete structures[J]. Jour-nal of Huazhong University of Science and Technol-ogy(Urban Science Edition)，2005，22(Suppl)：87-90(in Chinese).［3］王一功，杨佑发. 针对场地地震反应分析的ANSYS 二次开发[J]. 地震工程与工程振动，2004，24(2)：42-45.Wang Yigong，Yang Youfa. A secondary developmentof ANSYS for site earthquake response[J]. Earth-quake Engineering and Engineering Vibration，2004，24(2)：42-45(in Chinese).［4］刘光栋，王解君，何放龙. 空间梁单元的几何非线性刚度矩阵的分解形式[J]. 湖南大学学报，1992，19(1)：60-71.Liu Guangdong，Wang Jiejun，He Fanglong. Resolved-formulation of geometrical nonlinear stiffness matrixfor three-dimensional beam element[J]. Journal ofHunan University，1992，19(1)：60-71(in Chi-nese).［5］刘齐茂，燕柳斌. 基于 Newmark 法敏度计算的刚架结构动力优化[J]. 工程力学，2010，27(3)：145-154.Liu Qimao，Yan Liubin. Dynamic optimization offrame structures using sensitivity calculation based onnewmark 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 equiva-lent 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 es-timating nonlinear seismic responses of building struc-tures[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 en-ergy spectra of inelastic SDOF systems with equiva-lent linear method[J]. Building Structures，2007，37(8)：74-77(in Chinese).［9］白建方. 复杂场地土层地震反应分析的并行有限元方法[D]. 上海：同济大学土木工程学院，2007.Bai Jianfang. Parallel Finite Element Method forSeismic Response Analysis of Irregular Site[D].Shanghai：School of Civil Engineering，Tongji Uni-versity，2007(in Chinese).［10］张国栋. 土与结构相互作用体系随机地震反应分析[D]. 武汉：武汉大学土木建筑工程学院，2004.Zhang Guodong. Stochastic Seismic Response Analy-sis for Soil-Structure Interaction System[D].Wuhan：School of Civil Engineering，Wuhan Uni-versity，2004(in Chinese).(责任编辑：樊素英)。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
弹塑性矩阵推导考虑材料的塑性，其增量形式的本构关系可表达为p d σ=(D -D )d ε (1)式(1)中，D 为弹性矩阵，p D 为塑性矩阵。
弹性矩阵D 的形式为422000333242000333224000333000000000000000K G K G K G K G K G K G K G K G K G G G G ⎡⎤+--⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥=⎢⎥--+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦D (2) 体积模量3(12)E K μ=-，剪切模量2(1)E G μ=+。
在应变空间内，塑性矩阵可表达为1()T f fA ∂∂=∂∂p D D D σσ(3) 式中，()()T T p f f f f A B ∂∂∂∂=--∂∂∂∂D D σσσσ(4)f为屈服函数；p σ为塑性应力，p p =σD ε；1/2(()())T p T pT p f f f f B f f f ωθε∂∂⎧⎪∂∂⎪∂∂⎪'=⎨∂∂⎪∂∂∂⎪⎪∂∂∂⎩σσI σσσ(5)[111000]T '=I ；p ω为塑性功；p θ为塑性体应变；p ε为等效塑性应变；κ为反映加载历史的参数。
当p κω=时当pκθ=时 当p κε=时对于Drucker-Prager 模型，其屈服条件为120f I J α== (6)1x y z I σσσ=++，22222221()2x y z xy yz zxJ S S S S S S =+++++，α为材料常数。
22f J α∂''=∂I σ (7)222Txyzxyyzzx S S S S S S '⎡⎤=⎣⎦S (8)22()32f K J J αα∂'''==+∂DD I I σ (9)222()()(3)92T T T f f A K K G J J ααα∂∂'''===+∂∂D I I σσ (10)1()T f f A ∂∂=∂∂p D D D σσ22222222222222112123113211212311321121231131121121121(3)(3)999(9)(9)(9)(T TT T T TK K K GJ J K G K G K G J K G J K G J m mn ml S m S m S m mnn nl S n S n S n ml nl l S l S l S l S m S n S l ααααααααβββββββββββββ''=++''''=++++++=p D I I I I S SS 211211222311221321231231231122231231132232113113113113212113223113)()()S S S S S S m S n S lS S S S S S m S n S lS S S S S ββββββββββββββββββββ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⋅⋅⎢⎥⎢⎥⋅⋅⎢⎥⋅⋅⎢⎥⎣⎦令43p K G =+，23q K G =- 弹塑性矩阵可表达为2112123113211212311321121231132112112112112112223112213123123123112223()(p m q mn q ml S m S m S mq mn p n q nl S n S n S n q ml q nl p l S l S l S l S m S n S l G S S S S S S m S n S lS S G βββββββββββββββββββββββ------------------=-=-----⋅-⋅----⋅-ep p D D D 21231132232113113113113212113223113)()S S S S m S n S l S S S S G S ββββββββββ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⋅⎢⎥----⋅-⋅-⎢⎥⎣⎦令122(9)K G J βα=+，229K Gβα=+，1112m S ββ=+，1222n S ββ=+，1332l S ββ=+ 2222222211222222111222299(9)(9)(9)(()(9)9x x x x K G S S S K G K G J K G J K G J S m K G J K Gαααααββαα⎡⎤=++⎣⎦++++==+=++p D22222221222222222221112122299(9)(9)(9)()((9)9(9)9()()y x x yx y K G S S S K G K G J K G J K G J K G J K GK G J K GS S mnαααααααααββββ⎡⎤=+++⎣⎦++++=+++++=++=p D2222222132222222221112133299(9)(9)(9)()((9)9(9)9()()z x x zx z K G S S S K G K G J K G J K G J S S K G J K GK G J K GS S mlαααααααααββββ⎡⎤=+⎣⎦++++=++++++=++=p D2221422111211211222222(9)(9)()()(9)(9)(9)xy x xyx xy G S S K G J K G J S S S S mK G K G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2221522111212312322222(9)(9)()()(9)(9)(9)yz x yzx yz G S S S K G J K G J S S S mK G K G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2221622111211311322222(9)(9)()()(9)(9)(9)zx x zxx zx G S S S K G J K G J S S S S mK G K G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D222221222222222111212229(9)(9)(9)()((9)9(9)9()()x y y xx y S S K G K G J K G J K G J K G J K GK G J K GS S mnααααααααββββ⎡⎤=+⎣⎦++++=+++++=++=p D 2222222222222222122222299(9)(9)(9)(()(9)9y y y y K G S S K G K G J K G J K G J S S n K G J K Gαααααββαα⎡⎤=++⎣⎦++++=+=+=++p D2222222232222222221222133299(9)(9)(9))((9)9(9)9()()z y y Zy z K G S S S K G K G J K G J K G J S K G J K GK G J K GS S nlαααααααααββββ⎡⎤=+++⎣⎦++++=++++++=++=p D2222422122211211222222(9)(9)()()(9)(9)(9)xy y xyy xy G S S K G J K G J S S S nK G K G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2222522122212312322222(9)(9)()()(9)(9)(9)yz y yzy yz G S S K G J K G J S S S S nK G K G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2222622122211311322222(9)(9)()()(9)(9)(9)zx y zxy zx G S S K G J K G J S S S nK G K G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2222222312222222221112133299(9)(9)(9)()((9)9(9)9()()x z z xx z K G S S S K G K G J K G J K G J S S K G J K GK G J K GS S mlαααααααααββββ⎡⎤=+⎣⎦++++=++++++=++=p D222232222222222122213329(9)(9)(9)((9)9(9)9()()y z z yy z S S S K G K G J K G J K G J S K G J K GK G J K GS S nlααααααααββββ⎡⎤=+++⎣⎦++++=++++++=++=p D2222222233222222133222299(9)(9)(9)(()(9)9z z z z K G S S K G K G J K G J K G J S S l K G J K Gαααααββαα⎡⎤=+⎣⎦++++=+=+=++p D2223422133211211222222(9)(9)()()(9)(9)(9)xy z xyz xy G S S K G J K G J S S S S lK G K G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2223522133212312322222(9)(9)()()(9)(9)(9)yz z yzz yz G S S K G J K G J S S S lK G K G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2223622133211311322222(9)(9)()()(9)(9)(9)zx z zxz zx G S S S K G J K G J S S S lK G K G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2224122111211211222222(9)(9)(()9(9)(9)xy xy xx xy G S S K G J K G J S S S S mK GK G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2224222122211211222222(9)(9)(()9(9)(9)xy xy yy xy G S S S K G J K G J S S S nK GK G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2224322133211211222222(9)(9)(()9(9)(9)xy xy zz xy G S S S K G J K G J S S S lK GK G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D222112244222()()(9)(9)xy xy S K G JK G J βαα⎡⎤===⎣⎦++p D1122232452222(9)(9)9xy yz xy yz S S S S S K G J K G J K G ββααα⎡⎤===⋅⎣⎦+++p D 1122132462222(9)(9)9xy zx xy zx S S S S K G J K G J K G ββααα⎡⎤===⋅⎣⎦+++p D 2225122111212312322222(9)(9)()()9(9)(9)yz yz xx yz G S S K G J K G J S S S mK GK G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2225222122212312322222(9)(9)()()9(9)(9)yz yz yy yz G S S S K G J K G J S S S S S nK GK G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2225322133212312322222(9)(9)()()9(9)(9)yz yz zz yz G S S K G J K G J S S S lK GK G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D1122232542222(9)(9)9xy yz xy yz S S S S K G J K G J K G ββααα⎡⎤===⋅⎣⎦+++p D222123255222()()(9)(9)yz yz S K G JK G J βαα⎡⎤===⎣⎦++p D1132232562222(9)(9)9yz zx yz zx S S S S S K G J K G J K G ββααα⎡⎤===⋅⎣⎦+++p D2226122111211311322222(9)(9)()()9(9)(9)zx zx xx zx G S S K G J K G J S S S mK GK G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2226222122211311322222(9)(9)()()9(9)(9)zx zx yy zx G S S K G J K G J S S S nK GK G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D2226322133211311322222(9)(9)()()9(9)(9)zx zx zz zx G S S S K G J K G J S S S S S lK GK G J K G J ααββββααα⎡⎤=+⎣⎦++==+⋅=+++p D1132122642222(9)(9)9zx xy zx xy S S S K G J K G J K G ββααα⎡⎤===⋅⎣⎦+++p D 1132232652222(9)(9)9zx yz zx yz S S S S S K G J K G J K G ββααα⎡⎤===⋅⎣⎦+++p D222113266222()()(9)(9)zxzx S S K G J K G J βαα⎡⎤===⎣⎦++p D。