基于有限元的区域分解方法在永磁聚焦系统仿真中的应用

谢 鹏 徐 立 尹俊辉 杨中海 李 斌

(电子科技大学电子科学与工程学院 成都 610054)

1 引言

微波管广泛应用于卫星通信、雷达系统、电子对抗以及科学研究领域[1–3]，其中永磁聚焦系统常作为外部磁场用于电子束聚焦[4]。目前，在大量的科学及工程应用中，有限元方法针对复杂结构仍然是主流的数值分析工具[5,6]。在永磁聚焦系统仿真设计领域，国外Maxwell[7]是最流行的3维有限元仿真设计商业软件，国内电子科技大学开发的微波管模拟器套装[8]里面的永磁聚焦模拟器[4]则是具有国家自主知识产权的3维磁场仿真设计软件。它们都拥有永磁聚焦系统仿真能力，但是当求解大规模和多尺度问题时，其有限元矩阵求解通常会花费大量的时间和内存，有时甚至由于缺少有效的预处理导致矩阵无法求解。

非重叠区域分解法采用“分而治之”的思想，将原求解域划分成若干个互不重叠的子区域进行求解，其具有天然的并行性，非常适合仿真大型复杂结构[9]。本文采用基于磁标势有限元的非重叠区域分解方法来进行永磁聚焦系统的仿真计算。归根结底，基于磁标势的永磁聚焦系统仿真属于求解泊松方程边值问题。在这类边值问题的基于有限元的区域分解方法中，目前的处理方式主要有两种，一种是拉格朗日乘子类型的区域分解法[10,11]；一种是基于内罚方式的区域分解法[12,13]。前者会有两种不同类型的推导过程：一是保留拉格朗日乘子的方式，这种方式会额外增加未知量个数，并且会生成一个对称不定的鞍点矩阵系统，不利于方程的求解；二是在推导过程中消去拉格朗日乘子，这种方式会生成一个对称正定的系数矩阵，但是这种方式有可能会因为巨大的计算代价而不能显式地计算系数矩阵[14]。基于内罚的区域分解方法则不需要引入诸如拉格朗日乘子类型的辅助变量，只需要将传输条件引入基于内罚方式的有限元弱形式推导过程中，目前采用的主要是Robin传输条件，由Lions[15]首次提出，但目前基于Robin传输条件的内罚区域分解法需要考虑法向偏导项的计算，并且最终形成的有限元矩阵是非对称的。

本文提出的区域分解方法同样是基于内罚方式的，但是引入了一种新型传输条件，其来源于接触热阻的定义。相比于现有的方法，该区域分解方法的主要优势包括：(1)不需要引入多余的未知量，使得有限元矩阵维数更少；(2)有限元矩阵集成过程更加简单，只需要考虑区域交界面上的物理量，而且不需要进行法向偏导项的计算，更重要的是最终产生的有限元矩阵满足对称正定性，矩阵性质更好，非常适合采用共轭梯度法进行求解。通过对多个永磁结构的仿真计算可以发现，相比于商业软件Maxwell，本文所提出的区域分解方法在保证求解精度的同时，可以更加高效地实现对微波管永磁聚焦系统的仿真。

2 基于有限元的区域分解方法

2.1 基于内罚的区域分解有限元弱形式

永磁磁场的磁标势有限元分析的边值问题为

其中， µ为磁导率，φ 为标量磁位，Hc为永磁材料的矫顽场强， n 为永磁边界的外法向矢量， Ω为求解域，Γv为 真空的边界，Γm为永磁材料的边界。

为了便于推导区域分解有限元弱形式，将求解域分成2个子区域，如图1所示，其中 Ω1, Ω2代表2个子区域， Γv1，Γv2为2个子区域的真空边界，Γm1, Γm2为2个子区域的永磁边界， Γ为2个子区域的交界面，n1, n2为交界面上的外法向矢量。

扩展单个区域的边值问题式(1)到2个子区域，可以得到

其中，式(4)和式(5)用来保证区域交界面上物理量的连续性，但是其收敛性很差，常用的Robin传输条件也是通过两式的线性变换得到的。本文抛弃了之前的传输条件构造方式，而是从接触热阻的定义[16]出发，开创性地提出了一种新型传输条件，其具体表达式为

其中，γ 表示区域交界面上物理量连续程度的物理量，理论上当 γ无穷大时，区域交界面上的物理量就会完全连续[17]，在实际应用过程中，只需要取一个比较大的值，106量级基本可以满足需求。

为了推导2个区域的磁标势有限元弱形式，用式(10)和式(11)代替式(4)和式(5)，可以得到残差表达式

图1 单个区域分成2个子区域示意图

首先定义体积分和面积分如式(20)

其中， u 是权函数， v 为残差项。由式(12)—式(19)可以得到如式(21)所示的线性组合方程

其中，c1, c2, c3, c4, c5和 c6为待定系数。

由格林公式，式(21)中的前两项可以写为

对于式(22)中第1类边界条件项，由于基函数 wi具有任意性，令wi=0，有

也就是说，第1类边界条件项在有限元弱形式推导过程中可以不考虑，但需要采用强强加的方式添加到最终的有限元矩阵方程里。接下来，对于区域交界面上的项有

对于永磁的边界项有

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

不难发现式(26)可以扩展到任意多个子区域的情形。

2.2 有限元离散

由于四面体单元在处理复杂边界时具有良好的适应性，同时为了使用较少的网格和自由度得到较高的计算精度，本文采用了基于四面体单元的2阶叠层标量基函数进行有限元离散。首先定义体积坐标，四面体内的体积坐标满足式(27)

2阶叠层标量基函数则从体积坐标出发，由1阶基函数构造出2阶基函数，包括了4个顶点基函数和6个边基函数[18]，构造形式为

用2阶基函数去离散式(26)，便可以得到式(29)所示的矩阵方程

由于式(29)的系数矩阵是对称正定的，本文采用了包括块雅可比和多波前块不完全楚列斯基分解[19]的两层预处理的共轭梯度算法来进行矩阵方程的求解，相比于传统有限元法，可以大幅提高求解效率和减少内存消耗。经过有限元分析得到标量磁位φ 的值，就可以根据式(30)得到磁感应强度 B的值

2.3 区域划分

本文使用METIS软件包[20]进行区域的划分，经过大量的对比分析，区域划分过程需要考虑以下几点：

(1) 区域划分的个数对并行效率的影响很大，随着子区域数目的增加，并行计算效率会逐步增加，虽然理论上子区域数目可以随意取值，但是实际上随着区域数目的进一步增加，线程之间的资源竞争会更加激烈，并且线程切换花销也随之增大，会使得并行效率降低。

(2) 考虑到求解过程中每个子区域矩阵都需要进行预处理，为了避免线程等待，划分区域时应尽量使得每个子区域大小相当。

(3) 划分区域时应尽量使得区域交界面数量少，可以加快矩阵求解过程收敛速度，从而提高计算效率。

3 仿真实例

本节通过仿真多个微波管永磁聚焦系统，并与商业软件 Maxwell对比，来验证所提出的基于有限元的区域分解方法的准确性和高效性。区域分解法中的因子γ 取值为108，多波前块不完全楚列斯基分解残差为10–4，预处理共轭梯度法收敛残差和Maxwell一样为10–6。所有的仿真计算都是在一台小型工作站(Windows 10, Intel Xeon 5122 3.60 GHz 3.59 GHz 双处理器，16 threads, 128 GB RAM)上完成的。

3.1 单周期结构仿真

如图2，选取了一个典型的单周期结构[4]，采用的区域划分方式为沿着Z轴方向并尽量使得每个子区域的大小相当。首先与商业软件Maxwell进行精度上的比较，图3绘制了本文提出的区域分解法和Maxwell软件轴切面上的磁感应强度云图分布，可以看到其磁感应强度云图分布趋势一致。由于磁钢、极靴与真空交界处磁场变化比较剧烈，此处的

图2 单周期结构计算模型及区域分解示意图

磁感应强度由于网格因素会产生奇异值，所以将显示范围固定为0～1 T。

如图4所示，将轴线上的磁场与Maxwell进行对比，可以看到其吻合情况很好，在两个峰值点处的相对误差分别为0.12%和0.06%。此外还与Maxwell进行了计算性能对比，如表1所示，随着子区域个数的增加，计算时间和内存相比于Maxwell的优势越来越大，当划分为20个子区域时，时间加速比达到了11.4倍，而内存只有Maxwell的53.5%. 这里需要注意的是：在线程数更多的计算机上，区域数的增加会带来更加优越的计算性能，这里划分到20个区域已经足以说明所提出的区域分解方法具有非常好的计算优势。

3.2 Wiggler结构仿真

本实例考虑两周期Wiggler结构[4]的仿真计算，磁钢材料为SmCo28，为了展示区域分解法的计算效率并兼顾区域划分的方便快捷，将计算模型固定划分为20个子区域，其计算模型和区域分解示意图如图5所示。首先进行磁感应强度云图的对比，如图6所示，可以看到两者轴切面上的云图分布趋势相同。为了进一步证明所提出的区域分解方法的准确性，选取了轴线上的磁感应强度与Maxwell进行对比，由图7可以看到其吻合程度非常好，并且其峰值处的相对误差最大不超过0.14%。

此外，如表2所示，与Maxwell对比了3组网格数目相当情况下的计算时间和峰值内存，可以看到3组不同实例下的时间加速比分别为3.7, 3.2和4.2，但是其峰值内存分别只有Maxwell的77%, 82%和73%，充分证明了所提出的区域分解法的高效性。

4 结束语

图3 区域分解法与Maxwell软件轴切面磁感应强度云图对比

图4 单周期结构轴线磁场B z 分布

表1 单周期结构区域分解法与Maxwell软件性能对比

图5 Wiggler计算模型和区域划分示意图

图6 区域分解法与Maxwell轴切面磁感应强度云图分布对比

图7 Wiggler结构轴线By分布和峰值相对误差曲线

表2 Wiggler结构区域分解法与Maxwell性能对比

本文针对微波管中的永磁聚焦系统仿真，提出了一种先进的基于有限元的区域分解求解技术，并对其理论进行了详细的描述，还给出了实际应用中区域划分的相关原则和技巧。通过对多个永磁结构的建模与仿真计算，并与商业软件Maxwell进行详细的对比，验证了本文提出的区域分解方法的准确性和高效性。本文给出的针对永磁聚焦系统仿真的区域分解求解技术后续有望集成到微波管模拟器套装[8]中，为管型设计师提供更好的仿真设计平台。