Particle Flows around an Intruder

Particle flows injected as a beam and scattered by an intruder are numerically studied. We find a crossover of the drag force from Epstein's law to Newton's law depending on the ratio of the speed to the thermal speed. These asymptotic laws can be reproduced by a simple analysis of the collision model between the intruder and the particle flows. We show the existence of turbulent-like behavior of the particle flows behind the intruder from the second invariant of the velocity gradient tensor and the relative mean square displacement for large speed regime. We also characterize the scattering of particles by the intruder using the angular distribution function of scattered particles.

Particle flows injected as a beam and scattered by an intruder are numerically studied. We find a crossover of the drag force from Epstein's law to Newton's law depending on the ratio of the speed to the thermal speed. These asymptotic laws can be reproduced by a simple analysis of the collision model between the intruder and the particle flows. We show the existence of turbulent-like behavior of the particle flows behind the intruder from the second invariant of the velocity gradient tensor and the relative mean square displacement for large speed regime. We also characterize the scattering of particles by the intruder using the angular distribution function of scattered particles.
Introduction.-To know the fluid flows around an intruder depending on the Reynolds number is a fundamental problem [1,2]. When the flow speed is low, the drag force acting on a spherical intruder obeys Stokes' law in which the drag is proportional to the speed, the viscosity of the fluid, and the radius of the intruder. Whereas the drag force satisfies Newton's law in the high Reynolds number, in which the drag is proportional to the square of the moving speed and the cross section.
We believe that the drag force still satisfies Stokes' law even for low speed molecules [3][4][5][6][7]. It is known that Stokes' law can be used only for systems in the zero Knudsen number limit, the ratio of the mean free path to the intruder size [8][9][10][11][12][13]. The correction to Stokes' drag for rarefied gases in the low Kundsen number is theoretically confirmed by the kinetic theory [8][9][10][11][12][13]. Whereas the drag force acting on slowly moving intruders in the large Knudsen number satisfies Epstein's law [14,15]. It is known that the drag law depends on the distance from the boundary of a container of particles through the simulation [3] and the theory of the fluid mechanics [16]. We also note that the Kármán vortices have already been observed in molecular dynamics (MD) simulations [7,17].
In this Letter, we numerically study the drag force acting on a spherical intruder in particle flows by controlling the ratio of the injected speed of the particles to the thermal speed which is proportional to the sound speed in terms of the MD. We characterize the turbulent-like be- * e-mail: takada@go.tuat.ac.jp havior of particle flows behind the intruder with the aid of the second invariant of the velocity gradient tensor and the relative mean square displacement. We also characterize the scattering of particles by the intruder using the angular distribution function of scattered particles.
Model.-Let us explain our model and setup of our simulation. The system consists of two parts: one is a fixed intruder and another is the mobile particles as shown in Fig. 1. The intruder is made by one core particle whose diameter is D c , and N s identical particles whose diameters are d on the surface of the core particle. We examine three sizes of the intruders as (D c /d, N s ) = (5, 144), (13,744), and (28,1600). Then, the diameter of the intruder is given by D = D c +2d. The intruder is fixed at the origin. We simulate the system containing N = 30000, 101250, and 810000 monodisperse mobile particles whose mass and diameter are m and d, respectively. The collisions are always elastic. We examine various initial volume fractions ranged from φ = 0.20 to 0.55. Before starting our simulation, the mobile particles are thermalized with the temperature T , and are moved with the translational speed V . The equation of motion of i-th particle is given by where we have introduced r ij = |r i − r j |, the spring constant k, and the step function Θ(x), i. e., Θ(x) = 1 for x ≥ 0 and Θ(x) = 0 otherwise. Initially we confine the mobile particles in tubes whose radii are R = 10d for D c = 5d and 13d and R = 30d for D c = 28d, respectively. We assume that the interactive force between the wall of the tube and the particles is identical to the interparticle force.  When the mobile particles collide with the intruder, the particles reflect at random with the temperature T w , where we set T w = T for simplicity. We examine two cases for the systems behind the intruder: one is a free scattering case, and another is a confined case, where the scattered particles are still confined in the tube. The used parameters are listed in Table I. In the following, all quantities with the superscript * represents the quantities nondimensionalized by m, d, and k. We use the time-averaged drag force in a steady state, and fix the speed V * = 0.1. We have verified that the results of our simulation for V * = 0.1 are consistent with those of the hard-core particles.
Results.- Figure 2 presents the results of the dependence of the drag force against the normalized colliding speed by the thermal speed v T ≡ 2T /m for both the free scattering and the confined cases. Both results agree well each other because the drag force is determined from the impulses by collisions of particles in front of the intruder. We note that our obtained drag forces are proportional to (D + d) 2 , i. e., the collision cross section [20,36]. The drag force is proportional to the square The former can be understood by a simple collision model. Because each momentum change in the flow direction for a collision is ∆p = 2mV cos 2 θ where θ is the angle between the intruder and the colliding particle as shown in Fig. 3, and the collision frequency is Ω c = (π/4)n(D+d) 2 V with the number density n = 6φ/(πd 3 ), the force acting on the intruder is given by which agrees well with the simulation for V > v T (see Fig. 2). Whereas the front of the beam of mobile particles diffuses before the collision for V /v T ≪ 1 as shown in Fig. 4(a). This expansion speed of the diffusive front is approximately given by V p = √ V v T as shown in Fig.  4 which agrees well with the simulation results for V v T as shown in Fig. 2.  (1) and (2), respectively. Let us characterize the particle flows behind the intruder for large V . We introduce the second invariant of the velocity gradient tensor [37]. Here, we adopt Einstein's rule for i and j where duplicated indices take summation over x, y, and z. Figure 5 shows the contour of Q = 0, where the field is coarse-grained with the scale w = 2d for visibility [38]. The vortex rich regions Q > 0 emerge behind the intruder. This behavior is similar to that observed in turbulent flows induced by an intruder Let us study how the particles are scattered after collisions with the intruder. We only focus on the relative motion of the particles which collide with the intruder almost simultaneously through the mean square displace- ment with , where we only select two particles (i and j) within the interval |t i − t j | < ∆t th ≡ 10 m/k with the collision times t i and t j with the intruder for i-th and j-th particles, respectively. Note that t c in Eq. (3) is larger time of t i and t j . Figure 6 shows the super-ballistic behavior ∆(t) ∼ t 2.2 . This behavior is analogous to the relative motion of two tracer particles in turbulent flows, which is known as Richardson's law ∆(t) ∼ t 3 [39][40][41]. We have checked that this result is insensitive to the choice of ∆t th in the range m/k ∆t th 20 m/k. The exponent 2.2 for the super-ballistic ∆(t) in our system seems to be much smaller than 3 for ∆(t) in turbulent flows. This is because the fluid flow behind the intruder in our system is not a fully developed turbulence but is a weak turbulence.
Next, we consider the scattering angle distribution of the mobile particles. The azimuthal angle is stored when For V /v T ≫ 1, the angular distribution of scattered particles has a sharp peak around the opening angle, which is the half of the apex angle of the cone of the beam scattered after collisions with the intruder [32,35] as shown in Fig. 8. We also note that this opening angle can be explained by a phenomenology as in Ref. [32]. Because we use repulsive and elastic particles, the opening angle is expressed as for D < 2R [32]. When we substitute D = 7d and R = 10d into Eq. (4), we obtain θ 0 = 0.50 (rad), which roughly agrees with the simulation results (see Fig. 7). For V /v T ≪ 1, on the other hand, the particles are scattered in various direction as shown in Fig. 7.
Discussion.-In our setup, we have not observed Stokes' flow around the intruder. This is because the mobile particles are packed in the tube as an initial condition. As already explained, the expansion wave exists and the drag force acting on the intruder, known as Epstein's law, is determined from this wave. To realize Stokes' flow, we should simulate the situation for V /v T ≪ 1 [3][4][5].
Conclusion.-We have numerically investigated the particle flows injected as a beam and scattered by a spherical intruder. We found the crossover from Epstein's law to Newton's law, depending on the ratio of the speed to the thermal velocity, which can be explained by a simple collision model. The turbulent-like behavior has been also observed, where the mean square displacement is super-ballistic. The scattering of particles by the intruder is also characterized by the angular distribution function of scattered particles.

Supplemental Materials: Particle Flows around an Intruder
Effect of the geometry of the intruder In this section, let us check the dependence of the boundary condition between the intruder and the mobile particles. Here, we refer the intruder introduced in the main text as "thermal and bumpy", because the mobile particles reflect at random when they collide with the intruder, and the small particles are attached on the surface of the core particle. We prepare three different types of the intruders: (i) "reflective and smooth," (ii) "reflective and bumpy," and (iii) "thermal and smooth." The first intruder consists of only one core particle, and the reflection between the intruder and the mobile particles is elastic. The second one consists of the core particle and the small particles on its surface, where we have used the boundary condition as that used in the first case. The third one consists of only one core intruder, where we have used the thermal boundary condition as used in the main text. Figure S1 plots the results of the drag forces under various boundary conditions. The results clearly indicate that the boundary condition on the intruder is not important for the drag force. This is because the force acting on the intruder is almost determined from the region in front of the intruder as explained in the main text.

Vorticity
In this section, we plot the vorticity induced by the scattering of the intruder. In the following, we convert the Cartesian coordinate (x, y, z) into the cylindrical coordinate (r, ϕ, z), in which the common z-axis is parallel to the flow direction. Here, the velocity components in the cylindrical coordinate (v r , v ϕ ) are expressed as where (v x , v y ) are the velocity components represented by the Cartesian coordinate. Because the twisted structure of the flow field is not observed in our simulation, we focus on the flow structure in rz-plane. Let us introduce the vorticity in rz-plane as (S2) Figure S2 shows

Legendre coefficients
In this section, we analyze the angular distribution in terms of the Legendre polynomials as where P ℓ (x) is the Legendre polynomial of degree ℓ [S2], and A ℓ is the coefficient defined by (S4) Figure S3 shows the coefficient for each degree ℓ. The coefficients converge slowly to zero for large ℓ. For V /v T ≫ 1, A ℓ has a peak at ℓ = 2.

Volume fraction dependence of the sound speed
In this section, we show how the sound speed depends on the volume fraction. Because the volume fraction of the mobile particles is finite, the equation of state deviates from that for the ideal gas. Thus, the equation of state at finite density is fitted by the Carnahan-Starling equation [S3] where p is the static pressure. For the adiabatic process, the first law of thermodynamics becomes where L 3 is the volume, C V = 3/2 is the heat capacity at constant volume. Substituting the equation of state (S5), the following quantity is conserved: Then, the sound speed in the adiabatic process is given by with f (φ) = 5 + 10φ − 3φ 2 − 24φ 3 + 37φ 4 − 22φ 5 + 5φ 6 3(1 − φ) 6 .
(S9) Figure S4 shows the volume fraction dependence of the coefficient f (φ). When we use the expansion speed V ′ p = c(φ)V , the drag force for the low V /c cannot be captured as shown in  [S1] See the ancillary file for the movie for V /vT = 1.