1 Introduction
1.1 Scope and motivation
The present work numerically studies the vorticity distribution in laminar vortex rings produced by a piston–pipe apparatus at Reynolds numbers high enough that an edge layer, thinner than the size of the ring, exists at the dividing streamline. The diagnostics to be presented focus on the vorticity distribution in the region interior to the dividing streamline, called the vortex ring bubble. A model shows that the boundary condition provided by the edge layer is needed to account for an extra Reynolds-number dependence, in addition to that accounted for by introducing a scaled time,
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D708}t$
. However, the detailed structure of the edge layer and wake is not considered in the present work. We hope that it is elucidated in the future. For that effort, high-Reynolds-number falling drops (Harper & Moore Reference Harper and Moore1968) should provide a useful analogy: their structure also consists of an advectively balanced state in the interior surrounded by an edge layer that sheds a wake.
One motivation for studying the laminar vorticity distribution is that azimuthal instabilities are sensitive to its precise form (Saffman Reference Saffman1978). Another motivation is that the non-dimensional energy,
$\widetilde{E}\equiv E/(I^{1/2}\unicode[STIX]{x1D6E4}^{3/2})$
, which determines the limiting stroke-to-diameter ratio according to Gharib, Rambod & Shariff (Reference Gharib, Rambod and Shariff1998), depends on the peakiness of the vorticity profile. (Here
$E$
is the kinetic energy,
$I$
is the impulse and
$\unicode[STIX]{x1D6E4}$
is the circulation; density is set to unity.) For example, Hill’s spherical vortex which has uniform
$\unicode[STIX]{x1D701}\equiv \unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}/r$
gives
$\widetilde{E}=0.16$
while the experimentally measured value for the limiting ring is a much larger
$\widetilde{E}\approx 0.33$
. In addition, the dividing streamline of Hill’s vortex has zero oblateness
$(a-b)/b$
while the thickest experimental ring has an oblateness of
$0.21$
as inferred from photographs in Gharib et al. (Reference Gharib, Rambod and Shariff1998),
$a$
and
$b$
being the semi-major and minor axes, respectively. The discrepancies in energy and oblateness are both consequences of the peakiness of
$\unicode[STIX]{x1D701}$
in actual rings.
The family of steady vortex rings (Norbury Reference Norbury1973) with uniform
$\unicode[STIX]{x1D701}$
, in which Hill’s vortex is the thickest member, has been used to model vortex ring formation (Mohseni & Gharib Reference Mohseni and Gharib1998; Shusser & Gharib Reference Shusser and Gharib2000; Linden & Turner Reference Linden and Turner2001). Since thinner rings have larger
$\widetilde{E}$
, the result in these studies is that the limiting ring corresponds to a thinner member of the family than Hill’s vortex. If one were to use a family of rings having a more realistic peaked vorticity distribution, the result would presumably be that the limiting experimental ring corresponds better with the thickest member of that family.
The physics of vortex sheet roll-up at a sharp edge implies a peaked distribution of
$\unicode[STIX]{x1D701}$
. Saffman (Reference Saffman1978) and Pullin (Reference Pullin1979) obtain the form of
$\unicode[STIX]{x1D701}$
for thin rings (and in the inner region of the core) by using Kirde’s (Reference Kirde1962) theory for two-dimensional vortex sheet roll-up at an edge. This theory gives the vorticity in terms of the hypergeometric function. For thicker rings, there is a calculation of their formation from axisymmetric vortex sheet roll-up (Nitsche & Krasny Reference Nitsche and Krasny1994), but no corresponding prediction of the vorticity distribution. We hope that this lack is addressed soon.
1.2 Present work
The vorticity equation is
$\text{D}\unicode[STIX]{x1D701}/\text{D}t=0$
for inviscid swirl-free axisymmetric flow. Hence, if the flow is steady in some uniformly propagating frame, we have
where
$\unicode[STIX]{x1D713}$
is the Stokes streamfunction in the ring frame. It is natural to wonder if viscous vortex rings generated in the laboratory obey (1.1) and, if so, what form
$F(\unicode[STIX]{x1D713})$
takes. Knowing this, one could solve for the main structure of the ring (i.e. apart from the wake and edge layer) by solving an elliptic free boundary problem (Eydeland & Turkington Reference Eydeland and Turkington1988). It was this possibility that motivated the present work.
The axisymmetric Navier–Stokes simulations performed here show that, throughout the vortex ring bubble (except very close to the dividing streamline), a state of
$\unicode[STIX]{x1D701}\approx F(\unicode[STIX]{x1D713},t)$
is reached after a period of relaxation. However, even after moderately large times,
$F(\unicode[STIX]{x1D713},t)$
is not a universal function and depends on Reynolds number,
$Re$
, and stroke-to-diameter ratio,
$L/D$
. The best that can be said is that
$F(\unicode[STIX]{x1D713},t)$
tends to an approximately linear function with a non-zero intercept, whose dependencies will be studied in the following.
For interpreting results, it is helpful to keep in mind some consequences of the vorticity equation. We adopt a reference frame moving with the ring at speed
$U(t)$
. Because the frame is decelerating, a fictitious force term,
$-\unicode[STIX]{x1D70C}\dot{U}\widehat{\boldsymbol{x}}$
, is present on the right-hand side of the equation for momentum per unit volume. However, for uniform density,
$\unicode[STIX]{x1D70C}$
, this term can be absorbed into the pressure,
$p\rightarrow p+\unicode[STIX]{x1D70C}\dot{U}x$
. Thus, the vorticity equation remains unchanged in the moving frame:
For axisymmetric swirl-free flow in cylindrical coordinates (Batchelor Reference Batchelor1967, p. 602), equation (1.2) becomes
where
define the Stokes streamfunction
$\unicode[STIX]{x1D713}$
. Note that for inviscid flow
$\unicode[STIX]{x1D701}\equiv \unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}/r$
is conserved following fluid elements; the
$r$
in the denominator accounts for vorticity stretching. Suppose that at a given instant we have
$\unicode[STIX]{x1D701}=F(\unicode[STIX]{x1D713},t)$
in a certain region of the flow. The simulations will show that this is approximately the case within the vortex ring bubble. Then by direct substitution into (1.3) one sees that the advection terms sum to zero, which is called advective balance. This does not imply that the subsequent evolution is viscous. The diffusion operator in (1.3) destroys advective balance when the ring is not thin; see paragraph (
$a$
) in § 8. It is suggested that maintenance of the condition
$\unicode[STIX]{x1D701}\approx F(\unicode[STIX]{x1D713},t)$
requires that nonlinear terms still be active. A solution for the time evolution of
$F(\unicode[STIX]{x1D713},t)$
must take advective effects into account and paragraph (
$a$
) in § 8 conjectures that this effect is shear dispersion, which leads to averaging of the vorticity field on a closed streamline. A similar suggestion was made by Rhines & Young (Reference Rhines and Young1983) for planar vortices.
1.3 Previous efforts
Berezovski & Kaplanski (Reference Berezovski and Kaplanski1987) obtained a time-dependent Stokes-flow vortex ring solution with a peaked vorticity distribution. This solution has zero thickness at
$\unicode[STIX]{x1D708}t=0$
and tends to the self-similar solution obtained by Phillips (Reference Phillips1956) as
$\unicode[STIX]{x1D708}t\rightarrow \infty$
. No advection effects are present aside from a spatially uniform drift, which was obtained by Kaplanski & Rudi (Reference Kaplanski and Rudi2001) using the procedures of Kambe & Oshima (Reference Kambe and Oshima1975) and Rott & Cantwell (Reference Rott and Cantwell1993a
,Reference Rott and Cantwell
b
). Kaplanski & Rudi (Reference Kaplanski and Rudi2005) obtain the value of
$\unicode[STIX]{x1D708}t$
at which the non-dimensional energy
$\widetilde{E}$
of the Stokes solution equals the value 0.33 of the experimental limiting ring. Interestingly, at this value of
$\unicode[STIX]{x1D708}t$
, the non-dimensional circulation and ring speed also roughly match experimental values. This suggests that the Stokes-flow solution does capture some of the overall properties of high-Reynolds-number rings.
A virtue of the Stokes solution is that it fully includes curvature effects for vorticity diffusion. However, some issues arise because vorticity advection is absent. (i) Since the vorticity field is not deformed by the curvature-induced strain, vorticity contours are not axially elongated as in actual high-Reynolds-number rings. To overcome this difficulty, Kaplanski, Fukumoto & Rudi (Reference Kaplanski, Fukumoto and Rudi2012) introduce axial elongation and radial compression factors into the vorticity field of the Stokes solution. The values of the two factors are selected by matching the non-dimensional energy and circulation of the model ring against values obtained from Navier–Stokes simulations. (ii) An edge layer, where advection and viscous terms nearly balance, is missing. (iii) Since viscosity enters the vorticity field of the Stokes solution only via the product
$\unicode[STIX]{x1D708}t$
, the solution cannot exhibit the extra Reynolds-number dependence observed in the present results. (iv) Since the viscous term acts to destroy advective balance, the Stokes ring cannot in general satisfy advective balance except at early times when the ring is thin.
Among other results, Fukumoto & Moffatt (Reference Fukumoto and Moffatt2000) obtained (from the Navier–Stokes equations) an asymptotic solution for the vorticity and streamfunction of a diffusing vortex ring to second order in the slenderness ratio. The solution was stated in terms of ordinary differential equations with respect to distance from the core. At second order, both vorticity and streamfunction contours become elliptical. In the future, it would be worthwhile to probe this solution from the point of view of advective balance.
Finally, we note the work of Couder & Basdevant (Reference Couder and Basdevant1986), who studied the generation of vortex pairs by two-dimensional von Kármán wakes. They briefly analysed the functional form of the vorticity in terms of the streamfunction by using a scatter plot, and found it to be linear, which corresponds to the Lamb dipole (Batchelor Reference Batchelor1967, p. 535).
2 Simulation set-up and post-processing
The axisymmetric and incompressible Navier–Stokes equations were solved using the commercial Fluent package. Advection terms were discretized using the quadratic upstream interpolation (QUICK) scheme of Leonard (Reference Leonard1979), which obtains values at cell faces using a linear combination of central and upwind interpolation. The iterative SIMPLE algorithm was selected to obtain, at the end of a time step, a consistent pressure–velocity pair that simultaneously satisfies the space–time discrete momentum equation and discrete mass conservation (see Ferziger & Perić Reference Ferziger and Perić2002). The implicit Euler scheme with linearization of the nonlinear term was selected for time discretization. Dissipative schemes require a grid refinement study to ensure reliability of the results, and this is described in appendix A.
Figure 1 shows the computational domain. It consisted of a pipe of diameter
$D$
and length
$3D$
. The number of grid points was
$N_{x}\times N_{r}=1103\times 251$
. The grid was non-uniform, with the highest density near the tube wall and exit plane (
$x=0$
). The region
$r/D\leqslant 1$
and
$-0.5\leqslant x/D\leqslant 8$
, which we call the ‘vortex region’, contained nearly all of the vorticity during the simulations and was chosen to have the highest resolution. A no-slip condition was applied at the pipe wall, and piston motion was simulated using a spatially uniform inlet velocity,
$U_{p}(t)$
, applied for
$t\in [0,T]$
. The form of
$U_{p}(t)$
was a trapezoidal pulse: the first and last 10 % of the pulse consisted of a linear acceleration and deceleration, respectively. The maximum piston velocity is denoted
$U_{0}$
. The resulting jet slug length was
A uniform pressure was specified at the outer boundaries. The Reynolds number is defined as
Computational domain. Not to scale.

The Stokes streamfunction was determined by solving
in the frame of reference moving with the ring. The solution for (2.3) was obtained in the neighbourhood of the vortex ring by interpolating
$\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}$
onto a regular grid and solving using a finite-difference scheme with second-order truncation error. To recast the velocity field (and hence the solution of (2.3)) in the frame of reference moving with the ring, the time-varying ring velocity was determined from the derivative of the axial position of the ring centroid, defined as the centroid of vorticity greater than 60 % of the peak vorticity. The time derivative was computed using a fourth-order finite-difference formula and smoothed using a Butterworth low-pass filter with a cutoff frequency of 10 % of the sample rate of the position data.
The procedure used to obtain
$\unicode[STIX]{x1D701}=F(\unicode[STIX]{x1D713},t)$
was to consider 100 level curves of
$\unicode[STIX]{x1D713}$
equispaced in the interval
$0\leqslant \unicode[STIX]{x1D713}\leqslant \unicode[STIX]{x1D713}_{max}$
. On each level curve, values of
$\unicode[STIX]{x1D701}$
at 100 points (equispaced in the angular direction and centred at the location of the peak in
$\unicode[STIX]{x1D701}$
) were obtained by spline interpolation from the solution grid. Mean and standard deviations of
$\unicode[STIX]{x1D701}$
on each curve were then calculated. This procedure, unlike the method of making a scatter plot of
$\unicode[STIX]{x1D701}$
versus
$\unicode[STIX]{x1D713}$
previously used in the literature, facilitates further analysis.
3 Overview and orientation
3.1 Overview
Simulation cases were:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}Re=1000\quad \text{and}\quad L/D\in \{1,4\};\\ Re=2000\quad \text{and}\quad L/D\in \{0.5,1,2,3,4,5,6\};\\ Re=4000\quad \text{and}\quad L/D\in \{1\}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
Quantities normalized by piston speed
$U_{0}$
and pipe diameter
$D$
are denoted by an asterisk:
A hat is used for
$\unicode[STIX]{x1D701}$
and
$\unicode[STIX]{x1D713}$
normalized by their peak values (which decay in time):
To diagnose the shape of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}})$
at different instants, coefficients of the least-squares cubic
fitted to the data are obtained. Only the mean value of
$\widehat{\unicode[STIX]{x1D701}}$
on each streamline was used for the fit. To avoid the edge layer, the last value included in the fit is
$\widehat{\unicode[STIX]{x1D713}}=0.99$
. To ensure that
$\unicode[STIX]{x1D701}=F(\unicode[STIX]{x1D713},t)$
to a good approximation, results of the fit and any quantity derived from it are plotted only if the standard deviation of
$\widehat{\unicode[STIX]{x1D701}}$
on every sampled streamline is less than
$0.04$
for
$0\leqslant \widehat{\unicode[STIX]{x1D713}}\leqslant 0.99$
, and the root-mean-square (r.m.s.) error of the cubic fit is less than
$0.005$
. After a referee’s comment, a least-squares fit was also performed to the form
In only a small fraction of the fields (0 %–16 % for any given case) was it found to give a smaller r.m.s. error than the cubic.
A diagnostic of interest will be the value of
$\widehat{\unicode[STIX]{x1D701}}$
as the dividing streamline is approached from the interior of the ring bubble:
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}}\rightarrow 1^{-})$
. The use of a limit acknowledges the presence of an edge layer governed by boundary-layer-type behaviour, which should match the interior flow. The limit should be thought of in the same way that one thinks of the limit of potential flow past a body as the surface is approached, which should match the outer limit of the boundary layer. The limit is obtained from the computations as
which represents a slight extrapolation since the last value included in the fit is
$\widehat{\unicode[STIX]{x1D713}}=0.99$
.
3.2 Orientation
(a) Contours of
$\unicode[STIX]{x1D701}^{\ast }$
at
$t^{\ast }=19.6$
(
$Re=2000,~L/D=4$
). White line:
$\unicode[STIX]{x1D713}=0$
contour in a frame moving with the ring bubble. (b) The corresponding
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}})$
function in the region enclosed by the dividing streamline. The width of the error bars is the standard deviation of
$\widehat{\unicode[STIX]{x1D701}}$
on each streamline and indicates the degree of advective balance on that streamline.

To orient the reader to typical flow behaviour, figure 2(a) shows contours of
$\unicode[STIX]{x1D701}^{\ast }$
for one case at a single instant. As is well known (Gharib et al.
Reference Gharib, Rambod and Shariff1998), for stroke ratios
$L/D$
below a limiting value
$(L/D)_{lim}$
, a single ring forms without leaving a trailing remnant. For
$L/D>(L/D)_{lim}$
, the ring cannot grow fast enough to accommodate all of the vorticity fed into it and the leading ring pinches off from the feeding layer, leaving behind a remnant trailing jet. In the present simulations
$(L/D)_{lim}$
is between 3 and 4. The stroke ratio,
$L/D=4$
, for the case shown in figure 2 is a little larger than
$(L/D)_{lim}$
and so there is a trailing remnant. The white line is the dividing streamline, the contour of
$\unicode[STIX]{x1D713}=0$
in the frame moving with the ring bubble. It separates fluid that is instantaneously moving with the ring from that flowing past. At the dividing streamline, there exists a viscous edge layer where outwardly diffusing vorticity is partly swept into the wake, and partly re-enters the vortex ring bubble. The viscous edge layer also includes the region within the bubble next to the symmetry axis. It is on the symmetry axis that total circulation of the upper half-plane is lost by vorticity diffusion.
Figure 2(b) shows the corresponding function
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}})$
inside the ring bubble. The width of the error bars is the standard deviation of
$\widehat{\unicode[STIX]{x1D701}}$
on each streamline; smaller values imply better advective balance. In this and future plots, the last plotted point is at
$\widehat{\unicode[STIX]{x1D713}}=0.99$
in order to avoid the edge layer where
$\widehat{\unicode[STIX]{x1D701}}\neq F(\widehat{\unicode[STIX]{x1D713}})$
.
4 Model: one-dimensional radial diffusion with an edge layer
A one-dimensional toy model based on vorticity diffusion in a cylindrical (and non-curved) vortex tube is presented here. It was developed to explain the extra Reynolds-number dependence observed for the vorticity on the dividing stream-surface, specifically, for the ratio
$\unicode[STIX]{x1D701}(\unicode[STIX]{x1D713}\rightarrow 0)/\unicode[STIX]{x1D701}_{max}(t)$
. The key ingredient in the model is a condition at the boundary of the vortex tube which accounts for the viscous layer at the dividing stream-surface of actual vortex rings. The model assumes that the edge layer is circular, which is not true even for thin rings. Hence the model is not a mathematical solution in any asymptotic limit and should be viewed as being merely illustrative.
For a thin vortex ring, i.e. one whose core radius
$\unicode[STIX]{x1D6FF}_{core}\ll R$
(the toroidal radius), variations in
$r$
across the cross-section of the tube can be neglected, that is
$r=R(1+O(\unicode[STIX]{x1D6FF}_{core}/R))$
. The vorticity dynamics is then locally (i.e. within the tube and in the co-moving frame) the same as for planar flow. Therefore, if the vortex initially has concentric circular streamlines and vorticity contours, they will continue to remain so. Therefore, it is reiterated that the model does not account for the effect of curvature-induced strain on the diffusion of vorticity in the interior of the vortex ring bubble.
A feature of circular streamline flow is that the advective terms are trivially zero and the azimuthal vorticity
$\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}$
obeys the planar diffusion equation
where
$s$
is radial distance measured from the centre of the core in a meridional plane, and
$s_{DS}(t)$
is the distance to the boundary of the vortex tube; it is allowed to grow in time as the ring expands. Equation (4.1) is solved numerically.
Next comes an assumption that makes the treatment an illustrative model (even for a thin ring) rather than a rigorous analysis: equation (4.1) will be applied all the way out to the dividing streamline (denoted using subscript DS), at which location neither local two-dimensionality holds nor are streamlines circular. The distance from the core centre to the notional dividing streamline is denoted as
$s_{DS}(t)$
, which represents roughly the average distance from the centre of the ring to the actual dividing streamline.
At the dividing streamline of an actual ring, there is an edge layer whose thickness
$\unicode[STIX]{x1D6FF}_{layer}$
(averaged along the length of the streamline, say) is
where
is the characteristic strain rate to which the layer is subject. The boundary condition applied to the diffusion equation (4.1) to model the edge layer is
where
$C_{model}$
is a modelling constant that allows one to change all
$\propto$
into
$=$
signs. Combining (4.2) and (4.3) one sees that
involves the instantaneous ring Reynolds number
$\unicode[STIX]{x1D6E4}(t)/\unicode[STIX]{x1D708}$
, a fact that will be important in the sequel.
The asymptotic analysis of Fukumoto & Moffatt (Reference Fukumoto and Moffatt2000) revealed that, to third order in core thickness, the radius of a diffusing ring grows linearly in time:
The rate of growth is different for the radius corresponding to the peak vorticity and the radius corresponding to the peak value of the co-moving streamfunction, the respective values of
$C_{FM}$
being 4.59 and 2.59. (Note that this difference implies loss of advective balance.) Tests with the model revealed that the value of
$C_{FM}$
influenced the decay rate of
$\unicode[STIX]{x1D713}_{max}^{\ast }(t)$
plotted in figure 4(e). We chose
$C_{FM}=2.9$
to obtain the best agreement between the model and simulation for
$\unicode[STIX]{x1D713}_{max}^{\ast }(t)$
. Ring expansion is not critical to the model and was included only to improve the postdiction of
$\unicode[STIX]{x1D713}_{max}^{\ast }(t)$
.
We set
$C_{model}=2.2$
and
$s_{DS}(t)/R(t)=0.50$
based on a visual matching of model and simulation results for the
$L/D=1,Re=1000$
case. The boundary condition (4.4) is of Robin type because it involves a linear combination of value and derivative; it is analogous to Newton’s law of cooling, which is used to represent a layer of advective cooling at the surface of a heat-conducting solid. Ring vorticity is undergoing a similar advective–diffusive process in the edge layer. Note that the boundary condition (4.4) is nonlinear in the vorticity due to the presence of
$\unicode[STIX]{x1D6E4}(t)$
, which is a functional of the vorticity.
The initial vorticity profile for the diffusion equation in terms of piston parameters is prescribed in appendix B. This initial profile, which also depends on
$\unicode[STIX]{x1D708}$
, uses the theory of self-similar vortex sheet roll-up and is valid for thin cores.
To obtain the streamfunction as a diagnostic and to plot profiles of
$\unicode[STIX]{x1D701}(\unicode[STIX]{x1D713})$
, we first calculate
$\unicode[STIX]{x1D701}\equiv \unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}/r$
, which
${\approx}\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}/R$
to consistent order. The circumferential velocity,
$u_{\unicode[STIX]{x1D703}}(s)$
, is then obtained from
The co-moving two-dimensional streamfunction,
$\unicode[STIX]{x1D713}_{2D}$
, is obtained by integrating
from which, finally, the Stokes streamfunction is
$\unicode[STIX]{x1D713}=R\unicode[STIX]{x1D713}_{2D}+C_{1}$
to consistent order. The integration constant
$C_{1}$
is chosen to make
$\unicode[STIX]{x1D713}=0$
at the notional dividing streamline
$s=s_{DS}$
.
5
$\unicode[STIX]{x1D701}(\unicode[STIX]{x1D713},t)$
for a thin Gaussian ring
A referee requested comparison of vorticity profiles from the simulation against those for a thin Gaussian ring (Saffman Reference Saffman1970). This solution has circular vorticity contours, which is true in the simulations only for small
$L/D$
at early times. Nevertheless, because it is a simple analytical model for a diffusing ring, its
$\unicode[STIX]{x1D701}(\unicode[STIX]{x1D713})$
profile provides a useful reference even outside its range of validity. It should be noted that, while the experiments of Cater, Soria & Lim (Reference Cater, Soria and Lim2004) obtained Gaussian-like vorticity profiles, they were thinner in the radial direction than the axial direction due to the effect of strain.
The relations necessary for plotting
$\unicode[STIX]{x1D701}(\unicode[STIX]{x1D713},t)$
for the Gaussian ring are presented below. Note that in making a comparison with simulations, the relations will be evaluated all the way to the dividing streamline which is outside their range of validity,
$s\ll R$
, even for thin rings (as before,
$s$
is the distance from the centre of the core in a meridional plane). The ring radius (
$R$
) and circulation (
$\unicode[STIX]{x1D6E4}$
) needed to evaluate the relations are, as for the model described in the previous section, obtained from the expressions in appendix B.
When (the core radius)
$\unicode[STIX]{x1D70E}(t)\ll R$
, the dynamics are locally two-dimensional and the two-dimensional Oseen (Gaussian) vortex is a valid local solution for
$\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}$
, and therefore
The speed of such a ring, needed later, was calculated by Saffman (Reference Saffman1970) as
Integrating the circumferential velocity of Oseen’s vortex gives its streamfunction:
where
$E_{1}$
is the exponential integral and the integration constant was chosen to non-dimensionalize the argument of the logarithm. The local Stokes streamfunction for the vortex ring is
The constant
$C_{2}$
is not disposable and is obtained by using the asymptotic matching rule described by Moore (Reference Moore1980) in his calculation of the speed of an elliptical core ring. The rule is that, in the region
$\unicode[STIX]{x1D70E}\ll s\ll R$
, equation (5.4) should match the streamfunction for a zero-thickness ring in the region
$s\ll R$
, which is given by Moore’s equation (2.27):
Matching gives
Equations (5.1) and (5.4) implicitly give
$\unicode[STIX]{x1D701}(\unicode[STIX]{x1D713})$
for the Gaussian ring. The slope of this function can be obtained explicitly at the core centre:
at
$s=0$
.
The behaviour of the core size is specified as
where
$\unicode[STIX]{x1D70E}_{0}$
is the core size at the end of the piston stroke. To obtain
$\unicode[STIX]{x1D70E}_{0}$
, we note that for the hypergeometric profile ((B 1) in appendix B) the radius of peak velocity is (Saffman Reference Saffman1978)
at the end of the piston stroke. On the other hand, the radius of peak velocity for the Gaussian vortex is
Equating (5.9) and (5.10) gives
Since the circulation of the Gaussian ring is assumed to be constant, growth of its radius would imply increasing impulse
$\propto \unicode[STIX]{x1D6E4}R^{2}$
. Hence its radius is kept fixed at its initial value.
6 Results
6.1 Reynolds-number dependence at
$L/D=1$
Profiles of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}})$
for
$L/D=1$
and three different Reynolds numbers increasing downwards. The simulation results are shown using error bars; the height of the error bars equals the standard deviation of
$\widehat{\unicode[STIX]{x1D701}}$
on each streamline. The left panels compare simulations with the model (lines). The right panels compare simulations with the Gaussian ring.
Figure 3 shows profiles of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}},t)$
for
$L/D=1$
and three Reynolds numbers. The left panels compare simulation results (error bars) with the model (lines); the height of the error bars equals the standard deviation of
$\widehat{\unicode[STIX]{x1D701}}$
on each streamline. The right panels compare simulations with the Gaussian ring. A small value for the standard deviation means that the vorticity is in approximate advective balance; this is achieved earlier and better for higher
$Re$
. Vorticity diffusion causes the profiles to become less curved with time, which happens slower with increasing
$Re$
, as expected. It is noteworthy that the value of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}}\rightarrow 1^{-})$
, i.e. as the dividing streamline is approached and ignoring the edge layer, is non-zero. The value of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}}\rightarrow 1^{-})$
increases with time and decreases with Reynolds number. At early times, the profiles for the model and Gaussian ring do not have as much vorticity in the middle region of the bubble (
$\widehat{\unicode[STIX]{x1D713}}\approx 0.5$
) as the simulation does; this reflects inaccuracy in the initial condition provided by the vortex sheet roll-up theory (appendix B). For the highest
$Re$
case, a small secondary vorticity peak, which represents an outer spiral turn, is observed in this region. The distance between adjacent spiral turns is larger than the diffusion length in this region (Pullin Reference Pullin1979).
The model initially lags the simulation in its development towards an approximately linear profile and then leads at later times. Overall, the Gaussian ring gives a similar degree of agreement with the simulation as the model does, its most obvious error being a significant overshoot in the value of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}}\rightarrow 1^{-})$
at late times. This occurs because there is no edge layer to sweep away diffused vorticity.
Reynolds-number dependence of profile parameters with time scaled by
$Re$
(
$L/D=1$
).

Figure 4 shows the evolution of various parameters of the vorticity profile with respect to
$t^{\ast }/Re=\unicode[STIX]{x1D708}t/D^{2}$
. Figure 4(a) displays the evolution of the cubic coefficients, which are plotted only after the r.m.s. error (figure 4
b) of the cubic fit becomes less than
$0.005$
. With time left unscaled, there was considerable difference in the curves (not shown) for the three
$Re$
cases. Scaling time by
$Re$
reduces this difference due to the important role of viscosity in the overall evolution of
$\unicode[STIX]{x1D701}(\unicode[STIX]{x1D713},t)$
. Some extra
$Re$
dependence of the coefficients remains, and this must ultimately be accounted for by Reynolds-number effects in (i) the ring formation process, (ii) the boundary condition at the edge layer, and (iii) the competition between viscous destruction of advective balance and its restoration by nonlinear terms.
As
$t^{\ast }/Re$
increases, the linear coefficient
$a_{1}$
tends to approximately unity, while
$a_{2}$
and
$a_{3}$
tend to approximately zero. This means that the profile becomes approximately linear with a small intercept equal to
$(1-a_{1})+a_{2}+a_{3}$
at the dividing stream-surface,
$\widehat{\unicode[STIX]{x1D713}}\rightarrow 1^{-}$
.
Figure 4(c) shows that the decay of the peak vorticity (
$\unicode[STIX]{x1D701}_{max}^{\ast }$
) for the three
$Re$
cases curves collapses quite well (but not perfectly) when time is scaled using
$Re$
. The three violet-coloured curves (which are nearly coincident) show that the one-dimensional model tracks the simulation results quite well; the maximum error is approximately 10 %.
Figure 4(d) shows that the value of
$\widehat{\unicode[STIX]{x1D701}}$
as the dividing streamline is approached has a significant extra Reynolds-number dependence. It was the need to explain this that prompted development of the model with an edge layer, whose results (violet lines) give the correct qualitative behaviour.
In the context of the model, the extra
$Re$
dependence could arise from the presence of
$\unicode[STIX]{x1D708}$
in both the initial condition (appendix B) and boundary condition (4.4). To ascertain which of the two was more important, control runs with the model were performed for the three Reynolds numbers in which the initial profile was fixed at the form given by the model for
$Re=2000$
. This resulted in only a very slight change compared to the model results in figure 4(d), indicating that most of the extra
$Re$
dependence comes from the edge-layer boundary condition rather than from changes in the initial profile. The Gaussian ring (cyan lines) exhibits an insignificant amount of extra
$Re$
dependence that is very different from that observed in the simulations and the model. This is because the Gaussian ring lacks an advective–diffusive edge layer and because
$\unicode[STIX]{x1D708}$
enters the solution (mostly) via the product
$\unicode[STIX]{x1D708}t$
; the small extra
$Re$
dependence is due to the presence of
$4\unicode[STIX]{x1D708}T$
in equation (5.11) for the initial core size. The behaviour of
$\widehat{\unicode[STIX]{x1D701}}(t)$
for the Gaussian ring at the dividing streamline never goes through an inflection point and therefore has the opposite sign for the curvature at large times.
Figure 4(e) shows the decay of the peak value,
$\unicode[STIX]{x1D713}_{max}^{\ast }$
, of the streamfunction. The best straight lines were obtained with a logarithmic abscissa. The model agrees very well with the simulations for
$Re=4000$
while having a small constant error for the other two Reynolds numbers. The small extra Reynolds-number dependence in
$\unicode[STIX]{x1D713}_{max}^{\ast }$
is not predicted by the model.
Figure 4(f) assesses the deviation from advective balance in the simulations using the standard deviation of
$\widehat{\unicode[STIX]{x1D701}}$
on each streamline averaged for
$0\leqslant \widehat{\unicode[STIX]{x1D713}}\leqslant 0.99$
. Note that the abscissa is measured from the end of the piston stroke at
$t=T$
. At higher
$Re$
, advective balance is achieved to a greater degree and is destroyed more slowly with respect to
$t^{\ast }/Re$
. In figure 4(f) the rate of destruction of advective balance with respect to
$t^{\ast }/Re$
is roughly proportional to
$Re^{-1}$
; hence the rate with respect to
$t^{\ast }$
is roughly proportional to
$Re^{-2}$
. The quadratic-in-
$\unicode[STIX]{x1D708}$
rate of destruction suggests that the advective term plays a role in restoring advective balance. One also observes that advective balance reaches the 0.02 threshold earlier (with respect to
$t^{\ast }/Re$
) with increasing
$Re$
. The increasing level of advective imbalance is likely a result of decreasing
$\unicode[STIX]{x1D6E4}(t)/\unicode[STIX]{x1D708}$
with time.
6.2
$L/D$
dependence at
$Re=2000$
Profiles of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}})$
showing the dependence on
$L/D$
(
$Re=2000$
).

Recall that the limiting value of
$L/D$
beyond which the ring trails a jet or sheds some vorticity during the formation process is
$3<(L/D)_{lim}<4$
in the simulations. Model results are presented for only
$L/D\leqslant 3$
: according to current understanding, there should be no
$L/D$
dependence for
$L/D\geqslant (L/D)_{lim}$
. The simulations will indicate that this is not the case for
$\unicode[STIX]{x1D713}_{max}^{\ast }$
.
Figure 5 shows
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}},t)$
profiles. The error bars indicate that for increasing
$L/D$
(thicker cores) it takes longer for profiles to reach advective balance near the dividing streamline. At late times, however, the smallest
$L/D$
case (
$L/D=0.5$
) is the least advectively balanced overall. We shall see below that this is because it has the smallest ring Reynolds number
$\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D708}$
. As observed in the Reynolds-number study presented earlier (figure 3), the development of model profiles initially lags and later leads the simulations. Apart from the overshoot at late times, the Gaussian ring profiles (figure 5
b,d) agree with the simulations about as well as the model does.
Figure 5(f) compares simulation profiles for the different
$L/D$
cases at a fixed late time of
$t^{\ast }=40$
. The profiles for
$L/D\geqslant 2$
are practically coincident. The
$L/D=1$
(dotted) case is a little more diffuse, and the
$L/D=0.5$
(solid) case somewhat more diffuse than the rest. This is because smaller
$L/D$
rings have a smaller
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
and are therefore more diffusive.
The
$L/D$
dependence of various parameters of the vorticity profile (
$Re=2000$
).

The cubic coefficients in figure 6(a) show that differences in the shape of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}})$
at early
$t^{\ast }/Re$
diminish later and the shape converges to a form that is weakly dependent on
$L/D$
; some variation and crossing of the curves at large times should be noted. The curve for
$a_{1}$
and
$L/D=0.5$
(solid red line) has a noticeable shift relative to the other curves, perhaps due to its low ring Reynolds number
$\unicode[STIX]{x1D6E4}(t)/\unicode[STIX]{x1D708}$
.
The convergence noted in the previous paragraph is reminiscent of the study by Stanaway, Cantwell & Spalart (Reference Stanaway, Cantwell and Spalart1988) of viscous vortex rings, which found (pp. 74–76 in their work) that, for different initial core thicknesses (but the same
$\unicode[STIX]{x1D6E4}_{0},\unicode[STIX]{x1D708}$
and
$R_{0}$
), plots of normalized ring speed versus
$\unicode[STIX]{x1D708}t/R_{0}^{2}$
all converged to approximately the same curve. In their plots,
$t$
is measured from a virtual origin when the ring had zero core thickness.
Figure 6(c) displays the decay in peak vorticity,
$\unicode[STIX]{x1D701}_{max}^{\ast }(t)$
. The model works best for
$L/D=1$
and has moderately large errors for the other
$L/D$
cases. These errors are present in the initial condition and this suggests improvement of the relations in appendix B.
Figure 6(d) shows
$\widehat{\unicode[STIX]{x1D701}}$
at the dividing streamline. At early times, a lower
$L/D$
produces a smaller value of this quantity. This is expected since smaller
$L/D$
leads to a thinner core. However, a cross-over occurs and eventually smaller
$L/D$
produces larger values of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}}\rightarrow 1^{-})$
with unabated growth at the final instant of the simulation. This is non-intuitive at first sight, and it is hypothesized that it is mainly due to smaller
$L/D$
rings having a smaller initial ring Reynolds number
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
. Since the trend of the dependence on
$L/D$
at large times is reproduced by the model, we can use it to test the hypothesis. We re-ran three
$L/D$
cases (
$L/D=0.5,~1$
and 2) using the model, but set the initial circulation for all cases to the value for
$L/D=2$
. Doing so reduced the range of variation in
$\widehat{\unicode[STIX]{x1D701}}_{model}(\widehat{\unicode[STIX]{x1D713}}=1,t^{\ast }/Re=0.025)$
from
$[0.091,0.152]$
to
$[0.091,0.097]$
. Hence a substantial portion of the variation with respect to
$L/D$
is due to changing
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
.
Next we hypothesize that this effect of
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
(in the context of the model) enters mainly through the presence of
$\unicode[STIX]{x1D6E4}(t)/\unicode[STIX]{x1D708}$
in the edge-layer boundary condition (4.5). To verify this, a test was performed in which
$\unicode[STIX]{x1D6E4}(t)$
in the model boundary condition was kept fixed in time at the
$\unicode[STIX]{x1D6E4}_{0}$
value for
$L/D=2$
. In this case the range of variation of
$\widehat{\unicode[STIX]{x1D701}}_{model}(\widehat{\unicode[STIX]{x1D713}}=1,t^{\ast }/Re=0.025)$
for the three
$L/D$
cases was very small, namely,
$[0.076,0.078]$
. We conclude from this that a substantial portion of the
$L/D$
dependence of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}}\rightarrow 1^{-})$
is due to the dependence of edge-layer physics on
$\unicode[STIX]{x1D6E4}(t)/\unicode[STIX]{x1D708}$
, where
$\unicode[STIX]{x1D6E4}_{0}$
depends on
$L/D$
.
To avoid clutter, results for the Gaussian ring have not been included in figure 6. Suffice it to say, it does not provide better predictions than the model.
Reynolds-number dependence at
$L/D=4$
of
$\widehat{\unicode[STIX]{x1D701}}$
at the dividing streamline.

An important feature of the simulations which the model is unable to reproduce is that, for
$L/D\geqslant 3$
, the temporal growth of
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}}\rightarrow 1^{-})$
saturates to a constant value. This constant value has a weak dependence on
$L/D$
for
$L/D\geqslant 3$
. This is because
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
loses its dependence on
$L/D$
for
$L/D>(L/D)_{lim}$
. The flat value for
$L/D=4$
decreases with Reynolds number, as shown in figure 7. Although the model is inaccurate at this value of
$L/D$
, we hypothesize on the basis of previous model results that this dependence mainly reflects the change in
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
which affects the edge layer.
Figure 6(e) shows that the peak value of the streamfunction increases with
$L/D$
, a trend that is also captured by the model. However, the value of
$\unicode[STIX]{x1D713}_{max}^{\ast }$
for
$L/D=3$
is overpredicted by a fair amount.
According to current thinking, the characteristics of the leading ring experimentally produced when
$L/D\geqslant (L/D)_{lim}$
should be independent of
$L/D$
; for example, the circulation becomes independent of
$L/D$
(Gharib et al.
Reference Gharib, Rambod and Shariff1998). However, figure 6(e) shows that
$\unicode[STIX]{x1D713}_{max}^{\ast }$
continues to increase with
$L/D$
. Since
$\unicode[STIX]{x1D713}^{\ast }=0$
at the dividing streamline,
$\unicode[STIX]{x1D713}_{max}^{\ast }$
represents the volumetric flux of recirculating flow in the bubble. The mechanism that causes this quantity to increase for
$L/D>(L/D)_{lim}$
should be investigated in the future.
Figure 6(f) assesses the degree of advective balance as a function of
$L/D$
. We observe that cases with smaller
$L/D$
relax faster. This is consistent with intuition: thinner rings have most of their vorticity in a region where local two-dimensionality holds and planar flows with circular streamlines remain in exact inviscid balance as they diffuse. However, at later times, the smaller
$L/D$
cases lose advective balance faster. This is because smaller
$L/D$
gives smaller
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
. There is very little difference in the loss of advective balance for cases
$L/D\geqslant 4$
because all these cases have approximately the same
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
.
The notion that the
$L/D$
dependence is mostly due to changes in
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
can also be approximately tested using the simulations in the following way. Substituting the expressions for piston stroke (2.1) and slug circulation (B 7) into the empirical equation (B 6) for
$\unicode[STIX]{x1D6E4}_{0}$
gives
where
$Re_{L}\equiv U_{0}L/\unicode[STIX]{x1D708}$
is the Reynolds number based on piston stroke. Equation (6.1) shows that matching
$Re_{L}$
approximately matches
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
. In particular, the two cases
$(Re=2000,L/D=0.5)$
and (
$Re=1000,L/D=1$
) have the same
$Re_{L}=1000$
and (6.1) gives their corresponding
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
as 856 and 702, respectively. Simulation results for the two cases are shown as the dot-dashed and solid black lines in figure 8. Comparing them, one observes a reasonable correspondence in
$\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}}\rightarrow 1^{-})$
,
$\unicode[STIX]{x1D701}_{max}^{\ast }$
, and the rate of loss of advective balance. The quantity
$\unicode[STIX]{x1D713}_{max}^{\ast }$
is not shown; it cannot be expected to correspond in the two cases since it is affected by the ratio of core size to radius.
Comparison of two cases that have similar
$\unicode[STIX]{x1D6E4}_{0}/\unicode[STIX]{x1D708}$
.

7 Summary
-
(i) After a period of relaxation, a state of approximate advective balance,
$\unicode[STIX]{x1D701}\approx \unicode[STIX]{x1D701}(\unicode[STIX]{x1D713},t)$
, is achieved within the region bounded by the dividing streamline. -
(ii) The function
$\widehat{\unicode[STIX]{x1D701}}=\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}},t)$
(where hats denote quantities normalized by peak values) evolves in time from being convex to being approximately linear, but with a non-zero intercept as the dividing streamline is approached. -
(iii) Introduction of the viscous time,
$t^{\ast }/Re=\unicode[STIX]{x1D708}t/D^{2}$
, captures some but not all of the Reynolds-number dependence in the evolution of
$\widehat{\unicode[STIX]{x1D701}}=\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}},t)$
. In particular, the value of
$\widehat{\unicode[STIX]{x1D701}}$
at the dividing streamline exhibits a strong extra Reynolds-number dependence. To explain this, a toy vorticity diffusion model with an edge layer at the boundary was developed. The model indicates that the extra Reynolds-number dependence arises primarily from the role of
$\unicode[STIX]{x1D708}$
in determining the thickness of the edge layer. The dependence on
$L/D$
at later times arises mostly from the ring Reynolds number
$\unicode[STIX]{x1D6E4}(t)/\unicode[STIX]{x1D708}$
, which determines the thickness of the edge layer. The conventional Gaussian profile, in which
$\unicode[STIX]{x1D708}$
enters (mostly) through the product
$\unicode[STIX]{x1D708}t$
, does not exhibit the extra
$Re$
dependence. -
(iv) A diagnostic of the function
$\widehat{\unicode[STIX]{x1D701}}=\widehat{\unicode[STIX]{x1D701}}(\widehat{\unicode[STIX]{x1D713}},t)$
was provided by the coefficients of its cubic fit plotted versus
$t^{\ast }/Re$
. They indicate a lessening with time (
$t^{\ast }/Re$
) of the
$L/D$
and (extra)
$Re$
dependence. Nevertheless, a non-negligible dependence on both parameters continues to remain, most prominently in the value of
$\widehat{\unicode[STIX]{x1D701}}$
as the dividing streamline is approached, for as long as we have run the simulations.
8 Suggestions for future work
-
(i) An important remaining question is how advective balance is maintained in the face of a diffusion term that acts to destroy it. Suppose that at some instant we have
$\unicode[STIX]{x1D701}=F(\unicode[STIX]{x1D713},t)$
. Then upon using the chain rule relations (8.1a,b )
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}x}=-u_{r}r\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}},\quad \frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}r}=u_{x}r\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}, & & \displaystyle\end{eqnarray}$$
(8.2)
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}^{2}F}{\unicode[STIX]{x2202}x^{2}}=(u_{r}r)^{2}\frac{\unicode[STIX]{x2202}^{2}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{2}}-r\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}, & \displaystyle\end{eqnarray}$$
(8.3)together with
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}^{2}F}{\unicode[STIX]{x2202}r^{2}}=(u_{x}r)^{2}\frac{\unicode[STIX]{x2202}^{2}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{2}}+\frac{\unicode[STIX]{x2202}(u_{x}r)}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}, & \displaystyle\end{eqnarray}$$
(8.4)the vorticity equation (1.3) becomes
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D701}=\frac{\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}}{r}=\frac{1}{r}\left(\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}x}-\frac{\unicode[STIX]{x2202}u_{x}}{\unicode[STIX]{x2202}r}\right), & & \displaystyle\end{eqnarray}$$
(8.5)where
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}=r^{2}q^{2}\frac{\unicode[STIX]{x2202}^{2}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{2}}-(r^{2}F-4u_{x})\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}, & & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D708}t$
and
$q^{2}=u_{x}^{2}+u_{r}^{2}$
. The coefficients
$r^{2}q^{2},r^{2}$
and
$u_{x}$
of the three terms on the right-hand side of (8.5) are not functions of
$\unicode[STIX]{x1D713}$
only. This should be contrasted with the planar case, where streamlines are circular and the diffusion operator is a function of
$\unicode[STIX]{x1D713}_{2D}$
. We conclude that the diffusion term in axisymmetric flow acts to destroy advective balance. How is it then restored? A reasonable conjecture is that disturbances to the advectively balanced state induced by the diffusion term are sheared by differential rotation and eventually average out along each closed streamline. Note that this is an advective effect. It is analogous to averaging by Taylor (or shear) dispersion of a passive scalar around closed streamlines studied for planar flow by Rhines & Young (Reference Rhines and Young1983), who conjectured that the same mechanism could occur for vorticity. When the Reynolds number,
$\unicode[STIX]{x1D6E4}(t)/\unicode[STIX]{x1D708}$
, becomes sufficiently small, destruction of advective balance by the diffusion term will dominate its restoration by differential rotation.
-
(ii) If streamline averaging accurately describes the restoring effect of advection, then the following algorithm could be used to compute the evolution of
$\unicode[STIX]{x1D701}$
.(i) Diffuse
$\unicode[STIX]{x1D701}$
for a short time
$\unicode[STIX]{x1D6FF}t$
, imposing the boundary condition dictated by the edge layer.(ii) To mimic the effect of advection, average
$\unicode[STIX]{x1D701}$
along each closed streamline. Upon combining the two steps into one, equation (8.5) becomes (8.6)within the vortex ring bubble. Here
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}=\langle r^{2}q^{2}\rangle \frac{\unicode[STIX]{x2202}^{2}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{2}}-(\langle r^{2}\rangle F-4\langle u_{x}\rangle )\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}} & & \displaystyle\end{eqnarray}$$
$\langle \cdot \rangle$
denotes an average around a closed streamline.
(iii) The final step of the algorithm would be to compute the resulting velocity field and streamfunction from the Poisson equation.
-
(iii) The most non-trivial aspect of the above procedure is obtaining the boundary condition on
$F(\unicode[STIX]{x1D713})$
where the flow in the interior of the bubble matches with edge layers. The present work ignored the detailed structure of the edge layer and wake. They should be studied in the future. Hints can be obtained from the work of Harper & Moore (Reference Harper and Moore1968) on falling drops. An analytical or semi-analytical treatment that matches the interior flow, edge layers and exterior potential flow would be most satisfying. An even more satisfying analysis would also quantify the deviation from advective balance at every location around a closed streamline. One aspect of this deviation is the forward lean of vorticity contours seen in figure 2(a) (B. Cantwell, private communication). This fore–aft asymmetry may be related to the fact that the windward side of the edge layer experiences positive strain along the dividing streamline, while the opposite is true on the leeward side. -
(iv) The Prandtl–Batchelor theorem asserts that
$F(\unicode[STIX]{x1D713})=\text{const.}$
for exactly steady axisymmetric flow with closed streamlines (Batchelor Reference Batchelor1956, pp. 186–187), a result obtained by integrating the momentum equation around a closed streamline. The critical assumption is that the rate of change of circulation around every closed streamline is zero. Even though this is not true in the present case, Batchelor’s integral analysis might still be useful. Note that the recirculation region of a wake behind an axisymmetric bluff body or the interior of a falling drop (Harper & Moore Reference Harper and Moore1968) can remain steady, in spite of viscous diffusion, because of a continual source of vorticity at the dividing streamline. For a bluff body, this source is the separated boundary layer of the body, and for a drop it is the boundary condition of stress and velocity continuity at the two-fluid interface. In the present situation, there is a sink of vorticity at the dividing streamline.
Acknowledgements
K.S. is grateful to T. Hartlep and A. Wray for their suggestions during the internal review. We also thank the referees for their useful suggestions.
Appendix A. Mesh and time-step refinement
Mesh refinement parameters.

Table 1 provides the grid sizes used in the refinement study and table 2 provides the results. The 2
$\times$
grid has twice the resolution of the baseline grid inside the tube and in the vortex region. This region is axially and radially shorter than the computational domain in order to keep the number of mesh points manageable. The quantity compared in table 2 is
$\unicode[STIX]{x1D701}^{\ast }$
at the same
$100\times 100$
points used to obtain the
$\unicode[STIX]{x1D701}(\unicode[STIX]{x1D713})$
profiles. This comparison was performed at three values of
$t^{\ast }$
indicated in the tables.
Change (baseline relative to 2
$\times$
grid) in
$\unicode[STIX]{x1D701}^{\ast }(\unicode[STIX]{x1D713})$
due to mesh and time-step refinement. Values within parentheses are for time-step refinement (
$\unicode[STIX]{x0394}t^{\ast }=0.02$
to
$\unicode[STIX]{x0394}t^{\ast }=0.01$
) while values outside parentheses are for mesh refinement. The columns for relative difference show the mean of the relative absolute difference and the maximum of the relative absolute difference. Here
$Re=2000$
.

Table 2 shows that, when a grid refinement is performed, the maximum of the relative change in
$\unicode[STIX]{x1D701}^{\ast }$
is 2.71 %; it occurs near the dividing streamline where
$\unicode[STIX]{x1D701}^{\ast }$
is small. When the time step is halved from
$\unicode[STIX]{x0394}t^{\ast }=0.02$
to
$\unicode[STIX]{x0394}t^{\ast }=0.01$
, this quantity is 4.12 %; this value occurs far from the centre of the vortex where the last portion of trailing jet flow is still being entrained and is not representative of errors at later times.
Appendix B. Initial vorticity profile for the one-dimensional model
For a sharp edge and constant (in time) velocity upstream of the edge, the theory of self-similar roll-up (Kirde Reference Kirde1962) gives the following profile for the circumferential velocity in the vortex at the end of the piston stroke,
$t=T$
:
where
$s$
is cylindrical radial distance from the core centre,
$s_{0}$
is the outermost radius of the vortex spiral at
$t=T$
, and
$M(a,b;x)$
is the confluent hypergeometric function (Abramowitz & Stegun Reference Abramowitz and Stegun1965, p. 503), and the constant
$C_{0}$
is
The diffusion model of the paper requires specification of the vorticity corresponding to (B 1), which is
with
$\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}=0$
for
$s>s_{0}$
. Kirde’s theory does not consider the profile in the outer region of the core, denoted as region I in Pullin (Reference Pullin1979); in terminating the velocity profile (B 1) at
$s=s_{0}$
and setting
$\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D719}}=0$
for
$s>s_{0}$
, we are following equation (2.16) in Saffman (Reference Saffman1978).
To specify the radius,
$s_{0}$
, of the vortex spiral and the initial ring radius,
$R_{0}$
, we use the following expressions from Saffman (Reference Saffman1978):
which use scalings from the theory of self-similar roll-up and constants from fitting experiments. For the 0.28 constant, see the fourth line after Saffman’s equation (3.6); and for the 0.11 constant, see the second line on p. 633 of Saffman’s paper. To prevent the spiral from overfilling the notional dividing streamline, we set
$s_{0}=\min [(s_{0})_{Saffman},s_{DS}]$
.
The circumferential velocity at the core boundary is obtained from the circulation, i.e.
where for the circulation
$\unicode[STIX]{x1D6E4}_{0}$
we use
which is equation (2.6) from Shariff & Leonard (Reference Shariff and Leonard1992) and represents a fit to experiments. For the trapezoidal piston velocity prescribed in the present work, the slug circulation is









