{'fraction_non_alphanumeric': 0.03203717290291025, 'fraction_numerical': 0.015488709545936251, 'mean_word_length': 4.348811859603226, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A numerical study of the motion of a neutrally buoyant cylinder in two dimensional shear flow

In this paper, we investigate the motion of a neutrally buoyant cylinder of circular or elliptic shape in two dimensional shear flow of a Newtonian fluid by direct numerical simulation. The numerical results are validated by comparisons with existing theoretical, experimental and numerical results, including a power law of the normalized angular speed versus the particle Reynolds number. The centerline between two walls is an expected equilibrium position of the cylinder mass center in shear flow. When placing the particle away from the centerline initially, it migrates toward another equilibrium position for higher Reynolds numbers due to the interplay between the slip velocity, the Magnus force, and the wall repulsion force.


Introduction
The problem of particle motions in shear flows is crucially important in many engineering fields such as the handling of a fluid-solid mixture in slurry, colloid, and fluidized bed. Segré and Silberberg [27,28] experimentally studied the migration of dilute suspensions of neutrally buoyant spheres in a tube Poiseuille flow. The particles migrate away from the wall and the centerline to accumulate at about 0.6 of the tube radius from the centerline. The experiments of Segré and Silberberg [27,28] have had a large influence on fluid mechanics studies of migration and lift of particles. Comprehensive reviews of experimental and theoretical works have been given by Brenner [2], Cox and Mason [6], Feuillebois [11] and Leal [18].
Among the theoretical studies of the neutrally buoyant particle migration in linear shear flow, Bretherton [3] found an expression for the lift force per unit length on a cylinder in an unbounded two-dimensional linear shear flow at small Reynolds number. Saffman's lift force [26] on a sphere of radius a in an unbounded linear shear flow with shear rate G is F s = 6.46ρV a 2 (Gν) 1/2 = 6.46ρνaV (Re p ) 1/2 where ν is the kinetic viscosity of the fluid, ρ is the density of the fluid, Re p = Ga 2 /ν is the particle Renolds number, and V is the slip velocity of the sphere. In a bounded linear shear flow, Ho and Leal [16] examined the motion of a rigid sphere with inclusion of the inertia effects at small Reynolds numbers by a regular perturbation method. The sphere reaches a stable lateral equilibrium position which is the midway between the walls. Vasseur and Cox [29] also obtained the same stable lateral equilibrium position. Ho and Leal require that Re p /κ 2 1 which is more restrictive than the one Re p /κ 1 required by Vasseur and Cox where κ = 2a/H is the confined ratio, H being the distance between two walls. Direct numerical simulations have been used for understanding particle motion in shear flows. Feng et al. [10] investigated the motion of neutrally buoyant and non-neutrally buoyant circular particle in plane shear and Poiseuille flows using a finite element method and obtained qualitative agreement with the results of perturbation theories and of experiments. The numerical results of a neutrally buoyant circular cylinder in a shear flow of Re p = 0.625 have been discussed in details. The cylinder migrates back to the midway between two walls due to the wall repulsion at the small Reynolds number. They have suggested that that three factors, namely the wall repulsion due to a lubrication effect, the slip velocity, and the Magnus type of lift, are possible responsible for the lateral migration. Ding and Aidun [8] studied numerically the dynamics of a cylinder of circular or elliptic shape suspended in shear flow at various particle Reynolds number. They obtained the transient from being rotary to stationary as the particle Reynolds number is increased for an elliptic cylinder. For the cases of the circular cylinder, the effect of the two walls on the rotation speed ω has been studied. For the confined ratio κ = 0.5, Ding and Aidun reported |ω|/G ∝ Re −0.28 p . Similar result,|ω|/G ∝ Re −0.25 p , has been observed experimentally by Zettner and Yoda [31].
In this paper, we first discuss the generalization of a distributed Lagrange multiplier/fictitious domain method (DLM/FD method) developed in [22] to non-spherical neutrally buoyant cylinders in two-dimensional shear flows and to the cases with zero angular velocity as a constraint. Similar DLM/FD methods for freely moving neutrally buoyant particle has been successfully applied, in [20,23,30], to simulate particulate flow in two and three dimensions with neutrally buoyant particles. But for the cases of a neutrally buoyant elliptic cylinder in two dimensional flows, the methodology has not been fully validated yet. We have validated the numerical method by comparing with the computational results in Ding and Aidun [8] for a cylinder of circular and elliptic shape and the experimental results in Zettner and Yoda [31] for a circular cylinder. On the wall effect on the angular velocity of the circular cylinder, we have obtained |ω|/G ∝ Re −0.2771 p for the confined ratio κ = 0.5 which is in a good agreement with the results obtained by Ding and Aidun [8] and Zettner and Yoda [31]. We have also studied the wall effect on the angular velocity of the elliptic cylinder which is more complicated due to the non-circular shape.
Concerning the equilibrium position of a neutrally buoyant circular cylinder in shear flow, recent studies by Cherukat, McLaughlin and Dandy [5] and Kurose and Komori [17] focus on lift and drag on a stationary sphere in unbounded linear shear flow. The equilibrium positions have not been studied in these works. Feng and Michaelides [9] have investigated the equilibrium positions of non-neutrally buoyant circular cylinders in two-dimensional shear flow. In their simulations, the density ratio between the solid and fluid is between 1.005 and 1.1. The equilibrium heights of their lightest circular cylinder (the density ratio of 1.005) are far below the centerline. We have obtained that the cylinder stays in the middle between two walls as expected when placing it there initially. But when the initial position of the mass center of a circular cylinder is away from the centerline, the equilibrium position depends on the particle Reynolds number Re p and the confined ratio κ. For those placing away from the centerline initially, the circular cylinder migrates back to the centerline for Re p < Re p,c where Re p,c is the critical value. For Re p > Re p,c , the equilibrium position is between the wall and the centerline. The critical particle Reynolds number is increasing when increasing the confined ratio. Concerning the Magnus lift effect on the equilibrium position, we have added a constraint, zero angular velocity, to the motion of a circular cylinder and obtained that the equilibrium position of the circular cylinder moving with zero angular velocity is lower than those of freely moving cylinder when both are away from the middle. These results show that the Magnus lift does play a role as expected. Also from the computed slip velocity of the circular cylinder, it shows that the circular cylinder lags the fluid, at least for Re p > Re p,c , which means that there is a force pushing the cylinder toward the wall (see Fig. 1 for the setup of the boundary conditions). Hence the equilibrium position of the cylinder is up to the interplay between the slip velocity, the Magnus lift and the wall repulsion force. The content of this paper is as follows: In Section 2 we introduce a fictitious domain formulations of the model problem associated with the neutrally buoyant long particle cases; then in Section 3 we discuss the time and space discretization and in Section 4 we present and discuss the numerical results.

A fictitious domain formulation of the model problem
A fictitious domain formulation with distributed Lagrange multipliers for flow around freely moving particles and its associated computational methods have been developed and tested in, e.g., [12,13,24,25]. For the cases of neutrally buoyant particles, similar methodologies have been developed in [20,22,23,30]. But for the cases of a neutrally buoyant elliptic cylinder in two dimensional flows, the methodol-  ogy has not been discussed and fully validated yet. In this paper, we first discuss the formulation for the case of a neutrally buoyant elliptic cylinder and then present the numerical results to validate the methodology.
Let Ω ⊂ IR 2 be a rectangular region filled with a Newtonian viscous incompressible fluid (of density ρ and dynamic viscosity µ) and containing a freely moving neutrally buoyant rigid particle B centered at G = {G 1 , G 2 } t of density ρ. The flow is modeled by the Navier-Stokes equations and the motion of the particle B is described by the Euler-Newton's equations. We define W g 0 ,p = {v|v ∈ (H 1 (Ω)) 2 , v = g 0 on the top and bottom of Ω and v is periodic in the x 1 direction}, W 0,p = {v|v ∈ (H 1 (Ω)) 2 , v = 0 on the top and bottom of Ω and v is periodic in the x 1 direction}, an inner product on Λ 0 (t) which can be the standard inner product on (H 1 (B(t))) 2 . For simple shear flow, we have g 0 = (U/2, 0) t on the top wall and (−U/2, 0) t on the bottom wall. Then the fictitious domain formulation with distributed Lagrange multipliers for flow around a freely moving neutrally buoyant particle of the elliptic shape is as follows where u and p denote velocity and pressure, respectively, λ is a Lagrange multiplier, V G is the translation velocity of the particle B, ω is the angular velocity of B, and θ is the angle between the horizontal direction and the long axis of the elliptic cylinder (see Fig.1). We suppose that the no-slip condition holds on ∂B. We also use, if necessary, the notation φ(t) for the function x → φ(x, t).
Remark 1. The hydrodynamical forces and torque imposed on the rigid body by the fluid are built in (1)- (7) implicitly (see [12,13] for details), thus we do not need to compute them explicitly in the simulation. Since in (1)-(7) the flow field is defined on the entire domain Ω, it can be computed with a simple structured grid.
Remark 2. In (3), the rigid body motion in the region occupied by the particle is enforced via the Lagrange multiplier λ. To recover the translation velocity V G (t) and the angular velocity ω(t), we solve the following equations Remark 3. To investigate the effect of the Magnus type of lift on the lateral migration of the cylinder, we have considered the cases of the cylinder freely moving in shear flow with zero angular velocity. For this special consideration, we have the following modified formulation with the modified multiplier space

Space approximation and time discretization
Concerning the space approximation of problem (1)-(7) by a finite element method, we have used P 1 -iso-P 2 and P 1 finite elements for the velocity field and pressure, respectively (like in Bristeau et al. [4]). We approximate then W g 0 ,p , W 0,p , L 2 and L 2 0 by the following finite dimensional spaces (17) and bottom of Ω and v is periodic at Γ in the (18) and bottom of Ω and v is periodic at Γ in the x 1 direction }, respectively; in (17)- (20), P 1 is the space of polynomials in two variables of degree ≤ 1.
At each time step we have a sequence of subproblems (32), (33), (35) and (39). The degenerated quasi-Stokes problem (32) is solved by a preconditioned conjugate gradient method introduced in [14], in which discrete elliptic problems from the preconditioning are solved by a matrix-free fast solver from FISHPAK by Adams et al. in [1]. The advection problem (33) for the velocity field is solved by a wave-like equation method as in [7,21]. Problem (35) is a classical discrete elliptic problem which can be solved by the same matrix-free fast solver. To enforce the rigid body motion inside the region occupied by the particles, we have applied the conjugate gradient method discussed in, e.g., [20,22,23].

The motion of an elliptic cylinder in linear shear flow
We have first considered the motion of a neutrally buoyant elliptical cylinder in linear shear flow studied by Ding and Aidun [8] for the validation purpose. Their results were obtained by the lattice Boltzmann equation. The domain of computation is Ω = [0, 5] × [0, 1], then the height is H = 1. The confined ratio is κ = 2a/H = 0.2, and the aspect ratio is AR = b/a = 0.5 where a is the length of the semi-major axis and b is the length of the semi-minor axis. The density of the fluid is ρ = 1 and the kinetic viscosity is determined by the specified value of the Reynolds number via ν = Ga 2 /Re p . The shear rate is fixed at G = 1 and U = −GH = −1 (so the moving directions of the two walls are the same as those in Ding and Aidun [8]). The mesh size is h = 1/320 and the time step is t = 0.001.
In Fig. 2, the computational results of the angle and angular velocity of the elliptic cylinder are in a good agreement with the Jeffery's solution [15] for Re p = 0 and the ones obtained by Ding and Aidun for Re p = 0.02 and 0.25, respectively. Fig. 3 shows that at various particle Reynolds numbers from 0 to 7.5, our results of angular velocity match very well with the ones obtained by Ding and Aidun [8]. The motion of the ellipse for Re p ≤ 7 is a periodic rotation with non-uniform angular velocity, while for Re p ≥ 7.5, the ellipse does not rotate at all; instead it takes a stationary orientation in shear flow (see, e.g., Ding and Aidun [8] for further details). An example of the velocity field of Re p = 8.25 around an elliptic cylinder with a stable orientation is shown in Fig. 4.
The minimum angular velocity decreases as the particle Reynolds number is increased with a nearly straight line relationship as shown in Fig. 5, where Re p,c = 7.25 is the critical particle Reynolds number above which the rotational motion is stopped. Fig. 5 shows that the period of rotation increases rapidly as the particle Reynolds number is increased close to Re p,c , and the results are in good agreements with those obtained by Ding and Aidun [8]. The further studies of the motion of a cylinder of axisymmetric shape in simple shear flow will be reported in a forthcoming paper.

The rotation of a cylinder in a simple shear flow
In, e.g., Aidun and Ding [8], the angular velocity of a circular cylinder suspended     in a simple shear flow has been studied via the direct numerical simulation. Experimental results have been also obtained in, e.g., Zettner and Yoda [31]. In this section, we focus on the wall effect of the rotation speed of a cylinder freely suspended in a simple shear flow with various confined ratios. We have first considered the cases of a circular cylinder suspended in the middle between two walls initially.  Fig. 6. At the zero Reynolds number, the rotation speed of a circular cylinder is G/2 from the Jeffery's solution. At the small particle Reynolds numbers, the ratio |ω|/G is about to converge to 0.5 when decreasing the confined ratio κ from 0.5 to 0.125 as shown in Fig. 7. When κ =0.5, the normalized angular speed |ω|/G is found to be about 0.420 for Re p ≤ 2. When κ =0.25, |ω|/G is about 0.482 for Re p ≤ 0.2. Both are very close to those obtained in Ding and Aidun [8]. When κ =0.125, |ω|/G = 0.496 is much closer to 0.5 for Re p ≤ 0.05 due to weaker effect from the walls.
When increasing the particle Reynolds numbers, the log-log plot of |ω|/G versus Re p in Fig. 7 shows that |ω|/G ∝ Re For an elliptic cylinder suspended in linear shear flow, we have studied the effect of the aspect ratio AR. The mass center of an elliptic cylinder also stays at the middle between two walls if it is placed there initially. Since the elliptic cylinder has zero angular velocity when the particle Reynolds number is larger than the critical value Re p,c , the behavior of the rotation of an elliptic cylinder is different from that of the circular cylinder. In Fig. 8, the minimal angular velocity decreases rapidly to about zero when the particle Reynolds number Re p is closer but less than Re p,c since the motion of the elliptic cylinder is about to transit into the one with a fixed orientation in linear shear flow. The minimal angular velocity decreases faster for the smaller aspect ratio shows that the slender shape is easier to reach a stable orientation in linear shear flow. But its maximal angular velocity shown in Fig. 8 have different behavior since the maximal angular velocity happens when the direction of the long axis is about perpendicular to the shear direction. For the small particle Reynolds numbers, the maximal and minimal values of the angular velocity are close to the Jeffery's solution as in Fig. 8.

The equilibrium position
In Ho and Leal [16] and Vasseur and Cox [29], they concluded that the sphere reaches a stable lateral equilibrium position which is the midway between the walls for small particle Reynolds numbers. In Feng et al. [10], the cylinder migrates back to the midway between two walls at Re p = 0.625 when placing it away from the middle between two walls. Feng et al. have suggested that that three factors, namely the wall repulsion due to a lubrication effect, the slip velocity, and the Magnus type of lift, are possible responsible for the lateral migration. In the previous section, we have obtained that the centerline is always (at least in the range we have studied in this paper) the equilibrium position for the cylinder of elliptic or circular shape when it is positioned there initially. When placing the mass center initially away from the middle between two walls, the mass center may not migrate back to the centerline.
For the cases of a circular cylinder of the confined ratio κ = 0.25 as in the previous section, the radius of the circular cylinder is a = 0.125, the density of the fluid is ρ = 1, and the kinetic viscosity is ν = 0.012. To check the effect of the length L of the channel, the initial height y init , the mesh size h and time step dt on the equilibrium height of the mass center, we have considered different sets of parameter values as indicated in Fig. 9. The final equilibrium heights in Fig. 9 are almost the same for each value of Re p considered here except those with zero angular velocity constraint. When Re p =1 and 2, the mass center of the freely moving cylinder migrates back to the middle between two walls. But for higher particle Reynolds numbers, we have   obtained an equilibrium height which is between the centerline and the wall. To find out the effect of the walls on the final equilibrium height, we have varied the distance H between two walls. The initial height y init is always 0.1 unit below the centerline. The length of the computational domain is L = 2. The equilibrium height versus the particle Reynolds number is shown in Fig. 10. As expected, the wall repulsion force is stronger for the larger confined ratio. Hence the cylinder migrates back to the midway between two walls for κ = 0.5 for the range of Re p considered here and the critical particle Reynolds number Re p,c is increasing when increasing the confined ratio. We believe that the symmetric breaking shown in Figs. 9 and 10 is not a numerical artifact. To further validate this phenomenon, we have compared ours with the results obtained in Feng and Michaelides [9] where only the non-neutrally Figure 13: A snapshot of the velocity field of Re p = 5 obtained by following the cylinder mass center which is moving to the left: freely motion (top) and zero angular velocity (bottom). circular cylinders were considered. For the results presented in Fig. 11, the radius of circular cylinder is a = 0.05, the initial height is y init = 0.1, the kinetic viscosity is ν = 0.05, and the length and the height of the computation domain are L = 2 and H = 1, respectively. The density ρ s of the non-neutrally circular cylinder is either 1.002 or 1.005. We have applied the DLM/FD method developed in, e.g., [12,13] to obtain the computational results for the cases of the non-neutrally circular cylinder. Our results are in a good agreement with the results obtained in Feng and Michaelides [9]. The transition of the equilibrium height from the one associated with a non-neutrally buoyant cylinder to that of a neutrally buoyant cylinder is consistent and support the existence of the symmetric breaking of the equilibrium height.
For the effect of the Magnus lift associated with the rotation, we have considered the cases of a circular cylinder of κ =0.25 moving with zero angular velocity. The equilibrium height of non rotating cylinder is lower than the one of the freely moving cylinder for Re p ≥ 2 as in Fig. 9. These results indicate that the Magnus lift does play a significant role as expected. The similar results concerning the Magnus lift have also been observed in Feng and Michaelides [9]. When Re p =1, the wall repulsion force is strong enough to push the cylinder back to the centerline even without the help from the Magnus lift.
We now focus on the effect of the slip velocity. Since the initial position of the cylinder is below the centerline, it is moving to the left in the lower region of the computational domain due to the given boundary condition (see Fig. 1). For getting the slip velocity, we first compute the fluid horizontal speed on the streamline through the mass center in front of the cylinder at the distance of the half of the computational domain width and then minus the horizontal speed of the disk to obtain the slip velocity. The negative slip velocity means that the cylinder speed to the left is slower than the fluid speed to the left since the both signs are negative. We then have that the cylinder lags the fluid. For Re p =1, the slip velocity of the non rotating cylinder becomes negative after a short initial transition period as in Fig. 12. The one associated with free moving circular cylinder remains positive for, at least, the first 110 time units and then oscillates about zero. The cylinder with no rotation lags the fluid. But both cylinders migrate back to the middle between two walls due to that the wall repulsion force dominates the very weak slip velocity effect. We can see freely moving circular cylinder moves toward the centerline faster than the one with no rotation does in Fig. 12. When Re p =5, the slip velocity becomes negative after the initial transition period as in Fig. 12. The cylinder for the both cases lags the fluid. The slip velocity of the one with no rotation is about two times larger than the other one. The rotating velocity of each case of free motion is about constant speed after a short initial transition period. Fig. 13 shows that the relative velocity field to the horizontal velocity of the cylinder mass center in which we can clearly see that the cylinder lags the fluid. Both cases have stronger slip velocity which creates force pushing the cylinder to the region with faster flow speed, which is the region next to the bottom wall. Without the extra help from the Magnus lift, the one without rotation is closer to the wall. Hence when the initial position is not at the middle of two walls, the balance between the effect of the slip velocity, the Magnus lift and the wall repulsion does play a role for determining the equilibrium position of the cylinder.

Conclusion
We have investigated the motion of a neutrally buoyant cylinder of circular or elliptic shape in two dimensional shear flow of a Newtonian fluid by direct numerical simulation. The numerical results are validated by comparisons with existing theoretical, experimental, and numerical results, including a power law of the normalized angular speed versus the particle Reynolds number. The rapid slow down of the normalized minimal angular speed of an elliptic cylinder is totally different from the behavior of the circular cylinder since the motion of an elliptic cylinder can transit from rotating to a fixed orientation when the particle Reynolds number is increased beyond the critical value. The midway between two walls is an expected equilibrium position of the cylinder mass center in shear flow. But when placing the particle away from the centerline initially, it migrates toward another equilibrium position between the wall and the centerline for higher Reynolds numbers which is caused by the interplay between the slip velocity, the Magnus force, and the wall repulsion force. The further study of the motion of a cylinder of axisymmetric shape in simple shear flow will be reported in a forthcoming paper.