Skin friction on a moving wall and its implications for swimming animals

Abstract Estimating the energetic costs of undulatory swimming is important for biologists and engineers alike. To calculate these costs it is crucial to evaluate the drag forces originating from skin friction. This topic has been controversial for decades, some claiming that animals use ingenious mechanisms to reduce the drag and others hypothesizing that undulatory swimming motions induce a drag increase because of the compression of the boundary layers. In this paper, we examine this latter hypothesis, known as the ‘Bone–Lighthill boundary-layer thinning hypothesis’, by analysing the skin friction in different model problems. First, we study analytically the longitudinal drag on a yawed cylinder in a uniform flow by using the approximation of the momentum equations in the laminar boundary layers. This allows us to demonstrate and generalize a result first observed semi-empirically by G.I. Taylor in the 1950s: the longitudinal drag scales as the square root of the normal velocity component. This scaling arises because the fluid particles accelerate as they move around the cylinder. Next we propose an analogue two-dimensional problem where the same scaling law is recovered by artificially accelerating the flow in a channel of finite height. This two-dimensional problem is then simulated numerically to assess the robustness of the analytical results when inhomogeneities and unsteadiness are present. It is shown that spatial or temporal changes in the normal velocity usually tend to increase the skin friction compared with the ideal steady case. Finally, these results are discussed in the context of swimming energetics. We find that the undulatory motions of swimming animals increase their skin friction drag by an amount that closely depends on the geometry and the motion. For the model problem considered in this paper the increase is of the order of 20 %.


Introduction
Over the past 75 years, there has been much interest and controversy about the energetic cost of swimming. These studies were initiated by the famous paper of Gray (1936) who suggested that the large velocities achieved by swimming dolphins could not be explained by the available muscle power unless drag reduction mechanisms U. Ehrenstein and C. Eloy were at play, a result known as the Gray's paradox (see the recent review by Fish (2006) for a fresh perspective). In fact, this early attempt to estimate the swimming energetics proved to be incorrect for two reasons. First, the velocities were observed on animals swimming along the side of a boat, and the dolphins are known to make use of the boat wake energy, a phenomenon called wave-riding (Williams et al. 1992;Weihs 2004). Second, dolphins were only observed during a short period of time and the measured velocities thus corresponded to burst or sprint speeds for which specific muscle fibres are employed (Weis-Fogh & Alexander 1977). These specific fibres have much higher power outputs than regular fibres or the human muscle fibres used by Gray (1936) as a model. Although Gray's paradox was flawed, it has inspired a considerable number of studies on the energetics of swimming and, in particular, on drag reduction mechanisms (Gero 1952;Fish & Hui 1991;Fish 2006). Several mechanisms have been suggested including the optimization of the animal shape (Hertel 1966), the use of compliant skin (Kramer 1960;Carpenter, Davies & Lucey 2000), the presence of ridges, folds or riblets (Aleyev 1977;Walsh 1990;Choi, Moin & Kim 1993;Bechert et al. 1997), the release of secretions (Rosen & Cornford 1971;Hoyt 1975), the heating of the boundary layer (Webb 1975) or its relaminarization induced by the swimming motion (Barrett et al. 1999).
In general, none of these mechanisms have proved to decrease the drag by more than 5 %, and the attempts to replicate them on bio-inspired vehicles have usually failed. The only exceptions seem to be the use of polymer additives to transport oil in pipelines or water in fire hoses (Bhushan 2009), and the swimsuits used by professional swimmers before their ban in 2009, but even in this case, it may just have been an effect of their buoyancy.
While many scientists focused on the drag reduction mechanisms employed by aquatic animals, Lighthill and others proposed that drag may actually be enhanced by the swimming motion. Using potential flow theory, Lighthill (1971) showed that the average thrust force generated by the animal's undulatory motion could be calculated from the kinematics of the caudal fin alone. Then, using kinematics data of Bainbridge (1963) on the dace (Leuciscus leuciscus), he showed that this calculated thrust was approximately four times larger than the drag on a rigid body of similar shape. The explanation proposed by Lighthill (1971), quoting discussions with Bone, is what is sometimes called the 'Bone-Lighthill boundary-layer thinning hypothesis'. The best way to introduce this hypothesis is to quote Lighthill (1971) (the notation has been modified slightly to be coherent with the present paper): 'The possibility that large drag augmentation may result from [lateral movements] is suggested by calculations of boundary-layer thicknesses on highly simplified assumptions. For example, a Blasius boundary layer (with uniform external stream velocity U ) has a 'frictional boundary-layer thickness' δ L (defined so that skin friction is µU /δ L ) equal to 3 νx/U at a distance x from the leading edge. We may compare with this value a frictional boundary-layer thickness δ L for a flat section of depth s moving perpendicular to itself with velocity U ⊥ , which is 0.6 √ νs/U ⊥ on the side towards which the section is moving'.
Following Lighthill (1971) and using the same methodology, other authors have shown, for different species, that the calculated thrust was between two and ten times larger than the drag of a motionless animal (Webb 1975;Alexander 1977;Videler 1981). Additionally, Weihs (1974) showed that fish could make use of this phenomenon by adopting a 'burst-and-coast' swimming gait.
However, the Bone-Lighthill hypothesis have always appeared somewhat mysterious since the calculation of the 'frictional boundary-layer thickness', δ L = 0.6 √ νs/U ⊥ , has never been detailed. One important point in this hypothesis is that the resulting enhanced drag is now proportional to √ U ⊥ , where U ⊥ is the normal velocity. It is remarkable that the same scaling was obtained by Taylor (1952) when he analysed semi-empirically the longitudinal drag on a yawed cylinder in uniform flow. But the link between these two results remains to be established. It is sometimes argued that the reactive model of Lighthill (1971) overestimates the thrust (Hess & Videler 1984;Anderson, McGillis & Grosenbaugh 2001;Shirgaonkar, MacIver & Patankar 2009) and that drag enhancement may not be as dramatic as a two-to 10-fold increase. Nevertheless, Anderson et al. (2001) have carefully measured the boundary-layer velocity profiles on swimming fish and confirmed that skin friction drag could be enhanced because of motion by a factor of 1.5-1.9 for scup and 3-5 for dogfish. Numerical simulations of Borazjani & Sotiropoulos (2008) have also shown increase of skin friction drag by swimming motion but with much smaller factors.
The Bone-Lighthill hypothesis of drag enhancement conflicts with Gray's paradox and the suggested mechanisms of drag reduction. This discrepancy is sometimes attributed to the fact that drag is ill-defined. Schultz & Webb (2002), Fish & Lauder (2006) and Shirgaonkar et al. (2009), among others, argued that when an aquatic animal is swimming at constant mean velocity, there is no general way to separate drag and thrust since they balance on average. In our opinion, this argument can be disputed. It may be true that pressure drag is difficult to define since thrust also arises from pressure forces, but skin friction drag can always be defined, if not measured, and knowing whether this drag is enhanced or reduced during swimming motion is of fundamental importance to evaluate the energetic costs of swimming.
In this paper, to address this question, we use combined theoretical and numerical approaches on simplified models. The theoretical aspects involve laminar boundarylayer theory, as described in the classical book of Schlichting (1979). The boundary-layer equations are based on a simplified locally parallel flow assumption for which scaling laws for the drag coefficient can be derived. The reliability of these theoretical prediction can be assessed by computing quasi-steady states for the full non-parallel Navier-Stokes system for the flow along a plate subject to uniform translational motion. The transient regime when the plate is set into motion is characterized by performing time-dependent Navier-Stokes simulations.
More specifically, the outline of this paper is as follows. In § 2 the model problem of a yawed elliptic cylinder in a uniform flow is addressed using a boundary-layer approach and a drag coefficient expression is derived. The plate with finite width as a limiting case of the cylinder is considered in § 3, and it is shown that this three-dimensional flow is equivalent to a two-dimensional problem with a plate moving in the vertical direction and accelerating outer flow. The scaling of the 'Bone-Lighthill boundary-layer thinning hypothesis' is retrieved and its reliability is assessed in § 4 by considering the full Navier-Stokes system, computing quasisteady states for different parameters associated with the two-dimensional problem as well as time-dependent solutions. The expression of the theoretical drag coefficient is applied in § 5 to a model of undulatory swimming and the equivalent twodimensional Navier-Stokes time integration is performed. The results are summarized and discussed in § 6 and a tentative interpretation of the results in the context of swimming energetics is provided. FIGURE 1. Sketch of the three-dimensional problem: (a) an elliptic cylinder is inclined with angle α in a uniform flow of velocity U ∞ ; (b) in the plane perpendicular to the cylinder axis, the boundary-layer problem is two-dimensional. The boundary layer of thickness δ is developing starting from the stagnation point until it separates at angle θ s . In the boundary layer, we define the local curvilinear coordinate system ξ -η. This method is used to obtain the tangential potential velocity Q e on the cylinder wall.

Three-dimensional problem
Consider the problem illustrated in figure 1 of a yawed elliptic cylinder in a uniform flow. This problem was originally analysed by Taylor (1952) and used to discuss the swimming gaits of elongated animals. To calculate the longitudinal component of the skin friction drag, classical methods of boundary layer theory (Schlichting 1979) are used here. First, the potential flow is obtained by conformal mapping of the flow around a circular cylinder. Then the boundary-layer equations are made dimensionless and this provides the natural scaling of the longitudinal drag. Finally, the full boundary-layer problem is solved by using the momentum equations and proper ansatz for the velocity profiles.

Potential flow
Decomposing the uniform velocity U ∞ onto its tangential and normal components, U and U ⊥ , as illustrated in figure 1(a), the potential flow problem can be solved as follows. In the tangential direction, the outer potential flow is unchanged and is equal to U . In the normal direction, the potential flow around the cylinder with elliptic cross-section is obtained by conformal mapping techniques as illustrated in figure 2. This calculation can be found in any good textbook (e.g. Batchelor 1967), and will be briefly recalled here.
In the Z-plane, the complex potential around the circular cylinder of radius c is given by (2.1) and the conformal mapping from the Z-plane into the ζ -plane is given by the transformation where c = (a+b)/2 and λ = √ a 2 − b 2 /2, with a and b the two semi-axes of the ellipse in the directions parallel and normal to the flow, respectively. The complex flow velocity in the ζ -plane is then obtained by differentiating the complex potential with respect to ζ (where ζ = y + iz) Finally, the potential flow velocity at the surface of the elliptic cylinder, Q e , is obtained by assuming that Z = c e iθ and taking the norm of (2.3), which yields where the symmetries with respect to the Oy and Oz axes appear as expected.
To solve the boundary-layer inner problem, we will use the coordinates ξ -η attached to the surface (see figure 1b). We will then need the stretching function dξ/dθ given by the norm of dζ /dθ when Z = c e iθ dξ dθ = dζ dθ = a 1 − (1 − b 2 /a 2 )cos 2 θ. (2.5) For a given elliptic cross-section, a typical length is defined as the diameter of the 'equivalent' circular cylinder, such that this equivalent cylinder has the same surface as the elliptic cylinder where E is the complete elliptic integral of the second kind.

Boundary-layer equations
Using the fact that the problem is independent of x, the boundary-layer equations can be written in the (ξ, η, x) coordinates as The problem can be made dimensionless by using U ⊥ and as characteristic velocity and length scales. Rescaling the variables as follows where the Reynolds number is the dimensionless boundary-layer equations become (2.10c) and the problem now only depends on the aspect ratio of the elliptical cross-section b/a through q e (ξ * ). Note that u * ξ and u * η do not depend on u * x , such that the twodimensional problem for (u * ξ , u * η ), given by (2.10a) and (2.10c), can be first solved in the ξ -η plane and then its solution can be used to solve the linear problem (2.10b) for u * x . From these dimensionless equations, it appears that the longitudinal drag force per unit length of cylinder is given by where C 3D is a longitudinal drag coefficient that depends only on the aspect ratio of the ellipse b/a and is given in dimensionless variables by This longitudinal drag can also be rewritten using the incident angle α (see figure 1a) as where the Reynolds number Re is now (2.14) Equation (2.13) makes apparent the similarity with the semi-empirical formula of Taylor (1952), numbered (2.22) in his paper. The boundary-layer equations (2.10a)-(2.10c) can be solved using the approximate method of the momentum equations (Sears 1948;Wild 1949;Schlichting 1979), as it is detailed in appendix A. The results of this calculation are illustrated in figure 3 where the drag coefficient C 3D is plotted as a function of the aspect ratio b/a. The coefficient C 3D takes its maximal value, C 3D = 1.792, for the aspect ratio b/a = 0.48, but surprisingly, varies very little, and is approximately equal to C 3D ≈ 1.8 on the whole range of aspect ratios. This value is 30 % smaller than the empirical value used by Taylor (1952), C 3D = 2.7, which he found by fitting the data of Relf & Powell (1917) on a yawed circular cylinder at Re ≈ 8000. The reasons for this discrepancy are not fully understood. It could be due to the outer velocity profile on the cylinder that is not the solution of the potential flow problem or contributions of the surface situated beyond the separation line which are neglected in the present analysis. It could also be errors in the measurements, effects of the turbulence level in the wind tunnel, or roughness of the cylinder surface. In addition boundary layers along a yawed cylinder are known to be unstable (Poll 1985). In a flow configuration above a critical Reynolds number, the instability is evidenced by an undulating separation line. This instability is however unlikely to affect significantly the skin friction along the axial direction, as long as the mean position of the separation line is not affected. As already pointed out by Taylor (1952), the present three-dimensional problem of a yawed cylinder in uniform flow is analogous to the problem of heat transfer of a cylinder in cross-flow. This analogy is explored in appendix B for the case of an elliptic cross-section.
The main result of this section is that the longitudinal drag per unit length is given by (2.11), with a longitudinal drag coefficient C 3D ≈ 1.8. Therefore, the skin friction drag will increase as the square root of the normal velocity component, a result of crucial importance to evaluate the energetic costs of undulatory swimming.

Two-dimensional problem
Consider now the three-dimensional problem illustrated in figure 4, which is simply the limit of the yawed cylinder analysed above in § 2 when the aspect ratio b/a tends to infinity, with a change of framework such that now the plate is moving at the normal velocity U ⊥ in a uniform parallel flow U . The potential flow U e at the surface of this plate is given by substituting a = 0 and b = H into (2.3), which yields . Three-dimensional problem of a plate of span s = 2H in a uniform flow U moving at velocity U ⊥ in the direction normal to its surface. Variation of the fluid-particle velocity along these streamlines as a function of the curvilinear coordinate X.
To illustrate this outer flow, its streamlines, as seen from above the plate, are drawn in figure 5(a) in the particular case U ⊥ = U . It is important to note that this flow is associated with an acceleration along the streamlines as fluid particles get around the plate (figure 5b). This acceleration leads to a favourable pressure gradient, a thinner boundary layer thickness and thus an enhancement of skin friction. An alternative way to understand this skin friction enhancement is to remark that, for this three-dimensional problem, the mid-plane Oxy is a plane of symmetry, in which the z-component of the velocity field vanishes. In this mid-plane, applying the conservation of mass in y = 0 and z = 0, one finds that the y-component of the outer potential flow scales as V e = U ⊥ (1 − y/H) near the plate (i.e. for y H). Hence, the fluid particles experience an acceleration towards the plate as they move around the plate, which then yields a compression of the boundary layer and an increase of skin friction drag. Naturally, these two interpretations of the skin friction enhancement either by invoking the acceleration of fluid particles along their streamlines, or by invoking the compression in the vertical direction cannot be dissociated as they are linked by the conservation of mass. These accelerations in the vertical and horizontal directions, which are instrumental in enhancing the skin friction, are inherently linked to the three-dimensional nature of the problem. Yet, an analogue two-dimensional problem can be formulated where these accelerations are recovered by confining the flow between the upward moving plate and an upper free boundary at height H ( figure 6). In this two-dimensional problem, we are looking for the skin friction on the bottom wall of a channel of height H when this wall is moving upward with a prescribed velocity. We will consider the case where this upward velocity can be written with ψ a dimensionless function of order one. The outer flow field, U e = U + ∇φ, can be found by solving a Laplace problem for the potential φ. This is classically achieved by use of the conformal mapping X + iY = e π(x+iy)/H which transforms the channel geometry of height H into the upper half-plane (cf. Batchelor (1967, §6.5)). Green's representation theorem can then be used in the transformed domain to provide the value of φ on the bottom wall for any graph ψ(x) as where the logarithm function is the Green function of the Laplace equation.
In particular, when ψ = Θ(x), with Θ the Heaviside step function, the streamwise component of the outer flow is found by taking the x-derivative of (3.3) and is equal to which yields the following approximation The analogy with the three-dimensional problem of figure 4 appears through the same vertical potential velocity V e as in the mid-plane of the three-dimensional problem, and through a similar acceleration along the streamwise direction. The analogue two-dimensional problem that we have proposed here can be solved following a method similar to that of § 2.2, where the momentum equations can be used to calculate the drag force (cf. Schlichting 1979, pp. 206-214). When the acceleration dU e /dx varies slowly compared with the typical y-scale, the drag force per unit surface on the bottom plate reduces to U. Ehrenstein and C. Eloy which gives, taking dU e /dx = U ⊥ /H, a drag force per unit x (assuming a span s = 2H as in figure 4) where Re ⊥ is given by (2.9), with = 4H/π here. Equation (3.7) is similar to (2.11), and therefore an analogous drag coefficient can be defined for the two-dimensional problem studied here, which depends on x in general In this approximation, this drag coefficient is C 2D = 2.12, close to the value C 3D ≈ 1.8 found for the three-dimensional case of § 2. Using U and δ 1 (the displacement thickness in x = 0) as the characteristic velocity and length to make the problem dimensionless, such that, in particular, the dimensionless shearing stress along the plate can also be calculated Note that the expression (3.6) can be rewritten as assuming again that the span is s = 2H. Here δ L is the 'frictional boundary-layer thickness' used by Lighthill (1971) and its expression is exactly the same as the one he used (see the quotation from Lighthill's paper in § 1). It is thus likely that Lighthill (1971) used the same two-dimensional problem to find the scaling of δ L in the case of a normal component of the velocity.

Two-dimensional Navier-Stokes solution
The theoretical predictions and the corresponding scaling for the skin friction (3.6) are based on the boundary-layer approximation, under the assumption of a locally parallel flow in the streamwise x-direction. To assess the reliability of these hypotheses, the nonlinear Navier-Stokes system is now considered in the two-dimensional domain of height H * and plate velocity U * ⊥ , given by 0 This domain is identical to the two-dimensional flow geometry considered in the previous section, for t = 0 (cf. figure 6). The domain is necessarily of finite length L and incoming boundary-layer flow is assumed at x = 0. In the following only the non-dimensional key parameters (3.9) are written with asterisks, which are omitted for all other variables for simplicity, the reference velocity still being U and the reference length scale being the displacement thickness δ 1 of the incoming flow at x = 0. The numerical procedure uses the mappinḡ which transforms the moving domain into the fixed geometry 0 x L, 0 ȳ H * . The gradient and the time derivative in the transformed system are Performing the mapping, the Navier-Stokes equations become with u = (u, v), the velocities u and v being the streamwise and wall-normal components, respectively, and Re the Reynolds number as defined by (3.9). The dimensionless Navier-Stokes equations (4.4a)-(4.4c) will be solved numerically in the following for different values of the channel height H * , the normal velocity U * ⊥ and the Reynolds number Re. These equations are supplemented with the following boundary conditions. At the lower boundary the fluid velocity field has to match the wall velocity, such that In the analogy between the three-and two-dimensional problems discussed in § 3, the conditions at the upper free slip boundary (cf. figure 6) are At inflow x = 0, an incoming Blasius boundary-layer profile (with displacement thickness equal to one) is prescribed whereas at outflow x = L the Neumann boundary condition ∇ x u = 0 is used to avoid upstream effects due to the finite length of the domain.

4.1.
Quasi-steady states The first numerical simulations address quasi-steady states. These equilibrium states of the system are sought by setting ∂u/∂t = 0 in (4.4a)-(4.4c). The solution procedure is similar to that described by Ehrenstein & Gallaire (2008) for a boundary-layer flow along a bump geometry, the system being discretized using Chebyshev-collocation and a grid-stretching is performed in y to take into account the boundary-layer structure near the wall. To avoid clustering of the points at inflow and outflow, an algebraic one-parameter coordinate transformation (cf. Peyret 2002, p. 306), is used to achieve a quasi-uniform distribution in x. In a manner similar to Ehrenstein & Gallaire (2008), the coupled nonlinear system (4.4) is then solved by using a quasi-Newton method known as the Broyden rank-one update procedure (Stoer & Bulirsch 1992). In this procedure, only one complete QR-decomposition of the nonlinear system's Jacobian matrix has to be performed when starting an iteration cycle: the operators for the subsequent linear systems are updated, using an affine approximation of the nonlinear function, so that the corresponding QR-decompositions can be computed with little extra cost (Stoer & Bulirsch 1992). It can be shown that this procedure converges superlinearly for initial guesses close to an equilibrium state.
A computational domain of length L ≈ 300 has been considered. In the previous section, the outer flow field (3.5) has been obtained assuming the ideal case in which the graph Ψ (x) is the Heaviside step function. However, for the numerical solution in the finite domain, Ψ (x) has to evolve continuously from zero to one within small buffer regions near the inflow and outflow. The graph that has been used is shown in figure 7(a). Equilibrium flow states are obtained starting with U * ⊥ = 0 and considering as initial guess the streamwise-independent Blasius profile: performing the quasi-Newton iterations, the system converges to the flat plate Navier-Stokes solution. The wall velocity U * ⊥ is then progressively increased to compute the successive equilibrium states.
For a given velocity U * ⊥ , the moving boundary ϕ(x, t) = U * ⊥ tΨ (x) can be at different heights at different times t. In the analogy between the three-and two-dimensional problems detailed in § 3, the domain height has been associated with the plate span s = 2H ( figure 4). Consequently, for a given incoming flow with displacement thickness δ 1 , the dimensionless distance between the plate and the top of the domain is to be constant and equal to H * . It means that the metric terms in the operators (4.3) can be taken at t = 0, which makes the problem both simpler and closer to the analysis of § 3.
We first present the results of numerical simulations corresponding to a domain height H * = 90 and Reynolds number Re = 200. For these computations, 200 collocation points have been used in x and 60 points in y. The wall velocity has been increased progressively up to U * ⊥ = 0.5. For the graph Ψ (x) considered, the theoretical potential outer solution can be computed integrating numerically (3.3) and the corresponding U * e is then obtained by taking the x-derivative. The graph is shown in figure 7(b) together with the Navier-Stokes streamwise velocity at y = H * . The idealized solution (3.5) for the Heaviside step function is depicted as well. As can be seen, besides transient regions near inflow and outflow of length approximately 30, the outer flow taking into account the graph Ψ fits well the approximation (3.5). In particular, the acceleration of the numerical outer flow far from the moving wall is seen to be almost identical with the theoretical prediction in the region 50 x 250.
The computed velocity profiles at x = 127 are compared for different wall velocities in figure 8. The streamwise velocity exhibits a gradient near the wall which increases with U * ⊥ leading to an enhancement of skin friction. The wall-normal velocity exhibits, outside a small region near the wall, the expected almost linear decrease to the imposed zero value at H * = 90. For the numerical steady states the dimensionless shearing stress along the plate is and the result for H * = 90 at small wall velocities U * ⊥ = 0.05, 0.1 is shown in figure 9. Near outflow for x > 280 there is an overshoot in the skin friction values when U * ⊥ = 0, which is a finite-domain effect associated with the gradient of Ψ (x). It appears that whenever U * ⊥ > 0 the flow experiences higher skin friction than for a fixed wall. In this latter case, using the theoretical expression of the boundary-layer displacement thickness δ 1 (x) = γ νx/U , with γ = 1.7208, for the Blasius similarity solution, the skin friction is τ (x) = µcU ∞ /δ 1 (x) with c ≈ 0.57 (Schlichting 1979). The local Reynolds number can be expressed in terms of the Reynolds number Re (based on the displacement thickness at inflow x = 0) used in the computations and the dimensionless skin friction as predicted by boundary-layer theory is where x is the dimensionless distance from inflow. In figure 9 the skin friction c f for U * ⊥ = 0 as resulting from the Navier-Stokes solution is shown as the solid line. The graph (4.8) is plotted as well (with circles) and the two curves are seen to be very close. The small difference is due to inflow effects, because the imposed Blasius profile at inflow is not an exact solution of the Navier-Stokes equations.

Skin friction scaling for the moving plate
The theoretical drag force prediction given by (3.6) is obtained when the potential outer flow is approximated by (3.5), which corresponds to a constant slope U ⊥ /H for x > 0 and the Heaviside step function for Ψ in (3.3). For the graph Ψ used in the numerical simulations (figure 7a), the integral in formula (3.3) is solved numerically and U e as well as its derivative U e can be computed (an example being shown in figure 7b). Three different domain heights H * = 30, 60 and 90 have been considered and the inverse of the normalized outer flow gradient (U e (x)/U ⊥ ) −1 is shown in figure 10 (the outflow region x > 280 being discarded). The curves with an almost singular behaviour near inflow progressively tend towards a constant value, which for the different heights is approximately H * . To take into account the transient region near inflow when comparing theory and the Navier-Stokes computation, in the theoretical drag force formula (3.6) we rather use dU * e /dx instead of U * ⊥ /H * , that is For the three heights considered, steady states have been computed for three different wall velocities U * ⊥ = 0.1, 0.3 and 0.5. Figure 11 shows the computed skin friction, in comparison with the theoretical result (4.9). The theoretical predictions meet the computed values at some distance from inflow which roughly scales with the domain height H * . For the smallest wall velocity, U * ⊥ = 0.1, the theory and computations fit almost perfectly, besides a transient region near inflow. For higher wall velocities (U * ⊥ = 0.3 and 0.5), the transient behaviour near inflow is underestimated by the theoretical curves (in particular, for the largest height H * = 90), whereas more downstream the shearing stress is slightly overpredicted. Overall, the theoretical formula (4.9) may be considered as reliable, in particular for lower heights of the domain.
Similar computations have also been performed at Re = 600, with the same heights H * and velocities U * ⊥ . To summarize the different results, the drag coefficient has been computed according to the formula (3.8), by writing τ = c f (x)ρU 2 with c f (x) the skin friction as resulting from the Navier-Stokes computations. Taking again dU * e /dx instead of U * ⊥ /H * , the two-dimensional drag coefficient is given by Re ∂u ∂y (x, 0). (4.10) Figure 12 shows the evolution of C 2D (x) along the plate. After a transient whose length scales like H * , the drag coefficient converges towards a mean value C 2D ≈ 2.1, for all values of H * and U * ⊥ , and whether Re = 200 or Re = 600. This clearly demonstrates the reliability of a constant drag-coefficient approximation and consequently the scaling (3.11) proposed by Lighthill (1971) for the idealized case in the absence of transient effects.

Time-dependent solution along the moving plate
In the previous section, the skin friction computations assumed the coupled fluid-plate system to be in a quasi-steady state. To characterize the transient regime when the plate is set into motion, a time-dependent Navier-Stokes simulation has also been performed. For this purpose, the time-dependent Navier-Stokes integration procedure described by Marquillie & Ehrenstein (2003) is used, and is adapted to account for moving boundaries, as in Gobert et al. (2010). In this latter work, the coupling between the fluid and an elastic wall was considered. Details of the procedure can be found in these two references. The discretization uses fourth-order finite differences in x, Chebyshev collocation in y, and a semi-implicit time marching. A fractional time step procedure ensures a divergence-free velocity field. The procedure has been designed for boundary-layer type computations and it uses a simplified mapping y =ȳ + ϕ(x, t) with a fixed lower boundaryȳ = 0 in the transformed domain 0 ȳ H * . This simplification is justified for boundary-layer flows, provided that the domain is sufficiently high to recover an outer flow field at H * independent of the vertical coordinate. This is true for the streamwise velocity u given the condition ∂u/∂y = 0 at the top of the domain and the gradient with respect to y of v is approximately U * ⊥ /H * , which for the case considered hereafter (H * = 60 and U * ⊥ = 0.3) is small. A large domain 0 x 920 has been considered and the moving boundary is again ϕ(x, t) = U * ⊥ tΨ (x) with a graph Ψ (x) similar to that shown in figure 7(a), such that Ψ (x) = 1 in the region 30 x 600 (in other words, a relatively large outflow buffer region of length 300 has been considered to avoid upstream effects during the time integration).
According to the mapping used for the time-dependent simulation, the transformed gradient and time derivative are (4.11) The moving boundary ϕ(x, t) = U * ⊥ tΨ (x) is uniform in the domain 30 x 600, but its streamwise gradient becomes increasingly stiff in the buffer region when progressing in time. This difficulty is overcome in the time-dependent computation by simply setting ∂ϕ(x, t)/∂x to zero in (4.11) and owing to the mapping y =ȳ + ϕ(x, t), the upper boundary thus remains at a constant height H * from the ascending plate. This is precisely the quasi-steady configuration considered in § 4.1. The timeintegration has been performed at Re = 200 with 6150 discretization points in x and 80 collocation points in y.
First, a time-integration has been performed for the flow along the flat plate with zero wall velocity until a steady state is reached. Using this boundary-layer solution as the initial condition, the plate is set into motion with the velocity U * ⊥ = 0.3 and the skin friction along the plate at successive times is shown in figure 13. The equilibrium 338 U. Ehrenstein and C. Eloy  state result (for the smaller domain 0 x 280) is shown as well and the transient regime in time can be estimated to last until t ≈ 60. For large times, due to the time-independent forcing term ∂ϕ/∂t = U * ⊥ Ψ (x) with a sharp gradient at inflow, the solution will ultimately exhibit an asymptotic instability during the time-integration procedure.
According to these time-dependent simulations, the increase in skin friction predicted by the steady-state analysis appears to be a lower bound reached only after a transient. The question is now to understand to what extend this skin friction enhancement prediction is retrieved for a more realistic swimming behaviour, such as an undulatory rather than uniform motion.

Skin friction for an undulatory swimming motion
Consider a flat rectangular plate of aspect ratio s/L = 1/5 performing an undulatory motion similar to a fish swimming at velocity U ( figure 14) y(r, t) = Ar cos (k(r − Vt)) , where r is the curvilinear coordinate along the plate length, A is the dimensionless amplitude of beating (A = 0.1 in this example), V is the wave speed of the undulations (V = 1.5 U) and k is the wavenumber (k = 2π/L). These values give a Strouhal number, St = AVkL/πU = 0.3, typical for a large number of fish and cetaceans with similar aspect ratios (Eloy 2012).
In the linear limit, the normal velocity of the plate with respect to the fluid is thus varying from ∼0.1 U to 0.3 U from head to tail. The longitudinal drag per unit length due to this normal velocity is given by (2.11), which can rewritten using the fact that = 2s/π for a flat plate and that C 3D ≈ 1.8, as where Re s = Us/ν is the Reynolds number based on the span. This drag force can be compared with the Blasius drag on the same flat plate when it is motionless These two drag forces are of same order when s/x ≈ 4 U ⊥ /U ≈ 0.4, which means that the drag will be dominated by the effects of the normal velocity on the last 60 % of the plate. A more precise comparison can be carried out by taking into account the nonlinear geometrical effects and by calculating the time-average drag, as illustrated in figure 15. In this figure, it is seen that indeed the effects of the normal component of the velocity are dominant on the last 56 % of the plate. The total drag can be estimated by taking the time and space average of the maximum between D Blasius and D . It results in a 23 % increase compared with the Blasius drag alone that would act on a motionless plate. Note that this drag enhancement is not only due to the 'thinning' of the boundary layer but also to the acceleration of the tangential component of the outer velocity U , as was originally remarked by Anderson et al. (2001).
To asses numerically the validity of this drag increase, we again consider the analogy between the three-and two-dimensional problem (cf. figures 4 and 6), the height of the two-dimensional domain H being half the plate span s. To mimic the aspect ratio between width and length of the plate, a two-dimensional domain of length 10H is to be considered for the Navier-Stokes computation. The displacement thickness δ 1 of the incoming Blasius profile is again considered as reference length and for the same domain as in § 4.3, that is H * = 60 and 0 x 920, the time-dependent Navier-Stokes solution has been computed for a dimensionless wall velocity with k = 2π/L * and L * = 10H * = 600, which corresponds to the case considered above in (5.1). Again a buffer region is introduced through the graph Ψ (x). The inhomogeneity in space associated with the plate motion is discarded to account for the analogy between the theoretical formula (5.3) and the two-dimensional unsteady computation, that is in the mapping (4.11), ∂ϕ/∂x = 0 and ∂ϕ/∂t = U * ⊥ (x, t). The period of the plate motion is T = 400 and the time-dependent Navier-Stokes integration has been performed, starting with the motionless flat plate steady state as initial condition, for the Reynolds number Re = 200. The equivalent Reynolds number Re L based on the plate's length is of order 10 5 for which three-dimensional Navier-Stokes computations for the undulatory geometry with finite end effects would be hardly possible. Note that this corresponds to a typical Reynolds number for swimming fishes with length 20-30 cm.
After a transient time a periodic regime is retrieved and the solution has been computed for several cycles. The plate velocity and the corresponding instantaneous skin friction are shown in figure 16 at two instants. As expected, skin friction enhancement as well as skin friction reduction (with respect to a motionless plate value), are observed along the plate. The time average of the skin friction c f has been computed for different periodic cycles and the result is shown in figure 17. The computations for two successive cycles almost superimpose, which confirms the periodic behaviour of this Navier-Stokes solution. For x > 150 the mean skin friction lies above the value for the motionless plate and a 35 % increase is predicted further downstream. Integrating the shear stress along the plate, one finds an increase with respect to the Blasius drag of roughly 20 % that confirms the scaling obtained above with a simple analytical approach based on formula (5.3).

Discussion
In this paper, we have addressed how the skin friction on an object is modified when the outer flow velocity is no longer parallel to its surface but has a normal component as well. A classical example is the calculation of the longitudinal drag on a yawed cylinder immersed in a uniform flow, first studied by (Taylor 1952). We have shown that, in this ideal case, the boundary layer thickness scales as √ νs/U ⊥ , where s is the span and U ⊥ the normal component of the outer velocity, as it was originally proposed by Taylor (1952). In this context, a drag coefficient has been defined, and it has been calculated to be approximately equal to C 3D ≈ 1.8 for elliptic cylinders. This longitudinal drag is intimately due to the acceleration of fluid particles as they move around the cylinder. As such, it is inherently a three-dimensional effect. However, an analogue two-dimensional problem can be conceived if the flow is artificially accelerated along the longitudinal direction. We have proposed such a flow, where the acceleration is achieved by a channel of finite height, and have shown that a similar scaling can be obtained. In this two-dimensional problem, a drag coefficient similar to the three-dimensional case can be calculated using boundary-layer theory and we have found C 2D = 2.1. These results confirm the 'Bone-Lighthill boundary-layer thinning hypothesis' (Lighthill 1971) in the ideal case of a uniform and steady problem.
These analytical results, based on the boundary-layer approximation with an outer potential flow varying slowly in the streamwise direction, have been shown to be relatively robust by comparing them with steady-state solutions of the full U. Ehrenstein and C. Eloy Navier-Stokes equations. Transient effects, due to non-uniform acceleration of the potential outer flow near the leading edge of the moving plate, are however underestimated by the theoretical prediction, in particular for the larger values of the domain height, corresponding to large spans in three dimensions. In addition, time-dependent computations have shown that transient effects during the upward motion against the incoming boundary-layer flow enhance skin friction, in comparison with the quasi-steady state which is reached only after some time. Generally speaking, Navier-Stokes simulations have shown that non-uniform wall velocities and transient effects tend to slightly enhance skin friction in comparison with the uniform and quasi-steady case assumed in the theoretical analysis.
In an attempt to generalize these results obtained for uniform plate motions to more realistic swimming motions, an undulatory plate velocity mimicking the swimming motion of a fish has been used in the theoretical formula and a 23 % increase of skin friction drag has been predicted. The corresponding two-dimensional numerical simulation exhibits instantaneous enhancement or reduction of the shearing stress during the undulating cycle, but, on average, an enhancement of the skin friction drag of the same order of magnitude is observed.
In summary, the simple three-dimensional example chosen and the two-dimensional numerical simulations have shown that one can expect a drag enhancement due to the motion of the swimming animal. This confirms the hypothesis of 'boundarylayer thinning' proposed by Bone and Lighthill (Lighthill 1971). In practice, this enhancement will closely depend on the particular geometry and on the motion. In our computations the theoretical value of skin friction increase for a model of undulatory motion is 23 % and the same order of magnitude is retrieved in the time-dependent two-dimensional computation. This increase is below the measurements of Anderson et al. (2001) on the scups where they find an increase by a factor of 1.5-1.9. Even though our model predictions are to be considered as a lower bound, skin friction increases by factors between 4 and 10, as proposed by Lighthill (1971), Webb (1975), Alexander (1977) and Videler (1981), seem however unlikely.

Appendix A. Momentum equation
The boundary-layer equations (2.10a)-(2.10c) are classically solved by using an integral form (Schlichting 1979). To do so, the displacement and momentum thicknesses are defined (all variables are now dimensionless and the asterisks have been dropped for simplicity) where T * = (T − T w )/(T ∞ − T w ) is the dimensionless temperature (with T w and T ∞ the temperatures on the cylinder wall and at infinity, respectively) and Pr = ν/α is the Prandtl number based on the thermal diffusivity α. With this dimensionless temperature, the boundary conditions are exactly the same as for the longitudinal velocity u * x , i.e. T * = 0 in η * = 0 and T * = 1 for η * → ∞. Therefore, for Pr = 1, the two problems are analogous.
One can define an analogue of the drag coefficient C 3D for the heat transfer problem which satisfies where Nu is the average Nusselt number based on the characteristic length and C 3D now also depends on Pr. There is a large amount of literature devoted to the correlation between Nu, Re and Pr for circular cylinder. The state of the art probably corresponds to the Churchill-Bernstein empirical equation (Churchill & Bernstein 1977) Nu = 0.3 + 0.62Re 1/2 Pr 1/3 1 + (0.4/Pr) 2/3 1/4 1 + Re 282000 5/8 4/5 , (B 3) which is valid over a wide range of Reynolds and Prandtl numbers, as long as PrRe > 0.2. For Pr = 1 and Reynolds numbers typical of swimming animals, 10 3 < Re < 10 6 , the Churchill-Bernstein equation simplifies to Nu ≈ 0.556Re 1/2 , (B 4) which gives C 3D = 1.75, remarkably close to the calculation detailed in appendix A that gives C 3D = 1.78 for b/a = 1. There have been relatively few studies of forced heat transfer on elliptical cylinder. Ota & Nishiyama (1984) andŽukauskas &Žiugžda (1985 performed experiments with air, for which Pr = 0.7, and the results are compared with the present boundarylayer calculation in figure 18. While there is a reasonable agreement for a circular cylinder (i.e. b/a = 1), the variation with the ellipse aspect ratio b/a is not recovered. We do not have any explanation for this discrepancy.
It should be noted that Khan, Culham & Yovanovich (2005) performed a boundarylayer calculation similar to what is presented here for b/a 1, and in the limit of asymptotically large Pr. In their analysis, they find that C 3D increases when b/a decreases, in agreement with the experiments of Ota & Nishiyama (1984) anď Zukauskas &Žiugžda (1985). However, their calculation should not be valid in the regime of the experiments for which Pr = 0.7, and the present analysis should in principle give a much more accurate result.  (1985) for which 10 4 < Re < 10 5 .