« Previous Up Next »Contents
Previous: 3.6 Void Evolution Top: Home Next: 4.2 Implementation in FEDOS

4 Numerical Implementation

Finding the analytical solution for partial differential equations in three dimensions is in general only for quite simple problems regarding

$\bullet$ geometries,

$\bullet$ boundary conditions,

$\bullet$ boundaries between domains with different material behavior, and

$\bullet$ material behavior
possible. General TCAD simulations are to complicated and numerical methods have to be used. Numerical methods produce approximations of the solutions. One of the methods mainly used and widespread for Partial Differential Equations (PDEs) in engineering is the Finite Element Method (FEM). Due to the availability of powerful computers since the mid $20^{\mathrm {t}\mathrm {h}}$ century this method got practicable. Since this method was used for the simulations in this work, the basics will be described in this, followed by a discussion of the implementation of the EM model in the two software packages used.

4.1 Finite Element Method

The finite element method is a common approach to solve PDEs [67]. Thereby the simulation region is filled with a mesh and the solution is approximated by a linear combination of weighted base functions, which is called the Ritz-Galerkin-approximation [17]. This approach enables the software to approximate the solution by solving a system of linear equations representable by matrix equations.

A second order PDE can generally be written as

$(4.1) \{begin}{gather} \displaystyle \sum _{i,j=1}^{d}\frac {\partial }{\partial x_{i}}a_{ij}(\mathrm {\b {x}})\frac {\partial u}{\partial x_{j}}+\sum _{i}^{d}a_{i}(\mathrm {\b {x}})\frac {\partial u}{\partial x_{i}}+a(\mathrm {\b {x}})u=f(\mathrm {\b {x}}),\qquad \mathrm {\b {x}}\in \Omega \mathref {(4.1)} \{end}{gather}$

for $\Omega \subset \mathbb {R}^{d}$ Without loss of generality it can be assumed that $a_{ij} = a_{ji}$ , which easily can be achieved due to the commutativity of the differentiation order.
Definition 1. The differential equation 4.1 is elliptic in $\mathrm {\b {x}}\in \Omega$ , iff all eigenvalues of the matrix $(\mathrm {a}_{ij})$ have the same sign and are not equal to zero.
Definition 2. The differential equation 4.1 is elliptic (in $\Omega$ ), iff it is elliptic in (almost) all $\mathrm {\b {x}}\in \Omega .$

For the further description $d$ will be set to three, as the simulations are carried out in the three-dimensional real space and the time derivatives are realized by time stepping to reduce the problem order by solving it iteratively approximated by a finite differences schema. Furthermore, the material parameters in EM are rotationally symmetric leading to the following simplified form of the differential equations, where the symbols $\nabla$ and $\triangle$ are used for the nabla and the Laplace operator, respectively.

$(4.2) \{begin}{gather} \nabla \cdot (\alpha (\mathrm {\b {x}})\nabla u(\mathrm {\b {x}}))+\pmb \beta (\mathrm {\b {x}})\ \cdot \nabla u(\mathrm {\b {x}})+\gamma (\mathrm {\b {x}})u(\mathrm {\b {x}})=b(\mathrm {\b {x}})\qquad \mathrm {\b {x}}\in \Omega \mathref {(4.2)} \{end}{gather}$

4.1.1 Weak Formulation of PDEs

In order to derive the weak formulation of (4.2) one first has to multiply (4.2) by a test function $\phi (\mathrm {\b {x}})$ and then integrate the result over the domain $\Omega$ giving [17]

$(4.3) \{begin}{gather} \displaystyle \int _{\Omega }\ (\nabla \ (\alpha (\mathrm {\b {x}})\nabla u(\mathrm {\b {x}}))+\pmb {\beta }(\mathrm {\b {x}})\ \cdot \nabla u(\mathrm {\b {x}})+\gamma (\mathrm {\b {x}})u(\mathrm {\b {x}})-b(\mathrm {\b {x}}))\phi (\mathrm {\b {x}})\mathrm {d}V=0. \mathref {(4.3)} \{end}{gather}$

For FEM based simulations only gradients of the solution $u(\mathrm {\b {x}})$ and of the test function $\phi (\mathrm {\b {x}})$ occur in the integral, allowing the usage of solution functions with lower smoothness. This makes the employment of Green’s first identity necessary.
Theorem 1. Divergence theorem: If $\b {F}$ is a vector field and $\b {n}$ is an outward pointing normal vector of the surface $\partial \Omega$ surrounding the volume $\Omega$ , then the volume integral of the divergence of $\b {F}$ can be converted to the flux integral through the surface by

$(4.4) \{begin}{gather} \displaystyle \int _{\Omega }\nabla \cdot \b {F} dV_{n}=\oint _{\partial \Omega }\b {F}\cdot \mathbf {n}dS_{n-1}. \mathref {(4.4)} \{end}{gather}$

Applying Theorem 1 to a vector field defined by $\b {F} = \alpha \psi \nabla \phi$ leads to Green’s first identity [70].
Theorem 2. Green’s first identity

$(4.5) \{begin}{gather} \displaystyle \int _{\Omega } (\psi \nabla \cdot (\alpha \nabla \varphi )+\alpha \nabla \varphi \cdot \nabla \psi ) dV=\displaystyle \oint _{\partial \Omega }\alpha \psi (\nabla \varphi \cdot \mathbf {n}) dS=\displaystyle \oint _{\partial \Omega }\alpha \psi \nabla \varphi \cdot d\mathbf {S} \mathref {(4.5)} \{end}{gather}$

Employing Theorem 2 to (4.3) results in the so called weak formulation of the problem, which is given by

$(4.6) \{begin}{gather} \displaystyle \oint _{\partial \Omega }\alpha (\mathrm {\b {x}})\phi (\mathrm {\b {x}})\nabla u(\mathrm {\b {x}})\cdot d\mathbf {S}+\notag \\ \displaystyle \int _{\Omega } (-\alpha (\mathrm {\b {x}})\nabla u(\mathrm {\b {x}})\cdot \nabla \phi (\mathrm {\b {x}})+\pmb \beta (\mathrm {\b {x}})\cdot \nabla u(\mathrm {\b {x}})\phi (\mathrm {\b {x}})+\gamma (\mathrm {\b {x}})u(\mathrm {\b {x}})\phi (\mathrm {\b {x}})-b(\mathrm {\b {x}})\phi (\mathrm {\b {x}}))\mathrm {d}V=0. \mathref {(4.6)} \{end}{gather}$

where the first part is a surface integral and the second integral is carried out over the whole volume. For the further discussion it is assumed that $a(\mathrm {\b {x}})$ is a material parameter being piecewise constant and therefore not depending on the position. For example in the drift-diffusion equation of the vacancies it is proportional to the diffusion coefficient (c.f. Section 3.2).

4.1.2 Meshing of the Simulation Domains

Before the actual solving of the PDE can be started, the simulation domain has to be meshed. The meshing procedure subdivides the simulation domain into many small elements, described by nodes, edges, and faces. The generated mesh can either be structured by simple subdivision rules or unstructured for arbitrary domains. Typical shapes used for FEM are triangles (cf. Figure 4. $1\mathrm {a}$ ) and quadrilaterals (cf. Figure 4. $1\mathrm {b}$ ) for two-dimensional (2D) domains and tetrahedrons (cf. Figure 4. $1\mathrm {c}$ ) and hexahedrons (cf. Figure 4. $1\mathrm {d}$ ) in 3D simulation domains [6]. Figure 4.2 shows an example triangular mesh for an $\mathrm {L}$ shaped 2D domain. The mesh is refined at locations of higher interest and smaller geometrical features. In this example this is the fillet, as the boundary is a quarter of a circle.

Figure 4.1.: Example meshing elements for 2D and 3D domains.

Based on the meshing and its element type a simple reference element can be chosen establishing a mapping. For triangular meshing this mapping is shown in Figure 4.3. Due to this mapping the integration of every element can be formulated as an integration

Figure 4.2.: An example mesh for an $\mathrm {L}$ shaped structure with a fillet.

$ \mathrm {L} $ — Figure 4.2.: An example mesh for an $\mathrm {L}$ shaped structure with a fillet.

over the reference element by employing the determinant $\mathcal {J}(\mathrm {\b {x}}(s,t))$ of the Jacobian matrix. If the mapping can be written as a linear combination of the vectors spanning the element, where the origin of the element is one of the corner coordinates, the integral melts down to the volume of the original element times an integral not depending on the original element giving the possibility to precompute the integrals of (4.6). For triangular meshing elements (cf. Figure 4.3) the transformation of the integral over an element and the reference element is given by [68]

$(4.7) \{begin}{gather} \displaystyle \int \limits _{T}f(\mathrm {\b {x}})\mathrm {d}V=\int \limits _{\substack {0\leq s+t\leq 1\\0\leq s,t\leq 1}} f(\mathrm {\b {x}}(s,t))\mathcal {J}(\mathrm {\b {x}}(s,t)) \mathrm {d}s\mathrm {d}t =\notag \\ \displaystyle \int \limits _{\substack {0\leq s+t\leq 1\\0\leq s,t\leq 1}} f(\mathrm {\b {x}}(s,t))|(\mathrm {\b {x}}_{1}-\mathrm {\b {x}}_{0},\mathrm {\b {x}}_{1}-\mathrm {\b {x}}_{0})| \mathrm {d}s\mathrm {d}t =2\displaystyle \mathrm {Vol}(T)\int \limits _{\substack {0\leq s+t\leq 1\\0\leq s,t\leq 1}} f(\b {x})\mathrm {d}s\mathrm {d}t. \mathref {(4.7)} \{end}{gather}$

Important for the transformation is the preservation of the orientation. Choosing for this meshing a test function space $V_{h}$ and, therefore, a solution space leads to the Ritz-Galerkin-approximation [17].

Figure 4.3.: Mapping between an arbitrary triangle and the reference triangle for 2D triangular meshing with the transformation law.

4.1.3 Ritz-Galerkin-Approximation

Let $V$ be $\mathrm {a}$ (real) Hilbert space and $a:V\times V\rightarrow \mathbb {R}$ a continuous bilinear form and $F\mathrm {a}$ bounded linear functional on $V$ . To solve a problem, given by [17]

$(4.8) \{begin}{gather} a(u,\ v)=F (v)\qquad \forall v\in V, \mathref {(4.8)} \{end}{gather}$

approximately, the infinite space $V$ is substituted by the finite space $V_{h} \subset V$ . The Ritz-Galerkin-Approximation $u_{h}\in V_{h}$ is defined as the solution of

$(4.9) \{begin}{gather} a(u_{h},\ v)=F (v)\qquad \forall v\in V_{h}. \mathref {(4.9)} \{end}{gather}$

This problem is finite-dimensional and with $\{\phi _{1},\ \phi _{N}\}$ being a complete basis of $V_{h}$ the solution can be written by the linear combination

$(4.10) \{begin}{gather} u_{h}=\displaystyle \sum _{i=1}^{N}u_{h,i}\phi _{i}. \mathref {(4.10)} \{end}{gather}$

(4.9) has to be valid for all $v\in V_{h}$ and to solve the problem it is enough to use only the basis functions $\phi _{i}$ transforming the problem to

$(4.11) \{begin}{gather} a(u_{h},\ \phi _{i})=F(\phi _{i}),\qquad i=1,\ N \mathref {(4.11)} \{end}{gather}$

and by substituting $u_{h}$ by (4.10) a system of $N$ linear equations

$(4.12) \{begin}{gather} \displaystyle \sum _{j=1}^{N}a(\phi _{j},\ \phi _{i})u_{h,j}=F(\phi _{i}),\qquad i=1,\ \cdots ,\ N \mathref {(4.12)} \{end}{gather}$

with $N$ unknowns is obtained. By defining the matrix $\mathbf {A}= (A_{ij}) = (a(\phi _{j},\ \phi _{i}))$ and the vectors $\mathbf {b}= (b_{i}) =(F(\phi _{i}))$ and $\mathbf {u}=(u_{j})$ the system can be written in matrix notation

$(4.13) \{begin}{gather} \mathbf {A u}=\mathbf {b}. \mathref {(4.13)} \{end}{gather}$

and reduces the implementation to the calculation of the solution of a system of linear equations.

Figure 4.4.: The hat test function for a single node located in the middle.

A typical test function used for the weak formulation (4.6) is the so called hat function shown in Figure 4.4 [6]. It is defined to be 1 at the meshing point considered and decreases linearly in direction of the neighboring meshing points to zero. For a single triangle the function is defined by

$(4.14) \{begin}{gather} f(\displaystyle \mathrm {\b {x}})=\frac {|(\mathrm {\b {x}}_{j}-\mathrm {\b {x}})-\left ((\mathrm {\b {x}}_{j}-\mathrm {\b {x}})\cdot \frac {\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{k}}{|\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{k}|}\right )\frac {\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{k}}{|\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{k}|}|}{|(\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{i})-\left ((\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{i})\cdot \frac {\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{k}}{|\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{k}|}\right )\frac {\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{k}}{|\mathrm {\b {x}}_{j}-\mathrm {\b {x}}_{k}|}|}. \mathref {(4.14)} \{end}{gather}$

For the reference triangle the coordinates of the points are given by

$(4.15) \{begin}{gather} \mathrm {\b {x}}_{i}=(0,0)^{T}. \mathref {(4.15)} \{end}{gather}$

$(4.16) \{begin}{gather} \mathrm {\b {x}}_{j}=\ (0,1)^{T}. \mathref {(4.16)} \{end}{gather}$

$(4.17) \{begin}{gather} \mathrm {\b {x}}_{k}=(1,0)^{T}. \mathref {(4.17)} \{end}{gather}$

$(4.18) \{begin}{gather} \mathref {(4.18)} \{end}{gather}$

and a general coordinate by

$(4.19) \{begin}{gather} \mathrm {\b {x}}=(s,t)^{T} \mathref {(4.19)} \{end}{gather}$

resulting in the hat function

$(4.20) \{begin}{gather} f_{i,(s,t)}(\mathrm {\b {x}})=1-s-t. \mathref {(4.20)} \{end}{gather}$

The same procedure can be applied to the other two meshing points of this triangle leading to the functions

$(4.21) \{begin}{gather} f_{j,(s,t)}(\mathrm {\b {x}})=s \mathref {(4.21)} \{end}{gather}$

and

$(4.22) \{begin}{gather} f_{k,(s,t)}(\mathrm {\b {x}})=t. \mathref {(4.22)} \{end}{gather}$

The weak formulation (4.6) can now be transformed in a summation of integrals over each of the $N$ triangles $T_{l}$ , leading to

$(4.23) \{begin}{gather} \displaystyle \sum _{l=1_{\partial T_{l}}}^{N}\oint _{\cap \partial \Omega } \alpha \phi (\mathrm {\b {x}})\nabla u(\mathrm {\b {x}}) \cdot d\mathbf {S}+\notag \\ \displaystyle \sum _{l=1_{T}}^{N}\int _{l} (-\alpha \nabla u(\mathrm {\b {x}})\cdot \nabla \phi (\mathrm {\b {x}})+\pmb {\beta }(\mathrm {\b {x}})\ \nabla u(\mathrm {\b {x}})\phi (\mathrm {\b {x}})+\gamma (\mathrm {\b {x}})u(\mathrm {\b {x}})\phi (\mathrm {\b {x}})-b(\mathrm {\b {x}})\phi (\mathrm {\b {x}}))\mathrm {d}V=0 \mathref {(4.23)} \{end}{gather}$

and by applying the approximation of all functions as a linear combination of the hat functions given by

$(4.24) \{begin}{gather} g(\displaystyle \mathrm {\b {x}})=\sum _{i=1}^{M}g_{i}f_{i}(\mathrm {\b {x}}) \mathref {(4.24)} \{end}{gather}$

for $M$ mesh points it follows

$(4.25) \{begin}{gather} \displaystyle \sum _{i=1}^{M}\sum _{j=1}^{M}u_{j}\phi _{i}\sum _{l=1}^{N} (\displaystyle \oint _{\partial T_{l}\cap \partial \Omega } \alpha f_{i}(\mathrm {\b {x}})\nabla f_{j}(\mathrm {\b {x}})\cdot d\mathbf {S}-\notag \\ \alpha \int _{T_{l}}\nabla f_{j}(\mathrm {\b {x}})\cdot \nabla f_{i}(\mathrm {\b {x}})\mathrm {d}V+\sum _{k=1}^{M}\pmb {\beta }_{k}\int _{T_{l}}f_{k}(\mathrm {\b {x}})\cdot \nabla f_{j}(\mathrm {\b {x}})f_{i}(\mathrm {\b {x}})\mathrm {d}V+\notag \\ \displaystyle \sum _{k=1}^{M}\gamma _{k}\int _{T_{l}}f_{k}(\mathrm {\b {x}})f_{j}(\mathrm {\b {x}})f_{i}(\mathrm {\b {x}})\mathrm {d}V)\ -\sum _{k=1}^{M}\sum _{l=1}^{N}b_{k}\int _{T_{l}}f_{k}(\mathrm {\b {x}})f_{i}(\mathrm {\b {x}})\mathrm {d}V=0. \mathref {(4.25)} \{end}{gather}$

As this result has to be valid for every function $\phi$ it is enough to choose only one hat function at a time.

$(4.26) \{begin}{gather} \displaystyle \sum _{j=1}^{M}u_{j}(\sum _{l=1}^{N}2\mathrm {V}\mathrm {o}\mathrm {l}(T_{l})(\oint _{\partial T_{l}\cap \partial \Omega } \alpha f_{i}(\mathrm {\b {x}})\nabla f_{j}(\mathrm {\b {x}}) \cdot d\mathbf {S}’-\notag \\ \alpha \int _{T_{l}}\nabla f_{j}(\mathrm {\b {x}})\cdot \nabla f_{i}(\mathrm {\b {x}})\mathrm {d}V’+\pmb {\beta }_{k}\int _{T_{l}}f_{k}(\mathrm {\b {x}})\ \nabla f_{j}(\mathrm {\b {x}})f_{i}(\mathrm {\b {x}})\mathrm {d}V’+\notag \\ \displaystyle \gamma _{k}\int _{T_{l}}f_{k}(\mathrm {\b {x}})f_{j}(\mathrm {\b {x}})f_{i}(\mathrm {\b {x}})\mathrm {d}V’)) =\displaystyle \sum _{l=1}^{N}b_{k}\int _{T_{l}}f_{k}(\mathrm {\b {x}})f_{i}(\mathrm {\b {x}})\mathrm {d}V’\quad \forall i\in \{1,\ \cdot \ \cdot \ \cdot \ ,\ M\} \mathref {(4.26)} \{end}{gather}$

The primed coordinate system is the coordinate system of the reference system and allows the calculation of every integral in advance. This leads to the matrix representation

$(4.27) \{begin}{gather} \mathbf {A}(\pmb \beta ,\pmb \gamma )\mathbf {u}=\mathbf {C}\mathbf {b}. \mathref {(4.27)} \{end}{gather}$

The elements of matrix $\b {A}$ are defined by

$(4.28) \{begin}{gather} A_{ij}=\displaystyle \sum _{l=1}^{N}2\mathrm {V}\mathrm {o}\mathrm {l}(T_{l})\left (\oint _{\partial T_{l}\cap \partial \Omega } \displaystyle \alpha f_{i}(\mathrm {\b {x}})\nabla f_{j}(\mathrm {\b {x}})\cdot d\mathrm {S}’-\alpha \int _{T_{l}}\nabla f_{j}(\mathrm {\b {x}})\cdot \nabla f_{i}(\mathrm {\b {x}})\mathrm {d}V’+\right .\notag \\ \displaystyle \left .\pmb \beta _{k}\int _{T_{l}}f_{k}(\mathrm {\b {x}})\cdot \nabla f_{j}(\mathrm {\b {x}})f_{i}(\mathrm {\b {x}})\mathrm {d}V’+\gamma _{k}\int _{T_{l}}f_{k}(\mathrm {\b {x}})f_{j}(\mathrm {\b {x}})f_{i}(\mathrm {\b {x}})\mathrm {d}V’\right ) \mathref {(4.28)} \{end}{gather}$

and of matrix $\mathbf {C}$ by

$(4.29) \{begin}{gather} C_{ik}=\displaystyle \sum _{l=1}^{N}b_{k}\int _{T_{l}}f_{k}(\mathrm {\b {x}})f_{i}(\mathrm {\b {x}})\mathrm {d}V. \mathref {(4.29)} \{end}{gather}$

and $b_{k}, u_{j}$ , and $\gamma _{k}$ is represented by the vectors $\mathbf {b}, \mathbf {u}$ , and $\gamma$ , respectively. The $\pmb \beta _{k}$ are joined into the matrix $\pmb {\beta }’$ These are sparse matrices as only those integrals in (4.26), where the $i, j$ and $k$ are edges of the triangle that is integrated over are not zero.

« Previous Up Next »Contents
Previous: 3.6 Void Evolution Top: Home Next: 4.2 Implementation in FEDOS