基于最小二乘無網(wǎng)格法的金屬變形過程模擬

周 研，雙遠(yuǎn)華，趙春江，茍毓俊

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

?

基于最小二乘無網(wǎng)格法的金屬變形過程模擬

周 研a,b，雙遠(yuǎn)華b,c，趙春江b,c，茍毓俊b,c

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

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

最小二乘無網(wǎng)格法；移動最小二乘近似；剛塑性；塑性變形

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

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

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

1 基本方程

為了使問題簡化，本文以平面應(yīng)變問題為例。

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

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

(1)

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

(2)

(3)

1.2 應(yīng)力-應(yīng)變率關(guān)系

假設(shè)材料為理想剛塑性,滿足體積不可壓縮條件,材料的應(yīng)力-應(yīng)變率關(guān)系為:

(4)

(5)

(6)

(7)

1.3 幾何方程

平面應(yīng)變問題的幾何方程為:

(8)

式中:

(9)

微分算子

(10)

1.4 控制方程

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

平衡方程:

(11)

應(yīng)力邊界條件:

(12)

速度邊界條件:

(13)

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

(14)

(15)

(16)

式中:

(17)

(18)

2 剛塑性最小二乘無網(wǎng)格法

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

(19)

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

(20)

(21)

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

(22)

改寫為矩陣形式

(23)

在迭代過程中,

(24)

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

3 數(shù)值算例

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

圖2 無網(wǎng)格節(jié)點與有限元單元Fig.2 Meshless nodes and FEM elements

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

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

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

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

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

A:0.35；B:0.40；C:0.45；D:0.50；E:0.55圖5 30%壓下量時無網(wǎng)格法與有限元所得等效應(yīng)變分布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%壓下量時無網(wǎng)格法與有限元所得等效應(yīng)變分布Fig.6 Contours of equivalent strain at 50% reduction by Meshless and FEM

4 結(jié)論

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

[1] 張雄,劉巖.無網(wǎng)格法[M].北京:清華大學(xué)出版社,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] 胡建華,雙遠(yuǎn)華,王付杰,等.斜軋穿孔過程的剛塑性無網(wǎng)格法數(shù)值模擬[J].塑性工程學(xué)報,2013,20(3):16-21.

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

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

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

[11] 張雄,胡偉,潘小飛,等.加權(quán)最小二乘無網(wǎng)格法[J].力學(xué)學(xué)報,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

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

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

TG316

A

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