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

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.


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 highresolution 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 N particles. 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 N body 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 w R and w D , 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 farfield 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 manybody 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., ε ¼ 1 where ε 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 particlefluid 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.

Dissipative hydrodynamics as unified particle dynamics
In this section, a unified view of preexisting particle dynamic method is discussed.

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 N molecules in volume V at temperature T. Newton's second laws of motion for this system is where F ij is a force exerted on particle j of mass m j from particle i of mass m i at time t. Position r i and velocity v i of molecule i are updated from time t to t þ δt, i.e., from its initial values of respectively, using the acceleration a i determined 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 surfacedeformation 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.

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].

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 m is the particle mass at position x, β is the drag coefficient, and f 0 is a random fluctuating force of zero mean: where T is the absolute temperature, k B is Boltzmann's constant, and δ t ð Þ is Dirac's delta function. Discarding the random force, one divides both sides of Eq. (2) by m to 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 a p is the radius of the primary particle [22].
Governing equation. A governing equation of DHD simulation is as follows: where M is a diagonal matrix of mass and moment of inertia; v and U are translational/rotational velocities of the particle and the fluid, respectively; Q p is the generalized interparticle and conservative force/torque vector; R is the grand resistance matrix; and B is the Brownian matrix of zero mean and finite variance: where δ t ð Þ is the Dirac-Delta function. And, dW is the Ito-Wiener process [23,24] having the following mathematical properties: Þ follows the normal distribution, and finally dW Á dW ¼ dt. The Brownian matrix B can be calculated by decomposing the grand resistance matrix such as where A is the decomposed matrix to be obtained and I is the identity matrix. Statistically, the the identity matrix can be expressed as where C is a vector with zero mean and unit variance, i.e., C h i ¼ 0 and and substituted into (10) to provide Therefore, B is obtained by calculating a square root of R matrix, which is equal to A matrix of Eq. (10). The identity relationship of Eq. (11) is not satisfied at specific time t but statistically by taking an average of C tr Á C for 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 0 ¼ dW=dt. Although A is deterministically calculated to satisfy Eq. (10), C tr Á C ¼ I is not valid at an instance but statistically. In the same sense, Eq. (11) is satisfied statistically.
Hydrodynamic tensors In Eq. (14), the generalized force requires a calculation of the grand resistance matrix R, which will allow to generate A. Consider particle i among N p particles in a given volume V, translating with a linear velocity v i and rotating with an angular velocity ω i at an instantaneous position r i t ð Þ. In the absence of particles, the fluid flow at the center of particle i can be represented as U ∞ r i ð Þ. At a point r ¼ x; y; z ð Þ∈ S i on surface S i of particle i from the particle center r i , the flow field is described as where U ∞ is the unidirectional velocity, Ω ∞ is the vorticity, and E ∞ is the rate of strain which are evaluated at r i . Because the rate-of-strain matrix is symmetric and traceless, the original nine components are reduced to 6

Advanced Computational Fluid Dynamics for Emerging Engineering Processes…
The disturbance velocity field at the particle surface S i is where u i and ω i are the translational and angular velocities of particle i, respectively. The translational/angular velocities and the rate of strain of particle i relative 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 j is moving with linear and angular velocities of u j and ω j under the influences of the ambient flow field characterized using U ∞ , Ω ∞ , and E ∞ , it experiences the hydrodynamic force F H and torque T H . The stresslet S H can 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, andF 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 N bodies, can be expressed as where S H is the hydrodynamic stresslet. The matrix μ ∞ (multiplied toF is called far-field grand mobility matrix. An inverse relationship of Eq. (25) is where the matrix R ∞ (multiplied to Δv j ; ÀE ∞ i Â Ã tr ) is the far-field grand resistance matrix as an inverse of μ ∞ , having the mathematical identity as where I is the identity matrix. Note that μ ∞ ð Þ À1 ¼ R ∞ and R ∞ ð Þ À1 ¼ μ ∞ for grand matrices but R ∞ Fv À Á À1 6 ¼ μ ∞ vF for sub-matrices. The grand resistance matrix R in Eq. (8) refers to R ∞ Fv of Eq. (26), and R ∞ FE Á E ∞ j is 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 forceF H by 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 r j with linear and angular velocities of v j and ω i , respectively, is the ambient flow field consisting of U ∞ ; ;Ω ∞ ; ÀE ∞ ð Þ at the particle location. Unlike MD and BD, DHD explicitly includes the sizes of individual spherical particles. In this case, r j indicates the position of the particle centroid or center of mass, and the fluid field is calculated at r j without 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.

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 U is the flow velocity, P is pressure, η and ρ are the viscosity and density of the fluid, respectively. In the adjacent space of particle j at r j , U can 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 U and 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 U; P ð Þare 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, r k , for k ¼ 1 À N p 2. To calculate the distance between r k and wall surfaces, if the particle is close to wall boundaries 3. To interpolate the flow field U ∞ ; ;Ω ∞ ; ÀE ∞ ð Þ at r k within 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.

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, p i , located at position r The volume of the tetrahedron can be calculated as where J is the Jacobian matrix given as Using this mathematical relationship, one can easily test whether a particle position r p ¼ x p ; y p ; z p is inside or outside the tetrahedron. If r p is 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: 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 r p within the tetrahedron, dividing both sides of Eq. (32) by V gives where, for i ¼ 1 À 4, which is solely determined by the internal position r p . A relationship between r p and ξ is where the inverse of the 4 Â 4 matrix can be analytically available. Eq. (36) assumes that r p is 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 Â 4 linear system, but Eq. (33) requires four times calculation of 4 Â 4 Jacobian matrices.
3.1.2 How to calculate a distance between wall surfaces to the position of my particle Using the four vertices of r 1 ; r 2 ; r 3 ; r 4 ð Þof a tetrahedron cell, one can identify four triangular surfaces of S 1 r 2 ; r 3 ; r 4 ð Þ , S 2 r 1 ; r 3 ; r 4 ð Þ , S 3 r 1 ; r 2 ; r 4 ð Þ , and S 4 r 1 ; r 2 ; r 3 ð Þ . For example, let us consider S 4 having vertices of r 1 ; r 2 ; r 3 ð Þ . 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 S 4 : The particle position relative to r 1 is r p=1 ¼ r p À r 1 . Then, the distance between surface S 4 and position r p is simply where the absolute value is necessary because n can direct inside or outside the tetrahedron volume, depending on choice of the reference position, r 1 . This wallparticle distance calculation is necessary when the surface S 4 is known as a wall boundary, so the distance d cannot be smaller than the particle radius a p . 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 k is found inside a tetrahedron (of index l), the flow field consisting of U ∞ ; ;Ω ∞ ; ÀE ∞ ð Þ needs 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 U ∞ i and ∂ j U i ∂U i =∂x j for i and j ¼ 1 À 3. (Note that pressure value is not required to calculate the generalized hydrodynamic forcesF H .) Four vertices of a cell can be represented as r l ¼ x l ; y l ; z l À Á for l ¼ 1 À 4. A scalar quantity of interest S at each vertex point can be denoted as S l . Then, the value of S at an arbitrary internal position r p ¼ x p ; y p ; z p can be calculated as a linear superposition of S l and ξ l : For DHD simulations, S p represents each element of the unidirectional velocity U ∞ , vortex Ω ∞ , and the rate of strain E ∞ , which are calculated using U ∞ i or its gradient ∂ j U i .

Collision rules
Suppose there are two particles colliding each other, which are particle i and j located at r i and r j , translating with v i and v j and rotating with ω i and ω j , respectively. Relative position of particle i with respect to particle j is defined as r ij ¼ r i À r j , and similarly the relative velocity of i to j is v ij ¼ v i À v j . A normal vector from j to i is denoted as After an instantaneous collision, the two particle have the following velocities [14]: 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.