1. Introduction
Suspensions of fibre-like objects dispersed in turbulent flows are encountered in a variety of natural and engineering problems, e.g. microplastics or non-motile microorganisms dispersal in aquatic environments, pulp production in papermaking, and other ecological or industrial processes (Lundell, Söderberg & Alfredsson Reference Lundell, Söderberg and Alfredsson2011; Guasto, Rusconi & Stocker Reference Guasto, Rusconi and Stocker2012; du Roure et al. Reference du Roure, Lindner, Nazockdast and Shelley2019). Several topics appear of particular relevance nowadays, such as the quantitative understanding of the fragmentation process of marine litter in the oceans (Cózar et al. Reference Cózar2014; Filella Reference Filella2015; Brouzet et al. Reference Brouzet, Guiné, Dalbe, Favier, Vandenberghe, Villermaux and Verhille2021). Other examples concern the locomotion of bacteria or planktonic organisms (Son, Guasto & Stocker Reference Son, Guasto and Stocker2013; Ardekani et al. Reference Ardekani, Sardina, Brandt, Karp-Boss, Bearon and Variano2017b; Michalec et al. Reference Michalec, Fouxon, Souissi and Holzner2017), the formation of algae aggregates (Verhille et al. Reference Verhille, Moulinet, Vandenberghe, Adda-Bedia and Le Gal2017), and the manufacturing of paper and composite materials (Parsheh, Brown & Aidun Reference Parsheh, Brown and Aidun2005; Mortensen et al. Reference Mortensen, Andersson, Gillissen and Boersma2008; Marchioli, Fantoni & Soldati Reference Marchioli, Fantoni and Soldati2010).
From a fundamental perspective, understanding the most salient features in this specific kind of multiphase flow represents a topic of active research. Compared to spherical particles, additional complexity is given by the anisotropic character of the dispersed objects, reflecting into peculiar phenomena such as the preferential alignment of the fibre's orientation with characteristic quantities of the flow, i.e. vorticity or strain rate principal directions (Parsa et al. Reference Parsa, Calzavarini, Toschi and Voth2012; Voth & Soldati Reference Voth and Soldati2017). Focusing on rigid fibres (i.e. rods), numerous efforts have been devoted to describe in detail how such alignment, and consequently the fibre's rotation rate, are qualitatively and quantitatively conditioned by the fibre's length (compared to the characteristic lengthscales of the turbulent flow, ranging from sub-Kolmogorov to inertial subrange) (Ni, Ouellette & Voth Reference Ni, Ouellette and Voth2014; Ni et al. Reference Ni, Kramel, Ouellette and Voth2015; Pujara, Voth & Variano Reference Pujara, Voth and Variano2019; Oehmke et al. Reference Oehmke, Bordoloi, Variano and Verhille2021; Pujara et al. Reference Pujara, Arguedas-Leiva, Lalescu, Bramas and Wilczek2021) and inertia (expressed by means of a representative Stokes number) (Bounoua, Bouchet & Verhille Reference Bounoua, Bouchet and Verhille2018; Kuperman, Sabban & van Hout Reference Kuperman, Sabban and van Hout2019).
In many applications, however, fibres are not rigid but flexible, and may experience large deformations under the action of the flow, resulting in non-trivial dynamics already when immersed in simple laminar flows (du Roure et al. Reference du Roure, Lindner, Nazockdast and Shelley2019; Żuk et al. Reference Żuk, Słowicka, Ekiel-Jeżewska and Stone2021). While a relatively extensive research has focused on suspensions of rigid fibres, fewer studies are available on the case of flexible fibres in turbulence. Besides the relevance in the framework of particle-laden and multiphase flows, the latter represents an intriguing problem in the context of fluid-structure interaction (FSI), started to be investigated only quite recently (Brouzet, Verhille & Le Gal Reference Brouzet, Verhille and Le Gal2014; Rosti et al. Reference Rosti, Banaei, Brandt and Mazzino2018). From the theoretical and numerical viewpoint, fibres are typically assumed to be inextensible, so that the attention is primarily devoted to the bending deformation (Allende, Henry & Bec Reference Allende, Henry and Bec2020), although some contributions focused on extensible objects as well, mostly related to suspensions of polymers (Picardo et al. Reference Picardo, Vincenzi, Pal and Ray2018, Reference Picardo, Singh, Ray and Vincenzi2020; Vincenzi et al. Reference Vincenzi, Watanabe, Ray and Picardo2021). The mechanical behaviour of flexible fibres in turbulence has been characterised with specific interest towards the buckling of short fibres, i.e. having sub-Kolmogorov length (Allende, Henry & Bec Reference Allende, Henry and Bec2018), and the spatial conformation and deformation of finite-size fibres, i.e. with length within the inertial subrange or comparable to the integral length scale (Brouzet et al. Reference Brouzet, Verhille and Le Gal2014; Gay, Favier & Verhille Reference Gay, Favier and Verhille2018; Dotto & Marchioli Reference Dotto and Marchioli2019; Sulaiman et al. Reference Sulaiman, Climent, Delmotte, Fede, Plouraboué and Verhille2019; Dotto, Soldati & Marchioli Reference Dotto, Soldati and Marchioli2020; Picardo et al. Reference Picardo, Singh, Ray and Vincenzi2020); further recent developments include the aforementioned modelling of fragmentation processes (Allende et al. Reference Allende, Henry and Bec2020; Brouzet et al. Reference Brouzet, Guiné, Dalbe, Favier, Vandenberghe, Villermaux and Verhille2021; Vincenzi et al. Reference Vincenzi, Watanabe, Ray and Picardo2021).
From a relatively complementary perspective, recent studies focused on the possibility of exploiting fibre-like objects to measure the properties of the turbulent carrier flow, both in the case of flexible (Rosti et al. Reference Rosti, Banaei, Brandt and Mazzino2018, Reference Rosti, Olivieri, Banaei, Brandt and Mazzino2020a) and rigid fibres (Brizzolara et al. Reference Brizzolara, Rosti, Olivieri, Brandt, Holzner and Mazzino2021). For flexible fibres, Rosti et al. (Reference Rosti, Banaei, Brandt and Mazzino2018, Reference Rosti, Olivieri, Banaei, Brandt and Mazzino2020a) unravelled the existence of different dynamical states that can be predicted by comparing the characteristic timescales involved in the problem. In particular, it was shown that in certain regimes the fibres behave as a proxy of turbulent eddies of comparable size, therefore enabling the measurement of two-point flow statistics based on the longitudinal velocity differences by simply tracking the motion of the fibre (or, more practically, only the two fibre's ends). More recently, Olivieri, Mazzino & Rosti (Reference Olivieri, Mazzino and Rosti2021) extended such findings to the case of a non-dilute suspension where the flow is substantially altered by the presence of the dispersed phase, showing that the same qualitative scenario is retained (although in the latter situation fibres do not measure the unperturbed flow anymore, but the fibre modulated one). In the case of rigid fibres, similar outcomes can be obtained when considering negligible inertia (i.e. vanishing Stokes number) and focusing on the transverse (instead of longitudinal) velocity differences: such findings have been recently corroborated in the laboratory framework by Brizzolara et al. (Reference Brizzolara, Rosti, Olivieri, Brandt, Holzner and Mazzino2021), leading to the development of a novel experimental technique, named ‘Fibre Tracking Velocimetry’, able to measure the properties of turbulence at a fixed length scale (i.e. the length of the fibre). Remarkably, exploiting this concept intrinsically overcomes the well-known issue of relative dispersion experienced with more traditional approaches based on tracer particles. Moreover, it was demonstrated that for sufficiently short fibres the proposed technique leads to the accurate measurement of the turbulent energy dissipation rate.
The majority of studies focusing on the dynamics of dispersed particles in turbulence are typically based on assuming that the suspension is dilute enough, so that the backreaction of the dispersed phase to the carrier flow can be safely ignored, and therefore dramatically simplifying the modelling of the problem by exploiting a one-way coupling approach (i.e. the fluid flow is not affected by the presence of the dispersed objects). However, one needs to consider the two-way coupled problem for relatively high concentrations, i.e. not only the fibres are transported and deformed by the flow, but the flow in turns gets altered by their backreaction due to the no-slip condition at the surface of the immersed objects (Balachandar & Eaton Reference Balachandar and Eaton2010; Maxey Reference Maxey2017; Brandt & Coletti Reference Brandt and Coletti2021). Modelling such dynamical feedback by means of suitable techniques, e.g. immersed boundary methods, poses significant demands which are becoming more affordable with the availability of more powerful computational resources. Besides, along with ensuring the mutual coupling between the two phases, such fully-resolved approaches are expected to give more accurate results compared with traditional models used for describing the dynamics of the dispersed phase, whose foundation rely on one-way coupling and often linear flow assumptions (Jeffery Reference Jeffery1922). This is especially the case when adopting the latter beyond its limit of applicability for anisotropic particles of finite-size, i.e. larger than the Kolmogorov length scale.
The modulation of turbulence in non-dilute conditions has been the subject of studies regarding several classes of particulate and multiphase flows, e.g. spherical and isotropic particles (Lucci, Ferrante & Elghobashi Reference Lucci, Ferrante and Elghobashi2010; Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2013; Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2018; Ardekani & Brandt Reference Ardekani and Brandt2019; Ardekani, Rosti & Brandt Reference Ardekani, Rosti and Brandt2019; Sozza et al. Reference Sozza, Cencini, Musacchio and Boffetta2020), anisotropic particles and fibres (Andersson, Zhao & Barri Reference Andersson, Zhao and Barri2012; Ardekani et al. Reference Ardekani, Costa, Breugem, Picano and Brandt2017a; Olivieri et al. Reference Olivieri, Brandt, Rosti and Mazzino2020b, Reference Olivieri, Mazzino and Rosti2021), as well as droplets or bubbles (Dodd & Ferrante Reference Dodd and Ferrante2016; Freund & Ferrante Reference Freund and Ferrante2019; Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019; Cannon et al. Reference Cannon, Izbassarov, Tammisola, Brandt and Rosti2021). Overall, when the concentration is large enough, this modulation effect typically causes a substantial departure from the classical phenomenology observed for a purely Newtonian fluid (e.g. the presence of a clear energy cascade for sufficiently large Reynolds number). Specifically, some common features are generally observed looking at the turbulent energy spectra despite the remarkable differences between these multiphase flows (Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2013; Dodd & Ferrante Reference Dodd and Ferrante2016; Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019; Olivieri et al. Reference Olivieri, Mazzino and Rosti2021): with respect to the reference single-phase case, it is typically found (i) a decrease of the energy content at the lowest wavenumbers (i.e. the largest, energy-containing scales), along with (ii) a relative enhancement of energy at higher wavenumbers (i.e. smaller scales). Nevertheless, both the qualitative and quantitative comprehension on the mechanisms underlying such complex energy redistribution process is far from being exhaustive. Furthermore, these effects are not only directly associated with the resulting flow dynamics, but are in turn potentially influential on the behaviour of the dispersed objects.
In this work, we investigate numerically the behaviour of a suspension of finite-size, inertial fibres, either rigid or flexible, dispersed in a homogeneous and isotropic turbulent flow in both dilute and non-dilute conditions. An illustrative example of the system under consideration is shown in figure 1. A vast parametric study is conducted over the main mechanical and geometrical properties of the fibres, i.e. linear density, length and bending stiffness, as well as their number, while neglecting the presence of gravity. In order to accurately describe the dynamics of fibres, whose length is well beyond the Kolmogorov length scale, and to take into account the backreaction of the dispersed phase to the carrier flow, a four-way coupled, direct numerical simulation (DNS) approach is employed, based on the combination of the finite difference and immersed boundary methods, and complemented by a fibre-to-fibre interaction model. The goal of the work is twofold: (i) to properly characterise the backreaction effect exerted by fibres to the turbulent flow, highlighting the role of the mass fraction as the representative control parameter and providing a scale-by-scale description of the energy transfer across the scales of the fluid motion; (ii) to describe various aspects of the dynamics of the dispersed fibres in the turbulent flow, including the existing flapping states, as well as the preferential concentration (i.e. tendency to clustering) and orientation (i.e. alignment with specific quantities of the flow). Note that this investigation partially recovers the results presented by Olivieri et al. (Reference Olivieri, Mazzino and Rosti2021), with the goal of substantially advancing such findings and offer a systematic and fully exhaustive analysis.
Snapshot from a direct numerical simulation of homogeneous isotropic turbulence where finite-size, flexible fibres (in red) are dispersed. In particular, we show a high-concentration case with very flexible denser-than-the-fluid fibres of intermediate length comparable to the inertial range of scales (i.e. 
$\Delta \tilde {\rho } = 1$, 
$c=0.5$, 
$\gamma =10^{-8}$ and 
$N=10^3$). The domain backplanes are coloured by the fluid velocity magnitude.
The rest of the paper is structured as follows: in § 2, we introduce the computational methodology used in our investigation along with providing information on the performed parametric study; in § 3, we present and discuss the results to address the two main objectives; finally, concluding remarks are given in § 4.
2. Methods
This section describes the methodology adopted for the present study, which is essentially the same adopted by Olivieri et al. (Reference Olivieri, Mazzino and Rosti2021). First, in § 2.1, we introduce the governing equations and main parameters of the problem. Then, in § 2.2, we describe how the former are solved numerically. Finally, in § 2.3, we provide details on the performed parametric DNS study.
2.1. Problem formulation
We consider an ensemble of 
$N$ fibres moving freely in homogeneous and isotropic turbulence (figure 1). The fibre's length 
$c$ is varied with respect to the characteristic length scales of the turbulent flow, ranging from being of the order of the dissipative (or Kolmogorov) scale to that of the integral scale. Fibres are modelled as one-dimensional slender objects, i.e. their length is assumed much larger than the diameter 
$d$, or equivalently we consider an aspect ratio 
$c/d \gg 1$. Moreover, we assume that the fibres are inextensible and focus solely on the elastic contribution due to the bending stiffness 
$\gamma$ (which for a homogeneous fibre is given by the product of the elastic modulus and the second moment of the area). Finally, 
$\Delta \tilde {\rho } = \tilde {\rho }_{s} - \tilde {\rho }_{f}$ is the difference in the linear density (linear and volumetric densities are denoted with and without the tilde, respectively) of the fibre and the fluid, i.e. the parameter controlling the inertia of the dispersed objects. In particular, in the limit 
$\Delta \tilde {\rho } \to 0$ we recover the iso-dense case. Here, we assume zero-gravity conditions in order to decrease the number of independent parameters. In particular, we note that the choice is justified when focusing on the limit of arbitrarily large Froude number.
The dynamics of each fibre evolves according to the Euler–Bernoulli beam equation for an elastic filament coupled with the inextensibility constraint, which read
$$\begin{gather} \Delta \tilde{\rho} \, \ddot{{\boldsymbol{X}}}= \partial_s \left( T \partial_s {\boldsymbol{X}} \right) - \gamma \partial^4_s {\boldsymbol{X}} - \boldsymbol{F}_{fs} + {\boldsymbol{F}}_{col}, \end{gather}$$
$$\begin{gather}\partial_s {\boldsymbol{X}} \boldsymbol{\cdot} \partial_s {\boldsymbol{X}} = 1, \end{gather}$$
where 
$\boldsymbol {X} ( s, t )$ is the position of a generic material point belonging to the fibre, 
$T ( s,t )$ is the tension enforcing the inextensibility constraint in the filament equation, 
$\boldsymbol {F}_{fs} ( s, t )$ is the fluid–solid interaction forcing and 
${\boldsymbol {F}}_{col} ( s, t )$ is a fibre-to-fibre collision forcing term (both forcings will be discussed in the following); all quantities are a function of the curvilinear coordinate 
$s$ and time 
$t$. Freely moving boundary conditions are imposed at both fibre's ends: 
$\partial _{ss} {\boldsymbol {X}} |_{s=0,c} = 0$, 
$\partial _{sss} {\boldsymbol {X}} |_{s=0,c} = 0$ and 
$T |_{s=0,c} = 0$. From a normal mode analysis in the case of such configuration (i.e. free–free or unsupported beam) we can readily obtain the natural frequency 
$f_{nat} = \alpha \, \sqrt {\gamma / (\Delta \tilde {\rho } c^4)}$, where 
$\alpha \approx 22.4/{\rm \pi}$. Dealing with finite-size objects, we neglect the Brownian contribution to the fibre dynamics, which is considered to be small compared with the hydrodynamic forcing.
Fibres are released within a cubic domain of size 
$L=2{\rm \pi}$ having periodic boundary conditions in all directions and where a homogeneous and isotropic turbulent flow is generated using the stochastic spectral forcing by Eswaran & Pope (Reference Eswaran and Pope1988), randomly injecting energy at the largest scales, i.e. within a low-wavenumber shell with 
$1\leqslant k \leqslant 2$ (with 
$k$ denoting the wavenumber). The carrier flow is governed by the incompressible Navier–Stokes equations for a Newtonian fluid, which read
$$\begin{gather} \partial_t {\boldsymbol{u}} + {\boldsymbol{u}} \boldsymbol{\cdot} {\boldsymbol{\partial}} {\boldsymbol{u}} ={-} {\boldsymbol{\partial}} p / {\rho_{f}} + \nu \partial^2 {\boldsymbol{u}} + {\boldsymbol{f}}_{tur} + {\boldsymbol{f}}_{fs}, \end{gather}$$
$$\begin{gather}{\boldsymbol{\partial}} \boldsymbol{\cdot} {\boldsymbol{u}} = 0, \end{gather}$$
where 
${\boldsymbol {u}} ( {\boldsymbol {x}},t )$ and 
$p( {\boldsymbol {x}},t )$ are the velocity and pressure fields, respectively, 
$\rho _{f}$ is the volumetric fluid density and 
$\nu$ is the kinematic viscosity; 
${\boldsymbol {f}}_{tur} ( {\boldsymbol {x}},t )$ is the external forcing used to sustain the turbulent flow (Eswaran & Pope Reference Eswaran and Pope1988) and 
${\boldsymbol {f}}_{fs} ( {\boldsymbol {x}},t )$ is an additional forcing mimicking the presence of the solid phase (as described in the following); lastly, 
${\boldsymbol {x}}$ denotes the position vector in the Eulerian framework.
The flow and fibre dynamics are coupled by the no-slip condition 
$\dot {{\boldsymbol {X}}} = {\boldsymbol {u}} ( {\boldsymbol {X}}( s, t ),t )$.
2.2. Numerical technique
To solve numerically the fully coupled FSI problem, we employ the IBM procedure originally proposed by Huang, Shin & Sung (Reference Huang, Shin and Sung2007) and subsequently modified by Banaei, Rosti & Brandt (Reference Banaei, Rosti and Brandt2020) and Olivieri et al. (Reference Olivieri, Akoush, Brandt, Rosti and Mazzino2020a,Reference Olivieri, Brandt, Rosti and Mazzinob, Reference Olivieri, Mazzino and Rosti2021). In particular, the mutual interaction between the fluid and solid phase is enforced indirectly by means of the singular force distributions 
${\boldsymbol {F}}_{fs}$ and 
${\boldsymbol {f}}_{fs}$ acting on the fibre and flow, respectively (Peskin Reference Peskin2002; Huang et al. Reference Huang, Shin and Sung2007). The fluid velocity at the position of the fibre Lagrangian point, 
$\boldsymbol {U} = {\boldsymbol {u}} ( \boldsymbol {X} ( s,t ),t )$, is obtained by interpolating the fluid velocity at the Eulerian nodes surrounding the Lagrangian point:
\begin{equation} \boldsymbol{U} \left( \boldsymbol{X} \left( s,t \right),t \right) = \int \boldsymbol{u} \left( \boldsymbol{x},t \right) \delta \left( \boldsymbol{x} - \boldsymbol{X}\left( s,t \right) \right) \mathrm{d}^3 \boldsymbol{x}, \end{equation}
where 
$\delta$ is the Dirac delta function. The interpolated velocity 
${\boldsymbol {U}}$ is used to compute the fluid–solid interaction forcing needed to enforce the no-slip condition
\begin{equation} {\boldsymbol{F}}_{fs}(s,t) = \beta \, ({\boldsymbol{U}} - \dot{{\boldsymbol{X}}}), \end{equation}
where 
$\beta$ is a large negative constant (Huang et al. Reference Huang, Shin and Sung2007). Finally, 
${\boldsymbol {F}}$ is spread to the surrounding fluid flow as
\begin{equation} {\boldsymbol{f}}_{fs} \left( \boldsymbol{x},t \right) = \frac{1}{{\rho}_{f}} \int_s {\boldsymbol{F}}_{fs} \left( s,t \right) \delta \left( \boldsymbol{x} - \boldsymbol{X}\left( s,t \right) \right) \mathrm{d} s. \end{equation}
Both the interpolation and spreading operations feature the Dirac operator. Numerically, this is transposed into the use of the regularised 
$\delta$ function proposed by Roma, Peskin & Berger (Reference Roma, Peskin and Berger1999). Consequently, in the numerical framework the fibre's diameter is also approximated to a finite value comparable to the grid spacing (note that the present method is not able to properly resolve the boundary layer around the individual fibre diameter); in the following we assume 
$d = 2 \, \Delta x$ as the effective diameter.
Each fibre is evenly discretised using 
$N_{L}$ Lagrangian points such that the spatial resolution 
$\Delta s = c / (N_{L} -1)$ is approximately equal to the Eulerian grid spacing 
$\Delta x$; consequently, the discretised curvilinear coordinate can be evaluated as 
$s_l = l \, \Delta s$, with 
$l = 1, \dots, N_{L}$. The numerical solution of (2.1) and (2.2) follows the scheme detailed in Huang et al. (Reference Huang, Shin and Sung2007) with the difference that the bending term is treated implicitly to allow for a larger timestep (Banaei et al. Reference Banaei, Rosti and Brandt2020; Olivieri et al. Reference Olivieri, Mazzino and Rosti2021). Furthermore, a fibre-to-fibre collision model is implemented to prevent that the distance 
$\boldsymbol {d}$ between Lagrangian points belonging to different fibres goes below the grid spacing 
$\Delta x$ and that fibres eventually cross each other. Specifically, we employ the minimal collision model proposed by Snook, Guazzelli & Butler (Reference Snook, Guazzelli and Butler2012), where a constant forcing 
${\boldsymbol {F}}_{col} = F_0 \, \hat {\boldsymbol {d}}$ is applied when 
$|\boldsymbol {d}| \leqslant \varDelta _{col}$. Here, 
$F_0$ is a free parameter and 
$\varDelta _{col} = O(\Delta x)$ is the critical distance below which the forcing is imposed. After several tests on the sensitivity of the results with respect to these parameters (both in simple laminar flow conditions with few fibres and in the turbulent case with a non-dilute suspension), we chose 
$F_0 = 1.0$ and 
$\varDelta _{col} = 3 \Delta x$. As a result, the influence of collisions turns out to be very limited, as shown in Appendix A: fibre-to-fibre interactions are extremely rare in the majority of the analyzed configurations and are always found to have a negligible effect on the numerical results, especially on the turbulence modulation and fluid–solid coupling mechanisms which are the main focus of this work.
Concerning the fluid flow, (2.3) and (2.4) are solved using the fractional step method on a staggered grid (Kim & Peskin Reference Kim and Peskin2007). The Poisson equation enforcing the incompressibility constraint is solved using a fast and efficient approach based on the fast Fourier transform (FFT). The numerical solution is based on the (second-order) central finite-difference method for the spatial discretisation and the (third-order) Runge–Kutta scheme for the temporal discretisation.
The described computational procedure is implemented in the in-house solver Fujin (https://groups.oist.jp/cffu/code). The code is parallelised using the MPI protocol and the 2decomp library for domain decomposition (http://www.2decomp.org). It has been extensively tested and employed in a variety of fluid dynamical problems (Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019; Rosti & Brandt Reference Rosti and Brandt2020; Rosti et al. Reference Rosti, Olivieri, Cavaiola, Seminara and Mazzino2020b, Reference Rosti, Cavaiola, Olivieri, Seminara and Mazzino2021), including in particular suspensions of fibres in both laminar and turbulent flows (Rosti et al. Reference Rosti, Banaei, Brandt and Mazzino2018, Reference Rosti, Olivieri, Banaei, Brandt and Mazzino2020a; Cavaiola, Olivieri & Mazzino Reference Cavaiola, Olivieri and Mazzino2020; Olivieri et al. Reference Olivieri, Akoush, Brandt, Rosti and Mazzino2020a,Reference Olivieri, Brandt, Rosti and Mazzinob, Reference Olivieri, Mazzino and Rosti2021; Brizzolara et al. Reference Brizzolara, Rosti, Olivieri, Brandt, Holzner and Mazzino2021). To further enrich the existing framework, in Appendix B we show a comparison of our results on the self-sustained flapping oscillation of a two-dimensional, hinged filament in a laminar flow with those originally reported by Huang et al. (Reference Huang, Shin and Sung2007).
2.3. Parametric study
The objective of this work is to perform an extensive parametric study over the characteristic properties of the fibre suspension. Therefore, we have performed DNS varying the fibre's linear density difference 
$\Delta \tilde {\rho }$, length 
$c$ and bending stiffness 
$\gamma$, along with the number of fibres 
$N$. Overall, this also leads to consider different concentrations of the suspension, i.e. from very dilute to non-dilute suspensions. Specifically, in a first baseline investigation we considered two different linear densities 
$\Delta \tilde {\rho }$, corresponding to essentially iso-dense and denser-than-the-fluid fibres, three different values of the fibre's length 
$c$, which are comparable to the dissipative, inertial and integral characteristic lengthscales of the turbulent flow, three different bending stiffnesses 
$\gamma$ so that we range from highly flexible to essentially rigid fibres. Finally, we varied the numbers of fibres 
$N$ to adjust the concentration of the suspension. Tables 1 to 3 in Appendix C report the parametric combinations along with the consequent relevant metrics for the concentration, e.g. number density or mass fraction, and the corresponding symbol used in the rest of the paper to refer to each case.
A reference configuration is chosen in the single-phase case (i.e. for 
$N=0$), where the turbulent flow is characterised by a micro-scale Reynolds number 
${Re}_\lambda \equiv u_{rms} \lambda / \nu \approx 120$, where 
$u_ {rms}$ is the root-mean-square velocity and 
$\lambda$ is the Taylor micro-scale. To simulate such configuration, the fluid domain is discretised into a uniform Eulerian grid using 
$256^3$ cells, having verified that the ratio between the Kolmogorov dissipative length scale and the grid spacing is 
$\eta /\Delta x = O(1)$. Moreover, we verified the convergence of the numerical results by doubling the spatial resolution in one representative case (not shown), observing negligible differences with respect to the adopted grid setting. After reaching the statistically stationary state in the single-phase configuration, the obtained flowfield is used as the initial condition for the multiphase flow simulations, i.e. fibres are added to the fully-developed turbulent carrier flow. The simulation is then evolved to reach the new stationary regime, the achievement of the latter being qualitatively assessed from the time histories of main statistical quantities for both phases (e.g. average kinetic energy). Hence, we disregard the transient and accumulate the statistics over a minimum of five integral timescales for all cases. The statistical convergence of our data concerning both the flow and fibre dynamics was verified by widely varying the number of samples and checking their sensitivity on the results, ensuring that all the reported differences are substantially larger than the statistical uncertainty.
3. Results
This section presents the results of our study. First, in § 3.1, we focus on the characterisation of the modulation of the turbulent flow due to the dispersed phase. Then, in § 3.2, we analyse in more detail several features of the fibre motion and transport.
3.1. Flow dynamics
3.1.1. Macroscopic behaviour
We begin our analysis by looking at the bulk properties of the turbulent flow when varying the parameters of the suspension. To start, figure 2(a–c) shows the dependence of the micro-scale Reynolds number 
${Re}_\lambda$ as a function of several quantities that can be used as an indicator to characterise the concentration of the suspension: (i) the (non-dimensionalised) number density 
$n \, c^3$, (ii) the volume fraction 
$\phi _V$ and (iii) the mass fraction 
$\mathcal {M}$. These three quantities are chosen because they are often used in the framework of fibre suspensions, or more generally particle-laden flows, to estimate (a priori) the importance of the backreaction of the dispersed phase on the carrier flow (Butler & Snook Reference Butler and Snook2018; Brandt & Coletti Reference Brandt and Coletti2021). From figure 2, it can be noted that for several cases the Reynolds number turns out to decrease with respect to the single-phase case (the value of the latter indicated by the horizontal dashed line in the figure); this appears as a robust effect that will be characterised in detail in the rest of the section. However, the ways in which the data appear when plotted as a function of the three different parameters are radically different.
Micro-scale Reynolds number as a function of the (a) (non-dimensional) number density, (b) volume fraction and (c) mass fraction (each of these indicator quantities are defined in the text; the dashed line in each plot indicates the reference value obtained in the single-phase case) and (d) turbulent kinetic energy 
$\mathcal{E}$ (normalised with the value of the single-phase case 
$\mathcal{E}_0$) as a function of the mass fraction. Colours and symbols are used consistently with the indication of tables 1 and 2: empty and filled markers denote the cases for iso-dense (
$\Delta \tilde {\rho } = 10^{-3}$) and denser-than-the-fluid (
$\Delta \tilde {\rho } = 10^0$) fibres, respectively. The colour denotes the fibre's length (green, 
$c=0.1$, i.e. short; red, 
$c=0.5$, i.e. intermediate; blue, 
$c=2.0$, i.e. long). The symbol indicates the bending stiffness (squares, 
$\gamma =10^{-8}$; triangles, 
$\gamma =10^{-4}$; rotated squares, 
$\gamma =10^{0}$). Finally, the brightness decreases with the number of fibres (light, 
$N=10^1$; medium, 
$N=10^2$; dark, 
$N=10^3$).
The (dimensional) number density, defined as 
$n = N / V_{f}$ (where 
$V_{f} = L^3$ is the volume of the domain), is a quantity widely used to investigate fibre suspensions, especially in the framework of rheological studies in low-Reynolds-number flows. Such quantity is typically made non-dimensional with the fibre length 
$c$, so that 
$n c^3$ represents a measure of the relative spacing between fibres and consequently of the degree of their mutual interaction. More specifically, the suspension is typically considered to be in the dilute regime if 
$n c^3 \ll 1$, whereas it is classified as semi-dilute for 
$1/c^{3} \leqslant n \ll 1/ (d \, c^2)$. If 
$n \gg 1/ (d \, c^2)$, i.e. the spacing between the fibres is below their diameter, we finally have the concentrated regime, where the fibre-to-fibre interactions are expected to be dominant in controlling the dynamics of the dispersed objects (Butler & Snook Reference Butler and Snook2018). We note that all the configurations considered in the present study fall either in the dilute or in the semi-dilute regime. As shown in figure 2(a), however, the number density does not appear to be a good indicator for describing the entity of the backreaction on the turbulent flow, without any systematic trend observed when plotting the numerical results with respect to such parameter. As will become clearer in the following, the reason for this is the incorrect scaling with the fibre length.
Let us now consider the same data but as a function of the volume fraction 
$\phi _V$, as reported in figure 2(b). Here, 
$\phi _V = V_{s} / V_{f}$ is defined as the ratio between the volume of the dispersed phase 
$V_{s} = N \, c \, {\rm \pi}d^2 / 4$ and that of the fluid domain 
$V_{f} = L^3$. Using this quantity, the results are outlined with a much more defined trend (compared with those for the number density), with 
${Re}_\lambda$ decreasing as 
$\phi _V$ is increased. However, it can be noted that such indicator misses to take into account the inertia of the fibres, here quantified by the linear density difference 
$\Delta \tilde {\rho }$. Indeed, the latter clearly appears as an additional key (and free) parameter, with a dramatic difference between the results for the iso-dense (empty symbols) and the denser-than-the-fluid (filled symbols) fibre cases.
In light of the evidence, it is therefore natural to consider, as another way to parametrise the backreaction effect, the mass fraction 
$\mathcal {M} = m_{s} / (m_{s} +m_{f})$, where 
$m_{s} = \rho _{s} V_{s}$ is the (total) mass of the solid dispersed phase and 
$m_{f} = \rho _{f} \, V_{f}$ is the (total) mass of fluid in the domain. As a result, from figure 2(c) it can be now appreciated how the mass fraction turns out to be the most representative parameter for describing the variation of the backreaction and consequent turbulence modulation at a macroscopic level. We highlight that, although this fact is well known for the case of point-like inertial particles (Brandt & Coletti Reference Brandt and Coletti2021), similar evidence is not reported for suspensions of finite-size fibres, such as those investigated in the present work. Focusing on figure 2(c), it can be observed that the Reynolds number starts to vary appreciably from the single phase value when the mass fraction becomes larger than around 
$10\,\%$, thus systematically decreasing with 
$\mathcal {M}$ and eventually reaching, for the highest mass fraction that was tested (
$\mathcal {M}\approx 89\,\%$), a value that is roughly half of that obtained in the single-phase case (
${Re}_\lambda \approx 60$).
For a more complete characterisation of the turbulence modulation, in figure 2(d) we also show how the turbulent kinetic energy 
$\mathcal{E}$ decreases with the mass fraction, following a trend that is qualitatively similar to what observed for the micro-scale Reynolds number.
From figure 2, a further comment can be made on a secondary yet systematic effect associated with the variation of the fibre's bending stiffness, which is especially evident for the longest fibres at the highest concentration (dark blue symbols). Remarkably, the Reynolds number is systematically found to decrease when the bending stiffness 
$\gamma$ is increased (the fibre is more rigid). This feature is particularly evident for the cases with the longest fibres and can be intuitively associated with the degree of compliance manifested by the elastic objects under the action of the flow, which reflects in a stronger friction exerted on the flow. Nevertheless, it can be also pointed out the huge variation (up to eight orders of magnitude) in the bending stiffness that was considered for the reported cases, in light of which we conclude that the effect associated with such parameter is overall largely subleading when compared with that caused by the inertia of the dispersed phase.
3.1.2. Energy spectra
To better understand the multiscale features of the backreaction of the dispersed phase on the carrier flow, in figure 3 we show the energy spectra computed from our simulations for the two linear densities 
$\Delta \tilde {\rho }$, i.e. iso-dense (a,c,e) versus denser-than-the-fluid fibres (b,d, f), and the three fibre lengths 
$c$ (corresponding to each row and increasing from top to bottom) investigated. In each plot, we vary the number of fibres 
$N$ (increasing from light to dark colour). Moreover, for the sake of comparison, the results from the single-phase case are also reported (black curve). Here, we retain the same intermediate value for the bending stiffness 
$\gamma$ because its influence on the turbulence modulation was previously shown to be subleading when compared with the combined variation of the other parameters determining the mass fraction of the suspension.
Energy spectra of the modulated turbulent flow for different fibre suspension: (a,c,e) iso-dense fibres; (b,d, f) denser-than-the-fluid fibres; (a,b) short fibres; (c,d) intermediate fibres; (e, f) long fibres. The number of fibres, and therefore the mass fraction, increases from light to dark. Colours and symbols are used consistently with tables 1 and 2. All cases are relative to the intermediate value of the bending stiffness. The black lines (without symbols) refer to the single-phase case.
Figure 3 shows that, as the concentration is increased, the energy spectrum is modified with the same qualitative behaviour in all the cases. The latter can be essentially described as the combination of an energy content depletion at the largest scales (i.e. low wavenumbers), along with a (relative) enhancement at the smallest scales (i.e. high wavenumbers), with respect to the single-phase configuration. Moreover, the energy distribution in the intermediate range of scales is also appreciably modified, giving rise to a rather different phenomenology compared with that of classical turbulence. Clearly, the effect is more pronounced when increasing the fibre density and/or length, consistently with the increase in the mass fraction as outlined from the macroscopic characterisation of the bulk flow properties.
To highlight that qualitatively similar flow conditions can be obtained using significantly different combinations of the suspension parameters, we can compare, e.g. the case of short, denser-than-the-fluid fibres in figure 3(b) with that of long, iso-dense fibres in figure 3(e) (considering for both cases the highest number of fibres). Although some quantitative differences are present (with the former showing a slightly stronger modulation), when compared with the single-phase case the resemblance between the two configurations is evident. Another observation can be inferred by looking at the enhanced energy content at the smallest cases, a robust feature manifesting already at relatively low mass fraction, and which can have potential relevance for mixing processes. On the other hand, it should be noted that the overall kinetic energy of the flow monotonically decreases with the mass fraction, consistently with the large-scale energy depletion previously discussed.
Lastly, we explore the effect of the bending stiffness on the energy spectrum for one representative configuration (i.e., denser-than-the-fluid fibres of intermediate length and maximum concentration), shown in figure 4. On one hand, it can be noted that the energy content at the large scales, i.e. low wavenumbers, is systematically decreasing as the bending stiffness is increasing, in agreement with what observed in terms of the overall energy depletion (cf. figure 2). On the other hand, in the high-wavenumber region the various curves are almost superimposed. From the physical viewpoint, we can argue that increasing the flexibility (i.e. decreasing 
$\gamma$) leads to a more relaxed flow–structure coupling (i.e. fibres are more compliant to the action of the flow), reflected in a less-pronounced backreaction if compared with the rigid case. Note that the saturation observed in figure 4 when increasing 
$\gamma$ corresponds indeed to the convergence towards the rigid configuration. In this regard, a more quantitative indication on the expected role of elasticity can be obtained in terms of the ratio between the structural and fluid characteristic timescales (later introduced in § 3.2), from which it can be deduced the extensive (non-dimensional) range that has been considered (cf. figure 7). Nonetheless, it has also to be pointed out that the variation with 
$\gamma$ appears overall very limited, confirming the secondary effect of flexibility in altering the carrier flow.
Effect of fibre's bending stiffness on the turbulent energy spectrum. Results are obtained for a suspension of 
$N=10^3$ denser-than-the-fluid fibres of intermediate length, varying the fibre's bending stiffness over eight orders of magnitude. The black curve refers to the single-phase case. Symbols are used accordingly to table 3.
3.1.3. Scale-by-scale energy transfer
In order to characterise the mechanisms underlying the outlined phenomenology, we can analyse the scale-by-scale energy transfer that for the present multiphase flow problem can be written as
\begin{equation} P(k) + \varPi(k) + \varPi_{fs}(k) + D(k) = \epsilon, \end{equation}
where the various terms appearing on the left-hand side are the production rate 
$P$ (here associated with the external forcing used to sustain the flow and acting only at the largest scales, i.e. low wavenumbers), the energy flux 
$\varPi$ associated with the nonlinear term, the energy flux 
$\varPi _{fs}$ associated with the fluid–solid coupling, and the viscous dissipation 
$D$. Note that the latter is defined such that the turbulent energy dissipation rate 
$\epsilon$ appears on the right-hand side and the overall balance can be visualised more conveniently. For a derivation of (3.1), see Appendix D. In the single-phase case 
$\varPi _{fs} \equiv 0$, and the dominant terms in the balance are the nonlinear term 
$\varPi$ and the viscous dissipation 
$D$ when approaching the lowest and highest wavenumbers, respectively (Pope Reference Pope2000).
In this analysis, we focus on the cases with 
$N=10^3$. However, the generality of the results is not restricted because the mass fraction, and consequently the entity of the backreaction, still varies considerably, thus ranging from configurations with negligible backreaction (i.e. essentially resembling the case without fibres) to configurations where the backreaction is strong and leads to a dramatic departure from the more classical scenario obtained for purely Newtonian fluid turbulence. Similarly, for the bending stiffness we select at first the same (intermediate) value as for the energy spectra reported in figure 3, focusing on the effect of the mass fraction.
Let us therefore start from considering the spectral power balance for the cases of iso-dense fibres (with different fibre lengths), shown in figure 5(a,c,e). Here, the mass fraction is always relatively small; in agreement with the previously discussed effects on the Reynolds number and energy spectrum, we therefore expect an overall negligible, or at most quite limited, alteration in the balance with respect to the single-phase configuration. Indeed, starting from the short fibres (figure 5a) it can be clearly observed that the results resemble those of the single-phase configuration (indicated by the black lines without symbols), with a predominance at low wavenumbers of the nonlinear term 
$\varPi$ and a negligible energy transfer associated with the fluid-structure coupling 
$\varPi _{fs}$. This situation is accompanied by the predominance at high wavenumbers of the dissipation term 
$D$ which eventually saturates to 
$\epsilon$ for increasing 
$k$. The same applies also to fibres of intermediate length (figure 5c), with only a very modest increase of 
$\varPi _{fs}$. On the other hand, when considering the case of long fibres (figure 5e) the situation starts to be different: at low wavenumbers, the nonlinear term is weakened and does not entirely control anymore the energy transfer, with a non-negligible contribution of the fluid–solid coupling which is maximised at an intermediate length scale; as 
$k$ increases, the dissipation progressively becomes dominant, although the saturation is less evident due to the enhanced small-scale energy transfer. Remarkably, the mass fraction 
$\mathcal {M}$ in this case is only around 
$1\,\%$, therefore indicating that a substantially different and complex dynamics of the turbulent flow may arise already in relatively dilute conditions.
Spectral power balance according to (3.1), showing the energy fluxes due to the nonlinear term (solid lines), fluid–solid coupling (dashed) and viscous dissipation (dotted), all normalised with the turbulent energy dissipation rate, for the cases with 
$N=10^3$ and intermediate bending stiffness: (a,c,e) iso-dense fibres; (b,d, f) denser-than-the-fluid fibres; (a,b) short fibres; (c,d) intermediate fibres; (e, f) long fibres. In addition, in (a) the nonlinear and dissipation terms of the single-phase case are also reported for comparison. Colours and symbols are used consistently with the indication of tables 1 and 2.
The alteration gets substantially more dramatic when considering the same three cases but for denser-than-the-fluid fibres, reported in figure 5(b,d, f). Already for short fibres (figure 5b) a similar behaviour can be noted, with a non-negligible fraction of the energy transfer that is associated to the fluid-solid coupling term. In fact, this case presented similarities already in terms of the energy spectrum with that with iso-dense long fibres, as mentioned previously. When increasing the mass fraction and considering the intermediate fibre length (figure 5d), the contribution of the nonlinear and fluid–solid coupling terms become comparable in magnitude. However, the former is dominant at lower wavenumbers, indicating that the energy from the largest scales is at first mainly transferred to the smaller scales by means of the classical hydrodynamic mechanism, although in a restricted subrange of scales compared with the single-phase case. Then, at sufficiently higher wavenumbers, such energy transfer is progressively replaced by another mechanism driven by the direct interaction between the carrier flow and the dispersed phase. We stress the fact that, as the mass fraction is increased, the nonlinear term decreases while the fluid–solid coupling is increased, whereas the viscous dissipation always eventually prevails at sufficiently large wavenumbers, i.e. small scales. Eventually, for the case of long denser-than-the-fluid fibres having the largest mass fraction (figure 5f), the nonlinear term is essentially suppressed with the fluid–solid coupling largely controlling the overall energy transfer, up to the dissipative scales.
As an overall observation, for all cases the energy transfer is eventually compensated by the viscous term, since both the nonlinear and fluid–solid coupling term vanish when computed over the full range of wavenumbers. This aspect concerns a remarkable difference between the case of elastic fibres, as those considered in the present work, and that of polymer suspensions and viscoelastic turbulence, where the polymer stresses both transfer and dissipate energy (De Angelis et al. Reference De Angelis, Casciola, Benzi and Piva2005; Valente, Da Silva & Pinho Reference Valente, Da Silva and Pinho2014). In our model, no internal dissipation mechanism is considered for the dispersed phase and therefore the fluid–solid coupling purely transfers energy from larger to smaller scales, in a strict analogy with the nonlinear term, with the total energy flux given by the sum of the two contributions.
To conclude the analysis, there remains to discuss in more detail the influence of the fibre's flexibility, quantified by the bending stiffness 
$\gamma$, on the backreaction and turbulence modulation. Figure 6 reports the energy flux from the nonlinear (a,c,e) and fluid–solid coupling (b,d, f) contributions for the cases with maximum and minimum bending stiffness (along with different density and length) in order to highlight the overall range of variation associated with this parameter. The plots confirm once more that the variation is always quite modest, if compared with that given by the other parameters (i.e. more universally expressed in terms of the mass fraction). Nevertheless, a systematic trend in all cases can be detected: as the bending stiffness is increased, the magnitude of the nonlinear term is found to decrease whereas the fluid–solid coupling contribution becomes more relevant. This effect, which could be intuitively expected, can be associated with the fact that rigid fibres are less compliant and do not easily adapt to the stresses applied by the turbulent flow in the same way as very flexible and more compliant fibres do. Remarkably, this mechanism is intrinsically different from the one governed by the inertia of the dispersed phase, and its potential relevance could be better outlined in very dense suspensions of iso-dense fibres.
Energy fluxes associated with the nonlinear term (a,c,e) and fluid-solid coupling term (b,d, f) for 
$N=10^3$ fibres with minimum versus maximum bending stiffness, along with different linear density and length. Colours and symbols are used consistently with the indication of tables 1 and 2: empty and filled markers refer to iso-dense and denser-than-the-fluid fibres, respectively. The symbol indicates the bending stiffness (squares, minimum; diamonds, maximum). The colour denotes the fibre's length (green, short; red, intermediate; blue, long).
3.2. Fibre dynamics
3.2.1. Flapping frequency
We now turn our attention into the dynamical behaviour of the fibres when freely dispersed in the turbulent flow, focusing at first on the main features of their individual motion and deformation. This topic has been investigated by Rosti et al. (Reference Rosti, Banaei, Brandt and Mazzino2018, Reference Rosti, Olivieri, Banaei, Brandt and Mazzino2020a) in the context of very dilute suspensions, and more recently by Olivieri et al. (Reference Olivieri, Mazzino and Rosti2021) with the aim of extending the analysis to non-dilute configurations. These previous studies highlighted the potential of using finite-size flexible fibres to measure relevant two-point statistics of turbulence, i.e. longitudinal velocity differences and structure functions, under a proper choice of the fibre's mechanical properties. More specifically, based on the comparison between the characteristic timescales that can be identified in the problem, they showed the existence of two main dynamical regimes, i.e. underdamped or overdamped, differing by the fact that the structural elasticity may manifest or not in the resulting dynamics. Furthermore, the classification could further be divided into different dynamical subregimes. However, only two qualitatively different flapping states were found to ultimately take place (Rosti et al. Reference Rosti, Olivieri, Banaei, Brandt and Mazzino2020a; Olivieri et al. Reference Olivieri, Mazzino and Rosti2021): (i) for sufficiently flexible and/or iso-dense fibres (i.e. 
$f_{nat}/f_{tur} \ll 1$) the dynamics is fully controlled by the flow, and therefore the fibres act effectively as a proxy of the turbulent eddies of comparable size; (ii) for sufficiently rigid and inertial fibres (i.e. 
$f_{nat}/f_{tur} \gg 1$) the characteristic response time is related to the natural frequency, and the fibre is not dynamically controlled by the turbulent forcing.
In this work, we complement our recent results that generalise such scenario in the case of non-negligible backreaction (Olivieri et al. Reference Olivieri, Mazzino and Rosti2021). As a starting point, figure 7 shows the behaviour of the flapping frequency 
$f_{flap}$ (identified as the dominant peak frequency from the time history of the fibre's end-to-end distance, as detailed in Rosti et al. Reference Rosti, Olivieri, Banaei, Brandt and Mazzino2020a) while varying the ratio between the natural and turbulent eddy frequency. As anticipated, for the denser-than-the-fluid fibres (corresponding to the underdamped regime) the latter indeed provides a very good indication to distinguish between the two flapping states, which differ from each other for the scaling law of the flapping frequency: for 
$f_{nat}/f_{tur} \ll 1$, the flapping frequency is controlled by, and therefore scales with, the turbulent eddy frequency (i.e. 
$f_{flap} = f_{tur}$), whereas for 
$f_{nat}/f_{tur} \gg 1$ the flapping frequency is given by the natural one (i.e. 
$f_{flap} = f_{nat}$). On the other hand, when the fibres are iso-dense (and thus correspond to the overdamped regime) they are always controlled by turbulence (i.e. 
$f_{flap} = f_{tur}$). We remark that the present results concern both dilute and non-dilute cases, with the entity of the backreaction depending on the corresponding mass fraction, as discussed in § 3.1. The effective qualitative and quantitative characteristics of the modulated flow can thus be appreciably different from the single-phase configuration. To account for this crucial effect, the turbulent eddy frequency is here evaluated as 
$f_{tur} = \alpha ' \, \sqrt {S_2} / c$, where 
$S_2 = \langle (\delta u_\parallel )^2 \rangle$ is the second-order structure function of the longitudinal velocity difference 
$\delta u_\parallel$ evaluated at the fibre's length scale and 
$\alpha ' = 3.0$ is an overall constant 
$O(1)$. As shown in figure 7, this choice leads to results in overall good agreement with the theoretical prediction even in the non-dilute configurations. Note however, that similar results are found when computing 
$S_2$ by means of Lagrangian fibre tracking (instead of using directly the fluid flow measured on the Eulerian grid), or if using the average end-to-end distance in place of the fibre length 
$c$ (with minimal variations that do not appreciably affect the observed scalings). If, on the other hand, one estimated the hydrodynamic frequency using the dimensional arguments based on Kolmogorov scaling in the inertial subrange, a more substantial departure from the reported findings would be obtained for the non-dilute cases where the backreaction is non-negligible (Olivieri et al. Reference Olivieri, Mazzino and Rosti2021). On the other hand, it can be pointed out that the alteration of 
$f_{tur}$ does not have a primary role in the resulting scaling laws (which remain the same already found for the dilute case).
Fibre's dominant flapping frequency as a function of the (free-response) natural frequency, both normalised with the effective eddy turnover frequency (evaluated at the fibre length scale). Dotted and dashed lines denote the expected scalings in the two flapping states. The inset reports the average elastic energy of the fibres for 
$N=10^3$ denser-than-the-fluid fibres of intermediate length and different bending stiffness. Colours and symbols are used consistently with the indication of tables 1 to 3.
A final comment can be made in this regard when considering the (average) elastic energy stored by the fibres during their deformation, shown in the inset of figure 7. Rosti et al. (Reference Rosti, Banaei, Brandt and Mazzino2018) showed that (for a single fibre and in the absence of any appreciable backreaction) this quantity exhibits a maximum when the natural frequency is equal to the turbulence frequency, because of a resonance condition between these two timescales. To extend this finding, we choose a representative non-dilute configuration and vary only the bending stiffness in order to substantially retain the same backreaction to the flow (cf. table 3). From the inset of the figure, it clearly appears that, also in the presence of strong turbulence modulation, the average elastic energy 
$\langle \mathcal {E}_{el} \rangle = ({1}/{N})\sum _{i=1}^N \int _0^c \frac {1}{2} \gamma \kappa _i^2 (s) \,{\rm d}s$, where 
$\kappa _i$ is the local bending curvature of the 
$i$th fibre, peaks when 
$f_{nat}/f_{tur} \approx 1$. Remarkably, the emergence of such peak confirms the idea of a resonance condition in the structural response occurring also in the non-dilute configuration, hence generalising the theoretical framework previously proposed in the case of negligible backreaction.
3.2.2. Maximum curvature
The maximum curvature experienced by the fibres under the action of the flow is another relevant observable when characterising the fibre's deformation, and we report it in figure 8. Remarkably, this quantity is of paramount importance in fragmentation processes involving, e.g. fibrous microplastics in oceanic turbulence (Allende et al. Reference Allende, Henry and Bec2018; Brouzet et al. Reference Brouzet, Guiné, Dalbe, Favier, Vandenberghe, Villermaux and Verhille2021). Therefore, we focus on understanding how 
$\kappa _{max}$ behaves as a function of the relevant mechanical parameters. We compute the maximum fibre's curvature 
$\kappa _{max}$ as follows: we first evaluate the maximum of 
$\kappa _i (s)$ for each fibre (at each time instant), and then average this quantity over the different fibres (and over different time instants).While the curvature is expected to be inversely proportional to the fibre's bending stiffness 
$\gamma$, less trivial is the prediction and corroboration of scaling laws (if any) to qualitatively describe the dependence from the bending stiffness as well as other parameters.
Maximum fibre curvature (normalised with the fibre length) as a function of the normalised bending stiffness, for (a) iso-dense and (b) denser-than-the-fluid fibres. The dashed and dotted lines indicate the scaling laws obtained when balancing the forcing term with the inextensibility bending contribution or the linear bending term, respectively. Colours and symbols are used consistently with the indication of tables 1–3.
To proceed, we adopt an approach similar to that proposed by Gay et al. (Reference Gay, Favier and Verhille2018) in the limit case of iso-dense fibres and dilute suspension, based on balancing different terms in the fibre's dynamical equation (2.1) estimated by means of dimensional analysis. However, differently from Gay et al. (Reference Gay, Favier and Verhille2018), here we directly compare these terms without introducing an effective elastic length nor focusing on an energy balance, and consequently obtain different scaling laws for the maximum curvature. For the sake of simplicity, the same assumptions, i.e. one-way coupling and inertial subrange scaling, are retained (and will be commented on later). For the forcing term we can write 
${\boldsymbol {F}} \sim \mu \, (\epsilon c)^{1/3}$, where 
$\mu$ is the dynamic viscosity, while the (linear) bending term can be dimensionally expressed as 
$\gamma \partial ^4_s {\boldsymbol {X}} \sim \gamma \kappa _{max} /c^2$. Moreover, as shown by Marheineke & Wegener (Reference Marheineke and Wegener2006) and Gay et al. (Reference Gay, Favier and Verhille2018), from the tension term 
$\partial _s ( T \partial _s {\boldsymbol {X}} )$ one can isolate the contribution 
$\partial _s ( \gamma \kappa ^2 \partial _s {\boldsymbol {X}} )$ ensuring the inextensibility constraint under bending deformations, which can be further split in two terms, one parallel to the fibre's axis and one parallel to the curvature vector; focusing solely on the latter (cubic) contribution, we have 
$\partial _s ( T \partial _s {\boldsymbol {X}} ) \sim \gamma \kappa _{max}^3$. Remarkably, the cubic term is subleading with respect to the linear one for 
$\kappa _{max} c \ll 1$, whereas the opposite occurs for 
$\kappa _{max} c \gg 1$. Therefore, we can predict the following two regimes: (i) for sufficiently rigid fibres, the maximum curvature is given by the balance between the forcing and the linear bending term, thus obtaining that 
$\kappa _{max} c \sim [ \gamma / (\mu \epsilon ^{1/3} c^{10/3})]^{-1}$; (ii) for sufficiently flexible fibres, the balance instead is between the forcing and the cubic bending contribution to the tension term, yielding 
$\kappa _{max} c \sim [ \gamma / (\mu \epsilon ^{1/3} c^{10/3})]^{-1/3}$.
Figure 8 shows the results from our numerical simulations along with the expected scaling laws, both for the iso-dense (a) and denser-than-the-fluid fibre cases (b). Overall, we observe a good agreement of the numerical results with respect to both proposed scalings, with only a limited deviation for the denser-than-the-fluid fibres approximately in the transition region between the two regimes. Nevertheless, some remarks on the underlying assumptions are needed to properly outline the limitations of the proposed model. First, using the inertial subrange scaling for the turbulence forcing is strictly justified only for sufficiently dilute suspensions and fibres of intermediate length. Second, although we note from the numerical results a systematic increase of 
$\kappa _{max}$ with the fibre's inertia (i.e. linear density difference), the influence of the latter is not accounted for in the derivation of the scaling laws. Consequently, it was not possible to obtain a unique master curve collecting all the data together. We note the difference with the scalings that were obtained for the flapping frequency, where the inertial term was indeed considered and had a crucial role to distinguish between the overdamped and underdamped cases. Here instead the same qualitative behaviour is obtained for both the iso-dense and denser-than-the-fluid fibres. Nevertheless, our predictions for the remaining set of parameters are corroborated not only for iso-dense, but also for the denser-than-the-fluid inertial fibres, thus being general for the widest parametric range.
3.2.3. Clustering
After the characterisation of the individual motion of the dispersed fibres, we consider some relevant aspects concerning their collective dynamics. We look in particular at the local concentration of the suspension to inspect the presence of clustering phenomena, a key feature in the analysis of turbulent particle-laden flows (Eaton & Fessler Reference Eaton and Fessler1994; Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007) and with a specific relevance in the context of fibres aggregation (Lundell et al. Reference Lundell, Söderberg and Alfredsson2011; Verhille et al. Reference Verhille, Moulinet, Vandenberghe, Adda-Bedia and Le Gal2017). Different metrics have been proposed to quantify the tendency of particles to preferentially accumulate along with sampling specific regions of the flow (Brandt & Coletti Reference Brandt and Coletti2021), among which is the radial distribution (or pair correlation) function (RDF) (Salazar et al. Reference Salazar, De Jong, Cao, Woodward, Meng and Collins2008; Saw et al. Reference Saw, Shaw, Ayyalasomayajula, Chuang and Gylfason2008; Olivieri et al. Reference Olivieri, Picano, Sardina, Iudicone and Brandt2014). Such a quantity indicates the probability of having a pair of particles at a given mutual distance, thus highlighting the presence of accumulation at particular lengthscales, while assuming a constant unit value when the particle distribution is locally uniform. Here, we compute the RDF using its classical definition for point-like particles and specifically considering the Lagrangian midpoints, being concerned with the mutual distance between the different anisotropic particles. Furthermore, we restrict the analysis to the cases with the largest number of fibres 
$N=10^3$ to ensure well-converged statistics.
Figure 9 shows the RDF obtained from the DNS results, from which several observations can be made. The left panels (figure 9a,c,e) refer to the iso-dense fibres of different length and bending stiffness, for which the clustering effect is always essentially negligible. Indeed, such evidence is consistent with the crucial role of inertia in the formation of clustering (Eaton & Fessler Reference Eaton and Fessler1994). On the other hand, the cases for denser-than-the-fluid fibres in the right panels (figure 9b,d, f) reveal a richer phenomenology. First, for all the investigated fibre lengths, it turns out that the accumulation increases when decreasing the bending stiffness. This systematic trend can be explained by the fact that, when the fibres are more flexible, they can adapt more easily to the local flow structure and thus sample more frequently the zones of lower vorticity similarly to what observed for point-like particles. Conversely, the fibres are prevented from doing so when increasing 
$\gamma$, because of the rigidity constraint, and consequently a weaker effect is found. As another prominent effect, we observe that the highest peak in the RDF is found for the long fibres (figure 9f), whereas less-remarkable variations are observed for the short and intermediate fibres (figure 9b,d). Such finding can be explained indeed in terms of the combined role of inertia (i.e. Stokes number) and flexibility (i.e. effective compliance), such that for the longest fibres clustering is more favoured.
Figure 10 provides a qualitative insight by collecting some instantaneous visualisations for the iso-dense cases of intermediate stiffness (a,d,g), and for both the most flexible (b,e,h) and most rigid (c, f,i) denser-than-the-fluid fibres. For iso-dense fibres clustering is indeed always negligible, without any peculiar role observed for the fibre's flexibility. Conversely, for denser-than-the-fluid fibres the influence of the latter becomes stronger as the length is increased. As a result, for long fibres a remarkable difference can be noted between the flexible and rigid case (figure 10h,i), with much more intense clusters that can detected for the former.
Two-dimensional snapshots of fibres of different length (a–c short; d–f intermediate; g–i long) and for three different configurations (a,d,g iso-dense with intermediate bending stiffness; b,e,h denser-than-the-fluid with lowest bending stiffness; c, f,i denser-than-the-fluid with highest bending stiffness).
As a final comment, we highlight that clustering does not appear to have a clear role in the turbulence modulation: as shown in § 3.1, it is the mass fraction the most relevant (and arguably unique) control parameter, without a specific effect found, e.g. for the fibre's length, in the backreaction mechanism.
3.2.4. Preferential alignment
As the final step of our investigation, we focus on a peculiar issue of anisotropic particles in turbulence (Voth & Soldati Reference Voth and Soldati2017): the statistical characterisation of the alignment between the fibre orientation and some specific quantities of the turbulent flow (i.e. vorticity and principal directions of the strain rate). The topic has been extensively investigated in the framework of fibres with infinitesimal length, and more recently for finite-size fibres (Pumir & Wilkinson Reference Pumir and Wilkinson2011; Ni et al. Reference Ni, Kramel, Ouellette and Voth2015; Pujara et al. Reference Pujara, Voth and Variano2019, Reference Pujara, Arguedas-Leiva, Lalescu, Bramas and Wilczek2021), pointing out the preferential alignment of sufficiently short fibres with the fluid vorticity and that of longer fibres with the most extensional eigenvector of the strain rate. We note, however, that such theoretical and computational studies still often rely on a one-way coupling assumption and are typically based on the classical Jeffery's model (the latter strictly holding only for fibres shorter than the dissipative length scale). Moreover, fibres are usually assumed to be iso-dense and rigid. Here, we tackle the problem using a four-way coupled simulation approach and span the wide parametric range in terms of fibre's inertia, length and flexibility.
We compare the fibre's local orientation (i.e. considering each segment connecting two adjacent Lagrangian points, and then averaging the results over different segments and fibres) with the local fluid flow to identify the existence of preferential alignment, and to explore how the latter varies with the main properties of the suspension. As in our four-way coupled model the fluid flow is locally perturbed by the presence of the fibre, when computing the vorticity and strain rate we employ a coarse-graining procedure to overcome such effect in the proximity of the Lagrangian points. Specifically, a stencil of 
$7^3$ Eulerian cells surrounding the Lagrangian point is used. To assess the validity of this choice, at first we employ this procedure to compute the alignment between the vorticity and strain rate, a well-known feature of homogeneous turbulent flows (Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987; Tsinober, Kit & Dracos Reference Tsinober, Kit and Dracos1992). Results obtained in a representative configuration with iso-dense fibres (for which the backreaction is overall negligible) are shown in figure 11(a), along with those evaluated using the fully Eulerian framework in the single-phase case: the same qualitative outcome is observed between the Lagrangian coarse-grained and fully Eulerian approaches, with the vorticity resulting mostly aligned with the intermediate eigenvector and orthogonal to the direction of maximum compression, i.e. the third eigenvector (Tsinober et al. Reference Tsinober, Kit and Dracos1992). On the other hand, only a minimal departure may be noted for the alignment with the first eigenvector (corresponding to the direction of maximum extension). Nevertheless, in light of the overall agreement the coarse-grained approach can be safely exploited for the following analysis. Furthermore, the same quantities are shown in figure 11(b) for a case with denser-than-the-fluid fibres which cause a significant modulation of the turbulent flow. In this case, the alteration in terms of vorticity-strain alignment appears not dramatic yet not completely negligible. Specifically, the main difference that can be observed is the mild attenuation of the alignment of vorticity with the intermediate eigenvector. Note that, this effect can play indirectly a role when observing the alignment of inertial fibres with the flow.
P.d.f.s of the alignment between the fluid vorticity and eigenvectors of the strain rate computed in the fully Eulerian framework for the single-phase case (black lines without symbols) and via a Lagrangian coarse-grained approach (brown lines with symbols) for two representative cases (at intermediate length):(a) iso-dense fibres (with lowest bending stiffness); (b) denser-than-the-fluid fibres (with intermediate bending stiffness).
Let us therefore consider the alignment between the fibre's orientation and the strain rate principal directions, starting from the iso-dense cases reported in figure 12. In particular, we first consider the results for short fibres (figure 12a–c) for which, regardless of the variation in the bending stiffness, a strong alignment with both the first and second eigenvectors is found, along with anti-alignment with the third one. Remarkably, this finding is in very good agreement with the experimental results reported by Ni et al. (Reference Ni, Kramel, Ouellette and Voth2015) (see figure 8 therein), specifically for what it concerns the predominant alignment with the first (rather than the second) eigenvector (conversely, as noted by Ni et al. (Reference Ni, Kramel, Ouellette and Voth2015) such qualitative agreement was not found when using the one-way coupling approach). We note the resemblance in the parameters (i.e. fibre length and density, micro-scale Reynolds number) between these cases and the experimental study by Ni et al. (Reference Ni, Kramel, Ouellette and Voth2015); in particular, here the fibre is longer, yet still comparable, than the Kolmogorov length scale. If, on the other hand, we increase the fibre length to the intermediate value (figure 12d–f), we observe that the same phenomenology is retained only for the most flexible case, whereas the alignment with the most extensional direction is progressively lost when increasing the bending stiffness. As a result, for essentially rigid fibres of intermediate length (i.e. belonging to the inertial subrange), the dominant alignment is found with the intermediate eigenvector. Finally, when the length is further increased so that the fibres are comparable with the integral length scale (figure 12g–i), several qualitative modifications are observed: the strongest alignment is always found with the intermediate eigenvector (further increasing with the rigidity of the fibres), whereas the p.d.f. relative to the third eigenvector gets more flattened (both compared with the shorter lengths and when increasing the bending stiffness). Furthermore, for flexible fibres the p.d.f. relative to the first eigenvector shows now a peak at an intermediate value. Such non-monotonic behaviour, which is not observed for other configurations, appears to be conditioned by having a relatively low bending stiffness, and a possible transition can be argued towards a more anti-aligned situation when further increasing the bending stiffness.
P.d.f.s of the alignment between the (local) fibre's orientation and the eigenvectors of the strain rate, for iso-dense fibres of different length (a–c short; d–f intermediate; g–i long) and bending stiffness (increasing from left to right panels). The solid, dashed and dotted lines indicate the alignment with the first, second and third eigenvector of the strain rate, respectively. Colours and symbols are used consistently with the indication of tables 1 and 2.
Figure 13 shows the results concerning the same configurations but for denser-than-the-fluid fibres. The resulting scenario is now dramatically different compared with the iso-dense cases, with two main effects that can be noted. First, with respect to the iso-dense cases, the preferential alignment (or anti-alignment) of fibres with the strain rate turns out to be substantially weaker. This is especially evident when considering short fibres (figure 13a–c) or more flexible fibres (figure 13d,g), where all curves are much more horizontal, with an approximately constant value around 
$1$, thus indicating a tendency towards the situation of random alignment. On one hand, this result could be intuitively expected because the dynamics for highly inertial particles is significantly delayed with respect to the forcing acted by the turbulent flow. As a result, at a given instant the local fibre's orientation and surrounding flow topology are substantially less correlated with a departure from the preferential alignment. However, when increasing the rigidity (i.e. bending stiffness) the situation systematically changes: for all the considered fibre lengths, we observe a stronger alignment with the intermediate eigenvector, along with a more pronounced anti-alignment with the third eigenvector. A difference can be noted instead in terms of the alignment with the first eigenvector between short fibres, which show a weak alignment, and intermediate or long fibres, for which a moderate anti-alignment is noted. Overall, it can be underlined that the inertia of the fibres and the rigidity constraint play an apparently opposite role, i.e. depleting versus promoting the preferential alignment of the dispersed particles with the flow.
P.d.f.s of the alignment between the (local) fibre's orientation and the eigenvectors of the strain rate, for denser-than-the-fluid fibres of different length (a–c short; d–f intermediate; g–i long) and bending stiffness (increasing from left to right panels). The solid, dashed and dotted lines indicate the alignment with the first, second and third eigenvector of the strain rate, respectively. Colours and symbols are used consistently with the indication of tables 1 and 2.
4. Conclusions and outlook
In this work, we have investigated the mutually coupled dynamics of dispersed fibres in homogeneous isotropic turbulence by means of four-way coupled DNS. The interaction between finite-size fibres and the turbulent flow has been tackled in the framework of an immersed boundary method to properly model the coupling between the two phases, thus going beyond the typical assumptions of infinitesimal length and one-way coupling. Fixing a reference single-phase condition, such methodology has been employed to perform an extensive numerical study (in zero-gravity condition) over the main parameters that control the dynamics of both the dispersed and carrier phase, i.e. the fibre's linear density, length and bending stiffness, as well as the number of fibres, in order to simulate and compare a variety of configurations (i.e. iso-dense versus denser-than-the-fluid, short versus long, flexible versus rigid and dilute versus non-dilute). Several outcomes can be outlined from our numerical results, concerning both the modulation of the turbulent flow and some relevant features of the suspension, such as the fibre's flapping states and deformation, as well as the possibility of clustering and preferential orientation.
To start, we have systematically described how the backreaction of the dispersed phase affects the macroscopic behaviour of the carrier flow, with an overall depletion of the turbulent kinetic energy that is reflected by the remarkable variation of the (decreasing) micro-scale Reynolds number. In agreement with reported evidence for other particle-laden flows, such macroscopic effect is controlled by the mass fraction of the suspension, and not by other reference quantities such as the number density or volume fraction. Indeed, we remark that the backreaction effect can be non-negligible already at number densities that correspond to the so-called dilute regime (Butler & Snook Reference Butler and Snook2018). The mass fraction appears to be the main control parameter also when analysing the modification of the energy spectra and the corresponding scale-by-scale energy transfer. When increasing the mass fraction, two robust features are observed from the energy spectra, i.e. an overall large-scale energy depletion along with a relative increase of the small-scale energy content. Such evidence can be better explained by looking at the energy fluxes in Fourier space, where a general tendency is observed as the mass fraction is increased, with the suppression of the nonlinear term along with the increasing relevance of the energy transfer directly associated with the fluid–solid coupling. Finally, we point out that the fibre's bending stiffness appears to have a minor yet systematic influence, with a weaker backreaction for more flexible fibres, as it could be intuitively argued.
Regarding the fibre dynamics, we have first characterised the fibre's deformation in terms of the dominant flapping frequency and possible flapping states, revisiting the phenomenological model of Rosti et al. (Reference Rosti, Banaei, Brandt and Mazzino2018, Reference Rosti, Olivieri, Banaei, Brandt and Mazzino2020a) for the case of non-negligible turbulence modulation. As a result, the same qualitative scenario is observed, with only two possible dynamical regimes: (i) for relatively flexible fibres (such that their natural frequency is lower than the hydrodynamic frequency of turbulence eddies of comparable size) or overdamped (such as iso-dense) fibres, the dynamics is fully controlled by the turbulence structures at the fibre's length scale, i.e. the flapping frequency is locked to that of turbulence, with potential application in exploiting the fibre to measure the two-point statistics of the flow; (ii) for relatively rigid (and inertial) fibres, the fibre's flapping frequency is conversely decoupled from the flow and manifesting the natural structural response to the fluid forcing. From a similar perspective, we have characterised the maximum curvature that fibres can experience while being transported and deformed by the carrier flow, outlining the existence of two different scaling laws whose physical explanation is in partial agreement with that proposed by a previous study dealing with iso-dense fibres (Gay et al. Reference Gay, Favier and Verhille2018). Next, we have looked at the local concentration of the dispersed fibres in order to detect the occurrence of clustering phenomena. Following an approach similar to that usually employed for particle-laden flows, based on evaluating the radial distribution function over pairs of fibres, we observed a systematic increase of clustering when increasing the fibre's flexibility (i.e. decreasing the bending stiffness) and, as expected, for inertial fibres. In particular, we found the strongest clustering for the longest and most flexible, denser-than-the-fluid fibres. Finally, we have assessed the preferential alignment between the fibre's (local) orientation and the strain rate principal directions, comparing our findings with those reported by previous studies in similar conditions, as well as presenting novel results for the flexible or inertial cases. Overall, we observe that both the alignment with the first or second principal direction (i.e. of maximum extension or intermediate extension/compression) and the anti-alignment with the third principal direction (i.e. of maximum compression) typically increase when the fibres become more rigid, whereas remarkably decrease when moving from iso-dense to denser-than-the-fluid fibres.
In conclusion, we outline possible developments of the present study that could provide the motivation for future work. On one hand, exploring the behaviour of finite-size fibre suspensions at substantially higher Reynolds numbers, in particular to assess whether the same or more complex turbulence modulation mechanisms are occurring. On the other hand, further efforts are deserved to better understand the specific features of this kind of suspensions in terms of clustering and, even more importantly, preferential alignment. For inertial fibres, another remarkable point is exploring how the dynamics is altered in finite Froude number conditions (i.e. when the gravitational forces are not negligible). Lastly, another framework deserving further research efforts is surely represented by wall-bounded turbulent flows, which are beyond the scope of the present investigation.
Acknowledgements
S.O. and M.E.R. acknowledge the computational resources provided on the Deigo cluster by the Scientific Computing and Data Analysis section of Research Support Division at OIST and on the Fugaku supercomputer by RIKEN through the HPCI System Research Project (Project IDs: hp210229 and hp210269).
Funding
The research was supported by the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan. A.M. acknowledges the financial support from the Compagnia di San Paolo, Project MINIERA No. I34I20000380007, and from the Project PoC – BUYT MAIH, No. C36I20000140006.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available from the corresponding authors upon reasonable request.
Appendix A. Evaluation of collision events
To account for the effective relevance of the collision model in the performed simulations, we have monitored over time the number of collision events (occurring when two Lagrangian points belonging to different fibres are found below the critical distance 
$\varDelta _{col}$). To properly estimate the relative importance of the fibre-to-fibre interactions, we normalise such number with the total number of Lagrangian point. In figure 14, such quantity is shown for cases with different fibre length, density and bending stiffness. Overall, it can be noted that the occurrence of collisions is always very limited, including cases with the highest concentration, such events typically involving less than 
$1\,\%$ of the ensemble of fibre material points.
Appendix B. Benchmark for the numerical technique
The numerical code used in the present analysis has been extensively validated in the past. For the sake of completeness, in this appendix we present the result of an additional test case. In particular, we consider a passively flapping filament in a uniform flow. For a description of the case set-up the reader is referred to Huang et al. (Reference Huang, Shin and Sung2007) and, in particular, figure 13 therein. Figure 15 shows the comparison between our results and those by the original paper, where one can note the good agreement both in the amplitude and frequency of the flapping oscillation. Such evidence, therefore, further confirms the validity of our numerical findings.
Benchmark results for the case of a flapping filament in uniform flow, showing the time history of the filament's trailing point position by Huang et al. (Reference Huang, Shin and Sung2007) (red circles) and that obtained with the code used for the present study (black line).
Appendix C. List of configurations investigated in the parametric study
This appendix provides additional information on the set-up chosen for our investigation and, in particular, the combinations of parameters that have been considered, which are listed in tables 1 to 3. Note that, the parameters 
$\Delta \tilde {\rho }$, 
$c$ and 
$\gamma$ are given in code units because in the case of finite-size fibres it is not straightforward to choose a representative non-dimensional form. The information is therefore complemented by that regarding the carrier flow: the kinematic viscosity is set to 
$\nu = 1.59 \times 10^{-3}$ and the fluid density is set to 1, whereas the parameters of the external forcing are such that the resulting turbulence is characterised in the single-phase case by a root mean square of the velocity fluctuations 
$u_{rms} \approx 5.14 \times 10^{-1}$ and average dissipation rate 
$\epsilon \approx 4.5 \times 10^{-2}$, and consequently in the Kolmogorov dissipative length scale 
$\eta = (\nu ^3/\epsilon )^{1/4} \approx 1.7\times 10^{-2}$ and the micro-scale Reynolds number 
${Re}_\lambda \approx 120$.
Table 1 reports the cases with almost iso-dense fibres (indicated with empty symbols) and table 2 the cases with iso-dense fibres (denoted instead with filled symbols); for the former, the resulting solid-to-fluid density ratio (which can only be estimated because the fibre's diameter is not explicitly controlled) is around 1 whereas for the latter it is approximately 500. In addition, for a particular configuration, i.e. iso-dense fibres with intermediate length and density, a refined study was performed over the bending stiffness, as detailed in table 3. To give a further overall indication, we can compare the three considered fibre lengths with the dissipative length scale and the domain length (used as a rough estimate for the integral lengthscale), yielding 
$c/\eta = \{ 6, 30, 116 \}$ and 
$c/L = \{ 1.6\times 10^{-2}, 8.0\times 10^{-2}, 0.3 \}$, respectively.
List of parametric combinations considered in the present work (where 
$\Delta \tilde {\rho }$ is the linear density difference between the fibre and the fluid, 
$c$ is the fibre's length, 
$\gamma$ is the fibre's bending stiffness, 
$N$ is the number of fibres dispersed in the triperiodic domain box; all values are given in code units). The corresponding number density 
$n\,c^3$, volume fraction 
$\phi _V$ and mass fraction 
$\mathcal {M}$ are also reported. Each case is denoted by the marker used in the figures throughout the paper. The present table lists the cases for iso-dense fibres, whereas iso-dense fibres are reported in table 2.
Same as table 1, but referring to the cases with iso-dense fibres.
Refined study over 
$\gamma$ for iso-dense fibres.
Appendix D. Derivation of the scale-by-scale energy transfer
In this appendix, we briefly show how to derive the scale-by-scale energy transfer discussed in § 3.1.3, along with identifying the various terms appearing in such balance. For a more detailed explanation, the reader is referred to classical textbooks on turbulence theory (e.g. Pope Reference Pope2000).
As the starting point, we perform the Fourier transform of (2.3) and (2.4),
$$\begin{gather} \partial_t \hat{{\boldsymbol{u}}} + \hat{\boldsymbol{G}} ={-} {\rm i} \boldsymbol{k} \hat{p}/\rho_{f} - \nu k^2 \hat{{\boldsymbol{u}}} + \hat{{\boldsymbol{f}}}_{tur} + \hat{{\boldsymbol{f}}}_{fs}, \end{gather}$$
$$\begin{gather}\boldsymbol{k} \boldsymbol{\cdot} \hat{{\boldsymbol{u}}} = 0, \end{gather}$$
where 
$\widehat{(\cdot)}(\boldsymbol {k},t) = \mathcal {F}\{ (\cdot )({\boldsymbol {x}},t) \}$ denotes the Fourier transform from physical to spectral space, 
$\boldsymbol {k}$ is the wavenumber vector and 
$k$ its magnitude, 
$i$ is the imaginary unit and 
$\boldsymbol {G}$ represents the nonlinear term appearing in the momentum equation. In addition, the same equations can be written for the complex conjugate 
$\hat {{\boldsymbol {u}}}^*$.
We now multiply (D1) by 
$\hat {{\boldsymbol {u}}}^*$, so that the pressure term drops due to the incompressibility constraint. The same applies in the momentum equation for 
$\hat {{\boldsymbol {u}}}^*$ when multiplying by 
$\hat {{\boldsymbol {u}}}$. Then, we sum the two equations for 
$\hat {{\boldsymbol {u}}}$ and 
$\hat {{\boldsymbol {u}}}^*$, thus obtaining an equation for the spectral kinetic energy 
$\hat {E}(\boldsymbol {k},t) \equiv \langle \hat {{\boldsymbol {u}}}^* \boldsymbol{\cdot} \hat {{\boldsymbol {u}}} \rangle / 2$ which reads
\begin{equation} \partial_t \hat{E} = \hat{T} + \hat{V} + \hat{F}_{tur} + \hat{F}_{fs}. \end{equation}Four different contributions can be identified in the equation above:
(i)
$\hat {T} = - \frac {1}{2} \, (\hat {\boldsymbol {G}} \boldsymbol{\cdot} \hat {{\boldsymbol {u}}}^* + \hat {\boldsymbol {G}}^* \boldsymbol{\cdot} \hat {{\boldsymbol {u}}})$ is the energy transfer associated with the nonlinear term;(ii)
$\hat {V} = - 2 \nu k^2 \hat {E}$ is the energy dissipation associated with the viscous term;(iii)
$\hat {F}_{tur} = \frac {1}{2} \, (\hat {{\boldsymbol {f}}}_{tur} \boldsymbol{\cdot} \hat {{\boldsymbol {u}}}^* + \hat {{\boldsymbol {f}}}^*_{tur} \boldsymbol{\cdot} \hat {{\boldsymbol {u}}})$ is the energy input associated with the turbulence forcing;(iv)
$\hat {F}_{fs} = \frac {1}{2} \, (\hat {{\boldsymbol {f}}}_{fs} \boldsymbol{\cdot} \hat {{\boldsymbol {u}}}^* + \hat {{\boldsymbol {f}}}^*_{fs} \boldsymbol{\cdot} \hat {{\boldsymbol {u}}})$ is the energy transfer associated with the fluid–structure coupling.
Next, by isotropically averaging (D3) over the generic sphere of radius 
$k$, we obtain a similar equation for the energy spectrum 
$E(k,t)$,
\begin{equation} \partial_t {E} = {T} + {V} + {F}_{tur} + {F}_{fs}. \end{equation}
Assuming a statistically stationary state, the left-hand side vanishes. To obtain the scale-by-scale energy transfer, we first integrate (D3) from 
$k$ to infinity, yielding
\begin{equation} 0 = \varPi + D' + P + \varPi_{fs}, \end{equation}where the various contributions appearing in the balance are
$$\begin{gather} \varPi(k) = \int_k^\infty T (k) \, \mathrm{d}k, \end{gather}$$
$$\begin{gather}D'(k) = \int_k^\infty V(k) \, \mathrm{d}k, \end{gather}$$
$$\begin{gather}P(k) = \int_k^\infty F_{tur}(k) \, \mathrm{d}k, \end{gather}$$
$$\begin{gather}\varPi_{fs}(k) = \int_k^\infty F_{fs}(k) \, \mathrm{d}k. \end{gather}$$
Again, these quantities are associated with the nonlinear, dissipative, external forcing and fluid–solid coupling terms, respectively. Remarkably, not only the nonlinear transfer 
$\varPi (k)$, but also the transfer associated with the fluid–solid coupling 
$\varPi _{fs} (k)$ are energy fluxes, without any net energy input/output, i.e. 
$\varPi (0) = \varPi _{fs}(0) = 0$, with the overall balance governed only by the external forcing and viscous dissipation, i.e.
\begin{equation} \int_0^\infty F_{tur}(k) \, \mathrm{d}k ={-} \int_0^\infty V(k) \, \mathrm{d}k = \epsilon. \end{equation}
Indeed, exploiting this result we recast the definition of the dissipation using 
$D(k) = - \int _0^k V(k) \, \mathrm {d}k = \epsilon + D'(k)$ in place of 
$D'(k)$, and we finally obtain (3.1). Note that, this choice is meant only for better describing the scale-by-scale energy transfer when plotting the various terms.
Open access
