基于流動理論的板材多步成形模擬方法

張 向 奎， 劉 偉 杰， 胡 平

( 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室, 遼寧 大連 116024 )

0 引 言

經(jīng)濟(jì)發(fā)展和環(huán)境保護(hù)日益增長的需求，使得生產(chǎn)更加輕量化、幾何形狀復(fù)雜和低消耗費用的車用機(jī)械零件變得日益迫切．通過冷熱沖壓工藝得到的大尺寸的金屬零件往往經(jīng)歷了多軸加載、高應(yīng)變率和高溫等復(fù)雜的加載條件．基于有限元的數(shù)值模擬方法可以對設(shè)計的零件提供零件成形性分析、起皺和破裂等缺陷估計，以輔助零件制造和優(yōu)化工藝參數(shù)．

多步成形模擬有限元方法，通過快速構(gòu)造若干的中間構(gòu)型，基于流動理論的本構(gòu)積分模型，以克服采用全量理論的一步法對最終零件的應(yīng)力估計精度不足的問題．其待解決的主要問題有：構(gòu)造滑移約束曲面、生成中間構(gòu)型初始解和中間構(gòu)型的平衡迭代．一步法是由Guo等[1]于1990年提出的采用全量理論并廣泛應(yīng)用于零件初始設(shè)計階段的一種快速有限元模擬方法．大量的數(shù)值驗證表明，該方法可以提供較好的應(yīng)變度量，但是由于僅考慮變形終了狀態(tài)而忽略了變形歷史導(dǎo)致應(yīng)力度量不準(zhǔn)確．多步成形模擬方法的工作最初見于Majlessi等[2-3]的文章，其采用基于最小勢能原理和比例構(gòu)造中間構(gòu)型的方法，用膜單元分析了簡單軸對稱零件的成形．隨后Lee等[4]在已知幾何中間變形狀態(tài)的情況下，應(yīng)用多步模擬方法對方盒和油底殼零件進(jìn)行了分析并優(yōu)化了初始板料形狀．Kim等[5-6]提出根據(jù)接觸節(jié)點構(gòu)造滑移約束曲面和截面線映射構(gòu)造初始解的方法分析了S-Rail零件的成形過程．Guo等[7]提出了通過求解空間曲面面積最小構(gòu)造滑移約束曲面的擬一步法，來提高應(yīng)力的度量精度．Huang等[8]提出了修改的截面線弧長法來獲得中間構(gòu)型的初始解．Li等[9]提出了一種直接標(biāo)量算法以實現(xiàn)快速更新多步模擬法構(gòu)型間的應(yīng)力．Tang等[10]將擬一步法擴(kuò)展到一般三維狀態(tài)下，并提出一種針對中間構(gòu)型間平衡迭代時的空間滑移策略．Halouani等[11]提出了自由曲面糾正的方法構(gòu)造中間構(gòu)型的初始解，并將擬一步法應(yīng)用到軸對稱冷軋成形分析中．

本文在一步法的基礎(chǔ)上提出多步成形模擬方法來提高對最終零件的應(yīng)力評估精度．

1 運動幾何關(guān)系

考慮如圖1所示的3個構(gòu)型：初始平板、中間構(gòu)型和最終零件間的非線性運動幾何關(guān)系．以初始平板構(gòu)型C0為參考構(gòu)型的C1構(gòu)型上的應(yīng)變增量度量幾何關(guān)系與一步逆成形法的參考平板構(gòu)型的非線性幾何關(guān)系是一致的．為引入變形后參考構(gòu)型的彎曲效應(yīng)，本文僅推導(dǎo)變形后參考中間構(gòu)型C1的最終零件構(gòu)型C上的應(yīng)變增量度量．

圖1 運動關(guān)系示意圖

考慮最終零件構(gòu)型的中性層上的物質(zhì)點p和沿其法向無限接近的點q：

(1)

式中：n和n1分別為物質(zhì)點p和p1處的法向，dn1不為0反映了參考中間構(gòu)型C1的彎曲程度．式(1)寫為如下的矩陣形式：

Fx=(o1+zn，xo2+zn，yn)=T(I+zb)

(z1n1，xz1n1，y0)

(2)

由此可以建立中間構(gòu)型與最終構(gòu)型間的變形梯度張量：

(3)

同樣的，左柯西-格林變形張量的逆可以表示為

(4)

式中系數(shù)如下：

a=A0(1+zr)2+2B0sz(1+zr)+C0(sz)2；

b=A0(1+zr)sz+B0(1+zr)(1+zt)+

B0(sz)2+C0sz(1+zt)；

c=A0(sz)2+2B0sz(1+zt)+C0(1+zt)2

(5)

和

A0=a0-z1n1，x(2(o1-up，x)+z1(n1，x))；

B0=b0-z1((o2-up，y)n1，x+(o1-up，x)n1，y+

z1n1，xn1，y)；

C0=c0-z1n1，y(2(o2-up，y)+z1(n1，y))

(6)

其中a0、b0和c0與一步逆成形法[1]中的系數(shù)一致．

從式(6)中可以看出，當(dāng)z1=0，物質(zhì)點位于中性面時，該應(yīng)變增量度量退化為一步逆成形法的幾何關(guān)系；當(dāng)z1≠0，物質(zhì)點位于厚度方向時，參考構(gòu)型是已變形彎曲的中間構(gòu)型C1，由于彎曲效應(yīng)此有限應(yīng)變增量度量不再等同于一步逆成形法的應(yīng)變度量方式．

2 本構(gòu)方程

本文采用厚向異性彈塑性耦合各向同性硬化(Power-Law硬化)材料模型，結(jié)合Hill48塑性勢、關(guān)聯(lián)流動法則和隱式返回算法更新應(yīng)力與應(yīng)變狀態(tài)．總變形梯度張量F根據(jù)乘法分解為彈性Fe和塑性Fp部分，為滿足客觀性要求，應(yīng)變與應(yīng)力等度量均在共旋系下操作，考慮到金屬成形中的小彈性應(yīng)變假設(shè)，應(yīng)變率可近似按加法

e+Dp分解為彈性部分Jaumann率

ε．

e和塑性部分Dp．

Hill48的二次塑性勢函數(shù)定義為

(7)

式中:σ是柯西應(yīng)力張量；P是各向異性矩陣；q是用于計算各向同性硬化的內(nèi)變量，通常取為累積塑性應(yīng)變．

根據(jù)關(guān)聯(lián)流動法則，塑性應(yīng)變率可表示為

Dp=λ．Pσ(σPσ)1/2

(8)

式中：λ為拉格朗日塑性乘子，并通過一致性條件

f．

(σ，q)=0求得

λ．=?fT?σCε．?fT?σC?f?σ-?f?q

(9)

利用上式可求得應(yīng)力率為

σ．=Cepε．=C-C?f?σ?fT?σC?f?σTC?f?σ-?f?q?è??????÷÷÷÷ε．

(10)

采用經(jīng)典彈性預(yù)測塑性返回算法以更新應(yīng)力狀態(tài)．(n+1)時刻第i次迭代的更新應(yīng)力可由下式計算：

(11)

多步成形模擬方法采用較少的中間構(gòu)型，使得構(gòu)型間的有限應(yīng)變增量取值較大，從而導(dǎo)致應(yīng)力更新的圖形返回算法迭代收斂困難．鑒于在一步法中廣泛采用的彈性厚向異性假設(shè)，可以得到泊松比與厚向異性系數(shù)存在的關(guān)系：

ν=r/(1+r)

(12)

3 中間構(gòu)型

3.1 滑移約束曲面

鑒于采用最小曲面作為目標(biāo)函數(shù)的優(yōu)化模型是極難以進(jìn)行求解的，作者等提出了采用擬最小面積法以獲得可以高效求解的優(yōu)化模型[12]．該方法假定滑移約束曲面是在凸模、凹模約束下的所有可能形狀中面積平方最小的曲面，則目標(biāo)函數(shù)修改為

(13)

經(jīng)過代數(shù)運算，單元面積的平方可表示為

(14)

其中ξe是單元的基本變量向量，He和Ce類似于單元的剛度矩陣和內(nèi)力向量．

同樣的，按組裝單元剛度矩陣和內(nèi)力向量的方式整合該面積坐標(biāo)矩陣和向量，則目標(biāo)函數(shù)式(13)變?yōu)?/p>

(15)

其中ξ、H和C分別對應(yīng)于結(jié)構(gòu)的基本變量向量、面積矩陣和面積向量．原問題即轉(zhuǎn)化為一個可在多項式時間內(nèi)求解的二次規(guī)劃問題．本文采用Matlab/quadprog( )進(jìn)行求解．

此時，人為假定容許誤差E，則任何滿足條件

ξ

(16)

的節(jié)點處于與凹模或凸模接觸的狀態(tài)．據(jù)此，可遍歷節(jié)點篩選得到與工具網(wǎng)格接觸的節(jié)點，基于這些接觸節(jié)點采用直接連接接觸邊界或Delaunay三角網(wǎng)格的方法重新構(gòu)造滑移曲面的側(cè)壁部分，以改善擬最小面積法構(gòu)造的側(cè)壁，本文稱為強(qiáng)化的PMA方法(EPMA方法)．

3.2 構(gòu)造初始解

由上一節(jié)構(gòu)造的中間構(gòu)型的滑移約束曲面僅為幾何曲面，本節(jié)則采用直接映射物質(zhì)點到滑移曲面上的方法生成初始的中間構(gòu)型，即中間構(gòu)型的初始解．其基本思想為構(gòu)造的中間構(gòu)型是成形過程中的臨時狀態(tài)，則其對應(yīng)的初始坯料與最終已知零件的初始坯料是一致的．由此可以求得中間構(gòu)型的初始解，其主要步驟為

(2)同樣求解最終構(gòu)型C得到初始平板構(gòu)型C0；

圖2 物質(zhì)點映射示意圖

(4)利用物質(zhì)點的唯一標(biāo)識，在滑移曲面構(gòu)型Ci上的單元編號ep內(nèi)，依據(jù)面積坐標(biāo)P(Ai，Aj，Ak)創(chuàng)建點P變形后對應(yīng)物質(zhì)點P′；

(5)循環(huán)投影所有物質(zhì)點，則得到了位于滑移約束曲面上的中間構(gòu)型的初始解．

該過程可進(jìn)行多次直至前后兩次初始解趨于一致，以避免過大的工藝補(bǔ)充面阻礙材料流動而影響一步法求解的坯料精度．

3.3 中間構(gòu)型的平衡迭代

通過以上擬最小面積法和初始解構(gòu)造方法，可快速生成若干真實的中間構(gòu)型，然而注意到初始解構(gòu)造過程中采用的一步法是基于全量理論的，因此要進(jìn)行采用流動理論更新應(yīng)力方法的中間構(gòu)型間的平衡迭代．

考慮一階線性三角形單元，其單元剛度矩陣和內(nèi)力向量為

(17)

式中：he和Ae為單元的厚度與面積，B為局部系下的單元應(yīng)變操作矩陣．

由于平衡迭代過程中，中間構(gòu)型上的物質(zhì)點被約束在滑移曲面上移動，通過罰函數(shù)的方法施加單元剛度矩陣約束：

Ke=BTCepBheAe+γNTN

(18)

式中：γ為經(jīng)試錯法確定的罰函數(shù)系數(shù)，N是單元的節(jié)點法向矩陣．

利用從全局系到局部系的單元坐標(biāo)變換矩陣Q，得到全局系下的單元剛度矩陣和內(nèi)力向量：

(19)

經(jīng)組裝可得到結(jié)構(gòu)的殘余力向量：

R(U(i))=Fext-Fint(U(i))=KΔU≠0

(20)

采用Newton-Raphson迭代求解：

U(i+1)=Ui+ΔU

(21)

4 數(shù)值模擬

某汽車彈簧盒支座的幾何形狀如圖3所示，材料參數(shù)如下：彈性模量207 GPa，泊松比0.28，屈服應(yīng)力154.31 MPa，冪指數(shù)硬化模量520.4 MPa，硬化參數(shù)0.232，平均各向異性系數(shù)1.653．成形工藝參數(shù)：壓邊力15 kN，摩擦因數(shù)0.15．

圖3 彈簧盒支座零件的CAD模型

圖4給出了擬最小面積法求解的中間構(gòu)型在凸模行程15 mm和25 mm時的比例系數(shù)ξ的分布．其中灰色區(qū)域表示滿足篩選條件式(16)(ξ<0.01或1-ξ<0.01)，處于與凸凹模接觸狀態(tài)的單元．

為評估本文提出的中間構(gòu)型構(gòu)造方法的有效性，圖5展示了由Ls-dyna增量法(INC)和多步成形模擬方法計算的各中間構(gòu)型上截面線AB(見圖6)的Z向坐標(biāo)值比較．

可以看出針對15 mm構(gòu)型時，PMA方法構(gòu)造的中間構(gòu)型在坐標(biāo)范圍-55～-40 mm(20～60 mm)內(nèi)低于、在坐標(biāo)范圍-75～-55 mm(60～80 mm)處高于INC的中間構(gòu)型．觀察此中間構(gòu)型的比例系數(shù)分布圖4(a)發(fā)現(xiàn)，只有曲面A和D上的節(jié)點處于接觸狀態(tài)(灰色區(qū)域)．結(jié)合接觸節(jié)點集合和Delaunay三角網(wǎng)格法，重新構(gòu)造的零件的側(cè)壁與INC方法得到的中間構(gòu)型相吻合，如圖5中的EPMA曲線．

對于25 mm的中間構(gòu)型，PMA方法的截面線僅在坐標(biāo)范圍-55～-40 mm(20～60 mm)內(nèi)低于INC的，觀察其構(gòu)型的比例系數(shù)分布可以發(fā)現(xiàn)，曲面A、C和D均處于接觸區(qū)域，而其他曲面處于非接觸區(qū)域，通過EPMA方法重新生成的側(cè)壁與INC的結(jié)果保持一致．

圖4 中間構(gòu)型比例系數(shù)ξ的分布

圖5 中間構(gòu)型上截面線AB的形狀

(a) 15 mm

(b) 25 mm

圖6 中間構(gòu)型的等效應(yīng)力分布

Fig.6 Distributions of equivalent stress for middle configurations

針對凸模行程35 mm時的中間構(gòu)型，PMA、INC和EPMA 3種方法的中間構(gòu)型趨于一致，這是由于35 mm時，曲面A、B、C和D全部處于接觸狀態(tài)，構(gòu)造的中間構(gòu)型的曲面之間連接的側(cè)壁部分趨同于凹模的側(cè)壁形狀，此時，該3種方法得到的中間構(gòu)型是沒有差別的．

由此可見，EPMA方法可有效改善PMA方法構(gòu)造的中間構(gòu)型的側(cè)壁．圖6給出了兩個中間構(gòu)型的等效應(yīng)力分布．

為表明該多步成形模擬方法在考慮了若干中間構(gòu)型對最終零件應(yīng)力估計精度的改善，圖7展示了由增量法、8個中間構(gòu)型的多步法(MSTP08)和一步法(OSTP)計算的截面線CD(見圖6)的X和Y方向上的應(yīng)力分量分布圖．

從圖7中可以看出，整體趨勢上多步法估計的應(yīng)力分量分布曲線與增量法計算的應(yīng)力分布保持同步的變化，而一步法不能提供經(jīng)歷復(fù)雜變形歷史位置點的有效應(yīng)力估計．在截面線X坐標(biāo)約25 mm處，增量法估計的X應(yīng)力分量存在明顯的局部峰值(約250 MPa)和谷值(約-100 MPa)，Y應(yīng)力分量有最小值(約-350 MPa)；而一步法估計該處的應(yīng)力分量分布均變化較平緩，分別對應(yīng)X應(yīng)力分量的約-15 MPa和10 MPa及Y應(yīng)力分量的-100 MPa；多步法預(yù)示了相應(yīng)處的X應(yīng)力分量局部峰值(約125 MPa)和谷值(約-75 MPa)，及Y應(yīng)力分量的最小值約-200 MPa．類似的對應(yīng)力分量的改善在X坐標(biāo)值約50 mm和200 mm等處觀察到．

可以發(fā)現(xiàn)，采用多步法顯著地改善了最終零件上的拉壓應(yīng)力狀態(tài)的估計精度，如修正了X坐標(biāo)25 mm處的X應(yīng)力分量的局部峰值分布，提高了Y應(yīng)力分量的應(yīng)力精度．

(a) X應(yīng)力分量

(b)Y應(yīng)力分量

圖7 最終構(gòu)型截面線CD上的應(yīng)力分量

Fig.7 Stress components of section lineCDon final configuration

5 結(jié) 論

(1)以擬最小面積法建立的二次規(guī)劃模型來高效生成滑移約束曲面，并根據(jù)設(shè)定閾值篩選接觸節(jié)點和Delaunay三角網(wǎng)格法以改善滑移約束曲面的側(cè)壁使之符合真實變形情況．借助直接映射物質(zhì)點的方法生成初始解構(gòu)型，并采用施加罰函數(shù)約束的方式使構(gòu)造的中間構(gòu)型達(dá)到平衡狀態(tài)．

(2)重新推導(dǎo)了基于中間構(gòu)型的有限應(yīng)變增量度量，以考慮已發(fā)生變形的參考構(gòu)型本身的彎曲效應(yīng)對中性層外物質(zhì)點的影響．

(3)通過對盒支座零件的分析驗證了強(qiáng)化的擬最小面積法生成的中間構(gòu)型質(zhì)量，表明了考慮8個中間構(gòu)型的多步法預(yù)測的應(yīng)力分布的精度，修正了一步法預(yù)估的最終零件上的錯誤的拉壓應(yīng)力和局部極值應(yīng)力分量分布，并得到了與增量法的應(yīng)力估計趨勢相一致的結(jié)果．

[1] GUO Y, BATOZ J, DETRAUX J,etal. Finite element procedures for strain estimations of sheet metal forming parts [J].InternationalJournalforNumericalMethodsinEngineering, 1990,30(8):1385-1401.

[2] MAJLESSI S A, LEE D. Further development of sheet metal forming analysis method [J].JournalofEngineeringforIndustry, 1987,109(4):330-337.

[3] MAJLESSI S A, LEE D. Development of multistage sheet metal forming analysis method [J].JournalofMaterialsShapingTechnology, 1988,6(1):41-54.

[4] LEE C H, HUH H. Three dimensional multi-step inverse analysis for the optimum blank design in sheet metal forming processes [J].JournalofMaterialsProcessingTechnology, 1988,80-81:76-82.

[5] KIM Seungho, KIM Seho, HUH H. Finite element inverse analysis for the design of intermediate dies in multi-stage deep-drawing processes with large aspect ratio [J].JournalofMaterialsProcessingTechnology, 2001,113(1/2/3):779-785.

[6] KIM S H, HUH H. Construction of sliding constraint surfaces and initial guess shapes for intermediate steps in multi-step finite element inverse analysis [J].JournalofMaterialsProcessingTechnology, 2002,130-131:482-489.

[7] GUO Y Q, LI Y M, BOGARD F,etal. An efficient pseudo-inverse approach for damage modeling in the sheet forming process [J].JournalofMaterialsProcessingTechnology, 2004,151(1/2/3):88-97.

[8] HUANG Ying, CHEN Yiping, DU Ruxu. A new approach to solve key issues in multi-step inverse finite-element method in sheet metal stamping [J].InternationalJournalofMechanicalSciences, 2006,48(6):591-600.

[9] LI Y M, ABBES B, GUO Y Q. Two efficient algorithms of plastic integration for sheet forming modeling [J].JournalofManufacturingScienceandEngineering-TransactionsoftheASME, 2007,129(4):698-704.

[10] TANG Bingtao, LI Yunjiang, LU Xiaoyang. Developments of multistep inverse finite element method and its application in formability prediction of multistage sheet metal forming [J].JournalofManufacturingScienceandEngineering-TransactionsoftheASME, 2010,132(4):041013-041019.

[11] HALOUANI A, LI Yuming, ABBES B,etal. Simulation of axi-symmetrical cold forging process by efficient pseudo inverse approach and direct algorithm of plasticity [J].FiniteElementsinAnalysisandDesign, 2012,61:85-96.

[12] 劉偉杰,胡 平,周 平,等. 板材多步逆成形中的滑移約束曲面優(yōu)化構(gòu)造方法[J]. 機(jī)械工程學(xué)報, 2012,48(4):21-25.

LIU Weijie, HU Ping, ZHOU Ping,etal. Research on intermediate configurations of multi-step inverse approach in sheet metal forming [J].JournalofMechanicalEngineering, 2012,48(4):21-25. (in Chinese)