風行電磁波模擬 (FonSinEM) - Free ElectroMagnetic Simulator

● 金屬-介電質-交界面 (Metal-Dielectric-Interface)

問題: 考慮一平面波入射一金屬層(材料為金)，觀察其反射頻譜 (100~1000 nm)並與理論解析解做比對

===================================================================

根據這個問題，上、下的位置是要加上吸收層(PML)防止波傳播至邊界時反彈干擾系統。那麼X方向呢？

這裡X的方向符合我們需要而能設定的有吸收層(PML)或週期性(Period)。若是加吸收層，則即是模擬一個我們所設定大小的介電質-金屬交界面結構。那為什麼可以用週期性？這是因為，我們要處理的交界面問題，在X方向都是相同的，也可以看成無窮遠延伸，因此週期的邊界條件也是適用的。

但兩者的觀測點位置是不一樣的

吸收層：觀測位置不能在吸收層裡(也就是若吸收層為16格，那觀測點就從至少17開始，請勿與吸收層位置重疊！)

週期性：觀測位置為整個系統X長度

如下圖

那麼，總尺寸的大小要設定多少？

Z方向的我們只要給予適當的大小即可(大於金屬膚深skin depth)。比如50,100網格都是可以的。只要擺的下結構與入射波源、量測的觀測點即可(建議波源與交界面之間還是有個10格的距離會較好，量測點再多個5格，因為金屬材料可能會有消逝波(evanecent wave)影響，怕觀測點受到干擾造成計算誤差)

X方向呢？

--若是X方向我們選擇用吸收層的話，那麼就必須要大於『兩端吸收層厚度』的總尺寸再加2格才可。也就是說若我們吸收層厚度為16格，那總尺寸的X方向，就至少要32+2=34格。

--若是週期性，因為週期邊界沒有厚度，而總尺寸任何方向最少都要2格網格，因此只要大於2格即可。

可以自行試看看，無論2格、10格,100格,1000格都會得到相同的結果，但網格數越大計算時間也將越長及記憶體使用將越大。

===================================================================

設定參數

A. 主要設定 (Main-Setup)

1. 定義總x,y,z長度相等於模型檔案，然後調整適當的各軸解析度

2. 按【建立網格 Cad to Grids (Create)】按鈕建立模擬網格

首先先選擇想要模擬的為維度，垂直入射平面波計算反射率的話、這裡選擇2D TM (ExEzHy)為例子。

設定總尺寸x=10,z=100; delta我們設5nm，因此系統總大小是x=50nm,z=500nm;

B. 邊界條件 (B.C.)

1. 定義邊界條件

2. 按【創建 (Create)】按鈕 => 建立總網格 (包含邊界條件&增加網格)

邊界條件我們X方向使用週期性，Z方向使用吸收層

C. 波源與觀測 (Wave & Ob.)

使用線性極化平面波由下方入射(theta=0，電場極化方向為Ey (psi=0)

反射頻譜觀測為-Z 方向

穿透頻譜觀測為+Z 方向

教學範例 : 觀測範圍 (DFT)

D. 結構設計 (Geo. Design)

1.勾選要使用的材料

2.輸入內建結構函數

3.執行區塊內的指令函數

4.檢查結構是否正確

5.確認後輸出結構

內建結構的寫法則如下

gdx是(grid) dx的大小, gdy = dy. gdz = dz;

ib= X grids, jb= Y grids, kb= Z grids;

icenter= ib/2; jcenter= jb/2; kcenter=kb/2;

內建結構函數矩形 Brick 參考 : 內建結構

Initializedgeometries

%=================矩形 Brick===========================

xstart=1*gdx;

xend=10*gdx;

ystart=1*dy;

yend=1*dy;

zstart=50*gdz;

zend=100*gdz;

nindex=1^2;

sigma=0;

choice='E_Model1'; %E_Iso,PEC,M_Iso,PMC,E_Model1,M_Model1,EM_Model1

gridtype=-1; %

Iso_Brick(choice,gridtype,nindex,sigma,xstart,xend,ystart,yend,zstart,zend)

%=================矩形 Brick===========================

Initializedgeometries : 初始化結構

gdx是grid dx的大小，這裡前面設定是5nm，而dy,dz則與dx相同

X從1到10格，Z從zstart=50*gdz至zend=100*gdz;

材料我們選擇'E_Model1'電性色散材料

確定無誤後，按下【輸出結構】即可輸出記憶體內描述結構的矩陣，好讓執行檔讀取

這時會產生Data_Materials.csv（因為只有色散材料，因此只需要此檔案即可）

E. 電,磁性材料模型(E,M Models)

設定色電性散材料的參數，我們選擇金，而風行有內建金材料的參數。

(若有用到色散材料則電性色散材料記得一定給參數，不然計算時會讀取色散參數失敗而錯誤。

F.頻譜 & 場 (Spectrum & Pattern)

1.設定頻譜分析 & 觀測波長範圍

2.檢察波源

G.模擬計算 (Calculation)

1.設定模擬時間長

2.存檔

3.計算

H.分析結果 (Result Analysis)

這個是FDTD所計算出來的結果，您可以多試試幾個看看，也可做三維、或者是用SplitField Method的計算試試。

理論解析解公式

clear ; clc;

%讀取電性色散參數檔案

a=textread('Data.Mats','%s','delimiter', 's');

Debyemodel1=str2num(char(a((1))));

Debyemodel2=str2num(char(a((2))));

Debyemodel3=str2num(char(a((3))));

Debyemodel4=str2num(char(a((4))));

Debyemodel5=str2num(char(a((5))));

Drudemodel1=str2num(char(a((6))));

Drudemodel2=str2num(char(a((7))));

Drudemodel3=str2num(char(a((8))));

Drudemodel4=str2num(char(a((9))));

Drudemodel5=str2num(char(a((10))));

Lorentzmodel1=str2num(char(a((11))));

Lorentzmodel2=str2num(char(a((12))));

Lorentzmodel3=str2num(char(a((13))));

Lorentzmodel4=str2num(char(a((14))));

Lorentzmodel5=str2num(char(a((15))));

Pademodel1=str2num(char(a((16))));

Pademodel2=str2num(char(a((17))));

Pademodel3=str2num(char(a((18))));

Pademodel4=str2num(char(a((19))));

Pademodel5=str2num(char(a((20))));

hDebyemodel1=str2num(char(a((21))));

hDebyemodel2=str2num(char(a((22))));

hDebyemodel3=str2num(char(a((23))));

hDebyemodel4=str2num(char(a((24))));

hDebyemodel5=str2num(char(a((25))));

hDrudemodel1=str2num(char(a((26))));

hDrudemodel2=str2num(char(a((27))));

hDrudemodel3=str2num(char(a((28))));

hDrudemodel4=str2num(char(a((29))));

hDrudemodel5=str2num(char(a((30))));

hLorentzmodel1=str2num(char(a((31))));

hLorentzmodel2=str2num(char(a((32))));

hLorentzmodel3=str2num(char(a((33))));

hLorentzmodel4=str2num(char(a((34))));

hLorentzmodel5=str2num(char(a((35))));

Pademodel1=str2num(char(a((36))));

Pademodel2=str2num(char(a((37))));

Pademodel3=str2num(char(a((38))));

Pademodel4=str2num(char(a((39))));

Pademodel5=str2num(char(a((40))));

epsinf1=str2num(char(a((71))));

sigmametal1=str2num(char(a((72))));

epsinf2=str2num(char(a((73))));

sigmametal2=str2num(char(a((74))));

epsinf3=str2num(char(a((75))));

sigmametal3=str2num(char(a((76))));

epsinf4=str2num(char(a((77))));

sigmametal4=str2num(char(a((78))));

epsinf5=str2num(char(a((79))));

sigmametal5=str2num(char(a((80))));

hepsinf1=str2num(char(a((81))));

hsigmametal1=str2num(char(a((82))));

hepsinf2=str2num(char(a((83))));

hsigmametal2=str2num(char(a((84))));

hepsinf3=str2num(char(a((85))));

hsigmametal3=str2num(char(a((86))));

hepsinf4=str2num(char(a((87))));

hsigmametal4=str2num(char(a((88))));

hepsinf5=str2num(char(a((89))));

hsigmametal5=str2num(char(a((90))));

gg=111;

%將各項係數分配入變數

for i=1:Debyemodel1

Detp1(i)=str2num(char(a((gg))));

gamabye1(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Debyemodel2

Detp2(i)=str2num(char(a((gg))));

gamabye2(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Debyemodel3

Detp3(i)=str2num(char(a((gg))));

gamabye3(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Debyemodel4

Detp4(i)=str2num(char(a((gg))));

gamabye4(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Debyemodel5

Detp5(i)=str2num(char(a((gg))));

gamabye5(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel1

Wd1(i)=str2num(char(a((gg))));

gamad1(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel2

Wd2(i)=str2num(char(a((gg))));

gamad2(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel3

Wd3(i)=str2num(char(a((gg))));

gamad3(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel4

Wd4(i)=str2num(char(a((gg))));

gamad4(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Drudemodel5

Wd5(i)=str2num(char(a((gg))));

gamad5(i)=str2num(char(a((gg+1))));

gg=gg+2;

end

for i=1:Lorentzmodel1

eps_p1(i)=str2num(char(a((gg))));

Wp1(i)=str2num(char(a((gg+1))));

gamap1(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

for i=1:Lorentzmodel2

eps_p2(i)=str2num(char(a((gg))));

Wp2(i)=str2num(char(a((gg+1))));

gamap2(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

for i=1:Lorentzmodel3

eps_p3(i)=str2num(char(a((gg))));

Wp3(i)=str2num(char(a((gg+1))));

gamap3(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

for i=1:Lorentzmodel4

eps_p4(i)=str2num(char(a((gg))));

Wp4(i)=str2num(char(a((gg+1))));

gamap4(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

for i=1:Lorentzmodel5

eps_p5(i)=str2num(char(a((gg))));

Wp5(i)=str2num(char(a((gg+1))));

gamap5(i)=str2num(char(a((gg+2))));

gg=gg+3;

end

cc=2.99792458e8;

muz=4.0*pi*1.0D-7;

epsz=1.D0/(muz*cc*cc);

lambda=linspace(100,1000,500)*1e-9;

twopic = 1.883651567308853e+09; % twopic=2*pi*c where c is speed of light

omega = twopic.*(lambda.^(-1)); % angular frequency of light (rad/s)

%=================================================

epsilon=epsinf1-sigmametal1./(j*omega*epsz);

%debye

for i=1:Debyemodel1

epsilon=epsilon-Detp1(i)./(1-j*omega*gamabye1(i));

end

%drude

for i=1:Drudemodel1

epsilon=epsilon-Wd1(i)^2./(omega.^2+j*omega*gamad1(i));

end

%lorentz

for i=1:Lorentzmodel1

epsilon=epsilon-eps_p1(i)^2./(omega.^2+j.*omega*gamap1(i)-Wp1(i)^2);

end

% epsilon_r 對 lambda

figure(1);hold on;

plot(lambda,real(epsilon),'r','linewidth',3);

plot(lambda,imag(epsilon),'b','linewidth',3);

hold off;

title('Au');

xlabel('\lambda (nm)');

ylabel('\epsilon r');

% 反射率頻譜

n=sqrt( (sqrt(real(epsilon).^2+imag(epsilon).^2)+real(epsilon))/2 );

k=sqrt( (sqrt(real(epsilon).^2+imag(epsilon).^2)-real(epsilon))/2 );

r=(n0-n1)./(n0+n1);

figure(1);

plot(lambda,abs(r).^2,'^b','linewidth',1);

title('R');

xlabel('\lambda (nm)');

ylabel('R');

Page updated

Google Sites

Report abuse