1. Introduction
In this work, the Stokes equation for the flow between two immobile cylinders is solved analytically. Such a configuration has similarities to printing processes, where the two rollers are rotating. The differentiation and interpretation of the solutions with homogeneous and inhomogeneous boundary conditions are central to the present case. The inhomogeneous boundary conditions, i.e. in particular the no-slip condition on the impermeable surface of the rotating rollers, together with the Stokes equation, lead to an analytical solution that describes the flow between the rollers. This case has been discussed at length in the literature and details can be found below. However, it is interesting to note that new and interesting solutions can be constructed for the Stokes equation that are based solely on homogeneous boundary conditions. These describe highly complex flow patterns, which form the core of the present work. The homogeneous solutions can be understood as the response to a disturbance in the far field, i.e. far away from the nip between the two rollers.
Jeffery (Reference Jeffery1922) may have been the first to transform the Stokes equation into a bipolar coordinate system, which allows the two roller surfaces to be represented as constant coordinate lines. He also constructed the first solutions. Another specific solution for the problem at hand has also been published by Müller (Reference Müller1942). There, inhomogeneous solutions for the flow induced by two equally large co- and counter-rotating cylinders are presented. By the superposition principle of linear systems, these can be combined to obtain a solution for arbitrary rotational speeds of the cylinders. For this reason, the discussion of these particular solutions is skipped within this work and we restrict our focus to the solutions of the homogeneous problem.
Pitts & Greiller (Reference Pitts and Greiller1961) analysed the flow between two cylinders partially immersed in liquid experimentally and then theoretically, neglecting the inertia term. They also found the position of the free stagnation points for fully immersed cylinders. A further specific solution for the Stokes flow in bipolar coordinates is also given by Wakiya (Reference Wakiya1975b ) in the second part of his three-part work. In it, he first provides a general solution for the behaviour of the Stokes flow in the vicinity of two cylinders. He starts from the solution of the biharmonic equation in the bipolar coordinate system given in Jeffery (Reference Jeffery1922). He also assumes that the cylinders cannot rotate independently. A more general approach to the problem is taken by Dorrepaal & O’Neill (Reference Dorrepaal and O’Neill1979), who solved the biharmonic equation of the streamfunction using ‘matched asymptotic expansions’. The following papers also consider the Stokes flow in the bipolar coordinate system, but with different boundary conditions. In the first part of Wakiya (Reference Wakiya1975a ) the ‘flow along a plane with a projection or a depression’ is considered and in the third part (Wakiya Reference Wakiya1978) the Stokes flow of steadily moving, rotating cylinders and of a cylinder eccentrically rotating within another cylinder is investigated. The latter is also investigated by Kazakova (Reference Kazakova2020). Finally, Meleshko & Gomilko (Reference Meleshko and Gomilko2000) analysed the Stokes flow within a semicircle.
The solution of the problem with homogeneous boundary conditions leads to an eigenvalue problem for the sought-after streamfunction, which can only be solved numerically. The resulting eigenvalues generate various associated flow patterns, which are characterised by vortices rotating alternately in opposite directions. Several results showing such viscous vortices are collected in the monograph by Shankar (Reference Shankar2007), including also the works on inhomogeneous boundary conditions of Jeffery (Reference Jeffery1922) and Dorrepaal & O’Neill (Reference Dorrepaal and O’Neill1979). A very well-known phenomenon of viscous vortex formation, and related to the present flow, is that in a corner with a moving wall, which was described by Moffatt (Reference Moffatt1964). Moffatt shows that a cascade of increasingly smaller counter-rotating vortices forms between two rigid walls at an opening angle of less than
$146^\circ$
within the corner. We will see later that the leading mode for the present flows is comparable to the Moffatt modes. Similar to Moffatt’s approach the emergence of viscous vortices in the cusp between two touching cylinders was investigated in Schubert (Reference Schubert1967). The recent work of Dormy & Moffatt (Reference Dormy and Moffatt2024) utilised the same approach to solve the flow in this cusp when the cylinders are rotating. These works differ from the present work in the utilised coordinate system, where we are not able to investigate the limiting case of touching cylinders directly and the coordinate transform of Schubert (Reference Schubert1967) and Dormy & Moffatt (Reference Dormy and Moffatt2024) is restricted to this case.
As mentioned above, in this work we employ the stationary Stokes equation via the biharmonic equation of the streamfunction in bipolar coordinates, which map the surfaces of the rollers to constant coordinate lines. By means of a further transformation, see e.g. Bluman & Gregory (Reference Bluman and Gregory1985), the streamfunction is converted into a linear differential equation with constant coefficients in the bipolar coordinate system. The general solution of this differential equation is derived in the Appendix A. When solving the two-roller printing problem, we distinguish, as indicated above, between an inhomogeneous and a homogeneous boundary value problem. For the inhomogeneous boundary conditions, the speeds at the two roller surfaces are specified independently of each other, as results for this case were already derived by Jeffery (Reference Jeffery1922) and Müller (Reference Müller1942), it is not repeated here. For the problem with homogeneous boundary conditions the rollers are virtually at rest, that is, immobile and impermeable. In § 3, the homogeneous problem is converted to an algebraic eigenvalue problem using the general solution from the Appendix A. The algebraic eigenvalue problem is solved numerically. The flow problem with homogeneous boundary conditions is then analysed and illustrated in detail in §§ 3.2 and 3.3, e.g. by plotting streamlines and by investigating three different types of stagnation points, so that the behaviour for varying eigenvalues and variable gap height becomes apparent. The analysis of the vortex spacing and decay of their relative intensity is conducted in § 3.4 Finally, the results are discussed in § 4 and classified in relation to the aforementioned works.
2. Formulation of the Stokes equation using bipolar coordinates
For the problem at hand, the use of a bipolar coordinate system (BCS), as shown in figure 1, is appropriate, since here the surfaces of the rollers are described by constant coordinate lines of the
$\eta$
coordinate. This simplifies the prescription of arbitrary boundary conditions (in our case homogeneous boundary conditions) on the roller surfaces and allows an exact solution to be found. Two different solution branches with respect to the
$\xi$
coordinate can be identified. These display a flow that either decays from left to right, i.e. in the positive
$\xi$
-direction, denoted with subindex
$\boldsymbol{l}$
, or from right to left, denoted with subindex
$\boldsymbol{r}$
. Through superposition of solutions, other combinations can be constructed, for example streamfunctions that are even (
$\boldsymbol{e}$
) or odd (
$\boldsymbol{o}$
) in
$\xi$
with respect to the nip at
$\xi =\pi$
, this form is also included in Appendix A table 3 but not further discussed in this work. The latter procedure can also be used for the
$\eta$
coordinate. These solution branches are labelled as symmetric (
$\boldsymbol{s}$
), with an odd streamfunction with respect to
$\eta$
, and asymmetric (
$\boldsymbol{a}$
). This results in a total of four solution types, which in principle can be superposed due to the linearity of the problem. In this section, the bipolar coordinates are first defined and then the Stokes equation is transformed into the BCS. Finally, the boundary conditions for the homogeneous solution are defined.
Figure 1. Coordinate lines of the BCS for an arbitrary
$c$
. The blue and red lines represent constant
$\eta$
and
$\xi$
lines respectively. The red crosses at
$(x,y) = (0,\pm c)$
mark the focal points, where all constant
$\xi$
lines meet. In black, constant
$\eta$
curves are shown that correspond to two equal sized rollers of radius
$R$
, separated by a gap of
$2h$
and whose midpoints are marked by black crosses.
2.1. Bipolar coordinates
In this work, the BCS is rotated
$90^\circ$
anti-clockwise compared with the usual definition, see Happel & Brenner (Reference Happel and Brenner1981, p. 498) or Moon & Spencer (Reference Moon and Spencer1971, p. 89), in order to comply with the usual convention of the flow direction from left to right. This results in the following definition in complex representation:
\begin{equation} x + \mathrm{i} y = -c\ \mathrm{cot}\left (\frac {1}{2} (\xi + \mathrm{i} \eta ) \right ), \end{equation}
where c is defined below. Here,
$z = x + \mathrm{i} y$
represent the Cartesian coordinates and
$\zeta = \xi + \mathrm{i}\eta$
denote the bipolar coordinates. The limits for the coordinates are
$-\infty \lt \eta \lt \infty$
and
$0 \leqslant \xi \lt 2\pi$
. At
$(x,y) = (0,c)$
and
$(x,y) = (0,-c)$
all constant
$\xi$
lines meet, at the same time these two coordinate points represent the positions
$\eta = \infty$
in the upper half and
$\eta = -\infty$
in the lower half. Equation (2.1) can be used to find explicit expressions for the Cartesian coordinates as a function of the bipolar coordinates
\begin{equation} x = \frac {- c\ \sin \xi }{\cosh \eta - \cos \xi }\ , \quad y = \frac {c\ \sinh \eta }{\cosh \eta - \cos \xi }. \end{equation}
The Jacobian determinant
$J$
(see Bluman & Gregory Reference Bluman and Gregory1985) plays an important role in the further course of this paper and is therefore already presented here
\begin{equation} J = \biggl | \frac {\rm d z}{\rm d \zeta } \biggr | = \frac {c}{\cosh \eta - \cos \xi }. \end{equation}
Half the gap width
$h$
and the radius of the rollers
$R$
serve as the specified length scales, see figure 1. The centre of the rollers is
$y_{{RC}} = h+R$
and the surfaces of the two rollers are each constant
$\eta$
coordinate lines, where
$\eta _0$
stands for the upper roller and
$\eta _1$
for the lower roller. This also defines the boundary conditions on constant
$\eta$
coordinate lines, which allows a comparatively simple construction of the solution. According to Happel & Brenner (Reference Happel and Brenner1981, p. 498), the centre point and the roller radius of a constant
$\eta$
coordinate line are
\begin{equation} y_{{RC}} = c \coth (\eta ) = h + R, \end{equation}
and
\begin{equation} R = \frac {c}{|\sinh (\eta )|}. \end{equation}
From (2.4) and (2.5) the constant
$\pm \eta _0$
for the surfaces of the rollers with same diameter can be determined as a function of the dimensionless gap ratio
$h/R$
, as well as the dimensionless focal point
$c/R$
\begin{equation} \eta _0 = \cosh ^{-1}\left (\frac {h}{R}+1\right ). \end{equation}
2.2. Transformation of the stokes streamfunction equation to bipolar coordinates
The stationary Stokes equation in dimensionless variables reads
\begin{equation} \boldsymbol{\nabla }p = \Delta \boldsymbol{u}, \end{equation}
and can be represented by the biharmonic equation of the streamfunction (see Happel & Brenner Reference Happel and Brenner1981, p. 60), whose curl determines the velocities. The connection between the velocity components and streamfunction in the two-dimensional BCS is
\begin{equation} u_\xi = \frac {1}{J}\ \frac {\partial \psi }{\partial \eta }\ , \quad u_\eta = - \frac {1}{J}\ \frac {\partial \psi }{\partial \xi }. \end{equation}
The biharmonic equation
\begin{equation} \varDelta ^2 \psi = 0. \end{equation}
is obtained by taking the curl of (2.7) to eliminate the pressure. In order to represent the biharmonic equation of the streamfunction, see (2.9), in the BCS as a linear differential equation with constant coefficients, it is first represented as the product of the Jacobian determinant, see (2.3), and a modified streamfunction
$\varPhi (\xi ,\eta )$
(see Bluman & Gregory (Reference Bluman and Gregory1985)
\begin{equation} \psi (\xi ,\eta ) = J \varPhi (\xi ,\eta ). \end{equation}
The biharmonic equation in terms of the modified streamfunction
$\varPhi (\xi ,\eta )$
is then obtained by substituting (2.10) into (2.9) and expanding the biharmonic operator in the BCS
\begin{equation} \varDelta ^2 \psi (\xi ,\eta ) = \frac {1}{J^3}\left [\left (\frac {\partial ^2}{\partial \xi ^2} + \frac {\partial ^2}{\partial \eta ^2}\right )^2 + 2 \left (\frac {\partial ^2}{\partial \xi ^2} - \frac {\partial ^2}{\partial \eta ^2}\right ) + 1\right ] \varPhi (\xi ,\eta ) = 0. \end{equation}
The general product ansatz solutions for the above (2.11) are listed in Appendix A.
2.3. Formulation of boundary conditions
In order to be able to determine the constants present in the general solution (see Appendix A), the following homogeneous boundary conditions are used as a basis. The kinematic boundary condition of the impermeable wall applies
\begin{equation} u_{\eta }(\xi ,\eta = \eta _0) = u_{\eta }(\xi ,\eta = -\eta _0) = 0, \end{equation}
as well as a no-slip condition with zero velocities
\begin{equation} u_{\xi }(\xi ,\eta = \eta _0) = u_{\xi }(\xi ,\eta = -\eta _0) = 0. \end{equation}
To determine the behaviour of the solution in
$\xi$
and since the BCS is
$2\pi$
-periodic, periodic boundary conditions could be used. However, as the focus is currently on the nip, we will not pursue this further here.
3. Solution of the homogeneous Stokes problem, eigenvalues and vortex solutions
Solutions to the Stokes problem with non-vanishing velocities of the rollers, i.e. a non-vanishing right-hand side in (2.13) can be superimposed with solutions of the homogeneous problem, which are presented within this section. The so-constructed patterns can be seen as the response of the flow to a disturbance in the far field, i.e. far away from the nip between the rollers. For this purpose, the solution path
$\textit{II}$
from table 3 of Appendix A is selected: it can be shown that the general solutions of the solution paths
$\textit{I}$
and
$\textit{II}$
coincide. Furthermore, the general solutions of the solution paths
$\textit{III}$
and
$\textit{IV}$
only provide the trivial solution based on the given homogeneous boundary (2.13). Accordingly, the following streamfunction is used with
$f_{\boldsymbol{l}/\boldsymbol{r}}$
and
$g_{\boldsymbol{a}/\boldsymbol{s}}$
as defined in table 3 in Appendix A:
\begin{equation} \psi (\xi ,\eta ;\alpha ) = J(\xi ,\eta ) \bigl [f_{\boldsymbol{l}}(\xi ;\alpha ) + f_{\boldsymbol{r}}(\xi ;\alpha )\bigr ] \bigl [ g_{\boldsymbol{a}}(\eta ;\alpha ) + g_{\boldsymbol{s}}(\eta ;\alpha ) \bigr ]. \end{equation}
For the eigenvalues still to be determined
$\alpha \in \mathbb{C} \setminus \{-\mathrm{i},-1,0,1,\mathrm{i}\}$
applies. The excluded values are required to construct the particulate solutions, which correspond to a non-zero rotation of the rollers.
3.1. Determination of the eigenvalues
Inserting the homogeneous boundary conditions from (2.13) into the streamfunction from (3.1) and rearranging leads to the following four equations:
\begin{align} 0 &= A (\alpha +1) \sinh ((\alpha +1)\eta _0) + C (\alpha -1) \sinh ((\alpha -1)\eta _0), \end{align}
\begin{align} 0 &= B (\alpha +1) \cosh ((\alpha +1)\eta _0) + D (\alpha -1) \cosh ((\alpha -1)\eta _0), \end{align}
\begin{align} 0 &= A \cosh ((\alpha +1)\eta _0) + C \cosh ((\alpha -1)\eta _0), \end{align}
\begin{align} 0 &= B \sinh ((\alpha +1)\eta _0) + D \sinh ((\alpha -1)\eta _0). \end{align}
It is immediately visible that this system of four equations can be decoupled into two two-equation systems, leading to two independent solutions. Equations (3.2a
) and (3.2c
) result in a relation between integration constants
$A$
and
$C$
and an equation to determine the discrete eigenvalues as
\begin{align} C &= - \frac {\cosh ((\alpha + 1)\eta _0)}{\cosh ((\alpha - 1)\eta _0)}\ A, \end{align}
\begin{align} 0&=\sinh (2\alpha \eta _0) + \alpha \sinh (2\eta _0) . \end{align}
The streamfunction of this solution branch (contribution
$g_{\boldsymbol{a}}$
in (3.1)) is even in the
$\eta$
coordinate and leads to an asymmetrical flow. Therefore, it is labelled with the index
$\boldsymbol{a}$
. On the other hand, the second solution branch, specified by (3.2b
), (3.2d
), result in
\begin{align} D &= - \frac {\sinh ((\alpha + 1)\eta _0)}{\sinh ((\alpha - 1)\eta _0)}\ B, \end{align}
\begin{align} 0&=\sinh (2\alpha \eta _0) - \alpha \sinh (2\eta _0) . \end{align}
For this solution the streamfunction (contribution
$g_{\boldsymbol{s}}$
in (3.1)) is odd in
$\eta$
and leads to symmetrical flow patterns. Therefore, it is labelled with the index
$\boldsymbol{s}$
. These two solutions lead to different eigenvalues and either of the branches can be superposed with any solution (in particular the trivial solution, i.e.
$B=D=0$
or
$A=C=0$
) for the other branch.
The eigenvalues
$\alpha =\alpha _r+i\alpha _i$
are determined numerically by first splitting (3.3b
) (or (3.4b
)) into its real and imaginary parts (for the determination of the
$\boldsymbol{a}$
eigenvalues the sign of the first term is positive, for the
$\boldsymbol{s}$
eigenvalues it is negative)
\begin{align} &\pm \sinh (2\alpha _r\eta _0)\cos (2\alpha _i\eta _0)+\alpha _r\sinh (2\eta _0) = 0, \end{align}
\begin{align} &\pm \cosh (2\alpha _r\eta _0)\sin (2\alpha _i\eta _0)+\alpha _i\sinh (2\eta _0) = 0, \end{align}
which can be reshaped to
\begin{align} \cos (2\alpha _i\eta _0) &= \mp \alpha _r\frac {\sinh (2\eta _0)}{\sinh (2\alpha _r\eta _0)} , \end{align}
\begin{align} \sin (2\alpha _i\eta _0) &= \mp \alpha _i\frac {\sinh (2\eta _0)}{\cosh (2\alpha _r\eta _0)} . \end{align}
Solving (3.6a
)
$^2+$
(3.6b
)
$^2=1$
reveals the same relation between the imaginary and real parts of the eigenvalue, valid for both solution branches and all modes
$n \in \mathbb{N}\setminus \{0\}$
, namely
\begin{equation} \alpha _{i,n} = \pm \frac {\cosh (2\alpha _{r,n}\eta _0)}{\sinh (2\eta _0)}\sqrt {1-\left (\alpha _{r,n}\frac {\sinh (2\eta _0)}{\sinh (2\alpha _{r,n}\eta _0)}\right )^2}. \end{equation}
However, for the numerical determination of the
$n$
th eigenvalue it proved to be more suitable to compute the relation between its imaginary and real part solely from (3.5a
), resulting in
\begin{align} \boldsymbol{a} : \alpha _{i,n} &= \frac {1}{2\eta _0}\left (\pm \cos ^{-1}\left (\frac {\alpha _{r,n}\sinh (2\eta _0)}{\sinh (2\eta _0\alpha _{r,n})}\right ) + 2(n-1)\pi + \pi \right ), \end{align}
\begin{align} \boldsymbol{s} : \alpha _{i,n} &= \frac {1}{2\eta _0}\left (\pm \cos ^{-1}\left (\frac {\alpha _{r,n}\sinh (2\eta _0)}{\sinh (2\eta _0\alpha _{r,n})}\right ) + 2n\pi \right ). \end{align}
This relation is then inserted into (3.5b
), which has exactly one root if the sign in front of
$\cos ^{-1}$
is chosen as
$+$
in (3.8a
) or (3.8b
) respectively, and solved numerically for
$\alpha _r$
. For any eigenvalue
$\alpha _n$
the conjugate and negative values
$-\alpha _n,\overline {\alpha }_n,-\overline {\alpha }_n$
are also valid eigenvalues. It can be shown that, in reference to the eigenvalue with
${\rm Re}(\alpha _n)\gt 0, \textrm{Im}(\alpha _n)\gt 0$
, these combinations lead to phase shifts, scalings and transformations from one solution branch to another. It is therefore sufficient to investigate only those eigenvalues positioned in the first quadrant. In table 1 the first three eigenvalues for both the asymmetric and symmetric solution branches are displayed for three different ratios of gap width to roller radius
$h/R$
, whose relation to
$\eta _0$
is given by (2.6).
Table 1. Numerically calculated eigenvalues.
3.2. Solutions for the four flow patterns
After one specific eigenvalue is determined and thereby the behaviour of the corresponding eigenmode in
$\eta$
as well, we still need to specified how to choose the constants
$\hat {E}$
and
$\hat {F}$
occurring in (3.1). Notice how
$f_{\boldsymbol{l}}$
decays exponentially for
$\textrm{Im}(\alpha )\gt 0$
, while
$\xi$
is passing from
$0$
to
$2\pi$
. Therefore,
$f_{\boldsymbol{l}}$
and
$f_{\boldsymbol{r}}$
control the behaviour of the eigenmode to either decay from left to right (index
$\boldsymbol{l}$
), in which case
$\hat {F} = 0$
, or to decay from right to left (index
$\boldsymbol{r}$
), in which case
$\hat {E} = 0$
. Their contribution to the solution is based on the location of the disturbance inducing the flow. Combining all multiplicative constants in a single complex value
$s$
the final form of the streamfunction is obtained as
\begin{equation} \psi (\xi ,\eta ;\alpha )_{\boldsymbol{l}/\boldsymbol{r},\boldsymbol{a}/\boldsymbol{s}} = s J(\xi ,\eta ) F_{\boldsymbol{l}/\boldsymbol{r}}(\xi ;\alpha ) G_{\boldsymbol{a}/\boldsymbol{s}}(\eta ;\alpha ), \end{equation}
which allows for four distinct solution branches, depending on which behaviour in
$\xi$
or
$\eta$
is displayed
\begin{align} F_{\boldsymbol{l}}(\xi ;\alpha ) &= e^{\mathrm{i}\alpha (\xi -\pi )}, \end{align}
\begin{align} F_{\boldsymbol{r}}(\xi ;\alpha ) &= e^{-\mathrm{i}\alpha (\xi -\pi )}, \end{align}
\begin{align} G_{\boldsymbol{a}}(\eta ;\alpha ) &= \frac {\cosh ((\alpha +1)\eta )}{\cosh ((\alpha +1)\eta _0)} - \frac {\cosh ((\alpha -1)\eta )}{\cosh ((\alpha -1)\eta _0)}, \end{align}
\begin{align} G_{\boldsymbol{s}}(\eta ;\alpha ) &= \frac {\sinh ((\alpha +1)\eta )}{\sinh ((\alpha +1)\eta _0)} - \frac {\sinh ((\alpha -1)\eta )}{\sinh ((\alpha -1)\eta _0)}. \end{align}
In the following sections further investigations of the solutions are carried out for
${\rm Re}(\psi _{\boldsymbol{l},\boldsymbol{a}})$
and
${\rm Re}(\psi _{\boldsymbol{l},\boldsymbol{s}})$
and restricted to eigenvalues from the first quadrant. The obtained conclusions can be transferred to the
$\boldsymbol{r}$
solution branches by reflection
$\xi =2\pi -\xi$
.
3.3. Characteristic flow patterns
Figure 2 shows the constant streamlines of
$\psi _{\boldsymbol{l},\boldsymbol{a}}$
for the first three asymmetric eigenvalues at a ratio
$h/R=0.01$
. The zero iso-contours are shown in black, and blue and red curves denote negative and positive iso-contours, respectively. For the first eigenvalue in figure 2(a), vortices with alternating directions of rotation form, which occupy the full channel width. The magnitude of the streamfunction decreases from left to right; the decay factor is investigated in § 3.4. Stagnation points only form on the walls.
Figure 2. Plots of constant streamlines of the streamfunction
$\psi _{\boldsymbol{l},\boldsymbol{a}}$
around
$\xi = \pi$
for increasing
$\alpha _n$
and
$h/R=0.01$
. (a)
$\alpha _1 = 8.0152 + 14.896i$
, (b)
$\alpha _2 = 11.028 + 37.902i$
, (c)
$\alpha _3 = 12.613 + 60.411i$
, (d) Decay of
$|\psi _{\boldsymbol{l},\boldsymbol{a}}|$
along
$\eta =0$
.
For the second eigenvalue in figure 2(b), two vortices of the same rotation direction form, which are stretched finger-like towards the right side of the channel and protrude into the neighbouring vortex structure. In addition, the vortex density, i.e. the number of vortices per unit length, increases compared with the first eigenvalue.
The third eigenvalue in figure 2(c) shows that the tendency towards a ‘finger-like’ flow continues. As with the second eigenvalue, strongly elongated vortices are formed in the same direction at the rollers. Towards the right side of the gap between the rollers, these vortices are deformed in a ‘stepped’ manner and push under the next vortex. Again, this results in ‘finger-like’ flow structures. The number of these steps remains constant along the horizontal axis. The vortex density continues to increase with increasing eigenvalue.
To provide a better understanding of how fast the intensity of these vortices decreases, the absolute value of the streamfunction
$\psi _{\boldsymbol{l},\boldsymbol{a}}$
along
$\eta =0$
is displayed in figure 2(d) for all three modes. The segments of alternating signs are displayed as red (positive) and blue (negative) lines (due to the logarithmic representation and finite resolution, the zero points are not explicitly shown). The slowest decay is observable for the first mode, but even here after the third zero crossing the relative intensity decreased by more than five orders of magnitude. Therefore, experimental observation appears to be possible only for a small number of these vortices. Furthermore, it is visible that the intensity of the higher modes is decreasing extremely fast and detection of these modes in experiments would be highly unlikely.
The same analysis is repeated in figure 3 for the constant streamlines of
$\psi _{\boldsymbol{l},\boldsymbol{s}}$
for the first three symmetric eigenvalues at a ratio
$h/R=0.01$
. First, it can be seen that there is a horizontal separation line in the flow. This causes several free stagnation points to occur in the flow along the horizontal dividing line. Above and below the horizontal dividing line, vortices form with alternating direction of rotation. The streamfunctions decrease in magnitude towards the right side of the gap. As observed for the asymmetrical solution branch the vortices are oriented towards the right.
Figure 3. Plots of constant streamlines of the streamfunction
$\psi _{\boldsymbol{l},\boldsymbol{s}}$
around
$\xi = \pi$
for increasing
$\alpha _n$
and
$h/R=0.01$
. (a)
$\alpha _1 = 9.8456 + 26.525i$
, (b)
$\alpha _2 = 11.909 + 49.181i$
, (c)
$\alpha _3 = 13.199 + 71.611i$
, (d) Decay of
$|\psi _{\boldsymbol{l},\boldsymbol{s}}|$
along
$\eta = \eta_0/2$
.
Figure 4. Depiction of three different kinds of stagnation points.
For higher eigenvalues, a finger-like flow results, similar to the asymmetric flow, whereby, in contrast to the asymmetric flow, no central vortex forms on the horizontal axis of symmetry. The vortex density increases with increasing eigenvalues. A schematic representation of the three different types of stagnation points arising in the flow is depicted in figure 4. Stagnation points occur at extrema of the streamfunction, which can be true minima/maxima (left), saddle points (middle) or saddle points at a zero iso-contour (right). In figures 5(a) and 5(b) a zoom into two distinct areas of figures 2(b) and 3(b) is visible. On the one hand, the transition of the wall-bound vortices into the elongated finger structure can be observed, and on the other hand, the three different types of stagnation points from figure 4 are clearly recognisable. In figure 5(a), the visible streamlines imply a stagnation point at the saddle points of co-rotating streamlines (approximately
$x\approx 0.006, y\approx \pm 0.0018$
), represented by round markers. In addition, the formation of the elongated, ‘finger-like’ vortex structures can be recognised, which encloses a single vortex, whose centre of rotation is marked by a cross. In addition to these stagnation points, the symmetric solution branch, shown in figure 5(a), exhibits an additional saddle point with a change of sign. At this location a free stagnation point forms in the flow, marked by a square marker. Finally, two crosses represent the approximate locations of the centres of rotation of the wall-bound vortices.
3.4. Spacing and decay of the vortex chain
To conclude the analysis of the identified flow patterns, this section investigates the distance between two consecutive vortices, as well as how fast their relative intensity decays. First, the
$\eta$
-velocity component, according to (2.8), and its derivative with respect to
$\xi$
are constructed
\begin{align} u_{\eta ;\boldsymbol{l},\boldsymbol{a}/\boldsymbol{s}}(\xi ,\eta ;\alpha _n) &= -s e^{i\alpha _n(\xi -\pi )}G_{\boldsymbol{a}/\boldsymbol{s}}(\eta ;\alpha _n)\left (i\alpha _n - \frac {\sin (\xi )}{\cosh (\eta )-\cos (\xi )}\right ), \end{align}
\begin{align} \frac {\partial u_{\eta ;\boldsymbol{l},\boldsymbol{a}/\boldsymbol{s}}}{\partial \xi }(\xi ,\eta ;\alpha _n) &= -s e^{i\alpha _n(\xi -\pi )}G_{\boldsymbol{a}/\boldsymbol{s}}(\eta ;\alpha _n) \\ \nonumber &\quad \times \left (-\alpha _n^2 - \frac {i\alpha _n\sin (\xi )}{\cosh (\eta )-\cos (\xi )} +\frac {1-\cosh (\eta )\cos (\xi )}{(\cosh (\eta )-\cos (\xi ))^2} \right ). \end{align}
The integration constant
$s$
and
$G_{\boldsymbol{a}/\boldsymbol{s}}(\eta ;\alpha _n)$
are merely constant complex prefactors when evaluating (3.11a
) and (3.11b
) along a constant
$\eta$
line, e.g.
$\eta =0$
. Hence, they merely cause a phase shift
$\epsilon$
and scaling by a real factor. For the symmetric solution branch the velocity vanishes along
$\eta =0$
and the analysis could be carried out e.g. at
$\eta =\eta _0/2$
, however, the overall spacing and decay behaviour is unaffected for higher modes and the symmetric branch. Additionally, in the gap region
$\pi /2\lt \xi \lt 3\pi /2$
and for sufficiently large eigenvalues
$|\alpha _{r,n}|, |\alpha _{r,n}| \gt 1$
the parenthesis in (3.11a
) and (3.11b
) is almost constant and the equations are dominated by the complex exponential function. We therefore obtain the following proportionalities:
\begin{align} u_{\eta ;\boldsymbol{l},\boldsymbol{a}/\boldsymbol{s}}(\xi ,0;\alpha _n) & \propto e^{i(\alpha _n(\xi -\pi )+\epsilon )}, \end{align}
\begin{align} \frac {\partial u_{\eta ;\boldsymbol{l},\boldsymbol{a}/\boldsymbol{s}}}{\partial \xi }(\xi ,0;\alpha _n) & \propto e^{i(\alpha _n(\xi -\pi )+\epsilon )}. \end{align}
Now, by finding the roots of (3.12a
) in
$\xi$
, the approximate locations
$\xi _0$
of vortex centres are obtained and the locations
$\xi _{+}$
of velocity maxima by repeating the same for (3.12b)
. The distance
$d$
between consecutive vortex centres or velocity maxima is then computed by forming the difference of subsequent locations. As the behaviour of both quantities is dominated by the complex exponential function, this results in the same value of
\begin{equation} d(\alpha _n) = |\xi _{0/+,m+1}-\xi _{0/+,m}|\approx \frac {\pi }{\alpha _{r,n}}. \end{equation}
Additionally, the decay rate
$k$
can be computed from the relative magnitude of the velocity at two consecutive velocity maxima. Using the approximate relation (3.12a
), this results in
\begin{equation} k(\alpha _n) = \left |\frac {u_{\eta }(\xi _{+,m+1},0;\alpha _n)}{u_{\eta }(\xi _{+,m},0;\alpha _n)}\right | \approx \left | \frac {e^{i(\alpha _n(\xi _{+,m+1}-\pi )+\epsilon )}}{e^{i(\alpha _n(\xi _{+,m}-\pi )+\epsilon )}} \right | =e^{-\pi \frac {\alpha _{i,n}}{\alpha _{r,n}}}, \end{equation}
for the approximate distance
$d$
between two consecutive vortices in
$\xi$
and the factor
$k$
, denoting the decay in their relative intensity.
The quantities displayed in (3.13) and (3.14), associated with the
$n$
th eigenmode, are functions of the corresponding eigenvalue, which ultimately depends on the gap ratio
$h/R$
. Following the procedure given in § 3.1, these eigenvalues can be determined in the limit
$h/R\rightarrow 0$
. In figure 6 this analysis is carried out for the decay factor
$k$
given by (3.14) for both the asymmetric and symmetric solution branches and the first three modes
$n=1,2,3$
. The most slowly decaying mode for all gap ratios is the first mode of the asymmetric solution branch. Furthermore, the decay factors of both branches and for all modes approach an asymptotic limit as
$h/R\rightarrow 0$
. As the first eigenmodes display the slowest decay, they dominate the behaviour of the solution. The two limiting values for the first mode in each branch are
\begin{equation} \lim _{h/R\rightarrow 0} \frac {1}{k(\alpha _{\boldsymbol{a},1})} \approx 357.7, \quad \frac {1}{k(\alpha _{\boldsymbol{s},1})} \approx 4952. \end{equation}
Figure 6. Decay factor
$k$
from (3.14) displayed as a function of
$h/R$
for the asymmetric and symmetric solution branches and the first three eigenvalues.
The
$n+1$
th mode decays approximately two orders of magnitude faster than the
$n$
th mode, in both solution branches.
Figure 7. Number of vortices present in the domain.
Figure 8. Iso-contours of
$\psi _{\boldsymbol{l,a}}$
at the critical radius and slightly below for phase
$\arg (s)=2.3672$
. (a)
$h/R$
= 0.5342, (b)
$h/R$
= 0.5289.
Similarly, one can examine the number of sign changes in (3.11a
), where each sign change corresponds to a vortex. Therefore, the total number of vortices
$N_{v\textit{ortex}}$
in the domain equals the number of sign changes. In figure 7
${\rm Re} (u_{\eta ;\boldsymbol{l},\boldsymbol{a}/\boldsymbol{s}}(\xi ,0;\alpha _{\boldsymbol{a},n}) )$
and
${\rm Re} (u_{\eta ;\boldsymbol{l},\boldsymbol{a}/\boldsymbol{s}}(\xi ,\eta _0/2;\alpha _{\boldsymbol{s},n}) )$
are evaluated numerically and the minimum number of sign changes, occurring over all phases
$\arg (s)\in [0,\pi ]$
, as
$\xi$
passes from
$0$
to
$2\pi$
, is counted. First, it becomes evident from the slopes in the limit
$h/R\rightarrow 0$
in figure 7 that the number of vortices increases as
\begin{equation} N_{v\textit{ortex}} \propto \left (\frac {h}{R}\right )^{-\frac {1}{2}}\quad \mathrm{or} \quad |\alpha _r| \propto \left (\frac {h}{R}\right )^{-\frac {1}{2}}. \end{equation}
Second, it is immediately visible that the first asymmetric mode possesses the fewest vortices inside the domain. Additionally, it is determined that there exists a critical ratio
\begin{equation} \left (\frac {h}{R}\right )_{\textit{crit}}\approx 0.5342, \quad \arg (s)_{\textit{crit}}\approx 2.3672 ,\end{equation}
below which more than one vortex arises for all modes and all solution branches. Up to the critical ratio the most dominant solution branch may exhibit one single global vortex, as displayed in figure 8(a). Below the critical ratio,
$h/R=0.99(h/R)_{\textit{crit}}$
in figure 8(b), an additional vortex arises on the right side of the gap. For decreasing ratios additional vortices will form.
4. Conclusion
In this paper, the stationary Stokes flow between two immobile cylinders is investigated using a streamfunction approach. Following Bluman & Gregory (Reference Bluman and Gregory1985), the biharmonic equation for the streamfunction is transformed into a linear differential equation with constant coefficients in the BCS. Then, the general solution for the streamfunction in bipolar coordinates is derived, see Appendix A.
The focus of the present work then lies in finding solutions to the Stokes problem using homogeneous boundary conditions, in which case the rollers are defined as stationary impermeable walls. This leads to an eigenvalue problem, which is constructed analytically and then solved numerically. The eigenvalue problem results in a division of the flow pattern into asymmetric and symmetric flow, describing the symmetry of the velocity fields around
$\eta = 0$
. Additionally, the streamfunction can be constructed so that it exhibits a decaying behaviour of the flow from left to right or right to left. Combining the behaviour in
$\xi$
and
$\eta$
a total of four different solution paths are identified, as summarised by (3.9).
All four solution branches demonstrate the formation of viscous vortices in the nip, similar to the corner vortices first observed by Moffatt (Reference Moffatt1964). Moffatt describes that, for smaller opening angles between two rigid plates, the number of vortices in the corner increases. This behaviour is similar to what is observed in this work for a reduction of the gap ratio, where the vortex density, i.e. the number of vortices per unit length, increases with decreasing gap ratio
$h/R$
and for higher eigenmodes
$\alpha _n, n \gt 1$
. Furthermore, for higher modes, finger-like structures appear in the flow, while the magnitude of the streamfunction decays faster.
Finally, an analysis of the vortex spacing and relative intensity revealed that the periodicity is determined by the real part of the eigenvalue and the relative intensity between two consecutive vortices is given by the quotient of the imaginary to real part of the eigenvalues. For the slowest decaying eigenmode, i.e. the most dominant one, a relative intensity of subsequent vortices of
$k^{-1}\approx 357.7$
is identified in the limit
$h/R\rightarrow 0$
, while the number of vortices tends towards infinity. This value is very close to the minimum decay factor of
$k^{-1}\approx 365$
identified for the case of viscous vortices between two rigid plates in Moffatt (Reference Moffatt1964).
With the equations presented here for the streamfunction, (3.9), it is not possible to represent the extreme case
$h/R=0$
, as then
$\eta _0 = 0$
and
$\eta _1 = 0$
apply and thus the roller surfaces coincide in the BCS. This interesting case requires a different solution approach, e.g. the matched asymptotic expansion chosen by Dorrepaal & O’Neill (Reference Dorrepaal and O’Neill1979) or an alternative coordinate transform as that employed in Schubert (Reference Schubert1967) and Dormy & Moffatt (Reference Dormy and Moffatt2024). Although Dorrepaal & O’Neill (Reference Dorrepaal and O’Neill1979) do not distinguish individual eigenvalues and the rollers in their work come into contact, their solutions are qualitatively similar to the symmetric flow for the first eigenvalues from this work, see figure 3(a).
Acknowledgements
We kindly acknowledge non-financial support by the Graduate School CE within the Centre for Computational Engineering at TU Darmstadt.
Funding
The work of M. Rieckmann is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 265191195.
Declaration of interests
The authors report no conflict of interest.
Data availability
At https://doi.org/10.48328/tudatalib-1820 scripts to reproduce all figures displayed in this work are available, as well as the figures themselves and additional illustrations of streamlines and pressure fields not explicitly shown in this print.
Appendix A. General product solution
In this section, the general product solution of the biharmonic equation in the BCS is derived. Based on the (2.10) and (2.11), the following product ansatz is chosen to construct the streamfunction:
\begin{equation} \varPhi (\xi , \eta ) = P(\xi )Q(\eta ), \end{equation}
where
$P$
depends only on
$\xi$
and
$Q$
only on
$\eta$
. Inserting (A1) into (2.11) results in
\begin{equation} P^{(4)}Q + 2\ P^{\prime\prime}Q^{\prime\prime} + \textit{PQ}^{(4)} + 2\ P^{\prime\prime}Q - 2\ \textit{PQ}^{\prime\prime} + \textit{PQ} = 0. \end{equation}
Dividing the above equation with
$\textit{PQ}$
results in
\begin{equation} \frac {P^{(4)}}{P} + 2\ \frac {P^{\prime\prime}Q^{\prime\prime}}{\textit{PQ}} + \frac {Q^{(4)}}{Q} + 2\ \frac {P^{\prime\prime}}{P} - 2\ \frac {Q^{\prime\prime}}{Q} + 1 = 0. \end{equation}
The following procedure follows the solution method shown by Doschoris (Reference Doschoris2012) for differential equations with mixed derivatives: taking derivatives of (A3) with respect to
$\xi$
and
$\eta$
we obtain
\begin{equation} \frac {d^2}{\rm d \xi d \eta }\left (\frac {P^{\prime\prime}Q^{\prime\prime}}{\textit{PQ}}\right ) = \frac {\rm d}{{\rm d} \xi}\left (\frac {P^{\prime\prime}}{P}\right ) \frac {d}{{\rm d} \eta }\left (\frac {Q^{\prime\prime}}{Q}\right )= 0, \end{equation}
from which it follows that one of the subsequent terms must be constant
\begin{equation} \frac {P^{\prime\prime}}{P} = \pm \ \alpha ^2 = \textit{const.}\ , \quad \frac {Q^{\prime\prime}}{Q} = \pm \ \beta ^2 = \textit{const.}\ , \quad \alpha , \beta \in \mathbb{C}. \end{equation}
The fourth derivatives in (A3) are also represented with the introduced constants
$\alpha$
and
$\beta$
by first rearranging the (A5) to
\begin{equation} P^{\prime\prime} = \pm \alpha ^2 P, \quad Q^{\prime\prime} = \pm \beta ^2 Q, \end{equation}
and then differentiating twice
\begin{equation} P^{(4)} = \pm \alpha ^2 P^{\prime\prime}, \quad Q^{(4)} = \pm \beta ^2 Q^{\prime\prime}. \end{equation}
If the equations from (A6) are inserted into the equations from (A7), this results in
\begin{equation} \frac {P^{(4)}}{P} = \alpha ^4 \ , \quad \frac {Q^{(4)}}{Q} = \beta ^4. \end{equation}
Finally, the differential (2.11) results in the differential equation systems given in table 2. Based on these equations, the solutions of the biharmonic equation in the BCS given in the tables 3 and 4 are obtained, which can be superposed due to linearity.
Table 2. Four sets of linear ordinary differential equations (ODE).
Table 3. General solutions to the differential equations in table 2.
Table 4. Special solutions of the differential equations from table 2.
Appendix B. Solution branches and pressure fields
In § 3 four distinct solution branches were found, and the streamlines displayed for two of them, namely the
$\psi _{\boldsymbol{l},\boldsymbol{a}}$
and
$\psi _{\boldsymbol{l},\boldsymbol{s}}$
solutions. To keep a compact structure within this work, streamlines for the other two branches are not shown here. They exhibit the same features as discussed in § 3.3, the main difference being that they decay from right to left. Graphical representation of these branches, for several ratios
$h/R$
are collected and available in the data repository https://doi.org/10.48328/tudatalib-1820.
The pressure gradients can be derived from (2.7), similar to Wakiya (Reference Wakiya1975a ). By inserting the four distinct solutions for the streamfunction in this connection
\begin{equation} \frac {\partial p}{\partial \xi } = \frac {\partial (\Delta \psi )}{\partial \eta }, \quad \frac {\partial p}{\partial \eta } = -\frac {\partial (\Delta \psi )}{\partial \xi }, \end{equation}
and integrating once, a solution for the pressure field can be obtained for each solution branch. For brevity the pressure solutions are not shown here, but several examples can be accessed in the data repository. Within the pressure solutions the same finger-like patterns and alternating signs (with respect to an arbitrary reference pressure
$p_0$
) can be observed. The most notable difference is that the pressure solution corresponding to a symmetric streamfunction (asymmetric flow) is asymmetric with respect to the
$\eta =0$
line.
Additionally, the repository contains an applet to interactively display the streamfunctions for arbitrary ratios
$h/R$
and phases
$\arg (s)$
. For the critical ratio and phase from (3.17), the repository contains two animated videos displaying the variation of the critical ratio at the fixed critical phase and vice versa.
Open access
