【Python】FDTD 法で z 方向を進む 1 次元の電磁波解析

真空中の z 方向へ進む電磁波を FDTD (Finite Difference Time Domain) 法で解析する。前回，FDTD 法で解析するコードを MATLAB で記載したが，今回は Python で記載した。なお，Python の解析コードは標準ライブラリのみを使用している。

【MATLAB】FDTD 法で z 方向を進む 1 次元の電磁波解析

解析コード

z 方向は 0 ～ 2 000 mm の解析領域を設ける。z = 1 000 mm の位置に，広がり 100 mm のガウシアンパルスの電界を波源として与える。このパルスが真空中を伝搬していく様を FDTD 法でシミュレーションする。

ライブラリの読み込み

解析に使用する Python のライブラリをインポートする。今回は，標準ライブラリ math のみを使用した。

import math

パラメータの設定

解析に用いるパラメータを設定する。

ep=8.85418782e-12#真空の誘電率
mu=4*math.pi*1e-7#真空の透磁率
cc=1/math.sqrt(ep*mu)#光速
dd=2e-3#メッシュサイズ 2 mm
NZ=1000#z 方向の格子数
CFL=0.5#CFL
dt=CFL*dd/cc#時間刻み
NT=1500#時間ステップ

ee=dt/ep/dd
hh=dt/mu/dd
cz=(cc*dt-dd)/(cc*dt+dd)

sig=100e-3#ガウシアンパルスの広がり

初期条件の設定～ z 軸，電界，磁界の配列～

z 軸，z 方向へ進む電磁界として Ex，Ey，Hx，Hy の配列を用意する。ここで，電界 Ex にガウシアンパルスを与えている。

zz=[]
Ex=[]
Ey=[]
Hx=[]
Hy=[]

for i in range(0,NZ+1):
    zz.append(i)
    Ex.append(i)
    Ey.append(i)
    Hx.append(i)
    Hy.append(i)
    zz[i]=i*dd
    Ex[i]=math.exp(-(zz[i]-NZ/2*dd)*(zz[i]-NZ/2*dd)*2/2/sig/sig)
    Ey[i]=0
    Hx[i]=0
    Hy[i]=0

FDTD 法

FDTD 法では，電界の更新，磁界の更新を順に行う。また，解析領域の境界には，Mur の 1 次吸収境界条件を適用した。

for t in range(0,NT+1):
    # Update Ex,Ey
    oldEx0=Ex[1]
    oldEy0=Ey[1]
    oldExN=Ex[NZ-1]
    oldEyN=Ey[NZ-1]
    
    for i in range(1,NZ):
       Ex[i]=Ex[i]-ee*(Hy[i]-Hy[i-1])
       Ey[i]=Ey[i]+ee*(Hx[i]-Hx[i-1])

    # Mur's first order absorbing boundary condition
    Ex[0]=oldEx0+cz*(Ex[1]-Ex[0])
    Ey[0]=oldEy0+cz*(Ey[1]-Ey[0])
    Ex[NZ]=oldExN+cz*(Ex[NZ-1]-Ex[NZ])
    Ey[NZ]=oldEyN+cz*(Ey[NZ-1]-Ey[NZ])

    # Update Hx,Hy
    oldHx0=Hx[1]
    oldHy0=Hy[1]
    oldHxN=Hx[NZ-2]
    oldHyN=Hy[NZ-2]

    for i in range(1,NZ-1):
       Hx[i]=Hx[i]+hh*(Ey[i+1]-Ey[i])
       Hy[i]=Hy[i]-hh*(Ex[i+1]-Ex[i])
    
    # Mur's first order absorbing boundary condition
    Hx[0]=oldHx0+cz*(Hx[1]-Hx[0])
    Hy[0]=oldHy0+cz*(Hy[1]-Hy[0])
    Hx[NZ-1]=oldHxN+cz*(Hx[NZ-2]-Hx[NZ-1])
    Hy[NZ-1]=oldHyN+cz*(Hy[NZ-2]-Hy[NZ-1])

解析結果

真空中の z 方向へ進む電磁波を FDTD 法で解析した結果を示す。電界 Ex，磁界 Hy の伝搬をアニメーションで示している。

まとめ

FDTD 法で z 方向を進む 1 次元の電磁波の伝搬を解析した。今回，用いた FDTD 法はシンプルなアルゴリズムであるが，正確にプログラムで記述するのに苦労した。

次回は Numpy を使って，FDTD 法のプログラムを記載したい。

アニメーションの作成

FDTD 法で解析した電界，磁界の伝搬のアニメーションを作成する方法について述べる。

ライブラリの読み込み

アニメーションを作成するときは，matplotlib を読み込む。

import matplotlib.pyplot as plt
import matplotlib.animation as animation

アニメーションの作成

matplotlib の animation を使って，アニメーションを作成した。

fig = plt.figure()
image_list = []

for t in range(0,NT+1):
    # FDTD 法のコード参照
    if t % 10 ==0:
        print(t,'/',NT)
        img = plt.plot(zz,Ex,color="blue")
        plt.xlim(0,dd*NZ)
        plt.ylim(-0.2,1.2)
        plt.xlabel('z [mm]')
        plt.ylabel('Ex [A/m]')
        image_list.append(img)

ani = animation.ArtistAnimation(fig, image_list, interval=50)
ani.save("FDTD1D_without_Numpy.gif", writer="pillow")