1 Introduction↩︎

Epoxy resins are widely used in industries such as aerospace, automotive and microelectronics, owing to their superior thermal stability, mechanical strength, and favourable strength-to-weight ratio. The exposure of epoxy networks to humid or wet environments degrades their mechanical and thermal properties, including reductions in Young’s modulus, tensile strength, bulk modulus, and glass transition temperature, accompanied by plasticization, swelling, microcracking, and chemical degradation [1], [2]. Indeed, their relatively high moisture permeability, compared to metals or ceramics, prevents the total exclusion of water from epoxy resins during or after their synthesis. However, the mechanisms by which water molecules alter the material properties of epoxy resins remain poorly understood.

The interaction between water molecules and epoxy resins has been the focus of several studies based on molecular dynamics (MD) simulations. Early on, the equilibrium organisation and dynamics of absorbed water was characterised by means of radial distribution functions (RDF) and mean squared displacements (MSD) of water molecules inserted in cured epoxy resins free volume [3], [4]. The dynamics of water molecules were found to be heterogeneous hinting that a distinction should be made between water molecules directly bound to the epoxy resin via hydrogen bonds and water molecules surrounded by other water molecules. In complementary FTIR spectroscopy and MD studies [5], hydrogen bonding was also found to be a key metric to distinguish different types or dynamics of water molecules confined in epoxy resins, and helped to determine that water molecules tend to disperse in epoxies rather than forming aggregates. At higher water content, reaching up to \(12\)%, however, cluster analysis revealed that water molecules finally tend to aggregate, with cluster size increasing exponentially with the water content [6]. This same study also hinted that at low hydration (\(3\)%) an increase of Young’s modulus could be found, while the same modulus started to decrease at higher hydration (\(3\)%). An inverse trend has been found recently using reactive potentials, and low water hydration, approximately \(1.5\)%, causing a decrease of the Young’s modulus of \(40\)% but an increase of the bulk modulus [7]. Overall, although it has been the subject of several papers, the effect of hydration on the elastic properties has not fully been determined. In general terms, there is a lack of studies in most reported MD studies, frequently with large and overlapping error bars [8]. The case of MD simulations of complex systems such as hydrated epoxy resins is particularly striking, whereby the mechanical properties at different hydration levels cannot be differentiated [9].

Overall, this calls for new investigations using MD simulations, this time, with much larger ensembles. In order to enhance the moisture resistance of epoxy systems, the incorporation of nanofillers such as graphene has attracted considerable attention [10]. Graphene offers outstanding mechanical strength and interfacial interactions with polymer matrices, making it a promising reinforcement for epoxy composites [11]. Using large-scale large-ensemble MD simulations, we recently demonstrated the mechanisms enabling the reinforcement of polymers by graphene [12]. When it comes to hydration, graphene is considered for multiple reasons: (i) its barrier effects, reducing water uptake, and (ii) to mitigate the decrease of mechanical properties of epoxy resins with increasing hydration. Early experimental results obtained with graphene-oxide supported such hopes [13], nonetheless clear mechanistic explanations remained missing. This led to several studies of hydrated epoxy-graphene nanocomposites using MD simulations. These confirmed that in graphene-oxide systems, water diffusion was found to be lower due to its affinity with the reinforcement [14]. The mechanical properties of graphene-epoxy nanocomposites were investigated observing the degradation with increasing hydration of the interface between the matrix and the nanoparticle [15]–[17]. Moreover, these studies left bulk behaviour unexplored and focused on only two hydration levels.

In this work, we employ MD simulations to systematically investigate the effects of water content and graphene nanoparticle inclusion on the properties of water and thermomechanical properties of crosslinked epoxy networks. We focus on water structure, density, glass transition and elastic properties. We model the common and commercially available epoxy resin model epoxy diglycidyl ether of bisphenol F (DGEBF) cured with 4,4’-methylene-bis-(2,6-diethylaniline) (MDEA) and provide insight into the molecular mechanisms by which graphene reinforcement modifies the properties of hydrated epoxy nanocomposites.

To make actionable predictions from molecular dynamics (MD) simulations, we ensure control over the intrinsic uncertainty of our predictions. In MD simulations, the stress tensor values, and by extension, the elastic moduli, calculated from finite deformations have significantly high variance and act chaotically [18]. We have shown in our previous studies that ensemble MD, where numerous independent replicas of the simulation system are run with different starting conditions, provides the means to quantify the aleatoric uncertainty through ensemble averaging, thereby creating reproducible estimates with known errors [19]. We have applied such an ensemble-based approach to determine the glass transition of epoxy resins [20].

In this study, we rigorously determine the parameters required to produce estimates of the elastic moduli of the epoxy resin systems, such as the number of replicas and the length of simulation. Indeed, we show that using a small number of replicas, or short simulation times, leads to unreproducible estimates of the macroscopic properties.

2 Methods↩︎

2.1 Epoxy-graphene nanocomposites↩︎

The investigated systems are cross-linked polymer chains consisting of the epoxy diglycidyl ether of bisphenol F (DGEBF) and the amine 4,4’-methylene-bis-(2,6-diethylaniline) (MDEA) (see figure 1.a). The chemical structures of the two monomers can be found in figure 1.c,d. In addition to a system of pure epoxy, a pristine non-oxidised grapheene flakes stack or tactoid is embedded inside the simulation box (see figure 1.b). This graphene tactoid comprises \(10\) layers of graphene with a radius of \(25\) Å. Due to the polymeric material being studied, the systems were parametrised with the class II polymer consistent force field (PCFF+) [21]–[23]. The pair-wise non-bonded interactions between the graphene sheets are also modelled using the all-atom PCFF+ parameters. It has been shown in prior studies that PCFF+ outperforms other all-atom force fields, such as CHARMM and Amber [24], when calculating material properties such as the glass transition temperature.

The modelling of the epoxy-graphene nanocomposites is proceeded as follows. First, one graphene tactoid is randomly inserted inside an empty isotropic simulation box. Then, the epoxy and amine monomers are randomly packed around the graphene tactoid. The epoxy amines were cross-linked using the MedeA [25] thermoset builder software. The process of building the thermoset alternates between energy minimization, molecular dynamics and topology modifications forming new cross-links between the epoxy and amine species whenever their proximity satisfies a defined distance criterion. The dynamics simulation of the cross-linking of the monomers is achieved in an NVT ensemble at \(700\) K to ensure that the monomers were free to move and form the necessary bonds. The simulation was stopped when \(90\) % of the epoxy groups had formed bonds with the amines. Due to the likelihood of ring catenation, a higher bond conversion was not used. The pure epoxy system is obtained from the graphene systems by removing the graphene sheets and re-equilibrating the system. After cross-linking, arriving at the final polymer, the system was then equilibrated at \(700\) K by \(1\) ns of NVT followed by \(10\) ns of NPT at \(1\) atm.

2.2 Hydration↩︎

Hydration is modelled inserting water molecules to achieve a given weight ratio (see figure 1.b). Water molecules are modelled using the rigid SPC/E water model run under the SHAKE algorithm [26]. Following the initial equilibration of the system, water molecules were inserted via random packing at \(700\) K to achieve specified weight percentages (wt%). To get a statistically robust sample of the epoxy and epoxy-graphene nanocomposite systems, for each water content weight, using different random seeds during the random insertion, ten replicas were produced where the water molecules were packed in different configurations, and the velocities of molecules were initialised with different random seeds drawn from a Maxwell-Boltzmann distribution [20].

2.3 Glass transition↩︎

Through simulated annealing, after the water molecules were inserted, a step-wise glass transition temperature \(T_g\) simulation was run from \(700\) K to \(300\) K. The \(T_g\) simulation used a step size of 10 K, where each step comprised \(0.5\) ns of NVT followed by \(2\) ns of NPT at \(1\) atm.

Using atomistic simulations, polymers have long relaxation times, requiring long simulations of longer than 100 ns of NPT to reach a state where the density reaches a maximum and is constant [20]. In addition to calculating the glass transition, the simulated annealing serves two purposes; not only does it lead to a series of ten statistically different starting configurations for the elastic simulations, but through a combined duration of 90 ns of NPT, the simulated annealing also directly ends on a value for the density of the polymer at 300 K close to equilibrium.

We have generated \(10\) replica during the hydration simulation, this means that \(10\) step-wise glass transition temperature simulations were performed for each wt%, one for each water configuration.

The glass transition temperature is calculated through a hyperbola fitted across all density measurements, using the method of [27]. The density is measured at \(26\) different temperatures spanning from \(700\) K to \(300\) K. Each density at a given temperature is an ensemble average of 10 replicas after 2 ns of sampling time. Denoting \(\rho\) to be the average density at a temperature point \(T\), \[\rho (T) = \rho_0 - a \left(T - T_0 \right) - b \left(~ \frac{1}{2}\left(T - T_0 \right)+\sqrt{\frac{\left(T-T_0\right)^2}{4}+e^c} ~\right) \label{patrone}\tag{1}\] where \(T_0,\rho_0,a,b,c\) are positive real constants and \(e\) is Euler’s constant. The glass transition temperature is given by \(T_0\), representing the hyperbola’s centre. The hyperbola is fit using the least squares algorithm.

2.4 Mechanical deformation↩︎

We calculate the elastic constant through finite deformations at a fixed temperature, averaged from deformations carried out positively and negatively in six directions. Let \(C_{jl}\in \mathbb{R}^{6\times6}\) be the 36-valued elastic constant tensor, then we calculate the modulus as follows: \[\begin{align} T_{11} &= \frac{C_{11}+C_{22}+C_{33}}{3}, \quad T_{12} = \frac{C_{12}+C_{13}+C_{23}}{3}\\ K &= \tfrac{1}{3}(T_{11}+2T_{12}) \\ G_1 &= \tfrac{1}{3}(C_{44}+C_{55}+C_{66})\\ G_2 &= \tfrac{1}{2}(T_{11}-T_{12}) \\ \nu &= (1+T_{11}/T_{12})^{-1} \\ E &= 3K(1-2\nu) \end{align}\] where \(K\) is the bulk modulus, \(G_1\) and \(G_2\) are the shear modulus, \(\nu\) is the Poisson’s ratio, and \(E\) is the Young’s modulus. The values \(T_{11}\) and \(T_{12}\) are intermediate values to aid in the calculation of the moduli. Due to the symmetry of the elastic constant, we average the non-diagonal component pairs of the elastic constant tensor when calculating the moduli.

The elastic deformation simulations are performed under the NVT ensemble at \(300\) K with a time step of \(1\) fs. A strain of amplitude \(1\) % is applied over \(10\) ns in the different directions, in compression and tension, and the resulting stress is sampled over \(1\) ns. The validation of the mechanical deformation simulation protocol is supported by the supplementary figure \(1\).

2.5 Ensemble-based simulations↩︎

In short, \(10\) replicas at the hydration stage, each replicated \(10\) times at the mechanical deformation stage (see figure 2).

The velocities of the system are then reinitialised at \(300\) K with different random seeds to generate \(10\) further replicas for that water configuration. This leads to a total of \(100\) different replicas for each material.

3 Results and discussions↩︎

3.1 Molecular structure of water↩︎

We investigate the structure of water molecules in the pure epoxy resin two ways (see figure 3). First, by means of the partial radial distribution functions (RDF) of the oxygen atoms in water molecules with respect to the oxygen atoms in water molecules or in epoxy groups. Second, by means of a water molecules cluster size analysis. All RDF become smoother with increasing hydration, as statistics improve with the increasing number of water molecules. The RDF of water molecules oxygen atom with themselves shows a sharp first peak around \(2.8\) Å which position remains stable with hydration. The height of the peak, however, decreases, and the right side of the peak tend decrease slower with the distance between oxygen atoms \(r\). As water molecules are being added, second and third peak, respectively at \(4.5\) Å and \(6-7\) Å in bulk water, appear. With increasing hydration, water molecules therefore start to aggregate and start to exhibit structural properties from bulk water. The RDF of water molecules oxygen and epoxy groups oxygen displays a similar trend, the first peak remains at \(2.8\) Å and its intensity decreases with hydration, yet, at longer distances no changes are observed with increasing hydration.

The cluster size distribution (see figure 3) reveals predominantly small water clusters (less than \(10\) molecules) at low hydration (\(1\) wt%). Water is mostly bound to high-affinity polymer sites. At higher hydration, much larger water clusters are found as the high-affinity polymer site become saturated with water, and interwater interactions begin to dominate. As described earlier, the observed cross-over from predominately small, dispersed water clusters at low water uptake to larger clusters can be explained by competing affinity and geometric effects. At low water concentrations, water molecules are efficiently solvated by polar sites in the cross-linked epoxy (espcially the hydroxyl groups), forming small hydrogen bonded complexes. As these sites approach saturation, additional water must occupy other environments and the water clusters grow.

We compare the water molecular structure at \(5\)% hydration between the pure epoxy resin and the graphene-epoxy nanocomposite (see figure 4). RDF between water molecules oxygen atoms and a variety of different atom types in the epoxy chains are not affected by the insertion of graphene. The water molecules are primarily hydrogen bonding to other water molecules or the epoxy groups on the polymer chains. This demonstrates that the presence of graphene does not affect the local hydration environment around the polar oxygen atoms of the epoxy polymer, nor that, as expected, the water molecules do not reside near the graphene sheets.

On the contrary, longer-range order observations, from the cluster size analysis, display a clear divergence, as the graphene containing systems contain a higher fraction of larger water clusters compared to the pure epoxy. We hypothesise that, although the graphene tactoid only occupies a relatively small volume, its high surface to volume ratio decreases the area available for water dispersion, thereby increasing the water-water interactions.

3.2 Density and glass transition↩︎

We first focus on the glass transition phenomena, which we investigate observing the evolution of the density with temperature in our molecular models (see figure 5). In all molecular models, the density decreases nonlinearly with the temperature, displaying an inflexion point associated with the glass transition. Note that the standard error of density measurements over the \(10\) replicas is displayed in figures 5.a,b, but is too small to be visible. The presence of water inside the epoxy resins and the epoxy-graphene nancomposites systematically decreases the density across the whole temperature range. This in emphasised at higher temperatures. Meanwhile, the addition of graphene tends to decrease the density of the material, independentely of the water content.

The hyperbola fitting of the evolution of the density with temperature (see equation 1 ) provides us with the glass transition temperature at different water contents (see figure 5.c). The error bars refer to the error in the least squares fit used to optimise the glass transition temperature parameter, \(T_0\), from the hyperbola in equation 1 . We compute glass transition temperatures in agreement with experimental data for pure dry epoxy resins, around \(410\) K [28]. The increase of the water content up to \(3\)% reduces the transition temperature of epoxy resins down to \(380\) K, then the transition temperature remains stable up to \(5\)% water content. These results are consistent with former molecular simulation studies [4]. Identical variations are observe when epoxy-graphene nanocomposites are considered. The inclusion of the graphene tactoid inside the molecular model of the epoxy resin does not modify the glass transition temperature.

We used a tactoid of ten sheets in our simulation to match the number of graphene sheets typically found in commercially available graphene. Our results show that the low interfacial area associated with this aggregated state and the lack of interaction between the graphene and the epoxy chains results in no discernable change in the glass transition temperature. Furthermore, the small difference in predicted glass transition temperatures between the epoxy polymers and graphene nanocomposites aligns well with results from comparable computational studies [29] of single replica simulations of atomistic dry epoxy polymer and dry graphene nanocomposite. Curiously, the absence of any effect from the insertion of graphene is contrary to what is expected. However, the potentially reduced mobility of the cross-linked epoxy chains in the immobilised interphase region near the graphene surface is indeed expected to yield an increase in the glass transition temperature.

Considering the effect of water molecules, we explain the decrease in glass transition temperature with increased mobility induced by the presence of water molecules. Firstly, penetation of water molecules into the epoxy matrix creates greater free volume for the epoxy chains to move. In addition, water molecules can compete with inter- and intra-chain hydrogen bonds, which will also allow the epoxy chains to move more freely. Since the glass transition temperature plateaus at \(3\)% water, the interrupting of hydrogen bonding within the epoxy chains is likely to be the dominant factor in the reduction at low water contents. At \(3\)% water content, the water molecules change from being mostly bound to the epoxy polar sites to becoming free water, thereby reducing the effective amount of chain-solvating water (see section 3.1), therefore additional water does not increase chain mobility. If the decrease in glass transition temperature was dominated by additional free-volume, the glass transition temperature would be expected to continue decreasing past \(3\)% water content as free-volume increases.

3.3 Mechanical properties↩︎

The various elastic moduli introduced in section 2.4 of the epoxy resin and the epoxy-graphene nanocomposite are computed as a function of the water content (see figure 6. Our baseline prediction of the Young’s modulus of the dry epoxy resin at approximately \(2.6\) GPa agrees with experimental data.

The Young’s modulus is similar for both systems, indicating minimal effect of the graphene sheets. As we discussed in our previous study, this is not unexpected as graphene sheets do not show significant enhancements to the Young’s modulus until they have lateral dimensions of order \(250\) nm, due to shear-lag effects and the lack of strong interactions between epoxy polymers and the graphene sheets [12]. The shear modulus, largely governed by the behaviour of the epoxy matrix, remains similiar between graphene and pure epoxy systems. However, we see significant differences in the bulk modulus and Poisson ratio with the addition of graphene. In the dry-state, the graphene-epoxy nanocomposite exhibits a higher bulk modulus than the neat epoxy. This arises because the graphene sheets act as incompressible inclusion that restrict volumetric deformation, increasing the effective bulk modulus. Overall, the relative behaviour of dry and hydrated systems remain consistent across the whole hydration range.

The evolution of mechanical properties with hydration show that both neat epoxy and graphene-epoxy remain stable upto \(2\)% water content. Beyond this threshold, the water molecules are no longer primarily bound to the polar sites on the epoxy chains; the larger water clusters plasticizes the matrix reducing the Young’s, bulk and shear moduli for both graphene containing and pure epoxy systems. Inversely, the Poisson ratio increases above 2 wt% water content, with the graphene system consistently higher. The increase in Poisson ratio for graphene inclusions can be rationalized via the relation: \(\nu = (3K - 2G) / (2(3K+G))\) ; the elevated \(K\) relative to \(G\) increases \(\nu\). The relation \(E=3K(1-2\nu)\) shows that the effects of higher \(K\) and \(\nu\) for graphene partially offset each other, leading to the Young’s modulus remaining unchanged.

4 Uncertainty quantification in heterogeneous molecular models↩︎

Polymers exhibit long relaxation times with respect to molecular dynamics accessible timescales, as such methods used for equilibration can lead to resulting molecular models with varying densities and structural configurations. To get a statistically representative sample of starting configurations, the starting structures of the mechanical testing simulations were the the result of an independent equilibration for over \(100\) ns. Although we used unconventionally large ensembles of replica (\(100\) replica) with very similar densities at \(300\) K, one striking result from our mechanical simulations is the large variability of the mechanical properties predictions (see error bars in figure 6). Therefore, it appeared essential to investigate in further detail the variability of our mechanical predictions across large ensembles. We now investigate the effect of ensemble size up to \(1000\) replicas for the highly hydrated epoxy (\(5\)% water content).

We computed the Young’s modulus for three ensemble sizes: \(20\), \(100\) and \(1000\) replicas (see figure 7). The ensembles distributions are primarily described by means of histograms which the bin size is calculated using the Freedman–Diaconis rule [30]. Quantile-quantile plots of the distributions are also provided to assess the deviation from a Gaussian distribution. The histogram for \(20\) replica appears to follow a bimodal distribution, with samples less likely to display the average Young’s modulus of the ensemble at approximately \(2.3\) GPa. Even for an ensemble of \(100\) replica, the most probable value of the replicas Young’s modulus exceeds the mean value of the ensemble. The distribution also appears skewed, biased toward higher values of the modulus. Only for a \(1000\) replica ensemble, the distribution of the Young’s modulus values resemble a Gaussian distribution, showing symmetry and centered around the mean value.

Other factors not covered in this study can influence the uncertainty of the simulation. The two most notable are the interatomic force field of the system and the size of the system [31]. The force field used was PCFF+, which has been shown to be significantly better than other classical all-atom force fields for modelling material properties such as the glass transition [24]. A former study demonstrated that the Young’s modulus of molecular models of pure epoxy only become independent of the model size when the simulation box lengths is above 40 Å [18]. In turn, the present molecular models, approximately \(100\) Å wide, should not display size-induced variability of the mechanical properties.

5 Conclusion↩︎

The present study aimed at investigating in depth the influence of hydration on epoxy-based materials, and in particular their density, glass transition and mechanical properties. We have built large ensembles of molecular model of hydrated epoxy resins and epoxy-graphene nanocomposites. From these models, a campaign of MD simulations enables us to predict the systems’ thermodynamic properties as well as the organisation of water molecular.

We have found that the main structural consequence of hydration was the growing size of water clusters in the resin and the nanocomposite. The addition of graphene reinforced the aggregation of the water molecules, with double the size of the largest water clusters produced in presence of graphene.

Meanwhile, the short-range order of water molecules was neither perturbed by the increasing water content nor the insertion of graphene. We have also found that embedding graphene sheets inside an epoxy resin affects the bulk modulus and the Poisson ratio of the result nanocomposites. The average predicted moduli of the nanocomposite were systematically higher than that of the epoxy resin, independently of the water content. Conversely, the addition of graphene did not affect the glass transition temperature. Nonetheless, and consistently with structural observations, hydration diminshes the mechanical properties of resins and nanocomposites in a similar fashion. The effect of hydration was particularly striking beyond \(2\)% to \(3\)% hydration. Last, through uncertainty quantification of very large ensembles, we demonstrated that hundreds of replicas are necessary to report converged distributions of elastic properties and, by extension, reproducible mechanical results.

Overall, our study emphasizes and clarifies the role of water molecules and hydration on the reduction of elastic properties of epoxy-based resins and nanocomposites. Our study also highlights the limited capability of graphene to mitigate the influence of hydration. Graphene-reinforced epoxy resins indeed display higher mechanical properties, but suffer from identical softening as their pristine counterpart.

6 Acknowledgements↩︎

This research was equally co-funded by the Engineering and Physical Sciences Research Council (EPSRC) and Hexcel. We acknowledge funding support from the European Union’s Horizon 2020 research and innovation program under Grant Agreement 800925 (VECMA project, www.vecma.eu), from the UK Engineering and Physical Sciences Research Council under Grant Agreements EP/W007711/1 (SEAVEA project, www.seavea-project.org), and the European Union’s Horizon 2020 Research and Innovation Programme under grant agreement 823712 (CompBioMed2, compbiomed.eu).
We are grateful to the United States Department of Energy (DOE) via the Argonne Leadership Computing Facility for providing access to and node hours on Aurora under a DOE INCITE award in the years 2024-2025.

References↩︎