On the shape functions and their requirements for a finite element analysis 1. Introduction:
This write up provides an insight into the meaning of shape functions and their requirements for a physically admissible solution.
Section 2 discusses in detail the definition of shape functions and section 3 gives the derivation of shape functions for a C0 and C1 continuous linear truss and flexural element.
2. What are shape functions?
Let us consider a two dimensional domain as shown in the Figure 1.We have a triangular finite element indicated in the domain.
Figure 1: Triangular finite element represented in the domain
Let us say Φ (x,y) is the field variable associated with the domain. The triangular element indicated in the figure has only external nodes. It may be noted that a node is specific point in the domain where the value of the field variable is explicitly computed. If the value of the field variable is explicitly computed at the nodes, how are the values at other points within a finite element computed? The values of the field variable computed at the nodes are used to approximate the values at non-nodal points (that is, in the element interior) by interpolation of the non-nodal values. For the three noded triangular element, the nodes are all at the exterior and at any other point within the element the field variable is described by the approximate relation;
Where;
Φ1, Φ2, Φ3, are the values of the field variables at the nodes and N1, N2, N3 are the interpolation functions also known as shape functions or blending functions. In the finite element analysis approach, the nodal values of the field variables are treated as unknown constants that are to be determined. The interpolation functions are most often polynomial forms of the independent variables, derived to satisy certain conditions at the nodes. The major point to be made here is that the interpolation functions are pre-determined, known functions of the independent variables; and these functions describe the variation of the field variable within the element.
3. Degree of continuity of shape functions
Filed quantity Φ is interpolated in a piecewise fashion over an element. That is, each interpolation piece is defined only within its element. So, while Φ can be guaranteed to vary smoothly within each element, the transition between elements may not be smooth. The symbol Cm is used to define the continuity of a piecewise field.
A field is Cm continuous if its derivatives upto and including the degree m are inter-element continuous.
Thus, in one dimension Φ = Φx is C0 continuous if Φ is continuous and its derivative is not, and Φ = Φx is C1 continuous if Φ = Φx and its first derivative is continuous while its second derivative is not.
The Cm terminology can be applied to element types. The bar element is C0 continuous and the beam element is C1 continuous.3.1 C0 interpolation: Shape functions of a linear truss element
Let us consider a field Φ which is C0 continuous. Thus, the field variable only is only inter-element continuous whereas its derivatives are not.
Expressing the variable Φ in terms of a linear interpolating polynomial, we have;
3.2 C1 interpolation: Shape functions of a flexural element
Recalling that the basic premise of finite element formulation is to express the continuously varying field variable in terms of finite number of values evaluated at the nodes, we note that, for flexural element the field variable of interest is the transverse displacement v(x) of the neutral surface away from its straight un-deflected position.
As seen in the figure below, the transverse displacement of the beam is such that the variation of deflection along the length is not adequately described by the displacements of the end points only;
As seen from the figure 2 a and b, the end points can be identical while the deflected shape of the two cases is quite different. Therefore, the flexural element formulation should take into account the slope of the beam as well as the end point displacement. In addition to avoiding the potential ambiguity of displacements, inclusion of beam element nodal rotations ensures compatibility of rotations at nodal connections between the elements, this precluding the physical unacceptable discontinuity depicted in figure 2c).
In light of these observations regarding rotations, the nodal variables to be associated with flexural element are depicted as shown in Figure 3 below. The nodal variables are denoted in positive direction and it is to be noted that the slopes are to be specified in radians.
Figure 4: Flexural element nodal variables shown in positive sense.
The displacement function is to be discretised such that;
Equation 16
subject to the boundary conditions;
Considering the four boundary conditions, we assume the displacement function as;
The choice of a cubic function to describe the displacement is not arbitrary because;
i. With the specification of four boundary conditions, we can determine not more than four constants in the assumed displacement function
ii. Considering the differential equation;
Equation 19
We see that, since the loads are applied only at the nodes, the second derivative of the displacement function is linear, hence the bending moment varies linearly along the length of the element.
Application of the boundary condition to Eq17 yields the following;
Solving Equation 20 simultaneously gives;
N1,N2,N3,N4 are interpolation functions which describe the distribution of displacement in terms of nodal values in the nodal displacement vector {δ}
4. Requirements of shape functions
4.1 Solution accuracy:
In finite elements analysis, the solution accuracy is judged in terms of convergence as the element “mesh” is refined. There are two major methods of mesh refinement:
i. h-refinement
ii. p-refinement
In the h-refinement, mesh refinement refers to the process of increasing the number of elements used to model a given domain, consequently reducing the individual element size.
In p-refinement, the element size is unchanged but the order of polynomials used as interpolation functions is increased. The objective of mesh refinement in either method is to obtain sequential solutions that exhibit asymptotic convergence to the value representing the exact solution.
3.2 C1 interpolation: Shape functions of a flexural element
Recalling that the basic premise of finite element formulation is to express the continuously varying field variable in terms of finite number of values evaluated at the nodes, we note that, for flexural element the field variable of interest is the transverse displacement v(x) of the neutral surface away from its straight un-deflected position.
As seen in the figure below, the transverse displacement of the beam is such that the variation of deflection along the length is not adequately described by the displacements of the end points only;
As seen from the figure 2 a and b, the end points can be identical while the deflected shape of the two cases is quite different. Therefore, the flexural element formulation should take into account the slope of the beam as well as the end point displacement. In addition to avoiding the potential ambiguity of displacements, inclusion of beam element nodal rotations ensures compatibility of rotations at nodal connections between the elements, this precluding the physical unacceptable discontinuity depicted in figure 2c).
In light of these observations regarding rotations, the nodal variables to be associated with flexural element are depicted as shown in Figure 3 below. The nodal variables are denoted in positive direction and it is to be noted that the slopes are to be specified in radians.
Figure 4: Flexural element nodal variables shown in positive sense.
The displacement function is to be discretised such that;
Equation 16
subject to the boundary conditions;
Considering the four boundary conditions, we assume the displacement function as;
The choice of a cubic function to describe the displacement is not arbitrary because;
i. With the specification of four boundary conditions, we can determine not more than four constants in the assumed displacement function
ii. Considering the differential equation;
Equation 19
We see that, since the loads are applied only at the nodes, the second derivative of the displacement function is linear, hence the bending moment varies linearly along the length of the element.
Application of the boundary condition to Eq17 yields the following;
Solving Equation 20 simultaneously gives;
N1,N2,N3,N4 are interpolation functions which describe the distribution of displacement in terms of nodal values in the nodal displacement vector {δ}
4. Requirements of shape functions
4.1 Solution accuracy:
In finite elements analysis, the solution accuracy is judged in terms of convergence as the element “mesh” is refined. There are two major methods of mesh refinement:
i. h-refinement
ii. p-refinement
In the h-refinement, mesh refinement refers to the process of increasing the number of elements used to model a given domain, consequently reducing the individual element size.
An example illustrating regular h-refinement as well as solution convergence is shown in the Figure 5 below, which depicts a regular elastic plate of uniform thickness fixed on one edge and subjected to a concentrated load on one corner:
The problem is modelled using rectangular plane stress elements and three meshes are used in sequence Figure 6.1 b- 6.1 d. Solution convergence is depicted in Figure 6.1 e in terms of maximal normal stress in the x-direction. For this problem, the exact solution is taken in terms of maximal normal stress in x-direction. For this problem, the exact solution is taken to be the maxiamal normal bending stress computed using elemental beam theory. The true exact solution is the plane stress solution is the plane stress solution from the theory of elasticity. For a general field problem, the field variable is expressed on an element basis in the discretised form:
4.2 Compatibility requirement
Along element boundaries, the field variable and its partial derivatives up to one order less than the highest order derivative appearing in the integral formulation of the element equations must be continuous. Given the discretised representation of equation 6.1, it follows that the interpolation functions must meet this condition, since these functions determine the spatial variation of the field variable. Similarly, the beam element includes he second derivative of displacement and compatibility requires continuity of both the displacement and the slope (first derivative) at the element boundaries.
Physical meaning of compatibility requirements
In addition to satisfying the criteria of convergence, the compatibility requirements give a physical meaning to the problem as well. In structural engineering problems, the requirement of displacement continuity along the element boundaries ensures that no gaps or voids develop in the structure as a result of modelling procedure.
4.3 Completeness requirement
In the limit as element size shrinks to zero in the mesh refinement, the field variable and its partial derivatives up to and including, the highest order derivative appearing in the integral formulation must be capable of assuming constant values. Again, because of discretisation, the completeness requirement is directly applicable to interpolation functions. The completeness requirement ensures that a displacement field within a structural element can take on a constant value (constant values within that element but not overall), representing a rigid body motion for example. Similarly, constant slope for a beam element represents a rigid body rotation for example while a state of constant temperature in thermal element corresponds to no heat flux through the element.
The problem is modelled using rectangular plane stress elements and three meshes are used in sequence Figure 6.1 b- 6.1 d. Solution convergence is depicted in Figure 6.1 e in terms of maximal normal stress in the x-direction. For this problem, the exact solution is taken in terms of maximal normal stress in x-direction. For this problem, the exact solution is taken to be the maxiamal normal bending stress computed using elemental beam theory. The true exact solution is the plane stress solution is the plane stress solution from the theory of elasticity.
For a general field problem, the field variable is expressed on an element basis in the discretised form:
Where; M is the number of element degrees of freedom, the interpolation functions must satisfy two primary conditions to ensure convergence during mesh refinement: the compatibility and completeness requirements discussed as below.
4.2 Compatibility requirement
Along element boundaries, the field variable and its partial derivatives up to one order less than the highest order derivative appearing in the integral formulation of the element equations must be continuous. Given the discretised representation of equation 6.1, it follows that the interpolation functions must meet this condition, since these functions determine the spatial variation of the field variable.
Considering the formulation of truss elements, the displacements must be continuous across element boundaries but none of the displacement derivatives is required to be continuous across such boundaries. Indeed, the truss element is a constant strain element so the first order derivative is, in general, discontinuous at the boundaries.
Similarly, the beam element includes he second derivative of displacement and compatibility requires continuity of both the displacement and the slope (first derivative) at the element boundaries.
Physical meaning of compatibility requirements
In addition to satisfying the criteria of convergence, the compatibility requirements give a physical meaning to the problem as well. In structural engineering problems, the requirement of displacement continuity along the element boundaries ensures that no gaps or voids develop in the structure as a result of modelling procedure.
Similarly, the requirement of slope continuity for the beam element ensures that no „kinks‟ are developed in the deformed structure.
In heat transfer problems, the compatibility requirement prevents the physically unacceptable possibility of jump discontinuity in temperature distribution.
4.3 Completeness requirement
In the limit as element size shrinks to zero in the mesh refinement, the field variable and its partial derivatives up to and including, the highest order derivative appearing in the integral formulation must be capable of assuming constant values. Again, because of discretisation, the completeness requirement is directly applicable to interpolation functions.
