劉永財(cái)，陳文亮，胡慶婉，鮑益東

(1. 南京航空航天大學(xué)機(jī)電學(xué)院，江蘇省精密與微細(xì)制造技術(shù)重點(diǎn)實(shí)驗(yàn)室，南京 210016； 2. 曲靖師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，曲靖 655011)

在板料成形有限元模擬領(lǐng)域，基于全量理論的一步成形有限元法[1]僅考慮初始構(gòu)形和最終構(gòu)形，計(jì)算效率高，能夠快速計(jì)算板料的初始厚度以及應(yīng)變等物理量分布，得到了廣泛的應(yīng)用。但是，一步成形有限元法也帶來最終零件的應(yīng)力精度估計(jì)不足的問題。為彌補(bǔ)這一不足，人們通過引入中間構(gòu)形的方式來表征加載路徑和變形歷史，提出了多步成形有限元法[2－3]，也稱為偽逆法(pseudo inverse approach)[4－5]，該方法被廣泛應(yīng)用于板料沖壓成形領(lǐng)域[6]。在多步成形有限元法中，如何快速準(zhǔn)確的求解出中間構(gòu)形是其中的一個(gè)關(guān)鍵問題?；谔撐灰圃韀7]，采用彈塑性有限元模型，中間構(gòu)形可以看作是一個(gè)關(guān)于節(jié)點(diǎn)位移的非線性平衡問題，這類問題通常采用Newton-Raphson 迭代法求解。由于迭代法的收斂性強(qiáng)烈依賴于初值與真實(shí)值的逼近程度，因此，一個(gè)快速準(zhǔn)確的中間構(gòu)形初始解預(yù)示算法，對(duì)于有限元多步法的收斂效率，具有重要的作用。

中間構(gòu)形初始解的預(yù)示算法得到了眾多學(xué)者廣泛的研究。Lee 等[8]首先提出了采用映射函數(shù)的反向映射法，該方法可以處理傾角較大的零件，然而對(duì)于復(fù)雜的成形零件卻難以找到對(duì)應(yīng)的映射函數(shù)。Kim 等[9]提出了網(wǎng)格映射法，這種方法需要將零件進(jìn)行分區(qū)處理，算法實(shí)現(xiàn)復(fù)雜，而且分區(qū)對(duì)最終結(jié)果影響較大，穩(wěn)定性較差。在截面線法的基礎(chǔ)上，吳建軍等[10]提出了雙展開平面映射法，該方法使用截面線法對(duì)兩個(gè)空間形狀進(jìn)行平面展開，簡(jiǎn)化了計(jì)算量，但是對(duì)于局部網(wǎng)格畸變程度大等情況，計(jì)算結(jié)果仍不理想。Guo 等[11]基于最小面積法的思想，使用序列二次規(guī)劃(SQP)算法求解一系列約束條件下的最小弧長(zhǎng)問題。Tang 等[12]將最小面積法擴(kuò)展到方盒件等形狀規(guī)則零件上的初始解預(yù)示應(yīng)用上。然而，最小面積法程序?qū)崿F(xiàn)繁瑣，通用性較差，對(duì)于復(fù)雜的成形零件難以得到理想的結(jié)果。張向奎等[13－14]以三角形單元面積平方和最小為目標(biāo)函數(shù)，提出了偽最小面積法，基于板料變形過程中單元節(jié)點(diǎn)不發(fā)生橫向位移的假設(shè)，將問題轉(zhuǎn)化為求解一系列的非線性最小二乘問題，但是該方法沒有考慮橫向位移外其他方向上的位移。Fu 等[15]通過網(wǎng)格單元位置信息，基于網(wǎng)格參數(shù)平面化的思想實(shí)現(xiàn)了柱坐標(biāo)下的初始解預(yù)示，并將其用于管材的有限元成形模擬中。陳偉等[16]研究了多步正成形算法中的中間構(gòu)形初始解預(yù)示算法，根據(jù)初始平板料和工具曲面展平板料之間的映射對(duì)應(yīng)關(guān)系，在兩個(gè)空間形狀間進(jìn)行插值計(jì)算，結(jié)合接觸判斷和修正算法，得到中間構(gòu)形的初始形狀，然而初始形狀強(qiáng)烈依賴于工具曲面的展平板料形狀。鮑益東等[17]使用張拉膜單元方法計(jì)算中間構(gòu)形初始解，并將其成功集成到自主開發(fā)的板料快速模擬成形QuickForm 軟件中。然而，該方法需要求解一系列復(fù)雜的非線性方程組，計(jì)算效率較低。

本文將利用解耦格式對(duì)中間構(gòu)形進(jìn)行分解計(jì)算，推導(dǎo)節(jié)點(diǎn)位移約束條件，首次利用構(gòu)造線性Laplace-Beltrami(LB)方程模型來計(jì)算中間構(gòu)形的初始解，該方法僅需少量的線性迭代就可以快速地生成非對(duì)稱復(fù)雜零件的性態(tài)良好初始解。數(shù)值算例表明，采用LB 方程的方法能快速準(zhǔn)確地生成中間構(gòu)形的初始解，加快迭代法的收斂速度，減少多步有限元法的運(yùn)行時(shí)間。

1 解耦格式

在多步有限元法中，將中間構(gòu)形分解為彎曲變形和拉伸變形的解耦格式，采用該格式求解能夠增加步長(zhǎng)以及加快收斂速度。

根據(jù)虛位移原理，在任一時(shí)刻，有：

式中：ε*是虛應(yīng)變；σ是應(yīng)力； Δu*是虛位移增量；f是外力；V是體積。注意到虛位移的任意性以及應(yīng)變與位移的幾何關(guān)系和應(yīng)力與應(yīng)變的本構(gòu)關(guān)系，則在t時(shí)刻下，可以得到的有限元系統(tǒng)的平衡條件為：

式中：tF表示t時(shí)刻外部作用節(jié)點(diǎn)力；tR表示t時(shí)刻單元應(yīng)力的節(jié)點(diǎn)力。

在t+Δt時(shí)刻下，有平衡條件：t+ΔtF-t+ΔtR=0，假設(shè)tR已知，令R=t+ΔtR-tR，并將該向量近似為R=tK·ΔU，則有：

將剛度矩陣Kt看成由膜單元?jiǎng)偠染仃嚭蛷澢鷦偠染仃噧刹糠纸M成的，因此，如圖1 所示，在節(jié)點(diǎn)所在的局部坐標(biāo)系中，可以建立如下方程組來求解時(shí)刻tt+Δ 的中間構(gòu)形[18]：

式中：1Δu和 2Δu為待求解節(jié)點(diǎn)的位移增量向量；0u為指定節(jié)點(diǎn)的位移向量；t t+ΔF為tt+Δ 時(shí)刻外力在切平面上引起的外部作用的摩擦力向量；tR1和tR2分別為t時(shí)刻指定節(jié)點(diǎn)的位移在法線方向上和切平面上所產(chǎn)生的力向量；tKmij和tKgij(i,j=1,2)分別為t時(shí)刻下的膜單元?jiǎng)偠染仃嚭蛷澢鷨卧獎(jiǎng)偠染仃?。利用坐?biāo)系之間的轉(zhuǎn)換關(guān)系，對(duì)式(4)進(jìn)行坐標(biāo)轉(zhuǎn)換。首先，將方程組從節(jié)點(diǎn)所在的局部坐標(biāo)系轉(zhuǎn)換到單元局部坐標(biāo)系中，再從單元局部坐標(biāo)系轉(zhuǎn)換到全局坐標(biāo)系中，從而得到全局坐標(biāo)系中的整體有限元方程組。在轉(zhuǎn)換的過程中，節(jié)點(diǎn)局部坐標(biāo)系下的單元?jiǎng)偠染仃嚤晦D(zhuǎn)換為全局坐標(biāo)系下的單元?jiǎng)偠染仃?，而膜單元矩陣中的元素值遠(yuǎn)遠(yuǎn)大于彎曲單元矩陣中的元素值，這將導(dǎo)致方程組高度病態(tài)，造成迭代無法快速收斂甚至迭代發(fā)散。

圖1 單元坐標(biāo)系統(tǒng) Fig.1 Element coordinate systems

為了解決上述問題，本文將薄板的中間構(gòu)形計(jì)算分解為彎曲變形和拉伸變形兩個(gè)獨(dú)立的解耦計(jì)算過程。它假設(shè)節(jié)點(diǎn)法向與其鄰接的有限元網(wǎng)格法向之間的夾角均可以忽略不計(jì)，那么，在法線方向和切平面上，式(4)將被解耦分解為式(5)和式(6)：

式(5)對(duì)應(yīng)薄板的彎曲變形，式(6)對(duì)應(yīng)薄板的拉伸 變形。

為了精確滿足tt+Δ 時(shí)刻下有限元系統(tǒng)的平衡條件，需要使用Newton-Rapshon 法[19]對(duì)式(5)進(jìn)行增量迭代求解，即：

式中，i= 1,2,3,…，初始條件為：

對(duì)式(6)的Newton-Rapshon 迭代法求解可以類似 得到。

可以看出，解耦格式在計(jì)算上能夠帶來兩方面優(yōu)勢(shì)：它將一個(gè)大型方程組分解為兩個(gè)小型方程組；兩個(gè)小型方程組單元?jiǎng)偠染仃囍械脑刂迪嗖畈淮?，條件數(shù)得到有效的改善，因此可以加快迭代算法的收斂速度。

綜上所述，中間構(gòu)形被分解為彎曲變形和拉伸變形兩個(gè)計(jì)算過程。其中，以彎曲變形為初始解，通過彈塑性流動(dòng)力學(xué)進(jìn)行迭代求解來計(jì)算拉伸變形，該結(jié)果為一個(gè)中間構(gòu)形。本文只研究計(jì)算彎曲變形，稱其為中間構(gòu)形初始解?？梢钥闯?，一個(gè)高效的中間構(gòu)形初始解構(gòu)造方法十分必要。由于Laplace-Beltrami 方程具有堅(jiān)實(shí)的數(shù)學(xué)理論和明確的物理意義，它已經(jīng)被成功應(yīng)用到數(shù)字幾何處理領(lǐng)域，受到了廣泛的關(guān)注和研究。下面將研究使用LB方程來構(gòu)造中間構(gòu)形初始解。

2 Laplace-Beltrami 方程

考慮一個(gè)嵌入三維空間的二維曲面，它定義了一個(gè)黎曼流形空間。假設(shè)g=(gij)2×2該流形空間上的度量張量，則標(biāo)量函數(shù)f的LB 算子[20]可以表示為：

式中，G=det(g),(gij) =(gij)-1，而gij是兩個(gè)協(xié)變張量基做內(nèi)積的結(jié)果。然而，對(duì)于一般的復(fù)雜曲面，則很難得到連續(xù)情形下LB 算子的顯式表達(dá)式。因此，在實(shí)際應(yīng)用中，通常將待處理的復(fù)雜曲面表示為離散的多邊形網(wǎng)格形式，并在網(wǎng)格上構(gòu)造離散LB算子的顯式表達(dá)式[21]。

這里以三角形網(wǎng)格為例來說明離散LB 算子表達(dá)式的構(gòu)造。首先考慮在某一頂點(diǎn)vi處的一階離散LB 算 子ci=-x(i+1)mod3+x(i+2)mod3,(i= 1,2,3)的 表達(dá)式，寫成如下形式：

式中：M表示三角形網(wǎng)格集合；f表示定義在M的函數(shù)；MΔ 表示M上LB 算子Δ 的離散逼近；n為M中頂點(diǎn)的個(gè)數(shù)；ijw為組合系數(shù)。從式(9)可以看出，頂點(diǎn)iv處的離散LB 算子可以表示所有頂點(diǎn)函數(shù)值的線性組合。

本文選取基于平均曲率計(jì)算理論推導(dǎo)的LB 算子的余切表達(dá)式[22]為：

圖2 Voronoi 區(qū)域和角度 Fig.2 Voronoi area and two angles opposite to an edge vivj

可以證明[23]，當(dāng)頂點(diǎn)存在6 個(gè)相鄰頂點(diǎn)且函數(shù)f充分光滑的條件下，如果使用上述式(10)的離散LB 算子ΔM，則它能夠以O(shè)(h2)的速度收斂于連續(xù)LB 算子Δ。其中，h表示三角形網(wǎng)格的尺寸，如圖2 所示。

接著，聯(lián)立所有頂點(diǎn)的一階離散LB 算子表達(dá)式形成一個(gè)線性方程組： 式中，W=diag(A1-mixed,A2-mixed, …,An-mixed)；b為右端項(xiàng)，LM為對(duì)稱系數(shù)矩陣；當(dāng)且僅當(dāng)頂點(diǎn)vi和頂點(diǎn)vj相鄰時(shí)，(LM)ij具有非零值，即：

類似地，可以推導(dǎo)出二階LB 方程為：

一階LB 算子對(duì)應(yīng)的是曲面第一基本形式的變化量，而二階LB 算子對(duì)應(yīng)的是曲面第二基本形式的變化量，為了綜合考慮這兩種形式的變化量。一般地，將一階和二階LB 算子相結(jié)合獲得LB 方程來表達(dá)網(wǎng)格曲面變形，即：

式中，參數(shù)1k和2k分別取為1 和10[24]。

LB 方程可以求解約束條件下的網(wǎng)格曲面變形問題。在求解過程中，取f(vi)依次取為節(jié)點(diǎn)vi處的空間位置坐標(biāo)時(shí)，則方程組的解為節(jié)點(diǎn)更新后的空間位置坐標(biāo)。此外，為確定出方程的邊界條件，需要將網(wǎng)格區(qū)域分為固定區(qū)域、約束區(qū)域和自由區(qū)域三部分。其中，固定區(qū)域可以通過網(wǎng)格與壓邊圈或凹模的接觸判斷來確定，而約束區(qū)域和自由區(qū)域?qū)⑼ㄟ^大位移小應(yīng)變理論來確定。

3 約束條件及方程求解

在薄板彎曲變形的過程中，盡管結(jié)構(gòu)的位移很大，但是結(jié)構(gòu)的應(yīng)變卻較小，因此它可以視為一個(gè)大位移小應(yīng)變問題[25]。此時(shí)，應(yīng)力和應(yīng)變之間具有線性關(guān)系，滿足式(15)：

式中：D為彈性矩陣；0σ為初應(yīng)力，而應(yīng)變和位移之間具有非線性關(guān)系，即：

式中：Lε和NLε分別為應(yīng)變的線性部分和非線性部分；a為節(jié)點(diǎn)的位移向量。

下面將給出LB和NLB的計(jì)算方法。在全局坐標(biāo)系下，上一個(gè)臨時(shí)構(gòu)形Ω和當(dāng)前臨時(shí)構(gòu)形Ω′上的變形前后兩個(gè)三角形分別記為ABC和A B C′ ′ ′，將它們分別映射到局部坐標(biāo)系下，并記為abc和a b c′ ′ ′，如圖3 所示。在局部坐標(biāo)系下，三角形單元變形可以視為一個(gè)平面問題，可得矩陣LB的表達(dá)式：

式 中，bi=y(i+1)mod3-y(i+2)mod3和ci=x(i+1)mod3-x(i+2)mod3,(i= 1,2,3)為局部坐標(biāo)系下三角形abc的三個(gè)節(jié)點(diǎn)的平面位置坐標(biāo)。根據(jù)變分原理，可推導(dǎo)出矩陣BNL的表達(dá)式[26]：

式中，矩陣A中的元素通過式(19)得：

式中，a是三角形三個(gè)節(jié)點(diǎn)的位移向量。

圖3 不同坐標(biāo)系下的三角形變形描述 Fig.3 Triangle mesh deformation in different coordinate systems

得到LB和NLB后，就可以通過式(16)來計(jì)算應(yīng)變。接著，通過式(15)計(jì)算出每個(gè)單元的應(yīng)力，進(jìn)而得到每個(gè)單元所受的內(nèi)力。根據(jù)靜力等效原則，通過單元內(nèi)力能夠計(jì)算出每個(gè)節(jié)點(diǎn)處的等效節(jié)點(diǎn)力?；诘刃Ч?jié)點(diǎn)力，可通過雙動(dòng)情形來建立LB方程的約束條件。

遍歷所有與凸模曲面網(wǎng)格發(fā)生接觸的單元節(jié)點(diǎn)，根據(jù)等效節(jié)點(diǎn)力的方向和法線方向之間夾角大小來判斷節(jié)點(diǎn)是否為約束節(jié)點(diǎn)。

1) 如果兩者之間的夾角小于90°，力將節(jié)點(diǎn)拉離接觸面，當(dāng)前構(gòu)形的約束節(jié)點(diǎn)在下一個(gè)臨時(shí)構(gòu)形中成為自由節(jié)點(diǎn)，則需要解除該節(jié)點(diǎn)的固定約束條件，如圖4 中的A點(diǎn)所示；

2) 如果兩者之間的夾角大于90°，力將節(jié)點(diǎn)繼續(xù)固定在接觸面上，則該節(jié)點(diǎn)在下一個(gè)臨時(shí)構(gòu)形中仍為約束節(jié)點(diǎn)，如圖4 中的B點(diǎn)所示。

可以看出，隨著迭代次數(shù)的增加，自由節(jié)點(diǎn)的個(gè)數(shù)逐漸增加，約束節(jié)點(diǎn)的個(gè)數(shù)逐漸減少，直至臨時(shí)構(gòu)形滿足相應(yīng)的平衡條件。

圖4 約束節(jié)點(diǎn)和自由節(jié)點(diǎn)示意圖 Fig.4 Constraint and free nodes diagram

從式(13)可以看出，矩陣W-1(LMW-1LM)不是對(duì)稱矩陣。為解決這個(gè)問題，在方程的兩側(cè)分別左乘W，可得式(20)：

由于LM為對(duì)稱矩陣，可知式(19)中的系數(shù)矩陣LMW-1LM為對(duì)稱正定矩陣。另外，從LM的稀疏性可知，矩陣LMW-1LM為稀疏矩陣，因此系數(shù)矩陣是一個(gè)稀疏對(duì)稱正定矩陣?？刹捎肔U 分解法、多重網(wǎng)格法或共軛梯度法等方法[27]對(duì)線性方程組進(jìn)行數(shù)值求解。而當(dāng)前時(shí)間步的中間構(gòu)形初始解是以上一個(gè)時(shí)間步中間構(gòu)形為基礎(chǔ)，通過一系列線性方程組迭代求解得到的。在整個(gè)迭代過程中，由于上一個(gè)時(shí)間步的中間構(gòu)形沒有發(fā)生改變，因而迭代涉及的系數(shù)矩陣保持不變。因此，當(dāng)使用LU 分解法時(shí)，系數(shù)矩陣可被分解為上三角矩陣和下三角矩陣的乘積，所以，在迭代過程中僅需要使用追趕法就能準(zhǔn)確求解。而對(duì)于共軛梯度法，當(dāng)方程組右端項(xiàng)發(fā)生改變的時(shí)候，它沒有利用系數(shù)矩陣不變的特性，需要重新迭代，增加了計(jì)算時(shí)間。多重網(wǎng)格法雖然計(jì)算時(shí)間較短，但是程序?qū)崿F(xiàn)較為困難。因此，結(jié)合算法的實(shí)現(xiàn)難易和效率好壞綜合考慮，本文采用LU 分解法來求解LB 方程。

4 算例

基于解耦算法的多步有限元算法，已經(jīng)被成功集成到自主開發(fā)的板料成形模擬軟件QuickForm中。其中，使用LB 方程來構(gòu)造中間構(gòu)形初始解。為驗(yàn)證算法的有效性，這里以板料成形商業(yè)有限元軟件DynaForm 的模擬結(jié)果作為標(biāo)準(zhǔn)值，該軟件因使用動(dòng)力顯式算法，因?yàn)槠淠M精度高、可靠性好，已經(jīng)被廣泛應(yīng)用到覆蓋件的板料成形的仿真分析中。

以NUMISHEET’1994 中某汽車翼子板標(biāo)準(zhǔn)零件的沖壓成形為例，來說明使用LB 方程中間初始解算法的有效性。該零件凹模的CAD 模型如圖5所示。初始坯料的尺寸為1920 mm×1530 mm，厚度為1.0 mm，材料應(yīng)力-應(yīng)變關(guān)系為σ= 648(ε+ 0.004)0.220MPa，彈性模量為E= 2.07GPa，泊松比為ν= 0.28，沖頭力的大小為Fpunch= 280 kN ，摩擦力系數(shù)為 0.12μ= 。坯料被劃分為11,737 個(gè)節(jié)點(diǎn)，23,040 個(gè)三角形單元網(wǎng)格，拉延深度為130.40 mm。使用測(cè)試的計(jì)算機(jī)處理器為 Inter(R)Core(TM) i5-7300，內(nèi)存為2.60 GHz。

圖5 零件凹模的CAD 模型 Fig.5 CAD model of a die part

圖6(a)為拉延深度為90 mm 時(shí)，采用LB 方程迭代6 次計(jì)算得到的初始中間構(gòu)形的離散三角形網(wǎng)格。此時(shí)，板料被劃分為43031 個(gè)節(jié)點(diǎn)、79923 個(gè)單元，相較于原始平板料的網(wǎng)格數(shù)有所增加，這是為了更好地表達(dá)曲率變化較大處零件特征而對(duì)三角形網(wǎng)格進(jìn)行加密處理的緣故，它的總面積為2945308 mm2?？梢钥闯?，算法構(gòu)造的臨時(shí)構(gòu)形網(wǎng)格質(zhì)量良好，沒有出現(xiàn)網(wǎng)格重疊等缺陷。圖6(b)給出了迭代1 次，迭代2 次以及迭代3 次的臨時(shí)中間形狀沿截面線截取的部分曲線形狀結(jié)果對(duì)比。從圖 6 可以看出，隨著迭代次數(shù)的增加，曲線上的約束節(jié)點(diǎn)個(gè)數(shù)越來越少，而自由變形的區(qū)域越來越大，曲線曲面的光滑性越來越好。

另外，圖6 中還分別給出了使用張拉膜單元方法和LB 算子方法沿著截面線得到的中間構(gòu)形初始解形狀?？梢钥闯?，兩種方法均可以獲得良好的中間構(gòu)形初始解形狀，它們得到的結(jié)果非常接近。但是，采用相同的迭代終止條件時(shí)，兩種方法在計(jì)算當(dāng)前步的中間構(gòu)形初始解所需時(shí)間卻是不同的，使用張拉膜單元方法計(jì)算時(shí)間為24.32 s，而使用LB算子方法計(jì)算時(shí)間為16.68 s。顯然，LB 算子方法的計(jì)算速度更快，計(jì)算效率更高，這是因?yàn)樗蠼獾木€性方程組維數(shù)更少的緣故。

圖6 初始中間構(gòu)形和沿截面線的曲線形狀 Fig.6 Intermediate configurations

圖7 分別給出了使用DynaForm、QuickForm 和AutoForm 三種軟件得到的合模狀態(tài)下的零件厚度分布圖，其中板料成形商業(yè)有限元軟件AutoForm占有率高，它使用靜力隱式算法，具有模擬速度快，計(jì)算準(zhǔn)確度高的特點(diǎn)。從圖7(a)～圖7(c)對(duì)比中可以看出，三者的厚度分布具有一致的趨勢(shì)。其中，最大減薄部分均位于零件的相同區(qū)域，見圖中標(biāo)記部分，它們分別為0.70 mm、0.79 mm 和0.82 mm，驗(yàn)證了多步法算法的有效性。

圖7 DynaForm、QuickForm 和AutoForm 三種軟件模擬的板料厚度分布圖 Fig.7 Thickness distributions by DynaForm, QuickForm and AutoForm

進(jìn)一步來分析動(dòng)力顯式算法和靜力隱式算法的模擬結(jié)果和運(yùn)行時(shí)間。合模狀態(tài)下零件在XOY平面上的邊界輪廓線如圖8 所示，三者得到的邊界線形狀基本一致。其中QuickForm 和DynaForm 兩者之間的最大距離為58.5 mm，和零件的總尺寸相比，這里產(chǎn)生的相對(duì)誤差小于3.12%。三者計(jì)算的單元總面積分別為2998960 mm2、3000394 mm2和2954358 mm2，可見靜力隱式算法和動(dòng)力顯式算法結(jié)果之間的相對(duì)誤差小于1.5%。然而，在相同的計(jì)算機(jī)配置下，使用基于動(dòng)力顯式算法的DynaForm軟件運(yùn)行的時(shí)間為6899 s，而使用基于靜力隱式算法 AutoForm 軟件和 QuickForm 軟件分別需要1203 s 和1086 s，可見隱式多步成形有限元法的效率大幅提升，大約僅為顯式有限元法的16%。進(jìn)一步，分析基于靜力隱式算法的 QuickForm 和AutoForm 的模擬結(jié)果和運(yùn)行時(shí)間。如圖8 所示，合模狀態(tài)下的兩者計(jì)算模擬的邊界輪廓線較為接近，而QuickForm 軟件較AutoForm 軟件的運(yùn)行時(shí)間減少了10.77%，說明本文提出的構(gòu)造中間解預(yù)示算法是快速且有效的。

以汽車某復(fù)雜零件的沖壓成形為例，來進(jìn)一步說明本文方法有效性。該零件凸模的CAD 零件圖如圖 9 所示。初始坯料的尺寸為 1680 mm× 1270 mm，厚度為1.0 mm，材料具體參數(shù)為：應(yīng)力應(yīng)變關(guān)系為σ= 784(ε+ 0.009)0.21MPa，彈性模量為 2.07E= GPa，泊松比為 0.28ν= ，摩擦力系數(shù)為。坯料被劃分為15 328 個(gè)節(jié)點(diǎn)，29 436 個(gè)三角形單元網(wǎng)格，拉延深度為85 mm。

圖8 最終構(gòu)形的外輪廓圖 Fig.8 Contours of final part

圖9 零件凸模的CAD 模型 Fig.9 CAD model of a punch part

合模狀態(tài)下零件在XOY平面上的邊界輪廓線如圖10 所示，三者得到的邊界線形狀基本一致。而使用基于動(dòng)力顯式算法的DynaForm 軟件運(yùn)行的時(shí)間為4 126 s，而使用基于靜力隱式算法AutoForm和QuickForm 軟件分別需要782 s 和663 s，可見隱式多步成形有限元法的效率大幅提升，大約僅為顯式有限元法的19%，而QuickForm 較AutoForm 的運(yùn)行時(shí)間減少了 15.38%，這說明了集成到QuickForm 軟件中的LB 算子方程方法來預(yù)示中間構(gòu)形是快速且有效的。

圖10 某復(fù)雜零件最終構(gòu)形的外輪廓圖 Fig.10 Contours of final part for a complicated part

5 結(jié)論

在多步成形有限元法中，本文將中間構(gòu)形的構(gòu)造分解為拉伸變形和彎曲變形的解耦格式；根據(jù)大位移應(yīng)變的理論，推導(dǎo)了中間構(gòu)形的節(jié)點(diǎn)位移約束，給出了利用LB 方程求解中間構(gòu)形的初始解預(yù)示算法。通過數(shù)值實(shí)驗(yàn)結(jié)果表明：

(1) 在多步有限元法中采用解耦格式，增加了每步的行進(jìn)步長(zhǎng)，有效地改善了方程組的條件數(shù)，具有更好地?cái)?shù)值求解穩(wěn)定性和收斂性，提高了有限元法的計(jì)算效率。

(2) 基于大位移應(yīng)變理論推導(dǎo)位移約束，能夠有效地模擬出彎曲變形過程中的節(jié)點(diǎn)位置變化移動(dòng)情況。

(3) 通過離散的LB 算子線性方程來進(jìn)行中間構(gòu)形的初始解預(yù)示，算法實(shí)現(xiàn)簡(jiǎn)單、通用性強(qiáng)、計(jì)算速度快。