基于h型自適應(yīng)有限元法在薄板沖壓成型中的應(yīng)用

劉習(xí)洲，王城璟，王 琥

基于型自適應(yīng)有限元法在薄板沖壓成型中的應(yīng)用

劉習(xí)洲，王城璟，王 琥

(湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國家重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410082)

在對薄板沖壓成型這一過程進(jìn)行有限元仿真分析時(shí)，難以精確分析場變量發(fā)生劇烈變化的應(yīng)力集中及應(yīng)變梯度大的區(qū)域，如何平衡精度和效率間的關(guān)系是沖壓成型仿真的關(guān)鍵。因此，基于非線性有限元大變形的相關(guān)理論，針對動(dòng)態(tài)仿真的網(wǎng)格自適應(yīng)關(guān)鍵技術(shù)，建立了自適應(yīng)分析模式下的薄板沖壓成型算法。為了提高計(jì)算精度，提出基于單元應(yīng)變能增量的能量誤差準(zhǔn)則以及基于板料成型幾何特征的幾何誤差準(zhǔn)則，并結(jié)合2類誤差準(zhǔn)則，建立了基于自適應(yīng)分析算法的誤差判斷準(zhǔn)則；為了提高計(jì)算效率，引入阻尼因子提出了阻尼子循環(huán)算法，將自適應(yīng)加密后的板料單元按時(shí)間步長劃分為若干個(gè)區(qū)域，并按照不同時(shí)間步長分別求解。結(jié)果表明，該算法有效提高了薄板沖壓成型有限元仿真的精度與效率。

非線性有限元分析；薄板沖壓成型；自適應(yīng)網(wǎng)格加密；誤差估計(jì)；子循環(huán)算法

薄板沖壓成型通過模具在金屬板料上施加壓力，板料在外力的作用下發(fā)生塑性變形，從而形成所需零件的形狀。在這一復(fù)雜的非線性大變形過程中，包含材料、幾何以及接觸三重非線性問題，在進(jìn)行有限元仿真分析時(shí)，難以精確分析場變量發(fā)生劇烈變化的應(yīng)力集中區(qū)域。為了提高仿真精度，需要對全域網(wǎng)格進(jìn)行加密，縮小網(wǎng)格尺寸，以達(dá)到提升精度的目的，因此也極大地增加了計(jì)算量，使得計(jì)算效率嚴(yán)重低下。在這種情況下，對局部網(wǎng)格進(jìn)行加密的自適應(yīng)技術(shù)[1]應(yīng)運(yùn)而生。自適應(yīng)技術(shù)通過利用中間變量作為誤差估計(jì)[2]來定位應(yīng)力集中的區(qū)域，并采取相應(yīng)的網(wǎng)格劃分方法進(jìn)行網(wǎng)格加密，通過局部提升場變量變化劇烈區(qū)域網(wǎng)格密度的方法達(dá)到分析精度的要求，同時(shí)減少了計(jì)算消耗，提高了計(jì)算效率。

影響自適應(yīng)分析結(jié)果的因素有誤差估計(jì)和自適應(yīng)策略。一般來說，誤差估計(jì)包括先驗(yàn)估計(jì)和后驗(yàn)估計(jì)2種。先驗(yàn)估計(jì)是基于解的光滑性等特性，并提供了當(dāng)自由度數(shù)趨于無窮時(shí)解的漸近收斂速度的定性信息。盡管先驗(yàn)估計(jì)對于某一類解的最大誤差是準(zhǔn)確的，如所有二階導(dǎo)數(shù)為常數(shù)的解，但其通常無法提供有限元模型的實(shí)際誤差。與之對應(yīng)的后驗(yàn)估計(jì)是利用求解過程中獲得的信息，以及一些關(guān)于解的先驗(yàn)假設(shè)。與不完整信息相關(guān)問題的抽象處理，可以提供離散化誤差的定量精確度量。在后驗(yàn)誤差估計(jì)中，基于殘差的誤差估計(jì)和基于恢復(fù)的誤差估計(jì)是最為常用的方法?；跉埐畹恼`差估計(jì)量是通過考慮平衡殘差恢復(fù)平衡通量，在每個(gè)單元中得到誤差估計(jì)，其包括顯式和隱式2種殘差估計(jì)。基于恢復(fù)的誤差估計(jì)量是通過使用從數(shù)值結(jié)果的后驗(yàn)處理得到恢復(fù)解，并且任何能減少誤差的恢復(fù)過程都會(huì)給出合理的誤差估計(jì)，如果恢復(fù)解的收斂速度比有限元解的收斂速度快，那么將得到更為接近精確的估計(jì)。這是由ZIENKIEWICZ和ZHU[3]首先提出，現(xiàn)已被廣泛應(yīng)用于有限元方法中；此外，JIANG等[4]提出一種新的基于單元能量投影技術(shù)的二維非線性有限元分析方法；CHEN等[5]使用基于恢復(fù)的后驗(yàn)誤差估計(jì)作為Allen-Cahn方程的誤差準(zhǔn)則。

在誤差準(zhǔn)則建立之后，需要采用適當(dāng)?shù)淖赃m應(yīng)策略應(yīng)用于自適應(yīng)分析過程中，其目的是使有限元計(jì)算的估計(jì)解逐步逼近于真實(shí)解。目前較為成熟且應(yīng)用廣泛的自適應(yīng)策略有4種類型，分別是型、型、型和型，如圖1所示。

圖1 自適應(yīng)策略類型

型自適應(yīng)策略基于誤差準(zhǔn)則通過對高誤差區(qū)域增加節(jié)點(diǎn)數(shù)量來減小單元的尺寸[6-7]；型通過提高單元形函數(shù)的階次[8]；型通過節(jié)點(diǎn)的重新定位[9]，但不改變單元的類型和數(shù)目，網(wǎng)格的連接性也保持不變[10]；型策略同時(shí)采用型和型自適應(yīng)策略[11]，結(jié)合了二者的優(yōu)勢，但具有難度大、不穩(wěn)定的特點(diǎn)，較難應(yīng)用到薄板沖壓成型中。

本文選擇型自適應(yīng)策略應(yīng)用到薄板沖壓成型過程中，并提出一種單元應(yīng)變能增量和板料變形幾何特征相結(jié)合的自適應(yīng)誤差準(zhǔn)則，可自動(dòng)識別局部誤差超過預(yù)定極限的應(yīng)變梯度及應(yīng)力集中區(qū)域并進(jìn)行標(biāo)記，利用自適應(yīng)網(wǎng)格加密策略對該區(qū)域網(wǎng)格單元進(jìn)行自適應(yīng)加密。為了提高計(jì)算效率，引入阻尼因子，提出了阻尼子循環(huán)算法，將自適應(yīng)加密后的板料單元按照時(shí)間步長劃分為若干個(gè)區(qū)域，每個(gè)區(qū)域按照時(shí)間步長的不同分別進(jìn)行數(shù)值積分。算例結(jié)果表明，該算法有效地提高了薄板沖壓成型過程的仿真精度和效率。

1 自適應(yīng)分析誤差準(zhǔn)則

在自適應(yīng)分析過程中將有限元模型某些變量的中間結(jié)果視作誤差估計(jì)來進(jìn)行誤差分析，通過對誤差的分析來探求誤差分布的真實(shí)情況，使其為網(wǎng)格和節(jié)點(diǎn)的密度設(shè)置提供指導(dǎo)，從而有效地提升有限元仿真的精度。因此，如何設(shè)計(jì)誤差準(zhǔn)則將直接影響到自適應(yīng)分析結(jié)果的好壞，在自適應(yīng)分析過程中設(shè)置一個(gè)可靠、穩(wěn)定的誤差準(zhǔn)則是非常關(guān)鍵。

誤差指計(jì)算過程中精確解和近似解之間的差值，即

然而用式(1)～(3)說明局部誤差是不穩(wěn)妥的。例如在點(diǎn)載荷作用下，位移和應(yīng)力的局部誤差為無限大。因此，需要引入表示某個(gè)積分標(biāo)量的能量范數(shù)來表示誤差。對于彈塑性問題，采用恢復(fù)解后位移的能量范數(shù)誤差估計(jì)形式為

誤差的精度采用效率指數(shù)為

有限元自適應(yīng)分析過程中需要滿足2個(gè)條件：

(1) 在整個(gè)離散域中，有

(2) 誤差需均勻分布到每一個(gè)單元上，即要求每個(gè)單元滿足

針對恢復(fù)解給出一個(gè)更高階的誤差，則

針對薄板沖壓成型這一非線性大變形過程，可以采用單元應(yīng)變能的增量作為局部網(wǎng)格加密的準(zhǔn)則。本算法中采用了Belytschko-Tsay四邊形殼單元，定義如圖2所示的局部隨動(dòng)坐標(biāo)系，表達(dá)式為

在隨動(dòng)坐標(biāo)系中計(jì)算每個(gè)單元的應(yīng)變能增量為

其中，和分別為節(jié)點(diǎn)和的編號，,=1,2,3,4且≠；與分別為節(jié)點(diǎn)和投影表面凸模單元法向量；為兩兩法向量之間形成的夾角；為臨界夾角值。

圖3 自適應(yīng)網(wǎng)格加密幾何誤差準(zhǔn)則

將單元應(yīng)變能增量的能量誤差準(zhǔn)則與板料成型幾何特征的幾何誤差準(zhǔn)則相結(jié)合作為自適應(yīng)誤差準(zhǔn)則，最終自適應(yīng)加密的單元選擇2個(gè)準(zhǔn)則標(biāo)記出的所有單元，從而可以提高自適應(yīng)分析的精確性。

2 自適應(yīng)分析算法

2.1 非協(xié)調(diào)性網(wǎng)格的處理

在采用型自適應(yīng)有限元法中對板料應(yīng)力集中及應(yīng)變梯度大的區(qū)域網(wǎng)格加密后，板料網(wǎng)格單元會(huì)產(chǎn)生疏密不同的情況，如何處理粗網(wǎng)格和細(xì)網(wǎng)格區(qū)域之間的協(xié)調(diào)性成為自適應(yīng)分析過程的關(guān)鍵。如圖4所示，在新生成的網(wǎng)格單元，和過渡網(wǎng)格單元之間的邊界上，采用插值的方法生成節(jié)點(diǎn)，該點(diǎn)是新生成網(wǎng)格單元和的節(jié)點(diǎn)，而與過渡網(wǎng)格單元無關(guān)，則稱其為非協(xié)調(diào)網(wǎng)格節(jié)點(diǎn)。在仿真計(jì)算中，若不對非協(xié)調(diào)節(jié)點(diǎn)進(jìn)行處理，會(huì)導(dǎo)致新生成單元，和過渡單元之間邊界不連續(xù)，無法得出正確結(jié)果。因此，本文采用DEVLOO等[12]提出的1-不規(guī)則原則對非協(xié)調(diào)節(jié)點(diǎn)進(jìn)行處理，即用多點(diǎn)約束法引入約束條件將非協(xié)調(diào)網(wǎng)格節(jié)點(diǎn)約束在邊12上

其中，為位移場，新生成的非協(xié)調(diào)網(wǎng)格節(jié)點(diǎn)稱作節(jié)點(diǎn)1和2控制的約束從節(jié)點(diǎn)，節(jié)點(diǎn)1和2則稱作主節(jié)點(diǎn)，起到控制非協(xié)調(diào)節(jié)點(diǎn)的作用。

圖4 多點(diǎn)約束法處理非協(xié)調(diào)網(wǎng)格

Fig. 4 The multi-point constraint method deals with the inconsistent element

虛功原理可將加密后板料單元約束從節(jié)點(diǎn)的節(jié)點(diǎn)力等價(jià)到其相鄰的主節(jié)點(diǎn)上，即

將式(23)代入式(22)，得

式(24)可將作用在從節(jié)點(diǎn)上的節(jié)點(diǎn)力等價(jià)到與之相鄰的主節(jié)點(diǎn)上，采取這種做法的原因是在進(jìn)行顯式計(jì)算中，控制方程不包含從節(jié)點(diǎn)，通過相對應(yīng)的主節(jié)點(diǎn)才能將從節(jié)點(diǎn)的作用表現(xiàn)出來。此外，在網(wǎng)格自適應(yīng)過程中，除了等效約束從節(jié)點(diǎn)的內(nèi)力和外力外，還需要對單元的質(zhì)量矩陣進(jìn)行修正，否則將出現(xiàn)整體板料質(zhì)量增加或減少過多的情況，這就要求在自適應(yīng)網(wǎng)格加密的過程中，約束從節(jié)點(diǎn)上的質(zhì)量矩陣需等效到主節(jié)點(diǎn)上，以解決網(wǎng)格非協(xié)調(diào)性的問題。

2.2 自適應(yīng)網(wǎng)格加密策略

如圖5所示，假設(shè)板料單元2經(jīng)過自適應(yīng)誤差準(zhǔn)則判斷后誤差不滿足分析精度要求，在單元4條邊的中點(diǎn)分別插入新的節(jié)點(diǎn)對其進(jìn)行加密，原一個(gè)單元被劃分為4個(gè)新單元。其中節(jié)點(diǎn)為約束從節(jié)點(diǎn)，原有的單元編號刪除，生成新單元的編號按照順時(shí)針排列。假設(shè)新生成的單元3，4，5和6誤差仍不滿足精度要求，程序?qū)⒗^續(xù)對這4個(gè)單元按照之前的步驟加密。在繼續(xù)加密后，邊界上出現(xiàn)了，和3個(gè)非協(xié)調(diào)節(jié)點(diǎn)，且不滿足1-不規(guī)則原則。因此，首先應(yīng)對單元1進(jìn)行加密，將節(jié)點(diǎn)轉(zhuǎn)化為主節(jié)點(diǎn)，再對單元3，4，5和6加密，以保證相鄰單元之間最多只存在一個(gè)約束從節(jié)點(diǎn)，流程圖如圖6所示。經(jīng)過上述處理，按照自適應(yīng)誤差準(zhǔn)則對誤差分布較大區(qū)域進(jìn)行多次加密，可提高求解精度。

圖5 自適應(yīng)網(wǎng)格加密策略

圖6 自適應(yīng)網(wǎng)格加密流程圖

2.3 阻尼子循環(huán)算法

經(jīng)過自適應(yīng)加密后的板料單元與其他未加密單元相比，其時(shí)間步長僅為未加密單元的1/2或1/4，在多次加密后會(huì)更小。如果對所有的板料單元均對最小的時(shí)間步長進(jìn)行積分，將會(huì)產(chǎn)生較大的計(jì)算量，增加求解難度，降低計(jì)算效率。為解決該問題，本文算法在常速度子循環(huán)算法的基礎(chǔ)上提出阻尼子循環(huán)算法。

如圖7所示，區(qū)域表示的區(qū)域時(shí)間步長較小，區(qū)域表示的區(qū)域時(shí)間步長較大，且區(qū)域是區(qū)域時(shí)間步長的二分之一，換言之，當(dāng)位于區(qū)域中的節(jié)點(diǎn)變量更新2次，區(qū)域則更新1次。子循環(huán)算法將板料單元節(jié)點(diǎn)根據(jù)時(shí)間步長劃分為若干個(gè)區(qū)域，每個(gè)區(qū)域采用不同的時(shí)間步長進(jìn)行數(shù)值積分，從而能夠提高模擬仿真的計(jì)算效率。

圖7 子循環(huán)算法數(shù)據(jù)更新

在薄板沖壓成型過程中，最常用的是常速度子循環(huán)算法。在按照時(shí)間步長大小分區(qū)之后，需要考慮各區(qū)域步長的倍數(shù)關(guān)系，倍數(shù)可以是整數(shù)，也可以是非整數(shù)。如果是非整數(shù)倍關(guān)系，需首先確定最小時(shí)間步長，各個(gè)區(qū)域的時(shí)間步長必須是這個(gè)最小時(shí)間步長的整數(shù)倍，但相鄰2個(gè)區(qū)域的時(shí)間步長不要求為整數(shù)倍關(guān)系。在進(jìn)行顯式計(jì)算時(shí)需采用時(shí)鐘驅(qū)動(dòng)的模式，其最小間隔就是最小時(shí)間步長，凡是時(shí)間在時(shí)鐘后的所有節(jié)點(diǎn)數(shù)據(jù)必須更新。對于各個(gè)區(qū)域時(shí)間步長為整數(shù)倍時(shí)，需按照2個(gè)原則進(jìn)行劃分：

(1) 最大時(shí)間步長是各個(gè)時(shí)間步長的整數(shù)倍；

(2) 如果一個(gè)節(jié)點(diǎn)屬于2個(gè)不同的時(shí)間步長區(qū)域，則大時(shí)間步長必須是小時(shí)間步長區(qū)域步長的整數(shù)倍。

(1) 當(dāng)主循環(huán)開始時(shí)(初始0狀態(tài))，區(qū)域和區(qū)域的節(jié)點(diǎn)數(shù)據(jù)同時(shí)更新，有

常速度子循環(huán)算法是假設(shè)大時(shí)間步長區(qū)域的節(jié)點(diǎn)在一個(gè)主循環(huán)按常速度運(yùn)動(dòng)，很顯然，這種設(shè)置忽略了加速度的存在，必然會(huì)帶來一定的誤差，為了消除其引起的誤差導(dǎo)致精度上的影響和不穩(wěn)定因素，本文提出了阻尼子循環(huán)算法，即在算法中引入阻尼因子，子循環(huán)的變量更新時(shí)有

阻尼因子一般取0.005～0.010。

3 算例分析

本文以具有代表性的杯突成型算例和空調(diào)蓋板算例為例，檢驗(yàn)了網(wǎng)格自適應(yīng)分析算法的有效性和穩(wěn)定性。

3.1 杯突成型

杯突成型又稱埃里克森杯突試驗(yàn)，是板料成型試驗(yàn)中最常用的模型, 可以檢驗(yàn)板料單元在變形過程中的穩(wěn)定性。算例初始建立的有限元模型共有4 164個(gè)單元，4 307個(gè)節(jié)點(diǎn)，自由度為25 842。模具(凸模、壓邊圈和凹模)單元數(shù)為6 512，節(jié)點(diǎn)數(shù)為6 778。杯突成型的計(jì)算模型如圖8所示，仿真參數(shù)見表1。

圖8 杯突成型計(jì)算模型

表1 杯突成型仿真參數(shù)

在薄板沖壓成型中，首先進(jìn)行壓邊圈運(yùn)動(dòng)將金屬板料夾緊，然后凹模固定，凸模朝著凹模方向做直線運(yùn)動(dòng)。當(dāng)凸模與板料接觸時(shí)，離散后的有限元網(wǎng)格逐漸發(fā)生變形。如圖9(a)所示，當(dāng)板料與凸模接觸后，可滿足本文自適應(yīng)誤差準(zhǔn)則，根據(jù)該算法對網(wǎng)格進(jìn)行了第1次加密，對加密后的杯突模型繼續(xù)進(jìn)行有限元數(shù)值模擬。凸模在閉合階段繼續(xù)行進(jìn)過程中，網(wǎng)格變形產(chǎn)生的誤差再次達(dá)到自適應(yīng)誤差準(zhǔn)則，則對網(wǎng)格進(jìn)行第2次加密，如圖9(b)所示。隨著凸模行程的增加，越來越多的板料單元被加密，在圖9(c)中加密單元的數(shù)量逐漸增加。杯突成型的最后時(shí)刻如圖9(d)所示，與初始板料尺寸(紅色虛線框)相比，經(jīng)過沖壓成型后的板料尺寸有所減小，與實(shí)際工程情況相符。

圖9 杯突成型自適應(yīng)加密過程((a)第1次自適應(yīng)加密；(b)第2次自適應(yīng)加密；(c)第3次自適應(yīng)加密；(d)最終成型)

表2通過對比杯突成型過程中有無自適應(yīng)網(wǎng)格加密的厚向云圖分布可以看出，無自適應(yīng)加密板料的最大厚度為0.928 mm，最小為0.734 mm；自適應(yīng)加密后的板料最大厚度為0.926 mm，最小為0.757 mm。有無網(wǎng)格自適應(yīng)加密的厚向云圖分布的規(guī)律基本一致，但是網(wǎng)格自適應(yīng)加密后，云圖顯示的變形區(qū)域更為細(xì)致，能夠更明顯地表現(xiàn)出有限元模型中產(chǎn)生大變形和應(yīng)力集中的局部區(qū)域。

表2 杯突成型厚向云圖分布

BELYTSCHKO等[13]率先在其算法中采用單元應(yīng)變能的增量作為是否進(jìn)行局部網(wǎng)格加密的準(zhǔn)則來解決殼體大變形的計(jì)算問題，因此將本文算法與Belytschko算法分別在圖9(a)中的點(diǎn)(= 144 mm，=110 mm)和點(diǎn)(=150 mm，=56 mm)處的自適應(yīng)網(wǎng)格節(jié)點(diǎn)的位移-時(shí)間曲線進(jìn)行了對比，得到結(jié)果如圖10所示。

通過對比可以發(fā)現(xiàn)，使用本文算法計(jì)算的和節(jié)點(diǎn)的位移-時(shí)間曲線規(guī)律和Belytschko算法計(jì)算的結(jié)果基本符合，而在相同時(shí)間內(nèi)，采用本文的誤差準(zhǔn)則，可以更迅速地標(biāo)記應(yīng)力集中和應(yīng)變梯度大的區(qū)域，從而能夠加密更多的單元，得到更準(zhǔn)確的結(jié)果。

3.2 空調(diào)蓋板成型

與杯突成型相比，空調(diào)蓋板模型更為復(fù)雜，具有更多的成型特征，模型尺寸也更大，對整個(gè)模擬過程的精度要求更高。此外，該算例還具有實(shí)際工程意義，可驗(yàn)證算法的實(shí)用性與穩(wěn)定性。空調(diào)蓋板的計(jì)算模型如圖11所示。

根據(jù)空調(diào)蓋板的對稱性，選取整個(gè)有限元模型的1/2進(jìn)行自適應(yīng)分析，空調(diào)蓋板的仿真參數(shù)見表3。

如圖12所示，空調(diào)蓋板在沖壓成型過程中分別進(jìn)行了3次網(wǎng)格自適應(yīng)加密。為了驗(yàn)證自適應(yīng)分析算法的準(zhǔn)確性，本文做出如下的對比實(shí)驗(yàn)：

圖10 位移-時(shí)間曲線對比((a)A節(jié)點(diǎn)；(b)B節(jié)點(diǎn))

首先，分別設(shè)置與第1次至第3次自適應(yīng)加密網(wǎng)格大小相同、板料相同的3個(gè)有限元沖壓模型，每一個(gè)模型都劃分成與對應(yīng)自適應(yīng)次數(shù)網(wǎng)格相同大小的均勻網(wǎng)格；其次，在各個(gè)板料的同一位置(=496 mm，=283 mm)處設(shè)置測試點(diǎn)，如圖12(a)點(diǎn)所示，測量和記錄板料在該點(diǎn)的厚度變化信息；然后，將3個(gè)模型設(shè)置成除了網(wǎng)格尺寸和密度不同以外，其他條件均相同，利用已經(jīng)成熟應(yīng)用的商業(yè)軟件LS-DYNA求解；最后，將求解的結(jié)果與本文自適應(yīng)分析算法求解的結(jié)果進(jìn)行對比。

圖11 空調(diào)蓋板成型計(jì)算模型

表3 空調(diào)蓋板成型仿真參數(shù)

圖12 空調(diào)蓋板成型自適應(yīng)加密過程((a)初始網(wǎng)格；(b)第1次自適應(yīng)加密；(c)第2次自適應(yīng)加密；(d) 第3次自適應(yīng)加密；(e)最終成型)

空調(diào)蓋板成型在不同網(wǎng)格密度時(shí)的板厚情況如圖13所示，金屬板料與凸模產(chǎn)生接觸后，板料在壓力的作用下發(fā)生流動(dòng)，均出現(xiàn)減薄現(xiàn)象。在本文自適應(yīng)算法中板料的減薄情況與全部劃分為3次自適應(yīng)加密的最小尺寸網(wǎng)格的減薄情況幾乎一致，效果好于其他兩種均勻尺寸網(wǎng)格，顯示出本文自適應(yīng)算法的優(yōu)越性。

此外，為了測試阻尼子循環(huán)算法的精度，本文在空調(diào)蓋板成型過程中還增加了10個(gè)應(yīng)變測試點(diǎn)，分別測量位于該點(diǎn)的主應(yīng)變和次應(yīng)變，測試點(diǎn)位置見圖12(e)。

圖13 不同網(wǎng)格密度板厚對比

圖14分別采用了阻尼子循環(huán)、常速度子循環(huán)及無子循環(huán)算法，且主應(yīng)變和次應(yīng)變的總體趨勢基本一致，其中無子循環(huán)算法和阻尼子循環(huán)算法的誤差較小，常速度子循環(huán)算法的誤差較大。

圖14 主應(yīng)變和次應(yīng)變對比((a)主應(yīng)變；(b)次應(yīng)變)

圖15對有無阻尼子循環(huán)算法的計(jì)算時(shí)間進(jìn)行了對比，采用阻尼子循環(huán)法后，有限元仿真求解效率提升顯著，杯突成型計(jì)算時(shí)間減少了19.70%，空調(diào)蓋板成型計(jì)算時(shí)間減少了26.04%，對于復(fù)雜零件的沖壓成型過程，阻尼子循環(huán)算法的計(jì)算效率會(huì)更高。

圖15 計(jì)算時(shí)間對比

4 結(jié) 論

綜上所述，本文建立了一種基于自適應(yīng)分析模式下的薄板沖壓成型算法：

(1) 提出了基于單元應(yīng)變能增量的能量誤差準(zhǔn)則以及基于板料成型幾何特征的幾何誤差準(zhǔn)則，并將2類準(zhǔn)則進(jìn)行結(jié)合，建立了基于自適應(yīng)分析算法的誤差判斷準(zhǔn)則，提高有限元仿真精度；

(2)根據(jù)自適應(yīng)分析求解的特性，引入阻尼因子提出了阻尼子循環(huán)算法，將自適應(yīng)加密后的板料單元?jiǎng)澐譃椴煌瑫r(shí)間步長的區(qū)域，對每個(gè)區(qū)域分別進(jìn)行求解，在減少數(shù)值積分計(jì)算量、提高仿真計(jì)算效率的同時(shí)，增強(qiáng)了算法的穩(wěn)定性。

最后，杯突成型和空調(diào)蓋板成型的仿真實(shí)驗(yàn)結(jié)果驗(yàn)證了算法的有效性和實(shí)用性。

[1] BABU?KA I, RHEINBOLDT W C. A-posteriori error estimates for the finite element method[J]. International Journal for Numerical Methods in Engineering, 1978, 12(10): 1597-1615.

[2] 張征, 吳化平, 柴國鐘, 等. 自適應(yīng)無網(wǎng)格法在生物涂層接觸問題中的應(yīng)用[J]. 力學(xué)與實(shí)踐, 2012, 34(3): 58-61, 42.

ZHANG Z, WU H P, CHAI G Z, et al. Adaptive meshless method and its application in the contact problem of the bioactive coating[J]. Mechanics in Engineering, 2012, 34(3): 58-61, 42 (in Chinese).

[3] ZIENKIEWICZ O C, ZHU J Z. A simple error estimator and adaptive procedure for practical engineering analysis[J]. International Journal for Numerical Methods in Engineering, 1987, 24(2): 337-357.

[4] JIANG K F, YUAN S, XING Q Y. An adaptive nonlinear finite element analysis of minimal surface problem based on element energy projection technique[J]. Engineering Computations, 2020, 37(8): 2847-2869.

[5] CHEN Y Y, HUANG Y Q, YI N. Recovery type a posteriori error estimation of adaptive finite element method for Allen-Cahn equation[J]. Journal of Computational and Applied Mathematics, 2020, 369: 112574.

[6] WANG Y H, PENG X J, XIE M, et al. High-order lattice Boltzmann framework and its adaptive mesh refinement in the neutron transport SP3solutions[J]. Progress in Nuclear Energy, 2020, 128: 103449.

[7] XIAO X F, ZHAO J P, FENG X L. A layers capturing type H-adaptive finite element method for convection- diffusion-reaction equations on surfaces[J]. Computer Methods in Applied Mechanics and Engineering, 2020, 361: 112792.

[8] BABUSKA I, GUO B Q. Direct and inverse approximation theorems for the p-version of the finite element method in the framework of weighted besov spaces. part I: approximability of functions in the weighted besov spaces[J]. SIAM Journal on Numerical Analysis, 2002, 39(5): 1512-1538.

[9] 徐崗, 王毅剛, 胡維華. 等幾何分析中的r-p型細(xì)化方法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2011, 23(12): 2019-2024.

XU G, WANG Y G, HU W H. R-p-refinement in isogeometric analysis[J]. Journal of Computer-Aided Design & Computer Graphics, 2011, 23(12): 2019-2024 (in Chinese).

[10] TABARRAEI A, SUKUMAR N. Adaptive computations using material forces and residual-based error estimators on quadtree meshes[J]. Computer Methods in Applied Mechanics and Engineering, 2007, 196(25-28): 2657-2680.

[11] ZANDER N, BOG T, ELHADDAD M, et al. The multi-level hp-method for three-dimensional problems: Dynamically changing high-order mesh refinement with arbitrary hanging nodes[J]. Computer Methods in Applied Mechanics and Engineering, 2016, 310: 252-277.

[12] DEVLOO P, ODEN J T, STROUBOULIS T. Implementation of an adaptive refinement technique for the SUPG algorithm[J]. Computer Methods in Applied Mechanics and Engineering, 1987, 61(3): 339-358.

[13] BELYTSCHKO T, WONG B L, PLASKACZ E J. Fission- fusion adaptivity in finite elements for nonlinear dynamics of shells[J]. Computers & Structures, 1989, 33(5): 1307-1323.

Application in sheet metal forming based on-adaptive finite element method

LIU Xi-zhou, WANG Cheng-jing, WANG Hu

(State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha Hunan 410082, China)

In the finite element simulation analysis of the sheet metal forming process, it is difficult to accurately analyze the stress concentration and strain gradient areas where the field variables change drastically. How to balance the relationship between accuracy and efficiency becomes the key to stamping simulation. Therefore, based on the related theory of nonlinear finite element large deformation and the key technology of mesh adaptation for dynamic simulation, an algorithm for sheet metal forming under the adaptive analysis mode was established. In order to improve the calculation accuracy, an energy error criterion based on the unit strain energy increment and a geometric error criterion based on the geometric characteristics of sheet metal forming were proposed. Combining these two types of error criteria, an error judgment criterion based on an adaptive analysis algorithm was established.In order to enhance the calculation efficiency, by introducing a damper factor, a damper subcycling algorithm was proposed, which divided the adaptive sheet element unit into several regions according to the time step, and each region was integrated separately according to the time step.The results show that the algorithm improves the accuracy and efficiency of finite element simulation of sheet metal forming.

nonlinear finite element analysis; sheet metal forming; adaptive mesh refinement; error estimation; subcycling algorithm

TP 391

10.11996/JG.j.2095-302X.2021060970

A

2095-302X(2021)06-0970-09

2021-01-18；

2021-03-11

國家自然科學(xué)基金項(xiàng)目(11972155，51621004)

劉習(xí)洲(1994–)，男，安徽明光人，碩士研究生。主要研究方向?yàn)橛邢拊Y(jié)構(gòu)仿真。E-mail：liuxizhou@hnu.edu.cn

王 琥(1975–)，男，湖南長沙人，教授，博士。主要研究方向?yàn)橛邢拊Y(jié)構(gòu)仿真、快速計(jì)算方法等。E-mail：wanghu@hnu.edu.cn

11 March，2021

18 January，2021；

National Natural Science Foundation of China (11972155, 51621004)

LIU Xi-zhou (1994–), male, master student. His main research interest covers FEM structure simulation. E-mail：liuxizhou@hnu.edu.cn

WANG Hu (1975–), male, professor, Ph.D. His main research interests cover FEM structure simulation, fast calculation method, etc. E-mail：wanghu@hnu.edu.cn