陳智源，周曉明，方 芳，何 敏，劉 浩,3，孫佳寧

（1. 電子科技大學(xué)機(jī)械與電氣工程學(xué)院，四川成都 611731；2. 中國電子科技集團(tuán)公司第十研究所，四川成都 610036；3. 中國電子科技集團(tuán)公司第二十九研究所，四川成都 610036；4. 電子科技大學(xué)廣西智能制造產(chǎn)業(yè)技術(shù)研究院，廣西柳州 545003）

引 言

熱力耦合是指結(jié)構(gòu)中溫度場(chǎng)的變化對(duì)結(jié)構(gòu)的影響以及結(jié)構(gòu)變形在外部載荷作用下對(duì)溫度場(chǎng)的影響，即結(jié)構(gòu)、溫度場(chǎng)通過交互作用而相互影響的物理現(xiàn)象。熱力耦合問題廣泛存在于航空航天[1-3]、機(jī)械[4-6]、半導(dǎo)體微電子[7-9]等領(lǐng)域，特別是隨著現(xiàn)代電子設(shè)備向著高度集成化的方向發(fā)展，電子設(shè)備內(nèi)部的熱力耦合效應(yīng)更為明顯，對(duì)電子設(shè)備進(jìn)行熱力耦合效應(yīng)研究的需求逐漸增加。

通過仿真工具對(duì)多物理場(chǎng)耦合現(xiàn)象進(jìn)行工程計(jì)算，是現(xiàn)代電子設(shè)備研究與設(shè)計(jì)中常用的方法[10]。目前能夠進(jìn)行熱力耦合仿真的多物理場(chǎng)耦合工具有Ansys Multiphysics，COMSOL Multiphysics，ABAQUS，ADINA，MATLAB等。ANSYS提供結(jié)構(gòu)力學(xué)、熱學(xué)等計(jì)算程序，程序之間有給定接口，主要耦合方式為順序耦合、算子分裂和固定點(diǎn)迭代，但核心代碼封裝，改動(dòng)耦合方式和求解算法的難度較大；與ANSYS不同，COMSOL Multiphysics是在同一框架內(nèi)的多場(chǎng)耦合，支持多種模塊，部分計(jì)算能夠?qū)崿F(xiàn)牛頓迭代的全耦合，主要耦合方式是順序耦合和固定點(diǎn)迭代；ABAQUS和ADINA可以全耦合求解熱電耦合；MATLAB是基于代碼的數(shù)學(xué)軟件，通過特定代碼可以實(shí)現(xiàn)各種耦合算法。

上述基于牛頓迭代的耦合計(jì)算工具在面對(duì)多場(chǎng)同時(shí)求解時(shí)可能產(chǎn)生大型稀疏矩陣難以收斂求解的問題，而其他基于多求解器并采用順序耦合、固定點(diǎn)迭代等的傳統(tǒng)耦合策略的仿真軟件在雙向耦合模式下進(jìn)行多場(chǎng)強(qiáng)耦合仿真時(shí)經(jīng)常難以收斂，而且多個(gè)求解器之間存在異構(gòu)網(wǎng)格映射、數(shù)據(jù)傳遞的問題。鑒于上述問題，亟需探索一種適用于電子設(shè)備熱力強(qiáng)耦合仿真的新方法。多物理場(chǎng)耦合仿真（Multiphysics Object-Oriented Simulation Environment，MOOSE）[11]作為新一代開源多物理場(chǎng)仿真框架，采用Jacobian-Free Newton-Krylov (JFNK) 方法[12]求解強(qiáng)耦合問題偏微分（Partial Differential Equation, PDE）方程組。本文開展基于MOOSE框架的熱力全耦合計(jì)算方法研究，有望實(shí)現(xiàn)高性能熱力耦合仿真，具有廣闊的應(yīng)用前景。

1 方法

1.1 基于MOOSE框架的耦合仿真應(yīng)用程序開發(fā)

實(shí)現(xiàn)熱力耦合首先需要在MOOSE框架的基礎(chǔ)上開發(fā)求解程序。如圖1所示，MOOSE框架分為3層。底層由libMesh有限元庫、PETSc等組件構(gòu)成。libMesh主要進(jìn)行網(wǎng)格讀入、有限元離散、殘差矩陣形成和結(jié)果輸出；PETSc求解器則可以求解殘差矩陣。中間層為框架通用接口，包括耦合接口、殘差估算和雅可比矩陣估算。物理層包括內(nèi)核系統(tǒng)、邊界條件和材料系統(tǒng)。工程師可在這3層基礎(chǔ)上結(jié)合需求開發(fā)應(yīng)用程序。

圖1 MOOSE框架結(jié)構(gòu)示意圖

MOOSE熱力耦合仿真應(yīng)用程序的核心是內(nèi)核(Kernel)，目前，MOOSE框架中內(nèi)置了熱學(xué)、化學(xué)、N-S方程、接觸、彈性力學(xué)、時(shí)間處理等多種內(nèi)核。MOOSE框架下的多場(chǎng)耦合基于多個(gè)單場(chǎng)模塊的內(nèi)核組合實(shí)現(xiàn)，因此為實(shí)現(xiàn)熱力耦合，需要構(gòu)建結(jié)構(gòu)力學(xué)和傳熱學(xué)模塊。

結(jié)構(gòu)力學(xué)模塊基于彈性力學(xué)方程構(gòu)建。典型彈性力學(xué)問題的基本方程如下：

平衡方程為

式中：σ為應(yīng)力；F為體積力，即F= [Fx Fy Fz]T；A為微分算子。A的矩陣形式為

幾何方程為

式中：ε為應(yīng)變；u為位移。對(duì)于小應(yīng)變的線性彈性問題，采用線性化小應(yīng)變理論，有

式中，?為那勃勒算子。

應(yīng)力應(yīng)變本構(gòu)方程為

式中，D為彈性矩陣，它完全取決于材料的彈性模量E和泊松比ν。

將幾何方程代入本構(gòu)方程，再代入平衡方程得：

傳熱學(xué)模塊構(gòu)建基于熱傳導(dǎo)功能，有瞬態(tài)熱傳導(dǎo)方程如下：

式中:cp為比熱容；ρ為密度；T為溫度；t為時(shí)間；k為導(dǎo)熱系數(shù)；x為空間坐標(biāo)系下的矢量；˙q為熱源。

瞬態(tài)熱方程弱形式為：

式中：φi為有限元試函數(shù)；Th為有限元結(jié)果；?n為邊界向外的法向量。

將上述理論模型寫成代碼，構(gòu)成了熱傳導(dǎo)內(nèi)核。熱力耦合中使用的是ADHeatCondution類內(nèi)核，此類內(nèi)核是熱傳導(dǎo)中熱擴(kuò)散方程在自動(dòng)微分框架下的實(shí)現(xiàn)，為前向微分。ADHeatConduction內(nèi)核實(shí)現(xiàn)了傅里葉定律給出的熱方程，其中熱通量q為：

在MOOSE中，在默認(rèn)情況下，熱內(nèi)核在位移網(wǎng)格上運(yùn)行，其變形由力學(xué)模型計(jì)算。其本質(zhì)是通過將熱傳導(dǎo)計(jì)算嵌入到網(wǎng)格可運(yùn)動(dòng)的力學(xué)計(jì)算過程中來求解熱力問題，熱和力的耦合通過熱膨脹本征應(yīng)變模塊施加。對(duì)于力學(xué)模塊的本征應(yīng)變，需要將力學(xué)模塊的應(yīng)力-應(yīng)變方程替換為下面的熱彈性體的應(yīng)力-應(yīng)變方程：

式中：ε0為溫度引起的應(yīng)變；α為熱膨脹系數(shù)；T0為參考溫度，即低于該溫度材料收縮，高于該溫度材料膨脹。

1.2 應(yīng)用程序校核

熱力強(qiáng)耦合現(xiàn)象中溫度場(chǎng)產(chǎn)生熱變形，熱變形又會(huì)影響傳熱，溫度場(chǎng)、變形場(chǎng)因此而耦合。本文選擇商業(yè)軟件COMSOL Multiphysics中簡單的基準(zhǔn)案例（即鋼材熱變形模型[13-14]）在所開發(fā)的應(yīng)用程序中進(jìn)行計(jì)算，所得結(jié)果與COMSOL Multiphysics結(jié)果進(jìn)行對(duì)比，完成程序的校核。

該模型為一長寬高均為100 mm的鋼塊，其比熱容為453 J/(kg·K)，導(dǎo)熱系數(shù)為45 W/(m·K)，楊氏模量為2.1×1010Pa，泊松比為0.3，密度為7 850 kg/m3。在Trelis Pro 16.5中對(duì)模型分網(wǎng)，采用尺寸為2 mm的正六面體網(wǎng)格，模型及網(wǎng)格如圖2所示。

圖2 用于校核的基準(zhǔn)模型及其網(wǎng)格

在MOOSE中按如下步驟編寫輸入文件：

1）在輸入文檔中寫入網(wǎng)格命令，并將網(wǎng)格命令導(dǎo)入外網(wǎng)格。

2）主變量規(guī)定為鋼塊溫度和三維形變量disp_x，disp_y，disp_z，且均為一階拉格朗日插值函數(shù)。

3）同時(shí)寫入彈性力學(xué)和熱傳導(dǎo)兩內(nèi)核。

4）設(shè)置邊界條件。固體外表面設(shè)置為狄利克雷邊界，在固體一角施加三向固定約束，使材料可以各向形變。

5）設(shè)置材料。設(shè)置鋼材料的楊氏模量、泊松比、導(dǎo)熱系數(shù)和熱膨脹系數(shù)。參考溫度設(shè)為293.15 K。

6）設(shè)置求解器。求解算法為Preconditioning Jacobian-Free Newton-Krylov(PJFNK)方法，非線性相對(duì)容差為1e-14，線性相對(duì)容差為1e-6，默認(rèn)預(yù)處理。

1.3 電子設(shè)備熱力耦合仿真

電子設(shè)備往往由于其中的高功率組件發(fā)熱而變形，這種問題涉及到熱力強(qiáng)耦合。為探索將開發(fā)的應(yīng)用程序應(yīng)用在復(fù)雜電子設(shè)備中的可能性，本文面向某機(jī)載相控陣天線陣面進(jìn)行熱力耦合仿真。

采用的天線對(duì)象來自文獻(xiàn)[15]設(shè)計(jì)的相控陣天線陣面。該天線的總體尺寸為316.8 mm×24 mm×0.404 mm，貼片及微帶線為銅質(zhì)材料，介質(zhì)板材料為Rogers RT/5880。其T/R組件尺寸為60.6 mm×56.4 mm，T/R 組件單通道總功率為45.88 W。天線高功率放大器的功率為44.80 W，極限熱功率為60.65 W，尺寸為6.77 mm×77 mm×2 mm，材料為GaAs。電源尺寸為100 mm×80.74 mm×8 mm，外殼為鋁合金。使用Trelies Pro 16.5工具對(duì)模型分網(wǎng)，采用正六面體網(wǎng)格，單元尺寸為2 mm，生成如圖3所示的網(wǎng)格。

圖3 某相控陣天線冷板及陣面模型

2 結(jié)果與分析

2.1 校核仿真結(jié)果與分析

設(shè)置不同熱邊界條件，在MOOSE框架下運(yùn)行編制的應(yīng)用程序，所得熱力耦合位移結(jié)果如圖4所示。

在COMSOL Multiphysics中調(diào)取上述基準(zhǔn)案例結(jié)果，對(duì)比MOOSE應(yīng)用程序仿真結(jié)果與COMSOL Multiphysics結(jié)果，見表1。

圖4 熱力耦合仿真校核案例仿真結(jié)果

表1 熱力耦合校核案例仿真結(jié)果

結(jié)果表明，MOOSE框架下的位移極值與商業(yè)軟件COMSOL Multiphysics基準(zhǔn)案例提供的位移數(shù)值的最大相對(duì)誤差小于0.1%，說明熱力耦合內(nèi)核能夠進(jìn)行正確、準(zhǔn)確的計(jì)算。

2.2 天線仿真結(jié)果與分析

在MOOSE框架下運(yùn)行所編制的應(yīng)用程序，所得熱力耦合結(jié)果如圖5和圖6所示。

圖5 天線熱力耦合仿真位移云圖

圖6 天線熱力耦合仿真溫度云圖

采集熱源中心點(diǎn)的溫度和熱應(yīng)變數(shù)據(jù)，結(jié)果見表2。由表2可知，天線陣面中心位置溫度高，靠近邊緣位置溫度低，同時(shí)天線面板上下兩側(cè)的熱變形較明顯，所得結(jié)果與工程經(jīng)驗(yàn)相符。

表2 熱力耦合結(jié)果溫度及熱位移表

3 結(jié)束語

本文面向電子設(shè)備熱力強(qiáng)耦合仿真需求，基于MOOSE框架內(nèi)置的結(jié)構(gòu)力學(xué)和傳熱學(xué)內(nèi)核，編寫了全耦合仿真程序，并且通過商業(yè)軟件的基準(zhǔn)案例對(duì)仿真程序進(jìn)行了校核，進(jìn)而將仿真程序應(yīng)用于工程案例。通過對(duì)某相控陣天線進(jìn)行熱力耦合仿真，得到了天線陣面溫度場(chǎng)和位移場(chǎng)的數(shù)據(jù)。本文的主要結(jié)論如下：

1）在校驗(yàn)案例仿真中，基于MOOSE框架的熱力耦合仿真應(yīng)用程序能夠?qū)崿F(xiàn)熱力全耦合仿真。通過在彈性力學(xué)內(nèi)核框架中嵌入傳熱學(xué)的熱應(yīng)變方程組，可以實(shí)現(xiàn)網(wǎng)格共享。仿真所得溫度場(chǎng)下的熱變形位移結(jié)果與基準(zhǔn)案例結(jié)果基本一致，表明所開發(fā)的程序可正確計(jì)算熱力強(qiáng)耦合問題。

2）相控陣天線陣面熱變形仿真結(jié)果表明，由于天線介質(zhì)板與冷板表面綁定接觸，熱變形量為0，因此脫離冷板約束后，天線熱變形量呈增加趨勢(shì)，最大變形量出現(xiàn)在中部上邊緣。所述結(jié)果符合預(yù)期，反映基于MOOSE框架的熱力強(qiáng)耦合仿真可應(yīng)用于復(fù)雜電子設(shè)備，有望解決當(dāng)前電子設(shè)備熱力強(qiáng)耦合的迫切需求，具有廣闊的應(yīng)用前景。