Subsections

• 3.1.1 Approximation of $E_{\mathrm{n}}$ and $\varepsilon$
• 3.1.2 Continuation of the Discretization
• 3.1.3 Boundary Conditions
  • 3.1.3.1 Dirichlet Conditions
  • 3.1.3.2 Requirements for an M-Matrix
  • 3.1.3.3 Neumann Conditions
  • 3.1.3.4 Cauchy Boundary Conditions

3.1 The Poisson Equation

The Poisson equation is fundamental for all electrical applications. The derivation is shown for a stationary electric field [26]. For the derivation, the material parameters may be inhomogeneous, locally dependent but not a function of the electric field. In this section, the principle of the discretization is demonstrated. In the case of a low field dependency, a recursive reinsertion of newly derived material parameters may result in stable solutions [54]. A highly nonlinear behavior has to be managed by linearization and Newton iterations, for instance [2][4][60].

The Maxwell equations for a stationary electric field reduce to

$${\operatorname{curl}{\mathrm{\bf E}} } {= 0,}$$ (3.1)

$${\operatorname{div}{\mathrm{\bf D}} } {= \varrho .}$$ (3.2)

${\mathrm{\bf E}}$ and ${\mathrm{\bf D}}$ are related by the permittivity $\makebox{\boldmath $\underline\varepsilon$}$

$${\mathrm{\bf D}} = \makebox{\boldmath$\underline\varepsilon$}\;{\mathrm{\bf E}}.$$ (3.3)

Physical properties require that the permittivity tensor $\makebox{\boldmath $\underline\varepsilon$}$ is symmetric and positive definite for lossless media (or a positive scalar for isotropic media). Equation (3.1) is satisfied implicitly by introducing the electric potential $\phi$

$${\mathrm{\bf E}} = - \operatorname{grad}\phi.$$ (3.4)

Combing (3.1)-(3.4), the Maxwell equations can be rewritten as

$$\operatorname{div}\;( \makebox{\boldmath$\underline\varepsilon$}\;\operatorname{grad}\phi) = - \varrho .$$ (3.5)

This can be transformed to its integral formulation. By applying the Gaussian integral formula (3.5) results in

$$\int\limits_V \operatorname{div}\;( \makebox{\boldmath$\underline... ...}\phi) \cdot \mathrm{d}{\mathrm{\bf A}}= - \int\limits_V \varrho \;\mathrm{d}v.$$ (3.6)

Equation (3.6) must be fulfilled within any arbitrary volume $V$, with $\partial V$ being the surface of this volume. While performing Box Integration, this formula must be satisfied in the Voronoi boxes of each grid point. In Figure 3.1 a typical Voronoi box of a point $p_i$, belonging to all the tetrahedrons connected to the point $p_i$, can be seen. The integration can be split into the integration over the different box parts, caused by all the tetrahedrons $\mathrm{tets_i}$ which share the same point $p_i$

$$\sum_{t \;\in\; \mathrm{tets}_{i}} \;\int\limits_{A_{t,i}} (\make... ...i) \cdot \mathrm{d}{\mathrm{\bf A}}= - \int\limits_{V_i} \varrho \;\mathrm{d}v.$$ (3.7)

Such a part of the Voronoi box within the tetrahedron $< p_ip_jp_kp_l >$ is shown in Figure 3.2. Each part around the point $p_i$ can also be split into the three area contributions that are caused by the different edges $< p_ip_j >$, $< p_ip_k >$ and $< p_ip_l >$ (see Figure 3.3)

$$\underset{\forall\;< p_ip_j >}{\underset{t \;\in\;\text{tets}_i,}... ...i) \cdot \mathrm{d}{\mathrm{\bf A}}= - \int\limits_{V_i} \varrho \;\mathrm{d}v.$$ (3.8)

Figure 3.1: Tetrahedrons around point $p_i$ with the resulting Voronoi box constructed by these tetrahedrons.

Figure 3.2: Box parts of the tetrahedron $< p_ip_jp_kp_l >$ with drawn outer-sphere $s_{ijkl}$ of the tetrahedron and outer-circle $c_{ijk}$ of the triangle $< p_ip_jp_k >$.

Figure 3.3: Different Voronoi faces which share the point $p_i$ of a tetrahedron.

Using Box Integration, an approximation of first order is usually assumed, which means that the electric potential $\phi$ within one tetrahedron is a linear function along the axes or, in other words, the electric field ${\mathrm{\bf E}}$ is constant within each tetrahedron. Material parameters are assumed to be constant within the tetrahedron, too. The electric charge density $\varrho$ on the right-hand-side is assumed to be constant within the whole Voronoi box as well and can therefore be analytically calculated and written as the electric charge $q_i$ of the Voronoi box of the point $p_i$, which results in

$$-\underset{\forall\;< p_ip_j >}{\underset{t \;\in\;\mathrm{tets}_... ...mathrm{d}{\mathrm{\bf A}}= - \int\limits_{V_i} \varrho _i \;\mathrm{d}v= - q_i.$$ (3.9)

As in the usual applications the permittivity tensor reduces to a scalar, which means the electric displacement and the electric field have the same direction, it follows

$$(\varepsilon \;{\mathrm{\bf E}}) \cdot \mathrm{d}{\mathrm{\bf A}}... ...E}}\cdot \mathbf{n}) \;\mathrm{d}A = \varepsilon \;E_{\mathrm{n}} \;\mathrm{d}A$$ (3.10)

and equation (3.9) can be simplified to

$$-\underset{\forall\;< p_ip_j >}{\underset{t \;\in\;\mathrm{tets}_... ...mathrm{n}} \;\mathrm{d}A = - \int\limits_{V_i} \varrho _i \;\mathrm{d}v= - q_i.$$ (3.11)

3.1.1 Approximation of $E_{\mathrm{n}}$ and $\varepsilon$

Usually the integral (3.11) is approximated by

$$-\underset{\forall\;< p_ip_j >}{\underset{t \;\in\;\mathrm{tets}_i}{\sum}} \;\varepsilon _{ij} \;E_{ij} \;A_{t,ij} = - q_i,$$ (3.12)

with

$$E_{ij} \approx \frac{\phi_i - \phi_j}{d_{ij}} \qquad \varepsilon _{ij} \approx \frac{\varepsilon _i + \varepsilon _j}{2}.$$ (3.13)

Figure 3.4: A detail of the tetrahedron $< p_ip_jp_kp_l >$ shown in Figure 3.3 which shows the coupling Voronoi region between the point $p_i$ and $p_j$. This region can be split into the two triangular components $A_{ij,k}$ and $A_{ij,l}$ which are spawned by the midpoint $m_{ij}$ of the edge $< p_ip_j >$, the midpoint $m_{ijk}$ of the outer-circle of the triangle $< p_ip_jp_k >$ and the midpoint $m_{ijkl}$ of the outer-sphere of the tetrahedron $< p_ip_jp_kp_l >$, or $m_{ij}$, $m_{ijl}$, and $m_{ijkl}$, respectively.

Figure 3.5: Potential distribution between two grid points $p_i$ and $p_j$ with different permittivities.

An estimation of the validity of this approximation based on the homogenous Poisson equation (Laplace equation) will be explained under the conditions shown in Figure 3.4. $(\xi ,\eta ,\zeta )$ spans a local Cartesian coordinate system with origin $p_i$. In the one-dimensional case and with constant permittivity the potential varies linearly with position as

$$\phi(\xi ) = \phi_i - E_{ij} \; \xi .$$ (3.14)

Whereas with non-constant permittivity we have at first order

$$\widetilde{\varepsilon }(\xi ) = \varepsilon _i + \frac{\Delta \varepsilon }{d_{ij}} \;\xi ,$$ (3.15)

with

$$\widetilde{\varepsilon }(0)=\varepsilon _i \qquad \widetilde{\varepsilon }(d_{ij})=\varepsilon _{j}$$ (3.16)

and therefore,

$$\Delta \varepsilon =\varepsilon _j-\varepsilon _i.$$ (3.17)

The continuity of the electric displacement

$$\widetilde{\varepsilon }(0) \;\widetilde{E}(0)= \widetilde{\varep... ...;\widetilde{E}(\xi )= \widetilde{\varepsilon }(d_{ij}) \;\widetilde{E}(d_{ij}),$$ (3.18)

with

$$\widetilde{E}=-\frac{\partial \widetilde{\phi}}{\partial \xi }$$ (3.19)

delivers the deformed distribution of $\widetilde{\phi}(\xi )$

$$\phi_i - \widetilde{\phi}(\xi ) = (\phi_i - \phi_j)\; \frac{\disp... ...)} {\displaystyle\ln \left(1+\frac{\Delta \varepsilon }{\varepsilon _i}\right)}$$ (3.20)

and therefore,

$$\widetilde{E}\left(\frac{d_{ij}}{2}\right) = \left. -\frac{\parti... ...vert _{\textstyle \frac{d_{ij}}{2}} \neq \frac{\phi_i - \phi_j}{d_{ij}}=E_{ij}.$$ (3.21)

The resulting potential distributions $\widetilde{\phi}(\xi )-\phi_i$ between $p_i$ and $p_j$ for different values of $\Delta\varepsilon$ are shown in Figure 3.5. The overestimation of the electric field $E_{ij}/\widetilde{E}(d/2)$ in the middle of the grid element in dependence of the permittivity change inside the element is shown in Figure 3.6. According to this figure, a limit for the validity of approximation (3.13) can be defined by

$$\vert\varepsilon _i - \varepsilon _j\vert < k \;\mathrm{min}(\varepsilon _i,\varepsilon _j) \qquad \forall \; i,j \;:\; \exists \; \;< p_ip_j >$$ (3.22)

or a limit for continuously defined permittivity

$$\vert\operatorname{grad}\varepsilon _i \vert \;d_{ij} < k\;\varepsilon _i \qquad \forall \; i,j \;:\; \exists \; \;< p_ip_j >.$$ (3.23)

With $k\approx 2$, the estimation error of $E_{ij}$ is limited to 10%.

Figure 3.6: Ratio of the approximated electric field $E_{ij}$ to the exact solution $\widetilde{E}(d_{ij}/2)$ at the center of the control volume.

On a three-dimensional grid (refer to Figure 3.4), the spatial derivatives of $\varepsilon$ cause an additional contribution to the discretization error. With linear interpolation of $\varepsilon$ and $\phi$, equation (3.11) results in

$$\begin{gather*}\begin{split}-&\underset{\forall\;< p_ip_j >}{\underset{t \;\in\;... ... }{d_{ij}} \;\mathrm{d} \eta \;\mathrm{d} \zeta = - q_i \end{split}\end{gather*}$$ (3.24)

and with $L H = 2 A_{ij}$

$$\begin{gather*}\begin{split}&\underset{\forall\;< p_ip_j >}{\underset{t \;\in\;\... ...\Delta\phi_\xi \frac{A_{ij}}{d_{ij}} \frac{H}{3} = -q_i.\end{split}\end{gather*}$$ (3.25)

The first addend of equation (3.25) equals the value given by discretization (3.12) and (3.13) whereas the remaining terms cause an additional discretization error. Often the derivative components of neighboring tetrahedrons compensate each other, but not always. Even the lengths of $L$ and $H$ are not limited to lie within the tetrahedron ranges (for instance, obtuse angled triangles for two dimensions) and a limitation for this approximation is difficult. For the whole simulation domain the estimation

$$\vert\operatorname{grad}\varepsilon _i\vert _{\mathrm{max}} \;d_{... ...max}} < k\;\varepsilon _{i, \mathrm{min}}\qquad \forall \; i,j \;:\; \exists \; \;< p_ip_j >$$ (3.26)

seems plausible, but is too pessimistic as in most cases $L$ and $H$ have the same order of magnitude as their based tetrahedrons. The previously derived estimation (3.23) is therefore also applicable to the three-dimensional case.

Another common approximation method for the effective permittivity $\varepsilon _{ij}$ is [27]

$$\varepsilon _{ij} \approx \frac{2 \;\varepsilon _i \;\varepsilon _j}{\varepsilon _i + \varepsilon _j}$$ (3.27)

which shows the behavior of a serial connection of two capacitors (or resistors for electrical flow fields). This approximation has the advantage that if the permittivity (more likely for the conductivity) is zero in a Voronoi box, no flux enters the box. Nevertheless, this zero flux leaves the electric potential of the included grid point undefined and may cause problems for the equation solver, applied later on.

3.1.2 Continuation of the Discretization

However, equation (3.11) is discretized as

$$-\underset{\forall\;< p_ip_j >}{\underset{t\;\in\;\mathrm{tets}_i... ...sum}} \varepsilon _{ij} \;\frac{\phi_{j} - \phi_{i}}{d_{ij}} \;A_{t,ij} =- q_i.$$ (3.28)

In compliance with the decomposition of the Voronoi regions in Section 2.4.3, $A_{t,ij}$ is the right signed value of the Voronoi area between $p_i$ and $p_j$ of the tetrahedron $t$. Which means in detail, area portions with

• outer-circle inside the appropriate triangle and outer-sphere inside the tetrahedron counted positive,
• outer-circle outside the appropriate triangle and outer-sphere outside the tetrahedron counted positive, and
• one of them inside and the other one outside counted negative.

Combining the portions of the tetrahedrons which share the point $p_i$ delivers the commonly used notation

$$\sum_{\forall\; j\; : \;\exists\; \text{edge}\;< p_ip_j >} \varepsilon _{ij} \;\frac{A_{ij}}{d_{ij}} \;(\phi_{i} - \phi_{j}) = q_{i}$$ (3.29)

with $A_{ij}$ the Voronoi area of the edge $< p_ip_j >$ (coupling area between $p_i$ and $p_j$)

$$A_{ij}=\underset{< p_ip_j > \;\in\; t}{\underset{t\;\in\;\mathrm{tets}_i}{\sum}} A_{t,ij}.$$ (3.30)

Equation (3.29) must be fulfilled for each grid point $p_i$ inside the simulation domain, which results in an equation system

$${\mathrm{\bf S}} \cdot {\mathrm{\bf x}} = {\mathrm{\bf B}}$$ (3.31)

with the system matrix ${\mathrm{\bf S}}=(s_{ij})$

$$\qquad\qquad\qquad\qquad s_{ij} = - \varepsilon _{ij} \;\frac{A_{ij}}{d_{ij}} \qquad\forall \;i,j \; :\;\exists\;\mathrm{edge}\;< p_i p_j >,\qquad\qquad$$ (3.32)

$$s_{ii} = \sum_{\forall\;j\;:\;\exists\;\mathrm{edge}\;< p_i p_j >} \varepsilon _{ij} \;\frac{A_{ij}}{d_{ij}}=-\sum_{\forall\;j\;:\;i \neq j} s_{ij} \qquad\forall \; i,$$ (3.33)

the unknowns

$${\mathrm{\bf x}}=(\phi_i),$$ (3.34)

and the right-hand side

$${\mathrm{\bf B}}=(q_i),$$ (3.35)

which is defined for all inner points $p_i$ of the domain.

3.1.3 Boundary Conditions

For solving this time-invariant elliptical partial differential equation, boundary conditions have to be applied to all boundaries. For each grid point lying at the boundary of the domain one of the two following additional equations or a linear-combination of them must be added for the discrete version, which fulfill the boundary conditions.

Each grid point of the boundary has to be defined by exactly one of the different conditions. Additionally, at least one grid point must be described by a Dirichlet condition to define a definite potential distribution. With this set of equations, the discrete Poisson equation can be solved.

3.1.3.1 Dirichlet Conditions

The first class of boundary conditions is an imposed potential $\Phi$

$$\phi=\Phi \partial B_\Phi$$ (3.36)

The discrete formulation of this boundary condition simply replaces the respective rows of equation (3.29) by

$$\phi_{i}=\Phi_i\qquad\forall \;i \;:\; p_i \in \partial B_\Phi$$ (3.37)

or within the matrix representation

$$s_{ij}=0\qquad\forall\;i,j \;:\; i\neq j,\;p_i \in \partial B_\Phi,$$ (3.38)

$$s_{ii}=1\qquad\forall\;i\;:\;p_i \in \partial B_\Phi,$$ (3.39)

$$b_i=\Phi_i\qquad\forall\;i\;:\;p_i \in \partial B_\Phi,$$ (3.40)

which at first delivers an asymmetric system matrix ${\mathrm{\bf S}}$. Since most equation solvers require a symmetric system matrix, it may be necessary to transform the matrix to a symmetric one, which can be simple obtained by summation of adequate rows.

3.1.3.2 Requirements for an M-Matrix

As shown in Section 2.3, for satisfying the maximum principle the system matrix ${\mathrm{\bf S}}$ has to be an M-matrix. Under the precondition of a valid Delaunay grid, the validity of the M-matrix will be shown.

As the grid is Delaunay, it follows that the coupling areas are positive

$$A_{ij}\geq 0$$

and therefore,

	(3.32)	$s_{ij}$	$\leq 0$		satisfies (2.4),
	(3.33)	$-\sum_{\forall\;j\;:\;i\neq j} s_{ij} = s_{ii}$	$\geq 0$		satisfies (2.5),
	(3.33)	$-\sum_{\forall\;j\;:\;i\neq j} s_{ij} = s_{ii}$	= $\sum_{\forall\;j\;:\;i\neq j} \vert s_{ij}\vert$		satisfies (2.6),

and finally by defining at least one boundary point $p_{\Phi,k}$ as a Dirichlet point

(3.38)-(3.39)

$1=s_{kk}$

> $\sum_{\forall\;j\;:\;k\neq j} \vert s_{kj}\vert = 0$

satisfies (2.7).

3.1.3.3 Neumann Conditions

Figure 3.7: A Voronoi box of the point $p_i$ at a boundary.

The second type of boundary condition specifies the imposed normal component ${\mathrm{\bf n}} \cdot {\mathrm{\bf D}} = D_{n}$ on the boundaries $\partial B_n$ (see Figure 3.7), which occurs between different segments. As shown in Figure 3.7, the normal vector ${\mathrm{\bf n}}$ may have various directions on different boundary pieces around $p_i$. This behavior could be handled by splitting the boxes (and also the point $p_i$ itself) into several pieces. Usually this approach is prevented and the normal derivatives are combined to one value. In this case, the boundary conditions result in a distortion of (3.29)

$$\sum_{\forall\;j\;:\;\exists\;\mathrm{edge} \;< p_ip_j >} \vareps... ...\phi_{j}) + D_{{n},i} \;A_i = q_{i}\qquad\forall\; i \;:\; p_i \in \partial B_n$$ (3.41)

with $D_{{n},i}$ the assigned normal component at the boundary at the grid point $p_i$. It must be distinguished between $A_{ij}$ which represents the area of the Voronoi box between $p_i$ and $p_j$, and $A_i$ as the assigned boundary area of the Voronoi box around the grid point $p_i$.

In the majority of cases, a special kind of this boundary condition is used, which shows a vanishing flux across the boundary (homogenous Neumann condition)

$$D_{\mathrm{\bf n}} = 0.$$ (3.42)

This type of condition specifies the artificial boundaries which terminate the simulation domains. As no distortion of equation (3.41) appears, these boundaries do not require a special treatment and the equations can be left unchanged as in (3.31)-(3.35).

3.1.3.4 Cauchy Boundary Conditions

The linear-combination

$$D_{\mathrm{\bf n}} + \alpha \;\phi= \beta$$ (3.43)

combines the two previous conditions and results in the discrete formulation

$$\sum_{\forall\;j\;:\;\exists\;\mathrm{edge} \;< p_ip_j >} \vareps... ...- \alpha_i\;\phi_i) \;A_i = q_{i}\qquad\forall \; i \;:\; p_i \in \partial B_c.$$ (3.44)