板料成形快速模擬的中間構(gòu)形構(gòu)造方法

劉永財(cái)， 鮑益東， 秦雪嬌， 劉玉琳， 陳文亮

(1. 南京航空航天大學(xué) 江蘇省精密與微細(xì)制造技術(shù)重點(diǎn)實(shí)驗(yàn)室， 南京 210016；2. 曲靖師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院， 云南 曲靖 655011； 3. 美國(guó)工程技術(shù)聯(lián)合公司， 南京 210016)

在金屬板料成形的模擬方法中[1]，為了彌補(bǔ)一步逆成形法模擬的應(yīng)力精度不高的問(wèn)題[2-3]，提出了引入中間構(gòu)形的多步成形有限元模擬方法并得到了廣泛應(yīng)用[4-5].其中，構(gòu)造中間構(gòu)形是多步成形法中的一個(gè)關(guān)鍵技術(shù).相關(guān)研究包括：Lee等[6]采用線彈性反向映射法計(jì)算中間構(gòu)形，并將中間構(gòu)形視為從初始板料構(gòu)形到最終板料構(gòu)形的線彈性變形過(guò)程的中間過(guò)程；Guo等[7]考慮軸對(duì)稱問(wèn)題，利用最小面積法構(gòu)造中間構(gòu)形，將三維中間構(gòu)形轉(zhuǎn)化為二維中間構(gòu)形；劉偉杰等[8-9]采用偽最小面積法構(gòu)造滑移約束面，以降低對(duì)零件形狀的依賴性；Shirin等[10]將最終零件的網(wǎng)格向初始坯料所在的平面投影并利用網(wǎng)格展開(kāi)進(jìn)行求解；陳偉等[11]采用比例插值方法在初始構(gòu)形與最終構(gòu)形之間進(jìn)行插值以構(gòu)造中間構(gòu)形.

本文提出了一種基于預(yù)應(yīng)力膜單元的板料中間構(gòu)形的平衡迭代計(jì)算方法，以適應(yīng)復(fù)雜的零件形狀，并將其集成到自主研發(fā)的基于大步長(zhǎng)靜力隱式有限元模擬分析軟件QuickForm中，通過(guò)典型汽車零件的應(yīng)用驗(yàn)證其有效性和準(zhǔn)確性.

1 彎曲與解耦計(jì)算的基本原理

對(duì)于金屬薄板成形來(lái)說(shuō)，由于剪切變形對(duì)整個(gè)板料成形結(jié)果的影響較小，所以板料變形被看作由彎曲變形和拉伸變形構(gòu)成.解耦是將板料的彎曲變形和拉伸變形過(guò)程分開(kāi)進(jìn)行求解，如圖1所示.采用解耦技術(shù)可將板料成形的模擬過(guò)程分為彎曲變形和拉伸變形來(lái)求解.前者是求解板料在當(dāng)前模具作用下的構(gòu)形，因此，可以采用一種相對(duì)簡(jiǎn)單的方法求解以提高計(jì)算效率；后者是在彎曲變形所求板料中間構(gòu)形的基礎(chǔ)上，將板料的拉伸變形限定在中間構(gòu)形上，并對(duì)包含彎曲效應(yīng)的膜剛度所引起的拉伸變形進(jìn)行求解.由于在拉伸變形的平衡迭代過(guò)程中，板料節(jié)點(diǎn)的運(yùn)動(dòng)被限定在由中間構(gòu)形確定的滑移約束面上，只需考慮滑移約束面上兩個(gè)方向的位移，所以降低了節(jié)點(diǎn)的自由度線性方程組的剛度系數(shù)矩陣的條件數(shù)，從而有效地提高了求解效率.

圖1 解耦計(jì)算的原理示意圖Fig.1 Diagram of principle of decoupling calculation

2 基于預(yù)應(yīng)力膜單元的中間構(gòu)形算法

2.1 算法的基本流程

假設(shè)已計(jì)算出上一個(gè)時(shí)間步的板料構(gòu)形，則在該構(gòu)形的基礎(chǔ)上計(jì)算當(dāng)前時(shí)間步的中間構(gòu)形.本文以雙動(dòng)沖壓成形為例進(jìn)行說(shuō)明.具體計(jì)算步驟如下：

(1) 判斷板料與壓邊圈以及凹模的接觸情況.尋找板料與壓邊圈相互接觸的節(jié)點(diǎn)，接觸節(jié)點(diǎn)應(yīng)附著在凹模上，并對(duì)其進(jìn)行固定約束，且在計(jì)算過(guò)程中節(jié)點(diǎn)的位置不變.

(2) 判斷板料與凸模的接觸情況.在雙動(dòng)沖壓過(guò)程中，凹模始終固定，凸模每行進(jìn)一個(gè)步長(zhǎng)，與板料就有可能發(fā)生穿透現(xiàn)象.本文采用節(jié)點(diǎn)接觸判斷算法將穿透凸模表面的節(jié)點(diǎn)通過(guò)投影操作拉至凸模表面.

(3) 確定非接觸區(qū).根據(jù)位移計(jì)算拉至凸模表面的板料節(jié)點(diǎn)的內(nèi)力，以確定節(jié)點(diǎn)力的狀態(tài).如果節(jié)點(diǎn)受到向外的拉力，則將其添加到非接觸區(qū)節(jié)點(diǎn)集合中；否則，將其約束在凸模表面，并添加到約束節(jié)點(diǎn)集合中.

(4) 計(jì)算中間構(gòu)形.采用迭代方法求解預(yù)應(yīng)力膜單元的平衡狀態(tài)方程，計(jì)算出非接觸區(qū)域板料中間構(gòu)形的形狀，而接觸區(qū)域板料中間構(gòu)形的形狀由凸?；虬寄Ｐ螤顏?lái)確定.

重復(fù)步驟(1)～(4)，直到中間構(gòu)形的所有單元都滿足預(yù)應(yīng)力膜單元的平衡條件，從而完成當(dāng)前增量步的中間構(gòu)形構(gòu)造，算法的流程如圖2所示.類似地，可將上述算法推廣到單動(dòng)沖壓成形的模擬.

圖2 中間構(gòu)形的計(jì)算流程Fig.2 Calculation flow of intermediate configuration

2.2 節(jié)點(diǎn)與模具接觸判斷算法

判斷板料節(jié)點(diǎn)與模具接觸時(shí)，首先采用搜索算法確定節(jié)點(diǎn)沿法線方向投影到模具表面的投影單元，進(jìn)而確定投影點(diǎn).在確定投影單元時(shí)，運(yùn)用劃分空間格的方法搜索與板料節(jié)點(diǎn)距離最近的模具節(jié)點(diǎn)，根據(jù)模具表面節(jié)點(diǎn)與單元的拓?fù)潢P(guān)系來(lái)確定投影單元的范圍，從而縮小搜索范圍，提高搜索速度.

確定投影點(diǎn)后，根據(jù)幾何關(guān)系判斷節(jié)點(diǎn)是否穿透模具表面，如圖3所示.具體計(jì)算公式為

g(xs)=(xs-xm)·n(xm)

(1)

式中：g為節(jié)點(diǎn)穿透量；xs為板料節(jié)點(diǎn)位置；xm為板料節(jié)點(diǎn)xs沿節(jié)點(diǎn)法線方向投影至模具表面的節(jié)點(diǎn)位置；n(xm)表示節(jié)點(diǎn)xm所在投影單元的外法線方向.

圖3 節(jié)點(diǎn)穿透判斷示意圖Fig.3 Diagram of penetrating judgment

進(jìn)行節(jié)點(diǎn)穿透判斷時(shí)，若g(xs)>0，則表示板料節(jié)點(diǎn)xs并沒(méi)有穿透模具單元；若g(xs)<0，則表示板料節(jié)點(diǎn)xs已經(jīng)穿透模具表面進(jìn)入模具內(nèi)部.此時(shí)，需要將該節(jié)點(diǎn)拉回到模具表面，使其成為附著在模具表面的固定點(diǎn).

經(jīng)過(guò)幾何處理后，所有穿透模具表面的節(jié)點(diǎn)都被拉回到模具表面.但是，實(shí)際上并非所有節(jié)點(diǎn)都附著在模具主表面，此時(shí)，需要進(jìn)一步確定節(jié)點(diǎn)狀態(tài).針對(duì)拉回到模具表面的節(jié)點(diǎn)，根據(jù)投影前后的位置變化來(lái)計(jì)算相應(yīng)的節(jié)點(diǎn)內(nèi)力、內(nèi)力方向與投影單元外法線的夾角.如果該夾角小于等于90°，則說(shuō)明節(jié)點(diǎn)受到向外的拉力，將其釋放為自由節(jié)點(diǎn)，并添加到非接觸區(qū)節(jié)點(diǎn)集合中；否則，該節(jié)點(diǎn)仍然為附著在模具表面的固定點(diǎn).

2.3 平衡方程的建立與求解

因?yàn)榻佑|區(qū)構(gòu)形可以根據(jù)接觸狀態(tài)和模具形狀來(lái)確定，所以本文采用基于預(yù)應(yīng)力膜單元建立平衡方程的方法來(lái)計(jì)算非接觸區(qū)構(gòu)形.

表示中性面的平面應(yīng)變，其中：

預(yù)應(yīng)力膜單元在變形過(guò)程中的幾何關(guān)系為

ε=εL+εNL=Ba=(BL+BNL)a

(2)

式中：ε為應(yīng)變；εL、εNL分別為線性應(yīng)變和非線性應(yīng)變；B為應(yīng)變與位移的轉(zhuǎn)換矩陣；BL、BNL分別為線性應(yīng)變和非線性應(yīng)變分析的矩陣項(xiàng)，前者與a無(wú)關(guān)，后者是a的函數(shù).

盡管預(yù)應(yīng)力膜單元在變形過(guò)程中的位移和撓度都較大，但由于結(jié)構(gòu)的應(yīng)變不大，所以單元具有大位移、小應(yīng)變的特點(diǎn)[12].根據(jù)這一特點(diǎn)，假設(shè)材料的應(yīng)力-應(yīng)變關(guān)系滿足如下線彈性關(guān)系：

σ=Dε+σ0

(3)

式中：D為彈性矩陣，其與材料屬性相關(guān)；σ0為膜單元的預(yù)應(yīng)力.

在單元所在的局部坐標(biāo)系下，根據(jù)幾何關(guān)系和本構(gòu)關(guān)系，利用虛功原理建立如下關(guān)系式：

(4)

式中：V為單元的體積；a*和ε*分別為虛位移和虛應(yīng)變；f為局部坐標(biāo)系下單元所受外力，在使用預(yù)應(yīng)力膜單元進(jìn)行有限元計(jì)算的過(guò)程中，f=0.

由虛位移a*的任意性所得非線性問(wèn)題的平衡方程為

(5)

本文采用Newton-Raphson迭代法進(jìn)行求解，即

(6)

式中：ω為迭代收斂松弛因子；an為第n迭代步的位移向量；Δan為第n迭代步的位移增量，n=0,1,2,…；Ψ(an)為an的函數(shù)；KT為切線剛度矩陣，且

KT=K0+Kσ+KNL

K0為小位移的線性剛度矩陣，與BL和D有關(guān)，KNL為大位移剛度矩陣，與BL、BNL及D有關(guān)，Kσ為幾何剛度矩陣，與應(yīng)力σ及形函數(shù)對(duì)坐標(biāo)的偏導(dǎo)數(shù)有關(guān).

在迭代過(guò)程中，采用位移收斂準(zhǔn)則與最大迭代次數(shù)相結(jié)合的方法進(jìn)行收斂判斷.若節(jié)點(diǎn)位移范數(shù)滿足位移收斂準(zhǔn)則，則求解結(jié)束；若迭代次數(shù)大于最大迭代次數(shù)時(shí)仍未滿足位移收斂準(zhǔn)則，則認(rèn)為迭代發(fā)散，求解不成功，此時(shí)，需要采用縮小步長(zhǎng)的方法重新構(gòu)造中間構(gòu)形.

3 實(shí)例驗(yàn)證

本文以NumiSheet 2005’標(biāo)準(zhǔn)考題的汽車前翼子板零件為例來(lái)說(shuō)明所提算法的有效性.其中，初始板料幾何尺寸如圖4(a)所示.圖4(b)為沖壓初始時(shí)刻板料與模具的裝配圖.其中，上層為凸模，中間為壓邊圈和板料，下層為凹模.板料的初始厚度為 0.7 mm，包括 1 313 個(gè)節(jié)點(diǎn)和 2 468 個(gè)三角形單元，所選板料的彈性模量為210 GPa，泊松比為 0.3，密度為 7.8 g/cm3，摩擦系數(shù)為 0.2，其單向拉伸的應(yīng)力與應(yīng)變關(guān)系為σ=545.0(0.004+ε)0.263，凸模的總行程為 262.30 mm.

為了驗(yàn)證本文提出的算法的有效性，將本文算法集成到自主研發(fā)的基于大步長(zhǎng)靜力隱式有限元模擬分析軟件QuickForm中，并將計(jì)算結(jié)果與顯示動(dòng)力算法的模擬分析軟件LS-DYNA的模擬結(jié)果進(jìn)行對(duì)比.取零件的截面線A-A′(如圖5所示)，并對(duì)比在A-A′ 線位置的中間構(gòu)形形狀.

由于重力和壓邊圈的作用，使得拉延深度小于150 mm的中間構(gòu)形形狀與初始重力作用的中間構(gòu)形形狀相似，所以本文只比較拉延深度為170和200 mm時(shí)中間構(gòu)形的計(jì)算結(jié)果.其中：當(dāng)拉延深度為170 mm時(shí)，LS-DYNA軟件中含 14 856 個(gè)節(jié)點(diǎn)和 26 732 個(gè)三角形單元，總面積為 1.015 m2；QuickForm軟件中含 6 842 個(gè)節(jié)點(diǎn)和 12 464 個(gè)三角形單元，總面積為 1.012 m2，兩者預(yù)測(cè)的總面積相當(dāng)接近，相對(duì)誤差為3%.拉延深度為200時(shí)，LS-DYNA軟件和QuickForm軟件計(jì)算的單元總面積分別為 1.029 和 1.025 m2，相對(duì)誤差為4%.沿著截面線A-A′比較拉延深度為170和200 mm時(shí)兩種軟件計(jì)算的xOz平面的中間構(gòu)形形狀如圖6所示.由圖6可見(jiàn)，在不同的拉延深度下，兩種算法預(yù)測(cè)的中間構(gòu)形形狀很接近，即兩種算法均能夠合理地計(jì)算當(dāng)前拉延深度下零件的成形形狀，從而驗(yàn)證了采用預(yù)應(yīng)力膜單元構(gòu)造中間構(gòu)形方法的有效性.圖7示出了拉延深度為200 mm時(shí)兩種軟件計(jì)算的中間構(gòu)形板料厚度分布情況.可以看出，兩者的厚度分布趨勢(shì)相近，板料的減薄區(qū)域和增厚區(qū)域相一致.

圖4 前翼子板尺寸及其與模具裝配圖Fig.4 Front fender module and assembly drawing

圖5 截面線示意圖Fig.5 Diagram of section line

在英特爾酷睿5處理器、頻率 3.20 GHz、內(nèi)存8 GB 的計(jì)算機(jī)上進(jìn)行計(jì)算，使用QuickForm軟件進(jìn)行沖壓成形模擬，耗時(shí)為 90.24 s，在合模狀態(tài)下的成形零件包含 11 768 個(gè)節(jié)點(diǎn)和 21 983 個(gè)三角形單元.同時(shí)，在相同的計(jì)算機(jī)上，利用LS-DYNA軟件進(jìn)行沖壓成形模擬，耗時(shí)為 426.68 s，圖8所示為采用兩種軟件計(jì)算的成形結(jié)果.由圖8可以看出，兩者的厚度分布趨勢(shì)一致，但LS-DYNA軟件比QuickForm軟件具有更多的板料厚度分布層數(shù)，這是由于前者采用了加密級(jí)數(shù)較多的網(wǎng)格的緣故.圖9示出了采用兩種軟件在合模狀態(tài)下計(jì)算的截面線A-A′上的厚度分布情況.可以看出，采用基于大步長(zhǎng)靜力隱式有限元算法(軟件QuickForm)與采用顯式動(dòng)力算法(軟件LS-DYNA)所獲厚度分布曲線吻合，兩者在增厚區(qū)域和減薄區(qū)域的厚度分布基本一致.

圖6 兩種拉延深度下中間構(gòu)形的截面線對(duì)比Fig.6 Comparison of section lines at 170 and 200 mm drawing depth

圖7 拉延深度200 mm時(shí)中間構(gòu)形板料的厚度分布情況Fig.7 Distribution of thickness in 200 mm drawing depth

圖8 合模狀態(tài)下成形板料的厚度分布情況Fig.8 Distribution of thickness in clamping status

圖9 合模狀態(tài)下截面線A-A′上的厚度分布情況Fig.9 Distributions of thickness with section line A-A′ in clamping status

4 結(jié)論

(1) 對(duì)于金屬薄板成形，采用預(yù)應(yīng)力膜單元的方法能夠快速準(zhǔn)確地得出基于大步長(zhǎng)靜力隱式有限元法的板料成形的中間構(gòu)形.

(2) 所構(gòu)造的中間構(gòu)形可作為板料成形快速模擬中拉伸變形計(jì)算過(guò)程的滑移約束面，使用基于大步長(zhǎng)靜力隱式有限元算法計(jì)算滑移約束面的拉伸變形，可以得到合理的成形模擬結(jié)果.