Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Modeling of dynamic fracture problems using AL finite element formulation Abdelgalil, Abdelgader I. 2002

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2002-749819.pdf [ 14.88MB ]
Metadata
JSON: 831-1.0090649.json
JSON-LD: 831-1.0090649-ld.json
RDF/XML (Pretty): 831-1.0090649-rdf.xml
RDF/JSON: 831-1.0090649-rdf.json
Turtle: 831-1.0090649-turtle.txt
N-Triples: 831-1.0090649-rdf-ntriples.txt
Original Record: 831-1.0090649-source.json
Full Text
831-1.0090649-fulltext.txt
Citation
831-1.0090649.ris

Full Text

MODELING OF DYNAMIC FRACTURE PROBLEMS USING ALE FINITE ELEMENT FORMULATION By Abdelgader I. Abdelgalil B.Sc., M.Sc. University of Garyounis, Libya, 1988, 1994 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY In THE F A C U L T Y OF G R A D U A T E STUDIES DEPARTMENT OF M E C H A N I C A L ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH C O L U M B I A May 2002 © Abdelgader I. Abdelgalil, 2002 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mechanical Engineering The University of British Columbia 2324 Main Mall Vancouver, Canada V6T 1Z4 Date: Abstract The problem of dynamic crack propagation is widely addressed in the literature. The few available analytical solutions are limited to simple and idealized geometries and loading conditions. On the other hand, major approximations and inconsistent assumptions exist in published numerical models. In this thesis, the problem of dynamic crack propagation is modeled using a fully coupled Arbitrary Lagrangian Eulerian (ALE) formulation. The A L E equilibrium equations are derived, discretized using isoparametric finite elements and implemented into an A L E dynamic fracture program (ALEFR), based on an implicit solution scheme. The advantage of the A L E formulation is that the computational grid (finite element mesh) may have an arbitrary motion with respect to the domain of the deformed body. Therefore, the complex nature of the developed boundary condition due to a propagating crack may now be modeled in a continuous and accurate manner. The process of creating new surfaces due to crack propagation is modeled by splitting material points. This allows for a more realistic representation of the actual physical process. The A L E boundary constraint is enforced on the free boundaries, including the continuously changing free crack surfaces, using a newly developed technique. The ii dynamic energy release rate is evaluated through the integration of material properties of Lagrangian grid material points. The developed formulations and techniques are then discretized and implemented into a finite element code. The developed code is tested by modeling dynamic stationary and propagating fracture problems. iii Table of Contents Abstract ii List of Figures viii Nomenclature x Acknowledgments xvii 1 Introduction and Background 1 1.1 Motivation 1 1.2 Objective 2 1.3 Background 3 1.3.1 Numerical Modeling of Dynamic Crack Propagation 3 1.3.2 The Arbitrary Lagrangian Eulerian (ALE) Formulation 5 1.3.3 History of the A L E 7 1.4 Scope of Work 10 1.5 Organization of Thesis 11 2 Formulation and Discretization of ALE Equations 13 2.1 A L E Governing Equations 13 2.1.1 Preliminaries 13 2.1.2 Quasi-Static Analysis 16 2.1.3 Dynamic Analysis 21 iv 2.2 Finite Element Equations 23 2.2.1 Isoparametric Finite Elements 23 2.2.2 Discretization of the Quasi-static A L E Equation 25 2.2.3 Discretization of the Dynamic A L E Equation 30 3 Dynamic Fracture Parameters 33 3.1 Introduction 33 3.2 General Mechanical Energy Balance 34 3.3 Energy Balance for Crack Tip Region 35 3.4 The dynamic Energy Release Rate 37 3.5 Domain Integrals for Energy Release Rate 38 3.6 Mixed Mode Stress Intensity Factors 40 3.7 Direction of Crack Propagation 42 4 Numerical Procedure 43 4.1 M E S H MOTION 43 4.1.1 Grid Displacement 44 4.1.2 Mesh Motion for Interior Nodes 44 4.1.3 Mesh Motion on Free Material Boundaries 45 4.2 Material Point Splitting (MPS) 51 4.3 Numerical Implementations 54 4.3.1 Solution of Equilibrium Equations 54 4.3.2 Numerical Integration of the Energy Release Rate 56 v 4.4 The A L E Dynamic Fracture Program 60 4.4.1 Program Structure 60 4.4.2 Characteristics of the Developed Program (ALEFR) 64 5 Numerical Examples 66 5.1 Stationary Mode I Crack Subjected to Single Step Pulse 66 5.1.1 Problem Description 66 5.1.2 Numerical Procedure 67 5.1.3 Numerical Results 68 5.1.4 Analysis of Contour Path Dependency 73 5.2 Stationary Mode I Crack Subjected to Symmetric Step Pulse 79 5.2.1 Problem Description 79 5.2.2 Numerical Procedure 79 5.2.3 Numerical Results 79 5.3 Propagating Mode I Crack Subjected to Single Step Pulse 82 5.3.1 Problem Description 82 5.3.2 Numerical Procedure 84 5.3.3 Numerical Results 85 5.4 Propagating Mode I Crack Subjected to Symmetric Step Pulse 89 5.4.1 Problem Description 89 5.4.2 Numerical Procedure 89 5.4.3 Numerical Results 89 5.5 General Comments about Mode I Fracture Problems 90 vi 5.6 Mixed Mode Fracture Examples 92 5.6.1 Curved Crack Growth 92 5.6.2 Dynamic Mixed Mode Crack Propagation 96 5.6.3 Numerical Results 97 6 Conclusions and Future Work 100 6.1 Summary of Accomplishments 100 6.2 Future Work 101 Bibliography 103 Appendix A Calculation of the Material Rate of Cauchy Stress 109 Appendix B Treatment of Convictive Terms 110 Appendix C Discretization of the Convective Terms I l l vii List of Figures (3.1) Crack tip contours 36 (4.1) The Displacement Coupling (DC) technique 46 (4.2) The True Boundary tracking (TBT) technique 48 (4.3) Nodal motion on free boundary: TBT vs DC technique 50 (4.4) Material Point Splitting (MPS) technique 53 (4.5) Motion of the contour and domain integrals due to dynamic crack Propagation 59 (4.6) Flowchart of the developed program (ALEFR) 63 (5.1) Horizontal edge crack on a thick plate, used in mode I fracture cases .. 67 (5.2) Snap shots of ay stress distribution due to the applied tensile stress pulse, in the stationary crack model 70 (5.3) Variation of the normalized strain energy release rate with the normalized time using the fine model with (801 x 321) Lagrangian grid, for the stationary crack model subjected to tensile stress pulse 72 (5.4) Variation of the dynamic stress intensity factor with the normalized time using the fine model with (801 x 321) Lagrangian grid, for the stationary crack model subjected to tensile stress pulse 73 (5.5) Path dependency of the contour and domain integrals corresponding to the normalized energy release rate at different time steps, using (801 x 321) Lagrangian grid 75 viii (5.6) Variation of the contour and domain integrals corresponding to the normalized energy release rate with the normalized time, using (801 x 321) Lagrangian grid 78 (5.7) Snap shots of ay stress distribution due to the applied symmetric tensile stress pulse, in the stationary crack model 80 (5.8) Variation of the normalized strain energy release rate with the normalized time, for the stationary crack model subjected to tensile stress pulse 82 (5.9) Snap shots of ay stress distribution in the dynamic stationary/ propagating crack model, subjected to a tensile stress pulse 87 (5.10) Variation of the normalized strain energy release rate with the normalized time, for the dynamic stationary/propagating crack model subjected to a tensile stress pulse 88 (5.11) Snap shots of ay stress distribution in the dynamic stationary/ propagating crack model, subjected to a symmetric tensile stress pulse .. 93 (5.12) Snap shots of ay stress distribution in the dynamic stationary/ propagating crack model, subjected to a symmetric tensile stress pulse ... 94 (5.13) Snap shots of ay stress distribution in the curved crack propagation model. True crack faces are shown as red dots 95 (5.14) The double notched plate subjected to impact loading 96 (5.15) Snap shots of von Mises stress distribution during mixed mode dynamic crack propagation 98 ix Nomenclature ( p u r e mode I auxiliary field quantity ( ) ( 2 ) pure mode II auxiliary field quantity X ) quantity occurs at time t t () quantity referred to time t Ai universal function of crack tip speed for mode I deformation An universal function of crack tip speed for mode II deformation 'a element nodal referential material acceleration vector a element nodal incremental referential material acceleration vector a0 initial crack length ' ai components of referential material acceleration vector ai components of incremental referential material acceleration vector ,BAi element shape function derivative matrix related to first A L E stiffness ,BA2 element shape function derivative matrix related to second A L E stiffness ,BLI element shape function derivative matrix related to Lagrangian material stiffness tBL1 element shape function derivative matrix related to Lagrangian geometric stiffness ' C element material constitutive matrix 'CM element first convective velocity-stiffness matrix due to A L E 'CA2 element second convective velocity-stiffness matrix due to A L E ' C A 3 element third convective velocity-stiffness matrix due to A L E x ' CM element fourth convective velocity-stiffness matrix due to A L E 'Cijkl components of material constitutive tensor 'Dkl components of rate of deformation tensor c d elastic dilatational wave speed c R elastic Rayleigh surface wave speed c s elastic shear wave speed 'dA elemental domain area 'dC elemental contour length 'dS elemental surface area 'dV elemental volume E Young's modulus Ep plastic modulus , e element strain vector , etj components of deformation tensor f * element equivalent nodal force vector 'f element internal force vector ' / arbitrary function ' / material time derivative of arbitrary function ' / ' grid time derivate of arbitrary function 'ft8 components of body force per unit mass 'ff components of the externally applied surface tractions xi G normalized dynamic energy release rate G dynamic energy release rate 8 universal crack tip velocity function Gaux dynamic energy release rate of an auxiliary field G c contribution of the contour integral in expression for the normalized dynamic energy release rate Gd contribution of the domain integral in expression for the normalized dynamic energy release rate Gi„t interaction dynamic energy release rate Gtotai total dynamic energy release rate (total plus auxiliary) H element interpolation matrix ha element shape function at node a I identity matrix/tensor K normalized mode I dynamic stress intensity factor K mode I dynamic stress intensity factor k crack tip university velocity function Ki mode I dynamic stress intensity factor in a mixed mode case Kiaux dynamic stress intensity factor due to pure mode I auxiliary field Ku mode II dynamic stress intensity factor in a mixed mode case Kuaux dynamic stress intensity factor due to pure mode II auxiliary field Km mode III dynamic stress intensity factor in a mixed mode case ' K* element equivalent stiffness matrix 'K" element convective stiffness matrix due to A L E xn 'KM element first convective stiffness matrix due to A L E ' KA2 element second convective stiffness matrix due to A L E 'KL element'Lagrangian stiffness matrix 'K L 1 element Lagrangian material stiffness matrix ' K L 2 element Lagrangian geometric stiffness matrix I element length ' M A element convective mass matrix due to A L E ' M l element Lagrangian mass matrix m unit vector normal to ro contour m. components of the unit vector normal to ro contour N number of element nodal points n unit vector normal to boundary ' ni components of the unit vector normal to boundary surface q arbitrary smooth function r polar coordinate 'SA] element stress matrix related to first A L E stiffness 'SL 2 element stress matrix related to Lagrangian stiffness 'T kinetic energy density 'f time derivative of kinetic energy density t time t normalized time 'U stress work density xiii time derivative of stress work density element incremental material displacement vector element incremental grid displacement vector element incremental Lagrangian displacement vector element nodal incremental material displacement vector element nodal incremental grid displacement vector element nodal incremental arbitrary displacement vector element nodal incremental Lagrangian displacement vector components of incremental material displacement vector components of incremental grid displacement vector components of nodal incremental material displacement vector components of nodal incremental grid displacement vector components of nodal incremental Lagrangian displacement vector element nodal incremental material velocity vector element nodal incremental grid velocity vector element nodal material velocity vector element nodal grid velocity vector prescribed impact velocity components of material point velocity vector components of grid point velocity vector components of instantaneous surface velocity xiv 'v,c components of instantaneous crack tip velocity V instantaneous crack tip velocity along the local xi direction vc steady crack tip velocity along the local x\ direction 'v. components of material point acceleration vector 'v\ components of material point referential acceleration vector 'Wex' external work ' x element coordinate vector rx element nodal coordinate vector Xf components of material point initial coordinates Xf components of grid point initial coordinates ' xt components of material and grid points position vector ' xf1 components of material point mapping function 'xf components of grid point mapping function 'xik components of nodal coordinate vector a mesh motion parameter vector due to the Coupled Displacement approach a* mesh motion parameter vector due to the True Boundary Tracking approach at mesh motion parameter vector B mesh motion parametric matrix R Newmark integration parameter P{j) components of mesh motion parameter vector ro Lagrangian crack tip contour xv T vanishing crack tip contour T+ portion of the vanishing crack tip contour on the upper crack face 'r~ portion of the vanishing crack tip contour on the lower crack face At time increment Atc critical time step referred to the Courant condition for finite elements 8 virtual or variation S(j components of Kronecker delta sl equivalent strain boundary curve y Newmark integration parameter v Poisson's ratio 6 angle of rotation of the local axes w.r.t. the global axes ' p material density cr0 applied normal stress ' a element stress vector a equivalent stress a initial yield stress 'ay components of Cauchy stress tensor '<jy components of Truesdell stress rate tensor xvi Acknowledgments I wish to express my deepest gratitude to my thesis supervisor, Professor M.S. Gadala, for his support, guidance and encouragement throughout the course of this research. It has been a great opportunity for me to be able to learn from his knowledge and experience. I am also grateful to Dr. H . Bayoumi for providing me with his source code and formulations, which constitute the cornerstone of this research. I am forever indebted to my parents in Libya for their lasting love and prayerful support. To them I dedicate this thesis. I also would like to express my appreciation to my wife, Aziza, for her patience, encouragements and all the help she offered me throughout the completion of this work. I would also like to thank the rest of my family in Libya for their prayers and emotional support. xvii Chapter 1 INTRODUCTION AND BACKGROUND 1.1 M O T I V A T I O N Most structural materials and mechanical components contain or will develop cracks, either as natural imperfections or as a result of fabrication and assembly processes. A short crack may grow in a stable manner through mechanisms such as fatigue, corrosion cracking, and so on, and might under certain conditions, become unstable. An unstable crack propagates dynamically in a structure, at velocities in the order of the materials elastic shear wave speed. This rapid separation of material would seriously damage the integrity and function of the structure, and may eventually lead to total structural failure. Interest in dynamic fracture problems has been increasing. This is because structural designs, which preclude fracture instability under all conditions, can be far too costly, and, in addition, there are applications where large scale unstable crack propagation would have catastrophic consequences. With dynamic fracture analysis, fracture failures of structures, such as: pressurized containers, nuclear reactors pressure vessels, transmission pipelines, off-shore oil production platforms, aircrafts and spacecrafts, welded ships, railroad tracks, bridge girders, and so on, may be prevented, and propagating cracks can be arrested before losing the structural integrity of the unit. 1 Chapter 1. Introduction and Background 2 Due to the complex nature of dynamic fracture problems, available analytical solutions have required the use of simplifying assumptions to render the problem more tractable. Since most realistic problems are much more complex, consistent numerical approaches are, therefore, required to model this class of problem. 1.2 OBJECTIVE In general, the problem of dynamic fracture is complex, due not only to the interaction between the reflected boundary waves and the crack tip asymptotic fields, but also to the complex nature of the boundary condition of a propagating crack. The continuously changing material domain of a body and its prescribed tractions and displacement boundary conditions, due to crack propagation, are neither a material-related type of boundary condition, i.e., they may not be consistently described in a Lagrangian formulation, nor a fixed frame-related type of boundary condition, i.e., they may not be consistently described in an Eulerian formulation. The objective of this work is to model the process of crack propagation utilizing a dynamic Arbitrary Lagrangian Eulerian (ALE) formulation, which may consistently describe the complex nature of the boundary condition involved. The resultant formulation is then to be implemented into a finite element modular program (ALEFR) or amended to the existing Lagrangian-based commercial codes. Chapter 1. Introduction and Background 3 1.3 BACKGROUND 1.3.1 Numerical Modeling of Dynamic Crack Propagation The problem of dynamic crack propagation has been widely addressed in literature. Analytical solutions are available for simple geometries and loading conditions, but they are generally of limited use in practical problems. Despite the fact that, fundamentally, the crack propagation problem is not a Lagrangian type of problem, Lagrangian models have been widely used in the numerical analysis of crack propagation. In particular, the node release method (Anderson, 1973) has enjoyed some popularity because it is very robust and easy to implement in commercial finite element codes. In this method, the crack remains stationary as a Lagrangian formulation is used to advance the solution over each time step. The crack is then advanced between time steps by splitting the old crack-tip node into two distinct nodes and advancing the crack-tip location by one element length at a time. In order to improve this discontinuous modeling, "gradual-node-release" techniques were proposed (Keegstra et al., 1978; Kobayashi, 1979; Hodulak et al., 1980; Caldis et al. 1979). In these models, the "holding-back" force was introduced to prevent discontinuous jumps of the crack-tip between two nodes. Mahanty and Maiti (1990) studied the stable crack growth using the node release method for mixed mode loading cases. This work involves the search for crack propagation direction prior to application of finite element analysis. However, the crack is then assumed to follow the same straight line along this direction. Chapter 1. Introduction and Background 4 The node release method often produces inaccurate results when explicit time integration is employed, since the prescribed nodal force is not generally in dynamic equilibrium with the state of stress in the adjoining finite elements. Additionally, the effect of impulse stress waves generated by the instantaneous application of the prescribed nodal force reduces the accuracy of the method. Further details are given in the critical study of Kanninen (1987). Bazant, et al. (1978) have proposed a special finite element procedure based on a moving coordinate system centered at the crack tip. The finite element mesh is translated as a rigid body to follow the crack tip motion. However, this procedure is restricted to semi-infinite strips whose surfaces are parallel to the direction of crack propagation. Shephard et al. (1985) presented a method that requires an entire structure remeshing at each crack increment. Valliappan and Murti (1985) introduced a window that moves with the crack tip, such that only the elements within this window are candidates for remeshing, however, this approach is limited to crack propagation through a region with a very regular mesh. Similarly, Wawrzynek and Ingraffea (1989) presented a technique where only local remeshing takes place at each crack increment. This technique was implemented in the finite element program F R A N C (FRacture ANalysis Code). An alternative approach to modeling running cracks is the moving element approach. Nishioka and Atluri (1986) have employed this concept using a singular hybrid element. Chapter 1. Introduction and Background 5 In their method, a special singular element is formulated using analytical asymptotic solutions for the singularities near the crack tip together with the variational statement that accounts for the motion of the singular element. The crack propagation process is simulated by remeshing the conventional finite elements around the crack tip at each increment. The main problem with the hybrid approach in general, is the complexity of formulation and the difficulty in implementation. In addition, the convergence behavior of such formulation is unpredictable. Most of the finite element methods that employ local remeshing use a Lagrangian kinematic model, in which the geometric discretization is selected a priori for material configuration, and the finite element mesh is required to follow the material motion during a time step. This process requires local remeshing every time step, and frequent overall remeshing; therefore, it induces errors due to the interpolation of the displacement and velocity fields in each remeshing. In addition, the numerical techniques used in these procedures are not compatible with the general Lagrangian-based finite element codes. 1.3.2 The Arbitrary Lagrangian Eulerian (ALE) Formulation In the field of computational mechanics, there are two classical formulations used to model continuum problems. These are: the Lagrangian (referential) formulation and the Eulerian (spatial) formulation. In the Lagrangian formulation, an initial or updated configuration is used as a reference frame. In a finite element method, based on this formulation, the computational grid is fixed to the material points of the deformed body. Chapter 1. Introduction and Background 6 Accordingly, the Lagrangian-based finite element method may consistently handle material-associated boundary conditions. On the other hand, the reference frame (the computational grid) in the Eulerian formulation is fixed in space. Eulerian formulations are generally suitable for flow type problems and steady state conditions. Since the process of crack propagation physically involves the creation of new surfaces by splitting material points, at the instantaneous location of the crack tip, the one-to-one mapping requirement in the Lagrangian formulation does not exist; that is, a single material point in a configuration before crack propagation may split into two material points in a configuration after crack propagation. Therefore, from a theoretical point of view, this type for formulation may not truly model the process of crack propagation. On the other hand, the only approach to model the process of crack propagation using the Eulerian formulation is to regard the process as a constrained flow of material. This nonphysical assumption precludes the use of this type of formulation in accurate modeling of crack propagation problems. Although neither of the above classical formulation are capable of simulating the crack propagation process, we believe that a combination of the two may consistently represent this process. This combination, which is known as the Arbitrary Lagrangian Eulerian (ALE) formulation, has emerged in recent years to alleviate many of the drawbacks of Chapter 1. Introduction and Background 7 the classical formulations. In the A L E formulation, it is not necessary for the reference frame (or the finite element mesh) to be adhered to the material or to be fixed in space, rather, it may move arbitrarily. Consequently, mesh distortion associated with large deformation, may be avoided by designing a proper mesh motion scheme. In this work, the merit of the Lagrangian formulation, in that it can easily describe the material-related boundary conditions, and the merit of the Eulerian formulation, in that it can easily model the process of the branching of a steady flow, are combined to consistently model the process of dynamic crack propagation. 1.3.3 History of the A L E The A L E formulation was first proposed to model fluid mechanics applications using finite difference (Noh, 1964; Hirt et al., 1974) and later introduced to finite element analysis (Hughes et al., 1981). The A L E was first introduced to the analysis of solid mechanics application by Huetink (1982). He modeled a quasi-static metal forming problem employing a technique referred to as operator split. In this technique, material deformation and convective effects are treated separately although they are coupled in the same equations. Thus, each time step may be split into two steps; a pure Lagrangian step, at which the computational grid moves with the material; followed by an Eulerian step, at which the Lagrangian solution to the reference grid and stresses are updated using convective effects. Due to the fact that the operator split technique is computationally efficient and easily implemented into the common Lagrangian finite element codes, it has been adopted in the majority of A L E analyses. In his work, Chapter 1. Introduction and Background 8 Huetink used an approximation method to update integration point variables. This was carried out first, by calculating nodal variables utilizing integration points from all elements sharing that node. A continuous field for each variable is obtained by interpolating nodal point values using element shape functions. Similarly, gradients of variables are evaluated using shape function derivatives. Haber (1984) introduced a different type of A L E formulation, termed Eulerian Lagrangian Description (ELD). He divided the total deformation in each increment into separate Eulerian incremental displacements and Lagrangian incremental displacements. An Eulerian deformation gradient defines the mapping from the initial configuration to the reference configuration, while a Lagrangian deformation gradient describes the mapping from the reference configuration to the current configuration. The product of these two gradients gives the total deformation gradient from the initial to the current configuration. The E L D was then used to model crack propagation (Koh et al., 1988). However, the existence of two sets of unknown displacements makes it difficult to relate other formulations and preclude it from implementation into the existing Lagrangian codes. In addition, crack propagation was modeled artificially by merely advancing the crack tip node, without treatment of the separated material points and the nodal motion on the newly created surfaces. The method therefore, is applicable only to self-similar crack propagation cases, at which the errors due to this inconsistent modeling are minimal. Chapter 1. Introduction and Background 9 Liu et al. (1987, 1988) derived implicit quasi-static A L E formulation. In their work, the operator split approach was replaced by a fully coupled or, 'unsplit', approach. In the fully coupled approach, the governing A L E equations, in which the material deformation and convective effects are coupled, are implemented and solved without decoupling. This coupled approach is more accurate and theoretically represents a true kinematical description. In this work, Liu et al. (1987, 1988) used mixed explicit and implicit calculations to handle the convective effects. Therefore, the resulting A L E equations are difficult to relate to the incremental form of displacement used in the common Lagrangian codes. In another development, Huerta and Casadei (1994) presented a dynamic A L E formulation based on a fully explicit calculation. The convective terms are treated in a manner similar to the one developed by Liu et al. (1987, 1988). However, explicit solution schemes in dynamic analysis generally suffer from the lack of generality of application due to the stringent stability conditions, which necessitate the use of very small time steps. A fully coupled A L E quasi-static formulation was developed by Wang and Gadala (1997). The convective terms are treated using an approach similar to that of Liu et al. (1987, 1988). However, the implementation of the developed formulation does not represent a strictly coupled approach in the sense that convective terms are not determined within the iteration of each load increment. Rather, the material associated properties are updated after convergence. This approach, which was originally Chapter 1. Introduction and Background 10 formulated to model metal forming processes, has been extended by Movahhedy (2000) to include thermal and contact modules to model metal cutting process. Bayoumi (2000) developed a fully coupled quasi-static and dynamic A L E formulation for metal forming applications. The convective terms are treated using a new approach derived from continuum mechanics laws, instead of the approximate approaches, which involve unjustified assumptions. In their work, an implicit solution scheme was used and the only independent variables are the displacements, as opposed to velocities, which is common practice in the A L E literature. 1.4 SCOPE OF WORK In this work, the dynamic fracture problem will be modeled by a fully coupled dynamic A L E formulation. The scope of work may therefore be summarized in the following: - the derivation of fully coupled A L E virtual work equations based on principles of continuum mechanics, through which, the crack propagation process may be consistently modeled; - the discretization of the A L E virtual work equations using isoparametric finite elements; - the modeling the process of crack propagation as a process of continuous separation of material points; Chapter 1. Introduction and Background 11 - the development of an efficient procedure to control the motion of the finite element nodes on free boundaries, including those that evolved due to crack propagation; - the development of a numerical procedure to continuously evaluate fracture parameters and control crack behavior; - the implementation of the discretized A L E equations and the above techniques into a modular 2-D finite element code; and - the testing of the developed code through the simulation of stationary and propagated dynamic fracture problems with known analytical solutions. 1.5 ORGANIZATION OF THESIS The derivation of the A L E equations is given in Chapter 2. This includes the A L E quasi-static equilibrium equation and the A L E dynamic equation of motion. The derivation of dynamic energy release rate integrals from the energy balance equations is shown in Chapter 3. The correlation between the dynamic stress intensity factors and the dynamic energy release rate is also discussed. The development of an ALE-based numerical techniques including: the Material Point Splitting (MPS), a technique for modeling the process of crack propagation; the True Boundary Tracking (TBT), a technique for imposing the A L E boundary constraint on the old and the developed free surfaces; and techniques to evaluate the dynamic energy Chapter 1. Introduction and Background 12 release rate integrals, are shown in Chapter 4. The implementation of these techniques into a finite element program is shown in Chapter 4. Testing of the developed code is covered in Chapter 5. Tests include static and dynamic Mode I crack propagation cases. A sample mixed mode case is also presented. Finally, the conclusions drawn from this numerical procedure, along with recommendations for future work, are presented in Chapter 6. Chapter 2 FORMULATION AND DISCRETIZATION OF A L E EQUATIONS 2.1 A L E GOVERNING EQUATIONS 2.1.1 Preliminaries Linearized virtual work equations for quasi-static and dynamic applications are derived to be solved in an implicit time-stepping approach. Throughout the derivation in this section, notations of time and configuration adopted are similar to those used by Bathe (1996). Left superscripts indicate the configuration in which the quantity occurs, whereas left subscripts indicate the configuration to which the quantity is referred. Omitted left subscripts imply that the quantity occurs in the same configuration in which it is measured. For an incremental quantity from time t to t + At, the left superscripts and subscripts are omitted. As in the standard indicial notation, right subscripts refer to tensor or vector components, repeated indices imply summation over the admissible Q r j range and the notation [ ];, means . dxi In the A L E description, the motion of particles is referred to an arbitrarily assigned reference system which is neither adhered to the material, as in the Lagrangian system, nor fixed in space, as in the Eulerian system. This means that the material configuration at any time t refers to the set of material particles, whereas the reference configuration 13 Chapter 2. Formulation and Discretization of the A L E Equations 14 consists of a set of arbitrarily moving grid points (i.e. finite element nodes) sharing a common boundary with the set of material particles. The material configuration is identified by a set of material point coordinates ' X ( m while the reference, or grid, configuration is identified by an independent set of grid point coordinates 'Xf . Let 'x™('XJ,t) and 'x , s ( 'XJ ,t) be the vector functions or the mappings that characterize the motion of the material point ' X J 1 and the grid point ' X J in space, respectively. The position of 'XJ at time t is given by .'*,.='x;('x;,o (2.i) The set of material particles is related to the set of grid points by requiring that the two configurations share the same space at all times. Any point within the common boundary is occupied by elements of the two sets. Thus, the position of the grid point ' X J that occupies the same point in space at time t as ' X J 1 is also given by 'xt as 'x^x?CX«,i) (2.2) The A L E formulation requires that the inverse of Equations (2.1) and (2.2) exist to ensure a one-to-one mapping between the two configurations. The material velocity 'v,, and the grid point velocity 'vf, at time t are given by Chapter 2. Formulation and Discretization of the A L E Equations 15 The boundary constraint, which ensures that the material and grid configurations have the same boundary at all times, may be expressed as (V'v/)'n,.| =0 (2.5) Ion the boundary where ' n{ is the unit normal to the boundary surface. The governing A L E equations involve the material time derivative of several quantities. The material derivative of an arbitrary function ' / is denoted by a superposed dot and is defined to be the rate of change of the function holding the material particle 'Xf fixed df 7 = a (2.6) <x? However, the grid configuration is the computational configuration that tracks the history of all quantities. Thus, it is convenient to define a grid time derivative, which is the time derivative of the function ' / holding the grid point 'Xf fixed, and denoting it by a superposed prime df f' = a (2.7) The relation between the two time derivatives is given by (Hughes, 1981) OX; The A L E formulation will be discretized using the isoparametric displacement based finite element method. The incremental material displacements and the corresponding Chapter 2. Formulation and Discretization of the A L E Equations 16 incremental grid displacements from time t to time t + At are respectively denoted by ui and uf . The position of the grid point in the configuration at time t + At may be expressed as (2.9) The local form of the conservation of mass, i.e., continuity, at time t is given by 'p=-'p-r- (2.10) d xt where 'p is the material density. Using Equation (2.8), the continuity equation with respect to an arbitrary moving grid point may be expressed as ' P ' = - ' p ^ - ( \ - ' v ! > ^ (2.11) oxi oxi 2.1.2 Quasi-static Analysis In processes like stable ductile crack growth or low-speed metal forming, the effects of inertia forces may be neglected. Since body configuration at time t + At is yet unknown, an approximate solution may be obtained by referring all variables to the grid configuration at time t and linearizing the terms of equilibrium equation. The solution is then refined by iterations within an implicit incremental approach. Stress and strain measures are chosen to accommodate large strain applications. Chapter 2. Formulation and Discretization of the A L E Equations 17 To express the equilibrium of the body at time t + At, the principle of virtual displacements is used. It may be expressed as j'^Vijd^ij '+A'dV = $'+Al p'+%BSUi ,+AtdV + j'+'f'Su, l+A,dS (2.12) where / + A ' c r is Cauchy stress tensor and ,+A(e,y is the conjugate strain tensor defined as The RHS of Equation (2.12) represents the external virtual work, S'+A'Wex', due to body force ,+A'ftB , and the applied surface traction ' + A ' / ; 5 . Linearization of virtual work equation Linearization of the terms of the virtual work equation is carried out by adopting an incremental approach. Variables at time t + At are assumed to be composed of their respective values at time t, plus an increment given by the grid time derivative of the variable multiplied by the time increment Ar. The Cauchy stress tensor may, accordingly, be decomposed into <rv+'o'vto (2.14) Using Equation (2.8), we get where, 'dy is the material rate of Cauchy stress calculated from the material constitutive Chapter 2. Formulation and Discretization of the A L E Equations 18 relation (Appendix A). The variation in the strain components is decomposed as S.^^S^+S^At (2.16) in which Sre'(i is the grid time derivative of df^ and is given by (Wang an Gadala, 1997) "J 2Kd\ d'Xj d'xk d'x, { } Substitution in Equation (2.16) yields The decomposition of material density, upon substituting in Equation (2. II), at time t + At, may take the form ,+*P=!p-'p^-{uk-ul)^- (2.19) Finally, the elements of volume and surface area are decomposed into (Malvern, 1969) ,+A'dV='dV+'dVAt = (l + ^ -)'dV (2.20) d'xk A.S 1 A.g A.g ,+A'dS='dSVdS'At = [l + ^ - - ( ^ + -^-ynm,nJdS (2.21) where 'n is the unit outward normal to the surface. Chapter 2. Formulation and Discretization of the A L E Equations 19 Upon substituting the above decomposed quantities into Equation (2.12) and ignoring the higher order terms, the linearized form of the virtual work equation may be written as \'aijSle'dV+ \'&ilAtSle'dV+ \ ^ - ! , e ' d V - \^L'a,^'dV J 'J i 'J J 'J i 'J J at 'J ' 'J J y at •V 'V 'V k 'V j k - \(uk-u*y^^'dV= \'p'+%BSUi'dV- l'p'+%\p--P-)&l'dV (2.22) >V k 'v 'V O Xk O Xk [ t + A t f B , t i SP c„ ( I T / , [t+At /-Sri , <^k , \t I l t I J O - I / , (uk — uk) — oui dV+ I / , + — _ ( _ _ + _ _ ) „ n]fa dS •i • dxk 's dxk 2 dxn 3xm The RHS of Equation (2.22) represents the external virtual work S'+A'WeM . Equation (2.22) has two convective terms that involve spatial derivatives of stresses and material density. These are the last term in the LHS and the third term in the RHS. Since, in finite element formulation, the stresses and material densities may not have the inter-element continuity, their gradients may not be reliably computed on the element level, when evaluating element matrices. Different convective term treatment approaches are available in the literature (Huetink, 1982; Wang, 1998; Liu et. A L , 1987). Most of these treatments are based on the simple interpolation of calculated quantities. Bayoumi and Gadala (2000) developed a consistent method to replace the convective terms by their equivalent non convective terms based on the fundamental continuum mechanics relations. This includes applying the divergence theorem to the integral containing the convective term and substituting the A L E boundary constraint, Equation (2.5) (Appendix B). After treatment of the convective terms, Equation (2.22) becomes Chapter 2. Formulation and Discretization of the A L E Equations 20 J ' < r ^ W + J ' a , ^ , ' d V + I-CT&.JLw-jffv„|^V •V 'V 'V UXk <VU i U k + \{uk-ufY^-^'dV = l'p'+AtfiBSu<dV+ \'p^-l^{uk-ut)8u<dV (2.23) 'V ^Xk 'v , 'V k + \'p-%\uk-u!)^L'dV+ [ ^ f ^ ^ - h ^ A n J n j S u ' d S >v d x k .\ dxk 2 8xn dxm Fully coupled A L E equilibrium equation By substituting the constitutive relations, Equation (B.l) to (B.3) in Appendix B, into Equation (2.23) the A L E equilibrium equation may be written as \'CijklleklSteij'dV+ l'<TvStrj0'dV 'V 'v ' + J K - u D ' A d V , J ( | L _ | L y f f ^ ^ (2.24) •i d x k i d x ; d*j d\ 'V i ty~. • . • . • 1 9ut du. where C,,w is the material constitutive tensor, ,77 = - — and 3 ' * 2 ffx, d'xj 8 ' + ^ ' = \'Pr%B+^-^(uk-u?)]5u-dV+ \<p'+%B(uk-ul)?P-'dV •i dxk >i sxk + ^ f n i + ^ - ^ + ^ ' n J n J S u / d S •s dxk 2 dxn dxm (2.25) Since Equation (2.24) is expressed in terms of both material and grid incremental displacements and uf, respectively, the grid may therefore, move independently relative to the material. This characteristic is very desirable in problems that involve large material deformation and/or change in boundary conditions. Chapter 2. Formulation and Discretization of the A L E Equations 21 The A L E formulation, Equation (2.24) may be seen as a general form for an incremental virtual work equation. If, for instance, the grid is adhered to the material, i.e., uf =ui, the equation will be reduced to the updated Lagrangian formulation. On the other hand, if the grid is chosen to be fixed in space, i.e., uf = 0, the equation will then represent the Eulerian formulation. 2.1.3 Dynamic Analysis In dynamic analyses, inertia forces are included as an extra term that should be added to the LHS of the virtual work equation, Equation (2.12). Using the relation between the grid time derivative and the material time derivative, Equation (2.8), the virtual work done by inertia forces at time t + At, may be written as The first term on the RHS of Equation (2.26) may be referred to as the referential inertia term, whereas the second term is referred to as the convective inertia term. Linearization of inertia terms The material point velocity ,+Atvi, the grid point velocity t+Atvf and the grid time derivative of the material point velocity l+Atai, at time t + At, may be related to their respective known values at time t as V V (2.26) t+At (2.27) Chapter 2. Formulation and Discretization of the A L E Equations 22 ,+*v*='v°+vf (2.28) ' + >='a.+a ; . (2.29) where the incremental quantities v(-, vf anda(. are controlled by the time integration method to be utilized. The incremental decomposition of the two inertia terms on the RHS of Equation (2.26) and the treatment of the convective terms, are handled in a manner similar to quasi-static analysis. The linearized referential and the convective inertia terms may be written as J ' + V + * « A ' + * < W = l'p'aiSui'dV+ j'pafa'dV f"*v 'v 'y (2.30) 5% P)t+At p)t+AtM \ '^PC\-<+\*)^±Sur'dV= \p^v-^v))0—^-8u.;dV vmv o Xj ,v O Xj , d,+Atv- dt+A'vs d'+Atv •v dxb 8xk dxj d ,dt+A'Vi (2.31) d xk dx} Fully coupled A L E equation of motion The form of the A L E equation of motion, obtained by adding the inertia terms, Equations (2.30) and (2.31) to Equation (2.24), may be written as Chapter 2. Formulation and Discretization of the A L E Equations 23 \'paiSui'dV+ [cmtekl5fij'dV+ fas^'dV + \'pCv-'v^^Su.'dV 'V 'V 'V 'V j + { V ( v J - v ? ) f i * , W + l i U t - u ! y A d v + j A - ^ y ^ f t w ,JV Sxj ,JV dxk ,JV dxj dx. dxk = S»*vr* - faS^'dV - I'p'afa'dV - \'p(!vk-'vl)^^-At'dV (2.32) 'V 'V 'V k •V U X i 'V k U Xk U X j 'V O Xk O Xj The first three terms on both sides of Equation (2.30) represent the updated Lagrangian formulation. The first term on the LHS of the equation corresponds to the Lagrangian mass matrix. The second, third, sixth and seventh terms on the LHS and the first and second on the RHS were already defined for quasi-static analysis. The fourth and fifth terms are the convective velocity stiffness matrices due to material/grid relative motion. The third term on the RHS corresponds to the Lagrangian inertia force vector, whereas the last four terms on the RHS are convective inertia force vectors due to material/grid relative motion. 2.2 FINITE ELEMENT EQUATIONS 2.2.1 Isoparametric Finite Elements In finite element discretization, element coordinates 'xi, and incremental displacements ui and uf , are correlated to their respective nodal values, !xik, uik and ufk, through the element shape function hk, i.e. Chapter 2. Formulation and Discretization of the A L E Equations 24 k=\ N k = l N (2.33) (2.34) (2.35) where i corresponds to the degree of freedom, k is the nodal point number and N is the number of nodal points in the element. For a two dimensional space, Equation (2.31) may be expanded as or x=H'x where 'x is the element coordinate vector given by y H is the element shape function matrix in the form H K o ; K \ 2x2N and ' x is the nodal coordinate vector given by Similarly, Equations (2.32) and (2.35) may be written in the form u = Hu (2.36) (2.37) (2.38) (2.39) (2.40) (2.41) Chapter 2. Formulation and Discretization of the A L E Equations 25 u* = Hu 4 (2.42) where u, u8, u and us are the incremental element and nodal material and grid displacement vectors, respectively. 2.2.2 Discretization of the Quasi-static ALE Equation The variables in the linearized A L E equilibrium equation, Equation (2.24), have to be discretized in order to develop the finite element equilibrium equation for quasi-static analysis. Considering the internal force term, the last term on the RHS of Equation (2.24), the stress may be arranged in a vector form as (2.43) Similarly, the incremental strain vector has the form dux d'x t exx duy teyy > — < d'y dux duy — - + — y -t^zz J d'y d'x . 'x (2.44) which upon using the element displacement expansion of Equation (2.32), gives Chapter 2. Formulation and Discretization of the A L E Equations 26 ,e = < , h&xU* h d'y yk ttKd'y xk d'x yk> N h Z J N Uxk 7=1 d'x 0 d'y K 2>A 0 d\ d'y d'x 0 (2.45) or ,e=B"u (2.46) where tBLX is the element shape function derivative matrix related to Lagrangian material stiffness and it may be written as ,BLl = d'x 0 d'y K 0 dK d'y d\_ d'x 0 (2.47) 4x2N The finite element form of the internal work term can be established by using Equations (2.43) and (2.46) as follows \'cjijSteij'dV = \(Sle)T 'o'dV 'V 'V = $(lBLldu)T'<i'dV 'V = (du)T \{tBu)T'v'dV 'V = (<5u)T'f (2.48) Chapter 2. Formulation and Discretization of the A L E Equations 27 where ' f is the internal force vector given by 'f = \(tBu)T'a'dV (2.49) •v The Lagrangian material stiffness virtual work term may be rewritten in matrix form as \cmtek]5fi;dV= l(Ste)T'C,e'dV 'V 'V = \(tBudu)T'C(tBL1u)'dV •v (2.50) = (du)T \{tBn)TtC,BLUdVu 'V where ' C is the material constitutive matrix and 'KLl is the Lagrangian material stiffness matrix given by (Bathe, 1982) 'Ku = §(tBu)T'C,BLUdV (2.51) 'v The Lagrangian geometric stiffness virtual work term may be rewritten as ^ ' d V ^ ' v ^ f ^ ' d V (2.52) V -v dxj dxi [(TyS^'dV = (<5uf \{rBL2)TtSL2tBL2tdVu •v -v (2.53) = (<Su)r'KL2u where Chapter 2. Formulation and Discretization of the A L E Equations 28 B L 2 =| d'x o% d'y 0 8K d'x d'y 0 (2.54) <SL2 = <?xx 0 <7 'a 0 xy yy 0 0 'a XX 0 0 'CJ xy 0 0 0 0 0 J5x2/V 0 " 0 0 0 '<x„ (2.55) 5x5 and the nonlinear geometric stiffness matrix ' K " is given by (Bathe, 1996) , K L I = J ( ( B " ) r ' S " , B " ' d V 'V The first convective stiffness virtual work term due to A L E may be discretized as d5£„ (2.56) \(uk-uD'cjy—^-'dv = (saf j(tBA1)TtsMn'dV(u -u«) ox, ,y = (Sa)T,KM{u-ug) 'V (2.57) where Chapter 2. Formulation and Discretization of the A L E Equations 29 ,B A] d\ d'x2 d\ d'y2 d\_ d'xd'y 0 0 0 0 (7 0 0 8% d'x2 d% d'y2 d'xd'y 1 dK K ) x d'x ' x2 'x d'y (2.58) J8x2W tSM = 0 0 0 yy yy xy 0 ' 0 CT . 2x8 and the first convective stiffness matrix due to A L E is given by , K A I = J ( ( B > Y ' S / U H ' d V •v (2.59) (2.60) The second convective stiffness virtual work term due to A L E may be discretized as r{a±_a£_ytj^k'dv = (saf [(tBA2)T'SL2BL2'dV(u-us) „J. dx, dx, A* J d'x, •v = ( < S u ) r ' K A 2 ( u - u s ) (2.61) where Chapter 2. Formulation and Discretization of the A L E Equations 30 d'x 0 0 d'x d'y 0 (2.62) 0 dhk d'y K 0 N 15x2 N and the second convective stiffness matrix due to A L E is given by 'KA2= \(tBA2f'SL2BL2'dV (2.63) Substituting all of the above discretized terms in Equation (2.24) the final A L E finite element equilibrium equation for quasi-static analysis may then be written as where; 'KLl and 'KL1 are the Lagrangian material and geometric stiffness matrices, ' K A 1 a n d 'KA2 are the A L E convective stiffness matrices and ,+AlVx' and ' f are the external and internal force vectors, respectively. 2.2.3 Discretization of the Dynamic A L E Equation The A L E equation of motion for dynamic analysis, Equation (2.32) is discretized to develop the finite element equation for dynamic analysis in a manner similar to that used to drive the quasi-static finite element equation. Thus, for isoparametric finite element discretization, velocity and acceleration vectors may be written as ('Ku+'KL2)u+('KM+'KA2)(u-ug)='+A'fex'-'f (2.64) Chapter 2. Formulation and Discretization of the A L E Equations 31 ' v = H v , V ' = H'v« v = H v , v* = H v s (2.65) a = H a , a = Ha where 'v, 'vg and 'a represent the nodal material velocity, grid velocity and referential material acceleration at time t, whereas v, v4' and a are the nodal incremental quantities of the same variables from time t to t + At. Using Equation (2.60), the Lagrangian mass term may be discretized as follows l'paiSui'dV= l'p(HSu)T(H&)'dV 'V 'V = (du)T \'pRTYl'dVa (2.66) •v = (<5u)r'Mz-a where '~ML is the Lagrangian mass matrix given by ' M L = \'pHTH'dV (2.67) •v The inertia force virtual work term is discretized as j'p'aiSul'dV= J'p(H<5uj7'(H'a)'dV' 'V 'V = (du)T j'pRTH'dV'a (2.68) •v = (<5u)r'ML'a The first inertia force virtual work term may be discretized as \'p(vk-'vgk)d{'a;SUi) At'dV = (Sa)T'MA'a (2.69) where 'MA is given by Chapter 2. Formulation and Discretization of the A L E Equations 32 'MA = 'MA 'M~A 1 Y 1 2i-l,2./-l I Y I 2i-\,2 j 'MA 'MA 2Nx2N in which i and j indicate node numbers from 1 to N, and >MA 1 Y 1 2i-1,2y-1 J ox ox •v k=\ dh 7)h N dy dy $^ 'MA ='MA =0 m2i-l,2j 1VI2i,2j-\ U (2.70) (2.71) (2.72) The discretization of the convective terms in the equation of motion is shown in Appendix C. Upon substituting all of the discretized terms into Equation (2.35), the dynamic A L E finite element equation may then be written as ' M ^ ' C V C ^ v - y ^ f K M ^ u + C K M ^ X u - u 8 ) = l + A , f e X , _ t f _(<ML+'MA)'a-( rCA 1+'CA 3+'CA 4)'v (2.73) Chapter 3 DYNAMIC FRACTURE PARAMETERS 3.1 I N T R O D U C T I O N According to linear electrodynamic fracture theory, the near-tip stress field is dominated by the singular term in the asymptotic expansion of the solution (Freund, 1990). For each fracture mode, determining this singular stress field may be reduced to determining the amplitude factors known as the dynamic stress intensity factors K,, Kn and Kw. These factors are crucial in assessing the subsequent crack behavior in terms whether or not the crack will propagate, in which direction and at what speed. This chapter introduces the dynamic fracture parameters for stationary and propagating cracks. The techniques and solutions in this work are primarily for stationary and propagating mode I (opening mode) crack problems. In addition, the in-plane mixed mode crack problem, that is, combined mode I and mode II (shearing mode), is briefly addressed. The dynamic stress intensity factors may be extracted from the energy release rate G defined in domain (area) and contour (line) integral forms. The energy release rate is derived from a mechanical energy balance in the area surrounding the crack tip. The procedure closely follows that of Freund (1990) and Moran and Shih (1987). 33 Chapter 3. Dynamic Fracture Parameters 34 3.2 GENERAL MECHANICAL ENERGY BALANCE In the absence of body forces, the Lagrangian form of equation of motion may be written as 'atJJ = 'p 'v. (3.1) By taking the inner product of Equation (3.1) with the material point velocity v, , the result is \ = 'P 'Vj 'V ; (3-2) Since, {'<Jy'vj) t= aiUi 'v i + 'ay 'vy._; , then Equation (3.2) may be rewritten in the form ^y'v.X^'cJy'v^'p'v./v, (3.3) Since the stress work density, 'U , and the kinetic energy density, 'T , at a material point are given by 'U= J'CT.. ''v..dt' (3.4) f'=0 'T = \'p'vi'vl (3.5) Equation (3.3) can be written as, ('ay 'vjl^'U+'f (3.6) Equation (3.6) is referred to as the general balance equation and is valid for any material response. Chapter 3. Dynamic Fracture Parameters 35 Equation (3.6) may be integrated over a general three-dimensional body of volume 'V bounded by the surface 'S with a unit outward normal m. Let 'vf be the instantaneous velocity of the surface 'S. This yields | ( ' f f / v J / d V = K'U+'t)'dV (3.7) 'V 'V Applying the divergence theorem to the LHS of Equation (3.7) and Reynolds transport theorem to the RHS, the mechanical energy balance may be written as J V,/v.m,. 'dS = — j('U+'T) 'dV-l('U+!T)'v?m, 'dS (3.8) •s  d t 'V 's 3.3 E N E R G Y B A L A N C E F O R C R A C K TIP R E G I O N Consider an arbitrary two-dimensional body, Figure (3.1) containing a crack and bounded by a Lagrangian (material associated) curve ro. The crack is assumed to be propagating with an instantaneous velocity 'vc along x;-axis. The crack tip region is bounded by a small contour 'r that starts at one of the traction free crack faces and ends at the other. This contour is fixed in shape and it translates with the crack tip at the same velocity. Integrating Equation (3.8) over the area 'A, bounded by the contours r„ , 'r, lT+ and *r~ and applying the divergence theorem to the LHS of the equation and Reynolds transport theorem to the RHS yields (Moran and Shih, 1987) Chapter 3. Dynamic Fracture Parameters 36 Figure (3.1) Crack tip contours. j'a^Vjm, 'dC = — 1('U+'t) 'dA- |('C/+'r) \m{ 'dC (3.9) where 'C = F„ U 'F+ U 'r~\J 'r. By assuming that crack faces are traction free, i.e., the stress vector 'cr.m. = 0, Equation (3.9) may be expanded as J'cr..' y. m, 'dC = — J('U+'t) 'dA - ^('U+'t) V m,. ] 'dC r0 D T - A T+UT- (3.10) - j[('C/+'r) V/ra1+'o-1/v(mI.] 'dC 'r The term on the LHS of Equation (3.10) is the rate of traction work being inputted into the body. The first term in the RHS represents the rate of increase of the internal energy, whereas the second term represents the rate of energy lost due to flux through crack faces. The last term on the RHS represents the rate of energy lost due to flux through 'T, which will be denoted by F. Since m = -n on 'F, F can be expressed as Chapter 3. Dynamic Fracture Parameters 37 F(T)= {[('[/+'r)'vcJ1,+'o-,/v.]n,. 'dC (3.11) •r 3.4 T H E D Y N A M I C E N E R G Y R E L E A S E R A T E The dynamic energy release rate G is defined as the energy released from the body due to a unit crack growth. The dynamic energy release rate may be related to the contour 'F, as it is shrinking into the crack tip, as G= lim J ^ - ^ L lim j - ? - f f'cr, V . +(U+'T) 'vcSu]nt 'del (3.12) In order for G to be a fracture parameter of fundamental significance, its value must be independent of the shape of the contour 'F. To examine the path dependency of the above contour integral, consider a closed loop formed by two crack tip contour integrals. Upon application of the divergence theorem to the energy flux integral, Equation (3.11), yields F(ri)-F(T2)= { 'An where 'AJ2 is the area within the closed path. In general, the above integrand is not necessarily zero and consequently the value of F will be path dependent. Consider, as a special case, a steady state crack growth ( V = v c = const.). In this case, any field quantity/^) depends on xj and t only through the combination %=xl-v°t. This causes the area integral of Equation (3.13) to vanish, which implies that the energy d'U dt •+ v d'U , d'u + p dx. fd2!u: , r d2'u: dt dt2 dtdx. 'dA (3.13) Chapter 3. Dynamic Fracture Parameters 38 release rate contour integral, Equation (3.12), is independent of the path f7"for steady state crack growth. The relation between material point velocity and crack tip velocity may be expressed as , a x _ v a y ( 3 1 4 ) dt fix. For steady state crack growth, the first term on the right hand side of Equation (3.14) vanishes. For the non-steady state conditions, this term is not necessarily zero; however, in areas close to the crack tip, where the displacement gradient is very large, the second term dominates and Equation (3.14) may be written as ' v , - V ^ (3.15) Using Equations (3.14) and (3.15), the dynamic energy release rate may be expressed as, G = \[(!U+,T)5u-tCTij 'u^-'dC (3.16) 'r where G is path-independent for steady crack growth and locally path-independent for non-steady crack growth. 3.5 D O M A I N I N T E G R A L S F O R E N E R G Y R E L E A S E R A T E In finite element dynamic analysis, it is difficult to reliably evaluate energy functions defined in path forms even for steady state problems. Therefore, Equation (3.16) will not provide reliable results if it is implemented in a finite element modeling of dynamic crack problems. Chapter 3. Dynamic Fracture Parameters 39 To overcome this difficulty, Moran and Shih (1987) applied the divergence theorem to Equation (3.16), considering the configuration in Figure (3.1) for a crack with straight faces, and taking the limit as lr—> 0. This alternative expression for G is expressed as G = - {[( 'U+'T)Su - '<r„ ] q,j 'dA- \[p \ %-'p 'a,. \ x ] q 'dA (3.17) 'A 'A where q is an arbitrary smooth function that is a unity on '/"and vanishes on r0. Organ (1996) added an extra term to Equation (3.17) to solve for general curved cracks, i.e., G= - J[( 'U+'T)SU - 'ay \ ] q j 'dA- \[p 'v, \ - ' p 'a,. \ ] q 'dA ' A ,/ \ (3.18) + jyU+'T)qm] 'dC T+\j'r-Other expression for G for straight cracks derived from Equation (3.16) include (Nishioka and Atluri, 1983) G= \lu+'T)n-'crymj X j 'dC- \\p \ \ - ' p \ \ ] 'dA (3.19) 'A and G = Hm | | 'p('v, c) 2«., uif] nx 'dC + J['t/n,-'cr.. m} 'uiX] 'dC + J|> 'a,. 'u{,] 'dA (3.20) ' r r 0 '/t In Equation (3.17) to (3.20), 'A represents the area bounded by the contour i~o. The first integral in Equation (3.20) is evaluated at the crack tip. Chapter 3. Dynamic Fracture Parameters 40 3.6 M I X E D M O D E STRESS INTENSITY F A C T O R S The dynamic form of Irwin's relationships that relates the energy release rate to mode I and mode II stress intensity factors in plane strain condition is expressed as G(vc,t)=l-^-[A, (vc)Kf+A„ (vc)Kf,] (3.21) where, vc is the crack tip speed and Ai (yc) and An (yc) are universal functions of crack tip speed and are given by .2 .2 A ( V ) = / X , A„(vc)= (3-22) (\-v)csD (l-v)csD where c-2 c2 I 1 V 1 V 2x2 (3.23) D = 4ad as-(\ + as) and cs and c</ are the elastic shear and dilatational wave speeds respectively. Equation is obtained by substituting the asymptotic crack field solution into Equation (3.16). For mode I problems, determining the dynamic stress intensity factor, Kh from the dynamic energy release rate G, is straightforward, i.e., substituting Ku = 0 in Equation (3.21) . However, in mixed mode problems, Equation (3.21) does not provide Ki and Ku as separate quantities. This may be resolved by using extraction formulas based on known auxiliary field solutions (Shih and Asaro, 1988). Any field that satisfies the equation of motion, Equation (3.1), may be superimposed on the actual field. Chapter 3. Dynamic Fracture Parameters 41 Considering a field equal to the first term in the asymptotic solution for pure mode I case with unit stress intensity factor, i.e., K,aux = 1, and Kuaux = 0, the energy release rate in Equation (3.21) becomes • : l-v2 E A {K,+KIUJ2+A„Kl, AD AU (3.24) where, G^al is the energy release rate of the total field (the actual field plus the auxiliary field), G ^ i s the energy release rate of the auxiliary field, and G-^ is termed as the interaction energy release rate. The superscript (1) implies that the auxiliary field is a pure mode I type. Using the interaction term of Equation (3.24), the mode I stress intensity factor may be expressed as l-v2 G (0 2A, (3.25) Similarly, by superimposing a pure mode case {Kiaux = 0, Kuaux = 1), the mode II stress intensity factor may evaluated as K, A2) rint Kl-v'JIA,, (3.26) Chapter 3. Dynamic Fracture Parameters 42 The interaction energy release rate may be evaluated by superimposing the auxiliary field solutions corresponding to a unit mode I stress intensity factor in Equation (3.16). The contour integral may then be converted into a domain form following the procedure of Section 3.5. For small strain linear elasticity, Shih and Asaro (1988) showed that = J[-('<7F F '4'> + p Vv<»)9, H'*? '«,,+'*<, '"?!)q.j 'A - (3.27) where y is linear strain and the superscript (1 ) refers to mode I auxiliary field quantity. The same expression can be obtained for G{m '(2) 'int ' 3.7 D I R E C T I O N OF C R A C K P R O P A G A T I O N Different theories are available in the literature to predict the direction of crack propagation. Erdogan and Sih (1963) proposed that cracks propagate in a radial direction corresponding to the plane perpendicular to the direction of the greatest tangential stress. Alternatively, Sih (1974) introduced a theory based on the strain energy density. It states that the crack propagation starts in a radial direction along which the strain energy density is a minimum. In this work, both theories will be considered in the numerical implementation. Chapter 4 NUMERICAL PROCEDURE This chapter discusses the implementation of the A L E and dynamic fracture formulations, described in Chapters 2 and 3, into the newly developed dynamic A L E finite element FRacture program (ALEFR). New procedures and their numerical implementation into the developed program are detailed in this chapter. 4.1 M E S H M O T I O N As indicated in previous chapters, the finite element grid points in A L E formulation may be moved arbitrarily to maintain a homogeneous mesh and to properly represent boundary conditions throughout the deformation process. Grid displacements are first related to material displacements through a set of assigned arbitrary mesh motion parameters. This simplifies the specification of pure Lagrangian, pure Eulerian or any arbitrary degrees of freedom. The choice of the arbitrary mesh motion parameters for interior degrees of freedom is handled by a special mesh motion scheme. However, special treatment for mesh motion on free material boundaries is necessary to satisfy the A L E boundary constraint, Equation (3.5). 43 Chapter 4. Numerical Procedure 44 4.1.1 Grid Displacement The grid displacement is related to the material displacement according to the following general form u s =a + Bu (4.1) where a and B are a vector and a matrix of mesh motion parameters, respectively. Vector a consists of appropriate grid displacements given by the mesh motion scheme, while matrix B consists of factors that allow the coupling of grid and material displacements. Setting a = 0 and B = I, produces a pure Lagrangian scheme, i.e., u* = u, whereas a pure Eulerian scheme, i.e., us = 0, may be obtained by setting a = 0 and B = 0 . Depending on the scheme of mesh motion on free boundaries, matrix B may have non-zero off-diagonal terms, which means that grid and material displacements are coupled at the same node. 4.1.2 Mesh Motion for Interior Nodes In this work, the transfinite mapping method (Haber et al., 1981) is used as the mesh motion scheme for the degrees of freedom interior to any mesh region bounded by any four specified boundary curves. This method provides a homogeneous mesh and matches the boundary of a given region at an infinite number of points. Another distinct advantage of the transfinite mapping method is that it allows the discrete representation of boundary curves; that is, the coordinates and displacements of boundary nodes may be Chapter 4. Numerical Procedure 45 used to find the optimum position of the nodes internal to the region. It also allows for discontinuities in the slope of boundary curves. 4.1.3 Mesh Motion on Free Material Boundaries To move nodes on free boundaries, three important aspects have to be considered. First, the true material boundary should be determined and tracked along the course of deformation. Second, the A L E boundary constraint, Equation (3.5), must always be consistently satisfied. Third, the boundary nodes should be properly spaced to ensure a good quality mesh within the region they bound. The description and analysis of two techniques for the treatment of boundary motion are given in this section. These are: the Displacement Coupling (DC) technique (Bayoumi and Gadala, 1999) and the True Boundary Tracking (TBT) technique (Abdelgalil and Gadala, 1999), developed in this work. The Displacement Coupling technique Bayoumi and Gadala (1999) controlled the A L E nodal motion on free boundaries, by using the general equation of grid motion, Equation (4.1). By considering point k on the free boundary, Figure (4.1), the application of Equation (4.1) yields Components of the vector a and the matrix B are determined such that the resultant coupling of degrees of freedom satisfies the A L E boundary constraint, Equation (3.5). (4.2) Chapter 4. Numerical Procedure 4b_ The free boundary is defined by a cubic spline interpolation of nodal coordinates. A local set of axes, x' and y' is defined at node k, such that x' is tangent to the boundary. This local set of axes makes an angle 6 with the global x and y axes. Components of the vector a and the matrix B are determined such that the resultant coupling of degrees of freedom satisfies the A L E boundary constraint, Equation (3.5). The free boundary is defined by a cubic spline interpolation of nodal coordinates. A local set of axes, x' and / is defined at node k, such that x' is tangent to the boundary. This local set of axes makes an angle #with the global x and y axes. y free material boundary Figure (4.1) The Displacement Coupling (DC) technique. Chapter 4. Numerical Procedure 47 An arbitrary displacement ax, is assigned to node k along the tangent axis x'. This displacement is set in order to control the nodal spacing on the boundary. Equation (4.1), may then be written as The above scheme satisfies the A L E boundary constraint requirements, but it does not truly track the location of the boundary node. This may be crucial, especially in the case of a moving crack with a curved or kinked propagation path. True Boundary Tracking (TBT) technique This technique is developed in this work, and it is based on the physical interpretation of the A L E boundary constraint. The objective is to accurately determine and track the true material boundaries, including the evolutional crack faces. Implementation of this technique into an A L E description is not straightforward, but it might be simplified, according to Figure (4.2), into the following steps: ax, cosd ax, sin 0 + sin 26 -sinr^cosf? \ux -sin6>cost9 cos2 0 [uy (4.3) The total nodal point displacement is resolved as; us =uL + uA (4.4) where u L is the Lagrangian component of the nodal displacement, that is, the displacement of a material point that instantaneously shares the location with the node, whereas uA is the arbitrary component of the nodal displacement. Chapter 4. Numerical Procedure 48 Initially, before load application, a number of Material Points (MPs) are created on each free boundary. The possession of these MPs to their corresponding elements is then determined. Due to the nature of the A L E motion scheme, the possession of MPs may change form one element to another; therefore, the search for the owner element is repeated each time step. L*-i (was MP£j) Ml* : Material point k possessed by element i Element i Material point ei : + : • : Finite element node Martial boundary Figure (4.2) The True Boundary Tracking (TBT) technique. Chapter 4. Numerical Procedure 49 After the application of the loading increment, each nodal point k is moved according to its respective Lagrangian displacement, uL. Then, for a material point j, possessed by the element q, M P J , the deformation is determined by the extrapolation ( n (4.5) where the subscripts i and k correspond to the degree of freedom and node number, respectively. The new locations of MP's constitute the new (true) material boundary. Nodal points may now be moved arbitrarily by a distance uA on the true material boundary in accordance with the A L E boundary constraint. The TBT motion scheme consistently satisfies the A L E boundary constraint in the sense that the arbitrary nodal motion occurs exactly on the true material boundary, whereas in the coupled displacement technique, this motion occurs in a virtual material boundary, created by curve fitting of the nodal point coordinates. The TBT scheme treats the material boundary and the finite element boundary as two separate entities. This gives it the advantage of moving boundary nodes on uneven or curved free boundary contours. Consider, as an example, the free body contour shown in Figure (4.3). For the sake of simplicity, assume that the body is unloaded, i.e., stress Chapter 4. Numerical Procedure 50 free. An arbitrary mesh motion should not, in principle, induce any stress or deformation on the body and its free boundaries. However, moving the boundary nodes according to coupled displacement scheme (in which, nodes are moved in a direction tangent to the cubic spline curve connecting the nodes) distorted the unstressed boundary. On the other hand, the way the TBT defines and tracks the material boundary, kept this boundary free from any distortion throughout the process of arbitrary nodal motion. Figure (4.3) Nodal motion on free boundary: the TBT versus the DC techniques To implement the TBT approach, it is more convenient to rewrite the general equation of relative motion, Equation (4.1), as u g = a* where, a* = \ uL Lagrangian 0 Eulerian Arbitrary value A L E (4.6) Chapter 4. Numerical Procedure 51 where the arbitrary value, in A L E description, is assigned by the transfinite mapping method for internal nodes, whereas the arbitrary motion of the free boundary nodes are assigned according to the TBT approach. 4.2 M A T E R I A L POINT SPLITTING (MPS) The most widely used technique in the modeling of crack propagation process in Lagrangian finite element codes, is the node release method (Anderson, 1973). In this method, the crack propagation is modeled by discontinuous jumps between fixed nodes in a straight line determined prior to the application of the finite element analysis. In each crack growth increment, the crack tip node is split into two nodes and the next node in line becomes the new crack tip node. Thus, the crack is forced to grow, in a self similar manner, a whole element length for each increment. These nonphysical assumptions introduce some errors, in the analysis of dynamic fracture in general and in modeling ductile material in particular, where the accuracy of determining the stress and strain histories affects the solution significantly. This drawback limited the application of node release approach to some particular classes of crack propagation problems. In this work, a new technique is developed to model the process of crack propagation in a continuous manner (Abdelgalil and Gadala, 1999). This method models the creation of new material surfaces (due to crack propagation) by splitting every material point that happens to lie on the crack propagation path into two material points, instead of splitting crack tip nodes. Therefore, this method is referred to as: Material Point Splitting (MPS). Chapter 4. Numerical Procedure 52 The MPS technique is incorporated with the TBT to enable the modeling of crack propagation using the A L E finite elements in a continuous and consistent manner. This implementation is simplified in Figure (4.4) for a self-similar crack propagation. According to this new approach, the crack tip is associated to a single finite element node throughout the course of crack propagation. The Crack Tip Node (CTN) has a material-point-generation function. If the crack is allowed to move, according to the propagation criteria, this function generates a new material point at the new crack tip location (MP C T ) . The old (MPCT) , directly behind the CTN, splits into two MPs. In the next time step, these split MPs are associated to the upper and lower crack faces. By splitting each MP, new upper and lower crack faces are created in each increment. The objective of the TBT technique is to control nodal motion on the upper and lower crack faces. This insures that the mesh is always kept uniform during crack growth, and it keeps the motion of boundary nodes consistent with the definition of the relative motion, with respect to the computational and material reference frames. Chapter 4. Numerical Procedure 53 At the end of time step t Lagrangian deformations are applied to all the nodes and MPs according to TBT. + M P Material surface 11 Node Element surface Arbitrary motion may be applied to all nodes except CTN. (No crack propagation yet) The CTN is released from the M P C T and allowed to advance leaving M P C T behind. The M P C T is then split into two MPs. At the end of time step f+1 the split MPs will move to the upper and lower crack faces according to their Lagrangian deformations. New M P C T N is created. Figure (4.4) Material Point Split (MPS) technique. Chapter 4. Numerical Procedure 54 4.3 NUMERICAL IMPLEMENTATIONS The implementation of the A L E formulations, Chapter 2, and dynamic fracture formulation and techniques, Chapters 3 and 4, into a dynamic A L E finite element fracture program (ALEFR) is outlined in this section. 4.3.1 Solution of Equilibrium Equations In order to solve the finite element equilibrium equations, Equations (2.64) and (2.73), using a modified Newton-Raphson iterative scheme, the two equations may respectively be rewritten as 'M l a l ' ) +'C A l v ( 0 + , C X 2 (v ( i ) -Y J , i ) )+('K l l + , K t 2 )u ( i ) +('KM+'KA2)(u(i) -u s ( 0 ) where u ( , ) , v ( , ) and a ( , ) are the corrections to the incremental material displacement, velocity and acceleration vectors in iteration i, respectively, and u s ( , ) and vA'0) are the corrections to the incremental grid displacement and velocity vectors, respectively. These correction vectors are then used to update their corresponding quantity as follows ('K^'+'K 1 2)^ +('K>ll+'KAJ)(u(') - u'(,))='+Af rt-'f(w) (4.7) and t+At (i) _t+At + U (<•) + V (/) a ( 0 = r + * a « - l ) + a ( 0 (4.9) t+At g(i) _t+At g(i-\) + u t + At yg(0 _t + Atyg(i-]) + V Chapter 4. Numerical Procedure 55 In dynamic analysis, time integration is carried out using Newmark implicit scheme (Bathe, 1996) in which ' + A ' u ( 0 ='u + At' v + At2 (| - ft)'* + At2p,+Ata(i) (4.10) '+A< v ( 0 =' v + At(\ - y)' a + Aty' + A ' a ( i ) (4.11) where (3 and ^  are parameters that control the accuracy and stability of integration. Using Equations (4.9), (4.10) and (4.11), v 0 ) , a< 0 and \ 8 ( , ) may be eliminated from Equation (4.8) by correlating them to u ( , ) and us(,) as follows v(o = _JLu<o Atfi • a ( 0 = — ^ - u ( 0 (4.12) At2fi Atfi Using Equation (4.12), Equation (4.8) may be rearranged as [-L-'ML + -^('CA]+'CA2)+('KL+'KA)]u(i)-[-^'CA2+'KA]u "' Atlfi Atfi Atfi (4.13) _t+At ^ext _t+Ai j(i-i) _(']yji'-)-']yj/'y+z1' a^ '"1' —(' c*' +'C'A^+'C'**)'+^1' v^ '""1' Equations (4.7) and (4.13) may be simplified further by eliminating u i , ( , ) . This may be carried out by rewriting both of the equations in the following general form ' K u ( 0 - ' K s u i , ( 0 =f ( ; ) (4.14) Chapter 4. Numerical Procedure 56 where ' K and ' K x are equivalent stiffness matrices corresponding to u ( , ) and u s ( , ) , respectively, while f ( , ) is the incremental load vector for iteration i. Then, using the A L E correlation between grid and material point displacement according to TBT, Equation (4.6), the general form, Equation (4.14) becomes ( ' K ) u ( / ) = f ( ' V K V (4.15) Equation (4.15) is used to eliminate the grid displacement on the element level. Conventional finite element assembly and elimination techniques may now be applied directly to solve for the unknown material displacements. 4.3.2 Numerical Integration of the Energy Release Rate In order to numerically evaluate the energy release rate G, and subsequently the dynamic stress intensity factor, the integrals derived in Sections 3.5 and 3.6 are implemented in the developed finite element program. These G integrals employ the stress work density, 'U, which, according to its definition in Equation (3.4), requires the integration of the product of stress and the velocity gradient, from the beginning until the current time t. It is not possible to employ the Gauss quadrature integration scheme, used in the standard finite element, because when the finite element mesh moves arbitrarily, due to crack propagation, the location of the Gauss points change. Therefore, the material point stress and velocity gradient histories are unknown, and the integration, Equation (3.4), cannot be performed. Chapter 4. Numerical Procedure 57 This problem is overcome by creating a Lagrangian grid of material points in the area around the crack tip where the G contours are applied. The grid should be created initially before load application. Each material point is possessed by an (owner) element, i.e., located within its boundaries. Due to the relative (ALE) motion between the finite element mesh and the Lagrangian grid in the event of crack propagation, the search for the owner element for each material point has to take place each time step. To accelerate the searching process, the search will first take place in the last owner element, then it is performed in the elements adjacent to it. A map of the neighboring elements for each element in the model, is determined by utilizing element connectivity at the beginning of the analysis. This map will not change throughout the course of crack propagation, because there is no remeshing or nodal release used in modeling of crack propagation. Once the owner element is located, material point properties are then updated. A square contour centered at the crack tip is constructed, as shown in the Figure (4.5), such that it starts at the lower crack face, passes through the MPs of the Lagrangian grid and ends at the upper crack face. The material points bounded by the contour, including the ones on the contour itself, are referred to as domain MPs, whereas the MPs on the contour are referred to as contour MPS. Each term of the contour and domain integrands of G integrals, is evaluated at each contour and domain MP, respectively. The one and two dimensional Trapezoidal rules are then applied to carry out the numerical integration. Chapter 4. Numerical Procedure 58 Once the crack is allowed to propagate at a certain speed, V , the crack tip node will then displace a distance d = V'At, where At is the current time step. Accordingly, the TBT, MPS and the transfinite mapping scheme arrange the elements topography, Figure (4.5). Holding the square contour while the crack is propagating causes contour eccentricity, i.e., the crack tip is shifted from the center of the square. Excessive eccentricity leads to inaccurate results (Organ, 1996), and eventually, causes the domain to disintegrate into two separate parts. Therefore, continuous contour reconstruction is needed as long as the crack is propagating. Figure (4.5) shows, as an example, a self-similar crack propagation, in which the contour was reconstructed by eliminating the contour MPs on the left and incorporating new contour MPs on the right. This process of contour reconstruction involves the search for the best combination of MPs to make the contour as square and eccentric as possible. Practically, the energy release rate is evaluated using more than one contour, and in this work, a large number of square and rectangular contours are used in each example to determine the sensitivity of these contours to the size and shape of the path of integration. Chapter 4. Numerical Procedure 59 4 ft f-node GDoma integral + + + -ft V 7* --+ f +• d = VAt + -+-• + + •4 +• t +-4 + i+ +-• if +-•4---4 -4-+ + 4-,-,-H + + Lagrangian grid Material point Finite element node Finite element mesh G contour integral Contour and domain integrals move with the crack tip Figure (4.5) Contour and domain integrals to determine G for dynamic crack propagation Chapter 4. Numerical Procedure 60 4.4 T H E A L E D Y N A M I C F R A C T U R E P R O G R A M The techniques developed in this work to model dynamic fracture problems, are implemented in a program (ALEFR). These techniques include: the Material Point Splitting (MPS), True Boundary Tracking (TBT) and techniques to determine the fracture parameters for stationary and propagating cracks. 4.4.1 Program Structure The developed program (ALEFR) is written in a modular scheme and it employs standard updated Lagrangian procedures, routines for general A L E analysis, and the A L E routines developed specifically to model fracture problems. The flowchart, Figure (4.6), shows the main routines of the developed program. Additional standard and developed subroutines are not listed here for brevity. The program starts by reading the variables needed for controlling the size of the arrays through the dynamic memory allocation. Next, all variables and arrays used in the program are initialized. The program then reads the standard finite element input data as well as fracture control data. Since the program is designed to handle the conditions of displacement controlled and load (stress) controlled fracture, the applied distributed loads are converted into their equivalent nodal values. A grid of material points is created within the domain of the body. A searching process is carried out within the whole domain, for the owner element of each MP. This grid is needed later in the computation of fracture parameters and in the TBT scheme. Chapter 4. Numerical Procedure 61 In each Newton equilibrium iteration, the fracture parameters have to be evaluated. This is carried out first by searching for the current owner element for each M P in the grid created at the beginning. Once all the owner elements are located, the M P properties are updated accordingly. The integration contours and domains are built on the M P grid and the strain energy release rate and the stress intensity factors are determined. According to the crack propagation criterion used, if a decision is made to propagate the crack, propagation speed and direction are subsequently determined and the MPS technique is used to create the new crack surface. The TBT and motion routines are called to control the A L E mesh motion by determining the motion parameter in the vector a* . For the free boundary nodes, including crack faces, this parameter is computed using the TBT routines, whereas, for internal nodes, it is determined by calling the transfinite mapping subroutines. In the condition of a load controlled fracture, the A L E motion due to crack propagation may alter the position of the nodes included in the calculation of the consistent load vector. This vector is re-evaluated in each iteration while the crack is propagating. When the boundary conditions are set, the standard predictor phase of the time integration algorithm is then preformed, and the elements lumped mass and stiffness matrices, including the extra terms due to A L E , are calculated. The global effective load vector is assembled from all the elements load vectors as well as the out of balance internal forces from previous iterations. The solution for the unknown displacement is carried out through a frontal solution scheme. A line search algorithm is used to improve convergence. The corrector phase of the time integration Chapter 4. Numerical Procedure 62 scheme is then performed to determine velocities and accelerations. Gauss point stresses are integrated and three convergence measures are determined. These measures are respectively based on: the iterative displacements, the residual (out of balance) forces and the residual energy. Convergence of iterations is then checked and the iterations are performed until convergence is achieved. At the end of each converged time step, nodal stresses are extrapolated from Gauss point stresses, and the results are written into output files. The program then moves to the next time step. Chapter 4. Numerical Procedure 63 ( START 7 ") Read control variables / CQNTRO / Initialize arrays and variables | INITIA | Input model and fracture data / INPUTD J Evaluate consistent nodal force vector Create grid of material points Determine grid point possession Update grid points properties Define contour and domain integrals Determine fracure parameters Determine propagation speed and direction Perform MPS Perform TBT for crack faces and free boundaries Set motion parameters Update consistent load vectors due to crack propagation Set boundary conditions Predictor phase Determine elements mass matrices including ALE terms Calculate elements stiffenss matrices including ALE terms Assemble elements effictive load vector including ALE terms Solve ALE equilibrim equations Line search algorithm Corrector phase Evaluate stresses at Gauss points Determine convergence parameters Extrapolate nodal stresses Write output data / OUTPUT / Figure (4.6) Flowchart of the developed program (ALEFR) Chapter 4. Numerical Procedure 64 4.4.2 Characteristics of the Developed Program (ALEFR) The characteristics of the developed implementation, as discussed above are highlighted in the following points: - Due to the developed true boundary tracking (TBT) technique, the motion of the free boundary nodes is accurately tracked, and the A L E free boundary constraint equation is satisfied. - The modeling of the crack propagation process using the material point splitting (MPS) technique (as opposed to the Lagrangian nodal release method), provides a closer representation of the real physical process, which involves the evolution of new material surfaces in a continuous manner. In addition, the MPS technique eliminates the need for dense mesh in the path of the propagating crack tip, as traditionally required in the nodal release method. - Unlike the nodal release method, the developed program is capable of modeling in-plane mixed mode crack propagation without the need to determine the direction of crack propagation prior to the start of analysis. Further, the developed program may easily model non self-similar crack propagation, i.e., curved or kinked cracks. - In the developed program, the choice of the time step At is independent of the crack tip velocity and element size. Variable crack tip velocities may, therefore, be modeled more naturally. Chapter 4. Numerical Procedure 65 - The strain energy release rate may be evaluated through integration over contours and domains of inhomogeneous materials, i.e., a single contour may pass though different materials. - The TBT and the MPS techniques are carried out within the iteration loop, i.e., before achieving convergence. This is due to the fact that the stresses calculated within the iterations already include convective effects. - A l l output data is written in formats, such that they can be readily read and animated by the data visualization software TECPLOT®. Chapter 5 NUMERICAL EXAMPLES In this chapter, the developed program (ALEFR) is tested in a series of numerical experiments that include stationary and propagating cracks. Test results are then compared to the benchmark solutions available in the literature. Although the developed program is capable of modeling linear and nonlinear material behavior, in all cases considered here, only linear elastic material is assumed. In addition, a mixed mode dynamic crack propagation problem is modeled to test the ability of the developed code in solving these kinds of problems. 5.1 STATIONARY MODE I CRACK SUBJECTED TO SINGLE STEP PULSE 5.1.1 Problem Description In this example, a semi-infinite crack in an infinite body is loaded with a tensile stress wave traveling perpendicular to the crack surface. Freund (1990) presented an analytical solution for the dynamic stress intensity factor, which starts to develop as soon as the wave reaches the crack tip; according to (5.1) 66 Chapter 5. Numerical Examples 67 where, v' is the uniform crack propagation speed ( v c = 0 in this example), v is Poisson's ratio, cr0 is the applied tensile stress pulse and cd is the dilatational wave speed. In the above equation, the measurement of time t started as soon as the wave reached the crack tip. Equation (5.1) shows that the stress intensity factor is a function of the square root of time t, and is valid as long as the reflected boundary waves do not arrive at the crack tip. 5.1.2 Numerical Procedure (TO H Y j ,-. T t Y Y ! A A A A -+X Figure (5.1) Horizontal edge crack on a thick plate used in mode I fracture cases. The plate shown in Figure (5.1), has a horizontal edge crack. The material properties are: £ = 211 GPa, and p = 7800 K g / m 3 . A step tensile stress plus a0 = 1000 N / m 2 , is applied on the upper surface of the plate creating a stress wave which arrives at the crack tip at time t = H/cd . Plate dimensions are H = 2 m and L = 10 m, and a horizontal edge crack Chapter 5. Numerical Examples 68 of length ao = 5 m are chosen such that the reflected boundary waves do not arrive at the crack tip within the time span reported in the analytical solution (about 1.0 ms). The plate is modeled by a coarse and fine model. In the coarse discretization, 1000 identical 4-noded quadrilateral plane strain elements in (50 x 20) arrangement, are used, whereas in the fine discretization, 4000 of the above element type are used in (100 x 40) arrangement. Different time steps ranging from 0.3 Atc to 1.0 Atc, where Atc is the critical time step according to the Courant condition (Isaacson and Keller, 1966) given by Atc = l/cd, where / refers to element length. Three Lagrangian grids of different densities have been used to mesh a square area (4 m x 4 m) bounded by the largest energy contour integral and centered at the crack tip. In these grids, material points are arranged in (101 x 41), (401 x 161) and (801 x 321) uniform meshes. To study the path dependency of the energy release rate integrals, 200 integration paths, centered at the crack tip, have been created, Figure (5.1). The shortest path is just 0.02 m away from the crack tip, while the longest path is 5 m away. These integration contours pass through material points of the Lagrangian grid. 5.1.3 Numerical Results Figure (5.2) shows snap shots of the ay stress distribution in the model at different time steps. The step pulse is applied at time t = 0, at which stress waves start to emanate from the upper plate surface and move downward, such that they reach the crack tip at Chapter 5. Numerical Examples 69 t-H/cd . The solution is carried out until time t = 3H/cd, the time at which the reflected waves from the upper and lower surfaces reach the crack tip. The snap shots also captured the reflected waves from the right and left free ends. These waves do not arrive at the crack tip within the specified solution period. Different combinations of discretization schemes, time steps and Lagrangian grid densities are utilized in this example. The finer discretization, denser Lagrangian grid and smaller time steps always give slightly smoother results; however, all of the results obtained are in very good agreement with the theoretical solution (Freund, 1973; Freund, 1990). The energy release rate is evaluated from the dynamic stress intensity factor using Equation (3.21). The change of the energy release rate (normalized byHa\ IE) with time (normalized by Hlcd) is shown in Figure (5.3). The above solution was obtained through the integration of Equation (3.19) over a path located 3 m away from the crack tip employing the finest model and the densest Lagrangian grid. A unit value of normalized time is equal to the time required for the wave emanating from the upper surface to reach the crack tip. Although the waves pass through all the Lagrangian grid points on the upper half before they reach the crack tip, all the terms of the domain and contours integrals of the energy release rate, Equations (3.19) and (3.20), vanish due to the instantaneous symmetry about the crack tip y axis during that period of time. Chapter 5. Numerical Examples 70 Figure (5.2) Snap shots of ay stress distribution due to the applied tensile stress pulse in the stationary crack model. Chapter 5. Numerical Examples 71 Time = 2 0 t = 0.663 ms 500 Figure (5.2) (continued) Snap shots of ay stress distribution due to the applied tensile stress pulse in the stationary crack model. -1000 Chapter 5. Numerical Examples 72 2.5 2 0.5 0 • G I j i L„.vx i i ,J uniiiiii mull &r i i i i i i i Figure (5.3) Variation of the normalized strain energy release rate with the normalized time using the fine model with (801 x 321) Lagrangian grid in the stationary crack model, subjected to a tensile stress pulse. After a unit normalized time t, the energy release rate starts to rise in an almost linear manner in very good agreement with the analytical solution until t < 2.6. After this time, oscillation of the results increases and it starts to deviate from the analytical solution. This may be attributed to the influence of the reflected boundary waves. Finally, the Chapter 5. Numerical Examples 73 dynamic stress intensity factor is extracted from the energy release rate using Equation (3.21), and its variation with the normalized time is shown in Figure (5.4). 2500 2000 1500 1000 500 K :a\) K ftheoretu - i - i I ; ! L i : : I 7° : : : : / " i i : ! 1° I ! i i ( j, i i i i i i i i ) 1 2 3 Figure (5.4) variation of the dynamic stress intensity factor with the normalized time using the fine model with (801 x 321) Lagrangian grid. 5.1.4 Analysis of Contour Path Dependency The integrals of Equations (3.18), (3.19) and (3.20) possess the fundamental characteristic that they should evaluate the energy release rate regardless of the size and Chapter 5. Numerical Examples 74 shape of the domain of integration used, i.e., they are path independent. To examine the path dependency of these integrals, as many as 200 paths of integration have been constructed. Each path starts from the lower crack face and ends at the upper one. Within 2 m of the crack tip, all the paths are square in shape, whereas larger contours are rectangular in shape. Domain and contour integration are carried out over the material points of the Lagrangian grid. Figure (5.5) shows the path dependency for each of the above integrals at different time steps starting at the normalized time t = 1.0, the time at which the stress wave arrives at the crack tip. Since, for straight cracks, Equation (3.18) contains domain integrals only, the variation of G values was very smooth and perfectly path independent for the contours further than 0.5 m from the crack tip. Equations (3.19) and (3.20) have two common interesting features: first, for all contour (and domain) integrals and at all the time steps, they both produce, up to the third decimal point, identical G values (therefore, the two equations appear as one curve in Figures (5.5)). Second, in each of these equations, after decomposing into contour and domain integrals, the domain integral is almost a mirror image of its contour integral, i.e., the contour and domain integrals complement each other to produce an almost path independent G. In both of these equations, the contour integral dominates in the area closer to the crack tip, while in Equation (3.19), the domain integral dominates afterwards. Chapter 5. Numerical Examples 75 time 1.0 2 <5 1 0 t = 0.331 ms G [Eq.3.18] G c [ E q . 3.19] G d [ E q . 3.19] G c [Eq. 3.20] G d [Eq. 3.20] G [Eq. 3.19 & 3.20] G (theoretical) I 0 1 2 3 4 5 Distance from crack tip [m] time 1.5 t = 0.497 ms| G [Eq.3.18] G c[Eq. 3.19] G d [Eq. 3.19] G c [Eq. 3.20] G d [Eq. 3.20] G [Eq. 3.19 & 3.20] G (theoretical) 2 O 1 0 --/ . • • . ? • / ;< — -i i 0 1 2 3 4 5 Distance from crack tip [m] Figure (5.5) Path dependency of the contour and domain integrals corresponding to the normalized energy release rate at different time steps using (401 x 161) Lagrangian grid. Chapter 5. Numerical Examples 76 time 2.01 = 0.663 msl (3 1 G [Eq. 3.18] G c [Eq. 3.19] G d [Eq. 3.19] G c [Eq. 3.20] G d [Eq. 3.20] G [Eq. 3.19 & 3.20] G (theoretical) ~ ~ t - T If ,-. A . A A . A . A A ' U U . 2 3 Distance from crack tip [m] time 2.5 t = 0.829 ms O 1 7 ^ i~ G [Eq. 3.18] G c [Eq. 3.19] G d [ E q . 3.19] G c [Eq. 3.20] G d [Eq. 3.20] G [Eq. 3.19 & 3.20] — G (theoretical) 2 3 Distance from crack tip [m] Figure (5.5) (continued) Path dependency of the contour and domain integrals corresponding to the normalized energy release rate at different time steps using (401 x 161) Lagrangian grid. Chapter 5. Numerical Examples 77 time 3.0 2 (3 1 0 t = 0.829 ms " . ..4 ^ ^ i S J\ A 1 J ^ W 1 r7! 1 *****/ 1 •**"* i. / / > A 7 . ,i ?:..!— ~ G [Eq.3.18] G c [ E q . 3.19] G d [ E q . 3.19] - G c [ E q . 3.20] G d [ E q . 3.20] G [Eq. 3.19 & 3.20] G (theoretical) -i i • / / '< / / — —. - .< " ': i 1 1 0 1 2 3 4 5 Distance from crack tip [m] Figure (5.5) (continued) Path dependency of the contour and domain integrals corresponding to the normalized energy release rate at different time steps using (401 x 161) Lagrangian grid. The variation of the normalized energy release rate, evaluated by integrating Equations (3.18), (3.19) and (3.20) over a contour located 3 m away from the crack tip, versus the normalized time, is shown in Figure (5.6). It may be noted that variations in domain integrals are generally smoother than those of contour integrals for the normalized time period t > 2.0. The increase in G, due to the domain integral part of Equation (3.19) for t > 2.2, was balanced by the fluctuating contour integral. The maximum error produced by all of the three equations is about % 11 at the later stages of modeling. Chapter 5. Numerical Examples 78 G [Eq. 3.18] G c[Eq. 3.19] Gd[Eq . 3.19] G c [Eq. 3.20] G d [Eq. 3.20] G [Eq. 3.19, 3.20] G (theoretical) Normalized time Figure (5.6) Variation of the contour and domain integrals corresponding to the normalized energy release rate with the normalized time, using (401 x 161) Lagrangian grid. Chapter 5. Numerical Examples 79 5.2 STATIONARY MODE I CRACK SUBJECTED TO SYMMETRIC STEP PULSE 5.2.1 Problem Description The plate shown in Figure (5.1) is simultaneously loaded by a step tensile pulse ao at both the upper and lower surfaces. Equation (5.1) is still applicable to this problem, and by superposition, the dynamic stress intensity factor may be expressed as 5.2.2 Numerical Procedure Since this problem possesses symmetry with respect to the crack plane only the lower half of the plate is considered for modeling. The model is discretized into 500 and 1000 4-noded identical square elements. Displacement boundary conditions are prescribed at the corresponding nodes on the axis of symmetry. Uniform time steps ranging from At -0.3 Atc to 1.0 Atc are used. Lagrangian grids of (41 x 21), (161 x 81) and (321 x 161) material point arrangement are employed. 5.2.3 Numerical Results Snap shots of av (= 'cr22) distribution at different time steps are shown in Figure (5.7). Variation of the dynamic stress intensity factor, normalized by (CT0^H/(1-V2) ), with the normalized time, is shown in Figure (5.8). It may be seen that the numerical results are in very good agreement with the theoretical solution throughout the allotted time interval. (5.2) Chapter 5. Numerical Examples 80 Figure (5.7) Snap shots of oy stress distribution due to the applied symmetric tensile stress pulse, in the stationary crack model. Chapter 5. Numerical Examples 81 Figure (5.7) (continued) Snap shots of ay stress distribution due to the applied symmetric tensile stress pulse, in the stationary crack model. Chapter 5. Numerical Examples 82 5.3 PROPAGATING MODE I CRACK SUBJECTED TO SINGLE STEP PULSE 5.3.1 Problem Description This problem is similar to that described in Section (5.1). However, in this problem the crack will , at a certain point in time, propagate at a prescribed velocity. Freund (1990) Chapter 5. Numerical Examples 83 obtained an analytical solution to determine the energy release rate for a horizontal edge crack in an infinite plate propagating in a self-similar fashion (mode I). For a propagating crack, Equation (5.1) may be written as (Freund, 1990) K,(ty) = k(vc)K,(t,0) = fc(v')gg°JC'(1 2 v ) t (5.3) l - v V n where, k(vc) is a universal function of crack tip speed, is equal to unity for stationary crack and decreases monotonically to zero when the crack speed approaches the Raleigh wave speed, cR, according to this approximation k(vc)= ) ~ v C ' C r (5.4) The energy release rate may then be related to the dynamic stress intensity factor using Equation (3.21). From Equations (3.21) and (5.1) to (5.4), it may be concluded that the energy release rate changes linearly with time at a slope of g{vc)<rl C, where C is a contact based on elastic material properties only and g(vc)is another universal function of crack tip velocity and is expressed as g(vc)=AI(vc)k2(vc) = l-vc/cR (5.5) Thus, the slope is a maximum for stationary cracks and it decreases with the increase of propagation speed v c . Equation (5.1) is applicable to mode I crack propagation provided that the reflected waves from boundaries do not reach the vicinity of crack tip. Chapter 5. Numerical Examples 84 5.3.2 Numerical Procedure The plate shown in Figure (5.1) is subjected to a single step pulse. The discretizations and time steps used in the stationary crack example (Section 5.1) is used to model this problem, that is., 1000 and 4000 elements and 0.3 Atc to 1.0 Atc time steps. Since the A L E motion due to crack propagation process causes a change to the dimensions of some elements of the model, the selection of time step size is made considering the shortest element length. Although the code is capable of changing the time step in an adaptive manner, only uniform time steps have been used. To capture the change in the energy release rate of the propagating crack, the whole model area is covered by the Lagrangian grid of the uniform densities of (101 x41), (401 x 161) and (801 x 321) material point arrangements. Integration of the energy release rate is carried out using 200 square and rectangular contour (and domain) integrals. The largest contour covers the entire model domain. To simplify the comparison of the numerical results with the analytical solution, the onset of crack propagation and the propagation velocity are fed to the code as prescribed values instead of computing them through a propagation criterion. The crack is, therefore, assumed to remain stationary until the normalized time t = 1.5, at which it is allowed to propagate at a uniform velocity v = 0.4 c, (~ 1290 m/s). Chapter 5. Numerical Examples 85 5.3.3 Numerical Results Animation of stress distribution during the crack propagation process shows that the process has been modeled in a very smooth manner. Snap shots from the animation of ay stress distribution are shown in Figure (5.9). The instantaneous loading and unloading processes, which take place in the front of, and behind, the propagating crack tip node respectively, are consistent with the motion of the crack tip node. This might be attributed particularly to the success of the MPS technique, which models the evolution of new surfaces due to crack propagation. The nodes on the free boundaries, including the original and the new crack faces, moved naturally on the true material boundary tracked by the TBT technique. Once the motion of the boundary nodes has been decided, the transfinite mapping technique continuously keeps the internal nodes in a uniform mesh. The integration of the energy release rate is carried out over the 200 contour and domain integrals composed by the material points of the Lagrangian grid. As far as the crack remains stationary, the constituency of the contour and domain integrals is kept unchanged. As soon as the crack begins to propagate, the integrity and centricity of these contours are checked every time step and, when needed, reconstruction of the contours is performed. Material point properties were updated by mapping them from the specific finite element that instantaneously 'owns' the material point. Ownership of the material points of the Lagrangian grid is tracked throughout the process of crack propagation. Chapter 5. Numerical Examples 86 The dynamic stress intensity factor, K, is extracted from the energy release rate in a way similar to that of the stationary crack shown in Section 5.1. The change of the dynamic stress intensity factor (normalized by CJ0^H/([-V2) ) with the normalized time is shown in Figure (5.10). For a normalized time t < 1.5, the period at which the crack remained stationary, the numerical and analytical solutions are in good agreement. The dotted curve in the figure represents the hypothetical value of the normalized K if the crack remained stationary throughout the solution period. The sudden drop of the normalized K at the normalized time t = 1.5, indicates that the start of the crack propagation process was captured by the numerical procedure. However, this drop is not as large as anticipated by the analytical solution. For the normalized time period 1.5 < t < 2.14, numerical results are smooth, however, they are about %32 higher than the analytically anticipated solution. For t > 2.14, numerical solutions show more fluctuation but the same average error is maintained. The solution very slowly converges to the analytical solution with increasing the density of the Lagrangian grid. Comparable results may not be reached without the use of an impractical, extremely dense Lagrangian grid. The reasons behind the inaccuracy of the numerical results are not known yet; however, they may be attributed to the accumulation of errors in the numerical techniques used to map material properties of the material points of the Lagrangian grid and the numerical techniques used to evaluate the energy domain and contour integrals. More investigation and future work is needed to clarify this point. Chapter 5. Numerical Examples 87 time = 1.5 t. ~ 0.497 ms Crack starts to propagate Figure (5.9) Snap shots of ay stress distribution in the dynamic stationary/ propagating crack model subjected to a tensile stress pulse. Chapter 5. Numerical Examples 88 1.5 -0.5 — K (numerical) K (theoretical) K stationary (theoretical) ^ -! A £.i .— / \ ^, Crack prop agating Stations xy crack ax i i : (UltVS i . i i i i Figure (5.10) Variation of the normalized dynamic stress intensity factor with the normalized time for the dynamic stationary/propagating crack model, subjected to a tensile stress pulse. Chapter 5. Numerical Examples 89 5.4 PROPAGATING MODE I CRACK SUBJECTED TO SYMMETRIC STEP PULSE 5.4.1 Problem Description This problem is similar to problem (5.3), however, this time the pulse is simultaneously applied to both the upper and lower plate faces. By applying the superposition principle, Equation (5.3) may be written as 5.4.2 Numerical Procedure Due to symmetry, only the lower part of the plate is modeled, and the necessary symmetric displacement boundary conditions are applied. Two models containing 500 and 1000 elements with a time step varying from 0.3 Atc to 1.0 Atc, are used. Lagrangian grids with (101 x 21), (401 x 81) and (801 x 161) material point arrangements covering the entire domain of the model are employed. In this problem, the crack is assumed to remain stationary until the normalized time t = 1.5, at which it is allowed to propagate at a uniform prescribed velocity v = 0.4 cs (~ 1290 m/s). 5.4.3 Numerical Results Snap shots of the animation of ay stress distribution are shown in Figure (5.11). In this example as well, the MPS technique produced a smooth crack propagation modeling. The free boundary nodes were moved continuously and consistently on the true material boundary during the crack propagation process. (5.6) Chapter 5. Numerical Examples 90 The integration of the energy release rate is carried out over the 200 domains and the dynamic stress intensity factor is subsequently extracted and plotted as a normalized value against the normalized time, Figure (5.12). The trend of the variation of the normalized K is similar to that of the previous example. The numerical results are in good agreement with the analytical solution as far as the crack remains stationary. The onset of crack propagation is evident by the sudden drop of K. Once the crack starts to propagate, the same deviation pattern as reported in Section (5.3) is shown. The same conclusion reached for the previous example is applied here as well. 5.5 G E N E R A L C O M M E N T S A B O U T M O D E I F R A C T U R E P R O B L E M S In modeling of the above dynamic stationary and propagating mode I crack problems, the following conclusions may be drawn. The modeling process may be divided into two parts: the simulation part and the dynamic stress intensity factor devaluation part. The simulation part is successful for both the stationary and the propagating problems. Stress waves are generated, propagated and reflected in the expected smooth manner. Stress field loading and unloading in front of and behind the propagating crack tip, respectively, takes place smoothly without discontinuities. The K results obtained for dynamically loaded stationary crack problems are in a good agreement with the available analytical solutions. In addition, the onset of crack propagation is accurately defined by a sudden drop in K. However, values of K obtained Chapter 5. Numerical Examples 91 for dynamic crack propagation differ from those expected analytically for the two cases considered. Very slow convergence of the numerical solution towards the analytical one is achieved by increasing the Lagrangian grid density and decreasing element sizes and time step. Further investigation is recommended to improve the convergence. In testing the path dependency, results obtained by all of the three formulae employed, Equations (3.18), (3.19) and (3.20), are nearly path independent for both stationary and propagation fracture cases. Sine Equation (3.20) does not include domain integrals; it produces more stable results compared to the other two equations, for both stationary and propagating crack cases. However, in propagation fracture models, for a normalized time t > 2.5, K values become path dependent. This may be associated to the dispersion of waves reflected from the boundary and to the numerical errors associated with the techniques used to evaluate K. The results obtained by the above three Equations are generally very close. Particularly, Equation (3.19) and (3.20) are computationally identical up to the third digit (of the normalized value). By examining all the components of the above two equations, we find that the contribution of the first integral in Equation (3.20) is negligible. We also noted that, for Equation (3.19), the second term of the contour integral, which include the kinetic energy density 'T, and the first term in the domain integral, are typically comparable in magnitude but are of opposite signs. Therefore, their net contribution is Chapter 5. Numerical Examples 92 negligible, and the above two equations may be simplified by the flowing form G = J['/Jn, - 'cr y m. \ , ] 'dC + ftp 'a. 'w., ] 'JA (5.7) ra 'A 5.6 M I X E D M O D E F R A C T U R E E X A M P L E S In order to explore the capabilities of the procedures and program in modeling the mixed mode crack propagation problem, we consider the modeling of a mixed mode problem in the following examples. 5.6.1 Curved Crack Growth A plate made of a ductile material with a horizontal edge crack is subjected to vertical tensile static loading. Typically, the crack is expected to grow in a self-similar (straight) manner for a relatively short distance before unstable crack growth and fracture occurs. In this test, however, the crack is not allowed to take a self-similar path; rather, it is forced to take an arbitrary curved path, as shown by the red curves in Figure (5.13). The dots in the figure constitute the updated true material boundaries tracked by the TBT technique. With each crack tip advance, a couple of material points are created by the MPS technique and subsequently tracked by the TBT technique. Chapter 5. Numerical Examples 93 time- 0.5 I <> lt>r> ins 100011 tim&~ 1.0t~ 0.331ms 1000 3500 3000 01* 0.994 m s Figure (5.11) Snap shots of ay stress distribution in the dynamic stationary/propagating crack model subjected to a symmetric tensile stress pulse. • -500 1 Crack propagated 0.645 m I B in 0.497 ms . Chapter 5. Numerical Examples 94 2.5 2 H 1.5 * 1 0.5 -0.5 K (numerical) — K fthenretirah < stationary (theoretical) a " 0 „*w 1 ::::::::::::2L:::::Z —* x •JT—\ ^tr. i \ i A "jf X /a i ^ :Jff 1 "jf r • y / . . ! / / ! ""/' i . . / ; fc::::!:::::: I IE... Crack propagating • Stationary crack at 1290 m/s : : : : : : [ : : : : r : : : : i : : : : : j 1 j i 1 i 0.5 1 1. t 2.5 Figure (5.12) Variation of the normalized dynamic stress intensity factor with the normalized time a propagating crack subjected to a symmetric step pulse. Chapter 5. Numerical Examples 95 Figure (5.13) Snap shots of ay stress distribution in the curved crack propagation model. True crack faces are shown in red dots. Chapter 5. Numerical Examples 96 By examining the motion of the nodes on the upper and lower crack faces, it is clear that throughout the course of crack propagation, these nodes always move on the updated true material boundaries, despite the sharp edges, at the original location of the crack tip, and the tight curved path created due to crack propagation. Hence, the A L E boundary constraint is truly satisfied through the employment of the TBT technique. 5.6.2 Dynamic Mixed Mode Crack Propagation: Problem Description Kalthoff and Winkler (1987) experimentally tested the problem of mixed mode crack propagation in a double notched plate subjected to side impact loading, Figure (5.14). Since the current code is capable of extracting only mode I stress intensity factor, Ki , from the energy release rate, G, modeling of the above problem is limited to the simulation of the crack propagation scenario. Due to the complexity of this type of problem, the available analytical solution, developed by Lee and Freund (1990), is only valid as long as the crack remains stationary. Figure (5.14) The double notched plate subjected to impact loading. Chapter 5. Numerical Examples 97 In the above experiment, a free plate is subjected to an impact loading between the two edge notches. The plate has the following dimensions: W = 4 m, H = 6 m, h = 3 m and the initial crack size a0 = 1 m. The plate material properties are: p = 7833 kg/m 3, E = 200 GPa and v = 0.25. The prescribed impact velocity vo = 16.5 m/s. 5.6.3 Numerical Results Due to symmetry, only the upper half of the plate is modeled. The model is discretized into 2400 identical linear quadrilateral elements, in (40 x 60) arrangement. The time step used is At = 0.56 Atc = 1.0 x 10 5 s, where Atc is the critical time step according to the Courant condition for finite elements. The snap shots, Figure (5.15), show the von Mises stress distribution at different time steps. The applied wave reached the crack tip at t = 0.18 ms and the crack is assumed to remain stationary until the time t = 0.25 ms, it then propagates at an angle of 70° with a velocity vc = 480 m/s (~ 0.15 c,). As may be noted from the figure, the crack propagation process took place in a smooth manner indicating that the TBT technique is working properly. The program did, however, terminate and stop at a time t = 2.5 ms due to element distortion. Careful investigation of the mesh at that time revealed severe distortion of elements and indicates the failure of the transfinite mesh to handle the case. Further investigation of the transfinite mesh motion technique and its shortcomings in solving the severe distortion Chapter 5. Numerical Examples 98 problems has been recently reported by Gadala et al. (2002). It should be mentioned that at some point in time, the overlap between crack faces exceed; o experimental notch o width. Therefore, a contact algorithm should have been considered for more realistic modeling. Figure (5.15) Snap shots of von Mises stress distribution during a mixed mode dynamic crack propagation. Figure (5.15) (continued) Snap shots of von Mises stress distribution during a mixed mode dynamic crack propagation. Chapter 6 CONCLUSIONS AND FUTURE WORK 6.1 S U M M A R Y O F A C C O M P L I S H M E N T S The aim of this work is to model the dynamic fracture problem in an efficient and reliable manner. The work accomplished in this regard may be summarized in the following points: - The A L E equations are formulated and implemented into a new finite element program (ALEFR), specifically developed to model dynamic fracture problems. The program is written in a structure consistent with the Lagrangian-based finite element method to ease the implementation into commercial codes. - In the A L E F R program, the process of crack propagation is modeled by a newly developed material point splitting (MPS) technique, instead of the traditional node release techniques widely used in the Lagrangian-based finite element methods. - The A L E boundary constraint is consistently satisfied using the True Boundary Tracking (TBT) technique, which was developed to control the motion of nodes on free boundaries, including the evolutional crack surfaces. 100 Chapter 6. Conclusions and Future Work 101 - The dynamic stress intensity factor and the energy release rate are computed from domain and contour integrals around the crack tip. The integration is carried out on a Lagrangian grid of material points. Properties of these material points are continuously updated. - The application of the developed techniques to problems of stationary cracks subjected to dynamic loading, demonstrates very good agreement with the analytical solution. - The process of mode I crack propagation has been successfully modeled. The propagating crack tip has been smoothly captured by a single node throughout the course of crack propagation. Therefore, the need for dense finite element mesh around the crack propagation path is eliminated and the choice of element size is independent of crack tip velocity or time step. - The onset of crack propagation was instantaneously captured by energy integrals and the trend of the solution afterwards is acceptable, however, the discrepancy relative to the analytical solutions is relatively large. The numerical results very slowly converge to the analytical solution with denser Lagrangian grids. 6.2 F U T U R E W O R K Future development may be directed towards the improvement of the current code, and the implementation and simulation of more dynamic fracture problems, as summarized in these points: - Further investigation and improvements are needed to enhance the convergence of the solution for the energy release rate and the dynamic stress intensity factors in dynamic crack propagation cases. This may include the improvement of the Chapter 6. Conclusions and Future Work 102 techniques used to update material properties and to evaluate the fracture measures. - The application to dynamic fracture models with nonlinear material behavior or inhomogeneous material may be investigated. Special treatment is needed when elements cross the boundary between materials, i.e., more than one material type exists in a single element. - Preliminary results obtained in modeling the in-plane mixed mode dynamic crack propagation are encouraging. Further work is, however, needed in the areas of extracting the dynamic stress intensity factors K] and Kn from the energy release rate integrals. This may include the use of modified forms of interaction integrals originally developed for static fracture cases. In order to maintain a more uniform mesh throughout the course of mixed mode crack propagation, the transfinite mapping scheme, which controls the motion of internal nodes, has to be improved or replaced by a more reliable mesh motion scheme, such as the isoparametric mapping. - A contact algorithm may be implemented in the current program. This will enhance the capability of the current program to include crack face closure cases, as has been seen in the mixed mode crack propagation example in Section 5.6.2, - Extension to 3-D analysis is possible since the developed A L E equations are general. This may include the use of 3-D brick and shell elements and the extension of the mesh motion scheme, the MPS and the TBT techniques and the procedures to evaluate the fracture measures to 3-D format. Bibliography Abdelgalil, A . I. and Gadala, M . S. (1999), "Modeling of Crack Propagation Using Arbitrary Lagrangian Eulerian Finite Elements", Proceedings of the 1st Canadian Conference on Nonlinear Solid Mechanics - CanCNSM'99, Victoria, B C , Canada, June 16, Vol . 2, pp. 356-365. Anderson, H. (1973), "The steadily Growing Elastic-Plastic Crack Tip in A Finite Element Treatment", International Journal of Fracture, Vol . 9, pp. 231-233. Atluri, S. N . and Kathiresan, K. (1980), "Influence of Flaw Shape on Stress Intensity Factors for Pressure Vessel Surface Flaws and Nossel Corners", International Journal of Pressure Vessel Technology, Vol . 102, pp. 278-286. Bathe, K. J. (1996), Finite Element Procedures in Engineering Analysis, Prentice-Hall. Bayoumi H.N. (2000), "Arbitrary Lagrangian-Eulerian Formulation for Quasi-static and Dynamic Metal Forming Simulation", Ph.D. Dissertation, The University of British Columbia. Bayoumi H.N. and Gadala M.S. (1999), "Finite Element Analysis of Large Strain Solid Mechanics Problems", Proceedings of the 1st Canadian Conference on Nonlinear Solid Mechanics - CanCNSM'99, Victoria, BC, Canada, June 16, Vol . 2, pp. 375-384. Bazant, Z. P., Glazik, J. I., Jr. and Achenbach, J. D. (1978), "Elastodynamic Fields Near Running Cracks by Finite Elements", Computers and Structures, Vol . 8, pp. 193-198. Caldis, E. S., Owen, D. R. J., Zienkiewicz, O. C. (1979), "Nonlinear Dynamic Transient Methods in Crack Propagation Studies", Nonlinear and Dynamic Fracture Mechanics, A S M E A M D , Vol . 35, pp. 1-17. 103 Bibliography 104 Erdogan, F. and Sih, G.C. (1963), "On the Crack Extension of Plates Under Plane Loading and Transverse Shear", Journal of Basic Engineering, Vol . 85, No. 4, pp. 519-527. Freund, L . B. (1973), "Crack Propagation in an Elastic Solid Subjected to General Loading- III Stress Wave Loading", Journal of the Mechanics and Physics of Solids, Vol. 21, pp. 47-61. Freund, L . B. (1990), Dynamic fracture Mechanics, First edition, Cambridge University Press. Gadala, M . S., Movahhedy, M . R. and Wang, J., (2002),"On the Mesh Motion for A L E Modeling of Metal Forming Processes", Finite Elements in Analysis and Design, Vol . 38, pp.435-459. Gadala, M.S. and Wang, J (1998), " A L E Formulation and Its Application in Solid Mechanics", Computer Methods in Applied Mechanical Engineering, Vol . 167, pp. 33-55. Gadala, M.S. and Wang, J. (1998), " A Practical Procedure for Mesh Motion in Arbitrary Lagrangian Eulerian Method", Engineering with computers, Vol . 14, pp. 223-234. Haber R., Shepard, M . S., Abel, J. F. Gallagher, R. H. and Greenberg, D. P. (1981), " A General Two-Dimensional, Graphical Finite Element Preprocessor Utilizing Discrete Transfinite Mappings", International Journal for Numerical Methods in Engineering, Vol . 17, pp. 1015-1044. Haber, R. B. (1984), " A Mixed Eulerian-Lagrangian Displacement Model for Large-Deformation Analysis in Solid Mechanics", Computer Methods in Applied Mechanics and Engineering, Vol . 43, pp. 277-292. Bibliography 105 Haber, R.B. (1984), " A Mixed Eulerian-Lagrangian Displacement Model for Large Deformation Analysis in Solid Mechanics", Computer Methods in Applied Mechanics and Engineering, Vol . 43, pp. 277-292. Hirt, C. W., Amsden, A . A. and Cook, J. L. (1974), "An Arbitrary Lagrangian-Eulerian Computing Method for A l l Flow Speeds", Journal of Computational Physics, Vol . 14, pp. 227-253. Hodulak, L. , Kobayashi, A. S. and Emery, A. F. (1980), " A Critical Examination of a Numerical Fracture Dynamic Code", Fracture Mechanics: twelfth Conference, A S T M STP 700, pp. 174-188. Huerta A . and Casadei, F. (1994), "New A L E Applications in Nonlinear Fast Transient Solid Dynamics", Engineering Computations, Vol . 11, pp. 317-345. Huetink, J. (1982), "Analysis of Metal Forming Processes Based on a Combined Eulerian-Lagrangian Finite Element Formulation", International Conference on Numerical Methods in Industrial Forming Processes, pp. 501-509. Hughes, T. J. R., Liu W. K. and Zimmermann, T. K. (1981), "Lagrangian-Eulerian Finite Element Formulation for Incompressible Viscous Flows", Computer Methods in Applied Mechanics and Engineering, Vol . 29, pp. 329-349. Hughes, T. J. R., Liu, W. K. and Zimmermann, T. K. (1981), "Lagrangian-Eulerian Finite Element Formulation for Incompressible Viscous Flows", Computer Methods in Applied Mechanics and Engineering, Vol . 29, pp. 329-349. Isaacson, E. and Keller, H . B. (1966), Analysis of Numerical Methods, John Wiley & Sons, New York. Bibliography 106 Kalthoff, J. F. Winkler, S. (1987), "Failure Mode Transition at High Rates of Shear Loading", International Conference on Impact Loading and Dynamic Behavior of Materials, Vol . l ,pp. 185-195. Kanninen, M . F. (1987), " A Critical Appraisal of Solutions Techniques in Dynamic Fracture Mechanics", Numerical Methods in Fracture Mechanics, Swansea, pp. 612-633. Keegstra, P. N . R., Head, J. L. and Turner, C. E. (1978), " A Two Dimensional Dynamic Linear Elastic Finite Element Program for the Analysis of Unstable Crack Propagation and Arrest", Numerical Methods in Fracture Mechanics, (Ed. Luxmoore, A . R. and Owen, D. R.), Swansea, pp. 634-647. Kobayashi, A . S. (1979), "Dynamic Fracture Analysis by Dynamic Finite Element Method - Generation and Propagation Analysis", Nonlinear and Dynamic Fracture Mechanics, A S M E A M D , Vol . 35, pp. 19-36. Koh, H. M . , Lee H . S. and Haber R. B. (1988), "Dynamic Crack Propagation Using Eulerian-Lagrangian Kinematic Descriptions", Computational Mechanics, Vol . 3, pp. 141-155. Lee, Y . J. and Freund, L .B . (1990), "Fracture Initiation Due to Asymmetric impact Loading of an Edge Cracked Plate", Journal of Applied Mechanics, Vol . 57, pp. 104-111. Liu, W. K., Chang, H. , Belytschko, T., and Chen, J. S. (1987), "Arbitrary Lagrangian-Eulerian Stress Update Procedures for Forming Simulations", Advances in Inelastic Analysis, A S M E - A M D , Vol . 88, pp. 153-175. Bibliography 107 Liu, W. K., Chang, FL, Chen, J. S. and Belytschko T. (1988), "Arbitrary Lagrangian-Eulerian Petrov-Galerkin Finite Elements for Nonlinear Continua", Computer Methods in Applied Mechanics and Engineering, Vol . 68, pp. 259-310. Malvern, L .E . (1969), "Introduction to the Mechanics of a Continuous Medium", Prentice Hall Inc., NJ, USA. Mahanty, D.K. and Maiti, S.K. (1990), "Experimental and Finite Elements Studies on Mode I and Mixed Mode (I&II) Stable Crack Growth", Engineering Fracture Mechanics, Vol . 37, No. 6, p. 1237-1275. Moran, B. and Shih C. F. (1987), "Crack Tip and Associated Domain Integrals From Momentum and Energy Balance", Engineering Fracture Mechanics, Vol . 27, No. 6, pp 615-641. Movahhedy, M . (2000), "ALE Simulation of Chip Formation in Orthogonal Metal Cutting Process", Ph.D. Dissertation, The University of British Columbia. Nashioka, T. and Atluri S. N . (1983), "Path-independent Integrals, Energy Release Rates, and General Solutions of Near-tip Fields in Mixed-mode Dynamic Fracture Mechanics", Engineering Mechanics, Vol . 18, No. 1, pp 1-22. Nishioka, T. and Atluri, S. N . (1980), "Numerical Modeling of Dynamic Crack Propagation infinite Bodies by moving Singular Elements, Part 1: Formulation", Journal of Applied Mechanics, Transactions of A S M E , Vol . 47, pp 570-576. Nishioka, T. and Atluri, S. N . (1986), "Computational Methods in Dynamic Fracture", In computational Methods in the Mechanics of Fracture (Edited by Atluri, S.N.), Elsevier, Amsterdam, pp. 335-383. Bibliography 108 Noh, W. F. (1964), " A Time Dependent Two-Space-Dimensional Coupled Eulerian-Lagrangian Code", Methods in Computational Physics, Vol . 3, pp. 117-179. Organ, D. J. (1996), "Numerical Solutions to dynamic Fracture Problems Using the Element-Free Galerkin Method', Ph.D. Dissertation, Northwestern University. Shepard, M.S., Yehia, N.A.B. , Burd, G.S. and Weidner,T.J. (1985), "Automatic Crack Propagation Tracking", Computers and Structures, Vol . 20, No. 1, pp. 211-223. Shih, C F . and Asaro, R.J. (1988), "Elastic Plastic Analysis of Cracks on Bimaterial Interfaces: Part I - Small Seal Yielding", Journal of Applied Mechanics, Vol . 55, pp. 299-316. Sih, G.C. (1974), "Strain Energy Density Factor Applied to Mixed Mode Crack Problems", International Journal of Fracture Mechanics, Vol . 10, pp. 305-321. Valliappan, S. and Murti, V. (1985), "Automatic Remeshing Technique In Quasi-Static And Dynamic Crack Propagation", Proceedings of Numeta'85 Conference, Swansea, U.K. , pp. 107-116. Wang, J. (1998), "Arbitrary Lagrangian Eulerian Method and Its Applications in Solid Mehanics", Ph.D. Dissertation, The University of British Columbia. Wang, J. and Gadala, M . S. (1997), "Formulation and Survey of A L E Method in Nonlinear Solid Mechanics", Finite Elements in Analysis and Design, Vol . 24, pp. 253-269. Wawrzynek, A . and Ingraffea, R. (1989), "An Interactive Approach to Local Remeshing Around A Propagating Crack", Finite Elements in Analysis and Design, Vol.5, pp. 87-96. Appendix A CALCULATION OF THE MATERIAL RATE OF CAUCHY STRESS The material rate of Cauchy stresses 'd\- is calculated from the material constitutive relation which is usually given in terms of an objective stress rate tensor such as the Truesdell stress rate tensor defined by (Bayoumi, 2000) w*<^-*Z«*-i%°> <A-1) The material constitutive relation in terms of the Truesdell stress rate is given by ' * J = ' < y D H (A.2) where ' Dtj is the rate of deformation tensor given by 1 d'v. d'y, D = — (——*- +—'-) (A.3) and 'Cijkl is the fourth order material constitutive tensor. 109 Appendix B TREATMENT OF CONVECTIVE TERMS The convective terms are treated according to the method developed by Bayoumi and Gadala (1999) and Bayoumi (2000). This method is based on fundamental A L E relations and it sidestep the computation of the spatial gradients of stresses. It, therefore, offers an accurate treatment of convective terms while maintaining the convenience of using displacement based finite elements. This method involves a transformation from volume integrals to surface integrals as offered by the divergence theorem. Use is also made of the boundary constraint in Equation (2.5). The last integral on the LHS of (2.22) may be rewritten as Applying the divergence theorem to the first integral in the RHS of Equation (B.l) and using Equation (2.5), we get d\(uk-ul)'cri.8te^l t , ode,:, dV-\(uk-ul)'cJij--^-'dV (B . l ) dV = \(uk -u*k)'cjySfy'nk'dS = 0 (B.2) 's Substituting in Equation (B.2) in Equation (B.l) yields 'V k 'v (B.3) 110 Appendix C DISCRETIZATION OF THE CONVECTIVE TERMS OF THE ALE EQUATION OF MOTION The convective terms of the A L E equation of motion (2.32) is discretized according to Bayoumi and Gadala (2000). The first convective velocity-stiffness virtual work term may be discretized as (C.l) where 'CM is given by L ' 2 i - l , 2 j - l l"2i-\,2j * - 2 i , 2 ; - l ^2i,2j (C.l) - 2Nx2N in which / and j indicate node numbers from 1 to N, and . dh N dh N ' C . 2 7 _ , = ' C £ 2 ; = J V ^ J ^ E ^ C ^ - ' ^ ) + ^ S ^ ( ' v ^ - ' v ^ ) ] W (C.3) ly OX k = \ 0y k=l 'rAl ='CM =0 (C4) The second convective velocity-stiffness virtual work term is discretized as (C.5) where 'CA2 is given by 111 Appendix C 112 ICA2 = 1 ("'2i-l,2;-l t^A2 ; ; i^A2 1 ^2i,2j-\ C 2 « , 2 y ! 2Nx2N in which z and j indicate node numbers from 1 to N, and 'V N k=i d'x ^dhktT7 , J tidy •v 'c£2H = j'ph.hjt^'dv ' V k=\ N d'x dhktT7 t. J T^dy 'V (C.6) (C.7) (C.8) (C9) ( C I O ) The second convective inertia force virtual work term may then be expressed as Slp(vr'v*)^Su'dV = (mT'CAUv •v dXj (CU) The third convective inertia force virtual work term may be discretized as J' •v where ' C A 3 is given by , , , 3 ' v , . 3 ' v f d'v. , T, ii, 'p(\-,v[)(^-2^)^SuiAt'dV = (dAi)T'CA3'v o xk o xk O Xj 'CA3 = \ L'2i-\,2j-\ tt-iATi ; '-21-1,27 ; j ts-<A3 • 1 l"2i,2H '-2/, 2j 1 (C.12) (C.13) 2Nx2N Appendix C 113 in which i and j indicate node numbers from 1 to N, and 'V dh N rlh  N dx k=i dy dh N dh N dyttdy M (C.14) L ' 2 / - l , 2 y — *-'2[,2 y - 1 (C.15) The fourth convective inertia force virtual work term may be discretized as f V ( ' v , - ' v f ) ( ' v , - ' v p - ^ - ( | ^ - ^ , . ) A ^ y = (du)T'CM' v J dxk oxj 'V where 'CM is given by , C A 4 = tsiA4 ts-,A4 ^2i-l 2 '-' 0 , H ^2i-\,2j t^A4 tpA4 J 2 / V x 2 J V in which / and j indicate node numbers from 1 to N, and 'cAi ='cAi = ['p\&^+h^)\t^H'v - » v « ) f ,v d'x d'x ' d'x' M 32 ^d'xd'y d'y d'x 1 d'xd'/h k xk V x k ) h k * yk dydy dy t f (C.16) (C.17) (C.18) * - '2 i - l ,2 j l-'2i,2y-1 "~ U (C19) 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0090649/manifest

Comment

Related Items