Mixed-Mode Fracture Behavior of Penta-Graphene: A Molecular Dynamics Perspective on Defect Sensitivity and Crack Evolution

1. Introduction

In recent years, carbon allotropes have become a cornerstone of modern materials science research due to their diverse structures and unique properties. These allotropes exhibit remarkable physical, chemical, and electronic characteristics, making them indispensable in various technological fields. Diamond, graphite, fullerenes, carbon nanotubes, and graphene have significantly influenced advancements in medicine, electronics, energy storage, and composite materials [1,2,3,4].

Among these materials, graphene, a two-dimensional carbon allotrope, has garnered extensive attention due to its exceptional properties, including superior mechanical strength [5,6], excellent electrical conductivity [7,8], and outstanding thermal conductivity [9,10]. These properties have enabled graphene to revolutionize applications in nanoelectronics, flexible electronics, energy storage devices, and high-performance composites. However, despite these advantages, graphene has inherent limitations that restrict its widespread implementation in certain fields. One of the primary drawbacks is its lack of an intrinsic bandgap, which hinders its direct application in semiconductor technologies without additional modifications [11,12]. Additionally, challenges related to large-scale production and consistent quality [13,14], along with its sensitivity to chemical degradation under specific conditions [15,16,17], have led researchers to explore alternative two-dimensional materials with improved or complementary properties [18,19,20,21].

One such promising material is penta-graphene, a theoretically proposed two-dimensional carbon allotrope composed entirely of pentagonal rings [22]. Unlike graphene, penta-graphene possesses a natural bandgap [22,23,24], making it highly suitable for semiconductor and nanoelectronic applications. Furthermore, penta-graphene retains excellent thermal properties [25,26,27] and exhibits remarkable structural robustness [28,29,30], making it an attractive candidate for next-generation electronic and optoelectronic devices. Its potential applications extend beyond electronics, as its unique structure offers promising implications in energy storage, sensor technologies, and mechanical reinforcement in composite materials.

The atomic configuration of penta-graphene, illustrated in Figure 1, showcases a unique hybridization of ${\mathrm{sp}}^2$ and ${\mathrm{sp}}^3$ carbon atoms, structured into three distinct layers—top, middle, and bottom. The outer layers consist solely of ${\mathrm{sp}}^2$ hybridized atoms, while the middle layer contains ${\mathrm{sp}}^3$ hybridized atoms, providing structural stability and tunability through processes such as hydrogenation—a feature not present in graphene [23,31,32,33]. This adaptability through external stimuli like strain engineering and chemical functionalization augments its potential for integration into tunable electronic devices and quantum materials research. Furthermore, penta-graphene demonstrates enhanced chemical reactivity relative to graphene, primarily due to the inherent strain within its pentagonal rings. This reactivity aids in functionalization and chemical modifications, expanding its applicability across various fields including nanoelectronics, energy solutions, and composite material development.

In addition to its electronic versatility, penta-graphene exhibits remarkable mechanical characteristics. Penta-graphene also exhibits intrinsic mechanical anisotropy due to its pentagonal symmetry and mixed hybridization. The interplay between ${\mathrm{sp}}^2$ and ${\mathrm{sp}}^3$ bonding results in an optimal balance of rigidity and flexibility, leading to distinctive mechanical characteristics [34,35,36,37,38,39]. These properties suggest that penta-graphene could serve as a reinforcement material in high-performance composites, aerospace applications, and protective coatings. One of penta-graphene’s most unique mechanical properties is its negative Poisson’s ratio [37,40], which means that under tensile strain, it expands laterally instead of contracting. This counterintuitive behavior, often referred to as auxetic behavior, arises from the flattening of its layered structure under applied strain [41]. Such behavior is advantageous for applications requiring enhanced mechanical stability and energy dissipation, including impact-resistant materials and biomedical scaffolds.

Despite its theoretical promise, penta-graphene’s mechanical and fracture properties remain poorly understood. Experimental validation is currently unfeasible due to the lack of successful synthesis methods, though theoretical studies predict penta-graphene to be metastable [22]. In the absence of experimental data, molecular dynamics (MD) simulations serve as a critical tool for investigating penta-graphene’s mechanical behavior at the atomic scale [34]. By leveraging computational approaches, researchers can explore the deformation mechanisms, failure modes, and mechanical responses of penta-graphene under various loading conditions.

This study employs MD simulations to examine the mechanical and fracture response of monolayer penta-graphene under different loading conditions. Specifically, we analyze the effects of load angles, crack defects, and point defects on the material’s failure mechanisms. Additionally, we use MD simulations to determine the critical stress intensity factor for penta-graphene under mixed mode I and mode II loading.

2. Atomic Structure of Penta-Graphene

Penta-graphene features a distinctive pentagonal lattice structure, diverging from the hexagonal symmetry characteristic of conventional graphene. This arrangement, based on the Cairo pentagonal tiling, forms a periodically repeating network of pentagons, as illustrated in Figure 1. The unique structural motif is expected to have a profound impact on the material’s mechanical, electronic, and chemical properties.

The atomic structure of penta-graphene, presented in both top-down and side-view projections in Figure 1, reveals a mixed hybridization state of carbon atoms. The material consists of both ${\mathrm{sp}}^2$- and ${\mathrm{sp}}^3$-hybridized carbon atoms, arranged in a three-layered configuration. The side-view projection highlights this structure, showing a central layer of ${\mathrm{sp}}^3$-hybridized carbon atoms (red) intercalated between two layers of ${\mathrm{sp}}^2$-hybridized carbon atoms (blue and gray). The vertical displacement of ${\mathrm{sp}}^3$-hybridized atoms relative to the ${\mathrm{sp}}^2$ layers results in an out-of-plane thickness of 1.2 Å.

Within this framework, the ${\mathrm{sp}}^2$-hybridized atoms form planar bonds within the pentagonal rings, whereas the ${\mathrm{sp}}^3$-hybridized atoms establish out-of-plane bonds. This layered configuration alleviates strain associated with non-hexagonal bonding, thereby enhancing structural stability. The combination of intrinsic stability, unique bonding characteristics, and three-dimensional layering imparts penta-graphene with both high stiffness and mechanical flexibility, making it a promising candidate for advanced structural and mechanical applications.

The unit cell, highlighted by a black square in Figure 1, features a lattice parameter of 3.64 Å, corresponding to the lengths of the unit vectors ${\mathit{a}}_{\mathbf{1}}$ and ${\mathit{a}}_{\mathbf{2}}$, which define the dimensions along the X and Y axes, respectively. The unit cell contains six carbon atoms. The atomic composition maintains an ${\mathrm{sp}}^3$:${\mathrm{sp}}^2$ ratio of 1:2, with two ${\mathrm{sp}}^3$ and four ${\mathrm{sp}}^2$ carbon atoms per unit cell. This specific stoichiometry plays a crucial role in determining the material’s behavior. The bond lengths, measuring 1.55 Å for ${\mathrm{sp}}^3$–${\mathrm{sp}}^2$ bonds and 1.34 Å for ${\mathrm{sp}}^2$–${\mathrm{sp}}^2$ bonds, reflect the distinct bonding environments within the lattice. The bond angle between ${\mathrm{sp}}^2$–${\mathrm{sp}}^3$–${\mathrm{sp}}^2$ is 134.2 °.

A summary of key structural and bonding parameters is provided in Table 1. Future research should focus on the synthesis and characterization of penta-graphene to validate these theoretical predictions and explore its potential applications.

3. Methods

To investigate the behavior of penta-graphene, tensile tests are conducted using the MD simulation package LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [42,43].

When selecting an empirical potential to model the interactions between carbon atoms in penta-graphene, Tersoff potentials are commonly utilized due to their capability to accurately simulate the complex interactions that are crucial for capturing the dynamic behavior of carbon-based materials under various stress conditions. Among various studies, it has been established that the Tersoff potential developed by Erhart and Albe in 2005 offers a suitable depiction of the mechanical and structural behaviors of penta-graphene [44]. This potential considers both the bond angle and the distance between atoms, which are vital for predicting the material’s behavior under tension, compression, and bending. The Tersoff potential models these interactions by considering the electronic environment of each atom, thereby allowing it to simulate how the material behaves under tension, compression, and bending. This is vital for accurately predicting fracture mechanics, as the potential can adapt to the changes in atomic configurations that occur near crack tips and other stress concentrators.

Specifically, the research by Winczewski et al. evaluated 14 different empirical potentials, including three Tersoff-type, two Stillinger-Weber-type, three REBO models, and three parameterizations of the REAX force field [45]. This comprehensive study highlighted that only the Tersoff-type potential by Erhart and Albe accurately reflects the characteristics of penta-graphene when compared to ab initio data, primarily because the potential’s parameters are derived from quantum mechanical calculations, ensuring that they are grounded in the fundamental physics of atomic interactions. Furthermore, the Tersoff potential has proven effective not only in modeling mechanical properties but also in calculating the stress intensity factor for various two-dimensional materials [46,47,48]. Due to its reliability in capturing mechanical and fracture behavior of carbon-based materials, we have employed the 2005 Erhart and Albe Tersoff potential for our simulations [44].

The simulation consists of two major stages: an equilibration stage and a deformation stage, each involving several specific steps to ensure an accurate representation of the system’s behavior. The equilibration stage begins with a minimization step to relax the system and eliminate any unfavorable interactions or strains. Following this, a constant pressure and temperature ensemble (NPT) is applied to the system for 10 picoseconds, with the temperature set at 1 Kelvin and the pressure set at 1 bar. This initial NPT phase allows the system to reach an equilibrium state under controlled pressure and temperature conditions. Subsequently, the temperature is gradually increased from 1 Kelvin to 300 Kelvin over 25 picoseconds, ensuring a smooth transition to the desired simulation temperature. The system is then held in an NPT ensemble for an additional 25 picoseconds to allow for further equilibration, during which the pressure and temperature are maintained constant. To stabilize the system, it is transferred to an NVE (constant number of atoms, volume, and energy) ensemble, where it is allowed to evolve for an additional 10 picoseconds, enabling any remaining fluctuations to dissipate.

In the deformation stage, the structure is subjected to tension under an NPT deformation. An engineering strain rate of 0.001ps⁻¹ is applied, causing the structure to deform and expand gradually. The deformation process is carried out for a total of 1000 picoseconds, during which the structure’s width is doubled relative to its original dimensions. The simulation may be stopped early upon reaching the point of total structural failure to save system resources. Data from the deformation stage is used to investigate the structure’s response to applied strain, providing insights into its mechanical properties.

Stress-strain curves can be obtained from this process by recording the stress in the material during loading. The modulus of elasticity can then be calculated from the initial linear elastic region of the stress-strain curve. The initial slope of this curve over the strain range from 0 to 0.1 is used to determine the modulus of elasticity, and the highest stress is used to determine the fracture strength of the material.

4. Impact of Defects on Mechanical Properties of Penta-Graphene

The mechanical properties of penta-graphene, like those of other two-dimensional materials, are highly sensitive to structural imperfections [49,50,51]. Among these, point defects such as vacancies and initial cracks play a particularly critical role. These localized disruptions can significantly reduce the material’s tensile strength, elasticity, and fracture toughness, despite the inherent robustness of its pentagonal lattice structure. Given that such defects are inevitable during synthesis, processing, and operational use, understanding their impact is crucial for both fundamental research and practical applications.

Voids, or atomic vacancies, are an inherent aspect of real materials, arising from imperfections during synthesis, fabrication, or operational wear. These missing atoms disrupt the local bonding environment, creating stress concentrations around defect sites and impairing the material’s ability to distribute mechanical loads uniformly. Such irregularities weaken the structural integrity of penta-graphene, reducing its load-bearing capacity and altering its mechanical response. The distribution and density of voids further influence the material’s elasticity and fracture resistance, making their characterization essential for predicting mechanical performance.

Similarly, initial cracks—whether microscopic or larger in scale—serve as critical stress concentrators, significantly affecting mechanical behavior. These preexisting flaws act as nucleation sites for crack propagation under applied loading, accelerating structural failure. Unlike vacancies, which primarily disrupt localized stress distribution, cracks compromise the overall integrity of the material by facilitating mechanical energy accumulation and fracture progression. The combined presence of voids and cracks exacerbates structural weaknesses, necessitating a comprehensive understanding of their impact for material optimization.

In this study, we investigate the influence of these defects on the stress-strain behavior, elastic modulus, and strength of penta-graphene. A series of simulations are conducted to systematically examine the effects of void percentage and distribution on mechanical properties. Additionally, we analyze how varying crack size and orientation impact the material’s response under mechanical loading.

4.1. Impact of Point Defects on the Mechanical Properties

In penta-graphene, the coexistence of ${\mathrm{sp}}^2$ and ${\mathrm{sp}}^3$ hybridized carbon atoms creates a complex and mechanically robust structure that combines planar and out-of-plane bonding characteristics. However, the introduction of point defects—such as vacancies, substitutions, and interstitial atoms—can disrupt this delicate balance and significantly alter the material’s mechanical behavior. It is essential to study point defects associated with both ${\mathrm{sp}}^2$ and ${\mathrm{sp}}^3$ hybridized carbon atoms because their effects on the mechanical integrity of penta-graphene differ due to the nature of their bonding.

Defects in ${\mathrm{sp}}^2$-hybridized regions primarily impact the planar bonding network, which can lead to reduced flexibility and increased brittleness. Since ${\mathrm{sp}}^2$ bonds are generally stronger but less flexible than ${\mathrm{sp}}^3$ bonds, a defect in this region could localize stress and act as a nucleation point for crack initiation and propagation. Conversely, defects in ${\mathrm{sp}}^3$-hybridized regions may introduce distortions in the three-dimensional stacking and interlayer interactions, potentially affecting the out-of-plane rigidity and load transfer mechanisms that contribute to penta-graphene’s mechanical resilience. These disruptions can alter how the material accommodates deformation, impacting its ability to distribute stress uniformly and leading to localized structural instabilities that may compromise its overall mechanical performance.

To thoroughly evaluate the impact of point defects in penta-graphene, void defects are independently introduced into both the ${\mathrm{sp}}^2$ and ${\mathrm{sp}}^3$ hybridized carbon atoms. Defect densities are systematically varied from 0% to 20% in increments of 5% for each carbon type. This approach results in a total of 25 unique simulations that explore various possible combinations of defect densities. For these simulations, a square sheet of side length 215.95Å is used, totaling 21,721 atoms. This provides a large enough sample of atoms to examine how void defects affect the mechanical properties of penta-graphene on average. By analyzing these scenarios, the study aims to assess the relative importance of defects in the ${\mathrm{sp}}^2$ and ${\mathrm{sp}}^3$ regions and provide insights into how they affect the mechanical resilience of penta-graphene.

Figure 2 presents select stress-strain curves from void defect simulations of the penta-graphene sheet, demonstrating the impact of varying ${\mathrm{sp}}^2$ and ${\mathrm{sp}}^3$ carbon defect concentrations. The results reveal that as the void defect density increases, the ultimate stress before fracture decreases. However, beyond a certain defect threshold, the material exhibits greater strain before failure, indicating enhanced ductility. When the proportion of ${\mathrm{sp}}^3$ defects exceeds this critical threshold, the material becomes even more flexible, delaying fracture to higher strain values while simultaneously reducing its maximum load-bearing capacity.

This behavior can be attributed to the distinct mechanical roles of ${\mathrm{sp}}^2$ and ${\mathrm{sp}}^3$ hybridized carbon atoms in penta-graphene. Initially, as void defects increase, the removal of bonds reduces the structural integrity of the material, leading to a decrease in ultimate stress. This is because fewer intact bonds are available to sustain the applied load, resulting in earlier failure. However, beyond a critical defect density, localized flexibility emerges due to disruptions in the ${\mathrm{sp}}^3$ hybridized network, allowing for increased deformation before fracture. The presence of voids in ${\mathrm{sp}}^3$ regions weakens the material’s out-of-plane rigidity, making it more adaptable under strain. Consequently, the failure point shifts to higher strain values, but the weakened atomic framework results in a reduced maximum stress before failure. This balance between the rigidity provided by ${\mathrm{sp}}^3$ bonding and the planar strength maintained by ${\mathrm{sp}}^2$ bonding governs penta-graphene’s mechanical response to void defects.

The Young’s modulus and maximum strength of the penta-graphene sheets as a function of defect percentage are shown in Figure 3. The graph demonstrates that the elastic modulus steadily declines with increasing void defects. This figure also indicates that void defects in ${\mathrm{sp}}^3$-hybridized carbon atoms have a more pronounced impact on reducing the material’s overall elasticity compared to defects in ${\mathrm{sp}}^2$-hybridized regions. An increase in defects from a pristine sheet to 20% ${\mathrm{sp}}^2$ voids leads to a drop in elasticity of roughly 60 N/m, while a similar increase in ${\mathrm{sp}}^3$ void defects reduces the elasticity 70 N/m. If 20 % ${\mathrm{sp}}^3$ defects are already present, ${\mathrm{sp}}^2$ defects have a negligible effect on the overall stiffness of the material. This suggests that the integrity of the central plane of ${\mathrm{sp}}^3$-bonded atoms play a critical role in maintaining the material’s stiffness. The exterior ${\mathrm{sp}}^2$ atoms also contribute to the stiffness of the structure, but require a higher defect percent to produce the same effect.

In contrast, the fracture strength data reveal that the difference in the impact of increasing ${\mathrm{sp}}^2$ versus ${\mathrm{sp}}^3$ defects is less pronounced. A reduction in the concentration of either carbon type is equally likely to disrupt these critical ${\mathrm{sp}}^2$–${\mathrm{sp}}^3$ bonds, leading to structural failure. This underscores the importance of the ${\mathrm{sp}}^2$–${\mathrm{sp}}^3$ bonds as key contributors to the mechanical strength of penta-graphene. Maintaining the integrity of these mixed-hybridization bonds is crucial for preserving the material’s load-bearing capacity and resistance to fracture.

4.2. Impact of Initial Crack Size on the Mechanical Properties

Understanding the impact of pre-existing cracks on the mechanical properties of penta-graphene is crucial for evaluating its structural reliability and suitability for practical applications. Initial cracks can significantly influence the material’s mechanical strength and fracture behavior, making it essential to study their effects for accurate failure predictions. Gaining insights into crack-induced weaknesses is key to enhancing the durability and mechanical resilience of penta-graphene in advanced technologies such as flexible electronics, nanoelectromechanical systems (NEMS), and high-strength composite materials.

To examine the relationship between crack size and mechanical performance, a series of simulations were conducted with varying crack widths. The cracks introduced into the pristine material range from 30.59Å to 73.78Å in width, increasing in increments of 7.2Å. This increment, equivalent to two unit cells, was chosen to maintain a consistent microstructural configuration of the crack. The crack opening is set at 7.2Å to minimize artifacts arising from interactions between free crack edges. Figure 4 shows a penta-graphene sheet with an initial crack with a crack opening of 7.2Å. In simulations of crack defects, we use a simulation domain of 431.91 Å by 215.95 Å, including 43,123 atoms. The longer dimension is aligned with the direction of the crack edge to help mitigate edge effects. This systematic approach enables a comprehensive analysis of how initial crack sizes influence fracture behavior and mechanical stability.

Figure 5 presents the stress-strain curves for penta-graphene with different initial crack sizes, while Figure 6a and Figure 6b show the Young’s modulus and fracture strength respectively. As expected, an increase in initial crack size results in a notable reduction in both overall strength and stiffness. The stress-strain curves exhibit a characteristic linear elastic region, followed by failure at higher strain levels. This behavior is consistent with brittle fracture mechanics, where the material deforms elastically before experiencing catastrophic failure.

A key observation is the progressive decline in peak stress at failure and the overall stiffness as crack size increases, highlighting a weakened mechanical response. Larger cracks serve as stress concentrators, amplifying localized stress and expediting fracture initiation. The reduction in stiffness arises from the disruption of the material’s bonding network, which diminishes its load-bearing capacity.

Interestingly, some simulations exhibited a multi-stage failure process, where the material does not fail catastrophically in a single event. This phenomenon is particularly evident for larger crack sizes, where an initial drop in stress is followed by partial structural recovery, allowing the material to endure higher strains before complete failure. This behavior can be attributed to the atomic restructuring that occurs during the initial fracture.

Penta-graphene’s rigid three-dimensional structure contributes significantly to this mechanical response. When oriented in the X–Y plane, many atomic bonds retain components in the Z-direction. Under tensile strain, these out-of-plane components undergo atomic de-wrinkling, progressively flattening the structure. However, due to its pentagonal tiling, which does not allow perfect tessellation in a two-dimensional plane, penta-graphene experiences internal compression. This compression increases transverse stress on the bonds, ultimately leading to bond rupture. If the initial fracture does not lead to immediate catastrophic failure, the material can adapt by restructuring its atomic configuration, forming carbon loops that contain 6 to 12 atoms, as noted in [41]. This restructuring process is depicted in Figure 7, where areas of atomic reconfiguration are highlighted in green.

The stress distribution and crack growth at different strain levels for a penta-graphene sheet with an initial crack length of $52.18$ Å are shown in Figure 8. As the crack propagates, larger loops of carbon atoms form behind the crack tip, keeping the top and bottom surfaces connected. This phenomenon helps increase the failure strain and prevents sudden failure after the initial maximum stress is reached. The ability of these atomic chains to redistribute stress enhances the structural integrity of penta-graphene, contributing to its overall fracture resistance.

This phenomenon is influenced by both applied stress and temperature. The formation of the larger carbon loops helps relieve localized stress, enabling the material to sustain higher strains after the initial fracture. Due to stress concentration around crack tips, these phase transitions frequently originate at the crack corners and propagate diagonally away from the crack site, redistributing stress throughout the structure. Such structural adaptability contributes to the material’s unique fracture behavior, making penta-graphene an intriguing candidate for applications requiring both high strength and flexibility.

4.3. Effect of Loading Direction on the Behavior of Penta-Graphene

The simulations of previous sections examined the mechanical behavior of penta–graphene sheets constructed using the unit cell depicted in Figure 1, under tensile loading along the ${\mathbf{a}}_{\mathbf{2}}$ direction. To understand the effects of loading direction on the mechanical properties of penta–graphene, we expanded our study to include uniaxial tensile tests on penta–graphene sheets oriented at various angles.

In this series of simulations, penta–graphene sheets were generated with orientations varying from 0° to 90°, using the unit cells illustrated in Figure 9. The dimensions of each unit cell for the respective orientations are listed in Table 2. To ensure uniformity across tests, the sheets are designed with a periodic structure in all directions, with dimensions of approximately 215.95 by 215.95 Å, mirroring the configurations used in the void defect simulations.

The stress-strain curves for different orientations, are displayed in Figure 10. Due to the rotational symmetry of the penta-graphene unit cell, the material exhibits similar mechanical properties at both 0° and 90°. Figure 11 a and Figure 11 b further delineate how the loading angle impacts the Young’s modulus and fracture strength of penta-graphene.

Analysis of the results reveals that the fracture strength is relatively invariant to changes in the loading orientation, yet the modulus of elasticity, derived from the linear portion of the stress-strain curves, demonstrates a noticeable dependence on the loading angle. The lowest elastic modulus occurs when the penta-graphene sheet is oriented at 45°, highlighting a directional anisotropy in mechanical stiffness. This trend is symmetric about the 45° angle, likely reflecting the 90° rotational symmetry characteristic of penta–graphene. In this orientation, the material shows a decreased resistance to deformation, necessitating a strain of about 0.02 before substantial load-bearing capacity is achieved, resulting in a shift of the entire stress-strain curve toward higher strain values. Contrary to the elasticity, the fracture strength does not display a consistent trend with the angle of rotation, remaining relatively stable around 30 N/m with a variation of ±5

Fracture initiation is consistently observed at the bonds between ${\mathrm{sp}}^2$ and ${\mathrm{sp}}^3$ carbon atoms. Initially, the rupture of a single bond combines two pentagonal (${\mathrm{C}}_5$) loops into an octagonal (${\mathrm{C}}_8$) loop, precipitating a brittle failure. This phenomenon is particularly evident in sheets oriented at 0°, as illustrated in Figure 7. The formation of ${\mathrm{C}}_8$ loops redistributes stress to adjacent bonds, leading to a sequential rupture of bonds and the formation of additional ${\mathrm{C}}_8$ loops, thereby exacerbating the propagation of the crack.

5. Surface Energy of Penta-Graphene

Surface energy ($\gamma$) in penta-graphene quantifies the energy required to create new surfaces within the material. This is calculated by comparing the total energy of a pristine sheet ($U_0$) with that of the same sheet after introducing free surfaces (U), using the equation:

$$\gamma ={\displaystyle \frac{U-U_0}{2lt}}$$

(1)

where l is the length of the sheet along the considered axis, and t is the thickness of the sheet. The factor of 2 accounts for the creation of two new surfaces.

In this study, we assume a penta-graphene thickness of 3.68 Å, as estimated from multilayer simulations with periodic boundary conditions in the thickness direction [52]. To determine $U_0$, the penta-graphene sheet is simulated with periodic boundary conditions in both in-plane directions (X and Y), while keeping the Z-direction fixed. The total energy U is then obtained by modifying the boundary conditions to introduce a free surface along the desired fracture direction. Specifically, for surfaces oriented along the ${\mathbf{a}}_{\mathbf{1}}$-direction (normal to the Y-direction), periodicity is maintained in the Y direction, with the X and Z directions fixed. Similarly, for surfaces oriented along the ${\mathbf{a}}_{\mathbf{2}}$-direction (normal to the X-direction), periodicity is maintained in the X direction, with the Y and Z directions fixed. By substituting these computed energy values into Equation 1, the surface energies of penta-graphene along the ${\mathbf{a}}_{\mathbf{1}}$ and ${\mathbf{a}}_{\mathbf{2}}$ directions are found to be 4.434 J/m² and 6.771 J/m², respectively. These values are lower than the surface energy of graphene, which is about 8 J/m² [53].

While several studies have investigated the critical stress intensity factors of two-dimensional (2D) materials using MD simulations [50,54,55], penta-graphene presents unique fracture behavior due to its distinct atomic configuration.

6. Molecular Dynamics Simulations for Mixed-Mode Fracture in Penta-Graphene

The critical stress intensity factor is a fundamental parameter in fracture mechanics that quantifies a material’s resistance to crack propagation under applied stress [56]. For emerging nanomaterials like penta-graphene, understanding the critical stress intensity factor is crucial for evaluating their mechanical reliability and fracture toughness, particularly in the presence of defects such as pre-existing cracks or voids.

Studying the critical stress intensity factor of penta-graphene offers insights into how the material manages localized stress, particularly near defects. An understanding of critical stress intensity factor is crucial for predicting whether a crack will remain stable or grow uncontrollably, which is vital for preventing catastrophic failure in real-world applications. Additionally, exploring how the critical stress intensity factor is affected by factors such as crack geometry and loading conditions can guide the development of more defect-tolerant penta-graphene based materials.

This section presents MD simulations aimed at analyzing the critical stress intensity factor and crack propagation path of penta-graphene under mixed-mode loading. A circular domain surrounding the crack tip is employed for the analysis. The MD simulations apply asymptotic crack-tip displacement fields to the outermost atomic layers, as illustrated in Figure 12.

The far-field behavior is assumed to be linear and isotropic due to the rapid decay of the crack-tip strain field, whereas the anisotropic nature of penta-graphene becomes significant only under large strain conditions. The asymptotic displacement fields at the crack tip for a linear isotropic material under mixed-mode I and II fracture are given by equations 2 and 3 [50,57]:

$$\begin{array}{cc}\hfill U_x=& {\displaystyle \frac{1+\nu }{E}}\sqrt{{\displaystyle \frac{r}{2\pi }}}[K_{I}cos{\displaystyle \frac{\theta }{2}}(\kappa -1+2{\sin }^2{\displaystyle \frac{\theta }{2}})+\hfill \\ & K_{II}\sin {\displaystyle \frac{\theta }{2}}(\kappa +1+2{\cos }^2{\displaystyle \frac{\theta }{2}})]\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill U_y=& {\displaystyle \frac{1+\nu }{E}}\sqrt{{\displaystyle \frac{r}{2\pi }}}[K_{I}sin{\displaystyle \frac{\theta }{2}}(\kappa +1-2{\sin }^2{\displaystyle \frac{\theta }{2}})-\hfill \\ & K_{II}\sin {\displaystyle \frac{\theta }{2}}(\kappa -1-2{\sin }^2{\displaystyle \frac{\theta }{2}})]\hfill \end{array}$$

(3)

Here, $U_x$ and $U_y$ represent displacements in the x- and y-directions, respectively. The parameter E denotes the elastic modulus, and $\nu$ is Poisson’s ratio. The variable r is the radius of the circular domain, while $\theta$ is the angle in polar coordinates. The Kolosov constant $\kappa$ is given by ${\displaystyle \frac{3-\nu }{1+\nu }}$ for plane stress. The terms $K_{\mathrm{I}}$ and $K_{\mathrm{II}}$ denote the mode I (tensile) and mode II (in-plane shear) stress intensity factors, respectively.

For mixed-mode loading, the effective stress intensity factor is determined as

$$K_{\mathrm{eff}}=\sqrt{K_{\mathrm{I}}^2+K_{\mathrm{II}}^2}.$$

(4)

The loading phase angle $\varphi$ is given by

$$\varphi ={tan}^{-1}\left({\displaystyle \frac{K_{\mathrm{I}}}{K_{\mathrm{II}}}}\right),$$

(5)

where $\varphi =0^{\circ }$ corresponds to pure mode II (shear) loading and $\varphi ={90}^{\circ }$ corresponds to pure mode I (tensile) loading. Based on our MD simulations, the Young’s modulus for penta-graphene is taken as $E=216$ GPa, and Poisson’s ratio is $\nu =-0.174$ [45]. To examine the effect of crack edge structure and crack-tip configuration on the critical stress intensity factor, two distinct crack configurations (configurations A and B) are analyzed, as illustrated in Figure 13. Both configurations share identical lengths and widths to ensure a valid comparison.

The boundary conditions are applied incrementally with $\Delta K_{\mathrm{eff}}=0.01\mathrm{MPa}·\sqrt{\mathrm{m}}$. Time integration is performed using the velocity-Verlet algorithm with a time step of 1 femtosecond (fs). Simulations are conducted at a temperature of 300 K. Following each loading increment, the boundary atoms are fixed while the internal atoms are relaxed for 60,000 time steps, corresponding to a strain rate of $5\times {10}^{-5}{\mathrm{ps}}^{-1}$.

As MD simulation results can be sensitive to loading rates and due to the lack of experimental data on the fracture properties of penta-graphene, this study primarily provides qualitative insights into fracture behavior and trends. However, additional MD simulations with an extended relaxation period of 500,000 time steps were conducted, showing no qualitative differences in the reported results.

Due to the asymmetry at the crack tip and along the crack edges, simulations are performed for phase angles ranging from $-{90}^{\circ }$ to ${90}^{\circ }$. The boundary conditions are periodic in the X, Y, and Z directions, and the Tersoff potential is used for the simulations, which are conducted at 500K.

The crack propagation paths for configurations A and B are presented in Figure 14. These figures demonstrate that, despite differences in crack tip geometry, the cracks in both configurations propagate in the same direction when subjected to identical loading phase angles. This indicates that the loading phase angle is a critical factor in determining the crack propagation direction. At a ${90}^{\circ }$ loading phase angle (pure mode I), the crack propagation in a self similar path while for all the other phase angles, the crack propagates by kinking. For positive loading phase angles, the crack propagates downward, whereas, for negative loading phase angles, it propagates upward.

A notable observation is that under mixed-mode loading, the penta-graphene sheet experiences out-of-plane deformation. The out-of-plane deformations of the sheet under pure mode I and mode II loading are illustrated in Figure 16. Similar deformations have also been reported in other two-dimensional materials [58,59]. These deformations contribute to the absorption of external energy, delaying crack propagation and thereby increasing the critical stress intensity factor.

Figure 15 illustrates the variation of the critical stress intensity factor $K_{\mathrm{cr}}$ as a function of the loading phase angle for configurations A and B. Both configurations exhibit nearly identical values of $K_{\mathrm{cr}}$ at specific loading phase angles, suggesting that differences in the crack tip or crack edge have no significant influence on the effective stress intensity factor. Furthermore, the graph demonstrates symmetry for both positive and negative phase angles.

Figure 15. Effective critical stress intensity as a function of loading phase angle (configuration A and B)

Figure 15. Effective critical stress intensity as a function of loading phase angle (configuration A and B)

Figure 16. The out-of-plane deformation of penta-graphene sheet a) pure mode I, b) pure mode II

Figure 16. The out-of-plane deformation of penta-graphene sheet a) pure mode I, b) pure mode II

From these plots, the critical stress intensity factor for pure mode I is determined to be $17\mathrm{MPa}·\sqrt{\mathrm{m}}$, while for pure mode II, it is $15\mathrm{MPa}·\sqrt{\mathrm{m}}$. The stress intensity factor for penta-graphene is higher than that of graphene. Experimental measurements indicate that the mode I SIF for graphene is approximately $4\mathrm{MPa}·\sqrt{\mathrm{m}}$[46,53], while molecular dynamics simulations report values ranging from $2.63\mathrm{MPa}·\sqrt{\mathrm{m}}$ to $4.21\mathrm{MPa}·\sqrt{\mathrm{m}}$ [55,60,61].

The elevated stress intensity factor (SIF) of penta-graphene, relative to graphene, can be primarily attributed to its unique atomic architecture and associated deformation mechanisms. While graphene features a flat hexagonal lattice composed exclusively of sp²-hybridized carbon atoms, penta-graphene adopts a buckled configuration formed by pentagonal carbon rings with a combination of sp² and sp³ hybridization. This non-planar structure introduces the capacity for pronounced out-of-plane deformation, which is an intrinsic toughening mechanism. Despite possessing lower surface energy than graphene, penta-graphene’s ability to dissipate mechanical energy through structural flexibility enhances its resistance to crack initiation and propagation, resulting in a higher critical SIF.

The role of out-of-plane deformation is particularly significant in fracture processes. Under applied stress, penta-graphene can bend and locally rearrange its atomic structure, facilitating stress redistribution around crack tips. This energy is partially absorbed through mechanisms such as lattice distortion and in-plane stretching, thereby reducing the effective energy available for crack advancement. As a result, crack propagation is delayed, and the material exhibits greater toughness and a higher threshold for fracture.

This fracture behavior aligns with observations in other two-dimensional materials, such as ${\mathrm{MoTe}}_2$ and ${\mathrm{MoS}}_2$, where out-of-plane deformation has been shown to enhance toughness [54,62,63,64]. In our simulations, penta-graphene displayed similar deformation-driven toughening, supporting earlier findings on strain-induced energy dissipation mechanisms in 2D systems [65].

6.1. Stress Distribution

To analyze the stress distribution in a penta-graphene sheet, the stress tensor per atom is determined using the virial formulation. This method is commonly used in atomistic simulations to compute local stress fields by incorporating both kinetic and potential contributions to stress at the atomic scale. For an atom a, the virial stress is computed using the following expression [51,66,67,68].

$${{\sigma }^a}_{ij}=-{\displaystyle \frac{1}{V^a}}\left[{\displaystyle \frac{1}{2}}{m}^a{v^a}_i{v^a}_j+\sum _{b=1}^n{r}_j^{ab}{f^{ab}}_i\right]$$

(6)

In this equation, ${{\sigma }^a}_{ij}$ represents the stress tensor for atom a. The term $V^a$ denotes the atomic volume, while $m^a$ represents the atomic mass. The velocity components of atom $\mathbf{a}$ in the $i^{th}$ and $j^{th}$ directions are given by ${v_i}^a$ and ${v_j}^a$, respectively. The parameter n corresponds to the number of neighboring atoms within the interaction domain of atom a. The distance vector component between atoms a and b along the $j^{th}$ direction is denoted by ${r_j}^{ab}$, while ${f_i}^{ab}$ represents the $i^{th}$ component of the interatomic force exerted between these two atoms. This formulation accounts for both the kinetic energy contribution (first term) and the potential energy contribution (second term), allowing for a detailed assessment of atomic-level stress distribution in penta-graphene under various loading conditions.

Figure 17 illustrates the stress distribution ${\sigma }_{yy}$ for a loading phase angle of ${30}^{\circ }$ for both configurations A and B. The system temperature is maintained at 500K. For Both the configurations, the $K_{\mathrm{eff}}\mathrm{is}171\mathrm{MPa}\sqrt{\mathrm{m}}$.

The results reveal that stress concentration is highest near the crack tip in both configurations, where intense local deformation occurs. Given the mixed-mode loading condition at the loading phase of ${30}^{\circ }$, the crack kinks to propagates diagonally, exhibiting a combined shear and tensile fracture mechanism. An important observation from the stress distribution analysis is that despite differences in crack edge morphology between configurations A and B, both exhibit nearly identical stress distribution patterns around the crack tip. This suggests that, for penta-graphene, the intrinsic atomic structure dominates the stress response, leading to a self-similar stress field irrespective of crack edge variations. This behavior highlights the robust mechanical properties of penta-graphene, where stress localization and crack evolution are primarily dictated by the loading conditions rather than the specific shape of the crack edge.

7. Conclusion

In this study, we conducted a molecular dynamics investigation into the mechanical and fracture behavior of penta-graphene, a novel two-dimensional carbon allotrope characterized by its non-planar pentagonal structure and mixed ${\mathrm{sp}}^2$–${\mathrm{sp}}^3$ hybridization. Our simulations revealed that penta-graphene exhibits significant mechanical resilience, even in the presence of defects such as atomic vacancies and initial cracks. The material demonstrates a distinct capacity for stress redistribution through structural transformations at the crack tip, where the original pentagonal rings reconfigure into larger stable loops, such as hexagons and octagons. This adaptive bond rearrangement serves as a strain-energy dissipation mechanism that delays fracture and contributes to penta-graphene’s robustness under mechanical loading.

Our analysis of void defects showed that disruptions in ${\mathrm{sp}}^3$-hybridized carbon atoms have a more pronounced effect on reducing stiffness compared to ${\mathrm{sp}}^2$ defects. Nonetheless, at higher defect concentrations, the material exhibits enhanced ductility, indicating a complex interplay between rigidity and flexibility dictated by its hybridized bonding network. Similarly, the presence of initial cracks was found to reduce both the fracture strength and modulus of elasticity, yet the formation of new carbon ring structures during crack propagation prolongs the failure process, highlighting penta-graphene’s damage tolerance.

Molecular dynamics simulations revealed that penta-graphene exhibits a high critical stress intensity factor, exceeding that of graphene. This superior fracture toughness is attributed to the material’s ability to undergo substantial out-of-plane deformation under stress. These deformations redistribute stress around the crack tip and allow for greater energy absorption before fracture occurs. Importantly, we found that the crack propagation path and critical SIF are largely governed by the applied loading conditions rather than the initial crack geometry, emphasizing the dominant role of the intrinsic atomic structure in determining mechanical response.

Together, these findings underscore the potential of penta-graphene as a mechanically robust material suitable for demanding applications in nanoelectromechanical systems, flexible electronics, and structural composites. Its combination of defect tolerance, fracture resistance, and directional anisotropy presents a promising foundation for the development of next-generation two-dimensional materials.