DEM analysis of some effects due to capillary forces in sands

. In this paper, the Discrete Element Method is used to study the e ﬀ ect of capillary forces, in equilibrium with the ambient relative humidity, on the response of irregular arrangements of equal-sized spheres, simulating a ﬁne uniformly graded quartz sand. The e ﬀ ect on the isotropic compression was investigated by applying an increase in total stress under constant relative humidity (wet and dry) and drying under di ﬀ erent constant total stress values. The e ﬀ ect of the capillary forces on the shear strength was evidenced by the simulation of the instability of a cylindrical sand column due to drying.


Introduction
The Discrete Element Method (DEM) is a very convenient computational tool to study the micromechanical (particle level) effects of capillary forces on soils. The capillary inter-particles forces due to water menisci that form at the contact points, including their dependence on relative humidity (RH), can be applied to the modelled soil particles and the effect of their variation on soil response analysed.
In this paper, the isotropic compression behaviour of sand in different wetting conditions and the crumbling of a drying initially wet sand column is modelled with the DEM.

Insertion of capillary forces in the DEM
The DEM simulations were made with the PFC3D program that models the movement and interaction of spherical particles, as described by Cundall and Strack [1]. The nonlinear Hertz-Mindlin contact model [2] was used to allow for the highly variable contact stiffness. The magnitude of the normal contact force, F, and the maximum elastic normal contact stress, σ max , are given by where E 1 and E 2 are the elastic moduli, ν 1 and ν 2 the Poisson's ratios, R 1 and R 2 the radius of each sphere and δ the depth of indentation. σ max is given by the classical solution for nonadhesive elastic contact. The water meniscus connecting two particles is under tensile stress and compress both particles together around the initial contact point. As shown in Figure 1, the appli- * e-mail: cpereira@lnec.pt ⇐⇒ Fcap Figure 1. Representation of the equivalence between the capillary forces applied at each contact between particles (left) and the resultant force of all the capillary forces of each particle applied at the centre of the particle (right). cation of the capillary forces at the contact points is equivalent to the application of the resultant force at the centre of each particle. Figure 2 presents the evolution of the capillary force, F cap , with the relative humidity, RH, obtained from the precise numerical solution of the combined Laplace-Young and Kelvin equations for capillarity. The capillary force is created by a water meniscus, in equilibrium with a specific RH value, between two rough spheres in contact, made of the same material, of equal radius R, separation distance d and wetting contact angle θ on the surface of the sphere. It is assumed that the nanoscale asperities of the particles' surface impose a separation between them and as such there cannot be a perfect contact even for highly polished surfaces. There can be significant capillary forces even for RH = 1 (same pressure on both sides of the surface of tension). In this case, the capillary force is produced by the surface of tension, which assumed a shape of a catenoid (surface with a null mean curvature at all points). Moreover, starting from RH = 1, the decrease of the relative humidity originates the increase of the capillary force until a maximum value and, after that, the cap- Fcap F a cap Fcap(RH = 1) Figure 2. Evolution of the capillary forces, given by a precise numerical solution of the combined Laplace-Young and Kelvin equations for capillarity, F cap , and by an approximate solution proposed in [3], F a cap , with the relative humidity, RH, between to equal radius spheres, R = 0.1 mm, with the separation d = 0.5, 1, 3 and 10 nm between them, the contact angle θ = 0 • and T = 20 ℃.
illary force decreases to zero. The maximum value of the capillary force occurs at values of RH close to 1.
Some approximate solutions can be found in the literature based on circular or parabolic approximations of the shape of the meniscus' surface [3,4]. These approximate solutions do not strictly enforce a constant mean curvature at all points of the meniscus' water-air interface. Butt and Kappl [3] suggested the following approximate solution for the capillary force, F a cap , where γ LG is the surface tension and λ k the Kelvin's length. For T = 20 ℃, γ LG = 72.43 × 10 −3 N/m and λ k = 0.539 nm. This approximate formulation gives good predictions for the evolution of the capillary force at relative humidities lower than the one corresponding to the peak force, as represented in Figure 2. Therefore, Eq. (3) should only be used for values of relative humidity lower than the value correspondent to the maximum capillary force. Taking into account that the maximum relative humidity used in the following calculations is 0.9, the approximate solution of the capillary force is used in all the following calculations. In this paper, the evolution of the capillary force is always a function of the relative humidity and the distance d (with the minimum value equal to 1 nm, defined by the size of the surface asperities), contrary to other works in the literature [5][6][7], where the water volume of the meniscus is used, and that, in general, are not in thermodynamic equilibrium. The effective stress tensor acting on an arrangement of spheres in static equilibrium can be computed by a homogenisation procedure as described in Cundall and Strack [8]. The effective stress tensor, σ, for a volume V with n p particles is given by where V i p and σ i p are respectively the volume and average stress tensor of particle i. The average stress tensor at a particle is obtained as where F j c is the force vector acting at the contact point j with position vector x j c , F p the resulting force vector due to all the water menisci on the particle, with position vector x p and n c the number of contacts on the particle. The applied capillary forces are those due to the water menisci, also contributing to F j c . Finally, the capillary stress tensor, σ cap , is given by where σ tot is the total stress tensor (due to the external forces applied by the boundary walls on the particles' arrangement).

Isotropic compression simulations
In this section, the influence of capillary forces in the isotropic behaviour of unsaturated sands is evaluated using the DEM with the characteristics described in section 2.
The properties of quartz were attributed to the spherical particles: shear elastic modulus G = 36.4 GPa, friction coefficient µ = 0.2 [9] and Poisson's ratio ν = 0.167. The spherical particles have R = 0.1 mm and intend to represent a uniformily graded fine sand. The separation distance due to roughness of the particles is assumed to be d = 1 nm, and the contact angle θ = 0 • (the lowest possible value). Low values of θ are associated with drying evolving quartz [10]. Other values, like 10 • or 20 • , do not change significantly the results obtained. The isotropic behaviour was analysed, considering both dry and wet conditions. The application of the mean total stress, p tot , at the boundaries was imposed by a servocontrol algorithm, which controls independently the velocity of each boundary wall in the three orthogonal directions (x, y, z). The boundary walls have no friction, and capillary bridges are not formed with the particles.
To avoid very large calculation times, the samples have approximately 600 particles, randomly generated inside an initial cubic volume (limited by the servocontrolled walls) of edge length 2 mm. No overlapping of particles in the initial assembly was allowed. The servocontrol algorithm stops when the arrangement of spheres can sustain the desired total normal stress imposed on the wall under a prescribed tolerance. There might be significant stress relaxation when the walls are stopped. The value of the relaxation parameter was chosen such that this relaxation was insignificant.    Figure 3 represents the results of the numerical simulation of three distinct mean total stress paths: isotropic loading of a dry sample (red lines), isotropic loading of a wet sample, with RH = 0.9 (blue lines), and the drying of the wet sample (RH = 0.9 → 0.34) (0.34 is the highest value of RH corresponding to a completely dry soil, when F cap = 0), at constant p tot = 0.1 kPa (green lines), 1 kPa (brown lines) and 10 kPa (gray lines). p tot is the mean total stress (applied by the sample boundary walls).
The initial void ratio of the dry sample is 0.828. It was not possible to obtain any arrangement with larger void ratios with enough stiffness to sustain the imposed p tot . This value is also quite similar to the maximum void ratios of uniformily graded dry or fully saturated sands, generally between 0.8 and 0.85. So, even for very small mean total stresses, it is still almost the same void ratio (Figure 3a). The increase of p tot reveals a very stiff arrangement, even for low values of the average coordination number C. The coordination number of a given particle C is the number of contacts. Even for p tot = 1 MPa, C is still lower than the correspondent value of the simple cubic arrangement (Figure 3b), which is C = 6. C increases with p tot .  Figure 4. Evolution of the stiffness K of the dry sample with p, considering the shear modulus of the quartz G quartz (red line) and 0.1G quartz (blue line). Also, the stiffness of two real sand samples, the Ottawa sand (D 50 = 0.2 mm) and the Sacramento River sand, are presented [11]. The mean values of the maximum stress on the spheres' surface at the contact between spheres, σ max , is presented for both numerical simulations. bulk modulus, K, was calculated by where ∆p is the variation of the mean effective stress (Eq. (4)). The stiffness of the spheres' sample with the elastic constants of quartz is much higher than the values for real sands, as represented in Figure 4, which includes the values of Ottawa (D 50 = 0.2 mm) and Sacramento River quartz sands, obtained from Yamamuro and Lade [11]. Assuming perfect spheres, the average maximum normal contact stress σ max is in the order of [GPa]. Possible contributing factors for the lower stiffness of real quartz sands can be that the actual shape of the particles has much higher curvatures at the points of contact, which would induce plastic yielding at the contacts. Another factor can be the yielding of the asperities at the contact. The variability of the elastic constants of quartz could also be a contributing factor. The result with the shear modulus equal to 0.1G quartz is also presented. This reduction in elastic constants produced bulk modulus' values closer to the measure ones. The initial wet sample was obtained with the same compaction procedure used for the dry sample. Additionally, it was assumed, during the servocontrol procedure, that a water meniscus in equilibrium with RH = 0.9 is formed whenever two particles are separated by a distance small enough for it to occur (the formation of water menisci can be very fast [12]). In order to represent the action of the water meniscus, a pair of equilibrated compression forces is added to the centre of each particle.
The wet sample has, at p tot = 0.1 kPa, a void ratio equal to 1.226. This value is substantially higher than the normal values of e max in sands. In practice, to achieve this high void ratio, it is necessary to remould by carefully mechanical manipulation the wet (not saturated) sand. This loose arrangement of spheres is extremely compressible (blue line of Figure 3a). The average coordination number increases with p tot to values higher than the values obtained with the dry sample (Figure 3b).
From the wet sample, three drying simulations were conducted at constant p tot . The decrease in RH implies a decrease of the effective stress p, due to the reduction of the capillary stress (Figure 3d), and the void ratio e (Figure 3c). This is not a stable configuration in the sense of Drucker's stability postulate becauseṗε p v < 0, under an effective stress decrease, while it is stable under an increase in the effective stress p (Figure 3a).ε p v is the rate of the volumetric plastic strain. This configures a metastable structure, only possible due to the capillary forces acting at the contact between particles. Without capillary forces, particles' arrangements with void ratios higher than the ones of dry arrangements are not achievable. Taking into account that the effective stress also decreases due to full saturation (wetting), the same type of instability will occur. During the drying process, the void ratio converges to the values of the dry sample. The difference between the void ratios at points A d , B d and C d and those of the dry sample (Figures 3a and b) are due to random variations in the arrangements. The maximum capillary stress is approximately 0.8 kPa and corresponds to C = 5.2. The capillary stress at point C w is higher than at points B w and A w because it has more menisci formed. The number of contacts decreases suddenly when the capillary force disappears (Figures 3b). Figure 5 represents the contacts' network of configurations A w and C w . The normal force between two spheres in contact is represented by a cylinder connecting both centres, with the radius being proportional to the magnitude of the force. The colour indicates the number of contacts C of each sphere. The magnitude of the shear force is represented by a thin disk positioned at the contact point with the excess radius relative to the cylinder representing the magnitude of the shear force. When the shear force is higher than 90% of the maximum shear force given by friction, the disk is coloured black. When the relative importance of the capillary stress decreases from A w to C w , there is a higher variability between different load carrying paths. The values of contact forces are more uniform in A w because the capillary forces are practically equal in all the contact points, and are much higher than the contribution of the mean total stress to the contact forces.

Crumbling of a drying wet sand column
A different aspect of the action of capillary forces can be illustrated by their influence in the increase of the shear strength, such that it is able to sustain a vertical sand column. The capillary forces can be removed by drying in which case the column crumble into a pile of dry sand with slopes inclined at the angle of repose, as illustrated in Figure 6, where a cylinder of compacted wet beach sand was dried under constant ambient temperature. The crumbling starts at the top edge of the column as the relative humidity decreases at a much faster rate closer to the boundary. The irregular shape is probably due to fluctuations in the void ratio and some residual salt. The last remaining lumps are extremely unstable. In this section, the results of the modelling of a suitable scale sand column are presented. For instance, to model the sand column with the same size as that represented in Figure 6, approximately 170 million particles would be necessary. As such, a rather smaller cylindrical column with a diameter of 2 mm, 2273 spherical particles, R = 0.1 mm and the properties of quartz, is modelled.
The assembly of the column model starts with the filling of a cylindrical container with randomly distributed spheres. Then, the gravity was turned on, causing the particles to settle at the bottom of the cylindrical container. The next step was the introduction of the capillary forces due to RH = 0.9. Finally, the container was removed, leaving the column standing by itself. In all these steps, the model was allowed to reach equilibrium. The final height was 5 mm and e = 0.67. The contact between particles and the container is assumed frictionless, and no capillary forces were connecting the particles and the wall.
Two different approaches for simulating the drying were considered in the simulations: one (A) was to assume that the RH decreases rapidly and uniformly in the column and the other (B) assumes, more realistically, that RH decreases much faster at the boundaries than inside the sample. Approach A is much faster and allows to proceed to the final equilibrium configuration of the pile of parti-cles. Approach B is necessary to simulate the experimentally observed crumbling, starting at the top edge of the column. Without the contribution of the capillary forces, the cylindrical column turns into a conical pile of particles The slope of the final conical shape of the pile of sand has an angle of approximately 23 • . This value is lower than the angle of the pile of real sand, shown in the Figure, that has an angle of approximately 30 • . The lower value can be due to different shapes or the nature of surface roughness of the particles.

Conclusions
In the void ratios commonly found in dry uniformly graded sands (less than approximately 0.8) the increase of the effective stress due to wetting is minimal to cause any significant decrease of the void ratio. This is due to the high bulk modulus of the particle arrangement. For an unusual very loose configuration, that is only possible to obtain by remoulding wet sand at very low stresses, a considerable decrease of the void ratio occurs when the effective stress is reduced by full drying or saturation, but only at very low isotropic total stresses not of engineering interest. Also, it was shown how to model capillary effects in isotropic compression of sphere arrangements with the DEM.
The approach used to introduce the capillary forces in the DEM, enable the simulation of the effect that the change of shear strength has on the stability of a wet sand column.