Self-diffusion in nanofluids of nonelongated particles in the dilute limit

The dynamic features of a dilute suspension of nanoparticles (nanofluid) are fully modified depending on the dominant particles slip mechanism acting in the suspension. Self-diffusion effects in highly sheared diluted suspensions (entrance conditions and microapplications) can lead to a particles distribution fully different from the bulk one. The reported investigation proposes a model to determine the self-diffusion of three-planes symmetric nonelongated particles inmersed in a sheared Stokes flow. The model is based on the real displacements between any pair of particles and an statistical approach to determine contact kinematic irreversibilities. According to the proposed model, the source of hydrodynamic irreversibility is closely related to the particles shape. This is clearly demonstrated through the application of the model to cubic particles. The main conclusion is that the particles shape plays a significant role in the dynamic behavior of the suspension and, as a result, in the self-diffusion coefficient. The reported results arising from the cubic particles trajectories in a Stokes flow are reasonably close to the ones reported by Brady and Morris (J Fluid Mech, 348:103–139, 1997) for suspensions under high Pe number, and Zarraga and Leighton (Phys Fluids 13(3):565-577, 2001).


Introduction
Nanofluids production and thermal performance have been intensely investigated and reported during the last twenty years [8,22]. Most of the investigations have stressed the convective heat transfer enhancement associated with nanofluids though, as expected, accompanied by a corresponding pressure drop increment. Their thermal performance has been related to a nanoparticles concentration field affected by different diffusion mechanisms [5]. It is interesting to note that though a number of investigations have obtained significant convective heat transfer increments, as suggested by [24], others have obtained slight enhancements and even decrements in the heat transfer rate under internal flow conditions of nanofluids [11]. In any case, it can be concluded that the nanoparticles concentration field is affected when the nanofluid is submitted to heat transfer from a heated wall, and, as a conclusion, thermal performance is closely related to this field [10]. The concentration field is significantly affected in the region close to the heated (or cooled) wall due to thermophoretic diffusion effects. Thus, it is certain that the mechanisms of nanoparticles transport with respect to the base fluid play an important role in determining the behavior of the nanofluid under heat transfer conditions.
Hwang et al. [15] clearly identified a convective enhancement higher than 8.0% for a water alumina nanofluid in the fully developed region of a 1.812-mm inner-diameter pipe. The small amount of particles in suspension, up to 0.3% in volume, does not seem to explain the reported enhancement in terms of the nanofluid modified transport properties due to the particles concentration field. They also stated that the detected convective enhancement depends on the distance Technical Editor: Monica Carvalho. from the inlet of the particular cross section along the fully developed region. Several mechanisms were proposed in [15] in order to explain the detected enhancement: (i) thermophoresis, (ii) self-diffusion and (iii) viscosity gradient.
In the same way, though under a numerical perspective, [19] developed a model inspired in the Buongiorno's one [5]. They reported significant enhancement under fully developed flow conditions for highly loaded nanofluids, up to 5.0% in volume, whereas for loads below 1.0% , as in [15], the enhancement became negligible. As a result, the enhancement observed in the [15] results must be explained in terms of different particles interaction mechanism, not considered in [19], which becomes macroscopically apparent.
Metzger et al. [18] developed an investigation to devise the apparent thermal conductivity of sheared suspensions. They reported an apparent increase in the thermal diffusion coefficient (up to three times) depending on the particles concentration and shear rate. Furthermore, they found a correlation between the apparent thermal diffusion and the measured self-diffusion coefficients. It must be stressed that self-diffusion cannot take place in the absence of kinematic irreversibilities [13], which are generated through two mechanisms: particles collisions [7] and hydrodynamic interactions between particles [23]. The investigation reported by [18] was developed for pure sheared suspensions, Couette laminar flow, with spherical PMMA (polymethylmethacrylate) particles of millimetric size. If the spheres surface is free of asperities in the dilute limit, [1] (which is not the case for nanofluids) cannot generate self-diffusion due to particles collisions. [18] experimentally proved the previous statement, even for their highly loaded suspensions tests, since their developed numerical results considering the particles hydrodynamics alone matched the experimental ones. Therefore, in [18], the effects of particles collisions are neglected or, at least, the collisional mechanism is not the main source of kinematic irreversibilities.
Self-diffusion, also designated as "shear diffusion," has been neglected for dilute suspensions of particles of sizes of the order of nanometers, though as will be demonstrated in this paper could play an important role under certain conditions. Its importance with respect to Brownian diffusion is characterized by the so-called Péclet number, Pe, defined for spherical particles by the following expression: where is the shear rate, a represents the radius of the equivalent sphere, stands for the fluid dynamic viscosity, k is the Boltzmann constant and T is the fluid temperature [K], respectively. According to Eq. 1, the Peclet number is the ratio between shear and Brownian diffusion. Thus, nanoparticles diffusion caused by shear effects should be (1) Pe = 6 a 3 kT disregarded as significant mechanism in nanofluids due to the low values of the Péclet number typical of these suspensions, except either in microsized applications [12], or in the case of nonspherical nanoparticles [21]. In the latter case, self-diffusion might be enhanced orders of magnitude with respect to the case of spherical particles, where there is a dominance of hydrodynamic interactions, due to the increment in the hydrodynamic irreversibility related to the interparticles forces (nonspherical particles) [2]. On the other hand, since characteristic dimensions of microapplications could attain values lower than 100 μ m (see, for example, [6] and [12]), shear diffusion would no longer be a negligible diffusion mechanism in those applications. In fact, for this range of channel sizes, the Péclet number could attain values significantly higher than the unity, Pe ≫ 1 , due to the shear rate which increases as the size of the channel diminishes. A literature survey of self-diffusion in suspensions has revealed two different approaches to its investigation: (i) the one that focuses on the particles trajectories and (ii) from a rheologic point of view, which focuses on the particles "induced stresses" related to the particles. The approach related to the particles trajectories is based on the wellknown Stokes flow, which is hydrodynamically reversible. However, shear or self-diffusion in suspensions is the result of particles induced irreversibility in the flow, a mechanism that was initially devised by [9] in their experimental investigation. Their study consisted in measuring the self-diffusion coefficient in neutrally buoyant spherical and diskshaped particles suspensions. They suggested that self-diffusion is the result of the action of two physical mechanisms, namely: (i) hydrodynamic interactions of a couple of approaching particles which lead to "multiple passing interactions overlap" and a net lateral streamlines displacement; (ii) particles interaction forces resulting from either particles collisions or other causes. Any of the two mechanisms would cause the so-called hydrodynamic irreversibility. In addition, [9] claimed that, for volumetric concentrations lower than 20% , the self-diffusion coefficient D s increases linearly with the volumetric concentration ( ) in such a way that D s a 2 ∝ . It is interesting to note that further investigation of the physical mechanism related to self-diffusion in suspensions by different researchers has confirmed the physical model proposed by [9].
Later on, [17], based on experiments performed in a typical concentric cylinders rheometer, reported self-diffusion results in both the shear plane direction and in the plane normal to it. [17] performed a thorough analysis of the [9] investigation, concluding that, contrary to what those authors claimed, possibly their data were affected by wall effects. As a result, according to [9,17] underestimated the shear diffusion coefficient, which instead they claim that should be proportional to the square of the volumetric concentration, its dimensionless expression being the following D s a 2 ∝ 2 . They also claim that their expression correlates very well [9] data for concentrations below 20%.
[2] developed an experimental investigation on selfdiffusion of dilute suspensions of PPMA (polymethylmethacrylate) particles of average diameter of the order of 771 m. Their experimental setup, consisting of a basic Couette viscometer, was similar to the one of the [17] investigation. The particles of their investigation were eccentric with an average aspect ratio of 1.19. Surface dimensionless roughness (referred to the average particles radius) of the particles was of the order of 1 × 10 −3 . The following are some of the main conclusions of this study: (i) the dimensionless self-diffusion coefficient depends linearly on the volumetric concentration; (ii) the observed particles diffusion was attributed to the breaking of the symmetry of two particles interaction; though the cause was not certain, the authors limited themselves to speculate about the possible mechanisms; (iii) the surface roughness of the particles could not account for its effect; (iv) the dimensionless self-diffusion coefficient obtained in the investigation was 20 times the one obtained by [7] for equally sized spheres and the same surface roughness; [2] claim that in order to obtain the same self-diffusion coefficient of their investigation, the average dimensionless roughness should be of the order of 0.07, that is, 70 times greater than the actual one, what led the authors to state that there must exist a kind of self-diffusion promoting mechanism related to the particles shape.
[23] developed a model aimed at determining the hydrodynamic perturbation caused by two particles over the flow streamlines and a third particle trajectory in the dilute limit of suspensions. By determining the streamlines and particles deflections in the velocity gradient and vorticity directions of a shear flow, they were able to conclude that the dimensionless self-diffusion coefficient is proportional to the square of the volumetric concentration, thus confirming the trend previously observed by [17]. However, their proportionality coefficient was somewhat lower than that proposed by [17]. Thus, it could be stated that [23] provided a theoretical support to the fact that self-diffusion is closely related to the hydrodynamic irreversibility caused by the hydrodynamic interaction of two particles over a third one.
An alternative approach to the hydrodynamic irreversibility was taken by [7], who proposed a model based on the collisions of spherical particles. They assumed that the surface of particles is rough, the roughness being characterized by the "asperity," , a dimensionless parameter defined in terms of the particle radius, a. They claimed that the closest dimensionless distance between the centers that two particles could attain is equal to (2+ ). This concept is closely related to the so-called excluded volume which is the volume surrounding a particle that is not penetrable by other particles, which in the case of the da Cunha and Hinch's model is the spherical ring with inner and external dimensionless radii, respectively, equal to 1 and (1+ 2 ). By determining the trajectories of colliding particles in a shear flow, [7] were able to obtain the self-diffusion coefficient in the velocity gradient and vorticity directions as a function of the particles asperity. According to their model results, the dimensionless self-diffusion coefficient in both directions increases with the asperity and is proportional to the volumetric concentration, that is, D s a 2 ∝ , a trend similar to the one proposed by [9].
The second approach to deal with self-diffusion in dilute suspensions is the one related to the "induced stresses." This procedure is based on the determination of the induced stresses from the rheological performance of the suspension. It is interesting to note that hydrodynamic irreversibility is closely related to the non-Newtonian behavior of the suspension. [4] proposed general expressions for the induced stresses based on the suspension resistance tensors. According to [4], three physical mechanisms are responsible for the induced stresses: (i) hydrodynamic interactions between particles; (ii) Brownian diffusion; and (iii) forces over the particles. The procedure proposed by [4] provides the tools needed for the induced tensions determination from the particular mechanism that originates them. In that respect, the procedure allows to investigate the hydrodynamic irreversibility resulting from the [23] and the [7] mechanisms, but, in this case, from the induced stresses point of view. Later on, [3], following a similar procedure, obtained self-diffusion coefficients for suspensions of hard spheres under different flow conditions. They were able to differentiate among the different mechanisms the one responsible for the observed hydrodynamic irreversibility. For a shear flow, they proposed the following expression for the self-diffusion coefficient in the velocity gradient direction: where c is half the minimum distance between particles centers, which, for the case where particles forces are dominant, does not coincide with the particle radius, a. One could argue that in a sense, a and c are closely related to the [7] asperity concept. In addition, [3] suggest that their results are only approximate in contrast to those from [7], which are based on "the actual pair trajectory and, therefore, in principle, have determined the correct diffusivity." Zarraga et al. [26] determined the self-diffusivity through the full set of approaches for non-Brownian suspensions of spheres. They used the concepts of the excluded volume radius, c, the hydraulic particles radius, a, and its relation to the da Cunha et al. asperity [7]. Zarraga et al. [26] compared their self-diffusivity results with [7]'s ones and obtained (2) D s = 16 45 c 2 c very good agreement for low asperities. Furthermore, they extended the da Cunha's model and reported results up to a value of c corresponding to an asperity equal to 10. They, even, compared this high asperity results with their proposed analytical solution "in the absence of hydrodynamics." Under this condition, they found a full agreement between the values reported by the da Cunha's model and their analytical one. However, Zarraga et al. [26] state that their results are different from the ones reported in [3], Eq. 2, by a factor of 4 3 and that they cannot "account for this discrepancy." It seems, therefore, that the particles trajectories method becomes a valid approach in order to calculate the self-diffusion of non-Brownian suspensions of spheres.
Several theoretical investigations have been reported on the self-diffusivity of hard spherical particles, but very few, [20] and [21], aimed at nonspherical particles, which are dominant in the field of nanofluids. There is, therefore, a lack of knowledge in the effect of particles shape (nanofluids particles are not spheres) both on self-diffusion and apparent thermal conductivity for sheared suspensions. Up to the authors knowledge, there is no previous work determining the full effect (hydraulic and thermal) on kinematic irreversibilities due to particles collisions for suspensions of nonspherical particles. In addition, the new insight reported herein might provide useful information to take advantage in some other disciplines, such as [14].
The present paper reports the results of an investigation whose main contribution is to demonstrate that the selfdiffusion coefficient in dilute suspensions depends strongly on the shape of the particles and, for this type of particles, is independent on the Péclet number. The development of a model applicable to hard particles with three symmetry planes will be presented along with its application in the determination of the self-diffusion coefficient for equally sized cubic particles. For this purpose, the authors present a methodology to calculate the self-diffusion coefficient for cubic particles suspensions in the dilute limit and prove that, under this condition, the collision irreversibilities become the main source mechanism responsible for self-diffusion in sheared flows. The proposed methodology can be easily adapted, in future investigations, to calculate the apparent thermal conductivity following, for example, the proposals reported in [16]. It is also shown that, for the case of cubic particles, the self-diffusion is clearly enhanced with respect to the typically expected ones reported, for example, in [2,7,27].

Self-diffusion model
The present paper aims at developing a procedure to determine the self-diffusion coefficient of suspensions of nonspherical nanoparticles, though not elongated, characterized by three symmetry planes. One can think of polyhedral particles with an "equivalent spherical particle" of radius a. The dynamic performance of the proposed model is developed in this section based on the assumption that the trajectories of the polyhedral-shaped particles are identical as the ones corresponding to the equivalent spheres. The relative motion of a pair of particles will be assumed as similar to the one proposed by [1]. Though the relative displacement of the pair of particles has been previously studied by [7], the relative rotation rate of the particles has not been treated yet. Thus, the present model includes the determination of both the relative linear and angular velocities for any relative position of the pair of particles.
The determined relative spatial positions and rotation rate will be used in order to estimate the collisions rate when the particles are in the neighborhood of each other. The "neighborhood" is defined as the region of the space where the particles can physically collide with each other. For a regular polyhedron particle, it is defined as the minimum possible distance between centers when both particles are oriented in such a manner that their vertices are aligned along the straight line that connects both centers. This concept is somewhat related to the "excluded volume" of [3] and the "asperity" of [7]. For the moment, the term "neighborhood" will be assumed as being the particles interaction region. Under the proposed model, the self-diffusion coefficient is fully affected by the kinematic irreversibility sources (collisions rate) inside the interaction region.
Finally, the deflection of the pair of particles trajectories will be used to determine the self-diffusion coefficient. In a reversible Stokes pure shear flow, both particles approach each other in the compressional quadrant and the relative trajectory of one particle with respect to the other (assumed at origin of coordinates) is fully symmetric with respect to the three coordinate planes. If kinematic irreversibilities take place in the interaction region, the trajectories are no longer symmetric but rather they are deflected. The pair of particles will be assumed as being immersed in a pure shear flow in the proposed model, with the relative positions of the particles varying from far upstream to far downstream of one of the particles with respect to the other. The model allows the determination of the trajectories deflection and finally the self-diffusivity associated with pure sheared suspensions. As an example, the model will be applied in the determination of the self-diffusion coefficient of cubic-shaped particles.

Relative motion description
The relative motion of a pair of polyhedral particles is assumed as in [1] for spherical particles. A schematic view of the particles is shown in Fig. 1. One of the particles is located in the origin of coordinates, whereas the other moves around. According to [1], the relative motion of a pair of spherical particles in a Stokes flow is described by the following expressions corresponding, respectively, to the velocity and rotational fields In the above relations, V is the particles velocity vector, r the particles relative position vector, ijk the Levi-Civita tensor, Θ j the rotational rate component of a pure sheared flow, E ij the straining rate of the sheared flow, and is the particles rotation rate. Note that, since the Stokes flow is a linear one, both the velocity and the rotational fields are linearly dependent on the strain rate, E ij . Functions A, B and C are known as "mobility functions." These functions depend on the relative positions of both particles, the flow field and the shape of the particles. For a pure shear flow, U = ( y, 0, 0) , and hard spherical particles [1] developed the following expressions for the components of the relative velocity between the particles: where the subscripts correspond to the components in the main directions, x, y and z are the inertial coordinates located at the center of particle 1, r is the distance between centers of the particles, that is, r 2 = x 2 + y 2 + z 2 . The same set of equations was used by [7] in their evaluation of the paths of the particles needed to determine the self-diffusion coefficient. In addition, in the present model, the relative rotational speed will be evaluated in order to determine the overall relative motion between the particles. For a pure shear flow, the components of the rotational motion, Eq. 4, can be written as: The present model assumes that the mobility functions are the ones applicable to a pair of equivalent spheres of radius a, proposed originally by [1] and later on by [25] and [7], with the following expressions: where L = −Ln(r − 2). The integration of the above set of equations (5 to 10) can be performed numerically. In a time interval t , the relative displacement of the particles results from the superposition of a linear displacement, r , with an angular one related to the relative rotational speed, . Both the relative displacement and rotational speed are clearly illustrated in the sketch of Fig. 1.
The collision rate inside the interaction region will be determined based on the relative slip of the surfaces of the equivalent pair of spheres. The relative slip, D , originates from the superposition of both relative motions: the angular, , and linear one, whose displacement is designated as r . In order to determine the collision rate, it is necessary to realize that not all the relative motions can lead to a collision. The only motions resulting in net slip between the surfaces of both particles constitute an effective source for collisions. In order to make it clear to the reader, the reference plane has been used in the analysis. It is a middistance plane normal to the segment joining the centers of the particles. Note that when the second particle enters the interaction region, only relative motions contained on the reference plane will promote an eventual collision.
In order to obtain the displacement D on the reference plane ( Fig. 1), it is necessary to focus on the projection of the angular velocity, , on this plane. This component of is the only responsible for the mentioned displacements, D . D is contained in the reference plane too but normal to the projection of . The rotational rate in the radial direction does not promote net displacements in the reference plane and will not be considered in the collision rate determination.
The net displacement vector, r , can be decomposed into two main directions: the projection over the reference plane, D| [(r× r)×r] r , and the radial one. Obviously, the collision rate depends only on the projected component. As a result, on the reference plane there are two main directions: (i) (r × r) × r which correspond to the projection of the vector r and (ii) the perpendicular one ( r × r ). The next subsection is intended to determine the relative slip components along these directions.

Slip displacements on the reference plane
In what follows, the relative displacement between particles for each t resulting from translation and rotation will be evaluated through their projection in the reference plane along the main directions shown in Fig. 1. Since the rotation rate vector, , promotes slips perpendicular to its direction, it is necessary to identify the components of along the main directions in the reference plane shown in Fig. 1.
The direction normal to the linear displacement projection in the reference plane, r × r , is where the symbol means finite difference and ê i is a versor. The modulus is needed to determine the vector in this direction.
The dimensionless angular velocity projection along this direction is: Since both particles are identical, it can be stated that they rotate at the same angular speed due to the symmetry of the C mobility function. Thus, the relative angular velocity must be equal to twice the particle rotational speed. The relative displacement of the second sphere surface due to rotation in the present direction can be written as (see Fig. 1): If the dimensionless time and rotational speed are defined as: t * = t and * = 1 2 , Eq. 18 can be easily manipulated to obtain the dimensionless slip along the direction (r × r) × r.
where D stands for the slip distance increment, equal to in its dimensionless form, t is time, means the relative rotational angle between particles and the superscript * refers to dimensionless slip distance and time.
Proceeding in a similar way as in the previous case, one can determine the projection of the slip displacement in the r × r direction, which in the reference plane is given by the rotational rate projection along the (r × r) × r direction.
The dimensionless angular velocity projection along this direction is: The dimensionless slip promoted by the rotation rate along the direction r × r: Finally the linear displacement in the direction of the projection of the linear displacement in the reference plane, (r × r) × r , promotes the last relative slip displacement component that can be easily determined from the following general expression: (r× r)×r |(r× r)×r| r whose expression in terms of the relative separation coordinates is: Note that Eq. 25 is expressed in its dimensional form. Its dimensionless form can be easily obtained simply by expressing all coordinates related to the radius of the particles (equivalent sphere), a. Note also that two of the above displacements (the first one and the last one) occur in the same direction ( (r × r) × r ) and must be added algebraically. Since the second displacement is in the normal direction to the other two, the net dimensionless displacement can be obtained from the following expression: It must be stressed at this point that the net dimensionless displacement corresponding to each t could be considered as a relative rotation of one particle with respect to the other since the reference plane is parallel to the tangent plane of the equivalent spherical particles.
Since the relative rotation of both particles is a significant information to determine the potential collision between them, it is interesting to develop a model to determine the maximum allowable rotation previous to a collision in order to compare with the real rotated angle, , given by Eq. 26. The model will be complete if the maximum allowable rotation is determined in terms of the distance between particles centers. This maximum allowable rotation function will be designated as a "mean free angle" (mfa), which obviously will be defined in the interaction region.

Mean free angle and mean deflection model
The objective now is how to determine the intersection of two polyhedrons (particles), the center of one of them, particle 1, is located at the origin of coordinates and stands still, whereas the other, particle 2, whose center is located at r , is free to rotate. The impact will occur when particle 2 gets into the neighborhood of particle 1. The procedure for the intersection of the two polyhedrons, corresponding to the particles impact, is a standard one and will not be introduced here for space saving purposes. The impact could be obtained by following the trajectory of particle 2 and checking if the impact might occur at each relative displacement corresponding to a time interval t . This procedure, though plausible, would require a significant computer CPU time. Instead, an approximate procedure has been devised in order to determine the impact condition and its effects over the relative deflection of the particles. The procedure involves (25) D| [(r× r)×r] r = r 2 x 2 + y 2 + z 2 − (x x + y y + z z) 2 the determination of the so-called mean free angle, mfa, and "mean deflection," md, parameters that can be determined in advance, previously to the integration of the trajectories set of equations, since they depend on the particles geometry. In the development of the procedures for the determination of mfa and md, particle 2 will be assumed in the neighborhood of particle 1, and thus, in a condition that potentially could lead to impact.
• Mean free angle (mfa). The relative position of the particles is determined by the position vector of particle 2, r , and its orientation vector, expressed as . The particles impact could be monitored by a flag function for contact detection, F, depending on r and , and defined as: For a given position of particle 2, a "free of impact angle" could be devised as being the angle that particle 2 would rotate without impacting into particle 1. Since the orientation (angular position) vector components are where p is the number of times that F occurs over a complete rotation, 0 ≤ i < 2 , and = 1 , 2 , 3 is the angular relative orientation between particles. The average extended to all possible angular orientations is: Two particles can approach each other either at the near end or at the beginning of the ⟨ (r)⟩ interval. In the first case, the impact would be imminent, whereas in the latter the particle would rotate ⟨ (r)⟩ for the impact. It can be stated that none of both cases are significant since the expected value cannot correspond to these two extreme cases. The present model assumes that the most probable (expected) free of impact interval is half ⟨ (r)⟩ . The mfa is the average of half ⟨ (r)⟩ extended to any spatial orientation, r . Any versor, r , is associated with its differential area element, dS p , of the surface of particle 1 and therefore: Since mfa is an average value, the associated rootmean-square could be determined in order to estimate the expected error. S p represents the superficial area of the shaped particle.
Having introduced the mean free angle, the next step is to check for the impact occurrence if the particle rotates an angle in a time interval t . Since is a fraction of mfa, the ratio mfa can be considered as the probability of impact in the time interval t . A "cumulative impact probability" could be defined as the sum of the impact probabilities over the relative displacement of the particles so that the impact would occur when the cumulative probability would attain the unity. After the particles impact, a deflection occurs, the accumulated probability is set to zero, and the procedure is resumed. • Mean deflection (md).
The particles impact will cause a deflection whose extent must be determined. A function is introduced for that purpose. It is defined as the distance that the particles must separate after the impact to reach a nonimpact condition. The previously defined flag function, F, will again be useful in the determination of the function . This can be done by keeping a given initial orientation constant and checking for possible impacts for successive r increments. The procedure continues up to the point where no impact is obtained. Thus, can be considered as the difference between the modulus of the final position vector and the initial one (just after impact) for a given and spatial, r , orientations. A function G(r, ) , similar to F, which reports the required increment in the separation distance to avoid the impact, is introduced and defined as: According to the proposed model, the "mean deflection" is defined as the average deflection, , extended to all possible orientations. As for the case of the mean free angle, in this case, for a complete rotation of the angular coordinate 3 , the average deflection, representing the average required increment in the separation distance in case of contact for a full impact period, is given by the following expression:

Page 9 of 18 392
The md function is obtained as the average extended to all possible spatial and angular orientations. The procedure is analogous to the one for mfa based on .

Application of the proposed model to cubic particles
The objective of this section is to present an example of application of the proposed procedure in the determination of the mean free angle and the mean deflection for cubic particles. In order to do so, the values of the functions F and will be determined based on a discrete number of relative positions.
Let particle 1 be centered at the origin of coordinates and the center of the particle 2 at a position r , as shown in Fig. 2a. Particle 1 is considered to be at rest with its faces oriented along the coordinates system, whereas particle 2 is free to rotate. Particle 2 angular orientation, , is related to a local coordinate system, XYZ, displaced r from the general one, xyz. As a result, each cube vertex can be located in space by a position vector, r , and an orientation one, .
Due to symmetry, the cubic particle full rotation ( 2 ) is not required in the function F determination. Instead, a rotation in the range ± 4 is enough in this case. The determination of the deflection, , is based on a geometric series iterative procedure in order to reduce the related uncertainty. If the impact has been confirmed, the minimum noncontact distance (final r) can only assume a value between the given radial distance r and the outer edge of the interaction region.
For a pair of cubic particles, the maximum r value is equal to 3 1 2 b, which corresponds to the particles in a position such that two opposite cube vertices are in contact. Obviously the direction in space,r , and the orientation, , must be kept constant over the contact check procedure. The minimum noncontact distance will be determined by checking the contact on the extreme of geometric series intervals. In the first iteration, contact occurs in the inner radius though not at the outer one. The next check is performed assuming the new check point as the mid one between the inner and outer radius of the previous iteration. Two possible scenarios will result at this point: (i) If contact does occur at this new location, it is assumed as the inner radius of the next iteration, whereas the outer one is the same as the previous one. (ii) In case of no contact, the new location is assumed as the outer radius, whereas the inner one remains the same of the previous iteration.
The iteration procedure will go on up to the point where the chosen precision is attained. It is interesting to note that the difference between radial distances after 9 or 10 iterations is of the order of 0.1% of the maximum possible range, 3 1 2 b. The number of iterations assumed in all the cases considered in the present study has been set equal to 9 unless otherwise specified.

Discretization of the angular, , and spatial r orientations.
The proposed procedure for the determination of mfa and md described in previous section requires integration over the orientation and separation vectors. This is done by a discretization procedure that will be described in what follows for cubic particles. Two particles close to each other are illustrated in Fig. 2a. It must be noted that due to symmetry of cubic particles, the direction can be set using just one of the faces of the still cube, for example, the upper surface. Again, due to symmetry, one can work with just one of the quadrants of this surface, for example, the positive one ( x ≥ 0 ; y ≥ 0 ). This quadrant is divided into I points in the x direction and J in the y one, as shown in Fig. 2b. Each of the points in each axis is designated, respectively, by i and j such that 1 ≤ i ≤ I and 1 ≤ j ≤ J . By joining a point on the surface, given by the pair (i, j), with the center of the cube different directions of the r vector are defined. Having defined the directions, it remains to characterize the modulus of the separation vector, r. It must be noted that in the socalled interaction region, the r varies in the range from b to 3 1 2 b. The separation distance r is discretized by dividing it in K points each one designated by k so that 1 ≤ k ≤ K . Thus, given the above discretization procedure, the separation vector can be expressed as: The dimensionless components of the position vector are given by the following expressions: where i, j, k are indexes related to the relative particles orientation discretization, I, J, K represent the total number of the relative orientation positions in the discretization, X, Y and Z are the local coordinates of the center of particle 2, and b is the side size of the cubic particle. In the case of the orientation vector, the rotation angle is divided in L points.
The procedure for mfa and md determination consists in evaluating the average of all directions, defined by I and J, and angular positions, L, for each value of k, corresponding to a given value of r. For each value of k, the average procedure is repeated.

Application of the discretization procedure to the F and « functions
As an example of application of the discretization procedure described above, functions F and will be determined as a function of the rotation of particle 2 for the X direction, X . The number of points in the upper face of particle 1 was set as I = J = 3 , so that only 9 pair of points are considered. The number of points in the separation vector has been set equal to 20, that is, K = 20 . The separation vector is thus determined from Eqs. 34, 35 and 36. As mentioned above, only the X component of the orientation vector of particle 2 will be considered, the other two being set equal to 0. This simple set of orientations allows to easily determine the particles relative location in space and their F and symmetries. The results of the discretization procedure are shown in Fig. 3 where F and functions are overlaid in the plot against X , for different pairs of points (i, j) of the upper surface of the still particle, characterizing different directions of r . The corresponding distances between centers for the different orientations of Fig. 3 are, respectively, the following: k = 4(i = 2, j = 1) ; k = 8(i = 1, j = 2) ; k = 8(i = 2;j = 1) ; and k = 9(i = 2;j = 2) . A total of 31 angular positions, corresponding to an angular step of 3 degrees, have been assumed for the results of Fig. 3. Results for directions corresponding to the edges of the superior plane, i = j = 3 , of particle 1 have not been displayed in this figure due to their straightforwardness given that the rotation of particle 2 is with respect to the X axis. The deflection function, , has  -The assumed symmetry with respect to 2 is confirmed. -As should be expected, the flag function and the deflection are asymmetric with respect to 0 degree rotation angle. -Depending on the particles relative vector position, the deflection function asymmetry does not follow the same pattern than that of the flag function. In that respect, it is interesting to note that the results for the pairs ( i = 1 , j = 2 ) and ( i = 2 , j = 1 ), for equal particles separation ( k = 8 ), present a shift in the pattern between them.

Mean free angle and deflection function results
F and have been determined according to the discretization procedure described in the preceding section. Thus, the evaluation of mfa and md will depend on the number of discrete points used in the determination of F and . The plot of Fig. 4 has been raised to check for the effect of the number of discrete points and, based on the obtained results, set this number aiming at obtaining reasonable precision in the mfa and md determination. In the results shown in Fig. 4, three different pairs of discrete points were used for the position vector, r , characterized by the values of I and J. As for the angular position vector, , three different values of L have been used to cover the range ± 4 . The number of shells, K, has been set equal to 10 in the results of Fig. 4, corresponding to the interaction region given by 1 ≤ r b ≤ 3 1 2 . From Fig. 4, it may be concluded that mfa and md do not depend on the chosen set of discrete points. Thus, in what follows, the number of discrete points will be set as I = J = 5 , K = 10 , and L = 9 , in order to optimize CPU time/memory. It could be argued that the mfa and md functions have lost their physical meaning since they are the result of an average procedure. However, the average value loses its meaning when the root-mean-square, rms, of the considered parameter is high. The plots of Fig. 5 have been raised to check for the adequacy of the mfa and md functions obtained from the proposed procedure through the root-mean-square of each of these parameters. The uncertainty bars of these plots correspond to the average value plus or minus the rms of each parameter, for the number of discrete points previously suggested for the present study. It can be noted that the mfa presents significant uncertainties in the mid of the range, whereas toward the ends the uncertainty tends to zero. In the case of the md, the uncertainty increases toward the inner limit.

Results
This section will focus on an analysis of the numerical results from the proposed model applied to cubic particles. The analysis will be focused initially in evaluating the model and its sensibility to certain intervening parameters. In the second part, model results will be applied in the determination of the self-diffusion coefficients of dilute suspensions of cubic particles. It must be stressed that the adopted procedure will be based on the one proposed by [7] for the determination of both the particles trajectories and the diffusion coefficients.  1 1.2 1.3 1.4 1.5 1.6 1.7 1 I=J=5 L=9  I=J=5 L=9  I=J=5 L=15  I=J=5 L=15  I=J=5 L=20  I=J=5 L=20  I=J=7 L=9  I=J=7 L=9  I=J=9 L=9  I=J=9

Model evaluation
As previously noted, the analysis involves the interaction of a couple of particles and their relative trajectories. The particles are located in a pure shear flow field of the kind U = ( y, 0, 0) with the initial relative vector position of the particles being the same as the one adopted by [7] in their analysis. The particles are assumed to depart from a far upstream plane in the compressional quadrant and perpendicular to the main velocity, U , located 10 dimensionless units from the origin of coordinates, r a = −10, y 0 , z 0 . In order to analyze the model performance, one single-particle trajectory departing from r = (−10, 0.1, 0.1) will be reported, as in [7]. The components of the position vector given above are related to the spherical radius, a, whose adequate specification will be described further on. The time step and the integration of the trajectories equations, by a fourth-order Runge-Kutta algorithm, will be the same as the ones of the paper by [7]. In the far field, r a > 2.5 , t * = 0.01 ; in the intermediate region 2.01 ≤ r a ≤ 2.5 , t * = 0.005 ; and in the lubrication region 2.0 < r a < 2.01 , t * = 0.0001 . Note, from the results shown, that the pair of particles separation are, in any case, above 2.4-2.5 [-], and therefore, they can never enter into the lubrication region.
The first set of results has been obtained for cubic particles of side equal to b with the corresponding spherical particle being that of radius a = b 2 , that is, the sphere is circumscribed by the cube. The corresponding functions mfa and md are the ones plotted in Fig. 4. The projections of the trajectories are reported in the plots of Fig. 6. It can be noted that the displacement in the y direction is higher than in the z one. It can also be noted that once the particles get in touch, successive impacts occur characterized by the ripples in the curves which continue up to the point of maximum deflection. Beyond that point, no contact is observed. Most of the impacts occur in the compressional quadrant, x < 0 , though a couple of collisions can be noted in the extensional quadrant, x > 0 . The latter collisions are clearer in the plots of 7, where the projections of the rotational speed on the reference plane, t and l , are plotted along with the dimensionless separation distance, r 5a , against the longitudinal position of the particles, x a . Note that t and l correspond to the components on the reference plane of the rotational speed in the (r × r) , Eq. 16, and (r × r) × r , Eq. 22, directions. The successive impacts are clearly shown in these plots. The most influential component of the rotational speed is t for the present particles. As a result, the component is also the dominant one in generating irreversibilities. With respect to the separation distance, two aspects are clearly noted in the plot of Fig. 7: (i) the seesaw behavior corresponding to the aforementioned successive impacts and (ii) the minimum separation distance between the particles, corresponding to the minimum of the curve, which in this case is equal to 2.5 radius of the circumscribed sphere. Regarding the results from Figs. 6 and 7, it can stated that the observed trends comply with the expected physical insight of the problem though the obtained trends are clearly different from the ones obtained by [7] for their highest asperity values.
A second sensibility analysis could be performed in connection with the mfa and md plots of Fig. 5 and their associated root-mean-square values. This analysis has to do with the maximum and minimum irreversibility generation. In that respect, as should be expected, the irreversibility not only increases with the observed deflection, characterized by the value of mean deflection, md, but also with the number of impacts of the particles which increase for lower values of the mean free angle, mfa. Thus, the results reported in Fig. 5 suggest two limit conditions, superior and inferior, with respect to the average, corresponding to the maximum and minimum irreversibility generation, respectively, whereas the average condition corresponds to the results of x/a [-]  The plots of Fig. 8a, b illustrate the effect of the variation of the pair (mfa, md) from its average to its superior and inferior values, respectively. The superior pair causes an increment in the deflection with the opposite effect being noted for the inferior pair, as should be expected. Similar effect is noted in the case of the separation distance, as shown in Fig. 8a. It is interesting to note that the deflection effects of Fig. 8b could be important in characterizing the uncertainty range in the determination of self-diffusion coefficients based on the proposed model.
A third sensibility analysis could be devised in relation with the effect of the mobility functions, A, B and C, since in the present case they were assumed as those for spherical particles from the study by [1]. Though not shown here, the obtained results clearly indicate that suppressing function C, in other words, setting it equal to zero, does not affect the relative deflection of the particles. This could be explained by the fact that the main rotation effect is the one related to the pure shear flow assumed in the present study. The effect of functions A and B will be investigated considering different relations between the sizes of cubic shape and the corresponding sphere. As an initial statement, it must be stressed that the use of the assumed mobility functions for nonspherical particles is questionable. However, given the unavailability of adequate functions, their use must be considered an approximation in the self-diffusion evaluation. Similar arguments could be used regarding the model proposed by [7], and even the one by [2], who assumed a spherical particle of radius a, with the so-called asperities being external to the spherical surface. From the conceptual point of view, this approach to the problem is also questionable since the fluid cannot penetrate into the sphere surface given the presence of the asperities. As previously mentioned, the effect of the mobility functions A, B and C will be evaluated by varying the relative size of the equivalent sphere with respect to the cubic particle. The procedure begins by assuming that the radius of the sphere is half the side of the cubic particle. Now, the side of the cubic particle will diminish with respect to the sphere radius from the condition where the sphere is circumscribed in the cube, a = b 2 , to the one where the cube is circumscribed in the sphere, a = 3 1 2 2 b = 0.866b . The former corresponds to a geometrical relation between the sphere particles and the cubic ones with a maximum irreversibility due to the number of particles collisions, whereas in the case of the latter relation there is no irreversibility since the particles are effectively spheres. The general conclusions previously drawn for the maximum irreversibility case will hold for cases for which the irreversibility diminishes down to the reversible case.
The trajectories for different a b ratios are overlaid in the plot of Fig. 9a, corresponding to the following particular cases: (i) a b = 0.5 , sphere circumscribed in the cube, the reference case (ii) a b = 0.62 , sphere of equal volume (iii) a b = 0.74 and (iv) a b = 0.866 , cube circumscribed in the sphere, the reversible case. Note that the coordinates of the plot are referred to the spherical particle radius which is constant for all cases, whereas the cube side is diminished for each case down to the reversible one. Using the terminology of [7], each case could be considered one of different asperity, characterized by the ratio a b , with the latter being the smooth surface sphere one. The location y a of the particles for each case far downstream x a = 15 tends to differ significantly between each other, with the one corresponding to the higher asperity a b = 0.5 being the most distant from the original location far upstream x a = −10 . This trend is a clear indication of the occurrence of irreversibility. The particles separation distance for each case is shown in the plot of Fig. 9b. From this plot, it can be concluded that the trajectory corresponding to the case a b = 0.74 is not physically realistic since the closest distance between particles for this case is almost twice the spherical particle radius. This would correspond to an approximation of the particles for which the lubrication regime would hold, a condition that is not physically attainable with cubic particles. The mobility functions are dominant in the flow in such cases, and as a result, they should be discarded as physically unsound. As a final note regarding the plots of Fig. 9, one must recognize that the sphere radius in the mobility functions has no physical meaning since it is just a reference size used in their evaluation. Thus, it would be interesting to check what the trajectories would look like if instead of being referred to the sphere radius, the coordinates were referred to the half side of the cube, b 2 . This is shown in Fig. 10 for the same cases of Fig. 9. By doing so, the result is that the trajectories tend to collapse into a single curve in the compressional quadrant, whereas their differences tend to diminish in the extensional quadrant, as can clearly be noted in the plot.
As previously suggested, results from cases (i) a b = 1 2 corresponding to a circumscribed sphere, and (iv) a b = 0.866 , for a circumscribed cube, illustrate limiting cases of the proposed model and do not reproduce realistic patterns. The two intermediate cases are (ii) a b = 0.62 and (iii) a b = 0.74 that, as shown in Figs. 9 and 10, display a different pattern. However, note that both figures show that, for a b = 0.74 , the minimum distance between particles is attained over a wide region, with almost no impacts, which does not seem realistic for cubic particles. On the other hand, note that in terms of dimensionless cube size (Fig. 10b) the case a b = 0.62 keeps almost the same distance between particles than for the a b = 0.5 case, though with a number of impacts significantly lower, which seems to constitute a more plausible . The x distance is referred to the cubic particle half size, b pattern. As a result, from the observed trends in Figs. 9 and 10, one can conclude that the optimal case in terms of the sphere/cube size ratio is the one corresponding to the equivalent sphere (equal volume) which is the one with a reliable physical meaning.
Based on the sensibility analysis, one can conclude that the coefficient of self-diffusion must be evaluated from particles trajectories obtained from mobility functions based on the equivalent spherical particles. The uncertainty in their evaluation must be estimated based on the root-mean-square values from Fig. 5.

Self-diffusion coefficients
Self-diffusivity is related to the net flux in a given direction. Thus, it is determined in terms of the streamline deflection of the relative trajectories of a pair of particles [7]. The center of mass of a pair of particles runs along a "streamline." The deflection of the mass center of a pair of particles in the direction i of space, H i , is related to the trajectories by the following expression: The superscript ±∞ denotes the far downstream and upstream trajectories coordinates, respectively. The trajectories are considered to be at x +∞ (far downstream) and x −∞ (far upstream) if x a ≥ 10 and x a = −10 , respectively. Consider a uniformly diluted suspension where different deflections are not correlated with each other. It can be stated that the collision rate between a particle at rest, located at the origin, and the entering particles far upstream at any given y −∞ is: n |y −∞ |dy −∞ dz −∞ , where n is the number of particles per unit volume. Therefore, the self-diffusivity can be expressed as [7]: Note that = nv p , where v p is the volume of a particle. Therefore, Eq. 38 dependence on n can be expressed in terms of . For the case of spherical particles, the dimensionless form of Eq. 38 assumes the following expression: where H * i sphere = H a , and y * sphere = y a . Similarly, the dimensionless self-diffusivity for a suspension of cubic particles can be written as: where H * i cube = 2 H b , y * cube = 2y b , and b is, again, the cubic particle side, ±∞ refers to far downstream and upstream, respectively, and H is the streamline deflection. Equations 39 and 40 can be related as: In the preceding section, it has been assumed that the equivalent spherical particle of a cubic one is the equal volume sphere. Therefore the 2a b relation is known, and as a result, Eq. 41 can be written as: The set of Eqs. 37, 39 and 42 allows the determination of the self-diffusivity (coefficient of self-diffusion for a given direction) of cubic particles in a dilute suspension. The procedure is directly related to the equivalent spheres suspension previously proposed by [7].
The streamlines at far downstream, x +∞ ≥ 10[−] , are depicted in Fig. 11 for departures from the plane x −∞ = −10[−] , previously indicated. The axis is in dimensionless form referred to the cubic particle side, b, for the three irreversibility cases corresponding to the average, superior and inferior cases. Figure 11a reports the results of the superior limit of the irreversibility generation, whereas Fig. 11b reports the inferior limit. Finally, the mean expected streamlines are shown in Fig. 11c. The displacement of the stream lines is clearly visible in this figure, with the displacement in the z direction being clearly inferior to the one in the y direction for the three cases considered.
The self-diffusivities are determined for the three cases shown in Fig. 11. The results provided for the superior, average and inferior kinematic irreversibilities are: 0.326; 0.233; 0.115 which correspond to dimensionless self-diffusion coefficient of D s y b 2 2 = 0.233 ± 45% . It must be noted that this is the self-diffusion coefficient, in the y direction, for a dilute suspension of cubic particles. It is interesting to observe that self-diffusion coefficient thus obtained is applicable to any shear rate and, as a result, is not associated with a particular Pe number. Furthermore, for a suspension of spherical particles in the absence of hydrodynamics, [3] found that the dimensionless self-diffusion coefficient, given by D s c c 2 = 16 45 , is equal to 0.113, a value that, though different, is unexpectedly close to the result of the present investigation. [26], for the case of absence of hydrodynamics, reports a value of 0.085, lower than the one proposed by [3]. The procedure proposed by [26] assumed that = 2 (c−a) a . Let the radius of the interaction region be the hard spere radius of [26], c = 3 The following three main conclusions can be drawn from the previous comparison: (i) The absence of hydrodynamics model for spherical particles can be correctly adapted to cubic particles in order to determine the self-diffusion coefficient; (ii) at least for cubic particles, it can be stated that the "equivalent" spherical particles size, a b = 0.62 , adopted to model the pair hydrodynamics seems to work properly to report the self-diffusion in the absence of hydrodynamics; and (iii) the collision pattern can be considered realistic since the self-diffusion is dependent on how, when and where the collisions take place.
Finally, it can be stated that the results provided by the proposed model demonstrate that the hard particles shape does play a role in the self-diffusion value up to the limit of the "absence of hydrodynamics."

Conclusions
A non-Brownian of hard cubic particles suspension in the dilute limit has been studied in order to obtain the selfdiffusion coefficient under pure shear rate. The proposed model considers the impacts between particles as sources (c) Average of kinematic irreversibilities in the Stokes flow. The collisions are studied in a statistical manner for the particular cubic shape through the so-called mean free angle and mean deflection functions and the slip mechanisms between particles. The following are the main conclusions drawn from the reported investigation.
-The mechanism responsible for the kinematic irreversibilities is the main rotational rate produced in the shear flow. Small changes in the separation of the particles inside the interaction region result in significant change in the number of impacts. -The model is quite sensitive to the root-mean-square of the mfa and the md functions. Displacements of mfa curve toward the inner or outer zones of the interaction region affect the self-diffusion values causing an increment of the uncertainty related to the root-mean-square of this parameter. -In a suspension of cubic particles, kinematic irreversibilities are generated by the particles shape alone. The selfdiffusion is strongly affected by the particles shape which can become the main source for the self-diffusivity. This issue has been previously proposed by [2]. [3] reported similar self-diffusivity for suspensions where the collisions between particles become dominant, though they linked this behavior to high values of the Péclet number, 1 ≪ Pe < ∞.
More research is under development to obtain a deeper insight into the nonspherical particles migration in dilute suspensions, considering particles of other formats. The reported results will be applicable to macroscopic CFD models which include self-diffusion as a potential particles diffusion in continuum mechanics [10].