UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 IMPLEMENTATION OF EXTENDED FINITE ELEMENT METHOD USING IMPLICIT BOUNDARY APPROACH By VISHAL HUNDRAJ HOTWANI A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DE GREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2011 PAGE 2 2 2011 Vishal Hundraj Hotwani PAGE 3 3 To my parents Deepa and Hundraj Hotwani, my sister Vinita Hotwani, Kamal Thadani and my advisor Dr. Ashok Kumar. PAGE 4 4 ACKNOWLEDGMENTS I want to tha nk my advisor Dr. Ashok Kumar for his support as it has only been possible because of his guidance and help. I am grateful that he was available at every point and was ready to offer his time and helpful suggestions. It was a great learning experience work ing with him and I am sure it will help me in every aspect of my life. I would also like to thank Dr. Nagaraj Arakere and Dr. Nam Ho Kim for their willingness to help me. I am grateful for their useful suggestion and inputs that helped me finish this disse rtation. I would also like to take this chance to thank all my family members who have helped me in every aspect of my life especially my grandparents Hariram, Rani, Madhavdas and Dhani. I also appreciate the support of my friends here at University of Flo rid a especially Natasha, Lucas, Kevin Taiki and Matthew PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 8 ABSTRACT ................................ ................................ ................................ ................... 10 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 12 Overview ................................ ................................ ................................ ................. 12 Goals and Objectives ................................ ................................ .............................. 16 Outline ................................ ................................ ................................ .................... 17 2 EXTENTE D FINITE ELEMENT ANALYSIS ................................ ............................ 18 Introduction ................................ ................................ ................................ ............. 18 General Formulation ................................ ................................ ............................... 21 Blending ................................ ................................ ................................ .................. 24 Application to Fracture Mechanics ................................ ................................ .......... 26 3 MESH INDEPENDENT FINTE ELEMENT ANALYSIS ................................ ........... 30 Developments in Mesh Independent and Meshless Methods ................................ 30 Implicit Boundary Finite Element Method ................................ ................................ 32 Gr id Generation and Application of EBC ................................ ................................ 33 4 IMPLEMENTATION OF MESH INDEPENDENT XFEM ................................ ......... 36 Ramped Heaviside Function ................................ ................................ ................... 36 Singularity Enrichment ................................ ................................ ............................ 38 Fixing of Nodes ................................ ................................ ................................ ....... 40 Derivation of Stiffness Matrix ................................ ................................ .................. 42 Implementation Scheme for Multiple Enrichments ................................ .................. 45 Integration of Enriched Elements ................................ ................................ ............ 48 Algorithm used for Implementation ................................ ................................ ......... 51 5 RESULTS ................................ ................................ ................................ ............... 55 Edge Crack Under Mixed Mode Loading ................................ ................................ 55 Convergence Study ................................ ................................ .......................... 61 SIF Computations ................................ ................................ ............................. 65 PAGE 6 6 Effect of Number of Enrichment Layers ................................ ............................ 66 Center Crack Under Uniform Far Field Mode I Loading ................................ .......... 68 Convergence Study ................................ ................................ .......................... 72 SI F Computation ................................ ................................ .............................. 74 Effect of Number of Enrichment Layers ................................ ............................ 76 6 Conclusion ................................ ................................ ................................ .............. 78 Summary ................................ ................................ ................................ ................ 78 Scope of Future Work ................................ ................................ ............................. 80 LIST OF REFERENCES ................................ ................................ ............................... 82 BIO GRAPHICAL SKETCH ................................ ................................ ............................ 86 PAGE 7 7 LIST OF TABLES Table page 5 1 Legend for enrichment functions. ................................ ................................ ....... 57 5 2 Enrichment scheme for cases 1 to 9 (Problem 1) ................................ ............... 57 5 3 Enrichment scheme for cases 10 to 17 (Problem 1) ................................ ........... 57 5 4 Enrichment s cheme for cases 18 to 22 (Problem 1) ................................ ........... 57 5 5 Maximum displacement values for case 1 to 4 (Problem 1) ............................... 58 5 6 Maximum displacement values for case 5 to 9 (Problem 1) ............................... 58 5 7 Maximum displacement values for case 10 to 13 (Problem 1) ........................... 59 5 8 Maximum displacement values for case 14 to 18 (Problem 1) ........................... 59 5 9 Maximum displacement values for case 19 to 22 (Problem 1) ........................... 59 5 10 SIF values for cas e 15 to 18 (Problem 1) ................................ ........................... 60 5 11 SIF values as a function of enrichment terms (Problem 1) ................................ 65 5 12 Effect of radius of enrichment ( Problem 1) ................................ ......................... 67 5 13 En richment scheme for different cases (problem 2) ................................ ........... 70 5 14 Maximum displacement values (problem 2) ................................ ....................... 70 5 15 Maximum displacement values (problem 2) ................................ ....................... 71 5 16 SIF (KI) values for case 5 and 7(problem 2) ................................ ....................... 71 5 17 SIF as a function of number of enrichment terms (problem 2) ............................ 75 5 18 Effect of number of enrichment layers (problem 2) ................................ ............. 76 PAGE 8 8 LIS T OF FIGURES Figure page 2 1 Demarcation of enriched region and blending region. ................................ ........ 19 2 2 Enrichment technique for model ing cracked domain. ................................ ......... 29 3 1 Model for meshless method showing nodes and boundary ................................ 30 3 2 Model representing mesh independent metho d ................................ ................. 32 4 1 Heaviside function for strong discontinuity along Y=0 ................................ ........ 37 4 2 Ramped step function with discontinuity along Y=0 ................................ ........... 37 4 3 Enrichment functions used for direct computation of SIF ................................ ... 39 4 4 Blending of Heaviside function using fixing of nodes. ................................ ......... 41 4 5 Blending singularity function using fixing of nodes. ................................ ............ 41 4 6 Chart showing inheritance of classes and formation of elements. ...................... 46 4 7 Chart showing inheritance of classes and formation of enriched elements. ....... 46 4 8 Formation of integration triangles for Heaviside f unction. ................................ ... 50 4 9 Crack opening due to special integration scheme. ................................ ............. 50 4 10 Geometry with a crack ................................ ................................ ........................ 51 4 11 Background grid generated in IBFEM. ................................ ................................ 52 4 12 Enriched elements are selected and boundary nodes are fixed. ........................ 53 5 1 IBFEM model for plate with edge crack. ................................ ............................. 55 5 2 Contour plots for displacement and Von Mises stress ................................ ........ 56 5 3 Enrichment layers for computing SIF ................................ ................................ .. 60 5 4 Crack tip opening convergence for Abaqus and singularity type I enrichment using fixing of nodes approach. ................................ ................................ .......... 61 5 5 Crack tip opening convergence using corrected XFEM approach. ..................... 62 5 6 Crack tip opening convergence using fixing of nodes. ................................ ........ 62 PAGE 9 9 5 7 Crack tip opening convergence using corrected XFEM. ................................ ..... 63 5 8 Crack tip opening convergence for comparing fixing of nodes against corrected XFEM ................................ ................................ ................................ .. 64 5 9 Crack tip opening convergence for comparing fixing of nodes against corrected XFEM ................................ ................................ ................................ .. 64 5 10 KI Convergence as a function of enrichment terms ................................ ............ 65 5 11 KII Convergence as a function of enrichment terms ................................ ........... 65 5 12 Figure showing number enrichment layers ................................ ......................... 66 5 13 Effect of radius of enrichment on crack tip opening (problem 1) ......................... 67 5 14 Plate with center crack under far field loading. ................................ ................... 68 5 15 IBFEM model for plate with center crack. ................................ ........................... 68 5 16 Contour plots for displacement and Von Mises stress ................................ ........ 69 5 17 Crack tip opening convergence for Singularity Type I enrichment using fixing of nodes ................................ ................................ ................................ .............. 72 5 18 Crack tip opening convergence for Singularity Type I enrichment usin g corrected XFEM ................................ ................................ ................................ .. 72 5 19 Crack tip opening convergence for Singularity Type II enrichment using fixing of nodes ................................ ................................ ................................ .............. 73 5 20 Crack ti p opening convergence for Singularity Type II enrichment using corrected XFEM ................................ ................................ ................................ .. 73 5 20 KI convergence as a function of enrichment terms ................................ ............. 75 5 21 Effect of radius of enrichment on crack tip opening (problem 2) ......................... 76 PAGE 10 10 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the R equirement s for the Degree of Master of Science IMPLEMENTATION OF EXTENDED FINITE ELEMENT METHOD USING IMPLICIT BOUNDARY APPROACH By Vishal Hundraj Hotwani December 2011 Chair: Ashok Kumar Cochair: Nagaraj Arakere Major: Mechanical Engineering XFEM or eXtended f inite element method is a very well known technique and is getting more popular due to its vast application domain. It is a modification of finite element method (FEM) where problems having a local phenomenon such as kinks, stress concentration, and singul arity in the solution are studied. XFEM has been most extensively applied to solve problems in solid mechanic s involving stress concentration at crack tip If FEM is used to solve such a problem then the crack would be represented by a mesh element edges h ave to align with crack XFEM is a mesh independent method and hence it allows crack to pass through the elements. IBFEM is a mesh independent method that uses a background mesh instead of a conforming mesh that represents the geometry of the object. A sch eme has been discussed to incorporate XFEM in IBFEM in order to obtain a truly mesh independent approach. Such an approach would make crack as well the geometry independent of mesh. XFEM locally enriches FEM solution by incorporating priori known analytic al solution in certain r egions. This enrichment should blend and merge with the regular finite element solution of the surrounding region. Blending of enriched solution structure PAGE 11 11 with regular finite element structure has always presented problems such as poor accuracy and affected convergence rate. The m ost popular app roach to deal with problem of blending is using a weighted function with a local support A new method suggested which in essence is a technique to bring the enriched sol ution to zero along the boundary of enriched region. In order to apply Dirichlet boundary conditions in regions where the solution is enriched, the solution must be shifted such that the solution at the nodes is equal to the nodal values of displacement A ramped Heaviside function is introduced for modeling discontinuity or crack within an element This in combination with fixing of nodes can entirely avoid the need for shifti ng the solution and the need for special blending elements. Exact analytical solu tion is used at crack tip element for obtaining SIF directly without any post processing or contour integ rals computation. This requires that all the enriched degrees of freedom of crack tip element be equal. to equate the values of enriched degrees of freedom of crack tip element which is usually achieved by using penalty method. PAGE 12 12 CHAPTER 1 INTRODUCTION Overview Finite element method has been used for several decades for solving various boundary value problem s. It is an important analysis tool for engineers and helps in solving complex boundary value problems. Several FEA programs are available commercially and are widely used in a lot of industries. Results obtained are now reliable and accepted by the engine ering world. Although it still has some limitations and certain modifications are required to increase its scope. This in turn has led to development of methods that are modification of FEA but are more specific to certain applications. Extended FEM (or X FEM) which was first discussed by Belytschko [6] is one such method that can be used for problems that have a priori known solution for local domains. Usually problems involving local phenomenon such as jumps, kinks and singularity are hard to model using regular FEA programs as normal polynomial functions used as trial and test functions would not capture the local phenomenon. XFEM helps to incorporate a known local phenomenon while solving the overall partial differential equation. Partition of unity met hod (PUM) by Babuska [3] was a step in the direction of locally altering the test and trial functions and incorporating the change in the overall finite element structure. It was a very important step that helped FEA to develop in lot of different areas. T he core idea was to use functions other than polynomials to approximate the local changes in solution in order to improve accuracy. PUM requires that in order to localize certain functions in finite element space they should be multiplied with shape functi ons whose sum at any point in unity. XFEM uses PAGE 13 13 P UM to solve problems involving k inks, jumps and singularities. One such important field of application is fracture mechanics that studies structures involving cracks. XFEM uses specific functions obtained fro m known analytical solution as enrichment function to incorporate priori known solution into FEA framework. Fracture mechanics studies propagation of cracks in material. It mostly uses experimental solid mechanics to characterize materials resistance to f racture. It is a very important field as it studies particular mode s of failure in different materials and offers a chance to predict it. Several different approaches have been presented to parameterize the fracture of materials and one such popular appro ach is to obtain the stress intensity factors associated with cracks. Different mathematical models have been developed that predict the crack growth based on the computed stress intensity factors. SIF completely characterize the state of stress near a cra ck tip in a linear elastic material. It is known that material fails when SIF reaches a critical value which is an alternate measure of fracture toughness. Thus it is very essential to be able to compute SIF values of models correctly to predict crack grow th. In this thesis, the progress made in formulation of XFEM for computing SIF was studied. Techniques were developed to implement XFEM in a mesh independent analysis approach called Implicit Boundary Finite Element Method (IBFEM) which uses a background structured mesh for the analysis. The geometry and the crack are modeled independent of the mesh. Although the implementation is limited to fracture mechanics the intention is to understand the idea behind the concept and suggest certain improvements based on results obtained. In order to model a crack it is essential to take care of two important factors namely discontinuity and singularity. Functions PAGE 14 14 suggested by Belytschko [6] for taking care of these factors are very versatile and produced good results. Although it could be seen that for kinked and curved cracks it requires mapping of co ordinates to calculate the parameters required to compute the integration functions which besides being sophisticated is computationally very expensive. Later papers hav e shown that use of signed distance function as enrichment function provides solution to this problem away from the crack tip and ensures discontinuity .Use of level set function for the same purpose has been advocated on the basis of two reasons. Firstly it helps identify the region where enrichment is required and secondly it is also used to define the enrichment function. Level set function was first proposed by Osher [47] where it was suggested that it is very useful in modeling moving interfaces. In th is thesis a method is presented that uses parametric equation of boundary of the model and avoids computation of level set function since this parametric equation can be used for identifying the region of enrichment as well computing signed distance functi on. This approach would be very beneficial in problems involving crack growth as the new surface formed would update its parametric equation in the process without any extra effort. Blending of enriched solution with regular finite elements requires a lot of attention as it affects the convergence and the accuracy of solution. Blending was addressed by suggested to solve this problem. An approach based on the same method was also d iscussed for blending of enrichment for discontinuity with singularity enrichment. This method is very effective but requires the solution to be shifted properly as well as it increases the order of enrichment function. A simple but effective method is pre sented PAGE 15 15 which takes care of problem of blending, is easy to implement and saves computational effort required. Nodes at the interface of enriched domain and finite element are identified and degrees of freedom corresponding to enrichment functions of those nodes are fixed. Thus at this interface the enriched part of solution is forced to be zero by means of pre assigning the values associated to enriched coefficients equal to zero. This not only maintains the order of enrichment function but it also reduces computational effort as the rows and columns from global stiffness matrix are reduced. The same strategy is applied for blending Heaviside function to singularity enrichment. g the value of either function from unity to zero, they are gradually brought down to zero inside same element b y fixing the nodal values To apply Dirichlet boundary conditions it is known that the solution has to be shifted in order to make sure as the enriched solution possess es kronecker property. IBFEM by Kumar [13] uses implicit boundary method (IBM) to apply essential boundary conditions which does not require the nodes to be on the boundary. In present method use of shifting and special blending elements has been completely avoided by fixing of nodal enrichment degrees of freedom. Ramped Heaviside function was used for enriching elements with a discontinuity. It is constructed in such a way that its value is zero at all the nodes of the enriched element. Thus use of these fu nctions would result in enrichment field with no blending elements and enrichment itself would die to zero along all four nodes. This thesis includes information necessary for effective implementation of XFEM in an already existing FEA program. The resu lt provided at the end act as benchmark to PAGE 16 16 prove the validation of the approach. Results are discussed in terms of displacement convergence and various suggestions made are compared and are enlisted in the conclusion. SIF are also directly obtained using t he enrichment functions suggested by X.Y. Liu [32] and studied as a function of number of enrichment terms used from the analytical solution. It is strongly recommended to use the s e functions as opposed to using the functions suggested by Belytschko [6 ] wh ich use only the first term in the analytical solution. Even if the intention is to model the crack and use contour integrals instead of directly obtaining the SIF without any post processing it is still advisable to use the later approach since it has a v ery good convergence rate compared to the other method. Goals and Objectives The most prominent goal of present research can be considered as an attempt to study and understand the developments made in modeling domains in solid mechanics with cracks. XFE M is implemented in IBFEM framework to make both crack as well as geometry mesh independent. The main objectives of this thesis are: To combine the advantages of IBFEM and XFEM approaches to model fracture mechanics problems To use parametric equations to represent cracks inste ad of level set functions as done in most XFEM implementation To development strategies for implementation of XFEM such that multiple terms from the analytical solutions can be used as enrichment and SIF can be calculated as a bypro duct of the analysis. To study alternate approaches for blending enrichment with the unenriched finite element solution. PAGE 17 17 Outline Th e remaining document can be summarized as follows: In C hapter 2 there is a detailed discussion of XFEM from its introducti on to its latest developments. Blending has been thoroughly explained and problems associated to it are listed. XFEM is then studied in context of cracked domains in fracture mechanics and hence there is a brief discussion of Linear Elastic Fracture Mechan ics (LEFM). Ch apter 3 explains the basics of m esh independent methods and their classification. Implicit boundary finite element method is discussed in detail and advantages of using it in synchronization with XFEM are documented. C hapter 4 has the det ails about the implementation scheme used for empowering IBFEM with capability of modeling fracture mechanics problems involving cracked domains. A step by step approach on implementation is discussed along with the advised changes for better performance i s discussed. Chapter 5 has the results obtained using the above implementation scheme in IBFEM software. These results are thoroughly understood using benchmark problems and computed values are compared to analytical solution as well as those obtained fr om Abaqus 6.10. Chapter 6 entails the conclusions obtained from the results and scope of future developments that need to be done for improving the efficiency and accuracy of the method. PAGE 18 18 CHAPTER 2 EXTENTED FINITE ELEM ENT ANALYSIS Introduction Finite el ement analysis is finding application in almost all branches of engineering. The increasing pool of applications has led to development of various methods that are more specific to a certain application XFEM is one such technique which focuses on capturin g local phenomenon such as weak and strong discontinuities in the solution field [6]. Early approaches in solving such problems involved use of polynomials as test functions but it required attention towards mesh refinement for obtaining reasonable results Improving mesh density has been termed as h refinement while the us e of higher order polynomials f or test functions has been named p version. J Fish [14, 15] tried to capture local phenomenon of high gradient in strai n field by using technique named as s version. A similar attempt was made by J Fish and Belytschko [16] in which a sub domain of interest with high gradient had a spectral approximation, and this corresponding mesh was overlaid on regular finite element mes h. A more recent method Multi scale Enrichment using Partition of Unity ( MEPU ) [18 ] exploiting benefit of s version and PUM was demonstrated which accounts for coupling coarse scale domain with a fine scale region without affecting sparsity of coarse field. Fish described methods such as XFE M/GFEM as sparse global enrichment method ( SGEM ) and argued that MEPU is a decent approach towards exploring use of the FEM in nanotechnology. In the last decade development in this context has led to use of analytical solutions directly as test function. Enrichment was realized using partition of unity which was first explained by Babuska et al as partition of unity method PUM [3]. This addition to regular FEM trial PAGE 19 19 he weak formulation. Whereas an alternative approach of replacing usual FEM shape functions with special functions to capture local phenomenon has been termed as domain with such func tions is computationally very expensive. This approach of using analytical solution locally has helped to capture exact solutions in models involving discontinuities. Following this path leads to division of domain into three different regions which have b most of the authors. Figure 2 1. Demarcation of enriched region and blending region In F igure 2 1 reg ion around crack colored as blue is composed of enriched eleme nts and that colored as green is composed of blending elements while rest of the elements are regular finite elements. Enriched domain is a collection of all elements that have all their nodes enriched. While blending elements are those that have some but not all enriched nodes. A clear definition for the same could be found in later section of Blending elements Enriched elements Crack PAGE 20 20 this report. Typical example where X FEM is applied can be found in f racture mechanics where displacement field shows asymptotic change while strain could be singula r at the crack tip. Several other problems such as interface problems involving fluids, contact stresses at join ts or multi material problems, s hear bands and dislocation models have also been modeled using this technique. Enriched elements have capabilit y to reproduce exact solutions depending on the type of enrichment used but blending elements do not have the same capability. Blending elements serve to be a means of ensuring compatibility betwee n enriched elements and thosewith regular shape function. P UM condition is not satisfied in blending elements hence they are unable to represent enriched function but instead they end up adding unwanted terms leading to impaired accuracy and convergence properties as discus sed by Chessa et al [21 ]. Fries presente d a method of global enrichment using a ramp function which was presented as corrected XFEM [1]. This method was further studied by Ventura who Fries approach. A J Fa wkes [4, 5] and colleagues tried to solve problems involving crack tip singularities using finite element method in early approaches. Extended finite element method presented by Belytschko [6, 7, and 8] uses analytical solutions for such solution fields ne ar crack tip which has been active area of development in past decade. Nagashima along with colleagues tried to model interface cracks between dissimilar materials using XFEM and proved its effectiveness for stress analysis and stress intensity factor calc ulation in bi material fracture problems [8]. Earlier papers [6 8, 9] used the leading or the lower order terms of asymptotic solution and calculations of PAGE 21 21 desired stress intensity factor were carried out using output results of analysis using complicated p ost processing. Karihaloo and Xiao [10, 11] made an attempt to incorporate higher order terms in enrichment function to improve accuracy of solution and also to develop technique to obtain SIF directly without any post processing. Results obtained by [10 ,11] shows it is highly accurate method compared to previous approaches and if we further limit the additional degrees of freedom due to crack tip enrichment to be equal, the solution assumes the asymptotic field and nodal enriched coefficient values of c rack tip element are the required SIF. General Formulation XFEM was first introduced with its application field concentrated on fracture mechanics. Lot of effort has been made in this direction in order to capture discontinuities in solution field which i n our case would be across the crack in a given material under loading. Discontinuities could be in form of a jump which is referred as as well as its gradient is discontinuous in case of strong discontinuities while only gradient is discontinuous in case of kinks. XFEM is also useful to reduce stiffness of approximation in the vicinity of high gradients which helps in improving accuracy without the need for mesh re finement. This reduction in sophistication of mesh generation as well as savings in computational effort is also a motivation towards XFEM. Regular FEM app roximates the solution field piece wise with polynomial shape functions. But in several cases these polynomial functions cannot approximate jumps and kinks resulting in stiff solution and degradation of accuracy. Examples of such events would be modeling of cracks for calculation of SIF, s hear bands, dislocations, PAGE 22 22 material and phase interface. In materi als possessing cracks we have a jump in displacement field across the crack while crack tip has a high gradient in stress intensity. In the case of global enrichment the entire domain is modeled with enriched elements and hence both the regular polynomial s and the enrichment functions would share same suppor t i.e. entire domain. T he shape functions used to maintain partition of unity may be of different order compare d to the regular FE approximation part. This is earlier. 2 1 Where Regular finite element shape functions. Partition of unity shape functions. Enrichment function. I Entire domain. Degrees of freedom associated with I due to enrichment. Equation 2 1 is the overall displacement obtained after enriching all the n odes in the model. Thus it can be seen that the total displacement is the summation of displacement due to regular nodal values and enriched nodal values. This approach may be simple to implement but is computationally very expensive as the entire region is enriched In most models the need of spec ial enrichment function is local for example discontinuity due to crack or stress concentration due to special feature of model. Thus in such a situation enriching entire domain makes little sense. So only a part of domain is enriched and such a formulatio achieved by identifying local nodes in need of treatment and hence only those nodes PAGE 23 23 are burdened with additional degrees of freedom. Such an approach would lead to the following formulation. (2 2) Where Enriched domain Degrees of freedom associated with I due to enr ichment. This was general structure given to XFEM approximation in [6] whereas definition of enrichment function and I depends on specific application. Also, This helps us to maintain partition of unity over entire domain which gives it capability to represent any enrichment function in domain I But this approximation loses the Kronecker property and application of essential bound ary conditions becomes a challe nge A shifting of enrichment function was proposed by Belytschko [7] which assumes following approximation of the displacement field. (2 3) is the enrichment function evaluated at node j, and if are nodal co ordinates at node j then This approximation helps us apply essential boundary conditions though the en richment might still be non zero on some part or boundary other than nodes. In case of multiple enrichment s the overall test function would take the following form. (2 4) PAGE 24 24 is the compact support for the K th enrichment. Some of the nodes are common between two or more enrichment s and hence they would pos s ess appropriate degrees o f freedom to represent the solution field. Blending O ne of the chief advantages of XFEM is that the enrichment is localized and is only effective in regions where the discontinuity is expected or known. At the borders between the enriched and normal FE d omain there is incompatibility due to additional degrees of freedom associated with the enriched domain. This leads to existence of special elements referred as blending elements that help in smooth transition between two different types of domains. The pr oblems associated with blending have been discussed in past [21, 22 and 2]. Figure 2 2 Enrichment f unction across blending element Blue nodes represent the enriched nodes while the green nodes represent the regular finite elemen t nodes. The Figure 2 2 represents weight function across the blending elements. Blending elements do not satisfy partition of unity and hence they are unable to represent enriched function but instead they act against and affect the convergence and accura cy. Most of the approaches developed in the past to handle this problem are specific to certain application and do not apply to arbitrary enrichment U( x) E1 E2 E3 E4 E5 E6 E7 E8 Blending element Regular finite element Enriched element PAGE 25 25 functions in general. A more general approach called as adopted in by Fries [22 ] that use s a smooth weight function with compact support usually in the form of a polynomial to bring enrichment down to zero at interface nodes This path leads to suppression of all parasitic terms associated with blending elements. A ramp function was multiplied with the enrichment function referred to as weight function blending by Ventura [2] This can be represented as follows. (2 5) Wherein is the modified enrichment function which can be described as being product of two functions i.e. the original enrichment function and its associated weigh t that localizes the enrichment The weighting function cho sen is such that it diminishes the enrichment away from t he region of interest. Fries [22 ] suggested using partition of unity functions for this purpose which could be described as follows, (2 6) In this mathematical expression x represents elements that ar e supposed to be enriched and I x represents nodes of element x Thus if an element has any common nodes with domain I which was defined in E quation 2 2, the n it should be considered an enriched element which results in new set of enriched nodes E The ne w enrichment function for element x would then be (2 7) PAGE 26 26 Figure 2 3. Avoiding formation of blending element s using corrected XFEM. Let W be the set of all nodes of element x so we define m = W I Thus in such a situation weight x would be one for all enriched elements as defined by regular XFEM and hence enrichment function is produced exactly. Whereas blending elements formulation cannot produce enrichment function exactly. T his new appr oach helps in gettin g rid of all the unwanted terms in the approximation of such blending elements. Thus this approach does not enable blending element to reproduce enrichment function but it helps in avoiding the parasitic effects of blending element on o verall solution and mergers the enriched solution smoothly with regular FE domain. Results obtained by Ventura [2] show that wei ght function blending improves accuracy and convergence properties. Gracie [21] tried to avoid formation of blending elements us ing Discontinuous Galerkin (DG) by coupling two regions and imposing continuity between the two domains using penalty approach. This method was compared to Assumed Strain (AS) approach discussed by Ch essa et al [21 ]. Application to Fracture Mechanics X FE M was first studied in detail to model discontinuities in a model containing crack by Belytschko [6]. XFEM in general is a technique to treat weak or strong discontinuities as well as kinks in the fiel d solution. Westergraad [25] wrote one of the U( x) E1 E2 E3 E4 E5 E6 E7 E8 Regular finite element Enriched element PAGE 27 27 most impo rtant paper s in fracture mechanics history which had the displacement and stress field of mode I crack problems which was further extended to mode II problems by Sih [26]. Westergraad argued that MKG (Muskhelishvili Kolosov Goursat) functions can be simpl ified and expressed in terms of a single function in case of special symmetry. We know that for two arbitrary functions and the Airy stress function generated using MKG functions will be bi harmonic. Here z is a complex number. (2 8) (2 9) Thus these stresses can be described in terms of as follows (2 10) (2 11) (2 12) Westergraad prescribed that for mode I shear stress should van ish along plane of symmetry xy (x, 0) = 0 and hence the following relation exists between two MKG functions (2 13) Using the above relation the general solution of stresses become, (2 14) (2 15) (2 16) PAGE 28 28 Consider an infinite plate with center crack of length 2a subjected to biaxial loading Westergraad functions given by and will satisfy the boundary condition for given state of stress and using equ ation for stresses given above and applying the appropriate boundary conditions we obtain the near crack tip stress by assuming the origin at individual crack tip for calculation of stresses at that crack tip (2 17) (2 18) (2 19) Also we define (2 20) (2 21) S imilarly for mode II we can derive correspondin g crack tip stress field. W e observe that in analytical solution the solution field has functions that are not polynomials Belytschko decided to enrich crack tip elemen t with following four functions (2 22) A nalytical so lution lies in the span of the above four functions and hence the error reduces and the analytical solution is better cap tured with fewer elements When a strong discontinuity is to be applied in solution field functions such as level set fun ction, sign distance function or Heaviside function have been u sed successfully in form of enrichment. In the present app roach Heaviside function is used as we need to assign PAGE 29 29 the ability to elements to let crack pass through them. If step enrichment is shifted it normally vanishes outside the enriched element and hence there are no blending elements. Even if approximation is not shifted, it does not affect much as the enrichment is constant in blending elements Figure 2 4 Enrichment technique for modeling cracked domain using structured grid. Nodes with red color represent singularity enrichment while nodes with yellow c olor represent Heaviside enrichment. It could be seen from the F igure 2 4 that only one element is enriched with S ingularity enrichment but as explained later a geometrical approach of enriching elements is used and a fixed region defined by a radius o f en richment is enriched for higher accuracy and improved convergence PAGE 30 30 CHAPTER 3 MESH INDEPENDENT FINTE ELEMENT ANALYSIS Developments in Mesh Independent and Meshless Methods All FEA programs represent geometry with help of a mesh. Hence the model itself is re presented as a set of elements that are aligned in such a way that they represent the shape of the object. This proves to be a challenge when a mesh is to be created for complex geometries and in such case s most automated mesh generation algorithm s generat e a poor mesh In modern day engineering sy stems are becoming more complicated and intricate which renders use of regular FEM difficult One such diffi culty is mesh distortion which means poorly shaped and grouped elements In the last two decades the foc us ha s been upon methods that avoid generating a mesh to represents the geometry of the structure. This is usually experienced when automated mesh generator algorithms are asked to generate mesh for a complex geometry. Mesh generation is particularly prob lematic for models involving singularities, h igh gradient, complex geometries, large deformation and nonlinear behavior. Mesh free methods were developed with a goal of avoiding complex mesh generation. The main objective of meshless methods was to generat e approximation based entirely on nodes without connecting them to form elements. Figure 3 1. Model for meshless method showing nodes and boundary Dispersed nodes insi de the boundary PAGE 31 31 One of the early methods in this direction was Smoothed Particle Hydrodynamics (SPH) [45, 46] It was a co mputational method for simulating fluid flow problems and it is a mesh free lagran gian method where the mesh moves with the flow. William G Hoover used Smoothed Particle Hydrodynamics for studying impact fracture in solid mechanics [27]. Element Free Gale rkin Method (EFGM) [28, 29] was developed as an alternate path that used moving least square (MLS) approximation. Nayroles et al [30] actually first used MLS function as set of approximation function for solving differential equation in Galerkin space an Babuska and M elenk [3] presented method known as Partition of Unity Method (PUM) which is a general approach to user defined domain with required approximation. Thus it could be said the SPH and EFGM are in a way instances of more general approach presented by PUM. This new structure helps formulation of variety of trail spaces to solve above mentioned problems involving higher order differential equation Mesh i ndependent methods use a back grou nd mesh or a structured grid. In such methods the mesh need not conform to the geometry. For this case various different approaches have been developed to apply boundary conditions. The Penalty boundary method (PBM) was presented by Clark and Anderson [40 41] I t explains use of a regular structured grid while essential boundary conditions are applied using a penalty method. Shapiro [42, 43] and Hollig [37, 38, and 39] proposed use of R functions to define the boundary of the problem domain and Dirichlet boundary conditions were satisfied using the solution structure proposed by Kantorovich and Krylov [31] XFEM was also developed as a mesh independent method that uses implicit equation of a PAGE 32 32 boundary to represent phenomenon such as cracks and boundary cond itions were applied using Lagrange multipliers. Figure 3 2 Mod el representing m esh independent method with a non conforming background mesh The Implicit boundary finite element method (IBFEM) was first presented by Padmanabhan [44] and Kumar [13] an d it uses Heaviside step function to construct a solution structure based on the technique proposed by Kantorovich and Krylov [31] In this C hapter there is a brief review of meshless methods and then emphasis is made on study of Implicit Boundary Method (I BM) [13] which was used in the implementation of XFEM for this thesis Implicit Boundary Finite Element Method trial function to obtain a C 0 continues solution for elastosta tic problems. This approach in FEA gives a discontinuous stress across the element boundaries which require use of various smoothing algorithm for graphics purpose. Various basis functions such as Hermite functions, Moving Least Square (MLS) or B spline fu nctions have been used to obtain higher order continuity such as C 1 or C 2 Although application of functions such Boundary value problem Background structured mesh PAGE 33 33 as B s pline functions is not simple as essential boundary conditions are difficult to apply. The solution structure so generated does not pos sess Kronecker and hence Dirichlet boundary conditions are not satisfied. Thus the use of B spline functions calls for a structur ed mesh or grid which limits its use to a certain type of geometry. Thus in order to apply boundary conditions in ca se when the grid does not confirm with the geometry one needs to come up with special scheme. Implicit boundary method [13] (IBM) is one such method that uses D functions or Dirichlet functions to apply boundary conditions using the implicit equation of t he boundary under conside ration. This approach makes sure that essential boundary conditions are applied correctly even though the geometry does not confirm with the grid. Thus a grid which does not confirm to geometry ca n be used which helps in using basi s functions such as B spline formulation or other high er order approximations that possess C 1 continuity. Grid Generation and Application of EBC In FEM essential boundary conditions are applied by specifying its value at nodes on the boundary But since IB M uses a structured grid, the boundary may or may not have a node. Thus special technique is used to apply essential boundary conditions in such a case. IBM uses implicit equations of boundary to impose Dirichlet boundary conditions. A solution structure i s constructed using Heaviside function which varies from zero to one over a narrow range of values This method was first presented by Kantorovich and Krylov [31] and is termed as Implicit Boundary Method (IBM) by Kumar et al [13]. A piece wise boundary v alue function is used in solution structure in case the essential boundary conditions are non homogenous. In case of problems dealing with PAGE 34 34 elasticity the essential boundary conditions are applied in terms of displacement and the solution structure U define d over domain R 2 or R 3 takes following form, (3 1) (3 2) If is implicit equation of boundary and u a (x) is the required boundary value function then this solution structure will enforce U(x) = u a (x) at the boundary and essential boundary conditions could be enforced These functions can be referred as weighting function or as D f unctions as mentioned in Kumar et al [13] as they help in imposing Diri chlet boundary conditions. The D functions do not directly use the im plicit equation of boundary. D fu nctions construct a step function using those implicit equations as follows, (3 3) Function D(x) is a step function which tends to be Heaviside function as The magnitude of 10 5 or smaller is used in the implementation Since D Function has value of one in all elements except boundary elements, the stiffness matrix of all internal elements are identical which increases computational efficiency The S a being the boundary for essential boundary condition and S t being boundary for traction is (3 4) PAGE 35 35 In E quation 3 4 } is virtual strain and is the virtual displacement. Thus a trial solution is constructed as (3 5) Where {u a } is a vector representing boundary value function, {u g } are grid variables while is a diagona l matrix constructed by Dirichlet functions such that variable part of solution representing grid variables vanish at boundaries. The same Dirichlet functions used to create the trial solution are used to create test functions. Thus stresses and strains be come (3 6) (3 7) Substituting test function in weak form gives (3 8) Now strain can be summarized as (3 9) (3 10) In the E quation 3 10 vector represents the nodal grid variables for element e and [B] is a strain displacement matrix. Thus the weak form can be represented in discrete form by assembling the respective equations in glob al stiffness matrix similar to finite element analysis. The system solves for the grid variables instead of displacement unlike regular finite element analysis and hence there is no limitation to use of shape functions. In such a formulation the shape func tions could have variety of basis functions such as Lagrange interpolation, B spline functions or other meshless formulation such as MLS. PAGE 36 36 CHAPTER 4 IMPLEMENTATION OF MESH INDEPENDENT XFEM The implementation scheme adopted and enrichment functions used ar e di scussed further in this C hapter An approach for avoiding blending elements by fixing of nodal variables similar to corrected XFEM [22] is proposed in this C hap ter. Corrected XFEM use s weighting function that ramps down the enriched solution gradually to zero Fixing of nodes requires the enriched d.o.f of nodes that are present at boundary of intersection of enriched and regular FE domain to be fixed to zero. Ramped Heaviside function is introduced which avoids the need for shifting and formation of bl ending elements While for singularity enrichment early approach by Belytschko [6] was used and was compared to solution obtained by X.Y. Liu and colleagues [32] for direct evaluation of mixed mode SIF. In order to obtain SIF directly using the later appro ach the enriched d.o.f of crack tip elements should be equal. X.Y. Liu used penalty method for equating the enriched d.o.f of crack tip element. In this thesis Instead of using a the connectivity of crack tip element. Ramped Heaviside Function As mentioned in Chapter 2 Heaviside function is used as step function to achieve discontinuity in the displacement solution. Signed distance function was used to determine if the Gauss poin t is inside or outside of crack boundary. Signed distance function is calculated as follows when coordinate of the closest point on crack boundary is is the sample Gauss point of interest and is the normal to the line. (4 1) PAGE 37 37 Thus depending on location of sample Gauss point with respect to crack it will either have a positive value or a negative value. The magnitude is not used but only the sign of the answer determines the value of Heaviside function as follows. (4 2) This function can be represented as follows in an iso parametric four node el ement with disconti nuity along Y =0 axis Figure 4 1. Heaviside function for strong discontinuity along Y=0 Ramped step enrichment can be summarized as follows (4 3) H is a Heaviside function whic h is described in E quation 4 2 and the weighting function can be chosen as follows. (4 4) Figure 4 2. Ramped step function with discontinuity along Y=0 PAGE 38 38 Where n is the number of nodes of the enriched element outside of crack boundary and N K is the value of k th shape function associated with that node Thus thi s function would be forced to zero at all the nodes that are inside the crack boundary. It could be seen that enr ichment dies down to zero at all the nodes and hence there no blending elements created using this formulation. Also it helps generating an element which does not require shifting. Singularity Enrichment Singularity enrichment is used to reproduce the depe ndency of analytical solution near crack tip on where is the radial distance of point with respect to crack tip. The analytical solution is given in C hapter 2 gives all the required functions to be represented in span of enrichment function. Belyts chko [6] in order to incorporate this solution field in inter polation scheme used functions given in E quation 2 22. Karihaloo and c olleagues [11 12 and 20 ] tried to use higher order terms in enrichment functions which can be written as follows (4 5) (4 6) Where are respective coefficients of nodal d.o.f associated to enrichment function The coefficients of first term represent mixed mode SIF s in Isotropic and homogenous materials i.e. represents Mode I and represents Mode II SIF PAGE 39 39 (4 7) These fu nctions when plotted assume the following shape. a) Function f 11 b) Function f 12 c) Overall assumed d) Function f 21 e) Funtion f 22 f) Overall assumed Figure 4 3. Enrichment fu nctions used for direct computation of SIF Thus ov erall approximation used after combining Heaviside function and singularity enrichment can be summarized as (4 8) is a Heaviside or sign distance function required to intr oduce a strong discontinuity in displacement field which enables crack to pass through an element. N is the number of truncated terms from the analytical solution field given in E quation 4 7. PAGE 40 40 Domain I is entire region I and I ** are the compact support of respective enrichment function The enrichment function is the exact analytical solution Therefore the nodal values of the enrichment for the crack tip element becomes the required SIF at that crack tip The nodal values are restricted to be equal so onl y one unique value of SIF is computed. T his would eliminate the extra effort required to calculate computationally costly J integral s. This approach also includes higher order terms from analytical solution which helps in achieving higher ac curacy compared to previous methods. Fixing of Nodes XFEM enables merging priori known analytical solution t o finite element solution locally for problems involving local phenomenon such as kinks or singularity The enrichment functions have compact support and hence the y are only active in certain region of model. S ome measure is to be taken in order to have a smooth transition of solution field from enriched region to regular finite element region. This is achieved by formation of blending el ements which were discussed in C hapter 2. Blending presents certain problems such as reduced accuracy and low convergence rate. In this section a technique to merge enriched solution with regular finite element solution without formation of blending elements is presented. In case of Heaviside enrichment a ll the elements that are cut by the crack are identified. E lement s con taining the crack tip are not to be enriched with step function since discontinuity in such an element is represented by singularity enrichment. Thus enrichment due to Heaviside function is brought down to zero in crack tip element. Fixing of nodes refers to constraining the value of degrees of freedom associated to enrichment function at those nodes to be zero. Equations 2 5 and 2 7 show that degree PAGE 41 41 of enrichment function is increased by degree of weight function used in corrected XFEM. Fixing of nodes doe s no t increase the degree of enrichment function as in corrected XFEM. Figure 4 4 Blending of Heaviside function using fixing of nodes. Figure 4 4 shows the location of crack boundar y and the crack tip element .F our nodes belonging to crack tip elemen t are being fixed which is represented by black squares There are in total six elements enriched with discontinuity enrichment but one of element which contains crack tip is removed from list by technique of fixing the nodes. This will ensure smooth conve rgence of discontinuity due to Heaviside function with discontinuity from Singularity enrichment functions Figure 4 5 Blending singularity function using fixing of nodes. Figure 4 5 shows the nodes that are fixed when both Heaviside and Singularity enr ichments are used. The nodes that are fixed are the nodes at the boundary of each PAGE 42 42 enrichment domain. The nodal values of enrichment are equated to zero at these fixed nodes and hence the corresponding rows and columns are removed from global stiffness matr ix. This ensures continuity as well as reduces computational effort Derivation of Stiffness Matrix In XFEM after adding additional nodes in enriched element s assembly of stiffness is done in same way as regular finite elements. Stiffness matrix could be d erived in following way. (4 9 ) In our case we have chosen = (4 10) Let be symmetric gradient operator (4 11) (4 12) (4 13) (4 14) Using this approximation in weak form of principle of virtual work we get, (4 15) (4 16) (4 17) (4 18) (4 19) PAGE 43 43 For 2D elasto static problem the displacement field assumed can be described as, (4 20) B ut for enriched element this assumed field becomes, (4 21) Strain in case of plain stress conditions can be written as follows for an enriched element. (4 22) (4 23) (4 24) PAGE 44 44 (4 25) (4 26) Whereas stiffness matrix can be derived from above as follows (4 27) Which will b e a 16 x 16 matrix in case of a four node element and it can be very easily derived using the a bove definition of B 1 and B 2 matrix. In case of singularity enrichment the number of enrichment functions is more than one. Hence the stiffness matrix of such as element should be derived accordingly. If we follow the above steps for more than one type of enrichment function on same element we would end up having following strain matrix. (4 28) PAGE 45 45 (4 29) Implementation Scheme for Multiple Enrichments It is know that some of elements will require more than one enrichment functions. Thus the model might have different number of enrichments in different elements. Therefore i t is advisable to have a versatile implementation framework that would allow us to have desirable number enrichmen ts on a particular element. If E quation 4 7 is used as the singularity enrichment, then number of enrichment functions used will depend on num ber of terms used from analytical solution. Number of enrichment functions in an element will in turn determine the number of d.o.f. In present approach this was achieved by having multiple number of elements with same d.o.f instead of having a single ele ment with multiple enrichments and large number of d.o.f. F or example an analysis using 2D four node elements will have multiple four node elements with same global co ordinates in order to create an enriched element with required d.o.f. The procedure for assembling required stiffness matrix can be understood from E quation 4 29. When element is enriched with first enrichment components without any c olor are assembled as shown in E quation 4 29. Second enrichment assembles components colored as blue. Similar ly j th enrichment will assemble components colored as green. Special elements are created that do not represent a specific enrichment function but are used to assemble missing components in the stiffness matrix that are colored as red. This assembly proced ure was achieved by generating appropriate connectivity for the elements created. IBFEM framework has PAGE 46 46 assembles components of elements stiffness at appropriate location in global stiffness matrix depending on their connectivity. In general an element in I BFEM is an object created by combining it with any of available analysis type and Interpolation type files. A new object of type AType i.e. ATEnrich was created along with a new interpolation type IEnrich which together would make enriched elements. This c ould be understood easily from following chart s Figure 4 6 Chart showing inheritance of classes and formation of elements. Figure 4 7 Chart showing inheritance of classes and formation of enriched elements. This scheme is very versatile and each en richment type can be individually altered to have different functions by changing the object file for that particular enrichment function Thus here object oriented programming helps to develop a general structure PAGE 47 47 for any type of enrichment function to be easily incorporated by overriding some basic function of super class. Some of the basic important functions that enrichment class must have could be listed as follows. Function to determine elements to be enriched with that particular enrichment function In this implementation the boundary ID for crack was used to determine elements through which an crack passes and hence those elements were enriched with Heaviside function. A radius was provided in input file which would help determine the elements enric hed with singularity enrichment. All nodes are checked if they are at a distance less than the radius entered and if they are all the nodes in that element are enriched. Function for avoiding formation of blending elements. In present implementation b len ding of solution is taken care by two means i.e. either by fixing the nodes or by using corrected XFEM. Fixing of nodes for blending singularity enrichment is achieved by using same radius that was used for enriching the elements. All the nodes that fall o utside that radius but are part of enriched element are collected in an array and enriched d.o.f associated to it are equated to zero. While corrected XFEM is implemented by using the shape functions of parent element or the underlying grid as weight funct ion. The value of weighting function is the summation of shape functions associated to nodes on which the enrichment function is supposed to be non zero. This is also achieved by using the radius of enrichment and identifying nodes that are outside the cir cle with that radius. While computing the weighting function at any point in an element the nodes of that element are checked to see if they belong to array of nodes that are out side the enriched radius. If nodes are found belonging to that array PAGE 48 48 then the value of shape function associated to that node at that point is subtracted from the sum of all shape functions at that point. Function to assemble the local stiffness and send it to global stiffness matrix when called for. This means assembly of B J matri x where j stands for the j th enrichment. H ence every enrichment class should have its own individual stiffness matrix which will be assembled in the element matrix as mentioned in above E quation 4 29 Read functions that help collect all data required by t he enrichment file to perform all the necessary steps in algorithm. This may depend on the enrichment type and hence is modified according to the enrichment function and its requirement. So Heaviside function would require the boundary ID number to identif y the crack location from the geometry file which is in VRML format exported from a CAD software. In order to generalize the enrichment all enrichment class are developed as being an object file of type enrichment which has all the functions that are commo n to all types of enrichment and only functions that are specific to a certain type of enrichment are overridden. One example of such a function would be definition of enrichment function which is different for all enrichments. Discontinuity enrichment wi ll have Heaviside function while singularity enrichment will have complex functions with singularity and trigonometric functions. Integration of Enriched Elements Most common technique to integrate functions in FEA is using Gauss quadrature. This techniqu e is very accurate in integrating polynomials if correct order of quadrature is applied. This method involves computing the function being integrated at certain fixed number of the points in domain at specified location s depending on order of integration a nd multiplying it with weights associated to that sample point. The value of these PAGE 49 49 w eights could be very easily found in literature as it is very popular method for integrating polynomials. Enriched functions usually involve analytical solution which most likely will not be polynomials. Special integration techniques have been used to integrate such functions depending on type of function being integrated. Certain steps have to be taken to assure good a ccuracy in solution as use of usual Gauss quadrature wi ll induce error s It has been noticed that if proper integration scheme is not adapted the error induced in the crack tip element is sometimes as much as total error induced in model by other elements. This is because crack tip element has singularity term along with other complicated trigonometric functions in its trial function Ventura [2] discussed a technique called into surface integration and hence is much more efficient and accurate. Although they have hinted that their method o nly works best for elements where all the nodes of element are enriched and d.o.f associated to those nodes are equal. Accuracy can be greatly increased by using polar integration [34, 35] technique w hich concentrates more sample Gauss points closer to crack tip. These methods require that element must be decomposed into integration triangles in a certain manner so as to have lot of Gauss points close to crack tip. Discontinuity enrichment needs the cr ack surfaces to be traction free. Thus the crack surface s should not induce any strain while they are being moved with respect to each other. In order to integrate step function the enriched element has to be sub divided into integration triangles in such a way that one of the edge s of the triangles lie along the crack edge or surface as shown in Figure 4 8 Since computing the value of integrand in one triangle does not affect the other triangle, these triangles can separate along the crack edge leaving th e crack surface stress free. PAGE 50 50 Figure 4 8 Formation of integration triangles for Heaviside function. Figure 4 9 Crack opening due to special integration scheme. The triangles are generated using points selected on both the crack edges Crack has t wo edges since it is modeled as actual crack with finite but very small width In the current implementation a Gauss quadrature of sixth order is used for singularity enrichment while elements enriched with Heaviside function do not need such high order of integration. Although reasonable accuracy is maintained in overall solution the PAGE 51 51 displacement s at crack tip element nodes do not converge if an accurate integration scheme is not use d Since the crack tip element displacement solution gives the required S IF directly it is essential that they are accurate enough. A quadrature of seventh order was used by X.Y. Liu [32] for obtaining the direct SIF for mixed mode fracture analysis and a very high accuracy was obtained. Implementation in IBFEM shows that the S IF computed by this manner are very sensitive to the integration technique. It is also to be noted that either increasing the order of integration or increasing t he number of integration triangles is computationally very expensive and hence it becomes a ve ry t ime consuming analysis Algorithm used for Implementation Incorporating XFEM for solving problems in solid mechanics with cracks in IBFEM framework involved following steps. Geometry is imported from any commercial CAD software in a VRML (Virtual Real ity Modeling Language). Geometry is re presented by set of parametric E quation that defines the bounding region of the model. Crack is represented by one of the boundaries in the geometry. Figure 4 10 Geometry with a crack exported from CAD software into IBFEM using a VRML type file. PAGE 52 52 A background mesh is generate d for the object The crack is represented using the parametric equation of the boundary of the model. The model generated from CAD system has a crack of finite width but the width is very small c ompare to dimensions of the model and hence it can be treated as a crack. One of the boundaries of crack is chosen as crack boundary whose parametric equation is used for further computation. Figure 4 11 Background grid generated in IBFEM. Using the pa rametric equation of the crack and the location of crack tip the elements to be enriched in the model are determined. Elements with discontinuity are determined by finding the elements though which the crack boundary passes. Elements with singularity enric hment are determined using an enrichment radius which is an input for the program. All the elements that have one or more nodes within the circle of this radius with center at the crack tip are enriched. Thus in this approach there are essentially no ble nding elements Blending of enrichment is then achieved by using constraints i.e. fixing the nodes. The nodes that are located on the interface between regular finite elements and the enriched elements are determined. The enriched degrees of freedom associ ated to these nodes are fixed or equated to zero. Crack represented by a boundary PAGE 53 53 Figure 4 12 Enriched elements are selected and boundary nodes are fixed. The nodes with black spot represent the nodes being fixed. The nodes of crack tip element are collapsed which means that conne ctivity of that element is altered such that the nodal values of enrichment are equal at all the nodes of the crack tip element. Example of crack tip element connectivity with collapsed nodes would be ( ). It could be seen tha t first four nodes correspond to regular polynomial shape function and last four nodes are enriched nodes. The last four nodes have the same node number as they are collapsed and hence all the stiffness associated to those four nodes would be added to sing le corresponding equation in global stiffness matrix. The nodal value corresponding to that collapsed node is the SIF for mixed mode fracture toughness. Local stiffness matrix of regular finite element as well as enriched elements is assembled. In case of elements with multiple enrichments local stiffness matrix is The local stiffness matrix is assembled into global stiffness matrix using the updated connectivity table. Th is is done in a similar manner as usual FEM and corresponding equation are assembled in the matrix form to be sent to solver. R egion enriched with singularity functions. Elements containing crack boundary are enriched with Heaviside function. PAGE 54 54 Special integration technique is used to calculate stiffness at Gauss points and also different approach is used in selecting int egration points. In case of elements enriched with Heaviside functions the elements are triangulated in such a manner that edges lie along the crack edges. This helps in ensuring traction free crack edges. While for singularity enrichment the element is di vided into small triangles and a 6 th order Gauss quadrature integration on each triangle. The integration techniq ue was explained earlier using Figure 4 12 and 4 13. PAGE 55 55 CHAPTER 5 RESULTS A brief discussion about implementation of XFEM in context of modeling mixed mode fracture analysis for 2D mode ls was discussed in Chapter 4 This chapter presents results that were obtained from the implementation which validates the implementation of code. The problems presented in this chapter have an alytical solution for their SIF and are checked against results obtained. Various graphs presented help us to study convergence using different possible enrichment functions Edge Crack Under Mixed Mode Loading A rectangular plate of dimension 7 x 16 is subjected to mixed mod e stresses. This problem is also used by X.Y. Liu and colleagues [32] as a verification problem. Height = 16 units. Width = 7 units. Crack length = 3.5 units. Shear stress = 1 unit. E = 100 units Figure 5 1. IBFEM model for plate with edge crack. This is a mixed mode problem with a shear force at the top edge while bottom edge of the plate is fixed. The elements used are 4 node 2D elements and mes h used is 11 x 23, 23 x 47, 47 x 95 and 71 x 161 These meshes are selected so as to generate PAGE 56 56 a pproximately square elements The deformation pattern and the stress distribution obtained from IBFEM are compared to one obtained from Abaqus a) b) c) d) Figure 5 2. Contour plots for displacement and Von Mises stress (a) Displacement plot from IBFEM (b) Von Mises stress plot from IBFEM (c) Displacement plot from Abaqus (d) Von Mises stress plot from Abaqus. Table 5 1 shows th e legend used for addressing a specific enrichment function It is to be noted that Singularity (Type I) refers to enrichment functions mentioned in E quation 2 22 and Singularity (Type II) refers to Singularity functions used by Karihaloo [18, 12] which a re also enlisted in E quation 4 7. PAGE 57 57 Table 5 1. Legend for enrichment functions. Enrichment Type Enrichment Number Heaviside 1 Ramped Heaviside 2 Singularity (Type I) 3 Singularity (Type II) 4 Table 5 2. Enrichment scheme for cases 1 to 9 (Problem 1) Ca se number 1 2 3 4 5 6 7 8 9 Enrichment Number 1,3 1,3 2,3 2,3 1,4[1] 1,4[3] 1,4[5] 1,4[7] 1,4[9] Shifting YES YES YES YES YES YES YES YES YES Weight function NO YES NO YES NO NO NO NO NO Table 5 3. Enrichment scheme for cases 10 to 17 ( Problem 1) C ase number 10 11 12 13 14 15 16 17 Enrichment Number 1,4[1] 1,4[3] 1,4[5] 1,4[7] 2,4[1] 2,4[3] 2,4[5] 2,4[7] Shifting YES YES YES YES YES YES YES YES Weight function YES YES YES YES NO NO NO NO Table 5 4 Enrichment scheme for cases 18 to 22 (Probl em 1) Case number 18 19 20 21 22 Enrichment Number 2,4[9] 2,4[1] 2,4[3] 2,4[5] 2,4[7] Shifting YES YES YES YES YES Weight function NO YES YES YES YES PAGE 58 58 Table 5 3 and Table 5 4 are used to discuss various permutations and combinations tried to study convergence of different enrichment functions and formulations used to model the same problem. Enrichment number is the legend used to specify a type of enrichment func tion used which is enlisted in T able 5 1. The number inside the square bracket [x] repre sents number of terms used from the expansion of the en richment function mentioned in E quation 4 7. This number is equal to the value of used in E quation 4 7 The following tables show the results obtained in various cases. Table 5 5 Maximum displacem ent values for case 1 to 4 (Problem 1) Case Mesh Abaqus 1 2 3 4 11 x 23 8.295 9.548 10.495 8.7 11 10.34 1 23 x 47 8.069 8.78 1 10.01 1 8.37 1 9.72 1 47 x 95 7.955 8.248 8.795 8.05 1 8.6 11 71 x161 7.924 8 .04 1 8.375 7.92 1 8.26 1 201 x 403 7.866 7.945 8.175 7.92 1 8.09 1 Table 5 6 Maximum displacement values for case 5 to 9 (Problem 1) Case Mesh 5 6 7 8 9 11 x 23 7.69 7.72 7.75 7.75 7.76 23 x 47 7.94 7.96 7.97 7.97 7. 97 47 x 95 7.93 7.94 7.95 7.95 7.95 71 x161 7.89 7.89 7.89 7.89 7.89 PAGE 59 59 Table 5 7 Maximum displacement values for case 10 to 13 (Problem 1) Case Mesh 10 11 12 13 11 x 23 7.87 7.92 7.97 8 23 x 47 8.04 8.04 8.07 8.0 8 47 x 95 7.97 7.98 7.99 8 71 x161 7.91 7.92 7.93 7.93 Table 5 8 Maximum displacement values for case 14 to 18 (Problem 1) Case Mesh 14 15 16 17 18 11 x 23 7.25 7.27 7.3 7.31 7.31 23 x 47 7.64 7.65 7.66 7 .67 7.67 47 x 95 7.76 7.76 7.77 7.77 7.77 71 x161 7.78 7.79 7.79 7.79 7.79 Table 5 9 Maximum di splacement values for case 1 9 to 22 (Problem 1) Case Mesh 19 20 21 22 11 x 23 7.4 7.45 7.51 7.54 23 x 47 7.71 7.76 7.77 7.78 47 x 95 7.79 7.8 7.83 7.83 71 x161 7.81 7.81 7.83 7.83 PAGE 60 60 Table 5 10. SIF values for case 15 to 18 (Problem 1) Case number 15 16 17 18 Mesh KI KII KI KII KI KII KI KII 11 x 23 32.48 5.68 33.04 6 .35 33.3 5.82 33.13 6 23 x 47 32.06 5.76 33.87 6.23 34.95 5.53 34.71 5.13 47 x 95 30.33 5.48 32.92 5.74 34.26 5.32 34.27 4.45 71 x161 30.03 5.33 33.07 5.48 34.53 5.14 34.8 4.38 It is to be kept in mind that the max displ acement whenever mentioned is obtained using only one layer of enriched elements around crack tip element while SIF are computed using two enriched layers of element a s shown in F igure 5 3. Figure 5 3. Enrichment layers for computing SIF Fixing no des to blend sing ularity function Fixing no des to blend Heaviside function PAGE 61 61 The squar e black spots around the crack tip element represent fixing of nodes for merging Heaviside enrichment while the second set of square black spots represent fixing of nodes corresponding to Singularity enrichment. I t could be seen that two layers of element s around crack tip element a re enriched in the F igure 5 3 The reason behind computing SIF after enriching two layers is that the accuracy is very low if only one layer is enriched Convergence S tudy Convergence of crack tip opening was studied for various cases and care was taken to just enrich the solution with one layer of enriched elements around the crack tip element in all the cases. The graphs below are obtained by comparing the change in the crack tip opening due changing of mesh. It is to be noted that in all graphs the difference in the converged value of displacements may look huge due to scale of the graph but it is 1% or less than 1% of displacement value calculated by Abaqus.. The only exceptions with 4 % difference are Case 2 and 4 and they us e Singularity Type I enrichment along with Ramped Heaviside function. Figure 5 4. Crack tip opening convergence for Abaqus and singularity type I enrichment using fixing of nodes approach. PAGE 62 62 Figure 5 5. Crack tip opening convergence for Abaqus and singu larity type I enrichment using corrected XFEM approach. Figure 5 4 and Figure 5 5 compare crack tip opening for Cases 1,2,3,4. Th ese cases use s ingularity Type I enrichment functions which is same as u sed by Abaqus But when we compare Cases 1and 3 that u ses approa ch of fixing the nodes to Case2 and 4 that use the corrected XFEM approach, it could be seen the fixing the nodes helps displacement to converge much faster. Figure 5 6. Crack tip opening convergence with singularity type II comparing Heavisid e function with ramped Heaviside function using fixing of nodes. PAGE 63 63 Figure 5 7. Crack tip opening convergence with singularity type II enrichment comparing Heaviside function with ramped Heaviside function using corrected XFEM. The data plotted Figure 5 6 a nd Figure 5 7 is obtained from Cases 8,13,17,22. They both use the Singularity enrichment (Type II) as enrichment scheme. This compar ison emphasizes fact the using s hifted Heaviside enrichment as discontinuity enrichment or using ramped Heaviside enrichmen t without shifting is same in context of convergence of solution. Thus shifting can be avoided for Heaviside enrichment completely by using the ramped Heaviside function as enri chment scheme. I t could be seen that the convergence property obtained by enri ching the elements with Singularity (Type II) enrichment i s excellen t There is a limitation to number of terms that can be used from analytical solution as enrichment functions as it has associated increase in d.o.f. I t was found from studies of other cas es that even if the Singularity (Type II) enrichment with first three terms of analytical solution is used, it has much better convergence compare to enriching elements with Singularity (Type I) enrichment. Hence it is beneficial to enrich with Type II enr ichment even though if it is not intended to obtain SIF directly as it gives better convergence than Type I enrichment scheme PAGE 64 64 Figure 5 8. Crack tip opening convergence for comparing fixing of nodes against corrected XFEM Figure 5 9. Crack tip opening convergence for comparing fixing of nodes against corrected XFEM Figure 5 8 and Figure 5 9 the convergence rates when the blending of domain is seen the converge nce obtained by fixing the nodes is equally good as that obtained by corrected XFEM. Fixing of nodes Corrected XFEM Corrected XFEM Fixing of nodes PAGE 65 65 SIF C omputations The stress intensity factors computed are enlisted in the Table 5 11 Table 5 11 SIF values as a function of enrichment terms (Problem 1) Number of enri chment terms KI KII KI /34 KII/4.55 1 24.72 5.23 0.73 1.15 3 30.03 5.33 0.88 1.17 5 33.07 5.48 0.97 1.20 7 34.53 5.14 1.02 1.13 9 34.8 0 4.38 1.02 0.96 Figure 5 10. KI Convergence as a function of enrichment terms Figure 5 11. KII Convergence as a function of enrichment terms PAGE 66 66 Figure 5 10 and Figure 5 11 show the convergence of values of KI and KII as a function of terms used from the analytical solution given in E quation 2 29. X.Y. Liu [32] has shown a better convergence for the SIF values for t he same problem. The probable cause that could result into these errors is that the SIF values are very sensitive to integration scheme used for the crack tip element. Thus in order to obtain a good efficiency for SIF a very sophisticated scheme for integr ating the crack tip element has to be adopted. Effect of Number of Enrichment L ayers There is a difference in solution when numbers of enrichment layers around crack tip element are changed Enriching a fixed geometric region of model has been referred be understood by means of figure 5 12 which explains the selection process of enriched nodes. Thus for a fixed radius of 0.12 for mesh density of 23x49 the set of enriched nodes could be seen in F igure 5 12 Figure 5 12. Figure showing number enrichment layers It could be seen that all the elements that are cut by the circle formed by the input radius are enriched. While the nodes that fall ou tside the circle are fixed so as t o blend PAGE 67 67 the solution In the Table 5 12 the results obtained by using a constant radius of enrichment are shown while th e mesh is changed in each step. Table 5 12 Effect of radius of enrichment (Problem 1) Crack Tip opening Mesh r=0.6 r=0.95 r=1.2 11 x 23 7.31 7.31 8.11 23 x 47 8.18 8.67 8.96 47 x 95 8.61 9.13 9.51 71 x161 8.74 9.31 9.7 0 101 x 203 8.8 0 9.39 9.77 121 x 243 8.82 9.41 9.82 It is clear from the T able 5 12 that solution is differen t for different radius of enrichment and it keep s on increas ing with increase in the radius E. Bechet [35] studied the converg ence with respect to error in energy norm while in this thesis it is being studied with respect to displacement. The effect of e nrichment radius on the crack tip opening for the problem under study can be better understood with Figure 5 13 Figure 5 13. Effect of radius of enrichment on crack tip opening (problem 1) PAGE 68 68 Center Crack U nder Uniform Far Field Mode I L oading Consider a s quare plate 10 x 10 with a center crack and far field tensile stress. The model can be shown in the F igure 5 14 Figure 5 14. Plate with center crack under far field loading. a = 1 unit E = 100 units. It could be seen that this problem is a pure Mode I problem. The abo ve problem was modeled in IBFEM as shown in F igure 5 15 Figure 5 15. IBFEM model for plate with center crack. In order to apply proper boundary condi tions the bo ttom edge was constrained to move in Y direction but it was free to move in X direction. While tensile pressure was PAGE 69 69 applied on top edge of required magnitude. The displacement and stress distribution obtained from the IBFEM program are compared with Abaqus a) b) c) d) Figure 5 16. Contour plots for displacement and Von Mises stress The same procedure as problem 1 was followed to study effect of changes in various formulation of trial solution. List of various enrichment functions used can be PAGE 70 70 obtained from Table 5 1 Various different cases studi ed for this example can be obtained from T able 5 13 Table 5 13 Enrich ment scheme for different cases (problem 2) Case number 1 2 3 4 5 6 7 8 Enrichment Number 1,3 1,3 2,3 2,3 1,4[7] 1,4[7] 2,4[7] 2,4[7] Shifting YES YES YES YES YES YES NO YES We ight function NO YES NO YES NO YES NO YES Only certain cases are studied as compared to previous example as only these few cases are used for obtaining the convergence graphs. Since all the cases were studied in last example and repeatability was obtaine d in their nature only those with significant importance are studied for plate with center crack under pure mode I loading. The results obtained for this example are shown in following tables Table 5 14 Maximum displacement values (problem 2) Case number Mesh Abaqus 1 2 3 4 33 x 67 0.1014 0.1075 0.114 0.107 0.1138 53 x 107 0.1003 0.1037 0.1065 0.1034 0.1063 73 x 147 0.0999 0.1022 0.1044 0.102 0.1043 93 x 187 0.0996 0.1014 0.103 0.1012 0.1029 103 x 207 0.0995 0.1011 0.1026 0.1009 0.1025 PAGE 71 71 Table 5 15 Maximum displacement values (problem 2) Case number Mesh 5 6 7 8 33 x 67 0.1011 0.1016 0.1009 0.1014 53 x 107 0.1 000 0.1006 0.0999 0.1005 73 x 147 0.0998 0.0997 0.0996 0.0996 93 x 187 0.0995 0.0997 0.0993 0.0996 103 x 207 0.0995 0.0996 0.0992 0.0995 Table 5 16 SIF (KI) values for case 5 and 7(problem 2) Case Mesh 5 7 33 x 67 1.9 0 1.9 0 53 x 107 1.87 1.89 73 x 147 1.82 1.87 93 x 187 1.79 1.79 103 x 207 1.79 1.77 Table 5 14 and Table 5 15 give the crack tip opening under different formulation and changing mesh. It could be seen that Case 6 and 8 only has displacement values even though they use singularity (Type II) enrichment scheme. The rea son being if you multiply a weight function with the analytical solution the trial solution does not represent the exact analytical solution anymore and the nodal variables do not correspond to actual SIF. Hence in order to obtain direct SIF it is essentia l that we do not use d Thus fixing of nodes was used in our implementation and it served PAGE 72 72 the purpose. Compare to previous example this example has a shorter crack length and it is also pure Mode I type problem hence solution has better conv ergence but it still follows the same trend as the last example. Convergence S tudy Figure 5 17. Crack tip opening convergence for Singularity Type I enrichment using fixing of nodes Figure 5 18. Crack tip opening convergence for Singularity Ty pe I enrichment using corrected XFEM Abaqus Heaviside Ramped Heaviside Abaqus Ramped Heaviside Heaviside PAGE 73 73 The above graphs show convergence for cases 1 to 4 all of which uses singularity (type I) enrichment scheme. It could be seen that all the cases have almost equal slopes and are very close to slope of Abaqus solution. Figure 5 19. Crack tip opening convergence for Singularity Type II enrichment using fixing of nodes Figure 5 20. Crack tip opening convergence for Singularity Type II enrichment using corrected XFEM The above graphs are plot of crack tip opening with changing mesh for cases 5 to 8 and all of them use singularity (type II) enrichment. It is very clear that compare to PAGE 74 74 cases 1 to 4 these graphs show a much better convergence. It is apparent that cases 5 and 6 wherein fixing of nodes approach is used conve rge better compared to cases 7, 8 and also Abaqus. Thus it can be concluded that enr ichment functions suggested in E quation 2 29 are better even for problems where SIF is to be calculated using contour integrals. Above two graphs also present the differenc e in fixing the nodes and using a weighting function for enforcing continuity. It could be seen that if we use approach of fixing the nodes once again the convergence is faster compared to other approach. Also to validate the use of ramped Heaviside enrich ment it could be seen that it has same convergence except the solution is offset by a small value. Thus shifting of trial solution can be avoided entirely by using a ramped Heaviside function and fixing the nodes along the domain border of singularity enri ched region. SIF C omputation The SIF can be computed when using the singularity (type II) enrichment directly without using any contour integrals. Analytical solution for the above problem could be found in Anderson TL. [36]. (5 1) (5 2) is the crack angle and is zero in this case. Hence the mode II stress intensity factor will be zero for this problem. While can be computed as follows for present problem. (6 3) The results were obtained from IBFEM using singula rity (type II) enrichment and a topological enrichment approach. Two layers of elements around the crack tip element PAGE 75 75 were enriched with singu larity enrichment an d compatibility was achieved using fixing of nodes. The number of terms used from analytical solution given in E quation 2 29 were changed keeping the mesh constant Table 5 17 SIF as a function of number of enrichment terms (problem 2) Number of enrichmen t terms KI KI /1.77 1 1.11 0.63 3 1.67 0.94 5 1.75 0.99 7 1.77 1.00 9 1.82 1.03 11 1.82 1.03 These values obtained were using a constant mesh of 103 x 207. A s in previous example it could be seen that when only one term is used in enrichment scheme the answer obtained is very erratic. It gets accurate from fifth term and stays accurate for higher number of terms used. Figure 5 20. KI convergence as a function of enrichment terms Thus it could be seen that by using Gauss quadrature with sixth ord er for integration we can get accuracy of 3% with use of five or more enrichment terms. PAGE 76 76 Effect of Number of Enrichment L ayers This section discuss es about the difference in topological enrichment scheme against geometrical enrichment. Similar approach as previous example was use d and the results obtained by varying radius of enrichment against mesh is studied. Results obtained by following a geometrical approa ch are listed in T able 5 18 Table 5 18. Effect of number of enrichment layers (problem 2) Max Di splacement Mesh r=0.1 r=0.15 r=0.2 33 x 67 0.099 0 0.1009 0.1009 53 x 107 0.0991 0.101 0 0.102 0 73 x 147 0.0997 0.1013 0.1023 93 x 187 0.1003 0.1014 0.1023 103 x 207 0.1003 0.1014 0.1023 133 x 267 0.1003 0.1014 0.1023 Figure 5 21. Effect of radius of enrichment on crack tip opening (problem 2) Three fixed radius of enrichment were used and mesh was changed to see the results. All three different radius used gave same convergence rates but the c rack tip PAGE 77 77 opening is offset by a small value with increase in the radius. Figure 5 21 is in confirmation with conclusions drawn from previous example where a similar study was conducted. In case of radius the first val ue is to be ignored since radiu s is too small for element size that only one element is enriched making it same as the first value for It is evident that the crack tip opening converges to a higher value with increase in the radius of enrichment. Heaviside function can be very e asily integrated with a very high accuracy with techniques mentioned in Chapter 4 and hence the discontinuity is very accu rately represented. But if there is no special integration technique available for singularity enrichment and usual higher order qua dr ature is being used then choosing value of radius of enrichment becomes a matter of discussion. PAGE 78 78 CHAPTER 6 CONCLUSION S ummary This thesis is description of an effort made towards combining XFEM with implicit boundary approach. The obvious motivation behind is that model as well as the crack would be entirely mesh independent. In most of FEA packages the geometry is represented by means of finite element mesh but in IBFEM geometry is represented by a set of equations in geometry file generated from standard C AD software. IBFEM uses a back ground structured mesh which is independent of geometry. In present implementation the crack was represented as a feature in the model and is represented by a boundary. Thus the crack has a parametric equation which is a part of CAD file imported into IBFEM. This helps in recognizing enriched elements as well as determining the value of Heaviside function at any point. This procedure followed eliminates computation of Level set function. Some of the major issues with XFEM suc h as Blending and need for Shifting were discussed. A method for handling blending of enriched solution to regular finite element Fixing of nodes was discussed in detail which proves to be very simple to implement technique wherein boundary conditions are applied to enriched degrees of freedom so as to bring the enriched solution down to zero alo ng the boundary of intersection between enriched and regular domain. In present implementation scheme it could be seen that essentially there are no blending elements which is same as corrected XFEM. Fixing of Nodes unlike corrected XFEM also does not incr ease the order of polynomial PAGE 79 79 Shifting of test and trial functions is done in order to f acilitate application of Dirichlet boundary conditions. Ramped Heaviside function along with implicit boundary method was used in order to avoid shifting completely. Shifting was completely avoided and results obtained show that solution has same convergen ce rate as compared to those obtained by shifting the solution. Thus it is beneficial to avoid shifting since it reduces the complications in the trial solution used and makes it computationally cheaper. Also a detailed documentation of all the difficulti es tackled for implementation of XFEM in an already existing FEA program are presented. IBFEM is a structured grid method and uses implicit boundary approach for applying essential boundary conditions. Implementation details of incorporating XFEM with IBFE M are discussed and benefits associated to it presented. SIF are calculated directly with 3% accuracy using direct analytical solution for singularity enrichment. This method has been studied in detail and is found very effective as it completely avoids c omputation of contour integrals for SIF computations. This method requires that enriched degrees of freedom associated to crack tip elements to be equal. This was previously done by using a penalty function. An alternate approach is used which is termed as collapsing the nodes and is implemented by changing the connectivity of crack tip element and collapsing all the enriched nodes associated to crack tip element to one single node. Also a detailed convergence study is done to compare displacement convergen ce for different enrichment functions and it is validated from results that using exact PAGE 80 80 analytical solution instead of just first term had better convergence properties. Even if direct computation of SIF is not considered it is beneficial to use exact anal ytical solution with three or more terms. Also an effort is made towards studying the impact of geometrical enrichment scheme for singularity enrichment. It was observed as a larger geometrical portion of domain was enriched the material at crack tip becam e softer and thus the crack tip opening increases with increase in radius of enrichment. This leads us to the discussion as what is the optimum region to be enriched with singularity enrichment. Integration scheme to be used for integrating special functio ns for crack tip enrichment was discussed. A conclusion on impact of integration technique used to the accuracy of this method was discussed. It is clear that in order to obtain high accuracy for computation of SIF directly without computing contour integr als it is essential to have a very accurate integration scheme. Using Gauss quadrature as high as up to sixth order gives reasonable accuracy and is to be compensated with use of large number of Gauss points which indeed is computationally challenging. S cope of Future Work The above implementation is applied for static problems but the advantages presented by technique will be really effective in problems involving crack growth If a structured grid method is used in combination of XFEM then no remeshing would be required at any stage. This is because the deformed object does not affect the back ground mesh which is not a part of the geometry. Although the enriched elements will have to be changed as the crack tip is relocated. While in all other approache s the entire mesh has to be regenerated and level set functions are redefined. Also concept of localized mesh refinement could be used in conjunction with XFEM to improve the accuracy in regions involving enrichment. Localized mesh refinement would result into PAGE 81 81 higher local mesh density in regions involving enrichment or stress concentrations. This automatic mesh refinement should be very easy to apply since IBFEM is a structured grid method and the resulting quality of mesh will be very good since all eleme nts are similar in shape and form. B spline elements could be used for XFEM which might result in increased accuracy. B spline elements are known to have a very high convergence rate and hence less number of elements should be needed to achieve same accura cy. PAGE 82 82 LIST OF REFERENCES [1] T.P Fries. A corrected XFEM approximation without problems in blending elements. Internat. J. Numer. Methods.Engrg., 75: 503 532, 2008. [2] G.Ventura, R.Gracie and T.Belytschko. Fast integration and weight function blending in extended finite element method. Internat. J. Numer. Methods Engrg., 77: 1 29, 2009. [3] I.Babuska, and J.M.Melenk. A partition of unity method. Internat. J. Numer. Methods Engrg., 40: 727 758, 1997. [4] Fawkes AJ, Owen DRJ, Luxmoore AR, An assessment of crack tp singulari ty models for use with isoparametric elements. Engineering Fracture Mechanics, 11:143 149, 1979 [5] Owen DRJ, Fawkes AJ, Engineering Fracture Mechanics : Numerical methods and application. Pineridge : Swansea 1983. [6] T. Belytschko, T Black. Elastic crack growth in finite elements with minimal remeshing. Internat. J. Numer. Methods.Engrg., 45: 601 620, 1999. [7] Belytschko T, Moes N, S Usui and C. Parimi. Arbitrary discontinuities in finite element Internat. J. Numer. Methods.Engrg., 46: 131 150, 1999. [8] Moes N, Dolb ow J, Belytschko T. A finite element method for crack growth without remeshing. Internat. J. Numer. Methods.Engrg., 46: 131 150, 1999. [9] Daux C, Moes N, Dolbow J, Sukumar N, Belytschko T. Arbitrary ranched and intersecting cracks with extended finite element method. Internat. J. Numer. Methods.Engrg., 48: 1741 1760, 2000. [10] Nagashima T, Omoto Y, Tani S. Stress intensity factor analysis of interface cracks using XFEM. Internat. J. Numer. Methods.Engrg., 56: 1151 1173, 2003. [11] Karihaloo BL, Xiao QZ. Direct evaluati on of accurate SIF using PUM. In Advances in Fracture Research ( Proceedings of ICF10), Ravi Chandar K et al. (eds). Pergamon: New York, 2001: 0177OR. [12] Xiao QZ, Karihaloo BL. Direct evaluation of accurate coefficients of the linear elastic crack tip asymptot ic field. Fatigue and Fracture of Engineering Materials and Structures 2003: 25:719 729. [13] Kumar A, Padmanabhan S, Burla R. Implicit boundary method for finite element analysis using non conforming mesh or grid. [14] Fish J. The s version of the finite element me thod. Computers and Structures 1992; 43:539 547. PAGE 83 83 [15] Fish J, Nath A. Adaptive and hierarchical modeling of fatigue crack propagation. International Journal for Numerical Methods in Engineering 1993; [16] Belytschko T, Fish J, Bayliss A. The spectral overlay on the finite element solutions with high gradients. Computer Methods in Applied Mechanics and Engineering 1990; 81:71 89. [17] Lee S H, Song J H, Yoon Y C, Zi G, Belytschko T. Combined extended and superimposed finite element method for cracks. International Journal for Numerical Methods in Engineering 2004; 59:1119 1136. [18] Fish J, Yuan Z. Multiscale enrichment based on partition of unity. International Journal for Numerical Methods in Engineering 2005; 62:1341 1359. [19] Fan R, Fish J. The rs method for material failure sim ulations. International Journal for Numerical Methods in Engineering 2008; 73:1607 1623. [20] B.L. Karihaloo and Q.Z. Xiao, Accurate determination of the coefficients of elastic crack tip asymptotic field by a hybrid crack element with p adaptivity, Engng. Frac t Mech 2001; 68:1609 1630. [21] J. Chessa, H. Wang, and T. Belytschko. On the construction of blending elements for local partition of unity enriched finite elements. Internat. J. Numer. Methods Engrg., 57:1015 1038, 2003.20. [22] T.P. Fries. A corrected XFEM app roximation without problems in blending elements. Internat. J. Numer. Methods Engrg., 75:503 532, 2008. [23] R. Gracie, H.Wang, and T. Belytschko. Blending in the extended finite element method by discontinuous Galerkin and assumed strain methods. Internat. J Numer. Methods Engrg., 74:1645 1669, 2008. [24] N. Moes, M. Cloirec, P. Cartraud, and J.F. Remacle. A computational approach to handle complex microstructure geometries. Comp. Methods Appl. Mech. Engrg., 192:3163 3177, 2003. [25] Westergaard H. M. Bearing pressu res and cracks, J. appl. Mech. 6:49 53, 1939. [26] Sih G. C. On the Westergaard method of crack analysis, J. Fract. Mech. 2:628 631, 1966. [27] Hoover, W. G. (2006). Smooth Particle Applied Mechanics: The State of the Art, World Scientific. [28] Belytschko T., Lu Y.Y., Gu L., 1994. Element free Galerkin methods. International Journal for Numerical Methods in Engineering 37(2) 229 256. PAGE 84 84 [29] Belytschko T., Krongauz Y., Organ D., Fleming M.,Krysl P., 1996. Meshless methods: an overview and recent developments. Computer Methods i n Applied Mechanics and Engineering 139(1 4) 3 47. [30] Nayroles, B., Touzot, G. and Villon, P. "Generalizing the finite element method: diffuse approximation and diffuse elements", Computational Mechanics 10, p. 307 318, 1992. [31] Kantorovich L.V., Krylov V.I., 19 58. Approximate Methods of Higher Analysis. Interscience, New York. [32] X.Y. Liu, Q.Z. Xiao, and B.L. Karihaloo. XFEMfor direct evaluation of mixed mode SIFs in homogeneous and bi materials. Internat. J. Numer. Methods Engrg., 59:1103 1118, 2004. [33] Wilson WK. Combined mode fracture mechanics. Ph.D. Thesis, University of Pittsburgh, 1969. [34] P. Laborde, J. Pommier, Y. Renard, and M. Sala¨un. High order extended finite element method for cracked domains. Internat. J. Numer. Methods Engrg., 64:354 381, 2005. [35] E. Be chet, H. Minnebo, N. Mo¨es, and B. Burgardt. Improved implementation and robustness study of the X FEM for stress analysis around cracks. Internat. J. Numer. Methods Engrg., 64:1033 1056, 2005. [36] Anderson TL. Fracture Mechanics Fundamentals and Applications (2nd edn). CRC: Boca Raton, FL, 1995. [37] Hllig K., Reif V., Wipper J., 2001. Weighted extended B spline approximation of Dirichlet problems. SIAM Journal on Numerical Analysis 39(2) 442 462. [38] Hllig K., Reif V., 2003. Non uniform web splines. Computer Aided G eometric Design 20(5) 277 294. [39] Hllig K., 2003. Finite element methods with B splines. SIAM, Philadelphia [40] Clark B.W., Anderson D.C., 2003a. The penalty boundary method. Finite Elements in Analysis and Design 39(5 6) 387 401. [41] Clark B.W., Anderson D.C., 200 3b. The penalty boundary method for combining meshes and solid models in finite element analysis. Engineering Computations 20(4) 344 36 [42] Shapiro V., 1998. Theory of R functions and applications: A primer. CPA Technical Report CPA88 3, Cornell Programmable Automation. Sibley School of Mechanical Engineering. PAGE 85 85 [43] Shapiro V., Tsukanov I., 1999. Meshfree simulation of deforming domains. Computer Aided Design 31(7) 459 471. [44] Padmanabhan, S., 2006. Analysis using non conforming structured grid and implicit boundary me thod. Doctoral Dissertation, University of Florida. [45] Lucy L.B., 1977. A numerical approach to the testing of the fission hypothesis. The Astronomical Journal 82(12) 1013 24. [46] Gingold R.A., Monaghan J.J., 1982. Kernel estimates as a basis for general particle methods in hydrodynamics. Journal of Computational Physics 46 429 53. [47] Osher, Stanley J.; Fedkiw, Ronald P. (2002). Level Set Methods and Dynamic Implicit Surfaces Springer Verlag. ISBN 0 387 95482 1. PAGE 86 86 BIOGRAPHICAL SKETCH Vishal Hotwani was born in Ahm edabad, India which is in state of Gujarat. He did Engineering from Gujarat University, L.D College of Engineering. He finished his Masters in Science in Mechanical engine ering from University of Florida, Gainesville, USA in 2011. His areas of interests are CAD/CAM and finite element analysis. 