劉永財， 鮑益東， 秦雪嬌， 劉玉琳， 陳文亮

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

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

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

1 彎曲與解耦計算的基本原理

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

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

2 基于預應力膜單元的中間構形算法

2.1 算法的基本流程

假設已計算出上一個時間步的板料構形，則在該構形的基礎上計算當前時間步的中間構形.本文以雙動沖壓成形為例進行說明.具體計算步驟如下：

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

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

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

(4) 計算中間構形.采用迭代方法求解預應力膜單元的平衡狀態(tài)方程，計算出非接觸區(qū)域板料中間構形的形狀，而接觸區(qū)域板料中間構形的形狀由凸模或凹模形狀來確定.

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

圖2 中間構形的計算流程Fig.2 Calculation flow of intermediate configuration

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

判斷板料節(jié)點與模具接觸時，首先采用搜索算法確定節(jié)點沿法線方向投影到模具表面的投影單元，進而確定投影點.在確定投影單元時，運用劃分空間格的方法搜索與板料節(jié)點距離最近的模具節(jié)點，根據(jù)模具表面節(jié)點與單元的拓撲關系來確定投影單元的范圍，從而縮小搜索范圍，提高搜索速度.

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

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

(1)

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

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

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

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

2.3 平衡方程的建立與求解

因為接觸區(qū)構形可以根據(jù)接觸狀態(tài)和模具形狀來確定，所以本文采用基于預應力膜單元建立平衡方程的方法來計算非接觸區(qū)構形.

表示中性面的平面應變，其中：

預應力膜單元在變形過程中的幾何關系為

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

(2)

式中：ε為應變；εL、εNL分別為線性應變和非線性應變；B為應變與位移的轉換矩陣；BL、BNL分別為線性應變和非線性應變分析的矩陣項，前者與a無關，后者是a的函數(shù).

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

σ=Dε+σ0

(3)

式中：D為彈性矩陣，其與材料屬性相關；σ0為膜單元的預應力.

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

(4)

式中：V為單元的體積；a*和ε*分別為虛位移和虛應變；f為局部坐標系下單元所受外力，在使用預應力膜單元進行有限元計算的過程中，f=0.

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

(5)

本文采用Newton-Raphson迭代法進行求解，即

(6)

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

KT=K0+Kσ+KNL

K0為小位移的線性剛度矩陣，與BL和D有關，KNL為大位移剛度矩陣，與BL、BNL及D有關，Kσ為幾何剛度矩陣，與應力σ及形函數(shù)對坐標的偏導數(shù)有關.

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

3 實例驗證

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

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

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

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

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

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

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

圖7 拉延深度200 mm時中間構形板料的厚度分布情況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 結論

(1) 對于金屬薄板成形，采用預應力膜單元的方法能夠快速準確地得出基于大步長靜力隱式有限元法的板料成形的中間構形.

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