Coarse-Grained Monte Carlo Simulations of Graphene-Enhanced Geopolymer Nanocomposite Nucleation

1. Introduction

The word “geopolymer,” introduced by Davidovits in 1978 [1], denotes materials formed by combining aluminosilicate powder with a precursor solution of potassium silicate or sodium silicate [2,3]. The initiation of the alkalinization process is significantly influenced by silicate solutions precursors. In the dissolution-precipitation geopolymerization reaction, hydroxyl (OH^-) groups in the alkaline solution first break the bonds between silicate and aluminate in the solid aluminosilicate powder material (metakaolin, MK) [4,5]. This is followed by a polycondensation reaction that leads to the creation of an aluminosilicate network [6]. This network comprises interconnected aluminate and silicate (Si) tetrahedra connected by oxygen bridging bonds [3]. Geopolymer binders demonstrate mechanical properties comparable to those of ordinary Portland cement, yet their production generates roughly 80–90% fewer global anthropogenic CO₂ emissions. This characteristic establishes geopolymer as a sustainable and eco-friendly substitute for conventional cement [7]. The core chemical and structural properties of geopolymers produced from MK [8], slag, and fly ash are analyzed to evaluate how the selection of raw materials impacts the characteristics of the resulting geopolymer composites [3]. Geopolymers exhibit superior resistance to acid compared to Portland cement, a characteristic frequently highlighted as one of the key advantages of geopolymers [9,10,11]. Additionally, geopolymers possess inherent fire resistance and exhibit exceptional thermal stability, surpassing that of conventional cements [12].

Nanoscale modeling of geopolymers has recently attracted significant interest among researchers [13,14]. Lolli et al. [15] developed an innovative framework for understanding the molecular structure of geopolymers, highlighting the significance of nanoscale interfaces between crystalline and amorphous regions in influencing the material’s mechanical behavior. Izadifar et al. [16] studied a polymerization mechanisms of alkaline aluminosilicate gels using a three-dimensional off-lattice coarse-grained Monte Carlo (CGMC) method, which also enabled the characterization of their nanostructural features, such as particle size distributions and pore size distributions. The input data for this analysis included the Gibbs free energy values of dimerization reactions for four different monomer species, obtained from the literature and computed using the DFT modeling approach as described by White et al. [7,17]. Moreover, Valencia et al. [18] lately employed CGMC simulation method with octree cells for geopolymer nucleation [19] at different pH values. Yang and White [20] employed alkali-activated materials (AAMs), specifically class F fly ash and MK, at the mesoscale using the CGMC simulation method, combined with the DFT computational method. Their findings revealed that in both H-activated fly ash and MK systems, the growth of the gel occurred via the development of medium-sized clusters in conjunction with the expansion of the largest particle. Subsequent to this process, certain small clusters persisted in the pore solution of the solidified gel, with their dimensions determined by the saturation level of pre-dissolved silicate concentrations within the system. Furthermore, Izadifar et al. [4,5] have recently calculated the enthalpy activation energy (ΔH^∗) under far-from-equilibrium conditions, utilizing the transition state theory (TST) to determine atomistic reaction rates for silicate tetrahedra dissolution in MK via the (DFT) computational approach [21]. White et al. [22] also employed quantum chemical-based interaction energies, computed via density functional theory computations. They utilized an on-lattice CGMC simulation to investigate the primary stages of gel/cluster formation in sodium silicate systems across various concentrations. Their model implemented on a cubic lattice comprising 125,000 sites in the canonical ensemble (NVT) with periodic boundary conditions, facilitated the assessment of structure evolution resulting from polymerization reactions through Monte-Carlo moves and total system energy minimization.

Materials with two-dimensional structures, like graphene [23], have gained remarkable attention in nanomaterials study because of their exceptional thermal [24,25,26,27], mechanical [28,29,30,31], and electrical properties [32]. The GRGN has been explored across various advanced material applications, including structural health monitoring, structural supercapacitors, hydrogen production, 3D printing, dye wastewater treatment, and concrete, showing promising outcomes [33]. Zhang et al. [34] synthesized a new graphene-bottom-ash-based geopolymer (GBAG) composite using alkaline-activated geopolymerization. The primary aim was to enhance the electrical conductivity of GBAG by integrating graphene into its otherwise non-conductive matrix. Zhong et al. [35] for the first time reported the extrusion based 3D printing of geopolymer/GO nanocomposite structures. Their study showed that incorporating GO significantly alters the rheological behavior of the geopolymer precursor, facilitating the 3D printing of structures with the GO/geopolymer (GOGP) blend, a capability not achievable with pure geopolymer. Thus, the primary aim of the current study is to perform a 3D off-lattice CGMC simulation method, both with (GRGN) and without a pristine graphene nanosheet, to investigate the nucleation mechanism of alkaline aluminosilicate gel in a silicate-activated system. The Gibbs free energy associated with the process of dimerization, derived from a study by White et al. [7], is used for four monomer species of Si(OH)₄, Al(OH)₄⁻.Na⁺, SiO(OH)₃⁻.Na⁺.3H₂O, and SiO₂(OH)₂²⁻.2Na⁺.6H₂O. Furthermore, we incorporated the adsorption energy (E_ad), including van der Waals (vdW) interactions between a pristine graphene nanosheet and four monomer species, as detailed in Table 1 of the recent study by Izadifar et al. [36], to investigate the nucleation process in this study. Each monomer type is modeled as a separate coarse-grained particle type to investigate the evolution of gel structure across varying numbers of iterations. We explore pH value of 11, employing an expanded simulation system with Octree cell expansion to enhance accuracy. The system’s total energy is calculated at various iterations, offering valuable insights into its behavior. Moreover, we analyze the formation of clusters and the dissolution of MK over 56,000,000 iterations. The final structure is analyzed to evaluate cluster size distribution and pore network features.

2. Simulation Model and Method

2.1. Atomistic Model Preparation

The off-lattice CGMC method begins by defining the number of particles for the four different monomer species, based on the selected activated solution and pH level. For the chosen silicate-activated system, the contributions of water as well as Na (68.4%), silicate in solution (10.6%), and MK (21.0%) have been acquired from the data (shown in Table 3) provided by White et al. [7]. For this, a cubic simulation box with a size of 100 Å was selected, and the percentage of each type of monomer specie was computed. For the silicate-activated system, the contribution of silicate in the solution and MK are equivalent to the simulation boxes with dimensions of 47.32, and 59.43 Å, respectively. Hence, the total count of silicate monomers in solution and MK particles (in a simplified crystal arrangement) within the simulation box, with average diameters of 3.3 and 4 Å, are determined as 2744 and 2744, respectively. Given the three types of silicate monomers in solution, namely Si(OH)₄, SiO(OH)₃^-.Na⁺.3H₂O, SiO₂(OH)₂^2-.2Na⁺.6H₂O, the distribution percentages of these three contributed silicate monomer types (2744 particles in total) to maintain the pH at 11 are computed as 5%, 90%, and 5%, respectively [37]. These percentages translate to particle numbers of 137, 2470, and 137 respectively, as per Figure 5 published by Sefcik and McCormick [37]. Furthermore, the proportion of particles in MK present as alumina (aluminum oxide) and silica (silicon dioxide) are determined to be 42% (1152 particles) and 58% (1592 particles), respectively [7]. Moreover, the adsorption energy (E_ad) incorporating van der Waals (vdW) dispersion forces between pristine graphene nanosheet and four monomer species of Si(OH)₄, Al(OH)₄⁻.Na⁺, SiO(OH)₃⁻.Na⁺.3H₂O, and SiO₂(OH)₂²⁻.2Na⁺.6H₂O extracted form Table 1, reported by Izadifar et al. [36].

2.2. Monte Carlo Approach: Implementation in MATLAB Code

A custom MATLAB code [38,39] was created to manage and analyze data within exclusive cells following an Octree pattern, as depicted in Figure 1. To include data from adjacent cells, a neighbor relationship was defined, where two sites were considered neighbors only if they shared a common face; adjacency via an edge or corner not being sufficient. This methodology ensured that all interactions between clusters and monomers across cells were captured. It is also worth mentioning that all the reported numbers from section 2.1 concerning silicate monomers in solution, as well as silica and alumina particles from MK must be divided by 8 partitions, each represents as a subsystem according to Octree pattern. Figure 2 provides a detailed side view where two subsystems share particles involved in the geopolymerization process. The software tracked detailed information during the simulation, covering monomer specifications, cluster dynamics, system energy fluctuations, and the complex MK dissolution process. Post-simulation analysis involved using the global scan method to extract cluster size distribution. A refined pore network model was constructed, assuming idealized spherical pores, with monomers and dimers considered as aqueous species, as reported by White et al. [7]. Particle structure digitization utilized a watershed algorithm and city-block distance transform function (Figure 3). Subsequently, the analysis of pore connectivity and size distribution was conducted based on the approach detailed by Izadifar et al. [16]. The system’s overall energy was evaluated using the Gibbs free energy associated with dimerization, as documented by White et al. [7].

2.3. Octree Cells Approach: MATLAB Program Development

Models based on atomistic theory often face challenges in capturing mesoscale phenomena in zeolites due to the high computational costs involved in simulating across the necessary length scales [2]. The principal aim is to explore particle behavior during polymerization processes by employing a simplified mechanical model.. This investigation is essential for understanding nanoparticle formation, as larger system sizes are required to produce statistically meaningful results. Furthermore, simplified models enable a more concentrated examination of particular system characteristics [40]. The simulation methodology used in this study builds upon our former research [16] with a system size of 100x100x100 ų. To facilitate this, an octree cell structure with an initial partitioning level was developed, enabling parallel simulation, which is commonly utilized in High-Performance Computing (HPC) environments. Specifically, eight identical simulations were executed simultaneously, each with 7,000,000 iterations, and at the end, their results were combined for a thorough structural analysis. (as detailed in Section 2.2). The use of Octree patterns considerably lowered the memory demands for the CGMC solver when simulating large-scale systems, balancing the reduction of artificial effects in smaller lattice simulations with achieving convergence in a rational time frame [18].

2.4. Density Functional Theory (DFT) Computational Modeling Method

Density functional theory (DFT) [41] computations modeling method was conducted to determine the geometrical parameters of the silicate and aluminate monomers (refer to Figure 4). The Vienna ab initio simulation package (VASP) [42,43,44,45,46] was utilized, employing the projected-augmented wave (PAW) method [47] and pseudopotential to model electron-ion interaction. The electron exchange and correlation functional were selected based on the generalized gradient approximation (GGA) with the Perdew−Burke−Ernzerhof (PBE) parametrization [48]. The Brillouin zone was sampled using a well-converged k-sampling equivalent, given by 1 × 1 × 1 Monkhorst-Pack k-points mesh size for the entire system [49]. A plane-wave cutoff energy of 400 eV was employed for structural relaxations, ensuring convergence. The electronic self-consistent cycles convergence criterion was set at 10^-6 eV, while ion relaxation continued until forces fell below 10^-3 eV/Å. Additionally, structural analysis was carried out using the three-dimensional visualization software VESTA [50].

3. Results and Discussions

The distances between bonded atoms and the angles formed by these bonds, computed through DFT calculations, are detailed in Table A1-2 from our recent study, respectively [18]. Based on tetrahedral formation, we calculated β angles of 135.12° and 138.28° for the dimerization reactions of Si3-O11-Al1 and Si1-O4-Si2 (as shown in Figure 4), respectively. For Si-O and Al-O bonds in tetrahedral monomer formation, we considered mean values of bond lengths of 1.65 and 1.76 Å, respectively, which define the radius of coarse-grained particles specific to each particle type. Additionally, mean angles of 108° and 107° were assumed for O-Al-O and O-Si-O angles in formation of tetrahedral monomer, respectively, specifying each particle pairing type. The process of polycondensation brings the monomers together, leading them to connect at a single point—the shared center of the bonding oxygen between both particles (as illustrated in Figure 4). At pH 11, images capturing the evolution of the two geopolymer systems including without and with the contribution of graphene nanosheet at a certain number of iterations of 0, 40,000, 80,000, and 56,000,000 are observed from Figure 5A-B, respectively. The coarse-grained particles representing various monomer building units are color-coded as follows: Si(OH)₄ particles are depicted in cyan, SiO(OH)₃^-.Na⁺.3H₂O particles in blue, SiO2(OH)₂^2-.2Na⁺.6H₂O particles in green, and Al(OH)₄^-.Na⁺ particles in red. The final equilibrium condition was reached after 56 million iterations. The percentage of total particles involved in gel formation (cluster formation) was found to be 59.48% without the contribution of a graphene nanosheet and 57.27% with the contribution of a graphene nanosheet. Figure 6 depicts the equilibrium conditions observed during the energy computation for the silicate-activated system at pH 11, both with and without the contribution of a graphene nanosheet, after 8 million iterations within the solution, excluding MK.

Indeed, this pre-equilibrium condition was conducted for dissolved activator silicates. It is notable that the system containing pristine graphene nanosheet exhibited a lower energy of -819 kJ/mol compared to the system without the contribution of graphene nanosheet, indicating a more favorable condition. As depicted in Figure 7, the moment when MK is introduced into the system is marked as iteration 0. In other words, pre-equilibration condition obtained from iteration 8,000,000 (as shown in Figure 6) has been set to iteration 0 with the MK involvement as shown in Figure 7. Therefore, the system’s energy at the simulation’s outset is not zero, indicating the pre-equilibration process, involving the silicate monomers in the solution as an activator, is examined.

Figure 7 indicates a swift reduction in energy, implying the dissolution of MK particles and the initiation of the polymerization process, where aluminum monomers from the MK particles play a crucial role. After 500,000 iterations, the energy trend begins to stabilize, and by 56 million iterations, the energy values for the systems, both with and without the presence of the graphene nanosheet converge completely to values of -13,997 kJ/mol, and 14,228 kJ/mol, respectively. Figure 8 depicts the variation in the percentage of aluminate and silicate monomers for two systems, one with graphene contribution and one without, at pH 11 throughout the simulation. At the initiation of MK dissolution, the concentration of aluminate monomers remains at zero for both systems with and without the contribution of a graphene nanosheet, as all aluminate particles are initially contained within the MK particle. In contrast, the total concentration of silicate monomers present in both systems from the pre-equilibrium condition at 8 million iterations is 36.42% (without graphene nanosheet) and 36.80% (with graphene nanosheet). These starting values represent the quantity of silicate monomers within the solution. Over time, the concentrations of aluminate and silicate monomers rise significantly due to the dissolution of the precursors. The highest amount of silicate monomers in the solution are 51.87% at iteration of 16,000 (with graphene nanosheet), and 51.00% (without graphene nanosheet) at iteration of 15,200.

Moreover, the highest amount of aluminate presented as monomers in the solution reach to 14.37% (with graphene nanosheet), and 14.27% (without graphene nanosheet) at iteration of 17,600. After 56 million iterations, the number of monomers in the solution decreased due to the cluster formation. It was observed that the remaining aluminate and silicate monomers in the solution was fewer in the system without graphene nanosheet compared to the system with graphene nanosheet. To enhance comprehension of the varying silicate contributions to the cluster formation over time, Figure 9 illustrates the total number of silicate monomer species over different iterations. The three different types of monomer species, namely Si(OH)₄, SiO(OH)₃^-.Na⁺.3H₂O, and SiO₂(OH)₂^2-.2Na⁺.6H₂O contributed 1.80, 32.81, and 1.80% (total 36.41%, as depicted in Figure 8) at the beginning of the simulation for the system without graphene nanosheet, respectively. In the case of graphene nanosheet contribution, the same previously monomer species contributed 1.82, 32.48, and 1.82% (total 36.12%, as depicted in Figure 8), respectively. It is also worth reporting that the amount of SiO(OH)₃^-.Na⁺.3H₂O in the solution after 8 million pre-equilibrium condition and at the beginning of MK dissolution is high (32.81% for system without graphene, and 32.48% for the system including graphene nanosheet) due to its high contribution equals to 90% (as explained earlier in section 2.1) in the silicate-activated system.

The highest amount of these three silicate monomer presented at 15,200 iterations due to MK dissolution, and at the end after 56 million iterations, the amount of SiO(OH)₃^-.Na⁺.3H₂O decreased to 30.96%, and 32.26% for both systems without and with the incorporation of a graphene nanosheet, respectively. Figure 10 illustrates the trend of total particles present as monomer and cluster during 56 million iterations for both systems, with and without the incorporation of a graphene nanosheet. It is essential to note that graphene particles have not been factored into the reported percentage of monomers and clusters contribution for the system including graphene nanosheet, as the role of graphene is crucial for adsorption of different species onto it. The point where both monomers and clusters achieved an equal percentage (50%) occurred at 466,400 iterations for the system with the graphene nanosheet and at 398,400 iterations for the system without it. After 56 million iterations, it is evident that the number of particles participating in cluster formation is lower in the system with graphene nanosheet (57.27%) compared to the system without it (59.48%). Conversely, the energy of the system including graphene nanosheet is more favorable (more negative) (-14,228 kJ/mol) than the system without graphene (-13,997 kJ/mol), as depicted in Figure 7. This discrepancy can be attributed to the higher adsorption energy between the graphene nanosheet and particles contributing to cluster formation.

Figure 11 illustrates the progression of existing MK particles during the dissolution process for both systems, with and without the presence of graphene nanosheet. The total dissolution of MK at various iterations has been observed for both systems, with and without the inclusion of a graphene nanosheet. The main reason of various computed iterations for MK dissolution obtained from both systems, can be attributed to the fact that the MK dissolution can be accepted or rejected. After a random MK particle dissolution (either alumina or silica), the system checks for overlaps with other particles. If no overlap is detected, the particle movement is accepted within the Monte Carlo (MC) framework if the total system energy after the move is lower than the energy before the move. If the energy increases (becomes less negative), the movement may still be accepted based on the probability X, which is compared to a randomly selected value between 0 and 1. The probability X is determined using the Boltzmann factor for the change in configuration, as shown in equation 1, where kB represents the Boltzmann constant, T is the temperature, and ΔE is the energy change.

$$X=e^{{\displaystyle \frac{-\Delta E}{kBT}}}$$

(1)

$$\Delta E=E\left(afterparticlemovement\right)-E\left(beforeparticlemovement\right)$$

(2)

If the movement is not accepted, the particle will return to its original position (MK), and the loop for the next step will begin. In the next step, a particle is randomly selected, and it is moved in any direction to a new position while ensuring that the connections between particles maintain the tetrahedral structure. The displacement of the particle is determined by generating a random number between zero and one, which is then multiplied by the radius of the moving particle. It is also worth mentioning that particles participate into cluster formation, when the distance between the moved particle and its nearest neighbor is within 1 Å. The system’s total energy is then recalculated to determine whether it has increased or decreased as a result of this change as described above. After every 20 iterations, the pH of system is also checked to make sure to maintain the pH of the system at 11. If the pH is changed, the pH must be accordingly controlled by 5% (Si(OH)₄), 90% (SiO(OH)₃^-.Na⁺.3H₂O), and 5% (SiO₂(OH)₂^2-.2Na+.6H₂O) [37]. Figure 12 also plotted to illustrate the distribution of pore sizes for both systems. The pore size distribution in both systems exhibits a similar trend. At a high probability density of 56.5%, pore diameters of 1.62, and 1.57 nm were observed for the systems with and without the incorporation of a graphene nanosheet, respectively.

4. Conclusions

This research explored the nucleation of aluminosilicate gel using the off-lattice coarse-grained Monte Carlo (CGMC) method, both with (GRGN) and without graphene nanosheets, at a system pH of 11. The tetrahedral geometry and binding energy parameters of aluminate and silicate monomers were derived from DFT simulations. To efficiently scale up the CGMC system sizes, a method employing the Octree cell approach was developed, significantly reducing computational time and taking full advantage of parallel high-performance computing. The results revealed that the proportion of particles involved in cluster formation was lower in the system with graphene compared to the one without it. In contrast, the system containing graphene displayed a more favorable energy state, which can be ascribed to the stronger adsorption energy between the graphene nanosheet and the particles participating in cluster formation. Additionally, the complete dissolution of MK required more iterations in the system with graphene than in the system without it. Regarding pore size distribution, both systems showed a similar trend, with pore diameters of 1.62 nm and 1.57 nm at a high probability density of 56.5% for the systems with and without graphene, respectively.