最近、オープンソースでFDTD法を用いて電磁界を計算してくれるopenEMSというものがあることを知った。openEMSは、OctaveまたはMatlabから利用できるようになっている。とりあえずインストールはできたので、使い方を覚えるためにチュートリアル Tutorials - openEMS をひとつずつ試していこうと思う。ソースは openEMS/matlab/Tutorials at master · thliebig/openEMS · GitHub からダウンロードできる。
目次
Tutorial: Parallel Plate Waveguide
まずは最初のチュートリアル「Parallel Plate Waveguide」から見ていこう。
Octaveスクリプト
OctaveのGUI。Octaveはスクリプトをファイルに書いて実行することもできるが、コマンドウィドウからひとつひとつ打ち込んで実行することもできる。
最初に、下のコマンドで実行環境をまっしろにしておく。
close all
clear
clc
まず、FDTDオブジェクトを作る。
FDTD = InitFDTD('NrTS',100,'EndCriteria',0,'OverSampling',50); FDTD = SetSinusExcite(FDTD,10e6); FDTD = SetBoundaryCond(FDTD,{'PMC' 'PMC' 'PEC' 'PEC' 'MUR' 'MUR'});
1行目でFDTDオブジェクトを初期化して、時間ステップの数(NrTS: Number of time steps)を100、計算終了の基準値(EndCriteria)を0(デフォルトは1e-5)、オーバーサンプリングを50に設定している。
2行目では、正弦波の波源(Excitation)の周波数を10 MHzとしている。
3行目では境界条件を { xmin xmax ymin ymax zmin zmax } の順で記述している。今回は「Parallel Plate Waveguide」なので、yminとymaxに完全導体(PEC: perfect electric conductor)を割り当てている。また、z方向に伝搬する電磁波を吸収するため、zminとzmaxにはMUR吸収境界を割り当てている。
次に、3Dデータを扱うためのCSXCADライブラリを初期化する。
CSX = InitCSX();
電磁界を数値解析するために領域をメッシュで分割する。デフォルトの単位は m(メートル)。ここでは、xy方向は-10 mから10 mまで、z方向は-10 mから30 mまでを1 m刻みで区切っている。メッシュの間隔は、解析する波長より十分小さい方が良い。
mesh.x = -10:10; mesh.y = -10:10; mesh.z = -10:30; CSX = DefineRectGrid(CSX, 1, mesh);
波源を追加する。
CSX = AddExcitation(CSX,'excitation',0,[0 1 0]); CSX = AddBox(CSX,'excitation',0,[-10 -10 0],[10 10 0]);
AddExcitation関数の第3引数は波源の種類を指定する。0だと「E-field soft excitation」という意味らしい。第4引数の [0 1 0] は、y方向に1 V/mの電界強度で励振するということのようだ。
AddBox関数は直方体を追加するための関数で、'excitation'という名前のプロパティに始点 [-10 -10 0] から終点 [10 10 0] までの長方形を割り当てている。
同じ要領で、電界を観測するためのプレーンを設定する。
CSX = AddDump(CSX,'Et'); CSX = AddBox(CSX,'Et',0,[-10 0 -10],[10 0 30]);
最後に、ここまでで作ったFDTDオブジェクトとCSXオブジェクトを、1つのxmlファイルに書き込む。
mkdir('tmp'); WriteOpenEMS('tmp/tmp.xml',FDTD,CSX);
3DデータをAppCSXCADで表示
Octaveスクリプトで定義した3D構造を可視化して確認することができる。
CSXGeomPlot('tmp/tmp.xml');
この関数を実行すると、AppCSXCADのウィンドウが自動的に開かれる。マウスでぐるぐる回転させたりできて便利だ。
OpenEMS シミュレーション
xmlファイルができたら、次のコマンドでopenEMSによる解析が始まる。時間ステップ数を100としているので、Et_0000000xxx.vtrみたいな名前のファイルが101個出力される。
RunOpenEMS('tmp','tmp.xml','');
結果をParaViewで表示
出力されたEt..vtrファイルはParaViewという別のソフトを使って表示できる。File->OpenからEt..vtrを選んで、左のPropertyタブにあるApplyをクリックして、ドロップダウンからE-Fieldを選択すればそれっぽい画像が表示される。再生ボタン▶を押すとアニメーションになって電磁波が伝搬する様子が見られる。