空间弯管仿真与回弹补偿

熊威甘忠

西北工业大学现代设计与集成制造技术教育部重点实验室，西安，710072

0 引言

作为一种重要的承载和输送气体、液体的轻质结构件，金属弯管零件在高科技领域(如航空航天、造船、汽车、能源和医疗)得到了越来越多的关注［1-2］。为了将管材弯成具有一定弯曲半径、角度和形状的零件，弯管工艺作为管材塑性成形的重要分支，是生产轻质结构的一种十分重要的制造技术［1，3］。

以前的弯管研究主要集中在平面弯管领域，而对空间弯管(在弯曲段中，管材轴线的挠率不恒等于零的管弯零件)较少涉及。近年来，土木工程和交通系统中对三维弯曲钢铝构件的需求日趋强烈，空间弯管因其在空间和流体力学上的优势，为设计者提供了新的选择［4］。例如美国设计生产的太阳能STERLING发电机中散热管系中的管型，管材轴线即为空间曲线［5］。文献［6］提出了一种自由弯曲方法(free－bending)，该方法将管材挤过一个可以在两个垂直方向偏转的弯曲模腔，几乎能够成形任意外形的弯管，但是弯曲半径要是管材半径的2.5倍，并且对材料参数的变化较为敏感。另一种新的空间弯管方法［4，7］是扭矩叠加空间弯曲工艺(torque superposed spatial bending)，它基于滚弯工艺，通过动态调整弯曲外形和施加扭矩，来达到弯曲不同半径和角度的目的。

一般而言，弯管成形的主要缺陷有失稳起皱、壁厚过度减薄、截面畸变、拉裂和回弹［8-9］。目前，对空间弯管的研究主要集中在对小直径管材的弯曲研究上［5-6，10］。相对于管材的破裂起皱和截面畸变，管材的回弹是主要的成形缺陷［10］。在平面弯管的研究中，如果使用理想弹塑性模型，并且假设在弯曲过程中管径保持不变，回弹可以使用弯曲梁理论进行计算［11-12］。但因为忽略了三维应力应变，故纯解析计算结果会和试验结果相差较大［1］，所以有限元分析成为弯管回弹的主要方法。在空间弯管的回弹研究方面，文献［5］将管材的空间轴线离散为若干段圆弧，并通过纯弯曲试验确定每段圆弧在回弹前的半径，最终在UG软件中拼接各段圆弧，得到绕弯模具型面。实质上是将空间弯曲用若干个密切平面上的平面弯曲近似处理。文献［10］发展了上述平面化的补偿方法，将空间轴线分段投影到密切平面和从切平面上，同样使用纯弯曲试验确定每段圆弧在回弹后的半径，再将对应的一对投影圆弧在空间中合成一段空间曲线，同样使用UG拼合，预测了空间弯管回弹后的形状。上述文献中没有论述平面化回弹预测和补偿的数学依据，并且没有涉及空间绕弯的数值仿真问题。

本文通过有限元仿真研究小直径空间弯管的回弹问题。提出了一种具有数学依据的空间弯管回弹补偿方法。

1 空间弯管数值模拟

1.1 有限元模型

空间弯管使用绕弯的成形方法，绕弯模具的几何模型(见图1)，是沿初始截面和终止截面之间的脊线放样而成的。两个截面分别位于XOY面和YOZ面上;脊线是一段空间曲线，起点为rs(0)=(x0，H －rm，0)，终点的Y坐标比起点的Y坐标小H，曲线的参数方程为

图1 模具的几何模型

空间弯管的有限元模型见图2，管材使用C3D8R网格，管材外半径rp=2.8mm，略小于型面半径rm=3.0mm，这是因为将模具型面的截面半径同管材半径设为一致，在进行网格划分时，两者的网格有可能相互干涉而导致分析失败。管材的内径为5.0mm，材料为钢材，使用线性强化模型，设定弹性模量为210GPa，泊松比为0.33，屈服应力σy=418MPa，硬化系数Hp=31 579GPa。模具和压板的弹性模量比管材的弹性模量大一个数量级，以保证在成形过程中不发生过大的变形。

压板是内径为6mm、厚1mm、旋转角为90°的环扇形，参考点设在圆心。压板比管材略长，这是出于有限元分析收敛性的考虑，长的压板可以控制自由端，防止产生不确定的运动，使成形过程为准静态。

图2 空间弯管有限元模型

模型使用ABAQUS/Standard求解，只设定一个分析步，时间为2s，在前一秒中压板参考点沿轨迹从XOY面运动到YOZ面，对管材加载;后一秒，压板参考点沿原路径返回XOY面，完成卸载。

1.2 用户幅值子程序

空间弯管使用绕弯的方式成形时，需要控制压板沿空间曲线运动，在弯管的过程中，压板的6个自由度需要相互协调运动，都是时间的函数。在ABAQUS中，需要使用用户幅值子程序(UAMP)来实现。在UAMP程序的接口中需要关注如下几个形参变量:

(1)ampName，用户指定的幅值名称，用于区分当前幅值所赋给的自由度。

(2)time，由当前分析步所完成的时间和全部分析步所完成时间组成的数组。因为模型中只设一个分析步，两个时间相等，所以统一记作t，用于判定加载和卸载。

(3)ampValueOld，记作tA，增量步开始时刻的幅值。

(4)ampValueNew，记作t+ΔtA，增量步结束时刻的幅值。

(5)dt，记作 Δt，时间增量。

根据压板和模具的装配关系，在加载过程中，压板参考点的运动轨迹为

幅值tA 有 6 个分量tAi，i=1，2，…，6，相应的t+ΔtA 也有6 个分量，为t+ΔtAi，i=1，2，…，6，其中前3个分量分别表示沿X、Y、Z轴的平动位移，后3个分量表示绕相应坐标轴的转动。

对于移动分量t+ΔtAi(i=1，2，3)，有

式中，irtr表示位置向量rtr的第i个分量。

对于转动分量t+ΔtAi(i=4，5，6)，有

在卸载的过程中压板参考点的运动轨迹为

同样可以使用式(3)和式(4)计算t+ΔtA。

1.3 仿真试验设计

在空间曲线局部建立Frenet标架，曲线同时具有曲率和挠率，在研究平面弯曲变形时，常使用纯弯曲试验(曲率为固定值)。对比这种做法，将模具的脊线设计为圆柱螺线(图3)，即

图3 圆柱螺线模具

表1 均匀设计因素及水平

圆柱螺线的曲率κ和挠率τ均为固定值，分别为

使用U6(6d)均匀设计表设计仿真试验，水平顺序见表1。

2 仿真数据处理

2.1 中心线提取

在研究管材的回弹问题时，为了便于描述管子的空间形状变化，可以以管子的中心线作为研究对象，因为根据弯曲理论，圆管的中心线位于弯曲的中心层上，可以认为长度不变。但是管子的中心轴线不是一个几何实体，在ABAQUS中不能设定Set集合，无法提取出节点坐标。可以提取出的节点坐标是管壁上的网格节点，所以需要根据管壁上的节点坐标计算出截面中心点的坐标，即进行中心线提取。目前关于中心线提取的研究［13-15］多是关于从CT图像中提取血管的中心线，关于空间节点集的中心线提取没有查获相关的算法，为此需要设计相应的算法。

管子外表面的部分节点如图4所示。要处理的节点在空间位置上有一定的分布规律，近似沿若干个截平面分布;相邻两个截平面之间的距离大体相同，后一个截平面相对于前一个截平面既有移动又有转动，移动距离大体相同，但是转动角度和方向未知;数据的有效位数有限，如果只采用少数点进行计算会带来较大的误差。由于在仿真的过程中，管子的一端是固定的，所以以固定端作为提取中心点的初始位置，初始位置的中心点是可以估算的。

图4 管子外表面节点

本文根据以上分析提出了一种中心线提取方法，基本思路如下:以前一个截面的中心点为球心，建立半径为r的球面，沿截面的法向量方向移动球面直到和下一个截面的节点接触为止，根据初步接触点拟合平面，查找在平面附近的节点，再次拟合平面，这个平面就是新球面的大圆截面，可以算出圆心点和截面的法线方向。具体算法流程如下。

(1)初始化。给定拟合误差限ε，给定终止计算时的节点个数NT，给定球心点C(0)、半径r和大圆截面初始法向量V(0)，设定试探移动步长h，次数k=0，扩张系数为λ。

(2)计算节点集中包含点的个数N(k)，如果N(k)＜ NT，停机，C(i)(i=0，1，…，k)即为要提取的截面中心点。

(3)设定试探球心Ct=C(k)，半径r⇒R，前球心 CF=C(k)+hV(k)，后球心 CB=C(k)－ hV(k)。

(4)计算点集中节点到球心的距离dj(j=1，2，…，N(k))，统计在球面内点的个数mI和落在球面上点的个数mO。

(5)如果 mI＞ 0，则令 CF=Ct，Ct=(CF+CB)/2，转步骤(4)。

(6)如果 mO=0，则令 CB=Ct，Ct=(CF+CB)/2，转步骤(4)。

(7)如果 mO＜ 3，λ R⇒R，否则，转步骤(9)。

(8)查找落在球面上点的个数mO，转步骤(7)。

(9)使用落在球面上的点拟合平面Pt，方向设为V(k+1)，V(k+1)和V(k)的夹角小于π，并计算节点集中的点到平面Pt的有向距离，查找所有距离小于ε的点，构成截面点集。

(10)使用截面点集拟合半径为r的球面，球面圆心设为C(k+1)。

(11)k+1⇒k，转步骤(2)。

完成中心点的提取后，令s0=0，对于C(i)(i=1，2，…，k)，使用 si=si－1+|C(i)－ C(i－1)| 作为参数，将 C(i)(i=1，2，…，k)3个方向上的分量分别对si拟合三次函数，得到中心线的参数方程r(s)。

2.2 仿真结果分析

中心线参数方程的曲率κ(s)和挠率τ(s)分别为

相应的模具的曲率和挠率可以按式(7)和式(8)计算。

建立回弹前后曲率、挠率的对应关系，如图5所示。可以发现两组数据在试验范围内都接近于线性关系。

图5 回弹前后曲率、挠率关系

分别拟合回弹前后曲率、挠率，得

式中，κp、τp分别为回弹后的曲率和挠率，mm－1;κm、τm为回弹前的曲率和挠率，mm－1。

3 回弹补偿

3.1 补偿算法

根据式(11)和式(12)可计算出目标曲线上的每一点在模具脊线上对应点的曲率和挠率。相当于模具脊线上每一点的曲率κm和挠率τm都是已知的。从微分几何的角度考虑，对于在空间中不含逗留点的两条曲线，如果能够适当地选择公共自然参数s，使得在对应点两条曲线有相同的曲率和挠率，那么两条曲线是可以通过一个空间刚体运动相互重合的。所以，如果不考虑空间位置的差异，根据目标曲线可以唯一地计算出一条模具脊线。

由曲率和挠率求解空间曲线方程实际上是自然方程的积分问题，即首先求解Frenet－Serret方程:

其中，α为单位切向量，β为单位主法向量，γ为单位副法向量，它们构成曲线上某点的基本三棱形。

得到α、β、γ后再求解微分方程:

考虑到简化UAMP子程序编写的需求，在此使用一种半数值的补偿方法。另外，因为在定义曲线时，为了直观起见使用的是一般参数，故需要进行参数转化才能使用上述公式。具体的补偿过程如下。

(2)初始化。计算目标曲线的初始位置rp(0)，以及在初始位置的三棱标架 α0=r'p(0)，β0=α'0/|α'0|，γ0=α0× β0作为方程(13)和方程(14)的边界条件。

(3)求解Frenet－Serret方程组。将区间［0，sm］等距划分，使用四阶Runge－Kutta方法求解方程(13)，得到一系列αi的值。

(4)拟合切向量。使用三次多项式函数在区间［0，sm］拟合 αi，得到方程(14)的近似表达式r'=α。

(5)位置矢量积分。结合步骤(2)中的边界条件rp(0)，对r'=α直接积分得到rm。

进一步分析上述对曲线基本公式积分的补偿方法，补偿只要求知道曲率和挠率方程，而和具体的成形过程和截面形状无关。而且如果将成形后零件的曲率和挠率表示成模具曲率和挠率的函数，对自然方程积分还可以求解回弹预测问题，例如对于挠率恒等于0的情况(金属窄条的二维弯曲［16］)，也可以使用上述方法求解，但因为不需要在三维空间中求解微分方程，补偿过程可以大大简化。

3.2 补偿验证

验证使用的目标弯管的中心线参数方程为

可以计算出式(15)的曲率和挠率方程:

从式(16)和式(17)可以看出，目标曲线的曲率和挠率都不是固定值。先使用3.1节的方法进行补偿计算，按照1.1节生成模具，建立有限元模型，按1.2节中所述的公式编写UAMP子程序。然后进行有限元仿真，对仿真结果使用2.1节的方法提取中心点，将得到的中心点数据导入CATIA，并删除延长段的数据点和噪点，拟合空间样条。最后将目标曲线也导入CATIA中，如图6所示。因为使用目标曲线的初始位置矢量和初始基本三棱形作为求解微分方程的边界条件，所以可以认为得到的曲线和目标曲线是已经进行了粗略配准的。使用CATIA的测量功能，得到两条曲线的距离最大值约等于4.649mm。外形基本吻合，与文献［10］中的结果对比，在文献［10］中，轴线在俯视方向上，最大误差为5.15mm，由此可见，补偿效果与文献［10］处于同一水平。

图6 回弹补偿效果检验

4 结论

(1)建立了空间绕弯的有限元模型，并编写了绕弯的UAMP子程序。通过圆柱螺线型模具上的绕弯仿真试验，建立了回弹前后曲率和挠率的函数关系。

(2)提出了一种提取管状节点集截面中心点的方法，使用该方法从仿真结果中提取出了弯管的中心线。

(3)提出了利用空间曲线基本方程对空间弯管进行回弹补偿的方法。仿真试验证明，上述方法可以合理计算出空间绕弯模具的形状。

［1］Yang He，Li Heng，Zhang Zhiyong，et al.Advances and Trends on Tube Bending Forming Technologies［J］.Chinese Journal of Aeronautics，2012，25:1-12.

［2］詹梅，杨合，江志强.管材弯曲成形的国内外研究现状及发展趋势［J］.机械科学与技术，2004，23(12):1509-1514.Zhan Mei，Yang He，Jiang Zhiqiang.State of the Art of Research on Tube Bending Process［J］.Mechanical Science and Technology，2004，23(12):1509-1514.

［3］杨合，孙志超，林艳，等.管成形技术发展基础问题研究［J］.塑性工程学报，2001，8(2):83-85.Yang He，Sun Zhichao，Lin Yan，et al.Advanced Plastic Processing Technology and Research Progress on Tube Forming［J］.Journal of Plasticity Engineering，2001，8(2):83-85.

［4］Hermes M，Chatti S，Weinrich A，et al.Three－dimensional Bending of Profiles with Stress Superposition［J］.International Journal of Material Forming，2008，1(S1):133-136.

［5］李雁鹏，吴建军.非平面弯管成形过程的回弹补偿研究［J］.锻压技术，2009，34(1):89-92.Li Yanpeng，Wu Jianjun.Study on Spring－back Compensation for Non－plane Tube Bending Process［J］.Forging ＆ Stamping Technology，2009，34(1):89-92.

［6］Gantner P，Herbert B，Harrison D K，et al.Free－Bending—A New Bending Technique in the Hydroforming Process Chain［J］.Journal of Materials Processing Technology，2005，167:302-308.

［7］Chatti S，Hermes M，Tekkaya A E，et al.The New TSS Bending Process:3D Bending of Profiles with Arbitrary Cross－sections［J］.CIRP Annals－Manufacturing Technology，2010，59(1):315-318.

［8］李恒.多模具约束下薄壁管数控弯曲成形过程失稳起皱行为研究［D］.西安:西北工业大学，2007.

［9］寇永乐，杨合，詹梅，等.薄壁管数控弯曲应变的网格法研究［J］.中国机械工程，2006，16(S1):31-34.Kou Yongle，Yang He，Zhan Mei，et al.Research on Strain Distribution of Thin－walled Tube in NC Bending Process Using Grid Method［J］.China Mechanical Engineering，2006，16(S1):31-34.

［10］张深，吴建军.空间弯管的回弹预测［J］.航空学报，2011，32(5):953-960.Zhang Shen，Wu Jianjun.Spring－back Prediction of Non－planar Tube Bending［J］.Acta Aeronautica et Astronautica Sinica，2011，32(5):953-960.

［11］Al－Qureshi H A.Elastic－plastic Analysis of Tube Bending［J］.International Journal of Machine Tools＆ Manufacture，1999，39(1):87-104.

［12］Al－Qureshi H A，Russo A.Spring－back and Residual Stresses in Bending of Thin－walled Aluminum Tubes［J］.Materials and Design，2002，23(2):217-222.

［13］王胜军，付玲，康雁，等.一种基于管状特征的冠脉中心线自动提取算法［J］.东北大学学报(自然科学版)，2011，32(1):27-31.Wang Shengjun，Fu Ling，Kang Yan，et al.Automatic Centerline Extraction Based on Tubular Feature for Coronary Arteries［J］.Journal of Northeastern U-niversity(Natural Science)，2011，32(1):27-31.

［14］许燕，胡广书，商丽华，等.基于Hessian矩阵的冠状动脉中心线的跟踪算法［J］.清华大学学报(自然科学版)，2007，47(6):889-892.XuYan， HuGuangshu， ShangLihua， etal.Adaptive Tracking Extraction of Vessel Centerlines in Coronary Arteriograms Using Hessian Matrix［J］.J.Tsinghua Univ.(Sci.＆Tech.)，2007，47(6):889-892.

［15］李颖超，刘越，王涌天.基于多尺度Hessian矩阵和Gabor滤波的造影图像冠脉中心线提取［J］.中国医学影像技术，2007，23(1):133-136.Li Yingchao，Liu Yue，Wang Yongtian.Algorithm for Centerline Extraction of Coronary Arterial Tree in Coronary Angiographic Projections Based on Hessian Matrix and Gabor Filter［J］.Chin.J.Med.Imaging Technol.，2007，23(1):133-136.

［16］王晓林，周贤宾.金属板弹塑性非圆弧弯曲回弹的计算［J］.塑性工程学报，1996，3(4):27-33.Wang Xiaolin，Zhou Xianbin.Springback Analysis of Varied Curvature Elastic－plastic Bending for Sheet Metal［J］.Journal of Plasticity Engineering，1996，3(4):27-33.