基于解耦計(jì)算的多步快速成形有限元法*

鮑益東，席 潔，秦雪嬌，劉永財(cái)

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

鈑金件有著強(qiáng)度高、重量輕、性能好、可大規(guī)模量產(chǎn)的優(yōu)點(diǎn)，并且成形鈑金件比常見(jiàn)的機(jī)械減材加工工藝具有更高的加工效率和材料使用率，因此被廣泛應(yīng)用于航空、航天、汽車(chē)、船舶、電子器件等工業(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è)計(jì)制造周期，節(jié)約研發(fā)成本，計(jì)算機(jī)技術(shù)和數(shù)值模擬技術(shù)逐漸應(yīng)用于生產(chǎn)當(dāng)中[1–2]。

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

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

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

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

一步快速成形有限元法基于全量理論，將零件成形過(guò)程簡(jiǎn)化為近似比例加載的過(guò)程，只考慮初始板料構(gòu)形及最終零件構(gòu)形兩種狀態(tài)，其基本思想為：從板料的初始構(gòu)形C0出發(fā)，C0上各節(jié)點(diǎn)P0的坐標(biāo)已知，在滿(mǎn)足一定邊界條件的情況下，通過(guò)有限元方法可以確定與初始板料上各節(jié)點(diǎn)P0對(duì)應(yīng)的最終構(gòu)形C上的節(jié)點(diǎn)P的坐標(biāo)，對(duì)比兩者坐標(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è)板料的變形滿(mǎn)足Kirchhoff假設(shè)，則初始板料和最終零件構(gòu)形上的任意質(zhì)點(diǎn)的幾何關(guān)系如圖2所示 （其中，Up為質(zhì)點(diǎn)p的位移函數(shù)，即質(zhì)點(diǎn)p和p0之間的位移向量；n為質(zhì)點(diǎn)p的單位法向量；n0為質(zhì)點(diǎn)p0的單位法向量；z為q點(diǎn)的厚向坐標(biāo)；z0為q0點(diǎn)的厚向坐標(biāo)；h為最終構(gòu)形的板料厚度；h0為初始板料厚度）。將初始狀態(tài)板料作為參考構(gòu)形，比較當(dāng)前構(gòu)形與初始構(gòu)形上對(duì)應(yīng)的質(zhì)點(diǎn) （如質(zhì)點(diǎn)q和質(zhì)點(diǎn)q0）的位置坐標(biāo)，即可得到相應(yīng)質(zhì)點(diǎn)的位移函數(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點(diǎn)和q0點(diǎn)的坐標(biāo)向量。

采用Lagrange 描述質(zhì)點(diǎn)的運(yùn)動(dòng)，則有

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

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

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

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

最終可以得到單元沿厚度方向任意一點(diǎn)變形的對(duì)數(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)變曲線(xiàn)的割線(xiàn)模量，即等效應(yīng)力與等效應(yīng)變的比值；矩陣[P]由面內(nèi)各向異性系數(shù)R0、R45、R90求得，即

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 算例驗(yàn)證

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

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

圖7 材料應(yīng)力–應(yīng)變曲線(xiàn)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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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