Transcript 数値計算法I - 九州工業大学

数値計算法
第12回
2次元ラプラス方程式
2006年度
九州工業大学工学部電気電子工学コース用講義資料
講師：趙孟佑
[email protected]
1
2次元ラプラス方程式の例
• Laplace方程式
– 電位f(x,y)の空間分布
 2f  2f
 2 0
2
x
y
f=0
等電位線
4角形の中の
電位分布
f=10
f=0
1.0
2.0
3.0
f=0
5.0
Dirichlet境界
f=0
2
2次元ラプラス方程式の応用
• 電界計算に多用
変電施設の地下化
電力機器の小型化
3
2次元ラプラス方程式の応用
• 電力機器の絶縁設計
高圧（電位~100kV)
接地（電位0）
4
2次元ラプラス方程式の応用
• 電力機器の絶縁設計
接地導体
高圧導体
等電位線
電力用遮断器の模式図
5
2次元ラプラス方程式の応用
• 電力機器の絶縁設計
• できるだけ電圧を高くしたい
– 損失を少なくするため
• できるだけ小さくしたい
– ビルの地下等限られたスペースに収めたい
• 電界=（電圧）÷（距離）
– 高電圧と小型化は相反する
• 導体の形状を変えて、電界を小さくする。
– どういう形状がよいかは、電界計算で決める。
6
2次元ラプラス方程式の応用
• 電力機器の最適絶縁設計
1. 電界が弱くなりそうな形状をインプット
2. Laplace方程式を解いて、領域内の電位
分布を計算
3. 電位分布から最大電界値を求める
4. 何回も繰り返して理想の形を探す
7
2次元Laplace方程式の差分化
• x方向をDx,y方向をDyで分割
 2f
x 2
 2f
y 2

i, j

i, j
fi 1, j  2fi, j  fi 1, j
Dx 2
fi, j 1  2fi, j  fi, j 1
Dy
2
(i-1,j)
(i,j+1)
(i,j)
(i+1,j)
(i,j-1)
8
2次元Laplace方程式の差分化
f f
 2 0
2
x
y
2
fi1, j  2fi, j  fi1, j
Dx
2
2

は、
fi, j1  2fi, j  fi, j1
Dy
2
0
とかける。
一般化して、
ai, jfi1, j  bi, jfi, j  ci, jfi1, j  di, jfi, j1  ei, jfi, j1  0
9
2次元Laplace方程式の差分化
ai, jfi1, j  bi, jfi, j  ci, jfi1, j  di, jfi, j1  ei, jfi, j1  0
において、
 1
1
1 
1
ai, j  2 ,bi, j  2  2  2  , ci, j  2 ,
 Dx
Dx
Dy 
Dx
di, j 
1
1
,
e

i, j
Dy 2
Dy 2
点(i,j)の値を求めるのに、
自分自身も含め、近傍の
5点の値を使っている。
(i,j+1)
(i-1,j)
(i,j)
(i+1,j)
(i,j-1)
10
2次元Laplace方程式の差分化
N+1
格子点状の電位Φを
計算する。
N
境界の値はDirichlet
条件で与えられている
Dy
N-1xN-1個の値を求める
Dx
2
1
1
2
N
N+1
11
2次元Laplace方程式の差分化
i=2~N-1,j=2,N-1の(N-1)x(N-1)個の点について、(N-1)x(N-1)の連立
一次方程式がかける。
行列の形にかける
非常に大きな行列になる。
消去法では時間がかかりすぎる。
反復法を使う
12
SOR(Successive Over-Relaxation)法
加速緩和法
n番目の繰り返しの値から、n+1番目の繰り返しの
値を以下の式によって求める

f (n1)  f (n)   a'i, j f (n)  b'i, j f (n)  c'i, j f (n)  d 'i, j f (n)  e'i, j f (n)
i, j
i, j
i1, j
i, j
i1, j
i , j1
i , j1
これを、i=2~N-1,j=2,N-1までの(N-1)x(N-1)個の点に
ついて行う。
13

SOR法
境界の点の電位(Dirichlet境界条件で与えられる）は、それをその
まま使う。

f (n1)  f (n)   a'2,2 f  b'2,2 f (n)  c'2,2 f (n)  d '2,2 f  e'2,2 f (n)
2,2
2,2
1,2
2,2
3,2
2,1
2,3
3
2
1
1
2
3
14

SOR法
領域内で、電位を指定されている点は、a=0,b=1,c=0,d=0,e=0とおいて

f (n1)  f (n)   f (n)  f0
i, j
i, j
i, j

i,j
15
SOR法
より一般的には、
ai, jfi1, j  bi, jfi, j  ci, jfi1, j  di, jfi, j1  ei, jfi, j1  fi, j
という式を解くと考えればよい。その場合、n+1ステップ目の
更新方法は、
f
(n1)
i, j

 f i, j
(n)
(n)
(n)
(n)
(n)
(n)
 a 'i, j f i1,

b
'
f

c'
f

d
'
f

e'
f
i, j i , j
i, j i1, j
i, j i , j1
i, j i , j1  f 'i, j
j
16

f
(n1)
i, j

SOR法
 f i, j
(n)
(n)
(n)
(n)
(n)
(n)
 a 'i, j f i1,

b
'
f

c'
f

d
'
f

e'
f
i, j i , j
i, j i1, j
i, j i , j1
i, j i , j1  f 'i, j
j
fi,jの値が指定されていない時、
1
1
2
2
Dx
Dx
,
,b 'i, j  1, c'i, j 
a 'i, j 
 1
 1
1 
1 
2 2  2 
2 2  2 
 Dx
 Dx
Dy 
Dy 
d 'i, j 
1
Dy 2
, e'i, j 
1
Dy 2
, f 'i, j  0
 1
 1
1 
1 
2 2  2 
2 2  2 
 Dx
 Dx
Dy 
Dy 
fi,jの値が指定されている時、
a 'i, j  0,b 'i, j  1, c'i, j  0, d 'i, j  0, e'i, j  0, fi, j  f0
こうしておけば、全ての点において同じ形の式を使える。
17

f
(n1)
i, j

SOR法
 f i, j
(n)
(n)
(n)
(n)
(n)
(n)
 a 'i, j f i1,

b
'
f

c'
f

d
'
f

e'
f
i, j i , j
i, j i1, j
i, j i , j1
i, j i , j1  f 'i, j
j
ω：加速緩和係数
1<ω<2
問題によって最適なωがある。
最適なωを求めるやり方もあることはある。但し、かなり
複雑。
18

f
(n1)
i, j

SOR法
 f i, j
(n)
(n)
(n)
(n)
(n)
(n)
 a 'i, j f i1,

b
'
f

c'
f

d
'
f

e'
f
i, j i , j
i, j i1, j
i, j i , j1
i, j i , j1  f 'i, j
j
各点での誤差は
解に近づけばゼロになる。
(n1)
(n)
 i,(n1)

f

f
j
i, j
i, j
収束判定は
 i,(n1)
  min
j
max
N-1xN-1点中の誤差
の最大値
で行う
予め自分で決めた収束限界
19

ラプラス方程式の例
• Laplace方程式
– 電位f(x,y)の空間分布
 2f  2f
 2 0
2
x
y
f=0
101
f=0
f=10
(55,55)
ω=1.5
f=0
(45,45)
Dirichlet境界
1
f=0
101
20
ラプラス方程式の例
• a’i,j,b’i,j,c’i,j,d’i,j,e’i,j,f’i,jを計算する。
• 境界条件、f1,j, fN+1,j, fi,1, fi,N+1を決定
• 領域内の電位を指定された点の
a’i,j,b’i,j,c’i,j,d’i,j,e’i,j,f’i,jを計算する。
• 1回目の電位の推定値fi,j(1)を代入
• 誤差が収束するまで繰り返し。
– n+1ステップ目の誤差
 (n1)  max f (n1)  f (n)
i, j
i, j
21
ラプラス方程式の例
102
101
error
100
10-1
10-2
10-3
10-4
10-5
0
200
400
600
800
1000 1200
iteration
収束の具合
22
Laplace方程式の例
電位分布の等高線
系列100
系列89
系列78
系列67
8-10
系列56
6-8
系列45
4-6
系列34
2-4
系列23
0-2
系列1
100
89
78
67
56
45
34
23
12
1
系列12
Excelで描かせたもの
23
Laplace方程式の例
c
c
c
c
c
c
c
c
c
implicit real*8 (a-h,o-z)
parameter(nmax=102)
dimension phi(nmax,nmax),phin(nmax,nmax)
dimension a(nmax,nmax)
dimension b(nmax,nmax)
dimension c(nmax,nmax)
dimension d(nmax,nmax)
dimension e(nmax,nmax)
dimension f(nmax,nmax)
dimension error(nmax,nmax)
dimension x(nmax)
dimension y(nmax)
n:分割数
f0:中心導体の値
ermin:誤差の最小値
omega:加速緩和係数
itermax:繰り返し最大数
dx:x方向刻み
dy:y方向刻み
c
c
n=100
f0=10
ermin=1.d-4
omega=1.5
itermax=10000
dx=1
dy=1
c
c
c
x座標、y座標の定義
do i=1,n+1
x(i)=(i-1)*dx
end do
do j=1,n+1
y(j)=(j-1)*dy
end do
c
c
c
差分の係数の定義
do i=1,n+1
do j=1,n+1
a(i,j)=1.d0/4.d0
b(i,j)=-4.d0/4.d0
c(i,j)=1.d0/4.d0
d(i,j)=1.d0/4.d0
e(i,j)=1.d0/4.d0
end do
end do
Laplace方程式の右辺の定義
do i=1,n+1
do j=1,n+1
f(i,j)=0
end do
end do
24
Laplace方程式の例
c
c* 中心導体の点
c
do i=45,55
do j=45,55
a(i,j)=0
b(i,j)=-1
c(i,j)=0
e(i,j)=0
d(i,j)=0
f(i,j)=-f0
end do
end do
c
c* 境界条件
c
do j=1,n+1
phi(1,j)=0
phi(n+1,j)=0
end do
c
do i=1,n+1
phi(i,1)=0
phi(i,n+1)=0
end do
c
c* 初期値
c
do i=2,n
do j=2,n
phi(i,j)=f0
end do
end do
c
c ここから繰り返し開始
c
iter=0
10 continue
c
c
SOR法による更新
c
do i=2,n
do j=2,n
phin(i,j)=phi(i,j)
#
+omega*(a(i,j)*phi(i-1,j)
#
+b(i,j)*phi(i,j)
#
+c(i,j)*phi(i+1,j)
#
+d(i,j)*phi(i,j-1)
#
+e(i,j)*phi(i,j+1)
#
-f(i,j))
error(i,j)=abs(phin(i,j)-phi(i,j))
phi(i,j)=phin(i,j)
end do
end do
c
c
各点での誤差の計算と最大値の決定
ermax=0
do i=2,n
do j=2,n
if(error(i,j).gt.ermax)ermax=error(i,j)
end do
end do
c
c
c
収束判定
write(6,*)iter,ermax,ermin
if(ermax.lt.ermin)goto 20
if(iter.gt.itermax)goto 30
iter=iter+1
goto 10
30 continue
write(6,*)' iteraction max exceeded'
20 continue
c
c
結果出力
c このファイルをそのままエクセルで開いて出力可能
c
open(unit=20,file='laplace.lis',status='unknown')
write(20,110)x(1),(x(i),i=1,n+1)
do j=1,n+1
write(20,110)y(j),(phi(i,j),i=1,n+1)
end do
close(20)
110 format(102(1x,f12.5))
25
stop
end