Introduction
A three-dimensional (3D) solid element can be considered to be the most general of all solid finite elements because all the field variables are dependent of x, y and z. An example of a 3D solid structure under loading is shown in Figure 9.1. As can be seen, the force vectors here can be in any arbitrary direction in space. A 3D solid can also have any arbitrary shape, material properties and boundary conditions in space. As such, there are altogether six possible stress components, three normal and three shear, that need to be taken into consideration. Typically, a 3D solid element can be a tetrahedron or hexahedron in shape with either flat or curved surfaces. Each node of the element will have three translational degrees of freedom. The element can thus deform in all three directions in space.
Since the 3D element is said to be the most general solid element, the truss, beam, plate, 2D solid and shell elements can all be considered to be special cases of the 3D element. So, why is there a need to develop all the other elements? Why not just use the 3D element to model everything? Theoretically, yes, the 3D element can actually be used to model all kinds of structurural components, including trusses, beams, plates, shells and so on. However, it can be very tedious in geometry creation and meshing. Furthermore, it is also most demanding on computer resources. Hence, the general rule of thumb is, that when a structure can be assumed within acceptable tolerances to be simplified into a 1D (trusses, beams and frames) or 2D (2D solids and plates) structure, always do so. The creation of a 1D or 2D FEM model is much easier and efficient. Use 3D solid elements only when we have no other choices. The formulation of 3D solids elements is straightforward, because it is basically an extension of 2D solids elements. All the techniques used in 2D solids can be utilized, except that all the variables are now functions of x, y and z. The basic concepts, procedures and formulations for 3D solid elements can also be found in many existing topics (see, e.g., Washizu, 1981; Rao, 1999; Zienkiewicz and Taylor, 2000; etc.).
Figure 9.1. Example of a 3D solid under loadings.
Tetrahedron Element
Strain Matrix
Consider the same 3D solid structure as Figure 9.1, whose domain is divided in a proper manner into a number of tetrahedron elements (Figure 9.2) with four nodes and four surfaces, as shown in Figure 9.3. A tetrahedron element has four nodes, each having three DOFs (u, v and w), making the total DOFs in a tetrahedron element twelve, as shown in Figure 9.3.
Figure 9.2. Solid block divided into four-node tetrahedron elements.
Figure 9.3. Atetrahedron element.
The nodes are numbered 1, 2, 3 and 4 by the right-hand rule. The local Cartesian coordinate system for a tetrahedron element can usually be the same as the global coordinate system, as there are no advantages in having a separate local Cartesian coordinate system. In an element, the displacement vector U is a function of the coordinate x, y and z, and is interpolated by shape functions in the following form, which should by now be shown to be part and parcel of the finite element method:
where the nodal displacement vector, de, is given as
and the matrix of shape functions has the form
To develop the shape functions, we make use of what is known as the volume coordinates, which is a natural extension from the area coordinates for 2D solids. The use of the volume coordinates makes it more convenient for shape function construction and element matrix integration. The volume coordinates for node 1 is defined as
where Vp234 and Vi234 denote, respectively, the volumes of the tetrahedrons P234 and 1234, as shown in Figure 9.4. The volume coordinate for node 2-4 can also be defined in the same manner:
Figure 9.4. Volume coordinates for tetrahedron elements.
The volume coordinate can also be viewed as the ratio between the distance of the point P and point 1 to the plane 234:
It can easily be confirmed that
since
It can also easily be confirmed that
Using Eq. (9.9), the relationship between the volume coordinates and Cartesian coordinates can be easily derived:
Equations (9.7) and (9.10) can then be expressed as a single matrix equation as follows:
The inversion of Eq. (9.11) will give
where
in which the subscript i varies from 1 to 4, and j, k and l are determined by a cyclic permutation in the order of i, j, k, l. For example, if i = 1, then j = 2, k = 3, l = 4. When i = 2, then j = 3, k = 4, l = 1. The volume of the tetrahedron element V can be obtained by
The properties of Li, as depicted in Eqs. (9.6) to (9.9), show that Li can be used as the shape function of a four-nodal tetrahedron element:
It can be seen from above that the shape function is a linear function of x, y and z, hence, the four-nodal tetrahedron element is a linear element. Note that from Eq. (9.14), the moment matrix of the linear basis functions will never be singular, unless the volume of the element is zero (or the four nodes of the element are in a plane). Based on Lemmas 2 and 3, we can be sure that the shape functions given by Eq. (9.15) satisfy the sufficient requirement of FEM shape functions.
It was mentioned that there are six stresses in a 3D element in total. The stress components are
To get the corresponding strains,
we can substitute Eq. (9.1) into Eq. (2.5):
where the strain matrix B is given by
Using Eq. (9.3), the strain matrix, B, can be obtained as
It can be seen that the strain matrix for a linear tetrahedron element is a constant matrix. This implies that the strain within a linear tetrahedron element is constant, and thus so is the stress. Therefore, the linear tetrahedron elements are also often referred to as a constant strain element or constant stress element, similar to the case of 2D linear triangular elements.
Element Matrices
Once the strain matrix has been obtained, the stiffness matrix ke for 3D solid elements can be obtained by substituting Eq. (9.18) into Eq. (3.71). Since the strain is constant, the element strain matrix is obtained as
Note that the material constant matrix c is given generally by Eq. (2.9). The mass matrix can similarly be obtained using Eq. (3.75):
where
Using the following formula [Eisenberg and Malvern, 1973],
we can conveniently evaluate the integral in Eq. (9.20) to give
An alternative way to calculate the mass matrix for 3D solid elements is to use a special natural coordinate system, which is defined as shown in Figures 9.5-9.7. In Figure 9.5, the plane of ξ = constant is defined in such a way that the edge P-Q stays parallel to the edge 2-3 of the element, and point 4 coincides with point 4 of the element. When P moves to point 1, ξ = 0, and when P moves to point 2, ξ = 1. In Figure 9.6, the plane of η = constant is defined in such a way that the edge 1-4 on the triangle coincides with the edge 1-4 of the element, and point P stays on the edge 2-3 of the element. When P moves to point 2, η = 0, and when P moves to point 3, η = 1. The plane of ζ = constant is defined in Figure 9.7, in such a way that the plane P-Q-R stays parallel to the plane 1 -2-3 of the element, and when P moves to point 4, ζ = 0, and when P moves to point 2, ζ = 1. In addition, the plane 1-2-3 on the element sits on the x-y plane.
Figure 9.5. Natural coordinate, where ξ = constant.
Figure 9.6. Natural coordinate, where η = constant.
Figure 9.7. Natural coordinate, where ζ = constant.
Therefore, the relationship between xyz and
‘ can be obtained in the following steps:
In Figure 9.8, the coordinates at point P are first interpolated using the x, y and z coordinates at points 2 and 3:
The coordinates at point B are then interpolated using the x, y and z coordinates at points 1 and P:
Figure 9.8. Cartesian coordinates xyz of point O in term of ξηζ.
The coordinates at point O are finally interpolated using the x, y and z coordinates at points 4 and B:
With this special natural coordinate system, the shape functions in the matrix of Eq. (9.3) can be written by inspection as
The Jacobian matrix between xyz and ξηζ is required, and is given as
Using Eqs. (9.26) and (9.27), the determinate of the Jacobian can be found to be
The mass matrix can now be obtained as
which gives
where Nij- is given by Eq. (9.21), but in which the shape functions should be defined by Eq. (9.27). Evaluating the integrals in Eq. (9.31) would give the same mass matrix as in Eq. (9.23).
The nodal force vector for 3D solid elements can be obtained using Eqs. (3.78), (3.79) and (3.81). Suppose the element is loaded by a distributed force f on the edge 2-3 of the element as shown in Figure 9.3; the nodal force vector becomes
If the load is uniformly distributed,
are constants,and the above equation becomes
where l3-4 is the length of the edge 3-4. Equation (9.33) implies that the distributed forces are equally divided and applied at the two nodes. This conclusion also applies to evenly distributed surface forces applied on any face of the element, and to evenly distributed body force applied on the entire body of the element. Finally, the stiffness matrix, ke, the mass matrix, me, and the nodal force vector, fe, can be used directly to assemble the global FE equation, Eq. (3.96), without going through a coordinate transformation.
