Intro to Computational Electrostatics: Method of Moments
Table of Contents
1 Method of Moments for Electrostatics
The below derivation and examples are originally from Electromagnetic Fields and Waves, an excellent book by Magdy F. Iskander. I have reworked and summarized them here.
Oftentimes you want to answer questions about an electrical conductor at rest (i.e., not under some dynamic perturbation). In other words, you want to characterize it.
Classically speaking, this involves answering a few simple questions: we want to know the electrical potential throughout the conductor, the electric field, and the capacitance. (From there we can start asking and answering questions about the group velocity of a wave passing through the conductor, current v. temperature, and so on.)
First, we need some equations.
1.1 Basic Equations
1.1.1 Faraday's Law
First, let's take Faraday's Law:
\begin{equation} \oint_c \boldsymbol{E} \cdot d\boldsymbol{l} = - \frac{d}{dt} \int_s \boldsymbol{B} \cdot d\boldsymbol{s} \end{equation}This says that an electric field through a closed shape (contour) is coupled to a time-varying magnetic field through the surface formed by that same closed shape, in a direction such that the magnetic field opposes the change (Lenz's law).
Now, a static electric field corresponds to no time-varying magnetic field (a static magnetic field may still be present, but decoupled):
\begin{equation} \oint_c \boldsymbol{E} \cdot d\boldsymbol{l} = 0 \end{equation}1.1.2 Boundary Conditions (possibly not needed)
We also need the equation that determines the boundary conditions across two materials:
\begin{equation} \boldsymbol{n} \cdot (\epsilon_0 \epsilon_{r_1} \boldsymbol{E}_1 \cdot \epsilon_0 \epsilon_{r_2} \boldsymbol{E}_2) = \rho_s \end{equation}1.1.3 Coulomb's Law
The main focus of this section will be Coulomb's law:
\begin{equation} \Phi = \frac{Q}{4\pi\epsilon_0(\epsilon_r)\boldsymbol{r}} \end{equation}Coulomb's law expresses the electric potential at a distance \(\boldsymbol{r}\) from a point charge whose value is \(Q\).
This law is additive, so, for example, you can sum up the charge of a bunch of point charges distributed in space:
\begin{equation} \Phi_{dist} = \sum^N_{k=1} \frac{Q_k}{4\pi\epsilon_0(\epsilon_r)\boldsymbol{r}_k} \end{equation}As you might expect, you can take another step and imagine a surface made up of infinitesimal point charges also:
\begin{align} \Phi_{surf} &= \int_s \frac{\rho_s d\boldsymbol{s}}{4\pi\epsilon_0(\epsilon_r)\boldsymbol{r}} &= \frac{1}{4\pi \epsilon_0\epsilon_r} \int_s \frac{\rho_s d\boldsymbol{s}'}{|\boldsymbol{r} - \boldsymbol{r}'|} \end{align}In the surface calculation, you have to include the effect of all of the charges present in the whole surface, for any given distance \(\boldsymbol{r}\). This formulation is a concise one that weights each of the charge terms in the integral by the distance from that point to the point being evaluated.
1.2 Motivating Example - Coaxial Conductor
Before we get into the computational angle, let's compute a motivating example.
Assume we have a coaxial conductor of length \(L\) and outer radius \(a\). The surface of this conductor is held at a constant voltage \(V\), while an inner grounding conductor of radius \(\delta\) is held at \(0\) volts. They are separated by a dielectric of relative permittivity \(\epsilon_r\).
We want to find the charge distribution on the surface, \(\rho_s\).
First, we track down which of Maxwell's equations mentions \(\rho_s\): Gauss' Law.
\begin{align} \oint_s \epsilon_0\epsilon_r \boldsymbol{E} \cdot d\boldsymbol{s} = \int_v \rho_s dv \end{align}The first thing we need to do is define the surface and volume elements in terms of the topology of the cylindrical conductor. Since we know that the cylinder has constant radius about its length, we can write double integrals for both sides:
\begin{align} \int^{2\pi}_{\phi = 0} \int^{L}_{ z = 0} \epsilon_0\epsilon_r \boldsymbol{E} \cdot d\boldsymbol{s} = \int^{2\pi}_0 \int^{L}_{0} \rho_s dv \end{align}We also know that the electric field only extends in the radial direction. So we just need to write in the differential surface element on the left side projected onto the radial vector, and the differential volume element on the right.
\begin{align} \int^{2\pi}_{\phi = 0} \int^{L}_{ z = 0} \epsilon_0\epsilon_r \boldsymbol{E}_\rho \boldsymbol{a}_\rho \cdot \rho d\phi dz \boldsymbol{a}_\rho = \int^{2\pi}_0 \int^{L}_{0} \rho_s a d\phi dz \end{align}Nearly everything in either equation is a constant, so we can quickly integrate.
\begin{align} \rho 2\pi L \epsilon_0\epsilon_r \boldsymbol{E}_\rho \boldsymbol{a}_\rho &= % 2\pi L \rho_s a \\ \boldsymbol{E}_\rho \boldsymbol{a}_\rho &= \frac{\rho_s a}{\rho \epsilon_0\epsilon_r} \end{align}Now we just need to relate the potential and the electric field.
\begin{align} V = - \int^{a}_{\delta} \boldsymbol{E}_{\rho} \cdot d\boldsymbol{l} \end{align}Again, we know that the potential is measured in the radial direction (inner and outer conductor), so we can set \(d\boldsymbol{l} = d\rho \boldsymbol{a}_\rho\).
\begin{align} V &= - \int^{a}_{\delta} \frac{\rho_s a}{\rho \epsilon_0\epsilon_r} \cdot d\rho \boldsymbol{a}_\rho \\ &= - \left[ \frac{\rho_s a}{\epsilon_0\epsilon_r} \right] \left[ \ln |a| - \ln |\delta| \right] \\ &= \frac{\rho_s a}{\epsilon_0\epsilon_r} \ln \left| \frac{\delta}{a} \right| \end{align}Assuming that we know the voltage, then the charge distribution is easily found as
\begin{align} \rho_s = \frac{V \epsilon_0\epsilon_r}{a \ln \left| \frac{\delta}{a} \right|} \end{align}Since the conductor is cylindrical, you can also compute the charge per unit length by nothing that \(\rho_l = 2 \pi a \rho_s\).
From here you can compute capacitance as needed, as well.
1.3 Deriving the Method of Moments
Those methods are fine if you want to solve the equations directly, but if you don't have an analytically tractable charge distribution, or in fact don't know the charge distribution at all, only its measured potentials, then you need another method.
1.3.1 Derivation
Imagine we have a cylindrical conductor of length \(L\), radius \(a\), held at a constant voltage of \(1\) V. Its potential according to Coulomb's law would be
\begin{equation} \Phi(\boldsymbol{r}) = \frac{1}{4 \pi \epsilon_0\epsilon_r} \int_s \frac{\rho_s(\boldsymbol{r}') d\boldsymbol{s}'}{|\boldsymbol{r} - \boldsymbol{r}'|} \end{equation}Instead of working out the integral for that topology directly, let's make the assumption that if we divide the conductor into small enough sections of length \(\Delta_{l_i}\), the charge distribution will be essentially constant. Specifically:
\begin{align} \Phi(\boldsymbol{r}) &\approx \frac{1}{4 \pi \epsilon_0 \epsilon_r} \sum^N_{i=1} \frac{\rho_{s_i} \Delta s_i}{|\boldsymbol{r} - \boldsymbol{r}'|} \\ \Delta s_i &= 2\pi a \Delta{l_i} \\ \Delta{l_i} &= \frac{L}{N} \end{align}This leads to \(N\) equations computed at the "center" of each section (since it doesn't matter what point you choose to compute the potential at in each constant section):
\begin{align} 1 &= \frac{1}{4\pi\epsilon_0\epsilon_r} \left[ \frac{\rho_{s_1} \Delta{l_1}}{|x_1 - x_1|} + \frac{\rho_{s_2} \Delta_{l_2}}{|x_1 - x_2|} + \frac{\rho_{s_3} \Delta_{l_3}}{|x_1 - x_3|} + \cdots + \frac{\rho_{s_N} \Delta_{l_N}}{|x_1 - x_N|} \right] \\ 1 &= \frac{1}{4\pi\epsilon_0\epsilon_r} \left[ \frac{\rho_{s_1} \Delta{l_1}}{|x_2 - x_1|} + \frac{\rho_{s_2} \Delta_{l_2}}{|x_2 - x_2|} + \frac{\rho_{s_3} \Delta_{l_3}}{|x_2 - x_3|} + \cdots + \frac{\rho_{s_N} \Delta_{l_N}}{|x_1 - x_N|} \right] \\ &\vdots \\ 1 &= \frac{1}{4\pi\epsilon_0\epsilon_r} \left[ \frac{\rho_{s_1} \Delta{l_1}}{|x_N - x_1|} + \frac{\rho_{s_2} \Delta_{l_2}}{|x_N - x_2|} + \frac{\rho_{s_3} \Delta_{l_3}}{|x_N - x_3|} + \cdots + \frac{\rho_{s_N} \Delta_{l_N}}{|x_N - x_N|} \right] \end{align}Note that the equations are singular because each of the self-compared terms have zero in their denominator. One way to get around this is to note that the charge is distributed on the conductor surface, and do \(N\) integral equations (for which generally we can choose a tractable approximate geometry) on the cylindrical section:
\begin{equation} \Phi = 1 = \frac{1}{4 \pi \epsilon_0 \epsilon_r} \int^{\frac{\Delta{l_i}}{2}}_{-\frac{\Delta{l_i}}{2}} \int^{2\pi}_0 \frac{\rho_{s_i} a d\phi dx'}{\sqrt{a^2 + x'^2}} \end{equation}To solve this equation, we need to use an Euler substitution.
- Euler Substitution
For any equation of the form
\begin{equation} \int \frac{1}{\sqrt{c^2 + x^2}} dx \end{equation}the Euler substitution is
\begin{align} \sqrt{c^2 + x^2} &= -x + t \\ x &= \frac{t^2 - c^2}{2t} \\ dx &= \frac{t^2 + c^2}{2 t^2} \\ \sqrt{c^2 + x^2} &= \frac{t^2 + c^2}{2t} \end{align}This integral is then
\begin{align} \int \frac{1}{ \sqrt{c^2 + x^2} } dx &\implies \int \frac{ \frac{t^2 + c}{2t^2} }{ \frac{t^2 + c}{2t} } \\ &= \int \frac{dt}{t} = \ln |t| + C = \ln \left| x + \sqrt{x^2 + c^2} \right| + C \end{align} - Simplifying the result
Thus, the cylindrical integral can be solved to get
\begin{equation} \frac{\rho_{s_i} a}{2 \epsilon_0} \ln \frac{ \frac{\Delta{l_i}}{2} + \sqrt{a^2 + \left( \Delta{l_i}{2} \right)^2 }}{ -\frac{\Delta{l_i}}{2} + \sqrt{a^2 + \left( \Delta{l_i}{2} \right)} } \end{equation}We can simplify this further if we assume that the radius \(a << \Delta l_i\). In the numerator, let's just throw out \(a\). In the denominator, we need to compute the Taylor expansion of the square root and then throw out higher order terms.
- Taylor expansion of square root\begin{equation} \sum^{\infty}{n=0} \frac{(-1)^n (2n)! d^n}{(1-2n)n!^2 4^n N^{2n -1}} = N + \frac{d}{2N} - \frac{d^2}{8N^3} + \frac{d^3}{16N^5} + \cdots \end{equation}
- Constructing the linear equations
With that in mind, the equation simplifies to:
\begin{align} 1 &= \frac{\rho_{si}a}{2\epsilon_0} \ln \frac{\Delta l_i}{-\frac{\Delta l_i}{2} + \frac{\Delta l_i}{2} + \left(\frac{1}{2}\right) \left(\frac{a^2}{\frac{\Delta l_i}{2}}\right)} \\ &= \frac{\rho_{si}}{2\epsilon_0} \ln \frac{\Delta l_i^2}{a} \end{align}Now we can build the equations:
\begin{equation} \begin{bmatrix} 4 \pi \epsilon_0 \epsilon_r \\ 4 \pi \epsilon_0 \epsilon_r \\ 4 \pi \epsilon_0 \epsilon_r \\ \vdots \\ 4 \pi \epsilon_0 \epsilon_r \end{bmatrix} = \begin{bmatrix} 4 \pi a \ln \frac{\Delta l_1}{a} & \frac{2\pi a \Delta l_2}{|x_1 - x_2|} & \frac{2\pi a \Delta l_3}{|x_1 - x_3|} & \cdots & \frac{2\pi a \Delta l_N}{|x_1 - x_N|} \\ \frac{2\pi a \Delta l_1}{|x_2 - x_1|} & 4 \pi a \ln \frac{\Delta l_2}{a} & \frac{2\pi a \Delta l_3}{|x_2 - x_3|} & \cdots & \frac{2\pi a \Delta l_N}{|x_2 - x_N|} \\ \frac{2\pi a \Delta l_1}{|x_3 - x_1|} & \frac{2\pi a \Delta l_2}{|x_3 - x_2|} & 4 \pi a \ln \frac{\Delta l_3}{a} & \cdots & \frac{2\pi a \Delta l_N}{|x_3 - x_N|} \\ \frac{2\pi a \Delta l_1}{|x_N - x_1|} & \frac{2\pi a \Delta l_2}{|x_N - x_2|} & \frac{2\pi a \Delta l_N}{|x_N - x_3|} & \cdots & 4 \pi a \ln \frac{\Delta l_N}{a} \end{bmatrix} \begin{bmatrix} \rho_{s_1} \\ \rho_{s_2} \\ \rho_{s_3} \\ \vdots \\ \rho_{s_N} \end{bmatrix} \end{equation}
1.4 Algorithm
The algorithm is now not too difficult. Here are the steps for finding the static charge distribution across a surface:
- Separate the solid into differential (surface) sections (for a
higher-dimensional case, repeat the derivation for a different \(\Delta{s_i}\)
and adjust the distance calculation appropriately) and formulate the
off-diagonal terms directly.
- Note: The distance calculation can be memoized and computed once based on the spacing of the grid centers.
- Also note that the diagonal terms are based on an integral of the chosen tesselated approximate surface (in the above case, it's a cylindrical slice).
- Formulate the diagonal terms by integrating the volume element of the differential sections directly.
- Fill in the matrix as above and invert.