PHLEXcrack: A Meshless Code for Dynamic Fracture Propagation

Introduction

The Finite Element Method applied to fracture mechanics to date has been limited to the following approaches: All of these approaches are very computationally expensive whenever a reasonably accurate prediction of the crack path is required. The existence of the finite element mesh is the main obstacle to obtaining efficient solution algorithms. Recently, a new class of analysis techniques, based on so called meshless methods [8], has been developed. Besides the ease of domain discretization, the ability to introduce cracks arbitrarily in the domain is the main advantage of these methods when applied to fracture mechanics problems.

As expected, freedom from meshing constraints comes at a price, and meshless methods are usually significantly more computationally expensive than traditional finite element algorithms. To address these difficult issues, a comprehensive research program was initiated at Computational Mechanics Co, Inc. (Currently Altair Engineering, Inc.). The objective of the project, funded partially by the Office of Naval Research (years 1996-2000), was to develop robust numerical methods that avoid both the mesh related shortcomings of the FEM and the performance disadvantages of traditional meshless methods. The outcome of this effort is a new variant of meshless methods, called the Element Partition Method (EPM). This method is particularly suitable for fracture mechanics, and is computationally at least as efficient (and often surpassing) traditional finite element algorithms. A discontinuous version of EPM, specifically developed for fracture mechanics applications, produces accurate solutions of the stress field (static or dynamic) while the crack arbitrarily passes through the mesh cells.

The EPM technology has been implemented into a specialized crack propagation program, called PHLEXcrack. The current version of the code can be used to solve three-dimensional problems with linear physics (small deformation, small strain theory, static or dynamic problems), and with arbitrary (user specified) fracture equations to predict the crack direction and speed. The code is designed to be extensible, allowing developers/users to change the physical equations of the material as well as use different crack prediction theories.

Overview of the PHLEXcrack code

PHLEXcrack is currently released, and is available for the researchers and commercial applications, although the fracture mechanics part is limited to linear elasticity model only with or without thermal stresses. The EPM implementation in PHLEXcrack allows for the solution of three-dimensional problems using tetrahedral and hexahedral cells, with arbitrary p-enrichment controlled manually or by automatic adaptation based on a residual error estimator. Linear, static and dynamic solution sequences have been implemented, and the current crack propagation algorithm uses a Freund model [2], combined with stress intensity coefficients evaluated using linear fracture theory. Linear EPM volume elements can be connected to traditional Finite Element entinties (beams, shells, etc.), esp. for modeling external forces and/or the structure outside of the crack area. A typical hybrid application would be the analysis of a ship hull or aircraft structure, where only the vicinity of the crack is modeled using EPM, and the remainder of the model uses traditional shell elements.

The PHLEXcrack analysis package is based on the proprietary ProPHLEX[7] finite element/meshless development environment. This development environment was designed at Altair specifically for handling computer simulations in an hp-adaptive context. The current basic capabilities of ProPHLEX include:

Most of these capabilities are available in PHLEXcrack, although not all are available in the crack propagation solution sequence. In particular, the nonlinear solid mechanics module has not been fully tested in connection with EPM and is currently available only in an alpha testing mode. PHLEXcrack includes preliminary integration with the Underwater Shock Analysis [3] ``USA'' code, for accurate simulation of boundary loads resulting from for instance underwater explosions near immersed or surface structures.

Theoretical formulation of EPM

The Element Partition Method is a variant of the Partition of Unity (PoU) Method [1,5], which serves as a base for several implementations of modern meshless methods. The Partition of Unity is a set of linearly independent approximation functions, which only need to be able to approximate a constant function, and are thus quite easy to build. In order to obtain solution approximations of a required accuracy, PoU base functions are ``enriched'' typically using a set of polynomials, but any set of functions can be used for the enrichment. In particular, asymptotic solutions of fracture mechanics problems have proven to be very effective in improving the local accuracy of the solution.

In EPM, the base partition of unity comes from the traditional finite element approximation. This approximation can subsequently be arbitrarily enriched using higher order polynomials (p-enrichment) or specifically tailored functions based on the locally asymptotic solution near crack tip. The combination of these two enrichments produces very effective and accurate methods for solving fracture mechanics problems.

This version of EPM has many similiarities to the Finite Element Method (FEM): the PoU approximation satisfies in general the same requirements as FEM approximations, and they both use the Galerkin method to generate discretized equations. As a result, most of the theoretical foundations of FEM, including convergence and error estimates, apply directly or can easily be extended to EPM. Virtually all problems which can be solved using FEM, can be also analysed using EPM. In particular, multiple materials, dynamics, nonlinear solid mechanics and other similiar problems can be accurately modeled, with efficiency comparable to or better than FEM. In addition, EPM provides for an easier implementation of adaptivity, and its accuracy is much less affected by the quality of the meshing. In our testing we have successfully used meshes with very badly conditioned tetrahedral cells, while for a similiar mesh the FEM was locking even after enrichment to order p=3 .

EPM for domains with cracks

For the cells of the mesh (``elements'') which are cut by the crack surface a different version of the partition of unity is used. It is formulated using a Shepard Partition of Unity [6] within each ``cracked'' cell, effectively producing independent displacement fields on both sides of the crack without mesh refinement of remeshing. Still, the resulting approximation is mathematically correct, preserving the required continuity along inter--element boundaries, and satisfying rigid body motions and patch test criteria. To illustrate the idea, lets consider part of an EMP domain with a single crack

Example of element partition discretization cut by a crack

Cell T1 is not intersected by the crack, so the approximation at point x will use PoU based on standard shape functions. For each integration point (Gauss point) within the cells cut by the crack surface (T2 and T3), the segments between the point and the cell vertices are tested for the intersection with the crack surface (e.g. segments 6-y, 7-y, 10-y, and 11-y in the cell T2). Subsequently the solution at point y is fully defined by the nodal DOFs of ``visible'' nodes only, here nodes 6 and 7, and is independent of the solution at nodes 10 and 11, which are located on the other side of the crack surface. The solution along the path between points r and y is fully C0 continuous, while an arbitrary jump in the displacement field can appear between points y and z, or r and s.

Crack handling in PHLEXcrack

PHLEXcrack offers unprecedented flexibility in modeling the crack geometry. This is because crack surface definitions are handled in PHLEXcrack using a powerful geometry engine, a customized version of the one used in the popular FEM pre/postprocessor Altair's HyperMesh. Crack data inside the geometry engine is maintained as STL surfaces (a collection of triangular facets). The computation of geometrical quantities at the crack front (normal direction, curvature), visibility calculations (intersection of segments with the crack), and tolerance control for the accuracy of triangular representation, are handled internally within the geometry engine. In addition, the geometry engine automatically detects intersections of multiple crack surfaces and intersections with the boundary of the domain. Thus, the crack can be completely enclosed within the body, or can start and propagate along the boundary surface of the domain.


Flowchart of PHLEXcrack




In the dynamic / crack propagation solution sequence,  PHLEXcrack switches between two phases of the solution:

The crack advancement phase consists of the following steps:
  1. calculation of geometrical data at each point along the active crack front (tangent and normal directions, crack front curvature, etc.),
  2. calculating stress intensity factors and other relevant quantities at each point along the active crack front,
  3. local prediction of crack advancement,
  4. modification to the geometry database (crack advancement, adding / deleting vertices along the crack front to maintain required accuracy, possible reduction of features of the crack surface(s) behind the crack front),
  5. modification of the EPM database: marking newly cut cells, changes to the enrichment pattern due to the crack advancement, and
  6. recalculation of the stiffness and mass matrices for cells which were modified due to the crack movements.
Steps 1 and 4 are handled internally within the geometry engine. Steps 2 and 3 are specific to the particular theory of fracture propagation (PHLEXcrack is designed to allow for the substitution of user supplied modules to accomodate different equations of the crack motion). The last two steps required non-trivial modification of the time integration code, and were possible primarily because PHLEXcrack is built using Altair's proprietary ProPHLEX environment, which was designed for adaptive simulations.

Examples

During the validation process, we solved several simple examples of static and dynamic fracture problems, including dynamic crack propagation problems, and compared the results with solutions available in the literature. (A detailed report [4]  and published paper [9] are available) We present here a contrived example showing some of the fracture mechanics related capabilities of the code: a welded T-section of a beam with a preexisting ``half-penny'' crack:


Welded T-section beam with a crack.

Both the web and the flange are made of the same linear elastic material with the following properties:

The Freund propagation model is used to advance the crack.

EPM discretization for welded T-section beam.

Note than only minimal clustering of the integration cells has been used in the model. In fact, due to the shape of the crack, less than 20 elements are initially cut by the crack. Although the initial crack seems to be aligned with the nodes, PHLEXcrack automatically detected this situation and distorted all cells in the crack vicinity. This process is repeated after every crack advancement to avoid cell vertices being too close to the crack surface.

Following figures show the crack advancement at two distinct points in time. Note that crack surface was preserved during the computations (no defeaturing behind the crack front) so the history of the crack advancement is visible in both figures.
 
 


Crack surface of welded T-section example at time 0.0015


Crack surface of welded T-section example at time 0.0030


Conclusions and Suggested Future Work

PHLEXcrack is a very powerful tool to analyze various fracture mechanics problems, including: Note that due to the use of the crack definition independent of the domain discretization, PHLEXcrack can be used to rapidly analyze various ``what if'' cases with different crack configurations, without the need to regenerate the model of the structure. This could significantly reduce the manual effort and timeframe for the analysis of lifetime prediction for aging parts and structures.

Possible expansions of the code capabilities may include:

PHLEX environment already has most of the capabilities needed for nonlinear analysis, so these extensions should be possible with relatively small amount of extra effort.
 

References

[1] I. Babuska and J.M. Melenk. The partition of unity finite element method. Technical Report BN-1185, Inst. for Phys. Sc. and Tech., University of Maryland, June 1995.

[2] L.B. Freund, Dynamic Fracture Mechanics, Cambridge University Press, Cambridge, 1990.

[3] J.A. DeRuntz, Jr, Underwater Shock Analysis code, Theory Manual. Unique Software Applications, Colorado Springs, Colorado, 1998.

[4] T.J. Liszka, C.A.M. Duarte, and O.N. Hamzeh. Hp-meshless cloud method for dynamic structure in fluid-structure interaction, SBIR-ONR-N00014-96-C-0329. The Computational Mechanics Company, 1999.

[5] J.M. Melenk, and I. Babuska, The partition of unity finite element method: Basic theory and applications, Computer Methods in Applied Mechanics and Engineering, 139:289-314, 1996.

[6] D. Shepard, A two-dimensional function for irregularly spaced data. In ACM National Conference, pp 517-524, 1968.

[7] T.J.Liszka, et al. ProPHLEX - An hp-adaptive finite element kernel for solving coupled systems of partial differential equations in computational mechanics,  Computer Methods in Applied Mechanics and Engineering, 150:251-271, 1997.

[8] T.J.Liszka, C.A.M. Duarte, hp-meshless cloud method,  Computer Methods in Applied Mechanics and Engineering, 139:236-288, 1996.

[9] C.A.Duarte, O.N.Hamzeh, T.J.Liszka, W.W.Tworzydlo, A generalized finite element method for teh simulation of three-dimensional dynamic crack propagation ,  Computer Methods in Applied Mechanics and Engineering, 190:2227-2262, 2001.