第四章-弹塑性有限元法基本理论及模拟方法

课程教学内容：

第一章绪论第二章弹性力学变分原理第三章有限元单元法第四章弹塑性有限元法基本理论与模拟方法第五章刚塑性有限元法基本理论与模拟方法第六章弹塑性有限元变形有限元基本方程第七章塑性加工中的传热问题第八章几种典型材料成形过程计算机模拟分析实例第四章弹塑性有限元法基本理论与模拟方法4.1非线性问题及分类在分析线性弹性问题时，假定：应力应变线性关系结构位移很小（变形远小于物体的几何尺寸）加载时边界条件的性质不变

如果不满足上述条件之一，就称为非线性问题非线性结构的基本特征：变化的结构刚度非线性问题可以分为三类：材料非线性：体系的非线性由材料的应力应变关系的非线性引起。如金属变形弹塑性行为、橡胶的超弹性行为等几何非线性：结构的位移使体系的受力状态发生了显著的变化。如板壳的大挠度问题——平衡方程必须建立于变形后的状态接触非线性：接触状态的变化所引起。如金属成形、跌落试验、多零件装配体等碰到障碍物的悬臂梁（端部碰到障碍物时，梁端部的边界条件发生了突然变化，阻止了进一步的竖向挠度。）板料的冲压成形接触非线性例子随着有限元算法理论、计算机硬件和软件技术的进步及实际工业的需求，CAE技术的应用逐步由线性模拟为主向非线性模拟为主快速发展。1969年，第一个商业非线性有限元程序——Marc诞生。目前几乎所有的商业有限元软件都具备较强的非线性问题的分析求解能力。非线性求解技术的先进性与稳健性已经成为衡量一个结构分析程序优劣的标准。非线性问题的有限元求解方法非线性方程(组)的求解方法直接迭代法Newton-Raphson迭代法修正的Newton-Raphson迭代法非线性问题通常采用增量法求解（追踪加载过程中应力和变形的演变历史。）每个增量步采用Newton-Raphson迭代法非线性问题有限元控制方程：非线性方程的迭代求解方法直接迭代法Newton-Raphson迭代修正的N-R迭代非线性方程组的迭代求解方法直接迭代法N-R迭代修正的N-R迭代第四章弹塑性有限元法基本理论与模拟方法非线性问题的增量法求解过程(1)将总的外力载荷分为一系列载荷段(2)在每一载荷段中进行迭代，直至收敛(3)所有载荷段循环，并将结果进行累加(1)将总的外力载荷分为一系列载荷段(2)在每一载荷段中进行迭代，直至收敛N-R迭代:(3)所有载荷段循环，并将结果进行累加4.2材料非线性问题及分类概念：由于材料的应力应变非线性关系引起的非线性。分类：不依赖时间的弹、塑性问题非线性弹性——橡胶弹塑性——冲压成形依赖于时间的粘（弹、塑）性问题蠕变——载荷不变，变形随时间继续变化松弛——变形不变，应力随时间衰减非线性弹性材料行为橡胶应力应变关系曲线弹塑性材料进入塑性的特征：载荷卸去后存在不可恢复的永久变形。应力应变之间不是单值对应关系，与加载历史有关。单轴应力状态下弹塑性材料行为单轴（一维）应力状态下材料的应力应变行为可以从拉伸试验中获得。单调加载硬化塑性理想弹塑性各向同性硬化:运动硬化:混合硬化:反向加载运动硬化各向同性硬化混合硬化

在简单拉伸的情况下，当材料发生塑性变形后卸载，此后再重新加载，则应力和应变的变化仍服从弹性关系，直至应力到达卸载前曾经达到的最高应力点时，材料才再次屈服（后继屈服）。这个最高应力点的应力就是材料在经历了塑性变形后的新的屈服应力。由于材料的强化特性，它比初始屈服应力大。第四章弹塑性有限元法基本理论与模拟方法为了与初始屈服应力相区别，我们称之为后继屈服应力。与初始屈服应力不同，它不是一个材料常数，而是依赖于塑性变形的大小和历史。后继屈服应力是在简单拉伸下，材料在经历一定塑性变形后再次加载时，变形是按弹性还是塑性规律变化的界限。第四章弹塑性有限元法基本理论与模拟方法

和简单应力状态相似，材料在复杂应力状态下同样存在初始屈服和后继屈服的问题。

材料在复杂应力状态下，在经历初始屈服和发生塑性变形后，此时卸载，将再次进入弹性状态（称为后继弹性状态）。第四章弹塑性有限元法基本理论与模拟方法

把复杂应力状态下，确定材料后继弹性状态的界限的准则就称为后继屈服条件，又称为加载条件。问题:

当材料处于后继弹性状态而继续加载时，应力（或变形）发展到什么程度材料再一次开始屈服呢？第四章弹塑性有限元法基本理论与模拟方法一般应力状态下弹塑性材料行为屈服准则（初始屈服条件）硬化法则（后继屈服函数、加载函数、加载曲面）流动法则加载、卸载准则屈服准则（初始屈服条件）在单向受力情况下，当应力达到材料的屈服强度时材料开始产生塑性变形。对于一般复杂的应力状态，应力状态由六个应力分量决定时，显然不能根据某个单独应力分量的数值作为判断材料是否进入塑性变形的标准。为此，引入以应力分量为坐标的应力空间，根据代表不同应力路径的实验结果，可以定出从弹性阶段进入塑性阶段的各个界限，即屈服应力点。在应力空间中，这些屈服应力点形成一个区分弹性和塑性的分界面——屈服面。描述这个屈服面的数学表达式就是我们所要寻求的一般应力状态下的屈服准则。常用的各向同性Von-Mises屈服准则：各向同性屈服准则：各个方向屈服应力相同各向异性屈服准则：不同方向屈服应力有差异三维主应力空间π平面上的屈服轨迹σ3=0平面上的屈服轨迹第四章弹塑性有限元法基本理论与模拟方法硬化法则塑性硬化法则规定了材料进入塑性变形后的后继屈服函数（又称加载函数或加载曲面）各向同性硬化运动硬化混合硬化运动硬化：该模型假设材料随塑性变形发展时，屈服面的大小和形状不变，仅是整体在应力空间作平动。

各向同性硬化：材料进入塑性变形以后，屈服面在各方向均匀地向外扩张，其形状、中心及其在应力空间的方位均保持不变。

材料的强化只与总的塑性变形功有关而与加载路径无关。

应力有反复变化时，等向强化模型与实验结果不相符合。

混合硬化：其实质就是将随动强化模型和等向强化模型结合起来，即认为后继屈服面的形状、大小和位置一起随塑性变形的发展而变化。该模型能够更好的反映材料的Bauschinger效应。各向同性硬化运动硬化第四章弹塑性有限元法基本理论与模拟方法流动法则塑性应变增量和应力分量的关系：塑性应变沿后继屈服面F=0的法线方向——是一正的待定系数，其具体数值和材料硬化准则有关加载、卸载准则对于硬化材料（当材料处于某一塑性状态）：4.3几何非线性问题及分类概念：由于大位移、大转动而引起的非线性。分类：大位移、大转动、小应变问题

——板壳的大挠度和后屈曲大位移、大转动、大应变问题

——薄板成形、弹性材料的受力比较：线弹性—几何非线性线弹性：小变形假设——假定物体发生的位移远小于物体本身的几何尺寸，应变远小于1。建立平衡方程时不考虑物体位置和形状的变化。几何非线性：物体发生有限变形——大位移、大转动的情况。建立平衡方程时必须考虑物体位置和形状的变化。4.4弹塑性矩阵

应力与应变的关系有各种不同的近似表达式和简化式。根据普兰特尔—罗伊斯（Prandtl-Reuss）假设和密赛斯屈服准则，当外作用力较小时，变形体内的等效应力小于屈服极限时为弹性状态。当外力增大到某一值，等效应力达到屈服应力，材料进入塑性状态，这时变形包括弹性变形和塑性变形两部分，即：

式中下脚e、p分别表示弹、塑性状态。

(4-10)在弹性阶段，应力与应变关系符合虎克定律。进入塑性状态后，符合Prandtl-Reuss假设。4.4.1弹性阶段

在弹性阶段，应力和应变的关系是线性的，应变仅取决于最后的应力状态，并且一一对应，而与变形过程无关，有下列全量形式：

式中为弹性矩阵。

(4-11)对于各向同性材料，由广义虎克定律可得：

或：(4-12)式中：是材料的弹性模量，是泊松比。对于各向同性材料，广义虎克定律：

是材料的剪切弹性模量式中：是材料的弹性模量，是泊松比。公式（4-12）的具体推导：

（2）+（3）有：将其带入（1）得：将其带入（1）得：同理可推得得表达式，写成矩阵的形式，就是：4.4.2弹塑性阶段

当材料所受外力达到一定值时，等效应力达到屈服极限，应力应变关系曲线由弹塑性矩阵决定，现推导弹塑性矩阵。

等效应力为：

对应力求导得：

式中为应力偏量，

(4-13)(4-14)公式（4-14）对应力求导的具体推导：

由普兰特尔—罗伊斯关系有：将式4-14代入式4-15得：写成矩阵形式为：(4-15)(4-16)(4-17)公式（4-16）的具体推导：

式中:

又因有：

写成矩阵乘积的形式为：(4-18)(4-19)设为硬化曲线上任一点的斜率，即将式4-20代入式4-19中得：

将式4-11写成增量形式为：

(4-20)(4-21)

(4-22)

再利用式4-10就可得到：

两边同乘以后可得:利用式4-17和式4-21，可将上式写成：(4-23)(4-24)(4-25)由此得：将式4-26代入式4-17得：(4-27)

(4-26)将式4-27代入式4-23得：(4-28)

由式4-14有：则:(4-29)(4-30)

将式4-12代入式4-30，并注意到：得：因有：

(4-31)

(4-32)

(4-33)

因为：所以：公式（4-31）的具体推导：

故：注意：因为：所以：公式（4-33）的具体推导：

即：那么：或者：故令：式4-28可写成：(4-34)

(4-35)

利用上述关系式可将式4-34表示成显式，即：(4-36)(4-31)公式（4-36）的具体推导：

对于平面应力状态，，则有：(4-38)(4-37)

公式（4-37）的具体推导：

对于平面应力状态，，则有：(1)(2)得：同理可推得的表达式:写成矩阵的形式，就是：按照上述同样方法可得：式中：在塑性区：(4-39)公式（4-39）的具体推导：

所以：

由普兰特尔—罗伊斯关系有：写成矩阵形式为：仿照前面，不难推得：而：其中：对于平面应变问题，有

只需从式4-12和式4-36中消去上述为零的分量，就可得到下列各式。(4-40)(4-42)(4-41)4.5变刚度法

变刚度法又称切线刚度法，它所采用的应力与应变关系见图3-1。在等效应力达到屈服极限后，应力与应变不再是线性关系，而是由下列关系式所确定。(4-43)

弹塑性矩阵[D]ep中含有应力，它是加载过程的函数。直接求解是困难的，通常采用增量形式来近似代替微分形式，这样使求解成为可能。计算中由于[D]ep在∆{σ}范围内变化不大，因此可假设在每一加载步中是一个常数，并以该加载步前的应力状态近似计算出[D]ep，即：

单元刚度矩阵[k]在一个加载步中也同样取作常数，即：

(4-44)(4-45)

在一个变形体中，不仅各点的应力状态是不相同的，而且随着加载而变化着，通常变形体受外力作用时，从一个区域到另一个区域，等效应力是逐渐地达到屈服极限，即进入塑性（弹塑性）状态。为了简化，本章所指进入塑性即为弹塑性状态，这就是说在变形体中，各单元的应力和应变状态是不一样的，随着加载又是变化的，且各有各的变化规律。变形体内的单元按状态可分为三类:

弹性单元塑性单元过渡单元各类单元有不同的本构关系和单元刚度矩阵。在加载过程中，各单元的状态是变化的，为此[K]也是变化的。在计算中，每增加一个载荷增量，就得重新计算一次整体刚度矩阵[K]，这也就是变刚度法名称的由来。式中[K]整体刚度矩阵；

n1、n2、n3分别为弹性单元、塑性单元和过渡单元的数量；

[k]e、[k]ep、[k]g

分别为弹性单元、塑性单元和过渡单元的单元刚度矩阵。对于整体来说，可用下列关系式表示：(4-46)整体刚度矩阵求得之后，就可根据下列载荷和位移的线性方程组求解出未知的节点位移增量。有了节点位移增量就能求得各单元的应变及应力增量。(4-47)(4-49)(4-48)4.5.1定加载法

定加载法又称等量加载法。它每次的加载量是预先给定的。这种加载法的加载量一般较大。由于每次加载量较大，每次加载中由弹性单元转变为弹塑性单元的过渡单元较多。过渡单元在加载步中达到屈服，式中m为加权系数，0≤m≤1，采用不同的加载方法，过渡单元的处理也有所不同。下面介绍几种加载方法。(4-50)

m的取值需要进行迭代来逼近，收敛性一般都很好，只需进行2～3次迭代就能达到满意的精确度。定加载法计算程序框图4.5.2变加载法这种方法又被称作r因子法。用这种方法计算，每次加载量是变化的，其大小是由计算结果来确定。计算开始时预先施加一个单位载荷增量，然后求出各单元在施加单位载荷增量后的等效应力增加量。根据这个增加量求出各弹性单元当达到屈服时所需要施加的增量值，最后取这些增量值中最小的一个增量值作为本次加载的加载量。各弹性单元的加载因子按下式进行计算：式中,i单元前次加载后的等效应力；

,i单元本次施加单位载荷增量后的等效应力；

,i单元达到屈服所需施加单位加载量的倍数。(4-51)为了加快计算步伐，常假设单元的等效应力接近屈服极限时就由弹性单元转变为弹塑性单元。一般可取单元等效应力，在下一次施加载荷增量的计算中，这单元就按弹塑性单元处理。采用这种处理方法，能保证每次加载后弹性单元中等效应力的最大者正好达到屈服。在下一次加载中该单元按弹塑性单元处理。这种方法能避免在每个加载步中单元由弹性转变为弹塑性所需要迭代计算m因子的过程，还能保证足够好的计算精度。变加载法的计算程序框图4.5.3位移法

对于有些问题如圆柱体镦粗，每次施加的增量不是控制加载力，而是控制压下量。假设工具为刚性体，在工具与坯料的接触面上，各节点的位移相同。而接触面上的压力分布是未知的。这样在计算中需对与接触表面上节点有关的方程进行处理。计算这类问题时，一般先假设接触表面上有一个已知的轴向位移增量，根据这个已知位移增量求出各单元的应力和应变增量。这种施加位移增量的方法又分为定位移增量法和变位移增量法。4.5.3.1定位移增量法定位移增量法每次施加的位移增量相同。计算过程先算出各单元的应力和应变增量，有了应力增量，累加后可得全应力和计算得到等效应力。根据等效应力的大小将单元分成弹性单元、塑性单元和过渡单元。对于过渡单元的处理与定加载法相同，为此这种方法在求取m因子时也需要进行迭代。

4.5.3.2变位移增量法

这种方法与变加载法相似，每次施加的位移增量由计算结果来决定。每步计算也是先施加一个单位位移增量，然后根据这位移增量计算应力增量，从而找出弹性单元中等效应力增长最快而又最先达到屈服极限的单元，以这个单元达到屈服所需要的位移增量为本次施加的位移增量。达到屈服极限的弹性单元在下次计算中按弹塑性单元处理。

4.6初载荷法

初载荷法是将塑性变形问题试图转化为弹性问题来求解，它把塑性变形部分视作初应力或初应变来处理。

在弹性有限元中，当弹性体的单元中存在初应变{ε}，如因温度而引起的应变，或有初应力，如残余应力，则应力和应变关系分别为：(4-52)(4-53)这样可得有初应变时的变形能为：即：

(4-54)(4-55)因：得：上式中最后一项为与节点位移无关的变形能，由此可得到与节点位移有关的变形能为：式中∑是对所有单元求和。

(4-56)(4-57)通过变分可得到初载荷时的有限元公式为：

同样可导出有初应力时的变形能为：这样，基本方程就比无初应变或无初应力时多一项{R}，{R}是作为载荷存在于基本方程中，称为载荷向量。它是由于有初应变或初应力而引起的，下面分别导出{R}的表达式。(4-58)(4-59)4.6.1初应力法

在小位移的弹塑性问题中，应力应变关系为：

为使问题线性化，对于区域中已达到屈服的单元，采用逐次加载法。如每次加载的载荷较小，可将上式的微分形式近似地写作增量形式，即：

(4-60)(4-61)得：

因有：将式:

代入上式，并将右边表示成矩阵乘积形式，变分后得到：

令：由前面可知：(4-64)(4-63)(4-65)(4-62)上式要表示成（*）的形式，必有：

写成增量形式为：式中：(4-66)(4-68)(4-67)由上式可以看出，初载荷{R}不仅与加载前的应力有关，而且与该次加载引起的应变增量有关，即式：的两边都含未知数。因此每此加载时，必须利用迭代法来求解。具体迭代过程：可先取，则得，求得初载矢量，在载荷下，解方程式：此时方程可写成：(4-69)解方程：得到后，再照上述方法依次迭代计算下去。写作一般迭代式为：

当相邻两次迭代所得初应力相差甚小时，可认为迭代结束。

求得后，再由求和初载荷矢量，然后进行第二次迭代计算。(4-70)(4-71)为了考虑过渡单元，运用前面引入的加权系数m，这样初载荷矢量为：当单元为弹性状态时，取m=1单元为塑性状态时，取m=0若为过渡单元，则取计算出的m值。(4-72)

计算步骤如下：

⑴施加全部载荷{P}

于结构，按线弹性计算；⑵算出各单元的等效应力，并取出最大值σmax

。若σmax≤σs

，则材料尚未达到或刚巧达到屈服，所计算结果就是最终结果。否则，令β=σs/σmax

，则β{P}是恰好使等效应力为σmax的单元达到屈服的载荷，同时存贮由β{P}

产生的应力、应变和节点位移。余下的载荷{P}−β{P}

如分n

次加载完，每次施加的载荷增量为：⑶再施加载荷增量∆{P}

于结构；⑷计算屈服单元和过渡单元的[D]p，应力取加载前的值，过渡单元取达到屈服时的应力值；⑸对载荷增量∆{P}

进行纯弹性计算，求得各单元的全应变增量∆{ε}；⑹用迭代所得的全应变增量，重新计算初应力；⑺求出相应的初载荷，与∆{P}一起作用，按弹性问题求解得节点位移增量和应变增量；⑻重复步骤⑹、⑺直到相邻两次计算所得的初应力非常接近时为止；

⑼

求出应力增量，把位移增量、应变增量、应力增量迭加到加载前的数值上；⑽输出本次加载后的计算结果；

⑾载荷若全部加完，则停机。否则，回到步骤⑶继续计算。

4.6.2初应变法

有初应变时，与节点位移有关的变形能为：将上式代入式4-64后，进行变分运算可得：(4-73)(4-64)

由此可得：

在小位移弹塑性问题中，应力应变关系可写成：把视为初应变，则与弹性问题中存在初应变的情况相似，只要在基本方程的右边加上由引起的初载荷矢量。由于使用的只是弹性矩阵，所以其相应的刚度矩阵计算与弹性时相同。

(4-74)(4-75)把塑性应变增量视为初应变，则:

由式4-17与式4-21可得：

写作增量形式为：(4-76)(4-78)(4-77)应用式4-76和式4-78的关系，式4-72可写成：(4-79)

这样，初载荷矢量不仅与加载前的应力有关，而且与这次加载的应力增量有关。因此在基本方程式4-67的等式两边都有未知数，必须用迭代的方法求解。迭代公式与初应力相似，即：(4-80)

当i=0时，∆{σ0}=0，{R0}=0，按纯弹性计算，当迭代中相邻二次计算所得的∆{ε}p相差很小时，则迭代完成。每次加载时，的计算可利用本次加载前的应力进行计算。以上表明:对于变刚度法，每次加载都必须重新计算刚度矩阵。对于初应力或初应变法，每一步加载只需求解一个具有相同刚度矩阵的问题。所以，变刚度法的计算量一般比初应力法或初应变法大一些。但是，后两种方法每加载一步都必须对初应力或初应变进行迭代，于是就存在迭代是否收敛的问题。对于应变强化的材料，初应力法一定收敛；一般初应变法也是收敛的。但对理想塑性