關于ABAQUS熱力耦合分析，經常會有一些問題匯總到這里，微辰三維今天就基于有限元軟件ABAQUS的順序熱力耦合分析進行平板焊接過程的數(shù)值模擬，希望能幫大家解決熱力耦合的問題。順序熱力耦合分析是ABAQUS中常用的熱-應力分析類型，這種分析類型適用于應力依賴于溫度場，而溫度場受應力影響較弱的情況。在分析過程中，需要完成兩次分析計算：首先進行熱傳導分析，獲得整個結構的溫度分布；隨后將計算得到的溫度以熱載荷的形式施加到結構上，并進行應力分析，獲得結構的應力分布。由于應力場和溫度場之間互不影響，因此相比于完全耦合的熱力分析，順序熱力耦合分析求解更為高效。

1. 平板焊接模型及焊接工藝參數(shù)

平板焊接采用的網格模型如圖1.1所示。平板的厚度為3mm，長寬均為100mm，由于分析具有對稱性，因此計算時采用1/2模型進行計算。為了保證計算精度，一般要求在焊接熔池長度內不能少于四個網格，因此焊縫區(qū)域的網格尺寸均小于1mm，而遠離焊縫的區(qū)域采用逐漸增大的網格進行過渡。

數(shù)值模擬選取的焊接工藝參數(shù)為：焊接總功率Q=1000W，焊接效率η=0.57，焊接速度v=5mm/s。

2. 焊接材料

本次焊接使用的材料為IN718。在熱傳導分析中，需要用到的材料參數(shù)有：導熱系數(shù)，比熱容，密度；在應力分析中，需要用到的材料參數(shù)有：彈性模量，泊松比，平均熱膨脹系數(shù)（線膨脹系數(shù)）和屈服強度。IN718相應的材料數(shù)據如圖2.1至圖2.3所示。

圖2.1 IN718導熱系數(shù)及密度隨溫度變化曲線

圖2.2 IN718比熱容及平均熱膨脹系數(shù)隨溫度變化

圖2.3 IN718楊氏模量及屈服應力隨溫度變化

圖2.3中給出了焊縫和母材對應的楊氏模量及屈服強度，本次計算時均取為母材對應的材料參數(shù)。盡管材料的泊松比也會隨著溫度變化，但將泊松比考慮為依賴于溫度的參數(shù)會增加計算收斂的難度，因此將IN718的泊松比取為0.3。

3. 焊接熱源模型

選取合適的熱源模型，對平板焊接殘余應力的計算至關重要。焊接熱源需要根據焊接類型進行選取，一般常用的焊接熱源有雙橢球熱源、高斯錐形熱源、均勻體熱源等。為了簡化起見，這里本文使用一個圓形面熱源進行計算，熱源的熱流密度分布可表示為：

其中η為焊接效率，Q為焊接總功率，r為圓形熱源的半徑，這里取為4mm。

在ABAQUS中，上述熱源可以通過用戶子程序DFLUX施加到有限元模型上。具體的計算原理這里不再贅述，可以參考基于DFLUX的焊接雙橢球熱源模擬，這里僅給出用到的Fortran程序。

需要注意的是，在上面程序中假設焊接熱源沿X的正向移動，并且初始時刻焊接熱源中心的坐標為(-r,0,0)，因此有限元網格也應與程序定義的方向保持一致。

4. 平板焊接溫度場計算

將有限元模型導入到ABAQUS中，定義材料IN718，并輸入計算需要用到的密度、導熱系數(shù)及比熱容，如圖4.1所示。

圖4.1 定義材料屬性

在Assembly模塊中插入平板對應的Part，如圖4.2所示。

圖4.2 平板有限元模型

注意圖4.2中全局坐標的中心與焊接線的起始節(jié)點是重合的，并且X正半軸指向焊接熱源移動的方向，這樣可以保證焊接熱源沿平板的中心移動。

在Step模塊中定義兩個分析步，分析類型為Heat transfer (Transient)，第一個分析步的計算時間為21s，由于平板長度為100m，而焊接速度為5mm/s，因此焊接熱源完全穿過平板所需的時間約為21s；第二個分析步的計算時間為1000s，用于計算焊接完成后的冷卻過程，如圖4.3所示。

圖4.3 創(chuàng)建熱傳導分析步

在第一個分析步中設定初始增量為0.1，最大增量為5，最小增量為1E-5，每個增量步允許的最大溫升設為500oC，如圖4.4所示。

圖4.4 熱傳導分析步的增量步設置

冷卻過程的熱傳導分析步中取初始增量步為0.001，最大增量為1000，最小增量為1E-5，每個增量步允許的最大溫升設為500oC。

在Interaction模塊中創(chuàng)建對流邊界條件Surfacefilm condition，對流換熱面選擇為除對稱面以外的所有單元面，如圖4.5所示。

圖4.5 選擇對流換熱面

根據牛頓冷卻定律，溫度高于周圍環(huán)境的物體向周圍媒介傳遞的熱量可表示為：

其中h為對流換熱系數(shù)，Tbc為周圍環(huán)境的溫度。本次計算中取對流換熱系數(shù)h為35W?m-2?K-1，環(huán)境溫度取為20oC，如圖4.6所示。

圖4.6 設置對流換熱邊界

除了對流散熱，物體還將以輻射的方式向周圍介質釋放熱量，由于在焊接過程中以熱輻射的形式釋放的能量較少，因此在本文在計算中沒有考慮熱輻射的影響。在實際計算中，為了提高計算效率，也可以通過調整對流換熱系數(shù)來考慮熱輻射的影響。

在Load模塊中創(chuàng)建一個面熱流（Surface heat flux），并將面熱流作用區(qū)域選擇為熱源作用的區(qū)域，如圖4.7所示。

圖4.7 選擇面熱流作用區(qū)域

需要注意面熱流作用區(qū)域只需要覆蓋熱源作用區(qū)域即可，盡管選擇平板整個表面不會對計算結果產生太大的影響，但由于在所有面熱流作用區(qū)域的單元積分點上ABAQUS都會調用DFLUX計算熱流密度，這顯然增加了不必要的計算量。

將面熱流密度的分布類型（Distribution）改為User-defined，并將幅值（Magnitude）設為1.

通過Predefined Field定義初始溫度場，并將工件初始溫度設定為20oC。

在Mesh模塊中，定義整個平板網格的單元類型為8節(jié)點線性熱傳導單元DC3D8。

在Job模塊創(chuàng)建熱傳導分析任務plate_thermal，并指定分析所使用的子程序文件，如圖4.8所示。

圖4.8 指定子程序文件

圖4.9給出了焊接14s時平板表面的溫度分布，取IN718材料的熔點為1240oC，圖4.9中灰色區(qū)域即為焊接熔池。

圖4.9 平板溫度場分布（t=14s）

提取沿焊接方向的熱循環(huán)曲線，如圖4.10所示。

圖4.10 焊接熱循環(huán)曲線

圖4.11給出了焊接熔池沿板厚方向的形貌。

圖4.11 焊接熔池形貌

從圖4.11中可以看出，由于焊接速度過快，因此平板在焊接過程中并沒有焊透，并且熔池呈現(xiàn)出扁平的形狀，而并非球形。

5. 平板焊接應力場計算

在彈塑性應力分析中，總的應變可以表示為：

其中εije為彈性應變，εijp為塑性應變，εijth為熱應變。

熱應變εth可表示為：

其中α為材料熱膨脹系數(shù)，Tref為參考溫度。

在順序熱力耦合分析中，通過熱傳導分析計算得到的溫度分布將會施加到應力模型上，從而產生相應的熱應變，并最終計算得到對應的焊接應力場。

應力分析使用的模型可以直接通過修改熱傳導分析模型獲得。在模型樹中復制熱傳導分析模型，并將其更改為plate_mechanic。在Property模塊添加應力分析需要使用的材料屬性，包括楊氏模量、泊松比、平均熱膨脹系數(shù)和屈服應力，如圖5.1所示。

圖5.1 定義材料屬性（應力分析）

計算時采用的彈塑性本構為等向強化模型，由于前面只給出了IN718的初始屈服應力，沒有給出相應的硬化曲線，因此這里假設為理想彈塑性材料。為了考慮材料熔化時的力學行為，在Plastic菜單的Suboptions中選擇Anneal Temperature，并將溫度值設為1240oC，即材料的熔點，如圖5.2所示。

圖5.2 設置材料熔點

當材料點的溫度超過Anneal Temperature時，ABAQUS會自動將材料點的等效塑性應變值（PEEQ）設為0，并假設材料為理想彈塑性，只有當溫度低于Anneal Temperature才會恢復硬化行為。相應的，在圖5.1中需要指定熔點對應的屈服應力值，這里假設與前一個溫度數(shù)據對應的屈服應力相同。

在Step模塊創(chuàng)建兩個分析步，分別代表焊接過程和冷卻過程，計算時間分別為21s和1000s，與熱傳導分析設定的時間保持一致，分析類型選擇Static, General。打開Nlgeom（大變形開關），設置初始增量步為0.01，最小增量步為1E-5，最大增量步分別為21和1000。

在Interaction模塊刪除熱傳導分析創(chuàng)建的對流換熱接觸。在Load模塊中刪除熱傳導分析定義的面熱流密度，約束對稱面上節(jié)點沿Y方向的位移自由度，選擇平板端部的節(jié)點，并約束X和Z方向的位移自由度，如圖5.3所示。

圖5.3 邊界條件設置

在Predefined Field中創(chuàng)建溫度場，選擇平板所有節(jié)點，將Distribution改為From resultsor output database file，即通過輸出數(shù)據定義溫度場。在File中選擇熱傳導分析計算得到的輸出數(shù)據庫文件（plate_thermal.odb），將Begin step設為1,Begin increment設為0，End step設為1，End increment留空，如圖5.4所示。

圖5.4 讀取焊接過程溫度歷程數(shù)據

圖5.4表示在焊接過程對應的分析步（分析步1）中讀取熱傳導分析中焊接過程的溫度歷程數(shù)據。同理，在冷卻過程的應力分析中選擇熱傳導分析在冷卻過程對應的分析步，如圖5.5所示。

圖5.5 讀取冷卻過程溫度歷程數(shù)據

在Mesh模塊中更改單元類型為帶有減縮積分的8節(jié)點線性六面體單元C3D8R。在Job模塊創(chuàng)建任務文件plate_mechanic，由于應力分析中并沒有調用DFLUX子程序，因此不再需要用到用戶子程序文件。

通過應力分析得到工件冷卻過后的Von Mises應力分布如圖5.6所示。

圖5.6 平板Von Mises應力分布（t=1021s）

從圖5.6中可以看出，在焊接過后平板出現(xiàn)了翹曲變形，焊接產生的最大殘余應力位于焊縫區(qū)域，約為440MPa。由于在計算中沒有考慮材料的硬化行為，因此計算得到的最大殘余應力與圖2.3中母材在常溫下的屈服應力非常接近。

圖5.7給出了焊接過程的等效塑性應變（PEEQ）分布。

圖5.7 焊接過程的等效塑性應變（PEEQ）分布（t=14s）

從圖中可以看出，在材料屬性中設置了Anneal Temperature后，熔池區(qū)域的PPEQ值為0，而熔池后方區(qū)域由于溫度逐漸下降至低于熔點，因此開始累積塑性應變值。

圖5.8給出了沿平板中心的殘余應力分布。

圖5.8 沿平板中心的殘余應力分布

從圖中可以看出，在焊縫區(qū)域，平板所受的縱向應力（S11）遠大于橫向應力（S22），并且當遠離焊縫之后，橫向應力出現(xiàn)了明顯下降的趨勢，這也表明焊接產生的殘余應力主要集中在焊縫區(qū)域，這主要是由于焊縫區(qū)域材料不均勻熱變形引起的。

圖5.9給出了沿焊接線上的殘余應力分布。

圖5.9 沿焊接線的應力分布

從圖5.9可以看出縱向應力和橫向應力在焊接線上表現(xiàn)出了相同的分布趨勢，只是應力水平有所不同。并且沿著焊接線的中部區(qū)域，出現(xiàn)了應力逐漸上升的趨勢。盡管平板是完全對稱的，但這里的應力分布卻并不是對稱的，本文認為這主要是由于在焊接過程中焊件的整體溫度逐漸上升，導致最終獲得殘余應力也出現(xiàn)了逐漸上升的趨勢。

參考文獻

[1] Tanner, David W. J., Life assessment of welded INCONEL 718 at high temperature[D].University of Nottingham, 2009.

原作者

tian