Skip to main content
Log in

Three-dimensional simulations of dilute and concentrated suspensions using smoothed particle hydrodynamics

  • Published:
Computational Particle Mechanics Aims and scope Submit manuscript

Abstract

A three-dimensional model for a suspension of rigid spherical particles in a Newtonian fluid is presented. The solvent is modeled with smoothed particle hydrodynamics method, which takes into account exactly the long-range multi-body hydrodynamic interactions between suspended spheres. Short-range lubrication forces which are necessary to simulate concentrated suspensions, are introduced pair-wisely based on the analytical solution of Stokes equations for approaching/departing objects. Given that lubrication is singular at vanishing solid particle separations, an implicit splitting integration scheme is used to obtain accurate results and at the same time to avoid prohibitively small simulation time steps. Hydrodynamic interactions between solid particles, at both long-range and short-range limits, are verified against theory in the case of two approaching spheres in a quiescent medium and under bulk shear flow, where good agreements are obtained. Finally, numerical results for the suspension viscosity of a many-particle system are shown and compared with analytical solutions available in the dilute and semi-dilute case as well as with previous numerical results obtained in the concentrated limit.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Institutional subscriptions

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5

Similar content being viewed by others

References

  1. Mewis J, Wagner NJ (2011) Colloidal suspension rheology. Cambridge University Press, Cambridge Books Online

  2. Maitland G (2000) Oil and gas production. Curr Opin Colloid Interface Sci 5(5):301–311

    Article  Google Scholar 

  3. Bian X, Ellero M (2014) A splitting integration scheme for the SPH simulation of concentrated particle suspensions. Comput Phys Commun 185(1):53–62

    Article  MathSciNet  Google Scholar 

  4. Brady JF, Bossis G (1988) Stokesian dynamics. Annu Rev Fluid Mech 20:111–157

    Article  Google Scholar 

  5. Sierou A, Brady JF (2001) Accelerated Stokesian dynamics simulations. J Fluid Mech 448:115–146

    MATH  Google Scholar 

  6. Kutteh R (2010) Rigid body dynamics approach to Stokesian dynamics simulations of nonspherical particles. J Chem Phys 132(17):174107

    Article  Google Scholar 

  7. Swan JW, Brady JF (2007) Simulation of hydrodynamically interacting particles near a no-slip boundary. Phys Fluids 19(11):113306

    Article  MATH  Google Scholar 

  8. Schaink HM, Slot JJM, Jongschaap RJJ, Mellema J (2000) The rheology of systems containing rigid spheres suspended in both viscous and viscoelastic media, studied by Stokesian dynamics simulations. J Rheol 44(3):473–498

    Article  Google Scholar 

  9. Ladd A, Verberg R (2001) Lattice-Boltzmann simulations of particle-fluid suspensions. J Stat Phys 104(5–6):1191–1251

    Article  MathSciNet  MATH  Google Scholar 

  10. Owen DRJ, Leonardi CR, Feng YT (2011) An efficient framework for fluid-structure interaction using the lattice Boltzmann method and immersed moving boundaries. Int J Numer Methods in Eng 87(1–5):66–95

    Article  MathSciNet  MATH  Google Scholar 

  11. Koch DL, Ladd AJ (1997) Moderate Reynolds number flows through periodic and random arrays of aligned cylinders. J Fluid Mech 349:31–66

    Article  MathSciNet  MATH  Google Scholar 

  12. Hoef MVD, Beetstra R, Kuipers J (2005) Lattice-Boltzmann simulations of low-Reynolds-number flow past mono-and bidisperse arrays of spheres: results for the permeability and drag force. J Fluid Mech 528:233–254

    Article  MathSciNet  MATH  Google Scholar 

  13. Pan C, Luo LS, Miller CT (2006) An evaluation of lattice Boltzmann schemes for porous medium flow simulation. Comput Fluids 35(8):898–909

    Article  MATH  Google Scholar 

  14. Uhlmann M (2005) An immersed boundary method with direct forcing for the simulation of particulate flows. J Comput Phys 209(2):448–476

    Article  MathSciNet  MATH  Google Scholar 

  15. Usabiaga FB, Pagonabarraga I, Delgado-Buscalioni R (2013) Inertial coupling for point particle fluctuating hydrodynamics. J Comput Phys 235:701–722

    Article  MathSciNet  Google Scholar 

  16. Vázquez-Quesada A, Usabiaga FB, Delgado-Buscalioni R (2014) A multiblob approach to colloidal hydrodynamics with inherent lubrication. J Chem Phys 141(20):204102

    Article  Google Scholar 

  17. Pan W, Caswell B, Karniadakis GE (2009) Rheology, microstructure and migration in Brownian colloidal suspensions. Langmuir 26(1):133–142

    Article  Google Scholar 

  18. Boek E, Coveney P, Lekkerkerker H (1996) Computer simulation of rheological phenomena in dense colloidal suspensions with dissipative particle dynamics. J Phys: Condens Matter 8(47):9509

    Google Scholar 

  19. Boek E, Coveney P, Lekkerkerker H, van der Schoot P (1997) Simulating the rheology of dense colloidal suspensions using dissipative particle dynamics. Phys Rev E 55(3):3124

    Article  Google Scholar 

  20. Moshfegh A, Jabbarzadeh A (2015) Dissipative particle dynamics: effects of parameterization and thermostating schemes on rheology. Soft Mater 13(2):106–117

    Article  Google Scholar 

  21. Monaghan JJ (2005) Smoothed particle hydrodynamics. Rep Prog Phys 68(8):1703

    Article  MathSciNet  MATH  Google Scholar 

  22. Sbalzarini I, Walther J, Bergdorf M, Hieber S, Kotsalis E, Koumoutsakos P (2006) PPM-A highly efficient parallel particle-mesh library for the simulation of continuum systems. J Comput Phys 215(2):566–588

    Article  MATH  Google Scholar 

  23. Español P, Revenga M (2003) Smoothed dissipative particle dynamics. Phys Rev E 67(2):026705

    Article  Google Scholar 

  24. Vázquez-Quesada A, Ellero M, Español P (2009) Consistent scaling of thermal fluctuations in smoothed dissipative particle dynamics. J Chem Phys 130(3):034901

    Article  Google Scholar 

  25. Grmela M, Öttinger HC (1997) Dynamics and thermodynamics of complex fluids. I. Development of a general formalism. Phys Rev E 56(6):6620

    Article  MathSciNet  Google Scholar 

  26. Öttinger HC, Grmela M (1997) Dynamics and thermodynamics of complex fluids. II. Illustrations of a general formalism. Phys Rev E 56(6):6633

    Article  MathSciNet  Google Scholar 

  27. Bian X, Litvinov S, Qian R, Ellero M, Adams NA (2012) Multiscale modeling of particle in suspension with smoothed dissipative particle dynamics. Phys Fluids 24(1):012002

    Article  Google Scholar 

  28. Kulkarni PM, Fu CC, Shell MS, Gary Leal L (2013) Multiscale modeling with smoothed dissipative particle dynamics. J Chem Phys 138(23):234105. http://scitation.aip.org/content/aip/journal/jcp/138/23/10.1063/1.4810754

  29. Bian X, Litvinov S, Ellero M, Wagner NJ (2014) Hydrodynamic shear thickening of particulate suspension under confinement. J Non-Newton Fluid Mech 213:39–49

    Article  Google Scholar 

  30. Morris JP, Fox PJ, Zhu Y (1997) Modeling low Reynolds number incompressible flows using SPH. J Comput Phys 136(1):214–226

    Article  MATH  Google Scholar 

  31. Hu X, Adams N (2006) Angular-momentum conservative smoothed particle dynamics for incompressible viscous flows. Phys Fluids 18(10):101702

    Article  Google Scholar 

  32. Müller K, Fedosov DA, Gompper G (2015) Smoothed dissipative particle dynamics with angular momentum conservation. J Comput Phys 281:301–315

    Article  MathSciNet  Google Scholar 

  33. Götze IO, Noguchi H, Gompper G (2007) Relevance of angular momentum conservation in mesoscale hydrodynamics simulations. Phys Rev E 76(4):046705

    Article  Google Scholar 

  34. Monaghan JJ (1994) Simulating free surface flows with SPH. J Comput Phys 110(2):399–406

    Article  MATH  Google Scholar 

  35. Kim S, Karrila SJ (1991) Microhydrodynamics: principles and selected applications. Butterworth-Henemann, Boston

    Google Scholar 

  36. Nguyen NQ, Ladd A (2002) Lubrication corrections for lattice-Boltzmann simulations of particle suspensions. Phys Rev E 66(4):046708

    Article  Google Scholar 

  37. Melrose J, Ball R (1995) The pathological behaviour of sheared hard spheres with hydrodynamic interactions. EPL (Europhysics Letters) 32(6):535

    Article  Google Scholar 

  38. Ball R, Melrose J (1995) Lubrication breakdown in hydrodynamic simulations of concentrated colloids. Adv Colloid Interface Sci 59:19–30

    Article  Google Scholar 

  39. Dratler D, Schowalter W (1996) Dynamic simulation of suspensions of non-Brownian hard spheres. J Fluid Mech 325:53–77

    Article  MATH  Google Scholar 

  40. Brady JF, Morris JF (1997) Microstructure of strongly sheared suspensions and its impact on rheology and diffusion. J Fluid Mech 348:103–139

    Article  MATH  Google Scholar 

  41. Brady JF, Bossis G (1985) The rheology of concentrated suspensions of spheres in simple shear flow by numerical simulation. J Fluid Mech 155:105–129

    Article  Google Scholar 

  42. Batchelor G, Green JT (1972) The hydrodynamic interaction of two small freely-moving spheres in a linear flow field. J Fluid Mech 56(02):375–400

    Article  MATH  Google Scholar 

  43. Jeffrey D, Onishi Y (1984) Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow. J Fluid Mech 139:261–290

    Article  MATH  Google Scholar 

  44. Einstein A (1906) Eine neue Bestimmung der Moleküldimensionen. Annalen der Physik 324(2):289–306

    Article  MATH  Google Scholar 

  45. Batchelor G, Green JT (1972) The determination of the bulk stress in a suspension of spherical particles to order \(\text{ c }^2\). J Fluid Mech 56(03):401–427

    Article  MATH  Google Scholar 

  46. Sierou A, Brady J (2002) Rheology and microstructure in concentrated noncolloidal suspensions. J Rheol 46(5):1031–1056

    Article  Google Scholar 

  47. Bertevas E, Fan X, Tanner RI (2010) Simulation of the rheological properties of suspensions of oblate spheroidal particles in a Newtonian fluid. Rheologica Acta 49(1):53–73

    Article  Google Scholar 

  48. Vázquez-Quesada A, Ellero M, Español P (2012) A SPH-based particle model for computational microrheology. Microfluid Nanofluid 13(2):249–260

    Article  Google Scholar 

  49. Chen S, Phan-Thien N, Khoo BC, Fan XJ (2006) Flow around spheres by dissipative particle dynamics. Phys Fluids 18(10):103605

    Article  MathSciNet  MATH  Google Scholar 

Download references

Acknowledgments

LRZ’s support for providing computing time on HPC system SuperMUC within the project “Multiscale Modelling of Particles in Suspension” (ID: pr58ye) is gratefully acknowledged. The authors also appreciate the useful discussions with Sergey Litvinov and his help calculating the Batchelor trajectories.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Adolfo Vázquez-Quesada.

Appendices

Appendix 1: Rigid body dynamics

Some pioneers of particle methods used

$$\begin{aligned} \hat{{\mathbf{v}}}_b = {\mathbf{V}}_{\alpha } + ({\mathbf{r}}_b - {\mathbf{R}}_{\alpha }) \times \varvec{\Omega }_{\alpha } \end{aligned}$$
(12)

to update the position of each solid-boundary-particle b at every time step. However, when the tangential component \(({\mathbf{r}}_b - {\mathbf{R}}_{\alpha }) \times \varvec{\Omega }_{\alpha }\) is multiplied by a finite time step \({\Delta } t\), the trajectory of each particle b can cause an artificially expanding solid structure, which causes numerical instability. Therefore, we adopt rotation matrix [49] to track the rotation of a solid particle in addition to its translation.

1.1 Rotation matrix

If \(\varvec{\Omega }_{\alpha }=\left[ \omega _x \, \omega _y \, \omega _z\right] \) is the angular velocity of the solid particle \(\alpha \), a rotation vector can be defined as

$$\begin{aligned} \varvec{\Phi } = {\Delta } t \varvec{\Omega }_{\alpha }. \end{aligned}$$
(13)

Suppose that \(\hat{{\mathbf{I}}}=[l_x\, l_y\,l_z]\) is the unit vector of \(\varvec{\Phi }\), namely, \(\hat{{\mathbf{I}}}=\frac{\varvec{\Phi }}{|\varvec{\Phi }|}=\frac{\varvec{\Phi }}{\phi }\), the resultant rotation matrix is defined as

$$\begin{aligned} \mathbf{RM}(\varvec{\Phi }) = \mathbf{RM}(\hat{\mathbf{I}}, \phi ) = \left[ \begin{array}{c@{\quad }c@{\quad }c} RM_{11}&{}RM_{12}&{}RM_{13}\\ RM_{21}&{}RM_{22}&{}RM_{23}\\ RM_{31}&{}RM_{32}&{}RM_{33} \end{array} \right] , \end{aligned}$$
(14)

where

$$\begin{aligned} RM_{11}= & {} l^2_x(1-\cos \phi ) + \cos \phi , \\ RM_{12}= & {} l_xl_y(1-\cos \phi ) - l_z\sin \phi , \\ RM_{13}= & {} l_xl_z(1-\cos \phi ) + l_y\sin \phi , \\ RM_{21}= & {} l_xl_y(1-\cos \phi ) + l_z\sin \phi , \\ RM_{22}= & {} l^2_y(1-\cos \phi ) + \cos \phi , \\ RM_{23}= & {} l_yl_z(1-\cos \phi ) - l_x\sin \phi , \\ RM_{31}= & {} l_xl_z(1-\cos \phi ) - l_y\sin \phi , \\ RM_{32}= & {} l_yl_z(1-\cos \phi ) + l_x\sin \phi , \\ RM_{33}= & {} l^2_z(1-\cos \phi ) + \cos \phi . \end{aligned}$$

We note that \(\mathbf{RM}^{-1}=\mathbf{RM}^t\).

The relative position of each boundary particle b to the solid particle center in fixed OXYZ coordinate at n time step will be updated as

$$\begin{aligned} \mathbf{rx}^{n}_b=\mathbf{RM}(\varvec{\Phi }^n)\mathbf{rx}^{n-1}_\mathbf{b}, \end{aligned}$$
(15)

where

$$\begin{aligned} \varvec{\Phi }^n = {\Delta } t \varvec{\Omega }^{n-1}_{\alpha }. \end{aligned}$$
(16)

The absolute position of each boundary particle b at n time step will be the sum of both rotation and translation:

$$\begin{aligned} {\mathbf{x}}^{n}_{b}= & {} \mathbf{rx}^{n}_{b}+(\mathbf{X}^{n-1}_{\alpha } + {{\Delta } \mathbf X}^{n}_{\alpha }), \nonumber \\= & {} \mathbf{rx}^{n}_b+\mathbf{X}^{n}_{\alpha }, \end{aligned}$$
(17)

where \({\mathbf{X}}_{\alpha }\) is the position of the solid particle center and \(\varvec{{\Delta } X}_{\alpha }\) is its translational displacement.

To record the history of rotation, we store the accumulative rotation matrix \(\mathbf{ARM}^{n}\) at n time step as

$$\begin{aligned} \mathbf{ARM}^{n} = \mathbf{RM}^{n} \mathbf{RM}^{n-1} \dots \,\mathbf{RM}^1 \mathbf{RM}^0, \end{aligned}$$
(18)

where \(\mathbf{RM}^{n}=\mathbf{RM}(\varvec{\Phi }^n)\). \(\mathbf{ARM}^{n}\) may be used to study the rotational behavior of the solid particle, such as the rational diffusion coefficient.

Since the moment of inertia tensor is real and symmetric, there exists a Cartesian coordinate system oxyz in which it is diagonal, having the form

$$\begin{aligned} \mathbf{I} = \left[ \begin{array}{c@{\quad }c@{\quad }c} I_x &{} 0 &{} 0 \\ 0 &{} I_y &{} 0 \\ 0 &{} 0 &{} I_z \end{array} \right] , \end{aligned}$$
(19)

where the coordinate axes are called the principal axes and the constants \(I_x,\, I_y\), and \(I_z\) are the principal moments of inertia.

We consider the body-attached frame that oxyz has its origin o at the center of mass of the solid particle and oxyz has the same orientation as the fixed frame OXYZ initially. Three axes of oxyz will be the principal axes during the rotation of the solid particle.

1.2 Rotation velocity

For a given torque \(\mathbf{T}_{\alpha }=\left[ \tau _x\,\tau _y\,\tau _z\right] \) on the solid particle, we applied the Euler equations to calculate the time derivative of the angular velocity \(\varvec{\Omega }_{\alpha }\) as

$$\begin{aligned} \begin{aligned} I_x\dot{\omega }_x + (I_z-I_y){\omega }_y{\omega }_z&= \tau _x,\\ I_y\dot{\omega }_y + (I_x-I_z){\omega }_z{\omega }_x&= \tau _y,\\ I_z\dot{\omega }_z + (I_y-I_x){\omega }_x{\omega }_y&= \tau _z, \end{aligned} \end{aligned}$$
(20)

where the torque is decomposed into three components along the principal axes. In a SPH simulation, \(\mathbf{T}_{\alpha }=\left[ \tau _x~\tau _y~\tau _z\right] \) is calculated from the summation of torque on all boundary particles, which constitute the solid particle. As we know the principal moments of inertia \(\mathbf{I}\) for a particular solid particle, we may use Eq. (20) to calculate \(\dot{\varvec{\Omega }}=[\dot{\omega }_x\,\dot{\omega }_y\,\dot{\omega }_z]\) and integrate \(\varvec{\Omega }\) thereafter.

1.3 Rotation frame

Equation (20) is for the case that oxyz frame coincides with the fixed OXYZ frame. In a real simulation, one has to use Eq. (20) on the rotating oxyz frame, which accumulates the history of all previous rotations, namely, \(\mathbf{ARM}^n = \mathbf{RM}^n \mathbf{RM}^{n-1} \ldots \mathbf{RM}^1 \mathbf{RM}^0\). Therefore,

$$\begin{aligned}&(o'x'y'z')=(oxyz)^n = \mathbf{ARM}^n (oxyz)^0 = \mathbf{ARM}^n (OXYZ), \nonumber \\\end{aligned}$$
(21)
$$\begin{aligned}&\varvec{\Omega }_{\alpha }^{\prime } = \mathbf{ARM}^n \varvec{\Omega }_{\alpha }, \end{aligned}$$
(22)
$$\begin{aligned}&\mathbf{T}_{\alpha }^{\prime } = \mathbf{ARM}^n \mathbf{T}_{\alpha }, \end{aligned}$$
(23)

where \(\varvec{\Omega }^{\prime }\) and \(\mathbf{T}^{\prime }\) are angular velocity and torque in \(o'x'y'z'\) frame and \(\varvec{\Omega }_{\alpha }\) and \(\mathbf{T}_{\alpha }\) are angular velocity and torque in OXYZ frame. Therefore in the rotating frame \(o'x'y'z'\)

$$\begin{aligned} \begin{aligned} I_x\dot{\omega }_x^{\prime } + (I_z-I_y){{\omega }_y^{\prime }}{\omega }_z^{\prime }&= \tau _x^{\prime },\\ I_y\dot{\omega }_y^{\prime } + (I_x-I_z){\omega }_z^{\prime }{\omega }_x^{\prime }&= \tau _y^{\prime },\\ I_z\dot{\omega }_z^{\prime } + (I_y-I_x){\omega }_x^{\prime }{\omega }_y^{\prime }&= \tau _z^{\prime }. \end{aligned} \end{aligned}$$
(24)

where we calculate \(\dot{\varvec{\Omega }}^{{\prime }n}_{\alpha }=[\dot{\omega }_x^{\prime }\,\dot{\omega }_y^{\prime }\,\dot{\omega }_z^{\prime }]^{n}\) at nth time step. Afterwards, we transfer back \(\dot{\varvec{\Omega }}^{{\prime }n}_{\alpha }\) into OXYZ frame by

$$\begin{aligned} \dot{\varvec{\Omega }}^{n}_{\alpha }=(\mathbf{ARM}^{n})^{-1}\dot{\varvec{\Omega }}^{{\prime }n}_{\alpha }=(\mathbf{ARM}^{n})^{t}{\varvec{\Omega }}^{{\prime }n}_{\alpha }. \end{aligned}$$
(25)

Appendix 2: Solution to the splitting system

In this appendix the solution to the linear system (10) is calculated. If the second equation is subtracted from the first we find

$$\begin{aligned} \tilde{\varvec{V}}_{\alpha \beta }\cdot \left[ \varvec{1} - A_{\alpha \beta }\varvec{e}_{\alpha \beta }\varvec{e}_{\alpha \beta }\right] = \varvec{V}_{\alpha \beta }^{\prime } \end{aligned}$$
(26)

where

$$\begin{aligned} A_{\alpha \beta }= & {} {\Delta } t_\mathrm{sweep} \left( \frac{1}{m_{\alpha }} + \frac{1}{m_{\beta }}\right) \nonumber \\&\times \left[ -6\pi \eta \left( \frac{a_{\alpha }a_{\beta }}{a_{\alpha } + a_{\beta }}\right) ^2\left( \frac{1}{s} - \frac{1}{s_c}\right) \right] \end{aligned}$$
(27)

The solution can be obtained by inverting the matrix of the system. This can be done by supposing that the inverse matrix is of the form \(\left[ Y_{\alpha \beta }\varvec{1} + Z_{\alpha \beta }\varvec{e}_{\alpha \beta }\varvec{e}_{\alpha \beta }\right] \), where \(Y_{\alpha \beta }\) and \(Z_{\alpha \beta }\) are the functions to be determined. Then, the solution to the system is

$$\begin{aligned} \tilde{\varvec{V}}_{\alpha \beta } = \varvec{V}_{\alpha \beta }^{\prime }\cdot \left[ \varvec{1} -\frac{A_{\alpha \beta }}{A_{\alpha \beta }-1} \varvec{e}_{\alpha \beta } \varvec{e}_{\alpha \beta }\right] . \end{aligned}$$
(28)

Given that the lubrication correction conserves linear momentum between pairs of particles, \(\tilde{\varvec{V}}_{\alpha }\) and \(\tilde{\varvec{V}}_{\beta }\) can be calculated as

$$\begin{aligned} \begin{aligned} \tilde{\varvec{V}}_{\beta }&= \frac{m_{\alpha }\left( \varvec{V}_{\alpha }^{\prime } - \tilde{\varvec{V}}_{\alpha \beta }\right) + m_\beta \varvec{V}_{\beta }'}{m_{\alpha } + m_{\beta }},\\ \tilde{\varvec{V}}_{\alpha }&= \tilde{\varvec{V}}_{\alpha \beta } + \tilde{\varvec{V}}_{\beta }. \end{aligned} \end{aligned}$$
(29)

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Vázquez-Quesada, A., Bian, X. & Ellero, M. Three-dimensional simulations of dilute and concentrated suspensions using smoothed particle hydrodynamics. Comp. Part. Mech. 3, 167–178 (2016). https://doi.org/10.1007/s40571-015-0072-5

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s40571-015-0072-5

Keywords

Navigation