電熱多物理耦合模擬的軟件設計模式研究*

龔 平，鄭宇騰，張愛清

北京應用物理與計算數(shù)學研究所，北京 100083

1 引言

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

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

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

2 電熱耦合理論方法

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

2.1 電場和溫度場方程模板

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

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

對于時諧電磁場，可由頻域麥克斯韋方程導出矢量電場波動方程。

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

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

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

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

2.2 電-熱耦合模型

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

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

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

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

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

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

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

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

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

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

3 數(shù)理方程構件化設計模式

電熱耦合模擬軟件的快速組裝要求仿真組件滿足可復用、可互操作、可定制三個特性?？蓮陀弥傅氖牵恍薷奈锢磉^程計算構件就可以直接在計算流程中調用；可互操作指的是，物理場數(shù)據(jù)在不同物理過程計算構件間的正確傳遞；可定制指的是，不修改物理過程計算構件，軟件開發(fā)者可以修改仿真模擬的流程和耦合過程。

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

3.1 可復用：封裝單物理場

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

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

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

（2）讀取計算必須的統(tǒng)一數(shù)據(jù)結構的物理場數(shù)據(jù)（輸入）；

（3）提供統(tǒng)一數(shù)據(jù)結構的物理場數(shù)據(jù)（輸出）。

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

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

Table 1 Equation component表1 方程構件

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

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

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

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

3.3 可定制：選擇方程構件和耦合構件、搭建計算流程

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

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

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

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

Fig.3 Computation flow圖3 計算流程

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

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

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

4.1 方程構件

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

以矢量電場波動方程構件為例，圖4為方程構件的基本結構。矢量電場波動方程構件包括的功能有讀取電學參數(shù)，并計算電場E。

4.2 耦合變量的管理

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

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

4.3 定制耦合構件與組合實現(xiàn)計算流程

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

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

Fig.4 Equation component structure圖4 方程構件結構

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

5 實驗驗證

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

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

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

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

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

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

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

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

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

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

5.3 可復用的快速組裝驗證

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

可復用軟件模塊包括上文介紹的電流連續(xù)性方程構件、矢量電場方程構件、穩(wěn)態(tài)熱平衡方程構件，以及數(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)耦合分析算例及溫度場分布對比

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

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

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

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)耦合分析算例及溫度場分布對比

6 結束語

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

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ù)