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
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 . 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 . 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) , 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 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 .
Langevin’s equation includes the hydrodynamic drag balanced by the random forces due to thermal fluctuation . Solving body Langevin equation is called Brownian dynamics (BD) . 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 and , 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 . In-depth comparative analyses of present particle dynamic methods can be found elsewhere . 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., 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 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.
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 molecules in volume at temperature . Newton’s second laws of motion for this system is
where is a force exerted on particle of mass from particle of mass at time . Position and velocity of molecule are updated from time to , i.e., from its initial values of and , respectively, using the acceleration 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 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:
Investigate the accelerating/decelerating motion of particles.
Satisfy the fluctuation-dissipation theorem for Brownian particles.
Include many-body hydrodynamic interactions.
Mimic inelastic collisions between spherical particles.
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 .
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 is the particle mass at position , is the drag coefficient, and is a random fluctuating force of zero mean:
where is the absolute temperature, is Boltzmann’s constant, and is Dirac’s delta function. Discarding the random force, one divides both sides of Eq. (2) by to have
which indicates that 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 is the radius of the primary particle .
Governing equation. A governing equation of DHD simulation is as follows:
where is a diagonal matrix of mass and moment of inertia; and are translational/rotational velocities of the particle and the fluid, respectively; is the generalized interparticle and conservative force/torque vector; is the grand resistance matrix; and is the Brownian matrix of zero mean and finite variance:
where is the Dirac-Delta function. And, is the Ito-Wiener process [23, 24] having the following mathematical properties: at , is continuous, follows the normal distribution, and finally . The Brownian matrix can be calculated by decomposing the grand resistance matrix such as
where is the decomposed matrix to be obtained and is the identity matrix. Statistically, the the identity matrix can be expressed as
where is a vector with zero mean and unit variance, i.e., and , respectively. The Brownian matrix is defined as
and substituted into (10) to provide
Therefore, is obtained by calculating a square root of matrix, which is equal to matrix of Eq. (10). The identity relationship of Eq. (11) is not satisfied at specific time but statistically by taking an average of for a much longer period than the particle relaxation time . The effective force acting on a particle of a swarm of many particles in a viscous fluid is then represented as
Hydrodynamic tensors In Eq. (14), the generalized force requires a calculation of the grand resistance matrix , which will allow to generate . Consider particle among particles in a given volume , translating with a linear velocity and rotating with an angular velocity at an instantaneous position . In the absence of particles, the fluid flow at the center of particle can be represented as . At a point on surface of particle from the particle center , the flow field is described as
where is the unidirectional velocity, is the vorticity,
and is the rate of strain
which are evaluated at . 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 is
where and are the translational and angular velocities of particle , respectively. The translational/angular velocities and the rate of strain of particle 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 .
When particle is moving with linear and angular velocities of and under the influences of the ambient flow field characterized using , , and , it experiences the hydrodynamic force and torque . The stresslet 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 , , and 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 bodies, can be expressed as
where is the hydrodynamic stresslet. The matrix (multiplied to ) is called far-field grand mobility matrix. An inverse relationship of Eq. (25) is
where the matrix (multiplied to ) is the far-field grand resistance matrix as an inverse of , having the mathematical identity as
where is the identity matrix. Note that and for grand matrices but for sub-matrices. The grand resistance matrix in Eq. (8) refers to of Eq. (26), and is an extra forcing term due to the rate of strain.
Note that at time , the right-hand side of Eq. (25) is known, and one can calculate the generalized hydrodynamic force 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 , located at with linear and angular velocities of and , respectively, is the ambient flow field consisting of at the particle location. Unlike MD and BD, DHD explicitly includes the sizes of individual spherical particles. In this case, indicates the position of the particle centroid or center of mass, and the fluid field is calculated at 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.
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 is the flow velocity, is pressure, and are the viscosity and density of the fluid, respectively. In the adjacent space of particle at , 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 and . 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 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:
To identify a cell of the constructed mesh grid, which contains the th particle’s position, , for
To calculate the distance between and wall surfaces, if the particle is close to wall boundaries
To interpolate the flow field at within a cell that contains th 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, , located at position for , where . The volume of the tetrahedron can be calculated as
where is the Jacobian matrix given as
Using this mathematical relationship, one can easily test whether a particle position is inside or outside the tetrahedron. If 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 :
where, for example, means the Jacobian matrix of Eq. (30) with the first column replaced by . Conversely, if is 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 within the tetrahedron, dividing both sides of Eq. (32) by gives
where, for ,
which is solely determined by the internal position . A relationship between and is
where the inverse of the matrix can be analytically available. Eq. (36) assumes that 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 linear system, but Eq. (33) requires four times calculation of Jacobian matrices.
3.1.2 How to calculate a distance between wall surfaces to the position of my particle
Using the four vertices of of a tetrahedron cell, one can identify four triangular surfaces of , , , and . For example, let us consider having vertices of . 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 :
The particle position relative to is . Then, the distance between surface and position is simply
where the absolute value is necessary because can direct inside or outside the tetrahedron volume, depending on choice of the reference position, . This wall-particle distance calculation is necessary when the surface is known as a wall boundary, so the distance cannot be smaller than the particle radius . 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 is found inside a tetrahedron (of index ), the flow field consisting of 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 , the basic quantities needed are and for and . (Note that pressure value is not required to calculate the generalized hydrodynamic forces .)
Four vertices of a cell can be represented as for . A scalar quantity of interest at each vertex point can be denoted as . Then, the value of at an arbitrary internal position can be calculated as a linear superposition of and :
For DHD simulations, represents each element of the unidirectional velocity , vortex , and the rate of strain , which are calculated using or its gradient .
3.2 Collision rules
Suppose there are two particles colliding each other, which are particle and located at and , translating with and and rotating with and , respectively. Relative position of particle with respect to particle is defined as , and similarly the relative velocity of to is . A normal vector from to is denoted as
After an instantaneous collision, the two particle have the following velocities :
where and are restitution coefficients in the normal and tangential directions, respectively, is a reduced mass defined as
and is 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 occurs, then the wall can be represented as a stationary spherical particle having infinite mass and radius, i.e., , , and :
and is the normal vector of the colliding wall surface inward to the liquid volume. During the collision of particle with the wall, the wall is not moving so that and .
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.1 The Oseen tensor and Faxen law
The Navier–Stokes equation for a laminar flow is from Eq. (28):
where is the Dirac-Delta function, which indicates
where is a volume enclosing the origin . The fundamental solution for and are
The Oseen tensor for the fluid velocity is given by
which is independent of fluid properties. The Oseen tensor for the pressure is
where is a constant at the ambient condition. The Faxen laws determine the hydrodynamic force and torque, especially, exerted on a sphere of radius , moving with the linear and angular velocities of and :
which indicates that the Stokes flow requires a quadrupole in addition to a monopole of , which is a drag on a sphere undergoing steady translation.