Uniform flow past a periodic array of cylinders

e


Introduction
Uniform flow past a circular cylinder is a fundamental flow in fluid mechanics, constituting a basic paradigm for aerodynamics, geophysics, and porous media flows, to name just a few application areas. The solution for inviscid, irrotational flow of speed U along the x-axis past a circular disc of radius a centred at the origin is often given in the form of the following complex potential for the uniform flow past N > 1 cylinders, in general configuration, has only recently been written down. In [3] a new "calculus" for constructing the complex potentials for general ideal fluid motions in multiply connected domains was described; this includes the potentials for uniform flow past an arbitrary finite assembly of circular cylinders [4]. The article [3] gives the complex potentials for such uniform background flows as first written down in [4] in terms of the so-called Schottky-Klein prime function [5].
Of no less interest in applications is the solution for uniform flow past a periodic array of cylinders, but this important problem is surprisingly resistant to known mathematical techniques. There are two canonical situations: uniform flow past the cylinders in tandem, or the cross-flow case. These are illustrated schematically in Figure 1. The linearity of the underlying boundary value problem means that the solution for a general uniform flow past this periodic cylinder grating follows immediately on superposition.
If one restricts consideration to the flow in a single period window then the restricted flow domain is doubly connected and the solution is available, in principle, by combining the aforementioned calculus [3] with some conformal mapping ideas. One option is to make use of recent results for constructing the conformal mappings to a doubly connected polycircular-arc domain from a preimage circular region [12]. This approach nevertheless involves solving an accessory parameter problem which is generally nonlinear and not completely straightforward. Another option is to use the symmetry of the geometrical configuration to reduce consideration to a quarter of the period window; this has the effect of rendering the domain of interest simply connected. However, use of a circular-arc conformal mapping to a simply connected canonical domain, such as the unit disc, is still needed and presents its own complications. Indeed it was precisely to avoid such matters that Richmond [17] devised an ingenious approximation to this conformal mapping, one that gives a good approximation to the required solution provided the radius of the cylinder is sufficiently small compared to the length of the period window. A third approach involves finding instead the inverse conformal mapping from the doubly connected flow domain to a preimage concentric annulus, say. Actually, the latter approach turns out to be a linearizable problem and we will say more about it later.
The purpose of this paper is to present an alternative approach to solving both the tandem and cross-flow problems without the need for any conformal mapping. The solutions are derived by means of a novel transform technique recently derived by the author [8,9]. The new technique itself is very gen- eral and applies to boundary value problems for harmonic fields in (possibly multiply connected) circular domains, i.e., domains whose boundaries are circular arcs and straight lines. The new transform method formulated in [8,9] is an extension (to circular, multiply connected domains) of a novel transform approach to harmonic boundary value problems in convex polygons first presented by Fokas & Kapaev [13]. All these transform techniques fit into a much broader collective known as the unified transform method pioneered by A. S. Fokas [14] (now often called the Fokas method) and which applies not only to linear boundary value problems but also to initial value problems for nonlinear integrable partial differential equations.
Since the general transform technique described in [8,9] has only recently been formulated it is instructive to fully document the details of its application to the canonical flow problems of Figure 1. The analysis of this paper is expected to point the way to application of the same transform method to many other flow scenarios, including harmonic boundary value problems in other physical science applications outside fluid mechanics.
Lamb [18] (see §64) has given an approximate solution for the velocity potential in the cross-flow problem for the case of small cylinders (relative to their periodic spacings): The quantity Λ, which has dimensions of length, is often referred to as the blockage coefficient. These blockage coefficients are important physically and they often emerge as the only quantities about these flows needed in applications (see [7,19,15] for a more general discussion, and additional references). The value of the blockage coefficient for the cross-flow problem furnished in (2) is its value in the dilute limit (small cylinders compared to their spacing); from our transform solution we are able not only to produce the blockage coefficients in arbitrary configuration but also to discern a convenient approximate formula for Λ, one that reduces to (2) in the dilute limit but which is accurate for a much larger range of cylinder spacings. An analogous approximate formula for the blockage coefficient in the tandem problem is also found. Concerning the previous literature on flows in this geometry, Balsa [1] has studied uniform irrotational ideal flow past cylinder arrays using the methods of matched asymptotic expansions, assuming the cylinders are small compared to their typical separation, and with particular interest in questions of stability. The cross-flow flow problem considered here can alternatively be viewed as flow past a cylinder between two plane channel wallsa problem that has been studied numerically, and at finite (rather than infinite) Reynolds numbers, by Chen, Pritchard & Tavener [2] and Zovatto & Pedrizzetti [20].

Mathematical formulation
Consider the complex potential where φ is the velocity potential and ψ is the streamfunction for the irrotational two-dimensional flow [16]. We require to find this function for both canonical problems (tandem and cross-flow).
The 2l-periodicity of the flows allows us to restrict to the domain D: a single period window comprising the channel region −∞ < x < ∞, |y| < l with |z| > 1. In both tandem and cross-flow problems the boundary of the cylinder in the period window must be a streamline. On the edges of the period window either φ is constant (tandem problem) or ψ is constant (cross-flow problem). The far-field boundary condition for the tandem flow with unit far-field velocity in the positive y-direction is that as |z| → ∞. For the cross-flow problem with unit far-field velocity along the as |z| → ∞.
To proceed it is convenient, when studying the cross-flow problem, to consider the modified complex potential Condition (6) implies the following far-field condition on h(z): which is now identical to (5). A comparison of the boundary value problem (restricted to a period window) for w(z) = φ + iψ in the tandem problem, and that for h(z) ≡ −iw(z) = Φ + iΨ in the cross-flow problem, reveals that the only difference between them is the nature of the boundary condition on the cylinder; compare Figures 2 and 4. Henceforth the aim is to find w(z) for the tandem problem and h(z) for the cross-flow problem.

Flow symmetries
The following arguments pertain to both problems, although we present in detail our treatment for the boundary value problem determining w(z) (Figure 2), with the minor modifications for the boundary value problem determining h(z) (Figure 4) stated afterwards.
The boundary conditions at the edges of the period window are   for some constant c and The symmetries of the flow are such that the velocities are invariant as z → −z. Hence we conclude that that is, w(z) is an odd function. From this property we see immediately that if, as x → ∞, for some constant d then, as x → −∞, implying that where the Schwarz conjugate is defined as w(z) ≡ w(z).

Motivated by these observations we introduceŵ(z) via the relation
where λ ∈ R (the blockage coefficient) and the functionŵ(z) remain to be determined. w(z) satisfies the far-field conditions (12) and (13), with d = iλ, provided thatŵ(z), taken to be analytic and single-valued in the flow domain, tends to zero as |z| → ∞.
On y = −l, where z = z + 2il, we have, from (9), On use of both (11) and (15) we find or, on substitution of (78),ŵ provided that we pick c = l. It is straightforward to confirm that the above considerations are consistent with the boundary condition (9) on y = l. Relation (19) means thatŵ(z) is 2il-periodic. Owing to this periodicity property the values ofŵ(z) at points on the two sides of the channel displaced by 2il are equal. For use later we introduce the function which denotes the value ofŵ(z) on the lower and upper boundaries.

Integral transform representation
We now import, without proof, specific results from a general transform scheme for harmonic boundary value problems in multiply connected circular domains recently presented by the author [8,9]; the reader is referred to those references for full details of the derivation of the following transform pairs. Here we simply show how to use those transform pairs to solve the problems at hand.
The domain D is a doubly connected circular region (its boundaries are made up of straight lines and circles). Since the correction termŵ(z) decays as |z| → ∞, it can be shown [8,9] that it is possible to write it, for z ∈ D, in the form of an integral representation: where the set of contours {L j |j = 1, 2, 3} are shown in Figure 5 [8,9]. The three functions of k appearing in this formula are defined to be and These are known as spectral functions. The remaining elements of the socalled spectral matrix ρ ij (k) [8,9] are given by and and ρ 31 (k) = ρ 11 (k) and ρ 13 (k) = ρ 33 (k). The functions appearing in the spectral matrix have additional analytical structure; they satisfy a set of so-called global relations. They are which are equivalent, and Both (26) and (27) must be analyzed to find the unknown spectral functions. While we have not given the derivation of the global relations (26) and (27) -interested readers should consult [9] for that -they are very natural and are, in fact, direct consequences of Cauchy's theorem. Inspection reveals that both global relations are consequences of the analyticity ofŵ(z) in D.
Of course, the spectral functions ρ 11 (k), ρ 22 (k) and ρ 33 (k) are not known a priori and must be found. In §6 and §7 the global relations (26) and (27) will be analyzed to determine them. Once found, they can be substituted into (21) to give the required integral representation ofŵ(z).
For the cross-flow problem we will write where we have changed the designation of the unknown blockage coefficient to Λ. An integral representation of the correction termĥ(z) analogous to (21) will then be found.

Tandem and cross-flow problems
Everything discussed so far pertains to the boundary value problems for both w(z) and h(z). We now discuss the differing conditions on |z| = 1.
For the tandem problem on |z| = 1 we know from (10) and (78) that Hence, on |z| = 1, we writê where H(z) is the unknown real part ofŵ(z) and Therefore, for the tandem problem we will write where the set of real coefficients {a n |n ≥ 1} are to be determined. The ansatz (32) satisfies (15). For the cross-flow problem on |z| = 1 we know from (10) and (78) that Hence, on |z| = 1, we writê where H(z) is now the unknown imaginary part ofŵ(z) and It should be noted that, for notational convenience in the next section, we have used H(z) to denote the unknown boundary data for both tandem and cross-flow problems. In contrast to (32), for the cross-flow problem we will write where the set of real coefficients {b n |n ≥ 1} are to be determined. The ansatz (36) satisfies (15).

Analysis of global relations (tandem)
In this section we analyze the global relations for the tandem problem. Analogous details for the cross-flow problem differ only slightly and are discussed in the next section.
The first of the global relations (26) gives The second independent global relation (27) takes the form where is recognized as the standard Fourier transform of G(x), and with It follows that The inverse Fourier transform implies that It is clear that we require in order to remove the simple pole singularity of the integrand at k = 0. The second global relation (A.15) implies that, for n ∈ N, It is easy to show that for n = 2m with m ≥ 1, (46) then implies that (49) On substitution of (A.20) for G(x), and on letting m → n, we find where An exercise in residue calculus reveals that for n ≥ 1, On substitution of (A.28) into (A.25), and after a change of integration vari-able so that all integrals are over the positive real k-axis, (A.25) becomes But so that The final system of equations for the coefficients {a n |n ≥ 1} and λ is where the first equation is the condition (45), and With λ and the coefficients {a n |n ≥ 1} determined as the solution to this system the spectral functions are given by These can be substituted into (21) thereby completing the solution.

Analysis of global relations (cross-flow)
It is easy to trace the minor changes of the preceding analysis required when treating the cross-flow problem: we find that the analogue of (55) acquires a minus sign: The only other changes are the redefinitions of R 1 (k), R 2 (k), L 1 (k) and L 2 (k) to the identical forms but with r(z) →r(z) and s(z) →s(z).
The final system of equations for the coefficients {b n |n ≥ 1} and Λ is then where the first equation is the condition (45), and Figures 7 shows typical streamline plots, computed by means of the integral transform solution, for the cross-flow problem in the case l = 1.2.

Uniform flow in general direction
The complex potential W (z) for a general uniform flow, where as |z| → ∞ for any given U, V ∈ R can now be obtained by a superposition of the two solutions just obtained. Indeed, or, alternatively, 9. Approximations to blockage coefficients Figure 8 shows graphs of the numerical solutions for |a 1 | and |b 1 | as functions of a = 1/l for a ∈ (0, 1). It is clear that for a-values less than 0.7 the values of |a 1 | and |b 1 | are uniformly small, indeed their maximum values occur at a = 0.7 and are given by a 1 = −0.0011 and b 1 = −0.0028. Even for a = 0.9 these values only rise to a 1 = −0.0063 and b 1 = −0.0534. The important thing, though, is the magnitude of these quantities relative to the size of the two other terms in the first equations in (A.31) and (60). It turns out that for all a-values less than around 0.9 they are relatively small so that, to a good approximation, we can set where, from (A.18) and (31) or (35), These contour integrals can be computed explicitly yielding the approximate formulas λ ≈ (π/l) 2 + (π/l) 2 /6 , In the limit 1/l → 0, to leading order we retrieve the dilute result (2) found by Lamb, i.e., However (69) provide continuations of Lamb's linear approximation into the nonlinear régime (i.e., 1/l not necessarily small). Figure 9 shows a graph of Λ for the cross-flow problem (solid line) with a superposed graph (shown dotted) of the approximation formula (69). For a in the interval (0, 0.9) formula (69) gives a very accurate approximation of Λ. Also shown in Figure 9 are the values given by an approximate conformal mapping method due to Richmond [17]. As reported in [15], Richmond's approximation is where β is a solution to the nonlinear equation Richmond's approximation, which clearly differs from ours, is also accurate for a good range of a-values although it starts deteriorating earlier at a ≈ 0.7. This is not surprising: Richmond [17] used an approximation of the conformal mapping to a quadrant of the period window exterior to the circular arc and showed that for a = 0.75 that mapping approximated a circular arc only to within an 11% error. Our approximation for Λ in (69) is an improvement on Richmond's and has the added virtue of being an explicit (rather than an implicit) formula.
Finally, it is worth mentioning that an exact mathematical analogy has recently been pointed out [7] between the blockage coefficients considered here and the so-called slip lengths relevant in many microfluidics applications for characterizing the slip properties of superhydrophobic surfaces. Our approximate formulas (69), and our general methods, will also be useful in those quite different fluid dynamical applications (see [9]).

Conformal mapping approach
We can provide a consistency check on the foregoing analysis using conformal mapping ideas akin to those underlying the "calculus" for ideal flows in general finitely connected domains advanced by the author in [3,4] and described in the Introduction. In that approach, if the conformal mapping z = z(ζ), say, from a circular preimage ζ-domain to a fluid domain of interest can be found, the complex potentials for an array of ideal flows in that domain can be written down explicitly in terms of a special transcendental function called the Schottky-Klein prime function [3,5,11]. For the doubly connected fluid region D studied in this paper it is natural to consider a conformal mapping z = z(ζ) from the concentric annulus in a parametric ζ-plane. The parameter ρ must be determined from the shape of D. The latter is an example of a doubly connected polycircular arc domain and, in principle, the method expounded by Crowdy & Fokas [12] could be used to find z = z(ζ). However, we will not pursue such an approach since the inverse function ζ = ζ(z) is more useful, as explained shortly. (Note that here we are adopting the convention of using z and ζ to denote both the coordinate in the respective complex planes and the conformal mapping function). The Schottky-Klein prime function associated with the annular domain (73) is (to within a normalization) P (ζ, ρ) satisfies the two functional relations The crosses show values given by an approximate conformal mapping method due to Richmond [17] as reported in [15].
which can be checked directly from (74). It can be shown that the complex potential for both uniform flow problems considered in this paper can be written explicitly in terms of the special function P (ζ, ρ) as follows: , (tandem).
(77) The cross-flow solution is just the logarithm of a relevant conformal slit mapping [3,11]; the tandem solution can be obtained by straightforward adaptation of the results in [10]. Regardless of their derivations, it is an easy exercise to use (75) to verify directly that (77) satisfy all the required conditions.
Formulas (77) are only useful if the conformal mapping from the annulus (73) to D is known or, even better, its inverse ζ = ζ(z), which is a conformal mapping from D to the annulus (73). If ζ = ζ(z) can be found then (77) gives explicit expressions for the required complex potentials. But, remarkably, the mapping function ζ = ζ(z) can itself be found using the new transform approach expounded in this paper. That derivation is given in the appendix. The resulting function ζ = ζ(z) can then be used, in combination with formulas (77), to provide a non-trivial check on the transform formulas for the complex potentials derived in §6 and §7.

Discussion
This paper has given a quasi-analytical mathematical solution to two canonical ideal flow problems: the uniform flow past a grating of cylinders in tandem or in cross-flow. By superposition, the solution for general uniform flow past the cylinder grating follows. The solutions are derived from a new transform technique for Laplace's equation in multiply connected circular domains described in [8,9]. The general approach is extendible to many other flow scenarios. Useful approximations for the blockage coefficients emerge from the analysis.
It is worth remarking that had the problem of interest involved just a boundary value problem for Laplace's equation in a channel it would have been natural to apply Fourier transform methods. On the other hand, if the problem had involved a boundary value problem around just a circular disc, with no channel walls, we might have used Taylor/Laurent expansions about the circle center or, perhaps, Mellin-type transforms [8]. But on placing a cylinder inside a channel we intuitively need some hybrid of these two methodologies. But this is exactly what our general transform scheme provides. The solution for w(z) has been derived in the form where we now rewrite (21) aŝ From this expression it is clear the our transform method produces exactly the hybrid scheme we seek: the integral representation (79) contains terms of "Fourier transform type", with associated spectral functions ρ 11 (k) and ρ 33 (k), as well as terms of "Mellin transform type", with associated spectral function ρ 22 (k). All these spectral functions are inter-related via the global relations (26) and (27). It is these global relations that we have analyzed here to find the spectral functions and, hence, to determine the solution.
Finally, an alternative approach to the problems considered here is via standard numerical boundary integral methods. That typically leads to integral equations for a set of unknown boundary data the solution of which produces integral expressions for the flow field variables akin to the integral representations for w(z) (in terms of spectral data) found here. The boundary integral method is, of course, intimately related to the transform method presented here with the place of the integral equations arising in the former scheme taken by our small linear systems for the unknown coefficients of the function we have, in both calculations, called H(z). The solution of those small linear systems requires almost no numerical effort and, in this respect, our approach is certainly competitive with alternative purely numerical schemes.
whereχ(z) → 0 as x → ±∞. The first term is simply the logarithm of the conformal mapping of the channel (without the hole) to a unit disc; the latter is easily derived using elementary considerations. Since the boundary conditions on χ(z) are then, on use of (A.2), the following boundary conditions onχ(z) pertain: Recall that ρ is not known in advance and must be found. By the symmetries of the proposed mapping between regions we expect that if a point ζ on the upper-half unit circle corresponds to z = x + il then the point ζ will correspond to z = x − il. This means that, for each x, implying the relation χ(z + 2il) = −χ(z). for some (purely imaginary) function G(x). On |z| = 1 we will writê where r(z) = log ρ + log coth and where the coefficients a 0 ∈ R, {a m ∈ C|m ≥ 1} are to be found. We also expect, on grounds of symmetry, implying that r(z) − iH(z) = r(z) + iH(z), or H(z) = −H(z). (A.12) This condition implies a 0 = 0 and a n = ib n , (A. 13) for some real set {b n }. We will make use of these facts later.
The functionχ(z) is single-valued and analytic in the fluid region D. It therefore has a transform representation of the form (21). As before, we now analyze the associated global relations. (26) and (27) give, respectively, It can be checked easily that E 0 is the only quantity that depends on the unknown log ρ. We can therefore solve (A.34) for the set of coefficients {b n |n ≥ 1} and then use (A.33) a posteriori to determine log ρ.