模拟滚弯非稳态过程的虚拟载荷法

郑子君，王 洪

(重庆理工大学 机械工程学院，重庆 400040)

1 引 言

滚弯是一种常见的型材和板料弯曲加工工艺。其涉及的机械结构简单，便于连续生产，使用一套辊轮就可以加工多种曲率的产品，因此有着独特的优势。以平面对称式三辊滚弯为例，其基本工艺流程为[1,2],将工件一段置于滚弯机中(图1(a))，在静态弯曲阶段，调整辊轮至设计位置，夹紧工件(图1(b))，然后在动态弯曲阶段，通过送料辊转动，带动工件在水平方向连续地通过滚弯机并发生塑性弯曲(图1(c))。经过一段时间后，滚弯机内各空间点处的工件构型将达到一个稳定状态，辊轮反力与输出曲率均趋于定值。

对于滚弯的稳定状态的研究，已有较多的理论模型，如三点弯模型[3，4]、相切圆模型[5，6]、曲梁模型[7]及曲率积分法[8,9]等。此外，有限元模拟也是重要的研究手段，如Feng等[10]模拟了非对称式三辊滚弯，并进行了实验验证；Ktari等[11]通过显式动力学模拟得到了稳态输出曲率和辊轮位置的关系；Fu等[4]采用平面应变单元来模拟稳态滚弯，在小变形情形下结果与三点弯模型的理论解吻合良好；Shin等[12]比较了平面单元和梁单元的模拟结果，指出平面单元表现更好；Kim等[5]采用板壳单元模拟了薄板滚弯，提出采用这种单元时可按15%对输出曲率进行修正。

对滚弯非稳态过程的研究，有助于控制直边效应和提高产品质量，主要通过有限元模拟进行，如Kagzi等[2]对圆锥滚弯过程中各辊轮受力的波动情况进行了仿真；Tran等[13]比较了有限元模拟和实际测量的板料表面应变演化规律；Groth等[14,15]模拟了辊轮位置发生调整后，输出曲率和工艺力出现波动并最终达到稳态的过程。

图1 对称式三辊滚弯工艺

上述研究在模拟滚弯时往往采用相似的有限元建模方案，即工件采用柔性体建模，辊轮采用刚体建模；工件和辊轮间定义接触；通过指定辊轮的运动，带动工件产生位移和塑性变形。然而实践中发现，这种方案的效率不高，通常明显慢于同等尺度的有模弯曲的模拟。原因主要有两点，(1) 滚弯过程中工件不断平移通过变形区(图 1(a)AC间)，必须在变形区外预留足够长的工件，导致网格数量过多；(2) 由于辊轮与工件间接触面积小，且有切向相对滑动，很容易发生穿透和接触状态突变，因此处理接触时需采用小时间步长、大接触刚度和大检索半径，降低了接触算法的效率[16,17]。

受计算流体力学的启发，本文基于欧拉观点描述三辊滚弯机变形区内的型材工件状态变化。把滚弯机内的变形区视作一维流场，该空间内的工件挠度和转角视为流场变量。其优点在于，工件未进入滚弯机的部分不需要建立网格，只需提供其在出入口处的状态；网格的拓扑形态和结点坐标均是固定的，使得结点和辊轮间的接触搜索得到简化。

然而，采用欧拉观点时，单元方程形式往往和通常的拉格朗日列式有明显的区别，增加了开发和使用的成本。因此本文针对小曲率滚弯工艺中工件移动的特点，提出一个虚拟载荷的概念，使得基于欧拉观点的滚弯控制方程可以直接套用较为成熟的拉格朗日梁单元列式。

2 虚拟载荷法

为简单起见，沿用文献[3,4]对工件平面滚弯模型作的假设，在变形区内的工件始终处于小挠度状态；辊间距L显著大于工件横截面的尺寸h；重力、剪力和惯性力对工件的变形影响可以忽略。此时工件适用欧拉-伯努利直梁模型。

建立如图1所示坐标系。记t时刻某空间坐标x处的曲率为κt(x)，材料和截面参数集合为向量Ct(x)，塑性内变量集合为向量Pt(x)。根据欧拉梁的本构关系，弯矩由变形程度、截面参数以及塑性内变量确定。由内外力增量的平衡得

ΔMe x(x)=M[κt + Δ t(x),Ct + Δ t(x),Pt + Δ t(x)]-

M[κt(x),Ct(x),Pt(x)]

(1)

设经过微小时间步长Δt，位于x处的材料截面由原本位于x+s处的材料截面受辊轮带动平移而来，即有Ct + Δt(x)=Ct(x+s)；又由步长Δt很小，可以假设任意材料点的应力路径是简单的，则计算弯矩时，只需要知道该时间步内材料截面的初始塑性内变量值Pt(x+s)。于是式(1)变为

ΔMe x(x)=M[κt + Δt(x),Ct(x+s),Pt(x+s)]-

M[κt(x),Ct(x),Pt(x)]

(2)

在小变形情形下，s为与坐标x无关的常数，物理意义为该时间步内的滚弯机的进料量。式(2)可进一步整理为

ΔMe x(x)={M[κt + Δt(x),Ct(x+s),Pt(x+s)]-

M[κt(x),Ct(x+s),Pt(x+s)]}+

{M[κt(x),Ct(x+s),Pt(x+s)]-

M[κt(x),Ct(x),Pt(x)]}

(3)

式(3)第二项的意义是材料水平流动后，为维持t时刻构型不变所需要的额外载荷，记为-ΔMa r t(x)。而式(3)第一项的物理意义是由变形引起的弯矩变化。若定义割线弯曲刚度ks为

ks(x)={M[κt + Δt(x),Ct + Δt(x),Pt + Δt(x)]-

M[κt(x),Ct + Δt(x),Pt + Δt(x)]}/

[κt + Δt(x)-κt(x)]

(4)

则平衡方程(3)可写为

ks(x)Δκ(x)=ΔMe x(x)+ΔMa r t(x)

(5)

可以看出，式(5)虽然基于空间坐标x得出，但与基于材料坐标的欧拉梁方程相比，差别仅是右端项除了实际弯矩增量ΔMe x(x)以外，还额外多了一项ΔMa rt(x)。它是空间构型不变时，由于材料流动而导致的内力不平衡量。由于该项在t+Δt时刻是已知的，故可以作为载荷项来处理。引入额外载荷项来考虑材料流动影响的方法可称为虚拟载荷法。从定义容易看出，如果梁的截面参数C和塑性内变量P在各处是均匀的，则虚拟载荷项自动消失。

3 有限元实施方案

考虑用有限元方法求解控制方程(5)。

(1) 网格方案。将两个底辊间的变形区空间(图 1中AC间)等分成若干单元。由于有限元程序通常只在高斯积分点处存储截面数据，为了提高精度，取工件每步进料长度s恰等于单元的尺寸l。则t+ Δt时刻各单元积分点的参数即为其上游单元积分点在t时刻的参数；而入口处单元的参数为从进口新流入工件的初始参数。

(2) 截面参数集合与本构关系。简单起见，材料采用双线性随动强化弹塑性模型。把横截面在厚度上分为若干层，则高斯积分点处的截面参数集合C包括初始曲率κ0，杨氏模量E，屈服强度σs，硬化模量Et，各层形心到中性轴的距离ym(m为层号)和惯性矩Im；塑性内变量集合P则包含各层形心处的轴向塑性应变εp m。则弯矩与曲率和截面参数的关系为

(6)

式(6)在线弹性时退化为弹性欧拉梁方程M=EIκ。

(3) 单元平衡方程。原则上而言，虚拟载荷法的单元列式可基于现有的任意梁单元方程推导[18-20]。简单起见，本文单元位移场采用两结点三次Hermite插值，即取挠度w和转角θ为结点变量，形函数N={1-3ξ2+2ξ3,(ξ-2ξ2+ξ3)l,3ξ2-2ξ3,(ξ3-ξ2)l}，则单元列式推导与经典梁单元[20]完全相同。因此本文不加推导地给出单元平衡方程的增量格式：

(7)

(8)

式中ξg和Wg为积分域为(0,1)的第g个高斯积分点和对应的权重。与结点外力载荷对应的右端项为

(9)

式中 Δq，Δf和Δm分别为单元受分布力、集中力和集中力偶的增量。而与虚拟载荷对应的右端项为

{M[κt(ξg),Ct(ξg),Pt(ξg)]-

(10)

(4) 方程的组装和约束。组装与约束过程和传统有限元方法完全相同。由于欧拉网格中结点不会随构型改变而移动，故接触处理可以简化。更进一步，如果忽略辊轮尺寸，即将模型简化为三点弯[3,4]，则模具对工件的作用简化为结点强制位移条件，可以通过矩阵代数变形直接施加[20]。

(5) 非线性方程的求解。当工件发生弹塑性变形时，式(8)的割线弯曲刚度不是常数。本文采用割线迭代法，对割线刚度进行迭代更新求解，直至位移增量收敛。迭代结束后，更新各高斯积分点处的塑性内变量Pt +Δt(ξg)，并进入下一时间步的求解。

在实际计算中，若遇到割线法迭代收敛困难，可以将每次平移产生的虚拟载荷分多个子步施加，并动态调整子步步长。

(6) 程序流程图。用虚拟载荷法模拟滚弯过程的流程如图2所示。

图2 采用虚拟载荷法模拟滚弯流程

4 算 例

取对称式三辊滚弯机的辊距L=400 mm，辊径Dt=Db=50 mm。中辊下压量固定为d=3.8 mm。工件的横截面为h=10 mm的矩形，初始时工件平直。材料采用理想弹塑性模型，杨氏模量E=200 GPa，屈服强度σs=200 MPa。总滚弯长度为600 mm。

采用Mathematica 11.0编写虚拟载荷法程序；底辊AC间的空间划分为 200个单元；单元内用两点高斯积分，每个积分点在厚度方向上等分为20层；动态弯曲阶段每步进料长度s=2 mm。作为对照，在ANSYS 18.0中分别用平面应力单元和板壳单元建立两个有限元模型，网格尺寸完全对照虚拟载荷法的模型；辊轮和工件间粗糙接触，左底辊轮主动转动带动工件进入滚弯机，其余辊轮被动转动；设置动态时间步长上限使工件每次平移不超过2 mm。此外，还通过曲率积分法[10]得出静态弯曲状态和最终稳定状态下的理论解作为参考。

各模拟方案均在CPU为i7-6830HQ，内存为16 GB的计算机上以串行模式进行。ANSYS平面应力模型耗时约575 min，板壳模型耗时约 61 min，而虚拟载荷法仅耗时15 min，明显快于传统方案。

求解得到的中辊的约束反力和出口处的曲率随滚弯长度S的变化如图3所示。可以看出各模型的结果曲线非常相似。大约在滚过1.5倍跨长后，各解基本达到稳态。对比曲率积分法给出的稳态解可知，虚拟载荷法与理论解最接近。板壳模型与其他模型相比，变形量明显偏小，下压力偏大，与文献[11,12]结论相符，这是因为板壳模型中存在弯曲面外的正应力，而其他模型均未计面外应力。平面应力模型的塑性变形偏大，因为该模型能考虑中辊下方的应力集中，材料更容易发生塑性屈服。

图3 中辊工艺力与输出曲率随滚弯长度的变化

图4给出了达到稳态后各模型在跨内的塑性应变分布。可以看出，平面应力模型受应力集中的影响，塑性应变关于中性层不对称，受压侧偏大。由于在计算中该模型的接触误差难以保持恒定，所得的结果始终有微小波动，故应变分布和输出曲线不光滑(图3和图4(b))。此外该模型还偶尔会出现接触穿透，如图4的箭头处，塑性应变突然下降。板壳模型的塑性应变分布虽然最为光滑，但数值明显低于其他模型。而虚拟载荷法从原理上杜绝了偶然穿透现象，能给出光滑的结果曲线，并且数值上与曲率积分法给出的理论解最为接近。

图4 稳定状态下变形区内轴向塑性应变的分布

达到稳定状态后变形区出口处的横截面上的轴向残余应力分布如图5所示。可以看到虚拟载荷法与理论解基本吻合；平面应力模型在工件顶部的残余应力明显大于理论解，在工件底部则略偏小；而板壳模型给出的残余应力则整体偏小。

图5 稳定状态下残余应力在厚度方向上的分布

在上述算例的基础上逐渐增加中辊下压量，得到稳态下的辊轮反力和曲率随下压量变化的关系如图6所示。随着变形区内的工件构型逐渐偏离小挠度情形，各有限元模型的误差增加，预测的输出曲率逐渐小于曲率积分得到的理论值。相比之下，平面应力模型得到的输出曲率误差最小，而虚拟载荷法得到的下压力与理论值吻合最好。

图6 稳定状态下中辊工艺力和输出曲率与下压量的关系

5 结 论

针对三辊式滚弯的数值模拟，本文基于弹塑性欧拉梁模型提出了一种虚拟载荷法。相比通常的有限元模拟方案，其具有如下特点。

(1) 方法基于欧拉观点，便于处理定义在空间坐标上的力和位移载荷。

(2) 只对变形区内的空间划分网格，减少了网格数量，计算效率高。

(3) 不需要建立模具网格，能避免接触面近似和偶然穿透导致的误差。

本文方法的主要限制是只适于小曲率情形。对大曲率滚弯，适用的梁板理论、材料流动追踪和接触处理都将变得更加复杂，仍需进一步的研究。