1、课程教学内容： 第一章 绪 论第二章 塑性成形分析的理论基础第三章 有限元法基本概念第四章 弹塑性有限元法基本理论与模拟方法第五章 刚塑性有限元法基本理论与模拟方法第六章 几种通用有限元分析软件介绍（ANSYS、MARC、ABAQUS）第七章 几种典型材料成形过程计算机模拟分析实例,刚塑性有限元法是在1973 年提出来的，这种方法虽然也基于小应变的位移关系，但忽略了材料塑性变形时的弹性变形部分，而考虑了材料在塑性变形时的体积不变条件。 它可用来计算较大变形的问题，所以近年来发展迅速，现已广泛应用于分析各种金属塑性成形过程。 刚塑性有限元法的理论基础是变分原理，它认为在所有动可容的速度场中，使泛
2、函取得驻值的速度场是真实的速度场。根据这个速度场可以计算出各点的应变和应力。,5.1 刚塑性有限元法及其变分原理介绍,泛函是函数的函数;在泛函进行变分时根据其有无附加条件而分为一般变分和 广义变分或条件变分；广义变分又分为不完全广义变分和完全广义变分。,对于实际的金属成形加工过程，弹性变形部分远小于塑性变形部分 ( 弹性应变与塑性应变之比通常在1/1001/1000 )，因而可忽略弹 性变形，将材料模型简化为刚塑性模型。 采用刚塑性模型可大大简化有限元列式和求解过程。 与弹塑性有限元法相比较，可采用较大的时间增量步长。在保证足 够的工程精度的前提下，可提高计算效率。,由于刚塑性有限元法采用率方
3、程表示，材料变形后的构形可通 过在离散空间对速度的积分而获得，从而避开了应变与位移之 间的几何非线性问题。 由于忽略了弹性变形，刚塑性有限元法仅适合于塑性变形区的 分析，不能直接分析弹性区的变形和应力状态，也无法处理卸 载和计算残余应力与变形。 由于刚塑性模型假设，对一般的体积不可压缩材料，因为其静 水压力与体积应变率无关，如要计算应力张量，还必须进行应 力计算的处理。,从数学的角度来讲，有限元法是解微分方程的一种数值方法。它的基本思想是：在整个求解区域内要解某一微分方程很困难(即求出原函数)时，先用适当的单元将求解区域进行离散化，在单元内假定一个满足微分方程的简单函数作为解，求出单元内各点的
4、解；然后，再考虑各单元间的相互影响，最后求出整个区域的场量。刚塑性有限元法的求解过程 (1) 离散化处理； (2) 单元分析的基础上集合成总体方程组； (3) 刚塑性有限元法集合成的总体方程组为一非线性方程组，还须 线性化处理并采用迭代方法求解。,刚塑性有限元法按照处理方法的不同分成如下5种： (1)流函数法；(2)拉格朗日乘子法；(3)罚函数法；(4)泊松系数v 接近0.5法；(5)材料可压缩性法。,5.1.1 刚塑性材料基本假设 对于大变形金属塑性成形问题，将变形体视为刚塑性体，即把变形中的某些过程理想化，便于数学上处理。此时，材料应满足下列假设： (1) 不计材料的弹性变形； (2) 材
5、料的变形流动服从 Levy- Mises 流动法则； (3) 材料是均质各向同性体； (4) 材料满足体积不可压缩性； (5) 不计体积力与惯性力； (6) 加载条件(加载面)给出刚性区与塑性区的界限。,5.1.2 第一变分原理 刚塑性材料的第一变分原理又称为马尔柯夫(Markov)变分原理, 其为： 在满足：(1) 速度-应变速率关系 (在 上)的一切动可容场 ， 中使泛函： 的变分为零，即： ，且取极小值的 ，必为本问题的真实解。,(2) 体积不可压缩条件,(3) 速度边界条件,证明：设真实解为 和 ，而许可解 由屈服条件和本构方程有： (a) 则有： 由最大塑性功原理，有： (b),由虚
6、功率原理得： (c) 将(c)式代人(b)式得：(d) 注意，在Su上。将(d)式代人(e)式有： (e) 将(a)式代人(e)式有： (f),即： 因此，泛函取最小值， 于是第一变分原理得证,5.1.3 完全广义变分原理在第一变分原理中，所选择的速度场必须满足(1)，(2)和(3)式，实际问题中，有些条件比较容易满足，而有些条件则不易满足。 为了容易选择速度场，应用条件变分的概念，引用拉格朗日(Lagrangian) 乘子 ， 和 ，将运动许可解所必须满足的条件引入泛函中， 则得到新的泛函： (*) 在任意选取的 、 中，真实解使(*)式的泛函取驻值，这就是刚塑性 完全广义变分原理。,第一变
7、分原理和完全广义变分原理对比第一变分原理所选择的 和 只要求满足运动许可条件，而静力许可条件是通过变分近似满足的。据广义变分原理预选的 和 不受任何约束，所有的方程均由变分近似满足。所以，由第一变分原理计算的近似解较广义变分原理得到的解更精确。但前者在预选满足运动许可条件的速度场时比后者困难。,5.2.3 不完全广义变分原理 在选取运动许可解 和 时，可将其应满足的三个条件中的任意两个或一个事先得到满足，而将其余的一个或两个，通过拉格朗日乘子引入泛函中，组成新的泛函，真实解使泛函取驻值，这就是不完全广义变分原理。在选择速度场时应变速率与速度的关系(1)式和速度边界条(3)式容易满足，而体积不可
8、压缩条件(2)式难于满足。因此，可以把体积不可压缩条件用拉格朗日乘子入引入到泛函中，得到新泛函：(*) 可以证明，在一切满足应变速率与速度关系和速度边界条件的 中，使泛函(*)式取驻值的 为真实解。,按照Markov变分原理求解时，面临速度场选取的困难。因而在实际求解时常采用不完全广义变分原理求解塑性变形过程。对刚塑性体和刚粘塑性体，按Markov变分原理确定的泛函为： (*),解决的问题是寻找某种方式将体积不可压缩条件(2)式引入泛函(*)中，构成新的泛函，使问题转变成对新泛函的无条件的驻值问题。通常采用拉格朗日乘子法、罚函数法及修正罚函数法来构造新的泛函。通过这样的方法将体积不可压缩条件引
9、入后，便能求静水压力 ， 从而解决了因忽略材料的弹性变形而带来的应力计算的困难。,5.2 刚塑性增量理论的广义变分原理 欲求解变形体在塑性变形时的场变量，首先要建立基本方程组。,5.2.1 基本方程 基本方程如下： 微分平衡方程或运动方程: (5-1) 速度与应变速率的关系： (5-2) 式中： 速度； 应变速率 列维密赛斯应力应变速率关系： (5-3),假设材料符合密赛斯屈服准则，即: 式中 k 是变形过程的函数，如材料是理想刚塑性体时， k=const式 5-3 两边平方后得： 将式5-4 代入式5-5 整理后得： (5-6) 将式5-6 代入式5-3 可得： (5-7) 这就是符合密赛斯
10、屈服准则的应力应变关系式。,(5-4),(5-5),体积不可压缩条件: (5-8) 边界条件: 边界条件分为力学边界条件和位移边界条件，分别为： (5-9) (5-10) 利用上述方程和边界条件，虽然在理论上是可以求解的，但实际上很 困难，只有在几种简单情况下才能求出解析解。,刚塑性有限元法借助于变分原理可求出近似解，对变形场的位能泛函进行变分，当变分取得驻值时，变形场满足平衡微分方程和力学边界条件。处理体积不变条件的方法有两种： 一是在假设初始速度场时，除了满足速度边界条件以外，还应严格满足体积不变条件，这种方法给假设初始速度场带来困难。另一种方法是假设初始速度场只满足速度边界条件，而对体积
11、不变引入约束条件，即拉格朗日乘子进行有条件变分。这种方法在运算中较易实现，目前已得到广泛应用，下面对这种方法进行详细论述。,5.2.2 不完全的广义变分原理 刚塑性有限元计算需要先选择初始速度场。在选择初始速度场时，速度边界条件容易满足，而体积不可压缩条件较难满足。因此，把体积不可压缩条件用拉格朗日乘子引入泛函中去。这种有条件的但并非将所有条件引入泛函的变分称之为不完全的广义变分，所建立的泛函为：,式中：Sp 变形体边界中应力边界部分； 克罗内克尔（Kronecker）符号。,(5-11),刚塑性不完全的广义变分原理认为：在所有满足速度-应变速率关系和速度边界条件的 vi 中，使泛函式 5-1
12、1 取得驻值的 vi 是真实解。 在忽略体力的情况下，式5-11 还可写成另一种形式，即,式中：,(5-12),注意在刚塑性有限元中，利用的屈服准则是密赛斯屈服条件，它的一阶导数是连续的，在计算中一般略去体力，并设外力在变形过程中不变。对于有硬化的材料，假设剪切屈服极限在一小段范围内是常数，采取台阶形硬化曲线来代替真实硬化曲线，这样处理可大大简化变分的运算。下面证明这一原理。,由式 5-12 变分得： 采用密赛斯屈服准则，有：,(5-13),(5-14),将式5-14 代入式5-13 得：,(5-15),且 ，在 上 ，由此得：,因为：,(5-16),将式 5-16 代入式 5-15 得：,因
13、 和 都是任意的，要使泛函的变分为零，即 取得驻值，必须满足下列等式： （在 表面上） (5-18) （在 V 体积内 ） (5-19) （在 V 体积内 ） (5-20) （在 V 体积内 ） (5-21),由式5-21 可知，泛函变分为零时，满足体积不变条件。由式5-20 可得： 拉格朗日乘子等于平均应力，即静水压力。这就是拉格朗日乘子的物理意义。 由式5-19 可知，泛函变分为零时，在整个体积内都满足运动方程，即平衡微分方程。 由式5-18 可以看到，泛函变分为零时，满足应力边界条件。 这就证明了在满足速度边界和应变速率与速度关系的速度场vi 中，当泛函变分为零时，满足所有的基本方程，所
14、以这个速度场就是真实解。,(5-22),5.3 拉格朗日乘子法为使有限元计算方便，将式 5-11 改写成矩阵形式如下： (5-23) 式中： 应变速率列阵； 速度列阵； 应力边界 Sp 上给定的表面力列阵； 矩阵记号， 体力列阵。,在计算中如材料有硬化作用时，采用阶段硬化曲线来代替真实硬化曲线，如下图所示。经这样的处理，变分时可将 k 视作常数，可从积分号中提出。计算中若忽略体积力，则泛函又可写成如下形式： (5-24),注意上式中：,阶段硬化曲线来代替真实硬化曲线,5.3.1 离散化 假设变形体被划分为 M 个单元，N 个节点，由此可知： (5-25) 对于一个单元而言，可建立下列泛函： (
15、5-26) 式中： V单元的体积； S单元的边界。,在单元内有： (5-27) (5-28) 其中 为单元节点的速度列阵。,将式5-27 和式5-28 代入式5-26 得到： (5-29),令： （这里 K不是刚度矩阵）,对于一个单元来说，节点的速度 和 都是定值。所以式5-29 可写成： (4-30),泛函 中只含单元的节点速度 和 ，未知数为 、 即：,集合成整体，得： 泛函变分为零（ ），即得： (5-32),由于变分， 和 是任意的独立变量，所以有： (5-33),式中： k 求解问题的维数。,(5-31),对于每个单元有：,由式 5-34 集合成的方程组是非线性的，求解时需先进行线性
16、化， 下面就此问题展开讨论。,(5-34),5.3.2 方程的线性化 求解非线性问题的一种常用方法是摄动法，这种方法是先假设一个初始解，根据这个解求出修正量，利用修正量修改原初始解，再由修正后的解求出新的修正解，如此这样通过反复迭代来逼近真解。采用这种求解方法就能把非线性方程组化为线性方程组来求解。,设有一个初始速度场u和相对应的速度增量u，则每次迭代之间的速度有如下关系：,将式5-35 代入式5-34 得： (5-36),现令：,(5-35),略去高次小量，并注意 是一个数，即：,把形如 的因子展成幂级数并取线性项得：,(5-37),因： 令： 则： 将上式代入式5-37 后得,(5-38)
17、,略去高次微量得： (5-39),再令： 得： (5-40),将式5-40 写成矩阵形式为： (5-41),其中：,(5-42),将单元分析得到的式5-41 代入式5-33，就得到刚塑性有限元法求解的矩阵方程组，即：,(5-43),求解真实速度场时采用迭代法，其迭代的收敛判据取范数比 ，当范数比小于某一定值，如0.00001 时，认为泛函已收敛，即：,(5-44),从以上论述可以看出： 引入拉格朗日乘子后，在假设初始速度场时，可以不满足体积 不可压缩条件，这对选择初始速度场有很大方便； (2)拉格朗日乘子有明确的物理意义，即收敛时的拉格朗日乘子就 是对应单元的静水压力。 (3)在变分运算中，假
18、设剪切屈服应力是常数。因此，在计算有加工硬化的材料时，每次取的移动量不能太大，特别对于硬化显著的材料要尽可能取较小的步长。 (4)在线性化中，采用了摄动法，并应用了牛顿（Newton）二项式展开，展开时假设u是小量，并略去了高阶微量，因此在计算中，每次的修正量要小，否则影响收敛性。,5.4 材料可压缩性法5.4.1 理论基础 刚塑性有限元法引入拉格朗日乘子后，可求得平均应力，但这样增加了许多未知数和方程数，计算量大大增加了。 从分析可知，一般刚塑性有限元法不能求解平均应力的原因在于屈服条件中没有考虑平均应力的影响。材料有可压缩的刚塑性有限元是假设屈服条件与平均应力有关，并写作为： (5-45)
19、 式中：g 为一个数值很小的正常数，一般取作0.01； 平均应力，,由式 5-45 可看出：当g=0 时，即为密赛斯屈服准则，屈服曲面在应力空间为一圆柱面；当g0 时，屈服面在应力空间为一椭球面，如下图所示。,屈服轨迹,密赛斯屈服条件下的势函数为： 同理，可压缩性材料的势函数为：,所以：,(5-46),(5-47),(5-48),又因：,体积变化速率为：,上式表明，可压缩性材料的体积变化速率不为零。,(5-49),(5-50),因为：,又因：,(5-52),(5-51),由式5-45 和式5-52 可得：,设等效应变速率为：,(5-53),(5-54),式5-53 可写作如下形式： (5-55
20、),将式5-56 代入式5-50 可得：,即：,(5-57),(5-56),将式5-56 代入式5-51 得：,(5-58),(5-59),上式写成矩阵形式为：,(5-61),(5-60),由以上可以看出，可压缩性材料的应力与应变之间的关系系数 与不可压缩材料有相同形式，区别在于对等效应力和等效应变有不同的定义，即：,(5-63),(5-62),5.4.2 系数 g 的取值 系数 g 的取值直接影响计算结果的精度，g 的取值可由实验来确定，表5-1 中给出了在单向压缩时，不同 g 值所对应的屈服应力近似值以及高度压缩率为10%时的体积变化。由表可看出，取g=0.01 时，式5-45与密赛斯屈服
21、条件相当接近，此时体积变化也很小。,5.4.3 求解方程的建立 对于变形体，可建立相同的泛函，即： (5-64) 为了计算方便，上式后面一项表示成矩阵形式为： (5-65),对于可压缩性材料发生塑性变形时，因塑性判据与平均应力有关，相当于隐含考虑了体积不可压。因此，对初始速度场的选择不需要严格满足体积不可压缩条件。,计算过程是先将变形体离散，划分成 M 个单元和 N 个节点。第 m 个单元的真实速度场的泛函为 ，假设位移模式的速度场所得的泛函值为,设:,则:,对泛函进行变分，当变分为零时，泛函取得驻值，这时的速度场是真实的速度场。,因 是任意的，所以有：,式中：k 为维数，对于三维问题，kN=
22、3N ，即可得3N 个方程。,(5-67),(5-66),(5-68),在单元内有: (5-69) (5-70) 将式5-54 代入式5-64 得： (5-71),上式写成矩阵形式为： (5-72),设: (5-73) 式中: (5-74),(5-76),(5-75),式5-72 又可写作： (5-77),式中 对式 5-77 求导得：,(5-78),将 按阶段硬化处理，将式 5-78 代入式 5-68 可得到一组方程，对于三维问题一共可得 3N 个方程，包含 3N 个未知数，即：,(5-79),所得的方程组式 5-79 是非线性的，不能直接求解，需进行线性化，采用牛顿-拉夫森 （Newton
23、-Raphson ）方法线性化，设第n+1次速度场 是第 n 次速度场增加一个 ，则有下列关系式： (5-80),将式 5-80 代入式 5-79 中，并在 处按台劳(Taylor)级数展开，略计高阶微小量，得到下列线性方程组，即：,由式5-81得到的线性方程组就能在假设一初始速度场的基础上进行反复迭代，最后收敛于真解。利用这种方法求出来的应力不再是应力偏量，而是全应力。,(5-81),对于轴对称问题，当单元采用等参四边形单元时，式5-81中的 有下列展开式：,式中：,(5-82),(5-83),(5-84),(5-85),式中：k2、k1、k、k+1为下标i、j、k、m 的代号，0为下标时表
24、示平均值。 B、C、D 分别由式2-28和式2-36中对应A 、B 、G 的给出。,5.5 罚函数法 5.5.1 求解方程的建立 刚塑性有限元法的一个基本假设是体积不变，罚函数法从这一点入手，引入一个很大的正数乘以体积应变速率的平方，即： (5-86),将此项添加到泛函中去，于是可将式 5-11 改写如下： (5-87),当 在每一点的变化率都接近零时，这个泛函就取得最小值，这时所对应的速度场就逼近真实速度场。可以按照处理拉格朗日乘子法同样的步骤导出罚函数法的计算公式。,罚函数法与拉格朗日乘子法的不同之处在最后一项，下面只对这一项进行分析。令这一项为 ，即：,将上式写成矩阵形式为：,其中因：
25、(5-90),(5-88),(5-89),将式 5-90 代入式 5-89 得：,对泛函变分后， 项为：,因 是一个数，故有 ，则,(5-91),(5-92),(5-93),是节点的速度，可将其从积分号中提出，即： (5-94),其中:,因泛函变分取得驻值时所建立的求解方程是速度的非线性函数，求解时需采用摄动法对方程进行线性化，其中体积变化率这一项为：,用式 5-95 代替代替式 5-40 中的 可得下式：,由前面推导有： (5-97),将式5-96 代入式5-97 得：,(5-96),(5-98),式中：,式 5-98 可写成如下的简单形式，即：,罚函数法对于三维问题一共只有 3 N 未知数
26、和 3N 个方程，这比拉格朗日乘子法要少 M 个方程和未知数，因此可节省内存和计算时间。这种方法还有一个特点就是收敛速度快，这给求解真实速度场提供了一个有效的办法。 罚函数法中的 是一个很大的正数，它的取值是否合理，直接影响到收敛速度。一般取 。 罚函数法的缺点是在求解时只能求出应力偏量，而平均应力不知道，为此不少人在这方面做了大量工作，下面介绍一种在罚函数法中求解平均应力的方法。,(5-99),利用已求出的应力偏量，可求出平均应力，推导如下。 应力张量可分解为应力偏量和球张量，即：,5.5.2 应力的求取,(5-100),式中： 为平均应力。,将式5-100 对坐标求偏导得：,(5-101)
27、,由平衡方程有：,由式 5-101 和式 5-102 得：,(5-102),(5-103),则:,因：,(5-104),(5-105),对上式两边积分,(5-106),得：,(5-107),有限元法将变形体划分成许多小单元，任一点的应力偏量可由它所在单元的各节点上的应力偏量来表示，即：,注意：在有限元计算中，一般都引进局部坐标系，对整体坐标微分时应采用复合微分，设局部坐标为Qj ，则：,(5-108),其中 表示 k 节点的应力偏量,(5-109),其中 为雅可比矩阵的逆矩阵。,又有 (4-110) 将式 5-109 和式 5-110 代入式 5-107 得： 这样，只要找到一个参考点，例如某
28、边界上一点的应力为已知，就可逐点积分来求出各点的平均应力值。,(4-111),5.6 三种应力计算方法的比较 拉格朗日乘子法由于引入了附加未知数 (乘子 )，使方程数增多和总刚度矩阵带宽增大，求解方程的计算时间要比另两种方法长；三种方法的总刚度矩阵均为对称的、稀疏的矩阵，但是拉格朗日乘子法的刚度阵不呈现带状，而罚函数法和修正罚函数法的刚度阵均呈带状。因此，后两种方法可节省贮存空间和求解方程时间。,刚度矩阵元素分布示意(a) 拉格朗日乘子法 (b) 罚函数法、修正的罚函数法,5.3.4 三种应力计算方法的比较 当初始速度场选取不当时，罚函数法和修正罚函数的方程组会成病态，此时 ， 或值很大，可能导致迭代过程不收敛。因此，罚函数法和修正罚函数法是对初始速度场敏感的。,
