To calculate the energy costs of swimming or flying, it is crucial to evaluate the drag force originating from skin friction. This topic seems not to have received a definite answer, given the difficulty in measuring accurately the friction drag along objects in movement. The incoming flow along a flat plate in a flapping normal motion has been considered, as limit case of a yawed cylinder in uniform flow, and applying the laminar boundary layer assumption it is demonstrated that the longitudinal drag scales as the square root of the normal velocity component. This lends credit to the assumption that a swimming-like motion may induce a drag increase because of the compression of the boundary layer, which is known as the ‘Bone–Lighthill boundary-layer thinning hypothesis’. The boundary-layer model however cannot predict the genuine three-dimensional flow dynamics and in particular the friction at the leeward side of the plate. A multi-domain, parallel, compact finite-differences Navier–Stokes solution procedure is considered, capable of solving the full problem. The time-dependent flow dynamics is analysed and the general trends predicted by the simplified model are confirmed, with however differences in the magnitude of the friction coefficient. A tentative skin friction formula is proposed for flow states along a plate moving at steady as well as periodic normal velocities.
There has been a considerable amount of studies on the energetics of swimming over the past decades and in particular on drag reduction mechanisms (for a fairly recent review, see ). While many investigations focused on the drag reduction mechanisms employed by aquatic animals, Lighthill and others proposed that drag may actually be enhanced by the swimming motion. The explanation proposed by Lighthill , quoting discussions with Bone, is what is sometimes called the ‘Bone–Lighthill boundary-layer thinning hypothesis’, which states that a plate of section s in an external stream velocity U∥ moving perpendicularly to itself at velocity U⊥ has a frictional boundary-layer thickness (on the side towards which the section is moving), such that the drag force per unit surface is τ≈μU∥/δL.
The drag enhancement formula being associated with simple uniform movements of the body in the fluid, it may apply to flapping-like motions, rather than to fish-like swimming . Free flapping wings or plunging aerofoils have for instance been considered in [4–7], to cite just a few studies. In , a rectangular wing flapping sinusoidally has been analysed and the observed loss of symmetry of the wake induced by the lateral edges has been related to unidirectional flight. Coherent motions as attracting states induced by flapping have also been reproduced numerically . The wake of a pinching foil in still environment has been analysed in , and the experimental as well as computational investigation of plunging aerofoils subjected to uniform flow is reported for instance in .
However, the skin friction along elongated bodies in swimming-like motion has found less attention, owing to the difficulty in measuring this quantity. The hypothesis of drag enhancement, as advanced by Lighthill , conflicts with suggested mechanisms of drag reduction [8,9]. This discrepancy is sometimes attributed to the fact that drag is ill-defined, given the difficulty of separating thrust and drag which balance on average when an animal is swimming at constant mean velocity [1,10,11]. While pressure drag is difficult to define since thrust also arises from pressure forces, there is however no doubt about the definition of skin friction drag. Careful measurements of boundary-layer velocity profiles on swimming fish reported in  confirmed that skin friction drag could be enhanced by factors of up to three to five for dogfish. Skin friction enhancement has also been reported in numerical simulations , with however smaller factors.
One important point of the Bone–Lighthill hypothesis is that the enhanced drag is proportional to . It is remarkable that the same scaling was obtained by Taylor  when he analysed semi-empirically the longitudinal drag on a yawed cylinder in uniform flow. In , the yawed cylinder problem has been readdressed, applying boundary-layer theory and a drag coefficient is derived. The plate with finite span is a limit case of this model problem and the scaling of the boundary-layer thinning hypothesis is retrieved. This skin friction enhancement can be understood as resulting from the acceleration of the fluid particles, and in  a two-dimensional model problem which takes into account this effect has been proposed, by confining the flow between the lower moving plate and a free upper boundary at height s/2. The factor 0.6 in the frictional boundary-layer thickness δL proposed by Lighthill is retrieved in this model and confirmed in  by two-dimensional numerical simulations of the Navier–Stokes system.
A full three-dimensional simulation, in the absence of reliable skin friction measurements along a moving plate, remains necessary to confirm the theoretical drag enhancement prediction. Here, a moving rectangular plate with vanishing thickness, that is without form drag, is immersed in a uniform flow. In most of the theoretical investigations on swimming or flying, the resistive forces are decomposed into pressure drag and viscous drag, as for instance in a recent work on the optimal design for undulatory swimming . This decomposition justifies one to analyse separately the skin friction as one component of the total drag. The numerical solution procedure must be capable of handling the plate's edges, which are singularities for the flow field, and the numerical method has to be sufficiently accurate as to provide reliable skin friction values. This is achieved by using a multi-domain approach together with a high-order compact finite-differences discretization, and full three-dimensional simulations have been undertaken in this work for different uniform plate velocities.
In §2 of this paper, the three-dimensional boundary layer model for the moving plate, which has previously been addressed in , is summarized. The three-dimensional numerical solution procedure is explained in §3 and validated for the fixed flat plate boundary layer. The simulation results for the flow around the moving plate are reported in §4. The predictions for different plate velocities are analysed in §5, addressing the question of a skin friction formula and a periodic plate velocity is considered as well. Some conclusions are drawn in §6.
2. Three-dimensional boundary layer model
A plate with span s in a uniform incoming flow U∥ and moving at normal velocity U⊥ is considered, the configuration being sketched in figure 1. The theoretical prediction of the longitudinal drag provided in  is obtained for a yawed elliptic cylinder in a uniform flow illustrated in figure 2, the plate problem being a limit case for an infinite aspect ratio of the elliptic cross section in the (y,z)-plane. In the following, we briefly summarize the results in . The uniform flow is decomposed onto its tangential and normal components, U∥ and U⊥, respectively, as illustrated in figure 2. The problem is considered to be independent of the tangential direction and the x-component of the potential flow is simply U∥. In the normal direction, the potential flow Qe around the cylinder with elliptic cross section is solved using conformal mapping techniques. To solve the boundary-layer inner problem around the elliptic boundary in the (y,z)-plane, coordinates ξ–η attached to the surface are used (figure 2). The boundary layer equations are written in the (ξ,η,x) coordinates which yields 2.1 2.2 2.3
In , a typical length l is defined such that πl is equal to the circumference of the ellipse (and hence πl=2s when the ellipse degenerates into the plate's cross section). The problem is made dimensionless, considering l in the direction ξ tangential to the boundary of the ellipse and a convenient boundary layer length scale is considered in the normal direction η (see  for general boundary-layer modelling), where the Reynolds number is Re⊥=U⊥l/ν. Accordingly, the reference velocities are U⊥ and in the ξ and η directions, respectively. The scaled equations equivalent to (2.1) and (2.2) are solved using the approximate solution of the momentum equations, details being provided in . Note that the developing boundary-layer profile uξ can only be determined as far as the flow is attached: hence, for each aspect ratio b/a, being the limit case of the plate's section parallel to the z-axis, there is a limiting angle θs, marked in figure 2b, at which the flow separates. This boundary-layer analysis, solving uξ and ux, provides a longitudinal drag coefficient C and the longitudinal drag force per unit length is given by 2.4 It is shown in  that C≈1.8 on the whole range of the elliptic cylinder's aspect ratios. For the forthcoming numerical analysis, it is convenient to use U∥ as reference velocity and the plate's span s as the length scale. Defining the Reynolds number 2.5 and given that l=2s/π, the theoretical prediction for the friction drag per unit length of the plate is 2.6 U*⊥=U⊥/U∥ being the dimensionless normal plate velocity. Note that this formula fails when , in which case the classical friction drag formula for a motionless plate in uniform flow U∥ has to be used instead . Formula (2.6) is therefore relevant only for wall velocities above a lower bound, which is likely to depend on the ratio between the plate's span s and length L.
3. Three-dimensional numerical simulation procedure
In order to assess the reliability of the theoretical predictions outlined in §2, the full three-dimensional problem is solved numerically, for a computational domain containing the plate with vanishing thickness. This numerical problem is particularly challenging, given the singularities associated with the leading and trailing edges as well as the lateral boundaries of the plate. Also, the procedure must be sufficiently accurate as to provide reliable skin friction results along the plate. A multi-domain approach has been used for the solution of the Navier–Stokes system (in the following the dimensionless variables are written without asterisks) 3.1 and 3.2
The partition is designed such that the edges of the plate coincide with contour lines of interfaces between subdomains (sketch in figure 3). The Reynolds number Re=U∥d/ν is formed with the incoming uniform flow velocity U∥ and a typical length scale d of the rectangular plate to be specified later. The main aspects of the solution procedure are summarized hereafter. A semi-implicit second-order backward-Euler time integration is used, the nonlinear terms being evaluated through an Adams–Bashforth scheme. A projection method is considered, that is a fractional-step method by solving at each time step tn=nΔt an intermediate pressure and velocity field followed by a pressure correction to ensure incompressibility, known as the Kim–Moin scheme (see  and for a review on projection methods ). Hence, at each time step a series of Helmholtz-type problems 3.3 for the velocity components, with σ=3 Re/(2Δt), and the pressure (with σ=0) have to be solved. The domain Ω=∪Ωk is partitioned into subdomains Ωk with interfaces Γij=Ωi∩Ωj (see the sketch in figure 3) and the Helmholtz problems in each subdomain are 3.4 where g is either an imposed boundary condition on the exterior of the whole computational domain, or a kinematic condition on the plate in the interior, depending on the specific subdomain considered. High-order compact finite differences schemes are considered for the discretization in the three space variables (x,y,z). The schemes are derived for non-uniform meshes : in particular, as shown in , a clustering of the points near the boundary is appropriate for the eighth-order scheme considered here, to avoid oscillations and which enables a boundary closure scheme of the same order as the interior. In a pre-processing step, the second derivative operators in each direction are diagonalized which gives rise to a fast direct solver of the Helmholtz problems in each subdomain during the time-stepping procedure. Continuity of the solution as well as of its normal derivative is required at the domain interfaces Γij and fields , are introduced such that 3.5 3.6 3.7 3.8 3.9
In this system, the right-hand side of equation (3.7), containing the explicit terms for the time discretization, is time-dependent; and at each time step, the boundary value λ on the interfaces has to be computed to fulfil the continuity of the normal derivatives (3.9). The algebraic formulation of this problem leads to a linear system, the solution of which providing the boundary condition between adjacent domains. This system involves the Schur complement matrix , also called influence matrix, and its internal block structure is determined consistently with the subdomain partition in a pre-processing stage. A parallel MPI algorithm has been designed using the Cluster IBM x3750 of the French computer centre IDRIS, a process being assigned to each subdomain. The Schur complement system is solved iteratively using the Portable, Extensible Toolkit for Scientific Computing (PETSc) computational environment  and more specifically the Krylov subspace package (KSP), using hierarchical GMRES options and Block Jocobi preconditioning . In each subdomain Ωk a 30×30×30 mesh has been used and the algorithm proved to scale almost linearly with the number (up to 120) of domains considered.
(a) Flat plate boundary-layer validation
Before addressing the flow along the moving plate, the steady boundary layer along the plate with finite edges has to be computed which subsequently will be used as the initial condition when the plate is set into motion. The edges of the plate, with vanishing thickness placed at y=0 (see sketch in figure 1), are singularities when the plate is in contact with an incoming uniform flow. This difficulty is overcome by construction using the multi-domain approach, the edges being border lines between adjacent domains and hence the singular values do not appear explicitly throughout the computations. A computational Cartesian domain has been considered, the rectangular plate with length L=36 and span s=6 being located in the y=0 plane with the leading edge at xl=6 and centred at z=0. Uniform flow (1,0,0) (the uniform flow U∥ at inflow being the reference velocity) at x=0 is considered and an advection outflow condition is used at x=60. The wall-normal and spanwise components of the flow velocity, respectively v and w, are supposed to vanish far from the plate at y=±8, whereas a far-field Neumann boundary condition is imposed for the streamwise component u. No-slip conditions for the three components of the velocity field are imposed on the plate. A Reynolds number Re=200 has been considered, that is Res=1200 when based on the plate's span s. The multi-domain partition used contains 120 subdomains, with (ndx,ndy,ndz)=(10,4,3) the number of domains in the three directions, that is the plate ranges over six domains in x and one domain in z. Starting with the uniform flow at inflow, the computations have been advanced in time with a time-step Δt=0.005 and at t=90 a quasi-steady flow field was reached. All variables are now dimensionless and the displacement thickness is a convenient lengthscale for the boundary layer along a flat plate. Figure 4a shows the displacement thickness at different spanwise locations. The value does not vary significantly along the span, besides the region close to the edge. The displacement thickness is seen to grow monotonically as expected by the theory , except in the region close to the trailing edge of the plate (with vanishing thickness) at xt=42, where the flow field has a singular behaviour. Note that the maximum value is δ(x)≈0.6 which yields a maximum Reynolds number based on the displacement thickness of Reδ≈120, that is the boundary layer is stable with respect to infinitesimal perturbations (the critical Reynolds number based on δ being ≈520 ). Also, note that the far-field boundary (with ) is sufficiently far away from the boundary-layer edge, the distance for which the boundary-layer profile recovers 99% of the uniform flow being ≈3δ.
The dimensionless friction drag force per unit surface, the skin friction, is computed as 3.10 τ being the shearing stress on the wall, and cf=0.57/Reδ(x) for the Blasius boundary layer along a spanwise infinite flat plate, when made dimensionless with the displacement thickness . This classical boundary layer formula applies for zero-pressure gradient flow as long as the flow remains attached. More involved asymptotics, such as the triple-deck structure of the flow field , have to be used to describe the behaviour near singular points such as leading and trailing edges. In this investigation, we focus on the flow along the plate and only the classical theory is considered for comparison with the numerical Navier–Stokes solution. Figure 4b shows the computed cf value for the flow state at the centre of the plate, which exhibits as expected a singular behaviour at the leading edge xl=6 and the trailing edge xt=42. Along the plate, the skin friction is close to the theoretical Blasius value depicted as the dashed line. The plate's singularities do not induce significant oscillations of the wall-normal velocity gradient and for this test case of a rectangular flat plate, the simulation procedure is seen to provide reliable skin friction values.
4. Flow over the moving plate
Once the steady flow is established, the plate is set into motion, the dimensionless and constant plate velocity U⊥ being from now on written without asterisk. The plate is initially located in the plane y=0 and its spatially uniform displacement is ϕ(t)=U⊥t. A mapping 4.1 with the computational fixed normal coordinate is considered. In the Navier–Stokes system (3.1), the time derivative has to be transformed accordingly and on the plate the kinematic condition applies, that is 4.2 In this procedure and according to the mapping, the far-field boundary, where the flow becomes uniform, remains at a constant distance from the plate throughout the time integration. For the discretization, 120 subdomains have been considered in the multi-domain procedure with the same 30×30×30 mesh per subdomain as for the boundary-layer computation described in §3. The plate with zero thickness, length L=36 and span s=6 forms a rectangle 6≤x≤42, −3≤z≤3 in the plane inside the overall computational domain Ω=[0,60]×[−8,8]×[−9,9].
The Reynolds number is Re=200, or equivalently Res=1200 when based on the plate's span. The system has been integrated in time (with a time-step Δt=0.005) for different plate velocities U⊥, starting with the flow velocity for the fixed plate as initial condition. The instantaneous flow structure around the plate at t=40 is illustrated in figure 5 for U⊥=0.1,0.2,0.3, the z=0 cut of the streamwise velocity field u in the neighbourhood of the plate (with zero thickness but made visible as a thin black line) being shown. For the smaller velocities U⊥=0.1,0.2, the effect of the motion is only visible near the leading edge and downstream of the trailing edge, the boundary-layer structure being qualitatively similar as for a motionless plate, the streamwise velocity component recovering its uniform value u=1 at a small distance from the plate's boundary. For the higher velocity U⊥=0.3, the flow however exhibits a separation at the leading edge which leads to the formation of a reversed flow region at the lower side, the plate being in an upward motion. The streamwise vorticity field ωx=∂w/∂y−∂v/∂z is depicted in figure 6 where a cut at x=L/3 from the leading edge is shown in the (z,y)-plane. Two opposite counter-rotating vortex structures form at the lateral edges of the plate as a consequence of its upward motion. The intensity of the vorticity increases with U⊥. For U⊥=0.3 some imperfect matching, the vorticity involving the gradients of the velocity field, is visible at lines, corresponding to subdomain boundaries, normal to the plate's edges. This is due to the error tolerance of the iterative procedure used to solve the Schur complement matrix system in this numerical problem.
Starting with the boundary-layer flow along the fixed plate and setting the plate into motion, the flow structure undergoes a transient regime and a crucial question is whether it converges to some quasi-steady state during the time integration. The dimensionless friction force per unit surface 4.3 at the lower and upper face of the plate, that is at y=0− and y=0+, respectively, for U⊥=0.1 at x=L/3 and at different times t=20,30,40 is shown in figure 7. It is seen that the flow at t=40 may be considered as to be in a quasi-steady state for this small plate velocity. Note that the lateral edges at z=±3 are singularities for the flow field and the skin friction is plotted except in the very vicinity of the plate's edges. The skin friction for the motionless plate is also shown as the dotted line, which is of course constant along the plate except in the region adjacent to the edges. The viscous friction enhancement is clearly demonstrated, already at this low plate velocity. The skin friction for a higher velocity U⊥=0.3 is shown in figure 8. Now, while at the upper side towards which the plate is moving the friction value shows a convergence behaviour, at the lower side the flow remains unsteady. Indeed, as shown in figure 5, the flow at U⊥=0.3 exhibits a relatively strong separation at the leading edge which in general is synonymous with an unsteady behaviour. Also, at the lower side, the skin friction exhibits two peaks, symmetric with respect to z=0, which are more pronounced for the higher wall velocity. It is likely that this local increase in friction drag is associated with the presence of the edge vorticity structures at the lower side induced by the upward motion and shown in figure 6.
5. Skin friction formula for the moving plate
Making the longitudinal friction drag (2.6) dimensionless using the span s yields 5.1 the dimensionless plate velocity being written without asterisk and the integration is to be taken along the upper and lower side of the span, omitting the plate's edges which are singular points in the numerical integration formula (a simple trapezoidal rule has been used). Whether a viscous drag coefficient can be defined is intimately related to the existence of a quasi-steady state. However, local features of the flow are likely to be unsteady at higher plate velocities, as shown in the previous section, owing to the strong separation of the flow at the leading edge and at the lateral edges. The highest plate velocity considered here is U⊥=0.4 and the spanwise integrated skin friction Cf has been computed up to t=80. The result is shown in figure 9, for t=40,60,80. While near the leading edge the behaviour is highly unsteady, a quasi-steady evolution for this quantity is seen more downstream. This gives some confidence that the viscous friction for different plate velocities can be compared at some fixed time, after the initial transient behaviour has disappeared. Results for U⊥=0.1,0.2,0.3,0.4 at t=40 are shown in figure 10. As expected, no consistent behaviour of the Cf-values is observed in the region close to the leading edge, but more downstream the curves are seen to be not far from parallel to each other. In figure 11, the quantity 5.2 is shown, starting at x=15, that is discarding one fourth of the plate length near the leading edge. While this quantity varies with x, a clustering of the curves, besides that for the lowest wall velocity U⊥=0.1, at a value around C3D≈1.8 is observed. This value is higher than the theoretical coefficient C3D=1.4 (see §2), which is not surprising, because the friction drag contribution beyond the separation line (the plate's lateral edges) is not taken into account in the theoretical model. Also, when deriving the friction drag formula, the boundary-layer structure in the spanwise direction is considered, assuming streamwise invariance of the flow and leading precisely to the scaling (see §2 and the detailed analysis in ). This scaling is of course modified by the streamwise boundary-layer evolution which leads to the observed streamwise dependence of C3D. Also, for low wall velocities, it is more questionable to focus mainly on the spanwise boundary-layer structure which explains that the result at U⊥=0.1 lies a little apart in figure 11.
(a) Periodic plate velocity
The wall motion in any swimming behaviour is periodic and in  it is shown that the normal body velocity for a large number of fishes and cetaceans typically varies from 0.1U∥ to 0.3U∥ from head to tail. In this model, no explicit spatial undulation of the plate is taken into account, but in order to address a periodic motion the wall velocity with A=0.3 and ω=0.06 has been considered. The maximum wall velocity is 0.3 and the displacement ϕ(t) of the plate varies between ±A/ω=±5, which is a rather large amplitude (compared to the plate's length L=36), at least with regard to typical undulatory swimming amplitudes. It would of course be hazardous to infer from a spatially uniform time-periodic motion of the plate the results one would get for a realistic undulatory motion. However, this model problem is likely to be considered as a kind of extreme case, with respect to normal plate velocity and amplitude of motion. The flow behaviour has been computed over two time periods 2 T, with T≈105, and the spanwise integrated friction value Cf is depicted in figure 12 at two positions (x=L/3,L/2) of the plate. This quantity is seen to inherit the periodicity of the plate's motion and as expected, after a transient initial time interval, the distance between two peaks or equivalently between two valleys of the curves is T/2≈52.
The time-averaged skin friction is shown in figure 13 and compared with the spanwise friction drag for the motionless plate. Integrating these curves in the range 12≤x≤36, that is discarding the portions of the plate near the leading and trailing edges, provides drag values of 0.34 and 0.58 for the motionless plate and moving plate, respectively, that is a drag increase of 70% for the plate with the periodic normal velocity. The dotted line in figure 13 shows the skin friction one would get with formula (5.1) (for C3D=1.8), that is , by considering the mean absolute value of the velocity 〈|U⊥|〉=2A/π=0.191. This Cf value is seen to be surprisingly close to the computed mean friction result, over two-thirds of the plate's length.
In , the theoretical prediction of the so-called ‘Bone–Lighthill boundary-layer thinning hypothesis’ had been strengthened by exploring a boundary-layer model along a plate moving at a normal velocity and considered as the limit case of a yawed cylinder configuration. The three-dimensional numerical simulations of this paper reinforce the theoretical prediction. These simulations remain a challenging problem and are particularly time-consuming and only one plate configuration with a length to span ratio L/s=6 has been considered, using a multi-domain Navier–Stokes solver, at a relatively small Reynolds number Res=1200, based on the incoming uniform velocity U∥ and the span s. The longitudinal drag (per unit length) formula is clearly reinforced, at least for wall-normal velocities U⊥ above some lower bound, by the numerical simulation results, with however a drag coefficient C3D slightly varying along the plate's streamwise direction. The computed coefficient is higher than the theoretical value of 1.4 and may roughly be estimated as 1.7<C3D<2 for the different plate's normal velocities considered. Interestingly, this result is not far from the semi-empirical value ≈2.1 used by Taylor . Although a spatially uniform motion of the plate is oversimplified, it however exemplifies the possibility of skin friction enhancement in swimming motion. In particular, a time-periodic spatially uniform motion with a maximum normal velocity U⊥=0.3U∥ of the plate, which is an upper bound regarding fish swimming , is seen to provide a mean skin friction increase, compared to a motionless plate, by roughly a factor of 1.7. Again, it has to be emphasized that the full three-dimensional numerical simulations are computationally involved and could only be performed for a limited set of parameter values. Spatial undulation of the plate will also have to be considered in the future.
Although based on simplified assumptions, our results lend credit to the conclusion that skin friction is enhanced through swimming motion. However, increases by factors between 4 and 10, as proposed among others in [2,26,27], are unlikely.
This work was granted access to the HPC resources of IDRIS under the allocation i20132a1741 made by GENCI (Grand Equipement National de Calcul Intensif).
One contribution of 15 to a Theme Issue ‘Stability, separation and close body interactions’.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.