2.1 General considerations

For a uniform through-surface flow, the force on a rigid body with surface $S_{B}$ is defined as

(Batchelor Reference Batchelor1967, p. 233), where $p$ is the pressure and ${\bf\tau}$ is the viscous stress tensor. (The momentum flux from the uniform through-surface component is $-\!\int _{S_{B}}{\it\rho}\boldsymbol{u}(\boldsymbol{u}\boldsymbol{\cdot }\hat{\boldsymbol{n}})\,\text{d}S=\mathbf{0}$ .) For an incident flow $\boldsymbol{U}=U_{\infty }\hat{\boldsymbol{x}}$ , the drag coefficient is defined as

where $A$ is the projected cross-sectional area (or width) of the body and is equal to $2a$ , ${\rm\pi}a^{2}$ for a cylinder/sphere, respectively. The pressure and viscous drag coefficient are defined as

The vorticity, ${\it\omega}$ , is non-dimensionalised as

The pressure and viscous drag coefficients are related to the surface vorticity, $\tilde{{\it\omega}}_{s}$ , at $\tilde{r}=r/a=1$ through

and

for a flat plate (Rayleigh Reference Rayleigh1911), cylinder (Dennis & Shimshoni Reference Dennis and Shimshoni1965) and sphere (Cliffe & Lever Reference Cliffe and Lever1985), respectively. ${\it\theta}$ is taken in an anticlockwise manner from the rear of the cylinder or sphere (see figure 2). The pressure coefficient is defined as

where $p_{\infty }$ is the pressure far from the rigid body. To obtain the pressure distribution on the rigid surface, the tangential momentum equation in polar and spherical coordinates (Batchelor Reference Batchelor1967, p. 602) is integrated to yield

and

for a cylinder and a sphere, respectively.

Figure 2.

(a) Schematic of the computational domain used for the planar flow calculation; for axisymmetric flow, the domain is cut in half through the cylinder. The cylinder (of radius $r=a$ ) is placed in a domain of width $2W$ and length $2L$ , where $W=1000a$ and $L=200a$ . In (b), the coordinate system used in this work is shown. Additionally the two circular regions of inner and outer radius (denoted by $R_{1}$ and $R_{2}$ ), radius $1.2a$ and $5a$ , concentric to the cylinder are shown (figure not to scale).

For steady flows, the drag force can be expressed more generally in terms of the volume flux in the wake and source flow. Hunt & Eames (Reference Hunt and Eames2002) showed that in general, the drag force on a body (with no net momentum flux through the body surface) is

where ${\it\rho}$ is density, $Q_{s}$ is the source flow generated by the body and $Q_{w}$ is the volume flux in the narrow downstream wake (the first term is usually described as the Lamb force – Lamb Reference Lamb1932). When $U_{b}=0$ , the source flow is equal to the volume flux in the wake (i.e. $Q_{w}=Q_{s}$ ) and (2.13) reduces to the classic result attributed to Betz (Reference Betz1925):

For a strong uniform sucking flow ( $-{\it\Lambda}\gg 1$ ), the boundary layer thickness for both cylinder and sphere scale as ${\it\nu}/|U_{b}|$ and the boundary layer is thin. As the flow is irrotational outside the thin attached boundary layer and downstream of the body, the wake volume flux $Q_{w}$ is negligible and the drag force tends to

where $A_{b}$ is the surface area of the rigid body. The integral approach provides details of the drag force but does not give any insight into the flow field. Expression (2.15) was derived by Pankhurst & Thwaites (Reference Pankhurst and Thwaites1953, appendix I) for flows with strong suction and no thrust.

Much less general information can be gleaned for the case of a strong blow flow ( ${\it\Lambda}\gg 1$ ). From the general discussion of forces, the source flux $Q_{s}$ is a combination of the wake and blowing flux, $Q_{s}=Q_{w}+Q_{b}$ , which combined with (2.13) gives

(also obtained by Cliffe & Lever (Reference Cliffe and Lever1985, equation (17)). As we shall see, a blowing flow causes the wake width to widen and the wake volume flux to increase, i.e. $Q_{w}\sim A_{b}U_{b}$ , so that the drag force is not as strongly affected by ${\it\Lambda}$ as compared to the dramatic effect of a sucking flow.

The kinematic effect of the rigid surface generates a source flow in the far field and a dipole contribution for a cylinder or sphere. The irrotational solution, $\boldsymbol{u}=\boldsymbol{{\rm\nabla}}{\it\phi}$ , is well known for a uniform flow past a flat plate (length $2a$ )

for a cylinder

and a sphere,

The slip velocity at the surface of the rigid body in a potential flow is $\partial {\it\phi}/\partial s|_{r=a}=-U_{\infty }$ , $-2U_{\infty }\sin {\it\theta}$ or $-(3/2)U_{\infty }\sin {\it\theta}$ , (where $s=a{\it\theta}$ for a cylinder and sphere). To examine the influence of this component, we decompose the velocity into a primary irrotational component and a rotational component to enforce the no-slip condition on the rigid surfaces, i.e.

The purpose of the discussion in this section is to examine the form and role of the rotational component of the flow, using a boundary layer analysis for short time, and a viscous description for blowing for long time. Much of the discussion will centre around the short blowing viscous ( $T_{{\it\nu}}$ ) and inertial ( $T_{b}$ ) time scales set by the blow velocity and the advective time scale ( $T_{a}$ ) defined as

2.2 Viscously dominated flows ( $Re\rightarrow 0$ and $Re_{b}\neq 0$ )

Klettner & Eames (Reference Klettner and Eames2015) analysed the low Reynolds number flow past a cylinder with uniform blowing or sucking. By assuming a symmetrical vorticity field (about ${\it\theta}={\rm\pi}/2$ ), a leading-order solution to the Oseen form of the vorticity equation was constructed:

where

and

The expression is valid for $Re\rightarrow 0$ and $-\infty <{\it\Lambda}<4/Re$ . The pressure and viscous drag coefficients are

such that the ratio of the viscous drag to the total drag coefficient is

For strong suction,

which agrees with (2.15). The flux of vorticity, defined as

is

Dukowicz (Reference Dukowicz1982) analysed the Stokes flow (low $Re$ ) past a blowing sphere, which, similar to the no-blowing case, admits a closed form expression for the flow field and force without the requirement of including an Oseen approximation. We repeat the main results of Dukowicz (Reference Dukowicz1982) (who includes advection and stretching of azimuthal vorticity caused by the blow flow) and show how they provide information for higher Reynolds number scenarios. Dukowicz’s analysis is strictly valid for $Re\ll 1$ and ${\it\Lambda}<0$ , such that the boundary layer is thin compared to the sphere radius. When ${\it\Lambda}>0$ , we require ${\it\Lambda}Re<1$ or ${\it\Lambda}<1/Re$ , which is similar to the cylinder case. Dukowicz (Reference Dukowicz1982) derived an expression for the drag on a rigid sphere for $Re\ll 1$ and finite values of ${\it\Lambda}Re$ and showed that the vorticity field is given by

(The Reynolds number defined by Dukowicz (Reference Dukowicz1982) is based on the sphere radius, and not the diameter. This is why there is the additional factor of $1/2$ in the expressions stated in this paper.) For strongly sucking flow ( $-{\it\Lambda}\gg 1$ ), the vorticity decays exponentially away from the surface and

which agrees with our later boundary layer analysis; while for blowing flows ( ${\it\Lambda}\gg 1$ ),

The drag coefficient is

For strongly sucking flows,

which agrees with (2.15). The fraction of the drag force that comes from the viscous stresses is

which at ${\it\Lambda}=0$ corresponds to a ratio of $2/3$ (for $Re\ll 1$ ); the ratio $C_{{\it\nu}}/C_{D}$ only depends on $Re{\it\Lambda}$ . For strongly sucking flows, $C_{{\it\nu}}/C_{D}=1$ , while for strongly blowing flows, $C_{{\it\nu}}/C_{D}\sim 8/Re{\it\Lambda}$ . The flux of azimuthal vorticity from the sphere surface is

Equation (2.37) can be written in dimensionless form

In the limit of large sucking flows, the vorticity flux tends to

2.3 Unsteady boundary layer analysis ( $Re\gg 1$ and $Re_{b}\gg 1$ )

At $t=0^{-}$ , the flow is irrotational and described by (2.17)–(2.19). When the flow is set into motion at $t=0$ , it must instantaneously adjust to the no-slip boundary condition applied at the rigid walls. For strongly sucking flows, the thickness of the boundary layer thickness scales as ${\it\delta}\sim {\it\nu}/|U_{b}|$ and as we shall see the flow is ultimately steady over a characteristic time scale $T_{{\it\nu}}$ given by (2.21). The boundary layer grows with time as it diffuses which is then balanced by convection resulting from the through-surface flow. When ${\it\Lambda}\gg 1$ , the boundary layer is blown away from the surface and initially ${\it\delta}\sim U_{b}t$ ( ${\it\delta}\ll a$ ).

To understand the development of the viscous and pressure forces with time, we first examine how the boundary layer develops in time. When the boundary layer thickness is much thinner than $a$ , the momentum equation tangential to the rigid surface for a flat plate ( $\mathscr{D}=0$ ), cylinder ( $\mathscr{D}=1$ ) and sphere ( $\mathscr{D}=2$ ) can be written, to leading order as

Previous work on an impulsively started cylinder (e.g. Bar-Lev & Yang Reference Bar-Lev and Yang1975) and sphere (e.g. Basset Reference Basset1888) focussed around using the Rayleigh solution for a flat plate and projecting this onto a cylinder or sphere. The Rayleigh solution, which involves neglecting ${\it\nu}\mathscr{D}r^{-1}(\partial u_{s}/\partial r)$ , can be used to explain the unsteady viscous component of the drag force but cannot account for the pressure drag component, which depends on the normal gradient of the surface vorticity and would otherwise require a second-order solution (see for example Bar-Lev & Yang Reference Bar-Lev and Yang1975). We can explain the dependence of $C_{P}$ and $C_{{\it\nu}}$ on time and ${\it\Lambda}$ by writing $\tilde{t}=t/T_{a}$ , ${\tilde{n}}=(r-a)/a$ and $F_{1}=u_{s}/U_{\infty }$ . Equation (2.40) can then be rewritten as

The dimensionless tangential velocity is required to match across the no-slip condition the inner region of the outer flow. The solution corresponding to $F_{1}({\tilde{n}},\tilde{t})\rightarrow {\it\lambda}\sin {\it\theta}$ as ${\tilde{n}}\rightarrow \infty$ , $F_{1}(0,\tilde{t})=0$ , $F_{1}({\tilde{n}},0)=1$ is

where ${\it\lambda}=1,2$ and $3/2$ for a flat plate, cylinder and sphere, respectively and $\text{erfc}(z)=2/\sqrt{{\rm\pi}}\int _{z}^{\infty }\exp (-t^{2})\,\text{d}t$ (Abramowitz & Stegun Reference Abramowitz and Stegun1972). The $\sin {\it\theta}$ dependence is dropped for a flat plate. The vorticity at the rigid boundary is related to $F_{1}$ by

The force coefficients can then be calculated from (2.7)–(2.9). These results are summarised in table 1 with the range of validity in table 2 and the full expressions are given in appendix A.

Table 1.

Summary of the short and long time results of the unsteady boundary layer analysis where the slip velocity is the difference between the primary irrotational flow and zero tangential velocity imposed on the rigid surface. The full expressions are given in appendix A. The authors cited in brackets studied the relevant case for ${\it\Lambda}=0$ .

Table 2.

Summary of the validity of the present analysis.

2.4 Sucking flow ( ${\it\Lambda}<0$ )

We explore in more detail the transient and steady force for strong sucking flow, $-{\it\Lambda}\gg 1$ , on these rigid bodies. For sucking flows, the vorticity in the boundary layer is to leading order

The leading-order expressions for vorticity for the flow past a sucking cylinder and sphere are described by (B 14) and (C 11) respectively; these in turn reduce to (2.44) when $-{\it\Lambda}\gg 1$ .

The flux of vorticity through the surface from ${\it\theta}=0$ to ${\rm\pi}$ (or over a unit length of the surface), can be calculated from (2.44) giving

for a flat plate, cylinder and sphere, respectively. The pressure gradient normal to the boundary layer is small and the surface pressure corresponds to the irrotational solution and does not contribute to the total drag force on the body. The pressure coefficient is

for a cylinder/sphere, respectively. The first-order analysis (2.43), while quite an accurate description of the vorticity field, it is not able to recover the $\sin ^{2}{\it\theta}$ dependence seen in (2.46). A second-order analysis results in a pressure coefficient distribution of

where the first and second term on the right-hand side are the viscous and inviscid contributions, respectively. This analysis is presented in appendix B.

In the limit of $-{\it\Lambda}\gg 1$ , the flow outside the boundary layer is irrotational and large gradients in the flow are confined to the boundary layer. This is a special case when the total dissipation in the entire body of fluid is finite and the energy equation can be applied to estimate the drag force. This dissipation technique is well known in the context of Levich’s (Reference Levich1949) estimate of the drag on a clean bubble (see Moore Reference Moore1963, Stone Reference Stone1992). From Landau & Lifshitz (Reference Landau and Lifshitz1966, p. 50), the unsteady form of the energy equation (in the frame of reference with the fluid at infinity at rest) is

where $S_{\infty }$ is a control surface far from the body. The viscous dissipation term ${\it\epsilon}$ (Batchelor Reference Batchelor1967, p. 174) is

and is evaluated over the whole flow domain ( $V$ ) excluding the body ( $V_{B}$ ). To estimate the drag force, we must evaluate (2.48) in the frame of reference with the ambient flow at rest and the body moving.

In the far field the flow is irrotational, so that $p+{\it\rho}u^{2}/2=p_{\infty }-{\it\rho}\partial {\it\phi}/\partial t$ and the viscous stresses are weak. Therefore the contribution from the integral over $S_{\infty }$ reduces to

since the flux from the surface of the body is independent of the shape of $S_{\infty }$ . As the body moves with a constant velocity and the flow is steady in the moving frame, the left-hand side of (2.48) is

Under these conditions, the right-hand side of (2.48) can be recast in terms of

Expressing the viscous stresses as ${\bf\tau}={\it\tau}_{n}\hat{\boldsymbol{n}}+{\it\tau}_{s}\hat{\boldsymbol{s}}$ and using the boundary conditions (2.1), (2.2), we can rewrite (2.48) in a general form:

For a clean gas bubble the tangential shear stress is zero ( ${\it\tau}_{s}=0$ ) and $U_{b}=0$ , which reduces to Levich’s result that $F_{D}={\it\epsilon}/U_{\infty }$ . In our case, there is no tangential slip velocity and the last term is zero. We can estimate the various terms on the right-hand side of (2.53) for the case of a cylinder and sphere, using the boundary layer analysis for the velocity field.

For a cylinder, the flow (in the frame in which the far field is at rest) is approximately

The dissipation term in polar coordinates is (Batchelor Reference Batchelor1967, p. 602)

The dissipation can be decomposed into contributions from the boundary layer region ( ${\it\epsilon}_{bl}$ ), the radial blowing component ( ${\it\epsilon}_{B}$ ), non-blowing irrotational component outside the boundary layer ( ${\it\epsilon}_{IR}$ ) and the nonlinear interactions. The exponential terms decay sufficiently rapidly that the various components are estimated as

which gives a total dissipation of

A scaling analysis allows for an estimation of the radial distance where 99 % of each contribution occurs. For a cylinder the boundary layer, dipole and blowing contributions to the dissipation extend to ${\sim}\!2{\it\delta}$ (from the cylinder surface), ${\sim}\!10^{1/2}a$ and ${\sim}\!10a$ (from the cylinder centre) respectively.

The integral of the normal viscous stress over the surface of the cylinder is

As we have already discussed, the irrotational flow outside the boundary layer determines the pressure field at the cylinder surface and since the velocity normal to the cylinder surface is continuous,

because the tangential velocity is reversed in sign across the boundary layer. Combined together, we find for a cylinder

For flow past a sphere, the flow (in the frame in which the far field is at rest) is

For axisymmetric flows the dissipation term is

The various components are estimated as

In combination, the dissipation term is to leading order

For a sphere, a scaling analysis allows for an estimation of where 99 % of each contribution occurs, similar as that for a cylinder. For a sphere the boundary layer, dipole and blowing contributions to the dissipation extend to ${\it\epsilon}_{bl}\sim 2{\it\delta}$ (from the sphere surface), ${\it\epsilon}_{IR}\sim 3a$ and ${\it\epsilon}_{B}\sim 5a$ (from the sphere centre) respectively.

Now, on the surface of the sphere the pressure distribution is

so the pressure contribution to the first integral on the right-hand side of (2.53), over the surface, is ${\rm\pi}a^{2}{\it\rho}U_{\infty }^{2}U_{b}$ . The viscous stress term is

In total, we have

which agrees with (2.15) for large $Re$ .

2.5 Blowing flow ( ${\it\Lambda}>0$ )

For a large blow velocity, we expect the vorticity to form detached shear layers that bear a similarity to the shear layers that occur in other types of flows, for instance, flows past groups of bodies (Nicolle & Eames Reference Nicolle and Eames2011) or porous sheets (Castro Reference Castro1971). Typically, at some distance downstream, the shear layers become unstable and generate a von Kármán vortex street. During the initial phase, the vortex sheets are well separated and the decay in the maximum vorticity is due to diffusion which is a relatively slow process. This can be seen in the flow field figures 16(e–g) and 18(e–g) for a cylinder and sphere, respectively. Depending on the Reynolds number of the flow, the transition to a faster rate of decay is unlikely to be observed since it requires the vorticity in adjacent vortex lumps to merge and annihilate (see Hunt & Eames Reference Hunt and Eames2002).

The downstream flow signature is intrinsically linked to the forces acting the body. Insight can be gained by analysing the viscous and pressure contributions to the force on the body. In the limit of large ${\it\Lambda}$ , viscous effects are weak near the surface of the cylinder and the pressure is dominated by the vortex force term ( $=\!{\it\rho}U_{b}{\it\omega}_{s}$ ). Although the flow will be unsteady far downstream, a steady analysis can identify some of the critical physics of the wake flow and its connection with the near field. Focussing on the two-dimensional case in the limit of ${\it\Lambda}\gg 1$ , the vortical component of the flow is largely advected in the near field by the irrotational flow. A velocity jump of magnitude $\int _{SL}{\it\omega}\,\text{d}\boldsymbol{s}$ is generated across the shear layer (SL), downstream of the cylinder. Here the integral is taken across the streamline. Since the flow in the wake is much less than $U_{\infty }$ , then in the region $y>0$ ,

which means that for ${\it\Lambda}\gg 1$ the vorticity flux through the cylinder surface is

The first term in the Fourier expansion of the vorticity field ( $\tilde{{\it\omega}}_{s}\cong c_{1}\sin {\it\theta}$ ) is

and this allows the pressure distribution and the drag force to be estimated. From (2.11), the surface pressure distribution is approximately

and the drag coefficient is

We now extend our discussion to a blowing sphere. Again, taking the first term in the Fourier expansion of $\tilde{{\it\omega}}_{s}$ (i.e.  $\tilde{{\it\omega}}=c_{1}\sin {\it\theta}$ ), then the flux relationship is,

The pressure distribution is

and the drag coefficient is estimated to be

The flux of azimuthal vorticity is not conserved in axisymmetric flows because vortex stretching is important. Instead we must look at ${\it\omega}_{{\it\phi}}/y_{SL}$ (where $y_{SL}$ is the distance from the centreline) which is proportional to the circulation of the vortical element (Hunt & Eames Reference Hunt and Eames2002). The flux of ${\it\omega}_{{\it\phi}}/y_{SL}$ is conserved providing there is no vorticity annihilation so that in the shear layer,

where the integral is taken across the shear layer. When $Re\gg 1$ , since the flux of circulation is conserved along the shear layer until vorticity annihilation occurs,

or

where $Y_{W}$ is the ultimate distance of the shear layer from the centreline. As with the cylinder case, the velocity jump across the shear layer corresponds to $U_{\infty }$ ; the ultimate wake width is estimated from the potential flow solution to be $Y_{W}/2a\sim {\it\Lambda}^{1/2}$ . Combined together, we estimate

This provides a closure to enable the leading-order terms to be estimated, i.e.

While these estimates are indicative and order of magnitude, they show that the vorticity flux and drag coefficient decrease with ${\it\Lambda}$ in contrast to the case of a cylinder. This is because the stretching of the azimuthal vorticity by the blowing flow means that a smaller surface vorticity is required to generate the step change of $U_{\infty }$ across the wake shear layer.

3.1 Geometry, solver and validation

Figure 2 shows a schematic of the finite computational domain, which for the two-dimensional problem has width $2W$ and length $2L$ ; for an axisymmetric flow the domain radius was $W$ . We chose $W/a=1000$ and $L/a=200$ . The challenge is that the force calculations are sensitive to the finite size of the computational domain. The sensitivity is related to the fact that the boundary conditions imposed on the computational domain are $U_{\infty }\hat{\boldsymbol{x}}$ on the inlet surface, kinematic and no-slip condition on the rigid cylinder, zero flux on the sidewalls and a constant pressure constraint on the outlet. In addition to making the domain large, the boundary layer is extremely thin for the case of strongly sucking flows and high $Re$ . Typically $10$ points are used to resolve the boundary layer whose thickness scales as $1/|Re_{b}|$ , which for $Re=100$ and ${\it\Lambda}=-3$ , gives elements which have sides of length $0.0006a$ . For a blowing flow, the boundary layer is much thicker and elements of the size of $0.006a$ were found to be sufficient for $Re=100$ and their size could be further increased at lower $Re$ . The computational domain in the vicinity of the cylinder and sphere were refined in two circular regions ( $R_{1}$ and $R_{2}$ ) of radius $1.2a$ and $5a$ . The domain was meshed using Gmsh (Geuzaine & Remacle Reference Geuzaine and Remacle2009). Between $300\,000$ and $1.5\times 10^{6}$ nodes were used in the computations. The other challenge was that, in contrast to the non-blowing case where the smaller elements tend to be in the slowly moving part of the domain, in this case, the faster moving region of the flow coincided with the smallest elements, meaning that the time stepping constraint was severe, tending to be $10^{-6}U_{\infty }/a$ for $Re=100$ and ${\it\Lambda}=3$ .

The continuity and momentum equations were solved using an in-house code called ‘ACEsim’ with the characteristic based split scheme (Zienkiewicz, Taylor & Nithiarasu Reference Zienkiewicz, Taylor and Nithiarasu2005) which has been shown to be effective at studying fundamental fluid mechanics problems (Nicolle & Eames Reference Nicolle and Eames2011; Klettner & Eames Reference Klettner and Eames2012). The defining equations are solved in three steps: involving first the calculation of an intermediate velocity field using the forcing by the viscous forces. The second step is the calculation of the pressure field, from the fact that the Laplacian of pressure is proportional to the divergence of the intermediate velocity. The pressure is applied to force the intermediate velocity field to be solenoidal, and the updated velocity field is calculated in the third step. An axisymmetric formulation was developed, tested and applied in this paper.

3.2 Effect of boundedness on drag force

Flow boundedness has a significant influence on the drag calculations. We first discuss the influence of a domain size on the two-dimensional results. The bounding side walls generate a symmetric flow perturbation to the flow in the vicinity of the cylinder, which has a negligible effect on the drag force; but as a consequence of the finite length of the domain, a streamwise flow perturbation is generated in the vicinity of the cylinder. This perturbation ( ${\rm\Delta}U$ ), can be estimated using the method of images to account for the kinematic effect of the channel walls and the inlet condition. This technique has been used by Eames & Duursma (Reference Eames and Duursma1997) for a different problem. The perturbation is given by

The cylinder ‘sees’ an effective streamwise flow $U_{\infty }+{\rm\Delta}U$ . So for ${\it\Lambda}<0$ , the drag on a cylinder is

The drag coefficient in a bounded flow $C_{D}$ is related to an unbounded flow $C_{D_{\infty }}$ , for a sucking flow, through,

When the domain is wide (e.g.  $L/W\gg 1$ ), the correction to the drag is $a{\rm\pi}^{2}{\it\Lambda}/L$ (approximately 30 % of this comes from the nearest image source term). For example, a sucking flow ${\it\Lambda}=-3$ , where $L/a=200$ , we expect that the numerical calculations give a drag coefficient 15 % lower than the unbounded estimate. This technique can also be used to correct the vorticity distribution on the surface of the cylinder. For a blowing flow, we anticipate that $F_{D_{\infty }}\sim {\it\rho}C_{D_{\infty }}U_{\infty }^{2}$ where $C_{D_{\infty }}$ has a weak dependence on $U_{b}$ . In this case, $C_{D}\sim (1+{\rm\Delta}U/U_{\infty })^{2}C_{D_{\infty }}$ and the effect of flow boundedness leads to the drag coefficient being over predicted by 30 % for ${\it\Lambda}=3$ . It is so strong that it changes the reported trend that $C_{D}$ decreases with ${\it\Lambda}$ for cylinders (see later).

For axisymmetric (or three-dimensional) flows through a pipe of radius $L$ , a similar argument of point sources can be applied. In this case, we anticipate ${\rm\Delta}U\sim a^{2}U_{b}^{2}/L^{2}$ and decays so quickly with $L$ that we expect that its influence on the drag coefficient to be negligible.

Figure 3.

The characteristics of the drag force on a cylinder at $Re=1$ and ${\it\Lambda}=-3$  (▵),  $-2~(+),-1$  (▹), 0 (○), 1 (▹), 2  $(+)$ and 3 (▵). In (a), both the pressure coefficient ( $C_{P}$ ) and drag coefficient ( $C_{D}$ ) are plotted as functions of $t$ , normalised by the advective time scale. In (b), the viscous component of the drag coefficient ( $C_{{\it\nu}}$ ) is plotted as a function of $t$ , normalised by the advective time scale. The red and blue symbols correspond to ${\it\Lambda}>0$ and ${<}0$ , respectively. In (a) the solid and dashed lines are the pressure drag coefficient (A 2) and the total drag coefficient, respectively. The dashed line in (b) is the viscous drag coefficient (A 3).

Figure 4.

The characteristics of the drag coefficient and viscous force on a cylinder at $Re=100$ and ${\it\Lambda}=-3$  (▵),  $-2~(+)$ , $-1$  (▹), 0 (○), 1 (▹), 2  $(+)$ and 3  $(\triangle )$ . In (a,b) the drag coefficient is shown as a function of $t$ , normalised by the advective time scale, for ${\it\Lambda}<0$ and ${>}0$ respectively. The viscous component of the drag force is plotted in (c,d); to highlight the different decay rates (for long and short time), the viscous force is plotted on a linear scale in (c) and logarithmic scale in (d). The dashed lines in (a,b) correspond to the total drag coefficient given by (A 2) and (A 3). The full lines in (c,d) are given by (A 3).

4.1 Cylinder

Figure 3(a,b) shows the variation of the pressure ( $C_{P}$ ), viscous ( $C_{{\it\nu}}$ ) and total drag ( $C_{D}$ ) coefficients with time for $Re=1$ . The blue symbols represent ${\it\Lambda}<0$ , black symbols ${\it\Lambda}=0$ and the red symbols ${\it\Lambda}>0$ . The total drag coefficient shows an initial weak dependence on ${\it\Lambda}$ . This is because the circulation across the thin boundary layer, ${\it\kappa}=\int _{0}^{\infty }{\it\omega}\,\text{d}n$ , depends only on the velocity outside the boundary layer, regardless of whether the flow is sucking or blowing, at least over the period that the boundary layer approximation is valid. The analysis in appendix C shows the variation of $C_{P}$ and $C_{{\it\nu}}$ for short time and confirms the $t^{-1/2}$ dependence observed numerically.

Figure 4 shows the results for $Re=100$ . The scalings for ${\it\Lambda}>0$ and ${<}0$ are sufficiently different that the drag coefficients must be plotted separately; for comparison the result for ${\it\Lambda}=0$ are included. The calculations for ${\it\Lambda}<0$ are run with a finer mesh to capture the thin boundary layer evolution and this results in a smaller time step than for ${\it\Lambda}>0$ . Note that to make the comparison easier, the results are normalised by the advective time scale $T_{a}$ . Figure 4(a,b) confirms that for short time, the drag coefficient decays as $t^{-1/2}$ and is weakly dependent on the blow velocity. For ${\it\Lambda}<0$ , $C_{D}$ has tended to a constant in a time that depends on $T_{{\it\nu}}$ . For ${\it\Lambda}>0$ , the drag force settles to a steady state over a time that depends on both $T_{a}$ and $T_{b}$ . The normalised viscous component of the force is plotted against time in figure 4(c,d) on linear and logarithmic scales. There is good agreement with (A 3) confirming that the drag coefficient scales as $2{\rm\pi}|{\it\Lambda}|$ and that this component decays over a time scale of $T_{{\it\nu}}$ . Figure 4(d) shows that the viscous component has an initial $t^{-1/2}$ dependence which is consistent with the boundary layer analysis. The later decay is exponential for the blowing case, up to when the viscous boundary layer becomes thick and streamwise advection becomes important. The pressure component to the total drag force is more difficult to interpret largely because the blowing time scale is comparable to the advective time scale.

4.2 Sphere

The results for the axisymmetric flow past a sphere with a through-surface flow, mirrors closely the observations for a cylinder. For $Re=1$ , figure 5(a,b) shows that the majority of the force has a viscous component with both the viscous and pressure components decaying as $t^{-1/2}$ . The pressure component in this case is one-third of the total drag force. This is consistent with the unsteady boundary layer analysis which is summarised in table 1 with the full expressions given in appendix A.

Figure 6 shows the variation of the total drag coefficient at $Re=100$ . Over a time scale of ${\sim}\!T_{{\it\nu}}$ the drag quickly decays to a steady state for ${\it\Lambda}<0$ and the initial decrease is $t^{-1/2}$ , in a similar manner to the cylinder. As with the case of the cylinder the initial development of the boundary layer structure largely explains the initial decrease of the viscous component of the force. The pressure component to the total drag force for ${\it\Lambda}<0$ is negligible and forms less than 2 % of the total force. For ${\it\Lambda}>0$ , the pressure component partially collapses for intermediate values of $t/T_{b}$ up to $5T_{b}$ , after which the decay for strong blowing decreasing as $t^{1/6}$ before asymptoting to a constant value.

Figure 5.

The characteristics of the drag, pressure and viscous coefficient on a sphere at $Re=1$ and ${\it\Lambda}=-3$ (▵), $-2~(+)$ , $-1$  (▹), 0 (○), 1 (▹), 2  $(+)$ and 3 (▵). In (a), both the pressure coefficient ( $C_{P}$ ) and drag coefficient ( $C_{D}$ ) are plotted as functions of $t$ , where $t$ is normalised by the advective time scale. In (b), the viscous component of the drag coefficient ( $C_{{\it\nu}}$ ) are plotted as a function of $t$ , normalised by the advective time scale. The red and blue symbols correspond to ${\it\Lambda}>0$ and ${<}0$ , respectively. In (a) the solid and dashed lines are the pressure drag coefficient (A 4) and the total drag coefficient respectively. The dashed line in (b) is the viscous drag coefficient (A 5).

Figure 6.

The characteristics of the drag force on a sphere at $Re=100$ and ${\it\Lambda}=-3$  (▵),  $-2~(+)$ , $-1$  (▹), 0 (○), 1 (▹), 2  $(+)$ and 3 (▵). The drag coefficient is shown as a function of $t$ , normalised by the advective time scale, for ${\it\Lambda}<0$ (a) and ${>}0$ (b). In (a,b) the dashed lines correspond to the decay $t^{-1/2}$ . To highlight the different decay rates (for short and long time), the viscous component of the drag force is plotted on a linear scale in (c) and logarithmic scale in (d). The dashed lines in (a,b) correspond to the total drag coefficient given by (A 4) and (A 5). The full lines in (c,d) are given by (A 5).

Figure 7.

The influence of ${\it\Lambda}$ on the drag coefficient for a cylinder and sphere are shown in (a,c) respectively for $Re=1$ ( $\times$ ), 10 (○), 100 (▹). The fraction of the total drag coefficient which is due to viscous effects as a function of ${\it\Lambda}$ is shown in (b,d) indicating that the force is predominately viscous dominated for sucking flows ( ${\it\Lambda}<0$ ) and pressure dominated for blowing flows ( ${\it\Lambda}>0$ ). The red line in (a) and (c) is the asymptotic expression $2{\rm\pi}|{\it\Lambda}|$ . The dashed curves in (a) and (b) correspond to (2.26) and (2.27) (from Klettner & Eames Reference Klettner and Eames2015). The dashed curves in (c) and (d) correspond to (2.34) and (2.36) (from Dukowicz Reference Dukowicz1982). The black curves in (c) and (d) correspond to the numerical results of Cliffe & Lever (Reference Cliffe and Lever1985). The numerical results of Bagchi (Reference Bagchi2007) are shown in (c) and (d) ( $+$ ).

Figure 8.

The steady drag coefficient for a cylinder at $Re=100$ is plotted as a function of ${\it\Lambda}$ as ▹. The drag coefficient is corrected for flow boundedness and plotted in (a,b) as a filled ▹ for ${\it\Lambda}<0$ and ${\it\Lambda}>0$ respectively and the corrected drag coefficient is $C_{D}=(1+{\rm\Delta}U/U_{\infty })C_{D_{\infty }}$ and $C_{D}=(1+{\rm\Delta}U/U_{\infty })^{2}C_{D_{\infty }}$ respectively.

5.1 Cylinder

The steady drag coefficients for $|{\it\Lambda}|\leqslant 3$ are shown in figure 7(a) for $Re=1,10$ and $100$ ; the fraction of the total force due to viscous stresses is shown in figure 7(b).

The influence of the sucking flow has a dramatic effect on the drag coefficient, as shown in figure 7(a). For $Re=10$ and $100$ , the drag coefficient shows an approximately linear increase with $-{\it\Lambda}$ , in accordance with the boundary layer analysis summarised in table 1. For a blowing flow, the numerical results show that the drag coefficient has a weak dependence on ${\it\Lambda}$ , with a small tendency of $C_{D}$ to increase with ${\it\Lambda}$ . In comparison Mathelin et al. (Reference Mathelin, Bataille and Lallemand2002) found that the drag coefficient increased by nearly 30 % for ${\it\Lambda}=2$ (relative to ${\it\Lambda}=0$ ) for a cylinder at $Re=3900$ . However, the domain width was only $10D$ and would have had a significant effect on the drag coefficient. The drag coefficient, in the present work, corrected for flow boundedness is shown for a cylinder and $Re=100$ in figure 8(a,b). The correction is as much as $15\,\%$ for ${\it\Lambda}=-3$ even though the flow domain is extremely large, tending to raise the drag coefficient for ${\it\Lambda}<0$ . In addition figure 8(b) shows that correcting for flow boundedness an increase of $C_{D}$ is not evident for ${\it\Lambda}>0$ .

In the sucking regime and $Re\geqslant 10$ , the drag force is largely dominated by the viscous component for ${\it\Lambda}<0$ and pressure component for ${\it\Lambda}>0$ , which is consistent with the mathematical models. For $Re=1$ , the variation of the viscous component with ${\it\Lambda}$ is 30 % over the range of ${\it\Lambda}$ considered, with $C_{{\it\nu}}/C_{D}\cong 1/2$ . For the strongly blowing case, (2.72) suggests $C_{D}\cong 1.55$ which is close to $C_{D}\cong 1.3$ (figure 8 b).

5.2 Sphere

The steady drag coefficient for $|{\it\Lambda}|\leqslant 3$ is shown in figure 7(c) for $Re=1,10$ and $100$ . The influence of flow boundedness is negligible for a sphere, which explains the excellent agreement between calculations and predictions for $Re=100$ . The results for $Re=1,10,100$ and ${\it\Lambda}>0$ agree with the curve from Cliffe & Lever (Reference Cliffe and Lever1985) and Bagchi (Reference Bagchi2007) confirming that for $Re<100$ , $C_{D}$ decreases as ${\it\Lambda}$ increases. The theoretical prediction by Dukowicz (Reference Dukowicz1982) is plotted for comparison, which gives a good agreement up to $Re=100$ .

The fraction of the drag force which arises from viscous stresses is shown in figure 7(d). For $Re=1$ , the viscous force is approximately two-thirds of the total force. At higher $Re$ , the drag force is largely dominated by viscous forces for ${\it\Lambda}<0$ and pressure forces for ${\it\Lambda}>0$ , consistent with the present mathematical model. The theoretical prediction (2.36) for $C_{{\it\nu}}/C_{D}$ by Dukowicz (Reference Dukowicz1982) agrees with the numerical results for $Re=1$ . For the strongly blowing case, (2.80) suggests $C_{D}\cong 0.9/{\it\Lambda}^{1/2}$ . For $Re=100$ and ${\it\Lambda}=3$ this gives a prediction of $C_{D}\cong 0.52$ which is close to the numerical simulation of $C_{D}\cong 0.5$ (figure 7 c).

Figure 9.

Vorticity distribution in the boundary layer of a cylinder at $Re=1$ , where ${\it\Lambda}=-3$  (▵),  $-2~(+)$ , $-1$  (▹), 0 (○), 1 (▹), 2  $(+)$ and 3 (▵). In (a), the vorticity through the boundary layer along the line of ${\it\theta}={\rm\pi}/2$ , normalised by its value at $\tilde{r}=1$ . In (b), the vorticity distribution over the surface of the cylinder. The black line in (b) is from Dennis & Shimshoni (Reference Dennis and Shimshoni1965) for ${\it\Lambda}=0$ .

Figure 10.

Vorticity distribution in the boundary layer of a cylinder at $Re=100$ , where ${\it\Lambda}=3$ (▵), 2  $(+)$ , 1 (▹), 0 (○),  $-1$  (▹),  $-2~(+)$ and $-3$  (▵). In (a,c) the vorticity through the boundary layer, along the line ${\it\theta}={\rm\pi}/2$ , is plotted; while in (b,d) the vorticity on the surface of the cylinder is plotted. The blue curve in (a,b) corresponds to (2.44).

Figure 11.

Vorticity distribution in the boundary layer of a sphere at $Re=1$ , where ${\it\Lambda}=-3$ (▵), $-2~(+)$ , $-1$  (▹), 0 (○), 1 (▹), 2  $(+)$ and 3 (▵). In (a), the vorticity through the boundary layer along the line of ${\it\theta}={\rm\pi}/2$ , normalised by its value at $\tilde{r}=1$ . In (b), the vorticity distribution over the surface of the sphere. The full curves correspond to the results of Cliffe & Lever (Reference Cliffe and Lever1985) for ${\it\Lambda}=0,1,3$ . (Note, ${\it\Lambda}=2$ was not calculated by Cliffe & Lever (Reference Cliffe and Lever1985).)

Figure 12.

Vorticity distribution in the boundary layer of a sphere at $Re=100$ , for ${\it\Lambda}=-1$ (▹), $-2~(+)$ , $-3$  (▵) and 3 (▵), 2  $(+)$ , 1 (▹) respectively. In (a,c) the vorticity through the boundary layer, along the line ${\it\theta}={\rm\pi}/2$ , is plotted; while in (b,d) the vorticity on the surface of the sphere is plotted. The blue curve in (a,b) corresponds to (2.44).

Figure 13.

Vorticity flux through the surface of (a,b) cylinder and (c,d) sphere for $Re=1$ ( $\times$ ), 10 (○), 100 (▹). On the left-hand side, the flow is sucking, while on the right-hand side the flow is blowing. In (a,c) the strong sucking prediction given by (2.45) respectively, are plotted with full lines. In (a,b) and (c,d) the dashed lines correspond to low Reynolds number predictions given by (2.30) and (2.38) respectively.

Figure 14.

Pressure distribution on the surface of a cylinder (a,b) and sphere (c,d) for ${\it\Lambda}=1$  (▹), 2  $(+)$ , 3 (▵) and $-3$  (▵),  $-2$   $(+)$ , $-1$  (▹) respectively. In (a,c), $Re=1$ , while for (b,d), $Re=100$ . The full blue curves in (b,d) correspond to the potential flow prediction (2.46). The full black lines in (a–d) correspond to the numerical calculations by Dennis & Shimshoni (Reference Dennis and Shimshoni1965), Berthelsen & Faltinsen (Reference Berthelsen and Faltinsen2008), Le Clair, Hamielec & Pruppacher (Reference Le Clair, Hamielec and Pruppacher1970) and Lee (Reference Lee2000) respectively.

6.1 Vorticity field

6.2 Vorticity flux

6.3 Surface pressure

7.1 Cylinder

Figures 15 and 16 show the vorticity distribution (with blue, green and red representing negative, zero and positive vorticity, respectively) and streamlines in the vicinity of a cylinder for $Re=1$ and $Re=100$ respectively. In all the figures depicting flow fields the flow is from left to right.

For $Re=1$ , a weak asymmetry in the flow is clearly observed for ${\it\Lambda}=0$ , indicating that inertia is still present. We also observe that the radial extent of the vorticity field increases with ${\it\Lambda}$ (for a blowing flow) and decreases for sucking flows (as confirmed in figure 15).

The influence of ${\it\Lambda}$ on the vorticity distribution is more pronounced at $Re=100$ . For ${\it\Lambda}<0$ , the boundary layer is extremely thin and is localised around the cylinder surface. Outside the boundary layer, the flow is irrotational. For ${\it\Lambda}>0$ , the vortex sheets that are formed by the blowing flow separate the positive and negative vorticity which is advected downstream. For $t<50T_{a}$ , a vortex bubble forms downstream of the cylinder which is ultimately unstable and is then shed downstream. The remnants of the vortex bubble (i.e. the large blue and red regions of vorticity between the shear layers) can be seen in figure 16(g) for ${\it\Lambda}=3$ . The vorticity in the shear layers decreases due to local cross-stream diffusion. Further downstream the vortex shear layers become unstable and form a von Kármán vortex sheet, as shown in figure 16(e,f) for ${\it\Lambda}=1,2$ respectively.

7.2 Sphere

Figures 17 and 18 shows the vorticity distribution downstream of a sphere for $Re=1$ and $Re=100$ respectively. The calculations are axisymmetric and reflected along the central axis.

Figure 17.

Vorticity field and streamlines for a blowing sphere at $Re=1$ . The images correspond to ${\it\Lambda}=$ (a) $-3$ , (b) $-2$ , (c) $-1$ , (d) $0$ , (e) $1$ , (f) $2$ and (g) $3$ . The vorticity is normalised so its value extends over the range $-1$ (blue) to 1 (red).

Figure 18.

Vorticity field and streamlines for a blowing sphere at $Re=100$ . The images correspond to ${\it\Lambda}=$ (a) $-3$ , (b) $-2$ , (c) $-1$ , (d) $0$ , (e) $1$ , (f) $2$ and (g) $3$ . The vorticity is normalised so its value extends over the range $-1$ (blue) to 1 (red).

Figure 19.

The variation of the components of (2.60) and (2.67) for different blowing velocities at $Re=100$ for (a) a cylinder and (b) a sphere respectively; $F_{D}U_{\infty }$  (▫);  ${\it\epsilon}$  (○);  $U_{b}\int _{S_{B}}(p+(1/2){\it\rho}u^{2}-p_{\infty })\,\text{d}S$  (▵);  $-U_{b}\int _{S_{B}}{\it\tau}_{n}\,\text{d}S$ ( $\ast$ ). Red and black symbols represent full unsteady numerical calculations and steady state calculations respectively. The lines represent the theoretical predictions (2.57)–(2.60) and (2.64)–(2.67) for a cylinder and a sphere respectively. The variation of ${\it\epsilon}$ with blow velocity for a (c) cylinder and (d) sphere for $Re=100$  (○) and $Re=1000$ (▹). The black line represents the combined contribution to the total dissipation from the boundary layer ( ${\it\epsilon}_{bl}$ ), the blowing flow ( ${\it\epsilon}_{B}$ ) and the irrotational field ( ${\it\epsilon}_{IR}$ ) for a cylinder (2.56) and a sphere (2.63).

For $Re=1$ the vorticity field is similar to that of a cylinder at $Re=1$ . For ${\it\Lambda}=0$ , there is a weak asymmetry in the vorticity field resulting from the inertia of the flow.

For $Re=100$ , two distinct shear layers are formed and as the flow calculations are axisymmetric, this stabilises the presence of the wake bubble that is initially seen behind a bluff body (in this case a cylinder and sphere). It takes a considerable amount of time for this wake bubble to grow and be advected downstream which results in the drag coefficient still evolving even after $t=400T_{a}$ (figure 6 a).