基于RVE-子模型法的多尺度封裝結構分析方法

孫國立 公顏鵬 侯傳濤 李堯 秦飛

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

0 引言

隨著集成電路發展進入后摩爾定律時代，集成電路產業已經從傳統的不斷減小結構尺寸的方法轉向通過先進封裝技術提升系統性能[1,2]。近年來，系統級封裝、異質集成、3D封裝等先進封裝技術得到越來越多的關注[3,4]。封裝后的芯片以封裝模塊的形式服役，目前，載荷環境愈發復雜[5]。當封裝模塊處在惡劣的工作環境時會導致模塊出現可靠性問題[6]，尤其對于彈載的電子產品，其可靠性對飛行安全至關重要[7]。因此對完整的封裝結構進行仿真建模分析尤為重要。但是，完整的封裝結構中會存在微米尺寸的RDL和TSV-Cu，同時也存在毫米尺寸的基板，甚至會在界面處形成了相互鑲嵌的微觀組織[8]。因此，封裝模型具有典型的結構多尺度特征[9]。雖然計算機硬件和軟件水平已有長足發展，但結構多尺度會導致計算規模迅速增加，耗費大量計算資源。同時計算過程中各類誤差不斷累積，不僅影響計算精度，計算的收斂性也面臨挑戰。因此，亟待發展一種多尺度分析方法，實現幾何多尺度結構的快速精確分析。

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

本文基于代表性體積單元法和子模型技術，將兩種方法耦合，發展了RVE-子模型法。建立多尺度封裝結構的等效模型，在保證計算精度的情況下，達到減小計算規模和降低計算資源的目的。基于該方法本文對三點彎仿真建立等效模型，通過對關鍵位置各應力、應變分量的對比，驗證了本文方法的有效性。

1 RVE-子模型法理論

多尺度結構建模的中心思想是尺度分離。對于周期性結構，利用直接平均理論計算單個RVE的本構，并將其賦給所建立的均質等效模型。再通過子模型技術，將重點關注的位置建立精細模型（子模型），得到關鍵區域的應力應變場。該方法可以對具有結構多尺度的非線性周期性結構進行分析。下面我們以TSV-Cu結構為例介紹本文算法。

1.1 RVE分析法

圖1為建立等效模型的示意圖，在實際的封裝結構中，TSV-Cu直徑在30-100 μm之間，在芯片或轉接板上呈周期性分布，宏觀結構可由單個胞體陣列得到。根據其結構特點，在某一區域內定義一個單胞作為RVE，并將RVE的力學參數作為整個周期性結構的力學參數。

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

RVE的邊界條件對獲得準確本構具有重要影響，研究表明施加周期性邊界條件使得單個RVE變形協調，相鄰RVE界面之間不會分離，可以保證相鄰RVE界面處的應力連續[18]。式（1）為任意RVE周期性邊界條件。

在ABAQUS中實現對RVE劃分周期性網格，并采用多點約束的方式施加周期性邊界條件。式（2）為施加的約束方程，使得相對面上節點在某一方向的位移差為一常數。計算中，需要對RVE每個面上相對節點之間自由度建立約束，由于節點數量過于龐大，本文通過 Python腳本實現線性多點約束。

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

對RVE進行單向拉伸，并利用公式（3）和（4）所示的直接平均方法得到如圖4所示的RVE的應力應變曲線，即等效模型的本構。

式中N為RVE模型的單元總數，Vm是RVE中第m個單元的體積，σm是RVE中第m個單元的應力，mε是RVE中第m個單元的應變。

1.2 子模型技術

子模型方法又稱特定邊界位移法，特定邊界就是子模型從全模型中切割出的邊界。在同一坐標系中對部分區域進行精細建模，并將其作為子模型，其邊界條件即為全模型邊界的位移值。子模型技術基于圣維南原理。下簡要說明子模型的有限元理論。式（5）為基于彈性力學最小位能原理的有限元表達式[19]。

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

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

將式（6）展開得到公式（7）

由（7）式可知，已知的位移向量[]1D轉變為載荷項，可求得未知的[D2]。

本文將等效模型得到的位移場作為子模型邊界條件，利用節點位移和單元應力，得到關注位置的應力應變場。

2 有限元求解

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

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

3 算例驗證

3.1 三點彎模型

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

圖3 三點彎模型Fig.3 Three-point bend model

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

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

3.2 等效模型的有效性驗證

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

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

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

圖6 關注區域各位移分量對比(ux、uy、uz、u)Fig.6 Comparisons of displacement (ux、uy、uz、u)

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

表2 不同模型的計算時間Table 2 Computation time for different models

圖7 關注區域各應力分量對比Fig.7 Comparison of stress components in the area of interest

圖8 關注區域各應變分量對比Fig.8 Comparison of strain components in the area of interest

結果證明本文提出的方法用于研究多尺度封裝結構的可行性。該方法可以降低計算規模，提高計算效率。

4 結論

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