弹塑性有限元方法

第三章 弹塑性有限元方法的实施§3.1 增量平衡方程和切线刚度矩阵 1、 分段线性化的求解思想塑性变形的特点决定了塑性本构关系的非线性和多值性，上面由塑性增量理论给出了塑性应力—应变关系{}{}ep d D d σε=⎡⎤⎣⎦其中 [][]{}{}[]{}[]{}Tep TFFD DD D FFA Dσσσσ∂∂∂∂=-∂∂+∂∂⎡⎤⎣⎦说明当前应力状态不仅与当前应变有关，而且和达到这一变形状态的路径（加载历史）有关。

这里包含了屈服准则、强化条件和加卸载准则。

由此，对物理非线性问题，通常采用分段线性化的纯增量法和逐次迭代的方法求解。

即将加载过程分成若干个增量步，选择其中任意一个增量步建立它的增量平衡方程并求解，对整个过程的求解有普遍意义。

2、 增量平衡方程和切线刚度矩阵设t 时刻（加载至i -1步终），结构（单元）在当前载荷（广义体力{}v f 和表面力{}s f ）的作用下处于平衡状态，此时物体内一点的应力、应变状态为{}{}σε、。

在此基础上，施加一个载荷增量{}{}v s f f ∆∆和，即从t t t →+∆时刻，则在体内必然引起一个位移增量{}u ∆和相应的{}σ∆、{}ε∆，只要{}{}v s f f ∆∆和足够小，就有{}{}ep D σε∆=∆⎡⎤⎣⎦。

倘若初始状态{}σ已知，加载过程已知，则ep D ⎡⎤⎣⎦可以确定（即pij d ε⎰可以确定，然后可在硬化曲线上得到1p ε所对应的硬化系数）于是上面的方程成为线性的。

在t t t →+∆这一增量过程中，应用于虚功原理可得到如下虚功方程：()()()0eeT T T V V s s V S f f u dV f f u dS σσδεδδ⎡⎤+∆-+∆∆-+∆∆=⎣⎦⎰⎰ （1）根据小变形几何关系u N q B q ε∆=∆∆=∆和，再由虚位移()q δ∆的任意性，并设()()eeT T v v s s VS P P N f f dV N f f dS +∆=+∆++∆⎰⎰，展开后，其中单元在t 时刻载荷等效节点力：eeT T v s VS P N f dV N f dS =+⎰⎰；t ∆内增量载荷的等效力eeT T v s VS P N f dV N f dS ∆=∆+∆⎰⎰。

这样，由方程（1）可得平衡方程： []{}{}eTV B dV P P σσ+∆=+∆⎰ （2）即： ()0eeT T t t V V F B dV B dV P P σσ+∆=+∆-+∆=⎰⎰因为t 时刻（第i 步终）结构处于平衡状态 0eTt V F P B dV σ=-=⎰（3）这样（2）式变为: eT V P B dV σ∆=∆⎰ 即：0eT V F B dV P σ∆=∆-∆=⎰ （4）将{}{}ep d D d σε⎡⎤=⎣⎦和{}B q ε∆=∆代入上式得增量平衡方程：eTep V B D BdV q P F ⋅∆-∆=∆⎰ （5）对增量位移求导：()()eT e ep t V d F B D BdV K d q ∆==∆⎰ （6）于是（5）式成为 et K q P ⋅=∆∆ （7）et K 为单元切向刚度矩阵。

集合所有单刚后得到结构总的增量平衡方程T K q P ⋅=∆∆ （8）方程（8）是线性的，可以直接求解。

3.2 硬化系数H '的数值表示根据单一曲线定理，对于一般稳定性硬化材料，在其简单加载过程中，σ和ε之间存在着一一对应的确定的函数关系()εσΦ=，这一关系可用单向拉伸实验来确定。

例如，对于Mises 各向同性硬化材料p d d H εσ/=' （8）在有限元分析中，作为初始参数应把这一曲线输入（用函数或数字的形式），在加载过程中弹塑性矩阵不断地修改，根据当前的应力或应变来确定。

目前，硬化曲线的输入格式有两种： 1） 解析表达式根据单一曲线定理，由单向拉伸试验曲线直接得出硬化曲线的解析式。

例如： （a ）Mises 各向同性线性硬化材料单向拉伸曲线有： 当 s σσ≤ εσE =当 s σσ≥ ()t s s E εεσσ-+= （9）则有 11111111t p p ett EE d d d H d d d d E E E Eσσσεεεε'=====--- （10） -----------------------------------------------------------------------------------------------------------------附：对于一般材料的硬化曲线的求法（求H '）如单拉曲线 则硬化曲线根据 11~εσ ===》 p εσ~===》⎰⎰==p ij p p d d εεε1 其中单拉时等效应变为 ()113232ευεεε+==ij ij 因为 εε=1，υεεε-==32，平均应变为 ()()132********ευεεεε-=++=所以10εεε∴=+ ，当 21=υ时 εε=1（b ）Mises 各向同性幂硬化材料单向拉伸曲线有： 当 s σσ≤ εσE =当 s σσ≥ m A σε= （11）由屈服点条件：mS SE A εε=得 1mSA E ε-= 据（8）式得 111111111111111p e m e d d EB H d d d E B d d Am E d d σσεεεσσεεε---'=====--⎛⎫⎛⎫-- ⎪ ⎪⎝⎭⎝⎭（12）其中： 1111)/(--==m S m Em Am B εεε2） 根据离散的单拉实验数据，采用样条插值计算H '（参看清华大学孟凡中教材：弹塑性有限变形理论和有限元方法） 3.3 过渡单元弹塑性矩阵的确定1． 三种变形状态弹塑性变形体中，在一个载荷增量步内可能有三种变形状态： 1）弹性区：加载前后均处于弹性状态，故采用弹性阵不变。

2）塑性区：加载前后均处于塑性状态，其弹塑性矩阵ep D ⎡⎤⎣⎦由塑性增量理论确定（与当前应力水平和塑性变形增量的总量有关）3）过渡区：加载前处于弹性状态，加载后进入塑性状态，所以，在这一过程中采用弹性矩阵[]D 或最终的ep D ⎡⎤⎣⎦都不合适，必须寻找一个合适的弹塑性矩阵ep D ⎡⎤⎣⎦。

2． 加权平均的弹塑性矩阵ep D ⎡⎤⎣⎦1）过渡单元在加载后的应力计算（以单拉状态为例）在t ∆时间步内施加一个增量载荷后，讨论某单元的应力应变状态。

设某单元加（卸）载前的应力状态0σ，相应的应变0ε（A 点）处于弹性状态（弹性区间O ’C ）。

加载后，按弹性计算得到应变增量ε∆，到达B 点。

显然B 点不是实际的应力状态，因为已经超过了C 点，进入了塑性变形阶段，假设实际应到达D 点。

该增量步的弹塑性矩阵ep D ⎡⎤⎣⎦是未知的，它的大小应该和该增量步内弹塑性应变所占比例有关，只能经过迭代试算得出。

因为: p e 111εεε∆+∆=∆ （13） 且设 11εε∆⋅=∆m p （14） 则 11)1(εε∆⋅-=∆m e （15） 显然，10≤≤m ，且m =0时是全弹性，m =1时是全塑性。

实际应力增量为111111111(1)(1)e p AC CD ep ep ep D m D m m D mD D σσσσσεεεε∆=∆+∆=∆+∆=-∆+∆⎡⎤=-+∆=∆⎣⎦ （16）推广到一般的应力状态{}σ∆，{}ε∆为()[]1ep ep m D m D D σεε⎡⎤⎡⎤⎡⎤∆=-+∆=∆⎣⎦⎣⎦⎣⎦ （17） ()[]1ep ep D m D m D ⎡⎤⎡⎤=-+⎣⎦⎣⎦（18） ep D ⎡⎤⎣⎦－加权平均弹塑性矩阵；m －比例系数。

2）比例系数m 的迭代公式已知A 点的0σ和0ε，同时由到达A 点的路径确定s ε'和sσ'spεεεεε'-∆+=∆=∆0 由定义：001p s s m εεεεεεεεε''+∆--∆===+∆∆∆ （19）3． 过渡单元m 的确定1）确定是过渡单元。

即在第（i -1）个增量步终（求解结束时）某单元是弹性的应力、应变状态0σ和0ε，且0sσσ'≤（或 0s εε'≤），进入第i 个增量步（t ∆内载荷增量 i P ∆），按弹性计算到达B 点，其应力0s σσ'≥，应变0s εε'≥，可以确定该单元在第i 个增量步内是过渡单元。

2）关于m 的迭代过程：按弹性矩阵[]D 计算该单元的切线刚度矩阵t k ，然后和其他单元集合成总刚T K ，代入结构的增量平衡方程并求解得总位移向量{})1(i q ∆，从{}1i q ∆中提取该单元的{})1(i q ∆，并求出{})1(i ε∆，{})1(i σ∆及{})1(i ε∆。

代入（19）式，计算出)1(m ，再将)1(m 代入（18）得()1ep D ⎡⎤⎣⎦， ep D ⎡⎤⎣⎦中的ep D ⎡⎤⎣⎦与当前应力和应变状态有关。

当前应力为： {}{}{}()11i i σσσ∆+=-，{}{}{}()11i i εεε-=+∆（a ）按第1次迭代的计算值()1ep D ⎡⎤⎣⎦代入该单元计算切线刚度阵，并与其他单元集合组装求解总的增量方程得{}()2iq ∆及相应的{}()2iε∆，{})2(iσ∆及{}{}{}{}()212i i iσσσ∆+=-。

此时0ε和s ε没有改变，再代入（19）式计算)2(m 。

将(2)m 和{}{}2σ（若是硬化材料，还要根据当时塑性应变总量确定H '的值）代入计算()2ep D ⎡⎤⎣⎦。

(b) 依次类推，求出 ()3ep D ⎡⎤⎣⎦；()4ep D ⎡⎤⎣⎦…………直至前后两次的m 值十分接近（到达给定的允许误差范围）停止迭代。

(c ）将迭代终止时的ep D ⎡⎤⎣⎦作为该单元的弹塑性矩阵，求单刚集合，解方程，求出{}n i q ∆及相应的{}n i ε∆，{}n i σ∆将其累加到上一步的终值上作为下一步的初值。

总位移 {}{}{}ni i q q q ∆+=-1{}{}{}ni i εεε∆+=-1{}{}{}ni i σσσ∆+=-1并记下每个单元的{}c σ和s ε' {}{}{}N i i i σσσ∆+=-1，以此作为（i +1）步的初始状态，继续加载。

3) 讨论上面采用的是最简单的纯增量法，并取其中一个增量步（i 步）内对m 值的迭代，最终确定加权平均的弹塑性矩阵ep D ⎡⎤⎣⎦（1） 采用加权平均的弹塑性矩阵ep D ⎡⎤⎣⎦，在同一增量步内，对过渡单元的m 值往往要迭代若干次，每次迭代都要重新计算单元的切线刚度阵，并重新组装总刚和解方程。