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

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

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

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

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

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

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

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

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

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

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

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

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

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