\section{Introduction}
\input{./collision_modeling/introduction}
\section{Governing equations}
\input{./collision_modeling/governing_eq}
\section{Contact detection}
\input{./collision_modeling/contact_detection}
\section{{\label{sec:collision_model}}Collision model}
\input{./collision_modeling/Collision_model}
\section{{\label{sec:lubrication_model}}Lubrication model}
\input{./collision_modeling/Lubrication_model}
\section{{\label{sec:sensitivity_analysis}}Sensitivity analysis}
\input{./collision_modeling/Sensitivity_analysis}
\section{Validation}
\input{./collision_modeling/Validation}
\section{Conclusions}
\input{./collision_modeling/conclusion}
\newcommand{\clearpage}{\clearpage}


\section*{Acknowledgements}
	All the simulations presented here were carried out at the Center for Information Services and High Performance Computing (ZIH), TU Dresden, Germany. The authors thank ZIH for providing the computational resources and support. RJ gratefully acknowledges the scholarship provided by the Land Sachsen, No. L-201535. TS was employed by the Institute of Fluid Mechanics, TU Dresden at the time this study was carried out.

\bibliographystyle{abbrvnat}


\subsection{Continuous phase}
The fluid field is described by the unsteady three dimensional Navier-Stokes equations for a Newtonian fluid of constant density, 
\begin{equation} 
 \frac{\partial{\vec{u}}}{\partial{t}} = - \nabla\cdot(\vec{u}\otimes\vec{u}) + \frac{1}{\rho_\mathrm{f}} \: \nabla\cdot\boldsymbol{\tau}+\vec{f}_\mathrm{v}+\vec{f}_\mathrm{IBM} \hspace{0.1cm},
\end{equation} 
\begin{equation}
\nabla\cdot\vec{u}=0     \hspace{0.1cm} .
\end{equation}
which are solved on an equidistant Eulerian grid in the domain $\Omega = \Omega^{\mathrm{f}} \cup \Omega^{\mathrm{p}}$, i.e. the union of the fluid domain and the particle interior (Fig. \ref{fig:sketch_domain}). This nomenclature differs from \cite{tschisgale2018}. Here, $\textbf{u}=(u,v,w)^{T}$, $p$, $\rho_\mathrm{f}$, and $t$ are the velocity vector along the Cartesian coordinates $x,y,z$, pressure, fluid density, and time, respectively. On the right-hand side, $\vec{f}_\mathrm{IBM}$ is the coupling force imposed to take into account the particle-fluid interaction, and $\vec{f}_\mathrm{v}$ represents the volume force applied to drive the flow. The latter is constant in space and adjusted in time to obtain the desired flow rate. The hydrodynamic stress tensor, $\boldsymbol{\tau}$, is 
\begin{equation}
 \boldsymbol{\tau}= -p \:\vec{\mathbb{I}} + \mu_\mathrm{f} \:(\nabla \:\vec{u} + \left(\:\nabla \:\vec{u} \:\right)^\top ) \hspace{0.1cm},
\end{equation}
with $\mu_\mathrm{f}$ the dynamic viscosity of the fluid, and $\vec{\mathbb{I}}$ the identity matrix.\\

\begin{figure}
	\hspace*{0.2\linewidth}
	\def0.7\linewidth{0.6\linewidth}
	\input{./collision_modeling/Figures/Particle.pdf_tex}
	\caption{An instantanous snapshot of an ellipsoid as a moving particle. The nomenclature is introduced in the text.}
	\label{fig:sketch_domain}
\end{figure}

\subsection{Disperse phase}
The solid phase is represented by a non-deformable ellipsoid defined by its longest axis, $a$, intermediate axis, $b$, and smallest axis, $c$, as shown in Fig. \ref{fig:sketch_domain}. Here, $\vec{x}_\mathrm{p}$ is the position of particle center, $\vec{x}_\mathrm{\Gamma}$ an arbitrary point on the particle surface $\Gamma$, and $\vec{r} = \vec{x}_\mathrm{\Gamma}-\vec{x}_\mathrm{p}$ the vector connecting these points. The translational momentum balance of a particle is given by
\begin{equation}\label{eq:part_lin}
m_\mathrm{p} \:\frac{\text{d}\vec{u}_\mathrm{p}}{\text{d} t} =  \int_{\Gamma} \!\boldsymbol{\tau} \cdot \vec{n} \:\text{d} S   + V_\mathrm{p}\:( \rho_\mathrm{p}-\rho_\mathrm{f} ) \:(\vec{g}+\vec{f}_\mathrm{v}) + \vec{f}_\mathrm{c} + \vec{f}_\mathrm{lub}  \hspace{0.1cm},
\end{equation}
Here, $V_\mathrm{p}$, $\rho_\mathrm{p}$, and $m_\mathrm{p}$ are the volume, mass density, and mass of the ellipsoid, respectively, while $\textbf{u}_\mathrm{p}$ is the translational velocity of the particle, $\textbf{n}$ the outward-pointing surface normal vector, and $\textbf{g}$ the gravitational acceleration. The collision force $\vec{f}_\mathrm{c}$ and lubrication forces $\vec{f}_\mathrm{lub}$ are provided in the later sections below. The hydrodynamic loads are integrated over the particle surface $\Gamma$, which envelopes the particle domain $\Omega^{\mathrm{p}}$. The temporal evolution of the center of mass of the particle $\textbf{x}_\mathrm{p}$ is obtained by integrating
\begin{equation}\label{eq:part_lin_xp}
\frac{\text{d}\vec{x}_\mathrm{p}}{\text{d} t} = \vec{u}_\mathrm{p} \hspace{0.1cm}.
\end{equation}
The angular momentum balance is formulated in a body-fixed and co-rotated frame ${\mathcal{P}}$. This local coordinate system is defined to coincide with the principal axes of the rigid particle (Fig.~\ref{fig:sketch_domain}), such that the tensor of inertia $\textbf{I}_\mathrm{p}$ is diagonal, i.e.
\begin{equation}\label{eq:inertia_tensor}
\textbf{I}_\mathrm{p}=\text{diag}\left(\text{I}_1,\text{I}_2,\text{I}_3 \right)  \hspace{0.1cm}.
\end{equation}
The relation between the global and the local co-rotated frame is established by a rotation in terms of a unit quaternion $\textbf{q}_\mathrm{p} \in \mathbb{S}^3$, where $\mathbb{S}^3$ is the set of unit quaternions \citep{Kuipers1999}. It represents the spatial rotation of the particle around its center of mass $\textbf{x}_\mathrm{p}$. The temporal evolution of $\textbf{q}_\mathrm{p}$ is connected to the local angular velocity of the particle $\boldsymbol{\Omega}_\mathrm{p}$ via
\begin{equation}\label{eq:omega_quat}
\frac{\text{d} \textbf{q}_\mathrm{p}}{\text{d} t} = \frac{1}{2} \: \textbf{q}_\mathrm{p} \circ \boldsymbol{\Omega}_\mathrm{p}  \hspace{0.1cm},
\end{equation}
where $\dots$$ \circ \dots$ is the multiplication of two quaternions. Any vector $\vec{v}$ can be rotated from one coordinate system to the other using 
\begin{subequations}
	\begin{align}
	\label{eqn:forward_rot}
	\textbf{v} &= \textbf{q} \circ \textbf{V} \circ \bar{\textbf{q}} \hspace*{1cm} \text{(forward rotation)} \\
	\label{eqn:backward_rot}
	\textbf{V} &= \bar{\textbf{q}} \circ \textbf{v} \circ \textbf{q} \hspace*{1.09cm} \text{(backward rotation)}
	\end{align}
\end{subequations}
with $\bar{\textbf{q}}$ indicating the conjugate of a quaternion ${\textbf{q}}$.
As a convention, capital letters, e.g. $\textbf{V} \in \mathbb{R}^3$, are given in the local co-rotated frame and small letters, $\textbf{v} \in \mathbb{R}^3$, in the global Eulerian frame. Consequently, the angular velocity in the global Eulerian frame can be expressed in terms of $\boldsymbol{\Omega}_\mathrm{p}$, by
\begin{equation}\label{eq:transformation}
\boldsymbol{\omega}_\mathrm{p} = \textbf{q}_\mathrm{p} \circ \boldsymbol{\Omega}_\mathrm{p} \circ \bar{\textbf{q}}_\mathrm{p}  \hspace{0.1cm}.
\end{equation}
The angular momentum balance in the local coordinate system, then, reads
\begin{equation}\label{eq:part_ang}
\textbf{I}_\mathrm{p} \! \cdot \! \frac{\text{d} \boldsymbol{\Omega}_\mathrm{p}}{\text{d} t} + \boldsymbol{\Omega}_\mathrm{p} \times \textbf{I}_\mathrm{p} \!\! \cdot \! \boldsymbol{\Omega}_\mathrm{p} =  \int_{\Gamma} \textbf{R} \times \left( \textbf{T} \cdot \textbf{N} \right) \:\text{d} S + \vec{M}_\mathrm{c} + \vec{M}_\mathrm{lub} \hspace{0.1cm},
\end{equation}
where $\textbf{N}$ the outward pointing surface normal vector in the frame ${\mathcal{P}}$, $\textbf{R}$ the vector of a surface point with respect to the center of mass of the particle, $\textbf{T}$ the hydrodynamic stress tensor, $\vec{M}_\mathrm{c}$ the torque due to the collision and $\vec{M}_\mathrm{lub}$ the torque due to the lubrication force.
\subsection{Coupling between the two phases}
The coupling of fluid and particle equations is thoroughly discussed in \citet{tschisgale2018}. The final discretized model is recalled here and is the starting point for the derivation of the collision model. The discretized linear equation of motion of an ellipsoid coupled with a fluid using the IBM presented in that reference is
\begin{equation}\label{eq:lin_par}
\textbf{u}_\mathrm{p}^n - \textbf{u}_\mathrm{p}^{n-1} = \:\left( m_\mathrm{p}+m_L \right)^{-1} \: \Delta{t}\left\{ \vec{f}_\mathrm{f} + \vec{f}_\mathrm{g} + \vec{f}_\mathrm{c} + \vec{f}_\mathrm{lub} +\vec{f}_\mathrm{v} \right\} \hspace{0.1cm},
\end{equation}
where $m_{L}$ is the mass of the fluid layer ${L}$ of thickness $d_{L}$ constituting the central element in this method with all details and justifications provided in the cited reference. The superscript $n$ denotes the time step, while 
\begin{subequations}\label{eq:f_fandf_g}
\begin{align}
\vec{f}_\mathrm{f} &=- \: \rho_\mathrm{f} \! \int_L \:\frac{\textbf{u}_\Gamma^{n-1} - \tilde{\textbf{u}}^{} }{\Delta{t}}\: \text{d} V + \frac{1}{\Delta{t}} \left[\int_{\Omega^{\mathrm{p}}} \rho_\mathrm{f} \:\textbf{u} \:\text{d} V  \right]^{\:n}_{\:n-1}  \hspace{0.1cm},\\
\vec{f}_\mathrm{g} &= \: V_\mathrm{p}\:( \rho_\mathrm{p}-\rho_\mathrm{f} ) \:\textbf{g} \hspace{0.1cm},
\end{align}
\end{subequations}
are the time-discrete fluid and gravitational force, respectively. 
Here, $\textbf{u}_\Gamma$ is the velocity vector at a point on the surface of the moving particle and $\tilde{\textbf{u}}$ is the fluid velocity without taking into account the disperse phase.

Analogously, the angular equation of motion can be formulated as
\begin{equation}\label{eq:ang_par}
\boldsymbol{\Omega}_\mathrm{p}^n - \boldsymbol{\Omega}_\mathrm{p}^{n-1} = - \left( \textbf{I}_\mathrm{p}+\textbf{I}_L \right)^{-1} \left\{  \int^{n}_{n-1} \!\! \boldsymbol{\Omega}_\mathrm{p} \times \textbf{I}_\mathrm{p} \!\! \cdot \! \boldsymbol{\Omega}_\mathrm{p} \: \text{d} t \right\} + \Delta{t} \left( \textbf{I}_\mathrm{p}+\textbf{I}_L \right)^{-1} \! \cdot \! \left\{  \vec{M}_\mathrm{f}  + \vec{M}_\mathrm{c} + \vec{M}_\mathrm{lub} \right\} \hspace{0.1cm},
\end{equation}
with the torque due to the fluid forces, 
\begin{equation}\label{eq:T_f}
\begin{split}
\vec{M}_\mathrm{f} = & - \: \rho_\mathrm{f} \! \int_L \textbf{R}\times \left(\frac{\textbf{U}_\Gamma^{n-1} - \tilde{\textbf{U}}^{} }{\Delta{t}} \right) \text{d} V + \frac{1}{\Delta{t}} \left[\bar{\textbf{q}}_\mathrm{p} \! \circ \int_{\Omega^{\mathrm{p}}} \rho_\mathrm{f} \:\textbf{r}\times \textbf{u} \:\text{d} V \:\: \circ \textbf{q}_\mathrm{p}  \right]^{\:n}_{\:n-1} \hspace{0.1cm}.
\end{split}
\end{equation}
Here, $\textbf{I}_L$ is the tensor of inertia introduced by the layer $L$ shown in Fig. \ref{fig:sketch_domain}.
Considering Eq. \eqref{eqn:backward_rot}, $\textbf{U}_\Gamma$ and $ \tilde{\textbf{U}}^{}$ are the same as $\textbf{u}_\Gamma$ and $\tilde{\textbf{u}}$, respectively, but in the co-rotated frame $\mathcal{P}$ of the respective particle. The vector joining a point on the surface of an ellipsoid to the center of that ellipsoid in the local frame of reference is denoted as $\vec{R}$ with $\vec{R} = \bar{\textbf{q}} \circ \vec{r} \circ \textbf{q}$ and $\vec{r}$ the same vector in the global frame.

\subsection{Normal collision of a spherical particle with a wall}

The benchmark rebound trajectories provided by \citet{gondret2002} are used to calibrate the model. Therefore, they can be reproduced fairly well, as shown in Fig. \ref{fig:lub_calib}a. Many further cases were investigated in the cited experiment to determine the dependency of the restitution coefficient for a normal collision, 
\begin{equation}
e_\mathrm{n} =  -\frac{\lvert \vec{u}_\mathrm{p,n,out} \rvert}{\lvert \vec{u}_\mathrm{p,n,in} \rvert} \hspace{0.1cm},
\end{equation}
as a function of the Stokes number. Here, $\vec{u}_\mathrm{p,n,in}$ and $\vec{u}_\mathrm{p,n,out}$ are the input and output particle normal velocities measured just before and after the lubrication zone, i.e. when the surfaces have a distance of $d_\mathrm{lub}$ identified with start and end of the collision process. Both velocities are normal to the colliding wall. The influence of the viscous fluid surrounding the particle can be highlighted very well when using the restitution coefficient of dry collision $e_\mathrm{d,n}$ as a reference, which is done in Fig. \ref{fig:validation_sphere}a for the data of \citet{gondret2002} and \citet{joseph2001}. The experimental behavior of $e_\mathrm{n}/e_\mathrm{d,n}$ versus $St_\mathrm{r}$ is reproduced very well using this model. All previously developed collision models have been compared with this curve in the literature. \\

The same situation was simulated with the new collision model using the same fluid and particle properties as mentioned in Tab. \ref{tab:sensitivity_analysis1} and a computational domain of size $13D_\mathrm{p}\times26D_\mathrm{p}\times13D_\mathrm{p}$ discretized with $256\times512\times256$ cells. The particle was dropped from a distance $y_\mathrm{p,0} = \SI{0.06}{m}$ to the wall. The Stokes number was varied by successively changing the gravitational acceleration from $g=0.02\times\SI{9.81}{m/s^2}$ to  $g=20\times\SI{9.81}{m/s^2}$, similar to \citet{biegert2017}. The parameter values, $d_\mathrm{lub}=2\Delta_x$, $D_\mathrm{p}/\Delta_x \approx 20$, $\mathit{CFL} = 0.5$ were used in these tests. Fig. \ref{fig:validation_sphere}a shows that the dependency of the restitution coefficient on $St_\mathrm{r}$ obtained in the simulations with the proposed model is captured very well. Observing the considerable experimental scatter a better agreement can hardly by expected. 
\subsection{Oblique collision of a spherical particle with a wall}
\label{sec:sph_obl_col}
\begin{figure*}[tp]
	\centering
	\setlength{\unitlength}{1cm}
	\begin{picture}(16,8)
	\put(0,0){\includegraphics{./collision_modeling/Figures/validation_sphere.pdf}}
	\put(8,0){\includegraphics{./collision_modeling/Figures/validation_oblique.pdf}}
	\put( 0. , 7.5 )     {a) }
	\put( 8. , 7.5 )     {b) }  
	\end{picture}
	\caption{Results of the simulations conducted to validate the collision model. a) Normalized restitution coefficient plotted against the particle Stokes number. The present simulation results are compared with the experiments of \citet{gondret2002} and \citet{joseph2001}. b) Present simulations without and with fluid to validate the collision model for an oblique wall collision. The experimental data of \citet{joseph2004} are used for reference.}
	\label{fig:validation_sphere}
\end{figure*}
The model also covers tangential forces so that it is also validated for oblique particle-wall collisions in a dry system and in the presence of a viscous liquid. A computational domain of $26D_\mathrm{p}\times13D_\mathrm{p}\times13D_\mathrm{p}$ with resolution $D_\mathrm{p}/\Delta_x \approx 20$ was used in these simulations. The particle was initially placed at $(X,Y,Z)=(4D_\mathrm{p},7D_\mathrm{p},7D_\mathrm{p})$ with its motion is driven by the acceleration 
\begin{equation}
\vec{g} = (\,g \sin(\phi_\mathrm{g}),-g\cos(\phi_\mathrm{g}),0 \,)^T \hspace{0.1cm}.
\end{equation}

Here, $\phi_\mathrm{g}$ is the angle between the wall-normal axis of the computational grid and gravity. Since the distance between the starting point of the trajectory is relatively short and the trajectory prior to collision is a straight line, so that the incidence angle of the particle with the wall-normal $\phi_\mathrm{in}=\phi_\mathrm{g}$ in these simulations. The magnitude of the imposed acceleration equals $g=20\times\SI{9.81}{m/s^2}$ to achieve the high Stokes numbers of the experiments. \\

The simulations are compared with the experiments of \citet{joseph2004} on oblique particle-wall collisions in air and aqueous solutions of glycerol.
\begin{table*} [tp]
	\centering
	\begin{tabular}{ccclclcc}
		\hline \hline
		$D_\mathrm{eq}\,[\mathrm{mm}]$ & $\rho_\mathrm{p}\,[\mathrm{kg\,m^{-3}}]$	& $\rho_\mathrm{f}\,[\mathrm{kg\,m^{-3}}]$ & $\nu_\mathrm{f}\,[\mathrm{m^2\,s^{-1}}]$	&  $e_\mathrm{d,n}$&$\mu_\mathrm{s}$ & $\mu_\mathrm{k}$ & $e_\mathrm{d,t}$\\
		\hline
		$2.5$ & $2540$ 	& - & -  & $0.97$	& $0.15$    & $0.15$ & $0.0$   \\
		$2.5$ & $2540$ 	& - & -  & $0.97$	& $0.15$    & $0.15$  & $0.39$ \\
		$2.5$ & $7800$ 	& $998$ & $1.007\times10^{-6}$  & $0.97$	& $0.02$    & $0.02$  & $0.34$ \\
	
		\hline
		\hline
	\end{tabular}
	\caption{Particle and fluid properties used in the simulations to validate the new model for oblique collisions. The parameters are similar to the ones used in the experiments of \citet{joseph2004}. Cases without fluid properties provided refer to dry collisions.}	
	\label{tab:oblique_collision}
\end{table*}
\citet{joseph2004} compared the rebound angle $\Psi_\mathrm{out}$ to the impact angle $\Psi_\mathrm{in}$, where
\begin{equation}
\Psi_\mathrm{in} = \frac{\lvert \vec{u}_\mathrm{c,t,in} \rvert}{\lvert \vec{u}_\mathrm{c,n,in} \rvert} \hspace{0.1cm},  \, \, \,
\Psi_\mathrm{out} = \frac{\lvert \vec{u}_\mathrm{c,t,out} \rvert}{\lvert \vec{u}_\mathrm{c,n,in} \rvert} \hspace{0.1cm},
\end{equation}
so that, e.g., $\Psi_\mathrm{in} = \tan(\phi_\mathrm{in})$. Following \cite{joseph2004}, the velocities at the contact point are employed here, instead of the particle center, as previously used in the definition of $e_\mathrm{n}$ above. Since in these tests the spherical particle is not rotating significantly before the collision the velocity of the contact point, $\vec{u}_\mathrm{c,in}$ is the same as the velocities at the particle center,  $\vec{u}_\mathrm{p,in}$. However, the tangential velocity of the contact point after the collision differs from the center point velocity and is calculated as $\vec{u}_\mathrm{c,t,out} = (\vec{u}_\mathrm{c}\!\cdot\vec{t})\,\vec{t}$. \\

Fig. \ref{fig:validation_sphere}b displays the results of the simulations and provides a comparison with the experiments for the same configuration. The dry collision was first simulated under the assumption of full stick at the contact point i.e. $e_\mathrm{d,t} = 0$. The results are represented with the diamond symbols. They match very well when the particle slides. The negative values in the rolling mode are not captured properly, however, as expected. This is obtained by introducing $e_\mathrm{d,t} = 0.39$ in Eq. \eqref{eq:edt}, as proposed by \citet{joseph2004}. The match between the results with the new model and the reference data is extremely good so that the validation is successfully accomplished.

\subsection{Oblique collisions of ellipsoidal particles} \label{ssec:Coll_NS}
\begin{figure}[t]
	\centering

	\def0.7\linewidth{0.45\linewidth}
	\input{./collision_modeling/Figures/particle_orientation.pdf_tex}
	\caption{Off-center collision of an ellipsoidal particle with a wall and definition of the particle orientation angle $\theta$.}
	\label{fig:par_orientation}	
\end{figure}

For non-rotating spherical particles colliding with a wall in viscous fluid the rebound only depends on the Stokes number i.e. $e_\mathrm{n}/e_\mathrm{d,n}$ = $f(St_\mathrm{r})$ \citep{gondret2002}. For ellipsoidal particles, even without rotation, other parameters such as the ratio of the particle axes and the particle orientation become pertinent in addition, i.e. $e_\mathrm{n}/e_\mathrm{d,n}$ = $\phi(St_\mathrm{r}, \, a/b, \, b/c, \, \theta/2\pi)$, as illustrated by Fig. \ref{fig:par_orientation}, so that the parameter space is substantially enlarged. If rotation is considered \textcolor{red}{in addition, as is done here}, $e_\mathrm{n}/e_\mathrm{d,n}$ is also a function of the particle angular velocity in case of ellipsoids because the magnitude of the normal velocity of the contact point is a sum of $\vec{u}_\mathrm{p} \cdot \vec{n}$ and $(\vec{\omega}_\mathrm{p} \times \vec{r}_{\mathrm{c}})\cdot \vec{n}$. The latter term vanishes for spherical particles. The orientation of an ellipsoidal particle just before collision, $\theta/2\pi$, and its angular velocity are difficult to control in an experiment. This is one of the reasons why, to the best of the authors knowledge, no experimental data similar to the ones presented for spherical particles by \cite{gondret2002} and \cite{joseph2001} are available with ellipsoids. \\

\begin{figure*}[tp]
	\centering
	\setlength{\unitlength}{1cm}
	\begin{picture}(16,8)


	\put(4,0){\includegraphics{./collision_modeling/Figures/validation_NS.pdf}}
	\end{picture}
	\caption{Ratio of the restitution coefficient of the normal collision,  $e_\mathrm{n}$, to the coefficient for dry collision,  $e_\mathrm{d,n}$, obtained for collisions of an ellipsoidal particle with an inclined wall using the present model. Diamond: Simulations with the present model and the Stokes number fixed at $St_\mathrm{r} \approx 210$, using different initial particle orientations, $\theta_0$. Square: Fixed initial particle orientation, $\theta_0 = \pi/4$, variable Stokes number, as indicated by the horrizontal axis. The experimental data of \citet{schmeeckle2001} for collisions of natural sediment particles with an inclined glass wall are plotted for reference. }
	\label{fig:validation_NS}
\end{figure*}

As a remedy collisions of natural sediment particles are considered. These have been studied intensively in the literature on sediment transport, as this issue constitutes a fundamental aspect of bed-load transport. In one such study, \citet{schmeeckle2001} explored the collision of natural sediment particles with an inclined glass wall in water. They reported a total of $67$ such collisions and observed that a collision can be elastic or inelastic depending upon the orientation of the colliding particle. Only in $18$ collisions of these a rebound of the particle center was noticed. The restitution coefficient, $e_\mathrm{n}$, of the rebounding particles varied between $0.07$ and $0.83$ depending on their orientation just before the collision. The particle orientation $\theta$ ranged between $17\pi/45$ and $\pi/2$ with an average of $\theta = 9\pi/20$. The Stokes number was also different in all these collisions with the majority of them, $12$ out of $18$, happening at $St_\mathrm{r} \approx 210$, as seen in Fig. \ref{fig:validation_NS}. The $49$ particles which did not rebound, hence exhibiting $e_\mathrm{n} = 0$, collided with the wall at a Stokes number between $St_\mathrm{r}=33$ and $St_\mathrm{r}=268$. For these non-rebounding collisions, the average particle orientation was $\theta = 7\pi/20$ in the experiment. The data of these experiments for natural sediment are plotted in Fig. \ref{fig:validation_NS} for later comparison. \\

The physical properties of the particle and the fluid used in the simulations reported here were selected to match the physical properties of natural sediment i.e. quartz, and water (Tab. \ref{tab:NS_collision}). Choosing the particle shape for the simulations requires some selection to narrow down the huge parameter space. The major and minor half-axes of the natural sediment investigated were measured from the photographs by \textcolor{red}{\citet{schmeeckle2001}}, so that the information of only two axes is available, with the ratio of the major axis to the minor axis $b/c$ being essentially in the range of $1.4$ to $2.5$. On this background, the natural sediment particles were represented by oblate ellipsoids, i.e. $a=b$, and an aspect ratio of $b/c = 2$. \textcolor{red}{Beyond the restricted availability of data the choice of an oblate shape was motivated by the measurements of the shape of quartz sand particles conducted by \citet{Smith_Cheung2002}. They found a much higher percentage of oblate-like particles ($42\%$) compared to prolate-like ones ($14\%$), blades ($9\%$), and sphere-like particles ($35\%$), so that an oblate shape seems to be a good choice.} \\

\begin{table*} [tp]
	\centering
	\begin{tabular}{cccclclcc}
		\hline \hline
		$D_\mathrm{eq}\,[\mathrm{mm}]$ & $b/c$ & $\rho_\mathrm{p}\,[\mathrm{kg\,m^{-3}}]$	& $\rho_\mathrm{f}\,[\mathrm{kg\,m^{-3}}]$ & $\nu_\mathrm{f}\,[\mathrm{m^2\,s^{-1}}]$	&  $e_\mathrm{d,n}$&$\mu_\mathrm{s}$ & $\mu_\mathrm{k}$ & $e_\mathrm{d,t}$\\
		\hline
		$6$ & $2$ & $2650$ 	& $998$ & $1.007\times10^{-6}$  & $0.97$	& $0.15$    & $0.15$  & $0.39$ \\
	
		\hline
		\hline
	\end{tabular}
	\caption{Particle and fluid properties employed for the simulations of an oblate particle colliding with an inclined glass wall.}	
	\label{tab:NS_collision}
\end{table*}

\begin{figure*}[tp]
	\centering
	\setlength{\unitlength}{1cm}
	\begin{picture}(16,11)
	\put(0,5.5){\includegraphics{./collision_modeling/Figures/Oblate_traj1.pdf}}
	\put(8,5.5){\includegraphics{./collision_modeling/Figures/Oblate_traj2.pdf}}
	\put(0,0){\includegraphics{./collision_modeling/Figures/Oblate_traj3.pdf}}
	\put(8,0){\includegraphics{./collision_modeling/Figures/Oblate_traj4.pdf}}
	\put( 0.5 , 10.5 )     {a) }
	\put( 8.5 , 10.5 )     {b) }  
	\put( 0.5 , 5. )     {c) }  
	\put( 8.5 , 5. )     {d) }  	
	\end{picture}
	\caption{Oblique collisions of an oblate ellipsoid ($b/c=2$) with a wall at $\phi_\mathrm{in} = \pi/6$. The Stokes number in all cases is $St_\mathrm{r} \approx 210$. Only the initial orientation on the particle was changed. a) $\theta = \pi/2$, b) $\theta = 4\pi/9$, c) $\theta = \pi/3$, and d) $\theta = \pi/4$. The dotted line represents the wall normal position of the particle center and the solid line shows the distance between the particle surface and the wall. Ellipses, not to scale, drawn with time interval of $\SI{0.01}{s}$ visualize the orientation of the particle. }
	\label{fig:oblate_traj}
\end{figure*}

A computational domain of size $\SI{0.04}{m} \times \SI{0.04}{m} \times \SI{0.04}{m}$ discretized with $256 \times 256 \times 256$ cells was used here. Since the majority of the collisions observed in the experiment occurred out at an incidence angle of $\phi_\mathrm{in} \approx \pi/6$, the same value was used in the simulations as well. In a first set of runs, the particle was dropped from a height adjusted in each case so as to reach the desired Stokes number, $St_\mathrm{r} \approx 210$. This height depends on the initial orientation and is about $y_\mathrm{p,0} \approx 1.15\,D_\mathrm{eq}$ (Fig. \ref{fig:oblate_traj}). A few other simulations (not shown here) were conducted with higher $y_\mathrm{p,0}$ but with the same Stokes number, orientation, and angular velocity before collision to understand the effect of the particle wake on the effective restitution coefficient. It was found that the ratio $e_\mathrm{n}/e_\mathrm{d,n}$ changed only very little. This is due to the short time interval between the states $\vec{u}_\mathrm{p,in}$ and $\vec{u}_\mathrm{p,out}$ not providing enough time for the wake to influence this process, even more for the Stokes numbers considered here. The rebound trajectories after the collision were different, though, but this is beyond the collision process to be modeled here and an entire topic on its own.

The simulations were conducted with different initial orientations $\theta_0$ of the particle and $\vec{\omega}_\mathrm{p,0} = 0$. It turned out that during the approach the particle orientation changed only very little, less than five degrees, so that $\theta \approx \theta_0$. The resulting restitution coefficients are reported in Fig. \ref{fig:validation_NS}. If the center of the oblate particle after the collision continues to approach the wall instead of rebounding, as the collision point, the restitution coefficient is negative. In such a case, $e_\mathrm{n}$ was set to zero here, as in the experiment. It was observed in the simulations that the particle does not rebound when $\theta \approx \pi/4$. To investigate this further, a second set of simulations was conducted with $\theta_0 = \pi/4$ varying the gravitational acceleration $g$ to achieve different Stokes numbers. The results are included in Fig. \ref{fig:validation_NS} as well. As a result of the simulations conducted here, it can be seen that the collision model reproduces the phenomenon observed in \cite{schmeeckle2001} very well. \\

The initial and boundary conditions of these simulations are well specified. To provide a reference for later investigations and an illustration of the computed solution four sample trajectories are shown in Fig. \ref{fig:oblate_traj} obtained for a Stokes number $St_\mathrm{r} \approx 210$ and different initial orientations $\theta_0 = [\pi/2, 4\pi/9, \pi/3, \pi/4]$. In every sub-figure, there are two curves. The variation of the distance between the surfaces over time is plotted with a solid line and the wall normal position of the particle center is shown using a dashed line. Additionally, the particle orientation is visualized every $\SI{0.1}{s}$ using ellipses. In Fig. \ref{fig:oblate_traj}a and \ref{fig:oblate_traj}b, a clear rebound of both the particle center as well as the contact point can be noticed. In case of initial orientation $\theta_0=\pi/3$, shown in Fig. \ref{fig:oblate_traj}c, though, the particle center rebounds, but the restitution coefficient is much smaller and equals $0.125$. If $\theta_0=\pi/4$, the center of the particle continues to approach the wall although there is a rebound of the contact point. The center of the particle, then rebounds after the second collision. 

\subsection{Normal collisions of oblate particles with a wall}
To explore the effect of the particle orientation and the ratio of the particle axes further, more simulations were conducted. Here, a setup similar to the experiments of \citet{gondret2002} is used with the parameters listed in Tab. \ref{tab:sensitivity_analysis1}. The same computational domain as in Sec. \ref{ssec:Coll_NS} was employed with $\phi_\mathrm{in} = 0$. The equivalent diameter of the ellipsoid equals $D_\mathrm{eq} = \SI{6}{mm}$ in the simulations. The change in the particle orientation before collision is negligible due to the high density ratio, $\rho_\mathrm{p}/\rho_\mathrm{f} = 7.8$, taken over from the experiments with spheres and the short distance, so that $\theta = \theta_0$. The particle shape was varied by changing the ratio $b/c$ while keeping $a/b=1$ and $D_\mathrm{eq}$ unchanged. Three different ratios were used, $b/c=1$, $b/c=1.5$, and $b/c=2.0$, with increasing flatness $(a+b)/2c$ according to the definition of \cite{Cailleux1945}. 

\begin{figure*}[tp]
	\centering
	\setlength{\unitlength}{1cm}
	\begin{picture}(16,7.5)
	\put(0,0){\includegraphics{./collision_modeling/Figures/validation_oblate.pdf}}
	\put(8,0){\includegraphics{./collision_modeling/Figures/Oblate_e_theta.pdf}}
	\put( 0. , 7.5 )     {a) }
	\put( 8. , 7.5 )     {b) }  
	\end{picture}
	\caption{Normal collision of particles of different shape with a wall showing the dependency of the relative restitution coefficient on the Stokes number, the particle orientation, and the particle shape. a) The Stokes number and the particle shape are varied for a constant orientation. b) The particle orientation and its shape are changed keeping the Stokes number unchanged, $St_\mathrm{r} \approx 75$, indicated by a dashed line in a). The vertical dashed lines in sub-figure b) are positioned at the values $\theta \approx 7\pi/36$ and $\theta \approx 10\pi/36$ mentioned in the text. The symbols give a hint on the particle orientation upon collision and the particle flatness.}
	\label{fig:oblate_flatWall}
\end{figure*}
Two sets of simulations were undertaken for the three shapes considered. One was conducted with $\theta = 0$ varying the Stokes number by changing $\lvert \vec{g}\rvert$. Fig. \ref{fig:oblate_flatWall}a shows the results. It is observed that for oblate particles with $\theta = 0$, in contrast to spherical particles, the relative restitution coefficient, $e_\mathrm{n}/e_\mathrm{d,n}$, increases with increasing Stokes number and saturates at some value below $1$. This maximum value decreases with increasing flatness of the particle, $e_\mathrm{n}/e_\mathrm{d,n}\approx 0.70$ for $b/c=1.5$ and $e_\mathrm{n}/e_\mathrm{d,n}\approx 0.60$ for $b/c=2.0$. This is in line with the expectation, since lubrication forces become larger with increasing radius of curvature. Furthermore, the added mass effect, generating a force opposite to accelerations, changes with particle shape and increases with aspect ratio when an oblate particle is accelerated along its shorter axis. 

Fig. \ref{fig:oblate_flatWall}b shows the change in the restitution coefficient when changing the particle orientation prior to the contact. This is investigated for the three particle shapes examined. Here, the Stokes number was kept constant at a value of about $St_\mathrm{r} \approx 75$. It can be seen that the particle orientation upon contact, $\theta$, has a huge influence on the restitution coefficient. For $b/c=1.5$, the value of $e_\mathrm{n}/e_\mathrm{d,n}$ decreases with increasing $\theta$ up to $\theta \approx 7\pi/36$. It stays roughly constant between $\theta \approx 7\pi/36$ and $\theta \approx \pi/4$ before increasing to its maximum value at $\theta \approx 17\pi/36$. Interestingly, for $\theta$ around $\pi/2$, the restitution coefficient attains values higher than for a spherical particle. This is due to reduced lubrication because of the smaller radius of curvature at the contact point and due to the reduced added mass effects for this orientation of the particle with respect to its trajectory.

The observed trends are similar in the case $b/c=2.0$ and are more pronounced. The restitution coefficient has negative values between $\theta \approx \pi/6$ and $\theta \approx 10\pi/36$. A negative value of $e_\mathrm{n}$ means that the velocity of the particle centre after the collision is negative and the particle center still continues to approach the wall, even if there is a rebound of the contact point. This leads to further collisions with the wall. After two or three rebounds, the particle comes into an enduring contact with the wall still rolling back and forth before eventually coming to rest. 

\textcolor{red}{Let us finally mention that a small number of simulations with prolate particles were conducted in a similar way as reported in Fig. \ref{fig:oblate_flatWall} showing the same qualitative behaviour.}



\subsection{Literature review}
When the distance between two approaching particles becomes small right before the particle surfaces actually touch, the fluid between them is squeezed out of the gap. This process takes place on small length and time scales but is important for the collision as a whole since it influences, e.g., the rebound height due to the related dissipation \citep{Davis1986}. In the IBM framework as used in this study the fluid motion in the liquid film between the particles is resolved as long as the gap between the particles is covered by at least two to three grid cells. When the gap between the particles becomes smaller, the fluid forces acting on the particles are no more resolved and are underpredicted \citep{kempe2012a, izard2014}. This is generally compensated by using an analytical expression for the lubrication force in addition to the fluid forces computed at the immersed boundary. This so-called lubrication model provides the correct inter-particle force, as done in the cited references. Here, the distance where this takes place is denoted $d^*_\mathrm{lub}$ with $d^*_\mathrm{lub} = 2\Delta_x$, as used by \citet{kempe2012a} .

Various lubrication models have been proposed in the literature to model the phases of the collision process with small surface distances before and after the phase of direct surface contact \citep{kempe2012a,  izard2014, biegert2017}. They are all based on the lubrication theory of \citet{brenner1961} and \citet{cox1967} and differ very little with respect to each other. For example, the model for spherical particles used by \citet{kempe2012a} reads
\begin{equation}\label{eq:lub_kempe}
\vec{f}^*_\mathrm{lub} = 
	\begin{cases} 
	0\hspace{0.1cm},   & \qquad \lvert\vec{d}\rvert > d^*_\mathrm{lub} \\
	\displaystyle -\frac{6\pi\mu_\mathrm{f}\vec{u}_{r,n}}{\lvert\vec{d}\rvert}\left( \frac{r_\mathrm{1}r_\mathrm{2}}{r_\mathrm{1}+r_\mathrm{2}}\right) ^2 \hspace{0.1cm}, & \qquad \eta \leqslant \lvert\vec{d} \rvert \leqslant d^*_\mathrm{lub} \\
	0 \hspace{0.1cm},   & \qquad \lvert\vec{d}\rvert < \eta \ \hspace{0.1cm},
	\end{cases}
\end{equation}
where $r_i$ is the radius of the $i$-th particle, $\vec{u}_{r,n}$ the normal component of the relative velocity, and $\mu_\mathrm{f} = \nu_\mathrm{f}\rho_\mathrm{f}$ the dynamic viscosity. The arguments of $\vec{f}^*_\mathrm{lub}$ are the independent variables for given fluid and particle pairing. For distances $\lvert \vec{d} \rvert > d^*_\mathrm{lub}$, the particle motion is resolved by the background Eulerian grid and no modeling of the lubrication forces is applied. Only if the surface distance is between $ \eta \leqslant \lvert\vec{d} \rvert \leqslant d^*_\mathrm{lub}$, the lubrication model sets in. The parameter $\eta$ constitutes a lower bound for the distance in Eq. \eqref{eq:lub_kempe}. It is motivated from a physical point of view by the roughness each and every surface in reality has. It is also needed from an algorithmic point of view since $\vec{f}^*_\mathrm{lub}$ diverges for $\lvert\vec{d} \rvert \to 0$.

The models used by \citet{izard2014}, and \citet{biegert2017} apply a constant force in the range $0 < \lvert\vec{d}\rvert < \eta$ instead of $\vec{f}^*_\mathrm{lub} = 0$ as in Eq. \eqref{eq:lub_kempe}. Moreover, the value of $\eta$ is different in each of these models. In the work of \citet{biegert2017} the parameter, $\eta$, is calibrated such that the rebound trajectory of a $\SI{6}{mm}$ steel particle after collision with a glass wall at Stokes number, $St_\mathrm{r} = 27$ matches with the experiments of \citet{gondret2002}. This results in a value of $\eta$ approximately $30$ times the experimental value $\eta = 10^{-4}D_\mathrm{eq}$ given in \citet{gondret2002} and used by \citet{kempe2012a}. \citet{jeffrey1982} proposed an extended version of the lubrication model which could be applied in the region $ \eta \leqslant \lvert\vec{d} \rvert \leqslant d^*_\mathrm{lub}$. This model includes extra terms to calculate the force. Nevertheless, this force is still inversely proportional to the separation distance.
\begin{figure*}[tp]
	\centering
	\setlength{\unitlength}{1cm}
	\begin{picture}(12,8)
	\put(0,0){\includegraphics{./collision_modeling/Figures/Lub_models_sample.pdf}}
	\end{picture}
	\caption{Lubrication force applied in traditional models compared to the new model proposed here. The value employed for $\eta$, not drawn to scale, is different in the two functions, used by \citet{kempe2012a} and \citet{biegert2017}. In this graph, it is taken to be the same for simplicity. The circles and triangles identify distances assumed in successive time steps in case of a CFL number equal to $0.5$, as discussed in the text.}
	\label{fig:lub_models_sample}
\end{figure*}
The inverse proportionality to the distance $\lvert\vec{d}\rvert$, as in Eq. \eqref{eq:lub_kempe}, leads to a substantial variation in the force applied on the particle depending on where the particle happens to land at a particular time step. The issue is highlighted by the triangles in Fig. \ref{fig:lub_models_sample} corresponding to the evaluation of the lubrication force in case of a constant time step and a constant particle velocity, such that the CFL number is $0.5$. This choice is made for illustration here and can be different in applications where smaller CFL number is usually employed. Similar cases are obtained by shifting the entire set of evaluations by some value in $\lvert \vec{d} \rvert$-direction. The evaluation at discrete times, i.e. distances, samples the curves. It is obvious that if the particle approaches such that the first triangle in the diagram is shifted slightly towards the left, the last triangle will move to a much larger value so that the sum of the forces becomes different if Eq. \eqref{eq:lub_kempe} is used. Furthermore, it is also obvious that the modification for $\lvert\vec{d}\rvert < \eta$, which is not drawn to scale here, does not have much impact, as it makes the body thicker by a tiny amount only ($\eta \approx 10^{-4} D_\mathrm{eq}$) in the \citet{kempe2012a}. As a solution to this problem, \citet{biegert2017} suggested a sub-stepping scheme in the lines of \citet{Kidanemariam2014a} computing the particle motion with smaller timesteps than the fluid motion by dividing a fluid timestep into $150$ substeps, for example.
\begin{figure*}[tp]
	\centering
	\setlength{\unitlength}{1cm}
	\begin{picture}(16,7.5)
	\put(0,0){\includegraphics{./collision_modeling/Figures/Lubrication_model.pdf}}
	\put(8,0){\includegraphics{./collision_modeling/Figures/calibration_constant.pdf}}
	\put( 0. , 7 )     {a) }
	\put( 8. , 7 )     {b) }  
	\end{picture}
	\caption{Calibration of the lubrication model to match the trajectories provided by \cite{gondret2002}. a) Particle trajectories after a normal collision with a wall in viscous fluid obtained in the experiments (symbols) and in the simulations (lines). Indigo: $St_\mathrm{r}=6$; brown: $St_\mathrm{r}=27$; green: $St_\mathrm{r}=100$; black: $St_\mathrm{r}=152$; blue: $St_\mathrm{r}=193$. Here, $t_\mathrm{c}$ is the time of collision. b) The lubrication model parameter $k$ as a function of $St_\mathrm{r}$. The color code of the symbols corresponds to the same values of $St_\mathrm{r}$ as in a). The continuous line is the fit according to Eq. \eqref{eq:k_st}.  }
	\label{fig:lub_calib}
\end{figure*}

However, in a sediment transport simulation, where tens of thousands of particles are colliding in one timestep, this method might prove inefficient and may increase the computational cost. \citet{Kidanemariam2014a} conducted sediment transport simulations with sub-stepping for the particle motion but used a simple collision model without any lubrication force. The coefficients of stiffness and damping in their soft-sphere collision model were calibrated to meet the experiments of \citet{gondret2002} and \citet{joseph2001}. Also, the dry coefficient of restitution was set to a value much smaller than the value for glass which is close to the one for quartz sand employed by \citet{gondret2002} and, thus, relevant for sediment transport. 

\subsection{New lubrication model}
To avoid the uncontrolled variation of $\vec{f}^*_\mathrm{lub}$ just before collision noticed by \citet{biegert2017} without sub-stepping, a different approach is proposed here. It consists of using a constant corrective force for all distances $ 0 < \lvert\vec{d} \rvert \leqslant d_\mathrm{lub} $, as shown in Fig. \ref{fig:lub_models_sample}, where $d_\mathrm{lub} $ is a predefined model parameter. It defines the distance range over which the lubrication force is applied and is much larger than $\eta$ in the models discussed above. The expression of this lubrication force is inspired \eqref{eq:lub_kempe}. This classical model contains factors reflecting geometry, approach velocity, viscosity, inverse distance, and a scalar constant. The first three are retained, $\lvert\vec{d} \rvert$ is replaced with $d_\mathrm{lub}$, and a new dependence on the Stokes number is introduced. The model is given by
\begin{equation}
	\vec{f}_\mathrm{lub} = 
	\begin{cases} 
		0\hspace{0.1cm},   & \qquad \lvert\vec{d}\rvert > d_\mathrm{lub} \\

		\displaystyle - k(St_\mathrm{r})\, \frac{\mu_\mathrm{f}\vec{u}_{r,n}}{d_\mathrm{lub}}\left( \frac{r_\mathrm{1}r_\mathrm{2}}{r_\mathrm{1}+r_\mathrm{2}}\right) ^2\hspace{0.1cm}, & \qquad \lvert\vec{d} \rvert \leqslant d_\mathrm{lub} \hspace{0.1cm}.
	\end{cases}
	\label{eq:lub_model}
\end{equation}
\cite{kempe2012a} argued that the fluid forces influencing the particle motion are well resolved as long as the distance between two particles, or between a particle and a wall is greater than $2\Delta x$. Following this recommendation, others have also applied the lubrication force in the region $ 0 \leqslant \lvert\vec{d} \rvert \leqslant 2\Delta_x$ \citep{izard2014,Ardekani2016,biegert2017}. The same value is used here setting $d_\mathrm{lub}=2\Delta x$. Furthermore, it will be reported in Sec. \ref{sec:sensitivity_analysis} below that the model is fairly insensitive to this parameter. 

The function $k(St_\mathrm{r})$ is determined such that the rebound trajectories match the experimental data of \citet{gondret2002}. The trajectories available in this reference, except those at very large $St_\mathrm{r}$, were used to calibrate the present model. This is accomplished using the relative Stokes number
\begin{equation}
	St_\mathrm{r} = \frac{1}{9}\frac{\rho_\mathrm{p}}{\rho_\mathrm{f}}\frac{\lvert \vec{u}_{r,n} \rvert}{\nu_\mathrm{f}} \left(\frac{r_\mathrm{1}r_\mathrm{2}}{(r_\mathrm{1}+r_\mathrm{2})/2}\right) \hspace{0.1cm}.
	\label{eq:rel_stokes}
\end{equation}
\textcolor{red}{Unlike the classical Stokes number for a single particle in a flow, which is based on the velocity difference between particle and fluid, this Stokes number is based on the difference between the velocities of the two particles, more precisely the component in the direction of the distance vector. This velocity scale was selected because it is the relevant scale for the lubrication forces. Also, the length scale $r_\mathrm{1}r_\mathrm{2}/((r_\mathrm{1}+r_\mathrm{2})/2)$, which is $r$ in the case of equal radii $r_\mathrm{1}=r_\mathrm{2}=r$, seems appropriate to address lubrication forces. In fact these scales already appear in the classical lubrication models looking like Eq. \eqref{eq:lub_kempe}. Finally, this choice is backed \textit{a posteriori} by the success of this parametrization achieved in covering the experimental data.} The dependence of the collision force on $St_\mathrm{r}$ is a new feature compared to previous models and introduces the particle density in the collision model.

The experiments of \citet{gondret2002} were conducted with spherical particles of different materials impacting on a glass wall in different fluids. In the present study a wall is represented as a sphere of radius $\SI{e+12}{m}$  and density $\SI{e+12}{kg/m^3}$. Replacing $r_\mathrm{1}$ with the radius of the colliding particle and $r_\mathrm{2}$ with $\SI{e+12}{m}$ in Eq. \eqref{eq:rel_stokes} yields the same definition of the Stokes number as used in the cited experiments. 
\subsection{Dependency of the new lubrication model on the Stokes number}

To provide a closed form expression of $k(St_\mathrm{r})$ to be used in Eq. \eqref{eq:lub_model}, a set of simulations of particle-wall collision were carried out. The physical parameters of particles and fluids as well as the dimensionless numbers characterizing the simulations are assembled in Tab. \ref{tab:lub_calibration}. The computational domain of size $6.6D_\mathrm{p}\times13.4D_\mathrm{p}\times6.6D_\mathrm{p}$ was discretized with $128\times256\times128$ cells in all the simulations. This gives a spatial resolution of $D_\mathrm{p}/\Delta_x=19$. Periodic boundary conditions were imposed in the wall-parallel directions and a no-slip and a free slip condition on the bottom and the top boundary, respectively. A fixed time step was used corresponding to $\mathrm{CFL}=0.1$ based on the terminal velocity of the particle. The initial condition was still fluid and a particle position at a vertical distance of $10D_\mathrm{p}$ above the wall.
\begin{table*} [htp]
	\centering
	\begin{tabular}{cccccllcl}
		\hline \hline
		Case	& $D_\mathrm{p}\,[\mathrm{m}]$	& $\rho_\mathrm{p}\,[\mathrm{kg\,m^{-3}}]$	& $\rho_\mathrm{f}\,[\mathrm{kg\,m^{-3}}]$ & $\nu_\mathrm{f}\,[\mathrm{m^2\,s^{-1}}]$	& $St_\mathrm{r}$ & $Re_\mathrm{p}$ & $e_\mathrm{d,n}$ 	& $k$\\
		\hline
		$1$ 	& $3\times10^{-3}$ 	& $7800$ 	& $965$ & $1.0363\times10^{-4}$ 	& $6$    & $6$   & $ 0.97 $ & $125$\\
		$2$ 	& $6\times10^{-3}$ 	& $7800$ 	& $965$ & $1.0363\times10^{-4}$ 	& $27$   & $30$  & $ 0.97 $ & $120$\\
		$3$ 	& $4\times10^{-3}$ 	& $7800$ 	& $953$ & $2.0986\times10^{-5}$ 	& $100$  & $110$ & $ 0.97 $ & $78$ \\
		$4$ 	& $3\times10^{-3}$ 	& $7800$ 	& $935$ & $1.0695\times10^{-5}$ 	& $152$  & $165$ & $ 0.97 $ & $37$ \\
		$5$ 	& $6\times10^{-3}$ 	& $7800$ 	& $953$ & $2.0986\times10^{-5}$ 	& $193$  & $212$ & $ 0.97 $ & $18$ \\
		\hline
		\hline
	\end{tabular}
	\caption{Material properties of particles and fluids used in the experiments of \citet{gondret2002} and in the present simulations. The values of the model parameter $k$ leading to matching of trajectories between simulation and experiment are listed as well. The particle Reynolds number, $Re_\mathrm{p}$, is defined using the particle velocity before collision and the particle diameter $D_\mathrm{p}$. }	
	\label{tab:lub_calibration}
\end{table*}

The experiments of \citet{gondret2002} were conducted such that the settling particle reached its terminal velocity before colliding. This would require a large enough computational domain or an artificial acceleration \citep{biegert2017}. In the simulations reported here, a moving frame methodology was used to accelerate the particle to its terminal velocity. Once the particle Stokes number of the experiment was attained, the moving frame was fixed and the particle was allowed to freely collide with the wall. 

The lubrication force was applied as discussed above. Each simulation was repeated varying the value of $k$ value until the rebound trajectory matched the experiment. The resulting trajectories are shown in Fig. \ref{fig:lub_calib}a. The value of $k$ resulting from this procedure is provided in Tab. \ref{tab:lub_calibration} and plotted in Fig. \ref{fig:lub_calib}b. It is found that $k$ decreases substantially with $St_\mathrm{r}$. \citet{izard2014} demonstrated that the effect of the lubrication force is negligible for $St_\mathrm{r} \geqslant 200$, supporting the present observation and indicating that the decreasing trend continues for values of $St_\mathrm{r}$ larger than the ones selected here. It turns out that the behavior of $k(St_\mathrm{r})$ can be fitted extremely well by a simple Gauss function which is also shown in Fig. \ref{fig:lub_calib}b. This exponential function is motivated by the correlation between Stokes number and the effective coefficient of restitution $e_\mathrm{n}$ conceptualized by \citet{izard2014}. These authors proposed a model to reproduce the dependency of the effective restitution coefficient on the Stokes number observed in experiments and numerical studies. The function used here is  
\begin{equation} {\label{eq:k_st}}
k(St_\mathrm{r}) = \alpha \exp\left( -\frac{1}{2\sigma^2} St_\mathrm{r}^2\right) \hspace{0.1cm}, 
\end{equation}
with $ \alpha = 125$, and $ \sigma = 100 $ obtained from the fit to the simulation data. 
\subsection{Extension to non-spherical particles}
\textcolor{red}{To model the lubrication force acting on non-spherical particles, particular geometrical properties of the contacting surfaces need to be considered in the model,  such as the two principal radii of curvature at the contact point \citep{Cox1974}. \cite{Janoschek2013}, proposed an accurate method to compute the lubrication corrections in a case of two approaching spheroids using a third-order description of the particular particle surfaces which requires the determination of the associated coefficients at each contact point. In a simulation with thousands of colliding ellipsoids, however, such an algorithm is computationally too expensive. As an alternative, \citet{Ardekani2016} proposed an expression of the lubrication force based on the Gaussian curvature of the spheroidal particles at the contact point, where the local particle surface is approximated by a sphere of the same curvature. Since this approach provides a good balance between accuracy and computational efficiency it is employed in the present work as well, just for the more general case of ellipsoids instead of spheroids.} The Gaussian curvature of an ellipsoid at a point $ (X,Y,Z) $ in a body-fixed frame is 
\begin{equation}
	\mathrm{G} = \displaystyle \frac{1}{a^2b^2c^2\left[ \frac{X^2}{a^4} + \frac{Y^2}{b^4} +\frac{Z^2}{c^4} \right]^2 } \hspace{0.1cm},
\end{equation}
and the radius of Gaussian curvature is given by \citep{gray1997}
\begin{equation}
R_\mathrm{G} = \frac{1}{\sqrt{\mathrm{G}}} \hspace{0.1cm}.
\label{eq:Rg}
\end{equation}
Along these lines the lubrication force between two non-spherical particles is determined here using Eq. \eqref{eq:lub_model} with the radius of the spherical particles, $r_\mathrm{1}$ and $r_\mathrm{2}$ replaced by $R_\mathrm{G,1}$ and $R_\mathrm{G,2}$ according to Eq. \eqref{eq:Rg}, respectively, which depend on the particular collision point. The Stokes number in Eq. \eqref{eq:rel_stokes} is also determined in this way, and the expression for $k$ in Eq. \eqref{eq:k_st} is employed with the same values of $\alpha$ and $\sigma$ as used for spherical particles.
\subsection{Single collision}
\label{sec:single_collision}
Once it is found that the two ellipsoids are in contact, the collision response of both particles needs to be computed. The associated collision model used for this purpose is based on the model developed by \citet{Tschisgale2018b} for the collision of rod structures in dry environment. In the present work, this algorithm is enhanced to model the collision of rigid ellipsoids colliding in viscous fluid. This section describes the basic ideas of the collision modeling and the derivation of the collision forces and torques, $\vec{f}_{\mathrm{c}}$ and $\vec{m}_{\mathrm{c}}$, respectively. 
Staring point of the derivation is the configuration of two colliding ellipsoids at the moment of their first contact, shown in Fig. \ref{fig:colliding_pair}. 

\begin{figure}
	\centering

	\def0.7\linewidth{0.7\linewidth}
	\input{./collision_modeling/Figures/collision.pdf_tex}
	\caption{Two ellipsoids in collision and the variables defining the system. The time index is dropped here for better readability. }
	\label{fig:colliding_pair}	
\end{figure}

To describe the time-discrete equations of motion of the particles, it is assumed that the instant shown in this figure coincides with the discrete time level $t^{n-1}$. Due to the high rigidity of the particles considered here, the timescale of contact is much smaller than the long-term behavior of the particles. Especially for perfectly rigid bodies, the collision time is infinitesimally small, so that at the moment of first touch, the velocities of the particles discontinuously change while the particle position and orientation remain unchanged. This implies that the collision force at the contact point has a Dirac-like impulsive nature. In the present time-discrete context, this collision force is approximated by a constant value acting over the time interval $\Delta t$ denoted $\vec{f}_{\mathrm{c}}$. The velocities of the particles after the collision process are then given at the new time level $t^{n}$.\\
The force $\vec{f}_{\mathrm{c}}$ can be derived taking advantage of the velocities at the contact point of both ellipsoids $\vec{u}_{\mathrm{c,1}}$ and $\vec{u}_{\mathrm{c,2}}$. After the collision, at time $t^{n}$, these are
\begin{equation}{\label{eq:velAtContact}}
\vec{u}^{n}_{\mathrm{c}} = \vec{u}^{n}_\mathrm{p} + \vec{\omega}^{n}_\mathrm{p} \times \vec{r}_{\mathrm{c}} \hspace{0.1cm},
\end{equation}
with the vector, $\vec{r}_{\mathrm{c}}$ joining the contact point $\vec{x}_\mathrm{c}$ and the particle center $\vec{x}_\mathrm{p}$. The index identifying the particle being dropped for ease of notation. In Eq. \eqref{eq:velAtContact}, the linear and angular velocity $\vec{u}^{n}_\mathrm{p}$ and $\vec{\omega}^{n}_\mathrm{p}$ can be replaced by utilizing the discretized equations of motion, \eqref{eq:lin_par} and \eqref{eq:ang_par}, so that
\begin{equation}{\label{eq:longlonglilong}}
\begin{split}
\vec{u}^{n}_{\mathrm{c}} & =  \textbf{u}_\mathrm{p}^{n-1} + \: m_\mathrm{t}^{-1} \: \Delta{t} \: \left\{ \vec{f}_\mathrm{f} + \vec{f}_\mathrm{g} + \vec{f}_\mathrm{c} + \vec{f}_\mathrm{lub} \right\} \\
& + \left[ \textbf{q}_\mathrm{p}^{n-1} \circ \left( \boldsymbol{\Omega}_\mathrm{p}^{n-1} + \Delta{t} \: \textbf{I}_\mathrm{t} ^{-1} \! \cdot \! \left\{  \vec{M}_\mathrm{f}  + \vec{M}_\mathrm{c} + \vec{M}_\mathrm{lub} \right\} - \textbf{I}_\mathrm{t} ^{-1} \int^{n}_{n-1} \!\! \boldsymbol{\Omega}_\mathrm{p} \times \textbf{I}_\mathrm{p} \!\! \cdot \! \boldsymbol{\Omega}_\mathrm{p} \: \text{d} t \right) \circ \bar{\textbf{q}}_\mathrm{p}^{n-1} \right] \times \vec{r}_{\mathrm{c}} \hspace{0.1cm},
\end{split}
\end{equation}
with $m_\mathrm{t} = m_\mathrm{p} + m_{L}$ and $\textbf{I}_\mathrm{t} = \textbf{I}_\mathrm{p} + \textbf{I}_{L}$. After rearranging, this equation can be expressed in more compact form as
\begin{equation}{\label{eq:afterRearrange}}
\vec{u}^{n}_{\mathrm{c}} = \vec{u}^{n-1}_{\mathrm{c}} + \vec{u}_{\mathrm{ex}} + \: m_\mathrm{t}^{-1} \: \vec{p}_\mathrm{c} 
+ \: \textbf{i}_\mathrm{t} ^{-1} \! \cdot \! \vec{l}_\mathrm{c} \: \times \vec{r}_{\mathrm{c}} \hspace{0.1cm},
\end{equation}
where $\vec{u}^{n-1}_{\mathrm{c}}$ is the velocity at the contact point before the collision. Furthermore, $\vec{p}_\mathrm{c} = \Delta t \, \vec{f}_\mathrm{c}$ is the linear momentum due to the collision and $\vec{l}_\mathrm{c} = \Delta t \, \vec{m}_\mathrm{c}$ the angular momentum. The inertia matrix in the global coordinate system is obtained by the transformation $\textbf{i}_\mathrm{t} = \mathcal{R} \cdot \textbf{I}_\mathrm{t} \cdot \mathcal{R}^{\mathrm{T}}$, with the rotation matrix $\mathcal{R}$ converting any vector $\vec{s}$ from the body-fixed frame into the global coordinate system, i.e. $\mathcal{R} \cdot \vec{s} = \textbf{q}_\mathrm{p} \circ \vec{s} \circ \bar{\textbf{q}}_\mathrm{p}$. Finally, comparing Eqs. \eqref{eq:afterRearrange} and \eqref{eq:longlonglilong} shows that
\begin{equation}{\label{eq:externalForce}}
\begin{split}
\vec{u}_{\mathrm{ex}}  & = \Delta{t} \: m_\mathrm{t}^{-1} \: \left\{ \vec{f}_\mathrm{f} + \vec{f}_\mathrm{g} + \vec{f}_\mathrm{lub} \right\} \\
                       & + \left[ \textbf{q}_\mathrm{p}^{n-1} \circ  \Delta{t} \: \textbf{I}_\mathrm{t} ^{-1} \! \cdot \! \left\{  \vec{M}_\mathrm{f}  + \vec{M}_\mathrm{lub} \right\}  \circ \bar{\textbf{q}}_\mathrm{p}^{n-1} \right] \times \vec{r}_{\mathrm{c}} \\
                       & - \left[ \textbf{q}_\mathrm{p}^{n-1} \circ \left( \textbf{I}_\mathrm{t} ^{-1} \int^{n}_{n-1} \!\! \boldsymbol{\Omega}_\mathrm{p} \times \textbf{I}_\mathrm{p} \!\! \cdot \! \boldsymbol{\Omega}_\mathrm{p} \: \text{d} t \right) \circ \bar{\textbf{q}}_\mathrm{p}^{n-1} \right] \times \vec{r}_{\mathrm{c}}
\end{split}
\end{equation}
which can be computed fully explicitly at $t^{n-1}$ before the collision since $\vec{f}_\mathrm{f}$ and $\vec{f}_\mathrm{g}$ are known from \eqref{eq:f_fandf_g}, provided that $\vec{f}_\mathrm{lub}$ and $\vec{M}_\mathrm{lub}$ are known. Such explicit expressions will be provided below. The quantity $\vec{u}_\mathrm{ex}$ in \eqref{eq:externalForce} represents the velocity induced by external forces including hydrodynamic, gravitational and lubrication forces acting during the collision time interval. It does not have a time index since $\vec{u}_{\mathrm{ex}}$ contains forces that have intermediate values and can not be associated to one time step. Note that the last term of Eq. \eqref{eq:externalForce} vanishes since
\begin{equation}
  \boldsymbol{\Omega}_\mathrm{p} \times \textbf{I}_\mathrm{p} \!\! \cdot \! \boldsymbol{\Omega}_\mathrm{p} = \bar{\textbf{q}}_\mathrm{p} \circ \frac{\text{d}\, \textbf{i}_\mathrm{p}}{\text{d}t} \vec{\omega}_\mathrm{p} \circ \textbf{q}_\mathrm{p}
\end{equation}
and ${\text{d}\, \textbf{i}_\mathrm{p}}/{\text{d}t} \to 0$ due to the assumption of infinitesimally short collision time \citep{Guendelman2003, Tschisgale2018b}.\\ 
In Eq. \eqref{eq:afterRearrange}, the angular momentum in the global Eulerian frame, $\vec{l}_\mathrm{c}$, is related to the linear momentum $\vec{p}_\mathrm{c}$ by 
$\vec{l}_\mathrm{c} = \vec{r}_{\mathrm{c}} \times \vec{p}_\mathrm{c} $ since there is only one contact point for colliding convex rigid bodies such as ellipsoids. This property allows to express Eq. \eqref{eq:afterRearrange} in the form
\begin{equation}{\label{eq:relVelFinal}}
\vec{u}^{n}_{\mathrm{c}} = \vec{u}^{n-1}_{\mathrm{c}} + \vec{u}_{\mathrm{ex}} + \: \vec{K} \cdot \vec{p}_\mathrm{c} \hspace{0.1cm},
\end{equation}
with the symmetric system matrix
\begin{equation}
\vec{K} =  m_\mathrm{t}^{-1}\mathbb{I} + [\vec{r}_{\mathrm{c}}]^\mathrm{T}_\mathrm{x} \cdot \textbf{i}_\mathrm{t} ^{-1} \cdot [\vec{r}_{\mathrm{c}}]_\mathrm{x} \hspace{0.1cm}.
\end{equation}
The term $[\vec{r}_{\mathrm{c}}]_\mathrm{x}$ is the skew matrix of the vector $\vec{r}_{\mathrm{c}}$ and is defined by the equality
\begin{equation}
[\vec{r}_{\mathrm{c}}]_\mathrm{x} \cdot \textbf{p}_\mathrm{c} = \vec{r}_{\mathrm{c}} \times \textbf{p}_\mathrm{c} \hspace{0.1cm}.
\end{equation}

The quantities in Eq. \eqref{eq:relVelFinal} differ for both colliding particles, except the shared collision load acting at the common contact point, so that
\begin{subequations}
	\begin{align}
	\label{eq:relVelFinal1}
\vec{u}^{n}_{\mathrm{c,1}} &= \vec{u}^{n-1}_{\mathrm{c,1}} + \vec{u}_{\mathrm{ex,1}} - \: \vec{K}_1 \cdot \vec{p}_\mathrm{c} \hspace{0.1cm} \\
  \label{eq:relVelFinal2}
\vec{u}^{n}_{\mathrm{c,2}} &= \vec{u}^{n-1}_{\mathrm{c,2}} + \vec{u}_{\mathrm{ex,2}} + \: \vec{K}_2 \cdot \vec{p}_\mathrm{c} \hspace{0.1cm}
	\end{align}
\end{subequations}
Here, the linear momentum transferred in a collision, $\vec{p}_\mathrm{c}$, has the same magnitude in both equations but is directed in opposite direction according to the principle of {\it actio et reactio}.
This allows to derive the collision load $\vec{p}_\mathrm{c}$ by subtracting Eq. \eqref{eq:relVelFinal2} from Eq. \eqref{eq:relVelFinal1} resulting in 
\begin{equation}{\label{eq:relVelFinal3}}
\vec{p}_\mathrm{c} = \vec{K}_{12}^{-1} \cdot (\vec{u}^{n}_{\mathrm{r}} - \vec{u}^{n-1}_{\mathrm{r}} - \vec{u}_{\mathrm{ex,r}}) \hspace{0.1cm},
\end{equation}
with the total system matrix $\vec{K}_{12} = \vec{K}_1 + \vec{K}_2$ and $\vec{u}^{n}_{\mathrm{r}} = \vec{u}^{n}_{\mathrm{c,2}}-\vec{u}^{n}_{\mathrm{c,1}}$ the relative velocity of the particles at the contact point after the collision at time $t^n$. Since $\vec{u}^{n}_{\mathrm{r}}$ is still unknown, the Poisson hypothesis \citep{Lubarda2009} is used here to close this equation.  
It relates the normal relative velocity $\vec{u}^n_\mathrm{r,n}$ before and after the collision via
\begin{equation}\label{eq:edn}
\vec{u}^n_\mathrm{r,n} = -e_\mathrm{d,n}\, \vec{u}^{n-1}_\mathrm{r,n} = -e_\mathrm{d,n}\, (\vec{u}^{n-1}_\mathrm{r}\cdot\vec{n}) \, \vec{n} \hspace{0.1cm} ,
\end{equation}
where $e_\mathrm{d,n}$ is the dry coefficient of restitution.
Analogously, the tangential velocities can also be related through a tangential restitution coefficient, $e_\mathrm{d,t}$, using
\begin{equation}\label{eq:edt}
\vec{u}^n_\mathrm{r,t} = -e_\mathrm{d,t}\, \vec{u}^{n-1}_\mathrm{r,t} = -e_\mathrm{d,t}\, (\vec{u}^{n-1}_\mathrm{r}\cdot\vec{t}) \, \vec{t} \hspace{0.1cm} 
\end{equation}
and the decomposition $\vec{u}^{n}_{\mathrm{r}} = \vec{u}^n_\mathrm{r,n}+\vec{u}^n_\mathrm{r,t}$. The unit tangent vector, $\vec{t}$, indicating the tangential projection of the particle motion is defined as $\vec{t}$ = ${(\vec{u}^{n-1}_\mathrm{r} - \vec{u}^{n-1}_\mathrm{r,n})}/{\lvert (\vec{u}^{n-1}_\mathrm{r} - \vec{u}^{n-1}_\mathrm{r,n})\rvert}$. 

Eqs. \eqref{eq:edn} and \eqref{eq:edt} allow to calculate the momentum transfer $\vec{p}_\mathrm{c}$ in Eq. \eqref{eq:relVelFinal3} from the initial state of the colliding particles at first touch, i.e. at time level $t^{n-1}$.
However, the tangential part of the relative velocity after the collision, $\vec{u}^{n}_\mathrm{r,t}$, remains an unknown quantity, since it depends upon whether the particles slide or stick to each other during the contact. In case of sticking, $\vec{u}^{n}_\mathrm{r,t}$ can be calculated using \eqref{eq:edt}, given that $e_\mathrm{d,t}$ is known {\em a priori}. In case of sliding, $\vec{u}^{n}_\mathrm{r,t}$ must be provided in a different way. 

In the present work, this is done by following the recommendations of \citet{Guendelman2003}.
As a first guess, it is assumed that the colliding particles stick fully at the contact point in tangential direction, i.e.~$e_\mathrm{d,t} = 0$. 
Using the Poisson hypothesis in normal direction with $e_\mathrm{d,n}$, the change in the relative velocity during a collision is
\begin{equation}	{\label{eq:Delta_u}}
	\Delta{\vec{u}} = \vec{u}^{n-1}_\mathrm{r} - \vec{u}^n_\mathrm{r} = \vec{u}^{n-1}_\mathrm{r} +\, e_\mathrm{d,n}\, (\vec{u}^{n-1}_\mathrm{r}\cdot\vec{n}) \, \vec{n} \hspace{0.1cm} ,
\end{equation}
so that Eq. \eqref{eq:relVelFinal3} turns into
\begin{equation}\label{eq:pc_stick}
    \vec{p}_\mathrm{c} = -\vec{K}_{12}^{-1} \cdot \Delta\vec{u} \; -\vec{K}_{12}^{-1} \cdot \vec{u}_{\mathrm{ex,r}} \hspace{0.1cm} .
\end{equation}
Here, the first term on the right hand side imposes the required velocity constraint of the Poisson hypothesis, commonly used in dry hard sphere collision models. The second term captures the momentum transferred by external loads during the collision process, and is usually not considered by common hard sphere models. 

In a second step of the collision modeling, $\vec{p}_\mathrm{c}$ is used to verify whether the particles stick or slide along each other during contact. 
According to the Coulomb friction model, the particle surfaces stick if $\lvert \vec{p}_\mathrm{c} \cdot \vec{t} \rvert \le \mu_s \, \lvert \vec{p}_\mathrm{c} \cdot \vec{n} \rvert$ with $\mu_s$ being the static coefficient of friction. Once this condition is fulfilled, $\vec{p}_\mathrm{c}$ according to Eq. \eqref{eq:pc_stick} is a good approximation of the collision response. Otherwise, $\vec{p}_\mathrm{c}$ has to be modified to allow sliding in tangential direction. Based on the Coulomb law of friction, $\vec{p}_\mathrm{c}$ during sliding can be expressed in the local coordinate system as
\begin{subequations}
	\begin{align}
	\vec{p}_\mathrm{c} &= p_\mathrm{n}\, \vec{n} + p_\mathrm{t} \,\vec{t}\\
		\label{eq:pc_sliding_2}
	&= p_\mathrm{n}\,(\vec{n} + \mu_\mathrm{k}\,\vec{t}) \, ,
	\end{align}
\end{subequations}
with the normal and tangential components $p_\mathrm{n}$ and $p_\mathrm{t}$, respectively. The factor $\mu_k$ is the kinetic coefficient of friction. It relates the tangential and normal part of the collision response via $p_\mathrm{t} = \mu_\mathrm{k} \, p_\mathrm{n}$.
The alteration of the tangential part of $\vec{p}_\mathrm{c}$ in Eq. \eqref{eq:pc_sliding_2} violates the Poisson hypothesis in normal direction. This is because $\vec{p}_\mathrm{c}$ in Eq. \eqref{eq:pc_stick} is derived by assuming that the particles stick to each other at the contact point, i.e. $\vec{u}^n_\mathrm{r,t} = 0$. Therefore, the normal momentum $p_\mathrm{n}$ in case of sliding has to be recalculated by inserting $\vec{p}_\mathrm{c}$ from Eq. \eqref{eq:pc_sliding_2} in Eq. \eqref{eq:pc_stick}
\begin{equation}
\left[\vec{K}_{12} \cdot p_\mathrm{n}\,(\vec{n} + \mu_\mathrm{k}\,\vec{t})\right] \cdot \vec{n} = \left[-\Delta\vec{u} -\vec{u}_{\mathrm{ex,r}}\right] \cdot \vec{n} \:  \hspace{0.1cm},
\end{equation}
so that
\begin{equation}
p_\mathrm{n} = - \; \frac{ \Delta\vec{u} \!\cdot\! \vec{n} \, + \, \vec{u}_{\mathrm{ex,r}} \!\cdot\! \vec{n} }{ \vec{K}_{12} \cdot \,(\vec{n} + \mu_\mathrm{k}\,\vec{t}) \cdot \vec{n} } \hspace{0.1cm} .
\end{equation}
Finally, the vectorial momentum transfer during a collision in case of tangential sliding is obtained from Eq. \eqref{eq:pc_sliding_2} and the associated torque via $\vec{l}_\mathrm{c} = \vec{r}_{\mathrm{c}} \times \vec{p}_\mathrm{c} $.

\subsection{Multiple simultaneous collisions}
\label{LCP}
When simulating sediment transport with thousands of particles, almost each single particle is in contact with several other particles at the same collision time. As a result, the collisional momentums, $\vec{p}_\mathrm{c}$ and $\vec{l}_\mathrm{c}$, must be an outcome of all the collisions. Using the approach in Section \ref{sec:single_collision}, the collision loads can only be calculated for a single collision between two particles. The derivation of multiple simultaneous collisions differs from the above and is considerably more extensive. To overcome this issue, a common assumption is to suppose that all particles collide in pairs, i.e. a single collision of two particles, one after another. This is also done here. Each single collision is modeled exactly as described above. Then, all single collisions are executed multiple times in an iterative procedure until $\vec{p}_\mathrm{c}$, acting between each single pair, prevents an interpenetration of the particles. This is the Signorini condition, i.e. $\vec{p}_\mathrm{c} \cdot \vec{n} \ge 0$. In each iteration step, the order of the collision pairs is varied at random to achieve the most realistic collision response of the sediment.
More details on the collision modeling of multiple simultaneous collisions can be found in \cite{Tschisgale2018b} including a pseudocode of the iterative procedure.

