A numerical and experimental investigation of rolling and sliding motion of rotating spheres inside a cylinder

Abstract Rotating spheres in cylindrical channels roll or slide along the channel depending on the physical and geometric conditions. For a thorough investigation of the phenomenon, finite-element modelling is utilized to obtain the resistance coefficients for the motion of a sphere in a cylindrical channel, with an emphasis on near-wall motion. Extracted coefficients are compared with the data in the literature and utilized in exploring the conditions for rolling versus sliding along the channel. Sliding occurs due to the pressure build-up in the nip region between the wall and the rotating sphere in small confinement ratios, whereas rolling occurs when the shearing forces on the sphere are dominant in larger ratios. According to numerical results, a flow reversal takes place in the nip region and reduces the shear as well. Rolling versus sliding is demonstrated in experiments by using magnetic spherical particles, which are rotated by means of an external magnetic field inside cylindrical channels filled with viscous fluids. Faster axial velocities are observed in narrow channels while sliding than in wider channels while rolling for the same rotation rate of the sphere. Experiment observations are compared with the velocities evaluated from the resistance coefficients, showing that the distance between the sphere and the wall, which is dominated by roughness, plays an important role in the velocity of the sphere.

fields, such as sedimentation (Jayaweera, Mason & Slack 1964;Batchelor 1972;Bungay & Brenner 1973;Herron, Davis & Bretherton 1975;Arigo et al. 1995;Zhang & Muller 2018), lubrication (O'Neill & Stewartson 1967;Barnocky & Davis 1989;Higdon & Muldowney 1995;Gopinath, Chen & Koch 1997;Marston, Yong & Thoroddsen 2010), microfluidics (Bhagat, Kuntaegowdanahalli & Papautsky 2009;Koklu, Sabuncu & Beskok 2010) and micro/nanorobotics (Avron, Kenneth & Oaknin 2005;Golestanian & Ajdari 2008;Silverberg et al. 2020). Limited analytical solutions are available under simplifying assumptions, especially at low Reynolds numbers. Basset (1888), Boussinesq (1903) and Oseen (1927) studied the motion of a sphere settling under the gravity force in a quiescent fluid. In such a fluid, disturbance to the flow occurs solely due to the settling motion of the sphere, which is of low Reynolds number, and this allows the deduction of the resulting fluid force on the sphere using the Stokes equations (Maxey & Riley 1983). Tchen (1947) included the effects of unsteady flows in his PhD thesis, which prompted an immense number of studies suggesting corrections to his equations. Among the notable corrections, Corrsin and Lumley's (1956) remark on the contribution of the pressure gradient on the net force acting on the particle, and Buevich's (1970) correction on the term suggested by Corrsin & Lumley (1956) should be listed as well. Soo (1975) and Gitterman & Steinberg (1980), on the other hand, offered their own solutions. Maxey & Riley (1983) gave the equation the form that is widely used to this day, with corrections by Auton, Hunt & Prud'Homme (1988) and Maxey himself. Brenner & Happel (1958) investigated the frictional drag on a confined sphere subjected to a Poiseuille flow using the method of reflections. They concluded that the drag is minimized at an optimal distance away from the cylindrical channel boundaries. However, their results are valid in asymptotic cases where the distance between the sphere and the channel wall is much larger than the sphere radius. Later, Brenner & Sonshine (1964) calculated the torque required to maintain steady rotation of a sphere inside a cylindrical conduit. Their data show that the resistance to rotation increases logarithmically as the confinement increases. Bungay & Brenner (1973) studied the motion of spherical particles in a tightly fitting cylindrical conduit and proposed an improvement on the existing lubrication theories, which is still widely used in the cases where the sphere and the channel wall are in close proximity.
The limitations of asymptotic models have been overcome only recently. Higdon & Muldowney (1995) used a spectral boundary element method to obtain translational resistance coefficients of torque-free spheres moving inside cylindrical conduits. They presented tabulated results for a range of confinement ratios at any distance from the channel wall. For the cases when the sphere is too close to the channel wall, they employed the lubrication theory. As zero torque conditions are applied, rotational resistance coefficients and coupling coefficients are not reported. The most recent and comprehensive study on the topic is presented by Bhattacharya, Mishra & Bhattacharya (2010) in which the authors presented a semi-analytical method called the basis transformed spectral method (BTSM) to calculate translational, rotational and coupling resistance coefficients for spheres at a wide range of radial positions and for various confinement ratios. In this method, reflection relations for separable solutions of the flow field, represented by a basis function expansion governed by the Stokes equations, at the surfaces of the spherical particle and the cylindrical channel wall are utilized. The study is the first to report the exact coupling coefficients between the rotation and translation of a spherical particle inside a cylindrical channel. The authors explain the transition from rolling to sliding through the change in the sign of the coupling coefficients, which come out from the opposing effects of the pressure and shear forces on the particle.
One of the earliest reports on rolling and sliding is by Goldman, Cox & Brenner (1967), where the authors deduce that a sphere should slip as it rolls near a boundary.
The phenomenon is demonstrated by Liu et al. (1993) experimentally. The authors found that when a sphere is dropped near a planar wall, depending on the nature of the fluid used (Newtonian versus non-Newtonian) and the angle of inclination of the wall, the sphere might perform normal or anomalous rolling. Anomalous rolling is defined as the sphere rolling in a direction against its direction of rolling in the dry rolling case, for example, rolling upwards as it falls through a vertical tube. It is similar to what we call sliding in this study, but the translation is induced not by the rotation of the sphere but by gravity. The sphere exhibits anomalous rolling in both Newtonian and non-Newtonian fluids, and it shies away from the wall when the wall is vertical. The researchers observed that the sphere transitions to normal rolling in Newtonian fluids once the inclination of the planar wall is beyond a critical angle. However, anomalous rolling persists in non-Newtonian fluids regardless of the inclination angle. Similar behaviour patterns are observed for spheres falling down cylindrical tubes as well (Humphrey & Murata 1992), and more studies reporting the behaviour of spherical particles approaching a boundary or falling near a boundary (Dreyfus et al. 2005;Takagi et al. 2014;Djellouli et al. 2017), and studies on the collective behaviour of multiple particles (Brenner 1961;Bico et al. 2009), are also available in the literature.
Anomalous rolling is attributed to shearing at the large space between the sphere and the wall (Humphrey & Murata 1992). Bhattacharya et al. (2010) highlight the effect of lubrication as the sphere gets closer to the boundaries, so much so that the coupling resistance changes its sign and the sphere exhibits rolling instead of sliding. One important consideration at close proximity becomes the surface morphology of the sphere as roughness elements start to affect the distance from the sphere to the boundaries. Smart, Beimfohr & Leighton (1993) investigate rough spheres rolling down planes and find that the change in the distance from the sphere to the plane changes the coefficient of friction of the sphere, which manifests itself as fluctuations in the sphere velocity. When a roughness element makes contact with the plane, the contact may initiate normal motion relative to the plane that would decrease the rotation and increase the slip. The authors also provide a theoretical model that is in quantitative agreement with their experimental results with rough spheres. A more rigorous model is developed by Galvin, Zhao & Davis (2001) for a sphere rolling down a tilted plane, in which they define two roughness scales for the sphere. Large roughness elements temporarily lift the sphere, and as it rotates, it moves away from the plane, and this causes the sphere to lose contact with the plane. The sphere is then pulled down by gravity to be lifted off from the surface once again, with an upcoming roughness element. The authors report that the sphere is in contact with the plane for a longer time at lower inclination angles than in higher inclination angles, but the hydrodynamic resistance appears to be greater at high inclination angles. This is explained by the sphere's faster rotation leading to more frequent contact of the large roughness elements with the plane. Based on the model of Galvin et al. (2001), Zhao, Galvin & Davis (2002 study the problem of a smooth sphere rolling down a rough plane with two different roughness scales again. Upon contact with a large bump, the translational velocity of the sphere decreases as the rotational velocity increases, then the sphere quickly loses contact with the bump and the translational velocity decreases further upon contact with the small bumps. The dimensionless translational velocity is much greater than the dimensionless rotational velocity, indicating that the sphere slips all this time. Upon contact with the second large bump, the velocities coincide and slipping stops. These observations are quite critical as the gap size between the sphere and the channel wall determines the mode of motion of the sphere. A comprehensive investigation of the motion of a rotating sphere in close proximity to the boundaries is necessary. Our study aims to understand the effects of geometric parameters and to elucidate the rolling and sliding of spheres in cylindrical channels. In that regard, we study the effects of the distance of the sphere from the channel boundaries, and the confinement ratio, which is the ratio of the radii of the channel and the sphere, numerically and experimentally. The transition between rolling and sliding is demonstrated with respect to the confinement ratio. To this end, we first introduce a finite-element method-based (FEM) numerical model to obtain the complete set of resistance coefficients for a sphere inside a cylindrical channel especially for the case where the sphere is very close to the channel wall. The resistance coefficients are systemically derived by evaluating the forces and torques on the sphere for given swimming velocities. Unlike in Higdon & Muldowney (1995), coupling and rotational resistances are included here alongside the translational resistance coefficients. Moreover, we verify that the FEM model is more efficient and accurate especially for the case when the sphere is very close to the channel wall compared to the semi-analytical model presented by Bhattacharya et al. (2010). Furthermore, we demonstrate the rolling and sliding of rotating spheres in cylindrical channels experimentally for the first time in the literature. In our experiments, magnetized spheres with considerable roughness are placed in viscous fluid-filled cylindrical channels and rotated with the help of a rotating magnetic field. Both the rolling and sliding cases are reported and characterized. Finally, the experiment results are confirmed with the velocities obtained from the resistance coefficients, showing that the sphere is rotating in close proximity to the channel boundaries but the exact proximity cannot be determined due to the roughness of the spheres used in the experiments and limitations in image processing capabilities.

Resistance coefficients
Consider a sphere with diameter D s rotating inside a viscous fluid-filled cylindrical channel with diameter D ch , as shown in figure 1. Inertial effects in low Reynolds number motion are generally negligible, which is why the forces and torques acting on the sphere are directly related to the linear and angular velocities of the sphere through a resistance matrix, R. The matrix is generally expressed in terms of its four subcomponents as In this equation, F is the viscous force acting on the sphere and τ is the torque. F tt is the translational resistance matrix, F rr is the rotational resistance matrix, and F tr and F rt are the coupling resistance matrices with F tr = F rt where the ' ' sign indicates the transpose. The linear velocity vector is U, and the angular velocity vector is ω. Two coordinate frames will be used in the text. One follows the notation in Bhattacharya et al. (2010): it is a cylindrical coordinate frame with the z axis placed alongside the long axis of the cylindrical channel. The radial direction is identified withρ, and the tangential direction is identified withβ -unit vectors in figure 1 -along with a global Cartesian frame.
The elements of the resistance matrix are obtained in the cylindrical coordinates by running a series of simulations in which one component of linear or angular velocities is set to unity and the rest are set to zero. Most of the off-diagonal entries are found to be infinitesimally small (10 5 times smaller than the parameters listed in table 1) so they are  assumed to be zero, resulting in the following explicit form for (2.1): ⎡

Inputs Outputs
(2. 2) The whole set of cases and the resulting equations from each one of the cases are listed in table 1. A total of six separate simulation runs are required to obtain all of the components in (2.2) for a given position of the sphere.

The finite-element model
Sphere motion in viscous fluids and at very small scales (Re 1) is governed by the Stokes equations. The equations are written in non-dimensional form as (2.3a,b) Here, u andp are the non-dimensional fluid velocity field and the pressure, respectively. The length scale for non-dimensionalization is the sphere radius R s , the time scale is the rotation frequency of the sphere f , and the mass scale is the fluid viscosity μ; u is non-dimensionalized with lf whilep is non-dimensionalized with μf . The COMSOL Multiphysics software package is used to solve the incompressible Stokes equations. No-slip boundary conditions are applied on the channel boundaries and the sphere surface. The sphere surface is modelled as a moving wall with a velocity profile expressed as where r is a position on the sphere surface S and r 0 is the position of the centroid of the sphere.
Taking advantage of the symmetries in the model, the computational domain is cut in half in circular and rectangular cross-sections of the cylinder through the sphere for computational efficiency, as shown in figure 2. A slip boundary condition is applied at the cut planes as a symmetry condition. The complete three-dimensional (3-D) geometry is required only to obtain the rotational resistance of the sphere in the radial direction, F rr ρρ , with the same meshing parameters used in the cut geometries.
The P2+P2 discretization of the fluid and MUMPS direct solver are employed for the simulations. Tetrahedral elements are used to mesh the fluid domain, and triangular surface mesh is applied to the sphere surface, with the same meshing applied to symmetric pairs of the sphere faces to improve the accuracy of the solutions. The coefficients in (2.2) are obtained from the scenarios listed in table 1 for a wide range of non-dimensional radial positions, defined asρ where ρ s is the dimensional radial position of the sphere and R ch is the radius of the cylindrical channel. The minimum distance from the sphere surface to the channel boundaries is identified by δ and is non-dimensionalized as We first demonstrate the mesh convergence for the configurationρ s = 0.8 and R ch /R s = 1.6. The narrowest channel size is considered here so as to demonstrate the convergence with respect to the densest meshing possible. The converged configuration will be used to obtain the resistance coefficients for 0 ≤ρ s < 0.9. The meshing strategy is slightly altered forρ s ≥ 0.9 as the convergence in this range is much more demanding.
The mesh convergence study over the domain element size shows that the system is relatively insensitive to this parameter (maximum element size ranging from 0.1R s to R s ), with a relative error of less than 1 % even at the coarsest meshing configuration (not shown). The meshing on the spherical surface appears to be more critical for convergence, with the results demonstrated in figure 3(a) for several key resistance coefficients. The relative percentile error, e, is defined as where {F tt zz , F rr zz , G} max indicate the values obtained at the maximum degrees-of-freedom that corresponds to the smallest mesh element size. The converged configuration results in around 4 × 10 5 degrees-of-freedom in cut geometries and takes up to 100 GB of random access memory (RAM) usage.
When the sphere is very close to the channel wall, special care must be taken to ensure converging results. In this work, we take 0.9 ≤ρ s ≤ 0.99 to be the close proximity range, corresponding to 0.006 ≤δ ≤ 0.2. The minimum element dimension in the mesh is adjusted to accommodate properly the small gap between the sphere and the channel. As the sphere gets closer to the channel boundaries, a large pressure gradient builds up across the nip region between the sphere and the channel wall when the sphere rotates in the β direction (azimuthal). An accurate solution of the pressure gradient is critical to obtain resistance coefficients with high accuracy. Hence the meshing density in the fluid surrounding the sphere is increased to match the density on the sphere (the regions coloured orange in figure 2). A mesh convergence study is carried out forρ s = 0.99 and R ch /R s = 1.6, corresponding to the tightest configuration in the scope of this work, and the convergence of F tt zz , F rr zz and G are displayed in figure 3(b). The results indicate that the relative error falls below 1 % as the degrees-of-freedom approach 3million. The maximum element size around the sphere is δ/2, and the minimum element size is δ/40 in the converged configuration.

Experiments
In experiments, radially magnetized nickel-plated sintered neodymium (NdFeB) spheres (SM Magnetics, Pelham, AL, USA) of diameters 1 mm and 1.9 mm are placed inside glass channels of diameters 1.6 mm, 3 mm and 5.7 mm. The roughness of the spheres, which is critical in their motility, is investigated using a Nanofocus 3-D surface metrology system. The results of the measurements are provided in Appendix A. We present a close-up image of the sphere with D s = 1 mm in figure 4(a), which shows the roughness of the spherical surface. Since the spheres are made of magnetic ceramics with nickel coating, it is very difficult to clean them free from pieces of chipped coating and other magnetic debris that accumulates at the surface.
The channels are filled with silicone oil mixtures with μ = 0.5 Pa s and μ = 1 Pa s. Removal of excessive air inside the liquid is critical in order to obtain matching results with the simulations. The tubes are placed into a vacuum chamber (0 PSIA, measured with Omega DPG5600B-30A PSIA) for degassing before experimentation. After the degassing procedure, the tubes are sealed tightly to prevent air from leaking back into the liquid.
The sealed tubes are placed horizontally inside the experiment set-up consisting of three orthogonal Helmholtz coil pairs to induce sphere rotation. The coil system is controlled with custom LabVIEW software via Maxon drivers connected to the computer controlling the experiments. The experiment set-up is drawn in figure 4(b) and pictured in figure 4(c).
The spheres are rotated with a rotating magnetic field to observe sliding and rolling trajectories. Two of the coil pairs, placed along the y and z directions, are excited with sinusoidal out-of-phase currents with amplitude I 0 to create a magnetic field rotating about the x direction in the global frame. The magnetic field applied by each coil is measured using Phidgets 1108 Magnetic Field Sensors to assure equal magnetic field strength in both directions. The magnetic spheres are actuated at different magnetic rotation frequencies (f ) ranging between 0.1 Hz and 20 Hz, and the trajectories are recorded using a digital microscope from above (refer to figure 4c). Gravity, denoted g, is acting in the −y direction.
As the rotating magnetic field is applied to the magnetic sphere, magnetic torque tends to align the sphere's magnetic dipole moment, m, with the direction of the applied magnetic field, B, so that the sphere rotates around the β direction in the cylindrical frame as the magnetic field rotates in the x direction in the global frame. The relationship to evaluate the induced magnetic torque, τ m , is given by the equation ( 2.8) As implied by the cross-product, actual magnetic torque acting on the sphere at any given instance depends on the sine of the angle, ζ , between the magnetic dipole moment vector of the sphere and the magnetic field vector, when the two vectors are co-planar as shown in figure 4(d). The angle depends on the viscous torque, which balances the magnetic torque assuming that the sphere rotates at the same rate as the rotating magnetic field. A schematic of the sphere motion due to the rotating magnetic field inside the channel is given in figure 4(d). The spheres are not able to rotate indefinitely faster as the magnetic torque rotating the sphere cannot overcome the viscous resistance beyond a certain f . The sphere rotation stutters at larger f , which is called step-out in the literature (Zhang et al. 2009).
A rotating sphere near the channel wall translates along the z axis of the cylindrical channel. The translation velocity, u z , is extracted from the experiment recordings using the image processing code reported in our previous work, which utilizes MATLAB's Image Processing Toolbox functions (Caldag, Acemoglu & Yesilyurt 2017).

Validation of the finite-element model
The results from the FEM model are compared with two different datasets from the literature for validation purposes.ρ s is varied from 0 to 0.99 in the FEM simulations. Selected R ch values are 1.6, 2 and 3, while R s is fixed at 1. For the sake of brevity, the verification results will be presented only for R ch = 2 in this section. The complete list of the resistance coefficients and comparisons with the data from the literature for R ch = 1.6 and R ch = 3 are provided in Appendices B and C.
The resistance coefficients are normalized as in Bhattacharya et al. (2010): (3.4) Figures 5 and 6 show the parameters obtained via the FEM model and the results from two studies in the literature. Translational and rotational resistance coefficients match quite well with the reported data for 0 <ρ s < 0.9. When we compared our results with those in Bhattacharya et al. (2010), we observed some discrepancies, especially asρ s → 1. The authors kindly provided their code for the re-evaluation of the resistance coefficients with higher accuracy at Λ max = 16, μ max = 10, l max = 12 and δ λ = 0.02 (refer to Bhattacharya et al. (2010) for the definitions of these parameters). Updated resistance coefficient values are shown in red in figures 5 and 6, and agree much better with the FEM results than the values reported in Bhattacharya et al. (2010), especially for near-wall values (ρ s > 0.9).
For the coupling resistancesḠ andḠ , shown in figure Bhattacharya et al. (2010) and the re-evaluated coefficients stem from the selection of the model parameters that are critical for convergence. The output of the Bhattacharya et al. (2010) model is a grand mobility matrix whose dimensions ideally go to infinity. The matrix is truncated to a certain dimension, denoted by q = 3l max (l max + 2), for matrix inversion, which is a required step in obtaining the friction matrix that gives out the resistance coefficients. Each term in the mobility matrix is also an approximation in itself as each term includes a truncation of an infinite summation and an infinite integration. The authors reported a convergence study over multiple model parameters, including Λ max , δ λ , μ max and l max for translational and rotational resistance coefficients, which tend to converge fast forρ s = 0.5 andρ s = 0.9, with relatively low computational requirements. The convergence of the coupling coefficients is omitted in that study; we find that they do not converge as fast, especially asρ s → 1. The results of a new convergence study with the Bhattacharya et al. (2010) code over l max forρ s = 0.99 (provided in the figure in Appendix D) show that the coupling coefficients barely converge for the largest l max tested. Increasing l max further had convergence issues in the model. The improvement in the matching of the results with the re-evaluated data from Bhattacharya et al. (2010) permits confidence in the high-resolution FEM results. One more point to note here is that it takes up to 24 hours for the Bhattacharya et al. (2010) code to finish running for l max = 18, whereas our FEM model with the densest meshing takes up to 10 hours on the same workstation (a minimum of six separate simulations are required, with around 1.5 hours of runtime for each) and provides high accuracy for most of the parameters at a much lower computational cost. Furthermore, one has to carry out a multi-parameter convergence study for the Bhattacharya et al. (2010) model by covering other parameters listed above, which would increase the overall computational cost even further. The FEM model is very useful for single-particle systems but it may become costly for solving systems involving multiple particles. In that case, BTSM can be utilized for a global solution, and the FEM model can be used to resolve local fields involving fewer particles. The convergence of BTSM in earlier work appears to be incomplete, especially in terms of the coupling resistances and at close proximities. A multi-parameter convergence on BTSM is necessary to fully benefit from this approach for spheres in close proximity to the channel boundaries.
3.2. Rolling and sliding Rolling and sliding are tied to the coupling and translational resistance coefficients. From (2.2), one can write where u β and u z are the rolling/sliding velocities in the respective directions. Utilizing the normalization in (3.3) and (3.4), we define a non-dimensional velocityū: (3.8) Figure 7 depicts the non-dimensional velocitiesū β andū z for all R ch /R s values tested forρ s > 0.9 with respect toδ. The rotation of the sphere around the z axis gives rise to the sliding motion in the β direction, as shown in figure 7(a). As the sphere rotates, a pressure gradient (with maximump = p H and minimump = p L ) develops between the fore and aft of the sphere that induces sliding motion. Previous experiments had shown the sliding behaviour along the channel boundaries as the sphere rotates around the channel's long axis (z axis), resulting in circular trajectories around the long axis of the channel without any translation in the z direction (Demir 2018). At sufficiently high rotation rates, the radius of the circular trajectory decreases and the sphere settles at the centre of the channel radially (Demir 2018). Here,ū β is positive for the cases depicted in the figure; however, it starts to decrease asδ → 0, especially at high values of the curvature (shown in figure 7b). AlthoughḠ keeps increasing asδ → 0,F tt ββ exhibits a logarithmic increase (also predicted by Higdon & Muldowney 1995) that leads to an overall decrease in the ratio.
Rotation around the β axis results in axial sliding or rolling motion along the channel, as shown in figure 7(c). As shown in figure 7(d),ū z values are mostly negative, indicating that the sphere slides. Also worth noting is the fact thatū z varies logarithmically withδ. Asδ → 0,ū z decreases and changes sign atδ = 0.02 for the largest curvature ratio R ch /R s = 3, which indicates that the force due to the pressure difference in the nip region is not large enough to overcome the shear force for sliding.
Rolling and sliding phenomena are associated with the relative dominance of lubrication (F v ) and pressure forces (F p ) acting on the sphere in the literature. Bhattacharya et al. (2010) report that the rapid increase of shear forces compared to the increase in pressure forces at close proximity (ρ → 1) leads to decreases in the magnitudes of these coefficients. However, the underlying mechanisms forḠ andḠ are not exactly the same.Ḡ is the term that relates the axial torque with tangential motion or the tangential force with axial rotation, whileḠ relates the axial force with tangential rotation or the tangential torque with axial motion. When the ratio of the magnitudes of the pressure and shearing forces in coupling resistances is plotted as in figure 8(a) forḠ, and figure 8(b) forḠ , it is observed that the pressure-induced forces remain dominant inḠ (|F p β /F v β | > 1, where the subscripts denote the direction) even asδ → 0, meaning that the sphere tends to slide along the boundary. On the other hand, forḠ , dominance of the pressure contribution lessens asδ → 0. As the sphere moves closer to the boundaries, it tends to roll along the channel boundary as dictated by the dominating shearing effects. The sphere slides as it  rotates around the z axis in most of the configurations, because the pressure-induced forces remain dominant.
The distributions of pressure and shear on the sphere help in understanding the dynamics of rolling and sliding. As the sphere rotates around the β axis (azimuthal direction), a pressure gradient develops between the fore and aft of the sphere (shown in figure 9a). Regions with large pressure amplitude are quite small but significant. The pressure profile along the bottom half of the sphere is drawn in figure 9(b). Note the dramatic change from positive to negative pressure values through the nip region between the sphere and the cylindrical channel. Negative pressure levels may be deemed an indicator for cavitation, but it should be noted that the zero pressure level is with respect to a faraway point inside the channel, meaning that the absolute values must be calculated with respect to the reference pressure. A simple calculation, provided in Appendix E, shows that there should be no cavitation. The difference between the maximum and minimum pressures, p =p H −p L , increases monotonically with respect toδ. As plotted in figure 9(c), p due to rotation of the sphere follows a monotonic trend withδ for all confinement ratios asδ → 0, and the slope in the logarithmic scale is −0.5, which is consistent with the lubrication theory (Higdon & Muldowney 1995).
Looking at the shear stress distribution along the bottom arc (the non-dimensional shear stress is denoted asτ ), plotted in figure 9(d), a striking dip is observed right at the bottom of the sphere. This could be explained by the flow reversal in the back and front of the sphere, as shown by the streamlines in figure 9(a). Shear rates go through a maximum at the edges of the nip region due to the flow reversal. The dip inτ disappears completely when the sphere is sufficiently far from the wall, atρ s = 0.2. The maximum shearτ max at lowδ (shown in figure 9e) exhibits a trend similar to that of p, albeit that the magnitude is an order of magnitude lower for a givenδ. Also note that the slope is smaller than the value of −0.5 observed for the pressure gradient; it comes out as −0.833.

Experiment results
This subsection reports the velocities of the magnetically rotated spheres from our experiments. The spheres are observed to be 'rolling', i.e. translating in the positive z direction as they are rotated counter-clockwise about the β axis when R ch /R s = 3, as sketched in figure 9(a). Translation in the opposite of the rolling direction, which is referred to as 'sliding', occurs as a response to counter-clockwise rotation about the β axis when R ch /R s = 1.6, also shown in figure 9(a).
Values of u z for the experiments with R ch /R s = 3 are plotted against the rotation frequency, f , in figure 10(a,b). Fluids with two different viscosities, μ = 0.5 Pa s and μ = 1 Pa s, and spheres with diameters D s = 1 mm and D s = 1.9 mm, are tested. When D s = 1 mm, D ch = 3 mm and μ = 1 Pa s, u z increases linearly with increasing f up to 8 Hz (figure 10a). As the rotation frequency is increased beyond this value, the sphere fails to rotate synchronously with the rotating magnetic field, thus the sphere velocity decreases with increasing f up to 20 Hz.
Step-out also causes large deviations in u z , as shown by the error bars in the plots, which denote the standard deviation values. The fluctuations in u z outside the step-out regime can be explained easily by the roughness of the spheres, which leads to a time-varying δ that alters the resistance coefficients (Smart et al. 1993). The step-out frequency and translation velocities are higher overall for the D s = 1.9 mm and D ch = 5.7 mm configuration (10b), owing to the stronger magnetization and increased weight of the sphere that improves the traction. u z values for μ = 1 Pa s and μ = 0.5 Pa s are more or less similar up to f = 3 Hz at both geometric configurations (shown in the insets), but they deviate at larger f . Values of u z for the experiment configurations with sliding (R ch /R s = 1.6) are displayed in figure 10(c,d). The red lines passing through u z = 0 highlight the transition from rolling (i.e. u z > 0) to sliding (i.e. u z < 0). The transition to sliding occurs at very low f (around 0.5 Hz) for μ = 0.5 Pa s, as shown in the insets. At low frequencies, having a considerable roughness at the surface, spheres establish contact with the wall and roll slowly due to traction. Note that such a motion occurs in only a very small number of experiments. Contact of the roughness elements with the channel boundary could induce a lift that would deter the traction and cause sliding, but the lift appears to be limited as the sphere is able to maintain rolling motion, whereas at higher rotation frequencies the traction is lost and the pressure difference leads to a sliding motion. The variations in u z with respect to f are mostly linear at the sliding region. Note the overall increases in the magnitudes of the maximum velocities attained before step-out in comparison to u z observed in rolling spheres. The increase is particularly notable as the viscous effects at narrower channels are expected to be more restraining against motion, as implied by the resistance coefficients reported in § 3.2. The enhanced swimming speeds are due to the large p that contributes to the sliding of the sphere as opposed to rolling, where p hinders the sphere motion.
The experiment results can also be compared with the velocities evaluated from the resistance coefficients with (3.6). This simple calculation means that several types of forces will be neglected. Unsteady forces, such as the history force and added mass forces, are known to play an important role in the swimming of micro-organisms. Jakobsen (2001) reports that Balonion comatum, a ciliate plankton, increases its velocity fivefold in a time period shorter than the time needed to advance the organism more than one body length. Such motions create unsteady disturbances in the flow field that can affect the velocity and trajectory of the swimmers even after the motion causing the disturbance ceases. However, in the scope of the study reported here, these unsteady forces can be neglected, as the density of the particle used in this study is not comparable to the density of the fluid used in the experiments (Van Aartrijk & Clercx 2010). Wang & Ardekani (2012) model a spherical unsteady swimmer and show that the Boussinesq-Basset history term and the added mass term can be neglected when the product of Strouhal (Sl) and Reynolds numbers is smaller compared to unity: Here, m s is the mass of the swimmer and m f is the mass of the fluid displaced by the swimmer. In this study, the highest Sl Re that occurs throughout the experiments is 0.1886, which is achieved when μ = 1 Pa s, f = 20 Hz and R s = 0.95 mm. However, as this rotation frequency is above the step-out frequency, above which the sphere rotation loses its synchronization with the rotating magnetic field, the rotation rate of the sphere does not reach 20 Hz at all. Therefore, the actual Sl Re value for this configuration is below 0.1886. Thus the effects of history and added mass can be discarded safely. A critical omission in this approach is the roughness of the sphere, which would bring about a time-varyingδ that would lead to time-varying resistance coefficients as reported in Galvin et al. (2001). The roughness causes non-continuous traction of the sphere on the surface of the channel, which is highly critical for rolling motion. Note δ in the experiments cannot be determined as accurately as needed due to the limitations in our image processing capabilities. Withδ being unknown, the calculation in (3.6) is carried out for multipleδ.
Another consideration would be the effect of the lift force on the sphere. There exist an extensive number of studies investigating the lift force on spheres at low Reynolds numbers (Saffman 1965;Cox & Brenner 1968;Ho & Leal 1974;Vasseur & Cox 1976;Cox & Hsu 1977;Drew 1988;McLaughlin 1993;Cherukat & McLaughlin 1995), but these studies either do not fit into our configuration or cannot be implemented due to their nonlinear nature. Other models for spheres swimming in bulk or at higher Reynolds numbers point to a strong correlation between an increase in lift force with increasing angular velocity. Lift force is reported to affect δ, and a relevant example would be the study by Bhattacharya, Gurung & Navardi (2013), where the authors report equilibrium radial positions (where the inertial lift is balanced by the rest of the hydrodynamic forces) for the spheres inside cylinders with respect to the curvature ratio. Nonetheless, the lift induced by the roughness elements on the spheres appears to be more significant as the inertial lift force should be very low considering the Reynolds number of the system. Experiment results and calculated u z values for multipleδ values are plotted in figure 11 with respect to f up to step-out frequencies for each case. Resistance coefficients for δ = 0.002 andδ = 0.0006 (whereρ s > 0.99) are evaluated with the FEM model (only the necessary terms), while the coefficients forδ = 0.01 andδ = 0.003 are interpolated using piecewise cubic spline interpolation. The results for the D s = 1 mm and D ch = 3 mm configuration (figure 11a) show thatδ is between 0.01 and 0.002 in experiments, whereas in D s = 1.9 mm and D ch = 5.7 mm configuration (figure 11b) the experiment values fall between the u z values forδ = 0.02 andδ = 0.01, which indicates that the gap for the larger sphere in the larger channel is higher. Either way, the sphere appears to be very close to the boundaries and rolls along with the help of the traction. For both of the sliding configurations (figure 11c,d), u z in experiments are between the results forδ = 0.006 and δ = 0.003. The sphere is still close enough to get traction, but p is so great that it results in sliding motion. Galvin et al. (2001) state that the average δ between the sphere and the boundary should be between the smaller and larger roughness sizes, so these fits give an estimate for the average roughness element sizes, giving 1-5 μm for the sphere with D s = 1 mm, and 9.5-19 μm for the sphere with D s = 1.9 mm, in dimensional terms. The scale of the values with respect to sphere radii is in agreement with what Smart et al. (1993) report. Also worth noting is that the values of u z fit to differentδ values at low and high f values, meaning that δ is higher at larger f . This is consistent with the Galvin et al. (2001) report where the authors state that a faster rotating sphere would have its large roughness elements making contact with the boundaries more frequently, which results in a larger δ in average. Finally, note that the estimated sizes of average roughness elements are of the order of the gap sizes reported for cavitation in Appendix E, therefore it is possible that cavitation might have occurred but we do not expect it to be significant compared to the effects caused by the roughness elements.
Sliding and rolling velocities of the sphere are higher in magnitude in the higher viscosity fluid than the ones in the lower viscosity fluid. This indicates that the sphere's distance to the wall, δ, is larger during sliding at the more viscous fluid, but also that δ is smaller during rolling at the more viscous fluid. It is not clear why such a dichotomy occurs in the results. Furthermore, the competing effects, pressure and viscous forces, scale linearly with μ at the Stokes regime, thus there should be no difference with respect to μ in sphere motion. For sliding, the roughness elements appear to deter the traction. As reported by Galvin et al. (2001), when a large roughness element on the sphere makes contact with the boundary, the sphere is temporarily lifted off from the surface. In a more viscous fluid, it would take more time for the gravitational forces to pull the sphere down to lower δ. That is why u z is higher in magnitude with the more viscous fluid and also why u z data for large f from the experiments fit better to the velocity profile for the larger δ. For rolling, the roughness elements appear to improve the traction of the sphere with increasing viscosity, thereby decreasing δ. A second possibility would be that in the high viscosity case, pressure may drop to a level that causes evaporation in the silicon oil, limiting the minimum pressure to vaporization pressure, which is very small for silicon oil. On the other hand, the high pressure in the nip region is not bounded, thus pressure forces may not cancel and a net pressure force may lead to the levitation of the sphere. However, the analysis in Appendix E indicates that no cavitation should occur. In our numerical analysis, we could not confirm the presence of a lift force on the rotating sphere with or without the inertial forces. Therefore, the trend appears to be related to the roughness elements on the sphere.

Conclusion
Motion of rotating spherical particles in cylindrical conduits is investigated here numerically and experimentally. For the first time, a unique experiment set-up is used to demonstrate the rolling and sliding behaviour of magnetically rotated spheres inside cylindrical channels filled with a viscous fluid. Elements of the resistance matrix are calculated using a validated finite-element model with a lower computational cost than analytical models, especially for the case when the motion of a single sphere very close to the channel wall is considered. The resistance coefficients are reported for different confinement ratios and a range of radial positions, with a special focus on the motion in close proximity to the channel boundaries.
Resistance and coupling coefficients in the resistance matrix are used to study the rolling and sliding phenomena. A near-wall rotating sphere slides along the channel for small confinement ratios due to the very high pressure gradient in the nip region between the channel and the sphere, whereas the shearing forces are dominant in rolling observed for larger confinement ratios. Another interesting finding is that a flow reversal causes a dip in the shear stress in the nip region.
In the experiments, radially magnetized spherical magnets are placed in viscous fluid-filled horizontal cylindrical channels. Due to their weight, spheres are observed to be resting on the channel wall. When the spheres are rotated in the azimuthal direction with the help of a rotating magnetic field, sliding motion is observed along the channel in the opposite direction to the rolling direction at high confinement ratios owing to the pressure buildup near the nip region between the sphere and the channel. It is shown that the pressure gradient has similar magnitudes in different confinement ratios for a given non-dimensional gap size between the sphere and the channel. This would imply that sliding should be possible at all confinement ratios, but the differences in resistance coefficients at different confinement ratios prevent this from happening.
In the experiments, sliding is observed for the confinement ratio 1.6, whereas rolling is observed for the confinement ratio 3 for both viscosity values. Sliding velocities are larger in magnitude than the rolling velocities for a given magnetic rotation frequency as the sphere has to overcome the pressure gradient in the rolling motion whereas the sphere is pushed by the pressure gradient in the sliding regime. Translation velocities evaluated from the resistance coefficients for the experiment configurations indicate that the sphere is very close to the channel boundaries in both swimming modes. The proximity to the wall cannot be predicted accurately due to both image processing limitations and the roughness of the spherical surface. In small confinement ratios, rolling can take place at very low rotation rates, but as the rotation rate increases, a transition from rolling to sliding occurs. Surface roughness of the sphere causes a lift-off from the channel boundaries that affects the average gap size between the channel wall and the sphere, leading to fluctuations in the axial velocity and loss of traction that is critical for rolling, and also larger axial velocities at more viscous fluids. Faster rotation rates of the sphere amplify the effects of the surface roughness, increasing the gap size as evidenced by the evaluated velocities from the resistance coefficients. Overall, the findings of this study are expected to improve understanding on the motion of spherical particles in cylindrical channels, which is of interest from different perspectives in fluid mechanics. For a specific example, results could prove useful to study magnetic spherical particles as micro robots in cylindrical conduits in microfluidic applications. ratio indicates that the surface has significant dips or peaks along the measured profile, and for this sphere the ratio comes out as 7.37. The sphere with D s = 1.9 mm exhibits lower roughness, with an Rt/Ra ratio of 3.68. A ratio of 1 would indicate that there are no significant dips and peaks that go further beyond the average roughness on the surface.

Appendix B. Complete list of resistance coefficients
The complete list of the resistance coefficients in normalized form is provided here.    may occur at closer distances under those conditions (μ = 1 Pa s, f = 20), atδ < 0.0027, or δ < 2.7 μm. For a smaller frequency, f = 10 Hz, for which we assume a synchronous rotation of the sphere with the magnetic field, the gap thickness for which cavitation is expected would be smaller, i.e.δ < 0.0019, or δ < 1.9 μm.