Grand Canonical Monte Carlo Simulation of Nitrogen Adsorption in a Silica Aerogel Model

In this paper, the Diffusion Limited Cluster Aggregation (DLCA) method is employed to reconstruct the three-dimensional network of silica aerogel. Then, simulation of nitrogen adsorption at 77 K in silica aerogel is conducted by the Grand Canonical Monte Carlo (GCMC) method. To reduce the computational cost and guarantee accuracy, a continuous-discrete hybrid potential model, as well as an adsorbed layer thickness estimation method, is employed. Four different structures are generated to investigate impacts of specific surface area and porosity on adsorptive capacity. Good agreement with experimental results is found over a wide range of relative pressures, which proves the validity of the model. Specific surface area and porosity mainly affect nitrogen uptake under low pressure and high pressure, respectively.


Introduction
Silica aerogels are materials with high porosity, high specific surface area (SSA) and low density.They are promising in many applications, such as in heat insulators [1], drug delivery [2], effluent disposal [3,4] and fuel storage.These applications are based on their extraordinary properties of low thermal conductivity and good absorbability, which is caused by high porosity and nano-scale pores.
Three factors contribute to the thermal conductivity of silica aerogels: solid conductivity, gas transport and radiation.Liu [5] found, in simulation, that gas transport exerts an important effect on temperature response.Lee [6] observed enhanced self-diffusion of adsorbed methanol in silica aerogel in an experiment.Sun [7] found that for gases that could be adsorbed, such as N 2 and CH 4 , flux caused by adsorption is on the order of the direct flux or larger during molecular dynamic simulation of molecular permeation through nanoporous graphene membranes.Adsorption could affect the heat insulation performance.Concerning various applications of silica aerogel, exact determination of the adsorption amount would be meaningful.However, there are few simulation studies on the adsorption in silica aerogels currently.On the other hand, just as in simulations that explored the design principles of metal-organic framework (MOF) [8,9] for certain purposes, molecular simulation may also guide us in surface modification of silica aerogels to achieve better performance in the fields mentioned above.
The grand canonical Monte Carlo (GCMC) method has been proved to be an effective way to simulate the adsorption isotherm and it has been widely used in predicting the uptake capacity of various materials [10,11].Nitrogen adsorption has been a useful method in characterizing the nanoporous materials for pore size distribution [12,13] and SSA estimation [14] and there is plenty of experimental data for comparison.MacElroy [15] studied properties of methane, such as the Henry constant, in microporous silica with a random micro-sphere system and vitreous state sphere structure, and found results in satisfactory agreement with experimental data, and concluded that with a proper microsphere structure the numerical method could give a reasonable prediction on properties of simple nonpolar gases in silica aerogels.Ulker [16] simulated a nitrogen adsorption isotherm in silica aerogel MOF composite with conventional GCMC and it agrees with the experimental data.However, he treated the composite as a crystal, which may not work for pure silica aerogel.Le [17] constructed a silica cylindrical channel and simulated propane adsorption through the MD method.Only qualitative agreement with experimental results was obtained.Gavalda [18] applied the GCMC method to nitrogen adsorption in carbon aerogels and the results are close to the experimental isotherm in a wide range of relative pressures.In his simulation, carbon aerogels are modeled as a pack of overlapping hard spheres.Considering the structural similarities between carbon and silica aerogels, where they are all random networks consisting of nano-scale particles, it is reasonable to study nitrogen adsorption in silica aerogels with the GCMC method and hard sphere system.
In this paper, an approximate structure of silica aerogel is generated.Then a GCMC simulation of nitrogen adsorption is conducted to verify the model in the hope that it may shed light on heat and mass transport in silica aerogel.

Physical Model
Silica aerogels are usually synthesized through the sol-gel process.Silica source solution with a catalyst undergoes chemical reactions and turns into a gel consisting of a three-dimensional network of SiO 2 .After aging and drying, primary particles form and then bond into secondary particles [19], which can be seen as a basic constructional element and constitute the skeleton.Ingredients, concentration and temperature can all influence the reaction rate and in turn affect the final skeleton shape.Thus, it is hard to trace the chemical reaction process to get the microstructure of aerogels.Meanwhile, the skeleton structure is vital to the study of adsorption behavior.Several models have been put forward.Zeng [20] proposed intersecting square rod, cylindrical rod and spherical structures to study thermal conductivity.However, characteristics of the microstructure are neglected for analytical convenience, which makes it not appropriate in adsorption simulation.Park [21] suggested four hard sphere models which are classified by overlappings and connections among these spheres.These models did not take characteristics of the formation process into consideration.Schaefer [22] found this network backbone to be fractal and it could be traced back to polymerization and colloid aggregation, which can be described by the diffusion limited cluster aggregation (DLCA) method [23].Good [24] modeled the micro-structure of silica aerogel using the DLCA method and found it consistent with experimental observation.So here we adopt the DLCA method to reconstruct the micro-structure of silica aerogels.
The process can be described as follows.Initially, a group of spheres with a uniform radius are well distributed in a 100 nm ˆ100 nm ˆ100 nm box.Random walk is conducted for each sphere and then collisions are checked.If the distance between two spheres is within a certain range, they merge into one cluster and move as a whole.The periodic boundary condition is adopted during movement.Finally, a connected none-overlapping sphere system is acquired as shown in Figure 1 and this matrix of spheres is considered to be rigid during the simulation.
The final structure can be influenced by three important factors: number of particles (porosity), sphere radius (specific surface area) and randomness.To investigate the impacts of these factors, four matrixes, named M1, M2, M3 and M4, are generated as shown in Figure 2. M1 and M2 are of the same particle number and sphere radius.They are generated to study whether structural differences caused by randomness during the formation process affect adsorption capacity.The radius of M3 is close to that of M1 and M2, but the porosity is lower.It is intended to investigate the impacts of porosity.M4 is of a different particle number and sphere radius.Since the radius of the secondary particles that constitute silica aerogels usually ranges from 2 nm to 10 nm, and we want the SSA to be close to the experimental conditions, the sphere radius is set between 4 nm and 5 nm.Detailed information is shown in Table 1.
Computation 2016, 4, 18 3 of 9 be close to the experimental conditions, the sphere radius is set between 4 nm and 5 nm.Detailed information is shown in Table 1.
(a) (b)   be close to the experimental conditions, the sphere radius is set between 4 nm and 5 nm.Detailed information is shown in Table 1.

GCMC Method
Both fluid-fluid and fluid-solid interactions are given by Lennard-Jones 12-6 potential which takes the following form: Considering the composition similarities, the model developed by Hou et al. [25] for ITQ-1, a purely silica zeolite, is adopted for silica aerogel.In this model, oxygen and silicon atoms are taken into account, while for nitrogen, the force field parameter given by Gavalda et al., which has been proved effective in their experimental and numerical study of nitrogen adsorption in carbon aerogels, is employed.In this model, the nitrogen molecule is considered as a single rigid Lennard-Jones (L-J) sphere.However, a realistic nitrogen molecule consists of two atoms bonded together.So for cautious reasons, a pairwise model of N 2 [26], in which nitrogen consists of two L-J sites at a distance of 0.1085 nm, is also adopted for comparison.Parameters are listed in Table 2.The Lorentz-Berthelot mixing rules are adopted to calculate interaction parameters between different atoms.In order to reduce computational cost, a continuum model [27] which describes interactions between a single molecular L-J site of nitrogen and a solid sphere is introduced.In this model, the aerogel nanosphere is considered as a continuum of L-J sites with a density ρ.Thus, the potential between a molecular L-J site and the sphere can be written as: where R is the radius of the sphere, d is the distance between the center of the sphere and the molecular L-J site, and u is the interaction between a molecular L-J site and a tiny part of the sphere defined by relative spherical coordinates pr, θ, ϕq.By using the L-J 12-6 potential, we get: Substitute Equation (2) into Equation (3) and integrate them, and we get: when d " R, Equation (4) turns into the following form: This continuum model may not work well when d is close to R due to the continuous L-J sites assumption of the solid sphere.So when d ´R is less than 3σ N , direct summation of discrete U ij between adsorbent atoms and adsorbate molecules is adopted.Thus, the final form of potential between a sphere and a nitrogen molecule takes the form: where No is the number of oxygen atoms and N Si is the number of silicon atoms in a sphere.
In molecule-solid interactions, the minimum image convention is used.As for interactions between nitrogen molecules, simple truncation is applied.The truncation distance in this paper is set as 20 Å, approximately 5σ N .
The GCMC method is used in determining the adsorption isotherm.For such a system, it is hard to get equilibrium, especially for points near saturated vapor pressure.To increase the insertion acceptance rate, the Frank-Halsey-Hill (FHH) equation is used to estimate adsorption layer thickness δ (Å): To get B and m, we assume each adsorbed molecule occupies the same volume V. Thus, we get: where N is the number of adsorbed molecules.So we get: By fitting experimental data, we get N. Then B and m could be acquired through Equation (10).For R = 4.05 nm, B = 12.02 Å and m = 0.4984.In order to accelerate equilibrium and avoid missing some space where molecules could be adsorbed, B = 20.0Å is used in the simulation.
Prior to these operations, we first divide the simulation box into 4 nm ˆ4 nm ˆ4 nm (two times the truncation distance) grids with coordinates of (x,y,z) (x, y, z are integers from 0 to 24; totally 25 ˆ25 ˆ25 grids in a simulation box).Then we divide each grid into eight 2 nm ˆ2 nm ˆ2 nm sub-grids with coordinates of (a,b,c) (a, b, c = ´1 or 1, totally 2 ˆ2 ˆ2 sub-grids in each grid).The creation operation follows the following steps: (1) Randomly choose a solid sphere and insert a molecule within the layer thickness; (2) Check whether the molecule overlaps any solid sphere.If yes, go back to step one.Otherwise, go to step 3; (3) Identify the grid (x,y,z) and sub-grid (a,b,c) location of the molecule and get the contiguous grids of the sub-grid where the molecule is.Then calculate interactions between the inserted molecule and the existing molecules in these eight grids (x ~x + a, y ~y + b, z ~z + c) and the solid spheres.For example, if the grid coordinate is (3,6,9) and the sub-grid coordinate is (´1,1,´1), then the eight grids we need are (2 ~3, 6 ~7, 8 ~9).Through this method, we just need to calculate eight grids instead of 27.The calculation amount is significantly reduced; (4) Decide whether the creation is accepted with probability of where N is the number of existing molecules, V is the volume of the simulation box, µ is the chemical potential, Λ is the de Broglie wavelength, k B is the Boltzmann constant, U is the potential and T is the temperature.If accepted, go to the displacement operation; otherwise, go back to step 1.
As to the displacement operation, we first randomly choose one from the list of existing molecules with equal probability.Record its original coordinate pr, θ, ϕq and move it to pr `∆r, θ `∆θ, ϕ `∆ϕq.After that, the spherical coordinate is changed to a rectangular coordinate.Coordinate transformation looks time-consuming.However, as we know, the potential of positions near the surface of the solid sphere is usually lower than that of positions that are far away.Moving in a spherical coordinate ensures the system moves towards a lower potential with more probability compared to a rectangular coordinate system.Then calculate the potential as described in step 3 of creation above and accept the movement with a probability of where s is the configuration before movement and s In operation removal, we first randomly choose a molecule from existing ones with equal probability, calculate the potential as described in step 3 of the creation operation, and accept removal with a probability of A typical step includes all three operations mentioned above and they are performed once respectively in a step.The adsorption amount under pressure from 0.103 to 0.923 is calculated.Then, 100 million steps when P/P 0 < 0.6, 300 million steps when 0.6 ď P/P 0 < 0.8 and 600 million steps when 0.8 ď P/P 0 are run to get equilibrium and another 100-200 million steps are run for the statistical average.The configuration of lower pressure is used as the primary configuration of the following higher pressure.The simulation is done through self-written code.

Results and Discussion
Results of the simulative adsorption isotherm compared to the experiment [28] are shown in Figure 3.All exhibit a type IV isotherm.Good agreement with experimental data is found during middle pressure.near the surface of the solid sphere is usually lower than that of positions that are far away.Moving in a spherical coordinate ensures the system moves towards a lower potential with more probability compared to a rectangular coordinate system.Then calculate the potential as described in step 3 of creation above and accept the movement with a probability of where s is the configuration before movement and s′ is the new configuration.If it is accepted, check whether the molecule is moved to another grid.If yes, add a new molecule at the end of the list of the new grid and delete it from the original grid; if not, record the new position.If it is not accepted, move it back.Then go to operation deletion.
In operation removal, we first randomly choose a molecule from existing ones with equal probability, calculate the potential as described in step 3 of the creation operation, and accept removal with a probability of A typical step includes all three operations mentioned above and they are performed once respectively in a step.The adsorption amount under pressure from 0.103 to 0.923 is calculated.Then, 100 million steps when P/P0 < 0.6, 300 million steps when 0.6 ≤ P/P0 < 0.8 and 600 million steps when 0.8 ≤ P/P0 are run to get equilibrium and another 100-200 million steps are run for the statistical average.The configuration of lower pressure is used as the primary configuration of the following higher pressure.The simulation is done through self-written code.

Results and Discussion
Results of the simulative adsorption isotherm compared to the experiment [28] are shown in Figure 3.All exhibit a type IV isotherm.Good agreement with experimental data is found during middle pressure.First, there is a little gap between the adsorption isotherms of the two nitrogen models.In the dimer model, electrostatic potential is not taken into consideration due to limited computational resources, which is responsible for the overall lower adsorption amount.The systematic underestimation of gas sorption is not caused by gas models.First, there is a little gap between the adsorption isotherms of the two nitrogen models.In the dimer model, electrostatic potential is not taken into consideration due to limited computational resources, which is responsible for the overall lower adsorption amount.The systematic underestimation of gas sorption is not caused by gas models.
There is little discrepancy between M1 and M2, which indicates that structural differences caused by randomness during the formation process do not notably affect the adsorption capacity.Simulation results of M1 and M2 are the closest to the experimental thanks to parameter likeness, while M3 and M4 have a lower adsorption amount compared to the experimental isotherm due to less specific surface area or less porosity.Several reasons might be responsible for differences between the isotherms of M1 and the experiment: (1) The first reason is microstructural differences of the sphere.There are various groups, such as methyl and oxhydryl, on the surface of realistic particles.However, only oxygen and silicon atoms are considered in the model.Additionally, surface atom density plays an important role in determining adsorptive capacity, but in the model it may not agree with the realistic surface density.Lastly, realistic particles may overlap nearby particles.However, they are simply modeled as non-overlapping spheres in the simulation.(2) The second reason is pore size distribution.Though the porosity is close, the pore size distribution may not be.(3) The third reason is neglected adsorption sites.This simulation only pays attention to adsorption on the surface of spheres and in the pores among them.Possible adsorption which takes place in the microchannels within the sphere is neglected.
Under low pressure, specific surface area affects the adsorption amount.A structure with more SSA takes in more nitrogen, but discrepancies are not significant.When it comes to high pressure (P/P 0 ě 0.72), porosity plays a dominant role in determining the adsorption capacity.The nitrogen uptake of M1 and M4 is similar due to similar porosity, though the SSA is different.In contrast, M1 and M3 are of similar SSA, but there are significant differences in their nitrogen load.
Snapshots of adsorption at different relative pressures in M1 are shown in Figure 4.At first, molecules are mainly absorbed on the surface of the spheres.As the relative pressure increases, multiple layers of adsorbed molecules form.Under high pressure, nitrogen molecules gradually fill the pores among the spheres.The number density distribution (Figure 5) shows that the adsorbed molecule exhibits an intermediate state between gas and liquid.
less specific surface area or less porosity.Several reasons might be responsible for differences between the isotherms of M1 and the experiment: (1) The first reason is microstructural differences of the sphere.There are various groups, such as methyl and oxhydryl, on the surface of realistic particles.However, only oxygen and silicon atoms are considered in the model.Additionally, surface atom density plays an important role in determining adsorptive capacity, but in the model it may not agree with the realistic surface density.Lastly, realistic particles may overlap nearby particles.However, they are simply modeled as non-overlapping spheres in the simulation.( 2) The second reason is pore size distribution.Though the porosity is close, the pore size distribution may not be.(3) The third reason is neglected adsorption sites.This simulation only pays attention to adsorption on the surface of spheres and in the pores among them.Possible adsorption which takes place in the microchannels within the sphere is neglected.
Under low pressure, specific surface area affects the adsorption amount.A structure with more SSA takes in more nitrogen, but discrepancies are not significant.When it comes to high pressure (P/P0 ≥ 0.72), porosity plays a dominant role in determining the adsorption capacity.The nitrogen uptake of M1 and M4 is similar due to similar porosity, though the SSA is different.In contrast, M1 and M3 are of similar SSA, but there are significant differences in their nitrogen load.
Snapshots of adsorption at different relative pressures in M1 are shown in Figure 4.At first, molecules are mainly absorbed on the surface of the spheres.As the relative pressure increases, multiple layers of adsorbed molecules form.Under high pressure, nitrogen molecules gradually fill the pores among the spheres.The number density distribution (Figure 5) shows that the adsorbed molecule exhibits an intermediate state between gas and liquid.

Conclusions
Through the comparison of simulative results and the experimental data, validity of the combination of the sphere system generated through the DLCA method and the continuous-discrete hybrid potential model is proved.The simulative isotherm corresponds to type IV.Adsorption takes place on the surface of spheres under low pressure.As pressure increases, nitrogen molecules gradually fill the pores among the spheres.Specific surface area affects the adsorption capacity under low pressure, while porosity controls the adsorptive capacity under high pressure.
The model did not take possible adsorption sites within the sphere into consideration, which might be the main reason for the discrepancies between the simulative and experimental data.The interior structure of the sphere needs to be improved for a better prediction.Surface atoms could be altered to study the impacts of surface modification on adsorbability.The thickness of the adsorbed layer grows rapidly with relative pressure and could reach a few nanometers, to the scale of pores, which could significantly affect the mass and energy transport in silica aerogel.Exact determination of the amount of the adsorbed gas molecules could help the multi-scale simulation of the transport process in the future.

Conclusions
Through the comparison of simulative results and the experimental data, validity of the combination of the sphere system generated through the DLCA method and the continuous-discrete hybrid potential model is proved.The simulative isotherm corresponds to type IV.Adsorption takes place on the surface of spheres under low pressure.As pressure increases, nitrogen molecules gradually fill the pores among the spheres.Specific surface area affects the adsorption capacity under low pressure, while porosity controls the adsorptive capacity under high pressure.
The model did not take possible adsorption sites within the sphere into consideration, which might be the main reason for the discrepancies between the simulative and experimental data.The interior structure of the sphere needs to be improved for a better prediction.Surface atoms could be altered to study the impacts of surface modification on adsorbability.The thickness of the adsorbed layer grows rapidly with relative pressure and could reach a few nanometers, to the scale of pores, which could significantly affect the mass and energy transport in silica aerogel.Exact determination of the amount of the adsorbed gas molecules could help the multi-scale simulation of the transport process in the future.
1is the new configuration.If it is accepted, check whether the molecule is moved to another grid.If yes, add a new molecule at the end of the list of the new grid and delete it from the original grid; if not, record the new position.If it is not accepted, move it back.Then go to operation deletion.