——科研到工程：Abaqus Goldak 雙橢球 + FROM FILE 實現可復現實驗結果（含 Goldak 熱源 DFLUX ）

適用人群：做焊接/鍵合殘余應力/變形預測、增材制造熱-力場分析的工程師與研究生

代碼環境：Abaqus/CAE 2019（Python 2.7），Abaqus/Standard（DFLUX Fortran 子程序）

本文提供 兩個腳本（Abaqus/CAE Python 自動網格建模腳本 + Fortran DFLUX 熱源子程序）梳理成一套可復用的 有限元計算流程：

• 從物理到實現 的清晰鏈路：能量輸入 → 傳熱 → 溫度–時間歷程 → FROM FILE 映射 → 彈塑性力學響應；
• 建模與求解流程：幾何、分區、網格、邊界、步長、輸出與文件命名；
• Goldak 雙橢球熱源與熱力耦合理論： 在 DFLUX 中的實現原理與關鍵參數；

目錄

- 用 Abaqus 做焊接/鍵合熱-力耦合的“一鍵批量建模與計算”
- 目錄
    - 1. 為什么要做焊接/鍵合熱-力耦合？
    - 2. Goldak 雙橢球熱源與能量守恒
    - 3. 總體流程與工程目錄
    - 4. 熱分析建模要點（Thermal）
    - 5. 力學分析建模要點（Mechanical）
    - 6. 自動化批量建模腳本（Python，最終版）
    - 7. DFLUX：Goldak 體熱源子程序（Fortran）
    - 8. 模型驗證
    - 9. 參考參數與推薦文獻

1. 為什么要做焊接/鍵合熱-力耦合？

焊接/鍵合是強非線性、強非穩態的多物理場過程：移動熱源瞬時把能量輸入到極小體積，熱擴散與對流/輻射把能量帶走，材料在不同溫度區間內經歷彈性–塑性–循環硬化乃至回復。順序耦合（先熱后力）是工業中最穩健的路線：

• 熱分析（Heat Transfer）：移動熱源 + 換熱邊界 → 得到溫度–時間歷程 T(x,y,z,t)；
• 映射（FROM FILE）：把熱場隨時間讀入力學模型；
• 力學分析（Static, nlgeom）：考慮 E(T)、σy(T)、α(T) 等 → 輸出殘余應力/變形。

為什么不能只做熱或只做力？

• 只做熱：沒有熱–力耦合的應力演化，無法預測殘余場；
• 只做力：沒有真實的溫度歷程驅動，熱應變與材料退化無從談起。

工程意義：

• 快速評估工藝窗口（功率/焊速/熱源形參）對峰溫、HAZ、殘余應力與翹曲的影響；
• 用自動化腳本把“手工建模”變成“可復用流程資產”，支撐 DOE/靈敏度/優化。

要讓結果可信，關鍵是：能量守恒、邊界換熱量級合理、材料熱物性/力學參數隨溫度變化；在力學側需最小約束消除剛體模態，并與熱網格一致以確保映射穩定。

2. Goldak 雙橢球熱源、能量守恒與熱力耦合

符號：坐標 ；熱源中心位置 ；半軸 ；有效功率 ；分配系數 （滿足 ）。

前半橢球（front，）

后半橢球（rear，$x

分段表達

能量守恒

軌跡（恒速 ，起點 ，起始時刻 ）

熱傳導控制方程（瞬態）

在域 、時間區間 內，溫度場 滿足瞬態能量守恒（不考慮相變）：

其中 為密度， 為定壓比熱， 為導熱系數， 為體熱源（W/m^3）。

初始條件：

邊界條件（三類任選/組合）：

1. 指定溫度（Dirichlet）：
2. 指定熱流（Neumann）：
3. 對流 + 輻射（Robin）：

其中 為對流換熱系數， 為表面發射率， 為 Stefan–Boltzmann 常數， 為環境溫度。

力學控制方程與熱應變

焊接/鍵合后的固體力學響應（小–中變形，幾何非線性可按需開啟）：

靜力平衡：

其中 為 Cauchy 應力， 為體力密度（可忽略）。

位移邊界：

彈性本構（各向同性）：

含溫度依賴的 。

塑性與硬化（示例：J2 + 組合硬化；僅示意）：

• 等效應力：， 為偏應力；
• 屈服函數：；
• 演化：，并可含各向/隨動硬化項；
• 一致性條件：。

順序耦合的本質：先用熱方程得到 ，再把它作為已知外場驅動固體力學問題（通過 與溫度依賴材料參數），時間上保持同一時間軸或可匹配的時間段。

熱-力耦合

熱傳導與移動熱源

其中 即 Goldak 雙橢球體熱源; 表面邊界含對流/輻射條件:

熱彈塑性平衡方程

耦合流程

順序耦合中，熱分析得到的 （或其在積分點/節點的離散值）通過 TEMPERATURE, FILE=... 輸入到力學模型。

1. 解熱方程得 ;
2. 在力學步中按時間步讀入 , 由 生成熱應變, 并以溫度退化的 與屈服準則/硬化規律推進塑性;
3. 在冷卻階段, 不可恢復的塑性應變與結構約束共同“鎖定”殘余應力。

要求：

1. 網格一致（最好共享相同拓撲與節點）；
2. 時間軸覆蓋（力學步的時間點應落在熱步范圍內，或可插值）；
3. 參考溫度一致（材料模型中的 與初始溫度設置一致）；
4. 邊界與最小約束合理，去除剛體模態。

3. 總體流程與工程目錄

• Python 腳本自動創建幾何/網格/材料/邊界與分析步，分別寫出熱（THinp）與力學（MEinp）輸入；
• DFLUX 計算體熱源 q(x,y,z,t)；
• Mechanical 通過 *TEMPERATURE, FILE=<thermal job> 讀取熱場。

工程目錄（由腳本在“當前工作目錄 CWD”下自動生成）：

<當前工作目錄 CWD>/
 ├─ THinp/   # 熱分析 .inp（文件名以 T- 開頭）
 ├─ MEinp/   # 力學分析 .inp（文件名以 M- 開頭）
 └─ UFLUX.for（或 DFLUX.for）

4. 熱分析建模要點（Thermal）

• 單元：一階 DC3D8；網格沿焊縫中心、端部與厚度方向加密；
• 邊界：除對稱面外的外表面施加對流+輻射；
• 步長：自動增量（初始 0.02 s，最小 ）；
• 初溫：20 ℃；
• 熱源：Body heat flux, USER（由 DFLUX 提供）。

5. 力學分析建模要點（Mechanical）

• 單元：一階 C3D8R；與熱網格一致；
• 約束：Y=0 對稱；兩端兩點最小約束去除剛體模態（右端約束 u3，左端約束 u1,u3）；
• 材料：E(T)，σy(T)，組合硬化+循環硬化，α(T)、cp(T)、k(T)；
• 步長：Static, nlgeom=ON，時間長度與熱分析一致；
• 溫度：*TEMPERATURE, FILE=<熱分析 job 名>，分 3 個階段（Welding/Cool1/Cool2）。

6. 自動化批量建模腳本（Python）

該腳本以當前工作目錄 CWD 為根，在 THinp/ 和 MEinp/ 寫入成對輸入文件；支持單值或區間 min-max-n 的參數掃參。

使用

• CAE GUI：File → Run Script... 選擇腳本；
• 命令行：abaqus cae noGUI=script.py -- 0.16 0.06 0.02 0.08 0.00227 1100
• 或區間：-- 0.16 0.06 0.02 0.06-0.10-3 0.0018-0.003-5 900-1300-3。

腳本（節選，僅突出關鍵處）

# -*- coding: mbcs -*-
from abaqus import *
from abaqusConstants import *
import section, regionToolset, part, material, assembly, step, load, mesh, job, sketch
import os, sys, numpy as np
from collections import OrderedDict

# 以“當前工作目錄”為根
BASE_DIR  = os.getcwd()
THERM_DIR = os.path.join(BASE_DIR, 'THinp')
MECH_DIR  = os.path.join(BASE_DIR, 'MEinp')

def ensure_output_directories():
    for d in (THERM_DIR, MECH_DIR):
        if not os.path.isdir(d):
            os.makedirs(d)

# ……（此處省略：材料、幾何、分區、網格、步驟、邊界、輸出等函數）……

# 主流程：循環參數 → 寫熱分析 inp → 寫力學 inp
if __name__ == '__main__':
    ensure_output_directories()
    # 解析參數（支持單值或  min-max-n 形式）
    # grid = {'length': [...], 'width': [...], ...}
    # ……

    for L in grid['length']:
        for W in grid['width']:
            for H in grid['height']:
                for BL in grid['bead_length']:
                    for VS in grid['arc_speed']:
                        for PWR in grid['heat_input']:
                            os.chdir(THERM_DIR)
                            th_job = make_thermal_model(L, W, H, BL, VS, PWR)
                            os.chdir(MECH_DIR)
                            me_job = make_mechanical_model(L, W, H, BL, VS, PWR, thermal_job_name=th_job)
                            os.chdir(BASE_DIR)

    print('All input files generated successfully.')

7. DFLUX：Goldak 體熱源子程序（Fortran）

      SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,JLTYP,
     &                 TEMP,PRESS,SNAME)
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION FLUX(1),TIME(2),COORDS(3)
      CHARACTER*80 SNAME

C ---- 用戶參數（需與 Python/工藝一致） ----
      Q      = 1100.0D0
      A1     = 4.5D-3
      A2     = 6.5D-3
      B      = 6.8D-3
      C      = 2.6D-3
      F1     = 1.0D0
      F2     = 1.0D0
      L      = 0.080D0
      V      = 2.27D-3       ! m/s
      T0     = 0.0D0
      X0     = -L/2.0D0

      PI     = 3.141592653589793D0
      JLTYP  = 1

C ---- 當前時間與位置 ----
      T  = TIME(1)
      X  = COORDS(1)
      Y  = COORDS(2)
      Z  = COORDS(3)

C ---- 熱源中心并限位 ----
      XNOW = X0 + V*(T - T0)
      IF ( XNOW .GT.  0.5D0*L ) XNOW =  0.5D0*L
      IF ( XNOW .LT. -0.5D0*L ) XNOW = -0.5D0*L

C ---- 規范化系數 ----
      
            .......................

      IF ( X .GE. XNOW ) THEN
         FLUX(1) = QF
      ELSE
         FLUX(1) = QR
      ENDIF
      RETURN
      END

若你用“線能量” （J/mm）與焊速 （m/s），則 Q = η * E_\ell * v 。

8. 模型驗證

1. 準備輸入

• Python 腳本：幾何尺寸 [L,W,H]、焊縫長度 BL、焊速 VS、有效功率 Q。
• Fortran 子程序：Goldak 形參 a1,a2,b,c 與分配 f1,f2，以及功率 Q、速度 v、長度 L、起點 X0 等。

1. 熱分析輸入（T-*.inp）

• 幾何與分區：長方體試樣；以 Y=0 為對稱平面；在 焊縫中心線與兩端 做加密分區（便于捕捉尖銳梯度與端部效應）。
• 網格：一階熱單元 DC3D8；全局較粗、焊縫區加密；
• 材料（熱物性）：Density、SpecificHeat(T)、Conductivity(T)；
• 邊界換熱：除 Y=0 對稱面外，其他外表面施加 Film + RadiationToAmbient；
• 分析步：Welding（持續 BL/VS），Cool1，Cool2，均為自動增量控制；
• 體熱源：*DFLUX（或輸入文件設為 USER），在 Welding 步生效；
• 輸出：NT/HFL/RFL 等場量。

1. 求解熱場

• 提交 T-*.inp + DFLUX.for（或 user=UFLUX.for）作業，得到 T-*.odb。
• 檢查代表性點溫度–時間曲線、峰值溫度、等溫線是否合理。

1. 力學分析輸入（M-*.inp）

• 幾何/網格：與熱側一致，一階 C3D8R；
• 約束：Y=0 對稱；兩端兩點 最小約束（去剛體模態：右端約束 u3，左端約束 u1,u3）；
• 材料（力學）：Elastic(T)、Plastic(Combined) + CyclicHardening(T)、Expansion(T)，并保留 cp/k(T) 以便溫度場讀取時內部一致；
• 分析步：Welding、Cool1、Cool2，與熱側 時間軸一致；
• 溫度映射：

*TEMPERATURE, FILE=<T-...>, OP=NEW
** Step-1: Welding   beginStep=1
** Step-2: Cool1     beginStep=2
** Step-3: Cool2     beginStep=3

• 輸出：S/U/LE/PEEQ/NT。

1. 求解力學場

• 提交 M-*.inp，Abaqus 將從 T-*.odb 讀取每一步對應的溫度場；
• 檢查殘余應力分布（縱向/橫向/厚向）、等效塑性應變、焊后翹曲。

1. 后處理與掃參

• 批量參數（幾何/工藝）→ 自動生成多組 T-*/M-*.inp → 批處理提交 → 統一提取峰值溫度、熔寬/熔深近似、最大殘余應力、變形等 → 建響應面或靈敏度分析。

FAQ（高頻問題）

• 能只用 FILM/RADIATE 做熱源嗎？ 不能，那是邊界換熱；DFLUX 是體熱源。
• 為什么不做全耦合？ 工程尺度下成本高、材料相依更復雜；順序耦合更穩健。
• FROM FILE 總失敗？ 先查網格一致與 Job 名，時間軸覆蓋，再看插值容差設置。
• Goldak 與高斯面熱源區別？ Goldak 是體熱源且前后不對稱；面熱源多用于薄板近表面近似

9. 參考參數與推薦文獻

• Goldak J., Chakravarti A., Bibby M., A new finite element model for welding heat sources, Metallurgical Transactions B, 1984.
• Lindgren L.E., Modelling for residual stresses and distortions, 2001–2014（系列綜述）。
• ABAQUS Analysis User’s Guide：Heat Transfer、DFLUX、TEMPERATURE FROM FILE.