## Abstract

In this paper, we demonstrate the usage of the Nernst–Planck equation in conjunction with mean-field theory to investigate particle-laden flow inside nanomaterials. Most theoretical studies in molecular encapsulation at the nanoscale do not take into account any macroscopic flow fields that are crucial in squeezing molecules into nanostructures. Here, a multi-scale idea is used to address this issue. The macroscopic transport of gas is described by the Nernst–Planck equation, whereas molecular interactions between gases and between the gas and the host material are described using a combination of molecular dynamics simulation and mean-field theory. In particular, we investigate flow-driven hydrogen storage inside doubly layered graphene sheets and graphene-oxide frameworks (GOFs). At room temperature and with slow velocity fields, we find that a single molecular layer is formed almost instantaneously on the inner surface of the graphene sheets, while molecular ligands between GOFs induce multi-layers. For higher velocities, multi-layers are also formed between graphene. For even larger velocities, the cavity of graphene is filled entirely with hydrogen, whereas for GOFs there exist two voids inside each periodic unit. The flow-driven hydrogen storage inside GOFs with various ligand densities is also investigated.

## 1. Introduction

Humans are facing an unprecedented environmental crisis in the twenty-first century arising from inhalation of tiny toxic particles. Moreover, the discharge of carbon dioxide into the atmosphere has significantly contributed to global warming. To tackle these problems, a gas separation technique could be used to filter these toxic and harmful gases before they are inhaled by individuals or are discharged into the atmosphere. Also, gases such as hydrogen, methane and natural gas could be stored inside fuel cells to generate clean energy with low or zero carbon emission.

Recent advances in nanoscience and nanotechnology have provided new nanomaterials that have profound implications for science and technology [1]. Typical examples of nanomaterials include but are not limited to nanotubes [2], nanorods [3], nanowires [4], nanoparticles [5], graphene [6], zeolite [7], metal–organic frameworks [8] and zeolitic imidazolate frameworks [8]. Some of the existing and prospective applications include the use of TiO_{2} nanoparticles to enhance the efficiency of solar cells [5,9], nanomaterials as drug delivery carriers [10–12], gas storage inside nanomaterials [13–15], ultrafiltration and seawater desalination using nanomaterials [16–20], nanomaterials as effective cathode materials for Li–ion batteries [21–24], the usage of nanomaterials for fuel storage [25–27–29], nanomotors and rotors for energy harvesting in micromachines [30], and enhancing the mechanical strength of host materials by incorporating appropriate nanomaterials into composites [31,32]. In particular, several nanomaterials such as graphene, metal–organic frameworks and zeolitic imidazolate frameworks have been proved to be better candidates for storing gases than conventional materials [25,28,33].

Numerous theoretical studies have been performed to investigate the mechanism of gaseous storage using nanomechanics. For instance, Dimitrakakis *et al.* [26] used *ab initio* and grand canonical Monte Carlo calculations to show that the maximum hydrogen uptake inside graphene is 2 wt.% at 77 K and 1 bar, which matches well with the experimental results of 1.2 and 1.7 wt.% given by Srinivas *et al.* [34] and Ghosh *et al.* [35], respectively. Thornton *et al.* [28] adopted a mathematical model to show that magnesium-decorated metal–organic frameworks accommodate more hydrogen than normal metal–organic frameworks. This might be explained by the extra molecular interactions arising from the decoration that occupies some of the physical volume. Dimitrakakis *et al.* [26] adopted molecular dynamics (MD) simulations to show that an envisaged three-dimensional-pillared graphene possesses the remarkable hydrogen uptake of 9 wt.% at 77 K and 1 bar. Adisa *et al.* [14] showed that the methane configuration inside carbon nanotubes depends on the nanotube radius. Chan & Hill [25] also showed that the molecular ligands between graphene sheets not only provide mechanical support for the graphene-oxide frameworks (GOFs) but also induce extra molecular forces to enhance hydrogen storage.

As a result of the technological importance in gas separation and fuel cells, there has been an enormous amount of theoretical work, but, to the best of our knowledge, none of these investigations takes into account the flow field which is the predominant mechanism to squeeze gas into nanomaterials. Unfortunately, it is numerically infeasible to incorporate the effects of flow into MD simulations for nanomodels, because the flow field acts at a macroscopic scale. Some studies have endeavoured to incorporate the diffusivity. However, only diffusion driven by the concentration and temperature gradient and the intermolecular interactions between gases and the host material have been considered [21,36]. Several particle-laden flows inside microchannels have been investigated. For instance, Yu *et al.* [37] investigated the transport of nanoparticles in microchannels mimicking capillaries, and Heim *et al.* [38] used a microreactor to mix a particle-laden gas stream with a clean gas stream. In such situations, understanding the flow-driven transport inside nanomaterials becomes crucial to optimize material selection and design for gaseous storage.

In this paper, we aim to address such problems using the principle of multi-scale modelling [39]. While the flow-driven transport of concentration inside nanomaterials is modelled by a macroscopic conservation equation, the interactions between the gas and the host nanomaterial are determined by a MD simulation along with mean-field theory. In previous work, we have successfully employed mean-field theory to investigate ultrafiltration using carbon nanotube membranes [16–18,40–43], molecular, energy and charge storage inside nanomaterials [25,44–46], the axial buckling of nanopeapods [47], nanoshuttle memory devices [48], hydrogen yield using functionalized nanomaterials [49] and several intriguing mechanisms occurring at the nanoscale [30,50]. An example of the hydrogen uptake inside graphene and GOFs will be investigated and serves as a demonstration for the present hybrid methodology. While graphene could be viewed as a two-dimensional sheet comprising carbon atoms that are arranged in a hexagonal lattice structure [6], graphene-oxide frameworks (GOFs) are made from graphene oxide through a chemical reaction involving the boronic acids, hydroxy groups [51] and benzenediboronic acid pillars, which we will refer to as ligands.

In §2, we introduce the theoretical and numerical background for the rest of the paper, and the corresponding numerical results are given in §3, where some basic results are presented in §3a followed by more detailed results in §3b. The conclusions are presented in the final section of the paper.

## 2. Theory and molecular dynamics simulations

Here, a general theory for the present topic is provided, where the evolution of gas driven by a flow with velocity *U* inside a nanomaterial is investigated (figure 1). We note that bold letters represent vectors throughout this paper. The movement of gas inside the nanomaterial is subject to diffusion across concentration gradients, advection by the flow field and the molecular interactions between molecules and between molecules and the host material. The conservation of mass inside the host material is reflected by the Nernst–Planck equation [52]
*c*, *t*, *D*, *k*_{B}, *T* and *D* and *T* are constant and that the flow is incompressible and only in the *z*-direction, equation (2.1) reduces to
*C*, ℓ, *x*_{i} and *m* denote the initial concentration, the material's length in the *z*-direction, {*x*,*y*,*z*} coordinates and the mass of a single molecule, respectively. Using this non-dimensionalization, equation (2.2) becomes
*Pe*=*U*ℓ/*D* denotes the usual Peclet number. We observe that, while the diffusion mechanism dominates when the flow field is low or the diffusion constant is large, the third term on the right-hand side in equation (2.4) dominates for large flow fields, strong molecular forces or low temperatures. Equation (2.4) could easily be solved using any standard numerical methods such as finite difference or finite-element methods.

The force field *et al.* [53,54] adopted a mean-field theory approach to approximate the pairwise interactions between molecules by certain surface integrals. Suppose that *V* (*ρ*_{ij}) denotes the non-bonded interaction between two molecules, where *ρ*_{ij} is the atomic distance between the *i*th and *j*th molecules, then the total pairwise interactions, *E*, between two molecules with surface areas *S*_{1} and *S*_{2} can be approximated by [53,54]
*η*_{1} and *η*_{2} denote the number density of atoms on the surface *S*_{1} and *S*_{2}, respectively. Given *E*, the force field between two molecules can be obtained by computing

## 3. Numerical result and discussion

In this section, the particle-laden flow of hydrogen molecules inside both the doubly layered graphene and GOFs is investigated, where the coordinate system adopted in this paper and the direction of flow are shown in figure 1. Because both nanomaterials admit the same periodic structure when taking ligands into consideration, we investigate only a repeated unit. It is worth noting that the results deduced here could be generalized to other gases and nanomaterials without any conceptual difficulties. A graph of the GOFs, where layers of graphene are connected and separated by benzenediboronic acid pillars, each comprising eight hydrogen, four oxygen, two boron and six carbon atoms [55], is shown in figure 1.

Referring to figure 1, the flow direction is parallel to the graphene sheets and is defined to be the *z*-direction. The *x*-axis is defined to be parallel to the graphene sheets but perpendicular to the flow direction, whereas the *y*-axis is defined to be perpendicular to the graphene sheets. We assume that hydrogen can be described as point masses and we smear carbon atoms on the graphene and hydrogen, oxygen, boron and carbon atoms on the ligands. We model the interaction between molecules by the 6–12 Lennard–Jones potential [56] and mean-field theory [53,54]. Under these assumptions, equation (2.5) gives the total force acting on a hydrogen molecule inside the cavity [25]
*i*, *η*, *A*, *B*, *S*, *L* denote the *i*th element on the ligands, the number density of carbon on the graphene, the attractive constant for hydrogen and carbon, the repulsive constant for hydrogen and carbon, the maximum base radius of the repeated unit and the separation between two graphene sheets, respectively. In addition, *A*_{i}, *B*_{i}, *η*_{i} and *r* denote the attractive constant for hydrogen and the *i*th element on the ligands, the repulsive constant for hydrogen and the *i*th element on the ligands, the number density for the *i*th element on the ligands and the radial distance in the plane of the graphene sheet given by *B* and *y*-direction, whereas the ligands generate forces in both the *x*- and *z*-directions. Numerical values for all the constants adopted in this paper are provided in table 1.

### (a) Basic results

Here, we investigate hydrogen uptake driven by a flow of strength 0.1 ms^{−1} inside a repeated unit of doubly layered graphene sheet and GOF120, where the base radius is *S*=9.84 Å and the separation of graphene is given by *L*=11 Å . Here, the number 120 indicates that there are 120 carbon atoms on each graphene sheet for a given repeated unit. Inspired by the form of the solution of the Navier–Stokes equations [59], we simply assume laminar flow in the *z*-direction, which could be written as
*U*, *X* and *Y* denote the maximum speed, the total length in the *x*-direction and the total length in the *y*-direction, respectively, and a numerical result is given in figure 2.

Owing to the confinement of the porous nanomaterials, the flow distribution follows a well-defined laminar flow in which the velocity strength declines to zero on the boundaries. We adopt the flow distribution given in figure 2 and substitute both the flow distribution and the force fields as given in equation (3.1) into equation (2.4) to determine the concentration evolution for both nanostructures, i.e. graphene and GOF120. We assume periodic boundary conditions and use a standard finite-difference approach. We further assume the initial concentration of 0.1 Å ^{−3} and, without loss of generality, the saturation concentration is assumed to be 1 Å ^{−3}, where equation (3.2) reaches its maximum due to the fact that the periodic boundary is assumed. Because the solution in the *x*, *z* plane is symmetric, only the results arising from the top-right corner (quarter) for both host materials are shown. The temperature is set to be 300 K. The numerical result of the equilibrium concentration inside the graphene in the (*x*, *z*) plane at the outlet is presented in figure 3.

Our numerical results indicate that hydrogen molecules diffuse rapidly towards the boundaries of the graphene and form a single layer almost instantaneously after 1.31 ns near the graphene surface leaving the central cavity of the graphene almost completely void. Owing to the repulsive force induced between the layer and the host material, this single layer oscillates near the graphene surface even if it reaches the steady state. To investigate the instability of the system under strong flow fields, using equation (2.4), and upon assuming equilibrium conditions and constant velocity *u*=*U*, we obtain
*y*. Because a small oscillation is observed in figure 3, this is feasible if and only if *E*^{2}−4*F*<0, from which we deduce the maximum velocity field for which the single layer is stable, which is given by

Now, the particle-laden flow inside GOF120 is investigated. Again, an external temperature of 300 K is assumed, and the maximum inflow velocity of 0.1 ms^{−1} is used; the numerical results with and without the consideration of H–H interactions are given in figure 4 for comparison. We note that, in contrast to the single layer induced by the graphene sheet, multi-layers are formed when the extra ligands are incorporated. Our numerical results suggest that the single layer is initially formed near the graphene surface, whereas the central cavity remains approximately void. Soon after that, the first layer emerges in the vicinity of the centre, which is parallel to the ligands. This layer is subsequently attracted towards the ligands. After the first layer reaches the proximity of the ligand, a second layer begins to form and then moves towards the first one. However, when the second layer approaches the first layer, they do not merge with each other but are separated by a certain distance. This is the case both with and without H–H interactions. Three layers are eventually formed as a result of the current numerical settings without H–H interactions, and hence the ligands not only support the GOFs but also induce stability to the original graphene system, resulting in more hydrogen uptake for extremely low-flow speed. To understand the effects of H–H interactions, equation (3.2) is considered, and we plot the results with and without H–H interactions in figure 4. It is worth noting that only two layers result in this case due to the fact that the H–H interactions maintain the correct distance between layers.

We now consider the same system when the inflow velocity magnitude is increased to 1 ms^{−1}. To investigate the corresponding changes in the concentration pattern for graphene and GOF120 both with and without H–H interactions, the numerical results are given in figure 5. Under this low velocity and for the case of graphene, multi-layers are formed in parallel with graphene with H–H interactions, where the H–H interactions between layers generate repulsive forces to prevent the merging of different layers, as happens without H–H interactions (figure 5*a,¡textit¿b*¡/textit¿). For the case of GOF120, without H–H interactions, a pattern of a single layer and multi-layers is formed in parallel with the graphene and the ligands, respectively. However, when H–H interactions are included, a dramatically different concentration pattern forms, where the resultant pattern is the superposition of the multi-layers arising from both graphene and ligands. In conjunction with the H–H interactions, a weaker layer is also formed in the proximity of the ligands (figure 5*d*).

We now consider the same system when the flow is further increased to 50 ms^{−1} to scrutinize the storage pattern under high-flow fields. The numerical results are given in figure 6. For graphene sheets, because the flow carries more hydrogen into the system, the cavity is filled almost entirely with hydrogen except in the vicinity of the graphene where the strong repulsive forces will push hydrogen away from the sheet. We note that higher velocity fields are needed to reach this configuration using the laminar flow in equation (3.3) than with a uniform flow for which the configuration can be achieved around 20 ms^{−1}. Therefore, the distribution of macroscopic flow becomes more significant for the storage pattern in the high-velocity limit. For GOF120, the central cavity is filled almost entirely, as in the case of graphene. However, as we approach the ligands, a peculiar thing happens. There exists a region where hydrogen is free to pass through. To illustrate this effect, H–H interactions are plotted using white arrows in conjunction with the hydrogen concentration in figure 6*b* for comparison. In the red region, H–H interactions exist and act against the molecular forces, so that the hydrogen reaches the equilibrium state, and the closer the hydrogen to the graphene, the larger these H–H interactions. However, in the void region, the H–H interaction is negligible, which implies that no molecular force exists in this area, resulting from the zero net interaction generated by the graphene and ligands.

### (b) Further investigations

In this section, we further investigate the effect of temperature and inflow velocity on the hydrogen concentration inside graphene and various GOFs based on the basic results given in the previous section. Here, simulations are run until the system reaches its equilibrium state. In figure 7, the average concentration is plotted against velocity for GOF120 and graphene for a temperature of 300 K. The inflow velocity is varied from 0.1 to 2 ms^{−1} with an increment of 0.1 ms^{−1} and then from 2 to 20 ms^{−1} with an increment of 1 ms^{−1}.

GOF120 admits a higher concentration than graphene in the extremely low-velocity limit, which could be partially explained by the fact that a single layer is formed parallel to the graphene, but multi-layers are induced by the ligands inside GOF120 (see also figures 3 and 4). However, for the higher velocity fields, graphene admits higher concentrations. There exists a jump around *U*=0.8 ms^{−1} which is caused by the sudden creation of the multi-layers induced by the graphene, and the concentration for both materials drops slightly afterwards due to the fact that certain parts of the layers are being destroyed by the higher flow field. When *U*>6 ms^{−1}, as the flow field squeezes more and more hydrogen into the systems, a higher concentration results in both cases. This concentration is approximately saturated when *U* is larger than 20 ms^{−1}, and the gap between the concentrations for both cases is caused by the void region in GOF120 (figure 6). It is worth noting that the full concentration, i.e. 1 Å ^{−3}, is never reached owing to the strong repulsive force existing in the vicinity of graphene. Now, three velocities—namely 0.1, 10 and 20 ms^{−1}—are fixed, and the temperature effect is investigated, which is shown in figure 8. For both materials, the change in temperature causes no effect on the concentration at *U*=0.1 ms^{−1} owing to the tight interactions between hydrogen and the host materials. On the other hand, the temperature generates only a small decay in the concentration when the cavity is almost saturated with hydrogen. It turns out that temperature changes cause more dramatic changes in the intermediate velocity regime. In this regime, the concentration for graphene decreases more rapidly than for GOF120 as the ligands inside GOF120 help to stabilize hydrogen storage.

Now, we investigate both the temperature and velocity effects on the graphene and various GOFs, including GOF120, GOF66 and GOF28, where the number indicates the total number of carbon atoms in each graphene sheet. Hence, the smaller the number, the higher the ligand density. Any GOFs with a ligand density higher than GOF28 will be deemed too dense for hydrogen storage. First, we fix the temperature at 300 K to investigate the velocity effect, then we fix the velocity at 10 ms^{−1}, which reveals the largest variation as shown above to investigate the temperature effect. The numerical results are given in figure 9. As can be seen in figure 9*a*, the concentration for all proposed GOFs is higher than that for graphene in the extremely low-velocity regime, the jump occurs for all GOFs and graphene around *U*=0.8 ms^{−1} owing to the formation of multi-layers by graphene, and the tiny drop in the concentration results for all cases when *U*≈1 ms^{−1}. These results are consistent with the numerical results given in figure 7. However, there exist significant differences in the concentration in the intermediate velocity regime, where the concentration of GOF28 outperforms that of other GOFs and graphene, which is in good agreement with another study conducted by Chan & Hill [25] without the consideration of any flow fields. However, when *U*>11 ms^{−1}, as more hydrogen is squeezed into these materials, the geometric effect becomes important and the graphene turns out to have the highest concentration. For the fixed velocity field as indicated by figure 9*b*, graphene possesses the highest concentration when *T*<220 K; however, its concentration decays more rapidly at higher temperatures. The decay of the concentration for GOF120 follows a similar pattern to that of graphene but at a slower rate. The concentration for GOF66 and GOF28 does not decay as fast as that of graphene and GOF120 owing to the extra stability generated by the denser ligands. When *T*>220 K, GOF28 has a higher concentration than graphene and the other GOFs owing to its appropriate ligand density.

In conclusion, the ligands in GOF28 not only support the graphene structure but also give an appropriate ligand density to encapsulate more hydrogen under moderate temperatures and velocity fields, and the graphene reveals better hydrogen storage capacity in the high-velocity limits where its cavity is completely filled.

## 4. Conclusion

In this paper, a multi-scale model has been proposed for solving the particle-laden flow inside nanomaterials, where the Nernst–Planck equation is used to model the concentration of particles and the force field between the particle and the host materials is captured by a MD simulation using mean-field theory. In particular, we demonstrate the validity of the model by investigating the concentration of hydrogen inside graphene and various GOFs. We find that graphene possesses the highest storage capacity in the high-velocity limit, and GOF28 possesses better storage capacity under moderate temperatures and flow fields. The present calculations are in good agreement with the existing literature, but the method is significantly more computationally efficient.

## Data accessibility

The sources of all data are provided in table 1.

## Authors' contributions

Y.C. carried out the theoretical study, participated in the numerical analysis and drafted the manuscript; J.J.W. carried out the theoretical study, provided some suggestions and proofread the article; L.X., Y.R. and Y-T.C. provided some suggestions and proofread the article.

## Competing interests

We have no competing interests.

## Funding

This research was supported by the Zhejiang Provincial Natural Science Foundation of China under grant no. Q15E090001, the Natural Science Foundation of China-Young Scientist Programme (NSFC51506103/E0605), the Qianjiang Talent Scheme (QJD1402009) and the Ningbo Natural Science Foundation (2014A610025), (2014A610172) and (2015A610281).

- Received April 26, 2016.
- Accepted July 15, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.