Adaptive finite element methods for differential equations

finite element methods for navier-stokes equations theory and algorithms and a finite element method for elliptic equations on surfaces
Dr.NaveenBansal Profile Pic
Dr.NaveenBansal,India,Teacher
Published Date:25-10-2017
Your Website URL(Optional)
Comment
7 FEM FOR TWO-DIMENSIONAL SOLIDS 7.1 INTRODUCTION In this chapter, we develop, in an easy to understand manner, finite element equations for the stress analysis of two-dimensional (2D) solids subjected to external loads. The basic concepts, procedures and formulations can also be found in many existing textbooks (see, e.g. Zienkiewicz and Taylor, 2000). The element developed is called a 2D solid element that is used for structural problems where the loading–and hence the deformation–occur within a plane. Though no real life structure can be truly 2D, experienced analysts can often idealize many practical problems to 2D problems to obtain satisfactory results by carrying out analyses using 2D models, which can be very much more efficient and cost-effective compared to conducting full 3D analyses. In engineering applications, there are ample practical problems that can be modelled as 2D problems. As discussed in Chapter 2, there are plane stress and plane strain problems, whereby correspondingly, plane stress and plane strain elements need to be used to solve them. For example, if we have a plate structure with loading acting in the plane of the plate as in Figure 7.1, we need to use 2D plane stress elements. When we want to model the effects of water pressure on a dam, as shown in Figure 7.2, we have to use 2D plane strain elements. y f y f x x Figure 7.1. A typical 2D plane stress problem. 129 “chap07” — 2002/12/14 — page 129 — 1130 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS z f x y x Figure 7.2. A typical 2D plane strain problem. Note that in Figure 7.1, plane stress conditions are usually applied to structures that have a relatively small thickness as compared to its other dimensions. Due to the absence of any off-plane external force, the normal stresses are negligible, which leads to a plane stress sit- uation. In cases where plane strain conditions are applied, as in Figure 7.2, the thickness of the structure (in the z direction) is relatively large as compared to its other dimensions, and the loading (pressure) is uniform along the elongated direction. The deformation is, there- fore, approximated to be the same throughout its thickness. In this case, the off-plane strain (strain components in the z direction) is negligible, which leads to a plane strain situation. In either a plane stress or plane strain situation, the governing system equation can be dras- tically simplified, as shown in Chapter 2. The formulations for plane stress and plane strain problems are very much the same, except for the difference in the material constant matrix. A 2D solid element, be it plane strain or plane stress, can be triangular, rectangular or quadrilateral in shape with straight or curved edges. The most often used elements in engineering practice are linear. Quadratic elements are also used for situations that require high accuracy in stress, but they are less often used for practical problems. Higher order elements have also been developed, but they are less often used except for certain specific problems. The order of the 2D element is determined by the order of the shape functions used. A linear element uses linear shape functions, and therefore the edges of the element are straight. A quadratic element uses quadratic shape functions, and their edges can be curved. The same can be said for elements of third order or higher. In a 2D model, the elements can only deform in the plane where the model is defined, and in most situations, this is the x–y plane. At any point, the variable, that is the displacement, has two components in the x and y directions, and so do the external forces. For plane strain problems, the thickness of the true structure is usually not important, and is normally treated as a unit quantity uniformly throughout the 2D model. However, for plane stress “chap07” — 2002/12/14 — page 130 — 27.2 LINEAR TRIANGULAR ELEMENTS 131 problems, the thickness is an important parameter for computing the stiffness matrix and stresses. Throughout this chapter, it is assumed that the elements have a uniform thickness of h. If the structure to be modelled has a varying thickness, the structure needs to be divided into small elements, where in each element a uniform thickness can be used. On the other hand, formulation of 2D elements with varying thicknesses can also be done easily, as the procedure is similar to that of a uniform element. Very few commercially available software packages provide elements of varying thickness. The equation system for a 2D element will be more complex as compared with the 1D element because of the higher dimension. The procedure for developing these equations is, however, very similar to that for the 1D truss elements, which is detailed in Chapter 4. These steps can be summarized in the following three-step procedure: 1. Construction of shape functions matrix N that satisfies Eqs. (3.34) and (3.41). 2. Formulation of the strain matrix B. 3. Calculation of k , m , and f using N and B and Eqs. (3.71), (3.75) and (3.81). e e e We shall be focusing on the formulation of three types of simple but very important elements: linear triangular, bilinear rectangular, and isoparametric linear quadrilateral ele- ments. Once the formulation of these three types of element is understood, the development of other types of elements of higher orders is straightforward, because the same techniques can be utilized. Development of higher order elements will be discussed at the end of this chapter. 7.2 LINEAR TRIANGULAR ELEMENTS The linear triangular element was the first type of element developed for 2D solids. The formulation is also the simplest among all the 2D solid elements. It has been found that the linear triangular element is less accurate compared to linear quadrilateral elements. For this reason, it is often thought to be ideal to use quadrilateral elements, but the reality is that the triangular element is still a very useful element for its adaptivity to complex geometry. Triangular elements are normally used when we want to mesh a 2D model involving complex geometry with acute corners. In addition, the triangular configuration with the simplest topological feature makes it easier to develop meshing processors. Nowadays, analysts are hoping to use a fully automated mesh generator to perform the complex task of analysis that needs repeated or even adaptive re-meshing. Most automated mesh generators can only create triangular elements. There are automated mesh generators that can generate a quadrilateral mesh, but they still use triangular elements as patches for difficult situations, and end up with a mesh of mixed elements. Hence, whether we like it or not, we still have to use triangular elements for many practical engineering problems. Consider a 2D model in the x–y plane, shown schematically in Figure 7.3. The 2D domain is divided in a proper manner into a number of triangular elements. The ‘proper’ meshing of a domain will be outlined in Chapter 11, where a list of guidelines is provided. In a mesh of linear triangular elements, each triangular element has three nodes and three straight edges. “chap07” — 2002/12/14 — page 131 — 3132 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS Figure 7.3. Rectangular domain meshed with triangular elements. 3(x , y ) y, v 3 3 (u , v ) 3 3 fsy fsx A 2(x , y ) 2 2 (u , v ) 2 2 1(x , y ) 1 1 (u , v ) 1 1 x, u Figure 7.4. Linear triangular element. 7.2.1 Field Variable Interpolation Consider now a triangular element of thickness h. The nodes of the element are numbered 1, 2 and 3 counter-clockwise, as shown in Figure 7.4. For 2D solid elements, the field variable is the displacement, which has two components (u and v), and hence each node has two Degrees Of Freedom (DOFs). Since a linear triangular element has three nodes, the total number of DOFs of a linear triangular element is six. For the triangular element, the local coordinate of each element can be taken as the same as the global coordinate, since there are no advantages in specifying a different local coordinate system for each element. Now, let us examine how a triangular element can be formulated. The displacement U is generally a function of the coordinates x and y, and we express the displacement at any point in the element using the displacements at the nodes and shape functions. It is therefore assumed that (see Section 3.4.2) h U (x, y)= N(x, y)d (7.1) e where the superscript h indicates that the displacement is approximated, and d is a vector e of the nodal displacements arranged in the order of    u   1     displacements at node 1     v  1       u 2 d = displacements at node 2 (7.2) e   v 2            u 3     displacements at node 3   v 3 “chap07” — 2002/12/14 — page 132 — 47.2 LINEAR TRIANGULAR ELEMENTS 133 and the matrix of shape functions N is arranged as N 0 N 0 N 0 1 2 3 N= 0 N 0 N 0 N 1 2 3 (7.3)    Node 1 Node 2 Node 3 in which N (i = 1, 2, 3) are three shape functions corresponding to the three nodes of the i triangular element. Equation (7.1) can be explicitly expressed as h u (x, y)= N (x, y)u + N (x, y)u + N (x, y)u 1 1 2 2 3 3 (7.4) h v (x, y)= N (x, y)v + N (x, y)v + N (x, y)v 1 1 2 2 3 3 which implies that each of the displacement components at any point in the element is approximated by an interpolation from the nodal displacements using the shape functions. This is because the two displacement components are basically independent from each other. The question now is how can we construct these shape functions for our triangular element that satisfies the sufficient requirements: delta function property; partitions of unity; and linear field reproduction. 7.2.2 Shape Function Construction Development of the shape functions is normally the first, and most important, step in developing finite element equations for any type of element. In determining the shape functions N (i = 1, 2, 3) for the triangular element, we can of course follow exactly the i standard procedure described in Sections 3.4.3 and 4.2.1, by starting with an assumption of the displacements using polynomial basis functions with unknown constants. These unknown constants are then determined using the nodal displacements at the nodes of the element. This standard procedure works in principle for the development of any type of element, but may not be the most convenient method. We demonstrate here another slightly different approach for constructing shape functions. We start with an assumption of shape functions directly using polynomial basis functions with unknown constants. These unknown constants are then determined using the property of the shape functions. The only difference here is that we assume directly the shape function instead of the displacements. For a linear triangular element, we assume that the shape functions are linear functions of x and y. They should, therefore, have the form of N = a + b x+ c y (7.5) 1 1 1 1 N = a + b x+ c y (7.6) 2 2 2 2 N = a + b x+ c y (7.7) 3 3 3 3 where a ,b and c (i = 1, 2, 3) are constants to be determined. Equation (7.5) can be i i i written in a concise form, N = a + b x+ c y, i = 1, 2, 3 (7.8) i i i i “chap07” — 2002/12/14 — page 133 — 5134 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS We write the shape functions in the following matrix form:   a   i   T N = 1 xy b = p α (7.9) i i    c i T p  α where α is the vector of the three unknown constants, and p is the vector of polynomial basis functions (or monomials). Using Eq. (3.21), the moment matrix P corresponding to basis p can be given by   1 x y 1 1   P= 1 x y (7.10) 2 2 1 x y 3 3 Note that the above equation is written for the shape functions, and not for the dis- placements. For this particular problem, we use up to the first order of polynomial basis. Depending upon the problem, we could use higher order of polynomial basis functions. The complete order of polynomial basis functions in two-dimensional space up to the nth order can be given by using the so-called Pascal triangle, shown in Figure 3.2. The number of terms used in p depends upon the number of nodes the 2D element has. We usually try to use terms of lowest orders to make the basis as complete as possible in order. It is also possible to choose specific terms of higher orders for different types of elements. For our triangular element there are three nodes, and therefore the lowest terms with complete first order are used, as shown in Eq. (7.9). The assumption of Eq. (7.5) implies that the displacement is assumed to vary linearly in the element. In Eq. (7.8) there is a total of nine constants to be determined. Our task now is to determine these constants. If the shape functions constructed possess the delta function property, based on Lemmas 2 and 3 given in Chapter 3, the shape functions constructed will possess the partition of unity and linear field reproduction, as long as the moment matrix given in Eq. (7.10) is of full rank. Therefore, we can expect that the complete linear basis functions used in Eq. (7.9) guarantee that the shape functions to be constructed satisfy the sufficient requirements for FEM shape functions. What we need to do now is simply impose the delta function property on the assumed shape functions to determine the unknown constants a ,b i i and c . i The delta functions property states that the shape function must be a unit at its home node, and zero at its remote nodes. For a two-dimensional problem, it can be expressed as  1 for i = j N (x ,y )= (7.11) i j j 0 for i = j For a triangular element, this condition can be expressed explicitly for all three shape functions in the following equations. For shape function N ,wehave 1 N (x ,y )= 1 1 1 1 N (x ,y )= 0 (7.12) 1 2 2 N (x ,y )= 0 1 3 3 “chap07” — 2002/12/14 — page 134 — 67.2 LINEAR TRIANGULAR ELEMENTS 135 This is because node 1 at (x ,y ) is the home node of N , and nodes 2 at (x ,y )and3at 1 1 1 2 2 (x ,y ) are the remote nodes of N . Using Eqs. (7.5) and (7.12), we have 3 3 1 N (x ,y )= a + b x + c y = 1 1 1 1 1 1 1 1 1 N (x ,y )= a + b x + c y = 0 (7.13) 1 2 2 1 1 2 1 2 N (x ,y )= a + b x + c y = 0 1 3 3 1 1 3 1 3 Solving Eq. (7.13) for a ,b and c , we obtain 1 1 1 x y − x y y − y x − x 2 3 3 2 2 3 3 2 a = ,b = ,c = (7.14) 1 1 1 2A 2A 2A e e e where A is the area of the triangular element that can be calculated using the determinant e of the moment matrix:     1 x y 1 1   1 1 1   1 x y A = P= = (x y − x y )+ (y − y )x + (x − x )y (7.15) e 2 2 2 3 3 2 2 3 1 3 2 1   2 2 2   1 x y 3 3 Note here that as long as the area of the triangular element is nonzero, or as long as the three nodes are not on the same line, the moment matrix P will be of full rank. Substituting Eq. (7.14) into Eq. (7.5), we obtain 1 N = (x y − x y )+ (y − y )x+ (x − x )y (7.16) 1 2 3 3 2 2 3 3 2 2A e which can be re-written as 1 N = (y − y )(x− x )+ (x − x )(y− y ) (7.17) 1 2 3 2 3 2 2 2A e This equation clearly shows that N is a plane in the space of (x,y,N) that passes through 1 the line of 2–3, and vanishes at nodes 2 at (x ,y )and3at(x ,y ). This plane also passes the 2 2 3 3 point of (x ,y , 1) in space that guarantees the unity of the shape function at the home node. 1 1 Since the shape function varies linearly within the element, N can then be easily plotted 1 as in Figure 7.5(a). Making use of these features of N , we can immediately write out the 1 other two shape functions for nodes 2 and 3. For node 2, the conditions are N (x ,y )= 0 2 1 1 N (x ,y )= 1 (7.18) 2 2 2 N (x ,y )= 0 2 3 3 and the shape function N should pass through the line 3–1, which gives 2 1 N = (x y − x y )+ (y − y )x+ (x − x )y 2 3 1 1 3 3 1 1 3 2A e 1 = (y − y )(x− x )+ (x − x )(y− y ) (7.19) 3 1 3 1 3 3 2A e “chap07” — 2002/12/14 — page 135 — 7136 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS y 3 (x , y ) 3 3 (a) N (x,y) 1 (u , v ) 3 3 N (x , y )=1 1 1 1 1 N (x , y )=0 A 1 2 2 2(x , y ) 2 2 N (x , y )=0 1 (u , v ) 1 3 3 2 2 1 (x , y ) 1 1 (u , v ) 1 1 x Shape function N 1 y (b) 3 (x , y ) 3 3 N (x,y) 2 (u , v ) 3 3 N (x , y )=0 2 1 1 N (x , y )=1 2 2 2 1 N (x , y )=0 A 2 3 3 2 (x , y ) 2 2 (u , v ) 2 2 1 (x , y ) 1 1 (u , v ) 1 1 x Shape function N 2 (c) 1 y N (x,y) 3 3 (x , y ) 3 3 N (x , y )=0 (u , v ) 3 1 1 3 3 N (x , y )=0 3 2 2 A A N (x , y )=1 3 3 3 2 (x , y ) 2 2 (u , v ) 2 2 1 (x , y ) 1 1 (u , v ) 1 1 x Shape function N 3 Figure 7.5. Linear triangular element and its shape functions. which is plotted in Figure 7.5(b). For node 3, the conditions are N (x ,y )= 0 3 1 1 N (x ,y )= 0 (7.20) 3 2 2 N (x ,y )= 1 3 3 3 and the shape function N should pass through the line 1–2, and given by 3 1 N = (x y − x y )+ (y − y )x+ (x − x )y 3 1 2 1 1 1 2 2 1 2A e 1 = (y − y )(x− x )+ (x − x )(y− y ) (7.21) 1 2 1 2 1 1 2A e “chap07” — 2002/12/14 — page 136 — 87.2 LINEAR TRIANGULAR ELEMENTS 137 which is plotted in Figure 7.5(c). The process of determining these constants is basically sim- ple, algebraic manipulation. The shape functions are summarized in the following concise form: N = a + b x+ c y (7.22) i i i i 1 a = (x y − x y ) i j k k j 2A e 1 b = (y − y ) (7.23) i j k 2A e 1 c = (x − x ) i k j 2A e where the subscript i varies from 1 to 3, and j and k are determined by the cyclic permutation in the order of i,j,k. For example, if i = 1, then j = 2,k = 3. When i = 2, then j = 3,k = 1. 7.2.3 Area Coordinates Another alternative and effective method for creating shape functions for triangular elements is to use so-called area coordinates L ,L and L . The use of the area coordinates will 1 2 3 immediately lead to the shape functions for triangular elements. However, we first need to define the area coordinates. In defining L , we consider a point P at (x,y) inside the triangle, as shown in Figure 7.6, 1 and form a sub-triangle of 2–3–P. The area of this sub-triangle is noted as A , and can be 1 calculated using the formula     1 xy   1 1   1 x y A = = (x y − x y )+ (y − y )x+ (x − x )y (7.24) 1 2 2 2 3 3 2 2 3 3 2   2 2   1 x y 3 3 y A A A k, 3 1 2 3 L = , L = , L = 1 2 3 A A A e e e A 1 P i,1 j, 2 x Figure 7.6. Definition of area coordinates. “chap07” — 2002/12/14 — page 137 — 9138 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS The area coordinate L is then defined as 1 A 1 L = (7.25) 1 A e Similarly, for L we form sub-triangle 3–1–P with an area of A given by 2 2     1 xy   1 1   A = 1 x y = (x y − x y )+ (y − y )x+ (x − x )y (7.26) 2 3 3 3 1 1 3 3 1 1 3   2 2   1 x y 1 1 The area coordinate L is then defined as 2 A 2 L = (7.27) 2 A e For L , we naturally write 3 A 3 L = (7.28) 3 A e where A is the area of the sub-triangle 1–2–P, calculated using 3     1 xy   1 1   A = 1 x y = (x y − x y )+ (y − y )x+ (x − x )y (7.29) 3 1 1 1 2 2 1 1 2 2 1   2 2   1 x y 2 2 It is very easy to confirm the unity property of the area coordinates L ,L and L . First, 1 2 3 they are partitions of unity, i.e. L + L + L = 1 (7.30) 1 2 3 that can be proven using the definition of the area coordinates: A A A A + A + A 2 2 3 2 2 3 L + L + L = + + = = 1 (7.31) 1 2 3 A A A A e e e e Secondly, these area coordinates are of delta function properties. For example, L will 1 definitely be zero if P is at the remote nodes 2 and 3, and it will be a unit if P is at its home node 1. The same arguments are also valid for L and L . 2 3 These two properties are exactly those defined for shape functions. Therefore, we immediately have N = L,N = L,N = L (7.32) 1 1 2 2 3 3 The previous equation can also be easily confirmed by comparing Eqs. (7.16) with (7.25), (7.19) with (7.27), and (7.21) with (7.28). The area coordinates are very convenient for constructing higher order shape functions for triangular elements. Once the shape function matrix has been developed, one can write the displacement at any point in the element in terms of nodal displacements in the form of Eq. (7.1). The next step is to develop the strain matrix so that we can write the strain, and hence the stress, at any point in the element in terms of the nodal displacements. This will further lead to the element matrices. “chap07” — 2002/12/14 — page 138 — 107.2 LINEAR TRIANGULAR ELEMENTS 139 7.2.4 Strain Matrix Let us now move to the second step, which is to derive the strain matrix required for computing the stiffness matrix of the element. According to the discussion in Chapter 2, T there are only three major stress components, σ =σ σ σ in a 2D solid, and xx yy xy T the corresponding strains, ε =ε ε ε can be expressed as xx yy xy ∂u ε = xx ∂x ∂v ε = (7.33) yy ∂y ∂u ∂v ε = + xy ∂y ∂x or in a concise matrix form, ε = LU (7.34) where L is called a differential operation matrix, and can be obtained simply by inspection of Eq. (7.33):   ∂/∂x 0   L= 0 ∂/∂y (7.35) ∂/∂y ∂/∂x Substituting Eq. (7.1) into Eq. (7.34), we have ε = LU= LNd = Bd (7.36) e e where B is termed the strain matrix, which can be obtained by the following equation once the shape function is known:   ∂/∂x 0   0 ∂/∂y B= LN= N (7.37) ∂/∂y ∂/∂x Equation (7.36) implies that the strain is now expressed by the nodal displacement of the element using the strain matrix. Equations (7.36) and (7.37) are applicable for all types of 2D elements. Using Eqs. (7.3), (7.22), (7.23) and Eq. (7.37), the strain matrix B for the linear triangular element can be easily obtained, to have the following simple form:   a 0 a 0 a 0 1 2 3   B= 0 b 0 b 0 b (7.38) 1 2 3 b a b a b a 1 1 2 2 3 3 It can be clearly seen that the strain matrix B for a linear triangular element is a constant matrix. This implies that the strain within a linear triangular element is constant, and thus so is the stress. Therefore, the linear triangular elements are also referred to as constant strain “chap07” — 2002/12/14 — page 139 — 11140 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS elements or constant stress elements. In reality, stress or strain varies across the structure, hence using linear triangular elements with a coarse mesh will result in a rather inaccurate stress or strain distribution. We would need to have a fine mesh of linear triangular elements in order to show an appropriate variation of stress or strain across the structure. 7.2.5 Element Matrices Having obtained the shape function and the strain matrix, the displacement and strain (hence the stress) can all be expressed in terms of the nodal displacements of the element. The element matrices, like the stiffness matrix k , mass matrix m , and the nodal force vector e e f , can then be found using the equations developed in Chapter 3. e The element stiffness matrix k for 2D solid elements can be obtained using Eq. (3.71): e       h T T T k = B cB dV = dz B cB dA= hB cB dA (7.39) e V A 0 A e e e Note that the material constant matrix c has been given by Eqs. (2.31) and (2.32), respec- tively, for the plane stress and plane strain problems. Since the strain matrix B is a constant matrix, as shown in Eq. (7.38), and the thickness of the element is assumed to be uniform, the integration in Eq. (7.39) can be carried out very easily, which leads to T k = hA B cB (7.40) e e The element mass matrix m can also be easily obtained by substituting the shape e function matrix into Eq. (3.75):     h T T T m = ρN N dV = dxρN N dA= hρN N dA (7.41) e V A 0 A e e e For elements with uniform thickness and density, we can rewrite Eq. (7.41) as   N N 0 N N 0 N N 0 1 1 1 2 1 3   0 N N 0 N N 0 N N 1 1 1 2 1 3      N N 0 N N 0 N N 0 2 1 2 2 2 3   m = hρ dA (7.42) e   0 N N 0 N N 0 N N 2 1 2 2 2 3 A e    N N 0 N N 0 N N 0 3 1 3 2 3 3 0 N N 0 N N 0 N N 3 1 3 2 3 3 The integration of all the terms in the mass matrix can be carried out by simply using a mathematical formula developed by Eisenberg and Malvern (1973):  mnp p m n L L L dA= 2A (7.43) 1 2 3 (m+ n+ p+ 2) A where L = N is the area coordinates for triangular elements that is the same as i i the shape function, as we have seen in Section 7.2.2. The element mass matrix m is e “chap07” — 2002/12/14 — page 140 — 127.3 LINEAR RECTANGULAR ELEMENTS 141 found to be   2 0 10 10   20 1 0 1     ρhA 2010   m = (7.44) e   201 12     sy. 20 2 The nodal force vector for 2D 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–3ofthe s element, as shown in Figure 7.4; the nodal force vector becomes      f sx T f = N  dl (7.45) e f 2–3 sy l If the load is uniformly distributed, f and f are constants within the element, so the sx sy above equation becomes   0         0       1 f x xf = l (7.46) e 2–3 f   2 y       f   x     f y where l is the length of the edge 2–3 of the element. 2–3 Once the element stiffness matrix k , mass matrix m and nodal force vector f have e e e been obtained, the global finite element equation can be obtained by assembling the element matrices by summing up the contribution from all the adjacent elements at the shared nodes. 7.3 LINEAR RECTANGULAR ELEMENTS Triangular elements are usually not preferred by many analysts nowadays, unless there are difficulties with the meshing and re-meshing of models of complex geometry. The main reason is that the triangular elements are usually less accurate than rectangular or quadrilateral elements. As shown in the previous section, the strain matrix of the linear triangular elements is constant, accounting for the inaccuracy. With advances in meshing algorithms, many models of complex geometry with sharp corners or curved edges can be modelled using quadrilateral elements. For the rectangular element, the strain matrix is not a constant, as will be shown in this section. This will provide a more realistic presentation in strain, and hence the stress distribution, across the structure. The formulation of the equations for the rectangular elements is simpler compared to the triangular elements, because the shape functions can form very easily due to the regularity in the shape of the rectangular element. The simple three-step procedure is applicable, and will be shown in the following sections. “chap07” — 2002/12/14 — page 141 — 13142 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS 7.3.1 Shape Function Construction Consider a 2D domain. The domain is discretized into a number of rectangular elements with four nodes and four straight edges, as shown in Figure 7.7. As always, we number the nodes in each element 1, 2, 3 and 4 in a counter-clockwise direction. Note also that, since each node has two DOFs, the total DOFs for a linear rectangular element would be eight. The dimension of the element is defined here as 2a×2b×h. A local natural coordinate system (ξ,η) with its origin located at the centre of the rectangular element is defined. The relationship between the physical coordinate (x,y) and the local natural coordinate system (ξ,η)isgivenby ξ = x/a, η= y/b (7.47) Equation (7.47) defines a very simple coordinate mapping between physical and natural coordinate systems for rectangular elements as shown in Figure 7.8. Our formulation can now be based on the natural coordinate system. The use of natural coordinates will make the construction of the shape functions and evaluation of the matrix integrations very much easier. This kind of coordinate mapping technique is one of the most often used techniques in the FEM. It is extremely powerful when used for developing elements of complex shapes. We perform the field variable interpolation and express the displacement within the ele- ment as an interpolation of the nodal displacements using shape functions. The displacement vector U is assumed to have the form h U (x, y)= N(x, y)d (7.48) e where the nodal displacement vector d is arranged in the form e    u   1   displacements at node 1     v   1       u   2   displacements at node 2   u 2  d = (7.49) e u   3     displacements at node 3   u   3        u   4     displacements at node 4 u 4 Figure 7.7. Rectangular domain meshed by rectangular elements. “chap07” — 2002/12/14 — page 142 — 147.3 LINEAR RECTANGULAR ELEMENTS 143 and the matrix of shape functions has the form N 0 N 0 N 0 N 0 1 2 3 4 N= 0 N 0 N 0 N 0 N 1 2 3 4 (7.50)     Node 1 Node 2 Node 3 Node 4 where the shape functions N (i = 1, 2, 3, 4) are the shape functions corresponding to the i four nodes of the rectangular element. In determining these shape functions N (i = 1, 2, 3, 4), we can follow exactly the same i steps used in Sections 4.2.1 or 7.2.2, by starting with an assumption of the displacement or shape functions using polynomial basis functions with unknown constants. These unknown constants are then determined using the displacements at the nodes of the element or the property of the shape functions. The only difference here is that we need to use four terms of monomials of basis functions. As we have seen in Section 7.2.2, the process is quite troublesome and lengthy. For many cases, one often constructs shape functions simply by some ‘shortcut’ methods. One of these is by inspection, and utilizing the properties of shape functions. Due to the regularity of the square element in the natural coordinates, the shape functions in Eq. (7.50) can be written out directly as follows, without going through the detailed process that we described in the previous section for triangular elements: 1 N = (1− ξ)(1− η) 1 4 1 N = (1+ ξ)(1− η) 2 4 (7.51) 1 N = (1+ ξ)(1+ η) 3 4 1 N = (1− ξ)(1+ η) 4 4 3 (1, +1) 4 (−1, +1) (a) (b) (u , v ) (u , v ) 4 4 3 3 y, v η 4 (x , y ) 3 (x , y ) 4 4 3 3 η (u , v ) (u , v ) 4 4 3 3 ξ f sy 2b 2b f sx ξ 2a 2 (x , y ) 1 (x , y ) 2 2 1 1 (u , v ) (u , v ) 2a 1 1 2 2 2 (1, −1) 1 (−1, −1) x, u (u , v ) (u , v ) 2 2 1 1 Figure 7.8. Rectangular element and the coordinate systems. (a) Rectangular element in physical system, (b) square element in natural coordinate system. “chap07” — 2002/12/14 — page 143 — 15144 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS It is very easy to confirm that all the shape functions given in Eq. (7.51) satisfy the delta function property of Eq. (3.34). For example, for N we have 3   1 N = (1+ ξ)(1+ η) = 0 3 at node 1 4 ξ=−1 η=−1   1 N = (1+ ξ)(1+ η) = 0 3 at node 2 4 ξ=1 η=−1 (7.52)   1 N = (1+ ξ)(1+ η) = 1  3 at node 3 ξ=1 4 η=1   1 N = (1+ ξ)(1+ η) = 0  3 at node 4 4 ξ=−1 η=1 The same examination of N ,N and N will confirm the same property. 1 2 4 It is also very easy to confirm that all the shape functions given in Eq. (7.51) satisfy the partition of unity property of Eq. (3.41): 4  N = N + N + N + N i 1 2 3 4 i=1 1 = (1− ξ)(1− η)+ (1+ ξ)(1− η)+ (1+ ξ)(1+ η)+ (1− ξ)(1+ η) 4 1 = 2(1− ξ)+ 2(1+ ξ)= 1 (7.53) 4 The partitions of unity property can also be easily confirmed using Lemma 1 in Chapter 3. Equation (7.53) should be called a bilinear shape function to be exact, as it varies linearly in both the ξ and η directions. It varies quadratically in any direction other than these two ξ and η directions. Denoting the natural coordinates of node j by (ξ ,η ), the bilinear shape j j function N can be re-written in a concise form: j 1 N = (1+ ξ ξ)(1+ η η) (7.54) j j j 4 7.3.2 Strain Matrix Using the same procedure as for the case of the triangular element, the strain matrix B would have the same form as in Eq. (7.37), that is B= LN   1− η 1− η 1+ η 1+ η − 0 0 0 − 0   a a a a   1− ξ 1+ ξ 1+ ξ 1− ξ   =  0 − 0 − 0 0    b b b b   1− ξ 1− η 1+ ξ 1− η 1+ ξ 1+ η 1− ξ 1+ η − − − − b a b a b a b a (7.55) “chap07” — 2002/12/14 — page 144 — 167.3 LINEAR RECTANGULAR ELEMENTS 145 It is now clear that the strain matrix for a bilinear rectangular element is no longer a constant matrix. This implies that the strain, and hence the stress, within a linear rectangular element is not constant. 7.3.3 Element Matrices Having obtained the shape function and the strain matrix B, the element stiffness matrix k , e mass matrix m , and the nodal force vector f can be obtained using the equations presented e e in Chapter 3. Using first the relationship given in Eq. (7.47), we have dx dy = ab dξ dη (7.56) Substituting Eq. (7.56) into Eq. (7.39), we obtain    +1 +1 T T k = hB cB dA= abhB cB dξ dη (7.57) e A −1 −1 The material constant matrix c has been given by Eqs. (2.31) and (2.32), respectively, for plane stress and plane strain problems. Evaluation of the integral in Eq. (7.57) would not be as straightforward, since the strain matrix B is a function of ξ and η. It is still possible to obtain the closed form for the stiffness matrix by carrying out the integrals in Eq. (7.57) analytically. In practice, we often use a numerical integration scheme to evaluate the integral, and the commonly used Gauss integration scheme will be introduced here. The Gauss integration scheme is a very simple and efficient procedure that performs numerical integral, and it is briefly outlined here. 7.3.4 Gauss Integration Consider first a one-dimensional integral. Using the Gauss integration scheme, the integral is evaluated simply by a summation of the integrand evaluated at m Gauss points multiplied by corresponding weight coefficients as follows:  m +1  I = f(ξ) dξ = w f(ξ ) (7.58) j j −1 j=1 The locations of the Gauss points and the weight coefficients have been found for different m, and are given in Table 7.1. In general, the use of more Gauss points will produce more accurate results for the integration. However, excessive use of Gauss points will increase the computational time and use up more computational resources, and it may not necessarily give better results. The appropriate number of Gauss points to be used depends upon the complexity of the integrand. It has been proven that the use of m Gauss points gives the exact results of a polynomial integrand of up to an order of n = 2m− 1. For example, if the integrand is a linear function (straight line), we have 2m− 1 = 1, which gives m = 1. This means that for a linear integrand, one Gauss point will be sufficient to “chap07” — 2002/12/14 — page 145 — 17146 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS Table 7.1. Gauss integration points and weight coefficients mξ w Accuracy n j j 10 2 1 √ √ 2 −1/ 3, 1/ 3 1,1 3 √ √ 3 − 0.6, 0, 0.65/9, 8/9, 5/95 4 −0.861136,−0.339981, 0.347855, 0.652145, 7 0.339981, 0.861136 0.652145, 0.347855 5 −0.906180,−0.538469, 0, 0.236927, 0.478629, 0.568889, 9 0.538469, 0.906180 0.478629, 0.236927 6 −0.932470,−0.661209,−0.238619, 0.171324, 0.360762, 0.467914, 11 0.238619, 0.661209, 0.932470 0.467914, 0.360762, 0.171324 give the exact result of the integration. If the integrand is of a polynomial of a third order, we have 2m− 1= 3, which gives m= 2. This means that for an integrand of a third order polynomial, the use of two Gauss points will be sufficient to give the exact result. The use of more than two points will still give the same results, but takes more computation time. For two-dimensional integrations, the Gauss integration is sampled in two directions, as follows:   n n y x +1 +1  I = f(ξ,η) dξ dη= w w f(ξ ,η ) (7.59) i j i j −1 −1 i=1 j=1 Figure 7.9(b) shows the locations of four Gauss points used for integration in a square region. The element stiffness matrix k can be obtained by numerically carrying out the integrals e in Eq. (7.57) using the Gauss integration scheme of Eq. (7.59). 2× 2 Gauss points shown in Figure 7.9(b) are sufficient to obtain the exact solution for the stiffness matrix given by Eq. (7.57). This is because the entry in the strain matrix, B, is a linear function of ξ or η. T The integrand in Eq. (7.57) consists of B cB, which implies multiplications of two linear functions, and hence this becomes a quadratic function. In Table 7.1, having two Gauss points sampled in each direction is sufficient to obtain the exact results for a polynomial function in that direction of order up to 3. Figure 7.9(a) and (c) show some other different and possible number of integration points in a square region. To obtain the element mass matrix m , we substitute Eq. (7.56) into Eq. (3.75) to obtain e     h T T T m = ρN N dV = dxρN N dA= hρN N dA e V A 0 A   +1 +1 T = abhρN N dξ dη (7.60) −1 −1 “chap07” — 2002/12/14 — page 146 — 187.3 LINEAR RECTANGULAR ELEMENTS 147 (b) η η (c) η (a) 4(–1,+1) 4(–1,+1) 3(1,+1) 3(1,+1) 4(–1,+1) 3(1,+1) ξ ξ ξ 1(–1, –1) 2(1, –1) 1(–1, –1) 2(1, –1) 1(–1, –1) 2(1, –1) Figure 7.9. Integration points for n = n = 1, 2 and 3 in a square region. x y Upon evaluation of the integral, after substitution of Eq. (7.50) into Eq. (7.60), the element mass matrix m is obtained explicitly as e   4 0 201020   40 2 0 1 0 2     402010     ρhab 40201   m = (7.61) e   4020 9     402     sy. 40 4 In obtaining element m in the mass matrix, the following integral has been carried out and ij repeatedly used:   +1 +1 m = ρhab N N dξ dη ij i j −1 −1   +1 +1 ρhab = (1+ ξ ξ)(1+ ξ ξ) dξ (1+ η η)(1+ η η) dη i j i j 16 −1 −1    ρhab 1 1 = 1+ ξ ξ 1+ η η (7.62) i j i j 4 3 3 For example, in calculating m , we use the above equation to obtain 33    ρhab 1 1 ρhab m = 1+ × 1× 1 1+ × 1× 1 = 4× (7.63) 33 4 3 3 9 In practice, the integrals in Eq. (7.60) are often calculated numerically using the Gauss integration scheme. “chap07” — 2002/12/14 — page 147 — 19148 CHAPTER 7 FEM FOR TWO-DIMENSIONAL SOLIDS The nodal force vector for a rectangular element can be obtained by using Eqs. (3.78), (3.79) and (3.81). Suppose the element is loaded by a distributed force f s on edge 2–3 of the element, as shown in Figure 7.8; the nodal force vector becomes     f  T sx f = N  dl (7.64) e f 2-3 sy l If the load is uniformly distributed within the element, and f and f are constant, the sx sy above equation becomes   0         0         f  x     f y f = b (7.65) e f  x       f  y        0     0 where b is the half length of the side 2–3. Equation (7.65) suggests that the evenly distributed load is divided equally onto nodes 2 and 3. The stiffness matrix k , mass matrix m and nodal force vector f can be used directly e e e to assemble the global FE equation, Eq. (3.96). Coordinate transformation is needed if the orientation of the local natural coordinate does not coincide with that of the global coordinate system. In such a case, quadrilateral elements are often used, which is to be developed in the next section. 7.4 LINEAR QUADRILATERAL ELEMENTS Though the rectangular element can be very useful, and is usually more accurate than the triangular element, it is difficult to use it for problems with any geometry rather than rectangles. Hence, its practical application is very limited. A much more practical and useful element would be the so-called quadrilateral element, that can have unparalleled edges. However, there can be a problem for the integration of the mass and stiffness matrices for a quadrilateral element, because of the irregular shape of the integration domain. The Gauss integration scheme cannot be implemented directly with quadrilateral elements. Therefore, what is required first is to map the quadrilateral element into the natural coordinates system to become a square element, so that the shape functions and the integration method used for the rectangular element can be utilized. Hence, key in the development of a quadrilateral element is the coordinate mapping. Once the mapping is established, the rest of the procedure is exactly the same as that used for formulating the rectangular element in the previous section. 7.4.1 Coordinate Mapping Figure 7.10 shows a 2D domain with the shape of an airplane wing. As you can imagine, dividing such a domain into rectangular elements of parallel edges is impossible. The job “chap07” — 2002/12/14 — page 148 — 20