Next: 3.2 The Diffusion Equation Up: 3. The Box Integration Previous: 3. The Box Integration

Subsections

# 3.1 The Poisson Equation

The Poisson equation is fundamental for all electrical applications. The derivation is shown for a stationary electric field [26]. For the derivation, the material parameters may be inhomogeneous, locally dependent but not a function of the electric field. In this section, the principle of the discretization is demonstrated. In the case of a low field dependency, a recursive reinsertion of newly derived material parameters may result in stable solutions [54]. A highly nonlinear behavior has to be managed by linearization and Newton iterations, for instance [2][4][60].

The Maxwell equations for a stationary electric field reduce to

 (3.1) (3.2)

and are related by the permittivity

 (3.3)

Physical properties require that the permittivity tensor is symmetric and positive definite for lossless media (or a positive scalar for isotropic media). Equation (3.1) is satisfied implicitly by introducing the electric potential

Combing (3.1)-(3.4), the Maxwell equations can be rewritten as

This can be transformed to its integral formulation. By applying the Gaussian integral formula (3.5) results in

Equation (3.6) must be fulfilled within any arbitrary volume , with being the surface of this volume. While performing Box Integration, this formula must be satisfied in the Voronoi boxes of each grid point. In Figure 3.1 a typical Voronoi box of a point , belonging to all the tetrahedrons connected to the point , can be seen. The integration can be split into the integration over the different box parts, caused by all the tetrahedrons which share the same point

 (3.7)

Such a part of the Voronoi box within the tetrahedron is shown in Figure 3.2. Each part around the point can also be split into the three area contributions that are caused by the different edges , and (see Figure 3.3)

 (3.8)

Using Box Integration, an approximation of first order is usually assumed, which means that the electric potential within one tetrahedron is a linear function along the axes or, in other words, the electric field is constant within each tetrahedron. Material parameters are assumed to be constant within the tetrahedron, too. The electric charge density on the right-hand-side is assumed to be constant within the whole Voronoi box as well and can therefore be analytically calculated and written as the electric charge of the Voronoi box of the point , which results in

As in the usual applications the permittivity tensor reduces to a scalar, which means the electric displacement and the electric field have the same direction, it follows

 (3.10)

and equation (3.9) can be simplified to

## 3.1.1 Approximation of and

Usually the integral (3.11) is approximated by

with

 and (3.13)

An estimation of the validity of this approximation based on the homogenous Poisson equation (Laplace equation) will be explained under the conditions shown in Figure 3.4. spans a local Cartesian coordinate system with origin . In the one-dimensional case and with constant permittivity the potential varies linearly with position as

 (3.14)

Whereas with non-constant permittivity we have at first order

 (3.15)

with

 and (3.16)

and therefore,

 (3.17)

The continuity of the electric displacement

 (3.18)

with

 (3.19)

delivers the deformed distribution of

 (3.20)

and therefore,

 (3.21)

The resulting potential distributions between and for different values of are shown in Figure 3.5. The overestimation of the electric field in the middle of the grid element in dependence of the permittivity change inside the element is shown in Figure 3.6. According to this figure, a limit for the validity of approximation (3.13) can be defined by

 edge (3.22)

or a limit for continuously defined permittivity

 edge (3.23)

With , the estimation error of is limited to 10%.

On a three-dimensional grid (refer to Figure 3.4), the spatial derivatives of cause an additional contribution to the discretization error. With linear interpolation of and , equation (3.11) results in

 (3.24)

and with

The first addend of equation (3.25) equals the value given by discretization (3.12) and (3.13) whereas the remaining terms cause an additional discretization error. Often the derivative components of neighboring tetrahedrons compensate each other, but not always. Even the lengths of and are not limited to lie within the tetrahedron ranges (for instance, obtuse angled triangles for two dimensions) and a limitation for this approximation is difficult. For the whole simulation domain the estimation

 edge (3.26)

seems plausible, but is too pessimistic as in most cases and have the same order of magnitude as their based tetrahedrons. The previously derived estimation (3.23) is therefore also applicable to the three-dimensional case.

Another common approximation method for the effective permittivity is [27]

 (3.27)

which shows the behavior of a serial connection of two capacitors (or resistors for electrical flow fields). This approximation has the advantage that if the permittivity (more likely for the conductivity) is zero in a Voronoi box, no flux enters the box. Nevertheless, this zero flux leaves the electric potential of the included grid point undefined and may cause problems for the equation solver, applied later on.

## 3.1.2 Continuation of the Discretization

However, equation (3.11) is discretized as

In compliance with the decomposition of the Voronoi regions in Section 2.4.3, is the right signed value of the Voronoi area between and of the tetrahedron . Which means in detail, area portions with
• outer-circle inside the appropriate triangle and outer-sphere inside the tetrahedron counted positive,
• outer-circle outside the appropriate triangle and outer-sphere outside the tetrahedron counted positive, and
• one of them inside and the other one outside counted negative.
Combining the portions of the tetrahedrons which share the point delivers the commonly used notation

with the Voronoi area of the edge (coupling area between and )

 (3.30)

Equation (3.29) must be fulfilled for each grid point inside the simulation domain, which results in an equation system

with the system matrix

 (3.32) (3.33)

the unknowns

 (3.34)

and the right-hand side

which is defined for all inner points of the domain.

## 3.1.3 Boundary Conditions

For solving this time-invariant elliptical partial differential equation, boundary conditions have to be applied to all boundaries. For each grid point lying at the boundary of the domain one of the two following additional equations or a linear-combination of them must be added for the discrete version, which fulfill the boundary conditions.

Each grid point of the boundary has to be defined by exactly one of the different conditions. Additionally, at least one grid point must be described by a Dirichlet condition to define a definite potential distribution. With this set of equations, the discrete Poisson equation can be solved.

### 3.1.3.1 Dirichlet Conditions

The first class of boundary conditions is an imposed potential

 on boundary . (3.36)

The discrete formulation of this boundary condition simply replaces the respective rows of equation (3.29) by

 (3.37)

or within the matrix representation

 (3.38) (3.39) (3.40)

which at first delivers an asymmetric system matrix . Since most equation solvers require a symmetric system matrix, it may be necessary to transform the matrix to a symmetric one, which can be simple obtained by summation of adequate rows.

### 3.1.3.2 Requirements for an M-Matrix

As shown in Section 2.3, for satisfying the maximum principle the system matrix has to be an M-matrix. Under the precondition of a valid Delaunay grid, the validity of the M-matrix will be shown.

As the grid is Delaunay, it follows that the coupling areas are positive

and therefore,

 (3.32) satisfies (2.4), (3.33) satisfies (2.5), (3.33) = satisfies (2.6),
and finally by defining at least one boundary point as a Dirichlet point
 (3.38)-(3.39) > satisfies (2.7).

### 3.1.3.3 Neumann Conditions

The second type of boundary condition specifies the imposed normal component on the boundaries (see Figure 3.7), which occurs between different segments. As shown in Figure 3.7, the normal vector may have various directions on different boundary pieces around . This behavior could be handled by splitting the boxes (and also the point itself) into several pieces. Usually this approach is prevented and the normal derivatives are combined to one value. In this case, the boundary conditions result in a distortion of (3.29)

with the assigned normal component at the boundary at the grid point . It must be distinguished between which represents the area of the Voronoi box between and , and as the assigned boundary area of the Voronoi box around the grid point .

In the majority of cases, a special kind of this boundary condition is used, which shows a vanishing flux across the boundary (homogenous Neumann condition)

 (3.42)

This type of condition specifies the artificial boundaries which terminate the simulation domains. As no distortion of equation (3.41) appears, these boundaries do not require a special treatment and the equations can be left unchanged as in (3.31)-(3.35).

### 3.1.3.4 Cauchy Boundary Conditions

The linear-combination

 (3.43)

combines the two previous conditions and results in the discrete formulation

Next: 3.2 The Diffusion Equation Up: 3. The Box Integration Previous: 3. The Box Integration

J. Cervenka: Three-Dimensional Mesh Generation for Device and Process Simulation