3.1.1 Evaluation of Element Matrices

To solve complex differential equation systems it is necessary to divide the domain of interest into a number of subdomains or elements. Once the concept of subdividing a physical domain into elements is accepted the problem of analyzing a variation of a variable across the whole domain becomes far easier, since the variation can be related with each single element. If such an approximation across an element is defined and since elements are interconnected, it is a simple matter that the variation across the whole domain has, in a piecewise manner, been effected.

The shape functions are a means of interpolation of the discrete nodal variables. In general they are set up in a local coordinate system ( $\xi$,$\eta$,$\zeta$) which implies an element transformation from global coordinates to local ones. In comparison to other discretization methods it seems to be an overhead to transform a general element located in the (x,y,z)-system to a standard element defined in a local coordinate system ( $\xi$,$\eta$,$\zeta$). At closer inspection this transformation turns out to enable necessary integrations over an element to be performed with ease. Furthermore, a generalized methodology can be defined offering the possibility to solve complex PDE systems on various elements with different interpolation orders and dimensions without necessity to change the underlying concept.

To perform finite element analysis the matrices defining element properties have to be found. They will always be of the form

$$\int_{\Omega}^{} \Omega$$ (3.6)

or in the example of the Laplace operator (3.2)

$$\int_{\Omega}^{} \Delta \psi \Omega \int_{\Omega}^{} \nabla^{2}_{} \psi \Omega \int_{\Omega}^{} \left(\vphantom{\frac{d^2\psi}{dx^2} +\frac{d^2\psi}{dy^2} +\frac{d^2\psi}{dz^2}}\right. {\frac{d^2\psi}{dx^2}} {\frac{d^2\psi}{dy^2}} {\frac{d^2\psi}{dz^2}} \left.\vphantom{\frac{d^2\psi}{dx^2} +\frac{d^2\psi}{dy^2} +\frac{d^2\psi}{dz^2}}\right) \Omega$$ (3.7)

Since the values of $\psi$ are only defined on discrete points of an element a variation of the function has to be chosen. This is normally done using (3.3) and introducing N as a shape function of choice (see Appendix B).

With these definitions (3.7) can be written as

$$\int_{\Omega}^{} \nabla^{2}_{} \psi \Omega \sum_{i=1}^{n} \left(\vphantom{ \int _\Omega w\cdot (\nabla ^2 N_i \cdot \psi_i) }\right. \int_{\Omega}^{} \nabla^{2}_{} \psi_{i}^{} \left.\vphantom{ \int _\Omega w\cdot (\nabla ^2 N_i \cdot \psi_i) }\right) \Omega$$ (3.8)

Galerkin's approach suggests to use weighting functions w that are set to be equal to the shape functions N (w_j = N_j with j = 1,2, ... ,m). Because N only depends on the global coordinates (x,y,z) (3.7) can now be transformed to

$$\int_{\Omega}^{} \nabla^{2}_{} \psi \Omega \sum_{j=1}^{m} \sum_{i=1}^{n} \psi_{i}^{} \left(\vphantom{ \int _\Omega N_j \cdot \nabla ^2 N_i }\right. \int_{\Omega}^{} \nabla^{2}_{} \left.\vphantom{ \int _\Omega N_j \cdot \nabla ^2 N_i }\right) \Omega$$ (3.9)

The resulting integral formulation has partial differential operators of second order, which can be eliminated using Green's formula

$$\int_{\Omega}^{} \nabla^{2}_{} \Omega \int_{\Omega}^{} \nabla \nabla \Omega \oint_{\Gamma}^{} {\frac{\partial v}{\partial n}} \Gamma$$ (3.10)

resulting in

$$\int_{\Omega}^{} \nabla^{2}_{} \psi \Omega \sum_{j=1}^{m} \sum_{i=1}^{n} \psi_{i}^{} \left(\vphantom{ - \int _\Omega \nabla N_j \cdot \nabla N_i + \oint _\Gamma N_j \cdot \frac{\partial N_i}{\partial n} \;d\Gamma }\right. \int_{\Omega}^{} \nabla \nabla \oint_{\Gamma}^{} {\frac{\partial N_i}{\partial n}} \Gamma \left.\vphantom{ - \int _\Omega \nabla N_j \cdot \nabla N_i + \oint _\Gamma N_j \cdot \frac{\partial N_i}{\partial n} \;d\Gamma }\right)$$ (3.11)

It is remarkable that two transformations are necessary to evaluate this integral over the domain of interest:

• because the shape functions N_i are defined in terms of local (even curvilinear) coordinates ( $\xi$,$\eta$,$\zeta$), it is necessary to devise some means of expressing the global derivatives of the type occurring in (3.11) by local derivatives.
• the element of volume (or surface in case of boundaries) over which the integration has to be carried out needs to be expressed in terms of the local coordinates with an appropriate change of integration limits.

To express these transformations, consider a set of local coordinates ( $\xi$,$\eta$,$\zeta$) and a corresponding set of global coordinates (x, y, z). By the rules of partial differential equations the derivatives of the shape functions have the form

$${\frac{\partial N_i}{\partial \xi}} {\frac{\partial N_i}{\partial x}} {\frac{\partial x}{\partial \xi}} {\frac{\partial N_i}{\partial y}} {\frac{\partial y}{\partial \xi}} {\frac{\partial N_i}{\partial z}} {\frac{\partial z}{\partial \xi}}$$ (3.12)

$${\frac{\partial N_i}{\partial \eta}} {\frac{\partial N_i}{\partial x}} {\frac{\partial x}{\partial \eta}} {\frac{\partial N_i}{\partial y}} {\frac{\partial y}{\partial \eta}} {\frac{\partial N_i}{\partial z}} {\frac{\partial z}{\partial \eta}}$$ (3.13)

$${\frac{\partial N_i}{\partial \zeta}} {\frac{\partial N_i}{\partial x}} {\frac{\partial x}{\partial \zeta}} {\frac{\partial N_i}{\partial y}} {\frac{\partial y}{\partial \zeta}} {\frac{\partial N_i}{\partial z}} {\frac{\partial z}{\partial \zeta}}$$ (3.14)

or in matrix notation:

$$\left(\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial \xi... ...N_i}{\partial \eta}\\ \frac{\partial N_i}{\partial \zeta} \end{array} }\right. \begin{array}{c} \frac{\partial N_i}{\partial \xi}\\ \frac{\partial N_i}{\partial \eta}\\ \frac{\partial N_i}{\partial \zeta} \end{array} \left.\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial \xi... ...N_i}{\partial \eta}\\ \frac{\partial N_i}{\partial \zeta} \end{array} }\right) \left(\vphantom{\begin{array}{ccc} \frac{\partial x}{\partial \xi... ...ial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta} \end{array} }\right. \begin{array}{ccc} \frac{\partial x}{\partial \xi} & \frac{\parti... ...rac{\partial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta} \end{array} \left.\vphantom{\begin{array}{ccc} \frac{\partial x}{\partial \xi... ...ial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta} \end{array} }\right) \left(\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial x}\... ...artial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z} \end{array} }\right. \begin{array}{c} \frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z} \end{array} \left.\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial x}\... ...artial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z} \end{array} }\right) \bf J \left(\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial x}\... ...artial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z} \end{array} }\right. \begin{array}{c} \frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z} \end{array} \left.\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial x}\... ...artial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z} \end{array} }\right)$$ (3.15)

The left-hand side of (3.15) depends only on the shape functions N_i that are specified in the local coordinate system. Furthermore, as (x, y, z) are given explicitly, the matrix J can be found in terms of local coordinates. This matrix is known as the Jacobian matrix.

Since the point of interest is the definition of derivatives in terms of locally defined shape functions the remaining task is to convert the Jacobian matrix J. As a result the derivatives of the shape functions in the general coordinate system are found to be:

$$\left(\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial x}\... ...artial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z} \end{array} }\right. \begin{array}{c} \frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z} \end{array} \left.\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial x}\... ...artial N_i}{\partial y}\\ \frac{\partial N_i}{\partial z} \end{array} }\right) \bf J^{-1}_{} \left(\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial \xi... ...N_i}{\partial \eta}\\ \frac{\partial N_i}{\partial \zeta} \end{array} }\right. \begin{array}{c} \frac{\partial N_i}{\partial \xi}\\ \frac{\partial N_i}{\partial \eta}\\ \frac{\partial N_i}{\partial \zeta} \end{array} \left.\vphantom{\begin{array}{c} \frac{\partial N_i}{\partial \xi... ...N_i}{\partial \eta}\\ \frac{\partial N_i}{\partial \zeta} \end{array} }\right)$$ (3.16)

Figure 3.2: Transformation of elements onto a local standard element within the limits zero and one

For integration purposes a general located element defined in the spatial coordinates (x, y, z) is mapped onto a standard element in the ( $\xi$,$\eta$,$\zeta$) coordinate system (Fig. 3.2). Again the variation approach can be used but now interpolating the coordinate vectors itself instead of a distributed value. The variation approach can be written as

$$\bf x^T \cdot N(\xi,\eta,\zeta)$$

$$\bf y^T \cdot N(\xi,\eta,\zeta)$$ (3.17)

$$\bf z^T \cdot N(\xi,\eta,\zeta)$$

where x,y and z are the coordinate vectors of the element points. Notice that (3.17) is the same formulation as (3.3) only applied to the global coordinate values. Now the Jacobian matrix can be written using a determinant style

$$\bf J \left\vert\vphantom{\frac{\partial(x,y,z)}{\partial(\xi,\eta,\zeta)}}\right. {\frac{\partial(x,y,z)}{\partial(\xi,\eta,\zeta)}} \left.\vphantom{\frac{\partial(x,y,z)}{\partial(\xi,\eta,\zeta)}}\right\vert \left\vert\vphantom{\begin{array}{ccc} \frac{\partial x}{\partial... ...ial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta} \end{array} }\right. \begin{array}{ccc} \frac{\partial x}{\partial \xi} & \frac{\parti... ...rac{\partial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta} \end{array} \left.\vphantom{\begin{array}{ccc} \frac{\partial x}{\partial \xi... ...y}{\partial \zeta} & \frac{\partial z}{\partial \zeta} \end{array} }\right\vert \left\vert\vphantom{ \begin{array}{ccc} {\bf x}^T \frac{\partial ... ...zeta}&{\bf z}^T \frac{\partial {\bf N}}{\partial \zeta}\\ \end{array} }\right. \begin{array}{ccc} {\bf x}^T \frac{\partial {\bf N}}{\partial \xi... ...partial \zeta}&{\bf z}^T \frac{\partial {\bf N}}{\partial \zeta}\\ \end{array} \left.\vphantom{ \begin{array}{ccc} {\bf x}^T \frac{\partial {\bf... ...}&{\bf z}^T \frac{\partial {\bf N}}{\partial \zeta}\\ \end{array} }\right\vert$$ (3.18)

and its inversion

$$\bf J^{-1}_{} \left\vert\vphantom{\frac{\partial(\xi,\eta,\zeta)}{\partial(x,y,z)}}\right. {\frac{\partial(\xi,\eta,\zeta)}{\partial(x,y,z)}} \left.\vphantom{\frac{\partial(\xi,\eta,\zeta)}{\partial(x,y,z)}}\right\vert \left\vert\vphantom{\begin{array}{ccc} \frac{\partial \xi}{\parti... ...tial \eta}{\partial z} & \frac{\partial \zeta}{\partial z} \end{array} }\right. \begin{array}{ccc} \frac{\partial \xi}{\partial x} & \frac{\parti... ...frac{\partial \eta}{\partial z} & \frac{\partial \zeta}{\partial z} \end{array} \left.\vphantom{\begin{array}{ccc} \frac{\partial \xi}{\partial x... ... \eta}{\partial z} & \frac{\partial \zeta}{\partial z} \end{array} }\right\vert$$ (3.19)

With the help of parameter comparisons and simple transformations the global derivatives can be extracted.

To transform the variables and the region with respect to which the integration is made a standard process will be used that also involves the determinant J. Thus a volume element becomes

$$\bf J \xi \eta \zeta$$ (3.20)

Assuming that the inverse of J can be found the problem of general evaluation of an integral (see also Section 3.1.2) is now reduced to element properties $\bf F$ of the form

$$\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \bf F \xi \eta \zeta \xi \eta \zeta$$ (3.21)

The integration is carried out within a standard element in the local coordinate system and not in the complicated distorted shape, thus accounting for the simple integration limits. Similar one- and two-dimensional problems will result in integrals with one or two coordinates within simple limits, since the complete concept is valid independent of the dimension used.

While the limits of the integration are simple, unfortunately the explicit form of $\bf F$($\xi$,$\eta$,$\zeta$) is not. Non standard algebraic integration usually defies the mathematical average skill, and numerical integration has to be resorted to. This, as will be seen in the next section, is not a severe penalty and has the advantage that the solution is not tied to a particular element or partial differential equation but can be written for various classes of problems.