MECSEM Electrostatic Field Solver

This repository implements a generalized version of the electrostaticFoam solver to model linear dielectrics. The solver is restricted so that there are no free charges in the volume.

Contributors: Cristian & David

Theory

The solver solves the equation

$$ \nabla\cdot(\textbf{k}\nabla\phi) = 0 $$

where:

• $\textbf{k}$ is the dielectric constant of the medium, i.e., a scalar field
• $\phi$ is the electric potential

The primary boundary conditions of the problem are set on the $\phi$ field. I.e., setting a fixed voltage for conductors or zeroGradient BC for atmosphere.

The electric field can then be obtained as

$$ \vec{\textbf{E}} = -\nabla\phi $$

For a linear dielectric, the electric displacement $\vec{\textbf{D}}$ is defined as

$$ \vec{\textbf{D}} = \epsilon_0\vec{\textbf{E}} + \vec{\textbf{P}} = \textbf{k} \epsilon_0 \vec{\textbf{E}} $$

where:

• $\epsilon_0$ is the vacuum permittivity
• $\vec{\textbf{E}}$ is the electric field
• $\vec{\textbf{P}}$ is the polarization

The polarization density is related to the bound charge density ($\rho_\textbf{b}$) through:

$$ \rho_\textbf{b} = -\nabla\cdot\vec{\textbf{P}} $$

We can rearrange the displacement equation:

$$ \begin{align*} \epsilon_0\vec{\textbf{E}} + \vec{\textbf{P}} &= \textbf{k} \epsilon_0 \vec{\textbf{E}} \\ \vec{\textbf{P}} &= (\textbf{k}-1) \epsilon_0 \vec{\textbf{E}} \\ \nabla\cdot\vec{\textbf{P}} &= \nabla\cdot\left[(\textbf{k}-1) \epsilon_0 \vec{\textbf{E}}\right] \\ \rho_\textbf{b} &= \epsilon_0\nabla\cdot\left[(1-\textbf{k}) \vec{\textbf{E}}\right] \end{align*} $$

Thus obtaining the bound charge density field.

Boundary conditions must also be imposed on the electric field. The tangential electric field is zero at the surface of all conductors. The dielectric constant is zeroGradient at all boundaries.

Test Case

A simple two plate capacitor is used to verify the solver. I have stumbled across an excellent Feynman lecture that covers the relevant analytical theory.

As shown in the lecture (10.6):

$$ E = \frac{\sigma_f - \sigma_b}{\epsilon_0} $$

Furthermore (10.9):

$$ \sigma_f = \epsilon_0 E k $$

Thus,

$$ \begin{align*} E &= \frac{\epsilon_0 E k - \sigma_b}{\epsilon_0} \\ \epsilon_0 E &= \epsilon_0 E k - \sigma_b \\ \sigma_b &= \epsilon_0 (k-1) E \end{align*} $$

But, for an 'infinite' pair of parallel plates at potential difference $V$ spaced $d$ units apart, $E = V/d$. Thus,

$$ \sigma_b = \frac{\epsilon_0 (k-1) V}{d} $$

Logically this makes sense. If $k=1$ (vacuum) there are no bound charges, but there are still some free charges built up on the conductor. If $k>1$, a dielectric is present and will have bound charges.

So, we can check the result of the simulation against this analytical result by looking at the bound charge density on the plates and confirming that the result is accurate.

Setup

For the test case, two 100x100mm plates, 1mm thick are spaced 4mm apart. A 1kV potential difference is applied to the plates. The dielectric comprises the 100x100x4mm volume within the two plates and is set to have a dielectric constant of 13. The entire setup is enclosed within a 200x60x200mm void. Approximately 30 million cells are used.

Results

Plots

1. Mesh Closeup
2. Electric Potential Plot
3. Potential Plot with Field Lines
4. Electric Field Magnitude Plot
5. Bound Volume Charge Density Plot

Analysis

The mean $\rho_b$ in the dielectric-facing surface cells of the plates is approximately $0.427937,\text{C/m}^3$. The height of these cells is $0.06207$ mm. Thus, the predicted bound surface charge density is $2.65621\cdot10^{-5},\text{C/m}^2$. The analytical result is $2.65626\cdot10^{-5},\text{C/m}^2$, implying an error of around 0.002%.