用于分析電子封裝結(jié)構(gòu)的有限元-邊界元耦合方法研究

馬沖 公顏鵬 侯傳濤 秦飛

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

0 引言

芯片封裝是指利用膜技術(shù)及微細(xì)加工技術(shù)，將芯片及其他要素在框架或基板上布置、粘貼固定及連接，引出接線端子并通過可塑性絕緣介質(zhì)灌封固定，構(gòu)成整體立體結(jié)構(gòu)的工藝[1]。封裝結(jié)構(gòu)涉及多種不同性能的材料，而且在界面處存在幾何不連續(xù)性。在服役過程中，溫度場、濕度場等惡劣環(huán)境會使封裝結(jié)構(gòu)產(chǎn)生多種失效問題，例如，空洞、開裂、翹曲和分層等。隨著封裝密度增加和功能的多樣化，具有多物理場、多尺度特征的電子封裝結(jié)構(gòu)的可靠性問題正成為設(shè)計的關(guān)鍵。

由于易實現(xiàn)、精度高、無環(huán)境限制等特點(diǎn)，數(shù)值模擬技術(shù)正成為電子封裝結(jié)構(gòu)設(shè)計領(lǐng)域的主流方案。有限元法（FEM）是一種常用的獲得各種工程問題近似解的數(shù)值分析技術(shù)。目前，商業(yè)有限元軟件ABAQUS、ANSYS等被廣泛用于電子封裝結(jié)構(gòu)的可靠性分析[2-7]。

實際上，由于電子封裝結(jié)構(gòu)的多尺度特征，有限元模型往往需要被離散成大量的單元，以保證數(shù)值結(jié)果的合理性。特別是對于一些具有多尺度特征的復(fù)雜模型，會耗費(fèi)大量計算時間和成本。邊界元法（BEM）是用以求解工程和科學(xué)問題的常用數(shù)值分析方法之一，其特點(diǎn)是降低求解問題的自由度數(shù)和只需要邊界離散化。有限元法和邊界元法各有特點(diǎn)，因此發(fā)揮各自特長的有限元-邊界元耦合方法一直受到研究者的重視。Liu等人[8]提出了一種用于解決動態(tài)彈塑性問題的有限元-邊界元耦合方法。Fernandes 等人[9]提出了一種分析非均勻材料板彎曲問題的有限元-邊界元耦合方法。Bonari 等人[10]提出了一種用于求解接觸問題的有限元-邊界元耦合方法，其中有限元用于處理宏觀模型、邊界元用于處理微觀粗糙模型。Dong 等人[11]利用邊界元法對力載荷作用下的集成電路封裝進(jìn)行應(yīng)力計算。基于模型幾何和加載的對稱性質(zhì)，得到相應(yīng)的基本解。Neto 等人[12]提出了一種將等幾何邊界元法（IGABEM）與有限元軟件相結(jié)合的方法。此研究工作為彈性力學(xué)中的三維IGABEM公式提出了網(wǎng)格自適應(yīng)方法，提供了精確的幾何表示和力學(xué)場描述，有助于BEM和CAD方案的耦合。Qin等人[13]提出一種新穎的有限元-邊界元耦合方法來研究多尺度結(jié)構(gòu)二維穩(wěn)態(tài)傳熱問題，其中邊界元用于處理不重要或線性區(qū)域。Gong等人[14]提出了一種能夠研究電子封裝中多尺度結(jié)構(gòu)傳熱問題的等幾何邊界元方法。數(shù)值結(jié)果與解析解和有限元解的比較驗證了此方法的有效性。雖然邊界元法可以減少單元的數(shù)量，但對于非線性、非均勻問題的分析，該方法還存在很多問題。

結(jié)合邊界元法和有限元法各自的優(yōu)勢，本文提出了一種用于分析電子封裝結(jié)構(gòu)的有限元-邊界元耦合方法，將自編寫的BEM代碼集成到商業(yè)軟件ABAQUS中。在耦合方法中，將整個模型域劃分為有限元（FE）域和邊界元（BE）域。有限元法（借助ABAQUS）用于分析具有非線性或非均勻特性的子域，而具有線性特征或不重要的域則由BEM代碼求解。邊界元域被視為有限元法的一個超單元，它的有效剛度和有效節(jié)點(diǎn)力可通過BEM代碼求解得到。通過子程序UEL將超單元的等效量裝配到整體模型的有限元公式。通過這種方法，實現(xiàn)用ABAQUS進(jìn)行電子封裝結(jié)構(gòu)的耦合方法分析。

1 有限元-邊界元耦合理論

有限元法與邊界元法的耦合方法需要把整個求解域分成兩部分，一部分以有限元法進(jìn)行離散，另一部分以邊界元法進(jìn)行離散。求解方法基本上分為兩類：一類是先求解邊界元域，把邊界元區(qū)域看作是有限元法的一個超單元，通過邊界條件和有限元方程耦合起來求解；另一類是先求解有限元域，然后將有限元解作為邊界元域的邊界條件，再用邊界元法求解[15]。本文采用第一類求解方法。

對于先求解邊界元域的方法，邊界元域的邊界面力應(yīng)轉(zhuǎn)化為有限元中的等效節(jié)點(diǎn)力。如圖1所示的模型，由一個BE域和一個FE域組成。對于BE域，其界面的位移{uc}和面力{tc}之間的關(guān)系可以表示為[15]

圖1 邊界元域與有限元域的界面力Fig.1 Interface forces between FE domain and BE domain

式中，{tc0}是所有界面位移為零時的面力向量，[KBE]是邊界元域的偽剛度矩陣。

對于FE域，界面位移{uc}和界面節(jié)點(diǎn)力{Fc}之間的關(guān)系為

式中，{Fc0}是界面位移為零時的節(jié)點(diǎn)力向量，[KFE]是有限元域界面節(jié)點(diǎn)的剛度矩陣。

在界面的某一節(jié)點(diǎn)上施加x方向的任意虛位移，則有限元域界面節(jié)點(diǎn)力做的功為

式中，n為單元e的節(jié)點(diǎn)數(shù)。根據(jù)守恒條件，節(jié)點(diǎn)力做的功應(yīng)與面力做的功相等

式中，Γ為積分域邊界（如圖1所示）。將面力tx和虛位移xuδ表示為插值形式

其中，Ni和Nj為界面處的形函數(shù)。

由公式(3)-(6)，可以得到

同理，施加y向虛位移，可得到

因此，對公式(1)左乘矩陣[M]，得到等效節(jié)點(diǎn)力的表達(dá)式

其中，由[Me]組合而成的矩陣[M]表示整個界面的變換矩陣，矩陣[Me]可以通過下式獲得

由公式(9)可知，可以通過轉(zhuǎn)換矩陣[M]建立有限元法與邊界元法的關(guān)系。將公式（9）帶入公式（2），即可形成完整的耦合方程。

2 耦合算法的實現(xiàn)

本文所開發(fā)的有限元法-邊界元法耦合算法是通過子程序UEL把邊界元程序集成到有限元軟件ABAQUS中。圖2展示了將自編寫的邊界元程序與ABAQUS耦合起來的過程：

圖2 耦合算法實現(xiàn)流程圖Fig.2 The flowchart of coupling scheme

1）根據(jù)模型的幾何信息、材料特性等參數(shù)將模型劃分為有限元子域和邊界元子域。其中有限元子域需通過ABAQUS /CAE處理，而邊界元子域通過邊界元程序進(jìn)行計算。

2）為了使邊界元程序能與ABAQUS軟件耦合在一起，需在ABAQUS/CAE中以標(biāo)準(zhǔn)操作建立兩個模型：有限元域模型和等效邊界元模型。等效邊界元部分是將邊界元域等效為與界面相重合的幾何模型， 即接觸界面。需要注意的是應(yīng)該保證等效邊界元部分所劃分的節(jié)點(diǎn)與有限元域界面部分的節(jié)點(diǎn)保持一致，然后兩個模型通過綁定約束“ Tie ”組裝起來。

3）ABAQUS/CAE完成兩個部件的前處理后，輸出存儲模型信息的初始input文件。

4）對初始input文件進(jìn)行修改，將等效邊界元部分修改為用戶自定義單元，重新定義其幾何信息、物理屬性和邊界條件等參數(shù)。

5）將修改后的input文件提交計算，同時調(diào)用包含自編寫邊界元程序的子程序UEL來實現(xiàn)邊界元程序與ABAQUS的耦合。

6）通過UEL得到用戶自定義單元（即邊界元域界面處）的等效剛度矩陣和等效節(jié)點(diǎn)載荷，并傳遞給ABAQUS，此時有限元域便可通過ABAQUS求解器進(jìn)行求解得到有限元域的所有物理量。

3 數(shù)值算例

本節(jié)通過兩個算例來實現(xiàn)該算法對常見結(jié)構(gòu)的分析并探索其在實際問題中的應(yīng)用，進(jìn)一步探究該算法的性能。

3.1 具有復(fù)雜界面的矩形板結(jié)構(gòu)

在這個算例中，我們選擇帶有復(fù)雜界面形狀的矩形板結(jié)構(gòu)，該模型由兩個部分組成，其模型示意圖如圖3所示。該模型的幾何參數(shù)為a= 1 mm，l1= 0.2 mm，l2= 1 mm。邊界條件為懸臂梁左側(cè)固定，右側(cè)有沿x方向分布的均布載荷q= 1000 N/mm。該模型兩個域的彈性模量為E= 1000 MPa，泊松比為v= 0.2。耦合方法對模型求解時，左側(cè)區(qū)域用有限元法離散，右側(cè)區(qū)域用邊界元法離散。為了驗證該耦合算法的精度，在此選用商業(yè)軟件ABAQUS的細(xì)化網(wǎng)格模型的計算結(jié)果作為參考解。耦合算法處理的模型被離散成1個自定義單元和800個四節(jié)點(diǎn)單元，而細(xì)化的有限元模型被離散成5200個四節(jié)點(diǎn)單元。針對該矩形板模型，從提交作業(yè)到分析完成的過程，耦合方法用時24 s，有限元法用時50 s。圖4 (a) 和 (b) 所示為該耦合方法和有限元法得到的左側(cè)部分的von Mises應(yīng)力云圖。從圖中可以看出該耦合算法對于彈性問題的求解有著較好的精度。

圖3 帶有復(fù)雜界面的多尺度結(jié)構(gòu)示意圖Fig.3 Multiscale structure with complex interfaces

圖4 兩種方法得到的左側(cè)結(jié)構(gòu)的von Mises應(yīng)力分布云圖Fig.4 von Mises stress contours computed by the present method and FEM

3.2 電子封裝TSV-Cu結(jié)構(gòu)

此模型為三維電子封裝中較為先進(jìn)的硅通孔（Through silicon via, TSV）結(jié)構(gòu)。TSV-Cu 技術(shù)是在硅片上通過干法刻蝕技術(shù)在硅基體中刻蝕孔洞，然后在其中電鍍填充銅來實現(xiàn)晶體管的機(jī)械支撐和電氣互連等目的。后續(xù)需要通過CMP技術(shù)去除沉積在硅表面的銅層，并實現(xiàn)晶圓背面TSV結(jié)構(gòu)的暴露及結(jié)構(gòu)正、背面的平坦化[16]。

在此算例中，對簡化的TSV結(jié)構(gòu)模型（如圖5所示）進(jìn)行力學(xué)分析。其中TSV-Cu 部分的直徑為10 μm，深度為100 μm，電鍍Cu淀積層的厚度為δ= 6μm；硅部分的幾何尺寸分別為a= 630 μm，b= 250 μm。硅部分的底部被固定，在淀積層的上表面施加剪切載荷t= 0.25 N/m。TSV-Cu部分的彈性模量ECu= 155 GPa，泊松比v= 0.3，屈服強(qiáng)度為47.91MPa；Si部分的彈性模量ESi= 131 GPa，泊松比v= 0.25。

圖5 TSV模型的示意圖Fig.5 Cross-section of TSV

在耦合算法計算過程中，TSV-Cu部分以有限元法劃分網(wǎng)格，而硅部分以邊界元法劃分邊界 網(wǎng)格，如圖6所示的網(wǎng)格示意圖，在TSV - Cu區(qū)域使用四節(jié)點(diǎn)四邊形單元，而在硅部分使用二次邊界單元進(jìn)行離散。耦合算法處理的模型被離散成1個自定義單元和7434個四節(jié)點(diǎn)單元，而細(xì)化的有限元模型被離散成161934個四節(jié)點(diǎn)單元。由于硅部分只需在邊界部分進(jìn)行離散，所以該耦合方法對于這種跨尺度結(jié)構(gòu)可以大幅度減少單元數(shù)量。對于硅通孔模型，從提交作業(yè)到分析完成的過程，耦合算法用時82 s，有限元法用時102 s。

圖6 TSV結(jié)構(gòu)的網(wǎng)格示意圖Fig.6 Meshes used in the TSV model

圖7展示了該耦合算法和有限元法計算得到的位于路徑L上的節(jié)點(diǎn)處的von Mises應(yīng)力結(jié)果，有限元法的計算結(jié)果作為參考結(jié)果。從圖中可以清晰的看出該耦合算法對于非線性問題的求解也有著較好的精度。

圖7 路徑L上各點(diǎn)的von Mises應(yīng)力分布值Fig.7 The von Mises stress at points along the contour L

4 結(jié)論

本文提出了一種將商業(yè)軟件ABAQUS和自編BEM代碼結(jié)合的耦合方法。在當(dāng)前的耦合方法中，根據(jù)結(jié)構(gòu)的幾何特征、材料屬性等，將整個模型分為BE域和FE域。自編 BEM 代碼的功能通過 ABAQUS的子程序 UEL 來實現(xiàn)，以此獲得用戶自定義單元的有效剛度和有效節(jié)點(diǎn)力。

數(shù)值算例說明了該算法在求解彈性問題和非線性問題時均有著較高的精度，能大幅降低單元數(shù)量并減少計算時間。