基于最小二乘無網格法的金屬變形過程模擬

周 研，雙遠華，趙春江，茍毓俊

(太原科技大學 a.機械工程學院，b.太原重型機械裝備協同創新中心，c.材料科學與工程學院，太原 030024)

?

基于最小二乘無網格法的金屬變形過程模擬

周 研a,b，雙遠華b,c，趙春江b,c，茍毓俊b,c

(太原科技大學 a.機械工程學院，b.太原重型機械裝備協同創新中心，c.材料科學與工程學院，太原 030024)

金屬塑性大變形過程可以視為典型的非線性剛塑性問題。給出了一種基于最小二乘無網格法的金屬塑性變形過程仿真方法。該方法使用移動最小二乘法構建未知場函數，利用加權最小二乘法直接由控制方程構建系統剛度矩陣。與有限元方法不同之處在于，該方法無需進行網格重分，并且不需要在節點及求解域內進行積分。最終的數值算例表明，該方法與有限元法結果基本吻合，驗證了該方法的正確性和有效性。

最小二乘無網格法；移動最小二乘近似；剛塑性；塑性變形

無網格法作為一種新興的數值計算方法，只需將問題域離散成場節點，利用近似函數求解節點未知數值，克服了有限元法對網格的依賴，在工程數值模擬應用中越來越受到人們的重視[1-2]。

國內外許多學者均在嘗試使用無網格方法解決金屬塑性成形問題。美國學者CHEN et al最早將RKPM方法應用于金屬環件延伸、冷墩粗和壓縮過程問題的研究[3-5]。胡建華、劉新等人將無網格Galerkin法應用于金屬塑性成形過程的數值仿真，該方法與有限元法在能量泛函的構建原理是相似的[6-7]。孫杰等基于無網格徑向點插值法(RPIM)對斜軋延伸過程進行了數值仿真，基于徑向基函數具有Delta函數性質，可直接施加本質邊界條件；該方法非常適合應用于金屬塑性成形過程的模擬[8]。

文獻[3-8]所使用的方法因其離散方案需要借助于背景網格進行數值積分，不屬于純無網格方法。溫宏宇等采用無網格配點法對金屬擠壓過程進行了數值模擬仿真；該方法采用核近似函數構建未知場函數，并對系統的控制方程直接進行求解[9]。配點法不需要背景網格積分，是真正的無網格法，但其解是不穩定的，需要采用特殊的穩定方案[10]。本文采用的加權最小二乘法首先由張雄等人提出并應用于求解線彈性問題[11]。本文利用控制方程殘量的加權平方和構建系統泛函，利用罰函數將體積不變條件、本質邊界條件引入系統泛函。最終的仿真結果表明，該方法與有限元法結果相吻合，驗證了本文給出方法的可行性和正確性。

1 基本方程

為了使問題簡化，本文以平面應變問題為例。

1.1 速度場移動最小二乘近似

求解域Ω內一點x的速度u(x)的移動最小二乘近似(MLS)為

(1)

式中,uI為節點速度向量;U為廣義速度向量;N為離散節點數目;MLS形函數矩陣為

(2)

(3)

1.2 應力-應變率關系

假設材料為理想剛塑性,滿足體積不可壓縮條件,材料的應力-應變率關系為:

(4)

(5)

(6)

(7)

1.3 幾何方程

平面應變問題的幾何方程為:

(8)

式中:

(9)

微分算子

(10)

1.4 控制方程

平面問題的控制方程為如下。

平衡方程:

(11)

應力邊界條件:

(12)

速度邊界條件:

(13)

將應力-應變率關系式(4)與幾何方程式(8)代入式(11)與(12),再將無網格近似函數式(1)代入,可得到式(11)至式(13)的矩陣形式如下:

(14)

(15)

(16)

式中:

(17)

(18)

2 剛塑性最小二乘無網格法

本文使用控制方程(11)至(13)殘量的加權平方和構建泛函Π,其表達式如下:

(19)

式中:N為系統離散節點總數;Nt為應力邊界Γt上的節點總數;Nu為速度邊界Γu上的節點總數;λt、λu為引入邊界條件的罰函數。將式(14)至(16)代入泛函式(19)并寫為矩陣形式:

(20)

(21)

該方程為非線性方程,需要使用Newton-Raphson迭代法進行求解。方程(21)經Taylor級數展開,忽略二階以上高階微分,取其線性部分得到

(22)

改寫為矩陣形式

(23)

在迭代過程中,

(24)

式中,n為迭代次數;β為衰減因子,β∈(0,1]。方程(23)可以經Newton-Raphson法迭代,最終得到穩定的速度場。

3 數值算例

根據上述算法,編制了相應的計算程序,對二維平面應變問題進行了模擬分析。模型的幾何尺寸見圖1。使用無網格法(Meshless method)與有限元(FEM)軟件Deform2D對同樣的工件進行了相同過程的模擬仿真,無網格離散節點模型與有限元單元見圖2。

圖2 無網格節點與有限元單元Fig.2 Meshless nodes and FEM elements

圖1中,工件尺寸20mm×20mm。模擬中使用的工件與模具間的摩擦系數μ為0.12;下模靜置,上模以1mm/s的速度勻速下壓,兩類模擬最大壓下量為總高的50%;材料為理想剛塑性,取σs=200MPa。如圖2,有限元單元數為256個,取單元節點為無網格法離散節點,節點總數共計289個;仿真過程中分別取罰函數α與λt為1×102與1×105。

圖3 壓下量30%下的幾何形狀Fig.3 Geometric shape at 30% reduction

圖3、圖4為使用加權最小二乘無網格法與有限元模擬得到的工件變形后的單元形狀與無網格節點位置,其壓下量分別為30%與50%。通過無網格節點位置與有限元單元形狀對比,兩者結果是非常接近的。從圖4中還可以觀察到,工件四角的單元形狀變化已經非常劇烈,當達到一定程度時有限元迭代將無法收斂,其將停止計算或進行網格重分；而無網格法避免了這一問題,計算將順利進行下去。

圖4 壓下量50%下的幾何形狀Fig.4 Geometric shape at 50% reduction

圖5、圖6為使用加權最小二乘無網格法與有限元模擬得到的工件變形后等效應變的等勢線圖,其壓下量分別為30%與50%。兩者計算結果曲線基本吻合。等效應變總體上從工件外部至內部逐漸增大,但最大值出現在工件的邊角區域;等效應變的變化趨勢從工件的邊界逐漸向內部減小。

A:0.35；B:0.40；C:0.45；D:0.50；E:0.55圖5 30%壓下量時無網格法與有限元所得等效應變分布Fig.5 Contours of equivalent strain at 30% reduction by Meshless and FEM

A:0.5；B:0.6；C:0.7；D:0.8；E:0.9；F:1.0圖6 50%壓下量時無網格法與有限元所得等效應變分布Fig.6 Contours of equivalent strain at 50% reduction by Meshless and FEM

4 結論

本文將加權最小二乘無網格法應用于金屬變形過程的模擬仿真,建立了金屬塑性變形過程的剛塑性無網格模型,推導了關鍵公式,并編寫了相關程序進行了金屬平面應變問題的過程仿真。與有限元結果對比表明,該方法建立的模型是正確、有效的。該方法與有限元及其他無網格法相比，有如下優勢:1)求解域節點離散，避免了有限元法由于單元劇烈變形而產生的迭代收斂困難及網格重分;2)直接由控制方程殘量的加權平方和構建系統泛函,避開了其他數值方法在構建能量泛函時使用背景網格積分及高斯積分的過程,使計算量減少,并且提高了求解的穩定性。

[1] 張雄,劉巖.無網格法[M].北京:清華大學出版社,2005.

[2]LiuGuirong.MeshFreeMethods:movingbeyondthefiniteelementmethod[M].BocaRaton,USA:CRCPress,2003.

[3]ChenJS,RoqueCMOL,PANChunhui,etal.Analysisofmetalformingprocessbasedonmeshlessmethod[J].JournalMaterialsProcessingTechchnology,1998,80-81:642-646.

[4]ChenJS,WANGHuiping,LiuWK.Meshfreemethodwithenhancedboundryconditiontreatmentformetalformingsimulation[C]∥The1999NSFDesign&ManufacturingGranteesConference,LongBeach,CA,USA,1999.

[5]YoonSP,WuCT,WANGHuiping,etal.Efficientmeshfreeformulationformetalformingsimulation[J].JournalofEngineeringMaterialsandTechnology,2001,123(4):462-467.

[6] 胡建華,雙遠華,王付杰,等.斜軋穿孔過程的剛塑性無網格法數值模擬[J].塑性工程學報,2013,20(3):16-21.

[7] 劉新,謝桂蘭,彭建新,等.基于無網格Galerkin法模擬金屬塑性成形的研究[J].機械科學與技術,2009,28(1):71-74.

[8] 孫杰,胡建華,雙遠華,等.斜軋延伸過程的無網格RPIM方法數值模擬[J].四川大學學報(工程科學版),2013,45(1):67-73.

[9] 溫宏宇,董湘懷,阮雪榆.基于配點型無網格法的金屬擠壓過程數值模擬[J].塑性工程學報,2006,13(2):57-59.

[11] 張雄,胡偉,潘小飛,等.加權最小二乘無網格法[J].力學學報,2003,35(4):425-431.

(編輯：張紅霞)

Simulation of Metal Forming Process Based on Meshless Weighted Least-square Method

ZHOU Yana,b,SHUANG Yuanhuab,c,ZHAO Chunjiangb,c,GOU Yujunb,c

(a.CollegeofMechanicalEngineering，b.CollaborativeInnovationCenterofTaiyuanHeavyMachineryEquipment，c.CollegeofMaterialScienceandEngineering,TaiyuanUniversityofScienceandTechnology,Taiyuan030024,China)

The analysis of metal forming is a nonlinear rigid-plastic problem. This paper presents a rigid-plastic meshless method for this kind of analysis based on meshless weighted least-square method (MWLS), which constructs the unknown function with moving least-square (MLS) approximation. A discrete form of the weighted-square of residuals in control equations is used to build the system of equations. Being different from the conventional finite element method (FEM), the method does not rely on the mesh generation and field interpolation. Finally the simulation results are in good accordance with the results obtained from the rigid-plastic FEM.

meshless least-square collocation method;moving least-square approximation;rigid-plastic;plastic deformation

1007-9432(2016)03-0294-05

2016-03-18

國家自然科學基金項目：薄壁管材高速旋壓工藝擬動力學特性及其可控機理研究(51375325)；山西省科技攻關項目：鎂合金管材可控張力熱連軋工藝與設備開發(20140321008-08)

周研(1983-)，男，湖南寧鄉人，博士研究生，主要從事無縫鋼管軋制工藝及設備、軋制過程的數值模擬研究，(E-mail)zy_harry@vip.163.com

TG316

A

10.16355/j.cnki.issn1007-9432tyut.2016.03.004