The microstructural change causing the failure of the Cox-Merz rule in Newtonian suspensions: experiments and simulations

Abstract

Newtonian non-Brownian concentrated suspensions show a mismatch between the steady state and the complex viscosity, whatever the strain amplitude imposed in the oscillatory flow. This result is counterintuitive in the two extreme cases of vanishing strain amplitude and very large one. In the first case, the oscillatory flow should not be able to alter the steady microstructure, as well as in the other opposite limit for which the strain amplitude is so high that the oscillatory flow resembles a steady flow reversal. If the microstructure is not altered with respect to the steady one, similarly the complex viscosity should be equal to the steady one. We here investigate experimentally and numerically the origin of the viscosities mismatch at any imposed strain amplitude. We focus on the first two or three cycles of oscillations and different particle concentrations. Experimental and numerical results agree and allow to prove that for intermediate amplitudes, the oscillatory shear induces the breakage of particle clusters and the microstructure modifies so to minimise particle collisions. For very small strain amplitudes, the oscillatory shear only induces the rotation of few couples of touching particles and the complex viscosity results slightly smaller than the steady one, while for very large strains, the oscillatory flow reshuffles the particles inducing a microstructure as clustered as the steady state one but with a different angular distribution function. We show that the vast majority of the microstructure rearrangement takes place in the first half cycle of oscillation.

Introduction

Industrial applications involving multiphase fluids are frequently encountered; often a characterization of the rheological properties of these fluids is required to design and control the process. However, this task is not always easy, as, e.g., it is often difficult to measure with rotational rheometers the steady-state viscosity at the high shear rates of industrial interest because of experimental artefacts as sample or edge fracture, secondary flows and so on. Similarly, fast in line viscosity measurements are highly desired when designing a control algorithm. Typically, the above-mentioned requirements are accomplished with dynamic experiments that are much faster than the steady state ones and that allow to safely rich very high frequencies. From the results obtained with the dynamic experiments, one can get information on the steady-state properties of the fluid under investigation in view of the renowned Cox-Merz rule. In 1958, Cox and Merz (1958) related the steady shear viscosity, η, and the modulus of the complex viscosity, $\left|{\eta }^{*}\right|$ (from here on referred to as ${\eta }^{*}$ for the sake of simplicity), through an empirical correlation stating their numeric equivalence, Eq. (1), if the shear rate, $\dot{\gamma }$, is expressed in s⁻¹, the angular frequency of the oscillatory shear, ω, in rad·s⁻¹ and the oscillatory experiments are executed in the linear regime.

$$\eta (\dot{\gamma })=\eta^*(\omega )$$

(1)

Much more recently, Winter (2009) proposed a modification of Eq. (1) such that the dynamic viscosity was considered as a function of the shear stress, while the complex viscosity of the complex modulus. The rule proposed by Winter provides additional insight into the sample viscoelasticity and may allow for a better comparison among different polymers. The practical impact of Cox-Merz rule is so important that several equations correlating the fluid elastic parameters in steady and oscillatory shear can be found in the literature, but none of them shows the same wide validity of Cox Merz rule. Cox and Merz themselves proposed a correlation for the elastic viscosity, the so-called Cox–Merz rule of consistency; however, this rule was not consistent with Eq. (1) (Booij et al. 1983), though experimentally it proved to work properly for a variety of polymers. Laun (1986) proposed an empirical rule relating the first normal stress difference to the storage and the loss modulus and, more recently, Sharma and McKinley (2012) introduced a more sophisticated empirical rule to predict the first normal stress difference either from dynamic or from steady viscosity data.

The Cox-Merz rule, Eq. (1), originally proposed for flexible polymers, finds a theoretical explanation for linear polymers whose dynamics can be described with the tube model corrected with non-linear mechanisms like the convective constraint release, as first proved by Marrucci (1996) and Ianniruberto and Marrucci (1996), and successively by Mead (2011) in a more consistent way. Snijkers and Vlassopoulos (2014) assessed the validity of the Cox-Merz rule also for polymers with different architectures, though a general theoretical interpretation for flexible polymers different from linear ones is still missing. The Cox-Merz rule has become such a powerful tool in rheology that its validity has been investigated also for complex fluids that exhibit deformation-dependent microstructure like, e.g., heterogeneous fluids and associating polymers. In these cases, Cox-Merz rule typically fails (Al-Hadithi et al. 1992; Larson 1999).

Despite this, though the suspensions are heterogeneous fluids, several attempts to find the range of validity of the rule for this class of materials, and eventually a valid extension of it, were done in the literature. The stress response to both small amplitude oscillatory shear (SAOS), i.e. in the limit of linear viscoelasticity, and large amplitude oscillatory shear (LAOS), necessary to test suspensions with a yield stress, was investigated. Following a chronological order, which leads us to discuss systems with a decreasing level of complexity, the response to LAOS attracted first the literature attention as the suspensions may often exhibit a yield stress. Doraiswamy et al. (1991) derived a new rule, the so-called Rutger-Delaware rule valid for shear thinning suspensions with yield stress, stating the correspondence of the plot of the complex viscosity vs. the maximum (or effective) shear rate (γ₀ ω) with that of the steady-state viscosity vs. shear rate:

$$\eta \left(\dot{\gamma }\right)=\eta^*\left({\gamma }_{0}\omega \right),$$

(2)

where γ₀ is the strain amplitude of the oscillatory shear. The Rutger-Delaware rule was validated experimentally by Mas and Magnin (1997) concluding that Eq. (2) can be used only for large strains, such that the stress resulting in the sample is larger than the yield value. Later on, Wang and Wang (1999) studied shear thinning fluids with yield stress and a measurable Newtonian plateau. In this case, Eq. (2) resulted valid in the shear thinning region, while failed in the Newtonian limit. Raghavan and Khan (1997) showed its validity for large strains also for shear thickening suspensions, both in the shear thinning and in the shear thickening region.

In the linear regime, starting from the observation that shear thinning suspensions under SAOS conditions obey Cox-Merz rule in the Newtonian plateau, while show ${\eta }^{*}>\eta$ in the shear thinning region, Gleissle and Hochstein (2003) derived a modified Cox-Merz rule, valid for this class of materials, applying the shear stress equivalent shear rate concept. They showed that a single parameter B, function of the solid volume fraction ϕ, is required to superimpose the shear stress vs. shear rate plot to that of the complex modulus vs. the angular frequency. The shear stress equivalent shear rate concept accounts for a horizontal shift of the stress curves, through the shift coefficient B, and a vertical shift of the complex moduli, again through B. In terms of viscosity, it states the necessity of a single horizontal shift of the dynamic viscosity plot (Eq. (3)).

$$\frac{\eta (B(\phi )\dot{\gamma })}{B(\phi )}=\frac{\eta^*(\omega )}{B(\phi )}.$$

(3)

Xu et al. (2008) showed that in a shear thinning suspension for particle volume fractions below a critical value, ϕ_c = 0.036, which corresponds to a percolation threshold, Cox-Merz rule held true, while above ϕ_c, a modified rule, actually coincident with Eq. (3), resulted verified.

Finally, Newtonian suspensions attracted less attention, although the modified rules reported in Eqs. (2) and (3) are both ineffective for such materials. Bayram et al. (1998) showed the validity of the Cox-Merz rule for Newtonian suspensions up to a particle volume fraction of 30%, while Kaully et al. (2007) showed its validity up to a particle volume fraction of 20–40%, depending on the particle size distribution. Both papers showed that for higher concentrations, the complex viscosity resulted smaller than the dynamic ones (${\eta }^{*}<\eta$), in contrast to the case of shear thinning suspensions.

To try to explain why Cox-Merz rule fails also for Newtonian suspensions under SAOS conditions, let us recall that non-Brownian concentrated suspensions may show a peculiar dynamic even when subjected to a linear flow like that of Stokes. In simple unidirectional shear, since the papers of Eckstein et al. (1977) and Leighton and Acrivos (1987), it is known that in the creeping flow limit, the particles of a non-Brownian suspension may diffuse because of hydrodynamics and interparticle interactions. In the same limit, under small amplitude oscillatory shear, Gadala-Maria and Acrivos (1980) observed particle re-suspension, driven by the shear-induced particle migration, and, later on, Pine et al. (2005) demonstrated experimentally the particle trajectory irreversibility if the strain amplitude overcomes a critical value, function of the concentration. Both phenomena are due to the suspension microstructure reorganization under oscillatory shear, driven by a diffusion-like mechanism that is suggested to be due to the chaotic nature of the hydrodynamic interactions (Drazer et al. 2002). More recently, based on a simplified numerical model, it was suggested that the origin of irreversibility and chaos in oscillating suspension are the contacts between particles (Metzger and Butler 2010; Metzger et al. 2013). Breedveld et al. (2001) showed that the particle self-diffusion revealed two different regimes, a long-term and a short-term one.

As regard the long-term regime, Bricker and Butler (2006) studied the rheological behaviour of concentrated suspensions in a long-term oscillatory shear flow as a function of the accumulated strain (γ_tot = 4 γ₀ n_c, with n_c the number of cycles, of the order of thousands). They showed that the regime complex viscosity, reached at large γ_tot, decreases with the applied strain amplitude, it passes through a minimum in correspondence of γ₀ ≈ 1 and then it increases again. At small amplitudes, the regime complex viscosity may also overcome the steady one (Martone et al. 2020). This behaviour was explained (Bricker and Butler 2007) with simulations based on the Stokesian dynamics (Brady and Bossis 1988) showing that at small strain amplitudes, the system self-organises in a crystal-like structure, at larger strain amplitudes, in correspondence of the minimum of the regime complex viscosity, a particle ordering into layers orthogonal to the velocity gradient direction is found, while at the largest applied strains, the formation of hydroclusters oriented along the compression direction of the shear flow field, which justifies the viscosity increase, is observed. Successively, in agreement with Breedveld et al. (2001), Lin et al. (2013) suggested that at the largest strain amplitudes, the final microstructure is immediately formed after practically a single oscillation, similarly to what happens in a flow reversal, and its dynamics is not as slow as those of the microstructural reorganizations at the smaller strains that are driven by particle diffusion. Recently, Martone et al. (2020) showed that, in the long-term regime, Newtonian suspensions may even exhibit a dependence on the frequency in oscillatory shear that can be predicted if a weak Van der Waals adhesive/repulsive force is considered in the equations of motion (Ge et al. 2020). This weak force, negligible in steady shear, is able to mediate the microstructure evolution in oscillatory shear, which will consequently depend on the imposed frequency.

As regard the short-term regime, Corté et al. (2008) and Park et al. (2011) argued that the microstructural change is due to the collisions between the particles that, at strains below a critical value, self-organise to avoid further collisions so to enter into an absorbing state where a reversible quasi-steady state is attained, while above the critical strain, particle collisions never cease and an irreversible quasi-steady regime is reached. Let us notice that in this regime the frequency dependence observed by Martone et al. (2020) has not enough time to show its effects. The short-term regime microstructural reorganization mechanisms are active only when the strain amplitude overcomes a small minimum critical value required to induce particle collisions.

The short-term and long-term mechanisms cannot be invoked to explain the failure of the Cox-Merz rule as they both require a large number of cycles to show their effects. In particular, the short-term mechanism requires few hundreds of cycles to reach a quasi-steady state (Lin et al. 2013); moreover, the minimum critical strain amplitude to be applied to induce the collisions between the particles decreases with the suspension particle volume fraction and, e.g., it is larger than 10% for ϕ = 0.4. In the SAOS tests usually run to obtain the complex viscosity frequency spectra, experimental points are collected within two or three cycles for each frequency, and the applied strain amplitude is typically much smaller than 10%.

The goal of the paper is to highlight the microstructure reorganization mechanisms that take place in a Newtonian suspension in no-more than three oscillation cycles. Particular relevant are the cases at both vanishing imposed strain amplitudes and very large ones. In the first case, one would expect that the shear should not be able to perturb the microstructure, while in the second case, the amplitude is so large that the oscillation resembles a steady flow reversal (Breedveld et al. 2001; Lin et al. 2013) and then, once again, the microstructure should not result altered by the oscillatory shear. As a consequence, in both cases, the steady viscosity should match the complex one, i.e. lato sensu Cox-Merz rule should hold. To this end, we experimentally study non-Brownian suspensions made of hollow rigid glass microspheres suspended in a Newtonian oil, spanning a wide range of volume fractions under small amplitude and large amplitude oscillatory shear. We then numerically investigate the behaviour of suspensions at three different concentrations made of mono-sized rigid spheres immersed in a Newtonian fluid using the Stokesian dynamics as implemented by Rexha and Minale (2011). We eventually individuate the mechanisms at the basis of the observed discrepancy between the complex and the dynamic viscosity of Newtonian suspensions, also in the two extreme cases of very small and large strain amplitudes.

Experiments

Materials and methods

The suspensions are made of inert glass hollow microspheres dispersed in Polyisobutene. The microspheres, generously supplied by Potters Industries (Sphericel® 110P8), have the volumetric distribution function shown in Fig. 1 with an average diameter of 15.4 μm. The mean density is equal to 1.1 g/cm³. Due to their dimensions, the microspheres can be considered non-Brownian, the Péclet number is indeed always much greater than one, and in the experimental conditions of interest for this paper, it is always larger than 10⁵.

Fig. 1

Particle distribution function measured with the Mastersizer 3000 (Malvern)

The Polyisobutene was provided by INEOS (Indopol® — grade H50). The fluid is Newtonian with viscosity equal to 15.8 Pa.s, at 21.2 °C. Its density is equal to 0.884 g/cm³. The very similar density value of the matrix and the suspended microspheres allows neglecting sedimentation phenomena during the experiments, as the Shields number is on average larger than 10³ (Martone et al. 2020), and as confirmed by the Newtonian behaviour of the suspension (Fig. 2a).

Fig. 2

Steady state and dynamic data for suspensions at different particle volume fraction: (a) viscosity and first normal stress difference (inset) vs. shear rate; (b) complex viscosity, storage (crosses) and loss (circles) moduli vs. angular frequency. The dynamic data are calculated by the rheometer SW at each angular frequency from a single cycle, either the second or the third, depending on the frequency and the quality of the signal. (Data taken with ARES G2 rheometer)

The investigated suspensions have a volume fraction of microspheres ranging from 2 to 45%. Each sample was prepared by hand mixing the desired amount of the two components up to the achievement of a macroscopically homogeneous “cream”, and the sample was then put under vacuum over night to remove the air bubbles entrapped during the mixing process.

Rheological tests were run with two rotational rheometers. A stress-controlled rheometer (NOVA, Reologica) equipped with a stainless steel cone and plate geometry with radius of 12.5 mm and cone angle of 4°. The cone tip truncation is such that the nominal gap is 150 μm and thus the average gap in the sample is 512 μm. The temperature was set to 21.2 °C and controlled with a forced convection oven (ETC) with a precision of ± 0.1 °C. The gap between the cone and plate is much larger than the particle diameter everywhere, also at the tip where it is about ten times the mean diameter; thus, the effects of particle confinements can be neglected. A strain-controlled rheometer (ARES G2, TA) was used with a cone and plate geometry with radius 12.5 mm, cone angle of 6.36° and truncation of 60 μm. In this case, the gap at the tip is about four times the mean particle diameter and the average gap is 727 μm; thus, also with this device, the confinement effects can be neglected. The temperature was set at 21.2 °C and controlled with a Peltier cell with a precision of ± 0.1 °C.

On each sample, both the steady-state viscosity and the complex viscosity were measured. The latter was measured in the linear regime with SAOS and outside it with LAOS. In agreement with the pertinent literature (Bricker and Butler 2006; Lin et al. 2013; Park et al. 2011), also outside the linear regime, a single synthetic parameter, the complex viscosity, was used to study the data so to have a unique whole picture of the suspension behaviour both in SAOS and LAOS conditions. The complex viscosity is typically calculated by the rheometer software from the main harmonic of the oscillatory response; it should be referred to as apparent complex viscosity in LAOS and SAOS conditions, if the response is not perfectly linear. For the sake of simplicity, we will refer to it as complex viscosity. The complex viscosity and the dynamic moduli are calculated in this paper within no more than three cycles, typically from a single cycle that is either the second or the third, depending on the frequency and the quality of the signal.

With preliminary tests, run with a plate-plate device at different imposed gaps, we checked that wall slip was negligible for all samples (Barnes 1995). Only at the highest volumetric fraction, 45%, wall slip sometimes starts to occur and the use of rough geometries could have been preferred (Carotenuto and Minale 2013; Carotenuto et al. 2015; Paduano et al. 2019). However, to avoid a change of the measuring geometry for the different investigated suspensions inevitably complicating their comparison, we preferred the use of smooth cones and plates for all cases and the data at 45% are retained only when wall slip was considered negligible.

Experimental results

The rheology of suspensions is sensitive to the microstructure and so it may depend on the initial conditions. To erase the flow history, we preconditioned each sample with a pre-shear at 1 s⁻¹ for 2 min. All the investigated suspensions essentially show a viscosity constant with the applied tangential stress, or shear rate, and a measurable negative first normal stress difference (N1) linear in the applied shear rate (Fig. 2a). Similarly, in the linear regime, the complex viscosity resulted constant with the applied angular frequency, the loss modulus, G”, linear with it and the storage one, G’, hardly measurable, but quadratic in the angular frequency (Fig. 2b). The occurrence of negative N1, typical of this kind of suspensions (Denn and Morris 2014; Tanner 2018), together with a positive G’ implies the failure of any empirical rule developed to predict N1 from dynamic data.

We checked the linearity regime of each sample by running stress sweep and strain sweep experiments at a constant frequency and we found that the higher the particle volume fraction, the smaller the critical strain that marks the transition from the linear to the non-linear regime. Figure 3 highlights that for volume fractions smaller than 20%, the linear regime extends up to the largest investigated strain. The non-monotonous trend shown in Fig. 3 for the 30%, 35% and 40% suspensions was already observed by Lin et al. (2013). It can be hypothesised that the oscillatory flow alters the initial microstructure, formed during the preconditioning steady state, proportionally to the applied strain amplitude until a limit value above which the strain is so high that the oscillatory shear may resemble a flow reversal after which the initial microstructure is reconstituted. For these large strains, the microstructure is thus similar to that of the smallest applied strains and so the complex viscosity increases back. In Fig. 3, the value of the viscosities of the steady shear preceding the oscillatory test is shown as horizontal segments and it is clear that the Cox-Merz rule holds until a volume fraction of about 30%, if the complex viscosity is measured well within the linear regime, or for very large strain amplitudes. For the two more concentrated suspensions, data of Fig. 3 show the mismatch between the complex viscosity and steady one, i.e. the failure of the Cox-Merz rule, both for the smallest strain amplitudes, 0.1%, and the largest ones (500%). To check the effect of the strain accumulated at each strain amplitude, the data were measured also after a steady preshear run before each single point and the results obtained were superimposable to those shown in Fig. 3.

Fig. 3

Strain sweep tests executed with the ARES G2 rheometer at ω = 1 rad/s on suspensions with different particle volume fractions. Data are calculated by the rheometer SW at each strain amplitude from the second cycle. Horizontal segments indicate the value of the viscosity measured in the pre-shear at 1 s⁻¹

In the remaining part of the paper, we will only show graphs of relative viscosities, i.e. of the ratio between the viscosity of the suspension and that of the suspending fluid. Figure 4 shows the measured relative steady state and relative complex viscosity as a function of the suspension volume fraction and the error bars are calculated as ± the standard deviation of 3 independent measurements, i.e. from data obtained from three independent loadings. The complex viscosity is obtained with an applied strain of 1%, i.e. in the linear regime for all the measured suspensions, and with an applied strain of 100%, which lays outside the linear regime for volume fractions larger than 20% (see Fig. 3), and it corresponds more or less to the minimum values of each curve of Fig. 3. Concerning the behaviour at the highest strains, let us consider that the values of complex viscosity measured at a strain amplitude of 1% also represent the largest strain data, see Fig. 3, and thus to avoid to excessively crowd Fig. 4, only the two mentioned strain amplitudes are considered. As expected (Raghavan and Khan 1997; Stickel and Powell 2005), ${\eta }_{r}$ and ${\eta }_{r}^{*}$ increase by increasing the particle volume fraction. They are almost coincident for small concentrations (ϕ < 0.2), while for larger values, they start to differ, obtaining always ${\eta }_{r}>{\eta }_{r}^{*}$ and ${\left.{\eta }_{r}^{*}\right|}_{{\gamma }_{0}=0.01}>{\left.{\eta }_{r}^{*}\right|}_{{\gamma }_{0}=1}$. Both inequalities increase by increasing the particle volume fraction. We interpolated the steady-state viscosity data with Krieger and Dougherty (1959) equation (Eq. (4)) obtaining ϕ_M = 0.671 ± 0.018 and a = 3.072 ± 0.059 with R² = 0.999. These parameters are consistent with the literature.

Fig. 4

Relative steady state (full symbols) and relative complex (empty symbols 1% strain, and black and white symbols 100% strain) viscosity versus volumetric concentration, at 21.2 °C. The line is Krieger and Dougherty interpolation (Eq. (4)) of the steady-state values. (Data taken with NOVA rheometer)

$${\eta }_{r}\left(\phi \right)={\left(1-\frac{\phi }{{\phi }_{M}}\right)}^{-a {\phi }_{M}}$$

(4)

The relative viscosity difference, normalised with the value of the relative steady-state viscosity, $\left({\eta }_{r}-{\eta }_{r}^{*}\right)/{\eta }_{r}$, henceforth referred to as normalised delta viscosity, is plotted vs. the particle volume fraction in Fig. 5. For volume fractions smaller than 20%, both normalised delta viscosities are smaller than 1%, and this value is of the order of the experimental error and thus the three viscosities can be considered coincident; consequently, in the graph in Fig. 5, we do not show these points. For ϕ > 20%, the normalised delta viscosity increases monotonically with the particle volume fraction reaching values larger than 30% at the highest concentrations.

Fig. 5

Normalised delta viscosity as a function of the volumetric fraction of microspheres at 21.2 °C. Lines serve to guide the eye

Numerical simulations

To individuate the possible mechanisms that lead to a complex viscosity smaller than the dynamic one, both in the linear and non-linear regime, we numerically investigated a suspension of non-Brownian spheres immersed in a Newtonian fluid undergoing simple shear flows, steady or oscillatory. The numerical approach is based on the Stokesian dynamics (Brady and Bossis 1988) with the Eulerian/Lagrangian scheme proposed by Rexha and Minale (2011). We first summarise the modified Stokesian dynamics and then present the numerical results.

Stokesian dynamics

The motion of the fluid is governed by the Stokes equations, and the particle motion is described by the coupled N-body Langevin equations, which can be written as follows

$${\bf m}\frac{d{\bf U}}{dt}={\bf F}^{H}+{\bf F}^{P}+{\bf F}^{B}$$

(5)

where m is the generalised mass/moment of inertia matrix of the suspended particles, U is the particles translational/rotational velocity vector, d/dt denotes the time derivative, F^H is the hydrodynamic forces vector acting on the particles due to their motion relative to the fluid, F^P is the deterministic non-hydrodynamic forces (which may be either interparticle or external) vector and F^B is the stochastic forces vector that gives rise to the Brownian particles motion; here it is set to zero since the suspensions of interest are non-Brownian. Particle inertia, i.e. the LHS of Eq. (5), is also set to zero being the flow Stokesian.

The following analytical expression, Eq. (6), for the hydrodynamic forces can be obtained if (i) the suspension undergoes a simple shear flow; (ii) the particle Reynolds number is small:

$${\bf F}^{H}=-{\bf R}_{FU}\cdot ({\bf U}-{\bf U}^{\infty })+{\bf R}_{FE}:{\bf E}^{\infty }$$

(6)

where U^∞ is the imposed unperturbed velocity, E^∞ is the unperturbed rate of deformation tensor and R_FU and R_FE are the configuration-dependent resistance matrices that give the hydrodynamic force/torque on the particles due to their motion relative to the fluid (R_FU) and due to the imposed shear flow (R_FE). The interested reader can find all the information required to compute both matrices in the literature (Arp and Mason 1977; Durlofsky et al. 2006; Jeffrey and Onishi 1984; Kim and Mifflin 1985).

The deterministic, nonhydrodynamic force vector F^P is arbitrary and may have almost any form. We here consider only short-range repulsive inter-particle forces required to prevent particle overlapping, as suggested by Dratler and Schowalter (2006), Eq. (7).

$${\bf F}_{\alpha \beta }^{*p}=\frac{\tau }{\dot{\gamma }}\frac{{e}^{-\tau h}}{1+{e}^{-\tau h}}\frac{{\varvec{r}}}{\left|{\varvec{r}}\right|}$$

(7)

where ${\mathbf{F}}_{\alpha \beta }^{*p}$ is the dimensionless force on particle α exerted by particle β, r is the centre-to-centre dimensionless separation vector directed from particle β to particle α, τ controls the spatial decay rate of ${\mathbf{F}}_{\alpha \beta }^{*p}$, $\dot{\gamma }$ is the nondimensional shear rate, such that the ratio $\tau /\dot{\gamma }$ tunes the relative importance of inter-particles forces with respect to shear forces and h is the dimensionless distance between the particles surfaces, $h=\left|{\varvec{r}}\right|-2$. It should be noted that the use of F^P is required only for the numerical stability of the code; thus, its precise form is unimportant provided that it decays rapidly to zero while increasing the particles separation. We, then, checked that in our calculations, (i) the contribution to the viscosity given by the inter-particle forces is always some orders of magnitude smaller than that of the hydrodynamic forces and (ii) that particle overlaps are absent. Once both the hydrodynamic and the inter-particles forces are known, the relative viscosity of the suspension can be computed (Brady and Bossis 1985). Let us underline that, being F^P negligible with respect to the other forces, the set of Eqs. (5) and (6) is linear in the imposed shear rate, or angular frequency in the case of oscillatory shear; thus, the corresponding dimensionless set of equations is rate independent. As a consequence, all the material functions are also predicted to be rate independent, in agreement with what observed in Fig. 2.

Due to the long computational time required to simulate an infinite 3D suspension, we focused on a monolayer of uniform rigid spherical particles, since it was shown that the bulk properties of a fully three-dimensional suspension are well caught by the monolayer, provided that the areal fraction of the spheres is opportunely scaled to obtain the corresponding volume fraction of the real 3D suspension (Brady and Bossis 1988). Moreover, to simulate a non-confined suspension, cell periodicity is introduced and we implemented a Eulerian/Lagrangian approach (Rexha and Minale 2011) following which we calculated the motion of a set of particles contained in a Eulerian box, but we evaluated the forces acting on each particle referring to a Lagrangian box centred in the particle itself. In Fig. 6, the computational Eulerian box, shadowed, is plotted together with its eight images, the test particle is shown as a hollow circle and the Lagrangian box, centred in it, is shown with a grey border. Finally, to reduce the computational time, we also implemented a variable time step algorithm. All details on the numerical scheme and on the values of the parameters required for the numerical integration can be found in Rexha and Minale (2011).

Fig. 6

Cell periodicity: The Eulerian box (shadowed) and the Lagrangian box (with grey border) centred in the test particle, a hollow circle, are shown

Numerical results

We investigated several configurations made of 9 rigid non-Brownian spheres immersed in a Newtonian fluid undergoing simple and oscillatory shear. The adopted areal fractions ϕ_a, and the equivalent volume fractions, are reported in Table 1. The choice of 9 spheres is appropriate since here we only aim at individuating the mechanism underlying the microstructure modification in the very firsts cycles of an oscillatory shear and the qualitative behaviour of the suspension is nicely caught already by a configuration with such a small number of particles, as demonstrated by Rexha and Minale (2011). However, by investigating samples with such a small number of particles, we are conscious that we may have difficulties to obtain a quantitative agreement with the data, but the gain in computational time is strong enough to justify our choice. We focused on volume fractions not smaller than 25% because we showed experimentally, Figs. 2 and 3, that Cox-Merz rule ceases to be valid for ϕ > 20%, and because we already numerically proved the validity of Cox-Merz rule for volume fractions up to 17% for the same kind of suspension (Rexha and Minale 2011).

Table 1 Areal fractions and the corresponding volumetric fractions of the samples investigated in the numerical simulations

The experimental complex viscosity is uniquely determined once the system is opportunely preconditioned with a steady pre-shear, thus, in the numerical simulations we first integrate in the time domain the set of Eqs. (5)–(7) with an imposed constant simple shear flow until steady state and then, starting from a few configurations representative of the steady state, the oscillatory shear is simulated. For each sample, we always checked to have reached a stable steady state. The results of the steady simple shear are shown in Fig. 7 for the 25% suspension, chosen as an example. At time zero, the initial configuration is a random one. The solid grey line represents the instantaneous relative viscosity, the dots are the relative viscosity mean values calculated starting from time 0 to the running time and their associated error bars are estimated as ± the standard deviation of the data. The horizontal line represents the final steady-state relative viscosity. Let us consider that the unique experimental value, determined on a sample with a practically infinite number of particles, can be seen as equivalent to that numerically obtained averaging the results coming from an infinite number of configurations of a finite number of particles. In steady state, this is done by averaging the results coming from a configuration evolving in the time domain, while in the oscillatory shear experiments, this approach cannot be followed. Indeed, each dynamic numerical simulation is based on a single configuration and strongly depends on it. To mimic the reality, it is necessary to average the results coming from a set of different initial configurations that must be chosen so that accurate results can be obtained with a set as small as possible. To properly choose the initial configurations set so that it nicely represents the steady state, we focus on the viscosity associated with each particle configuration and we require that the distribution function of these viscosities has the mean value coincident with that of steady state and, to improve the accuracy of the numerical predictions, we also require that its standard deviation be equal to that of the steady-state data, see Appendix. The stars of Fig. 17 individuate the viscosities of the set of the 9 initial different configurations complying with the discussed requirements, for the 25% sample. The configurations associated with the stars, shown in Fig. 7, are those forming the initial set used for the oscillatory shear simulations.

Fig. 7

Numerical simulations: Relative viscosity vs dimensionless time in steady shear. ϕ = 0.25. The initial configurations used for the numerical simulations in oscillatory shear are shown. Two configurations, A and B, are highlighted (see Fig. 11)

Figure 8 shows the mean complex viscosity, obtained averaging the values originated from the chosen sets of 9 different configurations, as a function of the applied strain amplitude for three different volume fractions. We simulated approximately the first two cycles for each strain amplitude, and we discarded the first half cycle in the calculation of the complex viscosity from the main harmonic of the periodic response. Analogously to the experimental data, it can be seen that the relative complex viscosity of the three investigated suspensions is always smaller than the relative steady state one, whose value is shown with a horizontal segment in the graph; the complex viscosity decreases as the applied strain amplitude increases, up to a strain of about 100%, and then it increases again. In agreement with what observed in Fig. 3, the larger the volume fraction, the more evident the effect of the applied strain on ${\eta }_{r}^{*}$, and the smaller the applied strain amplitude corresponding to the minimum complex viscosity. Let us notice that the numerical data of Fig. 8 qualitatively reproduce the experimental data of Fig. 3, and this confirms that the numerical simulations based on only 9 particles correctly catch the physics of the problem. The mismatch between the steady and the complex viscosity at small and large strain amplitudes is also predicted.

Fig. 8

Numerical simulations: Relative complex viscosity vs applied strain amplitude for suspensions with ϕ = 0.25, 0.30, 0.35. The horizontal segments show the relative viscosity of the preceding steady shear

In Fig. 9, we compare the relative steady state and complex viscosities, for γ₀ = 1% and 100%, as a function of the volume fraction confirming that the numerical simulations show a behaviour in qualitative agreement with the data plotted in Fig. 4: The larger the volume fraction, the larger the difference between the steady-state viscosity and the complex one, at both applied strain amplitudes.

Fig. 9

Numerical simulations: Relative steady-state viscosity and relative complex viscosity, at γ₀ = 0.01 and 1 vs. the suspension volume fraction. Lines serve to guide the eye

Mechanisms at the basis of Cox-Merz rule failure

In Fig. 10, we compare the experimental and simulation results restricting our attention to the volume fractions investigated numerically. We show the normalised delta viscosity as a function of the volume fraction. The numerical results over-predict the data at γ₀₀ = 100%, while under-predict them at γ₀ = 1%, however, this does not surprise since the choice of using a set of 9 particles can only lead, as previously underlined, to a qualitative agreement with experimental data, that is however fully caught by our simulations.

Fig. 10

Normalised delta viscosity, at γ₀ = 1% (circles) and 100% (squares) vs. the suspension volume fraction. Comparison between experiments (hollow symbols) and simulations (filled symbols). Lines serve to guide the eye

To understand the possible physical origin of the observed discrepancy between the complex and the steady-state viscosity, let us remind that the more aggregated the particle microstructure, the higher the viscosity of the suspension. Being the complex viscosity smaller than the dynamic one, we can then hypothesise that the oscillatory shear induces a structure dilation, i.e. a particle declustering, of the suspension; this means that on average after the first two oscillations, the particles are better dispersed in the matrix. This idea is strengthened by the results shown in Fig. 11 where the complex viscosity coming from two different configurations, namely A and B of Fig. 7, is plotted vs. the strain amplitude. The structure of configuration A is very dilated (declustered), being its corresponding steady-state viscosity smaller than the mean one, while that of configuration B is more clustered, being its corresponding steady-state viscosity larger than the mean one (Fig. 7). Starting from the dilated configuration A, the oscillatory shear (in only two cycles) tends to aggregate the particles and indeed the complex viscosity increases with the applied strain amplitude (Fig. 11); the opposite holds when starting from the clustered configuration B, at least until an applied strain amplitude equal to 100%, for larger strain amplitudes, the complex viscosity increases again, thus suggesting a re-clustering effect. Let us also notice that the complex viscosity of configuration B deviates from the initial value since the smallest applied strain amplitudes, while to induce an increase of the complex viscosity of configuration A, a much larger strain amplitude is required, > 20%. As a reference in Fig. 11, the mean complex viscosity is plotted vs. the applied strain amplitude (data of Fig. 8) as grey dots. Though the viscosities of configurations A and B differ from the mean value approximately of the same amount, at each applied strain amplitude, the dilation effect induced by the oscillatory shear is more effective than the clustering one thus leading to an overall reduction of the complex viscosity with the strain amplitude. This suggests that in the first two cycles, the oscillatory shear is much more effective to disaggregate than to aggregate the microstructure until a strain amplitude of 100%. For larger strain amplitudes, the data of configurations A and B are very similar because, as we already discussed, in these conditions, the oscillatory shear resembles a flow reversal whose steady value is independent from the initial conditions.

Fig. 11

Relative complex viscosity as a function of the applied strain for a dilated particle configuration (A) and a clustered one B with φ = 0.25. The horizontal lines correspond to the initial relative viscosities of the two configurations, see Fig. 7. As a reference, the relative suspension viscosity (dotted line) and the mean relative complex viscosity (grey dots, data of Fig. 8) are shown

To confirm that the reduction of complex viscosity is due to a microstructure disaggregation, we evaluated the radial pair distribution function immediately before the oscillatory shear and upon it for selected strain amplitudes and three volume fractions, Fig. 12. The radial pair distribution function, rdf, describes how particle density varies as a function of the distance r from a reference particle and is calculated following the Eulerian/Lagrangian approach implemented in our model. The reference particle is put in the centre of the Lagrangian box, concentric circles with progressively increasing radius r are drawn around it and the number of particles contained in each shell is calculated. This number is then divided by the area of each shell and by the average density of particles in the system. This calculus is repeated considering in turn as reference all the particles contained in the Eulerian box and all the results are then averaged to obtain the rdf of a single configuration. The final rdf is obtained after averaging the rdf of the nine configurations used for the simulations.

Fig. 12

Radial pair distribution functions (rdf) before the oscillatory shear (initial blue line) and upon it (final red line) for five selected applied strain amplitudes and three volume fractions. The horizontal line is for rdf = 1

For 10% ≤ γ₀ ≤ 100%, the radial pair distribution function upon oscillation is significantly smaller than that before the oscillation only at r = 2 for all the three volume fractions. This means that the probability to find a couple of particles touching each other is significantly reduced by the oscillatory shear; thus, the existing clusters generated by the steady shear are disrupted by the oscillatory flow. When the strain amplitude is larger than 100%, the picture changes with the particle concentration: At γ₀ = 200%, for ϕ = 25% and 30%, the radial pair distribution function upon the oscillatory shear is still significantly smaller than the initial one at r = 2; thus, a particle declustering is still visible, while for the suspension at ϕ = 35%, the two rdf’s are very similar at every r and thus the disruption of clusters cannot be appreciated; at γ₀ = 400%, for ϕ = 30%, we still measure a difference between the initial and the final rdf at r = 2, while this is no more visible for ϕ = 25%. What observed for γ₀ ≥ 10% clearly explains the results of Figs. 3 and 8: For intermediate strain amplitudes, the oscillatory shear is capable of inducing a particle declustering that leads to a complex viscosity significantly smaller than the steady one; for larger amplitudes, with a threshold depending on the volume fraction, the oscillatory shear starts to behave like a flow reversal and the final rdf tends to the initial one; consequently, the complex viscosity increases back towards the steady one, though it remains at the highest concentration always smaller than it. Finally, when γ₀ ≤ 1%, the two radial pair distribution functions are practically coincident and so the dilation action is absent. This holds for all the three concentrations and in these conditions, the numerical complex viscosity is similar to the steady one; nevertheless, at the higher concentrations, it remains smaller than it, even if the regime complex viscosity reached at large accumulated strain may be larger than the steady one.

As discussed above, the data of Fig. 8 at γ₀ ≤ 1% and γ₀ = 200% show that for the higher volume fractions, a residual difference between the steady state and the complex viscosity is still measurable. This happens even if in these cases the radial pair distribution functions of the samples before and after the oscillatory shear are practically superimposable. A similar discrepancy between oscillatory and steady data is observable also experimentally. In these conditions, the reduction of complex viscosity cannot be thus attributed to a particle declustering. We then hypothesise that the oscillatory shear at small strain amplitudes and at very large ones, due to its fore-aft asymmetry, may induce a rotation of anisotropic structures that mostly contribute to the suspension viscosity, i.e. clusters of couple of particles will orient to exert a smaller resistance to flow. We then plot in Fig. 13 the angular pair distribution function of the 35% suspension before and after exactly two cycles of oscillatory shear, for γ₀ = 1% and 200%, considering only practically touching particles, i.e. particles whose centres have a distance r ≤ 2.01. In the case of γ₀ = 200%, the two adf’s differ significantly, proving that a macroscopic rotation is induced to the vast majority of clusters of couple of particles. In the case of γ₀ = 1%, the two adf are much more similar and only few, but significant, differences are observed. These small differences may justify the small difference between the complex and the steady viscosity observed in these conditions. Let us underline that when γ₀ = 200%, the adf before and after the oscillatory shear remain different also for non-touching couple of particles, r > 2.01, while for γ₀ = 1%, they result always practically coincident. This confirms our interpretation that when the applied strain amplitudes are large enough, the oscillatory shear behaves similarly to a flow reversal and the final particle configuration results as clustered as the initial one, as revealed by the rdf, but differently oriented with respect to the flow, as shown by the adf. This difference of orientation may justify the observed difference in viscosity. Probably, a much larger strain amplitude is necessary to obtain the adf of the steady state. Conversely, in the case of small strain amplitudes, γ₀ ≤ 1%, the structure is practically unaffected by the oscillatory shear and only few touching particle results rotated by the oscillatory flow. Since this anisotropic structure largely contributes to the suspension stress, even such a small rotation may justify the observed complex viscosity reduction.

Fig. 13

Angular pair distribution functions of particles (r ≤ 2.01) before the oscillatory shear (initial blue line) and after exactly two oscillations (final red line) for the suspension with ϕ = 35%. 0° is the velocity direction of the steady shear flow. The inset defines the angle θ

Further insight into the microstructural reorganization mechanism

Going back to the structure evolution of a concentrated suspension under oscillatory shear flow, we already discussed in the Introduction that two mechanisms can be individuated: a short-term and a long-term one. According to the first (Corté et al. 2008, 2009; Lin et al. 2013), if the strain amplitude is larger than a critical value, the system enters into an absorbing state where particle collisions are suppressed and a reversible pseudo steady state is attained, with η* < η; if the applied strain is larger than a second higher critical value, an irreversible steady state is reached where particle collisions never cease and the complex viscosity is not very different from the steady state one. What observed in this paper, both experimentally and numerically in the first very few cycles of oscillation (< 4) at intermediate strain amplitudes, resembles much to what happens according to the short-term mechanism in the first few hundreds of cycles. To further investigate this, we experimentally run strain sweep hysteresis loops by first increasing and successively decreasing the applied strain amplitude at a fixed angular frequency. The first point of the hysteresis loop is always run after the conditioning pre-shear. The results are shown in Fig. 14 for the suspensions at 30, 35 and 40%. Both the modulus of the relative complex viscosity and the out-of-phase component of the relative complex viscosity, η″_r, are shown. The latter reflects the elasticity of the system that in non-Brownian suspensions is due to either collisions between particles, or interparticle interactions that become relevant when the particles get very close to each other.

Fig. 14

Forth and two back strain sweep experiments at 1 rad/s (a, c); Forth, back and forth strain sweep experiments at 10 rad/s, compared with the data at 1 rad/s (b, d). (Data taken with ARES G2 rheometer)

The first run is performed in the same way in which we obtained the data of Fig. 3, and indeed, at constant angular frequency, by increasing the applied strain amplitude, the relative complex viscosity of all three investigated suspensions first decreases towards a minimum at about γ₀ = 100%, and then it increases again. Contemporaneously, by increasing the applied strain amplitude, first η″_r is not measurable at the smallest strain amplitudes, then from γ₀ = 2–5%, which corresponds to the value at which η*_r starts to decrease, it monotonically increases suggesting that the elasticity due to collisions increases monotonically with the applied strain amplitude. Let us underline that η″_r is always at least one order of magnitude smaller than η*_r, this means that the system is mainly viscous, but few collisions between particles are always present. At a fixed angular frequency, the hysteresis loop is performed by stopping increasing the applied strain amplitude, at any chosen value, and from there, γ₀ is decreased at the same rate of the preceding “forth” branch. This “back” branch is started from two strain amplitudes: 500% and 100%. In the first case, the complex viscosity follows the “forth” branch until its minimum value reached at γ₀ ~ 100%; from there on, it remains more or less constant and a hysteresis shows up. In the second case, when the applied strain amplitude is decreased starting from γ₀ = 100%, the complex viscosity remains practically constant and not significantly different from the value at γ₀ = 100%. This holds at all the volume fractions. η″_r also shows the hysteresis (Fig. 14c) and, similarly to the complex viscosity, the “back” branch of η″_r detaches from the “forth” branch at γ₀ ~ 100%, from where it rapidly drops to a very small value that is in the limit of measurability of the instrument. This suggests that in the “back” branch, for γ₀ ≤ 100%, there are no more collisions between particles during the oscillatory shear and thus the microstructure does not change much anymore. This is in agreement with the numerical predictions of particle declustering induced by the oscillatory shear which led to a final configuration with particles distributed in a such way that they were no longer prone to collide or interact. To fortify this interpretation, we experimentally showed that the “back” branch is reversible; indeed, starting from any point of it, we increased again the strain amplitude and the “back” branch was retraced by both η*_r and η″_r (Fig. 14b, d). We also checked that the angular frequency had no effects on this phenomenon by running the same experiments at 10 rad/s and obtaining results perfectly superimposable to those at 1 rad/s (Fig. 14b, d).

To validate the observations of Fig. 14, we run a hysteresis loop for the 35% suspension also numerically (Fig. 15). The data of the “forth” branch are those of Fig. 8, while those of the “back” branch are obtained using as initial conditions for the simulations the configurations reached after the oscillatory shear at strain amplitude of 100%. Only three strain amplitudes are investigated: 100%, 20% and 5% and the qualitative agreement with the experimental data of η*_r is excellent, though the quantitative agreement is not perfect. This does not surprise since we simulated configurations with only nine particles, and this allows catching the physics of the problem but makes a quantitative comparison difficult. The out-of-phase component of the complex viscosity, Fig. 15b, is rather small and difficult to be predicted; however, it also shows the hysteresis: In the “forth” branch, an increase of η″_r with the strain amplitude is detected, while in the “back” branch, it immediately falls down towards values too small to be accurately numerically predicted. The numerical data clearly show that in the “back” branch, the elasticity, given by the particle interactions, is much smaller than that of the “forth” branch; thus, the particles are oriented by the oscillatory flow so to minimise their interactions.

Fig. 15

Numerical simulations: Relative complex viscosity (a) and out-of-phase component of the complex viscosity (b) vs. the applied strain amplitude for the suspension with ϕ = 0.35. Forth and back strain sweep experiments

What happens in the first cycles is then in agreement with the first steps of the structure evolution dictated by the short-term dynamics, that is however not yet completed as demonstrated in Fig. 14 where the data of the “back” branch still show an evolution of the complex viscosity. Indeed, the constant value attained in “back” branch is always smaller than the minimum found in the “forth” one. These results suggest that in the first two or three cycles, the microstructure evolution is very fast and pronounced leading always to a decrease of viscosity. This holds regardless the strain amplitude and the long-term regime complex viscosity, which can be larger or smaller than the steady-state viscosity. This fast microstructure reorganization is followed by evolutions with much slower dynamics leading to the short-term regime, eventually followed by the long-term one.

Conclusions

We have experimentally investigated the behaviour of suspensions made of rigid non-Brownian spheres immersed in a Newtonian fluid in steady and oscillatory shear flow both at small and large strain amplitudes focusing on the first two or three cycles. We showed that if the solid volume fraction is larger than 20%, the complex viscosity is systematically smaller than the steady state one, and the higher the concentration the larger the difference. The discrepancy between the steady and the complex viscosity increases with the applied strain amplitude at any concentrations until a maximum that is reached for amplitudes of about 100%, and for larger amplitudes, the complex viscosity increases again tending towards the steady one. To explain the mechanism at the basis of these differences, which imply the failure of the Cox-Merz rule, we numerically investigated the same kind of suspension at three different volume fractions both in steady simple shear and in oscillatory shear, focusing on the first two cycles. The qualitative agreement between data and simulations is very good and we numerically evidenced that for strain amplitudes 1% < γ₀ < 200%, the oscillatory flow induces a microstructure declustering. More precisely, if the oscillatory flow starts from a clustered initial configuration, it induces the particle declustering, while it induces the particle clustering, if the initial configuration is more dilated. On average, the particle declustering prevails and this leads to the observed complex viscosity reduction. This is proven by the comparison of the rdf of the microstructure as calculated before and after the oscillatory shear flow. Indeed, when the complex viscosity is significantly smaller than the steady one, a pronounced reduction of the probability to find couple of particles in contact is observed after the oscillatory shear. For strain amplitudes larger than 200%, the complex viscosity, though still smaller than the steady one, increases back. In these cases, the oscillations are large enough to remix the microstructure so cancelling the dilation effects, but not so large to cancel also the angular microstructural rotations. In the final structure, any couple of particles, no matter the distance between the centres, have the axis passing through their centres rotated, thanks to the flow irreversibility met by these systems, with respect to the initial configuration, in such a way that the final viscosity becomes slightly smaller than the initial one. In these conditions, the oscillatory shear behaves similarly to a flow reversal without having, however, enough time to reobtain the microstructure of the steady state.

When the applied strain amplitude is extremely small, γ₀ ≤ 1%, the oscillatory shear essentially does not modify the suspension microstructures, but because of the fore-aft asymmetry of flow, it manages to slightly rotate few couples of touching particles, and this is enough to justify the observed small reduction of the complex viscosity, regardless the values reached at regime by complex viscosity that in these conditions can exceed the steady one.

We showed that at intermediate strain amplitudes, 10% ≤ γ₀ ≤ 100%, in the first cycles, the microstructure evolves towards a configuration where the collisions between particles are unlikely. This agrees with the idea that in the first cycles of oscillation, the first steps of the structure evolution dictated by the short-term dynamics are taken. The structure evolution is very pronounced and fast in these first two cycles and it will be followed by an evolution with a much slower dynamics that will lead the system towards the short-term regime, eventually followed by the long-term one.

In Fig. 16, we show the stress vs. time as numerically calculated for a single characteristic initial particle configuration for γ₀ = 100%. It is clearly shown that the vast majority of the stress amplitude reduction is realised already in the first cycle, or even in the first half cycle, while the following evolution appears much slower. In the first cycle, the fore-aft asymmetry is indeed more evident because the flow is imposed on the microstructure deriving from the steady state.

Fig. 16

Numerical simulations: Stress vs. time in oscillatory shear for a single initial particle configuration. ϕ = 0.35, 0.30 and 0.25 and γ₀ = 100%. Horizontal lines show the amplitude of the best sinusoid fitting the data after the first half cycle

Appendix

In Table 2 in Appendix, the values of the mean viscosities, and their standard deviations, of the set of 5, 7 and 9 initial configurations used for the oscillatory shear simulations are reported for the three suspensions with volume fraction of 25%, 30% and 35%. As a comparison, also the mean viscosity and the standard deviation of the steady state are shown. The three different sets of initial configurations have the same mean viscosity of the steady state, while their standard deviation tends to that of the steady state as the number of the initial configurations increases, coinciding with it only for the set of 9 configurations, at any volume fraction.

Table 2 Mean value and standard deviation of the viscosities of the steady state and of the different sets of initial configurations used for the oscillatory shear simulations

In Fig. 17 in Appendix, we show the normalised delta viscosities, calculated for γ₀ = 1% and 100%, as a function of the volume fraction for the different sets of initial configurations and we compare them with those of the steady state. It is clear that for both applied strains, the agreement between data and simulation improves as the standard deviation of the set of the initial configurations tends to that of the steady state. In addition, the qualitative trend of the data is well caught only by the set with 9 configurations for both applied strains.

This suggests that the set of initial configurations must be chosen so to share with the steady state both the mean value of the viscosity and the standard deviation.

Fig. 17

Effect of the distribution function of the viscosities of set of initial configurations used for the simulations of the oscillatory shear. Normalised delta viscosity, at γ₀ = 1% and 100% vs. the suspension volume fraction. Lines serve to guide the eye