Continuum Modeling with Functional Lennard-Jones Parameters for Methane Storage inside Various Carbon Nanostructures.

Methane capture and storage are of particular importance for the development of new technology to reduce the effects of climate change and global warming. Carbon-based nanomaterials are among several porous nanomaterials proposed as potential candidates for methane storage. In this paper, we adopt a new continuum approach with functional Lennard-Jones parameters to provide interaction energies for methane inside carbon nanostructures, namely fullerenes, nanotube bundles, and nanocones. This study provides a significant improvement to previous continuum modeling approaches using the Lennard-Jones potential.

Introduction

Methane (CH4) is the smallest hydrocarbon molecule, comprising one carbon and four hydrogen atoms (Figure a), and is well-known as one of the main greenhouse gases in the atmosphere causing global warming and climate change.[1,2] Even though the concentration of carbon dioxide (CO2) is much higher than that of methane in air, the effect of methane on the rising earth temperature and global warming is found to be 20–28 times higher than that of CO2.[3] A 2021 report also shows the rapid rise of methane in the atmosphere, where methane concentration was found to have risen past 1800 ppb, tripling preindustrial levels.[4] Capturing methane and other greenhouse gases from the atmosphere is therefore an urgent and critical task which attracts research attention worldwide. However, capturing methane from the atmosphere is very challenging due to its dilution in air, its stable C–H bond that is hard to activate,[5] and the difficulty of separating it from nitrogen.[6] As a result, a number of techniques have been proposed to capture methane, which include membrane absorption,[7] adsorption to nanoporous materials,[8] and more advanced adsorption methods.[9−12]

Figure 1

Three representations of methane (CH4): (a) discrete stick-and-ball model, (b) semicontinuous (point–shell model), and (c) fully continuous heterogeneous sphere model.

Another related area of research is methane storage. This area has also received much research interest ranging from designing and synthesizing suitable molecular containers to understanding mechanisms of adsorption and interactions between methane and potential storage structures.[13−16] A number of porous nanomaterials have been proposed for methane storage including metal–organic frameworks, porous aromatic frameworks, covalent organic frameworks, zeolites, activated carbon, and other porous materials.[17−21] Further, carbon-based nanomaterials are also considered as potential candidates for methane storage due to their excellent mechanical strength and chemical and thermal stability.[22,23]

The adsorption of methane onto surfaces of carbon structures is typically governed by the physisorption mechanism.[24−26] As the forces present in physisorption are predominantly van der Waals forces, the interactions between methane molecules and carbon structures can be evaluated using the Lennard-Jones potential[27]which computes the potential energy, E, between two atoms at distance ρ apart. The coefficients A and B are the attractive and repulsive constants of the interaction, respectively. In the conventional discrete approach adopted in molecular dynamics studies, the potential energy between two molecules is simply the pairwise sum of the potentials of their constituent atoms. In order to reduce computational costs, as well as to produce mathematical expressions that are tractable and readily computable, continuous models are often adopted to replace the discrete atomic structures of the molecules.[28−35]

In the continuum approach, atoms on the interacting molecules are assumed to be distributed, typically uniformly, over the surfaces for shell-like molecules or the volumes for dense molecules. Thus, the interaction energy can be evaluated by computing integrals of the interatomic potential over the surfaces or volumes, scaled by their atomic densities.[29] While this approach is suitable for homogeneous structures comprising the same type of atoms (also known as homonuclear molecules), it seems to not work well for heterogeneous, or heteronuclear, molecules, such as methane and coronene, which comprise more than one type of atoms.[36,37] Specifically, coronene and methane comprise carbon atoms in the central region and hydrogen atoms at the perimeter. Smearing the entire surface (or volume) homogeneously implies the same interaction strength throughout the molecule, which seems counterintuitive as interactions involving hydrogen atoms are generally weaker than those involving carbon atoms. To address this shortcoming, a semicontinuous approach has been proposed that involves approximating a heteronculear molecule by multiple homonuclear pieces. For example, benzene or coronene uses a combination of concentric rings of carbon and hydrogen atoms, where the carbon or hydrogen atoms are uniformly distributed on each ring.[31,38] For methane, a semicontinuous approach can be thought of as a point–shell model, where the carbon atom is represented as a point centered inside a spherical shell over which the four hydrogen atoms are homogeneously smeared (see Figure b).

The point–shell approximation has been used to model the interactions of methane with various nanostructures, such as nanotubes[39] and nanotube bundles,[40] graphene,[41] nanocones,[42] and fullerenes,[41] in order to investigate their adsorption properties and their potential use as methane storage. While the semicontinuous model may give a reasonable approximation for the interaction energy for small molecules, such as coronene and methane, extending the model to larger size structures is computationally expensive. To address this issue, Stevens et al.[36] propose a new continuum model taking into account the effect of different types of atoms in heterogeneous structures. This model considers a molecule as one continuum structure for which the interaction strengths are varied as a continuous function depending on the location and type of the atoms comprising the molecule. For the structure of methane, we assume a dense sphere (see Figure c) such that the interaction strength is stronger in the central region (red) due to the carbon atom and weaker toward the spherical surface (blue) due to hydrogen atoms. Accordingly, interaction functions A(r) and B(r) are used in place of the interaction constants (A and B) in the Lennard-Jones potential.

For methane–nanotube and coronene–nanotube interactions, Stevens et al.[37] demonstrate that applying a homogeneous approximation to a heterogeneous molecule can lead to inaccurate results, particularly at small intermolecular distances. Furthermore, results from Stevens et al.’s[36] continuous model of methane encapsulated within a carbon nanotube show better agreement with molecular dynamics simulations than that of the semicontinuous approach used by Adisa et al.[39] As such, this paper applies the continuous model assuming methane as a whole sphere of varying interaction strength[36] to study the problems presented by Adisa et al.,[40−42] which are methane encapsulated in a fullerene, methane within a nanotube bundle, and methane entering an open nanocone. Our aim is to provide an alternative method for energy calculations that give better agreement to molecular dynamics simulations and is less computationally expensive compared to a semicontinuous model. We comment that it is important to continue the development of a continuum approach to produce more realistic results that can be used as a guideline for experiments prior to investing in highly computational methods, such as ab initio and density functional theory.

In section , we describe the whole sphere continuum model for methane and state analytical expressions for the potential energy of methane interacting with fullerenes, nanotube bundles, and nanocones. In section , numerical results are presented in comparison with those obtained from molecular dynamics simulations and semicontinuous models, which we reproduce based on Adisa et al.[40−42] Conclusions are given in section , and Appendix A, Appendix B, and Appendix C provide detailed derivations for the integrals K for the three problems. Finally, in Appendix D, we outline basic details for the molecular dynamics simulations.

Modeling Approach

In this section, we model the interactions between methane and three carbon nanostructures, namely fullerenes, nanotube bundles, and open nanocones.

Since the three carbon nanostructures comprise only carbon atoms and have shell-like structures, we approximate them as continuum surfaces where carbon atoms are uniformly distributed over their entire surfaces. For the methane molecule, we model it as a sphere of radius a. As shown in Figure c, the sphere has an inner core region (red) of carbon, an outer layer (blue) of hydrogen, and a transitionary layer (purple) that is a blend of the two.

According to Stevens et al.,[37] the interaction energy E between methane and a homogeneous nanostructure can be evaluated bywhere η is the atomic density of a given molecule in the interaction and K (n = 3, 6) is the integral defined bywhere ρ is the Euclidean distance between two molecules, dS is the surface element of the carbon nanostructure, and dV is the volume element of the methane molecule, which has the standard Cartesian representation (ar sin ϕ cos θ, ar sin ϕ sin θ, ar cos ϕ), where r ∈ [0, 1], θ ∈ [0, 2π], and ϕ ∈ [0, π] are three variables parametrizing the sphere. The functions f for n = 3 and 6 are interaction functions representing the attractive and repulsive coefficients of the Lennard-Jones potential, respectively. The interaction functions are introduced as a smooth transition of interaction strength between the inner carbon core and the outer hydrogen region, which are treated as discrete entities in the point–shell model. Using the interaction functions has the twofold benefit of accounting for the heteronuclear nature of methane and maintaining a fully continuous approximation of the molecule, allowing for a single expression to determine the potential energy between methane and the carbon nanostructure. As discussed in Stevens et al.,[37] we assume that f3(r) = A(r) and f6(r) = B(r) and that they follow a sigmoidal profile, where a stronger interaction is found in the inner region, due to the carbon–carbon interaction with the nanostructure, and a weaker interaction is toward the surface of the molecule, due to the hydrogen–carbon interaction. Here, we assume that the sigmoidal function has the form α arctan(m(r0 – r)) + β, where α and β for A(r) and B(r) can be found, respectively, from the conditions A(0) = AC–C, A(a) = AC–H, B(0) = BC–C, and B(a) = BC–H, where the attractive and repulsive constants for carbon–carbon (C–C) and carbon–hydrogen (C–H) are given in Table . The parameters r0 and m affect the profile of the sigmoidal function such that they control the location of the inflection point and steepness of transition, respectively. Here we use r0 = 0.565 and m = 20.[36] Numerical values for other constants used in our model, the semicontinuous model, and molecular dynamics simulations are given in Table .

Table 1

Numerical Values of Constants Used in This Paper

parameter	value
AC–C	560.44 kcal/mol·Å6
BC–C	1121755.66 kcal/mol·Å12
AC–H	129.67 kcal/mol·Å6
BC–H	91727.95 kcal/mol·Å12
methane radius	1.07 Å
mean volume density of methane	0.974 Å-3
mean surface density of nanotube	0.382 Å-2
radius of (8, 8) nanotube	5.428 Å
mean surface density of nanocone	0.382 Å-2
mean surface density of C60	0.379 Å-2
radius of C60	3.55 Å
mean surface density of C240	0.377 Å-2
radius of C240	7.12 Å
mean surface density of C540	0.390 Å-2
radius of C540	10.5 Å

Methane–Fullerene Interaction

Here we consider the interaction between a spherical fullerene of radius b interacting with a methane molecule, which is assumed to be encapsulated inside the fullerene (see Figure ). The derivation of the interaction energy for nested spheres is similar to that of two separated spheres, as outlined in Adisa et al.[41] However, here we have an additional radial integral to perform over the methane’s radius. Assuming the methane molecule is centered at (0, 0, δ), where δ is the distance from the center of the fullerene located at the origin, the integral K can be evaluated as

Figure 2

Diagram of methane–fullerene interaction.

The detailed derivation of this expression can be found in Appendix A.

Methane–Nanotube Bundle Interaction

Here we consider a methane molecule interacting with a bundle of N carbon nanotubes, each of which has radius b and is aligned parallel to and equidistant from a common axis, which we refer to as the bundle axis (see Figure ). The perpendicular distance from the bundle axis to the axis of each constituent nanotube is termed the bundle radius, denoted here by R. Thus, a general point on the ith tube in the bundle is given bywhere θ ∈ [0, 2π] and z ∈ (−∞, ∞) are variables parametrizing the ith nanotube. Here, we prescribe the coordinates for the center of the methane molecule to be (ε cos ψ, ε sin ψ, 0), where ε is the perpendicular distance between the methane’s center and the bundle axis.

Figure 3

Diagram of methane–nanotube bundle interaction.

Following Adisa et al.,[40] the evaluation of the integral K giveswhere

Full derivation of this expression can be found in Appendix B.

Methane–Nanocone Interaction

Here we consider the interaction of a methane molecule with an open carbon nanocone with apex angle 2ω, on which the coordinates of a typical point are given by (z tan ω cos θ, z tan ω sin θ, z), where z ∈ [L1, L2] and θ ∈ [0, 2π] are variables parametrizing the nanocone. As illustrated in Figure , L1 and L2 are the distances from the origin along the z-axis indicating the locations of the smaller and larger openings of the nanocone, respectively. For a methane molecule located at (0, 0, Z), the integral K can be found to be a function of Z, namelywhere I(z) and a detailed derivation of K are given in Appendix C.

Figure 4

Diagram of methane–open nanocone interaction.

To evaluate K, we need to determine the location of the cone’s smaller opening (L1) along the z-axis and the radius at this opening. This location is critical since it is where the methane molecule is radially closest to the cone. Hence, we expect the highest interaction energy at L1. On the other hand, we expect a negligible interaction between methane and the cone at L2. This is due to the size of the methane molecule compared to the size of the opening at L2; thus we use L2 = ∞ to simplify the integral K.

Here we only consider the largest carbon nanocone with angle ω ≈ 56.4%, which has been flattened as shown in Figure . The open cone can be made by removing layers of carbon atoms from the apex located at the origin. This is equivalent to removing the atoms within the circle shown in Figure . For example, we can remove all atoms inside the innermost circle for a small opening cone or all atoms inside the second or third circle to have a larger opening of the cone.

Figure 5

Flat representation of single-pentagon nanocone.

The circumference of these circles is given by , where and a0 is the bond length between two carbon atoms in the nanocone, taken here to be 1.42 Å.

According to the geometry shown, L1 is given byWe call this the “geometric” L1. In this paper we also consider an “adjusted” L1, where we extend L1,geom toward the origin by a certain length to account for the loss of atomic contribution at the edge of the opening. The adjusted form of L1 is given bywhere ηg is the atomic density of graphene and k is the number of atoms found at L1,geom. Adisa et al.[42] present another method for determining the location of the cone’s opening using the number of hexagons from the apex, K, which is given as .

Results and Discussion

Here we present the interaction energy between a methane molecule and three types of carbon nanostructures, namely, (i) spherical fullerenes of various radii (C60, C240 and C540), (ii) carbon nanotube bundles comprising six (8, 8) nanotubes with varying bundle radii, and (iii) open carbon nanocones with varying sizes of the smaller opening. The results obtained from our approach presented here are compared with molecular dynamics (MD) simulations (representing a fully discrete model) and those of Adisa et al.,[39−42] which consider methane as a point of carbon atom at the center of a spherical shell of four smeared hydrogen atoms. In the following, we refer to our approach as “a whole sphere model” and that of Adisa et al.[39−42] as “a semi-continuous model”.

All numerical evaluations for the interaction energies are performed with MAPLE, and all MD simulations are performed with LAMMPS.[43] As shown in this section, for a heterogeneous structure, such as methane, our improved continuum model gives excellent agreements with MD studies compared with a classical continuum approach previously used by Adisa et al.[39−42]

(i) The potential energies for a methane molecule inside spherical fullerenes C60, C240, and C540 are shown together in Figure . For C60, the potential energy is positive for all positions of methane inside the fullerene, indicating a repulsive force between the two structures. This implies that methane cannot be encapsulated inside a C60 fullerene. For larger fullerenes, we can see that the methane molecule can be accommodated inside C240 and C540 fullerenes due to the negative energy profile. From Figure , we obtain the equilibrium position of the methane molecule inside a fullerene, which is where the energy is at the minimum. It is clear that, in C240 and C540, the methane prefers to be closer to the inner surfaces of the fullerenes rather than locating at the centers of the fullerenes. We comment that the energy profile for methane inside a fullerene is repulsive when methane is close to the inner surface of the fullerene. As shown in Figure , the energy rises steeply when the methane position is beyond an equilibrium position (δmin). We find that, within approximately 2 Å from the fullerene’s shell, the repulsion strongly dominates the interaction.

Figure 6

Comparison between energy profiles from three models for the methane–fullerene interaction with varying fullerene radii: (a) 3.55 Å (C60); (b) 7.12 Å (C240); (c) 10.5 Å (C540). ×, simulation; —, whole sphere; − –, point–shell.

In Table , we summarize the equilibrium positions (δmin) and the corresponding minimum energies (Emin) for the methane molecule inside C240 and C540 fullerenes using the three approaches.

Table 2

Comparison of Results between Whole Sphere and Semicontinuous Methane Models against Simulation Values, for Varying Fullerene Radii

fullerene radius (Å)	model	δmin (Å)	Emin (kcal mol–1)
7.12 (C240)	MD simulation	3.40	–5.598
	whole sphere	3.38	–5.258
	semicontinuous	3.39	–8.070
10.5 (C540)	MD simulation	6.90	–4.201
	whole sphere	6.85	–4.176
	semicontinuous	6.81	–6.289

As shown in Figure , for all cases, results from all three approaches behave similarly. Although it is clear from Table that our results from the whole sphere model give much better agreement with the simulation values compared to the semicontinuous approach. To confirm this, we benchmark energy results (Emin) against the MD simulations and find that the absolute errors for the whole sphere model are much smaller than those for the semicontinuous model for both C240 and C540 fullerenes. For the C240 fullerene, the absolute errors for the whole sphere and semicontinuous models are εwhole = 0.34 kcal mol–1 and εsemi = 2.472 kcal mol–1, respectively. For the C540 fullerene, the absolute errors are εwhole = 0.025 kcal mol–1 and εsemi = 2.088 kcal mol–1. While the whole sphere model clearly outperforms the semicontinuous approach in determining energy minima, it is worth noting that both models perform equally well in locating the equilibrium positions. From Table , we can see that the absolute errors for the equilibrium positions (δmin) of methane inside C240 and C540 fullerenes are within 0.1 Å from the MD results.

We note that some significant energy fluctuations observed in the MD results presented in Figure arise as the methane molecule approaches the fullerene due to effects of atomic surface roughness of the fullerene and methane molecules. This phenomenon is not seen in the continuous models since the molecules’ surface is smeared out.

(ii) The methane–nanotube bundle interaction is investigated using six (8, 8) nanotubes in a bundle, with the methane molecule centered within the bundle. Figure shows the energy profiles using the three approaches, and Table depicts from each approach the critical bundle radius (Rmin) and the corresponding minimum energy (Emin). Again, compared to the MD simulations, our whole sphere model outperforms the semicontinuous model. This is evident from Figure , which depicts excellent agreement between the whole sphere model and the MD results, whereas the semicontinuous model overestimates the potential energy for the entire interval of the bundle radius considered. Benchmarking against the MD results, we find from Table that the absolute errors for Emin for the whole sphere and semicontinuous models are εwhole = 0.411 kcal mol–1 and εsemi = 1.63 kcal mol–1, respectively. However, we again note that both models perform equally well in determining the critical bundle radius (Rmin).

Figure 7

Comparison of energy profiles between three models for methane–nanotube bundle interaction, assuming six (8, 8) nanotubes in the bundle. Here, we fix the center of methane molecule at the origin while allowing the bundle radius R to vary. ×, simulation; red —, whole sphere; blue —, point–shell.

Table 3

Comparison of Results between Whole Sphere and Semicontinuous Models against Simulation Values for Methane inside a Bundle Comprising Six (8, 8) Nanotubes

model	Rmin (Å)	Emin (kcal mol–1)
MD simulation	9.08	–12.808
whole sphere	9.065	–12.397
semicontinuous	9.030	–14.438

(iii) The methane–nanocone interaction is investigated using a nanocone of apex angle 112.9° and various fixed opening locations (L1) as determined in section . In Figure , we show the relationship between the potential energy and the location Z of methane along the z-axis for four different sizes of the cone’s opening at L1, ranging from the smallest opening (removing carbon atoms inside the first circle) to the opening that is made by removing all atoms inside the fourth circle in the geometric model (see Figure ). Further, in Figure , we plot the energy using the two values of L1, namely L1,adj and L1,geom, and compare the results with MD simulations. As shown in Figure , using L1 = L1,adj gives results that are in closer agreement with the MD simulations for all four sizes of the openings. In addition, we note from Figure a that the energy is positive (repulsive) at the opening of this cone, which implies that the methane molecule cannot enter the nanocone due to the size restriction at this opening. On the other hand, Figure b–d indicates the acceptance of methane inside the cones.

Figure 8

Energy profiles from the whole sphere model using two values of L1 (L1,adj and L1,geom). Results for all four open sizes are compared with MD simulations. ×, simulation; —, geometrically determined L1; − · –, adjusted L1.

Using L1,adj, we show in Figure that the whole sphere model gives results that are much closer to those from the MD studies than the semicontinuous model used in Adisa et al.,[42] especially for large sizes of the open cones. The summary of the critical location Zmin and its corresponding minimum energy is given in Table for all the models presented here. Using the energy results from MD simulation as a benchmark, we obtain the absolute errors for the whole sphere model (with adjusted L1) to be decreasing as the opening size increases: εwhole,2 = 0.247 kcal mol–1, εwhole,3 = 0.05 kcal mol–1, and εwhole,4 = 0.009 kcal mol–1. We note that at first the absolute error between the semicontinuous model and the MD simulation is slightly smaller than that of the whole sphere model and then becomes comparatively larger as the opening size increases: εsemi,2 = 0.164 kcal mol–1, εsemi,3 = 0.187 kcal mol–1, and εsemi,4 = 0.065 kcal mol–1.

Figure 9

Comparison between energy profiles from the three approaches for the methane–nanocone interaction for four sizes of the cone opening. We note that we use L1,adj for the whole sphere model, while for the semicontinuous model we adopt the value L1,Adisa that is proposed in Adisa et al.[42] ×, simulation; ···, point–shell; − · –, whole sphere.

Table 4

Comparison of Results between Whole Sphere and Semicontinuous Methane Models against Simulation Values, for Varying Sizes of the Cone Opening at L1

opening size	model	Zmin (Å)	Emin (kcal mol–1)
second circle	MD simulation	3.250	–3.134
	whole sphere (geometric)	3.484	–2.154
	whole sphere (adjusted)	3.424	–2.887
	semicontinuous	3.694	–3.298
third circle	MD simulation	4.500	–0.742
	whole sphere (geometric)	4.671	–0.625
	whole sphere (adjusted)	4.556	–0.692
	semicontinuous	4.211	–0.929
fourth circle	MD simulation	6.070	–0.227
	whole sphere (geometric)	6.118	–0.209
	whole sphere (adjusted)	5.939	–0.236
	semicontinuous	5.598	–0.292

We comment that the magnitudes of the potential energies for the methane–nanocone interactions (e.g., Figure b–d) appear to be much smaller compared to other interactions considered in this paper. This is due to the growing distance between the methane molecule and the nanocone as the methane molecule travels along the z-axis toward the cone’s larger opening. Since the van der Waals forces dominate at short distances, the attraction is greater when the molecules are closer, which is the case for methane inside a fullerene and a nanotube bundle.

Conclusions

Carbon nanostructures have potential use for methane capture and storage. In this paper, we use a recently proposed continuum model for heterogeneous molecules to determine the interaction energies between methane and three types of carbon nanostructures. Results obtained are in closer agreement with MD simulations compared to those calculated based on a semicontinuous model used by Adisa et al.[40−42] We comment that the difference between the fully continuous and semicontinuous models is due to the assumption made in the space between the central carbon atom and the spherical shell of hydrogen atoms, as shown in Figure b,c. The semicontinuous approach restricts the interaction effect of atoms to the point and the shell, while in the fully continuous model the interaction effect of the atoms in methane is distributed across the molecule. We believe that the assumption made for the fully continuous model is closer to reality as each atom’s or molecule’s electrons do not exist in a fixed position; hence the fully continuous model gives results that are in close agreement with MD simulations. Further, we found for each type of structure the critical distance, size, or configuration that optimizes the binding energy with methane. These findings are crucial for the development of structures for methane capture and storage. It is also clear that future modeling of gas storage involving the Lennard-Jones potential and a continuum approach should aim to include as much detail as possible within continuous structures in order to minimize deviation from the discrete case. The methods presented here can be further applied to other gases, such as CO2, or other types of storage structures.