**Dissipative Dynamics of Granular Materials**

DOI: 10.5772/intechopen.69196

Dissipative Dynamics of Granular Materials

Albert S. Kim and Hyeon-Ju Kim

Additional information is available at the end of the chapter Albert S. Kim and Hyeon-Ju Kim

http://dx.doi.org/10.5772/intechopen.69196 Additional information is available at the end of the chapter

## Abstract

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.

Keywords: granular dynamics, least action principle, classical mechanics, Newton's law of motion, Hertz's law, inelasticity, compression, restitution coefficient, parallel particle dynamics, dissipative hydrodynamics

## 1. Introduction

## 1.1. Mechanics

Mechanics is 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),

© 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

which is closely related to engineering mechanics. Classical mechanics has two branches: statics and dynamics. Statics is concerned with the equilibrium of a body that is either at rest or moves with a constant velocity. Dynamics deals 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 t from their initial state at t = 0. Newton's laws for the motion of bodies can be summarized as follows:


$$F = ma \tag{1}$$

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 m can be interpreted as the object's resistance to the velocity change.)

which is closely related to engineering mechanics. Classical mechanics has two branches: statics and dynamics. Statics is concerned with the equilibrium of a body that is either at rest or moves with a constant velocity. Dynamics deals 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.

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 t from their initial state at t = 0. Newton's laws

1. 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.) 2. Newton's second law states that the net (unbalanced) force F on an object is equal to the

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

F ¼ ma ð1Þ

1.2. Principles in classical mechanics

for the motion of bodies can be summarized as follows:

product of the mass m and acceleration a of the object:

1.2.1. Governing equations

10 Granular Materials

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

$$\mathbf{M} = \mathbf{l}\alpha \tag{2}$$

where M is the torque or moment, I is 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:

$$
\begin{bmatrix} F \\ M \end{bmatrix} = \begin{bmatrix} M & 0 \\ 0 & I \end{bmatrix} \begin{bmatrix} \mathfrak{a} \\ \mathfrak{a} \end{bmatrix} \tag{3}
$$

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

$$\mathbf{a} = \ddot{\mathbf{r}} = \frac{d^2 \mathbf{r}}{dt^2} \quad \text{and} \quad \mathbf{a} = \ddot{\theta} = \frac{d^2 \theta}{dt^2} \tag{4}$$

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

$$
\sigma = \dot{r} = \frac{dr}{dt} \quad \text{and} \quad \omega = \dot{\theta} = \frac{d\theta}{dt} \tag{5}
$$

Where a and α have three components each so that [a, α] <sup>T</sup> in Eq. (3) is a vector of six components, where the superscript T indicates a transpose. For mathematical simplicity, a generalized coordinate q, generalized velocity q\_, and generalized force Q of an object are written as

$$\boldsymbol{\dot{q}} = \begin{bmatrix} \boldsymbol{r} \\ \boldsymbol{\theta} \end{bmatrix}, \quad \dot{\boldsymbol{q}} = \begin{bmatrix} \dot{\boldsymbol{r}} \\ \dot{\boldsymbol{\theta}} \end{bmatrix}, \quad \text{and} \quad \boldsymbol{\mathcal{Q}} = \begin{bmatrix} \boldsymbol{F} \\ \boldsymbol{M} \end{bmatrix} \tag{6}$$

The classical mechanics problems are usually reduced to solve for q(t) and q\_ðtÞ at time t under 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(\mathbf{x}) = m\frac{dv}{dt} \tag{7}$$

Because dx = vdt, we multiply dx with f(x) and vdf with mdv / dt to derive

$$f(\mathbf{x}) \, d\mathbf{x} = m\mathbf{v} \, d\mathbf{v} \tag{8}$$

and integrate each side from state 1 of (x1, v1) to state 2 of (x2, v2) to obtain

$$\int\_{x\_1}^{x\_2} f(\mathbf{x}) \, d\mathbf{x} = m \int\_{v\_1}^{v\_2} v d\mathbf{v} = \frac{1}{2} m v\_2^2 - \frac{1}{2} m v\_1^2 \tag{9}$$

Then, Eq. (9) can be rewritten as

$$
\Delta W = T\_2 - T\_1 \tag{10}
$$

where ΔW ¼ ð<sup>x</sup><sup>2</sup> x1 <sup>f</sup>ðxÞdx is a work done and <sup>T</sup> <sup>¼</sup> <sup>1</sup> <sup>2</sup> mv<sup>2</sup> is 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 x in Eq. (9) can be exact if the force depends on the particle position only. If so, this force is called conservative and 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(\mathbf{x}) = -\frac{dV(\mathbf{x})}{d\mathbf{x}} \quad \text{in 1-D} \tag{11}$$

$$f = -\nabla V(\mathbf{r}) \quad \text{in 3-D} \tag{12}$$

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

$$\int\_{x\_1}^{x\_2} f(\mathbf{x}) \, d\mathbf{x} = -\int\_{x\_1}^{x\_2} dV = -V(\mathbf{x}\_2) + V(\mathbf{x}\_1) \tag{13}$$

Substitution of Eq. (13) into Eq. (9) gives

Dissipative Dynamics of Granular Materials http://dx.doi.org/10.5772/intechopen.69196 13

$$T\_1 + V\_1 = T\_2 + V\_2 \tag{14}$$

$$E\_1 = E\_2 = \text{constant} \tag{15}$$

where total energy E defined 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<sup>1</sup> and x2, respectively. Note that the total energy is conserved only if the force depends on only particle location. In advanced classical mechanics, the total energy E is 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 L is defined as the difference between the kinetic and potential energies:

$$L = T - V\tag{16}$$

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 classical mechanics problems are usually reduced to solve for q(t) and q\_ðtÞ at time t under

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

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

fðxÞ ¼ m

ð<sup>v</sup><sup>2</sup> v1

vdv <sup>¼</sup> <sup>1</sup> 2 mv<sup>2</sup> <sup>2</sup> � <sup>1</sup> 2 mv<sup>2</sup>

the work done is equal to the kinetic energy change from states 1 to 2. The integration of force with respect to x in Eq. (9) can be exact if the force depends on the particle position only. If so, this force is called conservative and becomes the satisfactory condition for the energy conserva-

tion principle. A conservative f(x) can be expressed as a gradient of a scalar function V:

<sup>f</sup>ðx޼� dVðx<sup>Þ</sup>

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

ð<sup>x</sup><sup>2</sup> x1

fðxÞ dx ¼ �

Because dx = vdt, we multiply dx with f(x) and vdf with mdv / dt to derive

and integrate each side from state 1 of (x1, v1) to state 2 of (x2, v2) to obtain

fðxÞ dx ¼ m

ð<sup>x</sup><sup>2</sup> x1

<sup>f</sup>ðxÞdx is a work done and <sup>T</sup> <sup>¼</sup> <sup>1</sup>

ð<sup>x</sup><sup>2</sup> x1

Substitution of Eq. (13) into Eq. (9) gives

dv

dt <sup>ð</sup>7<sup>Þ</sup>

<sup>1</sup> ð9Þ

fðxÞ dx ¼ mv dv ð8Þ

ΔW ¼ T<sup>2</sup> � T<sup>1</sup> ð10Þ

<sup>2</sup> mv<sup>2</sup> is the kinetic energy. Eq (10) indicates that

dx in 1-D <sup>ð</sup>11<sup>Þ</sup>

dV ¼ �Vðx2Þ þ Vðx1Þ ð13Þ

f ¼ �∇VðrÞ in 3-D ð12Þ

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

than Newton's second law?"

12 Granular Materials

law in one-dimensional space:

Then, Eq. (9) can be rewritten as

ð<sup>x</sup><sup>2</sup> x1

where ΔW ¼

The most general formulation of the law governing the motion of mechanical systems is the principle of least action or Hamilton's principle [1]. According to this principle, a mechanical system is characterized using a definite Lagrangian function Lðq1, q2, …, qs, q\_ <sup>1</sup>, q\_ <sup>2</sup>, …, q\_ <sup>s</sup>, 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<sup>1</sup> and t2, particle positions are defined by two sets of the generalized coordinates, q(t1) and q(t2). The condition is that the system moves between these two positions, minimizing the integral

$$S = \int\_{t\_1}^{t\_2} L(q\_\prime \dot{q}\_\prime t) \mathbf{d}t \tag{17}$$

to the least possible value. The integral of Eq. (17) is called the action. Note that the Lagrangian contains generalized coordinates and velocities, q and 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 S is a minimum. This means that S changes when q(t) is replaced by

$$q(t) + \delta q(t) \tag{18}$$

where δq(t), called a variation of q(t), is a function, assumed to be small everywhere in the interval of time from t<sup>1</sup> to t2. Since Eq. (18) must include the values of q(t1) and q(t2), we can now conclude that

$$
\delta q(t\_1) = \delta q(t\_2) = 0 \tag{19}
$$

In this case, the change in S when q is replaced by q + δq is equal to

$$\int\_{t\_1}^{t\_2} L(q + \delta q, \dot{q} + \delta \dot{q}, t) - \int\_{t\_1}^{t\_2} L(q, \dot{q}, t) \tag{20}$$

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

$$
\delta \mathcal{S} = \delta \int\_{t\_1}^{t\_2} L(\eta, \dot{\eta}, t) = 0 \tag{21}
$$

or equivalently

$$\int\_{t\_1}^{t\_2} \left( \frac{\partial L}{\partial q} \delta q + \frac{\partial L}{\partial \dot{q}} \delta \dot{q} \right) \mathbf{d}t = 0 \tag{22}$$

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

$$
\delta S = \left[\frac{\partial L}{\partial \dot{q}} \delta q\right]\_{t\_1}^{t\_2} + \int\_{t\_1}^{t\_2} \left(\frac{\partial L}{\partial q} - \frac{\mathbf{d}}{\mathbf{d}t} \frac{\partial L}{\partial \dot{q}}\right) \delta q \mathbf{d}t \tag{23}
$$

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

$$
\left[\frac{\partial L}{\partial \dot{\eta}} \delta q\right]\_{t\_1}^{t\_2} = 0 \tag{24}
$$

and the remaining integral is

$$\int\_{t\_1}^{t\_2} \left( \frac{\partial L}{\partial q} - \frac{\mathbf{d}}{\mathbf{d}t} \frac{\partial L}{\partial \dot{q}} \right) \delta q \mathbf{d}t = 0,\tag{25}$$

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

$$
\frac{\partial}{\partial t} \left( \frac{\partial L}{\partial \dot{q}} \right) - \frac{\partial L}{\partial q} = 0 \tag{26}
$$

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

$$\frac{\partial}{\partial t} \left( \frac{\partial L}{\partial \dot{q}\_i} \right) - \frac{\partial L}{\partial q\_i} = 0 \quad (i = 1, 2, \dots, s) \tag{27}$$

where s is 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 Np objects, s, is equal to 6Np.

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

$$\frac{\partial}{\partial \mathbf{d}} \left( \frac{\partial L}{\partial \dot{q}\_i} \right) - \frac{\partial L}{\partial q\_i} = Q\_i^\dagger \quad (i = 1, 2, \dots, s). \tag{28}$$

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 L instead of force vectors.

## 2. Particle dynamics simulation methods

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

δqðt1Þ ¼ δqðt2Þ ¼ 0 ð19Þ

Lðq, q, t \_ Þ ð20Þ

Lðq, q, t \_ Þ ¼ 0 ð21Þ

dt ¼ 0 ð22Þ

¼ 0 ð24Þ

δqdt ¼ 0, ð25Þ

<sup>∂</sup><sup>q</sup> <sup>¼</sup> <sup>0</sup> <sup>ð</sup>26<sup>Þ</sup>

¼ 0 ði ¼ 1; 2;…, sÞ ð27Þ

δqdt ð23Þ

In this case, the change in S when q is replaced by q + δq is equal to ðt2 t1

terms. Then, the principle of least action may be written in the form

ðt2 t1

<sup>δ</sup><sup>S</sup> <sup>¼</sup> <sup>∂</sup><sup>L</sup> ∂q\_ δq � �<sup>t</sup><sup>2</sup>

or equivalently

14 Granular Materials

and the remaining integral is

Eq. (25) is zero. Thus, we have

Lðq þ δq, q\_ þ δq, t \_ Þ �

ðt2 t1

� �

∂L <sup>∂</sup><sup>q</sup> � <sup>d</sup> dt ∂L ∂q\_

t1

� �

δS ¼ δ

∂L ∂q δq þ ∂L ∂q\_ δq\_

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

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

ðt2 t1 ∂L <sup>∂</sup><sup>q</sup> � <sup>d</sup> dt ∂L ∂q\_

d dt

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

total degrees of freedom of a system containing Np objects, s, is equal to 6Np.

� <sup>∂</sup><sup>L</sup> ∂qi

∂L ∂q\_ i � �

d dt t1 þ ðt2 t1

∂L ∂q\_ δq � �<sup>t</sup><sup>2</sup>

� �

∂L ∂q\_ � �

which must vanish for all values of δq. This can be satisfied if and only if the integrand of

where s is 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

� ∂L

We expand the difference in powers of δq and δq\_ in the integrand and leave only the first-order

ðt2 t1

> 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 dt must 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 + δt can be expanded at time t using Taylor's series:

$$r(t + \delta t) = r(t) + \mathbf{v}(t)\delta t + \frac{1}{2}\mathbf{a}(t)\delta t^2 + \mathcal{O}(\delta t^3) \tag{29}$$

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

$$r(t - \delta t) = r(t) - \mathbf{v}(t)\delta t + \frac{1}{2}\mathbf{a}(t)\delta t^2 - O(\delta t^3) \tag{30}$$

The position r and velocity v at the future time t + δt can be calculated by adding and subtracting Eqs. (29) and (30) such as

$$\mathbf{r}(t+\delta t) = \mathbf{2}\mathbf{r}(t) - \mathbf{r}(t-\delta t) + \delta t^2 \mathbf{a}(t) \tag{31}$$

and

$$v(t) = \frac{r(t + \delta t) - r(t - \delta t)}{2\delta t} \tag{32}$$

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 a obtained 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:

$$m\frac{d\mathbf{v}}{dt} = -\xi\mathbf{v} + \mathbf{f}(\mathbf{r}) + \mathbf{f}'(t) \tag{33}$$

where – ξv is 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 <sup>0</sup> (t) is the random fluctuation force caused by adjacent, rapidly-moving solvent molecules. It is assumed that f <sup>0</sup> is a Gaussian process with infinitely small correlation time, which renders the famous fluctuation-dissipation theorem [17]

$$
\langle \mathbf{f}'(t\_1) \cdot \mathbf{f}'(t\_2) \rangle = 2\xi k\_B T \delta(t\_1 - t\_2) \tag{34}
$$

$$
\langle f'(t) \rangle = 0 \tag{35}
$$

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

$$
\left<\frac{d\boldsymbol{v}}{dt}\right> \simeq -\left<\frac{\boldsymbol{v}}{m/\xi}\right> \simeq -\frac{\langle\boldsymbol{v}\rangle}{\pi} \tag{36}
$$

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 dt in 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 πdpμ, where dp 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

$$D\_{\rm SE} = \frac{k\_B T}{\xi} = \frac{k\_B T}{\Im \pi \mu d\_p} \tag{37}$$

Although BD adopts the Oseen tensor, it still treats a particle as a point mass (like MD). From BD simulations, the hydrodynamic diameter dp 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

rðt þ δtÞ ¼ 2rðtÞ � rðt � δtÞ þ δt

<sup>v</sup>ðtÞ ¼ <sup>r</sup>ð<sup>t</sup> <sup>þ</sup> <sup>δ</sup>tÞ � <sup>r</sup>ð<sup>t</sup> � <sup>δ</sup>t<sup>Þ</sup>

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 a obtained 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

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

dt ¼ �ξ<sup>v</sup> <sup>þ</sup> <sup>f</sup>ðrÞ þ <sup>f</sup> <sup>0</sup>

where – ξv is 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

random fluctuation force caused by adjacent, rapidly-moving solvent molecules. It is assumed that f <sup>0</sup> is a Gaussian process with infinitely small correlation time, which renders the famous

form), f(r) is the conservative force derived from the potential energy function, and f <sup>0</sup>

m dv

developed for particles treated as point masses.

2.2. Brownian dynamics

Eq. (1) to Langevin's equation:

fluctuation-dissipation theorem [17]

and

16 Granular Materials

2

aðtÞ ð31Þ

<sup>2</sup>δ<sup>t</sup> <sup>ð</sup>32<sup>Þ</sup>

ðtÞ ð33Þ

(t) is the

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 i from other particles is written in the form of:

$$m\mathfrak{a}\_{i} = \sum\_{j=1, j \neq i}^{N\_{\mathbb{P}}} (\boldsymbol{\mathcal{F}}\_{ij}^{\boldsymbol{P}} + \boldsymbol{\mathcal{F}}\_{ij}^{\boldsymbol{D}} + \boldsymbol{\mathcal{F}}\_{ij}^{\boldsymbol{R}}) \tag{38}$$

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

$$\boldsymbol{F}\_{\vec{\boldsymbol{\eta}}}^{D} = -\gamma \boldsymbol{\omega}\_{D}(\boldsymbol{r}\_{\vec{\boldsymbol{\eta}}}) (\boldsymbol{\varrho}\_{\vec{\boldsymbol{\eta}}} \cdot \boldsymbol{\nu}\_{\vec{\boldsymbol{\eta}}}) \boldsymbol{v}\_{\vec{\boldsymbol{\eta}}} \tag{39}$$

$$\mathbf{F}^{\underset{\text{i}\dot{\mathbf{R}}}{\text{R}}}\_{\dot{\mathbf{i}}\dot{\mathbf{j}}} = \sigma \boldsymbol{\omega}\_{\mathbf{R}}(\mathbf{r}\_{i\dot{\mathbf{j}}}) \mathbf{e}\_{i\dot{\mathbf{j}}} \zeta\_{i\dot{\mathbf{j}}} \tag{40}$$

where rij = |r<sup>i</sup> � <sup>r</sup>j|, <sup>v</sup>ij <sup>=</sup> <sup>v</sup><sup>i</sup> � <sup>v</sup>j, <sup>e</sup>ij = (r<sup>i</sup> � <sup>r</sup>j)/rij, and <sup>σ</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffi 2kBTγ p . Importantly, ζij = ζji is a Gaussian white-noise ther of zero mean and unit variance. The stochastic differential equation (SDE) of DPD consists of

$$d\mathbf{r}\_i = \mathbf{v}\_i dt\tag{41}$$

$$md\nu\_{i} = \sum\_{j} [\mathbf{F}\_{i\bar{j}}^{P} + \mathbf{F}\_{i\bar{j}}^{D}]dt + \sum\_{j} \sigma \omega\_{\mathbb{R}}(r\_{i\bar{j}}) \mathbf{e}\_{i\bar{j}} d\mathcal{W}\_{i\bar{j}} \tag{42}$$

where dWij (= dWji) is independent increment and of the Wiener process, satisfying

$$d\mathcal{W}\_{\vec{\eta}}d\mathcal{W}\_{kl} = (\delta\_{ik}\delta\_{jl} + \delta\_{jl}\delta\_{jk})dt\tag{43}$$

i.e., dWij(t) is an infinitesimal of order 1/2 (i.e., proportional to ffiffiffiffi <sup>δ</sup><sup>t</sup> <sup>p</sup> ).

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 ffiffiffiffi dt <sup>p</sup> . Hydrodynamic resistances are presumed to be pairwise and described as (intuitively chosen) simple functions of the inter-particle distance r:

$$
\omega\_D = \omega\_R^2 = \begin{cases}
\left(1 - r/r\_c\right)^2 & \text{for} \quad r < r\_c \\
0 & \text{otherwise}
\end{cases}
\tag{44}
$$

where rc 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<sup>∞</sup> is built using the pairwise superposition of the two-body mobility matrix. Its product with the hydrodynamic force F<sup>H</sup> gives the relative velocities of particles to the fluid flow, ΔU[37–41]:

$$
\Delta \mathbf{U} = \left[ \mathcal{M}^{\circ \circ} \right] \cdot \mathbf{F}^{\mathrm{H}} \tag{45}
$$

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 nearfield lubrication such as

$$\left[\mathcal{R}\right] = \left[\mathcal{M}^{\circ}\right]^{-1} - \left[\Delta \mathcal{R}\_{2B}\right] \tag{46}$$

where <sup>½</sup>ΔR2B� indicates a pairwise addition of the exact two-body resistance matrix <sup>R</sup><sup>ð</sup>2<sup>Þ</sup> i,j subtracted by the inverse of two-body far-field mobility matrix, M<sup>ð</sup>2<sup>Þ</sup> ij :

$$\mathbb{E}\left[\Delta \mathcal{R}\_{2B}\right] = \sum\_{i,j} \left( \mathcal{R}\_{i,j}^{(2)} - [\mathcal{M}\_{ij}^{(2)}]^{-1} \right) \tag{47}$$

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 drifting motion 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 dt is 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

dr<sup>i</sup> ¼ vidt ð41Þ

dWijdWkl ¼ ðδikδjl þ δjlδjkÞdt ð43Þ

<sup>δ</sup><sup>t</sup> <sup>p</sup> ).

<sup>Δ</sup><sup>U</sup> ¼ ½M<sup>∞</sup>� � <sup>F</sup><sup>H</sup> <sup>ð</sup>45<sup>Þ</sup>

<sup>2</sup> for r < rc

0 otherwise

dt

<sup>p</sup> . Hydrodynamic resis-

ð44Þ

σωRðrijÞeijdWij ð42Þ

mdv<sup>i</sup> <sup>¼</sup> <sup>X</sup> j <sup>½</sup><sup>F</sup> <sup>P</sup> ij <sup>þ</sup> <sup>F</sup><sup>D</sup>

i.e., dWij(t) is an infinitesimal of order 1/2 (i.e., proportional to ffiffiffiffi

<sup>ω</sup><sup>D</sup> <sup>¼</sup> <sup>ω</sup><sup>2</sup>

the inter-particle distance r:

18 Granular Materials

2.4. Stokesian dynamics

particles to the fluid flow, Δ

field lubrication such as

U[37–41]:

tional to dt, and due to the random fluctuating, force is linear on ffiffiffiffi

ij �dt <sup>þ</sup><sup>X</sup> j

In Eq. (42), the momentum changes due to the conservative and dissipative forces are propor-

tances are presumed to be pairwise and described as (intuitively chosen) simple functions of

where rc 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 automati-

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<sup>∞</sup> is built using the pairwise superposition of the two-body mobility matrix. Its product with the hydrodynamic force F<sup>H</sup> gives the relative velocities of

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-

<sup>R</sup> <sup>¼</sup> <sup>ð</sup><sup>1</sup> � <sup>r</sup>=rc<sup>Þ</sup>

�

cally satisfied [17] so that the thermodynamics of the system are well described.

where dWij (= dWji) is independent increment and of the Wiener process, satisfying

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 Np 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

$$\mathbf{M} \cdot d\mathbf{u} = [\mathbf{Q}^\mathbf{v} - \mathbf{R} \cdot (\mathbf{u} - \mathbf{U})]dt + \mathbf{B} \cdot dW = \mathbf{Q}dt\tag{48}$$

where M is the mass/moment-of-inertia matrix of 6Np � 6Np dimension, u and U are translational/ rotational velocities of particles and the fluid, respectively, Q<sup>p</sup> is the conservative force/torque vector, R is the grand resistance matrix, B is the Brownian matrix of zero mean and finite variance, i.e., 〈B〉 = 0 and 〈Btr � <sup>B</sup>〉 = 2kBT, where kB is the Boltzmann constant, and <sup>d</sup><sup>W</sup> is the Ito-Wiener process of 6Np elements [45, 46]. The generalized force/torque is defined, similar to Eq. (42), as

$$\mathcal{Q} = \mathcal{Q}^{\flat} - \mathbf{R} \cdot (\mathbf{v} - \mathbf{U}) + \mathbf{B} \cdot \mathbf{w} dt \tag{49}$$

where <sup>w</sup> <sup>¼</sup> <sup>d</sup>W<sup>=</sup> ffiffiffiffi dt <sup>p</sup> . 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 i and j of mass mi and mj, respectively, undergo an inelastic collision, as shown in Figure 1. The relative velocity of particles i and j, denoted as gij, at the point of contact is determined by the translational and rotational particle velocities:

$$\mathbf{g}\_{i\dot{\jmath}} = \mathbf{v}\_{i\dot{\jmath}} - (a\_i \boldsymbol{\omega}\_i + a\_{\dot{\jmath}} \boldsymbol{\omega}\_{\dot{\jmath}}) \times \mathbf{n}\_{i\dot{\jmath}} \tag{50}$$

where vij = v<sup>i</sup> � v<sup>j</sup> is the relative velocity of the center of mass of particle i to j of velocity v<sup>i</sup> and vj, respectively. In general, particles are polydispersed, and ai and aj are the radii of particles i and j, respectively. The normal and tangential collision velocities given by the projections of gij are

$$\mathbf{g}\_{i\dot{\jmath}}^{n} = (\mathbf{g}\_{i\dot{\jmath}} \cdot \mathbf{n}\_{i\dot{\jmath}}) \mathbf{n}\_{i\dot{\jmath}} \tag{51}$$

and

$$\mathbf{g}\_{i\dagger}^{t} = -\mathbf{n}\_{i\dagger} \times (\mathbf{n}\_{i\dagger} \times \mathbf{g}\_{i\dagger}) \tag{52}$$

respectively. The coefficients of restitution in normal and tangential direction, E <sup>n</sup> and E t , are defined as

$$(\mathbf{g}\_{ij}^n)' = -\epsilon^n \mathbf{g}\_{ij}^n \tag{53}$$

$$(\mathbf{g}\_{ij}^{t})' = +\epsilon^{t}\mathbf{g}\_{ij}^{t} \tag{54}$$

where 0 ≤ E <sup>n</sup> <sup>≤</sup> 1 and �<sup>1</sup> <sup>≤</sup> <sup>E</sup> nt ≤ 1. The primed and unprimed variables indicate the pre- and postcollision quantities, respectively. Then, the velocities of particles i and j after the collision can be represented as functions of E n , E t , g<sup>n</sup> ij, g<sup>t</sup> ij, mi , and mj. Now, for particles i and j, we finally have

$$\mathbf{v}'\_{i} = \mathbf{v}\_{i} - \left(\frac{\mu\_{i\bar{\boldsymbol{\eta}}}}{m\_{i}}\right) \left[ (1 + \boldsymbol{\epsilon}'') \mathbf{g}^{\boldsymbol{n}}\_{i\bar{\boldsymbol{\eta}}} + \frac{1 - \boldsymbol{\epsilon}^{t}}{1 + 1/\bar{\boldsymbol{\beta}}} \mathbf{g}^{t}\_{i\bar{\boldsymbol{\eta}}} \right] \tag{55}$$

$$\mathbf{v}'\_{\dot{j}} = \mathbf{v}\_{\dot{j}} + \left(\frac{\mu\_{\dot{\imath}\dot{\jmath}}}{m\_{\dot{\jmath}}}\right) \left[ (1 + \epsilon^{\imath}) \mathbf{g}^{\imath}\_{\dot{\imath}\dot{\jmath}} + \frac{1 - \epsilon^{t}}{1 + 1/\dot{\jmath}} \mathbf{g}^{t}\_{\dot{\imath}\dot{\jmath}} \right] \tag{56}$$

Figure 1. Collision of sphere i and j.

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

In the dynamic motion of granular particles, ballistic collisions are one of the most fundamental and important interactions. Suppose two (spherical) particles i and j of mass mi and mj, respectively, undergo an inelastic collision, as shown in Figure 1. The relative velocity of particles i and j, denoted as gij, at the point of contact is determined by the translational and

where vij = v<sup>i</sup> � v<sup>j</sup> is the relative velocity of the center of mass of particle i to j of velocity v<sup>i</sup> and vj, respectively. In general, particles are polydispersed, and ai and aj are the radii of particles i and j, respectively. The normal and tangential collision velocities given by the projections of gij are

gn

respectively. The coefficients of restitution in normal and tangential direction, E

ðgn ijÞ <sup>0</sup> ¼ �E <sup>n</sup>g<sup>n</sup>

ðgt ijÞ <sup>0</sup> ¼ þE t gt

collision quantities, respectively. Then, the velocities of particles i and j after the collision can be

ð1 þ E <sup>n</sup>Þg<sup>n</sup> ij þ

ð1 þ E <sup>n</sup>Þg<sup>n</sup> ij þ

gt

gij ¼ vij � ðaiω<sup>i</sup> þ ajωjÞ � nij ð50Þ

ij ¼ ðgij � nijÞnij ð51Þ

ij ¼ �nij � ðnij � gijÞ ð52Þ

, and mj. Now, for particles i and j, we finally have

gt ij

gt ij

nt ≤ 1. The primed and unprimed variables indicate the pre- and post-

<sup>1</sup> � <sup>E</sup><sup>t</sup> <sup>1</sup> <sup>þ</sup> <sup>1</sup>=~<sup>J</sup>

<sup>1</sup> � <sup>E</sup><sup>t</sup> <sup>1</sup> <sup>þ</sup> <sup>1</sup>=~<sup>J</sup>

" #

" #

ij ð53Þ

ij ð54Þ

<sup>n</sup> and E t , are

ð55Þ

ð56Þ

to simulate a reasonably large numbers of particles.

2.6. Discrete element method

20 Granular Materials

rotational particle velocities:

and

defined as

where 0 ≤ E

<sup>n</sup> <sup>≤</sup> 1 and �<sup>1</sup> <sup>≤</sup> <sup>E</sup>

n , E t , g<sup>n</sup> ij, g<sup>t</sup> ij, mi

v0

v0 <sup>j</sup> ¼ v<sup>j</sup> þ

<sup>i</sup> <sup>¼</sup> <sup>v</sup><sup>i</sup> � <sup>μ</sup>ij

mi � �

μij mj � �

represented as functions of E

$$
\omega\_i' = \omega\_i + \left(\frac{\mu\_{\vec{\imath}}}{m\_l a\_i}\right) \left(\frac{1 - \epsilon^t}{1 + \tilde{\jmath}}\right) \mathfrak{n}\_{\vec{\imath}\vec{\jmath}} \times \mathfrak{g}\_{\vec{\imath}\vec{\jmath}}^t \tag{57}
$$

$$
\omega\_{\rangle}^{\prime} = \omega\_{\rangle} + \left(\frac{\mu\_{\stackrel{\scriptstyle \partial}{\dot{\phantom{a}}}}}{m\_{\stackrel{\scriptstyle \partial}{}}a\_{\stackrel{\scriptstyle \partial}{}}}\right)\left(\frac{1-\epsilon^{t}}{1+\tilde{\phantom{\phantom{a}}}}\right)\mathfrak{n}\_{\stackrel{\scriptstyle \partial}{}} \times \mathfrak{g}\_{\stackrel{\scriptstyle \partial}{}}^{t} \tag{58}
$$

where for <sup>k</sup> <sup>=</sup> <sup>i</sup>, <sup>j</sup>. Here, Jk is the mass moment of inertia and <sup>~</sup><sup>J</sup> <sup>k</sup> <sup>¼</sup> Jk=mka<sup>2</sup> <sup>k</sup> is its dimensionless form. A spherical particle has Jk = 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 k is represented as

$$T\_k^{(\prime)} = \frac{1}{2} \mathfrak{m}\_k \mathfrak{v}\_k^{(\prime)} \cdot \mathfrak{v}\_k^{(\prime)} + \frac{1}{2} J\_k \mathfrak{w}\_k^{(\prime)} \cdot \mathfrak{w}\_k^{(\prime)} \tag{59}$$

One can calculate the energy loss of particle k, ΔTk ¼ T<sup>0</sup> <sup>k</sup> � Tk, using the pre- and post-collision velocities of ðvk, ωkÞ and ðv<sup>0</sup> k, ω<sup>0</sup> <sup>k</sup>Þ, 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., ω<sup>D</sup> and ω<sup>R</sup> 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.


Table 1. 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" dt indicates that dt can be arbitrarily determined for accurate and fast calculation regardless of the particle relaxation time τ.

## 3. Granular dynamics: theory and simulations

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

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.

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., ω<sup>D</sup> and ω<sup>R</sup> 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

MD BD SD DPD DHD [44] DEM

Gov. Eq. F = ma Langevin Langevin SDE SDE various Acceleration Yes Possible Possible Yes Yes Yes Force Conservative Yes Yes Yes Yes Yes Yes

Time interval dt Any δt ≫ τ δt ≫ τ Any Any Any

arbitrarily determined for accurate and fast calculation regardless of the particle relaxation time τ.

Table 1. 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" dt indicates that dt can be

Brownian No Simple Rigorous Approx. Rigorous No need Hydrodynamics No Constant Accurate Pairwise Accurate No Constraint Yes Yes Yes Yes Yes [54] Yes Collision No No No No (possible) Yes

2.7. Overview

22 Granular Materials

DEM method.

#### 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 = \begin{cases} 0 & r > 2a \\ \Leftrightarrow & r < 2a \end{cases} \tag{60}$$

where a is the particle radius. The mathematical discontinuity of the wall potential at r = 2a indicates 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.,

$$
\xi\_{ij} \equiv a\_i + a\_j - |\mathbf{r}\_i - \mathbf{r}\_j| > 0 \tag{61}
$$

where ξij is the mutual compression of particles i and 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 i from particle j (conceptually denoted as i j) is described by

$$F\_{ij} = \begin{cases} F\_{ij}^{\
u} + F\_{ij}^{\
u} & \text{for} \quad \xi\_{ij} > 0\\ 0 & \text{for} \quad \xi\_{ij} \le 0 \end{cases} \tag{62}$$

where F<sup>n</sup> ij and F<sup>t</sup> ij are the normal and tangential components of the contact force, respectively. For simplicity, Eq. (62) can be rewritten as

$$F\_{\vec{\eta}} = (F\_{\ \vec{\eta}}^n + F\_{\ \vec{\eta}}^t)H(\xi\_{\ \vec{\eta}}) \tag{63}$$

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

$$H(\mathbf{x}) = \begin{cases} 1 & \text{for } \quad \mathbf{x} > 0 \\ 0 & \text{otherwise} \end{cases} \tag{64}$$

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

$$A \times (\mathbf{B} \times \mathbf{C}) = \mathbf{B}(A \cdot \mathbf{C}) - \mathbf{C}(A \cdot \mathbf{B}) \tag{65}$$

one can replace A and B by n and C by Fij to write

$$
\mathfrak{n} \times (\mathfrak{n} \times \boldsymbol{F}\_{\vec{\imath}\boldsymbol{j}}) = \mathfrak{n}(\mathfrak{n} \cdot \boldsymbol{F}\_{\vec{\imath}\boldsymbol{j}}) - \boldsymbol{F}\_{\vec{\imath}\boldsymbol{\not\jmath}} \tag{66}
$$

$$F\_{\vec{\eta}} = \mathfrak{n}(\mathfrak{n} \cdot F\_{\vec{\eta}}) - \mathfrak{n} \times (\mathfrak{n} \times F\_{\vec{\eta}}) \tag{67}$$

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

$$\boldsymbol{F}\_{\vec{\text{ij}}}^{n} = (\mathfrak{n} \cdot \boldsymbol{F}\_{\vec{\text{ij}}}) \mathfrak{n} \tag{68}$$

$$\boldsymbol{F}^t\_{\vec{\boldsymbol{\eta}}} = (\mathfrak{n} \times \boldsymbol{F}\_{\vec{\boldsymbol{\eta}}}) \times \mathfrak{n} \tag{69}$$

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

$$F^{\imath} = \frac{2Y\sqrt{a\_{\text{eff}}}}{3(1-\nu^2)} \left(\xi^{3/2} + A\sqrt{\xi}\frac{d\xi}{dt}\right) \tag{70}$$

where Y and ν are Young's modulus and Poisson's ratio, respectively, A is the dissipative constant being a function of material viscosity, and aeff is the effective radius that can be interpreted as a harmonic sum ai and aj, i.e., a�<sup>1</sup> eff <sup>¼</sup> <sup>a</sup>�<sup>1</sup> <sup>i</sup> <sup>þ</sup> <sup>a</sup>�<sup>1</sup> <sup>j</sup> . Parameter A explains 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 A needs 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

Dissipative Dynamics of Granular Materials http://dx.doi.org/10.5772/intechopen.69196 25

$$F\_{ij}^n = \frac{4\sqrt{a\_{\text{eff}}}}{3} \left(\frac{1-\nu\_i^2}{Y\_i} + \frac{1-\nu\_j^2}{Y\_j}\right)^{-1} \xi^{3/2} \tag{71}$$

for particles of different Y and ν values. Assuming the dissipative constant A is also particlespecific and additive, the most general expression of the normal force between two viscoelastic spheres may be [51]

$$F\_{\vec{\eta}}^{\imath} = \frac{4\sqrt{a\_{\text{eff}}}}{3} \left(\frac{1-\nu\_i^2}{Y\_i} + \frac{1-\nu\_j^2}{Y\_j}\right)^{-1} \left(\xi^{3/2} + \frac{A\_{\text{i}}+A\_{\text{j}}}{2}\dot{\xi}\sqrt{\xi}\right) \tag{72}$$

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 = tc, after which restoration begins. When the relative velocities between the two particles at t = 0 and t = tc are <sup>ξ</sup>\_ð0<sup>Þ</sup> and <sup>ξ</sup>\_ðtcÞ, respectively, then the normal restitution coefficient can be calculated by measuring velocities <sup>ξ</sup>\_ð0<sup>Þ</sup> and <sup>ξ</sup>\_ðtcÞ, i.e.,

$$
\epsilon^n = \dot{\xi}(t\_c) / \dot{\xi}(0) \tag{73}
$$

This indicates that E <sup>n</sup> depends on the relative collision velocity, unless A = 0. Theoretically, at least three experiments are required for collisions between particles i and j for pairs of (i, j), (i, i), and (j, j). Then, Ai and Aj 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

$$\boldsymbol{\nu}\_{\rm rel}^{t} = (\boldsymbol{\nu}\_{j} - \boldsymbol{\nu}\_{i}) \cdot \boldsymbol{\mathcal{e}}\_{ij}^{t} + a\_{i}\boldsymbol{\omega}\_{i} + a\_{j}\boldsymbol{\omega}\_{j} \tag{74}$$

which provides the tangential force of

compression can be interpreted as the overlapped surface-to-surface distance. The force acting

ij <sup>þ</sup> <sup>F</sup> <sup>t</sup>

<sup>H</sup>ðxÞ ¼ 1 for <sup>x</sup> <sup>&</sup>gt; <sup>0</sup>

In three-dimensional space, vector quantities between particles i and j can be decomposed into

�

0 othersies

A � ðB � CÞ ¼ BðA � CÞ � CðA � BÞ ð65Þ

n � ðn � FijÞ ¼ nðn � FijÞ � Fij ð66Þ

Fij ¼ nðn � FijÞ � n � ðn � FijÞ ð67Þ

ij ¼ ðn � FijÞn ð68Þ

ij ¼ ðn � FijÞ � n ð69Þ

<sup>j</sup> . Parameter A explains the depen-

ij for ξij > 0 0 for ξij ≤ 0

ij are the normal and tangential components of the contact force, respectively.

ijÞHðξijÞ ð63Þ

ð62Þ

ð64Þ

ð70Þ

ij <sup>þ</sup> <sup>F</sup> <sup>t</sup>

on particle i from particle j (conceptually denoted as i j) is described by

<sup>F</sup>ij <sup>¼</sup> <sup>F</sup> <sup>n</sup>

the normal and tangential directions. Using the mathematical identity

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

F t

<sup>F</sup><sup>n</sup> <sup>¼</sup> <sup>2</sup><sup>Y</sup> ffiffiffiffiffiffi

aeff p

F <sup>n</sup>

The contact force between elastic spheres was originally developed by Hertz [55] and was later

<sup>3</sup>ð<sup>1</sup> � <sup>ν</sup><sup>2</sup><sup>Þ</sup> <sup>ξ</sup><sup>3</sup>=<sup>2</sup> <sup>þ</sup> <sup>A</sup> ffiffiffi

where Y and ν are Young's modulus and Poisson's ratio, respectively, A is the dissipative constant being a function of material viscosity, and aeff is the effective radius that can be

dence 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 A needs 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

eff <sup>¼</sup> <sup>a</sup>�<sup>1</sup>

<sup>ξ</sup> <sup>p</sup> <sup>d</sup><sup>ξ</sup> dt

� �

<sup>i</sup> <sup>þ</sup> <sup>a</sup>�<sup>1</sup>

where F<sup>n</sup>

24 Granular Materials

ij and F<sup>t</sup>

For simplicity, Eq. (62) can be rewritten as

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

one can replace A and B by n and C by Fij to write

generalized for viscoelastic (damped) particles [56, 57] as

interpreted as a harmonic sum ai and aj, i.e., a�<sup>1</sup>

�

<sup>F</sup>ij ¼ ð<sup>F</sup> <sup>n</sup>

$$F^t = -\text{sign}(\boldsymbol{\nu}\_{\text{rel}}^t) \cdot \min(\boldsymbol{\gamma}^t | \boldsymbol{\nu}\_{\text{rel}}^t | \boldsymbol{\mu} | F^t) \tag{75}$$

where γ<sup>t</sup> 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 |Ft | ≤ μ|Fn |, 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 Fn is constant, Eq. (75) provides an inconsistent value of Ft because Fn (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]:

$$F\_{ij}^t = -\mu F\_{ij}^t \left(\frac{\zeta\_{ij}}{\zeta\_0} - \left\lfloor \frac{\zeta\_{ij}}{\zeta\_0} \right\rfloor \right) \tag{76}$$

where ζ is the relative tangential shift, ζ<sup>0</sup> is its macroscopic maximum value, and ⌊x⌋ denotes the integer of x. Typical values of ζ<sup>0</sup> / aeff range from 10�<sup>7</sup> to 10�<sup>3</sup> . Similar to Eq. (73), the tangential coefficient of restitution can be calculated as

$$
\epsilon^t = \dot{\check{\zeta}}(t\_c) / \dot{\check{\zeta}}(0) \tag{77}
$$

If the friction coefficient μ is known, experimental measurement of E <sup>t</sup> 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 tc, while tc 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 tc in compressing/ restoring phases, and much longer than tc 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 thickening indicates 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

<sup>1</sup> A fictional green substance in the Dr. Seuss book Bartholomew and the Oobleck.

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

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

E <sup>t</sup> <sup>¼</sup> \_ <sup>ζ</sup>ðtcÞ<sup>=</sup> \_

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 tc, while tc 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 tc in compressing/ restoring phases, and much longer than tc for collision events for fluidized granules of a low

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

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

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

If the friction coefficient μ is known, experimental measurement of E

concentration.

26 Granular Materials

3.1.2. Shear stress

1

2

ζð0Þ ð77Þ

<sup>t</sup> allows us to inversely

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

Figure 2. Simulation of granular damping to a sudden impact: an intruder (a) approaches with high speed, (b) collides with a few granules on the packed surface, and (c) penetrates the loosely packed inelastic (energy absorbing) granules. Snapshots are visualized using Visual Molecular Dynamics (VMD) [74–76]. Radii of the intruder and packed spheres are 2 and 1 mm, respectively, and their specific gravity values are commonly 2.75. The spheres and the walls have Poisson's ratio of 0.4 and 0.6, and Young's modulus of 4.0 107 Pa and 1.0 109 Pa, respectively. The introduced energy dissipation rate <sup>A</sup> was set at 2.5 <sup>10</sup><sup>5</sup> , for both spheres and walls [56].

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 vz = 13.00 <sup>10</sup><sup>3</sup> m/s and spins with angular velocity of ω<sup>z</sup> = 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 A reduces the impact of the one-tomany 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

Figure 3. Force chain visualization of the intruder using open visualization tool (OVITO).

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

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 vz = 13.00 <sup>10</sup><sup>3</sup> m/s and spins with

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 A reduces the impact of the one-tomany 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

angular velocity of ω<sup>z</sup> = 0.56 rad/s.

28 Granular Materials

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 i and j is

$$
\sigma\_{i\dot{\jmath}} = (\mathbf{r}\_i - \mathbf{r}\_{\dot{\jmath}})^2 - d\_{i\dot{\jmath}}^2 = \mathbf{0} \tag{78}
$$

where dij is the fixed distance between particles i and j, usually an average of two contacting rigid bodies of diameters di and dj, i.e., dij = (di + dj)/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 Np particles as

$$\Phi = \frac{1}{2} \sum\_{i=1}^{N\_p} \sum\_{j>i}^{N\_p} \lambda\_{ji} \sigma\_{ij} \tag{79}$$

where λji is a symmetric Lagrange's multiplier and <sup>1</sup> <sup>2</sup> in the front of the summation is by convention. Exchanging positions of particles i and 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 j can be derived as a negative derivative of the constraint potential:

$$\mathbf{F}\_{\dot{j}}^{\mathcal{C}} = -\nabla\_{\dot{j}} \Phi = \sum\_{i \neq j} \lambda\_{ji} (\mathbf{r}\_i - \mathbf{r}\_j) \tag{80}$$

where i runs from 1 to Np except for the case i = j (if so, λii = 0), or simply

$$\boldsymbol{F}\_{\dot{j}}^{\mathbb{C}} = \sum\_{i=1}^{N\_p} \lambda\_{\dot{i}!} (\mathbf{r}\_i - \mathbf{r}\_{\dot{j}}) (1 - \delta\_{\dot{i}\dot{j}}) \tag{81}$$

where δij is the Kronecker delta symbol, defined as

$$
\delta\_{ij} = \begin{cases} 1 & \text{for} \quad i = j \\ 0 & \text{otherwise} \end{cases} \tag{82}
$$

Among Np � 1 pairs made by particle j (excluding itself), if particle j does 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 j by k is

$$F\_{j \gets k}^{\mathbb{C}} = -\nabla\_{\rangle} \Phi = \lambda\_{jk} (\mathbf{r}\_{j} - \mathbf{r}\_{k}) \tag{83}$$

and, similarly, the same force on k by j is

$$F\_{k \gets j}^{\mathbb{C}} = -\nabla\_k \Phi = \lambda\_{k\not\!j} (\mathbf{r}\_k - \mathbf{r}\_j) \tag{84}$$

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

$$F\_{j \gets k}^{\mathbb{C}} = -F\_{k \gets j}^{\mathbb{C}} \tag{85}$$

which follows Newton's third law, the action and reaction principle. 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 j can be expressed as

$$\mathbf{a}\_{\rangle} = \mathbf{a}\_{\rangle}^{\dagger} + \mathbf{a}\_{\rangle}^{\mathrm{C}} \tag{86}$$

where a† <sup>j</sup> is the unconstrained acceleration and aC <sup>j</sup> is the constrained acceleration, represented as

$$\mathbf{a}\_{\dot{j}}^{\mathcal{C}} = \frac{1}{m\_{\dot{j}}} \boldsymbol{\lambda}\_{\dot{j}k} (\mathbf{r}\_k - \mathbf{r}\_{\dot{j}}) \tag{87}$$

Assume particle j evolves from rj(t) at time t to r(t + δt) at time t + δt. Then, the future position at t + δt can be decomposed into two parts:

$$r\_{\dot{\jmath}}(t+\delta t) = r\_{\dot{\jmath}}^{\dagger}(t+\delta t) + \frac{1}{2}a\_{\dot{\jmath}}^{\mathbb{C}}\delta t^2\tag{88}$$

where

body. A constraint embedded between the two particles keeps the inter-particle distance

where dij is the fixed distance between particles i and j, usually an average of two contacting rigid bodies of diameters di and dj, i.e., dij = (di + dj)/2. This type of geometrical constraint is called holonomic. To develop an inter-particle interaction to satisfy Eq. (78), one defines the

i¼1

Exchanging positions of particles i and 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

X Np

j>i

i6¼j

<sup>2</sup> � <sup>d</sup><sup>2</sup>

ij ¼ 0 ð78Þ

λjiσij ð79Þ

λjiðr<sup>i</sup> � rjÞ ð80Þ

ð82Þ

λjiðr<sup>i</sup> � rjÞð1 � δijÞ ð81Þ

<sup>j</sup> <sup>k</sup> ¼ �∇j<sup>Φ</sup> <sup>¼</sup> <sup>λ</sup>jkðr<sup>j</sup> � <sup>r</sup>kÞ ð83<sup>Þ</sup>

<sup>k</sup> <sup>j</sup> ¼ �∇k<sup>Φ</sup> <sup>¼</sup> <sup>λ</sup>kjðr<sup>k</sup> � <sup>r</sup>jÞ ð84<sup>Þ</sup>

<sup>k</sup> <sup>j</sup> <sup>ð</sup>85<sup>Þ</sup>

<sup>2</sup> in the front of the summation is by convention.

σij ¼ ðr<sup>i</sup> � rjÞ

<sup>Φ</sup> <sup>¼</sup> <sup>1</sup> 2 X Np

exerted on particle j can be derived as a negative derivative of the constraint potential:

i¼1

�

<sup>δ</sup>ij <sup>¼</sup> 1 for <sup>i</sup> <sup>¼</sup> <sup>j</sup> 0 otherwise

Among Np � 1 pairs made by particle j (excluding itself), if particle j does not have any constraint to particle k, then λjk = 0, and symmetrically vice versa. For a two-body case, the

<sup>j</sup> ¼ �∇j<sup>Φ</sup> <sup>¼</sup> <sup>X</sup>

F <sup>C</sup>

where i runs from 1 to Np except for the case i = j (if so, λii = 0), or simply

F <sup>C</sup> <sup>j</sup> <sup>¼</sup> <sup>X</sup> Np

F <sup>C</sup>

F <sup>C</sup>

F <sup>C</sup>

<sup>j</sup> <sup>k</sup> ¼ �F<sup>C</sup>

which follows Newton's third law, the action and reaction principle. Summation over all constraint particles will give a zero resultant force. This is because the constraint force is an

invariant. A translational constraint between particle i and j is

constraint potential for all Np particles as

30 Granular Materials

where λji is a symmetric Lagrange's multiplier and <sup>1</sup>

where δij is the Kronecker delta symbol, defined as

holonomic constraint force acting on j by k is

and, similarly, the same force on k by j is

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

internal force and the sum of internal forces is zero.

$$\mathbf{r}\_{\rangle}^{\dagger}(t+\delta t) = \mathbf{r}\_{\rangle}(t) + \mathbf{v}\_{\rangle}(t)\delta t + \frac{1}{2}\mathbf{a}\_{\rangle}^{\dagger}\delta t^{2} \tag{89}$$

is the evolved position at time t + δt in the absence of the constraint force. A similar equation for particle k can 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:

$$
\lambda\_{i\dot{\eta}} \delta t^2 = \frac{\left[\Delta \mathbf{r}\_{i\dot{\eta}}^{\dagger}(t + \delta t)\right]^2 - d\_{i\dot{\eta}}^2}{\mu\_{i\dot{\eta}}^{-1} \left[\Delta \mathbf{r}\_{i\dot{\eta}}^{\dagger}(t + \delta t)\right] \cdot \left[\Delta \mathbf{r}\_{i\dot{\eta}}(t)\right]} \tag{90}
$$

where

$$
\Delta \mathbf{r}\_{ij}^{\dagger}(t + \delta t) = \mathbf{r}\_i^{\dagger}(t + \delta t) - \mathbf{r}\_j^{\dagger}(t + \delta t) \tag{91}
$$

$$
\Delta \mathbf{r}\_{\overline{i}\overline{j}}(t) = \mathbf{r}\_{\overline{i}}(t) - \mathbf{r}\_{\overline{j}}(t) \tag{92}
$$

and

$$
\mu\_{\vec{\eta}}^{-1} = m\_i^{-1} + m\_j^{-1} \tag{93}
$$

is called the reduced mass of particles i and j. In Eq. (90), λij can be determined using an iterative method.


This iterative procedure will continue until the Lagrange multiplier converges within a preset tolerable error for the position evolution from time t to t + δt. So far, the constraint force modifies the particle position at time t + δt from its unconstrained position r† <sup>j</sup> ð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

$$(\mathbf{r}\_i - \mathbf{r}\_j) \cdot (\mathbf{v}\_i - \mathbf{v}\_j) = \mathbf{0} \tag{94}$$

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

$$\mathbf{v}\_{\dot{j}}(t') = \mathbf{v}\_{\dot{j}}^{\dagger}(t') + \mathbf{a}\_{\dot{j}}^{\mathbb{C}} \delta t \tag{95}$$

where

$$\mathbf{v}\_{\rangle}^{\dagger}(t') = \mathbf{v}\_{\rangle}(t) + \mathbf{a}\_{\rangle}^{\dagger} \delta t \tag{96}$$

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

$$
\Delta \mathbf{r}\_{\vec{\eta}}(t + \delta t) \cdot \Delta \mathbf{v}\_{\vec{\eta}}^{\dagger} = -\Delta \mathbf{r}\_{\vec{\eta}}(t + \delta t) \cdot \Delta \mathbf{a}\_{\vec{\eta}}^{\mathbb{C}} \delta t \tag{97}
$$

where

$$
\Delta \mathbf{a}\_{i\dot{j}}^{\mathbb{C}} = \frac{\kappa\_{i\dot{j}}}{\mu\_{i\dot{j}}} (\mathbf{r}\_{\dot{j}} - \mathbf{r}\_{i}) \tag{98}
$$

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

$$
\kappa\_{\vec{\imath}\vec{\jmath}}\delta t = \frac{\Delta r\_{\vec{\imath}\vec{\jmath}}(t+\delta t) \cdot \Delta r\_{\vec{\imath}\vec{\jmath}}^{\dagger}}{d\_{\vec{\imath}\vec{\jmath}}^2 \mu\_{\vec{\imath}\vec{\jmath}}^{-1}}\tag{99}$$

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

1. Calculate the unconstrained velocity at the next time step, v† (t + δt) for particles i and j, and calculate their difference Δv† ij <sup>¼</sup> <sup>v</sup>† <sup>i</sup> � <sup>v</sup>† <sup>j</sup> at time t + δt.


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

This iterative procedure will continue until the Lagrange multiplier converges within a preset tolerable error for the position evolution from time t to t + δt. So far, the constraint force

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

valid both at time t and t + δt. Similar to the position evolution, the translational velocity at

Þ ¼ <sup>v</sup>jðtÞ þ <sup>a</sup>†

is the updated velocity from time t without the holonomic constraint. Substitution of Eq. (95)

is the constraint acceleration for velocity correction. Here, κij plays a similar role of λij, but it is

<sup>κ</sup>ijδ<sup>t</sup> <sup>¼</sup> <sup>Δ</sup>rijð<sup>t</sup> <sup>þ</sup> <sup>δ</sup>tÞ � <sup>Δ</sup>v†

1. Calculate the unconstrained velocity at the next time step, v† (t + δt) for particles i and j,

<sup>j</sup> at time t + δt.

<sup>i</sup> � <sup>v</sup>†

d2 ijμ�<sup>1</sup> ij

ij ¼ �Δrijð<sup>t</sup> <sup>þ</sup> <sup>δ</sup>tÞ � <sup>Δ</sup>a<sup>C</sup>

ij

modifies the particle position at time t + δt from its unconstrained position r†

vjðt 0 Þ ¼ <sup>v</sup>† j ðt 0 Þ þ <sup>a</sup><sup>C</sup>

v† j ðt 0

<sup>Δ</sup>rijð<sup>t</sup> <sup>þ</sup> <sup>δ</sup>tÞ � <sup>Δ</sup>v†

Δa<sup>C</sup> ij <sup>¼</sup> <sup>κ</sup>ij μij

independently determined only to update the velocity, which is calculated as

ij <sup>¼</sup> <sup>v</sup>†

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

and calculate their difference Δv†

velocity after the constrained evolution of position is the same as before the evolution.

<sup>j</sup> , go to step 2, unless λij converges to a finite value. Otherwise, store

ðr<sup>i</sup> � rjÞ�ðv<sup>i</sup> � vjÞ ¼ 0 ð94Þ

<sup>j</sup> δt ð95Þ

<sup>j</sup> δt ð96Þ

ðr<sup>j</sup> � riÞ ð98Þ

ij δt ð97Þ

ð99Þ

<sup>j</sup> ðt þ δtÞ, but the

at that time: r†

32 Granular Materials

gives

where

where

into Eq. (94) gives

5. Having the updated r†

time t + δt is represented as

<sup>j</sup> ðt þ δtÞ rjðt þ δtÞ.

information at t + δt, and go to the next time step.

4. Replace v<sup>j</sup> (t + δt) by v† <sup>j</sup> ð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<sup>3</sup> = 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 (H2O) 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.

Figure 4. A dimer of bi-dispersed spheres of holonomic constraint only.

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

$$\mathbf{p}\_k = m\_k \mathbf{v}\_k \tag{100}$$

$$H\_k = \mathfrak{J}\_k \mathfrak{o}\_k \tag{101}$$

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

$$M\_a \mathcal{V}\_a = m\_1 \mathfrak{v}\_1 + m\_2 \mathfrak{v}\_2 \tag{102}$$

where Ma = m<sup>1</sup> + m<sup>2</sup> is the total mass of the two-body aggregate and V<sup>a</sup> is the translational velocity of the center of mass of the aggregate. The total angular momentum has, however, a slightly different formulation:

$$H\_d = (\mathbf{J}\_1 + m\_2 \mathbf{d}\_1^2) \boldsymbol{\omega}\_1 + (\mathbf{J}\_2 + m\_2 \mathbf{d}\_2^2) \boldsymbol{\omega}\_2 = \mathbf{J}\_a \boldsymbol{\Omega}\_a \tag{103}$$

where

$$J\_a = I\_1 + m\_1 d\_1^2 + I\_2 + m\_2 d\_2^2 \tag{104}$$

is the mass moment inertia of the aggregate. Here, dk is the shortest distance between the center of mass of particle k and the rotation axis of the aggregate, passing through the center of mass of the aggregate, rCM:

$$d\_k = \frac{(\mathbf{r}\_k - \mathbf{r}\_{\rm CM}) \cdot \boldsymbol{\Omega}\_u}{|\boldsymbol{\Omega}\_u|} \tag{105}$$

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

$$\mathbf{M} = \sum\_{\mathbf{j}} (\mathbf{r}\_{\mathbf{j}} - \mathbf{R}\_{\mathbf{CM}}) \times (\mathbf{F}\_{\mathbf{j}}^{\mathrm{Ex}} + \mathbf{F}\_{\mathbf{j}}^{\mathrm{C}}) \tag{106}$$

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

$$\frac{dH\_d}{dt} = M\tag{107}$$

$$\mathbf{H}\_{\mathbf{a}}(\mathbf{t} + \delta \mathbf{t}) = \mathbf{H}\_{\mathbf{a}}(\mathbf{t}) + \mathbf{M} \delta \mathbf{t} \tag{108}$$

$$
\Omega \Omega\_{\mathfrak{a}}(t + \delta t) = \Omega\_{\mathfrak{a}}(t) + f\_{\mathfrak{a}}^{-1} \mathbf{M} \delta t \tag{109}
$$

Finally, all associated particles to the aggregate have the same angular velocity Ω<sup>a</sup> 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, Ja.

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

Consider a perfectly inelastic collision event between two approaching particles i and j, moving with translational and rotational velocities of (vi, ωi) and (vj, ωj), respectively, before their

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

where Ma = m<sup>1</sup> + m<sup>2</sup> is the total mass of the two-body aggregate and V<sup>a</sup> is the translational velocity of the center of mass of the aggregate. The total angular momentum has, however, a

<sup>1</sup>Þω<sup>1</sup> þ ðJ<sup>2</sup> <sup>þ</sup> <sup>m</sup>2d<sup>2</sup>

is the mass moment inertia of the aggregate. Here, dk is the shortest distance between the center of mass of particle k and the rotation axis of the aggregate, passing through the center of mass

dk <sup>¼</sup> <sup>ð</sup>r<sup>k</sup> � <sup>r</sup>CMÞ � <sup>Ω</sup><sup>a</sup>

After the perfectly inelastic collision, two compounded particles will have the same angular

<sup>ð</sup>r<sup>j</sup> � <sup>R</sup>CMÞ�ð<sup>F</sup> Ex

dH<sup>a</sup>

Ωaðt þ δtÞ ¼ ΩaðtÞ þ J

Finally, all associated particles to the aggregate have the same angular velocity Ω<sup>a</sup> at any time unless they break apart of slide into each other. It is assumed that the mutual compressions ξ

<sup>j</sup> <sup>þ</sup> <sup>F</sup> <sup>C</sup>

�1

<sup>1</sup> <sup>þ</sup> <sup>I</sup><sup>2</sup> <sup>þ</sup> <sup>m</sup>2d<sup>2</sup>

p<sup>k</sup> ¼ mkv<sup>k</sup> ð100Þ

H<sup>k</sup> ¼ Jkω<sup>k</sup> ð101Þ

<sup>2</sup>Þω<sup>2</sup> ¼ JaΩ<sup>a</sup> ð103Þ

<sup>2</sup> ð104Þ

<sup>j</sup> Þ ð106Þ

<sup>j</sup>Ωa<sup>j</sup> <sup>ð</sup>105<sup>Þ</sup>

dt <sup>¼</sup> <sup>M</sup> <sup>ð</sup>107<sup>Þ</sup>

<sup>a</sup> Mδt ð109Þ

Haðt þ δtÞ ¼ HaðtÞ þ Mδt ð108Þ

MaV<sup>a</sup> ¼ m1v<sup>1</sup> þ m2v<sup>2</sup> ð102Þ

collision. For k = i, j, the linear and angular momenta are

<sup>H</sup><sup>a</sup> ¼ ðJ<sup>1</sup> <sup>þ</sup> <sup>m</sup>2d<sup>2</sup>

Ja <sup>¼</sup> <sup>I</sup><sup>1</sup> <sup>þ</sup> <sup>m</sup>1d<sup>2</sup>

velocity Ωa. The total torque acting on the aggregate is rigorously represented as

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

<sup>M</sup> <sup>¼</sup> <sup>X</sup> j

slightly different formulation:

of the aggregate, rCM:

where

34 Granular Materials

Granular dynamics simulations in a fluid medium must be an open problem in state-of-theart computational research. Since the length scale of the forces acting on touching granules

Figure 5. A spherical trimer of (a) holonomic constraint only and (b) holonomic and non-holonomic constraints. In each case, the center-to-center distances between three pairs of particles are fixed.

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 i from j has the same magnitude and opposite direction to that from particle j from 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.

## 4. Concluding remarks

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.

## Acknowledgements

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 i from j has the same magnitude and opposite direction to that from particle j from 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

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

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 signif-

questionable to predict that the granular state will be the fifth matter phase.

doubled.

36 Granular Materials

4. Concluding remarks

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

## Author details

Albert S. Kim<sup>1</sup> \* and Hyeon-Ju Kim2

\*Address all correspondence to: albertsk@hawaii.edu

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

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

## References


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

[9] Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. Journal

[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

[13] Einstein A. Zur Theorie der Brownschen Bewegung. Annalen der Physik. 1906;324

[14] Langevin P. Sur la théorie du mouvement brownien. Comptes-rendus de I'Académie des

[15] Chandrasekhar S. Stochastic problems in physics and astronomy. Reviews of Modern

[16] Chandrasekhar S. Brownian motion, dynamical friction, and stellar dynamics. Reviews of

[17] Kubo R. The fluctuation-dissipation theorem. Reports on Progress in Physics. 1966;29:

[18] Stokes GG. On the effect of internal friction of fluids on the motion of pendulums. Trans-

[19] Ermakova L, Sidorova M, Savina I, Aleksandrov D. Structural and electrochemical parameters of asymmetric membranes for reverse osmosis. Colloids and Surfaces A.

[20] Ermak DL, Buckholz H. Numerical integration of the Langevin equation: Monte Carlo

[21] Ermak DL, Mccammon JA. Brownian dynamics with hydrodynamic interactions. Journal

[22] Ermak D. A computer simulation of charged particles in solution. I. Technique and

[23] Ermak DL. A computer simulation of charged particles in solution. I. Technique and

[24] Chen JC, Elimelech M, Kim AS. Monte Carlo simulation of colloidal membrane filtration: Model development with application to characterization of colloid phase transition.

equilibrium properties. Journal of Chemical Physics. 1975;62(10):4189–4196

equilibrium properties. Journal of Chemical Physics. 1975;62(10):4189–4196

actions of the Cambridge Philosophical Society. 1851;9:1–106

simulation. Journal of Computational Physics. 1980;35:169–182

General Existence of Active Molecules in Organic and Inorganic Bodies

of Computational Chemistry. 1983;4(2):187–217

(2):371–381

38 Granular Materials

255–284

Sciences (Paris). 1908;146:530–533

Modern Physics. 1949;21(3):383–388

of Chemical Physics. 1978;69(4):1352–1360

Journal of Membrane Science. 2005;255(1):291–305

Physics. 1943;15(1):1–89

1998;142(2–3):265–274


[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

[41] Brady JF, Phillips RJ, Lester JC, Bossis G. Dynamic simulation of hydrodynamically

[42] Wang M, Brady JF. Spectral ewald acceleration of stokesian dynamics for polydisperse

[43] Sierou A, Brady JF. Accelerated stokesian dynamics simulations. Journal of Fluid Mechanics.

[44] Kim AS. Dissipative hydrodynamics of rigid spherical particles. Chemistry Letters.

[46] Ito M. An extension of nonlinear evolution equations of the K-dV (mK-dV) type to higher

[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

[49] Kim AS, Kang ST. Microhydrodynamics simulation of Single-collector granular filtration.

[50] Hoomans BPB. Granular Dynamics of Gas-Solid Two-Phase Flows. Universiteit Twente;

[51] Poschel T, Schwager T. Computational Granular Dynamics, Models and Algorithms.

[52] Verwey EJ, Overbeek JTG. Theory of the Stability of Lyophobic Colloids. Amsterdam:

[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

[54] Kim AS. Constraint dissipative hydrodynamics (HydroRattle) algorithm for aggregate

[55] Hertz H. Ueber die Verdunstung der Flussigkeiten, insbesondere des Quecksilbers, im

[56] Brilliantov NV, Spahn F, Hertzsch JM, Poschel T. Model for collisions in granular gases.

[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

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

orders. Journal of the Physical Society of Japan. 1980;49:771–778

settling. Journal of Colloid and Interface Science. 2004;271(1):110–119

Chemistry Letters. 2012;41(10):1288–1290

dynamics. Chemistry Letters. 2012;41(10):1285–1287

Physical Review E. 1996;53(5):5382–5392

luftleeren Raume. Annalen der Physik. 1882;253(10):177–193

viscoelastic spheres. Physical Review E. 1999;60(4):4465–4472

Enschede, Netherlands; 2000

Springer, Berlin; 2005

USSR. 1941;14:633–662

Elsevier; 1948

interacting suspensions. Journal of Fluid Mechanics. 1988;195:257–280

suspensions. Journal of Computational Physics. 2016;306:443–477

2001;448

40 Granular Materials

2012;41(10):1128–1130

