Explanation of parallel plates capacitor implementation with finite element methods and with in-homogenous domain.
조회 수: 7 (최근 30일)
이전 댓글 표시
function cap
close all;
hx = 0.002;
vx = 0:hx:0.1;
hy = 0.002;
h = hx;
vy = 0:hy:0.1;
nx = length(vx);
ny = length(vy);
n = nx*ny;
A = zeros(n);
b = zeros(n,1);
eps1 = 4
eps2 = 14
for ix = 1:nx
for jy = 1:ny
i = (jy-1) * nx + ix; % the global index of the ix,jy vertex
x = vx(ix); % the geometrical x coordinate
y = vy(jy); % the geometrical y coordinate
if ix == 1
% (1)
A(i,i) = 1;
b(i) = 0;
elseif ix == nx
% (2)
A(i,i) = 1;
b(i) = 10;
elseif jy == 1
if ix < (nx+1)/2
% (3)
A(i,i) = eps1*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps1*(1/hy^2);
elseif ix > (nx+1)/2
% (4)
A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps2*(1/hy^2);
else
% (5)
A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = (eps1+eps2)*(1/hy^2);
end
elseif jy == ny
% (6)
A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-nx) = eps2*(2/hy^2);
elseif jy < (ny+1)/2 && ix < (nx+1)/2
% (7)
A(i,i) = eps1*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i-nx) = eps1*(1/hy^2);
A(i,i+nx) = eps1*(1/hy^2);
elseif jy > (ny+1)/2 || ix > (nx+1)/2
% (8)
A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps2*(1/hy^2);
A(i,i-nx) = eps2*(1/hy^2);
elseif jy < (ny+1)/2
% (9)
A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = ((eps1+eps2)/2)*(1/hy^2);
A(i,i-nx) = ((eps1+eps2)/2)*(1/hy^2);
elseif ix < (nx+1)/2
% (10)
A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps2*(1/hy^2);
A(i,i-nx) = eps1*(1/hy^2);
else
% (11)
A(i,i) = (eps1+(eps2*3))*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps2*(1/hy^2);
A(i,i-nx) = ((eps1+eps2)/2)*(1/hy^2);
end
end
end
% solve the problem
u = A\b;
s = reshape(u,nx,ny);
[X,Y] = meshgrid(vx,vy);
surf(X,Y,s)
end
댓글 수: 1
답변 (1개)
Shubham Khatri
2021년 2월 16일
Hello,
As far as the code is concerned, Please remove the last 'end' from the code for it to work as in the code below.
close all;
hx = 0.002;
vx = 0:hx:0.1;
hy = 0.002;
h = hx;
vy = 0:hy:0.1;
nx = length(vx);
ny = length(vy);
n = nx*ny;
A = zeros(n);
b = zeros(n,1);
eps1 = 4
eps2 = 14
for ix = 1:nx
for jy = 1:ny
i = (jy-1) * nx + ix; % the global index of the ix,jy vertex
x = vx(ix); % the geometrical x coordinate
y = vy(jy); % the geometrical y coordinate
if ix == 1
% (1)
A(i,i) = 1;
b(i) = 0;
elseif ix == nx
% (2)
A(i,i) = 1;
b(i) = 10;
elseif jy == 1
if ix < (nx+1)/2
% (3)
A(i,i) = eps1*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps1*(1/hy^2);
elseif ix > (nx+1)/2
% (4)
A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps2*(1/hy^2);
else
% (5)
A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = (eps1+eps2)*(1/hy^2);
end
elseif jy == ny
% (6)
A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-nx) = eps2*(2/hy^2);
elseif jy < (ny+1)/2 && ix < (nx+1)/2
% (7)
A(i,i) = eps1*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i-nx) = eps1*(1/hy^2);
A(i,i+nx) = eps1*(1/hy^2);
elseif jy > (ny+1)/2 || ix > (nx+1)/2
% (8)
A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps2*(1/hy^2);
A(i,i-nx) = eps2*(1/hy^2);
elseif jy < (ny+1)/2
% (9)
A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = ((eps1+eps2)/2)*(1/hy^2);
A(i,i-nx) = ((eps1+eps2)/2)*(1/hy^2);
elseif ix < (nx+1)/2
% (10)
A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps2*(1/hy^2);
A(i,i-nx) = eps1*(1/hy^2);
else
% (11)
A(i,i) = (eps1+(eps2*3))*(-(2/hx^2)-(2/hy^2));
A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
A(i,i-1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));
A(i,i+nx) = eps2*(1/hy^2);
A(i,i-nx) = ((eps1+eps2)/2)*(1/hy^2);
end
end
end
% solve the problem
u = A\b;
s = reshape(u,nx,ny);
[X,Y] = meshgrid(vx,vy);
surf(X,Y,s)
Hope it helps
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Other Formats에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!