Open access peer-reviewed chapter
By Albert S. Kim and Hyeon-Ju Kim
Submitted: November 7th 2016Reviewed: April 13th 2017Published: September 6th 2017
DOI: 10.5772/intechopen.69196
Downloaded: 591
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.
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), 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 for the motion of bodies can be summarized as follows:
Newton’s first law states that an object will remain at rest or in linear motion unless acted upon by an external force. (This first law describes a constant velocity motion in the absence of an external force, which is a special case of zero acceleration in the second law below.)
Newton’s second law states that the net (unbalanced) force F on an object is equal to the product of the mass m and acceleration a of the object:
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.)
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
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:
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
respectively, where r and θ are the linear and angular positions of the object, of which time derivatives are the linear and angular velocities, respectively:
Where a and α have three components each so that [a, α]T 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
The classical mechanics problems are usually reduced to solve for q(t) and
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:
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
Then, Eq. (9) can be rewritten as
where
where V is called the potential energy function. Then, the force integration is simply
Substitution of Eq. (13) into Eq. (9) gives
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 x1 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:
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.
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
At time t1 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
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
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
where δq(t), called a variation of q(t), is a function, assumed to be small everywhere in the interval of time from t1 to t2. Since Eq. (18) must include the values of q(t1) and q(t2), we can now conclude that
In this case, the change in S when q is replaced by q + δq is equal to
We expand the difference in powers of δq and
or equivalently
Since
The condition of Eq. (19) shows that the integrated term in Eq. (23) is zero:
and the remaining integral is
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
When the system has more than one degree of freedom, then Eq. (26) becomes
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
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.
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.
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 δt4 order. The particle position at time t + δt can be expanded at time t using Taylor’s series:
and by replacing +δt with +δt, one obtains
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
and
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.
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:
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′(t) is the random fluctuation force caused by adjacent, rapidly-moving solvent molecules. It is assumed that f′ is a Gaussian process with infinitely small correlation time, which renders the famous fluctuation-dissipation theorem [17]
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
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
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].
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:
where FP is a conservative inter-particle force, FD is a dissipative force, and FR is a random fluctuation force. It is required that FD and FR are linear and independent of the momentum, respectively. A simple form of these force are hypothesized as
where rij = |ri − rj|, vij = vi − vj, eij = (ri − rj)/rij, and
where dWij (= dWji) is independent increment and of the Wiener process, satisfying
i.e., dWij(t) is an infinitesimal of order 1/2 (i.e., proportional to
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
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.
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
where the superscript ∞ indicates that the near-field lubrication forces are excluded. The grand resistance matrix is formed by inverting the grand mobility matrix and correcting the near-field lubrication such as
where
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].
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
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, Qp 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 ⋅ B⟩ = 2kBT, where kB is the Boltzmann constant, and dW is the Ito-Wiener process of 6Np elements [45, 46]. The generalized force/torque is defined, similar to Eq. (42), as
where
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:
Collision of sphere i and j.
where vij = vi − vj is the relative velocity of the center of mass of particle i to j of velocity vi 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
and
respectively. The coefficients of restitution in normal and tangential direction, ϵn and ϵt, are defined as
where 0 ≤ ϵn ≤ 1 and −1 ≤ ϵnt ≤ 1. The primed and unprimed variables indicate the pre- and post-collision quantities, respectively. Then, the velocities of particles i and j after the collision can be represented as functions of ϵn, ϵt,
where for k = i, j. Here, Jk is the mass moment of inertia and
One can calculate the energy loss of particle k,
Table 1 summarizes specific features of the above-discussed simulation methods. Acceleration can be included by all methods in principle, but is barely employed in BD and SD. This is mainly due to the time interval restriction. Instead, the fluctuation-dissipation theorem is intrinsically embedded and satisfied in BD, SD, DPD, and DHD formalisms. All the simulation methods can surely include effects of conservative, external force fields. Brownian motion and multi-body hydrodynamic forces are (most) accurately implemented in SD and DHD, but only hydrodynamics is important in the many-body motion of non-Brownian granules. The DHD is the only simulation method that can mimic many-body hydrodynamics without the time interval restriction. Constraint forces/torques can be easily applied to any dynamics methods to simulate compound particles or aggregates. SD has the same capability but the collision rules are not included. Due to typical ranges of particle sizes in the simulation methods, the inelastic collision (or contact) forces are included only in DEM. Since DPD allows particle deformation by using the relaxed, pairwise hydrodynamic resistance between two particles, i.e., ωD and ωR in Eq. (44), employing the collision force of rigid bodies in DPD must be out of its original scope. The current state-of-the-art DEM algorithms can be further improved to mimic complex phenomena of dry granules of various shapes [50, 51]. However, in our opinion, granular motion in a liquid phase can be accurately simulated by only using DHD or SD. To accurately simulate a large number of granules of various sizes in complex fluid environments, multi-body hydrodynamics must be rigorously implemented in the current DEM method.
| 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 |
| 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 | |
| Time interval dt | Any | δt ≫ τ | δt ≫ τ | Any | Any | Any | |
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 τ.
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.
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:
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.,
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
where
where H (x) is the Heaviside step function, defined as
In three-dimensional space, vector quantities between particles i and j can be decomposed into the normal and tangential directions. Using the mathematical identity
one can replace A and B by n and C by Fij to write
From Eq. (67), the normal and tangential force components can be expressed as
The contact force between elastic spheres was originally developed by Hertz [55] and was later generalized for viscoelastic (damped) particles [56, 57] as
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.,
for particles of different Y and ν values. Assuming the dissipative constant A is also particle-specific and additive, the most general expression of the normal force between two viscoelastic spheres may be [51]
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
This indicates that ϵn 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
which provides the tangential force of
where γt is a fitting coefficient, proportional to the tangential dissipative force in magnitude. In Eq. (75), the shear force is limited by Coulomb’s friction law of |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]:
where ζ is the relative tangential shift, ζ0 is its macroscopic maximum value, and ⌊x⌋ denotes the integer of x. Typical values of ζ0 / aeff range from 10−7 to 10−3. Similar to Eq. (73), the tangential coefficient of restitution can be calculated as
If the friction coefficient μ is known, experimental measurement of ϵt allows us to inversely calculate ζ0. Techniques to determine the coefficients of restitution of colliding viscoelastic spheres can be found elsewhere [58]. This approach can be readily applied to densely packed granular medium with small movements or vibrations such as the Brazilian bean problem [59–64]. However, time integration per collision, consisting of compression and restoration, should be from 0 to 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.
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 particle contacts, followed by inelastic overlaps in a series of collision events. This is because the lubrication force is logarithmically proportional to the surface-to-surface distance between two particles and thus it diverges when two particles start inelastically overlapping. Second, the order-disorder phase transition includes the changes in the flow structure between ordered and disordered phases, which yields increases in the drag force between particles [70]. These ordered states are different from fixed 3D structures found in solids, but similar to amorphous aggregates. Third, the dilatancy mechanism describes the shear-thickening such that the effective volume of particle packing increases under the shear. When particles are confined in local spaces or partially jammed, the shear force pushes particles toward the containing walls. Additional stress can be developed on the walls and backward responses may generate extra stress between particles and walls [71] in contact interfaces. Recently, discontinuous shear thickening (DST) was proposed to explain the shear thickening experiments [69], of which comprehensive review can be found elsewhere [72, 73]. To the best of our knowledge, the fundamental mechanism of the shear thickening has not been fully discovered.
Sample simulation: Figure 2 shows results of a sample simulation using the mechanisms described in Section 3.1.1, which can be further extended to DST simulations. The gray spheres are loosely packed, ideally forming a body-centered cubic structure. The top layer consists of 5 × 5 = 25 regularly packed spheres, under which 6 × 6 = 36 spheres are located. These spheres are closely located to each other but not touching, having very small surface-to-surface distance between the nearest neighbors. The two layers are packed five times vertically, and therefore the total number of gray particles is 366. This preliminary simulation aims to see the impact responses of the small packed spheres when hit by a big, fast intruder. The intruder particle initially moves with the downward velocity of vz = 13.00 × 10−3 m/s and spins with angular velocity of ωz = 0.56 rad/s.
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 A was set at 2.5 × 10−5, for both spheres and walls [56].
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-to-many collisions between the intruder and dissipating granules. The intruder’s impact is transmitted through dynamic chains of contacting spheres, but not all of the packed granules participate in the force transmission. There are same granules almost free from the intruder’s impact.
Calculation of forces and torques exerted on each particle and their visualization can help understand the transient behavior and design dynamic granular materials. Granules initially located near the container walls have much less spaces to move. The kinetic energy of the intruder is transmitted to these granules at the cul-de-sac and mostly dissipated on the wall surfaces. Returning force to the intruder is initiated from the boundary. Applications of this simulation method can be designed to include a large number of real applications for characterization and prediction of transient granular material properties. One important key issue is, as indicated above, to identify and visualize the transmission chains of force/torque, which dynamically form and disappear. Figure 3 shows the initial impact event when the intruder starts penetrating the loose granular packing, visualized using open visualization tool (OVITO) [77]. The top and bottom rows show the top and side views of the intruder collision. The left column shows particle positions as well as force vectors, and the right column shows only particle configurations where color indicates the magnitude of the net force. Impact from the intruder is somehow irregularly distributed around the top granules. On the left-column images, arrows and their colors indicate force directions and magnitudes. A number of downward vectors indicate that the gravitational force is dominant for non-touching granules. The upward arrow in the intruder implies that the contact force to the intruder, which is similar to the normal force developed on the interface between an object and a wall, exceeds the gravitational force exerted on the intruder. Even though the force direction is temporarily inverted, the intruder still goes down due to the inertia of the high initial velocity. Figure 3 clearly shows that only a partial number of granules participate in the force transmission from the intruder to the packed granules. These force chains are very transient, and more importantly, the magnitude of the transmitted force diminishes as time goes by due to the intrinsic inelastic nature of the granules. For the future, designs of smart damping materials can be achieved not only by understanding various mechanical properties of the granules, but also by controlling specific initial and boundary condition, which leads inelastic many-granule systems to behave as smart materials.
Force chain visualization of the intruder using open visualization tool (OVITO).
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
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
where λji is a symmetric Lagrange’s multiplier and
where i runs from 1 to Np except for the case i = j (if so, λii = 0), or simply
where δij is the Kronecker delta symbol, defined as
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
and, similarly, the same force on k by j is
Since λjk is symmetric (λjk = λkj), one can show that
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
where
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:
where
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 (δt4) and higher, the representation of the Lagrange multiplier:
where
and
is called the reduced mass of particles i and j. In Eq. (90), λij can be determined using an iterative method.
Initially, the unconstrained position
Insert
Calculate the constraint acceleration of particle j,
Update the particle position using Eq. (88), which is under the influence of constrained and unconstrained accelerations. Set this updated position as the unconstrained position at that time:
Having the updated
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
Velocity evolution: Differentiation of the holonomic constraint Eq. (78) with respect to time gives
valid both at time t and t + δt. Similar to the position evolution, the translational velocity at time t + δt is represented as
where
is the updated velocity from time t without the holonomic constraint. Substitution of Eq. (95) into Eq. (94) gives
where
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
To update κij, a similar iteration method can be used:
Calculate the unconstrained velocity at the next time step, v† (t + δt) for particles i and j, and calculate their difference
Calculate κij using Eq. (99) using previously determined ri(t + δt) and rj(t + δt).
Update
Replace vj (t + δt) by
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.53 = 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.
A dimer of bi-dispersed spheres of holonomic constraint only.
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.
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
respectively. After the two particles are permanently attached, the total linear momentum is
where Ma = m1 + m2 is the total mass of the two-body aggregate and Va is the translational velocity of the center of mass of the aggregate. The total angular momentum has, however, a slightly different formulation:
where
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:
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
and the time evolution equations related to the total angular momentum are
Finally, all associated particles to the aggregate have the same angular velocity Ωa at any time unless they break apart of slide into each other. It is assumed that the mutual compressions ξ between the associated particles are small enough to be neglected in calculating the aggregate’s mass moment of inertia, 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.
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.
Granular dynamics simulations in a fluid medium must be an open problem in state-of-the-art computational research. Since the length scale of the forces acting on touching granules is much smaller than the granular sizes, simulations seem to be very efficient if parallel algorithms are adequately used. Domain decomposition scheme is one of the most widely used parallel algorithms. The system is divided into several small sub-domains, and a dynamic simulation in a spatial sub-domain is conducted by an individual computing unit, such as a core. This decomposition method is very efficient if particles are evenly distributed in space, which is usual in equilibrium simulations of MD and Monte Carlo. As granules are highly subjected to the gravitational force and hydrodynamic interactions, spatial biases are almost inevitable in both their transient and stationary states. Granular dynamics has at least three length scales of different orders of magnitude. The mutual compression distance, i.e., inter-particle overlap distance, is at least three or four orders of magnitude smaller than the particle size. In a fluid medium, hydrodynamic interactions are long-ranged and quite significant when the surface-to-surface distance between nearest neighbors is about the particle diameter. Motion of heavy granules in a fluid flow may distort the ambient flow-field as well as the hydrodynamic forces exerted on adjacent particles. In granular dynamics, when the contact force exerted between two colliding particles, the force acting on particle 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.
In nature, phases of matter have been conventionally believed to be those of gas, liquid, and solid, in which specific phase transformations are possible between the states. A recent addition of the plasma state has increased the total number of material states from three to four. It is questionable to predict that the granular state will be the fifth matter phase.
On the other hand, the granules dynamically change the representative phase based on external influences. Stationary and compressed granules behave similar to amorphous solids or gel materials. Moving like a fluid, mud or sediments create their own pathways by minimizing hydrodynamic influences. A fast flow with granular materials, such as in streams and ocean, creates a dense turbidity flow, but a decelerating flow field initiates a granular phase change from a flowing liquid to a packed solid. Small dry granules behave similar to dust in the wind, for which standard gas transport theory can predict dynamic behaviors of granular gases. In our opinion, granules are chameleon materials, transitioning their phases dynamically. No equilibrium exists in a granular phase so that thermodynamic fluctuation of the granular state cannot provide material or phase constants. As granular dynamics in transient force/torque fields significantly removes steady-state behaviors, the initial and boundary conditions become extraordinarily important in analyzing many-body granular motion. To utilize specific behavior of granules, granules can be dynamically controlled in vibrating, oscillating, or swinging phases. The time-correlation scale of dissipating granules is not short enough to use the Markovian chain concept because granular paths in both the real and phase spaces significantly influence their fates.
This work was financially supported by the National R&D project of “Development of 1MW Ocean Thermal Energy Conversion Plant for Demonstration” (PMS3680) from the Korean Ministry of Oceans and Fisheries and used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575.
591total chapter downloads
1Crossref citations
Login to your personal dashboard for more detailed statistics on your publications.
Access personal reportingBy Hao Shi, Dalila Vescovi, Abhinendra Singh, Sudeshna Roy, Vanessa Magnanimo and Stefan Luding
Edited by Michael Sakellariou
We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.
More About Us