基于物質(zhì)點(diǎn)法金屬體積成形過(guò)程的仿真

謝桂蘭 于 超 龔曙光 田 杰 胡駿迪

湘潭大學(xué)，湘潭，411105

?

基于物質(zhì)點(diǎn)法金屬體積成形過(guò)程的仿真

謝桂蘭 于 超 龔曙光 田 杰 胡駿迪

湘潭大學(xué)，湘潭，411105

針對(duì)傳統(tǒng)有限元法開(kāi)展金屬體積成形分析過(guò)程中存在的網(wǎng)格畸變，以及在常規(guī)無(wú)網(wǎng)格法分析中存在的本質(zhì)邊界條件施加困難等問(wèn)題，建立了基于物質(zhì)點(diǎn)法金屬體積成形過(guò)程的仿真模型。因其背景網(wǎng)格積分采用了有限元形函數(shù)，從而有效地解決了邊界條件施加困難的問(wèn)題。通過(guò)對(duì)圓柱頂鐓和反向擠壓等金屬體積成形過(guò)程的仿真，并將仿真結(jié)果與實(shí)驗(yàn)和有限元結(jié)果進(jìn)行了對(duì)比分析，結(jié)果顯示物質(zhì)點(diǎn)法不僅能有效地消除網(wǎng)格畸變，而且在金屬發(fā)生大變形時(shí)其計(jì)算精度優(yōu)于有限元法的計(jì)算精度。所得結(jié)論可為物質(zhì)點(diǎn)法在彈塑性大變形分析中的應(yīng)用提供指導(dǎo)。

物質(zhì)點(diǎn)法；彈塑性大變形；金屬體積成形；數(shù)值模擬

0 引言

金屬體積成形是固態(tài)金屬材料成形的主要方式，其成形過(guò)程屬于幾何非線性和物理非線性的大變形問(wèn)題。目前有限元(finite element,FE)數(shù)值模擬技術(shù)是金屬體積成形問(wèn)題廣泛采用的分析方法，但由于有限元法依賴于網(wǎng)格，在開(kāi)展大變形分析時(shí)易產(chǎn)生網(wǎng)格畸變，從而造成計(jì)算失敗，故需要對(duì)其進(jìn)行網(wǎng)格重構(gòu)，同時(shí)新舊網(wǎng)格節(jié)點(diǎn)場(chǎng)量映射會(huì)降低其計(jì)算效率與精度[1-2]。針對(duì)有限元法的不足，近年來(lái)SPH(smooth particle hydrodynamics)[3]、RKPM(reproducing kernel particle method)[4]、EFG(element free Galerkin method)[5]等無(wú)網(wǎng)格法在計(jì)算金屬體積成形問(wèn)題中得到廣泛應(yīng)用，但它們存在本質(zhì)邊界條件施加困難[6]、粒子響域搜索費(fèi)時(shí)[7]、積分時(shí)間步長(zhǎng)隨粒子間距減小而減小[8]等不足。

Sulsky等[9]提出了一種新型的數(shù)值方法即物質(zhì)點(diǎn)法(material point method, MPM)。它將材料域離散成質(zhì)點(diǎn)集合，質(zhì)點(diǎn)攜帶所有的物質(zhì)信息如質(zhì)量、速度、應(yīng)力、應(yīng)變等，且用質(zhì)點(diǎn)的運(yùn)動(dòng)來(lái)描述材料的流動(dòng)與變形，同時(shí)采用規(guī)則背景網(wǎng)格來(lái)計(jì)算空間導(dǎo)數(shù)和求解動(dòng)量方程[10]。目前物質(zhì)點(diǎn)法已在高速碰撞[11]、爆炸[12]、動(dòng)態(tài)裂紋擴(kuò)展[13]等領(lǐng)域得到應(yīng)用。Wieckowski[14]借助MPM研究了顆粒材料流動(dòng)、金屬材料切削等問(wèn)題；Carter等[15]采用MPM模擬了沙粒崩潰現(xiàn)象；Nair等[16]針對(duì)廣義物質(zhì)點(diǎn)法提出了一種隱式時(shí)間積分方法，并應(yīng)用于超彈性問(wèn)題的求解。鑒于MPM具有拉格朗日和歐拉算法的優(yōu)勢(shì)，非常適合模擬涉及材料特大變形的斷裂破碎等問(wèn)題，本文開(kāi)展了基于MPM的金屬體積成形的仿真研究。

1 物質(zhì)點(diǎn)法(MPM)簡(jiǎn)介

MPM的基本思想是將連續(xù)體離散成若干個(gè)集中質(zhì)量的物質(zhì)點(diǎn)集合，并按照連續(xù)體的形狀分布在背景網(wǎng)格內(nèi)，在整個(gè)計(jì)算過(guò)程中，物質(zhì)點(diǎn)和背景網(wǎng)格節(jié)點(diǎn)之間應(yīng)用節(jié)點(diǎn)形函數(shù)完成兩次映射計(jì)算并采用顯式時(shí)間積分算法求解，以實(shí)現(xiàn)各物質(zhì)點(diǎn)應(yīng)力和應(yīng)變的更新。圖1為物質(zhì)點(diǎn)法離散示意圖。其理論推導(dǎo)見(jiàn)文獻(xiàn)[9-10]。

圖1 物質(zhì)點(diǎn)法示意圖

采用MPM求解時(shí)，連續(xù)體滿足下面的控制方程(不考慮熱量交換)。

動(dòng)量方程：

ρa(bǔ)=·σ+ρb

(1)

質(zhì)量守恒方程：

(2)

式中，ρ為密度；a為加速度矢量；v為速度矢量；σ為柯西應(yīng)力張量；b為單位質(zhì)量上的體積力。

任取試函數(shù)w代入式(1)得動(dòng)量方程的弱形式為

(3)

式中，Ω為當(dāng)前構(gòu)形；σs為比應(yīng)力；ts為比邊界面力。

當(dāng)材料域用質(zhì)點(diǎn)離散時(shí)，連續(xù)體的密度等于所有物質(zhì)點(diǎn)質(zhì)量之和，寫(xiě)成δ函數(shù)為

(4)

式中，Np為質(zhì)點(diǎn)總數(shù)；mj為質(zhì)點(diǎn)j的質(zhì)量；δ為Dirac Delta函數(shù)；xij為t時(shí)刻質(zhì)點(diǎn)j的坐標(biāo)。

將式(4)代入式(3)，積分形式轉(zhuǎn)變?yōu)榍蠛托问剑?/p>

(5)

式中，h為假想的邊界層厚度。

物質(zhì)點(diǎn)法在背景網(wǎng)格上計(jì)算動(dòng)量方程，通過(guò)背景網(wǎng)格節(jié)點(diǎn)的形函數(shù)NI(x)來(lái)實(shí)現(xiàn)質(zhì)點(diǎn)與背景網(wǎng)格之間的相互映射。對(duì)于三維問(wèn)題背景網(wǎng)格采用規(guī)則八節(jié)點(diǎn)六面體，其形函數(shù)為

(6)

I=1,2,…,8

其中,ξI、ηI、ζI分別是網(wǎng)格節(jié)點(diǎn)I的自然坐標(biāo)。由于形函數(shù)NI(x)兼具緊支性和插值性，故質(zhì)點(diǎn)xp的影響域是其所在網(wǎng)格單元，同時(shí)也消除了無(wú)網(wǎng)格法中本質(zhì)邊界條件施加困難的缺陷。

質(zhì)點(diǎn)的任意某參數(shù)cp可通過(guò)背景網(wǎng)格節(jié)點(diǎn)相應(yīng)的參數(shù)cI插值得到：

(7)

(8)

其中，Ng是背景網(wǎng)格單元節(jié)點(diǎn)總數(shù)，NIp=NI(xp)表示網(wǎng)格節(jié)點(diǎn)I的形函數(shù)在質(zhì)點(diǎn)p處的值，c參數(shù)可以代表位移、加速度、速度等。

將式(5)中質(zhì)點(diǎn)位移、加速度以及試函數(shù)由背景網(wǎng)格節(jié)點(diǎn)的插值替代，整理得

(9)

J=1,2,…,Ng

(10)

(11)

(12)

為了加快求解速度，采用集中質(zhì)量矩陣元素：

(13)

則動(dòng)量方程式(9)簡(jiǎn)化為

(14)

I=1,2,…,Ng

MPM通過(guò)在背景網(wǎng)格上求解動(dòng)量方程即式(14)，獲得網(wǎng)格節(jié)點(diǎn)運(yùn)動(dòng)信息并返回質(zhì)點(diǎn)，用來(lái)更新質(zhì)點(diǎn)的運(yùn)動(dòng)狀態(tài)以及應(yīng)力應(yīng)變等。因?yàn)楸尘熬W(wǎng)格每次參加計(jì)算后都可重新生成，所以變形不會(huì)累加，網(wǎng)格不會(huì)發(fā)生畸變。

2 金屬成形的MPM仿真模型

(15)

σijp=Sijp+σmpδij

(16)

采用J2流動(dòng)理論的徑向返回算法[17]對(duì)應(yīng)力進(jìn)行更新,返回映射法由彈性預(yù)測(cè)步和塑性修正步組成，在彈性預(yù)測(cè)步中有

(17)

(18)

將式(17)和式(18)代入式(15)得到應(yīng)力更新格式：

(19)

(20)

(21)

(22)

(23)

(24)

(25)

根據(jù)物質(zhì)點(diǎn)法的基本理論及金屬體積成形的仿真模型，利用FORTRAN語(yǔ)言編寫(xiě)了金屬體積成形MPM程序。下面采用編寫(xiě)的程序分別對(duì)棒材頂鐓和反向擠壓兩種工藝過(guò)程進(jìn)行仿真模擬。

3 棒材頂鐓工藝過(guò)程模擬

圖2為棒材頂鐓成形過(guò)程示意圖。選用文獻(xiàn)[18]中的幾何參數(shù)和材料模型，已知圓棒坯料初始尺寸為φ20 mm×49.69 mm，凹模型孔尺寸為φ20 mm×20 mm。坯料材料為AlMgF10，材料在室溫下的應(yīng)力應(yīng)變關(guān)系為σ=180.65×106ε0.183(Pa)。

圖2 圓棒頂鐓模型

凹模固定不動(dòng)，凸模向下運(yùn)動(dòng)，下降速度v=5 m/s。凸模和凹模簡(jiǎn)化為剛體，模具與工件的摩擦采用庫(kù)侖摩擦，摩擦因數(shù)取0.12。根據(jù)模型的對(duì)稱特點(diǎn)，本文采用1/4模型，棒材被均勻離散成11 067個(gè)質(zhì)點(diǎn)，如圖3a所示。通過(guò)仿真分析得到自由段高度方向的變形比(D=Δh/h)分別為40%、50%、60%的結(jié)果如圖3b～3d所示。圖4所示為文獻(xiàn)[18]的實(shí)驗(yàn)結(jié)果。

(a)D=0 (b)D=40% (c)D=50% (d)D=60%圖3 圓棒在不同變形程度下的結(jié)果

圖4 文獻(xiàn)[18]實(shí)驗(yàn)結(jié)果

對(duì)比圖3和圖4可知，本文MPM計(jì)算所得的工件變形與實(shí)驗(yàn)結(jié)果相吻合，質(zhì)點(diǎn)的排布也符合金屬流動(dòng)規(guī)律，圓棒下端由于凹模的固定作用而保持不變，上端材料沿徑向流動(dòng)，從而形成了頂鐓件的頭部。

為進(jìn)一步研究工件與模具在a和b接觸處的金屬流動(dòng)規(guī)律，設(shè)RS=Ra/R0，RI=Rb/R0，其中，Ra、Rb分別為圓棒在位置a和b處的半徑，R0為圓棒初始半徑。RS、RI隨下壓行程s的變化曲線如圖5所示。

圖5 RS和RI隨下壓行程的變化曲線

從圖5可看到，MPM得到的RS、RI曲線與實(shí)驗(yàn)結(jié)果相吻合，而FEM在下壓量小時(shí)結(jié)果也相吻合。但隨著下壓量增大，利用FEM得到的RS、RI曲線與實(shí)驗(yàn)結(jié)果產(chǎn)生了偏離，這主要是因?yàn)殡S著下壓量增大，坯料的變形程度增加，此時(shí)FEM網(wǎng)格發(fā)會(huì)生畸變，需要不斷地進(jìn)行網(wǎng)格重構(gòu)。

圓棒在頂鐓時(shí)頭部變形呈鼓形。壓下量不同時(shí)MPM與FEM計(jì)算結(jié)果的最大鼓形半徑對(duì)比見(jiàn)表1。從表1中得知，MPM結(jié)果與FEM結(jié)果吻合較好，最大誤差為0.62%。

表1 最大鼓形半徑

(a)MPM (b)FEM圖6 等效應(yīng)變分布云圖

(a)MPM (b)FEM圖7 等效應(yīng)力分布云圖

D=60%時(shí)MPM和FEM計(jì)算得到的等效應(yīng)變和應(yīng)力的云圖分別見(jiàn)圖6和圖7。比較圖6與圖7可知，MPM計(jì)算所得的應(yīng)力應(yīng)變與FEM計(jì)算所得到的結(jié)果相近，即工件的自由段由于缺少約束，其變形比較大，最大的應(yīng)變?cè)谧杂啥蔚闹虚g部位，鍛件下端由于凹模的固定作用，故其應(yīng)變很小。

4 反擠壓工藝過(guò)程模擬

金屬坯料三維反擠壓的成形模型及尺寸如圖8所示。工件材料參數(shù)如下：彈性模量E=117 MPa，初始屈服應(yīng)力σ0=400 MPa，硬化模量H=100 MPa，密度ρ=8930 kg/m3，泊松比μ=0.35。凸模和凹模簡(jiǎn)化為剛體，凹模固定，凸模以v=5 m/s的速度向下擠壓工件。采用庫(kù)侖摩擦，摩擦因數(shù)取0.12。為驗(yàn)證MPM仿真金屬大變形成形過(guò)程的有效性，暫不考慮變形過(guò)程中的材料失效。

圖8 金屬反擠壓成形過(guò)程的物理模型和尺寸

工件在不同壓下量Δh時(shí)，F(xiàn)EM與MPM的計(jì)算結(jié)果分別如圖9所示。

(a)FEM(Δh=2 mm) (b)MPM(Δh=2 mm)

(c)FEM(Δh=4 mm) (d)MPM(Δh=4 mm)圖9 不同Δh下反擠壓模擬結(jié)果

(a)Δh=4.5 mm (b)Δh=5.5 mm

從圖9可以看出，當(dāng)垂直壓下量Δh較小時(shí)，MPM和FEM均能得到很好的結(jié)果。但隨反向擠壓垂直壓下量Δh的增大，即當(dāng)壓下量為4 mm時(shí)，有限元法計(jì)算的網(wǎng)格已發(fā)生了嚴(yán)重畸變，需要進(jìn)行重構(gòu)才能繼續(xù)求解，而MPM的計(jì)算卻不受影響。若繼續(xù)增大壓下量，MPM依然能計(jì)算，圖10所示為后續(xù)擠壓時(shí)的一些變形結(jié)果，該結(jié)果與文獻(xiàn)[19]相同。

(c)Δh=6 mm (d)Δh=7 mm圖10 MPM法模擬反擠壓過(guò)程結(jié)果

從圖10可看到，即使工件變形非常大，MPM也能夠在不需要任何特殊處理的情況下完成整個(gè)反擠壓過(guò)程的模擬，并且其結(jié)果保持光滑，這說(shuō)明MPM在計(jì)算金屬大變形時(shí)具有較好的優(yōu)越性。

5 結(jié)論

(1)無(wú)論是頂鐓成形還是反擠壓成形，MPM均能得到非常好的仿真結(jié)果。

(2)由于MPM在背景網(wǎng)格的計(jì)算中采用有限元形函數(shù)，故MPM能克服傳統(tǒng)無(wú)網(wǎng)格法中施加本質(zhì)邊界條件困難的問(wèn)題。

(3)MPM能克服有限元法在模擬材料大變形時(shí)所出現(xiàn)的網(wǎng)格畸變，當(dāng)材料發(fā)生非常大的變形時(shí)MPM不需要任何特殊處理就能一次性完成整個(gè)計(jì)算過(guò)程，說(shuō)明MPM在模擬金屬塑性大變形問(wèn)題時(shí)具有較大的優(yōu)勢(shì)。

[1] 劉玉紅, 李付國(guó), 吳詩(shī)惇. 體積成形數(shù)值模擬技術(shù)的研究現(xiàn)狀及發(fā)展趨勢(shì)[J].航空學(xué)報(bào), 2002, 23(6): 547-551. Liu Yuhong, Li Fuguo, Wu Shichun. Numerical Simulation for Bulk Forming Process: Its Current Status and Future Progress[J]. Acta Aeronautica et Astronautica Sinica, 2002, 23(6): 547-551.

[2] 齊會(huì)萍, 李永堂, 華林, 等. 環(huán)形零件輾擴(kuò)成形工藝研究現(xiàn)狀與發(fā)展趨勢(shì)[J]. 機(jī)械工程學(xué)報(bào), 2014, 50(14): 75-80. Qi Huiping, Li Yongtang, Hua Lin, et al. Research Status and Developing Trends on the Rolling Forming Process of Ring Parts[J]. Journal of Mechanical Engineering. 2014, 50(14): 75-80.

[3] 朱蒙蒙，謝桂蘭，曹尉南，等. 無(wú)網(wǎng)格SPH法在金屬鐓粗成形中的應(yīng)用[J]. 熱加工工藝,2012,41(5): 20-23. Zhu Mengmeng，Xie Guilan，Cao Weinan，et al. Application of Smooth Particle Hydrodynamic in Metal Forging[J]. Hot Working Technology,2012,41(5): 20-23.

[4] 劉永輝，陳軍. 任意模具形狀的金屬三維體積成形過(guò)程剛塑性無(wú)網(wǎng)格RKPM分析方法[J]. 塑性工程學(xué)報(bào), 2008, 15 (3): 105-109. Liu Yonghui, Chen Jun. Rigid-plastic Reproducing Kernel Particle Method for Three-dimensional Bulk Metal Forming with Arbitrarily-shaped Dies[J]. Journal of Plasticity Engineering,2008, 15(3): 105-109.

[5] 康永林，朱國(guó)明，陳偉. LS-DYNA無(wú)網(wǎng)格伽遼金方法在軋制過(guò)程三維仿真分析中的應(yīng)用[J]. 塑性工程學(xué)報(bào), 2007, 14(4): 140-146 Kang Yonglin, Zhu Guoming, Chen Wei. EFG Method’s Application in the 3D Simulation of Rolling Process in LS-DYNA[J]. Journal of Plasticity Engineering,2007, 14(4): 140-146.

[6] 崔青玲，李長(zhǎng)生，劉相華，等. 無(wú)網(wǎng)格法及其在金屬塑性成形中的應(yīng)用[J]. 塑性工程學(xué)報(bào)，2005, 12(1): 38-42. Cui Qingling，Li Changsheng，Liu Xianghua，et al. Meshless Methods and Its Application in Metal Forming[J]. Journal of Plasticity Engineering，2005, 12(1): 38-42.

[7] Ma S, Zhang X, Qiu X M. Comparison Study of MPM and SPH in Modeling Hypervelocity Impact Problems[J]. International Journal of Impact Engineering，2009, 36(2): 272-282.

[8] Smolinski P, Palmer T. Procedures for Multi-time Step Integration of Element-free Galerkin Methods for Diffusion Problems [J]. Computers & Structures，2000, 77(2): 171-183.

[9] Sulsky D, Chen Z, Schreyer H L. A Particle Method for History-dependent Materials [J]. Computer Methods in Applied Mechanics and Engineering,1994, 118(1/2): 179-196.

[10] 廉艷平，張帆，劉巖，等. 物質(zhì)點(diǎn)法的理論和應(yīng)用[J].力學(xué)進(jìn)展,2013,43(2): 237-264. Lian Yanping,Zhang Fan,Liu Yan,et al.Material Point Method and Its Applications[J].Advances in Mechanics,2013,43(2): 237-264.

[11] 馬上，張雄，邱信明. 超高速碰撞問(wèn)題的三維物質(zhì)點(diǎn)法[J]. 爆炸與沖擊,2006，26(3): 273-278. Ma Shang，Zhang Xiong，Qiu Xinming. Three-dimensional Material Point Method for Hypervelocity Impact[J]. Explosion and Shock Waves,2006，26(3): 273-278.[12] 王宇新，陳震，孫明. 滑移爆轟問(wèn)題無(wú)網(wǎng)格MPM法數(shù)值模擬[J].力學(xué)與實(shí)踐,2007, 29(3): 20-25. Wang Yuxin，Cheng Zhen，Sun Ming. Numerical Simulation of Slippage Detonation by Material Point Method-MPM[J]. Mechanics in Engineering, 2007, 29(3): 20-25.

[13] Guo Y J, Nairn J A. Three-dimensional Dynamic Fracture Analysis Using the Material Point Method[J]. Computer Modeling in Engineering and Sciences,2006, 16(1): 141-155.

[14] Wieckowski Z. The Material Point Method in Large Strain Engineering Problems[J]. Computer Methods in Applied Mechanics and Engineering,2004, 193(39/41): 4417-4438.

[15] Carter M M, Pedro A, Peter M H, et al. Simulation Granular Column Collapse Using the Material Point Method[J]. Acta Geotechnina, 2014, 10(1):101-116.

[16] Nair A,Roy S.Implicit Time Integration in the Generalized Interpolation Material Point Method for Finite Deformation Hyperelasticity[J]. Mechanics of Advanced Material and Structures, 2012, 19(6): 465-473.

[17] 張雄,王天舒.計(jì)算動(dòng)力學(xué)[M].北京:清華大學(xué)出版社,2007.

[18] Xiong S, Li C S, Rodrigues J M C, et al. Steady and Non-steady State Analysis of Bulk Forming Processes by the Reproducing Kernel Particle Method[J]. Finite Elements in Analysis and Design,2005, 41(6): 599-614.

[19] Li G, Belytschko T. Element-free Galerkin Method for Contact Problems in Metal Forming Analysis[J]. Engineering Computations, 2001, 18(1/2): 62-78.

(編輯 陳 勇)

Simulation of Metal Bulk Forming Based on MPM

Xie Guilan Yu Chao Gong Shuguang Tian Jie Hu Jundi

Xiangtan University,Xiangtan,Hunan,411105

Aiming at the problems that mesh distortions in finite element method(FEM) and difficulties of applying the essential boundary conditions in the traditional mesh-less method when doing the analysis of metal bulk forming process, a simulation model of metal bulk forming process was presented based on MPM. Because the shape function of the finite element method was used on the background grid nodes, MPM might solve the difficulty problems in dealing with the boundary conditions effectively. Two metal bulk forming processes(heading of cylindrical billets and backward extrusion) were performed, compared the results of experiments with those of finite element method, it is shown that MPM may eliminate the mesh distortions, and may get more accurate results than that of the finite element method when large deformations occur. The conclusions may provide guidance for the MPM to be applied in the large elastic-plastic deformation analyses.

material point method(MPM); large elastic-plastic deformation; metal bulk forming; numerical simulation

2016-01-19

國(guó)家自然科學(xué)基金資助項(xiàng)目(51475403)

TG302

10.3969/j.issn.1004-132X.2016.22.019

謝桂蘭，女，1966 年生。湘潭大學(xué)機(jī)械工程學(xué)院教授。主要研究方向?yàn)樾虏牧狭W(xué)性能。發(fā)表論文 30 余篇。于 超，男， 1990年生。湘潭大學(xué)機(jī)械工程學(xué)院碩士研究生。龔曙光，男，1964年生。湘潭大學(xué)機(jī)械工程學(xué)院教授。田 杰，男, 1988 年生。湘潭大學(xué)機(jī)械工程學(xué)院博士研究生。胡駿迪，男，1990年生。湘潭大學(xué)機(jī)械工程學(xué)院碩士研究生。