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

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

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

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

2. 焊接材料

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

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

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

圖2.3 IN718楊氏模量及屈服應(yīng)力隨溫度變化

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

3. 焊接熱源模型

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

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

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

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

4. 平板焊接溫度場計(jì)算

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

圖4.1 定義材料屬性

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

圖4.2 平板有限元模型

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

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

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

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

圖4.4 熱傳導(dǎo)分析步的增量步設(shè)置

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

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

圖4.5 選擇對流換熱面

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

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

圖4.6 設(shè)置對流換熱邊界

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

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

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

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

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

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

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

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

圖4.8 指定子程序文件

圖4.9給出了焊接14s時平板表面的溫度分布，取IN718材料的熔點(diǎn)為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. 平板焊接應(yīng)力場計(jì)算

在彈塑性應(yīng)力分析中，總的應(yīng)變可以表示為：

其中ε_ij^e為彈性應(yīng)變，ε_ij^p為塑性應(yīng)變，ε_ij^th為熱應(yīng)變。

熱應(yīng)變ε^th可表示為：

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

在順序熱力耦合分析中，通過熱傳導(dǎo)分析計(jì)算得到的溫度分布將會施加到應(yīng)力模型上，從而產(chǎn)生相應(yīng)的熱應(yīng)變，并最終計(jì)算得到對應(yīng)的焊接應(yīng)力場。

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

圖5.1 定義材料屬性（應(yīng)力分析）

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

圖5.2 設(shè)置材料熔點(diǎn)

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

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

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

圖5.3 邊界條件設(shè)置

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

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

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

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

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

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

圖5.6 平板Von Mises應(yīng)力分布（t=1021s）

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

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

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

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

圖5.8給出了沿平板中心的殘余應(yīng)力分布。

圖5.8 沿平板中心的殘余應(yīng)力分布

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

圖5.9給出了沿焊接線上的殘余應(yīng)力分布。

圖5.9 沿焊接線的應(yīng)力分布

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

參考文獻(xiàn)

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

原作者

tian