openEMSを使ってみる。その３：Circular Waveguide

openEMS

Tutorial: Circular Waveguide

今回は「Circular Waveguide」（円形導波管）のチュートリアル Tutorial: Circular Waveguide - openEMS をやってみた。前回の「Rectangular Waveguide」と基本的に流れは同じなので、前回あまり説明しなかったSパラメータの計算方法についてもう少し詳しく書いていこうと思う。

Octaveスクリプト

いつものように実行環境を空にしておく。

close all; clear; clc

単位をmmとする。

physical_constants;
unit = 1e-3; %drawing unit in mm

導波管の長さを2 m(=2000 mm)、半径を35 cm(=350 mm)とする。

% waveguide dimensions
length = 2000;
rad = 350;     %waveguide radius in mm

着目する周波数帯域を300 MHzから500 MHzとする。

% frequency range of interest
f_start =  300e6;
f_stop  =  500e6;

メッシュの細かさを決める。今回は円筒座標系(ρ,φ,z)で与えている。

mesh_res = [10 2*pi/49.999 10]; %targeted mesh resolution

FDTDオブジェクトを初期化。座標系を円筒座標系(=1)としている。

FDTD = InitFDTD('EndCriteria',1e-4,'CoordSystem',1);
FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start));

境界条件を設定する。±Z方向をPML (perfectly matched layer) としている。

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

CSXCADを作成してメッシュを定義する。

CSX = InitCSX('CoordSystem',1); % init a cylindrical mesh
mesh.r = SmoothMeshLines([0 rad],mesh_res(1)); %mesh in radial direction
mesh.a = SmoothMeshLines([0 2*pi],mesh_res(2)); % mesh in azimuthal dir.
mesh.z = SmoothMeshLines([0 length],mesh_res(3));
CSX = DefineRectGrid(CSX,unit,mesh);

円形導波管のポートを追加する。モードはTE11。

start=[mesh.r(1)   mesh.a(1)   mesh.z(8)];
stop =[mesh.r(end) mesh.a(end) mesh.z(15)];
[CSX, port{1}] = AddCircWaveGuidePort(CSX,0,1,start,stop,rad*unit,'TE11',0,1);
 
start=[mesh.r(1)   mesh.a(1)   mesh.z(end-13)];
stop =[mesh.r(end) mesh.a(end) mesh.z(end-14)];
[CSX, port{2}] = AddCircWaveGuidePort(CSX,0,2,start,stop,rad*unit,'TE11');

電界分布を取得する領域を追加する。

CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','4,4,4');
start = [mesh.r(1)   mesh.a(1)   mesh.z(1)];
stop  = [mesh.r(end) mesh.a(end) mesh.z(end)];
CSX = AddBox(CSX,'Et',0,start,stop);

FDTDの設定とCSXCADデータをxmlに書き込んで、シミュレーションを実行する。

Sim_Path = 'tmp';
Sim_CSX = 'circ_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)

作成した3DモデルはAppCSXCADで確認できる。CSXGeomPlot関数でxmlファイルを指定するだけなので簡単である。

CSXGeomPlot([Sim_Path '/' Sim_CSX]);

f:id:neurois:20200712202824p:plain
AppCSXCAD

Sパラメータとインピーダンスの計算

ここの部分は前回説明を省いたので少し詳しく書いてみる。Sパラメータを計算するには以下の後処理が必要となる。着目している周波数のベクトルをcalcPort関数に渡している。この関数は、それぞれの周波数におけるポートの電圧と電流、伝搬定数、特性インピーダンスを計算してくれるようだ。

freq = linspace(f_start,f_stop,201);
port = calcPort(port,Sim_Path,freq);

Sパラメータは次のスクリプトで計算できる。ここで、「uf」は周波数に対する電圧、「if」は周波数に対する電流、「inc」は入射（incident）、「ref」は反射（reflected）、「tot」はtotalの意味である。

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;

S11 はポート1に入射する電力の波 a1 と 反射する電力の波 b1 の比で与えられるが、無損失の場合 a1 と b1 は電圧に比例するので、a1 はport{1}.uf.inc（ポート1に入射する電圧）、b1 はport{1}.uf.ref（ポート1で反射する電圧）に対応する。したがって上のような記述となる。


計算したSパラメータをプロットする。

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:20200712201739p: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,port{1}.ZL,'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:20200712201802p:plain
Waveguide impedance