弹塑性有限元法理论基础

1、弹塑性有 限元法理论基础金属橱料的塑性威型是一种相H复杂的禅塑性变形过程、它即包括儿何|:线 性、又包括物理IE线性，同时也具仃IE常复杂的边畀条件。仃限变形的禅塑性仃 限元法在求解金属塑性成型问题时能够比较正确地切合客观实际，能器有数的模 撩出金属塑性成型的全迁程对变电过程中的应力成变分布情况、质点的流动规 律以及卸载后残余应力、.应变的分们情况等都能够较好的与实际情况相符合。因 此、在金属塑性加一:领域中，弹塑性仃;艮兀法越来越受到研窕者的重视，利用禅 塑性仃限兀法来研究金属塑性.伐型规律的研窕报道也越来越条、获得的住果也越 来越显著。随着有限元技术和塑性力学理论的不断最展以及与计算机技术

2、的结合， 使得弹塑性仃琅元法的应用空间和应用前景更却广阁。22材料非线性本构关系22.1材料弹塑性行为的表述具仃弹塑性特性的材料在外力作用F发无塑性变形后，该部分塑性变形将在 外力卸载后保留卜米住为血久的变形、这说明这类材料在加裁刑卸载情况下的应 力应变之间存在圳线性关系到，如图1所湛O图21弹塑性材料应力和成亶关系示意图Fig.2.1 The re Latin juhip of gtre 拓-5tnkii of elastic-platie tnaletLaL对一般的硬化材料而育，都存在一个屈服应力知在低于应力的情况卜, 材料表现为弹性状态，这种状态下材料提外力作用I、发仁的变形是订逆的。咛

3、应力到达b祯后，此时材料则进入弹塑性状态，此后，材料发生的变形是不可逆的。 也就是说，材料在弹塑性状态下发生的塑性变形，在外力被卸载后将会被永久的 保留下来。对已发生过塑性变形的材料进行加载，材料将继续发生变形，同时应 力也随之增大，重新进入塑性变形的成力值将比初始应力值/。大，这就是所谓的 加工硬化现象。此时应力Z(0)为塑性应变、的函数，表示为b，=g(2.1)2.2.2塑性力学基本法则屈服准则对于各向同性材料，在某一应力状态下，当应力满足如下条件时材料发生屈 服，即F(气)=0(2.2)式中为应力张量分量：尸()=。的儿何意义代表了应力空间上的一个超曲而， 迪常被称为屈服而。实际应用中有

4、多种屈服准则，其中大多数金属材料采用的Von Mises屈服准则 可表示为：5 =房(2.3)式中为应力张量分量，为材料的初始屈服应力。强化准则强化法则用于描述屈服而随内变量的改变规律。LI前广泛采用的强化模型有 等向强化模型和随动强化模型。等向强化模型认为屈服而形状不变，只是作均匀的扩张，后继屈服而仅与一 个参数/有关，可表示为- f (a) - /(A:) = 0(2.4 )其中参数/是内变量的函数。随动强化模型则认为屈服而大小和形状均不变，仅是整体地在应力空间中作 平动，其后继屈服而可表示为=/(b-Qbp) = 0(2.5 )大多数材料的屈服而规律介于两者之间。如果应力空间中应力方向变

5、化不大, 等向强化与实际较符合。如果，当需要考虑循环载荷下能耗时，则需采用能够反 映包兴格(Bauschingcr)效应的随动强化模型。甘前，多数各向同性硬化材料采 用的硬化准则就是等向硬化准则o流动准则流动准则确定了塑性应变分量在塑性变化时的大小和方向。塑性理论认为材 料进入塑性状态后存在一个势函数g(b,%,k)=0,简称塑性势。当材料发生塑性变形时，其瞬时塑性增量可由使其发生塑性变形的势能推导出来，表示为d矿 2 皂(2.6)式中必是一个非负的尺度因子，d0表示加我，必=0表示其他情况， 舞气定义的向量方向与是沿着应力空间后继屈服而0的法线方向一致。加载与卸载法则对理想塑性材料来说，由于

6、屈服面的大小和形状不随内变量而改变，因此， 屈服而为厂同=0。当应力增量心增加时，有塑性变形山出现，此时，材料处于 塑性加我状态。加我时，材料从一个塑性状态变到另一个塑性状态，应力点在屈 服而上。如果& = (),反应是纯弹性的，此时材料处于塑性卸我状态。卸我时， 材料从塑性状态退回到弹性状态，应力点将离开屈服而。用公式表示理想塑性材料的加卸我准则为:0=0卸载加载(2.7)对于具有强化的塑性材料，除了加我。卸我外还存在一个中性变我状态，此时 d% =()，但应力点仍然在屈服而上。公式表示为(2.8)却栽中性变载加栽223应力应变增量关系具有各向同性硬化的弹塑性金属材料的本构方程是非线性的，本

7、构方程是应 变的函数，与材料本身的物理机械性能和加我历史均有关系。不论我们研究的是 材料的非线性问题还是儿何的非线性问题，或者两者的混合，经过有限元离散列 式后，最终都会归结为求解一个非线性代数方程组的解问题。在计算有限元方程 的解时，通常情况下是无法求得它们的精确解的。而是采用芥种数值离散方法来 代替真实条件，其实质是用一系列线性方程组的解去逼近所讨论非线性方程组的 解。通常情况下，首先会在微小的变形增量范用内对其进行线性化处理，然后再 通过迭代法进行求解，这使得在整个大变形过程中仍然保持它原有的非线性特性 65O弹塑性增虽本构方程的张量形式为:(2.9)23弹塑性加载卸载的处理方法在利用弹

8、塑性有限元法模拟金属塑性成型过程时，通常采用增量步长的加载 方式。在增量步长加我过程中，判断某一单元或单元体是否处于弹性状态、塑性 状态，还是处于弹性、塑性交织状态（即单元的某一部分发生弹性状态，另一部 分发生塑性状态），是通过检查单元内的高斯积分点上的应力、应变值的大小来确 的，并以此来确定加载和卸载过程。然后再对应力和应变项进行相应的调整，以 满足材料屈服准则和本构关系阿可。在计算某一单元内高斯积分点应力时会有以下情况：纯弹性加载。即在某一增量步载荷的加载前后，该单元内的高斯积分点均 为弹性，则该点的应力按照弹性规律来计算。塑性加载。即在某一增量步载荷的加载前后，该单元内的高斯积分点均为

9、塑性，则该点的应力按照塑性规律来计算。过渡点。即在某一增量步载荷的加载以前，该单元内的高斯积分点处于弹 性状态，而增量步载荷加载后，该点就进入了塑性状态。以下是具体计算高斯积分点应力的步骤：1）给变形休施加载荷增量AP作为笫次迭代的外加载荷，通过迭代求解有限 元方程组可以得到笫r次的应变增量4/和位移增量血L2）根据公式=DM得到应力增量的的大小。3）可用犬 5+奴累加可得到每个高斯点的应力,式中，式中表示第 -1次迭代后的高斯点应力。4）判断了 +其中为初始屈服应力槌为硬化系数，和7”分别为笫1此迭代后的等效应力和等效塑性应变。若检查的结果：步骤4）中不等式成立，说明该高斯点的应力在步增量载

10、荷加载之前巳经 超过材料的屈服应力即该点的材料己经屈服，再判断此时的等效塑性应力是否超 出了此前的等效应力，即判断淇了”是否成立。若成立，则表明该高斯点在此前 产生屈服后随着加载的继续应力还在继续增加，此时的高斯点应力己经脱离了屈 服而，超出矢量大小为；-了”，这时就必须通过调整缩减系数使其重新退回到图 2.2所示的屈服m （F=0）上。否则，高斯点处在弹性卸载阶段，转至步骤7）。步骤4）中不等式不成立，则说明步增量载荷施加之前高斯点没有达到屈 服状态，则需判断此时的等效塑性应力是否超出了材料的屈服应力。,，即判断 ；/是否成立。如果是，在施加本次迭代的载荷时，该高斯点己经出现屈服，超 出材料

11、屈服应力的那部分应力同样需要退回到屈服而（F=（）上，如图2.3所示。否 则，高斯点仍将处在弹性状态，转至步骤7）。图2.2弹塑性体内已屈服点上的应力增堡变化示意图Fig.2.2 Schematic diagram of lhe stress incremental change on the yield node in elaslic-plastic body图2.3弹塑性体内初始屈服点上的应力增量变化示意图Fig.2.3 Schematic diagram of lhe stress incremental changeon lhe initial yield node in elasli

12、c-plaslic bcxly对于巳经出现屈服的高斯点，只需计算满足屈服准则的区段+Q-R) db-对于已经产生屈服的单元，其应力点只能通过在屈服而上的移动来满足材 料的屈服准则和本构关系的条件。因此剩下的应力RM就必须等效的消除，从 而使应力终点A退回到屈服而上来。又图2.4可得：妃=以1&34(2.10)或者，b = b + do = b * + db： -dad”( 2.11)其中,儒回S儒如(2.12)(2.13)式中心为材料的塑性系数，心为材料的本构矩阵在单元弹塑性转化时的增量。式（2.11）给出了应力寸”产生增量时do：时总应力寸应满足的塑性条件。对 给定的一定大小的应力增量，可能

13、使对应的最终应力点。离开屈服面（F=0） o 当增量载荷取的足够小时，这种由此带来的偏离是可以忽略不计的。而且只需调 整矢量寸的大小，就可以使得D点退回到屈服而上来。图24应力在屈服而上的修正过#呈示氐图Fig.2.4 Ainendalory procedure restoring stress to the yield surface5(1设淇为应力/所对应的等效应力，当。位于屈服而上时有$ =/修正因子取（4 /；）/7,则进过修正后的退回屈服而上的应力为：(2.14)/ = b,（W吗b当增量载荷较大时，用上述方法引起的误差是较大的。如图2.4所示的成力点。和 D，将弹性应力按照应力这算

14、的方法退回到。点，然后再通过比例减少至屈服而上的D,很容易看出这样做的误差是很大的。为了提高计算的精度我们将超出的应力 部分成若干部分，此处分为3部分，然后依次对每个部分增量进行处理使其在三个 周期后退回到屈服而上，对应图中应力点矿 再通过适当的调整，就能得到最终应 力点E。从图2.4不难看出，两种处理方法得到的最终应力点E和D相差是很大的o 而后者得到的点相对于前者的D误差减少了。另一种修正方法是在退回的每个(2.17)(2.18)(2.19)(2.20)(2.21 )周期，都把应力点修正到屈服而上来。超出应力区段AB所划分的段数越多，就越 精确，但这也使得计算更为复杂和耗时。然而我们可以采

15、取一种折中的方法勺，即 将超出的部分应力RW：分成段，其中为整数且满足：rnw 8x（丝-%）+（2.15）匕,式中，X-/为超出的应力值，为初始屈服应力。7）利用公式=7*妃来计算高斯点应力的大小。2.4弹塑性有限元求解列式2.4.1虚功率方程按拉格朗日参数描述的虚功率方程可表示为：% 房,W =、职心 +2.16）式中，气为克希荷夫应力张量的分量：是虚格林应变速率张量的分量：P, 为作用在变形体外力作用而S上的表而力向量的分量：f；为物休单位休积力的分 量：$匕为物休内某质点的虚功速度的分量：V为变形体的休积。在外力作用表面S面以外的其余表而上，如果已知质点的速度分量为则可 得己知边界上的

16、速度边界条件为：_X = Vi2.4.2弹塑性问题的增量方程首先假设：材料为连续、均匀、各向同性材料。位移与应变仍为微小量，线性儿何关系成立。物体在平衡状态下，如果外力缓慢增加，其增力为： = M M，”7（体积力）P=此 M APzf（面积力）-相应地将产生位移、应变、应变增量为：&= 冬 S2” 编b=Ab、My m Ar 么增量应满足的系统方程为：1）儿何方程=囱叫2）平衡方程M =3）力边界条件（在Sp上）4）位移边界条件（在Su上）即问=印(2.22)在上而各式中,2WC?XO o o O 。-力 o Aaxo oab= f%，d%, 牛(2.23)回=七ny mJ为边界方向余弦2.

17、4.3弹塑性本构方程根据塑性流动理论，增量型的Prandll-Reuss（P - R）方程可写为:则:旗如4 匕标ki - f 1(1+|h(2.24)(2.25)式中，H=dq/d矿，n.为我荷性质判断因子。2.4.4弹塑性有限元格式及解法(2.26)利用虚功原理建立的弹塑性增量形式的有限元方程为：K0U=AP式中，K卜zJK1ZD为单元刚度矩阵组合为整体刚度矩阵： AK1 = J Bj DlBNU为单元刚度矩阵： AVAP =Z J AP,为整体节点裁荷增量： A AAP = J N了 ” dV + & J N了 AP； dS为单元节点裁荷增量。AVAS其中N为插值形函数：8|为成变形函数

18、：|切为弹塑性矩阵，在塑性区时它与当 前应力状态有关，因此，单元刚度矩阵可能与应力有关。对塑性区每次加载要重 新计算K1.称式（2.26）为变刚度矩阵方程，其计算过程如下：将外裁荷分为M步P“=P妇+印 m=123,.,M 2.27 )每步计算：K(a)-AUm=APm2.28)每步加载完成时计算下列各量Um=Um+AUm叽+L+W代2.29)同,” =”,+不断增加载荷直到加完为止。2.5接触问题有限元理论旋压件的成形主要琵旋轮进给时工作圆角部位对坯料施压使其产生逐点连续 塑性变形，从而获得满足一定形状的产品。由于旋轮进给过程中伴随着自身的旋 转运动，使得旋轮与坯料的接触过程变得极其复杂。其

19、接触问题是典型的边界非 线性问题。在利用有限元法解决边界非线性问题时，既可以采用全量迭代法求解, 也可采用增量方法处理。对接触条件的处理可采用直接法，引入接触条件并不断 迭代修正接触区域，直到最后满足所有边界条件为止：同时也可采用罚因子和数 学规划方法来处理接触条件。2.5.1增星接触边界条件设变形体A、B在外力作用下发生接触，接触区域为q,如图2.5 (a)所示。 在以上分布有接触力R，为了分析方便，先设置可能接触区，如图2.5 (b)所示。 建立增量接触条件。图2.5变形体A、B的接触示意图Fig.2.5 The sketch of contact between A and B设在f时刻

20、，接触体上己再我荷P作用，接触体表而局部方程为:(2.30)r1 少 g%*式中，土、与、勺为接触的局部更新坐标，它满足协调条件，即增加我荷AP后， 接触面的非嵌入条件为： +（2.31）式中，耳为表而位移增量。 若取B为参考，对上式进行展开，并取线性项后得：（/小：）计.+ 心 m2.32）式中， = 心AU:为增量位移向量：/,. = 等 芸 1,小（劄+（割 为接触外法向余弦：为更新的接触法向间隙。接触力的存在的条件为：(2.33)(2.34)r pb p B oI p； = - p：其中，1为接触力。当不考虑摩擦时，接触力为法向力。2.5.2增量接触有限元方程对A、B两接触体的分别建立

21、增量有限元方程如下： K： KtllL74l LF1 J 氐Ur* JUk； K曲1诚式中，AUJ为接触点的位移：Ai7为其余节点的位移：为节点外我荷：AR； 为可能接触节点上的接触力增量。在引入边界接触条件的情况下，即可对增量接触有限元方程进行求解了。2.6动态显式算法基础2.6.1动态显式分析概述利用有限元法求解答类非线性问题，最终都转化为求解非线性方程组的解的 问题。根据求解方程方法的不同，可以将有限元法分为动态显式有限元法和静态 隐式有限元法两种。隐式有限元法计算准确度高，能够较准确预测残余应力和回 弹变形，但其计算时间相对较长，对于处理具有复杂接触条件的三维金属成形成 型问题时较难实

22、现。然而，动态显式有限元法却能够很好的处理这类问题，可以采职多种方式对计算速度进行调节，大大节省了计算时间。动态显式有限元方法 常常被用于求解复杂接触问题、复杂曲而问题、高度非线性准静态问题、材料破 坏和失效等。本文中，旋压在加工锥-简形零件时，旋轮与坯料都在转动，旋轮通 过在一定的成型轨迹下的运动来实现工件的成形，是典型的动态分析，故本文分 析中采用动态分析有限元算法。2.6.2显式时间积分对于无阻尼系统而言，在某一时间段开始时，程序将利用式(2.35)来解动态 平衡方程，即利用节点质量矩阵M，乘以节点加速度】，等于节点的合力：Mu = P-F(2.35)式中，P为作用于节点上的外力：尸为单

23、元内力。设当前时刻为f，则该时刻的加速度为：u =M(PF)L(2.36)由于显式算法通常采用块状的或是一个对角的质量矩阵，这使得加速度的求 解比较简单，不需要求解联立方程。每一个节点的加速度完均由作用于节点上的 合力和节点的质量决定，使节点计算变得较为容易。(2.37)节点的加速度是通过中心差分法的时间积分得到的，即通过假定加速度为不 变的常量，来求解速度的变化，再利用求得的速度的变化值加上一个时间段中节 点的速度，以此来确定当前时间段的中心速度。用公式表示如下：U = U + :U2时间段结束时的位移等于该段时间内节点速度沿时间积分的值再加上该段时 间以前的位移。公式表达为：W +广y(2.38)这样，在研究的时间段开始时，提供了满足动力学平衡条件的加速度。得到以下是显式了加速度，即可