Effect of finite Weissenberg number on turbulent channel flows of an elastoviscoplastic fluid

Abstract Direct numerical simulations are carried out to study the effect of finite Weissenberg number up to $Wi=16$ on laminar and turbulent channel flows of an elastoviscoplastic (EVP) fluid, at a fixed bulk Reynolds number of $2800$. The incompressible flow equations are coupled with the evolution equation for the EVP stress tensor by a modified Saramito model that extends both the Bingham viscoplastic and the finite extensible nonlinear elastic-Peterlin (FENE-P) viscoelastic models. In turbulent flow, we find that drag decreases with both the Bingham and Weissenberg numbers, until the flow laminarises at high enough elastic and yield stresses. Hence, a higher drag reduction is achieved than in the viscoelastic flow at the same Weissenberg number. The drag reduction persists at Bingham numbers up to 20, in contrast to viscoplastic flow, where the drag increases in the laminar regime compared with a Newtonian flow. Moreover, elasticity affects the laminarisation of an EVP flow in a non-monotonic fashion, delaying it at lower and promoting it at higher Weissenberg numbers. A hibernation phenomenon is observed in the EVP flow, leading to large changes in the unyielded regions. Finally, plasticity is observed to affect both low- and high-speed streaks equally, attenuating the turbulent dissipation and the fragmentation of turbulent structures.


Introduction
The manuscript considers unsteady and turbulent channel flows of non-Newtonian viscoelastic and elastoviscoplastic fluids. The dynamics of non-Newtonian turbulence has attracted a growing interest due to numerous industrial applications dealing † Email address for correspondence: outi@mech.kth.se with non-Newtonian fluids. Real industrial and natural fluids often present several non-Newtonian effects at the same time, such as plasticity (yield stress) and elasticity. Turbulent flows of yield-stress fluids are found in several industries such as petroleum, paper, mining and sewage treatment (Hanks 1963(Hanks , 1967Maleki & Hormozi 2018). Highly inertial flows of elastoviscoplastic fluids can also be found in geophysical applications, such as mudslides and the tectonic dynamics of the Earth (Oishi, Martins & Thompson 2017).
Yield-stress fluids behave like solids below a local stress threshold, and flow like liquids above this threshold. Assuming that the material is a rigid solid at stresses lower than the yield stress, purely viscoplastic models are obtained, such as the Bingham model (Bingham 1922), where the solvent viscosity of the yielded fluid flow follows a Newtonian law, and the Herschel-Bulkley model (Herschel & Bulkley 1926), where the yielded fluid flow is shear thinning. However, many yield-stress fluids deform like elastic solids in the unyielded state and behave as viscoelastic liquids in the yielded state, displaying elastic (E), viscous (V) and plastic (P) properties.
A simple dynamic elastoviscoplastic (EVP) constitutive model that can be integrated with direct numerical simulations was proposed by Saramito (2007). The model has proven to capture viscoplastic and elastic effects (Cheddadi et al. 2011;Fraggedakis, Dimakopoulos & Tsamopoulos 2016) and to properly match experimental results and observations (Holenberg et al. 2012) (common materials used to study this type of EVP fluids are Carbopol solutions and liquid foams Firouznia et al. 2018;Zade et al. 2020). The model was extended by the same author to account for shear-thinning effects (Saramito 2009), combining the Oldroyd viscoelastic model with the Herschel-Bulkley viscoplastic model, with a power-law index that allows a shear-thinning behaviour in the yielded state. When the index is equal to unity, the model reduces to the one proposed in his previous work, i.e. Saramito (2007). Apart from the models proposed by Saramito, many other EVP models exist in the literature. The interested reader is referred to Crochet & Walters (1983), Balmforth, Frigaard & Ovarlez (2014), Saramito (2016) and Saramito & Wachs (2017) for a thorough review of models and numerical methods.
Turbulent flow of purely viscoelastic fluids has gained attention in the drag-reduction and flow control communities, since a tiny amount of polymer (parts per million) has proven efficient in reducing friction drag in pipe flows (Virk 1971). Drag reduction by polymers is related to their ability to modify coherent structures in wall-bounded turbulence (Dubief et al. 2004(Dubief et al. , 2005. Polymers influence the turbulent cycle in two ways: they attenuate near-wall vortices, but at the same time they also increase the streamwise kinetic energy of the near-wall streaks. The net balance of these two opposite actions leads to a self-sustained drag-reduced turbulent flow. More recently, Xi & Graham (2010, 2012a proposed that polymeric drag reduction is a time-dependent process: turbulent flows with polymers spend relatively more time in hibernating phases, where the turbulence is weak, and less time in active phases with strong turbulence. The turbulent coherent structures, streamwise vortices and streaks, were found to differ in appearance depending on whether the flow was in an active or a hibernating phase. Biancofiore, Brandt & Zaki (2017) examined the secondary instability of streaks in viscoelastic flows, showing that the streaks reach a lower average energy with increasing elasticity due to a resistive polymer torque that opposes the streamwise vorticity and, as a result, opposes the lift-up mechanism. Numerous studies have addressed the topic of polymeric drag reduction and all cannot be reviewed here for brevity; the interested reader is referred to White & Mungal (2008) for a thorough introduction to this subject. Turbulent flow of yield-stress fluids (plasticity) has been studied much less, despite its relevance for applications. The reason is twofold. Time-resolved measurements in yield-stress fluids are very challenging, even in laminar flows. Most real-life yield-stress fluids are opaque and hence do not provide optical access. The transparent laboratory fluid Carbopol, often used in these experiments, needs careful preparation to achieve well-controlled properties (Dinkgreve, Fazilati & Bonn 2018;Tabuteau, Coussot & De Bruyn 2018), and exhibits slip on smooth surfaces that needs to be avoided or carefully controlled (Firouznia et al. 2018). On the other hand, direct numerical simulations of yield-stress fluids have not been affordable until recently. Rosti et al. (2018a) first studied the turbulent channel flow of a yield-stress fluid in a near-viscoplastic limit using the Saramito EVP model (Saramito 2007), adding a tiny amount of elasticity to achieve numerical stability (the Weissenberg number was Wi = 0.01). When the Bingham number was gradually increased from zero, the flow became less turbulent and more correlated in the streamwise direction, until it completely relaminarised at Bi = 2.8. The velocity correlations revealed that the size and length of the near-wall streaks increased with the Bingham number. The friction factor decreased with the Bingham number in the turbulent regime, while it increased with the same in the laminar regime. Le Clainche et al. (2020) analysed further the simulation data of Rosti et al. (2018a) using high-order dynamic mode decomposition, and compared the modes with those in Newtonian fluids, and also in purely viscoelastic fluids. Their results indicated that elasticity and plasticity have similar effects on the coherent structures; in both cases, the flow is dominated by long streaks disrupted by rapid, localised perturbations. The Newtonian flow, on the other hand, displays short streaks and a more chaotic dynamics. Very recent experiments in a duct flow of Carbopol confirmed that the energy content at low wavenumbers and streamwise anisotropy were higher than in Newtonian turbulence (Mitishita et al. 2021), indicating that streamwise near-wall structures (streaks) were enhanced in Carbopol. These experiments at higher Reynolds numbers also confirmed the qualitative changes that viscoplasticity causes in mean flow profiles and Reynolds stresses, as observed by Rosti et al. (2018a).
In addition, many earlier direct numerical simulation studies addressed the turbulence in Bingham pseudoplastic fluids (obtained by regularising the Bingham model by a large viscosity) (Rudman & Blackburn 2006;Guang et al. 2011;Zhu et al. 2020), and their shear-thinning version, regularised Herschel-Bulkley fluids by Rudman et al. (2004) and Singh, Rudman & Blackburn (2017). Very recently, Zhu et al. (2020) studied the turbulence of a Bingham pseudoplastic fluid in a vertical channel with particles, at a higher bulk Reynolds number (Re b = 8000) than Rosti et al. (2018a). They also obtained a drag reduction with an increasing Bingham number, along with an increasing streamwise coherence of the flow structures and observed that the turbulent statistics were asymmetric with respect to the centreline at higher Bingham numbers.
Transitional flows of yield-stress fluids in pipes and channels have been addressed in many studies, both experimentally and computationally. An asymmetric mean flow profile is a characteristic feature observed in pipe flow experiments of yield-stress fluids (Escudier et al. 2005;Peixinho et al. 2005;Esmael & Nouar 2008;Guzel, Frigaard & Martinez 2009), and for fluids with a Herschel-Bulkley-type rheology (Escudier & Presti 1996). According to Esmael & Nouar (2008), this was caused by an asymmetric nonlinear coherent structure consisting of two counter-rotating vortices. Experiments on laminar steady flow of a yield-stress fluid in a duct laden with particles were performed recently by Zade et al. (2020). Computational studies on transitional flows have mainly focused on the linear and nonlinear stability of channel flows (Frigaard, Howison & Sobey 1994;Nouar & Frigaard 2001;Frigaard & Nouar 2003;Nouar et al. 2007;Metivier, Nouar & Brancher 2010;Nouar & Bottaro 2010;Bentrad et al. 2017). Most importantly, within the Bingham model, the central plug regions of a channel flow remain unyielded despite linear perturbations, and hence the flow is always linearly stable. In non-modal analysis, in contrast to Newtonian fluids, the optimal disturbance is found to be an oblique wave (Nouar et al. 2007), associated with the lift-up effect (Schmid 2007;Brandt 2014). However, Nouar & Bottaro (2010) found that a slightly perturbed mean flow profile supports exponential amplification of streamwise-travelling waves, indicating another possible transition scenario for the plane channel flow of a yield-stress fluid.
Concluding, all previous studies on EVP turbulence focused on either purely viscoelastic or near-viscoplastic fluids. The questions remain as to what happens when both elasticity and plasticity effects are finite and interact with each other. The recent experimental study by Mitishita et al. (2021) also raised the question of how the findings of Rosti et al. (2018c) would be affected by finite elasticity. The present study focuses on EVP fluids at finite Weissenberg numbers, i.e. fluids with finite viscoelasticity and yield stress.
In this work, we perform direct numerical simulations of turbulent channel flows of an incompressible EVP fluid and extend the one by Rosti et al. (2018c) by considering a wide range of Weissenberg numbers 0 ≤ Wi ≤ 16 and Bingham numbers 0 ≤ Bi ≤ 22.4, all at the bulk Reynolds number Re = 2800. The non-Newtonian flow is simulated by solving the full unsteady incompressible Navier-Stokes equations coupled with a modified version of the model proposed by Saramito (2007) for the evolution of the additional EVP stress tensor.
The manuscript is organised as follows. In § 2, we first discuss the flow configuration and the governing equations, and then present the numerical methodology used. The new model with some test cases is reported in § 2.1, while the new results are presented in § 3. In particular, we discuss the role of elasticity and plasticity on a turbulent channel flow. Finally, a summary of the main findings and conclusions are collected in § 5.

Formulation
We consider the laminar and turbulent flows of an incompressible EVP fluid through a plane channel with two impermeable rigid walls. Figure 1(a) shows a sketch of the geometry and the Cartesian coordinate system, where x, y and z (x 1 , x 2 and x 3 ) denote the streamwise, wall-normal and spanwise coordinates, while u, v and w (u 1 , u 2 and u 3 ) denote the respective components of the velocity field. The lower and upper stationary impermeable walls are located at y = 0 and 2h, respectively, where h represents the channel half-height.
The fluid motion is governed by the conservation of momentum and the incompressibility constraint where ρ is the fluid density and σ ij the total Cauchy stress tensor, which is written as σ ij = −pδ ij + 2μ f D ij + τ ij , where p is the pressure, μ f the fluid molecular dynamic viscosity (also called solvent viscosity), δ the Kronecker delta and D ij the strain rate tensor defined as D ij = (∂u i /∂x j + ∂u j /∂x i )/2. The additional EVP stress tensor τ ij accounts for the non-Newtonian behaviour of the fluid, here described by a modified version of the model proposed by Saramito (2007). In the original model, when the stress σ is below the yield stress τ 0 , the system predicts only recoverable Kelvin-Voigt viscoelastic deformations, while when the stress exceeds the yield value τ 0 , the fluid behaves as an Oldroyd-B viscoelastic fluid. Thus, the total strain rateε is shared between an elastic contributionε e and a plastic oneε p (Cheddadi et al. 2011). Since the Oldroyd-B model assumes infinitely stretched dumbbells, the range of application of the model is limited to low elasticity. To overcome this drawback, we instead use the finite extensible nonlinear elastic-Peterlin (FENE-P) model. Thus, the instantaneous values of all the components of the stress tensor τ ij are found by solving the following transport equation: where Here, λ is the relaxation time, μ m is an additional viscosity, L is the maximum polymer extensibility, τ 0 the yield stress and |τ d | represents the second invariant of the deviatoric part of the added stress tensor. The EVP parameters μ f , μ m , λ and τ 0 can be obtained by experimental data following the procedure detailed by Fraggedakis et al. (2016), based on the determination of the linear material functions, i.e. the storage modulus G and the loss modulus G . The above constitutive equation gives a Kelvin-Voigt solid in the limit F ≈ 1, which is the case in the unyielded state, and a behaviour like a FENE-P fluid in the yielded state. The previous equation can be rewritten in terms of the conformation tensor where the EVP stress tensor τ ij is related to the conformation tensor by the relation We use the log conformation method to overcome the well-known high Weissenberg number problem (Fattal & Kupferman 2004;Hulsen, Fattal & Kupferman 2005;Izbassarov & Muradoglu 2015;De Vita et al. 2018). In this approach, (2.4) is rewritten in terms of the logarithm of the conformation tensor through an eigen-decomposition, i.e. Ψ = log C, which ensures the positive definiteness of the conformation tensor. The core feature of the formulation is the decomposition of the gradient of the divergence free velocity field ∂u j /∂x i into two anti-symmetric tensors, Ω ij (pure rotation) and N ij , and into a symmetric tensor, B ij , which commutes with the conformation tensor, i.e.
The modified transport equation to be solved is thus Once this equation is solved, the conformation tensor can be obtained using the inverse transformation C = e Ψ . Finally, the full system of equations can be rewritten in a non-dimensional form as where we have used the same symbols to define the non-dimensional variables for simplicity.  Shahmardi et al. 2019). The governing equations are discretised with the second-order centred finite difference scheme on a staggered uniform grid, except for the advection term in (2.7) where the fifth-order WENO (weighted essentially non-oscillatory) scheme is adopted (Shu 2009;Sugiyama et al. 2011). The time integration is performed with a fractional-step method (Kim & Moin 1985) and a third-order explicit Runge-Kutta scheme except for the EVP stress term which is advanced with the Crank-Nicolson method.
In all the simulations, periodic boundary conditions are used in the streamwise and spanwise directions, while the no-slip and no-penetration boundary conditions are enforced on the solid walls. For all the turbulent cases considered hereafter with a bulk Reynolds number equal to Re b = 2800, the equations of motion are discretised by using 1728 × 576 × 864 grid points on a computational domain of size 6h × 2h × 3h in the streamwise, wall-normal and spanwise directions, with the resolution satisfying the constraint x + = y + = z + < 0.6, where the superscript + indicates the wall units. The spatial resolution has been chosen equal to the one used in a previous work (Rosti et al. 2018a) in order to properly resolve the turbulent scales as well as the unyielded plug regions which form intermittently in the domain, and verified a posteriori with a grid refinement study. In the low Reynolds fully laminar cases, the spatial resolution was relaxed and the domain size in the homogeneous directions reduced. Finally, all the turbulent cases are initialised with a fully developed channel flow with zero EVP stress (τ ij = 0); after the flow and stresses have reached a statistically steady state, the calculations are continued for an interval of approximately 500 bulk time units, during which the statistics are computed.

Validations of the FENE-P-based EVP model
The present implementation for single and multiphase flows of an EVP fluid has been extensively validated in the past; here, we report further test cases of the new FENE-P-based EVP model also to present its rheological behaviour.
First, we consider a three-dimensional uniaxial elongational flow, with a constant elongational rate˙ 0 : the Weissenberg number Wi = λ˙ 0 and extensibility parameter L 2 are varied in the range of 0.25 ≤ Wi ≤ 0.75 and 10 ≤ L 2 ≤ ∞, respectively, while the Bingham number Bi = τ 0 /(μ 0˙ 0 ) = 1 and the viscosity ratio β = 0. The fluid flow is assumed to have a constant velocity gradient ∂u i /∂x j = diag{1, −1/2, 1/2} and (2.2) can be simplified as with τ ii (0) = 0 for i = 1, 2 and 3. The time evolution of τ 11 − τ 22 (the normal stress difference) is reported in figure 2; we observe that, initially, the stress components grow linearly, but when the stress level is above a threshold, i.e. the yield stress, the Saramito model (L 2 = ∞) predicts an unbounded growth for Wi > 0.5, while the FENE-P-based model exhibits a limited growth with the stress difference reaching a plateau. This example easily explains the reason why the FENE-P-based model is chosen for the current work, where we aim to investigate the effect of finite elasticity. Next, we consider a simple constant shear flow, with shear rateγ 0 ; the Weissenberg number is varied Wi = λγ 0 in the range 0.01 ≤ Wi ≤ 100, while the Bingham number Bi = τ 0 /(μ 0γ ) assumes 4 distinct values and the viscosity ratio β = 0. We consider two cases, with L 2 = 100 and L 2 = ∞ to study the properties of the two models. For this two-dimensional problem, we have ∂u i /∂x j = 0 1 0 0 and (2.2) can be rewritten as The dimensionless steady shear viscosity μ s /μ 0 = β + τ 12 /γ as a function of Wi is reported in figure 3. As expected, both models exhibit a shear-thinning behaviour: in particular, for L 2 = ∞ the model tends to a plateau at both low (Wi 0.1) and high (Wi 10) Weissenberg numbers, while for L 2 = 100 the steady shear viscosity reaches a plateau at low Wi but it decreases monotonically at high Wi. Note that the shear-thinning effect is controlled by Bi, which modifies the plateau value at low Wi.

Laminar flow
Our first analysis considers the laminar flow of an EVP fluid. All laminar cases are fixed at β = 0.9, Re = 2800 and L 2 = 3600, while the Weissenberg and the Bingham numbers are varied in the range of 0.01 ≤ Wi ≤ 16 and 0.1 ≤ Bi ≤ 100, respectively. Note that the domain size in the homogeneous directions was kept very small to avoid the development of a turbulent flow and that no perturbations were added to the zero initial condition.
First, we consider the combined effects of the Bingham and the Weissenberg numbers on the frictional resistance of the flow quantified by the Fanning friction factor f , defined as 2τ w /ρU 2 b with τ w the total wall shear stress including both the viscous and EVP contributions. Figure 4(a) shows the Fanning friction factor f as a function of the Weissenberg number for various Bingham numbers. In particular, the cyan, brown, purple  and green lines are used for Bi = 0.1, 2.8, 11.2 and 100, respectively. The friction factor f decreases nearly linearly with Wi (in log scale). However, at Bi = 0.1 the factor f is constant, which is consistent with the viscoelastic flows in the laminar regime. Results clearly show that the elastic effects become more prominent with Bi, leading to a steeper decay with Wi. Next, we consider the reverse case, where we study the effects of Bingham number at a fixed Wi, as shown in figure 4(b). We find that f increases nonlinearly with the Bingham number Bi, and that the dependency on Bi decreases with the Weissenberg number Wi.
The increase of f can be explained by subtle changes in the mean streamwise velocity profile. The non-dimensional streamwise velocity profile u/U b is shown in figure 5(a) at Wi = 0.01 and Bi = 0.1, 2.8, 11.2, 100 as a function of the wall-normal distance y/h. It can be seen that the solid plug in the middle of the channel increases with the Bingham number. The plug moves with a uniform velocity, and its velocity decreases with Bi due to mass conservation, leading to an increase in the wall shear stress.
The variations of the non-dimensional mean velocity profile u/U b at Bi = 22.4 with the Weissenberg number, see figure 5(b), demonstrate that the solid plug in the middle of the channel decreases with Wi. Indeed, the total stress magnitude in a channel flow increases with elasticity, leading to earlier yielding (Chaparian & Tammisola 2019). To further quantify these observations, the volume of the solid region is shown in figure 6 for various combinations of Wi and Bi. The results are consistent with the observations in figure 4, since the solid volume Vol s is proportional to the friction factor f .

Turbulent flow
We examine turbulent channel flows of an EVP fluid, together with the baseline Newtonian (Bi = 0 and Wi ≈ 0), viscoelastic (Bi = 0) and almost viscoplastic (Wi ≈ 0) cases. All simulations have been performed until statistically steady state (at least 500 non-dimensional units). Hence, in the turbulent flow cases, the turbulence is always fully developed. However, for some parameter values the flow is not turbulent, because it has laminarised by increasing elasticity or plasticity.
The flow rate is constant in all simulations, so that the flow Reynolds number based on the bulk velocity is fixed, i.e. Re = ρU b h/μ 0 = 2800, where the bulk velocity U b is the average value of the mean velocity computed across the whole domain and μ 0 is the total zero-shear viscosity. In the turbulent regime, this bulk Reynolds number corresponds to a nominal friction Reynolds number Re τ = ρu τ h/μ 0 = 180 for a Newtonian fluid, u τ being the friction velocity. A wide range of elasticity and plasticity is investigated: the Weissenberg number Wi = λU b /h is varied in the range 0 ≤ Wi ≤ 16 to study the role of fluid elasticity, while the Bingham number Bi = τ 0 h/μ 0 U b is varied between 0 and 22.4. The viscosity ratio and the extensibility parameter are fixed to β = μ f /μ 0 = 0.9 and L 2 = 3600.
The viscous units used above are indicated by the superscript + , and are built using the friction velocity u τ as the velocity scale and the viscous length δ ν = ν/u τ as the length scale. Here, we define the friction velocity as where we have taken into account also the EVP shear stress at the wall. Note that the actual value of the friction velocity in our simulations is computed from the friction coefficient, obtained with the driving streamwise pressure gradient, rather than from its definition. We first discuss the qualitative flow characteristics for various combinations of Wi and Bi. A qualitative overview of the resulting flow regimes as functions of Wi and Bi is presented in figure 7.
The empty and filled symbols in figure 7 represent turbulent and laminar regimes, respectively. As expected, the flow becomes laminar and steady with an increasing Bingham number. The critical Bi c can be defined as the lowest value of Bi where laminarisation occurs, at each Wi. At Wi = 0.01, the flow is near viscoplastic and fully laminarises at the relatively low critical Bi c = 2.8 (Rosti et al. 2018c). (Note that results for viscoplastic cases used in this work are taken from our previous work Rosti et al. (2018c).) Laminarisation of the viscoplastic flow with increasing Bingham number is closely related to the increase of the core unyielded region which grows from the centreline towards to the walls, and damps out turbulent fluctuations in the core. As described in Rosti et al. (2018c), the near-wall streaks were intensified and became more coherent. The low-speed streaks, usually associated with positive wall-normal fluctuations, reach higher wall-normal distances than the high-speed streaks, thus inducing the flow to yield at higher wall-normal distances if the local stress reaches the yield-stress threshold. Indeed, the unyielded regions preferentially form above high-speed streaks. Overall, the flow becomes more and more correlated in the streamwise direction when increasing the Bingham number, with high levels of flow anisotropy close to the wall, similarly to what observed in other drag-reducing flows. Differently from the other flows, however, both the streamwise and spanwise correlations grow with the Bingham number also away from the wall, due to the growth of the unyielded region.
As the flow becomes more elastic, the critical Bi c shows a non-monotonic behaviour; it first increases with Wi, then decreases. The reason for this non-monotonic behaviour is attributed to the fact that the coherent structures are highly influenced by the relative width of the plug compared with the yielded regions around it, as will be described later.
The colours in figure 7 represent a measure of drag reduction compared with a Newtonian fluid. Panel (a) shows DR = [1 − (Re τ /Re τ,Wi=0 ) 2 ] × 100 %. For flows in the turbulent regime, DR increases with Bi for all Wi studied in this work. In the laminar regime, however, at low Wi = 0.01, the DR changes sign, showing drag increase, whereas for higher Wi ≥ 2 values, the DR only slightly decreases with Bi. This is consistent with the results for the Fanning factor shown in figure 8. It is well known that a finite viscoelasticity reduces drag compared with a Newtonian flow, and this has been exploited in many previous studies and in the industry. The figure seems to show that a combination of finite Wi and Bi results in a much greater drag reduction than a finite Wi alone. However, a closer inspection shows that the very high drag reduction DR > 70 % occurs when the flow laminarises, because we compare with the drag in the Newtonian flow in the turbulent regime. We wanted to eliminate the effect of laminarisation on drag reduction, and also try to exclude the effects of shear thinning. Therefore, in (b), we show another measure of drag reduction which excludes the influence of shear thinning (Housiadas & Beris 2003): where f EVP and f N denote skin friction for EVP and Newtonian cases, respectively. Following Dean (1978)  where Re w = ρU2h/μ w and μ w = τ w /γ w , and subscript w denotes wall quantities. Using DR 2 as a measure, the highest drag reduction is achieved either at high Wi, or in EVP flow when Wi and Bi are both moderate, in the same region where laminarisation is delayed to higher Bi. In general, the drag-reduction values remain very similar for these two measures in the turbulent regime, indicating that shear thinning has no major influence on drag reduction at the studied parameter values. Despite this, the viscosity profiles (figure 8) show shear thinning for VE, EVP and VP cases. Apparently, plasticity plays a strong role on shear thinning, as both VP and EVP show stronger shear thinning. For the EVP case at Wi = 4 and Bi = 5.6, we are approaching results of VP at Bi = 1.4, thus elasticity seems to decrease shear thinning. The main area responsible for the non-monotonic trends is the unyielded region of the flow. The unyielded regions show a non-trivial and complex trend with Bi and Wi, with significant differences also between the laminar and turbulent regimes. In order to better understand the flow, we start by showing visualisations of the instantaneous distributions of the yielded and unyielded regions in figure 9. In particular, the figure shows instantaneous cross-stream planes (x-y plane) where the unyielded regions are coloured in brown; in the yielded regions we report colour contours of the spanwise vorticity, ω z . At Bi = 0 and Wi ≈ 0 we recognise the classic vorticity field of turbulent channel flows, with high vorticity levels at the walls, and the footprints of the classical turbulent streaky structures. As Wi increases, the flow becomes smoother, the flow structures seem to be more elongated and the streamwise coherency of the flow increases, as is typical of drag-reducing flows. Note that, although the turbulence activity is reducing with Wi, the purely viscoelastic flow does not laminarise. On the other hand, for every Wi, the flow returns to laminar for a sufficiently high Bingham number (Bi > Bi c ). In particular, we observe that as Bi increases, the amount of fluid which is unyielded grows, eventually forming a fully connected plug spanning the whole streamwise and spanwise lengths, thus leading to the decay of turbulence and the return to a fully laminar regime. However, at intermediate Wi, we observe that a higher Bi c is required to laminarise the flow than at low or high Wi. At intermediate Wi, the equilibrium solutions indicate that the unyielded region is narrower than in the viscoplastic flow (Chaparian & Tammisola 2019), and we propose that this allows instabilities to be sustained in the yielded regions. In the present simulation at intermediate Wi, a plug was formed initially and the flow seemed to laminarise, but soon a large-scale roll-up motion developed at the edges of the unyielded region, recreating the unsteadiness and turbulence. On the other hand, at higher values (Wi > 4) this effect does not occur thus leading to lower Bi c . This phenomenon could be attributed to the fact that the unyielded region becomes too thin to create any instability. Moreover, the yielded region becomes smoother at high Wi enhancing laminarisation.
It is not obvious what drives the instability leading to the roll-up motion, because the laminar flow profile is not inflectional. A sharp viscosity contrast can lead instabilities in channel and shear flows. The unyielded region could play a similar role as a high-viscosity fluid in the centre of the channel, a configuration which can be unstable at high Reynolds numbers (Hooper & Boyd 1987;Govindarajan & Sahu 2013). On the other hand, a channel flow of Bingham fluid is known to be linearly stable at all Reynolds numbers. Although the linear stability of EVP flow has not been studied, the transition that we observe could be subcritical, like for Bingham fluids. However, the roll-up structure we observe does not resemble the optimal non-modal instabilities in Bingham fluid: streamwise-travelling waves (Nouar et al. 2007) or oblique waves (Nouar & Bottaro 2010).
For laminar flows (Bi > Bi c ), a further increase of Bi leads to a growth of the unyielded region, which induces an increase of the wall shear stress and a consequent drag increase,

Flow statistics
Here, we discuss the turbulent statistics for EVP flows. To facilitate comparisons, we provide results for the two limiting cases: near-viscoplastic and viscoelastic flows, and then compare these with the flows with finite elasticity and plasticity. The near-viscoplastic flow was presented in our earlier work (Rosti et al. 2018c), and was computed using the Saramito model at Wi = 0.01, β = 0.95, and varying Bi. The viscoelastic flow is computed at varied Weissenberg number 0 ≤ Wi ≤ 16 and fixed Bingham number Bi = 0, and the FENE-P model parameters are β = 0.9 and L 2 = 3600, unless indicated otherwise.
We first discuss the mean streamwise velocity profiles in relation to drag reduction shown in figure 7. In figure 10, the velocity profiles are shown in wall units. The profiles for viscoplastic flow in figure 10(a) are similar to the Newtonian case for Bi < 1.4. For Bi = 1.4, the difference becomes more significant with zero-shear region at the centreline. At Bi = 2.8, the flow becomes laminar with a large zero-shear region, leading to drag increase as observed in figure 7. We can observe that the flow remains unyielded mostly in the logarithmic and outer layer, while it is always yielded in the viscous sublayer.
For viscoelastic flows, DR% increases monotonically with Wi in the range of Wi = O(1) and start to saturate at Wi > 10. The maximum drag reduction is close to 40 % at Wi = 16. Figure 10(b) shows the mean streamwise velocity profiles in wall units. It can be seen that increasing Wi the velocity profiles collapse in the viscous sublayer, however, shift upward in the buffer region (thickening of the buffer layer) and are parallel to the Newtonian profile. Note that we do not reach Virk's maximum drag-reduction (MDR) asymptote (Virk 1975), which is out of scope of the present study. The shift in the velocity profile is consistent with results in the literature (Xi & Graham 2010;Shahmardi et al. 2019).
We next consider the EVP flow. The drag-reduction values and the corresponding mean streamwise velocity profiles for Wi = 4, β = 0.9, L 2 = 3600 and Bi = 0, 2.8, 5.6, 11.2 and 22.4 are shown in figure 11. Similarly to the viscoplastic flow, the DR% increases with Bi in the turbulent regime. Unlike pure viscoelastic case, plasticity enhances the DR going from low-drag reduction (LDR, DR 40) at 0 ≤ Bi < 5.6 to high drag-reduction (HDR, 40 < DR 60) at Bi = 5.6. Moreover, the laminar regime delays to larger Bi > 5.6, which may be attributed to the dramatic effects due to the polymer dynamics. Unlike the viscoplastic case, DR% only slightly decreases in the laminar case as shown also earlier in figure 7. The change of mean streamwise velocity profile with Bi is more significant than in figure 10(a) and more similar to figure 10(b) for Bi ≤ 5.6, i.e. the velocity profile collapse in the viscous sublayer, and shift upward in the buffer and log-law layer approaching Virk's asymptote. In contrast, for Bi > 5.6 the flow fully laminarises and becomes similar to the case in figure 10(a), with the zero-shear region at the centreline. Next, we analyse the wall-normal distribution of the Reynolds stress components u u and v v , shown in figure 12 together with the Newtonian case by Kim et al. (1987) (blue symbol). The profiles are normalised with u 2 τ . The results for the viscoelastic case (a,b) show that the u u increases while v v decreases with Wi. In general, the shape of the profiles is similar to the corresponding Newtonian one, but the magnitude changes with Wi. Moreover, all the peaks move towards the channel centre, in agreement with the thickening of the buffer layer observed in the mean velocity profile. A similar behaviour has been noticed in the literature for drag-reducing turbulent viscoelastic flows (White & Mungal 2008;Dallas, Vassilicos & Hewitt 2010). The VP case (Wi = 0.01) is depicted in (c,d) to examine the effects of the Bingham number in the range 0 ≤ Bi ≤ 2.8. We clearly observe that the peak of the streamwise component u u increases with Bi, while it decreases for v v . A similar trend is observed for the EVP case (Wi = 4) as depicted in (e,f ).
Although the overall trends are similar, the combined effects of elasticity and plasticity in the EVP flows result in larger deviations from the Newtonian flow than for the two limiting cases, particularly for the peak in the u u profile. It is also interesting to note that in the EVP flow the peak moves non-monotonically, unlike in the VE and VP cases: it first moves away from the wall as Bi increases (Bi ≤ 5.6), and then moves slightly back towards the wall. The stress components u v and w w show similar trends as v v , and are shown in figure 13 for completeness. Figure 14 shows P s , the probability of the fluid being unyielded in a given vertical position, and Vol s , the percentage of the volume that was unyielded on average, for VP (Wi = 0.01) and EVP (Wi = 4 and 8) cases. The probability P s is displayed against the  Kim et al. (1987). wall-normal distance y, whereas the volume percentage is depicted in the legend. The probability P s ranges from 0 (always fluid) to 1 (always elastic solid).
For the laminar case, as expected, P s sharply changes along the interface between the unyielded and yielded region. For turbulent cases, the probability P s increases smoothly across the channel with a maximum in the vicinity of the centreline and a minimum near the wall. This happens due to the time-dependent nature of the turbulent flow.
For the VP case, P s and Vol s increase with Bi, reaching the maximum (Vol s = 54 %) at Bi = 2.8, when the flow is fully laminar. The volume of the unyielded region further  Kim et al. (1987). increases when Bi increases in the laminar regime, as illustrated in figure 5(a). Let us now consider the EVP cases. A careful inspection of figure 14 reveals that the unyielded volume Vol s is smaller in EVP flows than in VP flows, in both laminar and turbulent states. Also, Vol s further decreases as Wi increases, from 30 % at Wi = 4 to 26 % at Wi = 8, as seen in figure 14(b,c). The reason is that the elastic stress in channel flows increases with Wi, resulting in a smaller unyielded volume. Moreover, the Vol s increases monotonically with Bi also for the EVP case, as long as Bi < Bi c . However, when the flow laminarises (Bi ≥ Bi c ), Vol s decreases owing to a narrow plug region. An explanation to the decrease can be found by comparing the laminar and turbulent mean flow profiles. In laminar VE  channel flows, the elastic stress is known to be proportional to both the local shear and Wi, and therefore, the elastic stress forms the largest part of the total stress at moderate Wi. In laminar flows, the shear is finite almost everywhere across the channel. In turbulent flows, however, the mean flow shear is small in the central region, and hence the elastic stresses are smaller on average, with the Reynolds stresses becoming more important in yielding the material. We now present the statistics of the polymer conformation tensor for VE, VP and EVP cases. The mean profiles of the polymer stretching Tr(C)/L are shown in figure 15(a) for the VE case at Wi = 2, 4, 8 and 16, where the dashed line in the figure indicates the case of coiled polymers (conformation tensor equal to identity). As expected, the polymer stretching increases monotonically with the Weissenberg number Wi. The mean polymer stretching first increases with maximum in the vicinity of the wall, then decreases with minimum at the centre. The observed near-wall peak is mostly due to the interaction between polymers and near-wall vortices (Xi & Graham 2010;Dubief, Terrapon & Soria 2013). Viscoelastic effects become more prominent throughout the domain when increasing Wi and the peak values move away from the wall.
The mean polymer extension for the VP cases at 0 ≤ Bi ≤ 2.8 is shown in figure 15(b). For Bi 2, corresponding to the turbulent regime, the stress profile is nearly constant except in the near-wall region where it is maximum. Note that, for the laminar case, the stress monotonically decreases along the wall-normal direction. Finally, the results pertaining the EVP flows at 0 ≤ Bi ≤ 11.2 and Wi =4 & 8 are presented in figure 15(c,d). For the turbulent EVP cases (Bi < 11.2), the stress profiles are qualitatively similar to the VE ones. The maximum values are in the vicinity of the wall and the minimum at the centreline (y = h). Moreover, the minimum values increase with Bi while peaks are not sensitive to Bi at Bi < Bi c . For the laminar flows, the profile is qualitatively similar to the viscoplastic one, with maximum and minimum values at the wall and centreline, respectively.

Energy and stress balance
To further characterise the flow, we study the energy and shear stress budget for the viscoelastic and EVP cases. We first consider how the energy stress budgets evolve with Weissenberg number, at zero Bingham number.
An overall view of the velocity fluctuations is obtained from the turbulent kinetic energy K = (u 2 + v 2 + w 2 )/2, normalised with u 2 τ (figure 16a). Following Dallas et al. (2010), the budget of turbulent kinetic energy integrated over the whole domain can be written as follows: y + Wi Figure 16. VE case: wall-normal profiles of (a) the turbulent kinetic energy K = (u 2 + v 2 + w 2 )/2 and of (b) the turbulent production P = −ρu v dū/dy, (c) the turbulent dissipation ε = μ∂u i /∂x j ∂u i /∂x j and (d) average contributions to the energy balance from P, ε and dissipation due to the polymers Π = ∂u i /∂x j τ ij , all normalised with the friction velocity u 3 τ . The symbols in the figure correspond to results of Kim et al. (1987) for the Newtonian case.
where V is the volume of the domain, P = −ρu v dū/dy is the turbulent production, ε = μ∂u i /∂x j ∂u i /∂x j is the turbulent dissipation and Π = ∂u i /∂x j τ ij is the dissipation due to the polymers.
The turbulent kinetic energy increases with Wi, and its profile is similar to the streamwise velocity fluctuation in figure 12(a), as the streamwise velocity is of larger magnitude than the other velocity components. On the other hand, the turbulent production by the Reynolds shear stress P = −u v dū/dy is continuously reduced when increasing Wi. The peak values occur in the buffer layer, moving towards the channel centre with increased DR consisted with thickening of the buffer layer. We also computed the turbulent dissipation ε = μ∂u i /∂x j ∂u i /∂x j of the fluctuating velocity field u i . The dissipation ε also attenuates with Wi, while the profile maintains its shape with maximum value at the wall and minimum at the centreline. Volume-averaged quantities of P dy, ε dy and dissipation due to polymers Π dy are shown in figure 16(d). It can be seen that the average turbulent production and viscous dissipation decrease with Wi, while the polymeric dissipation increases with Wi. The decrease in P and ε with Wi is consistent with the observed drag reduction. Note that Π first increases with Wi, but slightly Next, we show the mean viscous (τ V ) and viscoelastic stress (τ E ) profiles as a function of the Weissenberg number, see figure 17. The viscous stress profiles coincide close to the wall and in the vicinity of the channel core but increase slightly in the buffer layer. The viscoelastic stress profile has a peak in the vicinity of the wall and minimum in the centre. Theτ E changes non-monotonically, it first increases reaching maximum in the buffer layer and then decreases reaching a minimum at the centreline. With increase in Wi the peak moves towards the core region.
To gain further understanding, the shear stress budget is shown in figure 18 for the VE cases with Wi = 2 and Wi = 16, normalised with the corresponding wall stress. The total shear stressτ T can be written asτ T =τ V +τ E +τ R = (1 − y/h)τ w where the viscous stressτ V = μ dū/dy, the Reynolds stressτ R = −ρu v and viscoelastic stressτ E tensor. Note that, unlike the Newtonian case, the total wall shear stress is sum of the viscous and viscoelastic contributions. For Wi = 2, the behaviour is similar to the Newtonian case with the viscous stress dominating close to the wall and the Reynolds stress in the core. The viscoelastic effects become more pronounced at Wi = 16, in that the viscous and viscoelastic stresses increase, while the Reynolds stresses diminish.
Finally, we perform the same energy and momentum budget analysis to consider the combined effects of elasticity and plasticity. The energy budget is shown in figure 19 at fixed Weissenberg number Wi = 4, while varying the Bingham number in the range 0 ≤ Bi ≤ 11.2. Qualitatively, the effect of Bingham number on the energy budget is similar to the viscoplastic case in Rosti et al. (2018a), the peak of K increases with Bi, while the peak values for P and ε decrease. The volume-averaged viscous dissipation decreases, while EVP dissipation increases with Bi. Surprisingly, however, the turbulent production slightly increases with Bi at this finite Weissenberg number, in contrast with the viscoplastic flow. The shear components of the mean viscous and EVP stress tensors are shown in figure 20. In viscoplastic flows, we observed that the shear component of the EVP stress tensor increased everywhere across the channel. In the combined case, however, (figure 20b) its value at the wall remains constant in the turbulent regime. In the laminar . EVP case: wall-normal profiles of (a) the turbulent kinetic energy K = (u 2 + v 2 + w 2 )/2 and of (b) the turbulent production P = −ρu v dū/dy, (c) the turbulent dissipation ε = μ∂u i /∂x j ∂u i /∂x j and (d) average contributions to the energy balance due to P, ε and dissipation due to fluctuation of the EVP stresses Π = ∂u i /∂x j τ ij , all normalised with the friction velocity u τ . The Weissenberg number Wi and the extensibility parameter L 2 are fixed and equal to 4 and 3600, respectively (FENE-P model). The symbols in the figure correspond to results of Kim et al. (1987) for the Newtonian case.  Figure 21. EVP case: normalised shear stress balance across the channel, for (a) Bi = 0 and (b) Bi = 5.6. The total shear stress varies linearly across the channel height. The Weissenberg number Wi and the extensibility parameter L 2 are fixed and equal to 4 and 3600, respectively (FENE-P model).
regime, as expected, the shear component changes linearly in the liquid region and is zero in the solid region. The EVP stress continuously increases with Bi towards the centre of the channel. Finally, we report in figure 21 the shear stress budget for the cases with Bi = 0 and Bi = 5.6, normalised with the corresponding wall stress. With increasing Bi, the EVP and viscous stresses grow more across the channel, further reducing the Reynolds stresses, and hence the combined case deviates even more from the Newtonian flow than the viscoelastic case does.
3.5. Intermittent dynamics Finally, we compare the intermittent dynamics for purely viscoelastic, viscoplastic and EVP cases. Previously, Xi & Graham (2010, 2012b observed two distinct time intervals, denoted as hibernating and active, for viscoelastic turbulent flow in minimal channel geometry. During the hibernating intervals, drag reduction is high and turbulence is partly attenuated, whereas the active intervals represent active turbulence similar to the Newtonian turbulent flow.   figure 22(b). Note that velocity profiles are shown in inner units based on the instantaneous friction velocity. For convenience, data for the Newtonian case (Kim et al. 1987), represented by blue symbols, are also shown. As expected from previous work mentioned above, the profile changes between the Newtonian and the Virk profile, approaching the latter one at the lowest Re τ (instant v).
The time evolution of the friction Reynolds number for the viscoplastic flow is shown in figure 22(c). Viscoplastic flow does not hibernate, so we choose four time instants in the turbulent regime. A more limited range of friction Reynolds numbers is observed, i.e. 165 < Re τ < 185. Indeed, the instantaneous velocity profile figure 22(d) at the chosen times only slightly deviates from the Newtonian flow, represented by blue symbols in the figure.
Let us now consider the EVP flow where both elastic and plastic effects are significant. It is remarkable to see that the EVP flow (figure 22e) behaves like the VE flow, with significant drag reduction in the hibernation regime. Interestingly, the elastic characteristics of the intermittent turbulent dynamics are even stronger for the EVP flow than for the VE flow as also observed in the velocity profiles. The profiles near the active regime (instants i and v) are similar to LDR profiles, whereas the ones under hibernation are more like a HDR, with instant (iv) being very close to the Virk profile. Correspondingly, since the EVP flow hibernates longer periods, we observe a more significant drag reduction (27 % < DR < 58 %) than for the purely viscoelastic flow. It is worth noting that we observed and compared several hibernation cycles in VE and EVP flows, in addition to the ones shown in the figure.
We also examine the instantaneous profiles of the mean polymer extension, wall-normal and shear components of the conformation tensor for VE and EVP flow in figure 23. The polymer extension is maximal in the active regime (black and red lines), while it is minimal in the hibernation phase (orange and green lines). The same was observed in viscoelastic turbulent flow in a minimal channel (Xi & Graham 2010, 2012b. While the qualitative picture is the same for VE and EVP, they do present some differences. It appears that at hibernation, the EVP flow has larger polymer extension and higher values of the wall-normal component than the VE flow, while the shear component of the EVP flow is lower, which explains drag reduction. On the other hand, the VE flow presents higher values of both components in the active turbulence regime. (We recognise that these profiles are only snapshots and may not represent accurately the whole dynamics.) EVP stresses are lowest in the hibernation state, and consequently, the probability P s of the fluid to be unyielded is highest in the hibernation state, as shown in figure 24. The P s -curves neatly divide into two distinct groups in active vs hibernating state, as characterised by Re τ , and also by the wall-normal component of the confirmation tensor C 22 .

Flow structures
The EVP character of the flow affects the near-wall turbulent structures, and this is visually confirmed in figure 25. In the figure we identify the low-(blue) and high-speed (red) near-wall streaks with isosurfaces of the streamwise velocity fluctuations u corresponding to the levels u + = ±0.25U b for the values of Bi and Wi pertaining all the turbulent cases we have studied. In the Newtonian case (Bi = 0 and Wi = 0), we recognise in the shape of the streamwise velocity fluctuations the typical footprints of the near-wall streaks and quasi-streamwise vortices that characterise the near-wall turbulence and that are responsible for the wall cycle that can sustain turbulence in classical wall-bounded flows (Jiménez & Pinelli 1999). In VE flows (the leftmost column in figure 25), as the elasticity increases, the low-and high-speed streaks become less fragmented, are correlated over longer distances in the streamwise direction and grow in size. This behaviour has been recognised in the past by several authors (Warholic et al. 1999;Escudier, Nickson & Poole 2009;Shahmardi et al. 2019;Le Clainche et al. 2020) who consistently reported that the inclusion of the polymers in the flow strongly modifies the near-wall structures, by overall reduction of the turbulent dissipation previously observed. This results in a reduction of the friction Reynolds number Re τ , i.e. drag reduction, similarly to what observed in other drag-reducing flows, such as riblets and anisotropic porous walls (Choi, Moin & Kim 1993;Rosti, Brandt & Pinelli 2018b). Although the modifications of the flow structures with Wi and Bi are similar, some differences are present. Indeed, plasticity alone in VP flows reduces the wall-normal fluctuations less than elasticity because the hibernating turbulent phases are reduced; this results in stronger ejection events and in turbulent structures that penetrate deeper in the bulk of the channel for high values of Bi. This difference is related to the presence of the unyielded region in the plastic cases, which is inherently intermittent and induces larger values of EVP dissipation rates Π than in the viscoelastic cases (figures 16d and 19d). Because of this, the drag reduction due to the plasticity is able to laminarise the flow while the one of elasticity is not. The two effects mix in a non-trivial way when both finite elasticity and plasticity are considered. In the EVP cases, the flow structures continue to become less fragmented and more streamwise coherent, with their wall-normal size depending on the relative magnitudes of Wi and Bi. Indeed, as also shown in figure 9, the unyielded region shows a complex trend with Wi and Bi, and this might influence the size of the flow structures. The flow always laminarises at high enough Bingham numbers. However, the VP flow is laminar at Bi > 2.8, showing purely plastic turbulence, while the EVP flow at Wi = 2 becomes laminar only at Bi > 5.6, and hence the turbulence shows combined elastic and plastic characteristics. At high Wi (Wi = 16), the trend is reversed and the flow laminarises again at Bi = 2.8, and the effect of plasticity on the turbulent viscoelastic flow is not large due to an early relaminarisation. This non-monotonic trend can be observed in figure 7. The different effect of elasticity and plasticity on the turbulent flow can be better appreciated in figure 26 where we show wall-parallel contours of the instantaneous polymer extension Tr(C)/L at y ≈ 0.15h, together with the footprints of the low-and high-speed streamwise velocity streaks on the plane. In the top row of the figure we consider the viscoelastic cases with Bi = 0 and growing Weissenberg number Wi from 2 to 16. The figure confirms that as Wi grows, the polymer extension is enhanced, resulting in more elongated and streamwise coherent structures. Also, it is evident from the figure that the flow structures grow in size, especially the low-speed ones. The second row depicts the results for a fixed Weissenberg number (Wi = 2) and growing Bingham number from 0 to 5.6. Again, we note that as Bi increases in the EVP flow, the flow becomes more ordered and streamwise coherent, however, from the figure we can also appreciate the different effect of elasticity and plasticity on the EVP flow: plasticity reduces the turbulent dissipation and the fragmentation of the turbulent structures, acting on both the low-and high-speed streaks equally. Indeed, both appear more elongated in the streamwise direction as Bi increases, with only a small growth in their spanwise dimension. On the other hand, elasticity affects mostly the low-speed streaks, while the high-speed ones remain fragmented despite growing in size. A common feature between the two flows is the tendency of the low-speed streaks to appear in regions where the polymers are highly elongated, while the high-speed streaks mostly appear in regions where polymers are not stretched.

Conclusions
The combined effect of finite Weissenberg and Bingham numbers in EVP flows has been investigated by direct numerical simulations in this work, and compared with the two limiting cases of VE and VP flows, where VP refers to the near-viscoplastic limit (Wi = 0.01) analysed in our previous work Rosti et al. (2018a). We study both laminar and turbulent flows, where the turbulent case has Re τ = 180 in the Newtonian limit. A wide range of Weissenberg and Bingham numbers is investigated in this work: Wi = 0-16 and Bi = 0-24.
Previously, we observed that the VP flow has a moderate drag reduction compared with Newtonian flow at intermediate Bingham numbers, while at higher Bingham numbers the flow relaminarises and drag increases. However, the EVP flow at finite Weissenberg numbers achieves higher drag-reduction values, up to 70%, than both viscoplastic and viscoelastic flows at the same Re and Wi. Moreover, the drag of EVP flow keeps decreasing after relaminarisation at high Bingham numbers.
The VP flow became laminar at Bi ≥ 2.8. The EVP flow is still turbulent at Bi = 5.6 at moderate Wi, but becomes laminar at Bi = 5.6 for higher Wi (Wi ≥ 8). The moderate value of Wi triggers instabilities that reinstate turbulence, since the thickness of the unyielded central plug region is considerably reduced compared with the viscoplastic case. At higher Wi, the central plug region becomes too thin to create instabilities leading to laminarisation.
When analysing the flow statistics, the streamwise component of Reynolds stress tensor increases while the other components decrease compared with the Newtonian flow. This effect is seen for EVP, VE and VP flows. However, it is strongest for the EVP flow. At Wi = 4, Bi = 5.6 the components of the EVP flow deviate more from the Newtonian case than VE flow at Wi = 16, or any of the turbulent VP cases.
To characterise the intermittent dynamics relating to the drag-reducing property of polymeric flows (Xi & Graham 2010, 2012b, we divided each flow into hibernating time intervals and periods of active turbulence. The VE flow at Wi = 16.0 and Bi = 0, EVP flow at Wi = 8.0 and Bi = 2.8 and VP flow at Wi = 0.01 and Bi = 1.4 were considered. The VP flow showed hibernation phenomenon with small drag-reduction variation. However, the EVP flow hibernates more than any of the others, with the velocity profile reaching close to the Virk profile. All EVP stress components were highest in active turbulence and lowest in hibernation stages. The probability of the fluid to be unyielded was clearly divided into two stages: low probability in active stages and high probability in hibernation.
Finally, we analysed the spatial flow structures and related them to the local elongation of the polymers. As previously noted, both elasticity and plasticity increase the streamwise coherence in the flow. In the EVP flows, we observed that polymer elongation increases with increasing Bingham number, again showing that elastic effects are stronger at finite Bingham numbers than in a purely elastic flow. In both VE and EVP flows, low-speed streaks appear in regions where polymers are more stretched, and high-speed streaks in regions where polymers are less stretched. However, plasticity in EVP flows reduces the turbulent dissipation and hence the fragmentation of both the low-and high-speed streaks equally, while in VE flows high-speed streaks remain fragmented.