沖擊力學(xué)有限元計(jì)算中的數(shù)值積分

1、第39卷第3期Vol. 39,No. 3 2009年3月J OURNAL OF UNIVE RSITY OF SCIE NCE AND TECHNOLO GY OF CHINAMar. 2009文章編號:025322778(2009 0320331206收稿日期:2006211213; 修回日期:2007206228基金項(xiàng)目:國家自然科學(xué)基金(10632080,10572134 資助.作者簡介:王峰, 男,1976年生, 博士. 研究方向:計(jì)算力學(xué). E 2mail : 通訊作者:王肖鈞, 教授. E 2mail :. cn沖

2、擊力學(xué)有限元計(jì)算中的數(shù)值積分王峰1,2, 王肖鈞1, 卞梁1, 1, 勞1(1. 中國科學(xué)技術(shù)大學(xué)近代力學(xué)系, 安徽合肥230027;2. , 摘要:在Lagrange 有限元基礎(chǔ)上, . 單點(diǎn)積分方法具有較高的計(jì)算效率, , ; 采用2×2×2高斯積分可以避免, ; 而采用局部2×2×2高斯積分則同時(shí)具有兩2×2×2高斯積分能夠很好地控制沙漏變形, 并有較高的計(jì)算效率; ,2×2×2高斯積分比單點(diǎn)積分更加接近理論值. 這說明所述方法和所建程序的合理性和有效性, 它為侵徹貫穿過程的數(shù)值分析提供了一種實(shí)用和有效的手段

3、.關(guān)鍵詞:爆炸力學(xué); 數(shù)值模擬; 有限元; 單點(diǎn)積分; 高斯積分中圖分類號:O385文獻(xiàn)標(biāo)識碼:ANumerical integration method of f inite elementcomputation in impact mechanicsWAN G Feng 1,2, WAN G Xiao 2jun 1, B IAN Liang 1, L I Jian 2guang 1, L AO J un 1(1. Depart ment of Modern Mechanics , Universit y of S cience and Technology of China , Hef e

4、i 230027, China;2. Depart ment of Unmanned A erial V ehicle and Mechanical Engineering , PL A A rtillery A cadem y , Hef ei 230031, China Abstract :Based on t he analysis of Lagrange finite element met hod , numerical integration met hods to comp ute nodal force were briefly described. In order to c

5、ont rol hourglass deformations , anti 2hourglass nodal forces have to be used for single point quadrat ure , which has higher efficiency. The hourglass modes would be cont rolled and t he numerical accuracy raised effectively if 2×2×2Gaussian quadrat ure were adapted , but t his would enta

6、il heavy comp utation. However , partial 2×2×2Gaussian quadrat ure has t he merit s of two kinds of integration met hod. The numerical example of penet ration indicated t hat partial 2×2×2Gaussian quadrat ure could cont rol hourglass deformations effectively and has higher effici

7、ency. And t he numerical simulations of one 2dimensional strain wave also shown t hat 2×2×2Gaussian quadrat ure is much better t han single point quadrat ure. It is concluded t hat t he met hod discussed and t he program we developed are reasonable and effective , providing a usef ul met h

8、od for t he numerical st udy of penet ration and perforation.K ey w ords :mechanics of explosion ; numerical simulation ; finite element ; single point quadrat ure ;Gaussian quadrat ure0引言侵徹問題的理論與實(shí)驗(yàn)研究始終是新武器研制和高科技領(lǐng)域中的重要課題, 無論是進(jìn)攻或者防御, 可靠預(yù)測侵徹效果、提高侵徹能力始終是新武器研究中的核心. 傳統(tǒng)的解析方法以及基于實(shí)驗(yàn)研究的經(jīng)驗(yàn)或半經(jīng)驗(yàn)方法已經(jīng)遠(yuǎn)遠(yuǎn)不能滿足現(xiàn)代高技術(shù)武

9、器發(fā)展的需要. 數(shù)值模擬方法, 由于具有假設(shè)少、物理模型逼真、計(jì)算精度高等明顯的優(yōu)勢,. ,程, . 為了提高數(shù)值模擬侵徹過程的計(jì)算效率, 通常采用單點(diǎn)積分法, 但這種方法的直接后果是立方元被認(rèn)為是常應(yīng)力常應(yīng)變單元, 在某些情況下會導(dǎo)致網(wǎng)格相交, 甚至網(wǎng)格翻轉(zhuǎn), 計(jì)算中斷, 因此必須通過引入抗沙漏節(jié)點(diǎn)力的方法以控制沙漏變形. 采用多點(diǎn)積分, 如2×2×2高斯積分不但可以提高計(jì)算精度, 還可以有效避免沙漏變形. 然而采用2×2×2高斯積分所引起計(jì)算量的增大是侵徹貫穿數(shù)值計(jì)算中必須面對的現(xiàn)實(shí)問題. 大量沖擊力學(xué)數(shù)值計(jì)算的實(shí)例表明, 出現(xiàn)網(wǎng)格嚴(yán)重畸變, 導(dǎo)致沙

10、漏變形發(fā)展失控的單元往往只是部分單元, 甚至是個(gè)別單元, 因此我們可以采用局部2×2×2高斯積分, 將單點(diǎn)積分計(jì)算效率高和全高斯積分可以消除沙漏變形的優(yōu)點(diǎn)結(jié)合起來. 在數(shù)值模擬侵徹貫穿問題中, 對部分網(wǎng)格嚴(yán)重畸變單元采用2×2×2高斯積分, 而其他的絕大部分單元仍采用單點(diǎn)積分, 計(jì)算實(shí)踐表明, 這一方法是有效的.本文對八節(jié)點(diǎn)立方元分別采用不同積分方法對平頭彈侵徹靶板進(jìn)行了數(shù)值模擬計(jì)算, 表明局部2×2×2高斯積分能夠很好地控制沙漏變形, 精度接近2×2×2高斯積分, 且具有較高的計(jì)算效率. 還采用單點(diǎn)積分和2

11、15;2×2高斯積分分別計(jì)算了板中一維應(yīng)變波的傳播, 并與理論解進(jìn)行了比較, 結(jié)果表明2×2×2高斯積分的計(jì)算結(jié)果明顯優(yōu)于單點(diǎn)積分.1計(jì)算方法簡介1. 1單點(diǎn)積分的沙漏變形與沙漏控制目前, 在大多數(shù)沖擊力學(xué)有限元程序里, 為提高計(jì)算效率, 在計(jì)算節(jié)點(diǎn)力時(shí), 通常采用的是單點(diǎn)積分法, 對于八節(jié)點(diǎn)六面體立方元(圖1 而言, 即在自然 坐標(biāo)系里, 在進(jìn)行單元體積分時(shí), 積分點(diǎn)取為=0. 單點(diǎn)積分的直接后果是立方元被認(rèn)為是常應(yīng)力常應(yīng)變單元, 因此其對應(yīng)的位移場只能是線性位移場, 顯然, 單點(diǎn)積分無法描述非線性位移模態(tài). 當(dāng)單元變形出現(xiàn)如圖2時(shí), , 這就是所謂零. , f

12、 int i =展, , 甚至網(wǎng)格翻 .針對單點(diǎn)積分中出現(xiàn)的這一問題,Belyt schko 等3,4提出了一種通過引入抗沙漏節(jié)點(diǎn)力的方法控制沙漏變形. 其基本思想是根據(jù)沙漏變形模態(tài)引入相應(yīng)的沙漏變形率,q i =v iI I . (1式中, 角標(biāo)i 表示空間坐標(biāo)軸方向, I 表示節(jié)點(diǎn), 表示沙漏變形模態(tài); v iI 為節(jié)點(diǎn)速度分量, I 為沙漏變形分量, 定義為I =I -V B iI x iJJ . (2 式中, B iJ =5N I (x i /5x i , I 如表1所示. 由虛功原理出發(fā), 可以建立抗沙漏節(jié)點(diǎn)力,f H G iI =Q i I . (3233中國科學(xué)技術(shù)大學(xué)學(xué)報(bào)第39

13、卷式中, Q i 是與沙漏變形相對應(yīng)的沙漏應(yīng)力,Q i =k3V q i , (4k 為人工沙漏控制系數(shù), 和是拉梅常數(shù).1. 22×2×2高斯積分的計(jì)算過程單點(diǎn)積分所對應(yīng)的位移場是線性位移場, 實(shí)際上八節(jié)點(diǎn)立方元所設(shè)定的位移函數(shù)應(yīng)包含若干非線性項(xiàng), 或者說在自然坐標(biāo)系里, 形函數(shù)應(yīng)具有如下的形式,N I (, , =8(1+ +I 8+I 3I +1I +2I +3I +4I .(5 式中, I , I , I , I , 1I , 2I , 3I , 1I , 2I , 3I , 4I 分別表示與節(jié)點(diǎn)相對應(yīng)的自然坐標(biāo)取值, 如表1所示, 其中, 1I , 2I , 3I

14、 , 4I 刻畫了位移場的非線性模態(tài). 表1節(jié)點(diǎn)的自然坐標(biāo)取值T ab. 1N atural coordinates of nodesv ( x , t =N ( x u (t . (6 式中, N ( x 是由形函數(shù)構(gòu)成的3×24矩陣,N ( x = N 100N 800 0N 100N 80 00N 100N ;u (t 表示由節(jié)點(diǎn)速度構(gòu)成的24×1列向量, u (t =u 11u 21u 31u 12u 22u 32u 18u 28u 38T ,u ij 表示第j 個(gè)節(jié)點(diǎn)i 方向的速度分量. u (t 只是時(shí)間t 的函數(shù). 于是,v ( x , t =N ( x u

15、(t . (7 應(yīng)變率矢量可以表示成11223312=x 1000x200x321x 3x 2x 30x 1v 1(x , t 2(x , t 3(x , t =L v (x , t . (8 式中,ii 為法向應(yīng)變率,ii =5xi;ij 為剪應(yīng)變率,ij =xj+xi; L 是由微分算子構(gòu)成的矩陣, v i (x , t 表示x i 方向上的速度分量. 將式(6 代入式(8 , 有= L N ( x u (t = B ( x u (t . (9 式中, B 稱為梯度矩陣, 由形函數(shù)的空間導(dǎo)數(shù)構(gòu)成, B ( x = L N ( x .由式(5 可知, 形函數(shù)N i (, , 是由自然坐標(biāo)給出

16、的, 根據(jù)偏微分法則,5=5x 15+5x 25+5x 35. (10 對于其余兩坐標(biāo)(, , 可以類似得到, 寫成矩陣形式有5=555x 1x 2x 3=Jx 1 x 2 x 3 (11 式中, J 稱為雅可比矩陣, 可記為(, , . 在等參元里, 空間坐標(biāo)可以表示為x j =N i x ji =N 1x j 1+N 2x j 2+N 8x j 8; i =1, 2, , 8; j =1, 2, 3.(12 式中, x ji 表示第i 個(gè)節(jié)點(diǎn)j 方向的坐標(biāo)值. 利用式(12 , J 可以顯式地表現(xiàn)為自然坐標(biāo)的函數(shù),J (, , =333第3期沖擊力學(xué)有限元計(jì)算中的數(shù)值積分 x 1i x 2

17、i x 3i x 1i x 2i x 3i 5x 1i 5x 2i 5 x 3i= x 11x 21x 31x 12x 22x 32x .(13于是N i 對于, , x 3 示為x 1x 2x 3=J -1.(14式中, J -1是J 的逆矩陣.體積微元可以表示為d V =d (d ×d . (15式中 ,d =d i +d j +d k , d =d i +d j +d k , d =d i +d j +d k ; i , j , k 是直角坐標(biāo)的單位向量, 于是d V =d d d =|J |d d d . (16節(jié)點(diǎn)內(nèi)力可以表示為f inti (t =eB T i ( x (

18、 x , t d =VB T i (, , (, , , t d V =|J |1-11-11-1B Ti(, , (, , , t d d d .(17采用2×2×2高斯積分5后, 式(17 可寫為 f int i (t =|J |2j =12k =12m =1HjH k H m B T i (j , k , m (j , k , m , t .(18式中,Ti 10x 20x 3x 20x 1x 30x 30x 2x 11122331223=x 111+x 212+x 331x 222+x 112+x 323x 333+x 223+x 131. (19對于2×

19、2×2高斯積分, H j =H k =H m =1, j =k =m =±0157735(j , k , m =1, 2 , 將式(14 和式(18 代入式(17 可以計(jì)算得到節(jié)點(diǎn)內(nèi)力. 由于在每個(gè)坐標(biāo)方向的積分點(diǎn)數(shù)取為n =2, 數(shù)學(xué)上可以證明, 2×2×2高斯積分的結(jié)果可以達(dá)到2×n -1=3階的精度. 也就是說, 如果 B T i (, , (, , =i j k, 且i , j , k 3, 則式(18 將能給出積分的精確解.2數(shù)值模擬計(jì)算及結(jié)果分析2. 1平頭桿彈對靶板的三維侵徹計(jì)算為了考查沙漏變形對有限元計(jì)算過程的影響,我們對平頭桿

20、彈侵徹靶板進(jìn)行了數(shù)值計(jì)算. 彈體和靶板材料取為線性硬化體, 其靜水壓力與體積變形之間用如下的Gruneisen 狀態(tài)方程描述6:P =(K 1+K 22+K 33(1-0. 5 +E. (20式中, =V 0/V -1為壓縮比, 為Gruneisen 參數(shù),E 為比內(nèi)能.應(yīng)力偏量服從增量型彈塑性應(yīng)力應(yīng)變關(guān)系, 即d S ij =2G (d e ij -d e p ij. (21 材料服從Mises 屈服準(zhǔn)則, 屈服強(qiáng)度取為線性硬化函數(shù),Y =Y 0+(Y u -Y 0 P / u .(22式中, Y 0為初始屈服應(yīng)力, Y u 為極限屈服應(yīng)力, P433中國科學(xué)技術(shù)大學(xué)學(xué)報(bào)第39卷表2材料參數(shù)

21、T ab. 2P arameters of m aterial/(kg m -3G /GPaY 0/GPaY u /GPaP /GPa u狀態(tài)方程參數(shù)/GPaK 1K 2K 3782377. 51. 432. 501. 700. 011642945002. 圖3不同積分方法下的侵徹圖像Fig. 3Penetration conf igurations of different integration methods為等效塑性應(yīng)變, u 為極限應(yīng)變.計(jì)算中用到的有關(guān)物理參數(shù)6和計(jì)算參數(shù)如表2所示.圖3給出了平頭桿彈以初速度V 0=1500m/s 侵徹靶板t =2710s 時(shí)的應(yīng)變及網(wǎng)格變形圖像.

22、 可以看出, 加入抗沙漏節(jié)點(diǎn)力的單點(diǎn)積分和2×2×2高斯積分都能夠很好地控制沙漏變形, 避免網(wǎng)格畸變, 而取消抗沙漏節(jié)點(diǎn)力的單點(diǎn)積分則出現(xiàn)了網(wǎng)格的畸變和相交, 嚴(yán)重影響了有限元計(jì)算的正常計(jì)算過程, 這也進(jìn)一步說明了上述計(jì)算方法的正確性.盡管2×2×2高斯積分的計(jì)算結(jié)果可以更為準(zhǔn)確地給出桿彈靶板的物理圖像, 然而其計(jì)算量的增加依然是數(shù)值模擬高速侵徹和貫穿問題時(shí)不可忽視的因素. 事實(shí)上, 計(jì)算中我們發(fā)現(xiàn), 出現(xiàn)網(wǎng)格畸變嚴(yán)重導(dǎo)致沙漏變形模態(tài)發(fā)展的只是部分單元, 這部分單元幾乎全部集中在侵徹通道附近. 為此我們開展了一組只對侵徹通道附近靶板單元采用2×2×2高斯積分, 其余單元都采取不加抗