Explanation of parallel plates capacitor implementation with finite element methods and with in-homogenous domain.

6 views (last 30 days)
suketu vaidya
suketu vaidya on 7 Feb 2021
Commented: Shubham Khatri on 16 Feb 2021
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

Sign in to answer this question.

Answers (1)

Shubham Khatri
Shubham Khatri on 16 Feb 2021
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

Categories

Find more on Special Functions in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!