1 Introduction
Many fluids in nature and industrial applications exhibit a non-Newtonian behaviour, i.e. a nonlinear relation between the shear stress and the shear rate, such as shear thinning, shear thickening, yield stress, thixotropic, shear banding and viscoelastic behaviours. Moreover, several non-Newtonian features are often present simultaneously. Here, we focus on elastoviscoplastic fluids, i.e. complex non-Newtonian fluids that can exhibit simultaneously elastic, viscous and plastic properties. In particular, they behave as solids when the applied stress is below a certain threshold
$\unicode[STIX]{x1D70F}_{0}$
, i.e. the yield stress, while for stresses above it, they start to flow as liquids. In this context, the aim of this work is to explore and better understand the laminar and turbulent flow of an elastoviscoplastic fluid by means of numerical simulations. Indeed, turbulent flows of elastoviscoplastic fluids occur in many industrial settings, such as petroleum, paper, mining and sewage treatment (Hanks Reference Hanks1963, Reference Hanks1967; Maleki & Hormozi Reference Maleki and Hormozi2018).
1.1 Stability of yield stress fluids
Several studies have been devoted to the stability of yield stress fluids (Nouar & Frigaard Reference Nouar and Frigaard2001; Metivier, Nouar & Brancher Reference Metivier, Nouar and Brancher2005; Nouar et al. Reference Nouar, Kabouya, Dusek and Mamou2007; Nouar & Bottaro Reference Nouar and Bottaro2010; Bentrad et al. Reference Bentrad, Esmael, Nouar, Lefevre and Ait-Messaoudene2017). The first study on the stability of viscoplastic fluid flows was reported by Frigaard, Howison & Sobey (Reference Frigaard, Howison and Sobey1994), who studied the linear stability of a Bingham fluid in a plane channel flow. More recently, Nouar et al. (Reference Nouar, Kabouya, Dusek and Mamou2007) performed a modal and non-modal linear stability analysis of the flow of a Bingham fluid also in a plane channel; these authors showed that the flow is always linearly stable and that the optimal disturbance for moderate/high Bingham number is oblique, i.e. not aligned with the Cartesian coordinate axes as in Newtonian fluids. A key results arising from the linear stability analysis is that the regions where the stress is below the yield stress value remain unyielded for linear perturbations, a fact that can lead to interesting mathematical anomalies. For example, Metivier et al. (Reference Metivier, Nouar and Brancher2005) showed that the critical Reynolds number for linear stability is different when the Bingham number tends to zero, compared to a Newtonian fluid with a null value of Bingham number. Thus, the authors suggest that the passage to the Newtonian limit of a yield stress fluid is ill defined in terms of stability. Besides linear analysis, fully nonlinear (energy) stability results were derived in Nouar & Frigaard (Reference Nouar and Frigaard2001). These authors showed that the critical Reynolds number for transition increases with the Bingham number; however they also observed that the energy stability results are very conservative. Moreover, since for yield stress fluids the nonlinearity of the problem is not simply in the inertial terms, but also in the shear stress and in the existence of unyielded plug regions, the gap between linear and nonlinear theories is much wider than with Newtonian fluids. While in Newtonian fluids weakly nonlinear theories provide useful insights, in the case of viscoplastic fluids, these methods are algebraically more complicated and only Metivier, Nouar & Brancher (Reference Metivier, Nouar and Brancher2010) have performed this type of analysis for a Rayleigh–Bénard–Poiseuille flow, finding that the range of validity of an amplitude equation is fairly limited. Only a small number of studies on the stability of more complicated geometries exist: recently, nonlinear (energy) stability analysis has been extended to multi-layer flows of yield stress and viscoelastic fluids by Moyers-Gonzalez, Frigaard & Nouar (Reference Moyers-Gonzalez, Frigaard and Nouar2004) and Hormozi & Frigaard (Reference Hormozi and Frigaard2012). Recently, in order to identify possible paths to transition Nouar & Bottaro (Reference Nouar and Bottaro2010) perturbed the base flow slightly, and found that very weak defects are indeed capable of exciting exponentially amplified streamwise travelling waves. Finally, Kanaris, Kassinos & Alexandrou (Reference Kanaris, Kassinos and Alexandrou2015) performed numerical simulations of a Bingham fluid flowing past a confined circular cylinder to study the viscoplastic effects in the wake-transition regime.
1.2 Friction losses and drag reduction
Turbulent flows of generalized Newtonian fluids occur in many industrial process. Despite the numerous applications, it has not been possible to estimate the force needed to drive a complex fluid yet, while in a Newtonian flow the pressure drop can be accurately predicted as a function of the Reynolds number, both in laminar and turbulent flows (Pope Reference Pope2001), and for different properties of the wall surface, e.g. roughness (Orlandi & Leonardi Reference Orlandi and Leonardi2008), porosity (Breugem, Boersma & Uittenbogaard Reference Breugem, Boersma and Uittenbogaard2006; Rosti, Cortelezzi & Quadrio Reference Rosti, Cortelezzi and Quadrio2015; Rosti, Brandt & Pinelli Reference Rosti, Brandt and Pinelli2018b ) and elasticity (Rosti & Brandt Reference Rosti and Brandt2017). This is due to the complexity of such flows where additional parameters become relevant, such as the yield stress value (above which the material flows), the relaxation time, the ratio of the solvent to the total viscosity, …; each of these parameters may affect the overall flow dynamics in different and sometimes surprising ways. Some work has been done on measuring and trying to estimate the hydraulic pressure losses in practical applications (Ryan & Johnson Reference Ryan and Johnson1959; Hanks Reference Hanks1963, Reference Hanks1967; Hanks & Dadia Reference Hanks and Dadia1971), with the most popular phenomenological approach suggested by Metzner & Reed (Reference Metzner and Reed1955). These authors provide a closure for the pressure drop as a function of a generalized Reynolds number defined using the local power-law parameters, subsequently extended to yield stress fluids by Pinho & Whitelaw (Reference Pinho and Whitelaw1990) and Founargiotakis, Kelessidis & Maglione (Reference Founargiotakis, Kelessidis and Maglione2008). Rudman et al. (Reference Rudman, Blackburn, Graham and Pullum2004) performed numerical simulations of a turbulent pipe flow of shear-thinning fluids and compared their results with the pressure drop closure discussed above, finding a decent agreement although with some differences.
There exists a large literature on the turbulent flow with polymer additives, with the main focus being the drag reduction (Logan Reference Logan1972; Pinho & Whitelaw Reference Pinho and Whitelaw1990; Escudier & Presti Reference Escudier and Presti1996; Den Toonder et al. Reference Den Toonder, Hulsen, Kuiken and Nieuwstadt1997; Beris & Dimitropoulos Reference Beris and Dimitropoulos1999; Escudier, Presti & Smith Reference Escudier, Presti and Smith1999; Warholic, Massah & Hanratty Reference Warholic, Massah and Hanratty1999; Escudier & Smith Reference Escudier and Smith2001; Dubief et al. Reference Dubief, White, Terrapon, Shaqfeh, Moin and Lele2004, Reference Dubief, Terrapon, White, Shaqfeh, Moin and Lele2005; Escudier et al. Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005; Escudier, Nickson & Poole Reference Escudier, Nickson and Poole2009; Xi & Graham Reference Xi and Graham2010; Owolabi, Dennis & Poole Reference Owolabi, Dennis and Poole2017; Shahmardi et al. Reference Shahmardi, Zade, Ardekani, Poole, Lundell, Rosti and Brandt2018). The interested reader is referred to the works by Berman (Reference Berman1978) and White & Mungal (Reference White and Mungal2008) for a thorough review on the subject.
1.3 Elastoviscoplastic fluid
Despite the numerous studies performed to analyse viscoelastic turbulent flows, much less attention has been given to viscoplastic and elastoviscoplastic fluids. Indeed, very few numerical works exist on fully turbulent flows of an elastoviscoplastic fluid, and to the best of our knowledge the only direct numerical simulations of the effect of a yield stress on a turbulent non-Newtonian flow were performed by Rudman & Blackburn (Reference Rudman and Blackburn2006) and Guang et al. (Reference Guang, Rudman, Chryss, Slatter and Bhattacharya2011). These authors simulated a yield–pseudoplastic fluid using the Herschel–Bulkley model and compared the results with experimental measurements. Although qualitative agreement was found, the simulation results strongly overpredict the flow velocity, and the authors were not able to find the source of the discrepancy. Their numerical results suggest that, as the yield stress increases, the mean velocity profile deviates more and more from the Newtonian one, and that the turbulent flow will be fully developed only for low values of the yield stress.
Many materials used in experiments, such as Carbopol solutions (i.e. a conventional yield stress test fluid) and liquid foams, exhibit simultaneously elastic, viscous and yield stress behaviour. Thus, in order to properly predict the behaviour of such materials, it is essential to model them as a fully elastoviscoplastic fluid, rather than an ideal yield stress fluid (e.g. using the Bingham or Herschel–Bulkley model). Recently, Saramito (Reference Saramito2007) proposed a new constitutive equation for elastoviscoplastic fluid flows, which reproduces a viscoelastic solid for stresses lower than the yield stress, and a viscoelastic Oldroyd-B fluid for stresses higher than the yield stress. Furthermore, in order to describe the yielding process it uses the von Mises yielding criterion, which has been also experimentally confirmed (Shaukat et al. Reference Shaukat, Kaushal, Sharma and Joshi2012; Martinie, Buggisch & Willenbacher Reference Martinie, Buggisch and Willenbacher2013). Cheddadi et al. (Reference Cheddadi, Saramito, Dollet, Raufaste and Graner2011) simulated the inertialess flow of an elastoviscoplastic fluid around a circular object using the model proposed by Saramito (Reference Saramito2007); these authors were able to capture the fore–aft asymmetry and also the overshoot of the velocity (negative wake) after the circular hindrance, which was previously observed experimentally by Dollet & Graner (Reference Dollet and Graner2007) for the flow of a liquid foam and by Putz et al. (Reference Putz, Burghelea, Frigaard and Martinez2008) who related this behaviour to the rheological properties of the fluid. Note that the Bingham model always predicts fore–aft symmetry and the lack of a negative wake, which is in contradiction with the aforementioned experimental observations. Recently, the loss of the fore–aft symmetry and the formation of the negative wake around a single particle sedimenting in a Carbopol solution was captured by the numerical calculations in Fraggedakis, Dimakopoulos & Tsamopoulos (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016) using the constitutive law by Saramito (Reference Saramito2007); their results are in a quantitative agreement with experimental observations obtained with a Carbopol gel.
The model proposed by Saramito (Reference Saramito2007) was extended by the same author to account for shear-thinning effects (Saramito Reference Saramito2009). The new model combines 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 (Reference Saramito2007). Apart from the models proposed by Saramito, many others exist in the literature. The interested reader is referred to Crochet & Walters (Reference Crochet and Walters1983), Balmforth, Frigaard & Ovarlez (Reference Balmforth, Frigaard and Ovarlez2014) and Saramito & Wachs (Reference Saramito and Wachs2016), Saramito (Reference Saramito2016) for a thorough review of models and numerical methods.
1.4 Outline
In this work, we present the first direct numerical simulations of both laminar and turbulent channel flows of an incompressible elastoviscoplastic fluid. In the laminar regime, a wide range of Reynolds numbers is investigated, while in the turbulent regime, we consider the bulk Reynolds number
$Re=2800$
. The non-Newtonian flow is simulated by solving the full unsteady incompressible Navier–Stokes equations coupled with the model proposed by Saramito (Reference Saramito2007) for the evolution of the additional elastoviscoplastic stress tensor. In § 2, we first discuss the flow configuration and the governing equations, and then present the numerical methodology used. A validation of the numerical implementation is reported in § 2.2, while the results on the laminar and on the fully developed turbulent channel flows are presented in § 3. In particular, we discuss the role of some of the parameters defining the elastoviscoplastic fluid, i.e. the Bingham number
$Bi$
and the viscosity ratio
$\unicode[STIX]{x1D6FD}$
. Finally, a summary of the main findings and conclusions are presented in § 4.
2 Formulation
We consider the laminar and turbulent flows of an incompressible elastoviscoplastic 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.
(a) Sketch of the computational domain. (b) Sketch of the mechanical model of the elastoviscoplastic fluid proposed by Saramito (Reference Saramito2007) and used in the present work.
The fluid motion is governed by the conservation of momentum and the incompressibility constraint:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}u_{i}u_{j}}{\unicode[STIX]{x2202}x_{j}}=\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D70C}$
is the fluid density and
$\unicode[STIX]{x1D70E}_{ij}$
the total Cauchy stress tensor, which is written as 
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{ij}=-p\unicode[STIX]{x1D6FF}_{ij}+2\unicode[STIX]{x1D707}_{f}\unicode[STIX]{x1D60B}_{ij}+\unicode[STIX]{x1D70F}_{ij},\end{eqnarray}$$
where
$p$
is the pressure,
$\unicode[STIX]{x1D707}_{f}$
the fluid molecular dynamic viscosity (also called solvent viscosity),
$\unicode[STIX]{x1D6FF}$
is the Kronecker delta and
$\unicode[STIX]{x1D60B}_{ij}$
the strain rate tensor defined as
$$\begin{eqnarray}\unicode[STIX]{x1D60B}_{ij}=\frac{1}{2}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right).\end{eqnarray}$$
In (2.2),
$\unicode[STIX]{x1D70F}_{ij}$
is the additional elastoviscoplastic stress tensor which accounts for the non-Newtonian behaviour of the fluid, here described by the model proposed by Saramito (Reference Saramito2007). A one-dimensional schematic of the mechanical behaviour of the model is shown in figure 1(b): when the stress
$\unicode[STIX]{x1D70E}$
is below the yield stress
$\unicode[STIX]{x1D70F}_{0}$
, the friction element is rigid and the system predicts only recoverable Kelvin–Voigt viscoelastic deformation due to the spring
$\unicode[STIX]{x1D705}$
and the viscous element
$\unicode[STIX]{x1D707}_{f}$
. When the stress exceeds the yield value
$\unicode[STIX]{x1D70F}_{0}$
, the friction element breaks and an additional viscous element
$\unicode[STIX]{x1D707}_{m}$
activates; the fluid then behaves as an Oldroyd-B viscoelastic fluid. Thus, the total strain rate
$\dot{\unicode[STIX]{x1D700}}$
is shared between an elastic contribution
$\dot{\unicode[STIX]{x1D700}}_{e}$
and a plastic one
$\dot{\unicode[STIX]{x1D700}}_{p}$
(Cheddadi et al.
Reference Cheddadi, Saramito, Dollet, Raufaste and Graner2011). The following limits can be identified: the model reduces to the Oldroyd-B model for
$\unicode[STIX]{x1D70F}_{0}=0$
, the Bingham model is recovered for
$\unicode[STIX]{x1D706}=0$
and the fluid is Newtonian with a total viscosity
$\unicode[STIX]{x1D707}$
equal to
$\unicode[STIX]{x1D707}_{f}+\unicode[STIX]{x1D707}_{m}$
for
$\unicode[STIX]{x1D70F}_{0}=0$
and
$\unicode[STIX]{x1D706}=0$
. The instantaneous values of all the components of the stress
$\unicode[STIX]{x1D70F}_{ij}$
are found by solving the following objective and frame-independent transport equation
$$\begin{eqnarray}\unicode[STIX]{x1D706}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}u_{k}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{k}}-\unicode[STIX]{x1D70F}_{kj}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}}-\unicode[STIX]{x1D70F}_{ik}\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{k}}\right)+\max \left(0,\frac{|\unicode[STIX]{x1D70F}_{d}|-\unicode[STIX]{x1D70F}_{0}}{|\unicode[STIX]{x1D70F}_{d}|}\right)\unicode[STIX]{x1D70F}_{ij}=2\unicode[STIX]{x1D707}_{m}\unicode[STIX]{x1D60B}_{ij}.\end{eqnarray}$$
Here,
$\unicode[STIX]{x1D706}$
is the relaxation time,
$\unicode[STIX]{x1D707}_{m}$
is an additional viscosity,
$\unicode[STIX]{x1D70F}_{0}$
the yield stress and
$|\unicode[STIX]{x1D70F}_{d}|$
represents the second invariant of the deviatoric part of the added stress tensor, i.e.
$|\unicode[STIX]{x1D70F}_{d}|=\sqrt{1/2\unicode[STIX]{x1D70F}_{ij}^{d}\unicode[STIX]{x1D70F}_{ij}^{d}}$
. Note that, the first four terms in the left-hand side of the previous equation are the upper convected derivative of the elastoviscoplastic stress tensor (Gordon & Schowalter Reference Gordon and Schowalter1972). The elastoviscoplastic parameters
$\unicode[STIX]{x1D707}_{f}$
,
$\unicode[STIX]{x1D707}_{m}$
,
$\unicode[STIX]{x1D706}$
and
$\unicode[STIX]{x1D70F}_{0}$
can be obtained by experimental data following the procedure detailed by Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016), based on the determination of the linear material functions, i.e. the storage modulus
$G^{\prime }$
and the loss modulus
$G^{\prime \prime }$
.
The previous set of equations can be rewritten in a non-dimensional form as
$$\begin{eqnarray}\displaystyle & \displaystyle Re\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}u_{i}u_{j}}{\unicode[STIX]{x2202}x_{j}}\right)=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}(-p\unicode[STIX]{x1D6FF}_{ij}+2\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D60B}_{ij}+\unicode[STIX]{x1D70F}_{ij}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Wi\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}u_{k}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{k}}-\unicode[STIX]{x1D70F}_{kj}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{k}}-\unicode[STIX]{x1D70F}_{ik}\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{k}}\right)+\max \left(0,\frac{|\unicode[STIX]{x1D70F}_{d}|-Bi}{|\unicode[STIX]{x1D70F}_{d}|}\right)\unicode[STIX]{x1D70F}_{ij}=2(1-\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D60B}_{ij},\qquad & \displaystyle\end{eqnarray}$$
$Re$
, the Weissenberg number
$Wi$
, the Bingham number
$Bi$
and the viscosity ratio
$\unicode[STIX]{x1D6FD}$
. The Reynolds number is the ratio of inertia and viscous forces
$Re=\unicode[STIX]{x1D70C}UL/\unicode[STIX]{x1D707}_{0}$
, the Bingham number the ratio of the yield and viscous stresses
$Bi=\unicode[STIX]{x1D70F}_{0}L/\unicode[STIX]{x1D707}_{0}U$
, the Weissenberg number the ratio of the elastic and viscous forces
$Wi=\unicode[STIX]{x1D706}U/L$
(Poole Reference Poole2012) and the viscosity ratio
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D707}_{f}/\unicode[STIX]{x1D707}_{0}$
the ratio between the fluid viscosity
$\unicode[STIX]{x1D707}_{f}$
and the reference one
$\unicode[STIX]{x1D707}_{0}$
. In the previous definitions,
$U$
and
$L$
are a characteristic velocity and length scales of the flow,
$\unicode[STIX]{x1D70C}$
the fluid density and
$\unicode[STIX]{x1D707}_{0}$
a characteristic viscosity, set equal to the total viscosity, i.e.
$\unicode[STIX]{x1D707}_{0}=\unicode[STIX]{x1D707}_{f}+\unicode[STIX]{x1D707}_{m}$
. Note that, the choice of the characteristic viscosity is an open topic of discussion in the community, with the most common choice being the total viscosity
$\unicode[STIX]{x1D707}_{0}$
.2.1 Numerical discretization
The equations of motion are solved with an extensively validated in-house code (Picano, Breugem & Brandt Reference Picano, Breugem and Brandt2015; Rosti & Brandt Reference Rosti and Brandt2017, Reference Rosti and Brandt2018; Rosti, Brandt & Mitra Reference Rosti, Brandt and Mitra2018a
). Equations (2.1) and (2.4) are solved on a staggered uniform grid with velocities located on the cell faces and all the other variables (pressure, stress and material component properties) at the cell centres. All the spatial derivatives are approximated with second-order centred finite differences except for the advection term in (2.4) where the fifth-order WENO (weighted essentially non-oscillatory) scheme is adopted (Shu Reference Shu2009; Sugiyama et al.
Reference Sugiyama, Ii, Takeuchi, Takagi and Matsumoto2011). The time integration is performed with a fractional-step method (Kim & Moin Reference Kim and Moin1985), where all the terms in the evolution equations are advanced in time with a third-order explicit Runge–Kutta scheme except for the elastoviscoplastic stress terms which are advanced with the Crank–Nicolson method; moreover, a fast Poisson solver is used to enforce the condition of zero divergence for the velocity field. Note that the choice of an explicit time integration is typically preferred for high Reynolds number turbulent flows. In particular, to solve the system of governing equations, we perform the following steps (see also Min, Yoo & Choi Reference Min, Yoo and Choi2001; Dubief et al.
Reference Dubief, Terrapon, White, Shaqfeh, Moin and Lele2005): (i) the elastoviscoplastic stress tensor
$\unicode[STIX]{x1D70F}_{ij}$
is updated by solving (2.4); (ii) the Navier–Stokes equations (2.1) are advanced in time by first solving the momentum equation (prediction step), then by solving a Poisson equation for the projection variable and finally by correcting the velocity and pressure to make the velocity field divergence free (correction step).
2.2 Code validation
The present implementation for single and multiphase flows of an elastoviscoplastic fluid has been extensively validated in Izbassarov et al. (Reference Izbassarov, Rosti, Niazi, Sarabian, Hormozi, Brandt and Tammisola2018), where the details of the algorithm are discussed in further detail. Nonetheless, we report here two validation cases for the sake of completeness.
(a) Time evolution of
$\unicode[STIX]{x1D70F}_{11}-\unicode[STIX]{x1D70F}_{22}$
(red) and
$\unicode[STIX]{x1D70F}_{12}$
(blue) in a stationary shear flow. The stress components are normalized with
$\unicode[STIX]{x1D707}_{0}\dot{\unicode[STIX]{x1D6FE}}$
. (b) Time evolution of the shear stress
$\unicode[STIX]{x1D70F}_{12}$
in an oscillating shear flow with
$Bi=0$
(red) and
$Bi=300$
(blue). The stress components are normalized with
$\unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D6FE}_{0}\unicode[STIX]{x1D714}_{0}$
. In both panels, solid lines are used for our numerical results and symbols for the analytical solution reported by Saramito (Reference Saramito2007).
First, we consider a simple constant shear flow, with the shear rate
$\dot{\unicode[STIX]{x1D6FE}}_{0}$
: the Weissenberg number is fixed to
$Wi=\unicode[STIX]{x1D706}\dot{\unicode[STIX]{x1D6FE}}_{0}=1$
, the Bingham number
$Bi=\unicode[STIX]{x1D70F}_{0}/(\unicode[STIX]{x1D707}_{0}\dot{\unicode[STIX]{x1D6FE}})=1$
and the viscosity ratio
$\unicode[STIX]{x1D6FD}=1/9$
. The time evolution of
$\unicode[STIX]{x1D70F}_{11}-\unicode[STIX]{x1D70F}_{22}$
(the normal stress difference) and
$\unicode[STIX]{x1D70F}_{12}$
(the wall-normal shear stress) are reported in figure 2(a) with red and blue lines, respectively. We observe that, initially both the stress components grow linearly, but when the stress level is above a threshold, i.e. the yield stress, the growth stops and they reach a plateau, as expected in the yielded state. As shown in the figure, we find a very good agreement with the analytical results by Saramito (Reference Saramito2007) depicted with symbols of the same colours.
Next, we consider a time-periodic uniform shear flow, i.e.
$\unicode[STIX]{x1D6FE}_{0}\sin (\unicode[STIX]{x1D714}_{0}t)$
, where
$\unicode[STIX]{x1D6FE}_{0}$
is the strain amplitude and
$\unicode[STIX]{x1D714}_{0}$
the angular frequency of the oscillations. The Weissenberg number is
$Wi=\unicode[STIX]{x1D706}\unicode[STIX]{x1D714}_{0}=0.1$
and two Bingham numbers
$Bi=\unicode[STIX]{x1D70F}_{y}/(\unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D6FE}_{0}\unicode[STIX]{x1D714}_{0})$
are considered:
$Bi=0$
and
$300$
. Note that, when
$Bi=0$
, the material behaves like a viscoelastic fluid, and when
$Bi=300$
as an elastic solid. The viscosity ratio
$\unicode[STIX]{x1D6FD}$
is null in both cases, i.e.
$\unicode[STIX]{x1D707}_{f}=0$
. The evolution of
$\unicode[STIX]{x1D70F}_{12}$
is plotted in figure 2(b) for the two cases (red
$Bi=0$
and blue
$Bi=300$
) and compared with the analytical solution provided by Saramito (Reference Saramito2007), shown with symbols in the figure. Again, an excellent agreement is found.
2.3 Numerical set-up
For all the cases considered hereafter, the equations of motion are discretized by using
$1728\times 576\times 864$
grid points on a computational domain of size
$6h\times 2h\times 3h$
in the streamwise, wall-normal and spanwise directions. The spatial resolution has been chosen in order to properly resolve the turbulent scales, as well as the unyielded plug regions which form intermittently in the domain. In the high Reynolds number simulations at
$Re_{b}=2800$
, the resolution satisfies the constraint
$\unicode[STIX]{x0394}x^{+}=\unicode[STIX]{x0394}y^{+}=\unicode[STIX]{x0394}z^{+}<0.6$
, where the superscript
$^{+}$
indicates the wall units defined in the next section. In one of the simulations at
$Re_{b}=2800$
, a grid refinement study was performed using
$2160\times 720\times 1080$
grid points in the streamwise, wall-normal and spanwise directions (around
$25\,\%$
more in each direction); the difference in the resulting friction coefficient
$C_{f}$
was less than
$2\,\%$
. Note that, in the low Reynolds fully laminar cases, the spatial resolution was relaxed and the domain size in the homogeneous directions reduced.
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. All the turbulent flows are initialized with a fully developed channel flow with zero elastoviscoplastic added stress (
$\unicode[STIX]{x1D70F}_{ij}=0$
). After the flow has reached statistically steady state, the calculations are continued for an interval of
$500h/U_{b}$
time units, during which around
$100$
full flow fields are stored for further statistical analysis. To verify the convergence of the statistics, we have computed them using a different number of samples and verified that the differences are negligible.
3 Results
We study both laminar and turbulent channel flows of an elastoviscoplastic fluid, together with the baseline Newtonian cases. All the simulations are performed at a constant flow rate, so that the flow Reynolds number based on the bulk velocity is fixed, i.e.
$Re=\unicode[STIX]{x1D70C}U_{b}h/\unicode[STIX]{x1D707}_{0}$
, where the bulk velocity
$U_{b}$
is the average value of the mean velocity computed across the whole domain and
$\unicode[STIX]{x1D707}_{0}$
is the total viscosity, i.e.
$\unicode[STIX]{x1D707}_{0}=\unicode[STIX]{x1D707}_{f}+\unicode[STIX]{x1D707}_{m}$
. In the present work, consistently with choosing
$U_{b}$
as the characteristic velocity, we opt for enforcing the constant flow rate condition; hence, the necessary value of the instantaneous streamwise pressure gradient is determined at every time step. This choice facilitates the comparison between the non-Newtonian and Newtonian flows. In the laminar regime, the bulk Reynolds number is varied between
$0.1$
and
$2800$
, where the corresponding baseline Newtonian solutions are known analytically; in the turbulent regime, the bulk Reynolds number is fixed to
$2800$
, corresponding to a nominal friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F}}h/\unicode[STIX]{x1D707}_{0}=180$
for a Newtonian fluid, being
$u_{\unicode[STIX]{x1D70F}}$
the friction velocity defined later on. In the turbulent case, we compare our Newtonian solution with the seminal work of Kim, Moin & Moser (Reference Kim, Moin and Moser1987).
The properties of the elastoviscoplastic fluids are chosen as follows: the Weissenberg number
$Wi=\unicode[STIX]{x1D706}U_{b}/h$
is fixed in all the simulations to
$0.01$
in order to limit the role of fluid elasticity in this first study of elastoviscoplastic flows; the Bingham number
$Bi=\unicode[STIX]{x1D70F}_{0}h/\unicode[STIX]{x1D707}_{0}U_{b}$
is varied in the range between
$0$
and
$1000$
in the laminar cases (
$0$
,
$0.1$
,
$1$
,
$10$
,
$100$
and
$1000$
), and between
$0$
and
$2.8$
(
$0$
,
$0.28$
,
$0.7$
,
$1.4$
and
$2.8$
) in the turbulent cases, which are computationally significantly more expensive than Newtonian turbulence. Finally, all the cases have been studied for three different viscosity ratios:
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D707}_{f}/\unicode[STIX]{x1D707}_{0}=0.25$
,
$0.5$
and
$0.95$
. Overall, we have performed
$108$
laminar and
$15$
turbulent simulations.
Viscous units, used above to express the spatial resolution, will be often employed in the following; they are indicated by the superscript
$^{+}$
, and are built using the friction velocity
$u_{\unicode[STIX]{x1D70F}}$
as the velocity scale and the viscous length
$\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}=\unicode[STIX]{x1D708}/u_{\unicode[STIX]{x1D70F}}$
as the length scale. For a Newtonian turbulent channel flow, the dimensionless friction velocity is defined as
$$\begin{eqnarray}u_{\unicode[STIX]{x1D70F}}=\sqrt{\frac{1}{Re_{b}}\left.\frac{\text{d}\overline{u}}{\text{d}y}\right|_{y=0}},\end{eqnarray}$$
where
$\overline{u}$
is the mean velocity, and the derivative is taken at
$y=0$
, the location of the wall. When the fluid is non-Newtonian, equation (3.1) must be modified to account for the elastoviscoplastic shear stress that is in general non-zero at the wall. Similarly to previous works with polymers (Shahmardi et al.
Reference Shahmardi, Zade, Ardekani, Poole, Lundell, Rosti and Brandt2018), we define
$$\begin{eqnarray}u_{\unicode[STIX]{x1D70F}}=\sqrt{\left.\left(\frac{1}{Re_{b}}\frac{\text{d}\overline{u}}{\text{d}y}+\overline{\unicode[STIX]{x1D70F}}_{12}\right)\right|_{y=0}}.\end{eqnarray}$$
Note that, the actual value of the friction velocity in our simulations is computed from the friction coefficient, found by the driving streamwise pressure gradient, rather than from its definition, i.e.
$u_{\unicode[STIX]{x1D70F}}=\sqrt{-\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D70C}\,\text{d}\overline{p}/\text{d}x}$
.
3.1 Laminar flow
We start our analysis by considering the laminar flow of an elastoviscoplastic fluid. First, we consider the effect of the Bingham number on the frictional resistance of the flow quantified by the Fanning friction factor
$f$
, defined as
$2\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D70C}U_{b}^{2}$
being
$\unicode[STIX]{x1D70F}_{w}$
the total wall shear stress including both the viscous and elastoviscoplastic contributions. Figure 3(a) shows the Fanning friction factor
$f$
as a function of the Reynolds number in the case with
$\unicode[STIX]{x1D6FD}=0.95$
,
$Wi=0.01$
and for various Bingham numbers. In particular, the grey, orange, brown, purple, cyan and gold lines are used for
$Bi=0$
,
$0.1$
,
$1$
,
$10$
,
$100$
and
$1000$
, respectively. In the figure we also show the Newtonian analytical solution
$f=6/Re_{b}$
with a black line. The results clearly show that all the non-Newtonian fluids have the same slopes as the reference Newtonian case, but with increasing
$f$
as the Bingham number
$Bi$
increases. As expected, the case with
$Bi=0$
(grey line) is almost indistinguishable from the Newtonian flow, since the elastic effects are small for the low value of the Weissenberg number
$Wi$
chosen. The results shown here are consistent with the experimental measurements reported by Guzel, Frigaard & Martinez (Reference Guzel, Frigaard and Martinez2009), who also found a linear relation between the Reynolds number and the friction factor in a laminar pipe flow.
Fanning friction factor
$f$
as a function of (a) the bulk Reynolds number
$Re_{b}$
for
$\unicode[STIX]{x1D6FD}=0.95$
and (b) the viscosity ratio
$\unicode[STIX]{x1D6FD}$
for
$Re=1$
, for different Bingham numbers
$Bi$
. In all the elastoviscoplastic cases, the Weissenberg number is fixed to
$Wi=0.01$
. Grey, orange, brown, purple, cyan and gold colours are used for
$Bi=0$
,
$0.1$
,
$1$
,
$10$
,
$100$
and
$1000$
, respectively, and the black line is the Newtonian analytical solution. Note that, the lines in the graphs are simple connections between the available data points.
Figure 3(b) shows the effect of
$\unicode[STIX]{x1D6FD}$
on the Fanning friction factor
$f$
at
$Re_{b}=1$
. We find that
$f$
decreases nonlinearly with the viscosity ratio
$\unicode[STIX]{x1D6FD}$
, and that the dependency on
$\unicode[STIX]{x1D6FD}$
increases with the Bingham number
$Bi$
. The increase in the friction factor, due to the increase of wall shear stress, comes from the change of the laminar streamwise velocity profile
$u$
shown in figure 4(a) as a function of the wall-normal distance
$y$
, with
$u$
and
$y$
being normalized with the bulk velocity
$U_{b}$
and
$h$
, respectively. Again, we observe that the viscoelastic flow (
$Bi=0$
) almost perfectly overlaps with the Newtonian solution due to the very low Weissenberg number considered in this study. As expected, as the Bingham number
$Bi$
increases, we note the appearance of a region in the middle of the channel with a uniform velocity, i.e. a plug is formed away from the walls flowing with uniform velocity; this corresponds to the region where the fluid is not yielded and behaves as an elastic solid. Consequently, the centreline velocity
$U_{c}=u(y=h)$
reduces and the wall shear increases for mass conservation. The volume of the unyielded fluid, denoted
$Vol_{s}$
, grows with the Bingham number
$Bi$
from
$0\,\%$
for
$Bi=0$
up to
$87\,\%$
for
$Bi=1000$
, as shown in figure 4(b).
(a) Mean streamwise velocity profile
$\overline{u}$
as a function of the wall-normal distance
$y$
. (b) Percentage of the unyielded volume
$Vol_{s}$
as a function of the Bingham number
$Bi$
. The Reynolds number is equal to
$1$
, and the colour scheme is the same as in figure 3. Note that, the lines in the graphs are simple connections between the available data points.
Finally, we provide a fit to our numerical data for the Fanning friction factor
$f$
(figure 3). In general, the Fanning friction factor
$f$
of an elastoviscoplastic fluid in a channel flow is a function of inertia (
$Re_{b}$
), elasticity (
$Wi$
), plasticity (
$Bi$
) and viscosity ratio
$\unicode[STIX]{x1D6FD}$
, i.e.
$f={\mathcal{F}}(Re,Wi,Bi,\unicode[STIX]{x1D6FD})$
. In our study the Weissenberg number is fixed to a very low value (
$Wi=0.01$
), thus we drop its dependency. A very good agreement with our data is found when using the following expression
$$\begin{eqnarray}f=\frac{6+{\mathcal{C}}\sqrt{Bi}}{Re_{b}},\end{eqnarray}$$
where
${\mathcal{C}}$
is a fit parameter which depends on
$\unicode[STIX]{x1D6FD}$
:
${\mathcal{C}}=2.47$
for
$\unicode[STIX]{x1D6FD}=0.25$
,
$7.55$
for
$\unicode[STIX]{x1D6FD}=0.5$
and
$8.86$
for
$\unicode[STIX]{x1D6FD}=0.95$
. Equation (3.3) clearly recovers the Newtonian analytical solution for
$Bi=0$
, and provides an error below
$2\,\%$
for all the elastoviscoplastic results. It is worth noticing, that an analogous expression was found by De Vita et al. (Reference De Vita, Rosti, Izbassarov, Duffo, Tammisola, Hormozi and Brandt2018) for the flow of an elastoviscoplastic fluid through a porous media, with the same dependency on
$Re_{b}$
and
$Bi$
.
3.2 Turbulent flow
Next, we examine the turbulent flow cases, all at a fixed bulk Reynolds number
$Re_{b}=2800$
and Weissenberg number
$Wi=0.01$
. The turbulent purely viscoelastic flow with
$Bi=0$
has a Fanning friction factor
$f$
higher than its laminar counterpart, rising by
$300\,\%$
from
$0.002$
in the laminar case to
$0.008$
in the turbulent one, as shown in figure 5(a) for the case with
$\unicode[STIX]{x1D6FD}=0.95$
. As the Bingham number increases, the Fanning friction factor progressively decreases, with an opposite trend to the laminar elastoviscoplastic cases. The decrease of
$f$
is small for the two lowest
$Bi$
(
$-1.4\,\%$
for
$Bi=0.28$
and
$-3.1\,\%$
for
$Bi=0.7$
), moderate for the intermediate
$Bi$
(
$-7.7\,\%$
for
$Bi=1.4$
) and large for the highest
$Bi$
(
$-41\,\%$
for
$Bi=2.8$
). For the highest
$Bi$
considered here, the turbulence cannot be sustained and hence
$f$
reaches its laminar value.
(a) Fanning friction factor
$f$
as a function of the Bingham number
$Bi$
for
$Re_{b}=2800$
and
$\unicode[STIX]{x1D6FD}=0.95$
. The points on the black line correspond to laminar flows, while those on the grey line to turbulent flows. (b) The friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}$
as a function of the Bingham number
$Bi$
, with the vertical error bars measuring the variance of
$Re_{\unicode[STIX]{x1D70F}}$
. The dashed, dash-dotted and solid lines are used for
$\unicode[STIX]{x1D6FD}=0.25$
,
$0.5$
and
$0.95$
, respectively. The Reynolds number is equal to
$Re_{b}=2800$
for all cases. The blue, magenta, red, orange and green colours are used for the turbulent cases with
$Bi=0$
,
$0.28$
,
$0.7$
,
$1.4$
and
$2.8$
, respectively, while the colour scheme for the laminar results is the same as in figure 3.
The friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}$
is depicted in figure 5(b) as a function of the Bingham number. Again, we observe a progressive decrease of
$Re_{\unicode[STIX]{x1D70F}}$
with
$Bi$
, corresponding to a net drag reduction when compared to a Newtonian turbulent channel flow at the same flow rate (horizontal grey line). The error bar in the figure represents the variance of the value, which is initially small, then grows suddenly for
$Bi=1.4$
(as discussed later), and finally becomes null for
$Bi=2.8$
. The different line styles used in figure 5(b) correspond to different values of the viscosity ratio
$\unicode[STIX]{x1D6FD}$
. We observe that, in the turbulent regime the results are almost independent of the viscosity ratio, both in terms of mean and root mean square (r.m.s.) values. To summarize, we identify three different regimes: (i) for low Bingham numbers (
${\lesssim}1$
) the friction Reynolds number decreases slowly, approximately linearly, with approximately constant r.m.s. values; (ii) for intermediate values,
$Re_{\unicode[STIX]{x1D70F}}$
decreases more than linearly and its r.m.s. increases; (iii) for high Bingham numbers (
${\gtrsim}2$
) the flow becomes stationary and fully laminar. Interestingly, these three separate regimes are found to be independent of the value of the viscosity ratio
$\unicode[STIX]{x1D6FD}$
.
We can define a Bingham number in wall units, i.e.
$Bi^{+}$
, as the ratio between the yield stress
$\unicode[STIX]{x1D70F}_{0}$
and the wall shear stress
$\unicode[STIX]{x1D70F}_{w}$
as
$Bi^{+}=\unicode[STIX]{x1D70F}_{0}/\unicode[STIX]{x1D70F}_{w}$
: from the results of our simulations we found that
$Bi^{+}=0$
,
$0.025$
,
$0.064$
,
$0.135$
and
$0.425$
for the turbulent simulations at
$Re=2800$
with
$Bi=0$
,
$0.28$
,
$0.7$
,
$1.4$
and
$2.8$
. Based on the values of
$Bi^{+}$
we have available, we can infer that the first regime depicted above holds for yield stress values that are below
$6\,\%$
the wall shear stress value, while the third regime holds for yield stress values above
$42\,\%$
the wall shear stress.
(a) Time history of the streamwise pressure drop
$\text{d}p/\text{d}x$
for different Bingham numbers
$Bi$
. (b) Probability of the fluid being unyielded
$P_{s}$
as a function of the wall-normal distance
$y$
. The percentages reported in the legend show the mean unyielded volume
$Vol_{s}$
. The blue, magenta, red, orange and green colours are used for
$Bi=0$
,
$0.28$
,
$0.7$
,
$1.4$
and
$2.8$
, respectively. The viscosity ratio
$\unicode[STIX]{x1D6FD}$
is fixed equal to
$0.95$
.
The time history of the instantaneous pressure drop along the channel
$\text{d}p/\text{d}x$
is shown in figure 6(a). This quantity represents the forcing term needed to drive the flow, which in the turbulent regime oscillates around a mean value in order to maintain a constant flow rate in the domain. We observe that for the low Bingham cases (
$Bi=0$
,
$0.28$
and
$0.7$
) the time histories of
$\text{d}p/\text{d}x$
are very similar, with only slightly different mean values; on the other hand, for
$Bi=1.4$
the mean value is further decreased while the amplitude of the oscillations increases. Finally, for
$Bi=2.8$
the pressure drop smoothly decays from the turbulent value imposed as initial condition to the final laminar value. Thus, the figure clearly confirms the differences between the three regimes highlighted above.
Contours of the instantaneous spanwise vorticity
$-\unicode[STIX]{x1D714}_{z}$
in an
$x{-}y$
plane (a,c,e,g,i) and of the streamwise vorticity
$\unicode[STIX]{x1D714}_{x}$
in a
$y{-}z$
plane (b,d,f,h,j). Colour scale ranges from
$-3U_{b}/h$
(blue) to
$3U_{b}/h$
(red). The brown areas represent the instantaneous regions where the flow is not yielded. The Bingham number
$Bi$
increases from top to bottom (
$Bi=0$
,
$0.28$
,
$0.7$
,
$1.4$
and
$2.8$
) and the viscosity ratio
$\unicode[STIX]{x1D6FD}$
is fixed equal to
$0.95$
.
In order to understand the physical origin of the three regimes, we start by showing visualizations of the instantaneous distributions of the regions where the flow is yielded and not yielded, see figure 7. In the wall-normal and cross-stream planes in the figure we also report colour contours of the spanwise (left column corresponding to a
$x{-}y$
plane) and streamwise vorticity (right column displaying
$y{-}z$
planes),
$\unicode[STIX]{x1D714}_{z}$
and
$\unicode[STIX]{x1D714}_{x}$
. At
$Bi=0$
we recognize the classic vorticity field of turbulent channel flows, with high vorticity levels at the walls, and the footprints of the classical turbulent streaky structures. For non-zero Bingham numbers, we see the appearance of unyielded regions – shown in brown – around the centre of the channel and far from the walls. These regions are mostly disconnected and with a limited spanwise length for the two lowest
$Bi$
(corresponding to the first regime), while for
$Bi=1.4$
and
$2.8$
the unyielded region extends over the full streamwise and spanwise directions. The bottom row of the figure clearly shows that the flow is fully laminar.
Intermittency of the yield/unyield process (F/S) as a function of time. The four panels correspond to probes located at
$x=3h$
,
$z=1.5h$
and different wall-normal distances: (a)
$y\approx 0.25h$
, (b)
$0.49h$
, (c)
$0.74$
and (d)
$0.98h$
. The lines in every panel are the results with different Bingham number, with the colour scheme being the same as in figure 6.
Figure 8 shows the instantaneous state of the fluid: S indicates an unyielded fluid and F a yielded one. This information is extracted by various numerical probes at different wall-normal locations. In particular, we show in the figure the results for four different wall-normal distances, i.e.
$y\approx 0.25h$
,
$0.49h$
,
$0.74$
and
$0.98h$
. The intermittent nature of the flow is evident; also it is clear that the fluid close to the wall is preferentially yielded, while close to the centreline it is mostly not yielded, even at low Bingham numbers.
As clearly indicated by the previous pictures, the flow and the yield/unyield process are inherently unsteady, thus we need to analyse the phenomenon in statistical terms. Going back to figure 6(b), we display the probability to have an unyielded region
$P_{s}$
as a function of the wall-normal distance
$y$
, together with the mean percentage of the unyielded volume
$Vol_{s}$
reported in the legend. The value
$P_{s}=1$
indicates a location where the material behaves as an elastic solid throughout the computational time, whereas the material behaves uniquely as a viscoelastic (Oldroyd-B) fluid when
$P_{s}=0$
. For
$Bi=0$
, we clearly have no solid anywhere, while as
$Bi$
increases, the probability of the fluid to be not yielded increases around the centreline, while still remains null in the near-wall region. Finally, for
$Bi=2.8$
when the flow is fully laminar, the probability of being yielded or unyielded is either
$0$
or
$1$
, with
$54\,\%$
of the total volume being not yielded. Note that, even for the Bingham number
$Bi=1.4$
, representative of the second regime, the probability to be unyielded in the middle of the channel is not exactly equal to unity. Indeed, as also shown in the third row of figure 7, instantaneous region where the fluid is yielded can appear around the centreline, thus decreasing the overall percentage. In particular, for
$Bi=0.28$
the probability of the fluid being yielded at the centreline is approximately
$85\,\%$
, for
$Bi=0.7$
is
$40\,\%$
and for
$Bi=1.4$
is
$8\,\%$
. This effect contributes to the highly unsteady and intermittent behaviour discussed in relation to the pressure drop and friction factor.
Mean streamwise velocity profile
$\overline{u}$
as a function of the wall-normal distance
$y$
in bulk (a) and wall units (b). The colour scheme is the same as in figure 6, with the addition of the black symbols used for the Newtonian case, taken from the results by Kim et al. (Reference Kim, Moin and Moser1987). The viscosity ratio
$\unicode[STIX]{x1D6FD}$
is fixed equal to
$0.95$
.
We now proceed by presenting the main flow statistics. Figure 9 shows the mean streamwise velocity component
$\overline{u}$
as a function of the wall-normal distance
$y$
. In (a) (in bulk units) we can again find the three different behaviours described above. Up to
$Bi\approx 1$
, the profiles are quite similar, with only little reductions of the wall shear and an increase of the centreline velocity as
$Bi$
grows. The difference with the Newtonian case becomes more noticeable for
$Bi=1.4$
, where a region with zero shear appears at the centreline. Finally the profile for
$Bi=2.8$
clearly differs from the other ones, with a large zero-shear region occupying more than
$50\,\%$
of the channel. In (b) the same velocity profiles are shown in wall units. In most of the turbulent cases we can identify three regions in the velocity profile, similarly to those found for a Newtonian turbulent channel flow (black symbols): first, the viscous sublayer for
$y^{+}<5$
where the variation of
$\overline{u}^{+}$
with
$y^{+}$
is linear; then, the so-called log-law region,
$y^{+}>30$
, where the variation of
$\overline{u}^{+}$
versus
$y^{+}$
is logarithmic; finally, the region between
$5$
and
$30$
wall units is called buffer layer and neither law holds here. As
$Bi$
increases, the extension of the inertial ranges reduces eventually disappearing, thus indicating the absence of an equilibrium range. By comparing the three regions discussed above present in the mean velocity profile and the extension of the unyielded region shown in figure 6, we can observe that the flow remains unyielded mostly in the logarithmic and outer layer, while it is always yielded in the viscous sublayer.
Wall-normal profiles of the different components of the Reynolds stress tensor, normalized with
$u_{\unicode[STIX]{x1D70F}}^{2}$
. (a–c) Show the diagonal components
$u^{\prime }u^{\prime }$
,
$v^{\prime }v^{\prime }$
and
$w^{\prime }w^{\prime }$
, while (d) shoes the cross-term
$u^{\prime }v^{\prime }$
. The colour scheme is the same as in figure 6.
We continue our comparison between the turbulent channel of a Newtonian and elastoviscoplastic fluid by analysing the wall-normal distribution of the diagonal component of the Reynolds stress tensor; these are shown in figure 10 together with the data from Kim et al. (Reference Kim, Moin and Moser1987) for the Newtonian case, represented with the black symbols
$+$
. Also in the Reynolds stress profiles, we find the distinction between the three regimes previously mentioned. First, for low
$Bi$
, the fluctuations are only slightly affected, with the differences being noticeable only in the buffer layer, while the profiles in the viscous sublayer and in the inertial range still show a good collapse in wall units. Then, for high Bingham numbers (
$Bi=1.4$
) the profiles undergo strong modifications, which are not limited to the buffer layer, but extend to the inertial range and viscous sublayer as well. Finally, the full laminarization of the flow for
$Bi=2.8$
is proved again by showing the null values assumed by all the Reynolds stress components, in the whole domain.
We also observe a clear trend, although not linear, of the Reynolds stress components with the Bingham number: the streamwise component
$\overline{u^{\prime }u^{\prime }}$
increases, while the wall-normal
$\overline{v^{\prime }v^{\prime }}$
and spanwise
$\overline{w^{\prime }w^{\prime }}$
decrease. Also, all the peaks are displaced away from the wall, towards the centreline. In relative terms, the wall-normal and spanwise components are the most affected ones, decreasing by almost
$40\,\%$
. Finally, figure 10(d) depicts the wall-normal profile of the off-diagonal component of the Reynolds stress tensor
$\overline{u^{\prime }v^{\prime }}$
. Also this cross component is affected by the elastoviscoplasticity of the fluid in a similar fashion as the diagonal components. In particular, the maximum value decreases and moves away from the wall as the Bingham number increases. Nevertheless, the stress profiles still vary linearly between the two peaks of opposite sign close to each wall, but with different slopes (not shown here). The Reynolds stress modifications due to the elastoviscoplastic property of the fluid are similar to what observed in other drag reducing flows, such as the turbulent flow over riblets (García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011), the turbulent flow over anisotropic porous walls (Rosti & Brandt Reference Rosti and Brandt2018), and turbulent flows with polymers (Dubief et al.
Reference Dubief, White, Terrapon, Shaqfeh, Moin and Lele2004; Shahmardi et al.
Reference Shahmardi, Zade, Ardekani, Poole, Lundell, Rosti and Brandt2018). In general, the increased amplitude of the streamwise fluctuations, and the reduction of the other components, is usually associated with the strengthening of streaky structures above the wall, which is true also in the present case, as shown in the next paragraphs.
Wall-normal profiles of (a) the turbulent kinetic energy
${\mathcal{K}}=(u^{\prime 2}+v^{\prime 2}+w^{\prime 2})/2$
and of (b) the turbulent production
${\mathcal{P}}=-\overline{u^{\prime }v^{\prime }}\,\text{d}\overline{u}/\text{d}y$
, both normalized with the friction velocity
$u_{\unicode[STIX]{x1D70F}}$
. The colour scheme is the same as in figure 6, with the blue, magenta, red, orange and green colours used for the turbulent cases with
$Bi=0$
,
$0.28$
,
$0.7$
,
$1.4$
and
$2.8$
, respectively.
Wall-normal profiles of (a) the turbulent dissipation
$\unicode[STIX]{x1D700}=\unicode[STIX]{x1D707}\overline{\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}}$
in wall units and of (b) the shear effective viscosity
$\unicode[STIX]{x1D707}_{e}$
. The colour scheme is the same as in figure 6. The inset figure in (b) shows
$\unicode[STIX]{x1D707}_{e}$
as a function of the shear rate
$\dot{\unicode[STIX]{x1D6FE}}$
.
An overall view of the velocity fluctuations can be inferred by considering the turbulent kinetic energy
${\mathcal{K}}=(u^{\prime 2}+v^{\prime 2}+w^{\prime 2})/2$
, shown in figure 11(a), normalized by the friction velocity
$u_{\unicode[STIX]{x1D70F}}$
. As usual, the symbols represent the profiles from the direct numerical simulation of Kim et al. (Reference Kim, Moin and Moser1987) of a turbulent channel flow. Close to the wall and close to the centreline, all the profiles coincide. On the contrary, in the region where the maximum of
${\mathcal{K}}$
is located, i.e. the buffer layer, we observe a strong increase for
$Bi=1.4$
, and only a moderate one for the other values of Bingham number. Also, the peak is displaced to higher wall-normal distances
$y^{+}$
than its Newtonian counterparts. The increased value of the peak is mainly due to the increase of the streamwise component of the velocity fluctuations discussed above. An opposite behaviour is evident in figure 11(b), where the turbulent production
${\mathcal{P}}=-\overline{u^{\prime }v^{\prime }}\,\text{d}\overline{u}/\text{d}y$
is displayed. Indeed, although all the profiles of
${\mathcal{P}}$
still collapse at the wall and at the centreline, in the buffer layer the turbulent production decreases with the Bingham number, with differences noticeable in the viscous sublayer as well. Figure 12(a) shows the turbulent dissipation
$\unicode[STIX]{x1D700}=\unicode[STIX]{x1D707}\overline{\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}}$
of the fluctuating velocity field
$u_{i}^{\prime }$
. We observe that
$\unicode[STIX]{x1D700}$
has a maximum at the wall and then decreases moving towards the centreline where it reaches its minimum value which is the approximately same for the considered turbulent cases. The case with
$Bi=0$
shows a dissipation profile similar to the one of a Newtonian fluid (Rosti et al.
Reference Rosti, Cortelezzi and Quadrio2015), while increasing
$Bi$
the dissipation decreases monotonically, being null in the laminar case when
$Bi=2.8$
. The decrease of dissipation with
$Bi$
is consistent with the progressive decrease of
$Re_{\unicode[STIX]{x1D70F}}$
previously observed.
(a) Trace and (b) shear component of the mean elastoviscoplastic stress tensor
$\overline{\unicode[STIX]{x1D70F}}_{ij}$
as a function of the wall-normal distance
$y$
. The colour scheme is the same as in figure 6.
We now discuss in more details the elastoviscoplastic stress tensor
$\unicode[STIX]{x1D70F}_{ij}$
. Figure 13 shows the mean profile of the microstructure stress tensor trace
$\overline{\unicode[STIX]{x1D70F}}_{ii}$
(a) and the shear component
$\overline{\unicode[STIX]{x1D70F}}_{12}$
(b) as functions of the Bingham number. All the stress profiles have their maximum values at the wall (
$y=0$
) and their minimum absolute values at the centreline (
$y=h$
), with the trace being symmetric with respect to
$y=h$
and the shear component anti-symmetric. In the turbulent flows (
$Bi\lesssim 2$
), the normal stresses vary only slightly across the channel, except in the near-wall region where they rapidly grow. On the other hand, this trend is almost inverted for the laminar flow (
$Bi=2.8$
). A similar behaviour is shown by the shear stress component
$\overline{\unicode[STIX]{x1D70F}}_{12}$
, except that the almost uniform region is less wide, since in the middle of the channel the stress needs to vanish. Moreover, this region with an almost uniform elastoviscoplastic shear stress further reduces with the Bingham number
$Bi$
. We observe that the shear stress component and the trace of the stress tensor are approximately of the same magnitude, and that the values of the stress components approximately scale with
$Bi$
, but only when the flow is turbulent.
Normalized shear stress balance across the channel, for (a)
$Bi=0$
and (b)
$Bi=1.4$
. The dashed, dotted and dash-dotted lines are used for the viscous, Reynolds and elastoviscoplastic shear stress, respectively, while the solid line is the total shear stress, which varies linearly across the channel height.
To gain further understanding, we report in figure 14 the shear stress budget for the cases with
$Bi=0$
and
$Bi=1.4$
, normalized with the corresponding wall stress. For
$Bi=0$
, the additional elastoviscoplastic stress is very small and the behaviour is therefore similar to that of a standard Newtonian turbulent channel flow, with the viscous stress dominating at the wall and then rapidly decreasing towards the channel core; the Reynolds stresses, on the other hand, are zero at the wall and at the centreline and attain a maximum relatively close to the wall. Note that, although very small, the elastoviscoplastic stress is not null at the wall, thus the total wall shear stress is the sum of the two contributions, the elastoviscoplastic and viscous stresses. The situation differs in the flow at the Bingham number
$Bi=1.4$
. Here, the elastoviscoplastic stress increases across the whole channel, reaching approximately
$15\,\%$
of the total stress at the wall. Its increase is compensated by a changes of the other two stress components; in particular, the Reynolds stress peak reduces from
$70$
to
$55\,\%$
of the total, while the viscous stress peak from
$95$
to
$85\,\%$
. An overall picture of the total shear stress profile can be gained by studying the effective shear viscosity
$\unicode[STIX]{x1D707}_{e}$
normalized with the total viscosity
$\unicode[STIX]{x1D707}_{0}$
, reported in figure 12(b). This is defined as follows:
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D707}_{e}}{\unicode[STIX]{x1D707}_{0}}=\frac{\displaystyle \unicode[STIX]{x1D707}_{f}{\displaystyle \frac{\text{d}\overline{u}}{\text{d}y}}+\overline{\unicode[STIX]{x1D70F}}_{12}}{\displaystyle \unicode[STIX]{x1D707}_{0}{\displaystyle \frac{\text{d}\overline{u}}{\text{d}y}}}.\end{eqnarray}$$
The effective shear viscosity
$\unicode[STIX]{x1D707}_{e}$
grows with the Bingham number
$Bi$
and moving from the wall towards the centre of the channel. The inset of the figure shows the same quantity
$\unicode[STIX]{x1D707}_{e}$
as a function of the shear rate
$\dot{\unicode[STIX]{x1D6FE}}$
here defined as
$\unicode[STIX]{x1D707}_{0}\,\text{d}\overline{u}/\text{d}y$
, i.e. the denominator of (3.4); from the figure we can appreciate the shear-thinning behaviour of the elastoviscoplastic fluid described by the Saramito model (Saramito Reference Saramito2007).
Wall-normal profiles of the cross-correlations
$\unicode[STIX]{x1D70C}_{i}$
of the streamwise (dashed line), wall-normal (solid line) and spanwise (dash-dotted line) velocity component
$u_{i}$
and elastoviscoplastic contribution
$f_{i}$
, defined in (3.5), for (a)
$Bi=0$
and (b)
$Bi=1.4$
.
Finally, figure 15 shows the cross-correlation
$\unicode[STIX]{x1D70C}_{i}$
defined as
$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{i}={\displaystyle \frac{\overline{u_{i}^{\prime }\,f_{i}^{\prime }}}{u_{i}^{\prime }\,f_{i}^{\prime }}},\end{eqnarray}$$
where
$f_{i}$
is the elastoviscoplastic volume force, i.e. the contribution to the Navier–Stokes equation of the elastoviscoplastic stress tensor
$\unicode[STIX]{x1D70F}_{ij}$
defined as
$f_{i}=\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}/\unicode[STIX]{x2202}x_{j}$
. Note that there is no summation over the repeated indices in the previous relation. The left and right panels of the plot show the cases with
$Bi=0$
(a) and
$Bi=1.4$
(b), with the dashed, solid and dash-dotted lines in each plot corresponding to the streamwise, wall-normal and spanwise cross-correlations:
$\unicode[STIX]{x1D70C}_{1}$
,
$\unicode[STIX]{x1D70C}_{2}$
and
$\unicode[STIX]{x1D70C}_{3}$
. In general, when
$\unicode[STIX]{x1D70C}_{i}$
equals
$1$
or
$-1$
, the flow velocity and elastoviscoplastic force are perfectly correlated or anti-correlated, while when
$\unicode[STIX]{x1D70C}_{i}=0$
they are not correlated. We observe that for the case
$Bi=0$
all the cross-correlation components
$\unicode[STIX]{x1D70C}_{i}$
are negative in most of the channel, except in the near-wall region where
$\unicode[STIX]{x1D70C}_{1}$
equals
$1$
and
$\unicode[STIX]{x1D70C}_{2}$
equals
$0$
. In this case, the cross-correlation is almost uniform across the whole channel height, with a negative value equal to
$\unicode[STIX]{x1D70C}_{i}\approx -0.5$
: the elastoviscoplastic body force and velocity are anti-correlated in the bulk of the flow away from the walls, thus indicating that the elastoviscoplastic contribution is opposing the turbulent fluctuations. Close to the wall, however, the high positive values attained by
$\unicode[STIX]{x1D70C}_{1}$
suggest a role played by the viscoelastic stresses in the increase of the streamwise velocity fluctuations. Note that, this is similar to what was found by Dubief et al. (Reference Dubief, Terrapon, White, Shaqfeh, Moin and Lele2005) for a turbulent channel flow with a polymer solution. On the other hand, the high Bingham case (
$Bi=1.4$
) shows a similar trend only close to the wall (
$y\lesssim 0.5h$
), while in the bulk the cross-correlations
$\unicode[STIX]{x1D70C}_{i}$
interestingly go to zero. The fact that the flow velocity and the elastoviscoplastic stress tensor are not correlated around the centreline can be associated with the continuous cycle of yielding and unyielding processes.
(a,d,g,j,m) Isosurfaces and ((b,e,h,k,n) and (c,f,i,l,o)) contours of instantaneous streamwise velocity fluctuation
$u^{\prime }$
. The flow goes from left to right and the colour scale ranges from
$-0.25U_{b}$
(blue) to
$0.25U_{b}$
(red). The brown regions in all the figures represent the unyielded fluid. The two slices on the right are
$x{-}z$
planes at
$y=0.15h$
and
$y=0.44h$
. The Bingham number
$Bi$
increases from top to bottom (
$Bi=0$
,
$0.28$
,
$0.7$
,
$1.4$
and
$2.8$
), and the viscosity ratio
$\unicode[STIX]{x1D6FD}$
is fixed equal to
$0.95$
.
The elastoviscoplastic character of the flow affects the near-wall turbulent structures, and this is visually confirmed in figure 16. The left panels identify the low- (blue) and high-speed (red) near-wall streaks with isosurfaces of the streamwise velocity fluctuations,
$u^{\prime }$
, corresponding to the levels
$u^{\prime +}=\pm 0.25U_{b}$
, while the pictures in the middle and right columns show the footprints of these structures on the wall-parallel planes at
$y=0.15h$
and
$0.44h$
. It is evident that the structures in the buffer layer are less fragmented and more elongated in the streamwise direction as the Bingham number increases. Also, their spanwise extension increases, and consequently the number of streaks reduces. Indeed, the attenuation of the small-scale features is consistent with a picture where the larger coherent structures grow in size due to the reduction of the friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}$
, i.e. drag reduction. This effect – decreasing drag and wider and more coherent structures – is similar to what is found in other drag reducing flows, such as riblets, polymer suspensions and anisotropic porous walls, as already discussed previously. From the three-dimensional visualizations, we observe that the low-speed streaks penetrate to higher wall-normal distances than the high-speed ones; the former are usually associated with wall-normal velocity fluctuations
$v^{\prime }$
away from the wall, and their high-speed counterparts with wall-normal velocities towards the wall (Kim et al.
Reference Kim, Moin and Moser1987). This tendency interacts with the yield/unyield process; indeed, from the rightmost panels in the figure we note that the regions where the fluid is not yielded are mostly located in the positions above an high-speed streak, while almost all the fluid above the low-speed streaks remains fully yielded. These visual observations will be now quantified statistically by analysing the autocorrelation functions.
Stack of two-point velocity autocorrelation functions across the channel
${\mathcal{R}}_{ii}(y)$
. (a–d) Shows the streamwise autocorrelation function of the streamwise velocity
${\mathcal{R}}_{11}$
, while (e–h) the spanwise autocorrelation function of the wall-normal velocity
${\mathcal{R}}_{22}$
. The Bingham numbers increase from left to right (
$Bi=0$
,
$0.28$
,
$0.7$
and
$1.4$
). The solid and dashed lines correspond to positive and negative values of autocorrelation, ranging from
$-0.1$
to
$0.9$
with a step of
$0.2$
between two neighbouring lines. The colour scale ranges from
$-0.1$
(purple) to
$1.0$
(red).
The effect of the Bingham number
$Bi$
on the flow coherence is quantified by the two-point velocity autocorrelation functions, reported in figure 17 for the cases with
$Bi=0$
,
$0.28$
,
$0.7$
and
$1.4$
. The two-point autocorrelation function
${\mathcal{R}}_{ii}$
is defined here as
$$\begin{eqnarray}{\mathcal{R}}_{ii}(\boldsymbol{x},\boldsymbol{r})=\frac{\overline{u_{i}^{\prime }(\boldsymbol{x})u_{i}^{\prime }(\boldsymbol{x}+\boldsymbol{r})}}{\overline{u_{i}^{\prime 2}(\boldsymbol{x})}},\end{eqnarray}$$
where the bar denotes average over time and the two homogeneous directions, and the prime the velocity fluctuation. Note that, there is no summation over the repeated indices in the previous relation. The top row in the figure shows the distribution in the
$x{-}y$
plane of the streamwise velocity component autocorrelation along the streamwise direction
$x$
, while the bottom row the distribution in the
$z{-}y$
plane of the wall-normal velocity component autocorrelation along the spanwise direction
$z$
. In the case
$Bi=0$
(shown in the leftmost column), the correlations appear to be very similar to the baseline Newtonian case with highly elongated streaky structures that alternate at the canonical spanwise distance (i.e.
$\unicode[STIX]{x0394}z^{+}\simeq 100^{+}$
) (Kim et al.
Reference Kim, Moin and Moser1987). As the Bingham number
$Bi$
is increased (panels from left to right), the streamwise correlation length monotonically increases, thus indicating a higher level of coherency of the flow structures; in particular, this reveals that the velocity streaks are more elongated in the streamwise direction. The spanwise correlation also increases when increasing
$Bi$
as the velocity streaks become wider than in a Newtonian fluid, despite the existence of yielded regions in the channel core (see figure 16). Note also that, the increased correlation lengths in both the streamwise and spanwise directions are not limited to the near-wall regions occupied by the streaky structures, as in the other drag reducing flows cited above (Dubief et al.
Reference Dubief, White, Terrapon, Shaqfeh, Moin and Lele2004; García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011; Rosti & Brandt Reference Rosti and Brandt2018; Shahmardi et al.
Reference Shahmardi, Zade, Ardekani, Poole, Lundell, Rosti and Brandt2018), but extend up to the centreline. However, this increased spanwise correlation does not extend towards the centreline as much as the streamwise coherence. This difference originates from the fact that the flow at a high wall-normal distance
$y$
is still not yielded, thus behaving as a viscoelastic solid.
4 Conclusion
We present numerical simulations of laminar and turbulent channel flow of a non-Newtonian elastoviscoplastic fluid. The elastoviscoplastic flow is simulated with the model proposed by Saramito (Reference Saramito2007), where the incompressible Navier–Stokes equations are coupled with an additional equation for the evolution of the elastoviscoplastic stress tensor. In particular, the model predicts only recoverable Kelvin–Voigt viscoelastic deformation for stress below the yield stress value, while when the stress exceeds the yield value, the fluid behaves as an Oldroyd-B viscoelastic fluid. For both the laminar and turbulent regimes, we examine the flow behaviour when changing the values of the material plasticity (Bingham number) and viscosity ratio (
$\unicode[STIX]{x1D6FD}$
), while keeping the elasticity constant (Weissenberg number) to a small value, to more clearly identify the role of the plasticity.
For the laminar channel flow, we carried out a full parametric study and find that the friction factor increases with the Bingham number and decreases with the Reynolds number and the viscosity ratio. The drag increase due to the Bingham number originates from the increase of the portion of the channel where the stress is below the yield stress value and thus the fluid is not yielded. In these regions, forming initially at the centreline and then growing towards the wall as the Bingham number increases, the flow presents a flat velocity profile with zero shear. We propose an empirical correlation for the friction factor in a laminar channel flow, which is function of the Reynolds number, Bingham number and viscosity ratio. We show that the Fanning friction factor is inversely proportional to the Reynolds number and proportional to the square root of the Bingham number (a result interestingly found also for the flow of the same kind of fluid in a porous medium).
In the turbulent flows, the bulk Reynolds number is fixed to
$2800$
due to computational costs and the effects of different yield stress values and viscosity ratios are studied via both statistical data and instantaneous visualizations. Unlike the case of laminar flows, the Fanning friction factor is almost independent of the viscosity ratio and decreases with the Bingham number. We show that all the elastoviscoplastic flow configurations analysed are drag reducing, and since the Weissenberg number considered is very low, we have demonstrated that this is not an effect of the elasticity. We identify three different regimes depending on the value of the Bingham number: for low Bingham numbers, the turbulence is only slightly modified, except for a slowly progressive reduction of the friction Reynolds number. Next, for intermediate Bingham numbers, the flow becomes highly intermittent, with a continuous cycle of yielding and unyielding process in the centre of the channel which is responsible for the increased fluid oscillations. We also document strong streamwise velocity fluctuations, with the mean velocity profile departing from the usual log law and the loss of the inertial equilibrium range; all of the flow statistics are affected both in the buffer and logarithmic layers. Finally, for high values of the Bingham number, the flow fully laminarizes.
We show that the progressive increase of the amount of fluid which is not yielded with the Bingham number has a strong influence on the flow. Indeed, these regions grow from the centreline towards the walls as the Bingham number increases, similarly to the laminar regime, but introduce a strong unsteadiness in the flow when they extend over the full streamwise and spanwise dimensions.
In the elastoviscoplastic flow we observed an enhancement of the near-wall streak intensity and of the associated quasi-streamwise vortices, regions with localized high stress values. 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.
The analysis performed here assumed a very low level of elasticity of the flow. The present results can therefore be extended by introducing this additional effect and investigating how the dynamics described here changes. Furthermore, more complex flow configurations, e.g. separating and fully inhomogeneous flows, as well as the addition of a dispersed solid phase in this complex matrix deserve further consideration. Another interesting extension of the present work is the analysis of these flows at higher Reynolds numbers, investigating how the friction factor depends on the Bingham number and the absence of unyielded regions in the viscous sublayer.
Acknowledgements
M.E.R. and L.B. were supported by the European Research Council grant no. ERC-2013-CoG-616186, TRITOS and by the Swedish Research Council grant no. VR 2014-5001. D.I. and O.T. acknowledge financial support by the Swedish Research Council through grants no. VR2013-5789 and no. VR 2014-5001. S.H. acknowledges financial support by NSF (grant no. CBET-1554044-CAREER), NSF-ERC (grant no. CBET-1641152 Supplementary CAREER). The authors acknowledge computer time provided by SNIC (Swedish National Infrastructure for Computing).
