电磁场数值计算.docx
《电磁场数值计算.docx》由会员分享，可在线阅读，更多相关《电磁场数值计算.docx（21页珍藏版）》请在冰豆网上搜索。
电磁场数值计算
第九章边界元法
有限元法的优点很多，但也有不足，例如：
①处理开域问题是人工边界会带来误差；②通过位函数求场量，计算精度受影响；③求解瞬态问题计算量大，有时产生振荡，不收敛。
而边界元法恰好能弥补有限元的不足。
积分方程法分为：
①体积分方程法；②边界积分方程法，只需要在场域边界进行数值化离散，又称边界元法。
边界元法的优点：
能较好的处理开域问
题；降维降阶，存储单元少。
但方程的系数矩阵中各元素要用数值积分得到，且为满阵；多种媒质分界面需要专门处理。
因此许多问题都采用有限元一边界元耦合方法，充分利用两者的优点。
积分方程的建立有直接法和间接法两种，直接法是从物理问题出发，直接导出积分方程。
间接法是从描述物理问题的微分方程出发导出积分方程。
本章介绍间接法。
间接法可以通过格林公式，将微分方程表示的边值问题变换为一个包括边界条件的边界积分方程。
也可以利用加权余量法的思想，引入权函数，再利用格林第二定理，将微分方程归结为边界上的积分方程。
然后，借助于有限元的思想，将边界离散化化，引入插值函数，得到线性代数方程组。
前者有严格的数学推导，论据充分；后者简单易行，不但对自伴算子适用，对复杂的非自伴算子也适用。
§9-1边界积分方程
9.1.1基本解的概念
边界积分方程利用格林公式从微分方程边值问题导出。
电磁场边值问题：
L=-f
盾|
=qs（9-1）
cnls2
L.SiS
若L为拉普拉斯算子“I2”，则对应静电场、恒定电场和恒定磁场。
若L为亥
姆霍兹算子“I2•k2”，“'2k”，则对应于正弦交流稳态场问题。
满足方程L=-人r-r（9-2）
的解:
*r,r称为对应于L二—f的基本解。
可见，基本解具有单位点源的性质。
将边界上连续的点源看成离散的分布点源的极限，把它们产生的场叠加后，就是问题的解。
如：
以无限远为参考点，体电荷产生的电位
对于静（恒定）电场、磁场问题的基本解（也称为格林函数）为
书和••是两个调和函数，令
（格林函数，即基本解G二「*）
「GGfdV
设在边界E上—-q，则有
这是边界积分方程的原始公式。
下面以三维泊松方程为例，推导边界积分方程。
由式（9-2）可知
'2G…r—r
代入式（9-8），再根据「•函数的性质，当场点与源点重合时，有
2
:
i-7Gdv--<■r-躲dv--i
式中，\是V域内集中点电荷所在处结点i的「值。
式（9-8）可以写为
可见，只要求出边界上的物理量：
「（在边界S>上）和q
上），即已知边界上的等效（电荷、电流）源分布以及
参考位，便可得到V域内任一点的物理量\0
当结点i在边界上时，场点与源点重合（r=rJ
处，面积分会出现分母为零的奇异积分，处理方法如
二一（在边界S
下:
设S面光滑，以边界S上场点i为球心，半径r。
二；做半球图9-1求解边
界上的物理量面S，如图9-1，令；》0求相应面积分在点i的极限值。
这时,
半圆球面S•中结点i仍作为内点，故式（9-9）仍适用。
下面分别讨论：
（1）场点i位于面上（第一类边界条件）
当；》0时，&>0，式（9-10）第一项为
:
G
（9-11）
将基本解
代入式（9-10）中第二项，并且申ss—=%，得到
二一旦2边
4兀e4兀ew丿$Pi4兀4兀
Q是场点i对S；所张的立体角
（2）场点i位于S2面上
同上推导,
（3）场点位于V域外
此时可12G=0式6（r_r'），故J申》2GdV=0，有
从式（9-9）、（9-13）、（9-14）、（9-15）可得边界积分方程通式
vfGdV
：
：
G
Gqsds
「n
（9-16）
其中
Ci
Q
4
0
场点在V域内
场点在SSiS2上
场点在V域外
G—VfGdV-..s:
—-Gqds-n
（9-17）
若为二维问题，同上推导，但基本解和立体角不同：
11
维场基本解申=G=丄In-^，所张平面角Ot日
2r—r]
代入式（9-16），得到
佻&Gq杆-叽閃
_Gqsdr
（9-18）
式中
「sfGds
「守一Gqdr
场点在s域内
场点在S域外
若边界光滑，心二，则Ci
解的步骤：
（1）用G=0.5的公式建立方程，每一系数均要积分;
（2）解方程，得到S上的一值，以及S2上的「值；
cn
⑶用Ci=1的公式，用已知的s上的一值，以及S2上的「值
En
求场域中任一点的值。
9.1.3用加权余量法推导边界积分方程
设定解问题
（1）的近似解为~，它是某一线性无关的完备函数集合
n
〜八Gi
id
（9-19）
那么，余量可以表示为
R八2〜f
在V域内
R十〜-s
在边界s上
（9-20）
c~〜
Rs2-一-qs=q-qs
在边界s＞上
.:
n
引入权函数W，使之在平均意义上令余量的加权积分为零。
根据误差分布
原理，可以得到
W
vRWdvsRs2Wds,Rsids二0（9-21）
vs2si■n
式（9-21）左边可以表示为
vRWdv=v、2「fWdv
利用格林第二定理，式（9-21）左边还可以表示为
用式（9-23）代入式（9-22），得到
§9-2边界积分方程的近似求解
9.2.1边界元素
二维场域d的边界r是一维曲线，将其离散为N个边界元r,】2,……，-N规定：
单元序号（或结点序号）与边界r的走向一致，
D
所示。
插值元素（函数）可分为定常元素、线性元素、混合元素和高阶元素
9.2.2积分方程离散化
、定常元素
图9-3定常单元
定常元素是指每一单元上的:
和q均为常数，且等于单元中点处的值，如图9-3所示。
这时，单元中点（不是单元端点）作为结点，边界积分方程（设f=0）为
Ci：
i「SfGds-
£G
-:
n
将边界离散后，式（9-25）写为
N
由于在单元上的「和q均为常数（包括未知数），故可将其提出积分号
NFGN
（9-27）
Gi、：
j—d八qjGd
j丄-jnj_i-j
积分与结点i和单元j有关，即积分项中被积函数有r「rj项，令
Hj-d,Gj二Gd（9-28）
ljcnTj
Hj,Gj一般由数值积分求解，特别简单的也有解析解，式（9-27）可以写为
NN
CijHij八qj—ij（9-29）
jj二
当i=j，场点与源点重合时，Ci=0.5;当场点在边界上，但i=j，G=0；当场点在D域内，i=j，Ci=1。
故上式可以写为
Hij"j
Hj0.5i=j
（9-30）
式中
其矩阵形式
（9-31）
由于-1上有N1个q未知，丨2上有N2个'未知，M+N2=N，重新整理上
式，有
A妝11（9-32）
F为已知边界条件，X二hq……,qN1,;i,2……\2T，解方程（9-32），得到q
场中任一点的（i点在D域内）为
-GNN
—.G1d-「'Gd八qjGjjHj（9-33）
1cn1cnpj=1
场中任一点的场强（上式中只有G与结点i有关：
r^rj项，并且边界上的「与q均已知和已求得）为
-E
二
:
x
—qix
；：
G:
:
；x：
：
n
d：
-E
■y
=qiy
d：
空d,
一n_x
J
my
q和「都已知，G是
上式积分仍然要经过离散化处理。
由于边界上所有节点的基本解，是x，y的已知函数，进行离散时与前面是同阶的，故「与」—具
ex£y
有相同的精度，这是边界元法的特点。
二、线性元素
3
q和「呈
（9-35）
（9-36）
-1
图9-4线性单元
位置按同次变化，如q是电位的导数，它的插值函数的次数最好比「低一次,
这就是混合元素）
边界离散为N个单元，单元的端点作为元素的结点，每一元素内线性变化。
设=ab
=■时，「1;—「时，=1，解得
a」「,b」「j—j
22
代入式（9-35），得到
4―;：
j11「j=2^j2vj1
矩阵形式
（9-38）
q有同样的插值函数
代入边界积分方程的离散式（9-27）
式中
还可以写为
GiNf'q1q2qN
i=1,2,3N
（9-41）
CiiHi1Hi2……HiN“1「2……:
CGi1Gi2
式中Hij=hj■hij4,Gij-gijgijJ
（说明式（9-41）:
当j=2时，hl笃-hiV：
3
当j=3时，hi化+h泸4
当j=4时，讥+识
乞h1h%卜……hi]笃+（hi]+hi笄3+（hl+hi弈4十……）j丘'j1
上式可以简化为
、Hj<八Gjqj
j4j=1
当结点i不属于单元j时
后面，与定常单元得到的矩阵方程相同的处理，可得到"上的q,丨2上的
（第十七次课）
§9—3亥姆霍兹（Helmholtz）方程求解（不讲）
9.3.1用加权余量法推导边界积分方程
在均匀、各向同性媒质中的波动问题可用下列方程组表示
式中'二E表示电场问题，'二H表示磁场问题。
根据用加权余量法推导的边界积分方程式（9-21）
在二维场中，上式为
边界光滑时，当结点i与单元j任一端点重合时
取权函数为基本解，即W=G，利用下式
丁—可2®ds=q畀电dF—JjWNZGds
S-「nS-「nS
（9-46）
二人Ng-k^Gds=*G计d「—LLG®2®+k2®Gds
代入式（9-45）后，得到
-JG」-k2Gds「..GqSd「G「小-..亠d.*空d】
S-21〔n-1：
「n-1^n
再利用（9-46）
:
G
一SG「ds—「点小J2Gds
得到
s'2Gk2G：
ds=-「Gqsd-】Gqd厂亠!
厂门d「亠卜込d
由于基本解满足下式
根据7函数的性质，代入上式得到
i-,GqsdJ二」.八石山
整理后，与拉普拉斯方程同样的推导，得到类似的方程
可见，表达式与拉普拉斯对应的公式（9-33）相同。
9.3.2代数方程组的形成
用线性元素，离散化和代数方程的推导与上述相同，只是H1,G1中各元素的积分不同而已。
9.3.3波动方程的解在无穷远处的性质
边界元法能处理开放域问题，这时，一部分区域在无穷远处，边界由丨和
、构成，如图9-5所示。
这时的边界积分方程为
图9-5无限远边界
在边界-1和-2上的积分，由于法线方向相反，故积分之和为零。
所以
Ci%=j/g兰-98h+dG兰/皂h
1Vcnon丿0cn的丿
在，：
上可认为是半径趋于二的球面，认为：
：
=0。
可以证明，第二项积分为零，即无穷远边界可以不考虑，那
；:
n：
：
n
（9-48）
可以证明，无论是二维场还是三维场，都只需要考虑有限边界上的积分§9—4瞬态涡流冋题
9.4.1边界积分方程
低频瞬态涡流场A二Ar,t边值问题为
2亠
、、A0
ct
在边界s上
（9-49）
当1=0时，在区域V内
这是二阶抛物型方程（扩散方程），与椭圆型方程不一样，这类方程如用有限差分或有限元法，需要很小的时间步长，否则数值解会产生严重的振荡，不易收敛，费机时。
用边界元法解瞬态涡流问题，可以采用时变基本解，因而时间步长可以取的比较大。
如果在分析中始终用同一个时间步长，则系数阵只需计算一次，因为系数阵只与■:
t有关，与t无关）。
另外，若边界条件不随时间变化，则矩阵只需求逆一次，省时。
时变基本解满足方程
G」旦…r-r7-t（
ti,t是两个任意时刻，且0乞t^tj。
（微分都是对第一个变量进行，如：
r,ti
对于二维问题，其基本解为
r-r
二tn-t
exp
（9-51）
若是三维问题d=3,二维问题d=2。
可以看出，基本解表明在时刻t，空间r处的单位强度瞬变点源，经过时间tn-1后，在时刻tn，空间r处的响应
根据格林第二定理
（9-52）
口GA—ANGds=G二一A
S八內
对时间积分
将式（9-54）代入式（9-53），将式（9-49）代入式（9-53）
左边
二—sGAyds-JsGAt乂ds-
第二项是t=0时GA的初值。
第一项GAt土是t=tn时的GA值，从式（9-51）
可以看出，当t=tn时，只有r=r，积分才有效，否则是无穷大，第一项积分
后得到」Art土，所以式（9-53）为
9.4.2数值离散化过程
对于式（9-55）的数值解，必须在时间和空间域中分别进行离散。
将边界『离散为N个单元，空间插值函数用••表示。
内部区域S剖分为若干子域，空间插值函数用.表示。
时间步长•进二ti-ty，时间插值函数用恒定时间和线性时间两种。
以恒定时间为例,
即函数A及兰值在每一时间段中为恒定值，方程（9-55）
cn
可以离散为
HAn=GQnJPnBAn」
（9-56）
Gdtd】i、ij
丄：
：
n
式中，时间积分均可解析解出
0：
：
x：
：
'
ti
tGtn,tdt二——Eian-ti丄
2-R
利用式（9-56）可以求出第i次时间步长时边界上的A及值，并以此为下一步
占n
的初始值。
（1）式中的H、G、J、B均不变的前提是分析过程中只用同一个：
t，这些矩阵只需计算一次。
（2）式（9-56）可整理为
UIx^=怡｝
式中'X■表示第i次时间步长时边界结点待求的位函数值和法向导数值。
由于在氏不变时，A1在各个步长中不变，仅需求逆一次，即
因此，每个时间步长内，只要计算d，然后与Ai-1相乘即可。
1
若边界光滑，亠2二，则G二丄。
也可写为
2
升级会员 