Uniformly Charged Sphere: Force Exerted by Southern Hemisphere on Northern Hemisphere

Electric Field Inside Sphere

Specify the radius of the sphere and charge density.

Rad = 1;
rho = 3;

Find the total charge of the sphere.

Q = 4/3*pi*Rad^3*rho;

Create an femodel object for electrostatic analysis and include the geometry of a sphere.

emagmodel = femodel(AnalysisType="electrostatic", ...
                    Geometry=multisphere(Rad));

Plot the geometry with face labels.

pdegplot(emagmodel,FaceLabels="on",FaceAlpha=0.3)

Specify the vacuum permittivity value in the SI system of units.

eps0 = 8.8541878128E-12;
emagmodel.VacuumPermittivity = eps0;

Specify the relative permittivity of the material.

emagmodel.MaterialProperties = ...
    materialProperties(RelativePermittivity=1);

Specify the charge density.

emagmodel.CellLoad = cellLoad(ChargeDensity=rho);

Specify the electrostatic potential on the boundary.

emagmodel.FaceBC = faceBC(Voltage=10);

Generate the mesh. This assignment updates the mesh stored in the Geometry property of the model.

emagmodel = generateMesh(emagmodel,Hmax=0.1);

Solve the model.

R = solve(emagmodel);

Obtain the mesh that is stored with the solution.

msh = R.Mesh;

Centers and Normal Vectors of Triangles Forming Northern Hemisphere

Create the triangulation object from the obtained mesh and use it to find triangular facets forming the surface of the northern hemisphere. Then find the centers and normal vectors for these triangular facets.

Create a 3-D triangulation object for the entire sphere. First, find all mesh nodes. These nodes are the triangulation points.

pts = msh.Nodes';

Obtain the connectivity list for the triangulation object by extracting the vertices of the mesh elements.

cl = msh.Elements(1:4,:)';

Remove unreferenced vertices.

ind = unique(cl(:));
pp(ind) = 1:length(ind);
cl = pp(cl);
pts = pts(ind,:);

Create a triangulation object.

volume = triangulation(cl,pts);

Extract the triangular facets of the 3-D triangulation that cover the surface of the sphere.

[cl,pts] = freeBoundary(volume);

Create a 2-D triangulation object representing the surface of the sphere.

surf = triangulation(cl,pts);

Find the centers of each triangle.

centers = incenter(surf).';

Obtain only those centers that belong to the northern hemisphere.

indx = find(centers(3,:)>=0);
centers = centers(:,indx);

Extract the triangles that belong to the surface of the northern hemisphere.

surfConnectivity_north = surf.ConnectivityList(indx,:);

Create a vector for storing the areas of each of these triangles.

area = zeros(size(surfConnectivity_north,1),1);

Find the area of each triangle.

for i = 1:size(surfConnectivity_north,1)
    a =  surf.Points(surfConnectivity_north(i,:),:);
    area(i) = 0.5*norm(cross(a(2,:)-a(1,:),a(3,:)-a(1,:)));
end

Find unit normal vectors for each of these triangles.

nMat = faceNormal(surf,indx');

Now compute the normal vectors by multiplying the unit vectors by the corresponding areas of the triangles.

nMat =  nMat.*repmat(area,[1 3]);
nMat = reshape(nMat',3,1,[]);

Plot the triangulation along with the centers and normal vectors.

trisurf(cl,pts(:,1),pts(:,2),pts(:,3), ...
     FaceColor="cyan",FaceAlpha=0.8);
axis equal
hold on  
quiver3(centers(1,:),centers(2,:),centers(3,:), ...
     nMat(1,:),nMat(2,:),nMat(3,:),Color="red");

Maxwell Stress Tensor

Generate the Maxwell stress tensor at the mesh nodes and add it to the results object, R.

R = generateMaxwellStressTensor(R);

Interpolate the Maxwell stress tensor at the midpoints of the triangles forming the northern hemisphere.

MSTIntrp = interpolateMaxwellStressTensor(R,centers);

Force Exerted by Southern Hemisphere on Northern Hemisphere

Find the force acting on the northern hemisphere from the southern hemisphere. First, calculate the product of the Maxwell stress tensor at the centers and the normal vectors for each triangle on the surface of the northern hemisphere.

Force = pagemtimes(MSTIntrp,nMat);

Reshape it to the suitable form.

Force = reshape(Force,3,[],1);

Sum the forces to approximate the total force.

ForceTotal = sum(Force,2)

ForceTotal = 3×1
1011 ×

    0.0001
    0.0063
    1.7548

Compare this result with the expected magnitude of the force [1].

F = 1/(4*pi*eps0)*Q^2/(8*Rad^2)

F = 
1.7741e+11

The approximated force points mainly in the z-direction, with the z-component being orders of magnitude larger than the x- and y-components. The first two digits also match the expected magnitude of the force. You can increase the accuracy by refining the mesh.

References

[1] Griffiths, David J. Introduction to Electrodynamics. 4th ed., Cambridge University Press, 2017.