3. 热机电耦合控制方程与弱耦合/强耦合策略
TSV与RDL的热机电可靠性仿真，本质上不是“把温度算出来再加到结构上”这么简单。温度场会改变材料刚度与屈服强度，应力又会通过塑性耗散和接触热阻反过来影响温度；电流在铜柱中产生焦耳热，同时原子在电子风力与应力梯度共同作用下发生电迁移，而电迁移造成的空洞与凸起又会改变局部电流密度与应力分布。要正确预测这类问题，必须先把各物理场的控制方程写清楚，再决定它们之间以什么方式、按什么顺序交换耦合变量。
一、五个物理场的控制方程
热传导方程描述温度场演化。瞬态形式为 ρcp∂T/∂t = ∇·(k∇T) + qgen，其中 qgen 包含焦耳热 J·E 与塑性/蠕变耗散。稳态分析时左端项为零。TSV铜柱的 k 约为 380～400 W/(m·K)，而 SiO2 介质层仅约 1.4 W/(m·K)，这种三四个数量级的差异决定了热量主要沿铜柱轴向导出，横向散热很弱。
热弹性方程由平衡方程 ∇·σ + f = 0 与本构关系 σ = C:(ε − εth − εp − εc) 组成，热应变 εth = α(T − Tref)I。这里 Tref 是零应力参考温度，必须与工艺历史一致：若铜在 250 ℃ 退火后冷却到室温，取 Tref = 250 ℃ 比取 25 ℃ 更接近真实残余应力状态。
塑性方程采用 von Mises 屈服准则与等向/随动硬化。屈服函数 f = σeq − σy(εpeq, T) ≤ 0，塑性流动方向由关联流动法则给出 dεp = dλ·∂f/∂σ。铜的屈服强度对温度极敏感，25 ℃ 时约 200～300 MPa，200 ℃ 时可能降到 80～120 MPa，因此温度场必须先于塑性更新。
蠕变方程常用 Norton–Bailey 幂律：dεc/dt = A·σn·exp(−Q/(RT))。对电镀铜，典型参数为 n ≈ 4～8，Q ≈ 0.8～1.2 eV。蠕变是焊点与铜柱在高温保载阶段的主导变形机制，其速率同时依赖应力与温度，是典型的强耦合项。
电迁移方程描述原子通量：Ja = (C·D/(kT))·(Z*·e·ρ·j − Ω·∇σh)。第一项是电子风力，第二项是应力梯度驱动的背应力。当 j 超过临界值（铜中约 104～105 A/cm²）且温度较高时，通量散度 ∇·Ja 不为零，在阴极侧形成空洞、阳极侧形成凸起。
二、弱耦合与强耦合的判定
弱耦合（顺序耦合）指各场依次求解，把上一场的解作为下一场的载荷输入，不迭代或只做少量外层迭代。典型流程是：先解稳态热传导得到 T 分布，再把 T 作为体载荷施加到力学模型，求解热应力与塑性/蠕变。这种策略适用于温度场受力学变形影响很小的场景，例如芯片在恒定环境温度下的热循环分析。它的优点是每个子问题规模小、可分别用最合适的求解器，计算成本低。
强耦合（全耦合）指把所有场的自由度组装进同一个方程组同时求解，或在每个时间步内反复迭代直到各场残差同时收敛。当耦合项量级不可忽略时必须用强耦合，例如：焦耳热使温度升高，温度升高又改变电阻率 ρ(T) 从而改变电流密度分布，电流密度变化再改变焦耳热——这个正反馈回路若用弱耦合会低估热点温度。电迁移与应力的双向作用同样属于强耦合：应力梯度影响原子通量，原子重分布又改变应力场。
三、求解顺序与变量传递
实际工程中更常用的是“分区强耦合”：把热场与电场合并为一个子问题（因为它们通过焦耳热和电阻率直接耦合），力学场与电迁移合并为另一个子问题，两个子问题之间在每个时间步内做 2～4 次交替迭代。这样既控制了方程组规模，又保留了主要反馈路径。
变量传递时要注意三个细节。第一，网格不一致时的插值：热场网格通常较粗，力学场网格较细，把 T 从热网格映射到力学网格必须用守恒插值，避免在界面处产生虚假温度梯度。第二，时间步匹配：蠕变和电迁移是慢过程，特征时间可达数百小时，而热瞬态可能只有几秒，需要采用子循环或自适应时间步。第三，参考状态一致性：热应变、塑性应变、蠕变应变都相对于同一个 Tref 和零应力状态定义，切换求解顺序时不能重复计入或漏计某一分量。
以一个典型 TSV 热循环为例：温度从 25 ℃ 升到 125 ℃ 再降回 25 ℃，升降温速率 10 ℃/min，高温保载 30 min。若采用弱耦合，先算瞬态温度场，再把每个时刻的 T 映射到力学模型计算弹塑性；若采用分区强耦合，则在每个时间步内先更新 T 和 j，再迭代求解 σ 与 εc，直到蠕变应变增量收敛。后者在保载阶段的应力松弛预测上明显更接近实际，因为蠕变速率对当前应力而非上一步应力敏感。
总结来说，控制方程给出“每个场怎么变”，耦合策略决定“它们怎么互相影响、按什么顺序算”。对 TSV 与 RDL，推荐默认采用热-电分区强耦合、力-电迁移分区强耦合的两级迭代框架；只有在温升小、电流密度低、蠕变不显著的快速筛选场景下，才退化为纯顺序弱耦合以节省计算资源。