Fracture Toughness Determination on an SCB Specimen by Meshless Methods

Abstract

This work investigates fracture characteristics of a marble semi-circular bend (SCB) specimen with a pre-defined crack under a compressive loading condition. It aims at evaluating how the fracture toughness can be affected by the crack and span length variation. Numerically, the model is solved using meshless methods, extended to the linear elastic fracture mechanics (LEFM), resorting to radial point interpolation method (RPIM) and its natural neighbor versions (NNRPIMv1 and NNRPIMv2). Alternatively, to validate the meshless method results, the problem is resolved following the finite element method (FEM) model based on the standard 2D constant strain triangle elements. As a result, fracture toughness and the critical strain energy release rate are characterized following the testing method on the cracked straight through semi-circular bend specimen (CSTSCB). A comparison is drawn amongst the theoretical, meshless methods and FEM results to evaluate the capability of advanced numerical methods. Encouraging results have been accomplished leading to validate the supporting numerical methodologies.

1. Introduction

The capability to endure a substantial amount of damage is a demand for current engineering structures; therefore, it has become progressively imperative to improve methodologies to anticipate failure phenomenon in damaged components. Damage tolerance analyses can be fulfilled within the linear elastic fracture mechanics (LEFM) concepts where the stress intensity factor (SIF) plays a substantial role [1]. The fracture mechanics theorem in conjunction with crack growth laws, i.e., Paris’ law, is characteristically engaged to examine and forecast crack growth and fracture behavior of structural integrity components [2]. To study crack growth and to evaluate the remaining lifetime of a structural component, demanding computational analyses must be performed and the SIFs shall be thereby evaluated. There exist several developed formulations to determine the fracture parameters such as mode I and mixed-mode SIFs, strain energy release rate, T-stress, etc.

The fracture toughness, which is the critical SIF of a sharp crack where propagation of the crack suddenly becomes unlimited and fast, remains a very imperative specification of a material that can demonstrate the material resistance against the structural cracks governing by the work required to fracture a material, known as work of fracture. It meant that fracture toughness is a quantitative definition to express the material’s resistance to a brittle fracture phenomenon in the presence of cracks. If the material possesses high fracture toughness, it would be more disposed to the ductile fracture [3]. On the other hand, brittle fracture is a characteristic of materials with less fracture toughness. There exist general factors which can affect the material toughness including strain rate, temperature, presence of any stress concentration on the specimen surface and the relationship between the material strength and its ductility [4]. The critical strain energy release rate describes the resistance of the material against the crack growth.

In this regard, Almeida-Fernandes et al. [5] developed numerical damage models calibrated with experimentally based fracture toughness parameters on the pultruded glass fiber reinforced polymers. Taylor et al. [6] studied the fracture strength of brittle materials by the effect of stress concentrations through SIF calculation. Mínguez [7] studied the fracture toughness of a center cracked plate by the Finite Element Method (FEM) in which the critical SIF variation through specimen thickness was characterized.

Amongst all testing specimens in the fracture mechanics studies, the semi-circular bend (SCB) specimen received considerable attention emerging in 1984 to test brittle materials in particular rocks contributing to determining the mode I fracture toughness [8]. Chong et al. [9] proposed the SCB specimen in 1987, due to the method of rock and geometrical specimens’ extraction being in a circular shape, simple geometry and demanding a common loading configuration. Subsequently, the SCB was adopted and improved on distinct solicitations in the field of solid and fracture mechanics. The SCB specimens can be simply manufactured to be adopted in the experimental tests and numerical simulations. There are other suggested testing methods for mode I fracture toughness of rock by the International Society for Rock Mechanics (ISRM) such as short rod specimen method (SR) [10], cracked notched Brazilian disk (CCNBD) [11], chevron-bend specimen method (CB) [12] and cracked straight through semi-circular bend specimen (CSTSCB) [13].

In this study, several numerical methods are used to characterize the mode I fracture toughness near the crack tip due to the remote load. The FEM and meshless methods including radial point interpolation method (RPIM) and natural neighbor radial point interpolation method (NNRPIM) are employed to determine the fracture toughness to be compared with the calculated values by CSTSCB method [13].

The necessity of mesh refinement in the problems that the geometry changes or contains a large distortion, FEM shows some shortcomings that leads researchers to emerge new techniques. One of the earliest methods proposed by Belytschko et al. [14] so-called Element-Free Galerkin (EFG) is capable to simulate cracks propagating statically and dynamically. Furthermore, the EFG is able to give valid values for irregular nodal distribution, which is one of the advantages of this technique. Muthu et al. [15] proposed a variant of the EFG method integrating the application of the level set and diffraction methods to simulate multiple interacting cracks using a coarse meshless nodal discretization. This methodology was deployed to model kinked cracks that have knee singularity.

On the other hand, one of the main factors affecting the acceptance of a numerical technique is computational cost. Generally, the cost is measured for a prescribed accuracy for a problem. Therefore, by taking into account that denser discretizations would increase the computational cost, it is necessary to optimize the cost meticulously without a negative effect on accuracy. In this regard, meshless methods are flexible to be optimized and may be adjusted to produce reliable results. Besides, it is worth mentioning that meshless methods provide the opportunity to impose the displacement constraints on the stiffness matrix directly, but the EFG method employs the Lagrange-multipliers for this purpose which is not as straightforward as it is for meshless methods. Meshless methods have been successfully implemented to address fracture and damage characterizations for a variety of material behaviors, c.f. [16,17,18,19,20].

So far, the FEM formulation with quadratic elements has been employed to calculate mode I and mixed-mode SIFs for SCB specimens [9,13,21,22]. Nevertheless, meshless methods have never been employed to calculate fracture toughness and critical strain energy release rate for SCB specimens.

It must be mentioned that the meshless methods are advantageous due to their simplicity in the formulations. Besides, it is feasible to use the FE mesh background to generate the nodal discretizations, which may lead to reducing the computational costs with the absence of finite elements. Other privileges of meshless methods can be related to their independence to the element shapes/types and the element distortion that may occur during the finite element analyses. Despite FEM, meshless methods do not require FE remeshing which is undesirable in terms of crack propagation problems and solution convergences on the nonlinearity in the engineering structures with complex material behaviors such as damage mechanics or porous models, c.f. [16,17,19,20,23,24].

Nevertheless, owing to the aforementioned related works, it can be inferred that the literature lacks a comprehensive analysis of SCB specimens for fracture toughness characterization by meshless methods. Therefore, the main contribution of this work can be related to the development of meshless method formulations combined with LEFM theory. Therefore, it contributed to elastostatically simulating the cracked domain in order to attain the fracture toughness and critical strain energy release rate.

2. Overview on Meshless Methods

In this study, the 2D deformation theory in plane stress state was taken into account for all computational analyses. The RPIM and NNRPIM meshless method formulations were taken into account [25] extended to the LEFM theory. Generally, meshless methods discretize the problem domain and respective boundaries using a nodal distribution. This nodal set cannot be categorized as a mesh, because no previous information regarding the relationship between each node is required to build the interpolation functions for the unknown variational fields [26].

The RPIM is an interpolator meshless method assuring the nodal connectivity within the influence-domain conception, which is chiefly adopted by many meshless techniques because of its simplicity. Similar to the FEM, RPIM meshless method stands as a discrete computational approach. The main difference could be related to the discretization procedure. The elements and nodes discretize the Finite Element (FE) problem domain whereas meshless methods discretize the problem domain only by nodal points. Likewise, in the FEM, the nodal connectivity must be enforced in the pre-processing phase through a predefined FE mesh establishment. Hence, any element owns the intrinsic and adequate information to impose nodal connectivity. Despite the FEM, meshless methods do not possess predefined connectivity mesh, it is essential to generate one [27]. Therefore, generally, the nodal interdependency must be enforced by the influence-domain geometrical construction in meshless methods. It can be acquired once the nodal discretization is created [25]. This is an overlap associated with the influence-domains allowing setting up nodal connectivity. Considering an interest point in the problem domain, it is possible to establish its influence-domain by a radial search centered in that interest point. It led to collecting a definite number of nodes inside a fixed area or volume, the 2D and 3D problems, respectively. The literature proved that this approach is simple to figure out and to implement contributing to support the expansion of several meshless techniques.

Regarding RPIM meshless method, the literature [25] recommends that any 2D influence-domain must include roughly $n=\left[9,16\right]$ nodes. Besides, a background integration mesh is necessary to mathematically integrate the integro-differential equations governing the physical phenomenon [26,28] due to the adaptation of the Galerkin weak-formulation.

In the case of the NNRPIM, the natural neighbor’s concept has been taken into account to impose the nodal connectivity as recommended by Sibson [29]. This mathematical notion could be described by the Voronoï diagram associated with the discretized domain. In the NNRPIM, the problem domain$\mathsf{\Omega }\in {ℝ}^d$, bounded by a physical boundary$\mathsf{\Gamma }\subset \mathsf{\Omega }$, is discretized in numerous arbitrarily distributed points $\mathit{N}=\left\{n_0,n_1,\dots ,n_N\right\}$ scattered in the space domain:$\mathit{X}=\left\{x_1,x_2,\dots ,x_N\right\}\in \mathsf{\Omega }$. The Voronoï diagram of $\mathit{N}$ is the domain partition determined by $\mathsf{\Omega }$ in sub-regions ${\mathit{V}}_i$, closed and convex where any sub-region ${\mathit{V}}_i$ is attendant with the node $n_i$, in a way that any point in the interior of the ${\mathit{V}}_i$ is closer to $n_i$ than any other node $n_j$, where $n_j\in \mathit{N}\wedge j\ne i,{\mathit{x}}_I\in \mathsf{\Omega }$.

$${\mathit{V}}_i:=\left\{{\mathit{x}}_I\in \mathsf{\Omega } \subset {ℝ}^d:\Vert {\mathit{x}}_I-{\mathit{x}}_i\Vert <\Vert {\mathit{x}}_I-{\mathit{x}}_j\Vert ,\forall i\ne j\right\}$$

(1)

However, a set of Voronoï cells${\mathit{V}}_i$defines the Voronoï diagram,$\mathit{V}=\left\{V_1,V_2,\dots ,V_N\right\}$. In the literature [25,30,31], it is possible to find relevant works properly addressing the Voronoï construction procedure.

In relation to the nodal connectivity in NNRPIM, it is made by a set of nodes in the neighborhood of an interest point${\mathit{x}}_i\in \mathit{X}$, the NNRPIM influence-cell notion contributes to forming the node connectivity. It contributes to naturally controlling the influence domain associated with${\mathit{x}}_i$. Meanwhile, only the determination of the 2D influence-cell is given, although this concept is applicable to n-dimensional space. In general, there are two types of influence-cells: the first-degree influence-cell and the second-degree influence-cell. The former one is composed of the first natural neighbors of ${\mathit{x}}_i$ and the latter one includes the natural neighbors of the nodes belonging to the first-degree influence-cell of node ${\mathit{x}}_i$ in addition to the first natural neighbors of${\mathit{x}}_i$. This concept is well described in detail in [23].

The second-degree influence-cell enforces higher nodal connectivity when compared with the first-degree influence-cell. The literature [23] states that regardless of the studied phenomenon, higher degree influence-cells permit to achieve more accurate solutions.

The standard FEM, RPIM and NNRPIM formulations are comprehensively demonstrated in [25] leading to a linear equations system exhibited as $\mathit{K}\mathit{u}=\mathit{f}$. In which, $\mathit{K}$stands as the stiffness matrix, $\mathit{f}$ is the force vector and $\mathit{u}$ is the displacement field. In accordance with Hooke’s law, it is feasible to correlate a relationship between the strain and stress fields as,

$$\mathit{\sigma }=\mathit{C}\mathit{\epsilon }\iff \left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\tau }_{xy}\end{array}\right\}=\frac{E}{\left(1+\nu \right)\left(1-\nu \right)}\left[\begin{array}{ccc}1 & \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\gamma }_{xy}\end{array}\right\}.$$

(2)

In this study, all integration cells are quadrilateral comprising of around 9 nodes and $n_Q\times n_Q$ integration points inside, respecting the Gauss-Legendre quadrature scheme. The literature [23] suggests that this integration scheme aggrandizes the efficiency if $n_Q=3$.

3. Mathematical Formulation

As shown in Figure 1, a cracked straight through semi-circular bend specimen (CSTSCB) is considered as the case study following the testing method suggested by Kuruppu et al. [13]. It aims at the determination of the mode I fracture toughness for various crack growth stages.

From the theoretical point of view, Equation (3) is applicable for the computation of fracture toughness for mode I [13],

$$K_{Ic}=\frac{P_{ }\sqrt{\pi a}}{2RB}{Y}_I.$$

(3)

In which $P_{ }$ is the maximum load, $R$ is the specimen radius, $B$ is the specimen thickness, and $Y_I$ is the normalized SIF, which is dimensionless and dependent on the normalized crack length and normalized span. By considering a thin specimen, $K_{Ic}$ in Equation (3) provides the critical SIF value. Regarding the normalized term, $Y_I$, Lim et al. [32] used the classical displacement extrapolation method and obtained the following equation by FE analysis;

$$Y_I=S/R\left(2.91+54.39 \alpha -391.4 {\alpha }^2+1210.6 {\alpha }^3-1650 {\alpha }^4+875.9 {\alpha }^5\right).$$

(4)

Being$\alpha =a/R$ indicates the ratio of the crack length to the radius. The low accuracy of the technique used to acquire Equation (4) in comparison to J-integral, encouraged Kuruppu et al. [13] to propose an alternative formulation. They used least-square for curve fitting where a second-order function of $\alpha$ is adopted for the simple relationship expressed in Equation (5). This equation has provided more accurate solutions in comparison to the ones obtained by Lim et al. [32] and Tutluoghlu and Keles [33].

$$Y_I=-1.297+9.516 \left(\frac{S}{R}\right)+\left(-0.47-16.457\left(\frac{S}{R}\right)\right)\alpha +\left(1.071+34.401\left(\frac{S}{R}\right)\right){\alpha }^2.$$

(5)

It must be noted that the recent formulation is valid if$\alpha \ge 0.2$. Nevertheless, a fairly deep notch is necessary on the bending effect to provide a robust mode I stress field near the notch tip. Therefore, it is suggested to respect$0.4\le \alpha \le 0.6$. Furthermore, it is recommended to consider $0.5\le S/R\le 0.8$ in which it is desirable to consider $S/R$ approaching to the maximum value [13].

Particularly in plane stress condition, the value of $K_{Ic}$ defined in Equation (3) at the point of crack extension is called the critical SIF value and expressed as $K_{Ic}$ in the literature [34]. It then defines the onset of crack extension. It does not necessarily indicate a fracture of the specimen; it depends on the crack stability. It is usually regarded as a material property and can be used to characterize toughness.

However, in this work, Equation (5) together with Equation (3) was considered as a reference theoretical solution on $K_{Ic}$ to validate the results obtained by numerical studies resorting to FEM study by the standard 2D constant strain triangle elements and meshless methods.

Considering $G$ as the strain energy release rate per crack tip, that is, for a double-ended crack within an infinite solid, the rate of release in strain energy per crack tip for a linearly-elastic material is:

$$G=\frac{K_I^2}{E}.$$

(6)

Thus, critical strain energy release rate$,G_c,$ is then related to the mode I fracture toughness calculated through the following equation:

$$G_c=\frac{K_{IC}^2}{E}.$$

(7)

Being $E$ the material elasticity modulus. It must be noted that Equation (7) will be used to obtain the theoretical solution and the numerical solutions on the critical fracture energy.

Regarding the radius of the fracture process zone, $r_c,$ so-called a critical distance from the crack tip, Taylor et al. [6] adopted an equation to determine the size of $r_c$ for brittle materials as follows:

$$r_c=\frac{1}{2\pi }{\left(\frac{K_{IC}^{ }}{{\sigma }_t}\right)}^2.$$

(8)

in which ${\sigma }_t$ denotes the material tensile strength. This term will be used to define the gap,$g,$ for fracture toughness calculations through numerical modeling.

Computationally, to determine the $K_{Ic}$within the obtained data, a numerical iterative routine was developed to analyze the stress field derived from numerical meshless and FEM analyses. It means that the stress field together with the point coordinates referenced to the crack tip must be extracted on a predefined window after the crack tip. Then, this data will be considered as an input for the developed numerical algorithm. This algorithm links the SIF computation routine and the stress calculation following principal stresses in the vicinity of the crack tip considering mode I condition.

Assuming the plane problem of a homogeneous isotropic solid, Williams’ series expansion on the plane stress would be appropriate as stated in the following formulations. Notice that the significance of $K_{\mathrm{I}}$ is demonstrated in Equation (10), [35]. Considering the maximum applied load, this mode I SIF will be the fracture toughness.

$${\sigma }_{yy}={\displaystyle \sum }_{n=1}^{\infty }\left(A_{\mathrm{I}n}\frac{n}{2}\right)r^{\frac{n}{2}-1}\left\{\left[2-{\left(-1\right)}^n-\frac{n}{2}\right]cos\left[\left(\frac{n}{2}-1\right)\theta \right]+\left(\frac{n}{2}-1\right)cos\left[\left(\frac{n}{2}-3\right)\theta \right]\right\},$$

(9)

$$K_{Ic}=A_{\mathrm{I}1}\sqrt{2\pi }.$$

(10)

In which $\left(r,\theta \right)$ represents the polar coordinates of the points referenced to the crack tip located in the defined window after the crack tip.

4. Analysis, Results and Discussion

A pre-defined cracked SCB specimen made of marble was considered in this study, as schematically depicted in Figure 1. Regarding geometrical characteristics, the specimen radius was considered as $R=25 \left(\mathrm{m}\mathrm{m}\right)$ and thickness was defined as$B=1 \left(\mathrm{m}\mathrm{m}\right)$. The mechanical parameters for the studied marble specimen are considered according to a reference article [36], reported in

Table 1. Mechanical properties of the studied SCB specimen, marble.

Parameter	Value
Young’s modulus	$E=75,500.0 \left(\mathrm{M}\mathrm{P}\mathrm{a}\right)$
Poisson’s ratio	$\nu =0.3$
Tensile strength	${\sigma }_t=16.1 \left(\mathrm{M}\mathrm{P}\mathrm{a}\right)$

.

Besides, the crack and the span lengths were accordingly defined as $a$ and $S$. In this study, different crack lengths and spans were taken into account led to obtaining the numerical results using meshless method formulations, see

Table 2. Specification on the span and crack lengths studied in this work.

$\mathbf{Span}\mathbf{Length},\mathit{S}\left(\mathbf{m}\mathbf{m}\right)$	$\mathbf{Crack}\mathbf{Length},\mathit{a}\left(\mathbf{m}\mathbf{m}\right)$	$\mathit{S}/\mathit{R}$	$\mathit{\alpha }=\mathit{a}/\mathit{R}$
$S_1=12.5$	$a_1=7.5$	$0.5$	${\alpha }_1=0.3$
$S_2=15.0$	$a_2=10.0$	$0.6$	${\alpha }_2=0.4$
$S_3=17.5$	$a_3=12.5$	$0.7$	${\alpha }_3=0.5$

.

As mentioned before, the critical mode I SIF$, K_{Ic},$ would be the most important parameter in fracture mechanics determining the resistance of the material to the impending urge to propagate the initial crack. This analysis is implemented on half of SCB specimen 2D domain under 3-point-bending condition due to symmetric geometry, as shown in Figure 2. The specimen contains an initial crack located in the middle of the SCB specimen for which the calculation of fracture attributes is taken into account.

Regarding the natural boundary condition, a compressive concentrated load was considered with a magnitude of $P=1.0 \left(\mathrm{k}\mathrm{N}\right)$ applied to the model vertically aligned with the initial crack. The essential boundary conditions were applied to the model constraining the nodes in the horizontal direction right above the crack tip to the end where the concentrated load is applied. On the other hand, depending on the various span lengths,$S$, a nodal restriction was imposed to restrict the model in the vertical direction to satisfy the requirements of the SCB specimen. Therefore, (as shown in Figure 2a), to simulate the boundary conditions regarding the span, an individual node ($P_A$) was singularly selected to restrict $v\left(x,y\right)$ displacement field in the corresponding point. Furthermore, a vertical line was assumed from the crack tip, ($P_B$), to the point, ($P_C$), under force to be bounded to $u\left(x,y\right)$.

The geometry of the model and the FE mesh was built through FEMAP© software. Then, the coordinates of the nodes/elements with their connectivity were imported as an input to the developed numerical algorithms implemented in MATLAB©. The numerical algorithms were based on a FEM code and used the FE mesh background to generate the nodal discretization. It was then extended to the LEFM theory to calculate the fracture toughness and strain energy release rate.

All stress and displacement types of boundary conditions were fully defined in the numerical models as explained. Regarding the FEM study, the problem domain was discretized by the standard 2D constant strain triangle elements, type S3. Then, the same discretization was used remaining nodal points to be analyzed by meshless methods.

For meshless method analyses, although in previous works [37,38], the number of nodes within an influence domain was suggested to be between $n=\left[9,16\right]$ nodes, Farahani et al. [26] recommended that each influence-domain should possess 20 nodes in order to get reliable results.

In this work, two different versions of NNRPIM were used, named NNRPIMv1 and NNRPIMv2 referring to the first-degree and the second-degree influence-cell, respectively. Referring to the previous works conducted on meshless methods by the authors [16,19,26], the parameters and the coefficients influencing the NNRPIM analysis were considered the same. This study follows the Radial Basis Function (RBF) formulation and definition as described in the previously published works by the authors, c.f. [19,25].

As Figure 2 depicts, the whole problem domain contains 4756 S3 elements and 2436 nodes.

In all numerical analyses, the discretization density of the problem domain, FE mesh as well as the nodal distribution of meshless methods, may affect the solution accuracy. Denser discretizations defined, more accurate results acquired. Hence, an interest region was defined in the problem domain possessing a denser discretization. This interest region is identified as the cracked area started from ($P_B$), right after the crack tip of the first crack length $a_1=7.5 \left(\mathrm{m}\mathrm{m}\right)$. It is dimensioned as $5\times 15\left(\mathrm{m}{\mathrm{m}}^2\right)$. This interest region was defined based on the first example, $a_1=7.5 \left(\mathrm{m}\mathrm{m}\right)$ and remains the same for the other crack length examples. It is demonstrated as a dash rectangular region highlighted in red as presented in Figure 2. It must be noted that this interest area includes a total number of 1037 nodes and 1920 elements of type S3.

Such interest region consideration permits to guarantee that the other crack tip will be located in a region with a denser FE mesh/meshless nodal distribution.

Concerning the FE mesh properties and meshless nodal distribution, it must be noted that the Authors have already studied convergence studies [26,39] on different densities for similar problems and it was concluded that this discretization owns sufficient potential to produce accurate and robust results.

In this study, referring to Table 2, three different crack lengths, $\left(a\right),$ are considered. Besides, three various span lengths $\left(S\right)$ are defined as another parametric value due to the considerable role of span length on the magnitude of bending moment around span coordinates leading to capture different results on the fracture characterization.

Owing to the abovementioned description, the numerical analyses have been performed and

Table 3. ${\mathit{K}}_{\mathit{I}\mathit{c} }$ obtained from theoretical formulation and numerical analyses.

$\mathit{S} \left(\mathbf{m}\mathbf{m}\right)$	$\mathit{a} \left(\mathbf{m}\mathbf{m}\right)$	${\mathit{K}}_{\mathit{I}\mathit{c} }\left(\mathit{M}\mathit{P}\mathit{a}\sqrt{\mathbf{m}\mathbf{m}}\right)$
		Theoretical	FEM	RPIM	NNRPIMv1	NNRPIMv2
12.5	7.5	242.30	242.03	243.32	236.83	237.13
	10.0	325.65	325.43	323.40	330.11	329.69
	12.5	461.17	461.89	460.92	462.45	463.03
15.0	7.5	316.81	315.10	318.74	316.40	315.89
	10.0	420.24	419.15	421.71	417.18	423.17
	12.5	585.10	584.57	586.38	583.31	585.31
17.5	7.5	391.32	391.12	392.73	383.03	386.58
	10.0	514.82	513.49	514.71	513.33	514.63
	12.5	709.02	710.71	710.23	711.92	707.85

reports the obtained results for the calculation of fracture toughness using all numerical simulations. It must be noted that the theoretical solution was calculated based on Equation (3) in which $Y_I$ was derived from Equation (5).

As explained in Section 3, the numerical fracture toughness was calculated in a defined window dimensioned as $3\times 5 \left(\mathrm{m}{\mathrm{m}}^2\right)$ considering a gap of$g\cong 0.5 \left(\mathrm{m}\mathrm{m}\right)$ ahead of the crack tip as illustrated in Figure 3.

Moreover,

Table 4. Deviation of ${\mathit{K}}_{\mathit{I}\mathit{c} }$ results obtained from numerical studies compared to the theoretical results.

$\mathit{S} \left(\mathbf{m}\mathbf{m}\right)$	$\mathit{a} \left(\mathbf{m}\mathbf{m}\right)$	$\mathbf{D}\mathbf{e}\mathbf{v}\mathbf{i}\mathbf{a}\mathbf{t}\mathbf{i}\mathbf{o}\mathbf{n} \left(\mathbf{\%}\right)$
		Theoretical and FEM 1	Theoretical and RPIM 2	Theoretical and NNRPIMv1 3	Theoretical and NNRPIMv2 4
12.5	7.5	0.11	0.42	2.26	2.14
	10.0	0.07	0.69	1.37	1.24
	12.5	0.16	0.05	0.28	0.40
15.0	7.5	0.26	0.61	0.13	0.29
	10.0	0.26	0.35	0.73	0.70
	12.5	0.09	0.22	0.31	0.04
17.5	7.5	0.05	0.36	2.12	1.21
	10.0	0.26	0.02	0.29	0.04
	12.5	0.24	0.17	0.41	0.17

¹$100\times \left|\left(K_{Ic}^{TH}-K_{ICc}^{FEM}\right)/K_{Ic}^{TH}\right|$; ²$100\times \left|\left(K_{Ic}^{TH}-K_{Ic}^{RPIM}\right)/K_{Ic}^{TH}\right|$; ³$100\times \left|\left(K_{Ic}^{TH}-K_{Ic}^{NNRPIMv1}\right)/K_{Ic}^{TH}\right|$; ⁴$100\times \left|\left(K_{Ic}^{TH}-K_{Ic}^{NNRPIMv2}\right)/K_{Ic}^{TH}\right|$.

shows the deviation computed amongst the numerical results compared to the theoretical solution. For all numerical calculations on fracture toughness, it must be stated that to avoid the plastic effect that occurred at the crack tip, a gap of $g\cong 0.5 \left(\mathrm{m}\mathrm{m}\right)$ was considered after the crack tip. This gap was defined following the plastic radius, $r_c,$ computation derived from Equation (8) as calculated by Xie et al. [39] for the similar SCB specimen.

Besides,

Table 5. Critical strain energy release rate obtained for different methods.

$\mathit{S} \left(\mathbf{m}\mathbf{m}\right)$	$\mathit{a} \left(\mathbf{m}\mathbf{m}\right)$	${\mathit{G}}_{\mathit{C} }\left(\mathbf{M}\mathbf{P}\mathbf{a}.\mathbf{m}\mathbf{m}\right)$
		Theoretical	FEM	RPIM	NNRPIMv1	NNRPIMv2
12.5	7.5	0.78	0.78	0.78	0.74	0.74
	10.0	1.40	1.420	1.39	1.44	1.44
	12.5	2.82	2.83	2.81	2.83	2.84
15.0	7.5	1.33	1.32	1.35	1.33	1.32
	10.0	2.34	2.33	2.36	2.31	2.37
	12.5	4.53	4.53	4.55	4.51	4.54
17.5	7.5	2.03	2.03	2.04	1.94	1.98
	10.0	3.51	3.49	3.51	3.49	3.51
	12.5	6.66	6.69	6.68	6.71	6.64

presents the critical strain energy calculation derived from the numerical analyses calculated through Equation (7). In conclusion, $G_C$ raised if the crack and span lengths increased. Straightforwardly, the higher fracture toughness and $G_C,$ then the greater the material’s fracture resistance, which then possibly translates into higher damage tolerance.

Figure 4 depicts the fracture toughness variation in terms of various $\alpha$ for different normalized span $\left(S/R\right)$ obtained from theoretical and numerical solutions. In these graphs, a reasonable increment in fracture toughness is noticeable as the crack length grows. All numerical methods provided similar results to the ones obtained from the theoretical formulation. The acquired curves experienced a parabolic behavior as anticipated by the former studies [32,33].

The stress distribution in y-direction, ${\sigma }_{yy},$ was also obtained from all numerical studies for different crack lengths and various normalized spans$\left(S/R\right)$, as demonstrated in Figure 5, Figure 6 and Figure 7. It must be noted that the stress contours were extracted in a window dimensioned as $3\times 5 \left(\mathrm{m}{\mathrm{m}}^2\right)$ considering a gap of$g\cong 0.5 \left(\mathrm{m}\mathrm{m}\right)$ ahead of the crack tip as illustrated in Figure 3.

The smooth contours on the integration nodes of meshless methods and the FEM infer that the density and the area of the stress distribution enhance by increasing the crack length and span. The cracked area accounted for the maximum stress conforming to the LEFM theory. As the crack length raised, the stress singularities increased around the crack tip. In most cases, NNRPIMv1 and NNRPIMv2 depict a much similar stress distribution in comparison to the RPIM and FEM even though there is no significant difference amongst the profiles.

The verified results acquired from different methods in terms of accuracy would prove the potential capability of meshless methods to provide accurate solutions. Therefore, it can be inferred that the obtained results have shown promising accuracy leading to validate the numerical methodologies.

5. Conclusions

This study deals with the fracture characterization on a marble SCB specimen using advanced discretization techniques. The numerical model was prepared through the FEM in FEMAP© software and the essential data of FE mesh was exported to feed the numerical algorithms developed in MATLAB©.

The numerical elastostatic algorithm was implemented based on the FEM fundamentals and RBF formulation in which the meshless nodal discretization was generated based on the FE mesh background, which was built by the standard 2D constant strain triangle elements. Then, respecting the FE mesh property, maintaining the nodal points, the RPIM and NNRPIM meshless models were generated in which the latter one employed on both first- and second-influence cell versions. As an output, the fracture toughness$, K_{Ic},$ and the strain energy release rate$, G_c,$ were computed.

The obtained numerical results (from FEM and meshless models) on $K_{Ic}$and $G_c$ were compared to a reference theoretical solution proposed for the SCB specimen available in the literature. It aimed at describing how the fracture toughness of the studied specimen can be affected by the crack and span lengths defined in the geometry. However, the results infer that $K_{Ic}$, $G_c$ and stress gradient around the crack tip accounted for a greater value if the crack length grows or the span length increases. The acquired results imply that although the FEM solutions were closer to the reference solution, meshless methods results were sufficiently accurate and promising.

A general discussion on the meshless method privileges can be drawn as:

The meshless methods are capable to offer numerically advantageous models to evaluate the fracture characterization of the engineering structures since their formulations are based on the FE mesh background (nodes) and some difficulties that FEM models may face are already overcome. These issues can be identified as the dependency on the element type and shape, FE distortions, convergence study on the nonlinearity, high degrees of freedom and FE remeshing in terms of crack propagation problems. Besides, in terms of computational cost, it must be mentioned that the meshless methods would have the potential to reduce the computational costs since they are solving the problems only using the nodal distribution in the problem domain. Moreover, it will be feasible to allocate denser discretization in the interest region/influence domain of the corresponding model. This aspect is more significant in solving large structures with complex material behaviors, for instance, elastoplastic analyses, damage models and porous materials.

However, a simple model was studied in this work but it can be pointed out that the meshless methods formulations extended to the LEFM are practically capable to simulate the fracture problems providing sufficiently accurate results to validate the supporting methodology.