FDTD 法による電磁界シミュレーションを実践していきます。
これまでの記事はこちら↓
前回までで電界と磁界について空間と時間の差分化を実施しました。
今回は差分化したものをコード化していき、最終的には磁界について空間の伝搬を計算して可視化してきます。なお、3次元での計算は大変なため、今回は2次元での電磁界の伝搬について扱います。
二次元での FDTD 法
前回までの差分化の式は、3次元に対して立式していたため、まずは次元を1つ下げて差分方程式を立式します。
前回は立方体に対して式を立式しましたが、今回は立方体のうち1平面に着目して立式します。

ここで注目したいのが、電界はx,y成分のみを持ち、磁界はz成分のみを持つということです。
すなわち、$E_{z}=H_{x}=H_{y}=J_{0z}=0$であるということです。これを前回の記事の最後に立式した方程式に対して適用すると、
磁界は
\begin{align}
H_{x}^{n+\frac{1}{2}}\Biggl( i,j+\frac{1}{2},k+\frac{1}{2}\Biggr)=0\tag{1} \\
H_{y}^{n+\frac{1}{2}}\Biggl( i+\frac{1}{2},j,k+\frac{1}{2}\Biggr)=0\tag{2}\\
H_{z}^{n+\frac{1}{2}}\Biggl( i+\frac{1}{2},j+\frac{1}{2},k\Biggr)=H_{z}^{n-\frac{1}{2}}\Biggl(i+\frac{1}{2},j+\frac{1}{2},k\Biggr)-\frac{\Delta t}{\mu \Delta d}\Biggl(E_{y}^n\Bigl(i+1,j+\frac{1}{2},k\Bigr)-E_{y}^n\Bigl(i,j+\frac{1}{2},k\Bigr)-E_{x}^n\Bigl(i+\frac{1}{2},j+1,k\Bigr)+E_{x}^n\Bigl(i+\frac{1}{2},j,k\Bigr)\Biggr) \tag{3}\\
\end{align}
電界は
\begin{align}
E_{x}^{n+1}\Biggl( i+\frac{1}{2},j,k\Biggr)=AE_{x}^{n}\Biggl(i+\frac{1}{2},j,k\Biggr)-BJ_{0x}\Biggl( i+\frac{1}{2},j,k\Biggr)+\frac{C}{\Delta d}\Biggl(H_{z}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j+\frac{1}{2},k\Bigr)-H_{z}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j-\frac{1}{2},k\Bigr)\Biggr) \tag{4} \\
E_{y}^{n+1}\Biggl( i,j+\frac{1}{2},k\Biggr)=AE_{y}^{n}\Biggl(i,j+\frac{1}{2},k\Biggr)-BJ_{0y}\Biggl( i,j+\frac{1}{2},k\Biggr)+\frac{C }{\Delta d}\Biggl(H_{z}^{n+\frac{1}{2}}\Bigl(i+\frac{1}{2},j,k+\frac{1}{2}\Bigr)+H_{z}^{n+\frac{1}{2}}\Bigl(i-\frac{1}{2},j+\frac{1}{2},k\Bigr)\Biggr) \tag{5}\\
E_{z}^{n+1}\Biggl( i,j,k+\frac{1}{2}\Biggr)=0\tag{6}\\
\end{align}
なお、A, B, Cはそれぞれ以下の通りです。
\begin{align}
A&=\frac{1-\frac{\sigma \Delta t}{2\varepsilon}}{1+\frac{\sigma \Delta t}{2\varepsilon}} \\
B&=\frac{\frac{\Delta t}{\varepsilon}}{1+\frac{\sigma \Delta t}{2\varepsilon}} \\
C&=\frac{\frac{\Delta t}{\varepsilon}}{1+\frac{\sigma \Delta t}{2\varepsilon}}
\end{align}
これにて2次元の電界・磁界の方程式は準備完了です。
コーディングしてみる
今回はMatlabで電磁界の計算をしていきます。特別なモジュール等は使用しないため、Python等他の言語でもコーディング可能かと思います。
解析空間・時間の設定
まずは解析空間と解析時間の設定です。今回のプログラムでは、1m四方のエリアのうち、xy方向にそれぞれ1000のブロックを配置しました。
また
clear vars all
close all
% 解析空間サイズ
xSize = 1;
ySize = 1;
%各方向のセル数
bl=250;
%セルのサイズ
d=xSize/bl;
% 解析終了時間
tend = 3e-10;
%図示
fig = figure;
imag = imagesc(xPos(1,:),yPos(:,1)',zHNow);
axis xy
axis equal
axis tight
hold on
colormap jet
各種変数の設定
次に、計算に使用する変数を設定します。座標の設定方法は、電界・磁界の方程式を導出する際にセルの端を原点としていたため、同様の設定としました。
時間ステップについては、解析を収束させるための条件がありましたので、こちらを適用します。詳しくはCourantの安定条件で調べてみましょう。(他人任せですみません)
また、電界・磁界の現在と過去の変数が出てくるので、過去は”Pre”、現在は”Now”を変数の語尾につけて区分してます。
初期条件で必要となる変数も設定しました。
今回は周波数1GHzで、単発の正弦波を入力させて解析をしようと思います。
% 光速
c = 2.99792458e8;
% 時間ステップ
dx=xSize/(bl/2);
dy=ySize/(bl/2);
dt = 1/c/sqrt((1/dx)^2+(1/dy)^2);
t = 0:dt:tend;
% 空気の誘電率と透磁率
AirEps = 8.8541878128e-12;
AirMu = 1.25663706212e-6;
f=1e9;
tt= 1.0/f; %period
npulse = round(tt/dt/2); %number of pulse excitation steps
% 座標
[xPos, yPos] = meshgrid(0:d:xSize,ySize:-d:0);
% 電場
[xEPre,xENow,yEPre,yENow]=deal(zeros(size(xPos)));
% 電流密度
[xJ,yJ]=deal(zeros(size(xPos)));
% 磁場
[zxHPre,zxHNow,zyHPre,zyHNow,zHPre,zHNow]=deal(zeros(size(xPos)));
% 比誘電率,比透磁率
[Eps,Mu]=deal(ones(size(xPos)));
% 導電率
Sig= deal(zeros(size(xPos)));
% 初期電流
xJE=zeros(size(xPos));
電界・磁界の方程式
ここがコーディングの山場です。
電界と磁界の計算に入ります。高速化には向かないですが、わかりやすいのでforループで計算していきます。初期電流はx方向に正弦波を1/4周期印加させます。
電界と磁界については配列を使用して定義しています。
$E_{x}$については、座標$(i+\frac{1}{2},j,k)$を配列の$(idx,idy)$に格納し、
$E_{y}$については、座標$(i,j+\frac{1}{2},k)$を配列の$(idx,idy)$に格納し、
$H_{z}$については、座標$(i+\frac{1}{2},j+\frac{1}{2},k)$を配列の$(idx,idy)$に格納して計算しています。
なお、境界条件を設定する必要があるのですが、簡単のため全て完全導体としました。つまり、境界の接線成分の電界を0としています。
(今後、吸収境界条件なども試していきたいと思います。)
%磁界、電界の定数
AE=(2.*(AirEps.*Eps)-Sig.*dt)./(2.*(AirEps.*Eps)+Sig.*dt);
BE=(2*dt)./(2.*(AirEps.*Eps)+Sig.*dt);
CE=(2*dt)./(2.*(AirEps.*Eps)+Sig.*dt);
%磁界の定数
AH=dt./((AirMu*Mu));
for n = 1:300
%波源設定
yEmitIdx=round(ySize/d/2);
xEmitIdx=round(xSize/d/2);
xJ = sig(n,npulse);
yJ = 0;
xJE = zeros(size(xENow));
xJE(xEmitIdx, yEmitIdx) = xJ;
yJE = zeros(size(yENow));
yJE(xEmitIdx, yEmitIdx) = yJ;
%電界を計算
xENow = AE.*xEPre + BE.*(...
[zHNow(1,:); (zHNow(2:end,:)-zHNow(1:end-1,:))]./dx-xJE);
yENow = AE.*yEPre + BE.*(...
-[zHNow(:,1), (zHNow(:,2:end)-zHNow(:,1:end-1))]./dy-yJE);
% 外周をPECに設定
yENow(:,end) = 0;
yENow(:,1) = 0;
xENow(1,:) = 0;
xENow(end,:) = 0;
%磁界を計算
zHNow=zHPre-AH.*([yENow(:,2:end)-yENow(:,1:end-1), yENow(:,end)]./dx...
-[xENow(2:end,:)-xENow(1:end-1,:); xENow(end,:)]./dy);
%電界の更新
xEPre = xENow;
yEPre = yENow;
%磁界の更新
zHPre = zHNow;
% 描画
imag.CData = zHNow;
drawnow
end
%電流源
function y = sig(x,npulse)
if x <= npulse/2 %正弦波の1/4波長のみ出力
y = -sin(pi*x/npulse)^2; %平滑化するため正弦波を2乗
else
y=0; %continuous sine wave
%y = -sin(pi*x/npulse);
end
end
解析結果
上記プログラムを実行した結果になります。
波源は(0.5,0.5)に設定しました。波源を中心に同心円状に磁界が伝搬していることがわかります。また、境界(グラフの端部)において、磁界が完全反射しているため、境界条件の設定も問題なさそうです。
まとめ
今回はFDTD法を用いた電磁界解析を自力で実装してみました。電磁波工学を学んでも研究や製品開発をする方はほとんどシミュレーションを使用するため方程式を解くことなどほとんどないと思います。かくいう私も今回の取り組みまでシミュレーションに頼っていました・・・
電磁波の理解を深める良い機会になりましたので皆さんもよかったら一度計算して見ると良いと思います。
マクスウェルの方程式を自分で解いてみて 電磁波 の理解を深める
FDTD 法の基礎を習得する