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

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

0 引 言

汽車制造過程需要將大量的金屬?zèng)_壓成形。為節(jié)省成本，沖壓成形數(shù)值模擬被廣泛應(yīng)用于汽車設(shè)計(jì)中。板料沖壓成形數(shù)值模擬最常用的方法是基于流動(dòng)理論的增量法和基于塑性變形理論的全量法。

增量法對板料進(jìn)行沖壓成形仿真預(yù)測時(shí)，需要在每一增量步開始時(shí)輸入所有參數(shù)，雖然結(jié)果相對準(zhǔn)確，但是計(jì)算相當(dāng)耗時(shí)。因?yàn)樵隽糠ǖ娜秉c(diǎn)較多，GUO等首先提出Inverse Approach，該方法可以根據(jù)最終構(gòu)型計(jì)算得到初始構(gòu)型，LEE等提出單步逆和多步逆分析方法。逆成形法是從沖壓件的最終構(gòu)型出發(fā)的，所以只需要輸入沖壓件的最終幾何形狀CAD模型即可，計(jì)算精度低于增量法，但是計(jì)算耗時(shí)大大減少。

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

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

本文利用Bézier曲面重構(gòu)剪裁單元并計(jì)算剪裁單元的剛度矩陣，然后通過轉(zhuǎn)換矩陣將其轉(zhuǎn)換到NURBS曲面上，以實(shí)現(xiàn)剪裁單元的處理，進(jìn)而將此方法應(yīng)用到等幾何一步逆法中。

1 NURBS曲面的剪裁算法

1.1 剪裁交點(diǎn)的搜索算法

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

求取剪裁交點(diǎn)的程序?qū)崿F(xiàn)流程如下：

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

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

圖1 細(xì)化后剪裁曲線節(jié)點(diǎn)向量映射到被裁剪曲面的參數(shù)點(diǎn)

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

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

圖2 剪裁曲線節(jié)點(diǎn)向量在NURBS曲面的參數(shù)空間示意

(5)根據(jù)(,)和(+1,+1)點(diǎn)，令

=|-+1||-+1|

(1)

=+|+1-|

(2)

即可以得到NURBS曲面參數(shù)+1對應(yīng)的B樣條參數(shù)。

(6)利用NURBS反演函數(shù)計(jì)算相對應(yīng)的(,)。

(7)比較與+1的差值是否在誤差范圍內(nèi)：若差值在誤差范圍內(nèi)，即認(rèn)為求得交點(diǎn)；若與+1的差值不在在誤差范圍內(nèi)，則繼續(xù)計(jì)算(8)。

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

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

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

1.2 剪裁單元分類

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

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

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

圖4 常見的剪裁單元示意

此外，還有可能出現(xiàn)剪裁線穿過單元角點(diǎn)的情況，雖然在實(shí)際計(jì)算中這種情況非常罕見，但是編寫程序時(shí)要注意這種情況，以增強(qiáng)程序的魯棒性。

1.3 剪裁單元重構(gòu)

為方便計(jì)算，使用Bézier曲面對剪裁單元進(jìn)行重構(gòu)，使用最小二乘法逼近重構(gòu)曲面與被重構(gòu)曲面。為使重構(gòu)曲面的精度足夠高，重構(gòu)曲面的階數(shù)應(yīng)不低于剪裁曲面。

使用最小二乘法重構(gòu)剪裁曲面，需要分別在2個(gè)曲面的參數(shù)空間內(nèi)選取配點(diǎn)，配點(diǎn)要盡量均勻，同時(shí)因?yàn)锽ernstein基函數(shù)與NURBS基函數(shù)具有邊界插值性質(zhì)，所以配點(diǎn)必須在曲面邊界上選取。Bézier基函數(shù)描述的曲面與NURBS基函數(shù)描述的曲面是等價(jià)的，即

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

(3)

則個(gè)配點(diǎn)的等價(jià)公式為

·==·

(4)

其中，

(5)

(6)

(7)

(8)

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

利用最小二乘法計(jì)算，

=(·)··()·

(9)

由此可知，用于重構(gòu)的Bézier曲面的控制點(diǎn)可以通過原始NURBS曲面的控制點(diǎn)轉(zhuǎn)換而來，即

=·

(10)

式中：為轉(zhuǎn)換矩陣，

=(·)··()

(11)

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

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

剪裁單元重構(gòu)后，初始解部分的求解與未重構(gòu)部分相同，都有完整的基函數(shù)支撐，區(qū)別在于有一部分單元不再是部件的一部分，所以不再對其進(jìn)行求解。

根據(jù)一步逆初始解算法，重構(gòu)后的單元在局部坐標(biāo)系下初始單元節(jié)點(diǎn)的位移可以表示為

(12)

單元的應(yīng)力可以表示為

(13)

最后，得到整體坐標(biāo)系下的節(jié)點(diǎn)內(nèi)力為

(14)

迭代公式為

(15)

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

迭代上述公式，直到內(nèi)力達(dá)到平衡狀態(tài)，使用位移節(jié)點(diǎn)向量的范數(shù)判斷計(jì)算收斂，

(16)

式中：為板料的總節(jié)點(diǎn)數(shù)。

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

塑性迭代可計(jì)算材料沖壓時(shí)的應(yīng)力、應(yīng)變和厚度的變化，剛度矩陣計(jì)算參照文獻(xiàn)[19]，根據(jù)Newton-Raphson迭代法求解非線性方程，即

(17)

式中：為迭代步時(shí)的殘余力；()為節(jié)點(diǎn)外力；()為節(jié)點(diǎn)內(nèi)力。

位移場更新公式為

+1=+Δ

(18)

因?yàn)橐徊侥娉尚嗡惴ㄖ袥]有施加邊界條件，所以添加剛度阻尼(>0)使剛度矩陣非奇異，即

,+Δ=

(19)

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

3 數(shù)值算例

3.1 剪裁S梁算例

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

圖5 剪裁S梁示意

圖7 S梁等效應(yīng)力云圖對比, MPa

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

4.2 引擎蓋算例

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

圖8 剪裁引擎蓋NURBS曲面和控制網(wǎng)格

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

圖10 引擎蓋等效應(yīng)力云圖對比, MPa

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

表 2 汽車引擎蓋算例有限元法與等幾何法計(jì)算結(jié)果對比

4 結(jié)束語

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

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