1. Introduction
The dynamics of vortex rings has long been an intriguing subject for researchers studying fluid dynamics (Lamb Reference Lamb1932; Batchelor Reference Batchelor1967; Saffman Reference Saffman1992). Interest has been shown in establishing analytical models of an ideal vortex ring, such as those of Fraenkel (Reference Fraenkel1972) and Norbury (Reference Norbury1973). Experimental research on vortex rings is usually carried out by using the piston–cylinder apparatus. It consists of a piston moving through a tube, leading to the generation of an isolated starting vortex ring at the exit orifice (Dabiri & Gharib Reference Dabiri and Gharib2004). It is prescribed here that vortex rings with the same sense of rotation as this kind of starting vortex ring are positive vortex rings (PoVRs) and that vortex rings with the opposite sense of rotation are negative vortex rings (NeVRs).
Counter-rotating vortex ring pairs (VRPs) which consist of a pair of a PoVR and an NeVR are prevalent in biological unsteady locomotion. Flying insects leave a linked vortex chain formed by a series of combined vortex rings with opposite sense of rotation in their wakes (Brodsky Reference Brodsky1991). Such VRPs are also observed in high-resolution experiments at the ends of wing strokes in flying insects (see Dickinson et al. Reference Dickinson, Farley, Full, Koehl, Kram and Lehman2000; Thomas et al. Reference Thomas, Taylor, Srygley, Nudds and Bomphrey2004; Lehmann, Sane & Dickinson Reference Lehmann, Sane and Dickinson2005; Hubel et al. Reference Hubel, Riskin, Swartz and Breuer2010) and in the wakes of self-propelled swimming animals (see Dabiri et al. Reference Dabiri, Colin, Costello and Gharib2005; Drucker & Lauder Reference Drucker and Lauder2005; Bartol et al. Reference Bartol, Krueger, Stewart and Thompson2009; Stewart, Bartol & Krueger Reference Stewart, Bartol and Krueger2010). In addition, VRPs are also presented in high-resolution computational studies (for instance Aono, Liang & Liu Reference Aono, Liang and Liu2008; Sahin, Mohseni & Colin Reference Sahin, Mohseni and Colin2009; Hamlet, Santhanakrishnan & Miller Reference Hamlet, Santhanakrishnan and Miller2011).
Considering the widespread biological evidence above, it is essential to study the effect of a VRP on biological propulsion. These correlated vortex rings involve unique dynamical behaviours. Regarding vortex formation, some visualisation results have shown that far less circulation is shed in the upstrokes than in the downstrokes (Dickinson & Götz Reference Dickinson and Götz1996; Warrick, Tobalske & Powers Reference Warrick, Tobalske and Powers2005). The formation of an NeVR is constrained by the prior PoVR. In terms of wake dynamics, it has been conjectured that the correlated vortex rings cause an increased mass flow in the wake, thus generating a larger force (Srygley & Thomas Reference Srygley and Thomas2002; Dabiri et al. Reference Dabiri, Colin, Costello and Gharib2005).
Vortex ring pairs take on various vortical configurations, which embody the dynamics of biological propulsion, as locomotive force is generated from a change in the vortical impulse in the fluid. Different patterns of vortical configurations can exist in the wakes of animals, and the patterns may correspond to propulsion performance characterised by different features (high efficiency, rapid acceleration, etc.) (Dickinson, Lehmann & Sane Reference Dickinson, Lehmann and Sane1999; Lehmann Reference Lehmann2008; Stewart et al. Reference Stewart, Bartol and Krueger2010).
Experimental examples of NeVRs are generated by time-dependent kinematics of physical bodies (Dabiri et al. Reference Dabiri, Colin, Costello and Gharib2005; Weathers et al. Reference Weathers, Folie, Liu, Childress and Zhang2010; Liu et al. Reference Liu, Ristroph, Weathers, Childress and Zhang2012). The circulation of NeVRs created in this way can be controlled by setting the flow parameters (Weathers et al. Reference Weathers, Folie, Liu, Childress and Zhang2010; Liu et al. Reference Liu, Ristroph, Weathers, Childress and Zhang2012), but the overall configuration is fixed by the experimental set-up. In contrast to this, numerical simulations can use artificial set-ups to gain deeper insight into vortex dynamics (Pradeep & Hussain Reference Pradeep and Hussain2004; Bergdorf, Koumoutsakos & Leonard Reference Bergdorf, Koumoutsakos and Leonard2007; O’Farrell & Dabiri Reference O’Farrell and Dabiri2012).
Researchers are faced with two main challenges in studying VRPs: (i) time-dependent wake measurement and (ii) control of the wake configuration. Vortex ring pairs, which appear as highly complicated three-dimensional vortical structures (for instance, see Warrick et al. Reference Warrick, Tobalske and Powers2005), are difficult to measure accurately, making it challenging to reveal the unsteady mechanisms in the wake. In addition, it is hard to actively vary vortical configurations by using a single experimental set-up. Since the evolution of wake vortex configurations is related to the direction and magnitude of the wake impulse, control of the wake implies control of the transient forces. However, because of the experimental difficulties outlined above, we do not know how transient forces and the locomotion performance respond to vortical configurations, and whether biological locomotion can be evaluated from the perspective of the dynamics of VRPs.
In this paper, numerical simulation is performed to investigate VRPs based on the piston–cylinder apparatus. We model the NeVR of a VRP by making reasonable assumptions about their vorticity distributions and spatial locations, which are initially prescribed. The flow induced by the NeVR is then superimposed on the flow where a PoVR is just generated from the apparatus (the stroke length is temporarily quite small). The result of superimposition is processed through pre-simulation, and then treated as the initial condition of the simulation. After the simulation begins, the fluid ejected from the apparatus rolls up into a much bigger PoVR. The counter-rotating VRP composed of a PoVR and an NeVR evolves with time under the interaction. The transient forces, calculated from the vortical impulse in the entire wake, are compared with those in the isolated PoVR flow under the same velocity inlet conditions.
The current study aims to quantitatively reveal the effect of VRPs on transient changes in vortical impulse, which is referred to as the force due to VRPs. The force generation is expected to be correlated to certain vortex configurations. It is also intended to study stroke-averaged forces for different amounts of momentum input.
The remainder of the paper is organised as follows. Section 2 introduces the computational set-up, flow initialisation and modelling of NeVRs. Section 3 investigates transient force augmentation at the beginning stage of a stroke, including the effects of circulation, vortex spacing and vortex configurations. Section 4 studies force augmentation over the entire stroke. A normalised force coefficient characterising the stroke-averaged force augmentation is given. Section 5 shows fully unsteady simulation of a VRP, including physical generation of a PoVR. The results are compared with those in § 4. Section 6 gives concluding remarks and future work.
2. Numerical approach
2.1. Computational set-up
Axisymmetric flows are assumed throughout the current study. The governing equations in the numerical approach are the axisymmetric, unsteady, incompressible and laminar Navier–Stokes equations, given by
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\partial u_{x}}{\partial x}+\frac{1}{r}\frac{\partial }{\partial r}\left(ru_{r}\right)=0,\\ \displaystyle \frac{\partial u_{x}}{\partial t}+u_{x}\frac{\partial u_{x}}{\partial x}+u_{r}\frac{\partial u_{x}}{\partial r}=-\frac{1}{{\it\rho}}\frac{\partial p}{\partial x}+{\it\nu}\left[\frac{\partial ^{2}u_{x}}{\partial x^{2}}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{x}}{\partial r}\right)\right],\\ \displaystyle \frac{\partial u_{r}}{\partial t}+u_{x}\frac{\partial u_{r}}{\partial x}+u_{r}\frac{\partial u_{r}}{\partial r}=-\frac{1}{{\it\rho}}\frac{\partial p}{\partial r}+{\it\nu}\left[\frac{\partial ^{2}u_{r}}{\partial x^{2}}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{r}}{\partial r}\right)-\frac{u_{r}}{r^{2}}\right],\end{array}\right\}\end{eqnarray}$$
where
$u_{x}$
and
$u_{r}$
stand for the velocity components in the axial (
$x$
) and radial (
$r$
) directions respectively, with
$t$
being time;
$p$
,
${\it\rho}$
and
${\it\nu}$
are the pressure, density and kinematic viscosity of the fluid respectively.
Figure 1. A sketch of the computational domain (not to scale). The fluid enters the straight tube from the inlet, goes through the tube, exits from the orifice and rolls up into a PoVR. The velocity program at the inlet simulates the piston motion.
The computational domain is sketched in figure 1 (not to scale in order to clearly show different regions of the domain). Since axisymmetric flows are assumed, the domain above the axis is shown. A squared-off lip at the orifice is adopted to avoid singularity, which follows Rosenfeld, Rambod & Gharib (Reference Rosenfeld, Rambod and Gharib1998). In the domain, the straight tube between the inlet and the orifice is
$30D$
in length, where
$D$
is the inner diameter of the tube. The outer boundary is placed at a distance of
$4D$
from the axis. The horizontal distance between the wall and the orifice is also
$4D$
, while the downstream boundary is placed at a distance of
$30D$
from the orifice. The size of the domain is sufficiently large that the computational results are independent of it. The physical scale of the domain in the simulation is chosen to be similar to that in the experiment of Gharib, Rambod & Shariff (Reference Gharib, Rambod and Shariff1998), and the choice of the scale of the domain (normalised by
$D$
) is close to that of previous studies using numerical simulation (Rosenfeld, Katija & Dabiri Reference Rosenfeld, Katija and Dabiri2009; Fu, Qin & Liu Reference Fu, Qin and Liu2014).
The piston used in experiments is modelled as a velocity inlet condition during the simulation; thus, no rigid body simulates the piston in the computation. The velocity profile at the inlet is independent of the radial direction. A no-slip condition is set at the walls. On the outer and downstream boundary, the gauge pressure of the fluid is fixed as zero and flow entrainment is allowed.
The commercial finite volume package ANSYS Fluent 12.0 is employed. Second-order temporal and spatial schemes are used with pressure-implicit with splitting of operators (PISO) pressure–velocity coupling. For each time step, the relative convergence criterion is
$10^{-3}$
for every variable calculated in the solver.
The Reynolds number (
$Re$
) is defined as
$Re=uD/{\it\nu}$
, where
$u$
is the inlet velocity. In the current study, the Reynolds number ranges from 750 to 1500. In this range, various numerical and experimental studies have confirmed that laminar vortex rings are generated (Cater, Soria & Lim Reference Cater, Soria and Lim2004; Rosenfeld et al.
Reference Rosenfeld, Katija and Dabiri2009), which validates our assumption of laminar flows in the simulation. In addition, the current Reynolds numbers are in the same order of magnitude as those in previous studies (
$Re\sim 10^{3}$
). Previous studies (Dabiri & Gharib Reference Dabiri and Gharib2004; Linden & Turner Reference Linden and Turner2004; Dabiri Reference Dabiri2009) have found that the flow behaviour is reasonably independent of
$Re$
, for this range of values of
$Re$
.
The inlet velocity is constant in each stroke throughout this paper, whereas its value varies from stroke to stroke in §§ 4 and 5 in order to reach different effective stroke lengths (
$L$
) over the entire period. In the simulations, a normalised velocity is used, defined as
$u^{\ast }=u/u_{0}$
, where
$u_{0}$
is the reference inlet velocity. It is chosen to satisfy
$u_{0}T/D=4$
, where
$T$
is the period of the entire stroke. In the recent literature on vortex rings (Dabiri Reference Dabiri2009; O’Farrell & Dabiri Reference O’Farrell and Dabiri2014), the formation time (
$t^{\ast }$
) is often used to describe unsteady vortex formation, defined as
$t^{\ast }=ut/D=L(t)/D$
. When the stroke is over, the PoVR has a final stroke ratio
$T^{\ast }$
given by
$T^{\ast }=L/D=uT/D$
.
Based on the computational results showing vorticity contours, the boundary of a vortex is empirically determined as the closed curve inside which the magnitude of the vorticity is not smaller than 3 % of the local maximum (Rosenfeld et al. Reference Rosenfeld, Rambod and Gharib1998, Reference Rosenfeld, Katija and Dabiri2009; Fu et al. Reference Fu, Qin and Liu2014). The area enclosed by the boundary is referred to as the core of a vortex ring. The circulation is calculated by integrating the vorticity inside the defined vortex boundary. Other physical parameters associated with vortices (e.g. impulse) are also evaluated in a similar way in the rest of the paper.
Figure 2. Total circulation in the computational domain downstream of the orifice. The current numerical result is compared with the experimental result (Gharib et al. Reference Gharib, Rambod and Shariff1998) for (a) mesh dependence and (b) time step dependence.
2.2. Flow evolution and validation
The grid and time dependence are checked with the experimental result (Gharib et al.
Reference Gharib, Rambod and Shariff1998) on an isolated PoVR, shown in figure 2. The horizontal axis plots the formation time of the vortex ring; the vertical axis plots the non-dimensional circulation, where
${\it\Gamma}_{to}$
is the time-dependent total circulation integrated over the entire computational domain downstream of the orifice. In the validation, an isolated PoVR is generated with
$ut/D=6$
, where the inlet velocity is constant. After
$ut/D=6$
, the inlet velocity drops to zero instantly, with the vortex ring propagating and the shear layer rolling up into secondary vortices. The total number of grids (
$N$
) ranges from
$0.9\times 10^{4}$
to
$11\times 10^{4}$
in figure 2(a), where the time step (
${\rm\Delta}t^{\ast }$
) is fixed to be 0.01 (
${\rm\Delta}t^{\ast }=u{\rm\Delta}t/D$
, where
${\rm\Delta}t$
is the physical time step). It can be seen from figure 2(a) that good agreement with the experiment is obtained for the stages of vortex ring formation (
$ut/D\leqslant 6$
) and propagation (
$ut/D>6$
). In figure 2(b), the time step ranges from 0.005 to 0.02, while the total number of grids is fixed to be
$2.9\times 10^{4}$
. It can be seen that nice agreement is also obtained regardless of the time step employed. For
$ut/D>6$
, the total circulation in the simulation is still close to the experiment value, without a decaying behaviour, suggesting that the numerical viscosity is not high in the current solver. Considering computational efficiency, a grid size of
$2.9\times 10^{4}$
and a time step of 0.01 are chosen in the rest of the paper, so that the maximum difference relative to the experimental result is restricted to 5 %.
Figure 3 shows the axial velocity profile at the nozzle exit, where
$u_{e}$
is the axial velocity at the nozzle exit. For the total number of grids equal to
$0.9\times 10^{4}$
,
$2.9\times 10^{4}$
and
$11\times 10^{4}$
, there are 10, 20 and 40 grid points across the boundary layer near the wall at
$t^{\ast }=2$
respectively. The results are obtained with different grid sizes and compared with a previous numerical study (Rosenfeld et al.
Reference Rosenfeld, Katija and Dabiri2009) at different formation times (
$t^{\ast }=2$
and 4). Here, we use the same flow parameters as were used there. It can be seen in figure 3 that the agreement of the velocity profile is very good for different grid sizes and it shows a convergent behaviour, which suggests that we can resolve the flow well at the nozzle tip.
Figure 4 shows the rolling-up process of a vortex ring at the circular nozzle, and figures 4(a) and 4(b) are previous experimental and numerical results of Didden (Reference Didden1979) and Nitsche & Krasny (Reference Nitsche and Krasny1994) respectively. Here, we adopt the same values of the velocity program and diameter as they used. The piston starts from a stationary state, with its velocity increasing to a steady state, and suddenly stops. It can be seen that when the flow starts, the fluid goes out of the nozzle exit and rolls up into a small vortex ring. The ring increases in diameter and goes further downstream. The nozzle lip results in a small secondary vortex with opposite vorticity. After the piston stops, the leading vortex ring stops growing but the secondary vortex continues to grow.
Figure 3. Axial velocity profile at the nozzle exit: (a)
$t^{\ast }=2$
; (b)
$t^{\ast }=4$
. The black squares mark the results of Rosenfeld et al. (Reference Rosenfeld, Katija and Dabiri2009).
Figure 4. The rolling-up process of a vortex ring at a circular nozzle for different formation times: (a,e,i)
$=0.29$
; (b,f,j)
$=0.96$
; (c,g,k)
$=1.53$
; (d,h,l)
$=2.10$
. The black arrow indicates the flow direction. (a–d) Visualisation with dye filaments (Didden Reference Didden1979); (e–h) numerical simulation (Nitsche & Krasny Reference Nitsche and Krasny1994); (i–l) current results (white and black regions for positive and negative vorticity contours respectively).
We also compare quantitatively the circulation shedding rate
$\text{d}{\it\Gamma}_{i}/\text{d}t$
with experimental measurements (Didden Reference Didden1979) and numerical simulation (Nitsche & Krasny Reference Nitsche and Krasny1994);
${\it\Gamma}_{i}$
is defined in equation (6) of Didden (Reference Didden1979). The shedding rate is closely associated with the unsteady process of vortex ring formation, because the shear layer shed from the inner edge of the nozzle feeds vorticity into the vortex ring. The result is plotted in figure 5. It can be seen that the current results agree reasonably well with the experimental measurement (Didden Reference Didden1979). For
$ut/D<0.2$
, the current simulation captures the rapid increase of
$\text{d}{\it\Gamma}_{i}/\text{d}t$
. For
$ut/D>1$
, the change of circulation remains almost invariant, and the differences between the current simulation and previous results are minor. From these results, we can conclude that the current simulation is capable of obtaining reliable results for vortex ring formation and evolution, particularly the vorticity close to the nozzle tip.
Figure 5. Comparison of the circulation shedding rate
$\text{d}{\it\Gamma}_{i}/\text{d}t$
between the previous experiment (squares, Didden (Reference Didden1979)), numerical simulation (dash-dotted line, Nitsche & Krasny (Reference Nitsche and Krasny1994)) and the current results (solid line).
2.3. Flow initialisation
In general, the steps of setting up the initial conditions for simulation are given as follows.
(i) A small PoVR (
$L/D=0.2$
) is generated from the apparatus, and its velocity field is calculated based on the numerical set-up presented in § 2.1. The location of the small PoVR is regarded as the reference point so that we can prescribe the location of the NeVR based on the location of the small PoVR. It should be noted that a typical PoVR in biological flows is much larger (
$L/D>2$
, shown by Linden & Turner (Reference Linden and Turner2004)), so the response after the jet starts can still represent the true response without the small PoVR.
(ii) An NeVR is modelled with a vorticity distribution and assigned a spatial location in order to create the configuration of PoVR and NeVR as desired. The details of the modelling will be discussed in § 2.3.1. The velocity field induced by the NeVR is found using the Biot–Savart law.
(iii) The flow induced by the NeVR calculated in (ii) is superimposed on the flow induced by the PoVR calculated in (i).
(iv) The velocity field obtained in (iii) is imported into the numerical solver, and then the computation referred to as pre-simulation is performed. The computation lasts for one time step, with the same numerical set-up as presented in § 2.1. This step ensures that the flow satisfies the governing equations. This process will be discussed in more detail in § 2.3.1.
(v) The velocity field calculated by the pre-simulation process is regarded as the initial condition, and imported into the solver. From this initial condition, the numerical investigation actually commences, allowing the calculation of the transient force from the fluid.
In §§ 3 and 4, the inlet velocity in step (i) for generating the small PoVR and in step (iv) for pre-simulation is consistent with that in the subsequent actual simulation. The details of NeVR modelling and pre-simulation are presented below.
2.3.1. Modelling of NeVRs
We prescribe a vorticity distribution and a spatial location to an NeVR. The vorticity distribution of the NeVR is modelled using the expression (Pozrikidis Reference Pozrikidis2009)
$$\begin{eqnarray}{\it\omega}(x,r)=\left\{\begin{array}{@{}ll@{}}\displaystyle -\frac{{\rm\pi}{\it\Gamma}^{-}}{({\rm\pi}^{2}-4)a^{2}}\left[1+\cos \frac{{\rm\pi}q(x,r)}{a}\right], & 0<q<a,\\ \displaystyle 0, & q>a,\end{array}\right.\end{eqnarray}$$
where
$$\begin{eqnarray}q(x,r)=[(x-x_{c}^{-})^{2}+(r-r_{c}^{-})^{2}]^{1/2}.\end{eqnarray}$$
In the equations above,
${\it\omega}$
and
${\it\Gamma}^{-}$
are the vorticity and magnitude of circulation of the NeVR, with
$a$
being the radius of the vortex ring core. Vorticity is defined as positive when the vortex has counterclockwise rotational motion in the plane shown in figure 1. This definition is consistent with the definition of a PoVR. The spatial location of the NeVR is determined by
$x_{c}^{-}$
and
$r_{c}^{-}$
, which are the coordinates of the vortex core centre (indicating the local extreme value of vorticity) in the (
$x,r$
) plane, and
$r_{c}$
is defined as the radius of the axisymmetric ring. Through the assignment of the coordinates of the NeVR, the configuration of a VRP can thus be manipulated. In the rest of the paper, the NeVR is modelled with (2.2) and (2.3) unless otherwise stated.
We do not intend to discuss how to realise the prescribed configurations of VRPs in laboratory experiments. The realisation is related to active control of flow structures, which is beyond the scope of the current study. It appears to be a standard procedure to examine the vortex interaction numerically by giving a vorticity distribution as an initial condition, although the configuration may not be easily achieved in reality (e.g. Aref & Zawadzki Reference Aref and Zawadzki1991; Pozrikidis Reference Pozrikidis2009).
A physically relevant flow is also ensured by choosing a reasonable vorticity distribution for the NeVR. Direct numerical simulation of the flow generated by a circular tube has shown that the strongly correlated PoVR and NeVR share a similar vorticity distribution (Nitsche & Krasny Reference Nitsche and Krasny1994). Thus, in order to confirm the reasonableness of the current modelling, we compare the vorticity distributions of the PoVR and NeVR. Specifically, three cases are studied: (1) the NeVR modelled by (2.2) and (2.3); (2) the NeVR after pre-simulation; (3) the isolated PoVR. In case (3), the isolated PoVR is generated by constant inlet velocity to reach a given
$L/D$
(2 or 4, see below). For case (1), the NeVR is modelled with the same magnitude of circulation as the PoVR generated in case (3), while its spatial location is relative to the small PoVR indicated in step (i) of § 2.3.
The vorticity distributions are plotted in figure 6. In both figures 6(a) and 6(b),
${\it\sigma}=r-r_{c}$
, where
$r_{c}$
is the radial coordinate of the core centre for either the PoVR or the NeVR;
${\it\omega}_{ex}$
is the local extreme value of vorticity. For case (1),
$(x_{c}^{-}-x_{c}^{+})/D=0.4$
and
$r_{c}^{-}/D=0.6$
, where
$x_{c}^{+}$
stands for the
$x$
coordinate of the core centre of the small PoVR generated in step (i) of § 2.3;
$a/D=$
0.25 in figure 6(a) or 0.3 in figure 6(b), where the values are chosen consistent with the PoVR in case (1). For case (2), the normalised inlet velocity is
$u^{\ast }=0.5$
in figure 6(a) or 1.0 in figure 6(b). These values are also consistent with the formation time of the PoVR in case (1), based on the definition of normalised velocity.
It is indicated in both figures 6(a) and 6(b) that all of the vorticity distributions generally share a similar bell-shaped profile. For
${\it\sigma}/a>0$
, it can be seen that the model captures the vorticity distribution of the PoVR well. For
${\it\sigma}/a<-0.5$
, PoVRs have larger vorticity because of the shear layers. The vorticity distributions after the pre-simulation have smaller magnitude because of the interaction between the PoVR and the NeVR. Pozrikidis (Reference Pozrikidis2009) points out that the precise vorticity distribution of the vortex ring core has little effect on propagation of a vortex ring. Considering the similar bell-shaped vorticity profile, it is hence believed that the modelling can represent reasonably the dynamics of real vortex rings.
2.4. Calculation of the force due to VRPs
The hydrodynamic force (
$\boldsymbol{F}$
) exerting on a solid body inside an infinite flow field can be obtained from the vorticity field of the fluid and the velocity field of the solid body, as given by Wu (Reference Wu1981, Reference Wu2005),
$$\begin{eqnarray}\boldsymbol{F}=-\frac{1}{2}{\it\rho}\frac{\text{d}}{\text{d}t}\int _{V}\boldsymbol{x}\times {\bf\omega}\,\text{d}V+{\it\rho}\frac{\text{d}}{\text{d}t}\int _{V_{s}}\boldsymbol{v}_{s}\,\text{d}V,\end{eqnarray}$$
where
$\boldsymbol{x}$
and
${\bf\omega}$
are displacement and vorticity vectors, and
$\boldsymbol{v}_{s}$
is the velocity of the solid body immersed in the fluid. In the first term on the right-hand side of (2.4),
$V$
is the total volume of fluid (
$V_{f}$
) and the solid body (
$V_{s}$
), and
${\bf\omega}$
takes into the account the rotational motion of the solid body, equal to twice its angular velocity. It should be noticed that (2.4) is not valid for the domain shown in figure 1, because the tube extends to the boundary of the domain so that it does not include all of the vorticity.
On the other hand, (2.4) is valid for the hypothetical domain shown in figure 7, where the domain includes all of the vorticity and the moving solid body (piston). For this case, the first term on the right-hand side of (2.4) can be split into two terms, namely the volume integral for
$V_{1}$
and another one for
$V_{2}$
,
$$\begin{eqnarray}\boldsymbol{F}=-\frac{1}{2}{\it\rho}\frac{\text{d}}{\text{d}t}\left(\int _{V_{1}}\boldsymbol{x}\times {\bf\omega}\,\text{d}V+\int _{V_{2}}\boldsymbol{x}\times {\bf\omega}\,\text{d}V\right)+{\it\rho}\frac{\text{d}}{\text{d}t}\int _{V_{s}}\boldsymbol{v}_{s}\,\text{d}V.\end{eqnarray}$$
The grey box denotes the computational domain shown in figure 1. This includes a VRP downstream of the nozzle exit. It is important to note that, in terms of the vorticity field downstream of the nozzle exit, the numerical results (e.g. circulation, circulation shedding rates) driven by a velocity inlet condition agree well with the experimental results driven by a piston, as validated in § 2.2. Therefore, the values of the
$V_{2}$
term are consistent for the cases shown in figures 1 and 7.
Figure 7. The hypothetical domain to be analysed for force calculation with notation (not to scale). Two blocks enclosed by the dashed lines are denoted by
$V_{1}$
and
$V_{2}$
respectively. The area enclosed by the grey box corresponds to the computational domain shown in figure 1. The empty or filled circles denote vortices with positive or negative vorticity respectively.
In the current study, the time derivative of the
$V_{2}$
term is examined (not the total force acting on the solid body). The motivation is that VRPs are widespread in biological propulsion. We would like to study the contribution solely from a VRP and to remove the contributions from other vortical structures inside the tube. Although the time derivative of the
$V_{2}$
term is not equivalent to the transient force acting on the solid body, this term is of particular significance to propulsion because it is directly related to the wake-capture mechanism (Dickinson et al.
Reference Dickinson, Lehmann and Sane1999), and it is the major contribution to the enhancement of transient force given in (2.5).
The time derivative of the impulse of the VRP in
$V_{2}$
is given by
$$\begin{eqnarray}\frac{\text{d}\boldsymbol{I}_{V_{2}}}{\text{d}t}=-\frac{1}{2}{\it\rho}\frac{\text{d}}{\text{d}t}\int _{V_{2}}\boldsymbol{x}\times {\bf\omega}\,\text{d}V,\end{eqnarray}$$
where
$\boldsymbol{I}_{V_{2}}$
means the impulse of the VRP in
$V_{2}$
. Considering the assumption of axisymmetric flows, the impulse has only an axial component. After some simple vector calculation, we have
$$\begin{eqnarray}\frac{\text{d}I_{V_{2}}}{\text{d}t}=-\frac{{\it\rho}}{2}\frac{\text{d}}{\text{d}t}\int _{V_{2}}r{\it\omega}\,\text{d}V=-\frac{{\it\rho}}{2}\frac{\text{d}}{\text{d}t}\int _{S_{2}}r{\it\omega}\,2{\rm\pi}r\,\text{d}x\,\text{d}r=-{\rm\pi}{\it\rho}\frac{\text{d}}{\text{d}t}\int _{S_{2}}r^{2}{\it\omega}\,\text{d}x\,\text{d}r,\end{eqnarray}$$
where
$I_{V_{2}}$
stands for the axial component. Hereafter, the time derivative of the impulse of a VRP is referred to as the force due to a VRP. This force is calculated according to the rightmost expression of (2.7). The integration of the vorticity field is carried out above the axis in the entire wake downstream of the nozzle exit.
The force due to a VRP is compared with the force generated in the flow with only an isolated PoVR. The flow with the isolated PoVR shares the same velocity inlet condition as the current flow with a VRP. In this way, we ensure that the same impulse is delivered to the wake in both kinds of flows. To compare these two cases, we define an augmentation ratio of force (
$F^{\ast }$
), which is given by
$$\begin{eqnarray}F^{\ast }=\frac{\left(\text{d}I_{V_{2}}/\text{d}t\right)_{VRP}}{\left(\text{d}I_{V_{2}}/\text{d}t\right)_{PoVR}},\end{eqnarray}$$
where, under the same inlet velocity,
$(\text{d}I_{V_{2}}/\text{d}t)_{VRP}$
is the force for the current flow with a VRP and
$(\text{d}I_{V_{2}}/\text{d}t)_{PoVR}$
is the force for the flow with an isolated PoVR. It should be noted that the axial component of the force is considered, as the flow is axisymmetric, and that the ratio is time-dependent since the impulse of a VRP varies with time.
Since all vorticity in
$V_{2}$
is included for calculating the impulse change, the effect of the tip vortex at the nozzle tip just on the boundary of
$V_{2}$
is also taken into account. It should be recalled that this tip vortex, which has negative vorticity and stays at the nozzle lip, is present in isolated vortex ring flows, as shown in figure 4; therefore,
$(\text{d}I_{V_{2}}/\text{d}t)_{PoVR}$
includes the impulse change of the tip vortex. Figure 8 shows the force contribution from all of the vorticity (solid line) in
$V_{2}$
and solely from the tip vortex (dash-dotted line), where
$C_{F}$
is
$\text{d}I_{V_{2}}/\text{d}t$
normalised by
${\it\rho}u^{2}D^{2}$
. We use the same velocity program and nozzle diameter here as in figure 5.
Figure 8. The contribution to the change of hydrodynamic impulse from all of the vorticity in
$V_{2}$
(
$C_{F1}$
, solid line) and from the tip vortex (
$C_{F2}$
, dash-dotted line).
It can be found in figure 8 that the changing rate of total hydrodynamic impulse (
$C_{F1}$
) reaches a peak at
$ut/D\approx 0.2$
, goes down to a minimum at
$ut/D=0.45$
, and recovers to an almost constant value afterwards. On the other hand, the changing rate of the hydrodynamic impulse of the tip vortex (
$C_{F2}$
) has a negative contribution to the changing rate of the total hydrodynamic impulse for
$ut/D<0.8$
. Its negative contribution peaks at
$ut/D=0.45$
. On comparing the contribution from the tip vortex with the contribution from all of the vorticity in
$V_{2}$
, we can find in figure 8 that at the beginning of the simulation (
$L/D=0.2$
, as indicated in § 2.3),
$C_{F2}$
is not very significant (approximately 10 % of
$C_{F1}$
), but at
$L/D=0.45$
,
$C_{F2}$
can contribute to 40 % of
$C_{F1}$
, and afterwards it gradually decays. This trend is consistent with the measurement of the shedding rate of the outer edge which is responsible for the formation of negative vorticity at the nozzle lip (Didden Reference Didden1979). From the previous flow visualisation in figure 4(a,b), before the piston stops, the tip vortex cannot be clearly seen. However, in the current study we include the contribution from the tip vortex so that we take into account the effect of the boundary (nozzle tip) on the growth of a VRP downstream.
Apart from the tip vortex, the shear layer close to the nozzle tip is also taken into account in figure 8 (
$C_{F1}$
). Previous studies have shown that the shear layer shed from the solid body and the leading and trailing edge vortices (counterparts of the leading vortex ring here) are the main sources of unsteady force (e.g. figure 7(a) of Lee et al. (Reference Lee, Hsieh, Chang and Chu2012)). This is consistent with the conclusion in some review articles on propulsion performance for biological swimming and flying (Peng & Dabiri Reference Peng and Dabiri2008; Dabiri Reference Dabiri2009).
3. Force augmentation at the beginning stages of strokes
To simulate the formation of a PoVR, a constant normalised inlet velocity is chosen as
$u^{\ast }=0.75$
for both flow initialisation and actual simulation. A velocity with this value can generate an isolated PoVR of
$L/D=3$
when the entire stroke is over, according to our definition of normalised inlet velocity given in § 2.1. For modelling of an NeVR, the radius of its ring core is chosen as
$a/D=0.27$
in (2.2). This choice is consistent with the value obtained from the isolated PoVR. The augmentation ratio
$F^{\ast }$
is evaluated at
$t/T=0.05$
, and the contours are also extracted at this time instant.
3.1. Effects of circulation and vortex spacing
Thanks to the presence of the NeVR, augmentation may result from the interaction between the vortices. The interaction induces velocity and change of the vortical impulse in the flow field. According to the formulation of the Biot–Savart law, the induced effect is dependent on the circulation of the NeVR (
${\it\Gamma}^{-}$
) and the spacing (
$d$
) between the NeVR and the PoVR. The spacing is defined as the horizontal distance between the centres of vortex cores, i.e.
$d=x_{c}^{-}-x_{c}^{+}$
. Therefore, we examine the effects of the circulation of the NeVR and vortex spacing on transient force augmentation in this subsection, where the radius of the NeVR is set to be
$r_{c}^{-}/D=0.60$
as the initial condition.
For the study on the effect of the circulation of the NeVR, we assign it different values with the vortex spacing set to be
$0.76D$
as the initial condition. Figure 9(a) plots the augmentation ratio as a function of the normalised circulation of the NeVR. The circulation of the NeVR is normalised by the inlet velocity and the inner diameter of the tube. It can be seen from the squares in figure 9(a) that the augmentation ratio is approximately a linear function of the circulation of the NeVR. The linear fitting line shows a slope of 0.51.
It should be recalled that in § 2.3 it has been mentioned that the pre-simulation lasts for one time step, with non-zero inlet velocity after the generation of a small PoVR of
$L/D=0.2$
. Here, we provide validation cases for this approach. The procedure for the validation cases is as follows: (1) introduce an NeVR; (2) run a pre-simulation with multiple time steps and a slip boundary condition; (3) continue the simulation with a no-slip boundary condition and inlet flow to generate a PoVR. The purpose of calculating for multiple time steps is to let the flow relax to a normal state. The slip boundary condition is used to avoid the generation of secondary vortices close to the wall (Lim & Nickels Reference Lim, Nickels and Green1995). The validation test is conducted with a range of time steps (10, 20 and 30) used in the pre-simulation, as shown in figure 9(a). As can be seen in the figure, slight differences are observed among the results obtained by using a range of time steps for the pre-simulation. A converging behaviour is also seen as the number of time steps increases. The differences are probably due to the relaxation of the flow to a natural state. Given the fact that the differences between the results of pre-simulations with different numbers of steps are small, the current results are still useful for describing the general character of the flow. Thus, the procedure described in § 2.3 is acceptable and is followed in the rest of the paper.
Figure 9. (a) The augmentation ratio as a function of the normalised circulation of the NeVR. The circulation is extracted after the pre-simulation. The squares indicate results by using the set-up described in § 2.3. The triangles, circles and stars mark the results for 10, 20 and 30 time steps used in the pre-simulation respectively. The straight lines indicate the results of linear fitting, with slopes of 0.51 (squares), 0.48 (triangles), 0.46 (circles) and 0.45 (stars). (b) The augmentation ratio as a function of the normalised vortex spacing. The spacing is extracted after the pre-simulation, with each axis of the plot on a logarithmic scale. The squares indicate simulation results, while the straight line shows the power-law fitting. The fitting line indicates a power of
$-0.31$
.
For the study on the effect of the spacing between the vortices, we assign it different values with the normalised circulation of the NeVR set to be 1.5 as the initial condition. Figure 9(b) plots the augmentation ratio as a function of the normalised vortex spacing, with each axis on a logarithmic scale. Its power-law fitting shows
$F^{\ast }\sim (d/D)^{-0.31}$
. It shows that the force augmentation ratio (
$F^{\ast }$
), representing the interaction effect, decreases with increasing
$d/D$
. This indicates that the effect of interaction becomes weaker with the initial position of an NeVR going further downstream. It can be inferred from the result that the transient force augmentation decays slowly if the NeVR stays in the near downstream.
It can be concluded from figure 9 that the presence of an NeVR can give rise to remarkable force augmentation relative to the force generated from an isolated PoVR, at least at the beginning stage of a stroke.
3.2. Discussion
The vorticity contours computed in § 3.1 are shown in figure 10(a,b). This simulation is initialised with
${\it\Gamma}^{-}/uD=2.0$
and
$d/D=0.76$
. For comparison, figure 10(b) shows the contours of the isolated PoVR with the same velocity input as the case shown in figure 10(a). The coordinates are normalised with
$D$
. The normalised vorticity (
${\it\omega}^{\ast }$
) is defined as
${\it\omega}^{\ast }={\it\omega}D/u$
. As indicated in figure 10(a,b), the shape of the PoVR remains generally the same despite the presence of an NeVR downstream of the orifice in figure 10(a).
Figure 10. (a) The contour plot of the normalised vorticity for the flow composed of the PoVR and NeVR. (b) The contour plot of the normalised vorticity for the isolated PoVR flow with the same inlet velocity as in (a). (c) The contour plot of
$F_{s}^{\ast }$
corresponding to the flow shown in (a). (d) The contour plot of
$F_{s}^{\ast }$
corresponding to the flow shown in (b). In all figures showing contours hereafter, the orifice is located in the range of
$r/D<0.5$
at
$x/D=0$
. The flow comes from the left, and the fluid rolls up into a small vortex ring in the vicinity of the orifice edge (i.e.
$x/D\approx 0$
and
$r/D\approx 0.5$
). The structures are independent of azimuthal angles, so the cross-sections above the axis are shown hereafter.
We also calculate the rightmost expression of (2.7) in every cell in the domain. Similar approaches to find the sources of force have been taken by Lee et al. (Reference Lee, Hsieh, Chang and Chu2012) and Fu et al. (Reference Fu, Qin and Liu2014). Between two adjacent time frames, the rate of the impulse change in each grid cell is then regarded as the contribution of force from each fluid element. The force from every fluid element is normalised and is given by
$$\begin{eqnarray}F_{s}^{\ast }=-\frac{\text{d}I_{cell}/\text{d}t}{{\it\rho}u^{2}S_{cell}}\approx \frac{{\rm\pi}}{u^{2}}\frac{{\rm\Delta}(r^{2}{\it\omega})}{{\rm\Delta}t},\end{eqnarray}$$
where
$F_{s}^{\ast }$
stands for the source of force from each fluid element;
$S_{cell}$
and
$\text{d}I_{cell}/\text{d}t$
are the area of the grid cell and the force due to the fluid element in the grid cell given by (2.7) respectively. It should be noted that we add a minus sign in order to show positive values for the sources leading to thrust. The difference in
$r^{2}{\it\omega}$
between two time steps is evaluated in every grid cell. The contours of
$F_{s}^{\ast }$
are shown in figure 10(c,d) for the VRP and for the isolated PoVR respectively. It should be noted that, in the regions without vorticity,
$r^{2}{\it\omega}$
is simply zero, and therefore the expression for
$F_{s}^{\ast }$
given by (3.1) represents the force contribution from the impulse of the vortices.
It is indicated in figure 10(c) that the NeVR contributes to overall positive force generation. This is attributed to the interaction effect between the PoVR and the NeVR. With the evolution of the VRP, the PoVR increases in the magnitudes of circulation and impulse, while the NeVR decreases in them. According to (2.7), either the increase in the magnitude of the impulse of the PoVR or the decrease in the magnitude of the impulse of the NeVR contributes to positive force generation. Moreover, the NeVR induces a velocity vector pointing in the positive radial direction in the location of the PoVR. The overall effect is that the vortices increase in their radii and impulses, resulting in force augmentation. It is also reported by Fu et al. (Reference Fu, Qin and Liu2014) that the interacting vortex rings embody another time scale, independent of vortex ring formation, related to thrust enhancement due to multiple vortex rings.
3.3. Effect of vortex configurations
3.3.1. Relative positions between PoVRs and NeVRs
While §§ 3.1 and 3.2 discuss force augmentation for a generally fixed configuration (i.e. an NeVR downstream of a PoVR), the combination of the PoVR and the NeVR can take on various configurations. The configuration can be defined by
${\it\theta}$
, as shown in figure 11, assuming one NeVR. It is shown that
${\it\theta}$
is defined as the angle between the axis and the line connecting the centres of the vortex cores. The distance between vortex core centres is defined as
$l$
, where
$l^{2}=(x_{c}^{-}-x_{c}^{+})^{2}+(r_{c}^{-}-r_{c}^{+})^{2}$
. The vortex spacing defined in § 3.1 satisfies
$d=l\cos {\it\theta}$
.
We calculate the force augmentation ratio for six configurations (different values of
${\it\theta}$
), all of which use
${\it\Gamma}^{-}/uD=2.0$
and
$l/D=0.4$
as the initial conditions. With the values of
$l$
and
${\it\theta}$
, the spatial location of the NeVR (i.e.
$x_{c}^{-}$
and
$r_{c}^{-}$
) is determined. The results for different configurations varying in
${\it\theta}$
are shown in table 1. According to the rightmost expression of (2.7), the thrust is dependent on
$r^{2}$
. With the given vorticity distribution and circulation, the vortical impulse of the NeVR can therefore be represented by
$(r_{c}^{-}/D)^{2}$
, which is shown in the third row in table 1. In the last row, the ratio of
$(F^{\ast })/[(r_{c}^{-}/D)^{2}]$
, namely the force augmentation divided by its impulse, suggests the capability of utilising the impulse from the surrounding fluid to generate force enhancement, because a larger ratio means that the system can gain more force augmentation from the given impulse of the NeVR.
Figure 11. A sketch of a vortex configuration dictating the relative positions of the vortices. The empty circle stands for the PoVR, while the filled circle stands for the NeVR. The configuration varies by changing
${\it\theta}$
, which is defined as the angle between the axis and the line connecting the centres of the vortex cores.
Figure 12. A sketch of the special configuration. The empty circle stands for the PoVR, while the filled region stands for the negative vorticity. The jet comes from the tube exit and generates a PoVR surrounded by vorticity of opposite sense.
Table 1. The results for six configurations defined by
${\it\theta}$
. The first row shows the values of
${\it\theta}$
, defining the configurations. The values of
${\it\theta}$
and
$(r_{c}^{-}/D)$
are extracted after the pre-simulation.
It can be seen in table 1 that in general the augmentation ratio shows an upward trend with increasing impulse of the NeVR, indicated by
$(r_{c}^{-}/D)^{2}$
. In terms of the relative positions between PoVRs and NeVRs, the force augmentation ratio increases when the NeVR moves from the inner side of the PoVR to the outer side. The last row suggests that we can gain more force augmentation for
${\it\theta}<0$
than
${\it\theta}>0$
given the same amount of impulse of an NeVR.
3.3.2. Cross-sectional shape of an NeVR
A flexible shape of the cross-section (not always circular) of an NeVR allows closer spacing between the PoVR and the NeVR, as previously observed (Sahin et al. Reference Sahin, Mohseni and Colin2009). The force augmentation has been found to increase significantly for sufficiently small vortex spacing (Fu et al. Reference Fu, Qin and Liu2014). In the following case we consider a ‘special’ configuration, illustrated in figure 12, where the vortex spacing is very small. The jet comes from the tube exit and generates a PoVR surrounded by vorticity of opposite sense. This special configuration is found in some numerical simulations for jellyfish propulsion (e.g. see Sahin et al. Reference Sahin, Mohseni and Colin2009). For this configuration, the philosophy is that the negative vorticity is kept as close to the PoVR as possible. The negative vorticity is modelled as
$$\begin{eqnarray}{\it\omega}(x,r)=\frac{{\it\omega}_{ex}}{2}\cos \left(\frac{{\rm\pi}{\it\theta}}{{\it\theta}_{2}-{\it\theta}_{1}}\right)\left(1+\cos \frac{2q(x,r)-a_{1}-a_{2}}{a_{2}-a_{1}}{\rm\pi}\right),\end{eqnarray}$$
where
$a_{1}\leqslant q\leqslant a_{2}$
and
${\it\theta}_{1}\leqslant {\it\theta}\leqslant {\it\theta}_{2}$
. The value of
$q$
here is defined in a similar way to (2.3), but
$x_{c}$
and
$r_{c}$
refer to the core centres of the PoVR. The parameters
$a_{1}$
,
$a_{2}$
,
${\it\theta}_{1}$
and
${\it\theta}_{2}$
are chosen as
$a_{1}=0.27D$
,
$a_{2}=0.4D$
,
${\it\theta}_{1}=-5{\rm\pi}/6$
and
${\it\theta}_{2}=5{\rm\pi}/6$
. By changing
${\it\omega}_{ex}$
, we get different circulation of the negative vorticity. We use a cosine function to characterise the vorticity profile from the core centre of the PoVR to the boundary of the negative vorticity in (3.2). With increasing
${\it\theta}$
, the negative vorticity is also assumed to have a cosine distribution. The resemblance between the model and what appears in practice is improved by pre-simulation.
Figure 13. The augmentation ratio as a function of the non-dimensional circulation for the ‘special’ configuration. The circulation is extracted after the pre-simulation. The plot uses a linear scale for the horizontal axis and a logarithmic scale for the vertical axis. The squares indicate simulation results, while the line shows the exponential fitting. The line shows the expression given by
$F^{\ast }\sim \exp (0.5\,{\it\Gamma}^{-}/uD)$
.
The augmentation ratio is plotted in figure 13 as a function of the circulation of the negative vorticity. Surprisingly, the augmentation ratio for this special configuration shows an exponential behaviour as a function of the circulation. The simulation results are fitted with the exponential function, giving the relationship as
$F^{\ast }\sim \exp (0.5\cdot {\it\Gamma}^{-}/uD)$
. For
${\it\Gamma}^{-}/uD=2.0$
, the augmentation ratio is as high as 3.9. The exponential behaviour may be attributed to the sufficiently small spacing between the vortices for this special configuration.
The normalised vorticity contour plot of the special configuration is shown in figure 14(a), where the flow is initialised with
${\it\Gamma}^{-}/uD=2.0$
. It can be seen that the PoVR is surrounded by negative vorticity, whereas there are differences between the vorticity distribution shown and modelled. The contour plot of
$F^{\ast }$
is shown in figure 14(b). It can be found that the interacting vortical structure induces an area of force source located at
$x/D<0.2$
and
$0.5<r/D<0.8$
. As the negative vorticity is kept sufficiently close to the PoVR, the configuration allows stronger interaction (Fu et al.
Reference Fu, Qin and Liu2014).
Figure 14. (a) The contour plot of normalised vorticity for the special configuration. A PoVR is surrounded by negative vorticity. (b) The contour plot of force source,
$F^{\ast }$
, corresponding to the flow shown in (a).
4. Force augmentation over the stroke period
In biological flows, the impulse of an NeVR is related to momentum input. The NeVR has a larger impulse for larger momentum input that generates the PoVR, assuming that the NeVR has the same circulation as the PoVR (conservation of circulation). Quantitatively, for the isolated PoVR, the circulation generated by the piston–cylinder apparatus can be estimated by the slug model (Didden Reference Didden1979; Lim & Nickels Reference Lim, Nickels and Green1995; Dabiri Reference Dabiri2009). The model assumes that, within the boundary layer of thickness
${\it\delta}\ll D$
, the
$r$
component of velocity is much less than the
$x$
component. Moreover, the flow outside the boundary layer is assumed to be spatially uniform and equal to the constant
$u$
. Under these assumptions, the circulation of the PoVR can be approximated as
$$\begin{eqnarray}{\it\Gamma}\approx \int _{0}^{T}\int _{D/2-{\it\delta}}^{D/2}u{\it\omega}\,\text{d}r\,\text{d}t\approx \int _{0}^{T}\int _{D/2-{\it\delta}}^{D/2}u\left(-\frac{\text{d}u}{\text{d}r}\right)\text{d}r\,\text{d}t=\int _{0}^{T}\left(-\frac{u^{2}}{2}\right)_{D/2-{\it\delta}}^{D/2}\,\text{d}t=\frac{1}{2}u^{2}T.\end{eqnarray}$$
Some revisions have been made to this model (Krueger Reference Krueger2005; Gao & Yu Reference Gao and Yu2010), but the slug model still reflects the fundamental relationship between the circulation of the vortex ring and the impulse delivery. The NeVR is assumed to have the same magnitude of circulation as the PoVR. Symmetric oscillatory motion leads to counter-rotating vortices with roughly the same amount of circulation. Thus, this assumption would tell us the general situation of force augmentation due to a VRP. According to the slug model, it is estimated here that the circulation of the NeVR is proportional to the square of the inlet velocity in order to reflect the relationship between the velocity input and the defined NeVR.
In this section, the normalised velocity varies from 0.6 to 1.2 in different cases, and is constant in each individual case. According to the definition introduced in § 2.1, these velocities correspond to final stroke ratios ranging from 2.4 to 4.8. The initial position of the NeVR is chosen as
$d/D=0.5$
and
$r_{c}^{-}/D=0.6$
, and the core size is chosen as
$a/D=0.28$
(see (2.2)). In this way, the NeVRs are constructed to resemble biological flows (Dabiri et al.
Reference Dabiri, Colin, Costello and Gharib2005).
Figure 15. Augmentation ratios as functions of time normalised by the stroke period
$T$
for normalised inlet velocities equal to 0.6, 0.8, 1.0 and 1.2. (a) The ratio over the entire stroke. (b) A magnification of (a) for
$t/T<0.25$
.
The force augmentation ratios are plotted in figure 15, as functions of normalised time for normalised inlet velocities equal to 0.6, 0.8, 1.0 and 1.2. Figure 15(a) indicates that large variations in the ratios occur during the first half of the period, whereas there is no evident force augmentation during the second half. It is also shown that for all of the cases the augmentation ratio increases and then declines, stabilising at 1 eventually. Different trends of the relatively large augmentation ratio are clearly observed in figure 15(b). It can be seen that for larger velocity input, the augmentation ratio peaks at a larger value at an earlier time instant, but drops more rapidly.
We define a force coefficient (
${\it\varepsilon}$
) reflecting the force augmentation over the entire stroke, which is given by
$$\begin{eqnarray}{\it\varepsilon}=\frac{\overline{(\text{d}I_{V_{2}}/\text{d}t)_{VRP}-(\text{d}I_{V_{2}}/\text{d}t)_{PoVR}}}{{\it\rho}u^{2}A},\end{eqnarray}$$
where
$A$
is the area of the nozzle exit and the overbar means a stroke-averaged value. In the equation above, the numerator stands for force augmentation and the denominator stands for the inlet momentum flow. The values of
${\it\varepsilon}$
for the different normalised velocities are shown in figure 16. It is shown that the maximum
${\it\varepsilon}$
corresponds to
$u^{\ast }=1.1$
. This normalised velocity corresponds to a final stroke ratio of
$T^{\ast }=4.4$
. For smaller inlet velocity, the system can only generate quite limited force augmentation. With the delivered momentum increasing, the force augmentation rises significantly, at a larger rate than the delivered momentum. Thus, for this range of inlet velocity, the value of
${\it\varepsilon}$
still increases. For larger inlet velocity, the force coefficient reaches a plateau for
$u^{\ast }>1$
.
Figure 16. The force augmentation coefficient
${\it\varepsilon}$
as a function of the normalised velocity.
Figure 17. The temporal evolution of the normalised vorticity: (a)
$u^{\ast }=0.8$
at
$t/T=0.05$
; (b)
$u^{\ast }=0.8$
at
$t/T=0.2$
; (c)
$u^{\ast }=1.2$
at
$t/T=0.05$
; (d)
$u^{\ast }=1.2$
at
$t/T=0.2$
. It should be noted that in (d) different ranges for the
$x$
and
$r$
axes are used to show the vortices clearly.
The contours of the normalised vorticity are shown in figure 17. Figure 17(a,b) shows the contours for
$u^{\ast }=0.8$
, while figure 17(c,d) shows the contours for
$u^{\ast }=1.2$
. Figure 17(a,c) corresponds to
$t/T=0.05$
, while figure 17(b,d) corresponds to
$t/T=0.2$
. It is shown in figure 17(a,b) that for
$u^{\ast }=0.8$
, the vortices are getting closer during the beginning of the stroke. As the stroke continues, the radii of the vortices increase, arriving at a value of
$r_{c}/D\approx 0.9$
. In comparison, it is indicated in figure 17(c,d) that for
$u^{\ast }=1.2$
, the vortices are getting closer more rapidly after the stroke begins. This results in greater force augmentation at the very beginning of the stroke. Due to the larger momentum delivered to the fluid and the stronger induced effect during their evolution, the radii of the vortices increase more rapidly.
The momentum is not effectively delivered to the VRP or entrained as vortical impulse as the momentum input increases. The circulation of the PoVR is calculated after it pinches off from the trailing jet. The normalised circulation is 0.76, 1.09, 1.16 and 1.06 for
$u^{\ast }=0.6,0.8,1.0$
and 1.2 respectively. For the isolated vortex ring flow, if the injected momentum is completely absorbed into the ring, the normalised circulation should be increased with increased
$u^{\ast }$
(or formation time), as suggested by figure 2. On the other hand, it can be seen for the VRP flows that, for gradually increasing
$u^{\ast }$
, the normalised circulation increases slowly and then saturates. This trend suggests that, for relatively large
$u^{\ast }$
, the momentum input is not as effectively absorbed into the PoVR as for small
$u^{\ast }$
. We note that this trend is well beyond the numerical uncertainty resulting from using the threshold to determine pinch-off. As stated in § 2.1, the threshold of the vortex boundary is chosen as 3 % of the core vorticity (local maximum). We use this threshold to determine whether the PoVR has pinched off from the trailing jet. When we change the threshold, the trend in normalised circulation, however, does not change. Furthermore, one can find that the VRP has moved upstream by
$t/T=0.2$
for
$u^{\ast }=1.2$
in figure 17(d). This flow phenomenon corresponds to the thrust in the opposite direction to the jet. It cancels out some of the thrust gained by fluid ejected downstream. These explain that the augmentation ratios drop below 1 after
$t/T=0.2$
for
$u^{\ast }=1.2$
. However, for
$u^{\ast }=0.8$
, the moderate interaction between the correlated vortices enables the augmentation ratio to maintain relatively high values over the first half of the stroke.
5. Fully unsteady interaction between the vortices
In the previous sections, NeVRs are introduced to the simulated flows as initial conditions. The generation of these NeVRs is not taken into account in the previous simulation. In this section, fully unsteady interaction between the NeVR and the PoVR is directly simulated. The NeVR is generated by forcing a negative velocity (
$-x$
direction) at the inlet, and subsequently a positive velocity (
$+x$
direction) is given so that a PoVR is formed in the presence of the NeVR. The velocity program models a piston moving to the left (generating an NeVR) and then to the right (generating a PoVR). The velocity program is given by
$u^{\ast }=-u_{1}$
(
$t<0$
) and
$u^{\ast }=u_{2}$
(
$t>0$
), where
$u_{1}$
and
$u_{2}$
are both positive constants in each simulation. The values of
$u_{1}$
are determined by equating the circulation of the NeVR to the corresponding circulation of the proposed NeVR in § 4 for each normalised velocity
$u^{\ast }$
at
$t=0^{+}$
. The values of
$u_{2}$
are the same as those in the corresponding cases in § 4 for
$t>0$
. It is expected that we would have different vortex configurations from those proposed in § 4. The purpose here is to compare the force and the structures, as the result of fully unsteady interaction, with those as the result of ‘prescribed’ interaction in § 4.
Figure 18. The contour plots of the normalised vorticity with the normalised inlet velocity of 0.8: (a,b) interaction with the proposed NeVR, copied from figure 17(a,b); (c,d) fully unsteady interaction with the NeVR. Here, (a,c) correspond to
$t/T=0.05$
and (b,d) to
$t/T=0.2$
.
The contours of the normalised vorticity during the fully unsteady interaction are shown in figure 18(c,d). For comparison, the contours with the prescribed NeVRs are copied from figure 17(a,b) and shown in figure 18(a,b). The normalised inlet velocity is 0.8 in all of the panels. It can be seen in figure 18(c,d) that, by forcing a negative velocity at the inlet, the NeVR is generated at the inner side of the PoVR, and that it is closer to the PoVR compared with figure 18(a,b). In addition, the cross-section of the NeVR is stretched alongside the PoVR thanks to the strong induced effect during the fully unsteady interaction between the vortices. Since the NeVR can induce an axial component of the velocity, the vortex pair goes further downstream in figure 18(d) than in figure 18(b). It is worth mentioning that the vortex configuration proposed in § 4 (the NeVR is downstream of the PoVR) could still exist in biological locomotion, as sketched by Dickinson (Reference Dickinson1996).
Figure 19. (a) The complete time history of the force ratio; (b) the ratio during
$t/T<0.25$
. The results shown by the dashed lines are copied from figure 15 for comparison. The lines with the same symbols share the same inlet velocity.
The unsteady force ratio is shown in figure 19. The complete time history of the force ratio is indicated in figure 19(a). The ratio during
$t<0.25$
is shown in figure 19(b), where the results indicated by the dashed lines are copied from figure 15 for comparison. It is shown that in general the force ratio has a similar temporal trend in this figure to the one it has in figure 15. It shows remarkably high values at the beginning stage of the stroke, and then decays to 1 for
$t/T>0.5$
. However, the comparison in figure 19(b) indicates clearly that, in the current case, the force ratio is even higher at the beginning stage. This is not surprising since the NeVR is much closer to the PoVR in this case.
The force coefficient
${\it\varepsilon}$
is also plotted for the current cases in figure 20. The results marked by the squares are copied from figure 16 for comparison. We see that the coefficients are larger in the current cases than in the previous cases of § 4. This is attributed to the smaller spacing between the vortices and the resultant stronger interaction, as discussed in §§ 3.1 and 3.3.2. As we have observed in § 4, the force coefficient
${\it\varepsilon}$
grows rapidly for
$0.6<u^{\ast }<0.8$
, and then reaches a plateau for large
$u^{\ast }$
. The difference is that for fully unsteady interaction the force coefficient increases more rapidly and its peak occurs at
$u^{\ast }=0.8$
. The changing cross-sectional shape of the NeVR during the fully unsteady interaction may be a factor. This could be an interesting research topic for further study.
6. Conclusions
Biological swimming and flying usually involve counter-rotating VRPs through ejection of fluid and boundary layer separation at the edges of solid bodies. In this respect, the model flow generated by the piston–cylinder apparatus is a reasonable representation of the interaction between PoVRs and NeVRs in nature or in laboratory experiments. In this paper, we have demonstrated numerically that the force due to a VRP is considerably larger than the force due to an isolated PoVR. The augmentation is directly related to the vortical configuration. The stroke-averaged force augmentation coefficient grows rapidly but then reaches a plateau as the momentum input increases. These results are consistent with the simulation in which the generation of an NeVR is included.
Figure 20. The force coefficient
${\it\varepsilon}$
as a function of the normalised velocity. The circles indicate the results obtained from fully unsteady interaction between the vortices, and the results marked by the squares are copied from figure 16 for comparison.
Future efforts could be directed towards oscillatory three-dimensional wings with spanwise flows. It would be intriguing to study both the common features and the differences between the fluid ejected from the tube and that going past the wing. Detailed understanding of counter-rotating VRPs is also significant to engineering applications. Bio-mimic machines can be designed to realise favourable vortical structures for high-performance self-propulsion. Moreover, the Lagrangian skeleton of VRPs could be uncovered to reveal the process of fluid transport.
Acknowledgements
We appreciate the financial support from the Natural Science Foundation of China (NSFC-91441205), the State Key Development Program of Basic Research of China (2014CB744802) and the Innovation Practice Program of Shanghai Jiao Tong University (IPP-201110248067). The first author is indebted to Y. Jin and Y. Xiang for their programming of the user-defined functions for flow initialisation, and thanks A. Gregory and D. Chappell for proofreading the manuscript.
