孫國(guó)立 公顏鵬 侯傳濤 李堯 秦飛

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

0 引言

隨著集成電路發(fā)展進(jìn)入后摩爾定律時(shí)代，集成電路產(chǎn)業(yè)已經(jīng)從傳統(tǒng)的不斷減小結(jié)構(gòu)尺寸的方法轉(zhuǎn)向通過(guò)先進(jìn)封裝技術(shù)提升系統(tǒng)性能[1,2]。近年來(lái)，系統(tǒng)級(jí)封裝、異質(zhì)集成、3D封裝等先進(jìn)封裝技術(shù)得到越來(lái)越多的關(guān)注[3,4]。封裝后的芯片以封裝模塊的形式服役，目前，載荷環(huán)境愈發(fā)復(fù)雜[5]。當(dāng)封裝模塊處在惡劣的工作環(huán)境時(shí)會(huì)導(dǎo)致模塊出現(xiàn)可靠性問(wèn)題[6]，尤其對(duì)于彈載的電子產(chǎn)品，其可靠性對(duì)飛行安全至關(guān)重要[7]。因此對(duì)完整的封裝結(jié)構(gòu)進(jìn)行仿真建模分析尤為重要。但是，完整的封裝結(jié)構(gòu)中會(huì)存在微米尺寸的RDL和TSV-Cu，同時(shí)也存在毫米尺寸的基板，甚至?xí)诮缑嫣幮纬闪讼嗷ヨ偳兜奈⒂^組織[8]。因此，封裝模型具有典型的結(jié)構(gòu)多尺度特征[9]。雖然計(jì)算機(jī)硬件和軟件水平已有長(zhǎng)足發(fā)展，但結(jié)構(gòu)多尺度會(huì)導(dǎo)致計(jì)算規(guī)模迅速增加，耗費(fèi)大量計(jì)算資源。同時(shí)計(jì)算過(guò)程中各類(lèi)誤差不斷累積，不僅影響計(jì)算精度，計(jì)算的收斂性也面臨挑戰(zhàn)。因此，亟待發(fā)展一種多尺度分析方法，實(shí)現(xiàn)幾何多尺度結(jié)構(gòu)的快速精確分析。

目前，針對(duì)多尺度封裝結(jié)構(gòu)的分析，最簡(jiǎn)單常用的方法是體積等效法，該方法可以很好地估計(jì)兩相復(fù)合材料的等效熱膨脹系數(shù)[10-12]和等效彈性模量[10]。根據(jù)該均勻化方法，Lee研究發(fā)現(xiàn)當(dāng)銅通孔占總體積14%以下時(shí)，等效楊氏模量接近于硅，而當(dāng)銅的含量為88%時(shí)彈性模量接近于銅[13]。Cheng 和 Shen使用有限元計(jì)算的熱應(yīng)變來(lái)估計(jì)等效的熱膨脹系數(shù)[11]。這種方法只能粗略的計(jì)算平面內(nèi)的等效參數(shù)。Che等用同樣的方法根據(jù)TSV-Cu在晶圓上的分布，得到銅/硅胞體的線彈性力學(xué)參數(shù)[14]。但以上方法并沒(méi)有考慮TSV-Cu塑性和胞體的邊界條件的影響，建立的等效模型與真實(shí)模型仍存在一定誤差。一些學(xué)者通過(guò)理論推導(dǎo)的方式得到解析的等效參數(shù)。Chen等人對(duì) TSV轉(zhuǎn)接板的平面內(nèi)和平面外等效機(jī)械性能的理論進(jìn)行了詳細(xì)闡述并驗(yàn)證其有效性[15]。在Chen基礎(chǔ)上，Ye等人通過(guò)在計(jì)算銅和燒結(jié)銀的等效力學(xué)參數(shù)時(shí)引入幾何參數(shù)，對(duì)其計(jì)算理論進(jìn)行了修正，數(shù)值結(jié)果表明全網(wǎng)格模型與等效模型應(yīng)力分布一致[16]。但是，由于封裝結(jié)構(gòu)的復(fù)雜度越來(lái)越高，這種方法的效率較低。因此，亟需發(fā)展一種封裝結(jié)構(gòu)多尺度分析方法。由于三維編織材料的細(xì)觀結(jié)構(gòu)呈周期性排布，因此利用代表性體積單元（Representative Volume Element，RVE）研究三維編織材料是一種有效且常用的方法[17]。與編織材料類(lèi)似，幾何周期性也同樣普遍存在于封裝結(jié)構(gòu)中。

本文基于代表性體積單元法和子模型技術(shù)，將兩種方法耦合，發(fā)展了RVE-子模型法。建立多尺度封裝結(jié)構(gòu)的等效模型，在保證計(jì)算精度的情況下，達(dá)到減小計(jì)算規(guī)模和降低計(jì)算資源的目的?；谠摲椒ū疚膶?duì)三點(diǎn)彎仿真建立等效模型，通過(guò)對(duì)關(guān)鍵位置各應(yīng)力、應(yīng)變分量的對(duì)比，驗(yàn)證了本文方法的有效性。

1 RVE-子模型法理論

多尺度結(jié)構(gòu)建模的中心思想是尺度分離。對(duì)于周期性結(jié)構(gòu)，利用直接平均理論計(jì)算單個(gè)RVE的本構(gòu)，并將其賦給所建立的均質(zhì)等效模型。再通過(guò)子模型技術(shù)，將重點(diǎn)關(guān)注的位置建立精細(xì)模型（子模型），得到關(guān)鍵區(qū)域的應(yīng)力應(yīng)變場(chǎng)。該方法可以對(duì)具有結(jié)構(gòu)多尺度的非線性周期性結(jié)構(gòu)進(jìn)行分析。下面我們以TSV-Cu結(jié)構(gòu)為例介紹本文算法。

1.1 RVE分析法

圖1為建立等效模型的示意圖，在實(shí)際的封裝結(jié)構(gòu)中，TSV-Cu直徑在30-100 μm之間，在芯片或轉(zhuǎn)接板上呈周期性分布，宏觀結(jié)構(gòu)可由單個(gè)胞體陣列得到。根據(jù)其結(jié)構(gòu)特點(diǎn)，在某一區(qū)域內(nèi)定義一個(gè)單胞作為RVE，并將RVE的力學(xué)參數(shù)作為整個(gè)周期性結(jié)構(gòu)的力學(xué)參數(shù)。

圖1 建立等效模型示意圖Fig.1 Schematic diagram of the established equivalent model

RVE的邊界條件對(duì)獲得準(zhǔn)確本構(gòu)具有重要影響，研究表明施加周期性邊界條件使得單個(gè)RVE變形協(xié)調(diào)，相鄰RVE界面之間不會(huì)分離，可以保證相鄰RVE界面處的應(yīng)力連續(xù)[18]。式（1）為任意RVE周期性邊界條件。

在ABAQUS中實(shí)現(xiàn)對(duì)RVE劃分周期性網(wǎng)格，并采用多點(diǎn)約束的方式施加周期性邊界條件。式（2）為施加的約束方程，使得相對(duì)面上節(jié)點(diǎn)在某一方向的位移差為一常數(shù)。計(jì)算中，需要對(duì)RVE每個(gè)面上相對(duì)節(jié)點(diǎn)之間自由度建立約束，由于節(jié)點(diǎn)數(shù)量過(guò)于龐大，本文通過(guò) Python腳本實(shí)現(xiàn)線性多點(diǎn)約束。

式中，Cn為根據(jù)節(jié)點(diǎn)和參考點(diǎn)的位置確定的系數(shù)，n為線性多點(diǎn)約束方程的項(xiàng)數(shù)（本文取3），是位移變量，i表示方向（=1,2,3i），P表示節(jié)點(diǎn)的位置。

對(duì)RVE進(jìn)行單向拉伸，并利用公式（3）和（4）所示的直接平均方法得到如圖4所示的RVE的應(yīng)力應(yīng)變曲線，即等效模型的本構(gòu)。

式中N為RVE模型的單元總數(shù)，Vm是RVE中第m個(gè)單元的體積，σm是RVE中第m個(gè)單元的應(yīng)力，mε是RVE中第m個(gè)單元的應(yīng)變。

1.2 子模型技術(shù)

子模型方法又稱(chēng)特定邊界位移法，特定邊界就是子模型從全模型中切割出的邊界。在同一坐標(biāo)系中對(duì)部分區(qū)域進(jìn)行精細(xì)建模，并將其作為子模型，其邊界條件即為全模型邊界的位移值。子模型技術(shù)基于圣維南原理。下簡(jiǎn)要說(shuō)明子模型的有限元理論。式（5）為基于彈性力學(xué)最小位能原理的有限元表達(dá)式[19]。

式中，[K]為整體剛度矩陣，[D]為待求位移向量，[F]為施加的載荷。

由于子模型的邊界條件是整體模型切割位移值，因此部分位移解已知。假設(shè)位移向量[D]已知的部分為[D1]，待求解的為[D2]。則式（5）可以寫(xiě)為

將式（6）展開(kāi)得到公式（7）

由（7）式可知，已知的位移向量[]1D轉(zhuǎn)變?yōu)檩d荷項(xiàng)，可求得未知的[D2]。

本文將等效模型得到的位移場(chǎng)作為子模型邊界條件，利用節(jié)點(diǎn)位移和單元應(yīng)力，得到關(guān)注位置的應(yīng)力應(yīng)變場(chǎng)。

2 有限元求解

本文以ABAQUS軟件為實(shí)驗(yàn)平臺(tái)，對(duì)RVE進(jìn)行單軸拉伸，采用直接均勻理論得到每一個(gè)子步下平均應(yīng)力和平均應(yīng)變值[20]。圖2利用RVE-子模型法建立計(jì)算宏觀位移場(chǎng)和關(guān)注位置的應(yīng)力應(yīng)變相應(yīng)的流程圖。求解過(guò)程通過(guò)商業(yè)軟件ABAQUS和Python腳本實(shí)現(xiàn)。假設(shè)RVE為各向同性材料，且由RVE在x向和y向陣列得到的非線性復(fù)合板可以等效為均質(zhì)板。首先，對(duì)RVE施加周期性邊界條件并進(jìn)行單軸拉伸得到其應(yīng)力應(yīng)變曲線。并將獲得的RVE本構(gòu)賦給建立的等效模型，求解得到等效模型的位移場(chǎng)，然后，與全網(wǎng)格模型的結(jié)果進(jìn)行對(duì)比，若誤差較小則將等效模型的位移作為子模型的邊界條件，可得到關(guān)注位置的應(yīng)力應(yīng)變場(chǎng)。若誤差較大，則重新選取RVE后再次進(jìn)行求解。

圖2 RVE-子模型法計(jì)算流程Fig.2 RVE-Submodel method calculation process

3 算例驗(yàn)證

3.1 三點(diǎn)彎模型

如圖3所示，考慮一種非線性周期性結(jié)構(gòu)的三點(diǎn)彎模型，模型由15個(gè)單胞組成，每個(gè)單胞的尺寸均為111×× mm，中間支撐圓柱的半徑為0.08 mm，整個(gè)模型的最大和最小尺寸相差兩個(gè)數(shù)量級(jí)，具備一定的多尺度特征。為體現(xiàn)RVE-子模型在處理非線性問(wèn)題上的優(yōu)勢(shì)，設(shè)置復(fù)合板的基體為T(mén)SV-Cu，其本構(gòu)如圖4所示[21]，增強(qiáng)體為硅，彈性模量為131 GPa，泊松比為0.3[22]。由于基體材料為非線性本構(gòu)，因此本文假設(shè)等效模型的也為非線性結(jié)構(gòu)，通過(guò)第2節(jié)所述方法得到等效模型的非線性本構(gòu)。有限元模型的兩個(gè)支撐圓柱和中間壓頭圓柱均被約束為剛體，兩端的支撐為固定約束，中間壓頭施加向下的0.3 mm位移。三個(gè)圓柱壓頭與中間的復(fù)合板之間采用綁定約束。全網(wǎng)格模型和基于RVE-子模型法建立的等效模型參數(shù)對(duì)比表1所示。為排除計(jì)算機(jī)硬件的影響，采用同一臺(tái)計(jì)算機(jī)和相同的單元類(lèi)型（C3D8R）對(duì)兩種模型進(jìn)行劃分網(wǎng)格，得到的全網(wǎng)格模型的單元數(shù)為835165，而等效模型的單元數(shù)為371925，等效模型單元約為全網(wǎng)格模型的一半。本算例假設(shè)關(guān)注區(qū)域?yàn)閳D3所示的RVE區(qū)域，因此，需要對(duì)該位置建立精細(xì)模型并將等效模型的計(jì)算結(jié)果作為子模型的邊界條件。

圖3 三點(diǎn)彎模型Fig.3 Three-point bend model

圖4 模擬單軸拉伸下的應(yīng)力應(yīng)變曲線Fig.4 Stress-strain curve under uniaxial tensile

表1 有限元模型參數(shù)Table 1 Parameters used in the finite element model

3.2 等效模型的有效性驗(yàn)證

為驗(yàn)證本文RVE-子模型法所建立的等效模型的有效性。首先對(duì)比了相同載荷下兩種模型的宏觀位移場(chǎng)，其次對(duì)比了關(guān)注位置的全網(wǎng)格模型和子模型的計(jì)算結(jié)果，最后比較了兩種模型在相同計(jì)算機(jī)配置下的計(jì)算時(shí)間。圖5為模擬三點(diǎn)彎曲實(shí)驗(yàn)下實(shí)際模型和等效均質(zhì)模型位移云圖。由圖可知兩種模型得到的最大位移相差僅7.6 μm，并且兩種方法計(jì)算結(jié)果的趨勢(shì)一致。

圖5 全網(wǎng)格模型與等效模型的位移云圖對(duì)比Fig.5 Comparisons of displacement obtained by the full mesh model and present method

圖6 給出了本文RVE-子模型法得到的關(guān)注區(qū)域位移值與全網(wǎng)格模型得到的值對(duì)比情況，ux、uy、分別為關(guān)注區(qū)域x、y、z三個(gè)方向的最大位移值和最大合位移。從圖中可以看出兩種方法得到位移值基本一致，最大相差0.3 mm，結(jié)果吻合較好。

圖6 關(guān)注區(qū)域各位移分量對(duì)比(ux、uy、uz、u)Fig.6 Comparisons of displacement (ux、uy、uz、u)

圖7 為關(guān)注位置的應(yīng)力分量對(duì)比，由圖可知本文算法得到的最大值與參考值趨勢(shì)相同。其中子模型得到的關(guān)注位置最大von Mises應(yīng)力為1464 MPa，參考值為1114 MPa，兩值相差較小。因此，可為材料失效分析提供依據(jù)。x方向正應(yīng)力模擬值為-786.3 MPa，與參考值的誤差僅為5.02%。許多經(jīng)典疲勞理論是基于應(yīng)變理論的（如Manson-Coffin等），因此應(yīng)變對(duì)于結(jié)構(gòu)可靠性有重要影響。關(guān)注位置的應(yīng)變分量對(duì)比在圖8給出，本文數(shù)值解和參考值在同一水平，與圖7對(duì)比可知，應(yīng)力大的分量應(yīng)變也相應(yīng)大，這與力學(xué)理論一致。兩種模型的計(jì)算時(shí)間如表2所示。與全模型相比，等效模型在計(jì)算結(jié)果精度一致的情況下，計(jì)算時(shí)間為40.85 min，比全模型減少了78.57%。等效模型和子模型的計(jì)算時(shí)間之和也僅有全模型的五分之一，這進(jìn)一步體現(xiàn)了RVE-子模型法處理結(jié)構(gòu)多尺度的優(yōu)勢(shì)。

表2 不同模型的計(jì)算時(shí)間Table 2 Computation time for different models

圖7 關(guān)注區(qū)域各應(yīng)力分量對(duì)比Fig.7 Comparison of stress components in the area of interest

圖8 關(guān)注區(qū)域各應(yīng)變分量對(duì)比Fig.8 Comparison of strain components in the area of interest

結(jié)果證明本文提出的方法用于研究多尺度封裝結(jié)構(gòu)的可行性。該方法可以降低計(jì)算規(guī)模，提高計(jì)算效率。

4 結(jié)論

針對(duì)封裝結(jié)構(gòu)中存在大量多尺度結(jié)構(gòu)的問(wèn)題，本文將代表性體積單元法和子模型技術(shù)相結(jié)合發(fā)展了RVE-子模型法。利用ABAQUS可以解決封裝結(jié)構(gòu)中大量周期結(jié)構(gòu)多尺度建模的困難，實(shí)現(xiàn)宏觀和微觀的雙尺度交互。文章給出了該方法的具體計(jì)算過(guò)程，并通過(guò)三點(diǎn)彎仿真驗(yàn)證了方法的有效性。結(jié)果表明， RVE-子模型法節(jié)省了大量計(jì)算時(shí)間，加快了計(jì)算進(jìn)程，為解決先進(jìn)封裝結(jié)構(gòu)多尺度建模問(wèn)題提供一種新方法。