Impact-Induced Hardening in Dense Frictional Suspensions

By employing the lattice Boltzmann method, we perform simulations of dense suspensions under impact that incorporate the contact between suspended particles as well as the free surface of the suspensions. Our simulation for a free falling impactor on dense suspensions semi-quantitatively reproduces experimental results, where the impactor rebounds for high speed impact and high volume fraction shortly after the impact before subsequently sinking. We observe that the response depends on the radius of the impactor, which leads to fitting our simulation data to a phenomenological model based on the Hertzian contact theory. When the rebound takes place, a localized jammed region is formed by the frictional contacts between suspended particles. Furthermore, persistent homology analysis is used to elucidate the significance of the topological structure of the force chains, where the total persistence of connected components correlates with the force supporting the impactor.

By employing the lattice Boltzmann method, we perform simulations of dense suspensions under impact that incorporate the contact between suspended particles as well as the free surface of the suspensions. Our simulation for a free falling impactor on dense suspensions semi-quantitatively reproduces experimental results, where the impactor rebounds for high speed impact and high volume fraction shortly after the impact before subsequently sinking. We observe that the response depends on the radius of the impactor, which leads to fitting our simulation data to a phenomenological model based on the Hertzian contact theory. When the rebound takes place, a localized jammed region is formed by the frictional contacts between suspended particles. Furthermore, persistent homology analysis is used to elucidate the significance of the topological structure of the force chains, where the total persistence of connected components correlates with the force supporting the impactor.
Introduction.-A typical example of the non-Newtonian behavior of dense suspensions is the fact that water containing cornstarch behaves like a fluid, while running persons can stay afloat on top of it. The phenomenon where the fluid exhibits solid-like response under fast impact has also gained attention for practical applications, such as protective vests. Such responses are often associated with both the continuous (CST) and discontinuous (DST) shear thickening in dense suspensions under simple shear flow. Nevertheless, the simple shear flow differs from running on the fluid. The latter should be treated as an impact problem in which the dependence on the impact speed and the surface effect are important. Some efforts have been made to reproduce the impact behavior through experiments [1][2][3][4][5][6]. Waitukaitis and Jaeger indicated that the added mass effect in solidification after an impact is important [1]. The propagation of the dynamic jamming-front is also observed in a quasi-two-dimensional experiment for impact on a suspension [2]. Moreover, a series of constant speed penetration experiments [3][4][5] suggested that the added mass effect alone is insufficient in explaining the hugh stress exerted on the impactor. Instead, they showed that the stress on the impactor increases as it approaches the bottom boundary. The importance of the boundary characterized by depth of the suspension H has been suggested by a free-falling impactor experiment [6], where rebounding motion of the impactor was observed shortly after the impact. In order to extract microscopic information of two-dimensional dry granular materials in experiments, one approach is to use photoelastic disks to visualize the force acting on each grain [7]. Similarly, numerical simulation is an important tool to understand the microscopic mechanism behind exotic phenomena in suspensions since the motion of the suspended particles is not visible in three-dimensional experiments. One of the remarkable insights from simulations is that the frictional contacts between particles are important for ob-serving DST under simple shear [8][9][10][11][12]. Nevertheless, a particle-based simulation of a free-falling projectile hitting a suspension has not been reported so far because of the difficulty of simulating suspension with free surface. As far as we know, the first fluid-based simulation of suspensions under impact has been conducted recently in Ref. [13], where the authors successfully reproduced various interesting processes for suspensions under impact, such as the viscoelastic response of a dense suspension to a rotating wheel. Since, however, their fluid simulations with a constitutive equation cannot capture the particle dynamics, the mechanism behind impact-induced hardening on the microscopic level remains elusive. In this Letter, we utilize the lattice Boltzmann method (LBM) [14][15][16] to perform a simulation that can incorporate particle dynamics with the free surface of the suspensions [17][18][19].
Simulation method.-The hydrodynamic fields are calculated on nodes inside the cells of a fixed Cartesian grid. Due to this discrete nature of the LBM, we take the lattice unit ∆x = 0.2a min (a min is the radius of the smallest particle) for the hydrodynamic fields calculations. To simulate the free surface of the fluid, it is necessary to introduce interface nodes between the fluid and gas nodes, where we calculate the fluid density in a single cell m f and the liquid fraction λ [17][18][19]. To maintain a smooth surface and conserve m f , the excess m f is distributed equally among the interface nodes during time evolution. The detailed explanation on how to calculate the hydrodynamic fields, and the evolution equation for m f are presented in Ref. [20].
The force and the torque on particle i are given by Here, u i , ω i , m i , and I i are the translational velocity, angular velocity, mass, and the moment of inertia of particle i, respectively. F g i = −m i gẑ is the gravitational force acting on the suspended particles, where g is the gravitational acceleration andẑ is the unit vector in the vertical direction. Note that our LBM accounts for both the short and the long-range parts for the hydrodynamic force F h i and torque T h i [12,20,21], though some previous simulations contains only the short-range force [9,10]. The contact forces F c ij and torques T c ij are computed using the linear-dashpot model with Coulomb friction rules and friction coefficient µ [22]. We mainly use µ = 1 for our simulations but will discuss µ−dependence later. Finally, we also introduce the electrostatic repulsive force F r ij arising from the double layer to prevent particles from clustering [10,12].
About 600 suspended particles (bidisperse suspension with bidispersity ratio a max = 1.2a min where the radii of the large and small particles are a max and a min , respectively) are confined into a rectangular box (W × D × H) where we adopt the width W = D = 24a min , and height H = 12a min . We adopt the reflection rule as the the boundary condition on the walls. The impactor is a solid spherical object with radius a I = 3a min , and density ρ I = 4ρ f , where ρ f is the density of the suspended particles and the solvent. The force and torque on the impactor are given by F I = F h I +F c I +F g I and T I = T h I +T c I , respectively. F g I = −m I gẑ is the gravitational force acting on the impactor with mass m I . The contact force F c I is the force acting on the suspended particles, while the hydrodynamic force F h I is calculated using the bounceback rules [14,15]. In order to maintain numerical stability, we ignore the lubrication torque on the impactor (approximating the surface of the impactor as planar) as in Ref. [21]. However, the torque arising from the contact force is still considered explicitly. The impactor is released from various heights H 0 that correspond to the impact velocity as u I 0 = √ 2gH 0 , which also specifies the units of time in our simulation τ = a min /2g, units velocity u * = √ 2ga min , units of force F 0 = 4 3 πρ f a 3 min g, and units of stress σ 0 = F 0 /a 2 min . Impact induced hardening.-We plot the impactor speeds u I z (t)/u * against time for various volume fractions φ in Fig. 1(a), where we set the time t = 0 and height z = 0 at the moment of impact. One can observe the rebound of the impactor (u I z (t)/u * < 0) when φ ≥ 0.54 which agrees semi-quantitatively with the freefalling projectile experiment [6]. This rebound occurs due to the hardening of the suspension shortly after the impact. After the rebound, the suspension becomes soft, and the impactor sinks with a constant speed. Note that φ for rebounds is a little higher than that in the experiment [6] since rolling friction is ignored in our simulations [23]. When the impactor bounces, one can observe a peak in F I,z shown in Fig. 1(b), where we plot the total force exerted on the impactor for φ = 0.57 and u I 0,z /u * = 4.2 (the blue lines in Figs. 1(a) and 1(b)). Furthermore, we can see a peak in the contact force, which follows the peak by weak hydrodynamic contribution. The time difference between these two peaks is so small that they appear as a single peak in the total force. In an experiment with rod impactor [1], two peaks in the acceleration of the impactor is observed for deep suspensions. This is in contrast for the shallower suspension, in which rebound takes place, where the separation between the peaks is not detectable. Moreover, they also observed the second peak when the impact force is transmitted to the boundary. Here, the peak of the contact forces between particles in our simulation corresponds to the second peak in Ref. [1]. This indicates that the hardening takes place when the contact force network percolates from the impactor to the boundaries [24]. To clarify the importance of contact friction further, we simulate various friction coefficients µ. In Fig. 1(c), one can observe that the bouncing motion is weakened as µ decreases. This stands as the second proof that frictional contacts between particles are necessary to induce the resistance that makes the impactor rebounds. This µ-dependence is analogous to that for DST in dense suspensions under steady shear [8][9][10][11][12] and for impact in dry granular materials [25]. Figure 2(a) shows a clear dependence on the impactor radius a I . Therefore, we try to characterize the elasticity of suspensions by the viscoelastic Hertzian law for contact [26,27]. The equation motion for the deformation h for the Hertzian contact is written as where A and B are fitting parameters which correspond to the elastic modulus and viscosity, respectively. In Fig.  2(a), we plot the results of the simulation alongside with the solutions of Eq. (1). Although our model treats the suspension as a solid body, one can see that the results of the simulation agree with the model shortly after the impact and deviate afterwards. Therefore, this does not contradict the observation in Ref. [3], where the suspension has plastic response under impact since their impactor penetrates to the bottom boundary. The experimental observation that it is possible to run on suspensions are successfully reproduced qualitatively through our simulations, where we also observe that the impact-induced hardening depends on the impact speed u I 0,z , as shown in the diagram of Fig. 2(b). Although not observed clearly in the experiment of Ref. [6], the speed dependence also exists in the simulation of a rotating wheel driving over the surface of suspension [13]. Note that the highest rebound volume fraction φ b,max = 0.57 in our simulation is still below the frictional (µ = 1) jamming fraction φ µ=1 J ≈ 0.585 [10], whereas rebound takes place for φ b,min ≤ φ ≤ φ b,max , where φ b,min = 0.51. This range is similar to the observed volume fractions for both CST and DST under simple shear in numerical simulations [9,10,12]. However, one should recognize that two processes are different since impact-induced hardening is a heterogeneous and transient process, while shear thickening is a homogenous steady states process.
There is another difference between two processes as mentioned in Ref. [20] in which the shear stress is as large as the normal stress in DST while the shear stress is much smaller than the normal stress in impact processes. To understand the stress response of the suspension, we visualize the stress σ zz on each suspended particle in a slice through the middle of suspension immediately after the impact in Fig. 3(a). Here we observe a localized stress with distinctively high value of σ zz , which extends from the impactor to the boundary. On the other hand, we observe a uniformly weak magnitude of the shear stress [20]. Note that this localized response of the normal stress clearly distinguishes be-tween the phenomena of impact-induced hardening and DST. To see the origin of this mechanism on the microscopic level, we visualize the magnitude of the normal contact forces between particles scaled by the gravitational force |F c,n ij |/F 0 immediately after the impact in Fig. 3(b). The percolating force chains span from the impactor to the boundary. The regions of large σ zz and the force chains correspond to the dynamically jammed region in Fig. 3(a) and Refs. [1,[3][4][5], where we also observe a spanning region of high particle displacement ∆z [20]. As indicated in Refs. [2,3,13], the propagation speed of the jamming front depends on the impact speed. One can imagine that the vanishing of the stress exerted on the suspension by the impactor allows the suspension to relax, thus becoming soft, which in turn results in the impactor sinking. Persistent homology.-To elucidate the role of the force chains in impact-induced hardening, we analyze its topological structure by persistent homology [28]. In addition to successfully distinguishing the liquid, amorphous, and crystalline states of e.g. silicon dioxide [29], persistent homology allows us to quantify the structure of the force chains in granular materials [30,31], and in dense suspensions under simple shear [32]. Since no visible loops or higher dimensional structures are observed in the force network in Fig. 3(b), the relevant topological structure is only the connected component represented by the zeroth Betti number β 0 . The idea of persistent homology is to filter the force chains by increasing threshold θ f , where a link in a force chain appears when |F c,n ij |/F 0 ≤ θ f . We regard this as the birth of a connected component. As the threshold further increases, the structure grows in size as additional contacts are added. When connected components merge, the structure that is born later in the filtration (which has higher birth θ f ) dies. We record the birth θ f as θ f,b and the death θ f as θ f,d . This rule ensures that θ f,d > θ f,b . The algorithm for filtering chains is available in public domain [33,34]. Note that in Refs. [30][31][32], θ f,b is always larger than θ f,d , since they adopt filtration by reducing the threshold. We plot θ f,d against θ f,b for all connected components appearing in Fig. 3(b) in the persistence diagram ( Fig. 4(a)). The time evolution of the force chains and the persistence diagram can also be seen in Ref. [35]. Shortly after the impact, we observe more points far from the diagonal, representing the connected components which persist through the increments of the force threshold with the life span (θ f,d −θ f,b ). Intuitively, the contact force between particles cannot change abruptly. Therefore the only possible mechanism for the occurence of a long lifespan for some connected components is by forming a long chain. This argument shows that percolated force chains exists. A component with θ f,d = −1 has infinite persistence, i.e. it does not die until the filtration ends. Persistent homology ignores the effect of such contact forces.
The total persistence of the connected components T P 0 is the sum of all life spans in the persistence diagram, P D 0 , This allows us to describe the persistence diagram by a single number. Higher T P 0 means merging of longer percolated force chains, while T P 0 = 0 means that no the connected components are merged. Note that we ignore the components with infinite persistence since we are only interested in extracting the structural information. We plot T P 0 scaled by the number of suspended particles N against time in Fig. 4(b). It is remarkable that T P 0 reaches its peak at the same time as the corresponding contact force and that the shape of T P 0 is similar to that of the contact force. When the suspension becomes soft and the impactor sinks, we found T P 0 = 0, even though there are still contributions from the contact forces. These contributions come from not percolated force chains near the impactor. In other words, the peak of the contact force inducing the hardening of the suspension originates from the existence of percolated force chains. Our results support the findings in Refs. [1,[3][4][5][6] that indicate the importance of the boundaries in sustaining percolated force chains. On the other hand, for sheared suspensions, the total persistence of loops is more significant as it behaves similarly to the viscosity [32]. This distinction exists because for sheared suspension, the force chains are more structured and uniformly distributed than for suspension undergoing impact.
Conclusions and outlook.-We have reproduced the impact-induced hardening and the dynamically jammed region in dense suspensions through particle-based simulations. With the aid of persistence homology, we have confirmed the existence of percolated force chains when impact-induced hardening takes place.
Due to limitation of our computational resources, we have not analyzed the system size dependence. A sinking impactor in dense suspensions shows a distinct behavior, as it oscillates and has a stop-go cycle near the bottom [36]. Our simulation method and analysis can also be applied to the impactor rebound within the framework of the jamming transition as in Ref. [37], where the concepts of shear jamming and the DST are unified. Observation of a universal scaling law for impacts on dry granular media was recently reported [38]. It would be interesting to investigate whether such a scaling law also exists in the case of suspensions. These are targets of our next research.

Supplemental Material for "Impact-Induced Hardening on Dense Frictional
Suspensions" In this section, we review the outline of the lattice Boltzmann method (LBM) based on Refs. [S1-S6]. Due to the discrete nature of the LBM, one needs to discretize the unit of length into the lattice unit ∆x. In LBM, the hydrodynamic fields (density ρ f and velocity u f ) are calculated on nodes r inside cells of a fixed Cartesian grid as

Pradipto and Hisao Hayakawa
where c q is the lattice velocity of the direction q, and ∆c 3 is the volume element in the velocity space with ∆c = ∆x/∆t. f q (r) is the abbreviation of f q (r, t) which is the discrete distribution function and has the dimension of mass density. The evolution equation for f q (r, c q ) is f q (r + c q ∆t, t + ∆t) = f q (r, t) + ∆t(Ω q,c + Ω q,f ), where Ω q,c is the collision operator and Ω q,f is an additional operator if a volumetric force densityf acts on the system. We use the Bhatnagar-Gross-Krook approximation for the collision operator [S7], which relaxes the system to the equilibrium state f eq q as where τ r is the relaxation time relating to the kinematic viscosity ν as τ r = ∆t(1/2 + ν/c 2 s ), with c s is the lattice sound speed, c s = ∆c 2 /3. The equilibrium distribution function f eq q = f eq q (ρ f , u f ) is calculated as where w q is the lattice weight that depends on the configurations.
For Ω q,f , we employ [S8] Ω q,f ∆c As a result, the macroscopic velocity is changed so the second term in Eq. (S1) becomes Handling the free surface of the fluid To simulate the free surface, we need to implement the mass tracking algorithm [S9-S11]. First, we assign a type of nodes such as the fluid, interface, or gas node for each node, where the interface node exists between the fluid and gas nodes as in Fig. S1. Note that Eqs. (S1) and (S2) are only used in the fluid and interface nodes.
A gas node represents the cell which is not occupied by the fluid, hence f q = 0. An interface node expresses the interface between fluid and gas, where the streaming and collision of f q exist as in fluid nodes. Here, we introduce a variable m f , which represents the density of the fluid in a single cell, to track the evolution of the surface. The interface node turns into a fluid node if m f ≥ ρ * f or into a gas node if m f ≤ 0, where ρ * f is the unit density of the fluid. Therefore, the state of each node is characterized by the liquid fraction λ: if the node is gas, where m f = λρ f . The evolution of the m f is determined by the balance between the populations streaming into the node f q ′ (r + c q ′ ∆t, t) (q ′ = −q) and out of the node f q (r, t) where α q is a function of λ of the neighboring node (located at r + c q ′ ∆t).
if f q ′ (r + c q ′ ∆t, t) streams from an interface node, 1 if f q ′ (r + c q ′ ∆t, t) streams from a fluid node, 0 if f q ′ (r + c q ′ ∆t, t) streams from a gas node.

(S9)
When an interface node turns into a fluid node, the neighboring gas nodes turn into interface nodes. When an interface node turns into a gas node, the neighboring fluid nodes turn into interface nodes. Although the density in a continuum model must be conserved, the discrete model can contain small loss or gain of m f . The surplus (or shortfall including the possibility of negative density) of m f is then computed at every time step and is corrected to satisfy the conservation among all interface nodes.
Fixed-pressure boundary condition. As stated before, LBM equations are solved only in the liquid and interface nodes. This creates a problem in the implementation since the population streaming to the interface nodes from gas nodes which is necessary in Eq. (S2) is not welldefined. Assuming that the gas node is always in equlibrium and has the same velocity of the interface node u in f and a constant atmospheric density ρ a < ρ f , the incoming distribution function (first term on the right hand side of Eq. (S2)) is replaced by the equilibrium distribution function with u in f and ρ a [S5], This is analogous to applying a fixed-pressure boundary condition at the interface, and local symmetry conditions for the velocity.
FIG. S1. Illustration of the the division of the lattice nodes into fluid, interface and gas nodes.
C. Solid boundaries and the fluid-particle coupling We implement two coupling schemes to handle solid boundaries within our simulations. We use the bounceback rules for no-slip boundary condition on walls and the surface of the impactor, while we use the direct forcing scheme for suspended particles.
The bounce-back rule simply states that whenever a population is streaming towards a wall, this population is reflected and bounced back in the opposite direction. This rule can be expressed as in LBM notation. If the wall is moving, the reflection has to take into account the momentum transfer by an addititonal term [S3, S4] {f q (r, t)−f q ′ (r, t+∆t) where u w is the wall velocity. Here, u w is calculated as where u I and ω I are the translational velocity and the angular velocity on the surface of the impactor, respectively, and R I denotes the center of mass of the impactor.  The momentum exchange described in Eq. (S12) results in a force on each node on the impactor surfaceF (r) as The hydrodynamic force on the impactor F h I is the sum of the forces for all nodes in the surface as F h I = r∈surfaceF (r), while T h I = r∈surface (r − R I ) ×F (r) is the hydrodynamic torque.
Direct forcing. By using the immersed boundary method, we calculate the hydrodynamic force through an additional discretization of particles into a set of segments r cell . These particle segments are related to the fluid simulation by an interpolating function [S12]. We implement the simplified version [S5, S11], where the segments correspond to the lattice nodes of the LBM r cell = r. Since the volume of a cubic cell is unity, the hydrodynamic force on each cellF cell (r) can be computed directly from the velocity differences where u cell is the velocity of the particle cell u cell (r) = u + (r − R) × ω, where u, R, and ω are the translational velocity, center of mass, and angular velocity of the suspended particles, respectively. The resultant hydrodynamic force on each suspended particle F h is the sum of all forces on the cells inside the particle l as F h = r∈lF cell (r). Similarly, the torque is given by T h = r∈l (r−R)×F cell (r). Note that this method requires a contribution to the body force density of the fluidf . Therefore, we calculatef in Eqs. (S6) and (S5) asf (r) = −ρ f gẑ −F cell (r) ∆x 3 . (S17) Note that the first term comes from the gravity.

II. ADDITIONAL FIGURES
In this section, we present several figures that are not placed in the main text. First, we plot the force exerted on the impactor for φ = 0.54 and u I 0,z /u * = 1.61 in Fig.  S2(a), when the impactor does not bounce. In contrast to Fig. 1(c) in the main text, no distinct peak exists in F I,z in this case. Moreover, there is a visible time delay between the maximum hydrodynamic contribution and the maximum contact contribution. It is easy to imagine that this response is viscous not to make the impactor rebounds. The plots for the height (z/a min ) of the impactor for different values of volume fraction φ can be seen in Fig. S2(b).
Normal displacement. In Fig. S3(a), we visualize the particle displacement in normal (z-) direction ∆z, sliced in the middle of the simulation box. Here, one can observe the existence of localized region of high normal displacements, which corresponds to the dynamically jammed region. Note that our observation is similar to the displacement field observed in an experiment of rod impactor [S13].
Shear stress. In Fig. S3(b), we visualize the dimensionless shear stress on each suspended particle σ xz /σ 0 , with σ 0 = 4 3 πρ f a min g, sliced in the middle of the simulation box. Here, one can observe that the magnitude of the shear stress is almost uniform. This indicates that the shear stress plays a minor role in the impact problem because it is unrelated to the force chains ( Fig. 3(b) in the main text). We also visualize the absolute value of the ratio of the shear to normal stresses |σ xz /σ zz | in Fig.  S3(c), which shows that the shear stress is much smaller than the normal stress. This clarifies another difference between the impact problem and the DST in which the shear stress is as large as the normal stress.