Introduction
In each of these elements, the profile of the displacements is assumed in simple forms to obtain element equations. The equations obtained for each element are then assembled together with adjoining elements to form the global finite element equation for the whole problem domain. Equations thus created for the global problem domain can be solved easily for the entire displacement field.
The above-mentioned FEM process does not seem to be a difficult task. However, in close examination of the above-mentioned process, one would naturally ask a series of questions. How can one simply assume the profile of the solution of the displacement in any simple form? How can one ensure that the governing partial differential equations will be satisfied by the assumed displacement? What should one do when using the assumed profile of the displacement to determine the final displacement field? Yes, one can just simply assume the profile of the solution of the displacements, but a principle has to be followed in order to obtain discretized system equations that can be solved routinely for the final displacement field. The use of such a principle guarantees the best (not exact) satisfaction of the governing system equation under certain conditions. The following details one of the most important principles, which will be employed to establish the FEM equations for mechanics problems of solids and structures.
It may be common for a beginner to be daunted by the equations involved when formulating the finite element equations. Perhaps for a beginner, full understanding of the details of the equations in this topic may prove to be a challenging task. It is thus advised that the novice reader just understands the basic concepts involved without digging too much into the equations. It is then recommended to review this topic again after going through subsequent topics with examples to fully understand the equations.
Strong and Weak Forms
The strong form, in contrast to a weak form, requires strong continuity on the dependent field variables (the displacements u, v and w in this case). Whatever functions that define these field variables have to be differentiable up to the order of the partial differential equations that exist in the strong form of the system equations. Obtaining the exact solution for a strong form of the system equation is usually very difficult for practical engineering problems. The finite difference method can be used to solve system equations of the strong form to obtain an approximated solution. However, the method usually works well for problems with simple and regular geometry and boundary conditions.
A weak form of the system equations is usually created using one of the following widely used methods:
• Energy principles (see, e.g. Washizu, 1974; Reddy, 1984)
• Weighted residual methods (see, e.g. Crandall, 1956; Finlayson and Scriven, 1966;
Finlayson, 1972; Ziekiewicz and Taylor, 2000)
The energy principle can be categorized as a special form of the variational principle which is particularly suited for problems of the mechanics of solids and structures. The weighted residual method is a more general mathematical tool applicable, in principle, for solving all kinds of partial differential equations. Both methods are very easy to understand and apply. This topic will demonstrate both methods used for creating FEM equations. An energy principle will be used for mechanics problems of solids and structures, and the weighted residual method will be used for formulating the heat transfer problems. It is also equally applicable to use the energy principle for heat transfer problems, and the weighted residual method for solid mechanics problems, and the procedure is very much the same.
The weak form is often an integral form and requires a weaker continuity on the field variables. Due to the weaker requirement on the field variables, and the integral operation, a formulation based on a weak form usually produces a set of discretized system equations that give much more accurate results, especially for problems of complex geometry. Hence, the weak form is preferred by many for obtaining an approximate solution. The finite element method is a typical example of successfully using weak form formulations. Using the weak form usually leads to a set of well-behaved algebraic system equations, if the problem domain is discretized properly into elements. As the problem domain can be discretized into different types of elements, the FEM can be applied for many practical engineering problems with most kinds of complex geometry and boundary conditions.
In the following section, Hamilton’s principle, which is one of the most powerful energy principles, is introduced for FEM formulation of problems of mechanics of solids and structures. Hamilton’s principle is chosen because it is simple and can be used for dynamic problems. The approach adopted in this topic is to directly work out the dynamic system equations, after which the static system equations can be easily obtained by simply dropping out the dynamic terms. This can be done because of the simple fact that the dynamic system equations are the general system equations, and the static case can be considered to be just a special case of the dynamic equations.
Hamilton’s principle
Hamilton’s principle is a simple yet powerful tool that can be used to derive discretized dynamic system equations. It states simply that “Of all the admissible time histories of displacement the most accurate solution makes the Lagrangian functional a minimum.”
An admissible displacement must satisfy the following conditions:
(a) the compatibility equations,
(b) the essential or the kinematic boundary conditions, and
(c) the conditions at initial (ίχ) and final time (t2).
Condition (b) ensures that the displacement constraints are satisfied; and condition (c) requires the displacement history to satisfy the constraints at the initial and final times.
Mathematically, Hamilton’s principle states:
The Langrangian functional, L, is obtained using a set of admissible time histories of displacements, and it consists of
where T is the kinetic energy, Π is the potential energy (for our purposes, it is the elastic strain energy), and Wf is the work done by the external forces. The kinetic energy of the entire problem domain is defined in the integral form
where V represents the whole volume of the solid, and U is the set of admissible time histories of displacements.
The strain energy in the entire domain of elastic solids and structures can be expressed as
where ε are the strains obtained using the set of admissible time histories of displacements. The work done by the external forces over the set of admissible time histories of displacements can be obtained by
where Sf represents the surface of the solid on which the surface forces are prescribed (see Figure 2.2).
Hamilton’s principle allows one to simply assume any set of displacements, as long as it satisfies the three admissible conditions. The assumed set of displacements will not usually satisfy the strong form of governing system equations unless we are extremely lucky, or the problem is extremely simple and we know the exact solution. Application of Hamilton’s principle will conveniently guarantee a combination of this assumed set of displacements to produce the most accurate solution for the system that is governed by the strong form of the system equations.
The power of Hamilton’s principle (or any other variational principle) is that it provides the freedom of choice, opportunity and possibility. For practical engineering problems, one usually does not have to pursue the exact solution, which in most cases are usually unobtainable, because we now have a choice to quite conveniently obtain a good approximation using Hamilton’s principle, by assuming the likely form, pattern or shape of the solutions. Hamilton’s principle thus provides the foundation for the finite element methods. Furthermore, the simplicity of Hamilton’s principle (or any other energy principle) manifests itself in the use of scalar energy quantities. Engineers and scientists like working with scalar quantities when it comes to numerical calculations, as they do not need to worry about the direction. All the mathematical tools required to derive the final discrete system equations are then basic integration, differentiation and variation, all of which are standard linear operations. Another plus point of Hamilton’s principle is that the final discrete system equations produced are usually a set of linear algebraic equations that can be solved using conventional methods and standard computational routines. The following demonstrates how the finite element equations can be established using Hamilton’s principle and its simple operations.
FEM procedure
The standard FEM procedure can be briefly summarized as follows.
Domain Discretization
The solid body is divided into Ne elements. The procedure is often called meshing, which is usually performed using so-called pre-processors. This is especially true for complex geometries. Figure 3.1 shows an example of a mesh for a two-dimensional solid.
The pre-processor generates unique numbers for all the elements and nodes for the solid or structure in a proper manner. An element is formed by connecting its nodes in a pre-defined consistent fashion to create the connectivity of the element. All the elements together form the entire domain of the problem without any gap or overlapping. It is possible for the domain to consist of different types of elements with different numbers of nodes, as long as they are compatible (no gaps and overlapping; the admissible condition (a) required by Hamilton’s principle) on the boundaries between different elements. The density of the mesh depends upon the accuracy requirement of the analysis and the computational resources available. Generally, a finer mesh will yield results that are more accurate, but will increase the computational cost. As such, the mesh is usually not uniform, with a finer mesh being used in the areas where the displacement gradient is larger or where the accuracy is critical to the analysis. The purpose of the domain discretization is to make it easier in assuming the pattern of the displacement field.
Figure 3.1. Example of a mesh with elements and node properly numbered.
Displacement Interpolation
The FEM formulation has to be based on a coordinate system. In formulating FEM equations for elements, it is often convenient to use a local coordinate system that is defined for an element in reference to the global coordination system that is usually defined for the entire structure, as shown in Figure 3.4. Based on the local coordinate system defined on an element, the displacement within the element is now assumed simply by polynomial interpolation using the displacements at its nodes (or nodal displacements) as
where the superscript h stands for approximation, nd is the number of nodes forming the element, and d, is the nodal displacement at the i th node, which is the unknown the analyst wants to compute, and can be expressed in a general form of
where nf is the number of Degrees Of Freedom (DOF) at a node. For 3D solids, nf = 3, and
Note that the displacement components can also consist of rotations for structures of beams and plates. The vector de in Eq. (3.6) is the displacement vector for the entire element, and has the form of
Therefore, the total DOF for the entire element is nd x nf.
In Eq. (3.6), N is a matrix of shape functions for the nodes in the element, which are predefined to assume the shapes of the displacement variations with respect to the coordinates. It has the general form of
where Ni is a sub-matrix of shape functions for displacement components, which arranged as
where Nik is the shape function for the kth displacement component (DOF) at the ith node. For 3D solids, nf = 3, and often
Note that it is not necessary to use the same shape function for all the displacement components at a node. For example, we often use different shape functions for translational and rotational displacements.
Note that this approach of assuming the displacements is often called the displacement method. There are FEM approaches that assume the stresses instead, but they will not be covered in this topic.
Standard Procedure for Constructing Shape Functions
Consider an element with nd nodes at
for one dimensional problems,
for two-dimensional problems, and for three-dimensional problems. We should have nd shape functions for each displacement component for an element. In the following, we consider only one displacement component in the explanation of the standard procedure for constructing the shape functions. The standard procedure is applicable for any other displacement components. First, the displacement component is approximated in the form of a linear combination of nd linearly-independent basis functions Pi (x), i.e.
where uh is the approximation of the displacement component, Pi (x) is the basis function of monomials in the space coordinates x, and ai is the coefficient for the monomial Pi (x). Vector a is defined as
The pi(x) in Eq. (3.12) is built with nd terms of one-dimensional monomials; based on the Pascal’s triangle shown in Figure 3.2 for two-dimensional problems; or the well-known Pascal’s pyramid shown in Figure 3.3 for three-dimensional problems. A basis of complete order of P in the one-dimensional domain has the form
A basis of complete order of p in the two-dimensional domain is provided by
and that in three-dimensional domain can be written as
As a general rule, the nd terms of Pi (x) used in the basis should be selected from the constant term to higher orders symmetrically from the Pascal triangle shown in Figures 3.2 or 3.3. Some higher-order terms can be selectively included in the polynomial basis if there is a need in specific circumstances.
Figure 3.2. Pascal triangle of monomials (two-dimensional case).
Figure 3.3. Pascal pyramid of monomials (three-dimensional case).
The coefficients ai in Eq. (3.12) can be determined by enforcing the displacements calculated using Eq. (3.12) to be equal to the nodal displacements at the nd nodes of the element. At node i we can have
where di is the nodal value of
Equation (3.17) can be written in the following matrix form:
where de is the vector that includes the values of the displacement component at all the nd nodes in the element:
and P is given by
which is called the moment matrix. The expanded form of P is
For two-dimensional polynomial basis functions, we have
Using Eq. (3.18), and assuming that the inverse of the moment matrix P exists, we can then have
Substituting Eq. (3.23) into Eq. (3.12), we then obtain
or in matrix form
where N(x) is a matrix of shape functions Ni (x) defined by
where
is the ith column of matrix
and
In obtaining Eq. (3.23), we have assumed the existence of the inverse of P. There could be cases where P-1 does not exist, and the construction of shape functions will fail. The existence of P-1 depends upon (1) the basis function used, and (2) the nodal distribution of the element. The basis functions have to be chosen first from a linearly-independent set of bases, and then the inclusion of the basis terms should be based on the nodal distribution in the element. The discussion in this direction is more involved.In this topic, we shall only discuss elements whose corresponding moment matrices are invertible.
Note that the derivatives of the shape functions can be obtained very easily, as all the functions involved are polynomials. The lth derivative of the shape functions is simply given by
Note that there are many other methods for creating shape functions which do not necessarily follow the standard procedure described above. Some of these often-used shortcut methods will be discussed in later topics, when we develop different types of elements. These shortcut methods need to make use of the properties of shape functions detailed in the next section.
