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

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

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

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

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

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

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

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

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

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

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

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

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

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

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