Open Access is an initiative that aims to make scientific research freely available to all. To date our community has made over 100 million downloads. It’s based on principles of collaboration, unobstructed discovery, and, most importantly, scientific progression. As PhD students, we found it difficult to access the research we needed, so we decided to create a new Open Access publisher that levels the playing field for scientists across the world. How? By making research easy to access, and puts the academic needs of the researchers before the business interests of publishers.

We are a community of more than 103,000 authors and editors from 3,291 institutions spanning 160 countries, including Nobel Prize winners and some of the world’s most-cited researchers. Publishing on IntechOpen allows authors to earn citations and find new collaborators, meaning more people see your work not only from your own field of study, but from other related fields too.

Granules are inelastic particles, undergoing dissipative and repulsive forces on contact. A granular state consists of a conglomeration of discrete, non-Brownian particles in a combined state of solid, liquid, and gas. Modern theoretical physics lacks general theories for the granular states. Simulation methods for particle dynamics include molecular dynamics (MD), Brownian dynamics (BD), Stokesian dynamics (SD), dissipative particle dynamics (DPD), and dissipative hydrodynamics (DHD). These conventional methods were originally designed to mimic the small-particle motion being less influenced by the gravitational force. There are three reasons that a conventional method cannot be directly applied to investigate granular dynamics. First, volume exclusion forces between colliding particles are often disregarded due to strong repulsive forces between negatively charged colloids and nanoparticles. Second, the gravitational force is not significant as applied to small, light particles, and therefore it is often discarded in force/torque calculations. Third, energy conservation in an equilibrium state is not guaranteed for the granular system due to the inelastic and frictional nature of the granular materials. In this light, this chapter discusses the fundamentals of particle dynamics methods, formulates a robust theoretical framework for granular dynamics, and discusses the current applications and future directions of computational granular dynamics.

Civil and Environmental Engineering, University of Hawaii at Manoa, Honolulu, HI, United States of America

Hyeon-Ju Kim

Seawater Utilization Plant Research Center, Korea Research Institute of Ships and Ocean Engineering, Goseong-gun, Gangwon-do, Republic of Korea

*Address all correspondence to: albertsk@hawaii.edu

1. Introduction

1.1. Mechanics

Mechanicsis the investigation of physical bodies when they are subjected to forces and torques in Euclidean three-dimensional spaces. It is often referred to as classical mechanics (in physics), which is closely related to engineering mechanics. Classical mechanics has two branches: statics and dynamics. Staticsis concerned with the equilibrium of a body that is either at rest or moves with a constant velocity. Dynamicsdeals with the accelerated motion of a body, classified into two parts: kinematics, which treats only the geometric aspect of motion, and kinetics, which analyzes the forces causing the motion. Two representative objects in mechanics are a particle and a rigid body. A particle is the most basic unit of matter, which contains a mass of a negligible volume. Since the particle is small enough to be regarded as a point mass, its angular motion is completely discarded in analyzing its dynamics. The total energy of particles depends on their velocities and positions as influenced by external and internal forces. The mass of each particle is assumed invariant and therefore energy conservation is independent of mass conservation in classical mechanics. If particles of interest have sub-atomic sizes (such as hydrogens and electrons), classical mechanics fails to predict their intrinsic duality behaviors. Quantum mechanics explains matter’s simultaneous wave-like and particle-like properties, and unifies matter and energy as they converge at the level of Planck’s constant. A particle’s position and velocity cannot be measured accurately at a particular moment because if one measures the position accurately, then the particle’s momentum will be disturbed, and vice versa. This is called the Heidelberg uncertainty principle. In the macroscopic engineering world, in which humans observe objects with their naked eyes, quantum phenomena are extremely rare to observe. Although quantum mechanics includes classical mechanics as a sub-set, most conventional engineering phenomena are macroscopic enough to neglect sub-atomic effects. Based on the above-mentioned characteristics and classifications of mechanics, the granular dynamics in this chapter focuses on kinetics and hydrodynamics of multiple non-Brownian particles in locally confined three-dimensional (3D) spaces, often filled with fluid media.

1.2. Principles in classical mechanics

1.2.1. Governing equations

Classical mechanics in this chapter is narrowly defined as the investigation of the motion of systems of bodies under the influence of forces and torques. The problem is to determine the positions of all the particles at an arbitrary time tfrom their initial state at t= 0. Newton’s laws for the motion of bodies can be summarized as follows:

Newton’s first law states that an object will remain at rest or in linear motion unless acted upon by an external force. (This first law describes a constant velocity motion in the absence of an external force, which is a special case of zero acceleration in the second law below.)

Newton’s second law states that the net (unbalanced) force Fon an object is equal to the product of the mass mand acceleration aof the object:

F=maE1

where the acceleration is defined as the rate of change in velocity with respect to time. (This law indicates that the net force modifies the object’s velocity with respect to time. Then, the mass mcan be interpreted as the object’s resistance to the velocity change.)

Newton’s third law states that all the forces in the universe occur in equal (in magnitude) but opposite directions. For example, when one body exerts a force on a second body, the second body simultaneously exerts a force, on the first body, equal in magnitude and opposite in direction. These two forces cancel each other, and therefore, the net sum is always zero, even if particles are inelastic.

As noted above, particle size is neglected in describing its motion. The possibility of doing so depends on the actual size of the object and/or its distance from the observer. But, when a group of constrained particles forms a rigid body, its rotational motion is described using an equation similar to Newton’s second law of force, which is

M=IαE2

where Mis the torque or moment, Iis the mass moment of inertia, and αis the angular acceleration. A rigid body has six degrees of freedom: three for translation and the other three for rotation. One can combine Eqs. (1) and (2) to write the governing equation of motion of a rigid body in a simple matrix form:

[FM]=[M00I][aα]E3

where the zeros in the off-diagonal terms indicate that the medium in which particles are moving is not viscous, i.e., conceptually similar to vacuum. The linear and angular accelerations are defined as

a=r¨=d2rdt2andα=θ¨=d2θdt2E4

respectively, where rand θare the linear and angular positions of the object, of which time derivatives are the linear and angular velocities, respectively:

v=r˙=drdtandω=θ˙=dθdtE5

Where_{ }aand αhave three components each so that [a, α]^{T}in Eq. (3) is a vector of six components, where the superscript Tindicates a transpose. For mathematical simplicity, a generalized coordinate q, generalized velocity q˙, and generalized force Qof an object are written as

q=[rθ],q˙=[r˙θ˙],andQ=[FM]E6

The classical mechanics problems are usually reduced to solve for q(t) and q˙(t)at time tunder the influence of specified Q, using the initial conditions of q(t= 0) and q˙(t=0).

1.2.2. Total energy sum and difference

Fundamental questions from physicists are “What governs the motion of an object?” and “How can the motion be described and predicted mathematically?” Then, a fundamental question that naturally rises is “Is there a more fundamental principle that nature follows other than Newton’s second law?”

Let’s consider a particle, i.e., a point mass, found at position r(t) with velocity v(t) under the influence of force f(r), depending on the particle position only. We consider Newton’s second law in one-dimensional space:

f(x)=mdvdtE7

Because dx= vdt, we multiply dxwith f(x) and vdfwith mdv/ dtto derive

f(x)dx=mvdvE8

and integrate each side from state 1 of (x_{1}, v_{1}) to state 2 of (x_{2}, v_{2}) to obtain

where ΔW=∫x1x2f(x)dxis a work doneand T=12mv2is the kinetic energy. Eq (10) indicates that the work done is equal to the kinetic energy change from states 1 to 2. The integration of force with respect to xin Eq. (9) can be exact if the force depends on the particle position only. If so, this force is called conservativeand becomes the satisfactory condition for the energy conservation principle. A conservative f(x) can be expressed as a gradient of a scalar function V:

f(x)=−dV(x)dxin1‐DE11

f=−∇V(r)in3‐DE12

where Vis called the potential energy function. Then, the force integration is simply

where total energy Edefined as a sum of the potential and kinetic energies, i.e., E= T+ V, is derived as constant regardless of the path the particle takes. The potential energy difference depends on only the initial and final positions of x_{1} and x_{2}, respectively. Note that the total energy is conserved only if the force depends on only particle location. In advanced classical mechanics, the total energy Eis replaced by a Hamiltonian function: H= T+ V, and problems can be solved identically to applying Newton’s second law. Using the Legendre transformation of the Hamiltonian H, a new function called Lagrangian Lis defined as the difference between the kinetic and potential energies:

L=T−VE16

Instead of dealing with force vectors in Newton’s second law, Lagrangian mechanics uses the scalar Lagrangian function, which is assumed to contain all the information of the mechanical system.

1.2.3. Principle of the least action

The most general formulation of the law governing the motion of mechanical systems is the principle of least actionor Hamilton’s principle[1]. According to this principle, a mechanical system is characterized using a definite Lagrangian function L(q1,q2,…,qs,q˙1,q˙2,…,q˙s,t)or briefly L(q,q˙,t), and the motion of the system is such that a certain condition (discussed below) is satisfied.

At time t_{1} and t_{2}, particle positions are defined by two sets of the generalized coordinates, q(t_{1}) and q(t_{2}). The condition is that the system moves between these two positions, minimizing the integral

S=∫t1t2L(q,q˙,t)dtE17

to the least possible value. The integral of Eq. (17) is called the action. Note that the Lagrangian contains generalized coordinates and velocities, qand q˙only (not the higher derivatives such as q¨), as independent variables.

Let us now derive the differential equations that minimize the action integral of Eq. (17). For simplicity, the system is assumed to have only one degree of freedom. Let q= q(t) be the function for which the action Sis a minimum. This means that Schanges when q(t) is replaced by

q(t)+δq(t)E18

where δq(t), called a variation of q(t), is a function, assumed to be small everywhere in the interval of time from t_{1} to t_{2}. Since Eq. (18) must include the values of q(t_{1}) and q(t_{2}), we can now conclude that

δq(t1)=δq(t2)=0E19

In this case, the change in Swhen qis replaced by q+ δqis equal to

∫t1t2L(q+δq,q˙+δq˙,t)−∫t1t2L(q,q˙,t)E20

We expand the difference in powers of δqand δq˙in the integrand and leave only the first-order terms. Then, the principle of least action may be written in the form

δS=δ∫t1t2L(q,q˙,t)=0E21

or equivalently

∫t1t2(∂L∂qδq+∂L∂q˙δq˙)dt=0E22

Since δq˙=dδq/dt, we integrate the second term of Eq. (22) by parts to obtain

δS=[∂L∂q˙δq]t1t2+∫t1t2(∂L∂q−ddt∂L∂q˙)δqdtE23

The condition of Eq. (19) shows that the integrated term in Eq. (23) is zero:

[∂L∂q˙δq]t1t2=0E24

and the remaining integral is

∫t1t2(∂L∂q−ddt∂L∂q˙)δqdt=0,E25

which must vanish for all values of δq. This can be satisfied if and only if the integrand of Eq. (25) is zero. Thus, we have

ddt(∂L∂q˙)−∂L∂q=0E26

When the system has more than one degree of freedom, then Eq. (26) becomes

ddt(∂L∂q˙i)−∂L∂qi=0(i=1,2,…,s)E27

where sis the total degrees of freedom of the particles in the system. Eq. (27) is a set of required differential equations, called in mechanics Lagrange’s equations. If there is no constraint, the total degrees of freedom of a system containing N_{p}objects, s, is equal to 6N_{p}.

The zero on the right-hand side of Eq. (27) implies that the total energy is conserved because particle forces depend on their positions only. There are no forces/torques that dissipate the energy. If a number of particles are investigated under the influence of dissipative forces, such as hydrodynamic drag and frictional forces, Eq. (27) should be modified to

ddt(∂L∂q˙i)−∂L∂qi=Qi†(i=1,2,…,s).E28

where Q^{†} is the generalized non-conservative force. Therefore, we take Eq. (28) as the general governing equation of motion for multi-body granular dynamics in the regime of classical mechanics and microhydrodynamics [2]. It generalizes Newton’s second law and usually provides great mathematical simplicity by dealing with a scalar function Linstead of force vectors.

In science and engineering disciplines, dynamic simulations of particulate materials allow researchers to investigate microscopic many-body phenomena and further predict macroscopic material properties. In a liquid (aqueous) phase, rigorous and accurate simulations of particle dynamics must consider the repulsive volume-exclusion between polydispersed particles. Here, we briefly review the historical development of particle dynamics simulation methods in various scales.

2.1. Molecular dynamics

Conventional molecular dynamics (MD) treats particles (such as ions and molecules) as interacting point-masses and updates their present positions and velocities (at time t) to those in the future (at time t+ dt) [3]. Basic potential forms include the hard-sphere and Lennard-Jones potential [4]. A specific analysis of particle trajectories can provide measurable macroscopic quantities such as solute diffusivities in an aqueous medium and heat capacities in various thermal conditions. The choice of time internal dtmust be much smaller than the typical time for a molecule to travel a distance of the same order as its size. One of the most widely used methods to integrate the governing Eq. (1) is the Verlet algorithm [5], which provides a direct solution of the second-order ordinary differential equation with errors of δt^{4} order. The particle position at time t+ δtcan be expanded at time tusing Taylor’s series:

r(t+δt)=r(t)+v(t)δt+12a(t)δt2+O(δt3)E29

and by replacing +δtwith +δt, one obtains

r(t−δt)=r(t)−v(t)δt+12a(t)δt2−O(δt3)E30

The position rand velocity vat the future time t+ δtcan be calculated by adding and subtracting Eqs. (29) and (30) such as

r(t+δt)=2r(t)−r(t−δt)+δt2a(t)E31

and

v(t)=r(t+δt)−r(t−δt)2δtE32

respectively. As the primary objective of MD is to evolve the particle system, one of the advantages of the Verlet algorithm of Eq. (31) is that the velocity does not need to be calculated unless needed during simulations. This advantage does not exist if hydrodynamic drag forces are considered, the magnitude of which increases with respect to the relative velocity of particles to that of the fluid medium. The particle position at the next time step is calculated using those of the past and current time steps and the acceleration vector aobtained by the net force exerted on the particles. Velocity vectors can be additionally calculated during or after simulations if the kinetic energy needs to be calculated. Advanced integration algorithms over Verlet’s algorithm include the half-step leap-frog scheme [6] and velocity-Verlet algorithm [7]. Several leading MD simulation packages include assisted model building and energy refinement (AMBER) [8], chemistry at HARvard macromolecular mechanics (CHARMm) [9], Groningen machine for chemical simulations (GROMACS) [10], and nanoscale molecular dynamics (NAMD) [11]. The main difference between these MD simulation packages is how the force fields (functional forms and parameter values) are determined and used as parameters in potential functions. Applications of MD simulations are diverse and can be used in a wide variety of chemical and environmental applications. The pre-developed time evolution algorithms can be used in most simulations of particle dynamics, governed by Newton’s second law. On the other hand, a direct application of MD for granular dynamics simulation is limited to mimicking granular phenomena and properties because MD was originally developed for particles treated as point masses.

2.2. Brownian dynamics

Experimental and theoretical studies on Brownian motion were initiated by Brown [12], Einstein [13], Langevin [14], and Chandrasekhar [15, 16]. When solute motion in a solution is of greater interest, the motion of the tremendous number of ambient solvent molecules cannot be or does not need to be computed. Therefore, effects of solvent motion are replaced by random, fluctuating forces acting on solutes. Then, the governing equation is switched from Newton’s Eq. (1) to Langevin’s equation:

mdvdt=−ξv+f(r)+f′(t)E33

where – ξvis the hydrodynamic drag force acting in the opposite direction to the relative velocity of particle in the fluid medium, ξis the drag coefficient (to be expanded to a tensor form), f(r) is the conservative force derived from the potential energy function, and f′(t) is the random fluctuation force caused by adjacent, rapidly-moving solvent molecules. It is assumed that f′ is a Gaussian process with infinitely small correlation time, which renders the famous fluctuation-dissipation theorem [17]

⟨f′(t1)·f′(t2)⟩=2ξkBTδ(t1−t2)E34

⟨f′(t)⟩=0E35

After a sufficiently long time after the initial state which ensures Eq. (35) is true, one can approximate Eq. (33) in the absence of fas

⟨dvdt⟩≃−⟨vm/ξ⟩≃−⟨v⟩τE36

where ⟨⋯⟩ indicates averages over time. Eq. (36) allows us to determine the time scale of the decelerating motion, the so-called relaxation time, defined as τ= m/ ξ. The time step dtin Brownian dynamics (BD) should be much longer than the particle relaxation time, i.e., dt≫ τ. Only for a dilute solution, the drag coefficient ξcan be regarded as a constant. Stokes derived the drag coefficient of a spherical particle moving in a fluid medium as ξ= 3 πd_{p}μ, where d_{p}is the (hydrodynamic) radius of the particle and μis the absolute viscosity of the fluid medium [18]. In general, Stokes-Einstein diffusivity is defined as

DSE=kBTξ=kBT3πμdpE37

Although BD adopts the Oseen tensor, it still treats a particle as a point mass (like MD). From BD simulations, the hydrodynamic diameter d_{p}can be inversely determined by matching experimental data and simulation results. As a consequence, the volume-exclusion based on particle sizes is commonly disregarded in MD as well as BD simulations [19–23]. Specific features of MD and BD as applied to colloidal systems can be found elsewhere [24].

2.3. Dissipative particle dynamics

The restrictive condition of the time interval (dt≫ τ) is fundamentally resolved in dissipative particle dynamics (DPD) for finite-sized particles by incorporating the Fokker-Planck equation and the Ito-Wiener process [25–31]. The total force on particle ifrom other particles is written in the form of:

mai=∑j=1,j≠iNp(FijP+FijD+FijR)E38

where F^{P}is a conservative inter-particle force, F^{D}is a dissipative force, and F^{R}is a random fluctuation force. It is required that F^{D}and F^{R}are linear and independent of the momentum, respectively. A simple form of these force are hypothesized as

FijD=−γωD(rij)(eij⋅vij)vijE39

FijR=σωR(rij)eijζijE40

where r_{ij}= |r_{i}− r_{j}|, v_{ij}= v_{i}− v_{j}, e_{ij}= (r_{i}− r_{j})/r_{ij}, and σ=2kBTγ. Importantly, ζ_{ij}= ζ_{ji}is a Gaussian white-noise ther of zero mean and unit variance. The stochastic differential equation (SDE) of DPD consists of

dri=vidtE41

mdvi=∑j[FijP+FijD]dt+∑jσωR(rij)eijdWijE42

where dW_{ij}(= dW_{ji}) is independent increment and of the Wiener process, satisfying

dWijdWkl=(δikδjl+δjlδjk)dtE43

i.e., dW_{ij}(t) is an infinitesimal of order 1/2 (i.e., proportional to δt).

In Eq. (42), the momentum changes due to the conservative and dissipative forces are proportional to dt, and due to the random fluctuating, force is linear on dt. Hydrodynamic resistances are presumed to be pairwise and described as (intuitively chosen) simple functions of the inter-particle distance r:

ωD=ωR2={(1−r/rc)2forr<rc0otherwiseE44

where r_{c}is the cut-off distance. For two particles in close proximity, it was reported that the pairwise summation of hydrodynamic forces can be erroneous in estimating the many-body hydrodynamic forces/torques, represented using fourth-order tensors [32]. Interestingly, DPD was widely used to investigate the motion of macromolecules and polymers, which are soft and deformable, in liquid phases [25–31, 33–35]. The computational advantages of DPD include that, first, the tensor-wise hydrodynamic interactions are simplified to pairwise forces (although this approach has less fundamental rigor); second, the smooth functional form of Eq. (44) allows multiple soft particles to physically overlap each other, in which repulsive lubrication forces are disregarded; and third, the fluctuation-dissipation theorem is automatically satisfied [17] so that the thermodynamics of the system are well described.

2.4. Stokesian dynamics

The lubrication tensors in Stokesian dynamics (SD) may mimic the volume exclusion forces, which is logarithmically proportional to the surface-to-surface distance between two particles [36]. When two spheres are colliding or in contact with each other, the surface-to-surface distance converges to zero and the lubrication force diverges to infinity. In SD, a many-body, far-field, grand mobility matrix M∞is built using the pairwise superposition of the two-body mobility matrix. Its product with the hydrodynamic force F^{H}gives the relative velocities of particles to the fluid flow, ΔU[37–41]:

ΔU=[M∞]⋅FHE45

where the superscript ∞ indicates that the near-field lubrication forces are excluded. The grand resistance matrix is formed by inverting the grand mobility matrix and correcting the near-field lubrication such as

[R]=[M∞]−1−[ΔR2B]E46

where [ΔR2B]indicates a pairwise addition of the exact two-body resistance matrix Ri,j(2)subtracted by the inverse of two-body far-field mobility matrix, Mij(2):

[ΔR2B]=∑i,j(Ri,j(2)−[Mij(2)]−1)E47

SD provides the many-body diffusion tensor, which has a significantly higher accuracy than those of BD and DPD. SD is, however, fundamentally limited to rigid particles of no contact. Similar to BD, SD uses the Langevin equation of the drag force (Eq. (33)) and therefore is subjected to the intrinsic restriction of the time interval, dt≫ τ. Most of the SD simulations often deal with zero-inertia motion so that only driftingmotion of particles is investigated with no-acceleration under conservative and external force fields in a viscous fluid media. To generally mimic the complex hydrodynamic motion of many particles in a fluid medium, the zero-force assumption should be relaxed and the time step dtis arbitrarily chosen for computational efficiency under given physical conditions. Recent studies on SD include fast numerical inversion of the grand mobility matrix to accelerate the computational time, but the hydrodynamics is similarly mimicked [42, 43].

2.5. Dissipative hydrodynamics

Dissipative hydrodynamics (DHD) was recently developed to unify the above-mentioned particle dynamics methods [44]. DHD employed specific advantageous features from various simulation methods, specifically the many-body hydrodynamic tensors from SD and the stochastic differential equation from DPD. DHD can mimic the translational as well as rotational motion of N_{p}particles in a viscous fluid of temperature T. Importantly, it is free from the restriction of the particle relaxation time. The stochastic governing equations of DHD are represented as

M⋅du=[Qp−R⋅(u−U)]dt+B⋅dW=QdtE48

where Mis the mass/moment-of-inertia matrix of 6N_{p}× 6N_{p}dimension, uand Uare translational/rotational velocities of particles and the fluid, respectively, Q^{p}is the conservative force/torque vector, Ris the grand resistance matrix, Bis the Brownian matrix of zero mean and finite variance, i.e., ⟨B⟩ = 0 and ⟨B^{tr} ⋅ B⟩ = 2k_{B}T, where k_{B}is the Boltzmann constant, and dWis the Ito-Wiener process of 6N_{p}elements [45, 46]. The generalized force/torque is defined, similar to Eq. (42), as

Q=Qp−R⋅(v−U)+B⋅wdtE49

where w=dW/dt. A simple time integration using an intermediate time step is adopted from standard DPD algorithms [31]. It is worth noting that the SDE (Eq. (48)) relaxes the intrinsic time interval restriction of dt≫ τ, which BD and SD are subjected to. Therefore, DHD allows mimicking accelerated motion of particles of different sizes. Mathematical details of the DHD formalism can be found elsewhere [44]. Applications of DHD and its related work in environmental engineering include, but are not limited to, aggregate formation [47, 48] and single collector granular filtration [49]. DHD provides the most fundamental and accurate simulation algorithms for polydisperse particles from nanometer to millimeter sizes without arbitrarily turning on or off specific force and torque terms. But, DHD is computationally intensive to the same degree as SD so that high-performance parallel computation is inevitable to simulate a reasonably large numbers of particles.

2.6. Discrete element method

In the dynamic motion of granular particles, ballistic collisions are one of the most fundamental and important interactions. Suppose two (spherical) particles iand jof mass m_{i}and m_{j}, respectively, undergo an inelastic collision, as shown in Figure 1. The relative velocity of particles iand j, denoted as g_{ij}, at the point of contactis determined by the translational and rotational particle velocities:

gij=vij−(aiωi+ajωj)×nijE50

where v_{ij}= v_{i}− v_{j}is the relative velocity of the center of mass of particle ito jof velocity v_{i}and v_{j}, respectively. In general, particles are polydispersed, and a_{i}and a_{j}are the radii of particles iand j, respectively. The normal and tangential collision velocities given by the projections of g_{ij}are

gijn=(gij⋅nij)nijE51

and

gijt=−nij×(nij×gij)E52

respectively. The coefficients of restitution in normal and tangential direction, ϵ^{n}and ϵ^{t}, are defined as

(gijn)′=−ϵngijnE53

(gijt)′=+ϵtgijtE54

where 0 ≤ ϵ^{n}≤ 1 and −1 ≤ ϵ^{nt}≤ 1. The primed and unprimed variables indicate the pre- and post-collision quantities, respectively. Then, the velocities of particles iand jafter the collision can be represented as functions of ϵ^{n}, ϵ^{t}, gijn, gijt, m_{i}, and m_{j}. Now, for particles iand j, we finally have

v′i=vi−(μijmi)[(1+ϵn)gijn+1−ϵt1+1/J˜gijt]E55

v′j=vj+(μijmj)[(1+ϵn)gijn+1−ϵt1+1/J˜gijt]E56

ω′i=ωi+(μijmiai)(1−ϵt1+J˜)nij×gijtE57

ω′j=ωj+(μijmjaj)(1−ϵt1+J˜)nij×gijtE58

where for k= i, j. Here, J_{k}is the mass moment of inertia and J˜k=Jk/mkak2is its dimensionless form. A spherical particle has J_{k}= 2/5. The dimensionless mass moment of inertia for (rigid and spherical) particle is assumed to be constant. During the collision, the kinetic energy of particles is lost and transferred into thermal energy of the ambient fluid. At an arbitrary time before or after the collision, the kinetic energy of particle kis represented as

Tk(′)=12mkvk(')⋅vk(')+12Jkωk(')⋅ωk(')E59

One can calculate the energy loss of particle k, ΔTk=T′k−Tk, using the pre- and post-collision velocities of (vk,ωk)and (v′k,ω′k), respectively. Post-collision velocities of two unequal spheres are completely solved using Eqs. (55)–(58), but an underlying assumption is that the inelastic collision is instantaneous without spending any time. In reality, however, any granular collision takes a finite amount of time, even if it is much shorter than the traveling time of the particles. This time duration in which two particles are in contact is defined as contact time(or collision duration). Therefore, it is more accurate to see the collision of two granules as an impulse event. The collision rule well determines the post-collision states, if the granules are in a fluid-like state. But, when densely-packed granules are slowly moving under mechanical or gravitational compression, the collision rule fails to predict the transient granular states because it does not take into account the compression and restoration of the inelastic particles.

2.7. Overview

Table 1 summarizes specific features of the above-discussed simulation methods. Acceleration can be included by all methods in principle, but is barely employed in BD and SD. This is mainly due to the time interval restriction. Instead, the fluctuation-dissipation theorem is intrinsically embedded and satisfied in BD, SD, DPD, and DHD formalisms. All the simulation methods can surely include effects of conservative, external force fields. Brownian motion and multi-body hydrodynamic forces are (most) accurately implemented in SD and DHD, but only hydrodynamics is important in the many-body motion of non-Brownian granules. The DHD is the only simulation method that can mimic many-body hydrodynamics without the time interval restriction. Constraint forces/torques can be easily applied to any dynamics methods to simulate compound particles or aggregates. SD has the same capability but the collision rules are not included. Due to typical ranges of particle sizes in the simulation methods, the inelastic collision (or contact) forces are included only in DEM. Since DPD allows particle deformation by using the relaxed, pairwise hydrodynamic resistance between two particles, i.e., ω_{D}and ω_{R}in Eq. (44), employing the collision force of rigid bodies in DPD must be out of its original scope. The current state-of-the-art DEM algorithms can be further improved to mimic complex phenomena of dry granules of various shapes [50, 51]. However, in our opinion, granular motion in a liquid phase can be accurately simulated by only using DHD or SD. To accurately simulate a large number of granules of various sizes in complex fluid environments, multi-body hydrodynamics must be rigorously implemented in the current DEM method.

Comparison of particle dynamics simulation methods. The conservative force includes external forces such as electromagnetic and gravitational, and inter-particle force such as DLVO [52, 53]. Here, “any” dtindicates that dtcan be arbitrarily determined for accurate and fast calculation regardless of the particle relaxation time τ.

Granular materials consist of a large number of particles whose typical size ranges from micrometers to centimeters [50, 51]. These particles interact via short-range forces through only mechanical contacts and (external) long-ranged electromagnetic or gravitational forces. Granular dynamics mimics dynamic motion of granular particles in a transient state such as excited or granules in a fluid media. Large-scale phenomena in this category include land sliding and snow avalanches. Wet granules such as sand in beaches undergo hydrodynamic forces due to tidal currents. Soil granules in unsaturated sub-surfaces are partially dry or wet. Interstitial water layers between granules can significantly change inter-granular interactions, especially when they are in a stationary contact phase. The significance of hydrodynamic interactions between non-Brownian granules is paid less attention. In this section, we describe granular dynamics as a microscopic extension of DEM in a shorter time scale by investigating the collision phenomenon between two (spherical) particles, as shown in Figure 1, during an impulse event.

3.1. Individual contact forces

3.1.1. Normal and tangential forces

Small particles such as suspended solid, colloids, and nanoparticles are naturally negatively charged. Their electrostatic repulsive forces decay exponentially with respect to the distance from their surfaces to the bulk phase. When the surface-to-surface distance between two particles is much smaller than their average sizes, forces such as electrostatic and Born repulsion are strong enough to repel each other. These forces are, however, not dominant for large, non-Brownian particles such as granules of an order of O(0.1 − 10) mm. The dominant granular force stems from contact during collisions. In conventional statistical mechanics, a hard sphere is characterized by the wall potential:

V={0r>2a∞r<2aE60

where ais the particle radius. The mathematical discontinuity of the wall potential at r= 2aindicates the infinite force that completely prevents any overlap between particles.

As granules are inelastic, the fundamental wall potential must be modified before it is applied to a granular system. Two spherical granules are in a mechanical contact if the sum of their radii exceeds their center-to-center distance, i.e.,

ξij≡ai+aj−|ri−rj|>0E61

where ξ_{ij}is the mutual compressionof particles iand j. Note that ξ_{ij}is positive when two granules overlap and becomes larger when the granules come closer. Thus, the mutual compression can be interpreted as the overlapped surface-to-surface distance. The force acting on particle ifrom particle j(conceptually denoted as i← j) is described by

Fij={Fijn+Fijtforξij>00forξij≤0E62

where Fijnand Fijtare the normal and tangential components of the contact force, respectively. For simplicity, Eq. (62) can be rewritten as

Fij=(Fijn+Fijt)H(ξij)E63

where H(x) is the Heaviside step function, defined as

H(x)={1forx>00othersiesE64

In three-dimensional space, vector quantities between particles iand jcan be decomposed into the normal and tangential directions. Using the mathematical identity

A×(B×C)=B(A⋅C)−C(A⋅B)E65

one can replace Aand Bby nand Cby F_{ij}to write

n×(n×Fij)=n(n⋅Fij)−FijE66

Fij=n(n⋅Fij)−n×(n×Fij)E67

From Eq. (67), the normal and tangential force components can be expressed as

Fijn=(n⋅Fij)nE68

Fijt=(n×Fij)×nE69

The contact force between elastic spheres was originally developed by Hertz [55] and was later generalized for viscoelastic (damped) particles [56, 57] as

Fn=2Yaeff3(1−ν2)(ξ3/2+Aξdξdt)E70

where Yand νare Young’s modulus and Poisson’s ratio, respectively, Ais the dissipative constant being a function of material viscosity, and a_{eff} is the effective radius that can be interpreted as a harmonic sum a_{i}and a_{j}, i.e., aeff−1=ai−1+aj−1. Parameter Aexplains the dependence of the restitution coefficient on the approaching velocity between two spheres. If A= 0, then Eq. (70) converges to the original Hertz’s equation for elastic granules. Therefore, parameter Aneeds to be inversely calculated using an experimentally observed coefficient of restitution. If two elastic particles are heterogeneous, then Hertz’s equation may be extended to

Fijn=4aeff3(1−νi2Yi+1−νj2Yj)−1ξ3/2E71

for particles of different Yand νvalues. Assuming the dissipative constant Ais also particle-specific and additive, the most general expression of the normal force between two viscoelastic spheres may be [51]

Fijn=4aeff3(1−νi2Yi+1−νj2Yj)−1(ξ3/2+Ai+Aj2ξ˙ξ)E72

as a generalized extension from Eq. (70). Note that this viscoelastic force is finite while two spheres are being overlapped so that the mutual compression ξis positive. Let’s define t= 0 as the moment of the contact of two particles. The compression continues until t= t_{c}, after which restoration begins. When the relative velocities between the two particles at t= 0 and t= t_{c}are ξ˙(0)and ξ˙(tc), respectively, then the normal restitution coefficient can be calculated by measuring velocities ξ˙(0)and ξ˙(tc), i.e.,

ϵn=ξ˙(tc)/ξ˙(0)E73

This indicates that ϵ^{n}depends on the relative collision velocity, unless A= 0. Theoretically, at least three experiments are required for collisions between particles iand jfor pairs of (i, j), (i, i), and (j, j). Then, A_{i}and A_{j}can be inversely calculated by using the trial-and-error method for numerical fitting.

The relative velocity of the spheres at the point of contact results from the relative translational/rotational velocities. The contact-point velocity has the tangential component of

vrelt=(vj−vi)⋅eijt+aiωi+ajωjE74

which provides the tangential force of

Ft=−sign(vrelt)⋅min(γt|vrelt|,μ|Fn|)E75

where γ^{t}is a fitting coefficient, proportional to the tangential dissipative force in magnitude. In Eq. (75), the shear force is limited by Coulomb’s friction law of |F^{t}| ≤ μ|F^{n}|, where μis the friction coefficient. Although this approach is conceptually straightforward, the level of approximation is still on the collision rule, as discussed in Section 2.6. Unless F^{n}is constant, Eq. (75) provides an inconsistent value of F^{t}because F^{n}(if Eq. (72) is used) is a function of not only the mutual compression, but also the relative velocity. The tangential contact force is correlated to the normal force [56, 57]:

Fijt=−μFijn(ζijζ0−⌊ζijζ0⌋)E76

where ζis the relative tangential shift, ζ_{0} is its macroscopic maximum value, and ⌊x⌋ denotes the integer of x. Typical values of ζ_{0} / a_{eff} range from 10^{−7} to 10^{−3}. Similar to Eq. (73), the tangential coefficient of restitution can be calculated as

ϵt=ζ˙(tc)/ζ˙(0)E77

If the friction coefficient μis known, experimental measurement of ϵ^{t}allows us to inversely calculate ζ_{0}. Techniques to determine the coefficients of restitution of colliding viscoelastic spheres can be found elsewhere [58]. This approach can be readily applied to densely packed granular medium with small movements or vibrations such as the Brazilian bean problem [59–64]. However, time integration per collision, consisting of compression and restoration, should be from 0 to t_{c}, while t_{c}is a negligibly short time in the collision-rule approach. For efficient computations, regular time-integrations and event-based simulations should be efficiently combined in order to have variable δt, which should be much smaller than t_{c}in compressing/restoring phases, and much longer than t_{c}for collision events for fluidized granules of a low concentration.

3.1.2. Shear stress

Shear thickening:When non-Brownian granules are densely packed, volume fraction is around 50%, depending on their polydispersity and shape. The granules form a loosely connected material, which responds to the external shear stresses in an unconventional way. Depending on the magnitude of the external stress, granules temporarily switch their phases between the liquid-like and the solid-like states. In a (pure) fluid, viscosity is defined as the ratio of shear stress to the shear rate during a steady flow, which represents the energy dissipation rate by the fluid flow. This dissipation rate in some cases decreases with respect to the shear rate, which is known as shear thinning. It is particularly desirable for paints such that pigments flow easily when brushed, but does not drip when brushing stops. On the other hand, shear thickeningindicates that the energy dissipation rate increases as the shear rate increases. In other words, the fluid becomes much more viscous if the external stress is strong enough. For example, if a large amount of cornstarch or Oobleck is mixed with water in a (small) swimming pool [65, 66],^{1} the dense suspension acts like a liquid if it is at rest in the absence of external stresses applied. But, when the suspension is sheared, the flow resistance increases dramatically and the fluid becomes locally amorphous for a short amount of time. If a person is continuously stepping on the dense suspension, the person will be able to dynamically stay on top of the semi-liquid surface. When the person slows down or stops moving, then he or she will slowly sink into the pool.^{2} This phenomenon is not only interesting, but also has practical importance to engineering systems like automobile brakes [67].

Possible mechanisms of this shear-thickening include hydroclustering, order-disorder transition, and dilatancy. First, in hydroclustering, particles gather together into transient, reversible clusters under shear flow, and this rearrangement leads to increases in the lubrication drag forces [68, 69]. Local heterogeneity of particle suspension creates regions in which particles undergo less drag forces and tend to agglomerate. This results in narrow flow channels among the dynamic groups of particles. Application of Stokesian dynamics to mimic the hydroclustering intrinsically prevents particle contacts, followed by inelastic overlaps in a series of collision events. This is because the lubrication force is logarithmically proportional to the surface-to-surface distance between two particles and thus it diverges when two particles start inelastically overlapping. Second, the order-disorder phase transition includes the changes in the flow structure between ordered and disordered phases, which yields increases in the drag force between particles [70]. These ordered states are different from fixed 3D structures found in solids, but similar to amorphous aggregates. Third, the dilatancy mechanism describes the shear-thickening such that the effective volume of particle packing increases under the shear. When particles are confined in local spaces or partially jammed, the shear force pushes particles toward the containing walls. Additional stress can be developed on the walls and backward responses may generate extra stress between particles and walls [71] in contact interfaces. Recently, discontinuous shear thickening (DST) was proposed to explain the shear thickening experiments [69], of which comprehensive review can be found elsewhere [72, 73]. To the best of our knowledge, the fundamental mechanism of the shear thickening has not been fully discovered.

Sample simulation:Figure 2 shows results of a sample simulation using the mechanisms described in Section 3.1.1, which can be further extended to DST simulations. The gray spheres are loosely packed, ideally forming a body-centered cubic structure. The top layer consists of 5 × 5 = 25 regularly packed spheres, under which 6 × 6 = 36 spheres are located. These spheres are closely located to each other but not touching, having very small surface-to-surface distance between the nearest neighbors. The two layers are packed five times vertically, and therefore the total number of gray particles is 366. This preliminary simulation aims to see the impact responses of the small packed spheres when hit by a big, fast intruder. The intruder particle initially moves with the downward velocity of v_{z}= 13.00 × 10^{−3} m/s_{ }and spins with angular velocity of ω_{z}= 0.56 rad/s.

Due to the loose packing of the small spheres on the bottom of the container, the intruder tends to penetrate the packed granules. In the initial stage of the penetration, the intruder is surrounded by small spheres, which bounce away from their initial positions. If the small spheres are initially touching each other, the force propagation from the top layer to the bottom layer must be almost instantaneous. The intruder will experience a strong normal force and may rebound in the opposite direction after an initial penetration of a short depth. As the collective phenomena of many granules are mostly transient, not only the material properties of the intruder but also initial and boundary conditions of granules play significant roles in their macroscopic dynamic behaviors. High dissipation rate of Areduces the impact of the one-to-many collisions between the intruder and dissipating granules. The intruder’s impact is transmitted through dynamic chains of contacting spheres, but not all of the packed granules participate in the force transmission. There are same granules almost free from the intruder’s impact.

Calculation of forces and torques exerted on each particle and their visualization can help understand the transient behavior and design dynamic granular materials. Granules initially located near the container walls have much less spaces to move. The kinetic energy of the intruder is transmitted to these granules at the cul-de-sac and mostly dissipated on the wall surfaces. Returning force to the intruder is initiated from the boundary. Applications of this simulation method can be designed to include a large number of real applications for characterization and prediction of transient granular material properties. One important key issue is, as indicated above, to identify and visualize the transmission chains of force/torque, which dynamically form and disappear. Figure 3 shows the initial impact event when the intruder starts penetrating the loose granular packing, visualized using open visualization tool (OVITO) [77]. The top and bottom rows show the top and side views of the intruder collision. The left column shows particle positions as well as force vectors, and the right column shows only particle configurations where color indicates the magnitude of the net force. Impact from the intruder is somehow irregularly distributed around the top granules. On the left-column images, arrows and their colors indicate force directions and magnitudes. A number of downward vectors indicate that the gravitational force is dominant for non-touching granules. The upward arrow in the intruder implies that the contact force to the intruder, which is similar to the normal force developed on the interface between an object and a wall, exceeds the gravitational force exerted on the intruder. Even though the force direction is temporarily inverted, the intruder still goes down due to the inertia of the high initial velocity. Figure 3 clearly shows that only a partial number of granules participate in the force transmission from the intruder to the packed granules. These force chains are very transient, and more importantly, the magnitude of the transmitted force diminishes as time goes by due to the intrinsic inelastic nature of the granules. For the future, designs of smart damping materials can be achieved not only by understanding various mechanical properties of the granules, but also by controlling specific initial and boundary condition, which leads inelastic many-granule systems to behave as smart materials.

3.2. Constraint force: holonomic and non-holonomic

3.2.1. Holonomic potential for translational constraints

Irregular-shaped granules can be mimicked as a large collection of polydispersed spherical particles. For simplicity, we first consider two particles (point masses) moving together as one body. A constraint embedded between the two particles keeps the inter-particle distance invariant. A translational constraint between particle iand jis

σij=(ri−rj)2−dij2=0E78

where d_{ij}is the fixed distance between particles iand j, usually an average of two contacting rigid bodies of diameters d_{i}and d_{j}, i.e., d_{ij}= (d_{i}+ d_{j})/2. This type of geometrical constraint is called holonomic. To develop an inter-particle interaction to satisfy Eq. (78), one defines the constraint potential for all N_{p}particles as

Φ=12∑i=1Np∑j>iNpλjiσijE79

where λ_{ji}is a symmetric Lagrange’s multiplier and 12in the front of the summation is by convention. Exchanging positions of particles iand j(i↔ j) should not change the sign and the magnitude of Φ. To satisfy this condition, the Lagrange multiplier is symmetric, i.e., λ_{ij}= λ_{ji}. The constraint force exerted on particle jcan be derived as a negative derivative of the constraint potential:

FjC=−∇jΦ=∑i≠jλji(ri−rj)E80

where iruns from 1 to N_{p}except for the case i= j(if so, λ_{ii}= 0), or simply

FjC=∑i=1Npλji(ri−rj)(1−δij)E81

where δ_{ij}is the Kronecker delta symbol, defined as

δij={1fori=j0otherwiseE82

Among N_{p}− 1 pairs made by particle j(excluding itself), if particle jdoes not have any constraint to particle k, then λ_{jk}= 0, and symmetrically vice versa. For a two-body case, the holonomic constraint force acting on jby kis

Fj←kC=−∇jΦ=λjk(rj−rk)E83

and, similarly, the same force on kby jis

Fk←jC=−∇kΦ=λkj(rk−rj)E84

Since λ_{jk}is symmetric (λ_{jk}= λ_{kj}), one can show that

Fj←kC=−Fk←jCE85

which follows Newton’s third law, the action and reactionprinciple. Summation over all constraint particles will give a zero resultant force. This is because the constraint force is an internal force and the sum of internal forces is zero.

Evolution of positions:Translational acceleration of particle jcan be expressed as

aj=aj†+ajCE86

where aj†is the unconstrained acceleration and ajCis the constrained acceleration, represented as

ajC=1mjλjk(rk−rj)E87

Assume particle jevolves from r_{j}(t) at time tto r(t+ δt) at time t+ δt. Then, the future position at t+ δtcan be decomposed into two parts:

rj(t+δt)=rj†(t+δt)+12ajCδt2E88

where

rj†(t+δt)=rj(t)+vj(t)δt+12aj†δt2E89

is the evolved position at time t+ δtin the absence of the constraint force. A similar equation for particle kcan be written easily. Note that Eq. (78) should be valid at all times. Substitution of Eq. (88) into (78) gives, neglecting terms on the order of (δt^{4}) and higher, the representation of the Lagrange multiplier:

Calculate the constraint acceleration of particle j, ajC, using Eq. (87).

Update the particle position using Eq. (88), which is under the influence of constrained and unconstrained accelerations. Set this updated position as the unconstrained position at that time: rj†(t+δt)←rj(t+δt).

Having the updated rj†, go to step 2, unless λ_{ij}converges to a finite value. Otherwise, store information at t+ δt, and go to the next time step.

This iterative procedure will continue until the Lagrange multiplier converges within a preset tolerable error for the position evolution from time tto t+ δt. So far, the constraint force modifies the particle position at time t+ δtfrom its unconstrained position rj†(t+δt), but the velocity after the constrained evolution of position is the same as before the evolution.

Velocity evolution:Differentiation of the holonomic constraint Eq. (78) with respect to time gives

(ri−rj)⋅(vi−vj)=0E94

valid both at time tand t+ δt. Similar to the position evolution, the translational velocity at time t+ δtis represented as

vj(t′)=vj†(t′)+ajCδtE95

where

vj†(t′)=vj(t)+aj†δtE96

is the updated velocity from time twithout the holonomic constraint. Substitution of Eq. (95) into Eq. (94) gives

Δrij(t+δt)⋅Δvij†=−Δrij(t+δt)⋅ΔaijCδtE97

where

ΔaijC=κijμij(rj−ri)E98

is the constraint acceleration for velocity correction. Here, κ_{ij}plays a similar role of λ_{ij}, but it is independently determined only to update the velocity, which is calculated as

κijδt=Δrij(t+δt)⋅Δvij†dij2μij−1E99

To update κ_{ij}, a similar iteration method can be used:

Calculate the unconstrained velocity at the next time step, v^{†} (t+ δt) for particles iand j, and calculate their difference Δvij†=vi†−vj†at time t+ δt.

Calculate κ_{ij}using Eq. (99) using previously determined r_{i}(t+ δt) and r_{j}(t+ δt).

Update ΔaijCusing Eq. (98) and use it to calculate v_{j}(t+ δt) of Eq. (95).

Replace v_{j}(t+ δt) by vj†(t+δt)and go to step 2 unless κ_{ij}converges to a constant value within a tolerable error. Otherwise, store information at t+ δt, and go to the next time step.

So far, we first defined the constraint potential using unknown Lagrange’s multipliers, derived the constraint force and acceleration, and updated iteratively particle positions and velocities until all the constraints are satisfied. The velocity evolution under this constraint is similarly done by introducing a new independent Lagrange’s multiplier, which is to satisfy the orthogonal relationship between position and velocity variations, in Eq. (94).

Figure 4 shows a simple test of the holonomic constraint between two unequally-sized spheres. It is clear that the two spheres are in contact with each other during their translational motion. A camera is moving with the same velocity of their center of mass so that only the relative motion is shown. Both the particles are non-Brownian, and the red sphere is 1.5 times bigger than the blue one, while their specific gravity is 2.75. The blue sphere rotates in the clock-wise direction around the red sphere. This is because the red one is 1.5^{3} = 3.375 times heavier so that the total center of mass is closer to that of the red particle. Careful observation indicates that the blue sphere rotates in the clock-wise direction about its center of mass, and the red sphere rotates in the opposite direction about its center of mass. This relative rotation stems from the presence of only holonomic (translational) constraints, which allows their smooth surface-elements to slide relative to each other. This phenomenon must happen if weakly attractive particles form loose aggregates in a fast viscous flow.

3.2.2. Non-holonomic torque for angular constraints

Dynamics simulations with the holonomic constraints work perfectly within tolerable errors, especially when the particles sizes are smaller than center-to-center distances. This method was successfully used for molecules of fixed structures such as water (H_{2}O) and organic compounds. On the other hand, if two spherical particles (such as colloids) of finite volumes are attached by sticky surface forces, the translational and rotational motion of these two spheres are constrained as they move as a single compound body. For simplicity, we will consider a pair of compounded golf balls. If only the two constraints discussed in Section 3.2.1 are considered regarding the motion of the compounded golf balls, their rotations about their own centers of mass are still allowed even if surface friction exists. Mathematically, the two balls, if perfectly glued on their small shared surfaces, should have the same angular velocity. This type of constraint is based on (angular) velocity so it is called non-holonomic.

Consider a perfectly inelastic collision event between two approaching particles iand j, moving with translational and rotational velocities of (v_{i}, ω_{i}) and (v_{j}, ω_{j}), respectively, before their collision. For k= i, j, the linear and angular momenta are

pk=mkvkE100

Hk=JkωkE101

respectively. After the two particles are permanently attached, the total linear momentum is

MaVa=m1v1+m2v2E102

where M_{a}= m_{1} + m_{2} is the total mass of the two-body aggregate and V_{a}is the translational velocity of the center of mass of the aggregate. The total angular momentum has, however, a slightly different formulation:

Ha=(J1+m2d12)ω1+(J2+m2d22)ω2=JaΩaE103

where

Ja=I1+m1d12+I2+m2d22E104

is the mass moment inertia of the aggregate. Here, d_{k}is the shortest distance between the center of mass of particle kand the rotation axis of the aggregate, passing through the center of mass of the aggregate, r_{CM}:

dk=(rk−rCM)⋅Ωa|Ωa|E105

After the perfectly inelastic collision, two compounded particles will have the same angular velocity Ω_{a}. The total torque acting on the aggregate is rigorously represented as

M=∑j(rj−RCM)×(FjEx+FjC)E106

and the time evolution equations related to the total angular momentum are

dHadt=ME107

Ha(t+δt)=Ha(t)+MδtE108

Ωa(t+δt)=Ωa(t)+Ja−1MδtE109

Finally, all associated particles to the aggregate have the same angular velocity Ω_{a}at any time unless they break apart of slide into each other. It is assumed that the mutual compressions ξbetween the associated particles are small enough to be neglected in calculating the aggregate’s mass moment of inertia, J_{a}.

Figure 5 compares dynamics of a trimmer, i.e., three constrained bodies. The six small black particles per sphere are imaginary, showing how the particle rotates in its transient motion. These imaginary markers do not generate any forces or torques. Initially, three spheres associated with a trimmer make the ideal L-shape. The downward gravitational force causes the settling of the trimmer. As noted above, the camera is moving with exactly the same velocity of the trimmer’s center of mass. Therefore, one can observe only relative motion of three identical spheres with respect to their center of mass. The gap between the two closest spheres is equal to the diameter of the black marker. In Figure 5(a), three particles undergo only holonomic constraints such that the center-to-center distances are kept constant. As the outside surfaces of each particle experience higher hydrodynamic stress, all three particles try to rotate toward the center, as viewed from the top of the trimmer. On the other hand, Figure 5(b) shows a trimmer of three rigidly attached (glued) spheres. Relative rotation of a sphere with respect to its neighbors is prevented. This indicates that all the three spheres in Figure 5(b) have identical angular velocity, which is equal to that of the whole trimmer as a compounded rigid body. If a member of an aggregate can freely rotate, then the hydrodynamic stresses exerted on the outer surface of the compound body must be relaxed, allowing rotation of individual spheres, and therefore the net shear stress is reduced. If the same shear force is applied to the rigid trimmer, then the non-holonomic constraint strongly resists the external hydrodynamic stress and adjusts its position to minimize the external stress. The angles made by connecting three particle centers in the last snapshots in Figure 5(a) and (b) indicate the different responses of the settling trimmer to the hydrodynamic drags and stresses. In addition to these hydrodynamic forces and torques, inelastic properties of granules significantly influence their transient rotating patterns.

3.3. Parallel algorithms

Granular dynamics simulations in a fluid medium must be an open problem in state-of-the-art computational research. Since the length scale of the forces acting on touching granules is much smaller than the granular sizes, simulations seem to be very efficient if parallel algorithms are adequately used. Domain decomposition scheme is one of the most widely used parallel algorithms. The system is divided into several small sub-domains, and a dynamic simulation in a spatial sub-domain is conducted by an individual computing unit, such as a core. This decomposition method is very efficient if particles are evenly distributed in space, which is usual in equilibrium simulations of MD and Monte Carlo. As granules are highly subjected to the gravitational force and hydrodynamic interactions, spatial biases are almost inevitable in both their transient and stationary states. Granular dynamics has at least three length scales of different orders of magnitude. The mutual compression distance, i.e., inter-particle overlap distance, is at least three or four orders of magnitude smaller than the particle size. In a fluid medium, hydrodynamic interactions are long-ranged and quite significant when the surface-to-surface distance between nearest neighbors is about the particle diameter. Motion of heavy granules in a fluid flow may distort the ambient flow-field as well as the hydrodynamic forces exerted on adjacent particles. In granular dynamics, when the contact force exerted between two colliding particles, the force acting on particle ifrom jhas the same magnitude and opposite direction to that from particle jfrom particle i. Therefore, the number of contact force calculations can be reduced to a half, if Newton’s third law is implemented during parallel computing. In this case, the computational efficiency of granular dynamics simulation can be as much as doubled.

In nature, phases of matter have been conventionally believed to be those of gas, liquid, and solid, in which specific phase transformations are possible between the states. A recent addition of the plasma state has increased the total number of material states from three to four. It is questionable to predict that the granular state will be the fifth matter phase.

On the other hand, the granules dynamically change the representative phase based on external influences. Stationary and compressed granules behave similar to amorphous solids or gel materials. Moving like a fluid, mud or sediments create their own pathways by minimizing hydrodynamic influences. A fast flow with granular materials, such as in streams and ocean, creates a dense turbidity flow, but a decelerating flow field initiates a granular phase change from a flowing liquid to a packed solid. Small dry granules behave similar to dust in the wind, for which standard gas transport theory can predict dynamic behaviors of granular gases. In our opinion, granules are chameleon materials, transitioning their phases dynamically. No equilibrium exists in a granular phase so that thermodynamic fluctuation of the granular state cannot provide material or phase constants. As granular dynamics in transient force/torque fields significantly removes steady-state behaviors, the initial and boundary conditions become extraordinarily important in analyzing many-body granular motion. To utilize specific behavior of granules, granules can be dynamically controlled in vibrating, oscillating, or swinging phases. The time-correlation scale of dissipating granules is not short enough to use the Markovian chain concept because granular paths in both the real and phase spaces significantly influence their fates.

This work was financially supported by the National R&D project of “Development of 1MW Ocean Thermal Energy Conversion Plant for Demonstration” (PMS3680) from the Korean Ministry of Oceans and Fisheries and used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575.

References

1.Landau LD, Lifshitz EM. Mechanics. New York: Pergamon Press; 1976

2.Kim S, Karrila SJ. Microhydrodynamics, Principles and Selected Applications. Mineola, New York: Dover Publication, Inc.; 1991

3.Allen MP, Tildesley DJ. Computer Simulation of Liquids. New York, NY: Clarendon Press; 1987

4.Jones JE. On the Determination of molecular fields. II. From the equation of state of a gas. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. 1924;106(738):463–477

5.Verlet L. Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Physical Review. 1967;159(1):98–103

6.Hockney RW. The Potential Calculation and Some Applications. Methods in Computational Physics. 9. New York: Academic Press; 1970

7.Swope WC, Andersen HC, Berens PH, Wilson KR. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. The Journal of Chemical Physics. 1982;76(1):637–649

8.Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Ferguson DM, et al. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. Journal of the American Chemical Society. 1995;117(19):5179–5197

9.Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. Journal of Computational Chemistry. 1983;4(2):187–217

10.Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJC. GROMACS: Fast, flexible, and free. Journal of Computational Chemistry. 2005;26(16):1701–1718

11.Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, et al. Scalable molecular dynamics with NAMD. Journal of Computational Chemistry. 2005;26(16):1781–1802

12.Brown R. A Brief Account of Microscopical Observations made in the Months of June, July, and August, 1827, on the Particles contained in the Pollen of Plants; and on the General Existence of Active Molecules in Organic and Inorganic Bodies

13.Einstein A. Zur Theorie der Brownschen Bewegung. Annalen der Physik. 1906;324(2):371–381

14.Langevin P. Sur la théorie du mouvement brownien. Comptes-rendus de I'Académie des Sciences (Paris). 1908;146:530–533

15.Chandrasekhar S. Stochastic problems in physics and astronomy. Reviews of Modern Physics. 1943;15(1):1–89

16.Chandrasekhar S. Brownian motion, dynamical friction, and stellar dynamics. Reviews of Modern Physics. 1949;21(3):383–388

17.Kubo R. The fluctuation-dissipation theorem. Reports on Progress in Physics. 1966;29:255–284

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

19.Ermakova L, Sidorova M, Savina I, Aleksandrov D. Structural and electrochemical parameters of asymmetric membranes for reverse osmosis. Colloids and Surfaces A. 1998;142(2–3):265–274

20.Ermak DL, Buckholz H. Numerical integration of the Langevin equation: Monte Carlo simulation. Journal of Computational Physics. 1980;35:169–182

21.Ermak DL, Mccammon JA. Brownian dynamics with hydrodynamic interactions. Journal of Chemical Physics. 1978;69(4):1352–1360

22.Ermak D. A computer simulation of charged particles in solution. I. Technique and equilibrium properties. Journal of Chemical Physics. 1975;62(10):4189–4196

23.Ermak DL. A computer simulation of charged particles in solution. I. Technique and equilibrium properties. Journal of Chemical Physics. 1975;62(10):4189–4196

24.Chen JC, Elimelech M, Kim AS. Monte Carlo simulation of colloidal membrane filtration: Model development with application to characterization of colloid phase transition. Journal of Membrane Science. 2005;255(1):291–305

25.Hoogerbrugge P, Koelman J. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhysics Letters. 1992;19:155

26.Koelman J, Hoogerbrugge P. Dynamic simulations of hard-sphere suspensions under steady shear. Europhysics Letters. 1993;21:363

27.Klingenberg DJ. Simulation of the dynamic oscillatory response of electrorheological suspensions: Demonstration of a relaxation mechanism. Journal of Rheology. 1993;37(2):199–214

28.Espanol P, Warren P. Statistical mechanics of dissipative particle dynamics. Europhysics Letters. 1995;30(4):191–196

29.Groot RD, Warren PB. Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. Journal of Chemical Physics. 1997;107(11):4423–4435

30.Martys NS, Mountain RD. Velocity Verlet algorithm for dissipative-particle-dynamics-based models of suspensions. Physical Review E. 1999;59(3):3733LP-3736 LP

31.Nikunen P, Karttunen M, Vattulainen I. How would you integrate the equations of motion in dissipative particle dynamics simulations? Computer Physics Communications. 2003;153(3):407–423

32.Satoh A, Chantrell RW, Coverdale GN. Brownian dynamics simulations of ferromagnetic colloidal dispersions in a simple shear flow. Journal of Colloid and Interface Science. 1999;209:44–59

33.Ramezani M, Shamsara J. Application of DPD in the design of polymeric nano-micelles as drug carriers. Journal of Molecular Graphics and Modelling. 2016;66:1–8

34.Karatrantos A, Clarke N, Composto RJ, Winey KI. Topological entanglement length in polymer melts and nanocomposites by a DPD polymer model. Soft Matter. 2013;9(14):3877

35.Prhashanna A, Khan SA, Chen SB. Micelle morphology and chain conformation of triblock copolymers under shear: LA-DPD study. Colloids and Surfaces A: Physicochemical and Engineering Aspects. 2016;506:457–466

36.Kim S, Karrila SJ. Microhydrodynamics: Principles and Selected Applications. New York: Dover Publications; 2005

37.Bossis G, Brady JF. Dynamic simulation of sheared suspensions. I. General method. Journal of Chemical Physics. 1984;80(10):5142–5154

38.Bossis G, Brady JF. Self-diffusion of Brownian particles in concentrated suspensions under shear. Journal of Chemical Physics. 1987;87:5437

39.Durlofsky L, Brady JF, Bossis G. Dynamic simulation of hydrodynamically interacting particles. Journal of Fluid Mechanics. 1987;180:21–49

40.Brady JF, Bossis G. Stokesian dynamics. Annual Review of Fluid Mechanics. 1988;20:111–157

41.Brady JF, Phillips RJ, Lester JC, Bossis G. Dynamic simulation of hydrodynamically interacting suspensions. Journal of Fluid Mechanics. 1988;195:257–280

42.Wang M, Brady JF. Spectral ewald acceleration of stokesian dynamics for polydisperse suspensions. Journal of Computational Physics. 2016;306:443–477

43.Sierou A, Brady JF. Accelerated stokesian dynamics simulations. Journal of Fluid Mechanics. 2001;448

45.Wiener N. Differential Space. Journal of Mathematical Physics. 1923;58:31–174

46.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:771–778

47.Kim AS, Stolzenbach KD. The permeability of synthetic fractal aggregates with realistic three-dimensional structure. Journal of Colloid and Interface Science. 2002;253(2):315–328

48.Kim AS, Stolzenbach KD. Aggregate formation and collision efficiency in differential settling. Journal of Colloid and Interface Science. 2004;271(1):110–119

49.Kim AS, Kang ST. Microhydrodynamics simulation of Single-collector granular filtration. Chemistry Letters. 2012;41(10):1288–1290

51.Poschel T, Schwager T. Computational Granular Dynamics, Models and Algorithms. Springer, Berlin; 2005

52.Verwey EJ, Overbeek JTG. Theory of the Stability of Lyophobic Colloids. Amsterdam: Elsevier; 1948

53.Derjaguin BV, Landau L. Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solutions of electrolytes. Acta Physicochim USSR. 1941;14:633–662

55.Hertz H. Ueber die Verdunstung der Flussigkeiten, insbesondere des Quecksilbers, im luftleeren Raume. Annalen der Physik. 1882;253(10):177–193

56.Brilliantov NV, Spahn F, Hertzsch JM, Poschel T. Model for collisions in granular gases. Physical Review E. 1996;53(5):5382–5392

57.Brilliantov NV, Spahn F, Hertzsch JM, Poschel T. The collision of particles in granular systems. Physica A: Statistical Mechanics and its Applications. 1996;231(4):417–424

58.Ramiraz R, Poschel T, Brilliantov NV, Schwager T. Coefficient of restitution of colliding viscoelastic spheres. Physical Review E. 1999;60(4):4465–4472

59.Rosato A, Strandburg KJ, Prinz F, Swendsen RH. Why the Brazil nuts are on top: Size segregation of particulate matter by shaking. Physical Review Letters. 1987;58(10):1038

60.Hong DC, Quinn PV, Luding S. Reverse Brazil nut problem: Competition between percolation and condensation. Physical Review Letters. 2001;86(15):3423

61.Walliser H. Comment on “Reverse Brazil Nut Problem: Competition between Percolation and Condensation”. arXiv preprint cond-mat/0111559. 2001

62.Mohamed Khayetbius ME, Lauderdale BE, Nagel SR, Jaeger HM. Brazil-nut effect: Size separation of granular particles. Nature. 2001;414(6861):270–270

63.Canul-Chay G, Belmont P, Nahmad-Molinari Y, Ruiz-Suárez J. Does the reverse Brazil nut problem exist? Physical Review Letters. 2002;89(18):189601

64.Breu AP, Ensner HM, Kruelle CA, Rehberg I. Reversing the Brazil-nut effect: Competition between percolation and condensation. Physical Review Letters. 2003;90(1):014302

65.Geisel TS. Bartholomew and the Oobleck. Random House Books for Young Readers. New York; 1949

66.Roché M, Myftiu E, Johnston MC, Kim P, Stone HA. Dynamic fracture of nonglassy suspensions. Physical Review Letters. 2013;110(14)

67.Tian T, Nakano M. Design and testing of a rotational brake with shear thickening fluids. Smart Materials and Structures. 2017;26(3):035038

69.Seto R, Mari R, Morris JF, Denn MM. Discontinuous shear thickening of frictional Hard-Sphere suspensions. Physical Review Letters. 2013;111(21).

70.Hoffman RL. Discontinuous and dilatant viscosity behavior in concentrated suspensions. II. Theory and experimental tests. Journal of Colloid and Interface Science. 1974;46(3):491–506

71.Brown E, Jaeger HM. The role of dilation and confining stresses in shear thickening of dense suspensions. Journal of Rheology. 2012;56(4):875–923

72.Cui M, Emrick T, Russell TP. Stabilizing liquid drops in nonequilibrium shapes by the interfacial jamming of nanoparticles. Science. 2013;342(6157):460–463

73.Brown E, Jaeger HM. Shear thickening in concentrated suspensions: Phenomenology, mechanisms and relations to jamming. Reports on Progress in Physics. 2014;77(4):046602

74.Humphrey W, Dalke A, Schulten K. VMD: Visual molecular dynamics. Journal of Molecular Graphics. 1996;14(1):33–38

75.Sharma R, Zeller M, Pavlovic VI, Huang TS, Lo Z, Chu S, Zhao Y, Phillips JC, Schulten K. Speech/Gesture interface to a Visual-Computing environment. IEEE Computer Graphics and Applications. 2000;20:29–37

76.Eargle J, Wright D, Luthey-Schulten Z. Multiple Alignment of protein structures and sequences for VMD. Bioinformatics. 2006;22(4):504–506

77.Stukowski A. Visualization and analysis of atomistic simulation data with OVITO–the Open visualization tool. Modelling and Simulation in Materials Science and Engineering. 2010;18(1):015012

Notes

A fictional green substance in the Dr. Seuss book Bartholomew and the Oobleck.

Can You Walk on Water? (Non-Newtonian Fluid Pool) https://youtu.be/D-wxnID2q4A

Written By

Albert S. Kim and Hyeon-Ju Kim

Submitted: November 7th, 2016Reviewed: April 13th, 2017Published: September 6th, 2017