Subsections

• 5.3.1 Two-Dimensional Behavior
• 5.3.2 Extension to Three Dimensions
  • 5.3.2.1 Prerequisites:
  • 5.3.2.2 Methodology of the Capacitor Model
  
  
• 5.3.3 Calculating the Equipotential Surfaces of the Capacitor
• 5.3.4 Topological Cuboids
• 5.3.5 Proof of the Two-Dimensional Duality
  • 5.3.5.1 Solution with a Scalar Potential $\pmb\phi$
  • 5.3.5.2 Solution with an (Electric) Vector Potential ${\mathrm{\bf A}}^e$
  • 5.3.5.3 Summary of Both Solutions
  • 5.3.5.4 Field Lines of ${\mathrm{\bf E}}$
  
  
• 5.3.6 The Three-Dimensional Capacitor Model
• 5.3.7 Evaluation of the Electric Potential of the Capacitor
  • 5.3.7.1 Why Choosing a Conjugate Gradient Solver

5.3 Potential-Based Grid Generation

In the particular case of MOS structures, a point distribution with low point density along the channel and high point density towards the bulk seems feasible. Nevertheless, to fulfill the Delaunay criterion, the grid generator sometimes has to insert additional grid points. Therefore, a method for point insertion which does not violate the Delaunay criterion must be found.

For ortho-product grids, the aspect ratio of the cuboids can be varied arbitrarily. As mentioned before, an ortho grid can be transformed easily to a tetrahedral grid by splitting each cuboid into five or six tetrahedrons. The Delaunay criterion will never be violated, since all resulting tetrahedrons have the same outer-sphere and neighboring grid points can never lie within. A pleasant side effect of this method is that the coupling areas of the tetrahedrons inside the original parallelepiped will become zero, the resulting equation system will have fewer off-diagonal entries and the assembled equation system will be the same as the one obtained from the original ortho grid.

However, this approach is difficult when cuboids have to be fitted along geometry lines that are not in the same direction as the edges of the cuboids (refer Section 2.1). For maintaining the accuracy, the number of cuboids and, therefore, also the number of grid points will increase rapidly. In this case, the general use of tetrahedral grids is required. There are fewer geometry related limitations and a Delaunay grid can be found in most cases. But even an enhancement of the grid density will cause the tetrahedrons to become nearly equilateral (see Figure 5.5(b)) and a directional dependent grid density cannot be tuned arbitrarily. In the method developed here, the advantages of ortho and tetrahedral grids are combined.

5.3.1 Two-Dimensional Behavior

The basic idea of the newly developed method can be explained by considering a nearly cuboidal capacitor with two electrodes. The electric field outside the capacitor is neglected by assuming a high-k dielectric with a permittivity $\varepsilon$ much higher than the vacuum permittivity. A sketch of such a capacitor is shown in Figure 5.6.

Figure 5.6: Two-dimensional capacitor with the electrodes and drawn field and equipotential lines.

After applying a voltage between the electrodes, the equipotential lines and the field lines can be evaluated. The field lines are orthogonal to the potential lines. Additionally, either the equipotential or field lines are boundary-conform or orthogonal to the boundaries of the capacitor. If the spacing of selected potential and field lines of the capacitor is sufficiently small, two potential and two field lines lying next to each other approximately form rectangles, as shown in the figure. The spacing between the selected field and potential lines can be varied independently. To create a grid suitable for device simulation, these rectangles must be split into triangles by a mesh generator.

This method can be used for creating a set of grid points in a given domain. Four points at the boundary of this domain are placed at the corners of the capacitor. Between the corners, the upper and lower electrodes are located. After the evaluation of the field distribution, the grid points can be set at arbitrary intersections of the equipotential lines and field lines.

The field lines are not evaluated directly. Rather, they are evaluated by a duality, which can be applied if the following prerequisites are satisfied:

• Constant permittivity $\varepsilon$ inside the capacitor.
• Four corner points are placed on the boundary.
• Between each two corners the electrodes are placed where the voltage is connected (Dirichlet boundary conditions).
• The remaining opposite boundaries show a vanishing normal component of the electric field (homogenous Neumann conditions).

After applying a voltage between the electrodes, the electric potential distribution ( ${\mathrm{\bf E}}=- \operatorname{grad}\phi$) inside the capacitor can be calculated and the equipotential lines ( $\phi= \mathrm{const}$) can be evaluated.

The field lines are evaluated via a duality in the specified way:

• The dual capacitor, with congruent shape, replaces the original capacitor.
• Constant permittivity $\varepsilon$ inside the capacitor.
• The same four corner points on the boundary.
• The original boundaries with vanishing normal component of the field are replaced by electrodes -- with imposed voltage.
• The original electrodes are replaced by boundaries with vanishing normal component of the electric field.
• The field lines of the original capacitor are the equipotential lines of the dual capacitor.

The dual electric field ( ${\mathrm{\bf F}}=- \operatorname{grad}\psi$) and its associated equipotential lines ( $\psi = \mathrm{const}$) are calculated.

Now each point $(x,y)$ inside the capacitor can be represented by its potential representation $(\phi,\psi)$.

5.3.2 Extension to Three Dimensions

In the two-dimensional case, there are four electrodes at the boundary of the capacitor -- two opposite electrodes of the original capacitor and the remaining two opposite boundaries as the electrodes for the dual capacitor. The extension of this method to three dimensions is as:

5.3.2.1 Prerequisites:

• Constant permittivity $\varepsilon$ inside the capacitor.
• There are exactly six boundary groups -- ``sides'',
• Which must not overlap and
• In general, all six ``sides'' must cover the whole surface.
• Two opposite ``sides'' act as electrodes,
• The other four ``sides'' act as Neumann boundaries.

Figure 5.7: Capacitor with electrodes at top and bottom.

Such a capacitor can be seen in Figure 5.7. The original capacitor delivers the field distribution ${\mathrm{\bf E}}=- \operatorname{grad}\phi$ with the equipotential surfaces $\phi= \mathrm{const}$.

5.3.2.2 Methodology of the Capacitor Model

The electric field ${\mathrm{\bf F}}$ and its electric potential $\psi$, and as a consequence of this three-dimensional expansion, the electric field ${\mathrm{\bf G}}$ and the potential $\vartheta$, are derived as the resulting field inside a capacitor of the following shape:

• The new capacitor region is the same as the original one.
• Constant permittivity.
• Two opposing Neumann surfaces are replaced by electrodes where a voltage is connected.
• The former electrodes are replaced by homogenous Neumann boundaries.
• The remaining two surfaces are left as homogenous Neumann surfaces.

With this new placement of the electrodes, the electric potential ${\mathrm{\bf F}}=- \operatorname{grad}\psi$ can be derived. By an additional replacement of the remaining two ``sides'' by electrodes and all others as insulator surfaces, the third potential ${\mathrm{\bf G}}=-\operatorname{grad}\vartheta$ can be calculated.

5.3.3 Calculating the Equipotential Surfaces of the
Capacitor

Inside the capacitor, the Laplace equation has to be solved (as shown in Section 3.1). The boundary conditions are

$$\phi(x,y,z)=\cases {\Phi_0}{\text{for $(x,y,z)\in \partial V_1$}} {\Phi_1}{\text{for $(x,y,z)\in \partial V_2$},}$$ (5.16)

for the electrodes where the voltage is connected and vanishing normal component

$$\partial_n \phi(x,y,z)=\cases {0}{\text{for $(x,y,z)\in \partial ... ...V_4$}} {0}{\text{for $(x,y,z)\in \partial V_5$ or $(x,y,z)\in \partial V_6$},}$$ (5.17)

on the other two opposing boundary groups. With this set of boundary conditions, the equipotential faces of the electric field can be evaluated.

It can be seen that the elements defined by two appropriate equipotential lines of $\phi$, $\psi$, and $\vartheta$ are spanning the desired cuboids for the tetrahedrization. The different electrode placements with their assigned field distribution are illustrated in Figure 5.8 (original distribution), 5.9 (first electrode-replacement), and 5.10 (second electrode-replacement). Resulting equipotential surfaces of all three field calculations are shown in Figure 5.11. The full set of differential equations (assuming constant permittivity) with their assigned boundary conditions can be seen in Table 5.1.

Figure 5.8: Cuboidal capacitor, original electrodes, potential distribution.

Figure 5.9: Cuboidal capacitor, first electrode-replacement, potential distribution.

Figure 5.10: Cuboidal capacitor, second electrode-replacement, potential distribution.

Figure 5.11: Cuboidal capacitor, with the equipotential surfaces of the three electric fields that are spawning the cuboids.

Table 5.1: Differential equations and boundary conditions of the capacitor and its dual capacitors.

$\operatorname{div}{\mathrm{\bf E}}= 0$	$\operatorname{div}{\mathrm{\bf F}}= 0$	$\operatorname{div}{\mathrm{\bf G}}= 0$
$\operatorname{curl}{\mathrm{\bf E}}= {\mathrm{\bf0}}$	$\operatorname{curl}{\mathrm{\bf F}}= {\mathrm{\bf0}}$	$\operatorname{curl}{\mathrm{\bf G}}= {\mathrm{\bf0}}$
${\mathrm{\bf E}}= - \operatorname{grad}\phi\qquad$	${\mathrm{\bf F}}= - \operatorname{grad}\psi\qquad$	${\mathrm{\bf G}}= - \operatorname{grad}\vartheta \quad$
$\phi=\Phi_1$	$\partial_n \psi=0$	$\partial_n \vartheta =0$	: for $(x,y,z) \in \partial V_1$
$\phi=\Phi_2$	$\partial_n \psi=0$	$\partial_n \vartheta =0$	: for $(x,y,z) \in \partial V_2$
$\partial_n \phi=0$	$\psi=\Psi_3\qquad$	$\partial_n \vartheta =0$	: for $(x,y,z) \in \partial V_3$
$\partial_n \phi=0$	$\psi=\Psi_4\qquad$	$\partial_n \vartheta =0$	: for $(x,y,z) \in \partial V_4$
$\partial_n \phi=0\qquad$	$\partial_n \psi=0$	$\vartheta =\Theta_5\qquad$	: for $(x,y,z) \in \partial V_5$
$\partial_n \phi=0\qquad$	$\partial_n \psi=0$	$\vartheta =\Theta_6\qquad$	: for $(x,y,z) \in \partial V_6$

5.3.4 Topological Cuboids

In general, a realistic devices are more complex than such capacitors and therefore, the whole simulation area must be split into segments of usually different materials. Not all of them will look like cuboids, either. However, if the following conditions are fulfilled, the segments can be treated as warped cuboids. In this case, the calculation of the equipotential surfaces can be performed as described.

Like the sides of a cuboid, it must be possible to define such six ``sides'' for the segment -- an upper, lower, left, right, top, and bottom ``side''. Usually the semiconductor segment has such a shape: The simulation domain of the silicon wafer with planar surfaces, except the upper face with the oxide-interface, for instance. In this case it is easy to extract the six ``sides'' of the warped cuboid. Even other regions with four surrounding flat faces and two warped faces, each at the top and the bottom, can be handled easily. Additionally other criteria for finding the ``sides'' of the cuboids are conceivable, even though there might be problems when these ``sides'' have to be detected automatically.

Such a segment with warped six ``sides'' gives a topological cuboid. By applying a voltage at two of the opposing ``sides'', the electric field distribution can be evaluated. By connecting a voltage to each of the remaining opposite side pairs and also by calculating the potential distribution we get three potential distributions.

With this methodology, the topological cuboids of the grid can be calculated. Doing so, the distances between the selected potential ticks can be tuned arbitrarily and the aspect ratios of the resulting cuboids can be varied independently.

To achieve a valid Delaunay grid, the grid points produced by this method are inserted in the global point set of the geometry and are meshed by the grid generator. The demands on Delaunay grids are nearly fulfilled by the orthogonality of these cuboids. As a result of the given boundary points, a distortion of these cuboids may be caused.

5.3.5 Proof of the Two-Dimensional Duality

Figure 5.12: Domain for the electric field calculation.

First we start with the stationary field equations. With constant permittivity $\varepsilon$ the equations reduce to

$${\operatorname{div}{\mathrm{\bf E}}} {= 0}$$ (5.18)

$${\operatorname{curl}{\mathrm{\bf E}}} {= {\mathrm{\bf0}}}$$ (5.19)

which must be satisfied in the capacitor domain $B$, shown in Figure 5.12. The contacts of the capacitor are the boundaries

$$\mathcal{C}_{2i}, \qquad i=1(1)n,$$ (5.20)

which alternate with the Neumann boundaries

$$\mathcal{C}_{2i-1}, \qquad i=1(1)n.$$ (5.21)

The tangent and normal vectors ${\mathrm{\bf t}}$ and ${\mathrm{\bf n}}$ are defined as

$$\n =\t\times {{{\mathrm{\bf e}}_z}}$$ (5.22)

and

$$\t ={{{\mathrm{\bf e}}_z}}\times \n .$$ (5.23)

5.3.5.1 Solution with a Scalar Potential $\pmb\phi$

The first way for solving the partial differential equation system (5.18) and (5.19) is performed by a gradient field

$${\mathrm{\bf E}}= -\operatorname{grad}\phi,$$ (5.24)

which satisfies (5.19) implicitly and delivers the differential equation

$$\operatorname{div}\operatorname{grad}\phi= \Delta \phi= 0 B$$ (5.25)

with the boundary conditions

$$\phi = \Phi_{2i}=\mathrm{const}\qquad \mathcal{C}_{2i}, \quad i= 1(1)n$$ Dirichlet: (5.26)

$$\qquad \partial_{n}\phi = 0 \qquad \mathcal{C}_{2i-1}, \quad i= 1(1)n$$ (5.27)

The boundary problem (5.25), (5.26) and (5.27) has a well-defined solution. Additionally the following integrals will be defined:

$$\G _{2i}= \int\limits_{\mathcal{C}_{2i}} E_n \; \mathrm{d}s=-\int\limits_{\mathcal{C}_{2i}} \partial_{n}\phi\; \mathrm{d}s, \qquad i=1(1)n$$ (5.28)

and the trivial portions on the homogenous Neumann boundaries

$$\G _{2i-1}=\int\limits_{\mathcal{C}_{2i-1}} E_n \; \mathrm{d}s= -\int\limits_{\mathcal{C}_{2i-1}} \partial_{n}\phi\; \mathrm{d}s=0,\qquad i=1(1)n.$$ (5.29)

Integration of (5.25) over the whole region $B$ and applying Gauss' integral theorem delivers

$$0=\int\limits_{B} \Delta \phi\; \mathrm{d}A=\int\limits_{B} \oper... ...int\limits_{\partial B}\partial_{n}\phi\;\mathrm{d}s = -\sum_{i=1}^{n} \G _{2i}$$ (5.30)

5.3.5.2 Solution with an (Electric) Vector Potential ${\mathrm{\bf A}}^e$

Different to the previous way, the basic approach is a vector potential

$${\mathrm{\bf E}}= \operatorname{curl}{\mathrm{\bf A}}^e, \qquad {\mathrm{\bf A}}^e= {\psi}(x,y) \;{{{\mathrm{\bf e}}_z}}$$ (5.31)

which satisfies (5.18). Because of the two-dimensionality of the problem, the vector potential shows only a z-component. Transformation delivers

$$\qquad\qquad\qquad\qquad\qquad {\mathrm{\bf E}} = \operatorname{curl}({\psi}\;{{{\mathrm{\bf e}}_z}})=\nabla \times ({\psi}\;{{{\mathrm{\bf e}}_z}}) = - {{{\mathrm{\bf e}}_z}}\times \nabla {\psi}$$

$$= \operatorname{grad}{\psi}\times {{{\mathrm{\bf e}}_z}}= \partial_y {\psi}\;{\mathrm{\bf e}}_x - \partial_x {\psi}\;{\mathrm{\bf e}}_y$$ (5.32)

$$\t\cdot{\mathrm{\bf E}} = \t\cdot ( \operatorname{grad}{\psi}\times {{{\mathrm{\bf e}}_z}... ... e}}_z}}\times \t ) = - \operatorname{grad}{\psi}\cdot \n = -\partial_{n}{\psi}$$ (5.33)

$$\n\cdot{\mathrm{\bf E}} = \n\cdot ( \operatorname{grad}{\psi}\times {{{\mathrm{\bf e}}_z}... ...\bf e}}_z}}\times \n ) = \operatorname{grad}{\psi}\cdot \t = \partial_{t}{\psi}$$ (5.34)

$$\operatorname{curl}{\mathrm{\bf E}} = \nabla \times (\operatorname{grad}{\psi}\times {{{\mathrm{\bf e... ...ratorname{grad}{\psi}- {{{\mathrm{\bf e}}_z}}(\nabla \operatorname{grad}{\psi})$$

$$= \underbrace{\partial_z \operatorname{grad}{\psi}}_{\displaystyl... ...}(\nabla \nabla {\psi}) = -{{{\mathrm{\bf e}}_z}}\Delta {\psi}= {\mathrm{\bf0}}$$ (5.35)

$$\Longrightarrow \qquad \Delta {\psi}= 0 B.$$ (5.36)

The boundary conditions (5.26) and (5.27) result to

$$\begin{gather*}\begin{alignat}{1} &\mathcal{C}_{2i}:\qquad \phi=\mathrm{const}\q... ...ightarrow}\quad \text{\bf Dirichlet conditions}\notag \end{alignat}\end{gather*}$$

The previously defined integral is carried out as

$$\int\limits_\mathcal{C}E_n \;\mathrm{d}s= \int\limits_\mathcal{C}... ...s_A^E \mathrm{d}\psi= \left. {\psi}\right\vert _E - \left. {\psi}\right\vert _A$$ (5.37)

and further on by defining the arbitrary points

$$Q_{2i-1}(s_{2i-1}) \in \mathcal{C}_{2i-1}$$ (5.38)

follows

$$\int\limits_{Q_{2n-3}}^{Q_{2n-1}} E_n \;\mathrm{d}s=\int\limits_{... ...) = -\int\limits_{s_{2i-3}}^{s_{2i-1}} \partial_{n}\phi\;\mathrm{d}s= \G _{2i}.$$ (5.39)

or

$$\Psi_{2i-1}=\sum_{j=1}^{i-1} \G _{2j} + {\Psi}_1, \qquad i=1(1)n$$ (5.40)

where the value of ${\Psi}_1$ can be chosen arbitrarily.

5.3.5.3 Summary of Both Solutions

$$\begin{gather*}\begin{align}\text{I:}&&\Delta \phi&= 0 \qquad &&\text{in } B\ &... ...\text{on } \partial B_{2i-1}, \quad i=1(1)n\qquad\notag \end{align}\end{gather*}$$

and

$$\begin{gather*}\begin{align}\text{II:}&&\Delta \psi &= 0 \qquad &&\text{in } B\\... ...uad &&\text{on } \partial B_{2i-1}, \quad i=1(1)n\notag \end{align}\end{gather*}$$

which describe the same fields, if $\Psi_{2i-1}$ satisfy (5.28) and (5.43).

5.3.5.4 Field Lines of ${\mathrm{\bf E}}$

The field lines of ${\mathrm{\bf E}}$ are defined by the differential equation

$$\frac{\mathrm{d}y}{\mathrm{d}x} = \frac{E_y}{E_x}.$$ (5.41)

By insertion of (5.32) it follows

$$0=\mathrm{d}y\;E_x - \mathrm{d}x\;E_y = \mathrm{d}y\;\partial_y {... ...ial_x {\psi}= \mathrm{d} {\psi}\quad\Longrightarrow\quad {\psi}= \mathrm{const}$$ (5.42)

which are the equipotential lines of ${\psi}$. Additionally, as the field lines of ${\mathrm{\bf E}}$ are orthogonal to the equipotential lines of $\phi$, it follows that the set of equipotential lines $\phi= \mathrm{const}$ and $\psi = \mathrm{const}$ are orthogonal.

This derivation has shown that the duality of the field lines and equipotential lines within the given set of differential equations can be obtained by a replacement of Dirichlet (homogenous Neumann) by homogenous Neumann (Dirichlet) boundary conditions. The values of $\Psi_{2n-1}$ have to be evaluated by (5.28) and (5.43).

In this generalization of the capacitor problem, even more than two electrodes in the original constellation could be handled. In this case, if the original field problem has $n$ boundaries with imposed voltage, it requires also $n$ Neumann boundaries (electrodes for the dual system) which alternate with the Dirichlet boundaries along the whole surface.

As neither a bias of $\Psi_1$ nor an arbitrary constant scaling factor of $\psi$ change the shape of the field lines (equipotential lines of $\psi$), in case of exactly two contacts, the actual values of the two dual boundary voltages can be set to two different arbitrary values, most notably to values suitable for the computation.

5.3.6 The Three-Dimensional Capacitor Model

Unfortunately such a duality is not guaranteed in three dimensions. The theoretical expansion of the two-dimensional capacitor model delivers a three-dimensional capacitor with six sides where two of them act as electrodes. To determine the field lines, three potentials have to be calculated: one with the original electrode setting, one with two other opposite electrodes and the third one with the remaining two opposite electrodes.

One might expect that the expanded two-dimensional duality formulation can be expressed as follows:

The three fields build an orthogonal system. The field lines of each electric field lie within the equipotential surfaces of each of the remaining potentials and the conclusion is that the intersections of two equipotential surfaces deliver the field lines of the third potential.

Nevertheless, this might not hold under certain circumstances. In three dimensions, the duality formulation is not guaranteed, because the three potentials, which can be considered as a curvilinear coordinate system, usually do not represent an orthogonal system. Particularly, it is not even guaranteed that the field lines on the surface of the simulation domain are the equipotential lines of one of the other potentials. Additionally, the different electric fields are not perpendicular [70].

Figure 5.13: Spherical capacitor with electrodes placed at the top and bottom (hidden) cups. The electrodes for the second potential calculation are placed at the remaining cups at the left and right side. The electrodes for the third potential calculation are the remaining areas of the surface.

The orthogonality of the equipotential surfaces only exists, if the coordinate lines on the surface follow lines of curvature. This happens, when the coordinate lines follow surface lines of maximum or minimum curvature. However, this is only fulfilled for simple geometries and usually not satisfied [70].

The following example should clarify the situation. Consider a spherical device region as shown in Figure 5.13. The coordinates of the device are given in spherical coordinates by

$$\begin{pmatrix}x\ y\ z\end{pmatrix}=r \begin{pmatrix}\cos \alpha \cos \beta\ \sin \alpha \cos \beta\ \sin \beta\end{pmatrix}. \pagebreak[4]$$ (5.43)

Constant radius $r=r_0=\mathrm{const}$ delivers the surface of this domain. The meridians of the sphere are formed by $\alpha=\mathrm{const}$ and $\beta=\mathrm{const}$ and represent circles of constant height. The upper and lower contacts of the domain of the capacitor are cups of the sphere and will be limited by the circles at $\beta=\beta_0$ and $\beta=\beta_1$. The extracted potential is $\phi$. With this contact setting the solution is rotationally symmetric and it follows that the equipotential lines at the surface of the sphere are circles of constant height $\beta=\beta_i=\mathrm{const}$ and the field lines at the surface are meridians between $\beta_0$ and $\beta_1$ with $\alpha =\alpha_i=\mathrm{const}$.

Inside the device, the equipotential surfaces are perpendicular to the surface of the sphere and the field lines are also perpendicular to the former. Fortunately, a complete analytical solution of the potential distribution is not necessary. The major prediction is that the field lines at the surface of the sphere $\mathcal{C}$ are represented by meridians

$$\mathcal{C}=\mathcal{C}(\beta)\Bigl\vert _{\substack{r=r_0\ \alp... ...x}\cos \alpha_i \cos \beta\ \sin \alpha_i \cos \beta\ \sin \beta\end{pmatrix}$$ (5.44)

and the equipotential lines $\mathcal{P}$ at the surface

$$\mathcal{P}=\mathcal{P}(\alpha)\Bigl\vert _{\substack{r=r_0\ \be... ...\cos \alpha \cos \beta_i\ \sin \alpha \cos \beta_i\ \sin \beta_i\end{pmatrix}$$ (5.45)

are lines of constant height. At an intersection point, the tangential vectors of these curves are perpendicular which can be expressed by a vanishing inner product

$$\frac{\mathrm{d} \mathcal{C}}{\mathrm{d} \beta} \biggl\vert _{\su... ... \alpha_i \cos \beta_i\ \cos \alpha_i \cos \beta_i\ 0\end{pmatrix} \equiv 0 .$$ (5.46)

So far, there was no specification of the other four electrodes. However, any possible specification of them does not change the layout of the original field lines. The second potential $\psi$ is defined by the left and right electrodes as the cup of the remaining part of the sphere at $y<y_0$ and $y>y_1$. The front and back contacts are the remaining parts of the surface and deliver the potential $\vartheta$. The borderline of the left contact and the front contact is a line of constant y-coordinate

$$\mathcal{Q}=\mathcal{Q}(\beta)\Bigl\vert _{\substack{r=r_0\ y=y_... ... \beta\end{pmatrix}=\begin{pmatrix}f(\beta)\ y_0\ r_0 \sin \beta\end{pmatrix}$$ (5.47)

which is an equipotential line of $\psi$ as well as an equipotential line of $\vartheta$ and should be expected to be a field line of $\phi$. However, it can be seen that the curves $\mathcal{C}(r_0,\alpha_i,\beta)$ and $\mathcal{Q}(r_0,y_0,\beta)$ are not identical. The derivative

$$\frac{\mathrm{d}\mathcal{Q}}{\mathrm{d}\beta}(r_0,y_0,\beta)= \begin{pmatrix}f'(\beta)\ 0\ r_0 \cos \beta\end{pmatrix}$$ (5.48)

has no y-component and the inner product at an intersection point (remember $y_0=r_0 \sin \alpha_i \cos \beta_i$) usually does not give zero.

$$\begin{gather*}\begin{split}\frac{\mathrm{d} \mathcal{Q}}{\mathrm{d} \beta} \big... ...c{1}{2} \;r_0^2 \;\tan \alpha_i \sin 2 \beta_i \neq 0\ \end{split}\end{gather*}$$ (5.49)

Therefore, these lines are not orthogonal and cannot be field lines of $\phi$. In Figure 5.14(a) the equipotential lines and field lines of the original electrode placement are shown. Figure 5.14(b) shows the original equipotential lines and the equipotential lines of the switched electrode placement. These equipotential lines do not represent the field lines. However, the resulting quadrangles at the surface appear to be nearly orthogonal. Therefore, for three dimensions the field lines of $\phi$ cannot be formed by the intersections of the equipotential surfaces $\psi$ and $\vartheta$. Moreover, the equipotential surfaces of $\phi$, $\psi$, and $\vartheta$ are not orthogonal any longer.

Nevertheless, usually even these surfaces will be nearly orthogonal, within the geometries considered and due to the smoothness of the resulting potential lines arising from the Laplace equation. Therefore, nearly cuboidal elements and smooth grids will be produced by the equipotential surfaces. The equipotential surfaces conform to the boundaries. As a result all points inside the simulation area will be reached by the equipotential surfaces of $\phi= \mathrm{const}$, $\psi = \mathrm{const}$ and $\vartheta =\mathrm{const}$ -- a geometrical point $(x,y,z)$ inside the simulation domain can be described by a triple $(\phi,\psi,\vartheta )$. Only in rare cases, an ambiguous mapping, a folding of the curvilinear coordinate lines, may arise [36][70].

Figure 5.14: Potential and field lines of the spherical capacitor.

$\parbox{8.8cm}{\includegraphics[width=8.8cm] {picsconveps/mehrere}}$
Field lines and equipotential lines defined by the original electrode setting

        

$\parbox{8.8cm}{\includegraphics[width=8.8cm] {picsconveps/nureine}}$
Equipotential lines defined by the replaced electrodes. Additionally the equipotential lines of the original electrode placement are drawn.

5.3.7 Evaluation of the Electric Potential of the Capacitor

By selecting several potential value ticks of the three different potentials and calculating the intersection points of the equipotential surfaces, the point set for the grid generation is derived. With the selection of the tick values of the potentials, the grid density can be controlled in various ways.

To simplify the expressions, the boundary values of the potentials are set to 0 and $1$. Each geometrical point $(x,y,z)$ inside the capacitor has a potential representation $(\phi,\psi,\vartheta )$ within the potential range $([0,1]\times[0,1]\times[0,1])$. As the discrete maximum principle is satisfied, it seems plausible that the potential values raise monotonously along the field lines, going from one electrode to the opposite one. Therefore, except for folding it is a one-one mapping of geometrical points to potential coordinates. In this way, a selected potential triple $(\phi,\psi,\vartheta )$ delivers a distinct geometrical point $(x,y,z)$ inside the capacitor.

In a first step it is necessary to find an initial Delaunay tetrahedrization of the capacitor domain. This grid generation step does not need a special grid refinement. All material properties are linear, in particular a constant permittivity is used, and the Laplace equation as an elliptic partial differential equation is numerically well behaved. In addition, as a big advantage different from generating a grid for the whole device, only one capacitor segment has to be meshed. This grid can be produced regardless of other segments. This means that no additional grid points are induced by geometry and grid constrains of neighboring segments. A relatively crude grid can be used. Only for preserving the quality of the resulting cuboids, it is useful to introduce a refinement criterion based on the directions of the electric fields.

The solutions of the field equations is usually obtained by the Box Integration method (refer Section 3.1). The differential equation systems for solving the three potential distributions can be seen in Figure 5.1. The resulting equation system is of the form

$${\mathrm{\bf A}} \cdot \makebox{\boldmath$\phi$}= {\mathrm{\bf b}}.$$ (5.50)

with

$$\makebox{\boldmath$\phi$}=(\phi_i)$$ (5.51)

the unknown potential values on the grid points $p_i$.

The elements of the system matrix ${\mathrm{\bf A}}=(a_{ij})$ and the right-hand side ${\mathrm{\bf b}}=(b_i)$ are set to

$${a_{ij}} {=-\varepsilon _{ij}\;\frac{A_{ij}}{d_{ij}}}$$ (5.52)

$${a_{ii}} {=\sum\limits_{\forall\; i\neq j} a_{ij}}$$ (5.53)

$${b_i} {=0}$$ (5.54)

$\forall\; i,j$ : $\exists\; < p_i p_j >$ and $p_i$ lying inside $\partial V$ or on $\partial V_3$, $\partial V_4$, $\partial V_5$ or on $\partial V_6$,

$${a_{ii}} {=1}$$ (5.55)

$${b_i} {=\Phi_1}$$ (5.56)

$\forall\; i,j$ : $\exists\; < p_i p_j >$ and $p_i$ lying on $\partial V_1$,

$${a_{ii}} {=1}$$ (5.57)

$${b_i} {=\Phi_2}$$ (5.58)

$\forall\; i,j$ : $\exists\; < p_i p_j >$ and $p_i$ lying on $\partial V_2$ and

$${a_{ij}} {=0}$$ (5.59)

for all other $i,j$.

These equations are used for calculating for the potential values $\phi_i$ on the grid points $p_i$. For the potentials $\psi$ and $\vartheta$ an similar equation system is set up. Basically, all three potential evaluations use the same system matrix ${\mathrm{\bf A}}$, except the lines and rows resulting from the boundary points where the adequate boundary conditions have to be inserted. For the evaluation of the potential values, a Conjugate Gradient Solver is advisable [71].

5.3.7.1 Why Choosing a Conjugate Gradient Solver

The matrix $\bf A$, for the equation system ${\mathrm{\bf A}} \cdot \makebox{\boldmath $\phi$}={\mathrm{\bf b}}$, is a sparse $n \times n$ matrix. Applying a Gaussian Solver destroys the sparse system, because by eliminating the entries under the diagonal a lot of entries above the diagonal become non-zero. Therefore, the memory usage for storing the system matrix elements increases to $O\left( n^2\right)$. Additionally the elimination of a fully occupied row requires $O\left( n^2\right)$ arithmetic operations and eliminating all rows is of $O\left(n^3\right)$. Therefore the use of a solver, which preserves the sparsity of the system and does not show a complexity of $O\left(n^3\right)$ is a better choice. With simple modifications, the system matrix ${\mathrm{\bf A}}$ can be made suitable for a Conjugate Gradient solver. Since a CG-solver requires a symmetrical system matrix, the original matrix must be transformed to become symmetric. The symmetry of the system matrix is violated in rows and columns where the values of $\phi_i$ are imposed. One possible method is to pre-eliminate these rows and columns and move the adequate entries to the right-hand side. This method requires a row number to point number remapping and each potential system delivers a completely different system matrix.

Another way is to eliminate asymmetrical matrix entries $a_{ji}$ on the left-hand side by

$$a_{jk}' = a_{jk} - a_{ik}\;\frac{a_{ji}}{a_{ii}}\qquad \forall\; k,$$ (5.60)

$$b_{j}' = b_{j} - b_{i}\;\frac{a_{ji}}{a_{ii}},$$ (5.61)

$\forall\; i,j$ with $i<j, a_{ji}\neq a_{ij}$ and $p_i$ on $\partial V_1$ or on $\partial V_2$.

As all off-diagonal elements $a_{ik}$ with $i\neq k$ of the concerned rows are zero, the transformation does not change any other entries than $a_{ji}$ and simplifies to

$$a_{ji}' = a_{ji} - a_{ii}\;\frac{a_{ji}}{a_{ii}},$$ (5.62)

$$b_{j}' = b_{j} - b_{i}\;\frac{a_{ji}}{a_{ii}},$$ (5.63)

$\forall\; i,j$ with $i<j$ and $p_i$ on $\partial V_1$ or on $\partial V_2$.

For an efficient memory representation of the equations, it is necessary to use a sparse matrix representation. The chosen method is only practicable, if the matrix values are left unchanged during the whole solving process and direct access to the row entries can be omitted. This representation can be assembled easily by iteration over all points and insertion of the couplings of the edges that are connected with the iteration point. The storage representation of a matrix works with four arrays diag[], offdiag[], colidx[], startcolidx[], where the entries are defined as [49]:

Algorithm
5.3.1   Sparse representation of ${\mathrm{\bf A}}=\left( a_{ij} \right)$

$a_{ii}=$ diag[i]

startcolidx[i] <= k <= startcolidx[i+1]

$j=$ colidx[k]

$a_{ij}=$ offdiag[k]

With this method, matrix-vector products can be performed directly on this data structure. Direct access to the row indices can be omitted by using the following algorithm:

Algorithm
5.3.2   Matrix-vector product ${\mathrm{\bf A}}\cdot {\mathrm{\bf x}}={\mathrm{\bf c}}$

for j in all point indices

c[j]=diag[j] $\cdot$ x[j]

for k=startcolidx[j] to startcolidx[j+1]-1

c[j]=c[j] + offdiag[k] $\cdot$ x[colidx[k]]

No special sorting of the column indices is necessary. Matrix-matrix products, however, require direct access to the row entries within the above matrix representation and thus a search algorithm within the row indices becomes necessary. This kind of products has to be avoided. Therefore, a simple CG-algorithm without pivoting and without preconditioning is used. The algorithm looks as follows:

Algorithm
5.3.3   CG-Method to solve ${\bf A}\cdot {\bf x}={\bf b}$

1. set

   $${\bf r}_0={\bf b}-{\bf A} \cdot {\bf x}_0, \qquad {\bf d}_0={\bf r}_0$$

   
   

2. $k=0,1,...$ until r is sufficiently set

   $$\begin{gather*}\begin{split}\lambda_k&=\frac{{\bf r}^T_k \cdot {\bf r}_k}{{\bf d... ...\bf d}_{k+1}&={\bf r}_{k+1}+\beta_k {\bf d}_k.\nonumber \end{split}\end{gather*}$$

   
   

In this case, only vector-vector products ( $O\left(n\right)$) and matrix-vector products are needed. If only the non-zero entries of the matrix $\bf A$ are stored, the memory usage is of the same order as the number of existing edges $m$. The number of arithmetic operations for a matrix-vector product is of $O\left(m\right)$, too. With the use of exact arithmetic, the equation solver delivers the solution at least after $n$ recursion steps. With non-exact arithmetic a border for the residuum ${\mathrm{\bf r}}$ must be given and usually this border is reached with less than $n$ steps [20][63][71].

While the mean value of edges in a tetrahedral mesh is of the same order than the number of grid points (usually multiplied with a typical factor, but surely not of quadratic or higher order), the matrix-vector product is also of $O\left(n\right)$. Therefore the memory usage is of $O\left(n\right)$ and the amount of arithmetic operations for solving the equation system is of order $O\left( n^2\right)$, which is a huge improvement compared to a Gaussian solver.