On the fully coupled dynamics of flexible fibres dispersed in modulated turbulence

The present work investigates the mechanical behaviour of finite-size, elastic and inertial fibres freely moving in a homogeneous and isotropic turbulent flow at moderate Reynolds number. Fully-resolved, direct numerical simulations, based on a finite difference discretisation and the immersed boundary method, are performed to mutually couple the dynamics of fibres and fluid turbulence, allowing to account for the backreaction of the dispersed phase to the carrier flow. An extensive parametric study is carried out over the characteristic properties of the suspension, i.e., fibre's linear density (from neutrally-buoyant to heavy fibres), length (from short fibres comparable with the dissipative scale to long fibres comparable with the integral scale) and bending stiffness (from highly flexible to almost rigid fibres), as well as the concentration (from dilute to non-dilute suspensions). Results reveal the existence of a robust turbulence modulation mechanism that is primarily controlled by the mass fraction of the suspension (with only a minor influence of the fibre's bending stiffness), which is characterised in detail by means of a scale-by-scale analysis in Fourier space. Despite such alteration with respect to the single-phase case due to the non-negligible backreaction, fibres experience only two possible flapping states (previously identified in the very dilute condition) while being transported and deformed by the flow. Additionally, we show that the maximum curvature obeys to different scaling laws that can be derived from the fibre dynamical equation. Finally, we explore the preferential concentration and alignment of fibres within the flow, highlighting the peculiar role of inertia and elasticity.


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 et al. 2011;Guasto et al. 2012;du Roure et al. 2019). 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. 2014;Filella 2015;Brouzet et al. 2021). Other examples concern the locomotion of bacteria or planktonic organisms (Son et al. 2013;Ardekani et al. 2017b;Michalec et al. 2017), the formation of algae aggregates (Verhille et al. 2017), and the manufacturing of paper and composite materials (Parsheh et al. 2005;Mortensen et al. 2008;Marchioli et al. 2010).
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. 2012;Voth & Soldati 2017). 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 et al. 2014(Ni et al. , 2015Pujara et al. 2019Pujara et al. , 2021Oehmke et al. 2021) and inertia (expressed by means of a representative Stokes number) (Bounoua et al. 2018;Kuperman et al. 2019).
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. 2019;Żuk et al. 2021). 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 et al. 2014;Rosti et al. 2018). 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 et al. 2020), although some contributions focused on extensible objects as well, mostly related to suspensions of polymers (Picardo et al. 2018(Picardo et al. , 2020Vincenzi et al. 2021). 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 et al. 2018), and the spatial conformation and deformation of finite-size fibres, i.e., with length within the inertial subrange or comparable to the integral lengthscale (Brouzet et al. 2014;Gay et al. 2018;Sulaiman et al. 2019;Dotto & Marchioli 2019;Dotto et al. 2020;Picardo et al. 2020); further recent developments include the aforementioned modelling of fragmentation processes (Allende et al. 2020;Brouzet et al. 2021;Vincenzi et al. 2021).
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. 2018(Rosti et al. , 2020a and rigid fibres (Brizzolara et al. 2021). For flexible fibres, Rosti et al. (2018Rosti et al. ( , 2020a 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 et al. (2021) 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. (2021), leading to the development of a novel experimental technique, named "Fibre Tracking Velocimetry", able to measure the properties of turbulence at a fixed lengthscale (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 2010;Maxey 2017;Brandt & Coletti 2021). 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 1922). 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 lengthscale.
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 et al. 2010;Gualtieri et al. 2013;Capecelatro et al. 2018;Ardekani & Brandt 2019;Ardekani et al. 2019;Sozza et al. 2020), anisotropic particles and fibres (Andersson et al. 2012;Ardekani et al. 2017a;Olivieri et al. 2020bOlivieri et al. , 2021, as well as droplets or bubbles (Dodd & Ferrante 2016;Freund & Ferrante 2019;Rosti et al. 2019;Cannon et al. 2021). 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. 2013;Dodd & Ferrante 2016;Rosti et al. 2019;Olivieri et al. 2021): 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 turns 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 lengthscale, 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 fiber-to-fiber 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. (2021), with the goal of substantially advancing such findings and offer a systematic and fully exhaustive analysis.
The rest of the paper is structured as follows: in section 2, we introduce the computational methodology used in our investigation along with providing information on the performed parametric study; in section 3, we present and discuss the results to address the two main objectives; finally, concluding remarks are given in section 4.

Methods
This section describes the methodology adopted for the present study, which is essentially the same adopted by Olivieri et al. (2021). First, in section 2.1, we introduce the governing equations and main parameters of the problem. Then, in section 2.2, we describe how the former are solved numerically. Finally, in section 2.3, we provide details on the performed parametric DNS study.

Problem formulation
We consider an ensemble of fibres freely moving in homogeneous and isotropic turbulence (figure 1). The fibre's length 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 modeled as one-dimensional slender objects, i.e., their length is assumed much larger than the diameter , or equivalently we consider an aspect ratio / 1. Moreover, we assume that the fibres are inextensible and focus solely on the elastic contribution due to the bending stiffness (which for a homogeneous fibre is given by the product of the elastic modulus and the second moment of the area). Finally, Δ = s − f is the difference in the linear density † of the fibre and the fluid, i.e., the parameter controlling the inertia of the dispersed objects. In particular, in the limit Δ → 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 is the position of a generic material point belonging to the fibre, ( , ) is the tension enforcing the inextensibility constraint in the filament equation, F fs ( , ) is the fluid-solid interaction forcing and F col ( , ) 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 and time . Freely-moving boundary conditions are imposed at both fibre's ends: X| =0, = 0, X| =0, = 0 and | =0, = 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 nat = √︁ /(Δ 4 ), where ≈ 22.4/ . Dealing with finite-size objects, we neglect the Brownian contribution to the fibre dynamics, which is considered to be small compared to the hydrodynamic forcing.
Fibres are released within a cubic domain of size = 2 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 (1988), randomly injecting energy at the largest scales, i.e., within a low-wavenumber shell with 1 2 (with denoting the wavenumber). The carrier flow is governed by the incompressible Navier-Stokes equations for a Newtonian fluid, which read

3)
· u = 0, (2.4) where u (x, ) and (x, ) are the velocity and pressure fields, respectively, f is the volumetric † Linear and volumetric densities are denoted with and without the tilde, respectively. 6 fluid density and is the kinematic viscosity; f tur (x, ) is the external forcing used to sustain the turbulent flow (Eswaran & Pope 1988) and f fs (x, ) is an additional forcing mimicking the presence of the solid phase (as described in the following); lastly, x denotes the position vector in the Eulerian framework.
The flow and fibre dynamics are coupled by the no-slip condition X = u (X ( , ) , ).

Numerical technique
To solve numerically the fully-coupled FSI problem, we employ the IBM procedure originally proposed by Huang et al. (2007) and subsequently modified by Banaei et al. (2020) and Olivieri et al. (2020bOlivieri et al. ( ,a, 2021. In particular, the mutual interaction between the fluid and solid phase is enforced indirectly by means of the singular force distributions F fs and f fs acting on the fibre and flow, respectively (Peskin 2002;Huang et al. 2007). The fluid velocity at the position of the fibre Lagrangian point, U = u (X ( , ) , ), is obtained by interpolating the fluid velocity at the Eulerian nodes surrounding the Lagrangian point: where is the Dirac delta function. The interpolated velocity U is used to compute the fluid-solid interaction forcing needed to enforce the no-slip condition where is a large negative constant (Huang et al. 2007). Finally, F is spread to the surrounding fluid flow as Both the interpolation and spreading operations feature the Dirac operator. Numerically, this is transposed into the use of the regularized function proposed by Roma et al. (1999). 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 = 2 Δ as the effective diameter. Each fibre is evenly discretized using L Lagrangian points such that the spatial resolution Δ = /( L − 1) is approximately equal to the Eulerian grid spacing Δ ; consequently, the discretized curvilinear coordinate can be evaluated as = Δ , with = 1, . . . , L . The numerical solution of equations (2.1) and (2.2) follows the scheme detailed in Huang et al. (2007) with the difference that the bending term is treated implicitly to allow for a larger timestep (Banaei et al. 2020;Olivieri et al. 2021). Furthermore, a fibre-to-fibre collision model is implemented to prevent that the distance d between Lagrangian points belonging to different fibres goes below the grid spacing Δ and that fibres eventually cross each other. Specifically, we employ the minimal collision model proposed by Snook et al. (2012), where a constant forcing F col = 0d is applied when |d| Δ col . Here, 0 is a free parameter and Δ col = O (Δ ) 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 0 = 1.0 and Δ col = 3Δ . As a result, the influence of collisions turns out to be very limited, as shown in appendix A: fiber-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.

7
Concerning the fluid flow, equations (2.3) and (2.4) are solved using the fractional step method on a staggered grid (Kim & Peskin 2007). 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 finitedifference method for the spatial discretization and the (third-order) Runge-Kutta scheme for the temporal discretization.
The described computational procedure is implemented in the in-house solver Fujin (https://groups.oist.jp/cffu/code). The code is parallelized 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. 2019;Rosti & Brandt 2020;Rosti et al. 2020bRosti et al. , 2021, including in particular suspensions of fibres in both laminar and turbulent flows (Rosti et al. 2018(Rosti et al. , 2020aCavaiola et al. 2020;Olivieri et al. 2020aOlivieri et al. ,b, 2021Brizzolara et al. 2021). 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. (2007).

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 Δ , length and bending stiffness , along with the number of fibres . 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 Δ , corresponding to essentially iso-dense and denser-than-the-fluid fibres, three different values of the fibre's length , which are comparable to the dissipative, inertial and integral characteristic lengthscales of the turbulent flow, three different bending stiffnesses so that we range from highly flexible to essentially rigid fibres. Finally, we varied the numbers of fibres 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 = 0), where the turbulent flow is characterized by a micro-scale Reynolds number Re ≡ rms / ≈ 120, where rms is the root-mean-square velocity and is the Taylor micro-scale. To simulate such configuration, the fluid domain is discretized into a uniform Eulerian grid using 256 3 cells, having verified that the ratio between the Kolmogorov dissipative lengthscale and the grid spacing is /Δ = 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 5 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.

Results
This section presents the results of our study. First, in section 3.1, we focus on the characterisation of the modulation of the turbulent flow due to the dispersed phase. Then, in section 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 2a,b,c shows the dependence of the micro-scale Reynolds number Re as a function of several quantities that can be used as an indicator to characterise the concentration of the suspension: (i) the (nondimensionalized) number density 3 , (ii) the volume fraction , and (iii) the mass fraction M. These three quantities are chosen since 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 2018; Brandt & Coletti 2021). 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.
The (dimensional) number density, defined as = / f (where f = 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 nondimensional with the fibre length , so that 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 3 1, whereas it is classified as semi-dilute for 1/ 3 1/( 2 ). If 1/( 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 2018). 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 2a, 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 it will be 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 , as reported in figure 2b. Here, = s / f is defined as the ratio between the volume of the dispersed phase s = 2 /4 and that of the fluid domain f = 3 . Using this quantity, the results are outlined with a much more defined trend (compared with those for the number density), with Re decreasing as 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 Δ . 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 above, it is therefore natural to consider, as another way to parametrize the backreaction effect, the mass fraction M = s /( s + f ), where s = s s is the (total) mass of the solid dispersed phase and f = f f is the (total) mass of fluid in the domain. As a result, from figure 2c 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 2021), similar evidence is not reported for suspensions of finite-size fibres, such as those investigated in the present work. Focusing on figure 2c, 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 M and eventually reaching, for the highest mass fraction that was tested (M ≈ 89%), a value that is roughly half of that obtained in the single-phase case (Re ≈ 60).
For a more complete characterisation of the turbulence modulation, in figure 2d we also show how the turbulent kinetic energy 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 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 to that caused by the inertia of the dispersed phase.

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 Δ , i.e., iso-dense (left panels) vs denser-than-the-fluid fibres (right panels), and the three fibre lengths (corresponding to each row and increasing from top to bottom) investigated. In each plot, we vary the number of fibres (increasing from light to dark color). 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 since its influence on the turbulence modulation was previously shown to be subleading when compared to the combined variation of the other parameters determining the mass fraction of the suspension. Figure 3 shows that, as the concentration is increased, the energy spectrum is modified with the same qualitative behavior 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 3b with that of long, iso-dense fibres in figure 3e (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 to 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 ) 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 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 section 3.2), from which it can be deduced the extensive (nondimensional) range that has been considered (cf. figure 7). Nonetheless, it has also to be pointed out that the variation with appears overall very limited, confirming the secondary effect of flexibility in altering the carrier flow.

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 ( ) + Π( ) + Π fs ( ) + ( ) = , (3.1) where the various terms appearing on the left hand side are the production rate (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 Π associated with the nonlinear term, the energy flux Π fs associated with the fluid-solid coupling, and the viscous dissipation . Note that the latter is defined such that the turbulent energy dissipation rate appears on the right hand side and the overall balance can be visualized more conveniently. For a derivation of  14 equation (3.1), see appendix D. In the single-phase case Π fs ≡ 0, and the dominant terms in the balance are the nonlinear term Π and the viscous dissipation when approaching the lowest and highest wavenumbers, respectively (Pope 2000).
In this analysis, we focus on the cases with = 10 3 . However, the generality of the results is not restricted since 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 isodense fibres (with different fibre lengths), shown in figure 5a,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 Π and a negligible energy transfer associated with the fluid-structure coupling Π fs . This situation is accompanied by the predominance at high wavenumbers of the dissipation term which eventually saturates to for increasing . The same applies also to fibres of intermediate length (figure 5c), with only a very modest increase of Π 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 maximized at an intermediate lengthscale; as increases, the dissipation progressively becomes dominant, although the saturation is less evident due to the enhanced small-scale energy transfer. Remarkably, the mass fraction 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.
The alteration gets substantially more dramatic when considering the same three cases but for denser-than-the-fluid fibres, reported in figure 5b,d,f. Already for short fibres (figure 5b) a similar behavior 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 the one with iso-dense long fibres, as previously mentioned. 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, while 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  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 , on the backreaction and turbulence modulation. Figure 6 reports the energy flux from the nonlinear (left panels) and fluid-solid coupling (right panels) 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 the one 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.

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. (2018Rosti et al. ( , 2020a in the context of very dilute suspensions, and more recently by Olivieri et al. (2021) 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. 2020a;Olivieri et al. 2021): (i) for sufficiently flexible and/or iso-dense fibres (i.e., nat / tur 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., nat / tur 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 generalize such scenario in the case of non-negligible backreaction (Olivieri et al. 2021). As a starting point, figure 7 shows the behaviour of the flapping frequency 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. (2020a)) 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 nat / tur 1, the flapping frequency is controlled by, and therefore scales with, the turbulent eddy frequency (i.e., flap = tur ), whereas for nat / tur 1 the flapping frequency is given by the natural one (i.e., flap = 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., flap = 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 section 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 tur = √ 2 / , where 2 = ( ) 2 is the second-order structure function of the longitudinal velocity difference evaluated at the fibre's lengthscale and = 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 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 (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. 2021). On the other hand, it can be pointed out that the alteration of tur does not have a primary role in the resulting scaling laws (which remain the same already found for the dilute case).
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. (2018) 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 E el = 1 =1 ∫ 0 1 2 2 ( ) , where is the local bending curvature of the -th fibre, peaks when nat / tur ≈ 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 generalizing the theoretical framework previously proposed in the case of negligible backreaction.

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. 2018;Brouzet et al. 2021). Therefore, we focus on understanding how max behaves as a function of the relevant mechanical parameters. We compute the maximum fibre's curvature max as follows: we first evaluate the maximum of ( ) 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 , 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.
To proceed, we adopt an approach similar to that proposed by Gay et al. (2018) 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. (2018), 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., oneway coupling and inertial subrange scaling, are retained (and will be commented on later). For the forcing term we can write F ∼ ( ) 1/3 , where is the dynamic viscosity, while the (linear) bending term can be dimensionally expressed as 4 X ∼ max / 2 . Moreover, as shown by Marheineke & Wegener (2006);Gay et al. (2018), from the tension term ( X) one can isolate the contribution 2 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 ( X) ∼ 3 max . Remarkably, the cubic term is subleading with respect to the linear one for max 1, while the opposite occurs for max 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 max ∼ /( 1/3 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 max ∼ /( 1/3 10/3 ) −1/3 . Figure 8 shows the results from our numerical simulations along with the expected scaling laws, both for the iso-dense (top) and denser-than-the-fluid fibre cases (bottom). 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 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.

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 1994;Bec et al. 2007) and with a specific relevance in the context of fibres aggregation (Lundell et al. 2011;Verhille et al. 2017). Different metrics have been proposed to quantify the tendency of particles to preferentially accumulate along with sampling specific regions of the flow (Brandt & Coletti 2021), among which is the radial distribution (or pair correlation) function (RDF) (Saw et al. 2008;Salazar et al. 2008;Olivieri et al. 2014). Such 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 = 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 1994). 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 to do so when increasing , 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 ones ( 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 (left panels), and for both the most flexible (centre) and most rigid (right) 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. As a final comment, we highlight that clustering does not appear to have a clear role in the turbulence modulation: as shown in section 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.

Preferential alignment
As the final step of our investigation, we focus on a peculiar issue of anisotropic particles in turbulence (Voth & Soldati 2017): 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 2011;Ni et al. 2015;Pujara et al. 2019Pujara et al. , 2021, 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 lengthscale). 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. Since 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. 1987;Tsinober et al. 1992). Results obtained in a representative configuration with iso-dense fibres (for which the backreaction is overall negligible) are shown in figure 11a, 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. 1992). 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 11b 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.
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,b,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. (2015) (see Fig. 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. (2015) 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. (2015); in particular, here the fibre is longer, yet still comparable, than the Kolmogorov lengthscale. If, on the other hand, we increase the fibre length to the intermediate value ( figure 12d,e,f ), we observe that the same phenomenology is retained only for the most flexible case, while 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 lengthscale (figure 12g,h,i), several qualitative modifications are observed: the strongest alignment is always found with the intermediate eigenvector (further increasing with the rigidity of the fibres), while the PDF 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 PDF 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. Figure 13 shows the results concerning the same configurations but for denser-thanthe-fluid fibres. The resulting scenario is now dramatically different compared with the iso-dense cases, with two main effects that can be noticed. First, with respect to the isodense 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,b,c) or more flexible ones ( 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 since 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 noticed. Overall, it can be underlined that the inertia of the fibres and the rigidity constraint play an apparently opposite role, i.e., depleting vs promoting the preferential alignment of the dispersed particles with the flow.

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 vs denser-than-the-fluid, short vs long, flexible vs rigid, and dilute vs 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 2018). 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. (2018Rosti et al. ( , 2020a 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 lengthscale, 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. 2018). 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 particleladen 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 gravitional 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.  Table 1: List of parametric combinations considered in the present work (where Δ is the linear density difference between the fibre and the fluid, is the fibre's length, is the fibre's bending stiffness, is the number of fibres dispersed in the triperiodic domain box; all values are given in code units). The corresponding number density 3 , volume fraction and mass fraction M are also reported. Each case is denoted by the marker used in the figures throughout the paper. The present table shows the cases for iso-dense fibres, whereas iso-dense fibres are reported in table 2.
the single-phase case by a r.m.s. of the velocity fluctuations rms ≈ 5.14 × 10 −1 and average dissipation rate ≈ 4.5 × 10 −2 , and consequently in the Kolmogorov dissipative lengthscale = ( 3 / ) 1/4 ≈ 1.7 × 10 −2 and the micro-scale Reynolds number Re ≈ 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 since the fibre's diameter is not explicitly controlled) is around 1 while for the latter it is approximately 500. Additionally, 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 lengthscale and the domain length (used as a rough estimate for the integral lenghtscale), yielding / = {6, 30, 116} and / = 1.6 × 10 −2 , 8.0 × 10 −2 , 0.3 , respectively.

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 section 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 (2000).
As the starting point, we perform the Fourier transform of equations (2.3) and (2.4), u +Ĝ = − kˆ/ f − 2û +f tur +f fs , (D 1) k ·û = 0, (D 2) where( ·) (k, ) = F {(·)(x, )} denotes the Fourier transform from physical to spectral space, k is the wavenumber vector and its magnitude, is the imaginary unit and G represents the nonlinear term appearing in the momentum equation. Additionally, the same equations can be written for the complex conjugateû * .

(D 3)
Four different contributions can be identified in the equation above: = − 1 2 (Ĝ ·û * +Ĝ * ·û) is the energy transfer associated with the nonlinear term; = −2 2ˆi s the energy dissipation associated with the viscous term; tur = 1 2 (f tur ·û * +f * tur ·û) is the energy input associated with the turbulence forcing; fs = 1 2 (f fs ·û * +f * fs ·û) is the energy transfer associated with the fluid-structure coupling. Next, by isotropically averaging equation (D 3) over the generic sphere of radius , we obtain a similar equation for the energy spectrum ( , ), = + + tur + fs .
(D 4) Assuming a statistically stationary state, the left-hand-side vanishes. To obtain the scale-byscale energy transfer, we first integrate equation (D 3) from to infinity, yielding 0 = Π + + + Π fs , (D 5) where the various contributions appearing in the balance are Again, these quantities are associated with the nonlinear, dissipative, external forcing and fluid-solid coupling terms, respectively. Remarkably, not only the nonlinear transfer Π( ), but also the transfer associated with the fluid-solid coupling Π fs ( ) are energy fluxes, without any net energy input/output, i.e., Π(0) = Π fs (0) = 0, with the overall balance governed only by the external forcing and viscous dissipation, i.e.,