用弹塑性有限元素法解平压头下的压入问题.pdf_大学文库

D0I:10.13374/j.issm1001-053x.1982.01.015 北京铜铁学晓学报 1982年第1期用弹塑性有限元素法解平压头下的压人问题力学教研室乔端钱仁根捕要本文采用弹塑性有限元素法中的变刚度法和加权平均弹塑性矩阵，解光滑压头下平面应变的弹塑性硬化材料的压入问题。得到压力一一位移曲线，屈服压力和塑性区的边界，并与理想刚塑性的滑移线场理论结果相比较。同时计算了在平面应变条件下，工业纯铝压缩时的冷变形抗力。计算结果与试验结果基本相符。一、引论平面应变条件下的理想刚塑性材料，在光滑平压头下的压入问题，可用滑移线场理论求解。被压材料的厚度(2h)和压头宽(2a)之比>1时（图1）其解由R.Hi11得到〔1〕，当h/a<1的情况，由A.P.Green解压头得〔2)。但平面应变条件下，硬化材料在光滑平压头下的压入问题，由于其弹塑性交界面比较复杂，用弹塑性理论求其解析解是相当困难的。但运用弹塑性有限元素法可求得数值解。本文中利用弹塑性有限元素法中 2a- 的变刚度法，解平面应变条件下的弹塑性硬化材料，在光滑平压头下的压图1示意图入问题。得到压头的压下位移和作用于压头上总压力之间的关系、屈服压力以及塑性区随压头位移的增加而变化的情况。并与相应的滑移线场理论的结果相比较。利用压入问题的结果，可计算金属材料平面应变条件下的冷变形抗力。只奖在压入问题计算中，令h/a=1,此时的平均屈服压力为冷变形抗力〔3)。冷变形抗力在金属冷加上成形中是一个很重要的工艺参数。一般是用拉伸试验或“平面压缩”试验测定，或把拉伸屈服应力乘上1.155为平面应变条件下的变形抗力。但由于材料的加工硬化和正交异性等因素，所以把屈服应力的1.155倍作为平面应变时的变形抗力，显然不是完全正确的。另一方而，常 *本文1980年12月27日收到。 170

北京栩铁学映学报年第期用弹塑性有限元素法解平压头下的压人问题力学教研室弄端钱仁根摘要本文采用弹塑性有限元素法中的变刚度法和加权平均弹塑性矩阵，解光滑压头下平面应变的弹塑性硬化材料的压入问题。得到压力— 位移曲线，屈服压力和塑性区的边界，并与理想刚塑性的滑移线场理论结果相比较。同时计算了在平面应变条件下，工业纯铝压缩时的冷变形抗力。计算结果与试验结果基本相符。一、引论平面应变条件下的理想刚塑性材料，在光滑平压头下的压入问题，可用滑移线场理论求解。被压材料的厚度和压头宽之比》时图其解由得到〔〕，当的情况，由解得。但平面应变条件下，硬化材料在光滑平压头下的压入问题，由于其弹塑性交界面比较复杂，用弹塑性理论求其解析解是相当困难的。但运用弹塑性有限元素法可求得数值解。本文中利用弹塑性有限元素法中的变刚度法，解平面应变条件下的弹塑性硬化材料，在光滑平压头下的压压头阮，声，各厂 ‘ 一，气令习图示意图入问题。得到压头的压下位移和作用于压头上总压力之间的关系、屈服压力以及塑性区随压头位移的增加而变化的情况。并与相应的滑移线场理论的结果相比较。利用压入问题的结果，一可计算金属材料平面应变条件下的冷变形抗力。只要在压入问题计算中，令，此时的平均屈服压力即为冷变形抗力〔。冷变形抗力在金属冷加工成形中是一个很重要的工艺参数。一般是用拉伸试验或 “ 平面压缩 ” 试验测定，或把拉伸屈服应力乘上，为平面应变条件下的变形抗力。但由于材料的加工硬化和正交异性等因素，所以把屈服应力的倍作为平面应变时的变形抗力，显然不是完全正确的。另一方面，常本文年月日收到。 DOI ：10．13374／j ．issn1001—053x．1982．01．015

用的拉伸试验和“平面压缩”试验，都有一些难于克服的困难。如“平面压缩”试验，不能完全克服介于压头与试样间的摩擦，而且由于被压材料的厚度在变化，所以难于做到 h/a=1。而这些在有限元素法的计算中是很容易做到的。本文计算了平面应变条件下，不同积累压下率工业纯铝的冷变形抗力，而且分有外区 (材料长度比压头宽度大许多)和无外区（材料长度等于压头宽度）二种情况进行计算。计算结果与试验结果基本符合。二、基本原理由于压入问题的对称性，所以只取被压材料的1/4研究。取x,y,z轴，z为被压材料的宽度方向。因为其宽度较大，可以认为是平面应变问题，计算单位宽度。利用三角形元素把材料离散化。对每个元素，设其节点的位移列阵为{ò}°，其增量为 △{ò}°，节点载荷列阵为{F}·,其增量为△{F}°，应力列阵为{σ}，其增量为△{o},应变列阵为{ε}，其增量为△{}，联系节点位移与应变的矩阵为〔B),对于处于坪性状态的元素，其弹性矩阵为〔D。),元素刚度矩阵为〔k。),对于处于屈服状态的元素，其弹塑性矩阵为(D。p),元素刚度矩阵为〔ke)。由变刚度法可知：〔k。)=〔B)T〔D。)〔B)t·A (1,(a)) 〔k。p)=〔B)T(D。p)(B)·t·A (1,(b)) 其中：t为元素的厚度，在平面应变时取t=1,A为元素的面积。对于弹性元素， (k。)△{8}°=△{F}°， (2,(a)) 对于塑性元素，〔kep)△{8}°=△{F} (2,(b)) 现在的〔B)和〔D。)与弹性有限元素法中的完全相同。平面应变条件下的〔D。。)可由空间应力状态的〔De)〔4〕刘去与△z,△Txz,△TYz相应的行和列，得到〔5) 1-v 对称 1-2v (D。p)= (3) 其中：σx',oy'为应力偏量，v为泊松比，G为剪切弹性模量， s=号(1+)， (4) 等效应力 (5) H'为硬化曲线上一点的斜率，可由材料的应力一一应变曲线求得〔6)。 171

用的拉伸试验和 “ 平面压缩” 试验，都有一些难于克服的困难。如 “ 平面压缩 ” 试验，不能完全克服介于压头与试样间的摩擦，而且由于被压材料的厚度在变化，所以难于做到。而这些在有限元素法的计算中是很容易做到的。本文计算了平面应变条件下，不同积累压下率工业纯铝的冷变形抗力，而且分有外区材料长度比压头宽度大许多和无外区材料长度等于压头宽度二种，’ 况进行计算。计算结果与试验结果基本符合。二、器本原理由于压入问题的对称性，所以只取被压材料的研究。取，，轴，为被压材料的宽度方向。因为其宽度较大，可以认为是平面应变问题，计算单位宽度。利用三角形元素把材料离散化。对每个元素，设其节点的位移列阵为乙卜，其增量为 △笼节点载荷列阵为，其增量为 △王应力列阵为，其增量为△ 应变列阵为。，其增量为△ 。联系节点位移与应变的矩阵为〔〕，对于处于弹性状态的元素，其弹性矩阵为〔。〕，元素刚度矩阵为〔〕，对于处于屈服状态的元素，其弹塑性矩阵为〔。〕，元素刚度矩阵为。〕。由变刚度法可知。〕〔〕〔。〕〔〕一一，〔。〕〔〕。〕〔〕一，其中为元素的厚度，在平面应变时取，为元素的面积。对于弹性元素，〔〕△魂各 △ ，，对于塑性元素，〔，〕△ 乙 △ ，现在的〔〕和〔。〕与弹性有限元素法中的完全相同。平面应变条件下的〔。。〕可由空间应力状态的〔。〕〔〕划去与△ ， △ ， △ 丫丫相应的行和列，得到〔〕一一对称一，，，，一、一川一一一，， “ 一一二一一旦兰三兰竺一、一丫竺里卫一一一三兰竺、艺一了、产、少一 … 其中厂， ’ 为应力偏量，为泊松比，为剪切弹性模量，笋 · 斋等效应力万〔多一，，，、，〕 ‘ 为硬化曲线上一点的斜率，可由材料的应力— 应变曲线求得〔〕

显然，〔k:)与应力有关。在运算时，可用经第i一1次运算后得到的应力组成(ke), 进行第i次运算，余此类推。若把每个元素的刚度矩阵集合成结构刚度矩阵〔K),则最终的基本方是为： {K}△{6}=△{F} (6) 其中：△{6}，△{F}为整个结构的节点位移增量和节点截荷增量。利用边界条件，解方 (6)得△{8}和△{F}中的术知量。由〔B)△{δ}°=△{e} (7) 和对弹性元素，〔De)△{e}=A{o}, (8,(a)) 对塑性元素， (Dep)△{}=△{o}, (8,()) 可解得每个元素的应变增量和应力增量。若设第i次运算得到的应力增量为△{σ}，i一1次运算后得到的应力为{o}:-1,i次运算后的应力为{o},则， {o},={o}-1+△{o}:o (9) 塑性区旁边的元素，在加载过程中（不是在加载后）有可能进入屈服，这一类尤素称为过渡元素。对于它们用(D。)和〔De)进行计算均不妥，会带来不小的误差。因此，我们采用加权平均弹塑性矩阵〔D。p)=m〔De)-(1-m)(Dep)i (10) 其中：m为加权系数，它是过渡元素为达到屈服所需要的等效应交：量这次加载所引的等效应变增量的比值，显然0<m<1。计算时，迭代~3次即可得到比较止：的位。若元素不属于过渡元素，则处于弹性变形状态时，可取m=1,若处于亚性次，态时则取m=0。这样，(1，(a)),(i,(b),(2,(a),(2,(b),(3,(a),(3,(b)均可分别合与j为 (k)=〔B)T(De。)〔B)·A (11) (E)△{8}°=△{F?° (12) (D。)△{e}=△{o} (13) 三、计算程序计算？序的摇图如图2。下面按框图的次序逐一说明如下： 1.输入数据。输入压头宽度，被压材料的长度、厚度，长度方向和厚度方向划分元素所分的段数，压头下划分元素的段数，被压材料的力学性质，如应力一应变曲线等。 2.自动划分元素及节点编号，计算节点座标。元素按图3所示的形式刻分，元素中节点是节点序号中的最小者，按逆时针转向排列i,j,m。 3.给定压头每次压下的位移。由下面的计算结果可知，压下位移不同，其结果也不同。若每次压下位移较大，则同时屈服的元素较多，这样计算结果误差较大。为了得到较好的结果，要求每次压下位移越小越好，但这样势必会延长计算时间。我们在计算时一般采用次位移增量为1×10~‘毫米或1×10-5毫米。 4,计算元素刚度矩阵并集合成结构刚度矩阵。计算时利用加权平均弹塑性矩阵，并当元素处于弹性时，取m=1,处于塑性时m=0,若为过渡元素，则取计算所得的m值。 5.边界条件。对于上部边界中，压头下的节点，其垂直位移等于每次压下的位移值， 172

显然，〔。，〕与应力有关。在运算时，可用经第一次运算后得到的应力组成〔。，〕，进行第次运算，余此类推。若把每个元素的刚度矩阵集合成结构刚度矩阵〔〕，则最终的基本方刃为 △ 乙 △ 其中 △ ， △ 为整个结构的节点位移增量和节点截荷增量。利用边界条件，解方门得 △ 和 △福中的未知量。由〔〕△蓬各 △ 。和对弹性元素，〔。〕△王。么，，对塑比元素，〔〕△ 。 △ ，，飞可解得每个元素的应变增量和应力增量。若设第次运算得到的应力增量为△ ，，一次运算后得到的应力为卜一，次运算后的应力为 ‘ ，则，灌 ‘ ‘ ， △ ‘ 。塑性区旁边的元素，在加载过涅中不是在加载后有可能进入屈服，这一类元素称为过渡元素。对于它们用〔。〕和〔。，〕进行计算均不妥，会带来不小的误差。因此，我 ’ 采用加权平均弹塑性矩阵〔〕〔〕一一〔〕、。其中为加权系数，它是过渡元素为达到屈服所需要的等效应变加童这次加载所引起的等效应变增量的比位，显然。丈。计算时，迭代。一次即可得列比校止你， ’ 他。若元素不属于过渡元素，则处于弹性变形状态时，可取，若处于些性伏态时则取二。这样，，，，，，，，，，咬。均可分别合与少骊卜〔〕。〕〔〕 · ‘ · 〔石〕△ 二 △ 少 ” 〔石〕△ 。 △ 三、计算程序计算陛序的框图如图。下面按框图的次序逐一说明如下输入数据。输入压头宽度，被压材料的长度、厚度，长度方向和厚度方向划分元素所分的段数，压头下划分元素的段数，被压材料的力学性质，如应力— 应变曲线等。自动划分元素及节点编号，计算节点座标。元素按图所示的形式划分，元素中节点是节点序号中的最小者，按逆时针转向排列，，。给定压头每次压下的位移。由下面的计算结果可知，压下位移不同，其结果也不同。若每次压下位移较大，则同时屈服的元素较多，这样计算结果误差较大。为了得到较好的结果，要求每次压下位移越小越好，但这样势必会延长计算时间。我们在计算时一般采用每次位移增量为一 ‘ 毫米或丁毫米。计算元素刚度矩阵并集合成结构刚度矩阵。计算时利用加权平均弹塑性矩阵，并当元素处于弹性时，取，处于塑性时。，若为过渡元素，则取计算所得的位。边界条件。对于上部边界中，压头下的节点，其垂直位移等于每次压下的位移值

水平节，点力等丁零（无摩擦）。压头以外的节点力为零。对于下部边界，因为对称，所以开始水平节点力为零，垂直位移为零。右边界为节数据点力为零。左边界，因为对称，所以垂直节点自动划分元素，节点力为零，水平位移为零。编号，计算节点座标 6.解方。采用〔K)阵一维紧缩存给定正头每次压下的位移储，非用追赶法求解方程。 PP:=1 7.计算应变增量和加权系数。为得到较好的m值，迭代2次。 YES 8.计算△{σ}，△{F}等。由计算所得的 No △{e},利用(13)，(9)和(5)式得 QQ:=1 A{o},{o},。并利用△{8}、〔K)和(6) 让饰元素刚度矩碎式得A{F}。则第i次计算所得的总节点力为 {},={}:-1+△{F}:。由节点力就可以求集合成结构刚度加阵得压头总的压力，它除以压头面积，即得平均运用边界条件压力。解方程，得△{8 9.输出结果。输出每个元素的m值，用它来判别元素屈服与否。输出元素的应力，等 i计算A{c}.m 效应力，节点力，压头压力，平均压力和节 Q- No QQ:QQ+ 点位移等。 YES 10.压头每向下压一个位移，必须计算计算△{o,{o}, 一次，PP丧示计算的次数。听PP=1时，可可A{F},F) 把整个结构作为弹性体计算，得到每个元素的 E头正力半均团等效应力g,确定·max,令0max/0,=L, 则把计算所得的值，除以L,即可得到被压材料进入屈服时的位移，应力，应变节点力等。 No 全屈服> ppi=pp+1 11,片刚性区与边界相连，即达到全屈 YES 服状态时，此时的平均压力即为平面应变的变停机形抗力。图2计算程序框图四、压入问题计算结果本文计流了已具方106潮件交形的工业纯铝的E入问题。其屈服极限σ。=6.17i公斤/ 毫米2，啪性模屉E=7201.6公斤/毫米°，泊松比v=0.33。它的应力一应变关系如图8中曲线I所示。计算长度1=5mm。取以下几种滑况进行计第。 1) 0=1毫米， h=2毫米， ,'h=1 2。 2) n=1毫米 h=1.5毫米， a/h=2 ● 3 3) a=1毫米 h=1毫米 ,a/h=1。 4) a=1毫米， h=0.5毫米，a/h=2。 173

水平节点力件千零无摩擦。压头以外的竹点力为零。对于下部边界，因为对称，所以水平节点力为零，垂直位移为零。右边界为节点力与零。左边界，因为对称，所以垂直方点力为零，水平位移为零。解方程。采用〔〕矩阵一维紧缩存储，并用追赶法求解方程、卜算应变增。和加权系数。为得到较好的故，迭代次。计算 △ ， △ 等。由计算所得的。卜，利用，和式得，，。并利用 △ 乙、〔〕和式得州卜。则第次计算所得的总节点力为魂，王 ‘ 一， △ ‘ 。由节点力就可以求得仄失总的力，它除以压头面积，即得平均压力。输出结果。输出每个元素的值，用它来判别元素屈服与否。输出元素的应力，等效应力，节点力，压头压力，平均压力和节点位移等。压少协向下压一个位移，必须计算一次，友示计算的次数。当二时，可把整个结构作与弹性体计算，得到每个元素的等效应力，确定。、，令。、、，则把卜算所得的值，除以，即可得到被压材料进入屈服时的位移，应力，应变和节点力等。当塑性区与边界相连，即达到全屈服状态时，此时的平均压力即为平面应变的变形抗力。自动划分元素，节点编号，计算节点座标计算△ 魔》，笼，下八卜王卜图计算程序框图四、压入问题计算结果术之计牛已具汀丫塑性呀形工业纯铝的压入问题。其屈服极限了公斤毫米 “ ，弹性模鼠落公斤毫米巳，泊松比二。。它的应力— 应变关系如图中线所示。计算长度。取以下几种洁况进子一八一算。一一一 ‘ ‘ 川了比且幻幻勺 ‘ 二毫米毫米 ‘‘ 毫米二毫米毫米毫米，毫米毫米勺性八月、、产产、产、声

二毫米，毫米，毫米，毫米，其元素划分如图所示。当压头每次压下位移为一 ‘ 毫米时。其位移与塑性区的关系如图所示。求得全屈时的平均压力与屈服应力的比值。并与理想刚塑性材料时的结果比较列表于下。表一‘︸，口一︸了︸‘一一二弹塑性有限元素法解。， · 二 · 二，滑移线场理论解自 ‘ 曰口用翻司口甲以加曰州口口，口 ‘ 目州州目， ‘ 门铃由文献〔〕〔〕提供约曲线而得、声 ‘ 尸，户户产产，，魂夕夕价口口口夕曰叼口叼匕团匕匕团匕匕匕四八，川侧口口口口口口口口口区叼口因团匕匕匕凶勺一月夕夕口口门口门口口口口叼国叼叼叼叼叼区团曰网口汤阴口口夕口口口口口口叼叼区匕匕匕夕口夕口口门口口口口口口侧例叼团口匕匕匕国冈口口圈口例口口胭口口例例冈团团叼叼团团团口圈侧团团口口叼夕口曰叉叼团叼叼叼叼日团团日团团叼团团团叼团叼叼团叼团区叼叼户、、二，一，五 ” 。了、。了。，孚一 “ 一几，，，，矛生，沙夕节夕口夕口口区区区匕匕匕匕匕瞬曰夕口口口曰口口口口口区口斌区忆区匕区以团一气刁口夕口口口口口口口口口口泛比忆区叼匕匕团沉夕沙口伊使泛泛口口口口口区健区匕匕匕匕叹夕任「任反仅沙汤口口口口口口反区区区区叼区区口区叼叼叼叼口叼区叼叼团团叼几少几二川，元，赦一，少，乃，飞 ’ 二之于亡丰从少七准， ’ ，’ 点故为 ‘了。叮，下玫 ‘ ，。因因口口团口口区口口因团口区泛口口区口冈冈比又口口口口少口口口口口口仁习口夕口口夕誉口口夕甲之夕区口口口口口夕口口区口口叼回夕叼叼叼回叼团口叼团叼冈团叼冈团，二，五二『， △ 二。夕夕夕夕夕夕口口夕夕夕口夕夕夕夕夕夕夕夕夕夕夕夕夕夕口夕夕夕夕夕夕尸冈夕冈冈夕夕，，」例夕口曰口口口夕口夕口口口夕口口口口夕夕夕口口夕夕夕少夕夕口口夕夕刁口训曰曰口口网夕国夕日夕夕习夕夕夕夕口夕夕口口日夕夕口厕回夕口口口夕口口口夕夕叼国夕口口国夕目口口曰叼口巨叼国国口口曰曰冈曰夕叼叼叼团团曰团团口口曰夕夕口叼团叼口口口口叼曰，‘ 了丫兮夕夕匕，，，，，二匕匕匕匕匕匕团回叼团叼口口口夕夕夕夕夕夕夕口夕夕夕口口口口夕口口口口口口口口夕口，一曰四匕团匕匕团囚团团团口口口口口夕夕夕夕夕口圈夕网夕侧夕夕夕夕沙夕夕防沙目台目目目，四叼团团团因叼团冈团口冈冈冈叼网口网网夕夕冈侧夕网叼回夕夕夕夕夕夕口口夕圈阴例冈冈口口口口冈口团口冈冈口夕口冈口口冈叼冈曰网网冈网冈夕夕夕夕夕夕夕曰夕夕夕国附 · ，，，、。认，，，，，，，，，，，，之巨团叼区口口口口口夕口夕夕夕口夕曰夕口口夕口口口口门口口口口户户户尸，一，门阵门，减叼叼叼叼因口团口口口网夕夕夕口口叼侧夕夕例夕口夕侧臼份目目勺万州卜二州‘ 才酬口夕冈冈口冈冈夕夕叼冈夕夕夕夕冈夕曰夕夕夕洲洲川夕夕 ’ ’ 洲万万一丁刊仁一门气尸一一侈卞夕叼冈冈冈冈网冈曰团冈叼曰冈冈网曰目叼冈冈冈日口叼冈叼曰阅阴夕口二弓阵刁，五一，，。 “ ‘ 一、 · 元不从，。，图元素划分图