UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations   Help 
Material Information
Thesis/Dissertation Information
Subjects
Notes
Record Information

Full Text 
PAGE 1 1 MIXED FORMULATION US IN G IMPLICIT BOUNDARY FINITE ELEMENT METHOD By HAILONG CHEN A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2012 PAGE 2 2 2012 Hailong Chen PAGE 3 3 To my parents and wife PAGE 4 4 ACKNOWLEDGMENTS I would like to express my gratitude to my advisor and chairman of my supervisory committee, Prof Ashok V. Kumar, for his guidance, enthusiasm and constant support throughout my master research. I would like to thank him for the numerous insights he provided during every stage of the research. Without his assistance it would not have been possible to complete this thesis. I would like to thank the memb ers of my advisory committee, Prof. Loc Vu Quoc and Prof for providing help whenever required, for reviewing this thesis and valuable sug gestions provided. I would lik e to thank my undergraduate mentor, Na Li, for her numerous encouragement and support during my undergraduate study in China and graduate study in University of Florida, US. I would like to thank my wife and parents, for their constant love and support wit hout which this would not have been possible. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABL ES ................................ ................................ ................................ ............ 8 LIST OF FIGURES ................................ ................................ ................................ ........ 10 LIST OF ABBREVIATIONS ................................ ................................ ........................... 14 ABSTRACT ................................ ................................ ................................ ................... 15 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 17 1.1 Overview ................................ ................................ ................................ ........... 17 1.2 Goals and Object ives ................................ ................................ ........................ 19 1.3 Outlines ................................ ................................ ................................ ............. 19 2 MESHLESS AND MESH INDEPENDENT METHOD ................................ ............. 21 2.1 Traditional FEM ................................ ................................ ................................ 21 2.2 Meshless a nd Mesh Independent Method ................................ ........................ 22 2.3 Implicit Boundary Finite Element Method ................................ .......................... 23 3 MIXED FORMULATION FOR NEARLY INCOMPRESSIBLE MEDIA .................... 26 3.1 Overview ................................ ................................ ................................ ........... 26 3.2 Mixed Formulation ................................ ................................ ............................ 27 3.2.1 Matrix Decomposition ................................ ................................ .............. 27 3.2.2 Weak Form ................................ ................................ .............................. 29 3.3 2D Plane Strain ................................ ................................ ................................ 33 3.4 3D Stress ................................ ................................ ................................ .......... 35 3.5 Numerical Examples and Results ................................ ................................ ..... 38 3.5.1 Bracket (Plane strain) ................................ ................................ .............. 38 3.5.2 Beam (3D stress) ................................ ................................ ..................... 39 3.6 Concluding Remarks ................................ ................................ ......................... 39 4 CLASSICAL PLATE THEORIES ................................ ................................ ............ 43 4.1 Overview ................................ ................................ ................................ ........... 43 4.2 Classical (Kirchhoff) Plate Theory (CPT) ................................ .......................... 44 4.2.1 Assumptions ................................ ................................ ............................ 44 4.2.2 Strain displacement Relationship ................................ ............................ 45 4.2.3 Governing Equations ................................ ................................ ............... 46 PAGE 6 6 4.3 Mi ndlin Plate Theory (First order Shear Deformation Theory) (FSDT) .............. 46 4.3.1 Assumptions ................................ ................................ ............................ 46 4.3.2 Strain displacement Relationship ................................ ............................ 47 4.3.3 Governing Equations ................................ ................................ ............... 48 4.3.4 Constitutive Relationship ................................ ................................ ......... 51 4.4 Analytical a nd Exact Solution ................................ ................................ ............ 52 4.4.1 Cantilever Plate ................................ ................................ ....................... 53 4.4.2 Squa re Plate ................................ ................................ ............................ 54 4.4.3 Circular Plate ................................ ................................ ........................... 58 4.4.4 30 degree Skew Plate ................................ ................................ ............. 60 4.4.5 60 degree Skew Plate ................................ ................................ ............. 61 5 MIXED FORMULATION FOR MINDLIN PLATE ................................ ..................... 66 5.1 Mixed Form ................................ ................................ ................................ ....... 66 5.2 Discrete Collocation Constraints Method ................................ .......................... 69 5.3 Applying EBC Using Implicit Boundary Method ................................ ................ 73 5.4 Numerical Results ................................ ................................ ............................. 78 5.4.1 Cant ilever Plate ................................ ................................ ....................... 78 5.4.2 Square Plate ................................ ................................ ............................ 79 5.4.3 Circular Plate ................................ ................................ ........................... 79 5.4.4 30 degree Skew Plate ................................ ................................ ............. 80 5.4.5 60 degree Skew Plate ................................ ................................ ............. 80 5.4.6 Flange Plate ................................ ................................ ............................ 81 5.5 Concluding Remarks ................................ ................................ ......................... 81 6 MIXED FORMULATION FOR 2D MINDLIN SHELL ................................ ............. 101 6.1 Governing Equations ................................ ................................ ...................... 101 6.2 Mixed Formulation ................................ ................................ .......................... 103 6.3 Numerical Examples and Results ................................ ................................ ... 106 6.3.1 60 degree Skew Plate ................................ ................................ ........... 106 6.3.2 Square Plate ................................ ................................ .......................... 106 6.4 Concluding Remarks ................................ ................................ ....................... 107 7 CONCLUSION ................................ ................................ ................................ ...... 111 7.1 Summary ................................ ................................ ................................ ........ 111 7.2 Scope of Future Work ................................ ................................ ..................... 112 APPENDIX A VOLUMETRIC LOCKING AND SHEAR LOCKING ................................ .............. 114 A.1 Volumetric locking ................................ ................................ .......................... 114 A.2 Shea r locking ................................ ................................ ................................ .. 114 B EQUILIBRIUM EQUATIONS OF 3D ELASTOSTATIC CASE .............................. 116 PAGE 7 7 C DERIVATION OF SHEAR CORRECTION FACTOR ................................ ............ 118 D DERIVATION OF THE JACOBIAN MATRIX IN IBFEM ................................ ........ 121 E FORMULATION OF MINDLIN PLATE ELEMENTS ................................ ............. 123 E.1 Element Q4D4 ................................ ................................ ................................ 123 E.2 Element Q5D6 ................................ ................................ ................................ 126 E.3 Element Q8D8 ................................ ................................ ................................ 128 E.4 Element Q9 D12 ................................ ................................ .............................. 130 E.5 Element Q16D24 ................................ ................................ ............................ 132 LIST OF REFERENCES ................................ ................................ ............................. 136 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 138 PAGE 8 8 LIST OF TABLES Table page 5 1 Location of three interpolation variables and the associated count conditions for patch test ................................ ................................ ................................ ....... 82 5 2 Cantilever plate (Shear force applied at free end) ................................ .............. 83 5 3 Cantilever plate (Bending moment applied at free end) ................................ ...... 83 5 4 Uniformly loaded, clamped square plate [a/t = 10] ................................ ............. 83 5 5 Uniformly loaded, clamped square plate [a/t = 100] ................................ ........... 83 5 6 Uniformly loaded, simply supported square plate [a/t = 10] ................................ 84 5 7 Uniformly loaded, simply supported square plate [a/t = 100] .............................. 84 5 8 Uniformly loaded, clamped circular plate [a/t = 10] ................................ ............. 84 5 9 Uniformly loaded, clamped circular plate [a/t = 100] ................................ ........... 84 5 10 Uniformly loaded, simply supported circular plate [a/t = 10] ............................... 85 5 11 Uniformly loaded, simply supported circular plate [a/t = 100] ............................. 85 5 12 Uniformly loaded, clamped 30 degree skew plate [a/t = 10] ............................... 85 5 13 Uniformly loaded, clamped 30 degree skew plate [a/t = 100] ............................. 85 5 14 Uniformly loaded, simply supported 30 degree skew plate [a/t = 10] ................. 86 5 15 Uniformly loaded, simply supported 30 degree skew plate [a/t = 100] ............... 86 5 16 Uniformly loaded, clamped 60 degree skew plate [a/t = 10] ............................... 86 5 17 Uniformly loaded, clamped 60 degree skew plate [a/t = 100] ............................. 86 5 18 Uniformly loaded, simply supported 60 degree skew plate [a/t = 10] ................. 87 5 19 Uniformly loaded, simply supported 60 degree skew plate [a/t = 100] ............... 87 5 20 Uniformly loaded, arbitrary shape plate ................................ .............................. 87 6 1 Transverse deflection of 60 degree skew plate with one edge clamped .......... 107 6 2 In plane displacement of 60 degree skew plate with one edge clamped .......... 107 PAGE 9 9 6 3 Transverse deflection of square plate with one edge clamped ......................... 107 6 4 In plane displacement of square plate with one edge clamped ........................ 108 PAGE 10 10 LIST OF FIGURES Figure page 2 1 Conforming mesh in traditional FEM ................................ ................................ .. 24 2 2 Scattered nodes in meshless methods ................................ ............................... 24 2 3 Nonconforming structured mesh in IBFEM ................................ ......................... 24 2 4 Step function configuration in IBFEM ................................ ................................ 25 3 1 Geometry of the 2D bracket ................................ ................................ ................ 40 3 2 Transverse displacement dis tribution after deformation (Q9M 130x80 mesh density) ................................ ................................ ................................ ............... 40 3 3 Maximum transverse displacement w.r.t Pois ......................... 40 3 4 ......................... 41 3 5 Geometry of 3D beam ................................ ................................ ........................ 41 3 6 Tran sverse displacement distribution after deformation (H exa8M 65x10x10 mesh density) ................................ ................................ ................................ ..... 41 3 7 Maximum transverse displacement w.r.t Pois .................... 42 3 8 .................. 42 4 1 Configuration for CPT ................................ ................................ ......................... 62 4 2 Configuration for FSDT ................................ ................................ ....................... 63 4 3 Definitions of variables for plate approximations ................................ ................ 63 4 4 Geometry of cantilever plate ................................ ................................ ............... 64 4 5 Geometry of square plate ................................ ................................ ................... 64 4 6 Geometry of ci rcular plate ................................ ................................ .................. 64 4 7 Geometry of 30 degree skew plate ................................ ................................ ..... 65 4 8 Geometry of 60 degree skew plate ................................ ................................ ..... 65 5 1 A typical background mesh using 20x2 Q4 element ................................ ........... 87 PAGE 11 11 5 2 Distribution of transverse displacement after deformation for 100x10 Q4 element (L1/t = 100) ................................ ................................ ........................... 88 5 3 Convergence of total strain energy for cantilever when shear applied (L1/t = 10) ................................ ................................ ................................ ...................... 88 5 4 Convergence of total strain energy for cantilever when shear applied (L1/t = 100) ................................ ................................ ................................ .................... 88 5 5 Distribution of transverse displacement after deformation for 100x10 Q4 element (L1/t = 100) ................................ ................................ ........................... 89 5 6 Convergence of total strain energy for cantilever when bending moment applied (L1/t = 10) ................................ ................................ .............................. 89 5 7 Convergence of total strain energy for cantilever when bending moment applied (L1/t = 100) ................................ ................................ ............................ 89 5 8 A typical background mesh using 10x10 Q9 element ................................ ......... 90 5 9 A typical background mesh using 10x10 Q9 element ................................ ......... 90 5 10 Distribution of transverse displacement after deformation for 150x150 Q9 element (t = 0.1) ................................ ................................ ................................ 90 5 11 Distribution of transverse displacement after deformation for 225x225 Q9 element (t = 0.1) ................................ ................................ ................................ 91 5 12 Convergence of total strain energy for clamped square plate (a/t = 10) ............. 91 5 13 Convergence of total strain energy for clamped square plate (a/t = 100) ........... 91 5 14 Convergence of total strain energy for simply supported square plate (a/t = 10) ................................ ................................ ................................ ...................... 92 5 15 Convergence of total strain energy for simply supported square plate (a/t = 100) ................................ ................................ ................................ .................... 92 5 16 A typical background mesh using 10x10 Q4 element ................................ ......... 92 5 17 Distribution of transverse displacement after deformation for 150x150 Q4 element (t = 1) ................................ ................................ ................................ .... 93 5 18 Convergence of total strain energy for clamped circular plate (a/t = 10) ............ 93 5 19 Convergence of total strain energy for clamped circular plate (a/t = 100) .......... 93 5 20 Convergence of total strain energy for simply supported circular plate (a/t = 10) ................................ ................................ ................................ ...................... 94 PAGE 12 12 5 21 Convergence of total strain energy for simply supported circular plate (a/t = 100) ................................ ................................ ................................ .................... 94 5 22 A typical background mesh using 10x10 Q8 element ................................ ......... 94 5 23 Distribution of transverse displacement after deformation for 200x200 Q8 element (t = 0.1) ................................ ................................ ................................ 95 5 24 Distribution of bending moment using 200x200 Q4 element (t = 0.1) ................. 95 5 25 Convergence of total strain energy for clamped 30 degree skew plate (a/t = 10) ................................ ................................ ................................ ...................... 95 5 26 Convergence of total strain energy for clamped 30 degree skew plate (a/t = 100) ................................ ................................ ................................ .................... 96 5 27 Convergence of total strain energy for simply supported 30 degree skew plate (a/t = 10) ................................ ................................ ................................ .... 96 5 28 Convergence of total strain energy for simply supported 30 degree skew plate (a/t = 100) ................................ ................................ ................................ .. 96 5 29 A typical background mesh using 10x10 Q16 element ................................ ....... 97 5 30 Distribution of transverse displacement after deformation for 200x200 Q9 element (t = 1) ................................ ................................ ................................ .... 97 5 31 Distribution of bending moment using 200x200 Q4 element (t = 0.1) ................. 97 5 32 Convergence of total strain energy for clamped 60 degree skew plate (a/t = 10) ................................ ................................ ................................ ...................... 98 5 33 Convergence of total strain energy for clamped 60 degree skew plate (a/t = 100) ................................ ................................ ................................ .................... 98 5 34 Converg ence of total strain energy for simply supported 60 degree skew plate (a/t = 10) ................................ ................................ ................................ .... 98 5 35 Convergence of total strain energy for simp ly supported 60 degree skew plate (a/t = 100) ................................ ................................ ................................ .. 99 5 36 Geometry of flange plate ................................ ................................ .................... 99 5 37 A typical background mesh using 20x20 Q9 element ................................ ......... 99 5 38 Distribution of transverse displacement after deformation for 150x150 Q9 element (t = 0.1) ................................ ................................ ............................... 100 5 39 Conve rgence of total strain energy for arbitrary shape plate (a/t = 10) ............. 100 PAGE 13 13 5 40 Convergence of total strain energy for arbitrary shape plate (a/t = 100) ........... 100 6 1 Convergence of total strain energy for 60 degree skew plate (a/t = 10) ........... 108 6 2 Convergence of total strain energy for 60 degree skew plate (a/t = 100) ......... 108 6 3 Dist ribution of transverse deflection after deformation for 60 degree skew plate using 150x150 Q4 element (t = 1) ................................ ............................ 109 6 4 Geometry of the 60 degree skew plate ................................ ............................. 109 6 5 Convergence of total strain energy for square plate (a/t = 10) .......................... 109 6 6 Convergence of total strain energy for square plate (a/t = 100) ........................ 110 6 7 Distribution of transverse deflection after deformation for square plate using 100x100 Q9 element (t = 0.1) ................................ ................................ ........... 110 6 8 The geometry of the square plate ................................ ................................ ..... 110 B 1 Stresses notations and directions ................................ ................................ ..... 116 C 1 Distribution of transverse shear stress through the thickness ........................... 119 D 1 Global coordinate and Local coordinate in IBFEM ................................ ............ 122 E 1 Collocation constraints on a 4 node Lagrange element ................................ .... 123 E 2 Interpolation nodes for Q4D4 element ................................ .............................. 124 E 3 Collocation constraints on a 5 node Serendipity element ................................ 126 E 4 Interpolation nodes for Q5D6 element ................................ .............................. 127 E 5 Collocation constraints on an 8 node Serendipity element ............................... 128 E 6 Interpolation nodes for Q8D8 element ................................ .............................. 128 E 7 Collocation constraints on a 9 node Lagrange element ................................ .... 130 E 8 Interpolation nodes for Q9D12 element ................................ ............................ 131 E 9 Collocation constraints on a 16 node Lagrange element ................................ .. 132 E 10 Interpolation nodes for Q16D24 element ................................ .......................... 135 PAGE 14 14 LIST OF ABBREVIATION S C3D8 H 8 node h ybrid hexahedral element constant pressure C3D20H 20 node h ybrid hexahedral element linear pressure CPE4H 4 node h ybrid plane strain quadrilateral element constant pressure CPE8H 8 node h ybrid plane strain quadrilateral element linear pressure CPT Classical Plate Theory EBC Essential Boundary Condition FSDT First Order Shear Deformable Theory H27 27 node hexahedral element H27M 27 node m ixed hexahedral element H8 8 node hexahedral element H8M 8 node m ixed hexahedral element IBFEM Implicit Boundary Finite Element Method LHS Left Hand Side ODE Ordinary Differential Equation Q4 4 node quadrilateral element Q4M 4 node m ixed quadrilateral element Q8 8 node quadrilateral element Q9 9 node quadrilateral element Q9M 9 node m ixed quadrilateral element Q16 16 node quadrilateral element RHS Right Hand Side S4R 4 node doubly curved thin or thick s hell element with reduced integration hourglass control and finite membrane strains S8R5 8 node doubly curved thin s hell elem ent with reduced integration using 5 degrees of freedom per node PAGE 15 15 Abstract of Thesis Prese nted to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science MIXED FORMULATION US IN G IMPLICIT BOUNDARY FINITE ELEMENT METHOD By Hailong Chen May 2012 Chair: Ashok V. Kumar Major: Mechanical Engineering Mixed formulation analysis was proposed for the purpose of avoid ing locking phenomena that occurs in displacement based finite element analysis. In displacement based analysis, volumetric locking will inevitably happen when the material is almost incompressible and the which results in a n infinite bulk modulus S hear locking occurs in Mindlin plate formulation when the plate is very thin but the shear strain s in plate do not go to zero due to the limitation of the interpolation functions. Aside from mixed formulation, several other techniques have also been proposed i n last three decades, such as reduced integration or selective reduced integration method, assumed natural strain method. Implicit Boundary F init e Element M ethod (IBFEM) is a mesh independent finite element method which is motivated by the desire to avoi d mesh generation difficulties in the traditional finite element method ( FEM ) Instead of generating a conforming mesh, a background mesh that does not represent the geometry is constructed for interpolating or approximating the trail and test functions. T he geometry of the model is exactly represented using equations obtained from CAD software T he essential PAGE 16 16 boundary condition is imposed by using implicit boundary method, which uses equations of the boundary and does not need to have nodes on the boundary IBFEM has been demonstrated for 2D and 3D displacement based structural analysis In this thesis, the main goal is to extend this approach to structural analysis using mixed formulation to eliminate volumetric lo cking and shear locking A three field mix ed formulation for incompressible media analysis and a two field mixed formulation for Mindlin plate s are us ed in this thesis A Mindlin 2D shell, is also discussed can model in plane strains as well as bending and shear Several benchmark problems are ut i lized to evaluate the validity of th is approach. PAGE 17 17 CHAPTER 1 INTRODUCTION 1.1 Overview The Finite Element M ethod (FEM) is a widely used numerical method solving problems arising in the engineering analysis. Mesh generation is the first necessary step in t raditional FEM and mesh generation algorithms have been developed that work acceptably for 2D problems but are still unreliable for comp licated 3D geometries, often leads to poor or distorted elements in some regions. Mesh generation is therefore often the most challenging process in the analysis In simulation of failure processes, due to the propagation of cracks with arbitrary and compl ex paths mesh regeneration is needed in each step in traditional FEM. It becomes even more challenge because of the discontinuity and complicated growing path of the cracks Due to aforementioned disadvantages in traditional FEM, there are challenges in i ts application in other field s as well such as manufacturing processes and fluid mechanics In order to better overcome these disadvantage s of traditional FEM a number of meshless or mesh free analysis techniques have been p roposed in last three decades. Meshless method s use a scattered set of nodes for the analysis but the nodes are not connected to form elements (Figure 2 2) Based on the method used to construc t a meshless approximation for the trial and test functio ns various meshless m ethod s exist One of the popular meshless approximation schemes is based on moving least squares method. Some other approaches are also used, such as kernels method and partition of unity method etc. [1] Most methods used to represent trial functions for the meshless approach do not have Kronecker delta properties, which results in difficult y to apply boundary conditions precisely along the boundary. PAGE 18 18 An alternative approach to reduce mesh generation difficulties is to use nonconforming mesh, often a structured background mesh to interpolate or approximate functions in the analysis domain. This approach was first proposed by Kantorovich and Krylov [10] A typical solution structure for applying e ssential boundary conditions is where is the impl icit equation of t he boundary and is the prescribed essential boundary condition is the unknown function that is interpolated piecewise over a mesh Several approaches were used to construct the implicit equation. Rvachev and S hieko [15] have developed an R function to construct a single implicit equation All boundary conditions including essential, natural, and convection boundary conditions are guaranteed in the solution structures. Belytschko et al. [4] has proposed extended finite element method ( X FEM ) bas ed on a structured mesh and implicit boundary representation to remove mesh generation process In X FEM, a pproximat e implicit function of the model was constructed by fitting a set of sample points on the boundary. Radial basis function was used for the i mplicit equation construction. Clark and Anderson [5] have used the penalty method to satisfy the prescribed EBC Another mesh independent method, the Implicit Bo undary F ini te Element M ethod (IBFEM) [11] [13] also utilizes a structured background mesh to interpolate or approximate the trial and test fun ction The geometry of the model is exactly represented by the equatio ns as exported from CAD s ystem A solution structure, similar to the one developed by Kantorovich and Krylov [10 ] is constructed to guarantee the EBC Thi s approach has been tested to be v alid for 2D and 3D structural PAGE 19 19 displacement based analysis In this thesis, we extend this approach to three field mixed formulation for near ly incompressible media and two field mixed formulation for Mindlin plate theory, both pure bending and combination cases 1.2 Goals and Objectives The goal of this thesis is to implement mixed formulation using Implicit Boundary Finite Element Method, so as to avoid volumetric locking for nearly incompressible media and shear locking in thin Mindlin plate s Volumetric locking and shear locking are the most common numerical phenomena that occur in the displacement base d finite element analysis For last three decades, various finite element techniques have been proposed to take care of these problems, such as reduced integration or selective reduced i ntegr ation method, mixed/hybrid method, assumed natural strain method, enhanced assumed s train m ethod, etc. [8] [20] In this thesis, w e will employ the m ixed formulation to remove these locking phenomena using Implicit Boundary Finite Element Method The main objectives of this thes is are 1. Extension of IBFEM to three field mixed formulation for near incompressible media ; 2. Extension of IBFEM to mixed formulation f or Mindlin plate s, pure bending ; 3. Extension of IBFEM to 2D Mindlin shells including both bending and in plane stretching 1.3 Outlines The rest of this thesis is or ganized as follo ws: In Chapter 1, a brief overview about FEM is presented, and the goals and objectives of this thesis are clearly stated. In Chapter 2, we give some details and reference to the m eshless and mesh i ndependent finite element method s. Two PAGE 20 20 challenges while using Implicit Boundary Finite Element Method are described and scheme to solve these issues is presented. I n Chapter 3, a three field mixed formulation for plane strain and 3D stress using IBFEM so a s to remove the volumetric locking for near incompressible media is derived in details. Two examples are used to test the performance of this method In Chapter 4, before proceeding to m ixed formulatio n for Mindlin plates, we review the two classical plat e theories, Kirchhoff Love plate theory and Mindlin plate theory. The analytical or exact solution of severa l benchmark problems, which used later to evaluate t he performance of IBFEM are given in Chapter 4 too The two field mixed formulation using Mindl in plate t heory is detailed in Chapter 5. The focus is on the derivation of m ixed form and imposing EBC using IBFEM An extension of Mindlin plate including both in plane stretching and bending is presented in Chapter 6. Two example s have been used in Chapter 5, the same geometry but different boundary conditions, are employed two test the performance of IB FEM for 2D Mindlin shell case The co nclusion and future work is in Chapter 7. An appendix is also included in order to give more details about some specific topics and also aimed to make this thesis more self contained. PAGE 21 21 CHAPTER 2 MESHLESS AND MESH INDEPENDENT MET HOD 2.1 Traditional FEM Finite Element Method ( [1] [9] [21] ) is a well established numerical technique an d is widely used in solving engineering problems such as stress analysis, heat transfer, fluid flow and electromagnetics in academia as well as in industry. According to Fish and Belytschko [6] the traditional Finite Element Method consists of five procedure s : 1. Preprocessing: subdividing the problem domain into finite elements and approximating the domain by these finite elements mesh generation; 2. Element formulation: derivation of equations in the element level discretization; 3. Element combination: obtaining the equation system for the approximated model from the equations of individual elements assembly ; 4. Solving the equations: using Gauss elimination, Cholesky decomposition or iterative schemes like Gauss Siedel to solve the equation system ; 5. Postprocessing : determining quantities of interest, such as displacement and force resultant, and visualizing the results for future evalu ation analysis Among above five procedures, mesh generation is the most challenge one, and still most endeavors is spend on devising an effective automatic mesh generator for 3D complex geometry, for which the generated mesh is unreliable nowadays. For finite element analysis, the domain of interest is subdivided into small elements by mesh generation techniques and the resulting element mesh approximates the geometry. A typical mesh is shown in Figure 2 1 The mesh is also used to approximate the solution by piece wise interpolation within each element. In procedure 2, employed to convert the strong from into weak form, namely, alleviate the interpo lation sha pe function degree requirement. PAGE 22 22 In order to avoid the problems associated with mesh generation, several approaches have been pr oposed which falls into two cate gories: 1. Mesh less methods; 2. Mesh independent methods Implicit Boundary Finite Element (IBFEM), [11] [13] falls into the mesh indepe ndent method category which utilizes a structured background mesh only for the purpose of interpolation We will present some details about m eshless and mesh independent method in following two section s 2.2 Mesh l ess a nd Mesh I ndependent Method The object ive of meshless methods is to eliminate at least part of mesh generation by constructing the approximat ion entirely in terms of nodes. A set of scattered nodes is used to construct the trial and test function see Figure 2 2 Several schemes are developed to approximate these functions, such as moving least square method, kernel method and partition of unity method, etc. [3] Mesh independent analysis is motivated by the desire to utilize accurate geometric models presented by equations rather approximated by mesh while using a background mesh s olely for the purpose of piecewise approximation or interpolation of the trial and test function. A solution structure is proposed by Kantorovich and Krylov [10] where is the impl icit equation of the boundary and is the essential boundary condition is the unknown function that interpolated piecewise over a mesh Rvachev and S hieko [15] have developed an R function to construct a single implicit equation to represent the entire boundary of a solid. PAGE 23 23 2.3 Implicit Boundary Finite Element Method Comparing to traditional finite element method, Implicit Boundary Finite Element Method (IBFEM) is a mesh independent finite element method in which the geometry of the model is exactly presented by the equations as exported from CAD system Although ther e is no geometry approximation in IBFEM a background mesh ( Figure 2 3 ) is still employed but solely for the purpose of approximation or interpolation of the trial and test function. Contrary to the R function technique used in Rvachev and S hieko [15] an approximate step function was used in IBFEM to construct the implicit equation of the interested domain. The solution structure used in IBFEM is (2 1) where is the step function. A typical approximate step function being used in IBFEM is (2 2 ) where can be is the distance between point to the boundary lines in the normal direction, and is the transition width see Figure 2 4 PAGE 24 24 Figure 2 1. Conforming mesh in traditional FEM Figure 2 2. Scattered nodes in meshless methods Figure 2 3. Nonconforming structured mesh in IBFEM PAGE 25 25 Figure 2 4. Step function configuration in IBFEM PAGE 26 26 CHAPTER 3 MIXED FORMULATION FO R NEARLY INCOMPRESSIBLE MEDIA 3.1 Overview It has been frequently noted that in certain constitutive laws, such as those of viscoelasticity and associative plasticity, the material behaves in a nearly inc ompressible manner. The incompressibility of these media at certain critical stage, e.g., metals at the range of plastic deformations ( yielding ), results in a locking phenomenon, named volumetric locking or Poisson locking. This phenomenon result s from the fact that when the material is near incompressible status, tends to 0.5, which makes the bulk modulus t end to wards infinity and hence result in a ill conditioned stiffness matrix in the finite element model Several approaches have been proposed in last three decades to reduce or alleviate volumetric locking occurrence. Reduced or s elective reduced integration technique wa s the first successful irreducible form of solutions for volumetric locking problems, although in the beginning not directed speci ally towards volumetric locking. L ater, other formulations succeeded in using augmented functional, when compared to that obtained from displacement based approaches, incorporatin g additional fields into the formulation and leading to the onset of general mixed methods. In the 1990s, the enhanced strain m ethod was also applied to alleviate this locking phenomenon. Generally, there are two choices of mixed formulation on additional fields, displacement and mean stress which named two field formulation, and displacement mean stress and volume change termed three field formulation ( [1] [9] [21] ) Which of these should be employed may depend on the form of the constitutive equation used. For situations where changes in volume affect only the PAGE 27 27 pressure the two field form can be eas ily used. However, for problems in which the response may become coupled between the deviatoric and mean components of stress and strain the three field formulation leads to much simpler forms from which to develop a finite element model. In this thesis, w e will only focus on the later formulation, three field mixed formulation. The layout of Chapter 3 is as follows: In section 3.2, we derive the general mixed formulation, and the resultant equations using three field formulation; In section 3.3 and 3.4, we specify our discussion to 2D plan strain and 3D stress respectively A further modification of the mixed form s for 4 node and 9 node Lagrange elements in 2D case and 8 node and 27 node hexahedron elements for 3D case ; In section 3.5, we employ some benchm ark problems to test the performance of elements developed in IBFEM. The conclusion of Chapter 3 is made in section 3.6. 3.2 Mixed Form ulation As mentio ned in the introduction of Chapter 3 a three field mixed formulation is adopted in this section to develop some valuable elements both for 2D plane strain and 3D stress using IBFEM. In this section, we will develop the general mixed form for both cases. 3.2.1 Matrix Decomposition Generally, in most cases, the strain and stress matrices can be s p li t into the devi ator (isochoric) and mean parts. Since the separation slightly differs for 2D and 3D. Here, we utilize 3D case for illustration. That for 2D will be presented in section 3.3. Accordingly, the m ean stress pressure part, can be expressed as PAGE 28 28 (3 1) We use sum notation here. And And the d eviator ic part can be defined as ( 3 2 ) where is the Kronecker delta function, which when is one and else zero for all Similarly, the m ean strain v olume change, can be defined as (3 3 ) and the d eviator ic strain as (3 4 ) Have above definitions for mean part and deviatoric part of stress and strain, t he strain and stress may now be expressed in a mixed form as Strain (3 5 ) Stress ( 3 6 ) where the mean matrix operator (3 7 ) and deviator ic matrix operator PAGE 29 29 (3 8 ) the operator (3 9 ) is the identity matrix. is the set of stress obtained directly from the strain rates, depending on the particular constitutive model form. 3.2.2 W eak F orm The governing equations of this topic are the same as those equilibrium equations of static elastic problems, which only differs at which equation we should employ for 2D plane strain and 3D stress scenarios. For detailed derivation of these governing equation s one can refer to Appendix B. use principle of minimization of potential energy to achieve t he same purpose. The strong form is (3 10 ) (3 11 ) PAGE 30 30 (3 12 ) go with natural and essential boundary conditions that on (3 13 ) on (3 14 ) and simplifying the results we can have the weak form as (3 15 ) (3 16 ) ( 3 17 ) Introducing the finite element approximations of variables as (3 18 ) And similar approximati ons to virtual quantities as (3 19 ) Hence, t he strain in an element becomes ( 3 20 ) In which is the standard strain displacement matrix. Similarly, the stresses in each element may be computed by using (3 21 ) PAGE 31 31 Since strain and stress are the element variables, we will keep it as constant while constructing the weak form system. Su bstitute Equations 3 18 and 3 19 into the weak f orm E quations 3 15 to 3 17 we have F or Equation 3 15 : (3 22 ) F or Equation 3 16 : (3 23 ) F or Equation 3 17 : (3 24 ) Writing Equations 3 22 to 3 24 in matrix form, t he finite element equation s ystem become s (3 25 ) Where PAGE 32 32 Above equation system cannot be solve in a global sense, since the stress in is not a global variable and not directly approximated. If the pressure and volumetric strain approximations are taken locally in each element it is possible to solve the above second and third equation in each element individually. If we further make that in each element, the array is now symmetric positive definit e. W e can use the second equation in 3 25 to solve for and in each element as (3 26 ) (3 27 ) The mixed strain in each element may now be computed as (3 28 ) Where defines a mixed form of the volumetric strain displacement equation. After solving the first two Equations of 3 25 the first equation can be write in a n alternative form as (3 29 ) where is the general stress strain matrix before using mixed formulation. Equation 3 29 can be formatted in the same way as general displacement based formulation. For impleme ntation pu rpose, we rewrite Equation 3 29 as PAGE 33 33 where is the same as defined in Equation 3 25 We have detailed how to get the mixed weak form to get rid of volumetric locking for both 2D plane strain and 3D stress. More attention should be focused on the solution of the mixed equation system. Since an element variable, element stress is involve d in the first equation of 3 25 it is impossible to globally solve this equation system. A n e fficient approach is adopted in order to achieve solvability of the equation system, that we solving the se cond and th ird equation of 3 25 in each element individually. The last part of section 3.2.2 is focused on h ow to use this approach to solve these two equations, and finally solve the whole system. Note should be made here that alternatives is available to solve ab ove mixed equation syst em, 3 25 Presented is just one possible approach 3.3 2D Plane Strain In thi s section, we will provide detail s on th e element variable s needed for 2D plane strain case For 2D Plain Strain, the matrix can be expressed using the shape functions of displacements as PAGE 34 34 The shape functions for displacement is the traditional Lagrange shape functions. While, a ccording to [1] [9] and [21] for 4 node Lagrange element the shape function needed for volumetric strain and pressure are And that for 9 node Lagrange element are The matrix can be immediately obtained using above shape functions for mean strain and stress. And the matrix can be simplified as PAGE 35 35 The matrix only dif fers from that provided in Equation 3 7 in dimension. where The elastic stress strain matrix is The mixed stress strain matrix is 3.4 3D Stress In this section, we continue our elaboration for 3D stress case. PAGE 36 36 For 3D stress the matrix can be expended using displacement shape functions as The shape functions for displacement is th e traditional Lagrange shape functions. While, according to [1] [9] and [21] the shape functions for 8 node hexahedron element of the mean strain and stress are And f or 27 node hexahedron element, the shape functions are The matrix can be obtained as And the matrix is PAGE 37 37 The matrix and matrix are already given in Equations 3 7 and 3 8 as The elastic stress strain matrix for 3D stress is PAGE 38 38 Hence, the mixed stress strain matrix is 3.5 Numerical Examples and R esults In this section, two examples are used to test the performance of these elements using IBFEM. The first example is a fixed bracket, whi ch also has been used in [1] 3D c lamped beam is utilized for 3D stress analysis Also, we compare our results with element CPE4 H and CPE8 H for plane strain C3D8H and C3D20 H for 3D s tress, which are all from ABAQUS pack age. Instead of tabulating those data, we give related plots to help visualizing the volumetric locking and the difference between different elements 3.5.1 Bracket (Plane strain) The bracket is fixed at the two vertices, and a uniformly distributed load o f value Th e geometry of this bracket is given in Figure 3 1 The distribution of transverse displacement a fter de formation analyzed by Q9M is in Figure 3 2 The compa rison between Q4 element and Q4M PAGE 39 39 element, Q9 element and Q9M element is plotted in Figure 3 3 and Figure 3 4 respectively. From Figure s 3 3 and 3 4 volumetric locking is very severe in Q4 element, while with the incre ase of degree of shape function, using Q9 element, this phenomenon can be better improved, but still cannot be avoid ed See Figure 3 4 3.5. 2 Beam (3D stress) A beam cl amp ed at the left end is loaded with a distributed load of 6000 at the top The geometry of this beam is shown in Figure 3 5 The defor med shape of the beam with transverse displacement distribution is given in Figure 3 6 The comparison between different elements also given in this example, Figure 3 7 and Figure 3 8 From Figure s 3 7 and 3 8 volumetric locking is also drastic in H8 element, while with the increase of de gree of shape function, using H27 element, this phenomenon is better improved, but st ill cannot avoid. See Figure 3 8 3.6 Concluding R emarks A three field mixed formulation scheme has been adopted in Chapter 3 to remove or alleviat e the volumetric locking in 2D plane strain and 3D stress. The detailed derivation of such formulation is presented in section 3.1, 3.2 and 3.3 Two examples are employed to test the validity of the developed elements using IBFEM. The comparison between di splacement based formulation and the three field mixed formulation is also presented in these two examples. Elements from Abaqus package are also used in the comparison. From the results and plots given these IBFEM elements can better remove or alleviate volumetric locking phenomenon in bot h 2D plane strain and 3D stress, and performs the same, if not better than, as the Abaqus elements. PAGE 40 40 Figure 3 1. Geometry of the 2D bracket Figure 3 2. Transverse displacement dis tribution after deformation (Q9M 130x80 mesh density) Figure 3 PAGE 41 41 Figure 3 Figure 3 5. Geometry of 3D beam Figure 3 6. Transverse displacement distribution after deformation (H exa8M 65x10x10 mesh density) PAGE 42 42 Figure 3 Figure 3 PAGE 43 43 CHAPTER 4 CLASSICAL PLAT E THEORIES 4.1 Overview The plate element has attracted considerable amount of attentions due to its great application in engineering field. Among numerous plate theories that have been developed since the late nineteenth century, the first widely accepted and used in engineering was proposed by Kirchhoff and developed by Love, named Kirchhoff Love (KL) plate theory, which is an extension of Euler Bernoulli beam theory. In KL plate ormals remain straight and perpendicular to the mid plane of the plate before and after deformation, and the strain in the thickness direction is negligible. In other word, the strains are negligibly zero in KL plate theory. These assumptions render KL plate theory only applicable for thin structures and C 1 continuity is required for shape functions. Mindlin and Reissner have relaxed Kirchhoff assumption to include the shear flexibility in the plate theory, which called Mindlin Rei ssner (MR) plate theory or First order shear deformation theory. In MR plate theory, the normal not necessary perpendicular to the mid plane and the transverse shear strains should be taken into consideration. Inclusion of the shear flexibility in the MR p late theory makes it comparatively suitable for moderately thick plates and only C 0 continuity is required. However, MR plate elements exhibit a phenomenon termed shear locking when the thickness of the plate tends to zero. Over last three decades, extens ive efforts have been spent to device an effective finite element scheme to overcome shear locking phenomenon for MR plate element. A lot of finite element techniques have been developed, including but not limited to: Reduced Integration or Selective Reduc ed Integration method, Mixed/Hybrid method, PAGE 44 44 Assumed Natural Strain method, Enhanced Assumed Strain method, Discrete Shear Gap method ( [8] [20] ) Recently, a large number of endeavors were focused on developing new plate/shell elements u sing various new finite element techniques, such as iso geometric method, smoothed finite element method, discontinuous Galerkin finite element method Mesh free finite element method, implicit boundary finite element method, using twist Kirchhoff theory i n finite element method, etc. Implicit Boundary Finite Element Method (IBFEM) is a mesh independent Finite Element Method. In IBFEM, the Dirichlet boundary conditions are imposed using implicit boundary method where approximate step functions are used as weighting functions to construct solution structures that enforce the boundary conditions. This method can impose boundary conditions even when the boundary is not guaranteed to have nodes on them or when the shape functions used for the analysis do not sa tisfy Kronecker delta property. Previous work by Kumar et al. [11] [13] has illustrated its validity and potential value by employing convent ional Lagrange elements and B Spline elements in their analysis. In Chapter 5 we use the discrete shear collocation method and IBFEM to develop a family of Mindlin plate elements which can avoid shear locking in thin plates 4.2 Classical ( Kirchhoff ) Plate Theory (CPT) 4.2.1 Assumption s There are two assumptions for CPT: 1. Transverse normals remain straight and perpendicular to the mid plane before and after deformation; 2. Strain in the tra nsverse normal direction is negligible. PAGE 45 45 4.2.2 Strain displacement R elation ship As depicted in Figure 4 1 taking the plane of the Cartesian coordinate system to coincide with the mid plane of the plate, we can have the plate displacement components represented as: (4 1 ) (4 2 ) (4 3 ) where are regarded as the weighted average for the deflection, is the coordinate in the thickness direction. Using above assumed displacements, the in plane strains can be deduced by the elasticity definition of strain components as following: (4 4 ) Usually, the matrix is termed as bending curvature matrix, which actually a pseudo curvature matrix. Transverse shear strains, which actually are the appropriate averages, are expressed in terms of deflection as PAGE 46 46 (4 5 ) The is zero according to the assumption for Kirchhoff plate 4.2.3 Governing E quations Governing equations for CPT can be derived from equilibrium equat ions for static elastic system. Since these equations are the same for both Kirchhoff plate theory and Mindlin plate theory, we only provide the result here. Details on notation and derivation will be elaborated in Mindlin plate section. T he governing equations of CPT are : (4 6 ) (4 7 ) (4 8 ) So far, we have finished our brief discussion about Kirchhoff plate theory. From now on, we will focus on the Mindlin plate theory and the development of valuable Mindlin plate element using Implicit Boundary Method. Some content provided without any derivation will be detailed in next section. 4.3 Mindlin P late Theory (First order Shear D eforma tion T heory) (FSD T) 4.3.1 Assumption s All assumption made in Kir chhoff plate theory is applicable in Mindlin plate theory, except one, that a ll plane sections normal to mid plane remain plane, but no necessary PAGE 47 47 normal to the mid plane after deformation. The configuration of Mindlin plate is shown in Figure 4 2 4.3.2 Strain displacement R elation ship Assuming there is no in plane forces acting on the plate, the local displacement in the directions of the and axes are presented as (4 9 ) (4 1 0 ) (4 11 ) Where are the rotation of the normals in the and direction s respectively. And this is the very distinction between CPT and FSDT. Immediately the strains in the and axes are available as In plane strains: (4 12 ) Transverse strains: (4 13 ) (4 14 ) PAGE 48 48 This condition will be generally valid when we apply transverse load distributed over a large area such that the normal stress is negligible. However, in situations where a very small area compared to the plate thickness is applied a transverse load, surface indentations may be created and the assumption for transverse displacement will not be valid. Lo w velocity impact due to e xternal objects may incur such phenomenon. In above expression for strains, we have adopted a notation that 4.3.3 Governing E quation s We will detail on the derivation of governing equations for both Kirchhoff plate and Mindlin plate in this section. The differential equilibrium equations of a static small deformation solid with body force are: (4 15 ) (4 16 ) (4 17 ) In above expansion, we omit the term, since we have assumed at the very beginning that this stress in negligible small These three equations are not satisfied at every point of instead we make them satisfied in an average sense through the PAGE 49 49 thickness of the pla te. This is accomplished by requiring the equations and their moments be satisfied on both side of the equation on an average sense. The notation used is shown in Figure 4 3 Integrating the E quation 4 15 in the thickness direction, we have (4 18 ) Where Multiply E quatio n 4 15 with and i ntegrate in the thickness direction, we get (4 19 ) Where and Integrating the E quation 4 16 in the thickness direction, we obtain (4 20 ) PAGE 50 50 Where Multiply E quatio n 4 16 with and integrate in the thickness direction, we get (4 21 ) Where and Integrating the E quation 4 17 in the thickness direction, we obtain (4 22 ) Where Multiply E quatio n 4 17 with and integrate in the thickness direction, we get (4 23 ) PAGE 51 51 Where and Omitting the in plane terms in E quatio ns 4 18 to 4 23 we finally have the governing equations as (4 24 ) (4 25 ) (4 26 ) 4.3.4 Constitutive R elationship In this section, we will derive the relationship between the plate resultants and the displacement and rotations. As i n last section, the shell resultants have been defined as Using the str ess strain relationship and the strain displacement relation, we can write the resultants in term of displacement and rotations as (4 27 ) Where PAGE 52 52 (4 28 ) where Constant is added here to account for the fact that shear stresses is the exact for rectangular, homogeneous section and corresponds to a parabolic shear stress distr ibution. Details about how to derive this constant can be found in Appendix C In next section, the analytical solution using above governing equation for various Kirchhoff plate and Mindlin plate is presented. The purpose of this work is to compare the an alytical result with numerical result using developed Mindlin plate elements. Hence, evaluate the performance of those elements. 4.4 Analytical a nd Exact S olution Before proceeding to mixed formulation of Mindlin plate problem, the analytical or exact sol ution s for those benchmark problems are derived for the purpose of comparing results and hence testing the validity of our numerical solution. Since this thesis is mainly focus on the finite element analysis of various kinds of Mindlin plates we provide the analytical or exact solution without any further derivation. For more details, one can refer to standard textbooks on plates and shells, such as Timoshenko et al. [18] and Reddy et al. [19] The benchmark problems will be studied here are : cantilever plate, clamped and simply supported square plate, clamped an d simply supported circular plate and clamped and simply supported 30 degree and 60 degree skew plate s For simplicity, we utilize following non dimensional constants for all testing examples. PAGE 53 53 Uniform load p = 100 or bending moment m = 100 ; Notations for geometries are as following: Lx: Length of edge x; a: Edge length; r: radius; t: thickness. The thickness varies among 0.1 and 1, so as to test the validation of the element for different leng th/thickness ratio. 4.4.1 Cantilever P late Cantilever plate is studied here for the purpose of testing the performance of new elements when only shear force or bending moment is applied. The cantilever in Figure 4 4 is clamped at the left edge and constant shear force or bending moment, both with value of 100, is applied separately at the free end. Due to the loading condition and the length to width ratio, Cantilever plate can be viewed as 1 D plate undertaking cylindrical bending. In these cases, when shear force or bending moment applied separately, the transversal deflection is derived here using 1 D assumption. And the res ults from using Timoshenko beam theory are also provided here without derivation. For mo re details, one can refer to Timoshenko et al. [18] S hear force app lied Us ing plate theory: T he deflection of the cantilever is (4 29 ) where is the component of stiffness matrix in row one and column one. Using beam theory: The deflection of the beam using beam theory is PAGE 54 54 (4 30 ) where and is the m oment of inertia Bending moment applied A positive bending moment of 100 is applied at the free end. Using plate theory : T he deflection of the cantilever i s (4 31 ) Using beam theory: From Timoshenko beam theory, we can have the deflection of the beam is (4 32 ) 4.4.2 Square P late Square plate is conventionally used to assess the element performance. Both camped and simply supported cases are studied here. The plate is subjected to a uniformly distri buted load with value of 100 The dimension is shown in Fi gure 4 5 C lamped The general analytical solution for problem of this case is not easy to be obtained using any series solution. But the exact solution of deflection and bending moment s at the center can be found in numerous in literatures. Below, same valuable results are given both for thin and thick case. For details, one can refer to the reference s The transverse displacement and bending moment solution: Thin case : PAGE 55 55 Here, we present several solutions, analytical and exact so as to fully study the pe rformance of elements using IBFEM. Analytical solution by Timoshenko et al. [18] : The central deflection is (4 33 ) The central bending moments are (4 34 ) The exact solution by ZienKiewicz et al. [22] : The central deflection is (4 35 ) The central bending moments are (4 36 ) Exact solution by Taylor et al. [17] : The central deflection is (4 37 ) The central bending moments are (4 38 ) Thick case : Exact solution by ZienKiewicz et al. [22] : The central deflection is PAGE 56 56 (4 39 ) The central bending moments are (4 40 ) S imply supported The transverse displacement and bending moment solution: Thin case : From Timoshenko et al. [18] the deflection of Kirchhoff plate of this case is (4 41 ) where and flexural rigidity The bending moments are (4 42 ) An approximate r esult was also given by Timoshenko et al. [18] as The central deflection is (4 43 ) PAGE 57 57 The central bending moments are (4 44 ) The exact solution by ZienKiewicz et al. [22] : The central deflection is (4 45 ) The central bending moments are (4 46 ) Thick case: According to Reddy et al. [19] the deflection of Mindlin plate of this case is the (4 47 ) where is the Marcus moment of Kirchhoff plate. Of this case, the Marcus moment is given by (4 48 ) Using above deflection relationship between Kirchhoff plate and Mindlin plate of this case, we have the deflection for Mindlin plate as (4 49 ) PAGE 58 58 where The bend ing moments for thick case are (4 50 ) The solution also given by ZienKiewicz et al. [22] as The central deflection is (4 51 ) The central bending moments are (4 52 ) 4.4.3 Circular P late A circular plate under uniformly distributed load with clamped and simply validity of imposing essential boundary condition using Implicit Boundar y Method. The geometry of this circular plate is shown in Figure 4 6 Clamped The transverse displacement and bending moment solution: Thin case : In Timoshenko et al. [18] the deflection in the form as PAGE 59 59 (4 53 ) And the bending moment (4 54 ) Thick case : For simply supported circular plate, according to Reddy et al. [19] the relationship of deflection and bending moments between Mindlin plate and Kirchhoff plate is (4 55 ) where Simply supported The transverse displacement and bending moment solution: Thin case : Fr om Timoshenko et al. [18] the deflection of a simply supported plate is of the form as (4 56 ) The bending moment s are PAGE 60 60 (4 57 ) Thick case : The relation between Mindlin plate and Kirchhoff plate is the same as that in the clamped case: (4 58 ) where 4.4.4 30 d egree Skew P late Skew plate usually being utilized in testing the performance of distorted element. In order to test the validity of the Implicit Boundary Method for imposing essential boundary condition, skew plate is a valuable problem for this purpose. W e directly provi de some results for skew plate according to Liew [14] and Sengupta [16] For more details, one can refer to Liew [14] and Sengupta [16] The ge ometry of this plate is shown in Figure 4 7 C lamped The transverse displacement solution: Thin case : According to Sengupta [16] the deflection at the center is PAGE 61 61 (4 59 ) For thick case, there is no analytical solution or exact solution available in literature. Simply supported The transverse displacement solution: Thi n case: According to Liew [14] et al. the deflection at the center is (4 60 ) Thick case : According to Liew et al. [14] the de flection at the center is (4 61 ) 4.4.5 60 degree S kew P late The geometry of this plate is shown in Figure 4 8 C lamped Transverse displacement solution: Thin case : According to Sengupta [16] the deflection at the center is (4 62 ) For thick case, there is no analytical solution or exact solution available in literature. PAGE 62 62 Simp ly supported The transverse displacement solution: Thin case : According to Liew et al. [14] the deflection at the center is (4 63 ) Thick case : According t o Liew et al. [14] the deflection at the center is (4 64 ) Figure 4 1. Configuration for CPT PAGE 63 63 Figure 4 2. Configuration for FSD T Figure 4 3. Definitions of variables for plate approximations PAGE 64 64 Figure 4 4. Geometry of cantilever plate Figure 4 5. Geometry of square plate Figure 4 6. Geometry of circular plate PAGE 65 65 Figure 4 7. Geometry of 30 degree skew plate Figure 4 8. Geometry of 6 0 degree skew plate PAGE 66 66 CHAPTER 5 MIXED FORMULATION FO R MINDLIN PLATE We will present the derivation of mixed form for Mindlin plate and validate our IBF EM elements in Chapter 5 Details will focus on deriving the mixed weak form and imposing essential boundary condition ( EBC ) using Implicit Boundary Method. The governing equations for Mindlin plate ar e already derived in Chapter 4 So we begin the finite formulation with the matrix form in Chapter 5 5.1 Mixed Form Writing in matrix form, the governing equation system is (5 1 ) (5 2 ) (5 3 ) (5 4 ) Generally, is not included in the plate theory. Equations 5 1 to 5 4 constitute the governing equation s in the plate theory. Eliminating the bending moment we have the three governing equations as: (5 5 ) (5 6 ) PAGE 67 67 (5 7 ) form. (5 8 ) For Equation 5 5 : (5 9 ) is the prescribed shear on the boundary For Equation 5 6 : (5 10 ) is the prescribed moment on the boundary For Equation 5 7 : (5 11 ) We use following approximations for di fferent variables : PAGE 68 68 (5 12 ) Substitute above approximation 5 12 into Equations 5 9 to 5 11 the weak form system can be represented as (5 13 ) Where For fairly thin plate, and where is the shape function for shear As the thin plate limit is approached, there are certain necessary criteria should be satisfied in order to achieve the solution stability According to Zienkiewicz et al. [21] certain necessary but not sufficient count condition s must be satisfied in order to develop useful and robust elements. The count conditions are: (5 14 ) (5 15 ) PAGE 69 69 Where and are the number s of and parameters When this count condition s are not satisfied then the equation system will be either lock ed or singular. Several elements satisfy the se criterions are presented in Table 5 1 In Table 5 1 we can see that the 4 node element ( [2] ) fails the count condition for four element patch test and the 9 node element ( [7] ) fails the count condition for single element patch test but the margin is small. If we increase the element number for the patch test, this count condition w ill be satisfied. The 5 node element, depicted in above table, passes the count condition for both single element and four element patch test, but it results in singular stiffness matrix, which illustrate the fact that the count condition is just a necessa ry, but not sufficient condition to ensure the system solvable and stable. In the next section, we will employ discrete collocation method to impose the above count condition and finally arrive at the mixed weak form. 5.2 Discrete Collocation Constraints Method To some degree, using three field mixed formulation to achieve satisfactory performance is limited. But a different approach uses collocation constraints for the shear approximation on the element boundaries, thus limiting the number of parameters and making the count condition more easily satisfied. We will focus on this method to derive a useful family of Mindlin elements using IBFEM from now on. The shear strain is uniquely determined given and at all the nodes. And hence is also uniquely determined there. Using this relationship and imposing the count conditions, t he prescribed values of shear at certain nodes are PAGE 70 70 where and are two constant matrices for each interpolation type. For simplicity, we directly give this expression here. A detailed derivation is presented in Appendix E. We will re derive the two field mixed weak form by above discrete shear field collocation method in this section Instead of using we adopt the principle of minimum potential energy method to derive the weak form for the static system. The potential energy of the system is (5 16 ) Where (5 17 ) (5 18 ) In order to make the potential energy minimum and hence a stable status of the system we let PAGE 71 71 (5 19 ) Plug ging the relation for the shear term, into equatio n 5 19 and interpolate .u sing their shape function 5 12 respectively, the variation of potential energy becomes (5 20 ) (5 21 ) (5 22 ) In the above deduction, we only introduced the variation of total strain energy terms in t he first two expressions, 5 20 and 5 21 since for the external force term, there are no shear collocation terms, and naturally, it follows a traditional routine. For any arbitrary and above E quation 5 22 must satisfied. Hence we obtain the weak form as (5 23 ) (5 24 ) Rearranging terms in E quatio n s 5 23 and 5 24 and being sorted in matrix form we have the weak form for two field mixed formulation as (5 25 ) PAGE 72 72 Where So far, we have finished the derivation of two field mixed formulation using discrete shear field collocation method. Some post processing issues are also discuss ed in following content. Strains and Stresses After solving the equation system for and at each node, the associated strains and stresses are of our interest sometimes to obtain these quantities for the implementation. (5 26 ) We can calculate the in plane strains at any l ayer on the condition that Equation 5 26 is calculated ahead. The in plane stra ins on the top surface of the plate are (5 27 ) For bottom in plane strains, due to the linearity with respect to thickness it has the same magnitude as those of the top surface of the plate, but in the opposite direction. After computing the in plane strains the corresponding stresses on the top surface of the plate also can be calculated by constitutive relation ship as PAGE 73 73 (5 28 ) where is the isotropic elasticity matrix. The average transverse shear stresses can be computed as (5 29 ) Equation 5 29 for the transver se shear stresses is just in an average sense. Most of times, it is of little use to predict the stress loading condition. In order to get a more accurate transvers e stresses distribution along the thickness, the equilibrium equations of the static system is used for this purpose, just as what we did to derive the shear correction factor in Appendix C. For more details, one can refer to Appe ndix C. Those bending moments and shear forces also can be calculated by (5 30 ) Again, the shear forces are also in an average sense. One can calculate the exact shear force distribution once given the transverse stresses distri bution. 5.3 Applying E B C U sing Implicit Boundary Method Two important issues in developing IBFEM for mixed formulation s is how to constructing the step function and solution structure and imposing the essential boundary condition. Variou s ways of constructing these functions are already presented in the Chapter 1 One can refer to those references for details. In this section, we will focus on the Implicit Boundary Method utilized by Kumar et al. [11] [13] PAGE 74 74 S olution structure and step function The solution structure to ensure the imposed essential boundary condition f or the transverse displacement and rotations are (5 31 ) where is the step function s to be constructed, are the prescribed essential boundary condition, and are the nodal values. can be constructed as (5 32 ) where is the distance between point s to the boundary lines in the bi normal direction, and is the transition width in order to avoid computation problem of sharp change from 0 to1. For transverse displacement, for rotations, Total s train energy In order to impose EBC, the two fi e l d mixed form need s to be modified. W e use the principle of potential energy method to re derive the weak form using IBFEM The potential energy of the system is (5 33 ) In order to make the potential energy minimum hence make the system stable we let the variation of the potential energy to be zero. That is PAGE 75 75 (5 34 ) Rewriting above equation in a matrix form, the equation becomes (5 35 ) For convenience of elaboration we define (5 36 ) (5 37 ) Hence the total strain energy term of E quation 5 35 can be rewritten as (5 38 ) Knowing that and plugging in the trial solution and test function Equation 5 31 the integrant of 5 38 can be expressed as (5 39 ) Where (5 40 ) PAGE 76 76 results from the prescribed essential boundary conditions. For implementation purpose, we can move this term to the RHS, considered to be a force term together wit h the original external force terms. This fictitious force term is (5 41 ) can be rewrite as (5 42 ) T he L HS can be further simplified as follows : Separat ing into two parts as, (5 43 ) Also, we define here (5 44 ) (5 45 ) Above two terms also can be expressed in terms of nodal displacement as: (5 46 ) (5 47 ) Expanding we can have PAGE 77 77 (5 48 ) Apparently, is what we have used in the collocat ion method for inner elements. For boundary elements, there are two more additional terms, that and due to the essential boundary condition. Integration issue For terms not including the step function and we use volume integration. That (5 49 ) (5 50 ) (5 51 ) (5 52 ) For terms contain these two step functions, instead of using volume integration, we can use boundary integration, since the step functions are define as co nstant except the boundary transit region. (5 53 ) (5 54 ) PAGE 78 78 So far, we have successfully imposed the EBC using Implicit Boundary Method. And details about formulation and implementation are given. In next section, we will compare our numerical solution with analytical or exact solution provided in Chapter 4. 5.4 Numerical R esults Numerical results for the benchmark problems in section 4.3 will be studied in this section. S4R and S8R5 eleme nts in commercial package ABAQUS were also used in this analysis to test the performance of the developed elements. Comparison among numerical results, exact solutions and analytical solutions are presented. Within each example, the convergence of total st rain energy is al so studied 5.4 .1 Cantilever Plate A canti lever plate with background mesh is show n in Figure 5 1 Shear applied Shear is applied at the free end of the cantilever plate. T he relevant result of this case is provided in Table 5 2 And the distribution of the transverse displacement after deformation of thin case is presented in Figure 5 2 The convergence plots of the total strain energy for both thin and thick cases are provided in Figure 5 3 and Figure 5 4 respectively. Bending moment applied Only bending is applied at the free end in this case. The relevant result of this case is presented in Table 5 3 both thin and thick cases. The distribution of the transverse displacement after deformation for thin case is presented in Figure 5 5 And the total strain energy plots are in Figure 5 6 and Figure 5 7 PAGE 79 79 5.4. 2 Square Plate We utilize two kinds of background mesh for the interpolation for this very conventional benchmark problem, as shown in Figure 5 8 and Figure 5 9 The results obtained from both these two approximations using IBFEM converges to the exact solution from Zienkiewicz et al. [22] Data tabulated in Table 5 4 to Table 5 7 are collected using background pattern one ( Figure 5 8 ). Fully clamped The rele vant data of this case is given in Table 5 4 and Table 5 5 The deformed transverse displace ment s for both two mesh patterns are provided in Figure 5 10 and Figure 5 11 respectiv ely. And the convergence plots for the total strain energy are in Figure 5 12 and Figu re 5 13 Simply supported The relevant data of this case is given in Table 5 6 and Table 5 7 The convergence plots for the total strain energy is in Figure 5 14 and Figure 5 15 5.4. 3 Circular Plate A typical bac kground mesh for this circular plate is depicted in Figure 5 16 Table 5 8 to Table 5 11 give the relevant testing results. Fully clamped The relevant data of this case is given in Table 5 8 and Table 5 9 And the convergence plot for the total str ain energy is in Figure 5 18 and Figure 5 19 A graph for the transverse displacement of the thick case is given in Figure 5 17 Per iphery is simply supported The relevant data of this case is given in Table 5 10 and Table 5 1 1 And the convergence plot for the total strain energy is in Figure 5 20 and Figure 5 21 PAGE 80 80 5.4 .4 30 degree Skew P late A typical background mesh for this 30 degree skew plate is provided in Figure 5 22 Table 5 12 to Table 5 15 giv e the relevant testing results. Fully clamped We found that the analytical or exact solution for a comparatively thick plate, such as length/thickness ratio is 0.1, is very rare, even impossible to refer to in literature. Due to this limitation, we will not provide any analytical solution or exact sol ution for this case for comparison. We only compare our results with S4R and S8R5 elements. The relevant data are given in Table 5 12 and Table 5 13 The deformed transverse displacement distribution is presented in Figure 5 23 and the convergence plot for the total strain energy is in Figure 5 25 and Figure 5 26 Fully simply supported Since t here is bending moment singularity at the obtuse corners for simply supported thin plate, we choose not to provide the bending moment for this case here, only with one figure, Figure 5 24 to demonstrate the singularity of bending moment at the obtuse corn ers. The relevant data of this case are given in Table 5 14 and Table 5 15 And the conv ergence plot for the total strain energy is in Figure 5 27 and Figure 5 28 5.4 .5 60 degree Skew P late A typical background mesh using 10X10 Q16 element is shown in Figure 5 29 Fully clamped The analytical or exact solution for a comparatively thick skew plate is also very rare, even impossi ble to find in literature For this case, we only compare the result with S4R and S8R5 elements. The relevant data of this case are formatted in Table 5 16 and PAGE 81 81 Table 5 17 And the convergence plot for the total strain energy is in Figure 5 32 and Figure 5 33 The deformed transverse displace ment distribution is in Figure 5 30 Fully si mply supported The bending moment singularity also occurs in this case at small thickness. This singularity is demonstrated in Figure 5 31 The relevant data of this case are given in Table 5 18 and Table 5 19 And the convergence plot for the total strain energy is in Figure 5 34 and Figur e 5 35 5.4.6 Flange P late In addition to the above conventional problems, we also present the analysis result using IBFE M for an arbitrary shaped plate It also aims to test the performance of these elements for arbitrary shapes. The plate is clamped at a ll four corner holes. The geometry of this plate is provided in Figure 5 36 and a typical background mesh is presented in Figure 5 37 The result data tabulated in Table 5 20 is for the fo ur inner holes clamped scenario The convergence plot for the total st rain energy and transverse displacement plot are in Figure 5 39 Figure 5 40 and Figure 5 38 respectively. 5.5 Concluding Remarks Implicit Boundary Finite Elem ent method was employed in Chapter 5 to develop a family of Mindlin plate elements. Discrete shear field collocation method is used to make the count condition satisfied. Implicit Boundary Method is adopted to impose the essential boundary condition. Several benchmark problems and one arbitrary shape problem are tested in order to fully a ssess the performance of these elements. From the obtained results and comparison with the exact solution, analytical solution and that from S4R and S8R5 of commercial package, these developed elements are valid and robust. In conclusion, employing the Imp licit Boundary Finite Element Method can solve PAGE 82 82 plate problem for both thin and thick cases using a background mesh and equations of the boundary to avoid generating a conforming mesh Table 5 1. Location of three interpolation variables and the associated count conditions for patch test Element Clamped edges Relaxed edges One element patch Four element patch One element patch Four element patch (F) (F) : deflection interpolation node; : rotations interpolation node; : Shear collocation point; : Shear collocation point. PAGE 83 83 Table 5 2. Cantile ver plate (Shear force applied at free end) Mesh density a/t = 10 tip displacement ( ) a/t = 100 tip displacement ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10X1 3.558 3.604 3.583 3.583 3.581 3.583 3.530 3.575 3.536 3.539 3.554 3.538 20X2 3.575 3.588 3.584 3.584 3.585 3.586 3.530 3.552 3.539 3.539 3.542 3.537 50X5 3.582 3.585 3.585 3.585 3.584 3.587 3.536 3.542 3.539 3.539 3.539 3.538 100X10 3.584 3.585 3.586 3.585 3.585 3.588 3.538 3.540 3.539 3.539 3.539 3.539 Analytical 3.571 3.571 Table 5 3. Cantilever plate (Bending moment applied at free end) Mesh density a/t = 10 tip displacement ( ) a/t = 100 tip displacement ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10X1 5.324 5.342 5.339 5.342 5.348 5.342 5.324 5.342 5.321 5.326 5.348 5.325 20X2 5.337 5.342 5.342 5.343 5.347 5.345 5.320 5.334 5.325 5.325 5.332 5.323 50X5 5.341 5.343 5.343 5.344 5.344 5.344 5.323 5.327 5.326 5.326 5.326 5.325 100X10 5.343 5.344 5.344 5.344 5.344 5.344 5.325 5.326 5.326 5.326 5.326 5.326 Analytical 5.357 5.357 Table 5 4. Uniformly loaded, clamped square plate [ a/t = 10 ] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 1.457 1.642 1.482 1.467 1.470 1.466 2.335 2.539 2.368 2.298 2.278 2.375 20x20 1.464 1.515 1.470 1.467 1.468 1.467 2.324 2.375 2.332 2.315 2.310 2.333 50x50 1.467 1.475 1.468 1.467 1.467 1.467 2.321 2.329 2.322 2.319 2.318 2.322 100x100 1.467 1.469 1.467 1.467 1.467 1.467 2.320 2.322 2.320 2.320 2.320 2.321 150x150 1.467 1.468 1.467 1.467 1.467 1.467 2.320 2.321 2.320 2.320 2.320 2.320 Exact 1.467 2.320 Table 5 5. Uniformly loaded, clamped square plate [ a/t = 100 ] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 1.227 1.417 1.251 1.236 1.239 1.236 2.316 2.516 2.338 2.269 2.253 2.346 20x20 1.234 1.286 1.240 1.236 1.237 1.236 2.297 2.349 2.302 2.286 2.281 2.304 50x50 1.236 1.245 1.237 1.236 1.236 1.236 2.292 2.300 2.293 2.290 2.289 2.293 100x100 1.236 1.238 1.236 1.236 1.236 1.236 2.291 2.293 2.291 2.291 2.291 2.291 150x150 1.236 1.237 1.236 1.236 1.236 1.236 2.291 2.292 2.291 2.291 2.291 2.291 Exact 1.236 2.290 PAGE 84 84 Table 5 6. Uniformly loaded, simply supported square plate [a/t = 10] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 4.395 4.602 4.499 4.501 4.472 4.476 5.012 5.204 5.117 5.075 4.968 5.127 20x20 4.465 4.518 4.503 4.501 4.493 4.498 5.066 5.115 5.102 5.091 5.063 5.106 50x50 4.495 4.504 4.502 4.501 4.500 4.501 5.090 5.098 5.097 5.095 5.090 5.098 100x100 4.500 4.502 4.502 4.501 4.501 4.502 5.094 5.096 5.096 5.096 5.094 5.096 150x150 4.501 4.502 4.502 4.501 4.501 4.502 5.095 5.096 5.096 5.096 5.095 5.096 Exact 4.503 5.096 Table 5 7. Uniformly loaded, simply supported square plate [ a/t = 100 ] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 3.954 4.146 3.982 3.978 3.977 3.969 4.793 4.981 4.832 4.781 4.706 4.849 20x20 3.967 4.016 3.982 3.988 3.978 3.975 4.795 4.843 4.813 4.806 4.779 4.813 50x50 3.979 3.987 3.992 3.996 3.988 3.987 4.803 4.811 4.816 4.818 4.808 4.813 100x100 3.988 3.990 3.996 3.997 3.994 3.995 4.811 4.813 4.819 4.819 4.816 4.818 150x150 3.992 3.993 3.997 3.997 3.995 3.996 4.815 4.816 4.819 4.819 4.818 4.819 Exact 4.004 4.826 Table 5 8. Uniformly loaded, clamped circular plate [ a/t = 10 ] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 0.737 1.194 1.091 1.109 1.110 1.128 1.573 2.199 2.048 1.992 1.970 2.058 20x20 0.970 1.143 1.111 1.113 1.119 1.126 1.860 2.065 2.029 2.012 2.000 2.029 50x50 1.099 1.119 1.112 1.113 1.126 1.129 2.003 2.025 2.020 2.017 2.029 2.035 100x100 1.108 1.115 1.113 1.113 1.126 1.125 2.013 2.020 2.019 2.018 2.030 2.031 150x150 1.111 1.114 1.113 1.113 1.126 1.124 2.016 2.019 2.018 2.018 2.031 2.031 Analy tical 1.039 2.031 Table 5 9. Uniformly loaded, clamped circular plate [ a/t = 100 ] Mesh density Central displacement ( ) Central bending moment Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 4.745 9.100 9.017 9.313 9.374 9.502 1.348 2.075 2.033 1.985 1.970 2.043 20x20 7.027 9.535 9.371 9.399 9.471 9.504 1.731 2.049 2.026 2.011 2.002 2.028 50x50 8.895 9.459 9.409 9.410 9.535 9.537 1.959 2.025 2.019 2.017 2.029 2.032 100x100 9.253 9.428 9.413 9.413 9.538 9.539 2.000 2.020 2.018 2.018 2.030 2.032 150x150 9.360 9.419 9.413 9.413 9.538 9.539 2.012 2.019 2.018 2.018 2.031 2.031 Analy tical 9.522 2.031 PAGE 85 85 Table 5 10. Uniformly loaded, simply supported circular plate [ a/t = 10 ] Mesh density Central displacement ( ) Central bending moment Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 1.076 3.127 3.046 3.742 3.978 4.053 2.008 4.257 4.180 4.825 5.049 5.182 20x20 1.788 3.535 3.848 3.950 4.029 4.050 2.788 4.536 4.972 5.047 5.114 5.153 50x50 3.497 3.922 3.959 3.982 4.053 4.059 4.591 5.019 5.075 5.096 5.152 5.159 100x100 3.805 3.986 3.989 3.991 4.055 4.055 4.914 5.098 5.106 5.107 5.155 5.157 150x150 3.920 3.986 3.989 4.028 4.056 4.054 5.034 5.102 5.106 5.107 5.156 5.156 Analy tical 4.103 5.156 Table 5 11. Uniformly loaded, simply supported circular plate [ a/t = 100 ] Mesh density Central displacement ( ) Central bending moment Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 0.737 1.557 1.114 1.532 3.805 3.875 2.002 4.202 4.143 4.795 5.049 5.168 20x20 0.966 1.683 1.706 2.393 3.857 3.875 2.767 4.506 4.962 5.034 5.115 5.153 50x50 1.637 2.215 2.517 3.073 3.881 3.883 4.481 5.013 5.073 5.092 5.152 5.157 100x100 3.120 3.211 3.417 3.743 3.883 3.884 4.912 5.098 5.102 5.107 5.155 5.157 150x150 3.406 3.690 3.708 3.743 3.883 3.884 5.033 5.100 5.106 5.107 5.156 5.156 Analy tical 3.900 5.156 Table 5 12. Uniformly loaded, clamped 30 degree skew plate [ a/t = 10 ] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 0.806 2.576 1.717 1.726 1.642 1.767 2.545 7.479 5.183 4.635 4.684 5.087 20x20 1.320 2.106 1.729 1.729 1.695 1.735 3.711 5.767 4.999 4.802 4.626 4.928 50x50 1.662 1.812 1.728 1.728 1.724 1.730 4.624 5.003 4.894 4.859 4.861 4.868 100x100 1.723 1.755 1.728 1.728 1.727 1.728 4.847 4.895 4.877 4.868 4.861 4.880 150x150 1.726 1.742 1.728 1.728 1.728 1.726 4.859 4.878 4.873 4.869 4.868 4.885 Table 5 13. Uniformly loaded, clamped 30 degree skew plate [ a/t = 100 ] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 0.143 1.480 0.971 1.055 0.957 1.062 1.107 6.430 4.685 4.578 4.801 4.880 20x20 0.664 1.362 1.060 1.062 1.027 1.059 3.628 5.731 4.903 4.708 4.554 4.782 50x50 0.967 1.126 1.063 1.062 1.058 1.062 4.490 4.945 4.802 4.764 4.755 4.780 100x100 1.047 1.079 1.062 1.062 1.062 1.062 4.743 4.818 4.782 4.772 4.763 4.776 150x150 1.057 1.070 1.062 1.062 1.062 1.062 4.760 4.794 4.778 4.773 4.772 4.776 Exact 1.064 PAGE 86 86 Table 5 14. Unifo rmly loaded, simply supported 30 degree skew plate [a/t = 10] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 1.217 6.258 4.759 5.004 4.796 5.081 0.366 1.554 1.225 1.217 1.148 1.271 20x20 2.556 5.456 5.001 5.043 4.945 5.050 0.651 1.337 1.267 1.235 1.225 1.269 50x50 4.429 5.140 5.043 5.047 5.031 5.053 1.075 1.277 1.265 1.258 1.259 1.263 100x100 5.001 5.080 5.046 5.047 5.044 5.053 1.247 1.267 1.263 1.262 1.261 1.265 150x150 5.024 5.066 5.047 5.047 5.046 5.049 1.255 1.265 1.263 1.262 1.262 1.264 Exact 5.047 Table 5 15. Uniformly loaded, simply supported 30 degree skew plate [ a/t = 100 ] Mesh density Central displacement ( ) Q4 Q8 Q9 Q16 S4R S8R5 10x10 0.255 4.584 2.707 3.451 3.897 4.069 20x20 0.867 4.036 3.334 3.802 4.013 4.083 50x50 1.671 4.032 3.909 4.081 4.105 4.118 100x100 3.111 3.909 4.080 4.132 4.130 4.135 150x150 3.410 4.122 4.122 4.138 4.135 4.138 Exact 4.070 Table 5 16. Uniformly loaded, clamped 60 degree skew plate [ a/t = 10 ] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 5.358 11.31 9.446 9.355 9.315 1.033 1.257 2.031 1.751 1.581 1.635 1.764 20x20 8.793 9.916 9.339 9.320 9.302 1.034 1.624 1.766 1.692 1.650 1.665 1.727 50x50 9.237 9.413 9.305 9.304 9.299 1.034 1.669 1.688 1.677 1.670 1.673 1.717 100x100 9.288 9.331 9.300 9.300 9.298 1.034 1.674 1.677 1.675 1.674 1.674 1.715 150x150 9.294 9.315 9.299 9.299 9.298 1.034 1.674 1.675 1.674 1.674 1.674 1.720 Table 5 17. Uniformly loaded, clamped 60 degree skew plate [ a/t = 100 ] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 2.764 8.688 7.539 7.525 7.539 7.527 1.283 1.897 1.728 1.989 1.623 1.704 20x20 5.643 8.072 7.546 7.519 7.521 7.528 1.526 1.750 1.672 1.644 1.645 1.666 50x50 6.872 7.622 7.523 7.517 7.518 7.529 1.622 1.670 1.656 1.652 1.652 1.656 100x100 7.357 7.545 7.519 7.517 7.517 7.529 1.647 1.657 1.654 1.653 1.653 1.654 150x150 7.456 7.530 7.518 7.517 7.517 7.529 1.651 1.655 1.653 1.653 1.653 1.654 Analy tical 7.517 PAGE 87 87 Table 5 18. Uniformly loaded, simply supported 60 degree skew plate [ a/t = 10 ] Mesh density Central displacement ( ) Central bending moment ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 0.809 3.103 2.898 2.923 2.881 2.999 1.721 4.082 3.895 3.798 3.743 3.937 20x20 2.349 2.975 2.912 2.914 2.900 3.009 3.201 3.932 3.875 3.845 3.829 3.907 50x50 2.837 2.925 2.910 2.910 2.907 3.010 3.768 3.875 3.863 3.858 3.855 3.896 100x100 2.904 2.916 2.909 2.909 2.908 3.010 3.854 3.867 3.861 3.860 3.859 3.895 150x150 2.907 2.914 2.909 2.909 2.909 3.010 3.858 3.865 3.860 3.860 3.860 3.894 Exact 2.909 Table 5 19. Uniformly loaded, simply supported 60 degree skew plate [ a/t = 100 ] Mesh density Central displacement ( ) Q4 Q8 Q9 Q16 S4R S8R5 10x10 0.352 2.722 2.399 2.489 2.507 2.510 20x20 0.868 2.596 2.487 2.520 2.511 2.517 50x50 1.951 2.551 2.527 2.535 2.524 2.530 100x100 2.430 2.537 2.535 2.537 2.532 2.536 150x150 2.483 2.539 2.536 2.537 2.534 2.537 Exact 2.554 Table 5 20. Uniformly loaded, arbitrary shape plate Mesh density a/t = 10 Central displacement ( ) a/t = 100 Central displacement ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 20x20 5.606 7.350 6.927 7.310 7.540 7.379 3.235 4.539 4.234 4.551 4.770 4.698 50x50 5.754 7.440 7.275 7.458 7.423 7.356 3.635 4.668 4.610 4.604 4.739 4.698 100x100 7.260 7.502 7.423 7.475 7.368 7.356 4.542 4.749 4.671 4.701 4.707 4.697 150x150 7.070 7.467 7.430 7.475 7.356 7.346 4.505 4.754 4.738 4.704 4.699 4.697 Figure 5 1. A typical background mesh using 20x2 Q4 element PAGE 88 88 Figure 5 2. Distribution of transverse displacement af ter deformation for 100x10 Q4 element (L1 /t = 100) Figure 5 3. Convergence of total strain energy for cantilever when shear applied (L1 /t = 10) Figure 5 4. Convergence of total strain energy for c antilever when shear applied (L1 /t = 100) PAGE 89 89 Figure 5 5. Distribution of transverse displacement af ter deformation for 100x10 Q4 element (L1 /t = 100) Figure 5 6. Convergence of total strain energy for cantilever when ben ding moment applied (L1 /t = 10) Figure 5 7. Convergence of total strain energy for cantileve r when bending moment applied (L1 /t = 100) PAGE 90 90 Figure 5 8. A typical background mesh using 10x10 Q9 element Figure 5 9. A typical background mesh using 10x10 Q9 element Figure 5 10. Distribution of transverse displacement after deformation for 150x150 Q9 element (t = 0.1) PAGE 91 91 Figure 5 11. Distribution of transverse displacement after deformation for 225x225 Q9 element (t = 0.1) Figure 5 12. Convergence of total strain energy for cl amped square plate (a/t = 10) Figure 5 13. Convergence of total strain energy for clamped square plate (a/t = 100) PAGE 92 92 Figure 5 14. Convergence of total strain energy for simply supported square plate (a/t = 10) Figure 5 15. Convergence of total strain e nergy for simply supported square plate (a/t = 100) Figure 5 16. A typical background mesh using 10x10 Q4 element PAGE 93 93 Figure 5 17. Distribution of transverse displacement after deformation for 150x150 Q4 element (t = 1) Figure 5 18. Convergence of total strain energy for clamped circular plate (a/t = 10) Figure 5 19. Convergence of total strain energy for clamped circular plate (a/t = 100) PAGE 94 94 Figure 5 20. Convergence of total strain energy for simply supported circular plate (a/t = 10) Figure 5 21. Con vergence of total strain energy for simply supported circular plate (a/t = 100) Figure 5 22. A typical background mesh using 10x10 Q8 element PAGE 95 95 Figure 5 23. Distribution of transverse displacement after deformation for 200x200 Q8 element (t = 0.1) Figure 5 24. Distribution of bending moment using 200x200 Q4 element (t = 0.1) Figure 5 25. Convergence of to tal strain energy for clamped 30 degree skew plate (a/t = 10) PAGE 96 96 Figure 5 26. Convergence of to tal strain energy for clamped 30 degree skew plate (a/t = 100) Figure 5 27. Convergence of total strain energy for simply su pported 30 degree skew plate (a/t = 10) Figure 5 28. Convergence of total strain energy for simply supported 30 degree skew plate (a/t = 100) PAGE 97 97 Figure 5 29. A typical background mesh using 10x10 Q16 element Figure 5 30. Distribution of transverse displacement after deformation for 200x200 Q9 element (t = 1) Figure 5 31. D istribution of bending moment using 200x200 Q4 element (t = 0.1) PAGE 98 98 Figure 5 32. Convergence of total strain energy for clamped 60 degree s kew plate (a/t = 10) Figure 5 33. Convergence of total strain energy for clamped 60 degree skew plate (a/t = 100) Figure 5 34. Convergence of total strain energy for simply supported 60 degree skew plate (a/t = 10) PAGE 99 99 Figure 5 35. Convergence of total strain energy for simply supported 60 degree skew plate (a/t = 100) Figure 5 36. Geometry of flange plate Figure 5 37. A typical background mesh using 20x20 Q9 element PAGE 100 100 Figure 5 38. Distribution of transverse displace ment after deformation for 150x150 Q9 element (t = 0. 1) Figure 5 39. Convergence of total strain energy for arbitrary shape plate (a/t = 10) Figure 5 40. Convergence of total strain energy for arbitrary shape plate (a/t = 100) PAGE 101 101 CHPATER 6 MIXED FORMULATION FO R 2D MINDLIN SHELL In plate elements we are mainly concern ed about bending and transverse loading. Plate s subjected to both transverse loading and in plane s tretching is a special case of s hells Chapter 6 is intended to model this spe cial case, and hence prepare for future work on curved 3D s hell like structures. 6 .1 Governing E quations The derivation of the g overning equations for Mindlin plate has already been done i n Chapter 4. In Chapter 4, we assume that there is no in pl ane stret ching, therefore, terms associated with stretching are omitted. In Chapter 6 we will consider in plane stretching also The derivation is the same as in Chapter 4, except that two mor e equations are added into the g overning equations. T he governing equati ons are (6 1) (6 2) (6 3) (6 4) (6 5) Writing in matrix form, the governing equation system can be expressed as PAGE 102 102 (6 6 ) (6 7 ) (6 8 ) (6 9 ) (6 10 ) (6 11 ) For clarity, notations used in above equations are rewritten here. (6 12 ) (6 13 ) (6 14 ) PAGE 103 103 Equation s 6 6 to 6 11 constitute the governing equation system for Chapter 6 Substituting notat ion 6 12 to 6 14 into E quation s 6 6 and 6 11 respectively, we can eliminate the in plane forces and moment and finally obtain the four governing equations as: (6 15 ) So far, we have obtained the governing equation for plates including in plane stretching. In next section, we will formulate the mixed form for this plate. 6 .2 Mixed F ormulation We have already derived the governing equations for Mindlin plate including in plane stretching. Next, we will focus on the derivation of the mixed form. As we have mentioned in the last section, except the fi rst governing equation, the remaining three ar e the same as those in Chapter 5 And the weak mixed form for these three equations are still the same since there is no coupling between in plane stretching and transverse bending in our formulation. For si mplicity, we will only derive the weak form for the first equation. For the other three, one can refer to Chapter 5 for details. Using integration by part the weak form of 6 6 to 6 11 can be deduc ed as PAGE 104 104 (6 16 ) is the prescribed traction on the boundary The final weak system of the combined Mindlin plate is (6 17 ) Writing in matrix form, the mixed weak form can be expressed as (6 1 8 ) Where (6 19 ) Strains and Stresses A fter solving the equations for and at each node, we need to calculate the associated strains and stresses. (6 20 ) PAGE 105 105 Once the membrane strain and curvature are obtained, t he in plane strains on any layer can be immediately calculated. A t the height equal to the membrane strains on the top surface of the plate i s (6 21 ) Using constitutive relation ship we can obtain the correspond ing stresses on the top surface of the plate as (6 22 ) where is the stiffness matrix. The average transverse shear stresses can be computed as (6 23 ) We also can calculate the in plane forces, bending moments and shear forces as (6 24 ) PAGE 106 106 6 .3 Numerical E xamples and R esults The way to impose the EBC is almost the same as what have been done in Chapter 5. We will not redo the procedures for combined Mindlin plates. One can refer to Chapter 5 for details. Examples employed in Chapter 6 are the 60 degree skew plate and square plate that have been used in Chapter 5. The geometry and elastic constants are exactly the same as in Chapter 5. 6.3.1 60 degree Skew P late The plate is clamped only at the left edge. A transverse pressure and in plane stretching force in the x direction, both with value of 100, are applied to this plate. See Figure 6 4 Since there is no analytical or exact solution for this loaded plate, we only compare our results with S4R and S8R5 element from Abaqus. Both thin case and thick case are studied in this section. The testing date for transverse deflection is in Table 6 1 and Table 6 2 The convergence plot of total strain energy is in Figure 6 1 and Figure 6 2 A plot of the transverse deflection for thin skew plate is also given as Figure 6 3 6.3.2 Square P late The square plate is clamped at the left edge. See Figure 6 8 A transverse pressure and in plane stretching force in the x direction, both with value of 100, are applied on this plate. Test data for displacement are tabulated in Tabl e 6 3 and Table 6 4 The convergence plot of tota l strain energy is in Figure 6 5 and Figure 6 6 The deformed transv erse deflection is in Figure 6 7 From above testing data and convergence plots of two examples, we can see that the new element s performance the same as S4R element, but better than S8R5 element. PAGE 107 107 6 .4 Concluding R emarks We have successfully extended Mindlin plate theory to incl ude in plane stretching in Chap ter 6 Since we study the elastic static problem, there is no coupling between in plane stretching and bending in our formulation. The mixed weak form is a li near combination of plane stress and Mindlin plate theory. Examples validate th e performance of th ese IBFEM element s that we have developed The content of Chapter 6 can be viewed as a special case of shells problem. Future work is to extend this into a general shell analysis. Table 6 1. Transverse deflection of 60 degree skew plate with one edge cla mped Mesh density a/t = 10 Maximum d eflection ( ) a/t = 100 Maximum deflection ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 0.970 1.233 1.218 1.219 1.218 1.230 0.585 1.198 1.184 1.184 1.185 1.182 20x20 1.176 1.222 1.218 1.219 1.218 1.232 0.934 1.188 1.184 1.185 1.183 1.183 50x50 1.209 1.219 71.218 1.219 1.218 1.232 1.068 1.185 1.185 1.185 1.184 1.184 100x100 1.216 1.219 1.219 1.219 1.219 1.232 1.160 1.185 1.185 1.185 1.184 1.185 150x150 1.217 1.219 1.219 1.219 1.219 1.232 1.174 1.185 1.185 1.185 1.184 1.185 Table 6 2. In plane displacement of 60 degree skew plate with one edge clamped Mesh density a/t = 10 Maximum d isplacement ( ) a/t = 100 Maximum displacement ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x10 1.034 1.034 1.057 1.033 1.033 1.033 1.034 1.034 1.057 1.033 1.034 1.033 20x20 1.035 1.033 1.045 1.033 1.033 1.033 1.035 1.033 1.045 1.033 1.033 1.033 50x50 1.034 1.033 1.038 1.033 1.033 1.033 1.034 1.034 1.038 1.033 1.033 1.033 100x100 1.034 1.033 1.036 1.033 1.033 1.033 1.034 1.034 1.036 1.033 1.033 1.033 150x150 1.034 1.033 1.033 1.033 1.033 1.033 1.034 1.034 1.033 1.033 1.033 1.033 Table 6 3. Transverse deflection of square plate with one edge clamped Mesh density a/t = 10 Maximum deflection ( ) a/t = 100 Maximum deflection ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x1 0 1.279 1.295 1.281 1.281 1.281 1.281 1.258 1.276 1.260 1.259 1.259 1.259 20x2 0 1.280 1.285 1.281 1.281 1.281 1.281 1.258 1.264 1.259 1.259 1.259 1.259 50x5 0 1.281 1.282 1.281 1.281 1.281 1.281 1.259 1.260 1.259 1.259 1.259 1.259 100x10 0 1.281 1.281 1.281 1.281 1.281 1.281 1.259 1.259 1.259 1.259 1.259 1.259 PAGE 108 108 Table 6 4. In plane displacement of square plate with one edge clamped Mesh density a/t = 10 Maximum displacement ( ) a/t = 100 Maximum displacement ( ) Q4 Q8 Q9 Q16 S4R S8R5 Q4 Q8 Q9 Q16 S4R S8R5 10x1 0 8.855 8.861 8.862 8.863 8.854 8.863 8.855 8.861 8.862 8.863 8.893 8.863 20x2 0 8.861 8.863 8.863 8.864 8.862 8.864 8.861 8.863 8.863 8.864 8.862 8.864 50x5 0 8.863 8.864 8.864 8.864 8.863 8.864 8.863 8.864 8.864 8.864 8.863 8.864 100x10 0 8.864 8.864 8.864 8.864 8.864 8.864 8.864 8.864 8 .864 8.864 8.864 8.864 Figure 6 1. Convergence of total strain energy for 60 degree skew plate (a/t = 10) Figure 6 2. Convergence of total strain energy for 60 degree skew plate (a/t = 100) PAGE 109 109 Figure 6 3. Distribution of transverse d eflection after deformation for 60 degree skew plate using 150x150 Q4 element (t = 1) Figure 6 4. Geometry of the 60 degree skew plate Figure 6 5. Convergence of total strain energy for square plate (a/t = 10) PAGE 110 110 Figure 6 6. Convergence of total strain energy for square plate (a/t = 100) Figure 6 7. Distribution of transverse d eflection after def ormation for square plate using 100x100 Q9 element (t = 0. 1) Figure 6 8. The geometry of the square plate PAGE 111 111 CHAPTER 7 CONCLUSION 7 .1 Summary In summary we have formulated and ev aluated elements that use mixed formu lation to avoid locking and the implicit boundary method to apply boundary conditions. These elements avoid volumetric locking and shear locking phenomena in near incompr essible media and thin Mindlin plate respectively These elements enables mesh independent analysis where the geometry is represented as equations and a background mesh is used to construct the trial solution. The background mesh can consist of uniform ele ments and need not conform to the shape of the analysis domain because the elements developed here use implicit boundary method to apply boundary conditions which does not need nodes on the boundary to apply boundary conditions. For volumetric locking, since this phenomenon is totally due to the fact that bulk modulus tends to field mixed formulation is adopted for the plane strain and 3D stress in Chapter 3. From the convergence plot and transverse displacement d istribution plot of 3D stress we can concl ude that these elements perform very well and th e volumetric locking is removed Before proceed ing to the formulation of Mindlin plate theory, two classic al plate theories were reviewed in Chapter 4. The difference between these two th eories is clearly stated in Chapter 4 For future evaluation purpose, the analytical or exact solutions of several benchma rk problems are provided in Chapter 4 also Mixed for mulation of Mindlin plates using discrete shear field collocation method and the imposition of the essential boundary condition using implicit boundary method are detailed in Chapter 5 From the comparison of the results from analytical solution, PAGE 112 112 exact sol ution, IBFEM element solution and Abaqus element analysis, the performance of the IBFEM elem ents can be easily assessed. It is clear that these elements are useful an d robust, and at least as good as, if not better than, the elements fr om the commercial pa ckage ABAQUS 2D M indlin shell that includes in plane stretching and transverse pressure, is also studied in this thesis. From the results, these elements also perform well. The objective of this extension of Mindlin plate is to prepare the formulation of general curved s hell elements using IBFEM. In conclusion, we have successfully implemented the mixed formulation using IBFEM so as to avoid volum etric locking and shear locking phenomena 7 .2 S cope of Future Work The three field mixed formulation we have employed in this thesis is only valid for nearly incompressible media, since there is still bulk modulus term in the modified stiffness matrix. Hence, in order to remove volumetric locking for completely incompressible media, new finite element te chniques should be used for this purpose Two field ( ) mixed formulation is a possible way according to Hughes [9] and Bathe [1] The mixed formulation using discrete shear field collocation method and IBFEM is developed in this thesis only for elastic static Mindlin plates. T he validity of this metho d is already evaluated in Chapter 5. Future work can extend this method to nonlinear and dynamic analysis of Mindlin plates. The mixed formulation of 2D Mindlin shell including both in plane st retching and bending is derived in Chapter 6. This c an be vie wed a special 2D case of s hell PAGE 113 113 formulation. Hence, g eneralizing this method for 3D shell analysis can be a prosp ective future direction PAGE 114 114 APPDENDIX A VOLUMETRIC LOCKING A ND SHEAR LOCKING Volumetric locking and shear locking are the most common numerical phenomena in the displacement base finite element formulation. We will gi ve some details in order to better under stand these two frequent phenomena in FEM. A.1 Volumetric locking Volumetric locking is very obvi ous as we already presented at the beginning o f Chapter 3 This phenomenon results from the fact that when the material is near tends to 0.5, which makes the bulk modulus tend to infinity and hence a ill c onditioned stiffness matrix in the finite element model. A.2 Shear locking Shear locking is an error that occurs in finite element analysis. It can be explained as following Let us consider a beam with thickness under combined shear and bending deformations. As we have done in Chapter 5 the total strain energy can be expressed as ( A 1 ) Where PAGE 115 115 ( A 2 ) ( A 3 ) From (A 1 ), t he portion of the strain energy associated with the bending defo rmation is propor tional to and the portion related to shear def ormations is proportional to If we reduce the thickness the value of approaches to zero much faster than and as a result, the total strain energy of the beam will largely come from shear deformation. This is not correct because for a thin beam it is always the bending deformation which provide s most part of the strain energy. So any shear strain resulting from numerical errors and due to the use of low order shape functions can lead to wrong results for very thin beam small displacement PAGE 116 11 6 APPENDIX B EQUILIBRIUM EQUATION S OF 3D ELASTOSTATIC CASE Figure B 1. Stresses notation s and direction s Comparing to these stresses, the body force is of high order, so we neglect it in above figure B 1 But in the process of deriving the equilibrium equation, we should take stresses. The equation of equilibrium of force in the x direction is (B 1) Dividing by (B 2) If the block is become smaller and small er, i.e., (B 3) By definition of derivative, we have PAGE 117 117 (B 4 ) (B 5 ) (B 6 ) Plugging B 4 to B 6 into B 2 we have the first differential equilibrium equation as (B 7 ) Using the same manner, we can derive the left two equations as (B 8 ) (B 9 ) PAGE 118 118 APPENDIX C DERIVATION OF SHEAR CORRECTION FACTOR The transverse shear correction coefficient can be derived using equivalence of shear strain energy before and after using this factor. Equivalence of Shear strain energy The first two equilibrium equations of a static case without body force are (C 1 ) (C 2 ) Integrating C 1 and C 2 we can obtain the transverse stress at any point through the thickness as (C 3 ) (C 4 ) a nd are the shear stresses at the bottom of the plate and both are equal to zero. Given the in plane strains, the in plane stresses can be obtained using the stress strain relationship as (C 5 ) PAGE 119 119 Substituting Equation C 5 into E quation s C 3 and C 4 and performing the integration, we obtain the following results (C 6 ) (C 7 ) For rectangular plate, the transverse stresses parabolically vary through the thickness Figure C 1. Distribution of transverse shear stress through the thickness In E quation s C 6 and C 7 we have following relations (C 8 ) Thus, and can be further expressed as (C 9 ) PAGE 120 120 (C 10 ) From Equation C 8, C 9 and C 10 we finally have (C 11 ) (C 12 ) Generally, the total strain energy per unit area due to transverse shears is (C 13 ) Pluggi ng Equations C 11 and C 12 into Equation C 13 we have (C 14 ) In shear deformable plate theory, we assumed at the beginning that (C 15 ) Hence, we also have the total strain energy per unit area due to transverse shears as (C 16 ) Compar ing Equation C 14 and Equation C 16 for total strain energy due to the same transverse shears, we have (C 17 ) PAGE 121 121 APPENDIX D DERIVATION OF THE JA COBIAN MATRIX IN IBF EM The Jacobian matrix for 2D is in the format as following: (D 1) We can modified above expression for of 9 node element as (D 2) In IBFEM, we have following relations (D 3) Using these relations in above expression for we can get (D 4) PAGE 122 122 Figure D 1. Global coordinate and Local coordinate in IBFEM PAGE 123 123 APPENDIX E FORMULATION OF MINDL IN PLATE ELEMENTS E.1 Element Q4D 4 Figure E 1 Collocation constraints on a 4 node Lagrange element Shape function for displacement and rotations The shape functions for rotation are the same as those for the displacement with corresponding node. These shape functions for displacement are (E 1) The shape functions for rotations are (E 2) Shape function for shear fields PAGE 124 124 The shape functions for shear fields can be derived using the same structure of Lagrange polynomials. And these shape fu nctions are (E 3) Derivation of the prescribed shear force value Figure E 2. Interpolation nodes for Q4D4 element For above Q4D 4 element, the displacement and the rotation can be approximated as and Using (E 4) we can derive the prescribed values for shear force at certain nodes using prescribed and as PAGE 125 125 (E 5) Where and are the and components of and In IBFEM, we have where and are th e length of edges of an element in and direction, respectively. Writing in the matrix form, we have where (E 6 ) (E 7 ) PAGE 126 126 (E 8 ) T he same procedure is used to get the collocated shear term and its shape function and we will not give more details in the following content about the derivation for each element And also, the shape functions for the displacement are the same as traditional ones. The shape functions for rotations are the same as those of displacement at the corresponding nodes E.2 E lement Q5D 6 Figure E 3. Collocation constraints on a 5 node Serendipity element Shape function for shear fields (E 9 ) PAGE 127 127 Figure E 4. Interpolation nodes for Q5D6 element and matrices (E 10 ) (E 11 ) (E 12 ) PAGE 128 128 E.3 Element Q8 D 8 Figure E 5. Collocation constraints on an 8 node Serendipity element Shape functions for shear fields (E 13 ) Figure E 6. Interpolation nodes for Q8D8 element and matrices PAGE 129 129 (E 14 ) (E 15 ) (E 16 ) PAGE 130 130 E.4 El ement Q9D 12 Figure E 7. Collocation constraints on a 9 node Lagrange element Shape functions for shear fields (E 17 ) (E 18 ) and matrices PAGE 131 131 Figure E 8. Interpolation nodes for Q9D12 element (E 19 ) (E 20 ) PAGE 132 132 (E 21 ) E.5 Element Q16D24 Figure E 9. Collocation cons traints on a 16 node Lagrange elemen t Shape function for shear fields PAGE 133 133 (E 22 ) (E 23 ) (E 24 ) PAGE 134 134 (E 25 ) and matrices (E 26 ) (E 27 ) PAGE 135 135 (E 28 ) Figure E 10. Interpolation nodes for Q16D24 element PAGE 136 136 LIST OF REFERENCES [1] Bathe K.J. (1996). Finite Element Procedures 2nd ed., Prentice Hall, New Jersey. [2] Bathe K.J Dvorkin, E.N. (1985). A four node plate bending element based on Mindlin/Reissner plate t heory and a mixed interpolation. International Journal of Num erical Methods in Engineering, 21(2), 367 383 [3] Belytschko T. Krongauz, Y., Organ, D., Fleming, M., Krysl, P. (1996). Meshless Methods: An Ov erview and Recent Developments. Comput er Methods in Appl ied Mechanics and Engineering 139 3 47 [4] Belytschko T., Parimi C., Moe s, N., Usui, S., Sukumar, N. (2003). Structure d extended finite element methods for solids defined by implicit surfaces. International Journal for Numerical Methods in Engineering 56(4) 609 635 [5] Clark, B.W., Anderson, D.C. (2002). Finite Element Analysis in 3D Using Penalty Boundary Method. Proceedings of Design Engineering Technical Conferences, Montreal, Canada [6] Fish, J., Belytschko, T. (2007). A first course in finite ele ments, John Wiley and Sons [7] Hinton, E., Huang, H.C. (1986). A family of quadrilateral Mindlin plate elements with substi tute shear strain fields. Computers and Structures, 23 (3), 409 431 [8] Hrabok, M.M., Hrudey, T.M. (1984). A review and catalogue of plate be nding finite elements. Computers and Structures, 19(3), 479 495 [9] Hughes, T.J.R. (2000). The Finite Element Method: Lin ear Static and Dynamic Finite Element Analysis, 2nd e dition. Dover Editions, New Jersey [10] Kantorovich, L.W., Krylov, W.I. (1956). Nherungsmethoden der Hheren Analysis. VEB Deutscher Ve rlag der Wissenschaften, Berlin [11] Kumar, A.V., Burla, R., Padmanabhan, S., Gu, L. (2008). Finite element an alysis using nonconforming mesh. Journal of Computing a nd Information Science in Engine ering 8, 031005 1 [12] Kumar A.V., Lee, J. (2003). Step function representation of solid models and application to mesh free engineering analysi s, Journal of Mechanical Design. Transactions of the ASME 128(1) 46 56 [13] Kumar, A.V. Padmanabhan, S. Burla, R. (2008). Implicit boundary method for finite element analysis using n on conforming mesh or mesh International Jo urnal for Numerical Methods in Eng inee r ing, 74 (9), 1421 1447 PAGE 137 137 [14] Liew, K.M., Han, J. B. (1997). Bending analysis of simply support ed shear deformable skew plates. Engineering Mechanics, 123( 3), 214 221 [15] Rvachev, V.L., Shieko, T.I. (1995). R functions in Boun dary Value Problems in Mechanics. Applied Mechanics Review 48 151 188 [16] Sengupta, D. (1995). Performance study of a simple finite element in the analysis of skew rhombic plates. Computers and Structures 54 (6) 1173 1182 [17] Taylor, R.L., Govindjee, S. (2004. ) Solution of clamped rectangular problems. Communications in Numer ical Methods in Engineering. 20, 757 765 [18] Timoshenko, S., Woinowsky Kreiger, S. (1987). Theory of Plates and Shells, second edition. McGraw Hill, New York [19] Wang, C.M., Reddy, J.N., Le e, K.H. (2000). Shear Deformable Beams and Plates: Relationships with Classical S olutions. Elsevier, Oxford [20] Yang, T.Y., Saigal, S., Masud, A., Kapania, R.K. (2000). A survey of recent shell finite elements. International Journal for Numerical Methods i n Engineering, 47, 101 127 [21] Zienkiewicz, O.C., Taylor, R.L., Zhu, J.Z. (2005). The Finite Element Method for Solid and Structural Mechanics, sixth ed ition Butterworth Heinemann [22] Zienkiewicz, O.C., Xu, Z., Zeng, L.F., Sam uelson, Wilbert, N.E. (1993). Linked interpolation for Reissner Mindlin plate elements: Part I A simple quadrilateral. Int ernational J ournal for numer ical Methods in En gineering, 36, 3043 3056 PAGE 138 138 BI OGRAPHICAL SKETC H Hailong Chen was born and brought up in Jiangxi, China. He graduated wi th a Bachelor of Science in engineering degree in mechanical engineering from Shanghai Normal University, Shanghai, China in June 2010 with honors After graduation, he directly enrolled in the master program in Mechanical and Aerospace Engineeri ng at th e University of Florida in f all 2010 with Achievement Award for New Engineering Graduate S tudent s In 2011, he and Xiao Zhou got married in China. His area of interest includes application of FEM in plates and shells, numerical methods and uncertainty quantification. 