Next: 6. Simulation of Thermal Up: 5. Discretization with the Previous: 5.2 Discretization with Tetrahedrons

Subsections

# 5.3 Assembling and Solving

Regarding assembling the finite element method is based on the principle that the components of the local element matrices must be assembled to a global matrix for building up a global equation system. Only with the global system all unknown variables on the grid nodes in the discretized domain can be determined. Each global grid node is shared by a varible number of finite elements which all make a contribution to the solution of the unknown values on the involved nodes.

## 5.3.1 Principle of Assembling

The assembling procedure from a local element matrix to a global matrix has the same routines for two- and three-dimensional structures. Therefore, the assembling procedure is demonstrated on a simpler two-dimensional example as shown in Fig. 5.3. The dimension of the local matrix is always , where is the number of grid nodes on the finite element ( for triangles and for tetrahedrons) and is the number of unknown variables on a grid node. The dimension of the global matrix is always , where is the total number of grid nodes in the discretized domain.
If it is assumed that there is only one sought variable

 (5.92)

the dimension of is in the two-dimensional case and the dimension of is .

The local matrix uses the local node indexes which are 1, 2, and 3 for every finite element. The global indexes for these grid nodes are different. There must exists a transformation which projects the local indexes and of the components to the global indexes and . For example, the element 66 in Fig. 5.3 with its local nodes 1, 2, and 3 has the global nodes 42, 30, and 25 and the index transformation is

 (5.93)

This means that the components of the local matrix from element 66 are transformed to the global matrix in the way

 (5.94)

Another important aspect is that a global grid node is shared by a number of different finite elements. This fact is taken into account during assembling of the global matrix by the so-called superposition principle. This means that the components of the global matrix are summed up from the contributions of the local element matrices. For example, the global grid node 30 is shared by the five elements 63, 64, 65, 66, and 67 (see Fig. 5.3), which all make a contribution to the node 30. So the global matrix components are found with the help of index transformations by the way

 (5.95)

The assembling of the global matrix from all elements with the index transformation can be described in the form

 (5.96)

Sometimes there are more than one variable on the grid nodes. For example, with two variables the size of the local matrix (two dimensions) is and for the global matrix . If the variables are independent, the offset of the entries for the second variable is 3 in and in . For assembling the second variable from the element 66 the index transformation (5.94) must be modified by adding the respective offset

 (5.97)

## 5.3.2 Dirichlet Boundary Conditions

Through the Dirichlet boundary conditions the values on the surface grid nodes are already fixed with the so-called Dirichlet value. Therefore, it is not necessary and even not allowed to recalculate the values on these grid nodes from the global equation system , because it is impossible to obtain the same Dirichlet values by solving the equation system. These surface grid nodes must be treated differently with the Dirichlet value. If on the global node there is a Dirichlet boundary value , the global equation system must be changed to

 (5.98)

From the mathematical point of view the global equation system has pseudo-equations if there are grid nodes with Dirichlet conditions, after setting all rows and columns from Dirichlet grid nodes in to 0.

In practice it is more comfortable to multiplicate with a transformation matrix

 (5.99)

which sets all rows and columns for the Dirichlet grid nodes in to 0, instead of doing it componentwise by and for . In the beginning is a unit matrix ( and ), but for every Dirichlet grid node the components are reset to . Therefore, all rows and columns in can easily be set to 0 at once with .

## 5.3.3 Mechanical Interfaces

The simulated structures generally consist of several segments. Normally a segment is a continuous region of one material and so there are everywhere the same material characteristics. Therefore, the electrical and mechanical behavior can be described with the same parameters within a segment.

The used meshing module makes a separated grid for every segment. This means that for a global grid point at the interface between two different segments there exist two different indices, because each segment has its own global index, as shown in Fig. 5.4.

The two global indices for the same grid point lead to an important aspect for the mechanics with regard to finite elements, because in the global stiffness matrix there exist entries for two node indices for the same grid point. After solving the linear equation system for the mechanical problem there would be two different, but wrong displacement vectors for the same grid point.

The calculated displacements would be wrong, because the entries from segment A for this interface points in the global stiffness matrix do not take into account that there is also a stiffness from the other segment B due to the separate assembling of the segments. The same but reverse explanation is valid for the entries from segment B. Therefore, the stiffness in matrix and also the force in must be corrected in the way that the entries with second index are added to the first index of the grid point and then all entries with second index are set to 0.

In the following the correcting procedure in is demonstrated for the representative interface point with Index 35 in Segment A and Index 85 in Segment B (see Fig. 5.4).

After the addition of all entries with Index 85 these components are set to zero: and for . The other interface points are handled in the same way. It is also necessary to rearrange the force vector with the same concept so that

 (5.101)

Like in (5.99) and can be manipulated for all interface points in a faster and simpler way with the same transformation matrix , because the Dirichlet boundary conditions have higher priority and must always be fulfilled. With the mechanical system becomes

 (5.102)

For the representative interface point with Index 35/85 the matrix elements from must be set to and in order to get the required effect.

After solving the linear system for displacements (5.102), has the correct value and . In order to also have the correct displacement for the second Index 85 (segment B), must be set to .

## 5.3.4 Complete Equation System for Oxidation

The oxidation problem induces 5 unknown variables: the oxidation concentration , the normalized silicon concentration , and the three components , , and (in x-, y-, and z-direction) of the displacement vector . If the dopant redistribution is taken into account, than the unknown variables for the concentrations of the different species have also to be considered. In case the five-stream diffusion model these are the 5 variables for , , , , and (see Section 4.2).

The following explanation is focused on the pure oxidation problem. Therefore, the discretization of the simulation domain with grid nodes leads to 5 variables on each grid node. After assembling of the whole equation system there are in total 5 unknown variables and, as already described in Section 5.3.1, the dimension of the global matrix is 55.

As shown in Fig. 5.5 the global matrix consists of the two sub-matrices and . contains the coupled entries from the oxidant diffusion (5.65) and the -dynamics (5.74). Because of the unknown variables and the size of is 22. is the global stiffness matrix for the mechanics with the three displacement components , , and , and so its dimension is 33.

The structure of the local system matrix for one finite element is the same as for the global one (see Fig. 5.5). The only difference is that because of the 4 nodes of a tetrahedron. So the dimension of the local matrix is 2020 with the 88 and 1212 sub-matrices.

The first part of the equation system which describes the oxidation diffusion and the -dynamics, is non-linear because of the coupling between and in the form (see (5.65) and (5.74)). The variables in the non-linear system can not be calculated directly, only their increments and . The second part of the system which is responsible for the mechanical problem is linear, as marked in Fig. 5.5.

For a quasi-stationary time step there is no coupling between the --system and the mechanical system, because the equations for and are not functions of displacements. Although the normalized additional volume after a time step is calculated with and , the mechanical equations for the displacements also are not functions of and . Because of the missing coupling the off-diagonal sub-matrices in the global equation system are zero (see Fig. 5.5).

On the right-hand side of the non-linear subsystem are the residuals calculated with the results from the last Newton iteration as explained in the next section. The right-hand side of the linear subsystem contains the (internal) forces on the grid nodes.

## 5.3.5 Solving with the Newton Method

In contrast to the linear mechanical subsystem, where the displacements can be calculated directly with one of the various standard methods like Gaussian elimination, the solving of the non-linear part demands other routines like the used Newton Method [96].

The global non-linear subsystem can be written in the form

 (5.103)

where stand for the variables and stands for the variables .

With the Newton formula the solution vector for the actual time step becomes

 (5.104)

where is the solution from the previous time step , is the inverse Jacobian matrix, and are the residuals, both determined by .

Transforming the Newton formula to the form

 (5.105)

leads to a linear equation system for the increments (see Fig. 5.5).
After solving this linear system the values for the actual timestep can be determined by

 (5.106)

This equation shows that the Newton method demands start values on all grid nodes for the first iteration .

Because the Newton method is only a first order approximation of the solution, an iteration never provides the exact results. Therefore the right-hand side of the non-linear system (5.103) can not be 0, instead there always exist residuals .

The quality of the approximation is increased with each iteration, but the number of iterations must be limited with termination conditions. With these conditions also the accuracy of the approximation can be controlled. It makes sense to use the following termination conditions [97]

 (5.107) (5.108) (5.109)

where , , and are given tolerances. is the Euclidean norm

 (5.110)

(5.107) and (5.108) are failure criterions, which only work, if the sequence converges quickly to the exact solution. If this converging sequence is slow, can be small but the demanded accuracy is by far not reached and the Newton loop is terminated too early. Due to this fact it is recommended to use the additional residual criterion (5.109). The combination of both criterions ensures good terminating conditions.

For one finite tetrahedral element with its 4 nodes () the non-linear system (5.103) is

 (5.111)

Here come from the oxidant diffusion (5.65)

 (5.112)

and the other equations describe the -dynamics (5.74)

 (5.113)

The local Jacobian matrix is

 (5.114)

For calculating the Jacobian matrix and the residuals in the Newton iteration the results for and from the previous iteration are used.
For the equations (5.112) of the Jacobian matrix components are

The partial derivatives of the equations (5.113) result in

The local residuals used for the actual Newton iteration are

The entries of the global Jacobian matrix and global residual vector as needed for (5.105) are assembled from the local ones as explained in Section 5.3.1.

Next: 6. Simulation of Thermal Up: 5. Discretization with the Previous: 5.2 Discretization with Tetrahedrons

Ch. Hollauer: Modeling of Thermal Oxidation and Stress Effects