Fundamentals for Finite Element Method Part 3

Formation of FE Equations in Local Coordinate System

Once the shape functions are constructed, the FE equation for an element can be formulated using the following process. By substituting the interpolation of the nodes, Eq. (3.6), and the strain-displacement equation, say Eq. (2.5), into the strain energy term (Eq. (3.4)), we have

where the subscript e stands for the element. Note that the volume integration over the global domain has been changed to that over the elements. This can be done because we assume that the assumed displacement field satisfies the compatibility condition (see Section 3.3) on all the edges between the elements. In Eq.(3.69), B is often called the strain matrix, defined by

where L is the differential operator that is defined for different problems.For 3D solids, it is given by Eq. (2.7). By denoting

which is called the element stiffness matrix, Eq. (3.69) can be rewritten as

Note that the stiffness matrix ke is symmetrical, because

which shows that the transpose of matrix ke is itself. In deriving the above equation, c = cT has been employed. Making use of the symmetry of the stiffness matrix, only half of the terms in the matrix need to be evaluated and stored.

By substituting Eq. (3.6) into Eq. (3.3), the kinetic energy can be expressed as

By denoting

which is called the mass matrix of the element, Eq. (3.74) can be rewritten as

It is obvious that the element mass matrix is also symmetrical.

Finally, to obtain the work done by the external forces, Eq. (3.6) is substituted into Eq. (3.5):

where the surface integration is performed only for elements on the force boundary of the problem domain. By denoting

and

Eq. (3.77) can then be rewritten as

Fb and Fi are the nodal forces acting on the nodes of the elements, which are equivalent to the body forces and surface forces applied on the element in terms of the work done on a virtual displacement. These two nodal force vectors can then be added up to form the total nodal force vector fe :

Substituting Eqs. (3.72), (3.76) and (3.80) into Lagrangian functional L (Eq. (3.2)), we have

Applying Hamilton’s principle (Eq. (3.1)), we have

Note that the variation and integration operators are interchangeable, hence we obtain

To explicitly illustrate the process of deriving Eq. (3.84) from Eq. (3.83), we use a two-degree of freedom system as an example. Here, we show the procedure for deriving the second term in Eq. (3.84):

In Eq. (3.84), the variation and differentiation with time are also interchangeable, i.e.

Hence, by substituting Eq. (3.85) into Eq. (3.84), and integrating the first term by parts, we obtain

Note that in deriving Eq. (3.86) as above, the conditionhave    been used, which leads to the vanishing of the first term on the right-hand side. This is because the initial condition at t1 and final condition at t2 have to be satisfied for any de (admissible conditions (c) required by Hamilton’s principle), and no variation at t1 and t2 is allowed.

Substituting Eq. (3.86) into Eq. (3.84) leads to

To have the integration in Eq. (3.87) as zero for an arbitrary integrand, the integrand itself has to vanish, i.e.

Due to the arbitrary nature of the variation of the displacements, the only insurance for Eq. (3.88) to be satisfied is

Equation (3.89) is the FEM equation for an element, while ke and me are the stiffness and mass matrices for the element, and fe is the element force vector of the total external forces acting on the nodes of the element. All these element matrices and vectors can be obtained simply by integration for the given shape functions of displacements.

Coordinate transformation

The element equation given by Eq. (3.89) is formulated based on the local coordinate system defined on an element. In general, the structure is divided into many elements with different orientations (see Figure 3.4). To assemble all the element equations to form the global system equations, a coordinate transformation has to be performed for each element in order to formulate its element equation in reference to the global coordinate system which is defined on the whole structure.

The coordinate transformation gives the relationship between the displacement vector de based on the local coordinate system and the displacement vector De for the same element, but based on the global coordinate system:

Figure 3.4. Local and global coordinate system.

T is the transformation matrix, which has different forms depending upon the type of element, and will be discussed in detail in later topics. The transformation matrix can also be applied to the force vectors between the local and global coordinate systems:

in which Fe stands for the force vector at node i on the global coordinate system. Substitution of Eq. (3.90) into Eq. (3.89) leads to the element equation based on the global coordinate system:

where

Assembly of Global FE Equation

The FE equations for all the individual elements can be assembled together to form the global FE system equation:

where K and M are the globe stiffness and mass matrix, D is a vector of all the displacements at all the nodes in the entire problem domain, and F is a vector of all the equivalent nodal force vectors. The process of assembly is one of simply adding up the contributions from all the elements connected at a node.It may be noted here that the assembly of the global matrices can be skipped by combining assembling with the equation solving. This means that the assembling of a term in the global matrix is done only when the equation solver is operating on this term.

Imposition of Displacement Constraints

The global stiffness matrix K in Eq. (3.96) does not usually have a full rank, because displacement constraints (supports) are not yet imposed, and it is non-negative definite or positive semi-definite. Physically, an unconstrained solid or structure is capable of performing rigid movement. Therefore, if the solid or structure is free of support, Eq. (3.96) gives the behavior that includes the rigid body dynamics, if it is subjected to dynamic forces. If the external forces applied are static, the displacements cannot be uniquely determined from Eq. (3.96) for any given force vector. It is meaningless to try to determine the static displacements of an unconstrained solid or structure that can move freely.

For constrained solids and structures, the constraints can be imposed by simply removing the rows and columns corresponding to the constrained nodal displacements. We shall demonstrate this method in an example problem in later topics. After the treatment of constraints (and if the constraints are sufficient), the stiffness matrix K in Eq. (3.96) will be of full rank, and will be Positive Definite (PD). Since we have already proven that K is symmetric, K is of a Symmetric Positive Definite (SPD) property.

Solving the Global FE Equation

By solving the global FE equation, displacements at the nodes can be obtained. The strain and stress in any element can then be retrieved using Eq. (3.6) in Eqs. (2.5) and (2.8).

Static Analysis

Static analysis involves the solving of Eq. (3.96) without the term with the global mass matrix, M. Hence, as discussed, the static system of equations takes the form

There are numerous methods and algorithms to solve the above matrix equation. The methods often used are Gauss elimination and LU decompositions for small systems, and iterative methods for large systems. These methods are all routinely available in any math library of any computer system.

Analysis of Free Vibration (EigenvALUE Analysis)

For a structural system with a total DOF of N, the stiffness matrix K and mass matrix M in Eq. (3.96) have a dimension of N x N .By solving the above equation we can obtain the displacement field, and the stress and strain can then be calculated. The question now is how to solve this equation, as N is usually very large for practical engineering structures. One way to solve this equation is using the so-called direct integration method, discussed in the next section. An alternative way of solving Eq. (3.96) is using the so-called modal analysis technique (or mode superposition technique). In this technique, we first have to solve the homogenous equation of Eq. (3.96). The homogeneous equation is when we consider the case of F = 0, therefore it is also called free vibration analysis, as the system is free of external forces. For a solid or structure that undergoes a free vibration, the discretized system equation Eq. (3.96) becomes

This solution for the free vibration problem can be assumed as

where φ is the amplitude of the nodal displacement, ω is the frequency of the free vibration, and t is the time. By substituting Eq. (3.99) into Eq. (3.98), we obtain

or

where

Equation (3.100) (or (3.101)) is called the eigenvalue equation. To have a non-zero solution for φ, the determinate of the matrix must vanish:

The expansion of the above equation will lead to a polynomial of λ of order N. This polynomial equation will have N roots,called eigenvalues, which relate to the naturalfrequency of the system by Eq. (3.100). The natural frequency is a very important characteristic of the structure carrying dynamic loads. It has been found that if a structure is excited by a load with a frequency of one of the structure’s natural frequencies, the structure can undergo extremely violent vibration, which often leads to catastrophic failure of the structural system. Such a phenomenon is called resonance. Therefore, an eigenvalue analysis has to be performed in designing a structural system that is to be subjected to dynamic loadings.

By substituting an eigenvalue Xi back into the eigenvalue equation, Eq. (3.101), we have

which is a set of algebraic equations. Solving the above equation for φ, a vector denoted by φ{ can then be obtained. This vector corresponding to the ith eigenvalue Xi is called the ith eigenvector that satisfies the following equation:

An eigenvector φi corresponds to a vibration mode that gives the shape of the vibrating structure of the ith mode. Therefore, analysis of the eigenvalue equation also gives very important information on possible vibration modes experienced by the structure when it undergoes a vibration. Vibration modes of a structure are therefore another important characteristic of the structure. Mathematically, the eigenvectors can be used to construct the displacement fields. It has been found that using a few of the lowest modes can obtain very accurate results for many engineering problems. Modal analysis techniques have been developed to take advantage of these properties of natural modes.

In Eq. (3.101), since the mass matrix M is symmetric positive definite and the stiffness matrix K is symmetric and either positive or positive semi-definite, the eigenvalues are all real and either positive or zero. It is possible that some of the eigenvalues may coincide. The corresponding eigenvalue equation is said to have multiple eigenvalues. If there are m coincident eigenvalues, the eigenvalue is said to be of a multiplicity of m. For an eigenvalue of multiplicity m, there are m vectors satisfying Eq. (3.105).

Methods for the effective computation of the eigenvalues and eigenvectors for an eigenvalue equations system like Eq. (3.100) or (3.101) are outside the scope of this topic. Intensive research has been conducted to-date, and many sophisticated algorithms are already well established and readily available in the open literature, and routinely in computational libraries. The commonly used methods are (see, e.g. Petyt, 1990):

•    Jacobi’s method;

•    Given’s method and householder’s method;

•    the bisection method (using Sturm sequences);

•    inverse iteration;

•    QR method;

•    subspace iteration;

•    Lanczos’ method.