FEM FOR 3D SOLIDS

FEM FOR 3D SOLIDS
Dr.NaveenBansal Profile Pic
Dr.NaveenBansal,India,Teacher
Published Date:25-10-2017
Your Website URL(Optional)
Comment
9 FEM FOR 3D SOLIDS 9.1 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 ofx,y andz. 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 f f f f Figure 9.1. Example of a 3D solid under loadings. 199 “chap09” — 2002/12/14 — page 199 — 1200 CHAPTER 9 FEM FOR 3D SOLIDS 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 books (see, e.g., Washizu, 1981; Rao, 1999; Zienkiewicz and Taylor, 2000; etc.). 9.2 TETRAHEDRON ELEMENT 9.2.1 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 Figure 9.2. Solid block divided into four-node tetrahedron elements. w 4 4= l v 4 w 3 u 4 w 1 v 3 1= i k 3= u 3 v 1 f sz u z=Z 1 f w sy 2 f sx y=Y u 2 2= j x=X u 2 Figure 9.3. A tetrahedron element. “chap09” — 2002/12/14 — page 200 — 29.2 TETRAHEDRON ELEMENT 201 (u, v andw), making the total DOFs in a tetrahedron element twelve, as shown in Figure 9.3. 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: h U (x,y,z)= N(x,y,z)d (9.1) e where the nodal displacement vector, d ,isgivenas e   u   1       v   displacements at node 1 1       w   1       u   2       v displacements at node 2   2     w 2  d = (9.2) e u   3       v displacements at node 3   3       w   3       u   4       v displacements at node 4   4     w 4 and the matrix of shape functions has the form node 1 node 2 node 3 node 4       N 00 N 00 N 00 N 00 1 2 3 4 (9.3)   0 N 0 0 N 0 0 N 0 0 N 0 1 2 3 4 N= 00 N 00 N 00 N 00 N 1 2 3 4 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 V P 234 L = (9.4) 1 V 1234 whereV andV denote, respectively, the volumes of the tetrahedrons P234 and 1234, P234 1234 as shown in Figure 9.4. The volume coordinate for node 2-4 can also be defined in the same “chap09” — 2002/12/17 — page 201 — 3202 CHAPTER 9 FEM FOR 3D SOLIDS 4= l 1= i 3= k P z y x 2= j Figure 9.4. Volume coordinates for tetrahedron elements. manner: V V V P 134 P 124 P 123 L = ,L = ,L = (9.5) 2 3 4 V V V 1234 1234 1234 The volume coordinate can also be viewed as the ratio between the distance of the point P and point 1 to the plane 234: d d d d P−234 P−134 P−124 P−123 L = ,L = ,L = ,L = (9.6) 1 2 3 4 d d d d 1−234 1−234 1−234 1−234 It can easily be confirmed that L +L +L +L = 1 (9.7) 1 2 3 4 since V +V +V +V =V (9.8) P 234 P 134 P 124 P 123 1234 It can also easily be confirmed that  1 at the home nodei L = (9.9) i 0 at the remote nodesjkl Using Eq. (9.9), the relationship between the volume coordinates and Cartesian coordinates can be easily derived: x =L x +L x +L x +L x 1 1 2 2 3 3 4 4 y =L y +L y +L y +L y (9.10) 1 1 2 2 3 3 4 4 z=L z +L z +L z +L z 1 1 2 2 3 3 4 4 Equations (9.7) and (9.10) can then be expressed as a single matrix equation as follows:      1 1111 L     1           x x x x x L 1 2 3 4 2   = (9.11)   y y y y y L   1 2 3 4  3         z z z z z L 1 2 3 4 4 “chap09” — 2002/12/17 — page 202 — 49.2 TETRAHEDRON ELEMENT 203 The inversion of Eq. (9.11) will give      L a b c d 1     1 1 1 1 1           1 L a b c d x 2 2 2 2 2   = (9.12)   L a b c d y     3 6V 3 3 3 3         L a b c d z 4 4 4 4 4 where     x y z 1 y z j j j j j     a = det x y z ,b =− det 1 y z i k k k i k k x y z 1 y z l l l l l (9.13)     y 1 z y z 1 j j j j     c =− det y 1 z ,d =− det y z 1 i k k i k k y 1 z y z 1 l 1 l l 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. Wheni = 2, thenj = 3,k = 4,l = 1. The volume of the tetrahedron elementV can be obtained by   1 x y z i i i   1 1 x y z j j j   V = × det (9.14)   1 x y z 6 k k k 1 x y z l l l The properties of L , as depicted in Eqs. (9.6) to (9.9), show that L can be used as the i i shape function of a four-nodal tetrahedron element: 1 N =L = (a +b x+c y+d z) (9.15) i i i i i i 6V It can be seen from above that the shape function is a linear function ofx,y andz, 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, xx yy zz yz xz xy ε ε ε ε ε ε , we can substitute Eq. (9.1) into Eq. (2.5): xx yy zz yz xz xy ε = LU= LNd = Bd (9.16) e e “chap09” — 2002/12/17 — page 203 — 5204 CHAPTER 9 FEM FOR 3D SOLIDS where the strain matrix B is given by   ∂/∂x 00   0 ∂/∂y 0     00 ∂/∂z   B= LN= N (9.17)   0 ∂/∂z ∂/∂y     ∂/∂z 0 ∂/∂x ∂/∂y ∂/∂x 0 Using Eq. (9.3), the strain matrix, B, can be obtained as   b 00 b 00 b 00 b 00 1 2 3 4   0 c 00 c 00 c 00 c 0 1 2 3 4     1 00 d 00 d 00 d 00 d 1 2 3 4   B= (9.18)   c b 0 c b 0 c b 0 c b 0 2V 1 1 2 2 3 3 4 4     0 d c 0 d c 0 d c 0 d c 1 1 2 2 3 3 4 4 d 0 b d 0 b d 0 b d 0 b 1 1 2 2 3 3 4 4 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. 9.2.2 Element Matrices Once the strain matrix has been obtained, the stiffness matrix k for 3D solid elements e can be obtained by substituting Eq. (9.18) into Eq. (3.71). Since the strain is constant, the element strain matrix is obtained as  T T k = B cB dV =V B cB (9.19) e e V e 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):   N N N N 11 12 13 14     N N N N T 21 22 23 24   m = ρN N dV = ρ dV (9.20) e   N N N N 31 32 33 34 V V e e N N N N 41 42 43 44 where   N N 00 i j   N = 0 N N 0 (9.21) ij i j 00 N N i j Using the following formula Eisenberg and Malvern, 1973,  mnpq p q m n L L L L dV = 6V (9.22) e 1 2 3 4 (m+n+p+q+ 3) V e “chap09” — 2002/12/17 — page 204 — 69.2 TETRAHEDRON ELEMENT 205 we can conveniently evaluate the integral in Eq. (9.20) to give   20 0 1 00100100   2 0 01 001 0 01 0     2 0 01 001 0 01     200100100     20010010     ρV 2001001 e   m = (9.23) e   200100 20     20010     2001     sy. 200     20 2 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–4ofthe 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 4 = l  = constant 3 = k =1 1 = i Q =0 z P j 2 = y =1 x Figure 9.5. Natural coordinate, whereξ = constant. “chap09” — 2002/12/14 — page 205 — 7206 CHAPTER 9 FEM FOR 3D SOLIDS 4= l  = constant k 3= 1= i =1 =0 P z j 2 = y =0 x Figure 9.6. Natural coordinate, whereη= constant. 4= l =0  = constant R Q k 3 = 1= i =1 =1 P z 2 = j y =1 x Figure 9.7. Natural coordinate, whereζ = constant. plane 1–2–3 on the element sits on thex–y plane. 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: x =η(x −x )+x P 3 2 2 y =η(y −y )+y (9.24) P 3 2 2 z = 0 P The coordinates at point B are then interpolated using thex,y andz coordinates at points 1 and P: x =ξ(x −x )+x =ξη(x −x )+ξ(x −x )+x B P 1 1 3 2 2 1 1 y =ξ(y −y )+y =ξη(y −y )+ξ(y −y )+y (9.25) B P 1 1 3 2 2 1 1 z = 0 B “chap09” — 2002/12/14 — page 206 — 89.2 TETRAHEDRON ELEMENT 207 4 = l =0  = constant k O 3 =  = 1 1= i  = 1 =0  = 1 B =0 Px = (x − x )+x , y = ( y − y )+y ,0 P 3 2 2 P 3 2 2 =1 z Bx = (x − x )+x , y = ( y − y ) + y ,0 B P 1 1 B P 1 1 ξ=constant y =1 j 2 =  = constant =0 =1 x O x=(1− ζ)(x − x ) + x , y=(1−)( y − y ) + y , z =(1− ) z 4 B B 4 B B 4 Figure 9.8. Cartesian coordinates xyz of point O in term ofξηζ. The coordinates at point O are finally interpolated using thex,y andz coordinates at points 4 and B: x =x −ζ(x −x )=x −ζ(x −x )+ξζ(x −x )−ξζ(x −x ) 4 4 B 4 4 1 2 1 2 3 y =y −ζ(y −y )=y −ζ(y −y )+ξζ(y −y )−ξζ(y −y ) (9.26) 4 4 B 4 4 1 2 1 2 3 z=(1−ζ)z 4 With this special natural coordinate system, the shape functions in the matrix of Eq. (9.3) can be written by inspection as N =(1−ξ)ζ 1 N =ξηζ 2 (9.27) N =ξζ(1−η) 3 N =(1−ζ) 4 The Jacobian matrix between xyz andξηζ is required, and is given as   ∂x/∂ξ ∂x/∂η ∂x/∂ζ   J= ∂y/∂ξ ∂y/∂η ∂y/∂ζ (9.28) ∂z/∂ξ ∂z/∂η ∂z/∂ζ Using Eqs. (9.26) and (9.27), the determinate of the Jacobian can be found to be     ζx +ηζx ξζx −x +ξx +ξηx 21 31 31 41 21 31   2   detJ= ζy +ηζy ξζy −y +ξy +ξηy =−6Vξζ (9.29) 21 31 31 41 21 31     0 z 0 4 “chap09” — 2002/12/14 — page 207 — 9208 CHAPTER 9 FEM FOR 3D SOLIDS The mass matrix can now be obtained as     1 1 1 T T m = ρN N dV = ρN N detJ dξ dη dζ (9.30) e V 0 0 0 e which gives   N N N N 11 12 13 14    1 1 1   N N N N 21 22 23 24 2   m =−6V ρ ξζ dξ dη dζ (9.31) e e   N N N N 31 32 33 34 0 0 0 N N N N 41 42 43 44 where N is given by Eq. (9.21), but in which the shape functions should be defined by ij 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–3ofthe s element as shown in Figure 9.3; the nodal force vector becomes     f   sx  T  f f = N dl (9.32) e sy    l 3–4 f sz If the load is uniformly distributed,f ,f andf are constants, and the above equation sx sx sz becomes   0  3×1       0  3×1       f  sx        f  sy        f  sz     1 f  sx f = l (9.33) e 3–4  f  2 sy        f  sz       0   3×1       0   3×1       0   3×1     0 3×1 wherel is the length of the edge 3–4. Equation (9.33) implies that the distributed forces 3–4 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, k , the mass e matrix, m , and the nodal force vector, f , can be used directly to assemble the global FE e e equation, Eq. (3.96), without going through a coordinate transformation. “chap09” — 2002/12/14 — page 208 — 109.3 HEXAHEDRON ELEMENT 209 9.3 HEXAHEDRON ELEMENT 9.3.1 Strain Matrix Consider now a 3D domain, which is divided in a proper manner into a number of hexahedron elements with eight nodes and six surfaces, as shown in Figure 9.9. Each hexahedron element has nodes numbered 1, 2, 3, 4 and 5, 6, 7, 8 in a counter-clockwise manner, as shown in Figure 9.10. As there are three DOFs at one node, there is a total of 24 DOFs in a hexahedron element. It is again useful to define a natural coordinate system (ξ,η,ζ) with its ori- gin at the centre of the transformed cube, as this makes it easier to construct the shape functions and to evaluate the matrix integration. The coordinate mapping is preformed in a similar manner as for quadrilateral elements in Chapter 7. Like the quadrilateral element, shape functions are also used to interpolate the coordinates from the nodal P P P P Figure 9.9. Solid block divided into eight-nodal hexahedron elements.  (–1, –1, 1) 5 8(–1, 1, 1) 5 8 (1, –1, 1)6 f 7(1, 1, 1) 6 sz 4  7 1 (–1, –1, –1)1 4(–1, 1, –1) 0 z f sy f sx 2 (1, –1,–1) 2 0 y 3 (1, 1, –1) 3  x Figure 9.10. An eight-nodal hexahedron element and the coordinate systems. “chap09” — 2002/12/14 — page 209 — 11210 CHAPTER 9 FEM FOR 3D SOLIDS coordinates: 8  x = N(ξ,η,ζ)x i i i=1 8  y = N(ξ,η,ζ)y (9.34) i i i=1 8  z= N(ξ,η,ζ)z i i i=1 The shape functions are given in the local natural coordinate system as 1 N = (1−ξ)(1−η)(1−ζ) 1 8 1 N = (1+ξ)(1−η)(1−ζ) 2 8 1 N = (1+ξ)(1+η)(1−ζ) 3 8 1 N = (1−ξ)(1+η)(1−ζ) 4 8 (9.35) 1 N = (1−ξ)(1−η)(1+ζ) 5 8 1 N = (1+ξ)(1−η)(1+ζ) 6 8 1 N = (1+ξ)(1+η)(1+ζ) 7 8 1 N = (1−ξ)(1+η)(1+ζ) 8 8 or in a concise form, 1 N = (1+ξξ )(1+ηη )(1+ζζ ) (9.36) i i i i 8 where(ξ,η,ζ ) denotes the natural coordinates of nodeI. i i i From Eq. (9.36), it can be seen that the shape functions vary linearly in theξ,η and ζ directions. Therefore, these shape functions are sometimes called tri-linear functions. The shape functionN is a three-dimensional analogy of that given in Eq. (7.54). It is very i easy to directly observe that the tri-linear elements possess the delta function property. In addition, since all these shape functions can be formed using the common set of eight basis functions of 1,ξ,η,ς,ξη,ξς,ης,ξης (9.37) which contain both constant and linear basis functions. Therefore, these shape functions can expect to possess both partitions of the unity property as well as the linear reproduction property (see Lemmas 2 and 3 in Chapter 3). “chap09” — 2002/12/14 — page 210 — 129.3 HEXAHEDRON ELEMENT 211 In a hexahedron element, the displacement vector U is a function of the coordinatesx, y andz, and as before, it is interpolated using the shape functions U= Nd (9.38) e where the nodal displacement vector, d is given by e    d  displacement components at node 1 e1       d displacement components at node 2   e2       d displacement components at node 3   e3     d displacement components at node 4 e4 d = (9.39) e d displacement components at node 5   e5       d displacement components at node 6   e6       d displacement components at node 7   e7     d displacement components at node 8 e8 in which   u   1 d = v (i = 1, 2,..., 8) (9.40) ei 1   w 1 is the displacement at nodei. The matrix of shape functions is given by   N= N N N N N N N N (9.41) 1 2 3 4 5 6 7 8 in which each sub-matrix, N ,isgivenas i   N 00 i   N = 0 N 0 (i = 1, 2,..., 8) (9.42) i i 00 N i In this case, the strain matrix defined by Eq. (9.17) can be expressed as   B= B B B B B B B B (9.43) 1 2 3 4 5 6 7 8 whereby   ∂N /∂x 00 i   0 ∂N /∂y 0 i     00 ∂N /∂z i   B = LN = (9.44) i i   0 ∂N /∂z ∂N /∂y i i     ∂N /∂z 0 ∂N /∂x i i ∂N /∂y ∂N /∂x 0 i i As the shape functions are defined in terms of the natural coordinates,ξ,η andζ , to obtain the derivatives with respect to x, y and z in the strain matrix, the chain rule of partial “chap09” — 2002/12/14 — page 211 — 13212 CHAPTER 9 FEM FOR 3D SOLIDS differentiation needs to be used: ∂N ∂N ∂x ∂N ∂y ∂N ∂z i i i i = + + ∂ξ ∂x ∂ξ ∂y ∂ξ ∂z ∂ξ ∂N ∂N ∂x ∂N ∂y ∂N ∂z i i i i = + + (9.45) ∂η ∂x ∂η ∂y ∂η ∂z ∂η ∂N ∂N ∂x ∂N ∂y ∂N ∂z i i i i = + + ∂ζ ∂x ∂ζ ∂y ∂ζ ∂z ∂ζ which can be expressed in the matrix form     ∂N /∂ξ ∂N /∂x i i ∂N /∂η = J ∂N /∂y (9.46) i i     ∂N /∂ζ ∂N /∂z i i where J is the Jacobian matrix defined by   ∂x/∂ξ ∂y/∂ξ ∂z/∂ξ   J= ∂x/∂η ∂y/∂η ∂z/∂η (9.47) ∂x/∂ζ ∂y/∂ζ ∂z/∂ζ Recall that the coordinates, x, y and z are interpolated by the shape functions from the nodal coordinates. Hence, substitute the interpolation of the coordinates, Eq. (9.34), into Eq. (9.47), which gives   x y z 1 1 1   ∂N ∂N ∂N ∂N ∂N ∂N ∂N ∂N   1 2 3 4 5 6 7 8 x y z 2 2 2      ∂ξ ∂ξ ∂ξ ∂ξ ∂ξ ∂ξ ∂ξ ∂ξ x y z 3 3 3       ∂N ∂N ∂N ∂N ∂N ∂N ∂N ∂N 1 2 3 4 7 8 x y z 5 6 4 4 4    J= (9.48)    x y z ∂η ∂η ∂η ∂η ∂η ∂η ∂η ∂η 5 5 5       ∂N ∂N ∂N ∂N ∂N ∂N ∂N ∂N x y z 6 6 6 1 2 3 4 5 6 7 8     x y z 7 7 7 ∂ζ ∂ζ ∂ζ ∂ζ ∂ζ ∂ζ ∂ζ ∂ζ x y z 8 8 8 or     8 8 8 x∂N /∂ξ y∂N /∂ξ z∂N /∂ξ i i i i i i i=1 i=1 i=1      8 8 8   J= (9.49) x∂N /∂η y∂N /∂η z∂N /∂η i i i i i i  i=1 i=1 i=1     8 8 8 x∂N /∂ζ y∂N /∂ζ z∂N /∂ζ i i i i i i i=1 i=1 i=1 Equation (9.46) can be re-written as     ∂N /∂x ∂N /∂ξ     i i −1 ∂N /∂y = J ∂N /∂η (9.50) i i     ∂N /∂z ∂N /∂ζ i i which is then used to compute the strain matrix, B, in Eqs. (9.43) and (9.44), by replacing all the derivatives of the shape functions with respect tox,y andz to those with respect to ξ,η andζ . “chap09” — 2002/12/14 — page 212 — 149.3 HEXAHEDRON ELEMENT 213 9.3.2 Element Matrices Once the strain matrix, B, has been computed, the stiffness matrix, k , for 3D solid elements e can be obtained by substituting B into Eq. (3.71):     +1 +1 +1 T T k = B cB dA= B cB detJ dξ dη dζ (9.51) e V −1 −1 −1 e Note that the matrix of material constant, c, is given by Eq. (2.9). As the strain matrix, B, is a function ofξ,η andζ , evaluating the integrations in Eq. (9.51) can be very difficult. Therefore, the integrals are performed using a numerical integration scheme. The Gauss integration scheme discussed in Section 7.3.4 is often used to carry out the integral. For three-dimensional integrations, the Gauss integration is sampled in three directions, as follows:    n m l +1 +1 +1  I = f(ξ,η) dξ dη= w w w f(ξ,η ,ζ ) (9.52) i j k i j j −1 −1 −1 i=1 j=1k=1 To obtain the mass (inertia) matrix for the hexahedron element, substitute the shape function matrix, Eq. (9.41), into Eq. (3.75):     1 1 1 T T m = ρN N dV = ρN N detJ dξ dη dζ (9.53) e V −1 −1 −1 e The above integral is also usually carried out using Gauss integration. If the hexahedron is rectangular with dimensions ofa×b×c, the determinate of the Jacobian matrix is simply given by detJ=abc=V (9.54) e and the mass matrix can be explicitly obtained as   m m m m m m m m 11 12 13 14 15 16 17 18   m m m m m m m 22 23 24 25 26 27 28     m m m m m m 33 34 35 36 37 38     m m m m m 44 45 46 47 48   m = (9.55) e   m m m m 55 56 57 58     m m m 66 67 68     sy. m m 77 78 m 88 “chap09” — 2002/12/14 — page 213 — 15214 CHAPTER 9 FEM FOR 3D SOLIDS where    1 1 1 m = ρabcN N dξ dη dζ ij i j −1 −1 −1       N 00 N 00 1 1 1 i j    =ρabc 0 N 0 0 N 0 dξ dη dζ i j −1 −1 −1 00 N 00 N i j      N N 00 1 1 1 i j   =ρabc 0 N N 0 dξ dη dζ (9.56) i j −1 −1 −1 00 N N i j or   m 00 ij   m = 0 m 0 (9.57) ij ij 00 m ij in which   +1 +1 m =ρabc N N dξ dη dζ ij i j −1 −1   +1 +1 ρabc = (1+ξξ)(1+ξ ξ) dξ (1+η η)(1+η η) dη i j i j 64 −1 −1  +1 × (1+ζζ)(1+ζ ζ) dζ i j −1     ρhab 1 1 1 = 1+ ξ ξ 1+ η η 1+ ζ ζ (9.58) i j i j i j 3 3 3 8 As an example,m is calculated as follows: 33     ρabc ρabc 1 1 1 m = 1+ × 1× 1 1+ × 1× 1 1+ × 1× 1 = 8× (9.59) 33 3 3 3 8 216 The other components of the mass matrix for a rectangular hexahedron element are: 8ρabc m =m =m =m =m =m =m =m = 11 22 33 44 55 66 77 88 216 m =m =m =m =m =m =m =m =m =m =m 12 23 34 56 67 78 14 58 15 26 37 4ρabc =m = 48 216 (9.60) m =m =m =m =m =m =m =m =m =m =m 13 24 16 25 36 47 57 68 27 38 45 2ρabc =m = 18 216 1ρabc m =m =m =m = 17 28 35 46 216 “chap09” — 2002/12/14 — page 214 — 169.3 HEXAHEDRON ELEMENT 215 Note that the equalities in the above equation can be easily figured out by observing the relative geometric positions of the nodes in the cube element. For example, the relative geometric positions of nodes 1–2 is equivalent to the relative geometric positions of nodes 2–3, and the relative geometric positions of nodes 1–7 is equivalent to the relative geometric positions of nodes 2–8. If we write the portion of the mass matrix corresponding to only one translational direction, say thex direction, we have   842 4 4 2 1 2   84 2 2 4 2 1     841242     ρabc 82124   m = (9.61) e   8424 216     842     sy. 84 8 The mass matrices corresponding to only they andz directions are exactly the same as m . e The nodal force vector for a rectangular hexahedron element can be obtained using Eqs. (3.78), (3.79) and (3.81). Suppose the element is loaded by a distributed force f on s edge 3–4 of the element, as shown in Figure 9.10; the nodal force vector becomes    f    sx  T f = N f dl (9.62)  e sy 3–4  l f sz If the load is uniformly distributed,f ,f andf are constants, and the above equation sx sx sz becomes   0  3×1       0  3×1       f  sx        f  sy       f   sz     1 f   sx f = l (9.63) e 3–4 f   2 sy       f   sz       0   3×1       0   3×1       0   3×1     0 3×1 where l is the length of edge 3–4. Equation (9.63) implies that the distributed forces 3–4 are equally divided and applied at the two nodes. This conclusion suggests also to evenly distribute surface forces applied on any face of the element, and to evenly distribute body forces applied on the entire body of the element. “chap09” — 2002/12/14 — page 215 — 17216 CHAPTER 9 FEM FOR 3D SOLIDS 8 5 4 8 1 1 5 3 8 8 6 4 1 6 7 8 7 6 2 3 3 6 6 1 1 3 3 2 Figure 9.11. A hexahedron broken up into five tetrahedrons. 9.3.3 Using Tetrahedrons to form Hexahedrons An alternative method of formulating hexahedron elements is to make use of tetrahedron elements. This is built upon the fact that a hexahedron can be said to be made up of numerous tetrahedrons. Figure 9.11 shows how a hexahedron can be made up of five tetrahedrons. Of course, this is not the only way that a hexahedron can be made up of five tetrahedrons, and it can also be made up of six tetrahedrons, as shown in Figure 9.12. Similarly, there is more than one way of dividing a hexahedron into six tetrahedrons. In this way, the element matrices for a hexahedron can be formed by assembling all the matrices for the tetrahedron elements, each of which is developed in Section 9.2.2. The assembly is done in a similar way to the assembly between elements. 9.4 HIGHER ORDER ELEMENTS 9.4.1 Tetrahedron Elements Two higher order tetrahedron elements with 10 and 20 nodes are shown in Figures 9.13(a) and (b), respectively. The 10-node tetrahedron element is a quadratic element. Compared with the linear tetrahedron element (four-nodal) developed earlier, six additional nodes are added at the middle of the edges of the element. In developing the 10-nodal tetrahedron element, a complete polynomial up to second order can be used. The shape functions for “chap09” — 2002/12/14 — page 216 — 189.4 HIGHER ORDER ELEMENTS 217 8 5 4 1 6 7 8 2 3 8 5 4 6 7 6 1 4 2 3 2 Break into three 5 8 5 1 4 6 1 4 4 6 6 2 Figure 9.12. A hexahedron broken up into six tetrahedrons. (a) (b) 4 4 14 12 9 13 20 8 16 11 2 2 8 7 7 17 10 9 19 1 1 5 15 5 18 10 5 6 6 3 3 Figure 9.13. Higher order 3D tetrahedron elements. (a) 10-node tetrahedron element; (b) 20-node tetrahedron element. this quadratic tetrahedron element in the volume coordinates are given as follows: N =(2L − 1)L for corner nodesi = 1, 2, 3, 4 i i i  N = 4L L  5 2 3    N = 4L L  6 1 3   N = 4L L (9.64) 7 1 2 for mid-edge nodes N = 4L L  8 1 4    N = 4L L  9 2 4   N = 4L L 10 3 4 “chap09” — 2002/12/14 — page 217 — 19218 CHAPTER 9 FEM FOR 3D SOLIDS whereL is the volume coordinate, which is the same as the shape function for the linear i tetrahedron elements given by Eq. (9.15). The 20-node tetrahedron element is a cubic element. Compared with the linear tetra- hedron element (four-nodal) developed earlier, two additional nodes are added evenly on each edge of the element, and four-node central-face nodes are added at the geometry centre of each triangular surface of the element. In developing the 20-nodal tetrahedron element, a complete polynomial up to third order can be used. The shape functions for this cubic tetrahedron element in the volume coordinates are given as follows: 1 N = (3L − 1)(3L − 2)L for corner nodesi = 1, 2, 3, 4 i i i i 2  9 9 N = (3L − 1)L L N = (3L − 1)L L 5 1 1 3 11 1 1 4  2 2    9 9  N = (3L − 1)L L N = (3L − 1)L L  6 3 1 3 12 4 1 4 2 2     9 9  N = (3L − 1)L L N = (3L − 1)L L 7 1 1 2 13 2 2 4 2 2 for edge nodes 9 9 N = (3L − 1)L L N = (3L − 1)L L  8 2 1 2 14 4 2 4  2 2    9 9  (9.65) N = (3L − 1)L L N = (3L − 1)L L  9 2 2 3 15 3 3 4 2 2    9 9  N = (3L − 1)L L N = (3L − 1)L L 10 3 2 3 16 4 3 4 2 2  N = 27L L L 17 2 3 4    N = 27L L L 18 1 2 3 for centre surface nodes  N = 27L L L 19 1 3 4   N = 27L L L 20 1 2 4 whereL is the volume coordinate, which is the same as the shape function for the linear i tetrahedron elements given by Eq. (9.15). 9.4.2 Brick Elements Lagrange type elements The Lagrange type brick elements can be developed in precisely the same manner as the 2D rectangular elements described in Chapter 7. Consider a brick element with n = d (n+ 1)(m+ 1)(p+ 1) nodes shown in Figure 9.14. The element is defined in the domain of (−1≤ξ ≥ 1,−1≤η ≥ 1,−1≤ζ ≥ 1) in the natural coordinatesξ,η andζ . Due to the regularity of the nodal distribution along theξ,η andζ directions, the shape function of the element can be simply obtained by multiplying one-dimensional shape functions with respect to theξ,η andζ directions using the Lagrange interpolants defined in Eq. (4.82) Zienkiewicz et al., 2000: p 1D 1D 1D n m N =N N N =l (ξ)l (η)l (ς) (9.66) i I J K I J K Due to the delta function property of the 1D shape functions given in Eq. (4.83), it is easy to confirm that theN given by Eq. (9.66) also has the delta function property. i “chap09” — 2002/12/14 — page 218 — 20