Field of a grounded sphere - Scilab

[elodie_1]

Homework Statement

There is a grounded sphere of radius R in the origin of the coordinate system. In the distance L (L>R) from the sphere’s center there is a point charge Q. The electric field (both intensity and potential) should be computed in the area of radius r_g = 5L (in the plane containing the sphere’s center and the charge).

The solution should be delivered in Scilab.

Homework Equations

Poisson and Laplace equations solved numerically.

∇² $\varphi$ = 0 Laplace
∇² $\varphi$ = f Poisson

The Attempt at a Solution

Using the Taylor series expansion of derivatives, one can convert the differential
equation to the difference equation.

In the case of the region outside the sphere, where the Laplace equation is used, basically, the electric potential in a given point is computed as the arithmetic average from its four neighboring potentials, the north, west, east and south potential.

However, this is not true for the region belonging to the point charge, where I assume, the Poisson equation should be used. Can you help me finding a way to use the Poisson equation in the code? Can you take a look at the code to check if until know everything is correct?

Are my assumptions with Laplace and Poisson equations right?

My code so far:clear // clear all the variables present in the memory
xdel(winsid()) //close all the graphical windows
multiplier = 3; //used to adjust the density of the grid
density = 30 * multiplier; //grid's density
L = 8;
Rg = 2*L;
Rc = 0.5;
R = 5;
q = 10;
X0 = - Rg;
Xk = Rg;
Y0 = - Rg;
Yk = Rg;
x = linspace(X0, Xk, density);
y = linspace(Y0, Yk, density);
[X,Y] = ndgrid(x,y);
dx = x(2) - x(1);
dy = y(2) - y(1);

sphere = zeros(density, density);
r = sqrt(X.^2 + Y.^2);
charge = sqrt ((X-L).^2 + Y.^2);
sphere(find(r <= R)) = 0;
sphere(find(charge <= Rc)) = q;
//xset('colormap', jetcolormap(256))
//plot3d1(x,y,sphere)
dok=1e-3;
count=0;
jeszcze=%t;
while jeszcze do
count = count + 1;
nu=([sphere(:,$), sphere(:,1:$-1)]+[sphere(:,2:$),sphere(:,1)]+[sphere(2:$,:);sphere(1,:)]+[sphere($,:);sphere(1:$-1,:)])/4;
nu(find(r <= R)) = 0;
nu(find(charge <= Rc)) = q;
nu(:,1)= (4*nu(:,2)-nu(:,3))/3;
nu(:,$)= (4*nu(:,$-1)-nu(:,$-2))/3;
nu(1,:)= (4*nu(2,:)-nu(3,/3;
nu($,:)= (4*nu($-1,:)-nu($-2,/3;

delta = max(abs(sphere-nu));
sphere = nu;
if (delta < dok) then jeszcze = %f; end;
end;

xset('colormap', jetcolormap(64))
plot3d1(x,y,sphere)
xtitle(msprintf('%g loops were performed, max incerement = %g / %g',count,delta,dok),'x','y')

grx = -(sphere(3:$,2:$-1) - sphere(1:$-2,2:$-1))/2/dx;
gry = -(sphere(2:$-1,3:$) - sphere(2:$-1,1:$-2))/2/dy;
dlgr = sqrt (grx.^2+gry.^2);
maxdl=10;
wspdl=dlgr/maxdl;
wspdl(wspdl<1) = 1;
ngrx = grx./wspdl;
ngry = gry./wspdl;

scf();
contour2d(x,y,sphere,10)
xtitle("","x","y");

//scf();
co = 2;
champ(x(2:co:$-1), y(2:co:$-1), ngrx(1:co:$,1:co:$), ngry(1:co:$,1:co:$));
a =gca();
a.isoview = 'on';// check the correctness of the solution
ip=24*multiplier;
ik=26*multiplier;
jp=15*multiplier;
jk=22*multiplier;
plot2d(x([ip,ip,ik,ik, ip]), y([jp,jk,jk, jp, jp]), style = 5)
gr1 = -(sphere(ip:ik,jk+1) - sphere(ip:ik,jk-1))/2;
gr1(1) = gr1(1)/2; gr1($)=gr1($)/2;
gr2 = -(sphere(ik+1, jp:jk) - sphere(ik-1,jp:jk))/2;
gr2(1) = gr2(1)/2; gr2($)=gr2($)/2;
gr3 = (sphere(ip:ik,jp+1) - sphere(ip:ik,jp-1))/2;
gr3(1) = gr3(1)/2; gr3($) = gr3($)/2;
gr4 = (sphere(ip+1,jp:jk)-sphere(ip-1,jp:jk))/2;
gr4(1)=gr4(1)/2; gr4($)=gr4($)/2;
F1 = sum(gr1);
F2 = sum(gr2);
F3 = sum(gr3);
F4 = sum(gr4);
FLUX = F1+F2+F3+F4;
xtitle(msprintf('Number of divisions is %g x %g Net flux = %g \nF1=%g F2=%g F3=%g F4=%g',density,density,FLUX,F1,F2,F3,F4),'x','y')You can try it out on your own.

 

[krome]

That's really cool. I know nothing about coding so I can't help you. I'm just wondering if you've compared your result with the known exact result via the image charge method. It would be cool to see the quality of this lattice method.