Theory

The nonlinear stress-strain response in fine-grained graphite was characterized by the damage law (Eq. 1) given in Sect. 3.1; however, it did not include the failure criterion. The critical stress or critical strain were often employed as failure criteria in graphite [26, 28], but they may be insufficient to characterize its fracture behavior, specifically stable crack propagation and toughening mechanisms [49, 50]. Some simulations indicated that using fracture energy as the damage criterion ensured good consistency with the experimental data [40, 51]. Therefore, the fracture energy G_f from fracture mechanics was chosen as the failure criterion to supplement the previous damage model.

To implement finite element analysis, based on previous studies on the nuclear graphite damage evolution and fracture process [4, 11], the degradation process of the graphite properties under stress was defined according to the following three states (as shown in Fig. 7(a)) in the FEM:

Fig. 7Full-screen View

Schematic of the damage-fracture model: (a) the three defined states of the damage-fracture model; (b) bilinear TSC and ultimate relative displacement ωult

First state: The material is in the diffuse damage state when its maximum principal stress is less than its tensile strength (${\sigma }_1<f_{\text{t}}$). In this state, micro-cracks begin to appear and initial defects and pores enlarge in the material owing to stress. These changes in the original microstructure in graphite create a nonlinear stress-strain relationship [11]. In the FEM, this process is characterized by the damage law (Eq. 12) of Young’s modulus with strain. The stress increases in a nonlinear monotonic manner as the strain increases in this state.

Second state: The material is in the damage localization state when the stress level reaches its tensile strength (${\epsilon }_1\ge \frac{f_{\text{t}}}{E_s}$, where E_s is the secant value of Young’s modulus at tensile strength). In this state, microcracks coalesce and become interconnected with the pores and other defects, and the stress between graphite particles is transferred through frictional interactions and mechanical interlocking across the microcrack surfaces, otherwise known as the bridging effect. In terms of mechanics, it appears that the material’s ability to sustain stress decreases with increasing tensile strain. In FEM, this behavior is described by a tensile softening curve (TSC). To simplify it without losing the features [36], a bilinear form of the TSC is adopted (Fig. 7(b)).

Third state: As the microcrack localization area develops, the cohesive stress between microcracks gradually reduces to zero. At this point, the microcracks evolve into stress-free macro cracks and continue to propagate with the ongoing deformation. In the FEM, when the cumulative energy dissipation of an element reaches the fracture energy G_f, that element is deleted to simulate crack propagation. To reduce the influence of element size on the calculation of energy consumption, relative displacement is used to calculate the energy dissipation of an element instead of strain [41], as follows: $\omega =l_{\text{CELENT}}\cdot \epsilon$ (7) where l_CELENT denotes the characteristic length of the element. For linear two-dimensional elements, it equals the square root of the element area (Eq. 9), which is calculated using ABAQUS and supplied to the subroutine. $l_{\text{CELENT}}=\sqrt{Are{a}_{\text{element}}}$ (8) The ultimate relative displacement ω_ult of an element is defined as the relative displacement ω at which the element’s cumulative energy dissipation (shown by the enclosed area in Fig. 7(b)) achieves the critical value, G_f. Subsequently, by dividing this ultimate relative displacement ω_ult by the element’s characteristic length l_CELENT, the ultimate strain value ε_ult that VUMAT uses to determine whether to delete an element can be obtained, with the specific calculations shown in Eq. 9: ${\epsilon }_{\text{ult}}=\frac{{\omega }_{\text{ult}}}{l_{\text{CELENT}}}=\frac{2\left(G_{\text{f}}-\frac{f_{\text{t}}^2}{2E_{\text{s}}{l}_{\text{CELENT}}}\right)}{f_{\text{t}}\left(R_1+R_2\right)l_{\text{CELENT}}}$ (9) where E_s is the secant value of Young’s modulus at the point where the maximum principal stress reaches the tensile strength f_t, G_f is the fracture energy of the material, and R₁ and R₂ shape the TSC (Fig. 7(b)). E_s can be determined using Eq. 1, and G_f, R₁, and R₂ must be obtained from fracture tests.

The above process only considers the Mode I (opening type) fracture condition. In applications where Mode II fractures predominate, the stress criterion f_t and fracture criterion G_f proposed above should be adjusted to adapt to the characteristics of such fractures [12].

Determination of fracture properties

To acquire the fracture parameters necessary for the aforementioned model, i.e. the fracture energy G_f and the shape parameters R₁ and R₂ of the TSC (Fig. 7(b)), an FEM integrating the above theory was established by VUMAT. This FEM had the same boundary conditions as the three-point bending test in Sect. 2.2, and its simulated load-displacement response was compared with the experimental results to optimize the values of the fracture parameters.

As shown in Fig. 8, quadrilateral elements and a two-dimensional plane stress assumption were used in the FEM. In the predicted crack propagation area, square elements with a size of 0.2 mm were employed to refine the mesh and decrease the model’s sensitivity to the element shape and size [40, 52]. Prior to the stress reaching the tensile strength, the damage law of Young’s modulus adhered to the inverse results obtained by the method described in Sect. 3.1. Poisson’s ratio was set to 0.14, and the tensile strength was set to 27.6 MPa [53].

Fig. 8Full-screen View

(Color online) Finite element model of the center-notched beam in three-point bending

In FEM, applying IG11 graphite’s typical Young’s modulus value of 10 GPa could lead to inaccuracies in parameter optimization. This is because Young’s modulus can vary between individual graphite specimens owing to the significant spatial heterogeneity of the properties in bulk graphite [54, 55]. Therefore, to input a more accurate Young’s modulus into the model, the Young’s modulus of the specimen was calculated from the experimental P-CMOD curve using the formulas (Eq. 10) provided in Tada’s Stress Intensity Factor Handbook [56]. For the standard three-point bending notched beam with a span-to-depth ratio of 4, these formulas were as follows: $W_{\text{CMOD}}=\frac{24\text{Pa}}{WBE}g\left(\frac{a}{W}\right)$ (10) $\begin{array}{ll}g\left(\frac{a}{W}\right)\hfill & =0.76-2.28\left(\frac{a}{W}\right)+3.87{\left(\frac{a}{W}\right)}^2\hfill \\ \hfill & -2.04{\left(\frac{a}{W}\right)}^3+\frac{0.66}{{\left[1-\left(\frac{a}{W}\right)\right]}^2}\hfill \end{array}$ (11) where W and B represent the depth and thickness of the beam, respectively, a is the notch depth, P is the load, W_CMOD represents the crack mouth opening displacement, and E is the Young’s modulus, which can be back-calculated by substituting the experimental P-CMOD data into this formula.

Linear elastic model results

As expected [57], all the disc specimens failed owing to cracking at their center. The failure loads and tensile strengths of three disc specimens with different diameters are listed in Table 4. The findings indicate that the tensile strength was affected by size, with larger specimens exhibiting lower tensile strengths. This trend is consistent with previous research where the strength data were analyzed using the Weibull weakest link theory [58].

To assess the suitability of the linear elastic model for simulating graphite deformation, using a 15 mm-diameter specimen as a representative example, the principal strain fields at failure load were calculated using a linear elastic FEM with a Young’s modulus of 10.5 GPa and compared with the strain fields measured using DIC in experiments, as shown in Fig. 9. The comparison revealed qualitatively similar strain distributions; however, quantitative differences were observed. More specifically, strain levels measured using DIC were higher than those obtained using the FEM. For the maximum principal strain, the average DIC-measured value (0.0038 ± 0.0029) was approximately 39.47% greater than the average simulated value (0.0023 ± 0.0019). Similarly, for the minimum principal strain, the average DIC-measured absolute value (0.0115 ± 0.0061) was approximately 35.65% greater than the simulated value (0.0074 ±0.0035). The primary reason for this mismatch is the nonlinear stress-strain behavior intrinsic to graphite [11], although it is less apparent in fine-grained graphite than in coarser grains. Thus, simulating graphite using a linear elastic model may yield an underestimation of the real deformation.

Fig. 9Full-screen View

(Color online) Comparison of principal strain fields obtained from DIC and FEM: (a) maximum principal strain distribution; (b) minimum principal strain distribution

Damage evolution law

By applying the procedure in Sect. 3.1 to the experimental data from one specimen per size, the damage evolution in the graphite specimens was characterized by quantifying the change in Young’s modulus with increasing strain, as shown in Fig. 10. The Young’s modulus was found to reduce with increasing strain level, and this damage behavior varied for compressive and tensile strains. Specifically, the Young’s modulus reduced more in tension at the same strain level. Additionally, the initial Young’s modulus obtained had an average value of 10.97 GPa, which was close to the 10.5 GPa provided by the supplier. As expected, the size of the specimens did not significantly affect the inversion results because damage law is an inherent property of graphite.

Fig. 10Full-screen View

(Color online) Inverted damage laws of Young’s modulus from three sizes of specimens: (a) under tensile strain; (b) under compressive strain

To simplify the use of the identified damage law in finite element analysis, the segmented results in Fig. 10 were fitted by the continuous equation in Eq. 12. $E=E^{\text{ini}}\left(1-a{\left|\epsilon -{\epsilon }_{\text{e}}\right|}^b\right)$ (12) where a and b are parameters to be fitted, and ε_e is the end value of the initial linear elastic segment. Their values are presented in Fig. 10.

Comparison between damage FEMs and experiments

To verify the accuracy of the damage law obtained using the inversion analysis, the fitted damage law (Eq. 12) was used in the finite element analysis to simulate disc specimens of three sizes, and the FEM-calculated strains were compared with the DIC-measured strains at the failure load, as shown in Fig. 11. The strain distributions and values of FEMs and DIC were almost consistent, exhibiting no significant variations with the size of the specimen. This indicates that the FEMs with the damage law can effectively simulate non-linear deformation in graphite specimens. It should be noted that the edge information of the specimen was ignored in the DIC analysis, causing slight differences with the FEM simulation results in the edge region. To quantify the differences between the simulation and the experiment, the strain error values, calculated by subtracting the simulated values from the experimental values, were plotted in frequency distribution graphs (on the right of Fig. 11). They indicate that the strain errors fluctuate around zero and are well fitted by a Gaussian function. This fluctuation may be due to the microscopic heterogeneity of graphite, tolerance of the specimen and the fixture, and combined measured error. This phenomenon has also been observed in previous graphite research [17].

Fig. 11Full-screen View

(Color online) Strain distribution of FEM and DIC and the corresponding strain error frequency distribution. (a) Maximum principal strain of the disc specimen with a radius of 10 mm; (b) Minimum principal strain of the disc specimen with a radius of 10 mm; (c) Maximum principal strain of the disc specimen with a radius of 7.5 mm; (d) Minimum principal strain of the disc specimen with a radius of 7.5 mm; (e) Maximum principal strain of the disc specimen with a radius of 5 mm; (f) Minimum principal strain of the disc specimen with a radius of 5 mm

To visually depict the advantage of this damage model over the linear elastic model in simulating graphite, the linear elastic model and the DIC-measured results were compared in terms of the transverse deformation U, longitudinal deformation V, and the maximum principal strain at the center position in the specimens, as shown in Fig. 12. It was demonstrated that the damage model matched the experimental data quite well, whereas the linear elastic model significantly underestimated the actual deformation and strains of the specimens. The damage model can accurately simulate both tensile and compressive responses, which may be difficult for a single damage variable model.

Fig. 12Full-screen View

(Color online) Strain distribution of FEM and DIC and the corresponding strain error frequency distribution. (a) Specimen with a radius of 10 mm; (b) Specimen with a radius of 7.5 mm; (c) Specimen with a radius of 5 mm

Young’s modulus

By inserting the CMOD values measured by the clip gauge and the corresponding applied load into the Eq. 10, the global Young’s modulus of the center-notched beam specimen was determined, as shown in Fig. 13. Because the notch length input into Eq. 10 was always the initial length, once the notch extended, the calculated Young’s modulus no longer reflected the actual material properties, but it could still be used to represent the specimen’s overall stiffness. According to Fig. 13, the initial Young’s modulus of the specimen was 11.29 GPa, and this value was incorporated into the subsequent inversion of fracture parameters. Additionally, when the load reached slightly less than the peak load (93% of the peak load), a sharp decrease in the global Young’s modulus of the specimen was observed. This decline may be attributable to the propagation of the main crack, which weakens the overall stiffness of the specimen. In another three-point bending test with the center-notched IG11 graphite beam [59], it was found that the critical load for the onset of crack extension was slightly smaller than the peak load. Thus, Eq. 10 may help identify the initiation point of the crack in the specimen.

Fig. 13Full-screen View

The overall Young’s modulus of the specimen degrades with displacement loading

Fracture parameters

The fracture parameters G_f, R₁, and R₂ in the FEM were optimized by comparing the simulated and actual load-displacement curves in the three-point bending test, including P-CTOD, P-CMOD, and P-δ (deflection) data, to obtain the values of the fracture parameters that offer the best match between the simulation results and the experimental outcomes. The optimization results indicated that when G_f was set to 176 N/m, R₁ to 0.397, and R₂ to 0.242 (Eq. 9), the simulated load-displacement curves exhibited a high degree of consistency with the experimental curves, as shown in Fig. 14. These optimized fracture parameters were close to those obtained in previous studies on fine-grained graphite [36]. Ultimately, all parameters required by the damage-fracture model for fine-grained graphite were determined and are displayed in Table 5.

Fig. 14Full-screen View

Comparison between the experimental and optimized simulated load-displacement curves: (a) P-CTOD; (b) P-CMOD; (c) P-δ

The specific fracture energy ${G^{\prime }}_{\text{f}}$ was calculated as 189 N/m from the simulated P-δ curve (Fig. 14(c)) using Eq. 13; this differed slightly from the defined fracture energy G_f of 176 N/m that was input into the model, with their ratio being approximately 1.07. Another study on damage models with a bilinear TSC discovered that this ratio, which is commonly near 1.0, tends to increase with the brittleness of the material [52]. ${G^{\prime }}_{\text{f}}=\frac{1}{A_{\text{f}}}{\displaystyle \underset{0}{\overset{{\delta }_{\text{ult}}}{\int }}}P\text{d}\delta$ (13) where FV is the applied load, δ denotes the vertical displacement, δ_ult denotes the ultimate vertical displacement, and A_f denotes the fracture surface, which is 25 mm × 30 mm in this test.

Simulation of fracture behavior

Figure 15 shows the distributions of the tensile damage factors in the FEM at different loading stages. This factor reflects the degree of damage caused by tensile strain, with values ranging from 0 (no damage) to 1 (complete damage). When the tensile damage factor reached 1, the corresponding elements were removed in the simulation, creating blank areas in Fig. 15. The damage distributions can be divided into three regions: front zone, wake zone, and macro-crack zone (Fig. 15(c)). These regions correspond to the three states described in Sect. 3.2.1. The size of the entire damage region was measured (Fig. 15), indicating that the size of the region initially increased from small to fully developed and then gradually decreased. This process was consistent with the changes in the fracture process zone (FPZ) observed in many graphite tests [35, 51, 60, 61]. That is, as the stress level in the specimen increases, the FPZ gradually enlarges until it is fully developed; when the crack propagates close to the specimen edge, the FPZ starts to shrink owing to the boundary effects.

Fig. 15Full-screen View

(Color online) Tensile damage factor distributions and the damage regions size at (a) 76% P_max (pre-peak); (b) P_max (peak); (c) 35% P_max (post-peak); and (d) 5% P_max (post-peak)

Further observations (Fig. 15) indicated that when the load exceeded 76% P_max (Fig. 15(a)), the damage began to localize, and the wake zone formed. A previous study [62] on IG11 three-point bending beams found that the P-CTOD curve transitions from linear to non-linear at approximately 84% P_max, which may be explained by the occurrence of damage localization in the model. In addition, when the FPZ was fully developed (Figs. 15(b) and 10(c)), the observed wake zone length was approximately 2.6 mm, which was close to the value of approximately 2.5 mm observed in previous graphite single edge notched beams [63]. However, it should be noted that the size of the FPZ often varies depending on the shape of the specimen, so it should not be considered an inherent property of graphite. In summary, the damage-fracture model effectively simulated some fracture characteristics observed in the experiment, including stable crack propagation, FPZ changes with crack propagation, and the front zone and wake zone in the process zone.

In previous fracture behavior simulations of graphite notched specimens [14, 40], although the damage evolution prior to reaching tensile strength was ignored, the simulation results still closely matched the experimental load-displacement response. This could be because the notch’s sharp shape concentrated severe stress in a limited area around the notch tip, limiting the influence of diffuse damage on the specimen’s overall response. If the notch tip were blunted, the influence of the front zone might be enhanced. Furthermore, this omission could overlook the damage in the front zone, which seems to diverge from the strain field concentration around the notch tip observed in experiments [60].

Finite element model

To enable comparison between the simulation results and previous experimental data [20], FEMs were created for specimens with geometric dimensions identical to those of the L-shaped graphite specimens described in Ref. [20]; the specific dimensions are shown in Fig. 16(a). In the FEM, quadrilateral elements and a two-dimensional plane stress assumption were used, and the mesh at the corner was refined to 0.3 mm. By coupling the midpoint of the holes, the hole at the bottom was fixed for both lateral and vertical displacements, whereas the hole at the top was fixed for lateral displacement and subjected to a vertical displacement load. The mesh and boundary conditions are shown in Fig. 16(b).

Fig. 16Full-screen View

(Color online) 2D finite element model of the L-shaped specimen: (a) model dimensions; (b) model mesh and its boundary conditions

To investigate the diffuse damage influence on the model, each dimensional model was assigned two material property sets: one encompassing damage and fracture parameters as detailed in Table 5, and the other considering only fracture parameters with damage parameters set to zero. Both sets adopted a Young’s modulus and Poisson’s ratio of 10 GPa and 0.14, respectively, representing typical values for IG11 graphite.

Prediction of failure load

The failure loads (peak load) predicted by both groups of models are shown in Fig. 17, indicating an upward trend with an increase in fillet radius. This trend was attributed to the decrease in stress concentration as the fillet radius increased [20]. Moreover, the model accounting for damage and fracture exhibited a higher load-bearing capacity and was generally closer to the experimental values, whereas the model accounting for only fracture displayed an increasing discrepancy with the experimental values as the fillet radius increased. (Fig. 17). This is due to the diffuse damage in the fillet area that blunted the stress concentration, preventing the local strain energy from prematurely reaching the fracture energy and triggering the appearance of macro cracks in the models.

Fig. 17Full-screen View

Failure loads predicted by the FEM and obtained from experiments [20] with varying fillet radii

Simulation of damage behavior

From the tensile damage distribution of the damage-fracture model at the failure load (Fig. 18), it can be observed that the area of damage distribution expanded as the fillet radius increased, enhancing the blunting effect. This may explain the trend observed in Fig. 17 where the difference in failure load between the damage-fracture model and the fracture model increased with the increase in fillet radius.

Fig. 18Full-screen View

(Color online) Distribution of the tensile damage factor in the FEM at the failure load with fillet radii of (a) 0.8 mm; (b) 1.4 mm; (c) 2.0 mm; and (d) 3.0 mm

For local strains, taking the model with a fillet radius r = 1.4 mm as an example, the average tensile strain at the fillet (i.e., vertical strain) was extracted and compared with experimental data, as shown in Fig. 19. The results indicated that the predictions of the damage-fracture model aligned well with the measured experimental values, whereas the fracture model underestimated the actual strain within 50% to 99% of the failure load, and when the stress in the fracture model transitioned into the softening state, the strain values rapidly approached the experimental values. This indicates that the fracture model can simulate the specimen deformation well after softening occurs but exhibits errors between half the failure load and the onset of softening.

Fig. 19Full-screen View

(Color online) Variation in average tensile strain at the fillet with load (fillet radius of 1.4 mm)

It should be noted that the damage evolution parameters used in the damage-fracture model were obtained from the inversion in CDI-1D graphite, but the simulation results were consistent with the experimental results of IG-11 graphite [20]. This may indicate that graphite samples with similar particle sizes and forming methods will follow similar damage laws.

Simulation of fracture behavior

Both sets of models exhibited the failure mode of quasi-brittle materials [11], characterized by stable crack propagation and a descending load-displacement curve after the peak load, as shown in Figs. 20 and 21. From Fig. 20, the load-displacement curves of the models with a larger fillet radius were steeper and reached higher levels, indicating greater stiffness and load-bearing capacity. However, once the crack extended, the curves rapidly converged. To analyze this, Fig. 21 depicts the tensile damage distribution in the damage-fracture models during crack propagation. The varying fillet radius led to different damage distributions around the fillets; however, after the crack had extended by a certain length, the damage distributions around the stress concentration areas near the crack tips became nearly identical. To quantify the models’ damage, the macro crack length (ac) and the length and the width of the process zone (l_FPZ and w_FPZ, respectively) were measured, with the results shown in Table 6. The measurements indicated that the crack lengths and FPZ dimensions of models with different fillet radii were almost equivalent under identical loads. This may explain why the load-displacement curves ultimately converged. Furthermore, it was noted that the size of the FPZ near the specimen edges was reduced owing to boundary effects, consistent with patterns observed in previous simulations with graphite beams, but the dimensions of the FPZ varied because different boundary conditions.

Fig. 20Full-screen View

(Color online) Load-displacement curve simulated by FEM with varying fillet radii: (a) fracture model; (b) damage-fracture model

Fig. 21Full-screen View

(Color online) Distribution of the tensile damage factor in FEM with varying fillet radii at (a) P = 100 N (post-peak); (b) P = 20 N (post-peak)

The specified fracture energy was calculated by dividing the area under the load-displacement curves in Fig. 20 by the fracture area (Eq. 13). It was compared with the defined fracture energy that was input into the FEM, as shown in Fig. 22. It was noted that the specified fracture energy of the fracture model fluctuated slightly above the defined fracture energy; however, the specified fracture energy of the damage-fracture model increased as the fillet radius increased. This increment could be attributed to the fact that, in the damage-fracture model, the dissipated energy was used not only for removing elements to form fracture surfaces but also in part for the evolution of diffuse damage. Furthermore, a larger fillet radius corresponds to a wider range of diffuse damage (Fig. 21); therefore, more energy was required to extend the same fracture surface in the damage-fracture models with larger fillet radii. As the defined fracture energy only accounts for the energy necessary to remove elements, it is less than the specified fracture energy determined from the damage-fracture model and close to that of the fracture model.

Fig. 22Full-screen View

Specific fracture energy of fracture models and damage-fracture models with different fillet radii