基于物質點法金屬體積成形過程的仿真

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

湘潭大學，湘潭，411105

?

基于物質點法金屬體積成形過程的仿真

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

湘潭大學，湘潭，411105

針對傳統有限元法開展金屬體積成形分析過程中存在的網格畸變，以及在常規無網格法分析中存在的本質邊界條件施加困難等問題，建立了基于物質點法金屬體積成形過程的仿真模型。因其背景網格積分采用了有限元形函數，從而有效地解決了邊界條件施加困難的問題。通過對圓柱頂鐓和反向擠壓等金屬體積成形過程的仿真，并將仿真結果與實驗和有限元結果進行了對比分析，結果顯示物質點法不僅能有效地消除網格畸變，而且在金屬發生大變形時其計算精度優于有限元法的計算精度。所得結論可為物質點法在彈塑性大變形分析中的應用提供指導。

物質點法；彈塑性大變形；金屬體積成形；數值模擬

0 引言

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

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

1 物質點法(MPM)簡介

MPM的基本思想是將連續體離散成若干個集中質量的物質點集合，并按照連續體的形狀分布在背景網格內，在整個計算過程中，物質點和背景網格節點之間應用節點形函數完成兩次映射計算并采用顯式時間積分算法求解，以實現各物質點應力和應變的更新。圖1為物質點法離散示意圖。其理論推導見文獻[9-10]。

圖1 物質點法示意圖

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

動量方程：

ρa=·σ+ρb

(1)

質量守恒方程：

(2)

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

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

(3)

式中，Ω為當前構形；σs為比應力；ts為比邊界面力。

當材料域用質點離散時，連續體的密度等于所有物質點質量之和，寫成δ函數為

(4)

式中，Np為質點總數；mj為質點j的質量；δ為Dirac Delta函數；xij為t時刻質點j的坐標。

將式(4)代入式(3)，積分形式轉變為求和形式：

(5)

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

物質點法在背景網格上計算動量方程，通過背景網格節點的形函數NI(x)來實現質點與背景網格之間的相互映射。對于三維問題背景網格采用規則八節點六面體，其形函數為

(6)

I=1,2,…,8

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

質點的任意某參數cp可通過背景網格節點相應的參數cI插值得到：

(7)

(8)

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

將式(5)中質點位移、加速度以及試函數由背景網格節點的插值替代，整理得

(9)

J=1,2,…,Ng

(10)

(11)

(12)

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

(13)

則動量方程式(9)簡化為

(14)

I=1,2,…,Ng

MPM通過在背景網格上求解動量方程即式(14)，獲得網格節點運動信息并返回質點，用來更新質點的運動狀態以及應力應變等。因為背景網格每次參加計算后都可重新生成，所以變形不會累加，網格不會發生畸變。

2 金屬成形的MPM仿真模型

(15)

σijp=Sijp+σmpδij

(16)

采用J2流動理論的徑向返回算法[17]對應力進行更新,返回映射法由彈性預測步和塑性修正步組成，在彈性預測步中有

(17)

(18)

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

(19)

(20)

(21)

(22)

(23)

(24)

(25)

根據物質點法的基本理論及金屬體積成形的仿真模型，利用FORTRAN語言編寫了金屬體積成形MPM程序。下面采用編寫的程序分別對棒材頂鐓和反向擠壓兩種工藝過程進行仿真模擬。

3 棒材頂鐓工藝過程模擬

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

圖2 圓棒頂鐓模型

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

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

圖4 文獻[18]實驗結果

對比圖3和圖4可知，本文MPM計算所得的工件變形與實驗結果相吻合，質點的排布也符合金屬流動規律，圓棒下端由于凹模的固定作用而保持不變，上端材料沿徑向流動，從而形成了頂鐓件的頭部。

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

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

從圖5可看到，MPM得到的RS、RI曲線與實驗結果相吻合，而FEM在下壓量小時結果也相吻合。但隨著下壓量增大，利用FEM得到的RS、RI曲線與實驗結果產生了偏離，這主要是因為隨著下壓量增大，坯料的變形程度增加，此時FEM網格發會生畸變，需要不斷地進行網格重構。

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

表1 最大鼓形半徑

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

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

D=60%時MPM和FEM計算得到的等效應變和應力的云圖分別見圖6和圖7。比較圖6與圖7可知，MPM計算所得的應力應變與FEM計算所得到的結果相近，即工件的自由段由于缺少約束，其變形比較大，最大的應變在自由段的中間部位，鍛件下端由于凹模的固定作用，故其應變很小。

4 反擠壓工藝過程模擬

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

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

工件在不同壓下量Δh時，FEM與MPM的計算結果分別如圖9所示。

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

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

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

從圖9可以看出，當垂直壓下量Δh較小時，MPM和FEM均能得到很好的結果。但隨反向擠壓垂直壓下量Δh的增大，即當壓下量為4 mm時，有限元法計算的網格已發生了嚴重畸變，需要進行重構才能繼續求解，而MPM的計算卻不受影響。若繼續增大壓下量，MPM依然能計算，圖10所示為后續擠壓時的一些變形結果，該結果與文獻[19]相同。

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

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

5 結論

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

(2)由于MPM在背景網格的計算中采用有限元形函數，故MPM能克服傳統無網格法中施加本質邊界條件困難的問題。

(3)MPM能克服有限元法在模擬材料大變形時所出現的網格畸變，當材料發生非常大的變形時MPM不需要任何特殊處理就能一次性完成整個計算過程，說明MPM在模擬金屬塑性大變形問題時具有較大的優勢。

[1] 劉玉紅, 李付國, 吳詩惇. 體積成形數值模擬技術的研究現狀及發展趨勢[J].航空學報, 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] 齊會萍, 李永堂, 華林, 等. 環形零件輾擴成形工藝研究現狀與發展趨勢[J]. 機械工程學報, 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] 朱蒙蒙，謝桂蘭，曹尉南，等. 無網格SPH法在金屬鐓粗成形中的應用[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] 劉永輝，陳軍. 任意模具形狀的金屬三維體積成形過程剛塑性無網格RKPM分析方法[J]. 塑性工程學報, 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] 康永林，朱國明，陳偉. LS-DYNA無網格伽遼金方法在軋制過程三維仿真分析中的應用[J]. 塑性工程學報, 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] 崔青玲，李長生，劉相華，等. 無網格法及其在金屬塑性成形中的應用[J]. 塑性工程學報，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] 廉艷平，張帆，劉巖，等. 物質點法的理論和應用[J].力學進展,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] 馬上，張雄，邱信明. 超高速碰撞問題的三維物質點法[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] 王宇新，陳震，孫明. 滑移爆轟問題無網格MPM法數值模擬[J].力學與實踐,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] 張雄,王天舒.計算動力學[M].北京:清華大學出版社,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

國家自然科學基金資助項目(51475403)

TG302

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

謝桂蘭，女，1966 年生。湘潭大學機械工程學院教授。主要研究方向為新材料力學性能。發表論文 30 余篇。于 超，男， 1990年生。湘潭大學機械工程學院碩士研究生。龔曙光，男，1964年生。湘潭大學機械工程學院教授。田 杰，男, 1988 年生。湘潭大學機械工程學院博士研究生。胡駿迪，男，1990年生。湘潭大學機械工程學院碩士研究生。