MULTIPLE INPUT MULTIPLE OUTPUT PHONONIC SUBSURFACES FOR PASSIVE BOUNDARY LAYER TRANSITION DELAY

Abstract

A multi-input, multi-output phononic system including a first interface surface and a second interface surface that respond to at least one of a pressure gradient or a velocity gradient in a wave of a turbulent fluid flow or a laminar fluid flow, the pressure gradient or the velocity gradient associated with complex motion of the flow exhibiting a plurality of frequencies exerted on one or more of the interface surfaces; and a subsurface feature extending from the interface surfaces, the subsurface feature comprising a phononic crystal or locally resonant metamaterial adapted to receive one or more of the pressure gradient or the velocity gradient from the fluid flow via the interface surfaces and to alter one or more of a phase and an amplitude of a plurality of frequency components of the fluid flow.

Description

FIELD OF THE INVENTION

The present invention relates generally to systems and methods for reducing skin friction of a surface and, more particularly, to systems and methods for reducing skin friction of a surface using multi-input, multi-output phononic subsurfaces for passive boundary layer transition delay.

BACKGROUND OF THE INVENTION

Skin friction drag is caused by the frictional resistance of an object moving through a fluid, and it is responsible for billions of dollars in energy costs each year. Understanding fluid-solid interactions and developing methods to reduce skin friction drag would lead to a reduction in energy consumption for the gas, aerospace, and shipping industries. With a focus on aerospace applications, both active and passive flow control methods have generally been proposed. Some passive approaches include contouring airfoils, adding riblets or dimples, using super hydrophobic surfaces, including micro/nano morphologies, or coating with compliant surfaces. While passive techniques can be relatively simple to implement, their drag reduction effectiveness can be quite modest and limited to specific flight conditions.

Alternatively or additionally, active approaches can be implemented. Active approaches can significantly delay a transition from laminar to turbulent flow along a surface and can thus decrease fuel costs in a variety of flight conditions. But active methods may require some type of energy input and higher maintenance costs. These inputs and costs can outweigh any benefit received from the active approach. Some active approaches include acoustic wave generation, piezoelectric actuators, periodic surface heating and cooling, wall blowing or sucking, or plasma actuators. Because of the potential inefficiencies of these active methods, though, it is highly desired to find a low-maintenance passive flow control method that provides a discernible transition delay in many different operating environments.

Passive flow control methods using phononic crystals (PnC) have been described. PnC can be formed by periodically alternating layers of materials that have a high contrast in mechanical impedance. Bragg diffraction from the PnC architecture, causes the formation of phononic band gaps, which are frequency bands where elastic or acoustic waves cannot propagate within the PnC. Passbands exist between the band gaps, defining frequency bands where waves can propagate. By engineering the material properties and the geometric architecture of the PnC, it is possible to precisely design the dispersion of elastic/acoustic waves. With the ability to control waves, the PnC design was applied to the task of passively controlling wave-like sinusoidal fluctuations in the freestream velocity during subsonic flow, termed Tollmien-Schlichting (T-S) waves. The amplification of T-S waves can cause the boundary layer to transition from laminar to turbulent; so, being able to restrict the amplitude growth of the T-S waves with PnCs holds hope for delaying transition. This problem has been addressed by attaching a PnC to the fluid/solid interface of a channel flow. This arrangement has been termed a phononic subsurface (PSub). As shown in FIG. 1A, fluid can interact with the top surface of the PSub with the PSub attached in a manner that allows it to displace normal to the channel wall. Computationally, it has been demonstrated that a T-S wave incident on a PSub with an engineered band gap and truncation resonance, can cause an out-of-phase displacement of the PSub interface with the T-S wave forcing. The T-S wave forces the surface of the PSub, and the constitutive response of the PSub is to displace upwards into the flow and out-of-phase with the growing T-S wave. This PSub behavior can cause a local reduction in kinetic energy of the T-S wave at the location of the PnC. Because this passive PSub flow control occurs over a single interaction area of the fluid/solid interface, it can be termed a single input, single output (SISO) device. The single input being the T-S wave forcing and the single output the PSub displacement.

It is important to contrast the use of PSubs for flow control to that of a compliant surface, as it provides a significant motivation for further investigation. In 2001, it was reported that the use of compliant surfaces for aeronautical laminar flow control required that the magnitude of the wall and fluid inertias must match to have a significant effect on a T-S wave. To meet these criteria, air vehicle walls would be so weak that conventional forces of flight would destroy them. However, PnCs can be made of rigid materials (ABS plastic and aluminum) which are sufficiently strong to survive the forces of flight. In this work, inspiration is taken from previous research on compliant surfaces, and an approach is developed to realize transition delay with realistic materials for flight. It has been documented in literature that compliant patches can stabilize T-S waves via an irreversible energy transfer to the wall. It has also been shown that increasing the compliance of the wall is better for decreasing the growth of wave instabilities, but increasing the damping in the wall increases growth of the T-S waves. The problem with compliant patches is that the elastodynamic response of the patch is locked by the material properties. The benefit of using PSubs is that they can be composed of an assembly of parts, as shown in FIGS. 1B and 1C, rather than a single or multiple slabs of material, which gives a greater ability to tailor the elastodynamic response and optimize the compliance and damping parameters. When a T-S wave interacts with the surface of a PSub, the PSub responds with a displacement of its interaction surface relative to the T-S wave, which can alter the amplitude of the growing instability. Another benefit of PSubs is that adaptive materials can be utilized to dynamically tune their properties in response to the flow environment; however, this would require active control, and here the sole focus is on passive control.

SUMMARY OF THE INVENTION

The present invention overcomes the foregoing problems and other shortcomings, drawbacks, and challenges described herein. While the invention will be described in connection with certain embodiments, it will be understood that the invention is not limited to these embodiments. To the contrary, this invention includes all alternatives, modifications, and equivalents as may be included within the spirit and scope of the present invention.

According to one embodiment shown and described herein, a multi-input, multi-output phononic system includes: a first interface surface and a second interface surface that respond to at least one of a pressure gradient or a velocity gradient in a wave of a turbulent fluid flow or a laminar fluid flow, the pressure gradient or the velocity gradient associated with complex motion of the flow exhibiting a plurality of frequencies exerted on one or more of the interface surfaces; and a subsurface feature extending from the interface surfaces, the subsurface feature comprising a phononic crystal or locally resonant metamaterial adapted to receive one or more of the pressure gradient or the velocity gradient from the fluid flow via the interface surfaces and to alter one or more of a phase and an amplitude of a plurality of frequency components of the fluid flow.

According to another embodiment shown and described herein, a method of controlling a flow includes providing a first interface surface and a second interface surface in a fluid flow along a flow structure, the fluid flow exhibiting turbulent or laminar flow characteristics, the first interface surface and the second interface surface arranged within the flow structure such that they each are exposed to one or more of a pressure gradient and a velocity gradient associated with complex motion of the fluid flow; and providing a subsurface feature extending from the interface surfaces, the subsurface feature comprising a phononic crystal or locally resonant metamaterial adapted to receive one or more of the pressure gradient or the velocity gradient from the fluid flow via the interface surfaces and to alter one or more of a phase and an amplitude of a plurality of frequency components of the fluid flow; and passively altering one or more of a phase or amplitude of a plurality of frequency components of the flow via the subsurface feature.

According to yet another embodiment, a system for reducing skin friction of a surface using multi-input, multi-output phononic subsurfaces for passive boundary layer transition delay includes: a first interface surface and a second interface surface that respond to at least one of a pressure gradient or a velocity gradient in a wave of a turbulent fluid flow or a laminar fluid flow, the pressure gradient or the velocity gradient associated with complex motion of the fluid flow exhibiting a plurality of frequencies exerted on one or more of the interface surfaces; and a subsurface feature comprising a plurality of masses extending from the interface surfaces, the subsurface feature adapted to receive one or more of the pressure gradient or the velocity gradient from the fluid flow via the interface surfaces and to passively alter one or more of a phase and an amplitude of a plurality of frequency components of the fluid flow.

Additional objects, advantages, and novel features of the invention will be set forth in part in the description which follows, and in part will become apparent to those skilled in the art upon examination of the following or may be learned by practice of the invention. The objects and advantages of the invention may be realized and attained by means of the instrumentalities and combinations particularly pointed out in the appended claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments of the present invention and, together with a general description of the invention given above, and the detailed description of the embodiments given below, serve to explain the principles of the present invention.

In this document, relational terms such as first and second, top and bottom, and the like may be used solely to distinguish one entity or action from another entity or action without necessarily requiring or implying any actual such relationship or order between such entities or actions. The terms “comprises,” “comprising,” “includes,” “including,” “has,” “having,” or any other variations thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. An element preceded by “comprises . . . a” does not, without more constraints, preclude the existence of additional identical elements in the process, method, article, or apparatus that comprises the element.

Reference throughout this document to “one embodiment,” “certain embodiments,” “an embodiment,” “implementation(s),” “aspect(s),” or similar terms means that a particular feature, structure, or characteristic described in connection with the embodiment is included in at least one embodiment of the present invention. Thus, the appearances of such phrases or in various places throughout this specification are not necessarily all referring to the same embodiment. Furthermore, the particular features, structures, or characteristics may be combined in any suitable manner in one or more embodiments without limitation.

The term “or” as used herein is to be interpreted as an inclusive or meaning any one or any combination. Therefore, “A, B or C” means “any of the following: A; B; C; A and B; A and C; B and C; A, B and C.” An exception to this definition will occur only when a combination of elements, functions, steps or acts are in some way inherently mutually exclusive. Also, grammatical conjunctions are intended to express any and all disjunctive and conjunctive combinations of conjoined clauses, sentences, words, and the like, unless otherwise stated or clear from the context. Thus, the term “or” should generally be understood to mean “and/or” and so forth.

All documents mentioned herein are hereby incorporated by reference in their entirety. References to items in the singular should be understood to include items in the plural, and vice versa, unless explicitly stated otherwise or clear from the text.

Recitation of ranges of values herein are not intended to be limiting, referring instead individually to any and all values falling within the range, unless otherwise indicated, and each separate value within such a range is incorporated into the specification as if it were individually recited herein. The words “about,” “approximately,” or the like, when accompanying a numerical value, are to be construed as indicating a deviation as would be appreciated by one of ordinary skill in the art to operate satisfactorily for an intended purpose. Ranges of values and/or numeric values are provided herein as examples only, and do not constitute a limitation on the scope of the described embodiments. The use of any and all examples, or exemplary language (“e.g.,” “such as,” or the like) provided herein, is intended merely to better illuminate the embodiments and does not pose a limitation on the scope of the embodiments. No language in the specification should be construed as indicating any unclaimed element as essential to the practice of the embodiments.

For simplicity and clarity of illustration, reference numerals may be repeated among the figures to indicate corresponding or analogous elements. Numerous details are set forth to provide an understanding of the embodiments described herein. The embodiments may be practiced without these details. In other instances, well-known methods, procedures, and components have not been described in detail to avoid obscuring the embodiments described. The description is not to be considered as limited to the scope of the embodiments described herein.

In the following description, it is understood that terms such as “first,” “second,” “top,” “bottom,” “up,” “down,” “above,” “below,” and the like, are words of convenience and are not to be construed as limiting terms. Also, the terms apparatus and device may be used interchangeably in this text.

FIG. 1A presents a bimaterial rod-based SISO of the prior art, according to one or more embodiments shown or described herein.

FIG. 1B presents a coiled mass/beam-based SISO of the prior art.

FIG. 1C presents an extended mass/beam-based MIMO.

FIG. 2 presents skin friction coefficient at several generator amplitudes.

FIG. 3A presents a neutral curve for Blasius boundary layer in the V Rex-F plane.

FIG. 3B presents kinetic energy growth downstream of a generator (2D simulation, B_g/L=2×10⁻⁵).

FIG. 4 presents a configuration for active control simulations.

FIG. 5 presents ΔK vs. phase for two controller segments, B_c=5 μm

FIG. 6 presents spatial growth rate o as a function of Y^†.

FIG. 7A illustrates how the first and third masses of the MIMO PSub react to original T-S wave pressure field, coupled through upstream (−) and downstream (+) interaction surfaces by displacing, thereby inducing additional T-S waves.

FIG. 7B illustrates that the induced T-S waves are dependent on the dispersion curves.

FIG. 7C illustrates that the induced T-S waves are dependent on the component receptances.

FIG. 7D illustrates a frequency of the T-S waves.

FIG. 8A illustrates an upstream apparent collocated receptance.

FIG. 8B illustrates a downstream apparent collocated receptance. These receptances may determine if the induced T-S waves will constructively or destructively interfere with the original T-S wave. Horizontal dashed lines in the phase panels are the boundaries of stability from Table 4 determined by CFD.

FIG. 8C illustrates a zoomed-in view of the phase panels.

FIG. 8D illustrates a zoomed-in view of the phase panels. These figures emphasize that stable frequency bands of phase are realized. Note, at 250 Hz and k_TSL_s=90°, phases are ∠H−=17.967° and ∠H+=−21.913°.

FIG. 9 presents a representative configuration for passive control simulations.

FIG. 10A presents amplitude for FSI surfaces alongside equivalent 2D simulation.

FIG. 10B presents phase for FSI surfaces alongside equivalent 2D simulation.

FIG. 11A presents an integrated TKE profile near PSub surfaces.

FIG. 11B presents an integrated TKE profile near skin friction profile, C_f.

FIG. 12 presents a three-dimensional flow structure for rigid and controlled cases using an iso-surface of Q-criterion (Q=50).

FIG. 13 is a first illustration of a PSub system according to one or more embodiments shown and described herein.

FIG. 14 is a second illustration of the PSub system of FIG. 13.

FIG. 15 shows an additional view of the PSub system of FIG. 13.

FIG. 16A shows the PSub system of FIG. 13 in a flow surface system with a cover plate installed.

FIG. 16B shows the flow surface system of FIG. 15A with the cover plate removed, showing the PSub system of FIG. 13.

It should be understood that the appended drawings are not necessarily to scale and may represent a simplified depiction of various features, which may be illustrative of the basic principles of the invention. The specific design features and/or the sequence of operations as disclosed herein, including, for example, specific dimensions, orientations, locations, and shapes of various illustrated components, will be determined in part by the particular intended application and use environment. Certain features of the illustrated embodiments have been enlarged or distorted relative to others to facilitate visualization and clear understanding. In particular, thin features may be thickened, for example, for clarity or illustration.

DETAILED DESCRIPTION OF THE INVENTION

Phononic subsurfaces (Psub) may be used to develop a passive boundary layer transition delay strategy. High-order implicit large-eddy simulation (ILES) may be used to explore control requirements for boundary layer transition delay in the context of guiding passive flow control strategies. Positive phasing of surface displacement relative to forcing from a Tollmien-Schlichting (T-S) wave is shown to be effective at attenuating instability. Although difficult to achieve with a single input system, positive phasing can be realized with a multi-input, multi-output (MIMO) system. Positive phasing, and therefore a small degree of passive transition delay, may be demonstrated on a flat plate boundary layer. Controller efficacy may be improved with further optimization of subsurface properties and placement of multiple devices.

Nomenclature

• • A=interaction area
  • b=width of PSub beam
  • B=amplitude of surface motion
  • c=wavespeed
  • C=damping matrix
  • C_f=skin friction coefficient
  • C_p=pressure coefficient, C_p=(p−p_∞)/½ ρU²_∞
  • D=dynamic stiffness matrix
  • d=displacement vector
  • d_i=displacement vector element
  • E=Young's modulus of PSub beam
  • f=frequency
  • f=loading vector
  • f_j=loading vector element
  • F=non-dimensional frequency, F=(2πfv/U²_∞)×10⁶
  • h=thickness of PSub beam
  • H=receptance matrix
  • H_ij=receptance matrix element
  • I=identity matrix
  • k=wavenumber in the flow
  • K=turbulent kinetic energy integrated across the boundary layer
  • K=stiffness matrix
  • l=length of PSub beam
  • L=reference length
  • L_c=controller length
  • L_g=generator length
  • L_s=distance between upstream and downstream interaction surfaces, λ_TS (n+γ)
  • L=primal assembly matrix
  • m=concentrated rigid mass
  • M=Mach number
  • M=moment (finite element model)
  • M=mass matrix
  • n=parameter of Ls specifying family of equivalent coupling locations, {0, 1, 2, . . . }
  • p=pressure
  • P=mechanical power
  • Re_L=Reynolds number based on reference length, Re_L=ρ_∞U_∞L/μ_∞
  • t=time
  • T=Floquet-Bloch condition matrix
  • S=spanwise extent
  • S=shear force (finite element model)
  • U_∞=free stream reference velocity
  • u, v, w=Cartesian components of flow velocity
  • v=displacement of PSub (finite element model-perpendicular to x-direction)
  • x, y, z=coordinate system aligned with streamwise, normal, and spanwise directions
  • Y^†=interaction area-receptance product, Y^†=AH*
  • α=beam element half length
  • γ=parameter of L_s controlling forcing phase shift
  • γ=travel time of T-S wave from upstream to downstream interaction surfaces, L_s/c_TS
  • Δτ=non-dimensional time step
  • ι=imaginary unit, √−1
  • κ=wavenumber in the PSub
  • λ=wavelength
  • μ=dynamic viscosity
  • ρ=fluid density
  • σ=spatial growth rate of T-S wave
  • τ=non-dimensional time, τ=tU_∞/c
  • υ=density of beam
  • ϕ-phase of complex number in exponential form
  • Φ_c=phase of controller motion relative to generator
  • ω=angular frequency, ω=2πf

Subscripts

• • ω=free stream reference value
  • −=upstream PSub interaction location
  • +=downstream PSub interaction location
  • b=beam element
  • {tilde over (b)}=disconnected arrangement of beam elements
  • b=assembled system of beam elements
  • c=parameter associated with controller
  • f^˜=disconnected arrangement of unit cells (finite PnC)
  • g=parameter associated with generator
  • i=index of displacement vector elements, also active control surface index
  • j=index of loading vector elements
  • PSub=parameter associated with PSub
  • r=real part of complex number
  • TS=parameter associated with T-S wave
  • w=value at the wall
  • ũc=disconnected arrangement of masses and assembled beam element systems (unit cell)
  • uc=assembled system of mass and beam elements (unit cell)
  • ι=imaginary part of complex number

Operators

• • |⋅|=magnitude of complex number
  • ∠=phase of complex number
  • F⁻¹{⋅}=Inverse Fourier transform

2.1. Fluid Dynamics Solver

Equations and other problems associated with the fluid dynamics described herein may be solved using a fluid dynamics solver (e.g., an aerodynamics solver). The solver used may be, for example, a solver such as the FDL3DI. Such code may solve full sets of compressible Navier-Stokes equations. The solver may be a high-order solver and employ, for example, a sixth-order compact difference scheme to compute spatial derivatives. Simulations may be integrated in time using a second-order implicit approximately factorized and diagonalized scheme. Time-integration may be augmented using a Newton-like sub-iteration procedure to drive down residual error and maintain second-order temporal accuracy.

In some embodiments, an implicit large-eddy simulation (ILES) technique may be implemented enabling solution of transitional and turbulent flows. This approach can use a high-order low-pass spatial filter applied to the flow field solution after sub-iterations (e.g., each sub-iteration). The filter can selectively damp the highest wavenumber content not supported by the grid resolution in place of subgrid-scale and heat flux models used in standard large-eddy simulation techniques. Such an ILES implementation can be an effective alternative to subgrid-scale models. As grid resolution increases or Reynolds number decreases, such an ILES technique can approach direct numerical simulation.

2.2. Fluid-Structure Interaction

Implicit coupling between, for example, a fluid dynamics solver such as the FDL3DI and the MIMO PSub model can be accomplished via a sub-iteration procedure. Within each sub-iteration, surface stresses can be integrated over each PSub interaction surface and transferred to the corresponding PSub nodes. The surface pressure can be balanced by a cavity pressure (i.e., Δp=p−p_cavity) where the cavity pressure, p_cavity, can be set to the time-mean surface pressure at each PSub interface in order to maintain a zero mean deflection. The resulting surface displacements and velocities are returned to the aerodynamics solver and used to adjust the boundary condition and grid motion at each PSub interface. This interchange may be repeated within each time step thereby synchronizing the two procedures and preserving second-order temporal accuracy for both the fluid and structural dynamics models. It should be noted that while the solver (e.g., aerodynamics solver) may be non-dimensional, the PSub model can be dimensional. The PSub structure is represented by a finite element model (FEM) model composed of Euler-Bernoulli beam elements and rigid masses, simulated in the time domain by one or more simulation methods, for example, the Newmark method. Forces and displacements can be translated between their dimensional/non-dimensional counterparts during this data exchange. Grid deformation can be accommodated using a cubic blending polynomial over a small region of influence above the plate surface.

3. CONFIGURATION

Various configurations are possible to achieve the results described herein. One exemplary configuration may be based on, for example, the use of piezoelectrically-driven oscillating disc actuators for generating disturbances and for subsequently controlling the growth of those disturbances. Some configurations may use, for example, a thick plate in order to house a subsurface in various wind tunnel experiments. In some embodiments, a 6:1 leading-edge aspect ratio may be used to minimize a risk of separation (e.g., the risk of separation that may be present with thicker geometry). In some embodiments, the configuration may be two-dimensional and exemplary relevant parameters are described in Table 1.

TABLE 1
Problem Setup
	Property	Value
	Freestream velocity, U∞	16	m/s
	Freestream air density, ρ∞	1.225	kg/m3
	Mach number, M∞	0.04665
	Reynolds number,	1,081,642
	ReL = p∞V∞L/μ∞
	Plate thickness	0.0508	m
	Leading edge aspect ratio	6:1
	Reference length, L	1	m
	Generator length, Lg	0.00635	m
	Generator amplitude, Bg	1 × 10−5	m
	Generator frequency, fg	250	Hz
	Generator position, xg	0.273	m
	Spanwise extent, s/L	0.05

Table 1—Problem Setup

In one particular embodiment, T-S waves were generated by an oscillating surface positioned at x_g=0.273 m

$\begin{array}{cc}{d}_g=\{\begin{array}{cc}{B_g\left[1-{\left(\frac{x-x_g}{L_g/2}\right)}^2\right]}^2\sin \left(2\pi \frac{f_{g}L}{U_{\infty }}\tau \right)& 2❘x-x_g❘\le L_g\\ 0& 2❘x-x_g❘>L_g\end{array},& \mathrm{Eq}.\left(1\right)\end{array}$

where B_g is the generator amplitude. In this exemplary embodiment, the generator was set to oscillate at a single dimensional frequency of f_g=250 Hz.

In embodiments, the computational mesh can be created using a hyperbolic extrusion from the plate surface. A half-plate model can be created by replacing the lower half with a symmetry plane. Dimensions for the resulting grid are recorded in Table 2.

Maximum values for grid spacing in wall units, Δx+, Δy+, Δz+, were obtained from spanwise periodic simulations of a generator-only case, (B_g/L=1×10⁻⁵), and provided in Table 2. These values occur at the peak in C_f shown in FIG. 2. All values were well within the bounds recommended for wall-resolved LES. Previous experiments have demonstrated suitability of this mesh resolution for the turbulent regions of flow on nearly the same configuration.

The laminar flow regions, of specific importance to embodiments based on the principles of the technologies described herein, are particularly well resolved. In one embodiment, the grid can give 46 points across the λ_TS=0.023 m wavelength of the T-S wave. In this embodiment, there can be 112 and 137 grid points across the boundary layer at the generator and controller locations, respectively.

In one exemplary embodiment, boundary conditions may be prescribed to the computational mesh as follows. The no-slip condition can be enforced at the plate surface. Surface position and velocity at the generator may be directly applied according to Eq. (1) and controller surface displacement/velocity may be imposed according to subsurface motion. Third-order adiabatic and vanishing normal pressure gradient conditions may also be imposed at the surface. Freestream conditions may be specified at the far field boundary located approximately 50 L away from the plate surface. First-order extrapolation of the primitive variables was applied at the outflow boundary. Three-dimensional simulations can be completed using periodic conditions at the spanwise boundaries. These boundaries can include a 9-point overlap to maintain the high-order stencil of the compact scheme and implicit filter. The domain width of the spanwise-periodic simulations may be s/L=0.05 in one or more embodiments.

As noted previously, the lower half of the plate may be replaced by a symmetry plane. Symmetry can be imposed on the grid line extending outward from the leading edge. Although not shown here, full plate simulations can be tested in two-dimensional cases demonstrating the symmetry assumption may not affect the upper surface flow behavior.

TABLE 2
Grid Parameters
Property        Value
Grid dimensions 2150 × 272 × 209
Grid size       122,223,200
Δx/L            5.0 × 10−4
Δyw/L           2.0 × 10−5
Δz/L            2.5 × 10−4
Δτ              1.6 × 10−5
Δxmax+          25.0
Δyw, max+       1.0
Δzmax+          12.5

In one or more embodiments, two-dimensional cases may be run with a very small time step of Δτ=3.2×10⁻⁵ (Δτ=2.0×10⁻⁶ s) which may allow for 2000 steps per cycle at a disturbance frequency of f_g=250 Hz (F=90.7).

Spanwise-periodic simulations can require a smaller time step of Δτ=1.6×10⁻⁵ (Δτ=1×10⁻⁶ s) to maintain stability in the fine grid region at the leading edge of the plate. All simulations were initiated from a previously converged unperturbed boundary layer where the generator and PSub/controller were enabled simultaneously. Simulations were then run until a time-asymptotic state was achieved before time-mean and statistical data was collected.

In various embodiments, several generator amplitudes were considered for testing with the PSub. Skin friction profiles for spanwise-periodic simulations for each generator amplitude are shown in FIG. 2. The gray bands in this figure indicate the eventual location of passive control surfaces for reference. The lowest generator amplitude, B_g=10 μm was selected for all cases in this study to place boundary layer transition (as indicated by a rapid increase in C_f) sufficiently far downstream of the planned control surfaces and to ensure that the PSub interacts with instabilities within the linear regime.

This simulation is referred to as the rigid case and used for comparison in each of the following sections. Theoretical predictions for laminar and turbulent boundary layer C_f are provided for comparison.

FIG. 3(a) shows the neutral curve for the Blasius boundary layer. The location of the generator, x_g, and passive controller surfaces, x− and x+, are shown for the disturbance frequency of f_g=250 Hz. Stability theory predicts the flat plate boundary layer is unstable to the chosen frequency in the region of interest and this is indeed supported by FIG. 3B which shows the growth of TKE integrated across the boundary layer downstream of the generator. Rapid growth of the 250 Hz instability is maintained for the region of interest.

4. SUBSURFACE DESIGN REQUIREMENTS USING ACTIVE CONTROL

Active control using small-amplitude prescribed motion of multiple adjacent control surfaces, as depicted in FIG. 4, was explored first in order to better understand design requirements for passive T-S wave mitigation using distributed subsurfaces. Effects of various combinations of phase, amplitude, and number of controller segments on T-S wave amplification/attenuation were explored.

Relevant properties of this configuration are provided in Table 3. While a total of four control surfaces have been considered, the discussion will primarily focus on the case of two control surfaces.

TABLE 3
Parameters for active control simulations
	Property	Value
	Generator location, xg	0.273	m
	Generator amplitude, Bg	1 × 10−5	m
	Generator frequency, fg	250	Hz
	Controller length, Lc	0.012	m
	Controller amplitude, Bc	5	μm
	Control surface 1 location, x1	0.388	m
	Control surface 2 location, x2	0.400	m
	TS wavelength, λTS	0.023	m

In one or more embodiments, the first control surface can be centered at x₁=0.388 m with additional surfaces placed adjacently at intervals, for example, at intervals of 0.012 m downstream. In one or more embodiments, each control surface can be flat with a streamwise length of L_c=0.012 m. This is roughly half the T-S wavelength and may maximize the surface area exposed to the T-S wave for the passive simulations otherwise described herein. Although not shown, the results of this section may not be qualitatively affected if the controller length is reduced. Some blending with adjacent surfaces can be applied at the edges of each controller to preserve grid quality. In some embodiments, surface motion may be defined according to a traveling wave equation, given as

$\begin{array}{cc}{d}_i=B_c\sin \left(2\pi x_i/{\lambda }_{\mathrm{TS}}-2\pi f_{g}t+{\Phi }_c\right)& \mathrm{Eq}.\left(2\right)\end{array}$

where d_i is the vertical motion (displacement) of control surface i, x_i is the controller position, λ_TS=0.023 m is the measured T-S wavelength, f_g=250 Hz is the generator/T-S wave frequency, and Φ_c sets the phase of motion. Matching the parameters to those of the T-S wave can allow for consistent phasing of each control surface relative to the T-S wave. Changing the value of Pc can shift the phase of surface motion relative to the T-S wave. The actual phase of each surface displacement relative to the T-S wave, ϕ(d_i/−p_i), for each Φ_c can be measured from simulation data.

Phase can be measured using the ratio of Hilbert transform of the time-domain signals of displacement, di, and negative surface pressure measured at the center of the controller. This definition is representative of the relationship between displacement and loading on a passive surface, i.e., f_i≈−p_iA. Note that phasing can be found to be nearly identical across all control surfaces when two or more controllers are active.

In one or more embodiments, one, two, and four controller segments were tested at amplitudes of B_c=1, 5, and 10 μm at intervals of ΔΦ_c=20°. Only the results for the two-segment case at B_c=5 μm are shown here for brevity. Two quantities were of interest; the phase of surface motion relative to the pressure force on the surface due to the T-S wave, ϕ(d_i/−p_i) and the effect of this phasing on the T-S wave downstream. For the computational fluid dynamics (CFD) simulations, this was quantified by integrating turbulence kinetic energy (TKE) across the boundary layer along coordinate lines normal to the plate surface defined as,

$\begin{array}{cc}K\left(x\right)=\frac{1}{{\mathrm{LU}}_{\infty }^2}\int \frac{1}{2}\left(\stackrel{\_}{u^{\prime }{u}^{\prime }}+\stackrel{\_}{v^{\prime }{v}^{\prime }}+\stackrel{\_}{w^{\prime }{w}^{\prime }}\right)\mathrm{dy}& \mathrm{Eq}.\left(3\right)\end{array}$

where u′, v′, w′ are the fluctuations of velocity components. FIG. 5 shows the difference between the controlled and uncontrolled cases, ΔK=K_control−K_rigid, at a single position downstream of the control surfaces, x=0.5 m. The horizontal axis portrays the phasing of the surface motion relative to negative surface pressure, ϕ(d/−p). In general, these results show that a positive phase is desirable for attenuating T-S waves with peak attenuation near ϕ≈57°.

This observation was further supported by linear stability analysis. A linear framework was developed using the classical Orr-Sommerfeld equation where the moving/compliant surfaces are modelled using a dynamic boundary condition in the form of complex admittance, Y(x)={circumflex over (v)}/{circumflex over (p)}.

The Orr-Sommerfeld equation is then recast as a general eigenvalue problem to attain the spatial growth rate σ. In this case, the effect of forced surface motion is introduced through the parameter Y.

FIG. 6 shows the spatial growth rate as a function of

$\begin{array}{cc}{Y}^{†}=\frac{\hat{d}}{-\hat{p}}=❘Y^{†}❘e^{-{\mathrm{\iota \angle Y}}^{†}}& \mathrm{Eq}.\left(4\right)\end{array}$

which is simply an altered version of the wall admittance in the context of displacement and negative surface pressure for consistency with the present definition of phasing. The linear stability map predicts T-S wave attenuation for generally positive phases in qualitative agreement with the CFD simulations. Bounds for T-S wave attenuation as well as optimal phase values for each are provided in Table 4 which are again very close between the two approaches. The optimal from the linear stability identifies the phasing that provides complete stabilization at the smallest amplitude whereas ϕ_optimal from the CFD is the value that minimizes K at a fixed amplitude. Note, the negative sign appearing in the argument of Eq. (4) accounts for a sign difference in the definition of the traveling wave phase function used in the stability analysis (−ιωt), and the one used in the definition of the inverse Fourier transform (+ωt) of the CFD generated time-domain signal analysis, as well as the analysis of the PSub to follow. Note, ι is the imaginary unit, √−1.

TABLE 4
Phase bounds for T-S wave attenuation
	ϕmin	ϕoptimal	ϕmax
	Linear stability	−10°	90°	170°
	CFD	−20°	57°	161°

Although not shown, in some embodiments, this analysis was repeated for several amplitudes, frequencies, number of controller segments, and Reynolds numbers and found to be broadly consistent across all parameters considered in both the linear stability and CFD frameworks. Similar results were also found for a small amplitude traveling surface wave. The greatest sensitivity was observed for the single controller case at larger amplitudes in the CFD simulations (linear stability was unaffected by the number of controller segments). This behavior seems related to an upstream effect due to acoustic emissions from the control surface. This effect appeared to be alleviated when multiple adjacent control surfaces were operating out of phase thereby reconciling with the theoretical results.

These active control simulations can show that distributed micrometer-scale surface motion can be effective for boundary layer transition delay on a flat plate. The critical factor appears to be the phasing of surface motion relative to the T-S wave, ϕ(d/−p). This observation is applied to the design of a subsurface capable of providing a passive stabilizing response as described herein.

5. PASSIVE CONTROL

Observations of the active control simulations were used to design and implement a new passive subsurface control strategy. Specifically, the phase ranges listed in Table 4 are used in one or more embodiments as stability criteria of the MIMO Psub response surfaces. Theoretical development is described herein, as well as implementation of fully-coupled fluid structure interaction (FSI) CFD simulations demonstrating successful boundary layer transition delay.

5.1. Multi-Input Multi-Output Phononic Subsurface

Psubs are a class of PnC, which are embedded in a flow-bounded material/structure and allowed to interact through discrete coupling interfaces that transduce pressure/displacements between the Psub and exterior flow phenomenon. When a Psub interacts with a flow containing a T-S wave, the vibration of the Psub can cause a localized reduction in the TKE of the flow in the volume above the PSub interaction surface. Designs of some Psubs rely on displacement to force phasing of −180°, provided by a truncation resonance. Since the PSub displacement-over-forcing response phase is stable (i.e., the value changes slowly over a long range of frequency) inside a band gap, the PSub acts a robust passive actuator. The approximately sinusoidal T-S wave pressure imparts a loading force on the interaction surface and causes the PSub to respond by vibrating according to its dynamic compliance (i.e., receptance). The displacement of the PSub, ultimately causes the generation of additional time varying pressure fluctuations in the flow that can lead to a localized cancellation of the T-S wave, indicated by a drop in the TKE of the flow.

Otherwise discussed herein are investigations of active control via forced surface motion using linear stability analysis as well as nonlinear CFD simulations. Such investigation may demonstrate that positive phasing between displacement and force at the interaction surface can cause global flow quieting (i.e., flow stabilization at points after the interaction surface). Specifically, +90° is identified by stability analysis as being ideal, but there is an envelope that includes almost all positive and some slightly negative phases, as being stabilizing, as identified in Table 4. However, positive phasing may be unfeasible from a passive vibrating system with only a single point of interaction (i.e., single-input single-output (SISO) system) with the flow.

A general discretized vibrating system can be represented, in the frequency domain, by

$\begin{array}{cc}\widehat{f}=\left[K+\mathrm{\iota \omega }C-{\omega }^{2}M\right]\hat{d}=D\left(\omega \right)\hat{d}& \mathrm{Eq}.\left(5\right)\end{array}$

where M, C, K, D are the mass, damping, stiffness, and dynamic stiffness matrices, {circumflex over (d)} and {dot over (f)} are the nodal displacement and forcing vectors, and ω=2πf is the angular frequency while f is frequency in Hz. The system may be solved for the displacement due to a given load for each frequency by D⁻¹ {circumflex over (f)}=H {circumflex over (f)}={circumflex over (d)}, where His the receptance matrix, whose elements H_ij(ω) relate the displacement {dot over (d)}_i (an element of {circumflex over (d)}) to the applied force {circumflex over (f)}_i (an element of {circumflex over (f)}).

Given that the forced motion study identified an envelope of phases that lead to T-S wave stabilization, it can be beneficial that a PSub be designed with equivalent SISO PSub response phases within the bounds of Table 4. Analysis discussed herein demonstrates that such response phases can be achieved by allowing multiple points on the structure to interact with the flow. In the simplest case, two points of interaction are considered so that f has two nonzero elements at the upstream (US), indicated by subscript-, and the downstream (DS), indicated by subscript +, locations.

Thus, i, j∈{−,+}, and the displacement at each interaction surface is a superposition, written as

$\begin{array}{cc}{\hat{d}}_{-}=H_{--}{\hat{f}}_{-}+H_{-+}{\hat{f}}_{+}& \mathrm{Eq}.\left(6\right)\end{array}$ $\begin{array}{cc}{\hat{d}}_{+}=H_{+-}{\hat{f}}_{-}+H_{++}{\hat{f}}_{+}& \mathrm{Eq}.\left(7\right)\end{array}$

The response at each point is governed by its collocated receptance, and a response that is propagated from each point to the other by a cross receptance. Forcing at each interaction location is provided by the incident T-S wave pressure, which can be modeled as

$\begin{array}{cc}{p}_{-}=p_{\mathrm{TS}}\left(0,t\right)={ℱ}^{-1}\left\{\hat{p}\left(\omega \right)\right\}& \mathrm{Eq}.\left(8\right)\end{array}$ $\begin{array}{cc}{p}_{+}=p_{\mathrm{TS}}\left(L_s,t\right)={ℱ}^{-1}\left\{\hat{p}\left(\omega \right)e^{-\mathrm{\iota \omega \delta }}\right\}={ℱ}^{-1}\left\{\hat{p}\left(\omega \right)e^{-\mathrm{\iota \omega }{L}_s/c_{\mathrm{TS}}}\right\}={ℱ}^{-1}\left\{\hat{p}\left(\omega \right)e^{-\iota k_{\mathrm{TS}}{L}_s}\right\}& \mathrm{Eq}.\left(9\right)\end{array}$

where p_TS is the time-domain perturbation pressure associated with the T-S wave and k_TS=2 π/λ_TS is the T-S wavenumber. p+ is equivalent to the time delayed signal observed at p−, as is shown in FIG. 7A.

Allowing the two points of flow interaction with the MIMO PSub to occur at different positions along the x direction means that there is a well-understood phase shift between the forcing at US and DS, that corresponds to the time delay δ=L_s/c_TS between the two signals, where T-S wavespeed is c_TS=λ_TSf_g and the T-S center frequency (i.e., generator frequency) is f_g=250 Hz.

It should be understood that the PSub will respond to the bandwidth of frequencies that are imposed by the T-S wave forcing, whose center frequency is f_g, but in order to study the response at other frequencies, a general frequency f is also defined that is related to the angular frequency as ω=2πf. When studying frequencies away from f_g, λ_TS=ω/c_TS. Here, c_TS is assumed to be constant over the studied range of frequencies.

Setting the distance between US and DS positions as x+−x−=L_s, the forces transduced to the structure are {circumflex over (f)}=−A^{{tilde over (p)}}(ω) and {circumflex over (f)}=−A{tilde over (p)}(ω), where A is the coupling area at each interaction surface, and thus the phase between the pair of forces can be adjusted by changing the distance between them in x (i.e., L_s). The length between interaction surfaces is specified as L_s=λ_TS(n+γ)∀n∈{0, 1, 2, . . . } and 0<γ<½, which specifies a family of points with a phase shift range (as a function of γ) of 0<k_TSL_s<180° between US and DS forcing. In the design of FIG. 7A, n=1 and γ=¼ so that k_TSL_s=90°, but a family of points (i.e., those L_s with other n values) is allowed since it may be necessary to include multiple wavelengths between the US and DS to accommodate the length of the MIMO PSub. Limits on y are exclusive, since at the limit bounds, the ACR will lose the large frequency ranges of nearly constant stabilizing phase associated with the band gap. Note, modeling the T-S wave as a propagating wave does not include its growth.

Multiple structures could be designed to produce the required response phasing, but PnCs have been used in the past and have well understood behavior which simplifies the design process. In this work, the MIMO PSub design consists of rigid masses, whose rotations are locked, linked by flexible beams, as shown in FIG. 7A below the Flow-PSub interface line. Masses along the length of the MIMO PSub alternate in value between m₁ and m₂ (i.e., a diatomic beam-based PnC). Parameters of the MIMO PSub model for some embodiments are given in Table 5. PnCs are periodic structures which, exhibit pass bands (i.e., white shaded background regions where dispersion curves exist) and band gaps (i.e., gray shaded regions with no dispersion curves) as shown in FIG. 7B. The finite PSub system resonances, identified as peaks in the component collocated and cross receptances shown in FIG. 7C, are usually confined to pass bands, while only truncation resonances can exist in the band gap. Since, the dispersion curves of the PSub are designed by specifying elements of the unit cell (i.e., the repeated unit used to create the periodic structure of the PnC), the pass band and band gap widths can be specified prior to dealing with resonances of the finite structure. Most importantly, since there are few if any resonances in the band gap it is assured to be a frequency band of stable phase. Note, otherwise described herein is at least one procedure for obtaining the dynamic stiffness matrix of the finite MIMO PSub which can be used to construct the receptance matrix. Also described herein is a process for obtaining dispersion curves, shown in FIG. 7B, from the unit cell dynamic stiffness matrices discussed herein.

The apparent collocated receptance (ACR) at US and DS interaction surfaces can be written:

$\begin{array}{cc}{H}_{-}=\frac{{\hat{d}}_{-}}{{\hat{f}}_{-}}=H_{--}+H_{-+}\frac{{\hat{f}}_{+}}{{\hat{f}}_{-}}=H_{--}+H_{-+}{e}^{-\iota k_{\mathrm{TS}}{L}_s}& \mathrm{Eq}.\left(10\right)\end{array}$ $\begin{array}{cc}{H}_{+}=\frac{{\hat{d}}_{+}}{{\hat{f}}_{+}}=H_{+-}\frac{{\hat{f}}_{-}}{{\hat{f}}_{+}}+H_{++}=H_{+-}{e}^{-\iota k_{\mathrm{TS}}{L}_s}{H}_{++}& \mathrm{Eq}.\left(11\right)\end{array}$

Equations (10) and (11) are the effective responses when the component collocated receptances, shown in FIG. 7C, are combined taking into account the 90° forcing phase shift between US and DS positions (i.e. L_s=5/4λ_TS=(n+γ)λ_TS). If γ had been zero or ½, then k_TSL_s would be 0° or 180°, which would stop the ACR from taking on phase values other than 0° or 180°. Consequently, as long as the phase of the ACR about the T-S wave frequency is within the stability bounds established by the forced motion studies for a single surface, it is expected that the induced motion of the MIMO PSub coupled surfaces will generate canceling T-S waves as in FIG. 7A, where the original T-S wave (red) induces T-S waves 1 (magenta) and 2 (aqua) with phase shifts dictated by the ACR. The ACR computed at the US and DS locations are plotted for multiple values of the forcing phase shift, k_TSL_s, in FIGS. 8A and 8B, respectively. Also, zoomed-in views of the ACR phase curves are shown in FIGS. 8C and 8D to demonstrate that the phase of the ACR is stable in the band gap and falls approximately within the stabilizing range. At US and DS, the ACR phases are ∠H−=17.967° and/H+=−21.913° for 250 Hz with k_TSL_s=90°, respectively, which are consistent with the stability bounds found by CFD in Table 4. These results can be related to the stability analysis described herein by realizing that the ACR at US or DS can be written as H={circumflex over (d)}/{circumflex over (f)}={circumflex over (d)}/(−{circumflex over (p)}A)=(Y^†)*/A. Note, the complex conjugation corrects for the negative in the argument of Eq. (4). It is important to note that the US surface reduction in TKE, dominates the effect of the DS surface which has a marginal impact on TKE, and can be understood by review of FIG. 5.

Finally, a method has been presented for realizing robust T-S wave suppression using a MIMO PSub that is designed with ACR (at each interaction surface) with broad frequency ranges of phase within the envelope of stabilization (as specified by the previous analysis). The main trade off here is that there is typically low receptance amplitude inside the band gap, but in the presented MIMO PSub, the ACR amplitudes slowly increase across the band gap as frequency approaches the second pass band, as shown in FIGS. 8A and 8B, which results in a sufficient amount of amplitude to create a meaningful T-S wave cancellation.

5.2. Computational Demonstration of Transition Delay

The MIMO PSub concept was incorporated into the flat plate geometry as shown in FIG. 9. Similar to the active control simulations, in one or more embodiments, the length of each control surface was L_c=0.012 m (≈λ_TS/2) with some smoothing at the edges to preserve grid quality. Control surfaces were placed at x−=0.388 m and x+=0.4165 m. The distance between x+ and x− was determined by the subsurface design which specified a separation of (5/4) λ_TS to achieve proper phasing of control surface motion, as described herein. The interaction surfaces are FEM nodes 7 and 21, which correspond to the first and third masses of the MIMO Psub shown in FIG. 7A. Forces on each control surface were obtained by integrating pressure over the relevant surface points. As noted herein, in one or more embodiments, three-dimensional simulations were completed using spanwise-periodic conditions. Therefore, forces were computed per unit span and then scaled to the specified width of the interaction surface.

TABLE 5
Parameters of the MIMO PSub
	Property	Value
	Density, v	2700	kg/m3
	Young's Modulus, E	68.9	GPa
	Beam thickness, h	0.35	mm
	Beam width, b	1	mm
	Beam length, l	14.375	mm
	Mass 1, m1	4	g
	Mass 2, m2	0.72	g

Passive control using Psubs may best be achieved through a distributed array of numerous subsurfaces, each providing small contributions to boundary layer transition delay with a large effect in aggregate. The effect of a single MIMO Psub device is shown here. Therefore, a fairly large width of s=0.4 m was used in one or more embodiments to produce a strong response with a clearly demonstrated effect on boundary layer transition. Further optimization of the subsurface may drive down surface area requirements alongside application of subsurface arrays.

TABLE 6
Parameters for passive control simulations
	Property	Value
	Generator location, xg	0.273	m
	Generator amplitude, Bg	1 × 10−5	m
	Generator frequency, fg	250	Hz
	Controller length, Lc	0.012	m
	Controller width, wc	0.4	m
	Upstream interaction surface, x−	0.388	m
	Downstream interaction surface, x+	0.4165	m
	TS wavelength, λTS	0.023	m

In one or more embodiments, passive control simulations were initiated from a previously converged boundary layer plate solution. Generator and controller were enabled at the same time and simulations were run until achieving a time-asymptotic state. This required a run time of τ=9.6 or 600,000 iterations at the chosen time step. Time-mean and statistical data was then collected over three runs of τ=3.2 (200,000 iterations) for a total run time of τ=19.2 (1,200,000 iterations). Quantities of interest were then compared over the three consecutive runs to verify the ILES settled on a consistent response as demonstrated herein. MIMO PSubs may be tested by running computationally efficient two-dimensional simulations, which may provide a good prediction of the flow response up to onset of three dimensional instabilities.

Time-mean amplitude (amplitude and phase were extracted from the fast Fourier transform of displacement and load in this case.), |d|, is shown in FIG. 10A for each of the PSub control surfaces. In embodiments, amplitude may not change significantly between consecutive runs and may be consistent with those predicted by purely two-dimensional simulations. Phasing of the displacement relative to the input load is shown in FIG. 10B for both control surfaces. The shaded region indicates phasing predicted by linear stability to provide attenuation, see FIG. 6. In one or more embodiments, the time-mean phase of both surfaces was found to be consistent between consecutive runs and effectively identical to the response of the purely two-dimensional simulation. The upstream surface may operate at a phase well within the gray region, indicating an attenuating effect on the T-S wave, whereas the downstream surface responded with a phase at the edge of the gray region, indicating little effect on the T-S wave.

FIG. 11(a) shows the streamwise distribution of K near the two control surfaces in one or more embodiments. A drop in K can indicate a reduction in kinetic energy across the boundary layer and delayed growth of instability and vice versa for a rise in K. The two shaded regions of the figure portray the extent of each control surface as labeled. Comparison of K with the rigid case just downstream of each control surface is of primary interest. A marked drop in K appears just downstream of the first control surface indicating a significant stabilizing effect on the flow. Most of this stabilizing effect is maintained downstream of the second control surface as indicated by the reduction in K relative to the rigid case. This reduction in K at the control surfaces results in a shift of the skin friction coefficient in FIG. 11B, indicative of delay in boundary layer transition. Boundary layer transition delay is visually confirmed in FIG. 12 which shows a top down view of the Q-criterion iso-surface. This quantity reveals three-dimensional vortical structure in the flow. Onset of spanwise instabilities are clearly shifted by roughly 5% of the reference length.

As before, K and C_f may not change substantially over time indicating consistency in the subsurface effect. Furthermore, the K profile is well-predicted by two-dimensional simulations. This demonstrates two-dimensional laminar simulation as an accurate predictor for the relative effect of a particular subsurface on the flow. Computationally expensive three-dimensional simulations need not be applied until verification of transition delay is desired.

As noted previously, positive phasing required to achieve the strongest T-S wave attenuation implies an unrealistic negative damping coefficient for single input systems. Energy considerations shed further light on how the MIMO system bypasses this restriction. Mechanical power, P={dot over (v)}·f, was measured at both control surfaces where v is the surface displacement, the overdot indicates a derivative with respect to time, and f=∫−pdA is the forcing on the surface. Time-averaging yields P₋=−3.26×10⁻⁴ W and P+=3.26×10⁻⁴ W for the upstream and downstream surfaces, respectively. In some embodiments, the MIMO subsurface appears to extract energy from the flow at the downstream surface where T-S waves are stronger and redeposits that energy into the flow at the upstream surface.

Referring back to FIG. 10B, it can be assumed that most of the subsurface efficacy stems from the motion of the first control surface given its positive phasing and comparable amplitude. The downstream surface can be assumed to have little effect on the flow given its phasing falls at the edge of the range of phases predicted by linear theory to reduce T-S wave instability. The controller surface at x+ appears to be slightly detrimental to T-S wave mitigation. However, input from the flow at this location is necessary to elicit the favorable response at x−. Further improvement to MIMO PSub devices would benefit from a similar approach, i.e. a more neutral/smaller amplitude downstream surface to extract flow energy and a mechanical system designed to provide optimal phasing at the upstream surface.

6. CONCLUSIONS

Boundary layer control for flow transition delay was explored using high-order ILES on a flat plate with a 6:1 ellipse leading edge operating at a Reynolds number of 1,081,642 per meter length and a very low Mach number of M∞=0.04665. Disturbances were generated using an oscillating surface located at x_g=0.273 m downstream of the leading edge. T-S waves were generated at a frequency f_g=250 Hz which falls within the unstable regime. The amplitude of the generator was selected to place flow transition sufficiently far downstream of the controller interfaces.

Active control using several small-amplitude oscillating surfaces was explored first in order to determine phase and amplitude necessary for a distributed array of controllers to delay onset of flow transition. It was found that several control surfaces could significantly delay flow transition with micrometer-scale oscillations. The critical factor was the phasing of the control surface relative to the T-S wave. This was measured as the phase of displacement relative to input load, i.e., ϕ(d_i/f_i)≈ϕ(d_i/−p). Linear stability analysis predicted T-S wave attenuation for −10°<ϕ(d_i/−p)<170° while similar bounds were found in the CFD simulations.

These results were used to inform the design of a passive subsurface capable of providing amplitude and phase within the stabilizing regime. This proved a challenging task as positive ϕ(d_i/f_i) is difficult to achieve for a single input system. Therefore, a multi-input multi-output (MIMO) system was designed. Allowing the structure to interact with the flow at multiple points provided forcing at each exposed node with a known phase shift between each, which was leveraged to obtain an apparent collocated receptance (ACR) capable of exhibiting phases within the stabilizing regime.

The change in TKE for a given response (displacement/force) phasing obtained during the forced motion study, is consistent with phases of the passive ACR exhibited by each interaction surface of the MIMO PSub.

Fully-coupled passive control simulations were completed using the MIMO PSub. Positive phasing was achieved at the first control surface and slightly negative phasing at the second. Onset of flow transition, as determined by a rapid increase in C_f, was shifted downstream by roughly 5% of the reference length for this single subsurface.

Positive phasing implies the mechanical subsurface performs work on the flow which is an unrealistic expectation from a single input system. This restriction is bypassed by the MIMO design.

Energy appears to be extracted from the flow at the downstream surface where T-S waves are stronger and redeposited into the flow at the upstream surface to mitigate instability where the T-S waves are weaker. While most of the efficacy in the current system stems from the upstream surface, the downstream surface operates at the edge of the attenuation regime and in theory has little effect on the stability of the flow. Further optimization of the subsurface structure is likely to improve sensitivity and efficacy of the MIMO PSub by taking these factors into consideration.

Distribution of multiple subsurfaces may further extend boundary layer control through incremental contributions from each passive controller. Optimization of the MIMO PSub design may further delay transition. Furthermore, one or more embodiments show that computationally efficient two-dimensional simulations were sufficient to evaluate the relative effect of control surfaces. More expensive three-dimensional simulations were only necessary to confirm transition delay of final subsurface design.

Finite Element Model of Phononic Subsurface

This section gives a procedure for constructing the finite PSub FEM model from simpler components. This method consists of successive rounds of constructing block matrices of non-interacting subcomponents (i.e., the disconnected system), and then assembling these subcomponents by imposing associations between the displacements and forces of the subcomponents. In one or more embodiments, this may be achieved by expressing them in terms of a reduced set of displacements and forces (i.e., those of the assembled system). A primal assembly matrix is responsible for imposing the associations between the variables of the disconnected and assembled systems.

The PSub studied in this work is composed of flexible beams and rigid masses. The first step is to assemble a set of finite elements describing the beam dynamics to form an assembled beam system. In one or more embodiments, consistent mass and stiffness matrices for a beam element may be given as

$\begin{array}{cc}{M}_b=\frac{\mathrm{vA}\alpha }{105}\left[\begin{array}{cccc}78& 22\alpha & 27& -13\alpha \\ & 8{\alpha }^2& 13\alpha & -6\alpha \\ & & 78& -22\alpha \\ \mathrm{sym}.& & & 8{\alpha }^2\end{array}\right],& \left(A\text{.1}\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}{K}_b=\mathrm{EI}\left[\begin{array}{cccc}\frac{3}{2{\alpha }^2}& \frac{3}{2{\alpha }^2}& -\frac{3}{2{\alpha }^2}& \frac{3}{2{\alpha }^2}\\ & \frac{3}{\alpha }& -\frac{3}{2{\alpha }^2}& \frac{1}{\alpha }\\ & & \frac{3}{2{\alpha }^2}& -\frac{3}{2{\alpha }^2}\\ \mathrm{sym}.& & & \frac{2}{\alpha }\end{array}\right].& \left(A\text{.2}\right)\end{array}$

• • where Young's modulus is E, density is v, cross-sectional area is A=bh, the second moment of area is I=bh³/12 and the element half-length is α=⅛. The beam element displacement vector is d_b={{circumflex over (v)}₁, {circumflex over (v)}′₁, {circumflex over (v)}₂, {circumflex over (v)}′₂}^t, and the force vector is {circumflex over (f)}_ij={Ŝ₁, {circumflex over (M)}₁, Ŝ₂, {circumflex over (M)}₂}^t, where the subscripts are indices of the element nodes. Here, {circumflex over (V)} is displacement perpendicular to the x-direction, {circumflex over (v)}′=d{circumflex over (v)}/dx, Ŝ is the shearing force and {circumflex over (M)} is moment. Note, the superscript t refers to the transpose operation. The beam elements may be derived with cubic shape functions such that fewer elements are necessary to represent multiple wavelengths along the beam.

Each beam used in the unit cell may be formed of four beam elements (i.e., N_be=4), each having two nodes, with two nodal degrees of freedom. To constructing the total stiffness and mass matrices of the beam, the elemental mass and stiffness matrices are arranged in a block diagonal fashion as

$\begin{array}{cc}{\hat{f}}_b=\left(\left[\begin{array}{ccc}{K}_b^1& \dots & 0\\ & \ddots & \\ 0& \dots & \text{?}\end{array}\right]-{\omega }^2\left[\begin{array}{ccc}{M}_b^1& \dots & 0\\ & \ddots & \\ 0& \dots & \text{?}\end{array}\right]\right)\left\{\begin{array}{c}{d}_b^1\\ ⋮\\ \text{?}\end{array}\right\}=D_b\left(\omega \right){\hat{d}}_b,& \left(A\text{.3}\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

• • where D_i(ω)=(K_{{tilde over (b)}}−ω²M_{{tilde over (b)}}) is the disconnected dynamic stiffness matrix of the beam FEM model.

The nodal degrees of freedom of the disconnected beam elements in {circumflex over (d)}_{{tilde over (b)}} can be expressed in terms of the reduced set of degrees of freedom for a fully connected beam of the unit cell by

$\begin{array}{cc}{\hat{d}}_{\tilde{b}}=\left[\begin{array}{ccccc}\left\{\begin{array}{c}1\\ 0\end{array}\right\}& 0& 0& 0& 0\\ 0& I_{2\times 2}& 0& 0& 0\\ 0& I_{2\times 2}& 0& 0& 0\\ 0& 0& I_{2\times 2}& 0& 0\\ 0& 0& I_{2\times 2}& 0& 0\\ 0& 0& 0& I_{2\times 2}& 0\\ 0& 0& 0& I_{2\times 2}& 0\\ 0& 0& 0& 0& \left\{\begin{array}{c}1\\ 0\end{array}\right\}\end{array}\right]\left\{\begin{array}{c}{\hat{d}}_b^1\\ {\hat{d}}_b^2\\ {\hat{d}}_b^3\\ {\hat{d}}_b^4\\ {\hat{d}}_b^5\end{array}\right\}=L_b{\hat{d}}_b,& \left(A\text{.4}\right)\end{array}$

where L_b is a primal assembly matrix. Here, L_b has been constructed such that the rotational degrees of freedom of the first and last nodes of the disconnected beam element system are set to zero, thereby imposing the zero rotation condition at the end of each assembled beam system. Thus, the first and last degrees of freedom in the displacement vector of the assembled system are {circumflex over (d)}_b=v, while the motion of all other nodes are described by displacement subvectors {circumflex over (d)}_b={{circumflex over (v)}, {circumflex over (v)}′}^t. The assembled beam matrix equation is then written

${\hat{f}}_b=L_b^t\left(K_{\tilde{b}}-{\omega }^2{M}_{\tilde{b}}\right)L_b{\hat{d}}_b=\left(K_b-{\omega }^2{M}_b\right){\hat{d}}_b=D_b{\hat{d}}_b.$

In addition to beams, rigid masses are also a component of the unit cells comprising the finite PSub studied in this work. The relationship between the motion of a mass and a force applied to it is expressed as {circumflex over (f)}_m=−ω²m{circumflex over (d)}_m where m is mass. Additional mass at the nodes can be added by incorporating this relationship into the FEM. Since the rotation of the masses is restricted there is no need to include any such degree of freedom. The fully assembled beam FEM model and rigid mass relationship will serve as the simple subsystems to be connected in the next step to form the unit cell FEM model.

A dynamic stiffness matrix of the unit cell is formed by arranging the elements in a disconnected block matrix as

$\begin{array}{cc}\text{?}=\left(\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& K_b^1& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& K_b^2& 0\\ 0& 0& 0& 0& 0\end{array}\right]-{\omega }^2\left[\begin{array}{ccccc}{m}_1/2& 0& 0& 0& 0\\ 0& M_b^1& 0& 0& 0\\ 0& 0& m_2& 0& 0\\ 0& 0& 0& M_b^2& 0\\ 0& 0& 0& 0& m_1/2\end{array}\right]\right)\left\{\begin{array}{c}{\hat{d}}_m^1\\ {\hat{d}}_b^1\\ {\hat{d}}_m^2\\ {\hat{d}}_b^2\\ {\hat{d}}_m^3\end{array}\right\}=\text{?}& \left(A\text{.5}\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

where the disconnected displacement vector is, {circumflex over (d)}_ũc, is made up of displacement scalars/vectors whose superscripts are indices of each type, and subscripts refer to the types of elements comprising the model (i.e., masses or beams). Non-bold elements of Eq. A.5 are scalar quantities, while bold elements are vectors or matrices. Note the first and last mass have a value of m₁/2 such that the unit cell is defined symmetrically.

The disconnected displacement vector is linked to the assembled displacement vector of the unit cell as

$\begin{array}{cc}\text{?}& \left(A\text{.6}\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

Then, the connected dynamic stiffness matrix of the unit cell is

$\begin{array}{cc}\text{?}& \left(A\text{.7}\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

The global dynamic stiffness matrix of the finite PnC, as depicted in FIG. 7A, can be formed by first arranging three unit cells in a disconnected dynamic stiffness matrix form as

$\begin{array}{cc}\text{?}=\left[\begin{array}{ccc}{D}_{\mathrm{uc}}^1& 0& 0\\ 0& D_{\mathrm{uc}}^2& 0\\ 0& 0& D_{\mathrm{uc}}^3\end{array}\right]\left\{\begin{array}{c}{\hat{d}}_{\mathrm{uc}}^1\\ {\hat{d}}_{\mathrm{uc}}^2\\ {\hat{d}}_{\mathrm{uc}}^3\end{array}\right\}=D_f{\hat{d}}_f& \left(A\text{.8}\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

where superscripts refer to the unit cell indices. The individual unit cell dynamic stiffness matrices are assembled just as the individual beam and rigid masses were assembled to form the unit cell dynamic stiffness matrix. The displacement vector of the assembled finite PSub is related to the displacement vector of the disconnected unit cell system as

$\begin{array}{cc}{\hat{d}}_f=L\hat{d},& \left(A\text{.9}\right)\end{array}$

where {circumflex over (d)} is the displacement vector of the assembled finite PSub, and the associated dynamic stiffness matrix is D=L^tD_{{tilde over (f)}}L. The force vector of the assembled finite Psub system is then defined as {circumflex over (f)}=D{circumflex over (d)}. It should be noted that the primal assembly matrix L, is structured to impose clamped boundary shown in FIG. 7A, by setting the first and last degrees of

• • • • freedom of the disconnected unit cell system to be equal to zero.

In one or more embodiments shown or described herein, D=K−ω²M, but generally a damping matrix C can be included if necessary, as is done in Eq. (5). Additionally, the same mass, stiffness and damping matrices are used to compute time-domain solutions of the system, such as the FSI modeling procedures used in Section 5.2.

Dispersion Curve Calculation

Methods for calculating the dispersion curves of structure represented by a dynamic stiffness matrix may be widely used. Appropriate Floquet-Bloch boundary conditions impose a phase shift between the nodes located along the boundaries between unit cells, that would act as coupling surfaces in an infinite tessellation of the unit cell model. Floquet-Bloch conditions are imposed on the 1D MIMO PSub unit cell by relating the first and last nodal displacements and forces as

$\begin{array}{cc}{\hat{d}}_{\mathrm{uc}}^{N_{\mathrm{uc}}}=\text{?}& \left(B\text{.1}\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}\text{?}& \left(B\text{.2}\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

• • where N_uc is the number of degrees of freedom in the unit cell, and κ is the normalized Bloch wavenumber (i.e. product of Bloch wavenumber and unit cell length). These conditions are represented in matrix form by

$\begin{array}{cc}{\hat{d}}_{\mathrm{uc}}=\text{?}{\hat{d}}_{\mathrm{uc}}=T{\hat{d}}_{\mathrm{uc}},& \left(B\text{.3}\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}\text{?}& \left(B\text{.4}\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

where T(κ) is that matrix that encodes the Floquet-Bloch conditions, subscripts on the matrices indicate their row-column size, and the superscript * is the Hermitian transpose. Thus, {circumflex over (f)}_ǔc=T*D_ucT{circumflex over (d)}_uc=D_ǔc{circumflex over (d)}_ǔc. Finally, the dispersion curves can be obtained by solving det(D_ǔc(f, κ)=0 for eigen frequencies as a function of normalized wavenumber. The dispersion curves plotted in FIG. 7B were generated by this method.

The following examples illustrate particular properties and advantages of some of the embodiments of the present invention. Furthermore, these are examples of reduction to practice of the present invention and confirmation that the principles described in the present invention are therefore valid but should not be construed as in any way limiting the scope of the invention.

Referring to FIG. 13, a potential embodiment of a multi-input, multi-output phononic system 100 is shown. FIG. 13 shows a flow structure 102 including a cavity 103. The flow structure 102 can be, for example, a wing or other rigid horizontal structure. The flow structure 102 may be positioned in a fluid flow 124 with fluid flowing along either or both of a top surface 118 and a bottom surface 120 from a front edge 122 of the flow structure 102.

Inside the cavity 103, multiple masses may be linked together along a flexible beam structure 107 comprising one or more flexible beam segments 108 (denoted individually by an increasing number of apostrophes (′) after the number 108—unless an individual segment is specifically referred to, the use of the number “108” refers to each and all of the flexible beam segments 108-108″″″). The flexible beam segments 108 may be rigid, horizontal pieces of flexible material that support the weight of one or more masses at one or both ends of the flexible beams 108. The flexible beams 108 may be coupled to the flow structure 102 or other structure at a forward end 104 and an aft end 106 of the flexible beam structure (i.e., at a forward end of the first flexible beam segment 108 and at an aft end of the sixth flexible beam segment 108″″″). The number of flexible beam segments 108 is not limited to the number shown in the figures and any number of flexible beam segments may be used in embodiments of the phononic system 100. Generally, a mass will be positioned at a forward end of one segment and an aft end of another segment.

In some embodiments, the phononic system 100 includes a first set of masses 110, 110′, 110″ (referred to collectively as first set of masses 110 unless specifically referenced in the text), and a second set of masses 112, 112′ (referred to collectively as second set of masses unless specifically referenced in the text). The masses 110, 112 may be rigid masses. Each of the masses 110, 110′, and 110″ may be of equivalent mass and each of the masses 112, 112′ may be of equivalent mass. In embodiments, the masses 110 may be equivalent to or different than the masses 112. As alluded to above, the masses 110 and the masses 112 may be alternatively, individually coupled with one another via the flexible beam segments 108, which may enable relative motion between the masses 110, 112. The masses 110, 112 may be inhibited or prohibited from rotating about a longitudinal axis of the flexible beam structure 107 and/or individual ones of the flexible beams 108. The masses 110, 112 can be cubes or other volumetric shape (e.g., sphere, triangular pyramid, square pyramid, amorphous, etc.) In embodiments, one or more of the masses 110, 112 may be coupled to interface surfaces, such as the first fluid interface 114 and the second fluid interface 116.

The first fluid interface 114 may have a first interface surface 114a and the second fluid interface 116 may have a second interface surface 116a. The interface surfaces 114a, 116a may each be exposed to the fluid flow 124. The first fluid interface 114 and the second fluid interface 116 may be rigidly coupled to different ones of the masses. For example, the first fluid interface 114 may be rigidly coupled to the mass 110 and the second fluid interface 116 may be rigidly coupled to the mass 110′. Gradients in the velocity, pressure, and other aspects of the fluid flow 124 across the first interface surface 114a and the second interface surface 116a may cause the first fluid interface 114 and the second fluid interface 116 to move upward or downward. Motion (e.g., cantilevered motion) may be possible based on the flexibility of the flexible beam structure 107. Because the various segments of the flexible beam structure 107 are able to move separately form one another, the first fluid interface 114 and the second fluid interface 116 may move independently of one another, enabling passive flow control using the multi-input, multi-output phononic system based on the principles described herein.

FIGS. 14 and 15 shows an exemplary embodiment of a flow structure 202, which could serve as one embodiment of the flow structure 102 of FIG. 13. FIG. 14 shows the flow structure 202 including a first fluid interface 114 and a second fluid interface 216. The fluid interfaces 214, 216 may be positioned and the flow structure 202 arranged such that fluid flows across the fluid interfaces 214, 216 in a direction substantially parallel to the x-axis of the coordinates 224. The fluid interface 214 may be coupled to a mass 210 and the fluid interface 216 may be coupled to a mass 210′. FIG. 14 shows an embodiment in which a mass 212 precedes the mass 210 along the flow direction, each of the masses 210, 212 linked by a flexible segment 208 (not all of the flexible segments 208 are labeled for clarity). As shown in FIG. 15, the flexible segments 208 can extend between alternating masses in generally opposing directions as shown by the generally opposite directions of segment 208a′ and segment 208b′ and segment 208a″ and segment 208b″. In some embodiments, motion of the masses may be limited to one or more degrees of freedom (e.g., rotation of the masses may be prevented, they may be prevented from moving forward and aft as compared to the flow direction, or other degrees of freedom of movement may be prevented). As shown in FIG. 15, the motion of the mass 212 is limited to translation along an axis (i.e., in the +/−Z direction) by an axis 236. Other ones of the masses 210, 212 may be similarly constrained using similar means or other types of motion constraints.

FIGS. 16a and 16b show a phononic system 300 in a flow structure 302. The flow structure 302 includes a leading edge 322, a bottom surface 320, and a top surface 318 with a fixed portion 318a and a removable portion 318b. The removable portion 318b may fit within a cavity 303 of the flow structure 302 and may interface with the flow structure 302 along a rim 350 of the cavity 303. The removable portion 318b may include a first slit 315 and a second slit 317. The first slit 315 and the second slit 317 may provide a space for the first fluid interface 314 and the second fluid interface 316 to interact with the fluid flow 324 in contact with the flow structure 302. FIG. 16a shows the removable portion of the top surface 318 removed, exposing the phononic system 300. As fluid flows across the flow structure 302, it may interact with the first fluid interface 314 and the second fluid interface 316 based on the principles described herein enabling passive flow control using the multi-input, multi-output phononic system.

While the present invention has been illustrated by a description of one or more embodiments thereof and while these embodiments have been described in considerable detail, they are not intended to restrict or in any way limit the scope of the appended claims to such detail. Additional advantages and modifications will readily appear to those skilled in the art. The invention in its broader aspects is therefore not limited to the specific details, representative apparatus and method, and illustrative examples shown and described. Accordingly, departures may be made from such details without departing from the scope of the general inventive concept.

Claims

1. A multi-input, multi-output phononic system comprising: a first interface surface and a second interface surface that respond to at least one of a pressure gradient or a velocity gradient in a wave of a turbulent fluid flow or a laminar fluid flow, the pressure gradient or the velocity gradient associated with complex motion of the flow exhibiting a plurality of frequencies exerted on one or more of the interface surfaces; anda subsurface feature extending from the interface surfaces, the subsurface feature comprising a phononic crystal or locally resonant metamaterial adapted to receive one or more of the pressure gradient or the velocity gradient from the fluid flow via the interface surfaces and to alter one or more of a phase and an amplitude of a plurality of frequency components of the fluid flow.

2. The multi-input, multi-output phononic system of claim 1, wherein: the subsurface feature is a homogeneous, uniform elastic structure that couples the first interface surface and the second interface surface, the multi-interface phononic system adapted to receive the wave via the interface surfaces and to alter a phase of the wave, andthe set of interface surfaces is adapted to vibrate at a frequency, phase, and amplitude in response to the altered phase of the wave.

3. The multi-input, multi-output phononic system of claim 2, further comprising a plurality of masses linked by a flexible beam structure, wherein the rotation of the masses about a longitudinal axis of the flexible beam structure is inhibited.

4. The multi-input, multi-output phononic system of claim 3, wherein the plurality of rigid masses comprises a first set of masses, each of the masses in the first set of masses having an equivalent mass,a second set of masses, each of the masses in the second set of masses having an equivalent mass, whereinindividual ones of the first set of masses and individual ones of the second set of masses alternate along the length of the flexible beams.

5. The multi-input, multi-output phononic system of claim 4, wherein the first interface surface is rigidly coupled to a first mass in the first set of masses and the second interface surface is rigidly coupled to a second mass in the first set of masses, and a mass of the second set of masses is coupled between the first mass in the first set of masses and the second mass in the first set of masses.

6. The multi-input, multi-output phononic system of claim 5, wherein interaction between the first interface surface, the second interface surface, and the fluid flow enable a positive phasing.

7. The multi-input, multi-output phononic system of claim 6, wherein the subsurface feature is adapted to reduce kinetic energy within the flow.

8. The multi-input, multi-output phononic system of claim 7, wherein the phononic subsurface feature is adapted to reduce formation or development of an energy cascade characteristic of partially or fully developed turbulence.

9. The multi-input, multi-output phononic system of claim 1, wherein the fluid flows at least partially as a Tollmien-Schlichting (T-S) wave.

10. The multi-input, multi-output phononic system of claim 5, wherein interaction between the first interface surface, the second interface surface, and the fluid flow enable a negative phasing, between the displacement and the forcing induced by the flow, at one or more interaction surfaces.

11. The multi-input, multi-output phononic system of claim 10 wherein the subsurface feature is adapted to increase kinetic energy within the flow.

12. A method of controlling a flow comprising: providing a first interface surface and a second interface surface in a fluid flow along a flow structure, the fluid flow exhibiting turbulent or laminar flow characteristics, the first interface surface and the second interface surface arranged within the flow structure such that they each are exposed to one or more of a pressure gradient and a velocity gradient associated with complex motion of the fluid flow; andproviding a subsurface feature extending from the interface surfaces, the subsurface feature comprising a phononic crystal or locally resonant metamaterial adapted to receive one or more of the pressure gradient or the velocity gradient from the fluid flow via the interface surfaces and to alter one or more of a phase and an amplitude of a plurality of frequency components of the fluid flow; andpassively altering one or more of a phase or amplitude of a plurality of frequency components of the flow via the subsurface feature.

13. The method of claim 12, further comprising vibrating the set of interface surfaces at a phase and amplitude of a plurality of frequency components of the flow.

14. The method of claim 13 wherein the set of interface surfaces are physically coupled to one or more masses in the subsurface feature, and the masses are tuned such that the interface surfaces vibrate at a plurality of frequencies in response to the pressure and/or velocity gradients of the fluid flow.

15. The method of claim 14, wherein the first interface surface and the second interface surfaces are physically coupled to equivalent masses.

16. The method of claim 15, wherein the equivalent masses physically coupled to the first interface surface and the second interface surface separated by a mass of a different amount than the equivalent masses.

17. The method of claim 12, wherein the first interface surface is upstream of the second interface surface and there is a phase shift between the forcing of the first interface surface and the second interface surface.

18. A system for reducing skin friction of a surface using multi-input, multi-output phononic subsurfaces for passive boundary layer transition delay comprising: a first interface surface and a second interface surface that respond to at least one of a pressure gradient or a velocity gradient in a wave of a turbulent fluid flow or a laminar fluid flow, the pressure gradient or the velocity gradient associated with complex motion of the fluid flow exhibiting a plurality of frequencies exerted on one or more of the interface surfaces; anda subsurface feature comprising a plurality of masses extending from the interface surfaces, the subsurface feature adapted to receive one or more of the pressure gradient or the velocity gradient from the fluid flow via the interface surfaces and to passively alter one or more of a phase and an amplitude of a plurality of frequency components of the fluid flow.

19. The system of claim 18, wherein the plurality of masses comprises: a first set of masses comprising a first mass, anda second set of masses comprising a second mass, andthe first set and the second set are arranged in an alternating pattern and connected with one another in series via a plurality of flexible segment beams.

20. The system of claim 19, wherein the flexible segment beams are linked to form a flexible beam structure and the flexible beam structure is anchored to a flow structure at a forward end and an aft end and can move perpendicularly in a plane parallel with respect to a direction of fluid flow across the flow structure.