數值電磁 - 電容模擬器自己寫 (part2/3)
此篇延續[數值電磁 - 電容模擬器自己寫 (part1/3)]
計算電容的數值模擬器，此篇延伸單一介質(homogeneous)至多種介質(non-homogeneous)。
如下圖，中間為方形金屬
上圖為單一介質
下圖為多種介質
主要參考文獻，University of Utah - 教授JR Nagel
https://my.ece.utah.edu/~ece6340/LECTURES/Feb1/Nagel%202012%20-%20Solving%20the%20Generalized%20Poisson%20Equation%20using%20FDM.pdf
首先介紹相關方程
因為介電材料ε不再是單一常數，是會隨x,y改變的，所以寫成ε(x,y)。
Maxwell Gauss law且source free
🜄·(ε(x,y)E) = 0...(1)
接著定義電位φ和E的關係
E = -🜄φ...(2)
將式(2)帶入式(1)
🜄·(ε(x,y)🜄φ(x,y)) = 0
拆解🜄成偏微分，且為了方便先只討論二維x,y
🜄ε(x,y)·🜄φ(x,y) + ε(x,y)🜄🜄φ(x,y) = 0
[(მε(x,y)/მx)*(მφ(x,y)/მx)] + [(მε(x,y)/მy)*(მφ(x,y)/მy)] + [ε(x,y)*(მ²φ/მx²+მ²φ/მy²)] = 0
採用「中央差分方」把討人厭的微分幹掉
且離散x,y為xx,yy，並把xy的mesh都設定成一樣大△x=△y=h
[(ε(xx+1,yy)-ε(xx,yy))/h * (φ(xx+1,yy)-φ(xx,yy))/h]
+ [(ε(xx,yy+1)-ε(xx,yy))/h * (φ(xx,yy+1)-φ(xx,yy))/h]
+ [ε(xx,yy)*(φ(xx+1,yy)+φ(xx,yy+1)+φ(xx-1,yy)+φ(xx,yy-1)-4φ(xx,yy))/h²]
= 0
接著把φ(xx,yy)整理到左邊
a0*φ(xx,yy) = a1*φ(xx,+1yy) + a2*φ(xx,yy+1) + a3*φ(xx-1,yy) + a4*φ(xx,yy-1)...(4)
其中
a0 = ε(xx+1,yy) + ε(xx,yy+1) + 2*ε(xx,yy)
a1 = ε(xx+1,yy)
a2 = ε(xx,yy+1)
a3 = ε(xx,yy)
a4 = ε(xx,yy)
觀察上式(4)可發現
non-homogeneous和homogeneous的電位方程架構是雷同的，只是non-homogeneous多了ε的係數。
「x,y點的電位 = 左、右、上、下四個點電位的組成」
若ε永遠為1，則式(4)也能變為Homogeirous的形式(向下相容)。
上述係數a0, a1, a2, a3, a4是由筆者自行推導，所以和文獻中不同。
文獻是導入積分型的Maxwell Gauss law，推得
a0 = ε(xx,yy) + ε(xx-1,yy) + ε(xx,yy-1) +ε(xx-1,yy-1)
a1 = 0.5 * (ε(xx,yy)+ε(xx,yy-1))
a2 = 0.5 * (ε(xx-1,yy)+ε(xx,yy))
a3 = 0.5 * (ε(xx-1,yy-1)+ε(xx-1,yy))
a4 = 0.5 * (ε(xx,yy-1)+ε(xx-1,yy-1))
不論是筆者自行推導的係數，或是文獻的係數，筆者都已測試過是可行的。
筆者也將兩組係數都寫在Matlab程式碼(文章最底下)，方便大家驗證。
求解問題如下圖，同軸方形波導管。
施加1V在中央金屬M1，0V在外圈金屬M2，我們想求內外金屬之間的電容。
兩金屬之間塞入3種介質，側面圖如下，包含DK1, 4, 10。
下圖為計算結果，mesh h為5e-5 m
我們用「收斂條件」來決定電位φ疊代的次數
透過相減疊代前和疊代前電位，若差異小於1e-4，我們可判定電位已經收斂，則停止疊代。
自製Matlab模擬器(2D)，金屬內對外單位電容為261pF。
Cadence Clarity Capacitance(3D)，金屬內對外單位電容為257 pF。
Ansys Q3D(3D)，金屬內對外單位電容為257pF。
MATLAB模擬結果
電位分布圖
電位等位圖
此次模擬，花費18982次疊代才達到收斂條件∆ε<1e-4，收斂過慢。
下次會介紹在模擬器中添加類似無限大空間的邊界條件。
「MATLAB程式碼 - 方形金屬同軸波導管單位電容器，內圈1cm外圈2cm，多種介電材料。 」
clc;clear;close all;
t0 = clock;
%金屬2 cm切400份
meshnum = 400;
dx = 0.02/meshnum; % dx=2cm/400
% Inner Conductor Region
metalsize = round(0.01/dx); % inner metal 1 cm
start = 1+floor((meshnum-metalsize)/2);
endmetal = meshnum - start;
metalPositiveVolt = ones(meshnum);
metalPositiveVolt(start:endmetal,start:endmetal)=0;
% phi Initial guess
phi = rand(meshnum,meshnum); %隨機初始化空間電位
% initial E,D,phi,Epsilon
E = zeros(meshnum);
Ex = zeros(meshnum);
Ey = zeros(meshnum);
D = zeros(meshnum);
Dx = zeros(meshnum);
Dy = zeros(meshnum);
Epsilon = ones(meshnum);
Epsilon(1:round(meshnum/2),1:round(meshnum/2)) = 4; %空間部分dk為4
Epsilon(round(meshnum/2)+1:end,1:round(meshnum/2)) = 10; %空間部分dk為10
Epsilon(:,1:round(meshnum/8)) = 1; %空間部分dk為1
% plot strcture
structure = Epsilon.*metalPositiveVolt; % plot structure
figure;pcolor(structure');colorbar;shading interp;
title('Non-homogeneous');
subtitle({['0 is metal'],['>0 is dieletric material']});
for nn=1:1e5 %疊代次數
phi_before = phi;
for xx=2:1:meshnum-1
for yy=2:1:meshnum-1
% 自己推
a0 = (Epsilon(xx+1,yy) + Epsilon(xx,yy+1) + 2*Epsilon(xx,yy) );
a1 = Epsilon(xx+1,yy) ; % x+1
a2 = Epsilon(xx,yy+1) ; % y+1
a3 = Epsilon(xx,yy) ; % x-1
a4 = Epsilon(xx,yy) ; % y-1
% 講義
% a0 = Epsilon(xx,yy) + Epsilon(xx-1,yy) + Epsilon(xx,yy-1) +Epsilon(xx-1,yy-1);
% a1 = 0.5 * (Epsilon(xx,yy)+Epsilon(xx,yy-1)); % x+1
% a2 = 0.5 * (Epsilon(xx-1,yy)+Epsilon(xx,yy)); % y+1
% a3 = 0.5 * (Epsilon(xx-1,yy-1)+Epsilon(xx-1,yy)); % x-1
% a4 = 0.5 * (Epsilon(xx,yy-1)+Epsilon(xx-1,yy-1)); %y-1
% Jacobi iteration of phi for laplace euqation
phi(xx,yy) ...,
= 1/a0* ...,
( phi_before(xx-1,yy)*a3 ...,% x-1
+ phi_before(xx,yy-1)*a4 ..., % y-1
+ phi_before(xx+1,yy)*a1 ..., % x+1
+ phi_before(xx,yy+1)*a2 ..., % y+1
);
end
end
% set B.C.
% Inner Conductor phi=1
phi(start:endmetal,start:endmetal)=1;
% Outter Conductor phi=0
phi(1,:)=0;
phi(:,1)=0;
phi(end,:)=0;
phi(:,end)=0;
deltaPhi = max(max(abs(phi-phi_before)));
if deltaPhi<1e-4 %收斂條件, phi差異小於1e-4
break;
end
end
t1=clock;
runtime = etime(t1,t0);
% calculate E
Ex(1:end-1,:) = (phi(2:end,:)-phi(1:end-1,:)) /dx;
Ey(:,1:end-1) = (phi(:,2:end)-phi(:,1:end-1)) /dx;
E = sqrt(Ex.*Ex + Ey.*Ey);
% calculate D
Dx = Ex.*Epsilon*8.8541878128e-12;
Dy = Ey.*Epsilon*8.8541878128e-12;
D = sqrt(Dx.*Dx + Dy.*Dy);
% calculate surface charge density
charge = (sum(abs(Dx(start-1,start-1:endmetal+1))) ...,
+ sum(abs(Dx(endmetal,start-1:endmetal+1))) ...;
+ sum(abs(Dy(start-1:endmetal+1,start-1))) ...,
+ sum(abs(Dy(start-1:endmetal+1,endmetal)))...,
) *dx*1e12; %in pF
% plot potential
figure; contour(phi'); hold on;
quiverArrowSpacing = round(meshnum/20);
quiver(1:quiverArrowSpacing:meshnum,1:quiverArrowSpacing:meshnum, ...;
-Ex(1:quiverArrowSpacing:meshnum,1:quiverArrowSpacing:meshnum)', ...;
-Ey(1:quiverArrowSpacing:meshnum,1:quiverArrowSpacing:meshnum)'); hold off;
title(['Capacitance = ',num2str(charge), ' pF']);
subtitle({['Iteration = ',num2str(nn), ', deltaPhi(%) = ',num2str(deltaPhi)], ...;
['dx(m) = ',num2str(dx)], ...;
['runtime = ',num2str(runtime),'(sec)']} ...;
);
colormap jet;colorbar;
% plot potential contour
figure; pcolor(phi'); shading interp;hold on;
quiverArrowSpacing = round(meshnum/20);
quiver(1:quiverArrowSpacing:meshnum,1:quiverArrowSpacing:meshnum, ...;
-Ex(1:quiverArrowSpacing:meshnum,1:quiverArrowSpacing:meshnum)', ...;
-Ey(1:quiverArrowSpacing:meshnum,1:quiverArrowSpacing:meshnum)'); hold off;
title(['Capacitance = ',num2str(charge), ' pF']);
subtitle({['Iteration = ',num2str(nn), ', deltaPhi(%) = ',num2str(deltaPhi)], ...;
['dx(m) = ',num2str(dx)], ...;
['runtime = ',num2str(runtime),'(sec)']} ...;
);
colormap jet;colorbar;
% plot E
figure; pcolor(E'); shading interp;colorbar;
title(['E field mag']);
% plot D
figure; pcolor(D'); shading interp; colorbar;
title(['D field mag']);
留言
張貼留言