Compressive consolidation of strongly aggregated particle gels
Abstract
The compressive yield stress of particle gels shows a highly nonlinear dependence on the packing fraction. We have studied continuous compression processes, and discussed the packing fraction dependence with the particle scale rearrangements. The 2D simulation of uniaxial compression was applied to fractal networks, and the required compressive stresses were evaluated for a wide range of packing fractions that approached close packing. The compression acts to reduce the size of the characteristic structural entities (i.e. the correlation length of the structure). We observed three stages of compression: (I) elasticdominant regime; (II) singlemode plastic regime, where the network strengths are determined by the typical length scale and the rolling mode; and (III) multimode plastic regime, where sliding mode and connection breaks are important. We also investigated the way of losing the fractal correlation under compression. It turns out that both fractal dimension and correlation length start to change from the early stage of compression, which is different from the usual assumption in theoretical models.
I Introduction
Particle gels form in suspensions if shortrange attractions act between contacting surfaces of particles. Depending on the strength of the attraction, one can distinguish them as two types: strongly and weakly aggregated particle gels (Larson, 1999). In strongly aggregated particle gels, on which we focus in this work, particles once joined hardly separate at room temperature (Witten and Pincus, 2004). The attractions also prevent tangential displacements of contacting particles (Dinsmore, 2010), thus particles form samplespanning loose networks. Their remarkable and complex mechanical properties under external stress are essentially due to the response of a deformable disordered network formed with all the particles.
Under shear stress, one can sketch out deformation of a particle gel as follows (Nguyen and Boger, 1992):

the system is elastic under very low stresses. If the stress is released, the system comes back to the original shape;

for larger stress, a particle gel may undergo finite deformations, though it does not flow. This is the plastic regime, and it gives malleability; i.e., the system mostly keeps the deformed shape if stress is released;

beyond a yield stress, the particles network is ruptured and the particle gel flows.
Fluids exhibiting these behavior are called yield stress (or Bingham plastic) fluids (Nguyen and Boger, 1992; Denn and Bonn, 2011). Since the yield stress is the highest stress that can prevent flow, it can be considered as the strength of particle gels. Naturally, this strength relates to the strength of the particle connections. In a particle gel under shear stress, transmitting stresses break a portion of particle connections and the sample spanning network is torn apart. The yield stress of particle gels is higher when the packing fraction of particles is larger. In many cases, this variation appears to follow a power law with exponents between 4 and 5 (Channell and Zukoski, 1997).
Under compressive stress, the behavior of a particle gel looks different, due to confinement (de Kretser, Boger, and Scales, 2003; Stickland and Buscall, 2009). Compressive stress appears naturally in a number of applications, for example paste drying (Brown, Zukoski, and White, 2002) or ceramic processing (Lewis, 2000). Osmotic compression is a generic important case. In this kind of compression process, liquid removal tends to pack the particles into a smaller volume. As long as the applied stress is lower than the strength of the gel, it will retain the volume. If the stress exceeds a compressive yield stress , network restructuring leads to the loss of a part of the liquid phase, and particles reorganize into a smaller volume (Buscall et al., 1987). The deformation stops when the mechanical resistance of the network equals the stress. Thus, unlike the shear case, deformation does not continue with a constant applied stress (Buscall and White, 1987). The growth of mechanical resistance, i.e. compressive consolidation, also shows a highly nonlinear dependence on the packing fraction, typically powerlaw relations with exponents between 3.2 and 4.5 (Miller, Melant, and Zukoski, 1996; Green and Boger, 1997; Gisler, Ball, and Weitz, 1999; Madeline et al., 2007; Parneix et al., 2009).
The highly nonlinear dependency of the compressive yield stress on the packing fraction relates to the heterogeneity of the network structure of particles (Buscall et al., 1987). The fractal model (Brown, 1986; Shih et al., 1990; Potanin and Russel, 1996; Mewis and Wagner, 2011) successfully captures such packingfraction dependence of rheological proprieties by considering the characteristic structure of strongly aggregated particle gels. In the model, a particle gel is assumed to be an assembly of fractal flocs—such structural correlation is evidenced by scattering observations (Dietler et al., 1986; Poon and Haw, 1997). The entire sample is filled up with flocs of a fractal dimension , and the correlation length is equivalent or proportional to the average size of the flocs. Thus, the correlation length can be associated with the packing fraction by using the fractal dimension : , where and are the spatial dimension and proportional coefficient. Under this assumption, a gel of a higher packing fraction consists of a higher number of smaller flocs of the same fractal dimension . By deducing average mechanical responses of the flocs, the fractal model predicts a powerlaw dependence of the shear modulus on with an exponent , where is the fractal dimension of the backbone structure. Although this theoretical prediction has often been cited in the literature, we should pay attention to the scope of the model. Since fractal correlation yields in the formation process of particle gels, there is no clear reason for and to be kept constant during the compression process. Actually, some scattering observations indicate that the compression process changes the local structure (fractal dimension) (Madeline et al., 2007; Parneix et al., 2009). To make progress in macroscopic modeling of aggregated colloidal or particle gels (Flatt and Bowen, 2006; Lester, Rudman, and Scales, 2010), we need to settle such microstructural issues.
There is a recent twodimensional(2D) simulation work to support the application of the fractal model to compression processes (Gilabert, Roux, and Castellanos, 2008). Gels which were (i) directly prepared to the packing fraction , or (ii) prepared by compression from a lower packingfraction gel showed the same structural characteristics. However, the investigated range of packing fractions was relatively high (). Since the fractal model targets on strongly aggregated gels, we should compare simulations at lower packing fractions. Indeed, our 2D simulation investigates a wider range of packing fractions ().
This article reports results of numerical simulation for compression of particulate networks to discuss the consolidation behavior with restructuring. Experimentally observed time evolution of compressed gels should be affected by both the hydrodynamic and network resistances (Manley et al., 2005b; Kim et al., 2007; Brambilla et al., 2011). But, we here focus on studying the latter; our simulation work investigates the quasistatic structural evolution determined by the balance between external stresses and network strength at long time. This limitation to quasistatic aspects is common with the abovementioned fractal model. So, compression considered in this work is sufficiently slow, such as osmotic compression. A contact model is employed to simulate sticky and breakable connections between contacting particles. In granular physics, there are many prior studies to investigate cohesive granular system (Iwashita and Oda, 1998; Tomas, 2007). Their simulation methods are an extension of the discrete element method (Cundall and Strack, 1979); the rolling friction has been additionally included into the contact model. Though our target is a suspension system, the spirit of the modeling is almost the same. In granular physics, there are also some studies on compression of cohesive powders (Kadau et al., 2002; Bartels et al., 2005; Wolf et al., 2005; Gilabert, Roux, and Castellanos, 2007, 2008), which obtained many common insights with our study. 2D simulation results are reported here. One obvious benefit of studying 2D systems is easy recognition of the structure changes. Though the densitydensity correlation is an effective tool to capture special sorts of heterogeneous structures, it might not capture all features of random structures. 2D experimental system has also been studied to observe the microstructure change under shear (Masschaele, Fransaer, and Vermant, 2009). Though we do not expect qualitative difference between 2D and 3D systems, this point should be checked by 3D simulation of the same modeling, which we hope to report in future work.
Ii Method
ii.1 Initial configuration: fractal gel
We use fractal aggregation models to generate the initial configurations. Networks of particles are formed in a rectangle box of sides and with periodic boundary conditions. We combine two processes. The first step models a formation process of fractal clusters in the dilute limit; the reactionlimited hierarchical clustercluster aggregation model was employed (Jullien and Botet, 1987). Each cluster consists of particles. The second step is to build a network with them. of the generated clusters are randomly placed in the box while avoiding overlap with each other, and then they are moved with a random walk. Collisions of clusters always lead to coalescence, so a single cluster is formed in the periodic boundary space after a while (see, e.g., Figure 7 (a), below). Thus, we can generate random networks of a packing fraction , where is the radius of a particle.
One can deduce structural features of the formed networks by evaluating the densitydensity correlation function (See appendix A). We prepared 10 samples of particles in a box . The packing fraction is . Figure 1 shows the averages and standard deviations of their correlation functions . A powerlaw decay of correlation is seen in the finite range . The slope indicates fractal dimension , which is almost the same fractal dimension as the clusters generated by the aggregation model. Though a crossover behavior is seen in the range , the correlation decays out after that.
ii.2 Contact forces
Our contact model is designed to simulate strongly aggregated particle gels. The contact model used in the original discrete element method (DEM) (Cundall and Strack, 1979) takes into account friction forces for the sliding mode, which is essential to simulate dense granular systems. For cohesive granular systems, the rolling friction has been added (Iwashita and Oda, 1998). This extension is also required to reproduce the behavior of particle gels having chainlike local structures. Our implementation is similar to the model used for cohesive granular systems. Since the dynamic friction, which is modeled in typical contact models, is unclear for colloidal systems, this part is different in our model. Though a brief summary will be given below, the detailed descriptions of our implementation were given in another paper (Seto et al., 2012).
We assume that interactions acting among particles are only short range forces. No interparticle interaction acts between two particles before getting into contact. Once the gap between two particles becomes zero, a cohesive force starts to act between them. For simplicity, our simulation model assumes immediate bond generation, which is expected to give the upper limit of the consolidation behavior. When compressionrate dependence is problem, the condition of bond generation can be crucial. Finite compression rates can suppress bond generation, leading to a weaker consolidation.
The forces and moments are assumed to be Hookian: Let us consider the contact interactions between particles and . The centertocenter distance between them is , and the gap . The normal (separation/overlap) displacement is related to the normal element of the contact force,
(1) 
where for and otherwise. is the normal unit vector defined by . The second term of eq. (1) is added to avoid significant overlap between particles. The sliding displacement relates to the tangential part of the contact force,
(2) 
The sliding displacement is the tangential projection of the deviation vector between the original contact points and ,
(3) 
where the vectors and indicate the original contact points from the respective centers of the particles and are fixed with them. The bending angle ( is related to the moment acting on the rolling mode,
(4) 
These elastic relations are broken when the force or moment exceed threshold values. We assume simple criteria given with three threshold values: the critical normal force , critical sliding force , and critical rolling moment . When the pulling force exceeds the value; , the bond is broken and the two particles separate. When the tangential force or rolling torque exceeds the threshold value; or , both the sliding displacement and rolling angle are reset to zero, which releases the tangential parts of the stored energy. The maximum strains are related to the threshold forces and moment as follows: , , and .
It must be noted that in our simple model the bond strength of the tangential mode does not depend on the normal load. However, the enlarged contact area due to the normal load may change the tangential strengths (Dominik and Tielens, 1995). Our neglect of this effect can affect simulations at higher packing fractions.
ii.3 Equation of motion
A particle may be in contact with other particles . The total contact force and torque acting on the particle can be written as follows:
(5) 
Particles are immersed in a viscous fluid, so that hydrodynamic interactions also act on them. For micro or nanosized particles, the particle Reynolds number is extremely small, so the Stokes equations should be used to evaluate the hydrodynamic interactions among many particles (Kim and Karrila, 1991). However, we do not go into the detailed hydrodynamics here. Hydrodynamic interactions in the Stokes regime are scaled with the velocities of particles. If the slow compression limit is considered, flows induced by particle movements are negligible in comparison with the contact forces. They just work as energy dissipators. Indeed, some experimental observations suggest that a solvent flow does not contribute the rearrangements of particles in a slow sedimentation of colloidal gels (Brambilla et al., 2011).
As an approximation, the Stokes formula is used in this work; the force and torque acting on a particle are proportional to the translational velocity and angular velocity , respectively:
(6) 
where is the viscosity of the solvent.
Since the interparticle interactions are assumed to be sufficiently stronger than thermal agitations, Brownian forces are also neglected in this work. Consequently, Newton’s equations of motions are given as the sum of the contact forces (5) and hydrodynamic interactions (6),
(7) 
Since particle inertia is very small for colloidal systems, it is a good approximation to neglect the inertia terms (the left hand sides) and solve the force balance equations (Ball and Melrose, 1997). As explained later (sect. II.5), however, we took a different strategy to find the static equilibriums.
ii.4 Uniaxial compression and static equilibrium
The compression is externally controlled. For simplicity, we focus on examining the uniaxial compression along the axis in this work. The periodic boundaries at and of the simulation box introduced in section II.1 are replaced by two flat walls. These walls can move along the axis like pistons, and they are assumed to be semipermeable; liquid passes through them, while particles do not. When particles get in contact with the walls, they are rigidly fixed there; particles contacting with the wall move together as one object. The periodic boundaries along the axis are still used during the compression simulation.
The stresscontrolled compression, in which an external stress is applied to cause the compression, is similar to typical experimental situations, such as osmotic compression. However, due to the higher efficiency of the simulation, we use straincontrolled compression. This also allows us to have better statistics in our data. This work investigates the slow compression limit; the time scale of compression is much longer than the relaxation time of particle motions. In this case, the results obtained with the straincontrolled compression are essentially equivalent to the ones with the stresscontrolled compression.
The compressive stress balancing the network resistance is evaluated from the total forces acting on the walls from particles. Since particles directly connecting to the walls (called wall particles) no longer have any mechanical freedom, the forces between wall particles and particles contacting with them are summed up as follows:
(8) 
When we refer to “compressive stress” in this 2D simulation work, our simulation box is assumed to have a thickness of particle diameter . Therefore, the compressive stresses are given as and .
Our compression simulations are performed in a stepwise manner in order to obtain numbers of equilibriums by a single simulation run. Compression steps are given by two increment ratios and ; the initial compression from to is given as , and the th packing fraction as . The increment ratio is determined with the number of steps from to : . The straincontrolled compression is simply realized by translations of the walls with the velocity , so that the packing fraction increases continuously.
When it reaches a target packing fraction , the translations of the walls are suspended to await the equilibrium. Usually, the excess stress due to the quickness of the compression is relaxed by spreading the deformations into the entire system. The following criteria are examined every time steps of the simulation:

The relative difference between the top and bottom stresses is smaller than a threshold value : , where is the average, .

The relative change of the compressive stress for the last steps is smaller than a threshold value : , where is the one of steps before.

No bonds were broken for the last steps.
Thus, the compressive stress balancing the network of packing fraction can be determined. Once is determined, the compression continues again for the next target .
ii.5 Non realtime simulation for quasistatic simulation
In order to reduce finite size effects, the simulation box needs to involve a sufficient number of typical structures. The compressive stress to balance sparse networks is very low, and the most of particles are slightly stressed. In this case, the relaxation time is significantly long. Usually, the overdamped equations of motion, which are given by dropping the right hand sides of eqs. (7) due to the negligible inertia of particles, are used to simulate colloidal suspensions. However, we have taken a different approach for efficiency. As restricting our scope in the slow compression, we may be allowed to optimize some parameters for a faster simulation. In the slow compression limit, velocities of particles are also slow. In this case, the contact forces dominate the displacements of particles. This is why a reduced viscosity still leads to the same equilibrium configurations with a shorter time. The following viscosity gives the critical damped motion in terms of the rolling mode.
(9) 
Though it is not simple problem to find the best value, we used in the equation of motion (7). This arbitrary selection of the viscosity is justified only if the compression is considered slow enough. This point will be checked in sect. III.3.
Iii Results and discussion
iii.1 Parameters for contact model
The introduced 2D contact model (see sect. II.2) is relatively simple, but it still has 7 parameters to be given. Traditional experimental measurements to characterize the particleparticle adhesion determine the pulloff force (Larsen, 1958), which corresponds to in this work. In recent years the atomic force microscope became a reliable method, which provides more precise data for normal forcedisplacement relationships (Hodges et al., 2002). However, there are fewer experimental measurements to determine the parameters of the other modes (Heim et al., 1999; Ecke et al., 2001). The experimental measurement of the threepoint bending test by using optical tweezers is remarkable (Pantina and Furst, 2005; Furst and Pantina, 2007), which provides parameters for the rolling mode, and (or ). In any case, all parameters for a single system have never been determined so far. Thus our knowledge of the contact forces is still limited, so that we need to set the parameters by hand. But, we may have a certain level of confidence by reproducing a similar bending behavior with the threepoint bending measurements. Our selected parameters are intended to imitate the observed mechanical behavior: a small extent of elastic deformation and irreversible reorganization due to a critical moment.
We have selected two sets of the parameters called bond 1 and bond 2, which are shown in Table 1. Bond 1 has the same strengths for the three modes, while bond 2 is 5 times stronger for the normal and sliding mode than the rolling mode. Both have the same maximum strains; elastic deformations give 5% of the radius.
bond 1  bond 2  

Normal strength  1  5  
Sliding strength  1  5  
Rolling strength  1  1  
Normal strain  0.05  0.05  
Sliding strain  0.05  0.05  
Rolling strain  0.05  0.05  
Parameter in eq. (1) 
The simulation for the threepoint bending test consists of the following three steps:

Prepare a linear chain under no stress.

Apply the forces at the central particle and at both the ends, and wait for equilibrium.

Stop applying the forces, and wait for equilibrium.
The three snapshots i–iii in Figure 2 (a) and (b) show the configurations of the respective equilibriums of the bond 1 simulation. Though the data is not shown here, the results of bond 2 are similar. The two results were obtained with two slightly different forces: (a) , (b) . If the external force did not reach the threshold, the initial linear shape remained after the test (see (a) iii). On the other hand, if the bond stress exceeded the threshold value, some irreversible rearrangement took place. As a result, the initial shape was not recovered (see (b) iii). In the latter case, the bending deformation proceeds until the moment acting on the broken bond becomes smaller than the threshold due to the increase of the bending angle. This is why the small increment of the applied forces leads to the visible difference of the deformation (see ii of (a) and (b) in Figure 2). The slight bending deformation seen by (a) ii is close to the limit of elastic deformation of this cluster.
Thus, the parameters were chosen to simulate some brittle behavior of colloidal aggregates; the local structures remain unless bond breakups take place. In this test, the rolling mode of the contact model is responsible for the cluster’s deformation.
iii.2 Parameters of the compression simulation
The initial configurations (sect. II.1) are fractallike networks with a packing fraction , where and for bond 1 simulation. Since the calculation cost of the bond 2 simulation is higher, smaller system size, and , is simulated. The first target packing fraction is set to . The increment ratio and the number of compression steps are and , respectively. The compression finally reaches to . The relaxation time depends on the parameters of the bonds. To account for this, the wall velocities are set to and for the bond 1 and bond 2 simulations, respectively, where . In the equilibrium checks (sect. II.4) after every time integration steps, is used for both the simulations, and and for bond 1 and for bond 2 simulations, respectively.
iii.3 Uniformity of compression
As explained in sect. II.5, it is useful to introduce artificial parameters to speed up the compression simulation. Therefore we started to check that our method reproduced the quasistatic compression as expected. If compression is too fast, the regions near the wall are easily crushed. So, we examined the uniformity of the deformations along the compression axis. The uniformity of compression can be seen via the dependence of the packing fraction , which is defined as the packing fraction in the sliced ranges: (where was chosen). Due to the heterogeneous structure and finite system size, the uniformity cannot be judged with a single compression simulation. This is why we took the averages over 10 simulation runs with different initial configurations.
Figure 3 shows the results for the bond 1 simulations, where every second compression step is plotted. Although we can see an undesirable drop at the middle of the first compression (), the compression is satisfactorily uniform along axis after that. The steep drops at both the edges relate to the integrity of the walls; particles connecting to the walls are not allowed to move (see section II.4). Equivalent results were also obtained with the bond 2 simulation with a slower compaction speed (see sect. III.2).
iii.4 Data of compressive consolidation
Once the system reaches a static equilibrium under a compressive stress , further compression requires a finite increment: . Therefore, the balancing stress can be considered as the strength of the gel network, i.e., the compressive yield stress (Buscall and White, 1987). A simulation run of the stepwise compression (sect. II.5) finds a series of vs. relations. Figure 4 plots the results of the bond 1 and bond 2 simulations, where the unit of the stress is . The compression from to shows wide ranges of (about four decades). These results indicate significantly rapid compressivestrain hardening.
Simulation provides various information, which allow us to understand the compression. Some of data indicate that there are three regimes in the compression process, which will be seen in the following section.
Bond creation
Bond creation is essential to understand the irreversible consolidation of the compression. Due to heterogeneous structure of the networks, the compressive deformations are not affine. Consequently, two separated particles may come in contact and newly bond with each other. The mean contact number per particle reflects new bond creations. The result of bond 1 simulation is shown in Figure 5 (a). The compression continuously changes the value from 2 to about 4 with increasing rate, so we do not see distinct regimes in this.
Stored and dissipated energies
The bonds have finite strength, and they can be broken up or irreversibly reorganized by certain applied stress. Such events cause the plastic deformation of the networks, which contribute the irreversibility. The translations of the wall by external stress inject energy into the system. Part of the energy are stored in contact bonds first and dissipated as bond ruptures. We can see this energy flow in the successive increments of the total stored and dissipated energies for each compression step (Figure 5 (b)). The total dissipated energy exceeds the total stored energy at the early stage, . Accordingly, we can identify the boundary of two regimes: elastic dominant regime (I) and plastic regime (II).
Bond rupture
The accumulated numbers of the bond ruptures were recorded as distinguishing the stress modes causing the breakup, where the indicates ‘connection break’, ‘sliding’, or ‘rolling’. The rupture rates can be defined as the number of breakups normalized with the total particle number and relative compressive strain:
(10) 
Figure 6 shows the results of both the bond 1 and bond 2 simulations. At lower packing fractions, only rupture events occur due to rolling. This indicates that compression of loose networks is achieved by bending deformations of particulate chains. When the compaction reaches a certain point, further compression requires sliding or separation of contacting particles. The onset of this mixedrupture mode depends on the parameters of the contact forces. The connection breaks start between and in bond 1 simulations, and between and in bond 2 simulations. They split the plastic deformation into two parts: singlemode plastic regime (II) below and multimode plastic regime (III) above.
iii.5 Three stages of compression
The date shown in the previous section indicate that three distinct behaviors can be recognized in the compression process. Let us discuss them one by one here.
(I) Elastic dominant regime
There is no percolation path in the 10 initial configurations of (see Figure 7 (a)). Hence, no compressive stress is required to compress them in the slow compression limit. After some compaction, the first spanning path appears, and the network begins to resist the compressive stress. At the gel point, a singlelinking chain (percolation path) resists the external stress between top and bottom. This systemspanning cluster is deformed as a whole, so that strains of single bonds can be very small. Therefore, the dominant response to the external stress is elastic deformation; the work by compression is stored in the contact bonds (see Figure 5 (b)).
In this range, the powerlaw form expanding at the gel point is suitable to capture the  curve (de Gennes, 1976):
(11) 
In Figure 4, the dotted curve shows the fitting result with the function (11) for the bond 1 simulation data of . The fitting gives the extrapolated gel point as . The evaluated exponent represents some nonlinearity. The compression induces some new contacts, so that the numbers of percolation paths are increased. The rise of resistance due to such increasing paths can be a plausible interpretation of the nonlinearity in this elasticdominant regime.
It must be noted that our simulation results in this range might have large statistical fluctuations. The typical length scale of the structure is comparable to the system size (compare the drawn circle and the box size in Figure 7 (a)), and consequently only a few paths enter the averaging. Therefore, the gel point obtained in the above fitting has only limited significance.
(II) Singlemode plastic regime
By considering two data (Figure 5 (b) and Figure 6), we identified the singlemode plastic regime in the ranges: for bond 1 simulations and for bond 2 simulations.
This dependence of the consolidation behavior has been discussed with a scaling concept (Buscall et al., 1988; Shih et al., 1990). The powerlaw relationships of and have been also observed in experiments (Buscall et al., 1987; Madeline et al., 2007; Parneix et al., 2009):
(12) 
A similar powerlaw behavior can be seen in our simulation results (Figure 4). The least squares fitting gives the almost same exponents: and , for the bond 1 and bond 2 results, respectively.
Let us look more closely at this regime. Figure 7 (b) shows the configuration of particles at the beginning of this range () in a bond 1 simulation. At the gel point, a single linking path spanned the entire box. In contrast, multilinking structures span the system here. Since the maximum bending moment of a randomwalk chain is proportional to the inverse of the chain size, the supportable loads of larger loops are smaller. The variety of the mesh sizes indicates a heterogeneity of the local strengths. Though large empty spaces remain, the relatively dense meshes fill the space from the bottom to the top (Figure 7 (b)). Further compression requires destroying a part of them. As shown in Figure 6, the rollingmode ruptures are the only rupture events seen in this range. Therefore, the size of the responsible mesh structures determines the compressive yield stress .
Thus, in the singlemode plastic regime, the length scale is essential; the compression breaks larger and weaker structures into smaller and more robust structures. Fractal model was employed to give the length scale from the packing fraction (Brown, 1986; Shih et al., 1990) (see sect. III.6). According to the nature of a fractal, a portion of a fractal cluster is also a small fractal cluster. If the compression of fractal networks causes only fragmentation of fractal clusters into smaller pieces, it can be expected that the correlation length is decreased while keeping the same fractal dimension. This optimistic expectation leads us to employ the fractal model for the compression problem. We will check this point in the next section.
(III) Multimode plastic regime
The onset of connection breaks indicates that the compression enters a new regime. This transition point depends on the parameters of the bonds. We consider the data above and in bond 1 and 2 simulations as the multimode plastic regime. In these ranges, the compression curves look to follow a power law. However, the evaluated exponents are totally different with each other: and . It is by chance that the exponent of bond 2 simulations is similar with the one in the singlemode plastic regime. The exponent in this regime is easily affected by the parameters of the contact model, that implies absence of a simple scaling law.
Figure 7 (c) depicts the beginning of this regime of a bond 1 simulation (). There are still large loops remaining, but they look like voids due to the denser surroundings. Here, we can easily find dense meshes and trimers (a trimer can be considered as the building block of a close packing). If they are next to each other, the entire object can be considered as a lump. The same scenario of the singlemode plastic regime seems to work; denser mesh domains are piled up over the entire system. However, due to the finite size of particles, the mechanical response of denser meshes is not the same as looser ones. The crush of them involves sliding ruptures and connection breaks, which are seen in Figure 6. Since the criteria of destruction includes the all three modes, the compression curve depends on the details of the contact force.
This regime terminates due to the close packing, where rigidity of particles prevents reorganizations of particles. Since our simplified contact model is not accurate for that, we cannot specify the upper end of the range.
iii.6 Fractal correlation under compression
During compression, separated particles can newly come into contact. Consequently, mesh size decreases as is increased by compression. If we employ the fractal model (Brown, 1986; Shih et al., 1990) to explain the dependence of , the structure change needs to follow a special manner; the fractal dimension characterizing the local structures does not depend on the packing fraction , and the correlation length is a decreasing function of the packing fraction , .
The fractal dimension can be estimated from the densitydensity correlation function (see Appendix A). Figure 8 shows averages of the correlation functions of every second compression step, (), of the bond 1 simulation results. The results of the bond 2 simulations were almost the same. The correlation of the density is mostly decreasing functions of the twopoint distance . The correlation of the configuration before compression is the largest and reaches the farthest distance (solid line). The compression reduces this correlation in a continuous manner.
The fractal correlation appears as a straight slope in the loglog plot of , i.e., . Though the straight lines can be recognized in the initial configuration, the compressed networks show slightly curved lines. In order to make them easier to see, we introduce the fractal dimension profile from the local slope of the correlation function at each distance :
(13) 
The results are shown in Figure 9 (a). As mentioned in sect. II.1, the initial configurations (solid line) consist of three parts: a fractal plateau () indicating , crossover range (), and 2Dlike uniform network (). Thus, though the recognizable fractal plateau is limited, the profiles of the initial configurations agree with the expected behavior as fractal networks.
Once the networks are compressed, the fractal correlations are affected immediately. Indeed, the first compression step of our simulation () shows the entire change of : a narrower plateau indicating a rise of the fractal dimension, and shifts of the crossover and uniform ranges to shorter lengths. These results suggest that not only longdistance correlations, but shortdistance correlations can be affected by the compressive deformation at the early stage. This tendency continues for the further compression.
Though the profiles for the compressed networks indicate that monofractal assumption does not work, it may be worth determining the dependence of the correlation length . Since the slopes of in Figure 8 are not constant, it is ambiguous to determine from them. In our purpose, we may consider the correlation length as the length scale where the network looks uniform. The typical length scale for the uniformity can be determined by evaluating fluctuation of local densities (Masschaele, Fransaer, and Vermant, 2009). The local densities are defined by boxes of size , which are denoted by . If the size of box is larger than the correlation length , the local densities show little fluctuation; the system looks uniform. Thus, the standard deviation of the local densities normalized by the global density, , measures the uniformity. Due to the finite system size, we compare here the lengths showing a common finite fluctuation, . The actual values of the correlation lengths should be larger than the evaluated lengths , so a common arbitrary factor may be given to adjust. Figure 9 (b) shows the result of the lengths satisfying . We find a powerlaw relation between and with an exponent . The obtained exponent is a little bit smaller than the fractal assumption: . The evaluated correlation lengths, , are also shown as the diameter of the circles drawn in Figure 7.
Iv Conclusion
We focused on investigating compression of particle gels. Inspired by scaling theory for polymer gels, the fractal model has been introduced to explain the packing fraction dependence of particle gels’ rheology (Brown, 1986; Shih et al., 1990), and referred for the compression process (Buscall et al., 1988). The model assumes that the structure of a particle gel is fractal network. However, the simple question whether a particle gel under compression is fractal network or not, has left unanswered. We have aimed to answer this question by introducing a particlescale simulation model in the limit of strongly aggregated particle gels. In this case, the cohesive force among contacting particles is assumed to be so strong that Brownian forces do not affect the structural evolution of gels.
Our simulation reproduced the reasonable consolidation behaviors as vs. relationships in the case of quasistatic compression. By monitoring various particlescale information, i.e., mean contact number, stored/dissipated energy, and number of bond ruptures, we identified the three distinct regimes: (I) elastic dominant regime, (II) singlemode plastic regime, and (III) multimode plastic regime. At the gel point, there are many “dead end” branches, where the particles in such branches do not contribute to the strength of the network . Compression makes separated particles come into contact, so that the mesh size becomes smaller and smaller. This explanation is already known and qualitatively true for the most part. Such shrinkage of the typical length scale causes a transition of the local mechanical response from (II) to (III) in the plastic regimes. While the length scale is large enough, the rolling mode is the responsible for the deformation. At a certain length scale, the sliding and normal modes become more important for further compression. Thus, we confirmed that a single regime could never cover the wide range of packing fractions in spite of subtle appearance in the vs. relationships. The simulation results also provided an insight that heterogeneity in the local mesh sizes plays some roles in the way of reorganization under compression and overall network strength. Mechanical correlations (connectivity) determine which parts of the network are destroyed under an applied stress. If a percolation path consisting of robust domains is built, they support a large fraction of the applied stress. In this case, the weakest domain in the path is most probable to be destroyed. This restructuring yields a smaller and robuster domain, which results in the reinforcement of the percolation path. However, such mechanism is subtle because compression is not continuing deformation and tends to reduce the heterogeneity.
In order to see whether the fractal assumption remains or not in particle gels under compression, we evaluated the densitydensity correlation for the simulation results. The evaluated fractal dimension profile says that compression affects both the local structure and correlation length in the early stage. Thus, our results are contradictory to the application of the fractal model to compression processes. However, the fractal dimension profiles show a systematic change in the singlemode plastic regime. This result seems to indicate required ingredients to upgrade the fractal model for compression processes.
The simulation model of this work is simple; we neglect various factors in real particle gels to make clear the focus of study and to optimize the simulation speed in order to achieve large system sizes. For example, we employed a simple contact model, which depends on only configuration of particles and has no history dependence. However, cohesion between particles in real system may be more complex, e.g. the bonding or breakup may require certain time periods. In order to deal with such systems, our quasistatic simulation framework needs to be extended. This extension is relevant for the future work because there are many experimental works investigating on the time/ratedepending aspects of colloidal gels and yield stress fluids (Manley et al., 2005a; Bonn and Denn, 2009; Brambilla et al., 2011).
Acknowledgements.
The authors would like to thank Prof. Morton Denn for valuable suggestions on the manuscript, and RS would like to express his gratitude to Prof. Richard Buscall for fruitful discussions, valuable suggestions, and warm encouragement.Appendix A Densitydensity correlation function
In order to study colloidal gels, we need to have some tools to characterize random structured systems. The densitydensity correlation function is used for this purpose. The local density function is defined as follows:
(14) 
By taking the average over the entire space of the system, the densitydensity correlation function is defined as follows:
(15) 
When the system is isotropic, the correlation is given in terms of the distance, , which can be obtained by taking averages over all directions:
(16) 
The correlation function can be regarded as a conditional probability to find particles. Since it is normalized by the entire probability, the value 1 means no correlation.
It is known that strongly aggregated gels have structures involving fractal correlations. The local structures are equivalent to fractal flocs. For fractal flocs, the relation between the number of particles and the radius of gyration follows a powerlaw:
(17) 
The exponent is called fractal dimension. Therefore, the average packingfraction profile of the clusters is given as follows:
(18) 
where is the distance from the center of mass. Since fractal gels can be considered as clusterfilling networks, the local structure should coincide with this density profile (18). Consequently, the densitydensity correlation function has the same functional form: . In a fractal gel, the correlation reaches a finite distance, which is called the correlation length . For , the correlation does not remain:
(19) 
If the local densities are evaluated in boxes of size , the system should look uniform.
References
 Ball and Melrose (1997) Ball, R. C. and Melrose, J. R., “A simulation technique for many spheres in quasistatic motion under frameinvariant pair drag and brownian forces,” Physica A 247, 444–472 (1997).
 Bartels et al. (2005) Bartels, G., Unger, T., Kadau, D., Wolf, D. E., and Kertész, J., “The effect of contact torques on porosity of cohesive powders,” Granular Matter 7, 139–143 (2005).
 Bonn and Denn (2009) Bonn, D. and Denn, M. M., “Yield stress fluids slowly yield to analysis,” Science 324, 1401–1402 (2009).
 Brambilla et al. (2011) Brambilla, G., Buzzaccaro, S., Piazza, R., Berthier, L., and Cipelletti, L., “Highly nonlinear dynamics in a slowly sedimenting colloidal gel,” Phys. Rev. Lett. 106, 118302 (2011).
 Brown, Zukoski, and White (2002) Brown, L. A., Zukoski, C. F., and White, L. R., “Consolidation during drying of aggregated suspensions,” AIChE Journal 48, 492–502 (2002).
 Brown (1986) Brown, W. D., Ph.d thesis, Department of Physics, Cambridge University (1986).
 Buscall et al. (1988) Buscall, R., Mills, P. D. A., Goodwin, J. W., and Lawson, D. W., “Scaling behaviour of the rheology of aggregate networks formed from colloidal particles,” J. Chem. Soc. Faraday Trans. 84, 4249–4260 (1988).
 Buscall et al. (1987) Buscall, R., Mills, P. D. A., Stewart, R. F., Sutton, D., White, L. R., and Yates, G. E., “The rheology of stronglyflocculated suspensions,” J. NonNewtonian Fluid Mech. 24, 183–202 (1987).
 Buscall and White (1987) Buscall, R. and White, L. R., “The consolidation of concentrated suspensions: Part 1. the theory of sedimentation,” J. Chem. Soc. Faraday Trans 83, 873–891 (1987).
 Channell and Zukoski (1997) Channell, G. M. and Zukoski, C. F., “Shear and compressive rheology of aggregated alumina suspensions,” AIChE J. 43, 1700–1708 (1997).
 Cundall and Strack (1979) Cundall, P. A. and Strack, O. D. L., “A discrete numerical model for granular assemblies,” Geotechnique 29, 47–65 (1979).
 de Gennes (1976) de Gennes, P.G., “On a relation between percolation theory and the elasticity of gels,” J. Physique Lett. 37, 1–2 (1976).
 Denn and Bonn (2011) Denn, M. M. and Bonn, D., “Issues in the flow of yieldstress liquids,” Rheol Acta 50, 307–315 (2011).
 Dietler et al. (1986) Dietler, G., Aubert, C., Cannell, D. S., and Wiltzius, P., “Gelation of colloidal silica,” Phys. Rev. Lett. 57, 3117 (1986).
 Dinsmore (2010) Dinsmore, A. D., “Soft random solids: particulate gels, compressed emulsions, and hybrid materials,” in Experimental and computational techniques in soft condensed matter physics, edited by J. Olafsen (Cambridge, 2010) Chap. 3.
 Dominik and Tielens (1995) Dominik, C. and Tielens, A. G. G. M., “Resistance to rolling in the adhesive contact of two elastic spheres,” Phil. Mag. A 72, 783–803 (1995).
 Ecke et al. (2001) Ecke, S., Raiteri, R., Bonaccurso, E., Reiner, C., Deiseroth, H.J., and Butt, H.J., “Measuring normal and friction forces acting on individual fine particles,” Rev. Sci. Instrum. 72, 4164–4170 (2001).
 Flatt and Bowen (2006) Flatt, R. J. and Bowen, P., “Yodel: a yield stress model for suspensions,” J. Am. Ceram. Soc. 84, 1244–1256 (2006).
 Furst and Pantina (2007) Furst, E. M. and Pantina, J. P., “Yielding in colloidal gels due to nonlinear microstructure bending mechanics,” Phys. Rev. E 75, 050402(R) (2007).
 Gilabert, Roux, and Castellanos (2007) Gilabert, F. A., Roux, J.N., and Castellanos, A., “Computer simulation of model cohesive powders: Influence of assembling procedure and contact laws on low consolidation states,” Phys. Rev. E 75, 011303 (2007).
 Gilabert, Roux, and Castellanos (2008) Gilabert, F. A., Roux, J.N., and Castellanos, A., “Computer simulation of model cohesive powders: plastic consolidation, structural changes, and elasticity under isotropic loads,” Phys. Rev. E 78, 031305 (2008).
 Gisler, Ball, and Weitz (1999) Gisler, T., Ball, R. C., and Weitz, D. A., “Strain hardening of fractal colloidal gels,” Phys. Rev. Lett. 82, 1064–1067 (1999).
 Green and Boger (1997) Green, M. D. and Boger, D. V., “Yielding of suspensions in compression,” Ind. Eng. Chem. Res. 36, 4984–4992 (1997).
 Heim et al. (1999) Heim, L.O., Blum, J., Preuss, M., and Butt, H.J., “Adhesion and friction forces between spherical micrometersized particles,” Phys. Rev. Lett. 83, 3328–3331 (1999).
 Hodges et al. (2002) Hodges, C. S., Cleaver, J. A. S., Ghadiri, M., Jones, R., and Pollock, H. M., “Forces between polystyrene particles in water using the afm: Pulloff force vs particle size,” Langmuir 18, 5741–5748 (2002).
 Iwashita and Oda (1998) Iwashita, K. and Oda, M., “Rolling resistance at contacts in simulation of shear band development by dem.” J. Eng. Mech. 124, 285–292 (1998).
 Jullien and Botet (1987) Jullien, R. and Botet, R., Aggregation and fractal aggregates (World Scientific, Singapore, 1987).
 Kadau et al. (2002) Kadau, D., Bartels, G., Brendel, L., and Wolf, D. E., “Contact dynamics simulations of compacting cohesive granular systems,” Comp. Phys. Commun. 147, 190–193 (2002).
 Kim et al. (2007) Kim, C., Liu, Y., Kühnle, A., Hess, S., Viereck, S., Danner, T., Mahadevan, L., and Weitz, D. A., “Gravitational stability of suspensions of attractive colloidal particles,” Phys. Rev. Lett. 99, 028303 (2007).
 Kim and Karrila (1991) Kim, S. and Karrila, S. J., Microhydrodynamics: Principles and selected applications (ButterworthHeinemann, Boston, 1991).
 de Kretser, Boger, and Scales (2003) de Kretser, R. G., Boger, D. V., and Scales, P. J., “Compressive rheology: an overview,” in Rheology Reviews 2003 (2003) pp. 125–165.
 Larsen (1958) Larsen, R. I., “The adhesion and removal of particles attached to air filter surfaces,” Am. Ind. Hyg. Assoc. J. 19, 265–270 (1958).
 Larson (1999) Larson, R. G., The structure and rheology of complex fluids (Oxford University Press, New York & Oxford, 1999).
 Lester, Rudman, and Scales (2010) Lester, D. R., Rudman, M., and Scales, P. J., “Macroscopic dynamics of flocculated colloidal suspensions,” Chem. Eng. Sci. 65, 6362–6378 (2010).
 Lewis (2000) Lewis, J. A., “Colloidal processing of ceramics,” J. Am. Ceram. Soc. 83, 2341–59 (2000).
 Madeline et al. (2007) Madeline, J. B., Meireles, M., Bourgegette, C., Botet, R., Schweins, R., and Cabane, B., “Restructuring of colloidal cakes during dewatering,” Langmuir 23, 1645–1658 (2007).
 Manley et al. (2005a) Manley, S., Davidovitch, B., Davies, N. R., Cipelletti, L., Bailey, A. E., Christianson, R. J., Gasser, U., Prasad, V., Segre, P. N., Doherty, M. P., Sankaran, S., Jankovsky, A. L., Shiley, B., Bowen, J., Eggers, J., Kurta, C., Lorik, T., and Weitz, D. A., “Timedependent strength of colloidal gels,” Phys. Rev. Lett. 95, 048302 (2005a).
 Manley et al. (2005b) Manley, S., Skotheim, J. M., Mahadevan, L., and Weitz, D. A., “Gravitational collapse of colloidal gels,” Phys. Rev. Lett. 94, 218302 (2005b).
 Masschaele, Fransaer, and Vermant (2009) Masschaele, K., Fransaer, J., and Vermant, J., “Direct visualization of yielding in model twodimensional colloidal gels subjected to shear flow,” J. Rheol. 53, 1437001460 (2009).
 Mewis and Wagner (2011) Mewis, J. and Wagner, N. J., Colloidal Suspension Rheology (Cambridge University Press, 2011).
 Miller, Melant, and Zukoski (1996) Miller, K. T., Melant, R. M., and Zukoski, C. F., “Comparison of the compressive yield response of aggregated suspensions: pressure filtration, centrifugation, and osmotic consolidation,” J. Am. Ceram. Soc. 79, 2545–2556 (1996).
 Nguyen and Boger (1992) Nguyen, Q. D. and Boger, D. V., “Measuring the flow properties of yield stress fluids,” Annu. Rev. Fluid Mech 24, 47–88 (1992).
 Pantina and Furst (2005) Pantina, J. P. and Furst, E. M., “Elasticity and critical bending moment of model colloidal aggregates,” Phys. Rev. Lett. 94, 138301 (2005).
 Parneix et al. (2009) Parneix, C., Persello, J., Schweins, R., and Cabane, B., “How do colloidal aggregates yield to compressive stress?” Langmuir 25, 4692–4707 (2009).
 Poon and Haw (1997) Poon, W. C. K. and Haw, M. D., “Mesoscopic structure formation in colloidal aggregation and gelation,” Adv. Colloid Interface Sci. 73, 71–126 (1997).
 Potanin and Russel (1996) Potanin, A. A. and Russel, W. B., “Fractal model of consolidation of weakly aggregated colloidal dispersions,” Phys. Rev. E 53, 3702 (1996).
 Seto et al. (2012) Seto, R., Botet, R., Auernhammer, G. K., and Briesen, H., ‘‘Restructuring of colloidal aggregates in shear flow: Coupling interparticle contact models with Stokesian dynamics,” Eur. Phys. J. E 35, 128 (2012).
 Shih et al. (1990) Shih, W.H., Shih, W. Y., Kim, S.I., Liu, J., and Aksay, I. A., “Scaling behavior of the elastic properties of colloidal gels,” Phys. Rev. A 42, 4772–4779 (1990).
 Stickland and Buscall (2009) Stickland, A. D. and Buscall, R., “Whither compressional rheology?” J. NonNewtonian Fluid Mech. 157, 151–157 (2009).
 Tomas (2007) Tomas, J., “Adhesion of ultrafine particles—a micromechanical approach,” Chem. Eng. Sci. 62, 1997–2010 (2007).
 Witten and Pincus (2004) Witten, T. A. and Pincus, P., Structured Fluids – Polymers, Colloids, Surfactants (Oxford University Press, 2004).
 Wolf et al. (2005) Wolf, D. E., Unger, T., Kadau, D., and Brendel, L., ‘‘Compaction and flow of cohesive powders,” in Powders and Grains 2005, edited by R. GarciaRojo, H. J. Herrmann, and S. McNamara (2005).