Electrostatics
| Introduction | Boundary conditions in electrostatics |
| Overview example | Nomenclature |
| Electrostatic equation | References |
| Secondary quantities |
Introduction
This monograph uses partial differential equations to model and analyze time independent electric fields, electric potentials, and charge densities distributions. The main assumption in an electrostatic simulation is that the electric and magnetic fields are not coupled.
The electrostatics equation is used to model electric fields produced by stationary charges. As a simplification of Maxwell's equations it is only used to describe insulating, or dielectric materials.
Modeling electromagnetic devices with partial differential equations (PDEs) is not the only way to model electromagnetic devices. Other techniques include setting up ordinary differential equations (ODEs). This approach is followed by the Wolfram System Modeler. Roughly speaking the system modeler approach is more suitable for large systems of electromagnetic devices interacting, while the partial differential equation approach is more suitable for a fine grained analysis of a specific device. In some cases it is beneficial to use a combination of the two approaches.
The approach taken here is that in the introductory section a capacitor is used to introduce an electrostatic analysis and the functionality available. This will be followed by a more theoretical explanation of the underlying ideas and concepts. The theoretical background is much easier to understand once an intuition for the electrostatic analysis exist. After that the available boundary conditions are discussed.
The equation used in this monograph is the electrostatic equation:
The goal of this electric device analysis is to find the potential distribution 
in units of volts 
under specific constraints. A subsequent step then finds secondary fields, such as the electric field intensity 
, and values such as the capacitance 
. The analysis and interpretation of these physical quantities are useful to create a better quality engineering design of the device under consideration.
The modeling process as such results in a system of partial differential equations (PDEs) that can be solved with NDSolve.
Extended examples of electromagnetic modeling can be found in the Model Collection.
The electromagnetic device analysis is typically done in stages. First, for the object to be analyzed a geometric model needs to be created. The geometric model is typically created within a computer aided design (CAD) process. CAD models can either be imported or created in product. To import geometries common file formats like DXF, STL or STEP are supported. These geometries can be imported with Import. The alternative is to create the geometrical models in product with, for example, by using OpenCascadeLink.
Once the geometric model is made available some thought needs to be put into what what type of analysis is to be performed. The next step is the setup of boundary conditions and constraints. Materials to be used further specify the PDE model. Once the PDE model is fully specified, the subsequent finite element analysis will then compute the desired quantities of the device under investigation. These quantities are then post processed. Either by visualizing them or some derived quantities are computed. This notebook shows the necessary steps for everything except the CAD model generation.
Electromagnetic devices are typically three dimensional objects. For the models, special cases exist that result in simplified 2D models. In fact, 2D models have been known to solve 90% of the important problems in electrical engineering [Cardoso, 2018].
The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.
Overview example
To illustrate the usage of the finite element method in these subcases of electromagnetics it is instructive to present a simple example and give an overview of the setup, various analysis types and post processing steps possible.
The first example is to introduce the workflow of setting up a simple electromagnetics PDE model, which will cover an electrostatic analysis.
Consider the following setup where we have a 3D circular parallel capacitor, which is composed of two conductive plates with a dielectric in between.
Depiction of a 3D capacitor and the cross-section of the dielectric region, in blue, it's height 
and radius 
.
To keep the model simple, we are only going to model the dielectric region, in blue.
Creating an electromagnetics model always comprises the same steps:
Geometry
When a capacitor devices is modelled there are two possible ways to model it. One is to consider the device and the surrounding volume. In this approach one can simulate the field that appears at the edge of the plates and extends away outside from the device into the surrounding volume considered. This field is known as the fringing field. The second approach, used here, considers the dielectric region. For more details on which approach to use and what are the effects on the model, refer to the capacitance section.
One thing to keep in mind is the scale used in the geometric model. If the length of the boundary mesh is, for example, in units of meters then the material parameters will need to be specified in consistent units.
More information on generating or importing 3D geometric models can be found in the Using OpenCascadeLink tutorial.
Material parameters
The next step is to assign material parameters. Generally, all parameters for a electromagnetics model are collected in an Association pars that includes the necessary parameter values.
The default vacuum permittivity value 
is used throughout the whole monograph. The "VacuumPermittivity" parameter can also be specified if the constant should have a different value, like 1, to scale the equation.
In this model, the dielectric of the capacitor will have a relative permittivity of 
. This value is just representing a general dielectric material.
Specifying a material as an Entity directly is currently not possible.
Material parameters where the values are given as Quantity objects can be used. The exact property names to specify material properties needed can be found on the reference pages of ElectrostaticPDEComponent.
Geometries consisting of several material can also be used and an example of such a use case is presented in the section on multiple materials of electrostatics.
Units
Should the units of the geometry be different from the material units, then the material units can be scaled.
Internally, all material data units are converted to "SIBase" units. As a consequence the default unit of length is "Meters". If the units of the geometry are also in meters then nothing needs to be changed. If the units of the geometry are not in meters then either the PDE and material properties need to be scaled to the units the geometry or the geometry needs to be scaled to "Meters". To scale the units of the PDE and material parameters the parameter "ScaleUnits" can be given. If not explicitly stated otherwise, examples in this notebook use the default "SIBase" units.
Boundary conditions
Boundary conditions must be defined at physical boundaries to fully specify the problem. In this case, the device will be charged by applying a voltage difference between the plates. In essence, this is done by applying a DirichletCondition or NeumannValue at the terminals respectively. Various boundary conditions can be used and will be discussed in more detail in the Boundary conditions in electrostatics section. For the purpose of this overview an electric potential condition will be sufficient.
The purpose of this section is to establish the positions where the boundary conditions are to be applied.
A way to find the positions where the boundary conditions are applied is to visualize them together with the outline geometry defined previously in the geometry section.
The same applies for the lower plate.
On the side of the dielectric, we have a zero NeumannValue. A zero NeumannValue boundary condition is the default boundary condition if nothing else is specified, and we can thus omit that boundary condition.
For more complicated geometries a different technique to specify boundary condition predicates may be appropriate and it can be seen in the Spherical Capacitor tutorial.
Mesh generation
To perform a finite element analysis, the geometric region needs to be discretized into a mesh.
More information about the mesh generation process can be found in the ElementMesh generation tutorial. Another option is to import a mesh. Some common mesh file formats can be imported with the help of the FEMAddOns.
Electrostatic analysis
An electrostatic analysis will seek the electric potential and electric field distribution inside the dielectric due to the static charge distribution imposed on the model. The dependent variable is the electric potential 
.
An electrostatic analysis in effect solves:
where 
is the vacuum permittivity and 
is the relative permittivity. This equation is provided by the function ElectrostaticPDEComponent.
By specifying the relative permittivity 
, the function generates the equation based on this linear constitutive relation. Others parameters, such as a polarization vector 
, can be specified.
The result is an InterpolatingFunction object, which gives the electric potential distribution 
.
More information about the solution process and its options can be found in the NDSolve Options for Finite Elements tutorial.
Post processing
The primary solution of the electrostatics PDE model is the electric potential.
Also of interest are the electric field intensity 
and the electric flux density 
from which the capacitance value 
can be obtained. These fields are recovered from the electric potential. More details about how to compute capacitance values are shown in the capacitance section.
For a 3 dimensional model, the function Grad return a list of three InterpolatingFunction with the three independent variables 
, 
and 
each.
Various components of the electric field can be accessed by using Part.
With the electric field intensity 
computed, we can find the capacitance of the capacitor:
where 
in units of 
is the charge's magnitude held on one of the plates of the capacitor, 
is the potential difference between the plates, and 
is the normal component of the electric flux density at one of the plates.
One must first compute the voltage at the upper plate to compute the potential difference between the plates. The upper voltage is obtained by integrating 
over the disk area and dividing it by the same area.
By definition, at ground we have a potential of 0.
In this case, the inward normal component at the upper plate is the 
-component in the negative direction 
. The 
-component is given by: 
.
To calculate 
it is necessary to extract the vacuum permittivity value used in the model. To make sure that the same values are used throughout a computation, we use the default model parameter computation.
For this example an analytical solution of the value of the expected capacitance is available.
The numerical capacitance obtained is in good agreement with the analytical one.
In the above case we made use of the symbolic region representation Disk as the upper part of the cylindrical region 
. In real world applications a symbolic representation of the boundary in question may not always be readily available. Even if a boundary representation is available it may differ from the numeric representation used during the computation. This is because a mesh is only an approximation to the symbolic representation like the Disk. For this reason, it may in general be beneficial to directly use the numerical approximation of the region in question. This can be done by using FEMNBoundaryIntegrate.
When using FEMNBoundaryIntegrate one needs to specify the element marker of the desired boundary surface. In this case , the number 2 refers to the upper plate boundary. For more information on markers in a mesh see the section Markers in the element mesh generation tutorial.
While the error with the numerical mesh is larger then with the symbolic surface region it is still acceptable. The use FEMNBoundaryIntegrate is convenient because it is a generally usable workflow when no symbolic representation of the surface in question is available.
As yet another alternative method, one can compute the capacitance by means of the electrostatic energy 
:
where 
is given by a volume integral:
Here, we have used the mesh of the region 
.
Details about the methods employed to compute the capacitance are shown in the capacitance section.
Electrostatic equation
This section gives an introduction to modeling electrostatics models.
Electrostatics is the subfield of electromagnetics that deals with electric fields produced by static charges. So, the driving force behind the generated electric field are distributed charges or a specified potential difference 
in units of volts 
. The electrostatic equation which is derived from the Maxwell's equations, is used for modeling steady electric fields produced by static charges within an insulating or dielectric material. Here we show the most general form of the equation within a dielectric medium:
The dependent variable in this equation is the electric potential 
, which varies with position 
. The equation is made up of a diffusive term 
with the vacuum permittivity 
as the diffusion coefficient, and of a derivative term 
in which the variable 
is the polarization vector. The term 
in units of 
denotes a volume charge density within the domain, and it is explained in the Source Types section.
Note that metallic components in the domain can not be considered with this equation; this equation is for dielectric or insulating material only because for conducting material there is no potential difference. This means that electrodes, for example, have to be modeled as boundary conditions. Also, any metal inside the domain needs to be removed.
The electrostatic equation can be solved for 1D axisymmetric, 1D, 2D, 2D axisymmetric models, and 3D models.
Model setup
The inputs needed for an electrostatic model are:
Next, we give some examples for the equation setup.
For 2D models there exists a special case of symmetry where the electric potential varies only in the 
- and 
-directions and is constant in the 
-direction. With this symmetry a similar equation is solved with the exception that a thickness variable 
denoting thickness in the 
-direction is added:
The default value for the thickness is 
.
Just as there can be a symmetry when simulating in the 
-plane, a symmetry can exist in 1D models. In a 1D electrostatic model, a cross-sectional area 
, which represents the 
- and 
-direction, is added to the equation:
The default value for the cross-sectional area is 
.
The 1D and 2D axisymmetric models take into account another type of symmetry that is observed in 3D solids, which have an axis of revolution. For more details on these cases, see section axisymmetric models.
Note that this model definition uses inactive PDE operators. "Numerical Solution of Partial Differential Equations" has several sections that explain the use of inactive operators.
Electrostatics
As the name electrostatic suggests, this subcase deals with electric fields that do not vary in time, and thus the fields are called static. To simplify things we will start by assuming vacuum as the simulation domain.
An electrostatics case will be described by two of Maxwells's equations: Gauss's law:
and Faraday's law without the magnetic induction term:
Gauss's law (Eq. 1) states that the divergence of the electric field intensity 
at any point in space is proportional to the volume charge density 
at that point. In other words, the charge density influences the behaviour of 
in that region. The divergence for the electric field 
will tell us how much the electric field is spreading out or converging at any point in space.
On the other hand Faraday's law (Eq. 2) states that 
is curl free and therefore a conservative vector field. This equation allows us to derive the scalar electric potential 
.
Using the properties of field operators we can deduce the following equation:
which states that the curl of the gradient of any scalar function 
is zero. Thus comparing both equations and substituting the electric field 
in Eq. 3, we can deduce that 
can be represented as the gradient of a function 
.
The negative sign was placed by convenience and for historical reasons.
Now, by substituting Eq. 4 in Eq. 5 one can obtain the electrostatic equation in free space:
The above equation is valid in vacuum. The electrostatic equation in vacuum can be extended to consider dielectric materials. Dielectric materials are distinguished from conductors by not having any free electrons, instead, charges in dielectric material are bounded by forces. When an external electric field is applied to a dielectric material, dipole moments are induced because the electron cloud in the atoms are displaced in the opposite direction of the 
field and the nuclei are displaced from their equilibrium positions in the same direction as the 
filed points to.
This displacement of charges creates a bound volume charge density 
and in consequence a vector field 
, which expresses the density of permanent or induced electric dipole moments. Due to this polarization vector, the electric field inside a dielectric material will be different from the one in free space [Sadiku, 2011].
On the left: a nonpolar atom. To the right: polarization of a nonpolar atom when an electric field 
is applied.
The polarization vector 
and the bound volume charge density 
are related by the following equation:
which resembles Gauss's law. The negative sign arises due to the opposite signs on the charges in the aligned dipoles of the dielectric, which are induced by the external 
field.
Now, we would like to know how these bound charges of a dielectric interact with the external electric field that is applied to the dielectric. It is helpful to distinguish between the free charges and the bound charges that make up the dielectric itself. Free charges denoted as 
can be either introduced to or taken away from a conductor, like a plate of a capacitor. On the other hand, bound charges are not mobile in dielectric.
So, in a system where we have both free and bound density charges, one needs to modify the Gauss's law to include the total charge density:
By substituting Eq. 6 into Eq. 7, we have:
which can be further simplified by introducing a new quantity 
, which name is the electric flux density:
We define the electric flux density 
by:
This constitutive relation basically says that the application of an external electric field 
to a dielectric material causes the flux density to be greater than it would be in free space. In free space 
.
If we substitute Eq. 8 into Eq. 9, then we get the electrostatic equation:
For linear materials, 
can be described in terms of the electric field as 
, where the unit-less 
is the electric susceptibility of the material. So, for linear materials, the constitutive relations is:
where the unit-less 
is the relative permittivity given by:
and 
is the absolute permittivity. Note that for free space 
. On the other hand, a material is said to be nonlinear if 
does not varies linearly with 
.
The relative permittivity can depend on the space coordinates 
and vary in the region considered, i.e an inhomogeneous material, or vary with direction, i.e. an anisotropic material.
It is common to model different materials in the simulation domain. For example, if one wants to simulate the electric field that extends beyond a capacitor into the surrounding space. This field is known as the fringing field. To simulate the fringing field of a capacitor, one needs to include the dielectric region and the air surrounding it.
Here we model a cross-section of a parallel capacitor, assuming that the electrodes extend both ways along the 
-axis, which is out of the 
plane. This model will be helpful to illustrate energy, capacitance, and force calculations in future sections.
Cross-section of a parallel capacitor surrounded by vacuum.
Also, in this example, we will employ element markers to specify the relative permittivity value for each subregion, and to specify the boundary conditions. Element marker and their use in meshes is explained in detail in the section Element Marker in the Element Mesh Generation tutorial.
of the box, the width 
and thickness 
of the electrodes, and half the distance 
between them:Next, we specify region markers with the "RegionMarker" option of ToElementMesh. To do so a coordinate within the dielectric and the air region should be given, as well as an integer marker. Optionally, an additional maximum cell measure can be specified to refine a subregion. The most accurate and efficient method to deal with multiple material regions is by element makers as then the mesh will have a specific subregion for the material which will result in an accurate solution.
Remember that metals should be removed from the mesh, when they are inside the domain of analysis.
The air region has a relative permittivity of 
while the dielectric region has a value of 
.
See this note about the set up of efficient PDE coefficients.
In this capacitor the upper electrode will have an electric potential of 
, and the lower electrode will be grounded. In addition to these boundary conditions, a zero surface charge density is applied to the box boundaries, this will constraint the electric field to remain inside the box.
A zero surface charge density is equivalent to a Neumann zero value, and as these are the natural default boundary condition, these can be left out.
In this model, markers are used in the boundary condition specification. As these are DirichletCondition boundaries, point element markers are used.
The element marker for the lower and upper electrodes are 
and 
respectively.
Next, we are going to compute the electric field intensity 
and the electric flux density 
, to see their behaviour.
To correctly compute the electric field 
and 
, one has to take into account that at interfaces between different materials there can be discontinuous in one of the components. To see details about the boundary conditions at interfaces, refer to the section Boundary Conditions in Electrostatics.
To address this discontinuity we make use of DiscontinousInterpolatingFunction or EvaluateOnElementMesh. These functions will prioritize the value from one of the subregions. The default uses the order given in "MeshElementMarkerUnion". To specify a priority you can give an option "MarkerPriority".
At the dielectric / air interface, the default, in this case, will prioritize the dielectric region of the capacitor over the air. In this example, a discontinuity in 
can be seen at the interface.
Since the electric field intensity 
is the gradient of the electric potential 
, a DiscontinuousInterpolatingFunction will be used first in the electric potential function and then the derivative will be computed from the new discontinuous function.
The DiscontinuousInterpolatingFunction is very useful when one tries to compute values at the interface between two different materials, and then used them for plots or as arguments for other computations.
The output is a DiscontinuousInterpolatingFunction, where it first shows the domain bounds and then the marker priority.
The vector plots visualizes that the zero Neumann values at the boundaries restrain the field lines to remain within the domain. In other words, the field is forced to be tangential to the exterior boundaries. This constraint in the field has a direct influence in the central part of the field and consequently in the capacitance value that is computed and discussed in the capacitance section.
Since the model has two different values for the relative permittivity 
, we must make use of EvaluateOnElementMesh to correctly calculate the electric flux density 
in the different material regions. The electric flux also will be represented with a DiscontinousInterpolatingFunction.
EvaluateOnElementMesh returns a list of discontinuous interpolating functions. These interpolating function do not have the independent variables 
and 
attached. Adding them is trivial with the Through command.
To see the electric flux discontinuity we can create a DensityPlot of the electric flux density magnitude. Sometimes, it is necessary to extract the maximum and minimum values of the plotted variable. So, first, we will create a rule allowing us to extract the "ValuesOnGrid" of the fields and then apply MinMax.
Dielectrics can be either isotropic or anisotropic. Isotropic materials are materials that have the same electric properties in all directions. For an isotropic 
, 
and 
always have the same direction due to the relation 
. In anisotropic materials the relative permittivity becomes a tensor 
and the electric properties change with direction. Crystalline materials are, for example, anisotropic [Sadiku, 2011].
For anisotropic materials 
or 
are now tensors. For a 3D case, the relative permittivity is a three-by-three the tensor that has 9 components:
where 
is the relative permittivity tensor. 
and 
are called the principal permittivity coefficients and off-diagonal permittivity coefficients, respectively.
With a relative permittivity as tensor, the relation between 
and 
changes to 
. This has the consequence that 
and 
are no longer parallel.
In the following example an electro-optic modulator (EOM) is modeled. This EOM consists of a lithium niobate core (
), a silica cladding (
), and two gold electrodes, which work as capacitor, that induce an electric field in the core. EOMs are used widely in many different photonic devices.
The lithium niobate core is an optical waveguide where the light passes through. This material is an orthotropic material. The material is symmetric along the principal directions 
and 
, and the off-diagonal permittivity coefficients are zero. The permittivity tensor for this case is:
The design of this EOM consists of a trapezoidal 
waveguide, resting on an 
plate, and two golden electrodes, placed above the core, which are embedded in a silica cladding. This design is based on [Loucas, 2023].
Cross-section of an electro-optical modulator. The electrodes are in yellow, the silica cladding in gray, and the LiNbO3 plate and waveguide in purple.
In this model, the electric field 
induced by the electrodes is simulated. Due to axial symmetry around the 
-axis and with the assumption that the EOM is sufficiently long, the electric field, has no dependence on 
. This allows to simplify the model and just model the cross-sectional area of the EOM.
As two different materials are being modeled, a region marker should be used to define the subregions.
The potential difference between the electrodes is of 
. With 
on the boundaries of the right-positive electrode, and 
on the left-negative electrode.
In this model, markers are used in the boundary condition specification. As these are DirichletCondition boundaries, point element markers are used.
The element marker for the left and right electrode are 3 and 4 respectively.
For nonlinear materials, the relationship between the electric field intensity 
and the electric flux density 
changes to:
where 
is the polarization when no electric field is present, known as the remanent polarization vector.
This equation takes into account the linear relationship 
as well as the additional polarization caused by Pr. In simple terms, the remanent polarization contributes to the overall electric flux density even in the absence of an external electric field.
The use of the constitutive relation expressed in Eq. 10 is limited to nonlinear materials that don't show hysteresis effects, such as ferroelectrics.
An axisymmetric simulation can be performed when a 3D solid has an axis of revolution. An axisymmetric model uses a truncated cylindrical coordinate system in 2D with independent variables 
instead of the cylindrical coordinates 
. The cylindrical coordinate variable 
disappears because the system is rotationally symmetric about the 
-axis.
The 2D axisymmetric electrostatic equation for linear materials is given as:
An axisymmetric model has the advantage that the computational cost in both time and memory is much less than in the case of solving a full 3D model.
The ElectrostaticPDEComponent function can produce the axisymmetric form of the electrostatic equation. To do so the parameter "RegionSymmetry" is set to "Axisymmetric".
A 2D axisymmetric model also applies when using a polarization density 
.
In this case, the electric potential is constant in the 
-direction, which implies that the electric field is tangential to the 
-plane. One can check the Spherical Capacitor model to see a 2D axisymmetric model.
In 1D axisymmetric geometries, the electric potential is also constant in the 
-direction, in addition to the 
-direction. The 1D axisymmetric electrostatic equation for linear materials is given as:
where 
is the thickness variable denoting thickness in the 
-direction.
A cylindrical capacitor is modeled in a 1D axisymmetric domain. In a cylindrical capacitor, the electric potential will vary only in the 
-direction. This capacitor consists of a solid cylindrical conductor of radius 
and length 
surrounded by a coaxial cylindrical shell of radius 
. The solid conductor will have an electric potential of 
and the shell will be grounded.
The space between the solid and the shell will have a relative permittivity of 1. The length of the capacitor is in the 
-direction, so the length of the device is equivalent to specify a thickness 
.
To verify the model, we will compute the capacitance and compare it with the theoretical value.
When doing volume integrals in axisymmetric geometries, one needs to include the term: 
.
Electrostatics sources are charge densities that effect the electric field. Electrostatics sources are categorized based on their charge distributions. We differentiate between volume, line, and point sources. With the special case that point sources in 2D models are equivalent to line charges in 3D models.
The source term 
on the right hand side of the electrostatic equation is used as a source to generate an internal electric flux density field 
, acting as a source when 
, or as a sink when 
.
where 
denotes the value of the volume charge density.
It is important that the mesh conforms to the geometrical bounds of the source term 
. The best way to do this is by explicitly generating the mesh for them. The finite element method, which is used as the solution method, benefits from a mesh where the individual elements do not cross a material boundary. An alternative is to make use of the MeshRefinementFunction. In this case, elements elements will cross the material boundary, but they will be small, if the mesh is sufficiently well refined. Combining both approaches is also possible.
A volume charge density 
or volumetric source can be used to model an arbitrarily shaped source or sink within the domain. The volume charge density value is always specified in units of 
. The name comes from the 3D incarnation of the volume charge but is used in other dimensions as well.
Consider the 1D case, where we have 
. When you specify the volume charge density 
in units of 
that value will be multiplied by the cross-sectional area 
in units of 
and we end up with a unit of 
.
Similarly in a 2D domain, where we have 
the volume charge density 
in units of 
is multiplied by the thickness 
in units of 
and we end up with units of 
.
In the following 2D example an electric dipole is modeled with two circular volume sources 
with opposite values 
.
In this example, the function RegionMember is used to specify the source region 
. The alternative option is to specify the sources with element markers, which will result in a more accurate solution.
for both sources with a value of 
and 
with a RegionMember function respectively:See this note about the set up of efficient PDE coefficients.
In this plot one can see the equipotential surfaces of a dipole and its respective electric field 
, which is normal to the equipotential lines everywhere.
A point charge 
or a point source can be used to model an internal source, when 
, or sink, when 
. A point charge is considered to have no spatial extension. A point charge can be used in 3D models or in 2D axisymmetric models.
In a 3D model a point charge can be positioned at any place.
In 2D axisymmetric models, a point charge can only exist if it is specified on the axis of symmetry. A point off, of the axis of symmetry would imply a line charge once the rotational aspect of the axisymmetric model is taken into consideration.
Point changes are not used in 2D models, as that would imply an equivalent to an out-of-plane line charge. See section Line Charges for details on this.
A point charge, shown in red, specified on the axis of symmetry of the 2D axisymmetric domain, in gray. One can see that revolving the axisymmetric domain does not transform the point charge into a line charge.
The units of the point change 
are always in units of 
while the volume charge density is in units of 
. Generally speaking we have that
where the 
is a Dirac delta function. The 
have a unit of 
in each of the directions of the source location 
This leads us to an expression for the change density 
that we can use:
The Dirac delta function, however, poses a problem in numerical simulations as it can not be resolved in the discretized spatial domain. This is because the Dirac delta function is singular at the source location 
. A second problem is that the evaluation of coefficients always happens within mesh elements, never at the edges. Hence, an approximation to the Dirac delta function is needed. The process of approximating the Dirac delta function is called regularization.
There are various regularized delta functions 
available [Peskin, 2002][Bilbao & Hamilton, 2017]. In this tutorial we choose:
where 
is the regularization parameter that controls the support of the regularized delta functions 
. Typically 
should have a size comparable to the mesh spacing 
. 
represents the difference 
.
The concept is best illustrated with examples. Before doing that, we create a helper function to compute a regularized Dirac delta function.
In the following 3D example a point charge 
is situated at 
with a value of 
.
To utilize the regularized delta function 
, we choose the regularization parameter to be 9 times the value of the mesh spacing: 
.
The regularization parameter 
must be sufficiently large so an integration point falls into the source point region. The following visualization shows a ball with a radius of 
in the point source location within the mesh and its integration points.
The figure shows a point source represented by a red ball with radius 
located at 
. The blue points represent the integration points of each mesh element. The integration points are the coordinates at which the PDE coefficients are evaluated. The point source needs enough extend 
such that sufficient integration points are inside it such that the numerical integration works.
The concept for axisymmetric models is the same as before and the expression for the charge density volume is given as:
Here the Dirac delta functions provide the units of 
at the source location 
in the axis of symmetry.
A line charge or a layer source 
models a source or a sink that is too thin to have a thickness in the model geometry. A line charge can be made use of in 3D, 2D, and 2D axisymmetric domains.
The units of 
are always specified in units of 
.
In 3D domains, line charges are can be specified along edges of a geometry, including lines floating arbitrary in the geometry.
In 2D axisymmetric domains, two options are available: a line charge specified on the symmetry axis or specified as a point that get's its length through the rotation of the simulation domain around the axis of symmetry.
In 2D, a point also represents a line charge because of the thickness 
of the model in the out-of-plane direction:
On the left: a 3D line charge and a plane. To the right: The same line charge from a 2D perspective.
In 3D models, a line charge can be specified at any edge of a geometry:
and the variations 
and 
are possible.
Here 
is a Dirac delta function at the source location 
and provides the units of 
and 
is the value of the line charge.
To model a line charge out-of-plane it is necessary to specify a point as a source.
Since the point has no spatial extension in all directions, the Dirac delta function 
should be applied on each dimension (i.e. 
) of the modeling domain 
:
Here 
is a Dirac delta function at the source location 
in units of 
, 
is the value of the line charge, and 
is the thickness in units of 
.
In the following 2D example, a two-dimensional dipole consisting of two lines with charge density 
on the 
-axis at positions 
, respectively. The line charges have a value of 
.
Due to symmetry only a quarter of the geometry is considered. At the left boundary, a ground potential is applied. A ground potential, 
,is equivalent to specify an antisymmetry boundary condition. At that boundary, the tangential component of the electric field 
is zero. At the top and right boundary, also a ground potential will be applied. And, at the bottom boundary, a symmetry boundary condition will be applied.
To utilize the regularized delta function 
, we choose the regularization parameter to be one-quarter of the mesh spacing: 
.
For this problem, an analytical solution exists for the electric potential 
given by:
where 
is the distance between the charge densities 
.
In 2D axisymmetric models, a line charge can be specified on the axis of symmetry. The volume charge density 
is given by:
Here 
is a Dirac delta function in units of 
at the source location 
and 
is the value of the line charge.
Secondary quantities
Once the simulation is complete, one can derive secondary quantities that provide deeper understanding of the model. These secondary quantities include: capacitance, forces, torques and energy-related parameters. These quantities provide deeper insights into the behaviour of electrostatic systems and guide the design and analysis processes for a wide range of applications.
Electrostatic energy
The Electrostatic energy 
is the energy that is stored in the electric field of charged particles and it can be expressed by the following equation:
where the electrostatic energy density 
is
In the case for a linear dielectric material with 
the electrostatic energy will be given by:
In the above example of the parallel capacitor the fields EfieldC and DfieldC were computed. In this section, the energy of the parallel capacitor is computed.
Capacitance
The capacitance 
is the property of a device to store electric charge and it is expressed as the ratio of the magnitude of the charge to the potential difference seen in the device. In a capacitor the capacitance is given by:
where 
is the charge's magnitude held in one of the plates of the capacitor, and 
is the potential difference between the plates.
In order to calculate the capacitance of the parallel-plate capacitor modeled in the multiple materials section, one first needs to know the charge distributed on one of the plates. The potential difference is already a know variable.
The charge 
will be distributed as a surface charge density 
. So, to obtain 
from 
we need to do a surface integral:
Note, that 
does not need to coincide with the boundary of the region 
.
From Gauss's law in its integral form, one can deduce that the surface charge density equals the normal component of the displacement field: 
. Therefore, the integral to solve is the following:
The field 
leaves or enters the metal plates in the direction of the normal, which in this case is the 
-direction. To completely capture the field 
leaving the metal plate, we integrate 
over a line from 
to 
.
In the above example of the parallel capacitor the fields EfieldC and DfieldC were computed. In this section, the energy of the parallel capacitor is computed.
The computed charge is the charge held in the lower plate.
In other cases, for example, when one specifies a current flowing into the capacitor, the way to know the voltage 
at the electrode is by using FEMNBoundaryIntegrate or NIntegrate. The voltage is integrated over the boundaries of the electrode and then divided by the area of these.
In the above case we could have made use of a symbolic region representation. In real world applications a symbolic representation may not always be available and even if available if may differ from the numeric representation. This is because the mesh is only an approximation to a symbolic representation of a geometry. For this reason it may in general be beneficial to directly use the numerical approximation of the region in question. This can be done directly using FEMNBoundaryIntegrate.
When using FEMNBoundaryIntegrate one needs to specify the element marker of the desired boundary surface. In this case , the number 3 refers to the upper plate boundary. For more information on markers in a mesh see the section Markers in the element mesh generation tutorial.
The theoretical values for a parallel capacitor is given by:
where 
is the plate area and 
is the distance between the plates.
The first thing to notice is that the theoretical values is lower than the numerical one. This difference is because the theoretical equation does not consider the fringing fields surrounding the device, which in the model case we are. In addition, the constraints imposed on the boundaries affect the capacitance value.
Ideally,when the surrounding space of the device is considered, in order to obtain the most accurate capacitance value, an infinite space should be modeled. To solve this there exists two approaches: extend the domain until the capacitance values converges to a specific value or to use infinite domains.
An alternative way of calculating the capacitance value is to use the electrostatic energy equation of a capacitor, given by:
where 
is the capacitance and 
is the potential difference. Solving for 
one gets:
The energy method offers a better approach to calculate the capacitance values, because one does not need to know in which boundary to do the integral and specify it.
Electrostatic force
This section shows how to compute electrostatic forces between objects. Computing forces can be necessary for several reason in various engineering and scientific applications. For example, many devices rely on electrostatic forces for control and actuation. Sensors also rely on changes in electrostatic forces, so they require accurate modelling to ensure precision in their designs.
Two methods are shown here. The first method is specific to compute only forces between metallic surfaces using an equation derived from 
. This approach is a bit simpler then the second case. The second method is a more general one, which can be used to compute forces between dielectrics or metals. This second method uses the Maxwell stress tensor 
.
Just as there is a Coulomb force between charges, there is also a force between continuous charge distributions. For example, for a surface charge density, the force per unit area is:
By doing a line integral over the boundary one gets the total electrostatic force 
.
In this section, the 
-component of the electrostatic force between the plates of the parallel capacitor is computed. For this computation the normal electric flux density field will be the 
-component 
.
Alternatively, when no symbolic representation of the line in question is available one can use FEMNBoundaryIntegrate.
If the geometry of the model is more complicated one can calculate the normal component to any surface by using BoundaryUnitNormal.
For force calculations it is not recommended to use FEMNBoundaryIntegrate in combination with BoundaryUnitNormal, because computation and the orientation of the normals at internal boundaries is not well defined.
The numerical electrostatic force can be compare with the following estimate:
where 
is the area of the plate. For a parallel capacitor the equation reduces to:
where 
is given by 
and where 
is the distance between the plates. The relative permittivity of the capacitor has the value 
.
An alternative method to compute electrostatic forces is by employing the Maxwell stress tensor, which is a more general method that can be applied to dielectrics and also to magnetostatic forces.
With the Maxwell stress tensor, the force can generally be expressed as:
where 
is the Maxwell stress tensor, 
is the normal to a surface, 
is the domain boundary, 
is the domain volume, and 
is the Poynting vector.
In the electrostatic case, we have no time dependency and the equation simplifies to:
where 
is the unit outward normal vector from an object and 
the surface.
The Maxwell stress tensor, which has two indices, is defined as:
where 
is the magnitude of E squared and 
is the Kronecker delta, which is 1 if the indices are the same and zero otherwise.
The first index indicates the direction of the force and the second index indicates the direction of the surface.
is a force per unit area acting on a surface. To compute the stress for a given surface, one takes the dot product with the desired surface normal such that:
Then, 
is integrated over the surface of the object that the force acts on:
An equivalent formula [Ida & Bastos, 2013] is:
The electrostatic force between the plates is also computed here. As the force is in the 
-direction we only need 
.
In the next example we are going to compute the force exerted on a dielectric and metallic rod. We will use both methods and compare both of them. These rods are in a vacuum and exposed to a transverse parallel field. As a result of Newton's third law the forces must be equal and opposite.
Cross section of a metallic and dielectric rod in a parallel 
field.
The use of regions markers is employed again in this example.
As in this section forces are calculated at boundaries, a special mesh is created here to get the best accuracy in the forces computation. A refinement is wanted at the internal boundaries of the domain, for that a manually boundary mesh is generated. Furthermore, if a NumericalRegion is used, a high-quality second-order mesh with minimal elements is generated. To use NumericalRegion presents other advantages when integrating the forces expressions, it allows us to use NIntegrate with a better accuracy.
To manually generate the boundary mesh, the idea is to generate each part of the region separately an then combine them in a new boundary mesh.
Note that the metal rod is removed from the boundary mesh.
Details on this matter can be seen in the Element Mesh Generation tutorial in the Numerical Regions. section.
Boundary markers were assigned to the boundaries of the boundary mesh.
The markers 1,2,3 were assigned to the boundaries of the metal rod, the dielectric rod, and of the box respectively. The point markers follow the same ID numbers.
A ground potential is set at the metal rod.
The parallel field is set up using an ElectricFluxDensityValue. For more details see the Electric Flux Boundary condition section.
When forces are computed at an interface between a dielectric and another material, one should pay special attention to the computed fields. At the interface, the electric field intensity 
and the electric flux density field 
can show discontinuities.
To address this discontinuity we make use of DiscontinuousInterpolatingFunction or EvaluateOnElementMesh. These functions will prioritize the value from one of the subregions. The default uses the order given in "MeshElementMarkerUnion". To specify a priority you can give an option "MarkerPriority".
In this case, we will use both the default priority, which prioritizes the dielectric region, and also the "MarkerPriority"
which will prioritize the vacuum region.
:First, we will compute the 
-component of the force on the metal rod. As we are dealing with a metallic surface, one can use Eq. 12:
At the interface between the metal rod and the vacuum there are no discontinuities because the metal region is not considered in the domain, so it won't matter which field is used.
Alternatively, one can use the Maxwell stress tensor:
Now, since we would like to compute the force applied to the dielectric the Maxwell stress tensor method needs to be used.
To correctly compute the 
-component of the force on the dielectric at the dielectric-vacuum interface one needs to prioritize the vacuum values over the dielectric values.
The value is of opposite sign and is really close to the previous calculated value, as one should expect.
There is another formula [Küpfmüller, 1968] for calculating forces in dielectrics, which is the following one:
where the subscript 
refers to the dielectric and the 
to the vacuum. The 
refers to the tangential component of the field and 
to the normal component.
The advantage of using the maxwell stress tensor is that it can be applied to multiple dielectric materials or to one dielectric or to a non-dielectric.
As only the 
-component of the normal force is calculated, the expression just need to be multiplied by the direction cosine of the angle between the normal and the 
-axis 
, instead of the complete unit normal 
.
Boundary conditions in electrostatics
Boundary conditions for electrostatics application fall into one of two categories. One are Dirichlet conditions which specifies the potential on conductors. An example is a capacitor whose electrodes are fixed to a potential, typically one is fixed to ground potential 
and the other to a specific value different from zero.
On the other hand, we have boundaries that specify the value of the normal component of the displacement field 
. These boundaries are realized with Neumann value boundary conditions.
In addition to these two type of boundary conditions, there is also a third type known as periodic boundary conditions, which are useful when the geometry has repetitive sections, for example, a multipolar rotating electric machine, and one wants to reduce the extension of the domain. So, this type of boundary condition specifies the electric potential 
at one part of the boundary to be the same at another part.
For these three general boundary condition types, the following electrostatic specific boundary conditions are introduced:
At an internal boundary separating two different dielectric materials the following conditions must be satisfied:
where 
is a surface charge density placed deliberately at the interface.
Eq. 13 and Eq. 14 can be better explain if we denote a new nomenclature, where 
refers to the normal component and 
to the tangential component.
Eq. 15 can be written as follows:
which means that the normal component of the electric flux density will be discontinuous if a surface charge density is placed deliberately at the boundary. If not, then the electric flux will be continuous across the boundary.
Nevertheless, the electric field will not be continuous. If we substitute 
by 
, then
showing that the normal component of 
is discontinuous at the boundary.
With Eq. 16, we deduce the following expression:
which means that the tangential component of the electric field 
is continuous along the internal boundary. On the other hand, if we substitute 
by 
, then we deduce the following equation:
which means that the tangential component of the electric flux is discontinuous.
Electric potential boundary condition
The purpose of an electric potential boundary condition is to set a specific potential on some part of the boundary.
The electric potential boundary condition is generally applied at parts of boundaries that model metal contacts such as electrodes. In most examples one electrode will be fixed at a value 
while a second electrode be grounded at 
.
An electric potential boundary condition set to zero 
can be used as an antisymmetry boundary condition.
An antisymmetry boundary condition models a boundary where a function 
being considered exhibits antisymmetry with respect to an axis. This means, that the value of the function being studied on one side is equal to its negative on the other side of the axis.
In electrostatics, the electric fields can exhibit this antisymmetric behaviour: Take for example the case of an electric field from a dipole source. Assume we have a positive and a negative charge separated by a distance 
and aligned along the 
-axis. One can see from the figure that the electric field exhibits an antisymmetric behaviour with respect to the dashed red line indicating the axis of antisymmetry.
The figure shows an electric field 
from a dipole, where the blue and red dots are the positive and negative charges respectively. An antisymmetry exists between the right and left side of the line of antisymmetry. We have 
.
It is important to realize that the tangential component of the electric field will be zero at the line of antisymmetry and that the field will only have a normal component. Thus, the antisymmetry boundary condition is given by:
and is fulfilled if the potential 
is set to zero at the boundary 
, due to the fact that the gradient of a constant potential is zero:
Just like a symmetry boundary condition models a boundary with mirror symmetry along an axis and is used to reduce the model size, an antisymmetry boundary condition can also be used to reduce the simulation domain size, as it is shown the following picture.
Reduced simulation domain of an electric dipole, where just the positive charge is modeled. A ground potential is set at the boundary 
to force the antisymmetry behaviour 
.
Electric flux density boundary condition
The purpose of an electric flux density boundary condition is to define the normal component of the electric flux density field at a boundary. In other words, it is used to give a flux condition on the specified boundary.
The electric flux density boundary condition can be specified in two different ways: One is by specifying the surface charge density 
and the second option is to specify the boundary electric flux density 
.
With a surface charge density 
the boundary condition is given by:
where the definition of 
is given by 
and by inserting 
we get to 
. The minus sign is added to Eq. 17 in order to be able to insert the field 
into the equation.
This Neumann value is directly related to the divergence of 
seen in the electrostatic equation.
A positive value denotes the inward electric flux, and a negative value denote an outward flux. In other words the flux is specified opposite to the outward normal 
.
A surface charge density can be specified at an outer boundary or at an interior boundary between two nonconducting materials. In this second case, the boundary condition models:
With a boundary electric flux density 
and boundary unit normal 
, the boundary condition is given by:
If the electric flux density is specified as a vector of zeros value or the surface charge density specified with a zero value, then there is no surface charge density at that boundary. This is equivalent to a Neumann zero boundary condition. The zero electric flux is the implicit default boundary condition used if no other boundary condition is specified at a given boundary:
At interior boundaries, the zero electric flux value means that no electric flux can penetrate the boundary and that the electric potential is discontinuous across the boundary.
In the following example, a capacitor is charged by specifying a surface charge density on an upper plate, and a ground potential on a lower plate.
In the following example a glass rod of circular cross-section is exposed to an external field. A normal component of an electric flux density is applied on the right and left boundary of a box which surrounds the glass rod.
Cross-section of a glass rod exposed to an external 
field.
A point at 
was added to the mesh because a ground potential will be applied there. Also, region markers were specified in the mesh to assign different relative permittivities.
The glass rod will have a relative permittivity of 7, while the surroundings of 1.
A ground potential will be set in the point 
as a reference value.
Symmetry boundary condition
A symmetry boundary condition is used when the computational domain, and the expected electric potential, have mirror symmetry along an axis of the simulation domain.
The symmetry boundary condition is given by:
which states that the normal component of the electric field intensity is zero.
If on some part of the boundary no boundary condition is set an implicit Neumann zero boundary condition is used. A symmetry boundary condition is equivalent to a zero charge density boundary condition.
In the following example the right half of a dielectric region of a capacitor is modeled.
Periodic boundary condition
The purpose of a periodic boundary condition is to map the electric potential from one boundary to another in order to model periodicity of the domain.
Given a function 
that maps the electric potential 
from the periodic boundary 
to the targeted boundary 
, the periodic boundary condition can be written as:
As an example an array of circular electrodes are modeled in a domain of perfect vacuum. It is possible to perform the simulation with just modeling one pair of electrodes by using the periodic boundary condition. The array is of 
equally spaced electrodes. The upper electrodes have an electric potential of 
while the lower electrodes are grounded.
Array of circular electrodes in a vacuum. The boundaries of the upper electrodes, in red, have a potential of 
, while the lower boundaries, in blue, have a potential of 0.
In the reduced domain only a single pair of electrodes is modeled. At the left boundary of the box a periodic boundary conditions is applied to map the the boundary electric potential to the target boundary, which is at the right boundary of the reduced box.
The reduced domain and the locations of the periodic boundary conditions.
Remember that electrodes inside a domain are removed, and just its boundaries are used.
The electrodes are modeled by removing them from the domain and replacing them with a electric potential boundary condition.
Note that the periodic condition cannot be in any boundary where a Dirichlet condition has been applied.
The same procedure for calculating the capacitance applies here, except that the final result must be multiplied by the total number of sections. We will use the energy method, so we will use the following equations:
The capacitance is lower than what would be obtained if the analysis was done in the full domain, this is because the left and right ends of the electrode array were not taken into account.
Nomenclature
References
1. Bilbao, S., & Hamilton, B. (2017, October). Directional source modeling in wave-based room acoustics simulation. In 2017 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA) (pp. 121-125). IEEE.
2. Cardoso, J. R. (2018). Electromagnetics through the finite element method. Boca Raton: CRC Press.
3. Ida, N., & Bastos, J. P. (2013). Electromagnetics and calculation of fields. Springer Science & Business Media.
4. Küpfmüller, K. (1968). Einführung in die theoretische Elektrotechnik. Springer.
5. Loucas, A. (2023). Optimization of lithium niobate electro-optical modulators (Master's thesis). Stevens Institute of Technology, New Jersey, NJ.
6. Sadiku, M. N. O. (2011). Elements of electromagnetics (5th ed.). Oxford University Press.
7. Peskin, C. S. (2002). The immersed boundary method. Acta numerica, 11, 479-517.