The interaction between a solid particle and a turbulent flow

The interaction between a fixed solid spherical particle and stationary turbulence with zero mean flow is investigated numerically. The object diameter, D, lies in the inertial range (D≈0.6L≈0.9λ≈8η, where L, λ and η, respectively, denote the integral scale, the Taylor microscale and the Kolmogorov length) and the particle Reynolds number is close to 20. It is found that the turbulence statistics at different distances from the solid/fluid interface are modified by the presence of the object in a region that extends more than 10 times further than the viscous layer. This estimate is confirmed by the analysis of the correlation between the force and torque on the particle and the force and torque on spherical surfaces surrounding the particle, although the torque decorrelates somewhat faster with increasing distance from the object surface. The angular slip velocity of the particle, a quantity of crucial importance for the modeling of the turbulent transport of large objects, is also characterized.


Introduction
Solid-fluid disperse multiphase systems play an important role in a variety of natural and industrial processes. Sediment transport in rivers and at the ocean bottom, clouds, dust storms and pollen dispersion are examples of natural systems of this type, while fluidized bed reactors and combustors, suspensions and oil recovery can be cited as instances in industry and technology. In most cases of interest, the flow is turbulent, which increases severalfold the complexity of the situation. Their widespread occurrence and economic importance have made these systems the object of an intense experimental, analytical and numerical study.
The dynamics of particles suspended in a given turbulent flow depends on their size and density: the two relevant parameters are the ratio between the object's diameter and the Kolmogorov scale of the flow, D/η, and the ratio between the particle and fluid densities. Very small, neutrally buoyant particles behave like passive tracers: this property is frequently used experimentally to study one-phase flows (laser Doppler velocimetry (LDV), particle image velocimetry (PIV) and particle tracking velocimetry (PTV) techniques). Otherwise, inertia effects come into play and the object's dynamics are different from those of fluid particles. For instance, numerical and experimental investigations of the collective behavior of inertial pointparticles in turbulence have shown that their distribution is far from uniform: light particles are concentrated in vortex cores (see e.g. [1]), while heavier ones are expelled from these structures under the action of the centrifugal force (see e.g. [2]). The distribution of inertial particles in the flow is therefore highly nontrivial. In most cases, the fluid/particles density mismatch leads to a transport problem with some analogies to compressible flow [3].
The dynamics of particles whose spatial extent can be neglected (in practice, such that D < η) is now fairly well understood, especially in conditions such that hydrodynamic drag is the dominant force. In this situation, one can describe the fluid phase in a Eulerian fashion by means of the Navier-Stokes equations, whereas each particle is tracked individually in a Lagrangian fashion by integrating its dynamical equation in which the hydrodynamic force is expressed in a parameterized form (see e.g. [4]- [6]). A crucial quantity appearing in this parameterization is the particle-fluid 'slip velocity', which is written as the difference between the particle velocity and the fluid velocity in its neighborhood, assumed substantially uniform 3 over the particle scale. In this way, it is possible to simulate the dynamics of several millions of point particles, thus nearly approaching situations of practical interest.
The situation is considerably more complicated if one considers particles whose spatial extent cannot be neglected, namely such that D > η, as the applicability of the point-particle models becomes then questionable [7], although its extension to particle diameter of several η can be used [8]. An obvious difficulty stems from the fact that the very notion of slip velocity loses meaning if the flow cannot be approximated as uniform at the particle scale. Furthermore, in these cases, the back-reaction of the particle on the surrounding flow, which is included somewhat indirectly in the point-particle model, must be taken into account with greater fidelity.
Several experimental studies have been recently carried out to characterize the dynamics of particles larger than the Kolmogorov scale [9]- [14]. In particular, the acceleration statistics of the objects, a quantity that reflects the hydrodynamic forces acting on them, have been measured for different values of the parameters. The most surprising result is the fact that the probability distribution function of the acceleration components normalized by their variance does not seem to depend strongly on particle size or density. In particular, as D tends to the integral scale, these distributions remain strongly non-Gaussian. The rotation of spherical particles comparable to the integral scale has been found to be very intermittent [14].
In order to understand these phenomena, and thereby to improve the modeling of two-phase flows, it is essential to carry out fully resolved numerical simulations of flows with extended particles based on first principles. For this purpose, one has to integrate the Navier-Stokes equation accounting for the no-slip condition at the particle surfaces. A number of methods have been proposed for this purpose: finite-element approaches (e.g. [15,16]), the Chimera algorithm [17], pseudo-penalization (PP) [18], immersed boundary (IB) [19] and lattice-Boltzmann (LB) methods [20] and an interesting combination of the LB and IB approaches [21,22]. Other authors have focused on particular methods enabling one to simulate a turbulent flow around a single finite-size particle [23]- [25].
In this paper, we use the Physalis algorithm [26]- [29], specifically designed to integrate the Navier-Stokes equations in the presence of spherical solid objects. This method has been tested in various laminar flows [29]. It is based on first principles, and the zero-velocity condition at the fluid/solid interface is applied exactly, without interpolation or averaging. We use this method to simulate the interaction between a fixed solid spherical particle and stationary turbulence without mean flow. The object diameter lies in the inertial range and the particle Reynolds number, based on this lengthscale and on the velocity fluctuations, is ≈20. We will focus, in particular, on (i) the back-reaction of the object on the turbulence and (ii) the angular slip velocity.
The paper is organized as follows. In section 2, we describe the numerical method: the Physalis algorithm and its implementation in the case of fixed particles are recalled in sections 2.1 and 2.2, and the specificities for its application to turbulence are given in section 2.3. The results of our study are gathered in section 3. We first characterize the local modification of turbulence due to the presence of the object by investigating the flow statistics at different distances from it (section 3.1) and the hydrodynamic forces and torques acting on spherical surfaces of various radii (section 3.2). The angular slip velocity of the object is then characterized in section 3.3. Concluding remarks are finally given in section 4.

Numerical algorithm
The Physalis method is briefly described in the next two subsections. More details, including the method in the case of freely moving particles, can be found in [28,29]. The modifications of the algorithm for its application to turbulence are given in section 2.3.

The Physalis method in the case of fixed particles
Let us consider a fixed spherical particle surrounded by an incompressible viscous fluid animated by a velocity u. The fluid motion is described by the Navier-Stokes equations with u = 0 on the particle surface. Here p is the piezometric pressure, and ρ and µ are the density and dynamic viscosity of the fluid. The left-hand side of equation (1) equals zero on the particle surface. Therefore, by continuity, this quantity must be small in a thin region surrounding the particle and, as a consequence, the right-hand side of equation (1) is small near the particle as well. Therefore, locally, u and p approximately satisfy i.e. the Stokes equations. This result can be easily interpreted: due to the no-slip condition, the fluid in contact with the particle is exactly stationary. Equation (3) is simply a linearization about this state of rest. For any finite Reynolds number Re, there exists a region where (3) is a good approximation to (1), even if its thickness is expected to decrease as Re increases. The general solution of the Stokes equations in the presence of a spherical boundary was given in [30], in terms of harmonic potentials: 5 where R is the particle radius and ν is the kinematic viscosity of the fluid. The functions p n , φ n and χ n are regular harmonics representing the incident flow. For instance, p n = r R n n m=0 P nm cos(mϕ) +P nm sin(mϕ) P m n (cos θ), where (r, θ, ϕ) are spherical coordinates centered at the center of the particle, P m n is an associated Legendre function and P nm andP nm are dimensionless coefficients. The other regular harmonics φ n and χ n are written in a similar way in terms of dimensionless coefficients nm , nm and X nm ,X nm , respectively.
The functions p −n−1 , φ −n−1 and χ −n−1 are singular harmonics, representing the flow disturbance due to the presence of the particle. The form of these harmonics satisfying the zero-velocity condition on the particle surface is The condition of vanishing velocity on the particle surface is therefore imposed exactly on the solution, which can be finally written as an infinite series of known functions multiplied by unknown coefficients. The suitable determination of these coefficients enables us to represent any flow taking place in the immediate vicinity of the particle where equation (3) is applicable. At this stage it should be stressed that the series previously mentioned represents the most general linearized flow compatible with the boundary condition on the particle surface. No assumption about the particular form of the flow away from the particle has been introduced. The plan of the calculation consists in determining the coefficients of this local solution, valid in the neighborhood of each particle, by matching the fields given by equations (5) and (6) with those calculated numerically in the rest of the domain.
For this purpose, the computational domain is covered with a Cartesian fixed grid irrespective of the presence of the particles. In general, the problem arising from the use of regular grids is the lack of conformity of the particle geometry to the grid structure. To overcome this difficulty, we use the previous analytical form of the flow fields valid near the particle. Matching the analytical and numerical solutions on the grid nodes belonging to a cage surrounding the particle (see figure 1) has the effect of implicitly accounting for the correct boundary condition on the particle surface. As a consequence, there is no need to explicitly account for the presence of the particle by locally modifying the grid structure.
The hydrodynamic force and torque acting on the particle can be expressed directly in terms of the low-order coefficients appearing in these expansions: F = π µν(P 11 + 6 11 ;P 11 + 6˜ 11 ; P 10 + 6 10 ), T = 8π µν R(X 11 ;X 11 ; X 10 ). , pressure points; ×, vorticity points; , horizontal velocity points; , vertical velocity points. The same representation in a three-dimensional flow can be found in [29].
The computational method described above exhibits several advantageous features as follows.
• The no-slip condition at the particle surface is imposed exactly.
• As the number of degrees of freedom used to describe each particle increases, the error decreases exponentially, rather than algebraically.
• As equations (11) and (12) show, the force and torque acting on the particles are obtained directly: there is no interpolation of the stress to the particle surface, nor integration over it.
• The computation time mostly depends on the extent of the computational domain, and only weakly on how many particles are contained in it.

Implementation in the general case
We give here some details about the implementation of the method previously introduced.

Algorithm.
Briefly, at each time step, the calculation proceeds as follows.
1. An available numerical estimate of the pressure and vorticity fields (e.g. values at the previous time step) is used to determine the coefficients appearing in the analytical solution by matching the numerical and analytical solutions at the grid nodes surrounding each particle. 7 2. These coefficients are used to calculate analytically from (5) the velocity at the grid nodes adjacent to the particle. 3. The Navier-Stokes equations are numerically solved by using a standard projection method (see details in section 2.2.2) and these velocity values as internal boundary conditions. 4. The coefficients are updated by using this numerical solution; the procedure is repeated until convergence.
We give here some details of particular aspects of the method. More information can be found in [29].

Flow solver.
We use a standard staggered grid arrangement (see figure 1), with pressure at cell centers, velocity components at the midpoints of cell sides, and vorticity components at the midpoints of cell edges [29]. As explained before, the calculation progresses from time level t n to t n+1 = t n + t by matching iteratively the local analytical solutions valid near the particles with the finite-difference one. Let us denote by κ the counter of these iterations.
For each iteration, we first used the projection method described in [31] slightly modified for the present purpose [29]. It turned out that, in the present case, a method of first order in time (explicit method) was more efficient than the previous one (implicit method). We therefore use the following projection method (1st order in time and 2nd order in space).
A provisional estimate u * of the new velocity is calculated according to where ∇ h denotes the spatial discretization on the finite-difference grid with a mesh size h in each direction and t is the time step. The convective term (u · ∇ h u) κ+1/2 is calculated by using the second-order Adams-Bashforth method. Equation (15) is solved by imposing directly on u * the prescribed outer boundary conditions (periodic in the present case). At the nodes surrounding the particles (cage nodes), u * is set equal to the velocity estimated from the analytical solution.
The pressure is updated from where P h is a projector that singles out the cage nodes and n c is the unit normal directed outward from the cage. The role of the last term of equation (16) consists of enforcing, at convergence, a zero-normal-gradient condition on the cage nodes without recognizing them as internal boundaries [27]- [29], which enables us to use fast Poisson solvers. Since equation (16) is implicit, we solve it by iteration. The velocity at the end of each iteration step is calculated from

Truncation of the summation.
The expansions for the flow fields (equations (5) and (6)) obviously need to be truncated. Truncations of these series at order n = N c retains (N c + 1) 2 coefficients (P nm ,P nm ), N c (N c + 2) coefficients ( nm ,˜ nm ) and N c (N c + 2) coefficients (X nm ,X nm ), which gives a total of 3 N c (N c + 2) + 1 coefficients.

8
As already mentioned, these coefficients are updated by matching the numerical values of the pressure and of the vorticity to their analytical expressions on the cage nodes. This operation gives rise to a linear system of approximately 4(4π R 2 / h 2 ) equations, where h denotes the grid spacing. As a first approach, this should give a reasonable estimate of the truncation order N c we should use.
However, aliasing considerations do not permit to retain all these coefficients. The optimal value of N c is rather determined by numerical experience. The matrix of the linear system for the coefficients is finally rectangular, with many more rows than columns. The system is solved by the singular value decomposition algorithm, which is equivalent to a least-squares procedure when all the singular values are retained, as here [32,33].

Grid resolution near the particles.
As previously explained, the Physalis method relies on an approximate solution in the fluid regions between the particle and surrounding cage. This has the effect of introducing an error in the calculation, which can be reduced by refining the grid, thus putting the cage nodes closer to the surface of the particle.
In order to estimate the maximal grid spacing necessary for a good accuracy, it can be noted that the grid points should lie well inside the boundary layer (BL) for the Stokes solution to be valid. The BL thickness can be roughly estimated as R/ Re p , where Re p is the Reynolds number based on the particle diameter. Finally, the number of nodes per particle radius R/ h should be sufficiently larger than Re p . This criterion is the same as that applicable to a standard finite-difference solution.
The method has already been tested in numerous situations [26]- [29]. We address in the present paper the interaction between a fixed particle and stationary turbulence.

Application to a turbulent flow
2.3.1. Forcing. In order to maintain a stationary level of turbulence, a forcing must be incorporated into the Navier-Stokes equations. In the present situation, this forcing must satisfy two criteria: (i) it should act in the physical space (this is not a spectral simulation); (ii) it should vanish at the particle surface (otherwise, the forcing might induce spurious forces and torques on the object). These two constraints are satisfied by the linear forcing proposed in [34] and studied in [35]: F = Au, where F is a force per unit mass and A is a parameter kept fixed during the simulation. As shown in [35], the parameter A can be determined from a balance of kinetic energy for a statistically stationary state, which leads to where ε is the mean energy dissipation rate per unit mass and u rms is the root-mean square (rms) of the velocity fluctuations. Prescribing this parameter A and keeping it constant during a simulation are therefore equivalent to imposing a prescribed eddy turnover time scale. In order to account for this forcing, the numerical scheme described in section 2.2.2 must be slightly modified. Equation (15) is simply replaced by Since, with a stationary object as here, the forcing smoothly vanishes at the particle surface, it is consistent to apply it at all the fluid velocity nodes external to and on the cage.

2.3.2.
Procedure. The procedure we use for studying the interaction between a particle and a turbulent flow is the following.
1. The stationary turbulence is generated in the absence of particles. The initial condition is a solenoidal isotropic velocity field written as the sum of Fourier modes such that the energy spectrum is of the form E(k) ∼ k −5/3 . The Navier-Stokes equations with the linear forcing are then integrated until a statistically stationary turbulent state is reached. 2. The particle is then introduced in the flow. At this stage, the forcing cannot be applied since the velocity around the objects is a priori nonzero. The parameter A is therefore set to zero during ten time steps, roughly the time needed for the flow to adjust around the particles. 3. A is then set to the prescribed value it had at the beginning of the simulation. The statistically stationary turbulent regime is reached in the fluid after a transient of duration close to one eddy turnover time (T e). 4. To be on the safe side, the statistics are only gathered starting at 1.5T e.
A disadvantage of the linear forcing is its intrinsic instability, in the sense that any mean flow is amplified under its effect. In a one-phase flow, this problem can be easily overcome: by replacing the fluid velocity by the velocity fluctuation in the force expression (F = A(u − u), where u is the instantaneous space average of the fluid velocity in the computational domain), the mean flow remains strictly stationary. This property is very useful, since the value of the mean flow can then be imposed simply by choosing suitable initial conditions. However, this trick cannot be used in the presence of particles, since the forcing would not vanish at their surfaces. In our case, the introduction of the particle in the flow (step 2 of the procedure) leads to a symmetry breaking (most likely of purely numerical origin), which induces a small mean flow u = 0, which gradually grows under the effect of the forcing. Therefore, we stopped the simulations after a few T e to make sure that u remained small. In practice, u was always smaller than 4% of u rms , which appears to be small enough to have negligible effects. Steps 2-4 were repeated several times, starting from different initial turbulent flows, and the statistics were accumulated over about 40T e in total.
This has been possible thanks to the following optimization of the code. The idea is to avoid the coefficient iteration (see step 4 in section 2.2.1) by using an explicit procedure with small time steps. The coefficients are then calculated from the computed pressure and vorticity; they are used to generate velocity boundary conditions at the particle surface, which are used to advance the calculation by one step, and so on. This explicit procedure needs very small time steps for being accurate, which is also time consuming. A good compromise is to iterate (step 4 in section 2.2.1) every ten time steps only 5 . In the system considered here, this method was found to be as accurate as the standard one (iterations at every time step) and led to a reduction of the computation time by a factor close to 5.

Numerical and physical parameters.
We consider in the present paper the interaction between a single particle and a turbulent flow. The object is kept fixed (no translation and no rotation). As already mentioned, the mean flow can be reasonably considered as negligible. The computational domain is a cubic box of linear size = 16R (R is the particle radius). The grid resolution is 128 3 , which results in a number of nodes per particle radius of 8. With this grid, the  18 20 maximum distance from the particle surface at which the Stokes approximation is used is less than R/8. It is worth mentioning that, in the case of a uniform steady flow past a fixed sphere, it has been shown that the Stokes solution as applied here does provide a correct approximation of the velocity field at distances up to R/7 from the object for a Reynolds number, based on the particle diameter and the incident velocity, of 50 [29]. The mesh size h is close to 0.5η, and the Courant number based on the instantaneous maximal velocity is always <0.4. The order of truncation of the analytical solution, N c , is equal to 3, which corresponds to retaining 46 coefficients. This number may be compared with the number of cells covering the particle surface, which may be estimated as 4π(R/ h) 2 ∼ 800. This would be the order of the number of degrees of freedom needed to represent the particle with a comparable accuracy in a conventional finite difference or finite element method (e.g. an IB approach). Conversely, 46 degrees of freedom per particle in such a context would be equivalent to using about 2 nodes per radius, which would be totally insufficient to achieve even a low accuracy.
The particle is located at the center of the domain and periodic boundary conditions are imposed in each direction. We therefore study an infinite set of particles arranged on a cubic lattice rather than an isolated object. The particles interdistance is equal to 16R, which corresponds to a volume fraction α p ≈ 10 −3 . The value of A(R 2 /ν) is set to 1.
With these values of the parameters, the particle Reynolds number Re p ≡ 2Ru rms /ν is close to 20, and the particle diameter lies in the inertial range. The values of the turbulent lengthscales and of the Reynolds numbers are gathered in table 1.

Results
We investigate here the local modification of turbulence due to the presence of the particle (see sections 3.1 and 3.2) and characterize its angular slip velocity (section 3.3).

One-point statistics at different distances from the particle
We first investigate the statistics of the fluid velocity u, kinetic energy k and dissipation ε at different distances from the particle. This will give us an indication of the 'region of influence' (RI) of the object, the size of which will be compared to the thickness of the viscous layer adjacent to the particle to which we loosely refer as 'boundary layer' (BL). Strictly speaking, in the present unsteady and relatively low-Reynolds-number flow, it may be inappropriate to think of a BL in the conventional sense. Nevertheless, dimensional analysis suggests to estimate the thickness of this viscous layer by BL ∼ R/ Re p so that with the values of our parameters, we find that BL ≈ 0.22R. We may expect BL so defined to characterize the order of magnitude of the fluid layer within which viscous effects dominate over inertial ones. This result will be confirmed in section 3.2. We return to this point in section 4. The results will be given as a function of α, the relative distance from the particle surface in units of R, defined as where r is the distance from the particle center. The statistics for each value of α will be calculated by accounting for all the grid points belonging to the shell r ∈ [R(1 + α) − h/2; R(1 + α) + h/2]. We first calculate the statistics of the Cartesian components of the fluid velocity, u x , u y , u z . Since the flow is isotropic, these three components can be treated equally. The probability density functions (PDFs) of u i (i ∈ x, y, z) for different values of α are plotted in figure 2(a). As expected because of the no slip condition, the distributions get narrower as α decreases. For α 3, the one-phase flow statistics (dashed line) seem to be reasonably recovered. This is clearer in figure 2(b), which shows the centered moments of orders 2, 3 and 4 of the Cartesian components u i . The variance σ 2 u increases with α, extrapolating to 0 when α → 0 and reaching its asymptotic one-phase flow value for α ≈ 3. The skewness, (u i − u i ) 3 /σ 3 u i , is always equal to zero up to the expected numerical accuracy, reflecting the distributions' symmetry. The flatness, (u i − u i ) 4 /σ 4 u i , is close to 3 (Gaussian statistics) for α 3 and is larger close to the object. Thus, the turbulence is significantly modified by the particle up to α ≈ 3. It is tempting to compare these results with the analytical solutions available around fixed spherical objects at Re p = 0 (Stokes flow) and at infinite Re p (potential flow). If one defines in these cases the thickness of the RI of the particle as the distance to its surface at which the velocity components are equal to 99% of their values at infinity, one obtains St ∼ 80R for a Stokes flow and pot ∼ 3R for a potential flow. We stress that in this analysis one then deals with a stationary velocity, not with velocity fluctuations as in the present computations. Anyway, the extent of the RI of the particle that we find conforms much more closely to the latter behavior than to the former, which is expected in view of the importance of fluid inertia in the present conditions. Because of the spherical symmetry of the system, it is interesting to distinguish the statistics of the radial component of u, u r , from those of its tangential component, u t . Figure 3 shows the rms of these quantities as a function of α, σ u r and σ u t ≡ (σ 2 u θ + σ 2 u ϕ )/2 . Both tend to zero as α → 0 and are in agreement with the one-phase flow value away from the object. However, the two quantities behave differently for α ∈ ]0; 2]: in this region σ u r is always smaller than σ u t . This reflects the fact that the radial velocity component is more strongly impeded by the presence of the particle than the tangential one.
We now investigate the statistics of quadratic quantities, the kinetic energy k and the dissipation ε. Figure 4 shows the PDF of ln k at different distances from the object. For α 3, the distribution calculated in single-phase flow is satisfactorily recovered. This distribution is p(ln k) = (2/ √ π σ 3 ) exp (3/2) ln k − (1/σ 2 ) exp(ln k) , which can be deduced from the fact that the velocity fluctuations are Gaussian with zero velocity and variance σ 2 , by using the change of variables: p(ln k)dln k = p(u)du. The mean value of k decreases as one approaches the surface of the particle, which is also reflected in figure 2(b) (velocity variance), as expected. Figure 5 shows the dissipation statistics. At large distances from the object, the one-phase flow statistics are recovered ( figure 5(a)). These statistics are only approximately log-normal, probably because of the low Reynolds number of our simulations. The mean dissipation is drastically enhanced close to the particle (see figure 5(b)) owing to the large velocity gradients caused by the no-slip condition [24]. To summarize, the turbulence is modified around the particle up to ≈3R from its surface. This defines an RI of the object, that is more than 10 times thicker than the viscous layer. One might argue that our forcing tends to overestimate the size of the RI, because it is smaller close to the object. It is therefore interesting to compare our results with other observations. Burton and Eaton [24] have presented some results on the mean kinetic energy and dissipation at different distances from a fixed spherical solid particle in decaying turbulence. In their paper, the particle is much smaller than ours (R ∼ η), but its Reynolds number is comparable to ours (Re p 20). They found that for Re p = 20, the RI is of order 6R, a value larger than ours. The authors also  observe that the RI gets larger when Re p decreases as one would expect on the basis of the large RI of a particle at zero Reynolds number mentioned before. Since the Reynolds numbers of our calculation and that of [24] are comparable, the difference in the RI size must be due to the different particle radii in the two simulations. Anyway, the linear forcing does not seem to overestimate the RI thickness.

Hydrodynamic forces and torques on spherical surfaces
We investigate in this subsection the local disturbance of the fluid due to the particle by adopting another point of view: the hydrodynamic forces and torques acting on spherical surfaces are  calculated. In particular, it will be interesting to study the way these quantities tend to the force and torque acting on the object as one approaches its surface. For each value of α 0, we calculate F α , the force acting inward on the surface S α defined by r = R(1 + α), by integration of the pressure and viscous stresses: with where τ i, j = µ(∂ i u j + ∂ j u i ) is the viscous stress tensor andê r is the radial unit vector. The torque T α is calculated from Since the moments of pressure forces on spherical surfaces are equal to zero, T α is purely viscous. It will be recognized that F 0 and T 0 are the hydrodynamic force and torque on the particle, which are given directly from the coefficients of the analytical solution according to equations (11)- (14). The pressure on the surfaces r = R(1 + α) is evaluated by triquadratic interpolation of the pressure field. The velocity gradients on these surfaces are calculated by a trilinear interpolation scheme. Figure 6(a) shows the rms of the force components F α,i (i ∈ {x, y, z}) as a function of α. As in the case of the Cartesian velocity components, the three components of the force can be treated on an equal footing due to the isotropy of the system investigated. As expected, F α,rms , F (p) α,rms and F (v) α,rms continuously tend to F 0,rms , F (p) 0,rms and F (v) 0,rms as α tends to 0. Whereas near the particle the pressure and viscous contributions are comparable, which is qualitatively similar to Stokes flow behavior, as the distance from the particle increases, the pressure contribution becomes dominant, indicating once more the importance of inertial effects: for α,rms /10. It also appears that our estimate of the BL thickness as the distance from the particle where viscous effects become smaller than inertial effects, namely ∼0.22R, was reasonable. It is worth noting that the two contributions exhibit different behaviors: F (p) α,rms is an increasing function of α, whereas F (v) α,rms is non-monotonic, decreasing for increasing α ∈ [0; 1]. As a consequence, the rms of the total force is roughly constant in this region. Further away from the particle, one would expect the traction to be statistically uniform so that F (p) α and F (v) α would be proportional to r 2 , i.e. (1 + α) 2 . The figure shows that this expectation is indeed approximately verified although, very far away from the particle, the artificial periodicity of the system begins to make itself felt. Another quantity of interest is the correlation between the force components at different distances from the particle and their value at the surface of the object, C 0,α ≡ F 0,i F α,i / (F 0,i ) 2 (F α,i ) 2 . C 0,α is plotted as a function of α in figure 6(b). For the pressure component, this correlation is always positive and it becomes close to zero for α 3, i.e. outside the RI evidenced in the previous subsection. On the other hand, the viscous component quickly becomes anti-correlated as one moves away from the particle, with C 0,α reaching a value of about −0.8 for α ∼ 2. The torque statistics will provide us a physical interpretation of this behavior. Since the viscous contribution is subdominant for α > 1, the correlation of the total force is dictated by the pressure contribution one.
The statistics of the torque components T α,i are represented in figure 7. As shown in figure 7(a), T α,rms continuously tends to T 0,rms as α → 0. T α,rms is roughly constant for α ∈ [0; 0.75]. Further away from the particle, one would expect that T α,rms scales as r 3 (one-phase flow scaling). This scaling seems to hold for 1 + α > 2.5, but it is not clearly visible because of the effects of the system periodicity, felt for 1 + α > 4. The correlation of the torque components T α,i with their value on the particle surface is plotted as a function of α in figure 7(b). The torque components T α,i decorrelate very quickly as α increases, and are anti-correlated with the components T 0,i for α ≈ 2. Further away, the correlation tends to zero. As expected since the torque is purely viscous, this behavior is similar to the viscous force one (see figure 6(b)). However, the physical interpretation is clearer from the torque point of view. Indeed, the quasi-anti-correlation of T 2,i with T 0,i means that at a distance of 2 radii from the object the fluid rotates in the opposite direction with respect to the one at the particle surface. This length is actually ≈ L/2, the order of magnitude of the large-scale structures. A visualization of the flow suggests that when a side of these structures is in contact with the particle, it has the tendency to turn around the object, whereas its tail (approximately located at 2 radii from the particle surface) turns in the opposite direction.

Characterization of the angular slip velocity
We have already mentioned in the introduction section the problem of the definition of a meaningful slip velocity for an extended particle. This quantity appears in particular in the Maxey-Riley-Gatignol (MRG) equation of motion of a solid sphere of mass m p , which is widely used in Eulerian-Lagrangian simulations of disperse flows when the local particle Reynolds number is small [36,37]: Here v p is the particle's velocity, u is the fluid velocity at the object location and m f is the mass of fluid displaced by the particle. This equation is, in principle, valid if Re p 1 and if the nonuniformity of the flow at the particle scale is weak. The first term in the rhs of equation (24) is the added mass force, namely the force required to accelerate the portion of fluid carried along with the particle. The second term is the viscous drag. The third one is the history force, which accounts for the interaction between the object and its wake (this is actually the unsteady portion of the drag). The last two terms are the buoyancy force and the force due to the fluid stress around the particle. The Faxén corrections terms, ∝ ∇ 2 u [38], account for the non-uniformity of the flow. They would be negligible if the flow was locally approximately uniform. It would be desirable to use the results of simulations such as the present one to explore the generalization of the MRG equation beyond the range for which it has been derived. Unfortunately, this does not seem feasible in view of the complexity of the equation and the difficulty of separating the total hydrodynamic force given by the numerical simulation into the various components appearing in it. A similar difficulty would be found starting from model equations applicable to large Reynolds numbers, such as the ones developed in [39] or in [40], which contain added mass and lift terms.
A more modest objective can be pursued by focusing on the torque acting on a spherical particle, namely [37,41] where I p denotes the moment of inertia of the object, p is its angular velocity, I f is the moment of inertia of a sphere of radius R and of the same density as the fluid, and ω f /2 is the fluid angular velocity 'seen' by the particle. This equation is valid under the same hypotheses as the previous equation (24). The angular momentum is resisted by the fluid viscosity through a drag torque and two history torque components. The last term is the angular acceleration of the fluid. The torque expression is simpler than the force one, since there is no component associated with pressure, gravity or added mass: the hydrodynamic torque is a purely viscous effect. If the relative angular acceleration and that of the fluid are negligible, the hydrodynamic torque reduces to [42] In this limit, the torque acting on the object is proportional to the relative angular velocity, which is a consequence of the linearity of Stokes flow. We use this property to characterize the angular slip velocity of the particle. For this purpose, let us adopt the reasonable assumption that ω f can be approximated as the fluid vorticity averaged in a spherical shell surrounding the particle. We have two characteristic lengths at our disposal: (i) the BL thickness; (ii) the thickness of the RI (see section 3.1). We therefore plot in figure 8 the components of the torque acting on the object, T i (i ∈ {x, y, z}), versus the vorticity components averaged in the RI, ω i RI ( figure 8(a)), and the vorticity averaged in the BL, ω i BL ( figure 8(b)). There is no clear correlation between T i and ω i RI . On the other hand, T i is clearly proportional to ω i BL (correlation coefficient greater than 0.995), and the proportionality coefficient is the one of equation (26) (see the dashed line in figure 8(b)). The BL thickness is only defined in order of magnitude, but we have checked that the correlation coefficient was always greater than 0.995 if one averages vorticity over shells of thickness ±20% greater than the values used in the figure.
We thus conclude that, in the present range of parameters, the angular slip velocity can be defined in terms of the vorticity averaged over a shell concentric with the particle and extending to the edge of the viscous layer. This result would of course be trivial if the averaging had to be carried out very close to the particle surface where the Stokes flow approximation is applicable. What makes it interesting is that it seems to hold at a much larger distance from the particle than one might expect.

Conclusion
We have studied the interaction between a fixed solid spherical particle of size lying in the inertial range and stationary turbulence without mean flow. We have first evidenced an RI of the object on the fluid, which roughly lies in [R; 4R]. This RI is the one in which the turbulence statistics are modified by the presence of the object. It is more than ten times as thick as the BL, defined as the region in which viscous effects dominate over inertial ones. Another way of characterizing the RI is to measure the hydrodynamic force on spherical surfaces centered on the particle center. The pressure component of this force, dominant with respect to the viscous contribution outside the viscous layer, is decorrelated with that acting on the particle iff the surface is outside the RI. The viscous force and the (purely viscous) torque acting on such surfaces decorrelate much faster with increasing distance from the solid/fluid interface. In particular, at 2 radii from it the fluid roughly rotates in a direction opposite to the one it has in the vicinity of the particle. Finally, the fact that the object is fixed enabled us to isolate the viscous drag torque and thereby to define the effective angular slip velocity of the particle.
This work has direct implications in the modeling of the turbulent transport of finitesize objects. In particular, the definition of the angular slip velocity is a first step toward the derivation of equations for the rotational motion of particles in turbulence, which is coupled to their translational motion, for instance through the lift force. The measurement of the size of the RI, RI , is also relevant for the physics of particle-laden flows. Indeed, the comparison between RI and the mean distance between the objects suspended in the fluid P enables us to predict the flow regime (one-, two-or four-way coupling) that holds. For instance, if P < 2 RI , then the RI of different particles will overlap and the four-way coupling will hold. The major difference between our study and the real situation of a particle moving freely in a turbulent flow is the fact that no mean wake can develop in our case because of the absence of mean flow. This made the interpretation of our measurements clearer, but the same study must be performed in the case where the object is not fixed. This is the subject of ongoing work, and will enable a direct comparison between the numerical results and the experimental data now available for integral scale objects [14].
The Reynolds number of our simulations is rather low, but it might nevertheless be of practical significance as the Reynolds number of the relative motion of a particle transported in a turbulent flow might be moderate even in the case of a relatively large particle provided, for example, that the particle-fluid density is not too large. This statement of course assumes that a meaningful particle-fluid slip velocity can be defined, otherwise the very notion of particle Reynolds number loses meaning. In the absence of a systematic study varying the particle size and the Reynolds number, one can reasonably hypothesize that, ultimately, the average RI of the object cannot be reduced below that of a sphere in potential flow, a limit that seems already essentially attained at the Reynolds number of the present simulations. A significant feature that would appear with increasing Re p , however, is a wake behind the particle (be it stationary or moving), a new feature that we could not examine as explained before. This wake would be unsteady, which might induce complex flow phenomena such as induced vortex shedding by trailing particles, unsteady lift forces and others.
The generalization of the equations of motion of solid objects valid in flows at vanishing or infinite Reynolds number to turbulent situations is a difficult problem. The very notion of BL on which one might try to build such a generalization might not be relevant as a particle transported in a turbulent flow is at every moment in unsteady motion subject to strongly 19 variable acceleration. It appears that studies such as the present one are a necessary step in elucidating the physics of the interaction between finite-size particles and turbulence, which is of crucial importance for the improvement of the modeling of particle-laden flows.