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.
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.
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).
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.
(a) Triangle
(c) Tetrahedron
(b) Quadrilateral
(d) Hexahedron
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 shaped structure with a fillet.
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].
Figure 4.3.: Mapping between an arbitrary triangle and the reference triangle for 2D triangular meshing with the transformation law.
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.
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
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.