# [

###### Abstract

We present numerical simulations of laminar and turbulent channel flow of an elastoviscoplastic fluid. The non-Newtonian flow is simulated by solving the full incompressible Navier-Stokes equations coupled with the evolution equation for the elastoviscoplastic stress tensor. The laminar simulations are carried out for a wide range of Reynolds numbers, Bingham numbers and ratios of the fluid and total viscosity, while the turbulent flow simulations are performed at a fixed bulk Reynolds number equal to and weak elasticity. We show that in the laminar flow regime the friction factor increases monotonically with the Bingham number (yield stress) and decreases with the viscosity ratio, while in the turbulent regime the the friction factor is almost independent of the viscosity ratio and decreases with the Bingham number, until the flow eventually returns to a fully laminar condition for large enough yield stresses. Three main regimes are found in the turbulent case, depending on the Bingham number: for low values, the friction Reynolds number and the turbulent flow statistics only slightly differ from those of a Newtonian fluid; for intermediate values of the Bingham number, the fluctuations increase and the inertial equilibrium range is lost. Finally, for higher values the flow completely laminarises. These different behaviors are associated with a progressive increases of the volume where the fluid is not yielded, growing from the centerline towards the walls as the Bingham number increases. The unyielded region interacts with the near-wall structures, forming preferentially above the high speed streaks. In particular, the near-wall streaks and the associated quasi-streamwise vortices are strongly enhanced in an highly elastoviscoplastic fluid and the flow becomes more correlated in the streamwise direction.

[figure]style=plain,subcapbesideposition=top
Turbulent channel flow of an elastoviscoplastic fluid]
Turbulent channel flow of an elastoviscoplastic fluid
M. E. Rosti, D. Izbassarov, O. Tammisola, S. Hormozi and L. Brandt]Marco E. Rosti^{†}^{†}thanks: Email address for correspondence: , Daulet Izbassarov, Outi Tammisola,

Sarah Hormozi and Luca Brandt

## 1 Introduction

Many fluids in nature and industrial applications exhibit a non-Newtonian behavior, i.e., a non-linear relation between the shear stress and the shear rate, such as shear thinning, shear thickening, yield stress, thixotropic, shear banding and viscoelastic behaviors. 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 , 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, 1963, 1967; Maleki & Hormozi, 2017).

### 1.1 Stability of yield stress fluids

Several studies have been devoted to the stability of yield stress fluids (Nouar & Frigaard, 2001; Metivier et al., 2005; Nouar et al., 2007; Nouar & Bottaro, 2010; Bentrad et al., 2017). The first study on the stability of viscoplastic fluid flows was reported by Frigaard et al. (1994), who studied the linear stability of a Bingham fluid in a plane channel flow. More recently, Nouar et al. (2007) 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. (2005) 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 (2001). 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 et al. (2010) has performed this type of analysis for a Rayleigh-Benard-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 geometry exist: recently, nonlinear (energy) stability analysis has been extended to multi-layer flows of yield stress and viscoelastic fluids by Moyers-Gonzalez et al. (2004); Hormozi & Frigaard (2012). Recently, in order to identify possible paths to transition Nouar & Bottaro (2010) perturbed the base flow slightly, and found that very weak defects are indeed capable to excite exponentially amplified streamwise traveling waves. Finally, Kanaris et al. (2015) 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, 2001), and for different properties of the wall surface, e.g., roughness (Orlandi & Leonardi, 2008), porosity (Breugem et al., 2006; Rosti et al., 2015, 2018b) and elasticity (Rosti & Brandt, 2017). 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 (Hanks, 1963, 1967; Hanks & Dadia, 1971; Ryan & Johnson, 1959), with the most popular phenomenological approach suggested by Metzner & Reed (1955). 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 (1990); Founargiotakis et al. (2008). Rudman et al. (2004) 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, 1972; Pinho & Whitelaw, 1990; Escudier & Presti, 1996; Den Toonder et al., 1997; Beris & Dimitropoulos, 1999; Warholic et al., 1999; Escudier et al., 1999; Escudier & Smith, 2001; Dubief et al., 2004, 2005; Escudier et al., 2005, 2009; Xi & Graham, 2010; Owolabi et al., 2017; Shahmardi et al., 2018). The interested reader is refereed to the work by Berman (1978) and White & Mungal (2008) for a through review on the subject.

### 1.3 Elastoviscoplastic fluid

Despite the numerous studies performed to analyze 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 (2006) and Guang et al. (2011). 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 over-predict 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 behavior. Thus, in order to properly predict the behavior 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 (2007) proposed a new constitutive equation for elastoviscoplastic fluid flows, which reproduces a viscoelastic solid for stresses lower then the yield stress, and a viscoelastic Oldroyd-B fluid for stresses higher then 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., 2012; Martinie et al., 2013). Cheddadi et al. (2011) simulated the inertialess flow of an elastoviscoplastic fluid around a circular object using the model proposed by Saramito (2007); 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 (2007) for the flow of a liquid foam and by Putz et al. (2008) 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 et al. (2016) using the constitutive law by Saramito (2007); their results are in a quantitative agreement with experimental observations obtained with a Carbopol gel.

The model proposed by Saramito (2007) was extended by the same author to account for shear-thinning effects (Saramito, 2009). The new model combines the Oldroyd viscoelastic model with the Herschel-Bulkley viscoplastic model, with a power law index that allows a shear-thinning behavior 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 others exist in the literature. The interested reader is referred to Crochet & Walters (1983); Balmforth et al. (2014); Saramito & Wachs (2016); Saramito (2016) for a through 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 . The non-Newtonian flow is simulated by solving the full unsteady incompressible Navier-Stokes equations coupled with the model proposed by Saramito (2007) for the evolution of the additional elastoviscoplastic stress tensor. In section 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 section 2.2, while the results on the laminar and on the fully developed turbulent channel flows are presented in section 3. In particular, we discuss the role of some of the parameters defining the elastoviscoplastic fluid, i.e., the Bingham number and the viscosity ratio . Finally, a summary of the main findings and conclusions are presented in section 4.

## 2 Formulation

[] \sidesubfloat[]