v10 3. vi. 5' SN“. - I dn‘.‘ '1']. u n 03114.7” ('0. P]! (All: I . . I9 0 .u’ I! I, ' l . QC; . ‘l 1.13 1ng . w .140.“ Quin. unvuc ‘ .. . an 3. [vbbA’ (Ill. w V1.5. . I- 1‘ I v . ‘ LOWE. .taI..CID..vIlun..-vs I“ (vs 2111!} v 1."? ‘ . I LLLBEDDL lily). .. I : . .I: . I? . '37.. p. 1.; p I .»» »h ... .D but-rt. o ‘ 1‘ ‘ D it... ,A.Pn? . -k n. . f 5 III-- 7 .Ir H} O IH‘UITHTWI'Hifliifilfliimifilifiliflfiliflifiifl (We) 3 1293 01710 3817 This is to certify that the thesis entitled A STUDY OF PERFORMANCE OF FINITE ELEMENTS IN SHELL EIGENVALUE PROBLEMS presented by Akmal Shah has been accepted towards fulfillment of the requirements for MS degree in Mechanical Engineering WW Major professor Date HAY 3% 1C) 1V0 0—7639 MS U is an Affirmative Action/Equal Opportunity Institution LIBRARY Michigan State Universlty PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE 1/” WM“ A STUDY OF PERFORMANCE OF FINITE ELEMENTS IN SHELL EIGENVALUE PROBLEMS By Akmal Shah A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Mechanical Engineering 1996 ABSTRACT A STUDY OF PERFORMANCE OF FINITE ELEMENTS IN SHELL EIGENVALUE PROBLEMS By Akmal Shah A study is carried out to evaluate the performance of the shell finite element of Belytschko and Leviatahan (QPH) in eigenvalue computations. The QPH element is a four noded quadrilateral shell element with one point quadrature that is physically stabilized to take care of spurious modes that result from reduced integration. Both consistent and lumped mass matrices are used in the evaluation. The subspace iteration method is used for determining eigenvalues and eigenvectors. The problems of a cylinder, a hemispherical shell and a flat plate are analyzed in the study. All the problems are also solved using a commercial finite element package (IDEAS). Wherever possible, the results are also compared with analytical results. The problems are solved for a wide range of thickness values and various boundary conditions. Analysis using the QPH element compares well with other solutions for all the problems. To my wife and children. The joy of going back to them really inspired me to finish my work on time. ACKNOWLEDGMENTS I would first thank the God Almighty for providing me with the opportunity to undertake and complete my graduate studies. I would also thank Pakistan Air Force and National University of Sciences and Technology, Pakistan for enabling me to come to Michigan State University for pursuing my this study. I would like to express my gratitude and sincere thanks to my advisor Dr. Alejandro Diaz for his able guidance and support without which it may not have been possible to pursue and complete this study. Whenever I needed assistance, he would take out time from his very busy schedule and guide me from his vast knowledge and experience. I also want to thank Dr. R. C. Averill who taught me the finite element course in the most befitting manner that really enabled me to have a far better understanding with my advisor. I also thank him for his off and on help in my thesis. I can not forget to mention the assistance provided by Venkat Rao, a Ph.D. student in Department of Materials Science and Mechanics at MSU. Despite his own commitments of classes and research, he would spare his precious time to discuss the problems that I used to face. Eric Leaman, the finite element consultant at MSU, also deserves special thanks for his valuable guidance in learning IDEAS software. I am thankful to my wife for her sacrifices, patience and continuous prayers. This was a great source of inspiration and encouragement for me. I would also like to mention about my daughters Erum and Maryam, talking to whom on telephone would provide me with much needed pleasure and relief in the moments of grief and tension. I also thank my brothers and my in laws for their moral support and prayers. Their occasional letters from back home were really encouraging and motivating. Finally, I would like to thank Fida Muhammad Khan for his unforgettable moral and logistic support that helped me a great deal in accomplishing this task. TABLE OF CONTENTS LIST OF TABLES .......................................................................................................... vii LIST OF FIGURES ....................................................................................................... viii CHAPTER 1 INTRODUCTION .............................................................................................................. 1 General ......................................................................................................................... 1 Brief Details of Study .................................................................................................. 2 CHAPTER 2 THE FINITE ELEMENT MODEL ................................................................................. 5 Governing Equations .................................................................................................. 5 Geometry and Kinematics ........................................................................................... 6 Locking Spurious Modes and Stabilization ................................................................. 8 Formulation of Stiffness Matrix ................................................................................. 15 Consistent and Lumped Mass Matrices ..................................................................... 20 Solution of the Eigenvalue problem .......................................................................... 22 CHAPTER 3 RESULTS ......................................................................................................................... 24 Linear Static Analysis ................................................................................................ 24 Dynamic Analysis of Flat Plate ................................................................................. 24 Dynamic Analysis of a Cylinder ................................................................................ 32 Dynamic Analysis of a Hemisphere ........................................................................... 36 CHAPTER 4 DISCUSSION ON RESULTS AND CONCLUSIONS ............................................... 48 General ....................................................................................................................... 48 Effect of Thickness .................................................................................................... 49 Effect of Drill Degrees of Freedom Stiffness ............................................................ 53 Effect of Mass Lumping ............................................................................................ 55 Skipping of Modes ..................................................................................................... 58 Conclusion ................................................................................................................. 58 BIBLIOGRAPHY ............................................................................................................. 61 vi Table Table Table Table Table Table Table Table Table LIST OF TABLES 1 - Problem Parameters for Flat Plate ................................................................... 25 2 - Results of Natural Frequencies in Hertz for Flat Circular Plate with Clamped Edge (a / R = 0.01) ........................................................................... 29 3 - Results of Natural Frequencies in Hertz for Flat Circular Plate with Clamped Edge (a / R = 0.50) ........................................................................... 30 4 - Results of Natural Frequencies in Hertz for Square Plate with two Opposite Edges Clamped ................................................................................. 31 5 - Results of Natural Frequencies in Hertz for Square Plate with one edge clamped and an adjacent edge simply supported .................................... 31 6 - Problem Parameters for Cylinder ................................................................... 34 7 - Results of Natural Frequencies in Hertz for a Cylinder .................................. 35 8 - Effect of Drill Stiffness Parameter or on Natural Frequencies of Cylinder (a/R = 1/30) ....................................................................................... 35 9 - Problem Parameters for Hemisphere ............................................................... 37 Table 10 - Results of Natural Frequencies in Hertz for Hemisphere with Clamped Edge for various a/R Ratios .............................................................. 43 Table 11 - Results of Natural Frequencies in Hertz for Hemisphere with Quarter Edge Clamped for various a/R Ratios ................................................. 44 Table 12 - Results of Natural Frequencies in Hertz for Hemisphere with Four Points Simply Supported for various a/R Ratios ..................................... 45 Table 13 - Comparison of Natural Frequencies with Consistent(C) and Lumped(L) Mass Matrices for Clamped hemisphere with QPH element ........................... 46 vii Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure LIST OF FIGURES 1 - Midsurface of QPH element ........................................................................... 7 2 - Mesh for which Incompressibility dictates Zero Displacements .................... 9 3 - Standard Mesh ............................................................................................... 10 4 - Circular Plate with Clamped Edge ................................................................ 27 5 - Square Plate with Opposite Edges Clamped ................................................. 28 6 - Square Plate with One Edge Clamped and an Adjacent Edge Simply Supported ...................................................................................................... 28 7 - Cylinder with Rigid Diaphragm at Each End ............................................... 33 8 - Hemisphere with Circumferential Edge Clamped ........................................ 38 9 - Hemisphere with Quarter Circumferential Edge Clamped ........................... 39 10 - Hemisphere with Four Points on Circumferential Edge Simply Supported .......................................................................................... 40 11 - Mode Shape of Second Mode of Thin Circular Plate with Clamped Edge ...51 12 - Mode Shape of Third Mode of Thin Circular Plate with Clamped Edge ..... 52 13 - Effect of Drill Stiffness Parameter on 6th Natural Frequency of a Cylinder ................................................................................................. 54 14 - Consistent and Lumped Mass Matrices Result for Clamped Hemisphere with QPH Element (a/R = 0.001) ............................................. 56 15 - Consistent and Lumped Mass Matrices Result for Clamped Hemisphere with QPH Element (a/R = 0.100) ............................................. 57 viii Chapter 1 INTRODUCTION General Considerable research has been devoted in the recent past for developing shell finite elements that are not subjected to the effects of locking and kinematic modes. Reduced or selective integration techniques have been widely used but even these have not always been successful in eliminating locking. Moreover, an excessively reduced integration may trigger spurious kinematic modes and thus lead to a kinematically unstable finite element model. The kinematic modes resulting from the use of reduced integration can be controlled by using a stabilization matrix. Stabilization can be achieved by adding artificial generalized strains that are orthogonal to all linear fields. The magnitude of associated generalized stresses in this case are governed by user input hourglass control parameters. These parameters are usually chosen to be just large enough to prevent hourglassing and the results may be sensitive to these parameters. This approach is therefore often termed as perturbation hourglass control. Many approaches have been used to nullify the difficulty introduced by the ad hoc control parameters in perturbation stabilization. Some of the approaches are (1) Developing a stabilization procedure based on three field Hu-Washizu principle as demonstrated by Belytschko and Bachrach [1]. 2 (2) Modifying the approach given in (1) above for nonlinear problems by augmenting the rank of the under-integrated element by assuming an strain field that does not require any stabilization parameters as done by Belytschko and Bindeman [2]. (3) Using an assumed stress field for the membrane and bending fields and an assumed strain field for transverse shear field as shown by Englemann and Whirley [3]. Brief Details of Study Belytschko and Leviathan [4] have introduced a four noded quadrilateral shell element using one point quadrature with physical stabilization and have named it QPH. The shell formulation used for this element is a Mindlin-Reissner formulation which admits transverse shear but assumes that pseudo-normals remain straight. In classical thin plate theory, normals to the mid-surface before deformation are assumed to remain normal to the mid-surface after the deformation, implying that transverse shear effects are negligible. As a result the free vibration frequencies calculated by using the thin plate theory are higher than those obtained by the Mindlin plate theory in which transverse shear and rotary inertia effects are included. In this study, the performance of the QPH element in eigenvalue problems is evaluated. The results for linear static analysis using the QPH element are shown in [4] and are quite satisfactory. However, no results are shown in [4] for the dynamic analysis and this study is thus carried out to evaluate the performance of QPH element in eigenvalue problems. 3 A computer program is developed for the purposes of this study that uses the closed form expression of the element stiffness matrix of [4] and the consistent and lumped matrices for shell elements given by Hughes [5]. The mass matrices given in [5] are in integral form and are numerically integrated in the program. The shape function matrix given by Hughes [5] for consistent mass is simplified by assuming that the fiber coordinate system coincides with the local Cartesian coordinate system for the QPH element. Both consistent and lumped mass matrices are further simplified by assuming a constant thickness over the element domain. The subspace iteration method is used for determining eigenvalues and the subroutines given by Bathe [6] are used in the program. Skyline storage is used for stiffness and mass matrices as required by subspace iteration subroutines of [6]. This form of storage also minimizes the storage requirements. The problems of a flat plate, a cylinder and a hemisphere are selected for this study as these problems are a part of the set of standard test problems to test the performance of a finite element given by MacNeal and Harder [7]. All problems are also solved using a commercial finite element package, IDEAS [8] . The element used from IDEAS element library is a thin shell linear quadrilateral element that is formulated on the basis of Mindlin plate theory. The results using the QPH element are compared with the results from IDEAS. For a flat plate problem, the results are also compared with analytical results that are given by Blevins [9]. No method has been devised for plotting the mode shapes and the same are plotted through IDEAS as the natural frequencies obtained from QPH are almost same as those obtained from IDEAS. The QPH element performs well when compared with IDEAS solutions and with the analytical solution of [9] in the case 4 of a flat plate. The QPH element is, however, faster than the fully integrated elements as it retains the speed of one point quadrature family of elements due to closed form integration of the stabilizing forces. The theoretical background of the element including governing equations, geometry and kinematics, formulation of stiffness matrix, formulation of consistent and lumped mass matrices and the solution technique, is given in Chapter 2. The results for the three problems under different boundary conditions and over a range of thickness values are given in Chapter 3. This chapter also includes the plots of the mode shapes that are obtained through IDEAS for some specific examples. Results for both consistent and lumped mass matrices are also given there. In Chapter 4, the results are discussed and conclusions are drawn to conclude the study. Chapter 2 THE FINITE ELEMENT MODEL Governing Equations The equations used for the QPH shell element are developed by starting with the equations for a continuum and then applying suitable constraints. The governing equations are then the momentum equation 0' + bi = pvi (2.1) M where 0'- is the Cauchy stress, bi is the body force, p is the density and Vi is the I] acceleration. The deformation of the continuum is measured by the rate of deformation tensor n given by 1 nij = 5(Vij + Vj,i) (22) The hypoelastic constitutive equation relating stress with T] is 8 = CT] (2-3) A denoting a local coordinate system. 6 The formulation of the QPH element is a Mindlin-Reissner formulation that admits transverse shear but assumes that pseudo-nonnals remain straight. The formulation thus begins with a full three dimensional velocity field 8 V. (Eng) = ZN?D(§,11,C M (2.4) I=l where NI3D are the full three dimensional isoparametric shape functions. Since the pseudo-normals remain straight so 4 a 2_ ,_ v. (£311.) = ZN.(§,n)[v.. + EU-Vuw .. + V..w ,. )1 (2.5) I=l where a is the thickness. Functions N[(§,n) are the two-dimensional isoparametric shape fimctions and V,1 and V,2 are the two nodal base vectors which are orthogonal to the pseudo-normal at a node 1. Geometry and Kinematics The geometry of the 4-node shell element is defined by its midsurface, as shown in Figure 1, with coordinate denoted by X” ( i = 1, 2, 3 representing the three Cartesian components and I = 1, 2, 3, 4 representing the four nodes), and by its thickness 3. The local Cartesian coordinate system for the element is defined by unit vectors é, , A 6,, 63 where 6, and 62 are tangent to the midplane at the element center and 63 is orthogonal to the midsurface. A 0y! ’03 yl 4 C . / A ,. / A —; vxl ’0) x1 my 3 ——B D A §,x 2 l A Figure 1 - Midsurface of QPH element Thus A x e, = M (2.6) Is. >< g2] where g, and g2 are covariant base vectors such that 6x 8x g1 : —— ’ g2 = —— (2.7) ai an where x = ( x, , x2 , x3 ) are the spatial coordinates and g and n are the referential coordinates. In Cartesian component form (expressed in the global system), g1 takes the form _a_x,_ 1+33£2_ez+§x_333 (2.8) “6:8 an 6; 8 The base vectors 6, , e, are tangent to the midsurface at the center and are orthogonal to E3 and to each other. Thus only one can be arbitrarily selected. The selection for QPH element is e] = i (2.9) [all and e, = e, x e, (2.10) Five degrees of freedom are defined at each node for the QPH element, three translational velocities 0”, {721, {731, and two out of plane rotational velocities (3 x, , (f) y, . All the degrees of freedom are expressed in the local coordinate system defined at the element center. The finite element interpolation for the midsurface geometry is £1=N1(§,n)f,, (2.11) The velocity interpolation for QPH element in the corotational components is same as equation (2.5) and can be written in terms of indicial notation as 91(§m,C)= N,(§.n)1 2 too few incompressibility constraints r = 2 optimal r < 2 too many incompressibility constraints r S 1 locking In the development of plate bending elements, an important consideration is the number of shear strain constants engendered in the thin plate limit. The constraint ratio 12 r, defined in equation (2.13) is again used here, but now r1,q is the total number of displacement and rotation equations after boundary conditions have been imposed and nc is the total number of shear strain constraints. Again, the idea is that as the number of equations in the standard mesh of Figure 3 approaches infinity, r should approximate the ratio of number of equilibrium equations to number of constraints for the given system of partial differential equations (in the present case, 3 and 2, respectively). The ideal ratio of r here would thus be 3/2. Smaller value will indicate too many shear strain constraints and a potential for locking. A larger value would indicate too few shear strain constraints and may poorly approximate the Kirchoff limit. The effect of locking can be reduced by using uniform reduced integration or selective reduced integration. In uniform reduced integration for plate elements, both the bending and shear terms are integrated with the same rule, which is of lower order than the one required for exact integration. In selective reduced integration for plate, the bending term is integrated with the normal rule (for exact integration), whereas the shear term is integrated with a lower order rule. The reduced integration represents a considerable improvement over normal integration. For a four noded isoparametric rectangular element with sides aligned with global x- and y- axes, the element expansion may be written as Wh=Bo+B1x+Bzy+B3xy (2-14) 93 =rao+rmx+ra2y+m3xy (2.15) 13 where Bi and ya, , OS i S 3, are constants that depend upon the nodal parameters w"; andOZa , Is a S 4 representing the four nodes, respectively. The conditions 0 = 1’1 = -91’+wf; = (-vlo+Bi)-1nx+(-712+B3)y-713xy (2-16) 0 = 72 = -O'2'+w"; = (-y20+l32)+(-y21+03)x-y22y-723xy (2-17) impose eight constraints per element and are approximately in force as t —> 0 if exact integration of element stiffness matrix is performed, which is two-by-two Gauss integration in this case. In a large rectangular mesh, there are approximately three degrees-of-freedom per element, and thus the element tends to be overly constrained with r = 3/8, which is indicative of locking. To reduce the effect of locking, using one-point Gauss quadrature for element stiffness matrix may be considered. Clearly this results in only two constraints per element, and now there are more degrees of freedom than there are constraints. The value of r is now 3/2, which is the optimal value. Although reduced integration alleviates the locking problem, it has an adverse side effect of rank deficiency, that is there are zero energy modes present in excess of the rigid body modes. These zero energy modes are known as spurious modes and can result in oscillatory errors in certain cases and, occasionally, even in a singular global stiffness matrix. Spurious modes are the eigenvectors of the stiffness matrix, such that 14 [K] {(1),} = 0 (2.18) A spurious mode, {(1),}, differs from a rigid body mode, {(1)0}, in that the modal strain state {Es} = [B] {(135} (2-19) may not be zero everywhere. It is, however, zero at the integration points. Spurious modes occur if the number of strain evaluations at integration points is less than the number of independent strain states provided by the nodal displacement. Thus the linear combinations of strain states will produce zero strains at all integration points. The difference between the number of independent strain states and the number of strain evaluations provide a lower limit to the number spurious modes The selective reduced integration has the advantage of fewer spurious modes than for the corresponding uniform reduced integration. In uniform reduced integration all the strain components are evaluated at a reduced set of points. However, locking is a result of incorrect interpolation of some, but not all, components of strain. In selective reduced integration some, but not all, of the strain components are evaluated at a reduced set of points. This should eliminate locking and may retain a sufficient number of strain evaluations to prevent spurious modes. This method still has some undesirable side effects. In general, the boundary conditions render the assembled stiffness matrix positive-definite, and so the zero energy modes are not globally present. A guideline that has been successfully used for selective integration elements is that if the boundary conditions preclude the rigid mode from forming in one element, it is also precluded in the remainder of the mesh. Due to the potential problems created by rank deficiency, a 15 number of research efforts have been undertaken to develop Mindlin-Reissner-type elements which have correct rank and maintain accuracy in both thin and thick plate applications. The QPH element is also based on the Mindlin-Reissner theory. In selective reduced integration some or all of the economic advantage of reduced integration is given away. In situations where this advantage may be of particular importance, it would be highly beneficial to retain uniform reduced integration while finding some simple way to restrain or stabilize the resulting spurious modes. The basic idea behind spurious modes stabilization is to separate the element’s modes into low order modes which are evaluated at integration points and higher order modes whose stiffnesses are approximated analytically. Care must be taken to preserve the accuracy of lower order modes while assigning penalty stiffness to the higher order modes to prevent large responses of spurious modes. Spurious modes stabilization occupies a middle position between reduced integration and full integration with strain modification. It retains some of the economic advantage of reduced integration and avoids spurious modes. Its main disadvantage is that it places a burden on the user to select a level of penalty stiffness. A level that is too low may cause the spurious modes to be evident in the solution, while a level that is too high may allow locking to return. An acceptable level may not always be possible. Formulation of Stiffness Matrix In deriving the stiffness matrix for the linear static case, the linear strains are related to the displacements exactly as the rate of deformation is related to the velocity. 16 The following orthogonalities have been used in deriving the element stiffness matrix for QPH element:- (a) Orthogonalities between constant in plane fields (membrane and bending) and the non-constant ones. (b) Orthogonalities between membrane and other fields. (c) Orthogonalities between shear and other fields. ((1) Orthogonalities between the decomposed non constant shear fields. The element stiffness matrix in local coordinates can then be written as A A K=K+KM, am» where K, is obtained by one point quardature rule with the quadrature point §=n=0, A K is the element stiffness matrix in the local coordinate system that consists of 16 blocks of (K)U , each of dimension 5x5, with I and J varying from 1 to 4 and Km is the stabilization matrix. The constitutive law for the stabilization matrix in (2.20), given in the local coordinate system is 6=éé a2n where 6 and 5 are the 5x1 stress and strain vectors, respectively, and C is a 5x5 matrix that can be partitioned into a 3x3 in plane part and a 2x2 transverse shear part as v 0 1 O (2.22) 1 0 - 1 — v 2 ( )_ O<~— Cm=Cb:1 2 -V b l7 1 0 G 0 (3:3): =2 (2.23) s 121+v o 1 6 0 G where E is Young’s modulus, v is Poisson’s ratio, G is the shear modulus and 5/6 is shear correction factor for a rectangular cross section. The element stiffness matrices must be expressed in global Cartesian coordinate system while assembling into global stiffness matrix. The 6x6 transformation matrix R maps from global to local coordinates as follows: a, = Rd, (2.24) where d, and cll are the nodal displacement vectors at node I expressed in local and global coordinates, respectively. The transformation matrix R has the form ~ A 0 R : [ 3x3 3x3 ] (2.25) 03x3 A3x3 where A3,.3 is computed using the unit vectors 6, , é, , é3 in the local Cartesian coordinate system. The QPH element has 5 degrees of freedom per node. No stiffness is associated with the drilling degrees of freedom. This may result in a singularity or ill conditioning of the stiffness matrix. A value of torsional stiffness factor kd for these degrees of freedom is calculated as below k, = a2 [ [xzaCs(1,l) + yzaCs(2,2)]dA (2.26) or kd = a2 J-[xzaCsU ,1) + yzaCs(2,2)]di (2.27) o 18 where or is a parameter taken from Quint [10], j is the determinant of the jacobian and [d0 indicates integration over a 2x2 master element. The jacobian maps an arbitrary :1 element onto a 2x2 master element that is a square element with dimension of each side equal to two in referential coordinate system. g and 1] thus vary from -1 to +1 over this element. In equations (2.26) and (2.27), aCs is the transverse shear stiffness, Cs being the 2x2 matrix in equation (2.23), and x and y are as below x= €181(§)+mir(n)+h.fir(§n) (2.28) y = €191(§) +n191(n)+ h19.(§n) (2.29) In equations (2.28) and (2.29), hl is the product of referential coordinates 9;! and n, at a node 1. Equation (2.27) is numerically integrated to evaluate the value of kd which depends upon a drill stiffness parameter or that is specified as input. For thick shells, where thickness to radius ratio is more than 1/20, this is seen to be affecting the natural frequencies for higher modes. With the addition of these degrees of freedom the K” block of the element stiffness matrix becomes a 6x6 block that facilitates transformation of element stiffness matrix from local to global coordinates. The element equilibrium equation for linear elastostatics in local coordinates can be written as A (1') A (e) A (e) [K] {11} = {F} (2.30) When the displacement and force vectors are transformed to global coordinates using transformation of equation (2.24), the above equation becomes [Kl‘e’eruw = 1111111“ (2.31) Since nodal displacements and forces of each node are to be multiplied by R , the matrix R in equation (2.31) would be a 24x24 matrix, for the four noded element, with blocks of 3x3 matrix A of equation (2.25) along the diagonal and zeros in the off diagonal positions. The equation (2.31) can also be written as [R] [K]“’[R]{U }‘”’= {F }“" (2.32) and since transformation matrices are orthogonal, so [R]_l = [R]t and equation (2.32) can be written as [R] [K]“"[R]{U }“”= {F }“’ (2.33) The element equation in global coordinates can be written as [K]“’ {u}“” = {F}“” (2.34) Comparing (2.26) with (2.27), it can easily be noticed that [K]“’ = [R]‘[K]”’[R] (2.35) The element stiffness matrix in the global coordinates can thus be written as KU = R‘KUR (2.36) where K” is a 6x6 block of element stiffness matrix after the addition of drilling degree of freedom and I and J vary from 1 to 4 for the four nodes. 20 Consistent and Lumped Mass Matrices The consistent and lumped mass matrices for shell elements as described by Hughes [5] are used in this study. The matrices have been used for the specific case of constant thickness. The fiber coordinate system is assumed to coincide with the local Cartesian coordinate system for the QPH element. The shape function matrix given by Hughes [5] for the consistent mass matrix has thus been modified as N, o o 0 —N,z N, = o N, 0 N,z 0 (2.37) o 0 N, o 0 The consistent mass matrix is therefore given by m=[m,,] (2.38) , where mu = [ [N;N,pjdzda (2.39) 0-2V 2 The product N:N , , after integrating over the thickness, becomes pN,N,a 0 0 0 0 O NINJa 0 0 0 1 O 0 NINJa 0 O N,NJ = 1 3 (2.40) 0 0 0 —N,N,a‘ 0 12 0 0 O 0 éNINfi3 The above expression for N}NJ is used and numerical integration is carried out to evaluate m” of equation (2.39), assuming a constant density, to obtain the consistent mass matrix. 21 I01 For lumped mass matrix, the procedure given by Hughes [5] is adopted. m, , except drilling degrees of freedom, that is the mass associated with each rotational degree of freedom of node I is initially calculated as mf‘” = [prij1 (2.41) CI .. +1 where j = [de (2.42) -1 For constant thickness a, j becomes aj. No mass is assigned to the drilling degrees of freedom. For the QPH element m?” = mf" I = 1,2, ..... ,nen (2.43) where m?” is the mass associated with each displacement degree of freedom and nen is the number of element nodes. Proportional lumping is then done to construct the lumped mass matrix. The volume V of the element is numerically integrated using the following integral v = [360 (2.44) D For constant density, the mass of the element (M) is evaluated as M = p V (2.45) The masses associated with rotational and displacement degrees of freedom are then normalized as my“ +- (M / 1W" )m;°‘ (2.46) m?” <— (M /1’\”4disp)m;“w (2.47) 22 where M” = 2...: mi“ (2.48) H and WW = Z“ in?” (2.49) 1=l Summations in equations (2.48) and (2.49) are carried out for one degree of freedom only. Finally the rotational inertia terms are adjusted. For the QPH element with constant thickness, the z-coordinate at the midsurface at each node is zero and both the thickness at each node of the element and the average thickness of the element are equal to the thickness a. The adjustment to rotational inertia terms given by Hughes [5] is thus simplified as m?“ (— 011m?” (no sum on I) (2.50) 2 where or, = max L,é (2.51) 12 8 and A is area of the element that is evaluated as A = V / a (2.52) The values of m1 and mf‘s" in equations (2.50) and (2.47), respectively, are used at the corresponding diagonals of the lumped mass matrix. Solution of the Eigenvalue Problem The subspace iteration method is used for determining eigenvalues and eigenvectors. The stiffness and mass matrices are stored in skyline form. A Sturm sequence check is used to check if any mode is skipped. Subroutines given by Bathe[6] are used for the 23 subspace iteration method. No method is devised for plotting the mode shapes and the same are plotted through IDEAS after solving all the problems using IDEAS. Chapter 3 RESULTS Linear Static Analysis The problems of a pinched cylinder with a rigid diaphragm at both ends, a hemisphere with radial loads at four points of the free circumferential edge and a circular plate clamped at the circumference with a transverse load at the center are solved to confirm the results of Belytschko and Leviathan [4] and to check the accuracy of the program. The circular plate problem is solved for a thin plate (a/R = 0.01) and a thick plate (a/R = 0.50). The QPH element gives satisfactory results for all the problems and in particular performs well for a thick plate for which Kirchoff plate theory predicts poor results. Dynamic Analysis of a Flat Plate A circular and a square plate are analyzed in this section. The problem parameters are shown in Table 1. The circular plate problem is solved for two thickness values, one in thin plate region and the other in thick plate region. A thin circular plate is considered to see the possible effects of membrane and shear locking and a thick circular plate is considered to see the performance of the QPH element in a region where rotary inertia 24 25 TABLE 1 PROBLEM PARAMETERS FOR FLAT PLATE Length L Radius R Thickness a Ratio a/R or a/L Young’s Modulus E Poisson’s Ratio v Density p Boundary Conditions Circular Plate 10 0.1 and 5.0 0.01 and 0.50 3.0E6 0.3 7.317372E-4 Circumferential edge of plate is clamped Rectangular Plate 10 0.05 0.005 2.9993 8E7 0.29 7.317372E-4 1- Opposite edges clamped 2- One edge clamped and an adjacent edge simply supported 26 and shear deformation effects become important. The circumferential edge of the circular plate is considered clamped for both the thickness ratios of the circular plate as shown in Figure 4. The square plate is solved for two different boundary conditions that are shown in Figure 5 and Figure 6. The problems are also solved using IDEAS software with a thin shell linear quadrilateral element which is based on Mindlin plate theory as stated in [8]. The analytical results for both circular and square plates could be obtained from Blevins [9]. The QPH results for these two plate problems are thus compared with IDEAS results as well as with analytical results. The results are shown in Table 2 and Table 3 for two thickness values of circular plate and in Table 4 and Table 5 for two boundary conditions for the square plate. In the thin circular plate problem, the results of QPH element are in good comparison with IDEAS results, but both results display an error with the analytical results, as can be seen in Table 2. The error is small in lower modes but in higher modes the error is as high as 12%. For the thick circular plate, the QPH results as well as the results of IDEAS do not compare well with the theoretical results obtained from Blevins [9], as can be seen in Table 3. This is discussed later in Chapter 4. The QPH results compare well with IDEAS results but QPH seems to be skipping some modes in this problem which is also discussed in Chapter 4. The square plate is analyzed for one thickness ratio only which is in very thin limit (a/L = 0.005), L being the length of a side of the square plate. The QPH results for this plate are in very good comparison with IDEAS and theoretical results as can be seen in Table 4 and Table 5. This was an expected result as the classical thin plate theory, on 27 \ Clamped Edge Figure 4 Circular plate with clamped edge 28 Clamped Edge Clamped Edge Figure 5 Square Plate with Opposite Edges Clamped Clamped Edge Simply Supported Edge Figure 6 Square Plate with one Edge Clamped and an Adjacent Edge Simply Supported 29 TABLE 2 RESULTS OF NATURAL FREQUENCIES IN HERTZ FOR FLAT CIRCULAR PLATE WITH CLAMPED EDGE (a/ R=0.01) 211.55 THEORY IDEAS QPH 31.52 31.79 31.69 65.56 68.21 67.86 68.21 68.23 107.56 112.78 112.88 114.98 114.48 122.64 133.50 133.22 157.40 169.14 168.61 169.14 169.45 187.56 211.55 210.72 211.57 30 TABLE 3 RESULTS OF NATURAL FREQUENCIES IN HERTZ FOR FLAT CIRCULAR PLATE WITH CLAMPED EDGES (a/ R=0.50) THEORY IDEAS QPH 1575.8 1038.5 1035.9 3278.1 1771.4 1762.7 1771.4 1770.1 5378.2 2102.7 2094.2 2102.7 6132.2 2452.1 2492.4 2498.8 2499.0 2509.8 7870.0 2787.9 2784.8 9378.0 3255.4 3239.4 31 TABLE 4 RESULTS OF NATURAL FREQUENCIES IN HERTZ FOR SQUARE PLATE WITH TWO OPPOSITE EDGES CLAMPED THEORY IDEAS QPH 108.23 108.65 108.65 128.93 129.30 129.38 212.18 213.28 213.28 298.73 303.80 303.80 328.28 333.03 333.03 388.29 392.80 392.79 TABLE 5 RESULTS OF NATURAL FREQUENCIES IN HERTZ FOR SQUARE PLATE WITH ONE EDGE CLAMPED AND AN ADJACENT EDGE SIMPLY SUPPORTED 311.27 310.22 THEORY IDEAS QPH 26.07 26.07 26.13 93.16 93.46 92.53 120.38 119.77 121.10 209.89 208.12 211.75 257.57 261.56 255.35 318.39 32 which the analytical results of Blevins [9] are based, predicts good results in this thin region. Dynamic Analysis of a Cylinder A right circular cylinder having a rigid diaphragm at each end, as shown in Figure 7, is analyzed here. The problem parameters are shown in Table 6. The problem is also solved using IDEAS software with a four noded quadrilateral element. In the absence of any analytical results, the QPH results are compared only with IDEAS results and are shown in Table 7. The QPH results compare very well with the IDEAS results. Analyses are carried out for this problem for thickness to radius ratios of l/ 1200 to 1/30 and QPH results compare well with the IDEAS results except for high thickness to radius ratios. The QPH results in this region could be improved by adjusting the value of drill stiffness parameter or (see equation (2.19)) and results are then in good comparison with IDEAS results. The effect of mesh refinement is studied in this problem and a 16x16 and a 20x20 mesh are considered. The natural frequencies are seen to be decreasing as the mesh is refined from 16x16 to 20x20, irrespective of the thickness of the cylinder. The effect of thickness is also studied in this problem and the thickness of the cylinder is varied from a very thin cylinder to a rather thick one (a/R from 1/ 1200 to 1/30) for both 16 x 16 and 20 x 20 meshes of the cylinder. Throughout, the QPH results are comparing very well with the IDEAS results. 33 \ X \ \\ \ \ \ \ A / A R l 1 Rigid Diaphragm Rigid Diaphragm Figure 7 Cylinder with rigid diaphragm at each end 34 TABLE 6 PROBLEM PARAMETERS FOR CYLINDER Length L Radius R Thickness a Ratio a/R Young’s Modulus E Poisson’s Ratio v Density p Boundary Conditions 600 300 0.25 - 10.0 1/1200 - 1/30 3.0E6 0.3 7.317372E-4 Supported at each end by a rigid diaphragm ux=uy=02=0 35 TABLE 7 RESULTS OF NATURAL FREQUENCIES IN HERTZ FOR A CYLINDER WITH A RIGID DIAPHRAGM AT EACH END 16 x 16 Mesh 20 x 20 Mesh IDEAS QPH IDEAS QPH 4.52 4.51 4.46 4.45 4.78 4.78 4.76 4.76 7.30 7.28 7.11 7.09 11.70 11.68 11.21 11.17 11.72 11.70 11.69 11.68 13.05 13.02 12.76 12.70 TABLE 8 EFFECT OF DRILL STIFFNESS PARAMETER or ON NATURAL FREQUENCIES OF CYLINDER IN HERTZ (a/R=1/30) IDEAS QPH RESULTS WITH BELOW VALUES OF or RESULT 1E-4 5E-4 8E-4 lE-3 lE-2 5E-2 lE-l SE-l 7.26 7.22 7.23 7.24 7.25 7.32 7.37 7.55 11.47 11.79 11.78 11.78 11.78 11.78 11.80 11.83 11.92 14.34 13.01 12.96 12.97 12.99 12.99 13.07 13.11 13.21 15.89 23.00 22.90 22.97 23.05 23.02 23.01 23.12 23.18 24.84 23.27 23.08 23.22 23.08 23.1 1 23.39 23.56 24.06 34.22 23.29 37.16 32.64 23.16 23.18 23.53 23.71 24.22 35.68 36 The effect of drill stiffness parameter or is also considered in this problem. At higher thickness values (a/R > 0.03), QPH does not calculate the higher modes frequencies correctly if the value of drill stiffness parameter or is not carefully adjusted. When this adjustment is done, QPH results compare well with IDEAS results. The adjustment to the value of or is done to have a good comparison of QPH results with IDEAS results. In the absence of any results to compare with, no method is available to adjust the value of or. A possible strategy is, however, suggested in Chapter 4. In a cylinder with thickness to radius ratio of 1/1200 to 1/30, results that compare well with IDEAS results are obtained with the value of drill stiffness parameter or equal to 1.0E-3. The results tend to drastically deviate for higher modes if the value of drill stiffness parameter or is decreased down to 1.0E-4 or below. If the value of or is increased beyond 1.0E-3, the results remain stable till 1.0E—l but are erroneous for the all the modes when the value is increased beyond 1.0E-1. The value of drill stiffness parameter or can thus be selected in the range of 1.0E-3-1.0E-1 for this problem. The results for various values of drill stiffness parameter or are shown in Table 8. Selection of the drill stiffness parameter or is further discussed in Chapter 4. Dynamic Analysis of a Hemisphere In this section the natural frequencies of a hemisphere under different boundary conditions are determined and analyzed. This problem was the most general problem addressed in this thesis and is thus studied in most detail. Problem parameters are shown in Table 9. The problem is analyzed for three different sets of boundary conditions that 37 TABLE 9 PROBLEM PARAMETERS FOR HEMISPHERE Radius R 10 Thickness a 0.01 - 1.00 Ratio a/R 0.001 - 0.100 Young’s Modulus E 6.825E7 Poisson’s Ratio v 0.3 Density p 7.317372E-4 Boundary l-Circumferential edge clamped Conditions 2-One quarter of circumferential edge clamped 3-Four points on circumferential edge at interval of 90 degrees simply supported 38 Clamped Edge Figure 8 Hemisphere with Circumferential Edge Clamped V 39 \ Clamped Portion Figure 9 Hemisphere with Quarter Circumferential Edge Clamped V 40 sym sym L Simple Supports Figure 10 Hemisphere with four points on circumferential edge simply supported 41 are shown in Figure 8, Figure 9 and Figure 10. In the problem shown in Figure 10, one quarter of the hemisphere is shown as symmetry was used in this problem. For the whole hemisphere in this problem, four points on the circumferential edge of the hemisphere at an interval of 90 degrees are simply supported. All rigid body modes are not restrained in this problem, as the displacements in the plane of circumferential edge are not restrained anywhere. Both QPH and IDEAS encountered singularities when the full hemisphere was used for this boundary conditions set without using the frequency shift option which was available in IDEAS. The problem with this boundary condition set could not be solved using QPH element without making use of symmetry, as the program developed in this study does not have the frequency shift option, and was thus solved using symmetry boundary conditions at one quarter of the hemisphere. The full hemisphere is considered in the other two boundary condition sets and no symmetry is used. The problem is also solved for all three boundary conditions sets with IDEAS software using a thin shell linear quadrilateral element. The thickness of the hemisphere was varied from a very thin limit to a thick limit (afR from 0.001 to 0.100) for all the boundary conditions. The behavior is exactly same as in the example of cylinder. The natural frequencies increase as the thickness is increased and the QPH results compare well with the IDEAS results. At higher thickness values the natural frequencies of higher modes given by QPH fall considerably as compared to those obtained from IDEAS. This situation is also controlled by increasing the value of drill stiffness parameter or. The increased value of or is to be carefully selected as the 42 results are sensitive to the value of this parameter. The selection of the value of or is further discussed in Chapter 4. The results for the three boundary condition cases are shown in Table 10, Table 11 and Table 12. The results are shown for six different thickness to radius ratios, starting from a very thin shell and going into the thick shell limit. In the thick shell region, the value of drill stiffness parameter or had to be increased to a very high value as compared to its value in the thin shell limits. For low thickness to radius ratios a value of 1.0E-3 for or is found to be appropriate but for high thickness to radius ratios the value of or has to be increased in the range of 0.5E0 to 1.0E0 to get results comparable with those of IDEAS in cases where full or part of the circumferential edge is clamped. The results shown in Table 10 are for or =1 .0E-3 for thickness to radius ratios upto 0.06 and are for or = 1.0EO for a/R equal to 0.08 and 0.10. Similarly the results shown in Table 11 are for a value of a = 1.0E-3 for thickness to radius ratios upto 0.80 and are for 01= 5.0E-1 for a/R= 0.1. The results shown in Table 12 are for or = 1.0E-3 upto a thickness to radius ratio of 0.06, for or = 7.0E-2 for a/R = 0.08 and for or = 8.0E-2 for a/R = 0.1. The effect of mass lumping is also studied for this problem. The clamped boundary condition example is selected to study the effect of mass lumping. Symmetry boundary conditions are used as the results for the same with consistent mass matrix were already obtained in the early part of the study. The results of both mass matrices are shown in Table 13. The QPH results obtained with lumped mass compare very well with the IDEAS results with the lumped mass option. The natural frequencies obtained using lumped mass matrix are found to be lower than the natural frequencies obtained using 43 TABLE 10 RESULTS OF NATURAL FREQUENCIES IN HERTZ FOR HEMISPHERE 5263.9 a / R = 0.001 IDEAS QPH 2729.3 2778.7 3632.7 3648.0 4306.2 4324.3 4307.4 4324.3 4393.8 4401.7 4401.9 4414.3 2728.8 2734.1 4515.5 4549.7 4629.4 4651.0 4637.7 4651.1 a I R = 0.060 IDEAS QPH 2901.6 2903.2 2902.5 2903.2 3950.4 3956.2 4542.2 4545.5 4556.2 4562.0 4726.2 4734.5 4726.5 4734.5 5045.8 5059.3 5050.3 5059.3 5287.9 a / R = 0.012 IDEAS QPH 2778.7 2780.5 3725.5 3733.4 4375.3 4386.8 4376.4 4386.8 4412.2 4415.8 4426.9 4432.4 2777.9 2780.5 4587.2 4617.1 4678.8 4690.6 4683.0 4690.7 a / R = 0.080 IDEAS QPH 2945.2 2949.9 2946.1 2949.9 4041.4 4053.0 4634.2 4646.3 4647.6 4663.3 4943.1 4963.7 4943.8 4963.7 5282.2 5314.9 5286.9 5314.9 5677.9 5715.7 WITH CLAMPED EDGE FOR VARIOUS a/ R RATIOS a / R = 0.020 IDEAS QPH 2803.8 2805.2 3769.8 3776.5 4418.1 4428.0 4418.9 4428.6 4424.8 4428.6 4440.0 4445.0 2802.9 2805.2 4651.4 4681.4 4719.1 4730.6 4723.1 4730.6 a / R = 0.100 IDEAS QPH 2987.4 2991.8 2988.4 2991.8 4137.4 4148.6 4743.9 4756.6 4756.6 4773.7 5184.7 5203.5 5185.6 5203.5 5551.3 5587.2 5556.2 5587.2 6037.2 61 16.8 RESULTS OF NATURAL FREQUENCIES IN HERTZ FOR HEMISPHERE 44 TABLE 1 1 WITH QUARTER EDGE CLAMPED FOR VARIOUS a/ R RATIOS a / R = 0.001 IDEAS QPH 23.2 23.5 25.3 25.4 75.1 75.5 82.0 82.5 101.8 162.3 162.3 163.6 264.8 266.6 279.0 278.0 397.6 399.4 419.6 421.7 a I R = 0.060 IDEAS QPH 148.0 148.8 203.1 203.0 476.4 478.7 647.0 644.4 1126.9 1123.6 1371.0 1371.5 1947.7 1949.5 2081.4 2076.8 2907.1 2896.1 2997.2 2995.8 a / R = 0.008 IDEAS QPH 39.3 40.3 42.9 43.7 131.4 132.9 138.6 139.2 273.3 275.0 296.2 297.7 478.1 480.1 482.4 486.0 687.0 692.0 741.8 744.7 a / R = 0.080 IDEAS QPH 184.6 185.5 238.0 238.5 551.8 553.7 810.7 808.0 1394.2 1390.3 1549.0 1553.3 2263.4 2258.4 2474.5 2470.6 3425.9 3397.1 3440.0 3407.9 a l R = 0.020 IDEAS QPH 73.0 73.7 86.5 86.5 257.5 258.1 286.2 287.1 527.8 530.9 600.6 600.8 896.6 898.4 1059.9 1061.3 1421.2 1422.9 1549.0 1555.9 a / R = 0.100 IDEAS QPH 220.3 226.2 263.6 269.7 629.0 635.7 952.4 960.2 1629.7 1642.3 1663.7 1677.3 2574.7 2595.1 2763.9 2786.1 3785.2 3806.3 381 1.1 3839.9 45 TABLE 12 RESULTS OF NATURAL FREQUENCIES IN HERTZ FOR HEMISPHERE WITH FOUR POINTS SIMPLY SUPPORTED FOR VARIOUS a/ R RATIOS a/R=0.00l a/R=0.004 a/R=0.040 IDEAS QPH IDEAS QPH IDEAS QPH 38.5 36.6 126.2 126.8 895.3 896.5 90.1 85.9 296.2 297.0 1888.1 1900.3 165.6 159.2 550.9 550.1 2988.6 3021.3 253.0 244.2 849.9 844.5 3548.7 3571.7 367.4 354.3 1244.9 1232.7 4369.7 4368.4 497.0 476.8 1696.4 1675.7 4590.2 4609.1 667.1 642.5 2313.8 2286.1 4704.0 4700.3 861.1 838.0 2937.4 2928.8 5124.5 5151.3 1135.2 1108.9 3622.8 3629.1 5172.7 5204.6 1518.7 1483.0 4108.8 4128.4 5204.7 5209.4 a/R=0.060 a/R=0.080 alR=0.l00 IDEAS QPH IDEAS QPH IDEAS QPH 1213.1 1215.4 1471.0 1492.0 1673.8 1700.2 2259.6 2284.6 2431.8 2478.4 2524.7 2573.9 3248.5 3280.7 3430.5 3475.5 3617.9 3659.5 4207.4 4172.1 4649.6 4666.7 4766.0 4784.4 4541.3 4548.3 4775.5 4796.0 4972.7 4991.0 4708.8 4722.0 5044.9 5094.6 5248.8 5294.1 5059.5 4948.5 5166.5 5210.1 5745.1 5803.8 5600.4 5004.7 6134.7 6143.0 6715.6 6720.1 5659.1 5172.9 6234.5 6263.0 6850.2 6867.0 TABLE 13 COMPARISON OF NATURAL FREQUENCIES WITH CONSISTENT(C) AND LUMPED(L) MASS MATRICES FOR CLAMPED HEMISPHERE WITH QPH ELEMENT UNDER SYMMETRY BOUNDARY CONDITIONS a / R = 0.001 C L 3648.2 2795.4 4401.8 2834.9 4549.9 2839.0 4734.8 2840.1 4764.1 2924.6 4790.7 2928.3 4882.3 2950.1 4891.1 2999.0 4917.4 3001.8 4961.2 3032.4 a / R = 0.016 C L 3755.8 3718.1 4421.5 4368.0 4646.9 4508.2 4800.1 4632.6 4867.7 4678.0 5004.7 4743.7 5097.0 4766.2 5140.2 4869.7 5241.8 4910.3 5482.0 4953.8 a / R = 0.004 C L 3678.0 3629.9 4406.0 4068.5 4567.9 4099.3 4736.9 4126.7 4778.9 4128.6 4809.9 4143.2 4884.2 4149.3 4926.6 4193.1 4932.0 4211.3 5012.3 4230.2 a / R = 0.040 C L 3868.5 3840.0 4474.9 4429.0 4930.5 4793.0 5112.6 4957.7 5180.6 5055.1 5936.7 5596.6 6061.1 5724.9 6101.6 5807.6 6181.5 5865.9 7400.0 6827.0 a / R = 0.008 C L 3708.2 3672.9 4410.8 4357.5 4591.3 4452.3 4749.1 4556.4 4802.1 4564.8 4849.7 4578.9 4927.6 4592.8 4977.8 4604.2 5027.5 4617.6 5095.1 4665.8 a / R = 0.060 C L 3957.8 3930.7 4548.0 4503.5 5293.3 5155.1 5525.6 5366.2 5574.8 5442.7 6876.8 6568.5 7144.0 6826.4 7242.3 6846.1 7273.6 6920.9 8110.9 7854.6 alR=0.012 C 3733.4 4415.8 4617.1 4771.2 4832.8 4915.0 4999.7 5048.0 5137.8 5246.4 a/R= C 4144.6 4751.4 6106.5 6531.4 6549.5 7845.8 9434.2 9466.0 9555.5 9862.2 L 3702.4 4367.8 4494.6 4615.9 4641.0 4704.7 4706.9 4744.2 4803.1 4858.8 0.100 L 4123.5 4714.0 6003.0 6386.2 6419.3 7723.5 9078.7 9097.0 9190.2 9474.6 47 consistent mass matrix. The error is, however, found to be not more than 5% except in thick region. The error for the first few frequencies is even smaller. The results with very low thickness to radius ratios (a/R < 0.01), however, display much greater errors in both QPH and IDEAS solutions which is discussed in Chapter 4. Chapter 4 DISCUSSION OF RESULTS AND CONCLUSION General The QPH element is based on Midlin-Reissner theory that takes into account the rotary inertia and shear deformation effects and thus performs better for thick shells where classical plate theory displays poor results. The QPH element displays good results for all the problems that are addressed in this study. It performs well in very thin shell regions where many elements face the effect of locking which verifies the stabilization scheme used for QPH element. The element is sensitive to the drill stiffness parameter or, especially for thick shells where the value of a must be carefully adjusted in order to get desirable results. The element works well with lumped mass matrix also but in very thin plates the performance deteriorates when lumped mass matrix is used. No study of computation times has been carried out but the QPH element should be much faster in comparison to the fully integrated elements due to the closed form integration of the stabilizing forces. The mode shapes plots have been made through IDEAS, as no method was devised in the program for plotting the modes. 48 49 Effect of Thickness For a thick circular plate the QPH results compare well with the IDEAS results but both QPH and IDEAS results are far off from the analytical results of Blevins [9]. This is because Blevins [9] has based his results on classical thin plate theory and has neglected the rotary inertia and shear deformation effects. In the problems where the rotary inertia and shear deformation terms become significant, the analytical results of Blevins [9] for natural frequencies are expected to be higher than those predicted by Mindlin plate theory. The thick circular plate considered for analysis has thickness to radius ratio of 0.5 for which the rotary inertia and shear deformation effects would be very significant. The natural frequencies obtained from Blevins [9] are thus higher than those obtained from QPH or IDEAS. The comparison of QPH results with IDEAS results for this problem indicates that some modes are being skipped when QPH element is used, as can be seen in Table 3. This is discussed later in this chapter. For thin plate, where the rotary inertia and shear deformation effects are far less important and can easily be ignored, the theoretical results obtained from Blevins [9] conform to the QPH or IDEAS results. In the thin circular plate problem, the natural frequencies of higher modes, obtained from both QPH and IDEAS have an error of about 12% in comparison to the analytical results of Blevins [9] as shown in Table 2. The error may be reduced by using a finer mesh. This error can not be attributed to the classical plate theory on which Blevins [9] has based his results, because if the rotary inertia and shear deformation effects are present the classical plate theory will predict higher frequencies whereas the results 50 shown in Table 2 indicate lower theoretical frequencies. Further as the mesh is refined the natural frequencies are known to converge to a lower value and thus refining the mesh further is expected to bring both QPH and IDEAS results closer to theoretical results. A very good comparison of QPH results with IDEAS results for the same mesh further proves this point. It can also be seen in Table 2 that both QPH and IDEAS have repeated modes which are not repeated in the analytical results. This is because the repeated modes in QPH and IDEAS have a certain number of half waves along one of the axes in the plane of the plate in one mode and the same number of half waves along the other axis in the second repeated mode. In the analytical results of Blevins [9] these modes are combined. The number of half waves along the radial direction is considered by Blevins [9] which could be along any of the axes in the plane of the plate. The plots for the second and third modes for this problem, which are repeated modes with both QPH and IDEAS, are plotted through IDEAS and are shown in Figure 11 and Figure 12, respectively, to illustrate this point. The results for two boundary conditions of a square plate are in very good comparison with IDEAS and theoretical results. Only a thin plate was analyzed and classical plate theory based results of Blevins [9] are, as expected, very close to the Mindlin-Reissner theory based results of QPH. For the cylinder and hemisphere problems no analytical solutions could be found and the study of the effect of thickness is restricted to the performance of QPH element in thin and thick shells in comparison to the IDEAS solutions. For a thin hemisphere problem with clamped edge, the performance of the QPH element declines in comparison to the 51 lhe-o/uo/uhnnnllA/old.utl Ilnvbflr l- I.c. I.Ionl 1.DIIILICIIIIT_I IODIX I Pllfi: CI.3II1| - IAG lrl. 0.00I‘co In]: 3.49Ioon DI'CIIITIOII 1' Le. 3,-0DI 2,3!‘ILLCIIII?_2 '03.! I PIIQI ‘I.IIGIE Dll'lalcllllf * no III: (LODIOOO III: Fl... 0' II'! III VII-III O'TXOIIACWIL 1.40loofl :.acnood 1.99Iood I .10'000 l.l)|oao I .ISIOIO I .D‘II'OI 7.!7I-01 lrlII-Ol .3 Figure 11 Mode Shape of Second Mode of Thin Circular Plate with Clamped Edge 52 /ho-o/uo/IllnaI-l/otd.Itl IIIULTI: J- I.C. 11.0 I J.DXI'LICIIIIT-J IODII J VIIG: 5.13141! ' IIG Ill: 0.00.000 Ill! I.IDIODO DI'OIIITXOI: I- I.C. IyIODI l.DII'LlCIII-?_J IOIII I IIIQI U..310)6 DIIFLACIIIIY ' Ilfi III: 0.00.900 IA]: . VILDI O’TXOIIICTD‘L VIII. 0' III: P If I.ADI§D° 3.10.000 I .DIIQOO 1.7lloou l .‘DIOIIO r.allooo 9.97l-OI 1.47l-01 I .III~01 Figure 12 Mode Shape of Third Mode of Thin Circular Plate with Clamped Edge 53 IDEAS when lumped mass matrix is used. This is further discussed later in this chapter. Effect of Drill Degree of Freedom Stiffness The effect of the drill stiffness parameter or is found to be very important for the QPH element. In linear static analysis the value of or is required to be kept low, somewhere in the range of 1.0E-3. The increase in the value of or results in considerably lower deformations but reduction in the value of or does not have any significant effect. In dynamic analysis also the value of or is to be kept in the same range of 1.0E-3 for thin shells. However, when the thickness is increased enough to move into the thick region, the higher modes frequencies are not found correctly if the same value of or is used. In the hemisphere problem the value of or was increased to a value in the range of 0.5E0- 1.0E0 to obtain results that compare well with the IDEAS results. In the cylinder problem, results that compared well with IDEAS results are obtained with the value of or equal to 1.0E-3. If the value is decreased to 5.0E-4 or below, the natural frequencies of the higher modes tend to increase sharply. However, when the value is increased beyond 1.0E-1, all the natural frequencies show a sharp upward trend. The natural frequencies thus appear to depend upon or in a quadratic manner. A graph is plotted in Figure 13 for the natural fiequency of 6th mode for the cylinder problem against the values of or to illustrate this phenomenon. A higher value of or is noted to be required for higher thickness ratios indicating that the proper choice of or directly depends on the thickness to radius ratio of the shell. The results shown in Table 10 verify this dependence where a Natural Frequency (Hz) 54 25« 1s - 1o 1. 1) 0.0001 0.001 0.01 0.1 Drtll Stiffness Parameter Figure 13 Effect of Drill Stifiness Parameter on 6th Natural Frequency of a Cylinder 55 gradual increase in the value of drill stiffness parameter a was required as the hemisphere was made thicker. Effect of Mass Lumping The use of lumped mass matrices is very common and is highly desirable due to the advantages of less storage space requirements and high computation speed. The study would thus be incomplete without studying the effect of mass lumping. The QPH results obtained with lumped mass compare very well with the IDEAS results with the lumped mass option. The real effect of mass lumping can only be studied if the same problem with same boundary conditions and using same element is solved with consistent and lumped mass matrices and a comparison of the two results is made. Hence the QPH results for the clamped hemisphere problem with consistent and lumped mass matrices are compared in order to study the effect of mass lumping. The results are shown in Table 13. The natural frequencies obtained using lumped mass matrix are found to be lower than the natural frequencies obtained using consistent mass matrix. The error is found to be not more than 5% except in thick region where error is higher. The error in results for the natural frequencies of first few modes is even smaller. For very thin shells (a/R < 0.01), both DEAS and QPH display much greater errors when lumped mass matrix are used which appears to be the result of mass lumping. The results of natural frequencies with consistent and lumped mass matrices for a very thin shell (a/R =0.001) and for a thick shell (a/R = 0.100) are shown in Figure 14 and Figure 15, respectively. Natural Frequency (Hz) 5000 4500 e 4000 + 3500 » 3000.. 2500«« 2000~ 1500., 1000~ 500«e Figure 14 Consistent and Lumped Mass Matrices Result for Clamped 56 { -o— Consistent i l: Lumped { 2 4 s 3 Mode Number Hemisphere with QPH Element (a/R = 0.001) 10 Natural Frequency (Hz) 57 5000 -. 4000’ l —o— Consistent | ~ + Lumped l 1 2 3 4 5 e 7 a 9 1o Mode Number Figure 15 Consistent and Lumped Mass Matrices Result for Clamped Hemisphere with QPH element (a/R = 0.100) 58 Skipping of Modes The Sturm sequence check is used in the program to determine if a mode has been skipped. For the thick plate problem it reported 10 eigenvalues missing. However, in comparison to the IDEAS results in Table 3, it appears that QPH element is skipping only two modes. This indicates that IDEAS is also skipping some modes for this problem. In the absence of enough evidence, it could not be ascertained if a particular mode i.e. membrane or bending was being skipped by QPH. The Sturm sequence check was also used for the hemisphere problem with clamped edge. The Sturm sequence check reported skipping of modes with the QPH element for thickness to radius ratios of 0.004, 0.040 and 0.080 indicating that the tendency of mode skipping exists for both thin and thick shells. The QPH and IDEAS results compare well for the clamped hemisphere but Sturm sequence check reported some modes being skipped with QPH element. Probably the same modes have been skipped by IDEAS also. The skipping of modes by IDEAS was observed in thick plate problem as well where out of the 10 modes that were reported to have been skipped by QPH, 8 were also skipped by IDEAS. Conclusion To conclude this study it can be said that QPH element performs very well for very thin shells where many elements would encounter locking. The stabilization scheme used in QPH element by Belytschko and Leviathan [4] is thus quite effective. However, a possible after effect of this stabilization is observed in the thick region where better 59 performance is generally expected with the Mindlin-Reissner plate theory. Since thick shells having thickness to radius ratio of 1/10 or more are not very common, QPH can be regarded as a very efficient element for eigenvalue computations. The results for dynamic analysis using the QPH element depend upon a drill stiffness parameter or. The results for thick shells greatly depend on the drill stiffness parameter or which requires careful adjustment. A general rule can be stated about the parameter or that the natural frequencies depend upon the value of or in a quadratic manner with good results at the apex. Thus if the natural frequencies are higher than expected then the value of or should be adjusted and if the frequencies further increase instead of reducing then the value of or should be changed in the other direction. If no other results are available with which the results of QPH can be compared, then the QPH results can be obtained for two more values of or in the vicinity of the value of or initially used. The new two values should be selected on either side of the initial value of or. The results for these three values may be used to determine the accuracy of results and any firrther adjustment required in the value of or. A value of 1.0E-3 appears to be the best choice for or in most situations and if the results with this value are not satisfactory then the value of a needs to be increased in almost all situations. A tendency of skipping of modes has been noticed with both QPH element and IDEAS. The tendency is observed in both thin and thick shells. In the absence of enough evidence it can not be said that the modes being skipped with QPH element are membrane or bending modes. 60 The use of mass lumping is highly recommended as it predicts frequencies that are very close to those predicted by consistent mass matrix. The small error that it gives is more than acceptable for the gain in storage requirements and solution time, that we get when lumped mass matrix is used. The error is somewhat greater in higher modes but even there it is not too significant. Further, in most practical situations we are interested in the first few modes only and there lumped mass matrices give good results. The results in very thin regions (a/R < 0.01) display much greater errors in both QPH and IDEAS solutions when lumped mass matrix are used which appears to be an effect of mass lumping. The use of lumped mass matrix for very thin shells should thus be exercised with caution. BIBLIOGRAPHY fl BIBLIOGRAPHY . T. Belytschko and WE. Bacharach, “Efficient Implementation of quadrilaterals with high coarse-mesh accuracy”, Computer Methods in Applied Mechanics and Engineering, 54 (1986) 279-301. T. Belytschko and LP. Bindeman, “Assumed strain stabilization of the 4-node quadrilateral with l-point quadrature for nonlinear problems”, Computer Methods in Applied Mechanics and Engineering, 88 (1991) 311-340. . B. E. Engelmann and R.G.Whirley, “A new elastoplastic shell element formulation for DYNA3D”, Lawrence Livermore National Laboratory, Report UCRL-JC 104826,1990. T. Belytschko and I. Leviathan, “Physical stabilization of the 4-node shell element with one point quadrature”, Computer Methods in Applied Mechanics and Engineering, 113 (1994) 321-350. T. G. R. Hughes, “The Finite Element Method linear static and dynamic finite element analysis”, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1987. K. J. Bathe, Finite Eelement Procedures in Engineering Analysis, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1982. R. H. MacNeal and R.L. Harder, “A proposed standard set of problems to test finite element accuracy”, Finite Elements in Analysis and Design, 11 (1985) 3-20. Structural Dynamics Research Corporation, “I-DEASTM Element Library P - 10022”, Structural Dynamics Research Corporation, Milford Ohio, 1993. R. D. Belvins, “F orrnulas for natural frequency and mode shapes”, Robert E. Krieger Publishing Company, Malabar Florida, 1984. 10. Quint Corporation, Tokyo, Japan, Private Communication. 61 nICHanN srnTE UNIV. LIBRQRIES lll1111111"lllllllWilllllllllllllllllllllllllllllllll 31293017103817