it ro a, Accepted 26 March 2013 Available online xxxx ear such as Ladd and ALD methods, viol ate conservation law of the translational and rotational momentum mental methods are expensive and cannot obtain details of the tional velocities. Analytical methods are limited to systems with a small number of particles and simple geometri es. Due to limitations of the experimental and analytical methods , numerical methods have been used extensively which made a sig- nificant contribution to our understand ing of suspension mechan- ics and rheology. The most widely used method for simulating omplex boundary ds in simulating hree dimensions. e irregula on the te particle position. tracted much attention as a promising alternative for simulation of fluid flows with complex physics such as multiphase and multi com- ponent flows, simulation of particulate suspensions and colloid hydrodyn amics [2,3]. LBM is a method based on the solution of the Boltzmann equation on a lattice with a discrete velocity field. Several issues explain the usefulness of LB method in particulate suspensi ons [4]. First, LBM provides straightforwar d and easy-to- impleme nt boundary conditions necessary to couple the fluid to the suspended particles. Second, the method is computationall y ⇑ Corresponding author. Advanced Powder Technology xxx (2013) xxx–xxx Contents lists available at d w. E-mail address:
[email protected] (S. Jafari). dynamical behavior of each particle such as translationa l and rota- In the last two decades, lattice-Boltzma nn method (LBM) has at- lyzing and understanding the transport of solid particles in the fluid, movement and dynamics of these particles should be thor- oughly investigated . Experime ntal, analytica l and numerical methods were used for investigatin g dynamics of particles suspended in fluids. Experi- main are needed which is not a simple task for c condition s and geometries. Using these metho dense flows has limited success, especiall y in t The main difficulty is that reconstructions of th are necessar y at every simulation time step based 0921-8831/$ - see front matter � 2013 The Society of Powder Technology Japan. Published by Elsevier B.V. and The Society of Powder Technology Japan. All rights http://dx.doi.org/10.1016/j.apt.2013.03.018 Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lees–Edwards boundary conditions with smo othed profile-lattice Bolt method s to introd uce shear into particle suspension s, Advanc ed Powde r Technolo gy (2013), http://dx.doi.org/10.1 016/j.a pt.2013.03.018 r mesh mporal Particle suspensions play an important role in many industrial and biological situation s. Even though the particles may be sus- pended in a Newtonian fluid, the resulting suspensi on exhibits much more complex and non-Newton ian flow characteristics. Understandi ng the complex rheology of these suspensi ons is important in the optimization of industrial processes, reduction of cost and waste, and in the biological case, saving lives. For ana- are simultaneou sly solved at discrete time steps for all the particles in the domain. However, it is extremely difficult to deal with dense suspensi ons and suspensions consisting of non-spherical particles by means of Stokesian dynamics due to its complicated mathemat- ical structure s. In conventional numerica l methods, such as finite volume and finite element methods, a proper grid for the computational do- Keywords: Lees–Edwards boundary condition Lattice Boltzmann method Smoothed profile method Suspension Shear 1. Introductio n (Galilean invariance ). In the present study, Lees–Edwards boundary condition has been combined with smoothed profile method (SPM) intending to eliminate Gal ilean invariance errors. The combined method is validated by allowing a particle to cross a Lees–Edwards boundary. Moreover, third-order interpolation is used for particle distribution functions leaving the domain in the velocity gradient direction to elimi- nate bumps in the angular velocity of the particle when crossing the Lees–Edwards boundary. As another test case, two interacting circular cylinders placed in a sheared domain using Lees–Edwards boundary condition. Comparing results with the ones presented in the literature shows good agreement. � 2013 The Society of Powder Technolog y Japan. Published by Elsevier B.V. and The Society of Powder Technology Japan. All rights reserved. suspensi on flow at creeping flow conditions has been Stokesian dynamics [1] in which linear equations describing Stokes flow Received 7 December 2012 Received in revised form 15 February 2013 ties of suspensions. Employing Lees–Edwards boundary conditions as an alternative method, removes these effects. Earlier methods of solid–fluid interactio ns in the framework of lattice Boltzmann method, Original Research Paper Combining Lees–Edwards boundary cond profile-lattice Boltzmann methods to int suspensions Ebrahim Jahanshahi Javaran a, Mohammad Rahnama a Department of Mechanical Engineering, Shahid Bahonar University, Kerman, Iran b Department of Petroleum Engineering, Shahid Bahonar University, Kerman, Iran a r t i c l e i n f o Article history: a b s t r a c t Using walls to introduce sh Advanc ed Pow journal homepage: ww ions with smoothed duce shear into particle Saeed Jafari b,⇑ into a domain causes wall effect s in the calculation of rheological proper- SciVerse ScienceDi rect e r Technology elsevier .com/locate /apt reserved. zman n Po efficient and scales linearly with the number of particles, N. The local nature and the regular lattice grid facilitate the use of LB method on distributed-m emory clusters, in contrast to competit ive numerica l methods such as Stokesian dynamics. Also, the dense and transient nature of suspensions limits the benefits of non-uniform grid tech- niques that would be commonly used in numerical schemes such as finite-element analysis. In general the solid–fluid interface is usually subject to the no-slip boundary condition, which requires the local fluid velocity at the boundary of the particle to be the same as that of the fluid. Since in the LBM, the unknowns are particle distribution functions, the no-slip boundary condition needs to be expressed in terms of the particle distribut ion functions. The most widely used and sim- plest way to implement no-slip boundary condition for LBM is the standard bounce-back method [5,6]. The main drawback of this method when applied to particle–fluid interaction problems with complex shapes is the step-wise representat ion of the particle boundary which causes fluctuations on the computation of forces exerted on the particle. Consequently, significant efforts have gone toward improving the accuracy of the no-slip boundary conditions; these improvements are broadly split into two groups. The first group retains its basis on links crossing the solid boundary and re- lies on interpolation schemes [7] to increase the accuracy when the solid does not lie at the midpoint of the link (link-based methods ). However, with the use of the single relaxation time (SRT) collision operator, there also exists an error in the boundary location that is dependent on the viscosity which can be reduced by using two relaxation (TRT) or multiple relaxation time (MRT) lattice Boltzmann methods. In the second group, which are called node-bas ed or immersed boundary methods, fluid-particle interactions are simulated by using a combination of the immerse d-boundary method (IBM) and the LBM [8]. In these methods, a Lagrangian grid is used to track the particle’s surface in addition to the Eulerian LBM grid. The points on the Lagrangian grid are allowed to freely advect with the local fluid velocity. Then, a solid model is used to calculate a restorative force term on the Lagrangian nodes, which is typically a high-stiffne ss elastic model for the simulation of rigid particles. As in traditional IB techniques, the force term is applied to the fluid as a force density through an appropriate weight function. In the LB method, this force density is prescribed through a forcing term to the LB equation. To avoid complicated tasks of calculating the force density, a different direct forcing method is introduced based on the smoothed profile method (SPM) [9]. It uses a fixed Eulerian grid for the host fluid and represents the particles by certain smooth body forces in the Navier–Stokes equations instead of treating the particles as boundary condition s to the fluid. SPM solves a single set of fluid dynamics equation s in the entire domain including the particle volumes, without any internal boundary conditions. A smoothly spreading interface layer is used to repre- sent the particle boundaries to give a transition from the rigid- body motion to the fluid motion. Using this simple modification, regular Cartesian coordinates can be employed for many particle systems with any arbitrary particle shape, rather than boundary- fitted coordina tes. The solid–fluid interface has a finite volume supported by multiple grid points. Thus, a round particle shape can be treated in a fixed Cartesian grid without difficulty. Further- more, the computational demands for this method include sensi- tivity to the number of grid points (volume of the total system); however, it is insensitive to the number of particles. Due to the common feature of using Cartesian mesh in the LBM and SPM, Jafari et al. [10] combined these methods and applied it to partic- 2 E. Jahanshahi Javaran et al. / Advanced ulate suspension and reported its effectivenes s in such flows. In the present work, this method is used for simulation of particle suspension in a fluid. Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lee method s to introd uce shear into particle suspension s, Advanced Powder Techn By far, the most widely used flow to probe the rheology and structure of suspensions is shear flow [11]. To this end, the system undergoes a shear by placing parallel walls at the top and bottom of the domain, which move in opposite directions. As it is known, effects of the wall such as formation of a depletion zone near the wall and wall slip [12], on the rheological propertie s of the system is not prohibited unless the size of the system is sufficiently large which is not computationall y efficient. To reduce wall effects and obtain a feasible simulation, Lees–Edwards (LE) boundary condi- tions are used to introduce a shear into the domain. This method, proposed by Lees & Edwards (1972) for molecular dynamics simu- lations [13], has been extended to the LBM by Wagner and Pagonabar raga [14]. In LE boundary condition, the flow and vortic- ity directions are treated in the typical periodic manner; however, in the shear direction, periodic domains are shifted continuo usly in time. Lorenz and Hoekstra [15] applied Lees–Edwards boundary condition for the lattice Boltzmann suspension simulatio ns. They showed that standard methods for lattice Boltzmann simulations of suspended particles , which are based on the momentum exchange algorithm, might lack accuracy or violate Galilean invari- ance, especiall y when the particles cross the LE boundary. They proposed corrected momentum algorithm [16] and applied it for a particle going across a LE boundary. Zhou et al. [17] combined Lees–Edwards boundary condition with immerse d moving bound- ary (IMB) conditions and demonstrat ed that their combined meth- od conserves the law of the translationa l and rotational momentum . In all above mentioned studies, first-order interpola- tion was used for particle distribution functions leaving one LE boundary and reinserted from the other LE boundary. This causes some bumps in the angular velocity of the particle when crossing this boundary . In link-based suspension methods when the particle crosses the LE boundary, the particle is partly described in one reference frame and the other part in another reference frame with a different underlyin g velocity. This behavior is not present in node-bas ed suspensi on methods. In smoothed profile method, which is a node-bas ed method, in order to satisfy the no-slip boundary condi- tion, a body force is introduced inside the solid domain to enforce the fictitious fluid inside this particle to satisfy the rigid-body mo- tion constraint. Using combination of SPM, with local interaction between solid and fluid or without required link between the solid and surrounding liquid, and Lees–Edwards boundary condition conserves Galilean invariance. In combined SPM–LEBC, the part of the particle that leaves the upper boundary will be reinserted from the bottom boundary at the proper location where the parti- cle distribution functions and fluid properties are previously affected by the LE boundary condition as in the case for LE bound- ary condition for fluid-only systems. In addition, the properties of the leaving part of the particle are corrected before reinsertion from the bottom boundary to match the fluid properties at the bot- tom boundary at particle location. The effect of violation the Galilean invariant on the rheological properties of the suspensi ons has been investigated by Clausen and Aidun [21]. In the current study, smoothed profile method which is a node- based method is combined with the Lees Edwards boundary condi- tion to reduce or eliminate Galilean invariant errors arising in the link-base d suspension methods when the particle moving in the sheared domain and then crosses the LE boundary. Moreove r, the violation of conservation law of the translational and rotational momentum due to some bumps in the angular velocity of the LE boundary crossing particle which is observed in previous investi- gations is resolved by using third-order interpolation for particle distribut ion functions leaving this boundary. wder Technology xxx (2013) xxx–xxx This paper is organized as follows. In Section 2, we present a description of LBM and SPM. Section 3 is dedicated to the Lees–Edwards boundary condition for fluid-only and particle s–Edwards boundar y conditions with smoo thed profile-lattice Boltzman n ology (2013), htt p://dx.doi.org /10.1016/ j.apt.2013.0 3.018 controls the relaxation of the viscous stress in the fluid and is con- is introduced as this particle to satisfy the rigid-body motion constrain t which is zero outside the solid domain.The fluid–solid interaction force /fP acting on the solid particle nodes is given by [10]: /ðx; tnÞfPðx; tnÞ ¼ /nðuPðx; tnÞ � uðx; tnÞÞ=dt ð8Þ In the above equation, uPðx; tnÞ and uðx; tnÞ are the particle and fluid velocities at t ¼ tn and x, respectivel y. The resulting force act- ing on the fluid nodes located inside the solid particles is given by: fHðx; tnÞ ¼ �/ðx; tnÞfPðx; tnÞ ð9Þ As mentioned before, in SPM only one equation is solved in the entire domain and the effect of solid particles is accounted for using a body force term. In LBM, there are several methods for introducing a body force term into the evolution equation . The most common one is adding a term into the collision operator as follows: faðx þ eadt; tn þ dtÞ ¼ faðx; tnÞ � 1s ½faðx; tnÞ � f eq a ðx; tnÞ� þxadt c2s ðfH:eaÞ ð10Þ d Po /ðx; tÞ ¼ XNP i¼1 /iðx; tÞ ð5Þ nected to the kinemati c viscosity via m ¼ c2s dtðs� 1=2Þ in which cs is the speed of sound and is defined by c2s ¼ c2=3, where c = dx/dt and dx is the lattice constant. The most widely used lattice model in two dimensions is the D2Q9 model. In this model, the discrete velocity vectors are e0 = (0, 0), ea ¼ ðcos ða� 1Þ p2 � � ; sin ða� 1Þ p2 � �Þc for a = 1–4 and ea ¼ cos ð2a� 9Þ p4 � � ; sin ð2a� 9Þ p4 � �� � ffiffiffi 2 p c for a = 5–8. The stream- ing speed, c, is defined as dx/dt. The equilibrium distribution func- tion, f eqa , is chosen such that the Navier–Stokes equation s for a weakly compressible system are obtained and for the D2Q9 model is of the form f eqa ¼ xaq 1 þ ea:uð Þ=c2s þ ðea:uÞ2=2c4s � u2=2c2s h i ð2Þ where xa is a weight factor which is 4=9 for a ¼ 0, 1=9 for a = 1–4 and 1=36 for a = 5–8. In the discretize d velocity space, macrosco pic fluid density and velocity can be evaluated as q ¼ X8 a¼0 fa ð3Þ u ¼ 1 q X8 a¼1 eafa ð4Þ 2.2. Simulation of fluid–solid interaction by the smoothed profile method (SPM) In SPM, the surface of the particle is treated not as a sharp inter- face having no thickness, but instead, an interface is introduced having a width compara ble to the grid spacing [9,18]. In other words, SPM represents each particle by a smoothed profile, which equals unity in the particle domain, zero in the fluid domain, and varies smoothly between one and zero in the solid–fluid interfacial domain, see Fig. 1. Here, field quantities such as velocity and particle density field are defined on the entire computati onal domain, which includes not only fluid regions but also regions in which particles exist. In order to determine regions in which particles exist, a function, /, suspension systems. In Section 4 validation test and numerica l re- sults for one particle and two particles undergoing shear using LE boundary are presented. The paper is concluded in Section 5. 2. Numerical method 2.1. Fluid flow simulation by the lattice Boltzmann method In the current work, the two dimensional lattice Boltzmann method is used for fluid flow simulation. The discretized form of the Boltzmann equation with single relaxation time collision oper- ator is expressed as follows [2]: faðx þ eadt; t þ dtÞ ¼ faðx; tÞ � 1s faðx; tÞ � f eq a ðx; tÞ � � ; ð1Þ in which, fa; x; ea; s and dt are the particle distribu tion function, coordinat es of the grid points, discrete velocity vectors, dimension- less relaxatio n time and time step, respective ly. The relaxatio n time E. Jahanshahi Javaran et al. / Advance where /iðx; tÞ 2 ½0;1� is the density profile of the ith particle and Np is the total number of the particles in the computa tional domain. Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lee method s to introd uce shear into particle suspension s, Advanc ed Powde r Techn Several analytic al forms of the smoothed profile for the spherical particles exist [9]. In this work, the following function is used: /iðx; tÞ ¼ sðRi � jx � RiðtÞjÞ sðxÞ ¼ 0 x < �n=2 1 2 sin px ni þ 1 � � jxj 6 n=2 1 x > n=2 8>>< >>: ð6Þ in which, Ri;Ri and ni are the radius, position vector and interfacia l thickne ss of the ith particle, respectively. Based on this concentratio n field, the particle velocity field is constructed from the rigid motions of the Np particles as /ðx; tÞuPðx; tÞ ¼ XNP i¼1 /iðx; tÞ½ViðtÞ þxi � fx � RiðtÞg�; ð7Þ here, fVi;xig ði ¼ 1; . . . NPÞ are the translat ional and angula r veloci- ties of the ith particle, respec tively. Fluid nodes covered by the solid particles must have the same velocity as the solid particles. To this end, a body force is intro- duced inside the solid domain to enforce the fictitious fluid inside Fig. 1. Smoothed profile of a typical particle (solid line), dotted line is shown for comparison. wder Technology xxx (2013) xxx–xxx 3 As the particle nodes coincide with the fluid nodes, the fluid–so- lid interaction force is active and this force is zero outside the par- ticles’ domain. s–Edwards boundary conditions with smo othed profile-lattice Boltzman n olo gy (2013), http://dx.doi.org/10.1 016/j.a pt.2013.03.018 The hydrodyn amic force FHi and torque T H i exerted on each solid particle are obtained from the momentum conservation law as [9,10]: FHi ¼ Z 8pi q/nðuðx; tnÞ � upðx; tnÞd8pi ð11Þ THi ¼ Z 8pi ðx � Rni Þ � q/nðuðx; tnÞ � upðx; tnÞd8pi ð12Þ The translationa l and angular velocities of each particle are up- dated using: Vnþ1i ¼ Vni þ M�1pi Z tnþdt tn ðFHi þ Fci þ Fext i Þds ð13Þ xnþ1i ¼ xni þ I�1pi Z tnþdt ðTHi þ Text i Þds ð14Þ fects including the formation of a depletion zone near the wall and algorithm, the periodic boundary condition is applied on the flow (x-) and vorticity (z-) directions ; however, in the shear direction (y-), periodic domains are shifted continuously in time with a velocity equal to cLy, where c is the imposed shear rate and Ly is the domain length in the shear direction. Materials crossing the top shear border, both fluid and suspended solids, must undergo a shift in the flow direction position equal to �cLyt and a shift in velocity in the flow direction of �cLy before reappearing in the bot- tom of the simulation domain. Lees–Edwards boundary conditions were extended by Wagner and Pagonabar raga [14] to be used in LBM for fluid-only systems. In the proposed method, periodic boundary conditions are applied to the distribution functions leaving the domain in the directions perpendi cular to the velocity gradient (i.e. x- and z-directions). In addition, it is assumed that an image of the system is moving with velocity þ0:5ULE with respect to the system on the upper part and velocity �0:5ULE relative to the system on the lower part of the 4 E. Jahanshahi Javaran et al. / Advanced Powder Technology xxx (2013) xxx–xxx the resultant wall slip [12]. One way to deal with these issues is to increase the overall size of the simulatio n region. This causes numer- ical inefficiencies, especially when computational power is limited. Lees and Edwards [13] proposed a remarkable algorithm in 1972 for use in molecular dynamics, which consists of a simple modification of periodic boundary condition s that enables the sim- ulation of Couette flow without introduction of physical bound- aries. This type of introducing shear into the domain overcomes the problems and limitations of the previous method. In this tn In Eq. (13) Mpi is the particle mass, F c i and F ext i are the collision force and the external force on the ith particle, respectively. Text i in Eq. (14) is the external torque on the ith particle and Ipi is the particle moment of inertia tensor. Finally, the new particle posi- tions are calculated by: Rnþ1i ¼ Rni þ Z tnþdt tn Vids ð15Þ 3. Lees–Edwards (LE) boundary conditio n 3.1. LE boundary condition for fluid-only systems Planar Couette flow is frequently used to introduce a shear into a domain. A simulation where bulk behavior is desired is then largely contaminat ed with boundary effects. In other words, because of the presence of the moving walls, it is unavoidable to prevent wall ef- Fig. 2. Applying Lees–Edwards boundary conditions, periodic copies of the system Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lee method s to introd uce shear into particle suspension s, Advanced Powder Techn system (Fig. 2). Consequently, periodic domains are shifted contin- uously in time with a velocity equals �0:5ULE, and therefore a uniform shear rate c ¼ ULE=Ly is maintained. Implementa tion of LE boundary condition in the lattice Boltz- mann method is not as straightforw ard as molecular dynamics . Two steps in addition to the molecular dynamics are needed for its application in LBM. At first, since in LBM we are dealing with a distribution function, the way to implement the velocity jump when crossing the system boundary in the velocity gradient direc- tion will be to perform a Galilean transformation , i.e., if the distri- bution function in the system creates a momentum of qu, its image node will have a momentum of qðu � ULEÞ for upper-le aving parti- cles, and qðu þ ULEÞ for lower-leavi ng particles. In this way, we will generate a linear velocity profile with shear rate of c ¼ ULE=Ly. In order to get the momentum transfer, we need to calculate the dif- ference between a density fi at the system velocity and the one after the Galilean transformat ion has been applied. So, we have: Dfi ¼ fiðu � ULEÞ � fiðuÞ � f eqi ðu � ULEÞ � f eqi ðuÞ ð16Þ where Dfi is the change in the density moving in the i-direction due to Galilean transform ation. Note that the transformat ion will only be applied to the particle distribu tion functions crossing the LE boundar y. The above transformat ion can be rewritten as: f 0i ðu � ULEÞ ¼ f out i ðuÞ þ f eqi ðu � ULEÞ � f eqi ðuÞ ð17Þ where f 0i and f out i are post-collisio n distribution functions after the Galilean transfor mation and in the system, respective ly. are moving with respect to the system with a horizontal velocity of �0:5ULE . s–Edwards boundar y conditions with smoo thed profile-lattice Boltzman n ology (2013), htt p://dx.doi.org /10.1016/ j.apt.2013.0 3.018 As in the molecular dynamics, materials leaving the sheared do- main in the positive and negative y-direction must undergo a shift in the flow direction position equal to �ULEt before reappearing from the opposite direction. Since the displacemen t ULEt is a real number in general, so this displacemen t maps the particle distribution functions crossing top and bottom boundaries into off lattice positions. In the second step, we will have to interpolate the densities crossing a LE boundary from the off-lattic e positions to which they stream, onto the regular lattice. To this end, this displacemen t is split into an integer part ðDLEÞI and areal part ðDLEÞR, DLE ¼ ULEt ¼ ðDLEÞI þ ðDLEÞR ð18Þ Linear interpolation is used to transform the densities from the off-lattice positions to the lattice points. In more detail, particle distribution functions leaving the upper boundary and inserted back to the lower boundary are f2; f5; f6, for which, after linear inter- polation, we have: out 0 I shown in Fig. 4. LE boundaries in fluid-only systems is applied here to these particles . In more detail, the implementation of Lees–Edwards boundary condition for a particle crossing the top boundary is as follows: As soon as each part of the particle leaves the upper boundary, this part is entered from the lower boundary, but with two modi- fications. The first modification is the change of the x-component velocity ðvxÞ of that part of particle by �ULE and the second modi- fication is the change of the x-compon ent position ðxÞ of that part by �ULEt as v 0x ¼ vx � ULE v 0y ¼ vy x0 ¼ modððx � ULEtÞ; LxÞ y0 ¼ y � Ly ð21Þ where x0 and y0 are the shifted position in the flow and velocity gra- dient directions of that part of particle that leaves the upper bound- ary and reappea rs from the bottom boundary. In general, when the particle is close to a corner of the domain, this particle may be split in four parts (Fig. 3). In this figure, Xc;Yc and X 0 c; Y 0 c are the coordinates of the center of mass of the considered particle and the shifted particle, respectivel y. In this case, integrals in Eqs. (11) and (12) must be split in four parts and the hydrodynam ics force and torque on each part is sepa- rately calculated and the sum of them is exerted on the particle that is placed on the considered corner. Furthermore, the force acting on the fluid nodes located inside this solid particle (Eq. (9)) must be changed in fluid nodes covered by each of these four E. Jahanshahi Javaran et al. / Advanced Powder Technology xxx (2013) xxx–xxx 5 For lower leaving particle distribut ion functions, namely f4; f7; f8, linear interpolation results in: f out i ðx; tÞ ¼ S1f 0i ðx � modððDLEÞI; LxÞ � 1; tÞ þ S2f 0i ðx � modððDLEÞI; LxÞ; tÞ i ¼ 4;7;8 ð20Þ 3.2. Coupling Lees–Edwards boundary condition to the combined SPM–LBM In a suspensi on system sheared by employing Lees–Edwards boundary condition, as time goes on, it is possible for sus- pended particles to cross the upper and lower LE boundaries. Similar method to particle distribution functions crossing the fi ðx; tÞ ¼ S1fi ðx þ modððDLEÞ ; LxÞ þ 1; tÞ þ S2f 0i ðx þ modððDLEÞI; LxÞ; tÞ i ¼ 2;5;6 ð19Þ in which Lx is the length of the domain in the x-direction, S1 ¼ modðDLE;1Þ, S2 ¼ 1 � modðDLE;1Þ and the modulus, or remain- der, operator (mod) divides two numbers in its argumen t and returns only the remainder as result. The direction of above and fol- lowing distribution functions with respect to the geometry is Fig. 3. Schematic representation of a particle that has crossed the top boundary. In this case the particle is placed close to the top right corner and is split in four parts. Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lee method s to introd uce shear into particle suspension s, Advanc ed Powde r Techn parts. Similar procedure is applied to the particles that cross the lower boundary . In this case, leaving part of the particle reappears in the simulatio n domain from the top boundary, but with change in the x-component velocity by þULE and x-component position by þULEt as v 0x ¼ vx þ ULE v 0y ¼ vy x0 ¼ modððx þ ULEtÞ; LxÞ y0 ¼ y þ Ly ð22Þ Fig. 4. Computational grid with particle distribution functions leaving the upper boundary, being reinserted from the lower boundary after applying Galilean transformation. s–Edwards boundary conditions with smo othed profile-lattice Boltzman n olo gy (2013), http://dx.doi.org/10.1 016/j.a pt.2013.03.018 by the LE boundary condition versus simulation time for different suspension methods. ary for the first time. (b) The time evolution of the particle angular velocity. Powder Technology xxx (2013) xxx–xxx Fig. 5. (a) The vertical velocity (Vy) of a particle placed in a domain which is sheared At the time step of around 12,800, the particle center of mass crosses the LE bound 6 E. Jahanshahi Javaran et al. / Advanced For the particles that cross the domain in the periodic direction x, each part of the particle that leaves the periodic boundaries from the positive and negative x directions are reinserted in the opposite side without any modifications. In all previous investigatio ns found in the literature, linear inter- polation was used for LE boundary -crossing particle distribution functions. This makes some bumps in the angular velocity of the particle when crossing the LE boundary. Such behavior is observed in the present study using linear interpolation which will be presente d in the results part of the paper. Since, we are looking for a method to conserve the law of the translationa l/rotational momentum; it is desirable to eliminate these bumps. Such behavior may be due to the linear interpolati on method used which intro- duces some errors (artificial diffusion) into the computati ons. In the present study, higher-order interpolations were used to check if this shortcomin g resolves. It is found that using third-order inter- polation makes these bumps to disappea r. The computati onal grid with four typical particle distribution functions for use in the third-order interpolation is shown in Fig. 4. For particle distribution functions crossing the top boundary and reappeari ng from the bottom boundary after third-order inter- polation we have: Fig. 6. Effect of the relaxation time on (a) particle vertical velocity, (b) particle angular velocity, using first-order interpolation for particle distribution functions leaving the domain in the shear direction. Fig. 7. Effect of relaxation time on the angular velocity of the particle using third- order interpolation for particle distribution functions leaving the domain in the shear direction. Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lees–Edwards boundar y conditions with smoo thed profile-lattice Boltzman n method s to introd uce shear into particle suspension s, Advanced Powder Technology (2013), htt p://dx.doi.org /10.1016/ j.apt.2013.0 3.018 f the d Powder Technology xxx (2013) xxx–xxx 7 Fig. 8. Effect of superposed velocity on (a) the angular velocity o E. Jahanshahi Javaran et al. / Advance f out i ðx; tÞ ¼ S2ðS2 þ 1ÞðS2 þ 2Þ 6 f 0i ðx þ modððDLEÞI; LxÞ; tÞ þ S1ðS2 þ 1ÞðS2 þ 2Þ 2 f 0i ðx þ modððDLEÞI; LxÞ þ 1; tÞ þ S1S2ðS2 þ 2Þ�2 f 0 i ðx þ modððDLEÞI; LxÞ þ 2; tÞ þ S1S2ðS2 þ 1Þ 6 f 0i ðx þ modððDLEÞI; LxÞ þ 3; tÞ i ¼ 2;5;6 ð23Þ For bottom-leaving particles, third-order interpolation results in: f out i ðx; tÞ ¼ S1S2ðS2 þ 1Þ 6 f 0i ðx � modððDLEÞI; LxÞ � 3; tÞ þ S1S2ðS2 þ 2Þ�2 f 0 i ðx � modððDLEÞI; LxÞ � 2; tÞ þ S1ðS2 þ 1ÞðS2 þ 2Þ 2 f 0i ðx � modððDLEÞI; LxÞ � 1; tÞ þ S2ðS2 þ 1ÞðS2 þ 2Þ 6 f 0i ðx � modððDLEÞI; LxÞ; tÞ i ¼ 4;7;8 ð24Þ Fig. 9. Effect of shear rate on (a) the angular velocity of the particle, (b) particle vertical v domain in the shear direction. Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lee method s to introd uce shear into particle suspension s, Advanc ed Powde r Techn particle, (b) particle vertical velocity ðc ¼ 2:500E � 5; s ¼ 1:00Þ. 4. Validation tests and numerical results 4.1. One particle going across the Lees–Edwards boundary For the simulations presente d in this section, a computati onal domain consisting of Lx � Ly ¼ 128 � 128 grid points was used. A particle with radius of R = 4.8 was placed in the domain which is elocity, using third-order interpolation for particle distribution functions leaving the Fig. 10. Simulation setup of two interacting circular particles in a sheared domain using LE boundary condition. s–Edwards boundary conditions with smo othed profile-lattice Boltzman n olo gy (2013), http://dx.doi.org/10.1 016/j.a pt.2013.03.018 velocity of this particle will be half of the imposed shear rate. With applicati on of a Galilean transformation the sheared flow is superimpos ed with a constant velocity (u(s,x), u(s,y)). Because the particle is placed at the center of the domain, its velocity after superimpos ition would be (u(s,x), u(s,y)). To satisfy the conserva- tive property of the translational and rotational momentum, the particle vertical and angular velocities should be kept at u(s,y) and half of the shear rate, respectively . This happens until the par- ticle crosses the LE boundary. At this time, some bumps are seen in the angular velocity. Although, the magnitudes of the bumps are very small, this bumps violate conservative property of the rotational momentum . In applying Galilean transformat ion to boundary crossing particle distribution functions, as was stated in the work by Wagner and Pagonabarra ga [14], there is an error term (Galilean transformation of the non-equilib rium part of the particle distribut ion function). This term is small and enters only into the pressure moments . The relaxation time in this term has a decreasing effect, namely by increasing the relaxation time this term decrease s. So, by increasing the relaxation time we expect that the magnitude of the bumps decreases. These errors are hid- two parallel walls moving in opposite directions (lines with symbols). Fig. 12. Effect of interfacial thickness on the particle vertical velocity. Po sheared using LE boundary condition. This particle will finally cross the LE boundary as explained below. Density ratio was set to qs=qf ¼ 10. 4.1.1. Validation test In order to validate the computati onal frame of LBM–LE meth- od, a horizontal shear rate of c ¼ 1:25 � 10�5 is introduced and maintained using Lees–Edwards boundary condition; the system is allowed to reach its steady state. Then a specific velocity of us ¼ ðus;x;us;yÞ ¼ ð0:025;0:005Þ is superimposed on the established steady-state flow velocity. At this time, a particle is placed at the mid-point of the simulatio n domain. The velocity of the center of mass of the particle would be equal to the superimposed velocity, vpðt ¼ 0Þ ¼ us as the steady state velocity at the center of the sheared velocity field is zero. The initial angular velocity of the particle is set to xðt ¼ 0Þ ¼ 0. After a short equilibration time, the near-particl e flow field is equilibrated and the particle moves freely along with the fluid and crosses the LE boundary. The y-com- ponent velocity of the center of mass of the particle is shown in Fig. 5a. The results with other suspension methods are also shown for comparison. Ideally, if the suspension method conserves the law of the translationa l momentum , the velocity should stay constant; a behavior which is clearly observed for the present com- putations using SPM–LE. However, deviation from a constant value is observed from the beginning for ALD method [19], while these deviations add up as time goes on. Ladd’s method conserves the constant speed perfectly until the particle crosses LE boundary. After the particle crossed LE boundary, its vertical velocity ap- proaches its initial value again. The reason for such behavior roots in non-Galilea n Invariance effects of the momentum exchange algorithm which is used in both ALD and Ladd’s method. Also shown in this figure are the results of the proposed corrected momentum algorithm for suspensions (CMES) by the authors of Ref. [16] which shows its preservation of Galilean Invariance. How- ever, as time goes on, the error due to violation of Galilean invari- ance grows. In addition, ideally, the angular velocity of the particle should be equal to one half of the shear rate in magnitude after the particle’s motion reaches equilibrati on. Fig. 5b depicts the time evolution of the angular velocity of the particle. This figure shows that the particle angular velocity reaches the expected value of c=2 after a short time; however there are bumps in the angular velocity when particle crosses the LE boundary. 4.1.2. Results In the first part of the computations , a particle is placed at the center of a shear flow with its translationa l velocity equal to the fluid velocity at that point and angular velocity equal to zero. It should be noted that the whole system is superpos ed with a veloc- ity us ¼ ðus;x;us;yÞ. Here, the effects of relaxation time, superpos ed velocity, shear rate and third-order interpolation on the angular velocity and particle vertical velocity are presented using SPM–LE method. Since relaxation time of s ¼ 1 is a special case and is believed to hide numerica l errors, at first the effect of the relaxation time is investigated and results of the vertical and angular velocities of the particle for different relaxation times are plotted in Fig. 6a and b respectivel y. Fig. 6a shows no tangible effect of s on the Vy while Fig. 6b depicts some bumps appeared for angular velocity of the particle crossing LE boundary. These bumps are reduced when larger relaxation times are used in the simulation. The max- imum height of the bump at relaxation time of s = 0.75 is about 0.04. This value decreases to 0.02 at s = 1.0 and 0.01 for relaxation times of 1.5 and 2.0. 8 E. Jahanshahi Javaran et al. / Advanced The reason for such behavior can be explained as follows. When a particle placed at the center of the sheared domain, reaches its equilibrium position with the surrounding fluid, the angular Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lee method s to introd uce shear into particle suspension s, Advanced Powder Techn Fig. 11. Flow trajectories of two interacting circular cylinders with radius R in a shear flow obtained using LE boundary condition (solid lines) and by employing wder Technology xxx (2013) xxx–xxx den when computing the force exerted on the particle because of the symmetry, while this errors present when computing the tor- que on the particle and angular velocity of the particle. In general, s–Edwards boundar y conditions with smoo thed profile-lattice Boltzman n ology (2013), htt p://dx.doi.org /10.1016/ j.apt.2013.0 3.018 d Po E. Jahanshahi Javaran et al. / Advance the errors in the angular velocity are truncation and numerical errors which can be decreased by increasing the relaxation time or by using higher-order interpolati ons. It should be mentioned that these bumps can also be seen in the results of the recently published work in this area in which first-order interpolation was used for the LE boundary crossing distribution functions [17]. In order to search for the cause of these errors, simulatio ns are done using third-order interpolation for particle distribution func- tions leaving the domain in the shear(y-) direction and it is found that these bumps are eliminated when third-order interpolation is used. In Fig. 7, time evolution of the angular velocity using third- order interpolati on is plotted in which the effect of higher-order interpolation on x is obvious as compared to Fig. 6b which shows it counterpart for the first-order interpolation. Therefore, here after, third-order interpolation is used in all simulations . In another simulation, the effect of superpos ed velocity us on the angular velocity of the particle for constant shear rate c ¼ 2:5 � 10�5 are explored and the results are shown in Fig. 8a. It is seen that after a short time, x reaches the expected value of c/2 with no effect of the superposed velocity on the angular veloc- ity. Fig. 8b shows the time evolution of particle vertical velocity for different superposed velocities. Although the flow field is initial- ized with different superposed velocities, no changes are found in the y-component of the particle velocity. It should be noted that, the y-component of the superposed velocity is kept constant in all simulations which causes the center of mass of the particle to cross the top LE boundary after the same times in all simulations ; its first-time-crossing occurs at the time step of t = 12800. Fig. 13. Effect of interfacial thickness on the angular vel Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lee method s to introd uce shear into particle suspension s, Advanc ed Powde r Techn wder Technology xxx (2013) xxx–xxx 9 Furthermore, the effect of shear rate c on the conservation law of the translational and rotational momentum of the method is explored . In this test, the superimposed velocity stays constant at us ¼ ð0:0125;0:005Þ. The results for the angular velocity and parti- cle vertical velocity are shown in Fig. 9a and b which shows similar behavior as those observed in previous simulations.One of the paramete rs which its variation may affect the computati onal results in SPM is interfacial thickness, n. In order to investigate depende ncy of results on the interfacia l thickness, particle vertical and angular velocities were computed for three interfacial thick- nesses; their results are shown in Figs. 12 and 13. As it is observed from these figures, using SPM with different interfacial thicknesses has no effect on the particle vertical velocity, while reduces small fluctuations in angular velocity of the particle. In addition, varia- tion of interfacial thickness has no effect on the equilibrium value of the angular velocity. It should be noted that these figures are based on third-order interpolati on of particle distribut ion functions crossing the LE boundary. The values of shear rate and relaxation time used in these computati on are selected as: c ¼ 0:625E � 5; s ¼ 0:75. 4.2. Two and many particles in a sheared domain using Lees–Edwards boundary condition The interaction between two circular cylinders subjected to simple shear flow obtained using Lees–Edwards boundary condi- tion is presented in this section to show the accuracy of the ocity of the particle (a) n = 0.0 (b) n = 1.0 (c) n = 2.0. s–Edwards boundary conditions with smo othed profile-lattice Boltzman n olo gy (2013), http://dx.doi.org/10.1 016/j.a pt.2013.03.018 10 E. Jahanshahi Javaran et al. / Advanced Powder Technology xxx (2013) xxx–xxx coupling LE boundary to combined SPM–LBM at near distances between suspended particles . The simulation set up is shown in Fig. 10. Two circular cylinders are placed in the domain at heights symmetr ical around the hori- zontal centerlin e; their centerlin es were placed at a distance of 19.25 R from the LE boundary and 0.75 R from the flow centerline. The initial horizontal distance between the two cylinders is 10 R where R is the radius of each particle. The particle based shear Rey- nolds number is Rep = 4cR2/t in which t is the kinematic viscosity of the fluid. In addition, when the gap width between two particles is in the order of one lattice spacing, the lubrication force is how- ever not exactly resolved with the lattice Boltzmann method. In the present work, the repulsive force proposed by based on the Fig. 14. A snapshot of the concentrated suspensions of 117 particles after the suspension reached a statistically steady-state. work of Feng and Michaeli des [8] is used to resolve this problem. Trajectories of these two interacting particles are shown in Fig. 11 along with those obtained by Kromkamp et al. [20] in which introducing shear into the domain was done by two parallel walls moving in opposite directions. They used a grid point distribution of 640 � 320 to eliminate wall effects. As is observed from this fig- ure, the present results are in good agreement with those presented in Ref. [20]. In order to show the applicability of the combined SPM–LEBC to concentrated suspensi ons, 117 particles are placed initially in a random manner in the domain while the system undergoes shear. A snapshot of these particles after the suspension reached a statis- tically steady-state is shown in Fig. 14. As it is observed in this fig- ure about 6 particles cross the top and bottom Lees–Edwards boundaries. The volume fraction of solid in this case is 0.25 and particle-bas ed Reynolds number is 0.31. 5. Conclusion s In order to solve a system of a large number of particles at a rea- sonable computational cost and remove wall effects, Lees–Edwards boundary condition s are used to introduce shear into the simula- tion region. It is important to set up a method which reduces the effect of numerical non-Galilea n invariance or violation of conser- vation law of the translational and rotational momentum . To this Please cite this article in press as: E. Jahanshahi Javaran et al., Combinin g Lee method s to introd uce shear into particle suspension s, Advanced Powder Techn end, in the present work, smoothed profile method (SPM), which is a node-based method, was combined with Lees–Edwards bound- ary conditions to simulate a suspension flow. A particle was al- lowed to cross the Lees–Edwards boundary and its vertical and angular velocities were recorded during the simulation. It was found that the vertical velocity stays constant which means the combined SPM–LE strictly conserves the Galilean invarianc e. Also, the angular velocity of the particle reaches the expected value of the half of the imposed shear rate, after a short time. Furthermore, by using first-order interpolation for the distribution functions, there were some bumps in the angular velocity of the particle when it crosses the LE boundary. These bumps were eliminated when third-order interpolation was used. Moreover, the flow tra- jectories of two interactin g particles placed in a sheared domain using LE boundary condition were shown and compare d with the ones found in the literature for which good agreement was ob- served. The last simulatio n was done for 117 particles placed ran- domly in a shear flow. A snapshot of these particles after the suspensi on reached a statistically steady-state showed about 6 particles cross the top and bottom Lees–Edwards boundaries. References [1] J.F.A. Brady, G. Bossis, Stokesian dynamics, Ann. Rev. Fluid Mech. 20 (1988) 111–157. [2] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond , Oxford University Press, Oxford, 2001 . [3] A.J.C. Ladd, R. Verberg, Lattice-Boltzmann simulation of particle–fluid suspensions, J. Stat. Phys. 104 (2001) 1191–1251. [4] C.K. Aidunand, J.R. Clausen, Lattice-Boltzmann method for complex flows, Ann. Rev. Fluid Mech. 42 (2010) 439–472. [5] A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation Part I. Theoretical foundation, J. Fluid Mech. 271 (1994) 285–310. [6] A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part II. Numerical results, J. Fluid Mech. 271 (1994) 311– 339. [7] M. Bouzidi, M. Firdaoussand, P. Lallemand, Momentum transfer of a Boltzmann-lattice fluid with boundaries, Phys. Fluids 13–11 (2001) 3452– 3459. [8] Z. Fengand, E.E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems, J. Comput. Phys. 195 (2004) 602–628. [9] Y. Nakayama, R. Yamamoto, Simulation method to resolve hydrodynamic interactions in colloidal dispersions, Phys. Rev. E 71 (2005) 036707 . [10] S. Jafari, R. Yamamoto, M. Rahnama, Lattice-Boltzmann method combined with smoothed-profile method for particulate suspensions, Phys. Rev. E 83 (2011) 026702 . [11] J.J. Stickle, R.L. Powell, Fluid mechanics and rheology of dense suspensions, Annu. Rev. Fluid Mech. 37 (2005) 129–149. [12] J. Kromkamp, D. van denEnde, D. Kandhai, R. van der Smanand, R. Boom, Lattice Boltzmann simulation of 2D and 3D non-Brownian suspensions in Couette flow, Chem. Eng. Sci. 61 (2006) 858–873. [13] A.w. Lees, S.F. Edwards, The computer study of transport processes under extreme conditions, J. Phys. C: Solid State Phys. 5 (1972) 1921–1929. [14] A.J. Wagner, I. Pagonabarraga, Lees–Edwards boundary conditions for Lattice Boltzmann, J. Stat. Phys. 107 (1/2) (2002). [15] E. Lorenz, A.G. Hoekstra, A. Caiazzo, Lees–Edwards boundary conditions for lattice Boltzmann suspension simulations, Phys. Rev. E 79 (2009) 036706 . [16] E. Lorenz, A. Caiazzo, A.G. Hoekstra, Corrected momentum exchange method for lattice Boltzmann simulations of suspension flow, Phys. Rev. E 79 (2009) 036705. [17] G. Zhou, L. Wang, W. Ge, Galilean-invariant algorithm coupling immersed moving boundary conditions and Lees–Edwards boundary conditions, Phys. Rev. E 84 (2011) 066701 . [18] X. Luo, M.R. Maxey, G.E. Karniadakis, Smoothed profile method for particulate flows: error analysis and simulations, J. Comput. Phys. 228 (2009) 1750–1769. [19] C.K. Aidun, Y. Lu, E. Ding, Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation, J. Fluid Mech. 373 (1998) 287 . [20] J. Kromkamp, D. van den Ende, D. Kandhai, R. van der Sman, R. Boom, Shear- induced self-diffusion and microstructure in non-Brownian suspensions at non-zero Reynolds numbers, J. Fluid Mech. 529 (2005) 253–278. [21] J.R. Clausen, C.K. Aidun, Galilean invariance in the lattice-Boltzmann method and its effect on the calculation of rheological properties in suspensions, Int. J. Multiphase Flow 35 (2009) 307–311. s–Edwards boundar y conditions with smoo thed profile-lattice Boltzman n ology (2013), htt p://dx.doi.org /10.1016/ j.apt.2013.0 3.018 Combining Lees–Edwards boundary conditions with smoothed profile-lattice Boltzmann methods to introduce shear into particle suspensions 1 Introduction 2 Numerical method 2.1 Fluid flow simulation by the lattice Boltzmann method 2.2 Simulation of fluid–solid interaction by the smoothed profile method (SPM) 3 Lees–Edwards (LE) boundary condition 3.1 LE boundary condition for fluid-only systems 3.2 Coupling Lees–Edwards boundary condition to the combined SPM–LBM 4 Validation tests and numerical results 4.1 One particle going across the Lees–Edwards boundary 4.1.1 Validation test 4.1.2 Results 4.2 Two and many particles in a sheared domain using Lees–Edwards boundary condition 5 Conclusions References