Utilizing cylindrical and cubical specimens with edge notch to determine size-independent fracture quantities of rock materials

In materials science, clay products such as bricks, natural stones, and concrete are generally classified as ceramic materials. From ancient times to the present day, bricks and natural stones such as limestone, granite, and sandstone, as well as Roman concrete, have been widely used in buildings, bridges, and retaining walls. These materials are categorized as quite brittle in materials science, meaning that while their compressive strength is reasonable, their tensile strength is poor. Therefore, in practice, ceramic materials are simulated using the Mohr-Coulomb criterion, one of the classical fracture theories, where materials are assumed to be flawless. Moreover, bricks, natural stones, and concrete contain many voids, both open and closed.
In the 1920s, experimental research on brittle materials such as glass revealed that no material was perfect. Subsequently, the principle of linear elastic fracture mechanics (LEFM) was introduced. Irwin [1] improved this theory in the 1950s and proposed two equivalent material quantities, fracture toughness, and the critical strain energy release rate, to simulate structures containing cracks. Irwin's approach was applied to concrete beams of different sizes in 1961 and observed that Irwin's material quantities were size-dependent for concrete. After the 1970s, experimental studies on the aforementioned ceramic materials revealed that these materials exhibit quasi-brittle failure, as a region containing a high density of microcracks, known as the fracture process zone (FPZ), can develop in the vicinity of the crack tip. In quasi-brittle materials, the initially sharp crack tip becomes blunted with increasing loading. As a result, failure in a cracked structure occurs only after the crack reaches a certain depth proportional to the FPZ, differing from truly brittle materials. Since 1976, various nonlinear fracture models, which employ at least two fracture quantities, have been developed to account for crack-tip blunting in quasi-brittle materials such as cementitious composites, rock, and asphalt [2-4].
The first applications of these fracture models focused on concrete and mortar materials. Until 1990, test results on rock material beams and compact tension (CT) specimens of various sizes revealed that the values of fracture toughness and the critical J-integral were independent of specimen width but dependent on specimen size. Considering this, a compliance technique commonly used in concrete fracture was first proposed for chevron bend specimens and short rod specimens to determine the nonlinear fracture toughness of rock materials by ISRM (International Society for Rock Mechanics) [5]. Bažant and Kazemi [3] successfully applied the size effect model (SEM) to experimental data obtained from beams and CT specimens of Hashida and Takahashi's Iidate granite. Bažant et al. [6] simulated the results of fracture tests on Indiana limestone beams of four different sizes using the SEM. Subsequently, Ouyang et al. [7] modeled Bažant and co-workers' [6] test results on Iidate granite using the two-parameter model (TPM).
Although beams and CT specimens containing cracks are widely used to determine fracture quantities of quasi-brittle materials and metals, the use of compact specimens, illustrated in Fig. 1, has recently become more common for cementitious materials, rock, and asphalt concrete. The wedge-splitting (WS) specimens and compact compression specimens, shown in Figs. 1a and 1b, were developed as alternatives to CT specimens and have been widely employed to determine the nonlinear fracture quantities of cementitious composites such as mortar and concrete. While splitting specimens in cylindrical and cubical forms have been widely used to indirectly evaluate the tensile strength of quasi-brittle materials, their centrally notched forms, shown in Fig. 1c, have also been applied to determine the nonlinear fracture quantities of cementitious materials. Ince [8] modeled edge-notched mortar cubes, which behave as splitting specimens,
using the TPM (Fig. 1d). The semi-circular bending (SCB) specimens illustrated in Fig. 1e are the most commonly used samples in fracture testing of rock and asphalt materials since they can easily be produced from cylindrical core samples. Many fracture tests on SCB specimens have been conducted to evaluate the critical J-integral value and the LEFM-based fracture toughness of asphalt composites. Zegeye et al. [9] modeled asphalt SCB specimens of various sizes at low temperatures using the SEM. Ince et al. [10] initially derived some LEFM relationships regarding crack blunting in SCB specimens and subsequently simulated various SCB tests on asphalt from the literature using the TPM in concrete fracture. While many studies on rock SCB specimens have been performed based on LEFM by several researchers, other investigators have modeled rock materials by considering the nonlinear fracture behavior of these materials. Wei et al. [11] modeled SCB specimens containing a notch and chevron-notched SCB specimens using the finite element method. They then conducted fracture tests on Changtai granite specimens of three different depths and determined the nonlinear fracture toughness parameters using their equivalent elastic fracture model. Guo et al. [12] determined fracture quantities based on the SEM for granite SCB specimens of various sizes at room and high temperatures. Recently, Ince [13] conducted two series of experimental studies with SCB and beam specimens made of concrete and mortar and discussed them based on the three most popular fracture models: the TPM, the SEM, and the double-K model. Subsequently, three sets of SCB rock tests from the literature were analyzed using the aforementioned concrete fracture models, emphasizing that the nonlinear fracture quantities of both rock and concrete materials can be effectively assessed using single-sized SCB specimens.
Tutluoglu and Keles [14,15] derived the stress intensity factor (SIF) values for straight-notched disk bending (SNDB) specimens with different height-to-diameter ratios (Fig. 1f), a simplified form of the short rod specimen proposed by ISRM [5]. Subsequently, they conducted bending fracture tests on SNDB specimens, SCB specimens, and cracked chevronnotched Brazilian disks composed of andesite and marble to comparatively determine mode I fracture toughness. They proposed a numerical model to compute the fracture process zone (FPZ) depth of these stone materials and concluded that the FPZ depth of the SNDB specimen is smaller than that of other specimens. As a result, deviations from the assumptions of linear elastic fracture mechanics (LEFM) are reduced. In another study, Tutluoglu and co-workers [16] performed extensive experiments on andesite using SNDB specimens, SCB specimens, and three-point bending plates. Their results revealed that the critical normal strain value could be considered a material property across all tested specimen types. The aforementioned studies primarily focused on compact specimens subjected exclusively to pure mode I loading. Ayatollahi and Aliha [17] modeled both SCB specimens and splitting cylinders under mixed-mode loading using the finite element method.
In this study, SNDB specimens and single-notch cube bending (SNCB) specimens, cubic forms of SNDB specimens that differ from the edge-notched cube shown in Fig. 1d only in boundary conditions, were investigated both experimentally and numerically. First, a two-dimensional crack analysis of SNCB specimens was conducted using the finite element method (FEM) in ANSYS, determining SIF values, crack blunting expressions, and Green's functions. A three-dimensional crack analysis of SNCB specimens was also performed to obtain SIF values. Additionally, the existing SIF formulation for SNDB specimens with a diameter-to-height ratio of 2, originally developed by Tutluoglu and Keles [14], was refined based on this 3D analysis. Other LEFM equations, including crack blunting expressions and Green's functions, were also derived for SNDB specimens. Using these LEFM formulations, test data from Tutluoglu and Keles [14] were simulated using both the compliance method and the peak load method based on the two-parameter model (TPM), and the results were compared with the numerical model proposed by Tutluoglu and Keles [14]. A significant statistical correlation was observed between the three methods. Furthermore, fracture tests were conducted on Elazig limestone cubes with different edge crack depths but identical specimen sizes. The experimental results were analyzed using the peak load method with the SNCB specimen expressions. Consequently, this study emphasizes that, unlike the compliance approach, the nonlinear fracture properties of rock materials can be efficiently determined using the peak loads of single-sized cylindrical and cubic specimens with different edge-notch depths, eliminating the need for complex testing equipment.

Building elements consisting of quasi-brittle materials, such as rock and cement-based materials, exhibit a three-stage behavior until failure occurs. As seen in a beam with initial crack length $a_0$$a_0$a_(0)a_{0}$a_0$, whose load-crack mouth opening displacement relationship is presented in Fig. 2, these stages include: (1) crack initiation, corresponding to the range between the initial state and the first cracking load $\left(P_i\right)$$\left(P_i\right)$(P_(i))\left(P_{i}\right)$\left(P_i\right)$; (2) stable crack growth, corresponding to the range between the first cracking load and the peak load $\left(P_c\right)$$\left(P_c\right)$(P_(c))\left(P_{c}\right)$\left(P_c\right)$; and (3) unstable crack growth after the peak load. The beam, which behaves almost linearly elastic due to defects in the material until the first cracking load, exhibits non-linear behavior after this point due to the formation of a fracture process zone (FPZ) behind the free crack (Fig. 2c).

To fully characterize this multiphase behavior, cohesive crack models, which take into account the cohesive nature of the FPZ, simulate quasi-brittle structural members using the finite element method and the boundary element method. On the other hand, equivalent elastic fracture models, which consider only the behavior up to the peak load of the structural member, aim to predict the crack extension at the peak load ( $\mathrm{\Delta }{a}_c$$\mathrm{\Delta }{a}_c$Deltaa_(c)\Delta a_{c}$\mathrm{\Delta }{a}_c$ ), which is proportional to the FPZ depth, as illustrated in

Fig. 2c [2-4]. It is emphasized that none of the above-mentioned fracture models employ a single fracture parameter, such as fracture toughness or fracture energy, contrary to classical LEFM.
Among the equivalent elastic fracture models, the two-parameter model (TPM) [2] predicts only the unstable crack growth at peak load in quasi-brittle structures by utilizing two critical quantities corresponding to the critical crack length $(a_c=a_0+\mathrm{\Delta }{a}_c)$$\left(a_c=a_0+\mathrm{\Delta }{a}_c\right)$(a_(c)=a_(0)+Deltaa_(c))\left(a_{c}=a_{0}+\Delta a_{c}\right)$(a_c=a_0+\mathrm{\Delta }{a}_c)$ : the fracture toughness $\left(K_{Ic}^s\right)$$\left(K_{Ic}^s\right)$(K_(Ic)^(s))\left(K_{I c}^{s}\right)$\left(K_{Ic}^s\right)$ and the crack tip opening displacement value $\left(CTO{D}_c\right)$$\left(CTO{D}_c\right)$(CTOD_(c))\left(C T O D_{c}\right)$\left(CTO{D}_c\right)$, which describes the blunting at the crack tip well, as defined above:

Here, ${\sigma }_{Nc}$${\sigma }_{Nc}$sigma_(Nc)\sigma_{N c}${\sigma }_{Nc}$ refers to the strength value computed according to the classical strength of materials theory for unnotched specimens, and it is commonly referred to as nominal strength. $E$$E$EE$E$ in Eqn. (2) represents Young's modulus of the material. In Eqns. (1) and (2), the normalized functions $Y,V_1$$Y,V_1$Y,V_(1)Y, V_{1}$Y,V_1$, and $M$$M$MM$M$ are based on crack lengths ( $a_0$$a_0$a_(0)a_{0}$a_0$ and $a_c$$a_c$a_(c)a_{c}$a_c$ ) and the specimen depth (d). These functions can be computed from crack analysis by employing numerical methods such as the finite element method (FEM) for the specimen considered or may be found in LEFM books for simple specimen types. As an example, the normalized functions employed for beam specimens with span/depth=2.5, which are also used in this study, can be given below:

The double- $K$$K$KK$K$ fracture model, based on equivalent elastic fracture models, simulates failure in quasi-brittle structures by utilizing two quantities: initiation and unstable fracture toughness ( $K_{Ic}^{un}$$K_{Ic}^{un}$K_(Ic)^(un)K_{I c}^{u n}$K_{Ic}^{un}$ and $K_{Ic}^{ini}$$K_{Ic}^{ini}$K_(Ic)^(ini)K_{I c}^{i n i}$K_{Ic}^{ini}$ ). Therefore, contrary to the TPM, both crack initiation and unstable crack growth in a quasi-brittle structure can be determined using Eqn. (1) through this model, where the values of $a_0$$a_0$a_(0)a_{0}$a_0$ and the first cracking load $\left(P_i\right)$$\left(P_i\right)$(P_(i))\left(P_{i}\right)$\left(P_i\right)$ are utilized to determine $K_{Ic}^{un}$$K_{Ic}^{un}$K_(Ic)^(un)K_{I c}^{u n}$K_{Ic}^{un}$.
As pointed out above, the only unknown points in Eqns. (1) and (2) are the critical (or effective) crack depth ( $a_c$$a_c$a_(c)a_{c}$a_c$ in the procedure for determining fracture quantities in the aforementioned fracture approaches. According to concrete fracture codes recommended by RILEM (International Union of Laboratories and Experts in Construction Materials, Systems and Structures), as illustrated in Fig. 2b, this critical quantity can be evaluated from the load-crack mouth opening displacement (P-CMOD) response of the test specimen (preferably a beam). However, to determine the P-CMOD response of the specimen, the use of a feedback-testing machine with CMOD control is necessary for the TPM, whereas no such requirement exists for the double- $K$$K$KK$K$ model.
In approaches based on RILEM, Young's modulus ( $E$$E$EE$E$ ) of the tested specimen is initially determined from the initial compliance value $(C_i=CMO{D}_i/P_i)$$\left(C_i=CMO{D}_i/P_i\right)$(C_(i)=CMOD_(i)//P_(i))\left(C_{i}=C M O D_{i} / P_{i}\right)$(C_i=CMO{D}_i/P_i)$, where $P_i$$P_i$P_(i)P_{i}$P_i$ is commonly assumed to be $P_c/2$$P_c/2$P_(c)//2P_{c} / 2$P_c/2$ due to difficulties in determining $P_i$$P_i$P_(i)P_{i}$P_i$ (Fig. 2b). Subsequently, assuming that $E$$E$EE$E$ remains constant for any pair of $P$$P$PP$P$ and $CMOD$$CMOD$CMODC M O D$CMOD$, the $a_c$$a_c$a_(c)a_{c}$a_c$ value of the notched specimen is computed using the secant compliance value ( $C_s$$C_s$C_(s)C_{s}$C_s$ ) in the double- $K$$K$KK$K$ model or the unloading compliance value ( $C_u$$C_u$C_(u)C_{u}$C_u$ ) in the TPM, as indicated in Fig. 2b. Although both models employ different test equipment and compliance techniques, experimental
studies on beams and wedge-splitting specimens by Xu et al. revealed that the fracture toughness and $CTO{D}_c$$CTO{D}_c$CTOD_(c)C T O D_{c}$CTO{D}_c$ values obtained using the TPM were in good statistical agreement with those of the double- $K$$K$KK$K$ model.
An alternative method, which utilizes only the ultimate loads of at least three specimens with different notch depths, was developed for the TPM by Tang et al. [18]. In this so-called peak-load method in concrete fracture, both fracture parameters ( $K_{Ic}^s$$K_{Ic}^s$K_(Ic)^(s)K_{I c}^{s}$K_{Ic}^s$ and $CTO{D}_c$$CTO{D}_c$CTOD_(c)C T O D_{c}$CTO{D}_c$ ) are treated as unknowns instead of $a_c$$a_c$a_(c)a_{c}$a_c$ values for the specimens. Consequently, this system, which consists of multiple equations but only two unknowns, can only be solved using an optimization-based approach. To achieve this, Tang et al. proposed that $K_{Ic}^s$$K_{Ic}^s$K_(Ic)^(s)K_{I c}^{s}$K_{Ic}^s$ and $CTO{D}_c$$CTO{D}_c$CTOD_(c)C T O D_{c}$CTO{D}_c$ curves corresponding to specific notch depths should first be determined for each specimen with a distinct $a_0$$a_0$a_(0)a_{0}$a_0$. The values of the sample standard deviation ( $s$$s$ss$s$ ) corresponding to the same $K_{Ic}^s$$K_{Ic}^s$K_(Ic)^(s)K_{I c}^{s}$K_{Ic}^s$ value are then calculated using Eqn. (6):

Here, $\overline{CTO{D}_c}$$\overline{CTO{D}_c}$bar(CTOD_(c))\overline{C T O D_{c}}$\overline{CTO{D}_c}$ is the mean value of $CTO{D}_{ci}$$CTO{D}_{ci}$CTOD_(ci)C T O D_{c i}$CTO{D}_{ci}$ values. Consequently, the $K_{Ic}^s$$K_{Ic}^s$K_(Ic)^(s)K_{I c}^{s}$K_{Ic}^s$ and $CTO{D}_c$$CTO{D}_c$CTOD_(c)C T O D_{c}$CTO{D}_c$ values, where $s\left(K_{Ic}^s\right)$$s\left(K_{Ic}^s\right)$s(K_(Ic)^(s))s\left(K_{I c}^{s}\right)$s\left(K_{Ic}^s\right)$ is at its minimum, are accepted as the fracture quantities according to the peak-load method.

As stated above, determining the initial cracking load $\left(P_i\right)$$\left(P_i\right)$(P_(i))\left(P_{i}\right)$\left(P_i\right)$ of the specimens is difficult. Therefore, Xu and Reinhardt [4] recommended a simplified approach, known as the inverse method in concrete fracture, as follows:

Here, $K_{Ic}^c$$K_{Ic}^c$K_(Ic)^(c)K_{I c}^{c}$K_{Ic}^c$ refers to the cohesive fracture toughness. As depicted in Fig. 3a, $K_{Ic}^c$$K_{Ic}^c$K_(Ic)^(c)K_{I c}^{c}$K_{Ic}^c$ corresponds to the toughness value consisting of cohesive stresses assumed to occur along the critical crack extension line ( $\mathrm{\Delta }{a}_c$$\mathrm{\Delta }{a}_c$Deltaa_(c)\Delta a_{c}$\mathrm{\Delta }{a}_c$ ). According to the inverse method, this cohesive stress distribution is trapezoidal, such that the maximum value of this distribution is supposed to correspond to the tensile strength of the material ( $f_t$$f_t$f_(t)f_{t}$f_t$ ) while its lower value, ${\sigma }_s\left(CTO{D}_c\right)$${\sigma }_s\left(CTO{D}_c\right)$sigma_(s)(CTOD_(c))\sigma_{s}\left(C T O D_{c}\right)${\sigma }_s\left(CTO{D}_c\right)$, is the stress value determined according to the chosen stress-crack tip opening displacement function ( $\sigma -w$$\sigma -w$sigma-w\sigma-w$\sigma -w$ ), as indicated in Fig. 3a. In many applications up to the peak load in concrete fracture, using a linear function, as shown in Fig. 3b, provides an adequate approach for $\sigma$$\sigma$sigma\sigma$\sigma$ $w$$w$ww$w$ [19]. The area under this function is equal to the fracture energy ( $G_f$$G_f$G_(f)G_{f}$G_f$ ), which is converted to fracture toughness using Irwin's G-K relation: $G_f={\left(K_{Ic}^{un}\right)}^2/E[1,3]$$G_f={\left(K_{Ic}^{un}\right)}^2/E[1,3]$G_(f)=(K_(Ic)^(un))^(2)//E[1,3]G_{f}=\left(K_{I c}^{u n}\right)^{2} / E[1,3]$G_f={\left(K_{Ic}^{un}\right)}^2/E[1,3]$. On the other hand, the toughness value created by such a trapezoidal stress distribution at the crack tip can be calculated by integrating the toughness values created by the individual forces ( $P=\sigma dx$$P=\sigma dx$P=sigma dxP=\sigma d x$P=\sigma dx$ ) shown in Fig. 3c:

Here, $x$$x$xx$x$ refers to the vertical distance from the crack mouth of the specimen, and $G$$G$GG$G$ is Green's function. For instance, Tada et al. [20] derived the following Green's function for an infinite strip as shown in Fig. 3c:

The results of the experiments on compact quasi-brittle samples by Ince [21] indicated that the fracture quantities ( $K_{Ic}^{un}$$K_{Ic}^{un}$K_(Ic)^(un)K_{I c}^{u n}$K_{Ic}^{un}$ and $K_{Ic}^{ini})$$\left.K_{Ic}^{ini}\right)${:K_(Ic)^(ini))\left.K_{I c}^{i n i}\right)$K_{Ic}^{ini})$ of the double-K model can also be evaluated using the peak-load method.

Two prominent energy methods based on the FEM, such as the J-integral and the crack closure integral (CCI) technique[22], which can also be successfully utilized for three-dimensional crack problems, were employed in this study. Although conventional methods commonly used in FEM crack problems, such as extrapolation methods based on stress matching, displacement matching, and the hybrid method, require a large number of fine elements near a crack for accurate results, the aforementioned methods do not have such refined mesh requirements. The J-integral is a contour integral evaluated along a path that encircles the crack tip, beginning and concluding at the surfaces of the crack. Since this integral is path-independent for any of the contours, it has the same value for any similar contour. The use of the CCI technique is not as widespread as the J-integral approach and is not embedded in many commercial FEM programs, but the implementation of the CCI technique is very simple compared to the J-integral method, as will be detailed below. The basic principles of the crack closure integral (CCI) technique can be derived using the fundamental principle of LEFM [1]. To determine the relationship between the energy release rate $G$$G$GG$G$ and the fracture toughness $K$$K$KK$K$, it is necessary to account for the work of the closing forces ( $P=\sigma dx$$P=\sigma dx$P=sigma dxP=\sigma d x$P=\sigma dx$ ), which corresponds to the initial state of the new crack profile when a crack of length $a$$a$aa$a$ increases by a small amount ( $\mathrm{\Delta }a$$\mathrm{\Delta }a$Delta a\Delta a$\mathrm{\Delta }a$ ), as shown in Fig. 4a. For a plate with a thickness of $b$$b$bb$b$, the energy released as a result of the small extension of the crack ( $Gb\mathrm{\Delta }a$$Gb\mathrm{\Delta }a$Gb Delta aG b \Delta a$Gb\mathrm{\Delta }a$ ) must be equal to the work done by the internal forces that closed it. Accordingly, the following equation can be written:

Here, $v$$v$vv$v$ is the crack opening displacement (COD). The above work expression can easily be applied by running two-step computer numerical approaches, such as the finite element method, where the elastic body model is assumed to be continuous only at the selected nodes. To calculate the crack closure work in the body modeled with the linear plane finite elements, as shown in Fig. 4b, the nodal force at the crack tip ( $P_i$$P_i$P_(i)P_{i}$P_i$ ) is determined for the original crack profile in the first step. Subsequently, the node at the crack tip is released, and the nodal displacement $\left(v_i\right)$$\left(v_i\right)$(v_(i))\left(v_{i}\right)$\left(v_i\right)$ is determined for the extended crack profile in the second step. On the other hand, when sufficiently small linear elements are used, the simple work equation, which can also be computed by utilizing the displacement value at the node just before the crack tip $\left(v_{i-1}\right)$$\left(v_{i-1}\right)$(v_(i-1))\left(v_{i-1}\right)$\left(v_{i-1}\right)$, could be obtained within acceptable limits by running a single step as follows $(\mathrm{W}(a)\approx \mathrm{W}(a+\mathrm{\Delta }a))$$(\mathrm{W}(a)\approx \mathrm{W}(a+\mathrm{\Delta }a))$(W(a)~~W(a+Delta a))(\mathrm{W}(a) \approx \mathrm{W}(a+\Delta a))$(\mathrm{W}(a)\approx \mathrm{W}(a+\mathrm{\Delta }a))$ :

Consequently, $G$$G$GG$G$ and $K$$K$KK$K$ can be determined for a cracked body modeled with linear plane finite elements by utilizing Eqns. (10) and (11) as follows:

Similarly, as depicted in Fig. 4c and Fig. 4d, in elastic crack problems where 8-node plane elements with or without singular are used, $G$$G$GG$G$ (and hence $K$$K$KK$K$ ) can be calculated as follows, respectively.

In this study, SNCB and SNDB specimens were simulated as deep beams subjected to three-point bending. While Tutluoglu and Keles [14] modeled SNDB specimens, as shown in Fig. 1f, they applied the external load to the specimen through rigid shell elements with an arch shape. However, Tang [23] revealed that the width of the distributed load has an insignificant effect on structures under three-point bending. Recently, Ince et al. [10] simulated semi-circular bending (SCB) specimens for asphalt concrete and concluded that this effect may also be ignored for SCB specimens. For instance, the relative error between concentrated loading and distributed loading was smaller than $0.7\mathrm{\%}$$0.7\mathrm{\%}$0.7%0.7 \%$0.7\mathrm{\%}$ in the computations of the SIFs of SCB specimens.
In this study, the SNCB specimen illustrated in Fig. 5a was initially analyzed using 8-noded plane elements. Because of the symmetry of the specimen, half of it, with a depth of $d=100\mathrm{mm}$$d=100\mathrm{mm}$d=100mmd=100 \mathrm{~mm}$d=100\mathrm{mm}$ and a width of $b=1\mathrm{mm}$$b=1\mathrm{mm}$b=1mmb=1 \mathrm{~mm}$b=1\mathrm{mm}$, was considered in the modeling, as shown in Fig. 5b. The specimen was discretized into 100 finite elements to minimize computation errors along the notch line. Furthermore, a total of eight quarter-point elements with six nodes were employed to consider the stress singularity around the notch. Note that Ince [21] previously utilized the same discretization along the notch line to model split-tension
cylinders containing central cracks, for which an analytical solution was available. Consequently, it was observed that this discretization provided an accuracy of $0.2\mathrm{\%}$$0.2\mathrm{\%}$0.2%0.2 \%$0.2\mathrm{\%}$.

In the finite element analysis (FEA) by using ANSYS Parametric Design Language (APDL), the elastic constants of the material, elasticity modulus ( $E$$E$EE$E$ ) and Poisson's ratio ( $v$$v$vv$v$ ), were selected as 1.0 MPa and 0.2 , respectively, and the concentrated load $(P)$$(P)$(P)(P)$(P)$ at the mid-span was set to 100 N . The relative notch depths ( ${\alpha }_0=a_0/d$${\alpha }_0=a_0/d$alpha_(0)=a_(0)//d\alpha_{0}=a_{0} / d${\alpha }_0=a_0/d$ ) were varied between 0.1 and 0.9 with a step interval of 0.1 , while span-to-depth ( $s/d$$s/d$s//ds / d$s/d$ ) ratios were chosen as $0.5,0.6,0.7$$0.5,0.6,0.7$0.5,0.6,0.70.5,0.6,0.7$0.5,0.6,0.7$, and 0.8 . The normalized function ( $Y$$Y$YY$Y$ ), based on the SIF of the SNCB specimen, was computed using both the J-integral and the CCI technique based on Eqn. (13) for each ${\alpha }_O$${\alpha }_O$alpha_(O)\alpha_{O}${\alpha }_O$ value. In this procedure, the value of the nominal strength ( ${\sigma }_N$${\sigma }_N$sigma_(N)\sigma_{N}${\sigma }_N$ ) in Eqn. (1) was taken as $P/(2bd)$$P/(2bd)$P//(2bd)P /(2 b d)$P/(2bd)$, which was previously used for the SCB specimens and SNDB specimens by Tutluoglu and Keles [14]. In above 2D FEA, the type of PLANE183 elements in ANSYS element library ( 8 -noded isoparametric element) were utilized except around the crack, where 8 quarterpoint elements (6-noded singular elements) were used to account for the stress singularity at the crack tip.
To simulate the SNDB specimen reliably, the SNCB specimen was also modeled in 3D. For this, a quarter of the SNCB specimen was considered due to the symmetry of the cubic geometry ( $b=100\mathrm{mm}$$b=100\mathrm{mm}$b=100mmb=100 \mathrm{~mm}$b=100\mathrm{mm}$ ), as indicated in Fig. 5c. The 3D SNCB specimen containing a notch was modeled using both 8 -noded (SOLID185) and 20-noded (SOLID186) solid elements without singular finite elements. In FEA, the elastic constants of the material were chosen as in the 2D analysis, while the concentrated load $(P=100\mathrm{N})$$(P=100\mathrm{N})$(P=100N)(P=100 \mathrm{~N})$(P=100\mathrm{N})$ at the mid-span was distributed along the specimen width, as depicted in Fig. 5c. The Y values of the 3D SNCB specimens were calculated using the J-integral and the CCI technique for each ${\alpha }_0$${\alpha }_0$alpha_(0)\alpha_{0}${\alpha }_0$ and $s/d$$s/d$s//ds / d$s/d$ ratio selected above. These values were determined for each finite element at the crack tip in the model with 20 nodes, whereas they were computed at the side surface of each finite element at the crack tip in the model with 8 nodes. In the crack analysis based on the CCI approach, Eqns. (12) and (14) were employed for the models with 8 -noded and 20 -noded solid elements, respectively.
In Fig. 6a, the $Y$$Y$YY$Y$ values computed for ${\alpha }_0=0.5$${\alpha }_0=0.5$alpha_(0)=0.5\alpha_{0}=0.5${\alpha }_0=0.5$ are illustrated for $s/d=0.5,0.6$$s/d=0.5,0.6$s//d=0.5,0.6s / d=0.5,0.6$s/d=0.5,0.6$, and 0.8 according to each solid element type. It is observed from Fig. 6a that there is a very close coincidence between both calculation methods, except for finite elements
at the boundary of the specimen. The mean values, excluding the boundary values of the CCI technique, were also plotted according to both crack analyses in Fig. 6a. Based on this, it can be concluded that the solutions with 8-node elements are sufficiently compatible with those of 20 -node elements.

Fig. 6b shows the FEA results of both 2D SNCB specimens with 8 nodes and 3D SNCB specimens with 20 nodes for all ${\alpha }_0$${\alpha }_0$alpha_(0)\alpha_{0}${\alpha }_0$ and $s/d$$s/d$s//ds / d$s/d$ values considered, comparatively. The equations of the regression line ( $Y=AX+C$$Y=AX+C$Y=AX+CY=A X+C$Y=AX+C$ ), the determination coefficients ( $R^2$$R^2$R^(2)R^{2}$R^2$ ) corresponding to these equations, and the maximum relative error (accuracy) computed for each data point were also reported in this figure. In an ideal model, the parameters $A,C$$A,C$A,CA, C$A,C$, and ${\mathrm{R}}^2$${\mathrm{R}}^2$R^(2)\mathrm{R}^{2}${\mathrm{R}}^2$ should attain values of 1,0 , and 1 , respectively. It can be seen from this figure that the results of the SNCB models based on 2D and 3D were satisfactorily consistent according to the J-integral way and the CCI technique.

Based on the information above, the SNDB specimens with only a depth-to-diameter ratio $(d/D)$$(d/D)$(d//D)(d / D)$(d/D)$ of 0.5 were simulated under the same conditions as the SNCB specimens for four different $s/d$$s/d$s//ds / d$s/d$ ratios: $0.25,0.30,0.35$$0.25,0.30,0.35$0.25,0.30,0.350.25,0.30,0.35$0.25,0.30,0.35$, and 0.4 . The FEM discretization of the quarter part of the SNDB specimen is illustrated in Fig. 5d.
In this study, the results of 2D FEA were used while deriving the following LEFM formulas for the SNCB specimens. According to the FEM discretization in Fig. 5b, in addition to the SIF, the vertical displacement values at each node along the crack were recorded to determine the COD profile of the SNCB specimen. This procedure was also performed along the cylinder axis for the SNDB specimen. Note that the maximum COD value is equal to the CMOD value.
Similar to the beam with span/depth=2.5 (Eqn. (3)), the following formulas were selected for the general form of $Y$$Y$YY$Y$ for the SNDB specimen and the SNCB specimen, respectively.

The normalized equations for $Y$$Y$YY$Y$ in Eqns. (15) and (16) were calculated by normalizing the FEM computations-based $K_I$$K_I$K_(I)K_{I}$K_I$ values with the values of $\left(P\sqrt{}\left(\pi a_0\right)/2bd\right)$$\left(P\sqrt{}\left(\pi a_0\right)/2bd\right)$(Psqrt()(pia_(0))//2bd)\left(P \sqrt{ }\left(\pi a_{0}\right) / 2 b d\right)$\left(P\sqrt{}\left(\pi a_0\right)/2bd\right)$ for each specimen geometry, $a/d$$a/d$a//da / d$a/d$, and span/depth ratio. In Tab. 1, the regression coefficients $\left(A_i\right)$$\left(A_i\right)$(A_(i))\left(A_{i}\right)$\left(A_i\right)$ computed using the least squares method are tabulated for each specimen geometry and span/depth ratio. The maximum relative error (accuracy) between the individual FEA analysis data and the above regression formulas was obtained as $0.2\mathrm{\%}$$0.2\mathrm{\%}$0.2%0.2 \%$0.2\mathrm{\%}$ and $0.3\mathrm{\%}$$0.3\mathrm{\%}$0.3%0.3 \%$0.3\mathrm{\%}$ for Eqns. (15) and (16), respectively.
According to the mechanics of structures containing cracks, the displacement expressions of a cracked structure, such as CMOD, depend on both the elastic constants ( $E$$E$EE$E$ and $v$$v$vv$v$ ) of the material and the characteristic size of the structure. Consequently, a general expression of CMOD that satisfies boundary conditions was chosen as follows:

Similar to beam (Eqn. (4)), the following two formulas were selected for the dimensionless $V_1$$V_1$V_(1)V_{1}$V_1$ function of the SNDB specimen and the SNCB specimen, respectively:

The regression coefficients ( $B_i$$B_i$B_(i)B_{i}$B_i$ ) corresponding to the above equations are reported for each specimen geometry and $s/d$$s/d$s//ds / d$s/d$ ratio in Tab. 1. Eqns. (18) and (19) are valid for $0.1\le {\alpha }_0\le 0.9$$0.1\le {\alpha }_0\le 0.9$0.1 <= alpha_(0) <= 0.90.1 \leq \alpha_{0} \leq 0.9$0.1\le {\alpha }_0\le 0.9$ within $0.5\mathrm{\%}$$0.5\mathrm{\%}$0.5%0.5 \%$0.5\mathrm{\%}$ and $0.4\mathrm{\%}$$0.4\mathrm{\%}$0.4%0.4 \%$0.4\mathrm{\%}$ accuracy, respectively. The dimensionless functions derived for cylindrical and cubical specimens are represented by the curved line for $0.1\le {\alpha }_0\le 0.8$$0.1\le {\alpha }_0\le 0.8$0.1 <= alpha_(0) <= 0.80.1 \leq \alpha_{0} \leq 0.8$0.1\le {\alpha }_0\le 0.8$ in Fig. 7, as the values for ${\alpha }_0=0.9$${\alpha }_0=0.9$alpha_(0)=0.9\alpha_{0}=0.9${\alpha }_0=0.9$ are too high. In the same figure, several values determined by Tutluoglu and Keles [14] are also plotted using + markers. It can be seen from this figure that the results of Eqn. (15) exhibit satisfactory accuracy. Note that $95\mathrm{\%}$$95\mathrm{\%}$95%95 \%$95\mathrm{\%}$ confidence interval (CI) values (CI-L lower and CI-U upper) were also presented for each coefficient in Tab. 1.

The COD/CMOD ratio depends on the structure's geometry and loading type but is independent of its size. When this ratio is known for any structure, the critical value of CTOD for the material can be easily determined. Consequently, Eqns. (20) and (21) were selected as the COD profiles for the SNDB specimen and the SNCB specimen, respectively:

Consequently, the $CTO{D}_c$$CTO{D}_c$CTOD_(c)C T O D_{c}$CTO{D}_c$ value for the material can be computed using the values $y=a_0,a=a_c$$y=a_0,a=a_c$y=a_(0),a=a_(c)y=a_{0}, a=a_{c}$y=a_0,a=a_c$, and $\alpha =a_c/d$$\alpha =a_c/d$alpha=a_(c)//d\alpha=a_{c} / d$\alpha =a_c/d$. For all chosen values, Eqns. (20) and (21) achieve an accuracy better than $3.3\mathrm{\%}$$3.3\mathrm{\%}$3.3%3.3 \%$3.3\mathrm{\%}$.
The validity of the above simulations performed was demonstrated in Fig. 8. Two mesh convergence plots for both 2D and 3D cases vs element size ( $1\mathrm{mm},2.5\mathrm{mm}$$1\mathrm{mm},2.5\mathrm{mm}$1mm,2.5mm1 \mathrm{~mm}, 2.5 \mathrm{~mm}$1\mathrm{mm},2.5\mathrm{mm}$, and 5 mm ) were illustrated for $a=0.1,0.3,0.5,0.7$$a=0.1,0.3,0.5,0.7$a=0.1,0.3,0.5,0.7a=0.1,0.3,0.5,0.7$a=0.1,0.3,0.5,0.7$, and 0.9 for SNCB specimens with $s/d=0.8$$s/d=0.8$s//d=0.8s / d=0.8$s/d=0.8$ in Fig. 8a. It is clearly seen from the figure that mesh structures with element sizes of 1 mm and 2.5 mm (corresponding to 100 and 40 elements along the crack line, respectively) gave very close results. In these comparisons, the maximum relative error calculated between 1 mm and 2.5 mm meshes was calculated as $1.3\mathrm{\%}$$1.3\mathrm{\%}$1.3%1.3 \%$1.3\mathrm{\%}$ for 2D mesh and $1.4\mathrm{\%}$$1.4\mathrm{\%}$1.4%1.4 \%$1.4\mathrm{\%}$ for 3D mesh, respectively, while these errors were $11.7\mathrm{\%}$$11.7\mathrm{\%}$11.7%11.7 \%$11.7\mathrm{\%}$ and $7.9\mathrm{\%}$$7.9\mathrm{\%}$7.9%7.9 \%$7.9\mathrm{\%}$ for the comparison between 1 mm and 5 mm meshes. As a result, it can be concluded that, as mentioned above, using 100 finite elements along the notch line is sufficient to minimize calculation errors. In Fig. 8b, the distributed load effect was investigated for the SNDB specimen. It is clearly seen from the figure that the relative errors between concentrated loading and distributed loading were smaller than $0.2\mathrm{\%}$$0.2\mathrm{\%}$0.2%0.2 \%$0.2\mathrm{\%}$ and $1.6\mathrm{\%}$$1.6\mathrm{\%}$1.6%1.6 \%$1.6\mathrm{\%}$ for ratios of the width of the distributed load to the specimen depth $b=0.05$$b=0.05$b=0.05b=0.05$b=0.05$ and 0.10 , respectively. As mentioned above, this effect may be ignored for bending specimens. Fig. 8c presents contour plots with various radii ( $r$$r$rr$r$ ) and corresponding $Y$$Y$YY$Y$ values computed to show path independence for the J-integral of the 2D-SNCB specimen with $a/d=0.5$$a/d=0.5$a//d=0.5a / d=0.5$a/d=0.5$. Since the differences for the Y values are very small here, it can be said that the calculated J-integrals are path independent.