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
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
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
In practice it is more comfortable to multiplicate
with a transformation matrix
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
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 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
.
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 5
5
.
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 2
2
.
is the global stiffness matrix for the mechanics with the three displacement components
,
, and
, and so its dimension is 3
3
.
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 20
20 with the 8
8
and 12
12
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.
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
With the Newton formula the solution vector for the actual time step becomes
Transforming the Newton formula to the form
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]
For one finite tetrahedral element with its 4 nodes () the non-linear system (5.103) is