基于最小二乘无网格法的金属变形过程模拟

周 研，双远华，赵春江，苟毓俊

(太原科技大学 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