Explanation of parallel plates capacitor implementation with finite element methods and with in-homogenous domain.
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentar
Antworten (1)
Shubham Khatri
am 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
0 Kommentare
Siehe auch
Kategorien
Mehr zu Other Formats finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!