Fundamental Rheology of Disperse Systems Based on Single-Particle Mechanics

A comprehensive review of the fundamental rheology of dilute disperse systems is presented. The exact rheological constitutive equations based on rigorous single-particle mechanics are discussed for a variety of disperse systems. The different types of inclusions (disperse phase) considered are: rigid-solid spherical particles with and without electric charge, rigid-porous spherical particles, non-rigid (soft) solid particles, liquid droplets with and without surfactant, bubbles with and without surfactant, capsules, core-shell particles, non-spherical solid particles, and ferromagnetic spherical and non-spherical particles. In general, the state of the art is good in terms of the theoretical development. However, more experimental work is needed to verify the theoretical models and to determine their range of validity. This is especially true for dispersions of porous particles, capsules, core-shell particles, and magnetic particles. The main limitation of the existing theoretical developments on the rheology of disperse systems is that the matrix fluid is generally assumed to be Newtonian in nature. Rigorous theoretical models for the rheology of disperse systems consisting of non-Newtonian fluid as the matrix phase are generally lacking, especially at arbitrary flow strengths.


Introduction
Understanding of the rheological behavior of particulate fluids (disperse systems) is important from a practical point of view.Particulate fluids form a large group of materials of industrial importance.Solids-in-liquid suspensions, emulsions (blends of two immiscible liquids), and foams (bubbly dispersions) are examples of particulate fluids.Both aqueous and non-aqueous particulate fluids are encountered in industrial applications [1][2][3][4][5].In the handling, mixing, storage, and pipeline transportation of fluidic particulate materials, knowledge of the rheological behavior is required for the design, selection, and operation of the equipment involved.
In this article, the fundamental rheology of dilute particulate fluids based on a single-particle mechanics is reviewed and discussed.The discussion covers particulate fluids composed of a variety of different types of particles, such as:

•
Rigid-solid spherical particles: uncharged and electrically charged   Ferromagnetic rigid-solid spherical and non-spherical particles

Bulk Stress and Bulk Rate of Strain in Particulate Fluids
Although particulate fluids are heterogeneous systems consisting of discrete particles of one phase suspended in a continuum of another phase, they can be regarded as effectively homogeneous systems when it comes to defining their bulk (macroscopic) rheological properties [6][7][8][9], provided that the length scale of the apparatus (L) is large in comparison with the average spacing ( ) between the centers of the adjacent particles.The bulk or average fields can then be defined as volume averages of the corresponding local fields over a representative sample element of volume V whose linear dimension, V 1/3 , is large compared with the average spacing between the particles , but small compared with the apparatus length scale L, that is, << V 1/3 << L. Under the condition << V 1/3 << L, the volume of the element V is large enough to contain a statistically significant number of particles but small relative to the characteristic macro-scale of the system.Thus, the bulk or volume-average stress tensor σ in particulate fluids can be defined as: where σ is the local stress tensor.The bulk rate of strain tensor E can be defined in a similar manner as: where E is the local rate of strain tensor defined as: In Equation (3), → u is the local velocity vector and superscript "T" refers to transpose of a quantity.

Rheological Constitutive Equation for Particulate Fluids
For a particulate fluid of identical particles (same shape, size, and internal constitution), the rheological constitutive equation can be expressed in the following general form [9]: where σ o is the stress tensor in pure matrix material under the condition corresponding to the imposed macroscopic rate of strain tensor E on the particulate fluid, n is the number density of particles in the particulate fluid, and S is the average value of dipole strength S for a reference particle.The dipole strength S of a particle could be interpreted as a measure of the increase in stress resulting from the replacement of matrix material by particle with the rate of strain held constant, that is: where V p is the volume of the particle, σ o is the local stress tensor in region V p when matrix material is present and σ is the local stress tensor when the matrix material is replaced by particle in the region V p , keeping the local rate of strain constant.The dipole strength S of a particle could also be expressed as [9,10]: where A p is the surface area of a particle, n is the unit outward normal to the particle surface, → u is the matrix fluid velocity at the particle surface, σ is the matrix fluid stress tensor at the particle surface, → X is the position vector with respect to a fixed origin, η c is the matrix fluid viscosity, and δ is a unit tensor.For a particulate fluid composed of identical particles of some shape, size, and internal constitution, the dipole strength of a reference particle is affected by two microstructural factors: (a) the orientation of the reference particle itself if the reference particle is non-spherical; and (b) the orientations and relative positions of the neighboring particles.However, for dilute particulate fluids, the interactions between the particles are negligible as the particles are far apart.The average value of the dipole strength S for the reference particle in dilute particulate fluids does not depend on the configuration of the other particles and therefore, S can be calculated assuming that the reference particle alone is situated in an infinite matrix medium.Thus, the rheological constitutive equation for dilute particulate fluids can be expressed as: (7) where S o is the mean dipole strength of the reference particle situated in an infinite matrix medium.
Note that in the case of non-spherical particles, S o depends on the orientation of the reference particle and therefore, S o is the mean value of different orientations.For dilute particulate fluids of spherical particles, the orientation effects are absent and consequently, the rheological constitutive equation could be expressed as: where φ is the volume fraction of particles and R is the particle radius.
It is clear from the preceding discussion that in order to develop the rheological constitutive equation for dilute particulate fluids, we need to determine the mean dipole strength S o of a reference particle situated in an infinite matrix medium.To evaluate the mean dipole strength, knowledge of the local flow fields (velocity, stress distribution) around the particle is needed (see Equation (6)).The local fields around the particle can be determined from the solution of the equations governing the motion of fluid around the particle, subject to the relevant boundary conditions.As an example, the local fields around a rigid-solid and electrically-neutral particle, suspended in a Newtonian incompressible fluid, can be determined from the solution of the Stokes equations given below subject to the following boundary conditions: (a) when r → ∞, → u → → u ∞ where → u ∞ is the fluid velocity far away from the particle; and (b) when r → R (that is, at particle surface), r where → ω is the angular velocity of the particle and → r is the position vector measured from origin fixed at the center of the particle:

Electrically-Neutral Particles
The unknown in the rheological constitutive equation of dilute particulate fluids of rigid-solid spherical particles (Equation ( 8)) is the dipole strength of a reference particle.As already pointed out,  6)) requires knowledge of the velocity and stress distributions around the particle.The velocity distribution around the reference particle is obtained by solving the Stokes equations (Equations ( 9)- (10)), subject to the relevant boundary conditions.Once the velocity distribution is known, the local rate of strain tensor is determined from Equation (3) and utilized in the following constitutive equation of the matrix fluid (assumed to be incompressible Newtonian fluid) to obtain the local stress field in the matrix fluid: where η c is the matrix fluid viscosity and E is the local rate of strain tensor.Upon substitution of the velocity and stress distributions at the particle surface in Equation ( 6), the following result is obtained for the dipole strength of a rigid-solid, electrically-neutral, spherical particle [9]: where E ∞ is the local rate of strain tensor far away from the particle.For a dilute suspension, the rate of strain tensor E ∞ far away from the particle can be equated to the bulk or imposed rate of strain tensor E on the suspension.From Equations ( 8) and (12), it follows that: where σ o is the stress tensor in pure matrix material (assumed to be incompressible Newtonian fluid) under the condition corresponding to the imposed macroscopic rate of strain tensor E on the particulate fluid, and P is the average pressure.To simplify the notation, the angular brackets < > can be dropped from Equation (13).Thus, the rheological constitutive equation of a dilute suspension of rigid-solid electrically-neutral spherical particles is as follows: where σ is the bulk stress tensor of the suspension and E is the imposed or bulk rate of strain tensor on the suspension.According to Equation (14), the suspension is a Newtonian fluid with an effective viscosity of η given as: Equation ( 15) is referred to as the Einstein equation as Einstein [11,12] was the first one to develop this equation.The relative viscosity (η r ) and intrinsic viscosity ([η]) of a suspension are defined as: [η] = Lim For a dilute suspension of solid, spherical, and electrically-neutral particles, η r and [η] are given as: [η] = 2.5 (19) Figure 1 shows comparison between the predictions of the Einstein equation (Equation (18)) and experimental data obtained for dilute suspensions of hard-spheres.Six sets of data are shown [13][14][15][16].The experimental data show excellent agreement with the predictions of the Einstein theory provided that φ ≤ 0.03.At higher volume fractions of the dispersed particles, the Einstein equation under predicts the relative viscosity.The deviation between the predicted and experimental values increases with the increase in φ.
Fluids 2016, 1, 40 5 of 52 For a dilute suspension of solid, spherical, and electrically-neutral particles, r  and ] [ are given as: [   (19) Figure 1 shows comparison between the predictions of the Einstein equation (Equation (18)) and experimental data obtained for dilute suspensions of hard-spheres.Six sets of data are shown [13][14][15][16].The experimental data show excellent agreement with the predictions of the Einstein theory provided that 0.03   .At higher volume fractions of the dispersed particles, the Einstein equation under predicts the relative viscosity.The deviation between the predicted and experimental values increases with the increase in  .It should be noted that the constitutive equation expressed in the form of Equation ( 14) is valid for dilute suspensions of rigid-solid and electrically-neutral spherical particles dispersed in a matrix of incompressible Newtonian fluid.When the matrix fluid is compressible in nature, the constitutive equation must also include the bulk viscosity of suspension as shown below: It should be noted that the constitutive equation expressed in the form of Equation ( 14) is valid for dilute suspensions of rigid-solid and electrically-neutral spherical particles dispersed in a matrix of incompressible Newtonian fluid.When the matrix fluid is compressible in nature, the constitutive equation must also include the bulk viscosity of suspension as shown below: where tr refers to trace of a tensor (trE = ∇ • → u ), η is the shear viscosity of suspension given by the Einstein equation (Equation (15)) and η k is the bulk viscosity of suspension given by the following expression developed by Brady et al. [17]: where η k c is the bulk viscosity of a matrix fluid.

Electrically-Charged Particles
The particulate fluids consisting of electrically-charged particles exhibit complicated rheology due to the presence of electric charge on the surface of the particles.When electrically charged particles are present in an electrolyte, they attract a cloud of counter-ions, as shown schematically in Figure 2. The charge on the particle surface shown in the figure is negative and therefore, the particle attracts positively charged ions from the dispersion medium.A monolayer of immobile positively-charged ions is formed just next to the negatively-charged particle surface.This layer is called Stern layer.Outside the Stern layer, the ionic cloud consists of mobile co-ions and counter-ions and the concentration of counter-ions is higher than that of co-ions.The mobile ions make up a diffuse outer layer referred to as Gouy or Debye-Huckel layer.The concentration of counter-ions in the diffuse layer gradually decreases toward the bulk concentration as one moves away from the particle.The charge cloud surrounding the particle surface is often referred to an electrical double layer.
charged ions is formed just next to the negatively-charged particle surface.This layer is called Stern layer.Outside the Stern layer, the ionic cloud consists of mobile co-ions and counter-ions and the concentration of counter-ions is higher than that of co-ions.The mobile ions make up a diffuse outer layer referred to as Gouy or Debye-Huckel layer.The concentration of counter-ions in the diffuse layer gradually decreases toward the bulk concentration as one moves away from the particle.The charge cloud surrounding the particle surface is often referred to an electrical double layer.
The suspensions of electrically charged particles are more viscous in comparison with the suspensions of neutral particles.The motion of fluid surrounding the charged particle distorts the electric charge cloud from its equilibrium state, resulting in additional stresses referred to as Maxwell stresses.The Maxwell stresses produced by the electrostatic attraction between the particle's surface charge and the surrounding charge cloud of counter-ions hinder the motion of fluid inside the double layer and consequently an extra dissipation of energy takes place.This effect whereby additional stresses are produced and extra dissipation of energy takes place in a particulate fluid due to distortion of electrical double layers around the particles caused by fluid motion is called the "primary electroviscous effect".The rheology of suspensions of charged particles is investigated by a number of researchers [18][19][20][21][22][23][24][25][26][27][28].The rheology of such suspensions is governed by a complex interplay of electrical, thermal, and viscous forces.The key factors to be considered are: the thickness of the electrical double layer with respect to the particle size, the distortion of the double layer from equilibrium due to fluid motion, and the resistance to fluid flow around the particle caused by electrical forces.For example, the viscosity of a suspension is expected to increase with the increases in the double layer thickness and electrical stress.Furthermore, shear-thinning is expected due to distortion of charge cloud around the particles caused by fluid motion.The suspensions of electrically charged particles are more viscous in comparison with the suspensions of neutral particles.The motion of fluid surrounding the charged particle distorts the electric charge cloud from its equilibrium state, resulting in additional stresses referred to as Maxwell stresses.The Maxwell stresses produced by the electrostatic attraction between the particle's surface charge and the surrounding charge cloud of counter-ions hinder the motion of fluid inside the double layer and consequently an extra dissipation of energy takes place.This effect whereby additional stresses are produced and extra dissipation of energy takes place in a particulate fluid due to distortion of electrical double layers around the particles caused by fluid motion is called the "primary electroviscous effect".
The rheology of suspensions of charged particles is investigated by a number of researchers [18][19][20][21][22][23][24][25][26][27][28].The rheology of such suspensions is governed by a complex interplay of electrical, thermal, and viscous forces.The key factors to be considered are: the thickness of the electrical double layer with respect to the particle size, the distortion of the double layer from equilibrium due to fluid motion, and the resistance to fluid flow around the particle caused by electrical forces.For example, the viscosity of a suspension is expected to increase with the increases in the double layer thickness and electrical stress.Furthermore, shear-thinning is expected due to distortion of charge cloud around the particles caused by fluid motion.
The dimensional analysis indicates that four dimensionless groups are important, namely the Hartmann number (Ha), the ion Peclet number (Pe), the relative thickness of the electrical double layer (1/Rκ), and the dimensionless surface potential ( ψ 0 ).The definitions of these dimensionless groups are as follows: .
Rκ = R e 2 I/εkT = particle − size double − layer thickness (24) where ψ 0 is surface potential (zeta potential), I is the bulk ionic strength defined as i is the number density of ionic species i far away from the particle, N is the total number of ionic species present in the fluid, D is the ion diffusivity, ε is the permittivity of the fluid, k is the Boltzmann constant, T is the absolute temperature, R is the particle radius, η c is the viscosity of fluid, e is the elementary electronic charge, .γ is the shear rate, and κ −1 is the Debye length defined as: The electric Hartmann number Ha is a measure of the ratio of electrical force to viscous force acting on the fluid.It represents rigidity of the charge cloud.Note that the diffusion of ions counteracts the Maxwell stresses and therefore, the Hartmann number includes the ratio of diffusion time to convection time of ions (see Equation (22)).The modification of the flow field around the particles caused by the electrostatic body force is small when Ha << 1.When Ha >> 1, the fluid adjacent to the particle experiences strong electrical forces and therefore, remains almost attached to the particle.
The ion Peclet number defined in Equation (23) measures the extent to which the motion of the fluid relative to the particle disturbs the charge cloud.For small Pe(Pe << 1) the ionic atmosphere is distorted only slightly from its equilibrium shape because the diffusion of ions is strong as compared to convection of ions.The dimensionless group Rκ defined in Equation ( 24) is the ratio of particle radius to Debye length κ −1 .As the double layer extends a distance of κ −1 from the surface of the particle, (Rκ) −1 is a measure of the thickness of the electric double layer, relative to the particle size.If Rκ >> 1, the electric double layer is thin and if Rκ << 1, the electric double layer is thick.
The dimensionless surface potential ( ψ 0 ) defined in Equation ( 25) is the ratio of electrical energy of ions to their thermal energy.When ψ 0 << 1, the number density of ions in the charge cloud differs only slightly from its value far away from the particle.In other words, the ion distribution in a fluid is affected only to a small extent near the particle surface, when a charged particle is introduced to the fluid.
When the distortion of charge clouds around the particles caused by fluid motion is small, that is, Pe << 1, the constitutive equation for a dilute dispersion of electrically charged spherical particles can be expressed as: where p is the primary electroviscous coefficient which is, in general, a function of surface potential ψ 0 , Hartmann number Ha, and double layer thickness (Rκ).The dispersion viscosity can be obtained from Equation (27) as: Equation ( 28) could be re-cast in terms of the intrinsic viscosity [η] as follows: where The primary electroviscous coefficient p at low surface potential ( ψ 0 << 1) and low Hartmann number (Ha << 1) is given as [24]: where f (Rκ) is a power series function of Rκ with the two limiting forms as follows: The form of f (Rκ) given in Equation ( 32) is valid for small Rκ(Rκ << 1), that is, thick double layers.The form of f (Rκ) expressed in Equation ( 33) is valid for large Rκ(Rκ >> 1), that is, thin double layers.For thin electrical layers (Rκ >> 1), the primary electroviscous coefficient p is given as: whereas for thick electrical layers (Rκ << 1), the coefficient p reduces to: Ohshima [29] has developed the following analytical expression for the primary electroviscous coefficient p, valid for arbitrary surface potentials (no restriction on ψ 0 ), large Rκ(Rκ >> 1), and Z-Z symmetrical type electrolyte where Z is the valence of the ion [29]: where In these equations, λ ± is the drag coefficient of ions, m ± is the scaled ionic drag coefficient, F is the Dukhin number [30], and ψ 0 = Zeψ 0 /kT.Using the relation λ = kT/D (where D is the ion diffusivity), it can be readily shown that: Ha For aqueous KCl solution at 25 • C, m − = 0.169 and m + = 0.176 [30].For low surface potential and large Rκ, Equation (36) reduces to: Assuming m − = m + = m, Equation (40) in conjunction with Equation (39) leads to the expression for p given by Equation (34).
Figure 3 shows the predictions of Equation (36) for different values of the double layer thickness Rκ.The electrolyte is KCl with m − = 0.169 and m + = 0.176.The primary electroviscous coefficient p decreases with the increase in Rκ (decrease in double layer thickness) for any given surface potential.Interestingly, the plot of p versus ψ 0 for any given Rκ exhibits a maximum at some intermediate value of ψ 0 due to relaxation effect [29].At high ψ 0 , p reaches a nonzero limiting value of 72m + (ln2/Rκ) 2 .It should also be noted that the value of ψ 0 where relaxation (peak) is observed increases slightly with the increase in Rκ.The electroviscous coefficient p could be determined experimentally using viscosity measurements of dilute suspensions.Upon re-arrangement of Equation ( 29), one can write: From the slope of The electroviscous coefficient p could be determined experimentally using viscosity measurements of dilute suspensions.Upon re-arrangement of Equation (29), one can write: From the slope of 1/η r versus φ data in the limit φ → 0 , one can determine the value of [η] (and hence the primary electroviscous coefficient p, see Equation (30)) experimentally.As an example, Figure 4 shows the plots of 1/η r versus φ data for polystyrene latex suspensions consisting of negatively charged particles of diameter 250 nm [27]   Figure 5 shows the variation of primary electroviscous coefficient p with the electrolyte KCl concentration for two negatively charged monodisperse polystyrene latex suspensions [27].The surface charge on the particles of more viscous latex suspension (top curve) was −6.85 μC/cm 2 and the particle diameter was 106 nm.The surface charge on the particles of less viscous latex suspension (bottom curve) was −0.20 μC/cm 2 and the particle diameter was 210 nm.With the increase in electrolyte concentration, the electrical double layer thickness decreases (  R increases) and therefore, the coefficient p decreases.For the same electrolyte concentration, the suspension of weakly-charged particles gives lower p values, as expected.Furthermore, the  R value of the suspension of weakly-charged particles is also larger (hence thinner double layer) due to large particle size.Figure 5 shows the variation of primary electroviscous coefficient p with the electrolyte KCl concentration for two negatively charged monodisperse polystyrene latex suspensions [27].The surface charge on the particles of more viscous latex suspension (top curve) was −6.85 µC/cm 2 and the particle diameter was 106 nm.The surface charge on the particles of less viscous latex suspension (bottom curve) was −0.20 µC/cm 2 and the particle diameter was 210 nm.With the increase in electrolyte concentration, the electrical double layer thickness decreases (Rκ increases) and therefore, the coefficient p decreases.For the same electrolyte concentration, the suspension of weakly-charged particles gives lower p values, as expected.Furthermore, the Rκ value of the suspension of weakly-charged particles is also larger (hence thinner double layer) due to large particle size.Thus far, our discussion about the primary electroviscous effect was largely restricted to small Pe ( Pe << 1).The ionic atmosphere around the particles was distorted only slightly from its equilibrium shape due to fluid motion.Russel [19] extended the work to arbitrary flow strengths (arbitrary Pe ) and developed the following constitutive equation for dilute dispersions of charged spherical particles valid under conditions of thin double layer (  R >> 1), small Hartmann number and low surface potential 0  << 1.The full constitutive equation developed by Russel [19] for suspensions of electrically charged spheres is as follows: where  is the deviatoric stress tensor and

Dt D/
is the Jaumann derivative.This constitutive equation predicts non-Newtonian shear-thinning behaviour of suspensions of electrically charged particles.Equation (42) gives the following expression for the suspension viscosity: Thus the viscosity of dilute suspensions of electrically charged particles decreases with the increase in ion Peclet number Pe .As an example, Figures 6 and 7  Thus far, our discussion about the primary electroviscous effect was largely restricted to small Pe(Pe << 1).The ionic atmosphere around the particles was distorted only slightly from its equilibrium shape due to fluid motion.Russel [19] extended the work to arbitrary flow strengths (arbitrary Pe) and developed the following constitutive equation for dilute dispersions of charged spherical particles valid under conditions of thin double layer (Rκ >> 1), small Hartmann number and low surface potential ψ 0 << 1.The full constitutive equation developed by Russel [19] for suspensions of electrically charged spheres is as follows: where τ is the deviatoric stress tensor and D/Dt is the Jaumann derivative.This constitutive equation predicts non-Newtonian shear-thinning behaviour of suspensions of electrically charged particles.Equation (42) gives the following expression for the suspension viscosity: Thus the viscosity of dilute suspensions of electrically charged particles decreases with the increase in ion Peclet number Pe.As an example, Figures 6 and 7 show the variations of intrinsic viscosity of suspensions of charged particles with ion Peclet number Pe.The plots are generated from Equation (43).As expected, the suspension exhibits shear-thinning in that the intrinsic viscosity decreases with the increase in the Peclet number.The [η] versus Pe plots exhibit with two distinct plateau regions: upper Newtonian plateau with high [η] value when Pe < 1 and lower Newtonian plateau with low [η] value when Pe > 1.With the increase in the electrical double layer thickness, that is, decrease in Rk, the intrinsic viscosity in the upper Newtonian region (Pe < 1) increases at a fixed value of Hartmann number Ha (see Figure 6).With the increase in the Hartmann number Ha, the intrinsic viscosity in the upper Newtonian region (Pe < 1) increases at a fixed value of the electrical double layer thickness (see Figure 7).6).With the increase in the Hartmann number Ha , the intrinsic viscosity in the upper Newtonian region ( 1  Pe ) increases at a fixed value of the electrical double layer thickness (see Figure 7).Equation ( 42) also predicts that the dilute suspensions of electrically charged particles are viscoelastic in nature.It leads to the following expressions for the normal-stress differences in steady shearing flow: Equation ( 42) also predicts that the dilute suspensions of electrically charged particles are viscoelastic in nature.It leads to the following expressions for the normal-stress differences in steady shearing flow: Fluids 2016, 1, 40 14 of 52  Figure 8 shows the plots of first normal-stress difference as N 1 / ε ψ 2 0 /R 2 versus ion Peclet number Pe for two different values of particle volume fraction.The first normal-stress difference increases initially with the increase in Pe and then levels off at high Pe.It also increases with the increase in particle volume fraction φ.The second normal-stress difference is plotted as −N 2 / εψ 2 0 /R 2 versus Pe in Figure 9.The second normal stress difference behaves in a manner similar to the first normal-stress difference.34) and ( 47) are compared for thin double layers (  R >> 1) and the predictions of Equations ( 35) and ( 48) are compared for thick double layers (  R << 1).Clearly, the bulk primary electroviscous coefficient is much higher than the shear primary electroviscous coefficient for the same electrical double layer thickness.34) and ( 47) are compared for thin double layers (  R >> 1) and the predictions of Equations ( 35) and ( 48) are compared for thick double layers (  R << 1).Clearly, the bulk primary electroviscous coefficient is much higher than the shear primary electroviscous coefficient for the same electrical double layer thickness.The constitutive equation developed by of Russel (Equation ( 42)) is valid for thin double layers (Rκ >> 1), small Hartmann number, and low surface potential.Lever [22] has extended the analysis of Russel to thick double layers (Rκ << 1) and conducted numerical calculations.The Hartmann number and surface potential were assumed to be small.His analysis shows that the suspension behaves as a non-Newtonian fluid with non-zero normal-stress differences and shear-thinning viscosity at high Peclet numbers when the charge cloud is deformed from equilibrium condition to a large extent.
Khair and Star [31] have considered the influence of particle surface charge on the bulk viscosity of dilute suspensions under the conditions of low surface potential, small Hartmann number, and gentle flow (small Peclet number).The bulk viscosity of a dilute suspension of electrically-charged spherical particles can be expressed as: Fluids 2016, 1, 40 where p k is the primary bulk electroviscous coefficient.For thin electrical layers (Rκ>> 1), the primary bulk electroviscous coefficient p k is given as: whereas for thick electrical layers (Rκ << 1), the coefficient p k reduces to: Figure 10 compares the bulk and shear primary electroviscous coefficients at a small Hartmann number of 0.1.The predictions of Equations ( 34) and ( 47) are compared for thin double layers (Rκ >> 1) and the predictions of Equations ( 35) and ( 48) are compared for thick double layers (Rκ << 1).Clearly, the bulk primary electroviscous coefficient is much higher than the shear primary electroviscous coefficient for the same electrical double layer thickness.

Rheology of Suspensions of Rigid-Porous Spherical Particles
The rheology of dilute suspensions of porous rigid spherical particles is strongly dependent on the particle permeability.Zackrisson and Bergenholtz [32] have shown that the dipole strength of a porous rigid spherical particle immersed in an infinite matrix of incompressible Newtonian fluid is given as: where  and    H are defined as: In Equation ( 50),  is the permeability of particles.For impermeable particles, and consequently, the dipole strength given by Equation ( 49) reduces to that of solid spherical particles given by Equation (12).For completely permeable particles,    , 0   , and consequently, the dipole strength given by Equation ( 44) reduces to zero.The full constitutive equation for dilute suspensions of electrically-neutral porous rigid spherical particles can be expressed as follows:

Rheology of Suspensions of Rigid-Porous Spherical Particles
The rheology of dilute suspensions of porous rigid spherical particles is strongly dependent on the particle permeability.Zackrisson and Bergenholtz [32] have shown that the dipole strength of a porous rigid spherical particle immersed in an infinite matrix of incompressible Newtonian fluid is given as: where β and H (β) are defined as: In Equation ( 50), α is the permeability of particles.For impermeable particles, α → 0 , β → ∞ , and consequently, the dipole strength given by Equation ( 49) reduces to that of solid spherical particles given by Equation (12).For completely permeable particles, α → ∞ , β → 0 , and consequently, the dipole strength given by Equation ( 44) reduces to zero.The full constitutive equation for dilute suspensions of electrically-neutral porous rigid spherical particles can be expressed as follows: where σ is the bulk stress tensor of the suspension and E is the imposed or bulk rate of strain tensor on the suspension.According to Equation ( 52), the suspension is a Newtonian fluid with effective viscosity of η given as: Equation ( 53) gives the following expression for the intrinsic viscosity of suspensions of electrically-neutral, porous, rigid and spherical particles: Figure 11 shows the plot of intrinsic viscosity [η] as a function of dimensionless permeability, defined as α/R 2 (note that α/R 2 = 1/β 2 ).The intrinsic viscosity decreases from a value of 2.5 corresponding to completely impermeable spheres ( α/R 2 → 0 ) to a value of zero, corresponding to completely permeable spheres ( α/R 2 → ∞ ).Natraj and Chen [33] have extended the analysis of the rheology of dilute suspensions of electrically-neutral porous spheres to electrically-charged spheres taking into consideration the primary electroviscous effect.According to the analysis presented by Natraj and Chen, the distortion of the charge cloud around the particles increases with the increase in the permeability of the particles leading to an increase in the intrinsic viscosity of suspension relative to intrinsic viscosity of suspension of uncharged porous particles.As an example, Figure 12 shows the plots of the ratio [η]/[η] neutral as a function of dimensionless permeability of particles for two different values of Rk (ratio of particle radius to electrical double layer thickness).For any given plot, the Hartmann number is kept constant.Clearly the suspension of charged porous particles is more viscous as compared with suspension of uncharged porous particles for any given particle permeability.The enhancement of suspension viscosity due to primary electroviscous effect becomes more intense with the increase in particle permeability.Equation ( 53) gives the following expression for the intrinsic viscosity of suspensions of electrically-neutral, porous, rigid and spherical particles: Figure 11 shows the plot of intrinsic viscosity ] [ as a function of dimensionless permeability, defined as ).The intrinsic viscosity decreases from a value of 2.5 corresponding to completely impermeable spheres ( ) to a value of zero, corresponding to completely permeable spheres ( ). Natraj and Chen [33] have extended the analysis of the rheology of dilute suspensions of electrically-neutral porous spheres to electrically-charged spheres taking into consideration the primary electroviscous effect.According to the analysis presented by Natraj and Chen, the distortion of the charge cloud around the particles increases with the increase in the permeability of the particles leading to an increase in the intrinsic viscosity of suspension relative to intrinsic viscosity of suspension of uncharged porous particles.As an example, Figure 12 shows the plots of the ratio as a function of dimensionless permeability of particles for two different values of Rk (ratio of particle radius to electrical double layer thickness).For any given plot, the Hartmann number is kept constant.Clearly the suspension of charged porous particles is more viscous as compared with suspension of uncharged porous particles for any given particle permeability.The enhancement of suspension viscosity due to primary electroviscous effect becomes more intense with the increase in particle permeability.

Rheology of Suspensions of Non-Rigid (Soft) Solid Particles
Particulate fluids consisting of deformable solid particles exhibit non-Newtonian behavior [9,34,35].In shear flow, they exhibit shear-thinning and normal-stress effects.In uniaxial elongational flow, the strain-thickening effect is observed.In simple shear flow, the initially spherical particle deforms into an ellipsoid of fixed dimensions and orientation at a given shear rate.The material inside the ellipsoidal particle undergoes continuous deformation and rotation (tank-treading motion) The rheology of suspensions of particles composed of deformable purely elastic material (Hookean solid), with elastic modulus of p G , depends on the following dimensionless group: Se is the ratio of viscous stress which tends to elongate and orient the particles in the direction of flow to elastic stress which tends to recover the original shape of the particles.When Se << 1, the suspension is more viscous as the elastic force dominates and the particles are nearly spherical.When Se >>1, the suspension becomes less viscous as viscous stresses now dominate and the particles become elongated and oriented in the direction of flow.
Goddard and Miller [34] calculated the dipole strength of soft purely elastic particles immersed in an infinite matrix of incompressible Newtonian fluid for weakly time-dependent flows.Their expression for o S could be expressed as follows: where d S refers to the symmetric and traceless part of the indicated tensor.From Equations ( 8) and (56), it follows that: Influence of particle permeability on the intrinsic viscosity of suspension of electrically-charged porous particles.

Rheology of Suspensions of Non-Rigid (Soft) Solid Particles
Particulate fluids consisting of deformable solid particles exhibit non-Newtonian behavior [9,34,35].In shear flow, they exhibit shear-thinning and normal-stress effects.In uniaxial elongational flow, the strain-thickening effect is observed.In simple shear flow, the initially spherical particle deforms into an ellipsoid of fixed dimensions and orientation at a given shear rate.The material inside the ellipsoidal particle undergoes continuous deformation and rotation (tank-treading motion) The rheology of suspensions of particles composed of deformable purely elastic material (Hookean solid), with elastic modulus of G p , depends on the following dimensionless group: Se is the ratio of viscous stress which tends to elongate and orient the particles in the direction of flow to elastic stress which tends to recover the original shape of the particles.When Se << 1, the suspension is more viscous as the elastic force dominates and the particles are nearly spherical.When Se >> 1, the suspension becomes less viscous as viscous stresses now dominate and the particles become elongated and oriented in the direction of flow.
Goddard and Miller [34] calculated the dipole strength of soft purely elastic particles immersed in an infinite matrix of incompressible Newtonian fluid for weakly time-dependent flows.Their expression for S o could be expressed as follows: where S d refers to the symmetric and traceless part of the indicated tensor.From Equations ( 8) and ( 56), it follows that: Note that for a dilute suspension, the rate of strain tensor E ∞ far away from the particle can be equated to the bulk or imposed rate of strain tensor E on the suspension.Equation ( 57) could be further approximated as [9]: Equation ( 58) predicts non-Newtonian behavior of dilute suspensions of deformable purely elastic solid particles.It gives the following expressions for the relative shear viscosity and reduced normal stresses when the suspension is subjected to steady shearing flow: 59) 60) where Se is the dimensionless group defined in Equation ( 55), N 1 is the first normal stress difference and N 2 is the second normal stress difference.From Equation ( 59), it follows that the intrinsic viscosity of a suspension of soft elastic solid particles is: It may be of interest to note that Roscoe [35] presented a similar set of equations for dilute suspensions of viscoelastic solid-like particles (Kelvin-Voigt type solid particles).In the case of purely elastic solid particles, the equations presented by Roscoe [35] reduce to the equations presented here (Equations ( 59)-( 62)).
Figure 13 shows the plot of intrinsic viscosity [η] as a function of Se.The intrinsic viscosity decreases with the increase in Se indicating that suspension of soft elastic particles is shear-thinning in nature.It should also be noted that at high Se(Se > 1), the intrinsic viscosity [η] becomes negative indicating that the viscosity of suspensions of soft elastic particles decreases with the increase in particle concentration at high Se.According to Equation ( 54), the relative viscosity η r of suspension becomes equal to (1 − 5φ/3) at high Se.The reduced normal stress differences are plotted in Figures 14  and 15 as N 1r /φ (Figure 14) and −N 2r /φ (Figure 15) versus Se.The reduced normal stress differences increase initially with the increase in Se, reach some maximum values, and then decrease with further increase in Se.The constitutive equation expressed in the form of Equation ( 58) also predicts strain-thickening behavior in uniaxial elongational flow, that is, the elongational viscosity ( E  ) increases with the increase in the elongation rate ( E  ) as indicated by the following equation: The constitutive equation expressed in the form of Equation ( 58) also predicts strain-thickening behavior in uniaxial elongational flow, that is, the elongational viscosity (η E ) increases with the increase in the elongation rate ( .γ E ) as indicated by the following equation: Finally, it should be realized that the Goddard and Miller analysis [34] and therefore the Equations ( 58)-( 63) assume infinitesimal deformation of particles.Thus the plots shown in Figures 13-15 are strictly valid for Se << 1.Little or no experimental data are available to verify these theoretical predictions.However, a number of interesting analytical and numerical studies have been published recently extending the works of Goddard and Miller [34] and Roscoe [35] to arbitrarily large deformations [36][37][38][39][40][41].Furthermore, Gao et al. [38] have considered the effect of the initial particle shape of elastic particles (usually assumed to be spherical) on the motion and rheology of elastic particles in a shear flow.While the initially-spherical elastic particles exhibit the well-known "tank-treading" motion in shear flow, the initially-ellipsoidal elastic particles exhibit so-called "trembling motion" in that they tend to oscillate without making a full rotation.As the initial shape of the particle becomes more elongated, the elastic particles undergo "tumbling motion" similar to that of rigid ellipsoidal particles in a shear flow.

Rheology of Emulsions
Emulsions are dispersed systems consisting of droplets of one liquid distributed throughout another immiscible liquid.Emulsions exhibit non-Newtonian rheology due to deformation and orientation of droplets in a flow field [9,[42][43][44][45].The shape and orientation of droplet with respect to the flow direction is governed by two dimensionless groups: the viscosity ratio λ, defined as the ratio of droplet fluid viscosity (η d ) to continuous-phase viscosity (η c ), and the capillary number Ca, defined as: where .
γ is the shear rate, R is the radius of the un-deformed droplet, and γ is the interfacial tension between the two immiscible liquids.There is competition between the hydrodynamic stress ∼ η c .γ which tends to stretch the droplet and the interfacial stress (∼ γ/R) which tends to contract the droplet into a sphere.With the increase in the capillary number, the droplet becomes more elongated and oriented with the flow direction.Due to this elongation and orientation of droplets emulsions exhibit shear-thinning behavior.Other non-Newtonian characteristics exhibited by emulsions are (i) non-zero normal stresses in shear flow and (ii) strain-thickening in extensional flows.
Assuming droplets to remain nearly spherical in shape, that is, Ca → 0 , Taylor [46] developed the following expression for the dipole strength of a droplet immersed in an infinite matrix of incompressible Newtonian fluid: where λ is the ratio of dispersed-phase viscosity to continuous-phase viscosity.The constitutive equation for the dilute emulsion in the limit Ca → 0 can be obtained from Equations ( 8) and (65) as follows: Thus the effective viscosity of a dilute emulsion in the limit Ca → 0 turns out to be Equation ( 67), referred to as the Taylor equation, reduces to the Einstein equation, Equation (15), in the limit λ → ∞ .From Equation (67), the intrinsic viscosity of an emulsion is as follows: Figure 16 shows the plots of relative viscosity η r versus droplet concentration φ for dilute emulsions for different values of the viscosity ratio λ.The plots are generated using the Taylor equation, Equation (67).With the increase in λ, the slope of the plot, that is [η], increases.When λ → ∞ , the slope [η] is 2.5 and the droplets behave as rigid particles.When λ → 0 , the slope [η] is 1.0 and the droplets behave as bubbles.Thus the relative viscosity of emulsions can be manipulated by varying the viscosity ratio while keeping all the other compositional variables constant.One possible way to manipulate the viscosity ratio, and hence the relative viscosity, in practical emulsions is to vary the temperature of the emulsion [47] provided that the temperature-dependence of viscosity is different for the two liquids involved.Figure 17 shows comparison between experimental viscosity data and theoretical predictions.
The emulsions tested are toluene-in-glycerol [48] with viscosity ratio  of 8.66 × 10 The constitutive equation of a dilute emulsion expressed in the form of Equation ( 66) does not predict non-Newtonian rheology of emulsions as deformation of droplets is neglected and 0  Ca .Schowalter et al. [42] and Frankel and Acrivos [43] extended the work of Taylor [46]   Figure 17 shows comparison between experimental viscosity data and theoretical predictions.The emulsions tested are toluene-in-glycerol [48] with viscosity ratio λ of 8.66 × 10 −4 .The emulsions are produced by mixing the two liquids in an agitated vessel at different agitator speeds.The intrinsic viscosity (Equation ( 68)) of these emulsions is 1.0013.The experimental data show good agreement with the Taylor viscosity equation when φ ≤ 0.05.At higher concentrations of dispersed-phase, the Taylor equation under predicts the emulsion viscosity.According to Frankel and Acrivos [43], Equation ( 70) could be recast as: where  is the deviatoric stress tensor and o  is the characteristic time constant of the emulsion given as: Equation ( 71) predicts shear-thinning in emulsions in steady-shear flows and leads to the following expression for the shear viscosity [9] of dilute emulsions: Based on Equation ( 73), the intrinsic viscosity of an emulsion can be expressed as: where and Ca is the capillary number (defined in Equation ( 64)).
Thus intrinsic viscosity of an emulsion depends on viscosity ratio  and capillary number Ca .
In the limit 0  Ca , Equation (74) reduces to Equation (68). Figure 18 shows the plots of intrinsic 0.9 The constitutive equation of a dilute emulsion expressed in the form of Equation ( 66) does not predict non-Newtonian rheology of emulsions as deformation of droplets is neglected and Ca → 0 .Schowalter et al. [42] and Frankel and Acrivos [43] extended the work of Taylor [46] by considering the first-order deformation of droplets and developed the following equation for the dipole strength of a droplet immersed in an infinite matrix of incompressible Newtonian fluid: From Equations ( 8) and ( 69), it follows that: where E ∞ has been equated to the imposed rate of strain tensor E on the emulsion.According to Frankel and Acrivos [43], Equation (70) could be recast as: Fluids 2016, 1, 40 where τ is the deviatoric stress tensor and τ o is the characteristic time constant of the emulsion given as: Equation ( 71) predicts shear-thinning in emulsions in steady-shear flows and leads to the following expression for the shear viscosity [9] of dilute emulsions: Based on Equation ( 73), the intrinsic viscosity of an emulsion can be expressed as: where h = (19λ + 16)(2λ + 3)/[40 (1 + λ)] and Ca is the capillary number (defined in Equation ( 64)).
Thus intrinsic viscosity of an emulsion depends on viscosity ratio λ and capillary number Ca.In the limit Ca → 0 , Equation ( 74) reduces to Equation (68). Figure 18 shows the plots of intrinsic viscosity generated from Equation (74).At low values of capillary number (Ca < 0.1), the intrinsic viscosity is always greater than or equal to 1, regardless of the viscosity ratio, indicating that the addition of droplets always increases the viscosity of emulsion due to interfacial tension effect.At high capillary numbers where the interfacial tension effect is negligible, the intrinsic viscosity is non-zero positive only if the viscosity ratio is larger than 1.When λ = 1, the intrinsic viscosity is zero at high capillary numbers indicating that the addition of droplets has no effect on the viscosity of emulsion.When λ < 1, the intrinsic viscosity is negative at high capillary numbers indicating that the addition of droplets results in a decrease in the viscosity of emulsion.
Fluids 2016, 1, 40 25 of 52 viscosity generated from Equation (74).At low values of capillary number ( Ca < 0.1), the intrinsic viscosity is always greater than or equal to 1, regardless of the viscosity ratio, indicating that the addition of droplets always increases the viscosity of emulsion due to interfacial tension effect.At high capillary numbers where the interfacial tension effect is negligible, the intrinsic viscosity is non-zero positive only if the viscosity ratio is larger than 1.When 1   , the intrinsic viscosity is zero at high capillary numbers indicating that the addition of droplets has no effect on the viscosity of emulsion.When 1   , the intrinsic viscosity is negative at high capillary numbers indicating that the addition of droplets results in a decrease in the viscosity of emulsion.versus Ca , keeping 1   .The normal stress differences increase initially with the increase in Ca due to interfacial-tension related elasticity of droplets.After reaching a certain maximum value, the normal stress differences fall off sharply with further increase in Ca .The decrease in the magnitude of the normal stress differences at high Ca is expected as viscous stresses dominate over the restoring interfacial stresses at high The Frankel and Acrivos equation, Equation (71), also predicts non-zero normal stress differences in shear flow of dilute emulsions.Equation (71) leads to the following expressions for the normal stress differences: Fluids The plots of normal stress differences for emulsions are shown in Figures 19 and 20.In Figure 19, the first normal stress difference is plotted as N 1r /φ versus Ca, keeping λ = 1.In Figure 20, the second normal stress difference is plotted as −N 2r /φ versus Ca, keeping λ = 1.The normal stress differences increase initially with the increase in Ca due to interfacial-tension related elasticity of droplets.After reaching a certain maximum value, the normal stress differences fall off sharply with further increase in Ca.The decrease in the magnitude of the normal stress differences at high Ca is expected as viscous stresses dominate over the restoring interfacial stresses at high Ca.According to the Frankel and Acrivos equation, Equation ( 71), the emulsions behave as strainthickening in uniaxial elongation flows as indicated by the following expression of elongation viscosity E  obtained from Equation (71): Finally it should be pointed out that the Frankel and Acrivos analysis and therefore the Equations ( 69)-(77) assume infinitesimal (first order) deformation of droplets.Thus the plots shown in Figures 18-20 are strictly valid for small deformation of droplets corresponding to 1  Ca .Wetzel and Tucker [39] have extended the small deformation analysis to arbitrary deformations for dispersions with unequal viscosities and zero interfacial tension.According to the Frankel and Acrivos equation, Equation ( 71), the emulsions behave as strainthickening in uniaxial elongation flows as indicated by the following expression of elongation viscosity E  obtained from Equation (71): Finally it should be pointed out that the Frankel and Acrivos analysis and therefore the Equations ( 69)-(77) assume infinitesimal (first order) deformation of droplets.Thus the plots shown in Figures 18-20 are strictly valid for small deformation of droplets corresponding to 1  Ca .Wetzel and Tucker [39] have extended the small deformation analysis to arbitrary deformations for dispersions with unequal viscosities and zero interfacial tension.According to the Frankel and Acrivos equation, Equation ( 71), the emulsions behave as strain-thickening in uniaxial elongation flows as indicated by the following expression of elongation viscosity η E obtained from Equation (71): Finally it should be pointed out that the Frankel and Acrivos analysis and therefore the Equations ( 69)-( 77) assume infinitesimal (first order) deformation of droplets.Thus the plots shown in Figures 18-20 are strictly valid for small deformation of droplets corresponding to Ca < 1. Wetzel and Tucker [39] have extended the small deformation analysis to arbitrary deformations for dispersions with unequal viscosities and zero interfacial tension.

Influence of Electric Charge Present on the Surface of Emulsion Droplets
Ohshima [29] has studied the influence of droplet surface charge on the rheology of dilute emulsions.The relative viscosity of dilute emulsions of electrically charged droplets at low capillary numbers can be expressed in the following modified-form of the Taylor viscosity equation (Equation ( 67)): where p is the primary electroviscous coefficient.Ohshima [29] has developed a general expression for p valid under the following conditions: arbitrary surface potential, arbitrary electric double layer thickness Rκ, weak flow (small Peclet number) such that the charge cloud around the droplet is disturbed from equilibrium state only slightly, and small capillary number so that the droplets remain spherical.For large Rκ ( Rκ → ∞ ) and arbitrary surface potential, the expression for p of emulsions reduces to [29]: where m ± is the scaled ionic drag coefficient, defined in Equation (37).In the limit ψ 0 → 0 , Equation ( 79) gives p = 0 as expected for the uncharged droplets.For very thin double layer ( Rκ → ∞ ) and high surface potential ψ 0 → ∞ , p approaches a nonzero limiting value given by: If we substitute this expression of into Equation ( 78), the Einstein viscosity equation, Equation (18), for dilute suspensions of uncharged rigid-solid spheres is recovered.This implies that emulsion droplets with thin double layers and very high surface potentials behave like uncharged rigid-solid spheres regardless of the viscosity ratio λ.This phenomenon is referred to as "solidification effect" in the literature [29].
Figure 21 shows the plots of p versus ψ 0 for different values of viscosity ratio λ.The plots are generated using Equation (79).The matrix is assumed to be an aqueous KCl solution.As expected, p increases initially with the increase in surface potential and then levels off at high surface potentials.With the increase in the viscosity ratio, the primary electroviscous coefficient p of emulsion decreases.

Influence of Surfactant on Emulsion Rheology
In the preceding discussion, it was observed that dilute emulsions of electrically-neutral droplets exhibit Newtonian behavior when As λ → ∞ , p → 0 , indicating that the influence of interface becomes less and less important with the increase in the viscosity ratio.

Influence of Surfactant on Emulsion Rheology
In the preceding discussion, it was observed that dilute emulsions of electrically-neutral droplets exhibit Newtonian behavior when Ca → 0 and the droplets are spherical.The emulsions exhibit shear-thinning and elastic behavior (non-zero normal stress differences) only when Ca > 0 and the droplets undergo significant deformation.It is interesting to note that dilute emulsions can exhibit non-Newtonian rheology even in the absence of any deformation, that is, when Ca → 0 and the droplets are spherical [49][50][51].The additional mechanism of non-Newtonian behavior in emulsions comes into play when surfactant molecules are present at the surface of droplets.The flow of continuous-phase fluid around the droplet results in non-uniform concentration of surfactant at the droplet surface.The resulting inhomogeneity in surfactant concentration at the droplet surface generates Marangoni stress due to interfacial tension gradient.The Marangoni stress results in a jump in tangential viscous traction across the droplet interface.Due to inhibition of the transmission of tangential stresses from the continuous-phase fluid to the droplet fluid, the droplets tend to behave more like solid particles.
The dimensionless group that determines the rheology of dilute emulsion of surfactant covered droplets when Ca → 0 is the Marangoni number Ma. Blawzdziewicz et al. [49] define Ma as: where ∆γ = γ o − γ e (γ o is interfacial tension in the absence of surfactant and γ e is equilibrium interfacial tension in the presence of surfactant at the droplet surface).At low shear rates (high Ma), the Marangoni stress dominates over the viscous stress and therefore, the emulsion droplets behave more like rigid particles.Consequently the emulsion viscosity is relatively high.At high shear rates (low Ma), the viscous stress dominates over the Marangoni stress resulting in lower emulsion viscosity.The Marangoni number defined in Equation ( 81) is relevant when the surfactant is non-diffusing type.When diffusion of surfactant is important, the dimensionless group that governs the rheology of emulsions is given as follows [51]: where E G is the Gibbs elasticity of the interface and D is the effective diffusion coefficient of surfactant.
As diffusion of surfactant molecules counteracts the Marangoni-stress, the dimensionless group m_Ma (modified Marangoni number) includes the ratio of diffusion time to convection time of surfactant molecules.When the diffusion of surfactant molecules is very fast, that is D is large, m_Ma is very small and therefore, the effect of interface elasticity on emulsion rheology can be neglected.
In gentle flows (small Peclet number, small capillary number), the effective viscosity of a dilute emulsion can be expressed as [52]: where ε r is the interfacial-rigidity parameter, defined as: The interfacial-rigidity parameter ε r varies from 0 to 1.When the interface is completely rigid, ε r = 1 and Equation ( 83) becomes the Einstein viscosity equation for suspensions of rigid-solid particles (Equation ( 15)).The interface is completely mobile (ε r = 0) when λ → 0 and m_Ma → 0 , that is, the droplet has zero viscosity and the interface has negligible elasticity.The interface is completely rigid and the interfacial-rigidity parameter ε r is 1, under the following conditions: (a) m_Ma → ∞ , that is, the interface possesses high Gibbs elasticity (note that when m_Ma → ∞ , ε r is 1 regardless of the droplet to matrix viscosity ratio) and (b) λ → ∞ , that is, the viscosity ratio is very high (note that when λ → ∞ , ε r is 1 regardless of the degree of elasticity of the interface).The plots of ε r versus m_Ma for different values of the viscosity ratio λ are shown in Figure 22.In the viscosity expression given in the form of Equations ( 83) and ( 84), it is assumed that the interface dilational and shear viscosities are negligible.In the presence of surfactant or other additives at the interface, the interfacial viscosities may not be negligible.In order to account for the interfacial viscosities, the interfacial-rigidity parameter ε r needs to be modified as follows [9]: where Bo s is the shear Boussinesq number and Bo k is the dilational Boussinesq number.The Boussinesq numbers are dimensionless groups defined as: where η is and η ik are interfacial shear and dilational viscosities, respectively.Figure 23 show the plots of ε r versus λ for different values of the interfacial viscosity term 4Bo s + 6Bo k .The interfacial elasticity term m_Ma is kept zero.Clearly the interfacial rigidity parameter ε r increases with the increase in the interfacial viscosities.The effect of interfacial viscosities on ε r is more at low values of the viscosity ratio λ.At high values of λ, ε r is 1 regardless of the magnitude of the interfacial viscosities.

Rheology of Bubbly Suspensions
In gentle flows (low Reynolds number, low capillary number), the rheological constitutive equation for a dilute suspension of spherical bubbles can be expressed as: where η k is the bulk viscosity of suspension and η is the shear viscosity of suspension.expressions for η and η k for dilute suspensions of electrically-neutral bubbles, in the absence of any additives present at the interface, are given as follows: When the bubbles are electrically-charged or when additives (such as surfactant) are present at the interface, the above expressions for η and η k need to be modified.

Influence of Electric Charge Present on the Surface of Bubbles
Both shear viscosity of suspension η and bulk viscosity of suspension η k are affected by the presence of electric charge on the surface of the bubbles.For dilute suspensions of electrically-charged bubbles in incompressible Newtonian fluids, the shear and bulk viscosities are given as [31]: where p is the primary shear electroviscous coefficient and p k is the primary bulk electroviscous coefficient.Under the conditions of very thin double-layers ( Rκ → ∞ ), Z-Z symmetrical type electrolyte, and arbitrary surface potentials (no restriction on ψ 0 ), the primary shear electroviscous coefficient p is given by the following expression [29]: At high surface potentials, Equation (93) gives p of 3/2.Under the conditions of low surface potential, small Hartmann number, and small Peclet number, the primary bulk electroviscous coefficient p k is given as [31]: It is interesting to note that both p and p k approach a zero value for suspensions of rigid particles with very thin double layers ( Rκ → ∞ ).However, in the case of bubbles or droplets with very thin double layers ( Rκ → ∞ ), both p and p k approach non-zero values at high surface potentials.One possible reason for this is that the relative fluid velocity on the surface of the rigid particles is vanishingly small (no-slip condition) as compared with the relative fluid velocity on the mobile surfaces of bubbles or droplets [31].The electroviscous stress in the thin double layers is expected to be small when the fluid velocity relative to the particles is small in that region.

Influence of Surfactant Present on the Surface of Bubbles
When surfactant is present on the surface of the bubbles, the interfacial shear and dilational viscosities may be significant enough to affect the bulk rheology of bubbly suspensions.In addition to interfacial viscosities, the interfacial elasticity (Marangoni effect) may also have a significant influence.The effective shear viscosity of dilute bubbly suspensions in the presence of non-zero interfacial viscosities and elasticity is given by Equation ( 83) with the following expression for the interfacial-rigidity parameter ε r : The bulk viscosity of a dilute bubbly suspension with non-zero interfacial dilational viscosity is given as [9,53]: where R is the radius of the bubble.

Influence of Capillary Number
The preceding discussion on the rheology of dilute bubbly suspensions was restricted to weak flows characterized by a small capillary number ( Ca → 0 ).The bubbles were assumed to be spherical in shape.At high capillary numbers, the bubbles become elongated and oriented with the flow direction.Consequently, the bubbly suspensions exhibit non-Newtonian rheology including shear-thinning and non-zero normal stress differences.
The full rheological constitutive equation for dilute bubbly suspensions, when bubbles undergo first-order deformation, is given as follows: In steady shear flow, Equation (98) leads to the following expressions for relative viscosity and reduced normal-stress differences [54,55]: In uniaxial elongation flows, the bubbly suspensions behave as strain-thickening as indicated by the following expression of elongation viscosity η E obtained from Equation (98):

Rheology of Suspensions of Capsules
Capsules are droplets of liquid covered with a thin, deformable, elastic-solid membrane [9,56].The motion and deformation of capsules and the rheology of dilute suspensions of capsules has been studied extensively by Barthes-Biesel and co-workers .In shear flow, capsules behave differently from rigid-solid particles and liquid droplets.At low shear rates, capsules act as rigid-solid particles and rotate as a rigid body regardless of the viscosity of the enclosed liquid.However, at high shear rates, the capsules behave more like liquid droplets in that they deform and attain a steady-state orientation.The main difference between the behaviors of capsules and droplets at high shear rates is that the thin elastic-solid membrane covering the surface of the capsules undergoes tank-treading motion whereas no such membrane exists on the surface of the droplets.The liquid present inside the capsules also undergoes internal circulation induced by tank-treading motion of the membrane.
The rheology of suspension of capsules is strongly dependent on the instantaneous shape and orientation of the capsules (assumed spherical initially) covered with thin elastic-solid membranes.The time evolution of the capsule shape and orientation is specified by two second-order symmetric and traceless tensors A and B where A specifies the radial displacement of the material points of the membrane and B specifies the tangential in-plane displacement of the material points of the membrane.Both A and B depend on the mechanical properties of the membrane as well as the imposed flow field.The mechanical behavior of the membrane material (assumed to be isotropic and purely elastic) is specified by a strain energy function W defined as internal strain energy per unit of initial membrane area.For small deformations, W can be expressed in the following form [58]: where E s is the membrane surface elastic modulus, λ 1 and λ 2 are extension ratios along the principal directions, and α i are material elastic coefficients depending on the type of membrane.In the literature on capsules, the following two types of membranes are often mentioned: MR (Mooney-Rivlin) type membrane and RBC (red blood cells) type membrane.The MR type membrane is made up of an infinitely thin sheet of isotropic, volume-incompressible elastomer, with the following values of α i : The RBC type membrane is made up of a lipid bilayer lined with a protein network that imparts elasticity to the membrane.The RBC type membrane is easily shearable but offers considerable resistance to any increase in local surface area, and has the following values of α i : The time evolution of the deformation tensors A and B is given as: (Ca) where Ca is the capsule capillary number defined as η c .
γR/E s , D/Dt is the Jaumann derivative, L and M are symmetric and traceless tensors which are linear functions of deformation tensors A and B, and a 0 , b i are viscosity-ratio dependent parameters.The tensors L and M as well as the parameters a 0 , b i are given below: The dipole strength S 0 of a single capsule located in an infinite matrix fluid has been determined by Barthes-Biesel and Rallison [58] by considering first-order (small) deformation of capsule.
The expression for S 0 is as follows [9,58]: Therefore, the rheological constitutive equation for a dilute suspension of identical capsules in a Newtonian matrix fluid can be expressed as (see Equation ( 8)): In steady simple shear flow, For MR type membranes, the expressions for A 12 and B 12 are as follows: where For RBC type membranes, the expressions for A 12 and B 12 are as follows: From Equation ( 116), it can be readily shown that the effective viscosity of a dilute suspension of capsules with MR type membranes is as follows: where Based on Equation ( 129), the limiting intrinsic viscosities of suspension of capsules (with MR type membranes) at low and high capillary numbers are as follows [9]: For dilute suspensions of capsules with RBC type membranes, Equation ( 116) leads to the following expression for the effective viscosity: Based on Equation ( 133), the limiting intrinsic viscosities of suspension of capsules (with RBC type membranes) at low and high capillary numbers are as follows [9]: The rheology of dilute suspensions of solid core-porous shell particles has been investigated by Zackrisson and Bergenholtz [32].The solid core-porous shell particles simulate sterically-stabilized soft particles consisting of solid core and a hairy polymer and/or surfactant coating.The dipole strength of a single electrically-neutral solid core-porous shell particle immersed in an infinite matrix of Newtonian incompressible fluid can be expressed as follows [9,32,77]: where a is the radius of the solid core, b is the radius of the composite core-shell particle, β is defined in terms of the permeability (α) of the porous shell as: β = b/ √ α, and L 1 (a/b, β) and L 2 (a/b, β) are given as: It should be pointed out that in the physical chemistry literature where sterically-stabilized particles (solid core-porous shell particles) are widely discussed [77], the friction coefficient f is used, instead of permeability α.The friction coefficient f is related to permeability α as f = η c /α and therefore β = b f /η c .
Using Equation ( 136), the constitutive equation for dilute suspensions of solid core-porous shell particles can be expressed as: From Equation (139), it follows that: [η] = 5 2 For impermeable shells, α → 0 and β = ∞ and therefore, the core-shell particles behave as hard spheres with L 2 /L 1 → 1.0 .When α → ∞ and β = 0, the porous shell vanishes and L 2 /L 1 → a 3 /b 3 ; in this case, the suspension consists of hard particles of radius a. Ohshima [77] has studied the influence of primary shear electroviscous effect on the viscosity of dilute suspensions of soft particles.Soft particles are solid core-porous shell particles where a solid core particle is covered with an ion-penetrable surface layer of polyelectrolytes (porous shell).The polyelectrolyte layer consisted of uniformly distributed charged groups of valence Z.The charge density of the immobile charges (fixed on the polyelectrolyte chains) is ZeN where N is the density of charged groups in the polyelectrolyte layer.The core particle surface is assumed to be uncharged.The effective viscosity of dilute suspensions of such soft particles is given as: Under the condition of low fixed charge density (low ZeN), the primary electroviscous coefficient p can be expressed as [77]: where z i is the valence of the ionic species i, λ i is the drag coefficient of ionic species i, λa = a f /η c = (a/b)β and ψ 0 is the surface potential at the outer surface of the porous layer (r = b) given as: The expression for L(κa, a/b, λa) under the restrictions κa >> Figure 24 shows the plots of L versus κa for two different values of a/b.The value of λa is fixed to be 75.Interestingly L (and hence p) exhibits a minimum at a certain value of κa and the minimum shifts to higher κa with the increase in a/b.No such minimum is observed in suspensions of electrically charged hard spheres. , the viscosity of the outer shell liquid be 2  , and the viscosity of the inner core liquid be 3  .

Solid Core-Liquid Shell
The rheology of dilute suspensions of solid core-liquid shell particles in another immiscible liquid was investigated by Davis and Brenner [78].Let the solid core be of radius a and the outer liquid shell be of thickness of b − a.Let the viscosities of the shell and matrix liquids be η 2 and η 1 , respectively.Assuming negligible deformation of the liquid shell, Davis and Brenner [78] obtained the following expression for the dipole strength S o of a single solid core-liquid shell particle immersed in an infinite matrix of Newtonian incompressible fluid: where ψ(a/b, λ) is a dimensionless function of viscosity ratio λ(= η 2 /η 1 ) and radii ratio (a/b), and E ∞ is the rate of deformation tensor far away from the particle.The function ψ is given as follows: Using the expression for the dipole strength (Equation ( 146)), the rheological constitutive equation for a dilute suspension of identical spherical particles of solid-core-liquid shell type can be expressed as: From Equation ( 148), the relative viscosity of dilute suspension is: When a = b, Equation (147) gives ψ = 1 and Equation ( 149) reduces to the Einstein viscosity equation for dilute suspensions.When a = 0, the solid core vanishes, ψ = λ+ (2/5)  1+λ and Equation (149) reduces to the Taylor viscosity equation for dilute emulsions.

Liquid Core-Liquid Shell
Liquid core-liquid shell particles are also referred to as double emulsion droplets.Let the inner core drop be of radius a and the outer liquid shell be of thickness b − a.Let the viscosity of the continuous-phase fluid be η 1 , the viscosity of the outer shell liquid be η 2 , and the viscosity of the inner core liquid be η 3 .
Assuming negligible deformation of double emulsion droplet, Stone and Leal [79] determined the velocity and pressure distribution interior and exterior to the globule.They also determined the following expression for the dipole strength S o of a single double-emulsion droplet located in an infinite matrix fluid: where I(λ 21 , λ 32 , a/b) is a dimensionless function of viscosity ratios λ 21 (= η 2 /η 1 ) and λ 32 (= η 3 /η 2 ), and radii ratio (a/b), given as: The coefficients δ i depend on a/b and are given as:  154) Using the expression for S o (Equation ( 150)), it can be readily shown that the rheological constitutive equation for a dilute suspension of identical spherical double-emulsion droplets is given by: From Equation (156), it follows that the relative viscosity of a dilute suspension of double-emulsion droplets is: When the outer shell liquid of double-emulsion droplets is very viscous ( λ 21 → ∞ ), I → 0 and consequently Equation ( 157) reduces to the Einstein viscosity equation for dilute suspensions of rigid spheres.When a = b, I → 1/(1 + λ 21 ) and therefore, Equation ( 157) reduces to the Taylor viscosity equation for dilute emulsions.It should also be noted that double emulsion droplets with very thin outer shells behave like rigid particles regardless of the viscosities of the inner core liquid and shell liquid.

Rheology of Suspensions of Rigid Non-Spherical Particles
The rheological behavior of suspensions of rigid non-spherical particles has been studied extensively [8][9][10].The rheology of dilute suspensions of non-spherical particles (assumed axially symmetric) is strongly dependent on the shape, aspect ratio (length of particle along its axis of symmetry to its length perpendicular to the axis of symmetry), and orientation distribution function of particles.Let → p be a unit vector parallel to the axis of symmetry of particle (and locked into the particle) that specifies the orientation of the particle.Let Ψ( → p , t) be the orientation distribution function of particles that specifies the orientation state of the particles present in the suspension at any time t.The orientation distribution function Ψ( → p , t) could be interpreted in two ways: (a) the probability that a particle is oriented in a neighborhood of Thus the calculation of the average dipole strength, and hence the bulk stress, of a dilute suspension of non-spherical particles requires the knowledge of the orientation distribution function Ψ( → p , t).Assuming that particles undergo rotary Brownian motion, the orientation distribution function of suspension can be obtained from the solution of the following Fokker-Planck equation in → p space: where div refers to divergence and D r is the rotary diffusivity of particles.
The average dipole strength of a dilute suspension of identical non-spherical particles, defined in Equation (158), is related to the rate of deformation tensor E of suspension as follows: where V p is the volume of a particle, A, B, C and F are shape factors dependent on the particle aspect ratio.The second and fourth moments of → p are defined as: Using Equations ( 7) and ( 160), the rheological constitutive equation of a dilute suspension of identical non-spherical particles can be written as [9]: For long thin rod-shaped particles ( r → ∞ ), the shape factors A, B, C and F are as follows: For flat circular disk-shaped particles ( r → 0 ), the shape factors A, B, C and F are as follows: Let us now consider two different cases: Case A where the axisymmetric non-spherical particles of the suspension are large enough to ignore the rotary Brownian motion, that is, the particles are non-Brownian (D r = 0); and Case B where the axisymmetric non-spherical particles are small enough that the rotary Brownian motion is significant, that is, the non-spherical particles are Brownian (D r > 0).
Case A: When non-Brownian axisymmetric particles are subjected to shear flow, they undergo rotational motion in fixed orbits, referred to as Jeffery orbits.Each particle rotates in a fixed orbit.The distribution of particles in different orbits is determined by the initial orientation state of the particles at t = 0. Assuming that all the particles of the suspension are identical, each particle possesses the same period of rotation T given as [9]: where r is the aspect ratio and .
γ is the shear rate.In this case, the orientation distribution function of particles Ψ is periodic in nature.As the bulk stress depends on the orientation state of particles, the rheological properties are also periodic functions of time.For example, the instantaneous intrinsic viscosity of suspension is small during that portion of the time period T when the particles are aligned with the flow field, and is large during that portion of the time period T when the particles are not aligned with the flow field, that is, the intrinsic viscosity oscillates with time period T. With the increase in shear rate .γ, the time period T decreases and the rheological properties oscillate at a faster rate.Interestingly, the rheological properties of a suspension of identical, non-Brownian and non-spherical particles, are indeterminate if the initial orientation state of the particles, specified by Ψ at t = 0, is not known.Note that Ψ = 1/4π at t = 0 if the particles are initially randomly oriented.
Case B: Small particles (particles with their largest dimension less than a few microns) undergo both translational and rotational Brownian diffusion, that is, random translational displacement and random angular displacement.However, it is only the rotary Brownian motion that affects the rheology of suspensions of non-spherical particles.The rotary Brownian motion affects the flow behavior of suspensions of non-spherical particles in two ways: (i) it alters the angular velocity of the particles.This is referred to as direct effect of Brownian motion on the rheology of suspension.The bulk stress, and hence rheological properties, of suspension are sensitive to the angular velocity of particles.Consequently a change in the angular velocity of particles due to rotary Brownian motion produces a corresponding change in the rheological properties of suspension.An important point to note is that the direct effect of rotary Brownian motion on suspension rheology is present even in the absence of imposed macroscopic motion, that is, even when E = 0, and is often referred to as diffusion stress; and (ii) the rotary Brownian motion also affects the orientation distribution function Ψ.It tends to randomize the orientations of the particles.This is referred to as indirect effect of Brownian motion on suspension rheology.Interestingly, the bulk stress (and hence rheological properties) of suspension in the presence of rotary Brownian motion is determinate under all conditions, regardless of the initial orientation state of the particles.If sufficient time is allowed, the orientation distribution Ψ reaches a unique steady state value that is independent of the initial orientation state of particles.It should be noted that the particles still rotate but not in fixed orbits.There is no equilibrium orientation of particles.Only Ψ is stationary in time indicating that for any orientation → p as many particles leave the neighborhood of → p as enter it.Interestingly, the dilute suspensions of non-spherical Brownian particles exhibit shear-thinning and normal-stresses when subjected to simple shear flow.The rotary Brownian motion tends to randomize the orientations of the particles whereas the viscous stresses tend to orient the particles such that they are least disruptive to flow.Thus the steady state value of the orientation distribution function Ψ varies with the rotary Peclet number, defined as Pe = .γ/D r .At low values of Pe, the viscosity of suspension is high as Brownian motion dominates and the particles are randomly oriented.At high values of Pe, the viscous stresses dominate over the Brownian motion and the particles become aligned with the flow resulting in lower viscosity of suspension.
Chen and Koch [131] have explored the influence of particle electric charge on the rheology of dilute suspensions of Brownian rod-shaped particles (charged fibers) with large aspect ratio under the conditions of small Hartmann number and thick double layer (κ −1 >> particle radius).The presence of electric charge on the particles leads to three additional contributions to the bulk stress of suspension: two direct effects and one indirect effect.One direct effect comes from the distortion of electrical double layer when the suspension is subjected to macroscopic motion.The motion of fluid surrounding the charged particle distorts the electric charge cloud from its equilibrium state, resulting in an electrical stress.Interestingly, the distortion of electrical charge cloud also results in an electric torque on the particle.The electric torque on the particle affects the angular velocity of the particles.This is the second direct effect of electric charge on the bulk stress of suspension.The electric torque on the particles also affects the orientation distribution function Ψ of particles.This is the indirect effect of electric charge on the rheology of suspensions.

Rheology of Suspensions of Ferromagnetic Particles
Suspensions of ferromagnetic particles are referred to as ferrofluids in the literature [132].A number of studies have been published on the mechanics, rheology, and applications of ferrofluids .The particles of ferrofluids are usually small in size (about 10 nm in diameter) and are made up of some ferromagnetic material such as iron oxide.The atoms of ferromagnetic materials possess a net magnetic dipole moment due to the presence of some unpaired electrons and hence incomplete cancellation of spin and orbital magnetic moments.Note that each electron in an atom possesses a magnetic moment due to two reasons: (1) orbital rotation of electron; and (2) spinning of electron.If the electrons are present in pairs, they cancel out each other's magnetic moments (both spin and orbital) as paired electrons spin and orbit in opposite directions with respect to each other.A large sized ferromagnetic material generally does not possess a net magnetic moment due to random orientation of many domains present in the material.However, due to the small size of the ferromagnetic nanoparticles, the particles behave as single-domain permanent magnets.
In shear flow of ferrofluids, the interaction of magnetic nanoparticles with the external magnetic field can result in some unusual rheological properties.In the absence of external magnetic field, the suspended nanoparticle rotates freely with the local vorticity of the fluid motion.However, when an external magnetic field is applied, the particle experiences a magnetic torque whenever its dipole vector is not aligned with the external field.Due to the presence of a magnetic torque, the particle is no longer able to rotate freely with the fluid.For example, the particle does not rotate at all when the magnetic field strength is very large and the magnetic field is perpendicular to the vorticity of the flow.In this case, the dipole vector of the particle aligns itself with the magnetic field and does not rotate with the fluid vorticity.The hindered rotation of magnetic nanoparticles increases the rate of mechanical energy dissipation in ferrofluid resulting in an enhancement of the fluid viscosity.This field-induced increase in the viscosity of ferromagnetic fluids is referred to as rotational viscosity and was first observed experimentally by McTague [133].The rotational viscosity is zero only in the special case when the external magnetic field is aligned with the vorticity of the flow.In this case the dipole vector of the particle aligns itself with the external field and the particle rotates freely with the local vorticity of the fluid motion.Interestingly, the bulk stress of suspension is no longer symmetric when the suspension is subjected to an external magnetic field.
Consider now a spherical magnetic nanoparticle of radius R suspended in an unbounded Newtonian liquid (viscosity η c ) and subjected to shear flow → u = (0, .γ x, 0).Let m be the strength of the magnetic dipole moment and → p be the unit vector locked into the particle in the direction of the dipole moment.When an external magnetic field → H is applied, the particle experiences a magnetic torque.The magnetic torque tends to align the particle such that its dipole vector → p is co-linear with the magnetic field → H whereas the hydrodynamic torque exerted by the fluid tends to align the particle such that its angular velocity vector → ω is co-linear with vorticity vector associated with the undisturbed flow.If the rotary Brownian motion of particles is absent, all the particles of the suspension attain a fixed orientation at steady state, except when γ = π/2 and 0 < λ < 1, where γ is the polar angle between the magnetic field → H and undisturbed vorticity vector of fluid 2 → ω 0 (see Figure 25), and λ is a dimensionless measure of the magnetic field strength with respect to the strength of hydrodynamic torque, defined as: The steady orientation of the particles, expressed in terms of the angles θ s and φ s (see Figure 25), is given as: The average dipole strength of a dilute suspension of magnetic spherical particles is related to the rate of deformation tensor E of suspension as follows [9]: Using Equations ( 7) and (177), the rheological constitutive equation for a dilute suspension of magnetic spherical particles can be written as: In the present situation where the dipole vectors of all particles are directed in the same direction the orientation distribution function  becomes the Dirac delta function: where  is the Dirac delta function and The average dipole strength of a dilute suspension of magnetic spherical particles is related to the rate of deformation tensor E of suspension as follows [9]: where Using Equations ( 7) and (177), the rheological constitutive equation for a dilute suspension of magnetic spherical particles can be written as: In the present situation where the dipole vectors of all particles are directed in the same direction the orientation distribution function Ψ becomes the Dirac delta function: It should be noted that in the special case where γ = π 2 and 0 < λ < 1, the particle dipole vector → p does not align in any fixed direction, but instead rotates forever in one of an infinite family of closed orbits with a period of rotation as 2π/ω 0 √ 1 − λ 2 where ω 0 is half the vorticity of the undisturbed shear flow.The axis of rotation of → p lies in a plane perpendicular to the direction of the applied magnetic field.The dipole vector → p sweeps over the surface of a circular cone as it rotates in a particular orbit.The distribution of particles in different orbits depends on the initial orientation state of the suspension.The orientation distribution function Ψ, and hence the bulk stress of suspension, is periodic in nature.Furthermore, the bulk stress is indeterminate if the initial orientation state of the suspension is not known.
In simple shear flow → u = (0, .γx, 0), Equation (181) predicts that the intrinsic viscosity of suspension varies from 2.5 to 4 depending on the steady state polar angle (θ s ) of dipole vector → p (see Figure 25) and is given as [135]: The polar angle θ s depends on magnetic field angle γ and dimensionless magnetic strength λ (see Equation ( 174)).For low values of λ(<<1), the hydrodynamic torque acting on the particle dominates over the magnetic torque and the particles are able to rotate freely with the fluid.Consequently the intrinsic viscosity [η] is 2.5 corresponding to the Einstein value for suspension of rigid spheres in the absence of any magnetic field.When λ >> 1, the particles are not able to rotate freely with the fluid and therefore, the intrinsic viscosity is higher than 2.5 provided that the angle γ between the magnetic strength and undisturbed vorticity vector of fluid is not zero.When γ = π 2 and λ >> 1, the particles do not rotate at all and the intrinsic viscosity is the highest, that is, [η] = 4.The constitutive equation, Equation (181), also predicts that there are no normal stresses in a suspension and that the stress tensor is not symmetric.
Let us now consider the effect of rotary Brownian motion on the rheology of suspension of magnetic particles.The rotary Brownian motion has a strong influence on the rheology of suspension of magnetic particles; it tends to randomize the orientation of the magnetic dipole moments ( → p vectors) of the particles and hence affects the orientation distribution function Ψ [136].The orientation distribution function Ψ is no longer a Dirac delta function; it is given by the Fokker-Planck equation in → p space (Equation (159)).If sufficient time is allowed, the orientation distribution Ψ reaches a unique steady state value that is independent of the initial orientation state of dipole vectors.Interestingly the bulk stress of a suspension, given by Equation (179), is now determinate under all conditions, regardless of the initial orientation state of particles.The presence of rotary Brownian motion eliminates the indeterminancy mentioned earlier in the special case where γ = π/2 and 0 < λ < 1.The rheological properties of suspension of Brownian magnetic particles when subjected to external magnetic field are governed by the following quantities: γ, λ and Pe, where γ is the angle between the magnetic field and undisturbed vorticity vector (ω 0 ), λ is a measure of magnetic strength relative to the strength of hydrodynamic couple (see Equation (173)), and Pe is the rotary Peclet number defined as Pe = .
γ/D r (some authors define the rotary Peclet number as Pe = ω 0 /D r = .γ/2D r ).For specified γ and λ the viscosity of the suspension decreases with the increase in Pe indicating shear-thinning behavior of suspension.There are no normal stresses present in the suspension and the stress tensor is antisymmetric.
Thus far in this section we have covered the rheology of suspensions of spherical dipolar (magnetic) particles, in the absence and presence of rotary Brownian motion.In what follows, the rheological behavior of suspensions of magnetic non-spherical particles is considered.Suspensions of magnetic non-spherical particles (such as rod-like and plate-like particles) are equally important from a practical point of view.Several authors [139,141,142,158,168] have investigated the mechanics and rheology of suspensions of non-spherical magnetic particles when subjected to external magnetic field in the absence and presence of rotary Brownian motion.In the absence of rotary Brownian motion, the motion of non-spherical magnetic particles is governed by hydrodynamic and magnetic torques acting on the particles.The hydrodynamic torque tends to orient the particle in the direction of flow whereas the magnetic torque tends to align the particle dipole with the magnetic field.Under certain conditions of external magnetic field orientation and field strength, the particles move in a family of closed orbits.Under other conditions of external field orientation and field strength, the particles reach a stable equilibrium orientation.For example, when the external magnetic field is perpendicular to fluid vorticity axis (polar angle of external magnetic field γ is π/2, see Figure 26), the motion of a non-spherical magnetic particle depends on the azimuthal angle β of the field orientation and the field strength parameter λ defined as: where µ is the mobility coefficient corresponding to rotation of particle about a transverse axis and m is the dipole strength of the particle.Note that this λ has the same significance as that of λ defined in Equation ( 173) for spherical dipolar particles.When β = 0 or π/2 and λ is small, the particle moves in a family of closed orbits.When 0 < β < π/2 and λ is small, the particle either reaches a stable equilibrium orientation or approaches a limit cycle on the equator [141].For large λ, the particle converges to a stationary stable orientation, regardless of the value of β.When the external magnetic field does not act in the plane of shear, that is γ = π/2, and cosγ ≈ order(1) particles approach the stable equilibrium orientation regardless of the value of the field strength parameter λ.It is expected that the bulk stress in a suspension of non-Brownian non-spherical magnetic particles is periodic when the particles move in fixed orbits and the orientation distribution function is periodic.Under such conditions, the bulk stress is determinate only if the initial orientation distribution of particles is known.
In those situations where the particles reach a stable orientation, the bulk stress is expected to reach a steady value regardless of the initial orientation of the particles.
In the presence of rotary Brownian motion, the bulk stress is determinate under all conditions regardless of the initial orientation state of particles.The particles do not rotate in fixed orbits.The orbits are stochastic in nature.If sufficient time is allowed, the orientation distribution function Ψ, given by the Fokker-Planck equation (Equation ( 159)), reaches a unique steady state valve independent of the initial orientation state of the particles.The suspensions of non-spherical Brownian magnetic particles, when subjected to magnetic field, exhibit non-Newtonian rheological characteristics such as: non-zero normal stresses and shear-rate dependent viscosity.Also the stress tensor is non-symmetric.Both Brownian and magnetic couples make a direct contribution to bulk stress by affecting the angular velocity of the particles.In addition to the direct contribution to bulk stress, the Brownian and magnetic forces also make an indirect contribution to the bulk stress by affecting the orientation distribution of particles.The steady-state rheological properties of suspensions of magnetic non-spherical Brownian particles depend on five parameters: polar angle γ of the magnetic field (see Figure 26), azimuthal angle β of the magnetic field, dimensionless magnetic field strength H defined as m H/kT, aspect ratio of particles, and rotary Peclet number defined as  Unlike suspensions of spherical magnetic particles, the rheological properties of non-spherical magnetic particles depend on both polar and azimuthal angles of the external magnetic field.Furthermore, the suspensions of non-spherical magnetic particles exhibit resonance phenomenon [139] at intermediate values of rotary Peclet number in strong magnetic fields with certain orientations.The intrinsic viscosity versus rotary Peclet number plot exhibits a maximum at resonance.For example, consider the flow of suspensions of magnetic oblate spheroids with aspect ratio of 0.4 [139].The suspension is subjected to simple shear flow given as

Concluding Remarks
The fundamental rheology of dilute disperse systems based on single-particle mechanics is reviewed.A variety of different types of inclusions are considered.Our current understanding of the theoretical rheology of dilute disperse systems is generally very good.However, some areas require further theoretical work.For example, the influence of primary electroviscous effect (shear and dilational) on the rheology of dilute disperse systems at arbitrary flow strengths and surface potentials is generally lacking.It should also be noted that most of the theoretical work on the rheology of dilute disperse systems is restricted to disperse systems consisting of Newtonian matrix fluid.Little progress has been made on the development of rigorous theoretical models for the rheology of dilute disperse systems consisting of non-Newtonian matrix fluid, especially at arbitrary flow strengths.In many cases, experimental data are lacking to validate the existing theoretical models and to determine the range of their validity.More experimental work is needed on the rheology of well-characterized dilute disperse systems in the following cases: (a) suspensions of core- Unlike suspensions of spherical magnetic particles, the rheological properties of non-spherical magnetic particles depend on both polar and azimuthal angles of the external magnetic field.Furthermore, the suspensions of non-spherical magnetic particles exhibit resonance phenomenon [139] at intermediate values of rotary Peclet number in strong magnetic fields with certain orientations.The intrinsic viscosity versus rotary Peclet number plot exhibits a maximum at resonance.For example, consider the flow of suspensions of magnetic oblate spheroids with aspect ratio of 0.4 [139].The suspension is subjected to simple shear flow given as → u = ( .γ y, 0, 0).The magnetic field orientation is specified as (0, 1, 0), that is, γ = π/2 and β = π/2 (the field is in xy plane, perpendicular to flow).The dimensionless filed strength H is 10.Under these conditions, the intrinsic viscosity versus rotary Peclet number plot shows a sharp peak at Pe ≈ 11.At low and high values of Pe, [η] ≈ 2.62 whereas at an intermediate Pe ≈ 11, [η] exhibits a peak value of 3.7.At low values of Pe, the particles are held in streamline configuration by the strong magnetic field.At high values of Pe, the shear flow is strong enough to keep the particles aligned with the flow.At intermediate values of Pe where resonance is observed, the particles spend significant time oriented in the direction of maximum strain resulting in high value of intrinsic viscosity.

Concluding Remarks
The fundamental rheology of dilute disperse systems based on single-particle mechanics is reviewed.A variety of different types of inclusions are considered.Our current understanding of the theoretical rheology of dilute disperse systems is generally very good.However, some areas require further theoretical work.For example, the influence of primary electroviscous effect (shear and dilational) on the rheology of dilute disperse systems at arbitrary flow strengths and surface potentials is generally lacking.It should also be noted that most of the theoretical work on the rheology of dilute disperse systems is restricted to disperse systems consisting of Newtonian matrix fluid.Little progress has been made on the development of rigorous theoretical models for the rheology of dilute disperse systems consisting of non-Newtonian matrix fluid, especially at arbitrary flow strengths.In many cases, experimental data are lacking to validate the existing theoretical models and to determine the range of their validity.More experimental work is needed on the rheology of well-characterized dilute disperse systems in the following cases: (a) suspensions of core-shell particles (solid core-porous shell, solid core-liquid shell, double emulsion droplets); (b) suspensions of porous particles and capsules; (c) suspensions of non-spherical particles; (d) suspensions of electrically-charged inclusions

Figure 1 .
Figure 1.Comparison of Einstein viscosity equation with experimental data for viscosity of suspensions of hard spheres.

)Figure 1 .
Figure 1.Comparison of Einstein viscosity equation with experimental data for viscosity of suspensions of hard spheres.

Figure 2 .
Figure 2. Charged particle surrounded by an ionic cloud [9].Variations of ionic concentration and electric potential are also shown.(a) Charged particle surrounded by an ionic cloud; (b) Variation of ionic concentration with distance from the particle surface; (c) Variation of electric potential with distance from the particle surface.

Figure 2 .
Figure 2. Charged particle surrounded by an ionic cloud [9].Variations of ionic concentration and electric potential are also shown.(a) Charged particle surrounded by an ionic cloud; (b) Variation of ionic concentration with distance from the particle surface; (c) Variation of electric potential with distance from the particle surface.

1 versusFigure 3 .
Figure 3. Variation of primary electroviscous coefficient p with dimensionless surface potential ψ 0 for different values of Rκ.
. The plots are shown for two different electrolyte (KCl) concentrations.As expected the experimental data exhibit a linear behavior with slopes of −[η].The intrinsic viscosity values for the data shown in Figure 4 are as follows: [η] = 6.9 at KCl concentration of 4 × 10 −5 M and [η] = 5.5 at KCl concentration of 1.2 × 10 −4 M. The corresponding p values (primary electroviscous coefficient) are: 1.76 and 1.20, respectively.

Figure 4 . 1 versus
Figure 4. Plots of r  / 1 versus  data for polystyrene latex suspensions consisting of negatively charged particles of diameter 250 nm suspended in KCl solutions.

Figure 4 .
Figure 4. Plots of 1/η r versus φ data for polystyrene latex suspensions consisting of negatively charged particles of diameter 250 nm suspended in KCl solutions.

Figure 5 .
Figure 5. Variation of primary electroviscous coefficient p with electrolyte concentration for electrically charged polystyrene latex suspensions.

2 Figure 5 .
Figure 5. Variation of primary electroviscous coefficient p with electrolyte concentration for electrically charged polystyrene latex suspensions.

Figure 6 .
Figure 6.Variation of [η] with Pe for suspensions of charged particles at different values of Hartmann number (Rk fixed at 8).

Figure 6 .
Figure 6.Variation of [η] with Pe for suspensions of charged particles at different values of Hartmann number (Rk fixed at 8).

Figure 7 .
Figure 7. Variation of [η] with Pe for suspensions of charged particles at different values of Rk (Ha fixed at 0.5).

Figure 8 numberversusFigure 7 .
Figure8shows the plots of first normal-stress difference as

Fluids 2016, 1 , 40 15 of 52 Figure 10
Figure10compares the bulk and shear primary electroviscous coefficients at a small Hartmann number of 0.1.The predictions of Equations (34) and (47) are compared for thin double layers ( 

Fluids 2016, 1 , 40 15 of 52 Figure 10
Figure10compares the bulk and shear primary electroviscous coefficients at a small Hartmann number of 0.1.The predictions of Equations (34) and (47) are compared for thin double layers ( 

Fluids 2016, 1 , 40 16 of 52 Figure 10 .
Figure 10.Comparison of bulk and shear primary electroviscous coefficients at a small Hartmann number of 0.1.

Figure 10 .
Figure 10.Comparison of bulk and shear primary electroviscous coefficients at a small Hartmann number of 0.1.

Figure 12 .
Figure 12.Influence of particle permeability on the intrinsic viscosity of suspension of electricallycharged porous particles.
Figure 12.Influence of particle permeability on the intrinsic viscosity of suspension of electrically-charged porous particles.

Figure 13 .
Figure 13.Intrinsic viscosity ] [ as a function of Se for suspensions of soft solid particles.

Figure 15 .
Figure 15.Reduced second normal-stress difference −N 2r /φ as a function of Se for suspensions of soft-solid particles.

Fluids 2016, 1 , 40 23 of 52 Figure 16 .
Figure 16.Relative viscosity r  versus droplet concentration  for dilute emulsions for different values of the viscosity ratio  .

− 4 .
The emulsions are produced by mixing the two liquids in an agitated vessel at different agitator speeds.The intrinsic viscosity (Equation (68)) of these emulsions is 1.0013.The experimental data show good agreement with the Taylor viscosity equation when 05 .0   .At higher concentrations of dispersed-phase, the Taylor equation under predicts the emulsion viscosity.
by considering the first-order deformation of droplets and developed the following equation for the dipole strength of a droplet immersed in an infinite matrix of incompressible Newtonian fluid:

69 )Figure 16 .
Figure 16.Relative viscosity η r versus droplet concentration φ for dilute emulsions for different values of the viscosity ratio λ.

Figure 17 .
Figure 17.Comparison between experimental viscosity data and theoretical predictions for emulsions.

Figure 17 .
Figure 17.Comparison between experimental viscosity data and theoretical predictions for emulsions.

Figure 18 .
Figure 18.Intrinsic viscosity of emulsion as a function of capillary number for different values of viscosity ratio  .The Frankel and Acrivos equation, Equation(71), also predicts non-zero normal stress differences in shear flow of dilute emulsions.Equation (71) leads to the following expressions for the normal stress differences:

1 versus 2 
The plots of normal stress differences for emulsions are shown in Figures19 and 20.In Figure19, the first normal stress difference is plotted as  r N Ca , keeping 1   .In Figure20, the second normal stress difference is plotted as  r N

Figure 18 .
Figure 18.Intrinsic viscosity of emulsion as a function of capillary number for different values of viscosity ratio λ.

Figure 20 .
Figure 20.Reduced second normal stress difference

Figure 20 .
Figure 20.Reduced second normal stress difference

Fluids 2016, 1 , 40 30 of 52 Figure 22 .
Figure 22.Plots of r  versus Ma m _ for emulsions for different values of the viscosity ratio  .

Figure 23 
Figure 23 show the plots of r  versus  for different values of the interfacial viscosity term

Figure 23 .Figure 22 .
Figure 23.Plots of r  versus  for different values of the interfacial viscosity term

Figure 22 .
Figure 22.Plots of r  versus Ma m _ for emulsions for different values of the viscosity ratio  .

Figure 23 
Figure 23 show the plots of r  versus  for different values of the interfacial viscosity term

Figure 23 .Figure 23 .
Figure 23.Plots of r  versus  for different values of the interfacial viscosity term

Fluids 2016, 1 , 40 38 of 52 Figure 24 .
Figure 24.Plots of L versus a  for two different values of b a / .

25 Figure 24 .
Figure 24.Plots of L versus κa for two different values of a/b.

→pp
at time t is given by Ψ( → p , t)d → p ; and (b) the fraction of the total number of particles in the suspension with orientations lying between → ) be the dipole strength of a reference particle with orientation → p , located in an infinite matrix fluid.Therefore, the average dipole strength of a dilute suspension can be expressed as:

)Fluids 2016, 1 , 40 43 of 52 Figure 25 .
Figure 25.Coordinate system.The magnetic field is assumed to be in X-Z plane with azimuthal angle of zero degree.The particles spin around the p  axis with an angular velocity of s p  is the dipole vector of the particles at steady state.Consequently, p  is the same as s p  and the rheological constitutive equation (Equation (179)) becomes:

Figure 25 .
Figure 25.Coordinate system.The magnetic field is assumed to be in X-Z plane with azimuthal angle of zero degree.

Figure 26 .
Figure 26.Coordinate system.The magnetic field has a polar angle of  and an azimuthal angle of  .
is in xy plane, perpendicular to flow).The dimensionless filed strength H ~ is 10.Under these conditions, the intrinsic viscosity versus rotary Peclet number plot shows a sharp peak at a peak value of 3.7.At low values of Pe , the particles are held in streamline configuration by the strong magnetic field.At high values of Pe , the shear flow is strong enough to keep the particles aligned with the flow.At intermediate values of Pe where resonance is observed, the particles spend significant time oriented in the direction of maximum strain resulting in high value of intrinsic viscosity.

Figure 26 .
Figure 26.Coordinate system.The magnetic field has a polar angle of γ and an azimuthal angle of β.
Pe .With the increase in the electrical double layer thickness, that is, decrease in Rk , the intrinsic viscosity in the upper Newtonian region ( 1  Pe ) increases at a fixed value of Hartmann number Ha (see Figure