采用試錯接觸算法的板料成形有限元模擬技術(shù)

1、采用試錯接觸算法的板料成形有限元模擬技術(shù)摘要：以金屬板料成形過程為對象，采用基于板殼理論的8節(jié)點殼單元，以Kirchhoff應(yīng)力張量和Green應(yīng)變張量作為應(yīng)力與應(yīng)變的度量，采用有限變形的Updated Lagrangian列式，研究有限元方法模擬成形過程的技術(shù)。在此基礎(chǔ)上，提出了一種直接試錯的接觸算法，將非線性的接觸邊界條件線性化處理，接觸處理與有限元計算相對獨立。通過試驗結(jié)果與計算結(jié)果的比較，說明采用的模擬算法是合理可行的。    敘詞：金屬板料  成形  有限元  模擬  接觸0 、前言板料沖壓成形是機械工業(yè)

2、中一種重要的加工方法，在航空、宇航、汽車等制造領(lǐng)域均有著廣泛的應(yīng)用。長期以來，保證成形零件的成形質(zhì)量、降低廢品率，一直是板料成形研究的目標(biāo)1，2。從力學(xué)角度而言，板料成形是一個同時涵蓋幾何非線性、材料非線性、邊界非線性的復(fù)雜的力學(xué)過程。以往，分析成形問題多采用基于塑性理論的解析方法，由于成形過程的物理復(fù)雜性，不得不作出較多的簡化和假設(shè)，這就使得解析方法只能進行很簡單的成形問題的分析。近些年發(fā)展起來的板料成形的有限元模擬技術(shù)，使得對復(fù)雜成形問題的模擬分析成為可能，在優(yōu)化成形工藝、提高成形質(zhì)量、降低產(chǎn)品開發(fā)成本等方面將發(fā)揮潛在的重要作用。然而，目前的板料成形模擬技術(shù)仍存在諸多問題亟待解決，因此仍是

3、國際上成形領(lǐng)域的研究熱點35。本文以金屬板料成形過程為對象，研究利用彈塑性有限變形增量有限元方法模擬成形過程的技術(shù)。采用基于板殼理論的8節(jié)點殼單元，以Kirchhoff應(yīng)力張量和Green應(yīng)變張量作為應(yīng)力與應(yīng)變的度量，建立了有限變形的Updated Lagrangian列式。在此基礎(chǔ)上，提出了一種直接試錯的接觸算法，在每一增量步進行接觸搜索，之后交替進行幾何協(xié)調(diào)處理與接觸力的協(xié)調(diào)處理，直至幾何協(xié)調(diào)與接觸力的協(xié)調(diào)同時滿足。通過試驗結(jié)果與計算結(jié)果的比較，說明采用的模擬算法是合理可行的。1、8節(jié)點等參殼單元由于板料實際上是一種板殼結(jié)構(gòu)，因此采用根據(jù)板殼理論建立的殼單元進行板料成形模擬，顯然是合適的。

4、作者在板料成形模擬中采用了一種8節(jié)點等參殼單元，該單元滿足板殼的兩個假設(shè)，其一是中面法線的物質(zhì)線元在板殼變形后仍為直線，其二是忽略與中面垂直的應(yīng)力分量所產(chǎn)生的應(yīng)變能，即認(rèn)為垂直于中面的應(yīng)力分量為零。在該單元中定義了4種不同的坐標(biāo)系：總體坐標(biāo)系(直角坐標(biāo)系)，節(jié)點坐標(biāo)、總剛度矩陣、位移以及外載力矢均定義于該坐標(biāo)系，坐標(biāo)用xi(i=1，2，3)表示，位移用ui表示，ei為3個坐標(biāo)軸方向的單位矢量；隨體的曲線坐標(biāo)系，如圖1所示，和是殼單元中面的兩個曲線坐標(biāo)，是厚度方向的線性坐標(biāo)；中面節(jié)點處的節(jié)點坐標(biāo)系vik，其中k指明此節(jié)點為單元的第k個節(jié)點；最后一種坐標(biāo)系為單元內(nèi)某一點處的局部坐標(biāo)系xi，x1取該

5、點處方向的切矢，由該點處方向的切矢與方向的切矢叉乘得到x3，x2由x3與x1叉乘得到，xi單位矢量化可得到局部坐標(biāo)系的坐標(biāo)基矢量ei。殼單元內(nèi)任意點的位置坐標(biāo)，可以由節(jié)點坐標(biāo)插值得到(1)式中m表示中面n殼單元的節(jié)點數(shù)k節(jié)點k處殼的厚度Nk(，)插值函數(shù)圖1  8節(jié)點殼單元極坐標(biāo)系在總體坐標(biāo)系中，殼單元的位移場由節(jié)點處法線的5個自由度來描述，即法線中點的3個位移umik和法線繞V1k和V2k的2個轉(zhuǎn)角k和k，如圖2所示。這樣，可由單元節(jié)點處的位移插值得到單元的位移場圖2  殼單元的位移模式(2)為了引入厚度方向應(yīng)力分量為零的假設(shè)，單元中的應(yīng)變在局部坐標(biāo)系xi中定義=x，y，

6、xy，yz，xzT(3)2、有限元求解列式物體的變形是連續(xù)的，從初始時刻到t+t時刻之間任一時刻的物體構(gòu)成都可以作為參考構(gòu)形建立有限元列式。實際應(yīng)用中常采用兩種參考構(gòu)形，一種以未變形時的初始構(gòu)形為參考構(gòu)形，稱作Total Lagrangian描述，另一種是以t時刻的平衡構(gòu)形為參考構(gòu)形來描述t+t時刻的變形，稱作Updated Lagrangian描述(或UL方法)。考慮到成形的特點，采用UL有限元列式。2.1  有限變形的有限元列式當(dāng)質(zhì)點在鄰域內(nèi)做剛性運動時，Kirchhoff應(yīng)力張量的各分量在固定于空間的坐標(biāo)系中保持不變。而Green應(yīng)變率張量也是與剛性轉(zhuǎn)動無關(guān)的量，并且Green

7、應(yīng)變ij與Kirchhoff應(yīng)力ij在能量上是共軛的。因此，在本構(gòu)方程中采用Kirchhoff應(yīng)力和Green應(yīng)變作為應(yīng)力和應(yīng)變的合理度量。以t時刻的構(gòu)形為參考構(gòu)形，t+t時刻的虛功方程為(4)式中V，At時刻構(gòu)形的體積區(qū)和受載表面區(qū)p,kt+t時刻定義在t時刻構(gòu)形上的面力b,kt+t時刻定義在t時刻構(gòu)形上的體力變分運算符由于以t時刻的構(gòu)形為參考構(gòu)形，t時刻位移和應(yīng)變均為已知，則式(4)可以寫成VijT(ij+ij)dV=AukTp,kdA+VukTb,kdV(5)(k1+kn)a=-Rs(6)式中k1小位移剛度矩陣kn由于大位移而引出的矩陣a節(jié)點位移增量t+t時刻外載荷的等效節(jié)點力矢量Rs應(yīng)

8、力場的等效節(jié)點力矢量2.2  本構(gòu)方程采用Euler描述時，屈服函數(shù)寫成一般形式為f(ij，E)=(ij)-s(E)=0(7)采用Hill的正交各向異性屈服準(zhǔn)則，并考慮到殼單元的假設(shè)，即3=0，于是得到2=a112+2a1212+a222+a3122+a4232+a5132(8)6個各向異性參數(shù)a1，a12，a2，a3，a4及a5可通過試驗確定。在t時刻，根據(jù)塑性流動法則，可導(dǎo)出本構(gòu)方程ij=Cijkl,tkl(9)上式寫成矩陣形式，則可得到正交異性的彈塑性物理陣(10)式中塑性流矢量彈性物理陣H塑性模量(單向拉伸試驗中應(yīng)力塑性應(yīng)變曲線的切線斜率)3、成形模擬中的接觸算法接觸問題一直

9、是成形模擬中的難點和研究熱點。經(jīng)過十幾年的發(fā)展，借助于有限元方法，產(chǎn)生了眾多用于板料成形模擬的接觸算法。本文中采用的是一種直接試錯的接觸算法，其主要特點是，有限元求解過程與接觸處理過程相對獨立，接觸條件以邊界條件的形式作用于有限元求解方程中。接觸處理過程如圖3所示，模具產(chǎn)生一個小的行程后，首先進行接觸搜索，找到穿透模具表面的板料節(jié)點；然后，進行幾何協(xié)調(diào)處理，根據(jù)這些節(jié)點的當(dāng)前接觸狀態(tài)，以相應(yīng)的方式將它們拉回到模具表面，之后，檢查這些節(jié)點上所受的接觸力是否協(xié)調(diào)，若不協(xié)調(diào)則進行接觸力調(diào)整。這個過程重復(fù)進行，直至同時達到幾何協(xié)調(diào)和接觸力協(xié)調(diào)。以上的接觸求解過程實際上包含兩層迭代，外層是將非線性接觸邊

10、界條件線性化的迭代，內(nèi)層則是根據(jù)給定的邊界條件進行的彈塑性有限變形有限元求解迭代。圖3  接觸處理過程3.1  接觸搜索方法由三角平面片組成的離散網(wǎng)格處理起來相對簡單，并可用于描述復(fù)雜形狀的模具，因此在接觸算法中采用離散網(wǎng)格法描述模具的幾何表面。接觸處理時，認(rèn)為模具是剛性的不發(fā)生變形，只須判別板料節(jié)點是否穿透模具表面，因此接觸搜索稱作單向接觸搜索。模具表面稱作主表面，板料節(jié)點稱作從節(jié)點。接觸狀態(tài)在從節(jié)點處考察，因此要在每個從節(jié)點x(s)處定義間隙函數(shù)g(x(s)=(x(s)-x(m).n(x(m)(11)式中x(m)x(s)在主表面上最近的投影點n主表面上x(m)處的單位外

11、法矢接觸搜索時，在每個從節(jié)點處檢查函數(shù)g的值，若g0則該節(jié)點與主表面是分離的未發(fā)生接觸，而g0則說明該節(jié)點穿透主表面進入了模具內(nèi)部，應(yīng)采取相應(yīng)的幾何協(xié)調(diào)處理。3.2  接觸中的幾何協(xié)調(diào)根據(jù)板料與模具之間的分離、附著和滑動這3種接觸狀態(tài)，將板料節(jié)點分為3類，即分離節(jié)點、附著節(jié)點和滑動節(jié)點。成形過程中，板料節(jié)點的狀態(tài)不斷變化，但必須滿足接觸條件：板料節(jié)點不能穿透模具表面。與模具接觸的板料節(jié)點同模具之間不能作用拉力。圖4  間隙函數(shù)的定義相應(yīng)于板料節(jié)點的分離、附著和滑動3種狀態(tài)，進行不同的幾何協(xié)調(diào)處理：(1)若從節(jié)點x(s)在當(dāng)前幾何協(xié)調(diào)處理前處于分離狀態(tài)，則將其拉回到它在主表面

12、上的投影位置x(m)，并認(rèn)為該節(jié)點當(dāng)前為滑動狀態(tài)。(2)若從節(jié)點x(s)在當(dāng)前幾何協(xié)調(diào)處理前處于附著狀態(tài)，則沿模具行程將其移動一個距離，并認(rèn)為該節(jié)點當(dāng)前仍為附著狀態(tài)。(3)若從節(jié)點x(s)在當(dāng)前幾何協(xié)調(diào)處理前處于滑動狀態(tài)，則將其拉回到它在主表面上的投影位置x(m)，并認(rèn)為該節(jié)點當(dāng)前仍為滑動狀態(tài)。3.3  接觸力協(xié)調(diào)處理由于將從節(jié)點拉回到模具表面的位置是近似的，因此得到的接觸力不一定是協(xié)調(diào)合理的。接觸力(模具作用于板料的力)存在兩種形式的不協(xié)調(diào)：法向接觸力為拉力，在這樣的從節(jié)點處應(yīng)將接觸力釋放掉，并令從節(jié)點成為分離節(jié)點。切向接觸力不協(xié)調(diào)，處于滑動狀態(tài)的從節(jié)點所受切向接觸力小于摩擦力，或

13、處于附著狀態(tài)的從節(jié)點所受切向接觸力大于摩擦力。切向接觸力不協(xié)調(diào)又分為兩種情況：(1)t時刻從節(jié)點的狀態(tài)為附著狀態(tài)，如圖5a所示。t+t時刻，若Fc,tFc，n，則應(yīng)在該從節(jié)點處施加調(diào)整力Fc=-(Fc,t-Fc，n)t(12)式中Fc,tt+t時刻的切向接觸力矢Fc,nt+t時刻的法向接觸力矢tFc,t的方向矢量松弛因子(01)同時，將該節(jié)點由附著狀態(tài)變?yōu)榛瑒訝顟B(tài)。(2)t時刻從節(jié)點的狀態(tài)為滑動狀態(tài)，如圖5b所示。t+t時刻，若Fc.u0則顯然違背了能量守衡，從節(jié)點不能滑動而應(yīng)將狀態(tài)轉(zhuǎn)變?yōu)楦街鵂顟B(tài)，并對從節(jié)點施加調(diào)整力Fc=-Fc(13)圖5  從節(jié)點所受的接觸力的變化4、試驗及模擬

14、算例4.1  方盒拉深零件的拉裂失穩(wěn)分析凹模尺寸為39mm×39mm，凸模尺寸為36.5mm×36.5mm，凸凹模圓角半徑均為5mm。毛料為91mm×91mm的正方形板料，厚度為1.0mm，材料為LY12M。試驗中，凸模行程13mm時零件靠近底部圓角處發(fā)生破裂，如圖6a。通過成形模擬，得到了成形深度為13mm1和2投放到成形極限圖上，如圖7所示，顯然其中某些點已超出了成形極限，這說明該模擬算法能夠用于預(yù)測拉裂失穩(wěn)。圖6  拉深成形的方盒零件圖7  應(yīng)變在FLD上的分布4.2  球頭形拉深件成形模擬凹模直徑為30mm，凸模球頭半徑為13.75mm，凹模圓角半徑為5mm。采用56mm×56mm的正方形毛料，厚度為0.8mm，材料為LY12M，凸模行程為15mm。分別用本文的成形模擬程序和非線形有限元軟件ABAQUS進行了成形模擬，對得到的成形力曲線進行了對比，圖8為成形力的模擬結(jié)果，顯然兩條曲線的趨勢是大體相同的。圖8  模