1. Introduction
The problem of a body immersed in a fluid, and which both translates and rotates while in contact with, or near, a wall is of practical interest in understanding and modelling sediment transport in rivers and coastlines, and improving the efficiency of various industrial processes involving particle-laden flows. One key aspect of this problem is to determine the forces and moments applied to the body, and hence, predict its motion along the boundary. Generally, the body is assumed to be offset from the wall by an effective hydrodynamic gap introduced by either surface roughness (Smart et al. Reference Smart, Beimfohr and Leighton1993; Galvin et al. Reference Galvin, Zhao and Davis2001; Thompson et al. Reference Thompson, Leweke and Hourigan2021; Houdroge et al. Reference Houdroge, Zhao, Leweke, Hourigan, Terrington and Thompson2023; Nanayakkara et al. Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024b ), cavitation (Prokunin Reference Prokunin2003; Ashmore et al. Reference Ashmore, Del Pino and Mullin2005), compressibility (Terrington et al. Reference Terrington, Thompson, Hourigan, Lei, Thornber and Armfield2022) or elastohydrodynamic effects (Rallabandi et al. Reference Rallabandi, Saintyves, Jules, Salez, Schönecker, Mahadevan and Stone2017; Bertin et al. Reference Bertin, Amarouchene, Raphael and Salez2022). While the introduction of a small hydrodynamic gap does not significantly affect the dominant flow structures, such as wake formation and vortex shedding (Stewart et al. Reference Stewart, Thompson, Leweke and Hourigan2010b ), it has a profound effect on the drag and moment coefficients, due to the large pressures arising in the narrow lubrication film between the body and the wall.
Analytical predictions of the forces and moments can be obtained for elementary geometries, such as spheres or infinite cylinders, translating and rotating near plane walls within the Stokes-flow regime. Expressions for the force and moment applied to an infinite circular cylinder in the Stokes regime were obtained by Jeffery (Reference Jeffery1922), Wakiya (Reference Wakiya1975) and Jeffrey & Onishi (Reference Jeffrey and Onishi1981), while O’Neill (Reference O’Neill1964) and Dean & O’Neill (Reference Dean and O’Neill1963) compute the corresponding force and moment applied to a sphere, also in the Stokes regime. Asymptotic expressions for the forces and moments, valid in the limit of a small gap, can also be obtained using the method of matched asymptotic expansions. Under this approach, an inner solution, valid in the thin interstice between the body and the wall, is obtained using lubrication theory. This is combined with an outer solution, which is independent of the gap between the body and the wall, to obtain the total force and moment. This approach was used by Goldman et al. (Reference Goldman, Cox and Brenner1967), O’Neill & Stewartson (Reference O’Neill and Stewartson1967) and Cooley & O’Neill (Reference Cooley and O’Neill1968) for spheres in the Stokes regime, and by Merlen & Frankiewicz (Reference Merlen and Frankiewicz2011) for an infinite circular cylinder in the Stokes-flow regime.
For inertial (finite-Reynolds-number) flows, however, analytical expressions for the forces and moments are generally not available, and numerical simulations are used to obtain predictions of the force and moment coefficients. For example, Stewart et al. (Reference Stewart, Hourigan, Thompson and Leweke2006, Reference Stewart, Thompson, Leweke and Hourigan2010b ), Rao et al. (Reference Rao, Stewart, Thompson, Leweke and Hourigan2011) and Houdroge et al. (Reference Houdroge, Leweke, Hourigan and Thompson2017, Reference Houdroge, Leweke, Hourigan and Thompson2020) perform numerical simulations for cylinders rotating and translating near a plane wall at finite Reynolds number, while Zeng et al. (Reference Zeng, Najjar, Balachandar and Fischer2009), Stewart et al. (Reference Stewart, Thompson, Leweke and Hourigan2010a ) and Houdroge et al. (Reference Houdroge, Leweke, Thompson, Hourigan, Ivey, Jones and Draper2016, Reference Houdroge, Zhao, Leweke, Hourigan, Terrington and Thompson2023) perform numerical simulations for spheres rotating and translating near a plane wall. There are several drawbacks to using numerical simulations to obtain the force and moment coefficients (Terrington et al. Reference Terrington, Thompson and Hourigan2023). First, a large number of extremely small elements must be used to resolve the inner lubrication flow, resulting in numerical stiffness and increased computational cost. This increased computational effort is essentially wasted, since analytical solutions to the inner flow are easily obtained using lubrication theory. Moreover, simulations must be performed for a range of gap-to-diameter ratios, despite the fact that this only affects the inner lubrication flow, which is already known from lubrication theory.
To counteract these problems, we proposed in a previous work (Terrington et al. Reference Terrington, Thompson and Hourigan2023) a combined numerical--analytical approach, where numerical simulations are performed to obtain the outer-flow solution only, while the analytical lubrication solution is used for the inner region. For the two-dimensional (2-D) flow over a circular cylinder, we demonstrated that the total force and moment coefficients can be decomposed into a gap-dependent part and a gap-independent part. The gap-dependent part was obtained analytically, using lubrication theory, while the gap-independent part was obtained using the numerical solution to the outer flow. The present paper applies this approach to the three-dimensional (3-D) flow over a finite-span cylinder that both translates and rotates near a plane wall.
Compared with the rolling sphere and infinite cylinder cases, the flow over a finite-span cylinder near a wall has received far less attention. Saintyves et al. (Reference Saintyves, Rallabandi, Jules, Ault, Salez, Schönecker, Stone and Mahadevan2020) has provided asymptotic expressions for the gap-dependent force and moment in the limit of a large aspect ratio. Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) extends this analysis, providing expressions for the gap-dependent force and moment valid for any aspect ratio, including the limiting cases of a large aspect ratio (thin rods) and small aspect ratios (thin disks). These studies do not compute the gap-independent contribution to the force and moment from the outer-flow region. The present paper presents a methodology to compute both the gap-dependent and gap-independent contributions to the force and moment coefficients for the inertial flow over a finite-span circular cylinder.
Previous numerical studies (Pirozzoli et al. Reference Pirozzoli, Orlandi and Bernardini2012; Javadi Reference Javadi2022) have considered the flow past a finite-length circular cylinder rolling while in contact with a plane wall. However, neither of these studies consider the significant influence of the hydrodynamic gap on the force and moment coefficients. Moreover, both papers report finite drag coefficients for cylinders in contact with the wall, while lubrication theory predicts an infinite drag for cylinders in contact with the wall. Therefore, accurate numerical predictions for the drag behind finite-span circular cylinders have not yet been obtained.
The flow over finite-span cylinders in contact with a wall is also relevant to the automotive and racing industries, as a simplified model for the aerodynamics of wheels (Diasinos et al. Reference Diasinos, Barber and Doig2015). Unlike the rigid circular cylinders considered in this work, pneumatic tyres will contact the road by a contact patch, rather than a singular contact point. Diasinos et al. (Reference Diasinos, Barber and Doig2015) report a significant variation in the computed drag coefficient with the size of this contact patch. We note that they report finite drag coefficients for isolated wheels in contact with the plane surface, in contrast with the infinite force predicted for a rigid cylinder. This may be due to the effects of tyre deformation, or some other aspect of the contact mechanics, which are outside the scope of the present work. It is unlikely that the assumption of a rigid cylinder with a small effective hydrodynamic gap is a reasonable approximation for the contact between a deformable tyre and the road. Additionally, the Reynolds numbers for studies on tyres are generally orders of magnitude higher than for the present study, which may lead to a reduced influence of the contact point on the drag coefficient.
A recent experimental investigation (Nanayakkara et al. Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a ) of finite-span circular cylinders rolling on an inclined plane found that the drag coefficient can be predicted by assuming an effective hydrodynamic gap equal to the peak surface roughness. However, they noted that expressions for the gap-independent force and moment coefficients are currently lacking in the literature, limiting the accuracy of their analysis. The present work aims to fill this knowledge gap by providing numerical computations of the gap-independent force and moment coefficients.
The structure of this paper is as follows. First, in § 1.1 we define the problem to be solved. Then, in § 2 we discuss the lubrication solution to the inner flow and the computation of the gap-dependent force and moment coefficients. Section 3 then discusses the numerical solution of the outer flow and the evaluation of the gap-independent force and moment coefficients. Finally, concluding remarks are made in § 4.
1.1. Problem description
The problem under consideration is as shown in figure 1. A circular cylinder of diameter
$d$
and span
$W$
travels along a plane wall with translational velocity
$U$
and rotational velocity
$\Omega$
, while maintaining a gap
$G$
between the wall and the cylinder. The cylinder is immersed in a fluid that has density
$\rho$
and dynamic viscosity
$\mu$
.

Figure 1. Problem considered in this work, showing (a) oblique, (b) side and (c) top views. A finite cylinder of diameter
$d$
and span
$W$
translates parallel to a plane wall with velocity
$U$
and angular velocity
$\Omega = 2kU/d$
, where
$k$
is the slip coefficient, while maintaining a gap
$G$
between the cylinder and the wall. Both Cartesian
$(x,y,z)$
and cylindrical
$(r,\theta, z)$
coordinate systems are used.
This flow is non-dimensionalised by
$d$
,
$U$
and
$\rho$
. Then, the flow is characterised by the following non-dimensional parameters: the Reynolds number
$Re = \rho U d/\mu$
, the gap-diameter ratio
$G/d$
, the aspect ratio
$A = W/d$
and the slip coefficient
$k = \Omega d/(2U)$
. We refer to this non-dimensionalisation as the ‘outer-flow’ non-dimensionalisation, or non-dimensionalisation by outer-flow variables, noting that we will introduce a new set of non-dimensional parameters for the inner lubrication flow in § 2.
The aim of the present work is to determine the functional dependence of the force and moment coefficients,

against
$Re$
,
$k$
,
$G/d$
and
$A$
. Here,
$D$
,
$L$
and
$M$
are the drag, lift and moment applied to the cylinder, respectively, and
$C_D$
,
$C_L$
and
$C_M$
are the corresponding force and moment coefficients.
We obtain expressions for the force and moment coefficients using the method of matched asymptotic expansions, which was previously applied to the 2-D cylinder flow by Terrington et al. (Reference Terrington, Thompson and Hourigan2023). Under this approach, we conceptually decompose the flow into separate inner and outer flows, which are matched asymptotically. The inner flow is the flow in the thin interstice between the cylinder and the wall, and solutions to the inner flow are obtained using lubrication theory. The outer flow represents the flow far from the interstice, and is approximately independent of
$G/d$
. In this work, the inner-flow solution is obtained by numerically solving the Reynolds equation in two dimensions, while the outer-flow solution is obtained using 3-D numerical simulations. The main benefit of this approach is that the influence of
$G/d$
can be treated semi-analytically, and the parameter space required for full 3-D simulations is reduced to three parameters (
$A$
,
$k$
and
$Re$
), greatly reducing the computational effort required to explore the parameter space.
Using this approach, we decompose the total force and moment coefficients into gap-dependent (
$G$
) and gap-independent (
$C$
) parts:

The gap-dependent parts of the force and moment coefficients are obtained using lubrication theory in § 2. We obtain expressions of the following form:



Here
$f(\hat {A})$
is determined using the numerical solution of the inner lubrication flow, and
$\hat {A} = A/\sqrt {G/d}$
is the aspect ratio normalised by inner-flow variables. The gap-independent parts of the force and moment coefficients are then obtained from the outer-flow solution in § 3.
2. Lubrication theory
In this section we obtain expressions for the gap-dependent part of the force and moment coefficients by considering the lubrication flow in the thin interstice between the cylinder and the wall. The outer flow is not considered until § 3. This approach was previously applied to the 2-D (or infinite) cylinder flow, by both the studies of Merlen & Frankiewicz (Reference Merlen and Frankiewicz2011) and Terrington et al. (Reference Terrington, Thompson and Hourigan2023). An approximate expression for the gap-dependent drag and moment for the finite-span cylinder was obtained by Saintyves et al. (Reference Saintyves, Rallabandi, Jules, Ault, Salez, Schönecker, Stone and Mahadevan2020), in the limit of a large aspect ratio. This analysis was extended to finite cylinders of any aspect ratio by Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022), including the limiting cases of large (thin rod) and small (thin disk) aspect ratios.
In this section we develop an alternative formulation of the gap-dependent force and moment coefficients, by considering the inner lubrication flow. The new solution is equivalent to Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) solution, differing only by a constant term. This new solution is easier to couple with the numerical simulation to the outer flow, which is considered in § 3. Additionally, our solution more explicitly reveals the dependence of the force and moment coefficients on the parameter
$\hat {A} = A/\sqrt {G/d}$
, while Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) solution appears to depend separately on
$A$
and
$G/d$
.
The structure of this section is as follows. In § 2.1 we introduce the equations and boundary conditions for the inner lubrication flow. Then, in § 2.2 we derive expressions for the gap-dependent force and moment coefficients. These expressions for the gap-dependent force and moment coefficients depend on a function
$f(\hat {A})$
, which is evaluated numerically in § 2.3. Finally, in § 2.4 we compare our results to Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) solution.
2.1. Lubrication solution
This subsection discusses the lubrication solution to the inner flow. The aim is to obtain solutions to the pressure in the inner region as a function of
$G/d$
,
$A$
,
$Re$
and
$k$
. The assumptions of lubrication theory are that fluid flow is confined to a thin film, in which inertial effects are negligible, pressure is constant across the film and cross-film derivatives of velocity are much larger than streamwise derivatives. These conditions are satisfied in the inner-flow region, so long as
$G/d$
is small. Combined with the assumption of steady flow in the inner region, one obtains the Reynolds equation, which is expressed in terms of outer-flow variables as

where
$\boldsymbol {\nabla }_{xz}$
is the 2-D gradient operator in the
$x$
-
$z$
plane,
$h$
is the thickness of the lubrication film and
$\boldsymbol {U}_1$
and
$\boldsymbol {U}_2$
are the velocities of the upper and lower walls, respectively.

Figure 2. Geometry of the inner lubrication region. Here,
$h$
is the film thickness, while
$\boldsymbol {U_1}$
and
$\boldsymbol {U_2}$
are the velocities of the wall and cylinder, respectively.
The geometry for the inner region is as shown in figure 2. The film thickness is approximated by a parabolic profile,

and this approximation is valid when
$x \ll 1$
. Near the contact point, the velocity of the wall is
$\boldsymbol {U_1} = -U \mathrm{\boldsymbol {\hat {e}}}_{x}$
, while the velocity of the cylinder is
$\boldsymbol {U_2} = -k U \mathrm{\boldsymbol {\hat {e}}}_{x}$
, where
$\mathrm{\boldsymbol {\hat {e}}}_x$
is the unit vector in the
$x$
direction. Substituting these expressions into the Reynolds equation gives the following equation for the lubrication pressure:

Equation (2.3) is non-dimensionalised by outer-flow variables, and the solution depends on the four parameters
$G/d$
,
$Re$
,
$k$
and
$A$
. This equation can be further simplified by using the following ‘inner-flow’ variable transformations (Terrington et al. Reference Terrington, Thompson and Hourigan2023):



This gives the following:

The boundary conditions for this equation are
$\hat {p} = 0$
at
$\hat {x} = \pm \infty$
and
$\hat {p} = 0$
at
$\hat {z} = \pm \hat {A}/2$
(Teng et al. Reference Teng, Rallabandi, Stone and Ault2022), since the pressure outside the lubrication region is negligible compared with the pressure within the lubrication region.
Note that the dimensionless pressure distribution
$\hat {p}(\hat {x},\hat {z})$
depends on a single dimensionless parameter,
$\hat {A}$
. In the limit
$\hat {A} \rightarrow \infty$
(i.e. 2-D flow) we can neglect
$\hat {z}$
derivatives, and (2.5) reduces to

The solution to this equation is the lubrication pressure for the 2-D cylinder flow (Jeffrey & Onishi Reference Jeffrey and Onishi1981; Merlen & Frankiewicz Reference Merlen and Frankiewicz2011)

We find it convenient to express the dimensionless pressure distribution for finite
$\hat {A}$
as

where
$\hat {g}$
is a modification to the 2-D pressure distribution due to finite-span effects. Equation (2.5) can be rearranged to give

Due to the pressure boundary conditions, we have
$\hat {g} = 0$
at
$\hat {z} = \pm \hat {A}/2$
. The boundary condition at
$\hat {x} = \pm \infty$
is less clear, since
$\hat {p}_{2D} = 0$
at
$\hat {x}=\pm \infty$
. However, Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) show that
$\hat {p}$
decays as
$\hat {x}^{-5}$
for large
$\hat {x}$
, and therefore,
$\hat {g}$
will decay as
$\hat {x}^{-2}$
. Therefore, we use the boundary condition
$\hat {g} = 0$
at
$\hat {x} = \pm \infty$
.
In this work we obtain profiles of
$\hat {g}(\hat {x},\hat {z})$
numerically for various values of the parameter
$\hat {A} = A/\sqrt {G/d}$
, using a finite-difference method implemented in MATLAB. Details of this implementation are provided in Appendix A.
Figure 3 presents profiles of
$\hat {p}$
computed using the 2-D finite-difference solution to the lubrication flow and from full 3-D numerical simulations (using the numerical scheme described in § 3.2) at
$\hat {A} = 10$
and
$\hat {A} = 100$
. Three-dimensional simulations are performed for
$Re = 50$
,
$A = 1$
and
$k = 1$
, and
$G/d= 10^{-4}$
for
$\hat {A}= 100$
and
$G/d = 10^{-2}$
for
$\hat {A} = 10$
. Excellent agreement is observed between subfigures (a) and (b), confirming that the lubrication solution is valid in the inner region, when the gap is small. Reasonable agreement is observed between subfigures (c) and (d), however, there are several small differences. Notably, the pressure distribution in subfigure (d) is not completely antisymmetric about
$\hat {x} = 0$
and the pressure does not decrease to zero at the boundaries. These differences are due to the larger gap-diameter ratio (
$G/d = 10^{-2}$
) used for subfigure (d), so the assumptions used to obtain the lubrication solution (inertia is negligible and the film thickness is parabolic) are less accurate. Regardless, the inner solution remains a reasonable approximation for the interstitial pressure distribution in this case.

Figure 3. Contours of dimensionless pressure
$\hat {p}$
for (a,b)
$\hat {A} = 100$
and (c,d)
$\hat {A} = 10$
, computed using (a,c) lubrication theory and (b,d) full 3-D numerical simulations. For the 3-D simulations, physical parameters were
$Re = 50$
,
$k = 1$
,
$A = 1$
and (b)
$G/d = 10^{-4}$
for
$\hat {A} = 100$
and (d)
$G/d = 10^{-2}$
for
$\hat {A} = 10$
.
The dimensionless pressure profiles presented in figure 3 are consistent with the results of Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022). Away from the ends of the cylinder, the dimensionless pressure
$\hat {p}$
is approximately equal to the 2-D solution
$\hat {p}_{2D}$
. There is a positive pressure peak just upstream of the point of minimum separation between the cylinder and the wall (which is situated at
$\hat {x} = 0$
), at
$\hat {x} = -1/\sqrt {3}$
, and a negative pressure peak just downstream at
$\hat {x} = 1/\sqrt {3}$
. While the pressure is approximately independent of
$\hat {z}$
away from the ends of the cylinder, it rapidly decreases to zero at the ends of the cylinder over a thin end correction region. The width of this end correction region, relative to the total width of the cylinder, decreases as
$\hat {A}$
is increased.
2.2. Force and moment coefficients
We now consider the computation of the gap-dependent part of the force and moment coefficients, using only the inner-flow solution. We assume that the inner and outer solutions are matched asymptotically, and that both inner and outer approximations are valid for some angular position
$\theta = \theta _0$
. Then, the inner solution is applied for
$\theta \leqslant \theta _0$
, while the outer solution is applied for
$\theta \geqslant \theta _0$
. Without directly computing the outer-flow solution, we obtain the
$G/d$
dependence of the total force and moment coefficients based on the assumption of asymptotic matching between the inner and outer solutions.
The structure of this section is as follows. First, in § 2.2.1 the force and moment coefficients are expressed as contributions from the inner and outer solutions. Next, in § 2.2.2, expressions for the inner-flow force and moment coefficients are obtained in terms of
$\hat {g}(\hat {x},\hat {z})$
. Then, in § 2.2.3 the
$\theta _0$
dependence of the outer-flow force and moment coefficients is deduced from the inner-flow solution, without directly computing the outer-flow solution. Finally, expressions for the
$G/d$
-dependent force and moment coefficients are obtained in § 2.2.4.
2.2.1. Contributions from the inner and outer regions
The total force and moment coefficients are expressed as



where
$C_{D,F}$
,
$C_{L,F}$
and
$C_{M,F}$
are the force and moment contributions from the end faces of the cylinder:



Here,
$(r,\theta, z)$
are cylindrical coordinates as defined in figure 1, overbars denote the mean across the cylinder span of quantities on the cylinder’s curved surface:

Here, subscript
$F$
denotes the sum of quantities on the cylinder’s front and back faces:

The total force and moment coefficients may be separated into contributions from the inner and outer regions, as

where a subscript
$[i]$
in terms such as
$C_{[i]}$
represents either the drag (
$C_D$
), lift (
$C_L$
) or moment (
$C_M$
) coefficient, as required. The outer-flow contributions are



while the contributions from the inner flow are



where the small angle approximation has been used in (2.16a
–
c
). Here,
$\theta _0$
is the boundary between the inner and outer regions, and
$\hat {x}_0 = \sin {\theta _0}/(2\sqrt {G/d}) \approx \theta _0/(2\sqrt {G/d})$
is this same boundary normalised by inner-flow variables. We note that both the inner and outer force and moment coefficients are functions of
$\theta _0$
. However, the total force and moment coefficients should be independent of
$\theta _0$
, assuming that the inner and outer solutions are matched asymptotically.
2.2.2. Force and moment contributions from the inner region
From (2.4a ), (2.8) and (2.12), the pressure in the inner region is expressed as

Under the assumptions of lubrication theory, the wall-shear stress in the inner region is (Ghosh et al. Reference Ghosh, Majumdar and Sarangi2014)

and after substituting (2.17) into (2.18), the mean stress is

while the
$y$
-shear stress is negligible.
Finally, substituting (2.17) and (2.19) into (2.16) gives the following expressions for the inner-flow force and moment coefficients:





2.2.3. Asymptotic matching to the outer flow
Without directly computing the outer-flow solution, we determine the
$\theta _0$
dependence of the outer-flow force and moment coefficient, assuming that the inner and outer solutions are matched at
$\theta = \theta _0$
.
First, we consider the asymptotic behaviour of
$\hat {p}$
for large
$\hat {x}$
. Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) provide an expression for this asymptotic pressure
$\hat {p}_o$
(their (2.25)), which we rewrite as a modification to the 2-D pressure
$\hat {g}_o$
:



Here,
$\hat {p}_o$
and
$\hat {g}_o$
are asymptotic expressions for the terms
$\hat {p}$
and
$\hat {g}$
(defined in (2.7) and (2.8)) in the limit
$\hat {x}\gg 1$
. Also,
$\hat {g}_o$
and
$\hat {p}_o$
are related as
$\hat {p}_o = \hat {g}_o \hat {p}_{2D,o}$
, where
$\hat {p}_{2D,o} = -1/\hat {x}^3$
is the asymptotic value of
$\hat {p}_{2D}$
for large
$\hat {x}$
. Averaging across the cylinder span gives the following expression for
$\bar {\hat {g}}_o$
:

Note that
$\hat {g}_o$
in (2.21) is a function of the parameters
$\hat {z}/\hat {A}$
and
$\hat {x}/\hat {A}$
, rather than of
$\hat {x}$
and
$\hat {z}$
alone. These parameters are equivalent to
$x/A$
and
$z/A$
and are independent of
$G/d$
. Therefore, the pressure far from the interstice is independent of
$G/d$
.

Figure 4. Profiles of
$\bar {\hat {g}}$
obtained using the finite-difference lubrication solution, along with the analytic outer-flow solution
$\bar {\hat {g}}_o$
against (a)
$\hat {x}$
and (b)
$x=\hat {x}\sqrt {G/d}$
. The aspect ratio is held constant at
$A = 1$
, and profiles of
$\bar {\hat {g}}$
are computed for a range of
$G/d$
.
Figure 4 presents profiles of both
$\bar {\hat {g}}_o$
obtained using (2.24) and of
$\bar {\hat {g}}$
computed using the finite-difference solution of (2.9) for a range of
$\hat {A}$
. Subfigure (a) presents the profiles against
$\hat {x}$
, which confirms that the inner-flow solutions (
$\bar {g}$
) approach the asymptotic solution
$\bar {\hat {g}}_o$
when
$\hat {x} \gg 1$
. Subfigure (b) presents the same profiles, but against
$x = \hat {x}\sqrt {G/d}$
(i.e. in outer-flow variables). The profiles of
$\hat {g}$
for different
$\hat {A}$
collapse onto a single universal outer profile given by
$\bar {\hat {g}}_o$
for large
$x$
, confirming that the pressure far from the interstice is independent of
$G/d$
.
The asymptotic value of the spanwise-averaged lubrication pressure for large
$\hat {x}$
can be recast into outer-flow variables as

where we have assumed that
$x \approx \theta /2$
due to the small angle approximation. Similarly, the asymptotic value of the wall-shear stress is

We propose that the outer-flow force and moment coefficients be decomposed as the sum of a
$\theta _0$
-dependent and constant term,

The terms
$C_{D,C}$
,
$C_{L,C}$
and
$C_{M,C}$
are independent of
$\theta _0$
, and are referred to as the gap-independent force and moment coefficients. These terms cannot be determined from the lubrication solution alone, and are evaluated using 3-D numerical simulations in § 3.
The terms
$C_{D,\theta _0}$
,
$C_{L,\theta _0}$
and
$C_{M,\theta _0}$
in (2.27) are referred to as the
$\theta _0$
dependent part of the outer-flow force and moment coefficients. These terms can be determined from the inner solution alone, up to an arbitrary constant, by requiring that the sum of the inner-flow and
$\theta _0$
-dependent force/moment coefficient be constant with respect to
$\theta _0$
. We define this arbitrary constant using the condition
$C_{[i],\theta _0}(\theta _0 = \infty ) =0$
, so that we have

and by substituting (2.16) into (2.28), we obtain the following expressions for the
$\theta _0$
-dependent force and moment coefficients:



Here, note that integrals have been simplified using symmetry. Finally, substituting (2.25) and (2.26) into (2.29) gives the following expressions for the
$\theta _0$
-dependent part of the outer-flow force and moment coefficients:




2.2.4. Gap-dependent force and moment coefficients
We define the gap-dependent force and moment coefficients as the sum of the inner flow and
$\theta _0$
-dependent part of the outer-flow force and moment coefficients:

Then, the total force and moment coefficients are the sum of the gap-dependent and gap-independent terms:

Specifically, the gap-independent terms are constant with respect to
$G/d$
, but may be functions of
$A$
,
$Re$
and
$k$
. These constant terms are determined using direct numerical simulations of the outer flow in § 3.
Using the condition
$C_{[i],\theta _0}(\theta _0 = \infty ) = 0$
, the gap-dependent drag can be expressed as




Using integration by parts on the final term in both
$\hat {f}_x$
and
$\hat {m}_z$
,

we obtain

where the factor of
$2\pi$
is introduced to simplify the notation, so that
$0\leqslant f(\hat {A}) \leqslant 1$
.
Finally, the gap-dependent force and moment coefficients are expressed as




2.3. Numerical evaluation of
$f(\hat {A})$
The function
$f(\hat {A})$
in (2.36) represents the effects of a finite aspect ratio on the gap-dependent part of the drag and moment coefficients. Specifically, for 2-D flow (
$\hat {A} = \infty$
), the value of this function is
$f(\infty ) = 1$
, and the force and moment coefficients equal the 2-D solutions (Merlen & Frankiewicz Reference Merlen and Frankiewicz2011)

For finite
$\hat {A}$
,
$f(\hat {A})$
is less than one, indicating a reduction in the force and moment coefficients compared with the 2-D solution.
In this study,
$f(\hat {A})$
is determined numerically. We first obtain profiles of
$\hat {g}$
for a particular value of
$\hat {A}$
by numerically solving (2.9) using the finite-difference approach described in Appendix A. Then, the function
$f(\hat {A})$
is evaluated at this
$\hat {A}$
by numerical integration of (2.36d
). We consider
$61$
logarithmically spaced values of
$\hat {A}$
between
$\hat {A} = 10^{-2}$
and
$\hat {A} = 10^{4}$
. For each of these
$\hat {A}$
, solutions of
$\hat {g}$
and
$f$
were obtained using three different grid resolutions as listed in table 1, and a Richardson extrapolation was performed to estimate the grid-independent solution. We note that the maximum relative difference between the predictions of
$f(\hat {A})$
on the finest mesh and the Richardson extrapolation is only
$1.34 \times 10^{-6}$
, and therefore, the meshes are fine enough to accurately compute
$f(\hat {A})$
.
Table 1. Resolution study for numerical computations of
$f(\hat {A})$
. Here
$n_x$
and
$n_z$
are the number of grid points in the
$x$
and
$z$
directions, respectively. The maximum relative error is computed with respect to the estimated grid-independent solution obtained using a Richardson extrapolation.


Figure 5. Profile of the function
$f(\hat {A})$
(defined in (2.36d
)), obtained from 2-D simulations of the lubrication flow. Markers indicate the numerical data, while the solid line is a cubic-spline interpolation between the data points.
A cubic-spline interpolation is used to obtain values of
$f(\hat {A})$
between the numerical data points. To ensure that a sufficient number of points are used in the spline interpolation, a second cubic-spline interpolation was created using half of the data points. The maximum difference between the two interpolations is only
$4.8\times 10^{-4}$
, and therefore, our spline interpolation is sufficiently accurate.
Figure 5 presents the function
$f(\hat {A})$
, with markers indicating the numerical data, and the solid line indicating the cubic-spline fit. As expected,
$f(\hat {A})$
approaches 1 as
$\hat {A}$
approaches
$\infty$
, corresponding to 2-D flow. For finite
$\hat {A}$
,
$f$
is between
$0$
and
$1$
, indicating a reduction in the gap-dependent force and moment coefficients compared with Stokes flow. The function
$f(\hat {A})$
is monotonically increasing. Increasing
$\hat {A}$
– either by increasing
$A$
or reducing
$G/d$
– results in force and moment coefficients closer to the 2-D theory, while reducing
$\hat {A}$
results in a greater relative difference between the 3-D and 2-D force and moment predictions.
Asymptotic expressions for
$f(\hat {A})$
for both small and large
$\hat {A}$
are obtained in Appendix B. For small
$\hat {A}$
, we obtain the asymptotic value

and this is confirmed in figure 6(a) by the close agreement between the numerical data and the asymptotic expression. For large
$\hat {A}$
, we obtain the expression

where
$\ln ()$
is the natural logarithm function, and the coefficient of the final term has been estimated by fitting to our numerical data. As shown in figure 6(b), the numerical data converge to this asymptotic expression for
$\hat {A}\gg 1$
. For
$\hat {A}\gt 10^3$
,
$f(\hat {A})\gt 0.99$
, so that 2-D lubrication theory holds to less than
$1\%$
error.

Figure 6. Asymptotic behaviour of the function
$f(\hat {A})$
for (a) small
$\hat {A}$
and (b) large
$\hat {A}$
, along with the asymptotic expressions obtained in Appendix B.
For the remainder of this work, we estimate
$f(\hat {A})$
as follows. Within the range
$10^{-2} \leqslant \hat {A}\leqslant 10^{4}$
, a cubic-spline interpolation fit to our numerical data is used. Outside this range, the asymptotic expression in (2.38) is used for
$\hat {A}\lt 10^{-2}$
, while (2.39) is used for
$\hat {A} \gt 10^{4}$
. A `.mat’ file containing the numerical data for
$f(\hat {A})$
is provided in the online supplementary material, along with MATLAB and Python scripts for estimating
$f(\hat {A})$
using the combined cubic-spline/asymptotic method.
2.4. Comparison to Teng et al
Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) have also computed the gap-dependent part of the force and moment coefficients for a finite-span cylinder rolling near a wall. In this subsection we compare our results to their solution. While the form of our equations (2.36) is quite different from their equations (2.17), we find that the two expressions differ only by a constant term with respect to
$G/d$
. Since the gap-dependent force and moment coefficients are only unique up to the addition of a constant term, this means that our solution is validated against the Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) solution.
Importantly for the present work, our derivation provides us with an expression for the
$\theta _0$
-dependent force and moment coefficients (2.30), which are used to compute the gap-dependent force and moment terms (
$C_{D,C}, C_{M,C}$
) in § 3. It would also appear that this form is simpler, since the effects of finite aspect ratio are described by a single parameter,
$\hat {A}$
, rather than depending separately on both
$A$
and
$G/d$
, as is the case for the Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) formulation.
Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) expressions for the force and moment can be expressed in terms of drag and moment coefficients as

where
$\mathcal I$
is the integral

Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) give the approximate solution for
$\mathcal{I}$
as

which is valid for
$G/d \ll 1$
. We emphasise that (2.40) does not consider the influence of the outer flow, and therefore, like our equation (2.36), describes only the gap-dependent part of the drag, due to lubrication effects in the inner region. These terms are only defined up to an additive constant, and therefore, (2.40) may differ from (2.36) by an additive constant.

Figure 7. Comparison between our numerically obtained
$C_{D,G}$
(see (2.36)), and Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) solution (see (2.40)) evaluated using both direct numerical simulations (see (2.41)), and the analytic approximation (see (2.42)). A range of
$G/d$
and
$A$
are considered, while the remaining parameters are fixed at
$Re = 50$
and
$k = 1$
. Note that curves for
$A\geqslant 10$
overlap in subfigure (a), while curves for
$A\leqslant 0.1$
overlap in subfigure (b).
Figure 7 presents the computed drag coefficients against
$G/d$
for a range of aspect ratios, using both Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) solution and our solution (2.36). For Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) solution, we have calculated the integral
$I$
both by using the analytic approximation (2.42) and by evaluating (2.41) directly using our numerically obtained solution for
$\hat {p}$
. For
$A \geqslant 0.1$
, the analytical approximation is in excellent agreement with the numerical computation, however, (2.42) is inaccurate when
$A\lt 0.1$
.
Figure 7(a) shows that our solution for the drag coefficient (2.36) displays the same
$G/d$
dependence as Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) solution, however, it differs by a constant with respect to
$G/d$
. This is more clearly shown in figure 7(b), which plots the difference between our solution and Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) solution (computed by directly evaluating (2.41) from the numerical solution). For a fixed
$A$
, the two solutions differ only by a constant term, which is a function of
$A$
, and likely also depends on both
$k$
and
$Re$
. This is expected, since expressions for the gap-dependent drag coefficients (
$C_{D,G}$
) are unique up to the addition of a term that does not depend on
$G/d$
.
In § 3 of this paper we obtain the gap-independent force and moment coefficients (
$C_{D,C}$
,
$C_{L,C}$
and
$C_{M,C}$
), which are added to the gap-dependent force and moment coefficients (
$C_{D,G}$
,
$C_{L,G}$
and
$C_{M,G}$
) to obtain the total force and moment coefficients. Since Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) expressions differ from ours by an additive constant, the corresponding gap-independent force and moment coefficients presented in § 3 need to be adjusted by this same additive constant if one uses Teng et al.’s (Reference Teng, Rallabandi, Stone and Ault2022) expressions for the gap-dependent force and moment coefficients, instead of our solution.
2.5. Applicability of the asymptotic matching assumption
We now consider the limitations of the approximations developed in this section, to inform the parameter space for which the model developed in this paper is applicable. First, the parabolic approximation (2.2) is valid for
$x\ll 1$
(or, alternatively,
$\theta \ll 1$
or
$\hat {x} \ll \sqrt {G/d}$
). The outer flow must be independent of
$G/d$
, which occurs for
$\hat x \gg 1$
, or equivalently,
$\theta \gg \sqrt {G/d}$
. The conditions of asymptotic matching require that
$\theta _0$
lies in a region where both inner and outer solutions are valid, which therefore requires
$\sqrt {G/d} \ll \theta _0 \ll 1$
. This condition is satisfied when
$\sqrt {G/d} \ll 1$
.
An additional consideration is that inertial effects are assumed negligible in the inner region. The film thickness at
$\theta = \theta _0$
is approximately
$h = 0.5 \theta _0^2$
, and the corresponding Reynolds number based on the film thickness is
$Re_h = 0.5 \theta _0^2 Re$
. The lubrication assumption is valid for
$Re_h \ll 1$
, giving a further restriction on
$\theta _0$
:
$\sqrt {G/d} \ll \theta _0 \ll 1/\sqrt {Re}$
, from which we obtain the condition
$G/d \ll 1/Re$
.
For high
$Re$
, we expect the gap-independent force and moment coefficients to remain of the order of
$O(1)$
, since these are dominated by inertial effects from the outer flow. Since the gap-independent force and moment coefficients scale as
$1/(Re\sqrt {G/d})$
(2.36), we obtain the following considerations at high
$Re$
: gap-dependent effects dominate for
$G/d \ll 1/Re^2$
, both gap-dependent and gap-independent effects will be significant for
$G/d = O(1/Re^2)$
, and gap-independent (outer-flow) effects dominate for
$G/d \gg O(1/Re^2)$
.
As a final consideration, Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) showed that, for a finite-gap height, the pressure does not decrease to zero at
$z = A$
, but extends past the ends of the cylinder by a distance proportional to
$G/d$
. As a result, the inner-flow solution is only accurate for
$G/d \ll A$
. For thin disks (
$A \ll 1$
), the
$G/d$
required to satisfy this assumption becomes very small.
3. Outer flow
In this section the outer flow in the region far from the gap is obtained using numerical simulations. Then, the gap-independent force and moment coefficients (
$C_{D,C}$
,
$C_{L,C}$
and
$C_{M,C}$
) are computed using the numerically obtained outer-flow solution. The structure of this section is as follows. First, in § 3.1 we discuss the numerical methodology used to obtain outer-flow solutions. Then, in § 3.2 we validate our approach against single-domain numerical simulations performed with a finite-gap ratio. Next, in § 3.3 we discuss the computation of the gap-independent contribution to the force and moment coefficients from the outer flow. Then, in § 3.4, flow visualisations are presented to discuss the effects of
$Re$
,
$A$
and
$k$
on the wake structures. Next, in § 3.5 we report the variation of gap-independent force and moment coefficients against
$Re$
,
$A$
and
$k$
. Then, in § 3.6 we compare the drag coefficient for a cylinder in contact with a wall to that of a cylinder in a uniform free-stream flow. Next, in § 3.7 we present the variation of total force and moment coefficients against
$A$
,
$Re$
,
$k$
and
$G/d$
. Finally, in § 3.8 we compare the results of the present study to experimental measurements.
3.1. Numerical set-up
This subsection discusses the numerical methods used to obtain solutions to the outer flow. The computational domain is as illustrated in figure 8. There are four external boundaries: the inlet at
$x = -20$
, the outlet at
$x = 50$
, the wall at
$z = 0$
and the far-field boundaries at
$y = \pm 50$
and
$z = 50$
, respectively.
Figure 8(b) shows the geometry near the cylinder. The outer-flow solution is independent of
$G/d$
and obtained under the assumption
$G/d = 0$
(Terrington et al. Reference Terrington, Thompson and Hourigan2023). To avoid the infinite pressure at the contact point, a small region near the contact point is cut out from the computational domain. The size of this cutout region is governed by the cutout angle
$\theta _c$
. The cutout introduces four inlet/outlet boundaries, on the front (
$I_f$
), back (
$I_b$
), inlet (
$I_i$
) and outlet (
$I_o$
) sides of the cylinder, respectively. Finally, the cylinder has three boundary surfaces, namely the front (
$C_f$
) and back (
$C_b$
) faces, and the curved surface
$C_c$
.

Figure 8. Computational domain for outer-flow simulations, showing (a) the outer boundaries and (b) a closeup of the cylinder, featuring a cutout region to avoid infinite pressures at the contact point. Solid arrows indicate the velocity boundary conditions at the inlet and on the wall.
Numerical simulations are performed in a moving reference frame, in which the translational velocity of the cylinder is zero. Then, the boundary conditions are as follows. The inlet is set to a constant velocity
$(u_x,u_y,u_z) = (1,0,0)$
, the outlet is set to a constant pressure
$p = 0$
and the far-field boundaries are set to shear-free walls. The no-slip condition (meaning no slip between the fluid and the solid boundary) is applied to both the plane wall (
$(u_x,u_y,u_z) = (1,0,0)$
) and the cylinder (
$(u_x,u_y,u_z) = (k \cos \theta, k\sin \theta, 0)$
on
$C_f$
,
$C_b$
and
$C_c$
). A constant pressure
$p = 0$
is applied to both
$I_f$
and
$I_b$
. Finally, the pressure at both
$I_i$
and
$I_o$
is given by (2.25), under the assumption that the outer-flow pressure approaches the outer-flow lubrication solution in the limit
$\theta \rightarrow 0$
(Terrington et al. Reference Terrington, Thompson and Hourigan2023).
Numerical simulations were performed using the commercial finite-volume solver ANSYS FLUENT. The momentum and continuity equations were discretised using the second-order upwind scheme for all spatial derivatives, and the second-order backwards time scheme was used for transient simulations. The SIMPLEC scheme was used for pressure--velocity coupling, and the resulting pressure and momentum equations were solved using the algebraic multigrid method. A time step of
$\Delta t = 0.01$
was used. While this resulted in a high Courant number in the inner region (
$Co\gt 25$
), due to the small cell size there, flow is steady in the inner region, so time accuracy errors due to a high Courant number are not an issue. The maximum Courant number in the outer region was
$Co \approx 0.63$
.

Figure 9. Illustration of the block-mesh scheme used in this study, showing the meshes for both (a) the outer region and (b) the inner region. Faces are coloured to indicate the plane wall (grey), the cylinder (green) and the interface between the two meshes (blue). The locations of the representative cell sizes (
$\Delta x_{i,min}$
,
$\Delta z_{i,min}$
,
$\Delta z_{i,max}$
,
$\Delta z_{o,min}$
,
$\Delta z_{o,max}$
) and the number of cells across various edges (
$N_{c,o}$
,
$N_{c,i}$
,
$N_y$
) used in table 2 are also indicated.
The computational domain was meshed with a block-structured approach, using the commercial software package ANSYS ICEM CFD. It was found that the region near the contact point required a much finer grid resolution than the rest of the computational domain, likely due to large pressure gradients near the contact point. Therefore, separate meshes were generated for the region near the contact point and the outer-flow region, as illustrated in figure 9. The two meshes were then coupled using a conformal mesh interface in ANSYS FLUENT.
The present numerical method was validated by performing numerical simulations of flow past a finite-span rotating cylinder in a uniform free-stream flow. Figure 10 presents a comparison between the present numerical scheme and the numerical study of Yang et al. (Reference Yang, Wang, Guo and Zhang2023), showing the variation of lift (
$C_L$
) and drag coefficients (
$C_D$
) against rotation rate
$k$
, at
$Re = 200$
and
$A = 1$
. Good agreement between the present solver and Yang et al. (Reference Yang, Wang, Guo and Zhang2023) are observed, with differences below
$2\,\%$
.

Figure 10. Comparison between the present numerical solver and the numerical results of Yang et al. (Reference Yang, Wang, Guo and Zhang2023) for flow over a rotating finite-span circular cylinder in a uniform free-stream flow at
$Re = 200$
and
$A = 1$
.
A mesh resolution study was performed to ensure the numerical results are sufficiently accurate. Table 2 lists the different meshes used, along with a range parameters used to characterise the grid resolution. These are as indicated in figure 9:
$N_i$
and
$N_o$
, the total number of cells used in the inner and outer meshes;
$\theta _c$
, the cutout angle;
$N_y$
, the total number of cells across the film thickness in the inner region;
$N_{c,i}$
and
$N_{c,o}$
, the number of cells along the cylinder circumference, in both the inner and outer regions; and
$\Delta x_{i,min}$
,
$\Delta z_{i,min}$
,
$\Delta z_{i,max}$
,
$\Delta z_{o,min}$
,
$\Delta z_{o,max}$
, which are various representative cell sizes on the cylinder surface.
Table 2. List of various meshes used in the grid resolution study, along with a variety of parameters used to characterise the grid resolution.

The resolution study is performed for
$Re = 200$
,
$k = 1$
and
$A = 1$
. For these parameters, the wake is unsteady and non-periodic (see figure 28
d). To obtain unsteady flow statistics, simulations were run for
$t = 200$
dimensionless time units prior to collecting statistics, to allow the flow to reach a statistically steady state. Then, simulations were run for a further
$500$
time units to obtain the mean and root-mean-square (r.m.s.) force and moment coefficients. To estimate the influence of the sampling window size, a single case (using mesh 5) was run for an additional
$500$
time units. The changes in the time-mean gap-independent drag and moment coefficients (
$\overline {C_{D,C}}$
and
$\overline {C_{M,C}}$
) were approximately
$0.1\,\%$
, while changes to the time-averaged gap-dependent lift coefficient (
$\overline {C_{L,C}}$
) were approximately
$1\%$
. Changes to the r.m.s. force and moment coefficient were more significant (
$4\%$
for
$C_{D,{rms}}$
,
$3.5\%$
for
$C_{M,{rms}}$
and
$8\%$
for
$C_{L,{rms}}$
). Therefore, while the time-mean force and moment coefficients are well resolved, the r.m.s. force and moment coefficients are not. We consider this to be acceptable, as the present study primarily focuses on the time-mean force and moment coefficients.
Table 3 lists the time-mean (
$\overline {C_{D,C}}$
,
$\overline {C_{L,C}}$
,
$\overline {C_{M,C}}$
) and r.m.s. (
$C_{D,{rms}}$
,
$C_{L,{rms}}$
,
$C_{M,{rms}}$
) gap-independent force and moment coefficients obtained for each mesh, which are computed using the methodology outlined in § 3.3. For meshes 1–4, there is an increase in grid resolution over the entire computational domain, with a relative cell spacing
$h_{rel}$
. A clear trend in convergence of the time-mean force and moment coefficients with increasing grid resolution is observed, and a Richardson extrapolation was used to estimate the grid-independent value of these quantities. Mesh 5 features the same outer-flow resolution as mesh 3, but with greater resolution in the inner region. The computation using this mesh predicts the time-averaged drag and moment coefficients with an error of below
$1\%$
, and predicts the time-mean lift coefficient to approximately a
$2\%$
error. Variation in the
$\mathrm{r.m.s.}$
force and moment coefficients between meshes 2–5 is approximately
$5\%$
, which is comparable to the error introduced by the sampling window size. Therefore, mesh 5 is considered sufficient for the remainder of this study.
Table 3. Resolution study for unsteady flow at
$Re = 200$
,
$k = 1$
and
$A = 1$
. A Richardson extrapolation is performed for the mean force and moment coefficients, using meshes
$1$
--
$4$
. The relative error compared with the Richardson extrapolation is given in parentheses.

Mesh 6 has a similar grid resolution to mesh 5, but with a larger cutout angle
$\theta _c = 0.05$
. Variation in the time-averaged drag and moment coefficients between meshes 5 and 6 are negligible (
$\overline {C_{D,C}}$
and
$\overline {C_{M,C}}$
change by less than
$0.1\%$
between meshes 5 and 6). Changes to the mean lift coefficient (
$\overline {C_{L,C}}$
) between meshes 5 and 6 are larger, at
$3.1\%$
. Therefore,
$\theta _c = 0.025$
is sufficiently small to accurately predict the drag and moment coefficients, and provides reasonable predictions of the lift coefficient.
3.2. Comparison to finite
$G/d$
simulations
To validate the asymptotic matching technique presented in this paper, a set of single-domain numerical simulations were performed for finite-gap--diameter ratios (
$G/d$
) between
$10^{-4}$
and
$10^{-1}$
, for
$Re = 200$
,
$A = 1$
and
$k = 1$
. The details of the numerical scheme are similar to the zero-gap outer-flow simulations discussed in § 3.1. However, the finite gap between the cylinder and the wall is included in the numerical domain. Therefore, these simulations consider a single numerical domain that includes both the inner and outer regions.
A grid resolution study was performed for the finite-gap computations. Three meshes were used, with statistics given in table 4. Meshes 1–3 here have a similar resolution in the outer-flow region to meshes 1–3 in table 2. However, these meshes use far more elements in the inner region, in order to resolve the lubrication flow in the interstice. In particular, the minimum cell size in the
$x$
direction (
$\Delta x_{min}$
) is an order of magnitude smaller for the finite-gap simulations, compared with the zero-gap outer-flow computations.
Table 4. List of meshes used in the grid resolution study for the finite-gap simulations, along with a variety of parameters used to characterise the grid resolution. Here,
$N_y$
is the number of cells across the film thickness in the gap-region.

Table 5 lists the time-mean force and moment coefficients obtained on each of the three meshes, for
$Re = 200$
,
$A = 1$
,
$k = 1$
and
$G/d = 10^{-4}$
. A Richardson extrapolation was performed to estimate the grid-independent value of the force and moment coefficients, and mesh 3 has errors of below
$1\%$
for
$C_D$
and
$C_M$
, and
$1.5\%$
for
$C_L$
compared with the grid-independent solution, and is therefore considered sufficient for this study.
Table 5. Resolution study for finite-gap simulations at
$Re = 200$
,
$A = 1$
,
$k = 1$
and
$G/d = 10^{-4}$
.

An estimate of the gap-independent drag and moment coefficients,
$C_{D,C}$
and
$C_{M,C}$
, is also provided in table 5. These are obtained by subtracting the gap-dependent force/moment coefficient (2.36) from the total force/moment coefficient. The gap-independent drag coefficient obtained using mesh 3 is accurate to within
$1\%$
, however, the error in the gap-independent moment coefficient has an error of approximately
$5\%$
. The finite-gap simulations produce relatively large errors in the gap-independent moment coefficient, since this term comprises only
$6\%$
of the total moment coefficient at
$Re = 200$
and
$G/d = 10^{-4}$
.
Comparing the estimated grid-independent values of the gap-independent drag and moment coefficients (tables 3 and 5), there is good agreement between the results of the finite-gap and zero-gap simulations. The zero-gap outer-flow simulation using mesh 5 required 69.1 hours (approximately 3 days) of wall-clock time using 64 CPU processors on a single node of the Setonix HPE located at the Pawsey Supercomputing Research Centre in Western Australia, while the finite-gap simulation using Mesh 3 required 139.8 hours (approximately 6 days). Therefore, the combined analytical--numerical approach is twice as fast as directly simulating the finite-gap case, in addition to more accurately predicting the moment coefficient.
Figure 11 compares the wake behind the cylinder at
$Re = 200$
,
$A = 1$
and
$k = 1$
, for (a,b) the zero-gap outer-flow computation, (c,d)
$G/d = 10^{-4}$
and (e,f)
$G/d = 10^{-2}$
single-domain finite-gap computations. Subfigures (a,c,e) show isosurfaces of time-mean streamwise vorticity
$\bar {\omega }_x = \pm 0.5$
, while subfigures (b,d,f) show isosurfaces of time-mean spanwise vorticity
$\bar {\omega }_z = \pm 0.5$
. The time-averaged flow structures are practically identical between the different values of
$G/d$
, confirming that the outer flow is approximately independent of
$G/d$
, when
$G/d$
is sufficiently small. Moreover, this figure confirms that the zero-gap outer-flow simulation is able to correctly predict the outer flow for finite
$G/d$
.

Figure 11. Isosurfaces of time-mean (a, c, e) streamwise-oriented vorticity
$\bar {\omega }_x = \pm 0.5$
and (b, d, f) spanwise-oriented vorticity
$\bar {\omega }_z = \pm 0.5$
, at
$Re = 200$
,
$k = 1$
and
$A= 1$
, and gap-diameter ratios (a,b)
$G/d = 0$
, (c,d)
$G/d = 10^{-4}$
and (e,f)
$G/d = 10^{-2}$
. Red and blue isosurfaces indicate positive and negative vorticity, respectively.

Figure 12. Profiles of the time-mean span-averaged pressure (
$\overline {p}$
) against angular position
$\theta$
on the cylinder circumferential face (
$C_c$
) for a range of gap-diameter ratios (
$G/d$
), at
$Re = 200$
,
$k = 1$
and
$A = 1$
. Subfigure (a) shows the pressure profiles in the outer region, while subfigure (b) shows the pressure profiles near the contact point. The asymptotic outer-flow pressure profile (
$\bar {p}_o$
) given by (2.24) is also included in (b). The gap region/contact point is located at
$\theta = 0,2\pi$
.
Figure 12 presents profiles of the time-mean span-averaged pressure (
$\bar {p}$
) on the cylinder’s curved surface for a range of
$G/d$
, at
$Re = 200$
,
$A = 1$
and
$k = 1$
. Subfigure (a) is scaled to show the pressure in the outer flow, away from the contact point (located at
$\theta = 0,2\pi$
). Subfigure (b) presents the same curve on a logarithmic axis to show the asymptotic behaviour of the span-averaged pressure near the contact point. For
$G/d \leqslant 10^{-2}$
, the pressure profiles for different
$G/d$
collapse to a single curve, confirming that the outer flow is independent of
$G/d$
, when
$G/d$
is small. The solid black line in subfigure (b) indicates the asymptotic profile for the outer flow obtained from lubrication theory (2.24). As
$\theta$
approaches zero, the outer-flow solutions approach this asymptotic limit. For a finite
$G/d$
, however, the pressure profiles deviate from the asymptotic limit when
$\theta \lessapprox \sqrt {G/d}$
, instead following the lubrication solution to the inner flow. This was confirmed by figure 3, which showed good agreement between the pressure profiles predicted from lubrication theory and those obtained using single-domain finite-gap numerical simulations. Finally, the outer-flow solution obtained for
$G/d = 0$
overlaps with the asymptotic solution for small
$\theta$
.
3.3. Gap-independent force and moment coefficients

Figure 13. Profiles of (a) the force and moment contributions from the outer region (
$C_{[D/M/L],O}$
), as well as the asymptotic profiles (
$C_{[D/M/O],\theta _0}$
), and (b) the difference between the outer flow and asymptotic force and moment coefficients (
$\Delta C_{[D/M/L]}$
) against
$\theta _0$
, for
$Re = 200$
,
$k =1$
and
$A = 1$
.
We now discuss the computation of the gap-independent force and moment coefficients using the numerical solution to the outer flow. Figure 13(a) presents profiles of
$C_{D,O}$
,
$C_{L,O}$
and
$C_{M,O}$
obtained for
$Re = 50$
,
$A = 1$
and
$k = 1$
. These profiles are obtained by numerically integrating (2.15) using the trapezoid rule for various values of
$\theta _0$
. As
$\theta _0$
approaches zero, both
$C_{D,O}$
and
$C_{M,O}$
approach infinity, while
$C_{L,O}$
approaches a finite value. The dashed lines in figure 13(a) indicate the
$\theta _0$
-dependent part of the outer-flow profiles (
$C_{D,\theta _0}$
,
$C_{L,\theta _0}$
and
$C_{M,\theta _0}$
) given by (2.30). The divergence of
$C_{D,O}$
and
$C_{M,O}$
towards infinity appears to be captured by the terms
$C_{D,\theta _0}$
and
$C_{M,\theta _0}$
. To confirm this, we compute the differences between these two terms,

which are plotted in figure 13(b). Each of these terms approaches a finite value as
$\theta _0\rightarrow 0$
.
When
$\theta _0$
is small, computing the terms
$\Delta C_{[i]}$
involves taking the difference of two large numbers, since the pressure approaches infinity near the contact point. This leads to some inaccuracies when calculating
$\Delta C_{[i]}$
for small
$\theta _0$
(Terrington et al. Reference Terrington, Thompson and Hourigan2023). To estimate the behaviour of these functions for small
$\theta _0$
, a polynomial fit is performed over the range
$0.2\leqslant \theta \leqslant 1$
. These polynomial fits are indicated by the black dashed curves in figure 13(b).
Using (2.27) and (2.32), the gap-independent force and moment coefficients are given by

which should be independent of
$\theta _0$
, at least when
$\theta _0$
is small. However, the profiles of
$\Delta C_{[i]}$
plotted in figure 13(b) are not constant with respect to
$\theta _0$
. This indicates that the assumptions used to derive (2.27) and (2.32), such as the assumption that inertial effects are negligible for small
$\theta$
, are not strictly satisfied. This leads to some additional errors in our estimate of the force and moment coefficients.
To estimate these errors, we assume that both
$\Delta C_D$
and
$\Delta C_M$
are approximately quadratic against
$\theta _0$
, while
$\Delta C_L$
is approximately linear in
$\theta _0$
, when
$\theta _0$
is small (Terrington et al. Reference Terrington, Thompson and Hourigan2023). Then, using
$\theta _0 \propto \sqrt {G/d}$
, we obtain the following expressions for the gap-independent force and moment coefficients (Terrington et al. Reference Terrington, Thompson and Hourigan2023):



Here the terms
$\Delta C_{[i]}(\theta _0= 0)$
are evaluated using the polynomial fit to
$\Delta C_{[i]}$
. Finally, the total force and moment coefficients are approximated as the sum of a gap-dependent and a gap-independent contribution:



The gap-dependent part of the force and moment coefficient is obtained using the 2-D lubrication solution (2.36a – d ), while the gap-independent part is obtained using the zero-gap outer-flow numerical simulation, using (3.3).

Figure 14. Variation of the time-mean predicted (a) drag, (b) moment and (c) lift coefficients against
$G/d$
, for
$Re = 200$
,
$A=1$
and
$k = 1$
. Markers indicate numerical data obtained using finite-gap simulations. Dotted lines indicate the gap-dependent force and moment predictions obtained from lubrication theory (2.33), and dashed lines indicate the force and moment coefficients predicted using (3.4).
To validate the proposed decomposition of the force and moment coefficients, figure 14 presents the variation of (a)
$C_D$
, (b)
$C_M$
and (c)
$C_L$
against
$G/d$
for
$Re = 200$
,
$A = 1$
and
$k = 1$
. This figure presents a comparison between the force and moment coefficients obtained directly from finite-gap numerical simulations, and the predictions of (3.4a
–
c
).
Figures 14(a) and 14(b) show that (3.4a
) and (3.4c
) provide an excellent estimate of the drag and moment coefficients,
$C_{D}$
and
$C_M$
, for gap-diameter ratios less than
$10^{-1}$
. The gap-dependent drag and moment coefficients (
$C_{D,G}$
and
$C_{M,G}$
) differ from the numerically obtained
$C_D$
and
$C_M$
by a constant factor (
$C_{D,C}$
and
$C_{M,C}$
), which can be computed using (3.3) using the numerically obtained outer-flow solution.
Figure 14(c) demonstrates that the lift coefficient predicted using (3.3) and (3.4) is less accurate than the predicted drag and moment coefficients, consistent with the error estimate in (3.3). Specifically, (3.4) predicts a constant lift coefficient, however, the finite-gap simulations display a significant variation against
$G/d$
. The variation of
$C_L$
against
$\sqrt {G/d}$
in figure 14(c) appears to be approximately linear for small
$\sqrt {G/d}$
, and the numerical data appears to approach the predicted
$C_{L,C}$
as
$\sqrt {G/d}$
approaches
$0$
, consistent with the error estimate in (3.4b
).
Therefore, the decomposition of the force and moment coefficients presented in this work can accurately predict the drag and moment coefficients for a wide range of gap ratios. The lift coefficient estimated by our method instead represents an upper bound on the
$C_L$
that occurs in the limit
$G/d \rightarrow 0$
. Similar behaviour was observed in our previous study of the 2-D flow over an infinite circular cylinder (Terrington et al. Reference Terrington, Thompson and Hourigan2023).
3.4. Flow visualisations
We now present the results of the numerical computations performed over a range of
$A$
,
$Re$
and
$k$
. The aspect ratio is varied between
$A = 0.1$
and
$A = 10$
, the slip coefficient is varied between
$k = 0$
(no rotation) and
$k = 1$
(no slip) and the Reynolds number is varied between
$Re = 10$
and
$Re = 200$
. This subsection presents flow visualisations to discuss the effects of
$A$
,
$k$
and
$Re$
on the flow structures and wake shedding. Then, the force and moment coefficients are presented in § 3.5.

Figure 15. Particle visualisations for steady flow past a cylinder with aspect ratio
$A = 1$
, and (a)
$Re = 100$
,
$k = 1$
; (b)
$Re = 100$
,
$k = 0$
; (c)
$Re = 50$
,
$k = 1$
and (d)
$Re = 200$
,
$k = 0$
. Three-dimensional particle trajectories are projected onto the
$x{-}z$
(left) and
$x{-}y$
(right) planes, respectively, and selected particle tracks are shown in colour for emphasis.
Particle tracking is used for flow visualisations, as it provides clear images of the wake structures, particle paths and comparison with dye visualisations from experiments. The discrete particle method from ANSYS Fluent was used, with massless particles. For cases where the flow is steady, particles are injected along two lines: one upstream of the cylinder, at (
$x = -1$
,
$y = 0.5$
), and one inside the recirculation region, at (
$x = 1$
,
$y = 1$
). For unsteady cases, particles are injected above the centre of the cylinder, along the lines (
$x = 0$
,
$y = 1.1$
), (
$x = 0$
,
$y = 1.2$
) and (
$x = 0$
,
$y = 1.3$
). Along each of these lines, particle injection points are spaced a distance of
$0.05$
apart in the
$z$
direction.
We first consider cases where flow is steady. Figure 15 presents particle visualisations for steady flow past a cylinder with aspect ratio
$A = 1$
, and for four combinations of
$Re$
and
$k$
. The majority of particle streams are plotted in a semi-transparent green, however, four selected particle streams are plotted in different colours (blue, red, purple and yellow, respectively) for clarity. The same general flow structure is observed for all cases, namely, a semi-circular recirculation region located behind the cylinder. The head of this structure is visible as a single recirculation zone from the side views (right), while the legs of this structure are visible as a pair of symmetric recirculation bubbles from the top views (left). For example, the red particle trace in subfigure (a) is entrained into the recirculation region at the head of this structure, and exits via one of the legs, while the purple particle trace in subfigure (d) is entrained into the recirculation region at one of the legs, and exits at the head.
Figure 15 demonstrates that increasing
$Re$
for a fixed
$k$
results in an increase in the size of the recirculation region. There also appears to be more fluid entrained into the recirculation region (e.g. the red particle trace in subfigure a and more fluid ejected from the recirculation region (the blue particle trace exits the recirculation bubble after fewer loops in subfigure a than in subfigure c). Similarly, increasing
$k$
for a fixed
$Re$
corresponds to an increase in the size of the recirculation region, and greater entrainment of fluid into the recirculation bubble. Finally, we note that, for
$k = 1$
, fluid exiting the recirculation region remains at approximately
$y = 1$
, while for
$k = 0$
, fluid exiting the recirculation region approaches the wall.

Figure 16. Particle visualisations for steady flow past cylinders with aspect ratio (a,c)
$A = 10$
and (b,d)
$A = 2$
, at
$Re = 10$
and
$k = 1$
. Three-dimensional particle trajectories are projected onto the
$x{-}z$
(left) and
$x{-}y$
(right) planes, respectively, and selected particle tracks are shown in colour for emphasis.

Figure 17. Particle visualisations for unsteady flow past cylinders with aspect ratio
$A = 1$
, slip coefficient
$k = 1$
and Reynolds number (a)
$Re = 120$
, (b)
$Re = 130$
, (c)
$Re = 160$
and (d)
$Re = 200$
. Left and right columns show top and side views, respectively.
Figure 16 presents particle visualisations for steady flow at aspect ratios
$A = 10$
(figures 16
a and 16
c) and
$A = 2$
(figures 16
b and 16
d), for
$Re = 10$
and
$k = 1$
. The structure of the recirculation region is similar between the two cases, however, the recirculation bubble is much larger for
$A = 10$
than for
$A = 2$
. Therefore, the size of the recirculation bubble increases as
$A$
is increased.

Figure 18. Particle visualisations for unsteady flow past cylinders with aspect ratio
$A = 1$
, slip coefficient
$k = 0.5$
and Reynolds number (a)
$Re = 150$
, (b)
$Re = 170$
, (c)
$Re = 180$
and (d)
$Re = 200$
. Left and right columns show top and side views, respectively.
We now consider unsteady cases. Figure 17 presents particle visualisations for unsteady flow at
$A = 1$
,
$k = 1$
and for four different Reynolds numbers. A transient animation is also provided in supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.329. At
$A = 1$
and
$k = 1$
, transition to unsteady flow was found to occur between
$Re = 110$
and
$Re = 120$
. At
$Re = 120$
, there is an unsteady oscillation in the recirculation bubble, but no clear vortex shedding (figure 17
a). Shedding of hairpin vortices is observed at
$Re = 130$
(figure 17
b), and the strength and intensity of vortex shedding increases as
$Re$
is increased to
$Re = 160$
(figure 17
c) and
$Re = 200$
(figure 17
d). In all cases, the oscillations are not symmetric about
$z = 0$
, which we refer to as the asymmetric mode.
At
$A = 1$
and
$k = 1$
, Nanayakkara et al. (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
) observe asymmetric oscillations in the recirculation region as low at
$Re = 79$
, and first observe asymmetric vortex shedding at
$Re = 116$
. Their experiments are performed with freely rolling cylinders, which potentially experience flow-induced vibrations, which might explain the earlier transition to unsteady flow. Additionally, environmental disturbances may be present in the experiments, leading to an earlier flow transition.
We now consider a reduced rotation rate of
$k = 0.5$
. Figure 18 presents particle visualisations for unsteady flow at
$A = 1$
and
$k = 0.5$
, for a range of
$Re$
. A transient animation is also provided in supplementary movie 2. For these parameters, transition to unsteady flow was found to occur between
$Re = 140$
and
$Re = 150$
. At
$Re = 150$
(figure 18
a), there is an unsteady oscillation in the recirculation region, but no vortex shedding. At
$Re = 170$
(figure 18
b), hairpin vortices are shed from the cylinder. For these cases, the flow is symmetric about
$z = 0$
, which we refer to as the symmetric mode. Between
$Re = 170$
and
$Re = 180$
, flow transitions to the asymmetric shedding mode. This transition corresponds to a reduction in the flow unsteadiness (as indicated by the r.m.s. force and moment coefficients; see § 3.5), and no hairpin vortices are observed at
$Re = 180$
(figure 18
c). As the Reynolds number is increased to
$Re = 200$
(figure 18
d), asymmetric hairpin vortices are shed from the cylinder.
The critical Reynolds number for transition to unsteady flow increases as
$k$
is decreased. This behaviour has also been observed for the flow over infinite-span cylinders (Rao et al. Reference Rao, Stewart, Thompson, Leweke and Hourigan2011; Stewart et al. Reference Stewart, Thompson, Leweke and Hourigan2010b
). For
$A = 1$
, the critical
$Re$
is between
$110$
and
$120$
at
$k = 1$
and between
$140$
and
$150$
for
$k = 0.5$
. For
$k = 0$
, the critical
$Re$
is above 200, and flow remains steady over the entire range of
$Re$
considered in this study.

Figure 19. Particle visualisations for unsteady flow past cylinders with aspect ratio (a,c,e)
$A = 10$
and (b,d, f)
$A = 3$
, at (a,b)
$Re = 100$
,
$k = 1$
; (c,d)
$Re = 200$
,
$k = 1$
and (e, f)
$Re = 100$
,
$k = 0$
.
Figure 19 presents particle visualisations for both (a,c,e)
$A = 10$
and (b,d,f)
$A = 3$
, for three combinations of
$Re$
and
$k$
. Transient animations for several different
$A$
at these same combinations of
$Re$
and
$k$
are also provided in supplementary movies 3, 4 and 5. For
$Re = 100$
and
$k = 1$
(figures 19
a and 19
b; movie 3), flow is steady for aspect ratios up to
$A=1$
, but becomes unsteady between
$A = 1$
and
$A = 1.5$
. For aspect ratios between
$A = 1.5$
and
$A = 7.5$
, symmetric hairpin vortices are shed from the cylinder (figure 19
b). However, between
$A = 7.5$
and
$A = 10$
flow transitions to the asymmetric vortex shedding mode (figure 19
a). For
$Re = 200$
and
$k = 1$
(figures 19
c and 19
d; movie 5), transition to unsteady flow occurs between
$A = 0.5$
and
$A = 0.75$
. At this combination of
$Re$
and
$k$
, asymmetric vortex shedding is observed for all aspect ratios above this critical aspect ratio. Finally, for
$Re = 100$
and
$k = 0$
(figures 19
e and 19
f; movie 4), the transition to unsteady flow occurs between
$A = 2$
and
$A = 3$
. Symmetric vortex shedding is then observed for all aspect ratios above the critical aspect ratio, at least up to
$A = 10$
.

Figure 20. Comparison between particle visualisations obtained in the present work (top), to experimental dye visualisations presented in Nanayakkara et al. (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
) (bottom) (reproduced with permission). Physical parameters are (a)
$Re = 160$
,
$k = 1$
and
$A = 1$
for the numerical simulation, and
$Re = 157$
,
$k = 1$
and
$A = 1$
for the experiment, and (b)
$Re = 100$
,
$k = 1$
and
$A = 10$
for the numerical simulation, and
$Re = 104$
,
$k = 1$
and
$A = 10.66$
for the experiment.
Figure 20 presents a comparison between our numerical results and dye visualisations presented in Nanayakkara et al. (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
). Subfigure (a) presents a comparison at approximately
$Re = 160$
,
$A = 1$
and
$k = 1$
, while subfigure (b) presents a comparison at approximately
$A = 10$
,
$Re = 100$
and
$k = 1$
. In both cases, asymmetric shedding of hairpin vortices is observed. A fair qualitative agreement is observed between the numerical particle visualisations and the experimental dye visualisations.
Summarising the results of this section, the size of the recirculation bubble behind the cylinder increases as any of
$Re$
,
$k$
and
$A$
are increased. This also corresponds to a reduction in the stability of the flow, so that the critical Reynolds number for the transition from steady to unsteady flow decreases as both
$A$
and
$k$
are increased. This is consistent with the experimental study of Nanayakkara et al. (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
), who found that the critical Reynolds number for unsteady flow decreases as
$A$
is increased, and with previous numerical studies of the flow over an infinite cylinder (Stewart et al. Reference Stewart, Thompson, Leweke and Hourigan2010b
; Rao et al. Reference Rao, Stewart, Thompson, Leweke and Hourigan2011), where the critical Reynolds number for transition to unsteady flow decreases as
$k$
is increased.
Three distinct flow regimes were identified: steady flow, symmetric wake shedding and asymmetric wake shedding. For some combination of physical parameters, flow transitions directly from steady flow to the asymmetric regime, while for other parameters, flow first transitions to the symmetric mode, then a second transition to the asymmetric mode is observed. This is consistent with previous studies. For
$A = 0.4$
and
$k = 1$
, Pirozzoli et al. (Reference Pirozzoli, Orlandi and Bernardini2012) observed steady flow at
$Re = 300$
, symmetric hairpin vortices at
$Re = 400$
and asymmetric hairpins at
$Re = 500$
. Nanayakkara et al. (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
), however, were unable to observe the symmetric mode for
$A = 1$
and
$k = 1$
, in agreement with the present study for these same parameters.
3.5. Force and moment coefficients
This subsection presents the variation of the gap-independent force and moment coefficients against the parameters
$Re$
,
$A$
and
$k$
. The predictions obtained in this section can be combined with the gap-dependent force and moment predictions obtained using lubrication theory (2.36), to obtain the total force and moment coefficients, as functions of
$G/d$
,
$Re$
,
$A$
and
$k$
.
3.5.1. Effect of Reynolds number
Figure 21 presents the variation of the time-mean (
$\overline {C_{D,C}}$
,
$\overline {C_{M,C}}$
,
$\overline {C_{L,C}}$
) and
$\mathrm{r.m.s.}$
(
$C_{D,{rms}}$
,
$C_{M,{rms}}$
and
$C_{L,{rms}}$
) gap-independent force and moment coefficients against Reynolds number, for
$A = 1$
and for three slip coefficients,
$k = 0$
,
$0.5$
and
$1$
. Error bars in this figure represent the uncertainty due to the duration of sampling interval for unsteady statistics. To estimate the uncertainties, reduced samples of 250 time units were taken from the full sampling window of 500 time units. The estimated uncertainty was then obtained as the difference between the minimum and maximum of various statistics obtained using the reduced samples.
Figures 21(d), 21(e) and 21( f) present the
$\mathrm{r.m.s.}$
drag, moment and lift coefficients, respectively. The r.m.s. force and moment coefficients reflect the different flow regimes discussed in § 3.4. For
$k = 0$
, flow is steady for all
$Re$
considered in this study, and therefore, the
$\mathrm{r.m.s.}$
force and moment coefficients are zero. For
$k=0.5$
, flow is steady for Reynolds numbers up to
$Re = 140$
. The symmetrical regime is observed between
$Re = 150$
and
$Re = 170$
, and within this regime the
$\mathrm{r.m.s.}$
force and moment coefficients increase approximately linearly with
$Re$
. Transition to the asymmetric regime occurs between
$Re = 170$
and
$Re = 180$
, and this is accompanied by a decrease in the
$\mathrm{r.m.s.}$
force and moment coefficient to nearly zero. The r.m.s. force and moment coefficients then increase approximately linearly with
$Re$
within the asymmetric vortex shedding flow regime. Finally, for
$k = 1$
, there is a transition from the steady flow regime to the asymmetric vortex shedding regime between
$Re = 110$
and
$Re = 120$
. Above this critical Reynolds number, the r.m.s. force and moment coefficients increase monotonically with increasing
$Re$
.

Figure 21. Variation of the time-mean and r.m.s. gap-independent force and moment coefficients (a)
$\overline {C_{D,C}}$
, (b)
$\overline {C_{M,C}}$
, (c)
$\overline {C_{L,C}}$
, (d)
$C_{D,{rms}}$
, (e)
$C_{M,{rms}}$
and (f)
$C_{L,{rms}}$
, against
$Re$
, for
$A = 1$
and for various slip coefficients. Markers indicate the numerical data, with error bars indicating the uncertainty due to the duration of time sampling. Solid lines indicate the best fit regressions listed in table 6.
Figures 21(a) and 21(b) present the variation of the mean gap-dependent drag and moment coefficients (
$\overline {C_{D,C}}$
and
$\overline {C_{M,C}}$
), respectively, against Reynolds number. The drag coefficient is positive and decreases monotonically as
$Re$
is increased. The moment coefficient is negative and the magnitude of the moment coefficient also decreases monotonically as
$Re$
is increased. There is no obvious effect of the different wake shedding regimes (steady, symmetric and asymmetric) on the trends of both
$\overline {C_{D,C}}$
and
$\overline {C_{M,C}}$
against
$Re$
. Empirical fits were obtained for
$\overline {C_{D,C}}$
and
$\overline {C_{M,C}}$
of the form

and the coefficients are reported in table 6. Equation (3.5) scales as
$a_1/Re$
in the limit
$Re \rightarrow 0$
, which was chosen so that the gap-independent drag and moment are finite, but non-zero, in the Stokes-flow limit. Specifically, forces and moments for Stokes flow are typically non-dimensionalised as
$F = F^*/(\unicode{x03BC} U d) \propto C_D Re$
, where
$F^*$
is the force in dimensional units, and therefore, force and moment coefficients for Stokes flow are proportional to
$1/Re$
.
Table 6. Empirical fits for the dependence of the mean gap-independent force and moment coefficients against Reynolds number for various
$k$
and for
$A = 1$
.

Figure 21(c) presents the variation of the time-mean lift coefficient,
$\overline {C_{L,C}}$
, against
$Re$
. Within the steady flow regime, the lift coefficients decrease monotonically with increasing
$Re$
. Empirical fits are obtained of the form

with the coefficients reported in table 6. Unlike the drag and moment coefficients,
$\overline {C_{L,C}}$
approaches a finite value as
$Re \rightarrow 0$
. This implies that the lift force for Stokes flow will be zero, when normalised by Stokes-flow variables, as is expected.
Above the critical Reynolds number, the lift coefficient begins to increase. Empirical fits for the lift coefficient in this regime are given as

with coefficients reported in table 6.
To understand the physical mechanisms and flow structures leading to the variation of gap-independent drag and moment coefficients with respect to
$Re$
, we decompose the gap-independent force and moment coefficients into three components: pressure forces on the curved surface (
$C_{D,p}$
), viscous stresses on the curved surface (
$C_{D,\tau }$
,
$C_{M,\tau }$
) and viscous stresses on the cylinder end faces (
$C_{D,F}$
,
$C_{T,F}$
). The terms
$C_{D,F}$
and
$C_{M,F}$
are computed according to (2.11). The remaining terms are evaluated as in § 3.3, but considering only pressure or viscous terms, respectively.
Figure 22 presents the variation of each of these contributions to the gap-independent drag and moment coefficients against
$Re$
, for
$A = 1$
and
$k = 1$
. The viscous terms (
$C_{D,\tau }$
,
$C_{M,\tau }$
,
$C_{D,F}$
and
$C_{M,F}$
), including contributions from the cylinder end faces, decrease monotonically with increasing
$Re$
, due to the reduced viscosity with increasing
$Re$
. The pressure term, however, increases monotonically with increasing
$Re$
, and becomes the dominant source of gap-independent drag at higher
$Re$
.

Figure 22. Contributions to the gap-dependent (a) drag and (b) moment coefficients (
$C_{D,C}$
and
$C_{M,C}$
) from pressure forces on the cylinder’s curved surface (
$C_{D,p}$
), viscous forces on the cylinder’s curved surface (
$C_{D,\tau }$
and
$C_{M,\tau }$
) and viscous stresses on the cylinder end faces (
$C_{D,F}$
and
$C_{M,F}$
), for different
$Re$
, at
$A = 1$
and
$k = 1$
.
To further explain these trends, we compute the contributions to the gap-independent drag and moment coefficient from different angular positions on the cylinder’s curved surface. These are evaluated as



where
$f_{p,o}$
is the local contribution to
$C_{D,p}$
,
$f_{\tau, o}$
is the contribution to
$C_{D,\tau }$
and
$m_{\tau, o}$
is the contribution to
$C_{M,\tau }$
. These expressions were obtained by subtracting the integrands of (2.30) from the integrands of (2.15), treating the viscous and pressure terms separately.
Figure 23 plots the distributions of
$f_{p,o}$
,
$f_{\tau, o}$
and
$m_{\tau, o}$
against
$\theta$
for various values of
$Re$
, for
$k = 1$
and
$A = 1$
. Subfigure (b) plots the contribution of viscous forces to the gap-independent drag coefficient (
$f_{\tau, o}$
), while subfigure (c) plots the contributions of viscous forces to the moment coefficient (
$m_{\tau, o}$
). In both cases, the general features are approximately independent of
$Re$
:
$f_{\tau, o}$
is positive for all angles, with a large peak just upstream of the top of the cylinder, while
$m_{\tau, o}$
is negative, but otherwise exhibits the same trends. The peak shear stress corresponds to the boundary layer on the upper surface of the cylinder, before it separates to form a shear layer above the recirculation region. Increasing
$Re$
decreases the magnitude of the viscous stresses, due to reduced viscosity, but does not significantly alter the pattern of shear stresses on the cylinder surface. Hence, the viscous drag and moment coefficients (
$C_{D,C}$
and
$C_{\tau, C}$
) decrease monotonically with increasing
$Re$
, as shown in figure 22.

Figure 23. Contributions to the gap-independent (a,d) pressure drag, (b) viscous drag and (c) viscous torque from different angular positions on the cylinder’s curved surface for different
$Re$
, at (a,b,c)
$A = 1$
and
$k = 1$
, (d)
$A = 1$
and
$k = 0.5$
.
Figure 23(a) plots the sectional contribution from pressure forces to the gap-independent drag coefficient. The pressure forces are largely dominated by three distinct regions. A positive pressure upstream of the cylinder, negative pressure on the upstream surface of the cylinder and a negative pressure downstream of the cylinder. These three contributions also occur for a rotating
$A=1$
cylinder in a uniform free-stream flow (Yang et al. Reference Yang, Wang, Guo and Zhang2023). Increasing
$Re$
increases the magnitude of positive pressure upstream of the cylinder, due to increased inertia of the fluid. The area of the negative pressure region behind the cylinder also increases, due to the increased size of the recirculation region. The negative contribution from negative pressure on the cylinder’s front face also decreases as
$Re$
increases. All of these effects produce the increasing pressure drag (
$C_{D,p}$
) as
$Re$
is increased.
There are only small changes to
$f_{p,o}$
between
$Re = 100$
and
$Re = 200$
, despite substantial changes in
$Re$
, and a change from steady to unsteady flow. Figure 23(d) plots
$f_{p,o}$
against
$\theta$
for
$A = 1$
and
$k = 0.5$
, for
$Re = 140$
(steady flow),
$Re = 150$
and
$Re = 170$
(symmetric shedding) and
$Re = 180$
and
$Re = 200$
(asymmetric shedding). The change from steady flow to symmetric vortex shedding produces no substantial change to the pressure distribution behind the cylinder. The change from symmetric to asymmetric shedding produces a small increase in the negative pressure behind the cylinder between
$\theta = \pi /2$
and
$\theta = \pi$
, which produces a small increase in the pressure drag.
3.5.2. Effect of slip coefficient

Figure 24. Variation of the time-mean and r.m.s. gap-independent force and moment coefficients (a)
$\overline {C_{D,C}}$
, (b)
$\overline {C_{M,C}}$
, (c)
$\overline {C_{L,C}}$
, (d)
$C_{D,{rms}}$
, (e)
$C_{M,{rms}}$
and (f)
$C_{L,{rms}}$
, against
$k$
, for
$A = 1$
and three different Reynolds numbers,
$Re = 10$
,
$Re=100$
and
$Re = 200$
. Markers indicate the numerical data, with error bars indicating the uncertainty due to the duration of time sampling. Solid lines indicate the empirical fits to the data listed in table 7.
Figure 24 presents the variation of the time-mean (
$\overline {C_{D,C}}$
,
$\overline {C_{M,C}}$
,
$\overline {C_{L,C}}$
) and r.m.s. (
$C_{D,{rms}}$
,
$C_{M,{rms}}$
$C_{L,{rms}}$
) gap-independent force and moment coefficients against slip coefficient, for
$A = 1$
and three different Reynolds numbers,
$Re = 10$
,
$Re = 100$
and
$Re = 200$
. Error bars indicate the estimated uncertainty due to the duration of sampling for unsteady statistics.
Table 7. Empirical fits for the dependence of the mean gap-independent force and moment coefficients against slip coefficient, for various
$Re$
and for
$A = 1$
.

Figures 24 (d), 24(e) and 24(f) present the r.m.s. force and moment coefficients. For
$Re = 10$
and
$Re = 100$
, flow is steady for all values of
$k$
considered, and the r.m.s. drag, lift and moment are all zero. For
$Re = 200$
, flow is unsteady for
$k\gtrapprox 0.3$
, and asymmetric vortex shedding occurs for all
$k$
above this value. The r.m.s. force and moment coefficients increase approximately linearly with
$k$
within this flow regime.
Figures 24(a) and 24(b) present the variation of the time-averaged gap-independent drag and moment coefficients, respectively, against slip coefficient. Both
$\overline {C_{D,C}}$
and
$\overline {C_{M,C}}$
vary linearly with
$k$
, and empirical fits of the form

were obtained, with coefficients provided in table 7. The drag coefficient is positive, while the moment coefficient is negative, and the magnitudes of both
$\overline {C_{D,C}}$
and
$\overline {C_{M,C}}$
increase as
$k$
is increased.
The linear variation of both
$\overline {C_{D,C}}$
and
$\overline {C_{M,C}}$
against
$k$
suggests that the drag and moment coefficients for any combination of
$k$
and
$Re$
can be obtained by linearly interpolating between the predictions of (3.5) for
$k = 0$
and
$k = 1$
, at least for
$A = 1$
,

In the above expression,
$\overline {C_{[D,M],C}}_{k=0}$
and
$\overline {C_{[D,M],C}}_{k=1}$
are the empirical fits for the force and moment coefficients against
$Re$
(3.5) obtained for
$k = 0$
and
$k = 1$
, respectively.
Figure 24(c) presents the variation of
$\overline {C_{L,C}}$
against
$k$
. The lift coefficient depends approximately quadratically against
$k$
, albeit with some scatter in the data. Empirical fits of the form

were obtained, with coefficients provided in table 7. The minimum lift coefficient occurs at approximately
$k = 0.8$
for all
$Re$
considered, although there are not enough data points to determine precisely the location of minimum
$\overline {C_{L,C}}$
.
Figure 25 plots the contributions to the gap-independent drag and moment coefficients (
$C_{D,C}$
and
$C_{M,C}$
) from pressure forces on the cylinder’s curved surface (
$C_{D,p}$
), viscous forces on the cylinder’s curved surface (
$C_{D,\tau }$
and
$C_{M,\tau }$
) and viscous stresses on the cylinder end faces (
$C_{D,F}$
and
$C_{M,F}$
) against
$k$
, for
$Re = 200$
and
$A = 1$
. At this
$Re$
, the drag coefficient is dominated by pressure forces, and the contributions from the cylinder end faces to the drag and moment are negligible. Both the viscous and pressure forces increase as
$k$
is increased. The increase is approximately linear, except between
$k = 0.2$
and
$k = 0.4$
, where there is a transition from steady flow to asymmetric vortex shedding, associated with an increase in the pressure drag.

Figure 25. Contributions to the gap-dependent (a) drag and (b) moment coefficients (
$C_{D,C}$
and
$C_{M,C}$
) from pressure forces on the cylinder’s curved surface (
$C_{D,p}$
), viscous forces on the cylinder’s curved surface (
$C_{D,\tau }$
and
$C_{M,\tau }$
) and viscous stresses on the cylinder end faces (
$C_{D,F}$
and
$C_{M,F}$
) against
$k$
, at
$A = 1$
and
$Re = 200$
.
Figure 26(a) plots the contributions to pressure drag (
$f_{p,o}$
) against
$\theta$
for
$A = 1$
,
$Re = 200$
and for various values of
$k$
. Increasing rotation rate produces a decrease in both positive and negative pressures from the front of the cylinder, and an increase in the contribution from negative pressure on the cylinder’s downstream face, due to an increasing size of the recirculation region. These general trends are in agreement with the rotating
$A = 1$
cylinder in a uniform free-stream flow (Yang et al. Reference Yang, Wang, Guo and Zhang2023). The transition from steady flow to asymmetric shedding between
$k = 0.2$
and
$k = 0.4$
results in a substantial increase in the negative pressure behind the cylinder, resulting in a jump in the pressure drag observed in figure 25.

Figure 26. Contributions to the gap-independent (a) pressure and (b) viscous drag from different angular positions on the cylinder’s curved surface for different
$k$
, at
$A = 1$
and
$Re = 200$
.
Figure 26(b) plots the viscous contribution to the drag coefficient for various values of
$k$
. The magnitude of the large peak (corresponding to the attached boundary layer on the cylinder’s front surface) is not significantly affected by the rotation rate. However, an increased rotation rate increases the viscous stresses on the remaining regions of the cylinder, resulting in an increased viscous drag coefficient. Similar trends occur for the moment coefficient (not shown for brevity).
3.5.3. Effect of aspect ratio

Figure 27. Variation of the time-mean and r.m.s. gap-independent force and moment coefficients (a)
$\overline {C_{D,C}}$
, (b)
$\overline {C_{M,C}}$
, (c)
$\overline {C_{L,C}}$
, (d)
$C_{D,{rms}}$
, (e)
$C_{M,{rms}}$
and (f)
$C_{L,{rms}}$
against
$A$
, for four combinations of
$Re$
and
$k$
. Markers indicate the numerical data, with error bars indicating the uncertainty due to the duration of time sampling. In subfigures (a)–(c), solid lines represent a cubic-spline interpolation of the numerical data. In subfigures (d)–(f), solid lines indicate the symmetric mode, while dashed lines indicate the asymmetric mode, with linear interpolation used between data points.
Figure 27 presents the variation in time-averaged and r.m.s. gap-independent force and moment coefficients against aspect ratio for four different combinations of
$Re$
and
$k$
. The first three cases,
$k = 1$
and
$Re = 10$
,
$k = 1$
and
$Re = 100$
and
$k = 1$
and
$Re = 200$
, were selected to demonstrate the effect of Reynolds number for different aspect ratios, while the fourth case,
$k = 0$
and
$Re = 100$
, was chosen to demonstrate the influence of slip coefficient for various aspect ratios. Finally, error bars indicate the estimated uncertainty due to the duration of sampling for unsteady statistics.
Figures 27(d), (e) and (f) present the r.m.s. drag, moment and lift coefficients, respectively. For
$Re = 10$
and
$k = 1$
, the flow is steady for all aspect ratios, and the r.m.s. force and moment coefficients are zero. For
$Re = 100$
and
$k = 0$
, the flow is steady up to
$A = 2$
and then transitions to the symmetric vortex shedding regime. In this regime, the r.m.s. force and moment coefficients are maximum at
$A = 3$
and then decrease as
$A$
is increased. For
$Re = 100$
and
$k = 1$
, the flow is steady up to
$A = 1$
, but transitions to the symmetric vortex shedding regime by
$A = 1.5$
. Within this regime, the r.m.s. force and moment coefficients increase to a maximum value at
$A = 3$
, then generally decrease as
$A$
is increased further. At
$A= 10$
, flow transitions to the asymmetric regime. This is accompanied by a significant increase in the error bars indicating the uncertainty due to the length of time sampling. Referring to the time histories of drag coefficient presented in figure 28, the asymmetric mode (subfigures c–f) are significantly less periodic than the symmetric modes (subfigures a and b), resulting in larger uncertainties on the unsteady statistics. Finally, for
$Re = 200$
and
$k = 1$
, flow is steady up to
$A = 0.5$
and then transitions directly to the asymmetric mode. Due to the large uncertainties, no clear trend is observed in the r.m.s. force and moment coefficients for this case.

Figure 28. Time histories of the gap-independent drag coefficient for a selection of cases.
Figure 27(a) presents the variation of time-mean drag coefficient against
$A$
. For small
$A$
, the drag coefficient increases significantly as
$A$
approaches zero. For
$Re = 10$
and
$k = 1$
, the mean drag coefficient decreases monotonically. The remaining cases reach a minimum value of
$\overline {C_{D,C}}$
at
$A = 1.5$
for
$Re = 100$
and
$k = 0$
, and at
$A = 2$
for both
$Re = 200$
and
$k = 1$
, and
$Re = 100$
and
$k = 1$
. The mean drag coefficient then increases slightly as
$A$
is increased above this value. For
$Re = 100$
and
$k = 1$
, the change of shedding mode between
$A = 7.5$
and
$A=10$
results in a decrease in
$\overline {C_{D,C}}$
. The dashed line in figure 27(a) indicates the value of
$C_{D,C}$
obtained for steady 2-D flow over an infinite-span cylinder (Terrington et al. Reference Terrington, Thompson and Hourigan2023), and the numerically obtained values of
$\overline {C_{D,C}}$
approximately approach this value as
$A$
is increased. Of course, some differences between the limiting value of
$\overline {C_{D,C}}$
for a finite-span cylinder and that obtained for steady 2-D flow are expected due to 3-D and unsteady effects.
Figure 27(b) presents the variation of time-mean gap-dependent moment coefficient against
$A$
. The moment coefficient is negative, and the magnitude of the moment coefficient increases significantly as
$A$
approaches zero. The magnitude of
$\overline {C_{M,C}}$
then decreases monotonically as
$A$
is increased, approximately approaching the value obtained for steady 2-D flow (dashed lines).
Finally, figure 27(c) presents the variation in the mean gap-independent lift coefficient against
$A$
. The lift coefficient increases monotonically as
$A$
is increased, approximately approaching the value obtained for steady 2-D flow (dashed lines) in the limit
$A\rightarrow \infty$
. For all cases with
$k = 1$
, the lift coefficient becomes negative when
$A$
is small. For
$k = 0$
and
$Re = 100$
, the lift coefficient remains positive, but close to zero, at least for
$A = 0.1$
.
Figure 29 plots the contributions to the gap-independent drag and moment coefficients (
$C_{D,C}$
and
$C_{M,C}$
) from pressure forces on the cylinder’s curved surface (
$C_{D,p}$
), viscous forces on the cylinder’s curved surface (
$C_{D,\tau }$
and
$C_{M,\tau }$
) and viscous stresses on the cylinder end faces (
$C_{D,F}$
and
$C_{M,F}$
) against
$A$
, for
$Re = 200$
and
$k = 1$
. The contributions from the cylinder end faces is negligible for large
$A$
, but rises sharply as
$A \rightarrow 0$
, due to the
$1/A$
factor in (2.11). Physically, this represents a larger drag from the end faces, relative to the cylinder width. The viscous drag force increases monotonically with
$A$
, while the viscous contribution to the moment decreases monotonically with increasing
$A$
. Finally, the pressure drag is minimum at
$A \approx 2$
and increases for either small
$A$
or large
$A$
.

Figure 29. Contributions to the gap-dependent (a) drag and (b) moment coefficients (
$C_{D,C}$
and
$C_{M,C}$
) from pressure forces on the cylinder’s curved surface (
$C_{D,p}$
), viscous forces on the cylinder’s curved surface (
$C_{D,\tau }$
and
$C_{M,\tau }$
) and viscous stresses on the cylinder end faces (
$C_{D,F}$
and
$C_{M,F}$
) against
$A$
, at
$k = 1$
and
$Re = 200$
.
Figure 30 plots the local contributions to the gap-independent drag and moment coefficients from viscous and pressure forces against angular position on the cylinder’s curved surface. The magnitude of viscous stress increases as
$A$
is reduced, however, the pattern of viscous stresses on the cylinder’s surface is not significantly affected. To explain this, consider that at the ends of the cylinder, momentum may be transported into the boundary layer in the spanwise direction, resulting in larger shear stresses near the ends of the cylinder. The relative contribution from this effect increases as
$A$
decreases. The contribution to the moment coefficient (subfigure c) is positive over the entire span of the cylinder, and therefore, the moment coefficient increases as
$A$
is decreased. However, the contribution from viscous forces to the drag force (subfigure b) includes regions of both positive and negative wall-shear stress, so the total viscous drag decreases as
$A$
is decreased.

Figure 30. Contributions to the gap-independent (a,d) pressure drag, (b) viscous drag and (c) viscous torque from different angular positions on the cylinder’s curved surface for different
$A$
, at (a,b,c)
$Re = 200$
and
$k = 1$
, (d)
$Re = 100$
and
$k = 1$
.
Figure 30(a) plots the local contributions of the pressure force to the gap-independent drag. Increasing
$A$
results in an increased force from the negative pressure region behind the cylinder, due to the increased size of the recirculation region. The pressure force on the front face is more complex. For small
$A$
, there is a positive pressure peak at
$\theta \approx 3\pi /2$
, which increases as
$A$
decreases, and drives the increase in pressure force for small
$A$
. However, for large
$A$
, there is a second pressure peak at
$\theta \approx 7\pi /4$
, which increases as
$A$
increases. This second peak, along with the increased contribution from the recirculation region behind the cylinder, produces the increasing pressure drag with increasing
$A$
. This second peak can be explained by considering that proportionally less fluid can flow around the sides of the cylinder as
$A$
increases (i.e. an increased blockage effect) resulting in a greater pressure to resist the fluid’s inertia.

Figure 31. Contours of time-averaged pressure, with the asymptotic pressure singularity removed (
$p - p_o$
), on the cylinder’s curved surface at
$Re = 100$
and
$k = 1$
, for (a)
$A = 1$
, (b)
$A = 2$
, (c)
$A = 3$
, (d)
$A = 5$
, (e)
$A = 7.5$
and (f)
$A = 10$
.
Finally, we consider the effects of changing wake modes on the drag. Figure 30(d) plots the contributions to pressure drag for both
$A = 7.5$
(symmetric mode) and
$A = 10$
(asymmetric mode), for
$Re = 100$
and
$k = 1$
. While the positive pressure peak upstream of the cylinder increases, consistent with the trends seen in figure 30(a), the downstream negative pressure peak decreases, despite the increase in aspect ratio. The loss of symmetry due to asymmetric vortex shedding evidently weakens the strength of the recirculation region for large
$A$
, resulting in a decrease in pressure drag.
Figure 31 plots contours of time-averaged pressure, with the asymptotic pressure singularity subtracted (
$\bar {p} - {p}_o$
), on the cylinder’s surface, for
$Re = 100$
and
$k = 1$
, to illustrate the spanwise contributions to the pressure force. The positive pressure peak upstream of the cylinder is located at the centre of the cylinder, consistent with the explanation of a blockage effect. The negative pressure behind the cylinder is strongest near the edges, particularly for large aspect ratios. This indicates flow separation from the ends of the cylinder provides the strongest contribution to the negative pressure behind the cylinder. The visualisations in figure 19 show that fluid flowing around the ends of the cylinder is entrained into the recirculation region, and the size of these recirculation cells increases as
$A$
is increased. For the symmetric shedding at
$A = 7.5$
, there is also an increased negative pressure at the centre of the cylinder (i.e. the symmetry plane). However, for asymmetric shedding at
$A = 10$
, this negative pressure region is lost, possibly a result of the loss of symmetry.
3.6. Comparison to 2-D and free-stream cylinders
Figure 32 presents a comparison between the gap-independent drag coefficient obtained in the present study, the gap-independent drag coefficient for 2-D flow over an infinite cylinder (Terrington et al. Reference Terrington, Thompson and Hourigan2023), the drag coefficient for a finite-span cylinder in a uniform free-stream flow (Yang et al. Reference Yang, Feng and Zhang2022, Reference Yang, Wang, Guo and Zhang2023), and the drag coefficient for the 2-D cylinder in the free stream (Qu et al. Reference Qu, Norberg, Davidson, Peng and Wang2013; Mittal & Kumar Reference Mittal and Kumar2003). Interestingly, there is a surprisingly close agreement between the gap-independent drag coefficient for the finite-span cylinder near a wall and the finite-span cylinder in a free-stream flow, despite the substantial influence of the wall on the wake.

Figure 32. Comparison between the gap-independent drag coefficient (
$C_{D,C}$
) for the finite-span cylinder near a wall and for the 2-D flow over an infinite cylinder near a wall (Terrington et al. Reference Terrington, Thompson and Hourigan2023); and the drag coefficient (
$C_{D}$
) for either a finite-span (Yang et al. Reference Yang, Feng and Zhang2022, Reference Yang, Wang, Guo and Zhang2023) or infinite-span (Mittal & Kumar Reference Mittal and Kumar2003; Qu et al. Reference Qu, Norberg, Davidson, Peng and Wang2013) circular cylinder in a uniform free-stream flow. Physical parameters are (a)
$A = 1$
and
$k = 0$
, (b)
$A = 1$
and
$Re = 200$
and (c)
$Re = 100$
and
$k = 0$
.
Figure 32(a) plots the variation in drag coefficient against
$Re$
for
$A = 1$
and
$k = 0$
. For all cases, there is a general trend of
$C_D$
decreasing as
$A$
is increased. However, the finite-span cylinders have a far greater variation compared with the infinite (2-D) cylinders, and exhibit a substantial increase in
$C_D$
as
$Re$
decreases. Figure 32(b) plots the variation in drag coefficients against
$k$
for
$A= 1$
and
$Re = 200$
. For both the 2-D and 3-D rolling cylinders, as well as the finite-span cylinder, the drag coefficient increases as rotation rate increases. This is not the case for the free-stream 2-D cylinder, for which drag decreases as
$k$
is increased. Finally, figure 32(c) plots the variation in drag coefficient against
$A$
for
$k = 0$
and
$Re = 100$
. The free-stream finite-span cylinder displays similar trends as the near-wall cylinder, i.e. a sharp rise in
$C_D$
as
$A$
decreases, which is likely attributed to the same general features: a rise in the pressure force and an increased contribution from viscous stresses on the cylinder side faces. For the wall-adjacent cylinder, there is a gradual increase in
$C_{D,C}$
towards the value for an infinite (2-D) cylinder as
$A$
is increased. This is not observed for the free-stream cylinder, likely because the maximum aspect ratio
$A = 2$
is too small to observe this trend.
3.7. Total drag and moment coefficients
Figure 33 plots the variation of total drag and moment coefficients (the sum of gap-dependent and gap-independent coefficients) against (a)
$Re$
, (b)
$k$
and (c)
$A$
, for three different values of
$G/d$
. Dashed and dash-dotted lines indicate the gap-dependent drag and moment coefficients, respectively. Decreasing
$G/d$
results in larger contributions from the gap-dependent drag and moment coefficients, so that the total drag increases substantially as
$G/d$
is decreased.

Figure 33. Variation of the total drag and moment coefficients against (a)
$Re$
, (b)
$k$
and (c)
$A$
, for three different values of
$G/D$
. Dashed and dash-dotted lines indicate the contribution from gap-dependent drag and moment coefficients, respectively. Unless otherwise stated, physical parameters are
$A = 1$
,
$k= 1$
and
$Re = 100$
.
Figure 33(a) plots the variation of drag and moment coefficients against
$Re$
. Both the gap-dependent and gap-independent drag and moment coefficients decrease monotonically with increasing
$Re$
, so this trend is also observed in the total drag and moment coefficients. Figure 33(b) plots the variation of drag and moment coefficients against
$k$
. The gap-dependent drag coefficient is constant with respect to
$k$
, while the gap-independent drag coefficient has a moderate linear increase as
$k$
increases. As a result, the total drag coefficient increases only weakly as
$k$
is increased. Both gap-dependent and gap-independent moment coefficients, however, increase linearly as
$k$
is increased, and the same trend is observed in the total moment coefficients.
Figure 33(c) plots the variation of drag and moment coefficients against
$A$
. For small
$A$
, the gap-dependent force and moment coefficients increase with increasing
$A$
, while the gap-independent terms decrease with increasing
$A$
. For
$G/d=10^{-2}$
, the increase in the gap-dependent term is less than the decrease in the gap-independent term, and the total drag increases as
$A$
is increased. For
$G/d = 10^{-3}$
, the changes in both gap-dependent and gap-independent terms are approximately balanced, and there is only a relatively mild variation in the total drag and moment coefficients. For
$G/d = 10^{-4}$
, the variation in the gap-dependent term dominates, and the total drag increases as
$A$
is increased.
For the range of parameters considered here, both gap-dependent and gap-independent contributions are significant. For
$A = 1$
and
$k = 1$
, the gap-independent term contributes
$5.1\%$
of the total drag at
$Re = 20$
and
$G/d = 10^{-4}$
,
$17.5\%$
at
$Re = 200$
and
$G/d = 10^{-4}$
,
$50\%$
at
$Re = 20$
and
$G/d = 10^{-2}$
, and
$63\%$
at
$Re = 200$
and
$G/d = 10^{-2}$
. For a sufficiently small
$G/d$
, the gap-dependent term will dominate over the gap-independent term. At
$Re = 200$
,
$A = 1$
and
$k = 1$
, the gap-independent term is less than
$1\%$
of the total drag for
$G/d \lt 2.5\times 10^{-7}$
. Likewise, the gap-independent term will dominate for a sufficiently large
$Re$
, assuming a fixed
$G/d$
. For
$G/d = 10^{-4}$
, and assuming that
$C_{D,C} \approx 1$
at large
$Re$
, the gap-dependent drag accounts for less than
$1\%$
of the total drag for
$Re \gt 1.26\times 10^5$
. For
$G/d = 10^{-2}$
, however, the gap-dependent contribution becomes negligible (below
$1\%$
) at
$Re \gt 1.26\times 10^{4}$
.
3.8. Application to experimental measurements of rolling cylinders
We now apply the results of our study to experimental measurements of finite-span cylinders rolling on a wall presented by Nanayakkara et al. (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
). Following the discussion in Nanayakkara et al. (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
), we assume that the cylinder contacts the wall via surface asperities. The corresponding hydrodynamic drag is estimated by assuming an effective gap approximately equal to the peak surface roughness. The total drag force also includes the contact force necessary to ensure the cylinder rolls without slip (
$k = 1$
), so that the effective drag coefficient is (Nanayakkara et al. Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
)

Figure 34 presents a comparison between Nanayakkara et al.’s (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
) experimental measurements and the theoretical predictions, with an effective gap determined by fitting to the experimental data. Subfigure (a) plots the variation in drag coefficient against
$Re$
for two different non-dimensional surface roughnesses (
$\xi = 5.9\times 10^{-4}$
and
$\xi = 2.0\times 10^{-3}$
), where
$\xi = R_p/d$
is the dimensionless peak roughness and
$R_p$
is the combined peak roughness of both the cylinder and panel. The solid line represents the theoretical predictions of (3.14) using the numerical and analytical results obtained in the present study, for values of
$G/d$
tuned to fit the experimental measurements. Good agreement between the numerical and experimental results is obtained, with an effective gap of the same order of magnitude as the peak roughness.

Figure 34. Comparison between Nanayakkara et al.’s (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a
) experimental measurements and the combined analytical/numerical predictions using (3.14) for the variation of effective drag coefficient against (a)
$Re$
and (b)
$G/d$
. The cylinder aspect ratio is
$A = 1$
.
Figure 34(b) plots the variation of effective drag coefficient against
$G/d$
for two Reynolds numbers (
$Re = 100$
and
$Re = 150$
). Experimental data for cylinders of various surface roughness are also provided, assuming that
$G/d = 0.58 \xi$
for
$Re = 100$
and
$G/d = 0.64\xi$
for
$Re = 150$
. The scaling factors for
$G/d$
and
$\xi$
were determined by fitting the experimental data. The results suggest that the drag applied to a rough cylinder rolling on a wall can be predicted by assuming a smooth cylinder and wall, with a gap equal to approximately
$0.6\xi$
.
4. Conclusions
We have computed the force and moment coefficients for a finite-span circular cylinder immersed in a fluid that travels parallel to a plane wall with a fixed rotation rate, while maintaining a fixed gap between the cylinder and the wall. The force and moment coefficients are a function of four dimensionless variables: the Reynolds number
$Re = \rho U d/\mu$
, the gap-diameter ratio
$G/d$
, the aspect ratio
$A = W/d$
and the slip coefficient
$k = \Omega d/(2U)$
. Here,
$d$
is the cylinder diameter,
$W$
is the cylinder span,
$U$
and
$\Omega$
are the translational and rotational velocities of the cylinder,
$\mu$
and
$\rho$
are the dynamic viscosity and density of the fluid, respectively, and
$G$
is the gap between the cylinder and the wall.
Using the method of matched asymptotic expansions, the total force and moment coefficients were expressed as the sum of gap-dependent and gap-independent terms. The gap-dependent force and moment coefficients were obtained by using lubrication theory to solve for the flow in the narrow gap between the cylinder and the wall. We expressed the gap-dependent force and moment coefficients analytically in terms of a nonlinear function
$f(\hat {A})$
, where
$\hat {A} = A/\sqrt {G/d}$
is the modified aspect ratio. The 2-D Reynolds equation was solved numerically for various values of
$\hat {A}$
to evaluate this function numerically.
The gap-independent force and moment coefficients were obtained by performing numerical simulations of the outer flow for various values of the parameters
$Re$
,
$A$
and
$k$
. We find that the combined predictions of lubrication theory for the inner flow, and numerical simulations of the outer flow, accurately predict the variation of the drag and moment coefficients against
$G/d$
. Our method does not predict the variation of the lift coefficient against
$G/d$
, but does predict the upper bound of
$C_L$
obtained in the limit
$G/d \rightarrow 0$
.
The critical
$Re$
for transition to unsteady flow was found to decrease as either
$A$
or
$k$
was increased. For
$A = 1$
, the gap-independent drag and moment coefficients both decrease in magnitude as
$Re$
increases and increase linearly in magnitude as
$k$
increases, and empirical fits for the force and moment coefficients were obtained. As
$A$
was increased, the gap-independent moment coefficient decreased in magnitude, while the gap-independent lift coefficient increases. The gap-independent drag coefficient is generally minimum near
$A = 2$
and increased as
$A$
was either increased or decreased from this value.
This work allows the drag coefficient for a finite-span rolling circular cylinder to be related to an effective hydrodynamic gap between the cylinder and the wall. Nanayakkara et al. (Reference Nanayakkara, Zhao, Terrington, Thompson and Hourigan2024a ) demonstrate experimentally that this effective gap is approximately equal to the sum of peak roughness of the cylinder and the wall. However, in their analysis, they assumed that the gap-independent force and moment coefficients were equal to the corresponding values for steady 2-D flow, since values for finite-span cylinders have not been obtained previously. The present study allows a more accurate prediction of the drag coefficient under the effective gap model, since the effects of finite aspect ratio on the outer-flow force and moment coefficients are included in the model.
By separating the flow into inner and outer regions, full 3-D numerical simulations are only required to obtain the outer-flow solution, which does not depend on
$G/d$
. This reduces the parameter space to be explored by numerical simulations from four to three, significantly reducing the computational effort required to predict the drag and moment applied to the cylinder. This method may also be applicable to other rolling bodies, such as spheres, ellipsoids or elliptical cylinders, which can be explored in future work.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.329.
Funding
This work was supported by the Australian Government through the Australian Research Council’s Discovery Projects funding scheme (projects DP200100704 and DP210100990), and by computational resources provided by the National Computational Infrastructure and Pawsey Supercomputer Centre (Merit Grant d71) under the National Computational Merit Allocation Scheme.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical computation of the lubrication flow
This appendix summarises the numerical methodology used to obtain solutions to (2.9). First, the infinite domain is transformed into a finite domain using a coordinate transformation. Then, the transformed version of (2.9) is discretised using the finite-difference method. Finally, the resulting matrix equation is solved using MATLAB’s mldivide (matrix left division) function.
We first transform the infinite domain
$-\infty \lt \hat {x} \lt \infty$
to a finite domain
$-1\lt \hat {t}\lt 1$
, using the following transformation:

The parameter
$k_t$
is a stretching parameter that controls the distribution of grid points in the physical domain. We have found
$k_t = 0.1$
to give suitable results.
For large
$\hat {A}$
,
$\hat {z}$
derivatives of
$\hat {g}$
are only significant in a thin end correction layer. Therefore, a second stretching function was used to concentrate a large number of grid points inside the end correction layer. We first estimate the thickness of the end correction layer. Assuming that
$\hat {z}$
derivatives are much larger than
$\hat {x}$
derivatives in this region, (2.9) becomes

which has solutions of the form
$\hat {g} \approx C \exp ( \pm \sqrt {6/H} \hat {z})$
Therefore, the width of the end correction region will approximately be proportional to
$\sqrt {H/6}$
, which has a minimum value of
$1/\sqrt {6}$
at
$\hat {x} = 0$
.
We use the following stretching function to transform the physical coordinate
$\hat {z}$
to a computational coordinate
$\hat {s}$
(Vafakos & Papadopoulos Reference Vafakos and Papadopoulos2020):

The parameters
$P$
and
$Q$
determine the distribution of grid points in the physical domain. Vafakos & Papadopoulos (Reference Vafakos and Papadopoulos2020) provide an algorithm to determine the values of
$P$
and
$Q$
required to place a particular percentage of grid points within a thin boundary layer. We use this algorithm to place
$25 \%$
of the grid points within the end correction layers (i.e. a distance of
$1/\sqrt {6}$
from the ends of the cylinder).
For
$\hat {A}\lt 8/\sqrt {6}$
, the end correction layers cover more than
$25\%$
of the domain. Therefore, we use a uniform transformation,

rather than the stretching function when
$\hat {A}\lt 8/\sqrt {6}$
.
Using the transformations A1, A3 and A4, the physical domain
$(\hat {x},\hat {z})$
is transformed to the computational domain
$(\hat {t},\hat {s})$
, where
$-1\lt \hat {t}\lt 1$
and
$-1\lt \hat {z}\lt 1$
. Furthermore,
$\hat {g}$
is symmetric about both
$\hat {s} = 0$
and
$\hat {t} = 0$
. Therefore, we consider the reduced domain
$0 \leqslant \hat {t} \leqslant 1$
and
$0 \leqslant \hat {s} \leqslant 1$
. Boundary conditions are
$\partial \hat {g}/\partial \hat {s} = 0$
at
$\hat s = 0$
,
$\partial \hat {g}/\partial \hat {t} = 0$
at
$\hat t = 0$
,
$\hat {g} = 0$
at
$\hat {t} = 1$
and
$\hat {g} = 0$
at
$\hat {s} = 1$
.
Finally, in transformed coordinates, (2.9) is expressed as

To ensure that each of the the terms remains finite, we also multiply the entire equation by
$\hat {x}/(1+\hat {x})$
,

The computational domain (
$\hat {t},\hat {s}$
) was covered with a uniform grid comprising
$n_t$
and
$n_s$
points in the
$\hat {t}$
and
$\hat {s}$
directions, respectively. Derivatives of
$\hat {g}$
in (A6) were discretised using second-order central finite differencing. The resulting matrix equation was solved using MATLAB’s mldivide function. The finite-difference matrix is a sparse, banded matrix, and therefore, MATLAB automatically selects the banded solver to perform the matrix inversion.
Appendix B. Asymptotic behaviour of
$f$
In this appendix we compute the asymptotic behaviour of the function
$f(\hat {A})$
, in the limits of both large and small
$\hat {A}$
. We first note that Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) have computed estimates for the asymptotic behaviour of
$\hat {p}$
in both of these limits.
We begin with the limit of small
$\hat {A}$
. From Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) equation (2.27), we have

and therefore,


Finally, we obtain

We now consider the limit of large
$\hat {A}$
. In this limit, Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) provide asymptotic solutions valid for both
$\hat {x}\lessapprox \hat {A}$
and for
$\hat {x}\gtrapprox \hat {A}$
. For
$\hat {x}\lessapprox \hat {A}$
, (2.22) from Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) gives

where
$\eta = (\hat {z} + \hat {A}/2)/\sqrt {H}$
is a non-dimensional distance from the end of the cylinder. Equation (B5) represents the effects of a single end of the cylinder, and assumes that the two ends of the cylinder are sufficiently separated to not significantly interact. Therefore, this solution is not valid for
$\hat {x}\geqslant \hat {A}$
. Noting that the opposite end of the cylinder is situated at
$\eta = \hat {A}/\sqrt {H}$
, and accounting for the end correction of both ends of the cylinder, we obtain the following expression for
$\bar {\hat {g}}$
:

Here the subscript
$(1)$
indicates the solution for
$\hat {x}\lessapprox \hat {A}$
, and we have used
$\sqrt {H}\approx \hat {x}$
. To estimate the asymptotic value of the integral
$f$
, we utilise the following series approximation for
$\bar {\hat {g}}_{(1)}$
, which converges for
$\hat {x}\leqslant \hat {A}$
:

For
$\hat {x}\gtrapprox \hat {A}$
, Teng et al. (Reference Teng, Rallabandi, Stone and Ault2022) provide the following solution for
$\hat {p}$
(their (2.25)):



Noting that, for
$\hat {x} \gg 1$
,
$\hat {g} \approx -\hat {p}/\hat {x}^3$
, and averaging across the cylinder span, we obtain

We then compute the total
$f$
as

For computing
$f_{(1)}$
, we use the asymptotic approximations of the integrals



to obtain an expression for
$f$
, up to order
$\hat {A}^{-1}$
:

Finally, evaluating the sum gives


The first term in
$f_{(2)}$
is:

Making the substitution
$x_b = \hat {x}/\hat {A}$
, the remaining terms of
$f_{(2)}$
are of the form

and therefore, we have
$f_{(2)}$
,

Therefore, the total
$f$
is

A comparison between the
$f$
predicted by (B22) and the
$f$
obtained using our numerical simulations is presented in figure 35(a). Here, we plot
$(1-f)\hat {A}$
against
$\log _{10}{\hat {A}}$
so that the asymptotic expression (B22) becomes linear. The slope of the line corresponds to the coefficient of the
$\ln (\hat {A})/\hat {A}$
term, while the
$y$
intercept corresponds to the
$\hat {A}^{-1}$
term. Equation (B22) correctly predicts the slope of
$f$
(the
$\ln (\hat {A})/\hat {A}$
term), however, it underpredicts the
$y$
intercept (the
$\hat {A}^{-1}$
term).

Figure 35. (a) Comparison between the asymptotic approximation for
$f$
(B22) and the
$f$
computed from our numerical simulations. Blue markers correspond to our numerical data, while red markers correspond to a direct numerical integration of (B12). The dashed line is (B22), while the solid line indicates the best fit approximation to our numerical data. (b) Comparison between the profiles of
$\bar {\hat {g}}$
obtained from our numerical simulations (
$\bar {\hat {g}}_{\mathrm{numerical}}$
), as well as the approximate solutions given in (B7) (
$\bar {\hat {g}}_{(1)}$
) and (B11) (
$\bar {\hat {g}}_{(2)}$
), for
$\hat {A}=100$
.
To account for this error, we note that the approximate solution
$\bar {\hat {g}}_{(1)}$
given by (B7) is only strictly valid for
$1 \ll \hat {x} \ll \hat {A}$
. This is confirmed by figure 35, which plots the numerically obtained profiles of
$\bar {\hat {g}}$
against the approximate solutions
$\bar {\hat {g}}_{(1)}$
and
$\bar {\hat {g}}_{(2)}$
. All three solutions are in excellent agreement for
$1 \ll \hat {x} \ll \hat {A}$
. As expected,
$\bar {\hat {g}}_{(1)}$
is inaccurate for
$\hat {x} \gtrapprox \hat {A}$
, while
$\bar {\hat {g}}_{(2)}$
is accurate for all
$\hat {x} \gg 1$
. Finally, we note that neither
$\bar {\hat {g}}_{(1)}$
nor
$\bar {\hat {g}}_{(2)}$
are accurate for
$\hat {x} \lessapprox 1$
. Specifically, both
$\bar {\hat {g}}_{(1)}$
and
$\bar {\hat {g}}_{(2)}$
approach
$1$
as
$\hat {x}$
approaches
$0$
, while
$\bar {\hat {g}}$
remains less than
$1$
. Evidently, this leads to an
$O(\hat {A}^{-1})$
error in the asymptotic approximation given by (B22).
Given that we do not have an accurate analytic expression for
$\bar {\hat {g}}$
valid when
$\hat {x} \lessapprox 1$
, we instead obtain the
$\hat {A}^{-1}$
term by fitting to our numerically obtained
$f$
, to give the following expression:

As shown in figure 35(a), this expression agrees with our numerical data for large
$\hat {A}$
.