1 Introduction

The interaction between elastic structures and fluid flows concerns several problems relevant for both biological and industrial applications, such as, e.g., flow control [1], passive locomotion of swimmers [2,3,4] and energy harvesting from flow-induced vibrations [5,6,7,8,9]. Moreover, fundamental studies have focused on the dynamics of constrained, elastic filaments in laminar flows, investigating the essence of flapping instabilities and related phenomena [3, 10,11,12].

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.

2 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

$$\begin{aligned} \frac{\partial u_i}{\partial t} + \frac{\partial u_i u_j}{\partial x_j}&= - \frac{1}{\rho } \frac{\partial p}{\partial x_i} + {\nu } \frac{\partial ^2 u_i}{\partial x_j x_j} + f_i^\mathrm {T} + f_i^\mathrm {F}, \end{aligned}$$
(1)
$$\begin{aligned} \frac{\partial u_i}{\partial x_i}&= 0, \end{aligned}$$
(2)

where \(u_i\) is the fluid velocity field, p the pressure, \(f_i^\mathrm {T}\) and \(f_i^\mathrm {F}\) two volume body forces used to sustain the turbulent flow and to account for the presence of the immersed fiber, respectively, and \(\rho\) and \(\nu\) the density and kinematic viscosity of the fluid (being \(\mu =\rho \nu\) the dynamic viscosity). The problem can be made non-dimensional by choosing reference length and velocity scales, \(U_*\) and \(L_*\), and by defining the Reynolds number as \(Re=\rho U_* L_*/\mu\). The equations of motion are solved in a triperiodic box, with periodic boundary conditions applied in all the three Cartesian directions.

The forcing \(f_i^\mathrm {T}\) 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_i^\mathrm {F}\) 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_i^\mathrm {F}\) 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

$$\begin{aligned} F_i=\frac{U_i^\mathrm {ib}-U_i^\varGamma }{\varDelta t}, \end{aligned}$$
(3)

where \(U_i^\mathrm {ib}\) is the interpolated fluid velocity on the Lagrangian points which does not satisfy the boundary condition on the immersed objects, \(U_i^\varGamma\) the desired velocity of the Lagrangian points, and \(\varDelta t\) 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.,

$$\begin{aligned} U_i^\mathrm {ib}=\int _V u_i~\delta \left( X_i-x_i \right) ~dV, \end{aligned}$$
(4)

where the integration is performed over the whole fluid domain V. Similarly, the spreading of the fluid–solid interaction force (Eq. (3)) from the Lagrangian grid to the Eulerian one is performed as

$$\begin{aligned} f_i^\mathrm {F}=\rho _1 \int _{c} F_i~\delta \left( X_i - x_i \right) ~ds, \end{aligned}$$
(5)

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]

$$\begin{aligned} \rho _1 \frac{\partial ^2 X_i}{\partial t^2}&= \frac{\partial }{\partial s} \left( T \frac{\partial X_i}{\partial s}\right) - \gamma \frac{\partial ^4 X_i}{\partial s^4} - F_i, \end{aligned}$$
(6)
$$\begin{aligned} \frac{\partial X_i}{\partial s} \frac{\partial X_i}{\partial s}&= 1, \end{aligned}$$
(7)

where T is the tension, \(\rho _1\) the the fiber linear density, \(\gamma\) 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,

figure a

= 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 \(\gamma\) in void, i.e., without fluid, and measure its main oscillation frequency \(f_\mathrm {osc}\). By standard normal mode analysis techniques, we can write that \(f_\mathrm {osc} = \alpha \sqrt{\gamma /(\rho _1 c^4)}\) (where \(\alpha \approx \pi /22.4\)), which can then be used to obtain the actual value of the fiber linear density due to its finite thickness.

2.1 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.1\pi\) 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.

Fig. 1
figure 1

Validation of the hanging filament under gravity: a sketch of the initial condition; b time evolution of the y-coordinate of the free-end of a filament due to gravity (the solid black line represents our numerical results, while the red dots those by Huang et al. [11]); c envelope of filament positions over time. (Color figure online)

Full size image
Fig. 2
figure 2

Validation of filaments in pulsating channel flow: a sketch of the case; b the streamwise displacement of the tip of the last filament at the right end of the row with respect to time. The solid black line is used to show our numerical results, while the red and blue dots experimental and numerical results from the literature [34]. (Color figure online)

Full size image

Next, we validate the fluid-structure interaction solver by considering a pulsating flow in a plane channel filled with glycerine (\(Re=u_\text {m} c/\nu = 60\) being \(u_\text {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.016 c/u_\text {m}\), matching the filaments natural frequency (the filaments Young modulus is \(2.05/\rho u_\text {m}^2\)). 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.

2.2 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.

Fig. 3
figure 3

A snapshot from our DNS of homogeneous isotropic turbulence in which a flexible fiber (line in dark red) is immersed. Isosurfaces of \(\mathcal {Q}\) (in cyan) depict the instantaneous vorticity field, and the three back planes are coloured according to the value of the enstrophy field. (Color figure online)

Full size image

The equations of motion are solved in a triperiodic box with size \(L=2\pi\), 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 \(\frac{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 \(\epsilon\), 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_\lambda =92\). Concerning the discretization of the fiber, the spacing between the Lagrangian points is taken equal to that of the Eulerian grid described before. Depending on the chosen length of the fiber, the number of Lagrangian points used in the simulations ranges from 16 to 41.

The results presented hereafter are obtained by solving the equation of motion and tracking the fibers dynamics for more than 40 large-eddy turnover times, with the statistical quantities averaged over at least 20 large-eddy turnover times.

Fig. 4
figure 4

Time history of the normalized velocity (red) and angular velocity (blue) of a filament translating and rotating in a quiescent fluid. The dashed lines are analytical expressions in the form of \(c_1 e^{t/t_0}\) plotted with the values of \(t_0\) that best fit our data and shifted upwards for major clarity. (Color figure online)

Full size image
Fig. 5
figure 5

a Time histories of the end-to-end distance and b corresponding spectra for the subcritical case \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} = 0.3\) (blue) and the supercritical case \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} = 50\) (green) in the under-damped regime. In b the dominant frequencies are marked with a circle and, for the sake of visualisation, values of the supercritical case are divided by \(10^3\). (Color figure online)

Full size image

3 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 \(\tau _{\gamma }\) by balancing the inertial and bending terms:

$$\begin{aligned} \tau _{\gamma }=\alpha \left( \frac{\rho _1c^4}{\gamma } \right) ^{1/2}, \end{aligned}$$
(8)

where c is the fiber length, \(\rho _1\) the linear density difference, \(\gamma\) the bending rigidity and the factor \(\alpha \approx \pi /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 \(\mathbf{F}\) with a simple viscous term \(\mathbf{F}=-\mu (\dot{\mathbf{X}}- \mathbf{u})\) [38], another characteristic time can be estimated by balancing the inertial term with the introduced damping:

$$\begin{aligned} \tau _{\mu }=\frac{2\rho _1}{\mu }. \end{aligned}$$
(9)

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:

$$\begin{aligned} \zeta =\frac{\tau _{\gamma }}{\tau _{\mu }}=\frac{\alpha c^2\mu }{2\rho _1^{1/2}\gamma ^{1/2}}. \end{aligned}$$
(10)

The critical condition \(\zeta =1\) represents the threshold between two different regimes: for \(\zeta <1\) (under-damped regime) the elasticity is expected to strongly affect the fiber dynamics, while for \(\zeta >1\) (over-damped regime) elastic effects are strongly inhibited. More details will be given in the next section.

Concerning the flow, from the well-known \(\frac{4}{5}\)th Kolmogorov law combined with dimensional analysis, the turnover time of turbulent eddies of size r is expected to behave as

$$\begin{aligned} \tau (r) \sim r^{2/3}{\epsilon }^{-1/3}, \end{aligned}$$
(11)

where \(\epsilon\) 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 \approx 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.

Fig. 6
figure 6

Superposition of instantaneous fiber configurations for the subcritical case \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} = 0.3\) (left) and the supercritical case \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} = 50\) (right) in the under-damped regime (same as Fig. 5). The red circle encloses the undeformed fiber configuration. (Color figure online)

Full size image

4 Fiber flapping states

4.1 Under-damped regime

For \(0<\zeta <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., \(\tau _{\gamma } \approx \tau (c)\), from which a critical value of the bending rigidity can be found:

$$\begin{aligned} \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}}\sim c^{8/3}\epsilon ^{2/3} \rho _1 \alpha ^2. \end{aligned}$$
(12)

Looking at the expression above, a further distinction can be made. In the limit of vanishing \(\gamma\) (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

$$\begin{aligned} St = \frac{t_0 u_0}{c_0}, \end{aligned}$$
(13)

where \(t_0\) is the Stokes time and \(c_0\) and \(u_0\) are characteristic length and velocity scales. Two different length and velocity scales are used: for the translational case the length scale is set equal to the fiber length \(c_0=c\) and the velocity scale is chosen equal to the rms fluid velocity \(u_0 = u'\), while for the rotational Stokes number the length scale is set equal to \(c_0=\pi c\) and the velocity scale to half the velocity difference \(u_0=\delta u_{\parallel }/2\) at the fiber scale. The two Stokes times are measured through numerical simulations where we consider a rotating and translating fiber with speed \(u_0\) in a quiescent fluid with the same viscosity and density used in the rest of the work. We then measure the decay time from our results, as shown in Fig. 4. The values of both these two non-dimensional parameters, computed using the expressions given above, are approximately \(\mathcal {O}(1)\). In particular, the translational Stokes number for \(c/L=0.16\) and \(\rho _1/\left( \rho _0 L^2 \right) = 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 \(\gamma\). 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-to-end 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 \(t u'/c \approx 6\), where \(u'\) 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/\tau _\gamma\), 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_{{\mathrm{turb}}} = 1/\tau (c)\).

Fig. 7
figure 7

The fiber oscillation frequencies (normalized by the turbulence frequency at the fiber length scale) as a function of the fiber bending rigidity (normalized by the corresponding critical value, given by Eq. (12) for the under-damped case (lower axis) and Eq. (15) for the over-damped case (upper axis). Blue: \(c/L=0.12\) and \(\rho _1/(\rho _0 L^2)=0.042\); red: \(c/L=0.16\) and \(\rho _1/(\rho _0 L^2)=0.042\); green: \(c/L=0.20\) and \(\rho _1/(\rho _0 L^2)=0.042\); orange: \(c/L=0.16\) and \(\rho _1/(\rho _0 L^2)=0.014\) or 0.125; gray: \(c/L=0.16\) and \(\rho _1/(\rho _0 L^2)=0\); brown triangles: \(c/L=0.16\) and \(\rho _1/(\rho _0 L^2)=0.042\), passive cases. (Color figure online)

Full size image

To confirm the theoretical predictions, we explore the parametric space by considering fibers with different combinations of density \(\rho _1\), length c and bending rigidity \(\gamma\). Figure 7 shows the ratio between the dominant frequency f and the turbulent frequency \(f_{{\mathrm{turb}}}\) as a function of the ratio between the fiber bending stiffness \(\gamma\) and its corresponding critical value \(\gamma ^{{\mathrm{ud}}}_{{\mathrm{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 \(\gamma =\gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}}\). For \(\gamma < \gamma ^{{\mathrm{ud}}}_{{\mathrm{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 \(\gamma > \gamma ^{{\mathrm{ud}}}_{{\mathrm{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 \(\gamma > \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}}\)) from a soft one occuring for \(\gamma < \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}}\).

Fig. 8
figure 8

Probability density function (PDF) of the longitudinal velocity increments for the subcritical case \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} = 0.3\) (left) and the supercritical case \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} = 50\) (right) in the under-damped regime. Comparison between the Lagrangian fiber measurement (bullets) and the Eulerian one (filled curve)

Full size image

A complementary analysis is performed by looking at the probability density function (PDF) of the longitudinal velocity difference \(\delta u_\parallel \equiv [\mathbf {u}(\mathbf {r},t)-\mathbf {u}(\mathbf {0},t)] \cdot \hat{\mathbf {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).

Fig. 9
figure 9

Second-order (solid line) and third-order (dashed line) velocity structure functions computed in the standard Eulerian way and the same measure obtained by means of Lagrangian fiber tracking for the subcritical, under-damped case with \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} = 0.3\) (blue bullets). Lengths and velocities are made dimensionless with the box size L and with the velocity root mean square, respectively. (Color figure online)

Full size image

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 third-order 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 \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} < 1\).

4.2 Over-damped regime

We now turn our attention to the case where \(\zeta >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:

$$\begin{aligned} \tau _{{\mathrm{re}}} \approx \frac{\mu c^4}{\gamma }. \end{aligned}$$
(14)

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., \(\tau _{{\mathrm{re}}} \approx \tau (c)\), so that the critical value of the bending rigidity for this regime can be estimated as

$$\begin{aligned} \gamma ^{{\mathrm{od}}}_{{\mathrm{crit}}}\sim \mu c^{10/3}\epsilon ^{1/3}. \end{aligned}$$
(15)

We shall discuss the expected behavior in the two limits also here. For \(\gamma /\gamma ^{{\mathrm{od}}}_{{\mathrm{crit}}}<1\), the relaxation is slower than the fluid forcing and thus the fiber is slaved to the turbulence. In the opposite case when \(\gamma /\gamma ^{{\mathrm{od}}}_{{\mathrm{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 \(\gamma /\gamma ^{{\mathrm{od}}}_{{\mathrm{crit}}}\).

Fig. 10
figure 10

Temporal spectrum of the end-to-end distance for the neutrally-buoyant, over-damped case with \(c/L = 0.16\) (gray line). The dominant frequency is marked with a circle

Full size image

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 (\(\rho _1 \sim 0\)), corresponding to \(\zeta \gg 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.

Fig. 11
figure 11

Second-order and third-order velocity structure functions for the neutrally-buoyant, over-damped case (gray bullets) and passive, under-damped cases (brown bullets). Lengths and velocities are made dimensionless with the box size L and with the velocity root mean square, respectively. (Color figure online)

Full size image

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).

Fig. 12
figure 12

Probability density function (PDF) of the velocity increments for the neutrally-buoyant, over-damped case with \(c/L = 0.16\). Comparison between the Lagrangian fiber measurement (bullets) and the Eulerian one (filled curve)

Full size image

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.

Fig. 13
figure 13

Temporal spectra of the end-to-end distance for subcritical case \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} = 0.3\) (brown) and supercritical case \(\gamma / \gamma ^{{\mathrm{ud}}}_{{\mathrm{crit}}} = 50\) (yellow) in the under-damped regime, with deactivation of the fiber-flow feedback (passive cases). The dominant frequencies are marked with a circle. (Color figure online)

Full size image
Fig. 14
figure 14

Probability density function (PDF) of velocity increments for the passive, a subcritical and b supercritical under-damped cases with \(c/L=0.16\). Comparison between the Lagrangian fiber measurement (bullets) and the Eulerian one (filled curve)

Full size image

5 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.

6 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 \(\gamma = EI = 6.76 \times 10^{-7} \mathrm {Nm}^2\) and the linear density \(\rho _1\) is \(0.86 \mathrm {g/m}\), for their Silicone II \(\gamma = EI = 1.74 \times 10^{-6} \mathrm {Nm}^2\) and \(\rho _1 = 0.41 \mathrm {g/m}\), and finally, for their Nylon (III) \(\gamma = EI = 2.75 \times 10^{-5} \mathrm {Nm}^2\) and \(\rho _1 = 0.16 \mathrm {g/m}\). We can also refer to the silk flexible filament used in Ref. [10]: it is a \(0.15 \mathrm {mm}\) diameter silk wire (Silk in the following) having \(\gamma = EI = 10^{-10} \mathrm {Nm}^2\) and linear density \(\rho _1 = 0.02 \mathrm {g/m}\). We can estimate the critical fiber length \(c_{{\mathrm{crit}}}\) to identify the threshold between the under-damped (\(c < c_{{\mathrm{crit}}}\)) and the over-damped regimes (\(c > c_{{\mathrm{crit}}}\)). Considering water as the working fluid, we have for Silicone I \(c_{{\mathrm{crit}}} \approx 0.15 \mathrm {m}\), for Silicone II \(c_{{\mathrm{crit}}} \approx 0.16 \mathrm {m}\), for Nylon III \(c_{{\mathrm{crit}}} \approx 0.26 \mathrm {m}\), and for Silk \(c_{{\mathrm{crit}}} \approx 0.0067 \mathrm {m}\). 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 \(\gamma < \gamma ^{{\mathrm{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 \(91 \mathrm {\mu m}\). This implies that, in water, the largest \(\epsilon\) is approximately \(50 \mathrm {m^2 /s^3}\). Assuming \(c = 13 \mathrm {cm}\) (that seems to be within the inertial range of the experiment) and considering the Silicone I fibers we get: \(\gamma ^{{\mathrm{ud}}} \approx 9.95 \times 10^{-7} \mathrm {Nm^2} > 6.76 \times 10^{-7} \mathrm {Nm^2}\), yielding a ratio \(\gamma /\gamma ^{{\mathrm{ud}}} \approx 0.68\). In 3D turbulence the under-damped regime is thus accessible.

7 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 simulations employing an immersed boundary technique for the fluid-structure interaction.

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.