直接多級有限元法在多尺度封裝結構中的應用

趙勝軍 公顏鵬 侯傳濤 秦飛

（1電子封裝技術與可靠性研究所，北京工業大學，北京 100124；2 可靠性與環境工程技術重點實驗室，北京強度環境研究所，北京 100076）

0 引言

電子封裝不斷向輕質量，多功能，高性能和高集成度方向發展[1]，其可靠性問題[2-4]一直是研究的熱點。封裝結構的特征尺寸一般會相差多個數量級[5]，表現出明顯的結構多尺度特征，例如，晶圓級封裝中的TSV結構[6]，RDL層[7]等。結構多尺度特征給封裝結構仿真分析帶來了極大的困難，這也對當前封裝結構多尺度仿真分析方法提出了新的需求和挑戰。

為了解決封裝結構多尺度仿真的問題，研究人員提出了各種分析方法。根據仿真模型是否為結構和載荷對稱的特點，將模型簡化為1/4，1/8模型或GPD模型[8-9]，但是經過這些簡化后，涉及到幾何多尺度的仿真模型時，其網格數量仍然很多，計算量仍然會大，甚至有時無法計算。

封裝結構中的異質材料通常會使用均勻化方法[10]將其等效成均質材料。秦飛[11]等根據復合材料力學分析方法，把三維硅通孔轉接板封裝結構中芯片與轉接板間的微凸點/下填料層以及轉接板與基板間的微焊點/下填料層等效為線彈性均質材料，對等效模型進行熱疲勞壽命的仿真分析。Cheng等[12]利用代表體單元（Representative Volume element，RVE）預測了復合材料的等效彈性參數。Omairey等[13]針對周期性RVE結構開發了一款Abaqus插件，用于預測等效彈性參數。以上研究只考慮了材料的線彈性參數，然而對于封裝結構中表現出非線性行為的材料可能會帶來較大誤差。這就需要非線性多尺度分析方法對異質材料結構進行分析。

Tan[14]提出了DFE法，它是一種非線性多尺度方法，能夠分析材料的非線性行為。該方法主要針對二維周期性分布的結構，無法直接應用于復雜的實際封裝結構多尺度結構的仿真。因此，本文將對DFE方法進行改進，以適用于復雜的封裝結構多尺度結構仿真分析。

本文基于DFE方法，提出了一種能夠用于封裝結構多尺度仿真的DFE-子模型方法。采用該方法對封裝結構進行了有限元仿真，將仿真結果與全模型結果、子模型和局部均勻化方法的結果進行對比，說明了該方法的準確性。

1 DFE理論及其在Abaqus中的實現

1.1 DFE理論

平衡方程的弱形式[12]:

其中，u是位移；σ是應力張量；b是體力；t是面力；V、S指計算域和邊界。

式(1)是虛功原理的表達式，即內力虛功 intWδ等于外力虛功extWδ。

有限元分析中，通常采用高斯積分法進行數值積分。因此，

其中，α是單元e的高斯積分點。Jα是雅可比行列式；ωα是高斯點的權重。

在DFE中，每個高斯點的應力是由相應RVE的體積平均應力計算得到的。因此方程(2)可以寫成：

其中，〈·〉α 表示與單元e內的高斯點α相關的RVE上的體積平均量。“～”用于表示微觀尺度計算中的量。

Hill-Mandel均質化條件要求，

聯立公式(3)和(4)得

把（5）式代入（1）式，得

式（6）的左邊是微觀的量，右邊是（1）式中的宏觀描述。

經有限元離散化，式（6）在有限元中表示為，

將(8)代入(7)并消去虛位移，得到

1.2 DFE在Abaqus中的實現

DFE方法在Abaqus中的實現可以通過兩個關鍵步驟完成。

1)通過Abaqus中的多點約束（multipoint constraints，MPC），建立微觀RVE邊界網格節點d和宏觀網格節點之間的關系。實現了從宏觀尺度到微觀尺度的過渡。

圖1 宏微觀單元Fig.1 Macro and micro scale element

RVE邊界AB和邊界CD節點之間的約束方程可以寫成：

為了約束 RVE 的剛體平移，需要施加額外的條件，

其中，0x指在宏觀單元內的高斯點，也是RVE中心的位置。

通過Abaqus中的MPC將約束方程施加到模型上，從而實現d和之間的聯系。

2）縮放RVE的剛度矩陣以獲得所有微觀RVE組成的整體剛度矩陣。

縮放系數可以直接確定。例如，對于二維DFE分析，使用2×2高斯積分點的矩形單元進行有限元分析時

為了更清晰的描述縮放系數與DFE有限元模型之間的關系，圖2給出了在宏觀單元內，不同縮放系數所對應的有限元模型。模型中RVE的體積Vα一般是固定的，因此，改變縮放系數的大小，根據公式（15）可知，宏觀單元體積也隨之變化。

圖2 不同縮放系數所對應的有限元模型Fig. 2 Finite element model for differentα

2 數值模擬

建立封裝結構中的焊盤和填充料（Cu-underfill）二維多尺度有限元模型并采用DFE方法進行仿真分析。均勻化方法和子模型技術也在算例中給出，用于與DFE方法的結果比較和分析。為了驗證計算結果的精度，本文用一個細化的全模型仿真結果作為參考解來計算不同方法的相對誤差。

2.1 DFE模型仿真分析

算例中Cu-underfill多尺度結構全模型有3200個RVE（40*80），RVE網格模型如圖3所示。RVE的長度和寬度為0.10.1mmmm× ；Cu材料在RVE的中間區域，Cu的中心與RVE中心重合，Cu的直徑為0.05mm；為了保證仿真結果的準確性，對全模型進行了網格收斂性驗證，當RVE的網格單元數量為457時，仿真結果趨于穩定，因此在后續的仿真中均采用該網格密度。

圖3 RVE網格模型Fig. 3 RVE mesh model

2.1.1 DFE仿真模型的建立

DFE仿真模型需要把RVE模型與宏觀尺度模型疊加在一起。宏觀單元中每個高斯積分點上均對應一個RVE，且RVE中心與高斯積分點重合。采用MPC的方式建立RVE邊界節點與宏觀單元節點之間的約束，所形成的DFE仿真模型如圖4所示，該模型的縮放系數為4。DFE模型的網格單元數量為365800，單元類型為CPS4R。

圖4 DFE仿真分析模型Fig.4 DFE model

2.1.2 材料參數

在宏觀尺度上，整個異質材料被離散化為均質的連續有限單元，宏觀尺度模型的計算不需要均勻的本構關系，這是因為整體剛度矩陣完全是由RVE的剛度矩陣乘以縮放系數得到。RVE 中Cu和underfill材料參數如表1所示，underfill為線彈性材料，Cu為非線性材料，屈服應力225MPa，切線模型6666MPa。由于宏觀尺度上不需要均勻的本構關系，因此，宏觀模型材料參數只需要設置較小的彈性模量即可。

表1 材料屬性Table 1 Material properties

2.1.3 施加邊界條件和載荷并求解

DFE模型中左邊界完全固定約束，在右邊界所有節點上施加單向拉伸載荷， 2xmm= 。A區域的von Mises 應力和等效塑性應變分布云圖如圖5所示，von Mises應力值為1001MPa，等效塑性應變為0.09699。仿真分析需要的運行時間為4.8min，所需內存576MB。將該模型仿真結果（S11、S22、S12、E11、E22和E12）與全模型相應結果進行了對比，如圖6和圖7所示。由圖可知，兩個模型的仿真結果較為吻合，說明了DFE方法是有效的。

圖5 DFE模型A區域仿真結果Fig.5 DFE model simulation results in area A

圖6 全模型與DFE模型不同應力結果對比Fig.6 The results of different stress for DFE and Full model

圖7 全模型與DFE模型不同應變結果對比Fig.7 The results of different strain for DFE and Full model

2.2 不同多尺度方法仿真結果的對比

采用子模型和局部均勻化方法[15]對2.1節中包含3200個RVE的Cu-underfill二維有限元模型進行仿真分析，RVE的網格數量、單元類型、邊界條件和載荷與DFE模型一致。所得結果用于與DFE方法的結果進行比較。把全模型的仿真結果作為參考解用于不同方法相對誤差的計算。不同方法的有限元仿真結果如表2所示，分析比較了仿真過程中的網格數量、內存、運行時間以及von Mises 應力和等效塑性應變的結果。不同多尺度方法對比分析的結果如下：

表2 不同多尺度方法仿真結果的對比Table 2 Comparison of simulation results for different multi-scale methods

a）子模型方法所需要的內存最少，結果的精度介于局部均勻化和DFE方法之間。

b）局部均勻化方法所得結果誤差最大，且網格數量仍然較多，所需內存也較多，在三種方法中效果最不理想。

c）DFE方法是預測精度最高，所需網格數量和運行時間最少的方法，且內存占比較少，在三種方法中效果最為理想。

3 DFE-submodel多尺度方法

DFE方法主要針對二維周期性分布的結構，不能夠直接應用于復雜的實際封裝結構多尺度的仿真。因此，基于DFE方法和子模型技術，發展了一種DFE與子模型結合（DFE-submodel）的多尺度方法。該方法首先利用DFE方法得到等效非線性參數，然后，對等效非線性均質化全局模型進行仿真分析，最后結合子模型技術得到模型中關鍵區域的應力應變值。

3.1 DFE法非線性等效參數的計算

首先計算所受拉力邊界上所有節點的位移值和節點反力，得到模型的載荷位移曲線；然后根據公式（16）和（17）得到名義應力應變曲線；最后根據公式（18）和（19）得到真實應力應變曲線。本算例，選擇不同的縮放系數對DFE方法進行仿真，采用細化的全模型結果作為參考解，所得應力應變曲線如圖6所示。根據公式（20）和得到的有限元結果可求得等效泊松比，本算例中，Cu-underfill的等效泊松比為0.3364。

由圖8可知，DFE模型與全模型結果吻合度較高，能夠較好的得到Cu-underfill等效應力應變關系，不同縮放系數的應力應變結果差別較小。根據應力應變關系可以得到等效的彈性模量為3.84GPa，將得到的彈塑性參數作為后續均質化模型的材料參數。

圖8 Cu-underfill等效應力應變曲線Fig.8 Equivalent stress-strain relationship for Cu-underfill

3.2 DFE-submodel仿真分析

均質化等效模型仿真結束后，選取均質化模型中關鍵區域進行子模型仿真。子模型結構中包含RVE模型，從而可以實現對模型中微觀區域仿真結果的提取。為了驗證該方法的有效性，采用DFE-submodel方法對Cu-underfill結構（2.1節中的全模型）進行仿真分析。對比了18個不同區域的最大von Mises應力結果，不同區域的位置如圖9所示。并且與全模型和局部均勻化模型的結果進行比較，結果如圖10所示。

圖9 均質化等效模型中不同區域的位置Fig.9 The locations of different areas in the homogenized equivalent model

圖10 不同區域von Mises應力結果的對比Fig.10 Comparison of results for von Mises stress in different areas

由圖10可知，DFE-submodel方法與全模型結果吻合度較高，相對誤差在1.97%-3.97%之間，比局部均勻化方法的精度高。

3.3 DFE-submodel方法在封裝結構中的應用

下面對簡化后的封裝結構進行仿真分析，仿真模型如圖11所示。該分析過程分為兩步，首先采用DFE方法對該結構的中間部分（10*10 RVE）進行均質化處理，然后對均質化模型進行有限元模擬。載荷為從25℃到150℃的升溫過程，RVE結構為Cu-underfill。有限元網格及邊界條件如圖12所示。

圖11 幾何模型Fig.11 The geometric model

圖12 有限元網格模型和邊界條件Fig.12 The mesh model and boundary condition

采用DFE-submodel方法對該模型模擬后，讀取Cu-underfill結構中Von Mises應力最大位置的仿真結果如圖13(a)所示，全模型結果在圖13(b)中也進行了展示，以驗證該方法的有效性。從圖中可以看出，DFE-submodel方法誤差較小。然后，對幾種方法的網格數量、所需內存、運行時間和不同方向上的應力應變都進行了對比。結果表明，本文所提出的方法能夠用于分析多尺度封裝結構的仿真。

圖13 有限元模擬結果Fig.13 The simulation results of FEM

4 結論

本文基于DFE方法得到了多尺度封裝結構的非線性等效參數，提出了一種能夠用于封裝結構仿真分析的DFE-submodel方法。對Cu-underfill多尺度結構的仿真分析表明，所建立的非線性等效應力應變關系與全模型的結果一致；采用提出的DFE-submodel方法的仿真分析結果與全模型結果吻合度較高。