Attenuation of Tollmien-Schlichting waves using resonating surface-embedded phononic crystals

A novel method for control of convective boundary layer instabilities using metamaterial concepts is investigated. Attenuation of Tollmien–Schlichting (TS) waves with surface-embedded one-dimensional phononic crystals (PCs) is theoretically and numerically modeled, capitalizing on the inherent frequency band stop of PCs. The PC is tuned to the targeted TS wave characteristics through the use of analytical models derived from transfer matrix and interface response theories, veriﬁed using a ﬁnite elements analysis. The interaction between TS waves and a single PC is investigated using coupled two-dimensional ﬂuid structure interaction simulations in the frequency domain. It is shown that TS waves are either ampliﬁed or attenuated depending on whether the PC free-face surface displacement and unsteady perturbation pressure at the wall are in-phase or out-of-phase, respectively. The perturbation pressure acts solely as the driver for the mechanical oscillation of the PC. The emerging hydrodynamic coupling between TS waves and the PC is found to be governed by a combination of the Orr mechanism and wall-normal velocity linear superposition near the wall. Finally, a metasurface comprised of an array of streamwise-distributed PCs is evaluated, resulting in an amplitude growth delay of 11.3% of the TS wavelength along the metasurface extent.


I. INTRODUCTION
The process of laminar to turbulent transition of two dimensional subsonic boundary layers under low external disturbance conditions is typically governed by the growth and eventual breakdown of Tollmien-Schlichting (TS) waves. 1 Suppression or mitigation of TS wave growth and, in consequence, transition delay has been the point of focus of numerous studies in the past decades due to direct implications on skin friction reduction and, hence, improvement of efficiency in a broad range of fluid flow applications.
The wavelike, two-dimensional, and initially linear behavior of TS modes motivated early active control attempts based on the principle of wave superposition.In these cases, a perturbation of equal amplitude but opposite phase is actively introduced in the boundary layer, 2,3 often in the vicinity or at the wall surface by means of devices, such as vibrating ribbons and plasma actuators.This approach has been successful in controlling artificially introduced monochromatic TS waves and set grounds for later development of advanced adaptive closed loop control methods, 4,5 which directly target non-deterministic, naturally occurring TS waves.These methods are typically either model-free and filter-based (e.g., Refs.6 and 7) or model-based (e.g., Refs.8-10) approaches.
A common trait in the aforementioned works is the use of active control schemes; thus, there is a requirement of sensors and externally powered actuators in order to assess the characteristics of the incoming TS waves and, through the appropriate control logic, to generate the required canceling disturbance.
At the same time, apart from the aforementioned active methodologies, numerous passive TS wave control approaches have been developed as they are generally characterized by simplicity and no requirement of energy expenditure.Passive control methods typically take advantage of the coupled dynamic fluidstructure interactions occurring at the interface between the boundary layer and the wall.An eminent example is compliant surfaces, 11,12 which function by manifestation of elastic waves along the surface caused by the TS waves themselves as their low surface stiffness gives rise to surface motions that interact with the flow.Yet, effective compliance is difficult to achieve, particularly when TS waves are considered in air, since the impedance disparity between the fluid and the flexible wall surface is very large.In addition, a critical limitation of compliant surfaces is the emergence of Rayleigh-type flutter modes that travel along the fluid-structure interface and which are known to accelerate transition.
The wavelike behavior intrinsic to boundary layer instabilities, brings forward new opportunities for control, by using local wallresonance phenomena.These can be deeply exploited through the concept of metamaterials, particularly those aimed at acoustic applications (e.g., Refs.13 and 14) Metamaterials have been proven to be excellent passive wave manipulators in various disciplines and, hence, are of potential interest when considering the wavelike behavior of convective TS wave instabilities.Furthermore, they constitute a relatively simple, passive control strategy that maintains considerable selfadaptability when it comes to varying forcing conditions.
To the authors' best knowledge, literature regarding the stabilization of convective flow phenomena using metamaterials or metasurfaces is very limited.Acoustic metasurfaces in the form of shallow cavities have been proposed for stabilizing hypersonic boundary layers. 15The surface achieves a near-zero impedance condition at the wall, hence, suppressing the growth of the second Mack mode.7][18] Therein, the frequency dispersion of the phonon is adjusted to form either a stop band or a passband for stabilization and amplification, respectively.
It must be emphasized that the particularities of convective boundary instabilities bring forward an important theoretical consideration on the operation of a metasurface.Specifically, in "classical" metamaterial systems such as acoustic manipulators, the propagation direction of the targeted wave is typically parallel to the axis of periodicity of the meta-atom (or phonon).The acoustic wave is, thus, propagating through the metamaterial.In contrast, TS waves propagate parallel to the solid wall (corresponding to the skin of an aerodynamic body, such as a wing) and in consequence normal to the axis of the wall-embedded metamaterial.
The salient differences in the mode of interaction between TS waves and metasurfaces are an important aspect in need of clarification.With this in mind, the aim of this work is to further assess the TS damping capabilities of one-dimensional multi-layered phononic crystals (PCs) arranged normal to the fluid/solid interface and at the same time identify the mechanisms of the interactions with the flow.Furthermore, the possibility of amplitude reduction by means of a metasurface that comprises of a series of streamwise distributed identical PCs is explored.Hereby, it is important to note that in the problem at hand, band stopping is only experienced by the elastic waves that propagate through the PC.Instead, TS waves are subjected to a wall perturbation at the interface whose amplitude and phase are governed by the mechanical resonance characteristics of the sub-surface; hence, the PC is perceived by the boundary layer as a local resonator.
Within the framework of this investigation, first, a methodology for the design of a single PC is established that is based on both theoretical modeling and finite element analysis.The interactions between the TS waves and a single PC as well as the PC metasurface are subsequently investigated by means of fluid-structure interaction (FSI) simulations, coupling the problems of Linearized Navier-Stokes (LNS) equations and elastic wave propagation in the frequency domain.

II. FLOW CONDITIONS AND NUMERICAL SETUP
A schematic representation of the problem, an overview of the computational domain, and the pertinent boundary conditions are provided in Fig. 1.TS waves develop on a 1 m long, zero pressure gradient, flat plate.The waves are artificially forced by a localized volume body force centered at a distance of 0:18 m from the leading edge ðx ¼ y ¼ 0Þ.As shown from linear stability analysis later in this section, at this location, TS waves are damped before amplifying, ensuring a stable perturbation.A one-dimensional, multi-layered, and multicelled PC is placed normal to the plate at x ¼ 0:65 m.Throughout this report, unless units are mentioned, quantities are scaled with characteristic quantities at the PC location, namely, the steady-state baseline boundary layer displacement thickness ðd Ã ¼ 1:2 mmÞ, freestream velocity ðU 1 ¼ 20:4 m=sÞ, and freestream static pressure ðp 1 ¼ 4:5 PaÞ.Hence, frequencies are scaled with d Ã =U 1 .
The computational domain comprises of two parts: a fluid and a solid domain.The fluid (air) temperature, pressure, and dynamic viscosity are set to 20 C; 1 atm, and 1:825 Â 10 À5 Pa s, respectively.The solid domain represents the PC and is described by the constitutive properties of elastic solids.The problem is discretized with the secondorder finite element basis functions in COMSOL multiphysics 5.6.The authors have validated this solver against experiments in similar conditions. 19Both fluid and solid domains are meshed with rectangular elements that are refined at the vicinity of the wall, the smallest elements being 6:3 Â 10 À5 m, approximately 20 times smaller than the displacement thickness at the PC location.The suitability of this mesh has been verified by a preliminary mesh convergence study, and sample results of which are shown in Table I.
Fluid structure interaction simulations are performed in two steps.First a steady-state simulation is run with which the baseline, time-invariant conditions are established.Hereby, the flow field is described by the steady, incompressible Navier-Stokes (NS) equations.

FIG.
1.An illustration of the geometry and simulation boundary conditions.Black and green font descriptors pertain to the boundary conditions for the steady state and frequency domain simulations, respectively.PML: perfectly matched layer.
In turn, the solid displacement of the PC is described by the steady balance of momentum for geometrically non-linear, isotropic, elastic solids, relating the second Piola-Kirchoff stress 20 to the Green-Lagrange strain. 21The two-way direct coupling of the fluid-structure interface is performed with an arbitrary Lagrangean-Eulerean method, 22 which ensures stress continuity between the fluid and the solid, i.e., pressure and viscous forces of the fluid are equated to the reaction force normal to the solid boundary, where n, p, I, l, and u f are the normal vector to the boundary, pressure, the identity matrix, fluid dynamic viscosity, and the fluid velocity vector, respectively.The second simulation step is performed in the frequency domain, solving the harmonic linearized incompressible Navier-Stokes equations, using the steady base flow solution extracted from the steady state simulation step as the mean (i.e., time-invariant) component.Perfectly matched layer (PML) inlet and outlet domains absorb convective perturbations (i.e., TS waves) by coordinate stretching. 23Regarding the solid domain, the same equations that describe the steady-state FSI are linearized and transformed into the frequency domain.A coupling that ensures continuity of velocities and stresses across the interface is employed, formulated as where u f and g s are vectors of fluid velocity and solid displacement, respectively.For both steady-state and frequency domain simulations, pressure gradient effects due to the leading edge are accounted for by means of a symmetry condition 24 extending upstream of the leading edge.Spatial linear stability theory (LST) formulated on the incompressible Orr-Sommerfeld equation 25 is employed on the baseline steady state solutions for identifying the unstable TS wave frequencies.Growth rates ða i Þ and amplification factors, Àa i dx: (3) x 0 being the neutral point, are shown in Fig. 2(a).A frequency of 0.0171 ð300 HzÞ is chosen for tuning the PC resonance mode.A TS wave of this frequency will have undergone substantial linear growth, achieving a sufficiently high amplification factor (N ¼ 2.5) when reaching the location of the PC.
In turn, as seen in Fig. 2(b), LST wavelength predictions for a TS wave with frequency of 0.0171 ð300 HzÞ range between 17 ð20 mmÞ and 21 ð25 mmÞ for the whole length of the flat plate (19.9, 23:2 mm at the PC).Given the aforementioned wavelength estimates of the expected TS waves and taking into account future experimental practicality, the width of the PC interface with the flow is set to 0.86 ð1 mmÞ, approximately one order of magnitude lower than the respective TS wavelength.

III. DESIGN AND FREQUENCY RESPONSE OF THE PHONONIC CRYSTAL
Notwithstanding the numerical methodology followed in this work, the design of the PC is performed provisioning for future experimental feasibility; hence, one of the primary goals is to obtain practical dimensions.Hereby, a one-dimensional design is considered; thus, the periodicity of the elastic properties and the displacement of the PC are prevalent along the y-direction (Fig. 1), under the assumption that the TS wave pressure load applied on the PC interface is much greater than local frictional loads.
Two aspects are considered for the design of the PC's unit cell.First, a frequency bandgap must be defined, within which elastic waves propagating through the one-dimensional PC are attenuated.This is based on the hypothesis that an out-of-phase surface response within the bandgap leads to destructive interference of TS waves at the fluidsolid interface. 16Second, the resonance frequency of the PC must lie TABLE I. Mesh convergence study.Amplitude of u 0 and v 0 at the TS wave peak for x ¼ 557:5 ð0:65 mÞ and f ¼ 0:0171 ð300 HzÞ.

Mesh
Elements ju 0 j ðm=sÞ j v 0 j ðm=sÞ within the bandgap as the displacement and, in consequence, the velocity of the PC interface is the largest; hence, the amplitude of interactions with the flow above the interface is maximized.
Given the above, a three-layer unit cell design is chosen since it allows to achieve low resonance frequencies while keeping the PC compact.The material properties of each layer and their dimensions are given in Table II.The bandgap of the three-layer unit cell can be analytically described via the transfer matrix analysis. 26,27In this case, the analysis yields where q is Bloch's wavenumber (i.e., the wavenumber of the elastic compression wave through the PC), D is the length of the unit cell, , and Z i ¼ q i c i , where x, d i , c i , and q i are the angular frequency, layer thickness, longitudinal speed of sound, and mass density of the ith layer (i ¼ 1, 2, 3), respectively.The longitudinal speed of sound can be expressed in terms of the Young's modulus (E) and Poisson's ratio () of the layer's material through the following expression: 28 In finding the bandgap with Eq. ( 3), the angular frequency and layer properties are given as inputs and q is solved for.The bandgap exists, where q is a complex quantity, hence, where the right-hand side of Eq. ( 3) is greater than unity.In turn, the resonance frequency within the bandgap is analytically derived following interface response theory (IRT), 13 which, for a three-layer unit cell, reads The material properties of each individual layer are input to Eq. ( 5), which is solved for the resonance frequency, x.The chosen layer thicknesses are 30 (35 mm) of rubber, 2.6 (3 mm) of aluminum, and 4.3 (5 mm) of steel, resulting in a total unit cell length of 36.9 (43 mm).The corresponding dispersion curve of the unit cell and the resonance frequency are shown in Fig. 3.The design of the PC is finalized by choosing the number of unit cells (N UC ).To this end, the finite element analysis is used to identify the effect of varying N UC on the PC eigenfrequencies.Two variations of boundary conditions are   examined.In the first case, the right and left sides of the PC are constrained in the x direction, while the top and bottom faces are left unconstrained.This set of boundary conditions (BC1) is used to verify Eq. ( 5).The second set of boundary conditions (BC2) is the same as the first, except that the bottom is fixed in both x and y directions.This set is a closer representation of practical experimental conditions and is, therefore, selected for the FSI analyses discussed.
Based on the IRT calculations, the first resonance frequency, which corresponds to the surface mode of the PC, is the N th UC eigenfrequency.The eigenmode corresponding to the N th UC eigenfrequency is verified to be a surface mode by inspection of displacement, as the deformation of the eigenmode is maximized on the surface of the PC.The surface mode resonance frequency for different values of N UC for the two sets of boundary conditions is given in Table III.Evidently, the surface mode resonance frequency predicted by IRT agrees very well with finite elements analysis (FEA) when BC1 is employed.As expected, due to the different problem formulation, the prediction of surface mode resonance frequency deviates from the IRT results when BC2 is used, as the latter inherently introduces non-uniformity in the response of the PC along the yÀdirection.This deviation is substantially reduced as N UC increases since the surface mode is gradually attenuated across the successive layers, i.e., converging to the IRT formulation for increasing number of layers.Based on the above results, a value of N UC ¼ 5 is used to finalize the PC design, noting that a three-layer design would also be acceptable (within 0.5% of the predicted response).
The complex frequency response of the PC to TS wave pressure disturbances of various wave frequencies at the fluid-structure interface pertaining to the LNS solution with BC2 is shown in Fig. 4. The amplitude response of the streamwise surface velocity, ju 0 =p 0 j, is found to be of the order of 10 À8 , far lower compared to TS wave velocities, which are of the order of 10 À3 to 10 À2 .Hence, the influence of the streamwise displacement of the PC can be considered negligible.Instead, the wall-normal surface velocity response, jv 0 =p 0 j, is of the order of 10 À3 , sufficient to impact the flow.
As shown in Fig. 4(b), the phase difference between the PC surface pressure and surface wall-normal displacement, /ðg 0 =p 0 Þ, is zero and Àp for frequencies below and above resonance, respectively.As expected, the corresponding phase difference between the surface pressure and surface wall normal velocity, /ðv 0 =p 0 Þ, follows the same trend, albeit with a p=2 lag, i.e., the phase difference between displacement and velocity.Given the above, for the remainder of this report, the response of the PC is classified into three regimes: "in-phase," "resonance," and "out-of-phase" when referring to below, at, or above the resonance frequency, respectively.

IV. INTERACTION OF A SINGLE PC AND A PC METASURFACE WITH TS WAVES
Following general approaches emerging from compliant surface studies (e.g., Ref. developing TS waves.These are the integrals along y of kinetic energy ðe k Þ, streamwise, and wall-normal velocity fluctuations.The variation of these variables along the x-direction at in-phase (0.0170, 298 Hz), resonance (0.0172, 300.71Hz), and out-of-phase conditions [0.0174, 304 Hz, Fig. 4(b)] is shown in Fig. 5.In all cases, it is observed that the influence of the PC is local rather than global and, at the same time, it is relatively small except at resonance.This can be justified by considering the normal mechanical admittance at the surface. 11,30dmittance follows the same trend as the displacement frequency response seen in Fig. 4(b).Consequently, far from resonance and on either phase direction, admittance is greatly reduced resulting in a lower amplitude of displacement.
Considering the integral of kinetic energy [Figs.5(a)-5(c)], it can be seen that the influence of the PC extends to both upstream and downstream directions, with an increase and decrease in e k at in-phase and out-of-phase conditions, respectively.Hence, locally, the PC either decreases or increases the TS wave amplitude.In addition, it can be concluded that a single PC is not capable of achieving persistent downstream TS wave amplitude reduction.In turn, at resonance, there is an abrupt increase in e k just upstream of the PC, followed by a sharp decrease over the PC.This causes substantial amplitude increase that persists downstream, following a behavior closely reminiscent of roughness induced scattering. 19,31,32he aforementioned observations in e k are directly related to the behavior of ju 0 j and jv 0 j as shown in Figs.5(d)-5(f) and 5(g)-5(i), respectively.As mentioned earlier, the premise of control relies on the excitation of the PC at the interface by the pressure fluctuations incurred by the TS waves.In turn, the PC produces a surface displacement and velocity, which, if tuned properly, counters the TS wave velocity perturbations.At in-phase conditions, jv 0 j is reduced with a simultaneous increase in ju 0 j, resulting in an overall increase in e k .The opposite is observed for the out-of-phase conditions, where an increase in jv 0 j and a decrease in ju 0 j occur, leading to an overall decrease in e k .An important observation is that increase or decrease in TS wave FIG. 6.Amplification factor at x ¼ 685 as a result of replacing the PC interface with an oscillating boundary condition of variable phase at out-of-phase conditions (0.0174, 304 Hz).Dashed line: Baseline, N base ; bottom axis: phase difference between the oscillating surface vertical velocity component and the TS wave pressure; top axis: phase difference between the TS wave and the oscillating surface vertical velocity component.The red mark corresponds to the favorable phase difference achieved with the current PC design.

FIG. 7.
Comparison of shape functions (a)-(c) and phase angles (d)-(e) of u 0 ; v 0 , and p 0 between baseline, oscillating wall, and PC cases pertaining to out-of-phase conditions (0.0174, 304 Hz) at the PC location.The insets in (a)-(c) correspond to the difference with respect to the baseline case (e.g., D ¼ ju 0 j À ju 0 j base Þ.
amplitude at in-phase or out-of-phase conditions does not take place at the PC location where e k remains unchanged.Instead, the effects pertain to the immediate upstream and downstream regions.
In order to elucidate the mechanism behind the observed amplitude reduction, the optimal scenario for TS wave attenuation must first be identified.To do so, the PC is replaced by an active oscillating wall of the same width, modeled in the FSI simulation by an imposed boundary condition on v, excited at the out-of-phase frequency (0.0174, 304 Hz) shown in Fig. 4(b).The velocity amplitude at the wall is matched to the local baseline TS wave amplitude, and the phase is varied with respect to the TS wave over a full cycle.Figure 6 depicts the results in terms of amplification factor downstream of the location of forcing, relative to the baseline.Evidently, the maximum reduction of N is achieved when the phase difference between jv 0 j and jp 0 j at the PC surface is zero.An asymmetry of the amplification factor with respect to the baseline is also noted in the figure which is biased in favor of N reduction, i.e., in the interval of À0:65p to 0:6p.The current PC design is, therefore, within the favorable phase interval for out-of-phase conditions [p=2, equivalent to À3p=2 shown in Fig. 4(b)], however, that does not exhibit the optimal phase relation.
Further insight into the salient interactions leading to the observed attenuation of the TS wave can be obtained by considering the shape functions of u 0 ; v 0 , and p 0 at the location of the PC. Figure 7 compares the aforementioned shape functions with the baseline as well as the oscillating wall cases.Notably, the same trends are evident for all shape functions between the PC and oscillating wall cases.Hence, it can be concluded that the PC mode of interaction is equivalent to an oscillating wall acting in the wall-normal direction with interactions being restricted within the local displacement thickness.Consequently, an increase in jv 0 j is observed at the surface with a simultaneous reduction of ju 0 j [Figs.7(a) and 7(b)], concurring the trends of the out-of-phase integrals shown in Fig. 5. Given the twodimensionality of the problem at hand, this behavior can be regarded as a Reynolds stress perturbation, thus pointing to the existence of an Orr mechanism (e.g., Refs.33 and 34).At the same time, as seen in Fig. 7(e), the phase of v 0 near the wall is reduced with respect to the baseline, opposing the maximum of jv 0 j of the TS wave (see Fig. 6).Given the above considerations and the linearity that governs the problem, it is concluded that the mechanism of interaction of the PC with the TS wave is a combination of an Orr process and linear superposition of the vertical velocity components near the surface.In addition, since no changes in the phase of p 0 [Fig.6(f)] are identified, the latter is found to be purely responsible for driving the PC and does not play a role in the attenuation mechanism.
In order to investigate the possibility of sustained attenuation along the streamwise direction, a metasurface comprising of twenty identical PCs is constructed.The pacing between PCs is 1.72 ð2 mmÞ; hence, sub-wavelength is relative to the TS wave.Unlike compliant surfaces, the benefit of a metasurface is that it does not exhibit streamwise traveling surface modes, i.e., PCs only interact with each other through pressure fluctuations.Stability is hereby assessed based on the amplification factors of ju 0 j; jv 0 j, and e k as shown in Fig. 8.In line with the observations for the single PC, the trend of decrease in ju 0 j and increase in jv 0 j occurs throughout the whole metasurface resulting in an overall reduction of e k .Through the combined action of the metasurface, the effective amplitude growth of the disturbance is delayed by approximately 11.3% of the TS wavelength.While this decrease is relatively mild, the overall potential for improvement can be increased by adjustment of the spacing between the PCs, their width, and the streamwise extent of the metasurface.At the end of the metasurface, amplification is observed due to the abrupt change in the boundary condition to a solid wall.This, however, does not pose concerns regarding the practicality of using PCs for transition delay as the metasurface can either be extended up to the end of the geometry at hand, or be designed to extend from the linear up to the non-linear and turbulent regimes, where the depicted linear amplification factor definition no longer holds.

V. CONCLUDING REMARKS
In conclusion, a novel technique for attenuation of TS waves and transition delay is investigated, employing PCs which serve as the units FIG. 8. Amplification factors based on the maxima of (a) ju 0 j, (b) jv 0 j, and (c) e k , in out-of-phase conditions (0.0174, 304 Hz) with twenty identical PCs located at the shaded regions.
for a passive metasurface.A methodology of the design of the individual PC characteristics based on analytical models is provided, agreeing with FSI simulations.The results for a single PC show a low magnitude local amplitude reduction at the vicinity of the PC when the surface displacement and the TS wave pressure perturbation are outof-phase and amplification when they are in phase.Pressure fluctuations are shown to be solely responsible as the driver of the PC displacement without directly participating in the interaction mechanism.Instead, the simultaneous decrease in streamwise and increase in wall-normal fluctuations, in conjunction with the opposing phases between the PC surface and the TS wave wall-normal velocity, suggest that attenuation is the result of a near-wall Orr mechanism and velocity superposition.In turn, with the introduction of a metasurface comprised of an array of identical PCs of sub-wavelength spacing, attenuation is extended downstream and transition delay of 11.3% of the TS wavelength is achieved.

FIG. 4 .
FIG. 4. (a) Streamwise velocity response of the PC surface.(b) Wall-normal velocity response of the PC surface along with the phase differences of wall-normal displacement (solid line) and wall-normal velocity (dashed line) with respect to surface pressure.This result pertains to the LNS solution with BC2.From left to right, arrows indicate in-phase (0.0170, 298 Hz), resonance (0.0172, 300.71Hz), and outof-phase (0.0174, 304 Hz) conditions shown in this report, respectively.

TABLE II .
Properties of the PC layers.

TABLE III .
Estimation of f SM from finite elements and IRT.