- Email: [email protected]
Adaptive Delaunay triangulation with object-oriented programming for crack propagation analysis
Finite Elements in Analysis and Design 40 (2004) 1753 – 1771 www.elsevier.com/locate/!nel
Adaptive Delaunay triangulation with object-oriented programming for crack propagation analysis Sutthisak Phongthanapanich, Pramote Dechaumphai∗ Department of Mechanical Engineering, Faculty of Engineering, Chulalongkorn University, Bangkok 10330, Thailand Received 19 March 2003; received in revised form 19 December 2003
Abstract A !nite element method, with the adaptive Delaunay triangulation as mesh generator, is used to analyze two-dimensional crack propagation problems. This paper describes the Delaunay triangulation procedure consisting of mesh generation, node creation, mesh smoothing, and adaptive remeshing, all with object-oriented programming. The adaptive remeshing technique generates small elements around crack tips and large elements in other regions. The resulting stress intensity factors and simulated crack propagation behavior are used to evaluate the e4ectiveness of the combined method. The accuracy of the stress intensity factor prediction is evaluated in three test cases, a center cracked plate, a single edge cracked plate and a compact tension specimen. Then, crack growth trajectories in a single edge cracked plate with three holes and a single edge cracked plate under mixed-mode loading are simulated and results assessed. ? 2004 Elsevier B.V. All rights reserved. Keywords: Delaunay triangulation; Adaptive mesh; Object-oriented programming; Crack propagation
1. Introduction The Delaunay triangulation, based on the concept of the Voronoi diagram [1,2], is one of the automated mesh generation algorithms that has recently gained popularity due to the ability to generate meshes of arbitrary geometry for both simply connected and multi-boundary domains. The procedure that is capable to produce meshes with regular nodal density and triangular elements for arbitrary two-dimensional geometry was !rst introduced by Weatherill and Hassan [3] and programmed as object-oriented codes by Karamete et al. [4]. In this paper, triangular mesh construction by the Delaunay triangulation for crack propagation analysis is described in details. In addition, an adaptive ∗
Corresponding author. Tel./fax: +66-2-218-6621. E-mail address: [email protected] (P. Dechaumphai).
0168-874X/$ - see front matter ? 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.!nel.2004.01.002
1754
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
remeshing technique is developed and incorporated into the triangulation in order to improve the solution accuracy of the !nite element method. The technique generates an entirely new mesh based on the solution obtained from the previous mesh; elements in regions with large changes of solution gradients are re!ned and elements in areas with little changes of gradients grow larger. The object-oriented based pseudo-code procedures for the Delaunay triangulation algorithm, the automatic node creation procedure and the adaptive remeshing technique are presented. The stress intensity factor is a critical parameter in the prediction of fatigue crack growths in crack propagation problems. In this work, the standard six-node isoparametric elements are used to model test specimens. In order to improve the accuracy of the near-tip stress !elds, elements around the crack tip have their mid-side nodes displaced from their nominal positions to quarter points [5,6]. The stress intensity factors are then calculated from the extrapolated nodal displacements on the crack faces next to the crack tip [7]. For crack propagation simulation, both the global remeshing [8] and the local mesh re!nements [9,10], which simply reduce element sizes near the crack tip while leaving other regions unchanged, have been proposed. This paper employs a global remeshing scheme with automatic mesh adaptation near crack tips using the error equidistribution property to determine appropiate element sizes. Because an entirely new mesh is generated for the whole domain, the initial mesh does not a4ect the solution accuracy. Small elements are automatically placed so that they clustered near crack tips, whilst larger elements are generated in the other regions. Thus, the scheme provides solution accuracy at crack tips, and reduces the computational e4ort. Several examples under mode I and mixed mode loadings, are modeled to evaluate the accuracy and e4ectiveness of the combined procedure. In addition, the capability of the proposed procedure is further demonstrated by the simulation of crack propagation trajectories in a single edge cracked plate with three holes and a single edge cracked plate under mixed-mode loading. The predicted solutions of the test cases are assessed by comparing with experimental and numerical results, respectively.
2. Stress intensity factor, crack propagation and nite element method The important parameters used in the linear elastic fracture mechanics are the stress intensity factors in various modes. Several methods have been proposed to determine the stress intensity factors, such as the displacement extrapolation near the crack tip [7], the J -integral [11], and the energy domain integral [12]. In this paper, the !eld extrapolation near the crack tip is used [7,13]. The stress intensity factors are determined from E (vc − ve ) 2 KI = 4(vb − vd ) − ; (1a) 3(1 + v)(1 + ) L 2 2 E (uc − ue ) KII = 4(ub − ud ) − ; (1b) 3(1 + v)(1 + ) L 2 where E is the modulus of elasticity, is the Poisson’s ratio, is the elastic parameter de!ned by (3 − 4) for plane strain and (3 − )=(1 + ) for plane stress problems, and L is the element length. The u and v are the displacement components in the x and y directions, respectively; the subscripts indicate their positions as shown in Fig. 1.
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1755
Fig. 1. Quarter-point triangular elements around the crack tip.
Fig. 2. A quarter-point six-node triangular element.
Crack propagation normally occurs under mixed mode loadings. There are several methods use to predict the direction of crack growths such as the maximum circumferential stress theory [14], the maximum energy release rate theory [15], and the minimum strain energy density theory [16]. In the maximum circumferential stress theory, the direction of crack propagation may be computed from KI sin + KII (3 cos − 1) = 0:
(2)
Eq. (2) implies that the crack propagates at zero angle for the pure mode I, and at non-zero angle in a mixed mode loading. Finite element equations for determining nodal displacements and stresses are derived from the governing partial di4erential equations that represent the equilibrium conditions. These equations can be written in matrix form as [K]{} = {F};
(3)
where {} is the vector that contains unknowns of the element nodal displacements {F} is the load vector, and [K] is the sti4ness matrix [17]. To felicitate high resolution accuracy near the crack tip, the six-node triangular elements that form up a circular zone surrounding the crack tip have their mid-side nodes displaced from their nominal positions to quarter points of the tip [6] as shown in Fig. 2. The radius of the circular zone is speci!ed to be no larger than one-eight of the initial crack length, with roughly one element every 30◦ in the circumferential direction [13]. The crack growth simulation is based upon the maximum
1756
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
circumferential stress [14]; the increment of the crack length that increases during each step as the crack propagates is speci!ed by the user. 3. Delaunay triangulation for crack propagation analysis For a given set of points in space {Pk }, k = 1; : : : ; n, the regions {Vk }, k = 1; : : : ; n, are Voronoi boundaries assigned to each point Pk and represent the space closer to Pk than to any other points in the set. Therefore, these convex regions satisfy Vk = {Pi :|p − Pi | ¡ |p − Pj |; ∀j = i}:
(4)
The above de!nition implies that boundaries of the Voronoi diagram must lie half way between the two points on either side of the boundary. If all points which have some segments of a Voronoi boundary in common are joined, the resulting shape is a Delaunay triangulation as shown in Fig. 3. In graph theory, Delaunay triangulation can be de!ned by the empty circumcircle property in which any circle in the plane is said to be empty if it contains no vertex in its interior. 3.1. Mesh generation procedure The Delaunay triangulation algorithm for constructing boundary triangles is based on the in-circle criterion [1,2]. In this paper, the algorithm is extended to crack propagation analysis with special isosceles six-node triangles around the crack tip as depicted in Fig. 4. The details of the Delaunay triangulation algorithm are described in Refs. [18,19]. The algorithm is implemented by a computer program written in VB.NET and is summarized as a pseudo-code in the Algorithm 1 with variables and parameters de!ned in Section 3.4. Algorithm 1. DelaunayTriangulation (P; T; p) Let P0 be the collection of node objects; Let T 0 be the collection of mesh objects; P0. Initialize; T 0. Initialize; t ← T . FindTriangleContainNode (p); T 0 ← T . IncircleTriangles (t; p); P0 ← T 0. DestroyTriangles (); T . CreateNewTriangles (P0; p); T . AssignNeighborhoodTriangles; End; 3.2. Automatic node creation procedure The Delaunay triangulation algorithm in the previous section does not explain the node creation method inside the domain. So far, researchers have introduced several approaches for creating new
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1757
Fig. 3. Delaunay triangulation layer over a Voronoi diagram.
Fig. 4. Removal of the crack tip node and creation of rosette nodes around the crack tip.
nodes inside the domain by re!ning boundary triangles such that the set of boundary points guide new node placements [20,21]. The new node creation procedure for geometry with crack propagations presented in this paper is extended from that proposed by Refs. [3,4]. The shapes and sizes of triangles are controlled by two coeKcients, the Alpha and the Beta. The Alpha coeKcient controls node density by changing the allowable shape of the formed triangles while the Beta coeKcient regulates the triangulation uniformly by disallowing new node insertion within a speci!ed distance of each other in the same sweep of the triangles. The automatic node creation procedure searches for the triangle that conforms to both Alpha and Beta testing criteria and inserts a new node at the centroid of that triangle. New triangles can then be created by Delaunay triangulation algorithm as described in Algorithm 1. The key idea of the procedure is summarized as pseudo-code in the Algorithm 2 below; Algorithm 2. MeshRefinement (P; T; alpha; beta; iteration) Let P0 be the collection of node objects; For i=l To iteration { Do t ← T . NextTriangle { p ← t. ComputeTriangleCentroid ( ); p.dp← t. ComputePointDistributionFunction ( ); p.dm(1 : 3) ← t. DistanceCentroidToVertices ( ); p. Rejected = FALSE; For j=l To 3 {
1758
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 5. Point distribution functions and automatic node creation scheme (Algorithm 2).
If (p. dm(j)¡ (alpha∗ p.dp)) { p. Rejected = TRUE; Break; };
}; If (Not p. Rejected) { P0. Initialize; P0 ← T . FindInsertedNodeOfNearestTriangles; Do p1 ← P0. NextNode { If (distance (p; p1) ¡ (beta∗ p. dp)) { p. Rejected = TRUE; Break; }; }; }; If (Not p. Rejected) P. AddNodeAsInsertedNode (p);
}; Do p ← P. NextInsertedNode { }; Call DelaunayTriangulation (P; T; p); };
}; End;
The point distribution function, dpi for the node, pi , is computed from dpi =
M 1 |pj − pi |; M j=1
(5)
where node i is surrounded by M neighbouring nodes (see Fig. 5). Shapes and sizes of triangles formed from the previous step can be improved by applying a several passes of mesh Laplacian smoothing technique through the entire set of all interior nodes. The node repositioning formula is derived from the !nite di4erence approximation of the Laplace’s equation. Each interior node is successively moved to the centroid of the area which is formed by connecting
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
neighbouring nodes. The new node locations are computed from M M yi i=1 xi and yic = i=1 i = 1; 2; : : : ; M; xic = M M
1759
(6)
where xi and yi are the coordinates of the surrounding M nodes. A demonstration of a new node creation inside a triangle with automatic node creation procedure (Algorithm 2) and the formation of new triangles by Delaunay triangulation algorithm (Algorithm 1) are shown in Fig. 5. Fig. 6 shows the di4erent stages of the domain discretization for the compact tension specimen with an outside crack. Fig. 7 shows detailed rosette triangular elements around the crack tip in the !nal mesh. 3.3. Adaptive remeshing technique The adaptive remeshing technique generates an entirely new mesh based on the solution obtained from a previous mesh [22,23]. In this paper, the technique is modi!ed and corporated into the Delaunay triangulation and the !nite element method for crack propagation analysis. There are two main steps in the implementation of the adaptive remeshing technique, the determination of proper element sizes and the new mesh generation. To determine appropriate element sizes at di4erent locations in the domain, the von Mises stress & is used as the indicator for computing proper element sizes. As small elements are required in the region where changes in the von Mises stress gradients are large, the second derivatives of the von Mises stress at a point with respect to global coordinates x and y are needed. Using the concept of principal directions from a given state of stresses at a point, the principal quantities in the principal directions X and Y where the cross-derivatives vanish are determined: 2 2 9& 92 & 9& 0 9x2 9X 2 9x9y : (7) ⇒ 2 9& 92 & 92 & 0 9x9y 9y2 9Y 2 The maximum principal quantities are then used to compute the proper element size, hi , by requiring the error to be uniform for all elements [22], h2i )i = h2min )max = constant; where )i is the higher principal quantity of the element that is considered, 2 2 9 & 9 & )i = max 2 ; 2 : 9X 9Y
(8)
(9)
In Eq. (8), )max is the maximum principal quantity of all elements and hmin is the minimum element size speci!ed by users. The mesh regeneration, based on the concepts of the Delaunay triangulation, the mesh re!nement and the adaptive remeshing technique, constructs a new mesh using the information from the previous
1760
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 6. Re!nement and smoothing of a mesh for compact tension specimen.
mesh (background mesh). A node may be inserted in an element which contains the background node with smaller point distribution function than hmax and the average of the distance to the three vertices of that element. With this technique, the new mesh is composed of small elements in the regions with large changes in solution gradients, and large elements in other regions where the changes in solution gradients are small. The process of adaptive remeshing technique is summarized as pseudo-code below;
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1761
Fig. 7. Rosette triangular elements around the crack tip in the compact tension specimen.
Algorithm 3. AdaptiveRemeshing (P; T; P0; Hmin; Hmax; threshold) Do { Do p ← P0. NextInteriorNode { If (p.hi6 Hmax) { t ← T . FindTriangleContainNode (p); pq ← t. ComputeTriangleCentroid (); pq.dp← t. ComputePointDistributionFunction (); pq.dm(1 : 3) ← t. DistanceCentroidToVertices (); pq.Rejected = FALSE; For j=1 To 3 { If (pq.hi¿ pq. dm. Average Or pq.dm(j)¡ Hmin) { pq.Rejected = TRUE; Break; }; }; If (Not pq.Rejected) P.AddNodeAsInsertedNode (pq); }; }; Do p ← P.NextInsertedNode { Call DelaunayTriangulation (P; T; p) ; }; } Loop Until (P.InsertedNodes <= threshold); End; The above Algorithm 3 can be sequentially described as follows: 1. Let P0 be the set of nodes of the background mesh. Let P be the set of nodes and T be the set of triangles of the current mesh. 2. Specify the values of hmin and hmax , respectively, as the mesh re!nement and enlargement criteria. 3. Read next interior node, pi , of the background mesh from P0. 4. If hi ¿ hmax , then return to step 3. 5. Search for triangle ti in T which contains the node pi by starting from the triangle which was last formed and used Lawson’s algorithm [18,19] to march from one triangle to the next in the direction of pi . This algorithm searches for the path and removes the need to search through the
1762
6. 7. 8. 9. 10.
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
entire domain. Then calculate the centroid of the triangle ti , de!ne it as node pq , and compute the nodal distribution value of node pq by Eq. (5). Compute the distance dm , m = 1; 2; 3, from node pq to each of the three vertices of the triangle ti . If hi ¿ average of dm or hmin ¿ dm , m = 1; 2; 3, then return to step 3. Otherwise accept the node pq for insertion by the main engine of Delaunay triangulation algorithm (Algorithm 1) and add node pq to P. Repeat steps 3 to 8 until all interior nodes in P are considered. Perform the Delaunay triangulation of the inserted nodes in P by Algorithm 1.
3.4. Mesh generation implementation This section presents the main algorithm for implementing the mesh generation from the Delaunay triangulation (Algorithm 1), the mesh re!nement procedure (Algorithm 2), and the adaptive remeshing technique (Algorithm 3) together. The combined program is developed using the object-oriented programming concept that takes into account the advantages of the code reusable, inheritance, and polymorphism capability. In addition, the main algorithm has incorporated a scheme for mesh generation around the crack tips. The crack tip nodes are !rst removed from the node-list prior to the domain discretization. Such process avoids the modi!cation of the original Delaunay triangulation (Algorithm 1) and the mesh re!nement procedure (Algorithm 2) previously described. Then the rosette nodes around the crack tips, with speci!ed angle by the parameter angle increment, are generated and added to the node-list. After the domain has been discretized, the crack tip nodes are inserted to form the rosette elements around the crack tips. The implementation of the main algorithm is summarized in the Algorithm 4 below. Algorithm 4. Main (P; T; angle increment; alpha; beta; iteration; Hmin; Hmax; threshold; isadaptive) Let BP be the collection of boundary node objects that stored in sequence of counter-clockwise direction for all outside boundaries and clockwise direction for all inside boundaries; Let P0 be the collection of background node objects; Let P be the collection of node objects; Let T be the collection of mesh objects; Let angle increment be the isosceles triangle crack tip angle; Let alpha be the constant that controls shape of formed triangles; Let beta be the constant that controls regularity of the triangulation; Let iteration be the number of loops to refine meshes; Let Hmin and Hmax be the minimum and maximum element size, respectively; Let threshold be the number of minimum increasing nodes for each iteration; Let isadaptive be the flag to generate background or adaptive mesh; BP. Initialize; P0. Initialize;
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1763
P. Initialize; T . Initialize; If (isadaptive) { P0. ReadBackgroundNodes; BP. RediscretizeBoundaryNodes; }; Else { BP. ReadBoundaryNodes; }; BP. MarkCrackTipNodes; BP. CreateRosetteNodesAroundCrackTips (angle increment); BP. CreateConvexHull; P. AddNode (BP:p1; BP:p2; BP:p3; BP:p4); T . AddTriangle (t1; BP:p1; BP:p2; BP:p3); T . AddTriangle (t2; BP:p3; BP:p2; BP:p4); Do p ← BP. NextBoundaryNode { Call DelaunayTriangulation (P; T; p); }; Call MeshRefinement (P; T; alpha; beta; iteration); If (isadaptive) Call AdaptiveRemeshing (P; T; P0; Hmin; Hmax; threshold); T . LaplaceSmoothing; T . RemoveOutsideDomainTriangles; Do p ← BP. NextCracktipNode { T . CreateRosetteTriangularElements (p); }; End; 4. Algorithm evaluation The simulation of fracture mechanics problems with !nite element program is used to evaluate the eKciency of the combined Delaunay triangulation and the adaptive remeshing technique. The procedure !rst predicts the stress intensity factors for problems with analytical solutions or experimental data so that their results can be compared. The procedure is then employed to capture the crack trajectory by adapting the mesh automatically with the crack growth. 4.1. Determination of stress intensity factors Three well-known geometries used in the evaluation of the proposed procedure are: (1) the center cracked plate, (2) the single edge cracked plate, and (3) the compact tension specimen. The center cracked plate: The geometry of the center cracked plate and its !nal adaptive mesh are shown in Fig. 8. The plate has an initial crack length 2a = 100 units, and the thickness t = 1
1764
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 8. Problem statement and the !nal mesh of the center cracked plate.
Fig. 9. Problem statement and the !nal mesh of the single edge cracked plate.
unit. The stress intensity factor for this problem can be calculated from [24] √ KI = 1:334& a:
(10)
Both the full and a quarter models are used to analyze the problem. The computed stress intensity factor, KI , from either model is 16.7611 comparing to 16.7192 from Eq. (10) with the value di4erence of less than 0:25%. The single edge cracked plate: The geometry of the single edge cracked plate and its !nal adaptive mesh are shown in Fig. 9. The plate has an initial crack length a = 0:4 unit, the width b = 1 unit, and the thickness t = 1 unit. The stress intensity factor of 2.357 was calculated from [25] √ KI = F& a; (11)
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1765
Fig. 10. Problem statement and the !nal mesh of the compact tension specimen.
where F = 1:12 − 0:2311 + 10:5512 − 21:7213 + 30:3914
and
a 1= : b
The computed stress intensity factor from the adaptive mesh shown in Fig. 9 is 2.358, which is less than 0:05% di4erent from the closed-form solution. The compact tension specimen: The geometry of the compact tension specimen with an initial crack length a = 3 mm, the width w = 50:8 mm and the thickness t = 25:4 mm, and its !nal adaptive mesh are shown in Fig. 10. The stress intensity factor can be calculated from [26]
a
a 2
a 0:886 + 4:64 − 13:32 KI = P 2 + w w w
a 3
a 4 √ a 3=2 t w 1− + 14:72 − 5:6 : (12) w w w The computed stress intensity factor obtained from the adaptive model is 27.718 comparing to 27.935 from Eq. (12) with 0:75% di4erence. 4.2. Simulation of crack propagation Fig. 11 describes the Oow-chart for predicting the crack growth trajectory. The combined adaptive meshing technique and the crack propagation analysis according to the !gure can be described in order as follows: 1. Create geometry model with initial crack de!ned by users. 2. Discretize the domain into triangular elements. Then apply the boundary conditions (loads and constraints) onto the !nite element model. 3. Solve the system of equations using Eq. (3) for displacement solutions. Then compute the stress intensity factors using Eqs. (1)(a)–(b). 4. For the initial mesh or the stress intensity factors as compared to the previous values change signi!cantly, go to step 5. Otherwise, go to step 6.
1766
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 11. Crack trajectory prediction Oow-chart.
5. Determine new element sizes according to the mesh adaptation process using Eqs. (8)–(9). Delete the old mesh and generate a new mesh with the computed element sizes. Apply the boundary conditions (loads and constraints) onto the new !nite element model and go to step 3. 6. If the crack tip has reached the !nal target position or the stability condition [10] is met, then stop the process. Otherwise determine the crack direction using Eq. (2) and update the model with the new crack geometry. The crack propagation increment used in this paper is between 20% and 50% of the initial crack length, depending on the KII =KI ratio. Then go to step 2. The above procedure has been implemented as a computer program that integrates the Delaunay mesh generation algorithm and the adaptive remeshing technique with a !nite element method together. The performance of the combined procedure is evaluated by simulating the crack growth trajectory of a single edge cracked plate with three holes under mixed mode loading as shown in Fig. 12. The plate is simply supported near the lower corners, and subjected to a concentrated load at the center of the upper edge. The experiments were conducted by Bittencourt et al. [10] with two cases of the initial crack length, a, and its location, b, as shown by the table in Fig. 12. For the !rst case, the initial crack length, a, and its location, b, are 1.5 and 5.0 units, respectively. The results of the adaptive !nite element meshes and the crack growth trajectory are depicted in Fig. 13. The !gure shows that the crack growth trajectory passes just above the lower hole and ended at the
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1767
Fig. 12. Problem statement for the single edge crack plate with three holes.
Fig. 13. Adaptive !nite element meshes and the crack growth trajectory for the single edge cracked plate with three holes under mixed-mode loading (Case I).
1768
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 14. Adaptive !nite element meshes and the crack growth trajectory for the single edge cracked plate with three holes under mixed-mode loading (Case II).
Fig. 15. Problem statement and the !nal mesh of the initial crack for the single edge crack plate.
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1769
Fig. 16. Adaptive !nite element meshes and the crack growth trajectory for the single edge cracked plate under mixed-mode loading.
middle hole. For the second case with the initial crack length, a = 1:0 and location b = 4:0 units, respectively, the crack propagates toward the middle hole (Fig. 14). In both cases, the predicted crack growth trajectories very closely resemble the experimental results. Fig. 15 shows the problem statement and !nal mesh of a single edge cracked plate under mixed-mode loading. This plate is !xed at the bottom edge and is subjected to far-!eld shear stress 3 = 1 unit along the top edge. The plate has an initial crack length a = 3:5 units. The modulus of elasticity and the Poisson’s ratio are 30 × 106 units and 0.25, respectively. The plane strain condition is assumed in the analysis. The computed stress intensity factors KI and KII from the adaptive mesh are 34.10 and 4.52 comparing to the reference values of 34.00 and 4.55, respectively. Fig. 16 shows the predicted crack trajectory that closely resembles the numerical solution in Ref. [27]. 5. Conclusions The adaptive !nite element method using Delaunay triangulation for crack propagation analysis was presented. The concept of Delaunay triangulation for two dimensional mesh construction was explained. The mesh generation procedure with automatic node creation and mesh smoothing, which was implemented by object-oriented programming, were presented in details. The technique was combined with the !nite element method for analyzing crack problems under single and mixed mode loadings. The isoparametric six-node triangular elements, with mid-side nodes displaced from their nominal positions to quarter points of the crack tip, were employed to form up a circular zone surrounding the tip in order to better capture the stress !eld. The solution accuracy was further improved by incorporating an adaptive remeshing technique into the Delaunay triangulation algorithm. The adaptive remeshing technique places small elements around the crack tips and in regions with large changes of stress gradients. At the same time, larger elements are generated in other regions to minimize the total number of unknowns and the
1770
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
computational time. Several examples were employed to evaluate the accuracy of the combined Delaunay triangulation, the !nite element method, and the adaptive remeshing technique by predicting the stress intensity factors for several geometries with di4erent loadings, as well as capturing the crack propagation trajectory. These examples demonstrated the capability of the combined adaptive Delaunay triangulation with the !nite element method for solving crack propagation problems e4ectively. Acknowledgements The authors are pleased to acknowledge the Thailand Research Fund (TRF) for supporting this research work. References [1] A. Bowyer, Computing Dirichlet tessellations, Comp. J. 24 (2) (1981) 162–166. [2] D.F. Watson, Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes, Comp. J. 24 (2) (1981) 167–172. [3] N.P. Weatherill, O. Hassan, EKcient three-dimension Delaunay triangulation with automatic point creation and imposed boundary constraints, Int. J. Numer. Methods Eng. 37 (12) (1994) 2005–2039. [4] B.K. Karamete, T. Tokdemir, M. Ger, Unstructured grid generation and a simple triangulation algorithm for arbitrary 2-D geometries using object oriented programming, Int. J. Numer. Methods Eng. 40 (2) (1997) 251–268. [5] R.S. Barsoum, On the use of isoparametric !nite elements in linear fracture mechanics, Int. J. Numer. Methods Eng. 10 (1) (1976) 25–37. [6] R.S. Barsoum, Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements, Int. J. Numer. Methods Eng. 11 (1) (1977) 85–98. [7] S.K. Chan, I.S. Tuba, W.K. Wilson, On the !nite element method in linear fracture mechanics, Eng. Fract. Mech. 2 (1) (1970) 1–17. [8] T.D. de Araujo, D. Roehl, T.N. Bittencourt, L.F. Martha, Adaptive simulation of elastic-plastic fracture process, new trends and applications, in: S. Idelsohn, E. Onate, E. Dvorkin (Eds.), Fourth World Congress on Computational Mechanics, Buenos Aires, Argentina, 1998. [9] D.O. Potyondy, P.A. Wawrzynek, A.R. Ingra4ea, An algorithm to generate quadrilateral or triangular element surface meshes in arbitrary domains with applications to crack propagation, Int. J. Numer. Methods Eng. 38 (10) (1995) 2677–2701. [10] T.N. Bittencourt, P.A. Wawrzynek, A.R. Ingra4ea, J.L. Sousa, Quasi-automatic simulation of crack propagation for 2D LEFM problems, Eng. Fract. Mech. 55 (2) (1996) 321–334. [11] D.M. Parks, A sti4ness derivative !nite element technique for determination of crack tip stress intensity factors, Int. J. Fract. 10 (4) (1974) 487–502. [12] B. Moran, C.F. Shih, A general treatment of crack tip contour integrals, Eng. Fract. Mech. 35 (1987) 295–310. [13] G.V. Guinea, J. Planas, M. Elices, KI evaluation by the displacement extrapolation technique, Eng. Fract. Mech. 66 (3) (2000) 243–255. [14] F. Erdogan, G.C. Sih, On the crack extension in plates under plane loading and transverse shear, J. Basic Eng. 85 (1963) 519–527. [15] R.J. Nuismer, An energy release rate criterion for mixed mode fracture, Int. J. Fract. 11 (2) (1975) 245–250. [16] G.C. Sih, Strain-energy-density factor applied to mixed-mode crack problems, Int. J. Fract. 10 (1974) 305–321. [17] O.C. Zienkiewicz, R.L. Taylor, Finite Element Method, 5th Edition, Butterworth-Heinemann, Woburn, 2000. [18] C.L. Lawson, Software for C1 surface interpolation, in: J.R. Rice (Ed.), Mathematical Software III, Academic Press, New York, 1977.
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1771
[19] S.W. Sloan, A fast algorithm for generating constrained Delaunay triangulations, Comp. Struct. 47 (3) (1993) 441–450. [20] W.H. Frey, Selective re!nement: a new strategy for automatic node placement in graded triangular meshes, Int. J. Numer. Methods Eng. 24 (11) (1987) 2183–2200. [21] H. Borouchaki, P.L. George, Aspects of 2-D Delaunay mesh generation, Int. J. Numer. Methods Eng. 40 (11) (1997) 1957–1975. [22] J. Peraire, M. Vahdati, K. Morgan, O.C. Zienkiewicz, Adaptive remeshing for compressible Oow computations, J. Comput. Phys. 72 (2) (1987) 449–466. [23] P. Dechaumphai, Improvement of plane stress solutions using adaptive !nite elements, J. Chin. Inst. Eng. 19 (3) (1996) 375–380. [24] M. Isida, E4ect of width and length on stress intensity factors of internally cracked plates under various boundary conditions, Int. J. Fract. 7 (1971) 301–316. [25] Y. Murakami, Stress Intensity Factors Handbook, Pergamon Press, Oxford, 1987. [26] ASTM 1996, Annual Book of ASTM Standards, American Society for Testing and Materials, Philadelphia, 1996. [27] B.N. Rao, S. Rahman, An eKcient meshless method for fracture analysis of cracks, Comput. Mech. 26 (4) (2000) 398–408.
Adaptive Delaunay triangulation with object-oriented programming for crack propagation analysis Sutthisak Phongthanapanich, Pramote Dechaumphai∗ Department of Mechanical Engineering, Faculty of Engineering, Chulalongkorn University, Bangkok 10330, Thailand Received 19 March 2003; received in revised form 19 December 2003
Abstract A !nite element method, with the adaptive Delaunay triangulation as mesh generator, is used to analyze two-dimensional crack propagation problems. This paper describes the Delaunay triangulation procedure consisting of mesh generation, node creation, mesh smoothing, and adaptive remeshing, all with object-oriented programming. The adaptive remeshing technique generates small elements around crack tips and large elements in other regions. The resulting stress intensity factors and simulated crack propagation behavior are used to evaluate the e4ectiveness of the combined method. The accuracy of the stress intensity factor prediction is evaluated in three test cases, a center cracked plate, a single edge cracked plate and a compact tension specimen. Then, crack growth trajectories in a single edge cracked plate with three holes and a single edge cracked plate under mixed-mode loading are simulated and results assessed. ? 2004 Elsevier B.V. All rights reserved. Keywords: Delaunay triangulation; Adaptive mesh; Object-oriented programming; Crack propagation
1. Introduction The Delaunay triangulation, based on the concept of the Voronoi diagram [1,2], is one of the automated mesh generation algorithms that has recently gained popularity due to the ability to generate meshes of arbitrary geometry for both simply connected and multi-boundary domains. The procedure that is capable to produce meshes with regular nodal density and triangular elements for arbitrary two-dimensional geometry was !rst introduced by Weatherill and Hassan [3] and programmed as object-oriented codes by Karamete et al. [4]. In this paper, triangular mesh construction by the Delaunay triangulation for crack propagation analysis is described in details. In addition, an adaptive ∗
Corresponding author. Tel./fax: +66-2-218-6621. E-mail address: [email protected] (P. Dechaumphai).
0168-874X/$ - see front matter ? 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.!nel.2004.01.002
1754
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
remeshing technique is developed and incorporated into the triangulation in order to improve the solution accuracy of the !nite element method. The technique generates an entirely new mesh based on the solution obtained from the previous mesh; elements in regions with large changes of solution gradients are re!ned and elements in areas with little changes of gradients grow larger. The object-oriented based pseudo-code procedures for the Delaunay triangulation algorithm, the automatic node creation procedure and the adaptive remeshing technique are presented. The stress intensity factor is a critical parameter in the prediction of fatigue crack growths in crack propagation problems. In this work, the standard six-node isoparametric elements are used to model test specimens. In order to improve the accuracy of the near-tip stress !elds, elements around the crack tip have their mid-side nodes displaced from their nominal positions to quarter points [5,6]. The stress intensity factors are then calculated from the extrapolated nodal displacements on the crack faces next to the crack tip [7]. For crack propagation simulation, both the global remeshing [8] and the local mesh re!nements [9,10], which simply reduce element sizes near the crack tip while leaving other regions unchanged, have been proposed. This paper employs a global remeshing scheme with automatic mesh adaptation near crack tips using the error equidistribution property to determine appropiate element sizes. Because an entirely new mesh is generated for the whole domain, the initial mesh does not a4ect the solution accuracy. Small elements are automatically placed so that they clustered near crack tips, whilst larger elements are generated in the other regions. Thus, the scheme provides solution accuracy at crack tips, and reduces the computational e4ort. Several examples under mode I and mixed mode loadings, are modeled to evaluate the accuracy and e4ectiveness of the combined procedure. In addition, the capability of the proposed procedure is further demonstrated by the simulation of crack propagation trajectories in a single edge cracked plate with three holes and a single edge cracked plate under mixed-mode loading. The predicted solutions of the test cases are assessed by comparing with experimental and numerical results, respectively.
2. Stress intensity factor, crack propagation and nite element method The important parameters used in the linear elastic fracture mechanics are the stress intensity factors in various modes. Several methods have been proposed to determine the stress intensity factors, such as the displacement extrapolation near the crack tip [7], the J -integral [11], and the energy domain integral [12]. In this paper, the !eld extrapolation near the crack tip is used [7,13]. The stress intensity factors are determined from E (vc − ve ) 2 KI = 4(vb − vd ) − ; (1a) 3(1 + v)(1 + ) L 2 2 E (uc − ue ) KII = 4(ub − ud ) − ; (1b) 3(1 + v)(1 + ) L 2 where E is the modulus of elasticity, is the Poisson’s ratio, is the elastic parameter de!ned by (3 − 4) for plane strain and (3 − )=(1 + ) for plane stress problems, and L is the element length. The u and v are the displacement components in the x and y directions, respectively; the subscripts indicate their positions as shown in Fig. 1.
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1755
Fig. 1. Quarter-point triangular elements around the crack tip.
Fig. 2. A quarter-point six-node triangular element.
Crack propagation normally occurs under mixed mode loadings. There are several methods use to predict the direction of crack growths such as the maximum circumferential stress theory [14], the maximum energy release rate theory [15], and the minimum strain energy density theory [16]. In the maximum circumferential stress theory, the direction of crack propagation may be computed from KI sin + KII (3 cos − 1) = 0:
(2)
Eq. (2) implies that the crack propagates at zero angle for the pure mode I, and at non-zero angle in a mixed mode loading. Finite element equations for determining nodal displacements and stresses are derived from the governing partial di4erential equations that represent the equilibrium conditions. These equations can be written in matrix form as [K]{} = {F};
(3)
where {} is the vector that contains unknowns of the element nodal displacements {F} is the load vector, and [K] is the sti4ness matrix [17]. To felicitate high resolution accuracy near the crack tip, the six-node triangular elements that form up a circular zone surrounding the crack tip have their mid-side nodes displaced from their nominal positions to quarter points of the tip [6] as shown in Fig. 2. The radius of the circular zone is speci!ed to be no larger than one-eight of the initial crack length, with roughly one element every 30◦ in the circumferential direction [13]. The crack growth simulation is based upon the maximum
1756
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
circumferential stress [14]; the increment of the crack length that increases during each step as the crack propagates is speci!ed by the user. 3. Delaunay triangulation for crack propagation analysis For a given set of points in space {Pk }, k = 1; : : : ; n, the regions {Vk }, k = 1; : : : ; n, are Voronoi boundaries assigned to each point Pk and represent the space closer to Pk than to any other points in the set. Therefore, these convex regions satisfy Vk = {Pi :|p − Pi | ¡ |p − Pj |; ∀j = i}:
(4)
The above de!nition implies that boundaries of the Voronoi diagram must lie half way between the two points on either side of the boundary. If all points which have some segments of a Voronoi boundary in common are joined, the resulting shape is a Delaunay triangulation as shown in Fig. 3. In graph theory, Delaunay triangulation can be de!ned by the empty circumcircle property in which any circle in the plane is said to be empty if it contains no vertex in its interior. 3.1. Mesh generation procedure The Delaunay triangulation algorithm for constructing boundary triangles is based on the in-circle criterion [1,2]. In this paper, the algorithm is extended to crack propagation analysis with special isosceles six-node triangles around the crack tip as depicted in Fig. 4. The details of the Delaunay triangulation algorithm are described in Refs. [18,19]. The algorithm is implemented by a computer program written in VB.NET and is summarized as a pseudo-code in the Algorithm 1 with variables and parameters de!ned in Section 3.4. Algorithm 1. DelaunayTriangulation (P; T; p) Let P0 be the collection of node objects; Let T 0 be the collection of mesh objects; P0. Initialize; T 0. Initialize; t ← T . FindTriangleContainNode (p); T 0 ← T . IncircleTriangles (t; p); P0 ← T 0. DestroyTriangles (); T . CreateNewTriangles (P0; p); T . AssignNeighborhoodTriangles; End; 3.2. Automatic node creation procedure The Delaunay triangulation algorithm in the previous section does not explain the node creation method inside the domain. So far, researchers have introduced several approaches for creating new
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1757
Fig. 3. Delaunay triangulation layer over a Voronoi diagram.
Fig. 4. Removal of the crack tip node and creation of rosette nodes around the crack tip.
nodes inside the domain by re!ning boundary triangles such that the set of boundary points guide new node placements [20,21]. The new node creation procedure for geometry with crack propagations presented in this paper is extended from that proposed by Refs. [3,4]. The shapes and sizes of triangles are controlled by two coeKcients, the Alpha and the Beta. The Alpha coeKcient controls node density by changing the allowable shape of the formed triangles while the Beta coeKcient regulates the triangulation uniformly by disallowing new node insertion within a speci!ed distance of each other in the same sweep of the triangles. The automatic node creation procedure searches for the triangle that conforms to both Alpha and Beta testing criteria and inserts a new node at the centroid of that triangle. New triangles can then be created by Delaunay triangulation algorithm as described in Algorithm 1. The key idea of the procedure is summarized as pseudo-code in the Algorithm 2 below; Algorithm 2. MeshRefinement (P; T; alpha; beta; iteration) Let P0 be the collection of node objects; For i=l To iteration { Do t ← T . NextTriangle { p ← t. ComputeTriangleCentroid ( ); p.dp← t. ComputePointDistributionFunction ( ); p.dm(1 : 3) ← t. DistanceCentroidToVertices ( ); p. Rejected = FALSE; For j=l To 3 {
1758
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 5. Point distribution functions and automatic node creation scheme (Algorithm 2).
If (p. dm(j)¡ (alpha∗ p.dp)) { p. Rejected = TRUE; Break; };
}; If (Not p. Rejected) { P0. Initialize; P0 ← T . FindInsertedNodeOfNearestTriangles; Do p1 ← P0. NextNode { If (distance (p; p1) ¡ (beta∗ p. dp)) { p. Rejected = TRUE; Break; }; }; }; If (Not p. Rejected) P. AddNodeAsInsertedNode (p);
}; Do p ← P. NextInsertedNode { }; Call DelaunayTriangulation (P; T; p); };
}; End;
The point distribution function, dpi for the node, pi , is computed from dpi =
M 1 |pj − pi |; M j=1
(5)
where node i is surrounded by M neighbouring nodes (see Fig. 5). Shapes and sizes of triangles formed from the previous step can be improved by applying a several passes of mesh Laplacian smoothing technique through the entire set of all interior nodes. The node repositioning formula is derived from the !nite di4erence approximation of the Laplace’s equation. Each interior node is successively moved to the centroid of the area which is formed by connecting
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
neighbouring nodes. The new node locations are computed from M M yi i=1 xi and yic = i=1 i = 1; 2; : : : ; M; xic = M M
1759
(6)
where xi and yi are the coordinates of the surrounding M nodes. A demonstration of a new node creation inside a triangle with automatic node creation procedure (Algorithm 2) and the formation of new triangles by Delaunay triangulation algorithm (Algorithm 1) are shown in Fig. 5. Fig. 6 shows the di4erent stages of the domain discretization for the compact tension specimen with an outside crack. Fig. 7 shows detailed rosette triangular elements around the crack tip in the !nal mesh. 3.3. Adaptive remeshing technique The adaptive remeshing technique generates an entirely new mesh based on the solution obtained from a previous mesh [22,23]. In this paper, the technique is modi!ed and corporated into the Delaunay triangulation and the !nite element method for crack propagation analysis. There are two main steps in the implementation of the adaptive remeshing technique, the determination of proper element sizes and the new mesh generation. To determine appropriate element sizes at di4erent locations in the domain, the von Mises stress & is used as the indicator for computing proper element sizes. As small elements are required in the region where changes in the von Mises stress gradients are large, the second derivatives of the von Mises stress at a point with respect to global coordinates x and y are needed. Using the concept of principal directions from a given state of stresses at a point, the principal quantities in the principal directions X and Y where the cross-derivatives vanish are determined: 2 2 9& 92 & 9& 0 9x2 9X 2 9x9y : (7) ⇒ 2 9& 92 & 92 & 0 9x9y 9y2 9Y 2 The maximum principal quantities are then used to compute the proper element size, hi , by requiring the error to be uniform for all elements [22], h2i )i = h2min )max = constant; where )i is the higher principal quantity of the element that is considered, 2 2 9 & 9 & )i = max 2 ; 2 : 9X 9Y
(8)
(9)
In Eq. (8), )max is the maximum principal quantity of all elements and hmin is the minimum element size speci!ed by users. The mesh regeneration, based on the concepts of the Delaunay triangulation, the mesh re!nement and the adaptive remeshing technique, constructs a new mesh using the information from the previous
1760
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 6. Re!nement and smoothing of a mesh for compact tension specimen.
mesh (background mesh). A node may be inserted in an element which contains the background node with smaller point distribution function than hmax and the average of the distance to the three vertices of that element. With this technique, the new mesh is composed of small elements in the regions with large changes in solution gradients, and large elements in other regions where the changes in solution gradients are small. The process of adaptive remeshing technique is summarized as pseudo-code below;
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1761
Fig. 7. Rosette triangular elements around the crack tip in the compact tension specimen.
Algorithm 3. AdaptiveRemeshing (P; T; P0; Hmin; Hmax; threshold) Do { Do p ← P0. NextInteriorNode { If (p.hi6 Hmax) { t ← T . FindTriangleContainNode (p); pq ← t. ComputeTriangleCentroid (); pq.dp← t. ComputePointDistributionFunction (); pq.dm(1 : 3) ← t. DistanceCentroidToVertices (); pq.Rejected = FALSE; For j=1 To 3 { If (pq.hi¿ pq. dm. Average Or pq.dm(j)¡ Hmin) { pq.Rejected = TRUE; Break; }; }; If (Not pq.Rejected) P.AddNodeAsInsertedNode (pq); }; }; Do p ← P.NextInsertedNode { Call DelaunayTriangulation (P; T; p) ; }; } Loop Until (P.InsertedNodes <= threshold); End; The above Algorithm 3 can be sequentially described as follows: 1. Let P0 be the set of nodes of the background mesh. Let P be the set of nodes and T be the set of triangles of the current mesh. 2. Specify the values of hmin and hmax , respectively, as the mesh re!nement and enlargement criteria. 3. Read next interior node, pi , of the background mesh from P0. 4. If hi ¿ hmax , then return to step 3. 5. Search for triangle ti in T which contains the node pi by starting from the triangle which was last formed and used Lawson’s algorithm [18,19] to march from one triangle to the next in the direction of pi . This algorithm searches for the path and removes the need to search through the
1762
6. 7. 8. 9. 10.
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
entire domain. Then calculate the centroid of the triangle ti , de!ne it as node pq , and compute the nodal distribution value of node pq by Eq. (5). Compute the distance dm , m = 1; 2; 3, from node pq to each of the three vertices of the triangle ti . If hi ¿ average of dm or hmin ¿ dm , m = 1; 2; 3, then return to step 3. Otherwise accept the node pq for insertion by the main engine of Delaunay triangulation algorithm (Algorithm 1) and add node pq to P. Repeat steps 3 to 8 until all interior nodes in P are considered. Perform the Delaunay triangulation of the inserted nodes in P by Algorithm 1.
3.4. Mesh generation implementation This section presents the main algorithm for implementing the mesh generation from the Delaunay triangulation (Algorithm 1), the mesh re!nement procedure (Algorithm 2), and the adaptive remeshing technique (Algorithm 3) together. The combined program is developed using the object-oriented programming concept that takes into account the advantages of the code reusable, inheritance, and polymorphism capability. In addition, the main algorithm has incorporated a scheme for mesh generation around the crack tips. The crack tip nodes are !rst removed from the node-list prior to the domain discretization. Such process avoids the modi!cation of the original Delaunay triangulation (Algorithm 1) and the mesh re!nement procedure (Algorithm 2) previously described. Then the rosette nodes around the crack tips, with speci!ed angle by the parameter angle increment, are generated and added to the node-list. After the domain has been discretized, the crack tip nodes are inserted to form the rosette elements around the crack tips. The implementation of the main algorithm is summarized in the Algorithm 4 below. Algorithm 4. Main (P; T; angle increment; alpha; beta; iteration; Hmin; Hmax; threshold; isadaptive) Let BP be the collection of boundary node objects that stored in sequence of counter-clockwise direction for all outside boundaries and clockwise direction for all inside boundaries; Let P0 be the collection of background node objects; Let P be the collection of node objects; Let T be the collection of mesh objects; Let angle increment be the isosceles triangle crack tip angle; Let alpha be the constant that controls shape of formed triangles; Let beta be the constant that controls regularity of the triangulation; Let iteration be the number of loops to refine meshes; Let Hmin and Hmax be the minimum and maximum element size, respectively; Let threshold be the number of minimum increasing nodes for each iteration; Let isadaptive be the flag to generate background or adaptive mesh; BP. Initialize; P0. Initialize;
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1763
P. Initialize; T . Initialize; If (isadaptive) { P0. ReadBackgroundNodes; BP. RediscretizeBoundaryNodes; }; Else { BP. ReadBoundaryNodes; }; BP. MarkCrackTipNodes; BP. CreateRosetteNodesAroundCrackTips (angle increment); BP. CreateConvexHull; P. AddNode (BP:p1; BP:p2; BP:p3; BP:p4); T . AddTriangle (t1; BP:p1; BP:p2; BP:p3); T . AddTriangle (t2; BP:p3; BP:p2; BP:p4); Do p ← BP. NextBoundaryNode { Call DelaunayTriangulation (P; T; p); }; Call MeshRefinement (P; T; alpha; beta; iteration); If (isadaptive) Call AdaptiveRemeshing (P; T; P0; Hmin; Hmax; threshold); T . LaplaceSmoothing; T . RemoveOutsideDomainTriangles; Do p ← BP. NextCracktipNode { T . CreateRosetteTriangularElements (p); }; End; 4. Algorithm evaluation The simulation of fracture mechanics problems with !nite element program is used to evaluate the eKciency of the combined Delaunay triangulation and the adaptive remeshing technique. The procedure !rst predicts the stress intensity factors for problems with analytical solutions or experimental data so that their results can be compared. The procedure is then employed to capture the crack trajectory by adapting the mesh automatically with the crack growth. 4.1. Determination of stress intensity factors Three well-known geometries used in the evaluation of the proposed procedure are: (1) the center cracked plate, (2) the single edge cracked plate, and (3) the compact tension specimen. The center cracked plate: The geometry of the center cracked plate and its !nal adaptive mesh are shown in Fig. 8. The plate has an initial crack length 2a = 100 units, and the thickness t = 1
1764
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 8. Problem statement and the !nal mesh of the center cracked plate.
Fig. 9. Problem statement and the !nal mesh of the single edge cracked plate.
unit. The stress intensity factor for this problem can be calculated from [24] √ KI = 1:334& a:
(10)
Both the full and a quarter models are used to analyze the problem. The computed stress intensity factor, KI , from either model is 16.7611 comparing to 16.7192 from Eq. (10) with the value di4erence of less than 0:25%. The single edge cracked plate: The geometry of the single edge cracked plate and its !nal adaptive mesh are shown in Fig. 9. The plate has an initial crack length a = 0:4 unit, the width b = 1 unit, and the thickness t = 1 unit. The stress intensity factor of 2.357 was calculated from [25] √ KI = F& a; (11)
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1765
Fig. 10. Problem statement and the !nal mesh of the compact tension specimen.
where F = 1:12 − 0:2311 + 10:5512 − 21:7213 + 30:3914
and
a 1= : b
The computed stress intensity factor from the adaptive mesh shown in Fig. 9 is 2.358, which is less than 0:05% di4erent from the closed-form solution. The compact tension specimen: The geometry of the compact tension specimen with an initial crack length a = 3 mm, the width w = 50:8 mm and the thickness t = 25:4 mm, and its !nal adaptive mesh are shown in Fig. 10. The stress intensity factor can be calculated from [26]
a
a 2
a 0:886 + 4:64 − 13:32 KI = P 2 + w w w
a 3
a 4 √ a 3=2 t w 1− + 14:72 − 5:6 : (12) w w w The computed stress intensity factor obtained from the adaptive model is 27.718 comparing to 27.935 from Eq. (12) with 0:75% di4erence. 4.2. Simulation of crack propagation Fig. 11 describes the Oow-chart for predicting the crack growth trajectory. The combined adaptive meshing technique and the crack propagation analysis according to the !gure can be described in order as follows: 1. Create geometry model with initial crack de!ned by users. 2. Discretize the domain into triangular elements. Then apply the boundary conditions (loads and constraints) onto the !nite element model. 3. Solve the system of equations using Eq. (3) for displacement solutions. Then compute the stress intensity factors using Eqs. (1)(a)–(b). 4. For the initial mesh or the stress intensity factors as compared to the previous values change signi!cantly, go to step 5. Otherwise, go to step 6.
1766
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 11. Crack trajectory prediction Oow-chart.
5. Determine new element sizes according to the mesh adaptation process using Eqs. (8)–(9). Delete the old mesh and generate a new mesh with the computed element sizes. Apply the boundary conditions (loads and constraints) onto the new !nite element model and go to step 3. 6. If the crack tip has reached the !nal target position or the stability condition [10] is met, then stop the process. Otherwise determine the crack direction using Eq. (2) and update the model with the new crack geometry. The crack propagation increment used in this paper is between 20% and 50% of the initial crack length, depending on the KII =KI ratio. Then go to step 2. The above procedure has been implemented as a computer program that integrates the Delaunay mesh generation algorithm and the adaptive remeshing technique with a !nite element method together. The performance of the combined procedure is evaluated by simulating the crack growth trajectory of a single edge cracked plate with three holes under mixed mode loading as shown in Fig. 12. The plate is simply supported near the lower corners, and subjected to a concentrated load at the center of the upper edge. The experiments were conducted by Bittencourt et al. [10] with two cases of the initial crack length, a, and its location, b, as shown by the table in Fig. 12. For the !rst case, the initial crack length, a, and its location, b, are 1.5 and 5.0 units, respectively. The results of the adaptive !nite element meshes and the crack growth trajectory are depicted in Fig. 13. The !gure shows that the crack growth trajectory passes just above the lower hole and ended at the
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1767
Fig. 12. Problem statement for the single edge crack plate with three holes.
Fig. 13. Adaptive !nite element meshes and the crack growth trajectory for the single edge cracked plate with three holes under mixed-mode loading (Case I).
1768
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
Fig. 14. Adaptive !nite element meshes and the crack growth trajectory for the single edge cracked plate with three holes under mixed-mode loading (Case II).
Fig. 15. Problem statement and the !nal mesh of the initial crack for the single edge crack plate.
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1769
Fig. 16. Adaptive !nite element meshes and the crack growth trajectory for the single edge cracked plate under mixed-mode loading.
middle hole. For the second case with the initial crack length, a = 1:0 and location b = 4:0 units, respectively, the crack propagates toward the middle hole (Fig. 14). In both cases, the predicted crack growth trajectories very closely resemble the experimental results. Fig. 15 shows the problem statement and !nal mesh of a single edge cracked plate under mixed-mode loading. This plate is !xed at the bottom edge and is subjected to far-!eld shear stress 3 = 1 unit along the top edge. The plate has an initial crack length a = 3:5 units. The modulus of elasticity and the Poisson’s ratio are 30 × 106 units and 0.25, respectively. The plane strain condition is assumed in the analysis. The computed stress intensity factors KI and KII from the adaptive mesh are 34.10 and 4.52 comparing to the reference values of 34.00 and 4.55, respectively. Fig. 16 shows the predicted crack trajectory that closely resembles the numerical solution in Ref. [27]. 5. Conclusions The adaptive !nite element method using Delaunay triangulation for crack propagation analysis was presented. The concept of Delaunay triangulation for two dimensional mesh construction was explained. The mesh generation procedure with automatic node creation and mesh smoothing, which was implemented by object-oriented programming, were presented in details. The technique was combined with the !nite element method for analyzing crack problems under single and mixed mode loadings. The isoparametric six-node triangular elements, with mid-side nodes displaced from their nominal positions to quarter points of the crack tip, were employed to form up a circular zone surrounding the tip in order to better capture the stress !eld. The solution accuracy was further improved by incorporating an adaptive remeshing technique into the Delaunay triangulation algorithm. The adaptive remeshing technique places small elements around the crack tips and in regions with large changes of stress gradients. At the same time, larger elements are generated in other regions to minimize the total number of unknowns and the
1770
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
computational time. Several examples were employed to evaluate the accuracy of the combined Delaunay triangulation, the !nite element method, and the adaptive remeshing technique by predicting the stress intensity factors for several geometries with di4erent loadings, as well as capturing the crack propagation trajectory. These examples demonstrated the capability of the combined adaptive Delaunay triangulation with the !nite element method for solving crack propagation problems e4ectively. Acknowledgements The authors are pleased to acknowledge the Thailand Research Fund (TRF) for supporting this research work. References [1] A. Bowyer, Computing Dirichlet tessellations, Comp. J. 24 (2) (1981) 162–166. [2] D.F. Watson, Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes, Comp. J. 24 (2) (1981) 167–172. [3] N.P. Weatherill, O. Hassan, EKcient three-dimension Delaunay triangulation with automatic point creation and imposed boundary constraints, Int. J. Numer. Methods Eng. 37 (12) (1994) 2005–2039. [4] B.K. Karamete, T. Tokdemir, M. Ger, Unstructured grid generation and a simple triangulation algorithm for arbitrary 2-D geometries using object oriented programming, Int. J. Numer. Methods Eng. 40 (2) (1997) 251–268. [5] R.S. Barsoum, On the use of isoparametric !nite elements in linear fracture mechanics, Int. J. Numer. Methods Eng. 10 (1) (1976) 25–37. [6] R.S. Barsoum, Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements, Int. J. Numer. Methods Eng. 11 (1) (1977) 85–98. [7] S.K. Chan, I.S. Tuba, W.K. Wilson, On the !nite element method in linear fracture mechanics, Eng. Fract. Mech. 2 (1) (1970) 1–17. [8] T.D. de Araujo, D. Roehl, T.N. Bittencourt, L.F. Martha, Adaptive simulation of elastic-plastic fracture process, new trends and applications, in: S. Idelsohn, E. Onate, E. Dvorkin (Eds.), Fourth World Congress on Computational Mechanics, Buenos Aires, Argentina, 1998. [9] D.O. Potyondy, P.A. Wawrzynek, A.R. Ingra4ea, An algorithm to generate quadrilateral or triangular element surface meshes in arbitrary domains with applications to crack propagation, Int. J. Numer. Methods Eng. 38 (10) (1995) 2677–2701. [10] T.N. Bittencourt, P.A. Wawrzynek, A.R. Ingra4ea, J.L. Sousa, Quasi-automatic simulation of crack propagation for 2D LEFM problems, Eng. Fract. Mech. 55 (2) (1996) 321–334. [11] D.M. Parks, A sti4ness derivative !nite element technique for determination of crack tip stress intensity factors, Int. J. Fract. 10 (4) (1974) 487–502. [12] B. Moran, C.F. Shih, A general treatment of crack tip contour integrals, Eng. Fract. Mech. 35 (1987) 295–310. [13] G.V. Guinea, J. Planas, M. Elices, KI evaluation by the displacement extrapolation technique, Eng. Fract. Mech. 66 (3) (2000) 243–255. [14] F. Erdogan, G.C. Sih, On the crack extension in plates under plane loading and transverse shear, J. Basic Eng. 85 (1963) 519–527. [15] R.J. Nuismer, An energy release rate criterion for mixed mode fracture, Int. J. Fract. 11 (2) (1975) 245–250. [16] G.C. Sih, Strain-energy-density factor applied to mixed-mode crack problems, Int. J. Fract. 10 (1974) 305–321. [17] O.C. Zienkiewicz, R.L. Taylor, Finite Element Method, 5th Edition, Butterworth-Heinemann, Woburn, 2000. [18] C.L. Lawson, Software for C1 surface interpolation, in: J.R. Rice (Ed.), Mathematical Software III, Academic Press, New York, 1977.
S. Phongthanapanich, P. Dechaumphai / Finite Elements in Analysis and Design 40 (2004) 1753 – 1771
1771
[19] S.W. Sloan, A fast algorithm for generating constrained Delaunay triangulations, Comp. Struct. 47 (3) (1993) 441–450. [20] W.H. Frey, Selective re!nement: a new strategy for automatic node placement in graded triangular meshes, Int. J. Numer. Methods Eng. 24 (11) (1987) 2183–2200. [21] H. Borouchaki, P.L. George, Aspects of 2-D Delaunay mesh generation, Int. J. Numer. Methods Eng. 40 (11) (1997) 1957–1975. [22] J. Peraire, M. Vahdati, K. Morgan, O.C. Zienkiewicz, Adaptive remeshing for compressible Oow computations, J. Comput. Phys. 72 (2) (1987) 449–466. [23] P. Dechaumphai, Improvement of plane stress solutions using adaptive !nite elements, J. Chin. Inst. Eng. 19 (3) (1996) 375–380. [24] M. Isida, E4ect of width and length on stress intensity factors of internally cracked plates under various boundary conditions, Int. J. Fract. 7 (1971) 301–316. [25] Y. Murakami, Stress Intensity Factors Handbook, Pergamon Press, Oxford, 1987. [26] ASTM 1996, Annual Book of ASTM Standards, American Society for Testing and Materials, Philadelphia, 1996. [27] B.N. Rao, S. Rahman, An eKcient meshless method for fracture analysis of cracks, Comput. Mech. 26 (4) (2000) 398–408.