Scaling and intermittency in turbulent flows of elastoviscoplastic fluids

Non-Newtonian fluids have a viscosity that varies with applied stress. Elastoviscoplastic fluids, the elastic, viscous and plastic properties of which are interconnected in a non-trivial way, belong to this category. We have performed numerical simulations to investigate turbulence in elastoviscoplastic fluids at very high Reynolds-number values, as found in landslides and lava flows, focusing on the effect of plasticity. We find that the range of active scales in the energy spectrum reduces when increasing the fluid plasticity; when plastic effects dominate, a new scaling range emerges between the inertial range and the dissipative scales. An extended self-similarity analysis of the structure functions reveals that intermittency is present and grows with the fluid plasticity. The enhanced intermittency is caused by the non-Newtonian dissipation rate, which also exhibits an intermittent behaviour. These findings have relevance to catastrophic events in natural flows, such as landslides and lava flows, where the enhanced intermittency results in stronger extreme events, which are thus more destructive and difficult to predict. Elastoviscoplastic fluids combine solid- and liquid-like behaviour depending on applied stress. Simulations of elastoviscoplastic fluids at high Reynolds number now show that plasticity plays a key role in the turbulent flows seen in these systems, leading for example to intermittency.


INTRODUCTION
Many fluids in nature and industry exhibit a non-linear relationship between shear stress and shear rate, which is referred to as non-Newtonian behaviour.Several non-Newtonian features can exist, and they are often present simultaneously.Here, we focus on the so-called elastoviscoplastic (EVP) fluids, which are fluids with elastic, viscous, and plastic properties.EVP materials combine solid-like behaviour and fluid-like response depending on the value of the applied stress: they behave like a solid when the applied stress is below a critical value known as the yield stress, and flow like a liquid otherwise [1].The elastic nature of these materials is present in their solid as well as liquid states [2].Such fluids are common in everyday life (e.g.toothpaste, jam, cosmetics, mud), and turbulent flows of EVP fluids are found in many industrial processes, including sewage treatment, crude oil transportation, concrete pumping, and mud drilling [3][4][5], and they are also found in nature as landslides and lava flows [6,7].
A great deal of work has been done in the past to properly characterize the viscoelastic behaviour of a fluid in both laminar and turbulent flows [8][9][10][11][12][13], while the e↵ect of plasticity has been studied mainly in low-Reynolds number of the applied stress: they behave like a solid when the applied stress is below a critical value known as the yield stress, and flow like a liquid otherwise [1].The elastic nature of these materials is present in their solid as well as liquid states [2].Such fluids are common in everyday life (e.g.toothpaste, jam, cosmetics, mud), and turbulent flows of EVP fluids are found in many industrial processes, including sewage treatment, crude oil transportation, concrete pumping, and mud drilling [3][4][5], and they are also found in nature as landslides and lava flows [6,7].
A great deal of work has been done in the past to properly characterize the viscoelastic behaviour of a fluid in both laminar and turbulent flows [8][9][10][11][12][13], while the effect of plasticity has been studied mainly in low-Reynolds number laminar conditions [1,14,15].Little is known about the plastic behaviour of an EVP fluid in turbulence; Rosti et al. [16] studied for the first time a turbulent channel flow of an EVP fluid, finding that the shape of the mean velocity profile controls the regions where the fluid is unyielded, forming plugs around the channel centreline that grow in size as the yield stress increases, similar to what is observed in a laminar condition.However, the presence of the plug region has an opposite effect on drag for laminar and turbulent flow configurations, resulting in drag reduction in the turbulent case and drag increase in the laminar one; the turbulent drag behaviour is due to the tendency of the turbulent flow to relaminarize, overall leading to a strongly non-linear relation between yield stress and drag coefficient.The simulation results were then employed by Le Clainche et al. [17] using high-order dynamic mode decomposition to study the near-wall dynamics, comparing them to those in Newtonian and viscoelastic fluids.Their work revealed that both elasticity and plasticity have similar effects on the near-wall coherent structures, where the flow is characterized by long streaks disturbed for short periods by localized perturbations.A recent experimental study by Mitishita et al. [18] on a turbulent duct flow of Carbopol solution de-facto verified the numerical results obtained by Rosti et al. [16] on the effect of plasticity on the mean flow profile and Reynolds stresses.Additionally, they observed an increase in the energy content at large scales and a decrease at small scales, when compared to a Newtonian fluid.Mitishita et al. reported a −7/2 scaling in the energy spectra at high wavenumbers during Carbopol flows compared to −5/3 scaling in the case of water flows.The newly observed scaling was attributed either to the decrease in the inertial effect in the presence of Carbopol solutions that shrinks the inertial range of scales, since the Reynolds numbers are much lower than in water flows, or to the elastic effects that become significant at large wavenumbers where the fluid experiences high frequencies.Moreover, the shear-thinning effects that Carbopol solutions exhibit affected the anisotropy and the overall flow behaviour.The elastic and shear thinning effects are rheological features of Carbopol solutions and can not be eliminated experimentally.
Homogeneous and isotropic turbulent flows have long been a focus of turbulence research for their simple theoretical analysis and the generality of their results.To this end, as has been extensively done in the past for viscoelastic flows, we study the tri-periodic homogeneous flow, where the celebrated K41 theory by Kolmogorov [19], can be directly applied to a classical Newtonian fluid.In this work, we study for the very first time a homogeneous isotropic turbulent (HIT) flow of an elastoviscoplastic fluid at high Reynolds number, as shown in Fig. 1.We aim to answer the following fundamental question: how does the Kolmogorov theory change when the fluid is elastoviscoplastic?We will mainly focus on its plastic behaviour and investigate how the yield stress affects the multiscale energy distribution and balance, and how the turbulent energy cascade is altered by the fluid's plasticity.Our results show profound modifications of the classical picture predicted by the K41 theory for Newtonian fluids, with the emergence of a new scaling range, the dominance of the non-Newtonian flux and dissipation at small and intermediate scales, and enhanced intermittency of the flow.

RESULTS
To investigate the problem at hand, we perform massive three-dimensional direct numerical simulations (DNS) of HIT where we solve the flow equations fully coupled with the constitutive equation of the EVP fluid, within a tri-periodic domain of size L, using 1024 grid points per side, as described in more detail in the Methods section.The flow is controlled by four main parameters: the Reynolds number Re Λ , the Weissenberg number W i Λ , the viscosity ratio α, and the Bingham number Bi Λ , all based on the root mean square velocity fluctuations u and Taylor's microscale Λ.We use the definitions Re Λ ≡ ρu Λ/µ t , W i Λ ≡ λµ t /ρΛ 2 0 , α = µ n /µ t , and Bi Λ ≡ τ y Λ 0 /µ t u 0 , where ρ is the fluid density, µ t ≡ µ f + µ n is the total dynamic viscosity with µ f being the fluid viscosity and µ n the non-Newtonian one, λ is the relaxation time, and τ y is the yield stress, and subscript 0 denotes quantities from the Bi Λ = 0 case.The Reynolds number describes the ratio of inertial to viscous forces, and we limit our analysis to high-Reynolds number flows, achieving a Taylor micro-scale Reynolds number Re Λ ≈ 435 for the Newtonian flow, at which statistics of the flow have been found to be universal and exhibiting a proper scale separation, with an extensive inertial range of scales extended to almost two decades of wavenumbers.The Reynolds number explored here is the highest ever reached in DNS of HIT of non-Newtonian fluids.The Weissenberg number describes the ratio of elastic to viscous   3: Scale-by-scale energy balance for BiΛ = 0 in black, BiΛ = 2.5 in dark green, BiΛ = 12.5 in light green, and BiΛ = 25 in orange.We plot the energy flux of the non-linear convective term Π using dashed lines, the solvent dissipation D using dotted lines, and the non-Newtonian contribution N using solid lines.Each term is normalized by the total dissipation rate t .N grows at intermediate and small scales when BiΛ is increased, eventually becoming the dominant contribution.forces, and here we limit the analysis to W i Λ 1, (i.e., W i Λ ≈ 10 −3 ), to ensure that elastic effects are sub-dominant and all the observed changes are due to plasticity.We also fix the value of α = 0.1 to represent a dilute concentration of polymers, in accordance with prior works on the subject [16,20].Thus, the key control parameter we vary is Bi Λ , which describes the ratio of the yield stress to the viscous stress, and thus correlates with the prevalence of unyielded regions.
Fig. 2 depicts the turbulent kinetic energy spectra of the cases analysed.The Bi Λ = 0 case is similar to the Newtonian case shown in Fig. S1 of the Supplementary Information, confirming that the effect of elasticity is subdominant and can be ignored.A clear E ∼ κ −5/3 range is visible for more than one decade, showing Re Λ is high enough to achieve scale separation, with the spectra exhibiting an inertial range of scales followed by a dissipative range.As Bi Λ increases, the inertial range is limited to the large scales (small wavenumbers κ), with the energy increasing at the large scales while decreasing at the small scales.A clear deviation from the Kolmogorov scaling becomes noticeable for Bi Λ > 1, resulting in the emergence of a new apparent scaling of E ∼ κ −2.3 that is shown more clearly by plotting compensated energy spectra (as shown in Fig. S3 of the Supplementary Information).The difference in scaling between the experimental work (−7/2) [18] and the current study (−2. the cases where Bi Λ < 1, Re Λ remains relatively unaltered, with Φ always close to zero, whereas when Bi Λ further increases, the micro-scale Reynolds number Re Λ and the volume Φ of the unyielded regions rapidly increase with a similar trend.
To fully characterize the change in the energy spectra, we study the turbulent kinetic energy balance, which in wavenumber space can be expressed as where F inj is the turbulence production introduced by the external forcing (injected at the largest scale κ L ≡ 2π/L); Π, D, and N are the non-linear energy flux, the fluid dissipation, and the non-Newtonian contribution, respectively.In addition to the classical bulk fluid dissipation rate f , here we have a non-Newtonian dissipation n which is the rate of removal of turbulent kinetic energy from the flow due to the non-Newtonian extra stress tensor (see the Supplementary Information for a derivation of this equation).Fig. 3 shows the turbulent kinetic energy balance for a few representative values of Bi Λ .When comparing with Fig. S1b of the Supplementary Information, the Bi Λ = 0 case closely follows the classical Newtonian turbulent flow, wherein energy is carried by Π from the large to small scales before being dissipated by the fluid viscosity D. The contribution of the non-linear convective term Π, which appears as an almost horizontal plateau at relatively large scales, progressively decreases with Bi Λ and shrinks towards larger scales, consistently with the reduction of the extension of the inertial range observed in Fig. 2. The reduced energy flux with Bi Λ is also accompanied by a decrease of the fluid dissipation D, which are instead compensated by the increase of non-Newtonian contribution N .At the small scales (large κ), the relative importance of the non-Newtonian contribution increases with Bi Λ , becoming comparable to the fluid dissipation for Bi Λ ≈ 2.5 and eventually becoming the dominant term for Bi Λ 12.5, corresponding to the emergence of the new scaling in the energy spectrum shown in Fig. 2; indeed, the non-Newtonian contribution can be interpreted as a combination of a pure energy flux (giving rise to the new scaling region) and a pure dissipative term, as recently suggested by Rosti et al. [21].Regarding the direction of energy flux, Fig. S4 in the Supplementary Information shows that we have a direct cascade of energy from large to small scales for all Bi Λ [22,23].
We extend the analysis done in the spectral domain, by computing the longitudinal structure functions defined as S p (r) = (∆u(r)) p , where p is the order of the structure function and ∆u(r) = u(x + r) − u(x) is the difference in the fluid velocity across a length scale r, projected in the direction of r.According to K41, S p (r) ∼ ( t r) p/3 ; however, when the structure functions are displayed as a function of r, as shown in Fig. 4a, they deviate from the K41 prediction as p increases.This phenomenon is thought to be due to the intermittency of the flow, i.e., extreme events which are localised in space and time, and thereby break Kolmogorov's hypothesis of self similarity in the inertial range [24].Intermittency results in the scaling exponent of r being a non-linear concave function of p (instead of p/3) [25].For the EVP fluid, two scaling regions appear at large Bi Λ , with scaling consistent with those from the energy spectra, and with intermittency present in both scaling regions.The role of intermittency in the scaling exponents can be better appreciated when the structure functions are displayed in their extended self-similarity form, obtained by plotting one structure function against another [26].In Fig. 4b, S 4 and S 6 are plotted against S 2 for all Bingham numbers considered.We note a clear power-law scaling, which deviates from Kolmogorov's prediction, even for the Bi Λ = 0 case shown in black.The departure from Kolmogorov's prediction progressively grows as the plasticity of the fluid increases, suggesting that the flow becomes more intermittent due to its plasticity.This becomes more obvious when we plot S n compensated by the intermittency correction at Bi Λ = 0 against S 2 (see Fig. S5 in the Supplementary Information).Also, intermittency appears to act equally in the two scaling regions present at large Bi Λ .
Intermittency originates from the multifractal nature of the turbulent dissipation rate [24].For Newtonian fluids, this can be quantified by the multifractal spectrum of the energy dissipation rate f [24,27], which we report in the inset of Fig. 4b.The graph demonstrates that F (α) is nearly identical for all Bi Λ cases except for minor variations at small and large values of α.This implies that the fluid dissipation rate is not the cause of the enhanced intermittency observed in the extended self-similarity analysis.
In the present flow, the turbulent kinetic energy is dissipated by two different terms f and n seen in Fig. 1; hence, we investigate their respective behaviour by looking at their probability distribution functions in Fig. 5.We name the non-Newtonian contribution n a "dissipation" because on average it removes energy from the flow, giving rise to the positive-skewed distributions in Fig. 5b; however, unlike the fluid dissipation, it can take positive or negative values at particular locations in space and time.Fig. 5a shows that the distribution of f narrows as Bi Λ increases [28]; on the other hand, from Fig. 5b, we see that that the distribution of n significantly broadens as Bi Λ increases.Since the non-Newtonian dissipation becomes dominant for the largest Bi Λ (as shown in Fig. 3), we can thus infer that the extreme values of n are indeed the source of the enhanced intermittency observed from the structure function analysis in Fig. 4.

DISCUSSION
By means of unprecedented high-Reynolds-number DNS of an elastoviscoplastic fluid, we have shown that plastic effects significantly alter the classical turbulence predicted by the Kolmogorov theory for Newtonian fluids.
We have proved that the non-Newtonian contribution to the energy balance becomes dominant at intermediate and small scales for large Bingham numbers, inducing the emergence of a new intermediate scaling range in the energy spectra between the Kolmogorov inertial and the dissipative ranges, where energy spectrum decays with a −2.3 exponent.Interestingly, this exponent has been recently found for turbulence of viscoelastic fluids at large Reynolds and Weissenberg number [21,29], suggesting a possible similarity among plastic and elastic effects on the turbulent cascade.This similarity in the scaling behaviour of the two cases could be attributed to a similar interaction mechanism in the Navier-Stokes equation between the convective and extra stress terms.It is also worth noting that in the context of viscoelastic flows at high Weissenberg number, an exponent less than or equal to -3 has been widely reported in the past [8]; however, this is only found at relatively lower Reynolds number than investigated here or explored in recent experimental and numerical work [21,29].The present work appears to be the first to report the −2.3 scaling in turbulent flows of highly plastic EVP fluids, and further studies on the size and distribution of the unyielded regions could shed more light on the origin of the newly found scaling.
We have also shown that the flow in the presence of plastic effects is more intermittent than in a Newtonian fluid, due to the combination of the classical intermittency originating from the multifractal nature of the turbulent dissipation rate, which remains substantially unaltered, and a new plastic contribution which instead grows with the Bingham number.A direct consequence of this result is that intermittency corrections for an elastoviscoplastic fluid are non-universal and dependent of the flow configuration, differently from viscoelastic flows.These results are relevant for several catastrophic natural flows with high plasticity, e.g., lava flows and landslides [30].Our findings explain why such flows are usually found to be intermittent and frequently aggressive, resulting in more damage.The non-universality of the flow intermittency in elastoviscoplastic fluids reflects also in an increased difficulty in their modelling.
The flow under investigation is governed by a system of a scalar, a vector and a tensorial equation, these are the incompressibility constraint, the conservation of momentum, and the constitutive equation for the non-Newtonian extra stress tensor, respectively.The incompressibility constraint and the momentum conservation equations can be written as where u is the fluid velocity, p is the pressure, ρ is the density, and µ f is the fluid dynamic viscosity.The term f inj represents the external force used to sustain turbulence; here we consider the Arnold-Beltrami-Childress (ABC) flow with forcing where i, j, k are the Cartesian unit vectors, A, B, and C are real parameters, and the flow has periodicity L in x, y, and z.In our simulations, we choose A = B = C and use an appropriate value of µ f to give a micro-scale Reynolds number Re Λ ≈ 435 for the Newtonian flow.The last term in equation 3 is defined as f evp ≡ ∇ • τ , where τ is the non-Newtonian extra stress tensor of the EVP fluid.We adopt the constitutive model proposed by Saramito [31] to express the evolution of the extra stress tensor which can be written as where ( ∇ .) denotes the upper convected derivative, i.e., µ n is the non-Newtonian dynamic viscosity, τ d is the magnitude of the deviatoric part of the stress tensor τ d ≡ τ − tr(τ )I/3, and I is the identity tensor, i.e., τ d = 1 2 (τ d : τ d ).Before yielding, i.e., τ d ≤ τ y , the model predicts only recoverable Kelvin-Voigt viscoelastic deformation, while after yielding, i.e., τ d > τ y , it predicts Oldroyd-B viscoelastic behaviour.This transition occurs in a continuous manner.There are other EVP models that take into account shear-thinning [32] or thixotropic behaviour [33]; however, we chose the one described above for its simplicity and the least number of involved parameters.Also, this model proved able to capture the main flow characteristics in a turbulent channel flow [16,18].

Numerical method
We use the in-house flow solver Fujin (https://groups.oist.jp/cffu/code) to solve the governing equations numerically on a staggered uniform Cartesian grid; velocities are located on the cell faces, while pressure, stresses, and the other material properties are located at the cell centres.The second-order central finite difference scheme is used for spatial discretisation except for the advection term that comes from the upper convective derivative in Eq. ( 5) where the fifth-order WENO (weighted essentially non-oscillatory) scheme is adopted [34].The second-order Adams-Bashforth scheme coupled with a fractional step method [35] is used for the time advancement of all terms except for the non-Newtonian extra stress tensor, which is advanced with the Crank-Nicolson scheme.To enforce a divergence-free velocity field, a fast Poisson solver based on the Fast Fourier Transform (FFT) is used for the pressure.The domain decomposition library 2decomp (http://www.2decomp.org)and the MPI protocol are used to parallelize the solver.The evolution equation of the extra EVP stress is formulated and solved using the log-conformation method [36] to ensure the positive-definiteness of the conformation tensor.The fluid domain is a periodic cubic box of length L discretized using 1024 grid points per side, resulting in a large grid resolution sufficient to represent the fluid properties at all the scales of interest with η/∆x = O(1), where η is the Kolmogorov length-scale, and ∆x is the grid spacing.

= 25 FIG. 1 :FIG. 1 :
FIG. 1: Instantaneous colormaps of the turbulent fluid dissipation ✏ f in homogeneous isotropic turbulence of an EVP fluid at di↵erent Bingham numbers.Yielded regions are shown with the black-red-yellow colorscale, while unyielded regions with the black-gray-white one.FIG.1: Instantaneous colourmaps of the turbulent fluid dissipation f in homogeneous isotropic turbulence of an EVP fluid at different Bingham numbers.Yielded regions are shown with the black-red-yellow colourscale, while unyielded regions with the black-gray-white one.

FIG. 2 :
FIG.2: Turbulent kinetic energy spectra of EVP flows with various Bingham numbers, plotted in colours from dark to light; BiΛ = 0, 0.0025, 0.025, 0.25, 2.5, 12.5 and 25 are plotted in black, purple, dark blue, light blue, dark green, light green, and orange respectively.The expected Kolmogorov scaling for a Newtonian fluid is shown by a grey dashed line, while the dash-dotted line represents an apparent new non-Newtonian scaling E ∼ κ −2.3 which emerges at large BiΛ.The inset of the figure reports how the mean values of the micro-scale Reynolds number ReΛ (plotted using squares on the right axis) and the volume fraction of the unyielded regions Φ (plotted using diamonds on the left axis) vary as a function of BiΛ.Error bars report the standard deviation of ReΛ in time, measured using 10 3 samples.Plastic effects start to appear for BiΛ 1, suggesting that Λ is the relevant length scale of the problem.
FIG.3: Scale-by-scale energy balance for BiΛ = 0 in black, BiΛ = 2.5 in dark green, BiΛ = 12.5 in light green, and BiΛ = 25 in orange.We plot the energy flux of the non-linear convective term Π using dashed lines, the solvent dissipation D using dotted lines, and the non-Newtonian contribution N using solid lines.Each term is normalized by the total dissipation rate

3 F
FIG.4:(a): Dependence of the longitudinal velocity structure functions S2 (pluses), S4 (squares), and S6 (hexagons) on the separation distance r.Symbol colour denotes BiΛ and is the same as in Fig.2.Dashed lines show scalings predicted by K41, and dash-dotted lines show scalings predicted using the new non-Newtonian scaling E ∼ κ −2.3 .(b): The two scalings collapse onto a single line when we plot the structure functions S4 and S6 in extended self-similarity form, i.e., against S2.To easily see changes in gradient, we have normalized the structure functions by their values at r ≈ Λ.The dotted line shows a best fit through the data for BiΛ = 0, which deviates from the K41 prediction (dashed line) due to intermittency.Increasing BiΛ further increases this deviation.The inset in (b) reports the multifractal spectrum of the energy dissipation rate carried out by the fluid f .

FIG. 5 :
FIG.5: Probability distribution function (PDF) of (a) the fluid dissipation rate f and of (b) the non-Newtonian dissipation rate n averaged over time, and normalised by 0 the total dissipation of the BiΛ = 0 flow.As BiΛ increases, the PDF of the fluid dissipation rate f slightly narrows, while the PDF of the non-Newtonian contribution n significantly widens.