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
[B]% now here this seems like  I am calculating wrongly [/B]
[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, :);

[B]% and the capacitance [/B]
% 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)');

 

[haruspex]

phi(i+1, j) * (1 + h / (2 * i))

I am not seeing the logic behind this expression. Please explain it.
E.g., i is a dimensionless number but h has units of metres, so 1 + h / (2 * i) is dimensionally inconsistent.

 

[kirito]

I am not seeing the logic behind this expression. Please explain it.
E.g., i is a dimensionless number but h has units of metres, so 1 + h / (2 * i) is dimensionally inconsistent.

you are correct I have preciously defined r=hi so if I have r_2=2h r_3=3h and so on then I said I will substitute it instead and forgot to delete the h from the upper term making it into 1+ 1/(2)

 

[kirito]

I did improve the code by a lot since I had a lot of confusion In application especially the relaxation method it is a bit hard for me visualise how it exactly works In a table I will add an updated version now

 

[kirito]

I am not seeing the logic behind this expression. Please explain it.
E.g., i is a dimensionless number but h has units of metres, so 1 + h / (2 * i) is dimensionally inconsistent.

with the dimension inconsistency dealt with ,really shows how important it is to check units it took me until yesterday morning to notice it
[CODE lang="matlab" title="an improved code with comments for a bit of clarity"]% Parameters
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
epsilon0 = 8.854187817e-12; % Vacuum permittivity

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

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

% Relaxation method
max_iterations = 1000000;
tolerance = 1e-6;
for iter = 1:max_iterations
old_phi = phi;
for i = 2:nr-1
for j = 2:nz-1
if j == ceil(d/(2*h)) && i <= ceil(R/h)
continue
end
phi(i, j) = (1/4) * (old_phi(i+1, j) * (1 + 1 / (2 * i)) + ...
old_phi(i-1, j) * (1 - 1 / (2 * i)) + ...
old_phi(i, j+1) + old_phi(i, j-1));
end
end
phi(1,:) = phi(2,:);
% Check for convergence
if max(max(abs(phi - old_phi))) < tolerance
fprintf('Converged after %d iterations\n', iter);
break;
end
end

% Compute the electric field using central difference approximation
Er = zeros(nr, nz);
Ez = zeros(nr, nz);
for i = 2:nr-1
for j = 2:nz-1
Er(i, j) = -(phi(i+1, j) - phi(i-1, j)) / (2 * h);
Ez(i, j) = -(phi(i, j+1) - phi(i, j-1)) / (2 * h);
end
end

% Compute the surface charge density at the edge of the plate
%r_edge = ceil(R / h);
%sigma = epsilon0 * (Er(r_edge+1, :) - Er(r_edge-1, / (2 * h);
% Compute the surface charge density at the edge of the plate
% Compute the surface charge density at the plate (at z = d/2)
z_plate = ceil(d / (2 * h)); % The index of z at the plate
sigma = epsilon0 * (Ez(:, z_plate+1) - Ez(:, z_plate-1)) / (2 * h);
% Compute the total charge on the plate by integrating sigma
% 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
fprintf('The total charge computed by integrating over the surface charged density:%e F\n', Q)
% Compute the capacitance
C = Q / V0;
fprintf('Computed capacitance: %e F\n', C);

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

% 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)');

% 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:nr-1) * h, sigma);
title('Surface Charge Density \sigma at the Edge of the Plate');
xlabel('r (m)');
ylabel('\sigma (C/m^2)');

% Numerical differentiation example
f = @cos;
x = 1;

d_f =@(x) -sin(x);
exact = d_f(x);
approx_1 = @(f, x, h) (f(x + h) - f(x)) ./ h;
approx_2 = @(f, x, h) (f(x + h) - f(x - h)) ./ (2 * h);
h = logspace(-1, -15, 100);

% Plot numerical differentiation errors
figure;
loglog(h, abs(exact - approx_1(f, x, h)), '*');
hold on;
loglog(h, abs(exact - approx_2(f, x, h)), '*');
hold off;
xlabel("log(h) from x=1");
ylabel("log(of the difference) from -sin(1)");
legend("normal approx", "symmetrical approximation");

[/CODE]

 

[haruspex]

This still looks wrong to me:

old_phi(i+1, j) * (1 + 1 / (2 * i)) + ... old_phi(i-1, j) * (1 - 1 / (2 * i)) + ... old_phi(i, j+1) + old_phi(i, j-1)

Why the 2*?

You seem to have fixed 0 potential at r=0 and r=rmax. Wouldn’t it be better to have it as linearly Vmin to Vmax at r=0, and maybe also at rmax?
Also, you set V0/2 at one plate, so shouldn’t you set -V0/2 at the other?

 

[kirito]

This still looks wrong to me:

Why the 2*?

You seem to have fixed 0 potential at r=0 and r=rmax. Wouldn’t it be better to have those as linearly Vmin to Vmax? Or maybe that at r=0 and uniformly Vavg at rmax?
Also, you set V0/2 at one plate, so shouldn’t you set -V0/2 at the other?

the 2* because i was asked to take r as ih which is the step taken in the r direction
Why did not i set the potential on the other plate to -vo/2 because i only looked at the upper plate taking it only into account

 

[kirito]

the 2* because i was asked to take r as ih which is the step taken in the r direction
Why did not i set the potential on the other plate to -vo/2 because i only looked at the upper plate taking it only into account

Though i may i have applied it incorrectly but when i asked a teacher for a hint to simplify the calculations since the run time is long , he told me a hint is to only look at the upper plate and try to find symmetry for the problem in the r theta coordinates, so i looked at the upper plate and only calculated the potential in a quarter of the space assuming that from the symmetry for rotation the potential i would get in each of the 4 quarters is the same

 

[kirito]

Though i may i have applied it incorrectly but when i asked a teacher for a hint to simplify the calculations since the run time is long , he told me a hint is to only look at the upper plate and try to find symmetry for the problem in the r theta coordinates, so i looked at the upper plate and only calculated the potential in a quarter of the space assuming that from the symmetry for rotation the potential i would get in each of the 4 quarters is the same

This still looks wrong to me:

Why the 2*?

You seem to have fixed 0 potential at r=0 and r=rmax. Wouldn’t it be better to have it as linearly Vmin to Vmax at r=0, and maybe also at rmax?
Also, you set V0/2 at one plate, so shouldn’t you set -V0/2 at the other?

I still have trouble understanding what potential to apply and for which condition, but i will reread this and try to understand how to apply what you are suggesting, i think it all comes down to the fact that though i know how to calculate relaxation I don’t understand how to apply Laplace conditions to most problems in question , i will be going to the office hours to see the code the lecturer wrote later on since the assignment time was due a while back , and i will read more about Laplace today

 

[haruspex]

the 2* because i was asked to take r as ih which is the step taken in the r direction

For a given circle (r value and z value) you have four adjacent circles: (r+h,z), (r-h,z), (r,z-h), (r,z+h).
The contributions of these to the potential at (r,z) at the next time step must be weighted according their "volumes": r+h, r-h, r, r. Normalising, 1+h/r, 1-h/r, 1, 1. But you have 1+h/(2r) etc.

Why did not i set the potential on the other plate to -vo/2 because i only looked at the upper plate taking it only into account

Yes, sorry, I forgot you were modelling only the upper half.
But what about the potentials at r=0 and r=rmax? It is not right for those to be fixed at 0.

 

[kirito]

For a given circle (r value and z value) you have four adjacent circles: (r+h,z), (r-h,z), (r,z-h), (r,z+h).
The contributions of these to the potential at (r,z) at the next time step must be weighted according their "volumes": r+h, r-h, r, r. Normalising, 1+h/r, 1-h/r, 1, 1. But you have 1+h/(2r) etc.

Yes, sorry, I forgot you were modelling only the upper half.
But what about the potentials at r=0 and r=rmax? It is not right for those to be fixed at 0.

The reason for the two because the formula i got when using symmetrical first derivative to approximate Laplace in spherical coordinates was this

 

[kirito]

View attachment 348616
The reason for the two because the formula i got when using symmetrical first derivative to approximate Laplace in spherical coordinates was this

For r=0 i thought its the condition so potential at z=0 is zero
For rmax i though i should take it as 0 since at infinite distance we consider potential zero

 

[haruspex]

For r=0 i thought its the condition so potential at z=0 is zero

Surely at r=0 the potential should increase approximately linearly from 0 to V0/2.

For rmax i though i should take it as 0 since at infinite distance we consider potential zero

It depends whether rmax is sufficiently large in relation to the disc radius. You only have it 25% larger. Maybe experiment with that later when you have it basically working, i.e. see whether increasing rmax makes a significant difference.

 

[haruspex]

Surely at r=0 the potential should increase approximately linearly from 0 to V0/2.

Better would be to set each potential at r=0 at t+1 from the potential at r=h, same z, at time t.

 

[kirito]

Better would be to set each potential at r=0 at t+1 from the potential at r=h, same z, at time t.

if you don't mind further explaining what this means , not sure if I get it , to set the potential at t+1 from r=h is bit unclear to me

 

[haruspex]

if you don't mind further explaining what this means , not sure if I get it , to set the potential at t+1 from r=h is bit unclear to me

At each time step you set the potential at each (r,z) coordinate from its neighbours' potentials at the previous time step, but not for r=0. I'm saying set it at (0,z) from the value at (h,z) at the previous time step.
Or you could do an average involving (0,z+1), (0,z-1) too, but I don't think that changes much.

 

[kirito]

At each time step you set the potential at each (r,z) coordinate from its neighbours' potentials at the previous time step, but not for r=0. I'm saying set it at (0,z) from the value at (h,z) at the previous time step.
Or you could do an average involving (0,z+1), (0,z-1) too, but I don't think that changes much.

your explanation is much appreciated, thank you