Search Continuum Mechanics Website

The usual FE notation will be used. In particular, \({\bf u} = (u,v,w)\) where \(u, v\), and \(w\) are the \(X, Y\), and \(Z\) components of the displacement of any given node.

As usual, we'll start with 2-D examples to explain the concepts before going to 3-D.

\[ \begin{eqnarray} u(X,Y) & = & {1\over 4}(1-X)(1-Y)u_1 + {1\over 4}(1+X)(1-Y)u_2 + {1\over 4}(1+X)(1+Y)u_3 + {1\over 4}(1-X)(1+Y)u_4 \\ \\ \\ v(X,Y) & = & {1\over 4}(1-X)(1-Y)v_1 + {1\over 4}(1+X)(1-Y)v_2 + {1\over 4}(1+X)(1+Y)v_3 + {1\over 4}(1-X)(1+Y)v_4 \\ \end{eqnarray} \]

where \(u(X,Y)\) is the displacement in the x-direction as a function of \(X\) and \(Y\), \(v(X,Y)\) is the displacement in the y-direction, and \(u_1, v_1, u_2, v_2\) etc., are the \(x\) and \(y\) displacements at nodes 1, 2, etc.

The equations are often written in shorthand as

\[ \begin{eqnarray} u(X,Y) & = & \phi_1(X,Y) u_1 + \phi_2(X,Y) u_2 + \phi_3(X,Y) u_3 + \phi_4(X,Y) u_4 \\ \\ v(X,Y) & = & \phi_1(X,Y) v_1 + \phi_2(X,Y) v_2 + \phi_3(X,Y) v_3 + \phi_4(X,Y) v_4 \end{eqnarray} \]

where \(\phi_1(X,Y) = (1 - X)(1 - Y)/4\) and \(\phi_2(X,Y) = (1 + X)(1 - Y)/4\), etc. These are called

This figure shows Shape Function #4, for example. This shape function is \(\phi_4(X,Y) = (1 - X)(1 + Y)/4\). It equals 1 at (-1,1) and zero at all other nodes.

Evaluating these equations for several special cases can help to better understand them. For starters, evaluate them at the center of the element: \(X=0, Y=0\).

\[ \begin{eqnarray} u(0,0) & = & { u_1 + u_2 + u_3 + u_4 \over 4 } \\ \\ \\ v(0,0) & = & { v_1 + v_2 + v_3 + v_4 \over 4 } \end{eqnarray} \]

Since the center is equidistant from all four nodes, each one gets equal weighting.

\[ \begin{eqnarray} u(0,-1) & = & { u_1 + u_2 \over 2 } \\ \\ \\ v(0,-1) & = & { v_1 + v_2 \over 2 } \end{eqnarray} \]

And setting \(X = 0.8\) and \(Y = 0.8\) calculates the displacement close to node #3.

\[ \begin{eqnarray} u(0.8,0.8) & = & 0.01 u_1 + 0.09 u_2 + 0.81 u_3 + 0.09 u_4 \\ \\ \\ v(0.8,0.8) & = & 0.01 v_1 + 0.09 v_2 + 0.81 v_3 + 0.09 v_4 \end{eqnarray} \]

The closest node gets the greatest weighting and the farthest node gets the least weighting.

\[ 0.01 + 0.09 + 0.81 + 0.09 = 1 \]

\[ \begin{eqnarray} u_1 = 3 \qquad u_2 = 2 \qquad u_3 = 3 \qquad u_4 = 2 \\ v_1 = 2 \qquad v_2 = 2 \qquad v_3 = 2 \qquad v_4 = 2 \end{eqnarray} \]

So the displacement of a point at \({\bf X} = (0.75,0.50)\) is

\[ \begin{eqnarray} u(0.75,0.50) & = & {1\over 4}[(1\text{-}0.75)(1\text{-}0.50)u_1 \text{+} (1\text{+}0.75)(1\text{-}0.50)u_2 \text{+} (1\text{+}0.75)(1\text{+}0.50)u_3 \text{+} (1\text{-}0.75)(1\text{+}0.50)u_4] \\ \\ \\ v(0.75,0.50) & = & {1\over 4}[(1\text{-}0.75)(1\text{-}0.50)v_1 \text{+} (1\text{+}0.75)(1\text{-}0.50)v_2 \text{+} (1\text{+}0.75)(1\text{+}0.50)v_3 \text{+} (1\text{-}0.75)(1\text{+}0.50)v_4] \\ \end{eqnarray} \]

which simplifies to

\[ \begin{eqnarray} u(0.75,0.50) & = & 0.0313 u_1 + 0.2188 u_2 + 0.6563 u_3 + 0.0938 u_4 \\ \\ v(0.75,0.50) & = & 0.0313 v_1 + 0.2188 v_2 + 0.6563 v_3 + 0.0938 v_4 \end{eqnarray} \]

It's clear that the displacement of node #3 gets the highest weighting because the coordinates \((0.75, 0.50)\) are closest to it.

Substituting in the u's and v's gives

\[ \begin{eqnarray} u(0.75,0.50) & = & 0.0313 * 3 + 0.2188 * 2 + 0.6563 * 3 + 0.0938 * 2 \\ \\ & = & 2.688 \\ \\ \\ v(0.75,0.50) & = & 0.0313 * 2 + 0.2188 * 2 + 0.6563 * 2 + 0.0938 * 2 \\ \\ & = & 2 \end{eqnarray} \]

So \({\bf u} = (2.688, 2)\). The x-component is closer to 3 than 2 because the point itself is close to Node #3, which had \(u_3 = 3\). The y-component is 2 because all the nodes move vertically that same amount.

\[ {\bf F} = {\bf I} + {\partial {\bf u} \over \partial {\bf X}} \]

The individual components of \({\bf F}\) are

\[ \begin{eqnarray} F_{11} & = & {1 \over 4} [-(1-Y)u_1 + (1-Y)u_2 + (1+Y)u_3 - (1+Y)u_4] + 1 \\ \\ F_{12} & = & {1 \over 4} [-(1-X)u_1 - (1+X)u_2 + (1+X)u_3 + (1-X)u_4] \\ \\ F_{21} & = & {1 \over 4} [-(1-Y)v_1 + (1-Y)v_2 + (1+Y)v_3 - (1+Y)v_4] \\ \\ F_{22} & = & {1 \over 4} [-(1-X)v_1 - (1+X)v_2 + (1+X)v_3 + (1-X)v_4] + 1 \\ \end{eqnarray} \]

Note that the deformation gradient varies throughout the element. It is not constant. This is actually a big advantage of quadrilaterals (and bricks in 3-D) over linear triangles (and tetrahedrals in 3-D). The deformation gradient in a linear triangle or tetrahedral does not vary from location to location. This tends to make them respond with artificially (and incorrectly) elevated stiffnesses.

\[ \begin{eqnarray} u_1 = 3 \qquad u_2 = 2 \qquad u_3 = 3 \qquad u_4 = 2 \\ v_1 = 2 \qquad v_2 = 2 \qquad v_3 = 2 \qquad v_4 = 2 \end{eqnarray} \]

The deformation gradient at point \({\bf X} = (0.75,0.50)\) is

\[ \begin{eqnarray} F_{11} & = & {1 \over 4} [-(1-Y)u_1 + (1-Y)u_2 + (1+Y)u_3 - (1+Y)u_4] + 1 \\ \\ & = & 0.25 [ -(0.50)(3) + (0.50)(2) + (1.50)(3) - (1.50)(2) ] + 1\\ \\ & = & 1.250 \\ \\ \\ F_{12} & = & {1 \over 4} [-(1-X)u_1 - (1+X)u_2 + (1+X)u_3 + (1-X)u_4] \\ \\ & = & 0.25 [ -(0.25)(3) - (1.75)(2) + (1.75)(3) + (0.25)(2) ] \\ \\ & = & 0.375 \\ \\ \\ F_{21} & = & {1 \over 4} [-(1-Y)v_1 + (1-Y)v_2 + (1+Y)v_3 - (1+Y)v_4] \\ \\ & = & 0.25 [ -(0.50)(2) + (0.50)(2) + (1.50)(2) - (1.50)(2) ] \\ \\ & = & 0.000 \\ \\ \\ F_{22} & = & {1 \over 4} [-(1-X)v_1 - (1+X)v_2 + (1+X)v_3 + (1-X)v_4] + 1 \\ \\ & = & 0.25 [ -(0.25)(2) - (1.75)(2) + (1.75)(2) + (0.25)(2) ] + 1\\ \\ & = & 1.000 \end{eqnarray} \]

So \({\bf F}\) equals

\[ {\bf F} = \left[ \matrix{ 1.250 & 0.375 \\ 0.000 & 1.000 } \right] \]

This makes sense because the upper right corner of the element is being stretched horizontally (\(F_{11} = 1.25\)) and is also shearing (\(F_{12} = 0.375\)).

\[ \begin{eqnarray} X_1 = -1 & \quad & X_2 = 1 & \quad & X_3 = 0 & \quad & X_4 = -1 \\ Y_1 = -1 & \quad & \, Y_2 = 0 & \quad & \, Y_3 = 2 & \quad & \, Y_4 = \;\;\;1 \\ \end{eqnarray} \]

Since the element is not aligned with the global coordinate system, the simple approach from above will not work. Instead, it is necessary to introduce two new local coordinates, \(r\) and \(s\), that span from -1 to +1.

The position within the element can be described in terms of the node coordinates and \(r\) and \(s\) as

\[ \begin{eqnarray} X(r,s) & = & {1\over 4}(1-r)(1-s)X_1 + {1\over 4}(1+r)(1-s)X_2 + {1\over 4}(1+r)(1+s)X_3 + {1\over 4}(1-r)(1+s)X_4 \\ \\ \\ Y(r,s) & = & {1\over 4}(1-r)(1-s)Y_1 + {1\over 4}(1+r)(1-s)Y_2 + {1\over 4}(1+r)(1+s)Y_3 + {1\over 4}(1-r)(1+s)Y_4 \\ \end{eqnarray} \]

These are the same shape functions as before, only in \(r\) and \(s\) this time. So these equations can also be written as follows

\[ \begin{eqnarray} X(r,s) & = & \phi_1(r,s) X_1 + \phi_2(r,s) X_2 + \phi_3(r,s) X_3 + \phi_4(r,s) X_4 \\ \\ Y(r,s) & = & \phi_1(r,s) Y_1 + \phi_2(r,s) Y_2 + \phi_3(r,s) Y_3 + \phi_4(r,s) Y_4 \end{eqnarray} \]

where \(\phi_1(r,s) = (1 - r)(1 - s)/4\) and \(\phi_2(r,s) = (1 + r)(1 - s)/4\), etc. Once again, each shape function equals 1 at the location of its node, and zero at all other node locations.

And the displacements in this example are

\[ \begin{eqnarray} u_1 = 3 & \quad & u_2 = 2 & \quad & u_3 = 4 & \quad & u_4 = 3 \\ v_1 = 2 & \quad & \, v_2 = 1 & \quad & \, v_3 = 1 & \quad & \, v_4 = 2 \end{eqnarray} \]

Just as before, the displacements within the element are interpolated using the same shape functions, but with \(r\) and \(s\).

\[ \begin{eqnarray} u(r,s) & = & {1\over 4}(1-r)(1-s)u_1 + {1\over 4}(1+r)(1-s)u_2 + {1\over 4}(1+r)(1+s)u_3 + {1\over 4}(1-r)(1+s)u_4 \\ \\ \\ v(r,s) & = & {1\over 4}(1-r)(1-s)v_1 + {1\over 4}(1+r)(1-s)v_2 + {1\over 4}(1+r)(1+s)v_3 + {1\over 4}(1-r)(1+s)v_4 \\ \end{eqnarray} \]

\[ \begin{eqnarray} X(r,s) & = & \phi_1(r,s) X_1 + \phi_2(r,s) X_2 + \phi_3(r,s) X_3 + \phi_4(r,s) X_4 \\ \\ Y(r,s) & = & \phi_1(r,s) Y_1 + \phi_2(r,s) Y_2 + \phi_3(r,s) Y_3 + \phi_4(r,s) Y_4 \end{eqnarray} \]

and the displacements \((u,v)\) can be calculated using the exact same shape functions

\[ \begin{eqnarray} u(r,s) & = & \phi_1(r,s) u_1 + \phi_2(r,s) u_2 + \phi_3(r,s) u_3 + \phi_4(r,s) u_4 \\ \\ v(r,s) & = & \phi_1(r,s) v_1 + \phi_2(r,s) v_2 + \phi_3(r,s) v_3 + \phi_4(r,s) v_4 \end{eqnarray} \]

The fact that the same shape functions are used for both coordinates and displacements is what makes this an

\[ \begin{eqnarray} {\partial u \over \partial X} = {\partial u \over \partial r} {\partial r \over \partial X} + {\partial u \over \partial s} {\partial s \over \partial X} \\ \\ \\ {\partial u \over \partial Y} = {\partial u \over \partial r} {\partial r \over \partial Y} + {\partial u \over \partial s} {\partial s \over \partial Y} \end{eqnarray} \]

and similar equations exist for \({\partial v \over \partial X}\) and \({\partial v \over \partial Y}\).

The complete set of equations can be written in matrix form as

\[ \begin{eqnarray} \left[ \matrix{ {\partial u \over \partial X} & {\partial u \over \partial Y} \\ \\ {\partial v \over \partial X} & {\partial v \over \partial Y} } \right] = \left[ \matrix{ {\partial u \over \partial r} & {\partial u \over \partial s} \\ \\ {\partial v \over \partial r} & {\partial v \over \partial s} } \right] \left[ \matrix{ {\partial r \over \partial X} & {\partial r \over \partial Y} \\ \\ {\partial s \over \partial X} & {\partial s \over \partial Y} } \right] \end{eqnarray} \]

The \({\partial u \over \partial r}\) and \({\partial u \over \partial s}\) terms, and the others in the first matrix are easy. They are

\[ \begin{eqnarray} {\partial u \over \partial r} & = & {1 \over 4} [-(1-s)u_1 + (1-s)u_2 + (1+s)u_3 - (1+s)u_4] \\ \\ {\partial u \over \partial s} & = & {1 \over 4} [-(1-r)u_1 - (1+r)u_2 + (1+r)u_3 + (1-r)u_4] \end{eqnarray} \]

However, the \({\partial r \over \partial X}\) and \({\partial s \over \partial X}\) terms in the second matrix are trickier because \(r\) and \(s\) are not explicitly functions of \(X\) and \(Y\). It's the other way around. Fortunately, the needed terms are just the inverse of a certain matrix. It is

\[ \begin{eqnarray} \left[ \matrix{ {\partial r \over \partial X} & {\partial r \over \partial Y} \\ \\ {\partial s \over \partial X} & {\partial s \over \partial Y} } \right] = \left[ \matrix{ {\partial X \over \partial r} & {\partial X \over \partial s} \\ \\ {\partial Y \over \partial r} & {\partial Y \over \partial s} } \right]^{-1} \end{eqnarray} \]

And this matrix is easy to compute and then take the numerical inverse of. So the entire process can be written as

\[ \begin{eqnarray} \left[ \matrix{ {\partial u \over \partial X} & {\partial u \over \partial Y} \\ \\ {\partial v \over \partial X} & {\partial v \over \partial Y} } \right] = \left[ \matrix{ {\partial u \over \partial r} & {\partial u \over \partial s} \\ \\ {\partial v \over \partial r} & {\partial v \over \partial s} } \right] \left[ \matrix{ {\partial X \over \partial r} & {\partial X \over \partial s} \\ \\ {\partial Y \over \partial r} & {\partial Y \over \partial s} } \right]^{-1} \end{eqnarray} \]

\[ \begin{eqnarray} X_1 = -1 & \quad & X_2 = 1 & \quad & X_3 = 0 & \quad & X_4 = -1 \\ Y_1 = -1 & \quad & \, Y_2 = 0 & \quad & \, Y_3 = 2 & \quad & \, Y_4 = \;\;\;1 \\ \end{eqnarray} \]

and the displacements are

\[ \begin{eqnarray} u_1 = 3 & \quad & u_2 = 2 & \quad & u_3 = 4 & \quad & u_4 = 3 \\ v_1 = 2 & \quad & \, v_2 = 1 & \quad & \, v_3 = 1 & \quad & \, v_4 = 2 \\ \end{eqnarray} \]

Let's compute the deformation gradient at the centroid of the element, where \(r = s = 0\). First, the centroid itself is at

\[ \begin{eqnarray} X_{cg} \; & = & \; {1 \over 4} (X_1 + X_2 + X_3 + X_4) \; & = & \; {1 \over 4} (-1 + 1 + 0 - 1) \; & = & \; -0.25 \\ \\ Y_{cg} \; & = & \; {1 \over 4} (Y_1 \, + Y_2 \, + Y_3 \, + Y_4) \; & = & \; {1 \over 4} (-1 + 0 + 2 + 1) \; & = & \; \;\; 0.50 \end{eqnarray} \]

And we can calculate \(u\) and \(v\) in order to determine where the centroid ends up. This step is not really necessary, but we'll do it anyway.

\[ \begin{eqnarray} u_{cg} \; & = & \; {1 \over 4} (u_1 + u_2 + u_3 + u_4) \; & = & \; {1 \over 4} (3 + 2 + 4 + 3) \; & = & \; 3.0 \\ \\ v_{cg} \; & = & \; {1 \over 4} (v_1 + v_2 + v_3 + v_4) \; & = & \; {1 \over 4} (2 + 1 + 1 + 2) \; & = & \; 1.5 \end{eqnarray} \]

So \(x_{cg} = -0.25 + 3 = 2.75\) and \(y_{cg} = 0.5 + 1.5 = 2.0\).

The next step in calculating \({\bf F}\) is to compute \({\partial u \over \partial r}\), \({\partial u \over \partial s}\), and all the other terms at the centroid.

\[ \begin{eqnarray} \left. {\partial u \over \partial r}\right|_{0,0} \; & = & \; {1 \over 4} (-u_1 + u_2 + u_3 - u_4) \; & = & \; {1 \over 4} (-3 + 2 + 4 - 3) \; & = \;\; 0.0 \\ \\ \\ \left. {\partial u \over \partial s}\right|_{0,0} \; & = & \; {1 \over 4} (-u_1 - u_2 + u_3 + u_4) \; & = & \; {1 \over 4} (-3 - 2 + 4 + 3) \; & = \;\; 0.5 \\ \\ \\ \left. {\partial v \over \partial r}\right|_{0,0} \; & = & \; {1 \over 4} (-v_1 + v_2 + v_3 - v_4) \; & = & \; {1 \over 4} (-2 + 1 + 1 - 2) \; & = -0.5 \\ \\ \\ \left. {\partial v \over \partial s}\right|_{0,0} \; & = & \; {1 \over 4} (-v_1 - v_2 + v_3 + v_4) \; & = & \; {1 \over 4} (-2 - 1 + 1 + 2) \; & = \;\; 0.0 \end{eqnarray} \]

So the matrix of partial derivatives of the displacement field is

\[ \begin{eqnarray} \left[ \matrix{ {\partial u \over \partial r} & {\partial u \over \partial s} \\ \\ {\partial v \over \partial r} & {\partial v \over \partial s} } \right] = \left[ \matrix{ 0 & 0.5 \\ \\ -0.5 & 0 } \right] \end{eqnarray} \]

The partial derivatives of the node coordinates are

\[ \begin{eqnarray} \left. {\partial X \over \partial r}\right|_{0,0} \; & = & \; {1 \over 4} (-X_1 + X_2 + X_3 - X_4) \; & = & \; {1 \over 4} ( 1 + 1 + 0 + 1) \; & = \;\; 0.75 \\ \\ \\ \left. {\partial X \over \partial s}\right|_{0,0} \; & = & \; {1 \over 4} (-X_1 - X_2 + X_3 + X_4) \; & = & \; {1 \over 4} ( 1 - 1 + 0 - 1) \; & = -0.25 \\ \\ \\ \left. {\partial Y \over \partial r}\right|_{0,0} \; & = & \; {1 \over 4} (-Y_1 + Y_2 + Y_3 - Y_4) \; & = & \; {1 \over 4} ( 1 + 0 + 2 - 1) \; & = \;\; 0.50 \\ \\ \\ \left. {\partial Y \over \partial s}\right|_{0,0} \; & = & \; {1 \over 4} (-Y_1 - Y_2 + Y_3 + Y_4) \; & = & \; {1 \over 4} ( 1 - 0 + 2 + 1) \; & = \;\; 1.0 \end{eqnarray} \]

The matrix of partial derivatives of the coordinates is

\[ \begin{eqnarray} \left[ \matrix{ {\partial X \over \partial r} & {\partial X \over \partial s} \\ \\ {\partial Y \over \partial r} & {\partial Y \over \partial s} } \right] = \left[ \matrix{ 0.75 & -0.25 \\ \\ 0.50 & \;\; 1.00 } \right] \end{eqnarray} \]

And the inverse is

\[ \begin{eqnarray} \left[ \matrix{ {\partial X \over \partial r} & {\partial X \over \partial s} \\ \\ {\partial Y \over \partial r} & {\partial Y \over \partial s} } \right]^{-1} = \left[ \matrix{ \;\;\; 1.1429 & 0.2857 \\ \\ -0.5714 & 0.8571 } \right] \end{eqnarray} \]

So

\[ \begin{eqnarray} \left[ \matrix{ {\partial u \over \partial X} & {\partial u \over \partial Y} \\ \\ {\partial v \over \partial X} & {\partial v \over \partial Y} } \right] = \left[ \matrix{ \;\;\;0 & 0.5 \\ \\ -0.5 & 0.0 } \right] \left[ \matrix{ \;\;\; 1.1429 & 0.2857 \\ \\ -0.5714 & 0.8571 } \right] \end{eqnarray} = \left[ \matrix{ -0.2857 & \;\;\;0.4286 \\ \\ -0.5715 & -0.1429 } \right] \]

And adding \({\bf I}\) to get \({\bf F}\) finally gives

\[ {\bf F} = \left[ \matrix{ 0.7143 & 0.4286 \\ \\ -0.5715 & 0.8572 } \right] \]

And a polar decomposition can be applied to find

\[ {\bf R} = \left[ \matrix{ \;\;\; 0.844 & 0.537 \\ \\ -0.537 & 0.844 } \right] \qquad \qquad {\bf U} = \left[ \matrix{ \;\;\; 0.909 & -0.099 \\ \\ -0.099 & \;\;\; 0.953 } \right] \]

So this central portion of the element rotates

\[ \theta \quad = \quad \text{Sin}^{-1}(-0.537) \quad = \quad -32.5^\circ \]

And the horizontal part of the element is compressed ~9%, the vertical part is compressed ~5%, and it is sheared -20%. These results may be intuitive whenever the element is originally aligned with the global coordinate system, however that does not seem to be the case when the element is randomly oriented.

\[ \begin{eqnarray} u(r,s,t) & \; = \; & {1\over 8}(1-r)(1-s)(1-t)u_1 \; + \; {1\over 8}(1+r)(1-s)(1-t)u_2 \; + \\ \\ & & {1\over 8}(1+r)(1+s)(1-t)u_3 \; + \; {1\over 8}(1-r)(1+s)(1-t)u_4 \; + \\ \\ & & {1\over 8}(1-r)(1-s)(1+t)u_5 \; + \; {1\over 8}(1+r)(1-s)(1+t)u_6 \; + \\ \\ & & {1\over 8}(1+r)(1+s)(1+t)u_7 \; + \; {1\over 8}(1-r)(1+s)(1+t)u_8 \end{eqnarray} \]

The equations for the \(v\) and \(w\) displacements are the same, as well as the coordinate interpolation equations for \(X\), \(Y\), and \(Z\).

The shape functions for triangles are

\[ \begin{eqnarray} X(r,s) & \; = \; & r X_1 \; + \; s X_2 \; + \; (1 - r - s) X_3 \\ \\ Y(r,s) & \; = \; & r Y_1 \; + \; s Y_2 \; + \; (1 - r - s) Y_3 \\ \end{eqnarray} \]

for the coordinates and

\[ \begin{eqnarray} u(r,s) & \; = \; & r u_1 \; + \; s u_2 \; + \; (1 - r - s) u_3 \\ \\ v(r,s) & \; = \; & r v_1 \; + \; s v_2 \; + \; (1 - r - s) v_3 \\ \end{eqnarray} \]

for the displacements.

The derivatives are easy in this case. Recall that

\[ \begin{eqnarray} \left[ \matrix{ {\partial u \over \partial X} & {\partial u \over \partial Y} \\ \\ {\partial v \over \partial X} & {\partial v \over \partial Y} } \right] = \left[ \matrix{ {\partial u \over \partial r} & {\partial u \over \partial s} \\ \\ {\partial v \over \partial r} & {\partial v \over \partial s} } \right] \left[ \matrix{ {\partial X \over \partial r} & {\partial X \over \partial s} \\ \\ {\partial Y \over \partial r} & {\partial Y \over \partial s} } \right]^{-1} \end{eqnarray} \]

Applying the derivatives to the triangle shape functions results in

\[ \begin{eqnarray} \left[ \matrix{ {\partial u \over \partial X} & {\partial u \over \partial Y} \\ \\ {\partial v \over \partial X} & {\partial v \over \partial Y} } \right] = \left[ \matrix{ u_1 - u_3 & u_2 - u_3 \\ \\ v_1 - v_3 & v_2 - v_3 \\ } \right] \left[ \matrix{ X_1 - X_3 & X_2 - X_3 \\ \\ Y_1 - Y_3 & Y_2 - Y_3 \\ } \right]^{-1} \end{eqnarray} \]

The key revelation here is that the result is independent of \(r\) and \(s\), which means that it is constant throughout the element. This constancy is the problem, for it tends to make triangles react much too stiff in FE analyses. Because of this, they should be avoided (just the linear ones, not the quadratics, which are fine) whenever possible. Although it is recognized that a few are almost always necessary anyway.

The shape functions for linear tetrahedrons are

\[ \begin{eqnarray} X(r,s,t) & \; = \; & r X_1 \; + \; s X_2 \; + \; t X_3 \; + \; (1 - r - s - t) X_4 \\ \\ Y(r,s,t) & \; = \; & r Y_1 \; + \; s Y_2 \; + \; t Y_3 \; + \; (1 - r - s - t) Y_4 \\ \\ Z(r,s,t) & \; = \; & r Z_1 \; + \; s Z_2 \; + \; t Z_3 \; + \; (1 - r - s - t) Z_4 \\ \end{eqnarray} \]

for the coordinates and

\[ \begin{eqnarray} u(r,s,t) & \; = \; & r u_1 \; + \; s u_2 \; + \; t u_3 \; + \; (1 - r - s - t) u_4 \\ \\ v(r,s,t) & \; = \; & r v_1 \; + \; s v_2 \; + \; t v_3 \; + \; (1 - r - s - t) v_4 \\ \\ w(r,s,t) & \; = \; & r w_1 \; + \; s w_2 \; + \; t w_3 \; + \; (1 - r - s - t) w_4 \\ \end{eqnarray} \]

for the displacements.

The derivatives are also easy in this case. The result is

\[ \begin{eqnarray} \left[ \matrix{ {\partial u \over \partial X} & {\partial u \over \partial Y} & {\partial u \over \partial Z} \\ \\ {\partial v \over \partial X} & {\partial v \over \partial Y} & {\partial v \over \partial Z} \\ \\ {\partial w \over \partial X} & {\partial w \over \partial Y} & {\partial w \over \partial Z} \\ } \right] = \left[ \matrix{ u_1 - u_4 & u_2 - u_4 & u_3 - u_4 \\ \\ v_1 - v_4 & v_2 - v_4 & v_3 - v_4 \\ \\ w_1 - w_4 & w_2 - w_4 & w_3 - w_4 \\ } \right] \left[ \matrix{ X_1 - X_4 & X_2 - X_4 & X_3 - X_4 \\ \\ Y_1 - Y_4 & Y_2 - Y_4 & Y_3 - Y_4 \\ \\ Z_1 - Z_4 & Z_2 - Z_4 & Z_3 - Z_4 \\ } \right]^{-1} \end{eqnarray} \]

Once again, the result is independent of \(r\) and \(s\), and now \(t\), which means that it is constant throughout the element. Not good!

Also, please consider visiting an advertiser above. Doing so helps to cover website hosting fees.

Bob McGinty

bmcginty@gmail.com

Click here to see a sample page in each of the two formats.

Click here to see a sample page in each of the two formats.

Copyright © 2012 by Bob McGinty