三維網(wǎng)格重劃分技術(shù)及其在焊接數(shù)值模擬中的應用

黃輝，趙耀，袁華，王波

（華中科技大學船舶與海洋工程學院，湖北武漢 430074）

0 引言

焊接是現(xiàn)代機械工業(yè)中的重要連接工藝，是結(jié)構(gòu)產(chǎn)生殘余變形和應力的主要來源。從20世紀50年代開始，人們開始焊接變形和應力的數(shù)值計算，取得了一系列的研究成果。近10多年來，隨著計算機的迅猛發(fā)展，使焊接過程的熱彈塑性計算成為模擬焊接變形、應力、組織轉(zhuǎn)變等復雜現(xiàn)象的主要手段。焊接過程中材料在高溫條件下發(fā)生分解、相變，而且不同溫度下的材料性能具有很大差異，使焊接局部區(qū)域經(jīng)歷高度非線性的過程，要模擬整個結(jié)構(gòu)的焊接變形和殘余應力，需要進行熱機耦合非線性有限元分析。由于焊縫區(qū)溫度和應力應變的梯度很大，進行有限元建模分析時，網(wǎng)格劃分必須足夠精細，并且進行時間增量步的分析。對于大型復雜結(jié)構(gòu)，其焊縫長度和焊接時間都很長，有限元網(wǎng)格規(guī)模和計算時間步勢必急劇膨脹，進行完整的熱彈塑性有限元分析將十分困難。

網(wǎng)格重劃分技術(shù)可以處理網(wǎng)格畸變及單元變形過大的問題，使熱軋、巖土等大變形分析得以進行下去，并保證計算結(jié)果的精度和可靠性。S．Brown和H．Song［1］結(jié)合網(wǎng)格重劃分和動態(tài)子結(jié)構(gòu)技術(shù)，采用殼體單元分析了平板激光焊接后的變形和應力，所使用的時間與完整的熱彈塑性計算相比，僅為后者的七分之一。L．－E．Lindgren［2］基于一種可變節(jié)點的體單元，開發(fā)了自適應網(wǎng)格技術(shù)，并應用三維有限元方法成功模擬了大型銅罐的電子束焊接，大幅縮短了計算的時間。本文開展了三維全六面體網(wǎng)格重劃分技術(shù)研究并將其應用于焊接數(shù)值模擬，對平板和圓筒的焊接模型分別進行了焊接熱力學仿真。

1 網(wǎng)格重劃分

焊接過程中，熱影響區(qū)的溫度、應力應變梯度較大，具有高度的非線性特征，有限元模擬中單元應當劃分得很細;然而，遠離焊縫的絕大部分區(qū)域仍然處于線彈性狀態(tài)，單元可以劃分得相對較粗。如果不斷對焊接模型進行網(wǎng)格重新劃分，使細密的網(wǎng)格隨著熱源一起移動，就能節(jié)省大量的計算資源并保證計算結(jié)果的精度。在網(wǎng)格重劃分前后，必須保證新舊網(wǎng)格的節(jié)點占據(jù)相同幾何空間，以順利地進行應力應變等變量的映射，使分析求解得以連續(xù)進行下去。下面介紹節(jié)點遷移和變量映射2個步驟。

1.1 節(jié)點遷移

由于網(wǎng)格疏密程度在網(wǎng)格重劃分前后變化較大，不得不完全改變節(jié)點和單元的編號，重新構(gòu)建有限元模型的網(wǎng)格。為了匹配新網(wǎng)格節(jié)點在舊網(wǎng)格中的位置，必須對網(wǎng)格單元一一尋址，逐一判斷節(jié)點是否處于某單元內(nèi)部［3］。節(jié)點遷移的基本問題就是已知單元8個節(jié)點的坐標和位于該單元中1個節(jié)點的坐標，反求這1個節(jié)點在該單元對應的等參元中的坐標。

對于六面體單元，其形函數(shù)為:

根據(jù)等參元的性質(zhì)，單元中任一點處的整體坐標

將表1中坐標值代入式（1）和式（2），可得如下方程組:

式中:ai，bi，ci為已知量xi，yi，zi的函數(shù)，它們的值均為常數(shù)。上述高階齊次方程組的求解方法有牛頓法、割線法、布朗方法等，本研究采用了牛頓法，對于比較規(guī)則的網(wǎng)格單元，進行為數(shù)不多的幾次迭代其結(jié)果就能收斂。求得節(jié)點在局部坐標系中的坐標后，可以用形函數(shù)對節(jié)點位移進行插值，從而得到下一步分析時的節(jié)點在整體坐標系中的坐標值。

1.2 變量映射

確定新網(wǎng)格的節(jié)點坐標后，需要將舊網(wǎng)格中求解的應力應變場映射到新網(wǎng)格上，以保證分析的連續(xù)性和準確性。采用了如下映射算法［4］:

1）采用插值技術(shù)，將上一步網(wǎng)格中積分點上的變量（應力、應變等）外插到每一個單元的節(jié)點上，然后對所有鄰接同一節(jié)點的單元的變量取平均。

2）新網(wǎng)格的積分點的位置由舊網(wǎng)格積分點來確定，分兩步:①先找到該點所屬舊網(wǎng)格的單元，這樣就可以確定該點在單元中的位置（該步驟假定新網(wǎng)格中所有積分點位于舊網(wǎng)格邊界內(nèi)）。②變量由舊網(wǎng)格的節(jié)點上內(nèi)插到新網(wǎng)格的積分點上。

2 計算實例

分析了平板表面堆焊和圓筒環(huán)焊縫的焊接過程，分別比較了完整模型與網(wǎng)格重劃分模型的計算結(jié)果。平板堆焊模型的尺寸為500 mm×300 mm×6 mm（圖1）。圓筒的內(nèi)徑為150 mm，外徑為156 mm，長度為400 mm（圖2）。平板和圓筒的焊接條件見表2。SS400的材料熱物性與熱力學性能參數(shù)曲線見文獻［5］。

2.1 熱傳導計算

采用順次耦合的熱彈塑性有限元計算方法，分別對平板模型和圓筒模型進行焊接熱力學仿真分析。首先進行焊接過程的熱傳導計算，把得到模型的瞬態(tài)溫度數(shù)據(jù)信息存儲在1個文件中，然后讀入每一時刻的溫度作為模型的載荷，進行焊接模型的應力應變計算。由于高斯熱源模型能很好地反映大多數(shù)焊接過程中的熱流分布特征，其在焊接熱傳導的數(shù)值仿真中應用甚廣。高斯熱流分布示意如圖4所示，函數(shù)表達式為:

圖3 焊縫中線處的圓筒橫剖面Fig.3Cross section of pipe model

式中:q（r）為半徑r處的表面熱流，W/m2;q（O）為熱源中心處的熱流量最大值;c為熱源集中系數(shù);r為距熱源中心的距離。

圖4 高斯熱源熱流分布示意圖Fig.4Heat flux distribution of Gauss heat source

熱傳導計算中考慮了相變潛熱，取為270 kJ/kg，固相線溫度為1 465℃，液相線溫度為1 544℃;環(huán)境溫度設(shè)為25℃，冷卻條件為空冷，表面換熱系數(shù)取常數(shù)10 W/（m2·℃）;考慮輻射散熱，發(fā)射率取為0.5。加熱過程中時間步長最大取0.5 s，冷卻過程取自動時間步長。為了兼顧計算的效率和精度，在焊縫附近區(qū)域網(wǎng)格劃分得很細，在遠離焊縫區(qū)劃分得較粗。沿焊縫方向根據(jù)熱源的位置，在網(wǎng)格重劃分模型中分4步實現(xiàn)焊接過程模擬，平板焊接模型的有限元網(wǎng)格如圖5所示。完整網(wǎng)格模型的有限元網(wǎng)格與網(wǎng)格重劃分模型第4步的網(wǎng)格相同，沿縱向均勻劃分，沿橫向由密變疏。圓筒焊接模型的網(wǎng)格劃分方法與重劃分步驟與平板焊接模型相似，在此不贅述。為了簡化，以FM（Full Model）表示完整模型，RM（Rezoning Model）表示網(wǎng)格重劃分模型。平板完整模型、平板網(wǎng)格重劃分模型、圓筒完整模型、圓筒重劃分模型分別表示為FMⅠ，RMⅠ，F(xiàn)MⅡ，RMⅡ。平板和圓筒焊接過程中的瞬時溫度場分布云圖如圖5和圖6所示。

2.2 應力應變計算

應力應變計算所使用的有限元模型與熱傳導計算模型的節(jié)點和單元完全相同，單元類型由熱單元改為力學單元，材料特性由熱物性參數(shù)變?yōu)闊崃W性能參數(shù)，同時增加了力學邊界。在平板有限元模型中只限制了結(jié)構(gòu)的剛體位移，在圓筒有限元模型中設(shè)置了一端簡支的約束條件。將熱傳導計算的節(jié)點瞬時溫度讀入有限元模型，開始應力應變的增量步計算。圖7為圓筒焊接后的Mises應力分布云圖。

圖7 圓筒焊后Mises應力云圖Fig.7Mises Stress distribution of pipe model after welding

3 計算結(jié)果

3.1 溫度場計算結(jié)果

平板焊縫長度中心處橫截面上3點的瞬時溫度結(jié)果見圖8，測點分布于平板表面，與焊縫的距離分別為0 mm，20 mm和44 mm。從圖中不難看出，焊接過程中各點經(jīng)歷了加熱和冷卻2個階段。距離焊縫越近的點，溫度上升過程越快，達到的最高溫度也越高。在冷卻階段，各點的溫度逐漸趨于一致。網(wǎng)格重劃分模型與完整模型的計算結(jié)果完全一致，說明采用網(wǎng)格重劃分技術(shù)在計算熱傳導過程時滿足精度要求。

3.2 變形場計算結(jié)果

平板焊接模型的變形結(jié)果見圖9，曲線表示了橫向收縮、角變形沿焊縫方向的分布情況。從圖中可以看出，橫向收縮在焊縫長度中間段有最大值，角變形隨焊縫長度增加略有變大，但整體上看，2種變形沿焊縫方向都較為均勻。

=圓筒焊接模型中A和B兩點（標示于圖3）的瞬態(tài)徑向位移示于圖10。從圖中可以明顯看出，當點所在區(qū)域處于加熱階段時，焊縫金屬發(fā)生膨脹，徑向位移變大，在熱源遠離該處，焊縫因冷卻而收縮，這時徑向位移減小。由于A，B兩點在環(huán)焊縫中受熱2次，可在位移圖中看到2個波峰。

由圖9和圖10可知，平板和圓筒的完整模型與網(wǎng)格重劃分模型計算結(jié)果相一致，說明網(wǎng)格重劃分模型預測的焊接變形結(jié)果沒有精度損失。

3.3 應力場計算結(jié)果

平板表面堆焊模型的殘余應力結(jié)果以等值線的方式繪制見圖11和圖12。由圖11可見，縱向殘余應力的分布情況，焊縫區(qū)域存在較大的拉應力，遠離焊縫區(qū)則呈現(xiàn)較低的壓應力以使構(gòu)件保持自身的平衡。圖12為平板焊后的橫向殘余應力結(jié)果，在焊縫兩端存在較大的壓應力，焊縫區(qū)的橫向應力水平較低。圖11和圖12表明，完整計算模型與網(wǎng)格重劃分模型計算的焊接殘余應力結(jié)果非常吻合。

圓筒焊接模型的0°和180°縱向剖面的殘余應力結(jié)果見圖13和14，圖中反映了圓筒焊后內(nèi)外側(cè)的應力分布情況。焊縫起始位置處（0°）的環(huán)向應力峰值較低，而焊縫中間段（180°）的環(huán)向應力峰值則超過了材料的屈服強度;同時，圓筒內(nèi)側(cè)的應力要高于外側(cè)的應力水平。圖14中的結(jié)果表明，圓筒內(nèi)側(cè)的軸向應力為拉應力，外側(cè)的軸向應力為壓應力，0°位置處的應力峰值低于180°位置處的應力水平。從圖中結(jié)果可以看出，完整模型和網(wǎng)格重劃分模型的橫向殘余應力和縱向殘余應力完全吻合，而且有限元分析預測的環(huán)向應力和軸向應力的分布規(guī)律與Burdekin［6］關(guān)于低碳鋼圓筒環(huán)焊縫的研究結(jié)果相一致。

3.4 有限元計算時間

平板與圓筒焊接模型的有限元網(wǎng)格節(jié)點數(shù)和計算時間如表3所示，2種模型在使用網(wǎng)格重劃分技術(shù)后，計算時間均減少了約30%。由于網(wǎng)格重劃分模型中的節(jié)點數(shù)一般少于完整模型的節(jié)點數(shù)，也即模型中總的自由度數(shù)前者少于后者，內(nèi)存的需求也必將因網(wǎng)格重劃分技術(shù)而得以減少。對大型結(jié)構(gòu)而言，計算時間和內(nèi)存大小是決定數(shù)值模擬能否實現(xiàn)的主要因素。計算結(jié)果表明，網(wǎng)格重劃分技術(shù)可以解決這方面的問題。本文暫時沒有研究熱源后端的網(wǎng)格粗化問題，若能選擇合理的粗化準則并優(yōu)化網(wǎng)格設(shè)計將可能大幅提高計算效率，同時保證結(jié)果的可靠性。

4 結(jié)語

1）焊接熱源影響下的區(qū)域具有高度非線性特征，遠離熱源區(qū)域為弱非線性，可以采用網(wǎng)格重劃分技術(shù)進行模擬。

2）使焊縫區(qū)細密網(wǎng)格隨熱源不斷移動，與完整的有限元模型相比，減少了節(jié)點自由度，提高了計算效率并保證了應力和變形結(jié)果的精度。

3）通過網(wǎng)格重劃分模型與完整模型的計算結(jié)果可知，提出的網(wǎng)格重劃分方法在充分保證結(jié)果精度的前提下，大約減少了30%的計算時間。

［1］BROWN SB，SONGH．Rezoninganddynamic substructuring techniques in FEM simulations of welding processes［J］．ASME Journal of Engineering for Industry，1993，155:415－423．

［2］LINDGREN L E，HAGGBLAD H A，et al．Automatic remeshing for three-dimensional finite element simulation 0f welding［J］．Computer Methods in Applied Mechanics and Engineering，1997，147（3）:401－409．

［3］MURTI V，VALLIAPPANS．Numericalinverseisoparametricmappinginremeshingandnodalquantity contouring［J］．Computers and Structures，1986，22（6）: 1011－1021．

［4］HIBBITT D，KARLSONB，SORENSONP．ABAQUS analysisuser'smanual［M］．ABAQUSVersion6.5 Documentation．

［5］WANG Rui，SHERIF R，et al．Numerical and experimental investigations on welding deformation［J］．Trans．JWRI，2008，37（1）:79－90．

［6］BURDEKIN F M．Local stress relief of circumferential butt welds in cylinders［J］．British Welding Journal，1963，10 （9）:483－490．