电磁场数值计算4-西安交通大学电气工程学院解读

第四章
二维电磁场有限元分析
二维电磁场：平行平面场、轴对称场的基本概念 §4-1 电磁场的微分方程及泛函
常用的泊松方程、拉普拉斯方程和亥姆霍兹方程是方程 f u y u y x u x y x =+⎪⎪⎭
⎫
⎝⎛∂∂∂∂-⎪⎭⎫ ⎝⎛∂∂∂∂-
βαα （4-1） 的几种特殊形式，它适用于非线性、非均匀、各向异性材料的静、稳态和简谐电磁场问题。

在扩散方程中ωβj =，在波动方程中2ωβ-=。

边界条件 01
u u
Γ=
q u y u x u 3
n y y x x =+⋅⎪⎪⎭⎫
⎝⎛∂∂+∂∂Γγααe e e —电场 （4-2a ） q u x u αy u αΓy y x x =+⨯⎥⎦⎤
⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂-∂∂3
e e e γn —磁场（z A u =） （4-2b ） 0n
u =∂∂对称面
E 线与对称面平行0=n E ，或B 线与对称面垂直0=t B
0=1
Γu
E 线与对称面垂直0=t E ，或B 线与对称面平行0=n B
式（4-1）对应的泛函求极值问题为
()()0
22222221u u dS fu d Γqu u γ dS βu y u αx u αu I S ΓS y x ==-⎪⎭⎫ ⎝⎛-+⎪⎪⎭
⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂=Γ⎰⎰⎰⎰⎰min （4-3）
（1） 若存在第一类边界条件，需要专门处理。

（2） 若存在媒质分界面，变分()0u I =δ中将自动满足
（3） 不论γβαα,,,y x 是实数还是复数，式（4-3）均成立。

如果这些参数
只是实数则可用另一泛函
()()
()
min =+---+
⎪⎪⎭
⎫ ⎝⎛+∂∂+∂∂=⎰⎰⎰⎰⎰*dS fu u f d Γqu u q u γ dS
u βy
u αx u αu I S
**Γ*S y x 212121222
22 （4-4）
若为电场问题： , ,u εαϕ==
y x y
x e e E ∂∂-∂∂-
=-∇=ϕ
ϕϕ （4-5） 若为磁场问题： ,1
,A u μ
α=
=或 , ,u m μαϕ==
y x x
A y A e e A
B ∂∂-∂∂=
⨯∇= （4-6） 若为非线性问题，()⎪⎪⎪⎭⎫ ⎝
⎛
⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂==22y A x A B μμμ §4-2 有限元分析 4.2.1 离散化（前处理）
将面区域离散为许多小面单元：三角形（三节点、六节点）、四边形（四节点、八节点） 1、离散的基本要求
* 单元之间既没有间隔，又没有重叠，单元通过它们的顶点相连； * 尽量避免较小内角的单元产生，可以证明，有限元解的误差反比与最小内角的正弦。

最理想的三角形是等边三角形；
* 一个单元的3个节点不能同时在一个边界上，若64123是第一类边界，则单元①的剖分法是错误的，这将导致24边等位；
* 单元越小，数值解结果越好，但未知量多会增加计算机内存和计算量，会带来计算误差。

因此，解变化剧烈的地方剖分可以较密些。

2、编码要求
为了压缩内存和简化程序，单元和节点编码应该注意： * 单元用一组整数编码（0e ,2,1e =）；
* 节点编码用另一组整数，它反映两个信息：①局部编码，表示节点在相应单元中的位置；②全局编码，表示该节点在全域中的位置。

两个功能用1个二维数组来表示，如()e ,i n ，其中，i =1,2,3反映三角单元3个顶点m ,j ,i 局部位置，
0e ,2,1e =是单元编码，()e ,i n 的值是节点的全局编码。

例：写出4单元6节点的节点编码 注：局部编码要逆时针进行。

* 节点坐标数组：() n 1,2,i ,y ,x xy i i = * 每个单元的材料参数及源密度 4.2.2 单元插值及基函数
以三角形3节点单元为例，采用线性插值函数逼近e 单元内的未知函数
()y c x b a y ,x u e e e e ++= （4-7）
e e e c ,b ,a 为待定常数，e 单元上3个节点的未知数分别为
()()()m e m e e m m e
m
j e
j e e j j e j i
e i e e i i e i y c x b a y ,x u y c x b a y ,x u y c x b a y ,x u ++=++=++= （4-8） 从中解出e e e c ,b ,a ，再代入式（4-7），得到
()()e j 3
1
j e j e
u y ,x N y ,x u ∑== （4-9）
式中， ()
3,1,2j y c x b a 21N e
j e j e j e
e j =++∆
= （4-10） 其中，
e i
e j e m e j e i e m e i e j e j e i e m e m e i e j e i e m e j e m e i e i e m e j e
j
e m e i e m e j e i e j e m e m e j e i x x c y y b y x y x a x x c y y b y x y x a x x c y y b y x y x a -=-=-=-=-=-=-=-=-= （4-11）
单元面积 ()
e
i e j e j e i e
m
e m e j e j e
i e i e c b c b y x y x y x Δ-==
2
111121 由此可见，基函数e N 是坐标y x,的线性插值函数，只与三角形的形状及节点分布有关，故又称之为形状函数。

基函数具有性质：
（1） 定义域上任意点p 的坐标代入基函数，有
()() e p
0e p p N p N e i e
i
⎩⎨⎧=上点不在单元上点在单元
体现方程具有稀疏矩阵的特性。

（2） 若p 点是单元节点()021e , ,i p i =，则()e i i e u p u =，基函数为
()k ,1,2,j i, j i j i
p N j e
i
=⎩
⎨⎧≠==01
表明：① 场点在节点i 上，()e
i i i e u ,y x u =；
② 当场点在节点i 的对边上时，()0=x,y N e i 。

所以，一个单元边上的
()y ,x u e 值与它的相对节点处的u 值无关，它仅由该边两端点处的u 值决定。

这就保证了单元两侧解的连续性。

（实际上，()e e jm e i x,y N ∆∆=/，()e e im e j x,y N ∆∆=/，()e e ij e m x,y N ∆∆=/，是面积坐标
()()[]
()
e e i e e i e
i e
j e m e m e j e e j e m e e m e j e
m
e m e j e j e
e e jm y c x b a y x y x y x x x y y y x y x y x Δ++=
-+-+-==
2
12
111121 一般而言，离P 点最近的节点的电位值所占比例最大。

与式（4-10）对照可知
()e e jm e i x,y N ∆∆=/。

当P 点落在j 、m 连线上，0=∆e pjm ，所以，()0=x,y N e i 。

当P 点落在i 节点上，e e pjm ∆=∆，()1=x,y N e i ）
这样，每个三角形单元的()y ,x u e 由该单元上各顶点的值的线性插值去逼近。

另外，由于两个相邻单元公共边和公共节点上的值相等，就有可能保持线性插值函数的连续性。

若为电场，根据式（4-5），电场强度的数值解为
()
e
j 3
1
j y e j x e j e
e
u c b 21
∑=+∆
-=e e
E （4-12）
显然，单元内的电场强度是一常矢量，为了提高精度，常采用高阶插值。

若为磁场，根据式（4-6），磁场强度的数值解为
()
e
j 3
1
j y e j x
e j
e
e
u b c 21
∑=-∆
=e e
B （4-13）
同理，单元内的磁感应场强度是一常矢量。

4.2.3 建立有限元方程
不论是基于变分原理的有限元法，还是伽辽金有限元法，最终的有限元方程是一样的
F KU =
其中，系数矩阵中的元素及右端列向量中的元素分别为
⎰=⎪⎪⎭
⎫ ⎝⎛+∂∂∂∂+∂∂∂∂=e s e j e i e j e i y e j e i x e ij 1,2,3j i, dxdy N N y N y N x N x N K βαα （4-14）
1,2,3i dxdy fN F e S
e i e i ==⎰⎰ （4-15）
设单元中媒质参数不变，为e e e
y e x f , , ,βαα（若不是常数，可取相应的单元的平
均值），通过三角形上的数值积分公式
() ！c b a a!b!c!dxdy N N N e e
S c b a
e 12
223
2
1
∆=∆+++=⎰⎰ （4-16） 其中，a,b,c 分别为1，1，0。

若,,021==c N N 上式积分等于6e ∆
根据基函数的性质，即积分单元区域中，只有与节点i ，或节点j 关联的单元参加积分，其余为零，可以得到式（4-14）、式（4-15）的解析表达式
()
()ij e e
j e i e y e j e i e x e
e
ij
δβΔc c αb b αΔ
K +++=11241 （4-17） e
e e
i f 3
F ∆= （4-18）
矩阵形式为（电场：εαα==y x ，且0=β）
[]
⎥⎥⎥⎥⎦
⎤⎢⎢
⎢⎢⎣⎡++++++∆=2m 2m j m j m i m i m 2
j 2j i j i j 2i 2i e e
c b c c b b c c b b c b c c b b
c b 4K 对称ε （4-19） []
⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡∆=1113
1F e c e e
ρ （4-20） 以下图为例，K 为6⨯6阶方阵
① 6,2,1j ,1i ==
11K —与节点1，
1关联的单元为（1），()11111K K = 12K —与节点1，
2关联的单元为（1），()11212K K = 13K —与节点1，3关联的单元为0，0K 13=
14K —与节点1，4关联的单元为（1）
，()
11414K K = 15K —与节点1，5关联的单元为0，0K 15= 16K —与节点1，6关联的单元为0，0K 16=
② 6,2,1j ,2i ==
21K —与节点2，1关联的单元为（1）
，()
12121K K = 22K —与节点2，2关联的单元为（1）、（2）、（3），()()()322
22212222K K K K ++= 23K —与节点2，3关联的单元为（3），(3)2323K K =
24K —与节点2，4关联的单元为（1）、（2），()()224
12424K K K += 25K —与节点2，5关联的单元为（2）、（3），325
(2)2525K K K += 26K —与节点2，6关联的单元为0，0K 26= 以此类推，最终有
∣----------- D -------------∣
()()()
()()()()()()()()()
()()()
()
()()()()()()()()()()()()()()()()()
()()()⎥⎥
⎥⎥
⎥⎥⎥
⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡++++++++ ++++= K K K 0 0 0K K K K K K K K K 0 K K K K K K 0 K K K 0 K 0 K K 0 0 K K K K K K K K K 0 0 K 0 K K K 466465464456455355255454254353352252446445245444244144242142141335333332325225224124323322222122121114112111 ()()()()()()
()()()()()()⎥⎥⎥⎥
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡++++++=464535254424143332221211F F F F F F F F F F F F F 系数矩阵K 的4个特点：
（1） K 是对称阵，ji ij K K = （2） 对角线元素占优，即
ij ii ii K K ,0K >>，且为正定阵；
（3） K 是稀疏矩阵，第i 行的n 个元素中，只有与节点i 关联的节点上是非
零元素，其余均为零，所以存储量小；
（4） K 中非零元素的位置呈现带状，带宽为D=1+H , H 指三角形单元节
点编码最大差值的绝对值，如图中单元①或②中节点的最大差值为H = 4-1=3或H = 5-2=3，带宽D=4。

不同的求解方法有不同的存储技术，大体分为：
（1） 等带宽存储：下三角带宽内的各元素依次（按行或按列）存储； （2） 变带宽存储：下三角带宽内从每一行第一个非零元素开始依次存储； （3） 非零元存储：下三角带宽内按行只存储非零元素。

前两者用于直接解法，后者用于迭代解法，如ICCG 法。

4.2.4 边界条件的处理 1、第一类边界条件处理
设节点3，5，6位于第一类边界上，其值分别为653p ,p ,p ，首先强加33p u =
⎥⎥
⎥
⎥⎥⎥⎥
⎥
⎦⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢
⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡65432165432166656456
55545352464544424135333225242322 2114121 1F F F F F F u u u u u u K K K 0 0 0K K K K K 0K K K 0 K K 0 K 0 K K 0 0 K K K K K 0 0 K 0 K K 做法一：令333i 33p F 1,2,4,5,6,i ,0K ,1K ====，这样将导致系数矩阵不对称。

做法二：令1,2,4,5,6
i ,0K K ,p K F F 3i i33i3i i ===-→, 令3333p F ,1K ==，并以此类推，方程为
⎥⎥⎥⎥
⎥⎥⎥⎥
⎦⎤
⎢⎢
⎢⎢⎢
⎢⎢⎢⎣⎡---------=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡65
64654534343
6265253232
6165153131654321444241242221141211100000010000000000100000000p p p K p K p K F p p K p K p K F p K p K p K F u u u u u u K K K K K K K K K 显然，方程可以降阶，且依然对称
⎥⎥⎥
⎦⎤
⎢⎢⎢⎣⎡----=⎥⎥⎥⎦⎤⎢
⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡64654545253232
1421444241242221141211p K p K F p K p K F F u u u K K K K K K K K K
2、第三类边界条件处理
第二类边界条件是第三类边界条件的特殊情况，故只讨论后者。

若为齐次边界条件，不需要增加任何项。

若为非齐次边界条件，需要增加一项
()d Γqu γu u I Γb ⎰⎪⎭
⎫
⎝⎛-=321 （4-21）
设边界离散后的单元边有s M 个，上式近似表示为每一线段积分的叠加
()()
∑==s
M 1
s s s b b u I u I （4-22）
每一线段上的未知函数可以近似为
∑==2
1j s j s j s
u N u （4-23） 式中，ς-=1N s 1，ς=s
2N ，s l /l =ς，其中ς是沿线段从节点1到节点2的归
一化距离，即节点1处，1u u s == ,0ς，节点2处，21u u ςs == ,，两个节点之间s u 按线性规律变化，s l 是线段长度。

将式（4-22）代入式（4-21），并对s i u 求导，得
∑⎰⎰=-=∂∂21
j 1010s s i s
s j s i s j s i s b d l qN d l N N u u I ςςγ （4-24） 写成矩阵形式
[]{}{}
s
s s s s b F u K u I -=⎭
⎬⎫⎩⎨⎧∂∂ （4-25） 式中系数阵是对称的
1,2j i,
d l N N K 1
s
s j s i s ij
==⎰ςγ （4-26） ⎰==1
s s i s i 1,2i d l qN F ς (4-27)
如果q ,γ在线段上是常数，分别用s s q ,γ表示，则可用解析解表示
()ij s
s s ij 1l 61K δγ+=
（4-28） s s s i l q 2
1
F = （4-29）
式中，⎩⎨⎧≠==j i
0j i
1ij δ
将式（4-27）与式（4-17）合并，式（4-28）与式（4-18）合并，得到最终方程
F KU = （4-30）
§4-3 轴对称场
电场与磁场的轴对称问题处理方法略有不同，分别介绍。

4.3.1 电场中的轴对称问题
设Z 轴为对称轴，泊松方程在圆柱坐标下的表达式为
c z z 1ρϕερϕερρρ=⎪⎭
⎫
⎝⎛∂∂∂∂-⎪⎪⎭⎫ ⎝⎛∂∂∂∂-
等式两边同乘以半径ρ，便可得到
ρρϕερρϕερρc z z =⎪⎭⎫ ⎝⎛∂∂∂∂-⎪⎪⎭
⎫ ⎝⎛∂∂∂∂-
（4-31） 做下面替换
ρρϕϕc , ======ερ , f αz , αρ , y x y x （4-32）
便可得到与平行平面场中泊松方程形式完全相同的方程。

也就是说，前面建立的有限元方程可以直接应用于该问题的求解。

轴对称情况下的开域问题可用点源公式设置人工边界
ϕϕr 1r -≈∂∂ （如r q πεϕ4= ， ϕπεϕr
r q r r 141-=-=∂∂） (4-33) 式中，22z r +=ρ。

设ϕ的源仅限于0r =附近，当∞→r 时，ϕ的渐进式为
r
A
≈
ϕ （4-34） 与式（4-14）、式（4-15）及式（4-17）、式（4-18）比较有
321,,i,j ρdρdz z N z N ερN ρN εK e S e j e i e j e i e ij =⎪⎪⎭
⎫ ⎝⎛∂∂∂∂+∂∂∂∂=⎰⎰
1,2,3i dz d N F e S
e i c e i ==⎰⎰ρρρ
[]
⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡++++++++∆=2m 2m j m j m i m i m 2
j 2j i j i j 2i 2i m j i e e
c b c c b b c c b b c b c c b b c b 34K 对称ρρρε
[]⎥⎥
⎥⎦⎤
⎢⎢⎢⎣⎡++∆=1113
3
1F m j i e c e e ρρρρ m j i ρρρ、、是三角形的3个顶点的半径。

i i c b ,等中的z y x →→,ρ替代。

电场强度(三角形3节点) ()
∑=+-=3
1
e e
21
E j y e j x
e j
e j
e
c b Δ
ϕe
（4-35）
（第七次课）
4.3.2 磁场中的轴对称问题
在磁场中，磁矢位满足的是矢量泊松方程，
z 2
z 2222222
A A A 2A A A 2A ∇+⎪⎪⎭
⎫ ⎝⎛-∂∂+∇+⎪⎪⎭⎫ ⎝⎛-∂∂-∇=∇e e e A ρφρρφρφρφφρφρρ 轴对称场中，电流只有φ方向的分量φφe J J =，那么()φφρe ,A z A =，考虑非线性磁介质，泊松方程为
()φφφμρρρμρJ z A 1z A 1=⎪⎪⎭
⎫
⎝⎛∂∂∂∂-⎥⎦⎤⎢⎣⎡∂∂∂∂- （4-36） 变形后
()()φφφρμρρρρμρJ z A 1z A 1=⎪⎪⎭
⎫ ⎝⎛∂∂∂∂-⎥⎦⎤⎢⎣⎡∂∂∂∂-
（4-37） 做下面替换
φφρρμ
ααρJ f ,A u ,1
z ,y ,x y x ===
=== （4-38）
电场中的边界条件与平行平面场的表达形式相同，不再重述。

下面讨论φρA 作为变量的边界条件是否相应于二维场讨论的一般形式 1、第一类边界条件
可以证明，轴对称场中，磁感应线方程为 C A φ=ρ，因此第一类边界条件
11
f A φ
=Γ → 11
f A φ
ρρ=Γ （4-39）
2、 第二类边界条件
若已知边界上切向磁密t B （z B 或r B ，即已知面电流K ），第二类边界条件是非齐次的，可以写为
()φe -B e A e t n n B =⨯=⨯∇⨯ （4-40）
由于 ()z A 1z A e e A ρ
ρρφρφ
∂∂+∂∂-=⨯∇
其中 z A ρz
A B φφ∂∂-
=∂∂-
=ρρ1 ()ρ
ρA ρB φz ∂∂=1 所以得到（n e 是roz 平面边界上的法向分量）
()()()()()z n φρn e ,e 1e e ,e 1e A e sin sin n ρ
ρA ρz ρA ρφφφ
∂∂-∂∂-=⨯∇⨯ 由于 ()()()()()()z n n n ,cos z
A ,cos A A n
A e e e e e ∂∂+
∂∂=
⋅∇=∂∂φρφφφρρ
ρρρ
又因
()()z n ,cos sin e e e ,e ρn = ， ()()z n e ,e e e sin ,cos n =ρ
所以 ()()φφφρρe e A e t n B n A 1-=∂∂-=⨯∇⨯ 所以，边界条件为
()t B n
A 1=∂∂φρρ 或 ()ρρφt
B n A =∂∂ （4-41） 显然，这是第二类边界条件。

若为齐次边界条件
()01=∂∂n
A φρρ
3、 分界面衔接条件
分界面连续边界条件 2n 1n B B =，2t 1t H H =，用磁矢位表示为
21A A φφρρ= ，
()()n
A 1n A 12211∂∂=∂∂φφρρμρρμ （4-42）
这样，以φφρA 'A =为变量的所有边界条件和连续性衔接条件均与ϕ的条件相对应，有限元方程自动满足分界面衔接条件。

单元磁感应强度 ()
∑=-=
3
1
e e
21
B j z e j ρ
e j
e φj
e
b c 'A Δ'e
（4-43） 实际单元磁感应强度 ()
∑=-=
=3
1
e e 21
B B j z e
j ρe j e φj e
e e e e
b c 'A ρΔρ' （4-44） 磁感应线就是等'A φ线 const A 'A ==φφρ （4-45） 节点的磁矢位 ρ
φφ'
A A = （4-46）
§4-4 非线性方程组解法
磁场多为非线性问题，材料参数是场量的函数，()B μμ=，而()A B B =，因此，有限元方程是非线性方程。

需要采用专门的解法。

4.4.1 非线性材料特性的数值逼近方法
基本磁化曲线（分析恒定磁场）、磁滞回线（分析磁滞损耗）、等效磁化曲线（分析时变场）都是非线性曲线。

它表明，如果激励源是正弦电流，那么，通过有限元法计算出来的磁矢位是正弦函数，由此取旋度得到的磁感应强度也是正弦函数，但通过μB H =得到的磁场强度一定不是正弦的。

如果激励源是正弦电压，那么，由于电感是非常数、非线性的梯度电感，因此电流一定不是正弦的，由电流计算得到的磁矢位、磁感应强度都不是正弦的。

在有限元法计算中，需要将B-H 曲线用函数形式表达，常用插值法和拟合法，每一种方法又有多种表达形式：
插值法：分段线性插值、拉格朗日插值、高次多项式插值、样条插值等； 拟合法：最小二乘法、分段线性拟合、相应表面法等。

以分段插值为例：先将B-H 曲线分成m 个子区间，每个子区间两个端点分别为() B ,H i i 及() B ,H 1i 1i ++，用线性插值，在此区间上任一点的 B 为
1i i
1i i
i 1i i 1i B H H H -H B H H H -H B ++++-+-=
（4-47）
这种插值方法使每一子区间段上的插值函数只依赖于本段上的节点值，与其它节点无关，优点是计算稳定。

但由于总体光滑度不高，导致计算精度不高。

拉格朗日插值是二阶插值函数，计算精度高于线性插值。

样条插值要求节点的一阶导数连续，光滑度高，计算精度也高。

4.4.2 线性化迭代算法
将非线性方程组分解为一系列线性方程组的迭代格式，然后逐次逼近非线性方程组的解。

迭代步骤如下
线性化迭代法的特点是：方法简单，有良好的收敛性，但收敛速度较慢，尤其当磁场进入饱和区时，收敛很慢，因此，当方程阶数不是很高时用此方法比较合适。

否则用迭代法。

4.4.3 牛顿—拉夫逊（Newton-Raphson）算法
该方法也属于线性化方法，但迭代格式不同，适用于非线性方程组求解。

尤
其对饱和磁场问题收敛性能好，比线性化迭代法更为优越。

1、一个未知量的N-R 法
先从一个未知量的非线性方程的解法来说明N-R 法。

设
F KA A f ==)( （4-48） 该式的解是图4-1(a )中焦点对应的A ，设A 的近似初解为)(0A ，)(A f 用泰勒级
f (A
图4-1牛顿—拉夫逊（Newton-Raphson ）算法示意图
数在)(0A 点展开
()()()()
+-+=)()()('000A A A f A f A f （4-49）
忽略第三项以后各项，
()()()()()()
)()()()()()('000000A A J A f A A A f A f A f -+=-+≈ （4-50）
)(0J 是)(A f 在)(0A 处的斜率，称为雅可比系数，是常数。

移项后有
()
()()()
)()()()(0000A f F A f A f A A J -=-=- （4-51）
当该式等于零时，得到解。

若不为零，求得的解作为第一个近似解)(1A ，上式方
程为
()()
)()()()(1121A f F A A J -=-
以此类推，有通式
()()
)()()()(k k k k A f F A A J -=-+1 （4-52）
当满足 ()[]
ε≤-2
)
(k A f F 时，)(k A 是最终的近似解。

这是一种变步长的方法，利用梯度确定步长,如图4-1(b)所示。

这种方法也有不收敛的情况。

当)(A f 是非单调的，可能出现)('A f 在相邻两次迭代中符号相反的情况，如图4-1(c)所示，不收敛。

形成雅可比矩阵很费时，提出一种修正的N-R 法。

2、修正的N-R 法
每次迭代都用)(0J ，一般的迭代公式为
()()
)()()()(k k k A f F A A J -=-+10 （4-53）
物理意义如图4-1(d)所示，但迭代次数增多。

3、对n 个未知函数的N-R 法
对有限元方程 F KA = 采用修正的N-R 法 []{}
{}
(){}{}
{}{}
)()
()
()
()
(k k k k KA F f F A A J -=-=-+10
非线性方程组的迭代格式为 []
()
{}(){}[]()[]()(){}(){})(k k k k k P A K J F A J =-+=+10 （4-54）
其中，()0J 是n n ⨯阶方阵。

一个三角形单元中
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=⎥⎥⎥⎥⎥⎦⎤
⎢⎢⎢⎢⎢⎣⎡=e m e m e j e m e i e
m
e m e
j e j e j e i e j e m e i e j e i e i e i e
mm e mj e mi e jm e jj e ji e im e ij e ii e A f A f A f A f A f A f A f A f A f J J J J J J J J J J （4-55）
而 e e m j i e mm e mj e mi e
jm e jj e ji e
im e ij e ii m j i A K A A A K K K K K K K K K f f f =⎥⎥⎥⎦⎤
⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡
其中，e
nq
e e nq D v K =，()
m j,i,q n, ,c c b b 41D e
q e n e q e n e
e nq =+∆
= 下面逐一推导。

由于ij K 中有μ，而μ是A 的函数，所以
()
()
e
m
e im e j e ij e i e ii e
i e e e e ii
e
e
m
e
i e im e j e i e ij e i e i e ii e ii
e
m
e im e j e ij e i e ii e i
e i e i e ii
A D A D A D A
B B v D v A A K A A K A A K K A K A K A K A A f J ++∂∂∂∂+=∂∂+∂∂+∂∂+=++∂∂=∂∂=
① 由于()
2e e e e e e e
e
e e B H B H B 1B H B B v -∂∂=⎪⎪⎭
⎫ ⎝⎛∂∂
=
∂∂，取分段线性插值函数
(
)()()
e
k
e k e
k e k e k
e
e
k
e
B B
H H B
B H H ---+=++11 则有 ()()
e k
e 1k 2e e 1k e k e k e 1k e e B B B B H B H B v --=∂∂+++ ② 因为2
2
y A x A B ⎪
⎪⎭
⎫
⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂=，前面已给出单元磁感应强度的离散式 ()
∑=-∆
=3
1
j y e
j e j x e j e j e
e
A b A c
21e e B
其模值为 2
12
31j e j e j 231j e
j e j e
e
A b A c 21
B ⎥⎥⎦
⎤
⎢⎢⎣
⎡⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛∆
=∑∑== 因此，
()()[]()()(
)[]
()
e
m
e im e j e ij e i e ii e e e
m e m e i e m e i e j e j e i e j e i e i e i e i e i e i e e e e
i
e m e m e j e j e i e i e i e m e m e j e j e i e i e
e e i e A D A D A D B
1 A b b c c A b b c c A b b c c 41B 1 b A b A b A b 2c A c A c A c 2B 41A B ++∆=+++++∆∆=+++++∆=∂∂
最终
()
e
m e im e j e ij e i e ii e e e e e ii
e
e ii
A D A D A D
B 1B v D v J ++∆
∂∂+=
同理可得
()
e
m e jm e j e jj e i e ji e e e e e jj
e
j j
e jj
A D A D A D
B 1B v D v A f J ++∆
∂∂+=∂∂=
()
e m e mm e j e mj e i e mi e e e e e mm e m m e mm
A D A D A D
B 1B v D v A f J
++∆
∂∂+=∂∂= 以及
()e
m
e jm e j e jj e i e ji e m e im e j e ij e i e ii e e e e e ij e j i e ji
e
ij
A D A D A D A D A D A D
B 1B v D v A f J J +++++∆
∂∂+=∂∂==(
)
e m
e mm e j e mj e i e mi e m e im e j e ij e i e ii e e e e e im e m i e mi
e im
A D A D A D A D A D A D
B 1B v D v A f J
J
+++++∆
∂∂+=∂∂==()
e m e mm e j e mj e i e mi e m e jm e j e jj e i e ji e e e e e jm
e
m
j e mj
e
jm
A D A D A D A D A D A D
B 1B v D
v A f J
J
+++++∆
∂∂+=∂∂=
=用通式表示
m j,i,t A D A D B 1B v D v J
t e t e
qt t e t e nt e e e e e nq
e
e nq
=⎪⎭
⎫ ⎝⎛⎪⎭⎫ ⎝⎛∆∂∂+=∑∑ （4-56）
三角形单元的右端项矩阵为
⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡e m e j e i e mm e mj e mi e
jm e jj e
ji e
im e ij e ii e e c e e c e e c e m e j e i A A A S S S S S S S S S ΔJ ΔJ ΔJ P P P 31 （4-57） 式中， A D A D B 1B v K J S t e t e jt t e t e it e e e e e nl
e nl
e
nl
⎪⎭
⎫ ⎝⎛⎪⎭⎫ ⎝⎛∆∂∂=-=∑∑
其中，m j,i,t
m ,j ,i l ,n ==，那么 i,j,m q,n A S ΔF P e n n
e qn e e e q =+=
∑31 （4-59） 最后，总雅可比矩阵J 由单元雅可比矩阵组合而成，只需将下标相同的相应项叠加即可
∑==0
e 1
e e
ij
ij J J （4-60） 雅可比矩阵J 与刚度矩阵相似，具有对称性和稀疏性，因此可以采用与K 一样的存储技术。

牛顿—拉夫逊法的特点是收敛快，远优于线性化迭代法。

缺点是： （1）要有较好的初值，否则可能迭代不收敛。

如图4- 3 所示。

（2）计算量大，每解一次方程都要重新形成雅可比矩阵，因为每一次的
e e B v ∂∂及e B 都变了。

为了加快运算速度，出现了许多的改进牛顿—拉夫逊法，如，形成一次J ，以后每次迭代都用相同的斜率J ，但K 随磁导率变化。

这样做的代价是增加迭代次数，但节省了每次计算的CPU ，因此，总的CPU 会节省。

可以采用两种方法相结合的方法：开始使用线性化方法，进入饱和区后再使用牛顿—拉夫逊法。

（第八次课）
4.4.4 预处理共轭梯度法（ICCG 法）
ICCG 法是解无约束优化问题中的一种基本优化方法。

那么，泛函极值问题
()⎰⎰=-⎥⎥⎦
⎤
⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂=S 22
S min JAdxdy dxdy y A x A v 21A I （4-61） 就可以作为无约束优化问题，未知函数()y ,x A 是优化变量，()A I 是目标函数。

ICCG 法的特点：
（a ）只需存储非零元素，节约内存；
（b ）降低系数矩阵的条件数（条件数的大小主要取决于主对角线元素数值之间的差值，差值越大，条件数越大，计算误差也越大）；
（c ）解大型稀疏矩阵方程时，CPU 时间比牛顿—拉夫逊法的时间少。

ICCG 法是在梯度法的基础上发展起来的，因此，先介绍梯度法的基本步骤： （1）第k 次迭代点用()k A 表示，设初值()0A —优化变量
（2）从()k A 出发，选择搜索方向()k P ，在该方向上寻找使目标函数下降的新点
（3）在()k P 方向上找出新点()()()()k k k 1k P A A α+=+，其中，()k α称为()k P 方向
上的步长因子
（4）如果()1k A +处的目标函数满足
()()()()
εA I A I k k <-+1 或 ()ε<∇+1k A I
迭代结束，()1k A +为优化解。

由于优化解是使目标函数取极小值的解，而0I =∇是极值解的必要条件，故可将其作为迭代终止的判据。

上述基本步骤可见，迭代过程的两个关键步骤是 ① 搜索方向；② 步长因子。

不同的()k P 和()k α将构成不同的计算方法。

1、最速下降法
该方法是将搜索方向定为目标函数下降的方向，即泛函在迭代点()k A 处的负梯度方向，即()
)(k A I 下降的方向
()
)(k k A I P -∇= （4-62） 步长因子选择 ()
()()()()()()k
T
k
k
T
k
k KP
P P P =α
（4-63）
可以证明，()k P 与()1k P +正交 ()
[]
()0=⋅+k T
k P P 1
这意味着两次搜索方向相互垂直，搜索过程呈现“锯齿”状。

图4-2 最速下降法物理含义
其物理含义如图4-2所示。

函数()y x I ,在几何上是一空间曲面，它与xoy 面相切的点即是它的零极小值点，如图4-2(a)所示。

对于空间曲面()y x I I ,=，如果用一系列平行于xoy 的平面C I =相截之，可以得到一族平面曲线，将他们投
影于xoy 面，得到如图4-2(a)所示的曲线族，称为曲面的等高线。

从初值点()00y x ,出发，沿着使I 值下降的方向逐步地下降I 值，直到降到它的零极小值。

该方法的优点：① 方法简单、直观，计算存储量小；② 开始收敛速度快，对任意初值均可收敛。

缺点：越接近极值点收敛越慢，对于条件数大的方程情况尤其严重。

通常与其他方法组合使用，如与牛顿—拉夫逊法组合使用，先用最速下降法，后用牛顿—拉夫逊法。

2、共轭梯度法（CG ）
该方法的第一步用最速下降法，即用目标函数在()0A 处的负梯度作为搜索方向，其余各步均用共轭方向作为搜索方向。

而共轭方向是用泛函的梯度来构造的。

共轭方向的定义：设K 为n n ⨯阶正定对称矩阵，P 为1n ⨯阶矩阵，若满足
()()
()1-n ,0,1,2,j i, j ,i 0KP P j T
i
=≠= （4-64）
则称()i P 和()j P 为对刚度矩阵K 共轭，且()i P 和()j P 为线性独立。

上式表明：第j 次的搜索方向()j P 经过K 矩阵变化后才与()i P 正交，这是与最速下降法区别之处。

基本步骤为：
（1）设初值()0A ，计算()()00KA F r -=，令第一次搜索方向()()00r P = （实际上()()
()000KA F A I r -=-∇=)(，一开始余量不为零，这是最速下降法）
（2）如果
()()
()ε<k T
k
r r ，则()k A 为解，程序停止。

否则按（3）进行
（3）按下面计算
第k 次步长因子 ()
()()()()()()k
k
k
T
k
k
KP
P r r =α k =0，1，2，…… 计算新点 ()()()()k k k 1k P A A α+=+ 第k +1次梯度 ()()()()k k k 1k KP r r α-=+ 计算系数 ()
()()()()()()()()
k
T
k
1k T 1k k r r r r ++=β
第k +1次搜索方向 ()()()()k k 1k 1k P r P β+=++
1k k += 转入步骤（2）。

该方法优点：①只需存储非零元素，节约内存；②从理论上讲，在无舍入误差的前提下，对N 个变量迭代N 次后必然收敛。

缺点：计算量大，且对条件数大的方程难以收敛。

3、预处理共轭梯度法（ICCG ）
ICCG 法是将刚度矩阵做不完全分解，得到乔列斯基的T LDL 阵作为预处理阵（L 为K 的下三角阵，D 为K 的对角阵），做这样处理的系数阵可以降低条件数，将处理后的系数阵放入CG 法的迭代公式中去，令()
K LDL 1
T -替代CG 法中的K ，
就形成了ICCG 法。

4.4.5 磁性参数的提取
从有限元方程中得到磁矢位后，可以通过下式得到磁感应强度
平行平面场 ()
∑=-∆
=3
1
q z e
q e q e q e q e
e
A b A c
21e e B φρφ （4-65）
轴对称场 ()
∑=-∆==3
1
q z e
q e q e q e
q e
e e e e
A b A c
21
'
e 'e 'B B φρφρ
ρ （4-66）
求出其模值e B ，通过查磁化曲线，或通过插值函数得到e H ，单元的磁阻率
e e e H B v =，若为轴对称场，e e e /v 'v ρ= §4-5 等参元有限元法
在有限元法中，形状函数有两个作用：
（1）用于未知函数的展开，即已知单元节点上的位函数值，采用形状函数便可计算单元内任一点的位函数值。

称之为形状函数1。

（2）用于局部坐标与整体坐标之间的变换，即将整体坐标系中形状复杂的单元通过形状函数（与（1）中的形状函数的概念不同）变换到局部坐标系中的标准单元，称之为形状函数2。

这样，单元刚度矩阵的各元素可以在局部坐标系
下，用前面论述的方法建立e
ij K 。

当形状函数1用于作用（1）时，根据计算精度的要求，取低阶或高阶插值的形状函数。

当形状函数2用于作用（2）时，整体坐标中单元的节点数可以与局部坐标
中的单元节点数不同，或者说两者的插值函数的阶数可以不同：
形状函数1的阶数与形状函数2的阶数相同时，称为等参元有限元法，本节介绍内容；
形状函数1的阶数小于形状函数2的阶数时，称为亚参元有限元法； 形状函数1的阶数大于形状函数2的阶数时，称为超参元有限元法。

4.5.1坐标系
整体坐标系：表示场域中所有节点的位置； 局部坐标系：表示单元节点在单元中的位置，在单元之外无任何意义。

局部坐标的表示方法是用面积比来表示的，故又称为面积坐标系。

以三角形单元为例， P 点是三
图4-3 面积坐标定义 角形单元中任一点，定义面积坐标：
∆∆=∆∆=∆∆=/L /L /L 332211 （6-67）
其中，∆是三角形单元面积，如图4-3所示。

面积坐标特点：
（1）面积坐标均小于等于1。

（2）平行于棱边12的任一直线，如直线
1'2'上的任意节点1P 都具有相同的3L 值，因为
三角形3∆的底和高相同；同理，对直线''''21上的任意节点2P 都具有相同的3L 值。

如图4-4
所示。

图4-4 面积坐标系
（3）将图4-1中的P 点拉到节点1，即可得到节点1的面积坐标，显然此时2∆和3∆都等于零，而1∆与∆相等，同理对节点2和节点3，因此有
节点1的面积坐标 0L 0L 1L 321=== （4-68） 节点2的面积坐标 0L 1L 0L 321=== （4-69） 节点3的面积坐标 1L 0L 0L 321=== （4-70） 三角形重心（中线的交点）的面积坐标 3
1
L L L 321=
== （4-71）
可见，任一节点对应的面积坐标有3个，但其中只有两个是独立的，因为
1L L L 321=++ （4-72）
所以只要用任意两个即可，即任一节点对应
的面积坐标有2个。

建立实际单元局部坐标系
''-ξη如图4-5所示，指定0L 3=的直线为'ξ轴，显然，当0L 3≠的那些直线恰是'η等于定值的直线。

指定0L 2=的直线为'η轴，显然，当0L 2≠的
那些直线恰是'ξ等于定值的直线。

因此有
132L ''1 L ' L '=--==ξηηξ 图4-5 实际单元局部坐标系
因此，在''-ξη坐标系中，只需用两个量即可确定单元中各点的位置。

将实际单元局部坐标系''-ξη转换成正交标准单元局部坐标系ξη-，如图4-5所示。

两坐标之间的关系为：
' 'ηηξξ==
在标准单元局部坐标系ξη-中，直角边边长均为 1。

对于平面上的曲边三角形，如图4-7所示。

图4-6 标准ζ-η 坐标系 图4-8所示6节点曲边三角形的标准单元坐标系，ξη-如图4-9所示。

图4-7 曲边三角形 图4-8 六节点曲边三角形局部坐标系
图4-9 标准ζ-η 坐标系
图4-10 四边形局部坐标系
图4-11 标准ζ-η 坐标系
4.5.2 整体坐标与局部坐标之间的转换关系
在图4-2中，1∆的面积为
()()()[]2332233233221x x y y y x y x y x 21
y x 1y x 1y
x
1
2
1
-+-+-==
∆
用前面定义过的符号可以看出：
[]y c x b a 2
1
1111++=
∆ 将其代入面积坐标可以得到局部坐标（归一化后的坐标均小于1）
[][][]y c x b a 21
L y c x b a 21
L y c x b a 21
L 333332222211111++∆
=
∆∆=
++∆
=
∆∆=++∆
=
∆∆= （4-73） 用矩阵形式表示为
⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡∆=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡y x 1c b a c b a c b a 21L L L 33
3
222111
321 （4-74）
反之，也可得到用面积坐标表示整体坐标的关系：
y
L y L y L y x L x L x L x 322113
32211++=++= （4-75）
比较式（4-71）～（4-73）、与前面讲的形状函数可知
e
3
3e 22e 11N L N L N L === 所以式（4-73）可以写成
∑∑===++==++=3
1
i i
e i 3e
32e 21e 13
1i i
e i 3e 3
2e 2
1
e
1y N y N y N y N y x N x N x N x N x （4-74）
可以证明，整体坐标与局部坐标变换后，①单元顶点的坐标经变换后保持不变；②单元之间不会重叠，也不会分离；③分界处保持连续。

4.5.3 等参元有限元法
根据前面所述已知
⎰=⎪⎪⎭
⎫ ⎝⎛+∂∂∂∂+∂∂∂∂=e s e j e i e j e i y e j e i x e ij 1,2,3j i, dxdy N N y N y N x N x N K βαα
1,2,3i
dxdy fN F e S
e i e i ==⎰⎰ 由于e i N 是ηζ,的函数，而()()ηζηζ, ,,y y x x ==，反推：由于
ζζζ∂∂∂∂+
∂∂∂∂=∂∂y y N x x N N e i e i e i η
ηη∂∂∂∂+
∂∂∂∂=∂∂y
y N x x N N e i e i e i 写成矩阵形式
⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡∂∂∂∂=⎥⎥⎥⎥
⎦⎤
⎢⎢⎢
⎢⎣⎡∂∂∂∂⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂y N x N J y N x N y x y x N N e i e i e i e i e i e i ηη
ζζηζ J 成为雅可比矩阵，它表示了坐标之间的转换关系。

其中各元素等于
i i e i x N x ∑=∂∂=∂∂31ζζ i i e i x N x ∑=∂∂=∂∂31ηη i i e i y N y ∑=∂∂=∂∂31ζζ i i e i y N y ∑=∂∂=∂∂3
1η
η 由此可得
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡∂∂∂∂=⎥⎥
⎥⎥⎦⎤
⎢⎢⎢⎢⎣⎡∂∂∂∂-ηζe i e i e i e i N N J y N x N 1 其中
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡∂∂∂∂-∂∂-
∂∂=-ζη
ζηx x y y
J J 11
η
ηζζ∂∂∂∂∂∂∂∂=y x y
x
J 可以证明，ηζd d J ds dxdy ==。

这样
3211
10
,, =⎪⎪⎭
⎫ ⎝⎛+∂∂∂∂+∂∂∂∂=⎰
⎰
-i d d J N βN y N y N αx N x N αK e j
e i e j e i y e j e i x e ij
ηζη
()3211010
,,
==⎰
⎰-i d d J fN F e
i
e i ηζη
最后得到有限与方程
F KA =
等参元法的优点：①经过坐标变换后，单元在局部坐标中非常“标准”（等腰直角三角形、正方形），容易写出形状函数的解析表达式。

②对于边界无法用折线或用折线误差大时，可以用这种方法来提高计算精度。

缺点：需要求雅可比矩阵的逆阵。

§4-6 二维涡流场
导体放入时变场中将会感应出涡流（交变磁场—感应电场—导体中的涡流—产生去磁磁场，导致集肤效应）。

稳态涡流场可以用复数形式表示，只需离散空间变量，有限元法建立方程与直流情况相同，但解方程要用特殊的解复数方程的方法求解。

瞬态涡流场用有限元法离散空间变量，用差分法离散时间变量，解方程与直流情况相同。

4.6.1 平行平面场问题
边值问题 ⎪⎪
⎪⎪
⎩
⎪⎪⎪⎪⎨⎧⎪⎪⎭⎫ ⎝⎛=⎪⎭⎫ ⎝⎛+∂∂=∂∂=-=-⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂ΓΓΓq A n A H n A 1A A J A j y A 1y x A 1x 221
t 0
S λαμωγμμ
（4-76） 约束条件 0e =⋅∇J （ 涡流是闭合线）
e J 是涡流，S J 是激励源电流，通常源电流区可以不计涡流引起的集肤效应（如多匝线圈组成的源区），只存在源电流密度。

而涡流区不存在源电流，只有感应电流。

因此
S J y A μy x A μx -=⎪⎪⎭
⎫ ⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂11 源电流区 e J y A μy x A μx -=⎪⎪⎭
⎫ ⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂11 涡流区 E J γ=s ，而e J 是未知的，可用A 表示：t
e ∂∂-=A
J γ
，说明如下： 根据Maxwell 方程可以得到感应电场
ϕ∇-∂∂-
=t
A
E 其中ϕ是标量电位，可以看作静电场中定义的标量电位在涡流场情况下的推广。

（如何消去ϕ？使场量e J 只和A 发生关系）
在平行平面场中，由于ϕ沿z 轴方向无变化，即()y ,x ϕϕ=， 所以
()0z
z
=∂∂=
∇ϕ
ϕ 而感应电场只有z 方向分量z z E e E =，所以 ()t
t t z ∂∂-=∇-∂∂-=∇-∂∂-=A A A E ϕϕ 因此，微分方程可以合并写为
S J A j ωωy A μy x A μx -=-⎪⎪⎭
⎫ ⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂11 γ和s J 通常按区域分别给出，一般不会同时出现在同一个区域。

求出z z A e A =后，
()y y x x z z B B A e e e A B +=⨯∇=⨯∇=
t
e ∂∂-=A J γ
注意：三维场就不能实现这样的简化，场量将与A 、ϕ同时发生关系。

运用式（4-3）前先作如下代换（以下A A
→ ） t s H q ,0 ,J f ,j ,1
,A ===→→→λωγβμ
αϕ
式（4-76）对应的泛函极值问题
()⎪⎪
⎩⎪⎪⎨⎧==--+⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂=⎰⎰⎰⎰⎰⎰⎰02
221
22121A A Adxdy J Ad ΓH dxdy A j ωdxdy y A x A μA I ΓS s Γt S S min γ （4-77）
式中的第一部分和第四部分离散结果与前述相同，分别为
{}[]()e
e
T
e
e 1A K A 2
1I =
和 {}{}e
T
e
e 4
A F I
=
其中， ()⎰⎰+∆=⎪⎪⎭
⎫ ⎝⎛∂∂∂∂+
∂∂∂∂=e S
j i j i e e e j e i e j e i e
ij
c c b b 41dxdy y N y N x N x N 1K μμ 1,2,3j i, J 3
F e
S e i =∆=
式中第二部分表示消耗在场域中的热能，它由涡流产生
⎰⎰
⎰⎰
=
=
e
e
S T S 2e 2Adxdy A 2
j dxdy A 2
j I ωγωγ
将[]{}
e e A N A =代入上式
[]{}{}[]{}{}{}[]{}e
e
T
e S
e
e
T
e
e
e
2A T A 2
1
dxdy A N A N 2j I e
==
⎰⎰ωγ
式中
[][][]{}
d x d y N N N N N N N N N N N N N N N N N N j d x d y N N j T
e e
S
e 3e 3e 2e 3e 1
e 3e
3e 2e 2e 2e 1e 2e
3e 1e 2e 1e 1e 1S e T
e e
⎰⎰⎰⎰⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==ωγωγ 应用积分
1,2,3
i 61
dxdy N N e S e
e
i e i =∆=
⎰⎰ 1,2,3i j i 121
dxdy N N e S e e j e i =≠∆
=⎰⎰
得到
[]
⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡∆=2 1 11 2 1 1 1 212j T e
ωγ 其中
()1,2,3j i, 112
j T T ij e
ji
e ij =+∆
==δωγ 这是虚数，反映损耗能量的参量。

式中第三部分由第二类边界条件产生
[]{}{}{}e
T
e
e e t e e t e e 3A B A H l H l 0I ==
对{}A 求极值后为{}B ，e l 是边界单元的边界线段。

0e 个单元叠加后，全场域有
4321I I I I +=+
变分后 []{}[]{}{}{}{}P B F A T A K =+=+
变分问题 [][](){}{}⎪⎩⎪⎨⎧==+Γ0A A P A T K 1
这是复数方程。

解方程后获得复磁矢位。

各剖分单元重心上的涡流密度为
e c 3
1i i
e
e
e A j 3
A j E J ωγωγγ-=-==∑
= ∑=∆
=∂∂=3
1
i e i
i e
e e
x A
c 21y A B
∑=∆
-=∂∂-=3
1
i e i
i
e
e e y
A
b 21x A B
式中，()
e
3e 2e 1e c A A A 3
1A ++=。

指出：①边界上0A 取值不同，不会影响磁密B 的大小，但会影响涡流的大小。

若满足约束条件 0e =⋅∇J ，或 0d L
=⋅⎰l J ，表示0A 取值正确。

否则调整0
A ，
重新计算。

②画场分布时，应使用各量的瞬时值，而不是向量的幅值。