三維網格重劃分技術及其在焊接數值模擬中的應用

黃輝，趙耀，袁華，王波

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

0 引言

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

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

1 網格重劃分

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

1.1 節點遷移

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

對于六面體單元，其形函數為:

根據等參元的性質，單元中任一點處的整體坐標

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

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

1.2 變量映射

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

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

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

2 計算實例

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

2.1 熱傳導計算

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

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

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

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

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

2.2 應力應變計算

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

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

3 計算結果

3.1 溫度場計算結果

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

3.2 變形場計算結果

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

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

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

3.3 應力場計算結果

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

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

3.4 有限元計算時間

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

4 結語

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

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

3）通過網格重劃分模型與完整模型的計算結果可知，提出的網格重劃分方法在充分保證結果精度的前提下，大約減少了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．