Flowing fibers as a proxy of turbulence statistics

The flapping states of a flexible fiber fully coupled to a three-dimensional turbulent flow are investigated via state-of-the-art numerical methods. Two distinct flapping regimes are predicted by the phenomenological theory recently proposed by Rosti et al. [Phys. Rev. Lett. 121, 044501, 2018]: the under-damped regime, where the elasticity strongly affects the fiber dynamics, and the over-damped regime, where the elastic effects are strongly inhibited. In both cases we can identify a critical value of the bending rigidity of the fiber by a resonance condition, which further provides a distinction between different flapping behaviors, especially in the under-damped case. We validate the theory by means of direct numerical simulations and find that, both for the over-damped regime and for the under-damped one, fibers are effectively slaved to the turbulent fluctuations and can therefore be used as a proxy to measure various two-point statistics of turbulence. Finally, we show that this holds true also in the case of a passive fiber, without any feedback force on the fluid.

On the other hand, in the field of particle-laden flows attention has been devoted to the dynamics of fiber-like objects dispersed in laminar or turbulent flows. In this latter framework, a further distinction can be made depending on the size of the fiber compared to that of the existing flow scales. A significant number of contributions have considered fibers with length smaller than the Kolmogorov scale [13,14], revealing the nature of effects such as preferential alignment for neutrally-buoyant rods [15] and preferential concentration when their inertia enters into play [16,17].
Fewer studies have addressed the case of finite-size fibers whose length is comparable to scales belonging to the inertial range. Among these, experimental analyses have outlined the potential of using dispersed rigid rods as a measurement tool for turbulent flows [18], since their tumbling rate is found to be approximately equal to the characteristic time of turbulent structures of comparable size (provided inertia does not affect the fiber dynamics [19,20]).
The case of flexible fibers has been recently addressed both experimentally [21,22] and numerically [23,24]. One of the main findings of Ref. [24] is the identification of a flapping regime where the fiber deformation is slaved to turbulent fluctuations, enabling to quantify their statistical properties by measuring only the distance and velocity difference between the fiber free ends.
Further attempts of such a Lagrangian description have seen the employment of other kinds of particles for evaluating both two-point (limited to distances between particles smaller than the Kolmogorov viscous scale) and single-point quantities [25,26], paving the way for new strategies to investigate turbulent flows. Overall, such efforts are needed in order to increase our understanding of turbulence and establish, in particular, a connection between scaling laws and spatial structures, e.g., vortex filaments [27,28]. In passive scalar turbulence, this connection leads to a complete understanding of the meaning of intermittency and anomalous scaling [29,30].
The goal of this work is to present an exhaustive description of the dynamical phenomenology associated to a long flexible fiber (i.e., having a rest length falling within the turbulence inertial range of scales) freely moving in a controlled turbulent flow. More specifically, we aim to categorize the different regimes that can be predicted theoretically by combining a simple structural model for the fiber with a widely studied turbulence model: the so-called homogeneous, isotropic and stationary turbulence model. We will present and discuss the different flapping states that may occur depending on the choice of the characteristic parameters, supporting our physical intuitions with evidence from fully-resolved numerical simulations. Furthermore, the underlying hypothesis of passive fiber, on which the phenomenological model that will be introduced in this work relies, will be reviewed by considering numerical experiments in which the feedback exerted by the fiber on the flow is deactivated; This will clarify whether this is crucial to capture the essential dynamics.
The idea of using a fiber to measure two-point turbulence statistics was recently proposed by Rosti et al. [24]. Our goal here is to give a detailed presentation on how to exploit a flexible fiber as a proxy of turbulence statistics and to present new results.
The present work is structured as follows: Sect. 2 presents the numerical method along with its validation. In Sect. 3 we introduce our phenomenological model while in Sect. 4 we discuss the different dynamical regimes and provide corroborations from DNS. Section 5 concerns the case of passive fiber and in Sect. 6 we propose parameters for possible experiments. Finally, Sect. 7 draws some conclusions.

Numerical method
We consider the fully coupled dynamics of a flexible fiber, governed by the Euler-Bernoulli equation, immersed in an incompressible three-dimensional turbulent flow field, governed by the Navier-Stokes equations.
In an inertial, Cartesian frame of reference the equations of momentum and mass conservation for the incompressible flow read as where u i is the fluid velocity field, p the pressure, f T i and f F i two volume body forces used to sustain the turbulent flow and to account for the presence of the immersed fiber, respectively, and q and m the density and kinematic viscosity of the fluid (being l ¼ qm the dynamic viscosity). The problem can be made nondimensional by choosing reference length and velocity scales, U Ã and L Ã , and by defining the Reynolds number as Re ¼ qU Ã L Ã =l. The equations of motion are solved in a triperiodic box, with periodic boundary conditions applied in all the three Cartesian directions.
The forcing f T i is used to generate and sustain a fully turbulent flow with homogeneous, isotropic, and stationary statistics. To do so, we use the spectral forcing scheme proposed by Eswaran and Pope [31], which involves the addition of energy to the Fourier modes of the velocity at wavenumbers inside a low wavenumber shell. The injected energy is obtained from a formulation based on Uhlenbeck-Ornstein processes and the resulting flowfield displays neither anisotropy nor unreasonably high correlation times [31].
The fluid-solid coupling force f F i is obtained by the Immersed Boundary Method (IBM), a technique used to simulate the flow past solid bodies, first developed by Peskin [32] to simulate blood flow inside a heart. The main feature of this method is that the numerical grid does not need to conform to the geometry of the object, which is replaced by the body force distribution f F i which mimics the presence of the body on the fluid and restores the desired velocity boundary conditions on the immersed surfaces. Two sets of grid points are needed: a fixed Eulerian grid x i for the fluid and a moving Lagrangian grid X i for the structure. The body force is found by first computing the fluid-solid interaction force as where U ib i is the interpolated fluid velocity on the Lagrangian points which does not satisfy the boundary condition on the immersed objects, U C i the desired velocity of the Lagrangian points, and Dt the time step. The interpolation from the Eulerian grid to the Lagrangian one of the fluid velocity is performed using a smooth Delta function, i.e., where the integration is performed over the whole fluid domain V. Similarly, the spreading of the fluidsolid interaction force (Eq. (3)) from the Lagrangian grid to the Eulerian one is performed as where s is the curvilinear coordinate along the fiber. The update of the Lagrangian points is achieved by solving a separate equation that describes the dynamics of the flexible fiber; in our simulations we use the Euler-Bernoulli beam equation together with the inextensibility condition [33] oX i os where T is the tension, q 1 the the fiber linear density, c the bending rigidity, and F i the fluid-solid interaction force. This model is fully justified as far as the ratio between the filament length and its diameter, = c/d, is much larger than unity. The fiber is free to move in the flow, thus, zero force, torque and tension boundary conditions are enforced at the two ends of the fiber. Gravitational effects are neglected, i.e., the Froude number is always much larger than unity. The fluid equations are solved numerically on a staggered grid, with pressure points located at the cell center and velocity components at the cell faces, using a second order finite difference code. Equations (1) and (2) are advanced in time by a fully explicit fractional step-method, where the third-order Runge-Kutta method is used, and the Poisson equation is solved by Fast Fourier Transform. To solve the fiber equation we follow the explicit two-step method proposed by Huang et al. [11].
Note that, the exact relation between the fiber volume and linear density is not clearly defined in the method described above due to the uncertain definition of the shape and cross-section of the fiber, which, although being mono-dimensional in the theory, has a finite thickness due to the spreading operation (Eq. (5)), with the Dirac-delta function having a support of 4 grid points. The problem can be solved as follows: we first simulate a free fiber with a prescribed bending rigidity c in void, i.e., without fluid, and measure its main oscillation frequency f osc . By standard normal mode analysis techniques, we can write that f osc ¼ a ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c=ðq 1 c 4 Þ p (where a % p=22:4), which can then be used to obtain the actual value of the fiber linear density due to its finite thickness.

Validation
The numerical code used in this work has been extensively validated in the past for turbulent flow simulations [35][36][37]. Here, we provide one more comparison with literature results for the specific case of a flexible filament.
First, we validate the structural solver by studying the oscillations of a hanging filament under gravity, as done by Huang et al. [11]. The filament is initially held stationary with an angle of 0:1p from the vertical (Fig. 1a) and then, after being released, starts swinging due to the gravity force. The unit filament is discretised in our simulation with 100 Lagrangian points. Figure 1b shows the time history of the free-end position obtained by our simulation (solid line) and by Huang et al. [11] (red dots); a good agreement with the literature data is evident.
Next, we validate the fluid-structure interaction solver by considering a pulsating flow in a plane channel filled with glycerine (Re ¼ u m c=m ¼ 60 being u m the maximum velocity). A row of 5 filaments is hinged vertically on the bottom wall (Fig. 2a). The pulsating frequency of the channel is 0:016c=u m , matching the filaments natural frequency (the filaments Young modulus is 2:05=qu 2 m ). Figure 2b shows the streamwise displacement of the tip of the last filament with respect to time: our results (solid line) are compared with the experimental measurements (red dots) and with the simulations (blue dots) reported by Pinelli et al. [34]. Both the frequency of oscillation and the magnitude of the displacement match the literature results.

DNS setup
Details are given here on the numerical simulations whose data are reported in the remaining part of the paper. We consider a flexible fiber immersed in an incompressible three-dimensional turbulent flow field; in Fig. 3 we report an instantaneous visualization of the turbulent flow along with the dispersed fiber to give a qualitative insight of the resulting scenario.
The equations of motion are solved in a triperiodic box with size L ¼ 2p, discretized on a Cartesian uniform mesh using 128 grid points per side. With this grid resolution, the resulting turbulence two-point statistics is consistent with the 4 5 th Kolmogorov law in the inertial range of scales (as one can notice for the Eulerian quantities in Fig. 11), with negligible differences if doubling the number of nodes in all directions. The resulting turbulent dissipation rate , made dimensionless with the cube of the velocity root-mean square divided by the size of the box, is about 2.6 and the Reynolds number at the Taylor microscale is about Re k ¼ 92. Concerning the discretization of the fiber, the spacing between the Lagrangian points is taken equal to that of the Eulerian grid described before. dynamics for more than 40 large-eddy turnover times, with the statistical quantities averaged over at least 20 large-eddy turnover times.

Phenomenological model
In this section, we start by reviewing the arguments presented in Ref. [24]. Our goal is to model the coupling between the fiber and the flow in a simple, yet effective way. For this purpose, we begin by estimating the different timescales of the problem.
Recalling the fiber dynamical equation (6) and neglecting at first the forcing from the flow, one can define the fiber elastic time s c by balancing the inertial and bending terms: where c is the fiber length, q 1 the linear density difference, c the bending rigidity and the factor a % p=22:4 results from a normal mode analysis, corresponding to the first natural frequency of the beam in the unsupported (free-free) case. If we model now the fluid-solid interaction force F with a simple viscous term F ¼ Àlð _ X À uÞ [38], another characteristic time can be estimated by balancing the inertial term with the introduced damping: It has to be noted that the expression chosen for the drag does not account for anisotropic effects (as done, e.g., for fibers in low-Reynolds flows [39]), since in our case the flow has an essentially isotropic nature and no preferential alignment is expected. An analogy can be drawn between our model and a damped harmonic oscillator, so that we can derive an expression for the equivalent damping ratio: The critical condition f ¼ 1 represents the threshold between two different regimes: for f\1 (underdamped regime) the elasticity is expected to strongly affect the fiber dynamics, while for f [ 1 (overdamped regime) elastic effects are strongly inhibited. More details will be given in the next section.
Concerning the flow, from the well-known 4 5 th Kolmogorov law combined with dimensional analysis, the turnover time of turbulent eddies of size r is expected to behave as where is the turbulent dissipation rate of kinetic energy. Let us now assume that the flow structures (i.e., eddies) effectively acting in the fluid-structure coupling are those with the same size of the fiber, i.e., r % c. This also implies that the activated oscillation mode of the fiber is essentially the first one (whose associated timescale is supplied by Eq. (8)). Having collected these relations, we are able to make predictions regarding the flapping states of the elastic structure depending on the physical parameters.

Under-damped regime
For 0\f\1 (under-damped regime), we expect that the fiber response shows an oscillatory behavior. Therefore, we impose a resonance condition between the elastic time and the hydrodynamic one, i.e., s c % sðcÞ, from which a critical value of the bending rigidity can be found: Looking at the expression above, a further distinction can be made. In the limit of vanishing c (subcritical case), the fiber can be thought to be slaved to the flow due to the relatively faster forcing compared with its response, therefore flapping at the eddy frequency. Such a slaved dynamics is expected to hold true when both the rotational and translational Stokes number are sufficiently small. The rotational and translational Stokes numbers are defined as computed using the expressions given above, are approximately Oð1Þ. In particular, the translational Stokes number for c=L ¼ 0:16 and q 1 = q 0 L 2 ð Þ¼0:042 is equal to 7 while the rotational one equals 0.3. The fact that the translational Stokes number is not smaller than unity, as it is for the rotational Stokes number, suggests that the possibility of using a fiber as a proxy of two-point statistics of turbulence is controlled by the (smallness of the) rotational Stokes number rather than by the translational Stokes number. Further analysis of this issue are however worth investigating to arrive at a firm conclusion.
In the opposite supercritical case, where the elastic time is much smaller than the hydrodynamic one, the fiber reaction is expected to be far more rapid than the fluid forcing.
To corroborate these expectations, we consider results from DNS of two cases, both belonging to this under-damped regime but with different c. We start by looking at how the fiber end-to-end distance varies in time (Fig. 5a). The two curves look clearly different: in the subcritical case, finite-size and continuous variations of the signal are found, with a mean end-toend value of around 0.55; conversely, in the supercritical case, the fiber remains almost unbent (the mean end-to-end distance is about 0.99) and isolated bursts are observed (such as that occurring from tu 0 =c % 6, where u 0 is the velocity root-mean square), when interactions occur with energetic eddies. These observations are confirmed by the pictorial views in Fig. 6 where the overlapped fiber configurations at different time instants are collected: we observe appreciable deformations in the subcritical case, while the fiber essentially behaves like a rigid rod in the supercritical regime.
For a deeper understanding, the corresponding temporal spectra are analysed and reported in Fig. 5b. A substantial difference can be noticed between the two cases: in the supercritical case, the dominant frequency is the natural frequency of the fiber, i.e., f ¼ 1=s c , with a clear peak in the spectra, while for the subcritical one the peak of the spectrum is less distinct and found at the turbulent frequency f turb ¼ 1=sðcÞ.
To confirm the theoretical predictions, we explore the parametric space by considering fibers with different combinations of density q 1 , length c and bending rigidity c. Figure 7 shows the ratio between the dominant frequency f and the turbulent frequency f turb as a function of the ratio between the fiber bending stiffness c and its corresponding critical value c ud crit . As one can observe, two clearly different regimes are found, defining two distinct flapping states of the fiber, with a well defined threshold around c ¼ c ud crit . For c\c ud crit , i.e., the subcritical cases, all the points lay on a horizontal line characterized by an oscillation frequency equal to the turbulent one; on the contrary, for c [ c ud crit , supercritical regime, all the points collapse on an inclined line with slope 0.5, indicating that the oscillation frequencies grow with the bending rigidity and are larger than the turbulent frequency. These findings thus confirm the idea that the subcritical case is fully governed by turbulence, unlike the supercritical one which exhibits a structural response behavior.
The net demarcation between the two regimes also allows us to quantitatively separate a hard regime of oscillation (for c [ c ud crit ) from a soft one occuring for c\c ud crit . A complementary analysis is performed by looking at the probability density function (PDF) of the longitudinal velocity difference du k ½uðr; tÞ À uð0; tÞ Ár sampled at the free ends of the fiber, compared with the same quantity computed in the Eulerian frame. Results for the two fibers considered previously are shown in Fig. 8: the PDF of the Eulerian data shows a non-symmetric bell shape, while two very different curves are found for the two fibers. Indeed, while we notice a good agreement in the subcritical case between the Eulerian and Lagrangian data, a significative difference is found for the supercritical one. This result supports one more time the idea that the fiber dynamics in the subcritical regime is dominated by turbulence, while in the supercritical one the fiber response comes from its structural elasticity. Concerning the first case, we ascribe the small mismatch occurring between the Lagrangian and Eulerian measurements to the different number of samples (about two order of magnitudes lower for the former).
Focusing on the subcritical situation and aiming at exploiting the capability of flexible fibers acting as a proxy of turbulent eddies, the consequent step is to use fibers of different length to obtain the flow velocity structure functions and estimate the associated scaling laws. To this end, Fig. 9 shows the second and thirdorder longitudinal velocity structure functions S p ðrÞ (with p ¼ 2; 3), proposing a comparison between the two approaches analogous to what was presented for the PDF. Notice that the plot abscissa refers to the time-averaged value of the end-to-end distance and not to the fiber rest length, this quantity being more representative of the actual fiber length scale (see Fig. 6a). For both the reported quantities, the Lagrangian and Eulerian measurements reveal a close resemblance with differences well within the statistical error. This leads to the conclusion that, in the under-damped regime, it is possible to measure the two-point statistics of the flow (e.g., PDF, structure functions) by means of dispersed flexible fibers tracked in time, provided that c=c ud crit \1.

Over-damped regime
We now turn our attention to the case where f [ 1 (over-damped regime), where dissipation becomes dominant. Once deformed, the fiber response is characterized by a relaxation timescale that can be estimated by balancing the elastic and viscous terms, yielding: Note that this relaxation process has exponential behavior, without oscillations. Indeed, the fiber equation becomes first-order in time.
In this case, a balance condition can be imposed between the relaxation time and the eddy turnover one, i.e., s re % sðcÞ, so that the critical value of the bending rigidity for this regime can be estimated as We shall discuss the expected behavior in the two limits also here. For c=c od crit \1, the relaxation is slower than the fluid forcing and thus the fiber is slaved to the turbulence. In the opposite case when c=c od crit [ 1, the fiber is appreciably deformed only by large strains, similarly to the under-damped regime. However, here elastic oscillations are not possible, and the dominant frequency in this case is expected to be the turbulent one. The fiber motion in the over-damped regime is therefore always expected to be slaved to turbulence, independently of c=c od crit . To prove this, we resort again to numerical experiments. When dealing with the subcritical, over-damped case, however, issues arise since the fibers are excessively flexible, leading to numerical instability. Therefore we consider only the zero-mass case (q 1 $ 0), corresponding to f ) 1. As done for the under-damped regime, the time trace of the end-to-end displacement is acquired and processed, extracting the peak of the frequency as shown in Fig. 10. These data are reported in Fig. 7 (gray symbols), to show that for all computed cases the frequency ratio is approximately 1, demonstrating that the fiber flapping is locked-in to the flow.
The comparison in terms of two points turbulent statistic presented for the under-damped regime is repeated for this regime as well. First, we consider the computed velocity structure functions and the results are shown for one representative case in Fig. 11 (along with results from the passive cases that will be introduced in Sect. 5).
Further, Fig. 12 depicts the resulting PDFs for the same case. Both observables measured by the Lagrangian fiber are in good agreement with the Eulerian data.
As shown for the under-damped regime, also the predictions for the over-damped case are confirmed by the results from numerical simulations.

Passive fiber
In the numerical framework considered so far, the presence of the fiber modifies locally the flow. In our phenomenological model, on the other hand, we have assumed that the former has an essentially passive behavior. The question that rises is therefore: are our findings confirmed if the fiber-flow feedback in the numerical method is deactivated? To address this point, simulations are performed neglecting the feedback to the flow.
The obtained picture is the same of what found for fibers with active feedback. The spectra from which the main flapping frequency was extracted are reported in Fig. 13. The data are added in Fig. 7 with brown symbols for the subcritical cases, yellow for the supercritical ones, showing that the fiber dynamics is consistent with the full model. Focusing on the subcritical case, we examine also the velocity structure functions (Fig. 11, brown bullets), for which again a good agreement with the Eulerian counterpart is found. Lastly, the PDF of the longitudinal velocity difference for c=L ¼ 0:16 is presented in Fig. 14, again confirming the conclusions above. In light of these evidences, it appears that the action of the fiber on the flow can be neglected when modeling the dynamics of single fibers in turbulence.

Suggestions for possible experiments
In this section we propose some estimations for planning laboratory experiments. To identify some possible materials for the fiber, we can refer to Ref. [21]. For their Silicone I, the resulting bending rigidity is c ¼ EI ¼ 6:76 Â 10 À7 Nm 2 and the linear density q 1 is 0:86g=m, for their Silicone II c ¼ EI ¼ 1:74 Â 10 À6 Nm 2 and q 1 ¼ 0:41g=m, and finally, for their Nylon (III) c ¼ EI ¼ 2:75 Â 10 À5 Nm 2 and q 1 ¼ 0:16g=m. We can also refer to the silk flexible filament used in Ref. [10]: it is a 0:15mm diameter silk wire (Silk in the following) having c ¼ EI ¼ 10 À10 Nm 2 and linear density q 1 ¼ 0:02g=m. We can estimate the critical fiber length c crit to identify the threshold between the under-damped (c\c crit ) and the over-damped regimes (c [ c crit ). Considering water as the working fluid, we have for Silicone I c crit % 0:15m, for Silicone II c crit % 0:16m, for Nylon III c crit % 0:26m, and for Silk c crit % 0:0067m. For the first 3 cases, fibers with lengths of a few centimeters (say, up to 10 cm) thus fall in the under-damped regime. For the Silk, with a centimeter-size fiber one falls in the over-damped regime, hence the fiber is expected to behave as a proxy of turbulence. Finally, for the first three fibers, we now need to estimate whether with such mechanical properties of the fiber we are in the region c\c ud . To estimate whether or not this is the case, we need to know something on the turbulence field: in Ref. [21], their Kolmogorov scale ranges between 12 and 91lm. This implies that, in water, the largest is approximately 50m 2 =s 3 . Assuming c ¼ 13cm (that seems to be within the inertial range of the experiment) and considering the Silicone I fibers we get: c ud % 9:95 Â 10 À7 Nm 2 [ 6:76 Â 10 À7 Nm 2 , yielding a ratio c=c ud % 0:68. In 3D turbulence the under-damped regime is thus accessible.

Conclusions
This study concerns the dynamics of flexible fibers dispersed in homogeneous, isotropic turbulence. Based on simple resonance conditions between different characteristic timescales, we proposed a phenomenological model able to classify the flapping regimes experienced by the fiber. The predictions have been corroborated by fully-coupled direct numerical For those regimes fully slaved to the flow, the fiber may be viewed as a Lagrangian tracker of turbulent eddies, exploitable for evaluating not only their characteristic time but also two-point statistical quantities such as, e.g., scaling exponents of velocity structure functions. We believe that this concept has potential applications in experimental methods for turbulence measurement paving the way to a new experimental strategy, a sort of ''Fiber Image Velocimetry'', suitable to study turbulence two-point statistics or multipoint statistics (once a single fiber is replaced by more complex elastic structures). A substantial difference between the principle of ''Particle Image Velocimetry'' and the one of the proposed ''Fiber Image Velocimetry'' is that while the relative distance between a pair of tracer particles is not maintained constant in time due to the celebrated Richardson law (the relative distance between particles grows in time as t 3=2 ), on the other hand, this requirement is intrinsically satisfied when considering an inextensible fiber.
While in this investigation we considered only the behavior of a single, isolated fiber, the described strategy can be employed seeding the flow with several elastic objects, potentially of different lengths in order to measure eddies of different size. Increasing the concentration of the dispersed phase, however, would determine an increase of the importance of the fiber-flow feedback, so that the assumption leading to the prediction for the critical values of the bending rigidity must be modified to account for the modulation of the flow statistics caused by the fiber feed-back.