Yield stress fluids and fundamental particle statistics

Abstract

Yield stress in complex fluids is described by resorting to fundamental statistical mechanics for clusters with different particle occupancy numbers. Probability distribution functions are determined for canonical ensembles of volumes displaced at the incipient motion in three representative states (single, double, and multiple occupancies). The statistical average points out an effective solid fraction by which the yield stress behavior is satisfactorily described in a number of aqueous (Si₃N₄, Ca₃(PO₄)₂, ZrO₂, and TiO₂) and non-aqueous (Al₂O₃/decalin and MWCNT/PC) disperse systems. Interestingly, the only two model coefficients (maximum packing fraction and stiffness parameter) turn out to be correlated with the relevant suspension quantities. The latter relates linearly with (Young’s and bulk) mechanical moduli, whereas the former, once represented versus the Hamaker constant of two particles in a medium, returns a good linear extrapolation of the packing fraction for the simple cubic cell, here recovered within a relative error ≈ 1.3%.

---

Yield stress fluids form a particular state of matter,¹ displaying non-linear and novel visco-plasto-elastic flow dynamics upon different boundary conditions. As their name says, they don’t flow until a certain load, the so-called yield stress (or point, τ₀), is applied. This value may be generally interpreted as a shear stress threshold for the breakage of interparticle connectivity.² Furthermore, as it initiates motion in the system, it is connected to mechanical inertia³ and particle settling, i.e. it is a terse summary of buoyancy, dynamic pressure, weight, viscous and yield stress resistances.⁴ For prototype systems such as colloids dispersed in a liquid, yield points sensibly depend on the mechanism by which the solid phase tends to interact or aggregate.^5–8 The macroscopic constitutive equations they obey, such as the Herschel–Bulkley model, were shown to correspond, over a four-decade range of shear rates, to the local rheological response.⁹

From the side of an experimenter, however, unambiguously defining a yield stress may not always be straightforward. It can be affected by the experimental procedure adopted, always considering a measurement or some extrapolation technique with the limit of zero shear. Conversely, unyielded domains may be defined by areas where the shear stress second invariant falls below the yield value, plus some small semi-heuristic constant.¹⁰ In addition, theoretically, the meaning of notions like τ₀ and rheological yielding were questioned to be only qualitative or even to stand for an apparent quantity.¹¹ The dependence they generally show on timescales characteristic of the applied (mechanical) disturbance, also suggested an intimate relationship¹² between yield stress and dispersion thixotropy.¹³ On the other hand, assigning a hydrodynamic or mechanical state below the yield point to a material that is not flowing seems not to be scientifically sound. Experimental values are normally obtained by extrapolation of limited data, whereas careful measurements below the yield point would actually imply that flow takes place.¹⁴

At any rate, the analysis of properly defined τ₀ concepts forms the subject of interesting investigations and is still a powerful tool in many applications, including macromolecular suspensions,¹⁵ gels, colloidal gels and organogels,^16–18 foams, emulsions and soft glassy materials.¹⁹ It allows for effective comparisons between the resistances which fluids initially oppose to the shear perturbation, somehow specifying a measure of the particle aggregation states taking place in a given dispersant. Electrorheological materials, for instance, exhibit a transition from liquid-like to solid-like behaviors, which is often examined by a yield stress investigation upon a given fluid model (e.g. the Bingham model or the Casson model).^20,21 The combination of yield stress measurements with AFM techniques can be used to well-characterize the nature of weak particle attractions and surface forces at nN scales.⁸ Further issues of a more geometrical nature, which naturally connect to τ₀, are rheological percolation²² and its differences from other connectivity phenomena, such as the onset of electric²³ or elastic percolation.^24,25 In granular fluids, it relates with the theory of jammed states,²⁶ originally pioneered by Edwards.²⁷

In nanoscience as well, the stability control and characterization in single and mixed dispersions or melts is an important and complex step.^28,29 Carbon nanotube suspensions,³⁰ for example, can be prepared in association with other molecular systems, like surfactants and polymers^31–33 or by (either covalent or non-covalent) functionalization of their walls with reactive groups, which increases the chemical affinity with dispersing agents.³⁴ As a consequence of large molecular aspect ratios and significant van der Waals’s attractions, the nanotube aggregation is highly enhanced, giving rise to strongly anisotropic systems of crystalline ropes and entangled network bundles, which are difficult to exfoliate, suspend or even characterize.³⁵ Stable CNT dispersions of controlled molecular mass may also exhibit polymeric behavior, and be quantitatively studied by equations taken from the well-established science of macromolecules.^36,37

This paper puts forward a basic approach, mostly focused on equilibrium arguments, to devise a yield stress law connected with particle statistics. By conjecturing an ensemble of effective volumes ‘displaced’ at the incipient state of motion, a statistical mechanics picture of τ₀ is proposed. This affords a phenomenological hypothesis that can be developed with reasonable simplicity. The derived relations are applied to typical disperse systems in colloid science and soft matter, such as aqueous and nonaqueous suspensions of ceramic/metal oxides and nanoparticles.

Defining an ensemble of volumes at the incipient motion

Yield stresses can be generally written as a sum of pairwise bonding contributions, where each particle pair is assigned a larger volume than the juxtaposition of the two initial units,^38,39 or by averaging the geometric part of the Hamaker expression over a representative pairwise cell.⁶ This analysis relies instead on a statistical definition of an effective volume concept (V_d) in a thought experiment. We suppose that, in yielding a dispersion, a canonical ensemble of volumes is displaced from the rest configuration and fulfils well defined statistical laws in thermal equilibrium, at a constant particle number and dispersion volume. The energy perturbation at the onset of motion will diminish with an increasing solid fraction in a representative large cell V_d, reflecting an increase in τ₀.

The disperse system will be conjectured to consist of an ensemble of elementary cells, each containing or not containing at least an aggregated solid unit (onwardly referred to as the “cluster” = cell + aggregates). An equilibrium total cluster number then can be identified with a conserved sum of aleatory variables (N_k):

whose values specify the aggregation statistics in each cell k. In the simplest limit cases, one may assign (A) a two-valued set, N_i = {0,1}, meaning that the i-th cell will be respectively empty or occupied by an aggregate, or (B) it may take any integer value, N_i = {0,1, … ∞}, defining in the extreme case, a cell that would be capable of hosting an open aggregate number. We complete this definition by a second relationship for the portions of fluid displaced in the experiment, conceived for simplicity as a set of discrete terms obeying the following combination law of volumes:where ρ_k = N_k/V_k is the k-th cluster density, and V_k is the k-th contribution to V_d. Yield strength will reflect the configuration statistics over the ensemble. Any cell containing aggregates will contribute to it, otherwise it won’t affect τ₀ (see Fig. 1). In ‘simple’ yield stress fluids, rest interactions are known to prevent the aggregation structure from breaking as a consequence of thermal agitation, and slow flows display a plastic behavior at very large deformations, without irreversible structural variations.^1,40 However, the equation systems (1) and (2) will be supposed to generally hold and be adopted irrespective of specific fluid dynamics properties (e.g. thixotropy or shear-thinning, pseudo-plasticity, etc.). Note that the density concept which the first two equations refer to does not coincide with the average dispersion density. The aim is to define an ad hoc (i.e. apparent, quasi-static) solid fraction value by which the yield stress response can be obtained by a statistical mechanics approach at the incipient state of motion.

Fig. 1 Scheme of the fluid cluster structure at the incipient motion. Cell volumes aren’t necessarily equal, as depicted here. Those in grey, with cell occupancy larger than zero, contribute to eqn (2).

A criterion assigning a sum over aggregation states is required for a thermodynamic framework, with this being promptly done by the usual partition function concept:

To give a suitable representation of the statistical issue, we define P_i = P_i(v,V_i) to be the probability of finding in V_i an empty (liquid) portion of volume v or, equivalently, an empty volume fraction ψ_i ≡ v/V_i. As we will tackle enough concentrated systems, this choice is suitable for an application of Poisson’s statistics. For spherically symmetric units and negligible excluded volumes, an expression like:

P_i ≈ e^−εvρ_i (4)

may be adopted,⁴¹ implying:with normalized canonical probabilities, constrained to eqn (1) and (2):

The quantity ε is introduced as a measure of cluster–cluster interaction strength,⁴² and regarded for simplicity as a homogeneous stiffness parameter, independent of {ρ_i}. The volume fraction ψ_i is a characteristic of the implied aggregation state.

We will relate the average occupancy in a generic single-particle state to the cluster interaction extent, and thus to τ₀. Particle indistinguishability, which is commonly a non-classical feature (but not only, see e.g. ref. 43), will be retained in this framework. From the rules of statistical mechanics, one therefore gets:

that, in light of eqn (2), may be rewritten as:Here, Z_i denotes the partition function with state i omitted:i.e. corresponding to the restricted sum:

Developing eqn (8) returns a distribution of (A) Fermi–Dirac or (B) Bose–Einstein type, and the proof of such a formal analogy is resumed in Appendix 1 upon mapping:

E_k = εk_BTψ_k (11)

with k_BT being the Boltzmann thermal energy. This phenomenological equation redefines an effective volume fraction in the displaced fluid (Fig. 2). It has an intuitive significance, as energy perturbations (E_k) will increase with increasing temperature, stiffness parameter and the velocity perturbation at the onset of motion, reflected by a larger liquid fraction (ψ_k), i.e. conditions normally implying a smaller yield stress. In the next section, some observations on fractional statistics, lying between (A) and (B), are reported as well the introduction of an intermediate state (C).

Fig. 2 Sketch of the canonical volume ensemble. Total volume V + V*, energy E + E*, and individual particle numbers are held constant, with n ≪ n*.

Finally, an equivalence between volume and energy was formerly introduced by Edwards and Oakeshott, although in a different context.²⁷ Their pioneering theory of powders faced the issue of a statistical mechanics analysis of non-thermal systems, like granular fluids.

Yield stress and cluster statistics

The most general pressure equation is written as a microscopic average of the stress tensor over the statistical distribution of particle states, provided here by eqn (7).⁴⁴ As yield stresses are likewise expected to be linear combinations of 〈N_i〉 with unknown (tensor) coefficients, we specialize the calculation to the mean cluster occupancy in the overall representative state conjectured by eqn (2):i.e.:

τ₀ ∼ 〈N(ψ)〉 (13)

The proportionality constant in this relationship will be removed by forming the experimental quantity τ* = τ₀/τ_M, i.e. dividing by the largest stress within a class of homogeneous and comparable measurements.

In (A), with a two-valued statistics, the arrangement of particle clusters takes the form (Appendix 1):

〈N_i〉 = (e^a+εψ_i + 1)⁻¹ (14)

with a being a characteristic fluid property, expressible as:

a = −εψ_a (15)

for some aggregation state denoted by ψ_a. It describes a liquid volume fraction at which the distribution function may either show a phase-like transition or a (rheological) percolation point.²³ Correspondingly, ε gives a measure of the average rate at which τ₀ is changing near ψ_a. The yield stress value as a function of ψ thus reads:

τ_0,1(ψ) = [e^{ε(ψ−ψ_a)} + 1]⁻¹ (16)

and, as the solid fraction θ complementing ψ obeys:

ψ + θ = ψ_a + θ_a ≡ 1 (17)

the normalized distributions of values in θ will scale as:still with , with θ_a being a critical threshold/maximum packing density, and:

ξ(θ) = e^{−ε(θ−θ_a)} (19)

In (B), the expression for an infinitely-valued statistic is regained as:

〈N_i〉 = (e^a+εψ_i − 1)⁻¹ (20)

so that:

Note that, as this distribution approaches infinity for ξ → 1⁺, it can no longer be normalized to unity. The idealized situation in which cells are permitted to be indefinitely occupied corresponds to the steepest liquid–solid transitions or percolation points contemplated by the model in its present form.

When the maximum cell occupancy is neither unitary nor infinite, it can be proven that:⁴⁵

upon validity of the implicit relation w_α^α(1 + w_α)^1−α = ξ(θ), α ∈ [0,1] ∩ Q. This parameter, α ≤ 1/〈N_r〉, gives a terse summary of the fractional interaction state, the extreme limits of which return type A (α = 1) and B (α = 0) laws. In between, an intermediate ‘statistical interaction’ follows, being attractive or repulsive depending respectively on whether or .⁴⁶ A general treatment of α-states would raise tough formal and numerical issues, falling beyond the experimental verification purposes of this work. In the next section, we thus limit ourselves to the neutral state (type C), which can be solved explicitly (e.g. ) to return .⁴⁵ The semi-heuristic position w_α′(ξ) ≡ ξ, now with α′ ∈ [−1,1], would be rather convenient for an intuitive data analysis, but is formally inadequate for a model comparison, and is therefore disregarded.

Results and discussion

Yield stress measurements are normally reported against solid fractions like the volumetric (ϕ) or the mass concentration. To get the dependence ξ = ξ(ϕ), consider a fluid layer bounded on one side by the plane surface xz. For a flow perturbation implied by a constant velocity δU_ex exerted at the surface, the velocity component x at a distance y from the plane will obey the profile , in which ν is the kinematic viscosity and δt is a small time interval.⁴⁷ The complementary error function, denoted by Φ⁻, quantifies the velocity fraction transmitted to the fluid upon shearing, and thus can be adopted for the evaluation of θ illustrated in Fig. 3 and explained in detail in Appendix 2.

Fig. 3 Fluid dynamics scheme for evaluating θ.

If the velocity profile is averaged over an incompressible liquid flow undergoing free momentum diffusion (no pressure gradient or volume forces, nor hydrodynamic perturbations, see the remarks in Appendix 2), the solid fraction can be estimated as:

θ ≈ 〈Φ⁻〉 (23)

i.e.:with (i = d, dispersion; l, liquid) and P being a Gaussian function. This integral turns out to be well defined, as it is time-independent, and only depends upon ϕ through the ratio k_d/k_i:Here [small nu, Greek, macron] ≡ ν_d/ν_l = [small eta, Greek, macron]/[small rho, Greek, macron], with shear viscosity and mass density being expressible as η_d = η_l[small eta, Greek, macron](ϕ) and ρ_d = ρ_l[small rho, Greek, macron](ϕ), i.e. the product of pure liquid properties times a function of the solid fraction. The reduced quantity [small nu, Greek, macron] = [small nu, Greek, macron](ϕ) generally increases with increasing ϕ, since viscosity changes should dominate over ρ values in the concentration regimes of interest. Correspondingly, a larger fluid inertia leads to a reduction of ψ (eqn (16)), as expected. Note that a critical threshold/maximum packing density (ϕ_a) cannot be identified from the last equation, being . A boundary/cutoff value θ_a ≡ ϕ_a needs to be set independently, completing the definition of ξ = ξ(θ) in the previous relationships, with the requirement θ < ϕ_a, or:

[small nu, Greek, macron](ϕ) < tan(πϕ_a) (26)

representing a model constraint for the viscosity, density and the critical volume fraction. Applications of the new equation family will be conducted in conformity with it, as shown in Fig. 7. While eqn (26) always depends on the viscosity model, it is more selective for B and C statistics (α < 1), describing steeper liquid–solid transitions.

Eqn (25) redefines an effective volume fraction, ϕ → θ(ϕ). To test its validity, experimental measurements from aqueous and non-aqueous systems were taken from the literature. In the first case, the yield stress of ceramic and metal oxide dispersions of Si₃N₄,^5,6 α-Ca₃(PO₄)₂,^5,6 ZrO₂,⁴⁸ and TiO₂ (ref. 49) (anatase) were regarded, whereas Al₂O₃/decalin⁵⁰ and multiwalled carbon nanotubes/polycarbonate⁵¹ (MWCNT/PC melts) were regarded as non-aqueous materials. In these systems, interaction mechanisms were mostly London–van der Waals with absent or irrelevant Vold’s effect.⁵² The first three were at/near their isoelectric points, where double-layer surface charges are negligible.⁵³ The structure of the anatase colloids was governed too by attractive van der Waals forces, which represents the main interaction mechanism in the last system as well.⁵⁴ Compressive yield stress in the fifth system, alumina in decalin, was still of the van der Waals type, modulated by a steric interparticle repulsive barrier of ∼0.7 nm of propionic acid, with no further electrostatic or structural energies coming into play. With the obvious exception of MWCNT/PC melts (T ≃ 533 K), all measurements were conducted at room temperature. The adopted extrapolation laws were Casson or Bingham, unambiguously written in the yield stress, viscosity, and shear rate, with no heuristic constants. In Al₂O₃ systems, compressive yield stresses were numerically inferred from the measured variation of volume fraction solids with elevation. Any further physical chemistry details may be found in the references.

Every plot of vs. ϕ was best fitted by means of each particle statistic, 1/α = 1 (A), ∞ (B), 2 (C) and the extracted model parameters are depicted in Table 1. Shear viscosity data vs. solid concentration were available for the system MWCNT/PC⁵⁵ and fulfilled Eiler’s law, (e.g. ref. 56), where a large intrinsic viscosity value ([η] = 165) seems to be a feature of other carbon nanotube suspensions.⁵⁷ Quemada’s model, [small eta, Greek, macron](ϕ) = (1 − ϕ/ϕ_a)⁻², was generally adopted in the other systems, still producing a good agreement with yield stress profiles. Examples of the suitability of each model is shown in Fig. 4–6, while Fig. 7 shows how the model constraint in eqn (26) is fulfilled in these cases. Obviously, a similar agreement is also found for Ca₃(PO₄)₂/H₂O, ZrO₂/H₂O, and Al₂O₃/C₁₀H₁₈.

Table 1 Yield stress model parameters

Chemical system (s/ l)	α	ϕa	ε
Si3N4/H2O		0.452	31.7
Ca3(PO4)2/H2O	0	0.509	10.5
ZrO2/H2O		0.393	26.1
TiO2/H2O	0	0.523	18.6
Al2O3/C10H18	0	0.527	30.9
C/(–O–(C=O)–O–)n	1	0.443	69.8

---

Fig. 4 Example of particle statistics α = 0. Reduced yield stress versus solid volume fraction for the anatase system TiO₂/H₂O (densely dashed line, α = 0; dashed line, α = 1; solid line, ). Model parameters for α = 0 are in Table 1. Best fits with similar quality were also met in Ca₃(PO₄)₂/H₂O and Al₂O₃/C₁₀H₁₈ systems.

Fig. 5 Example of particle statistics . Reduced yield stress versus solid volume fraction for Si₃N₄/H₂O. Lines and symbols are as in Fig. (4), and model parameters for are in Table 1. Best fits with similar quality were also obtained in ZrO₂/H₂O.

Fig. 6 Example of particle statistics α = 1. Reduced yield stress versus solid volume fraction for MWCNT/PC. Lines and symbols are as in Fig. 4, and model parameters for α = 1 are in Table 1.

Fig. 7 Representation of eqn (26) in each of the instances in Fig. 4–6 with [small rho, Greek, macron] = 1 + (ρ_s/ρ_l − 1)ϕ. From the top, TiO₂/H₂O (ρ_s/ρ_l ≈ 4.23, short-dashed), Si₃N₄/H₂O (ρ_s/ρ_l ≈ 3.17, medium-dashed) and MWCNT/PC (ρ_s/ρ_l ≈ 1.17, long-dashed). As expected, each of the functions θ–ϕ_a is always positive in the respective experimental domain, ϕ ≤ 0.12 (TiO₂), ϕ ≤ 0.32 (Si₃N₄), and ϕ ≤ 3 × 10⁻² (MWCNT).

Single-valued distribution (A) is the only one predicting a plateau, as is the evident case here of MWCNT/PC melts. A similar situation arises from kaolin colloids in water, paraffinic oil and liquid polybutadiene rubber, for which an S-shaped functional form analogous to eqn (18) was proven to hold.^58,59 Unluckily, the sparseness of physical chemistry properties of kaolin powders, especially when commercially supplied, do not allow for a data comparison with the other chemical systems, in Tables 1 and 2. Furthermore, the present model turns out to be fundamentally non-linear. Since chemical compositions of kaolin materials are markedly heterogeneous, to average over distinct solid components would require an extension of this approach to the framework of multi-phase media. Finally, on increasing α, (B) and (C) distributions imply steeper behaviors, as they are usually indicative of stronger interactions. Improving the agreement in Fig. 6 likely requires the combination of further interaction states.

The maximum solid loading is known exactly in two cases, the first of which, ϕ_m = 0.54 (Al₂O₃/C₁₀H₁₈),⁵⁰ agrees with ϕ_a ≈ 0.53 in Table 1. Volume corrections due to propionic acid layers were estimated to be less than 1%. The second, ϕ_m = 0.146 (TiO₂/H₂O), was highly affected by porosity (nanoparticles forming an interlinked porous network in the liquid) and thus is not directly comparable with ϕ_a.⁴⁹ The fraction ϕ_a was reported at the third digit, as the model was sensitive to it.

To interpret the model parameters by the extent of particle interactions, relevant energetic quantities at mesoscopic/macroscopic scales were collected in Table 2, i.e. the Hamaker constant for particles in the dispersant (A_H), which reflects the difference between the polarizabilities of solid and liquid molecules, and Young’s/bulk moduli of the solid phase (E, K), defining the elastic response in homogeneous isotropic materials. While A_H is known to relate with τ₀,^6,48 elastic (and loss) moduli, detected from oscillatory tests at ∼1 Hz, were recently found to be proportional to the yield stress of gel (Carbopol) solutions.⁹

Table 2 Particle–particle interactions and elastic solid constants

Chemical system (s/l)	AH [zJ]	E [GPa]	K [GPa]
Si3N4/H2O	46	320	245
Ca3(PO4)2/H2O	23	104	76
ZrO2/H2O	88	241	181
TiO2/H2O	37	167	192
Al2O3/C10H18	20	400	255
C/(–O–(C=O)–O–)n	8	740	442

---

The aqueous Hamaker constant was determined by the knowledge of optical spectra of Si₃N₄ (ref. 60) and ZrO₂ (ref. 61) from Lifshitz theory. For TiO₂, it was deduced from linearly correlating yield stress data with the square zeta potential.⁶² As in the case of colloidal Al₂O₃ in decalin,⁵⁰ the Tabor–Winterton approximation was adopted to evaluate A_H of tricalcium phosphate particles. Upon neglecting retardation effects, it gives:⁶³

with refractive indices in the visible spectrum set to n_s ≈ 1.6, and n_l = 1.33 and static dielectric constants ε_s ≈ 15.4, and ε_l = 80. The UV absorption frequency should range here in (ref. 63) (3–4) × 10¹⁵ Hz so that a mean v_e ≈ 3.5 × 10¹⁵ Hz was assumed. Because data for α-Ca₃(PO₄)₂ are scarce, n_s is inferred from optical spectra of calcium phosphate glasses⁶⁴ in the proportion [CaO] : [P₂O₅] = 1 : 3 and ε_s is that of hydroxyapatite, as it was successfully employed in a DLVO theory for amorphous particles.⁶⁵ The Hamaker constant for MWCNT/PC was also evaluated, from the combining law⁶⁶ for the pure materials values A_sH ≈ 100 zJ (ref. 67) and A_lH = 50.8 zJ:⁶⁸

Elastic features were taken from rather recent literature on each of the solid compounds listed in the tables. For Si₃N₄, an average E value among three families of samples was considered, which is then representative of Ceradyne Ceralloy.⁶⁹ Near-theoretical density values were adopted for Al₂O₃ (ref. 69) and K of β-Si₃N₄,⁶⁹ which was the most abundant silicon nitride phase (∼84.4% wt) in the original fluid samples of Fig. 5.⁷⁰ The Young’s modulus for α-Ca₃(PO₄)₂ was experimentally measured, with the bulk one being calculated instead using ab initio density functional theory calculations (DFT),⁷¹ i.e. the same numerical framework employed to get the elastic constants of TiO₂ (anatase)⁷² and ZrO₂.⁷³ Average values between the predictions from generalized gradient (GGA) and local density approximations (LDA) were regarded for the monoclinic ZrO₂ phase. Concerning MWCNT, E was taken from an atomistic potential simulation,⁷⁴ whereas the diamond value was used for K.⁷⁵

To proceed, the quantitative insight now is to link A_H to ϕ_a, and E, K to ε. On increasing the extent of the repulsive forces, particles exhibit an effective size (r_e) and volume fraction (ϕ_e), relating with unperturbed values as⁷⁶ ϕ/ϕ* = (r/r*)³.

Accordingly, by representing ϕ_a = ϕ_a(A_H), one can verify a meaningful linear correlation that extrapolates, in the limit A_H → 0⁺, a good estimate for the packing fraction of the simple cubic cell, ϕ_a (0⁺) = 0.517 against π/6 ≈ 0.524 (Fig. 8). Regarding ε = ε(E, K), it still shows a reasonable linear trend with both elastic constants (Fig. 9 and 10). As expected upon E, K → 0⁺, the stiffness parameter ε(0⁺) → 0 as well. The behaviors in Fig. 8–10 did not vary appreciably upon changing the yield stress data from Bingham to Casson, which were both available for anatase.⁴⁹

Fig. 8 Packing fraction versus the Hamaker constant, best fitted by ϕ_a = 10⁻³ × (517.2 − 1.2A_H). Data from Table 2.

Fig. 9 Stiffness parameter versus the Young’s modulus, best fitted by ε = 9.3 × 10⁻²E [GPa]. Data from Table 2.

Fig. 10 Stiffness parameter versus the bulk modulus, best fitted by ε = 14.0 × 10⁻²K [GPa]. Data from Table 2.

Such correlations, especially with the Young’s modulus, are suggestive of a linear law between ε and the interatomic spring constant of the solid compounds, consistent with a microscopic definition of the stiffness parameter, i.e.:

where V = V(r) is the bond energy and its curvature is calculated at the equilibrium interatomic distance r = r₀. The constant term on the right identifies an intrinsic elastic modulus, as it would promptly follow from the knowledge of potential energy details (e.g. Lennard-Jones) combined with the best fit of ε vs. E.

In conclusion, it is noteworthy to point out a couple of formal issues met in this analysis. First, the equations for τ₀ cannot be reduced to a power law in ϕ where, as in fractal-like structures, the exponent is expected to account for the interaction rate.⁷⁶ This doesn’t mean, evidently, that (B) or (C) distributions would not yield good numerical approximations for such cases as well. Second, the model inadequacy to reproduce in the limit of infinite dilution is anyway recovered upon ε → ∞. The actual values of stiffness parameters should guarantee a rather fast convergence since, for (A) and (B) distributions with and :

the numbers in Table 1 reasonably return when ϕ → 0⁺.

Conclusive remarks

Solid concentration affects the properties of complex fluids, and we put forward the (re)definition of an effective volume fraction (θ) for the yield stress behavior, here evaluated in terms of a reduced kinematic viscosity. Aggregation clusters contributing to τ₀ are modelled by canonical ensembles of (displaced) volumes, with particle statistics determined by the occupancy number.

This conjecture points out an average cluster that is representative of the incipient state of motion, with given liquid and solid fractions. A class of statistical mechanics laws defines yield stress in terms of two coefficients, the maximum packing threshold (ϕ_a) and the particle stiffness parameter (ε), e.g.:

1/τ₀(ϕ) ∼ e^{−ε[θ(ϕ)−ϕ_a]} ± 1 (31)

which turn out to relate to the Hamaker constant and elastic moduli of the solid phase.

Conflicts of interest

There are no conflicts to declare.

Appendixes

Appendix 1

We resume here some of the well known results for eqn (8) in (A) and (B).⁷⁷ In comparison to the original theories, the energy of the level k(E_k) is mapped initially into ψ_k, the Boltzmann energy (k_BT) into 1/ε and the chemical potential (μ) into ψ_a. In fact, the energy representation is regained upon eqn (11).

Consider thus the restricted partition function specified by eqn (9) and (10). In (A), the mean particle number in a state i reads:

where, to work it out, the following expansion can be used upon Δn/n ≪ 1:

lnZ_i(n − Δn) ≈ lnZ_i(n) − a_iΔn (A.2)

and, for a sum over many states, the following approximation:implies:

Replacing this result into eqn (A.1) proves eqn (14). In (B), let’s use the former logarithmic expansion:

to get:the numerator of which is expressible through:

Obviously, after noting the infinite geometric series:

Eqn (20) is recovered at once, since:

Remember that eqn (A.3) upon eqn (11) gives a = −μ/(k_BT). We accordingly interpret the ratio −a/ε in the corresponding relationship to assign a characteristic aggregation state of the disperse system (a percolation-like point ψ_a) in every statistic here regarded.

Appendix 2

For the incompressible flow in Fig. 3, with δu_y = δu_z = 0, with no pressure gradient and volume forces, the balanced equation for the momentum along x brings us to:where the constants, ρ and η, denote the mass density and shear viscosity coefficient. This is the basic equation to get the complementary error function and the Gaussian probability distribution in eqn (24).

To evaluate the fraction θ in eqn (17), we form the quantity:

with each displacement at the numerator being time-dependent, with δX_i = δx_i (0). Since no perturbation develops along y and z, one has δy/δY = δz/δZ = 1, and:where δx/δu_x = δX/δU = δt. To average the solid fraction and eliminate the dependencies on time and space, we limit ourselves to the free momentum distribution (P) in the liquid flow and disregard perturbations from long-ranged hydrodynamic interactions among particles.⁷⁸ They are expected to be negligible upon decreasing speed at the incipient motion⁷⁹ and increasing dilution of the dispersed units.⁸⁰ Hydrodynamic back-flows may be effectively screened as well by charged particles,⁸¹ while in systems like complex fluid interfaces they are generally not.⁸²

Therefore, still with the same notations of eqn (24), it turns out that (x_i = x,y,z):

which, since:and:reduces, after some mathematical developments, to eqn (25).

Acknowledgements

Jordi Sancho-Parramon and Prabhu Nott are acknowledged for their useful remarks. This work started years ago and was partly supported by a former Italian MURST (40%).

References

1. P. Coussot, Rheol. Acta, 2017, 56, 163–176.
2. P. Coussot, Rheol. Acta, 2018, 57, 1–14.
3. R. R. Huilgol, Phys. Fluids, 2002, 14, 1269–1283.
4. A. Wachs and I. A. Frigaard, J. Non-Newtonian Fluid Mech., 2016, 238, 189–204.
5. S. A. Mezzasalma, J. Colloid Interface Sci., 1997, 190, 302–306.
6. S. A. Mezzasalma, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 3134–3141.
7. G. V. Franks, J. Colloid Interface Sci., 2002, 249, 44–51.
8. E.-J. Teh, Y.-K. Leong, Y. Liu, V. S. J. Craig, R. B. Walsh and S. C. Howard, Langmuir, 2010, 26, 3067–3076.
9. J. Boujlel and P. Coussot, Soft Matter, 2013, 9, 5898–5908.
10. A. R. Koblitz, S. Lovett and N. Nikiforakis, Physical Review Fluids, 2018, 3, 023301.
11. H. A. Barnes, J. Non-Newtonian Fluid Mech., 1999, 81, 133–178.
12. P. Coussot, Q. D. Nguyen, H. T. Huynh and D. Bonn, Phys. Rev. Lett., 2002, 88, 175501.
13. M. Jenny, S. Kiesgen de Richter, N. Louvet, S. Skali-Lami and Y. Dossmann, Physical Review Fluids, 2017, 2, 023302.
14. B. Coleman, H. Markovitz and W. Noll, Viscometric Flows of Non-Newtonian Fluids: Theory and Experiment, Springer, Berlin, 1966.
15. P. Samyn and H. Taheri, J. Mater. Sci., 2016, 51, 9830–9848.
16. F. Yang, K. Paso, J. Norrman, C. Li, H. Oschmann and J. Sjöblom, Energy Fuels, 2015, 29, 1368–1374.
17. C. Zhao, G. Yuan, D. Jia and C. C. Han, Soft Matter, 2012, 8, 7036–7043.
18. A. Chakrabarty and U. Maitra, J. Phys. Chem. B, 2013, 117, 8039–8046.
19. G. Ovarlez, Q. Barral and P. Coussot, Nat. Mater., 2010, 9, 115.
20. H. J. Choi and M. S. Jhon, Soft Matter, 2009, 5, 1562–1567.
21. S. A. Mezzasalma and G. J. Koper, Colloid Polym. Sci., 2002, 280, 160–166.
22. A. K. Kota, B. H. Cipriano, M. K. Duesterberg, A. L. Gershon, D. Powell, S. R. Raghavan and H. A. Bruck, Macromolecules, 2007, 40, 7400–7406.
23. D. Stauffer and A. Aharony, Introduction to Percolation Theory, 1985.
24. S. Arbabi and M. Sahimi, Phys. Rev. Lett., 1990, 65, 725–728.
25. M. Sahini and M. Sahimi, Applications Of Percolation Theory, CRC Press, London, 1994.
26. P. Ball, Nat. Mater., 2016, 15, 1227.
27. S. Edwards and R. Oakeshott, Phys. A, 1989, 157, 1080–1090.
28. C. Zhao, G. Yuan, D. Jia and C. C. Han, Soft Matter, 2012, 8, 7036–7043.
29. M. T. Zafarani-Moattar and R. Majdan-Cegincara, Fluid Phase Equilib., 2013, 354, 102–108.
30. E. K. Hobbie and D. J. Fry, J. Chem. Phys., 2007, 126, 124907.
31. C. A. Mitchell and R. Krishnamoorti, Macromolecules, 2007, 40, 1538–1545.
32. M. J. Kayatin and V. A. Davis, Macromolecules, 2009, 42, 6624–6632.
33. Y. Dror, W. Pyckhout-Hintzen and Y. Cohen, Macromolecules, 2005, 38, 7828–7836.
34. V. Georgakilas, K. Kordatos, M. Prato, D. M. Guldi, M. Holzinger and A. Hirsch, J. Am. Chem. Soc., 2002, 124, 760–761.
35. R. Bandyopadhyaya, E. Nativ-Roth, O. Regev and R. Yerushalmi-Rozen, Nano Lett., 2002, 2, 25–28.
36. M. S. P. Shaffer and A. H. Windle, Macromolecules, 1999, 32, 6864–6866.
37. M. L. Usrey, A. Chaffee, E. S. Jeng and M. S. Strano, J. Phys. Chem. C, 2009, 113, 9532–9540.
38. R. J. Flatt and P. Bowen, J. Am. Ceram. Soc., 2006, 89, 1244–1256.
39. R. J. Flatt and P. Bowen, J. Am. Ceram. Soc., 2007, 90, 1038–1044.
40. P. Coussot, L. Tocquer, C. Lanos and G. Ovarlez, J. Non-Newtonian Fluid Mech., 2009, 158, 85–90.
41. W. Strieder and R. Aris, Variational Methods Applied to Problems of Diffusion and Reaction, Springer-Verlag, Berlin Heidelberg, Germany, 1973.
42. A. Coniglio, U. D. Angelis and A. Forlani, J. Phys. A: Math. Gen., 1977, 10, 1123.
43. A. Bach, Indistinguishable Classical Particles (Lecture Notes in Physics, Monograph 44), Springer-Verlag, Berlin Heidelberg, New York, Germany, 1997.
44. L. Durand, Am. J. Phys., 2004, 72, 1082–1094.
45. G. Joyce, S. Sarkar, J. Spalek and K. Byczuk, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 990–993.
46. Y.-S. Wu, Phys. Rev. Lett., 1994, 73, 922–925.
47. R. Bird, W. Stewart and E. Lightfoot, Transport Phenomena, John Wiley & Sons, New York, US, 1960.
48. Y. Leong, P. Scales, T. Healy, D. Boger and R. Buscall, J. Chem. Soc., Faraday Trans., 1993, 89, 2473–2478.
49. W. J. Tseng and K.-C. Lin, Mater. Sci. Eng., A, 2003, 355, 186–192.
50. L. Bergström, C. Schilling and I. Aksay, J. Am. Ceram. Soc., 1992, 75, 3305–3314.
51. S. Jin, D. Choi and D. Lee, Colloids Surf., A, 2008, 313–314, 242–245.
52. N. Dzuy and D. Boger, J. Rheol., 1985, 29, 335–347.
53. R. Pradip and S. Malghan, J. Am. Ceram. Soc., 1996, 79, 2567–2576.
54. S. Pande, A. Chaudhary, D. Patel, B. Singh and R. Mathur, RSC Adv., 2014, 4, 13839–13849.
55. D. Dijkstra, M. Cirstea and N. Nakamura, Rheol. Acta, 2010, 49, 769–780.
56. S. A. Mezzasalma, Colloid Polym. Sci., 2001, 279, 22–32.
57. K. Donovan and K. Scott, J. Chem. Phys., 2013, 138, 244902.
58. D. T. Beruto, A. Lagazzo, R. Botter and R. Grillo, Particuology, 2009, 7, 438–444.
59. D. Beruto, A. Lagazzo and R. Botter, Colloids Surf., A, 2012, 396, 153–160.
60. H. D. Ackler, R. H. French and Y.-M. Chiang, J. Colloid Interface Sci., 1996, 179, 460–469.
61. L. Bergstrom, A. Meurk, H. Arwin and D. Rowcliffe, J. Am. Ceram. Soc., 1996, 79, 339–348.
62. A. Gòmez-Merino, F. Rubio-Hernàndez, J. Velàzquez-Navarro, F. Galindo-Rosales and P. Fortes-Quesada, J. Colloid Interface Sci., 2007, 316, 451–456.
63. J. N. Israelachvili, Intermolecular and Surface Forces, 1991.
64. E. Lee and E. Taylor, Opt. Mater., 2006, 28, 200–206.
65. T. Morgan, T. Goff and J. Adair, Nanoscale, 2011, 3, 2044–2053.
66. J. N. Israelachvili, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1972, 331, 39–55.
67. M.-F. Yu, T. Kowalewski and R. Ruoff, Phys. Rev. Lett., 2001, 86, 87–90.
68. R. French, J. Am. Ceram. Soc., 2004, 83, 2117–2146.
69. D. de Faoite, D. Browne, F. Chang-Díaz and K. T. Stanton, J. Mater. Sci., 2012, 47, 4211–4235.
70. S. Mezzasalma and D. Baldovino, J. Colloid Interface Sci., 1996, 180, 413–420.
71. L. Liang, P. Rulis and W. Ching, Acta Biomater., 2010, 6, 3763–3771.
72. X. Ma, P. Liang, L. Miao, S. Bie, C. Zhang, L. Xu and J. Jiang, Phys. Status Solidi B, 2009, 246, 2132–2139.
73. X.-S. Zhao, S.-L. Shang, Z.-K. Liu and J.-Y. Shen, J. Nucl. Mater., 2011, 415, 13–17.
74. J. Uddin, M. Baskes, S. Srinivasan, T. Cundari and A. Wilson, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 104103.
75. J. Huntington, Solid State Physics, ed. F. Seitz and D. Turnbull, Academic, New York, 1958.
76. N. Willenbacher and K. Georgieva, Rheology of Disperse Systems, Wiley-Blackwell, 2013, pp. 7–49.
77. F. Reif, Fundamentals of Statistical and Thermal Physics, Mc Graw-Hill Book Company, New York, US, 1965.
78. S. A. Mezzasalma, Macromolecules in Solution and Brownian Relativity, Academic Press - Elsevier, Oxford, 2008.
79. P. Coussot, Rheophysics, Matter in all its States, Springer International Publishing, Switzerland, 2014.
80. C. Reichhardt and C. J. Olson Reichhardt, Phys. Rev. Lett., 2006, 96, 028301.
81. D. O. Riese, G. H. Wegdam, W. L. Vos, R. Sprik, D. Fenistein, J. H. H. Bongaerts and G. Grübel, Phys. Rev. Lett., 2000, 85, 5460–5463.
82. C. Klein, A. Theodoratou, P. A. Rühs, U. Jonas, B. Loppinet, M. Wilhelm, P. Fischer, J. Vermant and D. Vlassopoulos, Rheol. Acta, 2019, 58, 29–45.

---