基于理想形變理論的有限元逆算法在板料成形中的應(yīng)用

朱光成，吳建軍

（西北工業(yè)大學(xué) 機電學(xué)院，陜西 西安 710072）

基于理想形變理論的有限元逆算法在板料成形中的應(yīng)用

朱光成，吳建軍

（西北工業(yè)大學(xué) 機電學(xué)院，陜西 西安 710072）

基于理想形變理論，研究了金屬板料成形的有限元逆算法，并開發(fā)了計算程序。采用線性三角形膜單元和厚向異性的剛塑性材料模型，計算了一個帶凸緣的方盒形件，并與Dynaform的一步逆向法和增量法的計算結(jié)果進行了比較。實例分析結(jié)果表明，在工程允許精度范圍內(nèi)，本文方法能夠有效分析零件的成形性能。由于該方法計算速度快，所以可用于零件的早期設(shè)計。

機械制造；板料成形；理想形變理論；有限元

1 引言

隨著計算機技術(shù)的迅速發(fā)展和有限元理論的日趨成熟，有限元方法在坯料計算及板料成形中的應(yīng)用越來越廣泛。目前，常用的板料成形數(shù)值模擬方法主要有兩種：基于增量理論（又稱流動理論）的有限元增量法和基于全量理論（又稱形變理論）的有限元逆算法。增量法可以全面考慮各種影響因素，是目前比較精確的方法。但在產(chǎn)品及工藝設(shè)計的早期階段，通常只有產(chǎn)品形狀而無模具信息，因而，使用增量法進行數(shù)值模擬比較困難。對設(shè)計人員來講，最重要的是如何根據(jù)產(chǎn)品或已經(jīng)完成工藝補充的沖壓件幾何形狀來快速預(yù)測它的毛料展開形狀和零件的厚度分布等信息，以此來預(yù)測零件的可成形性和工藝的可加工性。因此，K.Chung和O.Richmond將形變理論與極值功路徑相結(jié)合，提出了理想形變的概念[1]。本文將其應(yīng)用于板料成形過程，開發(fā)了有限元逆算法程序。采用線性三角形膜單元和厚向異性的剛塑性材料模型，計算了一個帶凸緣的方盒形件，驗證了本文方法的有效性。

2 理想形變理論[2]

理想形變的基本假設(shè)為：變形沿著最小塑性功路徑進行，材料為剛塑性，滿足R.Hill’79屈服準(zhǔn)則。整個變形過程達(dá)到最小塑性功路徑的條件有兩個：塑性變形的主伸長對應(yīng)于固定的主物質(zhì)線；塑性變形的主自然應(yīng)變比保持不變。

理想形變理論認(rèn)為變形體在整體塑性功取得相對極值的條件下變形較為均勻。變形取得相對極值的條件為：在成形的最后一刻，變形體在邊界上僅受法向力作用而處于平衡狀態(tài)，且變形分布較為均勻，并以變形過程中消耗的整體塑性功最小作為應(yīng)變優(yōu)化的基礎(chǔ)。文獻[1]中已經(jīng)對該假設(shè)做了證明。將變形過程看作一個勻速的過程，考慮到塑性變形體積不可壓縮的假設(shè)，包括不均勻變形的塑性功可表示為：

式中：σ——Cauchy應(yīng)力；

D——應(yīng)變速率；

σˉ——等效應(yīng)力；

εˉ——等效應(yīng)變；

V0——初始體積。

當(dāng)物質(zhì)微元的變形路徑確定之后，塑性功就只由位移確定，可看作是位移的函數(shù)W=W（U）。

3 理想形變理論的有限元表達(dá)及實現(xiàn)

3.1 變形分析[3][4]

在Lagrange描述法中，物質(zhì)坐標(biāo)系在初始構(gòu)形中與整體坐標(biāo)系重合。變形前，物質(zhì)坐標(biāo)系為直角坐標(biāo)系，各坐標(biāo)軸之間的夾角為直角，且基矢量為單位矢量。變形后，物質(zhì)坐標(biāo)系不再是直角坐標(biāo)系。根據(jù)Kirchhoff假設(shè)，變形前的中面法線在變形后仍然垂直于中面。用P0和Pt分別表示物質(zhì)點P在變形前后的空間位置，我們可以得到變形后物質(zhì)坐標(biāo)系在P點的基矢量t→1、t→2和t→3的表達(dá)式為：

由于已知毛坯形狀為平板，所以單元的右Cauchy－Green變形張量[C]可表示為：

[C]的特征值為單元面內(nèi)兩個主伸長λ1、λ2的平方，則有：

根據(jù)塑性變形中體積不可壓縮的假設(shè)，可得板料上各點的厚向伸長量λ3為：

右Cauchy－Green變形張量[C]也可由坐標(biāo)旋轉(zhuǎn)矩陣[M]與主伸長λi（i=1，2，3）表示為：

其中，坐標(biāo)旋轉(zhuǎn)矩陣[M]為：

θ是λi與最終構(gòu)形的局部坐標(biāo)系x軸的夾角，求解（6）式得：

由右Cauchy－Green變形張量[C]容易得到對數(shù)主應(yīng)變表達(dá)式為：

那么單元沿著厚度方向任意一點的大變形對數(shù)應(yīng)變表達(dá)式為：

3.2 塑性功及其極值[5]

將理想形變應(yīng)用于板料成形設(shè)計時，考慮到零件的最終構(gòu)形和坯料的初始輪廓（一般是平面）是已知的，成形過程的邊界條件未指定，最終構(gòu)形上的物質(zhì)點必須滿足幾何約束條件x3=x3（x1，x2）,再根據(jù)U= x-X（X，x分別是物質(zhì)點在初始構(gòu)形和最終構(gòu)形上的節(jié)點坐標(biāo)），取得整體塑性功極值的靜力平衡方程：

采用Newton－Raphson迭代法，構(gòu)造方程組，求{R（X）}=0，有：

定義切線剛度矩陣為：

求解方程組式（13），就可得到初始構(gòu)形上節(jié)點坐標(biāo)的調(diào)整量{d Xi}，調(diào)整后的節(jié)點坐標(biāo)為：

式中：β——減速因子，取值一般為0～1，本文取為1。

以上各式的具體表達(dá)形式可以非常容易的求得，故這里不再詳細(xì)介紹。

3.3 程序流程圖及毛料初始解的獲取

在Visual C++6.0環(huán)境下開發(fā)有限元逆算法程序，程序的計算流程圖如圖1所示。

有限元逆算法需要一個毛料的初始猜測值。一個合理的毛料初始解將直接影響到求解的收斂性和準(zhǔn)確性。本文采用幾何映射法[6]獲得該初始猜測值。

3.4 收斂準(zhǔn)則及收斂性的提高[7][8]

采用Newton－Raphson迭代法必然存在兩次迭代之間滿足什么準(zhǔn)則可以終止本次迭代的問題。根據(jù)有限元逆算法的實際計算情況，本文選擇位移準(zhǔn)則來進行收斂性判斷，計算節(jié)點位移向量的范數(shù)：

式中：n——節(jié)點總數(shù)；

ρ——給定的收斂判斷因子，對位移收斂準(zhǔn)則ρ一般取10－3。

當(dāng)網(wǎng)格中出現(xiàn)畸變單元或單元大小突變比較嚴(yán)重時，可能使切線剛度矩陣呈現(xiàn)病態(tài)，程序算法無法收斂。為了提高逆算法的網(wǎng)格處理能力，添加阻尼矩陣來改善切線剛度矩陣的病態(tài)性。

由式（13）和式（14）可得：

圖1 程序計算流程圖

在切線剛度矩陣中加入阻尼矩陣δi[I]，[I]為單位矩陣，δi[I]為對角矩陣，則方程組的解為：

在本文程序算法中取：

式中：Kjj——第i迭代步中切線剛度矩陣第j行第j列的對角線元素值。

4 有限元逆算法應(yīng)用實例

圖2 劃分有限元網(wǎng)格的產(chǎn)品最終構(gòu)形

以尺寸為50mm×50mm×20mm、各處圓角半徑均為6mm的帶凸緣的方盒形件為例進行計算。圖2是劃分三角形單元后產(chǎn)品的最終構(gòu)形。單元尺寸為2.5mm，盒形件劃分為2384個單元和1237個節(jié)點。材料的應(yīng)力應(yīng)變曲線σˉ=520.4εˉ0.232MPa，厚向異性指數(shù)為r= 1.65，毛料的初始厚度為t=1.0mm。采用本文方法、Dynaform的一步逆向法和增量法的計算結(jié)果見3所示，三種方法計算的零件厚度的云圖分布規(guī)律基本上是一致的，只存在個別區(qū)域的差異。本文方法與Dynaform一步逆向法相比：增厚最大誤差為2.60%，減薄最大誤差為 0.82%；本文方法與Dynaform增量法相比：增厚最大誤差為12.84%，減薄最大誤差為6.77%。采用Dynaform增量法用時1976s，而本文程序迭代6次就收斂，用時僅493s，本文方法明顯比增量法快很多。本例計算是在Windows XP系統(tǒng)上完成的，CPU是AMD Processor 5050e/2.61GHz，內(nèi)存為2GB。

在板料成形中，零件的厚度變化是一重要的因素，因為變薄嚴(yán)重會有破裂的可能，增厚嚴(yán)重會有起皺的可能，所以算例中主要對零件的厚度進行了分析和比較。

5 結(jié)論

根據(jù)理想形變理論結(jié)合板料成形的具體情況，給出了有限元逆算法的數(shù)學(xué)公式和有限元表達(dá)，開發(fā)了工程分析模塊。同Dynaform的一步逆向法和增量法的計算結(jié)果進行對比，本文方法計算的零件厚度分布情況同實際情況更加吻合。由于本文方法在計算過程中簡化了工藝參數(shù)對板料成形的影響，這樣的結(jié)果在工程分析中是可以接受的。通過實例可以證明，有限元逆算法能夠較好的滿足設(shè)計的需要，是一種高效的數(shù)值模擬工具。雖然逆算法模擬精度相對于增量法低了一些，但隨著算法的不斷完善，其模擬精度還可以得到進一步的提高。

圖3 零件的厚度分布

[1]K.Chung，O.Richmond.Idealforming-Ⅱ.Sheetforming withoptimum deformation[J].Int.J.Mech.Sci.，1992，34（8）:617-633.

[2] 蘭 箭，董湘懷，李志剛.用有限元逆算法計算板料成形毛坯形狀和應(yīng)變分布[J].塑性工程學(xué)報，2001，8（2）:60-62.

[3] 吳建軍，陳衛(wèi)彬，李順平.復(fù)雜形狀拉深件快速展開與成形模擬[J].機械科學(xué)與技術(shù)，2004，23（3）:370-375.

[4] 陳衛(wèi)彬.基于反向模擬法的拉深成形性分析技術(shù)研究[D].西安:西北工業(yè)大學(xué)出版社，2003.

[5] K.Chung，S.Y.Lee，F(xiàn).Barlat.Finite elementsimulation ofsheetforming based on a planaranisotropic strain-rate potential[J].InternationalJournalofPlasticity，1996，12（1）:93-115.

[6] 吳建軍，楊漢平.考慮外法向的鈑金零件展開方法研究[J].中國機械工程，2006，17（15）:1546-1549.

[7] Y.Q.Guo，J.L.Batoz，H.Naceur.Recent developments on the analysis and optimum design of sheet metal forming parts using a simplified inverse approach[J].Computers and Structures，2000，78:133-148.

[8] 徐國艷，高 峰，杜發(fā)榮，張立玲.基于快速有限元分析的沖壓件毛料展開[J].機械科學(xué)與技術(shù)，2006，25（4）:386-389.

Application of finite element inverse approach in sheet metal forming process based on ideal forming theory

ZHU Guangcheng，WU Jianjun
（Northwestern Polytechnical University, School of Mechatronics, Xi'an 710072，Shanxi China）

Based on ideal forming theory, the finite element inverse approach for sheet metal forming has been explored, and the computer program has been implemented. By use of the linear triangular membrane element and the thickness anisotropy rigid-plasticity material model，a flanged square box has been analyzed and its result has been compared with the results obtained from MSTEP and incremental approach of Dynaform software. The results show that the finite element inverse approach is an effective approach for sheet metal forming only if its error is allowable during the engineering calculation. Because of its rapid computing speed, this method can be used at the early design stage of sheet metal forming parts.

Sheet metal forming；Ideal forming theory；Finite element

TG386

A

1672－0121（2011）03－0049－03

國家863計劃資助項目（2008AA04Z120）

2011－02－22

朱光成（1985－），男，碩士在讀，主攻板料成形數(shù)值模擬技術(shù)