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

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.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),

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}$ ,

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

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

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

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.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.3. Effect of vortex configurations