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

5 vues (au cours des 30 derniers jours)
suketu vaidya
suketu vaidya le 7 Fév 2021
Commenté : Shubham Khatri le 16 Fév 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

Connectez-vous pour répondre à cette question.

Réponses (1)

Shubham Khatri
Shubham Khatri le 16 Fév 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

Catégories

En savoir plus sur Animation dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by