龔 平，鄭宇騰，張愛清

北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100083

1 引言

電熱多物理場(chǎng)耦合數(shù)值模擬是大規(guī)模集成電路設(shè)計(jì)的重要手段。特別是隨著設(shè)計(jì)及制造工藝朝著三維系統(tǒng)級(jí)封裝、異質(zhì)異構(gòu)集成的方向演進(jìn)，導(dǎo)致功率密度顯著增加，熱與應(yīng)力可靠性問(wèn)題更加突出[1]，并且與電磁問(wèn)題相互耦合、相互影響，成為影響與制約電路性能與可靠性的核心因素，因此電、磁、熱、力等多物理場(chǎng)耦合分析是大規(guī)模集成電路設(shè)計(jì)的一個(gè)關(guān)鍵環(huán)節(jié)[1-4]。而大規(guī)模集成電路結(jié)構(gòu)精細(xì)，工況多樣，物理機(jī)理復(fù)雜，涉及電-熱、電-熱-力等各種不同的多物理場(chǎng)耦合，以及靜態(tài)、瞬態(tài)、時(shí)諧等多種分析類型[2-4]。面向這種多樣化的多物理場(chǎng)耦合分析需求，如何快速研制批量電熱多物理場(chǎng)耦合模擬軟件，是一個(gè)挑戰(zhàn)性的問(wèn)題。

基于算子分裂方法，復(fù)用單一物理場(chǎng)求解代碼，實(shí)現(xiàn)多物理場(chǎng)耦合計(jì)算，是快速開發(fā)多物理耦合軟件的有效途徑。但以往國(guó)內(nèi)外科研學(xué)者通常專注于物理模型和數(shù)值算法的研究，不關(guān)注代碼復(fù)用，一般會(huì)將物理建模細(xì)節(jié)，甚至工程建模細(xì)節(jié)，硬編碼在數(shù)值求解子程序中，導(dǎo)致幾乎沒(méi)有代碼重用的可能性，耦合軟件都只能從頭開發(fā)。例如Jin等人[5]對(duì)直流饋電網(wǎng)絡(luò)的電壓降及熱問(wèn)題的研究，Shao等人[6]對(duì)高功率集成電路的電-熱耦合仿真的研究，是眾多從頭開發(fā)電熱多物理場(chǎng)耦合模擬軟件的兩個(gè)普通例子。因此為滿足多樣化的多物理場(chǎng)耦合模擬需求，需要一種高效率、高可復(fù)用、能夠快速組裝計(jì)算流程的多物理耦合軟件開發(fā)模式，實(shí)現(xiàn)單一物理場(chǎng)求解代碼一次開發(fā)，多次利用。

針對(duì)電熱多物理場(chǎng)耦合模擬需求，本文提出了支持快速組裝、高可定制的并行模擬軟件設(shè)計(jì)模式。該模式的主要思路是基于可復(fù)用的電場(chǎng)和磁場(chǎng)方程模板，開發(fā)可復(fù)用的方程求解構(gòu)件庫(kù)，然后應(yīng)用基于方程求解構(gòu)件組裝耦合計(jì)算流程，從而達(dá)到快速研發(fā)多物理場(chǎng)耦合計(jì)算軟件的目的。

2 電熱耦合理論方法

芯片工作時(shí)，電磁場(chǎng)損耗產(chǎn)生焦耳熱，焦耳熱作為溫度場(chǎng)的熱源，導(dǎo)致芯片溫度上升；而溫度升高則又會(huì)改變芯片電學(xué)參數(shù)，進(jìn)而影響電磁場(chǎng)分布。因此在芯片的數(shù)值模擬中，至少需要考慮電場(chǎng)和溫度場(chǎng)間的耦合效應(yīng)[9]，涉及的電磁數(shù)值求解包括靜電問(wèn)題和時(shí)諧電磁問(wèn)題，分別對(duì)應(yīng)電流連續(xù)性方程、矢量電場(chǎng)波動(dòng)方程，熱數(shù)值模擬為穩(wěn)態(tài)熱平衡，對(duì)應(yīng)的方程為穩(wěn)態(tài)熱平衡方程，物理場(chǎng)數(shù)據(jù)耦合需要考慮電磁場(chǎng)產(chǎn)生的熱源引起的溫度場(chǎng)變化，同時(shí)根據(jù)溫度場(chǎng)的分布來(lái)更新電學(xué)參數(shù)。

2.1 電場(chǎng)和溫度場(chǎng)方程模板

由電流連續(xù)性方程可以導(dǎo)出靜電問(wèn)題所要求解的方程[10]：

其中，σ(r)為材料電導(dǎo)率；U(r)為電位；Ub(r)為第一類邊界上的電壓值。對(duì)于電位未知量，使用結(jié)點(diǎn)有限元單元離散。

對(duì)于時(shí)諧電磁場(chǎng)，可由頻域麥克斯韋方程導(dǎo)出矢量電場(chǎng)波動(dòng)方程。

其中，μ(r)為材料磁導(dǎo)率；ε(r)為材料介電常數(shù)；ω為角頻率；E(r)為電場(chǎng)；J(r)為電流源分布。對(duì)于電場(chǎng)未知量，使用棱邊有限元單元離散。

穩(wěn)態(tài)熱平衡方程是：

其中，κ(r,T)、c、ρ分別為材料的熱導(dǎo)率、熱容、密度；T(r)為溫度分布；T為第一類邊界上的溫度值；h為第三類邊界上的熱對(duì)流系數(shù)；Tambient為環(huán)境溫度；f(r)為熱源項(xiàng)。對(duì)于溫度場(chǎng)未知量，使用結(jié)點(diǎn)有限元單元離散。

在電場(chǎng)、電位、溫度場(chǎng)未知量離散后，使用伽遼金方法建立有限元矩陣方程[11-12]，并通過(guò)基于Krylov子空間的迭代方法[13]求解，最終得到場(chǎng)分布。

2.2 電-熱耦合模型

電-熱耦合分析涉及正向和反向兩種耦合模式，分別是電磁場(chǎng)產(chǎn)生熱源的正向耦合模式和溫度場(chǎng)改變電學(xué)參數(shù)的反向耦合[14]。

對(duì)于正向耦合，包括靜電場(chǎng)產(chǎn)生的熱源，計(jì)算公式為：

f(T,t)=σ|?U(t)|2

以及時(shí)諧電磁場(chǎng)產(chǎn)生的熱源，計(jì)算公式為：

f(T,ω)=σ|E(ω)|2

對(duì)于反向耦合，溫度將改變電學(xué)參數(shù)。上述方程（1）、（2）中的電學(xué)參數(shù)將演變?yōu)闇囟群涂臻g的函數(shù)，其形式如下：

σ(r)=σ(T(r),r)

μ(r)=μ(T(r),r)

ε(r)=ε(T(r),r)

本文的耦合處理方案是松耦合求解。每次耦合迭代中，依次求解電磁場(chǎng)方程，計(jì)算電磁場(chǎng)產(chǎn)生的熱源，求解熱平衡方程，根據(jù)溫度場(chǎng)更新電學(xué)參數(shù)。通過(guò)溫度場(chǎng)分布計(jì)算殘差，判斷是否終止迭代。

3 數(shù)理方程構(gòu)件化設(shè)計(jì)模式

電熱耦合模擬軟件的快速組裝要求仿真組件滿足可復(fù)用、可互操作、可定制三個(gè)特性?？蓮?fù)用指的是，不修改物理過(guò)程計(jì)算構(gòu)件就可以直接在計(jì)算流程中調(diào)用；可互操作指的是，物理場(chǎng)數(shù)據(jù)在不同物理過(guò)程計(jì)算構(gòu)件間的正確傳遞；可定制指的是，不修改物理過(guò)程計(jì)算構(gòu)件，軟件開發(fā)者可以修改仿真模擬的流程和耦合過(guò)程。

而根據(jù)第2章所展示電-熱場(chǎng)的控制方程與物理場(chǎng)的數(shù)據(jù)耦合，單物理場(chǎng)求解是多物理場(chǎng)耦合分析的共性，對(duì)應(yīng)的模擬個(gè)性則是根據(jù)計(jì)算流程定制的數(shù)據(jù)耦合處理方法。因此，使用統(tǒng)一數(shù)據(jù)結(jié)構(gòu)的并行編程框架，開發(fā)共性的單物理場(chǎng)計(jì)算構(gòu)件及個(gè)性定制的物理場(chǎng)間耦合器，便可快速開發(fā)高效的耦合軟件。

3.1 可復(fù)用：封裝單物理場(chǎng)

基于統(tǒng)一的數(shù)據(jù)結(jié)構(gòu)封裝共性，即將單物理求解流程封裝為求解單控制方程的計(jì)算構(gòu)件（簡(jiǎn)稱方程構(gòu)件），在仿真模擬軟件開發(fā)過(guò)程中實(shí)現(xiàn)復(fù)用。

如圖1所示，封裝后的方程構(gòu)件具有以下功能：

（1）配置控制參數(shù)及求解操作；

（2）讀取計(jì)算必須的統(tǒng)一數(shù)據(jù)結(jié)構(gòu)的物理場(chǎng)數(shù)據(jù)（輸入）；

（3）提供統(tǒng)一數(shù)據(jù)結(jié)構(gòu)的物理場(chǎng)數(shù)據(jù)（輸出）。

Fig.1 Coupling components in electric-thermal coupling simulation圖1 電-熱耦合模擬中的耦合構(gòu)件

對(duì)于本文求解的電磁-熱耦合問(wèn)題，封裝的方程構(gòu)件分別是電流連續(xù)性方程構(gòu)件、矢量電場(chǎng)方程構(gòu)件、穩(wěn)態(tài)熱平衡方程構(gòu)件，方程構(gòu)件的功能和輸入輸出如表1所示。

Table 1 Equation component表1 方程構(gòu)件

3.2 可互操作：物理場(chǎng)間數(shù)據(jù)交互

在電-熱耦合模擬過(guò)程中，電磁場(chǎng)在導(dǎo)體或介質(zhì)中損耗，產(chǎn)生焦耳熱，焦耳熱改變溫度場(chǎng)；溫度場(chǎng)的改變又導(dǎo)致材料的電學(xué)參數(shù)發(fā)生變化。電磁場(chǎng)對(duì)溫度場(chǎng)的作用通過(guò)改變焦耳熱實(shí)現(xiàn)，溫度場(chǎng)對(duì)電磁場(chǎng)的作用通過(guò)改變電學(xué)參數(shù)實(shí)現(xiàn)，因此焦耳熱和電學(xué)參數(shù)被稱為耦合變量；對(duì)于僅存在于一個(gè)物理過(guò)程中的變量（如電位）稱為非耦合變量。物理場(chǎng)間的數(shù)據(jù)交互指的就是耦合變量的交互。

Fig.2 Equation component in electro-thermal coupling simulation圖2 電-熱耦合模擬中的方程構(gòu)件

對(duì)于非耦合變量，由于其不需要在物理場(chǎng)間交互，便由所在方程構(gòu)件單獨(dú)管理；相對(duì)照的，耦合變量會(huì)被多個(gè)物理過(guò)程計(jì)算模塊訪問(wèn)和修改，因此需要使用統(tǒng)一的數(shù)據(jù)耦合構(gòu)件創(chuàng)建和管理，在其他的方程構(gòu)件中需通過(guò)數(shù)據(jù)耦合構(gòu)件，創(chuàng)建耦合變量和訪問(wèn)耦合變量。

3.3 可定制：選擇方程構(gòu)件和耦合構(gòu)件、搭建計(jì)算流程

在3.1 節(jié)中抽象和封裝單物理場(chǎng)控制方程的共性，形成可復(fù)用的方程構(gòu)件，而對(duì)于仿真模擬中的個(gè)性部分，包括耦合變量的創(chuàng)建與管理部分，則會(huì)被封裝成數(shù)據(jù)耦合構(gòu)件；根據(jù)特定計(jì)算場(chǎng)景或模型分析的需求，通過(guò)組合方程構(gòu)件與數(shù)據(jù)耦合構(gòu)件，定制計(jì)算流程。

在電-熱耦合模擬中，物理場(chǎng)數(shù)據(jù)耦合部分包括計(jì)算焦耳熱、更新材料的電學(xué)參數(shù)，如圖1。對(duì)于不同的仿真模擬應(yīng)用，可以通過(guò)定制數(shù)據(jù)耦合構(gòu)件，滿足不同的耦合需求。

以矩形導(dǎo)體棒中的電-熱耦合問(wèn)題為例，涉及的控制方程有電流連續(xù)性方程與穩(wěn)態(tài)熱平衡方程，如圖2 所示。電流連續(xù)性方程構(gòu)件從數(shù)據(jù)耦合構(gòu)件中讀取電學(xué)參量作為輸入，計(jì)算電場(chǎng)分布，并作為輸出量傳遞給數(shù)值耦合構(gòu)件；穩(wěn)態(tài)熱平衡方程，方程構(gòu)件從數(shù)據(jù)耦合構(gòu)件中讀取熱源作為輸入，計(jì)算溫度場(chǎng)分布，并作為輸出量傳遞給數(shù)值耦合構(gòu)件。

最終個(gè)性化軟件定制時(shí)，如圖3 所示，只需要在主控函數(shù)中組合使用方程構(gòu)件與數(shù)據(jù)耦合構(gòu)件。首先調(diào)用電磁方程構(gòu)件求解電流連續(xù)性方程或矢量電場(chǎng)方程；然后調(diào)用物理場(chǎng)數(shù)據(jù)耦合構(gòu)件更新焦耳熱，穩(wěn)態(tài)熱平衡方程構(gòu)件求解溫度場(chǎng)；物理場(chǎng)數(shù)據(jù)耦合構(gòu)件更新電學(xué)參數(shù)，最后判斷計(jì)算是否結(jié)束，如果未結(jié)束開始下一輪循環(huán)。

Fig.3 Computation flow圖3 計(jì)算流程

4 基于JAUMIN的并行實(shí)現(xiàn)

并行自適應(yīng)非結(jié)構(gòu)網(wǎng)格應(yīng)用支撐軟件框架（J parallel adaptive unstructured mesh applications infra-structure，JAUMIN），提供超大規(guī)模并行計(jì)算和網(wǎng)格自適應(yīng)加密等功能[15]。

本章將介紹基于JAUMIN 框架如何實(shí)現(xiàn)電熱耦合并行仿真模擬的快速開發(fā)。

4.1 方程構(gòu)件

根據(jù)JAUMIN 框架對(duì)軟件結(jié)構(gòu)的要求，物理過(guò)程計(jì)算構(gòu)件包括網(wǎng)格片層次結(jié)構(gòu)時(shí)間積分算法、網(wǎng)格層時(shí)間積分算法、網(wǎng)格片計(jì)算策略類。在網(wǎng)格層時(shí)間積分算法、網(wǎng)格片計(jì)算策略類中配置單個(gè)物理場(chǎng)的參數(shù)和計(jì)算，通過(guò)網(wǎng)格片層次結(jié)構(gòu)時(shí)間積分算法調(diào)用網(wǎng)格層時(shí)間積分算法、網(wǎng)格片計(jì)算策略類，實(shí)現(xiàn)自動(dòng)并行。

以矢量電場(chǎng)波動(dòng)方程構(gòu)件為例，圖4為方程構(gòu)件的基本結(jié)構(gòu)。矢量電場(chǎng)波動(dòng)方程構(gòu)件包括的功能有讀取電學(xué)參數(shù)，并計(jì)算電場(chǎng)E。

4.2 耦合變量的管理

在JAUMIN 框架中，使用變量管理器統(tǒng)一管理物理場(chǎng)數(shù)據(jù)，通過(guò)變量ID索引和訪問(wèn)。因此方程構(gòu)件的數(shù)據(jù)輸入輸出可以通過(guò)傳遞變量ID實(shí)現(xiàn)分享數(shù)據(jù)的目的。

耦合變量的創(chuàng)建和管理由數(shù)據(jù)耦合構(gòu)件負(fù)責(zé)，電流連續(xù)性方程構(gòu)件、矢量電場(chǎng)方程構(gòu)件和穩(wěn)態(tài)熱平衡方程構(gòu)件通過(guò)接口可獲得目標(biāo)數(shù)據(jù)變量ID。方程構(gòu)件與耦合構(gòu)件間交互數(shù)據(jù)變量的偽代碼如下所示。

4.3 定制耦合構(gòu)件與組合實(shí)現(xiàn)計(jì)算流程

依照J(rèn)AUMIN 框架要求，耦合構(gòu)件接口和結(jié)構(gòu)設(shè)計(jì)如圖5所示。與方程構(gòu)件相同，耦合構(gòu)件包括網(wǎng)格片層次結(jié)構(gòu)時(shí)間積分算法、網(wǎng)格層時(shí)間積分算法、網(wǎng)格片計(jì)算策略類，主要的功能接口包括獲取耦合變量ID、計(jì)算焦耳熱和更新電學(xué)參數(shù)。

計(jì)算流程定制在main函數(shù)通過(guò)組合方程構(gòu)件和耦合構(gòu)件實(shí)現(xiàn)，以靜電-熱穩(wěn)態(tài)耦合分析為例，main函數(shù)偽代碼形式如下所示。

Fig.4 Equation component structure圖4 方程構(gòu)件結(jié)構(gòu)

Fig.5 Data coupling component structure圖5 數(shù)據(jù)耦合構(gòu)件結(jié)構(gòu)

5 實(shí)驗(yàn)驗(yàn)證

根據(jù)本文提出的快速組裝的并行電熱耦合模擬軟件設(shè)計(jì)模式，編寫靜電分析、頻域分析以及熱穩(wěn)態(tài)分析模塊，在計(jì)算流程中根據(jù)需要組合使用電-熱模塊和可定制的耦合模塊搭建計(jì)算流程，就可以快速完成多物理耦合模擬軟件開發(fā)。

為驗(yàn)證本文提出的軟件模式可靠性、復(fù)用性，對(duì)于電-熱耦合問(wèn)題，開發(fā)了以下兩款軟件。

5.1 靜電-熱穩(wěn)態(tài)模擬軟件

仿真程序的流程圖如圖6 所示。驗(yàn)證模型為矩形截面的導(dǎo)體棒，模型尺寸長(zhǎng)、寬、高為50 mm、25 mm、10 mm，材料為銅，波導(dǎo)兩端設(shè)為恒溫邊界，溫度為295 K。

在仿真計(jì)算過(guò)程中調(diào)用靜電分析以及熱穩(wěn)態(tài)分析，與原有的程序溫度場(chǎng)分布對(duì)比如圖7，其中圖7（a）是定制程序仿真結(jié)果，圖7（b）是原有程序仿真結(jié)果。從圖中可以看出，溫度場(chǎng)數(shù)值是一致的，定制程序的仿真結(jié)果正確。

5.2 頻域-熱穩(wěn)態(tài)模擬軟件

驗(yàn)證模型為濾波器[16]，如圖8 所示，內(nèi)部材料包括空氣、硅和銅，問(wèn)題的求解頻率為30 GHz，目標(biāo)整體尺寸長(zhǎng)、寬、高分別為6.88 mm、4.12 mm、1.37 mm，濾波器上下表面與空氣的對(duì)流系數(shù)為30 W/(m2·K)。

Fig.6 Flow graph for electrostatic and steady-state thermal coupling analysis圖6 靜電-熱穩(wěn)態(tài)耦合分析流程圖

仿真程序的流程圖如圖9所示。

在仿真計(jì)算過(guò)程中調(diào)用頻域分析以及熱穩(wěn)態(tài)分析，與原有的程序溫度場(chǎng)分布對(duì)比如圖10，其中圖10（a）是定制程序仿真結(jié)果，圖10（b）是原有程序仿真結(jié)果，從圖中可以看出溫度場(chǎng)數(shù)值是一致的，定制程序的仿真結(jié)果正確。

5.3 可復(fù)用的快速組裝驗(yàn)證

本文通過(guò)對(duì)比代碼復(fù)用率[17]，也就是程序中可復(fù)用軟件模塊與不可復(fù)用軟件模塊代碼行數(shù)，證實(shí)了本文提出的軟件設(shè)計(jì)模式實(shí)現(xiàn)代碼的快速?gòu)?fù)用。

可復(fù)用軟件模塊包括上文介紹的電流連續(xù)性方程構(gòu)件、矢量電場(chǎng)方程構(gòu)件、穩(wěn)態(tài)熱平衡方程構(gòu)件，以及數(shù)值離散模塊、材料管理器。靜電-熱穩(wěn)態(tài)分析程序、頻域-熱穩(wěn)態(tài)分析程序中，各模塊的代碼行數(shù)（去掉空行和注釋后）如表2、表3所示。

Fig.7 Example of electrostatic,steady-state thermal coupling analysis and comparison of temperature distribution圖7 靜電-熱穩(wěn)態(tài)耦合分析算例及溫度場(chǎng)分布對(duì)比

Fig.8 Frequency-domain,steady-state thermal coupling analysis model structure圖8 頻域-熱穩(wěn)態(tài)耦合分析仿真模型結(jié)構(gòu)

靜電-熱穩(wěn)態(tài)分析程序中，全部程序代碼行數(shù)為7 219，可復(fù)用代碼行數(shù)占全部代碼行數(shù)為87.6%；頻域-熱穩(wěn)態(tài)分析程序中，全部程序代碼行數(shù)為7 888，可復(fù)用代碼行數(shù)占全部代碼行數(shù)為88.7%。

兩個(gè)程序中，數(shù)值離散模塊的代碼行數(shù)最多，復(fù)用率最高，其次是矢量電場(chǎng)波動(dòng)方程構(gòu)件（表中Electric-FieldFluctuation 構(gòu)件），復(fù)用率最低的是穩(wěn)態(tài)熱平衡方程構(gòu)件（表中SteadyStateHeatConduction 構(gòu)件）；不可復(fù)用代碼占全部代碼比例不足15%，預(yù)計(jì)后續(xù)的其他仿真模擬應(yīng)用開發(fā)的工作量將會(huì)有效減少。

Fig.9 Flow graph for frequency-domain and steady-state thermal coupling analysis圖9 頻域-熱穩(wěn)態(tài)耦合分析流程圖

Fig.10 Example of frequency-domain,steady-state thermal coupling analysis and comparison of temperature distribution圖10 頻域-熱穩(wěn)態(tài)耦合分析算例及溫度場(chǎng)分布對(duì)比

6 結(jié)束語(yǔ)

本文研究多物理仿真模擬中，電-熱耦合模塊的快速組裝。實(shí)驗(yàn)結(jié)果證明該方法在保證計(jì)算結(jié)果正確的基礎(chǔ)上，實(shí)現(xiàn)了電子芯片仿真模擬中的電-熱耦合模擬的快速組裝開發(fā)。本文提出的軟件復(fù)用模式是否適用于更多更復(fù)雜的非線性多物理耦合或多尺度多物理耦合應(yīng)用，還有待進(jìn)一步研究和驗(yàn)證。

Table 2 Number of code lines of coupled electrostatic,steady-state thermal analysis software表2 靜電-熱穩(wěn)態(tài)耦合模擬軟件代碼行數(shù)

Table 3 Number of code lines of coupled frequency-domain,steady-state thermal analysis software表3 頻域-熱穩(wěn)態(tài)耦合模擬軟件代碼行數(shù)