坯料随机局部弯曲对滚弯成形结果影响的蒙特卡洛分析

郑子君，陶裕梅

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

滚弯是对型材、板材进行弯曲加工的一种常见工艺。以对称式三辊滚弯机为例，其工作方式为将坯料置于底辊和中辊之间，中辊以给定的下压量压紧坯料，通过辊轮或送料辊转动，带动坯料连续地通过滚弯机并发生塑性弯曲。当加工过程开始一段时间后，滚弯机两底辊之间的坯料将逐渐达到一个“定常”状态，此时辊轮上的反力与流出出口的产品曲率不再随着时间变化。

由于滚弯工艺属于无模成型，回弹量大，实践中要得到目标曲率往往需要对工艺参数进行大量的试错[1]。为减少工艺开发成本，已有大量文献对坯料参数、工艺参数和产品构型之间的关系进行了探索。在理论研究方面，Basset 等[2]、Fu 等[3]基于三点弯曲模型推导了辊轮位置；Kim 等[4]、黄世军等[1]提出采用圆弧来近似滚弯机内的坯料构型；王安恒等[5]提出在大曲率滚弯中应考虑中性层的偏移；刘志芳[6]、张子骞等[7]采用曲率积分的方式求解变形区内的坯料构型。在数值模拟方面，Fu 等[3]采用平面单元来模拟稳态滚弯，结果与三点弯模型得到的理论解吻合良好；Kim 等[4]采用板壳单元模拟薄板的滚弯，提出可按15%对输出曲率进行修正；Feng 等[8]模拟了非对称式三辊滚弯，并进行了实验验证；Ktari等[9]通过显式动力学分析得到设计形状和辊轮位置的关系图谱；Shin 等[10]比较了平面单元和梁单元的模拟结果；Kagzi 等[11]对圆锥滚弯过程中非定常阶段各辊轮受力的波动情况进行了仿真；Tran 等[12]比较了有限元模拟和实验测量的板料表面应变演化规律；Groth 等[13]模拟了辊轮位置发生调整后，输出曲率和工艺力的变化过程，并对调整速率提出了建议。

前述研究都假设坯料初始是平直的。实际生产中，由于制造精度的限制和偶然因素的影响，坯料虽然宏观上整体平直，但若取很短的长度来观测，则常常有随机的局部曲率，这种随机性在通常的研究中被忽略了。这样的近似处理是否会对分析结果产生影响则未见讨论。

要分析随机参数对力学过程的影响，将有限元模拟与蒙特卡罗法结合起来是一种直观而有效的办法，其基本思路是：根据模型输入参数的统计分布规律，随机生成大量的有限元模型，求解后再对结果进行统计分析，从而得到输出的分布规律。如朱健等[14]在碳纤维布加固的厂房结构中考虑了尺寸、材料强度、载荷的随机误差对抗震能力的影响；陈力波等[15]在简支梁桥中考虑结构参数的随机误差，分析了其在地震中的易损性；金路等[16]分析了钢架中梁柱的初始侧移和直线度对整体刚度的影响；杨智勇等[17]对土质边坡的多种失效模式进行了概率分析，指出次风险滑面也应予以重视；Rafiee 等[18]、Abebe 等[19]和吴永强[20]采用蒙特卡洛方法实现了机械制造的误差六西格玛分析和鲁棒性设计。然而，基于有限元法的蒙特卡洛模拟，始终受到计算效率低的困扰[14−16,21−23]。通过简化有限元模型[15]，建立代理模型[15,18−21]，优化抽样方案[16]，引入矩方法近似计算[22]，以及利用模型对参数的梯度信息[23]等方式，可以一定程度上提高分析的效率。

对坯料有随机局部曲率的滚弯过程作蒙特卡洛模拟，可先建立整个坯料的模型并为其每个代表性单元长度(representative elementary length, REL，即曲率不发生明显改变的长度)随机生成初始曲率，然后进行有限元模拟，最后分析流出滚弯机的工件形状。当坯料足够长时，相当于蒙特卡洛模拟的样本数量足够多。注意到本文考虑的是在局部(单元级别)的随机初始曲率，模型参数的数量远远多于文献中的情形，难以使用代理模型。而传统的基于拉格朗日观点的有限元模型进行超长坯料的滚弯模拟时，时间代价很高。这是因为：一方面，模型需要为在滚弯机外的坯料划分网格，使得单元个数过多；另一方面，辊轮和坯料的接触面积过小，使得接触搜索算法效率不高；此外，初始构型的计算也增大了建模的难度。

近年大量文献尝试了采用欧拉或任意拉格朗日—欧拉网格进行金属塑性加工的模拟，发现欧拉网格在接触判断、网格数量和反畸变方面具有优势[24−28]。但欧拉网格在塑性加工中的应用主要集中在锻造[25]、挤压[26]、厚板轧制[27−28]等体积成型领域，采用三维或平面实体单元。而对于以弯曲变形为主的滚弯过程，显然基于梁/板理论建立欧拉观点的结构单元计算效率更高，但目前未见相关研究。

本文从计算流体力学受到启发，提出基于欧拉观点的滚弯模拟方案，并在经典梁单元技术的基础上，引入一个附加载荷项来处理材料在单元间流动带来的影响，从而提高有限元模拟的效率。在此基础上采用蒙特卡洛法研究正态分布的局部初始曲率对对称式三辊滚弯工艺的输出形状的影响。

1 基于欧拉网格的滚弯模拟方案

模拟塑性加工过程时，通常采用拉格朗日网格，这种网格会随着材料点的运动而移动，便于追踪材料点的受力路径。然而，对于待加工坯料很长的滚弯过程，实际新增塑性变形的却只有滚弯机内的一小部分，采用拉格朗日网格是不经济的。若将滚弯机内的空间视为一维流场，左、右底辊视为流场的入口和出口，坯料在滚弯机内的挠度视为流场变量，则可以采用流体力学中常用的欧拉网格，在任意时刻只分析滚弯机内的部分坯料。此外，欧拉网格中节点没有水平位移，坯料节点与辊轮的可能接触位置是确定的，接触搜索变得容易；在接触点处，辊轮对坯料的推动作用可简化为欧拉网格节点上的指定位移，用代数方法直接施加。

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

图1 对称式三辊滚弯工艺的示意图Fig.1 Schematic sketch of symmetrical three-rolled bending

设在t+∆t时刻，位于x处的材料截面由t时刻位于x+∆s处的材料截面被辊轮带动平移而来(图2(a)、图2(b))，于 是有C t+∆t(x)=C t(x+∆s)。又当 ∆t很短时，可认为任意材料点的应力路径是简单的，于是计算t+∆t时刻x处的弯矩时只需知道t时刻该材料截面的塑性内变量P t(x+∆s)，于是式(1)变为：

图2 基于欧拉网格的滚弯模拟原理示意图Fig.2 Schematic sketch of Eulerian-grid-based simulation of roll bending process

可以发现形式上仅仅相差了一项∆Ma，因此对式(4)进行有限元离散时，只需要在拉格朗日列式的基础上增加与内力不平衡量∆Ma对应的载荷项即可。如果采用经典的2节点4变量Hermite 插值，则式(4)的离散形式可写为：

为了实现简单，实际计算时可取单元长度l恰等于每步进料长度 ∆s，并忽略辊轮的尺寸。此时具体的模拟步骤如下：

1)初始化：将两个底辊间的变形区空间等分成若干单元，初始化各单元高斯积分点处的塑性内变量P和截面参数C。

2)模拟调辊过程：约束底辊处的两节点的垂向位移，为中辊下方节点施加指定的位移量，取附加载荷项为零，组装式(6)并进行非线性求解，更新塑性内变量P。

3)模拟材料流动：从左到右各单元依次用右侧相邻单元的内变量P和截面参数C覆盖自身值(图2(a)、图2(b))；滚弯机入口处的单元采用流入坯料的初始参数；将出口处单元的数据写入输出文件(图1)。

4)平衡附加载荷：按式(9)计算单元内附加载荷项，并组装得到整体附加载荷向量；约束三个辊轮处的节点垂向位移，组装式(6)并进行非线性求解，更新变形区的构型和各单元的塑性内变量P(图2(c))。

5)进入下一时间步：返回第3)步继续计算，直至总滚弯长度达到预定目标。

在进行每一时间步的非线性方程的求解时，本文采用了割线迭代法更新kb的方式[29]。求解方案采用Mathematica 10.0编程实现。

2 模拟方案的验证算例

为验证提出的基于欧拉网格的模拟方案，同时作为后续蒙特卡洛模拟的参照，考虑如下算例。

算例1中取对称式三辊滚弯机的底辊距L=300 mm，各辊直径Dt=Db=50 mm，中辊的下压量为d=2.7 mm；坯料初始时完全平直，其横截面为边长10 mm的正方形，材料为理想弹塑性，弹性模量为E=200 GPa ，屈服强度YS=200 MPa ；总滚弯长度为600 mm。

使用本文方案时，底辊间的空间等分为150个单元；单元内用两点高斯积分，每个积分点在厚度方向上等分为20层；每步进料长度∆s=2 mm(恰等于单元长度)。进行非线性求解时，收敛标准取为连续两次迭代的位移差向量的二阶范数，小于中辊下压量的10−7。作为对照，在ANSYS 19.0中建立平面应力模型(图3)，坯料长度和厚度方向的网格尺寸分别为2 mm 和1 mm；设辊轮和坯料间粗糙接触，左底辊主动转动带动坯料进入滚弯机，其余辊轮被动转动；取动态时间步长上限使得坯料每步平移不超过2 mm。此外，还通过曲率积分法得出最终定常状态下的理论解作为参考。

图3 ANSYS中的滚弯平面应力模型Fig.3 Plane stress model of roll bending process in ANSYS

算例在CPU i7-6820HQ，内存32 GB的工作站上以串行方式计算。ANSYS平面应力模型的模拟耗时近10 h，而本文方法仅耗时10 min，计算效率提升明显。

中辊处的支座反力和滚弯机出口处的坯料局部曲率随着滚弯长度变化关系如图4所示。由于平面应力模型的输出曲率是用曲率圆拟合每5个相邻节点的方式求得，这使曲率结果趋于平均，因此在滚弯长度150 mm处的峰值低于本文模型。除此之外本文方法和ANSYS平面模型得到的结果曲线几乎重合。在滚弯开始时，受到直边效应的影响，曲线有较大的波动，而当滚弯长度超过400 mm后，两种模型的结果曲线均不再有明显变化，指示均已基本达到定常状态。这时两种模型的结果均与曲率积分法得到的理论解吻合，误差小于2%。

图4 中辊反力和输出曲率随滚弯长度的变化Fig.4 Top roller reaction forceand output curvature versus roll bending length

提取定常状态下产品横截面上的残余应力分布如图5所示。两种有限元模型得到的残余应力均与理论解吻合良好，其中本文模型的吻合程度更高。这是因为本文模型和理论解均基于欧拉—伯努利梁理论，而ANSYS平面应力模型并不遵循该理论的假设。

图5 定常状态下输出坯料的残余应力在厚度方向上的分布Fig.5 Distribution of residual stress through the thickness of the crosssection in steady state

此外，由图4还可以看出ANSYS模型的辊轮反力和输出曲率曲线始终有微小波动，这主要是由于接触算法的数值误差引起的。这种误差具有一定的随机性，因此可能在蒙特卡洛模拟中与坯料局部曲率导致的结果波动相混淆；但要缩小接触误差将导致接触算法收敛困难。而本文模型能从原理上避免这种误差，输出曲线非常光滑，因此更适于进行蒙特卡洛模拟。

3 坯料有随机局部弯曲时产品的曲率半径分布

机械制造中，包括初始曲率半径在内的坯料误差常被认为服从正态分布。那么运用本文方法进行滚弯的蒙特卡洛模拟时，只需在前述模拟流程的基础上，在每个时间步为入口处的REL生成正态分布的伪随机初始曲率κ0即可。为研究模型参数对产品的曲率半径分布的影响，在算例1的坯料模型参数和网格设置的基础上，选取了4组不同的辊轮位置和REL 长度，并将四种情形分别记为算例2～算例5，如表1所示。对每个算例，考虑了4种偏差水平(σκ0=0.005 m−1∼0.02 m−1)的零均值正态分布随机初始曲率。对每种偏差水平，进行滚弯长度为2 00 m的蒙特卡洛模拟(即样本数量为5万～10万个)，并在统计时去掉初始的1000个样本以避免直边效应的影响。

表1 各数值算例的参数设置Table 1 The parameter sets of the numerical examples

先考虑算例2。当σκ0=0.02 m−1时，随机生成的坯料初始曲率与由此重构得的初始构型如图6(a)所示。可以看到虽然坯料局部初始曲率分布范围较大，但宏观上看200 m 长度的坯料仅有0.3 m 浪高，平均每米浪高仅0.14 mm，在工程上属于直线度很高的坯料，在通常的有限元模拟中可以当作完全平直来处理。而用蒙特卡洛模拟计算输出产品的各单元曲率和重构得的构型如图6(b)所示，各截面的输出曲率虽然很不一致，但产品整体上仍可以用一个圆形很好地拟合。由于该拟合圆的半径就是工程实践中观测到的产品半径，会直接影响到零件的安装和配合，不妨称其为宏观曲率半径；与之对应，由各单元位移场求得的曲率半径称为局部曲率半径；又称当坯料平直时(σκ0=0)滚弯定常状态得到的产品构型为目标形状，其半径为目标曲率半径。

图6(b)中显示考虑了微小的随机初始曲率后，宏观曲率半径略小于目标曲率半径，这并非偶然现象。图7给出了不同偏差水平下的局部曲率半径统计分布。可以看到局部曲率半径基本符合正态分布，且随着初始曲率偏差增大，拟合分布的均值变小，标准差近似线性增长；而宏观曲率半径则逐渐下降，偏离目标曲率半径。

图6 当 σκ0=0.02 m−1时算例2的工件形状Fig.6 Theconfigurations ofthework-piece in Case2when σκ0=0.02m−1

图7 算例2产品的局部曲率半径分布Fig.7 The distributionsof theoutput local radiusin Case 2

采用相同的方法分析了其余算例，其偏差水平在σκ0=0.02 m−1时的输出曲率半径分布如图8所示，可见各算例的输出半径均近似满足正态分布，都有宏观半径小于目标半径，但分布参数不同。

图8 当 σκ0=0.02 m−1时算例3～算例5的局部曲率半径分布Fig.8 The distributions of the output local radius in Case 3～Case 5 when σκ0=0.02 m−1

各算例的统计分布参数随偏差水平的变化如图9所示。从图9(a)看到相比曲率半径的算术平均值，宏观曲率半径的下降趋势更有规律，并且宏观曲率半径更符合工程使用的需要，更适合作为描述结果分布规律的位置参数。图9(b)显示产品的曲率半径标准差与坯料曲率的标准差基本成线性关系。图9(c)显示当坯料曲率标准差由小变大时，分布逐渐由负偏斜变为正偏斜；峰度先减小后增大。在本文考虑的几个算例中，偏度系数和峰度系数都很接近正态分布的对应值(0和3)，因此，输出的曲率半径分布可以用正态分布很好地近似。

图9 产品的局部曲率半径分布参数随坯料曲率标准差的变化规律Fig.9 Distribution parametersof output local radiusversus standard deviation of input curvature

对比算例2和算例3发现，两者的目标曲率半径相近，此时在辊距相差50%的情况下，曲率半径的分布规律和分布参数非常相似。这表明对于给定的目标曲率半径，局部输出曲率的分布与辊轮位置参数几乎无关。

而算例4具有和算例2相同的跨度，但目标曲率半径增大了近30%，结果宏观曲率半径下降幅度增大了近1倍，局部曲率半径的分布更加分散，标准差增大近90%。这表明随着设计目标曲率半径增大，产品受到坯料初始随机曲率的影响逐渐变大。

算例5和算例3只有代表性单元长度REL不同。这里REL 表征随机局部曲率的变化剧烈程度。当坯料曲率的标准差不变时，该长度越大，局部就越均匀，坯料的初始弯曲也就越明显。对比结果发现，当REL 增大时，产品曲率半径的分布更加分散，而宏观曲率半径则没有明显变化。

4 结论

在用有限元模拟研究滚弯过程时，通常假定坯料是平直的或者具有一致的曲率。这种假设是否对模拟结果有偏向性的影响，则未见讨论。为此本文在小曲率、细长梁假设的基础上，采用一种基于欧拉网格的滚弯模拟方案，使得超长长度的滚弯模拟成为可能，在此基础上运用蒙特卡洛方法，对坯料有零均值正态分布的随机局部曲率的滚弯过程进行了研究，结果发现：

(1)产品宏观曲率半径随着初始曲率偏差的增大而减小；这说明要得到理想的产品形状，坯料局部曲率波动较大时，应补偿性地减小下压量。

(2)产品的局部曲率半径大致服从正态分布，随着初始曲率的标准差增大，输出曲率半径分布的标准差线性增加。

(3)目标曲率半径相同时，辊轮位置参数对局部曲率半径的分布没有明显影响。换言之，若不更换坯料，无法通过调整辊轮位置参数来使得输出曲率更加集中到目标曲率附近。

(4)想要得到的产品半径越大，代表性单元长度越长，滚弯结果越容易受到坯料局部初始弯曲的影响。