有限元数值分析中的剪力锁闭(shearlock)

2016-08-30 by:CAE仿真在线来源:互联网

在进行有限元数值分析时有时候会遇到计算结果与理论分析不一致,或者计算结果不合理的情况,虽然也选择了合理的单元类型和材料本构模型,但得到的结果却就是无法解释,这时候可能就需要考虑所选用的单元插值函数是否合理了。这里来解释一下有限元数值分析中剪力锁闭出现的原理及处理方法。

一、剪力锁闭的概念
以受弯曲变形的梁来做为分析对象,在材料力学中有这样的假定,梁在受弯矩作用下发生纯弯曲变形时,变形前后平行于中截面的面仍然保持相互平行,且在变形前与中截面垂直的面在变形后仍然与中截面保持垂直。从连续梁中取一微小微元体分析,变形前后的图形示意如图1所示

有限元数值分析中的剪力锁闭(shear <wbr>lock)

图1 连续梁微元体域受弯曲材料的变形

从图1中可以看出,在受到弯矩作用下,梁顶面受拉而伸长为水平拉应力,底面受压而缩短为水平压应力,在变形前后竖直方向各线段长度不变,即竖直方向就没有应变,也就没有应力;此外,水平线与竖直线在变形前后仍然垂直,也就没有剪切变形,剪应力为0,这解释了材料力学的假定。需要注意一个重要现象就受拉压的上下边界在变形后是弯曲的。有限元数值分析中的剪力锁闭(shear <wbr>lock)

图2 有限元分析得到的受弯曲材料的变形

当用有限元对连续介质进行分析时,需要把连续介质离散为很多小单元和结点,当选取单元型式为线性单元时,每个边由两个结点构成,单元在变形前后边仍然保持为直线。假如以四结点四边形单元为例,受纯弯荷载之后变形情况如图2所示。从图2可以看出(图中虚线可以理解为通过高斯积分点的线),变形之后水平线仍然保持平行,水平线上面受拉应力而拉伸,下面受压应力而压缩。由于变形之后边界仍然要保持直线,竖直线与水平线不再保持垂直,而是有了增大或减小,这说明产生了剪切变形,而这个剪应变是不应该出现的,由于这个剪切变形要消耗一定的变形能,这就会导致梁不能再发生弯曲或者弯曲变形要小了,产生的弯曲绕度也要减小。这种现象就称为剪力锁闭现象,通常出现在用常规单元模拟以弯曲变形为主要变形的结构中。这也是为什么我们用平面四结点四边形单元模拟受纯弯或者以弯曲变形为主的构件时计算得到弯曲变形较理论值要小的原因。

二、剪力锁闭现象有限元验证

上面对剪力锁闭现象做了个解释,那有限元分析中是否真的存在如上分析的剪力锁闭呢。来考查一个长1m,截面尺寸为0.1m*0.1m,两端简支的梁,中间受一个集中荷载(F=100kPa)的算例,结构示意图如图3所示,这样的例子材料力学是有理论解的(y=pl^3/(48EI)=2.5cm)。(当然这个问题是可以用梁单元直接模拟的,这里为了说明剪力锁闭现象,采用平面实体单元来分析)

材料参数:E=1GPa,mu=0.2

图3 简支梁受集中荷载结构示意图

我们用平面应力问题,采用四结点四边形等参单元来离散,采用不同的网格划分密度进行分析,在长度方向上剖分10份,在厚度方向上分别剖分1份、2份、4份和8份,其中对厚度剖分8份的情况还进行了水平向剖分20份和50份的分析,有限元网格图如图4所示,然后取梁中间最下面一个结点的的竖直向位移与材料力学的理论解进行对比。

图4 梁结构有限元网格图

对上述不同有限元网格模型进行分析,采用PLANE182单元默认的积分方式,得到梁中间下面结点的竖直向位移如表1所示,表中ndivX为水平向剖分单元份数,ndivY为梁厚度方向剖分的单元份数。从表中可以看出,只有当梁采用非常密的网格时才能得到与理论解比较接近的值。这说明:1、对这种细长结构,采用实体单元进行离散是非常不经济的,而采用线性形状的梁单元只需要很少的单元就能够得到精确的结果;2、当单元比较少时,由于出现的剪力锁闭现象导致计算结果偏小,这是由于剪力锁闭消耗了一定的剪切能量,弯曲变形就小了。

表1 全积分计算结果(单位:cm)

通过上面的分析,我们可以看出在用有限元进行分析时,剪力锁闭现象的确是存在的。如果我们不去很好的分析结构受力情况,而只是按照教程一步一步的做,然后把这个过程套用到其它研究分析中,可能也计算得到了结果,但很有可能结果就是错误的。因此,一个问题就来了,那剪力锁闭现象究竟在什么时候才会存在?通常情况下,所分析问题中是否会有剪力锁闭现象,大概与两个方面有关:一是单元形状有关,当离散后的单元长比与短边之比越大,单元越狭长,就越可能出现剪力锁闭;二是与单元受力有关,在前面一个条件满足时,当单元再承受垂直与长边的荷载(如上面的梁承受法向荷载),使单元的变形以弯曲变形为主,这时出现剪力锁闭的可能性就大。如上面的分析中,水平方向分10份时,由于单元都比较狭长,因此,由于剪力锁闭的存在导致结果都失真;当水平方向分50份时,单元形状接近正方形,不再是狭长单元,计算结果就接近真值。

三、消除剪力锁闭现象的方法

那么,当单元形状和所受荷载不可避免剪力锁闭时,有没有好的解决办法?当然是有的,现在很多商业软件如ANSYS、ADINA、ABAQUS等都提供了相应的解决方法,大体上有两种,一是采用积分减缩,二是采用强化应变单元模式。这里以ANSYS来介绍一下几种方法。

3.1采用减缩积分来消除剪力锁闭

用常规的线性实体单元,软件默认采用全积分格式,即在每一个方向存在两个高斯积分点,这样两点确定一条直线,在单元变形后每条边仍然是直线,就可能造成剪力锁闭。如果对同一个单元只采用一个积分点,是不是就可以避免剪力锁闭?这或许是可能的,但一个高斯积分点会不会降低计算精度呢?通过计算分析来说明,仍然采用上面的分析模型,同样的计算网格。单元类型仍然取PLANE182号平面四结点等参单元,只是单元特性中将K2选项参数做些修改,改为Reduced Integration即减缩积分,如图5(在其它商业软件里也应该有类似的设置)。

图5 ANSYS中设置单元减缩积分模式

通过分析的结果如表2中减缩积分一列数值。从表2中可以看到,减缩积分模式下ANSYS的计算结果与理论解判别还是挺大,这个仅是针对本算例的结果,ANSYS帮助手册及其它理论文献中也都讲减缩积分对消除剪力锁闭是有帮助的。ABAQUS入门手册中也给出了相应的分析数值,减缩积分后与理论解非常接近。

表2 消除剪力锁闭的计算结果(单位:cm)

这里值得一提的是当梁厚度方向只有一层单元时,减缩积分的结果居然出现63.36cm的值,完全失真。这是为什么?Why?这是因为减缩积分会导致另外一个奇怪的现象,就是传说中的沙漏现象(hourglassing)。再次考虑受弯单元减缩积分时的变形情况,如图6所示。单元由于采用了减缩积分模式,只有一个高斯积分点,通过单元高斯积分点的水平和垂直两条线代表了单元变形后的情况。在单元弯曲变形之后,由于两条线长度没变,相互之间的夹角也没有变,这就是说即没有产生正应变,也没有剪应变,因此,应变能为0,单元的弯曲变形成为一个零能量模式。这种单元模式下没有刚度,不能够抵抗变形。当网格比较粗时,如上面的只有一层单元,就会产生一个无意义的结果(从理论上讲,由于沙漏现象的出现,完全没有刚度是无法计算的,但商业软件都采用了一些技术手段保证能有一个计算结果)。

当然,商业软件在处理这种没有刚度的减缩积分单元时会采取一些技术手段,如引入少量的人工“沙漏刚度”,让其限制沙漏模式的扩展。但是,即便如此,如果单元数量不够多的话,这种限制沙漏模式扩展的作用也不大,仍然会出现不理想的结果。只有当单元数量达到一定量后,才会起到明显的作用。

图6 受弯单元减缩积分时的变形

3.2采用强化应变单元模式(非协调单元)来消除剪力锁闭

虽然减缩积分可以在一定程度上消除剪力锁闭现象,但有时候效果不佳,有时候还会出现沙漏现象。这里再来介绍另外一种方法,就是强化应变单元模式(enhanced strain),也称为非协调单元。常规的线性单元位移沿一个方向是线性变化的,应变沿一个方向是常量,这就导致当单元受弯曲变形时容易出现剪力锁闭。那么,如果采取某种方式让应变沿一个呈线性变化是不是就可以解决剪力锁闭?是的,的确如此(二次单元剪力锁闭就比较难出现,采用二次单元也是消除剪力锁闭的方法之一)。基于这种思路,就有人提出了强化应变的概念,在单元位移模式上做文章,通过增加一些虚拟的附加自由度,让单元内部应变模式为线性变化,如图7所示。由于这种增加变形梯度完全是在单元内部,与单元节点无关,因此,即不增加求解结构的整体自由度数,也可以保证在边界上位移仍然是连续的。

图7 非协调单元模式(左)和常规单元模式(右)

在ANSYS中平面的PLANE182、三维的SOLID185等单元都是支持强化应变模式的,是通过修改单元的选项值来实现的。对本算例,采用PLANE182单元,通过修改其K3选项值为2来实现,如图5所示。计算结果如表2最后一列的强化应变单元。从表2中的结果可以看出,采用强化变形模式后有限元数值分析的结果与理论解非常接近,这说明强化应变单元的效果非常不错。

当然,强化应变单元也有它本身的限制和弱点,比如当单元形状比较畸形时计算结果会非常差,这里所说的畸形是指单元每条相邻两条线的夹角如果太大或太小,就可能导致非常不好的结果;若相邻两条边的夹角接近90度,则是采用强化应变单元,那是再好不过的。

3.3采用高次单元消除剪力锁闭现象

通过前面的分析,知道了剪力锁闭产生的本质原因是由于以弯曲变形为主要变形的结构在用有限元分析时,由于采用线性单元离散,在边界上仍然是线性变形。那么增加单元的形函数阶次,就可以在一定程度上消除剪力锁闭现象。这当然是可以的,但问题就是由于单元节点增加,计算工作量就提高不少。

四、结语

剪力锁闭是否会产生,与结构受力、单元形状、单元模式选取等多因素相关,当剪力锁闭不可避免时,可以采用减缩积分、强化应变单元、高次单元等方法去尽可能的消除这种现象,但也要了解每一种方法的适用范围、优缺点等,选取合适、正确的方法去进行分析,只有这样,才能得到合理的计算结果。

==================附本分析中用到的ANSYS命令流========================

finish
/prep7

len=1.0
b=0.1
h=0.1
ndivX=10 !水平向份数
ndivY=8 !竖直向份数
force=-100e3

mp,dens,1,density !定义材料
mp,ex,1,1e9
mp,nuxy,1,0.2

ET, 1, plane182,0,,0 ! 182号单元平面应力,182后面那个参数为0:全积分;1:减缩积分;2:强化应变

rect,0,len,0,h

lsel,s,loc,x,0.1,len-0.1
lesize,all,,,ndivX
lsel,inve
lesize,all,,,ndivY
allsel,all

mshape,0,2d
mshkey,1
amesh,all

/pnum,mat,1
/number,1

save
finish

/solu
antype,static

allsel,all
nsel,s,loc,x,0
nsel,r,loc,y,0
d,all,all,0
allsel,all
nsel,s,loc,x,len
nsel,r,loc,y,0
d,all,uy,0
allsel,all
nsel,s,loc,x,len/2.
nsel,r,loc,y,h
f,all,fy,force

save
solve
finish

/post1
set,last
nsel,s,loc,x,len/2.
nsel,r,loc,y,0
prnsol,u,y