鮑益東，席 潔，秦雪嬌，劉永財

（1.南京航空航天大學(xué)，南京 210016；2.常州工學(xué)院，常州 213032）

鈑金件有著強度高、重量輕、性能好、可大規(guī)模量產(chǎn)的優(yōu)點，并且成形鈑金件比常見的機(jī)械減材加工工藝具有更高的加工效率和材料使用率，因此被廣泛應(yīng)用于航空、航天、汽車、船舶、電子器件等工業(yè)產(chǎn)業(yè)。隨著現(xiàn)代飛機(jī)產(chǎn)量需求的加大和制造要求的提高，飛機(jī)零部件的生產(chǎn)愈發(fā)精密復(fù)雜，并且生產(chǎn)周期短、質(zhì)量要求高。鈑金零件是飛機(jī)機(jī)體構(gòu)成的重要零件，在飛機(jī)零件總量中可占30%～50%，沖壓成形是飛機(jī)鈑金件制造的關(guān)鍵技術(shù)之一，飛機(jī)鈑金零件的生產(chǎn)效率和成形精度直接影響著飛機(jī)的生產(chǎn)周期和整機(jī)質(zhì)量，因此鈑金零件的生產(chǎn)也需要保證高質(zhì)量、高精度、短周期和低成本。為了縮短鈑金件的設(shè)計制造周期，節(jié)約研發(fā)成本，計算機(jī)技術(shù)和數(shù)值模擬技術(shù)逐漸應(yīng)用于生產(chǎn)當(dāng)中[1–2]。

有限元法是板材成形數(shù)值模擬的主要分析方法，主要分為基于流動理論的增量法和基于形變理論的全量法。增量法的主要優(yōu)點在于模擬過程中考慮了加載歷史的影響，模擬結(jié)果精度相對較高，但計算時間較長，難以滿足設(shè)計階段快速評價產(chǎn)品可制造性的需求[3]。法國的Batoz[4]和Naceur[5]等首先提出了基于全量理論的一步逆成形有限元法，基本原理是簡化設(shè)置整個成形過程，將零件的最終形態(tài)和初始板料形態(tài)進(jìn)行對比，簡單且快速地求出零件最終形態(tài)的應(yīng)力應(yīng)變分布以及初始板料形狀。Li 等[6]提出在一步逆成形有限元法中用小應(yīng)變方程解決大變形問題。我國對一步逆成形有限元法的研究從20世紀(jì)90年代開始，上海交通大學(xué)的沈啟彧[7]和Wang[8]等對方盒件的成形使用一步逆成形有限元算法進(jìn)行了預(yù)示。周平等[9]對一步逆成形法的求解方法進(jìn)行了改進(jìn)，提出一種基于切平面的解耦列式法，使一步逆成形法的計算精度得到提高。在一步逆成形有限元法的基礎(chǔ)上，具有一定工程實用性的一步正成形有限元法發(fā)展起來，與一步逆成形有限元法相比，一步正成形有限元法的求解過程更加符合零件的實際成形過程，根據(jù)規(guī)則的初始坯料計算得到成形后零件的形狀。一步成形有限元法 （包括一步逆成形和一步正成形）計算效率雖然很高，但是由于計算過程中采用了過多的簡化處理與條件假設(shè)，忽略了變形歷史，使得應(yīng)力計算的精度無法達(dá)到預(yù)期。為了改善這一現(xiàn)狀，研究人員多采用引入中間構(gòu)形的方式來模擬成形過程中的變形歷史。Lee 等[3]在單步逆成形法基礎(chǔ)上提出了多步逆成形法，通過線彈性映射法計算成形過程中的中間構(gòu)形，并根據(jù)數(shù)值結(jié)果驗證了可行性，但對于復(fù)雜零件難以找到對應(yīng)的映射函數(shù)進(jìn)行計算，具有一定的局限性。Guo 等[10]提出偽逆成形法，采用最小面積法構(gòu)造成形加載過程中的中間構(gòu)形，但只考慮軸對稱問題，對于非對稱問題并不適用。唐炳濤等[11]采用雙截面法和單元面積最小的方法求解中間構(gòu)形初始解，可以求解非對稱問題，但是容易出現(xiàn)局部網(wǎng)格嚴(yán)重畸變的問題。陳偉等[12]對多步正成形算法中的中間構(gòu)形初始解預(yù)示算法進(jìn)行研究，將初始平板料投影到最終構(gòu)形上后，采用比例插值的方法計算出中間構(gòu)形，并對中間構(gòu)形進(jìn)行接觸修正。

目前，國內(nèi)外對中間構(gòu)形基本理論的研究尚未成熟，在構(gòu)造復(fù)雜形面的中間構(gòu)形時仍然存在一些問題，對于構(gòu)形間應(yīng)力更新算法方面的研究相對更少，因此要做到高效準(zhǔn)確地預(yù)示中間構(gòu)形，從而提高零件成形模擬的精度仍然需要更加深入透徹的研究。

本文在一步快速成形有限元法（一步正成形有限元法）的理論基礎(chǔ)上，研究多步快速成形有限元法，提出了一種基于解耦格式的中間構(gòu)形預(yù)示算法，將中間構(gòu)形的求解過程解耦為彎曲變形和拉伸變形兩個階段，解決了迭代求解復(fù)雜非線性方程組容易出現(xiàn)病態(tài)矩陣的問題，提高了求解效率，并且可以適應(yīng)復(fù)雜零件形狀的求解。最后，通過實際數(shù)值算例驗證了該算法的有效性和準(zhǔn)確性，可以用于指導(dǎo)鈑金件的實際生產(chǎn)。

1 一步快速成形有限元法基本理論

一步快速成形有限元法基于全量理論，將零件成形過程簡化為近似比例加載的過程，只考慮初始板料構(gòu)形及最終零件構(gòu)形兩種狀態(tài)，其基本思想為：從板料的初始構(gòu)形C0出發(fā)，C0上各節(jié)點P0的坐標(biāo)已知，在滿足一定邊界條件的情況下，通過有限元方法可以確定與初始板料上各節(jié)點P0對應(yīng)的最終構(gòu)形C上的節(jié)點P的坐標(biāo)，對比兩者坐標(biāo)的變化可以獲得成形后零件的應(yīng)力應(yīng)變狀態(tài)以及厚度分布情況，如圖1所示。

圖1 一步快速成形有限元法示意圖Fig.1 Diagram of one-step quick forming finite element method

1.1 幾何關(guān)系

一步快速成形有限元法假設(shè)板料的變形滿足Kirchhoff假設(shè)，則初始板料和最終零件構(gòu)形上的任意質(zhì)點的幾何關(guān)系如圖2所示 （其中，Up為質(zhì)點p的位移函數(shù)，即質(zhì)點p和p0之間的位移向量；n為質(zhì)點p的單位法向量；n0為質(zhì)點p0的單位法向量；z為q點的厚向坐標(biāo)；z0為q0點的厚向坐標(biāo)；h為最終構(gòu)形的板料厚度；h0為初始板料厚度）。將初始狀態(tài)板料作為參考構(gòu)形，比較當(dāng)前構(gòu)形與初始構(gòu)形上對應(yīng)的質(zhì)點 （如質(zhì)點q和質(zhì)點q0）的位置坐標(biāo)，即可得到相應(yīng)質(zhì)點的位移函數(shù)Uq，即

圖2 一步快速成形有限元法的幾何關(guān)系Fig.2 Geometric relations of one-step quick forming finite element method

式中，Uq=[U1U2U3]T；Xq=[XqYqZq]T；Xq0=[Xq0Yq0Zq0]T。Xq和Xq0分別為q點和q0點的坐標(biāo)向量。

采用Lagrange 描述質(zhì)點的運動，則有

式中，[F]為變形梯度張量，由此可得左柯西–格林 （Cauchy–Green）矩陣[B]的逆矩陣為

對于當(dāng)前構(gòu)形材料內(nèi)部任一點來說，可以用3 個主方向上的主伸長量來描述這一點的變形。通過計算[B]–1的特征值可以得到板料上各點的面內(nèi)主拉伸λ1、λ2，然后通過體積不變的原理可以得到單元厚度方向的主拉伸λ3為

單元應(yīng)變的主方向是由[B]–1的特征向量決定的，即

式中，[M]為單元應(yīng)變方向向量矩陣；n1、n2、n3分別為單元應(yīng)變3 個主方向的單位向量。

最終可以得到單元沿厚度方向任意一點變形的對數(shù)應(yīng)變?yōu)?/p>

1.2 本構(gòu)關(guān)系

一步快速成形有限元法的本構(gòu)方程是采用比例加載條件下的全量形變理論描述的，并將材料的各向異性特性考慮在內(nèi)，本構(gòu)方程為

式中，{σ}為應(yīng)力矩陣，{ε}為應(yīng)變矩陣；Es為材料應(yīng)力–應(yīng)變曲線的割線模量，即等效應(yīng)力與等效應(yīng)變的比值；矩陣[P]由面內(nèi)各向異性系數(shù)R0、R45、R90求得，即

1.3 有限元方程的建立與求解

根據(jù)虛功原理在最終構(gòu)形上建立平衡方程，即

式中，W為單元虛功；e為單元編號；{U*}為全局坐標(biāo)系下的虛位移；為節(jié)點內(nèi)力向量；為節(jié)點外力向量。

將所有的節(jié)點內(nèi)力向量和節(jié)點外力向量組裝在一起即可獲得整體殘余力向量，即

一步快速成形有限元法具有強烈的非線性特征，一般采用牛頓–拉普森 （Newton–Raphson）方法進(jìn)行迭代求解，即

式中，{Ui}為第i個迭代步的近似初始解；ωi為引入的迭代松弛系數(shù)，以保證求解能夠順利收斂。

根據(jù)一步快速成形有限元法的實際計算情況，選擇位移準(zhǔn)則和殘余力準(zhǔn)則進(jìn)行收斂性的綜合判斷，若迭代步小于設(shè)定閾值且迭代計算過程滿足任何一個收斂準(zhǔn)則，則認(rèn)為迭代收斂結(jié)束計算。

2 基于解耦計算的多步快速成形有限元法

一步快速成形有限元法只考慮初始板料構(gòu)形和最終零件構(gòu)形，忽略了變形的中間狀態(tài)，計算效率得以大幅提高，可以快速獲得成形后零件的應(yīng)力應(yīng)變狀態(tài)和厚度分布情況，但是對中間加載歷史的忽略也會導(dǎo)致最終零件分析精度不足。為解決這一問題，提出了多步快速成形有限元法，通過引入中間構(gòu)形的方式考慮加載路徑和變形歷史從而提高計算精度，每一個中間構(gòu)形都作為當(dāng)前時間步的最終構(gòu)形通過一步快速成形有限元法進(jìn)行迭代計算。因此，多步快速成形有限元法的關(guān)鍵在于準(zhǔn)確地構(gòu)造每個時間步的中間構(gòu)形，以及各中間構(gòu)形之間的應(yīng)力迭代更新。

對金屬薄板成形來說，變形主要由彎曲變形和拉伸變形構(gòu)成，由于剪切變形對整個板料成形結(jié)果影響比較小，可以忽略其影響效應(yīng)。在多步快速成形有限元法中，中間構(gòu)形的計算被分解為彎曲變形和拉伸變形的解耦形式，其中對彎曲變形起主要作用的是彎曲剛度，對拉伸變形起主要作用的是膜剛度。通過解耦技術(shù)，首先進(jìn)行彎曲變形過程的求解，獲得板料在當(dāng)前模具作用下的構(gòu)形，如圖3（a）所示；其次是進(jìn)行拉伸變形過程的求解，在上一步求解出的板料構(gòu)形的基礎(chǔ)上，限定材料的非線性流動在該構(gòu)形上，通過彈塑性流動力學(xué)進(jìn)行迭代求解來獲得當(dāng)前時間步的中間構(gòu)形，如圖3（b）所示。

圖3 解耦計算原理示意圖Fig.3 Diagram of principle of decoupling calculation

2.1 中間構(gòu)形滑移約束面的構(gòu)造

板料彎曲變形的求解過程就是板料中間構(gòu)形滑移約束面構(gòu)造過程，假設(shè)已計算出上一個時間步的板料構(gòu)形，則在該構(gòu)形的基礎(chǔ)上計算當(dāng)前時間步的中間構(gòu)形。其首要問題是準(zhǔn)確地得到板料和模具的接觸狀態(tài)，其中包括確定凸模、凹模、壓邊圈等工具與板料的接觸狀況，只有正確地判斷出與模具相接觸的板料節(jié)點，才能為中間構(gòu)形的預(yù)示算法提供計算依據(jù)。由于板料和模具的接觸狀態(tài)僅被用來計算板料的形狀，所以只需判斷哪些板料節(jié)點與工具發(fā)生接觸，而不用進(jìn)一步判斷板料節(jié)點在工具表面是發(fā)生滑動還是固定不動。在判斷板料節(jié)點與模具接觸狀態(tài)時，首先要通過搜索算法確定板料節(jié)點沿法線方向在模具表面上的投影單元，然后進(jìn)一步確定具體的投影點。在接觸搜索階段，為了縮小搜索范圍，加快接觸搜索的速度，采用劃分空間格的方法搜索距離板料節(jié)點最近的模具節(jié)點，根據(jù)節(jié)點與單元之間的拓?fù)潢P(guān)系確定模具表面的投影單元范圍。

確定了板料節(jié)點在模具表面的投影點后，根據(jù)幾何關(guān)系判斷節(jié)點是否穿透模具表面，如圖4所示。節(jié)點穿透量計算方程式為

圖4 節(jié)點穿透判斷示意圖Fig.4 Diagram of node penetration judgment

式中，g（xs）為節(jié)點xs的穿透量；xs為板料上某一節(jié)點；xm為節(jié)點xs在模具表面的投影點；n（xm）為節(jié)點xm所在投影單元的外法線方向。

當(dāng)g（xs）＞0 時，表明板料節(jié)點xs并未穿透模具表面；當(dāng)g（xs）＜0 時，表明板料節(jié)點xs穿透模具表面。判斷出所有發(fā)生穿透的板料節(jié)點后，將其沿節(jié)點法線方向直接投影到工具表面，假定這些發(fā)生穿透的板料節(jié)點就是可能與工具發(fā)生接觸的節(jié)點，根據(jù)發(fā)生接觸的板料節(jié)點會與工具表面產(chǎn)生作用力的特點排除其中一些誤判的板料節(jié)點，剩下的發(fā)生穿透的板料節(jié)點則是真正與工具發(fā)生接觸的板料節(jié)點。最后根據(jù)接觸算法，搜索出與工具發(fā)生接觸的板料節(jié)點信息，并確定出接觸區(qū)域的板料構(gòu)形。

在彎曲變形過程中，板料與凸凹模具接觸的區(qū)域受到模具面的約束作用，未接觸區(qū)域面內(nèi)會存在膜面的初始應(yīng)力而不承受其他力的作用，因此非接觸區(qū)域的構(gòu)形采用基于預(yù)應(yīng)力膜單元建立平衡方程的方法進(jìn)行計算。預(yù)應(yīng)力膜單元頂點的位移向量為a=[uvw]T（其中u、v和w分別為x、y和z方向的位移分量）。預(yù)應(yīng)力膜單元在變形過程中的幾何關(guān)系為

式中，ε為應(yīng)變；εL、εNL分別為線性應(yīng)變和非線性應(yīng)變；B為應(yīng)變與位移的轉(zhuǎn)換矩陣，BL、BNL分別為線性應(yīng)變和非線性應(yīng)變分析的矩陣項，前者與位移向量a無關(guān)，后者則是a的函數(shù)。

預(yù)應(yīng)力膜單元具有大位移、小應(yīng)變的特點，因此假設(shè)材料的應(yīng)力–應(yīng)變關(guān)系滿足如下線彈性關(guān)系，即

式中，D為彈性矩陣；σ0為膜單元的預(yù)應(yīng)力。

在單元所在的局部坐標(biāo)系下，根據(jù)幾何關(guān)系和本構(gòu)關(guān)系，基于虛功原理建立平衡方程為

式中，V表示體積，對體積進(jìn)行積分；B和均為位移向量的函數(shù)，=BL+2BNL。

采用Newton–Raphson 方法進(jìn)行迭代求解，并根據(jù)位移準(zhǔn)則和殘余力準(zhǔn)則進(jìn)行收斂性的綜合判斷，若迭代步小于設(shè)定閾值且迭代計算過程滿足任何一個收斂準(zhǔn)則，則認(rèn)為迭代收斂結(jié)束計算。否則，需要采用縮小步長的方法重新計算板料非接觸區(qū)域中間構(gòu)形。

重復(fù)上述過程，直至中間構(gòu)形的所有單元都滿足預(yù)應(yīng)力膜單元的平衡條件，從而完成當(dāng)前時間步的中間構(gòu)形滑移約束面的構(gòu)造。

2.2 基于滑移約束面的中間構(gòu)形計算

在拉伸變形求解過程中，為了確保板料節(jié)點在平衡迭代過程中的移動始終被限定在板料中間構(gòu)形上，將2.1 節(jié)計算出的中間構(gòu)形作為板料拉伸變形計算過程的滑移約束面，約束板料的運動和變形?；萍s束面即節(jié)點的切平面區(qū)域，滑移約束面如圖5所示?？梢钥闯觯?jié)點坐標(biāo)系定義在節(jié)點切平面上，軸為節(jié)點法線方向，單元坐標(biāo)系（s，t）則定義在單元所在平面上。單元上的剛度矩陣ke和殘余應(yīng)力都將轉(zhuǎn)換到節(jié)點坐標(biāo)系下，轉(zhuǎn)換公式如下：

圖5 滑移約束面示意圖Fig.5 Diagram of slip constraint surface

經(jīng)過2.1 節(jié)中的幾何迭代判斷，滑移約束面上各節(jié)點與模具表面之間不會發(fā)生穿透，但在模具表面上圓角較小的地方，滑移約束面的網(wǎng)格單元會出現(xiàn)穿透模具表面的情況，如圖6（a）所示，因此，中間構(gòu)形上的節(jié)點雖然被約束在滑移約束面上，但在基于滑移約束面進(jìn)行迭代計算時，節(jié)點在滑移約束面上的移動會導(dǎo)致中間構(gòu)形上部分節(jié)點與模具表面發(fā)生穿透，需要對迭代求解的中間構(gòu)形進(jìn)行修正。首先按照2.1 節(jié)中的方法將約束在滑移約束面上的中間構(gòu)形節(jié)點與模具進(jìn)行幾何接觸判斷，搜索到入侵至模具內(nèi)部的節(jié)點，判斷出中間構(gòu)形上所有發(fā)生穿透的節(jié)點后，將這些節(jié)點拉回到模具表面，使其成為附著在模具表面的固定點。通過直接投影的方法將這些節(jié)點沿法線方向直接投影到工具表面，從而可以獲得修正后的中間構(gòu)形，如圖6（b）所示。

圖6 中間構(gòu)形修正過程示意圖Fig.6 Diagram of intermediate configuration modification process

3 算例驗證

本文所提出的基于解耦計算的多步快速成形有限元法算法已集成于自主開發(fā)的QuickForm 板料成形模擬軟件中。為驗證該算法的有效性，將計算結(jié)果與現(xiàn)有的基于隱式算法求解的板料成形商業(yè)有限元軟件AutoForm 軟件以及動力顯式求解的Dynaform 軟件的模擬結(jié)果進(jìn)行對比分析。翼子板具有造型復(fù)雜多變、曲率變化大，表面質(zhì)量和尺寸精度要求高等特點，因此成形難度大，容易出現(xiàn)開裂、起皺以及棱線滑移等缺陷，并且回彈量大難以準(zhǔn)確預(yù)測[13]。以翼子板的沖壓成形為例進(jìn)行分析，可以充分地驗證本文所提出算法的可行性。

零件所用材料為熱鍍鋅合金鋼板，材料牌號為GD53D–ZF–V22，材料參數(shù)為：抗拉強度370.4 MPa、屈服強度165.5 MPa、各向異性系數(shù)1.61、硬化指數(shù)0.22，應(yīng)力–應(yīng)變關(guān)系式為：σ-=502.3（0.00644+ε-）0.22MPa，曲線如圖7所示。初始板料厚度為0.65 mm，初始板料劃分三角形網(wǎng)格后如圖8所示，初始單元數(shù)量為21105 個，初始節(jié)點數(shù)量為10756個，板料成形過程中可根據(jù)板料接觸區(qū)域的工具表面曲率變化情況自適應(yīng)加密板料網(wǎng)格尺寸，加密級數(shù)為3 級，板料與工具之間的摩擦系數(shù)為0.15。板料與模具的裝配如圖9所示，從下到上依次是凸模、壓邊圈、板料以及凹模，凹模向下運動，當(dāng)凹模與壓邊圈接觸后開始拉延過程，總行程為424.23 mm，壓邊圈行程為85 mm。采用的計算機(jī)配置為英特爾Core i5–10400@2.90 GHz 六核/12 線程的CPU，32 GB 內(nèi)存。

圖7 材料應(yīng)力–應(yīng)變曲線Fig.7 Material stress–strain curve

圖8 初始板料網(wǎng)格示意圖Fig.8 Diagram of initial sheet grid

圖9 板料與模具裝配圖Fig.9 Assembly drawing of sheet and die

拉延過程中，工具向下運動接觸板料后，板料受到力的作用開始發(fā)生變形，該過程采用自適應(yīng)時間步長控制，時間步長即工具向下移動的距離。有效的自動時間步長控制可以確保接觸約束只在從一個增量到下一個增量之間適度改變。因此要選擇合適的時間步長，既可以提高計算效率也可以減少對計算精度的影響，不會出現(xiàn)難以收斂的問題。本算法設(shè)置的拉延過程最大時間步長為5 mm，初始時刻步長為2 mm，發(fā)生的變形較小，時間步長自動增加為5 mm，當(dāng)即將合模時，為避免迭代發(fā)散，時間步長減少至0.5 mm。每一個時間步長都對應(yīng)一個經(jīng)過彎曲和拉伸變形求解之后的中間構(gòu)形，中間構(gòu)形的變化反映了成形過程中的材料的流動情況，板料初始狀態(tài)至合模狀態(tài)之間若干時間步長的中間構(gòu)形如圖10所示。

圖10 板料變形過程示意圖Fig.10 Diagram of sheet deformation process

上下模距離為12 mm 時的中間構(gòu)形網(wǎng)格如圖11（a）所示，構(gòu)形上特征明顯的區(qū)域網(wǎng)格經(jīng)過多次加密，節(jié)點數(shù)目較多，特征不明顯的區(qū)域網(wǎng)格沒有被加密，節(jié)點數(shù)目較少。為驗證中間構(gòu)形預(yù)示算法的有效性，在圖11（a）所示中間構(gòu)形上做縱向及橫向的截面線，與AutoForm 對應(yīng)位置的構(gòu)形截面線進(jìn)行對比。圖11（b）為沿截面線AA'上分布的曲線形狀，可以發(fā)現(xiàn)QuickForm 和AutoForm軟件計算的曲線形狀及拉延深度整體比較接近，但是局部區(qū)域存在一些差別。圖11（c）為沿截面線BB'上分布的曲線形狀，兩個軟件計算的曲線形狀及拉延深度十分接近，說明了本文提出的中間構(gòu)形預(yù)示算法是有效且可行的，能夠較好地反映材料的流動情況。

圖11 上下模距離為12 mm 時的中間構(gòu)形形狀結(jié)果Fig.11 Intermediate configuration shape results when the upper and lower die distance is 12 mm

合模狀態(tài)下成形性和厚度分布對比見圖12和13，其中一步快速成形有限元法與多步快速成形有限元法的最終模擬結(jié)果如圖12（a）和（b）、圖13（a）和（b）所示。一步快速成形有限元法沒有中間構(gòu)形的預(yù)示過程，因此模擬所得的零件成形性稍差，零件上特征復(fù)雜的區(qū)域成形不充分，最大厚度與最小厚度也偏大，精度稍差。一步快速成形有限元法的計算時間為29 s，多步快速成形有限元法 （QuickForm）的計算時間為216 s，QuickForm 在保證了模擬精度的情況下仍可在幾分鐘內(nèi)完成計算，保證了求解效率。

從圖12可以看出，QuickForm、AutoForm 二者成形性的整體分布情況接近，安全、變形不足，起皺趨勢和趨勢區(qū)的位置基本一致。從圖13可以看出，QuickForm、AutoForm 模擬的板料厚度整體分布趨勢一致，且最大增厚和最大減薄的位置一致。最大厚度QuickForm 的計算結(jié)果為0.691 mm，AutoForm 為0.673 mm，相對誤差為2.7%；最小厚度QuickForm的計算結(jié)果為0.526 mm，AutoForm為0.499 mm，相對誤差為5.4%。相較于QuickForm 和AutoForm 兩種軟件，Dynaform 模擬的零件整體成形性較差，變形不足的區(qū)域較大，且最大增厚和最大減薄的位置與另外兩種軟件相比有所差別，出現(xiàn)差距的主要原因是Dynaform 軟件求解采用的是動力顯式算法，該算法會引起慣性效應(yīng)，從而使拉延結(jié)果受到影響。以上模擬結(jié)果說明了本文所提出的算法計算結(jié)果是合理可行的。在求解時間上，QuickForm 的計算總時間為216 s，AutoForm計算總時間為165 s，Dynaform計算總時間為39 min+48 s，QuickForm 計算效率要遠(yuǎn)高于Dynaform 但是略低于AutoForm 軟件。QuickForm 的算法目前是單核計算，尚未實現(xiàn)并行求解，因此計算效率低于AutoForm，但是考慮到未來實現(xiàn)本文提出算法的并行計算版本后，計算效率有很大的提升空間。

圖12 合模狀態(tài)下成形性分布對比Fig.12 Comparison of formability distribution under clamping state

圖13 合模狀態(tài)下厚度分布對比Fig.13 Comparison of thickness distribution under clamping state

基于上述仿真模擬的工藝參數(shù)，進(jìn)行零件的拉深成形試驗，獲得成形后的產(chǎn)品如圖14所示，翼子板外觀無開裂、起皺等缺陷，成形質(zhì)量較好。將試驗制備的最終成形零件的截面厚度分布、外輪廓線等結(jié)果與QuickForm 自研算法的計算結(jié)果進(jìn)行對比，如圖15所示，可以發(fā)現(xiàn)試驗結(jié)果與QuickForm 的模擬結(jié)果基本一致，模擬精度可以滿足實際生產(chǎn)需求。

圖14 成形后的產(chǎn)品Fig.14 Finished product

圖15 最終成形零件與QuickForm 的計算結(jié)果對比Fig.15 Comparison of final formed parts and calculation results of QuickForm

4 結(jié)論

（1）在一步快速成形有限元法的基礎(chǔ)上，引入了中間構(gòu)形的多步快速成形有限元法，并提出了一種基于解耦計算的中間構(gòu)形構(gòu)造方法，將板料成形過程中的變形解耦為彎曲變形和拉伸變形兩個過程，提高了求解效率，并通過算例的模擬結(jié)果對比驗證了本文提出的中間構(gòu)形構(gòu)造方法的有效性。

（2）以典型鈑金件翼子板為案例進(jìn)行模擬分析，并將計算結(jié)果與AutoForm、Dynaform 軟件的計算結(jié)果進(jìn)行對比，驗證了本文所提出算法對鈑金零件成形預(yù)測的準(zhǔn)確性以及高效率；在相同工藝參數(shù)下進(jìn)行成形試驗，試驗結(jié)果進(jìn)一步表明本文所提出算法的可行性，模擬精度符合實際生產(chǎn)要求，可以在鈑金件生產(chǎn)設(shè)計初期進(jìn)行參數(shù)優(yōu)化，從而降低模具制造成本，縮短產(chǎn)品試制周期。