DOI:10.19476/j.ysxb.1004.0609.2002.06.007

金属三维挤压成形过程数值模拟的若干关键技术

陈军 张向 阮雪榆

  上海交通大学模具CAD国家工程研究中心  

  上海交通大学模具CAD国家工程研究中心 上海200030  

摘 要：

推导了采用六面体等参单元网格划分的金属三维挤压成形过程数值模拟中的摩擦力功率泛函在局部坐标系的一阶和二阶偏导数的表达式 , 提出了有效的三维六面体网格划分和网格质量控制技术处理方法。以此为基础 , 建立了三维刚塑性有限元数值模拟系统和六面体网格再划分系统。方坯料反向挤压过程数值模拟的结果与实验非常吻合 , 验证了所推导的公式和开发的系统 ;圆坯料侧向挤压过程的模拟验证了提出的网格质量控制技术的第三种方法。

关键词：

挤压工艺;数值模拟;三维刚塑性有限元法;

中图分类号： TG371

收稿日期：2002-03-05

Research of couple of key techniques in numerical simulation of 3D metal extrusion processes

Abstract：

The derivatives of the first order and the second order of the calculus of the function of the frictional power induced by the boundary friction force were deduced, which will be the foundation in simulating 3-D extrusion processes. In addition, the effective methods for initial meshing and mesh quality control were presented. The self-developed simulation system and remeshing system coupled with the commercial system DEFORM3D were used to simulate several extrusion processes, and the simulation results fit well with the experimental results, which verified the reliability of the proposed methods.

Keyword：

extrusion; numerical simulation; 3-D rigid-plastic finite element formulations;

Received： 2002-03-05

冷挤压技术是一种高效的少无切削加工技术, 具有较高的生产效率和材料利用率, 而且成形后可以获得较高的综合力学性能和较好的制品表面质量, 因而该技术在运载工具制造和机床制造业中具有广泛的应用。 众所周知, 冷挤压时金属的变形抗力非常大, 因而传统的模具设计和试模, 既耗时耗费, 又具有很大的危险性。 基于刚塑性有限元法的金属挤压成形过程数值模拟, 可以大幅度减少设计的反复次数, 同时也大大降低了试模的危险性 ^{[1,2,3,4,5,6]} , 因此具有非常广阔的应用前景。

本文作者针对采用综合计算精度和效率比较高的三维六面体单元数值模拟挤压成形过程所涉及到的摩擦功率泛函和网格划分问题, 推导了边界摩擦力引起的功率泛函的一阶和二阶偏导数的表达式, 并提出了有效的六面体网格划分/再划分技术。 综合作者提出的其它技术, 建立了三维刚塑性有限元数值模拟系统和六面体网格再划分系统, 并对几种挤压成形过程进行了数值模拟和实验验证。

1刚塑性有限元数值模拟的理论基础

金属挤压成形过程中塑性变形远大于弹性变形, 因此常采用刚塑性本构关系。 在用于刚塑性有限元数值模拟的几种方法中, 修正的罚函数法在保证计算精度的条件下, 具有较高的计算效率 ^[7,8] 。 基于罚函数法的泛函可以按照下列模式描述。

设刚塑性/刚粘塑性体的体积为V, 表面积为S, 在S_F面上给定面力F_i, 在S_V面上给定速度u_i, 则在满足几何条件、 体积不可压缩条件和速度边界条件的一切运动许可的速度场中, 真实解使泛函式 (1) 取极值:

$\begin{array}{l}\pi ={\displaystyle \int {\displaystyle \int {\displaystyle {\int }_V\stackrel{?}{\sigma }}}}?\stackrel{?}{\epsilon }\text{d}v+\frac{\alpha }{2}{\displaystyle \int {\displaystyle \int {\displaystyle {\int }_V(}}}\stackrel{?}{\epsilon }{}_{ij}?\delta {}_{ij}){}^2\text{d}v-\\ {\displaystyle {\int }_{S{}_{\text{F}}}F}{}_{i}u{}_i\text{d}s\text{?}\text{?}\text{?}(1)\end{array}$

式中 α—惩罚因子。

基于六面体等参单元网格划分的刚塑性有限元法求解时, 将上式采用矩阵表达式, 然后线性化, 通过摄动法求解:

[S]{u}={R} (2)

式中 [S]—刚度矩阵, {u}—节点的速度或速度增量列阵, {R}—节点的载荷列阵 。

2摩擦力功率泛函及其一阶、 二阶偏导数表达式

式 (1) 中的最后一项中, 外力F_i即为边界摩擦力。 反正切函数摩擦力模型, 既考虑了分流点两侧金属所受摩擦力的连续变化, 而且二阶可导, 因而对摩擦力的描述最为合理和完备。 反正切函数摩擦力模型的数学表达式为 ^[8]

$f=-mkt=-\frac{2}{\pi }mk\arctan (\frac{w{}_{\text{s}}}{A})\text{?}\text{?}\text{?}(3)$

式中 m—摩擦因子, k—最大剪应力, w_s —工具与材料接触面上的相对滑动速度, A—比模具速度小几个数量级的正常数, 一般取10^-4～10^-5。

当采用反正切函数模型时, 摩擦力与速度有关, 因此泛函式 (1) 中的摩擦项应为

π_f=-∫∫_SFf_iv_ids (4)

对于三维成形过程的数值模拟, 当六面体边界单元的边界表面与模具的表面接触时, 由于二者存在相对滑动而产生摩擦力。 图1所示为仅考虑4个边界节点同时接触模具表面的情况。

设模具的速度分量为{v_d1, v_d2, v_d3}, 单元触模节点的速度分量为{u_1i, u_2i, u_3i}, (i=1, 2, 3, 4) 。 该节点的局部坐标系相对于整体坐标系的方向余弦矩阵 (3×3) 为[T]_i, (i=1, 2, 3, 4) , 则该节点在与模具该位置的切平面上有两个方向的滑动速度分量u_i′和v_i′。 因此, 接触面上任意一点的滑动速度可以表示为

图1 边界单元与模具边界的接触示意图

Fig.1 Contact of boundary element to die boundary

$w{}_{\text{s}}={\displaystyle \underset{i=1}{\overset{4}{\sum }}{\rm N}}{}_{i}w{}_{\text{s}i}\text{?}\text{?}\text{?}(5)$

式中 ${\rm N}{}_i=\frac{1}{4}(1+\xi \xi {}_i)(1+\zeta \zeta {}_i)(i=1,2,3,4),w{}_{\text{s}i}=\sqrt{u{}_{\text{s}i}^2+v{}_{\text{s}i}^2}$ 。 其中,

$\{\begin{array}{l}u{}_{\text{s}i}=u{}_{i^{\prime }}-u{}_{d^{\prime }}\\ v{}_{\text{s}i}=v{}_{i^{\prime }}-v{}_{d^{\prime }}\end{array}$

而

$\{\begin{array}{l}u{}_{i^{\prime }}={\displaystyle \underset{j=1}{\overset{3}{\sum }}u}{}_{ji}?{\rm T}{}_{i2j}\\ v{}_{i^{\prime }}={\displaystyle \underset{j=1}{\overset{3}{\sum }}u}{}_{ji}?{\rm T}{}_{i3j}\end{array}\text{?}\text{?}\{\begin{array}{l}u{}_{d^{\prime }}={\displaystyle \underset{j=1}{\overset{3}{\sum }}v}{}_{dj}?{\rm T}{}_{i2j}\\ v{}_{d^{\prime }}={\displaystyle \underset{j=1}{\overset{3}{\sum }}v}{}_{dj}?{\rm T}{}_{i3j}\end{array}$

式中 u_d′, v_d′—节点所对应的模具上的点在切平面内的两个滑动速度分量; T_i2j —节点i局部坐标系中第二主方向的3个余弦分量; T_i3j—节点i局部坐标系中第三主方向的3个余弦分量 (第一主方向为法向) 。

由上面的式子可得

由此可得摩擦泛函π_f关于u′_i和u′_j的一阶偏导数和二阶偏导数表达式分别为

$\{\begin{array}{l}\frac{?\pi {}_f}{?u{}_{i^{\prime }}}=-{\displaystyle \int {\displaystyle {\int }_{S{}_{\text{F}}}m}}k\frac{2}{\pi }\arctan (\frac{w{}_{\text{s}}}{A}){\rm N}{}_i\frac{u{}_{i^{\prime }}-u{}_{d^{\prime }}}{\sqrt{u{}_{\text{s}i}^2+v{}_{\text{s}i}^2}}\text{d}s\\ \frac{?\pi {}_f}{?v{}_{j^{\prime }}}=-{\displaystyle \int {\displaystyle {\int }_{S{}_{\text{F}}}m}}k\frac{2}{\pi }\arctan (\frac{w{}_{\text{s}}}{A}){\rm N}{}_j\frac{v{}_{j^{\prime }}-v{}_{d^{\prime }}}{\sqrt{u{}_{\text{s}j}^2+v{}_{\text{s}j}^2}}\text{d}s\end{array}\text{?}\text{?}\text{?}(8)$

$\frac{?{}^2\pi {}_{f^{\prime }}}{?u^{\prime }{}_i^2}=-{\displaystyle \int {\displaystyle {\int }_{S{}_{\text{F}}}m}}k\frac{2}{\pi }\frac{A}{A{}^2+w{}_{\text{s}}^2}({\rm N}{}_i\frac{u{}_{i^{\prime }}-u{}_{d^{\prime }}}{\sqrt{u{}_{\text{s}i}^2+v{}_{\text{s}i}^2}}){}^2\text{d}s+$

$\begin{array}{l}{\displaystyle \int {\displaystyle {\int }_{S{}_{\text{F}}}m}}k\frac{2}{\pi }\arctan (\frac{w{}_{\text{s}}}{A}){\rm N}{}_i^2\frac{(u{}_{i^{\prime }}-u{}_{d^{\prime }}){}^2}{\sqrt{(u{}_{\text{s}i}^2+v{}_{\text{s}i}^2){}^3}}\text{d}s-\\ {\displaystyle \int {\displaystyle {\int }_{S{}_{\text{F}}}m}}k\frac{2}{\pi }\arctan (\frac{w{}_{\text{s}}}{A})\frac{{\rm N}{}_i}{\sqrt{u{}_{\text{s}i}^2+v{}_{\text{s}i}^2}}\text{d}s\\ \text{?}\frac{?{}^2\pi {}_{f^{\prime }}}{?u{}_{i^{\prime }}?v{}_{j^{\prime }}}=-{\displaystyle \int {\displaystyle {\int }_{S{}_{\text{F}}}m}}k\frac{2}{\pi }?\\ \frac{A}{A{}^2+w{}_{\text{s}}^2}\frac{{\rm N}{}_i{\rm N}{}_j}{\sqrt{u{}_{\text{s}i}^2+v{}_{\text{s}i}^2}\sqrt{u{}_{\text{s}j}^2+v{}_{\text{s}j}^2}}?\\ (u{}_{i^{\prime }}-u{}_{d^{\prime }})(v{}_{j^{\prime }}-v{}_{d^{\prime }})\text{d}s\\ \text{?}\frac{?{}^2\pi {}_{f^{\prime }}}{?v^{\prime }{}_i^2}=-{\displaystyle \int {\displaystyle {\int }_{S{}_{\text{F}}}m}}k\frac{2}{\pi }\frac{A}{A{}^2+w{}_{\text{s}}^2}({\rm N}{}_i\frac{v{}_{i^{\prime }}-v{}_{d^{\prime }}}{\sqrt{u{}_{\text{s}i}^2+v{}_{\text{s}i}^2}}){}^2\text{d}s+\\ {\displaystyle \int {\displaystyle {\int }_{S{}_{\text{F}}}m}}k\frac{2}{\pi }\arctan (\frac{w{}_{\text{s}}}{A}){\rm N}{}_i^2\frac{(v{}_{i^{\prime }}-v{}_{d^{\prime }}){}^2}{\sqrt{(u{}_{\text{s}i}^2+v{}_{\text{s}i}^2){}^3}}\text{d}s-\\ {\displaystyle \int {\displaystyle {\int }_{S{}_{\text{F}}}m}}k\frac{2}{\text{π}}\arctan (\frac{w{}_s}{A})\frac{{\rm N}{}_i}{\sqrt{u{}_{\text{s}\text{i}}^2+v{}_{\text{s}\text{i}}^2}}\text{d}s\text{?}\text{?}\text{?}(9)\end{array}$

对以上诸式进行辛普森数值积分时, 应构造一个平面积分区域来近似代替空间曲面积分区域, 作者根据4个节点的位置和法矢, 构造了一个虚拟接触直纹面近似代替实际接触的模具表面 ^[9] 。

3三维六面体网格划分和质量控制技术

三维挤压成形过程数值模拟时, 六面体网格的划分和网格质量控制至关重要。 合理的初始网格划分可以确保可靠的模拟精度和效率。 对于截面为方形或圆形的坯料, 作者提出4种截面网格划分模式, 如图2所示, 同时还可以根据需要实现局部网格加密。

众所周知, 挤压过程伴随大的局部塑性变形, 因此初始有限元网格可能产生畸变从而使模拟的局部精度降低。 为此, 作者提出了3种网格质量控制方法:1) 基于原有网格拓扑结构的局部网格质量的光滑处理, 2) 改变局部网格拓扑结构的局部网格再划分, 3) 基于表面偏置理论的六面体网格全局再划分技术。 根据模拟时的具体情况, 综合用于六面体网格质量的控制, 并且建立了基于表面偏置理论的六面体网格全局再划分系统 ^[10] 。

4 计算实例

金属三维挤压成形过程的数值模拟还涉及到其它具体问题的处理, 如三维模具几何模型表示、 初始速度场的自动生成和速度场的迭代收敛和控制技术、 动态边界处理技术和刚塑性区的判断准则等, 采用作者在文献 [ 11] 中提出的技术和文献 [ 8] 中的技术, 在CAD/CAM软件Unigraphics平台上建立了三维有限元数值模拟系统。 采用该系统模拟了方坯料的反向挤压成形过程, 如图3所示。 其中 (a) 为变形网格, (b) 为对称截面上的变形网格实验结果与模拟结果的对照, (c) 为对角线截面上的变形网格实验结果与数值模拟结果的对照 (模拟条件和实验结果取自文献 [ 12] ) , 二者吻合很好, 证明了所推导的摩擦力功率泛函的偏导数的表达式和提出的局部网格加密技术是有效的。

图4所示为圆柱形坯料侧向挤压成形过程的模拟结果, 模拟分析器采用商业化软件DEFORM3D, 网格再划分系统采用作者自主开发的基于曲面偏置的六面体网格生成系统。 图4中 (a) 为模拟的畸变网格, (b) 为网格再划分后自动生成的六面体网格, 新生成的网格具有较好的形状, 说明了提出的网格再划分技术是有效的。

图2 常见坯料形状的二维截面网格划分模式

Fig.2 2D sectional meshing patterns for general billet shapes

图3 方坯料反向挤压的数值模拟结果与实验结果

Fig.3 Numerical simulation results and experimental results of square billet reverse extrusion

图4 圆柱形坯料侧向挤压成形过程数值模拟结果

Fig.4 Numerical simulation results of side extrusion (a) —Before remeshing; (b) —After remeshing

参考文献

[1]　AltanT , ThomasW , VazquezV , etal.Simulationofmetalformingprocesses-applicationsandfuturetrends[A ].Proceedingsofthe6thICTP [C ].Germany:Nuremberg, 1999.210.

[2]　LeeSK , KoDC , KimBM .Optimaldieprofiledesignforuniformmicrostructureinhotextrudedproduct[J].InternationalJournalofMachineTools&Manufacture, 2000, 40:14571478.

[3]　LeeCM , YangDY .Athree dimensionalsteady statefiniteelementanalysisofsquaredieextrusionbyusingau tomaticmeshgeneration[J].InternationalJournalofMa chineTools&Manufacture, 2000, 40:3347.

[4]　TangJui peng, WuWT , WaltersJ.Recentdevelop mentandapplicationsoffiniteelementmethodinmetalforming[J].JournalofMaterialsProcessingTechnology, 1994, 46:117126.

[5]　LiGuo ji, WuWT , ArrindA , etal.3Dfiniteelementanalysesofindustrialmetalforming processes[A].5thICTP[C].Columbus, Ohio, USA , 1996.479484.

[6]　张晋明, 王学文, 阮雪榆.挤压模CAD中的数值模拟系统[J].锻压技术, 1992 (6) :4144.　　ZHANGJin ming, WANGXue wen, RUANXue yu.NumericalsimulationsysteminextrusionCAD[J].Forg ing&Stamping, 1992 (6) :4144.

[7]　王祖唐, 关廷栋, 肖景容, 等.金属塑性成形理论[M].北京:机械工业出版社, 1989.　　WANGZu tang, GuanTing dong, XIAOJing rong, etal.TheoryofMetalForming[M ].Beijing:MechanicalIndustryPress, 1989.

[8]　KobayashiS , OhSI, AltanT .MetalFormingandtheFiniteElementMethod[M ].NewYork:OxfordUniver sityPress, 1989.

[9] 　陈　军.虚拟模具制造及其金属成形过程三维仿真技术研究[D].上海:上海交通大学, 1996.　　CHENJun.VirtualDie&MouldManufacturingandResearchofTheThreeDimensionalMetalFormingPro cessSimulation[D].Shanghai:ShanghaiJiaotongUni versity, 1996.

[10]　陈　军, 张　向, 阮雪榆.基于曲面偏置的六面体有限元网格再划分算法[J].上海交通大学学报, 2002, 36 (4) :2124.　　　CHENJun, ZHANGXiang, RUANXue yu.Surfaceoffset basedremeshingalgorithmofgeneratinghexahe dralelements[J].JournalofShanghaiJiaotongUniver sity, 2002, 36 (4) :2124.

[11]　陈　军, 彭颖红, 卫原平, 等.体积成形过程的三维刚塑性有限元模拟技术研究[J].上海交通大学学报, 1996, 30 (9) :127132.　　　CHENJun, PENGYing hong, WEIYuan ping, etal.Aresearchofmetalvolumeformingprocessbythree di mensionalfinitesimulation[J].JournalofShanghaiJiaotongUniversity, 1996, 30 (9) :127132.

[12]　LeeNK , YoonJH , YangDY .Finiteelementanalysisoflargedeformationbyautomaticrenodingasaweakremeshingtechniques[J].InternationalJournalofMe chanicalScience, 1992, 34 (4) :255273.