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

7 visualizaciones (últimos 30 días)
suketu vaidya
suketu vaidya el 7 de Feb. de 2021
Comentada: Shubham Khatri el 16 de Feb. de 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

Iniciar sesión para responder a esta pregunta.

Respuestas (1)

Shubham Khatri
Shubham Khatri el 16 de Feb. de 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

Categorías

Más información sobre Other Formats en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by