有限元數(shù)值分析中的剪力鎖閉(shearlock)

     在進行有限元數(shù)值分析時有時候會遇到計算結果與理論分析不一致,或者計算結果不合理的情況,雖然也選擇了合理的單元類型和材料本構模型,但得到的結果卻就是無法解釋,這時候可能就需要考慮所選用的單元插值函數(shù)是否合理了。這里來解釋一下有限元數(shù)值分析中剪力鎖閉出現(xiàn)的原理及處理方法。

一、剪力鎖閉的概念
     以受彎曲變形的梁來做為分析對象,在材料力學中有這樣的假定,梁在受彎矩作用下發(fā)生純彎曲變形時,變形前后平行于中截面的面仍然保持相互平行,且在變形前與中截面垂直的面在變形后仍然與中截面保持垂直。從連續(xù)梁中取一微小微元體分析,變形前后的圖形示意如圖1所示

圖1 連續(xù)梁微元體域受彎曲材料的變形

圖2 有限元分析得到的受彎曲材料的變形

    當用有限元對連續(xù)介質進行分析時,需要把連續(xù)介質離散為很多小單元和結點,當選取單元型式為線性單元時,每個邊由兩個結點構成,單元在變形前后邊仍然保持為直線。假如以四結點四邊形單元為例,受純彎荷載之后變形情況如圖2所示。從圖2可以看出(圖中虛線可以理解為通過高斯積分點的線),變形之后水平線仍然保持平行,水平線上面受拉應力而拉伸,下面受壓應力而壓縮。由于變形之后邊界仍然要保持直線,豎直線與水平線不再保持垂直,而是有了增大或減小,這說明產生了剪切變形,而這個剪應變是不應該出現(xiàn)的,由于這個剪切變形要消耗一定的變形能,這就會導致梁不能再發(fā)生彎曲或者彎曲變形要小了,產生的彎曲繞度也要減小。這種現(xiàn)象就稱為剪力鎖閉現(xiàn)象,通常出現(xiàn)在用常規(guī)單元模擬以彎曲變形為主要變形的結構中。這也是為什么我們用平面四結點四邊形單元模擬受純彎或者以彎曲變形為主的構件時計算得到彎曲變形較理論值要小的原因。

二、剪力鎖閉現(xiàn)象有限元驗證

    上面對剪力鎖閉現(xiàn)象做了個解釋,那有限元分析中是否真的存在如上分析的剪力鎖閉呢。來考查一個長1m,截面尺寸為0.1m*0.1m,兩端簡支的梁,中間受一個集中荷載(F=100kPa)的算例,結構示意圖如圖3所示,這樣的例子材料力學是有理論解的(y=pl^3/(48EI)=2.5cm)。(當然這個問題是可以用梁單元直接模擬的,這里為了說明剪力鎖閉現(xiàn)象,采用平面實體單元來分析)

材料參數(shù):E=1GPa,mu=0.2

圖3 簡支梁受集中荷載結構示意圖

圖4 梁結構有限元網格圖

   對上述不同有限元網格模型進行分析,采用PLANE182單元默認的積分方式,得到梁中間下面結點的豎直向位移如表1所示,表中ndivX為水平向剖分單元份數(shù),ndivY為梁厚度方向剖分的單元份數(shù)。從表中可以看出,只有當梁采用非常密的網格時才能得到與理論解比較接近的值。這說明:1、對這種細長結構,采用實體單元進行離散是非常不經濟的,而采用線性形狀的梁單元只需要很少的單元就能夠得到精確的結果;2、當單元比較少時,由于出現(xiàn)的剪力鎖閉現(xiàn)象導致計算結果偏小,這是由于剪力鎖閉消耗了一定的剪切能量,彎曲變形就小了。

表1  全積分計算結果(單位:cm)

    通過上面的分析,我們可以看出在用有限元進行分析時,剪力鎖閉現(xiàn)象的確是存在的。如果我們不去很好的分析結構受力情況,而只是按照教程一步一步的做,然后把這個過程套用到其它研究分析中,可能也計算得到了結果,但很有可能結果就是錯誤的。因此,一個問題就來了,那剪力鎖閉現(xiàn)象究竟在什么時候才會存在?通常情況下,所分析問題中是否會有剪力鎖閉現(xiàn)象,大概與兩個方面有關:一是單元形狀有關,當離散后的單元長比與短邊之比越大,單元越狹長,就越可能出現(xiàn)剪力鎖閉;二是與單元受力有關,在前面一個條件滿足時,當單元再承受垂直與長邊的荷載(如上面的梁承受法向荷載),使單元的變形以彎曲變形為主,這時出現(xiàn)剪力鎖閉的可能性就大。如上面的分析中,水平方向分10份時,由于單元都比較狹長,因此,由于剪力鎖閉的存在導致結果都失真;當水平方向分50份時,單元形狀接近正方形,不再是狹長單元,計算結果就接近真值。

三、消除剪力鎖閉現(xiàn)象的方法

    那么,當單元形狀和所受荷載不可避免剪力鎖閉時,有沒有好的解決辦法?當然是有的,現(xiàn)在很多商業(yè)軟件如ANSYS、ADINA、ABAQUS等都提供了相應的解決方法,大體上有兩種,一是采用積分減縮,二是采用強化應變單元模式。這里以ANSYS來介紹一下幾種方法。

3.1采用減縮積分來消除剪力鎖閉

     用常規(guī)的線性實體單元,軟件默認采用全積分格式,即在每一個方向存在兩個高斯積分點,這樣兩點確定一條直線,在單元變形后每條邊仍然是直線,就可能造成剪力鎖閉。如果對同一個單元只采用一個積分點,是不是就可以避免剪力鎖閉?這或許是可能的,但一個高斯積分點會不會降低計算精度呢?通過計算分析來說明,仍然采用上面的分析模型,同樣的計算網格。單元類型仍然取PLANE182號平面四結點等參單元,只是單元特性中將K2選項參數(shù)做些修改,改為Reduced Integration即減縮積分,如圖5(在其它商業(yè)軟件里也應該有類似的設置)。

    通過分析的結果如表2中減縮積分一列數(shù)值。從表2中可以看到,減縮積分模式下ANSYS的計算結果與理論解判別還是挺大,這個僅是針對本算例的結果,ANSYS幫助手冊及其它理論文獻中也都講減縮積分對消除剪力鎖閉是有幫助的。ABAQUS入門手冊中也給出了相應的分析數(shù)值,減縮積分后與理論解非常接近。

表2 消除剪力鎖閉的計算結果(單位:cm)

     這里值得一提的是當梁厚度方向只有一層單元時,減縮積分的結果居然出現(xiàn)63.36cm的值,完全失真。這是為什么?Why?這是因為減縮積分會導致另外一個奇怪的現(xiàn)象,就是傳說中的沙漏現(xiàn)象(hourglassing)。再次考慮受彎單元減縮積分時的變形情況,如圖6所示。單元由于采用了減縮積分模式,只有一個高斯積分點,通過單元高斯積分點的水平和垂直兩條線代表了單元變形后的情況。在單元彎曲變形之后,由于兩條線長度沒變,相互之間的夾角也沒有變,這就是說即沒有產生正應變,也沒有剪應變,因此,應變能為0,單元的彎曲變形成為一個零能量模式。這種單元模式下沒有剛度,不能夠抵抗變形。當網格比較粗時,如上面的只有一層單元,就會產生一個無意義的結果(從理論上講,由于沙漏現(xiàn)象的出現(xiàn),完全沒有剛度是無法計算的,但商業(yè)軟件都采用了一些技術手段保證能有一個計算結果)。

     當然,商業(yè)軟件在處理這種沒有剛度的減縮積分單元時會采取一些技術手段,如引入少量的人工“沙漏剛度”,讓其限制沙漏模式的擴展。但是,即便如此,如果單元數(shù)量不夠多的話,這種限制沙漏模式擴展的作用也不大,仍然會出現(xiàn)不理想的結果。只有當單元數(shù)量達到一定量后,才會起到明顯的作用。

圖6 受彎單元減縮積分時的變形

      雖然減縮積分可以在一定程度上消除剪力鎖閉現(xiàn)象,但有時候效果不佳,有時候還會出現(xiàn)沙漏現(xiàn)象。這里再來介紹另外一種方法,就是強化應變單元模式(enhanced strain),也稱為非協(xié)調單元。常規(guī)的線性單元位移沿一個方向是線性變化的,應變沿一個方向是常量,這就導致當單元受彎曲變形時容易出現(xiàn)剪力鎖閉。那么,如果采取某種方式讓應變沿一個呈線性變化是不是就可以解決剪力鎖閉?是的,的確如此(二次單元剪力鎖閉就比較難出現(xiàn),采用二次單元也是消除剪力鎖閉的方法之一)?；谶@種思路,就有人提出了強化應變的概念,在單元位移模式上做文章,通過增加一些虛擬的附加自由度,讓單元內部應變模式為線性變化,如圖7所示。由于這種增加變形梯度完全是在單元內部,與單元節(jié)點無關,因此,即不增加求解結構的整體自由度數(shù),也可以保證在邊界上位移仍然是連續(xù)的。

圖7   非協(xié)調單元模式(左)和常規(guī)單元模式(右)

     在ANSYS中平面的PLANE182、三維的SOLID185等單元都是支持強化應變模式的,是通過修改單元的選項值來實現(xiàn)的。對本算例,采用PLANE182單元,通過修改其K3選項值為2來實現(xiàn),如圖5所示。計算結果如表2最后一列的強化應變單元。從表2中的結果可以看出,采用強化變形模式后有限元數(shù)值分析的結果與理論解非常接近,這說明強化應變單元的效果非常不錯。

    當然,強化應變單元也有它本身的限制和弱點,比如當單元形狀比較畸形時計算結果會非常差,這里所說的畸形是指單元每條相鄰兩條線的夾角如果太大或太小,就可能導致非常不好的結果;若相鄰兩條邊的夾角接近90度,則是采用強化應變單元,那是再好不過的。

3.3采用高次單元消除剪力鎖閉現(xiàn)象

   通過前面的分析,知道了剪力鎖閉產生的本質原因是由于以彎曲變形為主要變形的結構在用有限元分析時,由于采用線性單元離散,在邊界上仍然是線性變形。那么增加單元的形函數(shù)階次,就可以在一定程度上消除剪力鎖閉現(xiàn)象。這當然是可以的,但問題就是由于單元節(jié)點增加,計算工作量就提高不少。

四、結語

     剪力鎖閉是否會產生,與結構受力、單元形狀、單元模式選取等多因素相關,當剪力鎖閉不可避免時,可以采用減縮積分、強化應變單元、高次單元等方法去盡可能的消除這種現(xiàn)象,但也要了解每一種方法的適用范圍、優(yōu)缺點等,選取合適、正確的方法去進行分析,只有這樣,才能得到合理的計算結果。

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

finish
/prep7

len=1.0
b=0.1
h=0.1
ndivX=10  !水平向份數(shù)
ndivY=8   !豎直向份數(shù)
force=-100e3

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

ET, 1, plane182,0,,0      ! 182號單元平面應力,182后面那個參數(shù)為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