Dense bidisperse suspensions under non-homogeneous shear

We study the rheological behaviour of bidisperse suspensions in three dimensions under a non-uniform shear flow, made by the superimposition of a linear shear and a sinusoidal disturbance. Our results show that (i) only a streamwise disturbance in the shear-plane alters the suspension dynamics by substantially reducing the relative viscosity, (ii) with the amplitude of the disturbance determining a threshold value for the effect to kick-in and its wavenumber controlling the amount of reduction and which of the two phases is affected. We show that, (iii) the rheological changes are caused by the effective separation of the two phases, with the large or small particles layering in separate regions. We provide a physical explanation of the phase separation process and of the conditions necessary to trigger it. We test the results in the whole flow curve, and we show that the mechanism remains substantially unaltered, with the only difference being the nature of the interactions between particles modified by the phase separation.


I. INTRODUCTION
Dense suspensions of hard particles immersed in a Newtonian solvent are common in many natural and industrial applications, such as the pharmaceutical industry (e.g. the process of blending powders for tablets production), food industry (e.g.powdery products), powder metallurgy (e.g. the processes of compacting and sintering blended pulverized metals), and sediment transport (e.g. transportation of nutrients and landscape shaping) [1][2][3][4][5], with applications reaching biofluidics and bacterial suspensions [6].These materials, when subject to a shear-rate γ, often show peculiar rheological behaviours that can increase or decrease the fluidity of the suspension.Among these behaviours, shear-thickening, i.e. a reduction of the fluidity properties of the suspension when the latter is subject to an increasing γ, is of particular interest.Shear thickening is perhaps the most astonishing and most studied non-Newtonian behaviour of dense suspensions, and until few years ago it was far from being understood.The reduction of the fluidity is now attributed to the modification of the internal microstructure of the fluid [7]; a recent theory experimentally proved that shear-thickening is triggered by the close interactions between the particles that appear in the form of frictional contacts (due to the microscopic asperities on the surface of the particles) and constrain their relative movements, causing an enhancement of the viscosity of the suspensions [8][9][10][11][12][13][14][15][16][17][18].The intensity of the shear-thickening can be quantified as an increase of the effective viscosity of the suspension η r and becomes stronger as the volume fraction of the suspension ϕ increases, until it displays an abrupt enhancement, behaviour known as discontinuous shear thickening [19,20].The dependence on γ and ϕ of the rheological properties generally characterizes the mechanics of the suspensions within the shear-thickening regime.When the particles are suspended in a simple homogeneous shear flow, Boyer et al. [21] demonstrated that such dual dependence can be reduced to a single parameter in analogy to granular flows [22].In fact, specifying the macroscopic friction coefficient µ = σ 12 /Π uniquely sets the volume fraction and the shear-rate, and general constitutive laws η r = η r (µ) and ϕ = ϕ(µ) can be formulated [21,23].In the previous relation, σ 12 is the shear component of the stress tensor σ and Π is the pressure, related to the trace of σ.This powerful general outcome, however, fails when dealing with spatial or temporal inhomogeneity within the flow.In fact, phenomena such as subyielding and overcompaction [24][25][26], particles migration [27][28][29] or flow instabilities that can lead to banding [30][31][32][33] and segregative phenomena [34] cannot be captured by such constitutive-law, since it would require a two-phase description [23,35].The route towards a complete model able to describe real suspensions that involve the aforementioned phenomena requires more studies that shed some light onto the elusive physical mechanisms activated by the inhomogeneities.
In this direction, we propose a three-dimensional numerical study that aims at analysing non-trivial phenomena generated when bidisperse particles are suspended in an inhomogeneous Newtonian flow characterized by the linear super-imposition of a plain shear flow and a sinusoidal velocity profile, i.e., with zero mean but a non-uniform shearrate.Similar numerical experiments were carried out on two-dimensional granular flows: e.g., Saitoh and Tighe [36] studied a Kolmogorov-like flow of soft athermal disks close to the jamming transition and found that stress profiles can be reproduced by nonlocal constitutive relations that account for fourth-order gradients [see also 37].
The importance of studying non-uniform profiles is suggested e.g., by the experimental shear-rheometer, where the condition of plain shear is hard to reach due the presence of the solid walls that inevitably introduce local inhomogeneity.Also, the present configuration is a quite general configuration of any pressure driven flows in channels, which are characterised by a non-uniform mean shear-rate [38].In particular, the goal of the work is to trigger a mechanism of particles migration that results in a demixing of the two phases.In this sense, many works dealt with arXiv:2206.07916v2[cond-mat.soft]8 Oct 2023

II. METHODOLOGY
We carried out three-dimensional numerical simulations of a dense suspension of rigid, spherical particles, not subject to Brownian motions, immersed in a Newtonian fluid with viscosity η 0 .The suspension is driven by a sum of a plain shear-flow, u 1 = γx 2 , characterised by the shear-rate γ (all the quantities • within the manuscript refer to their mean value along the shearwise direction y, while the quantities •(y) to their local value) and a sinusoidal disturbance in the velocity of the form where A is the amplitude of the disturbance, a 0 is the reference length-scale of the problem (typically the radius of the smallest particle of the suspension), κ 0 = 2π/L is the fundamental wavenumber of the signal (being L the size of the computational domain), n an integer that sets the wavenumber of the signal and x j the direction of the wave.In particular, the subscripts i and j span the spatial three-directions in a Cartesian frame of reference, where x 1 , x 2 and x 3 (sometimes also referred to as x, y and z) are adopted to identify the streamwise, shearwise and spanwise directions, and u 1 , u 2 and u 3 to identify the corresponding velocity components (u, v and w).In Equation (1), the constraint i ̸ = j is imposed.
The suspensions considered are composed of N = 2 16 particles and are binary, i.e. the particles have two different sizes, with radii a 1 = a 0 and a 2 = 3a 0 .The two phases are dispersed with relative volume V 2 /V 1 = 0.25.The computational box containing the particles is a cube with size designed to contain the desired volume fraction ϕ = 0.50.The uniformity of the shear-flow at the edges of the domain is preserved through the adoption of the 3D-periodic Lees-Edwards boundary conditions [58].
The numerical investigation is performed using a validated and publicly available software, CFF-Ball-0x [59].The The numerical investigation is performed using a validated and publicly available software, CFF-Ball-0x [60].The software tackles the translational and rotational dynamics of the particles by solving the Newton-Euler equations, where the subscript i indicates the particle i 2 [1, N], being N the number of particles.The right-hand side of System (2) lists the forces and torques that are applied to the centre of mass of the ith particle, with mass m i and inertia tensor I i and cause a variation in the translational and angular velocities of the particles, here denoted by the symbols u i and !i , respectively.Such forces and torques result from the particle-fluid and particle-particle interactions of the ith particle with its jth neighbour.The nature of the interactions is indicated by the superscripts H, C and E, i.e. hydrodynamics, inelastic contacts and electrochemical potentials, respectively.
Here, we give a brief overview of the models used to support the three aforementioned contributions.For hydrodynamics, we utilized the approach developed by Ball and Melrose [61] and Mari et al. [11], which considers the e↵ects of a dense suspension of rigid particles in a low-Reynolds-number flow.These particles experience a Stokes drag and a pair-wise lubrication force [11] that arises from the relative motion of nearby particles, leading to the squeezing of fluid between them.To account for these forces and torques, we employed a linear relationship between velocities (angular velocities) and forces (torques) [11,61].This relationship is expressed as F H = R(u U 1 ), where R is a resistance matrix that takes into account only the dominant near-field divergent elements resulting from the squeeze, shear, and pump modes, as described in the work by Mari et al. [11].The far-field e↵ects are neglected.
The stick-and-slide method [62] is used to model the contact contribution, which simulates inelastic contacts with spring-dashpot systems that become active when two spheres overlap.The force is expressed as F C = (k n + n ˙ )n + k t ⇠t, with the spring-dashpot systems oriented along the normal (centre-to-centre) n and tangential t directions, FIG.1: Top row, panels a-c: sketch of the sinusoidal disturbance analysed in the respective columns.Bottom row, panels d-f : time-history in terms of strain γ of the relative viscosity η r .The colours distinguish the different cases, with the black solid lines referring to the values of η r for the particles suspended in a plain shear-flow.
software tackles the translational and rotational dynamics of the particles by solving the Newton-Euler equations, where the subscript i indicates the particle i ∈ [1, N ], being N the number of particles.The right-hand side of System (2) lists the forces and torques that are applied to the centre of mass of the ith particle, with mass m i and inertia tensor I i and cause a variation in the translational and angular velocities of the particles, here denoted by the symbols u i and ω i , respectively.Such forces and torques result from the particle-fluid and particle-particle interactions of the ith particle with its jth neighbour.The nature of the interactions is indicated by the superscripts H, C and E, i.e. hydrodynamics, inelastic contacts and electrochemical potentials, respectively.
Here, we give a brief overview of the models used to support the three aforementioned contributions.For hydrodynamics, we utilized the approach developed by Ball and Melrose [60] and Mari et al. [11], which considers the effects of a dense suspension of rigid particles in a low-Reynolds-number flow.These particles experience a Stokes drag and a pair-wise lubrication force [11] that arises from the relative motion of nearby particles, leading to the squeezing of fluid between them.To account for these forces and torques, we employed a linear relationship between velocities (angular velocities) and forces (torques) [11,60].This relationship is expressed as where R is a resistance matrix that takes into account only the dominant near-field divergent elements resulting from the squeeze, shear, and pump modes, as described in the work by Mari et al. [11].The far-field effects are neglected.
The stick-and-slide method [61] is used to model the contact contribution, which simulates inelastic contacts with spring-dashpot systems that become active when two spheres overlap.The force is expressed as F C = (k n δ + γ n δ)n + k t ξt, with the spring-dashpot systems oriented along the normal (centre-to-centre) n and tangential t directions, characterized by the spring and dashpot constants (k n ,γ n ) and (k t ,γ t = 0), respectively.The values δ and ξ indicate the normal and tangential displacements of the springs, and δ represents the normal projection of the overlapping velocity.The Coulomb's law, |F C t | ≤ µ|F C n |, applies to the tangential contribution, where µ is the frictional coefficient, set to µ = 0.5 for all the scenarios examined in this work, as in our previous study [18] where we compared our results with experimental data.The last contribution considered in this work is the electrochemical properties of the suspension, modelled as a sum of an inter-particle, distance-decaying repulsion force and an attraction force in van der Waals form [16,62].The resulting force is F E = F R + F A , where the where the superscripts E, R, A stand for electrochemical, repulsive, and attractive, respectively.The repulsive contribution is given by F R = −F 0 e −d/Ls n, where F 0 is the magnitude of the force, L s is the screening length, and d is the particle-particle surface distance, while the attractive force can be written as F A = H A a/12(d 2 + ε 2 )n, where H A is the Hamaker constant, a is the harmonic mean radius of the two particles involved, and ε is a smoothing term to eliminate the singularity when the two particles touch, i.e. d = 0.
The physics of the suspension described by the system of equations ( 2) is governed by several parameters.By analyzing the translational equation of (2) and combining the dimensional quantities, we can use Buckingham's Π theorem to identify a set of non-dimensional groups that can be adjusted to control the desired dynamics of the suspension.When using the hydrodynamic scales as the fundamental independent quantities, we obtain three nondimensional groups.The first is the Stokes number, which compares the time-scale of the particles with that of the hydrodynamics and arises from the inertial term.In the equation above, ρ p is the density of the particles (equal to the density of the carrier fluid in the cases we considered), a 0 is the typical length scale of the system (the radius of the particles of the smallest phase of the suspension), η 0 is the viscosity of the carrier fluid and γ is the shear-rate applied to the suspension.The constraint St ≪ 1 is applied to enforce the inertialess regime.The second number is the non-dimensional The second is the non-dimensional stiffness which evaluates the importance of the contacts contributions relative to the hydrodynamic term.The constraint k ≫ 1 is imposed to ensure that the particles behave as rigid [11].The normal spring constant k n is considered as the dominant term of the spring-dashpot systems, and additional constraints are applied to enforce this assumption.
In particular, we impose k t = 2/7k n and γ n γ/k n ≪ 1 ∼ O(10 −7 ), where the latter is the relaxation time of the spring-dashpot system.The third non-dimensional group is the equivalent shear-rate which measures the time-scale introduced by the electrochemical contribution and can be modified to impose a shearrate dependent rheology on the suspension.A table that recaps the parameters adopted in this work is reported in the Supplementary Material.

III. RESULTS
We start by considering a perturbation described by Equation ( 1) with A = 20 and n = 1.Thus, the particles are driven by a non-uniform shear-rate, due to the sinusoidal disturbance, but under the same mean total shear-rate.The effects of the perturbed shear-rate on the rheology of the suspensions are reported in Figure 1.In particular, the figure shows the relative viscosity (i.e., the shear-stress ratio) η r = σ 12 /η 0 γ, where σ 12 is the shear component of the suspension stress tensor σ, as a function of the non-dimensional time, i.e. the strain γ = γt, for all the possible configurations compatible with the boundary conditions (i.e., the waves u 2 , u 3 ∼ sin (x 1 ) are discarded since they are not consistent with the Lees-Edwards boundary conditions).In particular, the figure is split into two rows, where the top one sketches the velocity disturbance and the bottom one shows the time-history of η r , keeping the same colour map for clarity reasons.In the same figure, the black solid lines refers to the reference case, i.e. the particles immersed in a plain shear-flow.Note that, all the cases start from the same initial configuration, that corresponds to a bidisperse suspension with total volume fraction ϕ = 0.50, and dispersion ratio λ = a 2 /a 1 = 3.The suspensions are sheared with a high mean shear-rate, i.e., γ/ γ0 = 1 (a description of the latter quantity can be found in the Supplementary Material), and such configuration can be located within the shear-thickening region of the flow curve [11,59].From Figure 1, we observe that three out of four cases do not show significant alteration in terms of η r , apart from a small increase of the relative viscosity; on the other hand, the remaining case manifests substantial modifications.This case corresponds to the sinusoidal disturbance acting in the shear plane x − y, i.e. u = 20 γa 0 sin (κ 0 y) (red solid line) and leads to a reduction of the effective viscosity by approximately 50%.
In the following part of this section, we will investigate the cause of the viscosity drop for the latter case mentioned.The investigation will be carried out via a parametric study aiming at unravelling the effect of the amplitude and of the wavenumber of the disturbance.shown in the legend within the panels).The panels in the top row, i.e. panel a and (b), show the trends of ⌘ r varying the amplitude A of the disturbance with n = 1.The panels in the bottom row, i.e. panel c and d, show the trends of ⌘ r as a function of the wavenumber n, with the amplitude of the sinusoidal wave fixed at A = 20.The colorbars on the left of panel a and c indicate the colour-scheme adopted for the amplitudes (top) and wavenumbers (bottom) analysed; the black solid lines refer to the reference case of the suspension subject to a plain shear-flow.
In the following part of this section, we will investigate the cause of the viscosity drop for the latter case mentioned.The investigation will be carried out via a parametric study aiming at unravelling the e↵ect of the amplitude and of the wavenumber of the disturbance.
Figure 2 reports the trends of the relative viscosity obtained varying the amplitude A (top row) and the wavenumber  = n 0 (bottom row) of the disturbance.At first, we analyse the e↵ect of the amplitude A; in particular, in the top row of Figure 2, we show the time-history of the relative viscosity ⌘ r (left panel), and we report the mean values of ⌘ r collected at statistically steady-state, together with their contributions breakdown (right panel).Starting from the zero-amplitude case, i.e. uniform shear flow (black line), we observe that the e↵ective viscosity monotonically decreases by increasing the amplitude, eventually saturating for A > 20.The reduction rate is strongly non-linear, with very little changes of ⌘ r for A  10, followed by a sudden decrease.This is clearly visible from the mean values of ⌘ r at statistically steady-state (right panel, total value in the legend).Therefore, as a first result, the present behaviour suggests the existence of a threshold value which activates the mechanism of the reduction of the relative viscosity.Note that, the final value of the e↵ective viscosity is reached after a long transient, i.e. ⇠ 15, and that increasing the amplitude of the disturbance contributes to shorten the convergence time of the rheological response of the suspension.
Next, we report the trends of the relative viscosity obtained varying the wavenumber  = n 0 of the disturbance.Figure 2 reports the trends of the relative viscosity obtained varying the amplitude A (top row) and the wavenumber κ = nκ 0 (bottom row) of the disturbance.At first, we analyse the effect of the amplitude A; in particular, in the top row of Figure 2, we show the time-history of the relative viscosity η r (left panel), and we report the mean values of η r collected at statistically steady-state, together with their contributions breakdown (right panel).Starting from the zero-amplitude case, i.e. uniform shear flow (black line), we observe that the effective viscosity monotonically decreases by increasing the amplitude, eventually saturating for A > 20.The reduction rate is strongly non-linear, with very little changes of η r for A ≤ 10, followed by a sudden decrease.This is clearly visible from the mean values of η r at statistically steady-state (right panel, total value in the legend).Therefore, as a first result, the present behaviour suggests the existence of a threshold value which activates the mechanism of the reduction of the relative viscosity.Note that, the final value of the effective viscosity is reached after a long transient, i.e. γ ∼ 15, and that increasing the amplitude of the disturbance contributes to shorten the convergence time of the rheological response of the suspension.
Next, we report the trends of the relative viscosity obtained varying the wavenumber κ = nκ 0 of the disturbance.As done for the amplitude A, in the bottom row of Figure 2 we show the time-history of the relative viscosity η r (left panel), and we report the mean values of η r collected at statistically steady-state, together with their relative As done for the amplitude A, in the bottom row of Figure 2 we show the time-history of the relative viscosity ⌘ r (left panel), and we report the mean values of ⌘ r collected at statistically steady-state, together with their relative contributions (right panel).Similarly to the previous cases, starting from the uniform shear flow (black line), we observe that the e↵ective viscosity non-linearly decreases when a disturbance with higher harmonics is applied to the suspension, reaching a minimum value between n = 4 and n = 8.For higher harmonics (e.g., n = 16), the relative viscosity starts increasing again, thus suggesting the presence of an optimum harmonic that minimizes the value of the relative viscosity of the suspension.This tendency is also confirmed by the mean values of the relative viscosity at statistically steady-state.Another outcome suggested by the time-history of the relative viscosity is that the harmonics with higher wavenumbers significantly reduce the convergence time (as done by the higher amplitudes), nevertheless increasing the overall level of fluctuations, and showing the presence of an optimum value (i.e., faster convergence rate) for wavenumbers included between n = 4 and n = 8.
To understand the cause of the viscosity reduction, we give a closer look to the several stress contributions to the contributions (right panel).Similarly to the previous cases, starting from the uniform shear flow (black line), we observe that the effective viscosity non-linearly decreases when a disturbance with higher harmonics is applied to the suspension, reaching a minimum value between n = 4 and n = 8.For higher harmonics (e.g., n = 16), the relative viscosity starts increasing again, thus suggesting the presence of an optimum harmonic that minimizes the value of the relative viscosity of the suspension.This tendency is also confirmed by the mean values of the relative viscosity at statistically steady-state.Another outcome suggested by the time-history of the relative viscosity is that the harmonics with higher wavenumbers significantly reduce the convergence time (as done by the higher amplitudes), nevertheless increasing the overall level of fluctuations, and showing the presence of an optimum value (i.e., faster convergence rate) for wavenumbers included between n = 4 and n = 8.
To understand the cause of the viscosity reduction, we give a closer look to the several stress contributions to the effective viscosity.In panel b and d of Figure 2, we compare the reference case (plain shear) to the suspensions driven by the non-uniform shear flow varying the amplitude A of the sinusoidal disturbance (panel b, with n = 1) and its wavenumber n (panel d, with A = 20).The three components of the stress η r are shown: the hydrodynamic contribution which includes the Stokes drag, the lubrication and the Newtonian component of the carrier fluid; the interactions caused by the electrochemical potentials and the stresses due to the particle-particle contacts.The the interactions caused by the electrochemical potentials and the stresses due to the particle-particle contacts.The suspensions analysed in the parametric study are located within the shear-thickening region of the flow curve, i.e., ˙ / ˙ 0 = 1, thus, the contributions to the shear-stress related to the contacts and the hydrodynamics are dominant (with the former being more important), while the electrochemical potentials play no e↵ective role.Within the parameters chosen in this work, we observe that the reduction of ⌘ r (for the cases with A 20 and n 1) is mostly due to a strong decrease of the contacts contribution to the shear-stress, while the hydrodynamics remains almost unchanged.
Worth noticing is the case with A = 20 and n = 16, where also the hydrodynamics contribution is significantly reduced compared to the other cases analysed.
To highlight the origin of the relative viscosity reduction, we show in Figure 3, panels a and b, the instantaneous snapshots of the suspensions studied when a statistically steady-state is reached.In particular, the two panels are organised in a 5 ⇥ 2 (or 6 ⇥ 2) matrix as follows: from left to right, the amplitude (or the wavenumber) increases, from A = 0 (n = 0) to A = 40 (n = 16); in the top row, the distribution of the larger dispersed phase is shown while, in the bottom row, the location of the smaller dispersed phase is displayed.The two panels qualitatively show the e↵ect of the sinusoidal disturbance on the suspension: the large particles tend to separate from the small ones, forming parallel layers normal to the shearwise (i.e.y) direction.The number of layers of large particles seem to be correlated with the wavenumber of the disturbance, while the role of the amplitude appears to be related to a threshold that activates suspensions analysed in the parametric study are located within the shear-thickening region of the flow curve, i.e., γ/ γ0 = 1, thus, the contributions to the shear-stress related to the contacts and the hydrodynamics are dominant (with the former being more important), while the electrochemical potentials play no effective role.Within the parameters chosen in this work, we observe that the reduction of η r (for the cases with A ≥ 20 and n ≥ 1) is mostly due to a strong decrease of the contacts contribution to the shear-stress, while the hydrodynamics remains almost unchanged.
Worth noticing is the case with A = 20 and n = 16, where also the hydrodynamics contribution is significantly reduced compared to the other cases analysed.
To highlight the origin of the relative viscosity reduction, we show in Figure 3, panels a and b, the instantaneous snapshots of the suspensions studied when a statistically steady-state is reached.In particular, the two panels are organised in a 5 × 2 (or 6 × 2) matrix as follows: from left to right, the amplitude (or the wavenumber) increases, from A = 0 (n = 0) to A = 40 (n = 16); in the top row, the distribution of the larger dispersed phase is shown while, in the bottom row, the location of the smaller dispersed phase is displayed.The two panels qualitatively show the effect of the sinusoidal disturbance on the suspension: the large particles tend to separate from the small ones, forming parallel layers normal to the shearwise (i.e.y) direction.The number of layers of large particles seem to be correlated with the wavenumber of the disturbance, while the role of the amplitude appears to be related to a threshold that activates the demixing of the phases, as previously postulated from the mean values of the relative viscosity (i.e., Figure 2).
Quantitatively, this becomes clear when showing the small and large particles concentration along the y-coordinate, as pictured in Figure 4 (note that the symmetry of the profiles is discussed in the Supplementary Material).The figure shows, within each row, the volume fraction of the small particles ϕ s (y), that of the large particles ϕ l (y) and the total volume fraction ϕ(y), as a function of the amplitude of the disturbance with wavenumber n = 1 (top row), and as a function of the wavenumber n, fixing the amplitude to A = 20 (bottom row).Note that the local volume fractions oscillate around the values ϕ s = 0.4, ϕ l = 0.1 and ϕ = 0.5, highlighted by the dashed lines; therefore, the labels on the abscissas simply indicate the aforementioned values shifted by an offset to show all the curves on a single the demixing of the phases, as previously postulated from the mean values of the relative viscosity (i.e., Figure 2).Quantitatively, this becomes clear when showing the small and large particles concentration along the y-coordinate, as pictured in Figure 4 [note that the symmetry of the profiles is discussed in the Supplementary Material 64].The figure shows, within each row, the volume fraction of the small particles s (y), that of the large particles l (y) and the total volume fraction (y), as a function of the amplitude of the disturbance with wavenumber n = 1 (top row), and as a function of the wavenumber n, fixing the amplitude to A = 20 (bottom row).Note that the local volume fractions oscillate around the values s = 0.4, l = 0.1 and = 0.5, highlighted by the dashed lines; therefore, the labels on the abscissas simply indicate the aforementioned values shifted by an o↵set to show all the curves on a single graph.From the top row of Figure 4, we observe the progressive formation of two regions with high large particles the demixing of the phases, as previously postulated from the mean values of the relative viscosity (i.e., Figure 2).Quantitatively, this becomes clear when showing the small and large particles concentration along the y-coordinate, as pictured in Figure 4 [note that the symmetry of the profiles is discussed in the Supplementary Material 64].The figure shows, within each row, the volume fraction of the small particles s (y), that of the large particles l (y) and the total volume fraction (y), as a function of the amplitude of the disturbance with wavenumber n = 1 (top row), and as a function of the wavenumber n, fixing the amplitude to A = 20 (bottom row).Note that the local volume fractions oscillate around the values s = 0.4, l = 0.1 and = 0.5, highlighted by the dashed lines; therefore, the labels on the abscissas simply indicate the aforementioned values shifted by an o↵set to show all the curves on a single graph.From the top row of Figure 4, we observe the progressive formation of two regions with high large particles graph.From the top row of Figure 4, we observe the progressive formation of two regions with high large particles concentration ϕ l , with the peaks expanding monotonically by increasing the amplitude of the disturbance.On the other hand, in the bottom row of Figure 4, as n grows we observe an increase of the number of layers where the large particles are collected, thus suggesting that the number of accumulation regions is controlled by the wavelength of the disturbance with the amplitude controlling only the amount of the relative concentration.This quantitative analysis corroborates what observed in the snapshots shown in Figure 3. Panel d and e of Figure 4 also show that, for small wavenumbers, i.e., n ≤ 8, the large particles accumulate into well-separated layers while, for large wavenumbers, i.e., n > 8, the small ones do.This is a consequence of the size of the particles compared to the wavelength of the perturbation: in the former case, the large particles have a diameter comparable to or smaller than the wavelength of the sinusoidal wave, defined as Λ = 1/(nκ 0 ), while in the latter case, the small ones do.Indeed, when the wavelength of the disturbance is much smaller than the particle size, its effect is filtered out by the particle size, with the particles feeling the disturbance as noise in the limit of Λ ≪ a 1 , a 2 .This becomes clear when looking at the local distributions of the large particles for n = 8 and n = 16 in panel e of Figure 4: the two curves show the same number of peaks ( 16), ration l , with the peaks expanding monotonically by increasing the amplitude of the disturbance.and, in the bottom row of Figure 4, as n grows we observe an increase of the number of layers w rticles are collected, thus suggesting that the number of accumulation regions is controlled by the wa isturbance with the amplitude controlling only the amount of the relative concentration.This qua corroborates what observed in the snapshots shown in Figure 3. Panel d and e of Figure 4 also show avenumbers, i.e., n  8, the large particles accumulate into well-separated layers while, for large wave 8, the small ones do.This is a consequence of the size of the particles compared to the waveleng ation: in the former case, the large particles have a diameter comparable to or smaller than the wave soidal wave, defined as ⇤ = 1/(n 0 ), while in the latter case, the small ones do.Indeed, when the wa isturbance is much smaller than the particle size, its e↵ect is filtered out by the particle size, with the he disturbance as noise in the limit of ⇤ ⌧ a 1 , a 2 .This becomes clear when looking at the local dist rge particles for n = 8 and n = 16 in panel e of Figure 4: the two curves show the same number of pe that the e↵ect of the wavelength on the large particles had reached a maximum and saturated.An aximum value of the wavenumber felt by particles of radius a i (with i = 1, 2) can be easily compute or the large particles we consider gives n max ⇡ 7, comparable to n = 8.The wavenumber n max is tween L and d because of the points of accumulation of the particles, which are located in the proxim ar values, as it can be seen in meaning that the effect of the wavelength on the large particles had reached a maximum and saturated.An estimate of the maximum value of the wavenumber felt by particles of radius a i (with i = 1, 2) can be easily computed as which for the large particles we consider gives n max ≈ 7, comparable to n = 8.The wavenumber n max is half the ratio between L and d because of the points of accumulation of the particles, which are located in the proximity of the zero-shear values, as it can be seen in Figure 5.In particular, panel a of the figure shows a snapshot at statistically steady-state of the large particles for the case A = 20 and n = 4, while in panel b the mean local volume fraction of the large particles (left side) together with the shape of the local shear-rate (right side) are plotted.From panel a and b, we can see how the large particles accumulate in parallel layers around the locations of the zero-shear regions (dotted horizontal lines in panel b).On the other hand, when considering disturbances with high wavenumbers, the smaller particles are reorganizing in parallel layers, as can be observed in panel c and d of Figure 5.The two panels mirror panel a and b for the case A = 20 and n = 16.For this case, it is clear that the small particles collect around the zero-shear regions.Next, to understand the mechanism of the particles accumulation, in Figure 6 we show the relevant components of the local relative mean stress tensor for the case A = 20 and n = 4 at statistically steady-state.In particular, panel a and b show the local relative viscosity (i.e., the local shear-stress ratio) and the local relative normal components of the stress tensor, respectively, along the shearwise direction y.Note that the integral value of the local shear-stress ratio along the shear-wise direction y corresponds to the mean relative viscosity.In these panels, the dotted horizontal lines outline the zeros of the total shear-rate profile, i.e., γ(y) = 0.The profile of the local shear-stress ratio promptly highlights that the suspension is sheared in two opposite directions (opposite signs of η r (y)) following the shape of the imposed shear-rate, with regions of zero shear-stress in the proximity of the large particles collection points (i.e., γ(y) = 0).Focusing on the normal stresses, instead, it is clear that, in the regions where the large particles accumulate, there are peaks of high-compressive stresses along the three directions.For σ 11 (y) (red line) and σ 22 (y) (green line) those peaks are enclosed by two adjacent local maxima, whose distance is comparable with the diameter of the large particles.A zoom on one of the regions enclosing the three aforementioned peaks (one minimum and two maxima) is shown in panel c of Figure 6.In particular, only the normal stress along the shearwise direction σ 22 (y) is reported.Furthermore, in the figure we also sketch the force f y that the gradients of the normal stress σ 22,y generate on the particles: we notice how the gradients enclosing the minimum peak tend to push away the particles from the zero-shear regions, while the two local maxima have an opposite effect.From the sketch, it is clear that small particles are effectively pushed out of the zero-shear regions, while the large particles are trapped by the two small stabilising peaks of σ 22 (y).Indeed, the effect of the stronger peak in the proximity of the zero-shear regions is filtered out by the large particles due to its size.This explains the mechanism of the particles migration for this configuration: the of the large particles.A zoom on one of the regions enclosing the three aforementioned peaks (one minimum and two maxima) is shown in panel c of Figure 6.In particular, only the normal stress along the shearwise direction 22 (y) is reported.Furthermore, in the figure we also sketch the force f y that the gradients of the normal stress 22,y generate on the particles: we notice how the gradients enclosing the minimum peak tend to push away the particles from the zero-shear regions, while the two local maxima have an opposite e↵ect.From the sketch, it is clear that small particles are e↵ectively pushed out of the zero-shear regions, while the large particles are trapped by the two small stabilising peaks of 22 (y).Indeed, the e↵ect of the stronger peak in the proximity of the zero-shear regions is filtered out by the large particles due to its size.This explains the mechanism of the particles migration for this configuration: the small particles migrate away from the zero-shear regions and, in turn, push the large particles in the accumulation points, eventually locking them there by acting with compressive forces on the extrema of the large particles.One of the reasons that trigger the phase separation and the demixing process can be sought in the strength of the forcing term that imposes positive and negative shear-rates to the suspension.To investigate this, we carry out a further parametric study where we choose the case with the streamwise wave having n = 4 (that showed an e cient demixing) and we vary the amplitude of the wavy disturbance.Panel a of Figure 7 reports the values of the relative viscosity collected at statistically steady-state for the cases with amplitude A = [0.5, 2.5, 5, 20], compared to the reference case, i.e. the suspension subject to a plain shear flow (black marker with A = 0).The figure shows a similar non-linear trend seen for the cases with n = 1, with a strong reduction of the relative viscosity between A = 2.5 and A = 5.In panel b, instead, we report the profiles of the applied shear-rates along y.Note that the case with A = 5 is the first one to cross the zero-shear line (dashed black line).This suggests that, when the disturbance is high enough to produce negative shear-rates, the demixing is activated causing a sudden reduction of ⌘ r .Note that, being the total shear-rate ˙ (y) = ˙ + A ˙ a 0 n 0 cos (n 0 y), either a disturbance with a strong amplitude or one with small particles migrate away from the zero-shear regions and, in turn, push the large particles in the accumulation points, eventually locking them there by acting with compressive forces on the extrema of the large particles.One of the reasons that trigger the phase separation and the demixing process can be sought in the strength of the forcing term that imposes positive and negative shear-rates to the suspension.To investigate this, we carry out a further parametric study where we choose the case with the streamwise wave having n = 4 (that showed an efficient demixing) and we vary the amplitude of the wavy disturbance.Panel a of Figure 7 reports the values of the relative viscosity collected at statistically steady-state for the cases with amplitude A = [0.5, 2.5, 5, 20], compared to the reference case, i.e. the suspension subject to a plain shear flow (black marker with A = 0).The figure shows a similar non-linear trend seen for the cases with n = 1, with a strong reduction of the relative viscosity between A = 2.5 and A = 5.In panel b, instead, we report the profiles of the applied shear-rates along y.Note that the case with A = 5 is the first one to cross the zero-shear line (dashed black line).This suggests that, when the disturbance is high enough to produce negative shear-rates, the demixing is activated causing a sudden reduction of η r .Note that, being the total shear-rate γ(y) = γ + A γa 0 nκ 0 cos (nκ 0 y), either a disturbance with a strong amplitude or one with a large wavenumber (or a combination of the two) may cause the demixing of the dispersed phases.
The final question of this work concerns the extension of the demixing process seen within the manuscript for the other rheological regimes of the dense suspensions, i.e. for any point belonging to the whole flow-curve (η r as a function of the shear-rate).To introduce the shear-rate dependency, we tune the amplitude of the electrochemical as in Equation (5) [see also 11].To determine the Hamaker constant and the intensity of the repulsive forces separately, we set the ratio |F R |/|F A | = 9 when two particles are at contact d ij = 0. Panel a of Figure 8 shows the values of the relative viscosity collected at statistically steady-state for the suspensions immersed in a plain shear flow (markers with the black edge) and for the suspensions subject to the sinusoidal disturbance with A = 20 and n = 4 (markers with the pink edge), varying the shear-rate γ/ γ0 .Panel b and c, instead, show the decomposition of the relative viscosity in the three main contributions (as seen in Figure 2 already) for the suspension subject to the plain shear-flow and the sinusoidal shear-flow, respectively.As it can be seen, the relative viscosity is always lower when a sinusoidal shear is imposed to the suspension for any γ/ γ0 , meaning that the demixing is always activated, despite the different nature of the contributions that govern the rheology of the suspensions in the various regimes, as it can be seen from the histograms.In particular, while for high γ/ γ0 the reduction of the relative viscosity is due to a reduction of the contacts, as discussed before, for low γ/ γ0 this is due to a reduction of the electrochemical interactions.This suggests that the segregation phenomenon in this case is driven by the response of the suspension to the high shear introduced by the sinusoidal wave, that in turns activates the selective process carried out by the normal stresses.The effect of the different contributions controlling the particle dynamics in the various regimes, however, can be spotted from the instantaneous snapshots of the large particles at statistically-steady state, as shown in the bottom row, i.e. panels d-f, of Figure 8.In particular, panel d-f show the instantaneous distributions of the large particles for γ/ γ0 = 10 −5 , γ/ γ0 = 10 −2 and γ/ γ0 = 10 0 , respectively.These three particular scenarios have been selected based on panel c of Figure 8 to show the effect of the different contributions: in fact, in the configuration of panel d, the dominant contribution is given by the electrochemical potentials (i.e., forces interactions along the centre-to-centre direction); in that of panel e, the dominant contribution is given by the hydrodynamics, while in the condition of panel f by the normal and frictional contacts.From the figure, it can be seen that when the contributions to the stress tensor are dominated by forces along the centre-to-centre direction (such as in panel d), the large particles accumulate in the region with positive shear-rate (see the pink profile in panel b of Figure 7) while, when forces with tangential direction start appearing (gradually increasing from panel e to panel f ), the large particles tend to stabilise around the zero-shear regions only.

IV. DISCUSSION
We have performed numerical simulations of dense binary suspensions with a large dispersion ratio, under a nonuniform shear flow, consisting of the combination of a linear shear and a sinusoidal disturbance.First, we find that the only disturbance altering the rheology of the suspension is the one acting on the streamwise component of the velocity in the shear-plane.By varying the amplitude and wavenumber of the perturbation, we discover that the demixing process can be triggered when the local shear-rate becomes negative.When the phases separate, the relative viscosity reduces in the whole flow curve.We explain the process by analysing the full stress tensor of the suspension: we find that large particles are locked in the accumulation points by compressive forces, while the small ones are pushed out by the gradients of the normal stresses.With this analysis, we highlight the importance of the full stress tensor to unravel the complex rhelogical behaviour of dense suspensions.
The phenomenon of particle separation was examined throughout the entire rheological curve under optimal conditions (A = 20 and n = 4).Regardless of the shear rate applied to the suspension, the suspension undergoes phase-separation, although the relative significance of the factors contributing to shear stress have a different nature.On the right side of the flow curve, contact interactions dominate, while in the shear-thinning region, electrochemical potentials take precedence.This demixing effect leads to the accumulation of larger particles (for n ≤ 4 in our study) in different regions along the direction of shear: the zeros of the shear-rate for conctact-dominated suspesions, positive shear-rate for electochemical potentials-dominated interactions.
It is important to underline that, to properly understand the demixing, considering the finite size of the particles is essential to filter out and select the proper information from the gradients of the normal stress tensor.Therefore, caution should be used when employing point-size particles models to explain or analyse such phenomena.
The outcomes of our analysis have far-reaching implications within several scientific disciplines, specifically rheology, microfluidics, biofluidics, and granular flows.By delving into these areas, our research contributes to a deeper understanding of the behaviour of various materials and fluids under different conditions.
In the field of rheology, our findings provide a foundational basis for developing a comprehensive rheological model that accurately captures the complexities of real-world situations.Traditionally, rheological models have overlooked structural inhomogeneity or natural disturbances, which can have substantial effects on material behaviour.However, our study addresses this limitation by considering these factors, thus offering a more realistic depiction of rheological phenomena.Moreover, our research extends its impact to the realm of microfluidics, where the manipulation and control of fluids at the microscale are essential.By incorporating our insights into microfluidic systems, researchers and engineers can gain a better understanding of how spatial disturbances influence the fluid flow and behaviour, enabling them to design more efficient and precise microfluidic devices.Additionally, our research contributes to the study of granular flows.By considering spatial disturbances and their effects on granular materials, we offer a more comprehensive framework for modelling and analysing granular flow behaviour, thereby facilitating improved process design and optimisation.
An exciting future direction for our study involves exploring the coupling of spatial disturbances with the more commonly studied temporal disturbances.By incorporating both spatial and temporal factors into our model, we can take a significant stride towards simulating and understanding real-world conditions more accurately.This advancement would enable researchers to investigate how disturbances propagate and interact over time, providing a more holistic understanding of complex systems and further bridging the gap between theoretical models and real-world scenarios.In this appendix, we report the expressions of the forces acting on the particles considered in this study, and provide more details on the numerical method.In particular, in our model we consider the hydrodynamics, inelastic contacts and electrochemical potentials, and indicate them by the superscripts H, C and E,respectively.
Concerning the hydrodynamics, dense suspensions of rigid particles immersed in a low-Reynolds-number flow are subject to a Stokes drag and a pair-wise, short-range lubrication force [1] caused by the reaction of the fluid squeezed in the narrow gaps between two neighbouring particles.This contribution is modelled in our software as a linear relationship between the forces and the velocities, F H = −R(u − U ∞ ), where R is a resistance matrix that is obtained by neglecting the far-field effects and considering only the dominant near-field (threshold gap between particles δ lub < a 0 ) divergent elements coming from the squeeze, shear and pump modes [1,2].More details on the coefficients of the resistance matrix adopted can be found in Monti et al. [3].
The contacts, instead, are replicated with the stickand-slide method [4], that mimics inelastic, frictional contacts between a pair of slightly overlapping particles with spring-dashpot systems placed in the normal (centre-to-centre) and tangential directions.Considering two particles, i and j, the stick-and-slide model can be written as where the vectors n and t, together with the subscripts n and t, indicate the directions of the force (i.e.normal and tangential) and the parameters k n , k t and γ n are the spring constants and the normal damping constant, respectively.The velocity u n,ij is the projection of the relative velocity vector between the particles ith and jth projected in the normal direction, while δ ij and ξ ij represents the displacement of the normal and tangential spring.The latter is an integral length that considers the history of the contact [4] and its value is such that the tangential force satisfies the Coulomb's law , where µ C is the friction coefficient (in this work, µ C = 0.5 for all the suspensions simulated).
Finally, the contribution from the electrochemical potentials is modelled as a combination of a distancedecaying repulsive force (with the distance-decay length scale set by the screening length) and an attractive force expressed in the van der Waals form [5,6], In Equation (2a), F 0 is the amplitude of the repulsive force, a = 2a i a j /(a i + a j ) the harmonic radius and L S is the screening length that controls the decay of the force.
In Equation (2b), instead, H A is the Hamaker constant and ε is a small regularisation term (generally ϵ ≈ 0.01a) introduced to avoid the singularity at contact, when the surface distance between the particles d ij = 0.All these contributions induce mechanical stress in the suspension and the stress tensor σ can be evaluated through the stresslet theory [7].From σ, all the rheological quantities (relative viscosity η r , first and second normal stress differences N 1 and N 2 and pressure Π) can be extrapolated.
The physics of the suspension described by the Newton-Euler equations is dominated by three nondimensional groups obtained by applying the Buckingham Pi theorem and choosing as fundamental quantities the properties that describe the hydrodynamic force, i.e.F H ∝ η 0 a 2 0 γ.The non-dimensional groups are listed in the main manuscript.In equation (4) of the manuscript, k n approximates the equivalent stiffness of the spring-dashpot systems, with the additional constraints k t ≪ k n and γ n γ/k n ≪ 1 [1].In equation ( 5) of the manuscript, we defined a new shear-rate, γ0 = F E /η 0 a 2 0 , being in F E , the dominant contribution is given by the repulsive forces, with the choice F R = 0.9F E and F A = 0.1F E .This non-dimensional group represents the additional time-scale introduced by the electrochemical contribution and it can be thought of as the equivalent of the Péclet number for the Brownian-suspensions; arXiv:2206.07916v2[cond-mat.soft]8 Oct 2023 thus, tuning the module of the repulsive force F R and the Hamaker constant appropriately, a shear-rate dependency can be imposed to the suspension.Note that in the manuscript, at times we omitted the subscript ij for the sake of readability.
From a computational viewpoint, the governing equations (2) in the manuscript, are discretized and advanced in time with the modified velocity-Verlet explicit scheme [8].The modified velocity-Verlet scheme is second-order accurate in time and, being explicit, it requires a timestep ∆t able to capture the shortest dynamics of the systems, typically established by the stiffness of the contacts; a good indicator of the time-step required can be obtained by building the time constant from the normal spring-dashpot system, i.e. ∆t γ = γ n γ/k n .Finally, as already mentioned above, the computational domain is a cubic box of size L 3 .The shear-rate γ is applied along the y direction, with the Lees-Edwards boundary conditions [9] preserving the ideality of the flow, removing the effect of the walls.
The set of parameters for carrying out these simulations are listed in Table I.Note that, for the sake of simplicity and completeness, a 0 = a 1 = 1 L, γ = 1 T −1 and ρ p = 1 M L −3 .Here L, T and M are the typical scales for length, time and mass, respectively.From Table I, all the parameters used in the simulations can be retrieved: k n = kη 0 a 0 γ; k n ; γ n = 10 −7 k n γ ; F E = (6π)η 0 a 2 0 γ0 ; therefore, Hamaker constant; Note that, for F E we multiply the non-dimensional group for the value 6π, that is the value that premultiplies the viscous force in the Stokes drag of a sphere [see 1, where the Stokes drag of a sphere has been chosen as reference value].
SUPPLEMENTARY MATERIAL TO: RESULTS

Symmetry of the local volume fraction profiles
Concerning the possible lack of symmetry of the profiles in Figure 4 of the main manuscript, we prove here that the lack of averaging and the finite number of particles (only 600 large particles) are to blame.In particular, we considered three initial conditions to the 8 cases with the sinusoidal perturbation, obtained by sampling the simulation with no perturbation (i.e.A = 0 and n = 0) at different time instants.The cases with identical amplitude and wavenumber are then ensemble averaged, together with 900 instantaneous fields per run sampled after convergence.Figures 1 and 2 compare the profiles with and without ensemble average respectively.The improvement of the symmetry is clearly visible.Therefore, the lack of perfect symmetry in the volume faction profiles is related to the lack of averaging only, together with the finite number of particles.
Finally, we also prove the correctness of our implementation of the Lees-Edwards boundary conditions.We selected the run with A = 20 and n = 1 (red curve in Figures 1 and 2) and we averaged its volume fraction distributions with an identical set-up whose streamwise velocity profile is shifted of a value u = 30a 0 γ.In case of a wrong implementation of the boundary conditions, the resulting local volume fraction should differ, while the red profiles in Figure 2

FIG. 1 :
FIG.1: Top row, panels a-c: sketch of the sinusoidal disturbance analysed in the respective columns.Bottom row, panels d-f : time-history in terms of strain of the relative viscosity ⌘ r .The colours distinguish the di↵erent cases, with the black solid lines referring to the values of ⌘ r for the particles suspended in a plain shear-flow.

FIG. 2 :
FIG. 2: Panel a and c: time-history in terms of strain of the relative viscosity ⌘ r for the disturbance u = A ˙ a 0 sin (n 0 y).Panel b and d: contributions to the mean relative viscosity ⌘ r accumulated at statistically-steady state of the contacts, hydrodynamics and electrochemical potentials (the respective colours are shown in the legend within the panels).The panels in the top row, i.e. panel a and (b), show the trends of ⌘ r varying the amplitude A of the disturbance with n = 1.The panels in the bottom row, i.e. panel c and d, show the trends of ⌘ r as a function of the wavenumber n, with the amplitude of the sinusoidal wave fixed at A = 20.The colorbars on the left of panel a and c indicate the colour-scheme adopted for the amplitudes (top) and wavenumbers (bottom) analysed; the black solid lines refer to the reference case of the suspension subject to a plain shear-flow.

FIG. 2 :
FIG. 2: Panel a and c: time-history in terms of strain γ of the relative viscosity η r for the disturbance u = A γa 0 sin (nκ 0 y).Panel b and d: contributions to the mean relative viscosity η r accumulated at statistically-steady state of the contacts, hydrodynamics and electrochemical potentials (the respective colours are shown in the legend within the panels).The panels in the top row, i.e. panel a and (b), show the trends of η r varying the amplitude A of the disturbance with n = 1.The panels in the bottom row, i.e. panel c and d, show the trends of η r as a function of the wavenumber n, with the amplitude of the sinusoidal wave fixed at A = 20.The colorbars on the left of panel a and c indicate the colour-scheme adopted for the amplitudes (top) and wavenumbers (bottom) analysed; the black solid lines refer to the reference case of the suspension subject to a plain shear-flow.

FIG. 3 :
FIG. 3: Panel a: instantaneous snapshots of the suspensions at = 20 varying the amplitude of the sinusoidal perturbation, fixing n = 1.Panel b: instantaneous snapshots of the suspensions at = 20 varying the wavenumber of the sinusoidal perturbation, fixing A = 20.The top row of the subfigures shows the distribution of the large particles, while the bottom row the respective distribution of the small particles of the suspensions.The colour-scheme reflects the one adopted in the top row of Figure 2.

FIG. 3 :
FIG. 3: Panel a: instantaneous snapshots of the suspensions at γ = 20 varying the amplitude of the sinusoidal perturbation, fixing n = 1.Panel b: instantaneous snapshots of the suspensions at γ = 20 varying the wavenumber of the sinusoidal perturbation, fixing A = 20.The top row of the subfigures shows the distribution of the large particles, while the bottom row the respective distribution of the small particles of the suspensions.The colour-scheme reflects the one adopted in the top row of Figure 2.

FIG. 4 :
FIG. 4: Profiles of the local volume fractions along the shearwise direction y.Panel a and d show the local volume of the small particles s (y); panel b and e show the local volume of the large particles l (y); panel c and f show the total local volume fraction (y).The panels in the top row, i.e. panel a-c, refer to the sinusoidal disturbance with n = 1 and increasing amplitude A, as shown in the respective legend on the left of the figure; the panels in the bottom row, i.e. panel d-f, refer to the cases with A = 20 and increasing n, as shown in the respective legend on the left of the figure.The black profiles in each panel refer to the reference case with no sinusoidal disturbance.Finally, the dashed lines refer to the values s = 0.4, l = 0.1 and = 0.5.Note that, the profiles for di↵erent A and n have been shifted horizontally by 0.2 for visual purpose.Within the Supplementary Material [64], we discuss the reasons of the lack of symmetry of the profiles.

FIG. 4 :
FIG. 4: Profiles of the local volume fractions along the shearwise direction y.Panel a and d show the local volume of the small particles ϕ s (y); panel b and e show the local volume of the large particles ϕ l (y); panel c and f show the total local volume fraction ϕ(y).The panels in the top row, i.e. panel a-c, refer to the sinusoidal disturbance with n = 1 and increasing amplitude A, as shown in the respective legend on the left of the figure; the panels in the bottom row, i.e. panel d-f, refer to the cases with A = 20 and increasing n, as shown in the respective legend on the left of the figure.The black profiles in each panel refer to the reference case with no sinusoidal disturbance.Finally, the dashed lines refer to the values ϕ s = 0.4, ϕ l = 0.1 and ϕ = 0.5.Note that, the profiles for different A and n have been shifted horizontally by 0.2 for visual purpose.

FIG. 5 :FIG. 6 :
FIG. 5: Panel a and c: instantaneous snapshot of one of the two dispersed phases of the suspension; panel b and d: plots of the local volume fraction along the shearwise direction y of the respective phase, alongside the local profile of the shear-rate applied to the suspension (right side).Panel a and b: case A = 20 and n = 4, with the snapshot of the distribution of the large phase in the suspension; panel c and d: case A = 20 and n = 16, with the snapshot of the distribution of the small phase in the suspension.The two plots in panel b and d are split by the vertical, black solid line in the middle of the figure; the other vertical straight lines represent, from left to right: the mean value of the volume fraction of the particles ( l = 0.1 or s = 0.4, pink/violet dashed line); ˙ (y) = 0 (blue solid line) and the mean value of the total shear-rate applied to the suspension (pink/violet dashed line).The dotted horizontal lines outline the y location of ˙ (y) = 0. Finally, the pink/violet circle in the graphs sketches a large/small particle of the suspension.

FIG. 5 :FIG. 5 :FIG. 6 :
FIG. 5: Panel a and c: instantaneous snapshot of one of the two dispersed phases of the suspension; panel b and d: plots of the local volume fraction along the shearwise direction y of the respective phase, alongside the local profile of the shear-rate applied to the suspension (right side).Panel a and b: case A = 20 and n = 4, with the snapshot of the distribution of the large phase in the suspension; panel c and d: case A = 20 and n = 16, with the snapshot of the distribution of the small phase in the suspension.The two plots in panel b and d are split by the vertical, black solid line in the middle of the figure; the other vertical straight lines represent, from left to right: the mean value of the volume fraction of the particles (ϕ l = 0.1 or ϕ s = 0.4, pink/violet dashed line); γ(y) = 0 (blue solid line) and the mean value of the total shear-rate applied to the suspension (pink/violet dashed line).The dotted horizontal lines outline the y location of γ(y) = 0. Finally, the pink/violet circle in the graphs sketches a large/small particle of the suspension.

FIG. 6 :
FIG. 6: Panel a: local distribution of the mean shear-stress ratio η r (y) along the shearwise direction y.The vertical dotted line highlights η r (y) = 0; the two circles sketch the two types of particles dispersed in the suspension.Panel b: local distribution of the mean normal stresses, σ 11 (y) (red line), σ 22 (y) (green line) and σ 33 (y) (blue line), along the shearwise direction y.Panel c: zoom of σ 22 (y) in one of the regions of panel b bounded by the dashed lines.The location ŷ/a 0 corresponds to the zero-shear location γ(y) = 0.The two circles represent the particles of the two phases of the suspension.The arrows show the directions of the gradients of the normal stress σ 22,y and the respective forces they generate on the particles.In all the panels of the figure, the horizontal dotted lines indicate the regions at γ(y) = 0.The case considered has a sinusoidal perturbation with A = 20 and n = 4.

7 :
Wavenumber n = 4. Panel a: mean relative viscosity ⌘ r at statistically steady-state for the distu ˙ a 0 sin (4 0 y), varying the amplitude A of the disturbance.Panel b: profiles of the shear-rate applie ions along the shearwise direction y.The colour-scheme indicates the amplitude imposed and reflect used in the colorbar.The black dashed line is the zero shear-rate.

Figure 5 .FIG. 7 :
FIG.7: Wavenumber n = 4. Panel a: mean relative viscosity η r at statistically steady-state for the disturbance u = A γa 0 sin (4κ 0 y), varying the amplitude A of the disturbance.Panel b: profiles of the shear-rate applied to the suspensions along the shearwise direction y.The colour-scheme indicates the amplitude imposed and reflects the one used in the colorbar.The black dashed line is the zero shear-rate.

FIG. 8 :
FIG. 8: Top row, panels a-c.Panel a: mean relative viscosity ⌘ r at statistically steady-state for suspensions driven by a plain shear-flow (markers with a black edge) and for the suspension subject to the disturbance u = 20 ˙ a 0 sin (4 0 y) (markers with a pink edge) as a function of the shear-rate ˙ / ˙ 0 .Panel b and c: contributions to the mean relative viscosity ⌘ r accumulated at statistically-steady state of the contacts, hydrodynamics and electrochemical potentials (the respective colours are shown in the legend within the panels), as a function of the shear-rate ˙ / ˙ 0 , for the suspension driven by a plain shear-flow (panel b) and for the suspension subject to the sinusoidal disturbance (panel c).Bottom row, panels d-f : instantaneous snapshots of the large particles of the suspensions at statistically-steady state varying ˙ / ˙ 0 .Panel d: ˙ / ˙ 0 = 10 5 ; panel e: ˙ / ˙ 0 = 10 2 ; panel f : ˙ / ˙ 0 = 10 0 .The colours in panel a and panels d-f refer to the colormap in the figure.Note that the colour pink can be changed to black in panel a to show the non-perturbed case.The suspensions reported are subject to the sinusoidal perturbation with A = 20 and n = 4.

FIG. 8 :
FIG. 8: Top row, panels a-c.Panel a: mean relative viscosity η r at statistically steady-state for suspensions driven by a plain shear-flow (markers with a black edge) and for the suspension subject to the disturbance u = 20 γa 0 sin (4κ 0 y) (markers with a pink edge) as a function of the shear-rate γ/ γ0 .Panel b and c: contributions to the mean relative viscosity η r accumulated at statistically-steady state of the contacts, hydrodynamics and electrochemical potentials (the respective colours are shown in the legend within the panels), as a function of the shear-rate γ/ γ0 , for the suspension driven by a plain shear-flow (panel b) and for the suspension subject to the sinusoidal disturbance (panel c).Bottom row, panels d-f : instantaneous snapshots of the large particles of the suspensions at statistically-steady state varying γ/ γ0 .Panel d: γ/ γ0 = 10 −5 ; panel e: γ/ γ0 = 10 −2 ; panel f : γ/ γ0 = 10 0 .The colours in panel a and panels d-f refer to the colormap in the figure.Note that the colour pink can be changed to black in panel a to show the non-perturbed case.The suspensions reported are subject to the sinusoidal perturbation with A = 20 and n = 4.

FIG. 1 :FIG. 2 :
FIG. 1: Profiles of the local volume fractions along the shearwise direction y.Panel a and d show the local volume of the small particles ϕ s (y); panel b and e show the local volume of the large particles ϕ l (y); panel c and f show the total local volume fraction ϕ(y).The panels in the top row, i.e. panel a-c, refer to the sinusoidal disturbance with n = 1 and increasing amplitude A, as shown in the respective legend on the left of the figure; the panels in the bottom row, i.e. panel d-f, refer to the cases with A = 20 and increasing n, as shown in the respective legend on the left of the figure.The black profiles in each panel refer to the reference case with no sinusoidal disturbance.Finally, the dashed lines refer to the values ϕ s = 0.4, ϕ l = 0.1 and ϕ = 0.5.Note that, the profiles for different A and n have been shifted horizontally by 0.2 for visual purpose.Version with no ensemble average.

TABLE I :
Set of parameters used for simulating the behaviour of dense suspensions varying the shear rate.Here, a 0 = a 1 is the reference length scale.