电磁场数值方法1

h4
h1
4
φ4 (i, j −1)
(i +1, j −1)
i增 加 的 方 向
⎪⎧⎪(
∂φ
∂x
⎨ ⎪⎩⎪(
∂φ
∂x
)0 )0
≈ ≈
φ1 − φ0
h1
φ0 − φ3
h3
+ O(h2) + O(h2)
×
使差分格式误差最小
待定系数的方法
α
(φ1
−
φ0
)
+
β
(φ3
−
φ0
)
=
(
∂φ
∂x
)0
(α
h1
−
β
h3
)
2
h2
3
A有源
0
1
B 无源
h3
ε1
h4
ε2
4
h1
φa1 + φa2 + φa3 + φa4 − 4φa0 + h2Wa = 0 φb1 + φb2 + φb3 + φb4 − 4φb0 = 0
φa2 = φb2 = φ2 ,φa4 = φb4 = φ4 ,φa0 = φb0 = φ0
ε1
(
∂φ1
∂n
(1-7)
对偏导数，可仿照上述方法，将 ∂u 表示为： ∂x
∂u ≈ u(x + h, y, z) − u(x, y, z) （1-8）
∂x
h
同样，二阶偏导数可表示为：
∂ 2u ≈ u(x + h, y, z) − 2u(x, y, z) + u(x − h, y, z)
∂x 2
h2
（1-9）
§2 差分方法的求解步骤
有限差分（FDM）
在电磁散射计算方法中，有限差分法自上世纪 五十年代以来得到了广泛的应用，该方法概念清晰 ，方法简单，直观。虽然其与变分法相结合所形成 的有限元法更有效，但有限差分还是以其固有特点 在数值计算中有其重要地位。
为求解由偏微分方程定解问题所构造的数 学模型，有限差分法是将定解区域（场区）离 散化为网格离散节点的集合。并以各离散点上 函数的差商来近似该点的偏导数，使待求的偏 微分方程定解问题转化为一组相应的差分方程 。根据差分方程组解出各离散点处的待求函数 值—离散解。
+
1 2!
∂2φ
∂x2
(α
h12
+
β
h32
)
+
L
α h12
+
β h32
=
0
=> α
=
−
h32 h12
β
∂φ
( ∂x
)0
≈
α (φ1
− φ2 ) α h1
+ −
β (φ3 β h3
−φ0 )
=
h32 (φ1
−φ0 ) −
h1h3 (h1
h12 (φ3
+ h3 )
−φ0 )
∂φ
( ∂x
)0
≈
φ1 − φ3
，第三类边界条件。
介质不连续处还要增加连接条件
⎧⎪φ1 |G = φ2 |G
⎨⎩⎪ε1
∂φ
∂n
|G
−ε 2
∂φ
∂n
|G
=
σ
j增 加 的 方 向
(i −1, j +1)
3
h3
φ3 (i −1, j)
(i −1, j −1)
φ2 (i, j +1) (i +1, j +1)
2
φ0
0 (i, j)
1(φi1+1, j) h2
h2
4
h3
h1
+
1 4
Rh2Wa
=
0
边界条件的处理(续)
• 5、与节点不重合的边界
– 应用不等间距差分格式
∇2φ
= 2[ h3 (φ1 − φ0 ) + h1(φ3 − φ0 ) +
h1h3 (h1 + h3 )
h4 (φ2 − φ0 ) + h2 (φ4 − φ0 ) ] =
h2h4 (h2 + h4 )
f0
−φ0 )
+ O(h3)
∇2φ
=
2[ h3 (φ1 − φ3 ) + h1(φ3 − φ0 ) +
h1h3 (h1 + h3 )
h4 (φ2 − φ0 ) + h2 (φ4 − φ0 ) ] =
h2h4 (h2 + h4 )
f
(x0 , y0 ) =
f0
h1 = h3 = hx , h2 = h4 = hy
h1 = ph, h3 = qh, h2 = h4 = h
A区
φ1
2 p(p + q)
+φ2
+φ3
2 q(p + q)
+φ4
− 2φ0 (1+
1 ) − h2 pq
f0
=
0
2
h2
3 h3
0 h1
1 B区
4h4
利用边界条件和不等间距差分格式消去虚元
φb1
2 p(p+
pR)
+φ2
+φa3
2R q(q+ pR)
+ h3 )
− φ0 )
3 h3 h2 h1 1
r0
h4
∂2φ
( ∂y2
)0
=
2
h4 (φ2
−φ0 ) +
h2h4 (h2
h2 (φ4
+ h4 )
− φ0 )
4
r
φ0
(
2 h2h4
+
2r0 + h3 − h1 ) h1h3r0
=
2 h2 (h2 +
h4 ) φ2
+
2 h4 (h2 +
h4 ) φ4
§1 差分与差商
设函数 f (x) 的自变量 x 有一小增量 Δx = h ，则 f (x) 的增量为 Δf (x) = f (x + h) − f (x) （1-1）
Δf (x) 为函数 f (x) 的一阶差分。当增量 h 足够小，差分 Δf 与微分 df 之间的差才足够小。 一阶差分 Δf 是自变量 x 的函数。按式（1-1）计算 Δf (x) 的差分 Δ2 f (x) 称二阶差分，且 Δ2 f (x) = Δf (x + h) − Δf (x) （1-2）
=
f0
1 hx2
(φi+1, j
− 2φi, j
+ φi−1, j ) +
1 hy2
(φi, j+1 − 2φi, j
+ φi, j−1) =
fi, j
i增 加 的 方 向
hx = hy = h
五点格式
φi+1, j + φi−1, j + φi, j+1 + φi, j−1 − 4φi, j = h2 fi, j
有限小差分 Δx 的商，称为差商。
一阶导数 f ' (x) 还可表示为：
向后差商
df ≈ Δf (x) = f ( x) − f (x − h)
dx Δx
h
向前差商
（1-4）
df ≈ Δf (x) = f (x + h) − f (x)
dx Δx
h
中心差商
（1-5）
df ≈ Δf ( x) = f ( x + h) − f ( x − h) （1-6）
1+ R
+ φ2
+ φa3
2R 1+ R
+ φ4
−
4φ0
=
0
▲
分界面为场域边界，满足第二类边界条件
∂φ
( ∂n )0
=
φ1 − φ3
2h
=
K
φ3 = φ1 − 2hK
B：2φb1 + φ2 + φ4 − 2hK − 4φ0 = 0
A：2φa3 + φ2 + φ4 − 2hK − 4φ0 + h2Wa 2 = 0
3
0
1
不会引起原方程变化
β < 0.5π
1，2，3，4点都在媒质B区中
Q
β μ1,
r J1
4
μ2 P
B区
Ab1 + Ab2 + Ab3 + Ab4 − 4 A0 = 0
(1)
A区 R
M
α =π
Ab1
+
Ab 2
+
R( Aa3
+
Aa4 ) −
2(1 +
R) A0
+
1 2
Rh2Wa
=
0
(2)
( 12) ×[(1) + (2)]
)0
=
ε
2
(
∂φ2
∂n
)0
=> ε1(φa1 − φa3 ) = ε 2 (φb1 − φb3 )
ε1 ε2 = R
2φb1
1+ R
+ φ2
+ φa3
2R 1+ R
+ φ4
−
4φ0
+
R 1+ R
h2Wa
=
0
▲
R = 1,
φa1
+ φ2
+ φa3
+ φ4
−
4φ0
+
1 2
h2Wa
=
0
▲
两边均无源，2φb1
+φ4
−2φ0(1+
1 )− pq
Rh2 q+
f0 p pR
=
0
边界条件的处理(续)
• 6、曲线边界的情形
– 第一类边界条件的处理
9 直接转移法 φ0 ≈ φ1
h2
9 线性插值法
3
h4
若x方向最靠近0点
φ0
=
h3φ1
h3
+ h1φ3
+ h1
A区
2
0
1
B区
若y方向最靠近0点
φ0
=
h4φ2
h4
+ +
h2φ4
+ φa4 )
φbx
=
1 2
(φ b1
+ φb2 )
φby
=
1 2
(φ b 3
+ φb4 )
y
B区
4 A区
2(φb1 +φb2)+2R(φa3 +φa4)−4(1+R)φ0 +Rh2Wa =0
边界条件的处理(续)
• 3、边界平行于网格，但有拐点
– 无法引入虚构点
L
N2
– 引入辅助线
α
– 0.5π < α < 1.5π ,0 < β < 0.5π
网格划分方式离散化场域
给出相应的差分计算格式
求解
§3 二维泊松方程和拉普拉斯方程的 有限差分法
• 差分格式的建立
∇2φ
=
∂2φ
∂x2
+
∂2φ
∂y 2
=
f
(x,
y)
φ |G = g( p)
，第一类边界条件；
∂φ ∂n
|G
=
g(
p)
，第二类边界条件；
φ
|G
+ g1( p)
∂φ ∂n
|G =
g2( p)
边界条件的处理(续)
• 2、边界不平行于网格，但是边界无拐点
– 边界旋转
q
2
p
3
1
0
B区
r
4
s
A区
φbp
2 1+ R
+φq
+φar
2R 1+ R
+φs
−4φ0
+R 1+ R
h'2 Wa
=
0
φbp
2 1+ R
+ φq
+ φar
2R 1+ R
+ φs
− 4φ0
+ 2R 1+ R
h2Wa
=
0
– 采用边界条件重新推导
Ab1
+
Ab 2
+
1 2
(1 +
R)( A3
+
A4 ) −
(3 +
R) A0
+
1 4
Rh2Wa
=
0
ϕ ~ A,ε ~ 1 μ
R = μ2 μ1
A是铁磁物质 R ≈ 0
2
铁磁物质
3
铁空磁气
μ物1,质Jr1
A区
0 空气 1
μ2
4
B区
Ab1
+
Ab 2
+
1 2
( A3
+
A4 )
−
3 A0
=
0
B是铁磁物质
R≈∞
2hx
+ O(h3)
h1 = h3 = hx
使差分格式误差最小
待定系数的方法
α
(φ1
−
φ0
)
+
β
(φ3
−
φ0
)
=
(
∂φ
∂x
)0
(α
h1
−
β
h3
)
+
1 2!
∂ 2φ
∂x2
(α
h12
+
β
h32
)
+
L
α h1
−
β h3
=
0
=> α
=
h3 h1
β
∂2φ
( ∂x2
)0
≈
2 α (φ1
−φ0 ) α h12
+ +
轴对称场差分格式
• 球、柱、回旋体等形状 （3D->2D）
柱坐标
∇2φ
=
∂2φ
∂r 2
+
1 r
∂φ
∂r
+
∂2φ
∂z 2
=
0
∂ 2φ
( ∂x2
)0
=
2
h3 (φ1
− φ0 ) +
h1h3 (h1
h1 (φ3
+ h3 )
− φ0 )
z
2
( ∂φ
∂x
)0
=
h32 (φ1
−φ0 ) −
h1h3 (h1
h12 (φ3
函数 f (x) 的一阶导数 f ' (x) 为：
lim f ' (x) = df =
Δf (x)
dx Δx ⎯⎯→ 0 Δx
应用差分， f ' (x) 可表示为
f ' (x) ≈ Δf (x) = f (x + h) − f (x) (1-3)
Δx
h
故 f ' (x) 可ቤተ መጻሕፍቲ ባይዱ示为差分 Δf (x) 除以
β (φ3 β h32
−φ0 )
=
2
h3 (φ1
−φ0 ) +
h1h3 (h1
h1 (φ3
+ h3 )
−φ0 )
+ O(h3)
∂ 2φ
( ∂x2 )0
=
φ1
−
2φ0
hx2
+ φ3
+ O(h3)
h1 = h3 = hx
∂2φ
( ∂x2
)0
=
2
h3 (φ1
−φ0 ) +
h1h3 (h1
h1 (φ3
+ h3 )
(i−1, j+1) φ2 (i, j+1) (i+1, j+1)
3
h3
φ3 (i−1, j)
2
φ0 0 (i, j)
1(φi1+1, j) h2
h4
h1
4
(i−1, j−1) φ4 (i, j−1) (i+1, j−1)
j增 加 的 方 向
φ1 − 2φ0
hx2
+ φ3
+ φ2
− 2φ0
hy2
+ φ4
+
2r0 + h3 r0h1(h1 + h3
)
φ1
+
r0
2r0 − h1 h3 (h1 + h3
)
φ3
h1 = h2 = h3 = h4 = h
r0 = (i −1)h0 ,i > 1
z
2
3 h3 h2 h1 1
r0
h4
4φ0
= φ2
+ φ4
+ [1+
1 2(i −
1)
]φ1
+ [1−
1 2(i −
1)
]φ3
？r0 = 0
1 lim( r r −>0
∂φ )
∂r
=
(∂φ )'
lim
r −>0
∂r r'
=
∂2φ
( ∂r 2 )r=0
4
r
∂ 2φ
2 ∂r2
+
∂2φ
∂z 2
=
0
★★
整个场域内点的差分格式共有两种!
φ1 = φ3
6φ0 = φ2 + φ4 + 4φ1
边界条件的处理
• 1、不同介质平面分界面的情形
dx Δx
2h
在上面三种差商形式中，中心差商的精度最高。
函数 f (x) 的二阶导数 f '' (x) 为
d2 f dx2
= 1 ( df Δx dx
− df x+Δx dx
)
x
≈
1 h
⎡ ⎣⎢
f
(x + h) − h
f
(x)
−
f
(x) − f (x − h) ⎤
h
⎦⎥
= f (x + h) − 2 f (x) + f (x − h) h2
A : φa1 + φa2 + φa3 + φa4 − 4φ0 + h2Wa = 0
B :φb1 + φb2 + φb3 + φb4 − 4φ0 = 0
2
边界：ε1(φax − φay ) = ε 2 (φbx − φby )
3
x 1
φax
=
1 2
(φ a1 + φ a 2 )
φay
=
1 2
(φ a 3
( A3 + A4 ) − 2 A0 = 0
r (J1 = 0)
边界条件的处理(续)
• 4、网格成对角线边界时的角形区域边界
– 可以用边界平行于网格，但有拐点的情形的处理方 法同样处理
2
3
0
1
A区
B区
4
Ab1
+
1 2
(1 +
R) A3
+
1 4
(3 +
R)( Ab2
+
Ab4 )
− (3 +
R) A0