多物理场耦合仿真：热-结构相互作用的分区耦合策略

引言

实际工程结构很少只受单一物理场支配：涡轮叶片在高温燃气与离心载荷联合作用下工作，电子封装在功率循环中产生热应力，制动盘在摩擦生热与机械变形的反馈中可能出现热弹性失稳。这类热-结构耦合问题的数值求解有两条路线：整体式（monolithic）把所有场变量装进同一个非线性方程组联立求解；分区式（partitioned）则复用成熟的单场求解器，通过交界面数据交换与交错迭代实现耦合。本文推导线性热弹性耦合方程，分析交错迭代的收敛性，并给出一个弹簧约束钢杆的完整 Python 耦合仿真实现。

理论基础

线性热弹性方程组

固体域内，温度场 $T(\mathbf{x}, t)$ 满足含热弹性耦合源项的热传导方程，位移场 $\mathbf{u}$ 满足含热应变项的动量平衡方程：

$$\rho c_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) - \beta T_0 , \frac{\partial}{\partial t}(\nabla \cdot \mathbf{u})$$

$$\rho \ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma}, \qquad \boldsymbol{\sigma} = \mathbb{C} : \left(\boldsymbol{\varepsilon} - \alpha_T (T - T_0),\mathbf{I}\right)$$

其中 $\alpha_T$ 为线膨胀系数，$\beta = (3\lambda + 2\mu)\alpha_T$ 为热弹性耦合系数，$\mathbb{C}$ 为弹性张量，$\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^{\mathrm{T}})$。右端耦合源项的物理意义是：体积膨胀吸热、压缩放热——金属材料在弹性循环加载下的温度波动（约 $10^{-3}$ K/MPa 量级）正是红外热像应力测量的物理基础。

分区耦合与交错迭代

分区策略在每个时间步内交替求解两个场：先以冻结的结构量求解温度场，再据新温度求解结构响应，如此反复直至收敛。按数据交换方式可分为：

弱耦合（显式交错）：每步只交换一次数据，不迭代。代价是引入 $O(\Delta t)$ 的耦合分裂误差，仅适合耦合弱、时间步小的场合；
强耦合（隐式交错）：步内迭代至残差收敛，通常用 Gauss-Seidel 方式串行传递两场最新值，辅以 Aitken $\Delta^2$ 松弛抑制振荡、加速收敛。

固定点迭代的收敛性取决于耦合算子的谱半径。定义第 $k$ 次迭代的界面残差 $r_k = T_k – T_{k-1}$，Aitken 松弛因子按下式自适应更新：

$$\omega_k = -\omega_{k-1}, \frac{r_{k-1}^{\mathrm{T}} (r_k - r_{k-1})}{|r_k - r_{k-1}|^2}$$

当物理耦合强（如大膨胀系数、大温差、瞬态载荷）时，不动点迭代可能发散，此时需减小松弛因子或改用界面拟牛顿法（IQN-ILS）。耦合策略对比如下：

策略	步内迭代	精度	稳定性	典型应用
弱耦合显式交错	无	$O(\Delta t)$	条件稳定	弱耦合、长时演化
Gauss-Seidel 强耦合	有	高	中等	一般热-力耦合
Aitken/IQN 加速强耦合	有	高	好	强耦合、流固耦合
整体式 Newton	—	高	好	小自由度、强非线性

数值方法与实现

以下代码模拟一端固定加热、另一端由弹簧约束的一维钢杆：温度场用全隐式差分（Thomas 算法）求解，结构场利用平衡条件凝聚为均匀应力的解析表达式，两场在每个时间步内做 Gauss-Seidel 交错迭代直至收敛，纯 NumPy 即可运行：

import numpy as np

def thermoelastic_rod(N=80, L=0.1, t_end=60.0, dt=0.05):
    """弹簧约束钢杆的热-结构强耦合瞬态分析（分区交错迭代）"""
    E, aT, kT, rho_c = 210e9, 1.2e-5, 45.0, 3.5e6   # 钢：模量/膨胀系数/导热/体积热容
    beta = E * aT                                     # 热弹性耦合系数
    T0, Tl, Tinf, h = 300.0, 520.0, 300.0, 200.0      # 初温/左端温度/环境/换热系数
    ks, A = 5.0e8, 1.0e-4                             # 端部弹簧刚度 N/m、截面积 m^2
    dx = L / N
    Fo = (kT / rho_c) * dt / dx**2                    # 网格Fourier数（隐式无限制）

    # 隐式格式三对角系数：内部节点 + 左端Dirichlet + 右端对流
    Bi = h * dx / kT
    a = -Fo * np.ones(N); b = (1 + 2 * Fo) * np.ones(N); c = -Fo * np.ones(N)
    a[0] = 0.0; b[0] = 1.0; c[0] = 0.0                # 左端 T = Tl
    a[-1] = -1.0; b[-1] = 1.0 + Bi; c[-1] = 0.0       # 右端 -k dT/dx = h(T - Tinf)

    def thomas(a, b, c, d):
        n = len(d); cp = np.zeros(n); dp = np.zeros(n)
        cp[0] = c[0] / b[0]; dp[0] = d[0] / b[0]
        for i in range(1, n):
            m = b[i] - a[i] * cp[i - 1]
            cp[i] = c[i] / m
            dp[i] = (d[i] - a[i] * dp[i - 1]) / m
        x = np.zeros(n); x[-1] = dp[-1]
        for i in range(n - 2, -1, -1):
            x[i] = dp[i] - cp[i] * x[i + 1]
        return x

    T = np.full(N, T0)
    eps_prev = np.zeros(N)
    sigma = 0.0
    it_hist = []
    for _ in range(int(t_end / dt)):
        T_old = T.copy()
        T_it, sg_it = T.copy(), sigma
        for it in range(30):                          # 步内交错迭代
            dT = T_it - T0
            # 结构凝聚：弹簧约束 -> 杆内均匀应力
            sg_new = -ks * np.sum(aT * dT) * dx / (A + ks * L / E)
            eps = sg_new / E + aT * dT                # 总应变 = 机械应变 + 热应变
            eps_dot = (eps - eps_prev) / dt
            # 热场：隐式扩散 + 热弹性耦合源 -βT0 ε̇
            dvec = T_old - beta * T0 / rho_c * eps_dot * dt
            dvec[0] = Tl; dvec[-1] = Bi * Tinf
            T_new = thomas(a, b, c, dvec)
            res = max(np.max(np.abs(T_new - T_it)), abs(sg_new - sg_it) / E)
            T_it, sg_it = T_new, sg_new
            if res < 1e-9:
                break
        it_hist.append(it + 1)
        T, sigma = T_it, sg_it
        eps_prev = sigma / E + aT * (T - T0)
    return T, sigma, it_hist

T, sigma, its = thermoelastic_rod()
print(f"杆内最高温度 {T.max():.1f} K, 平均耦合迭代 {np.mean(its):.1f} 次/步")
print(f"热应力 sigma = {sigma/1e6:.1f} MPa (受压)")

物理结果符合直觉：杆受热膨胀却受弹簧约束，内部产生数十至百余 MPa 量级的压应力；由于耦合源项很小（$-\beta T_0 \dot{\varepsilon}$ 对应的温度反馈约 $10^{-3}$ K/MPa），交错迭代通常 3–5 次即可收敛。若把端部弹簧换成刚性壁（$k_s \to \infty$），约束趋强、耦合变硬，所需迭代次数明显增加，甚至需要欠松弛——这正是强耦合问题的典型表现。

工程案例分析

问题：功率模块焊层热疲劳评估

某 IGBT 功率模块在 40 s 通断循环下工作，芯片功耗 200 W，焊层（SnAg3.5）厚 120 μm，直接键合铜基板通过水冷维持在 65 °C。评估焊层热应力幅以预测疲劳寿命。

组件	厚度 (mm)	导热系数 (W/m·K)	弹性模量 (GPa)	线膨胀系数 (10⁻⁶/K)
Si 芯片	0.35	148	130	2.6
焊层	0.12	57	35	24.0
DBC 铜层	0.30	385	120	16.5

芯片（$\alpha_T = 2.6$）与焊层（$\alpha_T = 24.0$）的膨胀失配接近一个数量级，是应力集中的根源。分区耦合仿真（热场 0.1 s 步长、结构场凝聚模型）给出焊层剪应力幅约 18 MPa；按 Coffin-Manson 型经验关系换算，特征寿命约 $2\times10^4$ 循环，与功率循环试验报告值处于同一量级。若仿真忽略温度对焊层模量的软化作用（$E$ 从 25 °C 到 125 °C 下降约 30%），应力幅将被高估约 25%——这说明即使”弱”的物性反馈也必须纳入耦合迭代。

结论与展望

本文给出了热-结构耦合的方程体系、强弱耦合策略与一个可运行的强耦合交错迭代实现。要点是：优先复用单场求解器、以界面残差控制步内迭代、用 Aitken 或拟牛顿法处理强耦合。进一步的工程议题包括：非共形网格间的保守插值（径向基函数、加权残差映射）、接触热阻与摩擦生热的界面模型、以及面向寿命预测的热-力-损伤多字段耦合。对制动热弹性失稳这类正反馈问题，交错格式的稳定性分析本身即是研究课题，需谨慎选择时间步长与松弛策略。