熊威甘忠

西北工業(yè)大學(xué)現(xiàn)代設(shè)計與集成制造技術(shù)教育部重點實驗室，西安，710072

0 引言

作為一種重要的承載和輸送氣體、液體的輕質(zhì)結(jié)構(gòu)件，金屬彎管零件在高科技領(lǐng)域(如航空航天、造船、汽車、能源和醫(yī)療)得到了越來越多的關(guān)注［1-2］。為了將管材彎成具有一定彎曲半徑、角度和形狀的零件，彎管工藝作為管材塑性成形的重要分支，是生產(chǎn)輕質(zhì)結(jié)構(gòu)的一種十分重要的制造技術(shù)［1，3］。

以前的彎管研究主要集中在平面彎管領(lǐng)域，而對空間彎管(在彎曲段中，管材軸線的撓率不恒等于零的管彎零件)較少涉及。近年來，土木工程和交通系統(tǒng)中對三維彎曲鋼鋁構(gòu)件的需求日趨強烈，空間彎管因其在空間和流體力學(xué)上的優(yōu)勢，為設(shè)計者提供了新的選擇［4］。例如美國設(shè)計生產(chǎn)的太陽能STERLING發(fā)電機中散熱管系中的管型，管材軸線即為空間曲線［5］。文獻［6］提出了一種自由彎曲方法(free－bending)，該方法將管材擠過一個可以在兩個垂直方向偏轉(zhuǎn)的彎曲模腔，幾乎能夠成形任意外形的彎管，但是彎曲半徑要是管材半徑的2.5倍，并且對材料參數(shù)的變化較為敏感。另一種新的空間彎管方法［4，7］是扭矩疊加空間彎曲工藝(torque superposed spatial bending)，它基于滾彎工藝，通過動態(tài)調(diào)整彎曲外形和施加扭矩，來達到彎曲不同半徑和角度的目的。

一般而言，彎管成形的主要缺陷有失穩(wěn)起皺、壁厚過度減薄、截面畸變、拉裂和回彈［8-9］。目前，對空間彎管的研究主要集中在對小直徑管材的彎曲研究上［5-6，10］。相對于管材的破裂起皺和截面畸變，管材的回彈是主要的成形缺陷［10］。在平面彎管的研究中，如果使用理想彈塑性模型，并且假設(shè)在彎曲過程中管徑保持不變，回彈可以使用彎曲梁理論進行計算［11-12］。但因為忽略了三維應(yīng)力應(yīng)變，故純解析計算結(jié)果會和試驗結(jié)果相差較大［1］，所以有限元分析成為彎管回彈的主要方法。在空間彎管的回彈研究方面，文獻［5］將管材的空間軸線離散為若干段圓弧，并通過純彎曲試驗確定每段圓弧在回彈前的半徑，最終在UG軟件中拼接各段圓弧，得到繞彎模具型面。實質(zhì)上是將空間彎曲用若干個密切平面上的平面彎曲近似處理。文獻［10］發(fā)展了上述平面化的補償方法，將空間軸線分段投影到密切平面和從切平面上，同樣使用純彎曲試驗確定每段圓弧在回彈后的半徑，再將對應(yīng)的一對投影圓弧在空間中合成一段空間曲線，同樣使用UG拼合，預(yù)測了空間彎管回彈后的形狀。上述文獻中沒有論述平面化回彈預(yù)測和補償?shù)臄?shù)學(xué)依據(jù)，并且沒有涉及空間繞彎的數(shù)值仿真問題。

本文通過有限元仿真研究小直徑空間彎管的回彈問題。提出了一種具有數(shù)學(xué)依據(jù)的空間彎管回彈補償方法。

1 空間彎管數(shù)值模擬

1.1 有限元模型

空間彎管使用繞彎的成形方法，繞彎模具的幾何模型(見圖1)，是沿初始截面和終止截面之間的脊線放樣而成的。兩個截面分別位于XOY面和YOZ面上;脊線是一段空間曲線，起點為rs(0)=(x0，H －rm，0)，終點的Y坐標(biāo)比起點的Y坐標(biāo)小H，曲線的參數(shù)方程為

圖1 模具的幾何模型

空間彎管的有限元模型見圖2，管材使用C3D8R網(wǎng)格，管材外半徑rp=2.8mm，略小于型面半徑rm=3.0mm，這是因為將模具型面的截面半徑同管材半徑設(shè)為一致，在進行網(wǎng)格劃分時，兩者的網(wǎng)格有可能相互干涉而導(dǎo)致分析失敗。管材的內(nèi)徑為5.0mm，材料為鋼材，使用線性強化模型，設(shè)定彈性模量為210GPa，泊松比為0.33，屈服應(yīng)力σy=418MPa，硬化系數(shù)Hp=31 579GPa。模具和壓板的彈性模量比管材的彈性模量大一個數(shù)量級，以保證在成形過程中不發(fā)生過大的變形。

壓板是內(nèi)徑為6mm、厚1mm、旋轉(zhuǎn)角為90°的環(huán)扇形，參考點設(shè)在圓心。壓板比管材略長，這是出于有限元分析收斂性的考慮，長的壓板可以控制自由端，防止產(chǎn)生不確定的運動，使成形過程為準(zhǔn)靜態(tài)。

圖2 空間彎管有限元模型

模型使用ABAQUS/Standard求解，只設(shè)定一個分析步，時間為2s，在前一秒中壓板參考點沿軌跡從XOY面運動到Y(jié)OZ面，對管材加載;后一秒，壓板參考點沿原路徑返回XOY面，完成卸載。

1.2 用戶幅值子程序

空間彎管使用繞彎的方式成形時，需要控制壓板沿空間曲線運動，在彎管的過程中，壓板的6個自由度需要相互協(xié)調(diào)運動，都是時間的函數(shù)。在ABAQUS中，需要使用用戶幅值子程序(UAMP)來實現(xiàn)。在UAMP程序的接口中需要關(guān)注如下幾個形參變量:

(1)ampName，用戶指定的幅值名稱，用于區(qū)分當(dāng)前幅值所賦給的自由度。

(2)time，由當(dāng)前分析步所完成的時間和全部分析步所完成時間組成的數(shù)組。因為模型中只設(shè)一個分析步，兩個時間相等，所以統(tǒng)一記作t，用于判定加載和卸載。

(3)ampValueOld，記作tA，增量步開始時刻的幅值。

(4)ampValueNew，記作t+ΔtA，增量步結(jié)束時刻的幅值。

(5)dt，記作 Δt，時間增量。

根據(jù)壓板和模具的裝配關(guān)系，在加載過程中，壓板參考點的運動軌跡為

幅值tA 有 6 個分量tAi，i=1，2，…，6，相應(yīng)的t+ΔtA 也有6 個分量，為t+ΔtAi，i=1，2，…，6，其中前3個分量分別表示沿X、Y、Z軸的平動位移，后3個分量表示繞相應(yīng)坐標(biāo)軸的轉(zhuǎn)動。

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

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

對于轉(zhuǎn)動分量t+ΔtAi(i=4，5，6)，有

在卸載的過程中壓板參考點的運動軌跡為

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

1.3 仿真試驗設(shè)計

在空間曲線局部建立Frenet標(biāo)架，曲線同時具有曲率和撓率，在研究平面彎曲變形時，常使用純彎曲試驗(曲率為固定值)。對比這種做法，將模具的脊線設(shè)計為圓柱螺線(圖3)，即

圖3 圓柱螺線模具

表1 均勻設(shè)計因素及水平

圓柱螺線的曲率κ和撓率τ均為固定值，分別為

使用U6(6d)均勻設(shè)計表設(shè)計仿真試驗，水平順序見表1。

2 仿真數(shù)據(jù)處理

2.1 中心線提取

在研究管材的回彈問題時，為了便于描述管子的空間形狀變化，可以以管子的中心線作為研究對象，因為根據(jù)彎曲理論，圓管的中心線位于彎曲的中心層上，可以認(rèn)為長度不變。但是管子的中心軸線不是一個幾何實體，在ABAQUS中不能設(shè)定Set集合，無法提取出節(jié)點坐標(biāo)。可以提取出的節(jié)點坐標(biāo)是管壁上的網(wǎng)格節(jié)點，所以需要根據(jù)管壁上的節(jié)點坐標(biāo)計算出截面中心點的坐標(biāo)，即進行中心線提取。目前關(guān)于中心線提取的研究［13-15］多是關(guān)于從CT圖像中提取血管的中心線，關(guān)于空間節(jié)點集的中心線提取沒有查獲相關(guān)的算法，為此需要設(shè)計相應(yīng)的算法。

管子外表面的部分節(jié)點如圖4所示。要處理的節(jié)點在空間位置上有一定的分布規(guī)律，近似沿若干個截平面分布;相鄰兩個截平面之間的距離大體相同，后一個截平面相對于前一個截平面既有移動又有轉(zhuǎn)動，移動距離大體相同，但是轉(zhuǎn)動角度和方向未知;數(shù)據(jù)的有效位數(shù)有限，如果只采用少數(shù)點進行計算會帶來較大的誤差。由于在仿真的過程中，管子的一端是固定的，所以以固定端作為提取中心點的初始位置，初始位置的中心點是可以估算的。

圖4 管子外表面節(jié)點

本文根據(jù)以上分析提出了一種中心線提取方法，基本思路如下:以前一個截面的中心點為球心，建立半徑為r的球面，沿截面的法向量方向移動球面直到和下一個截面的節(jié)點接觸為止，根據(jù)初步接觸點擬合平面，查找在平面附近的節(jié)點，再次擬合平面，這個平面就是新球面的大圓截面，可以算出圓心點和截面的法線方向。具體算法流程如下。

(1)初始化。給定擬合誤差限ε，給定終止計算時的節(jié)點個數(shù)NT，給定球心點C(0)、半徑r和大圓截面初始法向量V(0)，設(shè)定試探移動步長h，次數(shù)k=0，擴張系數(shù)為λ。

(2)計算節(jié)點集中包含點的個數(shù)N(k)，如果N(k)＜ NT，停機，C(i)(i=0，1，…，k)即為要提取的截面中心點。

(3)設(shè)定試探球心Ct=C(k)，半徑r?R，前球心 CF=C(k)+hV(k)，后球心 CB=C(k)－ hV(k)。

(4)計算點集中節(jié)點到球心的距離dj(j=1，2，…，N(k))，統(tǒng)計在球面內(nèi)點的個數(shù)mI和落在球面上點的個數(shù)mO。

(5)如果 mI＞ 0，則令 CF=Ct，Ct=(CF+CB)/2，轉(zhuǎn)步驟(4)。

(6)如果 mO=0，則令 CB=Ct，Ct=(CF+CB)/2，轉(zhuǎn)步驟(4)。

(7)如果 mO＜ 3，λ R?R，否則，轉(zhuǎn)步驟(9)。

(8)查找落在球面上點的個數(shù)mO，轉(zhuǎn)步驟(7)。

(9)使用落在球面上的點擬合平面Pt，方向設(shè)為V(k+1)，V(k+1)和V(k)的夾角小于π，并計算節(jié)點集中的點到平面Pt的有向距離，查找所有距離小于ε的點，構(gòu)成截面點集。

(10)使用截面點集擬合半徑為r的球面，球面圓心設(shè)為C(k+1)。

(11)k+1?k，轉(zhuǎn)步驟(2)。

完成中心點的提取后，令s0=0，對于C(i)(i=1，2，…，k)，使用 si=si－1+|C(i)－ C(i－1)| 作為參數(shù)，將 C(i)(i=1，2，…，k)3個方向上的分量分別對si擬合三次函數(shù)，得到中心線的參數(shù)方程r(s)。

2.2 仿真結(jié)果分析

中心線參數(shù)方程的曲率κ(s)和撓率τ(s)分別為

相應(yīng)的模具的曲率和撓率可以按式(7)和式(8)計算。

建立回彈前后曲率、撓率的對應(yīng)關(guān)系，如圖5所示?？梢园l(fā)現(xiàn)兩組數(shù)據(jù)在試驗范圍內(nèi)都接近于線性關(guān)系。

圖5 回彈前后曲率、撓率關(guān)系

分別擬合回彈前后曲率、撓率，得

式中，κp、τp分別為回彈后的曲率和撓率，mm－1;κm、τm為回彈前的曲率和撓率，mm－1。

3 回彈補償

3.1 補償算法

根據(jù)式(11)和式(12)可計算出目標(biāo)曲線上的每一點在模具脊線上對應(yīng)點的曲率和撓率。相當(dāng)于模具脊線上每一點的曲率κm和撓率τm都是已知的。從微分幾何的角度考慮，對于在空間中不含逗留點的兩條曲線，如果能夠適當(dāng)?shù)剡x擇公共自然參數(shù)s，使得在對應(yīng)點兩條曲線有相同的曲率和撓率，那么兩條曲線是可以通過一個空間剛體運動相互重合的。所以，如果不考慮空間位置的差異，根據(jù)目標(biāo)曲線可以唯一地計算出一條模具脊線。

由曲率和撓率求解空間曲線方程實際上是自然方程的積分問題，即首先求解Frenet－Serret方程:

其中，α為單位切向量，β為單位主法向量，γ為單位副法向量，它們構(gòu)成曲線上某點的基本三棱形。

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

考慮到簡化UAMP子程序編寫的需求，在此使用一種半數(shù)值的補償方法。另外，因為在定義曲線時，為了直觀起見使用的是一般參數(shù)，故需要進行參數(shù)轉(zhuǎn)化才能使用上述公式。具體的補償過程如下。

(2)初始化。計算目標(biāo)曲線的初始位置rp(0)，以及在初始位置的三棱標(biāo)架 α0=r'p(0)，β0=α'0/|α'0|，γ0=α0× β0作為方程(13)和方程(14)的邊界條件。

(3)求解Frenet－Serret方程組。將區(qū)間［0，sm］等距劃分，使用四階Runge－Kutta方法求解方程(13)，得到一系列αi的值。

(4)擬合切向量。使用三次多項式函數(shù)在區(qū)間［0，sm］擬合 αi，得到方程(14)的近似表達式r'=α。

(5)位置矢量積分。結(jié)合步驟(2)中的邊界條件rp(0)，對r'=α直接積分得到rm。

進一步分析上述對曲線基本公式積分的補償方法，補償只要求知道曲率和撓率方程，而和具體的成形過程和截面形狀無關(guān)。而且如果將成形后零件的曲率和撓率表示成模具曲率和撓率的函數(shù)，對自然方程積分還可以求解回彈預(yù)測問題，例如對于撓率恒等于0的情況(金屬窄條的二維彎曲［16］)，也可以使用上述方法求解，但因為不需要在三維空間中求解微分方程，補償過程可以大大簡化。

3.2 補償驗證

驗證使用的目標(biāo)彎管的中心線參數(shù)方程為

可以計算出式(15)的曲率和撓率方程:

從式(16)和式(17)可以看出，目標(biāo)曲線的曲率和撓率都不是固定值。先使用3.1節(jié)的方法進行補償計算，按照1.1節(jié)生成模具，建立有限元模型，按1.2節(jié)中所述的公式編寫UAMP子程序。然后進行有限元仿真，對仿真結(jié)果使用2.1節(jié)的方法提取中心點，將得到的中心點數(shù)據(jù)導(dǎo)入CATIA，并刪除延長段的數(shù)據(jù)點和噪點，擬合空間樣條。最后將目標(biāo)曲線也導(dǎo)入CATIA中，如圖6所示。因為使用目標(biāo)曲線的初始位置矢量和初始基本三棱形作為求解微分方程的邊界條件，所以可以認(rèn)為得到的曲線和目標(biāo)曲線是已經(jīng)進行了粗略配準(zhǔn)的。使用CATIA的測量功能，得到兩條曲線的距離最大值約等于4.649mm。外形基本吻合，與文獻［10］中的結(jié)果對比，在文獻［10］中，軸線在俯視方向上，最大誤差為5.15mm，由此可見，補償效果與文獻［10］處于同一水平。

圖6 回彈補償效果檢驗

4 結(jié)論

(1)建立了空間繞彎的有限元模型，并編寫了繞彎的UAMP子程序。通過圓柱螺線型模具上的繞彎仿真試驗，建立了回彈前后曲率和撓率的函數(shù)關(guān)系。

(2)提出了一種提取管狀節(jié)點集截面中心點的方法，使用該方法從仿真結(jié)果中提取出了彎管的中心線。

(3)提出了利用空間曲線基本方程對空間彎管進行回彈補償?shù)姆椒ā７抡嬖囼炞C明，上述方法可以合理計算出空間繞彎模具的形狀。

［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］詹梅，楊合，江志強.管材彎曲成形的國內(nèi)外研究現(xiàn)狀及發(fā)展趨勢［J］.機械科學(xué)與技術(shù)，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］楊合，孫志超，林艷，等.管成形技術(shù)發(fā)展基礎(chǔ)問題研究［J］.塑性工程學(xué)報，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］.鍛壓技術(shù)，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］李恒.多模具約束下薄壁管數(shù)控彎曲成形過程失穩(wěn)起皺行為研究［D］.西安:西北工業(yè)大學(xué)，2007.

［9］寇永樂，楊合，詹梅，等.薄壁管數(shù)控彎曲應(yīng)變的網(wǎng)格法研究［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］張深，吳建軍.空間彎管的回彈預(yù)測［J］.航空學(xué)報，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］.東北大學(xué)學(xué)報(自然科學(xué)版)，2011，32(1):27-31.Wang Shengjun，F(xiàn)u 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］.清華大學(xué)學(xué)報(自然科學(xué)版)，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］.中國醫(yī)學(xué)影像技術(shù)，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］.塑性工程學(xué)報，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.