Adaptation of quadtree meshes in the scaled boundary finite element method for crack propagation modelling

Accepted Manuscript Adaptation of Quadtree Meshes in the Scaled Boundary Finite Element Method for Crack Propagation Modelling E.T. Ooi, H. Man, S. Natarajan, C. Song PII: DOI: Reference:

S0013-7944(15)00363-X http://dx.doi.org/10.1016/j.engfracmech.2015.06.083 EFM 4742

To appear in:

Engineering Fracture Mechanics

Received Date: Revised Date: Accepted Date:

9 September 2014 22 June 2015 27 June 2015

Please cite this article as: Ooi, E.T., Man, H., Natarajan, S., Song, C., Adaptation of Quadtree Meshes in the Scaled Boundary Finite Element Method for Crack Propagation Modelling, Engineering Fracture Mechanics (2015), doi: http://dx.doi.org/10.1016/j.engfracmech.2015.06.083

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ADAPTATION OF QUADTREE MESHES IN THE SCALED BOUNDARY FINITE ELEMENT METHOD FOR CRACK PROPAGATION MODELLING E. T. Ooi*, H. Man†, S. Natarajan‡, C. Song† June 30, 2015 *School of Science, Information Technology and Engineering, Federation University, Ballarat, VIC 3353, Australia †School of Civil and Environmental Engineering, The University of New South Wales, Sydney, NSW 2052, Australia ‡Department of Mechanical Engineering, Indian Institute of Technology, Madras, Chennai, 600036, India

Abstract A crack propagation modelling technique combining the scaled boundary finite element method and quadtree meshes is developed. This technique automatically satisfies the compatibility requirement between adjacent quadtree cells irrespective of the presence of hanging nodes. The quadtree structure facilitates efficient data storage and rapid computations. Only a single cell is required to accurately model the stress field near crack tips. Crack growth is modelled by splitting the cells in the mesh into two. The resulting polygons are directly modelled by the scaled boundary formulation with minimal changes to the mesh. Four numerical examples demonstrate the salient features of the technique. Keywords: scaled boundary finite element method; quadtree; fracture; crack propagation

1

Introduction

The quadtree decomposition [1] is a hierarchical-type data structure, in which each parent is recursively divided into four children. Quadtree algorithms derived from this quarternary tree structure are particularly efficient for data 1

storage and retrieval. The efficiency of these algorithms make them particularly popular in applications that demand high speed and efficiency such as computer graphics [2, 3] and image processing [4]. In computational mechanics, quadtree algorithms are usually used in large scale simulations typical in the modelling of earthquake and ground motions [5], flood [6] and tsunami [7]. In solid mechanics, where the finite element method is most widely used, the quadtree decomposition is an intermediate step that is usually employed in the automatic mesh generation of triangular meshes [8, 9] and quadrilateral meshes [10, 11]. Mesh generators employing quadtree algorithms are fast, efficient and are capable of achieving rapid and smooth transitions of element sizes between regions of mesh refinement. Adaptation of quadtree meshes for direct computational analyses within the framework of the finite element method, however, is not widespread. The primary reason for this is the presence of hanging nodes (shown as the filled circles in Fig. 1) between two adjacent elements of different sizes. The presence of hanging nodes destroys the displacement compatibility between the adjacent elements, prohibiting them to be directly used with quadtree meshes. Several methods have been proposed in the literature to resolve the displacement incompatibility introduced by the hanging nodes so that quadtree meshes can be used in finite element computations.These include constraining the displacements of the hanging nodes through Lagrange multipliers or penalty methods e.g. [5], constraining the hanging nodes to the corner nodes and adding temporary triangular or rectangular elements [12]. These approaches can only be used when the quadtree mesh satisfy the 2 : 1 rule i.e. the maximum difference between two adjacent elements is less than one [8]. Constraint equations have to be properly constructed to treat the displacement incompatibility whenn the 2 : 1 rule is not satisfied e.g. Ainsworth et al. [13]. An alternate approach that enables quadtree meshes to be directly used in finite element computations is by developing special transitional finite elements for the cells in a quadtree mesh that contain hanging node(s) e.g. [14, 15, 16]. Gupta [14] derived the shape functions for quadrilateral finite elements containing hanging nodes that are located at the midpoints of the element edges by modifying the shape functions of the corner nodes of the standard quadrilateral elements. Transition finite elements can also be derived using the assumed stress hybrid method of Pian [17]. Within this approach, several element designs are possible; as outlined by Lo et al. [15] and Sze and Wu [16]. Alternative approaches that can be used to construct transition elements include using incompatible modes [18]. The use of quadtree meshes within the framework of the extended finite element method [19] to model problems involving fracture and discontinuous surfaces has also been reported in the literature. Fries et al. [20] reported two approaches to adapt quadtree meshes into the extended finite element method. The first approach uses the transition elements developed by Gupta [14]. The second approach constrains the displacements of the hanging nodes to be the average of the neighbouring corner nodes and is similar to the approach adopted by Bielak et al. [5]. 2

Sukumar and Tabarraei [21] outlined a technique to develop conforming polygon elements having arbitrary number of sides. The shape functions of the polygon elements can be derived, for example, from Wachspress interpolants, mean value coordinates and natural neighbour interpolants. By employing polygons that can have arbitrary number of sides, the cells in a quadtree mesh containing hanging nodes can be treated simply as a polygon. Application of these polygon elements within the context of quadtree meshes was reported by Tabarraei and Sukumar [22] for stress analysis and by Folch et al. [23] for phase induction computations. Tabarraei and Sukumar [24] also developed their polygon-based elements for fracture analysis by coupling with the extended finite element method. The scaled boundary finite element method (SBFEM) is a semi-analytical method developed by Song and Wolf [25] to solve boundary value problems. It is known to be very efficient when applied to modelling fracture e.g. [26, 27]. When applied to fracture problems, stress singularities of any kind can be accurately modelled as their corresponding stress modes are analytically represented in the SBFEM solution of stresses [26]. The SBFEM is also sufficiently flexible such that it can be formulated on polygons with arbitrary number of sides [28, 29]. This feature makes it particularly suitable to be adapted for computations with quadtree meshes via the approach of Tabarraei and Sukumar [24]. Each cell in a quadtree mesh, regardless of the presence of hanging nodes, is treated as a generic polygon. This enables the structure of the quadtree to be exploited for efficient computations. The ability to assume any number of sides also enables the SBFEM to better discretize curved boundaries. The cells that are trimmed by the boundaries are directly treated as polygons. No excessive refinement using square cells is required. Motivated by these advantages, an approach that combines a quadtree mesh with the SBFEM to model crack propagation is developed in this study. When a crack propagates, it is only necessary to split each cell cut by the crack into two. These newly formed cells can be regarded as polygons having arbitrary number of sides. They are directly modelled by the SBFEM and no further remeshing or enrichment is necessary. Changes to the mesh are minimal. This, coupled with the computational efficiency derived from the quadtree data structure, leads to an elegant technique to model crack propagation. This paper is structured as follows. In Section 2 the fundamentals of the SBFEM are briefly reviewed. Section 3 discusses the adaptation of quadtree meshes into the SBFEM and its efficient implementation. In Section 4, the application of the SBFEM in fracture analysis is reviewed briefly. This is followed by a description of the crack propagation methodology that is developed in this study. In Section 5, the salient features of the proposed technique are demonstrated using four numerical examples and discussed in detail. The major conclusions of the developed technique are summarised in Section 6.

3

2

Scaled Boundary Finite Element Method

In this section, the fundamentals of the SBFEM will be briefly reviewed. A complete discussion on this subject has been reported in numerous instances in the literature [e.g. 25, 30]. The SBFEM can be implemented on polygons with arbitrary number of sides [28]. The geometry of a polygon is required to satisfy a so-called “direct visibility criterion”: there exists a point within the polygon where the entire boundary is directly visible from it. This point is called the scaling centre. It is to be noted that this geometric requirement is also referred to as the “star convexity” criterion. The scaled boundary finite element coordinate system of a generic polygon is shown in Fig. 2. In each polygon, the scaled boundary coordinate system (ξ, η) is established. ξ is the radial coordinate and is 0 at the scaling centre and 1 at the polygon boundary. η is the local coordinate defined for each one-dimensional finite element that discretizes the polygon boundary. The equilibrium condition in each polygon can be derived for example, starting from the method of weighted residuals, which result in the following first order differential equation { } { } u(ξ) u(ξ) ξ =−Z (1) q(ξ) ,ξ q(ξ) where u(ξ) are displacement functions, q(ξ) is the equivalent internal force vector and Z is a Hamiltonian matrix with [ ] T E−1 −E−1 0 E1 0 Z= (2) −1 T E1 E−1 0 E1 − E2 −E1 E0 The coefficient matrices E0 , E1 and E2 are computed on each line element discretising the polygon boundary and assembled in a manner similar to the finite element method [30]. The matrix Z can be decoupled into pairs of eigenvalues (−λi , λi ) using an eigenvalue decomposition or a block-diagonal Schur decomposition [30] as ZΨ =ΨS

(3)

where Ψ is a transformation matrix with independent column vectors. In the case where an eigenvalue decomposition is adopted, S is a diagonal matrix with its diagonal corresponding to the eigenvalues. When a block diagonal Schur decomposition is adopted, S is in real Schur form and is a quasi-upper triangular matrix consisting of 2 × 2 blocks on the diagonal corresponding to complex conjugate eigenvalues or 1 × 1 blocks corresponding to real eigenvalues. The eigenvalues in S are sorted in ascending order of their real parts. Ψ and S can be partitioned conformably as [ ] Ψun Ψup Ψ= (4) Ψqn Ψqp [ ] Sn S= (5) Sp 4

In Eq. (5), Sn and Sp contain the eigenvalues with negative and positive real parts resulting from an eigenvalue or a block-diagonal Schur decomposition of the Hamiltonian matrix Z. The stiffness matrix of each polygon Kp can be evaluated as Kp =Ψ−1 qn Ψun

(6)

For bounded polygons considered in this paper, the displacement field of a sector covered by a line element on the polygon boundary is [30] u(ξ, η) =N(η)Ψun ξ −Sn c

(7)

where N(η) is the shape function matrix and c are the integration constants that depend on the boundary conditions. They are determined from the nodal displacements on the polygon boundary ub = u(ξ = 1) as c =Ψ−1 un ub

(8)

The stress field in a line sector covered by a line element on the polygon boundary is σ(ξ, η) =Ψσ (η)ξ −Sn −I c where I is an identity matrix and Ψσ (η) = [ Ψσxx (η) Ψσyy (η) are the stress modes Ψσ (η) =D(−B1 (η)Ψun Sn + B2 (η)Ψun )

(9) Ψτxy (η) ]T (10)

In Eq. (10), D is the material constitutive matrix, and B1 (η) and B2 (η) are the standard scaled boundary finite element strain-displacement matrices [25].

3

Adaptation of Quadtree Meshes in the Scaled Boundary Finite Element Method

A cell in a quadtree mesh can be directly modelled by the SBFEM. Each cell is treated as a polygon having arbitrary number of sides similar to the approach of Tabarraei and Sukumar [22]. Unlike transitional finite elements, no distinction is made between a corner node and a hanging node. All the cells in the quadtree mesh are directly modelled by a unified formulation within the framework of the SBFEM. From the point-of-view of finite element mesh generation algorithms, this approach is particularly appealing given that a quadtree decomposition is efficient and usually employed prior to triangularisation of a mesh. The ability to use the quadtree mesh directly eliminates this triangularisation process and reduces the computational resources for mesh generation. The SBFEM can be applied to quadtree meshes without adhering to the 2 : 1 rule. However, the 2 : 1 rule offers some appealing features that facilitate more efficient computations. When the 2 : 1 rule is enforced, the cells in the mesh 5

are limited to only the patterns shown in Fig. 3. Considering all orthogonal rotations of each of the cells, the total number of different cells that can be present in the quadtree mesh is 24. This can be reduced to 16 for isotropic homogeneous materials considering that rotation does not have any effect on the 4-node and 8-node cells. For the 6-node cells of the first type (the top one in Fig. 3), only two rotations are necessary. If the mesh is purely made up of cells of the type shown in Fig. 3, only the stiffness matrices of these 16 cells need to be computed in a linear elastic analysis. These can be pre-calculated and stored for quick data retrieval. These 16 square cells shall be refered to as the “master cells” throughout the remainder of this manuscript.

3.1

Hybrid quadtree mesh generation

In the present study, a hybrid quadtree mesh is adopted for mesh generation. The hybrid quadtree mesh adopted in this study is predominantly made up of the master cells shown in 3. Irregular polygons are used to better fit the geometry of the boundary. This section describes in brief the procedure for the generation of the hybrid quadtree mesh. Suppose, it is desired to generate a generic domain as that shown in Fig. 3.1a. A background rectangular domain ΩR is first defined as is shown in Fig. 3.1b. It is discretized entirely by square cells. The geometry of a generic domain, Ω ∈ ΩR , to be analysed is defined using the signed distance function dΩ . Within Ω, the signed distance function associated with a point x ∈ R2 , which is a subset of R2 is defined as dΩ (x) =sΩ (x) min ||x − y|| y∈∂Ω

(11)

where ||x − y|| is the Euclidean norm in R2 with y ∈ ∂Ω. The sign function sΩ (x) is equal to −1 when x lies inside the domain and is equal to 1 otherwise. The boundary of the domain ∂Ω and is defined by the zero level set of dΩ . Complex geometries can be derived by combining various distance functions of simple geometries through Boolean operators [31]. In Fig. 3.1b, the signed distance function dΩ used to define the boundary are indicated in bold. Along the boundary ∂Ω, an arbitrary number of seed points ngeo can be defined. The seed points ngeo is input manually and divides the boundary ∂Ω into ngeo − 1 segments and is similar to the seed function in the commercial finite element software ABAQUS. ngeo therefore, allows the user to control the desired resolution of the boundary ∂Ω. The mesh density along the boundary can be controlled by specifying the maximum number of seed points allowable in a square cell nall . The value for nall is input manually. Depending on the specified value of ngeo , the square cells in Fig. 3.1b may contain a number of seed points, ncell . The number of seed points inside a cell is compared with nall to ensure that ncell < nall . Any square cell that does not meet this criterion is then refined whilst adhering to the 2 : 1 rule. In Fig. 3.1c, the cells around the boundary of the smaller circle is refined to meet the abovementioned criterion. Both ngeo and nall controls the final resolution of the quadtree discretization. A 6

high value of ngeo together with a small value of nall will lead to a dense mesh along ∂Ω and vice versa. The hybrid mesh of Ω is then obtained by trimming the square cells in ΩR that are cut by ∂Ω. This process may result in irregular shaped quadrilaterals or sometimes, polygons of more than 4 sides. Hereafter, any quadrilaterals or polygons that fall into these categories are referred to as irregular cells. The irregular cells can be directly modelled by the SBFEM without additional subdivision. As a result, the cells are square only in regions of the mesh that are not in direct contact with the boundary ∂Ω. The final hybrid quadtree mesh of the geometry shown in Fig. 3.1a obtained after this trimming procedure is shown in Fig. 3.1d.

3.2

Comparison of a hybrid quadtree and pure quadtree mesh

The hybrid quadtree mesh offers some advantages over a purely quadtree based mesh in discretising non-regular geometries. The primary reason for this is due to the fact that a purely quadtree mesh is composed of only squares having horizontal and vertical lines. When modelling curved boundaries, which is a generic feature in real life structures and components, the regions near these boundaries require further mesh refinement in order to reduce the discretization error. This may result in an overly refined mesh in these regions. Despite this, the boundary may still not be smooth and may result in unrealistically high stresses. Fig. 3.2a and Fig. 3.2b show two typical geometries that are discretized with the hybrid quadtree meshes. Their corresponding pure quadtree meshes are also shown for comparison. In the hybrid quadtree meshes shown in Fig. 3.2, the master cells shown in Fig. 3 are plotted in grey. The black lines in the mesh, which are predominantly visible along the external and internal boundaries of the domain are made up of irregular cells. A visual comparison of both meshes reveals that the boundary of the domain is more accurately modelled by a hybrid quadtree mesh compared with a purely quadtree mesh. Moreover, the mesh density required by a hybrid quadtree mesh is significantly reduced compared with a mesh of purely quadtrees. The pure quadtree mesh shown in Fig. 3.2a is made up of 4501 square cells. Its corresponding hybrid quadtree mesh, however, require only 308 polygons, out of which, 170 are the master cells shown in Fig. 3 and 138 of these are irregular cells. Fig. 3.2 compares the quality of the discretization near the boundary of both the hybrid quadtree and the pure quadtree meshes for Geometry 1 in Fig. 3.2a. It is interesting to note that the boundary is still not accurately represented by the square cells despite the high mesh density. In the hybrid quadtree mesh, the boundary is more accurately represented due to the flexibility of the SBFEM to assume polygons of any number of sides. At this juncture, it is worth mentioning that if desired, the boundaries can be discretized with more elements or high-order elements to better fit the curved boundaries. The construction of high-order shape functions for polygons with 7

arbitrary number of sides within the framework of the SBFEM is rather straightforward and does not require, for example, the solution of an optimisation problem when the maximum entropy shape functions are used [32]. During numerical execution of the algorithm, only the stiffness matrices of the polygons having arbitrary number of sides in addition to the 16 master cells have to be computed. This reduces the number of stiffness matrix computations in the mesh drastically. For the quadtree-polygon mesh shown in Fig. 3.2a, the amount of stiffness matrix computations is reduced by 50%. Savings in computational time can be derived from the reduced time required for the stiffness matrix computations. The extent of reduction in the number stiffness matrix computations depends on the geometry that is being modelled. For example, in the hybrid quadtree mesh shown in Fig. 3.2b, there are 574 cells, out of which, 266 are irregular. The total reduction in the number of stiffness matrix computations for this mesh is 50.9%. In order to obtain a better grasp of the amount of savings in computational time due to the reduction of stiffness matrix computations, the computational times required to calculate and assemble the stiffness matrices by (i) pre-computing the stiffness matrices of the master cells in Fig. 3 and (ii) computing the stiffness matrices of all the cells on the fly are investigated. For this purpose, the hybrid quadtree meshes shown in Fig. 3.2a and Fig. 3.2b will be considered. Table 3.2 shows the comparison of the computational times required for both approaches. All the computations are performed using an Intel(R) Core(TM) i7-4500U CPU @ 1.8GHz laptop. The percentage of savings in computational times in both meshes amounts to approximately 47%.

Table 1. Comparison of computational times for hybrid quadtree meshes.

On the fly Geometry 1 Geometry 2

4 4.1

0.9223 1.1495

Computational Times (s) Pre-computing basic cells Master cells Irregular Total 0.0382 0.3685 0.4067 0.0360 0.5428 0.5788

% Savings 44.09 49.64

Fracture Analysis Computation of stress intensity factors

A crack is sufficiently modelled by the SBFEM using a single polygon as shown in Fig. 4.1. The crack surfaces are not discretized. The line elements discretising the polygon boundary do not form a close loop. The stress solution (Eq. (9)) in a crack polygon contains two eigenvalues, λi , i = 1, 2 satisfying −1 < Re(λi ) < 0 in the matrix Sn . By inspection of Eq. (9), it is found that these eigenvalues lead to a singular stress field as ξ → 0. This enables the SBFEM to accurately model 8

the stress singularity at the crack tip without the need for local mesh refinement and local enrichment with asymptotic functions. This feature of the SBFEM makes it particularly appealing when applied to fracture problems [26, 28, 27] Using the two modes contributed by the eigenvalues λi , a singular stress field σ (s) (ξ, η) can be obtained from Eq. (9) as −Sn σ (s) (ξ, η) =Ψ(s) σ (η)ξ

(s)

−I (s)

(12)

c

(s)

where Sn is a block-diagonal matrix containing the two eigenvalues λi in the case where the Schur decomposition is used. When an eigenvalue decompsotion (s) is applied to Eq. (3), Sn is a diagonal matrix containing the two eigenvalues (s) λi . c(s) are the integration constants corresponding to Sn . The singular stress (s) (s) (s) (s) T mode Ψσ (η) = [ Ψσxx (η) Ψσyy (η) Ψτxy (η) ] in Eq. (12) is (s) (s) (s) Ψ(s) σ (η) =D(−B1 (η)Ψun Sn + B2 (η)Ψun )

(13) (s)

where Ψ(s) un are the modal displacements in Ψun corresponding to Sn . The stress intensity factors are directly evaluated from their definitions. For a crack that is aligned with the Cartesian coordinate system as shown in Fig. 4.1, the stress intensity factors are defined as { } { √ } KI 2πrσyy |θ=0 √ = lim (14) KII r→0 2πrτxy |θ=0 Substituting the stress components in Eq. (12) into Eq. (14) and using the relation ξ = r/LA , where LA is the distance from the scaling centre to the boundary along the direction of the crack at θ = 0 (see Fig. 4.1) leads to { } { } √ (s) Ψ(s) KI σyy (η|θ=0 )c = 2πLA (15) (s) KII Ψ(s) τxy (η|θ=0 )c

4.2

Modelling of a crack in a quadtree-polygon mesh

In order to use the hybrid quadtree mesh fracture analysis with the SBFEM, the mesh in the crack tip region has to be modified so that it takes the form of a cracked cell as is shown in Fig. 4.1. The procedure required for this is depicted in Fig. 4.2. Consider for instance, the mesh in the region near a crack as shown in Fig. 4.2a. The crack path is shown in bold and the crack tip is indicated by the cross. To generate the cracked cell, it is first necessary to identify the cells surrounding the crack tip. These cells are depicted Fig. 4.2b by the shaded area. Once identified, these cells are refined so that they have the same size as the cell containing the crack tip. This results in the mesh shown in Fig. 4.2b and avoids the case where a crack tip is too close to the edges of a cell, which may affect the mesh quality and the solution accuracy [28].

9

Next, the cells around the crack tip are merged to form a cracked cell. The intersection point between the edge of the resulting cell and the crack is computed to define the crack mouth. The other cells cut by the crack path are split into two by creating a pair of vertices and edges at the intersection point. This procedure is similar to that used in splitting the cells around the boundary of the domain described in Section 3.1 and will generate irregular cells along the crack path. However, this does not pose a problem in the SBFEM as it makes no distinction between a square cell or a polygon. In the crack polygon, the creation of a pair of vertices at the crack mouth results in a polygon that does not form a close loop as shown in Fig. 4.1. This is also directly modelled by the SBFEM. The final mesh after merging the cells are shown in Fig. 4.2.

4.3

Crack propagation direction

The evolution of the crack geometry is determined from the stress intensity factors computed in Eq. (15). A modified expression of the maximum circumferential stress criterion outlined by Sukumar and Prevost [33] is used to determine the crack propagation direction θp ( ) −2KII (θc )/KI (θc ) −1 √ θp =2 tan (16) 1 + 1 + 8(KII (θc )/KI (θc ))2 where θc is the crack angle.

4.4

Automatic remeshing algorithm

Once the the crack propagation direction has been determined from Eq. (16), the mesh is updated to accommodate the new crack propagation length |f |. The remeshing process is shown in Fig. 4.4. It involves the following three major steps: (i) quadtree refinement at the location of the old and new crack tips; (ii) regenerating a new crack polygon at the new crack tip and; (iii) splitting the cells cut by the crack into two. The remeshing process commences by first refining the mesh around the old crack tip. The minimum size of the cells adjacent to the old crack polygon hmin is identified. The old cracked cell is refined so that it is composed of cells having the same size as hmin . It is to be noted that the the 2 : 1 rule is adhered to during this process. Next, the region around the new crack tip which requires mesh refinement is identified. This region is determined by plotting a circle Ωcir of radius rcir < 0.5|f | centered at the location of the new crack tip as shown in Fig. 4.4a. Any cells that are in contact with this circle will be refined. The extent of mesh refinement required is controlled by specifying a number of “seed points”, ncir , along the boundary of Ωcir . The value of ncir can be determined prior to an analysis. It suffices to use a value of ncir between 10 ≤ ncir ≤ 30. A quadtree decomposition is then recursively applied to these cells until each cell contains only one seed point. The 2 : 1 rule is then enforced. The region affected by 10

remeshing may therefore be much larger than Ωcir . Fig. 4.4b shows the crack tip mesh in which, the region encompassed by Ωcir has been refined (in grey). It is to be noted that the mesh refinement ensures that there are enough nodes on the boundary of the cracked cell for accurate computation of the stress field in the vicinity of the crack tip . It also ensures that the size of the crack polygon generated in the subsequent step is of reasonable size and does not interfere with the surrounding boundaries. This will be crucial when there are multiple cracks that are propagating close to each other. A cracked cell can then be generated at the location of the new crack tip using the procedure outlined in Section 4.2. This leads to the mesh shown in Fig. 4.4c.

5

Numerical Examples

Four crack propagation benchmarks are modelled to validate the developed technique. In the sections that follow, each edge of a cell or a polygon, which is used to discretize the computational domain, is modelled using one linear element. For a cracked cell, each edge is discretized using two linear elements. It is worth noting that discretising each edge on a cracked cell into two elements does not create irregular cells around it. The resulting cells are of the 5-node type master cell shown in Fig. 3. The simulations were executed on an Intel(R) Core(TM) i74500U CPU @ 1.8GHz laptop. An in-house serial code implemented in Matlab is used for mesh generation, solution and remeshing.

5.1

Crack path deflection in a plate due to a hole

The 15mm×20mm plate with an interior hole of radius r = 3.45mm and a crack of length 2.75mm on its vertical left edge as shown in Fig 5.1a is considered. The plate is subjected to a uniform tension σ on both its horizontal edges. The material properties of the plate are as follows: Young’s modulus E = 98GPa and Poisson’s ratio ν = 0.3. Plane stress conditions are assumed. This crack propagation benchmark was used by Rashid [34] to validate the arbitrary local mesh replacement technique that was proposed to model crack propagation. Bouchard et al. [35] used the same problem to validate his finite element code that combines an automatic local remeshing algorithm. 5.1.1

Crack propagation modelling

Fig. 5.1b shows the initial hybrid quadtree mesh of the plate. The mesh has 249 cells, of which 88 of these are irregular. The total number of nodes in the mesh is 384. Crack propagation is modelled by specifying an incremental crack length of |f | = 1mm. Fig. 5.1.1b compares the predicted crack trajectories of the proposed approach with that reported by Rashid [34] and by Bouchard et al. [35]. The crack paths predicted by the three methods agree very well with each other.

11

Due to the presence of the cavity, the crack path initially turns towards the hole. As the crack tip reaches the middle of the plate, the crack reorients itself and propagates horizontally until it reaches the vertical edge on the opposing side. The final hybrid quadtree mesh is shown in Fig. 5.1.1b. The final mesh has 390 cells of which 164 are irregular. 5.1.2

Investigation of effect of mesh density

The benchmark described in Section 5.1 will be used to investigate the effects of mesh density on the predicted crack propagation paths. For this purpose, the meshes shown in Fig. 5.1.2 are considered in addition to the coarse mesh in Fig. 5.1b. The medium mesh (Fig. 5.1.2a) has 613 cells and 881 nodes out of which, 197 are irregular cells. The fine mesh (Fig. 5.1.2b) has 1098 cells and 1524 nodes out of which, 327 are irregular cells. The crack propagation simulation is repeated using these meshes with the same crack incremental length |f | = 1mm. Fig. 5.1.2 compares the predicted crack paths using the three meshes: coarse, medium and fine. It is evident from Fig. 5.1.2 that the mesh density does not have any appreciable effect on the predicted crack paths. This observation is attributed to the accuracy of the SBFEM in modelling the stress singularity at the crack tip. The coarse mesh is sufficient to simulate the crack propagation in the plate.

5.2

Cracked beam with three holes

The cracked beam with three holes subjected to three-point bending shown in Fig. 5.2 is considered. The load P = 4.45N. This problem is a widely used numerical benchmark for crack propagation problems [e.g. 36, 37]. The beam is made of polymethylcrylate with the following material properties: Young’s modulus E = 29 × 106 kPa and Poisson’s ratio ν = 0.3. Plane stress conditions are assumed. Two test cases, I and II, having different values of a and b as shown in Fig 5.2 are considered. Fig. 5.2a and Fig. 5.2b show the initial scaled boundary finite element meshes of the beam for Cases I and II, respectively. The mesh for Case I has 320 cells of which 242 are the master cells shown in Fig. 3 and 78 are irregular cells. The total number of nodes for this mesh is 438. In Case II, the mesh has 476 cells of which 354 are master cells and 122 are irregular cells. The total number of nodes for this mesh is 661. In both Cases I and II, a crack incremental length of |f | = 0.5mm is specified for the simulation. Fig. 5.2a and Fig. 5.2b show the final hybrid quadtree meshes of the beam for Cases I and II, respectively. For Case I, the final mesh has 664 cells of which, 215 are irregular. For Case II, the final mesh has 829 cells of which, 264 are irregular. The difference in the crack length and location, as outlined in Fig. 5.2 for both the cases considered, is observed to affect the crack propagation path. In Case I, the crack initially curves towards the bottom hole. As it approaches the bottom hole, the crack turns towards the the middle hole. The crack path once again deflects towards the loading point as it approaches the

12

middle hole. In Case II, the crack propagation is a curved line evolving towards the middle hole. The crack trajectories in both cases obtained by the SBFEM are reminiscence of the experimental results reported by Bittencourt et al. [36] and the finite element simulation of Azocar et al. [37]. A comparison of the results is shown in Fig. 5.2a and Fig. 5.1.2b. Fig. 5.2 shows a 20mm × 10mm plate with two circular cavities of diameter d = 4mm and two edge cracks of length 1mm along its opposing vertical edges. The material properties of the plate are: Young’s modulus E = 200GPa and Poisson’s ratio ν = 0.3. Plane stress conditions are assumed. The propagation of the cracks in the plate subjected to uniformly distributed vertical displacements along the horizontal edges is modelled. The results of the present method is compared with the finite element solution of Bouchard et al. [38]. The plate is discretised by the hybrid quadtree mesh shown in Fig. 5.2b. There are 371 cells in the mesh, of which 132 are irregular cells. The mesh has 588 nodes. Crack propagation is modelled specifying a crack incremental length of |f | = 1mm. Fig. 5.2a compares the crack path predicted by the present method with the finite element solution of Bouchard et al. [38] and the adaptive finite element method with superconvergent patch recovery technique developed by Khoei et al. [39]. The crack propagation of both cracks are symmetric. They initially curve towards the adjacent circular cavities and reorient almost horizontally and propagate towards each other until they reach the middle of the plate. The crack paths then curve towards the circular cavity at the opposite end of the plate. Fig. 5.2b shows the final mesh of the plate. Compared to the initial mesh, the number of cells have increased to 770 of which, 365 are irregular.

5.3

Propagation of two cracks emanating from holes

Fig. 5.3a shows a square plate of dimensions 5 × 5 having two internal circular cavities with radius r = 5/64. A crack of length a = 0.1 emanates from each hole at an angle θ = π/4. The plate is subjected to a uniform tension σ = 2 on both its horizontal edges. The material properties of the plate are: Young’s modulus E = 200 and Poisson’s ratio ν = 0.3. Plane stress conditions are assumed. A set of consistent units is chosen. Fig. 5.3b shows the initial hybrid quadtree mesh of the plate. The mesh has 454 cells, of which 110 are irregular. There are 667 nodes in the mesh. A crack incremental length of |f | = 0.1 is specified. The final mesh has 950 cells and contains 363 irregular cells. Fig. 5.3a shows a close up view of the mesh in the region surrounding the crack path. Due to the nature of the problem both crack trajectories are symmetric. Both cracks initially propagate horizontally towards holes opposite to them. When the cracks reach the middle of the plate, they curve towards each other. Both cracks then curve away from each other and propagate in an almost horizontal manner as they approach the holes opposite to them. Fig. 5.3b compares the predicted crack path of the present method with the extended finite element method results of Moes et al. [19]. Good agreement is observed. 13

6

Conclusion

A crack propagation modelling technique has been developed by adapting quadtree meshes in combination with the SBFEM. The use of quadtree meshes complements the SBFEM in that the formulation makes no distinction between hanging nodes and corner nodes. Moreover, the SBFEM does not have to adhere to the 2 : 1 rule that is typically enforced by standard numerical approaches adopting quadtree meshes. However, by enforcing the 2 : 1 rule in the quadtree mesh, the cells in the mesh is limited to a few patterns. The stiffness matrices of these cells are pre-calculated and stored for quick retrieval during analysis. This results in huge reduction in the required computational time during analysis. In order to achieve more accurate discretization of the boundary, the proposed technique combines the quadtree meshes with polygons with arbitrary number of sides. This results in a hybrid quadtree mesh. Along the boundary, polygons with arbitrary number of sides, which can be directly modelled by the SBFEM are used so that the boundary is modelled more accurately. Therefore, only the cells that are not in direct contact the boundary satisfy the 2 : 1 rule. An analysis carried out on two generic meshes showed that this hybrid approach is able to reduce the number of required stiffness matrix computations by approximately 50% and the required computational time by approximately 47%. When applied to fracture problems, the scaled boundary finite element formulation enables accurate modelling of the stress singularity at the crack tip. It does not require local mesh refinement or local enrichment with asymptotic functions. Moreover, the stress intensity factors can be computed directly from their definitions. An efficient approach to model crack propagation was developed by exploiting the flexibility of the scaled boundary formulation in assuming polygons of arbitrary shapes, its acuracy in modelling stress singularities at the crack tip and the computational efficiency of the hybrid quadtree mesh. The remeshing process used to accommodate the evolution of the crack trajectories involves application of the quadtree decomposition in the region near the crack tips and splitting the cells cut by the crack into two. The developed technique was successfully validated using three numerical benchmarks. It was found that the predicted crack propagation paths are not affected by the mesh density. This is attributed to accuracy of the stress intensity factors computed by the SBFEM. The developed technique therefore, offers an interesting and viable alternative to the standard approaches reported in the literature to model crack propagation.

14

References References [1] H. Samet, The quadtree and related hierarchical data structures, ACM Computing Surveys 16 (1984) 187–260. [2] H. Samet, R. E. Webber, Hierarchical data structures and algorithms for computer graphics, i. fundamentals, IEEE Computer Graphics and Applications 8 (1988) 48–68. [3] H. Samet, R. E. Webber, Hierarchical data structures and algorithms for computer graphics, ii. applications, IEEE Computer Graphics and Applications 8 (1988) 59–75. [4] H. Samet, Neighbour finding techniques for images represented by quadtrees, Computer Graphics and Image Processing 18 (1982) 37–57. [5] J. Bielak, O. Ghattas, E. J. Kim, Parallel octree-based finite element method for large-scale earthquake ground motion simulation, Computer Modelling in Engineering and Sciences 10 (2005) 99–112. [6] Q. Liang, G. Du, J. Hall, A. Borthwick, Flood inundation modelling with an adaptive quadtree grid shallow water equation solver, Journal of Hydraulic Engineering 134 (2008) 117–132. [7] S. Popinet, Quadtree-adaptive tsunami modelling, Ocean Dynamics 61 (2011) 1261–1285. [8] M. Yerry, M. Shepard, A modified quadtree approach to finite element mesh generation, IEEE Computer Graphics and Applications 3 (1983) 39–46. [9] M. O. Freitas, P. A. Wawrzynek, J. B. Cavalcante-Neto, C. A. Vidal, L. F. Martha, A. R. Ingraffea, A distributed-memory parallel technique for twodimensional mesh generation for arbitrary domains, Advances in Engineering Software 59 (2013) 38–52. [10] X. Liang, M. S. Ebeida, Y. Zhang, Guaranteed-quality all-quadrilateral mesh generation with feature preservation, Computer Methods in Applied Mechanics and Engineering 199 (2010) 2072–2083. [11] F. B. Atalay, S. Ramaswami, D. Xu, Quadrilateral meshes with provable angle bounds, Engineering with Computers 28 (2012) 31–56. [12] P. Krysl, P. Grinspun, E. Schroder, Natural hierarchical refinement for finite element methods, International Journal for Numerical Methods in Engineering 56 (2003) 1109–1124.

15

[13] M. Ainsworth, L. Demkowicz, C. W. Kim, Analysis of the equilibrated residual method for a posteriori error estimation on meshes with hanging nodes., Computer Methods in Applied Mechanics and Engineering 196 (2007) 3493–3507. [14] A. K. Gupta, A finite element for transition from a fine to a coarse grid, International Journal for Numerical Methods in Engineering 12 (1978) 35– 45. [15] S. H. Lo, K. H. Wan, K. Y. Sze, Adaptive refinement analysis using hybridstress transition elements, Computers and Structures 84 (2006) 2212–2230. [16] K. Y. Sze, D. Wu, Transition finite element families for adaptive analysis of axisymmetric elasticity problems, Finite Elements in Analysis and Design 47 (2011) 360–372. [17] T. H. H. Pian, Development of element stiffness matrices by assumed stress distribution, AIAA Journal: The Americal Institute of Aeronautics and Astronautics 2 (1964) 1333–1364. [18] C. K. Choi, E. J. Lee, Nonconforming variable-node axisymmetric solid element, Journal of Engineering Mechanics 130 (2004) 578–588. [19] N. Moes, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering 46 (1999) 131–150. [20] T. Fries, A. Byfut, A. Alizada, K. W. Cheng, A. Schroder, Hanging nodes and xfem, International Journal for Numerical Methods in Engineering 86 (2011) 404–430. [21] N. Sukumar, A. Tabarraei, Conforming polygonal finite elements, International Journal for Numerical Methods in Engineering 61 (2004) 2045–2066. [22] A. Tabarraei, N. Sukumar, Adaptive computations on conforming quadtree meshes, Finite Elements in Analysis and Design 41 (2005) 686–702. [23] J. R. Folch, J. Perez, M. Pineda, R. Puche, Quadtree meshes applied to the finite element computation of phase inductances in an induction machine, in: S. Wiak, A. Krawczyk, I. Dolezel (Eds.), Intelligent Computer Techniques in Applied Electromagnetics, volume 119 of Studies in Computational Intelligence, pp. 117–124. [24] A. Tabarraei, N. Sukumar, Extended finite element method on polygon and quadtree meshes, Computer Methods in Applied Mechanics and Engineering 197 (2008) 45–438. [25] C. Song, J. P. Wolf, The scaled boundary finite-element method - alias consistent infinitesimal finite-element cell method - for elastodynamics, Computer Methods in Applied Mechanics and Engineering 147 (1997) 329–355. 16

[26] C. Song, F. Tin-Loi, W. Gao, A definition and evaluation procedure of generalised stress intensity factors at ccrack and multi-material wedges, Engineering Fracture Mechanics 77 (2010) 2316–2336. [27] W. Becker, S. Goswami, Computation of 3-d stress singularities for multiple ccrack and crack intersections by the scaled boundary finite element method, International Journal of Fracture 175 (2012) 13–25. [28] E. T. Ooi, C. Song, F. Tin-Loi, Polygon scaled boundary finite elements for crack propagation modelling, International Journal for Numerical Methods in Engineering 91 (2012) 319–342. [29] S. Natarajan, E. T. Ooi, I. Chiong, C. Song, Convergence and accuracy of displacement based finite element formulations over arbitrary polygons: Laplace interpolants, strain smoothing and scaled boundary polygon formulation., Finite Elements in Analysis and Design 85 (2014) 101–122. [30] C. Song, A matrix function solution for the scaled boundary finite-element equation in statics, Computer Methods in Applied Mechanics and Engineering 193 (2004) 2325–2356. [31] C. Talischi, G. H. Paulino, A. Pereira, I. F. M. Menezes, Polymesher: a general-purpose mesh generator for polygonal elements written in matlab, Structural and Multidisciplinary Optimisation 45 (2012) 309–328. [32] N. Sukumar, Quadratic maximum-entropy serendipity shape functions for arbitrary planar polygons, Computer Methods in Applied Mechanics and Engineering 263 (2013) 21–41. [33] N. Sukumar, J. H. Prevost, Modelling quasi-static crack growth with the extended finite element method. part i: Computer implementation, International Journal of Solids and Structures 40 (2003) 7513–7537. [34] M. M. Rashid, The arbitrary local mesh replacement method: an alternative to remeshing for crack propagation analysis, Computer Methods in Applied Mechanics and Engineering 154 (1998) 133–150. [35] P. O. Bouchard, F. Bay, Y. Chastel, I. Tovena, Crack propagation modelling using an advanced remeshing technique, Computer Methods in Applied Mechanics and Engineering 189 (2000) 723–742. [36] T. N. Bittencourt, P. A. Wawrzynek, A. R. Ingraffea, Quasi-automatic simulation of crack propagation for 2d lefm problems, Engineering Fracture Mechanics 55 (1996) 321–334. [37] D. Azocar, M. Elgueta, M. C. Rivara, Automatic lefm crack propagation method based on local lepp-delaunay mesh refinement, Advances in Engineering Software 41 (2010) 111–119.

17

[38] P. O. Bouchard, F. Bay, Y. Chastel, Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria, Computer Methods in Applied Mechanics and Engineering 192 (2003) 3887–3908. [39] A. R. Khoei, H. Azadi, H. Moslemi, Modelling crack propagation via an automatic adaptive mesh refinement based on modified superconvergent patch recovery technique, Engineering Fracture Mechanics 75 (2008) 2921– 2945.

Node

Hanging node

Figure 1: Hanging nodes in a generic quadtree mesh

Figure 1

18

η (x2,y2), η = 1

Line element

(x1,y1), η = −1

ξ=1

S

Scaling Centre

Radial Lines Similar curve to S

Figure 2: Scaled boundary finite element representation of a generic polygon

Figure 2

19

Node Hanging node

4-node cell

5-node cell

6-node cell

7-node cell

8-node cell

Figure 3: Cells in a quadtree mesh satisfying the 2:1 rule

Figure 3

20

Ω

Ω

(a)

(b)

(c)

(d)

R

Figure 4: Generation of a hybrid quadtree mesh: (a) a domain to be discretized; (b) background rectangular mesh and computational domain ΩR; (c) mesh refinement on boundary; (d) final hybrid quadtree mesh

Figure 4

21

hybrid quadtree

pure quadtree (a)

hybrid quadtree

pure quadtree (b)

Figure 5: Two generic meshes discretized by quadtree-polygon mesh and purely quadtree mesh: (a) Geometry 1; (b) Geometry 2 Figure 5

22

(a)

(b)

Figure 6: Close up view of generic domain discretized by quadtree-polygon mesh and pure quadtree mesh shown in Fig. 5

Figure 6

23

y

r

θ

x LA

A

Crack tip

Figure 7: Modelling of a crack using the SBFEM

Figure 7

24

(a)

(b)

(c) Figure 8: Generation of a crack cell in a hybrid quadtree mesh: (a) cells in the vicinity of crack tip; (b) subdivision of cells around crack tip; (c) construction of crack cell

Figure 8

25

Ωcir

refined regions

(b)

(a)

(c) Figure 9: Automatic remeshing procedure: (a) typical mesh at load step n; (b) quadtree refinement; (c) generation of new crack cell

Figure 9

26

σ (0,20)

(15,20)

(5,11.9)

R=3.45 (2.75,5.3)

(0,0)

(15,0)

σ

(a)

(b)

Figure 10: Crack path deflection in a plate due to a hole: (a) geometry; (b) hybrid quadtree mesh

Figure 10

27

Present Rashid (1998) Bouchard et al. (2000)

(a)

(b)

Figure 11: Predicted crack trajectory for the crack path deflection in a plate due to a hole problem: (a) comparison with typical results reported in the literature; (b) final mesh

Figure 11

28

(a)

(b)

Figure 12: Medium and fine quadtree-polygon meshes for crack path deflection due to a hole problem: (a) medium mesh; (b) fine mesh

Figure 12

29

Coarse Medium Fine

Figure 13: Comparison of predicted crack path by various mesh densities

Figure 13

30

P 1.25

r = 0.25

2.0 8.0

2.0 b 4.0

a

1.0

9.0

9.0

CASE

a

I

1.5

5.0

II

1.0

4.0

1.0

b

Figure 14: Cracked beam with three holes

Figure 14

31

(a)

(b) Figure 15: Initial meshes of the cracked beam with three holes: (a) Case I; (b) Case II

Figure 15

32

(a)

(b) Figure 16: Final crack paths for cracked beam with three holes: (a) Case I; (b) Case II

Figure 16

33

Present Bittencourt et al (1996) Azocar et al (2010)

Present Bittencourt et al (1996) Azocar et al (2010)

(a)

(b)

Figure 17: Comparison of predicted crack paths for the cracked beam with three holes with published results in the literature: (a) Case I; (b) Case II

Figure 17

34

F (0,10)

(20,10) (19,7.15)

(3,7) R=2

R=2

(1,2.85)

(17,3) (20,0)

(0,0)

F

(a)

(b) Figure 18: Multiple crack propagation in a plate with two holes: (a) geometry; (b) hybrid-quadtree mesh

Figure 18

35

Present Bouchard et al. (2003) Khoei et al. (2008)

(a)

(b) Figure 19: Predicted crack path for the multiple crack propagation in a plate with two holes problem: (a) comparison with finite element solution of Bouchard et al [35]; (b) final hybrid quadtree mesh

Figure 19

36

σ=2

2.5 r = 5/64

θ 2

1

2

2.5

(a)

(b)

Figure 20: Propagation of two cracks emanating from holes: (a) geometry; (b) hybrid quadtree mesh

Figure 20

37

(a)

Present Moes et al (1999)

(b) Figure 21: Predicted crack path for the multiple crack propagation in a plate with two holes problem: (a) close up view of crack path; (b) comparison with results published in the literature

Figure 21

38

Highlights     

The SBFEM is combined with quadtree meshes for fracture modelling The quadtree structure leads to rapid mesh generation and stiffness computations The SBFEM accurately models the asymptotic stress field at the crack tip Crack evolution is conveniently modelled by splitting quadtree cells The developed approach is accurate and efficient for crack propagation problems

Nomenclature

1

Nomenclature B D K N S Z c q u K Ψ Ψσ θ λ ξ η dΩ sΩ Ω ∂Ω |f |

scaled boundary finite element strain-displacement matrix material constitutive matrix stiffness matrix shape function matrix Schur/eigenvalue matrix Hamiltonian matrix integration constants internal force vector displacement vector stress intensity factors transformation matrix stress mode angle eigenvalue radial coordinate local coordinate signed distance function sign function domain boundary of domain crack propagation length