PSI - Issue 80

3

D.C. Gonçalves et al. / Procedia Structural Integrity 80 (2026) 443–450 Author name / Structural Integrity Procedia 00 (2019) 000–000

445

propagation, drawbacks such as the definition of the influence domains near the crack tip and interfaces require extra care. Reference (Zhuang et al., 2011) settles a general formulation of meshless methods and the visibility criterion for accurate fracture modelling using meshless techniques. In this work, the natural neighbours radial point interpolation method (NNRPIM) is implemented to predict crack propagation in adhesively bonded joint specimens. In opposition to conventional meshless techniques, the NNRPIM uniquely relies on the nodal discretization to construct the integration points based on the Voronoi diagram and Delaunay triangulation. Geometric crack propagation is explicit, and not hindered by the remeshing of local fixed elements, as in finite element methods. Numerical results demonstrate the suitability of the proposed NNRPIM-based crack propagation algorithm to simulate fracture propagation in adhesively bonded joints. 2. The Natural Neighbours Radial Point Interpolation Method The first step in the NNRPIM is to discretize the problem domain in a set of independent field nodes (Fig. 1a). The Voronoi diagram can then be constructed (Fig. 1b), allowing to define the natural neighbours of each node (Belinha, 2014). Then, Voronoi cells are divided in integration subcells (Fig. 1c), and integration points are distributed considering the constructed integration subcells (Fig. 1d). Differently from classic finite element techniques, nodal connectivity is established by the overlap of influence domains. In the NNRPIM, the influence domain of a given integration point is defined by the natural neighbours of the nearst field node. Due to the nature of the influence domain being the Voronoi cells, in the NNRPIM formulation the influence domains are usually described as influence cells. With the integration points and influence cells established, the following step is the shape functions calculation. The interpolation function ( ! ) is defined by the combination of a Radial Basis Function (RBF), ( ! ) , and a Polynomial Basis Function (PBF), ( ! ) : ( ! )=( " ( ! ) " ( ! ) # "$% +( & ( ! ) & ( ! ) ' &$% , (1) being the number of nodes in the influence domain of ! and the number of monomials in the PBF. Although several RBF are available in literature, this work considers the multiquadrics (MQ-RBF) (Hardy, 1990), "& =1 " ( & + ( 4 ) , (2) where the Euclidean distance, "& , is the only variable in the RBF; and and are shape parameters of the RBF. For the present NNRPIM formulation, these parameters were previously determined and optimised by solving several 2D tests (Belinha, 2014), being = 0.0001 and = 0.9999 the optimal parameters used in this work. The coefficients " ( ! ) and & ( ! ) , in equation (1), are determined by establishing the following equation for each node inside the influence domain of ! : { } < ? = * , (3) being * the vector of the corresponding nodal values. To obtain a set of + equations with + unknowns, an extra set of equations, = , is added to the system (Golberg et al., 1999), which yields: < CDED F ? < ? = < ? . (4) Solving equation (4) permits to define the interpolation function as: ( ! ) = { ( ! ) . ( ! ) . } /% < ? = { ( ! ) . ( ! ) . }< ? , (5) where ( ! ) . is a by-product of equation (5) that is negligible (Belinha, 2014); thus, ( ! ) = ( ! ) . In addition, the shape functions’ spatial partial derivatives can then be easily computed as detailed by Belinha (Belinha, 2014). The procedure to obtain the final set of discrete equations, which is established in the Galerkin weak form, is analogous to the FEM. The equilibrium equations are given by ∇ + = , with = ̅ on the natural boundary Γ 0 , and = R on the essential boundary Γ 1 , being the Cauchy stress tensor, the body force per unit volume vector, the unit outward vector normal to Γ 0 , ̅ the traction force on Γ 0 , and R the prescribed displacement on Γ 1 . The Galerkin weak form for the present linear-static problem is defined as: