4.1 Two-Dimensional Scalar Finite Element Method

In this section the principle of the finite element method is demonstrated by solving the Poisson equation (4.11) in the two-dimensional domain $\mathcal{A}$ . One part of the domain boundary $\partial\mathcal{A}$ represents the Dirichlet boundary $\mathcal{C}_D$ and the remaining part -- the Neumann boundary $\mathcal{C}_N$ , respectively ( $\partial\mathcal{A}= \mathcal{C}_D + \mathcal{C}_N$ ).

As discussed in Chapter 3 the unknown function $\varphi$ is approximated by $\tilde{\varphi}$ as in (3.5). The known basis functions will be notated with $\lambda_j$ instead of $N_j$ . The reason for this will became clear in the vector finite element chapter, where systems of partial differential equations appear, which require a combination of vector and scalar interpolation functions.

$$\varphi(\vec{r}) \simeq \tilde{\varphi}(\vec{r}) = \sum_{j=1}^nc_j\lambda_j(\vec{r}) + v(\vec{r}).$$ (4.13)

The function $v(\vec{r})$ complies with the Dirichlet boundary function on $\mathcal{C}_D$ . The set of known functions $\lambda_i, i\in[1,n]$ builds a fundamental function system vanishing on the Dirichlet boundary $\mathcal{C}_D$ . To find the numerically approximated solution of (4.11) for $\tilde{\varphi}$ it is necessary to determine the unknown multiplier coefficients $c_i$ . This is performed by substituting $\varphi$ in (4.11) by its approximation (4.13) and then weighting the resulting residual of (4.11) with the set of functions $\lambda_i$ in the domain $\mathcal{A}$

$$\int_{\mathcal{A}}[\vec{\nabla}\cdot(\utilde{\epsilon}\cdot\vec{\... ...\mathrm{d}A = \int_{\mathcal{A}}f(\vec{r}) \lambda_i \mathrm{d}A, i\in[1;n].$$ (4.14)

Equation (4.14), using the first scalar theorem of Green (3.15) analogously to (3.16), leads to

$$\int_{\partial\mathcal{A}}\lambda_i \vec{n}\cdot\utilde{\epsilon... ...\mathrm{d}A = \int_{\mathcal{A}}f(\vec{r}) \lambda_i \mathrm{d}A, i\in[1;n],$$ (4.15)

with the fact that all functions $\lambda_i (i\in[1;n])$ vanish on $\mathcal{C}_D$ to

$$\int_{\mathcal{C}_N}\lambda_i \vec{n}\cdot\utilde{\epsilon}\cdot... ...\mathrm{d}A = \int_{\mathcal{A}}f(\vec{r}) \lambda_i \mathrm{d}A, i\in[1;n].$$ (4.16)

The Neumann boundary condition on $\mathcal{C}_N$ can be written analogously to (3.18) as

$$-D_n(\vec{r}) = -\vec{n}\cdot\vec{D}(\vec{r}) = \vec{n}\cdot\utilde{\epsilon}\cdot\vec{\nabla}\tilde{\varphi}(\vec{r})$$ (4.17)

using $\vec{\nabla}\tilde{\varphi} = -\vec{E}$ and $\vec{D} = \utilde{\epsilon}\cdot\vec{E}$ . After substitution of $\tilde{\varphi}$ in (4.16) by the approximation (4.13) and using the Neumann boundary condition (4.17) the following linear equation system for the unknown coefficients $c_i (i\in[1;n])$ is obtained

$$\int_{\mathcal{A}}\vec{\nabla}\lambda_i\cdot\utilde{\epsilon}\cdo... ...bda_i \mathrm{d}A - \int_{\mathcal{C}_N}\lambda_iD_n \mathrm{d}s, i\in[1;n].$$ (4.18)

The above equation can be converted to

$$\begin{split}\sum_{j=1}^nc_j\int_{\mathcal{A}}\vec{\nabla}\la... ...mathcal{C}_N}\lambda_iD_n \mathrm{d}s, i\in[1;n], \end{split}$$ (4.19)

which complies with the following linear equation system

$$[K]\{c\} = \{d\}.$$ (4.20)

The matrix $[K]$ and the right hand side vector $\{d\}$ are given by the expressions

$$\begin{split}K_{ij} & = \quad \int_{\mathcal{A}}\vec{\nabla}\... ...mathrm{d}s, & \quad \quad i\in[1;n], j\in[1;n]. \end{split}$$ (4.21)

Let assume that $\mathcal{C}_N$ is a part of the outer boundary of the two-dimensional domain $\mathcal{A}$ . Furthermore $\mathcal{A}$ is assumed sufficiently large to allow that $D_n$ can be neglected on $\mathcal{C}_N$ . Since the corresponding simulations are performed in finite domains, this assumption will lead to systematic error which became smaller with increasing domain size. This issue will be discussed in Subsection 4.1.5. It is also assumed that there is no electric charge density distribution $\rho$ in the domain of interest. Thus the Neumann boundary term and the source term in (4.21) are set to zero. Otherwise, if $\vec{D}$ is perpendicular to $\vec{n}$ on $\mathcal{C}_N$ , the boundary condition on $\mathcal{C}_N$ will be zero independently from the domain size.

The function $v(\vec{r})$ , which satisfies the Dirichlet boundary $\mathcal{C}_D$ can be analogously to the solution approximation (4.13) written as a sum of known functions multiplied with coefficients

$$v(\vec{r}) = \sum_{j=n+1}^mc_j\lambda_j.$$ (4.22)

Now the coefficients $c_j$ in (4.22) where $j\in[n{+}1;m]$ are known values. They are obtained easily from $v(\vec{r})$ for $m{-}n$ points on the Dirichlet boundary $\mathcal{C}_D$ by the following linear equation system

$$v(\vec{r}_i) = \sum_{j=n+1}^mc_j\lambda_j(\vec{r}_i), i\in[1;m{-}n], \vec{r}_i\in\mathcal{C}_D.$$ (4.23)

With (4.22) the right hand side vector $\{d\}$ from (4.21) can be rewritten as

$$\begin{split}d_i & = - \int_{\mathcal{A}}\vec{\nabla}\lambda_... ...hrm{d}A = - \sum_{j=n+1}^mK_{ij}c_j, i\in[1;n]. \end{split}$$ (4.24)

In general the form functions $\lambda_i$ can be defined in the entire simulation domain. For example, this is the case by the weighted residual method. In the case of the finite element method the domain $\mathcal{A}$ is divided into smaller sub-domains $\mathcal{A}_e$ . This process is known as domain discretization and the resulting sub-domains $\mathcal{A}_e$ are called elements. The index $^e$ or $_e$ will be used in this work generally for quantities, which comply with an element. It can be considered as a numbering index. The shape functions $\lambda_i$ are non-zero only in a few neighboring elements and vanish in the remaining simulation area. $\lambda_i$ is a global shape function. Its local representation in each element, in which it is non-zero is termed as element shape function $\lambda_i^e$ . Now $^e$ is written as a superscript index, because there is already another numbering index available -- the subscript index $_i$ . In each element usually many element shape functions are defined, which are part of different global shape functions. As it will be shown later for this purpose low order polynomials are used. It is very important to notice the difference in the indexing (index $i$ ) between the global and the element basis functions. In the global basis function $\lambda_i$ $i$ complies with the number of global functions in the whole region and with the number of unknown coefficients $c_i$ , respectively. By the element basis function $\lambda_i^e$ the index $i$ corresponds to the number of the element basis functions in the element. In practice at first the element form functions are defined for each element in the domain and then the global ones are constructed by them.

In this work the area of interest is discretized on triangular elements. On these elements linear triangular element functions are employed.