## Abstract

The interaction between a fluid and a solid surface in relative motion represents a dynamical process that is central to the problem of laminar-to-turbulent transition (and consequent drag increase) for air, sea and land vehicles, as well as long-range pipelines. This problem may in principle be alleviated via a control stimulus designed to impede the generation and growth of instabilities inherent in the flow. Here, we show that phonon motion underneath a surface may be tuned to passively generate a spatio-temporal elastic deformation profile at the surface that counters these instabilities. We theoretically demonstrate this phenomenon and the underlying mechanism of frequency-dependent destructive interference of the unstable flow waves. The converse process of flow destabilization is illustrated as well. This approach provides a condensed-matter physics treatment to fluid–structure interaction and a new paradigm for flow control.

## 1. Background

The field of flow control may be traced back to prehistoric times when primitive tools and weaponry, such as spears and arrows, had features deliberately included to interact favourably with the surrounding flow [1]. These features emerged by empirical evolution rather than an understanding of the underlying physics. The scientific method for flow control appeared only at the start of the twentieth century with Prandtl's proposition of the boundary layer theory accompanied by an explanation of the physics of flow separation, and a demonstration of several experiments for boundary layer control [2]. Following this seminal contribution, the concept of flow control, especially for wall-bounded flows, emerged as a major research thrust in fluid mechanics and its development progressed over several stages, or eras, in which numerous passive and active approaches have been extensively investigated [3]. These approaches encompass streamline shaping of the surface, surface heating or cooling [4], injection of polymer additives into the flow [5], addition of ribs on the surface [6], suction and blowing [7–9] and coating of the surface with a compliant material [10–12], among others.

Regardless of the approach for stimulating a change in flow behaviour, successful intervention may be realized when particular undesirable flow structures and mechanisms are identified and clearly understood. Concerning flow transition and increase in surface (skin friction) drag, the nucleation and growth of unstable disturbances are profoundly detrimental. Owing to their wave nature, the impediment of growth of these disturbances (i.e. their stabilization) is possible and may be induced by wave cancellation or at least a degree of destructive interference [13]. Wave cancellation or superposition has been extensively investigated using active means. Results, however, have been modest especially when applied to complex conditions such as flow transition where primary disturbances give rise to residual disturbances at a variety of frequencies, phases and orientations, which render the control process intractable [14].

In a wall-bounded flow, there is a mutual dependence between the dynamical behaviour of the fluid and the solid. This dependence is shaped by the nature of the fluid–structure interaction at the interface. It is therefore conceivable, in principle, to use the fluid to affect the constitutive response of the solid and, conversely, to favourably alter the ‘character and disposition of a flow field’ [3] by tuning the elastodynamic response of the solid surface. As mentioned above, this latter notion has been explored in the literature by the utilization of surfaces with significantly compliant elastic properties. The concept was introduced in 1957 by Max O. Kramer after conducting an experiment on a motorboat in which a model with a dolphin-like skin was towed in the sea and shown to exhibit a more than 50% reduction in drag [10]. Although this result was later questioned due to lack of well-controlled experimental conditions, it has helped trigger much interest in the subject. Numerous investigations were conducted on the ensuing effects on phenomena as complex as laminar-to-turbulent transition and skin-friction drag (e.g. [11,12]). A compliant surface predominantly admits Rayleigh elastic waves along the surface [15],^{1} and due to its low stiffness allows for the possibility of large surface motion and hence significant interaction with the flow. Its main advantage stems from its passivity and simplicity, i.e. no active control devices, wires, ducts and slots, etc. are needed. This is an economically desirable benefit, in fact critical, as the energy consumed in operating active devices may often exceed the energy saved by altering the flow. The concept, however, suffers from several crucial drawbacks. The large wall motion is mostly undesirable as it increases the likelihood of surface instabilities (e.g. flutter) [16–18]. Furthermore, a considerably high compliance is generally not welcome in intense operating environments where load-bearing materials are needed. A more fundamental disadvantage is that there is no clear route to mechanistic tuning for precise frequency- and phase-dependent intervention with the flow.

## 2. Subsurface phonons

In this paper, we couple a surface-bounded flow to *subsurface phonons*,^{2} and in doing so introduce a new paradigm for passive flow control that is based on tuning the emerging coherent fluid–structure wave field. This approach, and the associated theory, draws from principles rooted in condensed matter physics. Our concept consists of introducing an elastic periodic medium, located at one or more points or regions of interest along the surface, and extending in a manner such that its spatial periodicity is along the depth, i.e. perpendicular to the surface. We will refer to this medium as a *phononic subsurface*. One realization is shown at the bottom of figure 1*a*, in which a segment of the bottom surface of a flow channel with otherwise all-rigid walls is replaced with a one-dimensional elastic phononic crystal consisting of alternating layers stacked perpendicularly underneath. This one-dimensional two-layer per unit-cell configuration is chosen for its simplicity; however, two-dimensional or three-dimensional unit cells with complex topologies such as those considered in [21,22] may be used in the future. The fundamental property that we use in a phononic material is that it possesses a frequency band structure, consisting of stop bands (band gaps) and pass bands, and that the phase of an incident wave travelling in the medium is altered based on whether it falls within a stop band or a pass band [23]. Specifically, a stop band admits out-of-phase waves and a pass band retains the phase. This characteristic is used in the proposed phononic subsurface concept to stabilize, or destabilize if desired, any flow disturbance wave interacting with the subsurface via the fluid–structure interface. Stabilization is accomplished within a stop band by inducing destructive interferences that lead to attenuation of unstable wave amplitudes in the flow. Conversely, flow destabilization is induced within a pass band by constructive interferences that amplify disturbance wave amplitudes.^{3} Both mechanisms take effect robustly due to the intrinsic correlation between band structure and phase within a phononic material, thus individual wave tracking is not needed.

### (a) Subsurface dynamical characteristics

We will consider a Tollmien–Schlichting (TS) [24,25] wave as an example of an unstable flow wave. A TS wave represents a streamwise instability that arises in a viscous boundary layer and also suitably predicts the evolution of a two-dimensional instability in parallel (and nearly parallel) shear flows. TS waves have been observed in laboratory experiments for boundary layers [26,27] and channel flows [28], in which the unstable waves were generated using a vibrating magnetic ribbon. We first consider a TS wave with a frequency *ω*_{TS}=1690 Hz. To stabilize this wave, we design the one-dimensional phononic subsurface unit cell to exhibit a stop band over a frequency range encompassing its frequency. The emerging unit cell consists of a layer of ABS polymer (90% volume fraction) and a layer of aluminium (10% volume fraction) and has a total length *L*_{UC}=40 cm. The cross-sectional area is assumed to be uniform and is therefore eliminated from the analysis. The resulting composite is relatively stiff and exhibits modest damping. The inertial, elastic and damping properties of the selected materials are provided in appendix A. The choice of unit-cell length is not constrained in the current investigation; future work will aim to shorten it by exploiting established band-gap-lowering mechanisms, such as using local resonators [22,29]. By employing Bloch's theorem [30] (see problem 1 in appendix A), we obtain dispersion curves for this unit-cell configuration, which appropriately comprise a stop band covering the range 1652≤*ω*≤2102 Hz, as shown in figure 1*b*. Pass band behaviour is observed below 1652 Hz and from 2102 Hz to 3376 Hz.

Next, we consider the phononic subsurface as an isolated finite structure, i.e. as shown in figure 1*a* except without the flow channel attached (see problem 2 in appendix A). This structure spans 10 unit cells and is fixed at the bottom end. Furthermore, the layering sequence is chosen such that the topmost layer is ABS polymer. We apply a harmonic force, *ω**≤3376 Hz (the position along the structure, time, force amplitude and displacement amplitude are denoted by *s*, *t*, *c* shows this information as a normalized frequency response function, *b*, except for an anomalous resonant peak at a frequency of *ω*_{Truncation}=1685.2 Hz which falls inside the first stop band. This resonance stems from the symmetry-breaking periodicity truncation [31].

To verify the correlation described above between the band structure and the phasing, we perform yet another analysis on the isolated 10-unit-cell phononic structure. Now an elastodynamic time simulation is executed (see problem 3 in appendix A). The same harmonic forcing is again applied at the top and the displacement response, now in transient form, is recorded at the same location. The time-averaged dot product of these two signals in normalized form provides a quantitative measure of the relative phase between the forcing and the response. This quantity is defined as
*d* for a range of excitation frequencies (considering a total simulation time of *t*_{T}=1.18 s). Beyond this time, no further changes in the time-averaged dot product is observed. We note that indeed, within the stop band, the two signals are predominantly out of phase (〈*ϕ*(*ω**)〉<0), and within the pass bands they are in phase (〈*ϕ*(*ω**)〉>0) between each resonance and its preceding anti-resonance. It is clear that the range of out-of-phase behaviour is larger within the stop band. This provides us with an unprecedented route for inhibiting the growth of an unstable TS wave via de-phasing, namely to ensure that 〈*ϕ*(*ω**)〉 is negative and has a high absolute value at *ω*_{TS}—these conditions are met at frequencies close to and to the right of *ω*_{Truncation}.

While operating at a stop band will provide the proper phase conditions for stabilizing a growing TS wave, it is also necessary for the surface to be in considerable motion for a significant effect to be transmitted to the flow field.^{4} The same can be said about the opposite function of pass-band-triggered destabilization (which is considered in more detail later). The displacement amplitude of the interface, |*η*(0,*t*)|, should therefore also be accounted for in the tuning scheme. On this basis, we establish a performance metric, defined as *c*,*d*, we obtain a prediction of the performance metric, *P*=*P*_{P}, which is plotted as a function of *ω** in figure 1*e*. It is notable that the value of *P* is strongly negative just to the right of the truncation resonance frequency and retains a relatively strong negative value over a substantial frequency range on that side of *ω*_{Truncation}. The left side, to the contrary, experiences immediate positive *P* values that rise sharply and subsequently oscillate as the frequency decreases. While truncation modes are usually undesirable in conventional phononic crystal applications, here they provide a favourable effect as they provide the opportunity of large-amplitude response within a stop band. If the interest is to have a large frequency range for stabilization, then the unit cell should be designed to (i) exhibit the largest stop band possible encompassing the TS waves' frequency range and (ii) have a truncation resonance located to the left of the TS waves' frequency range and furthest to the left within the stop band to cover as wide a range as possible. If the interest is in destabilization, then the unit cell should be designed to exhibit the highest possible value of *P*—in the positive domain—at or near the TS wave frequency (or frequency range).

It is notable that *P* is also strongly negative at very low frequencies and is most negative between the first resonance and anti-resonance frequencies. In principle, these low-frequency windows could also be chosen as operating ranges for stabilization. However, utilization of a truncation resonance in a stop band provides more robust conditions. This is because both the stop band and the truncation resonance are unit-cell properties [31] and are therefore insensitive to global changes to the finite structure. By contrast, all pass band resonances are global characteristics that are sensitive to structural alterations such as contact with other components and changes in the boundary conditions. Furthermore, a design principle that is focused on intrinsic properties is more systematic, efficient and general in scope.

### (b) Subsurface interaction with the flow

In this section, coupled fluid–structure direct numerical simulations are executed to examine the actual behaviour of the phononic subsurface and its performance in altering the flow response via the interface motion. This represents problem 4 in appendix A. The TS wave is incorporated as an evolving disturbance in a fully developed plane Poiseuille (channel) flow driven by a mean pressure gradient. The channel is formed by parallel walls at the bottom and at the top with periodic boundary conditions applied in the spanwise *z*-direction, and a buffer (sponge) layer is used to model the outflow [9]. In one arrangement serving as a reference, the two walls are rigid everywhere along the channel. In the main arrangement, the phononic subsurface is installed as described earlier. It is positioned approximately at the centre of the channel and fully covers the spanwise length. We also examine the extreme arrangements of installing a structure at the same location and of the same length and boundary conditions but consisting of either entirely ABS polymer or entirely aluminium. Figure 1*a* illustrates the overall fluid–structure set-up for all four arrangements.

In the fluid part of the coupled model, the base flow is an exact solution of the Navier–Stokes equation [32,33]. An unstable spatial solution (eigenfunction) of the Orr–Sommerfeld equation [34–36] is superimposed at the inflow boundary of the channel to perturb the parabolic base velocity, precisely modelling the conditions in typical laboratory experiments [28]. We consider water as the working fluid, a Reynolds number *Re*=7500, and a non-dimensional unstable frequency *ω*_{R}=0.25. Details on the fluid model, the numerical integration schemes for the fluid and solid domains, and the coupling scheme are provided in appendix A.

The results from these simulations are remarkable. The interface steady-state vibration amplitude, interface phasing quantity and actual performance metric, *P*=*P*_{A}, for four values of the TS frequency, *ω*_{TS}=1600,1690,1700 and 1800 Hz (cases A–D), are plotted as red dots in figure 1*c*,*d*,*e*, respectively, and are in excellent agreement with the response of the artificially forced isolated structure. This indicates that the phononic subsurface may be designed *a priori*, and independently, from the coupled fluid–structure simulations. This is a key advantage because it allows for the design process to be performed outside the simulation loop, thus providing orders of magnitude of savings in computational time and expense.

We now focus our attention on (i) case B, for which *ω*_{TS}=1690 Hz lies within a stop band and falls to the right and close to the truncation frequency *ω*_{Truncation} and, for comparison, (ii) case A, for which *ω*_{TS}=1600 Hz lies within a pass band and falls between an anti-resonance frequency and a resonance frequency. In selecting these two cases, the intention is to provide diverse examples. Other than the stop band versus pass band contrast, case B exhibits a high value of |*P*|, while case A exhibits a low value of |*P*|. We first observe that the coupled simulation displacement field inside the phononic subsurface obeys perfectly the expected stop band (attenuating) and pass band (propagating) behaviour that would appear in an otherwise stand-alone phononic crystal subjected to harmonic excitation [23] (see §3a for a detailed energy analysis).

Turning to the actual flow behaviour for the same two cases, we consider two quantities in the fluid: the space-averaged time-dependent perturbation kinetic energy, *T**_{f} [in units of J (joules)] and the time-averaged position-dependent perturbation kinetic energy, *K**_{f} [in units of J m^{−1} (joules per metre)]. The former quantity, which is based on a volume integration, is given by
**u**=(*u*,*v*,*w*) represents the velocity components in the streamwise, *x*, wall-normal, *y*, and the spanwise, *z*, directions, respectively; *ρ*_{f} is the fluid density; *L*_{x} and *L*_{y} are the channel-base dimensions; and *δ* is the channel half-height. The symbols (∧) and (*) are used to represent a fluctuating component of the velocity field and a dimensional quantity, respectively (see appendix A for the detailed velocity decomposition formulation). The term *t*. The volume of integration over the bottom half of the channel is bounded by its physical domain, i.e. the buffer domain is excluded from the integration. The *K**_{f}, on the other hand, is based on a time–area integration and is a function of the dimensional streamwise coordinate *x**. It is defined as
*t*_{beg} denotes the beginning of the time-averaging window and *t*_{end} denotes the end of the time-averaging window. In our calculations,

The *T**_{f} quantity is plotted in figure 2*a*,*b* and is observed to correlate, most stringently, with the *P* metric. In the stop band frequency case (figure 2*a*), the perturbation kinetic energy, while oscillating, is predominantly lower in value when a phononic subsurface is installed than with all-rigid walls. This is consistent with the anticipated performance metric value in figure 1*e*, which is a negative *P* with a high absolute value. By contrast, the two extreme arrangements of an all-ABS and all-aluminium subsurface display negligible effect. Figure 2*b* shows the corresponding results for the pass band frequency case where we observe a predominantly positive increase in the perturbation kinetic energy—consistent with a positive *P*—and generally a significantly smaller effect overall—consistent with a low |*P*|. Also of significance is that neither the flow stop band nor the pass band responses are decaying in time, thus the perturbation kinetic energy behaviour observed above takes effect in a steady-state and sustainable fashion. To explicitly examine the role of phase in the phenomenon observed, figure 2*c* monitors the relative phase between the flow perturbation vertical velocity right above the control segment at *y*=0.0108*δ* when a phononic subsurface is installed and that with all-rigid walls. We see that, in the stop band case, the introduction of the phononic subsurface alters the phase significantly compared with the pass band case, where the phase is preserved on average. This result confirms the existence of destructive interferences in the vertical velocity components of the flow near the interface for case B, and, conversely, the existence of constructive interferences for case A.^{5} Hence the influence of the band structure of the phononic subsurface is shown to have extended to the flow.

It is notable in figure 2*a*,*c* that the phononic subsurface case B signals appear to oscillate at a very low frequency of roughly 7 Hz. This is a beat frequency reflecting the close proximity of *ω*_{TS} and *ω*_{Truncation}. From the coupled model, *ω*_{Truncation}≈1683 Hz, which is within only 2 Hz of the uncoupled value (providing additional reassurance on the excellent accuracy of the numerical simulation).

The *K**_{f} quantity is plotted in figure 3*a*,*b* for cases B and A, respectively. Compared with the reference all-rigid-wall case, it is observed that *K**_{f} is lower for case B, and higher for case A, across the length of the phononic subsurface—both results being consistent with the *T**_{f} response and what is suggested by the performance metric. Furthermore, upstream effects are displayed which are particularly strong in figure 3*a*. In that region, the gradient of *K**_{f} with respect to position over the leading edge of the phononic subsurface is mostly negative for case B, and mostly positive for case A, indicating stabilization and destabilization, respectively. By installing a double phononic subsurface (i.e. two contiguous segments), the region of negative gradient roughly doubles for case B. At the downstream end of the phononic subsurface, a high positive *K**_{f} gradient is observed in case B. This is a manifestation of the underlying stabilization mechanism being spatially local. However, as demonstrated by the double-segment results, the control region is easily extendable by simply stretching the area covered by the phononic subsurface. A 13% and 24% maximum reduction in the instability kinetic energy is realized for case B when, respectively, a single or double phononic subsurface segment is installed. On the other hand, case A displays a corresponding 1% and 1.5% increase in *K**_{f}. The high percentage values in case B and the low percentage values in case A are due to a high |*P*| and a low |*P*| for each case, respectively. An analysis of the flux and production rate of the perturbation energy for each of these two cases is presented in §3b.

## 3. Further analysis

### (a) Phononic subsurface energetics

In this section, we analyse the time evolution and spatial distribution of the displacement field inside the phononic subsurface in the context of the coupled fluid–structure simulations described above. Out of the four cases of different TS frequencies considered in figure 1, here we again focus on case B and case A. The stop band nature of case B suggests that the response takes a standing-wave, attenuating form (i.e. decaying with depth), and, conversely, the pass band nature of case A suggests a propagating-wave, non-decaying profile. Furthermore, as seen in figure 1*e* and pointed out earlier, the value of the performance metric, both the predicted, *P*_{P} (black solid line), and the actual, *P*_{A} (red dots), indicate a high negative value for case B and a low positive value for case A. As such, we expect in case B to see a large overall response at the interface between the fluid and the phononic subsurface, and in case A a small overall response. In figure 4, we show average quantities of the total energy, defined as *a*,*b* presents a space-averaged total elastodynamic (kinetic and potential) energy for case B and case A, respectively, displaying oscillating patterns and a two-orders-of-magnitude larger overall response for case B, as anticipated. Also notable is that, in both cases, the motion of the phononic subsurface, including the top end which interfaces with the flow, is stable—i.e. there is no occurrence of flutter. Figure 4*c*,*d* presents a time-averaged total elastodynamic energy for case B and case A, respectively, for the top half of the structure and over the total simulation time. The peaks correspond to the regions occupied by aluminium layers, which exhibit a faster speed of sound. The energy profiles for cases B and A are attenuating and non-attenuating, respectively, which is perfectly consistent with the characteristic stop band/pass band behaviour, as predicted.

### (b) Flow energetics

From a fluid dynamics perspective, the fundamental mechanism that drives the dynamical phenomenon directly observed in figures 2 and 3 can be further elucidated by considering the equation for the perturbation energy flux integrated across the full channel height, written as
*ω*_{z} denotes spanwise vorticity, subscript ‘w’ refers to a quantity evaluated at the wall, and 〈〉 again expresses a time average. The symbol ()_{b} is used to represent a base component of the velocity field. Following [37], the first two terms on the left-hand side are grouped as the integrated mechanical energy flux. The most significant contributions to the energy budget are from the convection term (I) and flow work (II), which are balanced mainly by the kinetic energy production term (i) and, to a much lesser extent, by the irreversible pressure–work term (iii).

In figure 5, we present the streamwise evolution of these quantities over the phononic subsurface obtained from the coupled simulations. From this figure, it is clear that there is a substantial decrease in the mechanical energy (terms I + II), with a maximum decrease of about 7%, compared with the all-rigid-wall arrangement. Strong upstream and downstream effects are also evident, as mentioned earlier. These effects are observed to extend to regions slightly outside the range of the phononic subsurface. Pronounced upstream/downstream effects are observed in the compliant wall simulations of Davies & Carpenter [37] which they attribute to possible interferences between incident and scattered waves at the front and trailing edges of their compliant panel along the flow direction. This behaviour may also be extant in the present simulations. Furthermore, it is interesting to note that there is a strong indication of reduction of mechanical energy, but no discernible net decrease in the production term. However, a true comparison should be made between the axial rate of change of the mechanical energy term (energy flux) and the production term as these are the actual terms that appear in equation (3.1). For the present simulation with the phononic subsurface, the mechanical energy flux term is indicated by the black dotted line; it is apparent that this term follows the trend of the production term very faithfully, i.e. a negative region followed by a positive region such that the aggregate effect is approximately null flux, as opposed to the all-rigid-wall arrangement where the flux term (indicated by the slope of the dashed blue line) is always positive. Hence, the phononic subsurface instigates changes in the production term and reduces the disturbance kinetic energy by modifying its flux (as caused by changes in the relative phase between the

Finally, we consider the bottom-half integrated rate of production of perturbation energy, which is defined as

## 4. Conclusion

In conclusion, we have shown that phasing relations may be structurally encoded at the phononic subsurface unit-cell level to enable stability control in a wall-bounded flow. Our results reveal that, when the frequency of a TS instability wave is accommodated to the right of *ω*_{Truncation} within a stop band, the wave stabilizes, and when accommodated appropriately within a pass band, the wave destabilizes. The underlying mechanism, which in principle is applicable to other types of instability waves as well, is based on the phenomenon of extension of wave interferences from the phononic subsurface into the adjacent flow. By accounting for the phononic subsurface frequency response function, the intensity of the fluid–structure interaction is also controlled. Figure 7 shows that the predicted performance metric (which is obtained independent of the coupled fluid–structure simulations) correlates precisely with the actual performance metric and, strikingly, with the actual change in the perturbation kinetic energy, both as computed from the coupled simulations. This represents a mechanistic approach for passive flow stability control that is effective, efficient, predictable and repeatable. And, of practical importance, the approach does not depend on having a highly compliant or a highly damped surface, thus simultaneously maintaining surface stability while providing powerful and sustained flow control. While the current implementation uses a single (or double) phononic subsurface line segment, the dispersion characteristics of the control region may be varied in space with pointwise local resolution, i.e. the flow stabilization/destabilization activity may be effectively ‘sculptured’ along the surface. This spatial versatility stems from the fact that the phonon band structure dynamics is generated via elastic wave propagation perpendicular to the surface, and not parallel to it.

The subsurface phonon-induced processes demonstrated in this work represent a fundamental fluid–structure mechanism that draws on concepts from phonon physics, including Bloch's theorem and the Brillouin zone. In conjunction with the Orr–Sommerfeld linear stability theory characterizing shear flow waves, the present theory of subsurface phonons provides an integrated and unified fluid–structure wave dynamics framework and allows, for the first time, wave motion in the flow and within the interacting elastic solid to be fully synchronized.

## Data accessibility

All data are included in the manuscript.

## Funding statement

This research was funded by the National Science Foundation grant no. 1131802 under the following programmes: Interdisciplinary Research (IDR), Mechanics of Materials, Dynamical Systems, Fluid Dynamics and Biotechnology, Biochemical, and Biomass Engineering. A CU-Boulder Innovation Seed grant funded the earlier stages of the research.

## Authors contributions

M.I.H. and S.B. conceived of the proposed concept and designed the research; M.I.H., S.B., O.R.B. and A.K. performed the research and analysed the results; M.I.H. and S.B. wrote the paper. M.I.H. and O.R.B. contributed mostly to the phononics aspects; S.B. and A.K. contributed mostly to the flow aspects.

## Conflict of interests

The authors have no conflict of interests to declare.

## Acknowledgements

The authors are grateful to the National Science Foundation programme managers Dr Bruce Kramer and Dr Glaucio Paulino for their support.

## Appendix A. Models, methods and parameters

In this work, we consider four separate problems. The first (problem 1) focuses on a single phononic subsurface unit cell modelled as an infinite one-dimensional periodic medium for which we obtain a band structure diagram. This diagram describes the dispersion characteristics of elastic wave propagation in the phononic subsurface. We consider wave propagation along only the direction of the periodicity. Problem 2 is concerned with the steady-state forced dynamic response, for a range of harmonic excitation frequencies, of the phononic subsurface when viewed as a structure composed of a finite number of unit cells with free and fixed boundary conditions at the top and bottom ends, respectively. In problem 3, we simulate elastic wave propagation in the same finite structure and under the same forcing conditions considered in the second problem. Finally, in problem 4, we consider a full simulation of the coupled fluid–structure domains. This appendix covers the modelling and computation methods, as well as the model parameters, used to investigate these four problems. The results and discussions are presented in §2 and §3. All notation is defined here regardless of prior appearance in the main sections of the paper.

**(a) Solid models**

The governing equation of motion for longitudinal wave propagation in a one-dimensional heterogeneous linear elastic solid with uniform cross-sectional area is
*η*=*η*(*s*,*t*), *f*=*f*(*s*,*t*), *ρ*_{s}=*ρ*(*s*) and *E*=*E*(*s*) denote longitudinal displacement, external force, solid density and elastic modulus, respectively. The position coordinate and time are denoted by *s* and *t*, respectively. The operation (.)_{,s} denotes differentiation with respect to position, while a superposed dot denotes differentiation with respect to time.

**(i) Problem 1: phononic subsurface unit-cell dispersion model**

To obtain the band structure for a given phononic subsurface unit-cell configuration, we set *f*=0 in equation (A 1) and assume a Bloch solution [30] in the form *κ* is the wavenumber, *ω* is the frequency and *L*_{UC} is the unit-cell length.

**(ii) Problem 2: phononic subsurface 10-unit-cell structural dynamics model**

In problem 2, a boundary value problem following equation (A 1) is set up for a finite phononic subsurface structure composed of *n*_{c} repeated unit cells of the type considered in problem 1. The boundary conditions chosen are *η*_{,s}(0,*t*)=0 and *η*(*l*,*t*)=0 for the top and bottom ends, respectively, where *l*=*n*_{c}*L*_{UC}. A harmonic forcing function is applied only at the top end, i.e. *f*(*s*,*t*)=0 for *s*>0, where *ω** is the excitation frequency and *n*_{c}=10, although in principle a smaller number of cells may be used [23].

**(iii) Problems 3 and 4: phononic subsurface 10-unit-cell wave propagation model**

In the third and fourth problems, we re-examine problem 2 but now as an initial boundary value problem, whereby no assumptions are made on the temporal dependency of the displacement field. In problem 3, we consider the initial conditions *η*(*s*,0)=0 and apply a point harmonic force *t*≤*t*_{T}, where *t*_{T} is the total simulation time. As in problem 2, we set *f*(*s*,*t*)=0 for *s*>0. In problem 4, we consider the same initial conditions and do not apply any external forcing (the forcing is induced by the interacting flow).

**(iv) Solid-domain numerical schemes**

We numerically analyse the one-dimensional solid domain representing the phononic subsurface using the finite-element (FE) method using one-dimensional 2-node iso-parametric elements [23]. Upon discretization of the unit cell in problem 1, we obtain a mass matrix, *q*_{1} and *q*_{2} are damping constants. The dispersion band structure is obtained by solving the eigenvalue problem *κ*≤*π*/*L*_{UC} (see [38] for a solution scheme for this form of eigenvalue problem). The number of nodes in the unit cell is denoted by *n*_{m}.

In problems 2–4, we use the FE method to spatially discretize the entire 10-unit-cell structure (in problem 4, we also consider homogeneous structures for comparison). This yields a global mass matrix, **M**, and a global stiffness matrix, **K**. For problem 2, we obtain the matrix problem *s*=0. The number of nodes along the entire structure is *n*_{s}=*n*_{c}(*n*_{m}−1)+1.

For problems 3 and 4, the second-order Newmark scheme is used for the time integration. The FE matrix equations following the scheme are
**D**^{(i)}, **V**^{(i)}, **A**^{(i)} and **F**^{(i)} represent the displacement, velocity, acceleration and force vectors at the *i*th time-step, Δ*t* is the time-step increment and **T**=**M**+(Δ*t*/2)**C**. Here, **F**^{(i)} vector are set to zero. In figure 1, *n*=4×10^{6}.

The value of Δ*t* is chosen such that the simulation is both accurate and stable for a given forcing frequency *ω**. By selecting *β*=0, the scheme is rendered fully explicit with an accuracy of the order of *O*(Δ*t*^{2}). The solution procedure is as follows. At each time step *i*, the displacement at the following time step *i*+1 is first obtained, from equation (A 2). After which, equation (A 4) is solved for the unknown acceleration **A**^{(i+1)}. Because of the incorporation of damping, **T** is generally not a diagonal matrix and thus obtaining **T**^{−1} is computationally expensive and introduces numerical errors, especially for large size matrices. However, using the property that **T** is a tridiagonal matrix, we use the compact-storage Thomas algorithm, which is an efficient form of Gaussian elimination suited for tridiagonal systems. Finally, equation (A 3) is used to obtain the updated velocity. This process is repeated, allowing the simulation to march forward in time until it reaches *t*_{T}. It is noteworthy that, if **C**=**0**, equation (A 4) becomes
**M** is a diagonal matrix.

Postprocessing of the results is performed to obtain quantities such as the total elastodynamic energy, which are plotted in figure 4. The space-averaged combination of kinetic energy and potential energy over the entire structure is defined and computed as follows:

Similarly, the time-averaged combination of kinetic energy and potential energy is
*n*_{e}=*n*_{c}(*n*_{m}−1) is the total number of elements in the structure.

**(b) Fluid model**

**(i) Overview**

We consider spatially evolving disturbances in fully developed plane Poiseuille (channel) flow driven by a mean pressure gradient between two nominally rigid, parallel walls. The mean velocity for this laminar flow is an exact solution of the Navier–Stokes equation [32,33] and is used as the base flow in our mathematical model. The hydrodynamic stability of plane Poiseuille flow is governed by the Orr–Sommerfeld equation [34–36] derived by linearizing the Navier–Stokes equation with the normal mode assumption. Growing eigensolutions of the Orr–Sommerfeld equation are known as the TS [24,25] waves, which predict the evolution of two-dimensional instabilities in parallel (and nearly parallel) shear flows. TS waves have been observed in laboratory experiments for boundary layers [26,27] and channel flows [28], in which the unstable waves were generated using a vibrating magnetic ribbon. Excellent agreement between experimental measurements and Orr–Sommerfeld solutions lends very strong evidence to the accuracy and viability of the linear stability theory. In this work, we impose unstable spatial solutions (eigenfunctions) of the Orr–Sommerfeld equation at the inflow boundary of the channel to perturb the parabolic base velocity, precisely modelling the conditions in typical laboratory experiments.

**(ii) Coupled fluid–structure simulation model**

The mathematical model is based on the time-dependent, three-dimensional, incompressible continuity and the Navier–Stokes equations non-dimensionalized by the channel half-height *δ* and the centreline velocity *U*_{c}. These equations read, respectively,
*Re*=*U*_{c}*δ*/*ν*, *ν* is the kinematic viscosity, and the vector **u**=(*u*,*v*,*w*) represents the velocity components in the streamwise, *x*, wall-normal, *y*, and the spanwise, *z*, directions, respectively; also, *p* is the non-dimensional pressure and *t* (in this section) represents non-dimensional time. In equations (A 9) and (A 10), the velocity vector can be decomposed into a base component **u**_{b} and a fluctuating component *a*
*b*
*c*
where *A*_{2D} is the amplitude of the two-dimensional disturbance, *ω*_{R} is the real frequency and **u**_{e2D} is a function to be prescribed. Periodic boundary conditions are applied in the spanwise, *z*, direction. Equation (A 12*a*) is the exact solution of the governing equations for plane Poiseuille flow, where the dimensional channel height spans from 0→2*δ*. For the spatially evolving channel flow under consideration, the outflow boundary conditions are imposed in a non-reflective buffer domain appended to the physical domain [9,39,40]. In this buffer domain, the streamwise perturbations in the governing equations are smoothly brought to zero to ensure strictly outgoing waves.

In equation (A 12*b*), **u**_{e2D} is obtained from the solution of the Orr–Sommerfeld equation, given as
^{′} represents differentiation in the wall-normal direction. The complex wave speed is defined as *c*=−*ω*_{R}/*α*, where *α*=*α*_{R}+i*α*_{I} denotes a complex wavenumber.

For spatially evolving channel flow, given real circular frequency *ω*_{R} and *Re*, equation (A 13) is solved for the complex eigenvalue *α* [41], where −*α*_{I} is the amplification rate of the corresponding eigenfunctions *b*)), which, in turn, are imposed as initial disturbances at the inflow boundary in the simulations. Accordingly, disturbances will be unstable (growing in space) if −*α*_{I}>0.

In the coupled system, i.e. problem 4, a portion of the bottom wall of the channel is replaced by the structure of interest covering a streamwise distance spanning from *x*_{s} to *x*_{e}. The structure also extends in the spanwise *z*-direction across the whole span of the computational domain without any variations, consistent with the flow periodicity in that direction. At every time step, the dimensional wall pressure *ρ*_{f} is the fluid density) is calculated on the midpoint between *x*_{s} and *x*_{e} on the fluid–structure interface. This pressure acts on the structure as a force and the resultant displacement, *η*(0,*t*), and velocity, *a*
and
*b*
These boundary conditions are obtained assuming small displacements and allowing wall motion only along the wall-normal *y*-direction. Hence, it is critical that *η*≪*δ* must be maintained throughout the computations. Additionally, the spanwise velocity, *w*, is zero at the fluid–structure interface.

**(iii) Fluid-domain numerical scheme**

A time-splitting procedure [9,39,40] is implemented to numerically integrate the fluid-domain governing equations where the Crank–Nicolson implicit method is used for the wall-normal direction diffusion term; all the other terms are treated explicitly with the Adams–Bashforth method. Spatial discretization is done by fourth-order central differences on a rectangular mesh system staggered and stretched along the wall-normal direction to accurately capture the fast flow-field gradients. The time-splitting procedure first advances the velocity field, lagging the pressure gradient term. The corrector step then computes the pressure field satisfying continuity at the advanced time level. Using this scheme, we have obtained excellent agreement with the linear theory, with a maximum deviation of 0.05% in the predicted perturbation energy growth [9].

**(c) Model parameters**

**(i) Phononic subsurface**

The unit cell of the phononic subsurface configuration considered is composed of two layers, one made out of aluminium and the other of ABS polymer with a respective volume fraction of 9 : 1. The density, *ρ*_{s}, and Young's modulus, *E*, for each of these two constituent materials are as follows: *ρ*_{Al}=2700 kg m^{−3}, *ρ*_{ABS}=1040 kg m^{−3}, *E*_{Al}=68.8 GPa and *E*_{ABS}=2.4 GPa. For the chosen volume fraction, the average static Young's modulus of the composite is 2.66 GPa, which is significantly high compared with that of conventional compliant coatings. For the Rayleigh proportional damping model, we select *q*_{1}=0 for both materials and *q*_{2}=6×10^{−9} for the aluminium and *q*_{2}=6×10^{−8} for the ABS polymer (preserving physical units). This translates into damping ratios of the order of *ζ*=0.001 for the composite phononic subsurface material, which represents a relatively modest value and is consistent with experimental observations for common structural materials [42]. The unit-cell length (i.e. depth into the surface) is *L*_{UC}=40 cm; the phononic subsurface structure as a whole consists of *n*_{c}=10 cells with a total length of *l*=4 m.

In the coupled simulations, the unit-cell width covers the streamwise range *x*_{s}≤*x*≤*x*_{e}, where *x*_{s}=2.926×10^{−3} m and *x*_{e}=3.602×10^{−3} m when a single segment is installed, and *x*_{s}=2.589×10^{−3} m and *x*_{e}=3.939×10^{−3} m for a double segment.

The phononic subsurface unit-cell FE model used for the dispersion curves calculations consists of 200 elements. In the uncoupled and coupled frequency response function calculations and elastodynamic wave propagation simulations, 100 FEs are used per unit cell to discretize the structure. The time step for the all-solid-domain simulations is Δ*t*=2.95×10^{−7} s, leading to a total simulation time of *t*_{T}=1.18 s.

**(ii) Channel flow**

The flow simulations are based on a Reynolds number *Re*=7500 and a non-dimensional frequency *ω*_{R}=0.25. This yields the least-damped eigenvalue of the Orr–Sommerfeld equation (equation (A 13)) *α*=1.0004−i0.0062. This eigenvalue represents an unstable spatial mode with a growth rate of −*α*_{I}=0.0062. The real frequency *f*_{R} is obtained by
*U*_{c} or *δ* is prescribed for a given fluid viscosity *ν* and *Re*. All simulations are done for water with *ν*=1×10^{−6} m^{2} s^{−1} and using *δ*=4.2×10^{−4} m. The centreline velocity for the primary simulation (case B) is *U*_{c}=17.8 m s^{−1}, which yields the dimensional frequency *f*_{R}=1690 Hz. The value of *U*_{c} is adjusted accordingly to obtain *f*_{R}=1600 Hz (case A), *f*_{R}=1700 Hz (case C) and *f*_{R}=1800 Hz (case D), respectively. While the present analysis is based on water, the geometric scales and the physical fluid properties considered approximately correspond to a free stream velocity of 150 m s^{−1} in an air boundary layer developing on a wing of a typical airliner.

In all simulations, the dimensional fluid domain lengths were set as *L*_{x}=20*δ*, *L*_{y}=2*δ* and *L*_{z}=2*πδ* for the streamwise, wall-normal and spanwise directions, respectively. The buffer domain length is approximately 30% of the total channel length, effectively preventing reflections at the outflow boundary. The computations were performed on the Kraken supercomputer at the National Institute for Computational Sciences, Oakridge, TN, USA, using a 225(*x*)×65(*y*)×8(*z*) mesh in the entire domain availing 60 grid points per TS wavelength (λ_{x}=2*π*/*α*_{R}) in the streamwise direction. Also, the time step Δ*t* was chosen such that there are 2000 time steps for every period of the TS wave. This spatial and temporal resolution is adequate in predicting the growth rates of unstable waves governed by the linear theory. The chosen length of the phononic subsurface structure along the streamwise direction, *L*_{s}=1.6095*δ*, corresponds to roughly one quarter wavelength of the primary disturbance wavelength λ_{x}.

A frequently used criterion states that a surface is hydraulically smooth if its roughness Reynolds number is less than 25, i.e. *Re*_{k}=*U*_{k}*k*/*ν*<25, where *U*_{k} is the local undisturbed velocity, *k* is the maximum permissible height (i.e. *y*-value) of the roughness element and *ν* is the kinematic viscosity [43]. In our calculations, we find that the maximum displacement of the fluid–phononic subsurface interface is *k*=1.54×10^{−5} m. The interface displacement can therefore be considered sufficiently small and does not act as a surface roughness element, and hence does not promote or enhance transition.

## Footnotes

↵1 In some studies, the surface is modelled as a taut membrane admitting scalar transverse waves [12,16].

↵2 The proposed concept uses classical elastic waves in a periodic medium in analogy to phonon motion in a crystal lattice [19].

↵3 In a finite structure, the phase is retained within a pass band at alternating frequency windows correlating with the anti-resonance and resonance frequencies.

↵4 Considerable motion in a phononic subsurface is still within the bounds of small-amplitude (infinitesimal) vibrations, which is in contrast to the typically large-amplitude (finite) deformation experienced in compliant surfaces.

↵5 The wave interferences subsequently also affect the horizontal components of the flow due to the inherent coupling between the components.

- Received December 1, 2014.
- Accepted March 24, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.