謝 鵬 徐 立 尹俊輝 楊中海 李 斌

(電子科技大學(xué)電子科學(xué)與工程學(xué)院 成都 610054)

1 引言

微波管廣泛應(yīng)用于衛(wèi)星通信、雷達(dá)系統(tǒng)、電子對(duì)抗以及科學(xué)研究領(lǐng)域[1–3]，其中永磁聚焦系統(tǒng)常作為外部磁場(chǎng)用于電子束聚焦[4]。目前，在大量的科學(xué)及工程應(yīng)用中，有限元方法針對(duì)復(fù)雜結(jié)構(gòu)仍然是主流的數(shù)值分析工具[5,6]。在永磁聚焦系統(tǒng)仿真設(shè)計(jì)領(lǐng)域，國(guó)外Maxwell[7]是最流行的3維有限元仿真設(shè)計(jì)商業(yè)軟件，國(guó)內(nèi)電子科技大學(xué)開發(fā)的微波管模擬器套裝[8]里面的永磁聚焦模擬器[4]則是具有國(guó)家自主知識(shí)產(chǎn)權(quán)的3維磁場(chǎng)仿真設(shè)計(jì)軟件。它們都擁有永磁聚焦系統(tǒng)仿真能力，但是當(dāng)求解大規(guī)模和多尺度問(wèn)題時(shí)，其有限元矩陣求解通常會(huì)花費(fèi)大量的時(shí)間和內(nèi)存，有時(shí)甚至由于缺少有效的預(yù)處理導(dǎo)致矩陣無(wú)法求解。

非重疊區(qū)域分解法采用“分而治之”的思想，將原求解域劃分成若干個(gè)互不重疊的子區(qū)域進(jìn)行求解，其具有天然的并行性，非常適合仿真大型復(fù)雜結(jié)構(gòu)[9]。本文采用基于磁標(biāo)勢(shì)有限元的非重疊區(qū)域分解方法來(lái)進(jìn)行永磁聚焦系統(tǒng)的仿真計(jì)算。歸根結(jié)底，基于磁標(biāo)勢(shì)的永磁聚焦系統(tǒng)仿真屬于求解泊松方程邊值問(wèn)題。在這類邊值問(wèn)題的基于有限元的區(qū)域分解方法中，目前的處理方式主要有兩種，一種是拉格朗日乘子類型的區(qū)域分解法[10,11]；一種是基于內(nèi)罰方式的區(qū)域分解法[12,13]。前者會(huì)有兩種不同類型的推導(dǎo)過(guò)程：一是保留拉格朗日乘子的方式，這種方式會(huì)額外增加未知量個(gè)數(shù)，并且會(huì)生成一個(gè)對(duì)稱不定的鞍點(diǎn)矩陣系統(tǒng)，不利于方程的求解；二是在推導(dǎo)過(guò)程中消去拉格朗日乘子，這種方式會(huì)生成一個(gè)對(duì)稱正定的系數(shù)矩陣，但是這種方式有可能會(huì)因?yàn)榫薮蟮挠?jì)算代價(jià)而不能顯式地計(jì)算系數(shù)矩陣[14]?；趦?nèi)罰的區(qū)域分解方法則不需要引入諸如拉格朗日乘子類型的輔助變量，只需要將傳輸條件引入基于內(nèi)罰方式的有限元弱形式推導(dǎo)過(guò)程中，目前采用的主要是Robin傳輸條件，由Lions[15]首次提出，但目前基于Robin傳輸條件的內(nèi)罰區(qū)域分解法需要考慮法向偏導(dǎo)項(xiàng)的計(jì)算，并且最終形成的有限元矩陣是非對(duì)稱的。

本文提出的區(qū)域分解方法同樣是基于內(nèi)罰方式的，但是引入了一種新型傳輸條件，其來(lái)源于接觸熱阻的定義。相比于現(xiàn)有的方法，該區(qū)域分解方法的主要優(yōu)勢(shì)包括：(1)不需要引入多余的未知量，使得有限元矩陣維數(shù)更少；(2)有限元矩陣集成過(guò)程更加簡(jiǎn)單，只需要考慮區(qū)域交界面上的物理量，而且不需要進(jìn)行法向偏導(dǎo)項(xiàng)的計(jì)算，更重要的是最終產(chǎn)生的有限元矩陣滿足對(duì)稱正定性，矩陣性質(zhì)更好，非常適合采用共軛梯度法進(jìn)行求解。通過(guò)對(duì)多個(gè)永磁結(jié)構(gòu)的仿真計(jì)算可以發(fā)現(xiàn)，相比于商業(yè)軟件Maxwell，本文所提出的區(qū)域分解方法在保證求解精度的同時(shí)，可以更加高效地實(shí)現(xiàn)對(duì)微波管永磁聚焦系統(tǒng)的仿真。

2 基于有限元的區(qū)域分解方法

2.1 基于內(nèi)罰的區(qū)域分解有限元弱形式

永磁磁場(chǎng)的磁標(biāo)勢(shì)有限元分析的邊值問(wèn)題為

其中， μ為磁導(dǎo)率，φ 為標(biāo)量磁位，Hc為永磁材料的矯頑場(chǎng)強(qiáng)， n 為永磁邊界的外法向矢量， ?為求解域，Γv為 真空的邊界，Γm為永磁材料的邊界。

為了便于推導(dǎo)區(qū)域分解有限元弱形式，將求解域分成2個(gè)子區(qū)域，如圖1所示，其中 ?1, ?2代表2個(gè)子區(qū)域， Γv1，Γv2為2個(gè)子區(qū)域的真空邊界，Γm1, Γm2為2個(gè)子區(qū)域的永磁邊界， Γ為2個(gè)子區(qū)域的交界面，n1, n2為交界面上的外法向矢量。

擴(kuò)展單個(gè)區(qū)域的邊值問(wèn)題式(1)到2個(gè)子區(qū)域，可以得到

其中，式(4)和式(5)用來(lái)保證區(qū)域交界面上物理量的連續(xù)性，但是其收斂性很差，常用的Robin傳輸條件也是通過(guò)兩式的線性變換得到的。本文拋棄了之前的傳輸條件構(gòu)造方式，而是從接觸熱阻的定義[16]出發(fā)，開創(chuàng)性地提出了一種新型傳輸條件，其具體表達(dá)式為

其中，γ 表示區(qū)域交界面上物理量連續(xù)程度的物理量，理論上當(dāng) γ無(wú)窮大時(shí)，區(qū)域交界面上的物理量就會(huì)完全連續(xù)[17]，在實(shí)際應(yīng)用過(guò)程中，只需要取一個(gè)比較大的值，106量級(jí)基本可以滿足需求。

為了推導(dǎo)2個(gè)區(qū)域的磁標(biāo)勢(shì)有限元弱形式，用式(10)和式(11)代替式(4)和式(5)，可以得到殘差表達(dá)式

圖1 單個(gè)區(qū)域分成2個(gè)子區(qū)域示意圖

首先定義體積分和面積分如式(20)

其中， u 是權(quán)函數(shù)， v 為殘差項(xiàng)。由式(12)—式(19)可以得到如式(21)所示的線性組合方程

其中，c1, c2, c3, c4, c5和 c6為待定系數(shù)。

由格林公式，式(21)中的前兩項(xiàng)可以寫為

對(duì)于式(22)中第1類邊界條件項(xiàng)，由于基函數(shù) wi具有任意性，令wi=0，有

也就是說(shuō)，第1類邊界條件項(xiàng)在有限元弱形式推導(dǎo)過(guò)程中可以不考慮，但需要采用強(qiáng)強(qiáng)加的方式添加到最終的有限元矩陣方程里。接下來(lái)，對(duì)于區(qū)域交界面上的項(xiàng)有

對(duì)于永磁的邊界項(xiàng)有

令c1=?1, c2=?1, c5=?1, c6=?1，可以得到如式(26)所示的有限元弱形式表達(dá)式

不難發(fā)現(xiàn)式(26)可以擴(kuò)展到任意多個(gè)子區(qū)域的情形。

2.2 有限元離散

由于四面體單元在處理復(fù)雜邊界時(shí)具有良好的適應(yīng)性，同時(shí)為了使用較少的網(wǎng)格和自由度得到較高的計(jì)算精度，本文采用了基于四面體單元的2階疊層標(biāo)量基函數(shù)進(jìn)行有限元離散。首先定義體積坐標(biāo)，四面體內(nèi)的體積坐標(biāo)滿足式(27)

2階疊層標(biāo)量基函數(shù)則從體積坐標(biāo)出發(fā)，由1階基函數(shù)構(gòu)造出2階基函數(shù)，包括了4個(gè)頂點(diǎn)基函數(shù)和6個(gè)邊基函數(shù)[18]，構(gòu)造形式為

用2階基函數(shù)去離散式(26)，便可以得到式(29)所示的矩陣方程

由于式(29)的系數(shù)矩陣是對(duì)稱正定的，本文采用了包括塊雅可比和多波前塊不完全楚列斯基分解[19]的兩層預(yù)處理的共軛梯度算法來(lái)進(jìn)行矩陣方程的求解，相比于傳統(tǒng)有限元法，可以大幅提高求解效率和減少內(nèi)存消耗。經(jīng)過(guò)有限元分析得到標(biāo)量磁位φ 的值，就可以根據(jù)式(30)得到磁感應(yīng)強(qiáng)度 B的值

2.3 區(qū)域劃分

本文使用METIS軟件包[20]進(jìn)行區(qū)域的劃分，經(jīng)過(guò)大量的對(duì)比分析，區(qū)域劃分過(guò)程需要考慮以下幾點(diǎn)：

(1) 區(qū)域劃分的個(gè)數(shù)對(duì)并行效率的影響很大，隨著子區(qū)域數(shù)目的增加，并行計(jì)算效率會(huì)逐步增加，雖然理論上子區(qū)域數(shù)目可以隨意取值，但是實(shí)際上隨著區(qū)域數(shù)目的進(jìn)一步增加，線程之間的資源競(jìng)爭(zhēng)會(huì)更加激烈，并且線程切換花銷也隨之增大，會(huì)使得并行效率降低。

(2) 考慮到求解過(guò)程中每個(gè)子區(qū)域矩陣都需要進(jìn)行預(yù)處理，為了避免線程等待，劃分區(qū)域時(shí)應(yīng)盡量使得每個(gè)子區(qū)域大小相當(dāng)。

(3) 劃分區(qū)域時(shí)應(yīng)盡量使得區(qū)域交界面數(shù)量少，可以加快矩陣求解過(guò)程收斂速度，從而提高計(jì)算效率。

3 仿真實(shí)例

本節(jié)通過(guò)仿真多個(gè)微波管永磁聚焦系統(tǒng)，并與商業(yè)軟件 Maxwell對(duì)比，來(lái)驗(yàn)證所提出的基于有限元的區(qū)域分解方法的準(zhǔn)確性和高效性。區(qū)域分解法中的因子γ 取值為108，多波前塊不完全楚列斯基分解殘差為10–4，預(yù)處理共軛梯度法收斂殘差和Maxwell一樣為10–6。所有的仿真計(jì)算都是在一臺(tái)小型工作站(Windows 10, Intel Xeon 5122 3.60 GHz 3.59 GHz 雙處理器，16 threads, 128 GB RAM)上完成的。

3.1 單周期結(jié)構(gòu)仿真

如圖2，選取了一個(gè)典型的單周期結(jié)構(gòu)[4]，采用的區(qū)域劃分方式為沿著Z軸方向并盡量使得每個(gè)子區(qū)域的大小相當(dāng)。首先與商業(yè)軟件Maxwell進(jìn)行精度上的比較，圖3繪制了本文提出的區(qū)域分解法和Maxwell軟件軸切面上的磁感應(yīng)強(qiáng)度云圖分布，可以看到其磁感應(yīng)強(qiáng)度云圖分布趨勢(shì)一致。由于磁鋼、極靴與真空交界處磁場(chǎng)變化比較劇烈，此處的

圖2 單周期結(jié)構(gòu)計(jì)算模型及區(qū)域分解示意圖

磁感應(yīng)強(qiáng)度由于網(wǎng)格因素會(huì)產(chǎn)生奇異值，所以將顯示范圍固定為0～1 T。

如圖4所示，將軸線上的磁場(chǎng)與Maxwell進(jìn)行對(duì)比，可以看到其吻合情況很好，在兩個(gè)峰值點(diǎn)處的相對(duì)誤差分別為0.12%和0.06%。此外還與Maxwell進(jìn)行了計(jì)算性能對(duì)比，如表1所示，隨著子區(qū)域個(gè)數(shù)的增加，計(jì)算時(shí)間和內(nèi)存相比于Maxwell的優(yōu)勢(shì)越來(lái)越大，當(dāng)劃分為20個(gè)子區(qū)域時(shí)，時(shí)間加速比達(dá)到了11.4倍，而內(nèi)存只有Maxwell的53.5%. 這里需要注意的是：在線程數(shù)更多的計(jì)算機(jī)上，區(qū)域數(shù)的增加會(huì)帶來(lái)更加優(yōu)越的計(jì)算性能，這里劃分到20個(gè)區(qū)域已經(jīng)足以說(shuō)明所提出的區(qū)域分解方法具有非常好的計(jì)算優(yōu)勢(shì)。

3.2 Wiggler結(jié)構(gòu)仿真

本實(shí)例考慮兩周期Wiggler結(jié)構(gòu)[4]的仿真計(jì)算，磁鋼材料為SmCo28，為了展示區(qū)域分解法的計(jì)算效率并兼顧區(qū)域劃分的方便快捷，將計(jì)算模型固定劃分為20個(gè)子區(qū)域，其計(jì)算模型和區(qū)域分解示意圖如圖5所示。首先進(jìn)行磁感應(yīng)強(qiáng)度云圖的對(duì)比，如圖6所示，可以看到兩者軸切面上的云圖分布趨勢(shì)相同。為了進(jìn)一步證明所提出的區(qū)域分解方法的準(zhǔn)確性，選取了軸線上的磁感應(yīng)強(qiáng)度與Maxwell進(jìn)行對(duì)比，由圖7可以看到其吻合程度非常好，并且其峰值處的相對(duì)誤差最大不超過(guò)0.14%。

此外，如表2所示，與Maxwell對(duì)比了3組網(wǎng)格數(shù)目相當(dāng)情況下的計(jì)算時(shí)間和峰值內(nèi)存，可以看到3組不同實(shí)例下的時(shí)間加速比分別為3.7, 3.2和4.2，但是其峰值內(nèi)存分別只有Maxwell的77%, 82%和73%，充分證明了所提出的區(qū)域分解法的高效性。

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

圖3 區(qū)域分解法與Maxwell軟件軸切面磁感應(yīng)強(qiáng)度云圖對(duì)比

圖4 單周期結(jié)構(gòu)軸線磁場(chǎng)B z 分布

表1 單周期結(jié)構(gòu)區(qū)域分解法與Maxwell軟件性能對(duì)比

圖5 Wiggler計(jì)算模型和區(qū)域劃分示意圖

圖6 區(qū)域分解法與Maxwell軸切面磁感應(yīng)強(qiáng)度云圖分布對(duì)比

圖7 Wiggler結(jié)構(gòu)軸線By分布和峰值相對(duì)誤差曲線

表2 Wiggler結(jié)構(gòu)區(qū)域分解法與Maxwell性能對(duì)比

本文針對(duì)微波管中的永磁聚焦系統(tǒng)仿真，提出了一種先進(jìn)的基于有限元的區(qū)域分解求解技術(shù)，并對(duì)其理論進(jìn)行了詳細(xì)的描述，還給出了實(shí)際應(yīng)用中區(qū)域劃分的相關(guān)原則和技巧。通過(guò)對(duì)多個(gè)永磁結(jié)構(gòu)的建模與仿真計(jì)算，并與商業(yè)軟件Maxwell進(jìn)行詳細(xì)的對(duì)比，驗(yàn)證了本文提出的區(qū)域分解方法的準(zhǔn)確性和高效性。本文給出的針對(duì)永磁聚焦系統(tǒng)仿真的區(qū)域分解求解技術(shù)后續(xù)有望集成到微波管模擬器套裝[8]中，為管型設(shè)計(jì)師提供更好的仿真設(shè)計(jì)平臺(tái)。