Open access peer-reviewed chapter

A Coupling Algorithm of Computational Fluid and Particle Dynamics (CFPD)

Written By

Albert S. Kim and Hyeon-Ju Kim

Submitted: December 4th, 2018 Reviewed: May 16th, 2019 Published: August 2nd, 2019

DOI: 10.5772/intechopen.86895

Chapter metrics overview

783 Chapter Downloads

View Full Metrics


Computational fluid dynamics (CFD) and particle hydrodynamics (PHD) have been developed almost independently. CFD is classified into Eulerian and Lagrangian. The Eulerian approach observes fluid motion at specific locations in the space, and the Lagrangian approach looks at fluid motion where the observer follows an individual fluid parcel moving through space and time. In classical mechanics, particle dynamic simulations include molecular dynamics, Brownian dynamics, dissipated particle dynamics, Stokesian dynamics, and granular dynamics (often called discrete element method). Dissipative hydrodynamic method unifies these dynamic simulation algorithms and provides a general view of how to mimic particle motion in gas and liquid. Studies on an accurate and rigorous coupling of CFD and PHD are in literature still in a growing stage. This chapter shortly reviews the past development of CFD and PHD and proposes a general algorithm to couple the two dynamic simulations without losing theoretical rigor and numerical accuracy of the coupled simulation.


  • computational fluid dynamics (CFD)
  • computational fluid and particle dynamics (CFPD)
  • dissipative hydrodynamics
  • tetrahedron mesh
  • mesh interpolation

1. Introduction

The first simulations of a liquid were conducted at the Los Alamos National Laboratory in the early 1950s, using the Los Alamos computer, Mathematical Analyzer, Numerical Integrator, and Computer (MANIAC). This computer was developed under the direction of Nicholas Metropolis, who is the pioneer of the modern (Metropolis) Monte Carlo simulation [1, 2]. The first MC simulation was conducted using the Lennard-Jones potential to investigate the material properties of liquid argon [3]. In the MC simulations, the phase space was searched to find more probable thermodynamic states and calculate macroscopic material properties (i.e., experimentally observable) using averages of the same physical quantity over the micro-thermodynamic states. In principle, the MC method assumes that the system of interest is in a static equilibrium state, and therefore, the time evolution is replaced by the phase-space averaging. The Boltzmann factor is used as a transition probability between two energy states.

The first CFD paper was published by Hess and Smith [4]. Their method is known as panel method as the surface was discretized with many panels. More accurate CFD work of advanced panel method can be found elsewhere [5, 6]. In general, CFD research is categorized by Eulerian and Lagrangian approaches, which are grid-dependent and meshfree, respectively. In the Eulerian category, the CFD performance regarding numerical accuracy and computational speed depends on how to discretize the computational domain. Popular methods include finite volume method (FVM) [7], finite element method (FEM), finite difference method (FDM), spectral element method (SEM), boundary element (BEM), and high-resolution discretization schemes.

The deterministic molecular dynamic (MD) method treats a particle as a point mass, which has a finite mass and other physical properties but does not have any volume and shape. Given external and interparticle forces, positions and velocities in linear motion are predicted using Newton’s second law for Nparticles. When the motion of solute molecules in the solvent fluid is of concern, then the deterministic forces form solvent to solute molecules can be mathematically replaced by random fluctuating forces. These random forces have a zero mean over time, and its magnitude is determined by the dissipation fluctuation theorem [8].

Langevin’s equation includes the hydrodynamic drag balanced by the random forces due to thermal fluctuation [9]. Solving Nbody Langevin equation is called Brownian dynamics (BD) [10]. Although BD can include effects of particle sizes in the dynamic simulations, it is fundamentally limited to the low concentration of solutes due to the Oseen tensor (see Appendix for details). Brownian dynamics (BD) was initially developed to reduce computational load by replacing deterministic interactions between a solute molecule and many solvent molecules by randomly fluctuating, probabilistic interaction as a net driving force. BD presumes a dilute solution, which means that a mean distance between solutes is much longer than the size of the molecule. A large number of solvent molecules exist around solutes, which is enough to exert random forces due to a tremendous number of collisions. When BD is applied to the dynamic motion of multiple particles, the lack of hydrodynamic interactions may provide erroneous results because the fluctuation-dissipation principle is not quantitatively well balanced. This limitation of BD to a single particle or a dilute solution is at the equivalent level of Stokes’ drag coefficient, used to calculate the particle diffusivity.

Dissipative particle dynamics (DPD) is an updated version of BD, which specifically includes the interparticle hydrodynamic forces [11, 12]. Two functions, often denoted as wRand wD, are proposed to quantify the presumed relationship of pair-wise hydrodynamic forces/torques, determined by the dissipation fluctuation theorem. The proposed forms of the hydrodynamic interaction are vector wise as the real viscous forces are tensor-wise. Dissipative particle dynamics (DPD) updates BD by including an approximate form of tensor-wise hydrodynamic interactions as a pair-wise vector form. A force exerted on a particle has three types such as conservative, dissipative, and random forces, and DPD satisfies the fluctuation-dissipation theorem by balancing the dissipative and random forces.

Stokesian dynamics (SD) uses the grand mobility matrix to include the far-field hydrodynamic forces/torques [13]. In-depth comparative analyses of present particle dynamic methods can be found elsewhere [14]. This mobility matrix is inverted and updated by adding differences between near-field lubrication forces and far-field two-body interactions. As the lubrication force is proportional to the logarithm of the surface-to-surface distance between two particles, it diverges during events of particle collisions. Besides, if two particles are close to each other, i.e., the surface-to-surface distance is much smaller than the particle diameter, the SD formalism of near-field hydrodynamic contribution to the global hydrodynamic resistance (i.e., mobility inverted) often overestimates the real hydrodynamic repulsion forces. One of the primary reasons is that the lubrication theory assumes that particle surfaces are perfectly smooth, while in reality particle surfaces are rough and therefore surface friction also plays an essential role during collision events. Stokesian dynamics includes the most accurate estimation of the hydrodynamic tensors using the grand mobility matrix and pair-wise lubrication interactions. The fluctuation-dissipation theorem is better satisfied in SD by using many-body hydrodynamic interactions. In SD simulations, the fluid flow is given at any point in the computational domain as a combination of the unidirectional velocity, vortex velocity, and the rate of strain that is traceless and symmetric. The feedback interaction between a particle and fluid is already included in the expansion of Faxen rule. Due to the logarithmic characteristics of the lubrication forces, if two particles touch each other, they experience infinite repulsive forces. This characteristic of the lubrication force indicates the hard-sphere behavior of the colliding particles conceptually, implicitly assuming a perfectly elastic interparticle collision i.e., ε=1where εis the coefficient of restitution.

Dissipative hydrodynamics (DHD) overcomes significant limitations of the particle dynamic method discussed above. DHD is a generalized method of which special cases converge to MD, BD, DPD, and SD by turning on or off specific force mechanisms. Details of DHD can be found elsewhere [14, 15, 16]. CFD and PHD were developed and applied without strong mutual influences. Particle tracking method can be viewed as a reasonable way to investigate the hydrodynamic motion of particles under the influence of ambient fluid flow. But, it has several fundamental limitations by neglecting particle density, particle shapes and sizes, and particle-fluid interactions. More importantly, the basic two-body interactions due to collisions, viscous flow, and electromagnetic properties are not included, and dispersion forces were dropped. In this light, the particle tracking method does not track particles but fluid elements moving along streamlines. While theories of particle hydrodynamics are not rigorously applied to engineering processes, this chapter includes a possible method to couple CFD and DHD in a seamless, robust way.


2. Dissipative hydrodynamics as unified particle dynamics

In this section, a unified view of preexisting particle dynamic method is discussed.

2.1 Overview

In the deterministic simulations, particle dynamics can be classified based on sizes of objects of interests. The purposes of particle dynamics are to provide the exact solution of a complex problem, bridging the theory and experiments. Note that the fundamental principles often provide governing equations, which were proven to be valid. It is difficult to solve a governing equation of a problem, if it has complex geometry and coupled boundary conditions. Experimental observations show natural processes using quantified information. Dissipative hydrodynamics is a generalized algorithm that unifies most of the preexisting particle dynamic simulation algorithms [16, 17].

At microscale of an order of nanometers, molecular dynamics deals with the motion of many molecules in various phases such as gas, liquid, and solid. Suppose that a system contains Nmolecules in volume Vat temperature T. Newton’s second laws of motion for this system is


where Fijis a force exerted on particle jof mass mjfrom particle iof mass miat time t. Position riand velocity viof molecule iare updated from time tto t+δt, i.e., from its initial values of rit=0=ri0and vit=0=vi0, respectively, using the acceleration aidetermined using Eq. (1). Numerical evolution of Eq. (1) requires a specific algorithm for double integration [18, 19, 20, 21].

A macroscale of order of millimeters, granular dynamics often called the discrete element method (DEM) includes specifically collision rules using the restitution and friction coefficients. During inelastic collisions of nonrotating particles, the particle kinetic energy is continuously lost, and their motion is decelerated. For a collision of rotating spheres, the surface friction provides an effective torque (as action and reaction) of the same magnitude and opposite directions to two colliding spheres. Rotational motion of a non-touching particle in a fluid medium generates angular dissipation of kinetic energy. Considering granules, i.e., non-Brownian particles, in a gas phase often neglect solute molecules and approximate the system as multiple particles undergoing the gravitational force field in a vacuum phase. As granular particle mass is much higher than that of colloids or nanoparticles in an aqueous solution, in a stationary phase the gravitational force is often balanced by normal forces developed at interfaces of particles to touching neighbors or a rigid wall. Implementation of the hydrodynamic lubrication interactions to granular particles is a difficult task, which requires an in-depth understanding of microscopic surface-deformation phenomena, linked to macroscopic particle motion.

A universal simulation method that can seamlessly include forces/torques exerted on arbitrary particles is therefore of great necessity. The method, first of all, should be able to:

  1. Investigate the accelerating/decelerating motion of particles.

  2. Satisfy the fluctuation-dissipation theorem for Brownian particles.

  3. Include many-body hydrodynamic interactions.

  4. Mimic inelastic collisions between spherical particles.

  5. Apply constraint forces to form a nonspherical rigid body consisting of unequal spherical particles.

Dissipative hydrodynamics (DHD) has all the features required to be a universal simulation method for particle dynamics by taking specific advantages from MD, BD, SD, DPD, and DEM. A detailed review can be found elsewhere [14].

2.2 DHD formalism

Particle relaxation time. For a single particle motion in a viscous fluid, the governing equation can be in 1D space for simplicity:


where mis the particle mass at position x, βis the drag coefficient, and fis a random fluctuating force of zero mean:


where Tis the absolute temperature, kBis Boltzmann’s constant, and δtis Dirac’s delta function. Discarding the random force, one divides both sides of Eq. (2) by mto have


which indicates that m/βhas a dimension of time. This time scale is called the particle relaxation time, define as


which is a time scale that after a particle noticeably slows down after it starts moving with an initial velocity under the drag force. Stokes derived the drag coefficient


where μis the absolute fluid viscosity and apis the radius of the primary particle [22].

Governing equation. A governing equation of DHD simulation is as follows:


where Mis a diagonal matrix of mass and moment of inertia; vand Uare translational/rotational velocities of the particle and the fluid, respectively; Qpis the generalized interparticle and conservative force/torque vector; Ris the grand resistance matrix; and Bis the Brownian matrix of zero mean and finite variance:


where δtis the Dirac-Delta function. And, dWis the Ito-Wiener process [23, 24] having the following mathematical properties: Wk=0at t=0, Wktis continuous, dWkWkt+δtWktfollows the normal distribution, and finally dWdW=dt. The Brownian matrix Bcan be calculated by decomposing the grand resistance matrix such as


where Ais the decomposed matrix to be obtained and Iis the identity matrix. Statistically, the the identity matrix can be expressed as


where Cis a vector with zero mean and unit variance, i.e., C=0and C2=1, respectively. The Brownian matrix is defined as


and substituted into (10) to provide


Therefore, Bis obtained by calculating a square root of Rmatrix, which is equal to Amatrix of Eq. (10). The identity relationship of Eq. (11) is not satisfied at specific time tbut statistically by taking an average of CtrCfor a much longer period than the particle relaxation time τp. The effective force acting on a particle of a swarm of many particles in a viscous fluid is then represented as


where W=dW/dt. Although Ais deterministically calculated to satisfy Eq. (10), CtrC=Iis not valid at an instance but statistically. In the same sense, Eq. (11) is satisfied statistically.

Hydrodynamic tensorsIn Eq. (14), the generalized force requires a calculation of the grand resistance matrix R, which will allow to generate A. Consider particle iamong Npparticles in a given volume V, translating with a linear velocity viand rotating with an angular velocity ωiat an instantaneous position rit. In the absence of particles, the fluid flow at the center of particle ican be represented as Uri. At a point r=xyzSion surface Siof particle ifrom the particle center ri, the flow field is described as


where Uis the unidirectional velocity, Ωis the vorticity,


and Eis the rate of strain


which are evaluated at ri. Because the rate-of-strain matrix is symmetric and traceless, the original nine components are reduced to


The disturbance velocity field at the particle surface Siis


where uiand ωiare the translational and angular velocities of particle i, respectively. The translational/angular velocities and the rate of strain of particle irelative to the ambient flow field have then 11 degrees of freedom such as


For non-Brownian particles, the governing Eq. (8) is reduced back to that of Stokesian dynamics, which is Langevin’s equation with the constant drag coefficient β, replaced by the grand resistant matrix R.

When particle jis moving with linear and angular velocities of ujand ωjunder the influences of the ambient flow field characterized using U, Ω, and E, it experiences the hydrodynamic force FHand torque TH. The stresslet SHcan be obtained but does not directly contribute to the particle acceleration. The generalized velocity and force are related through the grand mobility matrix μ. Here, we use q, q̇, and F˜for generalized coordinates, velocities, and forces, respectively:


The generalized relative velocity is


for both translational and angular motion. Then, the hydrodynamic interactions, i.e., forces and torques exerted on Nbodies, can be expressed as


where SHis the hydrodynamic stresslet. The matrix μ(multiplied to F˜iHSiHtr) is called far-field grand mobility matrix. An inverse relationship of Eq. (25) is


where the matrix R(multiplied to ΔvjEitr) is the far-field grand resistance matrix as an inverse of μ, having the mathematical identity as


where Iis the identity matrix. Note that μ1=Rand R1=μfor grand matrices but RFv1μvFfor sub-matrices. The grand resistance matrix Rin Eq. (8) refers to RFvof Eq. (26), and RFEEjis an extra forcing term due to the rate of strain.

Note that at time t, the right-hand side of Eq. (25) is known, and one can calculate the generalized hydrodynamic force F˜Hby using an appropriate solver in numerical linear algebra [25, 26]. As noted in Eq. (15), the required information to evolve the motion of particle j, located at rjwith linear and angular velocities of vjand ωi, respectively, is the ambient flow field consisting of U;ΩEat the particle location. Unlike MD and BD, DHD explicitly includes the sizes of individual spherical particles. In this case, rjindicates the position of the particle centroid or center of mass, and the fluid field is calculated at rjwithout considering the presence of the particles. Therefore, the calculation of the ambient flow field is highly dependent on a mesh structure used for CFD simulations.


3. Computational fluid dynamics coupled with dissipative hydrodynamics

A general governing equation for fluid dynamics is Navier–Stokes equation, of which most general form for incompressible fluid is


where Uis the flow velocity, Pis pressure, ηand ρare the viscosity and density of the fluid, respectively. In the adjacent space of particle jat rj, Ucan be expanded as expressed in Eq. (15). Most fluid dynamic problems have at least three boundaries, which are the inlet, outlet, and side walls. For inlet and outlet surfaces, Neuman and Dirichlet or mixed boundary conditions are often used to set values or conditions of Uand P. On the wall, zero velocity and zero-gradient pressure are usually assigned. The former condition assumes that there are strong adhesion forces between solvent molecules and wall surfaces. The solvent molecules are fixed on the wall. The velocities of the wall-adsorbed solvent molecules are equal to those of the solid walls, which is zero for non-moving walls. Values of UPare calculated at grid points in internal spaces surrounded by the boundary surfaces. The simulation accuracy and numerical convergence highly depend on the mesh structures. When CFD is coupled with particle dynamics, which is in this study, DHD, there are additional requirements that should be satisfied to evolve the motion of multi-particles moving in a fluid flow:

  1. To identify a cell of the constructed mesh grid, which contains the kth particle’s position, rk, for k=1Np

  2. To calculate the distance between rkand wall surfaces, if the particle is close to wall boundaries

  3. To interpolate the flow field U;ΩEat rkwithin a cell that contains kth particle

Possible methods to satisfy the three requirements are dependent on available CFD solvers and flexibility of applying customized modifications. Here we suggest fundamental approaches to meet the requirement numerically.

3.1 Geometric calculations

3.1.1 How to determine a cell containing the centroid position of my particle

Computational grid cells have various structures such as hexahedron, wedge, prism, pyramid, tetrahedron, and tetrahedron wedge. Among them, hexahedron followed by tetrahedron structures is widely used to generate mesh structure of bulk (internal) spaces. Cubic and rectangular shapes are representative structures of hexahedrons, consisting of eight vertices (points) and six rectangular (or square) surfaces. On the other hand, a tetrahedron cell has only four vertices (as compared to eight in hexahedron) and three triangular surfaces. Each of these two cell structures has its advantages and disadvantages in CFD simulations. Formation of hexahedron meshes is straightforward, and numerical solutions are well converged to provide accurate results within a tolerable error, especially if edge lines are well aligned to the flow directions. However, if a computational domain includes complex and curved surface structures such as human faces or globes, the hexahedron meshes often provide unrealistic small exuberances instead of well-curved surfaces. As three triangular surfaces surround the tetrahedron volume, it can form well-fitted boundary layers of arbitrary shapes. As edges of tetrahedrons cannot be fully aligned on a straight line, the numerical convergence of tetrahedron meshes is often more sensitive to the fineness of generated mesh structures than that of hexahedron meshes. To overcome this limitation, one can make tetrahedron meshes often finer than that of hexahedron meshes to solve the same problem with similar accuracy. Nevertheless, there are unique mathematical advantages of using tetrahedron meshes for coupled simulations of CFD and particle hydrodynamics.

Location test using volume calculation. A tetrahedron have four points, pi, located at position rifor i=14, where ri=xiyizi. The volume of the tetrahedron can be calculated as


where Jis the Jacobian matrix given as


Using this mathematical relationship, one can easily test whether a particle position rp=xpypzpis inside or outside the tetrahedron. If rpis inside of the tetrahedron, then the total volume of the tetrahedron can be divided by four small pieces and their sum is equal to V:




where, for example, Jp,2,3,4means the Jacobian matrix of Eq. (30) with the first column 1x1y1z1trreplaced by 1xpypzptr. Conversely, if rpis outside of the tetrahedron, we have an inequality of


If a particle is very close to one of the triangular surfaces of the tetrahedron, then the inequality check of Eq. (33) may not be done accurately. In this case, this location check can be extended to the nearest neighbor cells, especially one that shares the triangular surface that the particle is closely located.

For rpwithin the tetrahedron, dividing both sides of Eq. (32) by Vgives


where, for i=14,


which is solely determined by the internal position rp. A relationship between rpand ξis


where the inverse of the 4×4matrix can be analytically available. Eq. (36) assumes that rpis located within the tetrahedron, but it can also use the particle location check as a better alternative to Eq. (33) because Eq. (36) requires one time solving of a 4×4linear system, but Eq. (33) requires four times calculation of 4×4Jacobian matrices.

3.1.2 How to calculate a distance between wall surfaces to the position of my particle

Using the four vertices of r1r2r3r4of a tetrahedron cell, one can identify four triangular surfaces of S1r2r3r4, S2r1r3r4, S3r1r2r4, and S4r1r2r3. For example, let us consider S4having vertices of r1r2r3. If we calculate relative position of vertices 2 and 3 with respect to vertex 1, denoted as


their cross product allows us to calculate a unit vector normal to surface S4:


The particle position relative to r1is rp/1=rpr1. Then, the distance between surface S4and position rpis simply


where the absolute value is necessary because ncan direct inside or outside the tetrahedron volume, depending on choice of the reference position, r1. This wall-particle distance calculation is necessary when the surface S4is known as a wall boundary, so the distance dcannot be smaller than the particle radius ap. After a single-step time evolution, if a particle is found overlapped with a wall surface, then collision rule will be applied to return the particle to the tetrahedron interior.

3.1.3 How to interpolate the flow field at the position of my particle

If a particle kis found inside a tetrahedron (of index l), the flow field consisting of U;ΩEneeds to be interpolated using values calculated adjacent locations. A proper choice of a set of these locations are the vertex points of the containing cell. Based on the definition of Ωand E, the basic quantities needed are Uiand jUiUi/xjfor iand j=13. (Note that pressure value is not required to calculate the generalized hydrodynamic forces F˜H.)

Four vertices of a cell can be represented as rl=xlylzlfor l=14. A scalar quantity of interest Sat each vertex point can be denoted as Sl. Then, the value of Sat an arbitrary internal position rp=xpypzpcan be calculated as a linear superposition of Sland ξl:


For DHD simulations, Sprepresents each element of the unidirectional velocity U, vortex Ω, and the rate of strain E, which are calculated using Uior its gradient jUi.

3.2 Collision rules

Suppose there are two particles colliding each other, which are particle iand jlocated at riand rj, translating with viand vjand rotating with ωiand ωj, respectively. Relative position of particle iwith respect to particle jis defined as rij=rirj, and similarly the relative velocity of ito jis vij=vivj. A normal vector from jto iis denoted as


After an instantaneous collision, the two particle have the following velocities [14]:


where 0εn1and 1εt1are restitution coefficients in the normal and tangential directions, respectively, μijis a reduced mass defined as


and gijis the relative velocity at the point of contact defined as


whose normal and tangential components are


If a collision between a wall surface and particle ioccurs, then the wall can be represented as a stationary spherical particle jhaving infinite mass and radius, i.e., mj, vj0, and aj:




and nis the normal vector of the colliding wall surface inward to the liquid volume. During the collision of particle iwith the wall, the wall is not moving so that ωj'=ωj=0and vj=vj=0.


4. Concluding remarks

Each of computational fluid dynamics and particle hydrodynamics is a challenging research topic, as applied to real engineering problems. Movements of particles (viewed as small solid pieces) in a moving fluid require rigorous interfaces to couple the two strongly correlated dynamic events. When the ambient flow pushes suspended particles in a liquid, dynamic responses of particles to exerting fluid change the fluid motion at the next time step, which returns to the particles with modified magnitude and direction. This fluid-particle (or fluid–solid) interaction is under the regime of Newton’s third law, i.e., action and reaction. Multi-body simulations including the fluid-particle interaction are generally a difficult task. In this chapter, we briefly reviewed the CFD and PHD literature and discussed a feasible method to simulate the coupled fluid-particle motion within a reasonable time duration. A tetrahedron-based mesh is proposed to take the mathematical advantages of tetrahedron structure, consisting of four vertices and four triangular surfaces. Advantages of tetrahedron meshes include the following features: first, to test a location of particle within or outside a tetrahedron mesh-cell; second, to calculate a distance between a particle surface and a wall surface; third, to interpolate a fluid velocity and its gradient using those of values given at vertex locations; and finally to track each particle from one cell to the other by using a pre-built list of nearest neighbor cells. As DHD uses the SD algorithm for hydrodynamics of non-Brownian particles and Ito-Weiner process for random fluctuating forces, it can be used as a general particle hydrodynamic simulation method when it is coupled with CFD using specific mesh structures. Hexahedron-based meshes can be used for the same purpose with the intrinsic advantages of aligning grid edges to estimated streamline directions. When particles move in a channel of complex geometry, boundary surfaces can be better constructed using tetrahedron meshes. Open-sourced meshes include gmsh, tetgen, and netgen, which can import structure files, generate meshes, and export them to a CFD solver package. The current coupling algorithm of CFD and DHD is limited to cases that particle Reynolds number does not exceed 1.0, but this restriction can be avoided by considering dominant forces/torques exerted on particles and simulation time intervals as compared to the particle relaxation time. The new coupling method covered in this chapter may provide a new foundation in a coupled simulation of CFD and DHD including DEM.



This work was financially supported by the National R&D project of “Development of 1MW OTEC demonstration plant (4/5)” (PMS4080) funded by the Ministry of Oceans and Fisheries of the Republic of Korea.


A. Appendix


A.1 The Oseen tensor and Faxen law

The Navier–Stokes equation for a laminar flow is from Eq. (28):


where δris the Dirac-Delta function, which indicates


where Vis a volume enclosing the origin r=0. The fundamental solution for vand pare


The Oseen tensor Gxfor the fluid velocity is given by


which is independent of fluid properties. The Oseen tensor for the pressure pis


where Pjis a constant at the ambient condition. The Faxen laws determine the hydrodynamic force and torque, especially, exerted on a sphere of radius a, moving with the linear and angular velocities of uand ω:


which indicates that the Stokes flow requires a quadrupole a2F2δxin addition to a monopole of 6πηau, which is a drag on a sphere undergoing steady translation.


  1. 1. Metropolis N, Ulam S. The Monte Carlo method. Journal of the American Statistical Association. 1949;44(247):335-341
  2. 2. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. The Journal of Chemical Physics. 1953;21(6):1087-1092
  3. 3. Wood WW, Parker FR. Monte carlo equation of state of molecules interacting with the lennard-jones potential. i. a supercritical isotherm at about twice the critical temperature. The Journal of Chemical Physics. 1957;27(3):720-733
  4. 4. Hess J, Smith A. Calculation of potential flow about arbitrary bodies. Progress in Aerospace Sciences. 1967;8:1-138
  5. 5. Rubbert P, Saaris G. Review and evaluation of a three-dimensional lifting potential flow computational method for arbitrary configurations. In: 10th Aerospace Sciences Meeting. Reston, VA: American Institute of Aeronautics and Astronautics; 1972
  6. 6. Maskew B. Prediction of subsonic aerodynamic characteristics: A case for low-order panel methods. Journal of Aircraft. 1982;19(2):157-163
  7. 7. Hyman JM, Knapp RJ, Scovel JC. High order finite volume approximations of differential operators on nonuniform grids. Physica D: Nonlinear Phenomena. 1992;60(1–4):112-138
  8. 8. Kubo R. The fluctuation-dissipation theorem. Reports on Progress in Physics. 1966;29:255-284
  9. 9. Langevin P. Sur la the’orie du mouvement brownien. Comptes Rendus de l'Académie des Sciences (Paris). 1908;146:530-533
  10. 10. Ermak DL, McCammon JA. Brownian dynamics with hydrodynamic interactions. The Journal of Chemical Physics. 1978;69(4):1352-1360
  11. 11. Hoogerbrugge P, Koelman J. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhysics Letters. 1992;19:155
  12. 12. Koelman J, Hoogerbrugge P. Dynamic simulations of hard-sphere suspensions under steady shear. Europhysics Letters. 1993;21:363
  13. 13. Brady JF, Bossis G. Stokesian dynamics. Annual Review of Fluid Mechanics. 1988;20:111-157
  14. 14. Kim AS, Kim H-J. dissipative dynamics of granular materials. In: Granular Materials. Rijeka: InTechOpen; 2017. pp. 9-42
  15. 15. Kim AS. Constraint dissipative hydrodynamics (HydroRattle) algorithm for aggregate dynamics. Chemistry Letters. 2012;41(10):1285-1287
  16. 16. Kim AS. Dissipative hydrodynamics of rigid spherical particles. Chemistry Letters. 2012;41(10):1128-1130
  17. 17. Kim AS, Kang S-T. Microhydrodynamics simulation of single-collector granular filtration. Chemistry Letters. 2012;41(10):1288-1290
  18. 18. Verlet L. Computer ‘experiments’ on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Physics Review. 1967;159(1):98-103
  19. 19. Martys NS, Mountain RD. Velocity Verlet algorithm for dissipative-particle-dynamics-based models of suspensions. Physical Review E. 1999;59(3):3733 LP-3733736
  20. 20. Hockney RW, Eastwood JW. Computer Simulation Using Particles. New York: Adam Hilger; 1988
  21. 21. Hockney RW. The potential calculation and some applications. In: Methods in Computational Physics. Vol. 9. New York: Academic Press; 1970
  22. 22. Stokes GG. On the effect of internal friction of fluids on the motion of pendulums. Transactions of the Cambridge Philosophical Society. 1851;9:1-106
  23. 23. Wiener N. Differential space. Journal of Mathematical Physics. 1923;58:31-174
  24. 24. Ito M. An extension of nonlinear evolution equations of the K-dV (mK-dV) type to higher orders. Journal of the Physical Society of Japan. 1980;49(2):771-778
  25. 25. Anderson E. LAPACK Users’ Guide. Philadelphia, PA: Society for Industrial and Applied Mathematics; 1987
  26. 26. Blackford LS, Choi J, Cleary A, D’Azevedo E, Demmel J, Chillon I, et al. ScaLAPACK User’s Guide. Philadelphia, PA: Society for Industrial and Applied Mathematics; 1997

Written By

Albert S. Kim and Hyeon-Ju Kim

Submitted: December 4th, 2018 Reviewed: May 16th, 2019 Published: August 2nd, 2019