基于分形映射与等效参数优化方法的 2.5D 封装高带宽存储器降阶热机械建模

《International Journal of Heat and Mass Transfer》：Reduced-order thermo-mechanical modeling of 2.5D package high-bandwidth memory based on fractal mapping and equivalent parameter optimization method

　　

•一种多级降阶建模框架，融合了分形理论与最小二乘参数反演方法。•复杂的重布线层（RDL）几何构型以自相似分形图案表示，保留关键力学特征。•弹性模量和热膨胀系数通过有限元仿真和最小二乘优化获取。•实现了高保真热-力学响应，均方根位移误差<0.02，热误差为10??。•该方法将仿真时间和内存使用量分别降低90%和80%，实现了高效的高性能仿真。

**1. 引言**

在以NVIDIA H200为代表的高性能计算单元（HPC/AI加速器）中，计算核心和高带宽存储器（HBM）通过先进封装技术和中介层结构实现高度集成。因此，单位面积的热流密度远高于传统通用计算平台[[1],[2],[3],[4],[5],[6]]。在AI训练和推理等长时间满载运行过程中，H200芯片核心区域产生的极高热流不仅对冷却系统构成严峻挑战，还通过中介层、封装基板和其他印刷电路板（PCB）结构将热量传导至周围的HBM单元。该过程产生显著的非均匀温升，从而诱发复杂的热膨胀驱动变形和力学响应[[7],[8],[9],[10],[11]]。从失效机制的角度来看，高性能计算单元产生的超高热流可在封装结构内产生显著的温度梯度和热应力集中，特别是在热膨胀系数（CTE）严重不匹配的区域，如硅-重布线层界面、硅通孔（TSV）-硅界面以及焊点-基板界面。长期的热-力学耦合可能导致多种可靠性问题，包括微凸点焊料接头的低周疲劳、重布线层的累积塑性变形以及TSV周围硅基体中的裂纹萌生[[12],[13],[14],[15]]。尽管这些损伤可能不会立即导致灾难性失效，但会逐渐降低信号完整性、增加接触电阻并提高存储单元的误码率。最终，这些效应表现为器件性能的渐进性退化和系统级故障频率的显著增加[[16],[17],[18]]。因此，在真实工况下准确评估PCB上HBM单元及其中介层和封装基板所承受的热应力，对于器件可靠性和长期稳定运行具有至关重要的意义。

然而，HBM中介层内存在的跨尺度微结构使得直接仿真和评估热-力学变形极为困难。在结构特征尺度方面，重布线层的线宽和间距、TSV的直径以及微凸点焊料的尺寸通常在10 μm量级，而中介层的整体横向尺寸和厚度则为毫米量级[[19],[20],[21],[22],[23],[24]]。如果采用完全解析的跨尺度建模方法并以1:1几何比例表示——显式建模每一条线路、TSV和焊点进行直接热-力学耦合仿真——则自由度数量将急剧增加。因此，计算资源消耗将呈指数级增长，仿真时间将变得不可接受，甚至可能无法实现稳定收敛[[25],[26],[27],[28],[29],[30]]。因此，这类高保真原始模型在芯片设计迭代和参数优化阶段几乎不切实际[31]。在工业研究与开发的实际工作流程中，封装和系统设计通常要求在有限的时间框架内进行多轮快速的结构和材料参数评估与优化。这提出了一个实际需求，即仿真模型需同时兼顾计算效率和物理精度。

仅依赖完全细化的模型不仅无法满足工程时间线要求，还可能由于数值噪声和局部几何细节的主导作用而掩盖结构层面的固有热-力学行为。因此，引入基于代表性体积元方法的降阶建模成为一种自然解决方案。该方法将复杂异质微结构均匀化为具有各向异性等效力学和热学性质的等效层。该策略在保持整体热膨胀变形和弯曲响应精度的同时，大幅减少了模型中的自由度数量[[32],[33],[34],[35],[36]]。近期研究为电子封装的跨尺度热-力学建模开发了高效的均匀化和数据驱动策略。Yoo等人[37]将自动铜线路图案分类与各向异性粘弹性均匀化相结合，实现了快速PCB翘曲预测；Yang等人[38]采用深度学习加速从复杂中介层布线图案中提取等效力学参数；Liu等人[39]进一步开发了针对TSV异质结构等效弹性表征的多级均匀化框架。然而，直接均匀化需要合适的代表性单元，而数据驱动方法通常需要足够大且具有代表性的训练数据集。

相比之下，与传统的RVE均匀化方法（直接以等效均匀材料替代异质微结构）不同，本文提出的框架在全解析几何与RVE均匀化模型之间引入了一个基于分形的中间参考模型。希尔伯特曲线作为一种经典的能在二维区域内以一维路径遍历同时允许控制线宽和间距的分形曲线，特别适合近似PCB中的线路图案[[40],[41],[42]]。首先以特征线宽和铜覆盖率作为代表性几何描述符，将不规则RDL图案映射为自相似分形几何结构。随后通过有限元响应拟合和最小二乘优化识别其等效热-力学参数，同时采用壳理论高效保留跨尺度薄RDL和通孔结构的力学贡献。这种分层策略使最终的分形-壳模型在保持计算可处理性的同时保留了关键微结构特征，从而为评估RVE选择和进一步均匀化提供了定量参考。因此，本工作的主要方法论贡献在于将分形几何降阶、基于响应的等效参数识别和壳表示整合为一个分层框架，用于验证高度异质电子封装的基于RVE的降阶建模。

**2. 方法**

本文提出的降阶建模框架遵循分层四阶段流程：(1) 从原始RDL图案中提取代表性几何描述符并构建自相似分形表示；(2) 通过有限元响应拟合和最小二乘优化识别等效热-力学参数（杨氏模量E和CTE）；(3) 采用薄RDL和通孔结构的壳表示构建分形-壳原始等效模型；(4) 对最终降阶模型进行RVE均匀化并验证其与原始等效模型的一致性。整体技术路线如图1所示。

**图1.** 融合分形映射、等效参数优化、壳建模和RVE均匀化的分层降阶建模框架技术路线图。

**2.1 基于有限元响应拟合和最小二乘参数反演的等效参数提取方法**

为获取具有复杂几何结构的原始重布线层（RDL）的等效力学和热膨胀参数，本研究提出了一种基于有限元响应拟合和最小二乘参数反演的等效参数提取方法。该方法通过比较在相同加载条件下原始分形单元与均匀化等效单元的位移响应来确定等效材料参数，随后通过最小二乘优化反演识别等效参数。所提出的参数提取框架可被解释为偏微分方程约束的反问题，其中通过最小化分形微结构模型与均匀化模型位移响应之间的差异来识别等效材料参数。该方法的数学框架如下所述：

**(1) 参数化有限元响应模型**

原始分形单元的真实结构定义在域Ω中。均匀化等效单元的材料参数向量定义为：

(1) **p** = [E, α]?

其中E表示等效杨氏模量，α表示等效线性热膨胀系数（CTE）。

对于给定的参数向量**p**，通过求解结构力学控制方程（结合上述有限元方法和线弹性本构关系），获得均匀化模型的位移场*u*(*x*;**p**)。由原始分形结构参考仿真获得的位移场记为*u*ref(*x*)。本研究仅选择位于位移提取曲线Γ上的离散节点作为响应对比点。

**(2) 离散位移响应向量的构建**

假设位移提取曲线Γ包含N个采样点，坐标为*x?*（*i* = 1, 2, …, *N*）。原始结构的参考位移向量定义为：

(2) **d**ref = [u ref z(*x*?), u ref z(*x*?), …, u ref z(*x*N)]?

均匀化等效模型的位移向量为：

(3) **d**(**p**) = [u_z(*x*?;**p**), u_z(*x*?;**p**), …, u_z(*x*N;**p**)]?

其中仅考虑z方向的位移分量来表征弯曲或热膨胀变形。

**(3) 最小二乘目标函数的构建**

为确保均匀化模型的整体变形行为尽可能接近原始结构，定义残差向量为：

(4) **r**(**p**) = **d**(**p**) ? **d**ref

采用最小二乘准则，构建目标函数为：

(5) *J*(**p**) = (1/2) Σ???? [u_z(**x?**;**p**) ? u ref z(**x?**)]2

目标函数的向量形式为：

(6) *J*(**p**) = (1/2) ‖**r**(**p**)‖2?

该目标函数代表均匀化模型预测的位移曲线与原始模型之间的整体误差。因此，优化问题可表述为：

(7) **p**\* = arg min_**p** *J*(**p**)

其中**p**\*表示最优等效参数向量。

**(4) 有限元-优化耦合求解**

由于位移响应*d*(**p**)由有限元求解器获得，目标函数对参数呈隐式依赖关系。因此，该问题属于偏微分方程约束优化问题。在每次参数迭代中，首先需要求解有限元离散方程：

(8) **K**(**p**)**U** = **F** + **F**th(**p**)

(9) **K**(**p**) = ∫_Ω **B**?**C**(E)**B** dΩ

其中**K**为参数相关的结构刚度矩阵。

获得节点位移向量**U**后，沿提取曲线的节点位移值被收集以构建**d**(**p**)，然后可评估目标函数*J*(**p**)。

**(5) 梯度计算与参数迭代**

在参数优化过程中，采用基于梯度的优化算法更新参数。目标函数的梯度可表示为：

(10) ?*J*(**p**) = **J**?**r** **r**(**p**)

(11) **J**_r = ?**r**/?**p**

其中**J**_r表示残差向量关于参数的雅可比矩阵。

参数更新公式为：

(12) **p**??? = **p**? ? λ(**J**?_r **J**_r)?1**J**?_r **r**

其中λ为步长因子。迭代过程持续进行，直至均匀化模型与参考模型的位移响应之间的残差平方和（RSS）满足预设收敛准则：

(13) RSS(**p?**) < ε

其中ε为预设的RSS容差。本研究采用ε = 0.1 mm2作为收敛阈值。RSS的详细定义见公式(14)。

**(6) 等效参数的物理诠释**

通过上述最小二乘优化过程，当均匀化单元的位移响应曲线在整体趋势上与原始分形单元最为接近时，所识别的等效参数被认为能够在宏观尺度上再现原始结构的全局力学响应。对于三点弯曲仿真，优化后的等效杨氏模量（E）代表重布线层在宏观弯曲条件下的有效刚度。对于热膨胀仿真，优化后的等效热膨胀系数（α）表征复杂线路结构在温度加载下的整体热膨胀行为。由于最小二乘方法本质上是最小化整体位移响应误差，因此该方法可以获得具有明确物理意义的等效材料参数，同时在优化过程中保持良好的数值稳定性。

为定量评估位移响应的吻合程度，采用残差平方和（RSS）作为拟合指标：

(14) RSS = Σ???? (u_hom_z,? ? u ref z,?)2

其中N为位移提取线上的采样点数量，u_hom_z,?和u ref z,?分别表示均匀化模型和参考分形模型在第i个采样点的z方向位移。等效材料参数通过迭代更新以最小化RSS。当RSS降至0.1 mm2以下时，认为优化收敛。该准则一致应用于等效杨氏模量和热膨胀系数的识别。

**2.2 Hill-Mandel能量一致性条件**

在基于代表性体积元的均匀化建模中，从微结构获得的等效材料参数必须满足Hill-Mandel能量一致性条件，以正确描述宏观结构的力学行为。该条件是微结构均匀化理论的理论基础，要求宏观应力-应变功率等于RVE上微观应力-应变功率的体积平均值。只有满足该条件时，从RVE中提取的等效参数才能确保宏观模型与真实微结构之间的能量等价性。

在连续介质力学的框架下，Hill–Mandel条件可以表示为：

(15) $\langle\sigma_{ij}\dot{\varepsilon}_{ij}\rangle=\Sigma_{ij}\dot{E}_{ij}$

其中，$\sigma_{ij}$ 和 $\varepsilon_{ij}$ 分别表示微观应力张量和微观应变张量，而 $\Sigma_{ij}$ 和 $E_{ij}$ 分别表示宏观应力和宏观应变。记号 $\langle \cdot \rangle$ 表示对代表性体积单元（RVE）的体积平均。

体积平均算子定义为：

(16) $\langle f\rangle=\frac{1}{V}\int_V f(\mathbf{x})dV$

其中 $V$ 是代表性体积单元的体积。相应地，宏观应力和应变可定义为对应微观场的体积平均：

(17) $\Sigma_{ij}=\langle\sigma_{ij}\rangle=\frac{1}{V}\int_V \sigma_{ij}(\mathbf{x})dV$

(18) $E_{ij}=\langle\varepsilon_{ij}\rangle=\frac{1}{V}\int_V \varepsilon_{ij}(\mathbf{x})dV$

Hill–Mandel条件要求在任意加载过程中，RVE内的平均应变能密度必须等于等效宏观材料的应变能密度，即：

(19) $\frac{1}{V}\int_V \sigma_{ij}\varepsilon_{ij}dV=\Sigma_{ij}E_{ij}$

该关系确保从微观结构得到的等效弹性刚度矩阵 $\mathbf{C}_{eff}$ 能正确表征真实材料在宏观尺度上的整体刚度。在有限元实现中，必须对RVE施加适当的边界条件，以使其与宏观应变场保持一致，从而满足Hill–Mandel条件。本文采用周期性边界条件。该方法要求RVE相对边界之间的位移差满足：

(20) $u_{+i}-u_{-i}=E_{ij}L_j$

其中 $L_j$ 表示RVE在第$j$方向的尺寸。周期性边界条件同时保证了位移的周期性和牵引力的连续性，因此在周期性微结构材料的等效参数计算中被广泛使用。所得的本构关系保证了复杂微观结构在宏观热-力学分析中的等效性。基于Hill–Mandel能量一致性条件提取的等效参数，使得复杂微观结构可以被均质化材料模型所替代，同时仍能保持精确的整体力学响应。这显著减少了跨尺度有限元模拟中的自由度和计算成本。该方法对于含有重布线层、硅通孔和微焊点阵列等多尺度结构的三维非均匀封装系统的热-力学仿真尤为重要。

### 3. 结果与讨论

#### 3.1 基于分形的重布线层图形几何简化

在高性能计算单元中，图形处理器（GPU）和高带宽存储器（HBM）通过2.5D封装技术高度集成于中介层上，中介层再通过微凸点焊料连接至封装基板。GPU与HBM通过中介层内密集分布的微米级重布线层（RDL）实现信号传输，而通过封装基板中的重布线层与外围电路进行通信（图2a）。本文仅选取HBM及其关联电路进行仿真建模。图2b展示了一个典型的四芯片堆叠存储器结构。五层方形中介层的边长为33 mm、厚度为1 mm，而十层方形封装基板的边长为50 mm、厚度为2.4 mm。

![图2]

**图2.** 高带宽存储器结构及重布线层降阶建模策略。（a）HBM 2.5D封装结构示意图。（b）四芯片HBM封装几何模型。（c）中介层重布线层的典型原始图形及其特征尺寸。（d）根据铜覆盖率与特征线宽，基于Hilbert图形的原始图形重布线层降阶表示。

图2c展示了中介层中的典型重布线层，其中黑色区域代表空隙区域，蓝色区域对应铜走线。根据空隙率，重布线层可划分为两个特征图形区域。黄色方框标记的图形对应高空隙区域，其最小铜线宽为18 μm；而白色方框标记的图形对应低空隙区域，最小铜线宽为120 μm。可以看出，重布线层图形具有强烈的不规则性和显著的多尺度特征。为应对这一几何复杂性，本文提出了一种基于分形理论的降阶建模方法，将复杂原始图形根据其特征线宽和单位面积铜覆盖率等效简化为分形图形。通过这种方式，复杂无序的几何图形被转化为自相似结构，从中可以方便地提取代表性几何特征。

对于图2c所示的高空隙区域，计算得到的特征线宽为23 μm（2 × 总图形面积 / 总图形周长），单位面积铜覆盖率为16.76%（总图形面积 / 边框面积）。对于低空隙区域，计算得到的特征线宽为146 μm，单位面积铜覆盖率为65.18%。Hilbert曲线作为一种经典的分形曲线，能够以一条一维路径遍历二维区域，同时允许控制线宽和间距，因此特别适合于近似具有方形轮廓的重布线层中的走线分布图形。相应地，开发了一个MATLAB程序来生成具有指定线宽和单位面积线面积比（图形轮廓包围的总面积为1 cm2）的分形曲线单元。如图2d所示，构建了线宽和线面积比与原始重布线层图形一致的分形单元。通过将分形单元以与原始图形相同的总面积排列成阵列，获得了等效的降阶Hilbert图形层。

由于Hilbert图形层的几何特征参数与原始图形一致，等效降阶图形层继承了原始重布线层的力学影响。此外，分形图形的自相似特性允许提取最小分形单元，从而将初始的复杂力学问题分解为单个单元问题。只需对最小分形单元提取等效力学参数，从而显著降低计算成本。进一步地，当将最小分形单元得到的等效力学参数赋予大规模阵列化分形图形层时，所得力学性能仍与原始图形等效。

#### 3.2 基于最小二乘优化的等效热-力学参数提取

在获得与原始重布线层几何等效的分形图案最小自相似单元后，提出了一种基于最小二乘法的参数优化方法来提取原始图案的等效弹性模量和热膨胀系数。主要步骤如下：

(1) 首先构建了一个与原始分形单元相同边长（1 mm）和厚度（240 μm）的均质等效单元（图3a）。在两个单元中，顶层材料（厚40 μm）被定义为铜以代表重布线层，底层材料（厚200 μm）被定义为FR4以代表基板。在原始分形单元中，顶层铜的分布以Hilbert曲线的形式设计，以重现原始重布线层的几何特征；而在均质等效单元中，顶层铜的分布被定义为均匀层，以表示等效重布线层几何。此外，原始分形单元的底层被建模为包含多个通孔的方形板，而均质等效单元的底层被建模为无通孔的均匀方形板。这样，原始基板中通孔的几何特征即可通过参数优化过程等效地由均质基板表示。上述两种单元中各组件的材料属性列于表1。

![图3]

**图3.** 基于三点弯曲和热膨胀仿真的等效热-力学参数提取方法。（a）原始非均匀单元与等效均质单元的三点弯曲仿真边界条件。（b）原始非均匀单元与等效均质单元的热膨胀仿真边界条件。（c）原始非均匀单元与等效均质单元的三点弯曲仿真应力云图。（d）原始非均匀单元与等效均质单元的热膨胀仿真应力云图。（e）基于参数优化方法提取等效弹性模量的算法。（f）基于参数优化方法提取等效热膨胀系数的算法。（g）原始非均匀单元与等效均质单元在相同三点弯曲边界条件下底部中心线的位移曲线对比。（h）原始非均匀单元与等效均质单元在相同热膨胀边界条件下界面中线位移曲线的对比。

**表1.** HBM及其外围电路的材料属性。

| 类型 | 材料 | 密度(kg/m3) | 弹性模量(GPa) | 泊松比 | 热膨胀系数(10?? K?1) | 导热系数(W/m·K) | 比热容(J/kg·K) |
|------|------|-------------|--------------|--------|---------------------|-----------------|----------------|
| HBM | 硅 | 2329 | 170 | 0.28 | 2.6 | 130 | 700 |
| 中介层 | — | — | — | — | — | — | — |
| 基板 | FR4 | 1900 | 22 | 0.15 | 18 | 0.3 | 1369 |
| RDL | 铜 | 8960 | 110 | 0.35 | 17 | 400 | 385 |
| 通孔 | VIA | — | — | — | — | — | — |

(2) 以上述两种单元作为几何模型，在相同条件下进行三点弯曲仿真（图3a）和热膨胀仿真（图3b）。三点弯曲仿真的边界条件如图3a所示。对每个单元的两侧底边施加固定约束，同时在底面中线施加线载荷 $F=n\,\text{N}$（$n=50, 100, 150$）。沿垂直于加载方向的底面线定义为位移提取曲线，以获取弯曲变形数据。热膨胀仿真的边界条件如图3b所示。固定两个单元的底面，并对整个结构施加50 °C的均匀温升。沿上下层界面定义位移提取曲线以提取热膨胀变形数据。

(3) 以三点弯曲仿真中原始单元的底部中线位移曲线为参考，不断调整均质单元顶层材料的弹性模量，直到均质单元的底部中线位移曲线与原始单元的匹配（图3e）。此时，两个单元在三点弯曲条件下的应力云图如图3c所示，其z方向底部中线位移曲线如图3g所示。尽管由于等效重布线层均匀分布铜的更强约束作用，均质单元的位移曲线在加载区域外不能完全与原始单元重合，但两者在加载区域内的良好吻合表明，均质单元的位移曲线能够近似重现原始单元在相同机械载荷下的最大变形特征。此时均质单元的弹性模量 $E=0.081\times110\,\text{GPa}$，即被取为原始重布线层的均质等效值。类似地，以热膨胀仿真中原始单元的界面中线位移曲线为参考，不断调整均质单元顶层材料的热膨胀系数，直到均质单元的界面中线位移曲线与原始单元匹配（图3f）。此时，两个单元的热膨胀应力云图如图3d所示，界面中线z方向位移曲线如图3h所示。尽管由于非均匀走线图形引起的局部应力集中，原始单元的界面位移曲线表现出更显著的峰值，导致两条曲线在细节上存在差异，但两者的整体吻合表明，均质单元的位移曲线能够重现原始位移响应的总体趋势。该近似足以在相同热膨胀条件下实现两个单元等效的整体z方向变形。此阶段均质单元的热膨胀系数 $\alpha=1.7\times10^{-5}\,\text{K}^{-1}$，即被取为原始重布线层的均质等效值。

将铜的实际弹性模量（$E=110\,\text{GPa}$）赋予以分形图形表示重布线层的单层中介层模型，有限元仿真得到的三点弯曲应力分布如图4a所示（底边载荷为$F=1\,\text{N}$）。使用相同的弹性模量（$E=110\,\text{GPa}$），构建第二个模型，其中重布线层以均质图形表示，并在相同加载条件下进行有限元仿真。

对比两个模型底部中线的 z 方向位移曲线，揭示出它们在弯曲变形行为上存在显著差异（图 4c）。随后，通过基于最小二乘法的参数优化算法获得的等效杨氏模量（E = 0.1×110 GPa）被赋予均质重布线层模型。相应的三点弯曲应力分布如图 4b 所示。在此阶段，两个模型底部中线的 z 方向位移曲线表现出极好的一致性（图 4d），表明基于最小二乘法的等效参数优化方法可以有效地用简化的均质表示取代复杂的微观结构图案，同时保持宏观力学等效。

**图 4. 基于分形模型和均质模型的单层中介层三点弯曲仿真。**
(a) 用分形图案表示重布线层的单层中介层建模方法及其对应的三点弯曲应力云图。
(b) 用等效均质图案表示重布线层的单层中介层建模方法及其对应的三点弯曲应力云图。
(c) 均质图案等效参数优化前两个模型三点弯曲位移曲线的对比。
(d) 均质图案等效参数优化后两个模型三点弯曲位移曲线的对比。

### 3.3. 基于壳的多层中介层和衬底降阶建模

在通过参数优化方法将原始单元中复杂的重布线层图案等效替换为均匀单元中的均质化几何形状后，复杂的原始模型可以被简化为简化的均质模型。如图 5a 所示，均质模型中的中介层和封装衬底均采用壳 - 板组合建模方法构建。蓝色平面代表降阶均质重布线层，而灰色板对应力学与原始模型尺寸相同的衬底层。中介层由五层“重布线层 - 衬底”复合材料层组成，其中重布线层厚度为 2 μm，衬底厚度为 200 μm。采用基于壳的表示方法，将超薄 RDL 的机械响应映射到二维边界上，同时保留其膜力和弯曲贡献。

**图 5. 基于壳理论的薄铜层和窄过孔建模方法。**
(a) 中介层和封装衬底中重布线层的壳建模。
(b) 中介层中硅通孔和封装衬底中过孔的壳建模。

这种基于壳的公式非常适合二维表面上的有限元离散化。与完整的三维实体建模相比，该方法在保持面内机械响应准确性的同时，显著减少了自由度数量，从而大大提高了数值模拟的效率和稳定性。通过这种方法，五层中介层的总厚度仅等于五层硅衬底的组合厚度（1 mm）。五层重布线层由硅衬底的上表面以及相邻硅层之间的界面定义。在 COMSOL 中，这些界面使用"Shell"（壳）物理接口定义为"Interposer RDL shells"（中介层 RDL 壳），壳厚度为 2 μm，位置指定为"mid-surface on boundary"（边界中表面）。通过添加"Solid–Thin Structure Connection"（固体 - 薄结构连接）多物理场耦合，"Interposer RDL shell"（中介层 RDL 壳）界面与相邻的硅衬底被定义为通过共享边界连接。同样，十层封装衬底使用相同的方法建模，其总厚度等于十层 FR4 衬底的厚度（2 mm）。十层重布线层由 FR4 板的上表面以及相邻 FR4 层之间的界面定义。在 COMSOL 中，这些界面通过"Shell"（壳）物理接口定义为"Package substrate RDL shells"（封装衬底 RDL 壳），壳厚度为 40 μm，位置设置为"mid-surface on boundary"（边界中表面）。再次引入"Solid–Thin Structure Connection"（固体 - 薄结构连接）多物理场耦合，使"Package substrate RDL shells"（封装衬底 RDL 壳）界面与其相邻的 FR4 层共享边界。这样，无需复杂的重布线层图案或跨尺度薄层结构，即可构建中介层和封装衬底的降阶模型。如图 5a 所示，虽然重布线层在几何模型中未被赋予真实厚度，但通过 COMSOL 中的基于薄壳理论的后处理，仍可获得具有实际厚度的 RDL 应力云图。

中介层中的 TSV 和封装衬底中的过孔也使用“壳方法”建模。TSV 是直径为 10 μm 的实心铜柱，穿透硅中介层；而封装衬底中的过孔是外径为 50 μm、内径为 30 μm 的空心铜柱，穿透衬底。这两种类型的铜柱都可以由在中介层和衬底中创建相应外径孔时生成的圆柱壁面定义。中介层中的孔壁面在 COMSOL 中通过"Shell"（壳）物理接口定义为"Interposer TSV copper shells"（中介层 TSV 铜壳），壳厚度为 5 μm，位置指定为"bottom surface on boundary"（边界底表面）。由于孔的外径为 10 μm，此壳厚度使弯曲的铜壳聚合成等效的实心铜柱。通过添加"Solid–Thin Structure Connection"（固体 - 薄结构连接）多物理场耦合，"Interposer TSV copper shell"（中介层 TSV 铜壳）界面与相邻的硅中介层被定义为通过共享边界连接。同样，封装衬底中的孔壁面通过"Shell"（壳）物理接口定义为"Package substrate VIA copper shells"（封装衬底过孔铜壳），壳厚度为 20 μm，位置设置为"top surface on boundary"（边界顶表面）。因为孔的外径为 50 μm，此壳厚度允许弯曲的铜壳形成内径为 30 μm 的空心铜柱。通过引入"Solid–Thin Structure Connection"（固体 - 薄结构连接）多物理场耦合，"Package substrate VIA copper shell"（封装衬底过孔铜壳）界面与相邻的封装衬底通过共享边界连接。这样，无需显式构建微观柱体或跨尺度弯曲薄层，即可实现过孔的降阶建模。如图 5b 所示，虽然 TSV/VIA 孔壁在几何模型中未被赋予真实厚度，但通过 COMSOL 中基于薄壳理论的后处理，仍可获得具有实际厚度的 TSV 和过孔的应力云图。如图 5b 所示，HBM 芯片与中介层以及中介层与封装衬底通过焊点连接。在 COMSOL 中，焊点和板在固体力学物理接口中定义为柔性连接。

### 3.4. RVE 均质化及等效弹性和热性质的提取

上述基于参数优化和壳建模的降阶建模方法为下文基于 RVE 理论的降阶仿真提供了参考。以封装衬底的基于 RVE 的降阶建模为例，本研究选择的 RVE 如图 6a 所示。从尺寸为 53.55 mm、总厚度为 2.4 mm 的十层封装衬底的中心区域提取了一个边长为 3.825 mm 的方形单元。该单元由厚度为 200 μm 的 FR4 层和线宽为 23 μm 的厚度为 40 μm 的重布线铜层组成。在 FR4 层的中心有一个空心铜柱，外径为 50 μm，内径为 30 μm。

**图 6. 基于双层 RVE 模型的等效热膨胀系数和等效弹性刚度矩阵提取方法。**
(a) 代表性体积单元的选取方法。
(b) 沿 x 方向施加单轴拉伸应变时 RVE 的应力分布。
(c) 沿 y 方向施加单轴拉伸应变时 RVE 的应力分布。
(d) 沿 z 方向施加单轴拉伸应变时 RVE 的应力分布。
(e) 在 XY 平面施加剪切应变时 RVE 的应力分布。
(f) 在 YZ 平面施加剪切应变时 RVE 的应力分布。
(g) 在 XZ 平面施加剪切应变时 RVE 的应力分布。
(h) 底面固定，温度均匀升高 50 °C 自由热膨胀时 RVE 的应力分布。
(i) 使用 RVE 方法获得的均质单元的等效弹性模量矩阵。

RVE 均质化是通过宏观尺度和微观尺度之间的一致性实现的，这意味着 RVE 的宏观应变能等于其微观结构的总应变能。通过等效弹性矩阵关系，RVE 的宏观平均应力 - 应变关系变得与均质材料等效。等效弹性矩阵是宏观均质材料的关键参数，定量描述了宏观应力 $\sigma$ 和宏观应变 $\varepsilon$ 之间的关系。对于三维各向异性材料，等效弹性矩阵采取对称 6×6 矩阵的形式：

$$ (21) \ \sigma = C \cdot \varepsilon $$

$$ (22) \ C = \begin{bmatrix} C_{11} & C_{12} & C_{13} & C_{14} & C_{15} & C_{16} \\ C_{12} & C_{22} & C_{23} & C_{24} & C_{25} & C_{26} \\ C_{13} & C_{23} & C_{33} & C_{34} & C_{35} & C_{36} \\ C_{14} & C_{24} & C_{34} & C_{44} & C_{45} & C_{46} \\ C_{15} & C_{25} & C_{35} & C_{45} & C_{55} & C_{56} \\ C_{16} & C_{26} & C_{36} & C_{46} & C_{56} & C_{66} \end{bmatrix} $$

其中 C 表示由单位 $j$ 方向应变引起的 $i$ 方向应力，例如，$C_{11}$ 是方向 1 上的拉伸刚度。

可以看出，基于 RVE 理论的降阶问题本质上是求解代表性单元的等效弹性矩阵的问题。从原始对象中提取的代表性单元中包含的微观结构的本构参数（密度、杨氏模量、泊孙比和热膨胀系数）是先验已知的。因此，在有限元仿真中，对代表性单元施加给定量和方向的指定应变会产生相应的应力响应。当应力和应变均已知时，可根据方程 (21) 和方程 (22) 直接确定等效弹性刚度。为了解决 6×6 等效弹性刚度矩阵，必须对 RVE 施加六个独立的宏观应变分量（三个正应变和三个剪切应变），并计算相应的微观应力响应。具体的加载方案如下：

对于图 6a 所示的代表性单元，在 COMSOL 的固体力学物理接口内，重布线层和 FR4 板定义为柔性连接，而空心铜柱和 FR4 板定义为刚性连接。然后通过将顶面和底面、左面和右面以及前面和后面定义为三对周期性边界，引入周期性边界条件。随后，施加表 2 中列出的六种应变状态。相应的平均应变矩阵由下式给出：

$$ (23) \ \varepsilon = \begin{bmatrix} \varepsilon_{11} & 0.5\varepsilon_{12} & 0.5\varepsilon_{13} \\ 0.5\varepsilon_{12} & \varepsilon_{22} & 0.5\varepsilon_{23} \\ 0.5\varepsilon_{13} & 0.5\varepsilon_{23} & \varepsilon_{33} \end{bmatrix} $$

**表 2. 提取等效弹性刚度矩阵的应变加载方案。**

| 加载顺序 | 应变分量 | 应变模式 | 物理含义 | 加载目的 |
| :--- | :--- | :--- | :--- | :--- |
| 1 | $\varepsilon_{11}$ | 正应变 | X 方向拉伸/压缩 | 求解 $C_{11}, C_{21}, C_{31}, C_{41}, C_{51}, C_{61}$ |
| 2 | $\varepsilon_{22}$ | 正应变 | Y 方向拉伸/压缩 | 求解 $C_{12}, C_{22}, C_{32}, C_{42}, C_{52}, C_{62}$ |
| 3 | $\varepsilon_{33}$ | 正应变 | Z 方向拉伸/压缩 | 求解 $C_{13}, C_{23}, C_{33}, C_{43}, C_{53}, C_{63}$ |
| 4 | $\varepsilon_{12}$ | 剪切应变 | X-Y 平面剪切 | 求解 $C_{14}, C_{24}, C_{34}, C_{44}, C_{54}, C_{64}$ |
| 5 | $\varepsilon_{13}$ | 剪切应变 | X-Z 平面剪切 | 求解 $C_{15}, C_{25}, C_{35}, C_{45}, C_{55}, C_{65}$ |
| 6 | $\varepsilon_{23}$ | 剪切应变 | Y-Z 平面剪切 | 求解 $C_{16}, C_{26}, C_{36}, C_{46}, C_{56}, C_{66}$ |

图 6b-d 展示了分别沿 x、y 和 z 方向对代表性单元施加单轴拉伸正应变时获得的应力云图。图 6e-g 展示了分别在 XY、YZ 和 XZ 平面上施加剪切应变时的应力云图。通过这些仿真，获得了相应的 6×6 等效弹性刚度矩阵，如图 6i 所示。将此弹性矩阵作为与原始封装衬底尺寸相同的均质衬底的力学性能，均质衬底的三点弯曲仿真可以重现原始衬底的应力响应和弯曲变形，从而完成基于 RVE 的三点弯曲问题的降阶仿真。同样，在热应力分析中，热膨胀系数描述了热应变与温度变化 $\Delta T$ 之间的关系。对于 3D 各向异性材料，CTE 表示为对称 3×3 矩阵：

$$ (24) \ \varepsilon_{th} = \alpha \cdot \Delta T $$

$$ (25) \ \alpha = \begin{bmatrix} \alpha_{xx} & \alpha_{xy} & \alpha_{xz} \\ \alpha_{yx} & \alpha_{yy} & \alpha_{yz} \\ \alpha_{zx} & \alpha_{zy} & \alpha_{zz} \end{bmatrix} $$

为了解决等效 3×3 矩阵，周期性边界条件也应用于 COMSOL 固体力学物理接口中的代表性单元，其中顶面和底面、左面和右面以及前面和后面被定义为三对周期性边界。

随后，对该单元施加自由热膨胀应变以执行有限元仿真，获得对应的热应力响应。图6h展示了在底部表面固定、整体升温50°C条件下代表性单元的应力云图。根据该应力结果以及所施加的温度升高，可以反推出等效热膨胀系数（CTE）矩阵为：(26) α = [[1.8×10??, 5.25×10?1?, -1.54×10?2?], [5.25×10?1?, 1.80×10??, -5.27×10?2?], [-1.54×10?2?, -5.27×10?2?, 1.54×10??]] 可以看出，上述等效CTE矩阵的非对角项趋近于零，表明该代表性单元表现为典型正交各向异性材料，这与PCB结构的分层材料分布特征相一致。将该等效CTE矩阵连同此前获得的等效弹性矩阵一并赋予与原封装基板相同尺寸的均质基板作为其力学性能，即可通过该均质基板的热膨胀仿真复现原基板的应力响应与膨胀变形。至此，基于RVE理论的降阶热膨胀仿真即告完成。

类似地，图7展示了另一种夹层型RVE模型。该单元由两层厚度各为1.2 mm的FR4层夹持一层厚度为40 μm的再布线铜层构成，并在FR4层中心引入一个内径30 μm、外径50 μm的过孔。按照与双层RVE模型相同的流程，可以提取夹层型RVE的等效弹性刚度矩阵和有效CTE。

下载：下载高分辨率图片（1MB） 下载：下载全尺寸图片
图7. 基于夹层RVE模型的等效热膨胀系数与等效弹性刚度矩阵提取方法。(a) 代表性体积元的选取方式。(b)-(d) RVE在x方向(b)、y方向(c)和z方向(d)施加单轴拉伸应变时的应力分布。(e)-(g) RVE在XY平面(e)、YZ平面(f)和XZ平面(g)施加剪切应变时的应力分布。(h) RVE在底部固定、整体升温50°C条件下自由热膨胀时的应力分布。(i) 利用RVE方法获得的均质单元等效弹性模量矩阵。

由夹层型RVE（图8a）和双层RVE（图8d）推导出的等效弹性刚度矩阵，分别被赋予边长50 mm、厚度2.4 mm的均质板模型，用于对单层封装基板进行三点弯曲仿真（边缘载荷F=100 N）。同时，在相同条件下对相同结构使用原始分形模型进行三点弯曲仿真作为参考。尽管基于夹层型RVE的均质模型所获得的应力分布与原始分形模型较为吻合（图8b），但z方向位移仍存在显著差异（图8c）。两种RVE构型的差异可归因于其不同的约束与载荷传递条件。在夹层型RVE中，RDL被两侧FR4层包裹，对其变形施加了强双边约束。这种构型抑制了图案化铜层的力学贡献，并改变了RDL与周围基板之间的载荷传递路径。因此，夹层型RVE推导出的等效刚度无法完全复现原始分形结构的变形特性。相比之下，双层RVE仅包含单一RDL-基板界面，避免了这种人为的双边约束，使RDL与基板之间的变形协调和载荷传递更接近于原始封装基板中相应的重复结构单元。因此，双层RVE能够更准确地保留RDL的有效力学贡献，使其在应力分布和位移响应两方面均与参考模型表现出显著改善的吻合度。结果发现，在相同边缘载荷下，均质模型的变形大于原始分形模型。相反，基于双层RVE的均质模型在整体应力分布（图8e）和z方向位移（图8f）两方面均与原始分形模型表现出一致性，表明双层RVE能更好地捕捉原始系统的基本微观结构特征，使得对应的均质模型能够高保真地复现原始分形模型的力学响应。因此，在所有后续基于代表性体积元法的降阶仿真中，均采用双层RVE选取策略。

下载：下载高分辨率图片（535KB） 下载：下载全尺寸图片
图8. 基于双层型和夹层型RVE模型的单层封装基板三点弯曲仿真。(a) 夹层型RVE模型。(b) 基于夹层型RVE的均质模型与原始分形模型的三点弯曲应力云图对比。(c) 基于夹层型RVE的均质模型与原始分形模型的三点弯曲位移曲线对比。(d) 双层RVE模型。(e) 基于双层RVE的均质模型与原始分形模型的三点弯曲应力云图对比。(f) 基于双层RVE的均质模型与原始分形模型的三点弯曲位移曲线对比。

**3.5 热力学载荷下降阶模型的精度与效率验证**

上述两种降阶建模方法随后被应用于相同的HBM 2.5D封装结构，在相同条件下进行三点弯曲和热膨胀仿真，以交叉验证其建模精度和计算资源需求。利用基于参数优化和壳体建模的降阶建模方法，可以构建原始等效模型（图9a）。在该模型中，中介层的五个硅层及层间TSV，以及封装基板的十个FR4层及其中的过孔，均按照实际结构进行显式建模。相邻板层之间的再布线层则被参数优化获得的均质壳模型所替代。利用基于代表性体积元理论的降阶建模方法，可以构建图9b所示的RVE降阶模型。在这种情况下，中介层和封装基板被建模为总厚度分别为1 mm和2 mm的均质单板。RVE分析获得的等效弹性刚度矩阵和等效热膨胀系数矩阵被赋予这些均质板作为其力学性能参数，从而考虑了原始电路中再布线层、TSV和过孔等微观结构的力学效应。如图9a、b所示，原始等效模型和RVE降阶模型均由包含HBM芯片、中介层和封装基板的三层结构组成。相邻层之间通过焊球连接，其尺寸和形态如图5b所示。在HBM芯片与中介层之间选取四个代表性焊球t1-t4，在中介层与封装基板之间选取另外四个焊球b1-b4，这八个焊球被选为应力提取点。

下载：下载高分辨率图片（938KB） 下载：下载全尺寸图片
图9. RVE降阶模型在三点弯曲和热膨胀条件下与原始等效模型的验证对比。(a) 原始等效模型的三点弯曲仿真边界条件。(b) RVE降阶模型的三点弯曲仿真边界条件及代表性焊球位置。(c) 原始等效模型与RVE降阶模型三点弯曲仿真应力云图对比。(d) 不同边缘载荷下原始等效模型与RVE降阶模型的底部中线位移曲线对比。(e) 原始等效模型与RVE降阶模型热膨胀仿真应力云图对比。

图9a、b展示了三点弯曲仿真的边界条件。对四芯片HBM 2.5D封装结构的横向和纵向截面施加对称边界条件，从而构建对应实物四分之一的几何模型。在纵向截面底部边缘施加线载荷F=n N（n=50, 100, 150）。与加载线平行的封装基板顶部边缘被定义为固定约束。横向截面底部边缘被定义为位移提取线，用于获取板的三点弯曲变形曲线。如图9c所示，基于原始等效模型和RVE降阶模型的三点弯曲仿真所得应力云图高度一致。图9d展示了两个模型在不同边缘载荷（F=50, 100, 150 N）下底部中线的z方向位移曲线。两条位移曲线之间的均方根误差在载荷为50、100和150 N时分别计算为0.00396、0.00792和0.0182，表明曲线吻合良好。

对于热膨胀仿真，原始等效模型和RVE降阶模型的底部表面被固定（图9a、b），并对整体施加均匀温度变化ΔT = m ℃（m = -50, 50, 125）。所得热膨胀应力云图如图9e所示，两个模型的应力分布也高度相似。图9f展示了不同温度条件下底部中线的z方向位移曲线。温度变化为-50、50和125°C时，两个模型位移曲线之间的均方根误差分别为1.95×10??、8.34×10??和2.92×10??，表明曲线吻合度优异。图9g展示了两种模型在不同边缘载荷条件下芯片底部焊球处的应力集中情况。提取这些焊球的体积平均应力，并以原始等效模型三点弯曲仿真中的焊球应力作为参考值。根据式(27)，可以计算RVE降阶模型三点弯曲条件下的焊球应力仿真精度。如图9h所示，b2、b3和t2焊球的应力预测误差均低于2%。表3列出了两种降阶模型三点弯曲仿真的计算时间和内存消耗。与原始等效模型相比，RVE降阶模型的求解时间和内存消耗分别减少了90.2%和78.1%。图9i展示了两种模型在不同温度条件下芯片底部焊球处的应力集中情况。提取这些焊球的体积平均应力，并以原始等效模型热膨胀仿真中的焊球应力作为参考值。根据式(27)，可以计算RVE降阶模型热膨胀条件下的焊球应力仿真精度。如图9j所示，b1、b4和t4焊球的应力预测误差均低于1%。表4列出了两种降阶模型热膨胀仿真的计算时间和内存消耗。与原始等效模型相比，RVE降阶模型的求解时间和内存消耗分别减少了84.0%和79.2%。本次验证旨在量化模型降阶所引入的误差，而非建立针对特定制造HBM封装的实验标定模型。因此，采用更高保真度的原始等效模型作为数值参考，以评估RVE降阶模型是否保留了主要的热力学响应。尽管两个建模层级之间的密切吻合证明了所提出降阶策略的有效性，但在对特定HBM产品进行定量预测之前，仍需要使用已制造的封装结构进行实验验证，因此这将构成未来工作的重要方向。

(27) Error_sim = (σ_RVE - σ_origin) / σ_origin × 100%

表3. 两种降阶模型三点弯曲仿真的计算时间与内存消耗对比。

| 三点弯曲仿真 | 求解时间(s) | 内存占用(GB) | 单元数 |
|---|---|---|---|
| 原始等效模型 | 231930 | 1.58467 | 059 |
| RVE降阶模型 | 基板127663 | 36646 | |
| | 中介层2522 | 18727 | 25 |
| 总计/最大值 | 75663 | 73913 | |
| 总和 | 22717 | 4.84 | / |

表4.

两种降阶模型热膨胀仿真计算时间与内存消耗对比

热膨胀仿真 | 求解时间(s) | 内存使用(GB) | 单元数
原始等效模型 | | |
总数 | 170829 | 1.984 | 67059

RVE降阶模型 | | |
基板 | 1406 | 0.643 | 36736
中介层 | 4025 | 4.872 | 725
总数/最大值 | 936 | 0.643 | 73913
合计 | 27317 | 1.91 |

(f) 不同整体温升下原始等效模型与RVE降阶模型底部中心线位移曲线对比。(g) 不同边缘载荷下原始等效模型与RVE降阶模型芯片底部焊球应力分布对比。(h) 不同边缘载荷下RVE降阶模型预测焊球应力的仿真误差，以相同条件下原始等效模型的焊球应力为参考值。(i) 不同整体温升下原始等效模型与RVE降阶模型芯片底部焊球应力分布对比。(j) 不同整体温升下RVE降阶模型预测焊球应力的仿真误差，以相同条件下原始等效模型的焊球应力为参考值。

4. 结论

针对HBM 2.5D封装结构跨尺度显式建模导致的热机耦合仿真中计算成本高、仿真时间长、收敛困难等挑战，本文提出了一种基于分形理论和壳模型的降阶建模方法。采用该方法构建的原始等效模型能够准确捕捉整体结构中的系统级行为，包括工作条件下温度分布引起的全局翘曲、局部弯曲和界面应力演化。同时，该方法显著缩短了计算时间并减少了计算资源消耗，使热机耦合仿真从一次性验证工具转变为迭代工程设计与优化的预测工具。此外，该方法弥补了基于全细节建模的参考仿真结果的缺失，为评估基于RVE理论的降阶仿真的计算精度和效率提供了定量基准。

以原始等效模型的仿真结果为参考，双层RVE降阶模型预测的三点弯曲位移曲线与参考曲线之间的均方根误差(RMSE)低于0.02。对于热膨胀仿真，位移曲线之间的RMSE低至10??数量级。在不同边缘载荷条件下，三点弯曲仿真中焊球应力的预测误差小于2%，而在不同温度条件下的热膨胀仿真中误差低于1%。此外，与原始等效模型相比，RVE降阶模型将三点弯曲和热膨胀仿真的求解时间分别减少了90.2%和84.0%，内存消耗分别降低了78.1%和79.2%。计算成本的降低促进了先进封装系统的重复设计迭代、参数研究和大规模可靠性评估。此外，更短的求解器运行时间和更低的内存需求可以减少相关电力消耗，从而有助于更节能和可持续的计算设计。

这些结果表明，原始等效方法和RVE降阶方法在仿真精度和模型降阶有效性方面相互补充和验证。尽管原始等效仿真比RVE降阶仿真需要更多的计算资源，但它能够提供再分布层和过孔层等微米尺度微观结构内部的详细应力分布。相比之下，RVE降阶仿真以显著更低的计算成本和时间消耗，实现了宏观板级整体应力和变形的快速预测。因此，本研究为评估极端AI工作负载下HBM单元的结构可靠性提供了可行途径，并为高性能计算时代先进异构集成系统的稳定运行和寿命设计提供了重要的理论和方法论支持。

CRediT作者贡献声明
Yanbo Du：原创草稿撰写、可视化、验证、方法论、形式分析、数据整理。Chenpeng Liu：可视化、数据整理。Wei Du：软件、调查。Daili Feng：项目管理。Yanhui Feng：监督、项目管理。Lin Qiu：审阅与编辑、验证、资源、项目管理、调查、经费获取、原创草稿撰写。

Yanbo Du | Chenpeng Liu | Wei Du | Daili Feng | Yanhui Feng | Lin Qiu