FEM: Getting Resulting Electric Force on each body from electric field

Question

Based on this excellent answer FEM: Electric Field between two arbitrary defined shapes I can compute the electric field ef between two conducting objects.

$$ F = qE$$ Now, I tried to compute the total resultant electric force on each object (acting at it´s geometric center), by simply integrating the electric field around the boundary of the object :

So according to the Mathematica Documentation, the correct way to specify a domain, is for example:

NIntegrate[1, {x, y, z} \[Element] 

ImplicitRegion[x^2 + y^2 == 1 [And] z == 0, {x, y, z}]]

Since

Region`RegionProperty[RegionBoundary[object1], {x, y}, 
         "FastDescription"][[1]][[2]]

gives me the implicit region of the contour of object1, I can find the force as follows:

NIntegrate[
 Evaluate[ef], {x, y} \[Element] 
  ImplicitRegion[
   Region`RegionProperty[RegionBoundary[object1], {x, y}, 
      "FastDescription"][[1]][[2]], {x, y}]]

Here is the full code to compute the electric field:

Needs["NDSolve`FEM`"];
(*Define Boundaries*)
air = Rectangle[{-5, -5}, {5, 5}];
object1 = Rectangle[{-2.5, 2.5}, {2.5, 2}];
object2 = Rectangle[{-2.5, -2.5}, {2.5, -2}];
reg12 = RegionUnion[object1, object2];
reg = RegionDifference[air, reg12]

mesh = ToElementMesh[reg, MaxCellMeasure -> 0.1];
mesh["Wireframe"]

eq = Laplacian[u[x, y], {x, y}]; V1 = 1; V2 = -2;
bc = {DirichletCondition[u[x, y] == V1, 
    Region`RegionProperty[RegionBoundary[object1], {x, y}, 
       "FastDescription"][[1]][[2]]], 
   DirichletCondition[u[x, y] == V2, 
    Region`RegionProperty[RegionBoundary[object2], {x, y}, 
       "FastDescription"][[1]][[2]]]};
U = NDSolveValue[{eq == 0, bc}, u, {x, y} \[Element] mesh];

ef = -Grad[U[x, y], {x, y}];



force1 = NIntegrate[
   Evaluate[ef], {x, y} \[Element] 
    ImplicitRegion[
     Region`RegionProperty[RegionBoundary[object1], {x, y}, 
        "FastDescription"][[1]][[2]], {x, y}]];



centroid1 = RegionCentroid[object1];



StreamDensityPlot[Evaluate[ef], {x, y} \[Element] reg, 
 ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
 FrameLabel -> {x, y}, StreamStyle -> LightGray, VectorPoints -> Fine,
  PlotRange -> Automatic, 
 Epilog -> {{Arrow[{centroid1, centroid1 + force1}], 
    Text["Force", centroid1 + force1/2, Background -> LightRed]}}]

I get:

My question is: Is this approach correct ? And if not: How would you do it ?

Answer 1

Electrostatic force acting on a conductor in electric field is given by $\vec {F}=\int{\frac {E^2}{8\pi}\vec {n}dS}$. In this case we have a plate with normal vector $\vec {n}=(0,\pm 1)$ for upper and lower surfaces respectively. To calculate electric field we take air more wide and compute sufficient mesh:

 Needs["NDSolve`FEM`"];
(*Define Boundaries*)
air = Rectangle[{-10, -10}, {10, 10}];
object1 = Rectangle[{-2.5, 2.}, {2.5, 2.5}];
object2 = Rectangle[{-2.5, -2.5}, {2.5, -2}];
reg12 = RegionUnion[object1, object2];
reg = RegionDifference[air, reg12];

mesh = ToElementMesh[reg, 
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    area > 0.003 (0.1 + 10 Norm[Mean[vertices]])]]
mesh["Wireframe"]
eq = Laplacian[u[x, y], {x, y}]; V1 = 1; V2 = -2;
bc = {DirichletCondition[u[x, y] == V1, 
    Region`RegionProperty[RegionBoundary[object1], {x, y}, 
       "FastDescription"][[1]][[2]]], 
   DirichletCondition[u[x, y] == V2, 
    Region`RegionProperty[RegionBoundary[object2], {x, y}, 
       "FastDescription"][[1]][[2]]]};
U = NDSolveValue[{eq == 0, bc}, u, {x, y} \[Element] mesh];

ef = -Grad[U[x, y], {x, y}];

force1 = -NIntegrate[Evaluate[ef.ef] /. y -> 2, {x, -2.5, 2.5}];
force1 = force1 + 
  NIntegrate[Evaluate[ef.ef] /. y -> 2.5, {x, -2.5, 2.5}];force=force1/(8 Pi);
centroid1 = RegionCentroid[object1];

StreamDensityPlot[Evaluate[ef], {x, y} \[Element] reg, 
 ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
 FrameLabel -> {x, y}, StreamStyle -> LightGray, VectorPoints -> Fine,
  PlotRange -> Automatic, 
 Epilog -> {{Arrow[{centroid1, centroid1 + {0, force1}}], 
    Text["Force", centroid1 + {0, force1/2}, 
     Background -> LightRed]}}]