剪裁CAD沖壓件的一步逆等幾何分析方法

張志功， 劉依娜， 張向奎， 王長生

(大連理工大學 汽車工程學院， 遼寧 大連 116024)

0 引 言

汽車制造過程需要將大量的金屬沖壓成形。為節省成本，沖壓成形數值模擬被廣泛應用于汽車設計中。板料沖壓成形數值模擬最常用的方法是基于流動理論的增量法和基于塑性變形理論的全量法。

增量法對板料進行沖壓成形仿真預測時，需要在每一增量步開始時輸入所有參數，雖然結果相對準確，但是計算相當耗時。因為增量法的缺點較多，GUO等首先提出Inverse Approach，該方法可以根據最終構型計算得到初始構型，LEE等提出單步逆和多步逆分析方法。逆成形法是從沖壓件的最終構型出發的，所以只需要輸入沖壓件的最終幾何形狀CAD模型即可，計算精度低于增量法，但是計算耗時大大減少。

可采用的數值計算方法有很多，例如有限差分法、有限元法和邊界元法等。增量法和逆成形方法一般都采用有限元法。先使用伽遼金等弱形式構建數學模型，再對CAD模型進行離散，即可求解所需要的未知量。然而，有限元法一般使用拉格朗日等基函數，由于拉格朗日基函數的插值性質，模型離散過程會改變模型的CAD數據，拉格朗日單元的前處理過程也相當耗時耗力。針對以上缺點，HUGHES等提出等幾何分析法，以解決有限元方法無法精確表達CAD模型的缺陷。一步逆成形法分可為初始解算法和塑性迭代算法。ZHANG等將IGA應用于一步逆成形方法中，開發板材沖壓一步逆成形等幾何分析初始解快速預示算法，其分析結果表明該方法在計算效率、求解精度和魯棒性等方面具有很好的效果。WANG等利用可塑性全變形理論，將一步逆成形算法與等幾何分析方法相結合，提出一步逆等幾何分析方法，該方法在不損失CAD數據的情況下具有很高的精度和準確性。WANG等借用TSA方法與Nitsche方法解決等幾何一步逆初始解算法無法計算復雜剪裁多片模型的限制，為進一步采用等幾何一步逆塑性迭代算法計算復雜多片剪裁模型奠定基礎。

根據張量積型NURBS曲面的特性，NURBS曲面表示復雜模型時需要進行裁剪。BEER等使用規則張量積曲面重構被剪裁曲面有效區域，此方法計算程序非常簡單，但是不能重構復雜的模型。XIA等使用三角形Bézier曲面轉換剪裁模型，將剪裁模型轉換成C0連續的小平面模型。WANG等使用準一致方法處理剪裁曲面，該方法的主要優點是可以直接采用單元的邊界曲線進行數值積分。KIM等提出TSA方法，其核心在于計算積分點，主要通過曲線插值和多次坐標變換完成。SCHMIDT等使用獨立曲面對剪裁區域進行局部重構，其核心在于曲面逼近，一般使用最小二乘法。

本文利用Bézier曲面重構剪裁單元并計算剪裁單元的剛度矩陣，然后通過轉換矩陣將其轉換到NURBS曲面上，以實現剪裁單元的處理，進而將此方法應用到等幾何一步逆法中。

1 NURBS曲面的剪裁算法

1.1 剪裁交點的搜索算法

一般使用B樣條對NURBS曲面進行剪裁。與計算機圖形學不同的是，等幾何分析需要求出B樣條和NURBS曲面在NURBS曲面參數空間中的交點，方便后續進行數值積分等工作。對于樣條曲線或者曲面來說，其擁有2個空間，一個是參數空間，一個是物理空間。B樣條曲線的空間為一維向量，參數為，而NURBS曲面的空間為二維平面，參數為(,)，與(,)相互獨立，不存在任何聯系。為使剪裁曲線B樣條與NURBS曲面處在同一個物理空間中，必須建立與(,)之間的聯系，即得到參數在物理空間的坐標(,,)，同時找到該坐標在NURBS曲面空間的(,)值，需要推導將曲線參數對應到物理空間坐標(,,)、再反演到NURBS曲面參數空間的標(,)的聯系。

求取剪裁交點的程序實現流程如下：

(1)對剪裁曲線進行節點插入細化，在參數空間中得到大量采樣點的值，將其命名為Line-t。

(2)采用NURBS曲面公式計算Line-t對應的物理空間的坐標(見圖1)，將其命名為Line-xyz。

圖1 細化后剪裁曲線節點向量映射到被裁剪曲面的參數點

(3)Line-xyz同在NURBS曲面上，利用NURBS曲面公式和Newton迭代法，將Line-反演到NURBS曲面上，得到相對應的NURBS參數空間的Line-和Line-。

(4)剪裁曲線節點向量在NURBS曲面的參數空間見圖2。綠色的方點為交點，紅色的圓點為采樣點，每一交點都被包含在(,)～(+1,+1)區間中，并且任意一個區間僅包含一個交點或者不包含交點，以保證正確求解所有交點。

圖2 剪裁曲線節點向量在NURBS曲面的參數空間示意

(5)根據(,)和(+1,+1)點，令

=|-+1||-+1|

(1)

=+|+1-|

(2)

即可以得到NURBS曲面參數+1對應的B樣條參數。

(6)利用NURBS反演函數計算相對應的(,)。

(7)比較與+1的差值是否在誤差范圍內：若差值在誤差范圍內，即認為求得交點；若與+1的差值不在在誤差范圍內，則繼續計算(8)。

(8)判斷(,)所在的網格區間，即是在+1左側還是右側：若是在左側，則將(,)賦值給(,)；若是在右側，則將(,)賦值給(+1,+1)。

(9)繼續循環(5)～(8)，得到誤差范圍內的所有交點，即認為交點求出。將得到的交點按照曲線節點向量大小排序，以便后續單元分類。

(10)若后續計算中發現得到的交點不全，則需要回到第一步，繼續加密剪裁曲線，直到正確求解所有交點且沒有遺漏。

1.2 剪裁單元分類

為方便剪裁單元的重構，需要人為地對剪裁單元進行分類，某剪裁后的單元分類示意見圖3。剪裁單元可按照形狀進行分類：紅色的三角形單元設定為類型A，綠色的四邊形單元設定為類型B，深藍色的五邊形單元為類型C，左側淺藍色網格為需要保留的單元，設定為類型S，右側沒有顏色的單元為丟棄的單元，設為類型D。

圖3 剪裁后的單元分類示意

由于Bézier曲線參數空間的節點向量是按照從小到大排列的，故最小節點位置通常被認為是曲線的起點，可以利用這個特點判斷要保留的單元。例如，可以人為設定剪裁曲線的左側為保留單元，右側單元不保留。常見的剪裁單元示意見圖4(單元的分類與前文一樣)。

圖4 常見的剪裁單元示意

此外，還有可能出現剪裁線穿過單元角點的情況，雖然在實際計算中這種情況非常罕見，但是編寫程序時要注意這種情況，以增強程序的魯棒性。

1.3 剪裁單元重構

為方便計算，使用Bézier曲面對剪裁單元進行重構，使用最小二乘法逼近重構曲面與被重構曲面。為使重構曲面的精度足夠高，重構曲面的階數應不低于剪裁曲面。

使用最小二乘法重構剪裁曲面，需要分別在2個曲面的參數空間內選取配點，配點要盡量均勻，同時因為Bernstein基函數與NURBS基函數具有邊界插值性質，所以配點必須在曲面邊界上選取。Bézier基函數描述的曲面與NURBS基函數描述的曲面是等價的，即

∑(,t,,t)=∑(,)

(3)

則個配點的等價公式為

·==·

(4)

其中，

(5)

(6)

(7)

(8)

式中：為剪裁NURBS單元的基函數矩陣；(,)為之前選取的配點所在參數點的基函數值；為剪裁NURBS單元的控制頂點；為重構曲面的Bernstein基函數矩陣；為Bézier曲面的控制頂點，是整個方程中的未知量。

利用最小二乘法計算，

=(·)··()·

(9)

由此可知，用于重構的Bézier曲面的控制點可以通過原始NURBS曲面的控制點轉換而來，即

=·

(10)

式中：為轉換矩陣，

=(·)··()

(11)

2 剪裁NURBS曲面與一步逆成形法的耦合

2.1 剪裁一步逆等幾何分析法初始解預示

剪裁單元重構后，初始解部分的求解與未重構部分相同，都有完整的基函數支撐，區別在于有一部分單元不再是部件的一部分，所以不再對其進行求解。

根據一步逆初始解算法，重構后的單元在局部坐標系下初始單元節點的位移可以表示為

(12)

單元的應力可以表示為

(13)

最后，得到整體坐標系下的節點內力為

(14)

迭代公式為

(15)

式中：Δ為節點位移增量；為松弛因子。是根據能量因子決定的，0<<1。殘余力矢量與Δ相乘即可求得能量因子。一般來說，能量因子越大，越小。

迭代上述公式，直到內力達到平衡狀態，使用位移節點向量的范數判斷計算收斂，

(16)

式中：為板料的總節點數。

2.2 剪裁一步逆等幾何分析塑性迭代求解

塑性迭代可計算材料沖壓時的應力、應變和厚度的變化，剛度矩陣計算參照文獻[19]，根據Newton-Raphson迭代法求解非線性方程，即

(17)

式中：為迭代步時的殘余力；()為節點外力；()為節點內力。

位移場更新公式為

+1=+Δ

(18)

因為一步逆成形算法中沒有施加邊界條件，所以添加剛度阻尼(>0)使剛度矩陣非奇異，即

,+Δ=

(19)

為加快迭代收斂速度，可使剛度阻尼系數取值較小。

3 數值算例

3.1 剪裁S梁算例

計算一個剪裁的S梁模型(見圖5)，此模型是NUMISHEET96的標準考題。原始S梁模型為雙二次張量積型NURBS曲面，共有控制點784個、單元676個。為驗證算法的準確性，只裁掉S梁的一小部分，非剪裁單元保留624個，剪裁單元個數為26個。用2種不同方法計算出S梁的厚度云圖和等效應力云圖，結果分別見圖6和7。2種方法求出的厚度與等效應力分布幾乎相同，等幾何方法求出的結果云圖更平滑一些。2種方法得到的厚度與等效應力的最大值、最小值以及二者的相對誤差見表1，厚度與等效應力結果相對誤差都在可接受的范圍內。

圖5 剪裁S梁示意

圖7 S梁等效應力云圖對比, MPa

表 1 S梁算例有限元法與等幾何法結果對比

4.2 引擎蓋算例

汽車引擎蓋是汽車最具有代表性的沖壓件之一。同時使用剪裁一步逆等幾何方法和有限元法計算一個引擎蓋算例，見圖8。該引擎蓋未剪裁的雙CAD數據為二次NURBS曲面，共有450個控制頂點、364個單元，剪裁后共有保留的剪裁單元為298個，刪除剪裁單元30個。2種方法計算得到的引擎蓋的厚度云圖和應力云圖分別見圖9和10，可見2種方法計算結果的一致性較高。

圖8 剪裁引擎蓋NURBS曲面和控制網格

圖9 引擎蓋厚度云圖對比, mm

圖10 引擎蓋等效應力云圖對比, MPa

2種方法計算得到的厚度和等效應力的最大值、最小值以及二者的相對誤差見表2。對于此算例，厚度和等效應力的相對誤差極小。

表 2 汽車引擎蓋算例有限元法與等幾何法計算結果對比

4 結束語

將剪裁算法和等幾何分析方法應用于基于一步逆方法的板材沖壓成形分析中，構建基于塑性形變理論的剪裁一步逆成形等幾何法。首先，預處理被剪裁曲面，包括求取剪裁交點、尋找剪裁單元、剪裁單元分類、剪裁單元重構、剪裁單元組裝單元剛度矩陣等；然后，利用初始解算法，展開NURBS模型的控制點網格；最后，利用Newton-Raphson迭代算法處理最終構建的非線性平衡方程，得到最終構型的厚度、等效應力和應變等物理量的分布。

選擇不同的模型算例，分別在基于剪裁NURBS曲面的一步逆等幾何方法程序和一步逆有限元法程序中進行計算，對厚度和等效應力計算結果及二者的相對誤差，驗證基于剪裁NURBS曲面的一步逆等幾何分析法的正確性和準確性。