Subsections

• 5.2.1 Shape Functions for a Tetrahedron
• 5.2.2 Coordinate Transformation
• 5.2.3 Differentiation in the Normalized Coordinate System
• 5.2.4 Discretization of the Oxidant Diffusion
• 5.2.5 Discretization of the ${\eta }$-Dynamics
• 5.2.6 Discretization of the Mechanics

5.2 Discretization with Tetrahedrons

An often used choice is to discretize the solution domain with tetrahedrons. On the one side a tetrahedron is a relative simple element, especially regarding meshing aspects, on the other side it is an efficient element to discretize structures with non-planar surfaces or complex geometries. As shown in Fig. 5.1, this element is limited by four triangles and has four vertexes which are in any case grid nodes in the mesh.

5.2.1 Shape Functions for a Tetrahedron

For using the weighted residual method the shape functions must be continuous on the transition from one element to its neighbor element. Within the elements they must be at least one-time differentiable. The shape functions will be defined locally on the tetrahedron. It should be noted that the global shape function $N_j(x,y,z)$ is assembled from the local shape functions of the elements which share the same node $j$.

If it is assumed that the discretization is carried out with linear shape functions, the four vertexes used are the four grid nodes on the element. This means that the shape functions must depend on the x-, y-, and z-coordinate linearly, so that in general form it can be expressed by [90]

$$N_i(x,y,z)=a_i + b_i x + c_i y + d_i z \qquad \mathrm{for} \qquad i=0,1,2,3.$$ (5.8)

$i$ are the numbers of the local grid node. Since for every grid node a separate shape function is needed, there are four shape functions on a four node element. The coefficients $a_i$, $b_i$, $c_i$, and $d_i$ must be determined in such a manner, that the respective shape function $N_i$ fulfills (5.4). This means, that for example the value of the shape function $N_0$(x,y,z) must be 1 in node P$_0$ with its coordinates ($x_0$,$y_0$,$z_0$), and 0 in all other nodes. With this information the following equation system with N$_0$ on the four nodes can be written:

$$a_0 + b_0 x_0 + c_0 y_0 + d_0 z_0=1$$

$$a_0 + b_0 x_1 + c_0 y_1 + d_0 z_1=0$$

$$a_0 + b_0 x_2 + c_0 y_2 + d_0 z_2=0$$ (5.9)

$$a_0 + b_0 x_3 + c_0 y_3 + d_0 z_3=0 \nonumber$$

With this equation system it is possible to determine the four unknown coefficients $a_0$, $b_0$, $c_0$, and $d_0$ for the shape function $N_0$(x,y,z).

Figure 5.1: Tetrahedral element in a global $(x,y,z)$-coordinate system

Figure 5.2: Tetrahedral element in a normalized $(\xi ,\eta ,\zeta )$-coordinate system

For calculation of the unknown coefficients Cramer's rule can be applied, which says: if there is an equation system $\mathbf{A}\cdot\vec{x}=\vec{b}$, the numbers

$$x_i=\frac{1}{\mathrm{Det}(\mathbf{A})} \left\vert \begin{array}{c... ... &\cdots &A_{n i-1} &b_n &A_{n i+1} &\cdots &A_{nn} \end{array} \right\vert$$ (5.10)

are the components of the solution $\vec{x}$.

With Cramer's rule, for example, the coefficient $a_0$ from the shape function $N_0(x,y,z)$ can be calculated by

$$a_0= \frac{1}{\mathrm{Det}(\mathbf{D})} \left\vert \begin{array}{... ...1 &y_1 &z_1 0& x_2 &y_2 &z_2 0& x_3 &y_3 &z_3 \end{array} \right\vert,$$ (5.11)

where the matrix $\mathbf{D}$ is

$$\mathbf{D}= \left[ \begin{array}{cccc} 1& x_0 &y_0 &z_0 1& x_1 &y_1 &z_1 1& x_2 &y_2 &z_2 1& x_3 &y_3 &z_3 \end{array} \right].$$ (5.12)

By deleting the row i and the column j of a n-rowed determinant, a new (n-1)-rowed sub-determinant $\alpha_{ij}$ with sign $(-1^{i+j})$ is constructed

$$\alpha_{ij}=(-1)^{i+j} \left\vert \begin{array}{cccccccc} a_{11} ... ..._{n2} &\cdots &a_{n j-1} &a_{n j+1} &\cdots &a_{nn} \end{array} \right\vert$$ (5.13)

To simplify the arithmetic, a n-rowed determinant can be calculated with the sum of n (n$-$1)-rowed sub-determinants (5.13). The determinant can be expanded from a row or column. If for example $k$ is the number of any column, than the determinant is

$$\left\vert \begin{array}{cccc} a_{11} &a_{12} &\cdots &a_{1n} a... ...n2} &\cdots &a_{nn} \end{array} \right\vert=\sum_{i=1}^n a_{ik} \alpha_{ik}.$$ (5.14)

With this rule (5.14) the coefficient $a_0$ (5.11) can be simplified, because only $a_{0,11}=1$, and the rest in the first column is 0. After expanding the determinant from the first row one obtains

$$a_0= \frac{1}{\mathrm{Det}(\mathbf{D})} \left\vert \begin{array}{cccc} x_1 &y_1 &z_1 x_2 &y_2 &z_2 x_3 &y_3 &z_3 \end{array} \right\vert.$$ (5.15)

In the same way all other coefficients of the shape function $N_0$ can be determined

$$b_0= \frac{-1}{\mathrm{Det}(\mathbf{D})} \left\vert \begin{array}... ...array}{cccc} 1 &x_1 &y_1 1 &x_2 &y_2 1 &x_3 &y_3 \end{array} \right\vert.$$ (5.16)

The coefficients for the other shape functions can be found in the same way, the only difference leis in the equation system (5.9). Because, for example, $N_1(x,y,z)=a_1 + b_1 x + c_1 y + d_1 z$ must be 1 in node P$_1$ and 0 in all other nodes, it can be formulated

$$a_1 + b_1 x_0 + c_1 y_0 + d_1 z_0=0$$

$$a_1 + b_1 x_1 + c_1 y_1 + d_1 z_1=1$$

$$a_1 + b_1 x_2 + c_1 y_2 + d_1 z_2=0$$ (5.17)

$$a_1 + b_1 x_3 + c_1 y_3 + d_1 z_3=0 \nonumber$$

With the above described procedure also the coefficients for $N_1(x,y,z)$ can be determined straightforwardly to

$$b_1= \frac{1}{\mathrm{Det}(\mathbf{D})} \left\vert \begin{array}{... ...ay}{cccc} 1 &x_0 &z_0 1 &x_2 &z_2 1 &x_3 &z_3 \end{array} \right\vert.$$ (5.18)

5.2.2 Coordinate Transformation

A coordinate transformation can help to simplify the calculation of integrals. For constructing the residual the calculation of the following element integral is frequently needed

$$I_e=\int_T N_i(x,y,z)N_j(x,y,z) dz dy dz \qquad \mathrm{where} \qquad i,j=0,1,2,3.$$ (5.19)

Here the multiplication of two (linear) form functions leads to a more complex polynomial which complicates the integration over the region. It is more practical to integrate over a normalized element $T^n$ (see Fig. 5.2). For this purpose, a tetrahedron with any location in the global (x,y,z)-coordinate system must be transformed into a normalized local $(\xi ,\eta ,\zeta )$-coordinate system.

Each point (x,y,z) of the tetrahedral element in the global coordinate system can be transformed to a corresponding point $(\xi ,\eta ,\zeta )$ in the normalized coordinate system with the following bijective projection rule

$$x=x_0+(x_1-x_0)\xi+(x_2-x_0)\eta+(x_3-x_0)\zeta,$$

$$y=y_0+(y_1-y_0)\xi+(y_2-y_0)\eta+(y_3-y_0)\zeta,$$ (5.20)

$$z=z_0+(z_1-z_0)\xi+(z_2-z_0)\eta+(z_3-z_0)\zeta.$$

This projection in matrix form leads to

$$\{r\}=\{r_0\}+\mathbf{J}\cdot\{\delta\}$$ (5.21)

and the conversion from the global to the normalized coordinates is

$$\{\delta\}=\mathbf{J}^{-1}\big(\{r\}-\{r_0\}\big),$$ (5.22)

with $\{r\}=(x,y,z)^T$ and $\{\delta\}=(\xi,\eta,\zeta)^T$.
$\mathbf{J}$ is the so-called Jacobian matrix which only depends on the global coordinates (x,y,z)

$$\mathbf{J}= \left[ \begin{array}{ccc} x_1 - x_0 & \;x_2 - x_0 & \... ...y_0 & \;y_3 - y_0 z_1 - z_0 & \;z_2 - z_0 & \;z_3 - z_0 \end{array} \right].$$ (5.23)

The element integral (5.19) calculated in the normalized coordinate system must be multiplied with the determinant of the Jacobian matrix

$$I_e=\mathrm{Det}(\mathbf{J})\int_{T^n} N_i^n(\xi,\eta,\zeta)N_j^n(\xi,\eta,\zeta) d\xi d\eta d\zeta \qquad \mathrm{where} \qquad i,j=0,1,2,3,$$ (5.24)

because the following relationship holds

$$\frac{\partial (\xi,\eta,\zeta )}{\partial (x,y,z) }=\mathrm{Det}(\mathbf{J}).$$ (5.25)

The shape functions for the normalized tetrahedron $T^n$ are simpler than those in the global coordinates, because they are reduced to [91]

$$N_0^n(\xi,\eta,\zeta) =1-\xi-\eta-\zeta$$

$$N_1^n(\xi,\eta,\zeta) =\xi$$ (5.26)

$$N_2^n(\xi,\eta,\zeta) =\eta$$

$$N_3^n(\xi,\eta,\zeta) =\zeta.$$

These shape functions lead to a simpler integrand in (5.24). A further advantage is that after normalization the lower integration limit is always 0, because as shown in Fig. 5.2, the tetrahedron $T^n$ starts in the origin of ordinates.

Also the upper limits can be found straightforwardly. As shown in Fig. 5.2, the element $T^n$ is bounded by a plane which goes through the points P$_1$(1,0,0), P$_2$(0,1,0), and P$_3$(0,0,1). This plane is described with the equation $\xi+\eta+\zeta=1$. The maximum on the $\xi$-axes is 1. The limit in the $\xi$-$\eta$-plane ($\zeta=0$) can be described with $\eta(\xi)=1-\xi$ and the limit in $\zeta$-direction is $\zeta(\xi,\eta)=1-\xi-\eta$.

With these limits the element integral of $T^n$ can be written in the form

$$I_e=\mathrm{Det}(\mathbf{J})\int_{\xi=0}^{1}\int_{\eta=0}^{1-\xi\... ...xtrm{for } i= j \frac{1}{120} \qquad\textrm{for } i\neq j \end{array} \right.$$ (5.27)

This is a real advantage of the normalized tetrahedron, because there exists a simple scheme for $I_e$. The result from integration over the element $T^n$ is either $\frac{1}{60}$ or $\frac{1}{120}$, only depending, if the two functions $N_{i}$ and $N_{j}$ are equal or not. So in fact, to calculate the integral $I_e$ in the normalized coordinate system, only the determinant of the Jacobian matrix must be calculated. This procedure is much easier than to find the element integral for a common tetrahedron in the global coordinate system, where each element has a different size. This means that each integral $I_e$ has a different result and must be calculated separately.

5.2.3 Differentiation in the Normalized Coordinate System

The differentiation of the whole projection (5.2.2), with respect to $x$ leads to

$$1=(x_1-x_0)\frac{\partial \xi}{\partial x}+(x_2-x_0)\frac{\partial \eta}{\partial x}+(x_3-x_0)\frac{\partial \zeta}{\partial x},$$

$$0=(y_1-y_0)\frac{\partial \xi}{\partial x}+(y_2-y_0)\frac{\partial \eta}{\partial x}+(y_3-y_0)\frac{\partial \zeta}{\partial x},$$ (5.28)

$$0=(z_1-z_0)\frac{\partial \xi}{\partial x}+(z_2-z_0)\frac{\partial \eta}{\partial x}+(z_3-z_0)\frac{\partial \zeta}{\partial x},$$

with respect to $y$ leads to

$$0=(x_1-x_0)\frac{\partial \xi}{\partial y}+(x_2-x_0)\frac{\partial \eta}{\partial y}+(x_3-x_0)\frac{\partial \zeta}{\partial y},$$

$$1=(y_1-y_0)\frac{\partial \xi}{\partial y}+(y_2-y_0)\frac{\partial \eta}{\partial y}+(y_3-y_0)\frac{\partial \zeta}{\partial y},$$ (5.29)

$$0=(z_1-z_0)\frac{\partial \xi}{\partial y}+(z_2-z_0)\frac{\partial \eta}{\partial y}+(z_3-z_0)\frac{\partial \zeta}{\partial y},$$

and with respect to $z$ leads to

$$0=(x_1-x_0)\frac{\partial \xi}{\partial z}+(x_2-x_0)\frac{\partial \eta}{\partial z}+(x_3-x_0)\frac{\partial \zeta}{\partial z},$$

$$0=(y_1-y_0)\frac{\partial \xi}{\partial z}+(y_2-y_0)\frac{\partial \eta}{\partial z}+(y_3-y_0)\frac{\partial \zeta}{\partial z},$$ (5.30)

$$1=(z_1-z_0)\frac{\partial \xi}{\partial z}+(z_2-z_0)\frac{\partial \eta}{\partial z}+(z_3-z_0)\frac{\partial \zeta}{\partial z}.$$

These derivatives can be also expressed in matrix form

$$\mathbf{I}=\mathbf{J}\times \left[ \begin{array}{ccc} \displaysty... ...partial y} &\displaystyle\frac{\partial \zeta}{\partial z} \end{array} \right],$$ (5.31)

and the following relationship for the partial differential operators in the normalized system can be found with

$$\mathbf{J}^{-1}= \left[ \begin{array}{ccc} \displaystyle\frac{\pa... ...partial y} &\displaystyle\frac{\partial \zeta}{\partial z} \end{array} \right].$$ (5.32)

The inverse of the Jacobian matrix (5.23) in the global coordinate system is given by

$$\mathbf{J}^{-1}= \frac{1}{\mathrm{Det}(\mathbf{J})} \left[ \begin... ... &L_{13} L_{21} &L_{22} &L_{23} L_{31} &L_{32} &L_{33} \end{array} \right],$$ (5.33)

where the components $L_{ij}$ of the inverse matrix are

$$L_{11}=-J_{23} J_{32}+J_{22} J_{33},$$

$$L_{12}=\phantom{-}J_{13} J_{32}-J_{12} J_{33},$$

$$L_{13}=-J_{13} J_{22}+J_{12} J_{23},$$

$$L_{21}=\phantom{-}J_{23} J_{31}-J_{21} J_{33},$$

$$L_{22}=-J_{13} J_{31}+J_{11} J_{33}$$ (5.34)

$$L_{23}=\phantom{-}J_{13} J_{21}-J_{11} J_{23},$$

$$L_{31}=-J_{22} J_{31}+J_{21} J_{32}$$

$$L_{32}=\phantom{-}J_{12} J_{31}-J_{11} J_{23},$$

$$L_{33}=-J_{12} J_{21}+J_{11} J_{22}. \nonumber$$

The components $J_{ij}$ of the Jacobian matrix (5.23) depend only on the location of the tetrahedron vertices in the global coordinate system. Because of

$$\frac{\partial \xi}{\partial x}=\frac{L_{11}}{\mathrm{Det}(\mathb... ...quad \frac{\partial \zeta}{\partial z}=\frac{L_{33}}{\mathrm{Det}(\mathbf{J})},$$ (5.35)

there exists a relationship between the partial differential operators in the normalized system and the coordinates of the four nodes from the global system.

If any continuous function $f(\xi,\eta,\zeta)$ in the normalized coordinate system is differentiated in respect to $x$, then the chain rule must be used so that

$$\frac{\partial f}{\partial x}= \frac{\partial f}{\partial \xi}\fr... ...partial x}+ \frac{\partial f}{\partial \zeta}\frac{\partial \zeta}{\partial x}.$$ (5.36)

The gradient of the function $f(\xi,\eta,\zeta)$ in the normalized system becomes to

$$\nabla f(\xi,\eta,\zeta)= \left[ \begin{array}{c} \displaystyle\... ...a} [2.8mm] \displaystyle\frac{\partial f}{\partial \zeta} \end{array} \right]$$ (5.37)

This means that the gradient operator $\nabla^n$ in the normalized system must be multiplied with the transposed of the inverse Jacobian matrix

$$\nabla=(\mathbf{J}^{-1})^\mathbf{T}\times\nabla^n= \frac{1}{\math... ...ta} [2.8mm] \displaystyle\frac{\partial }{\partial \zeta} \end{array} \right]$$ (5.38)

5.2.4 Discretization of the Oxidant Diffusion

In the continuum formulation from (3.2), the diffusion of oxidants through the oxide material is described with

$$D \Delta C = \eta k_{max} C.$$ (5.39)

When Galerkin's method is applied, it is multiplied with a weight function $N_j(x,y,z)$ and integrated over the domain $\Omega$, which leads to

$$\int_\Omega N_j(x,y,z) D \Delta C d\Omega = \int_\Omega N_j(x,y,z) \eta k_{max} C d\Omega.$$ (5.40)

If there are two functions $u(\vec{x})$ and $v(\vec{x})$ defined on a domain $\Omega$, Green's theorem says that

$$\int_\Omega \nabla u \nabla v d\Omega + \int_\Omega u \Delta v d\Omega= \int_\Gamma u \frac{\partial v}{\partial n} d\Gamma,$$ (5.41)

where $\Gamma$ is the boundary of the domain.

With Green's theorem the Galerkin formulation from (5.40) can be rewritten in the form

$$-D\int_\Omega \nabla N_j \nabla C d\Omega + \int_\Gamma N_j \frac{\partial C}{\partial n} d\Gamma= k_{max}\int_\Omega N_j \eta C d\Omega.$$ (5.42)

Here the diffusion coefficient $D$ and the maximal reaction rate $k_{max}$ do not directly depend on the location and, therefore, they do not need to be integrated over space and can stand outside of the integral. Furthermore, it is assumed that there is no flow of oxidants through the boundary surface and the boundary condition becomes

$$\int_\Gamma N_j \frac{\partial C}{\partial n} d\Gamma= 0$$ (5.43)

With this Neumann boundary condition the Galerkin formulation for the oxidant diffusion can be reduced to

$$-D\int_\Omega \nabla N_j \nabla C d\Omega = k_{max}\int_\Omega N_j \eta C d\Omega\qquad \mathrm{for} \qquad j=0,1,2,3.$$ (5.44)

With the finite element method it can be assumed that this equation is only valid on a single tetrahedral element $T$. Furthermore, the scalar function for the oxidant concentration $C(\vec{x},t)$ is here approximated linearly with

$$C(\vec{x},t=t_n)=\sum_{i=0}^3 c_{i}^{(t_n)} N_i(x,y,z),$$ (5.45)

where $c_{i}^{(t_n)}$ is the value of the oxidant concentration in node $i$ and at discrete time $t_n$. $N_i(x,y,z)$ is the respective shape function from this node.
The distribution of the normalized silicon $\eta(\vec{x},t)$ is approximated in the same way so that

$$\eta(\vec{x},t=t_n)=\sum_{i=0}^3 \eta_{i}^{(t_n)} N_i(x,y,z),$$ (5.46)

where $\eta_{i}^{(t_n)}$ is the value of the normalized silicon in node $i$ and at discrete time $t_n$. $N_i(x,y,z)$ is the linear shape function (5.8) for node $i$.

With the approximation for $C(\vec{x},t)$ and $\eta(\vec{x},t)$, the oxidant diffusion on a single element $T$ can be described with

$$-D\int_T \nabla N_j \nabla \Big(\sum_{i=0}^3 c_{i}^{(t_n)} N_i\... ... \sum_{i=0}^3 \eta_{i}^{(t_n)} N_i \sum_{i=0}^3 c_{i}^{(t_n)} N_i d\Omega.$$ (5.47)

In the approximation for $C(\vec{x},t)$ and $\eta(\vec{x},t)$ the shape function $N_i(x,y,z)$ is the same. Since the values $c_{i}$ and $\eta_{i}$ in the nodes do not depend on the spatial location, (5.47) can be rewritten as

$$-D \sum_{i=0}^3 \Big(c_{i}^{(t_n)} \int_T \nabla N_j \nabla N_i\... ...um_{i=0}^3 \Big(c_{i}^{(t_n)} \eta_{i}^{(t_n)} \int_T N_j N_i d\Omega\Big),$$ (5.48)

where it is assumed that $\sum \eta_i N_i \sum c_i N_i \approx \sum \eta_i c_i N_i$.
By substituting the integrals with

$$K_{ij}=\int_T \nabla N_i \nabla N_j d\Omega$$

$$M_{ij}=\int_T N_i N_j d\Omega$$ (5.49)

the discretized equation for the oxidant diffusion is simplified to

$$-D \sum_{i=0}^3 K_{ij} c_{i}^{(t_n)} = k_{max}\sum_{i=0}^3 M_{ij} c_{i}^{(t_n)} \eta_{i}^{(t_n)} \qquad \mathrm{for} \qquad j=0,1,2,3.$$ (5.50)

The components $M_{ij}$ were already calculated in Section 5.2.2. After integration of $\int_T N_i N_j d\Omega$ in the normalized coordinate system, it was found that

$$M_{ij}=\mathrm{Det}(\mathbf{J})\cdot \left\{ \begin{array}{ll} \f... ...xtrm{for } i= j \frac{1}{120} \qquad\textrm{for } i\neq j \end{array} \right.$$ (5.51)

For calculating the components $K_{ij}$ the integral is also transformed from the global to the normalized coordinate system, and with (5.25) follows

$$K_{ij}=\mathrm{Det}(\mathbf{J})\int_{\xi=0}^{1}\int_{\eta=0}^{1-\xi\;}\int_{\zeta=0}^{1-\xi-\eta}\nabla N_i^n \nabla N_j^n d\zeta d\eta d\xi,$$ (5.52)

where $N_{i,j}^n(\xi,\eta,\zeta)$ are the shape functions for a normalized tetrahedron $T^n$ from (5.26).
It was demonstrated in (5.38) that the transformation of the gradient operator $\nabla$ is carried out by a multiplication with the matrix $(\mathbf{J}^{-1})^\mathbf{T}$, so that the integrand from (5.52) becomes

$$\nabla N_i^n \nabla N_j^n =(\mathbf{J}^{-1})^{\mathbf{T}} \nabla^n N_i^n (\mathbf{J}^{-1})^{\mathbf{T}} \nabla^n N_j^n =$$

$$=\frac{1}{\mathrm{Det}(\mathbf{J})} \left[ \begin{array}{ccc} L_{... ...2.8mm] \displaystyle\frac{\partial N_j^n }{\partial \zeta} \end{array} \right]=$$ (5.53)

$$= \frac{1}{\mathrm{Det}(\mathbf{J})} \left[ \begin{array}{c} L_{1... ...partial \eta} +L_{33}\frac{\partial N_j^n }{\partial \zeta} \end{array} \right]$$

After the multiplication of the two vectors and rearranging of this scalar product, (5.2.4) is

$${\phantom{+}\frac{\partial N_i^n }{\partial \xi}}{\frac{\partial ... ...ac{\partial N_j^n }{\partial \eta}} (L_{11}L_{21}+ L_{12}L_{22}+ L_{13}L_{23})+$$

$$+ {\frac{\partial {N_i^n} }{\partial \xi}} \frac{\partial {N_j^n}... ...ac{\partial {N_j^n} }{\partial \xi} (L_{21}L_{11}+ L_{22}L_{12}+ L_{23}L_{13})+$$

$$+ {\frac{\partial {N_i^n} }{\partial \eta}}{\frac{\partial {N_j^n... ...\partial {N_j^n} }{\partial \zeta}} (L_{21}L_{31}+ L_{22}L_{32}+ L_{23}L_{33})+$$ (5.54)

$$+ {\frac{\partial {N_i^n} }{\partial \zeta}}{\frac{\partial {N_j^... ...l \zeta}}{\frac{\partial N_j^n }{\partial \zeta}} (L_{31}^2+L_{32}^2+L_{33}^2).$$

Here the $L_{xy}$-terms which only depend on the location of the nodes in the global coordinate system, can be replaced by six constant coefficients $G_A - G_F$

$$\phantom{xxxxxxxxxxxxx}G_A = L_{11}^2+L_{12}^2+L_{13}^2$$

$$G_B = L_{11}L_{21}+ L_{12}L_{22}+ L_{13}L_{23}$$

$$G_C = L_{11}L_{31}+ L_{12}L_{32}+ L_{13}L_{33}$$

$$G_D = L_{21}^2+ L_{22}^2+ L_{23}^2$$ (5.55)

$$G_E = L_{21}L_{31}+ L_{22}L_{32}+ L_{23}L_{33}$$

$$G_F = L_{31}^2+L_{32}^2+L_{33}^2.$$

With (5.55) the scalar product (5.54) can be simplified to

$${\phantom{+}\frac{\partial N_i^n }{\partial \xi}}{\frac{\partial ... ...l {N_i^n} }{\partial \zeta} \frac{\partial {N_j^n} }{\partial \xi}\Big) G_C +$$

$$+ \frac{\partial {N_i^n} }{\partial \eta}{\frac{\partial {N_j^n} ... ...\partial N_i^n }{\partial \zeta}}{\frac{\partial N_j^n }{\partial \zeta}} G_F.$$ (5.56)

For example, the simplified scalar product (5.56) for $i=j=0$, $(N_0=1-\xi-\eta-\zeta)$ is

$$G_A+2 G_B+2 G_C+G_D+2 G_E+G_F,$$ (5.57)

and for $i=0$, $j=1$ $(N_1=\xi)$ it is

$$-G_A-G_B-G_C.$$ (5.58)

For all combinations of $i,j=0,1,2,3$ the simplified scalar product can be written in the form

$$\mathbf{S_A} G_A+\mathbf{S_B} G_B+\mathbf{S_C} G_C+\mathbf{S_D} G_D+\mathbf{S_E} G_E +\mathbf{S_F} G_F,$$ (5.59)

Instead of finding the scalar product (5.56), and the components of the matrices $\mathbf{S_A}-\mathbf{S_F}$ for all combinations of $i,j=0,1,2,3$ by the way like in (5.57), it is more comfortable to use

$$\mathbf{N^n}= \left[ \begin{array}{c} 1-\xi-\eta-\zeta \xi \e... ...n{array}{c} -1 \phantom{-}0 \phantom{-}0 \phantom{-}1 \end{array}\right],$$ (5.60)

and get the matrices $\mathbf{S_A}-\mathbf{S_F}$ with

$$\mathbf{S_A}=\frac{\partial \mathbf{N^n} }{\partial \xi}\times\Bi... ...}0 \phantom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0 \end{array}\right],$$ (5.61)

$$\mathbf{S_B}= \left[ \begin{array}{c} -1 \phantom{-}1 \phant... ...ntom{-}0 &\phantom{-}0 &\phantom{-}0 &\phantom{-}0 \end{array}\right],\; \ldots$$ (5.62)

Because the derivatives of the linear shape functions $N_{i,j}^n(\xi,\eta,\zeta)$ can only result in the values $-1$, 0 or $+1$, the integral from (5.52) becomes

$$\int_{\xi=0}^{1}\int_{\eta=0}^{1-\xi\;}\int_{\zeta=0}^{1-\xi-\eta... ...a=0}^{1-\xi-\eta} \pm 1 d\zeta d\eta d\xi = \pm\frac{1}{6}, %$\mathbf{D}$$$ (5.63)

which means that all matrices $\mathbf{S_A}-\mathbf{S_F}$ must be weighted with $\frac{1}{6}$.

After finding $\nabla N_i^n \nabla N_j^n$ and integration over the (normalized) element, the coefficients $K_{ij}$ are

$$K_{ij}=\frac{1}{6}(\mathbf{S_A}_{,ij} G_A+\mathbf{S_B}_{,ij} G_... ...,G_D+\mathbf{S_E}_{,ij} G_E+\mathbf{S_F}_{,ij} G_F)/\mathrm{Det}(\mathbf{J}).$$ (5.64)

Galerkin's method assumes that the residual from the (discretized) equation of the oxidant diffusion is zero, and so (5.50) is rewritten in the form

$$\sum_{i=1}^4\Big(D K_{ij} c_{i}^{(t_n)}+ k_{max} M_{ij} c_{i}^{(t_n)} \eta_{i}^{(t_n)}\Big)=0,\qquad j=0,1,2,3.$$ (5.65)

This is a system with four equations, but with eight unknown variables $c_{i}^{(t_n)}$ and $\eta_{i}^{(t_n)}$. Because of more unknowns than equations it is impossible to solve this equation system in present form. In the next section the required 4 equations introduced.

5.2.5 Discretization of the $\boldsymbol{\eta}$-Dynamics

The dynamics of the normalized silicon concentration $\eta$ in the continuum formulation (3.5) is described by

$$\frac{\partial \eta} {\partial t} = - \frac{1}{\lambda}\eta k_{max} C/N_1,$$ (5.66)

After applying Galerkin's method with a weight function $N_j(x,y,z)$ one obtains in a domain $\Omega$ is

$$\int_\Omega N_j \frac{\partial \eta} {\partial t} d\Omega = -K_... ...th} \quad K_A=\frac{k_{max}}{\lambda N_1} \qquad \mathrm{and} \quad j=0,1,2,3.$$ (5.67)

Because of the time dependence in this equation an additional time discretization of the term $\frac{\partial \eta(\vec{x},t)} {\partial t}$ is necessary. This time discretization is performed with a simple backward-Euler method as

$$\frac{\partial \eta(\vec{x},t=t_n)} {\partial t} = \frac{\eta(\vec{x},t_n)-\eta(\vec{x},t_{n-1})} {\Delta t},$$ (5.68)

where $t_n$ and $t_{n-1}$ are two successive discrete times.

For an equation $\frac{\partial y(t)}{\partial t} = f(y,t)$ two Euler methods can be applied. The forward-Euler method is an explicit simple method, because the new value $y^{(t_{n+1})}=y^{t_n}+\Delta t\cdot f(t_{n},y^{t_{n}})$ is defined in terms of values that are already known [92]. The backward-Euler method comes from using $f(y,t)$ at the end of a time step, when $t=t_{n+1}$. It is an implicit method, because in order to obtain the new discrete value $y^{(t_{n+1})}$ a linear equation of the form $y^{(t_{n+1})}=y^{t_n}+\Delta t\cdot f(t_{n+1},y^{(t_{n+1})})$ must be solved [92], which requires additional computing time. But compared with the forward-Euler method the most important advantage of the backward-Euler method is that a much larger time step size $\Delta t$ can be used. The reason is that implicit methods are usually much more stable for solving a stiff equation. A stiff equation is a differential equation for which certain numerical methods for solving the equation are numerically unstable, unless the time step size is taken to be extremely small [93].

The spatial approximation for the oxdiant concentration $C(\vec{x},t)$ and the normalized silicon $\eta(\vec{x},t)$ for one finite element $T$ is the same as in Section 5.2.4

$$C(\vec{x},t=t_n)=\sum_{i=0}^3 c_{i}^{(t_n)} N_i(x,y,z),$$ (5.69)

$$\eta(\vec{x},t=t_n)=\sum_{i=0}^3 \eta_{i}^{(t_n)} N_i(x,y,z).$$ (5.70)

With the time discretization (5.68) and the spatial approximation for $C(\vec{x},t)$ and $\eta(\vec{x},t)$, the Galerkin formulation for the $\eta$-dynamics on a finite element $T$ becomes

$$\int_T N_j\frac{1}{\Delta t}\Big(\sum_{i=0}^3 \eta_{i}^{(t_n)}N_i... ... N_j \sum_{i=0}^3 \eta_{i}^{(t_n)}N_i \sum_{i=0}^3 c_{i}^{(t_n)}N_i d\Omega,$$ (5.71)

Because of the same shape function $N_i$ for $C(\vec{x},t)$ and $\eta(\vec{x},t)$ and no spatial dependence of $c_{i}^{(t_n)}$ and $\eta_{i}^{(t_n)}$, the last equation can be rearranged to

$$\frac{1}{\Delta t}\sum_{i=0}^3\Big(\big(\eta_{i}^{(t_n)} -\eta_{i... ...A\sum_{i=0}^3\Big(\eta_{i}^{(t_n)} c_{i}^{(t_n)} \int_T N_i N_j d\Omega\Big),$$ (5.72)

After substituting $M_{ij}= \int_T N_i N_j d\Omega$ from (5.51) the discretized equation for the $\eta$-dynamics is simplified to

$$\frac{1}{\Delta t}\sum_{i=0}^3 M_{ij}\big((\eta_{i}^{(t_n)} -\eta_{i}^{(t_{n-1})}\big)= -K_A\sum_{i=0}^3 M_{ij}\eta_{i}^{(t_n)} c_{i}^{(t_n)}.$$ (5.73)

In order to fulfill Galerkin's demand that the residual should be zero, the last equation can be rewritten as

$$\sum_{i=0}^3\Big(M_{ij}\big(\eta_{i}^{(t_n)} -\eta_{i}^{(t_{n-1})... ...i}^{(t_n)} c_{i}^{(t_n)} \Delta t\Big)=0\qquad \mathrm{for} \qquad j=0,1,2,3,$$ (5.74)

which is also a system with four equations and eight unknown variables $c_{i}^{(t_n)}$ and $\eta_{i}^{(t_n)}$. The values for $\eta_{i}^{(t_{n-1})}$ are already determined at the previous time step.

By combining the two equation systems (5.65) and (5.74), a non-linear but fully determined equation system for one finite element, with 8 equations and the 8 unknowns $c_{0}^{(t_n)}-c_{3}^{(t_n)}$ and $\eta_{0}^{(t_n)}-\eta_{3}^{(t_n)}$, is obtained. The system is non-linear because of the product $\eta_{i}^{(t_n)} c_{i}^{(t_n)}$ in (5.65) and in (5.74). The complete equation system can be solved (for example with the Newton method) at each time point $t_n$ and the values for $c_i^{(t_n)}$ and $\eta_i^{(t_n)}$ can be determined.

5.2.6 Discretization of the Mechanics

The main interest in the continuum mechanics is the deformation of a body by internal or external forces. The deformation is expressed by the displacements $d(x,y,z)$. The displacement of a point in a three-dimensional elastic continuum is defined by three displacement components $u(x,y,z)$, $u(x,y,z)$, and $u(x,y,z)$ in directions of the three coordinates $x$, $y$, and $z$, so that

$$\vec{d}(x,y,z)= \left\{ \begin{array}{c} u(x,y,z) v(x,y,z) w(x,y,z) \end{array}\right\}$$ (5.75)

In contrast to the previous differential equations (5.39) and (5.66) for the mechanics Galerkin's method is not needed. Instead the virtual work concept is used [94]. The displacement components $u$, $v$, and $w$ are directly discretized on a finite tetrahedral element

$$u(x,y,z)=\sum_{i=0}^3 u_{i} N_i(x,y,z),$$

$$v(x,y,z)=\sum_{i=0}^3 v_{i} N_i(x,y,z),$$ (5.76)

$$w(x,y,z)=\sum_{i=0}^3 w_{i} N_i(x,y,z),$$

where $u_i$, $v_i$, and $w_i$ are the displacement values in $x$-, $y$-, and $z$-direction on node $i$ and $N_i$ is the linear shape function from (5.8).

The strain components in the elastic case are first order derivatives of the displacement components

$$\tilde{\varepsilon}^e = \left\{ \begin{array}{c} \varepsilon_{xx}... ...B} \vec{d}^e=[\mathbf{B_0},\mathbf{B_1},\mathbf{B_2},\mathbf{B_3}] \vec{d^e}.$$ (5.77)

The element displacement is defined by the 12 displacement components of the 4 nodes as

$$\vec{d}^e = \left\{ \begin{array}{ccc} \vec{d_0} \vec{d_1} \v... ...begin{array}{ccc} u_0 v_0 w_0 \end{array} \right\} \qquad \textnormal{etc.}$$ (5.78)

The submatrix $\mathbf{B_i}$ of displacement derivatives for the node $i$ is [95]

$$\mathbf{B_i}= \left[ \begin{array}{ccc} \displaystyle{\frac{\part... ... d_i \;&0 \;&b_i \end{array} \right], \qquad \mathrm{where} \qquad i=0,1,2,3.$$ (5.79)

The coefficients $b_i$, $c_i$, and $d_i$ are the same as already presented in (5.16) and (5.18).

The entire inner virtual work on a continuous elastic body, and so also on a finite element is [94]

$$W_{inner}=\int_{V}\{\tilde{\varepsilon}^e\}^T \sigma^e dV, \vspace*{-2mm}$$ (5.80)

with the stress tensor $\tilde{\sigma}$ (3.15). Here it is assumed that there is no residual stress $\tilde{\sigma_0}$, and the stress tensor is

$$\tilde{\sigma} = \mathbf{D} (\tilde{\varepsilon} - \tilde{\varepsilon_0}).$$ (5.81)

In discretized form the transposed strain tensor is

$$\{\tilde{\varepsilon}^e\}^T = \vec{d^e}^T \mathbf{B^T}$$ (5.82)

After discretization the stress tensor (5.81) can be arranged as a function of the element displacement vector

$$\sigma^e = \mathbf{D} (\tilde{\varepsilon}^e-\tilde{\varepsilon_0}^e) = \mathbf{D} \mathbf{B} \vec{d^e}-\mathbf{D} \tilde{\varepsilon_0}^e.$$ (5.83)

Together with the transposed strain tensor (5.82), this stress tensor leads to the following discretized form of the equation for the inner virtual work on a finite element

$$W_{inner}=\vec{d^e}^T\int_{V}\big(\mathbf{B^T}\mathbf{D}\mathbf{B} \vec{d^e}-\mathbf{B^T}\mathbf{D} \tilde{\varepsilon_0}^e\big)dV.$$ (5.84)

The outer virtual work on a finite element, caused by the external nodal forces is

$$W_{outer}=\vec{d^e}^T\vec{f^e}^{ext}=0,$$ (5.85)

because it is assumed that during the oxidation process there are not external forces acting.

On any elastic body, and so also on a finite element, the inner work must be equal with the outer work

$$W_{inner}=\vec{d^e}^T\int_{V}\big(\mathbf{B^T}\mathbf{D}\mathbf{B} \vec{d^e}-\mathbf{B^T}\mathbf{D} \tilde{\varepsilon_0}^e\big)dV=0=W_{outer},$$ (5.86)

which can be simplified to

$$\int_{V}\mathbf{B^T}\mathbf{D}\mathbf{B} \vec{d^e} dV=\int_{V}\mathbf{B^T}\mathbf{D} \tilde{\varepsilon_0}^e dV.$$ (5.87)

Here the integrals can be substituted as sketched in [95]

$$\mathbf{K^e}=\int_{V}\mathbf{B^T}\mathbf{D}\mathbf{B} dV= \mathbf{B^T}\mathbf{D}\mathbf{B} V^e,$$ (5.88)

$$\vec{f^e_{int}}=\int_{V}\mathbf{B^T}\mathbf{D} \tilde{\varepsilon_0}^e dV=\mathbf{B^T}\mathbf{D} \tilde{\varepsilon_0}^e V^e,$$ (5.89)

where $\mathbf{K^e}$ is the so-called stiffness matrix and $\vec{f^e_{int}}$ can be declared as internal force vector. Since the integrands are not functions of the $x$-, $y$-, or $z$-coordinates, the integration over the volume is equal with its much more simpler multiplication. The volume of any tetrahedron in the global coordinate system can be calculated with the determinant of matrix (5.12) by

$$V^e=\frac{1}{6} {\mathrm{Det}(\mathbf{D})}$$ (5.90)

The most important fact is that the residual strain tensor $\tilde{\varepsilon_0}^e$ in (5.89) loads the mechanical system. Because the residual strain components $\varepsilon_{0,ii}$ are directly proportional to the normalized additional volume (3.37), there is a relationship between the volume expansion and the internal nodal forces.

With the integral substitutions (5.88) and (5.89), the balance equation (5.87) becomes a linear equation system for the mechanical problem on one finite element

$$\mathbf{K^e} \vec{d^e}=\vec{f^e}$$ (5.91)

The system is fully determined, because there are 12 equations and also 12 unknown displacement-components (three on each node) on the tetrahedron.