Local lubrication model for spherical particles within incompressible Navier-Stokes flows

The lubrication forces are short-range hydrodynamic interactions essential to describe suspension of the particles. Usually, they are underestimated in direct numerical simulations of particle-laden flows. In this paper, we propose a lubrication model for a coupled volume penalization method and discrete element method solver that estimates the unresolved hydrodynamic forces and torques in an incompressible Navier-Stokes flow. Corrections are made locally on the surface of the interacting particles without any assumption on the global particle shape. The numerical model has been validated against experimental data and performs as well as existing numerical models that are limited to spherical particles.

Figure 1

Sketch of two particles in interaction.

Figure 2

Contact of two particles with notations associated to the soft sphere model.

Figure 3

Sketch of the decomposition for the domain $D$ into a fluid domain $D_{\text{f}}$ and the solid particles $D_{\text{s}}$.

Figure 4

Sketch of a grid cell with the notations used for the spatial discretization of the governing equations.

Figure 5

Two-dimensional sketch of the IPC method. Velocity at the solid ghost point G is corrected using its symmetric F. The velocity at F is obtained by interpolation of the velocity of its four neighboring fluid points (eight neighbors in three-dimensional cases). The orthogonal projection of G at the particle surface is denoted B.

Figure 6

Representation of the two parameters ${\epsilon }_{\text{lub}}$ and ${\epsilon }_{\text{col}}$ of the lubrication correction model. Typically, $a{\epsilon }_{\text{lub}}\sim 2–3$ grid cells.

Figure 7

Sketch description of the numerical algorithm used to compute the state $n+1$ of the whole system from its state $n$. The algorithm starts at the estimation of the time step $\mathrm{\Delta }t$ which is the time elapse between the system state $n$ and $n+1$. The particle dynamics is solved in $n_t$ sub-time steps $\delta t=\frac{\mathrm{\Delta }t}{n_t}$. $m$ denotes the current substate of the system between the state $n$ and $n+1$ while the particle dynamics is computed. Dashed boxes contain the reference to the main equations computed at the given step.

Figure 8

Sketch of the two interacting particles with the notations used to evaluate ${\mathbf{F}}_{i,j}^{\text{lub}}$ and ${\mathbf{T}}_{i,j}^{\text{lub}}$.

Figure 9

Sketch of two cross-sections of the domain with its configuration and initial location of the particle.

Figure 10

Evolution of the vertical velocity of the particle versus the nondimensional gap size. Simulations using local (LLCM) and tabulated (CLM) lubrication correction are compared to Harada [38] experimental measurements and the model H. The LLCM and CLM are activated for $\epsilon$ smaller than the blue square and red star, respectively.

Figure 11

Total hydrodynamic force as a function of $\epsilon$ in the case of a single particle approaching a solid wall. Simulations results are compared against the analytical solution of Brenner [39] (dashed line). All simulations are run with $h=1/40$. The LLCM is activated on all $\epsilon$ smaller than ${\epsilon }_{\text{lub}}=2\frac{\mathrm{\Delta }y}{a}$.

Figure 12

Total hydrodynamic force, according to the $y$ direction as a function of $1/\epsilon$ during the approach phase. For each grid mesh resolution $h$, lubrication correction is activated for $\epsilon \le {\epsilon }_h$.

Figure 13

Total hydrodynamic force, according to the $y$ direction as a function of $1/\epsilon$ during the approach phase. For all curves, the lubrication correction is activated for $\epsilon \le {\epsilon }_{\text{lub}}$.

Figure 14

Representation of the global rate of convergence on velocity (a) and position (b) of the particle center of gravity according to the norms $L^2$ and $L^{\infty }$. The particle position $Y_h$ and velocity $U_h$ obtained with a grid mesh resolution $h$ are compared to the particle position $Y_{1/100}$ obtained with a grid mesh resolution $h=1/100$.

Figure 15

Total hydrodynamic force, according to the $y$ direction as a function of $1/\epsilon$ during the approach phase for the critical lubrication distances $a{\epsilon }_{\text{lub}}$ equal to $0.1\mathrm{\Delta }y,1\mathrm{\Delta }y,2\mathrm{\Delta }y,20\mathrm{\Delta }y,3\mathrm{\Delta }y$, and $5\mathrm{\Delta }y$ (by increasing $F/U$ at the highest $1/\epsilon$). Grid mesh resolution is $h=1/40$ for all curves.

Figure 16

Evolution of the total energy (energy potential and kinetic) of the system during an elastic (${\xi }_{\text{max,n}}=1$) or plastic (${\xi }_{\text{max,n}}=0.5$) collision of a steel particle with a wall. In both cases the impact Stokes number is ${\text{St}}_d\approx 6900$.

Figure 17

Distribution of normalized effective coefficient of normal restitution ${\xi }_{\text{n}}/{\xi }_{\text{max,n}}$ of a single particle impacting a wall in respect of the particle Stokes number at the impact. Filled markers represent results obtained using the LLCM using the same experimental setup than Joseph et al. [44] (data represented by hallow markers). The black curve represents a correlation made on experimental data proposed by Legendre et al. [45]

Figure 18

Sketch of the simulation domain and initial location of the particle. Inside the dashed box, particle characteristic velocities and angles are displayed before and after collision with the wall.

Figure 19

Comparison between the normalized incidence ${\mathrm{\Psi }}_{\text{in}}$ and rebound ${\mathrm{\Psi }}_{\text{out}}$ angles of steel and glass particles in aqueous solution (a) or in the air (b). Results of simulations are represented by filled dots and are compared to Joseph and Hunt [46] experimental measurements (by circles).