openEMSを使ってみる。その２：Rectangular Waveguide

openEMS

Tutorial: Rectangular Waveguide

今回は「Rectangular Waveguide」（方形導波管）のチュートリアル（Tutorial: Rectangular Waveguide - openEMS）をやってみた。

Octaveスクリプト

最初に実行環境を空にしておく。

close all; clear; clc

物理定数（光速とか誘電率とか）を読み込んで、単位をμmとする。

physical_constants
unit = 1e-6; %drawing unit in microns

WR42の規格の方形導波管の寸法を指定する。幅10.7 mm、高さ4.3 mm、長さ50 mm。

a = 10700;   %waveguide width
b = 4300;    %waveguide height
length = 50000;

使用する周波数帯域を20 GHzから26 GHzとする。4行目では中心周波数 f_0=24 GHz の波長をμmで求めている。

f_start = 20e9;
f_0     = 24e9;
f_stop  = 26e9;
lambda0 = c0/f_0/unit;

方形導波管の伝搬モードをTE10とする。方形導波管の電磁波の伝搬の仕方は複数あり、進行方向に対して電界が垂直なものをTEモード（Transverse electric modes）、磁界が垂直なものをTMモード（Transverse magnetic modes）と呼んでいる。さらに、電界・磁界の分布によってTEmn、TMmn（m,n=0,1,2,...）がある。

TE_mode = 'TE10';

メッシュの細かさを波長あたり30セルとする。

mesh_res = lambda0./[30 30 30];


FDTDオブジェクトを初期化する。SetGaussExcite関数は、時間領域でガウス関数を包絡線とするような波源を設定することができる。第2引数は中心周波数で、第3引数は20dBカットオフ周波数を指定する。ガウス関数のフーリエ変換はガウス関数なので、周波数の幅を決めると時間波形の幅が決まるっぽい。

FDTD = InitFDTD('NrTS',1e4,'OverSampling',5);
FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start));

境界条件を設定する。z方向に伝搬する波を吸収するために、導波管の両端をPML（perfectly-matched layers）に割り当てる。導波管の壁はPEC（perfect electric conductor）にする。0=PEC、3=PML。

BC = [0 0 0 0 3 3]; %pml in pos. and neg. z-direction
FDTD = SetBoundaryCond(FDTD,BC);


CSXCADのメッシュを設定する。SmoothMeshLines関数は、第1引数で指定した区間を、第2引数の分解能で分割してくれる関数である。

CSX = InitCSX();
mesh.x = SmoothMeshLines([0 a],mesh_res(1));
mesh.y = SmoothMeshLines([0 b],mesh_res(2));
mesh.z = SmoothMeshLines([0 length],mesh_res(3));
CSX = DefineRectGrid(CSX,unit,mesh);

次は導波管のポートを設定する。ポート1は下端から11セルの位置、ポート2は上端から13セルの位置に置く。ポートは直方体の形であり、励振（excitation）は始点の位置、電圧・電流のプローブは終点の位置となる。終点はポートの基準面となる。
AddRectWaveGuidePort関数の第10引数は振幅であり、ポート1はActive（振幅=1で励振）、ポート2はPassive（振幅を指定しない or 振幅=0）としている。

start=[mesh.x(1)   mesh.y(1)   mesh.z(11)];
stop =[mesh.x(end) mesh.y(end) mesh.z(15)];
[CSX, port{1}] = AddRectWaveGuidePort(CSX,0,1,start,stop,'z',a*unit,b*unit,TE_mode,1);
 
start=[mesh.x(1)   mesh.y(1)   mesh.z(end-13)];
stop =[mesh.x(end) mesh.y(end) mesh.z(end-14)];
[CSX, port{2}] = AddRectWaveGuidePort(CSX,0,2,start,stop,'z',a*unit,b*unit,TE_mode);

時間領域の電界分布（E-field）をキャプチャする領域を定義する。AddDump関数の中で'FileType'属性を指定しているが、これは保存するファイル形式（0=vtkファイル、1=hdf5ファイル）を決めている。

CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','2,2,2');
start = [mesh.x(1)   mesh.y(1)   mesh.z(1)];
stop  = [mesh.x(end) mesh.y(end) mesh.z(end)];
CSX = AddBox(CSX,'Et',0,start,stop);


最後に、FDTDオブジェクトとCSXCADオブジェクトをxmlファイルに書き込み、シミュレーションを実行する。

Sim_Path = 'tmp_mod';
Sim_CSX = 'rect_wg.xml';
 
[status, message, messageid] = rmdir(Sim_Path,'s');
[status, message, messageid] = mkdir(Sim_Path);
 
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
RunOpenEMS(Sim_Path,Sim_CSX)

結果

電界分布

下のスクリプトでE-fieldをプロットすることができる。まああってそう。

dump_file = [Sim_Path '/Et.h5'];
PlotArgs.slice = {a/2*unit b/2*unit 0};
PlotArgs.pauseTime=0.01;
PlotArgs.component=0;
PlotArgs.Limit = 'auto';
PlotHDF5FieldData(dump_file,PlotArgs)

f:id:neurois:20200705064820g:plain
E-field

結果をHDF5ファイルで保存したことで、ParaViewを使わずにOctaveだけで確認することができる。このHDF5というフォーマット、初めて知ったのだが一般に使われているのか。調べてみるとHDFViewというソフトで見ることができて、pythonとか他の言語でも読み書きできるっぽい。

Sパラメータ

Sパラメータを計算するためのスクリプト。「uf」は周波数領域の電圧、「if」は周波数領域の電流、「inc」は入射（incident）、「ref」は反射（reflected）、「tot」はトータルという意味、、、らしいがかなり分かりにくい。

freq = linspace(f_start,f_stop,201);
port = calcPort(port, Sim_Path, freq);
 
s11 = port{1}.uf.ref./ port{1}.uf.inc;
s21 = port{2}.uf.ref./ port{1}.uf.inc;
ZL = port{1}.uf.tot./port{1}.if.tot;
ZL_a = port{1}.ZL; % analytic waveguide impedance

で、プロットするためのスクリプト。

figure
plot(freq*1e-6,20*log10(abs(s11)),'k-','Linewidth',2);
xlim([freq(1) freq(end)]*1e-6);
grid on;
hold on;
plot(freq*1e-6,20*log10(abs(s21)),'r--','Linewidth',2);
l = legend('S_{11}','S_{21}','Location','Best');
set(l,'FontSize',12);
ylabel('S-Parameter (dB)','FontSize',12);
xlabel('frequency (MHz) \rightarrow','FontSize',12);

f:id:neurois:20200705071502p:plain
S-parameters

導波管のインピーダンスを理論値と数値計算で比較。一致しててすごい（すごい）。

figure
plot(freq*1e-6,real(ZL),'Linewidth',2);
hold on;
grid on;
plot(freq*1e-6,imag(ZL),'r--','Linewidth',2);
plot(freq*1e-6,ZL_a,'g-.','Linewidth',2);
ylabel('ZL (\Omega)','FontSize',12);
xlabel('frequency (MHz) \rightarrow','FontSize',12);
xlim([freq(1) freq(end)]*1e-6);
l = legend('\Re(Z_L)','\Im(Z_L)','Z_L analytic','Location','Best');
set(l,'FontSize',12);

f:id:neurois:20200705071504p:plain
Waveguide impedance