馬沖 公顏鵬 侯傳濤 秦飛

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

0 引言

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

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

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

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

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

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

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

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

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

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

式中，{Fc0}是界面位移為零時(shí)的節(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向虛位移，可得到

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

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值算例

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

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

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

圖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ù)是在硅片上通過(guò)干法刻蝕技術(shù)在硅基體中刻蝕孔洞，然后在其中電鍍填充銅來(lái)實(shí)現(xiàn)晶體管的機(jī)械支撐和電氣互連等目的。后續(xù)需要通過(guò)CMP技術(shù)去除沉積在硅表面的銅層，并實(shí)現(xiàn)晶圓背面TSV結(jié)構(gòu)的暴露及結(jié)構(gòu)正、背面的平坦化[16]。

在此算例中，對(duì)簡(jiǎn)化的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

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

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

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

圖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)的幾何特征、材料屬性等，將整個(gè)模型分為BE域和FE域。自編 BEM 代碼的功能通過(guò) ABAQUS的子程序 UEL 來(lái)實(shí)現(xiàn)，以此獲得用戶自定義單元的有效剛度和有效節(jié)點(diǎn)力。

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