# FEM Capacitor Not Converging with Parallel Plate Model

I am trying to model a 3D electrostatics finite element model of a capacitor to see if I can accurately model the behavior of a parallel plate capacitor before I alter the geometry of the electrodes and simulate a more complex scenario such as coplanar geometry.

First I define my boundary conditions in units of cm:

Needs["NDSolveFEM"];
(*Define boundaries*)
(*Volume of Mesh*)
axmin = -.9; axmax = .9; aymin = -.1; aymax = .1; azmin = -.9; azmax \
= .9;
(*Volume of dielectric*)
dxmin = -.9; dxmax = .9; dymin = -.09; dymax = .09; dzmin = -.9; \
dzmax = .9;
(*Electrode 1 geometry*)
e1xmin = -.9; e1xmax = .9; e1ymin = -.1; e1ymax = -.09; e1zmin = -.9; \
e1zmax = .9;
(*Electrode 2 geometry*)
e2xmin = -.9; e2xmax = .9; e2ymin = .09; e2ymax = .1; e2zmin = -.9; \
e2zmax = .9;

air = Cuboid[{axmin, aymin, azmin}, {axmax, aymax, azmax}];
object1 = Cuboid[{e1xmin, e1ymin, e1zmin}, {e1xmax, e1ymax, e1zmax}];
object2 = Cuboid[{e2xmin, e2ymin, e2zmin}, {e2xmax, e2ymax, e2zmax}];

reg12 = RegionUnion[object1, object2];
reg = RegionDifference[air, reg12];
(*Create mesh*)
bmesh = ToBoundaryMesh[reg, "MaxBoundaryCellMeasure" -> 0.0001];
mesh = ToElementMesh[bmesh, "MeshElementType" -> TetrahedronElement];
mesh["Wireframe"]


Then I set up the permittivity distribution using a piecewise function and Poisson's equation and solve using NDSolveValue:

eps0 = 8.85*10^(-14); (*F per cm*)
eps1 = 4;
eps[x_, y_, z_] :=
Piecewise[{{eps1*eps0,
dxmin <= x <= dxmax && dymin <= y <= dymax &&
dzmin <= z <= dzmax}}, eps0];
eq = Div[eps[x, y, z]*Grad[u[x, y, z], {x, y, z}], {x, y, z}] == 0;
V1 = 1; V2 = 0;
(*Boundary Conditions*)
bc = {DirichletCondition[u[x, y, z] == V1,
RegionRegionProperty[RegionBoundary[object1], {x, y, z},
"FastDescription"][[1]][[2]]],
DirichletCondition[u[x, y, z] == V2,
RegionRegionProperty[RegionBoundary[object2], {x, y, z},
"FastDescription"][[1]][[2]]]};
(*Solve Poission's Equation*)
U = NDSolveValue[{eq, bc}, u, {x, y, z} \[Element] mesh];


And then finally I calculate the total charge Q and find C:

Q = Quantity[
Total[NIntegrate[-eps[x, y, z] Grad[U[x, y, z], {x, y, z}], {x, y,
z} \[Element] mesh]], "Coulombs"];
cap = UnitConvert[Abs[Q/Quantity[V1 - V2, "Volts"]], "Femtofarads"]


Running this instance gives a Capacitance of 1147 fF while calculating the capacitance using the parallel plate model using the area of the electrodes and the distance between them:

a = Quantity[(Abs[e1xmin] + e1xmax)*(Abs[e1zmin] + e1zmax), (
"Centimeters")^2];
d = Quantity[Abs[e1ymax] + e2ymin, "Centimeters"];
UnitConvert[


gives a capacitance of 6372 fF. I am in the same order of magnitude, and maybe that is all I could hope for but I would think the results would converge as the plates move closer together. Moving the plates closer together gives:

Plotting the capacitance vs separation distance gives:

As you can see, the FEM model shows an increase in capacitance as separation increases, which should not be the case. I have tried decreasing the MaxBoundaryCellMeasure which did not have any effect. I am able to visualize both the potential gradient and electric field vectors using SliceContourPlot3D and VectorPlot3D, respectively, and I don't see any issues. I have limited the simulation volume to the edges of each electrode which should mitigate any edge effects but maybe I am wrong about that. If there are any glaring errors in my code or mathematical faux pas please point them out to me. I am a chemist and have not taken any electrostatics classes since undergrad. Below are the SliceContourPlot3D and VectorPlot3D graphics for the code shown:

In the parallel plate capacitor one plate has a positive electric charge $$+q$$, and another plate has a negative electric charge $$-q$$, therefore the total electric charge is zero. To avoid this zero effect we use absolute value $$q$$ to define $$C$$ as follows

Needs["NDSolveFEM"];
(*Define boundaries*)
(*Volume of Mesh*)
axmin = -.9; axmax = .9; aymin = -.1; aymax = .1; azmin = -.9; azmax \
= .9;
(*Volume of dielectric*)
dxmin = -.9; dxmax = .9; dymin = -.09; dymax = .09; dzmin = -.9; \
dzmax = .9;
(*Electrode 1 geometry*)
e1xmin = -.9; e1xmax = .9; e1ymin = -.1; e1ymax = -.09; e1zmin = -.9; \
e1zmax = .9;
(*Electrode 2 geometry*)
e2xmin = -.9; e2xmax = .9; e2ymin = .09; e2ymax = .1; e2zmin = -.9; \
e2zmax = .9;

air = Cuboid[{axmin, aymin, azmin}, {axmax, aymax, azmax}];
object1 = Cuboid[{e1xmin, e1ymin, e1zmin}, {e1xmax, e1ymax, e1zmax}];
object2 = Cuboid[{e2xmin, e2ymin, e2zmin}, {e2xmax, e2ymax, e2zmax}];

reg12 = RegionUnion[object1, object2];
reg = RegionDifference[air, reg12];
(*Create mesh*)
bmesh = ToBoundaryMesh[reg, "MaxBoundaryCellMeasure" -> 0.0001];
mesh = ToElementMesh[bmesh, "MeshElementType" -> TetrahedronElement]

eps0 = 8.85*10^(-14); (*F per cm*)
eps1 = 4;
eps[x_, y_, z_] :=
Piecewise[{{eps1,
dxmin <= x <= dxmax && dymin <= y <= dymax &&
dzmin <= z <= dzmax}, {1, True}}];
eq = Div[eps[x, y, z]*Grad[u[x, y, z], {x, y, z}], {x, y, z}] == 0;
V1 = 1; V2 = 0;
(*Boundary Conditions*)
bc = {DirichletCondition[u[x, y, z] == V1,
RegionRegionProperty[RegionBoundary[object1], {x, y, z},
"FastDescription"][[1]][[2]]],
DirichletCondition[u[x, y, z] == V2,
RegionRegionProperty[RegionBoundary[object2], {x, y, z},
"FastDescription"][[1]][[2]]]};
(*Solve Poission's Equation*)
U = NDSolveValue[{eq, bc}, u, {x, y, z} \[Element] mesh];


Visualization

DensityPlot[U[x, y, 0], {x, dxmin, dxmax}, {y, dymin, dymax},
PlotLegends -> Automatic, ColorFunction -> "Rainbow",
AspectRatio -> Automatic]


Total absolute electric charge is $$2 q=$$

-  eps0 eps1 NIntegrate[
Derivative[0, 1, 0][U][x, dymin, z] +
Derivative[0, 1, 0][U][x, dymax, z], {x, dxmin, dxmax}, {z, dzmin,
dzmax}, AccuracyGoal -> 5]

Out[]= 1.2744*10^-11


Hence, $$C=q/(V1-V2)= 6.372*10^{-12}$$. From the other side we have for the parallel plate capacitor

C = eps0 eps1 (dxmax - dxmin) (dzmax - dzmin)/(dymax - dymin)

Out[]= 6.372*10^-12

• Thank you for this, I think I am getting the picture here. If I were to change the geometry such that the dielectric does not occupy all of the space between each electrode, how might the situation change? Commented Apr 13, 2023 at 16:26
• @DallonPenney The capacitance still increases, while it should be decreasing. There are two major errors in your code, 1 the total charge formular is incorrect and 2 you're integrating over dielectric material where no charge is present (yet your code computes non-zero charge, because the formula isn't correct). One more thing, unless we're integrating at an enclosing boundary via divergence theorem, we need to think about how we should deal with inner volume of each plate, because FEM might have huge error at boundary interpolation. I'll share a bit more in the answer I'm going to write. Commented Apr 13, 2023 at 16:35
• @DallonPenney It depends on geometry. If you show your new geometry then we can try to compute capacitance. Commented Apr 13, 2023 at 16:45
• @ShinKim I'm interested to see your answer! Like I mentioned in the question, my ultimate goal is to model a coplanar capacitor where there is some 3D permittivity distribution representative of a dielectric material somewhere between both coplanar plates Commented Apr 13, 2023 at 20:50

When modeling capacitors, there are a few things that need to be considered. First, one must decide whether to model the volume and the dielectric region or just the dielectric.

By using the first approach one can simulate the field that appears at the edge of the plates and which extends away from the device inside the volume considered, known as the fringing field. The second approach, the one used here, just considers the field inside the dielectric region.

The first approach will need infinite domains to do it the right way. If there are not available, one should use specific boundary conditions at the boundaries of the surrounding domain. One will need to make bigger the surrounding volume if one wants to get a more accurate value. By making it bigger, the capacitance value will start to converge to a specific value.

One should expect that the capacitance value obtained from just modeling the dielectric will be lower because it does not consider the fringing fields.

This is the workflow that I used.

The mesh is the same as yours. If you are considering modeling different materials I could recommend you to use ElementMarkers.

Just to simplify things, I avoid using the Piecewise function, as were are just modeling the dielectric region.

eps0 = 8.85*10^(-14); (*F per cm*)
eps1 = 4;
eps = eps1*eps0;


Define variables and equation:

vars = {V[x, y, z], {x, y, z}};
eq = DiffusionPDETerm[vars, eps*IdentityMatrix[3]];


I also simplified the preds here, if you are going to model more complex geometries you could use ElementMarkers also.

V1 = 1; V2 = 0;
(*Boundary Conditions*)
bc = {DirichletCondition[V[x, y, z] == V1, y == dymin],
DirichletCondition[V[x, y, z] == V2, y == dymax]};


Solve the equation:

Vfun = NDSolveValue[{eq == 0, bc}, V, {x, y, z} \[Element] mesh];


Then, we compute the electric field $$E$$ and the electric displacement $$D$$:

Ex[v__] := -Derivative[1, 0, 0][Vfun][v];
Ey[v__] := -Derivative[0, 1, 0][Vfun][v];
Ez[v__] := -Derivative[0, 0, 1][Vfun][v];

Dx[v__] := eps*Ex[v];
Dy[v__] := eps*Ey[v];
Dz[v__] := eps*Ez[v];


The equation used for computing the capacitance comes from the fact that Gauss's law implies that: $$D_n =\sigma_s$$ where $$D_n$$ is the normal component of the displacement field $$D$$, which in this case is y-component $$D_y$$. Then, we know that to get the charge from the surface charge density $$\sigma_s$$ we need to do an integral surface:

$$Q=\int_S \sigma_s dS$$ where $$S$$ is the surface, which in this case is one of the electrodes. Here we use Rectangle[] to represent this surface.

Now that we have the charge $$Q$$ we can compute the capacitance using the following formula: $$C=Q/ \Delta V$$

In[183]:= C0 =
eps*NIntegrate[
Ey[x, dymin, z], {x, z} \[Element]
Rectangle[{dxmin, dzmin}, {dxmax, dzmax}]]/(V1 - V2)

Out[183]= 6.372*10^-12


The other alternative is to compute the electrostatic field energy given by:

$$W=\frac{1}{2}\int \int \int_V D \cdot E~dV$$

energy =
NIntegrate[
Sqrt[(Dx[x, y, z])^2 + (Dy[x, y, z])^2 + (Dz[x, y, z])^2]*
Sqrt[(Ex[x, y, z])^2 + (Ey[x, y, z])^2 + (Ez[x, y, z])^2]*0.5, {x,
y, z} \[Element] mesh]



Then, we know that the energy of a capacitor is given by:

$$W=C(\Delta V)^2/2$$

From the last equation, we can know the capacitance value is:

In[181]:= 2*energy/(V1 - V2)^2

Out[181]= 6.372*10^-12