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

geometries,

boundary conditions,

boundaries between domains with different material behavior, and

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 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

for Without loss of generality it can be assumed that , which easily can be achieved due to the commutativity of the differentiation order.
Definition 1. The differential equation 4.1 is elliptic in , iff all eigenvalues of the matrix have the same sign and are not equal to zero.
Definition 2. The differential equation 4.1 is elliptic (in ), iff it is elliptic in (almost) all

For the further description 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 and are used for the nabla and the Laplace operator, respectively.

##### 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 and then integrate the result over the domain giving [17]

For FEM based simulations only gradients of the solution and of the test function 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 is a vector field and is an outward pointing normal vector of the surface surrounding the volume , then the volume integral of the divergence of can be converted to the flux integral through the surface by

Applying Theorem 1 to a vector field defined by leads to Green’s first identity [70].
Theorem 2. Green’s first identity

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

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 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. ) and quadrilaterals (cf. Figure 4. ) for two-dimensional (2D) domains and tetrahedrons (cf. Figure 4. ) and hexahedrons (cf. Figure 4. ) in 3D simulation domains [6]. Figure 4.2 shows an example triangular mesh for an 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.

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

over the reference element by employing the determinant 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]

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

##### 4.1.3  Ritz-Galerkin-Approximation

Let be (real) Hilbert space and a continuous bilinear form and bounded linear functional on . To solve a problem, given by [17]

approximately, the infinite space is substituted by the finite space . The Ritz-Galerkin-Approximation is defined as the solution of

This problem is finite-dimensional and with being a complete basis of the solution can be written by the linear combination

(4.9) has to be valid for all and to solve the problem it is enough to use only the basis functions transforming the problem to

and by substituting by (4.10) a system of linear equations

with unknowns is obtained. By defining the matrix and the vectors and the system can be written in matrix notation

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

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

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

and a general coordinate by

resulting in the hat function

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

and

The weak formulation (4.6) can now be transformed in a summation of integrals over each of the triangles , leading to

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

for mesh points it follows

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

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

The elements of matrix are defined by

and of matrix by

and , and is represented by the vectors , and , respectively. The are joined into the matrix These are sparse matrices as only those integrals in (4.26), where the and are edges of the triangle that is integrated over are not zero.

Previous: 3.6 Void Evolution    Top: Home    Next: 4.2 Implementation in FEDOS