10
$\begingroup$

I was wondering how to do the following: I would like to compute the electrostatic field between two shapes using the FEM method. 

(*Define Boundaries*)
air = Rectangle[{-3, -3}, {3, 3}];
object1 = Disk[];
object2 = Rectangle[{2, 0}, {2.5, 2}];
Show[Graphics[{Blue, air}], Graphics[{Magenta, object1}],Graphics[{Green, object2}]]

enter image description here

Calculation of the electric field at every point {x,y} in 2D space:

enter image description here $r_i$ is the vector of the point charge; $r$ is the vector to the point in 2D (or also 3D) space where we want to calculate the electric field.

I make a Mathematica function out of it (for the moment I omit the constant term):

eField[x_, y_] := q Sum[({x, y} - pts[[i]])/Norm[{x, y} - pts[[i]]]^3, {i, n}]

where pts[[i]] are the boundary points of the charged object and x and y are coordinates of the "air" object.

How I would proceed:

  1. I calculate the electrostatic field of object 1 -> $E_1$

  2. I calculate the electrostatic field of object 1 -> $E_2$

  3. I use superposition to get the resultant electric field: $E_{Total} = E_1 +E_2$

    Needs["NDSolve`FEM`"];

r1 = RegionDifference[air, object1];
r2 = RegionDifference[air, object2];
mesh1 = ToElementMesh[r1];
mesh2 = ToElementMesh[r2];
mesh1["Wireframe"]
mesh2["Wireframe"]

enter image description here

I would really appreciate if someone could show me how to do it in Mathematica using finite elements (FEM).

EDIT: Basend on the excellent answer below, I would like to use the answer here Get Perimeter Region from Object to automate the finding of the region boundaries for the DirichletCondition: 

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

mesh = ToElementMesh[reg, {{-5, 5}, {-5, 5}}, 
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    area > 0.001 (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}];

DensityPlot[U[x, y], {x, y} \[Element] reg, 
 ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
 FrameLabel -> Automatic, PlotPoints -> 50, 
 PlotRange -> {{-4, 4}, {-4, 4}}]

StreamDensityPlot[Evaluate[ef], {x, y} \[Element] reg, 
 ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
 FrameLabel -> {x, y}, StreamStyle -> LightGray, VectorPoints -> Fine,
  PlotRange -> {{-3, 3}, {-2.5, 3}}]

enter image description here enter image description here

EDIT 2: Just for beauty: Parallel Plate Capacitor Use this answer to make it work: FEM Simulation: Meshing two Arbitrary objects in an "air" mesh

enter image description here

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Have you seen this or this $\endgroup$
    – user21
    Commented Apr 3, 2020 at 8:59
  • 2
    $\begingroup$ @james The solution depends on the electrical conductivity of the disk and the rectangle. There is one solution for a conductor, and another for a dielectric. $\endgroup$
    – Alex Trounev
    Commented Apr 3, 2020 at 10:33
  • $\begingroup$ The comment of @AlexTrounev is very important. Are the disk and the rectangle conductor ? $\endgroup$
    – andre314
    Commented Apr 3, 2020 at 10:41
  • 1
    $\begingroup$ @AlexTrounev here on MSE, there is a recurrent problem of confusion between 1) dielectric versus conductor 2) potential imposed at boundary versus presence of charges. I give up, though if you want to explain the things, I could do an effort to participate. The problem is not related to Mma but some people think the Mathematica+FEM will solve all this magically. $\endgroup$
    – andre314
    Commented Apr 3, 2020 at 10:54
  • $\begingroup$ One point more : The formula above giving the electrical field E=..r/r^3 is wrong in the 2D case. In the 2D case it is E=..r/r^2. I have already seen this error on MSE too. $\endgroup$
    – andre314
    Commented Apr 5, 2020 at 13:22

1 Answer 1

Reset to default
14
$\begingroup$

In the case of two metal objects, we can set the potential of each object as $V_1, V_2$. Then the code for a numerical solution in 2D is

Needs["NDSolve`FEM`"];
(*Define Boundaries*)air = Rectangle[{-5, -5}, {5, 5}];
object1 = Disk[];
object2 = Rectangle[{2, 0}, {2.5, 2}]; reg12 = 
 RegionUnion[object1, object2];
reg = RegionDifference[air, reg12];
mesh = ToElementMesh[reg, 
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    area > 0.001 (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, x^2 + y^2 == 1], 
   DirichletCondition[
    u[x, y] == 
     V2, (x == 2 || x == 2.5 && 0 <= y <= 2) || (y == 0 || 
       y == 2 && 2 <= x <= 2.5)]};
U = NDSolveValue[{eq == 0, bc}, u, {x, y} ∈ mesh];

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

Visualisation of solution 

{DensityPlot[U[x, y], {x, y} ∈ reg, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> Automatic, PlotPoints -> 50, 
  PlotRange -> {{-4, 4}, {-4, 4}}], 
 StreamDensityPlot[Evaluate[ef], {x, y} ∈ reg, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {x, y}, StreamStyle -> LightGray, 
  VectorPoints -> Fine, PlotRange -> {{-1, 3}, {-1, 3}}]}

Figure 1

Update 1. Next code is devoted to solve electrostatic problem for combination of dielectric and conducting objects (glass cylinder and metal strip). For dielectric we put electric charge $q_1$, and for metal we put potential $V_2$. Code:

Needs["NDSolve`FEM`"];
par = {eps1 -> 3.5, eps2 -> 1.0}; air = 
 Rectangle[{-5, -5}, {5, 5}];
object1 = Disk[]; q1 = 1; vol1 = 
 NIntegrate[1, {x, y} ∈ object1]; rho1 = q1/vol1;
object2 = Rectangle[{2, 0}, {2.5, 2}]; 
rho[x_, y_] := rho1 Boole[{x, y} ∈ object1];
eps[x_, y_] := 
 eps2 + (eps1 - eps2) Boole[{x, y} ∈ object1]; reg = 
 RegionDifference[air, object2];
mesh = ToElementMesh[reg, 
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    area > 0.001 (0.1 + 10 Norm[Mean[vertices]])]]
mesh["Wireframe"]
 V2 = -2; eq = 
 Inactive[Div][
   eps[x, y] Inactive[Grad][u[x, y], {x, y}], {x, y}] == -2 Pi rho[x, 
    y]; bc = 
 DirichletCondition[u[x, y] == V2, {x, y} ∈ object2];
U = NDSolveValue[{eq /. par, bc}, u, {x, y} ∈ mesh];

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

Visualisation 

{DensityPlot[U[x, y], {x, y} ∈ mesh, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> Automatic, PlotPoints -> 50, 
  PlotRange -> {{-4, 4}, {-4, 4}}], 
 StreamDensityPlot[Evaluate[ef], {x, y} ∈ reg, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {x, y}, StreamStyle -> LightGray, 
  VectorPoints -> Fine, PlotRange -> {{-1, 3}, {-1, 3}}], 
 StreamDensityPlot[Evaluate[ef], {x, y} ∈ reg, 
  ColorFunction -> Hue, FrameLabel -> {x, y}, StreamStyle -> Blue, 
  PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}}, PlotLegends -> Automatic]}

Figure 2

Improve this answer
$\endgroup$
2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.