有限元方法的应力后处理技术_节点应力计算与等值线绘制_李明瑞

第卷第期年月机械强度宝公比。有限元方法的应力后处理技术—节点应力计算与等值线绘制李明瑞鲍妍光北京农业工程大学扭共本文介绍几种快速高精度的将弹性与塑性有限元解表示成节点应力的计算方法和绘制等值线的方法,可绘制正应力、切应力、应力差应力和、等效应力、主应力等种等应力迹线,不仅可用以模拟光弹和弥补光弹之不足,并可表示塑性有限元计算中塑性区的分布及其发展。有限元计算得出的应力值一般均表示为在单元内若干高斯积分点上的值一维的杆、梁单元除外,但是高斯点的位置在有限元剖分中并不直接给出,人们往往更希望能知道在剖分节点上的应力值目前还缺乏一套良好的节点应力计算方法其次,工程技术人员不仅关心在个别点上的应力值,如在何处有最大应力等,更关心整体的应力分布情况,以便为强度合理设计提供依据,如选择最佳结构形状及配置加强筋板等目前,一般采用光弹方法,但其缺点是要用异质材料进行模型试验,不能直接测试出主应力,试验周期长等本文提出若干种节点应力计算方法基于对节点二,口,,二,的计算可转换成节点,,。,,。一,,。的计算其中侧二,一二,二,生是等效应力并进一步将节点应力插值到单元域内,给出了一套快速的以单元为主体的等值线绘制方法‘比起传统的逐点追踪‘二州划等值线方法具有速度高、不漏划等优点用这套方法编制的有限元软件九如切加呱城加全砚可任选绘制等值线的类型和条纹级数,使用灵活方硬,不仅可与光弹试验结果对比,还由于其灵活性和等值线类型的多样化,大大丰富和弥补了光弹实验之木足一节点应力的全域小二桑法遇近在用有限元方法求解时,一般来说单元间位移场是连续的,但其应力场却是不连续的,因此只能直接求出单元内的应力而得不到节点应力为求节点应力值必须采用某种间接方法,由于方法不同其结果也会有所差异设令为用有限元方法计算得出的各单元内的应力,口可分别是二,。,,,二,,是欲求得的节点应力,是单元的位移插值函数另定义单元内的假设应力,它可由节点应力。及插值函数表示为少刃一般而言,计算应力口与假设应力是不相同的。现在提出这样一条原则以确定口。选取气使口与的误差平方在全域内的积分值达到最小,即,乏五。界。,’“卜将代入后并将化为一组等价的五口,二。则可得条件尧亦即有艺一。。口“‘一子玉。‘口上式是一组仁厂为未知量的方程组,未知量个数等于节点总数,方程组的系数矩阵类似于一致质量矩阵,可由组集方法求得方程组的右端项对不同应力项二,,,二,是不同的,因此求解相当于求解有三个不同右端项的线性方程组问题在求得各节点的二,“,,二,的基础上可由此而构造各节点的。,,,,一。,。,。一办’瓦一口二二了共中一可与光弹实验得出的等应力差互相比较二节点应力的单元最小二乘法加权平均以上方法确可得出比较满意的节点应力值,但要求解一组较大的线性方程组为此可提出另一条原则选取。使与。的误差的平方在单元域内积分值达到最小,即上面的第二方案仍需求解多个小型线性方程组也可直接求出节点应力。对于位移法求解时可先求出应变再算出应力,这种方法的精度稍低些对于采用假设应力的混合元或杂交元而言,采用直接计算单元节点应力然后再用加权平均法来求出全域的节点应力的方案还是可取的。四塑性状态时的节点应力计算以上方案均基于节点应力与单元域内假设应力的双线性插值关系,因此只适用于线弹性计算当单元进入塑性状态时由于单元域内应力与节点应力不再满足双线性插值,必须对以上算法进行修正以同时适合弹塑性状态的节点应力计算首先仍用以上方案技线弹性关系计算出节点应力然后对每一节点的各应力分量二,,,,二,对照屈服准则进行验算,如果超过则可近似地将各应力分量按比例。,缩回其中一训砚不千,盆一二,二,是等效应力,,则是屈服极限。于是单元节点应力的等效应力就正好处于屈服状态以上假设隐含地引进了比例加载假设,但当采用不太大的加载步沃时,对每一给定的加载步长而言,这条假设是允许的,从计算实践来说其效果也是良好的、少,、,矛‘、‘又。一丁。万一,’““一‘。亦即相当于对每一单元求解一组低阶方程叉。二“云一“由此可求得每个单元的节点应力,但是在同一节点处的应力将随单元而异为此可对不同单元的节点应力采取加权平均法以求得总体的节点应力见七五等值线的挟速绘制当各节点的应力值求得后,采兀按单元扫描各级等值线的思想〔,并在此基础上对四乏砰’。‘节点与多节点单元提出了高精度快速绘制等值线的方法如图所示的四节点单元,以。,吕二二二一一一二牛宁二二二,一乙川’式中牙‘是可选取的权因子,一般可选取为对四节点单元而言,如采用双线性插值函数,则有限元计算精度达到。,其中是网格剖分的测度为此可取平‘扭‘,’是单元的面积二因为可以设想,面积较小的单元应力应该比面积较大的单元应力更精确些,权因子应更大些,而月‘则正比于本方案的计算工作量与存贮量均小于第一方案,其计算精度图‘表示四个节点的坐标,以。表示四个节点的应力在等参单元设计中用同一插值函数表示位移与坐标变换,现假设这套场值公式对应力也适用,特别对于雪刀的中也是满意的〔三单元节点应力的直雄计苏与加权平均心点。处,其坐标与应力值分别为。艺,八。艺,八。艺,八然后将该单元魁分成图示的四个三角形分别考虑对每一三角形的等值线画法以三角形。”为例,设令。为某一指定的应力值,可以断定当。一一一一一时,在边间必无值口。的等值线通过这是因为采用了四节点等参元的双线性插值函数,排斥了多次通过单元域内任一线段的可能性当不等式不成立时,则必有一条。等值线通过线段,其中也包括了该等值线通过。点或点的可能性等值线与线段交点的坐标可用线性擂值公式近似地表示为二二。一。一口。一口。。一口二。一犷。一。。一。,。一分然后再逐一判别。等值线与线段或。线段是否有交点由等值线的连续性可以断言只要有等值线进入一三角形子域,则还必定要走出该子域。由此可以肯定如果口等值线与三角形某一边相交则必定而且只与该三角形的另一边相交反之,如等值线与三角形的其中两边不相交,则必定不与第三边相交以上这些论断可应用于程序设计中以提高运行速度顺次对四节点单元的四个子三角形扫描一遍后就完成了该单元的某一等值线的绘制由此可得出一套绘制有限元等值线的算法求出各节点应力对所有单元扫描完成以下各步确定单元中心坐标与应力值对所有需要绘制的各级等应力线值扫描对每一单元的四个子三角形扫描,以直线段联接绘制出该段等值线。为了提高绘制精度曾有人建议将一个单元划分成多个子单元’,例如将一个四节点单元划分成个或‘个四边形子单元笼·笔者认为这种做法一般不妥,等值线绘制精度必须与有限元计算精度相配匹才是有意义的如果原单元是多节点、高精度插值单元时,才有可能进一步封分原单元成多个子域例如对于节点单元采用双二次插值时,则可合理地将原单元劫分成四个子四边形单元,而对普通的四节点单元再作进一步划分则是多余和不合理的。为简单起见上面介绍的等值线绘制法都是在单元的子三角形内,用直线段联接三角形两边与等值线的交点而成为提高绘制精度,可以在单元内部将所需要的点先予以暂时存贮,而后用曲线拟合法连续绘出不难看出,对于点进行拟合,从而大大提高计算速度,而使得在微机上绘制等值线成为现实可能。,八火、犯儿、、‘图选取子三角形中与该等值线无交点的第三边中点,联接线段,同样用线性插值方法求得等值线与之交点于是可用文〕所介绍的用已知三点画曲线的方法得到比较光滑的等值线应该指出,应力等值线的绘制精度主要取决于有限元计算精度,亦即网格剖分的精度,当有限元网格剖分较密时,即使采用最简单的直线段绘制法也是令人满意的六弹组性状态时妞性区城的绘制法在有限元计算中通常都是以单元内的高斯积分点处的等效应力值来判定该点是否已进入塑性状态在图形中虽然可用不同记号来表示高斯点的应力状态孰为弹性与塑性,但是这种以分散点来表示一个区域的应力状态的方法既不直观也不够正确利用以上这种灵活绘制等值线的方法,只要指定在屈服极限附近画出数较多的等效应力等值线,就可相当精确地勾画出结构的塑性区域在不同的加载等级下分别绘制这些塑性区域图,就可形象地表示塑性区域的发展情况出几个算例例一个内外径之比为,的圆环,上下对径受压这是一个典型的光弹实验例题,它的结果可在许多光弹实验书籍中找到在用有限元计算时利用对称’性只计算四分之一区域,按。网格剖分等值线的各级数值都‘可自动计算、打印及标志在图纸上为节省篇幅,本文只附上用打印机拷屏的图象图至图分别给出了计算绘制的图象与文中所给出的图象的对比。其中节点应力按单元最小二乘加权平均方法计算圆环在对径压缩下的理论等厚应力条纹图(a)电木圆环在径向集中载荷。2公斤作用下的光弹性应力条纹图。内直径2r二19.5斑米;夕值径Do3.毫米;厚度t二4.5毫米;模型的拉伸(或压缩)条纹值ZF=韶.T公斤/厘米君。(b)径向受压团环的等和线图5(七)径向受压圆环的等差线图4‘.)早里权肺压维跳目环的等翻线(左半)’左固右铰梁中央受载P二1。。。.时一。:冰,的称口I线(七)径向受压圈环的等倾线图6例2:一根简支梁,长宽比为6:l,在上部中心处受集中载荷,由对称性只计算一半区域,有限元网格剖分为20xlo.节点应力按全域最小二乘法计算.图7给出了拷屏所得的等差线图与光弹实验照片,1的对比.〔c)左固右锥架中央受旅P二10000时扭性度分布圈吕笋r一(.)左固右铰弟中央受徽护二12。。晰.知,的等。I线(a)跨长B借于高度的跨中受载愉支梁的应力条纹圈(b)左固右饮梁中央受城P二1200。时一。,之。,的娜。:线肾肾纂薰薰(b)简支梁的计算等差线图7例3:一根左端固支、右端饺支的梁,上部中点受集中载荷.分三次加载,载荷分别是10000,120,o,14000,其他参数是E=Zxlo.,,==0·3,长160,宽30,口,=2000(各种单位是相容的).图8、9、10给出了三种载荷情况下的大拉应力区(口:》。,),大压应力区(一。:口,)和塑性区a:),,的分市区域图.节点应力按单元节点应力直接计算并进行加权平均处理.当节点应力口:口,时按(四)中所述进行处理。(c)左固右锐璀中央受载护二任000时蟹性区分布圈9(.)左固右铰梁中央受载P,14000时。:沙,的等。I线(b)左固右铰粱中央受绷p二“叨瞬一诊:却,的娜。:线(t)左固右铰梁中央受载P,10000时.知,的等0.线(c)左固右铰梁中央受该P,1‘000时塑性区分布图10(下转第9页)各类优化计算结果的比较粼、寸{‘生一址二二仁二‘一口火垫巴U一一}一一终o!29.35!,4。b(em)}于片~一一厂一一一万蔺-一一!一一一不犷万不一~~-一~~二一-三---!一{~—上丫__l曰“.ov一止红旦业全里工}_15600000}156而06丁esl不百丽而丽厂32.974022。27-‘-~闷比份‘悦‘,一一一一l一l一25。881519。111050000078000006.19X10一,1050000078000006.09X10一了1050000078000005.99X10一了156000001050000078000005.93X10一了304624371901545413486514311012721240683414858524刃5.o6xlo一,}4.71xl0一,一1一3一婆一1.83肠!一23,9呱55xzo5N/emb甲=60em,cl20”N/em,e:取Z,二l”二6oem,b‘=O,~110火一o5N/em,e:‘=85x=2ez‘,eZ“=2cZ‘.BOXI二RUH.Okl’明:歇二说,梢下主主甘入“,c.刚刚性程序的运算过程见图3由计算结果看出,.1口,泛月阮三巴,七月民,卜」上到匕飞.」活书.怕,下.C一匀.C.,1入毛刀叻C川).谁计必旅二的蓦彝度习9.匀呱,通过优化计算主轴前端(C川)一么曰劣万忿7:二二理,,。:‘一。a。。。下降了99.9肠,即静刚度提高了可以认为实现了无穷大刚度设计.参考文献杨叔子:三支承主轴部件静刚度的分析与讨.0“曰尹口洲).es…口a一艺C月个4,一曰91口反活(1〕论,机床,19了9年第9期〔2)金属切削机床设计编写组:金属切削机床设计(修