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

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

1 引言

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

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

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

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

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

永磁磁場的磁標(biāo)勢有限元分析的邊值問題為

其中， μ為磁導(dǎo)率，φ 為標(biāo)量磁位，Hc為永磁材料的矯頑場強(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ū)域的邊值問題式(1)到2個(gè)子區(qū)域，可以得到

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

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

為了推導(dǎo)2個(gè)區(qū)域的磁標(biāo)勢有限元弱形式，用式(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)可以寫為

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

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

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

2.3 區(qū)域劃分

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

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

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

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

3 仿真實(shí)例

本節(jié)通過仿真多個(gè)微波管永磁聚焦系統(tǒng)，并與商業(yè)軟件 Maxwell對比，來驗(yàn)證所提出的基于有限元的區(qū)域分解方法的準(zhǔn)確性和高效性。區(qū)域分解法中的因子γ 取值為108，多波前塊不完全楚列斯基分解殘差為10–4，預(yù)處理共軛梯度法收斂殘差和Maxwell一樣為10–6。所有的仿真計(jì)算都是在一臺小型工作站(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)度云圖分布趨勢一致。由于磁鋼、極靴與真空交界處磁場變化比較劇烈，此處的

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

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

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

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)度云圖的對比，如圖6所示，可以看到兩者軸切面上的云圖分布趨勢相同。為了進(jìn)一步證明所提出的區(qū)域分解方法的準(zhǔn)確性，選取了軸線上的磁感應(yīng)強(qiáng)度與Maxwell進(jìn)行對比，由圖7可以看到其吻合程度非常好，并且其峰值處的相對誤差最大不超過0.14%。

此外，如表2所示，與Maxwell對比了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é)束語

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

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

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

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

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

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

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

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