20.4. The electric field in a capacitor

We model electrostatics. The field variables are

symbol	unit	meaning
\(\phi\)	V	electrostatic potential
\(\rho\)	\(\frac{As}{m^3}\)	charge density
\(E\)	\(\frac{V}{m}\)	electric field
\(D\)	\(\frac{As}{m^2}\)	displacement current density

where \(V\) is Volt and \(A\) is Ampere. The equations are

\[ E = - \nabla \phi \qquad D = \varepsilon E \qquad \operatorname{div} D = \rho, \]

The material parameter \(\varepsilon\) is called dielectric parameters, and is often written as \(\epsilon = \varepsilon_r \varepsilon_0\), with \(\varepsilon_0 = 8.854 \cdot 10^{-12} \frac{As}{Vm}\) is the dielectric constant of vacuum. Sign convention follows the tradition.

Combining the equations we obtain

\[ - \operatorname{div} (\varepsilon \nabla \phi) = \rho. \]

Possible boundary conditions are

• Dirichlet boundary conditions for the potential \(\phi\)

• Neumann boundary conditions for \(D_n = \varepsilon \frac{\partial \phi}{\partial n}\)

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

def MakeMesh():
    shape = MoveTo(0,0).RectangleC(20,20) \
        .MoveTo(0,1).RectangleC(5,0.5, "el1").Reverse() \
        .MoveTo(0,-1).RectangleC(5,0.5, "el2").Reverse() \
        .Face()
    shape.edges["el.*"].vertices.hpref=1
    mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=2))
    return shape,mesh

We model the geometry and generate the mesh for a plate capacitor:

geo, mesh = MakeMesh()
Draw (geo)
Draw (mesh);

Assemble and solve the finite element system. The potential at the electrodes is prescribed as Dirichlet boundary condition (+1 Volt and -1 Volt at upper and lower electrodes).

eps0 = 8.854e-12

def SolveProblem(mesh):
    fes = H1(mesh, order=5, dirichlet="el.*")
    u,v = fes.TnT()

    gfu = GridFunction(fes)
    gfu.Interpolate( mesh.BoundaryCF( {"el1":1, "el2":-1 }), mesh.Boundaries(".*"))

    a = BilinearForm(eps0*grad(u)*grad(v)*dx).Assemble()
    inv = a.mat.Inverse(freedofs=fes.FreeDofs())
    gfu.vec.data -= inv@a.mat * gfu.vec
    return gfu

gfu = SolveProblem(mesh)

ea = { "euler_angles" : [-70,0,-40]} 
Draw (gfu, deformation=True, scale=5, **ea);

The electric field is the gradient:

Draw (Norm(grad(gfu)), mesh, order=3, deformation=True, min=0, max=3, **ea);

To resolve the singularities at the corners we apply strong mesh refinement. We use RefineHP:

geo, mesh = MakeMesh()
mesh.RefineHP(3, 0.25)
gfu = SolveProblem(mesh)
Draw (gfu, deformation=True, scale=5, **ea)
Draw (Norm(grad(gfu)), mesh, order=3, deformation=True, min=0, max=3, **ea);

To avoid such singularities, we can smooth the sharp corners:

def MakeMeshRounded():
    square = MoveTo(0,0).RectangleC(20,20).Face()
    el1 = MoveTo(0,1).RectangleC(5,0.5).Face()
    el1 += MoveTo(2.5,1).Circle(0.25).Face()
    el1 += MoveTo(-2.5,1).Circle(0.25).Face()
    el1.edges.name="el1"
    el2 = MoveTo(0,-1).RectangleC(5,0.5).Face()
    el2 += MoveTo(2.5,-1).Circle(0.25).Face()
    el2 += MoveTo(-2.5,-1).Circle(0.25).Face()
    el2.edges.name="el2"
    geo = square - el1 - el2
    mesh = Mesh(OCCGeometry(geo, dim=2).GenerateMesh(maxh=2))
    return geo,mesh

geo, mesh = MakeMeshRounded()
Draw(geo)
Draw (mesh);

gfu = SolveProblem(mesh)
Draw (gfu, deformation=True, scale=5, **ea)
Draw (Norm(grad(gfu)), mesh, order=3, deformation=True, min=0, max=2, **ea);

We still see spikes at the vertices of the mesh. To avoid such home-made singularities we can curve the elements using the Curve function.

geo, mesh = MakeMeshRounded()
mesh.Curve(5)

Draw (mesh)
gfu = SolveProblem(mesh)
Draw (gfu, deformation=True, scale=5, **ea)
Draw (Norm(grad(gfu)), mesh, order=3, deformation=True, min=0, max=3, **ea);

20.4.1. Exercises:

• Compute the total energy in the electrostatic field \(\tfrac{1}{2} \int \varepsilon | E |^2\)

• Compute the total charge at one electrode \(\int_{\Gamma_{El1}} D_n\)

• Visualize the Neumann data, the surface charge density \(\rho_s = D_n\)

• Use symmetry of the problem to reduce the computation to a quarter of the domain

• Model the capacitor in 3D

Use the following trick to compute Neumann data: Define the linear functional on \(H^1\).

\[ r(\psi) = \int \varepsilon \nabla \phi \nabla \psi - \int \rho \psi \]

By properties of the solution, it vanishes for \(\psi \in H_0^1\). Perform integration by parts backwards to obtain

\[ r(\psi) = - \int_\Omega (\operatorname{div} \varepsilon \phi + \rho) \psi + \int_{\partial \Omega} \varepsilon \partial_n \phi \psi = \int_{\partial \Omega} D_n \psi \]