坯料随机局部弯曲对滚弯成形结果影响的蒙特卡洛分析
摘 要：采用蒙特卡洛有限元模拟方法研究了坯料局部的随机弯曲对滚弯成形结果的影响。为了提高滚弯模拟的效率，提出了基于欧拉网格和经典梁单元的滚弯模拟方案，并与传统有限元模型和理论解对比验证了该方案的正确性。在此基础上模拟了具有零均值正态分布的局部曲率的超长坯料的滚弯过程，并统计产品曲率半径的分布规律。结果表明：输出曲率半径分布近似满足正态分布，且随着坯料曲率标准差的增大，均值减小，方差增大，宏观上表现为产品半径减小。产品的目标半径越大，代表性单元长度越长，受初始弯曲的影响就越大;对于给定的目标形状，辊轮位置参数对实际输出半径的分布没有影响。
郑子君; 陶裕梅， 工程力学 发表时间：2021-08-02
关键词：滚弯;欧拉网格;蒙特卡洛法;欧拉梁;随机误差
滚弯是对型材、板材进行弯曲加工的一种常见工艺。以对称式三辊滚弯机为例，其工作方式为将坯料置于底辊和中辊之间，中辊以给定的下压量压紧坯料，通过辊轮或送料辊转动，带动坯料连续地通过滚弯机并发生塑性弯曲。当加工过程开始一段时间后，滚弯机两底辊之间的坯料将逐渐达到一个“定常”状态，此时辊轮上的反力与流出出口的产品曲率不再随着时间变化。
由于滚弯工艺属于无模成型，回弹量大，实践中要得到目标曲率往往需要对工艺参数进行大量的试错[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 所示坐标系。在欧拉观点中，所有变量是一维空间坐标 的函数。记 t 时刻 处的曲率为 ，将 处横截面的几何尺寸、材料属性、初始曲率等截面参数合并记为向量 ，塑性内变量集合记为向量 。根据欧拉梁的本构关系，弯矩可以由变形程度、截面参数以及塑性内变量确定，那么经过 的时间，弯矩增量为：
设在 时刻，位于 处的材料截面由 t 时刻位于 处的材料截面被辊轮带动平移而来 (图 2(a)、 图 2(b))， 于 是 有 。又当 很短时，可认为任意材料点的应力路径是简单的，于是计算 时刻 处的弯矩时只需知道 t 时刻该材料截面的塑性内变量 ，于是式 (1) 变为：
在小变形情形下， ∆s 是与坐标 x 无关的常数，其物理意义为该时间步内滚弯机的进料长度。式 (2) 进一步改写为：
显然，式 (3) 中第一项是由于空间 处的坯料构型改变引起的，形式上可以写为 ，其中为截面参数为 ，塑性内变量为时的割线弯曲刚度。实际计算时， 需要用迭代法确定。
而式 (3) 第二项是 处构型不变，只是材料发生流动时引起的，类似于流体力学中的迁移项。该项涉及的所有变量在 时刻已经求得，因此是已知的。记该项为 ，物理意义为 x 处的构型不变而截面参数发生变化时的内力不平衡量。则式 (3) 可写为： kb (x)∆κ (x)=∆Mex (x)+ ∆Ma (x) (4) 与材料坐标下的经典梁方程比较： kb (s)∆κ (s)=∆Mex (s) (5)
可以发现形式上仅仅相差了一项 ，因此对式 (4) 进行有限元离散时，只需要在拉格朗日列式的基础上增加与内力不平衡量 对应的载荷项即可。如果采用经典的 2 节点 4 变量 Hermite 插值，则式 (4) 的离散形式可写为： Ke∆u=∆Fex + ∆Fa (6) ∆u = {∆wI ,∆θI ,∆wJ ,∆θJ } Ke ∆Mex(x) ∆Fex 其中：单元位移向量 由节点处的挠度和转角的增量组成;单元割线刚度矩阵 和与 对应的单元载荷向量 的公式与经典梁单元中的推导完全相同，即[29]：
式中： 为单元局部坐标; 为形函数向量; 为单元长度; 、 、 分别是分布力集度、作用在 处的集中力和作用在 处的集中力偶的增量。而与 对应的单元附加载荷向量 为：
至此已得到式 (4) 的单元离散方程。在滚弯模拟中的每一时间步，总刚矩阵和载荷向量的组装、约束的施加与非线性方程迭代求解过程均与常规的有限元方法基本相同，不同之处是该步开始前，需更新单元参数，即令 ，以模拟材料在滚弯机内的水平流动。每一步模拟的原理性示意图如图 2 所示。
为了实现简单，实际计算时可取单元长度 恰等于每步进料长度 ，并忽略辊轮的尺寸。此时具体的模拟步骤如下：
1) 初始化：将两个底辊间的变形区空间等分成若干单元，初始化各单元高斯积分点处的塑性内变量 和截面参数 。
2) 模拟调辊过程：约束底辊处的两节点的垂向位移，为中辊下方节点施加指定的位移量，取附加载荷项为零，组装式 (6) 并进行非线性求解，更新塑性内变量 。
3) 模拟材料流动：从左到右各单元依次用右侧相邻单元的内变量 和截面参数 覆盖自身值 (图 2(a)、图 2(b));滚弯机入口处的单元采用流入坯料的初始参数;将出口处单元的数据写入输出文件 (图 1)。
4) 平衡附加载荷：按式 (9) 计算单元内附加载荷项，并组装得到整体附加载荷向量;约束三个辊轮处的节点垂向位移，组装式 (6) 并进行非线性求解，更新变形区的构型和各单元的塑性内变量 (图 2(c))。
5) 进入下一时间步：返回第 3) 步继续计算，直至总滚弯长度达到预定目标。
在进行每一时间步的非线性方程的求解时，本文采用了割线迭代法更新 的方式[29]。求解方案采用 Mathematica 10.0 编程实现。
2 模拟方案的验证算例
为验证提出的基于欧拉网格的模拟方案，同时作为后续蒙特卡洛模拟的参照，考虑如下算例。
算例 1 中取对称式三辊滚弯机的底辊距，各辊直径 ，中辊的下压量为 ;坯料初始时完全平直，其横截面为边长 的正方形，材料为理想弹塑性 ， 弹 性 模 量 为 ， 屈 服 强 度;总滚弯长度为 。
使 用 本 文 方 案 时 ， 底 辊 间 的 空 间 等 分 为 150 个单元;单元内用两点高斯积分，每个积分点在厚度方向上等分为 20 层;每步进料长度 mm (恰等于单元长度)。进行非线性求解时，收敛标准取为连续两次迭代的位移差向量的二阶范数，小于中辊下压量的 10−7。作为对照，在 ANSYS 19.0 中建立平面应力模型 (图 3)，坯料长度和厚度方向的网格尺寸分别为 和 ;设辊轮和坯料间粗糙接触，左底辊主动转动带动坯料进入滚弯机，其余辊轮被动转动;取动态时间步长上限使得坯料每步平移不超过 。此外，还通过曲率积分法得出最终定常状态下的理论解作为参考。
算例在 CPU i7-6820HQ，内存 32 GB 的工作站上以串行方式计算。ANSYS 平面应力模型的模拟耗时近 10 h，而本文方法仅耗时 10 min，计算效率提升明显。
中辊处的支座反力和滚弯机出口处的坯料局部曲率随着滚弯长度变化关系如图 4 所示。由于平面应力模型的输出曲率是用曲率圆拟合每 5 个相邻节点的方式求得，这使曲率结果趋于平均，因此在滚弯长度 处的峰值低于本文模型。除此之外本文方法和 ANSYS 平面模型得到的结果曲线几乎重合。在滚弯开始时，受到直边效应的影响，曲线有较大的波动，而当滚弯长度超过后，两种模型的结果曲线均不再有明显变化，指示均已基本达到定常状态。这时两种模型的结果均与曲率积分法得到的理论解吻合，误差小于 2%。
提取定常状态下产品横截面上的残余应力分布如图 5 所示。两种有限元模型得到的残余应力均与理论解吻合良好，其中本文模型的吻合程度更高。这是因为本文模型和理论解均基于欧拉 —伯努利梁理论，而 ANSYS 平面应力模型并不遵循该理论的假设。
此外，由图 4 还可以看出 ANSYS 模型的辊轮反力和输出曲率曲线始终有微小波动，这主要是由于接触算法的数值误差引起的。这种误差具有一定的随机性，因此可能在蒙特卡洛模拟中与坯料局部曲率导致的结果波动相混淆;但要缩小接触误差将导致接触算法收敛困难。而本文模型能从原理上避免这种误差，输出曲线非常光滑，因此更适于进行蒙特卡洛模拟。
3 坯料有随机局部弯曲时产品的曲率半径分布
机械制造中，包括初始曲率半径在内的坯料误差常被认为服从正态分布。那么运用本文方法进行滚弯的蒙特卡洛模拟时，只需在前述模拟流程的基础上，在每个时间步为入口处的 REL 生成正态分布的伪随机初始曲率 即可。为研究模型参数对产品的曲率半径分布的影响，在算例 1 的坯料模型参数和网格设置的基础上，选取了 4 组不同的辊轮位置和 REL 长度，并将四种情形分别记为算例 2~算例 5，如表 1 所示。对每个算例，考虑了 4 种偏差水平 ( ) 的零均值正态分布随机初始曲率。对每种偏差水平，进行滚弯长度为 的蒙特卡洛模拟 (即样本数量为 5 万~10 万个)，并在统计时去掉初始的 1000 个样本以避免直边效应的影响。
先考虑算例 2。当 时，随机生成的坯料初始曲率与由此重构得的初始构型如图 6 (a) 所示。可以看到虽然坯料局部初始曲率分布范围较大，但宏观上看 200 m 长度的坯料仅有 0.3 m 浪高，平均每米浪高仅 0.14 mm，在工程上属于直线度很高的坯料，在通常的有限元模拟中可以当作完全平直来处理。而用蒙特卡洛模拟计算输出产品的各单元曲率和重构得的构型如图 6 (b) 所示，各截面的输出曲率虽然很不一致，但产品整体上仍可以用一个圆形很好地拟合。由于该拟合圆的半径就是工程实践中观测到的产品半径，会直接影响到零件的安装和配合，不妨称其为宏观曲率半径;与之对应，由各单元位移场求得的曲率半径称为局部曲率半径;又称当坯料平直时 ( ) 滚弯定常状态得到的产品构型为目标形状，其半径为目标曲率半径。
图 6 (b) 中显示考虑了微小的随机初始曲率后，宏观曲率半径略小于目标曲率半径，这并非偶然现象。图 7 给出了不同偏差水平下的局部曲率半径统计分布。可以看到局部曲率半径基本符合正态分布，且随着初始曲率偏差增大，拟合分布的均值变小，标准差近似线性增长;而宏观曲率半径则逐渐下降，偏离目标曲率半径。
采用相同的方法分析了其余算例，其偏差水平在 时的输出曲率半径分布如图 8 所示，可见各算例的输出半径均近似满足正态分布，都有宏观半径小于目标半径，但分布参数不同。
各算例的统计分布参数随偏差水平的变化如图 9 所示。从图 9 (a) 看到相比曲率半径的算术平均值，宏观曲率半径的下降趋势更有规律，并且宏观曲率半径更符合工程使用的需要，更适合作为描述结果分布规律的位置参数。图 9 (b) 显示产品的曲率半径标准差与坯料曲率的标准差基本成线性关系。图 9 (c) 显示当坯料曲率标准差由小变大时，分布逐渐由负偏斜变为正偏斜;峰度先减小后增大。在本文考虑的几个算例中，偏度系数和峰度系数都很接近正态分布的对应值 (0 和 3)，因此，输出的曲率半径分布可以用正态分布很好地近似。
对比算例 2 和算例 3 发现，两者的目标曲率半径相近，此时在辊距相差 50% 的情况下，曲率半径的分布规律和分布参数非常相似。这表明对于给定的目标曲率半径，局部输出曲率的分布与辊轮位置参数几乎无关。
而算例 4 具有和算例 2 相同的跨度，但目标曲率半径增大了近 30%，结果宏观曲率半径下降幅度增大了近 1 倍，局部曲率半径的分布更加分散，标准差增大近 90%。这表明随着设计目标曲率半径增大，产品受到坯料初始随机曲率的影响逐渐变大。
算例 5 和算例 3 只有代表性单元长度 REL 不同。这里 REL 表征随机局部曲率的变化剧烈程度。当坯料曲率的标准差不变时，该长度越大，局部就越均匀，坯料的初始弯曲也就越明显。对比结果发现，当 REL 增大时，产品曲率半径的分布更加分散，而宏观曲率半径则没有明显变化。
4 结论
在用有限元模拟研究滚弯过程时，通常假定坯料是平直的或者具有一致的曲率。这种假设是否对模拟结果有偏向性的影响，则未见讨论。为此本文在小曲率、细长梁假设的基础上，采用一种基于欧拉网格的滚弯模拟方案，使得超长长度的滚弯模拟成为可能，在此基础上运用蒙特卡洛方法，对坯料有零均值正态分布的随机局部曲率的滚弯过程进行了研究，结果发现：
(1) 产品宏观曲率半径随着初始曲率偏差的增大而减小;这说明要得到理想的产品形状，坯料局部曲率波动较大时，应补偿性地减小下压量。
(2) 产品的局部曲率半径大致服从正态分布，随着初始曲率的标准差增大，输出曲率半径分布的标准差线性增加。
(3) 目标曲率半径相同时，辊轮位置参数对局部曲率半径的分布没有明显影响。换言之，若不更换坯料，无法通过调整辊轮位置参数来使得输出曲率更加集中到目标曲率附近。
(4) 想要得到的产品半径越大，代表性单元长度越长，滚弯结果越容易受到坯料局部初始弯曲的影响。