A fast and efficient tool to study the rheology of dense suspensions

The physics of granular particles suspended in a fluid has become of fundamental importance in the last few decades due to their heavy involvement in many natural and industrial processes. For instance, clarifying the behavior of blood flows is a high-priority topic since it can potentially lead to the diagnosis and prevention of cardiovascular diseases.¹ From a geophysical viewpoint, phenomena like sediment transport are highly attractive since dense suspensions contribute to transporting nutrients and providing ecological habitats while shaping the surrounding landscape.² Industries working with paints or adhesives, instead, are interested in the control of the flow to guarantee the required performance and quality of the final products approaching the market.³ 

The physics of the processes listed above is dominated by the interactions between the particles and fluid that determine the flowability properties of the granular suspensions. Unscrambling the cumbersome microscopic mechanisms of the interactions between the particles shaping those materials and their flowing behavior has become a high-priority topic, with a particular need for finding a physically meaningful constitutive law.⁴ Exceptionally challenging is the rheological behavior shown by dense suspensions, i.e., suspensions with a similar volume fraction of particles and liquid. In fact, under stress, dense suspensions may show rheological behaviors that span from a yielded behavior to discontinuous shear-thickening.^5–8 

Despite the practical relevance, a few analytical studies dealing only with very dilute⁹ and semi-dilute¹⁰ granular suspensions have been successfully carried out. The major mathematical difficulties stand in the lack of proper analytical formulations that deal with multi-body interactions, especially for very large volume fractions, where enduring direct contacts becomes fundamental.⁶ Conversely to theoretical approaches, numerical algorithms are ideal for tackling the numerous particles interacting with each other. Numerical techniques, spanning from molecular dynamics^11,12 to the solution of the linear Stokes equations for every suspended particle,¹³ stormed the field of granular suspensions. Recently, since the work done in the seminal paper by Cundall and Strack,¹⁴ the discrete (or distinct) element method (DEM) for studying dense suspensions and granular flows has become the standard for the flow analysis. Modern approaches to DEM reduce the description of the particle–fluid and particle–particle interactions to minimal models able to correctly address the macroscopic behavior of the flow without excessive loss of information. On this fashion, simplified relationships for short-range interactions have been proposed, e.g., for hydrodynamic interactions^15,16 and for adhesive–repulsive contributions (a satisfactory overview can be found in the works of Refs. 16–19). Concerning contacts, the most implemented technique is the soft-particle model, where contacts among several particles and adhesion effects are handled.¹⁸ In particular, Luding²⁰ developed a minimal model that, based on the original work by Cundall and Strack,¹⁴ addresses the surface roughness of the particles with tangential contacts replicating the sliding, rolling, and torsion resistances, with the option of reproducing both static and dynamic frictions. More advanced contact models contemplating non-linearity, hysteresis, and complex effects can be found in his work.²⁰ Many other models have been adopted in the field of dense suspensions and DEM;^8,21 the interested reader is referred to the reviews by Li et al.¹⁸ and Guo and Curtis²² for a complete overview.

The major drawback of DEM is the expensive computation of short-range interactions acting on every particle. The computational cost, if care is not taken, may quadratically increase with the number of particles considered, limiting the size of the domain analyzable or the length of the simulated process. To reduce the computational cost, algorithms of fixed-radius near neighbors search²³ have been adopted and become the standard in DEM implementations. A widely implemented method is related to the use of Verlet lists.¹² The Verlet list is a data structure that stores and tracks all particles within a cutoff distance of each other. The list can be naively built by checking the $N(N−1)/2$ distances [cost $O(N2)$] between particles or, with more modern techniques, through efficient cell-lists or tree-structures that largely reduce the computational cost [O(N) for cell-lists, $O(N log N)$ for tree structures].²⁴ In particular, the former consists of dividing the computational box into a Cartesian lattice (with spacing Δ that accommodates the particles within a pre-selected cutoff distance) and assigning the particles to the points of the lattice. After that, for each particle a distance-checking procedure is applied to all the other particles belonging to the 27 neighboring cells. This procedure, however, may bring to a high computational and memory cost when dealing with highly polydispersed particles, since the cutoff distance is usually calibrated with the largest particles; therefore, a cell may contain a large number of small particles that must be enquired for neighborhood. While computationally the extra-cost can be decreased by adopting a cutoff radius that suits the smallest particle of the set, and extending the size of the searching-stencil for the larger particles as done in LAMMPS, memory-wise the cost will remain of the order $O(N2)$⁠, since the cell-list has to spatially cover the whole simulation box.

To further reduce the computational cost, the Verlet-lists are updated every n time steps, calibrated such that within the n − 1 remaining steps the locations of the neighboring particles do not significantly vary.¹² This last point, however, is very delicate since updating the Verlet list too frequently is computationally expensive and, on the contrary, increasing the time-period of the update can introduce errors in the computation.

Here, following the work by Hoetzlein,²⁵ we implement a highly parallelisable method based on the cell-list approach that does not require the explicit computation of a Verlet list, thus zeroing its memory-cost, since the particles will be simply reordered within the Cartesian lattice in a new array by means of a very efficient counting-sort algorithm [computational cost O(N)]. Since the particles are ordered by the cell, it is straight-forward to look for neighbors within the adjacent cells. The counting-sort algorithm will be applied to the particles at every time iteration, avoiding any loss of information on the neighboring pairs. Within the framework of this manuscript, we propose and validate the new software that combines the cutting-edge algorithm just mentioned for the neighbors query with methods that satisfactorily deal with complex rheological behaviors observed in dense suspensions. For the latter problem, a rigorous mathematical formulation will be provided. Finally, the software is parallelized within the multi-threading architecture and is intended to be used without accessing the expensive computational power provided by supercomputers or costly Graphics Processing Unit (GPU) accelerators.

The manuscript is organized as follows. Section II deals with the mathematical formulation, focusing on the models adopted and detailing the parameters that govern the flow. In addition, the numerical integration scheme together with its implementation and parallel performance are described. In Sec. III, we show the validation of the code by studying the rheology of dense suspensions to ensure the reliability of the results. Some benchmark tests on the potentiality of the code are also proposed. Finally, in Sec. IV a summary of the code capabilities and novelty ends the manuscript.

The in-house software used in this work models a dense, non-Brownian suspension of quasi-inertialess, neutrally buoyant, quasi-rigid spherical particles in a shear-flow defined by the shear-rate, $γ̇$⁠, and the strain-rate tensor, $E∞$ (with the only non-zero quantities, $E12∞=E21∞=γ̇/2$⁠). The code tackles the Newton–Euler equations that govern the translational and rotational dynamics of the rigid particles,

where the subscript i indicates the particle $i∈[1,N]$⁠, N being the number of particles. In Eq. (1), $FiM$ and $TiM$ denote the force and torque applied to the center of mass of the ith particle, with mass m_i and inertia tensor $Ii$⁠. The translational and angular velocities are here denoted by the symbols $ui$ and $ωi$⁠, respectively. The superscript M refers to the nature of the forces applied to system (1), resulting from particle–particle and particle–flow interactions. In particular, here we consider the contributions arisen from hydrodynamics, inelastic contacts, and electro-chemical effects. Specifically, the right-hand-side of Eq. (1) can be written as

where the summations are performed on the number of particles neighboring the ith particle considered, and the superscripts H, C, and E stand for hydrodynamics, inelastic contacts, and electro-chemical effects, respectively.

Concerning the hydrodynamics of the system, dense suspensions of rigid particles immersed in a low-Reynolds-number flow are subjected to a Stokes drag and a pair-wise, short-range lubrication force.²⁶ The latter is caused by the relative motion of particles that squeeze the fluid flowing in the narrow gaps between them. Following the work of Ball and Melrose¹⁵ and Mari et al.,²⁶ to replicate those contributions we implement a linear relationship between forces and velocities,

where $n=(xj−xi)/|xj−xi|$ is the center-to-center unit vector (being $xi$ the Cartesian position of the particle i), and $n⊥,i(A)=2n×[(Ai+AiT)·n]$ is a vector perpendicular to n (being A_i the strain-rate tensor referred to the ith particle, $Ei∞$⁠, in this case), $ℝStokes$ is a diagonal matrix that models the Stokes drag forces and torques, including the Faxén laws for non-uniform shear-flows,²⁷ and $ℝLub$ and $ℝ′Lub$ are sparse matrices that approximate the far-field mobility matrices, dropping the low moment contributions and considering a pair resistance model with only divergent terms.¹⁵ In particular, those sparse matrices model the contribution of the lubrication through the squeeze, shear, and pump modes.²⁶ The detailed formulation of the matrices defining the hydrodynamics resistances is given in A. Figure 1 sketches the lubrication forces and torques acting on the ith particle (left) due to the presence of the jth particle (right). The dashed red circumference highlights the region of influence of the lubricant interactions, parameterized with δ_lub. Note that beyond this threshold, i.e., even at mid-range distances between the particles, the hydrodynamics is modeled only through the Stokes drag. We acknowledge that this is an approximation; however in the context of dense suspensions, it is acceptable since its contribution is minimal.

The contact contribution is modeled by the stick-and-slide method²⁰ that is able to describe frictional contacts. This method mimics the inelastic contacts with spring-dashpot systems that start operating when the surface distance between the spheres i and j, $dij=|xi−xj|−(ai+aj)$ (being a_i and a_j the particles radii), is negative (overlap). Considering two particles, i and j, the stick-and-slide model can be written as

where the Coulomb's law $|FijC, tan |≤μC|FijC,nor|$ must be satisfied (μ_C being the friction coefficient). The superscripts nor and tan indicate the directions of the force (i.e., normal and tangential) and the parameters k_n, k_t, and γ_n are the normal and tangential spring constants and the normal damping constant, respectively. The velocity $un,ij=n⊗n·(uj−ui)$ is the normal velocity vector, while $ξij$ represents the stretch of the tangential spring. The latter is computed following the procedure described by Luding:²⁰ 

where $T=I−n⊗n$ is the tangential projection tensor, $ut,ij=T·(uj−ui−(aiωi+ajωj)×n)$ is the tangential velocity vector, and t₀ is the time at the beginning of the collision between the particles i and j [with $ξij(t0)=0$]. A sketch of the model adopted in this work is shown in Fig. 2. The possibility of adding the rolling resistance to the model, as in Refs. 19 and 20, is described in B.

The last contribution considered in this work is forces of conservative nature that arise from the chemo-physical properties of the solid–fluid and solid–solid interactions and play a fundamental role in systems of dense suspensions.^19,28 In particular, the model, sketched in Fig. 3, addresses the contribution by adding an inter-particle, distance-decaying repulsion force,

In Eq. (6), F_ER is the module of the force, a₀ the reference radius, $a¯=2aiaj/(ai+aj)$ the harmonic radius, and λ_D is the Debye length that controls the decay of the force. In Eq. (7), instead, A is the Hamaker constant and ϵ is a small regularization term (generally $ϵ≈0.01a¯$⁠) introduced to avoid the singularity at contact, when d_ij = 0. To simplify the notation, below the repulsive and attractive effects may be gathered together in the electro-chemical contribution, $FE=FR+FA$⁠, marked by the E superscript.

The forcing contributions described above induce mechanical stress in the suspension that can be computed using the stresslets theory,²⁷ 

where $Σ$ is the bulk stress tensor, η₀ is the viscosity of the carrier fluid, and V is the volume of the simulated box. The stresslets due to the interactions, instead, are denoted by $S$ and are defined as

where the hydrodynamic contribution (H) has been separated into the Stokes (S) and lubrication (L) terms, i.e., $SijH=SiS+SijL$⁠. From the bulk stress tensor, $Σ$⁠, we can derive the quantities of rheological interest, i.e., the shear-stress σ, the normal stress differences N₁ and N₂, and the pressure Π, defined as $σ=Σ12, N1=Σ11−Σ22, N2=Σ22−Σ33$⁠, and $Π=−Tr(Σ)/3$⁠. Normalizing the shear-stress and the pressure with the properties of the carrier flow, we can define the relative viscosity and the dimensionless pressure as $ηr=σ/(η0γ̇)$ and $ηn=Π/(η0γ̇)$⁠, respectively.

Several parameters dominate the physics of the suspension described by System (1). We consider the translational equation of (1),

where we gathered the dimensional quantities (being the non-dimensional quantities marked with the superscript $*$⁠). In Eq. (10), $U=a0γ̇$ is the typical velocity in a shear-flow, T is the timescale introduced by the particles (in this case, the particles introduce no additional time-scales, therefore $T=1/γ̇$⁠), k_eq is the equivalent spring constant of the model represented in Fig. 2, and F is the module of the force introduced by the electro-chemical interactions. Applying the Buckingham Pi theorem, the parameters are reduced to a set of non-dimensional numbers that can be tuned to outline the model desired. Considering the dimensional quantities in Eq. (10), ${ρ,a0,η0,γ̇,keq,F}$ and choosing ${a0,η0,γ̇}$ as independent fundamental quantities, three non-dimensional groups emerge. The first one is the Stokes number,

where the constraint $St≪1$ is applied to remain in the inertialess regime. The non-dimensional stiffness, instead, governs the importance of the contacts contributions compared to the hydrodynamic term,

where $k̂≫1$ forces the particles to be rigid; note that, the equivalent stiffness constant,

has been approximated to k_n by adding additional constraints $kt≪kn$ and $γnγ̇/kn≪1$⁠.²⁶ Finally, the last non-dimensional group (the equivalent shear-rate) refers to the timescale introduced by the electro-chemical contribution,

where we defined a new shear-rate, $γ̇0=F(dij)/η0a02$⁠; thus, tuning F (i.e., adjusting F_ER and the Hamaker constant) appropriately, a shear-rate dependency can be imposed to the suspension.

Following the work by Ge et al.,^29,30 the governing equations (1) are advanced in time with the modified velocity-Verlet explicit scheme,³¹ 

where $xi(n), ui(n)$⁠, and $αi(n)$ denote the position (orientation), velocity (angular velocity), and acceleration (angular acceleration) vectors of the particle i, respectively, at time $t=nΔt$⁠, and $F$ resembles the right hand side of Eq. (1). The modified velocity-Verlet scheme is second-order accurate in time and, being explicit, it requires a time step $Δt$ that resolves the smallest timescale of the system (typically established by the stiffness of the contacts; a good approximation of the time step can be computed through the time constant that can be defined in the normal spring-dashpot system, $Δtγ̇=γnγ̇/kn$⁠). Different from other approaches, e.g., considering an over-damped system dropping the inertial terms,²⁶ this approach allows us to avoid the inversion of the resistance matrix of the lubrication forces.

The domain is a box of size $Lx×Ly×Lz$ (or L³ if cubic) with a shear-rate, $γ̇$⁠, applied along the y direction. At the edges of the domain, the Lees–Edwards boundary conditions³² are adopted to remove the effect of the walls,

The interaction between the particles represents the most expensive part of the algorithm from a computational point of view with a cost of $O(N2)$⁠, if a brute force algorithm is used. To overcome the problem, we adopted a highly parallelisable, fixed-radius near neighbors algorithm.^23,25 The algorithm is based on rearranging the particles on a coarser lattice (cell-list) by adopting an efficient counting sort algorithm²⁵ that requires a computational cost O(N) and few memory arrays of size N. The spacing of the cell-list Δ must suit the desired searching distance; a sketch of it is shown in Fig. 4. By considering the near points of the coarser grid in the search for neighbors, i.e., the blue-filled region in Fig. 4 (considering a 3 × 3 macro searching-cell) the computational cost dramatically decreases to O(cN), where $c≪N$ is the number of particles inquired as neighbors and depends on the lattice space Δ. The dependence of c on Δ is related to the number of particles that fit each macro-cell and an optimal choice of Δ is therefore related to the properties of the suspension. For instance, fixing $Δ=2amax$⁠, where a_max is the radius of the largest particle in the suspension, minimizes the number of macro-cells to be inquired for the search of neighbors. This is ideal for monodisperse suspensions where there may be at maximum one particle per cell and may look reasonable when considering polydispersed suspensions with not-very-high dispersion ratios. On the contrary, the choice may be detrimental for suspensions with $amin≪amax$ (a_min being the radius of the smallest particle), since it may result in the presence of a high number of particles within each macro-cell, with the limit of c = N when $2amax$ corresponds to the size of the domain. In this case, to reduce the number of particles within each searching cell, a better choice is to set $Δ<2amax$ and increase the stencil of neighboring cells.^24,25

As a drawback, the search shown in Fig. 4 implies that the forces and torques acting on the pair of particles (i, j) are computed twice, although this procedure guarantees a better load balance through the processes when the code is parallelized. To avoid the double counting of the forces keeping a good load distribution among the processes, we propose a symmetric search, where the neighbors are scanned only over half of the points of the macro searching-cell of the near Cartesian lattice, as sketched in Fig. 5, thus obtaining an additional significant speed-up by halving the computational operations. Note that this procedure does not require any particular symmetry of the problem.

Memory-wise, instead, the main problem arises from the computation of the tangential stretch, Eq. (5), where the history of all contacts must be traced, with a potential memory cost $O(N2)$⁠. To reduce the cost, we use a reduced matrix with size N × M, where M is an overestimation of the maximum number of contacts per particle based on the volume of the spherical shell surrounding the largest particle of the suspension, with thickness $Δa=2a0$⁠. The contacts, once established, are then recorded into a hash-map that charts them along with the indexes of the particles involved. Then, every array of the map (i, M) is sorted by the indexes, and the history of the contacts for the particle i is sought over the sorted array (i, M). This procedure has an actual memory cost O(N), adding little to the computational cost.

The parallelization of the solver is reached by implementing the multi-threading paradigm on a single node with a multiprocessing system, using the OpenMP directives. The high-scalability of the code, obtained on a 2-socket AMD EPYC 7702 machine (zen2 microarchitecture), is shown in Fig. 6. The scalabiltiy graph is obtained by fixing the domain properties and varying the number of processes used. The case considered is a bidispersed suspension with a high number of particles N = 65 536 to test also the memory-efficiency, volume fraction $ϕ=0.5$⁠, bidispersity ratio, and volume $a2/a1=3$ and $V2/V1=0.25$⁠, respectively. All the forcing contributions listed above in this section have been considered. Referring to the configurations shown in Fig. 6, the total time per iteration (on average) in seconds is reported in Table I. The table also shows the breakdown (in percentage) of the total time needed by the several parts sampled within the core of the software. In particular, t₁ corresponds to the time needed to update the kinematical properties of the particles with the modified velocity-Verlet scheme, t₂ to sort the particles within the corresponding cells (counting sort), t₃ to compute the interactions between neighbors, and t₄ to sort the hash-map of the frictional contacts list, while the amount of time needed for writing the outputs is negligible. The breakdown of the total time per iteration shows how the counting sort algorithm is convenient for reordering the particles within the searching cells, with its cost not exceeding the 6% of the total time in the configuration considered.

For readers interested in the algorithm and in benchmarking the software, the code is available at https://github.com/marco-rosti/CFF-Ball-0x.

In this section, the robustness and accuracy of the method proposed above for dense suspensions is investigated. We consider, first, a relatively small domain with particles suspended in a uniform, plain shear-flow, and we compare the data obtained with studies available in the literature varying the volume fraction and the shear-rate; we then increase the domain size to test the memory-efficiency of the software, and we study the behavior of a dense suspension subject to a wavy shear-flow, the wave acting on the stream-shear cross-plane. All initial conditions used in this work are generated with the sphere-packing software developed by Donev.³³ 

The reliability of the code has been tested both in jamming transition scenarios and in a shear-dependent behavior of rheological relevance.

For the jamming transition, we consider a cubic box containing a bidispersed suspension (⁠$a2/a1=1.4$⁠, in equal volumes $V2/V1=1$⁠) of 200 frictionless, rigid [see Eq. (12)] particles, with no electro-chemical interactions. The side length of the box accommodates the volume fraction $ϕ$ imposed, being $a0=a1=1$⁠. The parameters chosen echo the validation carried out by Ge and Brandt²⁹ and are listed in Table II. Figure 7 shows the comparison between the simulation results and the data available in the literature with a satisfying agreement between the computed relative viscosity, η_r, and the power-law fitting of Mari et al.²⁶ for frictionless particles (solid black curve).

We then activate the friction (⁠$μC=1$⁠), and we increase the number of particles to N = 1000 (see Table II). Figure 7 shows the comparison between the computed relative viscosity η_r and the data fitting obtained by Mari et al.²⁶ (solid red curve) and Cheal and Ness³⁴ (dashed red curve) with excellent agreement.

The shear-dependent behavior of the particles suspended in a uniform shear-flow is obtained by considering all contributions listed in Sec. II. At first, we will show the results obtained from a bidispersed dense suspensions (size ratio $a2/a1=1.4$⁠) with volume fraction $ϕ=[0.45,0.50,0.55]$⁠, to see if the software is able to catch the typical shear-dependent behavior.²⁶ The set of parameters for carrying out these simulations have been chosen to mimic those by Mari et al.²⁶ and are listed in Table III. Note that for these simulations, the electro-chemical contribution is purely repulsive (⁠$FR/FA=∞$⁠), as set in Mari et al.²⁶ Figure 8 (left panel) shows the trend of the relative viscosity η_r varying the shear-rate applied to the suspension $γ̇/γ̇0$⁠. In particular, we see that the three suspensions show the same qualitative behavior, with values of the relative viscosity substantially increasing (especially at higher shear-rates) with the volume fraction. From Fig. 8, we can notice that at low shear-rates the suspensions show the typical shear-thinning behavior with a similar slope in the three cases considered, as already observed by Mari et al.²⁶ At those shear-rates, the repulsive contribution dominates the mechanics of the suspensions (higher F^E, lower $γ̇/γ̇0$⁠), holding the particles off at a mid-range distance. On the other hand, at higher shear-rates the contacts become dominant, triggering the typical shear-thickening behavior, with a growing slope as the volume fraction is increased, eventually leading to a quasi-discontinuous shear-thickening at $ϕ=0.55$⁠. When considering dense suspensions, to characterize the strain-induced anisotropy of the fluid, two additional rheologically meaningful functions are usually measured: the two independent normal stress differences N₁ and N₂. Fig. 8 shows the first (middle panel) and the second (right panel) normal stress differences, respectively, of the shear-dependent suspensions considered (see Table III). In particular, the trend observed is similar to the one reported by Mari et al.,²⁶ proving the capability of the software to catch the strain-induced anisotropies of the suspensions.

Next, we consider a dense suspension with volume fraction $ϕ=0.46$⁠, and we compare our numerical results with experimental data obtained from quasi-neutrally buoyant silica colloidal particles, processed through the Stöber method,^35,36 suspended in a glycerol–water solution with relative concentration $80%−20%$⁠, respectively. The electro-chemical properties of the suspension were altered by adding sodium-chloride (NaCl) to the solution, as done in Rathee et al.³⁷ In particular, the numerical and experimental parameters are reported in Table IV. Figure 9 shows a very good agreement between the relative viscosity, η_r, obtained from the numerical and experimental data, for both salt concentrations considered. For the experimental results, the shear-rate $γ̇$ has been normalized to match the domain span of the numerical data. The two suspensions behave very differently at lower shear-rates, i.e., in the shear-thinning region. In that regime, the electro-chemical potential acting on the suspension dominates its mechanics. Increasing the salt concentration, the electro-chemical properties of the suspension are modified, with the contribution due to the van der Waals attractive forces becoming more important (lower $FR/FA$⁠). This induces the formation of large clusters that highly contribute to the higher resistance of the suspension to flow,³⁷ and this is reflected in Fig. 9, with the red curve showing much higher values of η_r. Conversely, at higher shear-rates, where the contacts between particles dominate, the behavior of the suspension does not depend on the salt concentration; therefore, the two suspensions have similar values of relative viscosities (see Fig. 9).³⁷ 

The typical size of the cluster of particles for the suspension with higher salt concentration is shown in Fig. 10. In particular, the two panels show the size distribution of the clusters at the extrema of the curve in Fig. 9, i.e., at shear-rate $γ̇/γ̇0=[10−5,100]$⁠, when the suspension reaches its steady-state regime. From Fig. 10, we can see that the distributions show two peaks, the leftmost one concerning small clusters and the rightmost large ones. Comparing the distributions of the latter at the two shear-rates considered, we can see that the peak is narrower in the region where the van der Waals forces dominate the behavior of the suspension (left panel), with the particles gathering in a spherical-shaped cluster. A snapshot of the structures for the case with a salt concentration of $0.4M$ is shown in Fig. 11 (left panel), where the two peaks seen in the distributions can be clearly recognized. In the right panel of Fig. 11, instead, a snapshot during the time evolution of the clusters is shown, with the large cluster fragmented in many smaller ones.

The examples discussed above show that the code is able to tackle particles suspended in a uniform shear-flow. We extend here the simulations to non-uniform shear-flows that activate the Faxén contribution introduced in the Stokes drag. In particular, we introduce a wavy shear-flow,

where $U$ is the amplitude of the wave and n imposes the number of spatial periods within the domain. This is a condition of relevant interest as it can be used for example to reproduce the wall effects that appear in a rheometer when attempting to replicate a uniform shear-flows or to take into account the back-reaction of the particles on the flow itself.

We consider a cubic box containing a bidispersed suspension with volume fraction $ϕ=0.5$⁠. The ratio between the radii of the two dispersed phases is $a2/a1=3.0$⁠, while the volume percentage ratio of the largest dispersed particles to the smallest one is $V2/V1=0.25$⁠. We set the number of particles to N = 65 536 to test the memory-efficiency of the software. The side length of the box accommodates the volume fraction $ϕ=0.5$⁠, being $a0=a1=1$⁠. The particles composing the suspension are considered to be rigid with a non-smooth surface (i.e., with a friction coefficient $μC=1$⁠), electrically charged with a screening length $λD=0.02a0$⁠. The equivalent shear-rate, Eq. (14), is set to $γ̇̂=1$ (i.e., within the shear-thickening region of Fig. 9), with a ratio $FR/FA=19$ when the pair of particles (i, j) are at contact. The wave was chosen to have a single spatial period, n = 1, with a large amplitude $U=20γ̇a0$ to amplify the effects of the wavy shear. Table V wraps up the set of parameters chosen and a snapshot of the domain and particles is shown in Fig. 12.

A wave with a large amplitude, such as the one applied in the case described above, significantly affects the local shear-rate, as shown in the left panel of Fig. 13. This local effect has a direct impact on the suspension as half of it is located in a low shear-rate region, while the other half feels a high shear-rate. Thus, the non-uniformity impacts the rheology of the suspension modifying the relative viscosity η_r of the suspension, as shown in the right of Fig. 13. Note that both cases analyzed started from an identical initial condition.

Next, in Fig. 14 we plot the shear-wise diffusivity (left panel) and the autocorrelation function of the shear-wise velocity fluctuations (right panel), respectively, defined as

and

where in both equations, the $⟨ ⟩$ brackets indicate the average over all the particles and t₀ is the initial time; y(t) and $v′(t)$ are the time-dependent shear-wise location and velocity fluctuations, respectively, of the particles. The figure shows how the particles, on average, diffuse more (twice more) when the wavy-shear rate is applied to the suspension, due to the non-uniformity of the shear rate. Moreover, the autocorrelation function of $v′$ shows a higher integral timescale for the suspension with the non-uniform shear-rate.

Finally, to understand the temporal behavior of the suspension, we show the distribution of the local volume fraction along the shear-wise direction. Figure 15 shows the distribution of the particles in the shear-wise direction separated by the size of the suspended phase. In particular the figure is organized as follows: in the left panel, the local distributions of the large (solid line) and the small (dotted line) particles of the suspension at the initial time are shown; in the middle panel, the same quantities are shown for the particles suspended in a uniform shear-flow at γ = 10; in the right panel, the local distributions are shown for the particles suspended in a wavy shear-flow at γ = 10. While the distribution of the particles in the case where a uniform shear-rate is applied to the suspension remains similar to the initial condition with the particles uniformly distributed along the shear-wise direction, in the case where the particles are suspended in a non-uniform shear flow, two clear peaks in the distribution of the large particles appear. This suggests that the large particles are migrating toward such region that, in particular, coincide with the location of zero-shear.

We have presented an efficient software (freely available at https://github.com/marco-rosti/CFF-Ball-0x) that tackles the dynamics of dense suspensions. First, the models of the contributions adopted to simulate the particle–particle and the particle–flow interactions have been presented. The models mimic the ones already available in the literature, echoing in particular, the implementation of Mari et al.²⁶ and Ge et al.^29,30 The core of the algorithm, i.e., the fixed-radius near neighbors search, follows the idea of Hoetzlein²⁵ with an improved symmetric search that guarantees a speed-up in the near neighbors inquiry and a good load-balance among processes when the parallelization is adopted. Particular attention has been also given to the memory efficiency aspect of the implementation as the representation of rheological suspension with a large number of non-smooth particles is of high interest, especially to cope with particle suspensions having a high size ratio, as a limited number of large particles would fit into small domains.

The aforementioned features of the implementation have been proved in the examples carried out in Sec. III. First, we verified the reliability of the method by comparing our results with data available in the literature obtained from other codes or experiments. Then, we carried out sample simulations with a high number of particles $N=216$⁠, comparing the rheology of suspensions driven by a uniform and a wavy shear-flow.

The present code can be extended in a series of ways. Models in consideration for future implementations concern the possibility of reproducing Brownian particles and the introduction of walls that constrains the particles in any of the three directions. The present code is meant to be used for studying the rheology of particle suspensions in several configurations, to unravel some of the complex phenomena arising when dealing with dense suspensions.³⁸ 

All authors gratefully acknowledge the support of Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan. The authors also acknowledge the computer time provided by the Scientific Computing section of Research Support Division at OIST. A. Q. S. acknowledges funding from the Joint Research Projects (JRPs) supported by JSPS and SNSF. V. R. also acknowledges financial support from the Japanese Society for the Promotion of Science (Grant No. 21K13895).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

The hydrodynamics of a dense suspension is dominated by contribution of the close approaching particles.¹⁵ In a Stokesian framework, this translates in dropping the full description of the flow and modeling the hydrodynamics contribution acting on the suspension with only close-range terms, described by the resistance matrices in Eq. (3).

Following the description in Mari et al.²⁶ and adding the Faxén laws,²⁷ the matrices for rigid spheres can be written as