Citation
Contact Modeling Using Implicit Boundary Finite Element Approach

Material Information

Title:
Contact Modeling Using Implicit Boundary Finite Element Approach
Creator:
Tupsakhare, Anand S
Place of Publication:
[Gainesville, Fla.]
Publisher:
University of Florida
Publication Date:
Language:
english
Physical Description:
1 online resource (106 p.)

Thesis/Dissertation Information

Degree:
Master's ( M.S.)
Degree Grantor:
University of Florida
Degree Disciplines:
Mechanical Engineering
Mechanical and Aerospace Engineering
Committee Chair:
Kumar, Ashok V
Committee Members:
Arakere, Nagaraj K
Kim, Nam Ho
Graduation Date:
5/5/2012

Subjects

Subjects / Keywords:
Boundary conditions ( jstor )
Contact loads ( jstor )
Cylinders ( jstor )
Flat plates ( jstor )
Geometry ( jstor )
Modeling ( jstor )
Shear stress ( jstor )
Step functions ( jstor )
Stress distribution ( jstor )
Symmetry ( jstor )
Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF
contact
Genre:
Electronic Thesis or Dissertation
born-digital ( sobekcm )
Mechanical Engineering thesis, M.S.

Notes

Abstract:
Several structured mesh based finite element analysis techniques have been developed recently to avoid mesh generation difficulties. Structured meshes are composed of regular shaped uniform elements which are easy to generate automatically as compared to mesh generation in traditional FE software. In the Implicit Boundary Finite Element Method, IBFEM, used in this work, solution structures for the trial and test functions are constructed using approximate step functions of the boundaries, such that allows boundary conditions to be imposed even when there are no nodes on the boundary. In this study, Implicit Boundary Finite Element Method is further extended to obtain solution for tied contact problems where the parts in contact stay in contact without sliding or separating during the deformation. The solution structure presented previously by Burla et al. (13) for treatment of material discontinuity in a structured mesh framework was improved to handle point contact between two bodies. In this approach, contacting bodies are modeled as an assembly of parts and separate mesh is constructed for each part. At the contact, elements of the mesh belonging to the contacting parts overlap. A solution structure is constructed such that the displacements are continuous at the contact. At the beginning of contact, non-conforming bodies with sufficiently large relative curvatures touch each other at a point or along a line. With increasing load, they deform and contact increases over a finite area which is usually extremely small in scale as compared to the overall dimensions of the contacting bodies. Contact stresses in this contact patch are significantly higher than the rest of the regions giving rise to very steep stress gradients at the boundary of the contact patch. In the current work, an attempt is made to capture these stresses as accurately as possible using sufficiently dense global mesh so as to model the contact patch width to near accurate dimensions as calculated from the analytical solutions. Ideally, for these analyses to be more efficient, localized refinement of the involved meshes is necessary. Here we focus only on the case in which identical linear elastic material properties exist for both the contacting bodies. Several tied contact problems including classical point and line contact problems are used to study the effectiveness of our approach used here in modeling contact using implicit boundary finite element method. ( en )
General Note:
In the series University of Florida Digital Collections.
General Note:
Includes vita.
Bibliography:
Includes bibliographical references.
Source of Description:
Description based on online resource; title from PDF title page.
Source of Description:
This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Thesis:
Thesis (M.S.)--University of Florida, 2012.
Local:
Adviser: Kumar, Ashok V.
Electronic Access:
RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2013-05-31
Statement of Responsibility:
by Anand S Tupsakhare.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright Tupsakhare, Anand S. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:
5/31/2013
Resource Identifier:
864716243 ( OCLC )
Classification:
LD1780 2012 ( lcc )

Downloads

This item has the following downloads:


Full Text

PAGE 1

1 CONTACT MODELING USING IMPLICIT BOUNDARY FINITE ELEMENT APPROACH By ANAND SHRIKANT TUPSAKHARE A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGRE E OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2012

PAGE 2

2 2012 Anand Shrikant Tupsakhare

PAGE 3

3 To my parents Sudha and Shrikant Tupsakhare my brother Swanand my fiance Ketki and my advisor Dr. Ashok Kumar

PAGE 4

4 ACKNOWLEDGMENTS I owe my deepest gratitud e to my ad visor, Dr Ashok Kumar whose guidance encouragement and support develop an d complete this thesis in a sensible and fruitful manner with deeper understanding of the subject s using the numerous insights he provi ded during every stage of this research. Without his assistance, this work would not have been possible. The great learning experience gained from him will surely help me build my career and every aspect of my future life well. I am also grateful to Dr. Na garaj Arakere and Dr. Nam Ho Kim for their willingness to provide me with help useful suggestions and inputs which really helped me complete this thesis successfully. Moreover, i t is an honor for me to thank all my family members who made this work possib le and helped in every aspect possible. And finally, it is a pleasure to appreciate the support of my lab mates and friends especially Jinsang, Vijaykrishna, Hailong, Sriram, Yongmin and Jinuk, for many fruitful discussions we had here at the Design and R apid Prototyping Laboratory of University of Florida

PAGE 5

5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ 4 LIST OF TABLES ................................ ................................ ................................ ........... 7 LIST OF FIGURES ................................ ................................ ................................ ........ 8 ABSTRACT ................................ ................................ ................................ .................. 11 CHAPTERS 1 INTRODUCTION ................................ ................................ ................................ ... 13 1.1 Overview ................................ ................................ ................................ ......... 13 1.2 Goals and Objectives ................................ ................................ ....................... 15 1.3 Outline ................................ ................................ ................................ ............. 16 2 OVERVIEW OF CLASSICAL CONTACT PROBLEMS ................................ .......... 17 2.1 Theory of Contact ................................ ................................ ............................ 17 2.1.1 Normal Contact of Elastic Solids ................................ ............................ 17 2.1.2 Fundament als of Hertz Theory ................................ ............................... 19 2.1.3 Growth of Contact Area, Deformation and Stresses with Increasing Applied Load ................................ ................................ ................................ 21 2.2 Plane Problems: Formulation and Analytical Solutions ................................ .... 23 2.3 Axi symmetric Problems Formulation ................................ .............................. 30 2.4 Traditional Methods for Numerical Contact Modeling ................................ ....... 35 2.4.1 Lagrange Multiplier Form ................................ ................................ ........ 37 2.4.2 Penalty Form ................................ ................................ .......................... 38 3 MESH INDEPENDENT FINITE ELEMENT ANALYSIS ................................ ......... 42 3.1 The Finite Element Method ................................ ................................ .............. 42 3.2 Mesh Independent an d Meshless Methods ................................ ...................... 42 3.3 Implicit Boundary Finite Element Method ................................ ......................... 44 3.4 Development of Modified Weak Form for Linear Elasticity in IBFEM ................ 47 4 IMPLEMENTATION OF POINT CONTACT USING IBFEM ................................ ... 49 4.1 Solut ion Structure for Contact ................................ ................................ .......... 49 4.2 Modified Weak Form for Contact Implementation ................................ ............ 50 4.3 Numerical Implementations for Contact ................................ ........................... 52 5 RESULTS ................................ ................................ ................................ .............. 58

PAGE 6

6 5.1 Example 1: Cantilever Loaded Through Contact with a Half Cylinder .............. 60 5.1.1 Stresses and Displacements ................................ ................................ .. 62 5.1.2 Converg ence Studies ................................ ................................ ............. 63 5.2 Example 2: Blocks in Contact ................................ ................................ .......... 64 5.2.1 Stresses and Displacements ................................ ................................ .. 64 5.2.2 Convergence Stu dies ................................ ................................ ............. 66 5.3 Example 3: Cylinder Contacting a Flat ................................ ............................. 66 5.3.1 Stresses and Displacements ................................ ................................ .. 67 5.3.2 Convergence Studies ................................ ................................ ............. 69 5.4 Example 4: Contact of Cylinders ................................ ................................ ...... 69 5.4.1 Stresses an d Displacements ................................ ................................ .. 70 5.4.2 Convergence Studies ................................ ................................ ............. 71 5.5 Example 5: Sphere Contacting a Flat ................................ .............................. 72 5.5.1 Stresses and Displacements ................................ ................................ .. 72 5.5.2 Convergence Stu dies ................................ ................................ ............. 74 5.6 Example 6: Contact of Spheres ................................ ................................ ....... 75 5.6.1 Stresses and Displacements ................................ ................................ .. 75 5.6.2 Convergence Studies ................................ ................................ ............. 76 6 CONCLUSION AND FUTURE SCOPE ................................ ............................... 101 6.1 Summary ................................ ................................ ................................ ....... 101 6.2 Scope of Future Work ................................ ................................ .................... 103 LIST OF REFERENCES ................................ ................................ ............................ 104 BIOGRAPHICAL SKETCH ................................ ................................ ......................... 106

PAGE 7

7 LIST OF TABLES Table page 2 1 Fundamental functions ................................ ................................ ...................... 39 5 1 Maximum stress (Von Mises) and displacement magnitudes for Example 1 ..... 78 5 2 Contact region stress and displacement for Example 2 ................................ ..... 78 5 3 Maximum contact region stresses for example 3 ................................ ............... 78 5 4 Maximum contact region stresses for Example 4 ................................ .............. 79 5 5 Maximum contact region stresses for Example 5 ................................ .............. 79 5 6 Maximum contact region stresses for Example 6 ................................ .............. 80

PAGE 8

8 LIST OF FIGURES Fig ure page 2 1 Non conforming surfaces in contact at point O ................................ .................. 39 2 2 Non conforming surfaces in contact at point O ................................ .................. 40 2 3 A half plane subject to line normal and tangential forces ................................ ... 40 2 4 Pressure distribution for contact of cylinders ................................ ..................... 41 2 5 Surface displacement in normal direction ................................ .......................... 41 3 1 Model representing mesh independent method with a non conforming background mesh ................................ ................................ .............................. 48 4 1 Approximate step function H 2 of part 2 ................................ .............................. 56 4 2 Division of computation domain for implementation of contact problems ........... 56 4 3 Tolerance limit criterion for detecting contact elements ................................ ..... 57 5 1 Comparison: Co ordinate system for reference [20] and IBFEM ........................ 80 5 2 Example 1: 2D Plane stress model of a cantilever loaded through contact ....... 80 5 3 Example 1: Schematic of a cantilever loaded wi th an overhang ....................... 81 5 4 Von Mises stress contours: Example with 4 node bi linear elements ................ 81 5 5 Displacement magnitude contours: Example 1 with 4 node bi linear elements ................................ ................................ ................................ ........... 81 5 6 Stress contours: Example 1 with 9 node bi quadratic elements ......................... 82 5 7 Displacemen t contours: Example 1 with 9 node bi quadratic elements ............ 82 5 8 Convergence of maximum Von Mises stress ................................ ..................... 82 5 9 Convergence of ma ximum displacement at loading point ................................ .. 83 5 10 Convergence of maximum displacement at the free end ................................ ... 83 5 12 Von Mises Stress contour s with 3D elements ................................ .................... 84 5 13 Von Mises stress contours: Example 2 with 3D 8Node elements ...................... 85 5 14 Normal stress contours: Example 2 with 3D 8Node elements ..................... 85

PAGE 9

9 5 15 Maximum shears stress contours: Example 8 with 3D 8Node elements ..... 86 5 16 Displacement magnitude contours: Example 2 with 3D 8Node elements .......... 86 5 17 Convergence of maximum Von Mises stress in contact of blocks using 8Node elements ................................ ................................ ................................ 87 5 18 Example 3: 2D Quarter symmetry contact of cylinder and plate ....................... 87 5 19 Stress contours: Example 3 with 4 node bi linear elem ents, a. Von Mises stress distribution, b. Zoomed in View ................................ ............................... 88 5 20 Stress contours: Example 3 with 4 node bi linear elements, a. Normal stress b. Shear stress ................................ ................................ ........... 88 5 21 Displacement contours: Example 3 with 4 node bi linear elements, a. magnitude distribution, b. Zoomed in view ................................ ......................... 89 5 22 Stress contours: Example 3 with 9 node bi quadratic elements, a. Von Mises stress distribution, b. Zoomed in View ................................ ..................... 89 5 23 Stress contours: Example 3 using IBFEM wit h 9 node bi quadratic elements, a. Normal stress b. Shear stress ................................ ........................... 90 5 24 Displacement contours: Example 3 using IBFEM with 9 node bi qu adratic elements, a. magnitude distribution, b. Zoomed in view ................................ .... 90 5 25 Convergence for maximum contact pressure in cylinder flat contact ................. 91 5 26 Example 5: 2D Quarter symmetry contact of cylinders ................................ ..... 91 5 27 Stress contours: Example 4 with 4 node bi linear elements, a. Von Mises stress distribution, b. Zoomed in view ................................ ................................ 92 5 28 Stress contours: Example 4 with 4 node bi linear elements, a. Normal stress b. Shear stress ................................ ................................ ........... 92 5 29 Displacement contours: Example 4 with 4 node bi linear elements, a. magnitude distribution, b. Zoomed in view ................................ ......................... 93 5 30 Stress contours: Example 4 with 9 nod e bi quadratic elements, a. Von Mises stress distribution, b. Zoomed in view ................................ ................................ 93 5 31 Stress contours: Example 4 with 9 node bi quadratic elements a. Normal stress b. Shear stress ................................ ................................ ........... 94 5 32 Displacement contours: Example 4 with 9 node bi quadratic elements ............ 94

PAGE 10

10 5 33 Convergence study for maximum contact pressure in cylindrical contact .......... 95 5 34 Example 5: a. A 3D symmetrical contact of sphere and plate ........................... 95 5 35 Stress contours: Example 6 using IBFEM with 8 node hexagonal elements, a. Von Mises stress distribution, b. Zoomed in view ................................ .......... 96 5 36 Stress contours: Example 5 using IBFEM wi th 8 node hexagonal elements, a. Normal stress b. Shear stress ................................ ........................... 96 5 37 Displacement contours: Example 5 with 8 node hexagonal elements, a. magnitude distribution, b. Zoomed in view ................................ ......................... 97 5 38 Convergence study for maximum contact pressure in sphere to plate contact .. 97 5 39 Example 7: A 3D symmetrical contact of spheres ................................ ............ 98 5 40 Stress contours: Example 6 with 8 node hexagonal elements, a. Von Mises stress distribution, b. Zoomed in view ................................ ................................ 98 5 41 Stress contours: Example 6 with 8 node hexagonal elements, a. Normal stress b. Shear stress ................................ ................................ ........... 99 5 42 Displacement contours: Example 6 with 8 node hexagonal elements, a. magnitude distribution, b. Zoomed in view ................................ ......................... 99 5 43 Convergence study for maximum contact pressure in spherical contact .......... 100 5 44 Variation of maximum contact pressure with applied total load in spherical contact 100

PAGE 11

11 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science CONTACT MODELING USING IMPLICIT BOUNDARY FINITE ELEMENT APPROACH By A nand Shrikant Tupsak hare May 2012 Chair: Ashok Kumar Major: Mechanical Engineering Several structured mesh based finite element analysis techniques have been developed recently to avoid mesh generation difficulties. Structured meshe s are composed of reg ular shaped uniform elements which are easy to generate automatically as compared to mesh generation in traditional FE software In the Implicit B oundary Finite Element M ethod (IB FE M) used in this work solution structures for the trial and test functions are constructed usin g approximate step functions of the boundaries such that allows boundary conditions to be imposed even when there are no nodes on the boundary. In this study, Implicit Boundary Finite Element Method (IBFEM) is further extended to obtain solution for tie d contact problems where the parts in contact stay in contact without sliding or separating during the deformation. The solution structure presented previously by Burla et al. fo r treatment of material di scontinuity in a structured mesh framework was impro ved to handle point contact between two bodies. In this approach, contacting bodies are modeled as an assembly of parts and separate mesh is constructed f o r each part. At the contact elements of the mesh belonging to the contacting parts overlap. A soluti on structure is constructed such that the displacements

PAGE 12

12 are continuous at the contact. A t the beginning of contact, non conforming bodies with suffici ently large relative curvatures touch each other at a point or along a line. With increasing load, they d eform and conta ct increases over a finite area which is usually extremely small in scale as compared to the overall dimensions of the contacting bodies. Contact stresses in this contact patch are significantly hig her than the rest of the regions giving ris e to very steep stress gradients at the boundary of the contact patch In the current work, a n attempt is made to capture these stresses as accurately as possible using sufficiently dense global mesh so as to model the contact patch width to near accurate dimensions as calculated from the analytical solutions Ideally, for the se analyses to be more efficient, localized refinement of the involved meshes is necessary Here we focus only on the case in which identical linear elastic material properties exist f or both the contacting bodies. Several tied contact problems including classical point and line contact problems are used to study the effectiveness of our approach used here in modeling contact using implicit boundary finite element method

PAGE 13

13 CHAPTER 1 INT RODUCTION 1.1 Overview Analytical solutions, though available for analyses of several engineering problems, are difficult to establish for complicated domains and hence numerical methods are needed to obtain better or accurate solutions in such cases The Finite Element Method (FEM) is one of the most widely used numerical methods in engineering analysis due to its applicability to solve a variety of physical problems as well as to handle complex geometry and boundary conditions. FEM subdivides the analysis domain into multiple finite elements which approximate its geometry. FEM approx imates the field variable on these finite elements using piecewise continuous approximations or interpolations. Even though FEM is a well established technique, it has certain limitations. D iscontinuity of gradients across the elements is common in FEM as approximation schemes used are typically continuous, which necessitates the smoothing of results. Commercial traditional FEA packages these days do have well developed a utomat ic mesh generation algorithms for 2D problems but those still can be unreliable for complex geometries and 3D geometries often resulting in po or or distorted elements in certain regions that can lead to large errors in the solution. It requires s ignifican t amount of user intervention to correct such problems making mesh generation the most time consuming step in the design process. These problems associated with mesh generation can be avoided w ith the use of a structured mesh A structured mesh is much eas ier to generate than a finite element mesh and all the elements in such a mesh have regular geometry (rectangles/cubo ids). The use of structured meshe s for

PAGE 14

14 analysis eliminates the need to generate a conforming mesh which is a significant advantage over tra ditional FEM. When structured meshes are used to approximate solution in the analysis domain, the boundaries of the domain need to be defined in dependent of the structured mesh and this can lead to several elements on the boundary which hav e partial volum e inside the analysis domain Special techniques are therefore needed to compute the stiffness matrices and load vectors in these partial elements. When nodes are not on the boundary, the traditional methods used in FEM for applying essential boundary cond itions cannot be used and new techniques are needed for applying e ssential boundary conditions. In Implicit Boundary Method (IBM) solution structures are constructed by using the approximate Heaviside step functions based on the implicit equations of the boundaries for exact imposition of essential boundary conditions. Implicit Boundary Method previously used by Burla et al. [13] for treatment of material discontinuity is extended in this work by adding a capability to model point and line contact, using approximate step functions based on the implicit equations of boundaries of the bodies in contact In addition, a detailed overview of the theory behind classical contacts is also presented. Deriving mathematical equations for these cases involve complex procedures which can be looked up from the theory of elasticity. The equations for contact patch deformation s pressure distribution, and contact stresses for the cases of static loading were originally developed and derived by Hertz in 1881, [21]. Referen ce [14] includes a n English translation of the same Later additions to th e understanding of this problem can be found in references [22], [23], [24], [25]. T he contact stresses are very high and localized only in a very tiny region of the contact

PAGE 15

15 patch of the contacting bodies as compared to the rest of the region s. Hence in this work we try to use very dense global mesh It was attempted that the elements in the contact region are smaller than the size of the contact half width of the contact patch. Ideal ly localized mesh refinement must be used for computational efficiency of the analyses. The performance and capability of the improved solution structure are studied using several problems involving point or line contact in both 2D and 3D with dense global meshes As stated above, ideally t hese require only local resolution of solution near r egions of stress concentrations Convergence analyse s are also performed using these ex amples to study the behavior of various element types 4Node bi linear and 9Node quadratic quadrilateral elements were used for 2D analyses. I t has been demonstrated that 9Node quadrilateral element s provide better and much more near accurate solutions with fewer number of elements when compared to 4Node bi linear elements as expected. It shall be noted that in this work, mesh and grid are considered as equivalent words to represent the background mesh. 1.2 Goals and Objectives The goal of this research is to extend Implicit Boundary Finite Element Method by implementing point and line contact in the originally developed solution structure for analysis of engineering boundary value problems i nvolving material discontinuity. The main objectives of this thesis are : D evelop a mesh independent analysis scheme for modeling contact between parts in an assembly ; S tudy the effectiveness of solution structure originally developed for modeling multi material systems in modeling contact and ensuring dis placement continuity at contact ;

PAGE 16

16 S tudy the accuracy of stress and strain distribution near cont act when modeled using this solution structure 1.3 Outline The rest of the thesis is organized as follows: Chapter 2 discusses He rtz theory of contact mechanics for classical contact problems. Formulations for plane problems and axi symmetric problems alo ng with expressions for surface displacement and stresses involved with contact problems are mentioned. A brief description about some of the numerical techniques for treatment of contact is also included Chapter 3 presents a su rvey on structured grid met hods and meshless methods. It also explains the main idea s behind the Implicit Boundar y Finite Element Method (IBFEM) Chapter 4 discusses the construction of solution structure and its implementation for contact modeling These techniques are implemented using the JAVA programming language in software called IBFEM. Chapter 5 presents results including s everal problems in classical contact mechanics involving localized stress concentrations. Convergence analyses are presented to show the behavior of the imp roved solution structure in terms of accuracy and efficiency. Finally in C hapter 6 conclusions are drawn from the nume rical experiments performed in C hapter 5 The advantages and shortcoming s of proposed solution structure are discussed. The scope for f uture work is presented at the end of the thesis

PAGE 17

17 CHAPTER 2 OVERVIEW OF CLASSICAL CONTACT PROBLEMS 2.1 Theory of Contact Initially when brought in contact, non conforming bodies hav ing sufficiently large relative curvatures touch each other at a point or along a line. As the load is increased, the bodies deform and contact increases ove r a finite area, extremely small in scale as compared to the contacting bodies. In the current chapter the underlying assumptions and conditions are explained under which the analytical solutions for classical contact problems are valid. These were used a s a benchmark to validate and compare the results obtained using IBFEM. At the end of this chapter, w e also try to describe briefly how c ontact is modeled traditionally in FEM using methods such as Lagrange multi plier method and penalty method. The data me ntioned below can be validated from reference [1]. 2.1.1 Normal C ontact of E lastic S olids We assume the origin of the Cartesian coordinate system to be at the point of first contact between the two bodies with the x y plane as the common tangent plane to t he two contacting surfaces T he z axis lying along the common normal to the surfaces, is directed positively into the lower solid, as shown in the Figure 2. We work with the following a ssumptions 1. Each surface is considered to be smooth at both the m icro scale and the macro scale T he two contacting surfaces are assumed to be elliptical and represented as and ( 2 1 )

PAGE 18

18 where and are the p rincipal radii of curvature of the s urfaces at the origin. These represent maximum and minimum values of radius of curvature of all possible cross sections of the profile. 2. The surfaces are assumed to be frictionless so that only normal pressure is transmitted between them. If the separation betw een the surfaces (h) is given as : I t can be shown that ( 2 2 ) w here, A and B are positive constants and and represent the principal relative radii of curvature. Also and must be en ough large to eliminate the high er order terms. Also, and ( 2 3 ) Convex surfaces are to have positive radius of curvature and concave to be negative. If the axes and of principal curvature of each surface, are at an angle to each other, then as shown in Appendix 2 of [1] ( 2 4 ) ( 2 5 ) An equivalent radius is defined a s : ( 2 6 )

PAGE 19

19 Eq uation 2 2 shows that the elliptical contours of cons tant gap h between the un deformed surfaces have len gth of axes in the ratio From symmetry of Equation 2 2 about the origin it is evident that the contact region extend s an equal distance on either side of the origin as the defo rmation progresses. Following cases can be considered when a normal compressive load is applied to the two solids in contact Solids of revolution with radii R 1 and R 2 (e.g Spheres) Both, the c ontours of constant separation between the contacting sur faces before loading and the contact patch area after loading, are circular and centered at the origin. Parallel (y axes) cylinders in contact with radii R 1 and R 2 In unlo aded state, c ontours of constant separation are straight lines parall el to the y axis. Under loading the contact patch is a narro w strip parallel to the y axis. Perpendicular (y axes) cylinders in contact with equal radii R In unloaded state, c ontours of constant separation are circles. This is identical to the case of a sphere of the same radius R in contact with a plane surface. 2.1.2 Fundamentals of Hertz Theory F igure 2 2 shows 2 solids undergoin g deformation under actio n of normal compressive load P. Equation 2 2 gives the separation between two corresponding surface points and b efore deformation begins During this process

PAGE 20

20 of compression and are the respective displacements for the distant points and towards O and parallel to the z axis. The overlapping dotted profiles are the case w hen the solids do not deform after load application. and represent the displacements of each body parallel to Oz due to contact pressure These are me asured positive into each body and relative to the points and A f ter deformation if and are coincident in the contact patch, then we can write: ( 2 7 ) Using Equations 2 1 and 2 2 for elastic displacements, with (x, y) = common coordinates of and projected on the x y plane: When and are in the contact region, ( 2 8 ) When and are outside the contact region, ( 2 9 ) T he normal surface displacements must satisfy the a bove conditions given by Equations 2 8 and 2 9 In the Hertz theory, to calculate f or local deformations each body is considered as an elastic half space loaded over a small elliptical region to treat the contact stresses separately from general stresses in the rest of the regions M ethods for solving boundary value problems involving half spaces were used to obtain solutions for contact problems In order for this simplification to be applicable following conditions must be met The significant dimensions of the contact area must be small compared to: Dimensions of each body Thi s ensure s that the boundaries of the infinite solid considered are far away from the regions of stress concentration while calculating the contact stress field.

PAGE 21

21 The relative radii of curvature of the surfaces This is to enforce that outside the contact re gion the surfaces shall approximate roughly to plane surfaces of half spaces and sufficiently small strains occur in the contact region for the problem to be in the elastic range. Half spaces Contacting bodies are considered as planar half spaces N ormal tractions at the interface act parallel to the z axis and tangential tractions are in the x y plane. Let us have the following a contact area half width, R Relative radius of curvature, R 1 R 2 Radii of the 2 bodies, Significant dimensions both laterally and in depth Then u sing section for normal contact of solids and above data in the current section assumptions made in the Hertz theory can be summarized as Contacting surfaces to be continuous and non conforming: The strains involved are small: Contacting solid s can be considered as elastic half space s : and F rictionless surfaces : representing zero tangential pressures The definition of the contact proble m as a problem in elasticity now can be stated as We need to solve for t he distribution of mutual pressure acting over an area S on the surface of tw o elastic half spaces that produces normal displacements of t h e surfaces satisfying Equation 2 8 within the contact zone S and Equation 2 9 outside it. 2.1.3 Growth of C ontact A rea, D eformation and S tresses with I ncreasing A pplied L oad Now, we try to stud y how the contact area, stress and deformation is expected to grow with increasing load along with their dependence on curvature and elastic moduli.

PAGE 22

22 We know that the contact grows equally in x and y direction. So we take as the variable coor dinate and write using F igure 2 1, and writing Equation 2 2 in non dimensional form: ( 2 10 ) Modifying above E quation 2 10 by and as the expression for deformation within the conta ct zone can be wri tten as: ( 2 11 ) Considering small deformations with the strain is characterized by ratio pressure acting mutually on each body, divided by the elas tic modulus. Hence ( 2 12 ) This means, contact pressu re and related stresses are direct ly proportion al to the linear dimension of the contact area. Contact of cylinders Here, l oad per unit axial length P ut in Equation 2 12, and ( 2 13 ) This implies that for cylinders, contact width and contact pressure are direc tly proportional to square root of the applied load P. Contact of spheres

PAGE 23

23 Here, compressive load Put in Equation 2 12, and ( 2 14 ) This implies that for spheres, radius of contact region and contact pressure, both are directly proportional to cube root of the a pplied load P. For 3D contact, since the compressions and are proportional to the local indentations and we have, ( 2 15 ) This shows that the approach of the bodies in contact region is proportional to 2.2 Plane P roblems : F ormulation and A nalytical So lutions The s olution procedure for plane problems in volves two steps (1) t o derive the expression for the contact pressure distribution ; ( 2 ) t o find int ernal stress and displacement fields A detailed generalized formulation is mentioned in t his section [2] which is applicable to various plane problems in volving contact. As shown in Figure 2 3, l et be the applied normal and tangential distributed traction on the sur faces in contact respectively. P, Q respectively be the applied normal and shear forces applied to brought the surfaces in contact. Relative displacement or the amount of overlap between the two bodies while penetrating freely into each other, between the two corresponding points on the contacting bodies in z direction can be used as: ( 2 16 )

PAGE 24

24 W here : The composite compliance ( 2 17 ) ( 2 18 ) ( 2 19 ) Here we deal with contactin g bodies of similar materi als and contact is frictio nless i.e. perfectly lubricated. This eliminates the tangential tractions Hence, the second term on r ight hand side of above Equation 2 16 vanishes and the transformed eq uation is then For in the contact patch ( 2 20 ) Inverting E quation 2 20 it is shown by Muskhe lishvili, Ilills, Nowell and Sackfield (1993 ) that for a given separation between the two bodies, the pressure d istribution can be written as: ( 2 21 ) The fundamental function depends upon the behavior of at the end points of contact patch as seen from T able 4 1 For consistency we must have: ( 2 22 ) C ontact of c ylinders : Notations used are given here Contact radius = a, p 0 = Peak pressure P = Total external load Here we briefly discuss the method for solutions to cylindrical contact as described in section 4.2 c, page 99 of [1].

PAGE 25

25 Consider two cylinders contacting each other over strip width 2a along line with force per unit l ength P and their axes being parallel to y axis stated in the co ordinate system defined previously at origin O located at the start point of contact The contact region is assumed to have length b (major axis of the ellipse) much greater than length a (mi nor axis of the ellipse) Using Equation 2 1, the separation between corresponding points on the unloaded cylinder s becomes: ( 2 23 ) and ( 2 24 ) For cont act region after loading, Equation 2 2 changes to: ( 2 25 ) Outside the contact region, using Equation 2 3 we can write: ( 2 26 ) To find the local contact stresses, differenti ating Equation 2 26 with respect to surface gradients are: ( 2 27 ) The surface gradient due to pressure on contact strip is ( 2 28 ) Here, s = distance along x axis from origin, q(x) = Tangential tractions Now, same pressure acts on both the surfaces, hence:

PAGE 26

26 ( 2 29 ) From Eq uation 2 27 and Eq uation 2 28 : ( 2 30 ) Details of solution procedure for above type of equation are found i n section 2.7 of reference [1 ]. We then get, ( 2 31 ) Until and unless, semi contact width is in relation with is not uniquely defined. Hence, it was attempted to develop a relation between and As pressure will be positive in contact region, ( 2 32 ) When will exceed RHS of Equation 2 3 2 we get: at This leads to infinite surface gradients just outside the loaded region ( pg.36 [1 ] ) The deformed profile given by this type of pressure distribution with such surface gradients are not in line with Eq uation 2 2 which says, contact should not occur outside the loaded area. Hence, we can safely say: OR ( 2 33 ) we here have and ( 2 34 ) Put P from Equation 2 33 in Equation 2 31 it gives:

PAGE 27

27 or ( 2 35 ) This pressure distribution is shown schematically in Figure 2 4. We n ote that the c ontact pressure becomes zero at the contact edge. Maximum pressure at, mean pressure ( 2 36 ) a. Full field stresses: [2] T he analytical solutions for stresses is given as follows ( 2 37 ) ( 2 38 ) ( 2 39 ) ( 2 40 ) ( 2 41 ) Superscript stresses developed only due to the normal traction applied. Tangential tractions are assumed to be absent due to frictionless contact or due to use of identical material properties for the contacting parts is the largest root of the quadratic equation: ( 2 42 )

PAGE 28

28 i.e. ( 2 43 ) In contact patch, as z and s tend to zero, y/s achieves a finite value giving ( 2 44 ) Outside the contact patch, as z tends to zero, s has a finite value giving ( 2 45 ) The expressions for surface stresses in the contact re gion are as follows ( 2 46 ) Exterior region to cont act patch has no normal surface stresses and only tangential stresses due to tangential shearing tractions, if present, will occur. b. Relative surface displacements: la w and integrating to solve for displacements [2]. It is worth noting that f ollowing displacement equations are exact if the two bodies in contact are elastically similar i.e. Ratio in which case becomes zero. They may also be applied to elastically dissimilar bodies with approximate correctness of the displacement values as they are ion of the Eq uation 2 16 above, given as:

PAGE 29

29 (For details, refer (Ch.2, [2]) ( 2 47 ) Here, for elastically similar bodies, becomes zero or for perfectly lubricated contact q(x) becomes zero and only first term in above equation remains on RHS. This transformed equation is further used in the contact patch to obtain equations for stresses in and out of the contact region as described in Section 2.3 of (Ch.2, [2]) b 1 Normal surface displacement ( ) ( 2 48 ) ( 2 49 ) where, ( 2 50 ) Determination of relation between are arbitrary, but are chosen such that displacements are continuous at Putting OR i n Equations 2 48 and 2 49 and then equating both, i.e. ( 2 51 ) With this relation, Equation 2 49 changes to following with any arbitrary value for C 0 as ( 2 52 ) b 2. Tangential surface displacement ( )

PAGE 30

30 is to be taken as( ) to assure anti symmetry of tangential displacements. ( 2 53 ) ( 2 54 ) ( 2 55 ) ( 2 56 ) where, Dunders Parameter = ( 2 57 ) with ( 2 58 ) ( 2 59 ) 2.3 Axi sym metric P roblems F ormulation Assumptions of Hertz theory are applicable for axi symmetric problems as well. Contacting solids are assumed as semi infinite bodies with similar materials or in perfectly lubricated contact so that shear tractions can be ignore d. The relationship between the contact pressure and the degree of overlap between the bodies is [ 2] ( 2 60 ) Where, as a transformed displacement is given as

PAGE 31

31 ( 2 61 ) This two stage procedure gives the contact pressure distribution as ( 2 62 ) Details can be found in Ch.2 of reference [2] Contact of s pheres t he H ertz p roblem : Notations used are given here for reference. C ontact radius = a, = Peak pressure, P = Total external load As the contact radius of the circular contact region is assumed to be much smaller than the radii of the contacting body, we can have the separation between corresponding point s on the two contacti ng bodies with = Relative curvature ( 2 63 ) and are the radii of the two contacting spheres. Using Equation 2 61 the transformed displacement is given by: ( 2 64 ) This further Equation 2 62 gives: ( 2 65 ) Equation 2 65 gives the most general solution for contacting spheres when and are specified. Note that it also depicts si ngularity at Now as the pressure shall fall off to zero continuously, we have

PAGE 32

32 and hence which further from Equation 2 65 gives ( 2 66 ) Using the equilibrium of normal forces P = Total External load = ( 2 67 ) This would give ( 2 68 ) where ( 2 69 ) Maximum pressure = ( 2 70 ) and Principal radii of curv ature of the surfaces at the origin, which are the maximum and minimum values of radius of curvature of all possible cross sections of the profile. Equation 2 2 within the contact region becomes ( 2 71 ) Where, = relative curvature The pressure distribution satisfying Equation 2 7 1 is give n by Equation 2 7 0 a bove. This pressure distribution is shown to give below mentioned n ormal displacements as per Equation 3 41a [1] (p 61),

PAGE 33

33 ( 2 72 ) The pressure acting on the second body is equal to that on the first, So that we can write from Equation 2 71 substituting equations for in above form, ( 2 73 ) Where, and ( 2 74 ) Further using this with where the contact circle radius is given by ( 2 75 ) This is the expression for Contact radius with The mutual approach of the distant points in the two solids is given by: ( 2 76 ) This is the expression for Approach of the two contacting bodies. This are the absolute value s for and P. Using the pressure distribution of the form, with n=1/2 for Hertz pressure in case of solids of revolutions such as contact of sphere s the closed form solution f or the normal displacement is given by: where is along

PAGE 34

34 the circumference of the loaded circle in contact region and is the angle of revolution (p 60, 61, [1]). With this, f ollowing results are obtained. a. S tresses: [2] Using expression for pressure distribution from, in presence of only normal tractions (details of de rivation in section 2.5. 1 of [2]), (alternate form in [ 20 ] is used for Ch 5 here) The stresses in the cont act patch ( 2 77 ) ( 2 78 ) ( 2 79 ) where ( 2 80 ) The stresses out of the contact patch ( 2 81 ) ( 2 82 ) ( 2 83 ) b. Relative surface displacements: b 1. Normal surface displacements:

PAGE 35

35 Schematic representation of surfa ce displacement in normal direction can be seen in Figure 2 5. Reference [1] has following details of deriv ations for surface displacement. D eriva tion on (p 60, 61, [1]) gives, N ormal surface displacement ( ), fo r points inside contact patch : ( 2 84 ) D eriva tion on (p. 61, [1]) gives, Normal surface displacement for points outside contact region: ( 2 85 ) b 2 Tangential surface displacements: D erivation on (pgs. 61, [1]) gives, Tangential displacement for points insid e contact region : ( 2 86 ) Derivation on (p 61, 62, [1]) gives Tangential displacement for p oints outside contact region : ( 2 87 ) 2.4 Traditional M ethods for N umerical C ontact M odeling To satisfy the impenetrability condition during contact, we need to enfo rce kinematic constraints which will prevent penetration of one boundary through the other

PAGE 36

36 and provide traction and displacement continuity between the bodies during persistent contacts. To solve a contact problem we need to have: (a) an algorithm to sea rch points on a boundary segment that are in interact ion with those on another boundary segment and (b) constraints to prevent penetration and transmit the traction correctly between the contacting bodies. There are several methods used for implementing contact constraints between two bodies. The most fundamental are Penalty and Lagrange Multiplier Approaches [3] We try to explain the underlying concept behind these two approaches using Node to Node contact. The contacting nodes are found by comparing the vertical position of each node pair. Denoting one part ( preferably t he upper part) as slave body contacting (lower) part as master body ( 2 88 ) T he contact constraint can then be realized as ( 2 89 ) Penetration can occur where the constraint condition is not imposed. These constraint conditions are used fo r any nodal pair or element in which the gap g is negative or zero. Zero is sometimes defined using specific tolerance values in some cases. Following are the ways to perform this function:

PAGE 37

37 2.4.1 Lagrange M ulti plier F orm Here we multiply the gap condition given in E quation 2 .89 by a multiplier. In a variational form: ( 2 90 ) Where is the surface traction, shows the pos ition on surface of the slave body, represents the position on the surface of the master body is a Lagrange multiplier force and is the gap given by Equation 2 8 9 First variation of is used in the variational equations to solve the co ntact problem which is given as : ( 2 91 ) This clears that Linearizing a bove equation we get a tangent matrix term used in a Newton solution process. Following Equation 2 92 shows the final tangent and residual f or nodal contact element. ( 2 92 ) Lagrang e multiplier method introduces additional v ariables to enforce the contact constraints directly and exactly and thus requires additional efforts to solve for multipliers. In addition, the equation of the Lagrange multipliers introduces zeros in the diagonal of the system equations posing additional difficulties in direct solution techniques to be used as it becomes necessary to avoid division by zero.

PAGE 38

38 2.4.2 Penalty F orm This serves as a way out from dealing with a zero diagonal from a Lagrange mu ltiplier method. Here the c ontact term is represented as follows : ( 2 93 ) Here, serves as the penalty parameter and is the constraint for contact establishment. Matrix equation for a nodal/ elemental pair is transformed as : ( 2 94 ) As clearly visible from above formulation of the Penalty approach, it avoids the need of additional variables by introducing an approximation of the constraint condition. An additional term enters in the weak form of the governing equations, which penalizes the dissatisfaction of the constraint condition by a large positive penalty parameter. Theoretically, as the penalty parameter tends to infinity, the constraint is enforced exactly. A disadvantage of penalty approach is that the resulting system of equations can become ill conditioned with increase in penalty parameter value. Thus, the choi ce of an appropriate penalty parameter becomes a balancing act between accuracy and stability of the penalty approach. A solution to disadvantages of above approaches is chosen as an Augmented Lagrangian Method which serves as a combin ation of the regula rizing effect of the penalty method and also carries the exact enforcement of constraints contrib uted from the Lagrange approach without introducing instability into the system.

PAGE 39

39 Table 2 1. Fundamental functions Behavior at end x = a x = a N S 0 S N 0 S S Non Zero N N 0 S: Singularity, N : Non Singu larity Figure 2 1. Non conforming surfaces in contact at point O

PAGE 40

40 Figure 2 2 Non conforming surfaces in contact at point O Figure 2 3. A half plane subject to line normal and tangential forces

PAGE 41

41 Figure 2 4. Pressure distribution for contact of cylin ders Figure 2 5. Surface displacement in normal direction

PAGE 42

42 CHAPTER 3 MESH INDEPENDENT FIN I TE ELEMENT ANALYSIS 3.1 The Finite Element Method Partial differential equations are used to demonstrate physical phenomen on such as solid mechanics, heat transfer, electrostatics, fluid flow, and so on These equations usually can be solved with minimum efforts to gain ana lytical or exact solutions only if the geometries are simple and the problem involves only certain restricted loading conditions Since complex geometry and loading conditions are common for real physical or engineer ing situations, alternative methods are necessary to solve for such problems. One of the most popular and well established method, both in academia and industry, which is highly dependent on domain discretization using a conforming mesh, is the finite elem ent method (FEM) serving as a numerical method for solving partial differential equations. But as a drawback, often mesh generators are not sufficiently reliable for 3D domain This inability contaminates the mesh with few highly distorted elements which f urther leads to inaccurate solutions 3.2 Mesh Independent and Meshless Methods Several meshless techniques have been developed with a goal of avoiding complex mesh generation [4]. Here dispersed nodes within a domain are used for interpolations without co nnecting them to form elements. These approximations do not possess the Kron ecker delta property. T he value of the approximation at a node is equal to the nodal value of the approximation. Lagrange multiplier approaches, penalty methods are used to deal wi th t he difficulties faced in applying essential boundary conditions in such analyses. Examples of meshless methods include techniques such as Hp clouds method by Duarte et al., reproducing kernel particle method (RKPM) by

PAGE 43

43 L iu et al., Belytschko et al. and so on Element free Galerkin method (EFGM) [5,6] used moving least square approximation method developed by Nayroles et al. for solving Babuska and Melenk [18 ] presented a Partition of Unity Method (PUM) It serves as a general approach to user defined domain with required interpolations Smoothed particle hydrodynamics (SPH) is one of the older meshless methods in literat ure (Lucy, 1977; Gingold, 1982). I t was first used to model astrophysical phenomenon and later was modified to include continuum mechanics problems (Monahan, 1982). The most desired advantage of mesh independent methods is in using a back ground mesh or a structured grid. In such methods, it is not required for the mesh to conform to the g eometry. This can be seen from F igure 3 1 Easy and automatically generated structured grids having regular and undistorted elements are another alternative to avoid traditional mesh generation which at times b ecomes heavily cumbersome. Here, geometry gets independently represented using implicit boundary equations, as used by Belytschko et al. [7]. It is also possible to use variety of interpolation combinations with structured grids. A solution structure prese nted by Kantorovich and Krylov using implicit equations of boundaries of analysis domain ensures that essential boundary conditions are guaranteed to be satisfied at these boundaries. Rvachev used R functions using Boolean operations between simple primiti ves to define implicit equations for domains. Shapiro and Tsukanov extended the R function approach to solve problem s with time varying geometries and made use of transfinite Lagrangian interpolation to apply essential boundary conditions [8]. Hollig et

PAGE 44

44 al proposed the Web method involving use of weighted extended B spline basis to g uarantee higher order continuity in solutions [9]. This method has been made use with R functions as implicit equations for constructing solution structures. In this study, the implicit boundary finite element method (IBFEM), developed by Kumar et al [10] [13], [14] is used which uses solution structures constructed using approximate step functions of the implicit equations of the boundaries. Approximate step functions possess a unit value inside the analysis domain and transition sharply to zero at the boundary within a very thin transition band It becomes advantageous to use step function in the solution structure as it makes possible all the internal elements having identica l element stiffness matrix. IBFEM has also been used with traditional FEM interpolations ( Lagrange interpolation s) as well as with uniform B spline approximations [12]. 3.3 Implicit Boundary Finite Element Method In the implicit boundary method, approxima te step functions of the boundaries are used to construct the solution structure for test and trial functions. These solution structures guarantee fulfillment of the essential boundary conditions at these boundaries. Consider a tr ial function which is bound to satisfy the boundary condition along the boundary which is a part or subset of the boundary of This allows us to write the trial function as: ( 3 1 ) Here is weighting function constructed such that it has zero value at the boundaries where essential boundary conditi ons are specified. It is worth noting that

PAGE 45

45 the weighting function referred to here as the Dirichlet function or D function, is not zero along the entire boundary of the domain. Instead it is assumed to be zero only at the boundar ies with essential boundary conditions. An advantage of such weighting function is that if are constructed as step functions, it becomes possible to define them locally within elements containing such bo undaries The solution structure for the t rial function defined in Equation 3 1 is guaranteed to satisfy the condition along these boundaries for any The boundary value function is to be defined such that at the boundaries it has the prescribed values The variable part of the solution structure is defined by piece wise interpolation over the elements of the corresponding grid. The Dirichlet function is a step function which is constructed such that at the boundary it takes a value of zero with a non zero gradient and has a unit value away from the boundary. T hese functions as sist in imposing th e Dirichlet boundary conditions hence they have been referred to as Dirichlet functions or D functions. As mentioned above, t hese functions are constructed as approximate step functions [15] using an implicit function of the boundary as follows : ( 3 2 ) Hollig [16] made use of similar functions for constructing solution structure for B spline FEM with large values of This allowed t he weighting function transition from

PAGE 46

46 zero to one over a wide band along the boundary. In the im plicit boundary method, the limit of is used. This makes the D function approximate the Heaviside step function having a unit value inside the domain where and is zero at the bo undary and outside The magnitude of 10 5 or smaller is used in the implementation. Moving boundaries were represented using a ppro ximate step functions in level set method [15] With use of approximate s tep functions as the D functions in implicit boundary method, only the boundary elements get affected by this function and all internal elements behave exactly similar to the traditional finite elements. Furthermore, this allows the D function to be define d locally within the boundary elements having boundaries with specified boundary conditions. All the internal elements have identical shape and size if the grid is uniform and therefore all these elements will have the same stiffness matrix. The variable p art of the solution structure, gets defined using finite element like piece wise interpolation or using smoother B spline or other approximations. Using finite element shape functions yields a piece wise interpolation which is continuous B spline basis can provide a high degree of continuity depending on the degree of the polynomials used (e.g. linear, quadratic, etc. ). Uniform B spline elements applicable to structured grids have been studied [17] to demo nstrate that the implicit boundary method can be used to apply boundary conditions even when the approximation used does not have Kronecker delta properties. This clears that even though the nodal values are not equal to the B spline approximation values a t the nodes, the solution structur e in Equation 3 1 guarantees that the essential boundary conditions are imposed.

PAGE 47

47 Special techniques are used to comput e the stiffness of boundary elements. This is because the D function in the solution structure can have steep gradients near the boundary. The terms in the volume integral for computing stiffness that contain the gradient of the D function are converted into surface integrals. The details of the stiffness computation technique for boundary elements can be f ound in [17]. 3.4 Development of M odified W eak F orm for L inear E lasticity in IBFEM Details about constructing the weak form and the stiffness matrix can be found in [11]. The modified weak form or principle of virtual work for a domain being the boundary with essential boundary condition and being boundary for traction application, is: ( 3 3 ) In Equation 3 3 is virtual strain and is the v irtual displacement. Thus a trial solution using the homogeneous portion (superscript s) and boundary value portion (superscript a) can be then constructed as: ( 3 4 ) H ere is a vector representing boun dary value function, are grid variables and is a diagonal matrix D matrix is constructed by Dirichlet functions in a way that makes the variable part of solution representing grid variables vanish at bou ndaries. The same Dirichlet functions used to create test functions as well The stresses and strains now can be written as: (Strain) ( 3 5 ) ( 3 6 )

PAGE 48

48 Here matrix C holds the elasticity coefficient. Substitut ing test function in weak form we then get: ( 3 7 ) Now strain and stress can be obtained as follows: and ( 3 8 ) In above Equation 3 8 vector has the nodal grid variables for element and [B] represents the strain displacement matrix which holds the derivatives of the shape functions This involves discretization of the weak form by assembling the respective equations in global stiffness matrix similar to traditional FEM. I n such a f ormulation variety of basis functions such as Lagrange interpolation, B spline functions can be used Figure 3 1. Model representing mesh independent method with a non conforming background mesh Arbitrary geometry Uniformly structured background mesh

PAGE 49

49 CHAPTER 4 IMPLEMENTATION OF POINT CONTACT USING IB F EM 4.1 Solution S tructure for C ontact Burla et al. (2009) [1 3 ] demonstrated a solu tion structure using IBFEM for micromechanics of structures which makes use of approximate step functions to satisfy the interface conditions. It is required to have continuous displacements and ability to allow discontinuous gradients at the material inte rface boundary while modeling multi material systems Unlike traditional FEM, while using struc tured meshes f or FE analysis, the nodes may not lie on the interface of the contacting objects posing a challenge in imposing interface conditions for material d iscontinuity. In the current work, we try to extend the application of this solution structure to model tied contact of two objects with two distinctly separate meshes overlapping with each other in the contact region The two meshes are named as grid 1 for part 1 and grid 2 for part 2, with displacements denoted by and interpolated piecewise within elements of grid 1 and grid 2 respectively. With this t he solution structure for displacement field for contact of two objects is represented as: where (1 H 2 = H 1 ) ( 4 1 ) In above Equation 4 1, is an approximate step function of part 2 which has a unit value within part 2 and transitions sharply to zero at the boundary of the part 2, as shown in F igure 4 1. Similarly it can be said that H 1 is t he approximate function of part 1 which behaves identical to the function (1 H 2 ) as shown in F igure 4 1. Following is the expression for H 2 which is identical to that of D(x) in Equation 3 2 is used

PAGE 50

50 W ith the step function a s unity inside part 2 the solution is from the grid 2 of part 2 and with as zero (or H 1 as 1) inside part 1 ; the solution is from the grid 1 of part 1 In the very thin banded region at the interface, varies from unity to zero and the solution in this region is the blend of the solutions from grid 1 and grid 2 i.e combination of displacements from both the parts. In this manner, t his solution structure ensures the displacement continuity of the solution throughout the analysis domain. In addition it also makes provision to allow for the normal strains to be discontinuous at the interface which can be clearly seen if gradient of the solution structure is considered [13]. 4.2 Modified W eak F orm for C ontact I mplementation We try to summarize the d etails of the modified weak form found in reference [13 ]. The finite element approximations for displacemen ts within grid elements are: for part 1: for part 2: where, are shape functions and is the i th component of displacement for j th node of the element U sing the solution structure from Equation 4 1, the displacement is: ( 4 2 ) Here, = Elemental nodal displacement vector = nodal values for grid 1, = nodal values for grid 2

PAGE 51

51 The matrices are expressed as follows : ( 4 3 ) ( 4 4 ) ( 4 5 ) ( 4 6 ) Th e virtual displacements are expressed as ( 4 7 ) Differentiating Equatio n 4 2 the strains for 2 D model are given as: ( 4 8 ) Combining the str ain displacement matrices of the two grids, the modified strain displacement matrix is given as: ,where ( 4 9 ) ( 4 10 ) Where, ( 4 11 ) and ( 4 12 )

PAGE 52

52 Using Equation 3 8 and Equations 4 1 to 4 12, substituti on into Equation 3 7 gives the modified weak form for contact implementation as follows ( 4 13 ) In the above E quation 4 13 s ummation is carried out over the number of elements ( NE ) and the number of boundary elements ( NBE ) in the grid used for analysis. A bove solution structure was implemente d simultaneously for 2D as well as 3D environments. 4.3 Nu merical I mplementations for C ontact T he elements in the grid are divided into three types: elements in the grid 1 ( part 1 ) elements in the grid 2 ( part 2 ) and those that contain the interface regio n of the contacting bodies as shown in Fi gure 4 2 As explained in the previous section, the modified strain displacement matrix is split into strain displacement matrix of each grid. Further these matrices were split into the parts depending on gradients of the shape functions and gradients of the approximate step function To make use of these decompositions for computational efficiency improvement as some of the components of the strain displacement matrix vanish in certain region of the analysis domain, the elements in the grids are divided into three types as follows: Mesh 1 for part 1 :

PAGE 53

53 For elements in mesh 1 only the degrees of freedom for mesh corresponding to part 1 are active, therefore the s tr ain displacement matrix becomes : because everywhere in part 1 or mesh 1. This leads to the following stiffness matrix for elements in side part 1 : ( 4 14 ) Mesh 2 for part 2 : For elements in mesh 2 only the degrees of freedom for mesh corresponding to part 2 are active, therefore the str ain displacement matrix becomes : because everywhere in part 2 or mesh 2. This leads to the following stiffness matrix for elements in side part 2 : ( 4 15 ) For both above cases, c o mponent is computed as an area integral for 2D or volume integral for 3D. Interface region: In this region, the approximate step function H 2 transitions from zero to unity. Interface elements are composed of overlapping e lements from both grids in the interface region whic h together contribute to the solution within the interface element s. These elements are common to both the grids and hence the corresponding elements of the original meshes always move together when displacement occurs. This ensures the displacement contin uity at the interface. Degrees of freedom from both meshes corresponding to part 1 and part 2 are active in overlapping elements, therefore the

PAGE 54

54 str ain displacement matrix becomes The stiffness matrix is then expressed as: ( 4 16 ) Each of the term from above matrix is decomposed as a sum of volume integral and a surface integral. We explain this procedure in detail for the component ( 4 17 ) Using the definition of from Equation 4 9 above Equation 4 17 is written as: ( 4 18 ) In Equation 4 18, the first term is non zero only outside part 2 where H 2 is zero. It is integrated by subdividing the eleme nt into triangles or tetrahedra. The remaining 3 terms and involve gradients of the approximate step function and are non zero only in the region of thin interfacial ba nd in which theses gradients are non zero. This allows conversion of the volume integral to surface integral with very small value of the range parameter in use. It is worth noting that even if is made very small, its contribution to the integration procedure is significant due to high gradients of the approximate step function in place.

PAGE 55

55 Criterion for d etecting c ontact e lements : The current work uses a specified tolerance limit on d istance of initial separation between entities of the contacting bodies. If any of the entity i.e. points, lines or surfaces of the contacting bodies are within this limi t of separation from each other then the corresponding elements on the contacting bodi es are designated as contact elements. This is explained using F igure 4 3 If either of the distances a or b is less than the tolerance limit specified then it is determined that a coincident boundary exists between the two contacting parts and the corresp onding overlapping element is design ated as a contact element. T he contact stresses are highly sensitive to the number of elements recognized in contact zone controlled by the above tolerance limit. This also define s the contact patch size depending upon t he region occupied by the contacting elements. The contact patch half width as calcu lated in examples presented in C hapter 5 are extremely small as compared to the overall dimensions of the contacting bodies. To model the contact patch very accurately, we need superfine mesh near the contact region. D ue to limitation of computer memory we could not model the contact patch with really small elements in most of the examples. This led to approximate modeling of the contact patch but nonetheless the technique u sed still gave much better, satisfactory and converging solutions and distributions for contact stresses even with the global mesh refinement in use. This approach used for modeling classic al contact problems using IBFEM, however, has some drawbacks. Thou gh, the maximum stresses along with their proper overall distribution occur almost at the exact location of where they should be as per the analytical solutions almost entire contact force is concentrated in the only element

PAGE 56

56 designated as the contact elem ent at the point of contact. Due to this, entire stress concentration occurs in the only element and the contact stress distribution is not exact though it looks roughly close to the actual distribution. Fluctuations in the stress magnitude are seen at the contact point for maximum contact stresses in the contact patch of the two bodies This however does no t hinder convergence of solutions for contact problems as seen from the convergence analysis in each of the examples in C hapter 5 Figure 4 1. A pproximate step function H 2 of part 2 Figure 4 2. Division of computation domain for implementation of contact problems Boundary ( ) Part 1 Part 2 1 H 2 H 2 Part 1 Part 2 H 1 = 1, H 2 = 0 H 1 = 0, H 2 = 1 Thin band over which H2 transitions from 0 to 1 Contact element at interface

PAGE 57

57 Figure 4 3. Tolerance limit criterion for detecting contact elements Overlapping interfacial element Boundary for part 1 Boundary for part 2 b a

PAGE 58

58 CHAPTER 5 RESULT S A brief discussion about implementation of contact using the implicit boundary f inite element approach in context of modeling contact problems involving local stress concentration for 2D and 3D mode ls was presented in Chapter 4 This chapter discusses the results that were obtained from th e implementation of contact in IBFEM, which va lidate the implementation of the generated finite element code. The results obtained using IBFEM are checked against the analytical solutions for identical models in most of the example below. In E xample 2 results obtained from ABAQUS were used as compari son to IBFEM results Various graphs and figures presented in this chapter help to study convergence of the technique using different element types. The first example of a simple cantilever beam loaded near its free end through a half cylindrical support demonstrates the validity of the solution structure for contact implementation for a problem in which the maximum stress does not occur at the point of contact but somewhere else. This problem has exact solution for elastic displacements and stress es as me ntioned below In order to demonstrate the performance of the solution structure and the implementation technique in the regions of stress concentrations, examples 2 to 6 are presented which involve classical contact problems and the convergence of the sol ution structure is studi ed Note that all t he examples presented here are tied contact examples Example 2 has line/ surface contact respectively in 2D/3D. Examples 3 t o 6 hav e a common characteristic of involving a point or line contact at the beginnin g in unloaded state. As the bodies a re loaded in due course of time;

PAGE 59

59 contact being inherently a non linear problem, shall bring different regions of the contacting bodies in cont act. But as seen from section 4 2, we have only implemented the tied contact of two bodies. These examples employ extra dense global structured mesh. Ideally, to obtain even better solutions for contact stresses and displacements, localized mesh refinement or XFEM or localized mesh refinement along with XFEM is needed mak ing these analyses truly efficient. But nonetheless, current use of dense global mesh also helps to explain the effectiveness of the improved solution structure for classical tied contact problems involving distinctly high stress gradients due to contact stresses in the contact patch of the two contacting bodies. Much o f the computational expenses will be saved with localized mesh refinement as the displacement solution does not change significantly outside the contact patch. The distribution of contact stresses how ever may not look very similar to the analytical solutions when low mesh densities are used in the contact patch, as this region is extremely tiny in the elastic range of the analysis. Nonetheless, even though we are mo deling these contact problems as line ar elastic models, gross distribution of the gross distribution of stresses and displacement is seen to be very close to that given by ABAQUS which uses non linear analysis to model contact It shall be noted that this distribution improves with very fine mesh used in the contact patch which shall be obvious as more elements are present to distribute the contact force and to capture the stress gradient in more satisfactory manner. Also, one has to manually manipulate the mesh density in ABAQUS and other com mercial FEA software the contact patch near the contact patch, so as to obtain satisfactory and convergent solutions.

PAGE 60

60 We used 4 Node bi linear and 9 Node quadratic elements for 2D analyses. In 2D, to start off with, the tolerance value to detect contactin g elements in the implemented JAVA code was adjusted so as to reach convergence while experimenting with the increasing global density for 4 Node elements. Then using the same value of the tolerance the mesh density required for 9 Node elements was found. Obviously fewer 9N elements wi ll be required than the number of 4 Node bi linear elements. For 3D, 8 Node hexagonal elements were used. Tolerance value adjusted in the implementation code for detecting contact elements in all the examples is 0 .000 2. For e xamples 3 through 6, analytical solutions are used from specific examples in reference [20[. The co ordinate system details for benchmark examples from reference [20] and IBFEM are as follows: As seen from the benchmark examples of reference [20], at the center of the contact patch along z axis. Note that in IBFEM, y axis is equivalent to the z axis used in corresponding examples from [20 ] i.e. we have correlation for principal directions as given by Figure 5 1. Hence, which is the numerical value of contact pressure. Also which is the numerical value of shear stress. 5.1 Example 1 : Cantilever L oaded T hrough C ontact w ith a H alf C ylinder Figure 5 2 shows the 2D plane stress mode l of a cantilever beam loaded at the overhanging end through contact with a half cylinder. The cross section of the beam is considered to be rectangular. A sufficiently dense global mesh was used to obtain satisfactory solutions for maximum stresses and di splacements at the fixed end of the beam. Note that in this case, the maximum stress will always occur at the fixed support of the beam and not at the contact patch of the beam surface and the cylinder top.

PAGE 61

61 F igure 5 2 is shown with a low mesh density to ge t a clear picture of the meshed geometry We note that there is a point contact between the bottom edge of the loaded cantilever and the top edge of the cylinder transferring it a unit traction which is applied at the bo ttom edge of the half cylinder. Al so the orign point in F igure 5 2 is fixed and the left edge is free in y direction and fixed only in x direction. Analytical solutions: Cantilever loaded at the free end with a overhang and with rectan gular cross section. Z = Section modulus of the beam cr oss section = ( 5 1 ) I = Moment of inertia about x axis for the cantil ever = ( 5 2 ) Where, B = Plane stress thickness of beam cross section and length of the cylinder = 1m d = depth of the beam = 1m D = Diameter of the cylinder = 1m A = A rea of the load application on the half cylin der = D x B = 1m 2 \ W = Applied total load = 1N, hence Applied traction at the cylinder bo ttom = W/A = 1 Pa Material properties (identical to both parts in contact): F or Figure 5 3 below are the analytical expressions for stresses, displacements. 1. Stress at a specific point ( 5 3 ) 2. Stress at the support ( with constant cross section) ( 5 4 )

PAGE 62

62 3. Deflection at any point between the support and load ( 5 5 ) 4. Deflection at any point beyond load ( 5 6 ) 5. Deflection at load ( 5 7 ) 6. Maximum deflection at load ( 5 8 ) Using Equations 5 3, 5 4, 5 7 and 5 8 we have the following analytical solutions to compare with IBFEM solutions for the identical problem: 1. Stress at the support ( with constant cross section) 60 Pa 2. Deflection at load 2.0E 8 m 3. Maximum deflection at load = 2.6E 8 m 5.1.1 Stresses and D isplacements Stresses a nd displacem ents were calculated in IBFEM using 4 Node bi linear and 9 Node bi quadratic element types and compared with the analytical solution. T able 5 1 provides the maximum stress and maximum displacement magnitude values for this example with differen t mesh densities and various element types. As expected the maximum for stresses occurs at the fixed support of the cantilever beam Also we see some stress developing at the point of load application due to contact forces being transferred at that very lo cation. But the stress value here will never be higher than that at the support. We expect the maximum stress and displacement values to improve and converge by approaching towards the analytical solution as the mesh density increases.

PAGE 63

63 Also it is obvious t hat higher order element will reach the expected value with much fewer elements. Analytical values: = 6.0E1 Pa, = 2.0E 8 m, = 2.6E 8 m F igures 5 4 to 5 7 show the most satisf actory maximum stresses and displacement contours chosen out of the above experiments. The mesh densities used here are: 4 N ode bi linear elements 100 x 100 9 Node bi quadratic elements 2 0 x 2 0 This indicate s that the maximum Von Mises stress always occurs at the fixed support and there is finite contact stress at the load application point which is much less than the maximum stress at the support. Also the deflection values at the support and at the free end along with the stress values match very c losely with the analytical solutions. 9 Node elements converge much faster to the accurate solution as compared to 4 Node elements providing improved computational efficiency. Also the time taken for a single analysis with maximum of the mesh densities men tioned in T able 5 1 is less than 2 3 seconds which also explains the effectiveness of the IBFEM based implementation technique. 5.1.2 Convergence S tudies Convergence was checked for maximum Von Mises stress values and displacements ( both at load applicat ion point and at the free end) magnitudes for the loaded cantilever, by comparing the values obtained using IBF EM with the analytical solutions mentioned above F igure s 5 8 to 5 10 show the behavior of the two element types with varying mesh densities I t c an be observed that the solution value s approach towards the analytical value s as the g lobal mesh density is increased. Also 9 node elements have

PAGE 64

64 faster convergence than 4 Node elements as they can capture the stress concentrations better than the 4 Nod e elements. T his example demonstrates that the use of IBFEM based contact implementation technique provides satisfactory answers for measuring the maximum stresses away from the contact point in cases where the maximum stresses do not occur in the contact patch but somewhere else. The conditions of equilibrium are accurately maintained and followed by the solution structure used for contact implementation and the forces are transmitted accurately through the contact patch. 5.2 Example 2 : Blocks in C ontact Figure 5 11 shows the 3 D model of blocks in contact with bottom surface of the lower block fixed. Traction in negative y direction was applied on the top surface of the upper block A globally dense mesh was use d to obtain satisfactory solutions in the con tact patch. Again the F igure 5 11 is shown with a low mesh density to get a clear picture of the meshed geometry. 5.2.1 Stresses and D isplacements Stresses were calculated in IBFEM using 8 node hexagonal element type with varying mesh density and were co mpared with traditional finite element solutions obtained using ABAQU S found in reference [27 ]. Identical details of the contacting blocks used as used in [23], are as follows: Bottom block: 8m x 5m x 8m Top block: 5m x 3m x 5m Downward normal unit tracti on (1Pa ) was applied on the top surface of the top block, Analysis type: 3D stress analysis Material properties (identical to both parts in contact):

PAGE 65

65 F igure 5 12 show s the comparison of the Von Mises stress obtained using IBFEM and ABAQUS for identical geometry and material parameters of block to block contact shown above. It can be clearly seen that the maximum Von Mises stress occurs at the corner points of the contact surfaces. Also the mesh densit y used for the upper block in region of stress concentration, for obt aining the stress seen in F igure 5 12 is similar for analyses both in IBFEM and ABAQUS. The mesh densities used here are: IBFEM: 4 x11 x 4, ABAQUS: 4 x 4 x 4 The plot of Von Mises s tress conoturs for comparison purpose from ABAQUS can be seen from reference [27]. Von Mises stress value is seen to be 1.411 Pa using ABAQUS from reference [27] for the above mentioned mesh density. For identical case, IBFEM achieves 1.439 Pa with given mesh density. Note that 11 elements used in y direction for IBFEM analysis help to better capture the stress distribution in the lower block but only 4 y directional rows contribute mainly towards getting stress values in upper block where the stress conce ntration occurs. This helps establish a valid comparison between results f rom both the software packages. T able 5 2 mentions the maximum stresses for this example with different mesh densities using 3D, 8 Node Hexagonal elements. As expected the maximum f or stresses occurs very near the corner contact points on the blocks. We expect these values to improve and converge as the mesh density increases. H igher order element will reach the expected value with much fewer elements. Figures 5 13 to 5 16 use mesh d ensity 3 5 x 35 x 35 We see that all principal stresses including normal stress corresponding to maximum pressure in the contact patch are maximum on the surface at the four corner

PAGE 66

66 points with heavy discontinuity and are compres sive Also note that the location of the maximum shear stress is subsurface. This distribution of normal and shear stresses eventually leads Von Mises stress to be maximum in the subsurface region corresponding to the contact zone involved. 5.2.2 Converge nce S tudies Convergence was checked for maximum contact pressure value in the contact patch by comparing the values obtained using IBF EM with the ABAQUS FEA solution for this statically loaded contact of blocks. F igure 5 17 shows the values of maximum con tact pressure for various densities of the global mesh and 3D hexagonal, 8 Node element s It can be observed that the maximum Von Mises stress value tends to converge as the g lobal mesh density is increased. Very close convergence could not be exercised du e to computational memory limitation with higher mesh densities. T his example demonstrates that the use of IBFEM based contact implementation technique provides satisfactory answers for the maximum contact pressure in the contact patch for a block to bloc k contact. 5.3 Example 3 : Cylinder C ontacting a F lat Figure 5 18 shows the 2D model of a cylinder contacting a flat plate where o nly a quarter of cylinder and plate is modeled and symmetry boundary conditions are used. Green edges in F igure 5 18 are fixed only in x direction. Similarly, red edge is fixed only in y direction. T raction in positive y direction was applied on the bottom edge of the plate A dense mesh was use d to obtain satisfactory solutions in the contact patch Even being computationally in efficient in absence of enrichments and local mesh refinement t his helps build a fairly good picture of contact stress and displacement

PAGE 67

67 distribution in the contact zone The following F igure 5 18 is shown with a l ow mesh density to get a clear picture of the meshed geometry Analytical solutions for cylinder to cylinder contact are mentioned in section 2 2 1 in contact of cylinders. While dealing with cylinder to flat contact, we make the radius of curvature for the flat plate to be infinity considering i t as a cylinder with infinite radius to simulate flatness of the plate. 5.3.1 Stresses and D isplacements Stresses and displacements were calculated in IBFEM using 4 Node and 9 Node element types The stresses were compared with analytical sol utions obtain ed from Hertz theory useds in example 7 2 from referece [20] (Ch.7, page 443) for statically loaded con tact of cylindrical wheel with a flat rail. F rom this reference, the details of the conta c ting bodies used are as follows: Cylindrical wheel radius, R1 = 6 inch = 0.1524 m Radius of the flat rail, R2 = Radial total load on the wheel = P = 5000 lb = 22241 N Plane strain analysis Length of the cylinder and plate perpendicular to paper = 0.875 inch = 0.022225 m The traction t o be applied on bottom edge of the plate is calculated as below: Load on half symmetry model of plate = 5000/2 = 2500 lb = 11120.55 N Area of the half symmetry plate loaded = (width in x direction of the plate) x ( length of the plate) = 5.25 in 2 = 0.00 338709 m 2 Traction value on bottom edge of the plate = 2500/5.25 = 476.2 psi = 3283283.41 Pa Material properties (identical to both parts in contact): E = 3.0E07 psi = 2.07E11 Pa Using above data with eq u a tions dev eloped using Hertz theory of elastic contact: The analytical solutions for this case are as follows:

PAGE 68

68 = 70243 psi = 4 84 3E8 Pa = psi = 2.712E8 Pa 21354 psi = 1.472E8 Pa at a depth of 0.041 in = 0.00104 m subsurface. Table 5 3 mention s the maximum stress and maximum displacement magnitude values for this example with different mesh densities and various element types We expect these values to improv e and converge towards the analytical solution as the mesh density increases. As mentioned at the beginning of this chapter, we first adjust the tolerance for detecting contacting elements using convergence study for 4 Node elements and then find the requi red density for 9 Node elements to achieve near y the same value for maximum contact stresses. Analytical values: = 4 84 3E8 Pa = 1.472E8 Pa Figures 5 19 to 5 24 show the most satisfactory contact str esses and displacement contour s chosen out of the above experiments We see that all principal stresses including normal stress corresponding to maximum pressure in the contact patch are maximum on the surface at y=0 and are comp ressive in nature. Also note that the location of the maximum shear stress is subsurface, approximately very close to the depth value given by the analytical so lution This distribution of normal and shear stresses eventually leads Von Mises stress to be m aximum in the subsurface zone. T he mesh densities used here are : 4 Node bi linear elements 5 00 x 5 00 9 Node bi quadratic elements 8 0 x 8 0 In this analysis, we noted that the number of contacting elements with the maximum mesh density mentioned in th e T able 5 3 for the two element types are as follows These approximately model the contact patch to the analytical value of the patch half width calculated. 4 Node : 9 contact elements 9 Node : 2 contact elements

PAGE 69

69 5.3.2 Convergence S tudies Convergence wa s checked for maximum contact pressure value in the contact patch by comparing the values obtained using IBFEM with the analytical solution from example 7 2 of reference [20] (Ch.7, page 443) for statically loaded contact of cylindrical wheels with a fla t rail. F igure 5 25 shows the values of maximum contact pressure for various densities of the global mesh and various element types It can be observed that the maximum contact pressure value approaches towards the analytical value as the global mesh den sity is increased and this progress is faster with increase in order of the element type used. Th is example demonstrates the use of IBFEM based contact implementation technique gives satisfactory answers for the maximum contact stresses in the contact pa tch for a cylinder to flat contact. 5.4 Example 4 : Contact of C ylinders Figure 5 26 shows the 2D model of cylinders in contact where only quarters of cylinders are modeled and symmetry boundary conditions are used. Green edges fixed only in x direction. Re d edges fixed only in y direction. Traction in positive y direction was applied on the bot tom edge of the bottom cylinder A dense mesh was used to obtain satisfactory solutions in the contact patch. T his example is also used to study contact stress and di splacement distribution in the contact zone when there exists a line contact Note that the F igure 5 26 is shown with a low mesh density to get a clear picture of the meshed geometry. Analytical solutions for cylinder to cylinder contact are mentioned in section 2 2 1 in contact of cylinders. Here we try to compare the stresses and displacement obtained

PAGE 70

70 from IBFEM based contact problem analysis with the analytical solution mentioned in section 2 2 1 using a specific example from reference [20]. 5.4.1 Stres ses and D isplacements Stresses and displacements were calculated in IBFEM using 4Node and 9Node element types and were compared with analytical solutions obtained using Hertz theory presented in example 7 2 of reference [20] (Ch.7, page 443) for staticall y loaded contact of cylinders. From this reference, the details of the contacting bodies used are as follows: C ylinder radi i = R1 = 6 inch = 0.1524m Plane strain analysis Radial total load on the cylinder = P = 5000 lb = 22241 N, Length of the cylinders perpendicular to paper = 0.875 inch = 0.022225 m The traction to be applied on bottom edge of the bottom cylinder is calculated as below: Load on half symmetry model of cylinders = 5000/2 = 2500 lb = 11120.55 N Area of the half symmetry model loaded = ( Ra dius of the cylinder ) x ( length of the cylinder ) = 5.25 in 2 = 0.00338709 m 2 Traction value on bottom edge of the plate = 2500/5.25 = 476.2 psi = 3283283.41 Pa Material properties (identical to both parts in contact): .07E11 Pa Using above data with equations developed using Hertz theory of elastic contact: The analytical solutions for this case are as follows: = 99450 psi = 6.856E8 Pa = 55692 psi = 3.839 E8 Pa Also analytical, 30233 psi = 2.084E8 Pa at a depth of 0.000728 m subsurface T able 5 4 mentions the maximum stress and displacement magnitude s

PAGE 71

71 Analytical values: = 6.856E8 Pa = 2.0841E8 Pa F igures 5 27 to 5 32 show the most satisfactory contact stresses and displacement contours chosen out of the above experiments. We see that all principal stresses including normal stress corresponding to maximum pressure in the contact patch are maximum on the surface at y=0 and are compressive in nature. Also note that the location of the maximum shear stress is subsurface, approximately very close to the depth va lue given by the analytical so lution This distribution of normal and shear stresses eventually leads Von Mises stress to be maximum in the subsurface zone. The mesh densities used here are: 4 Node bi linear elements 450 x 405 9 Node bi quadratic eleme nts 95 x 95 In this analysis, we noted that the number of contacting elements with the maximum mesh density mentioned in the T able 5 4 for the two element types are as follows These approximately model the contact patch to the analytical value of the pa tch half width calculated. 4 Node : 4 contact elements 9Node : 1 element 5.4.2 Convergence S tudies Convergence was checked for maximum contact pressure value in the contact patch by comparing the values obtained using IBFEM with the analytical solution from example 7 2 of reference [20] (Ch.7, page 443) modified for statically loaded contact of cylinders. F igure 5 33 shows the values of maximum contact pressure for various densities of the global mesh and various element types It can be observed that the maximum contact pressure value approaches towards the analytical value as the global mesh

PAGE 72

72 density is increased and this progress is faster with increase in order of the element type used T his example demonstrates the use of IBFEM based contact implem entation technique to get satisfactory answers for the maximum contact pressure in the contact patch for a cylindrical contact. 5.5 Example 5 : Sphere C ontacting a F lat Figure 5 34 shows the 3D model of a sphere contacting a flat plate whe re only a 1/8 th of sphere and quarter of the plate is modeled and symmetry boundary conditions are used. The left faces of both sphere and plate and also the right face of the plate are fixed only in x direction. Top face of the sphere is fixed only in y d irection. Back face of the sphere is fixed only in z direction. Traction in positive y direction was applied on the bottom surface of the plate. A dense mesh was used to obtain satisfactory solutions in the contact patch. 5.5.1 Stresses and D isplacements Stresses and displacements were calculated in IBFEM using 8 Node Hexagonal elements and were compared with analytical solutions obtained using Hertz theory presented in example 7 1 of reference [20] (Ch.7, page 439) for statically loaded contact of spheri cal ball with a flat race. From this reference, the details of the contacting bodies used are as follows: Sphere radius, R1 = 0.197 inch = 0.0050038 m Radius of the flat race, R2 = Radial total load on the sphere = P = 21.5 lb = 95.63 N 3D stress analysis Length of the plate perpendicular to paper = 0.197 inch = 0.0050038 m Width of the plate in x direction = 0.197 inch = 0.0050038 m The traction to be applied on bottom edge of the bottom cylinder is calculated as below:

PAGE 73

73 Load o n quarter symmetry plate = 21.50/4 = 5.375 lb = 23.9075 N Area of the half symmetry plate loaded = ( Width of the plate in x direction) x (length of the plate) = 0.197 x 0.197= 0.038809 in 2 = 0.0000250 m 2 Traction value on bottom edge of the bottom plate = 5.375/ 0.038809 = 138.5 psi = 954923.88 Pa Material properties (identical to both parts in contact): = 3.0E07 psi = 2.07E11 Pa Using above data with equations developed using Hertz theory of elastic contact: Th e analytical solutions for this case are as follows: = 305381psi = 2.105 E 9 Pa = 238197 psi = 1.642 E 9 Pa Also, 103083 psi = 7.107 E6 Pa at a depth of 0.00009 4 m subsurface and 44789 psi = 3.088E 6 Pa at edge of contact patch with tensile stress of same magnitude T able 5 5 mentions the maximum stress and maximum displacement magnitude values for this example with differen t mesh densities and I3DHexa8N ( a 3D 8 node hexagonal) element We expect these values to improve and converge towards the analytical solution as the mesh density increases. As mentioned at the beginning of this chapter, we first adjust the tolerance for detecting contacting elements using convergence study for 8Node elements Analytical values: = 2.105E9 Pa, = 7.107 E6 Pa Figures 5 35 to 5 37 show the most satisfactory contact stresses and displaceme nt contours chosen out of the above experiments. We see that all principal

PAGE 74

74 stresses including normal stress corresponding to maximum pressure in the contact patch are maximum on the surface at y=0 and are compressive in nature. A lso note that the location of the maximum shear stress is subsurface, approximately very close to the depth value given by the analytical so lution This distribution of normal and shear stresses eventually leads Von Mises stress to be maximum in the subsur face zone. The mesh densities used here are: 8 Node hexagonal elements 275 x 275 In the above analysis, we noted that only a single contacting element is seen to be present when the maximum mesh density mentioned in the Table 5 5 for the 3D 8Node hexag onal element type is used. This approximately models the contact patch to the analytical value of the patch half width calculated. 5.5.2 Convergence S tudies Convergence was checked for maximum contact pressure value in the contact patch by comparing the values obtained using IBFEM with the analytical solution from example 7 1 of reference [20] (Ch.7, page 439) for statically loaded contact of sphere and flat race. F igure 5 38 shows the values of maximum contact pressure for various densities of the glob al mesh and 3D 8Node hexagonal element type It can be observed that the maximum contact pressure value approaches towards the analytical value as the global mesh density is increased and this progress is expected to be faster with increase in order of th e element t ype. T his example demonstrates the use of IBFEM based contact implementation technique to get satisfactory answers for the maximum contact pressure in the contact patch for a sphere to flat contact.

PAGE 75

75 5.6 Example 6 : Contact of S pheres Figure 5 3 9 shows the 3D model of a sphere contacting a sphere where only 1/8 th of the spheres were modeled and symmetry boundary conditions are used. Traction in positive y direction was applied on the bott om surface of the bottom sphere Left faces of the spheres are fixed only in x direction. Back faces are fixed only in z direction. Top face of the top sphere is fixed only in y direction. Analytical solutions for sphere to sphere contac t are mentioned in section 2 3 1 in contact of spheres. 5.6.1 Stresses and D isplacements Stresses and displacements were calculated in IBFEM using 8 Node Hexagonal elements and were compared with analytical solutions obtained using Hertz theory presented in example 7 1 of reference [20] (Ch.7, page 439) modified for statica lly l oaded contact of spheres From this reference, the details of the contacting bodies used are as follows: Sphere radii, R1 = 0.197 inch = 0.0050038 m 3D stress analysis Radial total load on the sphere = P = 21.5 lb = 95.63 N, The traction to be applied o n bottom surface of the bottom sphere is: Load on 1/8 th symmetry model = 21.50/4 = 5 .375 lb = 23.9075 N Area of the 1/8 th symmetry sphere loaded = (1/4) x (cross sectional area of the sphere at half symmetry along y axis) = 0.25 x P i x 0.197 2 = 0.038 8 0 in ch 2 = 0.000025 m 2 Traction value on bottom surface of the bottom sphere = 5.375/ 0.03048 = 1 76.34 psi = 1215821.23 Pa Material properties ( Same for both parts ): E = 2.07E11 Pa

PAGE 76

76 The analytical solutions for this case are as follows: = 484715 psi = 3.341E9 Pa = 378077 psi = 2.606 E9 Pa Also, 163601 psi = 1.127E6 Pa at a depth of 0.0007366 m subsurface and 71092 psi = 4.9E8 Pa at the edge of contact patch with tensile stress of same magnitude Ta ble 5 6 mentions the maximum stress and maximum displacement magnitude values wit h different mesh densities and 3D 8 node hexagonal element. Analytical va lues: = 3.341E 9 Pa, = 1.127E 6 Pa Figures 5 40 to 5 42 show the most satisfactory contact stresses and displacement contours chosen out of the above experiments. We see that all principal stresses inc luding normal stress corresponding to maximum pressure in the contact patch are maximum on the surface at y=0 and are compressive in nature. Also note that the location of the maximum shear stress is subsurface, approximately ver y close to the depth value given by the analytical so lution This distribution of normal and shear stresses eventually leads Von Mises stress to be maximum in the subsurface zone. The mesh densities used here are: 8 Node hexagonal e le ments 275 x 275 In this analysis, we noted that only single contacting element is seen to be present when the maximum mes h density mentioned in the T able 5 6 for the 3D 8Node hexagonal element type is used. This approximately models the contact patch to the analytical value of the patch half width calculated. 5.6.2 Convergence S tudies Convergence was checked for maximum contact pressure value in the contact patch by comparing the values obtained using IBFEM with the analytical solution from

PAGE 77

77 example 7 1 of reference [20] (Ch .7, page 439) modified for statically loaded contact of sphere and flat race. Figure 5 43 shows the values of maximum contact pressure for various densities of the global mesh and 3D 8Node hexagonal element type It can be observed that the maximum cont act pressure value approaches towards the analytical value as the global mesh density is increased and this progress is expected to be faster with increase in order of the element t ype. Also the variation of the maximum contact pressure with the total loa d was compared with variation of analytical solution with applied total load. The analytical solution for maximum contact pressure in contact of spheres is given in reference [27],which is: A plot of contact pressure or Sigma YY versus total applied load was obtained as shown in F igure 5 44, using 300 x 300 x 1 mesh. This plot also help to see the close match bwtween solutions obtained using IBFEM and the analytical solution. T his example demonstrates that the use of IBFEM bas ed contact implementation technique yields satisfactory answers for measuring the maximum contact pressure in the contact patch for a sphere to sphere contact.

PAGE 78

78 Table 5 1. Maximum stress (Von Mises) and displacement magnitudes for Example 1 4Node Mesh IBFEM (Pa) IBFEM (m) IBFEM (m) 9 N ode Mesh IBFEM (Pa) IBFEM (m) IBFEM (m) 30 x 30 5.703 E1 1. 7 E 8 3.207 E 8 5 x 5 5.649 E1 2.00 E 8 3.502 E 8 50 x 50 5.899 E1 1.9 E 8 3.401 E 8 10 x 10 5.974 E1 2.01 E 8 3.519 E 8 70 x 70 5.963 E1 1.93 E 8 3.459 E 8 20 x 20 6.062 E1 2.02 E 8 3.521 E 8 100 x 100 6.002 E1 1.999 E 8 3.49 E 8 30 x 30 6.066 E1 2.02 E 8 3.521 E 8 150 x 150 6.024 E1 2.00 E 8 3.507 E 8 40 x 40 6.068 E1 2.02 E 8 3.521 E 8 Table 5 2 Contact region stress and displacement for Example 2 8Node Mesh IBFEM (Pa) IBFEM (Pa) IBFEM (m) 10 x 10 x 10 1.305 1.428 2.716E 7 20 x 20 x 20 1.934 2.126 2.834E 7 30 x 30 x 30 2.420 2.985 2.871E 7 35 x 35 x 35 2.977 3.850 2.883E 7 Table 5 3. Maximum contact region stresses for example 3 4Node Mesh IBFEM (Pa) IBFEM (Pa) 9 N ode Mesh IBFEM (Pa) IBFEM (Pa) 100 x 100 1 .90 E8 4.59 E7 20 x 20 1.35 E8 1.61 E7 200 x 200 2.89 E8 6.14 E7 30 x 30 2.00 E8 2.41 E7 400 x 400 4.60 E8 7.68 E7 40 x 40 2.66 E8 3.20 E7 450 x 450 4.71 E8 7.79 E7 50 x 50 3.31 E8 3.99 E7 460 x 460 4.82 E8 7.84 E7 60 x 60 3.54 E8 3.16 E7 475 x 475 4.85 E8 7.95 E7 70 x 70 4.62 E8 3.69 E7 500 x 500 4.89 E8 8.02 E7 80 x 80 4.70 E 8 4.21 E7

PAGE 79

79 Table 5 4 Maximum contact region stresses for Example 4 4Node Mesh IBFEM (Pa) IBFEM (Pa) 9 N ode Mesh IBFEM (Pa) IBFEM (Pa) 100 x 1 00 2.61 E8 6.788 E7 20 x 20 1.53 E8 1.65 E7 200 x 200 3.99 E8 9.42 E7 30 x 30 2.26 E8 2.46 E7 300 x 300 5.14 E8 1.10 E8 40 x 40 3.00 E8 3.27 E7 400 x 400 6.21 E8 1.16 E8 60 x 60 4.47 E8 4.90 E7 430 x 430 6.81 E8 1.29 E8 80 x 80 5.93 E8 6.52 E 7 450 x 450 6.92 E8 1.35 E8 95 x 95 6.88 E8 7.41 E7 Table 5 5 Maximum contact region stresses for Example 5 8Node Mesh IBFEM (Pa) IBFEM (Pa) 100 x 100 4.992 E7 1.648 E6 150 x 150 8.111 E8 4.1 32 E6 200 x 200 1.289 E9 6.369 E6 225 x 225 1.563 E9 7.29 E6 250 x 250 1.884 E9 8.45 E6 260 x 260 2.119 E9 8.891 E6 275 x 275 2.218 E9 9.514 E6

PAGE 80

80 Table 5 6 Maximum contact region stresses for Example 6 8Node Mesh IBFEM (Pa) IBFEM (Pa) 175 x 175 1.513 E9 4.303 E6 190 x 190 1.76 E9 4.506 E6 200 x 200 1.921 E9 4.62 E6 250 x 250 2.84 E9 4.643 E6 260 x 260 3.123 E9 4.673 E6 275 x 275 3.345 E9 3.743 E6 280 x 280 3.395 E9 3.751 E6 Figure 5 1 Comparison: Co ordinate system for reference [20] and IBFEM Figure 5 2. Example 1: 2D Plane stress model of a cantilever loaded through contact

PAGE 81

81 Figure 5 3. Example 1: Schematic of a cantilever loaded w ith an overhang Figure 5 4. Von Mises stress contours: E xample with 4 node bi linear elements Figure 5 5. Displacement magnitude contours: Example 1 with 4 node bi linear elements W=1N x y L =10m b=2m W WL

PAGE 82

82 Figure 5 6. Stress contours: E xample 1 with 9 node bi quadratic eleme nts Figure 5 7. Displacement contours: Example 1 with 9 node bi quadratic elements Figure 5 8. Convergence of maximum Von Mises stress

PAGE 83

83 Figure 5 9. Convergence of maximum displacement at loading point Figure 5 10. Convergence of maximum displacemen t at the free end

PAGE 84

84 Figure 5 11. Example 2: A 3D contact of blocks Figure 5 12 Von Mises Stress contours with 3D elements

PAGE 85

85 Figure 5 13. Von Mises stress contours: Example 2 with 3D 8Node elements Figure 5 14 Normal stress contours: Example 2 with 3D 8Node elements

PAGE 86

86 Figure 5 15 Maximum shears stress contours: Example 8 with 3D 8Node elements Figure 5 16 Displacement magnitude contours: E xample 2 with 3D 8Node elements

PAGE 87

87 Figure 5 17. Conve rgence of maximum Von Mises stress in contact of blocks using 8Node elements Figure 5 18. Example 3: 2D Quarter symmetry contact of cylinder and plate

PAGE 88

88 Figure 5 19. Stress contours: E xample 3 with 4 node bi linear elements a. Von Mises stres s distribution, b. Zoomed in View Figure 5 20. Stress contours: Example 3 with 4 node bi linear elements a. Normal stress b. Shear stress

PAGE 89

89 Figure 5 21. Displacement contours: Exampl e 3 with 4 node bi linear elements a. magnitude distribution, b. Zoomed in view Figure 5 22. Stress contours: E xample 3 with 9 node bi quadratic elements a. Von Mises stress distribution, b. Zoomed in View

PAGE 90

90 Figure 5 23. Stress contours: Example 3 using IBFEM with 9 node bi quadratic elements a. Normal stress b. Shear stress Figure 5 24. Displacement contours: Example 3 using IBFEM with 9 node bi quadratic elements a. mag nitude distribution, b. Zoomed in view

PAGE 91

91 Figure 5 25. Convergence for maximum contact pressure in cylinder flat contact Figure 5 26. Example 5: 2D Quarter symmetry contact of cylinder s

PAGE 92

92 Figure 5 27 Stress contours: E xample 4 with 4 node bi linear elements a. Von Mises stress distribution, b. Zoomed in view Figure 5 28. Stress contours: E xample 4 with 4 node bi linear elements a. Normal stress b. Shear stress

PAGE 93

93 Figure 5 29 Displacement contours: Example 4 with 4 node bi linear elements, a. magnitude distribution, b. Zoomed in view Figure 5 30. Stress contours: Example 4 with 9 node bi quadratic elements, a. Von Mises stress distribution, b. Zoomed in vi ew

PAGE 94

94 Figure 5 31. Stress contours: Example 4 with 9 node bi quadratic elements a. Normal stress b. Shear stress Figure 5 32. Displacement contours: Example 4 with 9 node bi quadrat ic elements

PAGE 95

95 Figure 5 33. Convergence study for maximum contact pressure in cylindrical contact Figure 5 34. Example 5 : a. A 3D symmetrical contact of sphe r e and plate

PAGE 96

96 Figure 5 35. Stress contours: Example 6 using IBFEM with 8 node hexagonal elements a. Von Mises stress distribution, b. Zoomed in view Figure 5 36. Stress contours: Example 5 using IBFEM with 8 node hexagonal elements a. Normal stress b. Shear stress

PAGE 97

97 Figu re 5 37. Displacement contours: E xample 5 with 8 node hexagonal elements a. magnitude distribution, b. Zoomed in view Figure 5 38. Convergence study for maximum contact pressure in sphere to plate contact

PAGE 98

98 Figure 5 39. Example 7 : A 3D symmetrical co ntact of sphe r es Figure 5 40. Stress contours: Example 6 with 8 node hexagonal elements a. Von Mises stress distribution, b. Zoomed in view

PAGE 99

9 9 Figure 5 41. Stress contours: Example 6 with 8 node hexagonal elements a. Normal stress b. Shear stress Figure 5 42. Displacement contours: E xample 6 with 8 node hexagonal elements a. magnitude distribution, b. Zoomed in view

PAGE 100

100 Figure 5 43. Convergence stud y for maximum contact pressure in spherical contact Figure 5 44. Variation of maximum contact pressure with applied total load in spherical contact

PAGE 101

101 CHAPTER 6 CONCLUSION AND FUTURE SCOPE 6.1 S ummary Most of the traditional FEA packages represent th e geometry by means of finite element mesh which acts as the object to be analyzed bringing in initial geometric inaccuracies IBFEM uses sets of mathematical equations extracted from geometry files generated using standard CAD software such as Pro Enginee r or Solid Works to represent the geometry This minimizes initial ina ccuracies in geometric modeling IBFEM g enerat es a uniform structured background mesh that encloses the geometry A structured mesh is easier to generate and re quires minimal efforts for regeneration or modification than a conforming mesh which is a major reason for heavy time consumption in contact analysis using traditional FEA software like ABAQUS. F urthermore in IBFEM all the elements used in the analysis are regular and undistorted unlike traditional FEM and hence problems associated with distorted elements are eliminated. As shown through several examples in this work models with much fewer higher order elements can be used for getting satisfactory solutions. T he implicit equation s of the essential boundaries are used to impose essential boundary conditions which eliminate the need of having nodes on the boundaries. Another advantage of this method is that all internal elements have identical stiffness matrix resulting in significa nt reduction in computation required to determine the global stiffness matrix. The primary motivation behind the current work is to be able to model tied contact problems using IBFEM It wa s desire d to obtain satisfactory solutions in the region of contac t stress concentration using structured mesh with minimized pre processing efforts Attempt was made to establish and maintain the state of static equilibrium along

PAGE 102

102 with continuity of displacement at the interface in a two part contact assembly model. The assembly model consists of two separate meshes, each one corresponding to the parts in contact. These meshes overlap in the region where the contact occurs between the two objects. A solution structure constructed previously by Burla et al. [13] using appr oximate step function for treatment of material discontinuity was use d here for test and trial function s. The solution structure is such that it ensures fulfillment of interface conditions on the contacting boundaries of the two parts in contact by blendin g the solutions from the two meshes in the region of contact. The contact region is defined using contact elements located in the overlapping region of the two meshes. C ontact element s are composed of two original structured elements, one from each mesh. T h ese contact element s bel ong to both the parts. The displacement in a contact element is a blend of solutions from corresponding original elements in the individual meshes. This makes possible displacement continuity at the interface. The stiffness of cont act elements gets determined using special technique as explained in detail in reference [13]. Displacements in the individual parts are governed by elements from corresponding individual meshes This solution structure was implemented both in two and thre e dimensions to detect the point and line contact. Note that c urrent work uses very dense mesh to model the contact patch which is computationally inefficient L ocalized mesh refinement or enrichment technique s, such as XFEM are needed Nonetheless, the convergence analyses performed show that the solution structure implemented for contact has good convergence capabilities for various element types and provides fairly good results It is able to captur e the contact stress es while modeling the contact betw een two bodies

PAGE 103

103 satisfactorily indicating that reasonable solutions can be obtain ed for tied contact problems using IBFEM with minimal pre processing efforts. 6.2 Scope of Future Work 1. Localized mesh refinement will allow contact analysis with fewer elements and nodes thus reducing the computational cost. 2. Ideally adaptive mesh refinement should be used in order to reduce the manual intervention in mesh generation and refinement. The h p and k adaptivity can be incorporated. 3. Another appr oach for locally improving the quality of the solution is to use XFEM approach which locally enriches the FEM solution using terms from the analytical solution. 4. The boundaries of the contacting bodies should be more accurately approximated near the contact when you have point or line contact. In the current implementation the geometry is approximated using line segments or triangles that are sufficiently small so that the geometry looks smooth but near the contact the geometry must be approximated even more accurately. 5. Contact is inherently a non linear problem. To achieve even more accurate and converging solutions along with improved stress and displacement distribution in the contact patch, non linear contact analysis needs to be implemented.

PAGE 104

104 LIST OF RE FERENCES [1] K L Johnson, Contact Mechanics 1985 [2] Hills, D. A., Mechanics of Fretting Fatigue [3] O C Zienkiewicz and R L Taylor, Finite Element Method for Solid and Structural Mechanics, 6th Edition [4] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, Meshless methods: An overview and recent developments," Computer Methods in Applied Mechanics and Engine ering, vol. 139, pp. 3 47, 1996 [5] Owen DRJ, Fawkes AJ, Engineering Fracture Mechanics: Numerical methods and application. Pineridge: Swansea 1983 [6] T. Be lytschko, T Black. Elastic crack growth in finite elements with minimal remeshing. Internat. J. Numer. Methods.Engrg. 45: 601 620, 1999 [7] T. Belytschko, C. Parimi, N. Moes, N. Sukumar, and S. Usui, "Structured extended finite element methods for solids defin ed by implicit surfaces," International Journal for Numerical Methods in Engineering vol. 56, pp. 609 635 [8] V. Shapiro and I. Tsukanov, Mesh free simulation of deforming domains," Computer Aided Design vol. 31, pp. 459 471, 1999 [9] K. Hollig, C. Apprich, and A. Streit, "Introduction to the Web method and its applications," Advances in Computational Mathematics vol. 23, pp. 215 237, 2005. 2003 [10] A. V. Kumar, R. Burla, S. Padmanabhan, and L. X. Gu, "Finite element analysis using nonconforming mesh," Journal of Co mputing and Information Science in Engineering, vol. 8, 2008 [11] A. V. Kumar, S. Padmanabhan, and R. Burla, "Implicit boundary method for finite element analysis using non conforming mesh or grid," International Journal for Numerical Methods in Engineerin g, vo l. 74, pp. 1421 1447, 2008 [12] R. K. Burla and A. V. Kumar, "Implicit boundary method for analysis using uniform B spline basis and structured grid," International Journal for Numerical Methods in Engineerin g, vol. 76, pp. 1993 2028, 2008 [13] R. K. Burla, A. V. Ku mar, and B. V. Sankar, "Implicit boundary method for determination of effective properties of composite microstructures," International Journal of Solids and Structure s, vol. 46, pp. 2514 2526, 2009

PAGE 105

105 [14] P. D. S. Periyasamy, "Finite element analysis of shell li ke structures using implicit boundary method," in Mechanical and Aerospace engineering, vol. Master of Science. Gainesville: Univ ersity of Florida, 2009, pp. 86 [15] Osher SJ, Fedkiw RP (2002) Level Set Methods and Dynamic Implicit Surfaces Springer Verlag, Ph iladelphia [16] Hllig K (2003) Finite element methods wit h B Splines. SIAM, Philadelphia [17] Burla RK, Kumar AV (2008) Implicit boundary method for analysis using uniform B spline basis and structured grid. International Journal for Numerical Methods in Engineerin g. 76(13): 19 93 2028. DOI: 10.1002/nme.2390 [18] I Babuska and J M Melenk, A partition of unity method, Internat. J. of Numer. Methods. in Engrg., 40: 727 758. 1997 [19] Tsukanov I, Shapiro V (2007) Adaptive multiresolution refinement with distance fields. Internat ional Journal for Numerical Methods in Engineering 72: 1355 1386 [20] Robert L Norton, 2006, Machine Design An I ntegrated Approach 3 rd Edition [21] H Hertz. On the Contact of Elastic Solids. J. Math., 92: pp. 156 171. 1881 (in German) [22] H. Hertz. Contact of Elastic So lids, in Miscellaneous Papers, P. Lenard, ed. Macmillan and Co. Ltd.: London, pp. 146 162, 1896 [23] H.L. Whittemore and S. N. Petrenko. Friction and Carrying Capacity of Ball and Roller Bearings, Technical Paper 201. National Bureau of Standards. Washinton, D. C., 1921 [24] H. R. Thomas and V.A. Horesch, Stresses due to the Pressure of One Elastic Solid on Another, Bulletin 212, U. Illinois Engineering Experiment Station, Champaign, IL, July 15, 1930 [25] E.I. Radzimovasky, Stress Distribution and Strength Condition of Tw o Rolling Cylinders, Bulletin 408, U. Illinois Engineering Experiment Station, Champaign, IL,Feb 1953 [26] J. O. Smith and C. K. Lui, Stresses due to Tangential and Normal Loads on Elastic Solids with Application to some Contact Stress Problems J. Appli. Mech Trans. ASME, 75:pp. 157 166, 1953 [27] Harry Harkness, George Ang, Prashanth Vijalpura and Dan Cozokaru Contact stress accuracy with robust and broadly applicable implicit contact algorithm NWC ( Nafems World Congress), 2011

PAGE 106

106 BIOGRAPHICAL SKETCH Anand Tups akhare was born in Aurangabad, India which is in state of Ma harashtra. He graduated with a b in m echanical e ngineering from College of Engineering, Pune, an autonomous institute of Government of Maharashtra in May 2008 He then worked as a design and execution engineer in the Department of Marine Equipments in Heavy Engineering division of Larsen and Toubro Limited, in Powai, Mumbai from July 2008 to July 2010. In f all 2010 he joined the m program at the Department of Mechani cal and Aerospace Engineering of University of F lorida, Gainesville, USA and graduated with the in mechanical engineering in spring 2012. His areas of interests are finite element analysis, contact mechanics and CAD