find the electric field between 2 finite discs

[kirito]

Homework Statement

there are two finite disc distance d from each other in the z axis potential difference between them V_zero find capacitance q and electric field between them

Relevant Equations

relaxation method , Q/C=V

I did make the problem simpler by looking at the the part from d/2 down the upper plate
here are my initial parameters I am making my size step be h since lowering it may make calculating harder
I am especially getting weird results for the field and capacitance



R = 0.1; % Radius of the plates in meters (10 cm)
d = 0.005; % Separation distance in meters (0.5 cm)
V0 = 1; % Potential difference in volts
h = d / 20; % Step size
rmax = R + 5 * d; % Maximum r
zmax = 10 * d; % Maximum z


% Create the grid
nr = ceil(rmax / h);
nz = ceil(zmax / h);
phi = zeros(nr, nz);


phi(1:R/h, d/(2*h)) = V0/2;% Set the potential on the positive plate
phi(end, :) = 0%; Boundary condition far from the plates
phi(:, 1) = 0 ;


% Relaxation method)
max_iterations = 10000;
tolerance = 1e-6;
for iter = 1:max_iterations


old_phi = phi;

for i = 2:nr-1
for j = 2:nz-1


if j == d/(2*h) && i < R/h
continue
end


phi(i, j) = (1/4) * (phi(i+1, j) * (1 + h / (2 * i)) + ...
phi(i-1, j) * (1 - h / (2 * i)) + ...
phi(i, j+1) + phi(i, j-1));
end
end
phi(:,1) = old_phi(:,2);
% Check for convergence
if max(max(abs(phi - old_phi))) < tolerance
fprintf('Converged after %d iterations\n', iter);
break;

end
end



% Plot the potential
figure;
imagesc((0:nr-1) * h, (0:nz-1) * h, flipud(phi'));
colorbar;
title('Electric Potential \phi');
xlabel('r (m)');
ylabel('z (m)');



















% Compute the electric field
% now here this seems like I am calculating wrongly
[Er, Ez] = gradient(-phi, h);


% Compute the surface charge density at the edge of the plate
r_edge = ceil(R / h) + 1;
sigma = -Er(r_edge, :);


% and the capacitance
% Note: We multiply by 2*pi*r to account for cylindrical coordinates
Q = 0;
for j = 1:nz-1
r = (j - 1) * h;
area = 2 * pi * r * h; % Area of the thin ring
Q = Q + sigma(j) * area;
end


% Compute the capacitance
% Compute the capacitance
C = Q / V0;
fprintf('Computed capacitance: %e F\n', C);


% Analytic capacitance for comparison
epsilon0 = 8.854187817e-12; % Vacuum permittivity
C_analytic = epsilon0 * pi * R^2 / d;
fprintf('Analytic capacitance: %e F\n', C_analytic);






% Plot the potential along z at r = 0
figure;
plot((0:nz-1) * h, phi(1, :);
title('Electric Potential \phi along z at r = 0');
xlabel('z (m)');
ylabel('\phi (V)');


% Plot the electric field
figure;
quiver((0:nr-1) * h, (0:nz-1) * h, Er', Ez');
title('Electric Field');
xlabel('r (m)');
ylabel('z (m)');


% Plot the surface charge density
figure;
plot((0:nz-1) * h, sigma);
title('Surface Charge Density \sigma at the Edge of the Plate');
xlabel('z (m)');
ylabel('\sigma (C/m^2)');