PythonでFDTD法で電磁界シミュレーションできるopenEMSを使う(1)例題にあるマイクロストリップラインのノッチフィルタ（スタブ）を動かして電磁界分布を動画で見てみる。CSXCADでモデルは確認できるし、ParaViewで電磁界分布が見られる。Sパラメータも計算できる。

オープンソースの電磁界シミュレータはいろいろある。

Open-Source Electromagnetic Simulation: FDTD, FEM, MoM

今回はその中でEC-FDTE法のopenEMSをPythonから使ってみる。

まず公式サイトはこちら。

https://www.openems.de/

WindowsでPythonを使うときのインストールはこちら。

https://docs.openems.de/python/install.html#windows

Pythonのバージョンは3.10か3.11でないと動かないっぽいので3.11を仮想環境をvenvで作って使う。

以下をpipでインストールして、

openEMS-0.0.36-cp311-cp311-win_amd64.whl

CSXCAD-0.6.3-cp311-cp311-win_amd64.whl

環境変数を設定。

setx OPENEMS_INSTALL_PATH インストールしたフォルダ

これで動く。

問題はPythonで使うときのドキュメントがものすごく不親切（もともとMatlab/Octaveから使うもの）で、両方見ないと分からん…
とりあえず例題のマイクロストリップラインのノッチフィルタ（スタブ）やってみよう。
PythonとMatlabを見比べないと分からない部分が多くある。

https://docs.openems.de/python/openEMS/Tutorials/MSL_NotchFilter.html

https://wiki.openems.de/index.php/Tutorial:_Microstrip_Notch_Filter.html

まずは必要なモジュールをインポート。

import os

import numpy as np

import matplotlib.pyplot as plt

from CSXCAD  import ContinuousStructure

from openEMS import openEMS

from openEMS.physical_constants import C0

設定。ここでPathは各自保存したいところに変更。 

Sim_Path = os.path.join("C:\\Users\\tomoh\\Documents\\Python Projects\\OpenEMS", 'NotchFilter')

post_proc_only = False

unit = 1e-6 # specify everything in um

MSL_length = 50000

MSL_width = 600

substrate_thickness = 254

substrate_epr = 3.66

stub_length = 12e3

f_max = 7e9

FDTDの設定。ここで励振はガウシアンで、

• f0: the center frequency of the pulse
• fc : 20dB cutoff frequency --> bandwidth is 2*fc

また境界条件は [xmin xmax ymin ymax zmin zmax]の順で、PML-8は完全吸収層で8セル使うもの、MURはムーアの吸収境界条件、
PECは完全導体でGNDを表している。x方向がPML-8なのはポートがつく方向で、終端したいので。

FDTD = openEMS()

FDTD.SetGaussExcite( f_max/2, f_max/2 )

FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'MUR', 'MUR', 'PEC', 'MUR'] )

メッシュの基本設定。電流分布はAddDumpで設定する。

CSX = ContinuousStructure()

FDTD.SetCSX(CSX)

mesh = CSX.GetGrid()

mesh.SetDeltaUnit(unit)

resolution = C0/(f_max*np.sqrt(substrate_epr))/unit/50 # resolution of lambda/50

third_mesh = np.array([2*resolution/3, -resolution/3])/4

Jt = CSX.AddDump('Jt', dump_type=3, file_type = 0)

Jt.AddBox([-MSL_length, -MSL_length, substrate_thickness], [MSL_length, MSL_length, substrate_thickness])

でここが一番面倒なところで、メッシュを細かくする部分を手動で記述する。HFSSのようなアダプティブメッシュに慣れきっていると
これがわからないところかも。基本は1/3ルール（2D 金属のエッジでは電界が強くなるため、FDTD では正しく計算するのが困難になります。

つまり、例えばマイクロストリップライン（MSL）の場合、線路インピーダンスと波動伝播はシミュレーションと実測で多少異なります。そのため、精密計算では、金属の1/3を内側、2/3を外側に配置するメッシュラインを配置することが推奨されます。）を使う。メッシュは円筒形もあり。

mesh.AddLine('x', 0)

mesh.AddLine('x',  MSL_width/2+third_mesh)

mesh.AddLine('x', -MSL_width/2-third_mesh)

mesh.SmoothMeshLines('x', resolution/4)

mesh.AddLine('x', [-MSL_length, MSL_length])

mesh.SmoothMeshLines('x', resolution)

mesh.AddLine('y', 0)

mesh.AddLine('y',  MSL_width/2+third_mesh)

mesh.AddLine('y', -MSL_width/2-third_mesh)

mesh.SmoothMeshLines('y', resolution/4)

mesh.AddLine('y', [-15*MSL_width, 15*MSL_width+stub_length])

mesh.AddLine('y', (MSL_width/2+stub_length)+third_mesh)

mesh.SmoothMeshLines('y', resolution)

mesh.AddLine('z', np.linspace(0,substrate_thickness,5))

mesh.AddLine('z', 3000)

mesh.SmoothMeshLines('z', resolution)

基板の設定。Rogersを使うこと想定。

substrate = CSX.AddMaterial( 'RO4350B', epsilon=substrate_epr)

start = [-MSL_length, -15*MSL_width, 0]

stop  = [+MSL_length, +15*MSL_width+stub_length, substrate_thickness]

substrate.AddBox(start, stop )

マイクロストリップラインポートの設定。ポートをx=0で交差させると勝手にスルーラインができるそうだ。

port = [None, None]

pec = CSX.AddMetal( 'PEC' )

portstart = [ -MSL_length, -MSL_width/2, substrate_thickness]

portstop  = [ 0,  MSL_width/2, 0]

port[0] = FDTD.AddMSLPort( 1,  pec, portstart, portstop, 'x', 'z', excite=-1, FeedShift=10*resolution, MeasPlaneShift=MSL_length/3, priority=10)

portstart = [MSL_length, -MSL_width/2, substrate_thickness]

portstop  = [0         ,  MSL_width/2, 0]

port[1] = FDTD.AddMSLPort( 2, pec, portstart, portstop, 'x', 'z', MeasPlaneShift=MSL_length/3, priority=10 )

 スタブを追加。

start = [-MSL_width/2,  MSL_width/2, substrate_thickness]

stop  = [ MSL_width/2,  MSL_width/2+stub_length, substrate_thickness]

pec.AddBox(start, stop, priority=10 )

計算。確認のためにCSXCADでメッシュも見てみる。

if 1:  # debugging only

    CSX_file = os.path.join(Sim_Path, 'notch.xml')

    if not os.path.exists(Sim_Path):

        os.mkdir(Sim_Path)

    CSX.Write2XML(CSX_file)

    from CSXCAD import AppCSXCAD_BIN

    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))

if not post_proc_only:

    FDTD.Run(Sim_Path, cleanup=False)

でSパラメータを計算して図示。

f = np.linspace( 1e6, f_max, 1601 )

for p in port:

    p.CalcPort( Sim_Path, f, ref_impedance = 50)

s11 = port[0].uf_ref / port[0].uf_inc

s21 = port[1].uf_ref / port[0].uf_inc

plt.plot(f/1e9,20*np.log10(abs(s11)),'k-',linewidth=2 , label='$S_{11}$')

plt.grid()

plt.plot(f/1e9,20*np.log10(abs(s21)),'r--',linewidth=2 , label='$S_{21}$')

plt.legend()

plt.ylabel('S-Parameter (dB)')

plt.xlabel('frequency (GHz)')

plt.ylim([-40, 2])

plt.show()

ParaViewで電磁界分布見てみる。ファイル開いてApply、リスケールするのと、表示をSolid ColorからRot-H Fieldに変更しないと見られない。
実はここにハマって時間消費した…

なるほどもっともらしい。メッシュを切る部分さえ慣れれば何でも計算できそうだ。

次はパッチアンテナやってみよう（続く）。