V

Colloidal gels are intermediates in the production of highly porous particle systems. In the production process, the gels are fragmented after their creation. These gel fragments consolidate to particles whose application-technological properties are determined by their size and porosity. A model of the consolidation process is proposed: The consolidation process of a gel fragment is simulated with the Molecular Dynamics (MD) method with the assumption of van der Waals forces in interplay with the thermal motion as driving forces for the consolidation. The simulation results are compared with experimental data and with a Monte Carlo (MC) simulation.


σ
(finite) distance at which inter-particle potential is zero r distance between the particles V interaction energy A cc

Hamaker constant for colloid-colloid interactions A cs
Hamaker constant for colloid-solvent interactions a radius of the colloidal particles fractal dimension n number of primary particles in the aggregate r agg radius of the aggregate r p radius of the primary particles

Introduction
Precipitated silica is used as a filling material in the manufacture of shoes, tires and general rubber goods.It is used in tires in order to reduce the rolling resistance, while the abrasion resistance and wet grip are increased [1].An alternative to precipitated silica is pyrogenic silica.It is produced by flame hydrolysis in the gaseous phase.Pyrogenic silica is a special form of amorphous silica with a surface area of 30-440 m²/g [2] and a primary particle size of 5-30 nm [3].Up to now, it has not been possible to form qualities of pyrogenic silica by a liquid route of manufacturing.However, because flame hydrolysis is expensive, it is desirable to find a way to produce highly porous silica by a liquid route, e.g., by a precipitation reaction.
In the industrial inorganic precipitation of silica, sodium silicate solution and sulfuric acid are fed into a water tank while the solution is stirred.The solution gelates after about 30 min.This gel is fragmented by the stirrer.The gel fragments consolidate to aggregates which are characterized by their size and porosity.
The simulation of aggregates makes the investigation of the development of the aggregate structure as a function of time possible.These insights into the aggregate structure can be used to tailor the aggregate structure by designing new processes in the future [4].
So far, the consolidation of aggregates has been modeled i.a. by the Stokesian Dynamics method [5] and by the Monte-Carlo method [6], in order to describe the development of the aggregate structure as a function of time.
The Stokesian dynamics method, as described in [5], allows the simulation of the motion of particles which are located in a fluid.For this purpose, the Langevin form of the 2nd law of Newton (F = m × a), which describes the particle motion in a Newtonian fluid, is solved.The Langevin form states that the sum of forces to which a particle is exposed in a fluid equals the mass of the particle times its acceleration [5].
Here, F H are the hydrodynamic forces to which the particle is exposed because of its movement relative to the fluid, F P stands for the non-hydrodynamic forces (interparticular and external forces) and F B stands for stochastic forces which cause Brownian motion of the particle [5].
First attempts of the modeling of the structure of colloidal aggregates can be found in [7,8].In [7], 2D-simulations of colloids with cohesive forces are described.The model described (Sticky Sphere Model) allows for some simplification of the Stokesian dynamics method.The hydrodynamic forces are strongly simplified.The particles are only affected by Stokes' frictional force (free draining approximation) [9]: F R is the frictional force in the direction of the particle velocity v, r is the particle radius and η is the dynamic viscosity of the fluid.Thus, all influences caused by flow fields which are generated by the particles, and which can influence the flow of neighboring particles, are neglected.The non-hydrodynamic forces are simplified in a manner that the cohesive particles can roll on each other.The bond between two particles breaks if a critical value of the normal force is exceeded.The Brownian motion is neglected [4,7].A model is described in [10] that looks at the hydrodynamic and interparticular forces in more detail.Upon the assumption that Stokesian friction takes effect on the particles, the frictional force is overestimated.This is especially the case for those particles which are completely surrounded by other particles.It is assumed that the hydrodynamic frictional force only takes effect on the particle surface which is directly exposed to the liquid [10].The interparticular forces are obtained by DLVO theory [11][12][13].The model according to [14,15] is used for particles which are in contact with each other.It has found wide use in the modeling of granular solids with the help of the Discrete Element Method (DEM) [16].By these assumptions, among other things, the deformation of three-dimensional aggregates is simulated, e.g., in [4,10,17,18].Furthermore, the consolidation of gel fragments is simulated by statistical methods, such as the Monte Carlo method.The model described in Schlomach (2007) [19] shows good agreement to the time-dependent behavior of a structural factor (fractal dimension) of the computed particle structures with experimental data.Nevertheless, this model is not based on the appropriate interaction potentials between the particles, e.g., according to van der Waals.These interactions are in interplay with the thermal motion of the particles.They are the driving forces for the consolidation [20,21].
MD simulation is used here in consideration of the argument above.Further arguments for this decision are given in Section 2.1.2.
The paper is organized as follows: In Section 2.1, a model of the consolidation of a gel fragment is introduced taking van der Waals interactions and thermal motion into account.In Section 2.2, simulation results of the consolidation behavior with the help of MD simulation are presented.The consolidation of a gel fragment is computed and compared with experimental data and a former MC simulation of Schlomach and Kind (2007).The temporal development of the structural factor that is measured in Schlomach (2004) [22] is only valid for a certain experiment.Since the experimental investigation of the consolidation of microscale gel fragments is difficult, an experiment for quantifying the consolidation of macroscale silica gels was developed by Sahabi and Kind (2011) [23].They suggest the oedometer test for the investigation of the consolidation of silica gels and measured consolidation curves for silica gels under various conditions (pH, T).Furthermore, in Section 2.2, results of the consolidation of a gel volume element by MD simulation are given and compared to experimental data.Finally, the findings are summarized in Section 2.3.

Modeling and Simulation
It is assumed that the driving forces for the consolidation are van der Waals forces [20].The shrinkage is finally stopped by the remaining repulsive forces [24].If there is a certain gel structure and van der Waals forces are applied to this structure, the structure will compact until further attraction is no longer possible because of the repulsive forces.The system is in equilibrium.Here, the sum of attractive (especially van der Waals) forces and repulsive forces is assumed to follow the Lennard-Jones (12-6) potential.In Figure 1, the Lennard-Jones (12-6) potential (V LJ ) [25] is shown qualitatively and the equilibrium condition (energy minimum) is obvious.In Figure 2, the motion of one particle is shown schematically.The Brownian motion, which is dependent on the temperature, causes the particle shown in Figure 2 to leave one minimum and reach another.

Statistical Physics
The properties (in physics: position and momentum) of the systems of many particles (e.g., atoms, molecules, particles) are analyzed in statistical physics.Statistical ensembles are thus an important tool in statistical physics.An ensemble is an idealization consisting of a large number of copies of a system.These are considered all at once; each of them represents a possible state of the real system [26,27].It is usual to investigate systems at constant pressure and constant temperature in an experimental situation.An isothermal isobaric ensemble (NPT ensemble) is an amount of systems of a given number of particles N, a given pressure p and a given temperature T [28].
The Ergodic hypothesis is the fundament of statistical physics.In the Ergodic hypothesis, it is assumed that the temporal average of a measured variable (e.g., of the pressure) is equal to the ensemble average: 2.1.2.MC vs. MD Monte-Carlo (MC) simulations are especially appropriate for the calculation of statistical average values of a certain quantity.The Monte-Carlo method is-in the form it was developed by Metropolis et al. [6] and extended in several ways [29,30]-a procedure for estimating average values for ensembles of constant temperature, e.g., the isothermal isobaric ensemble.The quantities, pressure and temperature are defined in advance [31].In molecular dynamics, Newton's equations of motion of N particles are numerically computed.An advantage of the MD method is that it delivers information about the temporal dependency and the extent of deviations of position and momentum, while MC is merely dealing with variables of position and does not provide any information about the temporal dependency of the deviations.It has been necessary to use the Monte-Carlo method in order to be able to fix the temperature and the pressure in a simulation, and abandons the possibility of gaining dynamical information in the same simulation [31].Andersen [31] describes an MD method for performing simulations of a fluid at constant pressure without the loss of a dynamical description of the fluid.This method is a hybrid of MD and MC [32].In Andersen's method, only changes in the volume of the MD simulation cell were possible, but not in its shape [33].Parrinello and Rahman [33,34] have extended Andersen's method so that changes both of the volume and the shape of the MD cell become possible.Moreover, they have shown the usefulness of the application of their new method for changes in the structure of solids [32].Nosé [32] proposes an MD method which can generate configurations belonging to the NPT ensemble by scaling time (with seconds) [35].Hoover [35] developed a slightly different set of equations, free of time scaling.They both are referred to in the following.

MD Model
The consolidation of gel fragments is simulated by the MD tool, LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [36].The time integration in LAMMPS is performed on Nosé-Hoover-style equations of motion [32,35], which are designed to generate positions and velocities sampled from the isothermal isobaric (NPT) ensemble.This is achieved by adding dynamic variables which are coupled to the particle velocities (thermostating) and simulation domain dimensions (barostating).The external pressure of the barostat can be specified as either a scalar pressure (isobaric ensemble) or as components of a symmetric stress tensor (constant stress ensemble) [37].The constant stress ensemble is used in the case of the simulation of the shrinkage of a gel fragment, and the isobaric ensemble is used for the simulation of the shrinkage of the particle network.The computations in the simulations here are carried out with dimensionless quantities.For the sake of being specific and clear about the equations used, they are repeated in the Appendix.

Structure of the Simulated Gel Fragment
A fractal aggregate according to Schlomach and Kind (2007) [19] consisting of 1,000 primary particles (see Figure 3) is used for the simulation of the consolidation of a gel fragment.An algorithm based on a model developed by Sutherland (1967) [38] is used to generate this cluster consisting of primary particles.The collision probability for the cluster-cluster aggregation is obtained from the diffusion limited aggregation constant [39].

Interaction Energies
The solvent environment in the simulation is achieved by a high density lattice of solvent molecules.The Lennard-Jones (12-6) potential V [25] is used for the interaction energy between the solvent molecules of size σ: V R describes Pauli/Born repulsion at short ranges due to overlapping electron orbitals, and V A describes attraction at long ranges (mainly van der Waals attractive forces).
A colloidal particle has a size of 2a > σ.Evraers and Ejtehadi [41] describe each colloidal particle as an integrated collection of Lennard-Jones particles of size σ.Their formula for the colloid-colloid interaction energy between two particles gives with a 1 = a 2 = a (monodisperse): where A cc is the Hamaker constant and a is the radius of the colloidal particles.
The colloid-solvent interaction energy is given by where A cs is the Hamaker constant and a is the radius of the colloidal particles.This formula is derived from the colloid-colloid interaction energy formula (see [41]), letting one of the particle sizes be zero.Hence, pairwise interactions between large colloidal particles and small solvent particles are computed using these three formulas [37] in case of Lennard-Jones (LJ) interactions between particles.The additional electrostatic repulsive force (V EL ) is added to the LJ potential for DLVO interactions (V DLVO ): Here, ε is the relative permittivity of the fluid, ε 0 is the electric constant, κ is the inverse screening length, and ψ is the surface potential [37].

Consolidation of a Single Cluster
A cluster consisting of 1,000 particles in the environment of ~400,000 solvent molecules is densified in LAMMPS as an isothermal isobaric (NPT) ensemble under the assumption of the interaction energies between the particles according to Equations (4-6) (LJ interactions).The isothermal isobaric ensemble was chosen in order to imitate the real situation where a gel cluster is confronted by a certain pressure in three dimensions in a stirred vessel.Thus, a certain temperature and a certain pressure can be allotted to the particle system.Figure 5 shows the computed model of a gel fragment before and after shrinkage.In the following, these simulation results are compared with the experiment of Schlomach and Kind [22].Figure 6 shows data from a stirred vessel experiment [22], along with the MC simulation of Schlomach (2006) [40] and the simulated data (LJ) computed here.The fractal dimension is plotted as a function of the time after the gelation point.The formula for the computation of the fractal dimension d f of the gel aggregates in the simulations is where n is the number of primary particles in the aggregate, and r agg and r p are the radii of the aggregate and the primary particles, respectively.The agreement of the MD simulation with the experimental data is not as good as with the MC simulation.The theoretical MD curve matches the experimental data at the beginning and at the end of the simulation; there are deviations in between.
The agreement shown of simulation data with the experimental data is admittedly only achieved by weighting the time step with a factor of 100 (1 min = 100 time steps).This indicates that the assumptions for the simulation are still incomplete.This concerns both the Lennard-Jones interaction energies between the particles as an approach for van der Waals and Pauli/Born interactions, and the use of an isothermal isobaric (NPT) ensemble.Apparently, the Nosé-Hoover thermostat as a deterministic method for the achievement of a constant temperature in the simulation (isothermal ensemble) is appropriate for the problem.Additionally, the equations used for the generation of an isobaric ensemble, and the assumption of van der Waals forces in interplay with the thermal motion as driving forces for the consolidation seem appropriate for the problem.Furthermore, the electrostatic repulsive force according to Equation ( 8) was taken into account to achieve DLVO interactions between the particles.An additional MD simulation leads to the data shown in Figure 6.The results show that DLVO interactions are not appropriate for this problem.(rhombi) along with the fractal dimension computed in the MC [40] resp.MD simulations (curves).

Consolidation of a Particulate Network
In order to simulate the shrinkage of an expanded gel sample, the highly porous cluster (left side of Figure 4) is stringed together regularly in the x-and z-direction (25 × 40 clusters of 1,000 particles).Thus, a particulate network consisting of 10 6 particles is generated.The kind of arrangement is shown in Figure 7 by means of 2 × 2 clusters.The simulation leads to the results shown in Figure 8 with the assumption of van der Waals attractive and Pauli/Born repulsive interaction energies between the particles, and consolidation in the z-direction at a constant temperature.A scalar pressure in the z direction is chosen in order to imitate the experimental situation of a gel sample in the oedometer test.The time step was weighted with factor 10 (1 min = 10 time steps) in order to bring the MD simulation results into reasonable agreement with experimental data.
The sharp bend of the consolidation ratio computed is because of the low resolution of the time step in the initial phase of this square-root-of-time plot.

Conclusions
The consolidation of gel fragments and of a particle network consisting of gel fragments have been investigated on the assumption of the interplay of van der Waals energies and thermal motion as driving forces for the shrinkage.Two molecular dynamics simulations have been performed with the tool LAMMPS: The fractal dimension of the MD simulation of the shrinkage of the gel fragment was compared to the results of a former MC simulation and also to experimental data.After time scale adjustment (factor 100) the simulation data could be aligned in agreement with the experimental data.Thus, the assumptions of the simulation have apparently been correct.Van der Waals forces in interplay with the thermal motion are considered to be driving forces for consolidation.
Furthermore, an MD simulation of the consolidation of a gel volume element was performed.The consolidation ratio of the consolidating particle network was compared to experimental data from oedometer tests.Good agreement of experiment and simulated data was also found here after time scale adjustment (factor 10).However, in order to describe the temporal course of the consolidation in more thoroughly, the fluid flow through the particle network needs to be taken into consideration.This has been neglected in the simulation so far.This coupling will, however, be investigated in the future.
The matrix Σ is defined as Here, t is the pressure difference which is applied to the system [42].The equations of motion above have the conserved quantity where G = h t h is the metric tensor.The sixth term on the right side of the above equation denotes the elastic energy due to the external pressure [42].The Jacobian of the coordinate transform is computed according to [43].
Thus, the partition function becomes [42] with If M = 1 is set in the equations of motion, a single Nosé-Hoover thermostat is obtained.For the simulations here, M = 1 is set.
The time integration schemes closely follow the time-reversible measure-preserving Verlet [45,46] and rRESPA [47] integrators derived by Tuckerman et al. [37,48].The form according to [42] is used as a factorizing schema of the time evolution operator Here, iL is the Liouville operator which is split into three components [42]: with

Figure 2 .
Figure 2. Schematic illustration of the particle motion.

Figure 5 .
Figure 5. Structure of the model gel fragment (cluster of 1,000 particles) before and after shrinkage.

Figure 6 .
Figure 6.Experimental data of the fractal dimension vs. time after gelation point[22]

Figure 8 .
Figure 8. Computed consolidation ratio (MD-simulation) of the particulate network and experimental data from [23] as a function of time.Thin continuous curves: experimental data; dotted curves: Terzaghi model; thick continuous curve: MD model.