{'fraction_non_alphanumeric': 0.02404446413468122, 'fraction_numerical': 0.007571791050787386, 'mean_word_length': 4.467958599427439, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Liquid migration in sheared unsaturated granular media

We show how liquid migrates in sheared unsaturated granular media using a grain scale model for capillary bridges. Liquid is redistributed to neighboring contacts after rupture of individual capillary bridges leading to redistribution of liquid on large scales. The liquid profile evolution coincides with a recently developed continuum description for liquid migration in shear bands. The velocity profiles which are linked to the migration of liquid as well as the density profiles of wet and dry granular media are studied.


Introduction
How does liquid spread in wet granular media under shear? We know from everyday experience that mixing large amounts of powders and liquid is a challenging task. Liquid bridges appearing in wet granular matter exert cohesive forces onto the grains resulting in substantial changes of the mechanical properties [1] highly enhancing the difficulty in obtaining a homogenous mixture of liquid and powder. The development of advanced techniques for the improvement of mixing processes requires gaining knowledge about the influence of the grain scale redistribution of liquid onto the large scale liquid migration in a granular medium.
The fundamental difference between ordinary cohesive and wet granular matter is that in the wet case the cohesive force acts only after liquid bridge formation exhibiting hysteresis, resulting in a rich variety of observed phenomena [1,2]. The interaction distance of two grains is determined by the rupture distance of a capillary bridge. The forces and rupture distances were experimentally measured in Ref. [3] where the authors provided well established formulas suitable for numerical simulations. Discrete element methods are widely used to model wet granular matter. The simplest model to account for hysteresis is the minimal capillary model [1]. Simplified models for capillary clusters were introduced in Ref. [4]. Most simulations implement geometry imposed or homogenous capillary bridge volumes [5]. Richefeu [6] implemented a redistribution scheme for liquid and investigated the mechanical properties of wet granular matter in direct shear but did not allow for larger changes in liquid concentrations. A very important problem which was rarely studied in the past is the transport of liquid in wet granular matter. It is still unclear which mechanisms are responsible for an altering liquid concentration. M. Scheel et al. [2] studied the equilibration of liquid bridges through thin films on the particles due to Laplace pressure differences. However, liquid is also redistributed due to relative motion of grains leading to steady rupture and creation of liquid bridges. In this paper, we study this redistribution process in detail and present a model for capillary bridges taking into account the redistribution of liquid after bridge rupture. As opposed to previous studies where the volume of the bridges is either homogeneous or given by the local geometry we explicitly account for the volume of individual capillary bridges given by the amount of liquid which was trapped by the roughness of the particles. The redistribution of liquid gives rise to locally changing liquid concentrations which will be studied in periodic simple shear between two rough walls [7,8].

Model
We model wet granular matter at low liquid content under shear. For this purpose we use Contact Dynamics to model spherical hard particles and develop a model for capillary bridges with adhesive forces.

Particle dynamics
Particles are modeled using Contact Dynamics originally developed by J. J. Moreau which is suitable for rigid particles [9,10]. Here, contact forces are calculated based on perfect volume exclusion and Coulombian friction. In the framework of Contact Dynamics, more contact laws have been established such as rolling friction, hydrodynamic lubrication and cohesion [11,12,13]. The fundamental difference to standard discrete element methods lies in the determination of the contact forces which are calculated based on perfect volume exclusion of the contacting particles. For instance, the contact normal force F n between two grains is calculated such that the two grains do not overlap at the next time step. A tangential constraint force F test t is subsequently calculated under the condition of no slip at the contact which requires the relative tangential velocity to vanish at the next timestep. However, |F t | must not exceed the threshold value µ|F n | where µ is the friction coefficient. If the previous calculation of |F test t | is greater than the threshold we set |F t | = µ|F n | which is a sliding, energy dissipating contact, else we set F t = F test t .

Capillary bridges
A wetting liquid forms a wetting layer on a grain surface. As soon as two grains touch, liquid from the wetting layer accumulates at the contact point to form a capillary bridge. This meniscus gives rise to an adhesive capillary force. The force was measured in detail in Ref. [3] and the authors derived an empirical formula for the capillary force given by where V is the bridge volume, r the average radius of curvature of the two spheres, Γ the surface tension of the liquid air interface and s the separation [3]. This empirical formula for the forces is well established and valid for small capillary bridges up to V ≈ 0.03r 3 . The contact angle θ is set to zero in this work. We see that the capillary force is independent of the volume when s = 0. However, roughness influences the attractive force as soon as the length scale of the liquid bridge is comparable to the length scale of the roughness [14]. In that case, the force can be much smaller. The authors distinguish three regimes for very small, for intermediate and for large amounts of liquid. The first two regimes are characterized by a monotonically increasing capillary force with volume before the force becomes constant in the third regime. We simplify the force law by choosing a threshold value below which the attraction vanishes. When two grains separate, the liquid bridge eventually ruptures. The above mentioned authors also provide an empirical formula for the rupture distance s c which is for zero contact angle and small volume. There are two relevant mechanisms leading to a redistribution of liquid. First, in case of rupture, the liquid is sucked back onto the grains very fast. At low water contents, the pressure in the bridges compared to the pressure in the wetting layers is small, such that the liquid is drawn into the already existing bridges [1]. The time scale for this redistribution process is expected to be of the same order as the time scale for bridge formation between two grains when the roughness is saturated with liquid. If the relative velocities of the grains are small enough this redistribution of liquid can be considered to take place instantaneously.
Secondly, as shown experimentally in Ref. [15], there is a flux of liquid between capillary bridges through the wetting layers on the bead surfaces or through the vapor phase driven by Laplace pressure differences. However, the equilibration of capillary bridges can be neglected if the lifetime of a capillary bridge is much smaller than the equilibration time scale. For a fixed equilibration time scale, this criterion is fulfilled for sufficiently high shear rates since the relative velocities of the grains become larger and the lifetimes of the capillary bridges become smaller with increasing shear rate. Therefore, a model in which the equilibration of liquid bridges is neglected and in which the liquid is instantaneously redistributed to neighboring bridges after bridge rupture is indeed valid in a certain range of shear rates: The relative motion of the grains must be slower than the redistribution of liquid after bridge rupture but faster than the equilibration of liquid bridges. The time scale for bridge formation is less than a millisecond for grains of one millimeter diameter [1,16]. On the other hand, typical equilibration times of liquid bridges are of the order one to five minutes [15] for glass beads of half a millimeter diameter which is a much larger time scale. Since we expect the time scale for liquid redistribution after bridge rupture to be of the same order as the time scale for bridge formation, the assumption of instantaneous liquid redistribution after bridge rupture and neglection of equilibration of liquid bridges will be valid in the range 10 −2 s −1 ≪γ ≪ 10 3 s −1 which spans a large range of applicable shear ratesγ.
Based on these considerations we propose the following model for redistribution: Particles can carry an arbitrary amount of liquid V f > V min in a wetting layer where V min is a fixed amount of liquid which is trapped in the roughness [2]. When two particles come into contact and at least one particle has V f > V min , a capillary bridge is instantaneously formed. Since for each grain, V min is fixed, the available liquid for bridge formation is V f −V min therefore, without loss of generality, we can use V min = 0 and the bridge volume V is given by the sum of the two involved films, i.e. after bridge formation If the separation of particles exceeds the critical distance s c , the bridge volume V rup is equally split between the two particles P j such that V j f = V rup /2, j ∈ {1, 2}: At the same time, if there exist further contacts, namely capillary bridges or dry contacts, all liquid is equally distributed among them such that the new volume of a neighboring contact C j i , i ∈ {1 · · · N j } of particle P j is given by where N j is the number of neighboring contacts of particle j. An example of this redistribution scheme is shown in fig. 1 where initially, all bridges have the same volume. If N j = 0 the liquid is kept in the wetting layer of particle j for the formation of a bridge at a later time. To prevent formation of liquid clusters via bridge coalescence [2] we use an upper threshold V max of the bridge volumes which must not be exceeded. In that case, all contacts of that particle with V < V max are filled first and if all are filled completely, the liquid remains in the film. This case however, happens very rarely since we consider relatively low liquid contents only. In our model, liquid mass conservation is ensured, whereas in suction controlled models, e.g. [5] the Laplace pressure is prescribed such that the total amount of liquid can change.

Simulated Setup
We simulate normal stress controlled periodic simple shear between two rough walls. We use particle radii uniformly distributed between 0.8 and 1 in units of the largest particle radius R. Fig. 2 shows the geometry used in our numerical experiments. The bottom wall is fixed at position z = 0 but moving at constant speed v shear along the x direction. The top wall is fixed in x direction and a pressure P is applied onto it such that the position of the wall z wall (t) fluctuates around z wall = z wall (t) . Gravity is neglected and we use periodic boundary conditions in x and y direction. The system dimensions are L x = 20R, L y = 12R, z wall ≈ 78R. The results will be presented in terms of the global shear strain γ, which is the shear displacement v shear t divided by the system height L z . As proposed by Rognon et al. [17] there are two dimensionless parameters which govern the mechanics of the system, the inertial number I =γR ρ/P wherė γ = ∂v x /∂z is the local shear rate and the cohesion number η = F max c /P R 2 = 2πΓ/P R which is the ratio of the largest cohesive force to the force exerted on a particle by the pressure P . We first study the simplest case η = 0 where shear flow is approximately homogeneous which can be realized in experiments by choosing a sufficiently high pressure. We subsequently study the case where η > 0 and in both cases, we fix I ≈ 0.008. Initially, we impose a Gaussian distribution of liquid bridge volumes. As a function of the bridge position z the volumes are initialized according to V b (z) = A exp(−(z − z wall /2) 2 /σ 2 0 ) with amplitude A and width σ 0 and we monitor the evolution of the liq-uid distribution during shear. Here, the shearing walls are assumed to be hydrophobic.
4 Results: η = 0 In this section, we present results in the limit of vanishing cohesion number. We show the obtained velocity profiles and the consequences for the liquid distribution.

Velocity profile
Velocity profiles have already been studied for noncohesive systems in detail [7,18]. It was shown that under high confining pressures the system is alternating between diffuse shear and shear banding occurring at random height z. These shear bands are however not stable, in fact they disappear and reappear constantly at different positions. However, when averaging the velocity profile over long time scales it becomes linear apart from a slight S-shape which would vanish in the limit I → 0 [18]. The time averaged velocity profile for our system is shown in the inset on the left of fig. 3.

Liquid profile
During shear the liquid profile changes. In fig. 2 we observe that the initial Gaussian liquid distribution (left picture) spreads towards the top and bottom walls (middle and right pictures). There are two relevant processes involved causing the spreading of liquid. It is known that in plane shear flows particles undergo a diffusive motion and therefore, also liquid which is carried by the menisci will diffuse in space [19]. Secondly, there is a transport of liquid associated to liquid bridge rupture. As explained in section 2.2 the liquid is redistributed to all neighboring liquid bridges after bridge rupture which means that locally, after a bridge rupture event, there is a liquid flux away from the rupture point. Fig. 3 shows the evolution of the liquid distribution at shear strains γ = 0, 30, 60 where Q is the average bridge volume per particle. The data was averaged over a small strain interval, five independent runs and in slices of width 2R. We observe the liquid distribution to remain Gaussian (solid lines) with an increasing width σ(γ) which is exactly the case for diffusive processes. The solution of the standard diffusion equation with Gaussian initial condition is a Gaussian for which σ 2 (t) increases linearly in time [20]. Since γ is proportional to t it can be expressed as σ 2 (γ) = σ 2 0 + 4Dγ where D is the diffusion coefficient. The inset in fig. 3 shows σ 2 as a function of γ (blue crosses) and an excellent linear fit to the data  After bridge rupture, liquid was redistributed to neighboring capillary bridges taking into account the Laplace pressure in the bridges as well as the distance from the rupture point to the neighboring bridges. A modified diffusion equation was derived by considering that liquid fluxes happen via bridge rupture events which are proportional to the volumes of the ruptured bridges and to the local rupture rate. The continuum equation was given bẏ where Q b is the average bridge volume, B is the bridge rupture rate and C is a constant. We notice that the equation is applicable to the plane shear geometry as well: Since the velocity profile is linear, we expect the bridge rupture rate to be constant such that we recover the ordinary diffusion equationQ b = D∂ 2 /∂z 2 Q b . Assuming that the amount of contacts per particle is constant, Q b can be replaced by Q.

Dependence on liquid content
The influence of the liquid content on the diffusion constants was investigated by varying the amplitude A of the Gaussian liquid distribution. Here, we used a smaller system for computational reasons. Fig. 4 shows the diffusion coefficients D for different amplitudes A of the Gaussian. It can be seen that the diffusion coefficients decrease with increasing liquid content. A certain dependence is indeed expected: Since the rupture distance of a capillary bridge scales as s c ≈ V 1/3 b the strain needed to rupture a capillary bridge increases with increasing bridge volume. The strain needed to break a contact in a dry granular medium is of order unity [17] such that the rupture time for dry granular media can be estimated by T sc=0 ∼ 1/γ. The additional strain γ r needed to rupture a capillary bridge can be estimated by γ r ∼ s c /R. Thus, the rupture time T sc=0 for dry contacts is increased by T sc ∝ T sc=0 + s c /Rγ ∝ 1 + s c /R such that we expect the diffusion constants to decrease with Q as D ∼ 1/(1 + bQ 1/3 ), where b is a fit parameter. Fig. 4 shows the diffusion coefficient D for different amplitudes of the Gaussian which can be well fitted by D(A) = a/(1 + bA 1/3 ) with fit parameters a and b.

Results: η > 0
Now we study the case of finite cohesion number η. We obtain a non-homogenous velocity profile and investigate its influence on the liquid migration.

The non-homogenous case: velocity and density profiles
As opposed to the case η = 0 here the capillary forces play a dominant role since weaker tangential stresses are less capable to rupture the bonds between the grains. Since wet grains stick together, the velocity profile develops a plateau and the shear rate drops inside the wet area. Examples for such plateaus can be seen in fig. 5 for two different cohesion numbers. The shearing of a mixture of two granular materials with different friction coefficients was studied recently by Unger [22]. The upper half of the shear cell was filled with material having a larger coefficient of friction than that of the lower part. He observed that in such a system, the upper part with larger friction coefficient undergoes shear hardening. The shear rate was found to be given bẏ where z c is the middle position of the upper part of the cell and therefore, the velocity profile is given by a sinh function. Since in our case friction is effectively enhanced due to the attractive force, the same arguments as in Ref. [22] are applicable to our system. Indeed, fig. 5 shows typical velocity profiles appearing during shear which can be fitted by E + A sinh(k(z − z c )/z wall ) inside the wet region. Since the liquid profile changes, we only show velocity profiles averaged over a small strain interval of 0.4. We see that increasing the cohesion number η leads to a non linearity in the velocity profile which becomes more pronounced with increasing η. The parameters E, A, k and z c fluctuate very much over time as shown exemplarily for the parameter k in fig. 6. We observe k to significantly increase with increasing cohesion number η. As pointed out by Unger, the amount of agitation a layer in the x-y plane receives from the next upper respectively lower layer decreases with increasing k. Here, the same behavior is found, i.e. the shear rate in the wet region drops more for larger cohesion numbers. The explanation for the appearance of a velocity profile given by a sinh function in Ref. [22] was based on the presence of a shear stress heterogeneity. However, we could not find noticeable stress heterogeneities in our systems such that we believe that in our case, a larger effective viscosity of the granular medium caused by cohesion is responsible for the locally decreased shear rate. Next, we study the density profiles for different cohesion numbers. The density φ of cohesive granular media in plane shear was found to decrease with increasing η [17], due to the formation of large stable pores in the medium during shear. Somewhat counterintuitively, we find that the density shown in fig. 7 for three different cohesion numbers η increases in the wet region despite the presence of cohesive forces. On the other hand, the shear rate respectively the local inertial number in the wet region is decreased which in fact inhibits the for- mation of larger pores. Our observation coincides with Ref. [22] where density differences originate from having two materials with different friction coefficients, although the microscopic origin for increased density is different in our work.

Liquid profile: The non-homogenous case
We showed that under shear the liquid migrates due to diffusion and liquid bridge rupture. Both processes are however, a function of the shear rate. Since the shear rate drops in the wet region, although the global shear rate is the same as for the case η = 0, we expect the liquid to migrate on larger time scales which reflects a frequently encountered problem in industrial applications when mixing granular matter with liquid. If the humidity profile was non-homogenous but the grains were wet everywhere, we would recover the linear velocity profile again which enhances the mixing properties. We can already deduce that mixing dry and wet powder is even more difficult than equalizing the liquid distribution of completely wet samples because of cohesion and the appearance of powder clumps. The evolution of the liquid profile for η = 0.47 is shown in fig. 8. As opposed to the case η = 0 we see a deviation from a Gaussian at the border between the wet and the dry region at z ∼ 50. It corresponds to the place where the shear rate profile starts to rapidly deviate from a constant. At this point, we cannot describe anymore liquid migration as a purely diffusive picture but it is necessary to take the shear rate profile into account. We apply the continuum descriptioṅ from Ref.
[21] to our system. To solve eq. (5) numerically we use the measured shear rate profile as an input to eq. (5) and the constant C was set to 0.475 in units of R 2 . The solid line in fig. 8 shows the numerical solution of eq. (5) which corrects the deviations from a Gaussian at the interface between dry and wet granulate. To have a better view on the accuracy of eq. (5) we plot ln(Q) as a function of the height in the inset on the right (filled symbols) in fig. 8 as well as the results for η = 0 for comparison. We see that the deviations from a Gaussian are much more pronounced for η = 0.47 than for η = 0 caused by the shear rate. This shows that in our model, the shear rate profile determines how liquid migrates within the system and must be taken into account in a diffusive description for liquid migration. The reason why we see deviations from a Gaussian for η = 0 is on the one hand due to the boundaries and on the other hand due to the slight non-constancy of the shear rate profile. The observation of purely diffusive liquid migration in a numerical experiment would be possible in larger systems and by choosing appropriate boundary conditions without walls, e.g. Lees-Edwards boundary conditions [23]. The inset in fig. 8 on the left shows the evolution of Q in slices corresponding to the data points in the main panel as a function of strain (red solid lines) which is followed by the numerical solution (black dashed lines) of eq. (5) very nicely. Since the local shear rate drives the liquid migration, improving mixing properties of grains and liquids can be achieved by decreasing the cohesion number. Experimentally, this requires increasing the pressure onto the top wall.
Since the liquid migrates during shear, the velocity profiles are expected to change in time. The shear rate can be fitted byγ(z) = c+a(z −z c ) 2 with fit parameters a, c and z c for sufficiently small |z − z c |. In our simulations, we found that c does not change much with strain for fixed η. Therefore, the parameter a is a measure for the width of the plateau which increases with decreasing a. Fig. 9 shows the parameter a as a function of γ for two different cohesion numbers. In both cases, a decreases consistently with strain. The fact that the velocity profile is affected by the liquid distribution could be used in experiments to indirectly measure the position of the interface between the dry and wet region by measuring the velocity profile of the granulate.  The parameter a of the fit of the shear rate is plotted as a function of the strain γ. Since liquid is migrating, the plateau in the velocity profile widens which is reflected by the decreasing fit parameter a with γ.

Conclusion
We studied liquid migration in unsaturated granular media in plane shear. We found that liquid migrates diffusively if the shear rate profile is constant along the sample in agreement with experimental observations made in process engineering [24,25]. However, a nonhomogeneous shear rate profile leads to a different liquid migration pattern in agreement with Ref. [21]. The influence of the cohesion number η on velocity profiles was studied and linked to the evolution of the liquid profiles. Homogenous mixtures of liquid and powder are achieved by increasing shear rates and stresses. Direct experimental measurements of the liquid distribution in sheared granulates are difficult, but would be feasible in indirect measurements of the velocity profiles since the liquid distribution affects the velocity profile.