THE EFFECT OF ADSORBED MOISTURE ON INTERPARTICLE FORCES

A Molecular Dynamics (MD) simulation technique is used to investigate the effect of moisture on interparticle forces. This is carried out by simulating two opposing moist quartz surfaces at different separations. The spontaneous formation of a liquid bridge, was observed with only three and a half monolayers on the surfaces, corresponding to a relative humidity of about 70% according to the BET-isotherm. In the article a brief account of the methodology is given. Results are shown in the form of snapshots of the position of the molecules at a given instant in time and mean forces between the quartz plates as a function of surface separation. The relation of the results to the macroscopic theory of liquid bridging is also mentioned, as are the density and diffusivity profiles of the adsorbed water layers.
KEY WORDS: Liquid bridging; molecular dynamics; inter particle forces; moisture; adsorption.

INTRODUCTION
In the processing industry design and optimization of unit operations involving the handling and/or processing of powders are almost universally based on bulk measures of powder flowability, which are determined, for instance, in the 'Jenike shear cell'. In the past few years numerical methods have emerged such as 'hard sphere simulation', in which the dynamics of individual particles are simulated, making assumptions about the interparticle and particle/wall interaction forces, (see, for instance, 1).

In order to progress in the modelling of powder behaviour, relating the bulk flowability indicators to the particle properties and the ambient conditions, or improving the predictive power of the 'hard sphere simulation' techniques, fundamental knowledge of the particle/particle and particle/wall interaction is required.

If the particle surface is completely pure and the surface morphology and elastic properties of the material are known, particle/particle interaction energies can be calculated approximately. In most real situations, though, interparticle forces are strongly affected by materials (gases, moisture) adsorbed on the particle surface. At high humidities 2 the effect of moisture on particle interaction can be described rather well by the macroscopic theory of liquid bridging 3-5, although some questions remain to be answered. A good understanding of the effect of adsorbed moisture at low and moderate relative humidities is still lacking. Although some pioneering work has been done (e.g.6) much is still based on conjecture. It remains an open question whether liquid bridging occurs at lower humidities, and how many monolayers must be adsorbed before the liquid is sufficiently mobile for the formation of bridges. This article presents a new approach for investigating the effect of adsorbed moisture (or other adsorbed agents) on interparticle interaction, which offers improved possibilities for fundamental understanding and for quantifying interparticle interaction in different systems.

The method is the technique of Molecular Dynamics (MD) simulation, which provides a method of investigating these phenomena on nanometer scale. The term Molecular Dynamics simulations is being used prolificly in the field of the modelling of granular dynamics by hard sphere simulation. The work in this article is concerned with Molecular Dynamics in the true sense: the simulation of the dynamics of individual molecules, which involves the evaluation of force fields very different from those in granular dynamics.

METHODS
MOLECULAR DYNAMICS
MD simulation 7-9 is a well established computational technique in different areas of science, for example in protein engineering, drug design, spectroscopy and the study of membranes and (mixtures of) liquids. In particle technology it is less used. A short introduction is, therefore, in order.

Performing MD simulations means numerically solving Newton's equations of motion for a large number of molecules, which may have internal degrees of freedom as well. The coupled differential equations are solved using a leap-frog algorithm 10, typically using a time step of 2 fs. The forces acting on the atoms at each time step are determined by the interaction potential, the 'force field', which consists of different (empirical) terms, such as harmonic or rigid bond and angle potentials, Coulomb interactions and Lennard-Jones potentials, mimicking Born repulsion and London dispersion. The values of the various parameters are in this case taken from the GROMOS/GROMACS force field 9. To save computational time, interactions are only calculated with neighbouring atoms within a certain cut-off radius Rc, since most forces level off very quickly with interparticle distance. Furthermore, only a finite number of molecules is simulated in a predefined box and, to avoid border effects, periodic boundary conditions are applied. This means that the simulation box is surrounded by its periodic images, and that the interactions of an atom are only calculated with its closest (periodic) neighbours. Simulations were carried out at constant temperature, using the weak coupling thermostat algorithm 11. The final file, the 'trajectory', containing positions, velocities and forces of all atoms at each time step may be used for further analysis.

THE STARTING STRUCTURE
Glass is a good first type of solid material to look at, since glass spheres are often used as a model powder in powder technology. The disadvantage of glass, however, is its amorphous polymer structure. Generation of a starting structure of amorphous glass for MD simulations would be a project in itself 12,13, with probably only little chance of improvement over a system built from a more simple compound with a regular crystal lattice: -quartz.



Figure 1 Two interacting surfaces with adsorbed layers of moisture.

The starting structure for the simulations, containing positions and velocities of all atoms at time zero, consists of two opposing quartz surfaces of finite thickness, separated on one side by water layers and on the other side by vacuo, see Figure 1. The quartz crystals were generated by multiplication of a single quartz unit cell along the (1 0 -1) Miller plane, which is a well studied surface 14. The surface oxygen atoms were hydrogen terminated as shown in Figure 2. The charges in Figure 2 were taken from the literature 15. During the equilibration and the simulations the bulk crystal atoms were kept fixed in space.



Figure 2. Crystal surface detail illustrating the applied topology. The crystal bulk atoms are fixed in space without harmonic bonds to neighbours. The surface atoms are connected with bonds as illustrated. In the left of the two tetraedals the different atom labels as used in the GROMOS/GROMACS force field are indicated, in the right one the atomic charges different from zero are indicated. The bold oxygen is used to define the distance between the two crystal surfaces.


After the surface had been generated, a layer of SPC water molecules 16 was added and the system mirrored in order to obtain two opposing surfaces (see Figure 3), the total system contained up to 10600 atoms. The structure was then allowed to equilibrate for 100 ps to remove bad interactions between the crystal surfaces and the water layers and between the water molecules themselves.



Figure 3 A start configuration. The single crystal dimensions are 5.49 * 4.91 * 1.45 nm.

THE CRYSTAL BULK CORRECTION (CBC)
Because of the strong anisotropy of the system, one has to correct for the neglected bulk part of the crystal-water interactions outside the cut-off radius by applying an extra potential to the water molecules: the Crystal Bulk Correction (CBC). This potential is obtained by integrating the Lennard-Jones pair potential:

(1)

over an infinite surface 17, where r is the interparticle distance and C6(IJ) and C12(IJ) are attractive and repulsive interaction parameters for atoms of type I and J. The final result for an atom of type I at a distance d from the surface is:

(2)

the parameters C6 and C12 are defined as:



Here rN is the number density of the crystal, xJ the mole fraction of the atom type J and the summation is taken over all atom types in the bulk.

COMPUTERS AND SOFTWARE
The simulations presented in this article are carried out using the GROMACS software package 8 on a Silicon Graphics Power Challenge and a Silicon Graphics O2 workstation.

RESULTS
Several simulations of 400 ps have been performed at 300 K, using a time step of 2 fs and a cut-off radius, Rc of 1.2 nm, with different numbers of water molecules (up to 1500) and different surface separations. Various properties of the adsorbed layers (e.g. density and diffusion profiles) have been calculated. The effect of the CBC on these properties was investigated as well 18. Detailed studies of the density profiles (see also the snapshots below) clearly show a layer packing of the water molecules close to the crystal surface up to three distinct layers. The diffusion profiles indicated that the first layer is almost immobile and that also the second layer is significantly less mobile than bulk water.

For the systems containing 1500 water molecules with about three and a half adsorbed monolayers, corresponding to relative humidities in the region of 70%, the force acting on the crystals was calculated for different separations of the surfaces.


Table 1 The time averaged forces per unit surface area with their statistical errors. The first column gives the separation between the crystal surfaces, D, the second column the average force on the single crystals. All forces are attractive, and have been averaged over the two opposing crystals as well.
D (nm) F (bar)
6.0 32±178
3.7 27±167
3.0 466±160



The evaluation of the total force is a simple summation of all forces acting on the individual crystal atoms. On a short time scale the force was strongly oscillating, and in order to obtain a statistically significant result averaging over a 'long' time was necessary. The results are shown in Table 1. Note the significant increase in the (attractive) force at 3.0 nm separation.

This large increase of the attractive force at 3.0 nm can be understood by inspecting the final structures of this simulation. In contrast to the final stable structures at 6.0 and 3.7 nm, which exhibit two distinct adsorbed water layers, at 3.0 nm the structure shown in Figure 4 is obtained, a liquid bridge-like structure.This structure already appears within the first 20 ps and remains stable during the rest of the simulation time. Note that the cavity has a cylindrical shape. Under the constraints of the periodic boundary conditions, this geometry leads to the minimal surface area for the case shown.



Figure 4. Front view of end configuration of the simulation with 3.0 nm separation distance between the crystal surfaces.

DISCUSSION, CONCLUSIONS AND OUTLOOK
These MD simulations show liquid bridging to be possible for conditions where the adsorbed layers are relatively thin. A phenomenological estimation of the force between two particles joined by a bridge, based on macroscopic theory 3, gives a result of the same order of magnitude to the force in Table 1 at 3.0 nm separation. This confirms the experimental result of Fisher and Israelachvili 19 that macroscopic theory is applicable to bridges with a curvature of molecular dimensions, although they found some deviation for water. Our results are, at this stage, not of a statistical significance high enough to permit a precise quantitative comparison like theirs.

At a given surface separation, the simulated system is a closed one under conditions of constant volume. When the available volume between the crystals exceeds the volume of the dense phase, the system contains two fluid phases in equilibrium. No gaseous H2O molecules are observed most of the time since the total volume available to the gas is only a fraction of the volume taken up by a single gaseous molecule. The 'gas phase' therefore consists of a molecule escaping the adsorbed layers once in a while, only to be caught up quickly again.

According to the BET-isotherm, three and a half monolayers correspond to a relative vapour pressure (or relative humidity) of around 70%. If the system is always in equilibrium one would expect bridge formation to take place at the surface separation when bridging becomes energetically advantageous, and the resulting bridge curvature should correspond to about the same relative humidity as the adsorbed layer thickness. On the other hand, bridge formation may be delayed by non-equilibrium effects. The adsorbed water layers will not join in a bridge until they 'feel' eachother, and there will then be a hysteresis effect in bridge formation and breakage. This is so in a physical system, and particularly so in the present system, where a cut-off radius of 1.2 nm is used. Once a bridge is formed, equilibrium (equal chemical potential) between the molecules in the bridge and those of the remaining adsorbed layers (clearly visible in Figure 4, the thickness of these has not yet been determined, but appears to be a couple of monolayers) should exist. The bridge curvature in Figure 4 is about 1.1 nm, and this would, according to the Kelvin Equation using the macroscopic value for the surface tension of water, correspond to a relative humidity of 62%. Taking into account the uncertainties involved, this shows that the observed phenomena are roughly consistent with what one would expect. It should be noted in this context that bridge formation in the system with other thicknesses of the adsorbed layers has not yet been studied, neither has the effect of changing the cut-off radius.

One issue, which it may be possible to investigate using this method in the future, is the question to what extent adsorbed moisture accumulates in valleys on an undulated particle surface, due to surface tension effects. It is likely that, if the surface is reasonably hydrophilic, some adsorbed moisture is present in contact points whatever the morphology of the particle surface.

A contact point will, if there is adsorbed moisture on the surface, look something like the sketch in Figure 5. A phenomenon akin to the bridging seen in Figure 4 will be present at the edges of the contact point, and at equilibrium the attractive force at the periphery will be balanced by a repulsion in the central region. There, in a region of high pressure, some water molecules will be 'squeezed out' (higher chemical potential due to the higher pressure), but some may remain 19 (their chemical potential being lowered due to proximity to the particle surfaces). In this region, deformation of the solid will also take place.





Figure 5 Sketch of an interparticle contact point.


Simulations of the sort of phenomena described above may, together with some of the experimental studies of intersurface interation being published (for instance 19,20) provide a key to understanding and estimating interparticle forces whether of a smooth or more irregular surface topology.

It is the longer term aim of this research to arrive at a description of a contact point, applicable to systems involving various adsorbed agents and various surfaces and thus gain quantitative information about particle/particle and particle/wall interaction energies.

Work will also be aimed at a quantitative thermodynamic description of the bridging seen, and the effect of the different parameters (model parameters as well as simulation parameters) on the formation of the bridge. The behaviour of adsorbed liquid as a function of surface properties, for example the charges on the surface, is another promising line of investigation.

A next step may be a series of simulations with crystal separations between 3.7 nm and 3.0 nm to check whether there is a significant increase of the calculated forces before the formation of a bridge. The results of these simulations may be useful to check the proposition of enhancement of interaction force by adsorbed layers without bridging. It is also possible to perform simulations with crystals moving relative to each other. This offers the possibility to investigate dynamic phenomena in particulates.

REFERENCES
1. Hoomans, B.P.B., Kuipers, J.A.M., Briels, W.J., van Swaaij, W.P.M., 1996, Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: a hard-sphere approach, Chem. Eng. Sci., 51: 99-118.

2. Turner, G.A. and Balasubramianian, M., 1974, Investigations of the contributions to the tensile strength of weak particular masses, Powder Technol., 10: 121-127.

3. Adams, M.J., and Perchard, 1985, The cohesive forces between particles with interstitial liquid, IChemE Symp. Ser., 92: 147-160.

4. Lian, G., Thornton, C., Adams, M.J., 1993, A theoretical study of the liquid bridge forces between two rigid spherical bodies, J. Coll. Int. Sc., 161: 138-147.

5. Gao, C., 1997, Theory of menisci and its applications, Appl. Phys. Lett., 71: 1801-1803.

6. Coelho, M.C., Harnby, N., 1978, Moisture Bonding in Powders, Powder Technol., 20: 201-205.

7. Allan, M.P., Tildesly, D.J, 1987, Computer Simulation of Liquids, Oxford University Press, Oxford.

8. van der Spoel, D., Berendsen, H.J.C., van Buuren, R., Apol, E., Meulenhoff, P.J., Sijbers, A.L.T.M. and van Drunen, R., 1996, GROMACS Manual user manual, BIOSON Research Institute 1996, internet: http://rugmd0.chem.rug.nl/ gmx1996.

9. van Gunsteren, W.F. and Berendsen, H.J.C., 1990, Computer simulation of molecular dynamics: Methodology, applications and perspectives in chemistry, Angew. Chem. Int. Ed. English 29: 992-1023.

10. Berendsen, H.J.C. and van Gunsteren, W.F., 1986, In: Molecular-Dynamics simulation of statistical-mechanical systems, Proceedings of the International School of Physics "Enrico Fermi", course 97. G. Ciccotti and W.G. Hoover (eds.), North-Holland, Amsterdam, p. 43-65.

11. Berendsen, H.J.C., Postma, J.P.M, van Gunsteren, W.F., Di Nola, A. and Haak, J.R., 1984, J. Chem. Phys. 81: 3684-3690.

12. Webb, E.B., Garofalini, S. H., 1994, Molecular Dynamics simulation of the approach and withdrawal of a model crystalline metal to a silica glass surface, J. Chem. Phys, 101 (11): 792-801.

13.Della Valle, R.G., Andersen, H.C., 1992, Molecular Dynamics simulation of silica liquid and glass, J. Chem. Phys. 97 (4): 2682-2689.

14. Jánossy, I. and Menyhárd M., 1971, LEED Study of Quartz Crystals, Surface Science, 25, 647-649.

15. Lee, S.H. and Rossky, P.J., 1994, A comparison of the structure and dynamics of liquid water at hydrophobic and hydophilic surfaces - Molecular Dynamics simulation study, J. Chem. Phys. 100 (4): 3334-3345

16. Berendsen, H.J.C., Postma, J.P.M., van Gunsteren, W.F. and Hermans, J., Interaction models for water in relation to protein hydration. In: Intermolecular forces, B. Pullman (ed.), p. 331-342. Reidel, Dordrecht.

17. Israelachvili, J.N., 1992, Intermolecular and surface forces, Academic Press, London.

18. Wensink E.J.W., Apol, M.E.F. and Hoffmann, A.C., 1997, Properties of adsorbed water layers and the effect of adsorbed layers on interparticle forces, submitted

19. Fisher, L.R. and Israelachvili, J.N., 1981, Direct measurement of the effect of meniscus forces on adhesion: a study of the applicability of macroscopic thermodynamics to microscopic liquid interfaces, Coll. Surf. 3: 303-319

20. Jarvis, S.P., Yamada, H., Yamamoto, S.I., Yokumoto, H. and Pethica, J.B., 1996, Direct mechanical measurement of interatomic potentials Nature 384: 247-249.

LIST OF SYMBOLS
C6 attractive parameter of CBC potential for type I
C6(IJ) attractive parameter of Lennard-Jones potential between types I and J
C12 repulsive parameter of CBC potential for type I
C12(IJ) repulsive parameter of Lennard-Jones potential between types I and J
d atom/surface distance
D surface/surface distance
F average force per unit area perpendicular to the surface
r distance between molecular or atomic particles
V Lennard-Jones pair potential
w crystal bulk correction potential
xI mole fraction of atom type I
rN number density