Finding the analytical solution for partial differential equations in three dimensions is in general only for quite simple problems regarding
boundaries between domains with different material behavior, and
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.
The finite element method is a common approach to solve PDEs . 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 . 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.
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 
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 .
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).
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 . 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 
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 .
Let be (real) Hilbert space and a continuous bilinear form and bounded linear functional on . To solve a problem, given by 
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 . 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
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.« PreviousUpNext »Contents