有限元法中間構形初始解預示的 Laplace-Beltrami 方程法

劉永財，陳文亮，胡慶婉，鮑益東

(1. 南京航空航天大學機電學院，江蘇省精密與微細制造技術重點實驗室，南京 210016； 2. 曲靖師范學院數學與統計學院，曲靖 655011)

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

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

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

1 解耦格式

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

根據虛位移原理，在任一時刻，有：

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

式中：tF表示t時刻外部作用節點力；tR表示t時刻單元應力的節點力。

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

將剛度矩陣Kt看成由膜單元剛度矩陣和彎曲剛度矩陣兩部分組成的，因此，如圖1 所示，在節點所在的局部坐標系中，可以建立如下方程組來求解時刻tt+Δ 的中間構形[18]：

式中：1Δu和 2Δu為待求解節點的位移增量向量；0u為指定節點的位移向量；t t+ΔF為tt+Δ 時刻外力在切平面上引起的外部作用的摩擦力向量；tR1和tR2分別為t時刻指定節點的位移在法線方向上和切平面上所產生的力向量；tKmij和tKgij(i,j=1,2)分別為t時刻下的膜單元剛度矩陣和彎曲單元剛度矩陣。利用坐標系之間的轉換關系，對式(4)進行坐標轉換。首先，將方程組從節點所在的局部坐標系轉換到單元局部坐標系中，再從單元局部坐標系轉換到全局坐標系中，從而得到全局坐標系中的整體有限元方程組。在轉換的過程中，節點局部坐標系下的單元剛度矩陣被轉換為全局坐標系下的單元剛度矩陣，而膜單元矩陣中的元素值遠遠大于彎曲單元矩陣中的元素值，這將導致方程組高度病態，造成迭代無法快速收斂甚至迭代發散。

圖1 單元坐標系統 Fig.1 Element coordinate systems

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

式(5)對應薄板的彎曲變形，式(6)對應薄板的拉伸 變形。

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

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

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

可以看出，解耦格式在計算上能夠帶來兩方面優勢：它將一個大型方程組分解為兩個小型方程組；兩個小型方程組單元剛度矩陣中的元素值相差不大，條件數得到有效的改善，因此可以加快迭代算法的收斂速度。

綜上所述，中間構形被分解為彎曲變形和拉伸變形兩個計算過程。其中，以彎曲變形為初始解，通過彈塑性流動力學進行迭代求解來計算拉伸變形，該結果為一個中間構形。本文只研究計算彎曲變形，稱其為中間構形初始解。可以看出，一個高效的中間構形初始解構造方法十分必要。由于Laplace-Beltrami 方程具有堅實的數學理論和明確的物理意義，它已經被成功應用到數字幾何處理領域，受到了廣泛的關注和研究。下面將研究使用LB方程來構造中間構形初始解。

2 Laplace-Beltrami 方程

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

式中，G=det(g),(gij) =(gij)-1，而gij是兩個協變張量基做內積的結果。然而，對于一般的復雜曲面，則很難得到連續情形下LB 算子的顯式表達式。因此，在實際應用中，通常將待處理的復雜曲面表示為離散的多邊形網格形式，并在網格上構造離散LB算子的顯式表達式[21]。

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

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

本文選取基于平均曲率計算理論推導的LB 算子的余切表達式[22]為：

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

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

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

類似地，可以推導出二階LB 方程為：

一階LB 算子對應的是曲面第一基本形式的變化量，而二階LB 算子對應的是曲面第二基本形式的變化量，為了綜合考慮這兩種形式的變化量。一般地，將一階和二階LB 算子相結合獲得LB 方程來表達網格曲面變形，即：

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

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

3 約束條件及方程求解

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

式中：D為彈性矩陣；0σ為初應力，而應變和位移之間具有非線性關系，即：

式中：Lε和NLε分別為應變的線性部分和非線性部分；a為節點的位移向量。

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

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

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

式中，a是三角形三個節點的位移向量。

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

得到LB和NLB后，就可以通過式(16)來計算應變。接著，通過式(15)計算出每個單元的應力，進而得到每個單元所受的內力。根據靜力等效原則，通過單元內力能夠計算出每個節點處的等效節點力。基于等效節點力，可通過雙動情形來建立LB方程的約束條件。

遍歷所有與凸模曲面網格發生接觸的單元節點，根據等效節點力的方向和法線方向之間夾角大小來判斷節點是否為約束節點。

1) 如果兩者之間的夾角小于90°，力將節點拉離接觸面，當前構形的約束節點在下一個臨時構形中成為自由節點，則需要解除該節點的固定約束條件，如圖4 中的A點所示；

2) 如果兩者之間的夾角大于90°，力將節點繼續固定在接觸面上，則該節點在下一個臨時構形中仍為約束節點，如圖4 中的B點所示。

可以看出，隨著迭代次數的增加，自由節點的個數逐漸增加，約束節點的個數逐漸減少，直至臨時構形滿足相應的平衡條件。

圖4 約束節點和自由節點示意圖 Fig.4 Constraint and free nodes diagram

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

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

4 算例

基于解耦算法的多步有限元算法，已經被成功集成到自主開發的板料成形模擬軟件QuickForm中。其中，使用LB 方程來構造中間構形初始解。為驗證算法的有效性，這里以板料成形商業有限元軟件DynaForm 的模擬結果作為標準值，該軟件因使用動力顯式算法，因為其模擬精度高、可靠性好，已經被廣泛應用到覆蓋件的板料成形的仿真分析中。

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

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

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

另外，圖6 中還分別給出了使用張拉膜單元方法和LB 算子方法沿著截面線得到的中間構形初始解形狀。可以看出，兩種方法均可以獲得良好的中間構形初始解形狀，它們得到的結果非常接近。但是，采用相同的迭代終止條件時，兩種方法在計算當前步的中間構形初始解所需時間卻是不同的，使用張拉膜單元方法計算時間為24.32 s，而使用LB算子方法計算時間為16.68 s。顯然，LB 算子方法的計算速度更快，計算效率更高，這是因為它求解的線性方程組維數更少的緣故。

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

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

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

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

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

圖8 最終構形的外輪廓圖 Fig.8 Contours of final part

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

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

圖10 某復雜零件最終構形的外輪廓圖 Fig.10 Contours of final part for a complicated part

5 結論

在多步成形有限元法中，本文將中間構形的構造分解為拉伸變形和彎曲變形的解耦格式；根據大位移應變的理論，推導了中間構形的節點位移約束，給出了利用LB 方程求解中間構形的初始解預示算法。通過數值實驗結果表明：

(1) 在多步有限元法中采用解耦格式，增加了每步的行進步長，有效地改善了方程組的條件數，具有更好地數值求解穩定性和收斂性，提高了有限元法的計算效率。

(2) 基于大位移應變理論推導位移約束，能夠有效地模擬出彎曲變形過程中的節點位置變化移動情況。

(3) 通過離散的LB 算子線性方程來進行中間構形的初始解預示，算法實現簡單、通用性強、計算速度快。