趙勝軍 公顏鵬 侯傳濤 秦飛

（1電子封裝技術(shù)與可靠性研究所，北京工業(yè)大學(xué)，北京 100124；2 可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京強(qiáng)度環(huán)境研究所，北京 100076）

0 引言

電子封裝不斷向輕質(zhì)量，多功能，高性能和高集成度方向發(fā)展[1]，其可靠性問題[2-4]一直是研究的熱點(diǎn)。封裝結(jié)構(gòu)的特征尺寸一般會(huì)相差多個(gè)數(shù)量級(jí)[5]，表現(xiàn)出明顯的結(jié)構(gòu)多尺度特征，例如，晶圓級(jí)封裝中的TSV結(jié)構(gòu)[6]，RDL層[7]等。結(jié)構(gòu)多尺度特征給封裝結(jié)構(gòu)仿真分析帶來了極大的困難，這也對(duì)當(dāng)前封裝結(jié)構(gòu)多尺度仿真分析方法提出了新的需求和挑戰(zhàn)。

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

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

Tan[14]提出了DFE法，它是一種非線性多尺度方法，能夠分析材料的非線性行為。該方法主要針對(duì)二維周期性分布的結(jié)構(gòu)，無法直接應(yīng)用于復(fù)雜的實(shí)際封裝結(jié)構(gòu)多尺度結(jié)構(gòu)的仿真。因此，本文將對(duì)DFE方法進(jìn)行改進(jìn)，以適用于復(fù)雜的封裝結(jié)構(gòu)多尺度結(jié)構(gòu)仿真分析。

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

1 DFE理論及其在Abaqus中的實(shí)現(xiàn)

1.1 DFE理論

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

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

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

有限元分析中，通常采用高斯積分法進(jìn)行數(shù)值積分。因此，

其中，α是單元e的高斯積分點(diǎn)。Jα是雅可比行列式；ωα是高斯點(diǎn)的權(quán)重。

在DFE中，每個(gè)高斯點(diǎn)的應(yīng)力是由相應(yīng)RVE的體積平均應(yīng)力計(jì)算得到的。因此方程(2)可以寫成：

其中，〈·〉α 表示與單元e內(nèi)的高斯點(diǎn)α相關(guān)的RVE上的體積平均量?！啊庇糜诒硎疚⒂^尺度計(jì)算中的量。

Hill-Mandel均質(zhì)化條件要求，

聯(lián)立公式(3)和(4)得

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

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

經(jīng)有限元離散化，式（6）在有限元中表示為，

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

1.2 DFE在Abaqus中的實(shí)現(xiàn)

DFE方法在Abaqus中的實(shí)現(xiàn)可以通過兩個(gè)關(guān)鍵步驟完成。

1)通過Abaqus中的多點(diǎn)約束（multipoint constraints，MPC），建立微觀RVE邊界網(wǎng)格節(jié)點(diǎn)d和宏觀網(wǎng)格節(jié)點(diǎn)之間的關(guān)系。實(shí)現(xiàn)了從宏觀尺度到微觀尺度的過渡。

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

RVE邊界AB和邊界CD節(jié)點(diǎn)之間的約束方程可以寫成：

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

其中，0x指在宏觀單元內(nèi)的高斯點(diǎn)，也是RVE中心的位置。

通過Abaqus中的MPC將約束方程施加到模型上，從而實(shí)現(xiàn)d和之間的聯(lián)系。

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

縮放系數(shù)可以直接確定。例如，對(duì)于二維DFE分析，使用2×2高斯積分點(diǎn)的矩形單元進(jìn)行有限元分析時(shí)

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

圖2 不同縮放系數(shù)所對(duì)應(yīng)的有限元模型Fig. 2 Finite element model for differentα

2 數(shù)值模擬

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

2.1 DFE模型仿真分析

算例中Cu-underfill多尺度結(jié)構(gòu)全模型有3200個(gè)RVE（40*80），RVE網(wǎng)格模型如圖3所示。RVE的長(zhǎng)度和寬度為0.10.1mmmm× ；Cu材料在RVE的中間區(qū)域，Cu的中心與RVE中心重合，Cu的直徑為0.05mm；為了保證仿真結(jié)果的準(zhǔn)確性，對(duì)全模型進(jìn)行了網(wǎng)格收斂性驗(yàn)證，當(dāng)RVE的網(wǎng)格單元數(shù)量為457時(shí)，仿真結(jié)果趨于穩(wěn)定，因此在后續(xù)的仿真中均采用該網(wǎng)格密度。

圖3 RVE網(wǎng)格模型Fig. 3 RVE mesh model

2.1.1 DFE仿真模型的建立

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

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

2.1.2 材料參數(shù)

在宏觀尺度上，整個(gè)異質(zhì)材料被離散化為均質(zhì)的連續(xù)有限單元，宏觀尺度模型的計(jì)算不需要均勻的本構(gòu)關(guān)系，這是因?yàn)檎w剛度矩陣完全是由RVE的剛度矩陣乘以縮放系數(shù)得到。RVE 中Cu和underfill材料參數(shù)如表1所示，underfill為線彈性材料，Cu為非線性材料，屈服應(yīng)力225MPa，切線模型6666MPa。由于宏觀尺度上不需要均勻的本構(gòu)關(guān)系，因此，宏觀模型材料參數(shù)只需要設(shè)置較小的彈性模量即可。

表1 材料屬性Table 1 Material properties

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

DFE模型中左邊界完全固定約束，在右邊界所有節(jié)點(diǎn)上施加單向拉伸載荷， 2xmm= 。A區(qū)域的von Mises 應(yīng)力和等效塑性應(yīng)變分布云圖如圖5所示，von Mises應(yīng)力值為1001MPa，等效塑性應(yīng)變?yōu)?.09699。仿真分析需要的運(yùn)行時(shí)間為4.8min，所需內(nèi)存576MB。將該模型仿真結(jié)果（S11、S22、S12、E11、E22和E12）與全模型相應(yīng)結(jié)果進(jìn)行了對(duì)比，如圖6和圖7所示。由圖可知，兩個(gè)模型的仿真結(jié)果較為吻合，說明了DFE方法是有效的。

圖5 DFE模型A區(qū)域仿真結(jié)果Fig.5 DFE model simulation results in area A

圖6 全模型與DFE模型不同應(yīng)力結(jié)果對(duì)比Fig.6 The results of different stress for DFE and Full model

圖7 全模型與DFE模型不同應(yīng)變結(jié)果對(duì)比Fig.7 The results of different strain for DFE and Full model

2.2 不同多尺度方法仿真結(jié)果的對(duì)比

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

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

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

b）局部均勻化方法所得結(jié)果誤差最大，且網(wǎng)格數(shù)量仍然較多，所需內(nèi)存也較多，在三種方法中效果最不理想。

c）DFE方法是預(yù)測(cè)精度最高，所需網(wǎng)格數(shù)量和運(yùn)行時(shí)間最少的方法，且內(nèi)存占比較少，在三種方法中效果最為理想。

3 DFE-submodel多尺度方法

DFE方法主要針對(duì)二維周期性分布的結(jié)構(gòu)，不能夠直接應(yīng)用于復(fù)雜的實(shí)際封裝結(jié)構(gòu)多尺度的仿真。因此，基于DFE方法和子模型技術(shù)，發(fā)展了一種DFE與子模型結(jié)合（DFE-submodel）的多尺度方法。該方法首先利用DFE方法得到等效非線性參數(shù)，然后，對(duì)等效非線性均質(zhì)化全局模型進(jìn)行仿真分析，最后結(jié)合子模型技術(shù)得到模型中關(guān)鍵區(qū)域的應(yīng)力應(yīng)變值。

3.1 DFE法非線性等效參數(shù)的計(jì)算

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

由圖8可知，DFE模型與全模型結(jié)果吻合度較高，能夠較好的得到Cu-underfill等效應(yīng)力應(yīng)變關(guān)系，不同縮放系數(shù)的應(yīng)力應(yīng)變結(jié)果差別較小。根據(jù)應(yīng)力應(yīng)變關(guān)系可以得到等效的彈性模量為3.84GPa，將得到的彈塑性參數(shù)作為后續(xù)均質(zhì)化模型的材料參數(shù)。

圖8 Cu-underfill等效應(yīng)力應(yīng)變曲線Fig.8 Equivalent stress-strain relationship for Cu-underfill

3.2 DFE-submodel仿真分析

均質(zhì)化等效模型仿真結(jié)束后，選取均質(zhì)化模型中關(guān)鍵區(qū)域進(jìn)行子模型仿真。子模型結(jié)構(gòu)中包含RVE模型，從而可以實(shí)現(xiàn)對(duì)模型中微觀區(qū)域仿真結(jié)果的提取。為了驗(yàn)證該方法的有效性，采用DFE-submodel方法對(duì)Cu-underfill結(jié)構(gòu)（2.1節(jié)中的全模型）進(jìn)行仿真分析。對(duì)比了18個(gè)不同區(qū)域的最大von Mises應(yīng)力結(jié)果，不同區(qū)域的位置如圖9所示。并且與全模型和局部均勻化模型的結(jié)果進(jìn)行比較，結(jié)果如圖10所示。

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

圖10 不同區(qū)域von Mises應(yīng)力結(jié)果的對(duì)比Fig.10 Comparison of results for von Mises stress in different areas

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

3.3 DFE-submodel方法在封裝結(jié)構(gòu)中的應(yīng)用

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

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

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

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

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

4 結(jié)論

本文基于DFE方法得到了多尺度封裝結(jié)構(gòu)的非線性等效參數(shù)，提出了一種能夠用于封裝結(jié)構(gòu)仿真分析的DFE-submodel方法。對(duì)Cu-underfill多尺度結(jié)構(gòu)的仿真分析表明，所建立的非線性等效應(yīng)力應(yīng)變關(guān)系與全模型的結(jié)果一致；采用提出的DFE-submodel方法的仿真分析結(jié)果與全模型結(jié)果吻合度較高。