有限元外推插值法在非线性磁场计算中的应用研究

1　引言

有限元外推插值法是近年来出现的一种计算电磁场的新方法^［1］。从它被引入到线性电磁场数值计算中的情况可知，它能在占用计算机时和内存量增加不多的情况下，显著提高数值解的整体逼近精度，而且方法易于为人们接受和掌握。那么，有限元外推插值法在人们十分关注的非线性场计算中是否也有类似的突出优点呢?基于这一想法，就有限元外推插值法计算非线性电磁场的可能性，位函数和场矢量的计算格式、计算精度、计算量和计算机内存占用等方面，我们做了初步的研究。此处介绍该方法形成和运用的基本情况。

2　数学基础

2.1　变系数椭圆边值问题有限元近似解的外推方法

电磁场的非线性边值问题在数学上属变系数椭圆边值问题，因而我们首先考虑变系数椭圆边值问题有限元近似解的外推方法。
对变系数椭圆边值问题

其中，,变系数a_ij、a≥0连续可微，i，j=1,2，…，n。
假定Ω的边界Ω分段属于三次连续可微，且在角点是局部凸的(包括凸多边形和光滑域)，则这个边值问题在Hilbert空间H¹₀∩H²有唯一解，且

‖u‖_2,Ω≤c‖f‖_0，Ω

令S^h为在Ω内的分段线性连续函数集合，即S^hH¹₀(Ω)。定义

对Ω进行三角单元剖分π_H，令M_H为与剖分相应的节点集，则原变系数椭圆边值问题的线性有限元解u_H∈S^h满足
B_h(u_H,ν)=f_h(ν)　ν∈S^h
即B_h(u-u_H,ν)=0ν∈S^h

细分π_H，使得细分后新增加的节点为π_H单元中对应边的中点，得到一新的三角剖分π_h，令M_h为细分后新增加的节点集，对应的线性有限元解u_h∈S^h/2，将细分后的节点集确定为M_h=M_h/M_H，则在三角单元顶点上定义有限元外推解

Eu_h=(4u_h-u_H)/3z∈M_h

假设z₁，z₂∈M_H为三角单元两相邻的顶点，且其中至少有一个顶点不在边界上，在它们的中点z₁₂=(z₁+z₂)/2∈M_h上定义有限元外推解

我们由此可得定理1：在以上的条件下，如剖分是没有内极点的分片强正规三角剖分，且u∈W^3，∞∩H¹₀，则如下的有限元外推解的误差估计式成立

有关的数学证明见文献［2］。

2.2　有限元外推插值的有关理论

在低阶有限元解的基础上，采用在单元内构造高阶插值函数的方法，来提高待求函数的精度，使待求位函数及场函数的整体逼近程度有进一步的提高。
在剖分π_H中将三角单元分成4个全等的小三角形而得到部分π_h，u_H为π_H的线性有限元解(即按一阶插值有限元计算所得的解)，i′_H表示π_H上的分片二次插值算子，取

则插值后，有限元解的导数有超收敛

即线性有限元解在三角单元内的导数本来只有一阶收敛，但对该有限元解进行二次插值后，导数便具有二阶超收敛。
将有限元外推法和有限元插值法相结合，我们称为有限元外推插值法，有定理2：对于分片强正规三角剖分，分片二次插值后导数有二阶超收敛，则有限元外推插值的导数有3阶超收敛，即

‖u^*-u‖₁≤ch³‖1nh‖

式中　u^*——有限元外推解
有关的数学证明见文献［1，3，4］。

3　计算非线性场的有限元外推插值法计算格式

将上述变系数椭圆边值问题的有限元外推插值法运用于二维非线性电磁场边值问题(此时x₁=x,x₂=y)。若在区域Ω上网格剖分是分片强正规的，当我们分别对初次剖分和加密一次细分场域基础上所得的非线性有限元方程进行计算，再以这些结果进行外推插值计算，就可得到非线性场有限元外推插值法的计算结果。
对图1所示的计算区域Ω，先采用较粗略的三角单元剖分π_H，得到E个单元，n个节点，加密剖分后为π_h，它有4E个单元，N个节点。设△123是π_H中的某一三角元，将△123的各边中点分别相连，便得到细分π_h的四个单元△165，△624，△456，△543，见图1。又设细分前△123的各节点位值的线性有限元解为u_Hi(i=1，2，3)，细分后各节点位值的线性有限元解为u_hi(i=1，2，…，6)，根据外推计算公式，可得单元△123中各节点位值的外推解u_i为

将其推广到整个域Ω，可得到求解全域位值外推解的矩阵形式

图1
Fig.1

式中　u——有限元外推解的N阶列矢量
u_H——粗分线性有限元解的n阶列矢量
u_h——为加密剖分后线性有限元解的N阶列矢量
O——n×(N-n)阶零矩阵
I——(N-n)阶单位阵
B——(N-n)×n阶矩阵
取I₁为n阶单位矩阵，则的元素为

粗剖分单元边的中点(如节点6)与该边的顶点(节点1，2)称为关联，很显然，粗剖分单元边中点仅与该边的两个顶点相关，而与其余节点不关联。因加密剖分新增加的节点为粗剖分单元边的中点，所以B矩阵中，每行只有两个非零元素。
对剖分π_H的单元运用二次插值，其插值函数值为

式中　u_i——单元上的有限元解(i=1，2，3，4，5，6)
N_i——插值基函数
若上式中的u_i用有限元外推解代入，便形成了有限元外推插值算法。
在计算非线性场问题时，对于含有非线性媒质的场域进行粗分和加密一次剖分，分别导出相应线性有限单元上的非线性有限元方程，然后用牛顿拉夫森迭代法计算这两个非线性方程组，再根据计算结果进行外推插值计算。

4　算例分析

设三种均匀媒质构成的二维静磁场区域Ω=Ω₀+Ω₁+Ω₂，其尺寸如图2所示，单位为cm。内部区域Ω₀为载流导体媒质，电流密度J=20.0A/cm²，磁阻率为ν₀=0.80×10⁴A/(T^.cm)；区域Ω₂为磁阻率ν₂=0.00008×10⁴A/(T^.cm)的线性媒质；区域Ω₁为非线性铁磁媒质，不饱和时磁阻率为ν=0.0005×10⁴A/(T^.cm)，其磁化曲线如表1。

图2
Fig.2

表1　磁化曲线