Subsections

• 5.3.1 Principle of Assembling
• 5.3.2 Dirichlet Boundary Conditions
• 5.3.3 Mechanical Interfaces
• 5.3.4 Complete Equation System for Oxidation
• 5.3.5 Solving with the Newton Method

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 $\mathbf{A^e}$ to a global matrix $\mathbf{A^g}$ 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 $\mathbf{A^e}$ is always $n k \times n k$, where $n$ is the number of grid nodes on the finite element ($n=3$ for triangles and $n=4$ for tetrahedrons) and $k$ is the number of unknown variables on a grid node. The dimension of the global matrix $\mathbf{A^g}$ is always $N k\times N k$, where $N$ is the total number of grid nodes in the discretized domain.
If it is assumed that there is only one sought variable $\varphi$

$$\mathbf{A^e} \varphi^e = b^e, \qquad \mathbf{A^g} \varphi^g = b^g$$ (5.92)

the dimension of $\mathbf{A^e}$ is $3\times 3$ in the two-dimensional case and the dimension of $\mathbf{A^g}$ is $N\times N$.

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 $T(\mathbf{A^e})$ which projects the local indexes $k$ and $l$ of the components $A^e_{kl}$ to the global indexes $i$ and $j$. 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

$$1 \to 42, \qquad 2 \to 30, \qquad \mathrm{and} \qquad 3 \to 25.$$ (5.93)

This means that the components of the local matrix $\mathbf{A^{66}}$ from element 66 are transformed to the global matrix $\mathbf{A^g}$ in the way

$$A^{66}_{1,1} \to A^{g}_{42,42}, \qquad A^{66}_{1,2} \to A^{g}_{42... ...66}_{1,3} \to A^{g}_{42,25}, \qquad A^{66}_{2,1} \to A^{g}_{30,42}, \quad\ldots$$ (5.94)

Figure 5.3: Part of a mesh with finite elements and grid nodes.

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 $A^g_{ij}$ of the global matrix are summed up from the contributions $A^e_{kl}$ 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 $A^g_{30,j}$ are found with the help of index transformations by the way

$$A^g_{30,25}= A^{65}_{2,1}+ A^{66}_{2,3}, \qquad A^g_{30,29}= A^{63}_{3,2}+ A^{64}_{1,2}, \quad\ldots$$ (5.95)

The assembling of the global matrix from all $N$ elements with the index transformation $T(\mathbf{A^e})$ can be described in the form

$$\mathbf{A^g}= \sum_{e=1}^N T(\mathbf{A^e}).$$ (5.96)

Sometimes there are more than one variable on the grid nodes. For example, with two variables the size of the local matrix $\mathbf{A^e}$ (two dimensions) is $6\times 6$ and $2N\times 2N$ for the global matrix $\mathbf{A^g}$. If the variables are independent, the offset of the entries for the second variable is 3 in $\mathbf{A^e}$ and $N$ in $\mathbf{A^g}$. For assembling the second variable from the element 66 the index transformation (5.94) must be modified by adding the respective offset

$$A^{66}_{3+1,3+1} \to A^{g}_{N+42,N+42}, \quad A^{66}_{3+1,3+2} \to A^{g}_{N+42,N+30}, \quad A^{66}_{3+1,3+3} \to A^{g}_{N+42,N+25}, \ldots$$ (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 $\mathbf{A^g} \varphi^g = b^g$, 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 $i$ there is a Dirichlet boundary value $\varphi_i = C_i$, the global equation system must be changed to

$$\left[ \begin{array}{ccccccccc} a_{1,1} &a_{1,2} &\cdots &a_{1,i-... ... b_2 \vdots b_{i-1} C_{i} b_{i+1} \vdots b_{N} \end{array} \right]$$ (5.98)

From the mathematical point of view the global equation system $\mathbf{A^g} \varphi^g = b^g$ has $m$ pseudo-equations if there are $m$ grid nodes with Dirichlet conditions, after setting all $m$ rows and columns from Dirichlet grid nodes $\varphi_i$ in $\mathbf{A^g}$ to 0.

In practice it is more comfortable to multiplicate $\mathbf{A^g}$ with a transformation matrix $\mathbf{T_b}$

$$\mathbf{T_b} \mathbf{A^g} \varphi^g = b^g,$$ (5.99)

which sets all rows and columns for the $m$ Dirichlet grid nodes $\varphi_i = C_i$ in $\mathbf{A^g}$ to 0, instead of doing it componentwise by $A^g_{ik}=0$ and $A^g_{ki}=0$ for $k=1\ldots N$. In the beginning $\mathbf{T_b}$ is a unit matrix ($T_{ii}=1$ and $T_{ij}=0$), but for every Dirichlet grid node $i$ the components $T_{ii}$ are reset to $T_{ii}=0$. Therefore, all $m$ rows and columns in $\mathbf{A^g}$ can easily be set to 0 at once with $\mathbf{T_b}$.

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.

Figure 5.4: Two segments with interface and its grid points.

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 $\mathbf{K^g}$ there exist entries for two node indices for the same grid point. After solving the linear equation system $\mathbf{K^g} \vec{d^g}=\vec{f^g}$ 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 $\mathbf{K^g}$ and also the force in $\vec{f^g}$ 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 $\mathbf{K^g}$ is demonstrated for the representative interface point with Index 35 in Segment A and Index 85 in Segment B (see Fig. 5.4).

$$\mathbf{K^g}= \left[ \begin{array}{ccccccccccccc} k_{1,1} &\cdots... ...s &k_{N,34} &k_{N,35}+k_{N,85} &k_{N,36} &\cdots &k_{N,N} \end{array} \right]$$

$$\mathrm{with}\qquad k_{SUM}=k_{35,35}+k_{35,85}+k_{85,35}+k_{85,85}$$ (5.100)

After the addition of all entries with Index 85 these components are set to zero: $k_{i,85}=0$ and $k_{85,i}=0$ for $i=1\ldots N$. The other interface points are handled in the same way. It is also necessary to rearrange the force vector $\vec{f^g}$ with the same concept so that

$$\vec{f^g}= \big[ f_{1} \;\;\;f_{2} \ldots f_{34} \;\;\;f_{35}... ...85} \;\;\;f_{36} \ldots f_{84} \;\;\;0 \;\;\;f_{86} \ldots f_{N} \big]^T$$ (5.101)

Like in (5.99) $\mathbf{K^g}$ and $\vec{f^g}$ can be manipulated for all interface points in a faster and simpler way with the same transformation matrix $\mathbf{T_b}$, because the Dirichlet boundary conditions have higher priority and must always be fulfilled. With $\mathbf{T_b}$ the mechanical system becomes

$$\mathbf{T_b} \mathbf{K^g} \vec{d^g}=\mathbf{T_b} \vec{f^g}.$$ (5.102)

For the representative interface point with Index 35/85 the matrix elements from $\mathbf{T_b}$ must be set to $T_{35,85}=1$ and $T_{85,85}=0$ in order to get the required effect.

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

5.3.4 Complete Equation System for Oxidation

The oxidation problem induces 5 unknown variables: the oxidation concentration $c$, the normalized silicon concentration $\eta$, and the three components $u$, $v$, and $w$ (in x-, y-, and z-direction) of the displacement vector $\vec{d}$. 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 $C_{A^+}$, $C_I$, $C_V$, $C_{AI}$, and $C_{AV}$ (see Section 4.2).

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

As shown in Fig. 5.5 the global matrix $\mathbf{A^g}$ consists of the two sub-matrices $\mathbf{J^g}$ and $\mathbf{K^g}$. $\mathbf{J^g}$ contains the coupled entries from the oxidant diffusion (5.65) and the $\eta$-dynamics (5.74). Because of the unknown variables $c_1\ldots c_N$ and $\eta_1\ldots \eta_N$ the size of $\mathbf{J^g}$ is 2$N\times$2$N$. $\mathbf{K^g}$ is the global stiffness matrix for the mechanics with the three displacement components $u_1\ldots u_N$, $v_1\ldots v_N$, and $w_1\ldots w_N$, and so its dimension is 3$N\times$3$N$.

Figure 5.5: Structure of the global equation system for oxidation.

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 $N=4$ because of the 4 nodes of a tetrahedron. So the dimension of the local matrix $\mathbf{A^e}$ is 20$\times$20 with the 8$\times$8 $\mathbf{J^e}$ and 12$\times$12 $\mathbf{K^e}$ sub-matrices.

The first part of the equation system which describes the oxidation diffusion and the $\eta$-dynamics, is non-linear because of the coupling between $c_i$ and $\eta_i$ in the form $c_i \eta_i$ (see (5.65) and (5.74)). The variables in the non-linear system can not be calculated directly, only their increments $\Delta c_i$ and $\Delta \eta_i$. 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 $c$-$\eta$-system and the mechanical system, because the equations for $c$ and $\eta$ are not functions of displacements. Although the normalized additional volume after a time step is calculated with $c$ and $\eta$, the mechanical equations for the displacements also are not functions of $c$ and $\eta$. 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

$$\begin{array}{c} f_1(x_1, x_2,\ldots, x_{2N})=0 f_2(x_1, x_... ...)=0 \vdots f_{2N}(x_1, x_2,\ldots, x_{2N})=0 \end{array}$$ (5.103)

where $x_1\ldots x_N$ stand for the variables $c_1\ldots c_N$ and $x_{N+1}\ldots x_{2N}$ stands for the variables $\eta_1\ldots \eta_N$.

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

$$\vec{x}^n=\vec{x}^{n-1}-\mathbf{J^g}^{-1}(x^{n-1}) \vec{R}(x^{n-1}),$$ (5.104)

where $\vec{x}^{n-1}$ is the solution from the previous time step $n-1$, $\mathbf{J^g}^{-1}$ is the inverse Jacobian matrix, and $\vec{R}$ are the residuals, both determined by $x^{n-1}_i$.

Transforming the Newton formula to the form

$$\mathbf{J^g}\vec{\Delta x}=\vec{R}\qquad\mathrm{with}\qquad \Delta x_i=x^n_i-x^{n-1}_i,$$ (5.105)

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

$$x^n_i =x^{n-1}_i+ \Delta x_i$$ (5.106)

This equation shows that the Newton method demands start values $x^{0}_i$ on all $N$ grid nodes for the first iteration $n=1$.

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 ${R}_i(x^{n-1}_i)$.

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]

$$\Vert\vec{x}^n-\vec{x}^{n-1}\Vert\leq\tau_{abs},$$ (5.107)

$$\Vert\vec{x}^n-\vec{x}^{n-1}\Vert\leq\tau_{rel}\Vert\vec{x}^n\Vert,$$ (5.108)

$$\Vert\vec{R}(x^n)\Vert\leq\tau_{f},$$ (5.109)

where $\tau_{abs}$, $\tau_{rel}$, and $\tau_{f}$ are given tolerances. $\Vert \Vert$ is the Euclidean norm

$$\Vert\vec{x}^n-\vec{x}^{n-1}\Vert=\sqrt{\sum\big(x_i^n-x_i^{n-1}\big)^2}\qquad\mathrm{and}\qquad \Vert\vec{x}^n\Vert=\sqrt{\sum \big(x_i^n\big)^2}.$$ (5.110)

(5.107) and (5.108) are failure criterions, which only work, if the sequence $\vec{x}^n$ converges quickly to the exact solution. If this converging sequence is slow, $\Vert\vec{x}^n-\vec{x}^{n-1}\Vert$ 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 ($N=4$) the non-linear system (5.103) is

$$\begin{array}{c} f_1(c_1^n, c_2^n, c_3^n, c_4^n, \eta_1^n,\et... ..., c_4^n, \eta_1^n,\eta_2^n, \eta_3^n,\eta_4^n)=0. \end{array}$$ (5.111)

Here $f_1\ldots f_4$ come from the oxidant diffusion (5.65)

$$f_j= \sum_{i=1}^4\big(D K_{ij} c_{i}^{n}+ k_{max} M_{ij} c_{i}^{n} \eta_{i}^{n}\big)=0 \quad\mathrm{for}\quad j=1,2,3,4, \vspace*{-0.1cm}$$ (5.112)

and the other equations $f_5\ldots f_8$ describe the $\eta$-dynamics (5.74)

$$f_{j+4}= \sum_{i=1}^4\big(M_{ij}\big(\eta_{i}^{n} -\eta_{i}^{{n-1... ..._{i}^{n} \Delta t\big)=0,\quad \mathrm{for} \quad j=1,2,3,4. \vspace*{-0.1cm}$$ (5.113)

The local Jacobian matrix is

$$\mathbf{J^e}= \left[ \begin{array}{cccccccc} \displaystyle\frac{\... ...n} &\displaystyle\frac{\partial f_8} {\partial \eta_4^n} \end{array} \right]$$ (5.114)

For calculating the Jacobian matrix and the residuals in the Newton iteration $n$ the results for $c_i$ and $\eta_i$ from the previous iteration $n-1$ are used.
For the equations $f_1\ldots f_4$ (5.112) of the Jacobian matrix components are

$$\phantom{xxxxxxxxxxxxxxxxxxxxxxx}J_{i,j} = \frac{\partial f_i} {\partial c_j^n}= D K_{ij}+ k_{max} M_{ij} \eta_{i}^{n-1},$$

$$J_{i,j+4} = \frac{\partial f_i} {\partial \eta_j^n}= k_{max} M_{ij} c_{i}^{n-1} \quad\qquad\mathrm{for}\quad i,j=1,2,3,4.$$ (5.115)

The partial derivatives of the equations $f_5\ldots f_8$ (5.113) result in

$$\phantom{xxxxxxxxxxxxxxxxxxxxxxx}J_{i+4,j} = \frac{\partial f_{i+4}}{\partial c_j^n}= K_A M_{ij} \eta_{i}^{n-1} \Delta t,$$

$$J_{i+4,j+4} = \frac{\partial f_{i+4}}{\partial \eta_j^n}= M_{ij}+ K_A M_{ij} c_{i}^{n-1} \qquad\mathrm{for}\quad i,j=1,2,3,4.$$ (5.116)

The local residuals used for the actual Newton iteration $n$ are

$$\phantom{xxxxxxxxxxxxxxxxxxxxxx}R_i= \sum_{j=1}^4\big(D K_{ij} c_{i}^{n-1}+k_{max} M_{ij} c_{i}^{n-1} \eta_{i}^{n-1}\big),$$

$$R_{i+4}= \sum_{j=1}^4\big(M_{ij}\big(\eta_{i}^{n-1} -\eta_{i}^{{n-2}}\big)... ...,\eta_{i}^{n-1} c_{i}^{n-1} \Delta t\big)\qquad \mathrm{for} \quad i=1,2,3,4.$$ (5.117)

The entries of the global Jacobian matrix $\mathbf{J^g}$ and global residual vector $\vec{R}$ as needed for (5.105) are assembled from the local ones as explained in Section 5.3.1.