Problems in fluid-structure interaction http://www.damtp.cam.ac.uk/user/pdm23/papers/smith.pdf

Paul D. Metcalfe

ii

Contents 1 Introduction 2 The basic problem 2.1 The Briggs–Bers method . . 2.2 Absolute instability . . . . . 2.3 Other effects of causality . . 2.4 Variants and extensions . . . 2.5 The elastic channel and pipe 2.5.1 Inviscid perturbation

3 . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

3 Ribbed elastic structures 3.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 The Green’s function . . . . . . . . . . . . . . . . . . . . 3.2 Reduction of Green’s function matrix . . . . . . . . . . . . . . . 3.2.1 Transfer matrix form . . . . . . . . . . . . . . . . . . . . 3.2.2 Reconstruction of solution, energy flux and other useful tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Branch line contribution . . . . . . . . . . . . . . . . . . 3.3 Periodically ribbed membranes . . . . . . . . . . . . . . . . . . 3.4 Disordered rib arrays . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Ribbed plates . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 The finite baffle and other problems 4.1 Numerical scheme . . . . . . . 4.2 Some results . . . . . . . . . . . 4.2.1 The flag problem . . . . 4.3 The lattice Boltzmann method .

1

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

5 . 6 . 7 . 8 . 9 . 11 . 13 . . . .

15 16 17 17 19

. . . . . .

19 20 21 22 26 27

. . . .

33 35 37 37 39

2

CONTENTS

Chapter 1

Introduction Many unusual phenomena are found when a fluid mechanical system with mean flow is coupled to a solid system; the typical example is flow in a system with elastic boundaries. Similar problems are a source of much current research interest, whether in biomechanics (Luo & Pedley, 2000) or in the vibrational problems which are the focus of this essay (Crighton, 1989). The basic problem in this field is that of (inviscid, irrotational, incompressible) flow over an infinite elastic plate (or membrane). Recent interest in the wave mechanics of this system dates back to numerical work by Brazier-Smith & Scott (1984), which was explained analytically by Crighton & Oswell (1991). Since this work many variations on this basic problem have been studied (Lucey & Carpenter, 1992; Peake, 1997; Lingwood & Peake, 1999; de Langre & Ouvrard, 1999; de Langre, 2000; Peake, 2000). Much of this work has been done to relax the (somewhat unphysical) assumptions of BrazierSmith & Scott (1984), to see if these phenomena can be found in more and more physically plausible systems. A related problem is that of disturbance propagation in a stiffened fluid loaded elastic structure. Many elastic structures have stiffeners or ribs attached to provide extra rigidity, and this problem has been studied by a number of workers (Rumerman, 1975; Crighton, 1984; Spivack, 1991; Cooper & Crighton, 1998a), but only in the case of zero mean flow. The inclusion of mean flow will be shown to lead to significant new effects. In this essay we first discuss our basic problem, flow over an elastic membrane. Much of this discussion has been adapted from work by Kelbert & Sazonov (1996). In §2.5 we then discuss some new variants of this basic problem. Changing focus slightly, in §3 we discuss disturbance propagation along a ribbed elastic structure. The inclusion of mean flow effects will be seen to permit some unusual new behaviour. This work is closely related to work in modern dynamical systems and solid state physics and these links will be briefly discussed. Finally, in §4 we discuss some largely numerical work on the response of a finite elastic plate and some other related problems. The work in §3 is original, as is the (largely incomplete) work in §4. Much of the work in §2 is not original, although §2.5 is my own. All the original work has been done as a part of my doctoral studies, and is as yet unpublished. I gratefully acknowledge the support of my supervisor, Nigel Peake, and of my industrial supervisor, Roger Kinns of BAe Systems. 3

4

CHAPTER 1. INTRODUCTION

Chapter 2

The basic problem and variations thereof As mentioned earlier, a number of interesting features appear in structural interaction with a nonzero mean flow. The simplest problem that exhibits all of them is the problem of flow over an elastic membrane which is driven by a line force on the line x = 0. The following exposition (§2.0–§2.4) has been freely adapted from work by Cairns (1979), Crighton & Oswell (1991) and Kelbert & Sazonov (1996). We suppose the flow to be incompressible, inviscid and irrotational, so that the fluid flow satisfies Laplace’s equation. Let the mean flow have speed U. On linearising the Bernoulli equation and the boundary conditions we get the set of equations mηtt − Tηxx = −p(y = 0) + δ(x)F(t) 2

∇ φ=0 p = −ρ(∂t + U∂x )φ φy (y = 0) = (∂t + U∂x )η,

(2.1) (2.2) (2.3) (2.4)

where η is the membrane displacement, m is the membrane density, T is the membrane tension, p is the perturbation fluid pressure, φ is the perturbation velocity potential, ρ is the fluid density and F is the time dependence of the line forcing. In order, later, to properly enforce causality, we suppose that F(t) = 0 for t < 0. We also have a boundary condition at infinity on φ, that ∇φ → 0 as y → ∞. These equations are nondimensionalised with a lengthscale m/ρ and a timescale (m/ρ)(m/T )1/2 . The lengthscale comes from finding the depth of fluid needed to balance the membrane mass, and the timescale comes from equating the tension and inertia terms in the membrane equation. On adopting these scalings we regain the same problem but with m = ρ = T = 1 and U 7→ U(T/m)1/2 . We now Fourier transform in streamwise co-ordinate x and time t, so that, for instance ZZ ˜η¯ (k, ω) = η(x, t)eıωt−ıkx dxdt. (2.5) 5

CHAPTER 2. THE BASIC PROBLEM

6

This is equivalent to making the transformations ∂t → −ıω and ∂x → ık. We can solve the resulting equations to get ¡ 2 ¢ ˜ k − ω2 − (ω − Uk)2 /γ η˜¯ = F, (2.6) where γ2 = k2 and Re γ ≥ 0. Therefore, Z Z ˜ 1 F(ω) η(x, t) = eıkx−ıωt dωdk, 2 4π Γk Γω D(ω, k) where

D(ω, k) ≡ k2 − ω2 − (ω − Uk)2 /γ

(2.7)

(2.8)

is the dispersion relation for free waves on an unforced membrane. The two contours of integration Γω and Γk are chosen to satisfy causality requirements. We start by letting Γk lie on the real axis. In order to satisfy causality and produce η = 0 for t < 0, we require Γω to pass above all the sheets of zeroes of D(ω, k) for k ∈ Γk . Then if t < 0 we can close the Γω contour in the upper half-plane and we find η = 0. Note, incidentally, that since Im ω > 0, we can assign a meaning to the temporal Fourier transform of H(t)e−ıω0 t (with ω0 real), as the analytic continuation of the integral Z1 ı eı(ω−ω0 )t dt = (2.9) ω − ω0 0 from Im ω > 0. We really want to consider time-harmonic forcing with F(t) = ˜ = ı(ω − H(t)e−ıω0 t , and in this case (2.9) tells us that we should put F(ω) −1 ω0 ) .

2.1 The Briggs–Bers method We now have an equation for the plate displacement η(x, t) =

ı 4π2

Z

Z D(ω, k)−1 Γk

Γω

eıkx−ıωt dkdω, ω − ω0

(2.10)

and we know that η ≡ 0 for t < 0. This expression, although an exact solution of the problem, is less than useful as a source of physical understanding. We need some way of evaluating it, preferably in terms of some wavelike forms. Supposing that we can deform Γω onto the real axis in the ω plane, we can simply close the Γω contour in the lower half plane. At large time, any poles in the lower half of the ω plane will give decaying waves, so we can asymptotically evaluate (2.10) as e−ıω0 t η(x, t) ∼ 2π

Z Γk

eıkx dk. D(ω0 , k)

(2.11)

This procedure fails if, in deforming the Γω contour, a pole approaches the Γk contour in the k plane. In this case we need to indent the Γk contour around the pole. Provided that we can keep indenting the Γk contour we can

2.2. ABSOLUTE INSTABILITY

7

continue to deform the Γω contour. If we can continue to deform the Γω contour onto the real axis then (2.11) still holds, and we have also determined the correct position of the Γk contour. Whether a root lies above or below the indented Γk contour determines whether it appears downstream or upstream (respectively). But if we cannot keep indenting the Γk contour then we cannot continue to deform the Γω contour and the asymptotic evaluation (2.11) is not valid. The only thing that can stop us indenting the Γk contour is if two sheets of roots, which start on different sides of the Γk contour when Im ω is large, merge at some value of ω. If this merging happens at (ω, k) = (ωp , kp ), the response will be dominated by a term of the form t−1/2 eıkp x−ıωp t . Since Im ωp > 0, this represents a response that is exponentially growing in time. This is absolute instability. Note that ωp and kp are independent of the form of the forcing and are solely dependent on the properties of the dispersion relation (2.8). Any causal forcing will excite this behaviour if it exists. It is simplest to think of the method this way, with a mapping determining roots in the k plane as functions of ω, and this is how this method was originally developed (Briggs, 1964). For technical reasons it is sometimes useful to determine roots in the ω plane as functions of k. A variant on the Briggs– Bers method which does this was given by Kupfer et al. (1987), but it is more complicated and is only tangential to this discussion.

2.2 Absolute instability for the membrane We now want to apply the Briggs–Bers method to the dispersion relation (2.8). It turns out that this dispersion relation has, in general, four roots — three with Re k > 0 and one with Re k < 0. (Actually, two of the three roots in the right half plane vanish when Re ω gets too large, but that is somewhat irrelevant.) The root in the left half plane, k1,− , comes from the lower half of the k plane as Im ω decreases to zero, does not interact with anything and all of the interesting behaviour happens for the three roots in the right half-plane. One of these roots, k2,− comes from the lower half of the k plane. The other two, k1,+ and k2,+ , come from the upper half of the k plane. The first thing we want to do is to find a criterion for absolute instability; and given the above results about the roots, we see that absolute instability can only occur as a result of a pinch between k2,− and either of k1,+ or k2,+ . It turns out that for a pinch with Im ωp > 0 we need to merge the roots k2,− and k2,+ , which happens for U greater than some critical value Uc . It is surprisingly easy to find a value for Uc , the speed at which we get marginal absolute instability with Im ωp = 0. When Im ω = 0, and Im k1,+ 6= 0, the two roots k1,+ and k2,+ form a complex conjugate pair, which at U = Uc happens to coincide with k2,− as Im k1,+ → 0. This gives p √us a triple root, with D, Dk and Dkk all equal to zero, and we find Uc = 6 3 − 9 ≈ 1.1800. The system is absolutely unstable for flows faster than this.

CHAPTER 2. THE BASIC PROBLEM

8

Figure 2.1: The closed contour of integration for x > 0

2.3 Other effects of causality What effects are possible if we are not in the absolutely unstable r´egime? In this case we only have to evaluate (2.11) and since the time dependence factors out we just need to consider Z 1 eıkx G(x, ω0 ) ≡ dk, (2.12) 2π Γk D(ω0 , k) where Γk is the indented contour as given by the Briggs–Bers method. The evaluation of (2.12) depends on the sign of x. For x ≥ 0 we close the k contour in the upper half plane, as in figure 2.1. Thus Z X eıkx 1 eıkx G(x, ω0 ) = ı − dk (2.13) Dk (ω0 , k) 2π down and up on ubc D(ω0 , k) k

for x > 0. The sum is over the modes which lie above the Γk contour and the remaining integral comes from the branch cut along the positive imaginary axis created by (k2 )1/2 . For x ≤ 0 we get Z X eıkx 1 eıkx G(x, ω0 ) = −ı − dk, (2.14) Dk (ω0 , k) 2π up and down on lbc D(ω0 , k) k

where the sum is over the modes which lie below the Γk contour and, again, there is a branch cut contribution. We see that, as mentioned in §2.1, the Briggs–Bers method determines which wave modes appear upstream and which appear downstream. It is the only rigorous way of doing so in problems with multiple energy sources. It also provides a way of dealing with a wave with complex group velocity, which is rather hard to interpret physically. Much of the behaviour of the system can be understood by looking at the neutral stability curve, which is qualitatively different for U < 1 and U > 1. For U < 1 there are essentially two regions, one with three neutral waves in the right half plane and one with only one neutral wave in the right half plane.

2.4. VARIANTS AND EXTENSIONS

9

When U > 1 an extra region appears, with only one neutral wave in the right half plane. This is shown in figures 2.2 and 2.3. Figure 2.4 is a close-up of the ‘nose’ in figure 2.3. We will mainly deal with the case 1 < U < Uc , so it will be useful to describe the response in a little detail. We start with 0 < ω0 < Ωl (U), where Ωl (U) is as defined graphically in figure 2.3. In this region there are two neutral waves and a complex conjugate pair. The neutral waves are observed upstream of the driver and the complex conjugate pair downstream. Note that this implies that the response is exponentially growing downstream — this is the region of convective instability. Now take the range Ωl (U) < ω0 < Ω(U). Here we have four neutral waves; the complex conjugate pair has moved onto the real axis, but this pair of waves are still observed downstream. Observe that near Ωl (U) there is a region where the mode k1,+ has negative group velocity, although we found that this mode always occurs downstream of the driver. Had we naively used the radiation condition to locate this mode, we would have put it upstream of the driver. We find, incidentally, that in this range the modes k1,+ , k2,+ and k2,− have negative activation energy, which is defined as the net energy required from an external agency to create a given wave mode from rest. An expression can be found for the activation energy of a neutral wave mode EA = −ω

∂D 2 |A| , ∂ω

(2.15)

where |A| is the wave amplitude (as in Cairns, 1979, although note that this paper uses a different definition of D which removes the minus sign in (2.15)). If a wave mode has negative energy, the excitation of this mode lowers the energy of the system; such waves are destabilised by membrane damping. In this range the modes k2,+ and k2,− transport energy towards the driver. The mode k1,+ transports energy away from the driver; it has negative group velocity and negative energy near the nose of the dispersion relation, but as ω is increased its group velocity and wave energy change sign and this mode becomes an ordinary positive energy wave. These unusual energy propagation results are found because there are two sources of energy in the system; the driver at x = 0 and the mean flow above the plate. When energy flows towards the driver, this simply means that the driver is providing a way of extracting energy from the mean flow. For ω > Ω(U) the modes k2,+ and k2,− form a complex conjugate pair, with Im k2,+ > 0. This region is entirely conventional. As ω increases further, the modes k2,+ and k2,− collide with the branch cut on Re k = 0 and vanish. This, as mentioned before, is somewhat irrelevant.

2.4 Variants and extensions The elastic plate (change Tηxx to −βηxxxx in (2.1)) behaves similarly. The dispersion relation becomes a quintic, which makes the solution somewhat harder, but it turns out that the qualitative behaviour of the plate system is roughly the same as that of the membrane system (Brazier-Smith & Scott, 1984).

CHAPTER 2. THE BASIC PROBLEM

10

0.9

0.8

0.7

0.6

ω0

0.5

0.4

0.3

0.2

Ω(U)

0.1

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

k

Figure 2.2: Neutral stability curve: U = 0.9

1.4

1.2

1

0.8

ω0 0.6

0.4

Ω(U) Ω l (U)

0.2

0 0

0.2

0.4

0.6

0.8

1

k

Figure 2.3: Neutral stability curve: U = 1.1

1.2

2.5. THE ELASTIC CHANNEL AND PIPE

11

0.24

0.23

0.22

ω 0 0.21

0.2

0.19

0.18 0.2

0.21

0.22

0.23

0.24

0.25

k

Figure 2.4: Close-up of neutral stability curve: U = 1.1

A number of other variants on this problem have been studied. The stability of flow through an elastic cylinder was considered by Peake (1997). The notable feature of this problem is that the hoop stress in the cylinder introduces a spring stiffness that pushes the absolute instability boundary to unphysical parameter r´egimes. The stability of Blasius flow over an elastic plate was considered by Lingwood & Peake (1999). For a thin shear layer, absolute instability only occurs at unphysical flow speeds, but as the shear layer thickens, absolute instability can be found at increasingly physically realisable speeds. Outside the absolutely unstable parameter r´egimes the behaviour with a realistic boundary layer was similar to the behaviour with a flat flow profile. Finally, Green & Crighton (2000) have examined the effects of point, rather than line, forcing. Absolute instability and the anomalous modes are still found.

2.5 The elastic channel and pipe Another variant on this problem is that of flow through an elastic channel. This problem has been previously examined by de Langre & Ouvrard (1999) using a lumped model of the fluid behaviour which is equivalent to plug flow in the limit of a thin channel, and by de Langre (2000), who studied a flat mean flow profile and potential flow. The next step in this (almost classical) problem is to consider Poiseuille flow in an elastic channel. Consider a channel of mean thickness 2H with flexible walls at y = H + ηT and y = −H + ηB through which flows a mean flow U(y) ≡ U∗max (1 − y2 /H2 ).

CHAPTER 2. THE BASIC PROBLEM

12

Linearise the Navier–Stokes equations about this mean flow to get ut + Uux + Uy v = −px /ρ + ν(uxx + uyy ) vt + Uvx = −py /ρ + ν(vxx + vyy ) ux + vy = 0.

(2.16) (2.17) (2.18)

Supposing that the channel walls are elastic plates with bending stiffness β and mass per unit length m gives the two equations βηTxxxx + mηTtt = p(y = H) βηB xxxx

+

mηB tt

(2.19)

= −p(y = −H).

(2.20)

The viscous normal stresses in the above actually evaluate to zero; if the fluid is viscous then we can use the tangential velocity condition and continuity; if the fluid is inviscid then the viscous stresses are by definition zero. Note that in using this plate model we neglect longitudinal plate displacements and tangential stresses. Finally we need the linearised boundary conditions v(y = H) = ηTt

(2.21)

ηB t

v(y = −H) = u(y = ±H) = 0.

(2.22) (2.23)

Now nondimensionalise with a vertical lengthscale H, a horizontal lengthscale L ≡ (ρ/m) and a timescale T ≡ (m/β)1/2 (m/ρ)2 . Let u = L/Tu 0 , U = L/TU 0 , v = H/Tv 0 , p = ρL2 /T 2 Hp 0 and η = Hη 0 . Drop the primes and Fourier transform in x and t to obtain the equations 2 2 −ıωu + ıkUu + Uy v = −ık²p + Re−1 fl (uyy − ² k u)

Re−1 fl (vyy

−ıωv + ıkUv = −py /² + ıku + vy = 0

− ² k v)

v(y = 1) = −ıωηT B

v(y = −1) = −ıωη u(y = ±1) = 0

k4 ηT − ω2 ηT = p(y = 1) 4 B

2 B

2 2

k η − ω η = −p(y = −1),

(2.24) (2.25) (2.26) (2.27) (2.28) (2.29) (2.30) (2.31)

where U(y) = Umax (1 − y2 ), Refl ≡ H2 /(νT ), ² ≡ H/L and Umax ≡ U∗max T/L. We can recombine these three parameters to get the fluid dynamic Reynolds number, Re ≡ Umax Refl /². The equations (2.24–2.31) only have nonzero solutions if ω and k satisfy a dispersion relation 0 = D(ω, k; Umax , ², Refl ) ≡ DS (ω, k; Umax , ², Refl )DV (ω, k; Umax , ², Refl ), (2.32) where DS is the dispersion relation for sinuous modes (ηT = ηB ) and DV the dispersion relation for varicose modes (ηT = −ηB ). We wish to trace the branches of the dispersion relations DS and DV .

2.5. THE ELASTIC CHANNEL AND PIPE

13

2.5.1 Inviscid perturbation Taking values for 5cm thick aluminium in water, we find that Refl À 1. We observe that Crighton & Oswell (1991) found instability in their problem at all nonzero flow speeds, so we consider the limit Refl → ∞, Umax ¿ 1 with 1/2 ² = o(Refl ). In taking this limit we lose the no-slip conditions (2.29), and derive Rayleigh’s equation vyy + (Uyy k/(ω − Uk) − k2 ²2 )v = 0,

(2.33)

subject to the boundary conditions (2.27, 2.28, 2.30, 2.31), where the pressure p is given by k²p = ı (ωvy /k + Uy v − Uvy ) . (2.34) This system can now be solved numerically, either by combining shooting with a root finder or by combining Chebyshev collocation with standard matrix eigenvalue software (Boyd, 1999; Anderson et al., 1999). Some care must be taken with critical layers (when ω = U(yc )k), these are avoided by continuing the integration into the complex plane (Lingwood & Peake, 1999). Some results are shown in figures 2.5 and 2.6. In figure 2.5 we observe convective instability; a wave mode that crosses the real axis as Im ω is decreased to zero. In figure 2.6 we observe anomalous propagation; a wave mode that crosses, and then circles back to, the real axis as Im ω is decreased to zero. This behaviour is very similar to that found by Brazier-Smith & Scott (1984) and Crighton & Oswell (1991) in their inviscid irrotational system. A full parameter search and a search for absolute instability has not yet been completed.

CHAPTER 2. THE BASIC PROBLEM

14

0.6

x

0.4

x 0.2 x

x

0 x

x

-0.2

-0.4

-0.6

x

-0.8 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

Figure 2.5: The paths of roots in the k plane as Im ω is decreased to zero, U = 7 × 10−2 , sinuous disturbance, ² = 1 and Re ω = 1.5 × 10−3 . × marks the start of the paths for large Im ω. 0.6

0.4

x x

0.2 x

x

0 x

x

-0.2

-0.4

-0.6

-0.8 -0.6

x

-0.4

-0.2

0

0.2

0.4

0.6

Figure 2.6: The paths of roots in the k plane as Im ω is decreased to zero, U = 7 × 10−2 , sinuous disturbance, ² = 0.95 and Re ω = 1.5 × 10−3 . × marks the start of the paths for large Im ω.

Chapter 3

Ribbed elastic structures in a mean flow In §2 we developed the Green’s function for a fluid-loaded membrane immersed in a mean flow and driven on one line. One application of this result is to study the properties of a ribbed membrane — an elastic membrane immersed in a mean flow and supported at certain ribs. This is of interest to naval engineers as a model of a ship’s hull, which has constraints applied to it at points along its length by its coupling to the superstructure. Actually, a ship’s hull is possibly better modelled as a plate, but this introduces complications that are avoided by considering the membrane. Both problems are discussed here. There are two different problems that we can consider in the context of the ribbed membrane or plate. The first, which is not discussed here, is wave scattering (Lawrie, 1989; Llewelyn Smith & Craster, 1999). The second, which is discussed here, is the problem of energy transmission. We suppose that one of the ribs in the array is oscillated, and ask whether the disturbance so generated can propagate along the array, or whether it is trapped near the driver. Previous work on this problem has only considered the case of static fluid, and in this case it is found that for periodic rib arrays there is a pass and stop band structure — certain frequencies can pass down the array without attenuation but others are trapped near the driver (Crighton, 1984; Spivack, 1991). Disordered arrays slightly delocalize the stop bands, but destroy the delicate phase matching needed for a pass band, resulting in exponential localization of the system’s response (Anderson, 1958; Hodges & Woodhouse, 1989; Sobnack & Crighton, 1994; Spivack & Barbone, 1994). This Anderson localization is relatively universal in wavebearing systems, and has been repeatedly observed both in acoustic experiments (Photiadis & Houston, 1999) and electromagnetic experiments (Dembowski et al., 1999). When the hydrodynamic component of the Green’s function is included, which provides a long range coupling term, the pass bands of the ordered system are not significantly altered (Cooper & Crighton, 1998b). However, the exponential localization found in the stop bands of the ordered system is destroyed, and the decay is only algebraic (Cooper & Crighton, 1998a). Similarly, 15

CHAPTER 3. RIBBED ELASTIC STRUCTURES

16

in disordered arrays, we find only an algebraic localization when the hydrodynamic part of the Green’s function is included (Cooper & Crighton, 1999; Baesens & MacKay, 1999). The inclusion of a mean flow is expected to provide interesting new behaviour, partly because it will break the symmetry of the zero mean flow problem, partly because it is a system which supports more free waves and partly because some of these free waves are actually instabilities. Any or all of these properties could lead to interesting results. This ribbed membrane problem, besides having its own intrinsic interest, is also very closely related to some areas of dynamical systems and solid state physics (Thouless, 1974; Flach & Willis, 1998). These links will be briefly discussed as and when they become relevant. In the rest of §3 we study the response of both periodically and aperiodically ribbed elastic structures under a mean flow. We will derive the governing equations of the system in §3.1, simplify them in §3.2 and discuss some results in §3.3, §3.4 and §3.5.

3.1 Governing equations Suppose that we have a membrane supported by ribs at points x1 , x2 , . . . , xN and driven on the Dth rib by a force FD . We may safely assume that x1 < x2 < · · · < xN . For simplicity, we will suppose that the non-driven ribs have infinite mechanical impedance and so must have a zero velocity. We want to find the forces F1 , F2 , . . . , FN on the rib array. We insist that we are not in the absolutely unstable r´egime, so that the time dependence of the Green’s function can be factored out and the difference between velocity and displacement is irrelevant. Now, let G(x) be the Green’s function of the membrane. By linearity, the plate displacement η is found to be η(x) =

N X

Fk G(x − xk ),

(3.1)

k=1

which gives us the set of N − 1 linear equations 0=

N X

Fk G(xj − xk )

for j 6= D.

(3.2)

k=1

We can view either the force applied to the Dth rib or the displacement of the Dth rib as an unknown. If we set the displacement of the Dth rib to be 1, we get the additional equation 1=

N X

Fk G(xD − xk ),

(3.3)

k=1

which just sets the scale of the system. This gives us the matrix equation G · F = η, where Gij = G(xi − xj ), which is a set of N linear equations in N unknowns and we now have a system we can solve to find the forces Fi .

3.2. REDUCTION OF GREEN’S FUNCTION MATRIX

17

It is in theory possible to numerically invert the matrix G by methods such as LU or QR factorisation, but unfortunately the matrix G is ill-conditioned — the neutral surface waves maintain a constant amplitude and do not decay. Thus the matrix G is far from diagonally dominant and one cannot make N particularly large before catastrophic loss of precision occurs in standard factorisation routines. The problem is even worse in the convectively unstable r´egime. We will see that it is possible to recombine these equations to give a numerically well-conditioned matrix. This rearrangement of the equations will also give us a means of understanding the qualitative form of the response of the system. Note that the non-driven ribs need not have an infinite impedance. This assumption can be trivially relaxed, as done by Cook (1998). One other point to note is that the Green’s function grows exponentially downstream at sufficiently low frequencies when U > 1. Thus it is not immediately obvious that the infinite rib array problem of Crighton (1984) is well-defined, and we must consider the finite array problem. Another point to note is that the infinite array problem can only be thought of as a model of the startup of forcing, when reflections from the ends of the array have not returned to the driver. We only have the Green’s function in the long-time limit, so we can only consider a fully developed solution and thus, necessarily, a finite array.

3.1.1 The Green’s function The Green’s function (from (2.12) of §2.3) is ± P R ıki,+ x 1 eıkx ı i Dke(ω,ki,+ ) − 2π dk x>0 ubc D(ω,k) R G(x) = P eıki,− x 1 eıkx −ı i Dk (ω,ki,− ) − 2π lbc D(ω,k) dk, x < 0,

(3.4)

where, as in §2.3, the integrals are taken along the appropriate branch cuts. We see that we can write the Green’s function matrix G as G = Gsw + Gbc ,

(3.5)

where Gsw contains the surface wave parts of G and Gbc almost all of the branch line contributions. As is usual (Crighton, 1984), we include the line admittances G(0) in Gsw , this means that the diagonal of Gbc is zero. We also observe that it is trivial to evaluate the Green’s function numerically. Previous studies of this problem have tended to use asymptotic results, in order to produce equations that can be manipulated analytically. In this problem there is little hope of doing any analytical manipulation and we will use numerically computed values for the Green’s function in all of the calculations in §3.

3.2 Reduction of Green’s function matrix The only true long range coupling in the Green’s function is that provided by the branch line component. This component of the Green’s function is weak, and is sufficiently strongly decaying to be trivially invertible. The surface wave components, although apparently a long range coupling, can be

CHAPTER 3. RIBBED ELASTIC STRUCTURES

18

recast as a purely local term. This derivation is not conceptually difficult, although it does require a little thought. ˆ ≡ S · G, where S is a pentadiagonal matrix with We consider the matrix G Si(i−2) = ai , Si(i−1) = bi , Sii = 1, Si(i+1) = ci and Si(i+2) = di , and observe ˆ sw pentadiagonal if ai , bi , ci and di satisfy the that we can make the bulk of G equations     ai −1 bi  −1    Mi ·  (3.6)  ci  = −1 di −1 for i = 3 . . . N − 2, where  ık1,+ (xi−2 −xi ) e eık2,+ (xi−2 −xi ) Mi ≡  eık1,− (xi−2 −xi ) eık2,− (xi−2 −xi )

eık1,+ (xi−1 −xi ) eık2,+ (xi−1 −xi ) eık1,− (xi−1 −xi ) eık2,− (xi−1 −xi )

eık1,+ (xi+1 −xi ) eık2,+ (xi+1 −xi ) eık1,− (xi+1 −xi ) eık2,− (xi+1 −xi )

 eık1,+ (xi+2 −xi ) eık2,+ (xi+2 −xi )  . eık1,− (xi+2 −xi )  eık2,− (xi+2 −xi ) (3.7)

The left end of the array gives the equations c1 eık1,− (x2 −x1 ) + d1 eık1,− (x3 −x1 ) = −1 c1 eık2,− (x2 −x1 ) + d1 eık2,− (x3 −x1 ) = −1 b2 eık1,− (x1 −x2 ) + c2 eık1,− (x3 −x2 ) + d2 eık1,− (x4 −x2 ) = −1

(3.8)

b2 eık2,− (x1 −x2 ) + c2 eık2,− (x3 −x2 ) + d2 eık2,− (x4 −x2 ) = −1 and the right end of the array gives aN−1 eık1,+ (xN−3 −xN−1 ) + bN−1 eık1,+ (xN−2 −xN−1 ) + cN−1 eık1,+ (xN −xN−1 ) = −1 aN−1 eık2,+ (xN−3 −xN−1 ) + bN−1 eık2,+ (xN−2 −xN−1 ) + cN−1 eık2,+ (xN −xN−1 ) = −1 aN eık1,+ (xN−2 −xN ) + bN eık1,+ (xN−1 −xN ) = −1 aN eık2,+ (xN−2 −xN ) + bN eık2,+ (xN−1 −xN ) = −1. (3.9) Note that b2 , c2 and d2 , and aN−1 , bN−1 and dN−1 are indeterminate. Although it might be useful to set d2 = aN−1 = 0, we observe that when the ribs are periodically spaced this choice will set c1 = b2 , d1 = c2 , bN−1 = aN and cN−1 = bN . Therefore the matrix S has zero determinant in this case and we ˆ It proves convenient to set d2 = aN−1 = 1. are not able to invert the matrix G. ˆ If we temporarily ignore Gbc , we see that this transformation produces the ˆ sw · F = S · η, which is pentadiagonal system of equations G αi,−2 Fi−2 + αi,−1 Fi−1 + αi,0 Fi + αi,1 Fi+1 + αi,2 Fi+2 = ai ηi−2 + bi ηi−1 + ηi + ci ηi+1 + di ηi+2

(3.10)

for i = 3 . . . N − 2, for some values αi,j which it is not useful to write explicitly. There are also four other equations from the endpoints of the rib array. It is clear that obvious variants of this reduction scheme are applicable to other similar wavebearing systems and will allow the reduction of a general wavelike (or exponential) coupling to a purely local effect. Unlike other such schemes (Spivack & Barbone, 1994; Photiadis, 1992), this trick is trivial to derive, understand and generalise.

3.2. REDUCTION OF GREEN’S FUNCTION MATRIX

19

3.2.1 Transfer matrix form We note that away from the driven rib the recurrence (3.10) can be cast in a transfer matrix form by writing     Fi+1 Fi  Fi      = Ti · Fi−1  , (3.11) Fi−1  Fi−2  Fi−2 Fi−3 where the ith transfer matrix Ti governs the propagation from the ith rib to the (i + 1)st rib and is simply given as  αi−1,1 − αi−1,2  1 Ti =   0 0

i−1,0 −α αi−1,2 0 1 0

− ααi−1,−1 i−1,2 0 0 1

 − ααi−1,−2 i−1,2  0 .  0 0

(3.12)

This formulation will prove to be useful later.

3.2.2 Reconstruction of solution, energy flux and other useful tools In this section we briefly derive some useful results that will be helpful later. Spivack & Barbone (1994) were able to produce some strong results in their system, and some of these will be seen to carry over to this problem. However, their brute force methods are much less useful here and we proceed from more general principles. Firstly, we note that we can reconstruct the displacement of the membrane between two ribs, given the forces on the ribs and neglecting the branch line component of the Green’s function. Suppose that xi−1 < x < xi , and that neither the (i − 2)th , (i − 1)th , ith or (i + 1)th ribs are driven. Then, by applying the ideas that led to equations (3.6,3.7), we find that η(x) can be written in the form T η(x) = B(x) · (Fi+1 Fi Fi−1 Fi−2 ) , (3.13) for some vector B(x) for which an expression can be found. Next we observe that, between two ribs, the pressure on the membrane is a linear functional of η(x), ¡ ¢ p(y = 0) = ω2 + ∂2x η, (3.14) and similarly we note that we can find the fluid velocity, φy (y = 0) = (−ıω + U∂x ) η.

(3.15)

Lastly, under the assumption that the branch line component of the solution is not present, we observe that we can write φ as a linear functional of η, ηx , ηxx and ηxxx , and we see that η, φy , φ and p are all linear in the state T vector (Fi+1 Fi Fi−1 Fi−2 ) .

20

CHAPTER 3. RIBBED ELASTIC STRUCTURES

To make all of these observations useful, we next note that the system allows us to define an energy flux (Crighton & Oswell, 1991), Z1 J(x, t) = −ηx ηt + φt φx dy + Uηφt . (3.16) 0

Since we work with time-harmonic solutions, we define an average energy flux µ ¶ Z1 J(x) = Im −η∗ ηx + Uφ∗ η +

φ∗ φx dy .

(3.17)

0

At the midpoint of a bay we can evaluate all of the functions in the mean enT ergy flux (3.17) as linear functionals of (Fi+1 Fi Fi−1 Fi−2 ) , and so we observe that conservation of energy implies the existance of a set of matrices Σ(i) such that Σ(i)† = −Σ(i) and  †   Fi+1 Fi+1  Fi    (i)  Fi    (3.18) Fi−1  · Σ · Fi−1  Fi−2 Fi−2 is independent of i away from the driven rib, which is equivalent to the statement T†i · Σ(i) · Ti = Σ(i−1) . (3.19) When the system is periodic, the matrices Σ(i) are constant and so the transfer matrix T preserves the alternating form (x, Σ · y). It is of use to investigate this special case. The main point to note is that if e is an eigenvector of T with eigenvalue λ, then Σ · e is an eigenvector of T∗ with eigenvalue 1/λ. Thus if λ is an eigenvalue of T so is 1/λ∗ , and if |λ| 6= 1 this construction provides a pair of eigenvalues of T.

3.2.3

Branch line contribution

Observing that left multiplication by S just consists of row operations, we see that S · Gbc has the same decay rate away from the central diagonal as Gbc . ˆ bc is accurate, if not cheap. Writing G ˆ bc = G ˆ (1) + Thus numerical inversion of G bc ˆ bc and G ˆ (2) is the remainder, we ˆ (1) is the pentadiagonal part of G ˆ (2) , where G G bc bc bc see that we have derived the set of linear equations ³ ´ ˆ sw + G ˆ (1) · F + G ˆ (2) · F = S · η. G (3.20) bc bc We now have two choices. If the number of ribs is not too large, we can diˆ using LU factorisation. This direct inversion is wellrectly invert the matrix G conditioned because we have localized the surface wave component; the surface wave inversion acts as a preconditioner for the whole system. For arrays of up to a few hundred ribs there is no need to resort to iterative methods and ˆ is the best way to proceed. For longer direct LU factorisation of the matrix G arrays an iterative scheme may be necessary, but an excessively complicated scheme should not be needed and the iteration ³ ´ ˆ sw + G ˆ (1) · F(n+1) = S · η − G ˆ (2) · F(n) , with F(0) = 0 G (3.21) bc bc

3.3. PERIODICALLY RIBBED MEMBRANES

21

(2)

(n) ˆ should be sufficient, since we expect G to be a small correction. This bc · F iteration also allows us to interpret results found by LU factorisation; we will ˆ sw ; the most see later that the dominant coupling is actually contained in G important part of the branch line contribution is contained in the line admittance G(0). This idea is equivalent to the manipulation of Cooper & Crighton (1998a), who took the iteration (3.21) as far as F(2) . Numerically, of course, there is no reason to do this and we might as well iterate (3.21) to convergence. This iterative scheme was found to be unnecessary for the sizes of problems studied here.

3.3 Periodically ribbed membranes As derived above, our simplifying transformation is applicable to both periodic and disordered rib arrays, but in this section we only consider periodic structures. Suppose that the inter-rib spacing xi+1 − xi is a constant; say h. In this case the coefficients αi,j and the transfer matrices Ti are found to be independent of i (and the i suffices on both these quantities will be omitted throughout the rest of §3.3). Similarly, the coefficients ai , bi , ci and di are constant for i = 3 . . . N − 2 and will be written a, b, c, d. ˆ using standard routines and results are shown We can invert the matrix G in figure 3.6 (on page 29). We observe a very distinctive pass and stop band structure. In the pass bands the response of the system is O(1) on the whole array. In the stop bands we observe exponential decay of the response away from the driven rib. Note also the massive asymmetry. This asymmetry, although due to the the mean flow, is not simply controlled by the mean flow. It is easy to find parameter values at which the array response upstream of the driver is large and the array response downstream of the driver is small. ˆ sw · F = η; This structure can be explained by the system of equations G the branch line parts of the Green’s function matrix just produce a change in the quantitative response; the form of the response is controlled by the recurrence (3.10) and the endpoints of the rib array. This can be seen clearly in figure 3.2. We find that of the four eigenvalues λ1 , . . . , λ4 of the transfer matrix T, |λ3 | À 1 and |λ1 | ≈ 1, with |λ3 λ4 | = |λ1 λ2 | = 1. Note that in §3.2.2 we predicted that eigenvalues off the unit circle come in pairs (λ, µ) such that λµ∗ = 1. This relation was satisfied with an error of 10−10 at most. In a pass band the eigenvalues λ1 and λ2 lie on the unit circle and move around it as the frequency is varied. When the system leaves a pass band, the two eigenvalues collide and move off the unit circle. When the system reenters a pass band, they collide again, attach themselves to the unit circle and begin to move around it. This can be observed in figure 3.1, which is a plot of |λ1 | and |λ2 | as functions of frequency. Note how closely figures 3.6 and 3.1 coincide; almost everywhere we see a localized response where the eigenvalues are off the unit circle and a delocalized response when they are on the unit circle. It is not difficult to produce an argument (Borland, 1963) which suggests that eigenmodes whose eigenvalues are not on the unit circle are localized

CHAPTER 3. RIBBED ELASTIC STRUCTURES

22 3.5

3

2.5

2

|λ| 1.5

1

0.5

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

ω

Figure 3.1: Plot of moduli of eigenvalues of transfer matrix of periodic system, with U = 1.1 and h = 10. near the driven rib and endpoints of the array. Thus only eigenmodes whose eigenvalues lie on the unit circle can pass along the whole of the array. This does not necessarily mean that they will do so with a large amplitude — we observe “pass bands” in which the amplitude in the bulk of the array is less than 1% of the maximum amplitude on the array. Compare, for instance, the low frequency regions of figures 3.6 and 3.1, where figure 3.1 shows that there are two eigenvalues on the unit circle, but figure 3.6 shows a localized response. This low amplitude response is simply because the boundary conditions and driver conditions can be satisfied by the localized modes without resort to the extended modes.

3.4 Disordered rib arrays When a periodic wavebearing system becomes disordered, we expect to observe the phenomenon of Anderson localization (Anderson, 1958), a universal behaviour which is to be expected when periodic structures become disordered. Generically, one expects the effect of Anderson localization to confine any response of the system to the region near the driver. In our case, this means that we expect the pass band response to become localized about the driven rib. For definiteness, we introduce disorder by displacing each rib from some mean position by a random amount, xi = h(i + σUi ),

(3.22)

where the random variables Ui are independent, each with uniform distribution on [−1, 1], and σ is a disorder parameter.

3.4. DISORDERED RIB ARRAYS

23

In §3.4 we neglect the branch line part of the Green’s function. This is not strictly necessary, but since the branch line part was found to have so little effect in §3.3 we make this simplification, which clarifies the problem and greatly speeds up the numerics. It will later be necessary to include the branch line component, but that is left as further work. ˆ sw · F = S · η and It is trivial to numerically solve the linear equations G some results are shown in figure 3.7 (on page 29). Figure 3.2 compares the response of ordered and disordered arrays, at a frequency which is in a pass band of the ordered system. Interestingly, we observe a massive delocalization in regions of parameter space where, although the recurrence (3.10) allows a pass band response in the periodic system, the boundary conditions at the ends of the rib array and the driver conditions prevent any extended modes from appearing. This can be seen clearly in figure 3.3. A physical explanation is simple; when the system becomes aperiodic the centre of the array cannot feel these end effects, wants to propagate and can now, to some extent do so. Although we expect the system response to become localized when the disorder is sufficiently large, it appears that there is a window in which disorder obscures the end effects but is not strong enough to kill off propagation altogether. This, then, is the “sensitivity to boundary conditions” used by Thouless (1974) as the characteristic of extended modes. This particularly dramatic manifestation of this sensitivity does not seem to have been previously observed in studies of similar systems (Anderson, 1958; Thouless, 1978; Flach & Willis, 1998). It is useful in this context to introduce the Lyapunov exponents of the system (Kottos et al., 1999), which give some measure of its localization behaviour, and the inverse of the smallest Lyapunov exponent (in absolute value) gives some kind of localization length of the system. If we introduce an ‘endto-end’ transfer matrix Ttot =

N−1 Y

Ti ,

(3.23)

i=4

with singular values σ1 , . . . , σ4 , we can define the Lyapunov exponents as ­ ® li = (N − 4)−1 log σi ,

(3.24)

where hfi denotes an average over realizations. There are various theorems which ensure that such an average exists. It is generally found that the addition of some kind of disorder to such a system causes the Lyapunov exponents to repel each other. This can be clearly seen in figures 3.4 and 3.5, which show the Lyapunov exponents as functions of frequency for the same parameter values as figures 3.6 and 3.7 respectively. Nevertheless, even with this repulsive effect we see that as we turn on the disorder there are frequency intervals in which the smallest Lyapunov exponents remain very small. This means that the system is still able to propagate signals even when disordered. In these r´egimes the effect of disorder is simply to obscure the end effects, which, in this system, can allow a signal to propagate.

CHAPTER 3. RIBBED ELASTIC STRUCTURES

24

1 ordered, no branch line ordered, branch line disordered, no branch line

0.9 0.8 0.7 0.6

force 0.5 0.4 0.3 0.2 0.1 0 0

10

20

30

40

50

60

70

80

90

100

rib number

Figure 3.2: Comparison of forces along a rib array for ordered and disordered systems, in pass band of ordered system. 1 ordered system disordered system 0.9 0.8 0.7 0.6

force 0.5 0.4 0.3 0.2 0.1 0 0

10

20

30

40

50

60

70

80

rib number

Figure 3.3: Delocalization caused by disorder.

90

100

3.4. DISORDERED RIB ARRAYS

25

4

3

2

1

0

-1

-2

-3

-4 0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

ω

Figure 3.4: Plot of finite length Lyapunov exponents for system of 100 ribs with U = 1.1, h = 10 and σ = 0. 4

3

2

1

0

-1

-2

-3

-4 0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

ω

Figure 3.5: Plot of finite length Lyapunov exponents for system of 100 ribs, with U = 1.1, h = 10 and σ = 5%.

CHAPTER 3. RIBBED ELASTIC STRUCTURES

26

3.5 Ribbed plates Since the Green’s function of the elastic plate under mean flow is qualitatively similar to that of the elastic membrane, and the plate problem is more physically plausible anyway, we would also like to apply these ideas to the problem of a ribbed fluid-loaded plate. The complication is that we must solve for both a force and a torque on each rib. Previous studies of the plate problem have either neglected the torque (Eatwell & Willis, 1982; Photiadis, 1992) or have used a dubious model of the fluid loading (Rumerman, 1975), or both, so it is of interest to do this problem taking both of these factors into account. Assuming that the plate is clamped at each rib (Guo, 1993), we find X X η(x) = Fj Gp (x − xj ) + Tj Gp0 (x − xj ), (3.25) j

j

where Gp is the plate Green’s function, Fj is the force on the jth rib and Tj is the torque on the jth rib. This gives us the set of equations X X 0= Fj Gp (xi − xj ) + (3.26) Tj Gp0 (xi − xj ) i 6= D j

1=

X

j

Fj Gp (xD − xj ) +

j

0=

X

X

Tj Gp0 (xD − xj )

j

Fj Gp0 (xi − xj ) +

j

X

i=D

Tj Gp00 (xi − xj ),

(3.27) (3.28)

j

where Eqs. (3.26,3.27) come from the plate displacement on the rib array and Eq. (3.28) comes from the clamped rib condition ηx (x = xi ) = 0. We can rewrite Eqs. (3.26–3.28) as a set of linear equations of the form ¶ µ ¶ µ ¶ µ 0 0 η + Gbc F Gsw + Gbc Gsw (3.29) = · 0 0 00 00 ηx Gsw + Gbc Gsw + Gbc T Next, we observe that the surface wave part of the plate Green’s function Gp (x) is qualitatively similar to the membrane Green’s function, that all the derivatives of the plate Green’s function have waves of the same wavelength and that the matrix S is only a function of wavelength and rib position. This means that we can apply the simplifying transformation to the plate equations by left-multiplying Eq. (3.29) by µ ¶ S 0 Sp ≡ , (3.30) 0 S to obtain

µ ˆ sw + G ˆ bc G 0 ˆ ˆ0 Gsw + G bc

0 ˆ sw ˆ0 G +G bc 00 ˆ ˆ 00 Gsw + G bc

¶ µ ¶ µ ¶ F S·η · = , T S · ηx

(3.31)

ˆ G ˆ 0 and G ˆ 00 are all pentadiagonal. Neglecting the branch line conwhere G, tributions, we can rewrite the simplified plate equations (3.31) as a matrix of bandwidth 11, which can be cheaply and accurately solved using standard library routines.

Results for the force and the torque on a ribbed plate, neglecting the branch line component of the Green’s function, are shown in figures 3.8 and 3.9 respectively. The frequency range shown covers the entire range of convective instability and anomalous propagation of the forced plate. We observe both pass and stop band behaviour, as in the ribbed membrane, although the distinction is much less sharp.

3.6 Discussion The pass and stop band structures found when the fluid is taken to be static are shown to persist as the mean flow speed is increased from zero. We see strong asymmetries at low frequency, when the effect of the fluid loading is important. An interesting result is that the branch line part of the membrane Green’s function has little effect on the response of a periodic array at the flow speeds we are interested in. This work needs to be extended to the plate problem and to the disordered array. These ribbed elastic structure problems provide yet another context for work on discrete breathers and band random matrices (Flach & Willis, 1998; Kottos et al., 1998), in which one is not restricted to the comparatively simple couplings that have been previously considered in the context of solid state physics. This more general interaction allows us to find, in a real system, massive delocalization of the response upon the introduction of disorder. Finally, we observe that we have a complete description of the response of a ribbed elastic plate or membrane in the time-harmonic steady state, and that this work could be trivially extended to consider, say, a ribbed elastic cylinder (Peake, 1997; Photiadis & Houston, 1999). The only major remaining problem is to study the transient response to forcing, which would allow us to investigate the absolutely unstable parameter r´egime.

27

Figure 3.6: Plot of th

|Fn |(ω) kFk1 (ω)

with U = 1.1 and h = 10 for system of 100 ribs,

driven at the 50 . The vertical axis is frequency, with a range 10−3 to 1, and the horizontal axis is rib number. The mean flow is from left to right.

Figure 3.7: Plot of

|Fn |(ω) kFk1 (ω) th

with U = 1.1, h = 10 and σ = 5%, for system of 100

ribs, driven at the 50 . The vertical axis is frequency, with a range 10−3 to 1, and the horizontal axis is rib number. The mean flow is from left to right.

Figure 3.8: Plot of force on rib array with U = 0.05 and h = 10 for system of 100 ribs, driven at the 50th . The vertical axis is frequency, with a range 10−3 to 8 × 10−3 , and the horizontal axis is rib number. The mean flow is from left to right.

Figure 3.9: Plot of torque on rib array with U = 0.05 and h = 10 for system of 100 ribs, driven at the 50th . The vertical axis is frequency, with a range 10−3 to 8 × 10−3 , and the horizontal axis is rib number. The mean flow is from left to right.

32

Chapter 4

The finite baffle and other numerical problems Although we are able to make a great deal of analytic progress on the infinite plate problem, the infinite membrane problem and the rib array problem, there are other problems in this area where the primary tool of investigation is numerical simulation. For instance, in §2 we discussed the problem of an infinite plate forced on one line, and it is of interest to see if the infinite plate problem has any relation to the plate of finite length. This, of course, is not a new problem and has been studied before, but with conflicting results. Lucey (1998) and Lucey & Carpenter (1992) did a full numerical simulation of the system using a boundary integral method with low order spline interpolation to discretize the surface. They reported a major difference between the finite plate and infinite plate; that the finite plate was unstable at all flow speeds. They were also unable to find the negative energy waves found by Crighton & Oswell (1991) on the infinite plate. Wright (2000) discretised the surface as a sum of a small number of Fourier modes and derived evolution equations for the mode amplitudes. Whilst this may make the analysis simple, it is not a numerically efficient way of simulating this system1 He reports that the finite plate is unstable only when the infinite plate is absolutely unstable, a conflict which can perhaps be resolved by noting that Wright (2000) used a rather crude numerical scheme and considered plates much shorter than those considered by Lucey (1998). Lucey attributes his instability to negative energy effects. An alternative explanation is to think of a convectively unstable wave generated as the forcing is switched on. It travels along the plate from the driver, absorbing energy from the mean flow as it goes. On reaching the far end of the driver other wave modes must be added in to the solution to fix the boundary conditions. These take energy back to the driver, or to the other end of the plate, where the convective growth process starts again. 1 In Wright’s problem, the plate displacement has its zeroth and second derivatives zero at the ends of the baffle. The (appropriately scaled) sine functions satisfy this condition, but also have all of their even derivatives zero at the endpoints. These extra conditions are not necessarily satisfied by the plate displacement. This increases the error from the exponentially small error that we might expect with a spectral discretisation to the algebraically small error that we get from finite differences (Boyd, 1999).

33

U

-a

a F

Figure 4.1: The finite baffle In general there is destructive interference between all this growth, but at some resonant frequencies there is constructive interference and we get temporal growth. This can be examined by considering a toy problem with one upstream wave mode and one downstream wave mode, fitted to zero boundary conditions at the endpoints of a domain and with a jump in the first derivative at the origin. Such a problem can only be solved with a timeharmonic response if resonance does not occur; it is this resonance to harmonics generated at startup which causes instability. This is a more complicated version of a problem examined by Tobias et al. (1998), who considered the effect of convective instability on solutions of the Ginzburg–Landau equation in finite domains. They noted that if it is possible for energy generated in convective growth to return to the body of the domain then such a feedback loop is possible. If it is not possible for this energy to return then this loop cannot establish itself and convective instability of the infinite system has no effect on the behaviour of the finite system. So, we consider an elastic plate of length 2a, Young’s modulus E, mass per unit area m with the moment of inertia per unit width of its cross-section I, attached to a hard wall at either end. Above it flows an (incompressible, inviscid, irrotational) fluid with density ρ at speed U1 , as sketched in figure 4.1. The natural length and time scales are the fluid loading ones; the lengthscale is (m/ρ) and the timescale is m3 /(ρ5 EI)1/2 . The lengthscale is natural; it is the depth of fluid that interacts with the plate. The timescale comes from equating the bending moment and inertia of a slice of plate. Choosing these units sets m = ρ = EI = 1, whilst U1 7→ U1 m2 /(ρ3 EI)1/2 and a 7→ aρ/m. Without much effort we obtain the linearised equations ∇2 φ = 0 φ→0

y>0

as y → ∞

φy = 0 |x| > a, y = 0 φy = (∂t + U1 ∂x )η |x| < a, y = 0 p = −(∂t + U1 ∂x )φ ηxxxx + ηtt = −p(y = 0) + F(x, t) |x| < a η(x = ±a) = 0,

(4.1) (4.2) (4.3) (4.4) (4.5) (4.6) (4.7)

where φ is the perturbation velocity potential, η is the plate displacement and 34

p is the perturbation fluid pressure. Note that we need two more conditions, which are boundary conditions on η at x = ±a. There are two possible choices, ηx (x = ±a) = 0 or ηxx (x = ±a) = 0. The second choice allows us to expand the plate displacement as a sum of sine functions, but as mentioned earlier, this is not numerically useful. Since Leppington et al. (1984) found the ηx = 0 boundary condition more consistent with experiment for the acoustic scattering problem, we will use that boundary condition. Another reason for this boundary condition is purely technical; it lessens the severity of corner singularities at the ends of the plate — see Garrad & Carpenter (1982) for more details.

4.1 Numerical scheme The obvious method of solving equations (4.1–4.7) is to use a boundary element method, and that is indeed what we will use. In these methods we convert the volume differential equations into surface integrals by use of the Green’s function of the system. This approach is standard, and should need little explanation. The next thing to note is that although Laplace’s equation and the linearised Bernoulli equation (4.1,4.5) imply ∇2 p = 0, this relationship is not maintained in the nonlinear r´egime. Although it would be simpler to solve for the pressure, we must solve for the velocity potential if we want the numerical scheme to generalise into the nonlinear r´egime. Actually, in the linearised r´egime it makes little difference what we do, since we end up producing the same integro-differential equation. But the derivation of this equation parallels the numerical scheme which must be adopted in the nonlinear r´egime. Let G(x, y, ξ, ζ) be the solution of ∇2x,y G = δ(x − ξ)δ(y − ζ)

(4.8)

Gy (x, 0, ξ, ζ) = 0 ∇x,y G → 0 as y → ∞,

(4.9) (4.10)

where the subscripts on the differential operators indicate the variables on which they work. This gives ¡ ¢ ¡ ¢¢ 1 ¡ log (x − ξ)2 + (y − ζ)2 + log (x − ξ)2 + (y + ζ)2 . 4π (4.11) Now, from Green’s identity we obtain Z Za φ(x, y)δ(x − ξ)δ(y − ζ) dxdy = G(x, 0, ξ, ζ)φy (x, 0) dx, (4.12)

G(x, y, ξ, ζ) =

y>0

x=−a

and evaluating this at ζ = 0 we find Za φ(x, 0, t) = 2 G(ξ, 0, x, 0)φy (ξ, 0, t) dξ,

(4.13)

x=−a

which we will write as φ(x, 0, t) = L[φy (x, 0, t)], which gives φ(x, 0, t) = L[ηt + U1 ηx ]. 35

(4.14)

In the linear r´egime, we can immediately derive the equation ηtt = −ηxxxx + (∂t + U1 ∂x )L[ηt + U1 ηx ] + F(x, t), where L[f(x, t)] = π

−1

(4.15)

Za log |x − ξ|f(ξ, t) dξ

(4.16)

ξ=−a

We now write η(x, t) =

N X

ηn (t)φn (x/a),

(4.17)

n=1

where the basis functions φn are yet to be determined. Traditional boundary element methods use basis functions with only a small region of support — the usual choice is a triangular function. This has the benefit that the necessary integrals become comparatively simple, but has the disadvantage that we cannot attain high smoothness on the interior of the domain, which limits the accuracy that we can obtain. We will require relatively high smoothness on the interior of the domain so that we can calculate the bending stiffness accurately. As an alternative we can use Chebyshev polynomials as basis functions, which will produce a hybrid method, coupling spectral accuracy with the simplicity of a boundary integral method. These methods form an obvious development of the boundary integral idea and have been applied to many problems of interfacial dynamics, such as the Kelvin–Helmholtz instability of a vortex sheet, the fingering instability in Hele–Shaw flow and the Mullins– Sekerka instability of crystal growth (Hou et al., 2000). One aeroacoustical application of these methods was by Llewelyn Smith & Craster (1999), who considered the numerical solution of zero mean-flow finite baffle problems in the frequency domain and developed a Galerkin spectral solver. Unfortunately, such work is not very helpful for mean flow problems; the point of the Briggs–Bers analysis is that one cannot naively work in the frequency domain. Another (minor) difficulty with this work is that their spectral basis functions are rather unusual and it is not clear that they are an optimal choice. Thus our basis function φn is approximately the (n + 2)th Chebyshev polynomial, modified so that φn (±1) = φn0 (±1) = 0 (Boyd, 1999). The boundary conditions at the end of the plate are thus implicitly used throughout the computations. On substituting the spectral series (4.17) into the plate equation (4.15) we immediately obtain a matrix equation for the spectral coefficients of the plate acceleration η¨ n that must be solved at every timestep (Baker et al., 1980) and provides the input for one’s ODE solver of choice. Note that we need to do integrals of the form Z1 I(x) ≡ f(ξ) log |ξ − x| dξ, (4.18) −1

which can be evaluated numerically by subtracting out the singularity: Z1 Z1 I(x) = f(x) log |ξ − x| dξ + (f(ξ) − f(x)) log |ξ − x| dξ −1

−1

(4.19)

= f(x) ((1 + x) log(1 + x) + (1 − x) log(1 − x) − 2) + J(x), where J(x) can be evaluated efficiently by Clenshaw–Curtis quadrature. 36

4.2 Some results Simulations were performed with the forcing term 2

F(x/a, t) = e−µ(x/a) sin2 (π(x/a + 1)/2)) H(t) sin(ωt)e−λt ,

(4.20)

where H(t) is the Heaviside step function. This forcing term (4.20) was used to model a point force switched on at t = 0. Typically, µ was taken to be 20 and λ was varied between 0 and about 2; the idea being to model a sudden excitation of the system and to allow it to reach its own equilibrium. When λ is set to 0 there are essentially two possible responses. The plate amplitude either rapidly attains some envelope and then grows in magnitude, slowly travelling downstream (figure 4.2), or it oscillates quasiperiodically (figure 4.3). One of the most interesting results is that for many parameter values the length of the plate increases enormously, so that linear elasticity can become invalid even if the system is not unstable (see figure 4.3, which shows a quasiperiodic oscillation of amplitude 4000 on a plate of length 40). One must therefore imagine a very weak forcing term if this theory is to have any validity. The addition of nonlinear elasticity will presumably cause a saturation to some stretched state in which the inertia of the system is relatively unimportant and the forcing term is balanced by the bending moment and the nonlinear tension due to stretching. Long plates are much more unstable than both short plates and infinite plates. In general, long plates become unstable at flow speeds lower than that required for absolute instability of the infinite plate; short plates need the flow speed to be higher. This is possibly because the resonant modes of the plate appear at lower and lower frequencies as the plate lengthens (Abrahams, 1981); so as the plate lengthens the resonances come closer to the low frequency r´egime in which we are interested. Another numerical complication is that the curvature of the plate at its endpoints can increase drastically, which eventually ends up polluting the whole solution. This could be overcome with an arclength parametrisation of the plate, but it is questionable whether or not this would be worth doing, since by this stage the response of the system has entered the nonlinear r´egime. Note (figure 4.3) that we do not obtain a single frequency response, and that transients generated at startup can persist. This becomes even clearer in a calculation with nonzero λ, as shown in figure 4.4.

4.2.1 The flag problem The finite baffle has been repeatedly studied. A more interesting version of this problem is that of a flag flapping in the breeze. The flag problem is complicated by the necessity to keep track of the vortex sheet which is shed at the free end of the flag. Hopefully, the finite baffle problem should give some insight into the correct numerical treatment of the flag problem. There may, however, be a much better way to simulate the system, allowing both geometric nonlinearity of the membrane and real fluid dynamics. It is briefly explained in the next section. 37

12000

10000

8000

6000

η 4000

2000

0

-2000

-4000 -100

-50

0

50

100

x

Figure 4.2: The finite baffle: surface response with a = 100, U = 0.01, ω = 10−3 , λ = 0. 4000

3500

3000

2500

η 2000

1500

1000

500

0 0

50000

100000

150000

200000

t

Figure 4.3: The finite baffle: maximum absolute deflection on plate with a = 20, U = 0.01, ω = 10−3 , λ = 0.

38

0.001 0.0009 0.0008 0.0007 0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0 0

200

400

600

800

1000

Figure 4.4: The finite baffle: maximum absolute deflection on plate with a = 10, U = 0.05, ω = 10−3 , λ = 2.

4.3 The lattice Boltzmann method The lattice Boltzmann method is a very powerful and comparatively new tool for the simulation of viscous fluid flows, based on the solution of the Boltzmann equation of thermodynamics (Chen & Doolen, 1998). Lattice Boltzmann codes are trivial to write and (possibly surprisingly) have a computational efficiency comparable to far more sophisticated methods. We can recover the Navier–Stokes equations from a solution of the Boltzmann equation (Huang, 1988) and so can produce a solution of the Navier–Stokes equations very simply. Lattice Boltzmann codes have been developed to handle multicomponent fluids with surface tension, and some recent work by Stelitano & Rothman (2000) has extended this approach to handle elastic interfaces. Unfortunately, this method produces an interface that is typically three lattice points thick and requires reconstruction at each timestep to produce the curvature. It is not clear that this method produces a surface on which the no-slip condition is satisfied and there are also some difficulties with an unphysical surface tension term. Since we need to know something about an exact interface anyway, another possible approach would be to advect the material surface around with the flow. We could then use the boundary conditions at the surface to recover the plate equation. This would have the advantage that the no-slip condition was guaranteed, and would also remove the need for interface reconstruction at each timestep. It is in principle easy to see how this could be done; to obtain the correct boundary conditions for the Navier–Stokes equations it 39

is generally necessary to insert fake lattice points inside the wall and to produce, by some means, a particle distribution on these points. This has been recently studied by Mei et al. (1999) and it should be possible to extend this work to handle stress jumps across the boundary. Further work in this area would be both interesting and relevant.

40

Bibliography Abrahams, I. D. 1981. Scattering of sound by a heavily loaded finite elastic plate. Proc. R. Soc. Lond. A 378, 89–117. Anderson, E. et al. 1999. LAPACK Users’ Guide. 3rd ed. Philadelpia, PA: Society for Industrial and Applied Mathematics. Anderson, P. W. 1958. Absence of diffusion in certain random lattices. Phys. Rev. 109(5), 1492–1505. Baesens, C. & MacKay, R. S. 1999. Algebraic localisation of linear response in a network with algebraically decaying interaction, and application to breathers in dipole-dipole networks. Helv. Phys. Acta 72, 23–32. Baker, G. R., Meiron, D. I. & Orszag, S. A. 1980. Vortex simulations of the Rayleigh–Taylor instability. Phys. Fluids 23(8), 1485–1490. Borland, R. E. 1963. The nature of the electronic states in disordered onedimensional systems. Proc. R. Soc. Lond. A 274, 529–545. Boyd, J. P. 1999. Chebyshev and Fourier spectral methods. Dover Publications. Brazier-Smith, P. R. & Scott, J. F. 1984. Stability of fluid flow in the presence of a compliant surface. Wave Motion 6, 547–560. Briggs, R. J. 1964. Electron-stream interaction with plasmas. Cambridge, MA: MIT Press. Cairns, R. A. 1979. The role of negative energy waves in some instabilities of parallel flows. J. Fluid Mech. 92, 1–14. Chen, S. & Doolen, G. D. 1998. Lattice Boltzmann method for fluid flows. Ann. Rev. Fluid Mech. 30, 329–364. Cook, P. D. 1998. Wave transmission along a ribbed fluid-loaded membrane. Ph.D. thesis, University of Cambridge. Cooper, A. J. & Crighton, D. G. 1998a. Transmission of energy down periodically ribbed elastic structures under fluid loading: algebraic decay in the stop bands. Proc. R. Soc. Lond. A 454, 1337–1355. Cooper, A. J. & Crighton, D. G. 1998b. Transmission of energy down periodically ribbed elastic structures under fluid loading: spatial periodicity in the pass bands. Proc. R. Soc. Lond. A 454, 2893–2909. 41

Cooper, A. J. & Crighton, D. G. 1999. Response of irregularly ribbed elastic structures, under fluid loading, to localized excitation. Proc. R. Soc. Lond. A 455, 1083–1105. Crighton, D. G. 1984. Transmission of energy down periodically ribbed structures under fluid loading. Proc. R. Soc. Lond. A 394, 405–436. Crighton, D. G. 1989. The 1988 Rayleigh medal lecture: fluid loading — the interaction between sound and vibration. J. Sound Vib. 133(1), 1–27. Crighton, D. G. & Oswell, J. E. 1991. Fluid loading with mean flow. I. Response of an elastic plate to localized excitation. Phil. Trans. R. Soc. Lond. 335, 557–592. Dembowski, C. et al. 1999. Anderson localization in a string of microwave cavities. Phys. Rev. E 60(4), 3942–3948. Eatwell, G. P. & Willis, J. R. 1982. The excitation of a fluid-loaded plate stiffened by a semi-infinite array of beams. IMA J. Appl. Math. 29, 247–270. Flach, S. & Willis, C. R. 1998. Discrete breathers. Phys. Rep. 295(5), 182–264. Garrad, A. D. & Carpenter, P. W. 1982. On the aerodynamic forces involved in aeroelastic instability of two-dimensional panels in uniform incompressible flow. J. Sound Vib. 80(3), 437–439. Green, M. R. & Crighton, D. G. 2000. Fluid loaded elastic plate with mean flow: point forcing and three-dimensional effects. J. Fluid Mech. (to appear). Guo, Y. P. 1993. Effects of structural joints on sound scattering. J. Acoust. Soc. Am. 93(2), 857–863. Hodges, C. H. & Woodhouse, J. 1989. Confinement of vibration by onedimensional disorder, I: Theory of ensemble averaging. J. Sound Vib. 130(2), 237–251. Hou, T. Y., Lowengrub, J. S. & Shelley, M. J. 2000. Boundary integral methods for multicomponent fluids and multiphase materials. J. Comp. Phys. (to appear). Huang, K. 1988. Statistical Mechanics. 2nd ed. John Wiley. Kelbert, M. & Sazonov, I. 1996. Pulses and other wave processes in fluids. Modern approaches in geophysics. Kluwer Academic Publishers. Kottos, T., Izrailev, F. M. & Politi, A. 1999. Finite-length Lyapunov exponents and conductance for quasi-1D disordered solids. Physica 131D, 155–169. Kottos, T., Politi, A. & Izrailev, F. M. 1998. Finite-size corrections to Lyapunov spectra for band random matrices. J. Phys.: Condens. Mat. 10, 5965–5976. Kupfer, K., Bers, A. & Ram, A. K. 1987. The cusp map in the complex-frequency plane for absolute instabilities. Phys. Fluids 30(10), 3075–3082. de Langre, E. 2000. Ondes variqueuses absolument instables dans un canal e´ lastique. C. R. Acad. Sci. Paris s´erie IIb 328(1), 61–65. 42

de Langre, E. & Ouvrard, A. E. 1999. Absolute and convective bending instabilities in fluid-conveying pipes. J. Fluids and Structures 13, 663–680. Lawrie, J. B. 1989. An infinite elastic shell with a finite number of ring constraints. J. Sound Vib. 130(2), 189–206. Leppington, F. G., Broadbent, E. G. & Heron, K. H. 1984. Acoustic radiation from rectangular panels with constrained edges. Proc. R. Soc. Lond. A 393, 67–84. Lingwood, R. J. & Peake, N. 1999. On the causal behaviour of flow over an elastic wall. J. Fluid Mech. 396, 319–344. Llewelyn Smith, S. G. & Craster, R. V. 1999. Numerical and asymptotic approaches to scattering problems involving finite elastic plates in structural acoustics. Wave Motion 30(1), 17–41. Lucey, A. D. 1998. The excitation of waves on a flexible panel in a uniform flow. Phil. Trans. R. Soc. Lond. 356, 2999–3039. Lucey, A. D. & Carpenter, P. W. 1992. A numerical simulation of the interaction of a compliant wall and an inviscid flow. J. Fluid Mech. 234, 121–146. Luo, X. Y. & Pedley, T. J. 2000. Multiple solutions and flow limitation in collapsible channel flows. J. Fluid Mech. 420, 301–324. Mei, R., Luo, L.-S. & Shyy, W. 1999. An accurate curved boundary treatment in the lattice Boltzmann method. J. Comp. Phys. 155, 307–330. Peake, N. 1997. On the behaviour of a fluid loaded cylindrical shell with mean flow. J. Fluid Mech. 338, 387–410. Peake, N. 2000. Nonlinear stability of a fluid-loaded elastic plate with mean flow. J. Fluid Mech. (to appear). Photiadis, D. M. 1992. Anderson localization of one-dimensional wave propagation on a fluid-loaded plate. J. Acoust. Soc. Am. 91(2), 771–780. Photiadis, D. M. & Houston, B. H. 1999. Anderson localization of vibration on a framed cylindrical shell. J. Acoust. Soc. Am. 106(3), 1377–1391. Rumerman, M. L. 1975. Vibration and wave propagation in ribbed plates. J. Acoust. Soc. Am. 57(2), 370–373. Sobnack, M. B. & Crighton, D. G. 1994. Anderson localization effects in the transmission of energy down an irregularly ribbed fluid-loaded structure. Proc. R. Soc. Lond. A 444, 185–200. Spivack, M. 1991. Wave propagation in finite periodically ribbed structures. Proc. R. Soc. Lond. A 435, 615–634. Spivack, M. & Barbone, P. E. 1994. Disorder and localisation in ribbed structures with fluid loading. Proc. R. Soc. Lond. A 444, 73–89. 43

Stelitano, D. & Rothman, D. H. 2000. Fluctuations of elastic interfaces in fluids: Theory, lattice-Boltzmann model, and simulation. Phys. Rev. E 62(5), 6667–6680. Thouless, D. J. 1974. Electrons in disordered systems and the theory of localization. Phys. Rep. 13(3), 93–142. Thouless, D. J. 1978. Percolation and localization. In R. Balian, R. Maynard & G. Toulouse, eds., Ill-condensed matter. North-Holland. Tobias, S. M., Proctor, M. R. E. & Knobloch, E. 1998. Convective and absolute instabilities of fluid flows in finite geometry. Physica 113E, 43–72. Wright, S. A. 2000. Aspects of unsteady fluid-structure interaction. Ph.D. thesis, University of Cambridge.

44