Evaluating the surface integral in an FEM (Finite Elements Method) procedure

Question

I want to evaluate the Surface force integral in an FEM procedure.

The basic reference tet is shown in the figure. The faces are numbered corresponding to the node opposite to them. For example the nodes 2-3-4 which forms a face opposite to node 1 is labeled as 1, and so on... I now transform the tet into the volume coordinate denoted by $(\xi_1, \xi_2, \xi_3, \xi_4 )$ which are also the shape functions of the linear element. The surface force vector on a face is denoted by $f_e$. I now need to evaluate the integral $$\int_{\Gamma_e} N^T f_e \ d\Gamma_e $$ where $\Gamma_e$ represent any face. $N^T$ represent the shape function matrix given by $$\begin{bmatrix} \xi_1 & 0 & 0 & \cdots & \xi_4 & 0 & 0 \\ 0 & \xi_1 & 0 & \cdots & 0 & \xi_4 & 0 \\ 0 & 0 & \xi_1 & \cdots & 0 & 0 & \xi_4 \end{bmatrix} $$ To evaluate the surface integral one takes the trace of the matrix on the particular face. For example on face one, $\xi_1=0$. I now need to use a Gaussian Quadrature and not sure how to evaluate the transformed Jacobians. I am actually implementing the above procedure on a 10 node element. So actually confused on how to implement the surface quadratures. Another source of doubt is that in case of the linear element one can do a node lumping i,e distributing the surface load equally among all the nodes and that gives a very simple expression for the surface vector. Can this be also be done for a higher element?

Answer 1

The particular surface integral you want to calculate is basically a specific case of integrating a function over a surface defined in terms of two parametric coordinates. Lets first consider this general case. A fundamental relation from differential geometry is (using your notation)

$$ d{\bf \Gamma_e} = \left\{\begin{array}{c} \frac{\partial x}{\partial\xi_1} \\ \frac{\partial y}{\partial\xi_1} \\ \frac{\partial z}{\partial\xi_1} \end{array}\right\} \times \left\{\begin{array}{c} \frac{\partial x}{\partial\xi_2} \\ \frac{\partial y}{\partial\xi_2} \\ \frac{\partial z}{\partial\xi_2} \end{array}\right\} $$ where $\xi_1$ and $\xi_2$ are the two parametric coordinates of the surface. An integral of a function, $f$ over the surface would be written as

$$\int f d\Gamma_e = \int f(\xi_1, \xi_2) ||d{\bf \Gamma_e}||_2 d\xi_1 d\xi_2$$ where $||d{\bf \Gamma_e}||_2$ is just the magnitude of the vector.

Now lets consider the specific case of faces of the ten-node tetrahedron. In an isoparametric formulation, $x, y,$ and $z$ are written in terms of the shape functions, $N_I$, e.g. $$ x = N_I(\xi_1, \xi_2) x_I $$ where $ x_I$ are the values of $x$ at the nodes of the triangular face. The expressions for $N_I$ on a face can be obtained from the element shape functions by specializing them for a particular face, or equivalently, by using the shape functions for the six-node triangle. From these equations, you can easily evaluate the vector cross product above.

The most general way to actually perform the integration is by using the appropriate Gauss-Legendre rule for the triangle.

Finally, you asked about the difference between formally evaluating the consistent load integral and simply lumping equal values at the nodes. The equivalent nodal loads are quite different from the simple lumped values. For a simple, straight-sided element, the equivalent nodal loads are non-zero only at the mid-side nodes. If you have a fine enough mesh, the simple lumped approach will give you the correct answer but it is not optimum in the Galerkin-FEM sense.

Answer 2

You may please refer chapter two of Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations by M. Asghar Bhatti.

Section 2.4.5 of the book deals with the evaluation of surface integrals. A workout example of evaluating the surface pressure loads with Mathematica code is also given in the book. A Matlab code is also given in the book site for 8 node brick element.