openEMSを使ってみる。その４：2D Cylindrical Wave

openEMS

Tutorial: 2D Cylindrical Wave

今回は「2D Cylindrical Wave」のチュートリアル Tutorial: 2D Cylindrical Wave - openEMS をやってみた。

Octaveスクリプト

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

close all; clear; clc

シミュレーションの半径を2.56 m(=2560 mm)とする。領域をr方向に5つに分割して、半径が境界を跨ぐと、θ方向のメッシュの数が2倍になるようにする。こうすることで、中心付近の電磁界のオーバーサンプリングを防ぐことができる。

physical_constants
mesh_res = 10;
radius = 2560;
split = ['80,160,320,640,1280'];
split_N = 5;
height = mesh_res*4;

周波数は1 GHzとする。波源は中心から1300 mm、45°の位置にずらして配置する。

f0 = 1e9;
excite_offset = 1300;
excite_angle = 45;

FDTDオブジェクトを初期化する。InitFDTD関数の中でMultiGridのオプションに先ほどの領域分割を与える。境界条件は+r方向をPMLに割り当てる。

FDTD = InitFDTD('NrTS',100000,'EndCriteria',1e-4,'CoordSystem',1,'MultiGrid',split);
FDTD = SetGaussExcite(FDTD,f0,f0/2);
BC = [0 3 0 0 0 0];             % pml in positive r-direction
FDTD = SetBoundaryCond(FDTD,BC);

内側の領域（r<80 mm）でのθ方向のメッシュラインの数を50として、外側の領域（r>1280 mm）のメッシュラインの数を計算する。+1してるのは植木算的な意味。

% 50 mesh lines for the inner most mesh
% increase the total number of mesh lines in alpha direction for all sub-grids
N_alpha = 50 * 2^split_N + 1;

CSXCADを初期化してメッシュを定義する。DefineRectGrid関数は円筒座標のグリッドにも使える。

CSX = InitCSX('CoordSystem',1);
mesh.r = SmoothMeshLines([0 radius],mesh_res);
mesh.a = linspace(-pi,pi,N_alpha);
mesh.z = SmoothMeshLines([-height/2 0 height/2],mesh_res);
CSX = DefineRectGrid(CSX,1e-3,mesh);

波源を置く。AddExcitation関数の第3引数で、波源の種類を1=「E-field hard excitation」に設定する。「hard excitation」では、波源での電磁界が強制的に（他の入射してくる波に関係なく）決まるのに対し、「soft excitation」は波源が作る電磁界と他の進行波との重ね合わせとなる。

start = [excite_offset excite_angle/180*pi-0.01 -20];
stop =  [excite_offset excite_angle/180*pi+0.01  20];
 
CSX = AddExcitation(CSX,'excite',1,[0 0 1]);
CSX = AddBox(CSX,'excite',0,start,stop);

フィールドをダンプする領域を定義する。今回は時間領域での結果をvtk形式で保存し、周波数領域（Frequency=f0）の結果をhdf5形式で保存する。AddDump関数のDumpTypeを0にすると「E-field time-domain dump」、10にすると「E-field frequency-domain dump」となる。

start = [mesh.r(1)   mesh.a(1)   0];
stop =  [mesh.r(end-8) mesh.a(end) 0];
 
% time domain vtk dump
CSX = AddDump(CSX,'Et_ra','DumpType',0,'FileType',0,'SubSampling','4,10,1');
CSX = AddBox(CSX,'Et_ra',0,start,stop);
 
% frequency domain hdf5 dump
CSX = AddDump(CSX,'Ef_ra','DumpType',10,'FileType',1,'SubSampling','2,2,2','Frequency',f0);
CSX = AddBox(CSX,'Ef_ra',0,start,stop);

FDTDとCSXをxmlに書き込んで、シミュレーションを実行する。

Sim_Path = 'tmp';
Sim_CSX = '2D_CC_Wave.xml';
 
[status, message, messageid] = rmdir(Sim_Path,'s'); % clear previous directory
[status, message, messageid] = mkdir(Sim_Path); % create empty simulation folder
 
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
RunOpenEMS(Sim_Path,Sim_CSX);


作成した3DモデルをCSXGeomPlot関数で確認してみる。円形にメッシュがあって、中心から離れて45°の位置にぽつんと青色のポートが見えているので良さそうだ。

CSXGeomPlot([Sim_Path '/' Sim_CSX]);

f:id:neurois:20200714003341p:plain
AppCSXCAD

周波数領域の電界分布をプロット

出力したhdf5形式のファイルを読み込み、電界分布をプロットする。

[field mesh_h5] = ReadHDF5Dump([Sim_Path '/Ef_ra.h5']);

r = mesh_h5.lines{1};
a = mesh_h5.lines{2};
a(end+1) = a(1);            %closeup mesh for visualization
[R A] = ndgrid(r,a);
X = R.*cos(A);
Y = R.*sin(A);

Ez = squeeze(field.FD.values{1}(:,:,1,3));
Ez(:,end+1) = Ez(:,1);      %closeup mesh for visualization

E_max = max(max(abs(Ez)));  %get maximum E_z amplitude

while 1
    for ph = linspace(0,360,41) %animate phase from 0..360 degree
        surf(X,Y,real(Ez*exp(1j*ph*pi/180)),'EdgeColor','none')
        caxis([-E_max E_max]/10)
        zlim([-E_max E_max])
        pause(0.3)
    end
end

f:id:neurois:20200714012845g:plain
E_z(f) animation