Particle–boundary interaction in a shear-driven cavity flow

The motion of a heavy finite-size tracer is numerically calculated in a two-dimensional shear-driven cavity. The particle motion is computed using a discontinuous Galerkin-finite-element method combined with a smoothed profile method resolving all scales, including the flow in the lubrication gap between the particle and the boundary. The centrifugation of heavy particles in the recirculating flow is counteracted by a repulsion from the shear-stress surface. The resulting limit cycle for the particle motion represents an attractor for particles in dilute suspensions. The limit cycles obtained by fully resolved simulations as a function of the particle size and density are compared with those obtained by one-way coupling using the Maxey–Riley equation and an inelastic collision model for the particle–boundary interaction, solely characterized by an interaction-length parameter. It is shown that the one-way coupling approach can faithfully approximate the true limit cycle if the interaction length is selected depending on the particle size and its relative density.


Introduction
Dispersed multiphase flows arise in a wide range of natural phenomena [1]. The problem is strongly characterized by the interaction between two fluid phases, a fluid and a particulate phase, or by both. Particle-laden flows are dealing with a continuously connected fluid phase and an immiscible dispersed phase made of particles. Examples for such flows are sand storms, debris flows, and transport of volcanic ash [2]. Particle-laden flows are also important in many technical processes like combustion [3], steel making [4], drug delivery and other biological applications [5]. Owing to the abundance of particle-laden flows their understanding, prediction and control has become an active research field for theoretical, experimental and computational fluid dynamics.
In experimental fluid mechanics, small particles suspended in the fluid are frequently used to visualize the flow [6]. The underlying assumption is a nearly perfect advection of the suspended particles. This condition is typically satisfied if the Stokes number is very small, the volume fraction of the particles is small and if the particles moves sufficiently far from the boundaries. While the application of this technique is not very complicated, a numerical simulation of the problem with even a single particle can be very demanding if all relevant length scales must be fully resolved. These scales range from the large scales of the global flow, over the medium scales for the flow around the particle, down to the scales of the flow in the lubrication film between the particle and the boundary. For a numerical solution of such problems, the use of efficient and accurate numerical methods is of crucial importance, in particular, for the simulation of the particle motion.
For very small particle volume fraction and low Stokes numbers, a numerical treatment by one-way coupling [3] seems to be justified. Within this approach, the feedback of the particles on the flow and particle-particle interactions are neglected. The uncoupling of the fluid flow from the particle motion makes this method extremely efficient, and a large number of problems have been treated successfully. However, a conventional one-way coupling approach may fatally fail to predict the behavior of particles if they experience repeated interactions with boundaries, even if the volume fraction and the Stokes number are small. Particle-boundary interactions become increasingly important in closed flows when the streamlines are compressed toward a smooth boundary. This situation typically arises in lid-driven cavities or surface-tension-driven flows. If the thickness of the layer of high volume flux adjacent to the boundary becomes of the order of the particle size the frequency of particle-boundary interactions for an individual particle can become very high in such recirculating closed flows. Since the domain of motion of the particle centroid in a typical one-way coupling approach is not restricted by the particle size, a dedicated theoretical model of particle-particle and particle-boundary interactions can make the difference between a realistic and a poor numerical prediction [7]. For this reason, the present study is devoted to the particle-boundary interaction via a fully resolved numerical simulation.
The finite particle size has been taken into account for the particle-boundary interaction in particular settings. These concern the settling of particles in a quiescent fluid in a gravity field (see, e.g., [8,9]) and particles moving near a vertical wall [10,11]. To the best of the authors' knowledge, the interaction between a particle and a boundary on which a shear stress is imposed in a recirculating flow has not been treated. Since particleparticle collisions are much more frequent than particle-boundary interactions in dispersed particle-laden flows, more attention has been paid to fully resolve particle-particle collisions. In corresponding simulations, the largest scales taken into account have been the typical particle scales [12][13][14]. Most of these studies served the purpose to derive models for particle-particle interactions mimicking colloidal forces in simulations in which the lubrication scales between the particles are not resolved for computational cost reasons (see, e.g., [15]).
In the following section, the problem of particle motion in a shear-driven cavity will be formulated mathematically. In Sect. 3 the numerical methods of investigation are introduced. In Sect. 4 the results of the numerical calculations are presented, focusing on the role of the particle-boundary interaction. Based on the comparison of the fully resolved simulations with results obtained by one-way coupling supplemented by the particle-surface interaction (PSI) model of Hofmann and Kuhlmann [16] (see also [7,17]), an improved particle-surface interaction model is suggested which closely approximates the particle-motion attractors obtained by the fully resolved simulations. The influence of the particle radius and the particle-to-fluid density ratio on the particle-boundary interaction is investigated for a constant shear Reynolds number Re = 1000. Finally, Sect. 5 is dedicated to summary and conclusions.

Problem formulation
We consider a shear-stress-driven square cavity filled with an incompressible Newtonian liquid with density ρ f and kinematic viscosity ν. The length of the cavity sidewalls is L. A plane flow is driven by a constant shear stress imposed along one side of the cavity. Figure 1 shows a sketch of the setup.
For the mathematical formulation of the problem, a viscous scaling is employed where x, t, u and p are the position vector, time, velocity field and pressure field, respectively. The caret indicates dimensional quantities. The non-dimensional Navier-Stokes and continuity equations read The flow field must satisfy the no-slip and constant-stress boundary conditions Re y x 1 a Fig. 1 Sketch of the two-dimensional shear-stress-driven square cavity. The particle sizes considered are shown to scale where the Reynolds number Re = (L 2 /ρ f ν 2 )τ 0 is a non-dimensional measure of the imposed shear stress τ 0 . The constant-shear-stress boundary condition is approximately realized in an open differential-heated low-Prandtl-number cavity due to the thermocapillary effect when the pressure due to the mean surface tension dominates all other normal stresses [18]. For a thermocapillary-driven flow the local shear stress on y = 1/2 is If the cavity is seeded with rigid particles, the above equations must be solved in the domain occupied by the fluid together with no-slip conditions on the surfaces of the particles. The motion of the particles (numbered by i) with mass M i and moment of inertia I i is governed by where F i and T i are the force and torque, respectively, acting on the i-th particle and V i and Ω i denote its translational and angular velocities, respectively. The coupling between the two phases results from the no-slip condition on the surfaces of the particles. Buoyancy forces acting on the particulate phase are neglected. This assumption holds true under weightlessness conditions or when the settling time of the particle is large compared to the characteristic time of the flow. The latter condition is met when the particles are sufficiently small as, e.g., in the experiments of Schwabe et al. [19].

Methods of investigation
The motion of a single finite-size particle is computed numerically, resolving all length and time scales ranging from the scales of the shear-driven flow over the scales of the particle to the scales of the lubrication gap between the particle and the boundaries. In addition, the flow is calculated in the absence of a particle and the particle motion is calculated by one-way coupling supplemented by a suitable particle-surface interaction (PSI) model when the particle moves close to any of the boundaries.

Fully resolved simulations
To solve the full problem, we use an Eulerian approach in which the flow-field equations are solved by a discontinuous Galerkin-finite-element method (DG-FEM) coupled to a smoothed profile method (SPM). In SPM particles are represented via a smooth concentration function φ, called smoothed profile. The smoothed profile method was originally proposed by Nakayama and Yamamoto [20] for simulating the interaction between fluid and particles in colloidal dispersions. The method has been devised to be used in combination with highorder methods and has successfully been employed for the simulation of particle-laden flows [13,14,[21][22][23], electro-charged colloids and electrophoresis of dense dispersions [24][25][26][27][28][29][30].
The space discretization is based on the discontinuous Galerkin-finite-element method. For all cases treated in the following, a nodal approach is adopted. Following Hesthaven and Warburton [31], triangular unstructured grids and a warp-blend distribution of nodes are employed. It has been shown by Luo et al. [32] that the convergence properties of this spectral/hp method are transferred to the SPM. This beneficial property provides a high flexibility via an element-wise mesh refinement (h-convergence) or a polynomial-order refinement providing exponential convergence (p-convergence). The details of the DG-FEM-SPM code are described in Romanò and Kuhlmann [33] including a comprehensive verification study.
Although we shall be concerned with a single particle only, the formulation will be given for a finite number N p of rigid particles. In the smoothed profile method, the particle and fluid subdomains are described by an indicator function φ. It takes the value 1 in the particle domain, 0 in the fluid domain, and it varies smoothly between these values in a transition zone between both phases. Following Nakayama and Yamamoto [20], the boundary between a particle and the fluid is smoothed by employing the smoothed profile for the single-particle indicator function where the index i enumerates all particles. The quantity d i is the signed minimum distance of x from the centroid P i of the i-th particle. It is defined such that it is positive outside and negative inside the particle. The thickness of the transition layer for the i-th particle is denoted ξ i . The profile (5) has been shown to yield very good results for a wide range of flows [33]. For non-overlapping particles, the global indicator field φ is defined as the linear combination of all single-particle indicator functions φ i . Thus, for N p rigid bodies, The rigid-body motion of a set of non-overlapping particles is then represented by a global particle velocity field u p which satisfies [32] φ( The total flow field u defined in both the fluid and the solid domain is written as a combination of the particle (u p ) and the fluid-phase velocity fields (u f ) using the indicator function as a weight It can be shown that the particle velocity field (7) is exactly solenoidal [20] and no-slip and no-penetration conditions are automatically satisfied by a smooth blending of both flow fields [20,32,33]. For the time integration of the global flow field, we employ the second-order stiffly stable splitting scheme of Karniadakis et al. [34]. The time discretization of (2) and (4) then leads to the following algorithm where represents the nonlinear term written in conservative form. The vectors r i and O i are the instantaneous position and rotation rate vectors, respectively, of the i-th particle. The perturbation-flow field caused by the particles is denoted u p , and the superscript * refers to the flow field obtained without updating the positions and velocities of the particles. The coefficients a k are derived from the backward-differentiation formula used for integrating the particle positions and velocities. The total pressure field can be reconstructed by p n+1 = p * + p p and the coefficients γ 0 , α 0 , α 1 , β 0 , β 1 are tabulated in Karniadakis et al. [34]. For small solid-to-liquid volume ratios φ V < 10 −3 particle-particle interactions can be neglected. For a fully resolved approach it is, therefore, sufficient to investigate the motion of a single circular particle only. In the following, the method is applied to the two-dimensional motion of a single circular particle in the shear-driven square cavity, because the two-dimensional case is the logical first step and greatly facilitates a variation of the particle parameters which could be extremely costly in three dimensions when the lubrication gap becomes very small. The fully resolving code has been extensively verified in two dimensions by Romanò and Kuhlmann [33] who have reproduced numerous two-dimensional benchmarks with high accuracy.

One-way-coupled simulations
In addition to the fully resolved simulations of circular particles in the two-dimensional shear-driven cavity, one-way-coupled simulations of spherical particles in a steady two-dimensional flow are carried out. An approximation to the steady flow is obtained by solving (9)c to (9)f until changes of the flow field have become negligible. A steady state is assumed to be reached once the convergence criterion was satisfied, where j is the number of grid points.

One-way-coupled particle motion
The motion of particles of sufficiently small size and volume fraction is frequently modeled by the Maxey-Riley equation [35]. Neglecting the Basset term, the Saffman and Faxén corrections and gravitational forces, the Maxey-Riley equation readsÿ where we have used the same scaling (1) as before, y is the position vector of the particle and D/Dt is the substantial derivative. In (11) we neglect finite-Reynolds-number corrections. This approximation is motivated in the "Appendix". The non-dimensional groups are the particle-to-fluid density ratio and the Stokes number defined as respectively, where a is the radius of a spherical particle. For a given steady flow field (11) is solved by a standard fourth-order Runge-Kutta method and time step Δt = 10 −6 . To verify the Maxey-Riley solver according to (11), we consider the motion in the Taylor-Green vortex flow [36] Figure 2 shows a comparison of the present results with those of Domesi and Kuhlmann [37] for a spherical particle whose centroid velocity is initially mismatched to the unperturbed fluid velocity. The agreement is very good in view of the sensitive dependence of the trajectory from the initial conditions due to the abundance of hyperbolic points in the flow. The simulation parameters are given in the caption of Fig. 2. The integration time was selected to demonstrate the high accuracy of the solver for a time window which is ten times larger than the typical integration time employed for integrating a particle in the investigated shear-stress-driven cavity flow. The Maxey-Riley equation is valid for small Stokes numbers St 1 and small particle Reynolds numbers Re p = a|ẏ − u|/ν 1. Moreover, the particle is assumed to move far from any boundary. Therefore, it is necessary to supplement (11) with a model which takes care of the particle-boundary interaction. This is particularly important in shear-driven flows, e.g., thermocapillary flows, because the highest velocity arises near the stress interface causing a very high volume flux near the liquid-gas interface. Therefore, suspended particles will be frequently transported very close to the interface. This necessitates a dedicated modeling of the boundary effect on the particle motion which will be presented in Section 4.2.1.

Results
We consider the motion of a circular particle near the stress boundary in a steady two-dimensional shear-stressdriven cavity flow introduced in Sect. 2. In such a closed vortex flow counteracting forces on the particle, due to inertia on the one hand side and due to the boundary effect on the other hand, can lead to a limit cycle for the particle motion. The parameter dependence of the limit cycle is investigated, and a model is proposed to take into account the boundary effect within a one-way coupling approach. For all fully resolved simulations, we employ ξ = 0.05a and select Δt according to Luo et al. [32]. This choice of parameters is intended to minimize the numerical error [33]. We employ fifth-order polynomial basis functions and a structure of elements which guarantees that all relevant flow scales are resolved. To meet with the latter requirement, two different grids are used. When the particle is moving near the stress boundary, the relevant scales to be resolved are the particle and the lubrication-gap length scales. Therefore, the computational mesh is designed to always have at least five elements (i.e., ∼100 grid points) between the surface of the particle and the shear-stress boundary and at least ten elements (i.e., ∼200 grid points) inside the rigid-body domain. When the particle moves further away from the shear-stress surface in the bulk of the liquid, a second computational grid is employed which only needs an accurate resolution of the particle scale. Passing from one grid to the other is accomplished by spectral interpolation. A representative example of the two computational grids is depicted in Fig. 3. The structure is motivated by the streamlines of the flow shown in Fig. 4.

Poincaré maps
The motion of a circular particle in a shear-driven cavity of Re = 1000 is considered. The main closed streamlines of the steady two-dimensional flow in the absence of the particle are shown in Fig. 4. The streamlines are bounded by a separating streamline along the wall which also separates the small and weak viscous corner eddies near (x, y) = (±1/2, −1/2) (not shown) from the main recirculation. The streamlines are crowded toward the shear-stress boundary at y = 1/2 where the highest velocities arise. Due to the no-slip boundary conditions on the remaining walls and the associated reduction in the mass flow rate, the streamlines are much wider near the solid walls than near the stress boundary. Therefore, the shear-stress boundary is much more important for the particle motion than the no-slip boundaries.
To distinguish between the effects on the particle motion caused by the shear-stress boundary and those caused by the solid walls including the bulk motion we consider Poincaré maps on y P = 0.3 (dashed line in Fig. 4). This plane subdivides the full domain into a region y P ≤ y < 1/2 near the shear-stress surface and its complement (bulk motion). The streamlines in y ≥ y P have a small curvature near x ≈ 0 and become more and more rectilinear the closer they are to the shear-stress surface. The streamlines in y < y P , on the other hand, are approximately circular. The small particles considered below have trajectories which do not deviate very much from the streamlines. The deviations from the nearly straight streamlines in y > y P are expected to be primarily caused by forces on the particle due to the proximity of the shear-stress surface, whereas for y < y P centrifugal forces are the primary cause for deviations from the streamlines. The choice of y P = 0.3 thus allows to approximately decouple the surface effects from the inertial effects and also guarantees a complete coverage of all trajectories of interest.
To analyze the particle motion the Poincaré section on y P is calculated for two Poincaré points when the particle is initiated at y sb s < y P , where the subscript s indicates the start position and the superscript sb refers to the motion in y > y P near the stress boundary. The Poincaré map P sb : x sb in → x sb out = P sb (x sb in ) defines a surface map which is primarily due to the particle-surface interaction (Fig. 5a). The bulk map with y b s > y P defines a bulk map (Fig. 5b). Due to the small Stokes numbers considered, the velocity of the particle relaxes from its initial value much faster than any other time scale involved [38]. Therefore, the initial velocity does not practically affect the mappings, in particular, if the particle is initialized velocity-matched and sufficiently far from y P . It turns out that initializing the particle's centroid velocity-matched at y sb s = 0.2 < y P is sufficient for the surface map. For the bulk map we select y b s = 0.4 > y P . The overturning motion for a full revolution of the particle is described by the oriented Poincaré map P : It is equivalent to the combined maps consisting of the bulk and the stress-boundary maps P sb (x) and P b (x), respectively (Fig. 5c).

Existence of a globally stable limit cycle
To reduce the parameter space, Poincaré maps are investigated for a heavy particle with particle-to-fluid density ratio = 2 and for Re = 1000. Three particle radii, a = 0.01, 0.03 and 0.05, are considered.
To establish the Poincaré maps, fully resolved simulations are performed for a set of initial conditions imposed at (x s , y s ) such that the first Poincaré point on y P is approximately covering the range x ∈ [0.2, 0.45]. Figure 6 shows the stress-boundary map P sb (open symbols) and the inverse bulk map [P b ] −1 (full symbols). The particle size is indicated by the type of symbol. A typical mapping along y P = 0.3 is illustrated by arrows for a = 0.05 (triangles). From x 0 = 0.255 (open diamond) the particle is mapped via P sb (x 0 ) = −0.39 to Fig. 5 Sketches illustrating the stress boundary, bulk and iterative maps P sb , P b and P, respectively, as well as the notation used. Poincaré points are indicated by small dots If the particle were only advected with the flow it would return to its initial position, i.e., x 1 = x 0 . The corresponding maps P sb ψ and [P b ψ ] −1 are indicated by a dotted curve in Fig. 6. Since the slope of these advection maps is less than one, |∂ x P sb ψ | = |∂ x [P b ψ ] −1 | < 1, the Poincaré points and thus fluid elements become compressed in x during the motion along the stress boundary, and dilated by the same amount during the bulk motion. This trivial effect is due to the asymmetry of the streamlines (Fig. 4).
We find the absolute value of the slope of the inverse bulk map for the particle motion always to be larger than that of the inverse map for the streamlines, |∂ ψ ] −1 |. Therefore, the dilatation of particle trajectories is less during motion in the bulk than that of the streamlines. This means that the particle trajectories are compressed in x relative to the streamlines. This behavior can be explained by the larger centrifugal displacement a particle experiences in the bulk (y < y P ) when moving closer to the vortex core as compared to a particle moving further away from the vortex center. For the motion along the stress boundary, particle trajectories become even more compressed in x than the streamlines, because |∂ x P sb | < |∂ x P sb ψ |. Obviously, this effect is caused by the increasing repulsion force a particle experiences the closer it moves to the shear-stress surface. Therefore, the distance between any two particle trajectories, measured by the distance |x n+1 − x n | of the corresponding Poincaré points, is monotonically shrinking with n, i.e., in time.
Interestingly, the graphs of P sb and [P b ] −1 intersect (Fig. 6) yielding a fixed point for the combined map P for a full revolution of the particle shown in Fig. 7. Therefore, the combined map P has a single fixed point x * for each particle size (Fig. 7). Since the combined Poincaré map P(x) is monotonic with slope 0 < ∂ x P(x) < 1 for all particle sizes, the corresponding limit cycle of the particle motion is stable and globally attracting. From the slopes P (x * ) (Fig. 7) the asymptotic rate of attraction is larger the larger the particle. Evidence for the existence of a stable, global limit cycle has been provided by [39] who found two-dimensional axisymmetric particle accumulation structures in cylindrical thermocapillary liquid bridges.
The existence of a stable limit cycle for heavy particles > 1 can also be heuristically understood by the competition of two effects. As long as the particle is transported far away from the stress boundary, the inertia effect dominates and the particle is centrifuged out of the center of the vortex. As the particle is transported to the outer streamlines of the vortex, it experiences a repulsive force from the stress boundary. The repulsion dominates for particles moving very close to the wall and the stress boundary. For a specific intermediate trajectory centrifugal forces and surface repulsion exactly balance in the mean along the limit cycle.

Parameter dependence of the limit cycle
The limit cycles for Re = 1000 and = 2 are shown in Fig. 8 for a = 0.01, 0.03, and 0.05. Close to the stress boundary the particle's trajectory is nearly parallel to the boundary, and the liquid film between the particle and the stress boundary is much smaller than the radius of the particle. A characteristic quantity of the limit cycle is the minimum distance of the particle's centroid from the stress boundary. We thus introduce the interaction   length Δ = a + δ where δ is the minimum film thickness between the surface of the circular particle and the shear-stress surface. Figure 9 shows the normalized interaction length Δ/a as a function of the particle radius (open squares). The normalized lubrication gap width δ/a = (Δ/a) − 1 exhibits a strong monotonic dependence on the particle radius a: the larger the particle the more it is capable, in proportion, of squeezing the lubrication gap (a 1 > a 2 ⇒ δ 1 /a 1 < δ 2 /a 2 ). In all cases investigated the lubrication gap is less than 25% of the particle radius. On the other hand, the absolute lubrication gap width δ (full squares in Fig. 9) takes a local minimum. Since centrifugal inertia increases with the particle radius like ∼a 2 (for the two-dimensional fully resolved simulations), the lubrication gap δ should decrease with a. This effect seems to hold for sufficiently small a. However, the equilibrium orbit also depends on a. As the particle radius increases, the local curvature of the closed orbit just before approaching the stress boundary and the particle velocity decrease. Both effects reduce the centrifugal inertia and could thus contribute to the increase of δ(a) beyond the minimum.
To establish the dependence of the limit cycle and the associated interaction length on the particle-to-fluid density ratio nine different density ratios ∈ [1.2, 2] have been considered for Re = 1000 and a = 0.05. As before, the limit cycles are obtained as fixed points of the Poincaré maps P shown in Fig. 10b obtained combining the stress-boundary and the bulk maps presented in Fig. 10a. The locus x * of the limit cycle depends sensitively on the density ratio. Since the slopes of the combined maps P (x * ) at the fixed points do not vary much, the asymptotic attraction rate to the limit cycle is nearly independent of the density ratio. Note, however, that the curvature P (x * ) > 0. Therefore, the approach of the limit cycle due to inertia from the bulk (for x < x * ) is slower than the attraction due to boundary repulsion from the surface (for x > x * ). This shows that the boundary effect acts on a faster time scale, for the present parameters, than the inertia effect. The stable periodic orbits are shown in Fig. 11. Due to the inertia of heavy particles with > 1 (centrifugal forces) the periodic orbits are larger the larger . We find the limit cycles for different to be continuously nested. Moreover, as the density becomes smaller the minimum distance of the limit cycle from the stress boundary is significantly larger than the particle radius a. The normalized interaction length Δ/a is shown in Fig. 12 as function of . As expected from Fig. 11, the lubrication gap increases significantly as the particle density is approaching the fluid density from above ( ↓ 1). To fit the data, we consider the limit of large and small relative density. For → ∞, the inertial forces will squeeze the lubrication gap completely. Therefore, Δ/a should admit an horizontal asymptote Δ/a = 1. Moreover, if = 1, the boundary repulsion will not be opposed by inertia. In this case, the particle is expected to move far away from the boundaries on a trajectory in the vicinity of the vortex core which is determined by some equilibrium between the repulsive actions due to the four boundaries of the cavity. In the limit a → 0 such a condition becomes a vertical asymptote for Δ/a. These considerations suggest the fit where A and s are constant fit parameters and B(a) is a function of the particle radius with lim a→0 B = 1. For a = 0.05, we find the least squares fit (full line in Fig. 12) The functional dependence of the interaction length Δ = f (a, ) on the particle radius a, relative particle density and, possibly, on the flow parameters can be exploited to construct a simplified model for the forces acting on the particle when its motion is approximated by one-way coupling.

Particle motion near a shear-stress boundary: one-way coupling
In view of the numerical effort required to simulate the flow and the motion of finite-size particles fully resolved on all length scales, a one-way coupling method is highly desirable for practical reasons in which the boundary effect on the particle motion can be approximated by a simple model. For no-slip boundary conditions along rigid walls such a model has been suggested by Brenner [40]. For dilute suspensions in a flow driven by a constant shear stress on a boundary, however, the particle-boundary interaction is much more important than the particle-wall interaction. Since we are not aware of any model valid near a shear-stress surface derived from first principles, we build on the inelastic collision model introduced by Hofmann and Kuhlmann [16].

Modeling the particle-boundary interaction
In order to model the strong repulsion forces a finite-size particle experiences during its motion near a free surface subject to thermocapillary shear stresses, Hofmann and Kuhlmann [16] have proposed the particlesurface interaction (PSI) model for one-way-coupled simulations. In their model, the particle motion is assumed to be governed by the Maxey-Riley equation up to the point at which the surface of the particle makes contact δ a a Δ release impact Fig. 13 Illustration of the particle-surface interaction (PSI) model. The solid line symbolizes the trajectory of a circular particle (its size is indicated by dashed circles) obtained using the Maxey-Riley equation disregarding the boundary effect. The dotted line represents the trajectory of the centroid of the same particle when the PSI model is applied. The interaction length Δ, comprising of particle radius a and the minimum lubrication gap width δ, defines the size of a prohibited region in which the particle's centroid cannot enter (dark gray area). Dots particle centroids, circles impact and release points, dashed line stress boundary with the free surface, assumed to be non-deformable. The lubrication gap was thus assumed to be vanishingly small. After an inelastic contact the particle is assumed to slide along the free surface, being pushed toward it by the outward-normal velocity component of the flow at the particle centroid, up to the release point at which the surface-outward-normal velocity component of the flow turns negative upon which the particle detaches from the free surface. Introducing the interaction length Δ, Mukin and Kuhlmann [17] pointed out that the minimum distance of the centroid of the particle from the free surface is necessarily larger than the particle radius and the lubrication gap must be included in Δ. Such PSI model is also justified by the observation that the true particle trajectory is nearly rectilinear and parallel to the stress boundary over quite some distance for small and heavy particles. The model can be understood as a contact force acting at a distance Δ from the boundary which instantaneously annihilates any component of the velocity which is directed outward normal to the boundary. Figure 13 depicts a sketch of the PSI model in a two-dimensional flow. At the impact point, which has a distance Δ from the stress boundary, the trajectory of the particle centroid (dotted line) is assumed to deviate from the trajectory predicted by the Maxey-Riley equation (solid line). Within the model, the particle centroid cannot enter the layer of thickness Δ on the boundary (shaded in dark gray). At the release point at which the surface-outward-normal velocity turns negative on y = 1/2 − Δ the particle centroid detaches from the prohibited region.
The results of the fully resolved coupled simulations presented in Sect. 4.1.3 suggest a refinement of the PSI model by including the lubrication gap width in the interaction parameter Δ. As an approximation to the true lubrication gap, we shall use the lubrication gap obtained for the limit cycles during the fully resolved simulations of the particle-boundary interaction process. This approximation insures that the minimum distance from the stress boundary of a periodic orbit obtained by the fully resolved simulation will be identical to that of the periodic orbit obtained from one-way coupling. However, the dynamics of the approach to the periodic orbit will still differ between both methods. Figure 14 shows a comparison of limit cycles from fully resolved simulations (full line with markers) with those obtained by solving the Maxey-Riley equation (11) together with the PSI model (full line without markers). In the PSI model the minimum distance of the particle centroid from the stress boundary is enforced to be the interaction length Δ = a + δ. For = 2, Re = 1000 and a = 0.05, we find a good agreement of the corresponding limit cycles. For a = 0.03, the agreement is still reasonable and for the smallest particle size the agreement between the limit cycles is good as well. We conclude that, for the range of parameters investigated, limit cycles obtained by fully resolved simulations are good approximations to the true limit cycles if the PSI model is used taking into account the lubrication gap. In addition, the trajectory of the particle centroid is shown as a dashed line for the case when the Maxey-Riley equation (11) is solved without the PSI model. The particle was initialized velocity-matched to the flow at (x, y) = (x * , y P ), where x * is the fixed point of the combined Poincaré map (Fig. 7). As can be seen, a limit cycle is absent when ignoring the particle-boundary interaction and the trajectory is asymptotically approaching the separating streamline, a result which is completely wrong.

Comparison of fully resolved simulations with one-way coupling approximations
A further comparison is made between the limit cycle of the fully resolved simulation and that of the Maxey-Riley equation combined with the PSI model for a = 0.05 and = 2. Figure 15 depicts the coordinate of the particle's centroid normal to the stress boundary (upper row), its normal velocity (middle row), and its normal acceleration (lower row) as functions of time during the approach, interaction and release from the stress boundary. The origin of time is the instant at which the particle centroid passes the line y p for x > 0. From the lower row, it can be seen that the deceleration of the particle in the fully resolved simulation (full line with triangles) reached a significant maximum just before t ≈ 0.01 during the approach of the stress boundary. In contrast, the particle moving according to the Maxey-Riley equation (full line without triangles) experiences a much smaller deceleration. Within the PSI model, this lack of deceleration is compensated by an impulsive acceleration at t ≈ 0.01 slightly after the maximum deceleration of the particle in the fully resolved simulation. Obviously, the impulsive acceleration and the PSI model lead to instantaneously transferring the particle practically to the same distance from the stress boundary as the particle in the fully resolved simulations attains somewhat later (upper row). As a result of the PSI model the velocity normal to the stress boundary returns to zero instantaneously (middle row), whereas the normal velocity of the particle smoothly tends to zero in the fully resolved simulation. In both cases, the normal velocity is exactly (or close to) zero for a considerable period of time. It is seen that the departure (release) of the particle from the boundary is slightly premature as compared to the fully resolved simulation. This hints at a possible improvement of the release condition within the PSI model which currently is defined by the zero of the normal acceleration of the fluid at the location of the particle (dotted curve in the bottom row). Comparing the dynamics obtained from the Maxey-Riley simulation without the PSI model [dashed lines, initiated velocity matched at (x * , y p )] leads to the following deficits as compared to both the fully resolved simulation and results of the Maxey-Riley equation with PSI modeling: (a) without PSI model the particle approaches the stress boundary much too closely (upper row). (b) the deceleration during the approach of the stress boundary is much too small compared to the fully resolved simulation (lower row). (c) The release from the stress boundary is too premature and occurs even earlier than in case of the PSI model. (d) The Maxey-Riley equation without PSI modeling does not exhibit a non-trivial limit cycle in a certain distance from the stress boundary. This shows that the PSI model in combination with the Maxey-Riley equation yields a much better approximation to the result of the fully resolved simulation than the Maxey-Riley equations without PSI modeling which completely fails for t → ∞.
The influence of the relative particle density on the limit cycles obtained using different approximations is shown in Fig. 16. Only the near-surface region is displayed. Shown are results obtained by fully resolved simulations (full lines with symbols), by one-way coupling (11) using the PSI model with Δ = a (dashed lines), and by one-way coupling employing the PSI model with Δ = a + δ including the lubrication gap (full lines without symbols). As before, Δ = a + δ yields a good approximation of the limit cycle for all density ratios. Disregarding the lubrication gap width in one-way coupling is a reasonable approximation for a large particle-to-fluid density ratio. However, as the density ratio decreases, the uncorrected interaction parameter Δ = a is no longer a good approximation to the distance of the limit cycle from the stress boundary which results in an appreciable error in the shape and location of the limit cycle. The inclusion of the lubrication gap  into the interaction parameter Δ is particularly important when the particle density becomes comparable to the fluid density, i.e., for ≤ 1.2 (see Fig. 16e).

Discussion and conclusions
The motion of finite-size particles in the closed recirculating flow inside a square cavity has been investigated. The flow was driven by a constant shear stress with Re = 1000 which was imposed along one of the four boundaries. The simulations were carried out for heavy circular particles of different radius and different relative densities. The flow has been resolved on all length scales including scales smaller than the thickness of the lubricating film between the particle and the stress boundary.
A first key result is the existence of a global limit cycle for the motion of heavy particles of finite size. It is created by opposing centrifugal and boundary-repulsion forces. The attraction is due to the dominance of centrifugal forces in the mean if the particle centroid moves in the flow region enclosed by the limit cycle (near the vortex center). Particles with centroids outside of the limit cycle are attracted to it by the dominance in the mean of the repulsive forces from the boundaries.
The limit cycle is characterized by its minimum distance Δ from the stress boundary, the repulsion from it being the dominant boundary effect. The interaction length Δ( , a) was found to depend significantly on the particle radius a and on its relative density . The functional dependence of Δ obtained from the fully resolved simulations was used to correct the particle-surface interaction (PSI) model of Hofmann and Kuhlmann [16] by including the lubrication gap of width δ between the particle and the stress boundary in the interaction parameter Δ = a + δ. Limit cycles in the square cavity by one-way coupling using Δ = a + δ were found to be in good agreement with the limit cycles by fully resolved simulations, whereas limit cycles obtained using one-way coupling with the uncorrected interaction length Δ = a deviated significantly from the reference Limit cycles for Re = 1000 and a = 0.05 close to the shear-stress surface for different density ratios = 2, 1.8, 1.6, 1.4, and 1.2 (from a to e). Lines with symbols: fully resolved simulations; full lines: one-way coupling using the PSI model with Δ = a + δ according to Fig. 12; dashed lines: one-way coupling using the PSI model with Δ = a. The dark and bright gray bars indicate the minimum distance of the respective one-way-coupled limit cycle from the stress boundary limit cycle. The smaller the particle radius and the smaller the density difference between particle and fluid with > 1 the more important, relative to the particle radius, becomes the lubrication effect on Δ and, thus, on the location and shape of the limit cycle. A second key result is the demonstration that one-way coupling combined with the PSI model can predict periodic orbits in two dimensions with good accuracy. It is concluded that the particle/stress-boundary interaction can be effectively modeled, for the parameters investigated, by one-way coupling using the PSI model if the lubrication-gap width is properly taken into account. A one-way coupling approach without any particle-surface interaction model is not able to reproduce the non-trivial limit cycle and will thus miss this essential feature of particle accumulation.
The PSI model for one-way coupling with Δ = a + δ can well describe the location and shape of the limit cycle. However, it does not correctly reproduce the dynamics of the attraction to the limit cycle. This applies particularly to particles which are initially located on the outer side of the limit cycle. Within the PSI model these particles are immediately transferred to the limit cycle upon a first interaction with the stress boundary, whereas the fully resolved simulations show the particles to only asymptotically approach the limit cycle. The attraction of particles from the interior of the limit cycles is expected to be very similar in the framework of both the PSI model and the full simulation, because their motion is primarily determined by inertia as these particles are centrifuged toward the limit cycle. Several extensions would be of interest. In order to broaden the database for the interaction parameter further, fully resolved simulations are required to study the effect of other governing parameters like the Reynolds number. Moreover, not only the particle/stress-boundary interaction has to be modeled in a one-way coupling approach, but also the lift forces which a shear-flow over a wall exerts on the particle when the particle moves close to it. Finally, particle-particle interaction can be of crucial importance to understand the dynamics and stability of particles densely occupying the limit cycle, after they have been attracted to it essentially by a single-particle mechanism.