Engineering with Computers https://doi.org/10.1007/s00366-017-0555-5

ORIGINAL ARTICLE

Simulating fracture propagation in brittle materials using a meshless approach J. Belinha1,2   · J. M. C. Azevedo1 · L. M. J. S. Dinis1 · R. M. Natal Jorge1 Received: 2 September 2016 / Accepted: 14 November 2017 © Springer-Verlag London Ltd., part of Springer Nature 2017

Abstract In the present work, a meshless method is combined with a crack growth algorithm, considering a linear elastic model to simulate the propagation of cracks in brittle materials. The proposed methodology is automatic, which means that it does not require the user intervention. In this study, the meshless method used is the Natural Neighbour Radial Point Interpolation Method (NNRPIM), an efficient discrete numeric method. Being a truly meshless method, the NNRPIM only requires an unstructured nodal distribution to fully discretize the problem domain. With the coordinates of the nodes, the NNRPIM formulation is autonomously capable: to define the nodal connectivity; to build the background integration mesh; and to construct the shape functions. The crack propagation algorithm here proposed simulates the crack growth by displacing iteratively the crack tip. Updating the crack tip location requires a local remeshing, which does not represent a numeric difficulty for the NNRPIM. The updated position of the crack tip depends on the estimated stress field in the vicinity of the crack tip. Thus, in each iteration, the stress field in the vicinity of the crack tip is estimated, and then, the direction of the crack propagation is calculated using the maximum circumferential stress criterion. Keywords  Meshless methods · Fracture · Crack growth · NNRPIM · Radial basis functions

1 Introduction The computational fracture mechanics have been extensively investigated over the recent years. This computational mechanic area deals with the development and enhancement of numerical techniques capable of simulating and analysing the numeric singularity created by a crack tip. In the literature, it is possible to find several numerical approaches capable to analyse fracture mechanics problems [1]. Nevertheless, most of the numerical methods developed to study the fracture problem are based on the finite element method (FEM). The first FEM works dealing with the crack propagation appear in the early 1970s [2], and since then, the FEM has

* J. Belinha [email protected] 1



Faculty of Engineering of the University of Porto, FEUP, Rua Dr. Roberto Frias, 4200‑465 Porto, Portugal



INEGI, Institute of Mechanical Engineering and Industrial Management, Rua Dr. Roberto Frias, 4200‑465 Porto, Portugal

2

been intensively used in the numerical analysis of fracture mechanic problems [3]. The simulation of the crack propagation is a non-linear procedure, requiring a recurrent modification of the problem physical domain. Consequently, the computational discretization (the element mesh) has to be frequently updated due to the crack growth. This remeshing task can be a burdensome and an expensive computational FEM task. Furthermore, the automatic remeshing could lead to highly distorted elements, worsening the performance of the FEM. Initial attempts to improve the FEM formulation for fracture mechanics included the hybrid-element approach [4, 5] and the fractal finite element method (FFEM) [6]. More recently, other FEM versions were developed, such as the extended finite element method (XFEM) [7], the node-based smoothed extended finite element method (NS-XFEM) [8] and edge-based smoothed finite element method (ESFEM) [9, 10]. The literature also shows that the strain smoothing in FEM and XFEM allows to achieve superior results in crack problems [11]. Nowadays, the XFEM is the most popular FEM version for computational fracture mechanics. The efficiency and

13

Vol.:(0123456789)



the accuracy of the XFEM is widely attested in the literature [12–14]. Nevertheless, all FEM versions available in the literature are mesh-based discrete approximation techniques. In the literature, it is possible to find advanced numerical techniques that do not require a structured mesh to discretise the problem domain—the meshless methods [15, 16]. Meshless methods permit to discretize the problem domain using only a set of computational points or nodes, without pre-establishing any sort of connectivity between them [17]. Alternatively to the “element” concept of the FEM, in meshless methods, the nodal connectivity is enforced with the “influence-domain” concept [15]. Then, the shape functions of meshless methods are constructed using the nodes inside each influence domain. Meshless methods using the weak form of Galerkin can be classified as approximation meshless methods [18–23] or as interpolation meshless methods [24–30]. When compared with the interpolation meshless methods, the approximation meshless methods show lower convergence rates [31] and are unable to directly impose the essential and natural boundary conditions [17]. Since interpolation meshless methods are capable to construct shape functions possessing the Kronecker delta function property, for these methods the essential boundary conditions can be directly imposed in the system of equations, as in the FEM. In this work, it is considered an interpolation meshless method—the natural neighbour radial point interpolation method (NNRPIM) [17, 28, 32]. The NNRPIM combines the radial point interpolators (RPI) [27, 33] with the natural neighbour concept [34]. When compared with other meshless method formulations, the NNRPIM formulation presents some notorious advantages [17]. For instance, the only computational mesh required by the NNRPIM formulation is the unstructured nodal distribution discretizing the problem domain. With the spatial information of the nodal discretization the NNRPIM is capable to autonomously (without the user intervention): obtain the influence-domain; build the background integration mesh; and construct the shape functions. Thus, the NNRPIM is a truly meshless method. In the literature, it is possible to find research works attesting the accuracy and efficiency of the NNRPIM in several computational mechanics fields, such as the analysis of laminated shell-like structures [35, 36], the analysis of functionally graded plates [37], the dynamic and transient analysis [38–40], the non-linear structural analysis [41, 42] and the non-linear biomechanical analysis [17, 43–46]. Meshless techniques seem more suitable to solve fracture mechanics problems. The stress field obtained with meshless methods is generally more accurate [17, 21] and more nodes can be easily inserted in the vicinity of the crack tip. The potential of meshless methods on fracture mechanics computation is evident ever since the early published works [21]. Thus, the development of crack opening algorithms

13

Engineering with Computers

combined with meshless methods is a topic of significant research activity [47, 48]. In the literature, it is possible to find several works showing the efficiency of meshless methods in the prediction of the crack growth, using a classic formulation [49, 50] or an enriched formulation [51–53]. For instances, in Rabczuk et al. [54] an enriched version of the element free Galerkin method (EFGM) was successfully used to model arbitrary crack initiation and crack growth in reinforced concrete structures using 3D geometrically non-linear formulation. Furthermore, in Rabczuk et al. [55] the EFGM is combined with a non-enriched formulation in which the crack is modelled by splitting particles located on opposite sides of the associated crack segments. Thus, using the visibility method, Rabczuk et al. were capable to predict 2D and 3D crack opening paths. This study aims to extend the NNRPIM to the simulation of the crack opening path using a linear elastic model. In some practical circumstances, there are cases in which is preferable to predict the crack growth using a linear elastic model, such are some cases of fatigue crack propagation and quasi-stable crack growth under compressive loading [56]. Therefore, in this work, an iterative crack opening procedure is adapted to the NNRPIM formulation. In each iteration of the analysis, the crack propagation direction is calculated using the maximum circumferential stress criterion [57]. Then, the crack is extended in a straight line segment and the nodal mesh is locally updated [58]. In the end, benchmark examples are studied considering the proposed meshless approach. The obtained solutions are compared with the crack opening paths predicted by other numerical approaches or with experimental results.

2 Meshless method In this work, to iteratively update the crack tip spatial position, the considered crack growth algorithm requires the stress field. Here, the stress field is obtained using a flexible and organic interpolation meshless method—the natural neighbour radial point interpolation method (NNRPIM) [28]. The NNRPIM uses mathematic concepts, such as the Voronoï diagram [59] and the Delaunay tessellation [60], to construct the influence cells (the basic structure of the nodal connectivity in the NNRPIM) and the background integration mesh. Both these numerical structures are completely dependent on the nodal mesh. The radial point interpolators (RPI) permit the construction of the NNRPIM interpolation functions, which afterwards are used as test functions in the Galerkin weak form. Additionally, the RPI function construction is simple and its derivatives are easily obtained. Since the RPI functions possess the delta

Engineering with Computers

Kronecker property, the boundary conditions can be directly imposed as in the FEM. When compared with the FEM, the NNRPIM shows a higher accuracy for the obtained variable field and it is more flexible in the remeshing process [35, 36]. In addition to the previously pointed advantages, the NNRPIM presents acceptable computational costs [17], which is relevant in non-linear studies, and allows to obtain high quality stress fields [17]. In this section, the main features of the NNRPIM are presented.

2.1 Nodal connectivity and integration mesh The natural neighbour concept [34] is used to enforce the nodal connectivity of the NNRPIM. The natural neighbour concept is materialized by the construction of the Voronoï diagram of the discretized domain. This geometric diagram is constituted by a set of Voronoï cells. To each one of these cells is associated one and only one node. Since it is simpler to represent graphically, only twodimensional examples are shown in this section. However, this geometric procedure can be applied to a d -dimensional space. Consider a domain Ω ⊂ ℝ2 bounded by a physical boundary, Γ ∈ Ω . Then, the domain Ω{ is discretized}with an unstructured nodal distribution: N = n0 , n1 , … , nN ∈ ℕ . The {Cartesian coordinates of N are represented as: } X = x0 , x1 , … , xN ⊂ Ω with xi ∈ ℝ2. The Voronoï diagram of N is defined by the partition of the domain Ω into Voronoï cells VI  , closed and convex. Each sub-region VI is associated with the node nI  , in a way that any point in the interior of the VI is closer to nI than any other node nJ ∈ N(J ≠ I) . The set of all Voronoï cells Vk , k = 1 , … , N  , defines the Voronoï diagram, as shown in Fig. 1a. In the literature, it is possible to find research works describing with detailed efficient methodologies to obtain the Voronoï diagram [17, 28, 29].

Using the Voronoï diagram, it is possible to enforce the nodal connectivity. Generally, in meshless methods, the nodal connectivity is imposed using the “influence-domain” concept [16], which corresponds approximately to the “element” concept of the FEM. Several techniques used to construct influence domains are well documented in the literature [16, 17]. The NNRPIM imposes the nodal connectivity using the “influence cell” concept [28]. The influence cells are organic influence domains that depend on the information from the Voronoï diagram. Considering the interest point xI ⊂ Ω , shown in Fig. 1a, the influence cell of xI is constituted by the interest point xI and the direct natural neighbour nodes of xI  , indicated in grey in Fig. 1a. This procedure only depends on the nodal discretization. Besides the nodal spatial position, no other information is required. Discrete numeric techniques using the Galerkin weak form require a background integration mesh to numerically integrate the integro-differential equations ruling the studied physical phenomenon. The NNRPIM is no exception. For the majority of meshless methodologies, the construction of the background integration mesh does not depend on the nodal distribution. However, in the NNRPIM, the integration mesh is obtained using directly and exclusively the nodal distribution [28]. It is possible to construct the background integration mesh using the information from the Voronoï diagram. First, using the Delaunay triangulation concept, the area of each Voronoï cell is subdivided in several subareas, as shown in Fig. 1(b). Afterwards, the area SI of the Voronoï cell VI of ∑k node xI is divided in k subareas SIj , being SI = j=1 SIj , as shown in Fig. 1b. Following the Gauss–Legendre quadrature rule, it is possible to distribute integration points inside each sub-area SIj . This procedure permits to obtain the integration mesh of the Voronoï cell VI  . Repeating this process for the N Voronoï cells from the Voronoï diagram, it is possible to determine the background integration mesh discretizing the problem

Fig. 1  a Voronoï diagram and influence cell of interest point xI; b subdivision of a Voronoï cell; c quadrilateral used to insert the integration points

13



Engineering with Computers

∑N ∑k domain Ω ⊂ ℝ2 . In the end: SΩ = I=1 j=1 SIj , being SΩ the area of the discretized domain. In this work, the integration mesh is constructed considering just one integration point per sub-area SIj , as shown in Fig. 1c, since previous research works on the NNRPIM shown that this simple integration scheme is sufficient to integrate accurately the integro-differential equations [17, 28, 32, 61]. This procedure is described with detail in the literature [17, 28, 32, 62, 63]. Since it is not possible to know in advance the most efficient integration scheme to accurately integrate a weak form whose approximation depends on the multiquadric (MQ) RBF function, in [28] four integrations schemes were tested. A first one, using one integration point per sub-area (order 0), and three others, in which the sub-area is divided in quadrilaterals and inside each quadrilateral are inserted 1 × 1 integration points (order 1), 2 × 2 integration points (order 2), and 3 × 3 integration points (order 3). The results of [28] show that the integration scheme of order 0 is more efficient, since it allows to obtain a solution very close to the one obtained using higher integration orders (such as order 3), and at the same time, it allows to save time. For the same analysis, using an integration scheme of order 0 allows to obtain the result: three times faster than integration scheme of order 1; 12 times faster than integration scheme of order 2; and 27 times faster than integration scheme of order 3.

2.2 NNRPIM shape functions In this work, the shape functions are constructed using the radial point interpolators (RPIs) [27]. The classical RPI requires a complete polynomial basis and a radial basis function (RBF). Nevertheless, previous works on the NNRPIM showed that the polynomial basis can be removed from the formulation and substituted by a constant unity basis [28, 32]. This simplification can be made if the shape parameters of the RBF are chosen carefully. With this innovation, the global NNRPIM performance was enhanced and the computational efficiency of the NNRPIM was increased. To understand the RPI function construction, consider a function u(x) defined in the domain Ω ⊂ ℝd  . Additionally, it is assumed that an interest point xI ∈ Ω possesses an influence cell containing n nodes: XI = {x1 , x2 , … , xn } ∈ Ω ∧ xi ∈ ℝd . Since the RPI function have contact support, only the nodes within the influence cell of the interest point xI will have effect on u(xI ) . Thus, for the generic point of interest xI  , the value of the function u(xI ) is defined by n ( ) ( ) ∑ [ ( ) ( )] u xI = Ri xI ⋅ ai xI + Cu ⋅ b xI i=1

( ) ( ) ( ) = R xI ⋅ a xI + Cu ⋅ b xI ,

13

(1)

where Ri (xI ) is the RBF and ai (xI ) are non-constant coefficients of Ri((xI )) . The parameter Cu is the unity basis, being Cu = 1 and b xI is a non-constant coefficient of Cu. The variable of the RBF is the distance rI i between the relevant node xI and the neighbour node xi , rI i ∶= |xi − xI | . The NNRPIM formulation uses (MQ) RBF ( the multiquadric )p function, Ri (xI ) = R(rIi ) = rIi 2 + c2  , in which c and p are two shape parameters. The MQ-RBF was initially proposed for surface fitting [64]. Previous works on the NNRPIM found that c should be close to zero, c ≅ 0 , and p should be close to one, p ≅ 1 , [28, 32]. However, these values cannot be c = 0 and p = 1 . The use of the exact integer value leads to a dependent system of equation, which does not permit to obtain a solution. Therefore, the literature [17] recommends the following values: c = 0.0001 and p = 1.0001 . These c and p values can be used in the construction of the interpolation functions, regardless the studied phenomenon [28, 32, 61]. Applying Eq. (1) to each node inside the influence cell, sequentially considering each node as the interest point, and ∑n including an extra equation i=1 ai (xI ) = 0 to guarantee the unique solution [65], it is possible to obtain the following equation system:

⎡ u1 ⎤ ⎡ R(r11 ) R(r12 ) ⋯ R(r1n ) 1 ⎤⎡ a1 (xI ) ⎤ ⎢ u ⎥ ⎢ R(r ) R(r ) ⋯ R(r ) 1 ⎥⎢ a (x ) ⎥ � � 21 22 2n ⎢ 2⎥ ⎢ ⎥⎢ 2 I ⎥ us ⋮ ⋱ ⋮ ⋮ ⎥⎢ ⋮ ⎥ ⇔ ⎢ ⋮ ⎥=⎢ ⋮ 0 ⎢ un ⎥ ⎢ R(rn1 ) R(rn2 ) ⋱ R(rnn ) 1 ⎥⎢ an (xI ) ⎥ ⎢0⎥ ⎢ 1 1 ⋯ 1 0 ⎥⎦⎢⎣ b(xI ) ⎥⎦ ⎣ ⎦ ⎣ � �� � R Cu a (2) = . b Cu T 0 The vector uS is constituted with the nodal function values uS = {u1 u2 … un }T of the nnodes inside the influence cell of interest point xI  . From Eq. (2), it is possible to obtain the non-constant coefficients a(xI ) and b(xI ) ,

[ ] [ ]−1 [ ] [ ] [ ] R Cu a us a −1 us = ⇒ =M . b b 0 0 Cu T 0

(3)

Thus, for an interest point xI  , the interpolation function 𝚽(xI ) is found by substituting in Eq. (1) the results obtained in Eq. (3), ( ) { } u xI = R1 (xI ) R2 (xI ) ⋯ Rn (xI ) Cu M−1

{

us 0

}

{

( ) = 𝚽 xI

us 0

}

.

(4) The last term of the interpolation function, 𝜓C , is associated with the constant term Cu and it does not possess any useful physical meaning. Therefore, only the n first components of 𝚽(xI ) = {𝜑1 (xI ), 𝜑2 (xI ), … , 𝜑n (xI ), 𝜓C } are relevant and constitute the interpolation function, [17]. The partial

Engineering with Computers

derivative of 𝚽(xI ) with respect to the spatial variables are defined as } { ( ) { 𝚽,x xI = R1 (xI ) R2 (xI ) ⋯ Rn (xI ) Cu ,x M−1 . ( ) { } (5) 𝚽,y xI = R1 (xI ) R2 (xI ) ⋯ Rn (xI ) Cu ,y M−1 Notice that since Cu is a constant scalar, then Cu,x = Cu,y = 0 . The partial derivatives of the MQ-RBF with respect to the spatial variables are obtained with

� �p−1 � � ⎧ � � 2 2 R r = 2p r + c xj − xi ij ⎪ ,x ij . ⎨ � � � �p−1 � � ⎪R,y rij = 2p r2 + c2 y − y j i ij ⎩

(6)

2.3 Linear elastic analysis using the NNRPIM Consider the body described by the domain Ω ⊂ ℝ and bounded by Γ , where Γ ∈ Ω ∶ Γu ∪ Γt = Γ ∧ Γu ∩ Γt = � , being Γu the essential boundary and Γt the natural boundary. The equilibrium equations governing the linear elastostatic problem are defined in Ω by ∇𝚲 + b = 0 . The symbol ∇ represents the gradient operator; the 𝚲 stands for the Cauchy stress tensor for a kinematically admissible displacement field u ; and the vector b represents the body force per unit volume. The boundary conditions are given by, 𝚲 n = ̄t on Γt and u = ū on Γu , where ū is the prescribed displacement on the essential boundary Γu , ̄t is the traction on the natural boundary Γt and n is the unit outward normal to the boundary of domain Ω . Using the Voigt notation and assuming the Galerkin procedure for linear elasticity, the weak form for the discrete problem can be written as, 2

∫Ω

𝛿𝜺T 𝝈 dΩ −

∫Ω

𝛿uT b dΩ −

∫Γt

𝛿uT ̄t dΓ = 0,

(7)

where 𝜺 is the strain vector, defined as 𝜺 = L u , and L the differential operator,

⎡ 𝜕 0 ⎤ ⎢ 𝜕x 𝜕 ⎥ L = ⎢ 0 𝜕y ⎥. ⎢ 𝜕 𝜕 ⎥ ⎣ 𝜕y 𝜕x ⎦

c=

(8)

It is possible to directly correlate the stress field with the strain field using the Hooke Law: 𝝈 = c𝜺 = c Lu . The square matrix c is the material constitutive matrix, and for a homogeneous and isotropic material, it depends on the following material properties: the elasticity modulus E and

⎡1 𝜐 0 ⎤ E ⎢ 𝜐 1 0 ⎥, ⎥ 1 − 𝜐2 ⎢ 0 ⎣ 0 1−𝜐⎦

(9)

and for the plane strain as,

c=

In the literature [17, 27–29], it is possible to find complete studies on the RPI interpolation functions: construction procedure and properties.

𝛿L =

the Poisson’s ratio 𝜐 . Thus, the matrix c can be defined for the plane stress as,

0 ⎤ ⎡1−𝜐 𝜈 E ⎢ 𝜈 1 − 𝜐 0 ⎥. ⎥ (1 + 𝜐)(1 − 2𝜐) ⎢ 0 1 − 2𝜐 ⎦ ⎣ 0

(10)

In the NNRPIM, the weak form has local support. Therefore, the discrete system of equations is developed first for each influence cell and later the local systems of equations are assembled to form the global system of equations. The NNRPIM trial function is given by Eq. (4). Thus, for each degree of freedom it is possible to write,

uh (xI ) =

n ∑

𝜑i (xI ) u(xi ) and vh (xI ) =

i=1

n ∑

𝜑i (xI ) v(xi ), (11)

i=1

where 𝜑i (xI ) is the NNRPIM interpolation function and u(xi ) and v(xi ) are the nodal parameters of the ith node belonging to the nodal set defining the influence cell of interest node xI  . Both expressions in Eq. (11) can be combined in one single equation, { h } } ]{ n [ ∑ u (xI ) 𝜑i (xI ) 0 u(xi ) h u (xI ) = = vh (xI ) 0 𝜑i (xI ) v(xi ) i=1

=

n ∑

(12)

Hi (xI ) u(xi ).

i=1

Thus, the linear relation between the strain field and the displacement field: 𝜺 = L u , can be expressed as,

𝜺(xI ) = L uh (xI ) = L

n �

Hi (xI ) u(xi ) =

i=1

⎡ n � ⎢ = ⎢ i=1 ⎢ ⎣

n �

Bi (xI ) u(xi )

i=1 𝜕𝜑i (xI ) 𝜕x

0

𝜕𝜑i (xI ) 0 𝜕y 𝜕𝜑i (xI ) 𝜕𝜑i (xI ) 𝜕y 𝜕x

⎤� � ⎥ u(xi ) ⎥ v(x ) . i ⎥ ⎦

(13)

Here, the deformation matrix is defined as B = LH. Additionally, using the Hooke law, it is possible to correlate the stress state in the interest point xI with the strain state, 𝝈(xI ) = c𝜺(xI ) . Thus, substituting the strain vector 𝜺(xI ) and the stress vector 𝝈(xI ) in the first term of Eq. (7) and the approximation function on Eq. (12) in the second and third term, it is possible to rewrite Eq. (7) for the interest point xI as,

13



Engineering with Computers

n n ∑ ∑

𝛿uTi

i=1 j=1

∫Ω ⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟

BTi cBj dΩ uj − (K)ij

∑

{

n

−

𝛿uTi

i=1

n ∑

} bx dΩ by ∫Ω ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ {

𝛿uTi

i=1

}

HTi

(f b )i

t HT x dΓt = 0. ∫Γt i ty ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟

(14)

(f t )i

The equation system can be represented in the matrix form: K u = f  . The square matrix K represents the stiffness matrix, the vector u is the displacement field vector and f = f b + f t is the vector of applied forces. Since the RPI shape function possesses the delta Kronecker property, the essential boundary conditions can be directly applied in the stiffness matrix [28]. After the determination of the displacement field, solving the linear equation system K −1 f = u , the strain in an interest point xI ∈ Ω can be obtained using 𝜺(xI ) = L u(xI ) , and then, considering the Hooke Law, 𝝈(xI ) = c(xI ) 𝜺(xI ) , it is possible to determine the stress field. Additionally, it is possible to obtain the ( two principal )stresses 𝜎(xI )i for each interest point xI , det 𝚲(x ( I ) − 𝜎(xI )i I = ) 0 , and the two principal directions n(xI )i : 𝚲(xI ) − 𝜎(xI )i I n(xI )i = 0 . Here, the Cauchy stress tensor is represented by 𝚲(xI ) and I is a 2 × 2 identity matrix.

3 Crack propagation with the NNRPIM The research work, here presented, allowed the development of the main numeric tools that enable the NNRPIM formulation to analyse the crack growth phenomenon considering a linear elastic model. Thus, in this section, the most important developed numeric tools are explained. First, the NNRPIM crack propagation iterative algorithm is presented, and then, the numeric techniques used to obtain the crack growth direction and the new location of the crack tip are described with detail.

3.1 Crack propagation algorithm In this work, to efficiently predict the crack propagation, an iterative algorithm was adapted and developed to take advantage of the unique features of the NNRPIM formulation. The schematic representation of the iterative algorithm is shown in Fig. 2. Thus, in step 1, the following tasks are performed: the solid domain is identified; the geometry of the structure is obtained; the material properties are established; the displacement constrains are fixed; and external forces are identified.

13

In the next step, step 2, the domain of the studied solid is discretized with a computational nodal distribution, which can be regular or irregular. Then, in step 3, using the procedure described in Sect. 2.1, the Voronoï diagram is built and the background integration mesh is constructed. Both these geometric structures are obtained uniquely from the initial computational nodal distribution and no other information is required from the user. In step 4, the initial crack is created in the numeric model. The procedure described in Sect. 3.4 is applied in this step. Hence, the initial Voronoï diagram (the nodal connectivity) and the integration mesh are adapted to the new domain configuration, with the crack opening. As Fig. 2 shows, there are no integration points between the crack opening edges. Notice that the crack opening edges became a domain boundary. After the pre-processing tasks are concluded (step 1 to step 4), in step 5 the NNRPIM shape functions are built. Afterwards, the local stiffness matrix is constructed for each interest point and then, the local stiffness matrices are assembled into the global stiffness matrix. Next, the global external force vector is calculated and the displacement constrains are imposed. As step 6 indicates, after defining the system of equations, the displacement field is obtained: u = K −1 f  , which permits to obtain the stress and strain fields, step 7. Then, as explained in Sect. 3.2, the weighted stress tensor in the vicinity of the crack tip is found. The principal stress tensor and the principal directions of the weighted stress tensor are obtained, as Sect. 2.3 specifies, and the direction of the crack propagation is determined, as Sect. 3.2 indicates. Afterwards, in step 8, the crack is extended in straight line segments, following the direction previously obtained, and domain discretization is updated. The crack tip update procedure is described with detail in Sects. 3.2 and 3.3. If the crack has reached its final length, the analysis ends. If not, the process moves to step 4 and  a new iteration begins.

3.2 Direction of the crack propagation The algorithm developed in this work assumes that the crack propagation direction is orthogonal to the first principal stress direction [57]. This assumption is in accordance with the maximum circumferential stress criterion [57]. This work uses this simple failure criterion because its main objective is to develop the most relevant numerical tools required to analyse the crack propagation phenomenon with the NNRPIM formulation. After the validation of the developed methodology, other failure criteria can be easily included in the algorithm. Thus, to obtain the crack propagation direction, first, the Euclidean distance between each integration point xI and the crack tip is obtained: d0I .

Engineering with Computers Fig. 2  Developed algorithm to predict the crack growth with the NNRPIM

Then, as Fig. 3a shows, the integration points satisfying the condition d0I < dref are selected. The scalar parameter dref should be defined in the pre-processing phase. The benchmark examples solved in Sect. 4 study the influence of this parameter in the final solution. Afterwards, the stress state of all integration points within the reference distance dref is considered to calculate the estimated stress state on the crack tip, 𝝈(x0 ),

nq ∑

𝝈(x0 ) =

w(x0 , xI ) ⋅ 𝝈(xI )

I=1 nq ∑

.

(15)

w(x0 , xI )

I=1

The number of integration points within the reference distance dref is defined with nq . The function w(x0 , xI ) represents a monotonic-decreasing weight function centred in n0 . The

13



Engineering with Computers

Fig. 3  a Radial search of the integration points. b Crack growth direction and elimination distance. c Update of crack tip and crack edge nodes

weight function w(x0 , xI ) has compact support, i.e. w(x0 , xI ) is positive and non-null for d0I < dref : { 1 − d0I ∕dref , if d0I ⩽ dref w(x0 , xI ) = . (16) 0, if d0I > dref After the determination of the estimated stress state on the crack tip, 𝝈(x0 ) , it is possible to determine at the crack tip the respective two-dimensional principal stresses ( 𝜎1 and 𝜎2 ) and principal directions ( v1 and v2). Since the maximum circumferential stress criterion postulates that the crack propagates in a direction orthogonal to the first principal stress direction, the crack propagation vector is defined as v2.

3.3 Explicit crack propagation To simulate the crack opening phenomenon, it is necessary { } to update the nodal distribution N = n1 , n2 , … , nN discretizing the problem domain Ω ∈ ℝ2. Once again, a graphic example is used to explain the procedure. Consider the discrete set-up represented in Fig. 3b. First, as explained in Sect. 3.2, the direction of the crack propagation is determined: v2 . Afterwards, all nodes ni satisfying the condition: d0i < delim are eliminated from the nodal set N  , Fig. 3b. In the former condition, the Euclidean distance between node ni and node n0 is represented as: d0i . The scalar parameter delim ∈ ℝ is fixed in the beginning of the analysis. It depends on how much the crack tip position moved in each iteration Δa , as shown in Fig. 3c. In this work, all studied examples use: delim ≅ Δa , which permits to obtain accurate and stable results. Next, as Fig.  3b shows, the old crack tip node n0 is divided in two other nodes: n3 and n4 . The Euclidean distance between nodes n3 and n4 is n3 n4 = da . The crack

13

opening distance, da , is defined in the beginning of the analysis. Being n1 and n2 the first nodes of the crack opening, then: da = n1 n2  . The coordinates of nodes n3 and n4 are obtained with the expressions: x3 = xold − (da ∕2) ⋅ v1 and 0   . The coordinates of the new crack tip x4 = xold + (d ∕2) ⋅ v a 1 0 old is determined with xnew . = x + Δa ⋅ v 2 0 0

3.4 Discretization of the crack opening To numerically study the crack propagation, first, it is necessary to discretize the opening edges of the initial crack. Consider a solid domain Ω ∈ ℝ2 bounded { by Γ ∈ Ω and } discretized with a nodal mesh: N = n1 , n2 , … , nN  . The {correspondent } Cartesian coordinates of N are: X = x1 , x2 , … , xN ∈ Ω  . Assume a crack opening defined by the contour ΓC ∈ Γ ∧ ΓC ∈ Ω . The crack opening is discretized with the nodal set NK ⊂ N  , being XK ∈ ΓC ∧ XK ∈ X the Cartesian coordinates of NK  . It is possible to divide the nodal set NK in two sub-sets: NK1 and NK2 . Each one of these sub-sets are associated with a physical side of the crack edge: NK1 ∪ NK2 = NK ∧ NK1 ∩ NK2 = �. This methodology is widely used in meshless methods [51, 66, 67] because it permits to easily identify to which side of the crack the edge nodes belong to. Since it is simple, assume the following practical example shown in Fig. 4a. Consider that the crack opening is defined by the nodal set: NK = {n0 , n1 , n2 , … , n6 } . At this point, the crack tip has moved from a previous location to a new location, n0 . The old crack tip node was divided in two new nodes: n1 and n2 . It is necessary to decide to which side n1 and n2 belong to. Thus, it was developed a simple algorithm to perform automatically this task. No user intervention is required. From the last iteration, the nodes {n3 , n5 } and {n4 , n6 } were already associated with NK1 and NK2 , respectively, {n3 , n5 } ∈ NK1 and {n4 , n6 } ∈ NK2.

Engineering with Computers Fig. 4  a Generic view of the relative position of line “s”; b–d remaining cases for the relative position of the line “s”

Hence, the following sequent nodes are selected: {n1 , n2 , n3 , n4 } , and a line s is defined passing through node n3 ∈ NK1 and node n2  , s ∶ y = a ⋅ x + b . Since the coordinates of node ni ∈ N are xi = {xi , yi } ∈ ℝ2 , the constants a and b can be defined as a = (y2 − y3 )∕(x2 − x3 ) and b = (x2 y3 − x3 y2 )∕(x2 − x3 ) . Afterwards, it is possible to test for nodes n1 and n4 the following conditions: { a ⋅ x1 + b − y1 < 0, condition 1 . (17) a ⋅ x4 + b − y4 < 0, condition 2 If n1 and n4 are in opposite sides, then, it is because condition 1 is false and condition 2 is true (or viceversa). Therefore, nodes n3 and n1 belong to the same side {n1 , n3 } ∈ NK2 , and consequently, {n2 , n4 } ∈ NK1. If both conditions are true or both conditions are false, it is because nodes n1 and n4 belong to the same side, i.e. {n1 , n4 } ∈ NK1 , and consequently, {n2 , n3 } ∈ NK2 . The four possible solutions are presented in Fig. 4a–d. Thus, with this process all nodes NK ⊂ N will be univocally associated with NK1 or NK2. Afterwards, it is necessary to update the natural neighbours set of each node ni ∈ N  . Since all nodes ni ∈ NK1 cannot belong to the influence cells of the nodes nj ∈ NK2 , all the nodes ni ∈ NK1 will be removed from the nodal set forming the influence cells of nodes nj ∈ NK2 . The same procedure is performed, now with respect to the influence cells of nodes ni ∈ NK1.

4 Examples of applications In this work, three benchmark examples are solved to show the performance of the NNRPIM formulation in computational fracture mechanics. The automatic crack propagation algorithm here developed requires the introduction of two sensitive parameters: the crack increment length Δa and the selection radius dref (both described with detail in Sect. 3).

Thus, in this section, the influence of the parameters Δa and dref in the accuracy of the simulated crack growth is studied. Both these parameters depend on a third parameter: the inter-nodal distance near the crack h . Considering the 2 problem { domain Ω }∈ ℝ discretized in a nodal distribution X = x1 , x2 , … , xN ∈ Ω , the inter-nodal distance near the crack can be expressed as,

‖ ‖ h = min‖xj − xi ‖, ∀{i, j} ∈ ℕ ∶ {i, j} ⩽ N ∧ i ≠ j, ‖ ‖

(18)

where || ⋅ || is the Euclidean norm. Since this work assumes a linear elastic model, the structural numerical analysis only requires the following material properties: Young’s modulus, E  , and the Poisson’s coefficient, 𝜐 . Thus, in all the examples are considered the following values: E = 1000 Pa and 𝜐 = 0.29.

4.1 Four point bending beam In this first example, as Fig. 5a schematically shows, a four point bending beam with a circular opening is studied. In Fig. 5a, it is possible to identify the geometric properties of the beam, and additionally, the displacement constrains and the applied forces. In the mid-span of the beam exists an initial vertical crack with the following length: a = 2.5 mm. To perform the simulation, the beam physical domain was discretized in a computational nodal distribution. In this example, two kinds of nodal discretization were used: a uniformly distributed nodal scattering, shown in Fig. 5b, d, and a regular nodal discretisation presenting a higher nodal density in the crack vicinity, as represented in Fig. 5c, e. This example was solved assuming the plain strain deformation theory. All the results obtain with the NNRPIM formulation are compared with the results obtained with the ESFEM formulation [9]. First, it is studied the influence of the nodal discretisation on the simulation of the crack propagation. Thus, this numerical example was studied considering the four

13



Engineering with Computers

Fig. 5  a The four point bending beam with circular opening and an initial crack. Nodal meshes used in the analysis: b 41 × 11 uniform mesh (461 nodes); c 41 × 11 locally refined mesh (769 nodes); d 81 × 21 uniform mesh (1845 nodes); e 81 × 21 locally refined mesh (3071 nodes)

0.025

0.0175

0.020

0.015 y [m]

0.0150

41x11 0.0125

41x11 (refined) 0.010

0.005

y [m]

Fig. 6  a Prediction of the crack growth of the four point bending beam, for different nodal meshes. dref = 8h and Δa = h . b Close-up detail near the circular opening

81x21

0.0100

81x21 (refined) Nguyen-Xuan et al.

0.000 0.030

0.035

0.040

0.045

0.050 x [m]

(a) nodal distributions represented in Fig. 5. The four analyses assumed the following parameters: dref = 8h and Δa = h . The parameter h varies with the nodal density and it is defined by Eq. (18) considering only the nodes in the vicinity of the crack. For the nodal distributions shown in Fig.  5, the h parameter are, respectively, h = 1.50 mm , h = 0.75 mm , h = 0.75 mm and h = 0.375 mm. The results obtained with the numerical approach presented in this work are shown in Fig.  6. As Fig.  6a

13

0.055

0.060

0.065

0.070

0.0075 0.0575

0.0600 x [m]

0.0625

(b)

shows, the iterative crack propagation initiates in point x0 = {0.0625 , 0.0025} m  . In the graphs presented in Fig. 6, only the iterative position of the crack tip node is shown. The close-up graph presented in Fig. 6b permits to visualise with more detail the differences between the four NNRPIM analyses and the reference solution [9]. From Fig. 6, it is possible to understand that the analyses performed with refined nodal distributions permit to obtain results closer to the ESFEM solution [9].

Engineering with Computers 0.025

0.0175

0.020

0.0150 0.015

y [m]

[0.5h] 0.010

0.005

0.0125

[h]

y [m]

Fig. 7  a Prediction of the crack growth of the four point bending beam, for different values of crack increment Δa , with dref = h . b Close-up detail near the circular opening

[2h]

0.0100

Nguyen-Xuan et al. 0.000 0.030

0.035

0.040

0.045

0.050 x [m]

0.055

0.060

0.065

0.0075 0.0575

0.070

(a) Since the solution obtained with the refined 81 × 21 nodal distribution is capable to produce continuous and smooth crack opening paths, without any abrupt changes in the direction, only the refined 81 × 21 nodal distribution, presented in Fig. 5e, will be used in the further numerical studies performed in this subsection. The next study aims to investigate how the crack increment length parameter Δa influences the prediction of the crack propagation. Thus, the solid domain was discretised with the nodal distribution presented in Fig. 5e. Here, it was assumed a much lower radial distance parameter: dref = h . Due to the configuration of the nodal distribution, the inter-nodal distance near the crack was considered as: h = 0.375 mm . Additionally, three crack increment length parameters Δa were separately considered: Δa = {0.5h, h, 2h}. The crack propagation paths obtained with the NNRPIM formulation are presented in Fig. 7. From Fig. 7a, it is perceptible that the three NNRPIM solutions are very close with the ESFEM solution. Nevertheless, a closer look,

Fig. 8  a Prediction of the crack growth of the four point bending beam for different values of dref  , with Δa = h . b Close-up detail near the circular opening

Fig. 7b, permits to verify that the solution obtained with Δa = h is much closer to the reference solution. Notice that the increase of the value of parameter Δa is directly associated with the increase of the number of eliminated nodes in the vicinity of the crack tip in each iteration, as shown in Fig. 3b. Therefore, assuming high values for Δa lead to the loss of accuracy of the stress field in the vicinity of the crack tip. The crack increment length parameters Δa should be near to Δa = h. In the last study regarding this benchmark example, the behaviour of the proposed numeric approach is analysed for distinct radial distance parameters dref  , as shown in Fig. 3a. As Sect. 3.2 explains, the parameter dref permits to control the number of integration points contributing to the estimated stress state on the crack tip 𝝈(x0 ) , as given in Eq. (15). Thus, in this study, four district radial distance parameters were assumed: dref = {0.5h, h, 2h, 4h} . The crack increment length parameter was kept constant: Δa = h . Since this example was analysed for each dref mentioned using the nodal discretisation presented in Fig. 5e, the value of the inter-nodal distance near the crack is also constant: 0.0175

0.020

0.005

0.0150

[0.5h]

y [m]

y [m]

0.010

0.0625

(b)

0.025

0.015

0.0600 x [m]

[h]

0.0125

[2h] 0.0100

[4h] Nguyen-Xuan et al.

0.000 0.030

0.035

0.040

0.045 0.050 x [m]

(a)

0.055

0.060

0.065

0.070

0.0075 0.0575

0.0600 x [m]

0.0625

(b)

13



h = 0.375 mm . Figure 8 presents the iterative position of the crack tip node of each NNRPIM analysis and the ESFEM solution [9]. As shown in Fig. 8, the variation of the parameter dref has a significant impact on the obtained crack propagation path. It is visible that the NNRPIM formulation using dref = h permits to obtain solutions very close with the ESFEM solution [9]. However, the results obtained with dref ⩾ 2.0h are smoother. To understand the efficiency of the NNRPIM, a similar FEM (2D three nodes constant strain triangle) code was developed to simulate the crack opening of the four point bending plate, as shown in Fig. 5a. The FEM code follows the same architecture of the NNRPIM code. First, the elements and the nodes are established and the shape functions are constructed. Then, the stiffness matrix and the external force vector are built. Afterwards, the essential boundary conditions are imposed using the direct imposition method. Later, the displacement filed is obtained, and consequently, strain and stress fields are achieved. To obtain the stress at the crack tip and the crack propagation, the same technique used for the NNRPIM is applied, see Sect. 3. The methodology to explicitly open the crack is the well-known “discrete cracked element” described in the literature [68], which allows to propagate cracks inside the finite elements. To establish the new finite element mesh (due to the splitting of the cracked element and the division of the adjacent elements ensuring the compatibility of the neighbouring finite elements), this technique requires the implementation of a local remeshing technique combined with an adaptivity method. The exposition of the FEM technique is beyond the present manuscript, therefore, the authors recommend to consult the literature [68]. The results obtained with the developed FEM code overlap the ones obtained by Nguyen-Xuan and coworkers [9]. Thus, such results are not depicted in Figs. 6, 7 and 8. To compare the performance of both FEM and NNRPIM codes, the several computational phases were identified and its individual computational time intervals were documented. The identified phases are indicated and briefly described in Table 1. In Fig. 9a, it is shown the computational time of each one the computational phases. Notice that the NNRPIM possess a significantly higher computational cost in the “preprocessing” and “shape-function” phases, a direct consequence of: the construction of the Voronoï diagram (computationally very expensive); the larger number of integration points; and the higher number of nodes used to construct the shape functions. Furthermore, since the NNRPIM shape functions are built using much more nodes than the FEM (2D three nodes constant strain triangle), it is understandable that the construction of the stiffness matrix (SM) and external

13

Engineering with Computers

force vector (NB) is much more time consuming that the corresponding FEM phases. Solving the equation system constructed using the NNRPIM takes more time than the solution of the FEM equations system (DF phase). This occurs because the bandwidth of the FEM equations system is smaller than the bandwidth of the NNRPIM equations system (once again, a direct consequence of the number of nodes used to construct the shape functions). Obtaining the strain and stress fields with the NNRPIM is more computational expensive, because in the FEM (recall that it was used the constant strain three nodes element) the strain/stress states are obtained at the element level and for the NNRPIM the strain/stress states are obtained for each integration point. For the same nodal discretization, there are much more integration points in the NNRPIM discretization than elements in the FEM discretization. To illustrate the statement, notice that a 2 × 2 nodal discretization leads to two triangular elements. However, the same 2 × 2 nodal discretization leads to eight integration points for the NNRPIM (see [17]). Regarding the MS phase, corresponding to the procedure described Sect. 3 (obtaining the stress state at the crack tip, the corresponding principal stresses and directions and defining the crack propagation directions and length), it is possible to visualise that the NNRPIM is capable to perform this task much faster. This is due to the Voronoï diagram. For the NNRPIM it is not necessary to search the integration points associated (directly or indirectly) to the crack tip, the influence cell concept facilitates the association. In opposition, for the FEM, it is necessary to search all the elements in the neighbourhood of the crack tip. The remeshing (RM) phase is also more time consuming for the FEM. Since a new crack tip node and two crack opening nodes are created at each step, new elements have to be created/split and their conformity/consistency has to be verified and validated. Alternatively, in this phase, for the NNRPIM it is only necessary to create for the new crack tip node and its two crack opening nodes, their influence cells, new integration points and the corresponding shape functions. In Fig.  9b, it shows the computational cost of the NNRPIM and the FEM analyses. At the abscissa, it is represented the number of the crack opening step and at the ordinate is shown the computational time required to achieve it. For both NNRPIM and FEM, it was solved the discretization represented in Fig. 5e and it was considered a crack increment length Δa = h and radial distance parameter dref = h . From Fig. 9b, it is perceptible that the NNRPIM possesses a higher computational cost. Notice in Fig. 9b that the NNRPIM final computational time is about two times higher than the FEM (the NNPRIM takes 12.8 s to finish the analysis and the FEM only takes 5.1 s). However, the slope of the both curves indicates that for larger analyses (with a higher number of crack increments) the FEM will possibly start to show higher total computational costs.

RM

DF SS MS

EB

NB

SF SM

NNRPIM

Establish the nodal mesh and the corresponding Voronoï diagram. Define the natural neighbours (influence cells) and construct the background integration mesh Shape functions Construct the shape functions for each element Construct the shape functions for integration Stiffness matrix Build the local stiffness matrix (element based). Assemble the local stiffness Build the local stiffness matrix (integration point based). Assemble the local matrix into a global stiffness matrix stiffness matrix into a global stiffness matrix Natural boundary Build the local external force vector (element based). Assemble the local Build the local external force vector (integration point based). Assemble the external force vector into a global external force vector local external force vector into a global external force vector Essential boundary Impose the essential boundary conditions (displacement constrains) in the Impose the essential boundary conditions (displacement constrains) in the system of equations system of equations Displacement field Obtain the displacement field Obtain the displacement field Strain/stress field Obtain the strain and stress at each element Obtain the strain and stress at each integration point As described in Sect. 3 Maximum stress Determine the stress state at the crack tip and determine the corresponding principal stresses and directions. Define the crack propagation directions and Determine the stress state at the crack tip and determine the corresponding principal stresses and directions. Define the crack propagation directions and length length ReMeshing Remesh the nodal/element mesh at the crack tip using local remeshing techRepeat phases PP and SF, just for the crack tip nodes nique combined with an adaptivity method

Establish the nodal/element mesh

PP

Pre-processing

FEM

Acronym Phase

Table 1  Acronyms and descriptions of the processing phases in both numerical methods: FEM and NNRPIM

Engineering with Computers

13

10 NNRPIM

10

FEM-3n

Computaitonal cost (s)

1 Computational cost (s)

Fig. 9  a Computational cost of each one of the computational phases. b Complete (accumulated) computational cost of the analysis for each crack opening step

Engineering with Computers

0.1

1

0.01

NNRPIM FEM-3n

0.001

PP

SF

SM

NB

EB

DF SS MS RM Computational phase

0.1

1

(a)

(b)

10 Crack opening step

4.2 Three point bending beam

Fig. 10  Generic representation of a three point bending beam with three circular openings and an initial crack

Fig. 11  Nodal meshes used in the analysis: a 41 × 17 uniform mesh (739 nodes); b 41 × 17 locally refined mesh (769 nodes); c 81 × 33 uniform mesh (2757 nodes); d 81 × 33 locally refined mesh (4797 nodes)

13

As the previous example has shown, the presence of holes in plates disturbs the stress/strain fields, leading to complex curvilinear crack growth trajectories. As the literature shows [12, 56, 69], the crack opening paths found experimentally vary significantly with the initial crack position and geometry. Thus, this section aims to analyse the example schematically represented in Fig. 10 for two distinct initial crack positions. As Fig. 10 shows, the crack initial position depends on the geometric parameters: a and b . The a parameter denotes the initial length of the crack opening. The parameter b indicates the position of the beginning of the crack opening along the Ox axis. Since the initial crack opening is perfectly aligned with the Oy axis, no other geometric information is required to fully characterize this example. In this section, all the studies are analysed assuming the plane strain deformation theory.

5

5

4

4

3

3

41x17

2

41x17 (refined) 81x33

1 0

y [cm]

y [cm]

Engineering with Computers

2 [0.5h] [h] [2h] Bittencourt et al.

1

81x33 (refined) Bittencourt et al. 3

4

5

6 x [cm]

7

8

0

9

3

4

5

6 x [cm]

7

8

9

Fig. 12  Prediction of the crack growth on the first specimen, for different nodal meshes. dref = 8h and Δa = h

Fig. 13  Prediction of the crack growth on the first specimen, for different values of crack increment Δa , with dref = 2h

The present study aims to continue the investigation regarding the influence on the predicted crack opening path of the following variables: the level of the nodal discretisation; the variation of the crack increment length parameter Δa ; and the variation of the radial distance parameter dref .

subsection, only the refined 81 × 21 nodal distribution, Fig. 11d, is considered. In the second study, the proposed problem is studied for three distinct crack increment length parameters: Δa = {0.5h, h, 2h} . The other parameters required by the crack opening algorithm are the following: dref = 2h and h = 0.125 cm. The NNRPIM solutions obt ained for each Δa = {0.5h, h, 2h} are shown in Fig. 13. The results show that using a crack increment length parameter Δa = h leads

In the first studied example, the crack initial geometric parameters indicated in Fig. 10 are considered as: a = 1 cm and b = 4 cm . The influence of the domain discretisation is studied first. Thus, the physical domain is discretized with the four nodal distributions as shown in Fig. 11. In this first study, the following parameters are considered: dref = 8h and Δa = h . The parameter h is obtained with Eq. (18). Hence, respectively, for each nodal mesh presented in Fig.  11: h = 0.50 cm , h = 0.25 cm , h = 0.25 cm and h = 0.125 cm. The results obtained with the crack propagation algorithm proposed in this work are presented in Fig. 12 for each nodal mesh represented in Fig. 11. The predicted crack opening paths are compared with the FEM solution proposed by Bittencourt et al. [56]. In the original paper of Bittencourt et al. [56], the crack propagation path predicted with the FEM is compared with an experimental test. Since the best FEM solution and the experimental solution are almost coincident [56], in this work, the crack opening path obtained with the NNRPIM formulation is compared only with the best FEM result in [56]. As Fig. 12 shows, the results obtained with the NNRPIM formulation are very close to each other, regardless the considered nodal distribution. For the following studies of this

5 4

y [cm]

4.2.1 Study 1

3 [h]

2

[2h] [4h]

1 0

[8h] Bittencourt et al. 3

4

5

6 x [cm]

7

8

9

Fig. 14  Prediction of the crack growth on the first specimen, for different values of dref  , with Δa = h

13

4.2.2 Study 2 The second study regarding the three point bending beam respects the same configuration presented in Fig. 10 and the following geometric parameters are assumed: a = 1.5 cm and b = 5 cm . The crack propagation path verified experimentally with this geometric configuration is much more complex than the crack opening path of the example in Sect. 4.2.1. This benchmark example can be found in the literature [12, 69]. The work of Ingraffea and Grigoriu [69] permits to reproduce with detail the experimental crack curve found for this example. Additionally, this example was analysed with the XFEM by Geniaut and Galenne [12]. Thus, the crack path predicted with the NNRPIM formulation is compared with these last reference solutions. Once again, first, the influence of the nodal discretisation is investigated. Thus, the problem domain Ω ∈ ℝ2 is discretised in the four nodal distributions presented in Fig. 15. The numeric prediction of the crack growth is a non-linear procedure, which depends on the accuracy of the stress field in the vicinity of the crack tip. Therefore, it is important Fig. 15  Nodal meshes used in the second specimen: a 61 × 25 uniform mesh (1527 nodes); b 61 × 25 (refined 1) locally refined mesh (2141 nodes); c 61 × 25 (refined 2) uniform mesh (3607 nodes); d 81 × 33 (refined 1) locally refined mesh (3961 nodes)

13

5

4

4

3.5

61x25

3

y [cm]

to a solution almost coincident with the crack propagation path presented in [56]. Next, the same nodal distribution presented in Fig. 11d is considered to analyse the variation of the distinct radial distance parameters dref  . Here, the crack increment is kept constant: Δa = h , with h = 0.125 cm , and the four district radial distance parameters are assumed: dref = {h, 2h, 4h, 8h} . The crack growths predicted with the NNRPIM formulation are shown in Fig. 14 and compared with the FEM simulation [56]. As Fig. 14 shows, the variation of the dref affects significantly the obtained solution. The results on Fig. 14 show that when dref = 2h is considered, the NNRPIM solution is coincident with the FEM result [56]. However, dref = 4h is capable to produce acceptable results.

Engineering with Computers

y [cm]



61x25 (refined 1)

3

61x25 (refined 2) 2

81x33 (refined 1)

2.5

Geniaut et al. Ingraffea et al. 1

5

6

x [cm]

(a)

7

8

2 5.25

5.75 x [cm]

6.25

(b)

Fig. 16  a Prediction of the crack growth on the second specimen, for different nodal meshes. dref = 2h and Δa = h . b Close-up detail near the circular opening

to increase the level of the nodal discretisation on the area of the expected crack opening path. The discretisation presented in Fig. 15c assures a nodal density four times higher than the rest of the domain. Additionally, notice that the nodal meshes shown in Fig. 15c, d possess approximately the same nodal density on the area of the expected crack path. However, globally, the nodal discretisation of Fig. 15c has less nodes that the nodal distribution shown in Fig. 15d, contributing to the efficiency of the analysis. Regarding the crack growth algorithm, in this study, the following parameters were considered: a radial distance dref = 2h ; a crack increment length Δa = h ; and a parameter h obtained with Eq. (18). Respectively, for each mesh presented in Fig.  15, the following inter-nodal distance near the crack were obtained: h = 0.333 cm , h = 0.167 cm , h = 0.125 cm and h = 0.083 cm.

5

4

4

3.5

4

3.5

3 [0.5h] [h]

2

3

3 [2h] [4h]

2

2.5

[2h]

Geniaut et al.

Ingraffea et al. 5

6

x [cm]

7

2 5.25

8

(a)

1

5.75 x [cm]

6.25

4

4

3.5

y [cm]

y [cm]

5

3 [0.5h] [0.75h]

3

2.5

[h] Geniaut et al.

1

Ingraffea et al. 5

6

x [cm]

(a)

7

8

2 5.25

5.75 x [cm]

Ingraffea et al. 5

6

x [cm]

(a)

(b)

Fig. 17  a Prediction of the crack growth on the second specimen, for different values of crack increment Δa , with dref = 2h . b Close-up detail near the circular opening

2

3

2.5

[8h]

Geniaut et al. 1

y [cm]

4

y [cm]

5

y [cm]

y [cm]

Engineering with Computers

6.25

(b)

Fig. 18  a Prediction of the crack growth on the second specimen, for different values of dref = {0.5h, 0.75h, h} , with Δa = h . b Close-up detail near the circular opening

The crack growth predictions obtained for each nodal discretisation are shown in Fig. 16. With the close-up graphs presented in Fig. 16b is easier to visualise the differences between the obtained solutions and the reference solutions [12, 69]. Figure 16b permits to observe that the nodal distributions presented in Fig. 15c, d are capable to capture with a higher accuracy the complex behaviour of the crack opening curve near the first circular opening. Additionally, with these two nodal distributions it is possible to approach the experimental crack propagation path obtained by Ingraffea and Grigoriu [69]. Thus, only the nodal discretization presented in Fig. 15d is assumed in the following analyses. In the second study, the proposed problem is studied for three distinct crack increment length parameters:

7

8

2 5.25

5.75 x [cm]

6.25

(b)

Fig. 19  a Prediction of the crack growth on the second specimen, for different values of dref = {2h, 4h, 8h} , with Δa = h . b Close-up detail near the circular opening

Δa = {0.5h, h, 2h} . The other parameters required by the crack opening algorithm are the following: dref = 2h and h = 0.083 cm . The NNRPIM solutions obtained for each Δa = {0.5h, h, 2h} are shown in Fig. 17. The results show that only with Δa = h it is possible to obtain a crack opening path similar with the reference solutions. The remaining values, Δa = {0.5h, 2h} , lead to inaccurate crack propagation paths. Thus, it is assumed that for the numerical approach proposed in this work, the value Δa = h can be used globally to predict the crack opening path. Next, a final study is performed. Here, the influence of the radial distance parameter dref on the obtained solution is investigated, Fig. 3a. First, the problem domain is discretized with the nodal mesh presented in Fig. 15d. Then, it is considered a crack increment Δa = h and an inter-nodal distance near the crack: h = 0.83 cm . Afterwards, the problem is analysed considering six district radial distance parameters: dref = {0.5h, 0.75h, h, 2h, 4h, 8h}. For each NNRPIM analysis, the iterative positions of the crack tip node are presented in Figs. 18 and 19. Once more, the NNRPIM solutions are compared with the experimental solution of Ingraffea and Grigoriu [69] and the XFEM solution of Geniaut and Galenne [12]. As in the previous examples, the results show that the numeric prediction of the crack growth is highly influenced by the variation of the radial distance parameter. From Figs. 18 and 19, it is possible to visualise that the NNRPIM solution considering dref = h permits to obtain a crack opening path very close with the experimental solution [69] and if dref = 8h is considered, the results are closer to the XFEM solution [12].

13



5 Conclusion In this work, an algorithm to predict the crack growth was developed and combined with the NNRPIM—a truly meshless method. The developed numeric approach is fully automatic. Being this a preliminary study, the algorithm considers a linear elastic model and the maximum circumferential stress criterion. First, an elastostatic analysis is performed and the correspondent displacement, strain and stress fields are obtained. Then, the stress field obtained in the vicinity of the crack tip is used to calculate an estimated stress state at the crack tip. With this approximated stress state, it is possible to obtain the respective principal stresses and principal directions. Afterwards, the direction of the crack propagation is defined as orthogonal to the principal stress direction. Subsequently, the crack length is increased in a straight line segment respecting the obtained direction and a nodal remeshing is performed in the vicinity of the crack tip. Then, a new elastostatic analysis is executed and the process restarts. As mentioned before, the new crack tip position depends on the stress field in the vicinity of the crack tip and the update of the location of the new crack tip requires a local remeshing. Therefore, the inclusion of the NNRPIM in a crack opening algorithm is a strong asset, because the NNRPIM is capable to perform automatically a nodal remeshing without the user intervention [17] and capable to produce high quality stress/strain fields [17]. The crack propagation algorithm here presented requires from the user the introduction of two parameters: the radius of selection of integration points dref and the crack increment length Δa . The predicted crack propagation path depends on the value of these two parameters. The parameter dref defines the number of integration points considered to calculate the propagation direction of the crack, as shown in Fig. 3a. If dref assumes high values, the estimated stress state at the vicinity of the crack tip will not be a correct representation of the actual stress state in that area, because the stress state of distant points will be included in the approximation. If dref takes very low values, incorrect oscillations appear in the crack path. From the examples studied in this work, it was found that the dref parameter should assume the following values: dref ∈ ]h, 2h[. The crack increment length Δa defines the number of nodes that discretize the crack opening edges. Thus, this parameter strongly influences the crack propagation path. Since the nodal discretisation influences significantly the quality of the constructed shape functions, the abrupt variation of the nodal distances are responsible for shape functions with low quality, resulting in inaccurate results. Thus, the crack increment length should assume values similar to the inter-nodal distance near the crack tip: h . This remark

13

Engineering with Computers

is corroborated with the obtained results. In all the studied examples, Δa = h permitted to obtain the most accurate crack propagation path. In this work, several challenging benchmark crack propagation examples were analysed with the developed meshless approach. The results obtained in this preliminary study permit to conclude that the numeric tools here presented are suitable to efficiently analyse fracture mechanic problems. Acknowledgements  The authors truly acknowledge the funding provided by Ministério da Ciência, Tecnologia e Ensino Superior—Fundação para a Ciência e a Tecnologia (Portugal), under Grants SFRH/BPD/75072/2010 and SFRH/BPD/111020/2015, and by project funding UID/EMS/50022/2013 (funding provided by the inter-institutional projects from LAETA). Additionally, the authors gratefully acknowledge the funding of Project NORTE-01-0145FEDER-000022—SciTech—Science and Technology for Competitive and Sustainable Industries, cofinanced by Programa Operacional Regional do Norte (NORTE2020), through Fundo Europeu de Desenvolvimento Regional (FEDER).

References 1. Aliabadi MH, Rooke DP (1991) Numerical fracture mechanics. Computational Mechanics, Southampton 2. Watwood VB (1970) The finite element method for prediction of crack behavior. Nucl Eng Des 11(2):323–332 3. Boljanović S, Maksimović S (2011) Analysis of the crack growth propagation process under mixed-mode loading. Eng Fract Mech 78(8):1565–1576 4. Tong P, Pian THH, Lasry SJ (1973) A hybrid-element approach to crack problems in plane elasticity. Int J Numer Methods Eng 7(3):297–308 5. Karihaloo BL, Xiao QZ (2001) Accurate determination of the coefficients of elastic crack tip asymptotic field by a hybrid crack element with p-adaptivity. Eng Fract Mech 68(15):1609–1630 6. Leung AYT, Su RKL (Aug. 1995) Mixed-mode two-dimensional crack problem by fractal two level finite element method. Eng Fract Mech 51(6):889–895 7. Moes N, Dolbow J, Belytschko T (1999) A finite element method for crack growth without remeshing. Int J Numer Methods Eng 46(1):131–150 8. Vu-Bac N, Nguyen-Xuan H, Chen L, Bordas S, Kerfriden P, Simpson RN, Liu GR, Rabczuk T (2011) A node-based smoothed extended finite element method (NS-XFEM) for fracture analysis. Comput Model Eng Sci 73(4):331–355 9. Nguyen-Xuan H, Liu GR, Nourbakhshnia N, Chen L (2012) A novel singular ES-FEM for crack growth simulation. Eng Fract Mech 84:41–66 10. Nguyen-Xuan H, Liu GR, Bordas S, Natarajan S, Rabczuk T (2013) An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order. Comput Methods Appl Mech Eng 253:252–273 11. Bordas SPA, Rabczuk T, Hung N-X, Nguyen VP, Natarajan S, Bog T, Quan DM, Hiep NV (2010) Strain smoothing in FEM and XFEM. Comput Struct 88(23–24):1419–1443 12. Geniaut S, Galenne E (2012) A simple method for crack growth in mixed mode with X-FEM. Int J Solids Struct 49(15–16):2094–2106 13. Fries T-P, Baydoun M (2012) Crack propagation with the extended finite element method and a hybrid explicit-implicit crack description. Int J Numer Methods Eng 89(12):1527–1558

Engineering with Computers 14. Mariani S, Perego U (2003) Extended finite element method for quasi-brittle fracture. Int J Numer Methods Eng 58(1):103–126 15. Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P (1996) Meshless methods: an overview and recent developments. Comput Methods Appl Mech Eng 139(1–4):3–47 16. Nguyen VP, Rabczuk T, Bordas S, Duflot M (2008) Meshless methods: a review and computer implementation aspects. Math Comput Simul 79(3):763–813 17. Belinha J (2014) Meshless methods in biomechanics—bone tissue remodelling analysis. Lecture notes in computational vision and biomechanics, vol 16. Springer, Netherlands 18. Gingold RA, Monaghan JJ (1977) Smoothed particle hydrodynamics—theory and application to non-spherical stars. Mon Not R Astron Soc 181:375–389 19. Monaghan J (2005) Smoothed particle hydrodynamics. Rep Prog Phys 68(8):1703–1759 20. Nayroles B, Touzot G, Villon P (1992) Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 10(5):307–318 21. Belytschko T, Lu YY, Gu L (1994) Element-free Galerkin methods. Int J Numer Methods Eng 37(2):229–256 22. Liu WK, Jun S, Zhang YF (1995) Reproducing kernel particle methods. Int J Numer Methods Fluids 20(8–9):1081–1106 23. Atluri SN, Zhu T (1998) A new Meshless Local Petrov-Galerkin (MLPG) approach in computational mechanics. Comput Mech 22(2):117–127 24. Sukumar N, Moran B (1998) Belytschko T, The natural element method in solid mechanics. Int J Numer Methods Eng 43(5):839–887 25. Sukumar N, Moran B, Semenov AY, Belikov VV (2001) Natural neighbour Galerkin methods. Int J Numer Methods Eng 50(1):1–27 26. Liu GR, Gu YT (2001) A point interpolation method for two-dimensional solids. Int J Numer Methods Eng 50(4):937–951 27. Wang JG, Liu GR (2002) A point interpolation meshless method based on radial basis functions. Int J Numer Methods Eng 54(11):1623–1648 28. Dinis LMJS, Jorge RMN, Belinha J (2007) Analysis of 3D solids using the natural neighbour radial point interpolation method. Comput Methods Appl Mech Eng 196(13–16):2009–2028 29. Belinha J, Dinis LMJS, Jorge RMN (2013) The natural radial element method. Int J Numer Methods Eng 93(12):1286–1313 30. Onãte E, Idelsohn S, Zienkiewicz OC, Taylor RL (1996) A finite point method in computational mechanics—applications to convective transport and fluid flow. Int J Numer Methods Eng 39(22):3839–3866 31. Gu YT (2005) Meshfree methods and their comparisons. Int J Comput Methods 2(4):477 32. Dinis LMJS, Jorge RMN, Belinha J (2008) Analysis of plates and laminates using the natural neighbour radial point interpolation method. Eng Anal Bound Elem 32(3):267–279 33. Wang JG, Liu GR (2002) On the optimal shape parameters of radial basis functions used for 2-D meshless methods. Comput Methods Appl Mech Eng 191(23–24):2611–2630 34. Sibson R (2008) A vector identity for the Dirichlet tessellation. Math Proc Cambridge Philos Soc 87(1):151–155 35. Dinis LMJS, Jorge RMN, Belinha J (2010) A 3D shell-like approach using a natural neighbour meshless method: isotropic and orthotropic thin structures. Compos Struct 92(5):1132–1142 36. Dinis LMJS, Jorge RMN, Belinha J (2010) Composite laminated plates: a 3D natural neighbor radial point interpolation method approach. J Sandw Struct Mater 12(2):119–138 37. Dinis LMJS, Jorge RMN, Belinha J (2010) An unconstrained thirdorder plate theory applied to functionally graded plates using a meshless method. Mech Adv Mater Struct 17(2):108–133

38. Dinis LMJS, Jorge RMN, Belinha J (2009) The natural neighbour radial point interpolation method: dynamic applications. Eng Comput 26(8):911–949 39. Dinis LMJS, Jorge RMN, Belinha J (2011) Static and dynamic analysis of laminated plates based on an unconstrained third order theory and using a radial point interpolator meshless method. Comput Struct 89(19–20):1771–1784 40. Dinis LMJS, Jorge RMN, Belinha J (2011) A natural neighbour meshless method with a 3D shell-like approach in the dynamic analysis of thin 3D structures. Thin Wall Struct 49(1):185–196 41. Dinis LMJS, Jorge RMN, Belinha J (2009) Large deformation applications with the radial natural neighbours interpolators. Comput Model Eng Sci 44(1):1–34 42. Dinis LMJS, Jorge RMN, Belinha J (2009) The radial natural neighbours interpolators extended to elastoplasticity. In: Ferreira AJM, Kansa EJ, Fasshauer GE, Leitão VMA (eds) Progress on meshless methods—computational methods in applied sciences. Springer, Netherlands, pp 175–198 43. Belinha J, Jorge RMN, Dinis LMJS (2012) Bone tissue remodelling analysis considering a radial point interpolator meshless method. Eng Anal Bound Elem 36(11):1660–1670 44. Belinha J, Dinis LMJS, Jorge RMN (2013) A meshless microscale bone tissue trabecular remodelling analysis considering a new anisotropic bone tissue material law. Comput Methods Biomech Biomed Eng 16(11):1170–1184 45. Moreira SF, Belinha J, Dinis LMJS, Jorge RMN (2014) A global numerical analysis of the ‘central incisor/local maxillary bone’ system using a meshless method. MCB Mol Cell Biomech 11(3):151–184 46. Belinha J, Dinis LMJS, Jorge RMN (2015) The mandible remodelling induced by dental implants: a meshless approach. J Mech Med Biol 15(4):1550059 47. Rabczuk T, Belytschko T (2007) A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput Methods Appl Mech Eng 196(29–30):2777–2799 48. Bordas S, Rabczuk T, Zi G (2008) Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Eng Fract Mech 75(5):943–960 49. Ventura G, Xu JX, Belytschko T (2002) A vector level set method and new discontinuity approximations for crack growth by EFG. Int J Numer Methods Eng 54(6):923–944 50. Rao BN, Rahman S (2001) A coupled meshless-finite element method for fracture analysis of cracks. Int J Press Vessel Pip 78(9):647–657 51. Duflot M, Nguyen-Dang H (2004) A meshless method with enriched weight functions for fatigue crack growth. Int J Numer Methods Eng 59(14):1945–1961 52. Duflot M (2006) A meshless method with enriched weight functions for three-dimensional crack propagation. Int J Numer Methods Eng 65(12):1970–2006 53. Li SC, Li SC, Cheng YM (2005) Enriched meshless manifold method for two-dimensional crack modeling. Theor Appl Fract Mech 44(3):234–248 54. Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H (2008) A geometrically non-linear three-dimensional cohesive crack method for reinforced concrete structures. Eng Fract Mech 75(16):4740–4758 55. Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H (2010) A simple and robust three-dimensional cracking-particle method without enrichment. Comput Methods Appl Mech Eng 199(37–40):2437–2455 56. Bittencourt TN, Wawrzynek PA, Ingraffea AR, Sousa JL (1996) Quasi-automatic simulation of crack propagation for 2D LEFM problems. Eng Fract Mech 55(2):321–334 57. Erdogan F, Sih GC (1963) On the crack extension in plates under plane loading and transverse shear. J Basic Eng 85(4):519

13

58. Sumi Y (1985) Computational crack path prediction. Theor Appl Fract Mech 4(2):149–156 59. Voronoï G (1908) Nouvelles applications des paramètres continus à la théorie des formes quadratiques. deuxième mémoire. recherches sur les paralléllo èdres primitifs. J fur die reine Angew Math 134:198–287 60. Delaunay B (1934) Sur la sphère vide. A la mémoire de Georges Voronoï. Bull Acad Sci USSR 6:793–800 61. Moreira S, Belinha J, Dinis LMJS, Jorge RMN (2014) Análise de vigas laminadas utilizando o natural neighbour radial point interpolation method. Rev Int Métodos Numéricos para Cálculo y Diseño en Ing 30(2):108–120 62. Belinha J, Dinis LMJS, Jorge RMN (2013) Analysis of thick plates by the natural radial element method. Int J Mech Sci 76:33–48 63. Belinha J, Dinis LMJS, Jorge RMN (2013) Composite laminated plate analysis using the natural radial element method. Compos Struct 103:50–67

13

Engineering with Computers 64. Hardy RL (1990) Theory and applications of the multiquadric-biharmonic method 20 years of discovery 1968–1988. Comput Math with Appl 19(8–9):163–208 65. Golberg MA, Chen CS, Bowman H (1999) Some recent results and proposals for the use of radial basis functions in the BEM. Eng Anal Bound Elem 23(4):285–296 66. Fleming M, Chu YA, Moran B, Belytschko T (1997) Enriched element-free galerkin methods for crack tip fields. Int J Numer Methods Eng 40(8):1483–1504 67. Belytschko T, Lu YY, Gu L (1995) Crack propagation by elementfree Galerkin methods. Eng Fract Mech 51(2):295–315 68. Mohammadi S (2008) Extended finite element method for fracture analysis of structures. Blackwell, Oxford 69. Ingraffea AR, Grigoriu M (1990) Probabilistic fracture mechanics: a validation of predictive capability. Cornell Univ Ithaca NY Dept Struct Eng ADA228877:1–155