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

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

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

(湖南大學(xué)汽車車身先進設(shè)計制造國家重點實驗室，湖南 長沙 410082)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在隨動坐標系中計算每個單元的應(yīng)變能增量為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

阻尼因子一般取0.005～0.010。

3 算例分析

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

3.1 杯突成型

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

圖8 杯突成型計算模型

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

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

圖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ū)域更為細致，能夠更明顯地表現(xiàn)出有限元模型中產(chǎn)生大變形和應(yīng)力集中的局部區(qū)域。

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

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

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

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

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

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

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

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

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

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

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

此外，為了測試阻尼子循環(huán)算法的精度，本文在空調(diào)蓋板成型過程中還增加了10個應(yīng)變測試點，分別測量位于該點的主應(yīng)變和次應(yīng)變，測試點位置見圖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)算法的計算時間進行了對比，采用阻尼子循環(huán)法后，有限元仿真求解效率提升顯著，杯突成型計算時間減少了19.70%，空調(diào)蓋板成型計算時間減少了26.04%，對于復(fù)雜零件的沖壓成型過程，阻尼子循環(huán)算法的計算效率會更高。

圖15 計算時間對比

4 結(jié) 論

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

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

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

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

[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é)與實踐, 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型細化方法[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 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é)基金項目(11972155，51621004)

劉習(xí)洲(1994–)，男，安徽明光人，碩士研究生。主要研究方向為有限元結(jié)構(gòu)仿真。E-mail：liuxizhou@hnu.edu.cn

王 琥(1975–)，男，湖南長沙人，教授，博士。主要研究方向為有限元結(jié)構(gòu)仿真、快速計算方法等。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