# Effective viscosity and elasticity in dense suspensions under impact

Pradipto<EMAIL_ADDRESS>Department of Mechanical Systems Engineering,
Tokyo University of Agriculture and Technology, Koganei, Tokyo 184-8588, Japan
111Present address Yukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawaoiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan  Hisao Hayakawa
Yukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawaoiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan

###### Abstract

The elastic response of dense suspensions under impact is investigated by
coarse-graining method of the coupled lattice Boltzmann method and discrete
element method (LBM-DEM) simulation. We are able to distinguish the viscous
force from the elastic force exerted on the impactor. We then propose an
empirical expression for the elastic force that can reproduce the behavior of
the LBM-DEM simulation. For the analytical treatment, we use the perturbation
method, which recovers the behavior of the impactor in a relatively short time
for the elastic force. We also discuss the impact of a foot-spring-body system
to mimic walking on the top of dense suspensions. The foot undergoes multiple
bounces and also multiple hops due to the spring force and the stiffness of
the suspensions. We also check the parameter dependencies of the hopping
motion and found that hopping and multiple rebounds are suppressed as the
spring stiffness increases.

††preprint: APS/123-QED

##  I Introduction

The phenomenon of being able to run on suspensions has attracted the interest
of both scientists and the general public [1, 2]. Such impact-induced
hardening of dense suspensions is often chosen as an example of discontinuous
shear thickening (DST) [1], but it has already been shown that the underlying
mechanism of impact-induced hardening is different from that of DST [3]. In
fact, impact-induced hardening is a transient process in which only normal
stress becomes large and the system is heterogeneous, whereas DST is a steady
process in which both shear and normal stresses become large and the system is
homogeneous.

Most physical studies of impact-induced hardening use a free-falling impactor
or a constant velocity penetrating intruder. Using a free-falling impactor,
Ref. [4] reported the existence of a localized rigid region under the
impactor, called the dynamically jammed region (DJR). As such a DJR grows in
size, Ref. [4] proposed the added-mass model, which treats the impact as an
inelastic collision between the impactor and the DJR that evolves with time.
Then, Ref. [5] visualized the flow field inside the dense suspension around
the penetrating intruder, and found that the strain rate peaked at the
boundary of the DJR. Inspired by this observation, Reference [6] proposed a
model based on the viscous force acting at the boundary of the DJR. However,
none of the above models can explain the existence of elastic response of
dense suspensions under impacts such as fracture [7], high stress near the
boundary [8] and rebound of the impactor [9]. In Ref. [8], a constitutive
model was proposed and the modulus of elasticity was measured when the DJR
spans from the impactor to the boundary. Then, the viscoelastic response of
dense suspensions under impact is captured using the floating + force chains
model [10], where the elastic force is only finite when the force chains of
contacting suspended particles percolate from the impactor to the bottom
boundary. However, such an analysis requires data on the position of the
suspended particles at a given time interval to resolve the force chains and
calculate the number of percolated force chains. In this paper, we propose a
model that does not depend on such simulation data to recover the motion of
the impactor. This can be achieved by studying the behavior of the viscous and
elastic forces acting on the impactor during the impact and constructing
models for them.

The motion of a running or walking person is more complicated than that of a
free-falling impactor or a penetrating intruder. An approach to directly
investigate the walking motion on the suspensions was described in Ref. [11].
They discussed the maximum penetration depth of a foot for different impact
velocities corresponding to walking, jogging, and running [8]. They also
showed that the added mass model is not sufficient to recover the response of
the suspensions under running motion. Nevertheless, little is known about the
dynamics of multiple impacts on dense suspensions, which are important for
running and walking motions. Some studies have attempted to reproduce
mechanical models for the locomotion of legged animals. One of the simplest
and most famous models is the spring-mass model inspired by biomechanical
observations [12]. Such a model has been realized as a one-legged hopping
robot [13]. In the spring-mass model, the human leg is represented by a spring
and the human body is simply represented by a point mass. Inspired by such a
model, this paper studies the motion of a foot-spring-body system on dense
suspensions.

((a))

((b))

Figure 1: Local quantities inside suspensions for $\phi=0.53$,
$u_{0}=4.5u^{*}$, $t=0.06t_{g}$: (a) Strain rate and (b) strain.

The structure of this paper is as follows. In Sec. II, we evaluate the
viscosity and elasticity around the impactor, starting from the setup for the
free-falling impactor simulation, the equation of motion of the impactor, the
coarse-graining technique to obtain relevant fields and evaluate the behavior
of the viscous and elastic force exerted on the impactor surface during the
impact process. Here, we also treat the elastic force as a perturbation to the
purely viscous floating model [10]. In Sec. III, we describe the simulation
setup for the foot-spring-body model and discuss the hopping motion of such a
system to clarify the criterion for the hopping motion. In Sec. IV, we
conclude our findings and discuss the future prospects of this study. In
Appendix A, we describe the details of the coupled lattice Boltzmann method
and discrete element method (LBM-DEM) used in our simulation. Finally, in
Appendix B, we compare our simulation results with relevant experiments.

##  II Evaluation of viscosity and elastic force around the impactor

###  II.1 Setup for a free-falling impactor simulation

We use the coupled LBM-DEM simulation as in Refs. [3, 10]. Details of the
simulation method are given in Appendix A. The simulation setup is as follows
[10]. A suspension with $N$ particles with a volume fraction
$\phi_{0}:=(4\pi/6)N(a_{\rm min}^{3}+a_{\rm max}^{3})$ is contained in a box
of size $W\times H\times D$, where $a_{\rm min}$ and $a_{\rm max}$ are the
radii of the smaller and larger suspended spheres, respectively. We only
analyze the case that the number of smaller spheres is equal to that of larger
spheres. Throughout this paper, we have assumed perfect density matching
between the solvent and the suspended particles, where the densities of the
particles and the solvent satisfy the relation $\rho_{p}=\rho_{f}$, where
$\rho_{p}$ and $\rho_{f}$ are the densities of the suspended particles and the
solvent fluid, respectively. In this section, a spherical impactor with
diameter $D_{I}$ (radius $a_{I}$) and density $\rho_{I}$ is released from
height $H_{0}$, which corresponds to the impact speed $u_{0}=\sqrt{2gH_{0}}$
with gravitational acceleration $g$. In our simulation $\rho_{I}$ and $D_{I}$
satisfy $\rho_{I}=4\rho_{f}$ and $D_{I}=6a_{\rm min}$ respectively, where
$a_{\rm min}$ is the radius of the smallest suspended particle. We also
introduce the time scale $t_{g}=\sqrt{a_{I}/2g}$, the velocity scale
$u^{*}=\sqrt{2ga_{I}}$ and the force scale
$F_{g}=\frac{4}{3}\pi\rho_{f}a_{I}^{3}g$. For most of the cases considered in
this paper, we use $\phi_{0}=0.53$, $H=3D_{I}$, $W=D=6D_{I}$ and $N=2200$.

###  II.2 Equation of motion of impactor

The equation of motion for a free-falling impactor along $z-$direction of the
center of mass $z_{I}$, and the deepest point of the impactor $z:=z_{I}-a_{I}$
onto dense suspensions can be written as

$m_{I}\ddot{z}=-m_{I}\tilde{g}-3\pi\eta_{\rm eff}\dot{z}z-F_{\rm el}(z),$ (1)

where $\tilde{g}$ is the effective gravity acceleration defined as
$\tilde{g}=g(\rho_{I}-\rho_{f})/\rho_{f}$, where $\rho_{I}$ and $\rho_{f}$ are
the density of the impactor and the suspensions respectively. Note that
gravity acts on the negative $z$ direction as shown in the first term on the
right-hand side of Eq. (1).

In order to propose a model that does not depend on the simulation data to
recover the motion of the impactor, we need to determine two unknown
variables: (i) the effective viscosity $\eta_{\rm eff}$ and (ii) the elastic
force $F_{\rm el}$.

###  II.3 Local quantities inside suspensions

((a))

((b))

((c))

Figure 2: (a) Visualization of the surface of the submerged impactor $S$, with
arrows represent the surface normals $\bm{n}$. (b) Plots of the time evolution
of effective volume fraction $\phi_{\rm eff}$ on $S$ for $\phi=0.53$ and
$u_{0}=4.5u^{*}$. (c) Plots of the time evolution of effective viscosity
$\eta_{\rm eff}$ on $S$ for $\phi=0.53$ and $u_{0}=4.5u^{*}$.

In this subsection, we will describe the method to obtain local quantities
inside suspensions that are relevant to elucidate the behavior of viscosity
and elastic force on the impactor during the impact process. The quantities we
are interested in are: (i) local volume fraction, (ii) local strain rate,
(iii) local strain and (iv) local stress. The approximate description of such
continuum fields from discrete particle data can be carried out using the
coarse-grain method that has been used for granular materials [14, 15]. The
local volume fraction $\phi(\bm{r})$ can be expressed with a coarse-grained
function $\Phi(\bm{r})$ as

$\phi(\bm{r},t):=\sum_{i}\Phi(\bm{r}-\bm{r}_{i}(t)),$ (2)

where $\bm{r}_{i}$ is the position of particle $i$. Here, we adopt

$\Phi(\bm{r}-\bm{r}_{i}):=\frac{1}{(w\sqrt{2\pi})^{3}}\exp\left[-\frac{(\bm{r}-\bm{r}_{i})^{2}}{2w^{2}}\right],$
(3)

where $\bm{r}=(x,y,z)$ and $\bm{r}_{i}=(x_{i},y_{i},z_{i})$ are the field and
particle positions respectively. For all the results presented here, we adopt
width $w=6a_{\rm min}$ and calculated all fields on a square grid of lattice
constant $0.5a_{\rm min}$. To satisfy the boundary conditions on the wall,
mirrored copies of the particle configurations are required on each side of
the wall before applying the coarse-graining methods [16].

The coarse-grained momentum density $\bm{p}(\bm{r},t)$ can be written as

$\bm{p}(\bm{r},t)=\sum_{i}m_{i}\bm{u}_{i}(t)\Phi(\bm{r}-\bm{r}_{i}(t)),$ (4)

where $\bm{u}_{i}$ and $m_{i}$ is the velocity and mass of particle $i$,
respectively. The velocity field $\bm{u}(\bm{r},t)$ is defined by
$\bm{u}(\bm{r},t)=\bm{p}(\bm{r},t)/\rho(\bm{r},t)$, with
$\rho(\bm{r},t):=\sum_{i}m_{i}\Phi(\bm{r}-\bm{r}_{i}(t))$. The stress tensor
$\overleftrightarrow{\sigma}$ contains the contributions from the contact
$\overleftrightarrow{\sigma}^{c}$ and the hydrodynamics
$\overleftrightarrow{\sigma}^{h}$

$\overleftrightarrow{\sigma}:=\overleftrightarrow{\sigma}^{h}+\overleftrightarrow{\sigma}^{c}.$
(5)

The contact contribution is calculated from the pairwise contact force

$\overleftrightarrow{\sigma}^{c}(\bm{r})=-\frac{1}{2}\sum_{i,j}\bm{F}^{c}_{ij}\otimes\bm{r}_{ij}\Phi(\bm{r}-\bm{r}_{i})$
(6)

where $\bm{F}^{c}_{ij}$ and $\bm{r}_{ij}$ are the pairwise contact force and
the interparticle distance between particles $i$ and $j$, respectively. Here
$\otimes$ denotes the tensor product. Meanwhile, the hydrodynamic contribution
is given by

$\overleftrightarrow{\sigma}^{h}(\bm{r})=\sum_{i}\overleftrightarrow{\sigma}_{i}^{h}\Phi(\bm{r}-\bm{r}_{i}),$
(7)

where $\overleftrightarrow{\sigma}_{i}^{h}$ is the hydrodynamic stress tensor
on each particle, obtained from the LBM and the lubrication stresslet [3].

The displacement field $\bm{U}(\bm{r})$ is defined by the particle
displacement from the equilibrium $\bm{U}_{i}(\bm{r}_{i})$ as

$\bm{U}(\bm{r})=\sum_{i}\bm{U}_{i}\Phi(\bm{r}-\bm{r}_{i}),$ (8)

where $\bm{U}_{i}$ is calculated as follows: For each time $t$, an additional
equilibration step is introduced where we freeze the motion of the impactor
and allow the suspended particles to equilibrate. Thus, the particle
configuration from the LBM simulation at each time $t$, $\bm{r}_{i}(t)$, is
considered as the initial condition in the equilibration process, i.e.
$\bm{r}_{i}(t,s=0)$. At each increment of the equilibration process $s+\Delta
s$, $\bm{r}_{i}(t,s)$ is updated considering only the hydrodynamic
lubrication, normal and tangential contacts, until the condition is reached
where the average overlap between particles $\langle\delta_{n}^{i}\rangle$ is
less than a threshold $\delta^{th}$, at $s=s_{th}$. Then the displacement from
the equilibrium of particle $i$ at time $t$ is defined as

$\bm{U}_{i}(t):=\bm{r}_{i}(t;s=0)-\bm{r}_{i}(t;s=s_{th}),$ (9)

where the second term on the right hand side of Eq. (9) is the equilibrated
position.

Once the flow field is obtained, the symmetric part of the local strain rate
tensor can be obtained $\overleftrightarrow{D}(\bm{r})$.

$\overleftrightarrow{D}(\bm{r}):=\frac{1}{2}(\nabla\bm{u}+\nabla\bm{u}^{T}).$
(10)

Meanwhile, the local strain tensor $\overleftrightarrow{L}(\bm{r})$ is defined
as

$\overleftrightarrow{L}(\bm{r}):=\frac{1}{2}(\nabla\bm{U}+\nabla\bm{U}^{T}).$
(11)

Let us introduce the scalar local viscosity $\eta(\bm{r})$ defined as [17,
18],

$\eta(\bm{r}):=\frac{1}{2}\frac{\overleftrightarrow{\sigma}(\bm{r}):\overleftrightarrow{D}(\bm{r})}{\overleftrightarrow{D}(\bm{r}):\overleftrightarrow{D}(\bm{r})},$
(12)

where $:$ is the scalar or double inner product. The local strain rate
$\dot{\varepsilon}$ is defined as

$\dot{\varepsilon}(\bm{r}):=\sqrt{2\overleftrightarrow{D}(\bm{r}):\overleftrightarrow{D}(\bm{r})}.$
(13)

A snapshot of the local strain rate field immediately after impact is shown in
Fig. 1(a). It can be seen that the position of the high rate region in our
simulation is reminiscent of that observed experimentally in Ref. [5],
although they used a constant penetrating intruder. Then the local viscous
stress $\sigma^{\rm(vis)}$ is simply given by [17, 18]

$\sigma^{\rm(vis)}(\bm{r})=\eta(\bm{r})\dot{\varepsilon}(\bm{r}).$ (14)

Finally, similar to the strain rate, the local scalar strain fields
$\varepsilon(\bm{r})$ are defined as

$\varepsilon(\bm{r}):=\sqrt{2\overleftrightarrow{L}(\bm{r}):\overleftrightarrow{L}(\bm{r})}.$
(15)

A snapshot of the local strain field immediately after the impact is shown in
Fig.1(b). Similar to eq. (12) the local rigidity can be introduced as

$G(\bm{r}):=\frac{1}{2}\frac{\overleftrightarrow{\sigma}(\bm{r}):\overleftrightarrow{L}(\bm{r})}{\overleftrightarrow{L}(\bm{r}):\overleftrightarrow{L}(\bm{r})}.$
(16)

The local elastic stress is then given by

$\sigma^{(\rm el)}(\bm{r})=G(\bm{r})\varepsilon(\bm{r}).$ (17)

###  II.4 Density around impactor: phenomenology for the viscosity

Once the local fields are obtained, we define a surface $S$ with normals
$\bm{n}$ corresponding to the surface of the submerged impactor (Fig. 2(a)).
First, the effective volume fraction around the impactor $\phi_{\rm eff}$ is
defined as

$\phi_{\rm eff}:=\frac{\int_{S}\phi(\bm{r})dA}{\int_{S}dA}.$ (18)

Similarly, the effective viscosity around the impactor $\eta_{\rm eff}$ is
defined as

$\eta_{\rm eff}:=\frac{\int_{S}\eta(\bm{r})dA}{\int_{S}dA}.$ (19)

The time evolution of $\phi_{\rm eff}$ can be seen in Fig. 2(b), where it can
be seen that $\phi_{\rm eff}$ increases immediately after the impact. We have
also plotted the time evolution of $\eta_{\rm eff}$, which can be seen in Fig.
2(c), where it also increases immediately after the impact. This suggests that
the effective viscosity $\eta_{\rm eff}$ can be estimated by the $\phi_{\rm
eff}-$dependent constitutive law for viscosity [19]

$\frac{\eta_{\rm eff}}{\eta_{0}}=\bigg{(}\frac{\phi_{\rm
eff}(\phi_{J}-\phi_{0})}{\phi_{0}(\phi_{J}-\phi_{\rm eff})}\bigg{)}^{2},$ (20)

where $\phi_{J}$ is the jamming volume fraction, $\phi_{\rm eff}$ is the
effective volume fraction around the impactor, and $\phi_{0}$ is the initial
volume fraction. In Fig. 2(c), we also plot Eq. (20) as dashed lines, where
one can see that the measurement (Eq. (19)) agrees with the model (Eq. (20))

###  II.5 Empirical model for the elastic force

((a))

((b))

Figure 3: (a) Plots of the force on the impactor from Eqs. (21) and (22),
black dashed lines represent the force measured directly on the impactor, blue
solid, blue dashed and blue dotted lines are obtained by Eqs. (22) and (23),
respectively. (b) Plots of elastic force against depth for various impact
velocities. Dashed lines represent the solution of Eq. (23) by numerically
obtaining $z_{\rm me}$ and $z_{\rm cut}$, and set $z_{\rm on}=0.18D_{I}$.

The force around the impactor can be obtained by integrating the stress in the
surface $S$. This way, we can also obtain the contribution from viscous and
elastic parts,

$\displaystyle F_{\rm vis}$ $\displaystyle=\int_{S}\sigma^{(\rm
vis)}\bm{n}dA,$ (21) $\displaystyle F_{\rm el}$
$\displaystyle=\int_{S}\sigma^{(\rm el)}\bm{n}dA,$ (22)

where $dA$ is the surface integral on the surface $S$. In Fig. 3(a), we plot
the total force obtained through this analysis (blue solid line) and compare
it to the one directly measured on the impactor (black dashed line). Although
the coarse-grained method based on Eqs. (22) and (23) cannot follow the sharp
pulses in short time as well as a little deviation between two measurements
can be detected after a peak, the agreement between two methods is reasonable.
Then, as one can see that viscous part dominates, yet the elastic part cannot
be ignored nonetheless. Also, the elastic force appears later than its viscous
counterpart.

To get a better understanding of the elastic force, we plot $F_{\rm el}$
against normalized depth $z/D_{I}$ for various $u_{0}$ with $\phi_{0}=0.53$ in
Fig. 3(b). Here, one can see that for various $u_{0}$, the onset of elastic
forces takes place at similar $z/D_{I}$, then grow linearly until reaching the
maximum value at certain $z_{\rm me}$ that depends on impact velocity $u_{0}$.

Based on these observations, we propose the following empirical expression for
the elastic force $F_{\rm el}$

$F_{\rm el}(z)=\begin{cases}0,\qquad z<z_{\rm on},\\\ k(z-z_{\rm on}),\qquad
z_{\rm on}\leq z\leq z_{\rm me},\\\ -k^{\prime}(z-z_{\rm me})+k(z_{\rm
me}-z_{\rm on}),\qquad z_{\rm me}<z\leq z_{\rm cut},\\\ 0,\qquad z>z_{\rm
cut},\end{cases}$ (23)

where $z_{\rm on}$, $z_{\rm me}$, $z_{\rm cut}$ are the position of the onset
of elastic force, the position of maximum elastic force, and the cut-off
position of the elastic force, respectively. Here, $k$ and $k^{\prime}$ are
fitting parameters that have the dimension of spring constant and are related
to the stiffness of the suspended particles.

The next task is to determine $z_{\rm on}$, $z_{\rm me}$, $z_{\rm cut}$. Using
the second expression for $F_{\mathrm{el}}(z)$, Eq. (1) can be rewritten as

$m_{I}\ddot{\tilde{z}}=-m_{I}\tilde{g}-3\pi\eta_{\rm
eff}\dot{\tilde{z}}\tilde{z}-k\tilde{z},\qquad t_{\rm on}\leq t<t_{\rm me}$
(24)

with $\tilde{z}:=z-z_{\rm on}$. It is natural that the maximum elastic force
takes place when the sign of impactor velocity switches ($\dot{\tilde{z}}=0$).
However, the full solution of Eq. (24) is complicated thus it is impossible to
determine $z_{\rm me}$ ($z$-position that corresponds to $\dot{\tilde{z}}=0$)
analytically. Thus, we solve Eq. (24) numerically and denote the rebound time
or time when $\dot{\tilde{z}}=0$ as $t_{\rm me}$, thus $z_{\rm me}=z(t_{\rm
me})$.

For $z_{\rm me}<z<z_{\rm cut}$, the equation of motion can be written as

$m_{I}\ddot{z}=-m_{I}\tilde{g}-3\pi\eta_{\rm eff}\dot{z}z-k(z_{\rm me}-z_{\rm
on})+k^{\prime}(z-z_{\rm me}).$ (25)

The mechanical energy of the system $E$ that consists of kinetic and elastic
energy is given as

$E=\frac{1}{2}m_{I}\dot{z}^{2}+\frac{1}{2}k^{\prime}(z-z_{\rm me})^{2}.$ (26)

Multiplying Eq. (25) with $\dot{z}$, one can get

$\displaystyle\frac{dE}{dt}$ $\displaystyle=-\dot{z}(m_{I}g+k(z_{\rm
me}-z_{\rm on}))-3\pi\eta\dot{z}^{2}z$
$\displaystyle+2k^{\prime}\dot{z}(z-z_{\rm me})$ (27)

Note that the elastic force is continuous at $z_{\rm me}$ and $E_{\rm
me}=E(t=t_{\rm me})=\frac{1}{2}kz_{\rm me}^{2}$. Then, $z_{\rm cut}$ takes
place when all of the remaining energy is completely dissipated at $t_{\rm
cut}$ i.e

$\displaystyle 0$ $\displaystyle=\frac{1}{2}kz_{\rm me}^{2}-\int_{t=t_{\rm
me}}^{t=t_{\rm cut}}\dot{z}(m_{I}g+k(z_{\rm me}-z_{\rm on}))$
$\displaystyle+3\pi\eta\dot{z}^{2}z-2k^{\prime}\dot{z}(z-z_{\rm me})dt.$ (28)

From Eq. (28), one can obtain $t_{\rm cut}$ numerically, then determine
$z_{\rm cut}=z(t_{\rm cut})$. The dashed lines in Fig. 3(b) are the solution
of Eq. (23) with $z_{\rm me}$ and $z_{\rm cut}$, leaving $z_{\rm on}$, $k$,
and $k^{\prime}$ as the remaining fitting parameters. One can see that this
approach yields reasonable agreement with the simulation.

###  II.6 Perturbation approach

((a))

((b))

Figure 4: Plots of the solution of the perturbation approach, compared with
the solution of the empirical model (Eqs. (1) and (23) with $z_{\rm
on}=0.18D_{I}$ and numerically obtained $z_{\rm me}$ and $z_{\rm cut}$) and
the simulation results for: (a) Depth of the impactor and (b) velocity of the
impactor.

In this section, let us attempt to solve Eq. (24) by treating the elastic
force as a perturbation to the dominant viscous force. Let us introduce the
dimensionless depth $\zeta=\tilde{z}/a_{I}$ and the dimensionless time
$\tau:=3\pi\eta_{\rm eff}a_{I}t/m_{I}$. Then, Eq. (24) can be written as

$\frac{d^{2}\zeta}{d\tau^{2}}=-G-\frac{d\zeta}{d\tau}\zeta-\epsilon\zeta,$
(29)

where $G=m_{I}^{2}\tilde{g}/9\pi\eta_{\rm eff}^{2}a_{I}^{3}$ and
$\epsilon=m_{I}k/(3\pi\eta_{\rm eff}a_{I})^{2}$. The perturbation solution for
Eq. (29)

$\zeta=\zeta_{0}+\epsilon\zeta_{1}+O(\epsilon^{2}).$ (30)

Plugging Eq. (30) into Eq. (29), ignoring higher order terms, and rearranging,
one can get up to the first order in $\epsilon$

$\displaystyle\frac{d^{2}\zeta_{0}}{d\tau^{2}}+G+\frac{d\zeta_{0}}{d\tau}\zeta_{0}$
$\displaystyle=0,$ (31)
$\displaystyle\frac{d^{2}\zeta_{1}}{d\tau^{2}}+\zeta_{0}\frac{d\zeta_{1}}{d\tau}+\frac{d\zeta_{0}}{d\tau}\zeta_{1}+\zeta_{0}$
$\displaystyle=0$ (32)

Equation (31) is equivalent to the floating model in Ref. [10] so the solution
with initial conditions $\zeta_{0}(0)=0$ and
$\zeta_{0}^{\prime}(0)=\tilde{u}_{0}$, with dimensionless impact velocity
$\tilde{u}_{0}:=u_{0}m_{I}/3\pi\eta a_{I}^{2}$, can be written in terms of
Airy functions

$\zeta_{0}(\tau)=\frac{\kappa[-\text{Ai}^{\prime}(\Phi)\text{Bi}^{\prime}(\Theta)+\text{Ai}^{\prime}(\Theta)\text{Bi}^{\prime}(\Phi)]}{\gamma[\text{Bi}(\Phi)\text{Ai}^{\prime}(\Theta)-\text{Ai}(\Phi)\text{Bi}^{\prime}(\Theta)]},$
(33)

where $\gamma=-G^{2/3}$, $\kappa=2^{2/3}G$,
$\Phi=(-G\tau+\tilde{u}_{0})/(2^{1/3}(-G)^{2/3})$, and
$\Theta=(\tilde{u}_{0})/(2^{1/3}(-G)^{2/3})$. Here, $\text{Ai}(x)$ is the Airy
function of the first kind, which is defined as
$\text{Ai}(x)=\int_{0}^{\infty}\cos(t^{3}/3+xt)dt/\pi$, and
$\text{Ai}^{\prime}(x)$ is its derivative. $\text{Bi}(x)$ is the Airy function
of the second kind, which is defined as
$\text{Bi}(x)=\int_{0}^{\infty}[\exp(-t^{3}/3+xt)+\sin(-t^{3}/3+xt)]dt/\pi$,
and $\text{Bi}^{\prime}(x)$ is its derivative.

The solution for Eq. (32) can be written in terms of the complementary and
particular solutions

$\zeta_{1}(\tau)=\zeta_{1,c}(\tau)+\zeta_{1,p}(\tau).$ (34)

The complementary solution $\zeta_{1,c}(\tau)$ can be solved by first solving
the homogeneous equation

$\frac{d^{2}\zeta_{1}}{d\tau^{2}}+\zeta_{0}\frac{d\zeta_{1}}{d\tau}+\frac{d\zeta_{0}}{d\tau}\zeta_{1}=0.$
(35)

Then, the complementary solution can be written as

$\displaystyle\zeta_{1,c}(\tau)$
$\displaystyle=c_{a}\zeta_{a}(\tau)+c_{b}\zeta_{b}(\tau),$ (36)
$\displaystyle\zeta_{a}(\tau)$
$\displaystyle=\frac{2^{\frac{1}{3}}\gamma\Lambda_{1}(\Theta,\Phi)+(Gt-
u_{0})\Lambda_{2}(\Theta,\Phi)-2^{\frac{4}{3}}\gamma\Lambda_{3}(\Theta,\Phi)+2^{\frac{1}{3}}\gamma\Lambda_{4}(\Theta,\Phi)}{G\left[\text{Ai}^{\prime}(\Theta)\text{Bi}(\Phi)-\text{Ai}(\Phi)\text{Bi}^{\prime}(\Theta)\right]^{2}},$
(37) $\displaystyle\Lambda_{1}(\Theta,\Phi)$
$\displaystyle=\text{Ai}^{\prime}(\Phi)^{2}\text{Bi}^{\prime}(\Theta)^{2},$
$\displaystyle\Lambda_{2}(\Theta,\Phi)$
$\displaystyle=[\text{Ai}^{\prime}(\Theta)\text{Bi}^{\prime}(\Phi)-\text{Ai}(\Phi)\text{Bi}^{\prime}(\Theta)]^{2},$
$\displaystyle\Lambda_{3}(\Theta,\Phi)$
$\displaystyle=\text{Ai}^{\prime}(\Theta)\text{Ai}^{\prime}(\Phi)\text{Bi}^{\prime}(\Theta)\text{Bi}^{\prime}(\Phi),$
$\displaystyle\Lambda_{4}(\Theta,\Phi)$
$\displaystyle=\text{Ai}^{\prime}(\Theta)^{2}\text{Bi}^{\prime}(\Phi)^{2},$
$\displaystyle\zeta_{b}(\tau)$
$\displaystyle=\frac{1}{\left[\text{Ai}^{\prime}(\Theta)\text{Bi}(\Phi)-\text{Ai}(\Phi)\text{Bi}^{\prime}(\Theta)\right]^{2}},$
(38)

where $c_{a}$ and $c_{b}$ are coefficients that will be determined later from
the initial conditions. Then, the particular solution $\zeta_{1,p}$ can be
obtained when we have finite Wronskian,
$W(\zeta_{a}(\tau),\zeta_{b}(\tau))\neq 0$ defined as

$W(\zeta_{a}(\tau),\zeta_{b}(\tau)):=\zeta_{a}(\tau)\zeta_{b}^{\prime}(\tau)-\zeta_{b}(\tau)\zeta_{a}^{\prime}(\tau).$
(39)

Plugging Eqs. (37) and (38), one can obtain

$W(\zeta_{a}(\tau),\zeta_{b}(\tau))=-\frac{1}{\left[\text{Ai}^{\prime}(\Theta)\text{Bi}(\Phi)-\text{Ai}(\Phi)\text{Bi}^{\prime}(\Theta)\right]^{2}}.$
(40)

The particular solution can be written as

$\displaystyle\zeta_{1,p}(\tau)=\zeta_{a}(\tau)$
$\displaystyle\int\frac{\zeta_{b}(\tau)\zeta_{0}(\tau^{\prime})}{W(\zeta_{a}(\tau^{\prime}),\zeta_{b}(\tau^{\prime}))}d\tau$
$\displaystyle-\zeta_{b}(\tau)\int\frac{\zeta_{a}(\tau)\zeta_{0}(\tau^{\prime})}{W(\zeta_{a}(\tau^{\prime}),\zeta_{b}(\tau^{\prime}))}d\tau^{\prime},$
(41) $\displaystyle=\zeta_{a}(\tau)$
$\displaystyle\left[\log\left(\zeta_{b}(\tau)\right)+C\right]$
$\displaystyle-\zeta_{b}(\tau)\int\frac{\zeta_{a}(\tau)\zeta_{0}(\tau^{\prime})}{W(\zeta_{a}(\tau^{\prime}),\zeta_{b}(\tau^{\prime}))}d\tau^{\prime}.$
(42)

Note that the integral in the second term of the r.h.s in Eq. (42) cannot be
calculated analytically. Plugging Eqs. (42) and (36) to Eq. (34), one can get

$\displaystyle\zeta_{1}(\tau)=$
$\displaystyle\zeta_{a}(\tau)\left[\log(\zeta_{b}(\tau))+C_{1}\right]$
$\displaystyle+\zeta_{b}(\tau)\left[C_{2}-\int\frac{\zeta_{a}(\tau)\zeta_{0}(\tau^{\prime})}{W(\zeta_{a}(\tau^{\prime}),\zeta_{b}(\tau^{\prime}))}d\tau^{\prime}\right],$
(43)

where $C_{1}$ and $C_{2}$ are constants that will be determined from the
initial conditions. Since the integral in the second term in the r.h.s of Eq.
(43) cannot be calculated analytically, a numerical evaluation for this
equation is necessary.

Let us discuss the appropriate initial conditions for this perturbation
problem. Note that the perturbative solution only exists ($\epsilon\neq 0$)
when elastic force exists after $\tau_{\rm on}$

$\zeta(\tau)=\begin{cases}\zeta_{0}(\tau),\qquad 0\leq\tau<\tau_{\rm on},\\\
\zeta_{0}(\tau)+\epsilon\zeta_{1}(\tau),\qquad\tau\geq\tau_{\rm
on}.\end{cases}$ (44)

Thus, the initial conditions for Eq. (29) are $\zeta(\tau_{\rm on})=\zeta_{\rm
on}$ and $\frac{d\zeta}{d\tau}(\tau_{\rm on})=\tilde{u}_{\rm on}$. Currently,
$\tau_{\rm on}=3\pi\eta_{\rm eff}a_{I}t_{\rm on}/m_{I}$ is another fitting
parameter and the value that correspond to $\zeta_{\rm on}=z_{\rm on}/a_{I}$
used in Fig. 3(b) is chosen. Then, $\tilde{u}_{\rm on}=u_{\rm
on}m_{I}/3\pi\eta a_{I}^{2}$ can be obtained from the solution of the floating
model after specifying $\tau_{\rm on}$. Then, after ignoring the higher-order
terms

$\displaystyle\zeta_{0}(\tau_{\rm on})+\epsilon\zeta_{1}(\tau_{\rm on})$
$\displaystyle=\zeta_{\rm on}$ (45)
$\displaystyle\frac{d\zeta_{0}}{d\tau}(\tau_{\rm
on})+\epsilon\frac{d\zeta_{1}}{d\tau}(\tau_{\rm on})$
$\displaystyle=\tilde{u}_{\rm on}.$ (46)

Eq. (46) must be valid for all $\epsilon$ close to zero. Nevertheless,
assuming some fluctuations in $F_{\rm el}$ one can allow
$d\zeta_{1}/d\tau(\tau_{\rm on})$ to have some finite but small value
$u_{1,\rm on}$. Therefore, the initial conditions can be written as

$\displaystyle\zeta_{0}(\tau_{\rm on})$ $\displaystyle=\zeta_{\rm
on},\quad\zeta_{1}(\tau_{\rm on})=0,$ (47)
$\displaystyle\frac{d\zeta_{0}}{d\tau}(\tau_{\rm on})$
$\displaystyle=\tilde{u}_{\rm on},\quad\frac{d\zeta_{1}}{d\tau}(\tau_{\rm
on})=u_{1,\rm on}.$ (48)

With this initial condition, $C_{1}$ and $C_{2}$ in Eq. (43) can be expressed
as

$\displaystyle C_{1}=$ $\displaystyle\frac{1}{\zeta_{a}(\tau_{\rm
on})}\Bigg{[}-\zeta_{b}(\tau_{\rm on})\int_{\tau_{\rm
on}}^{\tau}\frac{\zeta_{a}(\tau)\zeta_{0}(\tau^{\prime})}{W(\zeta_{a}(\tau^{\prime}),\zeta_{b}(\tau^{\prime}))}d\tau^{\prime}-C_{2}\zeta_{b}(\tau_{\rm
on})-\zeta_{a}(\tau_{\rm on})\log\zeta_{b}(\tau_{\rm on})\Bigg{]},$ (49)
$\displaystyle C_{2}=$ $\displaystyle\frac{\zeta_{b}(\tau_{\rm
on})\left[\frac{\zeta_{a}(\tau_{\rm on})\zeta_{0}(\tau_{\rm
on})}{W(\zeta_{a}(\tau_{\rm on}),\zeta_{b}(\tau_{\rm on}))}+\int_{\tau_{\rm
on}}^{\tau}\frac{\zeta_{a}(\tau)\zeta_{0}(\tau^{\prime})}{W(\zeta_{a}(\tau^{\prime}),\zeta_{b}(\tau^{\prime}))}d\tau^{\prime}\right]+\log\zeta_{b}(\tau_{\rm
on})\left[\zeta_{a}^{\prime}(\tau_{\rm on})-\zeta_{a}(\tau_{\rm
on})\right]}{\zeta_{b}(\tau_{\rm on})-\zeta_{b}^{\prime}(\tau_{\rm on})}$
$\displaystyle+$ $\displaystyle\frac{\zeta_{b}^{\prime}(\tau_{\rm
on})\left[\frac{\zeta_{a}(\tau_{\rm on})}{\zeta_{b}(\tau_{\rm
on})}-\int_{\tau_{\rm
on}}^{\tau}\frac{\zeta_{a}(\tau)\zeta_{0}(\tau^{\prime})}{W(\zeta_{a}(\tau^{\prime}),\zeta_{b}(\tau^{\prime}))}d\tau^{\prime}\right]-u_{1,\rm
on}}{\zeta_{b}(\tau_{\rm on})-\zeta_{b}^{\prime}(\tau_{\rm on})},$ (50)

where $\zeta_{a}^{\prime}=d\zeta_{a}/d\tau$ and
$\zeta_{b}^{\prime}=d\zeta_{b}/d\tau$. Then, we can solve Eq. (31) with
$\tilde{u}_{0}$ as its initial condition ($\tau=0$). Then at $\tau_{\rm on}$
we solve Eq. (30) with Eqs. (47) and (48) as initial conditions and stitch
this with the solution of Eq. (31). As one can see in Fig. (4), we get a good
agreement between the simulation results, Eq. (1) with Eq. (23), and the
perturbative approach until the impactor reach minimum velocity. After the
impactor reaches the minimum velocity, it starts sinking due to the relaxation
of suspensions. Indeed, the sinking/relaxation process is currently ignored in
Eq. (29).

##  III Foot-spring-body dynamics in dense suspensions

((a))

((b))

((c))

((d))

((e))

Figure 5: a) Schematic of the simulation setup of the foot-spring-body system.
The body (black circle) is connected with the foot (black rectangle) by a
massless spring (yellow tube). (b) Time evolutions of the velocities of the
foot (blues) and body (purples) in $z-$coordinate, where the black dashed line
expresses $u_{z}=0$, the solid lines express the solutions of Eq. (51), and
the symbols express the results of the simulations. (c) Time evolutions of the
positions of the foot (blues) and body (purples) in $z-$coordinate, where the
black dashed line expresses the suspension surface, the solid lines express
the solutions of Eq. (51), and the symbols express the results of the
simulations. All results here are obtained with $k_{s}=100m_{0}/(a_{\rm
min}t_{g}^{2})$ and $u_{0}=4u^{*}$. (d) Time evolutions of the foot position
in $z-$direction $z_{p}$ with various spring stiffness $k_{s}$ at
$u_{0}=4u^{*}$. (e) Time evolutions of the foot position in $z-$direction
$z_{p}$ with various initial velocity $u_{0}$ at $k_{s}=100m_{0}/(a_{\rm
min}t_{g}^{2})$.

Let us now introduce the foot-spring-body model as a model for expressing the
bouncing motion on suspended liquids. The foot in the foot-spring-body model
is represented by a rectangular plate impactor with volume $V_{p}=W_{p}\times
H_{p}\times D_{p}$ and mass $m_{p}=V_{p}\rho_{p}$, where $\rho_{p}$ is the
density of the foot plate and we assume $\rho_{p}=1.2\rho_{f}$. The body is
represented by a sphere with diameter $D_{b}$ and mass $m_{b}$. We take the
density of the body $\rho_{b}$ as $\rho_{b}=2\rho_{p}=2.4\rho_{f}$. The body
and foot are then connected by a massless spring with stiffness $k_{s}$ and
natural length $L_{0}$. The schematic of this setup can be seen in figure
5(a). Here we use $\phi=0.51$, $H=2D_{b}$, $W=D=4D_{b}$ and $N=618$. Note that
we are only interested in the vertical ($z$ direction) motion of the system.
We then modify Eq. (1) to the set of equations

$\displaystyle m_{b}\ddot{z}_{b}=$ $\displaystyle-
m_{b}g-k_{s}(z_{b}-z_{p}-L_{0})-\zeta_{s}\dot{z}_{b}$ $\displaystyle
m_{p}\ddot{z}_{p}=$ $\displaystyle-m_{p}g-3\pi\eta_{\rm
eff}z_{p}\dot{z}_{p}+F_{\rm el}(z_{p})$
$\displaystyle+k_{s}(z_{b}-z_{p}-L_{0})-\zeta_{s}\dot{z}_{p},$ (51)

where $m_{b}$ and $z_{b}$ are the mass and the vertical position of the body,
respectively, $z_{p}$ is the vertical position of the base of the plate
impactor, $\zeta_{s}$ is the damping constant,
$\zeta_{s}=\sqrt{k_{s}(m_{p}+m_{b})/2}$. We also use Eq. (23) to recover the
elastic force $F_{\rm el}(z_{p})$ in this case. A typical motion of the foot-
spring-body system can be seen in Figs. 5(b) and 5(c). It can be seen that the
solution of eq. (51) agrees well with the results of the simulation. Thus, our
empirical model Eq. (51) is the reasonable model to analyze the motion of the
foot-spring-body model. Initially, the foot experiences a strong deceleration
as in the free-falling impactor due to the interaction between the foot and
the suspensions. Meanwhile, the body continues to accelerate due to gravity.
The system then exhibits a damped oscillation. Due to the spring force and the
stiffness of the suspensions, the foot undergoes multiple bounces ($u_{p}<0$)
and also multiple jumps ($z_{p}>0$). This result suggests that composite
materials with elastic springs inside the body can maintain their position for
a while.

We then evaluate the multiple rebound motion of the foot in detail. First, we
check how the motion of the foot depends on the stiffness of the spring
$k_{s}$. The simulation results for different $k_{s}$ are shown in Fig. 5(d).
Here one can see a lower tendency to multiple bounces or hops for higher
$k_{s}$. Furthermore, for a rigid system ($k_{s}\rightarrow\infty$) the foot
only bounces once and then sinks. This is similar to the prediction of the
added mass model in Ref. [11], where they suggested that running on
suspensions is impossible for a perfectly stiff leg. We then examine the
dependence of the initial velocity ($u_{0}$) in Fig. 5(e). As expected, the
foot sinks and does not hop at low $u_{0}$, since the impact-induced hardening
is stronger at high $u_{0}$ [4, 9, 6, 3].

##  IV Conclusions and discussions

By using the coarse-grained function, defining the deformation as the
displacement of the suspended particles from equilibrium, and evaluating the
local quantities on the surface of the impactor, we are able to separate the
viscous and elastic forces acting on the impactor and evaluate their
contributions. We observe increases in viscosity and density around the
impactor immediately after impact. Then, by looking at how $F_{\rm el}$
depends on depth $z$, we propose an empirical expression with five fitting
parameters (onset of elastic force $z_{\rm on}$, position of maximum elastic
force $z_{\rm me}$, position $z_{\rm cut}$ where the initial mechanical energy
$E$ becomes completely dissipated, spring constants $k$ and $k^{\prime}$), two
of which ($z_{\rm me}$ and $z_{\rm cut}$) can be obtained from the equations
of motion and energy of the impactor. Alternatively, we also solve the
equation of motion of the impactor when the elastic force is linear with
depth, treating the elastic force as a perturbation to the dominant viscous
force, which agrees well with the simulation results before the suspensions
start to relax. To mimic the hopping motion, we discussed the impact of the
foot-spring-body system on the top of dense suspensions. Our model agrees well
with the simulation results, and we found a lower tendency for hopping and
multiple rebounds at higher spring stiffness $k_{s}$ and low initial velocity
$u_{0}$.

Our simulation method can be extended to capture the time evolution of the
sinking intruder in dense suspensions, where oscillations and slip-stick
motions have been observed [20]. However, such a sinking process is beyond the
scope of this paper. Relatedly, little is known about the relaxation process
of hardening suspensions after impact, where our current perturbation approach
fails and our toy model can only give a rough approximation ($F_{\rm el}=0$
when $z>z_{\rm cut}$) and ignores the actual relaxation process of the
suspensions [21, 22, 23]. Future studies should focus on this relaxation
phenomenon of dense suspensions under impact.

## Acknowledgements

We thank Ryohei Seto and Satoshi Takada for their useful comments. This work
is partially supported by the Grant-in-Aid of MEXT for Scientific Research
KAKENHI (Grant No. JP21H01006). All numerical calculations were carried out at
the Yukawa Institute for Theoretical Physics (YITP) Computer Facilities, Kyoto
University, Japan.

##  Appendix A  LBM + DEM with free surface

We employ the LBM involving suspensions and the free surface of the fluid. The
details of the LBM are explained in Ref. [3]. The suspended particles in LBM
are represented as a group of solid nodes, while the surrounding fluids are
represented by fluid nodes. The hydrodynamic field is calculated from the time
evolution of the discrete distribution function at each fluid node. We select
the lattice unit $\Delta x=0.2a_{\text{min}}$, where it gives sufficient
accuracy but is still not computationally expensive as shown in the previous
LBM for suspensions literature [24, 25, 26]. In addition, to simulate the free
surface of the fluid, it is necessary to introduce interface nodes between the
fluid and gas nodes [27, 28, 29, 3].

Equations of motion and the torque balance of particle $i$ are, respectively,
given by

$m_{i}\frac{d\bm{u}_{i}}{dt}=\bm{F}_{i}^{c}+\bm{F}_{i}^{h}+\bm{F}_{i}^{\rm
lub}+\bm{F}_{i}^{r}+\bm{F}_{i}^{g},$ (52)
$I_{i}\frac{d\bm{\omega}_{i}}{dt}=\bm{T}_{i}^{c}+\bm{T}_{i}^{\rm
lub}+\bm{T}_{i}^{h}.$ (53)

Here, $\bm{u}_{i}$, $\bm{\omega}_{i}$, $m_{i}$, and
$I_{i}=(2/5)m_{i}a_{i}^{2}$ (with $a_{i}$ the radius of particle $i$), are the
translational velocity, angular velocity, mass, and the moment of inertia of
particle $i$, respectively. $\bm{F}_{i}^{g}=-m_{i}g\hat{\bm{z}}$ is the
gravitational force acting on the suspended particle $i$, where $\hat{\bm{z}}$
is the unit vector in the vertical direction.

Note that our LBM accounts for both the short-range lubrication force
$\bm{F}_{i}^{\rm lub}$ and torque $\bm{T}_{i}^{\rm lub}$, as well as the long-
range hydrodynamic force $\bm{F}_{i}^{h}$ and torque $\bm{T}_{i}^{h}$ as in
Ref. [26, 30]. The long-range parts ($\bm{F}_{i}^{h}$ and $\bm{T}_{i}^{h}$)
are calculated using the direct forcing method [29, 3], while the lubrication
force $\bm{F}_{i}^{\rm lub}$ and torque $\bm{T}_{i}^{\rm lub}$ are expressed
by pairwise interactions as $\bm{F}_{i}^{\rm lub}=\sum_{j\neq
i}\bm{F}_{ij}^{\rm lub}$ and $\bm{T}_{i}^{\rm c}=\sum_{j\neq
i}\bm{T}_{ij}^{\rm lub}$, respectively [31, 32, 26, 30]. The explicit
expressions of $\bm{F}_{ij}^{\rm lub}$ and $\bm{T}_{ij}^{\rm lub}$ can be
found in Ref. [30].

We adopt the linear spring-dashpot version of the DEM [33] for the contact
interaction between particles, which involves both the normal and the
tangential contact forces. Note that we omit the dissipative part for the
tangential contact force. For the particle $i$, the contact force
$\bm{F}^{c}_{i}$ and torque $\bm{T}^{c}_{i}$ are, respectively, written as
$\bm{F}^{c}_{i}=\sum_{i\neq
j}(\bm{F}^{\text{nor}}_{ij}+\bm{F}^{\text{tan}}_{ij})$ and
$\bm{T}^{c}_{i}=\sum_{i\neq j}a_{i}\bm{n}_{ij}\times\bm{F}^{\text{tan}}_{ij}$,
where $a_{i}$ is the radius of particle $i$. The normal force is explicitly
expressed as

$\bm{F}^{\text{nor}}_{ij}=(k_{n}\delta_{ij}^{n}-\zeta^{(n)}u_{ij}^{(n)})\bm{n}_{ij},$
(54)

where $k_{n}$ is the spring constant, $\delta_{ij}^{n}$ is the normal overlap,
$\bm{n}_{ij}$ is the normal unit vector between particles, $u_{ij}^{(n)}$ is
the normal velocity difference of the contact point
$u_{ij}^{(n)}=u_{i}^{(n)}-u_{j}^{(n)}$, and $\zeta^{(n)}=\sqrt{m_{0}k_{n}}$ is
the damping constant, where $m_{0}$ is the average mass of the suspended
particles. If the tangential contact force is smaller than a slip criterion,
the tangential contact force is represented as

$\tilde{\bm{F}}^{\text{tan}}_{ij}=k_{t}\delta_{ij}^{t}\bm{T}_{ij},$ (55)

where $k_{t}$, assumed to be $0.2k_{n}$, is the tangential spring constant,
$\delta_{ij}^{t}$ is the tangential compression and $\bm{T}_{ij}$ is the
tangential unit vector at the contact point between particles $i$ and $j$. We
adopt the Coulomb friction rules as

$\displaystyle|\bm{F}_{ij}^{\text{tan}}|$
$\displaystyle=\mu|\bm{F}_{ij}^{\text{nor}}|\quad\text{if
}|\tilde{\bm{F}}^{\text{tan}}_{ij}|\geq\mu|\bm{F}^{\text{nor}}_{ij}|\quad\text{(slip)},$
(56) $\displaystyle|\bm{F}_{ij}^{\text{tan}}|$
$\displaystyle=|\tilde{\bm{F}}^{\text{tan}}_{ij}|\quad\text{if
}|\tilde{\bm{F}}^{\text{tan}}_{ij}|\leq\mu|\bm{F}^{\text{nor}}_{ij}|\quad\text{(stick)},$
(57)

whereas $\delta_{ij}^{t}$ is updated each time with relative tangential
velocity [33].

Finally, $\bm{F}_{i}^{r}$ is the electrostatic repulsive force, also expressed
by pairwise interactions as $\bm{F}_{i}^{\rm r}=\sum_{j\neq
i}\bm{F}_{ij}^{r}$. The explicit expression of $\bm{F}^{r}_{ij}$ is expressed
by the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory [34, 35, 36] for the
double layer electrostatic force as

$\bm{F}^{r}_{ij}=F_{0}\exp(-h/\lambda)\bm{n}_{ij},$ (58)

where
$F_{0}=k_{B}T\lambda_{B}\hat{Z}^{2}(e^{a_{\text{min}}/\lambda}/(1+a_{\text{min}}/\lambda))^{2}/h^{2}$
with the charge number $\hat{Z}$, the Bjerrum length $\lambda_{B}$ and the
Debye-Hückel length $\lambda$. Note that $\lambda_{B}$ can be expressed as
$\lambda_{B}=e^{2}/(4\pi\epsilon_{0}\epsilon_{r}k_{B}T)$ where $e$,
$\epsilon_{0}$, $\epsilon_{r}$, and $k_{B}$ are the elementary charge, the
vacuum permittivity, the dielectric constant, and the Boltzmann constant,
respectively [36]. Here, we adopt the Debye length
$\lambda=0.02a_{\text{min}}$. Our simulation ignores the Brownian force. Thus,
the electrostatic repulsion force is important to prevent the suspended
particles from clustering [30, 32].

For the impactor, the contact force $\bm{F}^{I,c}$ and torque $\bm{T}^{I,c}$,
which arise from the interactions with the suspended particles, are also
calculated by the DEM. The lubrication force $\bm{F}^{I,\rm lub}$ and torque
$\bm{T}^{I,\rm lub}$ are also calculated in a similar manner as used in
suspended particles. The long-range hydrodynamic force $\bm{F}^{I,h}$ and
torque $\bm{T}^{I,h}$ are calculated using the bounce-back rule which
satisfies the no-slip boundary condition between the fluid and the surface of
the impactor [24, 25]. In the bounce-back rule, the LBM discrete distribution
function that streams from fluid nodes to the boundary nodes is reflected.
Then, the hydrodynamic force on each node is calculated from the momentum
transferred in this reflection process. In our implementation, the bounce-back
rule is implemented by treating the surface of the impactor as boundary nodes.

##  Appendix B  Comparison with experiments

Figure 6: Plots of the velocity of the impactor for simulations (reds on the
top $x-$axis) and experiments of Ref. [9] (blacks on the bottom $x-$axis).

In this section, we compare the results of LBM-DEM simulations for the
velocity of a free-falling impactor into dense suspensions with a
corresponding experiment in the same setup [9]. In Fig. 6, we plot the time
evolution of velocities for the same impact velocity. We found a good
agreement between the simulation and experiment once we introduce the scaled
time and change the volume fraction. Simulations require a higher volume
fraction than the experiments since the interaction between particles in LBM-
DEM does not consider the rolling friction, which must exist in the actual
consider cornstarch particles. It has been known that the introduction of
rolling friction lowered the critical volume fraction of discontinuous shear
thickening [37]. Another important difference is on the timescale, whereas
simulation is twice as fast compared to the experiments. This is indeed caused
by the finite size effect in our simulation, as it was observed earlier in
Ref. [3], the time scale of impact can be scaled by the number of particles.

## References

  * Brown and Jaeger [2014] E. Brown and H. M. Jaeger, Shear thickening in concentrated suspensions: phenomenology, mechanisms and relations to jamming, Rep. Prog. Phys. 77, 046602 (2014).
  * Ness  _et al._ [2022] C. Ness, R. Seto, and R. Mari, The physics of dense suspensions, Annu. Rev. Condens. Matter Phys. 13, 97 (2022).
  * Pradipto and Hayakawa [2021a] Pradipto and H. Hayakawa, Impact-induced hardening in dense frictional suspensions, Phys. Rev. Fluids 6, 033301 (2021a).
  * Waitukaitis and Jaeger [2012] S. R. Waitukaitis and H. M. Jaeger, Impact-activated solidification of dense suspensions via dynamic jamming fronts, Nature 487, 205 (2012).
  * Han  _et al._ [2016] E. Han, I. R. Peters, and H. M. Jaeger, High-speed ultrasound imaging in dense suspensions reveals impact-activated solidification due to dynamic shear jamming, Nat. Commun. 7, 12243 (2016).
  * Brassard  _et al._ [2021] M. Brassard, N. Causley, N. Krizou, J. A. Dijksman, and A. H. Clark, Viscous-like forces control the impact response of shear-thickening dense suspensions, J. Fluid. Mech. 923, A38 (2021).
  * Roche  _et al._ [2013] M. Roche, E. Myftiu, M. C. Johnston, P. Kim, and H. A. Stone, Dynamic fracture of nonglassy suspensions, Phys. Rev. Lett. 110, 148304 (2013).
  * Maharjan  _et al._ [2018] R. Maharjan, S. Mukhopadhyay, B. Allen, T. Storz, and E. Brown, Constitutive relation for the system-spanning dynamically jammed region in response to impact of cornstarch and water suspensions, Phys. Rev. E 97, 052602 (2018).
  * Egawa and Katsuragi [2019] K. Egawa and H. Katsuragi, Bouncing of a projectile impacting a dense potato-starch suspension layer, Phys. Fluids 31, 053304 (2019).
  * Pradipto and Hayakawa [2021b] Pradipto and H. Hayakawa, Viscoelastic response of impact process on dense suspensions, Phys. Fluids 33, 093110 (2021b).
  * Mukhopadhyay  _et al._ [2018] S. Mukhopadhyay, B. Allen, and E. Brown, Testing constitutive relations by running and walking on cornstarch and water suspensions, Phys. Rev. E 97, 052604 (2018).
  * Blickhan [1989] R. Blickhan, The spring-mass model for running and hopping, J. Biomech. 22, 1217 (1989).
  * Raibert and Tello [1986] M. H. Raibert and E. R. Tello, Legged robots that balance, IEEE Expert 1, 89 (1986).
  * Zhang  _et al._ [2010] J. Zhang, R. P. Behringer, and I. Goldhirsch, Coarse-graining of a physical granular system, Prog. Theor. Phys. 184, 16 (2010).
  * Saitoh and Hayakawa [2013] K. Saitoh and H. Hayakawa, Quantitative test of the time dependent gintzburg-landau equation for sheared granular flow in two dimensions, Phys. Fluids 25, 070606 (2013).
  * Ries  _et al._ [2014] A. Ries, L. Brendel, and D. E. Wolf, Coarse graining strategies at walls, Comp. Part. Mech. 1, 177 (2014).
  * Giusteri and Seto [2018] G. Giusteri and R. Seto, A theoretical framework for steady-state rheometry in generic flow conditions, J. Rheol. 623, 713 (2018).
  * Giusteri and Seto [2021] G. G. Giusteri and R. Seto, Shear jamming and fragility of suspensions in a continuum model with elastic constraints, Phys. Rev. Lett. 127, 138001 (2021).
  * Boyer  _et al._ [2011] F. Boyer, E. Guazzelli, and O. Pouliquen, Unifying suspension and granular rheology, Phys. Rev. Lett. 107, 188301 (2011).
  * von Kann  _et al._ [2011] S. von Kann, J. H. Snoeijer, D. Lohse, and D. van der Meer, Nonmonotonic settling of a sphere in a cornstarch suspension, Phys. Rev. E 84, 060401(R) (2011).
  * Maharjan and Brown [2017] R. Maharjan and E. Brown, Giant deviation of a relaxation time from generalized newtonian theory in discontinuous shear thickening suspensions, Phys. Rev. Fluids 2, 123301 (2017).
  * Cho  _et al._ [2022] J. H. Cho, A. H. Griese, I. R. Peters, and I. Bischofberger, Lasting effects of discontinuous shear thickening in cornstarch suspensions upon flow cessation, Phys. Rev. Fluids 7, 063302 (2022).
  * Barik and Majumdar [2022] S. Barik and S. Majumdar, Origin of two distinct stress relaxation regimes in shear jammed dense suspensions, Phys. Rev. Lett. 128, 258002 (2022).
  * Ladd [1994a] A. J. C. Ladd, Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 1. theoretical foundation, J. Fluid. Mech. 271, 285 (1994a).
  * Ladd [1994b] A. J. C. Ladd, Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 2. numerical simulations, J. Fluid. Mech. 271, 311 (1994b).
  * Nguyen and Ladd [2002] N. Q. Nguyen and A. J. C. Ladd, Lubrication corrections for lattice-boltzmann simulations of particle suspensions, Phys. Rev. E 66, 046708 (2002).
  * Švec  _et al._ [2012] O. Švec, J. Skocek, H. Stang, M. R. Geiker, and N. Roussel, Free surface flow of a suspension of rigid particles in a non-newtonian fluid: A lattice boltzmann approach, J. Non-Newton. Fluid. 179-180, 32 (2012).
  * Leonardi  _et al._ [2014] A. Leonardi, F. K. Wittel, M. Mendoza, and H. J. Hermann, Coupled dem-lbm method for the free-surface simulation of heterogeneous suspensions, Comp. Part. Mech. 1, 3 (2014).
  * Leonardi  _et al._ [2015] A. Leonardi, M. Cabrera, F. K. Wittel, R. Kaitna, M. Mendoza, W. Wu, and H. J. Herrmann, Granular-front formation in free-surface flow of concentrated suspensions, Phys. Rev. E 92, 052204 (2015).
  * Pradipto and Hayakawa [2020] Pradipto and H. Hayakawa, Simulation of dense non-brownian suspensions with the lattice boltzmann method: shear jammed and fragile states, Soft Matter 16, 945 (2020).
  * Seto  _et al._ [2013] R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Discontinuous shear thickening of frictional hard-sphere suspensions, Phys. Rev. Lett. 111, 218301 (2013).
  * Mari and Seto [2014] R. Mari and R. Seto, Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions, J. Rheol. 58, 1693 (2014).
  * Luding [2008] S. Luding, Cohesive, frictional powders: contact models for tension, Granul. Matter 10, 235 (2008).
  * Derjaguin and Landau [1941] B. Derjaguin and L. D. Landau, Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solutions of electrolytes, Acta Physicochim. U.R.S.S 14, 633 (1941).
  * Verwey and Overbeek [1948] E. J. W. Verwey and J. T. G. Overbeek,  _Theory of the Stability of Lyophobic Colloids: The Interaction of Sol Particles Having an Electric Double Layer_ (Elsevier, 1948).
  * Israelachvili [2011] J. Israelachvili,  _Intermolecular and Surface Forces_, 3rd ed. (Academic Press, 2011).
  * Singh  _et al._ [2020] A. Singh, C. Ness, R. Seto, J. J. de Pablo, and H. M. Jaeger, Shear thickening and jamming of dense suspensions: The “roll” of friction, Phys. Rev. Lett. 124, 248005 (2020).