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

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

?

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

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

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

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

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

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

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

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

1 基本方程

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

1.1 速度場(chǎng)移動(dòng)最小二乘近似

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

(1)

式中,uI為節(jié)點(diǎn)速度向量;U為廣義速度向量;N為離散節(jié)點(diǎn)數(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)變問(wèn)題的幾何方程為:

(8)

式中:

(9)

微分算子

(10)

1.4 控制方程

平面問(wèn)題的控制方程為如下。

平衡方程:

(11)

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

(12)

速度邊界條件:

(13)

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

(14)

(15)

(16)

式中:

(17)

(18)

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

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

(19)

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

(20)

(21)

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

(22)

改寫為矩陣形式

(23)

在迭代過(guò)程中,

(24)

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

3 數(shù)值算例

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

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

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

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

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

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

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

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

4 結(jié)論

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

[1] 張雄,劉巖.無(wú)網(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)華,王付杰,等.斜軋穿孔過(guò)程的剛塑性無(wú)網(wǎng)格法數(shù)值模擬[J].塑性工程學(xué)報(bào),2013,20(3):16-21.

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

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

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

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

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

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

TG316

A

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