基于有限元的區(qū)域分解方法在永磁聚焦系統(tǒng)仿真中的應用

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

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

1 引言

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

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

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

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

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

永磁磁場的磁標勢有限元分析的邊值問題為

其中， μ為磁導率，φ 為標量磁位，Hc為永磁材料的矯頑場強， n 為永磁邊界的外法向矢量， ?為求解域，Γv為 真空的邊界，Γm為永磁材料的邊界。

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

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

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

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

為了推導2個區(qū)域的磁標勢有限元弱形式，用式(10)和式(11)代替式(4)和式(5)，可以得到殘差表達式

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

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

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

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

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

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

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

對于永磁的邊界項有

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

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

2.2 有限元離散

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

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

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

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

2.3 區(qū)域劃分

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

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

(2) 考慮到求解過程中每個子區(qū)域矩陣都需要進行預處理，為了避免線程等待，劃分區(qū)域時應盡量使得每個子區(qū)域大小相當。

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

3 仿真實例

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

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

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

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

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

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

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

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