Collective dynamics of dense hairy surfaces in turbulent flow

Flexible filamentous beds interacting with a turbulent flow represent a fundamental setting for many environmental phenomena, e.g., aquatic canopies in marine current. Exploiting direct numerical simulations at high Reynolds number where the canopy stems are modelled individually, we provide evidence on the essential features of the honami/monami collective motion experienced by hairy surfaces over a range of different flexibilities, i.e., Cauchy number. Our findings clearly confirm that the collective motion is essentially driven by fluid flow turbulence, with the canopy having in this respect a fully-passive behavior. Instead, some features pertaining to the structural response turn out to manifest in the motion of the individual canopy elements when focusing, in particular, on the spanwise oscillation and/or on sufficiently small Cauchy numbers.


Introduction
The coherent waving motion of seagrass meadows in marine currents or plant crops under the action of wind is a fascinating example of the complex interaction between fluid flows and a multitude of deformable structures, namely hairy surfaces or filamentous beds.More generally, the emergence of large-scale collective motion (e.g., the so-called honami or monami phenomena in terrestrial and aquatic canopy flows 1 , respectively) represents an important reason why this kind of flow-structure interaction (FSI) keeps attracting interest from both the fundamental and applied viewpoint [2][3][4][5][6] .Such intriguing features indeed manifest in a vast number of problems pertaining to various fields, ranging from aforementioned environmental processes (e.g., terrestrial or aquatic canopy flows) 2,5,[7][8][9] to engineering applications (e.g., energy harvesting and drag reduction) 10,11 and, in a broader sense, biological systems (e.g., ciliated organisms) [12][13][14] .
Environmental phenomena arguably represent the most classical example where collective motion can be observed and numerous studies have been devoted to understand the key features of the mutual coupling between the fluid and structural motion 8,15 .Along with the decreased drag due to the adaption of compliant stems to the action of the flow, the canopy flexibility is typically responsible for an attenuation of vertical mixing between the outer flow and the canopy 2,16 .Yet, several questions of paramount importance remain not fully understood.Firstly, the mathematical description of canopy flows is still the subject of theoretical studies aiming at providing comprehensive models able to predict the onset of the honami/monami-like instability 3,15,[17][18][19][20][21][22] .On the other hand, despite their importance for understanding the incipient cause of such phenomena, these models typically rely on linear approximations, and it remains elusive to predict and measure the characteristic properties of fully-developed honami/monami in the nonlinear stage.Furthermore, from the current literature it can be noticed conflicting evidence concerning the identification of the main features of their collective motion, such as its characteristic frequency and lengthscale 4,5,15 .It is also important to highlight the technical challenges for experimental investigations, in particular, for measuring the spatial conformation of the fluid flow in the intra-canopy region 5 , as well as isolating the effect of the various governing parameters (e.g., Cauchy and Reynolds number, solid-to-fluid density ratio, reduced velocity).Both issues, on the other hand, can be overcome by the support of computational studies based on direct or large-eddy numerical simulations [23][24][25][26] , with more recent advancements based on fully-resolved approaches where the stems/elements of the canopy are resolved individually and not modelled, e.g., as a continuous porous medium 4,6,11,[27][28][29] .
In this work, we aim at deciphering the highlighted issues by means of a high-fidelity computational approach to accurately and comprehensively access both the flow physics and canopy collective dynamics, as well as the motion of the individual stems.Overall, we aim at establishing whether the fluid and structural dynamics act competitively or not, and which of the two has a dominant role.Moreover, we explore the connection between the individual motion of the elements (e.g., filaments or blades) and the collective coherent motion.

Results
To investigate the essential features of the FSI between a turbulent flow and a flexible canopy, we have performed high-fidelity, direct numerical simulations (DNS) where the canopy elements are modelled individually and where the two-way coupling between the fluid and the solid dynamics is realized by means of an immersed boundary method (IBM).The considered setup is that of an open channel at constant flow rate with the flexible-canopy elements randomly distributed 28 and individually clamped at the bottom boundary, as shown in Fig. 1 for a representative case.See Methods section for a complete description of the problem setup and numerical method.
We fix the hydrodynamic parameters to a reference configuration with bulk Reynolds number Re = UH/ν = 5000 (where U is the imposed bulk flow velocity, H the height of the open channel and ν the kinematic viscosity), with a submergence level h/H equal to 0.25 (where h is the height of the canopy in the rigid vertical configuration) and with canopy solidity λ = h d/∆S2 ≈ 1.43 (with d being the stem diameter and ∆S the average spacing between adjacent stems).The latter choice corresponds (in the rigid case) to a configuration well within the dense regime, thus avoiding the transitional regime whose definition remains debated 8,28,30 .The selected solid-to-fluid density ratio is ρ s /ρ f = O(1) in order to consider the case of almost neutrally-buoyant filaments.Since we aim at isolating the essential flow-structure interaction mechanisms, we neglect the presence of gravity (i.e., the buoyancy and Galileo number are set to zero).
We have performed a set of simulations varying the Cauchy number Ca, such parameter representing the ratio between the forcing exerted by the flow and the elastic restoring force exerted by the filaments, defined as where γ is the bending stiffness of the filaments.Specifically, the following values of the Cauchy number have been chosen: Ca = {0, 1, 10, 25, 50, 100}, with Ca = 0 corresponding to the rigid canopy case where the stems cannot deform.Such range is expected to include the variety of possible dynamical regimes for flexible canopies, both in terms of the individual element motion (i.e., weak vs strong reconfiguration) 25 as well as collective motion (i.e., uncorrelated gentle swaying vs honami/monami) 31 .Note that the analysis can also be expressed in terms of a variation of the reduced velocity U r = U/( f n,1 h) ∼ √ Ca, where f n,1 is the natural frequency of the filaments.For the chosen configurations, the (average) canopy height expressed in wall units is overall O(10 2 ), monotonically decreasing with the Cauchy number from ∼ 280 in the rigid-canopy case (i.e., Ca = 0) to ∼ 95 in the most flexible case (Ca = 100) (note that the friction velocity used here is that evaluated at the virtual wall origin, i.e. the virtual wall seen by the outer boundary-layer 27 ).To give a qualitative indication, snapshots from the representative case at Ca = 50 are shown in Fig. 1, from which it can be observed the large-scale coherent motion undergone by the filaments (left panel), and superimposed the large-scale flow structures in the proximity of the canopy top (right panel).As it will be shown, the two are strictly related and represent the key features of the underlying physical mechanisms.A revised view of honami/monami Fig. 2 (top panels) shows some instantaneous configurations of the canopy for several values of Ca considered in our study.From the elevation of the filaments it can be observed the typical pattern of the honami/monami motion, which appears more pronounced as the Cauchy number is increased.To extract the dominant lengthscales of the canopy coherent motion, we look at the (spatial) spectra of the canopy top surface, reported in Fig. 2 (bottom panels), which is obtained by the surface enveloping the vertical position of the filament tips, Y (x, z,t).The first insight is that a peak can be observed for all cases, and specifically also at the smallest Cauchy number, i.e., Ca = 1.Note that in this condition, often referred to as the "gently swaying" regime, the motion of the filaments is typically claimed to be uncorrelated 31 .Secondly, the location of the peak is typically at O(H), both in the streamwise and spanwise wavelength component, and does not appear to appreciably vary within the investigated range.Increasing Ca, for the spanwise component λ z /H remains around a constant value of 0.6 (only for Ca = 1, it drops around 0.4), whereas a more pronounced variation is observed in the streamwise component with λ x /H approximately ranging between 0.7 and 2. Nevertheless, the absence of a systematic modification in the peak location with Ca suggests that the dominant wavelength is essentially imposed by the flow rather than the structure or that, in other terms, the coherent patterns of the canopy collective deformation can be seen as a signature of the large-scale turbulent structures 11,17 .
To better characterize the relevant lengthscales of the turbulent flow that are involved in the honami/monami phenomenon, Fig. 3 shows the magnitude of the fluid velocity spectra for all components (from top to bottom row) as a function of the streamwise wavelength λ x and wall-normal distance y.Similarly, Fig. 4 reports the same quantities but as a function of the spanwise wavelength λ z .From Fig. 3, a clear peak at λ x ≈ O(H) can be observed for all cases and velocity components.Specifically, for the streamwise component (top panels) the peak is found at λ x /H approximately between 1.2 and 1.8, while for the wall-normal component (middle panels) it is found at λ x /H ≈ 0.6 ± 0.1, and for the spanwise component (bottom panels) at λ x /H ≈ 0.8 ± 0.1.For all these observables, no systematic variation is reported that can be associated with the Cauchy number.Instead, the peaks in the spectra can be associated with the large-scale coherent structures that dominate the collective motion of the canopy.The wall-normal location of the peak appears to gradually decrease while increasing the Cauchy number (from left to right panels), indicating that the dominant turbulent structures remain relatively close to the canopy top when the filaments get bent by the action of the flow (with a maximum offset of about 0.1H observed for the wall-normal and spanwise velocity components at Ca = 100).Note that the observed trend consistently complements with experimental evidence recently reported for specific configurations 5 .Considering the content in the spanwise wavelength, shown in Fig. 4, similar observations can be made.Here, for the streamwise velocity component λ z /H ≈ 0.5 ± 0.1, for the wall-normal component λ z /H ≈ 0.35 ± 0.05 and for the spanwise component λ z /H ≈ 0.8 ± 0.1.Again, the variation in both the streamwise and spanwise wavelength of the peak is rather minimal and does not appear to follow a clear trend.This evidence further confirms that the canopy flexibility (and more generally the structural dynamics) does not have a crucial effect in what concerns the spatial features of the flow-structure interaction.Note that the premultiplied spectra shown in Figs. 3 and 4 are scaled by the bulk velocity U, in . Same as Fig. 3, but as a function of the spanwise wavelength λ z /H and wall-normal coordinate y/H.The grey levels range in: [0, 0.005] with a 0.0005 increment for the wall-normal velocity component; [0, 0.01] with a 0.001 increment for the streamwise and spanwise velocity components.Again, the peak of the spectrum is not substantially affected by the variation of Ca. order to better show the large-scale peak.An alternative normalization, using the local friction velocity to better visualize the flow structures in the intra-canopy region, is provided in Figs.S1 and S2.
A further confirmation of the marginal role of the Cauchy number on the characterization of the honami/monami comes from the joint probability density function (JPDF) of the streamwise and wall-normal velocity fluctuations components, computed on a plane parallel to the wall at the average canopy tip.As Fig. S3 shows, the shape of the JPDF remains substantially the same varying the Cauchy number, with more frequent weak ejections (peak in the second quadrant) and fewer strong sweep events (tails in the forth quadrant).The effect of the Cauchy number can only be noticed with an attenuation of the strong sweep events increasing the flexibility of the canopy stems, in agreement with the literature 16 .

Spatio-temporal properties of honami/monami
Expanding our analysis we now consider not only the spatial, but also the temporal content of the characteristic flow and structural features.Having previously shown that the large-scale spatial content of the flow is not appreciably altered by the canopy flexibility (i.e., Cauchy number), we choose as a representative configuration the one at Ca = 50 where the collective motion is clearly visible (and can be classified as honami/monami regime in the most classical sense).We can now look at the overall picture by comprehensively comparing, in Fig. 5, the spatio-temporal spectra of the flow field (probed in the proximity of the filament tips) and of the canopy tip envelope Y .A low-frequency (approximately ranging between 0.1 and 0.5 U/H, depending on the different velocity components), large-wavelength peak can be generally observed at a very similar location for both the flow and canopy motion.Remarkably, the frequency of the filament tip motion does not appear to agree with the natural frequency (indicated by the violet dashed lines in Fig. 5), thus suggesting the absence, at least for the chosen set of governing parameters, of any fingerprint of structural response in the collective dynamics and, in particular, of peculiar phenomena such as lock-in conditions 15,18,33 .

Individual stem dynamics
Lastly, we contrast the observation on the canopy collective dynamics by focusing on the motion of the individual filaments composing the canopy.Fig. 6 shows the power spectra computed from the time histories of the filament tip velocity (averaged over the different filaments), for both the streamwise (top panels) and spanwise (bottom panels) components.To verify the presence of possible scalings for the dominant peak frequency, the horizontal axis is normalized either with the (first mode) natural frequency f n,1 (left) or the bulk flow timescale H/U (right panels).In fact, it can be noticed that the latter works better 5/16 Figure 6.Temporal spectra of the filament tip velocity (top: streamwise component; bottom: spanwise component) for different Cauchy numbers (increasing with the curve brightness, Ca = 1, 10, 25, 50, 100).Left panels: normalization with the filament natural frequency; right panels: normalization with bulk flow quantities.The former normalization turns out to work better for the spanwise component (i.e., the filaments exhibit their free response behavior), while the streamwise oscillation does not show a similar behavior (except for Ca = 1) and is fully controlled by the turbulent flow.The black dashed lines indicate the slopes f 2 and f −5/3−2 which are obtained by modeling the filament as a harmonic oscillator forced by a fully-developed turbulent flow 32 , while the vertical magenta lines in the left panels indicate the first natural frequency.

6/16
for the streamwise component, while for the spanwise oscillation the dominant frequency scales well with the natural one.This qualitative difference (to the authors' knowledge so far not reported) can be explained by the presence of a net convective transport of the turbulent fluctuations in the streamwise direction.However, it can be observed that when decreasing the Cauchy number the dominant frequency eventually becomes comparable with the natural frequency, and it can be consequently argued that a qualitative change occurs at Ca = 1, with the flapping frequency being now characteristic of the structural response also in the streamwise component.Finally, looking at the whole oscillation spectra, it is worth mentioning the possibility of detecting the fingerprint of turbulence over a certain range of timescales.Indeed, we can model the filament as a harmonic oscillator subject to a fluctuating aerodynamic force that contains the major features of the incoming turbulent flow, and in particular its characteristic spectral content 32,34 .The resulting power spectra associated with the filament oscillation can thus be linked to that of the turbulent forcing by means of a response function dependent on the structural parameters, yielding for the investigated cases the emergence of the two power-law slopes reported in Fig. 6, f 2 and f −5/3−2 , in the lower and higher frequency subranges, respectively.Overall, our numerical results appear in good agreement with such reduced-order model.

Discussion
With the goal of unravelling the essential mechanisms in the mutual interaction between a turbulent flow and a dense hairy surface made of flexible filaments, we have performed high-fidelity direct numerical simulations over a range of the Cauchy number (or, equivalently, of the reduced velocity) in order to isolate the effects related to the canopy flexibility and fluid vs structural characteristic timescales.
Our results show that the so-called honami/monami instability and consequent coherent motion is always observed, also at the small Cauchy numbers to which it is commonly attributed the absence of such physical mechanism.When decreasing the Cauchy number, the deformation of the canopy clearly gets smaller as approaching the rigid limit.Yet, the signature of the large-scale flow structures over the canopy collective motion is still present, therefore suggesting that there is not a particular threshold between different dynamical regimes (e.g., gently swaying vs honami/monami).
The absence of a clear variation in the main features of the flow while varying the Cauchy number similarly supports the idea that the collective dynamics phenomena (i.e., honami/monami) are essentially driven by the resulting turbulent flow, with the flexible canopy having always a passive dynamical behavior in this respect 4,35 .Here, passive means that the coherent waving of the canopy does not manifest the characteristic structural response of the individual canopy elements and is fully driven by the flow.Conversely, certain features of the individual filament motion may (or may not) show the fingerprint of the structural response (i.e., natural frequency).This is typically the case for the individual spanwise oscillation, while for the streamwise oscillation it can be argued only for sufficiently rigid stems.Nevertheless, our conclusions on the fully-passive nature of the honami/monami collective dynamics are not affected by the particular choice of Cauchy number.
In conclusion, our findings provide novel physical insight on the nature of the fascinating collective dynamics that can be observed in a multitude of environmental settings as well as engineered textured surfaces with relevant applications to, e.g., drag reduction, mixing and flow control.Further work aiming at an extension of our results may incorporate the effect of gravity or consider the role of other parameters such as inertia or canopy geometry, although the main driving mechanism of the collective motion can be expected to remain those highlighted in the present investigation.

Direct numerical simulation of turbulent flows
The numerical simulations have been performed using the in-house code Fujin (https://groups.oist.jp/cffu/code).The fluid flow is assumed to be governed by the incompressible Navier-Stokes equations, where u (x,t) and p (x,t) are the velocity and pressure fields, respectively; ρ f is the volumetric fluid density, ν is the kinematic viscosity, and f fs (x,t) is an additional forcing mimicking the presence of the canopy structural elements.We consider a computational domain with size L x /H = 2π, L y /H = 1 and L z /H = 3/2π with periodic boundary conditions along the streamwise and spanwise directions (i.e., x and z), while the no-slip and free-slip conditions are imposed at the bottom (y = 0) and top (y = L y = H) boundaries, respectively 28 .The size of the domain is sufficient to contain the largest, energy-containing scales that are developed in the outer flow 4 .A uniform pressure gradient is imposed in the streamwise direction so that the flow rate is kept constant.Eq. ( 2) are solved numerically using the (second-order) central finite-difference method for the spatial discretization and the (second-order) Adams-Bashforth scheme for advancing in time.We employ a fractional step procedure where the Poisson equation for the pressure is solved using an efficient Fast-Fourier-Transform (FFT) based approach.The solver is parallelized using the MPI protocol and the 2decomp library (http://www.2decomp.org).

Flexible filament and canopy model
The flexible stems composing the canopy are modelled as elastic inextensible filaments by the set of equations where ∆ ρ is the linear density difference between the solid and the fluid, X (s,t) is the position of a generic material point of the filament at the curvilinear coordinate s and time t; T (s,t) is the tension enforcing the inextensibility constraint, F fs (s,t) is the fluid-solid interaction forcing and F col (s,t) is a filament-to-filament or filament-to-wall collision modeling term.The filament's lower end (s = 0) is clamped to the bottom wall while the other end (s = h) is free to oscillate, reflecting into the following set of boundary conditions: X| s=0 = X 0 , ∂ s X| s=0 = ϕ 0 , ∂ ss X| s=h = 0, ∂ sss X| s=h = 0 and T | s=h = 0. From the normal mode analysis for such clamped configuration, the natural frequency of the filament can be evaluated as f nat = (2π) −1 α γ/(∆ ρh 4 ), where α ≈ 3.5160.Both a filament-to-filament and filament-to-wall collision model are implemented which prevent that the stems cross each other or the bottom wall while deforming 36 .However, after extensive testing on different kinds of collision models and their calibration parameters, the influence of the filament-to-filament collision term was found to be very weak, whereas the filament-to-wall interaction model turned out to be necessary only for sufficiently large Ca.The numerical solution of Eq. ( 3) follows the scheme detailed in Ref. 37 with the difference that the bending term is treated implicitly to allow for a larger timestep [38][39][40] .

Simulation setup
In the present study, for all cases the fluid flow is described using n x = 1152, n y = 384 and n z = 864 (Eulerian) grid nodes along the streamwise, wall-normal and spanwise direction, respectively.A nonuniform distribution (with a finer and locally uniform resolution in the lower region containing the canopy) is used in the y-direction in order to better describe the sharper gradients in the proximity of the canopy top.The canopy is made of N x × N z = 15552 flexible stems, with N x = 144 and N z = 108 stems placed along the streamwise and spanwise direction, respectively.Subdividing the horizontal plane in regular tiles of size ∆S = L x /N x = L z /N z , each tile is occupied by a stem that is randomly positioned in order to prevent preferential flow channeling effects 28 .In the rigid canopy case (i.e., Ca = 0), each filament is discretized using N L = 81 (Lagrangian) points such that the spatial resolution ∆s = h/(N L − 1) is approximately equal to the Eulerian grid spacing ∆y in the wall-normal direction.For the flexible canopy case, using such Lagrangian resolution would impose an excessively small timestep when solving Eq. ( 3) and therefore we have set N L = 21, noting that relaxing the requirement on ∆s ≈ ∆y is acceptable when the stems are more compliant because of the decrease in the relative velocity between the stem and the flow.

Validation and convergence study
The numerical method has extensively been tested in a variety of problems [40][41][42][43][44][45] and, in particular, the employed solver for the filament dynamics has been validated in the past with several test cases 40,46 .Results for the representative FSI problem of a flapping filament in laminar flow are reported in Fig. S4.As another validation, we compare our numerical results with available experimental measurements for the case of a rigid-canopy flow at Re = 7070, h/H = 0.65 and λ = 0.83 47 .Results showing good agreement in both the mean velocity profile and Reynolds shear stress distribution are reported in Fig. S5.Finally, considering for the present parametric study the rigid-canopy case (i.e., Ca = 0) where the velocity gradients at the canopy top are more severe, we have assessed the convergence of our results with respect to the spatial resolution along the wall-normal direction ∆y: as shown in Fig. S6, negligible differences are observed when further refining the grid resolution.Note that the more evident variation for the Reynolds shear stress is arguably associated to a lack of statistical convergence rather than for the effect of the numerical resolution.

Figure 1 .
Figure 1.Snapshots of flexible-canopy filaments (left), colored by their vertical position, increasing from light to dark, in order to highlight the coherent waving motion, and with superimposed large-scale flow structures (right), depicted by isosurfaces of positive and negative streamwise velocity fluctuations at u = ±0.4U(in dark and light blue, respectively) in the vicinity of the canopy edge.Results are obtained from the case Ca = 50.

Figure 3 .
Figure 3. Magnitude of the premultiplied spectra of each fluid velocity component 2πΦ u u /(U 2 λ x ), where U is the bulk velocity and u the generic component of the velocity fluctuations, (top: streamwise; middle: wall-normal; bottom: spanwise) as a function of the streamwise wavelength λ x /H and wall-normal coordinate y/H.Results are shown in different columns as a function of the investigated Cauchy number (from left to right, Ca = 0, 1, 10, 25, 50, 100).The white marker denotes the peak of the spectrum, whereas the red horizontal dashed line indicates the averaged height of the filament tips.The grey levels range in: [0, 0.005] with a 0.0005 increment for the streamwise and wall-normal velocity components; [0, 0.01] with a 0.001 increment for the spanwise velocity component.Overall, the characteristic peak in the wavelength of the energy-containing scale appears always O(H) without an appreciable and systematic change associated with Ca.

Figure 5 .
Figure5.Premultiplied spatio-temporal spectra of the fluid velocity, 2π f Φ u u /(U 2 λ x,z ) (where U is the bulk velocity and u the generic component of the velocity fluctuations), probed at the canopy average height (from left to right: streamwise, wall-normal and spanwise component) and canopy top wall-normal displacement (rightmost panel), for a representative case with Ca = 50.Top and bottom panels show the spatial content in the streamwise and spanwise wavelengths, respectively.The vertical red line indicates the wavelength associated to the averaged vertical height of the filament tip, whereas the horizontal violet line corresponds to the first natural frequency of the filament.The white marker denotes the peaks of the spectra.The grey levels range in: [0, 2 × 10 −5 ] with a 2 × 10 −6 increment for the fluid probed at the canopy average height in the streamwise direction; [0, 2 × 10 −4 ] with a 2 × 10 −5 increment for the fluid probed at the canopy average height in the spanwise direction; [0, 10 −6 ] with a 5 × 10 −8 increment for the canopy top wall-normal displacement of the tips (rightmost panels).

Figure S4 .
Figure S4.Validation of the employed simulation procedure for a flapping filament in uniform flow at Re = 200 (for more information on the problem setup, see Huang et al.37 and, in particular, figure13 therein).Comparison of the time history of the filament's trailing point transverse position obtained with Fujin (black line) and that of Huang et al.37 (red circles).Note that the time is normalised with the inflow velocity U and the length of the filament l, while the position of the trailing point is normalised with the length of the filament.

Figure S5 .Figure S6 .
Figure S5.Validation of the employed simulation procedure for the rigid canopy case (R31) experimentally investigated by Shimizu et al.47 .Comparison between our numerical results (black circles) and the experimental measurements (red stars).Left: mean streamwise velocity profile; right: Reynolds shear stress distribution, as a function of the wall-normal location.