一种提高管材弯曲成型过程仿真效率的方法及仿真方法

摘要

本发明属于机械加工技术领域，公开了一种提高管材弯曲成型过程仿真效率的方法及仿真方法，该方法对用于描述弯曲成型过程的运动方程的时间域进行当量变换，将实际的弯曲时间域[0, T]变换为标准的时间域[0，1]，空间域保持不变；同时用时间域变换的变换系数对管材密度、初始速度和初始角速度进行变换，获得标准时间域下的管材密度、初始速度与初始角速度；再基于变换后获得的等价运动方程和变换后管材密度、初始速度和初始角速度，对管材弯曲成型过程进行仿真求解。本发明方法可缩短仿真计算的时间，可以适用于管材弯曲成型模具与工艺的设计中，具有很好的工程应用前景。

说明书

技术领域

本发明属于机械加工技术领域，涉及弯管成型过程仿真技术，关于一种关于管材成 型仿真的计算方法，特别指一种弯管绕弯成型过程仿真分析与提高仿真效率的方法。

背景技术

以管材为毛坯，通过塑性加工手段，制造管材零件的工艺称为管材塑性加工。管材 塑性加工，由于易于实现产品的轻量化、强韧化和满足低耗高效、精确制造等方面的要求， 并且能获得形状复杂的制品，具有精确稳定、无废料绿色环保等优点。管材塑性成型得到的 弯管零件具有轻量化、结构稳定、高刚度等优良的力学性能，而成为塑性加工技术发展的一 个重要方向。管材弯曲成型技术是管材塑性加工的重要组成部分。用管材制造的弯曲件，具 有重量轻、吸震力强、介质流通量大等一系列优良性能，广泛应用于普通气体、液体的输送 管路及汽车、航空航天、化工和其它高新技术领域。

在管材的弯曲成型过程中，容易出现回弹、起皱、截面畸变及减薄拉裂等问题。弯 曲成型过程中，管材的外侧壁受拉而减薄甚至拉裂，内侧壁受压而增厚起皱，内、外侧受力 不匀导致截面畸变，卸载后又存在回弹现象。这些问题一直是管材弯曲成型加工行业未能有 效解决的技术难题，也是当今国内外塑性加工学科研究的难点和热点。采用计算机辅助工程 分析，即CAE技术，对管材弯曲成型过程进行仿真分析可以有效预测弯曲成型过程中出现 的各种缺陷，是目前流行的设计手段与分析工具。

对管材弯曲成型问题采用计算机辅助工程分析的CAE技术研究弯曲模具结构及弯曲 成型工艺，涉及到管材的几何大变形、材料的弹塑性变形及管材与弯曲模具的接触等问题。 这是一个高度非线性的仿真模拟过程，成型过程仿真分析的计算量非常大。

为缩短仿真时间、提高分析效率，工程上采用一些特殊的方法来加速模拟过程，如 在板料冲压成形的模拟中，常用速度加速法，通过人为提高冲头的运动速度，即将几十毫米 /秒的实际冲压速度提高到3000mm/s左右(称为虚拟冲压速度)，来缩短弯曲过程仿真的时 间。又如质量缩放法，即将某些特定单元的质量密度放大来加大时间步长，但是缩放系数的 确定需要试探。这种单方面调整某一物理量的方法尽管提高了计算效率但容易产生误差，且 不够严谨、缺乏理论基础的支持。

发明内容

本发明目的是提出一种管材绕弯成型过程仿真分析与提高仿真效率的计算方法。该 方法可以对管材弯曲成型过程进行仿真计算，并且缩短了仿真计算的时间，因此可以有效应 用于管材弯曲成型模具与工艺的设计中，具有很好的工程应用前景。

基于上述目的，本发明首先提出一种提高管材弯曲成型过程仿真的效率的方法，该 方法对用于描述弯曲成型过程的运动方程的时间域进行当量变换，将实际的弯曲时间域[0,T] 变换为标准的时间域[0，1]，空间域保持不变；同时用时间域变换的变换系数对管材密度、 初始速度和初始角速度进行变换，获得标准时间域下的管材密度、初始速度与初始角速度； 再基于变换后获得的等价运动方程和变换后管材密度、初始速度和初始角速度，对管材弯曲 成型过程进行仿真求解。

具体地说，上述方法采用变换方程t₁＝αt对运动方程进行变换，其中取α＝1/T。管 材密度的变换公式为ρ₁＝ρ/T²，初始速度变换公式为v₁＝v₀T，初始角速度变换公式为 ω₁＝ω₀T＝v₀T/R。

作为改进，为了使求解递推公式的临界时间步长时间域变换后与变换前保持不变， 在仿真求解前，还需要选择管材或弯曲模具的某一零件，如芯棒，对其质量密度进行密度缩 放，密度缩放公式为ρ_i＝ρ₁/α²＝ρ₁T²。

基于上述方法，本发明进一步提出一种基于当量变换的管材弯曲成型过程仿真方法， 包括如下步骤：

步骤一，将用于描述管材弯曲成型过程的运动方程的时间域进行当量变换，将实际的弯曲时 间域[0,T]变换为标准的时间域[0，1]，空间域保持不变；

步骤二，用步骤1时间域变换中的变换系数对管材密度、初始速度、初始角速度进行等价变 换，获得标准时间域下的管材密度、初始速度和初始角速度；

步骤三，选择管材或弯曲模具的某一零件，如芯棒，对其质量密度进行缩放，使求解递推公 式的临界时间步长与时间域变换前保持不变；

步骤四，基于步骤1-3变换后的标准时间域下的等价运动方程、管材密度、初始速度、初始 角速度，对管材弯曲成型过程进行仿真求解。

上述方法中，步骤1采用变换方程t₁＝αt对运动方程进行变换，其中取α＝1/T。

步骤2中管材密度的变换公式为ρ₁＝ρ/T²，初始速度变换公式为v₁＝v₀T，初始角速 度变换公式为ω₁＝ω₀T＝v₀T/R。

步骤3中采用密度缩放公式ρ_i＝ρ₁/α²＝ρ₁T²对所选择的管材或弯曲模具的某一零件 的密度进行缩放。

下面分步对本发明方法做出详细的解释说明。

步骤1：弯曲成型基本方程变换。

管材弯曲成型问题属于结构动力学问题。结构动力学问题的控制方程可以通过运动

方程描述，它是借助牛顿(Newton)第二定律这一基本力学原理建立的数学表达式，描述

了结构中力与位移(包括速度和加速度)的关系。结构动力学问题的运动方程一般形式：

σ_ij,j-ρu_i,tt-μu_i,t+f_i＝0 t∈[0,T],i＝x,y,z∈Ω    (1)

其中，σ_ij是应力张量，u_i是位移矢量，ρ是材料密度，μ是阻尼系数，T是弯曲成型所需的 时间，u_i,t是位移矢量对时间的一阶导数即速度，u_i,tt是位移矢量对时间的二阶导数即加速度， μu_i,t是粘性力，ρu_i,tt是惯性力。

方程(1)边界条件中的初始条件为：

t＝0,u＝u₀,v＝v₀    (2)

其中，u₀是初始位移，v₀是初始速度。

方程(1)的求解一般使用数值方法，如采用有限元方法对方程(1)进行离散求解。 由于方程(1)包含空间坐标变量与时间坐标变量，因此有限元离散求解的做法是：空间域 网格划分建立有限元模型，时间域采用差分格式建立递推公式。

本发明的做法是对公式(1)所描述的弯曲成型基本方程的时间域进行变换，即时间 变量的变换表达式如下：

t₁＝αt    (3)

其中，α为变换因子，又称缩放系数。这样将实际弯曲成型过程的时间域[0,T]变换到了当量 仿真时间域[0，T₁]＝[0，αT]。很显然，若取α＝1/T，则有T₁＝1，即将实际弯曲成型时间T 变换到当量仿真时间1。α也可取其它小于1的值，这样可变换到更小的时间域。称当量仿 真时间域[0，1]为标准时间域。

时间域变换的同时对弯曲问题(1)中的空间变量(x,y,z)不作变换，这样(1)式 等价为

${\sigma }_{\mathrm{ij},j}-{\alpha }^2\rho u_{i,t_1{t}_1}-\mathrm{\alpha \mu }{u}_{i,t_1}+f_i=0,t_1\in [0,T_1],i=x,y,z\in \Omega ---\left(4\right)$

初始条件(2)也变换为：

t₁＝0,u₁＝u₀,v₁＝v₀/α    (5)

这样通过(3)式的时间域变换将(1)式所描述的管材弯曲成型问题等价变换为(4)描述 的不同时间域的管材弯曲成型问题。很显然问题(1)与问题(4)的管材及模具的尺寸与拓 扑关系一致，因此有限元网格模型也不变，由于变换而发生改变的是管材及模具的材料参数 与边界条件中的初始条件。这样的变换方便工程实际应用。

步骤2：参数计算。

经过第一步的时间变换后，用于管材弯曲仿真的主要参数，即材料密度与初始速度

可以从原始数据变换计算得到(取α＝1/T)：

1)材料密度：ρ₁＝ρ/T²    (6)

2)初始速度：v₁＝v₀T    (7)

使用本方法时须注意各个物理量的量纲，即量纲要统一。如采用国际单位制SI，基本量纲为： 毫米mm，千克kg，秒s，则力的量纲为N(牛顿)、材料密度量纲为kg/mm³，速度的量纲 为mm/s，弹性模量E与应力σ的量纲为MPa(N/mm²)。很显然，变换因子α是无量纲的， 为一放大缩小的缩放因子。

管材弯曲工艺中有时还包含弯曲模的转动角速度ω，可以按公式ω＝v/R(R弯曲模 半径)通过弯曲速度进行计算，即初始角速度的变换计算为：

3)初始角速度：ω₁＝ω₀T＝v₀T/R    (8)

步骤3：材料密度缩放。

根据结构动力学关于控制方程(1)数值求解的分析，为保证数值求解的收敛性对时 间域上求解递推公式的时间步长取值有个临界值，即临界时间步长。这个临界时间步长与材 料参数、网格尺寸有关，即由下式计算：

$\Delta t_{\min }=l_{\min }/\sqrt{E/(1-{\mu }^2)\rho }---\left(9\right)$

其中：l_min是单元最小尺寸，ρ是材料密度，E是材料的弹性模量，μ是泊松比。工程上较多 使用的有限元分析软件进行冲压、弯曲成型过程的仿真计算时都采用(9)式进行时间步长 的自动计算。

由公式(9)可以推出，仿真时间域变小后临界时间步长也变小，而且缩放的系数也 与时间变换的系数相同，即都是α倍。为了提高计算效率，应使求解所需的计算时间步长在 时间变换的前后不变。这样变换后时间域变小了，而步长不变，则总的求解步数减少了。而 公式(9)也表明在对时间坐标变换之后，可以对有限元网格模型中单元尺寸最小的单元 (即对应l_min的单元)的质量密度进行缩放来放大时间步长，从而保证计算时间步长不变。

因此，本算法在时间坐标变换之后再选择性的对管材及弯曲模具的材料密度进行缩 放，目的是使材料密度放大，从而放大时间步长。如何选择缩放的系数很关键，以往的“质 量缩放法”是单独对模具或管材的质量缩放，即对其材料密度缩放，但采用的做法是试凑， 要达到满意的效果很困难。本算法从公式(9)推导求解得到了最佳缩放的系数，即要求时 间坐标变换前后的临界时间步长不变。具体算法是将单元尺寸最小单元的质量密度放大，即：

ρ_i＝ρ₁/α²＝ρ₁T²    (10)

其中，i为需要密度缩放的单元号，即单元尺寸最小的单元；ρ₁是时间变换后的材料密度。

上式中的i也可以是需要密度缩放的零件号。这样有利于实际工程应用。因为在使用 仿真软件分析时，网格划分是自动完成的，而且网格数目很大，靠手工寻找单元尺寸最小的 单元很困难。但可以人为设定需要密度缩放的零件，并控制其网格划分的疏密，从而保证满 足单元尺寸最小的要求。

本算法提出了选择密度缩放因子的方法，即在时间变换的基础上选择时间变换因子α 的平方作为密度缩放因子，这样一方面使时间域缩小，另一方面保证了时间步长不变，使得 仿真时间缩短，达到提高仿真效率的效果。这样克服了以往的“质量缩放法”的缺点，即质 量缩放无规律可循，因而质量缩放很容易造成计算失败。

常用的管材弯曲成型仿真模型中有一种简化刚体模型是将弯曲模具作为刚体对管材 施加约束。由于只有弹性体可以进行密度缩放，因此简化刚体模型中只有管材和芯棒连接体 是弹性体，是可以进行密度缩放的。实际计算时可选择一个弹性体进行密度缩放，如选芯棒 连接体或选管材，按公式(10)密度缩放。

步骤四：仿真求解。

经过前面的变换后，仿真求解过程则与常规则的仿真求解一样了。同样是使用有限 元仿真分析软件进行仿真求解，常见的有限元仿真分析软件有Ls-dyna及Abaqus。

附图说明

图1管材弯曲成型过程仿真提高计算效率的当量法流程图。

图2为本发明实施案例1的弯管截面厚度变化的对比，其中，图2a为现有方法的结 果，图2b为本发明方法的结果。

图3为本发明实施案例1的弯管截面椭圆度的对比，其中，图3a为现有方法的结果， 图3b为本发明方法的结果。

图4为本发明实施案例2的304不锈钢薄壁圆管芯棒工具设置。

图5为本发明实施案例2的弯管壁厚变化情况计算对比，其中，图5a为现有方法的 结果，图5b为本发明方法第一种方案的结果，图5c为本发明方法第二种方案的结果。

图6为本发明实施案例2的弯管椭圆度计算对比，其中，图6a为现有方法的结果， 图6b为本发明方法第一种方案的结果，图6c为本发明方法第二种方案的结果。

具体实施方式

下面结合几个实例进一步说明本发明所述的管材弯曲成型过程仿真效率提高的当量 法计算方法的有效性。

例1，弯曲问题描述如下：

圆管材料为Q235、密度为、圆管半径为16.85mm、壁厚为3mm。弯管的弯曲角度为90度、 弯曲半径为130mm。弯曲工艺为：弯曲速度，弯曲时间T＝8秒，即时间域[0，T]＝[0，8]。

在有限元仿真分析软件平台，如Ansys LS-dyna上进行仿真求解，根据上述数据采用 常规方法进行弯曲成型过程仿真，其弯曲成型过程仿真的计算时间是623秒，得到的计算结 果称为原始数据计算结果。

以下是采用本发明的当量法对本例描述问题进行计算步骤。

步骤1：弯曲成型基本方程变换。

根据弯曲工艺的弯曲时间T＝8s，取缩放系数α＝1/T＝1/8，则时间变换公式为 t₁＝t/8，对方程(1)进行时间域的时间变量变换，将方程(1)表达的弯曲成型基本方程 转换为等价的弯曲成型基本方程，变换后的等价运动方程为：

${\sigma }_{\mathrm{ij},j}-\rho u_{i,t_1{t}_1}/64-\mu u_{i,t_1}/8+f_i=0,t_1\in [0,1],i=x,y,z\in \Omega$

其中当量仿真时间域为[0，1]。

步骤2：参数计算。

从原始数据变换计算获得材料密度ρ₁与初始速度v₁，依据公式(6,7)有：

ρ₁＝ρ/T²＝7.86×10^-6/8²＝0.12266×10^-6kg/mm³，

v₁＝v₀T＝26.939×8＝215.51mm/s。

步骤3：材料密度缩放。

根据该管材的璧厚推论无须使用芯棒工具即可满意弯曲成型，因此本例只有管材弹 性体可以进行密度缩放。依据公式(10)，对管材弹性体进行密度缩放，有： ρ_i＝ρ₁T²＝0.1226610^-6×64＝7.86×10^-6kg/mm³。

步骤4：仿真求解。

在Ansys LS-dyna平台上进行仿真分析，仿真计算时间为94秒。

弯管弯曲成型的成型质量指标主要是考虑弯管截面厚度变化的减薄率、增厚率及弯 管截面形状变化的畸变程度，即截面的椭圆度。将弯曲成型后弯管截面厚度减小的最大值与 弯曲前的厚度比定义为减薄率，将弯曲成型后弯管截面厚度增加的最大值与弯曲前的厚度比 定义为增厚率，将成型后弯管截面形状变化引起的最大直径减最小直径与弯曲前的管材直径 比定义为椭圆度。减薄率、增厚率及椭圆度的计算公式如下：

$\xi =\frac{(t-t_{\min })}{t}---\left(11\right)$

$\xi =\frac{(t_{\max }-t)}{t}---\left(12\right)$

式中：ξ—壁厚减薄率；t—管材原始壁厚；t_min—为管材弯曲后最小壁厚，t_max—为管材弯曲 后最小壁厚。—椭圆度；d_max—弯曲后管材同一截面的最大外径尺寸；d_min—弯曲后管材同 一截面的最小外径尺寸；d—弯曲前管材同一截面的外径尺寸。

对常规分析方法与当量法二种方法的仿真结果从减薄率、增厚率与椭圆度的指标进 行对比，即根据公式(11～13)进行计算，表1给出了相关计算数据。对比常规分析方法与 当量法可见，减薄率的差距为15％，增厚率的差距为0.6％，而椭圆度二者的差距为0.4％。 可见误差在可接受范围，但仿真计算时间大大缩短，只有原始数据仿真计算时间的15％。图 2给出了二种算法弯曲成型后弯管截面厚度变化的对比，图3给出了二种算法的弯管截面变 化的对比。

表1：例1相关计算结果

  减薄率/％ 增厚率/％ 椭圆度/％ 计算时间/s 原始计算 3.93 10.55 0.779 623 当量法 4.52 10.61 0.776 94

例2，弯曲问题描述如下：

这是材料为304不锈钢、外径为50mm、壁厚为0.6mm的薄壁圆管的弯曲，管材的密度为 ρ＝7.9×10^-6kg/mm³。弯管的弯曲角度90°，弯曲半径150mm。弯曲工艺为：弯曲速度0.4 弧度/秒，弯曲时间T＝3.93秒。

由于管材壁厚太薄，弯曲过程容易产生塌陷起皱的缺陷，因此使用芯棒工具，即在 管材内放入由若干芯头组成的芯棒工具。芯棒工具随管材一起运动弯曲，在弯曲过程中对管 材起到支撑防止管壁塌陷起皱的作用。见图5，使用了7个芯头才保证了管材顺利成型。这 样增加了管材与芯头的接触对数，使弯曲仿真的时间大大增加，采用常规方法仿真计算的时 间是17739秒，约5小时。

以下是采用本发明的当量法对本例描述问题进行计算步骤。

步骤1：弯曲成型基本方程变换。

根据弯曲工艺的弯曲时间T＝3.93s，取缩放系数α＝1/T＝1/3.93，则时间变换公式为 t₁＝t/3.93，对方程(1)进行时间域的时间变量变换，将方程(1)表达的弯曲成型基本方 程转换为等价的弯曲成型基本方程，其等价运动方程为：

${\sigma }_{\mathrm{ij},j}-\rho u_{i,t_1{t}_1}/15.44-\mu u_{i,t_1}/3.93+f_i=0,t_1\in [0,1],i=x,y,z\in \Omega$

其中当量仿真时间域为[0，1]。

步骤2：参数计算。

从原始数据变换计算获得材料密度ρ₁与初始速度v₁，依据公式(6,8)，有：

ρ₁＝ρ/T²＝7.9×10^-6/3.93²＝0.5115×10^-6kg/mm³，

ω₁＝ω₀T＝0.4×3.93＝1.572弧度/s。

步骤3：材料密度缩放。

在本例中，对芯棒结构采用超弹体连接模型进行分析建模，即各芯头的连接件用超 弹体，芯头与连接件之间是摩擦接触关系。这样的分析模型有两个对象可以进行质量缩放， 一个是管材，一个是连接件。因此分别对这二个对象零件进行质量缩放比较效果。

方案1，对芯棒采用密度缩放，有：

ρ_i＝ρ₁T²＝0.5115×10^-6×3.93²＝7.9×10^-6kg/mm³。

方案2，对管材密度缩放，即管材缩放后的密度为：

ρ_i＝ρ₁T²＝0.5115×10^-6×3.93²＝7.9×10^-6kg/mm³。

步骤4：仿真计算。

方案1，仿真的计算时间为8099秒(约2小时15分)。

方案2，计算时间为5048秒(约1小时24分)。

图5给出了弯管截面厚度变化的对比，图6给出了二种算法的弯管截面变化的对比。 根据公式(11～13)进行计算，方案1的减薄率最大差距为1.4％，方案2的减薄率最大差 距为0.9％。表2给出了相关数据。而仿真计算时间都大大缩短，尤其当量方案法2，计算效 率提高是非常明显的。

表2：例2相关计算结果

  减薄率/％ 增厚率/％ 椭圆度/％ 计算时间/s 原始计算 11.49 10.67 0.31 17739 方案1 11.33 10.52 0.17 8099 方案2 11.39 10.81 0.47 5048

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范 围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员 根据本发明构思所能够想到的等同技术手段。