4.2.3 Derivation of Equations over an Element

The derivation of algebraic equations in a finite element approach associated with the model DE, defined by equations (4.1) and (4.2), over a mesh element involves a three-step approach [125]:

  1. Construction of the weak formulation of the DE;
  2. Choice of the form of the approximate solution over the finite element;
  3. Derivation of the finite element equations in matrix form by substituting the approximate solution into the weak formulation.

In the following, the construction of the algebraic equations over the isolated element Ω1 with two nodes 1 and 2, whose endpoints have coordinates x=x1 and x=x1, is presented.

Step 1 - Weak Formulation

Different variational methods, such as the Ritz-Galerkin and weighted-residual, provide a background for the development of the algebraic equations over the element [116]. In all these methods, the basic idea is to seek a polynomial approximation φh1 of the solution to equation (4.1) within the finite element Ω1 in the form

\[\begin{equation} \varphi^{1}_\text{h}= \sum_{k=1}^{2} c^{1}_\text{k}u^{1}_\text{k}(x), \end{equation}\] (4.3)

where ck1 are the unknown nodal values of the solution φ(x)) at the nodes of the element Ω1 and uk1 are the basis functions over the element. The scope of any variational method is to determine the nodal values ck1 in order to find an approximate solution φh1 that best satisfies the DE at every point of the element Ω1 of interest. The strong solution to equation (4.1) must be as smooth as required by the DE. In the present case, the solution should be continuously differentiable until at least the second derivative. To avoid such a requirement, the weak formulation of the problem is introduced. The main advantage compared to the strong form is that the weak form lowers the demand on continuity of the solution to only the first derivative.

The construction of the weak formulation of the problem requires one to recast the DE in a weighted-integral form [125]. The weighted-integral statement of the DE is constructed as follows

\[\begin{equation} \int_{x_\text{1}}^{x_\text{2}} w_\text{j}(x)R(x,c^{1}_\text{k}) \text{d} x=0,\quad \text{for} \quad j=1,2,...,n, \end{equation}\] (4.4)

where

\[\begin{equation} R(x,c^{1}_\text{k})=-\frac{d}{dx}\left(a\frac{d{\varphi^{1}_\text{h}}}{dx}\right) +b\varphi^{1}_\text{h}-f. \end{equation}\] (4.5)

R(x,ck1) is called residual of approximation of the DE within the element Ω1 and wj(x) is the weight function. The integral statement of equation (4.4) allows one to choose n linearly independent functions wj(x) and obtain a set of n equations for ck1 by making R(x,ck1) equal to zero at the selected points of the element. It therefore requires that the basis function uk1 be such that the approximate solution φh1 is twice differentiable. Note that the weighted-integral statement of (4.4) is equivalent to the DE and does not include any boundary conditions. The boundary conditions will come into play subsequently.

Due to the different choices of the weight function wj(x) and/or the basis function uk1, the system of algebraic equations within the element will have different characteristics when different variational methods are used [125]. A general weighted-residual method is a variational method which requires only the weighted-integral statement of the DE to determine the nodal values ck1 of the approximation. The weight functions wj(x) can be chosen from an independent set of functions, but are required to be linearly independent, otherwise the set of n equations for ck1 will not be solvable due to the underdetermination of the system of equations [88]. The Ritz-Galerkin method is a weighted-residual method which uses a set of weight functions wj(x) equal to the set of basis functions uk1 such that the weighted-integral statement (4.4) reduces to

\[ \begin{equation} \int_{x_\text{1}}^{x_\text{2}} u_\text{j}(x)\left[-\frac{d}{dx}\left(a\frac{d{\varphi^{1}_\text{h}}}{dx}\right) +b\varphi^{1}_\text{h}-f\right] \text{d} x=0,\quad\text{for}\quad j=1,2,...,n. \end{equation}\] (4.6)

The weighted-residual approach provides residual minimization by integrating over the element. Using the integration by parts [148], the order of differentiation in φh1 is reduced as follows

\[\begin{equation} \int_{x_\text{1}}^{x_\text{2}}\left(a\frac{d{u_\text{j}(x)}}{dx}\frac{d{\varphi^{1}_\text{h}}}{dx}+bu_\text{j}(x)\varphi^{1}_\text{h}\right) \text{d} x=\int_{x_\text{1}}^{x_\text{2}}u_\text{j}(x)f \text{d} x+\left[au_\text{j}(x)\frac{d{\varphi^{1}_\text{h}}}{dx}\right]_{x_\text{1}}^{x_\text{2}}. \end{equation}\] (4.7)

In this way, it is necessary that the basis functions uk1 be such that the approximate solution φh1 is only once-differentiable, and not twice differentiable as in equation (4.6). Equation (4.7) represents the weak formulation of the problem. The term "weak" refers to the reduced or weakened continuity of uk1 such that φh1 is differentiable as many times as called for the original DE in the weighted-integral statement [125].

Step 2 - Approximate Solution

The approximate solution φh1 within the finite element Ω1 must be found in a way that the approximate solution of the problem φh is convergent to the actual solution φ as the number of elements increase. For this purpose it is required to satisfy the following conditions [125]:

Consider the element Ω1 composed of two nodes at x1 and x2, as depicted in (4.3).
Figure 4.3: Schematic view of the linear finite element approximation. The true solution φ(x) (dotted line) and the construction of the approximate solution φh1 (solid line), within the element Ω1, are presented.
femappchap4

Since only the first derivatives of φh1 are involved in equation (4.7), the solution φh1 within the element Ω1 is approximated using a linear polynomial as follows

\[\begin{equation} \varphi^{1}_\text{h}(x)=g^{1}_\text{1}+g^{1}_\text{2}x\quad\quad x_\text{1}\leq x\leq x_\text{2}, \end{equation}\] (4.8)

where g11 and g21 are unknown coefficients. The expression (4.8) fulfills the first two conditions of an approximation described above. The third condition is satisfied if the unknown coefficients g11 and g21 are expressed in terms of nodal values of the solution φ11(x1) and φ21(x2). The nodal values are given by

\[\begin{align} \varphi^{1}_\text{1}(x_\text{1})=c^{1}_\text{1}=g^{1}_\text{1}+g^{1}_\text{2}(x_\text{1}),\\ \varphi^{1}_\text{2}(x_\text{2})=c^{1}_\text{2}=g^{1}_\text{1}+g^{1}_\text{2}(x_\text{2}). \end{align}\] (4.9)

By solving the system of equations (4.9) and using h1=x2-x1, the unknown coefficients can be obtained as follows

\[\begin{align} g^{1}_\text{1}=\cfrac{c^{1}_\text{1}x_\text{2}-c^{1}_\text{2}x_\text{1}}{h_\text{1}},\\ g^{1}_\text{2}=\cfrac{c^{1}_\text{2}-c^{1}_\text{1}}{h_\text{1}}. \end{align}\] (4.10)

The substitution of equations (4.10) into equation (4.8) permits one to approximate the solution within the element as follows

\[\begin{align} \varphi^{1}_\text{h}(x)=\cfrac{c^{1}_\text{1}x_\text{2}-c^{1}_\text{2}x_\text{1}}{h_\text{1}}+\cfrac{c^{1}_\text{2}-c^{1}_\text{1}}{h_\text{1}}=\cfrac{x_\text{2}-x_\text{1}}{h_\text{1}}c^{1}_\text{1}+\cfrac{x-x_\text{1}}{h_\text{1}}c^{1}_\text{2}=u^{1}_\text{1}(x)c^{1}_\text{1}+u^{1}_\text{2}(x)c^{1}_\text{2}=\sum_{k=1}^{2}u^{1}_\text{k}(x)c^{1}_\text{k}, \end{align}\] (4.11)

where u11(x) and u21(x) are the basis functions or interpolation functions [116].

The interpolation property of these functions requires that uk1(x) is unity at the kth node and zero at the other node of the given element. Another property of the interpolation functions is that their sum is unity. In this way, the solution φh1(x) is interpolated using its nodal values c11 and c21, is approximated by a piecewise linear polynomial, and its gradient is constant within the element. The geometric representation of the approximate solution inside the element Ω1 is shown in (4.3).

Step 3 - Ritz-Galerkin Method

Once the form of the approximate solution φh1 within the element Ω1 is known, the application of the Ritz-Galerkin method allows one to derive the set of algebraic equations over the finite element. The substitution of equation (4.11) into equation (4.7) will give the set of equations among the nodal values of the element as follows

\[\begin{equation} \sum_{k=1}^{2}\left[\int_{x_\text{1}}^{x_\text{2}}\left(a\frac{d{u_\text{j}(x)}}{dx}\frac{d{u^{1}_\text{k}}}{dx}+bu_\text{j}(x)u^{1}_\text{k}\right) \text{d}x\right] c^{1}_\text{k}=\int_{x_\text{1}}^{x_\text{2}}u_\text{j}(x)f \text{d}x+\sum_{k=1}^{2}u_\text{j}(x_\text{k})\beta^{1}_\text{k}, \end{equation}\] (4.12)

where

\[\begin{equation} \beta^{1}_\text{k}=\left(a\frac{du^{1}_\text{k}}{dx}\right)\Bigg|_{x=x_\text{k}} \end{equation}\] (4.13)

are the unknown source terms. Usually, such a relation for the finite element Ω1 is presented as a matrix equation [161] of the form

\[\begin{equation} [K^{1}]\{c^{1}\}=\{f^{1}\}+\{\beta^{1}\}, \end{equation}\] (4.14)

where

\[\begin{equation} [K^{1}]=\int_{x_\text{1}}^{x_\text{2}}\left(a\frac{du_\text{j}(x)}{dx}\frac{du^{1}_\text{k}}{dx}+bu_\text{j}(x)u^{1}_\text{k}\right)\text{d}x, \end{equation}\] (4.15)

\[\begin{equation} \{f^{1}\}=\int_{x_\text{1}}^{x_\text{2}}u_\text{j}(x)f dx, \end{equation}\] (4.16)

and

\[\begin{equation} \{\beta^{1}\}=\sum_{k=1}^{2}u_\text{j}(x_\text{k})\beta^{1}_\text{k}. \end{equation}\] (4.17)

The matrix [K1], which is called coefficient matrix, and the vector {f1}, which is the source vector, can be determined for the given element and known functions a, b, and f. The relation (4.14) provides two equations for the finite element used to determine four unknowns: c11, c21, β11, and β21. For a linear element, the matrix form is written as follows

\[\begin{equation} \begin{bmatrix} K^{1}_\text{11} & K^{1}_\text{12} \\[2pt] K^{1}_\text{21} & K^{1}_\text{22} \\ \end{bmatrix} \begin{Bmatrix} c^{1}_\text{1} \\[2pt] c^{1}_\text{2} \\ \end{Bmatrix} = \begin{Bmatrix} f^{1}_\text{1} \\[2pt] f^{1}_\text{2} \\ \end{Bmatrix} +\begin{Bmatrix} \beta^{1}_\text{1} \\[2pt] \beta^{1}_\text{2} \\ \end{Bmatrix}, \end{equation}\] (4.18)

where the coefficients of the matrix (Kjk1) and the source vector (fj1) can be easily evaluated for the given element from equations (4.15) and (4.16), respectively, by using the values of the constants a, b, and f.




M. Rovitto: Electromigration Reliability Issue in Interconnects for Three-Dimensional Integration Technologies