Quantitative relationship between surface chemistry of graphene and compatibility with rubbers established by two-dimensional solubility parameters

Abstract

Hildebrand (δ_T) and Hansen (δ_D, δ_P, δ_H) solubility parameters are important indexes to predict the compatibility of components intuitively. Currently, almost all the experiments only measured the solubility parameters of the pristine graphene. Therefore, there is a lack of quantitative relationship between the surface chemistry of graphene and solubility parameters, resulting in no theoretical guidance for the surface modification of graphene. In this work, three-dimensional Hansen solubility parameters are converted to two-dimensional solubility parameters. Hildebrand and two-dimensional solubility parameters of six functionalized graphene as a function of grafting ratio are calculated by molecular dynamics (MD) simulation. Interestingly, if the functional group is at the edge of graphene, the δ_T decreases with the increase of the grafting ratio, whereas if the functional group is in the plane of graphene, the δ_T decreases first and then increases with grafting ratio. Two-dimensional solubility parameters are proved to be a good predictor of the compatibility between functionalized graphene and rubbers. The quantitative relationship between the surface chemistry of graphene and compatibility with rubbers based on two-dimensional solubility parameters is constructed. The optimum grafting ratio corresponding to the best compatibility is given. Finally, the effect of temperature on the compatibility behaviors of graphene/rubber mixtures is elucidated.

---

1 Introduction

Graphene is considered to be an ideal filler for the rubber industry due to its extremely high specific surface area, thermal and electrical conductivity.^1–3 The graphene–rubber interactions and the dispersion of graphene in the matrix are two key factors affecting the properties of graphene/rubber composites.⁴ From a thermodynamic point of view, the compatibility of graphene and rubber plays a key role in the above two factors.⁵ Since pristine graphene is neither hydrophilic nor oleophilic, the compatibility of pristine graphene and rubber is poor. Therefore, it is necessary to modify the surface chemistry of graphene to improve the thermodynamic compatibility of graphene and rubber. One of the most commonly used methods for the preparation of graphene is the oxidation–reduction method. The introduction of oxygen-containing groups (such as hydroxyl, carboxyl, epoxide groups, etc.) by the oxidation–reduction method is unavoidable.^6–8 These oxygen-containing groups provide the possibility for the modification of graphene. At present, there are some mature modification methods to prepare aminoacylated,⁹ carboxylated,¹⁰ sulfhydrated,^11,12 alkylated,¹³ and silanized graphene.¹⁴

The most commonly used index to predict the thermodynamic compatibility of two components is the solubility parameter proposed and developed by Joel Henry Hildebrand and Charles M. Hansen.^15–17 The thermodynamic theory of mixing shows that the miscibility of two components is related to their cohesion energy density (CED), and the CED is the cohesion energy per unit volume (E_coh/V).¹⁸ The Hildebrand solubility parameter (δ_T) is defined as the square root of CED. Charles M. Hansen divides cohesion energy into dispersion (E_coh,D), polarity (E_coh,P), and hydrogen bonding (E_coh,H) interactions:¹⁷

Therefore, the Hansen solubility parameters (δ_i, i = D, P, H) are proposed:

Then the relation between Hildebrand and Hansen solubility parameters is as follows:

δ_T² = δ_D² + δ_P² + δ_H² (3)

Generally, the Hildebrand solubility parameter is a good predictor for the compatibility of non-polar systems, while the Hansen solubility parameter is proved to be more suitable for the judgment of the compatibility of polar systems.¹⁷ Besides, in the past decade, some studies have reported the solubility parameters of carbon materials (graphene,¹⁹ carbon nanotube,^20,21 fullerenes,²² etc.), and the compatibility between carbon materials and polymers is well predicted by solubility parameters.^21–23 Such as, Hernandez et al.¹⁹ measured the solubility parameters of graphene through dissolution experiments: δ_T = 23 MPa^1/2, δ_D = 18 MPa^1/2, δ_P = 9.3 MPa^1/2, and δ_H = 7.7 MPa^1/2. To solve the problem of poor compatibility between graphene and polyvinylidene fluoride (PVDF), Zhu et al.²⁴ used polyimide (PI) and polyaniline (PANI) to modify reduced graphene oxide (RGO), respectively, and measured Hansen solubility parameters of functionalized RGO. Further research shows that the matching of Hansen solubility parameters is crucial for the dielectric and dynamic mechanical properties of RGO/PVDF composites. To the best of our knowledge, although some experiments have measured the solubility parameters of functionalized graphene, the quantitative relationship between functional groups and the solubility parameters of graphene is still lacking due to the limitations of modification and characterization methods.^25–27 Therefore, the current surface modification of graphene relies on trial and error method, and there is no new theory to guide the interface design of high-performance graphene/rubber composites.

Alternatively, molecular dynamics (MD) simulation based on Newton's laws of motion is already a powerful computational tool that can provide an important means for understanding and predicting the behaviors/properties of materials.^28–31 In the field of polymer composites, MD simulation can be used to study the mechanical, thermal, and electrical properties of polymers (e.g., natural rubber, styrene butadiene rubber, polyethylene, epoxy) composites filled with graphene.^32,33 Compared with experimental methods, MD simulation can easily and cheaply establish quantitative relationships between the structure and properties of materials. Such as the experimental measurement of solubility parameters requires dozens of solvents and consumes a lot of time,²³ but the research cost and time can be greatly reduced by MD simulation.³⁴ Especially, with the improvement of the force field precision, the solubility parameters of small molecules, polymers, and carbon materials calculated by MD simulation have been in good agreement with the experimental values.^34–36 Therefore, MD simulation is expected to be an effective method to explore the relationship between functional groups and solubility parameters of graphene. Besides, MD simulation is expected to be a powerful means to establish the quantitative relationship between surface chemistry of graphene and compatibility with polymers.

In this work, three-dimensional Hansen solubility parameters were converted to two-dimensional solubility parameters based on the force field method. Six functionalized graphene models were constructed by introducing six functional groups (phenyl isocyanate, 1-hexyl isocyanate, aniline, 1-octadedecylamine, 2-mercaptoethylamine, and 3-mercaptopropionic acid). Then the effect of grafting ratio on Hildebrand and two-dimensional solubility parameters of functionalized graphene was investigated by MD simulation. Whether two-dimensional solubility parameters can effectively predict the compatibility between graphene and rubbers with different polarities was studied based on the Flory–Huggins lattice model. The functionalization principle and the optimal grafting ratio of graphene for graphene/rubber mixtures were further explored. Additionally, the effect of temperature on the compatibility behaviors of graphene/rubber mixtures was elucidated.

2 Methodology

2.1 Energy and two-dimensional solubility parameters

In MD simulation, the ab initio Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) force field is one of the most widely used force fields for calculating solubility parameters of polymers and carbon materials.^28,37,38 The simulated results based on COMPASS force field have been widely proved to be in good agreement with the experimental values. Therefore, the solubility parameters of graphene and rubbers are also calculated using the COMPASS force field in this work. In the COMPASS force field, the potential energy of a system (E_total) can be expressed as a sum of valence (E_valence), cross term (E_crossterm), and nonbond (E_nonbond) interactions:³⁹

E_total = E_valence + E_crossterm + E_nonbond (4)

The valence interactions are generally accounted for by bond stretching, valence angle bending, dihedral angle torsion, out-of-plane interactions, and Urey-Bradley terms. The cross term is introduced to account for such factors as bond or angle distortions caused by nearby atoms so that higher accuracy can be achieved. The energy of interactions between nonbonded atoms (E_nonbond) is accounted for by van der Waals also known as dispersion (E_vdW) and electrostatic (E_elec) terms:

E_nonbond = E_vdW + E_elec (5)

E_vdW is calculated by the Lennard-Jones 9–6 function (eqn (6)), and E_elec is calculated by the CONST-ε function (eqn (7)).where R is the distance between two atoms, D₀ is the equilibrium well depth, and R₀ is the equilibrium distance.where R is the distance between atom i and atom j. q is the charge of the atom. ε is the relative dielectric constant and R is the distance between two atoms. C = 332.0647 kcal mol⁻¹ Å e⁻² is a unit conversion factor.

So far, no force field can accurately express hydrogen bond interactions. Some force fields (such as Dreiding), while giving a separate hydrogen bond term, are based on low-precision empirical or semi-empirical force field parameters.⁴⁰ Therefore, in order not to sacrifice computational accuracy, the COMPASS force field does not contain a separate hydrogen bond term, but instead merges the polar and hydrogen bond terms into the electrostatic term:

E_elec = E_coh,P + E_coh,H (8)

Then the simulated solubility parameters based on the COMPASS force field can be given by

δ_T² = δ_vdW² + δ_elec² (9)

δ_elec² = δ_P² + δ_H² (10)

Consequently, three-dimensional Hansen solubility parameters are converted to two-dimensional solubility parameters based on the force field method. We present a chart to show the difference between the two solubility parameters so that readers can better understand, as shown in Fig. 1. Gupta et al.²⁸ successfully predicted the miscibility of pharmaceutical compounds by similar two-dimensional solubility parameters from MD simulation. Lee et al.³⁴ determined two-dimensional solubility parameters of single-walled and double-walled carbon nanotubes using computational simulation. In this work, we calculate two-dimensional solubility parameters of graphene and rubbers by COMPASS force field. Also, whether two-dimensional solubility parameters can predict the compatibility of graphene and rubber is investigated.

Fig. 1 The definition of three-dimensional and two-dimensional solubility parameters.

2.2 Model and simulation details

The functionalization of graphene is often achieved through the reaction of oxygen-containing groups on the surface, such as the hydroxyl group in the graphene plane and the carboxyl group at the edge.⁴¹ Therefore, according to the reports in the literature,^9–13 we selected the modifiers containing isocyanate groups (phenyl isocyanate and 1-hexyl isocyanate) that can react with hydroxyl groups and the modifiers containing amino groups (aniline and 1-octadedecylamine) that can react with carboxyl groups, as shown in Fig. 2a–d. The isocyanate group reacts with the hydroxyl group to form the urethane group, and the carboxyl group reacts with the hydroxyl group to form the ester group. In addition, another popular modification method is to functionalize graphene by click chemistry, so we select modifiers containing mercaptan groups (2-mercaptoethylamine and 3-mercaptopropionic acid) to perform a click chemical reaction with C=C on graphene, as shown in Fig. 2e and f. Notably, it's not how important these groups are, but these chemical modification methods are mature. Therefore, we selected these representative compounds containing isocyanates, amino groups, and thiols to modify graphene. The purpose of our study is not to emphasize how important these groups are for the properties of graphene/rubber composites but to construct a quantitative relationship between functional groups and the compatibility of functionalized graphene and rubbers.

Fig. 2 Chemical structure of (a) phenyl isocyanate, (b) 1-hexyl isocyanate, (c) aniline, (d) 1-octadedecylamine, (e) 2-mercaptoethylamine, (f) 3-mercaptopropionic acid, (g) nature rubber (NR), (h) butadiene rubber (BR), (i) neoprene rubber (CR), (j) nitrile butadiene rubber (NBR), and (k) styrene butadiene rubber (SBR). Models of (l) pristine graphene, (m) functionalized graphene in the plane, and (n) functionalized graphene at the edge. Amorphous cells of (o) BR and (p) graphene/BR composites. The grey, white, blue, red, and yellow balls represent C, H, N, O, and S atoms, respectively.

To investigate the general law for the compatibility of graphene surface chemistry with rubber, five different polarities of rubber were constructed: natural rubber (NR), butadiene rubber (BR), neoprene rubber (CR), nitrile butadiene rubber (NBR), and styrene butadiene rubber (SBR), as shown in Fig. 2g–k. The mass fraction in the MD simulation is the same as that of a commercial NBR (N220S, Japan Synthetic Rubber Co., Ltd) with 41 wt% of acrylonitrile and 59% of butadiene.⁴² The simulated chemistry structure of SBR with 21.0 wt% of styrene, 47.4 wt% of 1,2-polybutadiene, and 31.6 wt% of 1,4-polybutadiene, which is consistent with that of a commercial SBR (2466, China Petroleum & Chemical Corporation).⁴³ The model of a single graphene with a length of 4.182 nm and a width of 3.690 nm was constructed. The unsaturated boundaries of graphene were terminated by adding hydrogen atoms. The pristine graphene and functionalized graphene are shown in Fig. 2l–n.

Amorphous cells (Cell-1) containing two graphene or functionalized graphene sheets were built to simulate the solubility parameters of graphene or functionalized graphene, as shown in Fig. S1 (ESI).† Amorphous cells (Cell-2) of rubber consisting of 10 chains with 50 repeating units are constructed to calculate the solubility parameters, as shown in Fig. 2o. Amorphous cells (Cell-3) consisting of one graphene or functionalized graphene sheet and 10 rubber chains are constructed to calculate the mixing energy (E_mix) and Flory–Huggins parameter (χ), as shown in Fig. 2p. This Cell-3 is constructed according to the following algorithm: the construction algorithm resembles that of a polymerization process, where molecules are “grown” into an empty box fragment by fragment. A fragment corresponds to the substructure of atoms between two rotatable bonds. Amorphous cell automatically detects rotatable bonds along the backbone of the input model. In adding the next fragment, the algorithm probes a number of potential placements. The energy of each placement is evaluated from which the probability of each placement is calculated. Based on the probabilities, one of the candidates is selected, and the fragment is given its final placement. The torsion degrees of freedom in the backbone of input molecules will be determined automatically. The initial temperature is set to 298 K and the initial density is set to 1 g cm⁻³. The number of segments to consider in growing the chains into the cell is 1. The maximum number of attempts made to load all molecules into the cell is 1000. A test on spearing and catenation is carried out whilst growing the chains into the cell. The number of random positions and orientations to be sampled for the first segment of each molecule is 10. The number of torsions to be sampled for all torsions of each molecule is 10. The number of molecules in the amorphous cell and the initial conditions for the construction of the cell are listed in Table 1.

Table 1 The number of molecules in the amorphous cell and the initial conditions for the construction of the cell

	Graphene sheet	Rubber chains	Initial density (g cm^−3)	Initial temperature (K)	Lattice type
Cell-1	2	0	1	298	Cubic
Cell-2	0	10	1	298	Cubic
Cell-3	1	10	1	298	Cubic

---

In MD simulation, three-dimensional periodic boundary conditions are used. The amorphous cell is first geometrically optimized by the Smart method with an energy convergence tolerance of 10⁻⁴ kcal mol and a force convergence tolerance of 0.005 kcal mol⁻¹ Å⁻¹. The resulting cell is then subjected to 1 ns NVT (constant number of atoms, volume, and temperature) ensemble simulation at 298 K, followed by 2 ns NPT (constant number of atoms, pressure, and temperature) ensemble simulation at 298 K and 10⁵ Pa, and finally, the equilibrium structure is obtained. The average of the last five frames was calculated to obtain the values of δ_T and χ. In the simulation, the time step is 1 fs and the cutoff radius of the van der Waals interaction is 1.25 nm. The Ewald method was used to calculate the electrostatic interactions with an accuracy of 1.0 × 10⁻⁴ kcal mol⁻¹ and a buffer width of 0.5 Å. The Berendsen barostat was used to control pressure, and the Nosé–Hoover thermostat was used to control temperature. The pressure is applied isotropically to the systems. The value used to scale the fictitious mass, Q, of the Nosé–Hoover thermostat is 0.01. The coupling constant is 1 ps for both Berendsen barostat and Nosé–Hoover thermostat. All modeling and simulations were operated using Materials Studio software with COMPASS force field.

Cell-1 is the amorphous cell of graphene or functionalized graphene. Cell-2 is the amorphous cell of rubbers. Cell-3 is the amorphous cell of graphene/rubber mixtures.

To verify the reliability of the simulation method and results, the simulated density (ρ) and solubility parameters were compared with the experimental values, as shown in Table 2. It can be seen that the simulated ρ and δ_T of rubbers are in good agreement with the experimental values. However, the simulated δ_T = 29.2 MPa^1/2 shows some difference from the experimental δ_T = 23.0 MPa^1/2 for graphene. It is worth noting that we construct a monolayer of graphene without any oxygen-containing group and defect in simulation. Nevertheless, the dissolution experiment of graphene found that even if the best solvent (N-methyl-2-pyrrolidone, NMP) was used to dissolve graphene, the number fraction of single-layer graphene was only 29%, and the maximum number of layers of graphene was 16.¹⁹ Besides, the graphene prepared by existing methods is almost impossible to avoid the existence of defects and functional groups. Therefore, we speculate that the reasons for the difference between the simulated and experimental δ_T are as follows: (1) defects, (2) functional groups, and (3) number of graphene layers.

Table 2 Simulated and experimental δ_T and ρ of rubbers and pristine graphene^a

	δT,s (MPa^1/2)	δvdW,s (MPa^1/2)	δelec,s (MPa^1/2)	δT,e (MPa^1/2)	ρs (g cm^−3)	ρe (g cm^−3)	Ref.
a Subscript s represents the simulated value and subscript e represents the experimental value.
NR	16.3	15.8	1.1	16.5	0.93	0.92–0.95	17
BR	16.6	15.4	0.9	16.5	0.90	0.90–0.92	17
SBR	16.9	16.2	2.2	16.7	0.97	0.95–1.00	17
NBR	20.8	18.5	8.0	21.1	1.04	1.06	17
CR	19.3	17.8	5.2	19.0	1.22	1.23–1.25	17
Pristine graphene	29.2	28.1	0.02	∼23.0			19

---

The best-known theory of the thermodynamics of mixing and phase separation in binary systems is the Flory–Huggins lattice model.^44,45 In the Flory–Huggins lattice model, the χ of the binary system can be calculated by MD simulation to characterize the compatibility of graphene and rubber:

where E_mix is the mixing energy; that is, the difference in free energy due to interaction between the mixed and the pure state. k_B is the Boltzmann constant. In the Flory–Huggins model, each component occupies a lattice site. For a lattice with coordination number Z, the E_mix iswhere E_ij (i, j = b or s) is the binding energy between a unit of component i and a unit of component j. For molecules, the binding energies are regarded as averages over an ensemble of molecular configurations.

3 Results and discussion

3.1 Validation of models and methods

To further verify the validity of the model and method, different simulation results are compared in different cell sizes, simulation times, and initial conditions. To verify whether the system reached convergence after the simulation of 2 ns, we extended the simulation time to 10 ns. The resulting δ_T of NR and pristine graphene as well as the χ value of graphene/NR mixture are presented in Fig. S2 (ESI).† The results indicate that the δ_T and χ do not change even if the time is extended from 2 ns to 10 ns. That is, these systems have converged after a 2 ns simulation. Besides, the energy, density, and temperature as a function of time within 10 ns are given for natural rubber, graphene, and graphene/NR mixture, as shown in Fig. S3–S5 (ESI).† Generally, if the fluctuation of energy, density, and temperature is within 10%, the system is considered to have reached convergence. As a result, the systems reach convergence after the simulation of the NPT ensemble for 2 ns.

The different initial conditions are listed in Table S1 (ESI)† and the simulation results are listed in Table S2 (ESI).† For example, in initial condition 1, the initial density is set to 0.6 g cm⁻³, and Berendsen barostat and Nose–Hoover thermostat are employed. The simulated results are much lower than the experimental values. Compared with the experimental δ_T, it is found that the simulation δ_T using initial condition 2 is inferior to that using initial condition 3, so the simulation in this manuscript is based on initial condition 3.

To consider the system size effect, we doubled the number of molecules in the system listed in Table 1 and kept the other simulation conditions unchanged. The number of molecules in the extended system is listed in Table S3 (ESI),† and the simulated results are shown in Fig. S6 (ESI).† The solubility parameters of pristine graphene and rubbers are constant even though the number of molecules has doubled. Although the χ value of 3-mercaptopropionic acid functionalized graphene/NR mixtures increases with the number of molecules, the changing trend of χ is the same. In summary, the size effect does not affect the trend.

3.2 Solubility parameters of functionalized graphene

The solubility parameters of six functionalized graphene as a function of grafting ratio are shown in Fig. 3. The grafting ratio is defined as the percentage of the number of functional groups to the number of carbon atoms (574) on the graphene. According to reports in the literature,⁴⁶ the oxygen content of graphene oxide prepared by the modified Hummers' method is between 25% and 45%. In the simulation, the maximum grafting ratio is 26.8% and the highest oxygen content is 27.2%. Therefore, the simulated grafting ratio and oxygen content are reasonable. Strikingly, if the functional group (aniline and 1-octadedecylamine) is at the edge of graphene, the δ_T decreases with the increase of the grafting ratio, whereas if the functional group (phenyl isocyanate, 1-hexyl isocyanate, 3-mercaptopropionic acid and 2-mercaptoethylamine) is in plane of graphene, the δ_T decreases first and then increase with grafting ratio. The change of δ_T and δ_vdW is consistent with the grafting ratio. Besides, the δ_ele increases with the increase of grafting ratio. We present the δ_elec on the same ordinate scale for a more intuitive comparison of the magnitude, as shown in Fig. S7 (ESI).† The magnitude of the δ_ele is as follows: 3-mercaptopropionic acid > 2-mercaptoethylamine > phenyl isocyanate > aniline > 1-octadedecylamine > 1-hexyl isocyanate. To save space, functionalized graphene is described here with modifiers. Besides, current experimental means may be able to prepare several modified graphene with different grafting ratios, but it is not possible to precisely prepare all the modified samples shown in Fig. 3. Since the existence of oxygen-containing functional groups in graphene is inevitable, it is difficult to exclude the effect of such oxygen-containing functional groups when we study the quantitative relationship between grafting ratio and solubility parameters. As a result, we cannot quantitatively obtain the quantitative relationship between grafting ratio and solubility parameters by experimental means. MD simulation can accurately construct functionalized graphene with different groups and different grafting ratios and exclude the inevitable effect of oxygen-containing functional groups in experiments. Therefore, our study aims to theoretically provide the quantitative relationship between different functional groups/grafting ratios and solubility parameters through MD simulation. This is also the significance of our study. It is hoped that the effectiveness of our simulation can be verified after the experimental preparation and characterization methods are more mature.

Fig. 3 Solubility parameters of aniline, 1-octadedecylamine, phenyl isocyanate, 1-hexyl isocyanate, 3-mercaptopropionic acid and 2-mercaptoethylamine functionalized graphene with grafting ratio: (a) Hildebrand (δ_T), (b) van der Waals (δ_vdW) and (c) electrostatic (δ_ele) solubility parameters.

To explain the variation trend of δ_vdw and δ_elec, we give the E_vdW and E_elec of 2-mercaptoethylamine and aniline functionalized graphene as a function of the grafting ratio, as shown in Fig. 4a and c. In addition, the volume changes of the two systems with grafting ratio are shown in Fig. 4b and d. For the 2-mercaptoethylamine functionalized graphene, the E_vdW first decreases and then increases. For the aniline functionalized graphene, the E_vdW decreases first and then increases slowly. Both E_elec and V of the two systems increase with grafting ratio. The decrease of E_vdw and the increase of V cause the δ_vdW of two systems to decrease within a certain grafting ratio. Above a certain grafting ratio, the δ_vdw of 2-mercaptoethylamine functionalized graphene increases, while the δ_vdw of aniline functionalized graphene decreases. The variation of δ_vdw can be inferred from the following equation:

where δ_vdW² is the δ_vdw at a high grafting ratio and δ_vdW¹ is the δ_vdw at a low grafting ratio. V₂ and V₁ are the volume at high and low grafting ratio. E_vdW,2 and E_vdW,1 are the E_vdW at high and low grafting ratio. According to the eqn (13), if δ_vdW²/δ_vdW¹ > 1, δ_vdw increases, and if δ_vdW²/δ_vdW¹ < 1, δ_vdW decreases. As can be seen from Fig. 4, although the E_vdW of 2-mercaptoethylamine functionalized graphene increases in the range of 2.6–12.2% grafting ratio, δ_vdw decreases in this range. This is because the increased multiple in the E_vdW is lower than the increased multiple in volume, resulting in a decrease in δ_vdw. As the grafting ratio exceeds 12.2%, the increased multiple of E_vdW is greater than the increased multiple of volume for 2-mercaptoethylamine functionalized graphene, increasing δ_vdw. For aniline functionalized graphene, even if the E_vdW increases, the increased multiple of the E_vdW is lower than the increased multiple of the volume, resulting in a continuous decline in δ_vdw. Besides, since the increased multiple of E_elec is always greater than the increased multiple of volume, δ_elec is always increasing.

Fig. 4 (a) van der Waals interactions (E_vdW), electrostatic interactions (E_elec), and (b) volume (V) of 2-mercaptoethylamine functionalized graphene with grafting ratio. (c) E_vdW, E_elec, and (d) volume of aniline functionalized graphene with grafting ratio.

Reasonably, the π–π interactions between graphene sheets decrease as the grafting ratio increases due to the continuous conversion of sp² hybrid carbon atoms to sp³ hybrids. We would like to give the law of π–π interactions with grafting ratio, but no molecular force field has been able to calculate the π–π interactions so far. Therefore, the reason why the solubility parameter changes with the grafting ratio is also based on the speculation of van der Waals and electrostatic interactions. We analyze the reasons for the solubility parameter variation with grafting ratio as follows: In the graphene system, the E_vdW consists of the π–π stacking interactions (E_π–π) and other van der Waals interactions besides the E_π–π (which we call the regular van der Waals interaction, E_vdW,r).

E_vdW = E_π−π + E_vdW,r (14)

With the increase of grafting ratio, E_π–π decreases, and E_vdW,r increases. At a low grafting ratio, the decreasing effect of E_π–π dominates, so E_vdw decreases. Above a certain grafting ratio, with the increasing number of atoms and the continuous breaking of sp² hybrid carbon atoms, the increasing effect of E_vdW,r rather than the decreasing effect of E_π–π is dominant, so the E_vdw increases. This explains the phenomenon in Fig. 3.

It is worth noting that the minimum δ_T of functionalized graphene is a very important parameter, which determines the closeness of the solubility parameter of graphene and rubber. The minimum δ_T and corresponding grafting ratio for different functional groups are shown in Fig. 5. It can be seen that the solubility parameters of phenyl isocyanate, 1-hexyl isocyanate, aniline, 1-octadedecylamine, 2-mercaptoethylamine, and 3-mercaptopropionic acid functionalized graphene decreased by 17.8%, 52.7%, 45.9%, 50.0%, 50.0% and 46.2%, respectively, compared to the solubility parameter of pristine graphene. This shows that the effect of functional groups on solubility parameters is significant. Besides, the minimum δ_T of functionalized graphene and the δ_T of rubbers are ordered from high to low value, as listed in Table 3. This order determines the thermodynamically compatible behavior of functionalized graphene and rubber, which will be further elaborated in the next section.

Fig. 5 The minimum δ_T and corresponding grafting ratio of functionalized graphene.

Table 3 The minimum δ_T of functionalized graphene and the δ_T of rubbers

Samples	δT (MPa^1/2)
Phenyl isocyanate	24.0
NBR	20.8
CR	19.3
SBR	16.9
BR	16.6
NR	16.3
Aniline	15.8
3-Mercaptopropionic acid	15.7
2-Mercaptoethylamine	14.9
1-Octadedecylamine	14.9
1-Hexyl isocyanate	13.8

---

3.3 Compatibility of functionalized graphene and rubbers

Two-dimensional solubility parameters of rubber and graphene have been obtained. Next, we want to know if two-dimensional solubility parameters can predict the thermodynamic compatibility of graphene and rubber. According to the thermodynamics of mixing, the compatibility of two components is determined by the Gibbs free energy of mixing (ΔG_M):

ΔG_M = ΔH_M − TΔS_M (15)

where ΔH_M and ΔS_M are enthalpy and entropy of mixing, respectively. T is the absolute temperature. Only if ΔG_M is negative can the two components be compatible. If ΔG_M is positive, then the two components are incompatible. ΔS_M is usually positive when the two components are mixed, so keeping ΔH_M as small as possible ensures that ΔG_M is negative. In the Flory–Huggins lattice model, ΔH_M can be expressed:^47,48where φ_A and V_m,A are volume fraction and molecular volume (lattice site volume in the lattice theory) of component A as solvent, respectively. When studying the heat of mixing, Hildebrand and Scratchard proposed an expression for the variation of heat after mixing two components:^47,48

ΔH_M = (δ_T,A − δ_T,B)²φ_A(1 − φ_A)V_M (17)

where δ_T,A and δ_T,B are Hildebrand solubility parameters of components A and B, respectively. V_M the total volume of the mixture. This expression states that the closer δ_T,A and δ_T,B are, the smaller ΔH_M is, making it more possible for ΔG_M to be negative. The relationship between χ and δ_T can be obtained by eqn (16) and (17):

This formula states that the closer δ_T,A and δ_T,B are to each other, the smaller χ is. Charles M. Hansen replaced the Hildebrand solubility parameter with the Hansen solubility parameter as follows:¹⁷

Then, an R₀ parameter is proposed to describe the closeness of the solubility parameters:

The smaller the R₀ value, the closer the solubility parameters, and the better the compatibility of A and B components. In our study, two-dimensional solubility parameters were substituted for three-dimensional Hansen solubility parameters. Similarly, an R parameter is proposed to describe the closeness of the two-dimensional solubility parameters:

Therefore, to explore whether two-dimensional solubility parameters are good predictors of graphene and rubber compatibility, R and χ need to be compared. If R and χ have the same trend, it shows that two-dimensional solubility parameters are a good predictor of compatibility, but the opposite is not.

We investigated the correlation of R and χ with grafting ratio for 30 graphene/rubber mixtures (six functionalized graphene and five rubbers), as shown in Fig. 6 and S8–S12 (ESI).† It was found that R and χ had the same variation trend with grafting ratio, which proved that two-component solubility parameters can predict the compatibility of functionalized graphene and rubber well. Besides, interestingly, the curve of R with grafting ratio shows two different shapes (V-shape and W-shape). That is, the graphene and rubber mixtures exhibit two thermodynamically compatible behaviors, as listed in Table 4. For example, for aniline and 1-octadedecylamine functionalized systems (the δ_T decreases with grafting ratio), the δ_T of functionalized graphene is closest to that of rubber at a certain grafting ratio. Therefore, there is a minimum value for the curve of R with grafting ratio, that is, V-shape. For phenyl isocyanate, 1-hexyl isocyanate, 3-mercaptopropionic acid, and 2-mercaptoethylamine functionalized systems, the δ_T decreases first and then increases with grafting ratio. Therefore, If the minimum δ_T of graphene is greater than the δ_T of rubber, then the minimum δ_T is the closest to the δ_T of rubber. That is, there is only one minimum value for the curve of R with grafting ratio, and the curve is V-shaped. If the minimum δ_T of graphene is less than δ_T of rubber, there are two grafting ratios at which the δ_T of graphene and the δ_T of rubber are the closest. Therefore, there are two minima for the curve of R with grafting ratio, that is, W-shape.

Fig. 6 R and χ of 3-mercaptopropionic acid functionalized graphene/rubber mixtures with the grafting ratio: (a) NR, (b) BR, (c) SBR, (d) NBR, and (e) CR.

Table 4 Compatibility behavior of functionalized graphene and rubbers^a

	NR	BR	SBR	NBR	CR
a V means the curve of R or χ as a function of grafting ratio is V-shaped and W means the curve is W-shaped.
Aniline	V	V	V	V	V
1-Octadedecylamine	V	V	V	V	V
Phenyl isocyanate	V	V	V	V	V
1-Hexyl isocyanate	W	W	W	W	W
3-Mercaptopropionic acid	V	V	W	W	W
2-Mercaptoethylamine	W	W	W	W	W

---

The conclusions in Fig. 6 and Table 4 illustrate the following two points: (1) the trend of R and χ values with the grafting ratio is consistent, indicating that the two-component solubility parameters can well predict the compatibility of graphene and rubbers; (2) different graphene/rubber mixtures exhibit different compatibility behaviors with the grafting ratio, which may be V-shape or W-shape.

According to the eqn (11), the χ originates from the E_mix. That is, exploring the reason for the variation of χ with grafting ratio (V-shape or W-shape) is to investigate the reason for the variation of E_mix with grafting ratio. The nonbond interactions in the COMPASS force field are the sum of the electrostatic and van der Waals interactions. Therefore, E_mix is also divided into van der Waals components (E_mix,vdW) and electrostatic components (E_mix,elec).

E_mix = E_mix,vdW + E_mix,elec (22)

To show why there is such a trend of V-shape or W-shape, the E_mix and its corresponding E_mix,vdW and E_mix,elec with grafting ratio for 1-hexyl isocyanate functionalized graphene/NR mixture, 3-mercaptopropionic acid functionalized graphene/NR mixture and 3-mercaptopropionic acid functionalized graphene/NBR mixture are given, as shown in Fig. 7. The change of E_mix with grafting ratio is W-shaped for 1-hexyl isocyanate functionalized graphene/NR mixture. The change of E_mix with grafting ratio is V-shaped for 3-mercaptopropionic acid functionalized graphene/NR mixture. The change of E_mix with grafting ratio is W-shaped for 3-mercaptopropionic acid functionalized graphene/NBR mixture. It can be seen that for either mixture, E_mix,vdW is positive and E_mix,vdW increases with grafting ratio. E_mix,elec is negative, and the smaller the E_mix,elec, the greater the electrostatic mixing energy. With the increase of grafting ratio, the electrostatic mixing energy first increases and then decreases. At a certain grafting ratio, the electrostatic mixing energy reaches the maximum. This shows that too high grafting ratio will cause the polarity difference between the functionalized graphene and the rubber to be too large, resulting in a decrease in electrostatic mixing energy. It is due to the different increasing (or decreasing) tendencies of E_mix,vdW and E_mix,elec that E_mix appears V-shaped or W-shaped.

Fig. 7 The variation of mixing energy (E_mix) and van der Waals (E_vdW) and electrostatic (E_elec) components of mixing energy with grafting ratio for (a) 1-hexyl isocyanate functionalized graphene/NR mixture, (b) 3-mercaptopropionic acid functionalized graphene/NR mixture and (c) 3-mercaptopropionic acid functionalized graphene/NBR mixture.

3.4 Functionalization principle of graphene

From the above analysis, the compatibility between graphene and rubber is optimal when corresponding to the smallest R or χ. Then, the grafting ratio corresponding to the best compatibility can be called the optimum grafting ratio. The optimum grafting ratios of different graphene/rubber mixtures are listed in Table 5. It can be seen that the optimum grafting ratio may be the same or different for different mixtures, depending on their compatible behavior. Many experiments have also found that the existence of an optimum grafting ratio makes the graphene/rubber composites the best performance.^2,49,50 Our work provides a theoretical basis for the explanation of these phenomena.

Table 5 Optimum grafting ratios of different graphene/rubbers mixtures

	NR (%)	BR (%)	SBR (%)	NBR (%)	CR (%)
Aniline	23.2	23.2	23.2	15.7	19.7
1-Octadedecylamine	12.2	10.5	15.7	7.0	8.7
Phenyl isocyanate	8.7	8.7	8.7	8.7	8.7
1-Hexyl isocyanate	15.7	15.7	15.7	23.2	19.7
3-Mercaptopropionic acid	10.5	10.5	10.5	23.2	19.7
2-Mercaptoethylamine	8.7	8.7	8.7	23.2	19.7

---

From a thermodynamic point of view, the properties of graphene/rubber composites are determined by two aspects: (1) the inherent properties of graphene and rubber, and (2) the compatibility of graphene and rubber. The introduction of functional groups into graphene will improve the compatibility of graphene and rubber; however, this will deteriorate the properties of graphene itself (electrical and thermal conductivity, etc.). Therefore, the functionalization principle of graphene is to make the functionalized graphene and rubber have good compatibility while not destroying the structure of graphene to the greatest extent. That is, the R-value should be as low as possible at a low grafting ratio. The R-values of graphene/rubber mixtures as a function of grafting ratio are compared in Fig. 8. Therefore, according to the functionalization principle of graphene, the priority of functional groups for NR, BR, and SBR is as follows: 1-octadedecylamine > 1-hexyl isocyanate > 3-mercaptopropionic acid > 2-mercaptoethylamine > aniline > phenyl isocyanate. The priority of functional groups for NBR is as follows: 3-mercaptopropionic acid > 1-octadedecylamine > 1-hexyl isocyanate > 2-mercaptoethylamine > aniline > phenyl isocyanate. The priority of functional groups for CR is as follows: 3-mercaptopropionic acid > 1-octadedecylamine > 1-hexyl isocyanate > 2-mercaptoethylamine > phenyl isocyanate > aniline. For NR, BR, and SBR with low polarity, the 1-octadedecylamine with long alkane chains and low polarity is the most preferred group. For high-polarity rubber (NBR and CR), the functional group with high polarity (3-mercaptopropionic acid) is preferred.

Fig. 8 R as a function of grafting ratio for different functional groups: (a) NR, (b) BR, (c) SBR, (d) NBR, and (e) CR.

3.5 Temperature dependence of compatibility

It is widely known that temperature has an important effect on the thermodynamic compatibility of two components.^51,52 Therefore, it is of significance to explore the effect of temperature on the compatibility for the processing of materials. The processing temperature of rubber is usually less than 200 °C. Above 200 °C, the higher the temperature, the less important it is for the actual processing. To investigate the compatibility between graphene and rubber at different temperatures, the equilibrated structures after 2 ns MD simulation were annealed. The annealing procedure with five annealing cycles and 25 heating ramps per cycle from an initial temperature of 273 K to a midcycle temperature of 500 K was then performed. The dynamics steps per ramp are 10⁶ and the total number of steps is 5 × 10⁷. The annealing procedure was carried out with NVE (constant number of atoms, volume, and energy) ensemble with a time step of 1 fs. The last annealing cycle was used to calculate E_mix from 273 K to 500 K. The E_mix of graphene/rubber mixtures as a function of temperature at optimum grafting ratio (see Table 5) is shown in Fig. 9.

Fig. 9 E_mix as a function of temperature for functionalized graphene/rubbers mixtures at optimum grafting ratios: (a) NR, (b) BR, (c) SBR, (d) NBR, and (e) CR.

According to the Flory–Huggins lattice model, the lower the E_mix, the better the compatibility. It can be seen that there are three types of E_mix change with temperature: (1) E_mix decreases continuously with the increase in temperature. That is, as the temperature increases, the compatibility of the two components becomes better. (2) E_mix increases first and then decreases with the increase in temperature. That is, the compatibility first gets worse and then better with temperature. (3) With the increase in temperature, E_mix first decreases, then rises, and finally decreases. That is, the compatibility of the two components first becomes better, then becomes worse, and finally becomes better with temperature. To explain the above three phenomena, we select three representative systems for structural analysis, namely, (1) 8.7% 2-mercaptoethylamine/NR (the first type), (2) 8.7% 2-mercaptoethylamine/SBR (the second type) and (3) 23.2% 2-mercaptoethylamine/NBR (the third type). We calculate the variation of the fractional free volume (FFV) with temperature for the three systems. The FFV is defined as follows:

where V_f is the free volume and V_t is the total volume of the system. The FFV can characterize the stacking ability of molecular chains, and the smaller the FFV, the closer the molecular chains are stacked. As shown in Fig. 10, the FFV and E_mix of the three systems have a good agreement with temperature. Therefore, we speculate that the FFV is the essential cause of E_mix variation. For the 8.7% 2-mercaptoethylamine/NR mixture, the FFV decreases as the temperature increases. This is because the increase in temperature improves the fluidity of the polymer, and the rubber chains are better able to dissolve the graphene so that the graphene and rubber are packed more closely and the compatibility is improved. For the 8.7% 2-mercaptoethylamine/SBR mixture, the increase in temperature may lead to a decrease in the interactions between rubber and graphene, allowing the rubber and graphene stack to become loose, and thereby reducing FFV. As the temperature continues to increase, the fluidity of the rubber chains increases, and the rubber molecules are better able to dissolve graphene, so that graphene and rubber are packed more closely. For 23.2% 2-mercaptoethylamine/NBR mixture, the initial increase in temperature is conducive to tight packing of graphene and rubber. Because NBR is a polar rubber, graphene and NBR will form hydrogen bonds. Therefore, when the temperature continues to rise, the hydrogen bonds will be dissociated, resulting in a decrease in the interactions between graphene and rubber. When the temperature exceeds a certain value, the increase in the fluidity of rubber molecules is conducive to the dissolution of graphene, so that rubber and graphene are tightly packed. In summary, the competing effects between the fluidity of rubber chains and interactions cause E_mix to exhibit three different changes with temperature. The above conclusions can be used as a theoretical basis for regulating temperature to improve the compatibility between functionalized graphene and rubber.

Fig. 10 Fractional free volume (FFV) as a function of temperature for (a) 8.7% 2-mercaptoethylamine/NR mixture, (b) 8.7% 2-mercaptoethylamine/SBR mixture and (c) 23.2% 2-mercaptoethylamine/NBR mixture.

4 Conclusions

In this work, Hildebrand and two-dimensional solubility parameters of six functionalized graphene as a function of grafting ratio were calculated by MD simulation. The quantitative relationship between the surface chemistry of functionalized graphene and compatibility with rubber was constructed based on the Flory–Huggins lattice model. This work will provide theoretical insights at the molecular level for the functionalization of graphene and the preparation of high-performance graphene/rubber composites. The main conclusions are as follows:

(1) Three-dimensional Hansen solubility parameters are converted to two-dimensional solubility parameters based on the force field method. If the functional group (aniline and 1-octadedecylamine) is at the edge of graphene, the δ_T decreases with the increase of the grafting ratio, whereas if the functional group (phenyl isocyanate, 1-hexyl isocyanate, 3-mercaptopropionic acid and 2-mercaptoethylamine) is in plane of graphene, the δ_T decreases first and then increase with grafting ratio. The change of δ_T and δ_vdW for graphene is consistent with the grafting ratio. The δ_ele increases with the increase of grafting ratio.

(2) Two-component solubility parameters can predict the compatibility of functionalized graphene and rubber with different polarities well. The functionalized graphene and rubber mixtures exhibit two thermodynamically compatible behaviors with grafting ratio (V-shape and W-shape). The grafting ratio corresponding to the best compatibility between functionalized graphene and rubber is given.

(3) According to the functionalization principle of graphene, the priority of functional groups for NR, BR, and SBR is 1-octadedecylamine > 1-hexyl isocyanate > 3-mercaptopropionic acid > 2-mercaptoethylamine > aniline > phenyl isocyanate. The priority of functional groups for NBR is 3-mercaptopropionic acid > 1-octadedecylamine > 1-hexyl isocyanate > 2-mercaptoethylamine > aniline > phenyl isocyanate. The priority of functional groups for CR is 3-mercaptopropionic acid > 1-octadedecylamine > 1-hexyl isocyanate > 2-mercaptoethylamine > phenyl isocyanate > aniline. The compatible behavior of different graphene/rubber mixtures as a function of temperature at optimum grafting ratios is given.

Our study aims to theoretically provide the quantitative relationship between different functional groups/grafting ratios and solubility parameters through MD simulation. The experimental verification is one direction we want to work on. It is hoped that the effectiveness of our simulation can be verified after the experimental preparation and characterization methods are more mature.

Data availability

The authors confirm that the data supporting the findings of this study are available within the article [and/or] as its ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 51903122).

Notes and references

1. A. K. Geim and K. S. Novoselov, Nat. Mater., 2007, 6, 183–191.
2. A. Murugesan, J. L. Saraswathy and R. Chandran, Int. Polym. Process., 2022, 37, 505–522.
3. J. Wang, S. Li, L. Yang, B. Liu, S. Xie, R. Qi, Y. Zhan and H. Xia, Molecules, 2024, 29, 1009.
4. Y. L. Leong, H. N. Lim and I. Ibrahim, Mol. Syst. Des. Eng., 2023, 8, 1229–1251.
5. A. M. Adnan and J. Wang, Constr. Build. Mater., 2023, 402, 133038.
6. W. Yu, L. Sisi, Y. Haiyan and L. Jie, RSC Adv., 2020, 10, 15328–15345.
7. V. Georgakilas, J. N. Tiwari, K. C. Kemp, J. A. Perman, A. B. Bourlinos, K. S. Kim and R. Zboril, Chem. Rev., 2016, 116, 5464–5519.
8. S. Eigler and A. Hirsch, Angew. Chem., Int. Ed., 2014, 53, 7720–7738.
9. S. Rani, M. Kumar, R. Kumar, D. Kumar, S. Sharma and G. Singh, Mater. Res. Bull., 2014, 60, 143–149.
10. C. Xu, X. Wang, J. Wang, H. Hu and L. Wan, Chem. Phys. Lett., 2010, 498, 162–167.
11. K. Zhao, Y. Wang, W. Wang and D. Yu, J. Mater. Sci., 2018, 53, 14262–14273.
12. R. Salvio, S. Krabbenborg, W. J. M. Naber, A. H. Velders, D. N. Reinhoudt and W. G. van der Wiel, Chem.–Eur. J., 2009, 15, 8235–8240.
13. A. B. Bourlinos, D. Gournis, D. Petridis, T. Szabó, A. Szeri and I. Dékány, Langmuir, 2003, 19, 6050–6055.
14. Z. Xu, J. Zhang, M. Shan, Y. Li, B. Li, J. Niu, B. Zhou and X. Qian, J. Membr. Sci., 2014, 458, 1–13.
15. A. Mahmood, Y. Sandali and J.-L. Wang, Phys. Chem. Chem. Phys., 2023, 25, 10417–10426.
16. A. F. Barton, Chem. Rev., 1975, 75, 731–753.
17. C. M. Hansen, Hansen Solubility Parameters: A User's Handbook, CRC Press, 2007.
18. X. Fang, B. Li, J. C. Sokolov, M. H. Rafailovich and D. Gewaily, Appl. Phys. Lett., 2005, 87, 094103.
19. Y. Hernandez, M. Lotya, D. Rickard, S. D. Bergin and J. N. Coleman, Langmuir, 2010, 26, 3208–3213.
20. S. D. Bergin, Z. Sun, D. Rickard, P. V. Streich, J. P. Hamilton and J. N. Coleman, ACS Nano, 2009, 3, 2340–2350.
21. S. Ata, T. Mizuno, A. Nishizawa, C. Subramaniam, D. N. Futaba and K. Hata, Sci. Rep., 2014, 4, 1–8.
22. C. M. Hansen and A. L. Smith, Carbon, 2004, 42, 1591–1597.
23. T. Dooher and D. Dixon, Polym. Compos., 2011, 32, 1895–1903.
24. W. Zhu, J. Ma, X. Nan, P. O. Lartey and Y. Yang, Colloid Polym. Sci., 2019, 297, 213–224.
25. Y. Luo, R. Wang, W. Wang, L. Zhang and S. Wu, J. Phys. Chem. C, 2017, 121, 10163–10173.
26. J. Ma, J. Liu, W. Zhu and W. Qin, Colloids Surf., A, 2018, 538, 79–85.
27. P. May, U. Khan, J. M. Hughes and J. N. Coleman, J. Phys. Chem. C, 2012, 116, 24390–24391.
28. J. Gupta, C. Nunes, S. Vyas and S. Jonnalagadda, J. Phys. Chem. B, 2011, 115, 2014–2023.
29. Y. Luo, X. Chen, S. Wu, S. Cao, Z. Luo and Y. Shi, Langmuir, 2020, 36, 9291–9305.
30. M. M. Hayder, T. M. Moumita, S. Chowdhury and K. A. Rahman, Mol. Simulat., 2024, 1–12.
31. Y. Keramati, R. Ansari, S. Haghighi and M. Eghbalian, Mol. Simulat., 2024, 1–13.
32. X. Huang, C. Zhi, Y. Lin, H. Bao, G. Wu, P. Jiang and Y.-W. Mai, Mater. Sci. Eng. R Rep., 2020, 142, 100577.
33. G. G. Vogiatzis and D. N. Theodorou, Arch. Comput. Methods Eng., 2018, 25, 591–645.
34. K. Lee, H. J. Lim, S. J. Yang, Y. S. Kim and C. R. Park, RSC Adv., 2013, 3, 4814.
35. H. Moradi, H. Azizpour, P. Khadiv-Parsi and R. G. Nia, Chem. Eng. Technol., 2023, 46, 2167–2174.
36. S. Chen, H. Yang, K. Huang, X. Ge, H. Yao, J. Tang, J. Ren, S. Ren and Y. Ma, Polymers, 2021, 13, 3778.
37. X. Chen, C. Yuan, C. K. Y. Wong and G. Zhang, J. Mol. Model., 2012, 18, 2333–2341.
38. Q. Chen, Z. Zhang, W. Huang, J. Qu, Q. Zhang, X. Wu, L. Zhang and J. Liu, Polym. Int., 2024 DOI:10.1002/pi.6647.
39. A. V. Savin and M. A. Mazo, Phys. E, 2020, 118, 113937.
40. S. Hoyas, V. Lemaur, Q. Duez, F. Saintmont, E. Halin, J. De Winter, P. Gerbaux and J. Cornil, Adv. Theory Simul., 2018, 1, 201800089.
41. G. Zhao, X. Li, M. Huang, Z. Zhen, Y. Zhong, Q. Chen, X. Zhao, Y. He, R. Hu, T. Yang, R. Zhang, C. Li, J. Kong, J.-B. Xu, R. S. Ruoff and H. Zhu, Chem. Soc. Rev., 2017, 46, 7469.
42. M. Song, M. Wang, C. Wang, J. Song, Y. Li, F. Cao, G. Yu and Q. Qin, Macromol. Theory Simul., 2023, 32, 202200072.
43. K. Luo, G. You, S. Zhang, W. Zheng and S. Wu, Polymer, 2020, 188, 122143.
44. S. Thakral and N. K. Thakral, J. Pharm. Sci., 2013, 102, 2254–2263.
45. S. Dhamankar and M. A. Webb, ACS Macro Lett., 2024, 13, 818–825.
46. D. W. Boukhvalov and M. I. Katsnelson, J. Am. Chem. Soc., 2008, 130, 10697–10701.
47. J. Chen, J. Li, J. Zhu, X. Chen, Y. Liu, Y. Gao, Z. Luo, S. Wu and Y. Luo, J. Phys. Chem. B, 2021, 125, 5420–5433.
48. Y. Luo, X. Chen, H. Liu, H. Zhang, M. Song, J. Liu and Z. Luo, Comput. Mater. Sci., 2021, 186, 8663.
49. Z. Tang, L. Zhang, W. Feng, B. Guo, F. Liu and D. Jia, Macromolecules, 2014, 47, 8663–8673.
50. X. Liu, L.-Y. Wang, L.-F. Zhao, H.-F. He, X.-Y. Shao, G.-B. Fang, Z.-G. Wan and R.-C. Zeng, Polym. Compos., 2018, 39, 1006–1022.
51. L. Zhang, Y.-K. Huang, L.-N. Yue, L. Xu, J.-Y. Qian and X.-D. He, Carbohydr. Polym., 2022, 296, 119951.
52. K. Shitajima, N. Karyu, K. Takada, S. Fujii, Y. Nakamura and Y. Urahama, Polym. Polym. Compos., 2015, 23, 121–128.

---

Footnote

† Electronic supplementary information (ESI) available: Validation of models and methods. R and χ of functionalized graphene/rubbers mixtures. See DOI: https://doi.org/10.1039/d4ra08062a

---