Molecular dynamics of C–S–H production in graphene oxide environment

2.1 Principles of molecular dynamics

When a simulated system contains more than one molecule or atom, the energy of the system includes the sum of the total kinetic energy and total potential energy of all the atoms. The total potential energy, U, is a function of the position of the individual atoms. Specifically, it can be subdivided into two main components: the non-bonded van der Waals interaction between atoms, U _vdw, and the intramolecular potential energy, U _int. Specifically,

U _vdw can be approximated as the sum of the van der Waals interactions between pairs of atoms, whereas U _int is the sum of the potential energies for each type of coordinate (bond length stretching energy and bond angle bending energy). According to classical mechanics, the force F on any atom i is the gradient of the potential energy, and the acceleration of the atom is

By integrating Newton’s equations of motion with respect to time, the velocity and position of atom i after time t can be obtained:

where r i → represents the position of atom i, v i → represents the velocity of atom i, the superscript 0 represents the initial moment, and t represents the moment. With t = δt representing a very short time interval, the position and velocity of the ith atom after the interval δt can be calculated. The acceleration at this point is then recalculated according to the new position, and the position and velocity at the next moment are then recalculated to obtain information on the position and velocity of all atoms in the system as a function of time.

2.2 Model building

The early generation of C–S–H was studied using reaction molecular dynamics with reference to the simulation of the silica oxide gel process [15]. During the polymerization reaction, the main components of silica oxide colloids, namely C₃S and C₂S, are hydrolyzed to Si(OH)₄, which then undergoes polymerization aging [16]. To simplify the reaction process, Si(OH)₄ molecules were added directly to the simulation. Considering the model size effect and computational limitations, the number of Si(OH)₄ molecules was determined to be 6³, i.e., 216; H₂O/Si = 3 based on the influence of H₂O on the polymerization rate [17]. The calcium:silica ratio of C–S–H produced after the hydration reaction was 1.75 on average, and combined with the active energy of C–S–H and its thermodynamic stability, the simulated Ca:Si ratio was determined to be 1.5 [18].

The simulations were performed in large-scale atomic/molecular massively parallel simulator (LAMMPS; Sandia National Laboratories) with a reactive force field (ReaxFF), and the trajectories of the atoms were calculated using the Verlet algorithm with a time step of 0.25 fs [19]. The parameters of the ReaxFF were derived from a fit to quantum calculations [20]. First, the model was simulated kinetically using a regular system synthesis (NVT) so that the initial system was chilled for 100 ps at 300 K. The system was then chilled for 500 ps at 300 K and 0 pressure in an isothermal isobaric system (NPT). After two chilling processes, the architecture was optimized, and a steady state condition of energy and volume was achieved [21]. A temperature of 2,000 K was used to keep the model reacting at a constant temperature under NVT synthesis, and the number of simulation steps was set to 2,000,000 with a total duration of 500 ps. Finally, the simulation was concluded by gradually cooling the gel-generating system to 300 K and chilling it for 100 ps under NPT system conditions at 0 pressure to obtain C–S–H.

The initial model of graphene was derived from the graphite crystal of the system, and a graphene sheet layer with a = 39.36 Å and b = 39.36 Å was obtained by replication in the x and y directions; a carboxyl group (–COOH) was added to the graphene surface to form a GO sheet layer [22]. Two types of GO sheet layer with 25 and 50% oxygen content were simulated and constructed to study the effects of different oxygen contents of GO on the cementitious materials. Slits with different spacing of 40 and 10 Å were constructed to study the effect of different doping amounts of GO. The same simulation procedure and system parameters were used to obtain early C–S–H gel crystal models by placing molecular reactant solutions that are consistent with the above model into different GO sheet layers and different spacing slits [23]. The model construction parameters are listed in Table 1, and the initial system model is shown in Figure 1.

Figure 1

Initial system model (C, O, H, Ca, and Si elements are shown in gray, red, white, green, and yellow, respectively). (a) Control, (b) GCOOH40-25, (c) GCOOH40-50, and (d) GCOOH10-25.

3.1 Structural evolution of C–S–H

First, the changes in the reactions involving oxygen and silicon atoms in the models were investigated to compare the differences in the early C–S–H crystal structures in the different system models. The number of bridging oxygen atoms (BO, i.e., one oxygen atom linked to two silicon atoms) in the structure and those of non-bridging oxygen atoms (NBO, i.e., oxygen atoms linked to only one silicon atom) can be calculated by determining the number of partial coordination sites per atom [24]. The overall degree of polymerization in silicon oxide chain-structured colloids can be represented by the ratio of bridging oxygen atoms to silicon atoms (BO/Si), and the reaction rate of polymerization can be obtained from the BO/Si curve with time [25]. As shown in Figure 3, the ratio of BO to Si increased with an increase in the reaction time, reflecting a continuous increase in the overall polymerization of the system, which gradually stabilized in the later stages of the aging reaction [26]. A comparison of the variation curves of the groups indicates that the GO-free reference group started the reaction with a BO number of 0. The reaction process was stable and equilibrium was achieved in a very short time. In the model with the GO slit, a small amount of BO was formed at the end of the chilling process and the BO/Si ratio was approximately 0.1 at the beginning of the reaction, which also indicates that the C–S–H polymerization barrier was lowered and the chain structure was more easily formed under the effect of the site resistance of the nanoslit. Figure 4 represents the final content of BO at the bridge position. The GCOOH40-25 and GCOOH40-50 systems formed the highest amount of BO, with a BO/Si ratio of approximately 0.6. In contrast, GCOOH10-25 formed the least amount of BO and the BO/Si ratio was only 0.28, which is the lowest of all the systems. The initial slope of the BO/Si content curve corresponds to the polymerization rate of the aging reaction, which is closer to the GCOOH40-25 and GCOOH40-50 reaction rates in the early comparison group. The modulating effect of the nanoslits in the middle of the reaction allows the polymerization reaction to continue, with higher amounts of BO formed.

Figure 3

Variation curves of Bo/Si for the different systems.

Figure 4

Final bridging oxygen atom content.

The silicon atoms were classified into Q ₀–Q ₄ by counting the number of silicon atoms bonded to BO. The change curve of Q _n in each system is shown in Figure 5, whereas the final Q _n content distribution is shown in Figure 6. It can be observed that the reference group only contained Q ₀ at the beginning of the reaction, whereas Q ₁, Q ₂, and a small amount of Q ₃ were gradually formed in subsequent reactions. In the GO slit system, a small amount of Q ₁ was already present at the beginning of the reaction, as in the case of BO. The number of Q ₁ units in the GCOOH40-25 and GCOOH40-50 systems increased rapidly as the reaction proceeded and decreased after reaching peak values, whereas the number of Q ₂ and Q ₃ units increased, indicating that the silica atoms are more likely to produce highly polymerized unit structures in the GO environment [27]. Ultimately, Q ₂ and Q ₃ units stabilized at 19.5 and 6.5% and 23.6 and 8.8%, respectively, which were much higher than those of the reference group. Furthermore, some Q ₄ units were produced in both the systems. In the GCOOH10-25 system, the number of Q ₁–Q ₃ units was the lowest of all systems, and that of Q ₄ units was almost zero, as in the reference group.

Figure 5

Variation curves of Q _n content. (a) Q ₁, (b) Q ₂, (c) Q ₃, and (d) Q ₄.

Figure 6

Distribution of Q _n content in each system.

A comprehensive comparison of the structural evolution of the early C–S–H crystal-formation process in each system indicates that the systems in the GO slit environment formed a certain number of chain structures in the chirality phase before the start of the reaction because of the site-blocking effect of the GO lamellae and the influence of the reactive oxygen-containing groups. During the reaction, the reaction rate in the system did not change significantly when the slit length was 40 Å. However, the number of silicon atoms in the Q ₁ unit increased rapidly early in the reaction process and then decreased. The total number was essentially the same as the reference group when the reaction reached equilibrium, whereas the number of BO and number of silicon atoms in the Q ₂ and Q ₃ units increased significantly; a significant portion of Q ₄ units was produced, indicating that the molecular structure tended to converge in the GO environment with a slit of 40 Å, and the overall polymerization of the system increased effectively. A comparison of the content change curves of GCOOH40-25 and GCOOH40-50 indicates that under the same slit spacing conditions, the effect of 50% oxygenated GO lamellae on the system was slightly higher in the high polymerization units, such as BO quantity, Q ₃, and Q ₄ than in the 25% oxygenated GO system, although the difference between the two is not significant. This indicates that in the simulation, GO oxygenation affects the structure of C–S–H in the system evolution, and the chain structure is more easily polymerized in the GO environment with high oxygen content [28].

When the slit distance was set to 10 Å, a comparison of the model curves shows that although the chain structure was generated during the initial stages, the system reached equilibrium in a very short time as the reaction proceeded; the BO content and Q ₂ and Q ₃ units were all significantly lower than in the simulated system at 40 Å, and lower than in the reference group without GO. The proportion of Q ₀ units in the final Q _n unit distribution was more than 55%, i.e., more silicon atoms did not polymerize effectively, and the proportion of Q ₁ units with low polymerization was much higher than the other units in the Q _1–4 units, indicating that even if polymerization did occur, it was only a short-range chain structure of a small size, with the lowest overall polymerization performance. This suggests that excessive GO doping and a nanoslit that is too small will restrict the reaction space in the solution and affect the full reaction of the molecules, preventing the C–S–H from early polymerizing effectively.

3.2 Analysis of the RDF

The RDF is a function representing the change in the density of a target particle within a certain distance from a reference atom. This function can reflect the change in the odds (g(r)) of a particular particle appearing within a certain range of any atom in the model and the arrangement, and is a common method of analysis in molecular dynamics [29]. The RDFs obtained through X-ray small-angle diffraction tests are mutually Fourier-dependent and can be analyzed against simulation and test results [30].

By analyzing the final C–S–H model obtained for each group, the RDFs for the four types of atomic pairs, namely, Si–O, Si–Si, Ca–O, and Ca–Ca, can be calculated, as shown in Figure 7. It can be observed from Figure 7a that the Si–O bond function curves for all the models show a peak point at 1.6 Å, but the peak intensities are significantly different. The peak of the 40 Å slit model shows a significant increase compared with the control group, whereas that for the GCOOH40-50 model is slightly higher than that for GCOOH40-25. The GCOOH10-25 peak point is the lowest, and the differences in the other groups are more pronounced, indicating that the number of Si and O atoms bonded in the model increased and the density increased significantly in the GO environment with a slit of 40 Å. However, the Si–O bonding decreased with a slit of 10 Å. The steamed bun peak at 4 Å of the function curve corresponds to the non-bonded Si and O atoms in the model system. The peak magnitudes did not differ significantly between the models, with the GCOOH10-25 model, which has a lower aggregation ratio, exhibiting a relatively small peak. Figure 7b shows the RDF of Si–Si. The Si atoms are linked to each other by bridge-site oxygen atoms, and the main peak of the curve occurs at 3.2 Å. The peak intensities of the GCOOH40-50 and GCOOH40-25 models are significantly higher than those of the other two curves, whereas the peak of GCOOH10-25 is closer to that of the reference group. This indicates that the polymerization of Si atoms is effectively enhanced by the increase in the oxygen atoms in the bridge position under the effect of nano-GO. The Si–Si function curve shows a bun peak at 5 Å. The peaks of the curves are still relatively close to each other, and the values for the GCOOH10-25 group are, in contrast, slightly higher. It can be inferred from the Q _n distribution that this may be related to the uneven dispersion of Si atoms when the slit is too small, and the lack of bonding of the locally concentrated Si atoms causes the Si–Si function curve to yield a higher intensity at the bun peak. Figure 7c shows the RDF of Ca–O. The high degree of polymerization of the GCOOH40-50 and GCOOH40-25 models results in a higher peak at 2.4 Å compared with the other two groups, with a significant increase in the Ca–O bond density. It can be observed from the Ca–Ca RDF curve in Figure 7d that the GCOOH10-25 model has a lower intensity in the Ca–O function curve; however, the intensity of the main Ca–Ca peak is closer to those of the other models, indicating that the effect of GO still resulted in significant aggregation of Ca atoms in the solution even at low polymerization levels.

Figure 7

Atomic radial distribution function. (a) Si–O, (b) Si–Si, (c) Ca–O, and (d) Ca–Ca.

3.3 Distribution of chemical components

The distribution of Ca and Si atoms in the model in the z-axis was determined as shown in Figure 8, allowing the differences in the C–S–H density generated in the model system to be analyzed. A comparison of the distribution curves of Ca and Si in the reference group is shown in Figure 8a; it indicates that the atoms are more evenly distributed in the z-axis, and the Ca/Si ratio at each position remains essentially the same (1.5) in the initial system, with a peak in the number of atoms distributed in the middle of 25–30 Å range. This region is the C–S–H enriched region, which indicates that the C–S–H generated in the system tends to aggregate and form local clusters. Because there is a GO lamella of approximately 3 Å in the GCOOH40-25 model, the peaks of the distribution curves of Ca and Si in the 4–37 Å z-axis range in Figure 8b show a clear shift toward one side of the GO lamella, with the enriched region of Ca atom distribution appearing in the z-axis range of 7–12 Å, and the enriched region of Si atoms in the 10–13 Å range. This indicates that GO has an adsorption effect on both Ca and Si atoms. The adsorption of Ca atoms is more obvious. The distribution curves of Ca and Si for GCOOH40-50 in Figure 8c are more similar to the GCOOH40-25 model, whereas the peak number of Ca atoms further increases, creating a region of higher Ca/Si ratio locally. It can be concluded that the oxygen atoms in GO can effectively adsorb Ca and Si atomic groups, and the higher the amount of oxygen contained, the better the adsorption effect. Simulations of C–S–H gel polymerization with different Ca/Si ratios show that the C–S–H production potential varies with different Ca/Si ratios, with the enthalpy of C–S–H production decreasing as this ratio increases, whereas the activation energy decreases and then increases as the ratio increases, reaching a minimum at a ratio of 1.8 [31]. At this point, if the number of Ca atoms continues to increase, more NBOs are induced within the system, decreasing the overall polymerization of the colloid [32]. Therefore, the enrichment of suitable Ca atoms in the GO-containing model increases the local Ca/Si ratio of the system [33], at which Ca atoms catalyze the early gel polymerization reaction and act to isolate the nano-space, thereby allowing the enrichment of Si(OH)₄ monomers, increasing the chance of collision, and increasing the rate of polymerization [33]. It can be observed from the Ca and Si distribution curves of GCOOH10-25 in Figure 8d that when the GO slit was too small, there was a non-uniform distribution of atoms in the system, with Ca atoms enriched at 3.8 and 9.2 Å and Si atoms enriched at 5 and 8 Å. The Ca/Si ratio varied significantly in different regions, resulting in the lowest degree of C–S–H polymerization.

Figure 8

Distribution of Ca and Si atoms in z direction. (a) Control, (b) GCOOH40-25, (c) GCOOH40-50, and (d) GCOOH10-25.

3.4 Distribution of Q _n

The atomic structure of the gel largely determines the chemical composition of C–S–H; therefore, by exploring the distribution of Q _n units along the z-axis, it is possible to characterize the structural changes of C–S–H in the system. Figure 9a shows the distribution of Q _n units in the control group. It can be observed that the Q ₁ and Q ₂ units appear to aggregate in larger numbers around 25 Å, whereas the Q ₃ units are fewer and randomly distributed, indicating that C–S–H is mostly a short-chain structure with clear convergence locations. There are also more Q ₀ units in the system, indicating that the system contains C–H. This demonstrates that the early C–S–H system had a low degree of polymerization and the system was not homogeneous. Figure 9b shows the distribution of Q _n units in the GCOOH40-25 model, where the Q ₀ units are reduced compared with the reference group. The overall distribution of Q ₁ units is relatively homogeneous with a slight increase in number at 10 Å; more Q ₂ and Q ₃ units are aggregated at this position, which is consistent with the silica-rich region in Figure 8b, indicating that the highly linked C–S–H structure is likely to occur in the calcium-rich region with a high Ca/Si ratio. The presence of Q ₄ units in the system, however, indicates that GO enhanced early C–S–H aggregation. Figure 9c shows the distribution of Q _n units in the GCOOH40-50 model. Similar to the GCOOH40-25 model, the number of Q ₁ units is further reduced, and Q ₂, Q ₃, and Q ₄ units all exhibited content peaks at 8–10 Å. The system is more polymerized and has more highly linked structures. In the GCOOH10-25 model shown in Figure 9d, the Q _n unit content distribution appears quite different from the other systems and essentially maintains a high consistency with the Si atom distribution in Figure 8d, with Q ₀ and Q ₁ units obviously converging at 4.5 and 8 Å, respectively. The atomic uniformity of the system is poor, and the Q ₂ unit content is already very low, with no Q ₄ units generated.

Figure 9

Q _n element distribution in z direction. (a) Control, (b) GCOOH40-25, (c) GCOOH40-50, and (d) GCOOH10-25.

The above analysis shows that when the GO slit spacing is moderate, i.e., when the GO doping amount is ideal, the adsorption and limiting effect of GO lamellae can effectively improve the polymerization of C–S–H and facilitate early C–S–H production. However, when the GO slit spacing is too small, i.e., when GO addition is significantly increased, the polymerization of C–S–H shows a significant weakening, which is consistent with the results of macroscopic tests where adding too much nanomaterial does not contribute to further improvement in the properties of the cementitious materials [34]. In most cases, this is due to the agglomeration effect of the nanomaterials, similar to macro-fiber addition, where there is a small area of entanglement and aggregation in the area of excessive fiber density at high addition levels [35]. In this area, the transition zone between the cement paste and the fibers is dense, whereas the anchoring effect of the paste on the fibers is weakened and the bridging effect of the fibers is not present, such that this area within the matrix becomes weak and less durable. Therefore, the material admixture is also usually too high, resulting in difficult or uneven dispersion, which is the main reason for the poor modification effect [36].

In addition to the agglomeration effect described above, doping by GO has another effect based on the analysis of molecular reaction process simulations. When the slit is set to 10 Å, the Ca and Si distribution curves appear significantly different, with Ca atoms aggregating in response to GO; however, the Ca/Si ratios in different regions demonstrate strong randomness. This resulted in the system reaching equilibrium in a very short time, although the chain structure was generated earlier in the early stage of the aging process with a high content of monomeric Q ₁ units, whereas the BO, Q ₂, and Q ₃ units all had a low content, and essentially no Q ₄ units with high polymerization were formed. Moreover, the products in the system were mostly short-range chain structures of a small size, with extremely low overall polymerization performance. This indicates that when GO addition is high, the spacing between the nanolayers is too small and the difference in the adsorption of ions by GO during the hydration process becomes more obvious. This limiting effect also causes the ions within the lamellae to accumulate and not move sufficiently, resulting in a lower overall polymerization of C–S–H. This microscopic fragmentation effect tends to break up and unitize the gel, reducing the overall strength of the material. In other words, when the amount of nanomaterials added is too high, even if there is no agglomeration and the dispersion effect is extremely homogeneous, this fragmentation effect may become more obvious, resulting in poor modification of the composite material.