Adaptive simulation of fracture processes based on spatial enumeration techniques

To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

Copyright © 1997 Elsevier Science Ltd Int. J. Rock Mech. & Min. Sci. Vol. 34, No. 3-4, 1997 ISSN 0148-9062 To cite this paper: Int. J. RockMech. &Min. Sci. 34:3-4, Paper No. 188

A D A P T I V E SIMULATION OF F R A C T U R E P R O C E S S E S BASED ON SPATIAL E N U M E R A T I O N T E C H N I Q U E S T. D. P. de A r a f i j o l; j . B. C a v a l c a n t e Net04; M . T. M . de C a r v a l h 0 2 ; T. N. B i t t e n c o u r t 3 ; L. F. M a r t h a 4

1 Departamento de Engenharia Civil, Pontificia Universidade Catdlica do Rio de Janeiro, (PUC-Rio), Rio de Janeiro, RJ, 22453-900, Brazil z TeCGraf-Grupo de Tecnologia em Computa~go Grfifica, Pontificia Universidade Catdlica do Rio de Janeiro, (PUC-Rio), Rio de Janeiro, RJ, 22453-900, Brazil 3 Laboratdrio de Mec~nica Computacional, Escola Polit6cnica da Universidade de Sgo Paulo, (EPUSP), Sgo Paulo, SP, 05508-900, Brazil 4 Departamento de Engenharia Civil, Pontificia Universidade Catdlica do Rio de Janeiro, (PUC-Rio), Rio de Janeiro, RJ, 22453-900, Brazil TeCGraf-Grupo de Tecnologia em Computa~go Grfifica, Pontificia Universidade Catdlica do Rio de Janeiro, (PUC-Rio), Rio de Janeiro, RJ, 22453-900, Brazil ABSTRACT

An interactive graphics computational system with self-adaptive, integrated, two-dimensional finite element analysis capabilities is described in this work. This system is able to handle both standard structural and fracture mechanics problems. The self-adaptive strategy is based on recursive spatial enumeration techniques: a binary tree partition for the boundary and the crack-line curves definition, and a quadtree partition for domain mesh generation. The 'apriori'refinement of the curves has the advantage of generating good transition meshes at the boundary regions. The system integrates different tools: a geometric modeler to create the model geometry, a pre-processor for mesh generation and attribute assignment, a numerical analysis module to evaluate the finite element response, and a post-processor for result visualization. The system is capable of deciding where to refine an initial mesh, of redoing the analysis, and of repeating this procedure until a pre-defined convergence criterion is achieved. Cracks can be introduced arbitrarily by the user at any position in the model. The system regenerates the meshes automatically taking into account the new created crack surfaces. For linear elastic analysis, quarter-point elements are inserted around the crack tips. The self-adaptive procedure is also considered in the crack propagation process. This procedure takes into account the arbitrarily generated crack geometry and the finite element error estimation analysis. Copyright © 1997 Elsevier Science Ltd

KEYWORDS Fracture • Self-Adaptive • Computer Graphics Programs • Crack Propagation • Mesh Refinement • Automatic Mesh Generation • Error Estimation • Delaunay Triangulation • Recursive Spatial Enumeration

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

INTRODUCTION The purpose of this work is to present a reliable, robust and efficient self-adaptive strategy that considers, adequately, the physical and geometric attributes of a structural model. The system implementing this strategy provides a regular mesh refinement for the free-boundary curves and for the interfaces between different materials. The adaptive process is based on 'a posteriori' error estimation. An h-refinement type is utilized in this process. An interactive graphics 10re-processor is used to generate the initial mesh information and the boundary conditions of the finite clement model. After analyzing the structural model, the system decides automatically where to refine the mesh. If it is necessary, the system refines the mesh considering the initial boundary conditions. A new analysis is then performed. This procedure is repeated until a pre-defined convergence criterion is achieved. The criterion adopted is a relative error in the global energy norm. It is important to mention that in each step of the cycle the user can interfere in the process and modify the specified relative error. The user can also visualize the responses using a graphics post-processor. The error estimator proposed by Zienckiewiecz, Zhu 1987 is employed here. This estimator considers a norm of the energy to measure the magnitude of the error. A nodal averaging scheme among adjacent elements is used to compute smoothed nodal stresses. This scheme is implemented in the numerical analysis module for all types of triangular or quadrilateral elements. The purpose of this work is not to question the accuracy of the chosen error estimator. Its selection has two main reasons: it is a well-known method and it is simple to implement. Different types of estimators could be implemented without any loss of generality of the proposed adaptive strategy. The mesh refinement is guided by a characteristic size of each clement, predicted according to a given error rate and the degree of the clement interpolation function. A recursive spatial enumeration procedure is used to perform the refinement in two levels: boundary refinement and domain refinement. The boundary refinement is based on a binary tree data structure while the domain refinement makes use of a quadtree data structure. The domain mesh generation combines the advantages of a quadtree based algorithm with a boundary contraction triangulation based on Delaunay properties. The quadtree triangulation is used in the entire domain except for a narrow region close to the boundary, which uses the Delaunay triangulation. To improve computational efficiency, the adjacency information provided by the quadtree data structure is also used in the boundary contraction algorithm. Many authors adopted a similar quadtree strategy for the self-adaptive process (Baehmann 1989). The advantage of the 'a priori'refinement of the boundary curves is that it insures mesh compatibility among distinct regions of the model. This also results in a mesh of better quality along the boundaries of the model. In crack propagation simulation, the geometry changes in each propagation step, therefore the global mesh has to be updated each time a crack moves. There arc some works in the literature that solve this problem with a local remeshing scheme, as Wawrzynek, Ingraffea 1987, Potyondy et al. 1995, and Bittencourt et al. 1996. In this work, however, a global remeshing scheme is employed. This global scheme relics on the efficiency of the mesh generation algorithm. Cracks can be introduced at any position and at any time, i.e., either in the initial structural model or at any stage of the adaptive process. After the introduction of a crack, the mesh is updated considering the new geometry.

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

R E C U R S I V E SPATIAL E N U M E R A T I O N A N D T R I A N G U L A T I O N Recursive spatial enumeration algorithms subdivide a particularly region of the space (two- or three-dimensional) into a number of smaller regions of a similar shape. The process is then repeated an arbitrary number of times. This information can be stored in a tree-type data structure. In the two-dimensional case, the basic shape is usually a square and each node of the tree has either four children or none. This data structure is called a quadtree (Samet 1984). The undivided regions are called cells. A cell size can be determined easily from its level in the tree. Algorithms which make use of the quadtree for mesh generation are computationally efficient, due to the properties of the quaternary tree. In addition, they provide a good transition between regions with different refinement orders. The main disadvantage presented by this technique, in its original version (Baehmann et al. 1987), is the boundary discretization in an irregular manner, creating new vertexes in non-controlled positions. This occurs because the boundary discretization is guided by the same tree. For that reason, the mesh generated by this technique cannot be compatible with a, previously supplied, fixed discretization for the boundary of the considered region. This fact makes it difficult to combine different mesh generation algorithms, especially in problems with local mesh modification. In the meshing algorithms based on boundary contraction, not only the boundary discretization is guaranteed but also the boundary points themselves are used as input data. There are many published boundary contraction algorithms, for example, Shaw 1978 and Lo 1989. The mesh generation process is developed in two phases: first, interior nodes are generated and second, a triangulation is performed. The triangulation is usually based on Delaunay criteria. However, the performance of the process is not very high because it involves a complexity of quadratic order in the worst case. There are also other algorithms that employ Delaunay triangulation schemes for regions with restrictions (no convex regions and with holes) whose performance is O(nlogn) (Chew 1989; De Florani et al. 1992). In this work, the finite element mesh is generated combining the quadtree procedure, the boundary contraction technique, and a Delaunay criterion. This combination tries to incorporate the efficiency and the capability of handling transitions of the quadtree technique with the advantages provided by the contraction algorithm (Cavalcante 1994).

SELF-ADAPTIVE STRATEGY Initially, the adaptive process requires a usually coarse finite element mesh for the model, the boundary and geometric description of this model, and their attributes. The geometric description is realized through curves that represent the boundary of the model and that define each model region (distinct materials, for example). The boundary condition description (support and/or loading) is associated to the curves and the domain parameters (material, properties, thickness, etc.) are associated to the regions. The process is started performing an analysis of the initial mesh. The results of this analysis, in addition to the description of regions, curves and attributes, are used for mesh refinement. A numerical analysis module is required to evaluate the discretization error. A finite element code, based on object oriented programming concepts (Martha et al. 1996), is employed for the implementations described in this work. The subsequent steps of the adaptive process start from the mesh of the previous step. One of the most important characteristics of the present self-adaptive meshing strategy is that the boundary refinement is enforced independently of the domain refinement. In fact, the domain

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

discretization requires 'apriori'boundary discretization. This enforces a better mesh gradation along the boundaries. In this boundary refinement, each boundary curve is discretized in its own parametric space. Therefore, the algorithm is generic for all classes of geometric curves. The algorithm used to refine each boundary curve is a unidimensional version of the algorithm that is used to refine the domain, which is based on a quadtree technique. Each curve is decomposed through a binary tree technique. The idea consists of recursively subdividing the curve into segments whose sizes are computed taking into account the characteristic sizes of the finite elements adjacent to the curve. These sizes are dictated by the error estimation analysis of the previous step in the adaptive mesh refinement. At the end of this phase, when all curves of the boundary are refined, the boundary conditions are reapplied to the model in a consistent way. For example, if a curve has a uniformly distributed load, the new equivalent nodal forces are calculated consistently with the new refinement. After the curves are discretized, the new mesh is generated as described in the previous section. Up to this point the error estimation analysis has been only considered in the boundary refinement. Therefore, the interior quadtree decomposition is further refined to also consider the error in the domain interior. This is performed visiting each element of the current mesh and testing its middle point with respect to the generated quadtree cells. If the cell which contains the element middle point has a size which is larger than the predicted characteristic element size, this cell is subdivided into four cells. This process is repeated until the cell size is less than the characteristic size. Following the traditional quadtree meshing approach (Baehmann et al. 1987), finite elements are generated in the interior cells using templates. These templates require that the difference in depth level for two adjacent cells is in the maximum one. The templates for quadrilateral elements were devised in such a way that the cells have twice the size of the cells for generating triangular elements. In the present case, this is considered to refine the quadtree cells taking into account the boundary discretization and the error estimation analysis. A major difference between the present algorithm and the traditional quadtree meshing is that only interior cells are considered for generating elements based on the quadtree decomposition. The narrow areas between the interior cells and the domain boundary are meshed in an unique process. In the present case, a boundary contraction procedure generates the mesh in the remaining area. No new interior node is generated in this process. A property of the Delaunay triangulation is used for the creation of triangular elements. Given a segment of the current boundary, the selection of a boundary node for the creation of a triangle is based on the maximum included angle. Because the boundary is not necessarily convex, additional checks are needed to avoid triangle overlapping (Shaw, Pitchen 1978). As mentioned previously, a problem with the boundary contraction technique is the quadratic complexity with respect of the number of mesh nodes. In the present case, the algorithm is enhanced by taking into account the existing quadtree domain decomposition. The selection of nodes for triangle or quadrilateral creation exploits the quadtree data structure to avoid testing of nodes that are not in the vicinity of the boundary segment in consideration. This is certainly one of the most important gains of the present algorithm when compared to previous boundary contraction procedures. The final step of the mesh generation is a node coordinate smoothing by averaging the coordinates of adjacent nodes.

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

C R A C K INITIATION AND PROPAGATION In the case of fracture simulation, a crack is arbitrarily introduced in the model in any phase of the process. A new geometry is generated and the previous mesh is deleted. A new mesh is generated utilizing the same algorithm described previously, using the boundary nodes of the previous mesh as input boundary information. In the moment, only triangular (T6) elements are used, although the adaptive process also uses quadrilateral elements. The two crack surfaces are considered geometrically coincident. The crack trajectory is treated like any geometry curve and before the system regenerates the mesh, the curve is subdivided and the subdivision nodes are added to the input boundary description together with the previous mesh nodes. To insure the generation of well-shaped elements at each crack tip, a uniform rosette pattern of quarter-point singular elements is inserted around each tip. The meshing algorithm just considers the rosette outline and its nodes are considered at the region boundary (Wawrzynek, Ingraffea 1987). In the mesh generation, crack curve nodes are considered as any other boundary node, i.e., there is only one instance of each subdivision node along a crack curve. In the boundary contraction procedure each crack curve is traversed twice: one for each direction of the curve. This provides a consistent boundary information for the meshing algorithm. For the finite element analysis module, the crack curve nodes are duplicated, except the crack tip nodes. For the analysis, the generated elements consider the appropriate nodes on each side of a crack curve. In the adaptive process, during the boundary refinement, the curve defining the crack is discretized considering the adjacent elements in each side of the crack. The numerical error is not considered for the quarter-point elements, neither to discretize the boundary nor the domain. The size of these elements are kept basically the same during the whole process.

EXAMPLE The present strategy was implemented in an interactive graphics system that supports simulations of discrete crack propagation along with the adaptive process. The system is totally interactive, so the user can interfere in the process at any time. In an adaptive simulation, the user can modify the control error parameter and, in crack propagation simulation, the user also controls the number of subdivisions of the crack curves. In both types of simulation, the mesh is generated automatically. All phases of the simulation strategy described in this work are illustrated by means of the following example. Figure 1 shows the model and its initial finite element mesh. This mesh, the geometric description and attribute information are taken as an input data set for the adaptive system. After interpreting the initial results, the user decides where to insert a crack and does it interactively (Figure 2). Defined the crack curve, the system asks for the number of segments that the user wants to subdivide it. The number of segments guide the definition of the size of the rosette elements (Figure 3). This definition is automatically performed by the system. After the crack curve and the rosette elements are defined, the system remeshes the model. To generate the new mesh, the previous mesh boundary nodes, the crack curve nodes, and the rosette nodes are used as input boundary information. The new mesh, in the example, is depicted in Figure 4. Generated the new mesh, a numerical analysis is required to produce results and also to estimate the

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

error. From these parameters, the system refines the mesh utilizing the strategy described in the self-adaptive process. Figure 5 shows the refined mesh after one step of the adaptive process, i.e., the error obtained with this mesh is not smaller than the one furnished by the user. Figure 6 shows the deformed model after one step of analysis. Crack propagation follows similar procedures. The system supports any number of cracks. Figure 7 illustrates a case where an edge crack and an internal crack are inserted simultaneously into the model. The deformed configuration of the model was obtained after analyzing the structure. CONCLUSIONS An interactive adaptive strategy for structural analysis of discrete crack propagation and of standard structural problems was described. This strategy is currently applicable to two-dimensional linear elastic finite element models. The resulting self-adaptive system possesses the following characteristics: robustness, computational efficiency, and reliability. As the geometry of the model changes after each step of crack growth, the finite element mesh has to be modified accordingly. The global mesh generation scheme has proven applicable in crack propagation simulations of large and complex models used in practice. The mesh generation algorithm can be applied to any type of geometry, with any number of cracks. In the case of fracture mechanics simulations only triangular meshes are currently generated. However, an easy extension can be made in order to generate quadrilateral or mixed meshes, as it is already available for standard structural problems.

FIGURES

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Mm. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

Paper 188, Figure 1. I1..'.

., ...

I I:

<

• .':

i

I ...

.:.

.....

~1:

Figure 1. Initial mesh of the model.

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Mm. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

Paper 188, Figure 2.

t, tt-'lt

Figure 2. Crack curve defined by the user.

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Mm. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

Paper 188, Figure 3.

,--.,-.~ 1. ~ .11. ~ ~"

1 Figure 3. Inserting crack curve and quarter-point elements

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Mm. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

Paper 188, Figure 4.

Figure 4. New mesh after the crack is inserted.

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Mm. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

Paper 188, Figure 5.

Figure 5. Refined mesh by the self-adaptive process.

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Mm. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

Paper 188, Figure 6.

Figure 6. Deformed configuration after analysis.

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

Paper 188, Figure 7.

Figure 7. Deformed configuration with multiple cracks.

References References Baehmann P. L., Wittchen S. L., Shepard M. S. et al. 1987. Robust geometrically based, automatic two-dimensional mesh generation. International Journal for Numerical Methods in Engineering, 24, 1043-1078. Baehmann P. L. 1989. Automated finite element modeling and simulation, Ph.D. Thesis, Rensselaer Polytechnic Institute, Troy-New York. Barsoum R. S. 1976. On the use of isoparametric elements in linear fracture mechanics. International Journal for Numerical Methods in Engineering, 10, 25-37. Bittencourt T. N., Wawrzynek R A., Ingraffea A. R., Sousa J. L. 1996. Quasi-automatic simulation of crack propagation for 2D LEFM problems. Engineering Fracture Mechanics, 55:2, 321-334. Cavalcante Neto J. B. 1994. Simula~fio self-adaptativa baseada em enumera~fio espacial recursiva de modelos bidimensionais de elementos finitos, Msc. Thesis, Pontificia Universidade Cat61ica do Rio de Janeiro. Chew L. R 1989. Constrained Delaunay triangulation. Algorithmic', 4, 97-108. De Florani L., Puppo E. 1992. An on-line algorithm for constrained Delaunay triangulation. Computer Vision, Graphics, and Image Processing, 54, 290-300.

ISSN 0148-9062

To cite this paper: Int. J. Rock Mech. & Min. Sci. 34:3-4, paper No. 188.

Copyright © 1997 Elsevier Science Ltd

Henshell R. D., Shaw K. G. 1975. Crack tip elements are unnecessary. International Journal for Numerical Methods' in Engineering, 9, 495-509. Lo S.H. 1989. Delaunay triangulation of non-convex planar domains. International Journal for Numerical Methods' in Engineering, 28, 2695-2707. Martha L. F., Menezes I. F. M., Lages E. M., Parente Jr. E., Pitangueira R. L. 1996. An OOP class organization for materially nonlinear FE analysis. Proceedings ofXVII CILAMCE, Venice, Italy, 229-232. Potyondy D. O., Wawrzynek R A., Ingraffea A. R. 1995. An algorithm to generate quadrilateral or triangular element surface meshes in arbitrary domains with applications to crack propagation. International Journal for Numerical Methods' in Engineering, 38, 2677-2701. Samet H. 1984. The quadtree and related hierarchical data structures. ACM Computer Surveys, 16:2, 187-260. Shaw R. D., Pitchen R. G. 1978. Modifications to the Suhara-Fukuda method of network generation. International Journal for Numerical Methods' in Engineering, 12, 93-99. Wawrzynek RA., Ingraffea A. R. 1987. Interactive finite element analysis of fracture processes: an integrated approach. Theo~ & AppL Frac. Mech., 8, 137-150. Zienkiewicz O. C., Zhu J. Z. 1987. A simple error estimator and adaptive procedure for practical engineering analysis. International Journal for Numerical Methods' in Engineering, 24, 337-357.

ISSN 0148-9062