Continuous discontinuous approach for the modelling of ductile fracture

UNIVERSIDADE DO PORTO Continuous-Discontinuous Approach for the Modelling of Ductile Fracture Mariana Rita Ramos Seabra Faculdade de Engenharia Departamento de Engenharia Mecˆanica October 2012 “Crave for a thing, you will get it. Renounce the craving, the object will follow you by itself.” Swami Sivananda (1887-1963) Abstract Ductile metals, such as aluminium or steel, are present in all sorts of industries, from aeronautical to automotive, construction to household goods. A deeper understanding of the fracture mechanisms associated to these type of materials is essential to optimize production processes as well as avoid catastrophic failures. Traditional fracture theories generally rely on a single energetic parameter to trigger crack propagation, which may not be adequate for ductile metals, as there is a substantial amount of plastic deformation prior to failure. On the other hand, Continuum Mechanics theories successfully handle large straining and describe successfully the plastic-hardening and plastic-softening stages of material behaviour. Nevertheless in the final stages of failure, a discontinuous methodology is essential to represent surface decohesion and macro-crack propagation. The aim of this work is to build a model for ductile fracture gathering the advantages of both continuous and discontinuous approaches. The Extended Finite Element Method (XFEM) is combined with the Lemaitre model for ductile damage in a way that crack initiation and propagation are governed by the evolution of damage. The model was built under a finite strain assumption and a non-local integral formulation is applied to avoid pathological dependency of results on the spatial discretisation. Special attention is given to the necessary adaptations of the XFEM to be used in a ductile fracture context, such as selection of the enrichment functions and numerical integration strategies. The energy consistency issue during the transition from damage to fracture is also addressed. The effects of the introduction of a cohesive law in the global ductile fracture model are analysed and a strategy to approximate the fully damaged material condition is proposed. Throughout the work, the efficiency of the proposed model is evaluated through various numerical examples and general conclusions are outlined. To finalize this thesis, some comments on the overall performance of the model are given together with some suggestions for future work. Resumo Os metais d´ucteis, como o alum´ınio ou o a¸co, est˜ao presentes nas mais diversas ind´ustrias, da aeronautica `a autom´ovel, da constru¸c˜ao aos bens de consumo. O aprofundamento dos conhecimentos relativos aos mecanismos de fractura associados a este tipo de metais ´e essencial para optimizar processos de produ¸c˜ao e evitar falhas catastr´oficas. As teorias cl´assicas da fractura baseiam-se num ´unico parˆametro energ´etico que pode n˜ao ser adequado a metais d´ucteis, pois estes apresentam deforma¸c˜oes pl´asticas significavas anteriores `a propaga¸c˜ao de fendas no material. Por outro lado, as teorias de Mecˆanica dos Meios Cont´ınuos descrevem adequadamente as fases de endurecimento e amaciamento pl´asticos. Contudo, nas fases finais de degrada¸c˜ao do material, uma formula¸c˜ao descont´ınua ´e fundamental para descrever a separa¸c˜ao entre as superf´ıcies de uma fenda em propaga¸c˜ao. Este trabalho tem como principal objectivo o desenvolvimento de um modelo de fractura d´uctil, reunindo as vantagens dos modelos cont´ınuos e descont´ınuos. O M´etodo dos Elementos Finitos Estendidos (XFEM) ´e associado ao modelo de dano d´uctil de Lemaitre, de forma a descrever a inicia¸c˜ao e propaga¸c˜ao. O modelo considera grandes deforma¸c˜oes e uma formula¸c˜ao n˜ao-local integral, de forma a evitar a dependˆencia patol´ogica dos resultados relativamente `a discretiza¸c˜ao espacial. Para conseguir uma transi¸c˜ao energeticamente consistente entre o dano e a fractura, analisam-se os efeitos da introdu¸c˜ao de uma lei coesiva no modelo de fractura d´uctil global, o que permite o desenvolvimento de uma estrat´egia para aproximar a condi¸c˜ao de total degrada¸c˜ao do material. Ao longo do trabalho, a eficiˆencia do modelo proposto ´e avaliada atrav´es de v´arios exemplos num´ericos. Encerra-se esta tese com importantes coment´arios gerais e sugest˜oes de desenvolvimentos futuros. Agradecimentos Ao meu orientador, Prof. Jos´e C´esar de S´a, Professor Catedr´atico da Universidade do Porto, por me ter acompanhado ao longo desta fase da minha vida, prestando todo apoio a n´ıvel cient´ıfico, com muita amizade e generosidade. Tamb´em por me proporcionar a oportunidade de participar em importantes congressos na ´area da Mecˆanica Computacional. ` A Funda¸c˜ao para a Ciˆecia e Tecnologia, pelo importante apoio financeiro para o desenvolvimento deste trabalho atrav´es da bolsa de doutoramento SFRH/BD/43798/2008. To the Slovene team: my colleague and friend Primoˇz ˇ Suˇstariˇc for the precious help programming in the AceGen environment and for the valuable comments and suggestions; Prof. Tomaˇz Rodiˇc which provided all the conditions for this cooperation; to Prof. Joˇze Korelc for kindly allowing me to use his AceGen system and for all the support provided. Najlepˇsa hvala! Aos meus amigos e tambm colegas de doutoramento Filipe Andrade, F´abio Reis, Ana Neves e Jaime Rodrigues, n˜ao s´o pelas valiosas discuss˜oes e sugest˜oes, mas tamb´em pelos momentos de descontra¸c˜ao, sem os quais teria sido muito mais dif´ıcil a realiza¸c˜ao deste trabalho. A toda a comunidade da FEUP por me proporcionar excelentes condi¸c˜oes de trabalho, num ambiente tranquilo e familiar. ` A minha fam´ılia, em especial aos meus pais e ao meu irm˜ao, por me amarem e apoiarem em tudo. Ao meu Grande Mestre do Yoga, Eng. Jorge Veiga e Castro e `a minha professora, Eng. Catarina Ferreira, por toda a sabedoria que me transmitiram, em especial as t´ecnicas que me permitiram aumentar a concentra¸c˜ao e melhorar o meu desempenho profissional em geral. ix List of Figures xvi 4.16 Regular Gaussian quadrature rules applied to an element crossed by a discontinuity: a) four Gauss points; b) nine Gauss points; c) sixteen Gauss points; d) nine Gauss points. . . . . . . . . . . . . . . . . . . . . . . 58 4.17 Application of the Schwarz-Christoffel conformal mapping to build an integration rule for elements crossed by a discontinuity: a) cracked element; b) polygon that will be mapped to the unit disk; c) unit disk containing the Gauss point positions; d) final Gauss point distribution on the upper partofelementa). ............................... 59 4.18 Integration rules built using the Swchwarz-Christoffel conformal mapping: a) four Gauss points; b) eight Gauss points; c) twelve Gauss points; d) twelve Gauss points; e) twenty-four Gauss points. . These figures were created making use of the MATLAB SC Toolbox [1]. . . . . . . . . . . . . 60 4.19 Integration rule for a regular 4-nodes element: 1 Gauss point for the reduced rule and 4 Gauss points for the complete rule; b) Integration rule for an element crossed by a discontinuity: 1 Gauss point per triangle as a reduced rule and 3 Gauss points per triangle as a complete rule. . . . . 62 4.20 Test examples for single element containing a crack. . . . . . . . . . . . . 64 4.21 Cracked plate with respective boundary conditions. . . . . . . . . . . . . . 65 4.22 a) Coarse finite element mesh. b) Fine finite element mesh. . . . . . . . . 66 4.23 a) Complete integration rule; b) reduced integration rule. . . . . . . . . . 66 4.24 Deformed configurations for different integration schemes: a) complete integration rule for all the elements; b) reduced integration rules for all the elements; c) B-bar methodology applied only to the regular finite elements; d) B-bar methodology for both regular and enriched finite elements. 67 4.25 Deformed configurations for different integration schemes: a) complete integration rule for all the elements; b) reduced integration rules for all the elements; c) B-bar methodology applied only to the regular finite elements; d) B-bar methodology for both regular and enriched finite elements. 68 4.26 Crack opening for integration schemes c (line) and d (dashed line). . . . . 69 4.27 Cook’s membrane problem. . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.28 a) Example of a mesh with the crack modeled through enrichment. b) Example of a mesh containing a crack explicitly modeled. . . . . . . . . . 70 4.29 Vertical displacement of the top right corner node, as a function of the number of elements per side, for the explicitly meshed crack and enriched crack. ...................................... 71 4.30 Comparison of the vertical displacement of the top right corner node, as a function of the number of elements per side, for different tip enrichment functions. .................................... 72 4.31 Vertical displacement of the top right corner node, as a function of the number of elements per side, obtained for the explicitly meshed crack of figure 4.28(b), using the F-bar methodology and the Enhanced strain approach used by Dolbow et al. . . . . . . . . . . . . . . . . . . . . . . . 74 4.32 Comparison of the vertical displacement of the top right corner node, as a function of the number of elements per side, for explicitly meshed and enriched approximations using the F-bar methodology and the Enhanced Strain approach employed by Dolbow et al. . . . . . . . . . . . . . . . . . 74 4.33 a) Maximum value of the deviation from the incompressibility condition, (J−1), in the element crossed by the crack marked in the mesh b). . . . 75 4.34 Definition of interpolatory shape functions used in the variable transfer. . 76 List of Figures xvii 4.35 Group of Gauss points (in black) that influence a particular one. When a Gauss point is located near a discontinuity, only the points in the same side of the discontinuity is selected. . . . . . . . . . . . . . . . . . . . . . . 77 4.36 Group of Gauss points (in black) that influence a particular one, located nearacracktip. ................................ 77 5.1 Ductile fracture process: a) inclusions in a metallic matrix; b) void nucleation; c) void growth; d) strain localization and necking between voids; e) void coalescence and formation of a fracture surface. . . . . . . . . . . 80 5.2 Set of control points to define the interpolatory B-spline patch of the centralelement.................................. 83 5.3 Search for the point with highest damage at the boundary of the domain: a) Selection of the element with the highest damage value and b) selection of the adjacent segments and respective nodes. . . . . . . . . . . . . . . . 84 5.4 Ductile crack growth: a) void nucleation ahead of the crack tip; b) void growth; c) voids connect to the main crack causing its propagation. . . . 85 5.5 Selection of points to determine the crack growth direction. . . . . . . . . 86 5.6 Iterative scheme for equilibrium recovery during a crack growth step. . . . 90 5.7 Time-stepping scheme before and after crack initiation/propagation. . . . 91 5.8 Pre-crackedplate................................ 92 5.9 Mesh refinement for the plate with an initial crack: a) 20 ×31, b) 40 × 61 and c) 50 ×75elements. ......................... 93 5.10 Damage contours obtained for the local case, for an applied displacement of 0.05mm for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75 elements. 93 5.11 Damage contours obtained for lr= 0.8mm, for an applied displacement of 0.05mm, for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75 elements. 93 5.12 Damage contours obtained for lr= 1.6mm, for an applied displacement of 0.05mm, for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75 elements. 94 5.13 Crack path obtained using lr= 0.8 mm, for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75elements........................ 95 5.14 Crack path obtained using lr= 1.6 mm, for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75elements........................ 96 5.15 Reaction force in function of the applied displacement lr= 1.6mm. . . . . 97 5.16 Crack length evolution lr= 1.6mm....................... 98 5.17 Axisymmetric Specimen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.18 Mesh refinement for the axisymmetric notched specimen, mesh density of a) 11, b) 31 and c) 51 elements per side. . . . . . . . . . . . . . . . . . . . 100 5.19 Relation between the reaction force and the displacement applied to the top edge of the notched specimen for the different meshes, considering lr= 1.6 mm and DC= 0.2. ..........................101 5.20 Relation between the reaction force and the displacement applied to the top edge of the notched specimen for the different meshes, considering lr= 2.0 mm and DC= 0.5. ..........................101 5.21 Comparison between the force-displacement curves obtained for the intrinsic lengths lr= 1.6 mm and lr= 2.0 mm. Critical damage has the value DC= 0.5..................................102 5.22 Damage contour prior to crack insertion for mesh a) 11, b) 21 and c) 31 elementsperside.................................102 List of Figures xviii 5.23 Damage contour and crack, for an applied displacement of 0.9 mm, for mesh a) 11, b) 21 and c) 31 elements per side. . . . . . . . . . . . . . . . . 103 5.24 Crack length evolution considering lr= 2.0 mm and DC= 0.5. . . . . . . 103 5.25 Relation between the residual force and the displacement applied to the top edge of the notched specimen for the mesh with 31 nodes per side, considering lr= 2.0 mm and different values for DC. ............104 5.26 Plane strain specimen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.27 Mesh refinement for the plane strain specimen, mesh density of a) 11, b) 21, c) 31, d) 41 elements per side. . . . . . . . . . . . . . . . . . . . . . . . 105 5.28 Reaction force as a function of the applied displacement of the top nodes of the plane strain specimen. . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.29 Damage contour and final crack for mesh a) 11, b) and 21 c) 31 elements perside......................................106 5.30 Different crack growth steps for the mesh with 41 elements per side. . . . 107 5.31ShearSpecimen.................................107 5.32 FEM meshes used in the double notched specimen: a) mesh with 16 nodes per side, b) mesh with 23 nodes per side. . . . . . . . . . . . . . . . . . . 108 5.33 Crack path obtained for a) mesh a, b) mesh b. . . . . . . . . . . . . . . . 108 5.34 Damage contour and final crack for mesh a) and mesh b). . . . . . . . . . 109 6.1 a) Traction-free crack b) Cohesive crack. . . . . . . . . . . . . . . . . . . . 112 6.2 a) Body containing a cohesive crack. b) Detail of the cohesive crack . . . 113 6.3 a) Bilinear cohesive law. b) Polynomial cohesive law. . . . . . . . . . . . . 114 6.4 a) Linear cohesive law. b) Exponential cohesive law. . . . . . . . . . . . . 115 6.5 Changes in the cohesive response in function of a) Fnb) ω0.........117 6.6 a) Point pon the crack surface. b) Images of point p. ...........118 6.7 Gauss points used in the integration of a cohesive law. . . . . . . . . . . . 119 6.8 a) Plate with initial crack b) FEM mesh. . . . . . . . . . . . . . . . . . . 120 6.9 Influence of Fnin the response of the cracked plate for constant ω0= 0.5 mm−1.....................................121 6.10 Influence of ω0in the response of the cracked plate for constant a) Fn= 200 MPa b) Fn=100MPa. ..........................121 6.11 Influence of a cohesive law in the reaction force during crack propagation. 122 6.12 Influence of a cohesive law in the evolution of the crack length. . . . . . . 122 6.13 a) Plane strain specimen b) FEM mesh. . . . . . . . . . . . . . . . . . . . 123 6.14 Reaction force-applied displacement curves. . . . . . . . . . . . . . . . . . 124 6.15 Reaction force-applied displacement curves and respective influence of the cohesivelaw. ..................................124 6.16 Influence of the cohesive law in the evolution of the crack length. . . . . . 125 6.17 Reaction force-applied displacement curves and respective influence of the cohesivelaw. ..................................128 6.18 Strain energy as a function of the critical damage. Correlation factor: R2= 0.985. ...................................128 6.19 Strain energy as a function of the critical damage. Correlation factor: R2= 0.967. ...................................129 6.20 Reaction force-applied displacement curves and respective influence of the cohesivelaw. ..................................129 6.21 a) Plane strain specimen b) FEM mesh. . . . . . . . . . . . . . . . . . . . 130 List of Figures xix 6.22 Strain energy as a function of the critical damage. Correlation factor: R2= 0.999. ...................................131 6.23 Strain energy as a function of the critical damage. Correlation factor: R2= 0.999. ...................................131 6.24 Reaction force-applied displacement curves and respective influence of the cohesivelaw. ..................................132 List of Tables 3.1 Calculation of the deformation gradient by deriving the displacement field approximation. The SMSFreeze command denotes the manual introduction of an intermediate variable, which in this situation is necessary to correctly define the dependencies in the automatic derivation procedure. . 36 4.1 A possible basis for the incompressible deformation space of an element totally crossed by a discontinuity. . . . . . . . . . . . . . . . . . . . . . . . 53 4.2 A possible basis for the incompressible deformation space of an element containing a discontinuity tip. . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.3 A possible basis for the incompressible deformation space of an element totally crossed by a discontinuity. . . . . . . . . . . . . . . . . . . . . . . . 56 4.4 A possible basis for the incompressible deformation space of an element totally crossed by a discontinuity. . . . . . . . . . . . . . . . . . . . . . . . 57 4.5 Deformed configurations for different integration rules. . . . . . . . . . . . 65 4.6 Maximum value of the deviation from the incompressibility condition, max |J−1|.................................... 75 5.1 Material properties of the plate with an initial crack. . . . . . . . . . . . . 92 5.2 Material properties used in the numerical examples. . . . . . . . . . . . . 99 6.1 Material properties of the cracked plate. . . . . . . . . . . . . . . . . . . . 120 6.2 Strain energy for Dc= 0.4−0.5........................127 6.3 Strain energy for Dc= 0.6−0.7........................127 6.4 Material properties of the axisymmetric notched specimen. . . . . . . . . 131 xxi Dedicated to my family and friends xxiii Chapter 1 Introduction Ductile fracture has been regarded as a crucial issue in many industrial areas. Although it is usually associated to failure, as most of the metals, at room temperature, fail according to this mechanism, it may be an integrating part of production processes, such as machining, blanking or cutting. Numerous techniques have been developed to optimize operations and products, nevertheless, a deeper understanding of this phenomenon would lead to more efficient designs, which, in a first stage, prevents the use of unnecessary material and energy and, in a second stage, avoids catastrophic failures. More than ever, the responsible use of resources is essential, not only from the financial point of view but also for a sustainable world. Up to the nineteenth century the design of structures and mechanical parts was mostly empirical and understanding of material behaviour quite poor. Nevertheless, with the massive production of iron and steel demanded by the Industrial Revolution came the need of defining feasible design laws. Thus, by the beginning of the twentieth century, the study of ductile metals gained a mathematical approach with the pioneering works of Tresca [2], Huber [3], von Mises [4] and Hencky [5], where there was an attempt of developing material laws able to reproduce the experimental observations. The main outcome of these contributions was the prediction of the onset of plastic behaviour through the establishment of yield functions. In these approaches, the material properties are described at the meso-scale and explained by a few energy mechanisms. As the material is regarded as a continuous media, these works later gave birth to the so-called Continuum Mechanics. Although the dimensioning of structures using the yield criteria became quite popular in the engineering media between the World Wars, catastrophic failures of structures such 1 Chapter 1. Introduction 8 •M. Seabra, P. ˇ Suˇstariˇc, J. C´esar de S´a, T. Rodiˇc, Damage driven crack initiation and propagation in ductile metals using XFEM, Computational Mechanics (accepted for publication). •M. Seabra P. ˇ Suˇstariˇc, J. C´esar de S´a, T. Rodiˇc, Towards energetically consistent transition from damage to fracture using XFEM. (in preparation). •M. Seabra P. ˇ Suˇstariˇc, J. C´esar de S´a, T. Rodiˇc, Automatic generation of XFEM codes and application to crack propagation, Technical note (in preparation). Chapter 2 Standard Continuum Damage Mechanics Ductile metals, which include iron and steel, aluminium and copper alloys, are characterized by the presence of moderate to large plastic deformations prior to failure. Plastic strains may result from crystalline slip through dislocation motion or from degradation of the mechanical properties and therefore both mechanisms are strongly coupled in this type of materials. Classically, the material properties may be inferred from the study of the microstructure of the material, the so-called Micro Mechanical theories, or the material behaviour may be described in phenomenological way, based on a thermodynamic framework, referred as the Continuum Damage Mechanics (CDM). In the first category, the probably most extensively used model was proposed by Gurson [16] and later extended by Tvergaard and Needleman [18], which considers the material as a porous media, where voids nucleate and coalesce. On the other way, the model proposed by Lemaitre [25–27, 29] is the most representative of the CDM models, serving as a base for many other material laws. Here material progressive degradation is treated as an internal variable, acting in an averaging sense within a representative volume, corresponding to a scale where the material may be regarded as homogeneous. 9 Chapter 2. Standard Continuum Damage Mechanics 10 In this work, the Lemaitre model is used to describe the degradation of the mechanical properties of a solid prior to the formation of macro cracks, capturing the elastic, plastichardening and plastic-softening stages of material behaviour. Therefore, in this chapter the basic concepts of CDM and constitutive modelling are reviewed, along with the Lemaitre model for isotropic damage, in both local and non-local formulations. Rather then presenting an intensive research on CDM models, this chapter intends to provide the basis theory for the broader model for ductile fracture built throughout this thesis. Among the existing literature on the subject, interested readers may consult the books [62–64]. This chapter is organized as follows. In the first section the basic kinematics of deformable bodies is introduced, followed by the main strain and stress measures. Next the balance equations and the principles of thermodynamics are presented. In section 2.3 the constitutive theory is derived, in particular the Lemaitre model for ductile damage, making use of the previously introduced concepts. This chapter finalizes with the re-formulation of the Lemaitre model following a non-local integral approach and some brief conclusions. 2.1 Kinematics Kinematics provides the basic tools for motion and deformation description. In this section the fundamental quantities employed in the constitutive description of ductile metals are defined. Motion Let Bbe a deformable body occupying the domain Ω ⊂ R3, with the positions of its particles defined by X, referred as the reference configuration. The body may undergo amotion or deformation, within a certain time, corresponding to the transformation X:B→ R3 x=X(X, t) (2.1) where xis the position at the current configuration, at time t. Assuming that the motion Xis invertible, the position Xat a time tmay be recovered as: X=X−1(x, t) (2.2) Chapter 2. Standard Continuum Damage Mechanics 11 and the displacement field may be introduced as u(X, t) = X(X, t)−X(2.3) The time derivative of the function X(X, t), defined as V(X, t) = ∂X(X, t) ∂t (2.4) is the velocity of a particle at the position X, and consequently, the velocity of a particle at the current configuration may be determined by applying equation 2.2: v(x, t) = V(X−1(x, t), t) (2.5) Deformation gradient The derivative of the deformation, F, is a second order tensor termed deformation gradient, relating Xand xas follows: F(X, t) = ∂X(X, t) ∂X =∂x ∂X=∇XX(X, t) (2.6) where ∇Xis usually termed material gradient operator. Alternatively, the deformation gradient may also be defined in the current configuration, making use of the spatial gradient operator,∇xas F(x, t)=[∇xX−1(x, t)]−1(2.7) The displacement field may also be employed in the definition of the deformation gradient, leading to F=I+∇Xu(2.8) where Iis the second order identity tensor. From the definitions presented so far, it can be shown that line, surface and volume changes respectively transform according to: dx=FdX da=JF−TdA dv =JdV (2.9) In the above equations, J, termed Jacobian of the deformation, is the determinant of the deformation gradient J= det F(2.10) Chapter 2. Standard Continuum Damage Mechanics 12 It follows that a volume-preserving or isochoric deformation is characterized by: J= 1 (2.11) Moreover, any deformation can be locally decomposed into an isochoric deformation, Fvol and a volumetric deformation, Fiso, [65–67] as F=FisoFvol (2.12) where Fiso = (det F)−1 3F=J−1 3F(2.13) Fvol = (det F)1 3I=J1 3I(2.14) and Iis the second order identity tensor. Alternatively, a motion may be split in pure stretches and pure rotations, making use of the polar decomposition, that is, F=RU =V R (2.15) where Uand Vare unique, positive definite, symmetric tensors, termed, respectively right and left stretch tensors.Ris an unique proper orthogonal tensor, called rotation tensor. The tensors Uand Vmay also be introduced as: C=U2=FTF,b=V2=F F T(2.16) were Cand bare, respectively, the right and left Cauchy-Green strain tensors. Strain measures The deformation gradient is the fundamental kinematic quantity characterizing the changes of each material point during motion. Nevertheless, to quantify the relative distance change between two material points, that is, to characterize the straining, proper strain measures have to be defined. As the strains are not measurable quantities (such as the displacements) but a conceptual quantity to simplify analysis, there are numerous choices of strain tensors. An important Chapter 2. Standard Continuum Damage Mechanics 13 family of strain measures are the Lagrangian strain tensors: E(m)=(1 m(Um−I), m 6= 0 ln U, m = 0 (2.17) Similarly, another family of strain tensors, called the Eulerian strain tensors, may be obtained using the left stretch tensor: e(m)=(1 m(Vm−I), m 6= 0 ln V, m = 0 (2.18) Velocity gradient and strain rates The velocity gradient, l, is the spatial field defined as l=∇xv.(2.19) Alternatively, lmay be expressed as a function of the deformation gradient as l=∂ ∂t∂X ∂X ∂X ∂x=˙ F F −1(2.20) The velocity gradient may also be split into a symmetric and skew part, originating two important tensors as follows: d=1 2(l+lT) (2.21) w=1 2(l−lT) (2.22) The tensor dis called the rate of deformation and is associated with the straining, while the tensor lis termed spin tensor and is associated with rigid body velocity. Stress measures In line with the previous sections, for each type of strain measure, it is possible to define a stress measure. One of the most important stress tensors is the Cauchy or true stress tensor, σ, as: t=σn (2.23) where tis the surface traction and nis the associated outward normal vector. The Cauchy stress tensor may be split in a deviatoric,s, and a pressure,p, contributions, Chapter 2. Standard Continuum Damage Mechanics 14 as follows σ=s+pI(2.24) where p=1 3tr[σ] (2.25) Another stress measure frequently used is the Kirchhoff stress tensor, given by τ=Jσ(2.26) Other convenient stress tensors are the First and the Second Piola-Kirchhoff tensors, respectively introduced as P=τF −T(2.27) and S=F−1τF −T(2.28) 2.2 Conservation principles In the previous sections, the quantities used in the mathematical description of motion, straining and stress were introduced. In this section, this concepts are related in some fundamental conservation principles, which are essential to formulate a continuum mechanics problem. The detailed derivation of this principles is out of the scope of this thesis and may be found, for instance, in the references [62–64, 68]. Conservation of mass The first principle regards the mass conservation and may be expressed as ˙ρ+ρdivx˙ u= 0 (2.29) where ρis the mass density at the deformed configuration. Chapter 2. Standard Continuum Damage Mechanics 15 Momentum balance The momentum balance equation describes the equilibrium between internal and external forces in a body, as follows divxσ+b=ρ¨ u(2.30) where bdenotes the body force vector at the deformed configuration. These equation is frequently referred as strong form of the equilibrium equation and has to fulfil the following boundary condition: ¯ t=σn (2.31) where ¯ tis a traction vector applied to the boundary of the body. First and second laws of thermodynamics The first law of thermodynamics postulates that the energy must be conserved and may be mathematically expressed as: ρ˙e=σ:d+ρr −divxq(2.32) where erepresents the specific internal energy, rrepresents the density of heat production and qthe heat flux. Throughout this thesis, only purely mechanical processes will be addressed, thus equation 2.32 reduces to ρ˙e=σ:d(2.33) This equation may be rewritten in terms of the Kirchhoff stress tensor, as follows ¯ρ˙e=τ:d(2.34) where ¯ρ=Jρ is the reference mass density. Second law of thermodynamics refers to the irreversibility of the entropy production, which may be expressed as ρT ˙s+ divxq−ρr ≥0 (2.35) where srepresents the entropy and Tthe temperature. Similarly to the first law, for isothermal processes equation 2.35 reduces to ρT ˙s≥0 (2.36) Chapter 2. Standard Continuum Damage Mechanics 16 Combining the first and the second laws and introducing the Helmholtz free energy concept, ψ,defined as ψ=e−Ts (2.37) the so-called Clausius-Duhem inequality is obtained τ:d−¯ρ˙ ψ≥0 (2.38) 2.3 Finite strain damage mechanics The laws presented are valid independently of the material considered. Therefore, to fully characterise the response of a solid body and distinguish between different materials, a constitutive model has to be formulated and boundary conditions have to be prescribed. In this thesis, the chosen constitutive model, the Lemaitre model for ductile damage [25–27, 29], relies in the thermodynamics with internal variables, which means that at any instant, the thermodynamic state at any point is completely defined by the values of a given set of state variables. In order to successfully capture all the stages of material behaviour, the choice of the set of internal variables is crucial. In the particular case of ductile materials, the material degradation, that is damage occurs simultaneously with plastic deformations, and therefore, the internal energy function should reflect this dependency, as it will be shown throughout this section. 2.3.1 The damage variable Failure in ductile metals initiates with the nucleation cavities and micro-cracks, which may grow and coalesce, eventually leading to macro-crack generation and propagation [27, 69]. The quantification of the micro-defects at the meso-scale level, may be taken in account in a averaging sense, defining the damage variable, Dˆn, as: Dˆn=∂SD ∂S (2.39) where ∂S is the area of intersection of a plane, with normal ˆn, with the representative volume element (RVE) and ∂SDis the area that effectively intersects cavities or microcracks, as can be observed in figure 2.1. The RVE of a certain material is defined in such a way that all the properties are represented by homogenized variables (around 0.1mm3 Chapter 2. Standard Continuum Damage Mechanics 17 ∂S ∂SD ˆn Figure 2.1: Representative volume element. for metals [27]). When only ductile damage is considered, the micro-defects density is often a weak function of the plane orientation. Therefore, Dmay be considered isotropic and represented by a scalar variable, ranging from 0 to 1, corresponding to undamaged to fully damage material, respectively, which is the case of this thesis. Nevertheless, a second or fourth order tensor representation may be adopted for other applications [27, 70]. 2.3.2 Effective stress concept and strain equivalence principle In the previous sections some strain and stress measures were introduced to be used in the constitutive equations that characterize a given material. These quantities were developed for plain material, however, in the presence of damage the resistant area is smaller and an the effective Cauchy stress tensor may be introduced as: ¯ σ=σ 1−D(2.40) where σis the Cauchy stress tensor for undamaged material. From this definition it follows the strain equivalence principle, which states that a damaged material is governed by the constitutive laws of the plain material, replacing the true stress by the effective stress [27]. 2.3.3 Multiplicative plasticity framework In ductile materials damage and plasticity are strongly connected, therefore the concepts associated with finite strain elasto-plasticity are now introduced. Chapter 3 Finite Element Distretization and Automatic Code Generation The description of ductile materials following the damage model described in chapter 2 results in a complex system of differential equations. As only a very restricted set of analytical solutions are available, the use of numerical methods becomes imperative. In particular, the Finite Element Method (FEM) stands out as the most widely used numerical method for solving solid mechanics problems for its effective predictive capabilities. The basis of the FEM are well established in the literature as, for instance, in the reference works of Zienkiewicz [36], Hinton and Owen [80], Hughes [81] or Crisfield [82]. Hence the method will be briefly reviewed in the first part of this chapter, which initiates by re-formulating the balance laws presented in chapter 2 in order to define a Initial Boundary Value Problem, suitable to be solved with the FEM. Next, the displacementbased discretisation is introduced, along with the non-linear incremental finite element equation. The second part of this chapter will focus on the automatic generation of FEM codes. The computer science advances, in both software and hardware, had a significant impact in generalizing the use of the FEM, allowing the simulation of problems with higher number of variables, described by increasingly complex mathematical models . Moreover, the combination of modern Computer Algebra Systems with an automatic derivation technique, allows a faster and more efficient implementation of FEM formulations. In this chapter, special attention is paid to one of these systems, the AceGen system, which was used to develop a considerable amount of computer codes used in the context of 25 Chapter 3. Finite Element Distretization and Automatic Code Generation 26 this thesis. 3.1 The quasi-static IBVP In Chapter 2 the fundamental equilibrium laws were presented. However, to solve a problem using the FEM, these laws are often formulated as variational principles. Although this methodology is not mandatory, approximate solutions are often obtained through the weak forms of field equations. Therefore, recalling the momentum balance equation: divxσ+b=ρ¨ u(3.1) and considering the quasi-static case, as throughout this thesis the inertial effects will be neglected, we may write: divxσ=0(3.2) Equation 3.2 may be multiplied by a virtual displacement, η, and integrated over the domain, as: ZX(Ω) (divxσ)Tηdv = 0 (3.3) After some algebra manipulation and making use of the divergence theorem, equation 3.3 may be written in the following form: ZX(Ω) σ:∇xηdv −ZX(∂Ω) (σ.n)Tηda = 0 (3.4) As introduced in chapter 3, the following boundary condition has to be fulfilled together with equation 3.1 ¯ t=σn (3.5) where ¯ tis a prescribed traction at the boundary of the body. Therefore, replacing this relation in equation 3.4, the weak form of the equilibrium equation of obtained: ZX(Ω) σ:∇xηdv −ZX(∂Ω) ¯ t·ηda = 0 (3.6) The weak form of the equilibrium equation is equivalent to the strong form and allows the definition of the weak form of the Initial Boundary Value Problem (IBVP) associated with a certain deformation process, as follows. Find a kinematically admissible displacement function, u∈ K such that the equation Chapter 3. Finite Element Distretization and Automatic Code Generation 27 ZX(Ω) σ:∇xηdv −ZX(∂Ω) ¯ t·ηda = 0 (3.7) is satisfied for all t∈[t0, tn]and for all η∈ Vt.The set of virtual displacements and the set of kinematically admissible displacements, at time tare respectively defined as Vt={η: Ω −→ X(Ω)|η=0∈∂Ω(t)}(3.8) K={u: Ω −→ X(Ω)|¯ u(X, t), t ∈[t0, tn],X∈∂Ω}(3.9) The Cauchy stress, at each point of the body is expressed as σ(t) = σ(F(t),α(t)) (3.10) where αrepresents the set of internal variables associated with the material and Fis a given prescribed deformation gradient history, defined as F(t) = I+∇Xu(X, t) (3.11) The solution of the IBVP reflects how a solid body will behave when subjected to certain boundary conditions. Material description Equation 3.6 reflects the spatial version of the Principle of Virtual Work, nevertheless, depending on the particular problem, the material version may provide a more efficient implementation. The first member of equation 3.6 is termed internal work,Wint, and may be written in the reference configuration, recalling the definition of the first Piola- Kirchhoff stress tensor: P=JσF −T(3.12) resulting in Wint =ZΩ P:∇XηdV (3.13) The first Piola-Kirchoff tensor is not completely related to the material configuration but it may be replaced by the second Piola-Kirchoff strain tensor as: P:∇Xη=S:FT∇Xη=S: (FT∇Xη+∇XηTF) = S:∂E(3.14) Wint =ZΩ S:∂EdV (3.15) Chapter 3. Finite Element Distretization and Automatic Code Generation 28 where ∂Eis the first variation of the Green-Lagrange strain tensor. The second member of equation 3.6, the external work,Wext may be written in the material description simply by taking the applied traction in the initial configuration, ¯ t0 Wext =ZX(∂Ω) ¯ t.ηda =Z∂Ω ¯ t0.ηdA (3.16) Either in the material or spatial version, the IBVP is usually non-linear and therefore linearisation and discretisation are necessary to produce accurate solutions using the FEM, as described in the upcoming sections. 3.2 Displacement-based finite elements The FEM for the numerical solution of the IBVP problem described in the previous section consists in replacing the sets Vtand K, defined in equations 3.8 and 3.9, respectively, by discrete subsets, Vh tand Khobtained from a finite discretization of the domain. In the case of displacement-based finite elements, the interpolated field variables are the displacements, which within a given element, e, assume the form: u(x) = nnode X i=1 N(e) i(x)ui(3.17) where N(e) iis the shape or interpolation function associated with node iand nnode is the number of nodes of the element. Similarly, the interpolation function may be defined over the entire approximate domain, which is constituted by a total number of npoin nodal points, as u(x) = npoin X i=1 Ng i(x)ui(3.18) In equation 3.18, urepresents the global vector of nodal displacements, which, in a problem of dimension ndim, is given by uh= [u1 1, . . . , u1 ndim, . . . , unpoin 1, . . . , unpoin ndim ] (3.19) and Ngis the global interpolation matrix defined as: Ng i= [ diag[Ng 1(x)] diag[Ng 2(x)] ··· diag[Ng npoin(x)] ] (3.20) Chapter 3. Finite Element Distretization and Automatic Code Generation 29 where diag[Ng 1] represents a ndim ×ndim diagonal matrix as follows diag[Ng 1] =        Ng i0··· 0 0Ng i··· 0 . . .. . ..... . . 0 0 ··· Ng i        (3.21) Equation 3.18 may be reformulated as: uh(x) = Ngu(3.22) Similarly, the field of virtual displacements may be written as: ηh(x) = Ngη(3.23) Using the presented notation the discretized sets Vh tand Khmay be defined as: Vh t=ηh(x) = npoin X i=1 Ng iηi|ηi=0if xi∈∂Ω(3.24) Kh=uh(x) = npoin X i=1 Ng iui|ui=¯ u(xi) if xi∈∂Ω(3.25) Now, equation 3.6 may be replaced by its discretized version, as follows ZX(Ωh) σ:∇x(Ngηh)dv −ZX(∂Ωh) ¯ t·(Ngηh)da = 0 (3.26) In classical FEM implementations it is usual to define the global discrete symmetric gradient matrix,B, which for a plain strain problem has the form: Bg=    Ng 1,10Ng 2,10··· Ng npoin,10 0 0Ng 1,20Ng 2,2··· 0Ng npoin,2 Ng 1,2Ng 1,1Ng 2,2Ng 2,1··· Ng npoin,2Ng npoin,1     (3.27) where the following notation was employed (·)i,j =∂(·)i ∂xj (3.28) Finally, the discretized virtual work expression can be rearranged to ZX(Ωh) [σTBgηh−b·Ngηh]dv −ZX(∂Ωh) ¯ t·Ngηhda = 0,∀ηh∈ Vh t(3.29) Chapter 3. Finite Element Distretization and Automatic Code Generation 30 For completeness, when using the material version of the virtual work equation, the discretised form of the variation of the Green-Lagrange strain tensor is provided herein: ∂E=1 2 npoin X i=1 [FT(ηh i⊗∇XNi)+(∇XNi⊗ηh i)F)] (3.30) where F= npoin X i=1 xi⊗∇XNi(3.31) The non-linear incremental finite element equation The equilibrium equation is, in general, non-linear. For the particular case of material non-linearities, such as in the Lemaitre model where the Cauchy stress is dependent of the history of strains to which the solid has been subjected, a suitable temporal discretisation is required. In this thesis, a pseudo-time discretisation between the time increments [tn, tn+1] will be considered, in a fully implicit scheme. The function, ˆ σ, defined as σn+1 =ˆ σ(Fn+1,αn) (3.32) is assumed to exist and is associated with an integration algorithm that delivers the behaviour for a given deformation gradient, Fn+1, and a set of internal variables, αn, which are assumed to be constant within the one increment. The Newton-Raphson algorithm is particularly attractive for the solution of this type of problems due to its robustness and quadratic rates of asymptotic convergence, and thus was used in this work. During a solution procedure in which equilibrium is not yet satisfied, there is a residual Rbetween the internal work and the external work. As in the problems considered in this thesis, the applied traction is independent of the deformation, the individual element contribution to the residual may be defined as: Re(un+1) = Ren+1 =ZΩe Sn+1 ∂En+1 ∂uen+1 dΩ (3.33) where Ωeis the domain of the element. Following the general procedures of the Newton method, a system of equations of the the form Rn+1 =0(3.34) is obtained. However equation 3.34 is not satisfied unless convergence has occurred. Since the current global residual Rn+1 depends on the global displacement vector of the Chapter 3. Finite Element Distretization and Automatic Code Generation 31 previous time step, un, and on the current global displacement vector, un+1, that is, Rn+1 =Rn+1(un+1,un) (3.35) the (k)th iteration step may be written as Kk Tn+1∆uk n+Rk n+1 =0(3.36) where KTis the global tangent stiffness matrix given by Kk Tn+1 =∂Rk n+1 ∂uk n+1 (3.37) Finally, the displacements are updated as follows: uk+1 n+1 =uk n+1 + ∆uk n+1 (3.38) Numerical integration In a finite element implementation, the exact integrals are replaced by a numerical integration procedure. In general, standard Gaussian quadrature will be used for this purpose. The integral of a generic function fover a domain Γ is given by ZΓ f(ξ)dξ ≈ ngaus X i=1 ωif(ξi)(3.39) where ngaus is the number of approximation points, with coordinates ξiand weights ωi. At this stage, it should be remarked that in the following chapters of this thesis, some alternative methods for numerical integration will be addressed, in particular in the presence of discontinuities. 3.3 Automatic code generation The implementation of a FEM code may be quite complicated and time consuming, especially when dealing with non-linear material models. In general, the discretisation and linearisation procedures described in the previous sections have to be performed prior to the actual code writing. The complexity of the required calculations increases the probability of error and, consequently, inefficient codes may result. Chapter 3. Finite Element Distretization and Automatic Code Generation 32 From another perspective, a new formulation represents essentially the coding of a different element kinematics and/or a different material model, which suggests that a large number of operations related to FEM programming could be automated. In this context, the use of automatic differentiation tools and modern symbolic and algebraic computer systems may represent a considerable advantage. However, in most of these systems, the classical symbolic derivation leads to the exponential growth of expressions and inefficient codes. In order to overcome this limitation, the AceGen System [83] employs the Simultaneous Stochastic Simplification of Numerical Code approach [84, 85], which combines an automatic differentiation technique with the simultaneous optimization of symbolic expressions, allowing the automatic generation of efficient numerical codes. In the following subsections the main features of this system will be described and the application to the FEM will be illustrated. 3.3.1 Automatic code generation features For an efficient code generation, the AceGen system gathers some important features characterized in the following paragraphs. Differentiation As shown in the beginning of this chapter, the development of a finite element formulation involves the calculation of numerous derivatives. Therefore a proper automatic differentiation tool, in which the dependence on intermediate quantities is correctly accounted for is essential. The derivative of a general function, f, which depends on a set of mutually independent variables aand a set of mutually independent intermediate variables b(that is, bdepend on a) may be expressed as follows: ∇f=∂f(a,b(a)) ∂(a)∂(b) ∂(a)=M(3.40) The total derivatives of the variables bwith respect to the variables amay be written in the matrix form, where Mis the respective Jacobian matrix. During the differentiation procedure, when there exists an explicit algorithmic dependence between band a, the derivatives may be obtained automatically by the chain rule without user intervention, Chapter 3. Finite Element Distretization and Automatic Code Generation 33 simply by replacing its value for the corresponding entrance of matrix M. Nevertheless, for some situations, the direct application of the chain rule leads to erroneous results. In particular, when function fonly depends explicitly on b, ∇f=∂f(b) ∂(a)∂(b) ∂(a)=M(3.41) or the dependency between band ahas to be neglected, ∇f=∂f(a,b(a)) ∂(a)∂(b) ∂(a)=0(3.42) a proper exception has to be considered by the automatic differentiation tool [84, 85]. A typical example of the first case is a differentiation that involves a transformation of coordinates, for instance, from a reference isoparametric finite element, to a global coordinate finite element. The introduction of the additional condition for the proper derivation of the deformation gradient is illustrated in section 3.3.2. Simplification of large mathematical expressions After setting up a proper differentiation tool, the methodology to simplify mathematical expressions is of major importance. Most of the existing symbolic systems, such as Mathematica or Matlab, only search for common sub-expressions after all the formulae have been derived, using a pattern-matching technique. However this methodology is insufficient to obtain efficient codes for calculating typical FEM quantities like the internal force vector or the stiffness matrix, due to the excessive swell of the mathematical expressions [86, 87]. To overcome this problem, AceGen employs the Simultaneous Stochastic Simplification of numerical code [84, 85]. In this methodology the search of common sub-expressions is performed after each automatic derivation step, optimizing the introduction of intermediate variables. In addition, the equivalence between expressions is not only searched symbolically but also numerically [88]. It is a fact that the correctness of the simplified expressions can only be determined with a certain probability, nevertheless, when dealing with relatively smooth functions, as it is the case of the FEM, this can be neglected. Chapter 4. The Extended Finite Element Method 40 4.1 Formulation and basic numerical implementation 4.1.1 Displacement approximation and choice of the enrichment functions The XFEM is a methodology which allows for the insertion of special features, such as discontinuities or interfaces, in a FEM problem, independently from the mesh, by enriching the displacement field, u, with one or more functions ψ, as follows: u(x) = n X i=1 Ni(x)ui+ nenr X j=1 Nj(x)ψ(x)ψj(4.1) where Nirepresents the element shape functions and ψjare the extra degrees of freedom associated to the function ψin each one of the nenr nodes. The function ψ, termed enrichment function, contains the description of the desired feature and satisfies the partition of unity principle, which may be stated as, X j Nj(x)ψ(x) = 1.(4.2) It should be noticed that considering Njas the standard FEM shape functions, this property is automatically satisfied, though this choice is not mandatory. The mathematical expression of ψvaries greatly, as this method may be applied to very diverse scientific areas. In particular case of this thesis the XFEM is used in the context of fracture and therefore the enrichment function should be capable of capturing the decohesion between the two surfaces of a crack. Discontinuous functions such as the following Heaviside function, H(ˆη), defined in terms of the coordinate ˆηperpendicular to the crack plane (figure 4.1), are particularly suitable for this purpose. H(ˆη) = (1: ˆη≥0 −1: ˆη < 0(4.3) The displacement approximation then reads: u= n X i=1 Niui+ nsplit X j=1 NjHaj(4.4) where ajare the degrees of freedom associated to the Heaviside function. When a finite element is totally cut by a crack is usually termed split element and therefore, has a number of nsplit enriched nodes. An element only partially crossed by a crack, that is, Chapter 4. The Extended Finite Element Method 41 ˆη (a) ˆη (b) Figure 4.1: Definition of the coordinate perpendicular to the crack in a) 2D problem b) 3D problem. ˆ ξ ˆη η ξ α lc ξ= 1ξ=−1 η= 1 η=−1 Figure 4.2: Crack tip coordinates for a 2D element or a cut along ˆ ζ= 0 for a 3D element. containing a crack tip is named tip element. A local coordinate system, ˆ ξ= (ˆ ξ, ˆη, ˆ ζ), centered at the crack tip, may be defined for these type of elements, as illustrated in figure 4.2, for a plane element or a cut along ˆ ζ= 0 for a 3D element. In 2D problems ˆ ξ= (ˆ ξ, ˆη), where ˆ ξis in the crack direction and ˆηis perpendicular to the crack line. In 3D problems ˆ ξ= (ˆ ξ, ˆη, ˆ ζ), in which ˆ ξand ˆ ζare in the crack plane, with ˆ ξ= 0 defining the crack tip line at the crack plane, and ˆηcoordinate is normal to the crack plane. Applying directly the displacement approximation of equation 4.4, in general u(xi)6= ui. This problem may be easily overcome by shifting the approximation around the node of interest, as follows: u(x) = n X i=1 Niui+ nsplit X j=1 Nj[H(x)−H(xj)]aj(4.5) Moreover, the shifted version of equation 4.5 prevents the enrichment of spreading to the neighbour elements, which contain standard and enriched nodes (figure 4.3), usually designated blending elements. Chapter 4. The Extended Finite Element Method 42 Figure 4.3: Crack inserted in the mesh through the XFEM. The nodes belonging to split elements are marked with circles, while the nodes belonging to tip elements are marked with triangles. Elements containing both enriched and standard nodes are designated blending elements. Enriching finite elements only with the Heaviside function always causes the extension of the crack to the edge of an element. Nevertheless, to represent cracks truly mesh independent, it is desirable to define a strategy to deal with the closure of a crack away from the element edge, as can be observed in figure 4.3. Elements containing crack tips, may then be enriched by a functio, which has to ensure that the enrichment vanishes exactly at the crack tip. Some alternative definitions for this function are proposed and presented below, in terms of crack tip coordinates and the crack length inside the tip element, lc. The first proposed form for the function, here named as R, is a Heaviside type function, defined as: R(ˆ ξ) = (1 for ˆ ξ≤0 0 otherwise (4.6) and determines that only Gauss points located behind the crack front should be enriched. The second proposed form is a linear type function, defined as: R(ˆ ξ) = (−ˆ ξ lcfor ˆ ξ≤0 0 otherwise (4.7) A higher order function, defined in reference [96], is also adopted here as: R(ˆ ξ) = (3( ˆ ξ lc)2+ 2( ˆ ξ lc)3for ˆ ξ≤0 0 otherwise (4.8) In this work, these functions are considered to be constant along the ˆ ζcoordinate, in the case of 3D elements. Chapter 4. The Extended Finite Element Method 43 Besides ensuring that the enrichment disappears exactly at the crack tip, the three forms of the Rfunction are particularly suitable for ductile fracture problems, where plasticity is not only confined to the region around the crack tip but widely spread. Furthermore, as in the case of the Heaviside function, these polynomial enrichments vanish outside the elements containing discontinuities, which is not the case if other types of enrichments are employed [94, 97, 98]. The full displacement approximation may then be written as: u(x) = n X i=1 Niui+ nsplit X j=1 Nj[H(x)−H(xj)]aj + ntip X k=1 Nk[R(x)−R(xk)][H(x)−H(xk)]bk (4.9) Remark 4.1.In the framework of linear elastic fracture mechanics, the most used enrichment functions to deal with crack tips are based in the asymptotic displacement fields derived from analytical solutions [54, 90, 99]. Therefore, the following set of four functions, ψi tip, is employed: ψi tip ={√rsin(θ 2),√rsin(θ 2) sin θ, √rcos(θ 2),√rcos(θ 2) sin θ}(4.10) These functions are widely used for brittle fracture problems, as they reproduce a singularity at the crack tip. Nevertheless, as previously referred, in a context of ductile fracture they may not be adequate. Finally, for a successful implementation of the method it is essential to determine the split nodes and the tip nodes, which may be done through the Level Set Method, described in the next section. 4.1.2 The level set method The Level Set Method [100, 101] is the most frequently used complement of the XFEM to accurately determine the position of the discontinuities, that is, the set of nodes to be enriched [102, 103]. Chapter 4. The Extended Finite Element Method 44 Ω Ω− Γ Ω+ d x xΓ Figure 4.4: Domain Considering a domain Ω, divided in the regions Ω+and Ω−by the interface Γ, the level set function, φ(x), is defined in a such a way that φ(x) = (1: x∈Ω+ 0: x∈Γ −1: x∈Ω− (4.11) A common choice for φ(x) is the signed distance function, defined as φ(x) = ±d=±kx−xΓk(4.12) where dis the distance between point xand the closest point to xlaying in the interface, xΓ, as depicted in figure 4.4. The sign of ddepends on which side of the interface is located point x. The function φ(x) determines if a point lies above or bellow a crack. Nevertheless, as usually a crack does not cut the entire domain, a second function,ϕ(x), has to be introduced to deal with the active crack tips. When a crack is defined by straight line segments, which is the case in this thesis, a function ϕi(x) may be defined by each active crack tip as follows: ϕ(x) = (x−xi).ˆ ti(4.13) where ˆ tis a unit vector tangent to the latest crack segment and xiis the location of the ith crack tip, as can be observed in figure 4.5. Although, functions φand ϕihave to be updated at each crack growth step, the search for newly enriched elements is only required around the previous crack tip. Furthermore, in a ductile fracture problem the direction and magnitude of a crack increment are determined by the material model, as it will be described in the following chapters. Chapter 4. The Extended Finite Element Method 45 φ < 0 φ > 0 ϕ2<0 ϕ2>0 ϕ1>0 ϕ1<0 ˆ t1 ˆ t2 Figure 4.5: Level set functions. Amax Amin (a) (b) (c) Figure 4.6: The crack is very close to a) a node b) an edge and totally contained in one element in c), possibly leading to ill-conditioned stiffness matrices. Consequently, the level set functions are updated by simple geometrical considerations. The level set function described may also describe crack propagation for planar cracks in 3D problems. A more detailed description of the procedure may be found in the references [104–107]. 4.1.3 Basic implementation The implementation of an XFEM formulation consists fundamentally in modifying the standard FEM displacement approximation and create a structure to deal with the extra degrees of freedom, introduced by the enrichment functions, when performing tasks such as the assembly of the internal force or the stiffness matrix. Some linear dependencies or ill-conditioning may arise due, essentially, to elements where the ratio of the areas/volumes between the two sides of the crack is very high, Amax/Amin, as illustrated in figures 4.6(a) and 4.6(b), or when the discontinuity lies only inside one element (figure 4.6(c)). The first situation is usually prevented by setting a tolerance from which the crack should be converted into a inter-element crack while in the second situation, the mesh may be refined [94, 108]. Chapter 4. The Extended Finite Element Method 46 Another key feature in the implementation of the XFEM is numerical integration. Due to its high importance the next section is fully devoted to this subject. Other features that had to be developed to deal with the use of the XFEM with nonlinear problems, such as transfer of history variables and non-locality, will be also addressed in the end of this chapter. 4.2 Numerical integration and the incompressibility issue Numerical integration has to be carefully addressed when implementing a XFEM code, as the results produced by regular Gaussian quadrature rules are not accurate enough, due to the presence of discontinuous functions. Two main techniques were developed since the early stages of XFEM. Probably the most common of those techniques consists of dividing the elements containing discontinuities in sub cells conforming the discontinuity, usually triangles or tetrahedrons, for 2D and 3D problems respectively, and applying a regular rule to each sub cell [90, 91, 109]. Alternatively, some authors [110–112] prefer to keep fixed Gauss point positions, but largely increase their number inside an element to compensate for not adjusting the discontinuity pattern. As computationally this last method may become extremely expensive, other approaches were created to suppress the need of element subdivision, starting with the works of Iarve [113] and Ventura [114] who derived continuous polynomial functions which, integrated over an element, reproduce the integrals of the discontinuous XFEM functions. This method was extended for regularized discontinuities by Benvenuti et al. [115], but it is only exact for triangular or tetrahedral elements. Another integration rule was proposed by Holdych [116] for the same type of elements, where the drawback of having fixed positions for the Gauss points is compensated by adopting variable weights. The relation between these last two methods was established by Belytschko et al. [117] and, subsequently, the same authors developed a boundary integration method [118], applicable to quadrangular elements, in which all nodes are equally enriched. Recently, Natarajan et al. [119] proposed an integration technique based on the Schwarz-Christoffel conformal mapping valid for any type of 2D finite elements, and Mousavi et al. [120] made use of the node elimination technique to produce accurate integration rules with a reduced number of Gauss points. As, so far, most of the XFEM applications are within the scope of brittle fracture, the accuracy of the methodologies described above is typically measured by the ability to reproduce the stress intensity factors, for particular problems that have an available analytical solution. Nevertheless, the direct extension of those methodologies to problems Chapter 4. The Extended Finite Element Method 47 involving ductile fracture may not be adequate. In fact, when ductile fracture is the issue, plastic deformation is widely spread and not just confined to the nearby zones around the crack, as it happens in brittle fracture. Furthermore it is well known that incompressibility of plastic deformation may lead to locking of the numerical solution, especially when low order finite elements are employed [121–126]. Consequently, in order to avoid it, adequate formulations or integration techniques must be used. Depending on the problem addressed solutions may include, for example, selective reduced integration techniques [122], mixed formulations [65], B-bar methods [123], F-bar methods [63, 127] or enhanced assumed strain methods [128], to name a few. Therefore when XFEM is used to model ductile fracture some questions emerge as important topics, namely those related with the discussion on the adequate number of Gauss points that should be used in the sub cells, when an element contains or is crossed by a crack, and the ability to satisfy the incompressibility constraints. In the following sections, having in view the application of XFEM to ductile fracture problems, the efficiency of different integration rules is investigated, focusing on their ability to satisfy volumetric incompressible constraints in the context of both infinitesimal and finite strains. As a consequence, the conventional formulations for B-bar and F-bar are adapted to include the XFEM enrichment functions. Other alternatives to deal with incompressibility in XFEM elements include the adaptation of the assumed enhanced strain method by Dolbow et al [129] and the mixed formulation technique by Legrain et al [130]. Here, the efficiency of the proposed techniques is evaluated with some numerical examples and the performance of the F-bar method for XFEM is compared with the method developed by Dolbow et al [129]. 4.2.1 Methodology to assess the ability to alleviate locking The methodology here adopted to assess the tendency or the ability of a particular method to depict or to alleviate locking will follow closely some previous work that may be found in references [124–126]. This methodology is based on the analysis of the underlying sub-space of incompressible modes embedded in the approximation method adopted, rather than relying on mere constraint indexes which may be misleading, as they are defined independently of the boundary conditions. The aforementioned analysis may give information on how the dimension of this subspace can be enlarged or reduced as a function of the integration rule adopted. Consequently it may give precious hints on the ability of a particular method to reproduce incompressible deformations, on how Chapter 4. The Extended Finite Element Method 48 spurious modes may be prevented or on how over integration rules may be avoided. Recalling equation 2.11 from chapter 2, an isochoric deformation is characterized by: J= 1 (4.14) Taking into account that the material time derivative of the volume ratio may be written as a function of the rate of deformation tensor or of the velocity field: ˙ J=Jtrd=Jdivv(4.15) and considering equation 4.14 and its time derivative, the incompressibility condition may alternatively be stated as: trd= 0 (4.16) For a small strain problem, disregarding second and higher order terms, this condition may be simplified to: divu= 0 (4.17) where urepresents the displacement field. Consequently, when imposing the incompressibility condition in a FEM framework via the Penalty Method or the Lagrange Multiplier technique [124], a local discretized form of equation 4.16 or equation 4.17 may be established at each numerical integration point of a particular element with the form Qwh= 0 (4.18) where whis the vector of element nodal variables from which the rate of deformation tensor or the displacement field are approximated and Qis a matrix i×n, in which iis the number of integration points and nthe number of unknowns in wh. Consequently, element wise, an admissible incompressible solution whshould belong to the null space of matrix Q, which defines locally the space of incompressible deformations, Ih Ih={wh∈U:Qwh= 0}(4.19) Locking occurs, when for a given set of boundary conditions, the expected solution cannot be reproduced by a linear combination of a given basis of Ih[124–126]. Having applications on ductile fracture problems, in which extensive plastic deformations may be present, the analysis of this underlying sub-space and how it may be affected by adopting or altering a certain integration rule will be addressed in the following sections of this chapter. Chapter 4. The Extended Finite Element Method 49 4.2.2 Incompressible deformation modes in elements containing XFEM enrichments A basis of Ih(equation 4.19), the space of admissible incompressible deformations, must be determined to study the vulnerability of the enriched elements to volumetric locking. It should be noticed that for the XFEM approximation, the nodal variables whused in the approximation of the strain or the strain rate encompasses both the displacements uhand the extra variables ahor bhor their pseudo-time derivatives. Therefore, for simplicity and in line with the previous works of Cesar de Sa et al. [124–126], a basis of Ihof a particular element, in which global coordinates xand natural coordinates ξare coincident, may be determined by obtaining a basis of the null space of matrix Q. For this particular case equation 4.19 reads: Qwh= [ ∂Ni ∂ξ ∂Ni ∂η ∂Ni ∂ζ ∂(Niψ) ∂ξ ∂(Niψ) ∂η ∂(Niψ) ∂ζ |i=nnodes ]{uh ψh}= 0 (4.20) where the number of lines in Qequals the number of Gauss points, nGauss, used in the finite element, ψrepresents the enrichment function associated with the extra variables ψh, which may be ahor bhwhether it is a split element or a tip element. To determine the admissible incompressible deformations of an element containing a discontinuity, equation 4.20 has to be particularized for the displacement approximations of equation 4.5. For an element totally crossed by a discontinuity, only the Heaviside enrichment is considered leading to the displacement approximation, expressed locally as: u(ξ) = n X i=1 Niui+ nsplit X j=1 NjHSh iaj(4.21) Therefore, the terms associated with the function Rvanish and each line, g, of matrix Qmay be written as: Qg={∂Ni ∂ξ ∂Ni ∂η ∂Ni ∂ζ |i=nnodes ∂Ni ∂ξ HSh i∂Ni ∂η HSh i∂Ni ∂ζ HSh i|i=nnodes }(4.22) where HSh idenotes the shifted version of the Heaviside function as: HSh i=H(ξ)−H(ξi), i= 1, nnodes (4.23) Chapter 4. The Extended Finite Element Method 56 Table 4.3: A possible basis for the incompressible deformation space of an element totally crossed by a discontinuity. Incompressible deformation modes reproduced by all the integration rules m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 (a) (b) Figure 4.14: Possible crack configurations containing crack kinks, with integration rules obtained considering subdivision into six triangles, four in one side of the crack and two in the other side. (a) (b) Figure 4.15: Examples of integration rules obtained for cracked hexahedrons: a) four points per tetrahedron; b) one point per tetrahedron. conforming the crack plane, in which a regular Gaussian rule is applied, as shown in figure 4.15. For example, for an hexahedral element with 1 ×12 Gauss points, corresponding to an element subdivision of twelve tetrahedrons, there are thirty seven linearly independent admissible incompressible deformation modes, whereas considering 4 ×12 Gauss points, for the same subdivision in twelve tetrahedrons, it is only possible to reproduce thirty four linearly independent modes. A complete basis for Ihis represented in table 4.4. Chapter 4. The Extended Finite Element Method 57 Table 4.4: A possible basis for the incompressible deformation space of an element totally crossed by a discontinuity. Incompressible deformation modes reproduced by all the integration rules m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14 m15 m16 m17 m18 m19 m20 m21 m22 m23 m24 m25 m26 m27 m28 m29 m30 m31 m32 m33 m34 Incompressible deformation modes only reproduced by a reduced integration rule m35 m36 m37 Chapter 4. The Extended Finite Element Method 58 (a) (b) (c) (d) Figure 4.16: Regular Gaussian quadrature rules applied to an element crossed by a discontinuity: a) four Gauss points; b) nine Gauss points; c) sixteen Gauss points; d) nine Gauss points. Once more, we conclude that is possible to expand the basis of incompressible deformation modes by using a reduced integration rule, which will be crucial to control volumetric locking. 4.2.2.2 Regular Gauss point distribution The simplest strategy to avoid the subdivision of the finite elements for integration purposes consists of placing a large number of Gauss points whose position is fixed in the element, and therefore does not adjust to the crack path [110–112], as illustrated in figure 4.16. Although this technique is relatively simple, the characterization of the underlying space of incompressible deformations associated with each Gauss point distribution may not be as systematic as in the case of element subdivision. For example, in the configurations given in figure 4.16 it was expected that Ihcould be characterized by the same basis for the cases in figure 4.16(b) and 4.16(d). Nevertheless the dimension of Ihis higher in the second case. In fact, in this case, the number of linearly independent incompressible modes depends not only on the number of Gauss points but also on the particular crack configuration. In figure 4.16(a), the integration rule yields the same solution as the one in which only one Gauss point per triangle is used in the strategy defined in section 4.2.2.1, and the twelve modes represented in table 4.1 are reproducible; for higher number of Gauss points like in figures 4.16(b)-4.16(c), or even twenty-five or sixty-four Gauss points, only the first ten modes are included in the basis of Ih. Nevertheless, when the crack configuration is slightly changed, as illustrated in figure 4.16(d), the number of linearly independent incompressible deformation modes increases to eleven, when only ten modes were expected. Chapter 4. The Extended Finite Element Method 59 (a) (b) (c) (d) Figure 4.17: Application of the Schwarz-Christoffel conformal mapping to build an integration rule for elements crossed by a discontinuity: a) cracked element; b) polygon that will be mapped to the unit disk; c) unit disk containing the Gauss point positions; d) final Gauss point distribution on the upper part of element a). These results suggest that when the integration rule does not conform the discontinuity, more likely spurious deformation modes may arise and the construction of strategies to avoid locking may not be so systematic as in the case of element subdivision. In addition, this strategy can become computationally more expensive than the referred subdivision due to the much larger number of Gauss points that has to be used to obtain the same accuracy. 4.2.2.3 Schwarz-Christoffel conformal mapping Another possible strategy used to avoid the subdivision of finite elements to perform numerical integration, based on the Schwarz-Christoffel conformal mapping was proposed by Natarajan et al. [18]. The Schwarz-Christoffel conformal mapping is an angle preserving transformation, which maps the complex upper half-plane, z={ζ∈C: Imζ > 0} to the interior of a polygon. Mathematically it may expressed as a function f(ζ), ζ ∈C, as follows: f(ζ) = ZζK (w−a)α/π(w−b)β/π(w−c)γ/π . . .dw(4.29) where Kis a constant, α,β,γ,. . . are the interior angles of the polygon and a,b,c,. . . are the values, along the real axis of the z-plane, of the vertices of the polygon. To apply this transformation to enriched XFEM elements, each polygon belonging to an element cut or intersected by a crack is mapped onto a unit disk (figure 4.17) through the inverse of the described function. The unit disk contains equally spaced Gauss points on the complex plane, which are mapped back into the original polygon by the transformation function [119, 132]. Chapter 4. The Extended Finite Element Method 60 (a) (b) (c) (d) (e) Figure 4.18: Integration rules built using the Swchwarz-Christoffel conformal mapping: a) four Gauss points; b) eight Gauss points; c) twelve Gauss points; d) twelve Gauss points; e) twenty-four Gauss points. . These figures were created making use of the MATLAB SC Toolbox [1]. The Schwarz-Christoffel conformal mapping technique may be a applied to an element totally crossed by a crack to construct the integration rules represented in figure 4.18. There is a certain freedom to distribute the Gauss points depending on the subdivision of the unit circle, as illustrated by figures 4.18(c) and 4.18(d), where twelve Gauss points are distributed in a different way. In line with the previous sections, equation 4.20 was applied to the integration rules of figure 4.18 to characterize the space of admissible incompressible deformations. Results exhibit the same structure of section 4.2.2.1, that is, twelve linearly independent incompressible deformation modes for four Gauss points (figure 4.18(a)), which corresponds to one Gauss point per triangle, when element sub division is used, and ten linearly independent incompressible deformation modes when more than four Gauss points are employed (figures 4.18(b)-4.18(e)). A possible basis for Ihis also represented in table 4.1. Therefore, we conclude that the incompressible deformation modes associated to integration rules built with the Schwarz-Christoffel conformal mapping are exactly the same as the ones obtained with the sub division in triangular sub cells technique, as long as the number of Gauss points above and below the crack is equal for both methods. Although this methodology suppresses the need of sub cells, it does not bring any advantage in terms of incompressible problems. Moreover, it requires solving two differential equations for each element totally crossed by a crack, which is computationally more expensive than to subdivide the element in triangular sub cells, and, as additional disadvantage, it is not readily extendable to 3D problems. Chapter 4. The Extended Finite Element Method 61 4.2.3 Finite element formulations for incompressibility The results previously obtained have shown, as expected, that it is possible to expand the space of admissible incompressible deformations, when an integration rule with a reduced number of Gauss points is used in the constraint equation which defines the incompressibility condition. Nevertheless these rules may introduce spurious deformation modes that should be controlled. This fact suggests that classical methods used to avoid locking in the framework of FEM, namely the B-bar method [123] and the F-bar method [127], should be used in conjunction with XFEM elements in order to prevent volumetric locking whilst preventing spurious deformation modes. 4.2.3.1 B-bar method The B-bar method [123] is based on the additive decomposition of the matrix of the shape functions derivatives, B, into a volumetric component, Bdil and a deviatoric component, Bdev, as follows: Bdev =B−Bdil (4.30) To avoid volumetric locking, Bdev is integrated with a complete Gaussian rule, whereas Bdil is integrated with a reduced Gaussian rule, which is able to reproduce a larger number of admissible incompressible deformations than the complete counterpart. In 2D problems, the B-bar formulation is usually applied to low order quadrilateral elements and Bdev is calculated using 4 Gauss points and Bdil is calculated only for one Gauss point in the center of the element, as shown on figure 4.19(b), while in 3D problems, the complete integration rule consists of eight Gauss points and the reduced integration rule consists of one Gauss point in the center of the hexahedral element. Here we propose to extend the B-bar methodology to elements crossed by discontinuities, based on the element subdivision procedure of section 4.2.2.1. Thus, in these finite elements the deviatoric part of the Bmatrix is calculated using a 3 points per triangle rule or four points per tetrahedron in 2D and 3D problems, respectively, while the volumetric contribution of each Gauss point is calculated in the center of corresponding triangle (figure 4.19) or tetrahedron. As the numerical integration rule consisting of one point per triangle is able to reproduce more incompressible deformation modes than the higher order rules, this methodology efficiently alleviates the volumetric locking in infinitesimal strain problems, as will be later illustrated with some numerical examples. Chapter 4. The Extended Finite Element Method 62 (a) (b) Figure 4.19: Integration rule for a regular 4-nodes element: 1 Gauss point for the reduced rule and 4 Gauss points for the complete rule; b) Integration rule for an element crossed by a discontinuity: 1 Gauss point per triangle as a reduced rule and 3 Gauss points per triangle as a complete rule. Although the extension of this methodology to XFEM is straightforward, it should be noticed that for cracked elements, the Bmatrix contains some extra terms due to the enrichment functions. 4.2.3.2 F-bar method For problems involving finite deformations, one of the simplest ways of taking advantage of the expansion of the incompressible deformation space, when a reduced integration rule is used, is through the F-bar method [63, 127]. While B-bar is based on the additive decomposition of the Bmatrix, the F-bar is based on the multiplicative decomposition of the deformation gradient[65–67], as described in chapter 2: F=FdevFvol (4.31) where Fdev = (det F)−1 3F=J−1 3F(4.32) Fvol = (det F)1 3I=J1 3I(4.33) The F-bar deformation gradient is simply obtained by replacing the regular volumetric part by the volumetric part, F0, calculated at the centroid if the element (for a 4-node quadrilateral or a 8-node hexahedron), reading: ¯ F=Fdev(F0)vol = (det F0 det F)1 3F= (J0 J)1 3F(4.34) The application of the F bar method to enriched elements follows the same principle as the B-bar method, hence, constructing the numerical integration rules in a sub triangulation base, for each Gauss point of the complete rule (typically 3 Gauss points per triangle/4 Gauss points per tetrahedron), the volumetric part of the deformation Chapter 4. The Extended Finite Element Method 63 gradient, F0, will be calculated at the centroid of the respective triangle/tetrahedron (figure 4.19). As in the B-bar formulation, it should be noticed that the deformation gradient of the elements crossed by discontinuities contains some extra terms associated with the enrichment of the displacements, which is clear when Fis defined in terms of the displacement field: F= (I+∇u) (4.35) Therefore, considering a general enrichment function, ψ, and respective associated degrees of freedom, ψ, the deformation gradient can be written as: F=I+ n X i=1 ui⊗∇XNi+ m X j=1 ψj⊗∇XNiψ(4.36) recalling that ∇Xdenotes the gradient with respect to the reference configuration, usually implying a coordinate transformation to the isoparametric mapping. Departing from the enriched form of the deformation gradient it is easy to adapt a finite strain code to include XFEM enrichment functions. Readers may consult references [133–135] for details. The presented extension of the B-bar and F-bar methodologies to XFEM is expected to alleviate locking in fracture problems where the satisfaction of the incompressibility constraint plays a major role. Therefore, the performance of these techniques will be illustrated in the next section with some numerical examples involving linear and nonlinear materials. 4.2.4 Numerical examples In this section, the conclusions regarding the numerical integration strategies and their ability to reproduce incompressible states will be illustrated through some numerical examples, not only focusing on each type of enriched element by itself, but also analysing their influence when inserted in a mesh. In the examples incompressibility in a XFEM framework is addressed in both small and finite strain problems. Chapter 4. The Extended Finite Element Method 64 (a) (b) (c) Figure 4.20: Test examples for single element containing a crack. 4.2.4.1 Single Element Analysing table 4.2 from section 4.2.2.1, it is prominent the absence of the deformation modes m6 and m7 for a complete integration rule, while the mode m4 may be reproduced by both complete and reduced integration rules. This fact suggested the single element test illustrated in figure 4.20, where a 4-node element in a plane strain condition is submitted to 3 different loading scenarios. The material considered is linear elastic and the nearly incompressible behavior is obtained setting the value of the Poisson coefficient to ν= 0.499999. The Young0s Modulus has a value of E= 103and the load is unitary. All the quantities are supposed to be expressed in a consistent unit system. The element contains a crack segment starting at one quarter of an edge and finishing at the element center. The element edges have length 2. In this analysis, the complete integration rule consists of eighteen integration points placed on six sub-triangles and the reduced integration rule consists of a single integration point per triangle. The B-bar technique is applied for the combination of these two rules as previously described. The tip enrichment is the one as defined in equation 4.8. In the table 4.5 it is possible to observe the deformed configuration obtained as a function of the loading and integration rule used. As it was expected, using a complete integration rule, the element exhibits volumetric locking for the cases b) and c) in figure 4.20. Differently, both reduced integration and B-bar methodologies are able to reproduce the correct deformations and do not present any locking behaviour. Nevertheless, the reduced integration may encompass some drawbacks when used on a finite element mesh, which will become clearer in the following numerical examples. Chapter 4. The Extended Finite Element Method 65 Table 4.5: Deformed configurations for different integration rules. Integration Rule Problem a b c Complete Reduced B-bar 9 9 Figure 4.21: Cracked plate with respective boundary conditions. 4.2.4.2 Incompressible behavior of XFEM elements in a mesh The problems presented in this section intend to analyse the influence of enriched elements within a mesh under an incompressibility constraint. A similar example as the previous one is tested with two types of mesh, containing regular and enriched elements as shown in figure 4.21. The material properties are the same as in the previous example. The plate is discretized with a coarse mesh and with a fine mesh containing regular and enriched finite elements, as can be observed in figures 4.21 and 4.22. Chapter 4. The Extended Finite Element Method 72 Figure 4.30: Comparison of the vertical displacement of the top right corner node, as a function of the number of elements per side, for different tip enrichment functions. Finite Strain Analysis In this section the same example but at finite strains is used to evaluate the performance of the proposed extension of F-bar formulation to XFEM and compare it with the enhanced strain formulation for XFEM proposed by Dolbow et al. [129]. A hyperelastic material response governed by the following stored energy function, W, is considered: W=µ 2[J−2 3tr[C]−3] + κ 2(J−1)2(4.37) where µrepresents the shear modulus and κthe bulk modulus. The nearly incompressible behavior is obtained setting the material parameters as κ= 40.0104Pa and µ= 80.0 Pa. In this type of constitutive behaviour, the relevant physical quantities can be derived directly from the energy function making it appealing for automatic differentiation and automatic code generation [84, 85]. Therefore, the previously described AceGen and AceFEM systems will be employed. In particular, for this hyperelastic problem, we will consider the volumetric/deviatoric split of the stored energy function [62, 138, 139], as follows: W=Wdev +Wvol (4.38) Wdev =µ 2[J−2 3tr[C]−3] (4.39) Wvol =κ 2(J−1)2(4.40) As outlined in chapter 3, using a Newton-Raphson procedure, the internal force, Re, may be obtained integrating, over the domain, the derivative of the strain energy in Chapter 4. The Extended Finite Element Method 73 respect to the displacement field. Consequently, the tangent stiffness matrix, KTis obtained as the derivative, in respect to the displacement field, of the internal force vector: Re=ZΩ ∂W ∂udΩ (4.41) KT=∂Re ∂u(4.42) Here, the derivative of the strain energy potential is presented in terms of the deviatoric part of second Piola-Kirchhoff tensor, S, and the hydrostatic pressure p: ∂W ∂u=1 2Sdev :Cdev ∂u−p∂J ∂u(4.43) with Sdev = (2∂Wdev ∂C)dev = 2∂Wdev ∂Cdev −2 3tr[Cdev ∂Wdev ∂Cdev ]C−1 dev (4.44) and p=−∂Wvol ∂J (4.45) where Cdev corresponds to the deviatoric part of the right Cauchy-Green tensor, depending on the deviatoric part of the deformation gradient, ¯ F: Cdev =¯ FT dev ¯ Fdev =J−2 3C(4.46) With the AceGen system all the presented derivatives are calculated automatically and optimized. It should be noticed that the global vector of unknowns includes the displacements as well as the extra degrees of freedom introduced by the XFEM. The analysis is performed using AceFem. For the analysis of the finite strain version of the Cook0s Membrane, a load of 1 Nis applied according to figure 4.27 during 30 load steps. Only the Heaviside enrichment is employed to model the crack. In line with the infinitesimal strain version, the vertical displacement of the top right corner is plotted as a function of the number of finite elements per edge. In addition the incompressibility constraint, J= det F= 1, is checked for all the Gauss points in the mesh. Initially, in order to evaluate the model implementation in the AceGen system, the solution for the example with the explicitly meshed crack, illustrated in figure 4.28(b), is compared with the one obtained in reference [129]. As can be seen in figure 4.31 a Chapter 4. The Extended Finite Element Method 74 Figure 4.31: Vertical displacement of the top right corner node, as a function of the number of elements per side, obtained for the explicitly meshed crack of figure 4.28(b), using the F-bar methodology and the Enhanced strain approach used by Dolbow et al. Figure 4.32: Comparison of the vertical displacement of the top right corner node, as a function of the number of elements per side, for explicitly meshed and enriched approximations using the F-bar methodology and the Enhanced Strain approach employed by Dolbow et al. very good correlation is obtained. Next the F-bar methodology applied to the XFEM is assessed. In figure 4.32 its results are compared with those obtained by Dolbow et al. [129] using the enhanced strain method and having as a reference solution the one with the crack explicitly meshed. Although, as expected, it converges slightly slower than the enhanced strain formulation developed by Dolbow et al., it reproduces better the incompressibility condition as can be observed in table 4.6. The F-bar methodology for XFEM approximates the incompressibility condition better than the Assumed Enhanced Strain approach developed by Chapter 4. The Extended Finite Element Method 75 Table 4.6: Maximum value of the deviation from the incompressibility condition, max |J−1|. Explicitly Meshed Enriched Elements per side F-bar Enhanced Dolbow et al F-bar for XFEM Enhanced Dolbow et al 5 0.000030 0.0004 0.00048 0.018 10 0.000047 0.0007 0.00047 0.024 20 0.000075 0.0012 0.00046 0.025 30 0.000093 0.0016 0.00044 0.030 -0.214e-5 0.1549e-4 J-1 (a) (b) Figure 4.33: a) Maximum value of the deviation from the incompressibility condition, (J−1), in the element crossed by the crack marked in the mesh b). Dolbow et al. Nevertheless, the value of Jslightly depends on the mesh, namely, on the size of the triangles obtained during the sub cell division for integration. In general, the higher deviation from J= 1 occurs around the crack tip, in the smaller sub cells, as depicted in figure 4.33. Besides the good performance exhibited, the F-bar method is relatively simple to implement on computer programs, making the proposed approach for XFEM suitable for simulating crack propagation in problems involving large strains where incompressibility plays a major rule, such as in elasto-plastic materials. 4.3 Interaction with non-linear material models In the previous section, suitable strategies for numerical integration in XFEM elements in the presence of incompressibility constraints were developed. Nevertheless, in a ductile fracture problem the crack path is not defined a priori but the discontinuities are Chapter 4. The Extended Finite Element Method 76 ξ η ξi=1 √3 ξi=−1 √3 ηi=−1 √3 ηi=1 √3 Ni=1 4(1 + 3ξiξ)(1 + 3ηiη) Figure 4.34: Definition of interpolatory shape functions used in the variable transfer. rather introduced at some stage of material degradation. This means that the positions and weights of the Gauss points of a certain finite element may change during the problem analysis. As the Lemaitre model for damage presented in chapter 2 is history dependent, a proper variable transfer strategy has to be defined. Here a local smoothing technique, consisting of a bilinear extrapolation of the quantities from the existing Gauss points to the new, using shape functions, is employed. This technique is used in reference [80] to transfer quantities stored in the Gauss points to the nodes. Basically, any variable in the new Gauss points, αnew, may be obtained using the following expression: αnew = ngpt X i Ng iαi(4.47) where αirepresents the values of the variable αat each one of the ngpt old Gauss points. The interpolatory shape functions are defined at the set of old Gauss points according to figure 4.34. In the particular case of the damage model considered, the variables which reflect the history of the deformation are the plastic multiplier, the strain and the damage. The stresses may be recovered using the return-mapping equations. It should be noticed that this technique differs from remeshing, as only the quantities stored in Gauss points are transferred inside the same finite element; the mesh, as well as the quantities stored in the nodes, such as the displacements, remain unchanged. Remark 4.2. Depending on the particular problem, different strategies to recover equilibrium, after the variable transfer may be required. These strategies are developed and discussed in chapter 5. Chapter 4. The Extended Finite Element Method 77 lc lc Figure 4.35: Group of Gauss points (in black) that influence a particular one. When a Gauss point is located near a discontinuity, only the points in the same side of the discontinuity is selected. P1 P2 P3 Figure 4.36: Group of Gauss points (in black) that influence a particular one, located near a crack tip. In chapter 2, a non-local formulation to prevent mesh pathologies, when implementing the damage model was also presented. In this methodology, the characteristic length, lc, determines which points should influence the history of a particular point x(figure 4.35). Nevertheless, when a crack is inserted in the mesh, its surfaces no longer interact, especially in the case of a traction-free crack. Therefore, the Gauss points near the interface should only be influenced by other Gauss points laying in the same side of the crack, as illustrated in figure 4.35. If a point is located close to the crack tip, different selection criteria may be adopted. Considering only a normal level set function, all the points P1, P2 and P3, shown in figure 4.36, are excluded from the group of Gauss points influencing the red point. In the other hand, considering an additional tangent level set function, P1 may be included in the set. Despite the different options, the influence in the results is nearly unnoticeable. Chapter 4. The Extended Finite Element Method 78 4.4 Conclusions In this chapter the general formulation of the XFEM was presented, focusing in the enrichment functions suitable to describe a ductile fracture process. The basic implementation steps were outlined, followed by a comprehensive study of integration strategies, which is one most relevant parts of this chapter. While in most of the existing studies on integration rules for XFEM, the accuracy of the integration schemes is based on the ability to reproduce the stress intensity factors in problems featuring known analytical solutions, here we favoured an approach based on the analysis of the underlying sub-space of incompressible modes embedded in the XFEM approximation. Several integration rules were compared in terms of the deformation modes that can actually be captured. In general, considering the sub division of the elements crossed by a crack it is possible to expand the underlying sub-space of incompressible deformations using a reduced integration rule over each sub cell, when compared with full integration schemes. This fact motivated the extensions of B-bar and F-bar methodologies for XFEM, which somehow rely on the selective integration rules, with encouraging results. In the final part of the chapter some more features required for the combination of the XFEM with non-linear materials, namely, the variable transfer and the non-local interaction, were also addressed. All these techniques are useful in the study of the final stage of failure of materials which undergo large strain processes without noticeable volume changes, as it will be described in the following chapters. Chapter 5 Ductile Fracture Model: Traction-Free Cracks Material failure is often related to the development of micro-defects at the microstructural level, in particular to damage mechanisms. Damage is characterized by the nucleation growth and coalescence of voids, which results in a loss of mechanical properties and leads to fracture [27]. In the case of ductile metals, the so-called ductile damage is strongly related with large plastic straining. The nucleation of voids occurs after a certain threshold of plastic strain and their growth and coalescence are governed by plastic instability, as schematically illustrated in figure 5.1 [27, 69]. This type of metals is usually constituted by second phase particles or inclusions immersed in a matrix (figure 5.1(a)). When sufficient stress is applied, matrix-particle debonding or particle cracking may occur, leading to the nucleation of voids (figure 5.1(b)). Subsequently these voids may grow (figure 5.1(c)) and, after reaching a certain size, they start interacting with the neighbouring voids. Plastic strain tends to concentrate along a preferential direction, causing local necking instabilities (figure 5.1(d)). Finally, voids coalesce and a fracture surface is created (figure 5.1(e)). The Lemaitre model for ductile damage presented in chapter 2 is able to describe material degradation by the introduction of a damage variable, D, which accounts for the loss of stiffness due to the presence of voids in an averaging way. Nevertheless, this material model is unable to capture the last failure stage corresponding to crack initiation and propagation. 79 Chapter 5. Ductile Fracture Model: Traction-free Cracks 80 (a) (b) (c) (d) (e) Figure 5.1: Ductile fracture process: a) inclusions in a metallic matrix; b) void nucleation; c) void growth; d) strain localization and necking between voids; e) void coalescence and formation of a fracture surface. A complete ductile failure model may be constructed combining the damage material model with the XFEM discretisation, introduced in chapter 3. The stage illustrated in figure 5.1(d) may be associated to a critical damage value, Dc, which is determined through the continuous theory. When this damage level is reached, a discontinuity surface may be inserted in the corresponding region using the XFEM. Therefore, in the first section of this chapter a damage based criterion for transition from damage to fracture is proposed. The same ideas are applied to crack propagation and therefore the two phenomena are dealt with in an unified way. In terms of numerical implementation the continuous-discontinuous model requires some further developments discussed in section 5.2. Next, in sections 5.3 and 5.4, the efficiency of the proposed model is illustrated through various numerical examples. This chapter finalizes with some general conclusions regarding the full ductile fracture model. 5.1 Transition criterion from damage to fracture The transition criterion from Damage to Fracture may be simply stated as a crack is inserted when critical damage is reached in the corresponding region of the domain. The Chapter 5. Ductile Fracture Model: Traction-free Cracks 81 main question is how to determine the crack characteristics, namely the initiation point, direction and length, from the critically damaged area. In a problem discretised through the FEM, the values of the damage variable are stored in each Gauss point, from which the damage distribution pattern follows directly. This information may be used to determine the characteristics of a crack, but to meet a good accuracy and flexibility an interpolation strategy to calculate the damage value at any point of the domain must be defined. A simple way to obtain the damage value at an arbitrary point is to use bilinear extrapolation with Lagrange polynomials. The damage value is determined making use of the values at the Gauss points of the element in which the point is contained, similarly to the variable transfer technique described in chapter 4. Nevertheless, to increase the accuracy of the crack definition it is interesting to gather the information of a patch of elements and, alternatively, the approximation based on B-spline functions [140–142] may be employed. The B-spline basis functions are defined recursively from a set of control points, [ξ1, ξ2, . . . , ξn, using the following recursion formula [143, 144] Ni,0(ξ) = (1 if ξi≤ξ≤ξi+1 0 otherwise (5.1) Ni,p(ξ) = ξ−ξi ξi+p−ξi Ni,p−1(ξ) + ξi+p+1 −ξ ξi+p+1 −ξi+1 Ni+1,p−1(ξ) (5.2) A B-spline curve, C, in Rdis built as a linear combination of B-spline basis functions as: C(ξ) = n X i=1 Ni,pPi(5.3) where idenotes the i-th control point and not each one of its coordinates. Unlikely Lagrange polynomials, functions constructed in this way are monotone and consequently do not introduce new maxima in the distribution. Moreover, a B-spline object of dimension dis itself a B-spline of dimension d−1, which means that a variable at a point at the boundary of a B-spline surface/volume may be interpolated with the same functions as a point lying inside the domain. These properties may bring important advantages in the treatment of a damage distribution, insuring that damage does not grow artificially due to the interpolation technique and allowing an unified treatment of crack initiation Chapter 5. Ductile Fracture Model: Traction-free Cracks 88 following section. 5.2.2 Equilibrium recovery algorithm Regarding numerical implementation, the equilibrium recovery required when new elements are enriched is the most critical part of the model. Each time new elements are enriched due to crack initiation and/or propagation, the number of degrees of freedom increases, as well as the number of Gauss points inside the enriched elements and hence the equilibrium between internal and external forces has to be re-established. In fact, when an element is enriched the positions and weights of the Gauss points change and consequently the internal variables have to be transferred from the old Gauss points to the new ones, which is done by applying a local smoothing technique as explained in chapter 4. It should be noticed that this process differs from remeshing as only the quantities stored in the Gauss points, namely the plastic strain and the damage variable, are transferred in a process occurring only inside each element; the mesh, as well as the quantities stored in the nodes remain the same. In a second stage, a load increment with null amplitude is performed. However, due to the highly non-linear nature of the problem, this procedure may be inefficient and a more complex algorithm had to be developed, as described in the following paragraphs. The key feature for equilibrium recovery is the constancy of the nodal quantities. When a discontinuity is introduced in a particular element, the continuous problem should be equivalent to the continuous-discontinuous problem. Even though the approximation function for the element changed from: u(x) = n X i=1 Niui(5.10) to u(x) = n X i=1 Niui+ nsplit X j=1 Nj[H(x)−H(xj)]aj(5.11) in elements which became totally cut by a crack, or to u(x) = n X i=1 Niui+ ntip X j=1 NjR(x)[H(x)−H(xj)]bj(5.12) Chapter 5. Ductile Fracture Model: Traction-free Cracks 89 in elements, which became tip elements, after solving the system of equations, the nodal displacements should remain the same. Therefore, when an element is about to be enriched, its displacements, unenrich, are stored. Then, after introducing the discontinuity, the problem is solved imposing the stored displacements, unenrich, as prescribed displacements, resulting in a set of residual forces, FRes, at the enriched nodes. Subsequently, the nodes are released and the problem is solved for progressively reduced residual forces, according to the scheme in figure 5.6. To finalize the algorithm description, it should be noticed that, although the crack length step is calculated using the methodology proposed in section 5.1.2, sometimes it is divided in several steps to facilitate equilibrium recovery. Moreover, during the analysis time stepping is usually not constant. When critical damage is reached, the program goes back one time step and continues with a fraction of the first time step according to the scheme in figure 5.7. In this way the crack initiation stage is captured with more detail. As crack propagation in ductile metals is usually stable it may be successfully simulated using a quasi-static formulation. Now that the fracture model is fully presented, its effectiveness will be illustrated through various numerical examples, presented in the remaining of this chapter. Chapter 5. Ductile Fracture Model: Traction-free Cracks 90 Crack growth step Select nodes to enrich Store the displacements of newly enriched nodes, unenrich Perform Newton iteration with null multiplier convergence Yes Continue analysis No set unenrich as prescribed displacements and solve the system of equations system of equations convergence Yes No Reduce crack growth step Resultant set of residual forces, FRes Release the imposed displacements on the enriched nodes Reduce progressively the residual forces until convergence is obtained for FRes < tol by a certain percentage α FRes = (1 −α)FRes Newton iteration convergenceYes No FRes < tol α≤αmin Yes No No Yes α=α+αReduce α Continue analysis Figure 5.6: Iterative scheme for equilibrium recovery during a crack growth step. Chapter 5. Ductile Fracture Model: Traction-free Cracks 91 Figure 5.7: Time-stepping scheme before and after crack initiation/propagation. 5.3 Numerical examples 1: Plate with an initial crack Before evaluating the total fracture model, involving crack initiation and propagation, a specimen with an initial crack is studied. This numerical example is intended to evaluate the effects of the non-local damage formulation, when applied in conjunction with the XFEM. In a second phase, damage based crack propagation is inspected as well. Following this line, this example consists of a cracked plate illustrated in figure 6.9. The material properties and dimensions are summarized in table 5.1. Three different FEM meshes are considered, as shown in figure 5.9, under plane strain assumption. For the analysis, besides the material parameters defined in table 5.1, the proposed model requires a value for the non-local regularization length, lr. This parameter cannot be measured directly from experiments but is usually obtained through inverse analysis [77, 79]. As the presented non-local formulation has proven to be able to alleviate pathological mesh dependence, lrmay be regarded as a numerical regularization parameter to be prescribed by the user. Therefore, three different situations, corresponding to three values of lrare analysed: lr= 0 (the local case), lr= 0.8mm and lr= 1.6mm. The damage contours obtained in each case, for an applied displacement of 0.05mm, are Chapter 5. Ductile Fracture Model: Traction-free Cracks 92 a b 2h Figure 5.8: Pre-cracked plate Table 5.1: Material properties of the plate with an initial crack. Property Value Elastic modulus E= 206.9 GPa Poisson’s ratio ν= 0.29 Damage exponent s= 1.0 Damage denominator r= 1.25 MPa Hardening function τy(R) = 450 + 129.24R+ 265(1 −e−16.93R) MPa Dimension h= 15 mm Dimension b= 20 mm Initial crack length a= 20/3 mm displayed in figures 5.10 - 5.12. Results demonstrate a clear mesh dependence pathology in the local case. Instead of providing better resolution in the crack tip vicinity, mesh refinement causes spurious damage localization. As expected, the non-local formulation alleviates the mesh dependence, as can be observed in figures 5.11 and 5.12, where the damage contours do not vary significantly in the different meshes. Nevertheless, for lr= 0.8mm, there is still some excessive damage concentration at the crack tip, which suggests that the larger lris a better choice. In fact, when analysing the crack growth, lr= 0.8mm will prove to be insufficient. As a final note on the use of the Lemaitre model under a non-local formulation with discontinuous cracks, it should be noticed that the maximum damage values region is connected Chapter 5. Ductile Fracture Model: Traction-free Cracks 93 (a) (b) (c) Figure 5.9: Mesh refinement for the plate with an initial crack: a) 20 ×31, b) 40 × 61 and c) 50 ×75 elements. (a) (b) (c) Figure 5.10: Damage contours obtained for the local case, for an applied displacement of 0.05mm for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75 elements. (a) (b) (c) Figure 5.11: Damage contours obtained for lr= 0.8mm, for an applied displacement of 0.05mm, for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75 elements. Chapter 5. Ductile Fracture Model: Traction-free Cracks 94 (a) (b) (c) Figure 5.12: Damage contours obtained for lr= 1.6mm, for an applied displacement of 0.05mm, for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75 elements. with the crack tip, and therefore correlates with micro mechanical approach described previously in this chapter: the triaxiality state at the crack tip favours the nucleation and growth of voids, which eventually connect to the main crack causing its propagation. Ductile Crack Growth The analysis of the damage-based crack growth is now presented. One more model parameter, the critical damage value, Dc, which triggers the crack advance must be defined. Following the work of Lemaitre [27], in real parts and structures the critical damage value is not 1, which would correspond to theoretical fully damaged material, but is rather located between 0.2 and 0.5. At this stage, the values of 0.2 is employed and further discussion is left for the full model involving crack initiation and propagation. In this study, the same three meshes and non-local lengths were employed. Results are depicted in figures 5.13 and 5.14, where the final crack paths and damage contours may be observed. For the regularization length of 0.8mm there is still some damage localization at the crack tip. In the coarse mesh (figure 5.13(a)) just a few Gauss points in the vicinity of the crack tip reach the critical damage. Consequently the convergence of the iterative scheme is affected during the crack growth simulation. The finer meshes (figures 5.13(b) and 5.13(c)) provide more resolution, however only with the later is possible to simulate the expected crack path. This case illustrates one of the disadvantages of the proposed Chapter 5. Ductile Fracture Model: Traction-free Cracks 95 (a) (b) (c) Figure 5.13: Crack path obtained using lr= 0.8 mm, for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75 elements. Chapter 5. Ductile Fracture Model: Traction-free Cracks 96 (a) (b) (c) Figure 5.14: Crack path obtained using lr= 1.6 mm, for thee meshes a) 20 ×31, b) 40 ×61 and c) 50 ×75 elements. Chapter 5. Ductile Fracture Model: Traction-free Cracks 97 Figure 5.15: Reaction force in function of the applied displacement lr= 1.6mm. damage criterion: due to its geometric nature, an accurate calculation of the crack direction, requires that critical damage is reached in a certain number of Gauss points. Nevertheless, this disadvantage is usually overcome by selecting an adequate non-local length, as it will be demonstrated in the remaining examples. In particular, for the length of 1.6mm, crack growth can be simulated successfully for all the meshes. Nevertheless, when the specimen is close to complete rupture, some numerical instabilities may arise causing a slight change of the crack direction. In figure 5.15 the reaction force-applied displacement curves are displayed. The curves are convergent upon mesh refinement, which demonstrates the efficiency of the proposed crack growth criterion, based on damage evolution. To finalize this example, the crack length evolution is depicted in figure 5.16. The Chapter 5. Ductile Fracture Model: Traction-free Cracks 104 Figure 5.25: Relation between the residual force and the displacement applied to the top edge of the notched specimen for the mesh with 31 nodes per side, considering lr= 2.0 mm and different values for DC. calibration of the proposed methodology is dependent on the characterization of the continuum model parameters, regardless of the good performance exhibited so far. 5.4.2 Plane strain specimen In this example, the performance of the proposed methodology under the plane strain condition is tested using the specimen depicted in figure 5.26. Figure 5.26: Plane strain specimen. For this analysis, the internal length, lrhas the value of 1.6 mm and the critical damage value is 0.5. The analysis is performed for the four meshes represented in figure 5.27, prescribing adequate symmetry conditions and applying a displacement to the top and Chapter 5. Ductile Fracture Model: Traction-free Cracks 105 (a) (b) (c) (d) Figure 5.27: Mesh refinement for the plane strain specimen, mesh density of a) 11, b) 21, c) 31, d) 41 elements per side. bottom edges. In figure 5.28 the reaction force-displacement curves obtained are represented. As in the previous example, all the stages of material behaviour, including hardening and softening up to failure, are clearly represented. Results are also convergent upon mesh refinement. Under a plane strain condition, the non-local integral model also avoids pathological mesh dependence and spurious damage localization, as may be observed in the damage distribution contours illustrated in figure 5.29. Finally, in figure 5.30 different crack growth steps for the mesh with 41 elements per side are illustrated, illustrating, once more, the close relation between damage evolution and crack growth in ductile metals. Chapter 5. Ductile Fracture Model: Traction-free Cracks 106 Figure 5.28: Reaction force as a function of the applied displacement of the top nodes of the plane strain specimen. (a) (b) (c) Figure 5.29: Damage contour and final crack for mesh a) 11, b) and 21 c) 31 elements per side. 5.4.3 Double notched specimen The main objective of this example is to compare the crack path obtained with different meshes. The previous examples are not fully illustrative, as the crack always progresses in a straight line. For that purpose, a double notched specimen, as represented in figure 5.31, is chosen and is loaded so that a shear-like failure mode would occur. The analysis is performed for the two meshes represented in figure 5.32, considering a critical damage value of 0.5, a regularization length of 1.6 mm and applying a displacement of 1.5 mm. Chapter 5. Ductile Fracture Model: Traction-free Cracks 107 (a) (b) (c) (d) Figure 5.30: Different crack growth steps for the mesh with 41 elements per side. rc= 1.0 mm r1= 2.0 mm r2= 2.5 mm a= 10 mm Figure 5.31: Shear Specimen Chapter 5. Ductile Fracture Model: Traction-free Cracks 108 (a) (b) Figure 5.32: FEM meshes used in the double notched specimen: a) mesh with 16 nodes per side, b) mesh with 23 nodes per side. (a) (b) Figure 5.33: Crack path obtained for a) mesh a, b) mesh b. In figure 5.33 the crack paths obtained for the two different meshes are represented. It can be observed that the paths are nearly the same for both cases, that is, the crack path can be considered as mesh independent. The non-local integral model avoids spurious localization of damage, resulting in similar damaged areas for different mesh refinements (figure 5.34). The subsequent insertion of the crack through the XFEM respects these damage contours independently from the mesh as well, indicating that this methodology is adequate to predict failure in specimens of arbitrary shape. Chapter 5. Ductile Fracture Model: Traction-free Cracks 109 (a) (b) Figure 5.34: Damage contour and final crack for mesh a) and mesh b). 5.5 Conclusions Phenomenologically, the initiation of a crack in ductile metals is connected to the evolution of damage, which may be described by a continuum model. Nevertheless, for a complete description of the failure process, a discontinuity should be inserted in the domain, once a critical damage value is met. In this chapter, a model able to represent all the stages of material behaviour was presented. It was shown that ductile fracture may be successfully dealt with by combining the XFEM with a plastic-damageable material model. Furthermore, the proposed methodology features the following advantages: . Crack characteristics are determined directly from the continuous model and therefore there is no need of identifying additional material parameters for fracture. Moreover, crack initiation and propagation are dealt with in an unified way. . Crack initiation locus does not need to be known in advance and crack progression is independent of the mesh, saving computational resources when compared with competing approaches such as remeshing. . Results are mesh independent upon a certain mesh refinement. . The Lemaitre model in its non-local integral formulation ensures similar damage patterns in coarse and fine meshes. Therefore crack advance velocity and crack patterns are also accurately determined, even in coarse meshes. Chapter 5. Ductile Fracture Model: Traction-free Cracks 110 The possible drawbacks of the proposed methodology are related to parameter identification, such as the non-local intrinsic length and the critical damage value, which may be determined by inverse analysis using a proper combined numerical-experimental approach. From the theoretical point of view, there is an energy gap when performing a transition from damage to a traction-free crack at damage values lower than 1 and therefore the introduction of a cohesive crack may regarded as an improvement of the model and will be addressed in the following chapter. Nevertheless, we believe that the approximation developed in this chapter is good enough to describe ductile fracture in a wide range of industrial processes. Chapter 6 Ductile Fracture Model: Cohesive Cracks In the previous chapter a model for ductile fracture, governed by damage evolution, was presented. When damage reaches a critical value a crack is initiated and subsequently propagates following the damage pattern. In theory, to ensure thermodynamical consistency in the model in use, the transition from damage to fracture should occur when the material is fully degraded, that is, for Dc= 1. At this point the damage energy release would be equivalent to the energy necessary to create a crack surface [145]. However, setting the damage value to 1 leads to a singularity in the continuum equations, with its natural numerical consequences. In the previous chapter it was also outlined that experimentally the critical damage value is, in general, located between 0.2 and 0.5,and therefore there is still an energy gap which should be fulfilled. One possible solution to this problem is to add a cohesive law to the model. In the previous model it was considered that once a crack is introduced, the newly formed surfaces no longer interact, that is, the cracks were considered traction-free. A cohesive law is a traction-displacement relation which models the interaction between the two surfaces of a crack [12, 13, 22, 146, 147]. From the micro-mechanical point view, it may be interpreted as material which is not fully damaged and, consequently there is still some connection between the two crack surfaces, as schematically illustrated in figure 111 Chapter 6. Ductile Fracture Model: Cohesive Cracks 112 (a) (b) Figure 6.1: a) Traction-free crack b) Cohesive crack. 6.1. In terms of macro-modelling, the transition from damage to fracture for critical damage values lower than one could be compensated by the cohesive law. This concept has been previously applied by Cazes et all. [58, 59], in a formulation that does not require a predefined shape for a cohesive law but in which the location of the crack must be known in advance. As additional disadvantage, the model was only fully developed for 1D cases. In this work, a shape for the cohesive law will be assumed but its parameters will be fitted following energetic considerations. Therefore, this chapter initiates with a brief review of of the basic equations and implementation of a cohesive law. Then the influence of the cohesive zone in the model developed in chapter 5 is investigated. Encouraged by the results of this study, a strategy to achieve an energetically consistent transition from damage to fracture is proposed and illustrated through some numerical examples. Finally, the main conclusions of the developed approach are outlined. 6.1 Basic equations and implementation of a cohesive law 6.1.1 Variational formulation In a domain containing a cohesive crack, the stress field must be related not only to the external loading but also to the cohesive tractions, which are active in the cohesive zone. An additional condition for the cohesive interface must be satisfied together with the equilibrium equation, as described in the following paragraph. Chapter 6. Ductile Fracture Model: Cohesive Cracks 113 ∂Ωcoh Ω ∂Ωt ∂Ωu (a) t+ con t+ cot t− cot t− con n+ n− (b) Figure 6.2: a) Body containing a cohesive crack. b) Detail of the cohesive crack Considering the body containing a cohesive crack, represented in figure 6.2, the equilibrium equation for the static case is given by: divxσ= 0 in Ω (6.1) and the applied tractions, ¯ t, must satisfy ¯ t=σn at ∂Ωt(6.2) In the cohesive zone, considering normal and tangential interactions as represented in figure 6.2(b), the cohesive tractions, t− co and t+ co, are given by: t+ con =σn+=−t− con =σn−(6.3) and t+ co = (t+ con +t+ cot) = −t− co =−(t− con +t− cot) (6.4) where n−and n+are the normal vectors to the crack surfaces, as illustrated in figure 6.2(b). The cohesive tractions are considered to be a function of the crack opening, ω, defined as: ω= (u−−u+) at ∂Ωn(6.5) Equation 6.5 may be split in a normal component and a tangential component as follows: ωn=ω·n ωt=ω·t(6.6) where tis a vector perpendicular to nand, consequently: tcon =tcon(ωn) tcot =tcot(ωt)(6.7) Chapter 6. Ductile Fracture Model: Cohesive Cracks 120 a b 2h (a) (b) Figure 6.8: a) Plate with initial crack b) FEM mesh. Table 6.1: Material properties of the cracked plate. Property Value Elastic modulus E= 206.9 GPa Poisson’s ratio ν= 0.29 Damage exponent s= 1.0 Damage denominator r= 1.25 MPa Hardening function ty(R) = 450 + 129.24R+ 265(1 −e−16.93R) MPa Critical damage Dc= 0.5 Non-local length lr= 1.6 mm that the effects of the cohesive zone are noticeable in the final response of the component. The influence of the cohesive law in the reaction force-applied displacement curve of the plate, prior to crack growth is represented in figures 6.9 and 6.10. The model exhibits a higher sensitivity to Fnthan to ω0. As expected, increasing the value of Fnthe strain under the reaction force-applied displacement curve increases significantly. On the other hand, it terms of ω0, the effects of the cohesive law nearly vanish for ω0>1 mm−1, as from this value the slope of the cohesive law becomes appreciable and, consequently the crack becomes nearly traction-free. This preliminary study is finalized inspecting the crack growth phase. The effect of the cohesive law in the reaction force-applied displacement curve and in the crack length evolution are displayed in figures 6.11 and 6.12, respectively. It can be observed that the cohesive law delays slightly the propagation of the crack and increases the global strain energy, which suggests that the cohesive law could have the same effect as increasing Chapter 6. Ductile Fracture Model: Cohesive Cracks 121 Figure 6.9: Influence of Fnin the response of the cracked plate for constant ω0= 0.5 mm−1. (a) (b) Figure 6.10: Influence of ω0in the response of the cracked plate for constant a) Fn= 200 MPa b) Fn= 100 MPa. Chapter 6. Ductile Fracture Model: Cohesive Cracks 122 Figure 6.11: Influence of a cohesive law in the reaction force during crack propagation. Figure 6.12: Influence of a cohesive law in the evolution of the crack length. the critical damage value. The main objective of this problem was to check the sensibility of the damage model to the cohesive law before moving to the transition problems. However, the importance given to the magnitude of the cohesive parameters is relative, as in this work they are regarded as numerical parameters, rather than material characteristics. Chapter 6. Ductile Fracture Model: Cohesive Cracks 123 (a) (b) Figure 6.13: a) Plane strain specimen b) FEM mesh. 6.2.2 Transition from damage to fracture In the previous example it was possible to get a first estimate of the magnitude of the cohesive law parameters: ω0≤1 mm−1and Fn≈(1 −Dc)τy. Now an example considering the full model of transition from damage to fracture is inspected. It consists of a plane strain specimen, illustrated in figure 6.13, which was previously analysed in the context of traction-free cracks. In the same figure it is represented the finite element mesh employed, which has a 2205 elements. The material properties are the same used in the previous example, with the exception of the critical damage value and the parameters of the cohesive law. In the first analysis the critical damage has the value of Dc= 0.4. The reaction forceapplied displacement curve of a traction-free crack is compared with the one of a cohesive crack. The cohesive law parameters are Fn= 200 MPa and ω0= 0.5 mm−1. Results are displayed in figure 6.14. The cohesive law ensures a smoother transition from damage to fracture and a smoother crack propagation, as a result the strain energy under the reaction force-applied displacement curve increases.Comparing the results with the case of a traction-free crack where Dc= 0.5, it can be observed in figure 6.15 that the curve corresponding to Dc= 0.5 may Chapter 6. Ductile Fracture Model: Cohesive Cracks 124 Figure 6.14: Reaction force-applied displacement curves. Figure 6.15: Reaction force-applied displacement curves and respective influence of the cohesive law. be approached by the curve Dc= 0.4 enhanced by the cohesive law. This fact suggests that the curve corresponding to the traction free crack Dc= 1 could be approached by setting the cohesive parameters at the right values. To finalize this example, the crack evolution is illustrated in figure 6.16, however, the effect of the cohesive law is not so significant. The crack growth speed is still higher for Dc= 0.4 combined with the cohesive law than for Dc= 0.5. Depending on the particular application, it may be more important to have an accurate prediction of the load level/applied displacement at which the failure of the component happens or an Chapter 6. Ductile Fracture Model: Cohesive Cracks 125 Figure 6.16: Influence of the cohesive law in the evolution of the crack length. accurate prediction of the crack length evolution. The important conclusion is that simulations with high critical damages may be approached by simulations with lower critical damages, combined with a cohesive law. Therefore, in the next section an energy relation between the traction-free model and the cohesive model will be proposed, with the objective of fitting the cohesive law parameters. By setting these parameters to the right values, it will be possible to approach the transition from damage to fracture at Dc= 1. 6.3 Energy balance and calibration of the parameters of the cohesive law The early crack propagation models, based in LEFM concepts, assume that there is a certain amount of dissipated energy, which is material-specific and is responsible for a certain crack extension. Fracture is triggered by a single parameter such as the Griffith fracture energy or the stress intensity factor. Even extensions of this theories to include crack tip plasticity still rely on a single parameter and consider that the total amount of dissipated energy in a deformation process is used to propagate cracks [147, 152]. In materials which exhibit substantial deformations prior to crack growth at least two energy consuming processes should be considered, one related to plastic deformation and another related to progressive material degradation, that is, related to damage. In this context, Mazars and Pijaudier-Cabot [145] developed the equivalent crack concept, Chapter 6. Ductile Fracture Model: Cohesive Cracks 126 in which the damage energy release rate is equivalent to the crack surface energy, as follows: ZV−Y˙ DdV =−GF˙ A(6.19) In equation 6.19 Vis the overall volume of the structure, Yis the damage energy release rate (as defined in chapter 2), Dis the damage variable, GFis the fracture energy and Ais the area of the crack. One of the main limitations of this approach is that requires the value of the fracture energy, which is a parameter derived from the LEFM and its applicability in a ductile fracture context may be questionable. Alternatively, Cazes et al. [58, 59], propose a formally identical expression, where GFis replaced by the area under the traction curve of a cohesive model. The equivalence of damage to fracture is also performed locally and the model is only fully developed for 1D problems. The work developed in this chapter is conceptually related to these approaches. In fact, the results presented so far indicate that a cohesive law may actually compensate for a transition from damage to fracture before for critical damage values lower than one. Nevertheless, in opposition to the described models, we believe that the energy balance between damage and fracture should be evaluated in a global way. A local balance may be ambiguous in terms of the elements which should contribute for the energy balance: all the elements which reached critical damage in a certain region, or only the elements which will contain the crack? In addition, the direct energy transfer from the damaged volume to the crack surface admits that the micro-mechanical damage processes in that volume stop evolving and that all the energy stored in the micro-structure is transmitted to the dominant macro-crack, which is not necessarily true. Here, having in mind that the strain energy may be given by Γ = ZX(Ω) σ:ddv(6.20) a function, Γ = Γ(Dc), will be constructed for traction-free cracks, in order to approximate the value of Γ(Dc= 1). Then, departing from a certain value of Dc<1, a cohesive law will be added in order to meet the condition: Γ(Dc<1 + cohesive law) = Γ(Dc= 1) (6.21) Chapter 6. Ductile Fracture Model: Cohesive Cracks 127 In the next section, the efficiency of this methodology will be illustrated through some numerical examples. 6.4 Numerical examples 6.4.1 Plane strain specimen This section starts by going back to the plane strain specimen of section 6.2.2 and analyse the strain energy involved in each one of the cases presented: traction-free crack with Dc= 0.4, traction-free crack with Dc= 0.5 and cohesive crack with Dc= 0.4. The respective values of the strain energy are displayed in table 6.2. Table 6.2: Strain energy for Dc= 0.4−0.5. Case Strain energy - J traction-free Dc= 0.4 101.4 cohesive Dc= 0.4 121.7 traction-free Dc= 0.5 121.8 In terms of energy balance, the introduction of a cohesive law is indeed equivalent to rise the critical damage value. To emphasise this fact, the influence of the cohesive law for Dc= 0.6 and Dc= 0.7 is illustrated in figure 6.17. The respective values of the strain energy may be found in table 6.3. Choosing the values Fn= 80 MPa and ω0= 1.0 mm−1, the curve corresponding to a traction-free case with Dc= 0.7 is very well approximated. In terms of energy balance, the choice Fn= 150 MPa and ω0= 0.5 mm−1also produces good results, but the crack propagation phase is oversmoothed, suggesting that approximations using lower values of the strain energy are preferable to approximations using higher values of the strain energy. Table 6.3: Strain energy for Dc= 0.6−0.7. Case Strain energy - J traction-free Dc= 0.6 132.4 traction-free Dc= 0.7 139.6 cohesive Dc= 0.6, Fn= 80MPa, ω0= 1mm−1136.4 cohesive Dc= 0.6, Fn= 150MPa, ω0= 0.5mm−1140.1 Chapter 6. Ductile Fracture Model: Cohesive Cracks 128 Figure 6.17: Reaction force-applied displacement curves and respective influence of the cohesive law. Figure 6.18: Strain energy as a function of the critical damage. Correlation factor: R2= 0.985. To achieve the energetically consistent transition from damage to fracture, one must determine the strain energy for Dc= 1. In figures 6.18 and 6.19 the strain energy as a function of the critical damage value is represented, together with two alternative fittings: an exponential function of the type y=aebx and a polynomial function as y=ax3+bx2+cx +d. The coefficients defining each function were determined using the least-squares method. Using the exponential function the strains energy value is Γ(Dc= 1) = 188.7 J, while using the polynomial function Γ(Dc= 1) = 186.4 J. The proximity of the two values suggests that the real value of Γ(Dc= 1) should be of this order of magnitude. Therefore, departing from a critical damage value of Dc= 0.8, a cohesive law was added in order to meet Γ(Dc= 0.8 + cohesivelaw)≈187 J. Chapter 6. Ductile Fracture Model: Cohesive Cracks 129 Figure 6.19: Strain energy as a function of the critical damage. Correlation factor: R2= 0.967. Figure 6.20: Reaction force-applied displacement curves and respective influence of the cohesive law. The best fit is obtained setting the parameters of the cohesive law to Fn= 110MPa and ω0= 0.01mm−1. The strain energy obtained is Γ(Dc= 0.8 + cohesivelaw) = 185.1 J, which is quite close to the expected value. In figure 6.20 it is possible to observe the difference between the traction-free and the cohesive reaction force-applied displacement curves. Throughout this chapter it has been mentioned that a cohesive law increases the total strain energy under the reaction force-applied displacement curve, which has been verified through various numerical examples. However, in the context of this work, the cohesive law parameters are not a characteristic of the material, but numerical parameters which adapt to the critical damage value. Physically they may be interpreted as