Intro to Computational Electrostatics: Finite Element Method
Table of Contents
1 Finite Element Method for Electrostatics
The below derivation and examples are originally from Numerical Techniques in Electromagnetics, by Matthew Sadiku. I have reworked and summarized them here.
This document introduces the finite element method for solving Laplace's equation. Since the theoretical background is the same as the finite difference method, refer to that document as needed.
The finite element method is a mechanical/civil engineering technique mathematically treated in 1943, but found use in electromagnetics in 1968. According to Sadiku, FD and MM are "conceptually simpler" and "easier to program" than FEM, but FEM is more powerful when it comes to complex geometries and inhomogeneous media. I agree with that.
So, the goal here is to deal with the conceptual complexity part.
1.1 Basic Algorithm
The algorithm is basically the following steps:
- Derive equations that characterize the quantity of interest in each of the elements.
- Discretize the region of interest into a finite number of elements.
- Assemble the elements into a complete system of equations.
- Solve the complete equations.
In this document, we'll consider Laplace's equation in two dimensions, and limit the discussion to triangular elements. However, it's worth noting that there are variations for each of these four steps that can be chosen independently of the others. I'll mention them as we go.
1.1.1 Characterizing The Potential
Let's assume that we've discretized the region of interest into \(N\) elements that are disjoint. In that case, we should write equations that are nonzero outside of their individual regions, and that are continuous across their boundaries with the adjacent elements.
An immediate consequence of that stipulation is that the approximate solution is a sum of \(N\) element equations:
\begin{equation} V(x, y) \approx \sum^N_{e=1} V_e(x,y) \end{equation}Next, let's assume the element equation takes the form of a polynomial:
\begin{equation} V_e(x,y) = a + bx + cy \end{equation}This linear equation implies that the electric field is uniform inside the element:
\begin{equation} \boldsymbol{E}_e = - \nabla V_e = - (b \boldsymbol{a}_x + c \boldsymbol{a}_y) \end{equation}You could use higher-order terms to describe the potential, or consider a different geometry instead of a triangle.
1.1.2 Discretizing the Region
Assuming that we've gone with the equations from the previous section, we now need to write them out for each node in a triangular element \(e\).
\begin{equation} \begin{bmatrix} V_{e1} \\ V_{e1} \\ V_{e3} \end{bmatrix} = \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} \end{equation}Note that these are just the barycentric coordinates of a triangle. Since we need to determine the coefficients as part of the solution method, let's expose them algebraically:
\begin{equation} \begin{bmatrix} a \\ b \\ c \end{bmatrix} = \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{bmatrix}^{-1} \begin{bmatrix} V_{e1} \\ V_{e1} \\ V_{e3} \end{bmatrix} \end{equation}We can write out the solution explicitly using Cramer's rule:
\begin{equation} V_e = \frac{1}{2A} \begin{bmatrix} 1 & x & y \end{bmatrix} \begin{bmatrix} (x_2 y_3 - x_3 y_2) & (x_3 y_1 - x_1 y_3) & (x_1 y_2 - x_2 y_1) \\ (y_2 - y_3) & (y_3 - y_1) & (y_1 - y_2) \\ (x_3 - x_2) & (x_1 - x_3) & (x_2 - x_1) \end{bmatrix} \begin{bmatrix} V_{e1} \\ V_{e1} \\ V_{e3} \end{bmatrix} \end{equation}Note that \(2A\) is the determinant of the coordinate matrix, equal to twice the area of the element:
\begin{equation} 2A = (x_1 y_2 - x_2 y_1) + (x_3 y_1 - x_1 y_3) + (x_2 y_3 - x_3 y_2) \end{equation}Here there's a convention to follow: shape nodes should be numbered counterclockwise to ensure that \(A\) is positive, if \(x\) increases to the right and \(y\) increases upward. (In other words, if using the image processing convention of the upper left corner being the origin, shape nodes should be numbered clockwise.)
Let's do one last thing before we move on. The equations need to be packaged up into something we can name succinctly in future steps. Since the geometry is constant, we can take everything to the left of \(\begin{bmatrix} V_{e1} \\ V_{e1} \\ V_{e3} \end{bmatrix}\) and package it up as a function. Writing out the product, we get:
\begin{equation} V_e = \sum^3_{i=1} \alpha_i (x,y) V_{ei} \end{equation}The \(\alpha\) functions here are just the previous matrices multiplied out:
\begin{align} \alpha_1(x,y) &= \frac{1}{2A}\left[(x_2 y_3 - x_3 y_2) + (y_2 - y_3)x + (x_3 - x_2)y \right] \\ \alpha_2(x,y) &= \frac{1}{2A}\left[(x_3 y_1 - x_1 y_3) + (y_3 - y_1)x + (x_1 - x_3)y \right] \\ \alpha_3(x,y) &= \frac{1}{2A}\left[(x_1 y_2 - x_2 y_1) + (y_1 - y_2)x + (x_2 - x_1)y \right] \end{align}These functions have the interesting property of being \(1\) at the node they correspond to, and \(0\) at the other nodes. In other words:
\begin{equation} \alpha_i(x_j, y_j) = \begin{cases} 1 & i = j \\ 0 & i \neq j \end{cases} \end{equation}1.1.3 Synthesizing the Equations
- The element coefficient matrix
Recall that we were trying to solve Laplace's equation before we got bogged down in all that algebra. Since we're looking for an electrostatic solution, we can note that there should not be any "extra" energy in the system.
\begin{equation} W_e = \frac{1}{2} \int \epsilon |\boldsymbol{E}_e|^2 dS = \frac{1}{2} \int \epsilon |\nabla V_e|^2 dS \end{equation}Unpacking that with our definition from the previous section, we get:
\begin{equation} W_e = \frac{1}{2} \sum^3_{i=1} \sum^3_{j=1} \epsilon_0 V_{ei} \left[ \epsilon_r^{(e)}\int \nabla \alpha_i \cdot \nabla \alpha_j dS \right] V_{ej} \end{equation}The integrand can be packed up into a matrix,
\begin{equation} C_{ij}^{(e)} = \epsilon_r^{(e)} \int \nabla \alpha_i \cdot \nabla \alpha_j dS \end{equation}which lets us write the functional as
\begin{equation} W_e = \frac{1}{2} \epsilon_0 \left[ V_e \right]^\top \left[ C^{(e)} \right] \left[ V_e \right] \end{equation}\(C^{(e)}\) is the "element coefficient matrix", and it describes the coupling between nodes \(i\) and \(j\), as we'll see in an example shortly.
At any rate, the coupling expressions are a function of the basic geometry, and their form is clearly determined in advance:
\begin{align} C_{11}^{(e)} &= \frac{1}{4A} \left[ (y_2 - y_3)^2 + (x_3 - x_2)^2 \right] \\ C_{12}^{(e)} &= \frac{1}{4A} \left[ (y_2 - y_3)(y_3 - y_1) + (x_3 - x_2)(x_1 - x_3) \right] \\ C_{13}^{(e)} &= \frac{1}{4A} \left[ (y_2 - y_3)(y_1 - y_2) + (x_3 - x_2)(x_2 - x_1) \right] \\ C_{21}^{(e)} &= C_{12}^{(e)} \\ C_{22}^{(e)} &= \frac{1}{4A} \left[ (y_3 - y_1)^2 + (x_1 - x_3)^2 \right] \\ C_{23}^{(e)} &= \frac{1}{4A} \left[ (y_3 - y_1)(y_1 - y_2) + (x_1 - x_3)(x_2 - x_1) \right] \\ C_{31}^{(e)} &= C_{13}^{(e)} \\ C_{32}^{(e)} &= C_{23}^{(e)} \\ C_{33}^{(e)} &= \frac{1}{4A} \left[ (y_1 - y_2)^2 + (x_2 - x_1)^2 \right] \end{align}
- Global element coefficient matrix
The global element coefficient matrix is formed by a slightly intricate procedure that follows some simple rules.
- Label all of the nodes with global indices.
- Define a square matrix \(C\) with as many rows as nodes.
- For each element of \(C_{ij}\), write an expression made up of a sum of all of
the local \(C_{ij}^{(e)}\), including all of the local nodes that are reached
along a path from global node \(i\) to to global node \(j\).
- If no direct path exists from global node \(i\) to \(j\), write 0.
The resulting matrix will have some interesting properties:
- It will be symmetric.
- It will grow sparse for large element count.
- It is singular.
- It also takes the same form regardless of the local element numbering.
The procedure for generating this matrix may seem confusing, but it may become more clear in the example.
1.1.4 Solving the equations
Since we stated we want to minimize the energy (static case), let's write out that statement's consequences: all of the derivatives of the energy function have to equal \(0\).
In general, this means \(\frac{\delta W}{\delta V_k} = 0 = \sum^N_{i=1} V_i C_{ik}\). This produces one equation for each \(V_i\).
There are at least two algorithms for solving this problem, which are introduced in the section on algorithms for calculation.
1.2 Algorithms for Calculation
1.2.1 Iterative Algorithm (Jacobi)
- Set unknown \(V_k\) to an initial value like \(0\) or \(\frac{1}{2}(V_{max} + V_{min})\).
- Calculate the potential \(V_k\) for all \(k\) according to that equation.
- Replace the potentials used in step 1 with the newly calculated potentials.
- Repeat the above two steps until \(|V_k \textrm{(previous)} - V_k \textrm{(current)}| < \delta\), where \(\delta\) is a convergence criterion.
1.2.2 Block Matrix Algorithm
Recall that the nodes whose value is set by the initial conditions are "fixed" or "prescribed" nodes, while the nodes whose value we want to solve for are "free" nodes. If we sort the terms corresponding to free nodes (with subscript \(f\)) and the terms corresponding to prescribed nodes (with subscript \(p\)) separately, we can write the energy functional in terms of block matrices:
\begin{equation} W = \frac{1}{2} \epsilon_0 \begin{bmatrix} V_f & V_p \end{bmatrix} \begin{bmatrix} C_{ff} & C_{fp} \\ C_{pf} & C_{pp} \end{bmatrix} \begin{bmatrix} V_f \\ V_p \end{bmatrix} \end{equation}Note that \(V_p\) is composed of constants. Thus, the derivative here is \(\frac{\delta W}{\delta V_f}\) –we can ignore the fixed terms. That leads to:
\begin{align} \begin{bmatrix} C_{ff} & C_{fp} \end{bmatrix} \begin{bmatrix} V_f \\ V_p \end{bmatrix} &= 0 \\ \begin{bmatrix} C_{ff} \end{bmatrix} \begin{bmatrix} V_f \end{bmatrix} &= - \begin{bmatrix} C_{fp} \end{bmatrix} \begin{bmatrix} V_p \end{bmatrix} \\ \implies \begin{bmatrix} V_f \end{bmatrix} &= \begin{bmatrix} C_{ff} \end{bmatrix}^{-1} \left[ - \begin{bmatrix} C_{fp} \end{bmatrix} \begin{bmatrix} V_p \end{bmatrix} \right] \end{align}Thus, the equations can be solved using matrix inversion techniques.
1.3 Complete Example - Trapezoid (3 Elements)
Let's work through a simple triangular element example. (For ease of calculation, we leave out \(\epsilon\) from the coefficient matrix.) Take a look at the trapezoid below:
The numbers on the outside of the trapezoid are the global node indices, while the numbers on the inside corners of each triangle are the local node indices. Notice that the global node indices don't have to be ordered counter-clockwise.
Let's give some coordinates and initial values to the nodes:
| Node | Coordinate | Value \((V_i)\) |
|---|---|---|
| 1 | 1,1 | 100 |
| 2 | 0,0 | 100 |
| 3 | 2,1 | * |
| 4 | 1,0 | * |
| 5 | 3,0 | * |
1.3.1 Writing the element coefficient matrix
Since the coupling matrix is symmetric, we just need to set about writing the upper-right (or lower-left) diagonal:
\begin{align} C_{11} &= C_{11}^{(1)} + C_{11}^{(2)} \\ C_{12} &= C_{13}^{(1)} \\ C_{13} &= C_{12}^{(2)} \\ C_{14} &= C_{12}^{(1)} + C_{13}^{(2)} \\ C_{15} &= 0 \\ % C_{22} &= C_{33}^{(1)} \\ C_{23} &= 0 \\ C_{24} &= C_{23}^{(1)} \\ C_{25} &= 0 \\ % C_{33} &= C_{22}^{(2)} + C_{11}^{(3)} \\ C_{34} &= C_{23}^{(2)} + C_{13}^{(3)} \\ C_{35} &= C_{12}^{(3)} \\ % C_{44} &= C_{22}^{(1)} + C_{33}^{(2)} + C_{33}^{(3)} \\ C_{45} &= C_{32}^{(3)} \\ % C_{55} &= C_{22}^{(3)} \end{align}In matrix form, this is:
\begin{equation} C = \begin{bmatrix} C_{11}^{(1)} + C_{11}^{(2)} && C_{13}^{(1)} && C_{12}^{(2)} && C_{12}^{(1)} + C_{13}^{(2)} && 0 \\ && C_{33}^{(1)} && 0 && C_{23}^{(1)} && 0 \\ && && C_{22}^{(2)} + C_{11}^{(3)} && C_{23}^{(2)} + C_{13}^{(3)} && C_{12}^{(3)} \\ && && && C_{22}^{(1)} + C_{33}^{(2)} + C_{33}^{(3)} && C_{32}^{(3)} \\ && && && && C_{22}^{(3)} \end{bmatrix} \end{equation}The lower triangle is omitted for clarity, since the values are the same.
- Figuring the element areas
Element \(4A\) \(1\) \(-2\) \(2\) \(-2\) \(3\) \(-4\)
- Calculating the element coefficients
\(\textrm{Element} / (4A)\) \(11\) \(12\) \(13\) \(22\) \(23\) \(33\) \(1\) \(-0.5\) \(0.5\) \(0\) \(-1\) \(0.5\) \(-0.5\) \(2\) \(-1\) \(0.5\) \(0.5\) \(-0.5\) \(0\) \(-0.5\) \(3\) \(-1\) \(0.5\) \(0.5\) \(-0.5\) \(-0.5\) \(-0.25\)
- Calculating the complete coefficient matrix
\begin{equation} C = \begin{bmatrix} -1.5 & 0 & 0.5 & 1 & 0 \\ & -0.5 & 0 & 0.5 & 0 \\ && -1 & 0 & 0.5 \\ &&& -1.75 & -0.5 \\ &&&& -0.5 \end{bmatrix} \end{equation}A good thing to check here is for weak diagonal dominance. Diagonal dominance is the property of a matrix (usually referring to rows, although it doesn't matter in this case) wherein the magnitude of the diagonal element in each row is larger than the summed magnitude of all of the other row elements:
\begin{equation} \forall i : \left| C_{ii} \right| \geq \sum^5_{j = 1, j \neq i} \left| C_{ij} \right| \end{equation}The above matrix is diagonally dominant, which means that the iterative algorithm will work.
1.3.2 Iterative calculations
At this point, all we need to do is repeat the following equations:
\begin{align} V_3 &= - \frac{1}{C_{33}} \sum^5_{i = 1, i \neq 3} V_i C_{3i} \\ V_4 &= - \frac{1}{C_{44}} \sum^5_{i = 1, i \neq 4} V_i C_{4i} \\ V_5 &= - \frac{1}{C_{55}} \sum^5_{i = 1, i \neq 5} V_i C_{5i} \end{align}- Iterating
To check your work, you can compare with the below calculated table.
Iteration \(V3\) \(V4\) \(V5\) \(0\) \(0\) \(0\) \(0\) \(1\) \(50\) \(150\) \(0\) \(2\) \(50\) \(150\) \(-100\) \(3\) \(0\) \(200\) \(-50\) \(4\) \(25\) \(175\) \(-100\) \(5\) \(-50\) \(225\) \(-75\) \(6\) \(12.5\) \(112.5\) \(-137.5\) \(7\) \(-18.75\) \(81.25\) \(-50\) \(8\) \(25\) \(125\) \(-50\) \(9\) \(25\) \(125\) \(-50\)