Turbulent Duct Flow with Polymers

(1)

IN

DEGREE PROJECT MECHANICAL ENGINEERING, SECOND CYCLE, 30 CREDITS

STOCKHOLM SWEDEN 2018,

Turbulent Duct Flow with Polymers

ARMIN SHAHMARDI

KTH ROYAL INSTITUTE OF TECHNOLOGY SCHOOL OF ENGINEERING SCIENCES

(2)

www.kth.se

(3)

Turbulent Duct Flow with Polymers

Author: Armin Shahmardi Supervisor: Luca Brandt Examiner: Luca Brandt

Master’s degree project

Stochholm, Sweden, January 2018 Department of Mechanical Engineering KTH Royal Institute of Technology

(4)

Abstract

Direct numerical simulation of the turbulent ow in a square duct with polymers has been performed and the results have been compared with the laboratory experiments done at KTH Mechanical engineering department. The polymer suspension is simulated with the FENE-P model, and the numerical results are used to elucidate the mechanism that provides drag reduction and the eect of polymers on the secondary motion which is typical of the turbulent

ow in ducts. Experiments are used to support and validate the numerical data, and to discuss the Reynolds number dependency of the obtained drag reduction.

The study shows that the Prandtl's secondary ow is modied by the polymers:

the classical 8 regions, in the cross section with high vorticity, are bigger in the polymer ow than those in the Newtonian case, and their centers are displaced towards the center of the duct away from the wall. In plane uctuations are reduced and streamwise coherence of the ow enhanced in the presence of polymers.

(5)

Sammanfattning

Direkt numerisk simulering av det turbulenta ödet i en kvadratisk kanal med polymerer har utförts och resultaten har jämförts med de laboratorieexperiment som gjorts vid KTH: s maskintekniska avdelning. Polymersuspensionen simule- ras med FENE-P-modellen och de numeriska resultaten används för att belysa mekanismen som ger dragreduktion och eekten av polymerer på sekundärrö- relsen som är typisk för det turbulenta ödet i kanaler. Experiment används för att stödja och validera de numeriska data och för att diskutera Reynolds bero- endet av den erhållna dragreduceringen. Studien visar att Prandtls sekundära

öde modieras av polymererna: de klassiska 8 regionerna i tvärsnittet med hög vorticitet är större i polymerödet än de i det newtonska fallet och deras centra är förskjutna mot centrum av kanalen bort från väggen. I planuktuationer re- duceras och strömningsförstärkt sammanhängning av ödet förbättras i närvaro av polymerer.

(6)

Introduction

1.1 Aim and objectives

Near-wall turbulence is responsible for significant drag penalties in many flows of engineering relevance, and because of that, many researchers are studying various ways to be able to properly control the flow (Choi et al., 1993; Dubief et al., 2004; Breugem et al., 2006; Orlandi &

Leonardi, 2008; García-Mayoral & Jiménez, 2011; Rosti et al., 2015). Among the many control strategies, the use of polymers has been demonstrated to be very efficient to reduce drag in pipelines (Virk, 1971). While flows through axisymmetric geometry have been extensively studied (Virk et al. , 1967; Cho & Harnett, 1982; Toonder et al., 1997; Escudier et al., 1999;

Ptasinski et al., 2003; Resende et al., 2006), less attention has been given to more complex geometries, such as square ducts. In this context, the aim of this work is to explore and bet- ter understand the interactions between the turbulent flow in ducts with square section and polymer additives.

1.2 Duct flow

Turbulent flow in a plane channel has been extensively studied in the past (Kim et al., 1987) while the flow in a duct with a square cross-section has received much less attention (Uhlmann et al., 2007; Pinelli et al., 2010; Samanta et al., 2015; Vinuesa et al., 2015; Owolabi et al., 2016), although its geometry is only mildly more complex. Peculiar features of the turbulent flow in a duct are that the mean velocity profile presents a non-uniform distribution of the skin friction coefficient along the edges, and the related appearance in the cross-sectional plane of secondary motions of the second kind, as classified by Prandtl (1926), which is a mean flow effect induced by gradients of turbulence fluctuations, e.g., a breaking of axisymme- try. The existence of such secondary mean motion in this geometrical configuration is well known since the experiments by Nikuradse (1926), which was followed by further experimental measurements (Brundrett & Baines, 1964; Gessner, 1973; Melling & Whitelaw, 1976) as well as direct and large eddy simulations (Madabhushi & Vanka, 1991; Gavrilakis, 1992;

Uhlmann et al., 2007; Pinelli et al., 2010) all of which subsequently extended our knowledge of this flow. Uhlmann et al. (2007) simulated the case of the marginal Reynolds number regime and showed that the buffer layer coherent structures play a crucial role in the appearance of the secondary flow and in the deformation of the mean streamwise velocity profile. In- deed, they proposed that the deformation of the mean streamwise velocity profile is due to the presence of preferential positions of quasi-streamwise vortices and velocity streaks.

1

(9)

2 CHAPTER 1. INTRODUCTION

Moreover, in such a marginal regime, short-time averaged velocity fields are found to ex- hibit a 4-vortex state instead of the usual 8-vortex secondary flow pattern found at higher Reynolds numbers. This feature was explained with the relation between coherent structures and secondary flow: if the dimension of the cross-section, in wall units, is below the one needed to accommodate a complete minimal turbulent cycle on all four walls (Jimenez

& Moin, 1991), then just two facing walls can alternatively give rise to a complete turbulent regeneration mechanism, while the other two faces remain in a relative quiescent state.

Such 4-vortex states at marginal Reynolds number have also been observed experimentally ( Owolabi et al., 2016). Pinelli et al. (2010) studied turbulent duct flows at higher Reynolds numbers, where the length of the edge of the square cross-section, expressed in viscous units, is larger, therefore allowing for the simultaneous presence of multiple pairs of high- and low- velocity streaks. These authors further proved the idea that the secondary flow is a footprint of the coherent motions embedded in the turbulent flow, and proposed a scenario in which the position of the center of the mean secondary vortices is determined by the preferential positioning of the quasi-streamwise vortices.

1.3 Polymer additives

Polymer addition is a very efficient strategy employed for drag reduction in wall-bounded turbulent flows, since drag reduction (DR) up to 80% have been achieved with only few ppm concentrations. Two distinct regimes are usually identified (Warholic et al., 1999), often called low and high drag reduction, LDR and HDR, respectively. The former (LDR) exhibits a log-law region of the mean velocity profile parallel to that of the Newtonian flow but with an upward shift associated with the amount of drag reduction; streamwise fluctuations are increased while wall-normal and spanwise components reduced, as well as the shear stress.

The regime characterised by high drag reductions (HDR), more than 40%, shows a strong increase in the slope of the log-law region and low levels of Reynolds shear stress (Ptasinski et al., 2003). The drag reduction is eventually bounded by a maximum (MDR) (Virk et al., 1970). Several works (Warholic et al., 1999; Ptasinski et al., 2003) observed low Reynolds shear stress and a consequent deficit in the stress balance equation, which has been interpreted as the input of energy from the polymers to the flow which ultimately sustains the asymptotic MDR.

The drag reduction mechanism is quite complex due to its multiscale nature (small amount of microscopic polymer molecules is needed to achieve significant drag reduction in the bulk flow) and its full physical understanding is still incomplete. Several explanations have been proposed to gain further insight in the mechanism of polymer drag reduction: Dimitropou- los et al. (2001) showed that streamwise enstrophy is inhibited by the extensional viscosity generated by polymers stretched by turbulence; Ptasinski et al. (2003) found a shear shelter- ing effect in the near-wall region which decouples the outer and inner layer vortices; Min, Yoo & Choi (2003a) and Min et al. (2003b) observed the transport of elastic energy from the viscous sublayer to the buffer and the log region, and proposed a criterion for the onset of drag reduction.

Dubief et al. (2004, 2005) studied the intermittency of polymers in turbulent flows, showing that the action of polymers is as intermittent as the near-wall vortices, and that the drag- reducing property of polymers is closely related to coherent turbulent structures. Polymers dampen near-wall vortices but also enhance streamwise kinetic energy in near-wall streaks.

The net balance of these two opposite actions leads to a self-sustained drag-reduced turbu-

(10)

CHAPTER 1. INTRODUCTION 3

lent flow: the polymers reduce turbulence by opposing the downwash and upwash flows generated by near-wall vortices, while they enhance streamwise velocity fluctuations in the very near-wall region. Recent studies by Xi & Graham (2010, 2012a,b) provided new insight on the mechanism by which polymer additives reduce the drag. These authors suggested that a turbulent flow is characterised by an alternate succession of strong and weak turbulence phases. The first are characterised by flow structures showing strong vortices and wavy streaks, the latter weak streamwise vortices and almost streamwise-invariant streaks. In the Newtonian flow, the so-called active turbulence dominates; with increasing viscoelasticity, on the contrary, active intervals becomes shorter while the so-called hibernating intervals are unaffected. Also, it is shown that during these hibernating turbulence intervals, the turbulent dynamics resemble MDR turbulence in both Newtonian and viscoelastic flows (Li et al., 2006; White et al., 2004).

Polymers can also alter flow instabilities and transition to turbulence. Indeed, Biancofiore et al. (2017) recently examined the secondary instability of streaks in a viscoelastic flow, 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.

1.4 Polymer solutions in square duct flows

Extensive work has been done experimentally to try to characterise polymer solution be- haviour in duct flows. In particular, Rudd (1972) and Logan (1972) reported limited measurements of the mean-flow and turbulence structure for flow of drag-reducing polymers through square tubes, with sufficiently low polymer concentrations, and provided no infor- mation about the secondary flow. More recently, Gampert & Rensch (1996) have reported the results of a systematic study of the influence of polymer concentration on near-wall turbulence structure for the flow through a square duct, but again have provided no direct infor- mation on the secondary flow. Also, they argued that there are two flow regimes in which the properties of the polymer structure, and hence the turbulent flow-field, differ significantly, depending on the polymer concentration. Escudier & Smith (2001) presented both global data (friction-factor versus Reynolds number) and detailed mean axial- and secondary-flow- velocity and turbulence- field data for fully-developed flow through a square duct of various polymer solutions. They found a reduction in the intensity of turbulent-velocity fluctuations transverse to the mean flow and a strong reduction in the secondary-flow velocities. More recently, Escudier et al. (2009) provided a comprehensive experimental database of previous experiental works ( Rudd , 1972; Gampert & Yong , 1990; Escudier & Smith , 2001); the authors pointed out that several of these studies have limited turbulence data, and considers mainly relatively low polymer concentrations. Escudier et al. (2009) gave special emphasis on the quantification of turbulence anisotropy and showed that with polymer additives, the flow displays a tendency towards the axisymmetric-turbulence limit. Also, there is a marked decrease in anisotropy with distance from the near-surface peak in all cases but this tendency progressively reduces with increasing concentration/drag-reduction level.

1.5 Outline

In this work, Direct Numerical Simulations (DNS) of a turbulent duct flow with polymer additives at moderate Reynolds number is presented. The viscoelastic fluid is modelled by the

(11)

4 CHAPTER 1. INTRODUCTION

constitutive FENE-P (Finite Extensible Nonlinear Elasticity-Peterlin) closure. In chapter 2, first, the flow configuration and governing equations is discussed, and then the implemented numerical methodology is presented. The results of simulation have been also compared to experimental study perfomed by Sagar Zade at KTH. The experimental apparatus and the techniques used to study the problem at hand are explained in the next section. The features of the turbulent duct flow are documented in chapter 6, where the results with and without the polymer additives have been compared, to elucidate their effect on the secondary motion and on the flow in general at this moderate Reynolds number. First and second order statistics are analysed, together with energy and vorticity budgets. Finally, a summary of the main findings and some conclusions are drawn in chapter 7.

(12)

Chapter 2

Mathematical Formulation

Studying non-Newtonian viscoelastic fluids such as polymer solutions, liquid crystal solutions, fiber suspensions, etc., has been of great interest during the last decades. Among a lot of strategies to reduce friction drag, polymer additive is the most efficient (Dubief et al., 2005). The stress endured by the viscoelastic element is governed by the qualitative and quantitative deformation that the element has been experienced. As a result, a coupling between the molecular configuration and rheological response is inevitable. Therefore, to study these kind of flows, a multiscale method is necessary to deal with both microscopic scale of kinetic theory and macroscopic scale of continuum mechanics (R.Keunings , 1997; Du et al., 2006). One of the most efficient multiscale models is FENE( Finite Extensible Nonlinear elasticity) which describes the dilute solution of a polymer by considering a single dumbbell (two beads connected by a elastic spring). In this model, the configuration variable (which is introduced later) represents the vector connects two mentioned beads.

2.1 FENE model

Using the FENE model, the solution of a polymer (now considered as a dumbbell) is achieved by solving a set of coupled PDEs involving the incompressible Navier-Stokes equation (macroscopic scale) and the Fokker- Plank equation (microscopic scale). The Fokker- Plank equation describes the Probability Density Function (PDF) of the dumbbell orientation (in microscopic level) (R.Keunings , 1997; Du et al., 2006):

∂~u

∂t +

~ u. ~∇

~

u + ~∇P = ~∇.~τc+ ν∆~u (2.1)

∂f

∂t +

~ u. ~∇

f + ∇_Q_~.

∇~u ~Qf

= 2 ξ∇_Q_~.

∇_Q_~ψ( ~Q)f

+2kT

ξ ∆_Q_~f (2.2) Where ~u is the velocity vector, p is the hydrostatic pressure, ν is the kinematic viscosity, τp is polymer contribution to stress, ~Q is polymer configuration variable, ξ is the friction coefficient of the dumbbell beads, T is the temperature,k is the Boltzmann constant and ψ( ~Q) is the elastic spring potential.

The second term in the left hand side of equation (2.2) corresponds to polymer convection (which needs especial consideration while the equation is discretized) and the third term represents polymer stretching by the macroscopic flow. The first right hand side term models the inner force of the dumbbell according to spring potential, and the last term in the right

5

(13)

6 CHAPTER 2. MATHEMATICAL FORMULATION

hand side (the diffusion term) illustrates the random collisions of the solvent particle with the polymer. Polymer contribution to stress is represented as follow (τp can be obtained by the least action principle):

τp= λ Z

(∇_Q_~ψ( ~Q)O ~Q)f (~x, ~Q, ~t) (2.3)

Where f is Fokker- Plank PDF and λ is the polymer density constant.

The Fokker- Plank equation with SDE reads (R.Keunings , 1997; Du et al., 2006):

d ~Q + ~u. ~∇ ~Qdt =

∇~u ~Q − 2 ξ∇_Q_~ψ

dt +

s 4kT

ξ d ~Wt (2.4)

Where ~Wt is the standard Brownian motion and acts as a noise term. To derive a simple constitutive equation for polymer stress from Fokker- Plank, the Hookean law is used as the main idea. The spring potential using the Hookean’s law is presented in the following equation:

ψ(Q) = HQ²

2 (2.5)

Where H is the elasticity constant and Q = | ~Q|. The FENE model includes a constraint for the extensibility, therefore there is no exact constitutive equation for polymer in this model.

According to the FENE model, spring force law reads:

∇_Q_~ψ = H ~Q

1 − (^Q_L)² (2.6)

Where L is the maximum dumbbell extension. Using the FENE model, the system of equation would be more complex and need more consideration in the case of boundary conditions and taking care of singularity of coefficients near the boundary of configuration domain. For simplicity, and without loss of accuracy, Du et al. (2006) assumed that the convective term (~u. ~∇)f and ∇~u as constant matrices.

Finally, a system of Navier-Stokes equation and the constitutive equation for polymer should be solved. To solve the coupled system of equation, the configuration space of ~Qshould be resolved or sampled at each point in space which is a difficult task. Several strategies have been introduced to tackle the SDE corresponds to Fokker-Plank equation, among which the Mont Carlo sampling is one the most efficient. Du et al. (2006) also devised the so-called moment closure approach to tackle Fokker-Plank equation.

By multiplying the two dimensional Fokker- Plank equation by < ~QNQ >~ and integrating on the probability space, and introducing configuration tensor as C =< ~QNQ >~ , where the brackets are the ensemble average over the space and is defined as follow:

< ~QO ~Q >=

Z Z

|Q|<Q0

Q~O ~Qf (~x, ~Q, t) d ~Q (2.7)

The Fokker- Plank equation reads:

∂C

∂t + (~u. ~∇) − (~∇~u)C − C( ~∇~u)^T = −4 ξ

Z H ~Q

1 − (^Q_L)²f (~x, ~Q, t) d ~Q + 4KT

ξ (2.8)

(14)

CHAPTER 2. MATHEMATICAL FORMULATION 7

To close the system, and reduce the computational costs of resolving the probability space, it is necessary to approximate the ensemble average term in the right hand side of the last equation, < ^{H ~}^Q^N^Q^~

1−(^Q

Q0)² >

2.2 FENE-P model

Various models have already been presented for an accurate approximation which agrees quantitatively with the original FENE equation. One of the best closure models is FENE-P model. FENE- P model is based on the following closure assumption(R.Keunings , 1997; Du et al., 2006):

< H ~QNQ~ 1 − (_Q^Q

0)² >= H < ~QNQ >~ 1 − (_Q^Q

0)² = HC

1 −^tr(C)_L2

(2.9) FENE-P model quantitatively differs from the original FENE models (even in steady cases).

An approach to tackle this problem is to clarify the class of equations which are suitable forms of PDF solutions to Fokker- Plank equation. As a result, the more the approximation for PDF is capable to resolve actual solution, the less unphysical solutions for the FENE-P model.

The final proposed equations for FENE-P model which is used in the present study are (Du- bief et al., 2005, 2004):

∂C_ij

∂t + uk

∂C_ij

∂x_k = Ckj

∂u_i

∂x_k + Cik

∂u_j

∂x_k − τ_ij (2.10)

τij = 1 W e

Cij

1 −^C_L^kk2

− δ_ij

!

(2.11) Once again, the first term on the right hand side represents the polymer stretching re- garding to the straining actions of the flow, τij is the stress tensor through which the polymer extract or release energy, and W e is the so-called Weissenberg number and is the ratio between the polymer relaxation time (λ ) and the time scale of the flow (tf).

W e = λ tf

(2.12)

tf = h Uc

(2.13) Where h is the reference length (half height of channel or duct) and Ucis the reference velocity (centerline velocity). Sometimes the Weissenberg number is defined based on the time scale of shear layer near the wall ( tf = ν/u²_τ).

W e = λu²_τ

ν (2.14)

Where uτ =q_τ

w

ρ is the friction velocity.

(15)

8 CHAPTER 2. MATHEMATICAL FORMULATION

2.3 Navier-Stokes equation

The dynamics of the incompressible fluid motion are governed by The Navier-Stokes equation together with the continuity equation:

∂

∂x_i (ρui) = 0 (2.15)

∂ui

∂t + uj

∂ui

∂x_j = −∂P

∂x_i +∂τ_ij^f

∂x_j + fi (2.16)

where uiare the velocity vectors, P is the hydrodynamic pressure, τ_ij^f is the fluid shear stress tensor and fi represents body forces, which hereafter will be neglected.

Genearlly, using the FENE model, the solution of a polymer (now considered as a dumbbell) is achieved by solving a set of coupled PDEs involving the incompressible Navier- Stokes equation (macroscopic scale) and the Fokker- Plank equation(microscopic scale) (R.Keunings , 1997; Du et al., 2006).

The modified Navier-Stokes equation, including the polymer effects can be written as:

∂u_i

∂t + u_j∂u_i

∂x_j = −∂P

∂x_i + β Re

∂

∂x_j

∂u_i

∂x_j

+1 − β Re

∂τij

∂x_j (2.17)

Here, β is the ratio between fluid viscosity and the solution viscosity. By solving the polymer’s evolution equation, equation 2.10, together with the Navier-Stokes equation, (equation 2.17) the solution will be achieved.

(16)

Chapter 3

Numerical Implementation

In this chapter, the numerical method which is used to solve the system of equations governs the dynamic of the flow of a polymer solution is presented.

Recent studies show that the numerical solution of viscoelastic fluids, regardless to the choice of the constitutive equation, discretization scheme, and time integration method encoun- tered with numerical breakdowns in the case of high Weissenberg numbers. The difficulties in these simulations are the main reason for the absence of the exact reason for the flow dynamics and the mentioned numerical breakdown. However, one possible reason is that since the constitutive equation does not include any diffusion term, any disturbance in numerical scheme amplifies and lead the code to diverge.

One of the earliest solution to this problem was introducing an artificial diffusivity(AD) to the constitutive equations or using discretization methods including dissipation errors.

Adding the global AD smears out the gradients of polymer stresses and the obtained flow field is completely different from the physical one. As a modification, global AD was replaced by local AD which has been using in many researches and provide physical results.

Spectral methods also have been used to solve the laminar flow of viscoelastic fluids. Al- though the results were successful in laminar flow, in the case of Direct Numerical Simula- tion of turbulent flow, once again the numerical breakdown occurred. To solve the problem, an AD was added locally and the physical solution was attained (Dubief et al., 2005).

Despite the fact that adding AD or using discretization method with dissipative errors mit- igate the mentioned problem, the polymer stress tensor is still smeared and the solution would deviate from the physical solution. The oscillations in the Cij tensor in high Weis- senberg number, which results in numerical breakdown can be explained by considering different terms of constitutive equation.

In equation 2.10, Ckj∂ui/∂x_k+ C_ik∂uj/∂x_k and −τij are the stretching and relaxation terms, respectively. Both the terms are local terms with respect to Cij. Therefore, no spatial diffusion is expectable in these terms. On the other hand, the advection term uk∂Cij/∂xkin the left hand side, includes spatial discretization. one can write the equation in the general form of an advection equation with a source term.

∂Cij

∂t + a∂Cij

∂x_k = f (Cij) (3.1)

Here, a is a constant with respect to conformation tensor. The source term consists of the stretching and relaxation terms. Since the stretching term is made up of multiplication of

9

(17)

10 CHAPTER 3. NUMERICAL IMPLEMENTATION

the velocity gradients and the conformation matrix, stretching occurs at all the scales, from the Kolmogrov scale to the scale of the largest energetic eddies. This means that in order for polymers to be relaxed and the drag reduction to occur, the relaxation time should be larger than Kolmogrov time scale. Therefore, according to the definition of Weissenberg number, it should be larger than or equal to one.

The Kolmogorov length η, velocity vη, and time tηscales are:

η = ν³

¹₄

, vη = (ν)¹⁴ , tη =

ν

¹

4 (3.2)

Where is dissipation of turbulent energy and ν is the kinematic viscosity. The Kol- mogrov time scale in the viscous sub layer is of order of ν/uτ and in the buffer layer is of order of ν/uτ2. Recalling the definition of Weissenburg number in both region, the last two kolmogrov scales corresponds to W e = 1. But even in the simulation of large Weissenburg number, the time scale should be larger or equal to Kolmogorov scales. In the absence of diffusion terms in the constitutive equation, this results in some sharp gradients in the flow including shock like solutions (according to the advection equation). Similar condition could happen for a passive scalar with low diffusivity. In this case, the smallest scale is Batchelor scale (Dubief et al., 2005):

Θ = √ηk

Sc (3.3)

Where Sc is Schmidt number and is defined as the ration of kinematic viscosity to molecular diffusion. It was shown that the smallest scale of Cij is of order of Batchelor scale.

A common approach to solve the advection equation is upwind scheme. Although this scheme reduces the risk of oscillations, it also increase the artificial dissipation.

To solve the problem, in this study, global artificial diffusivity has been considered in early steps of simulation. By approaching to the statistically steady state solution, the AD was gradually reduced so that it was totally removed at the the final steps. Moreover, a fifth order WENO scheme with dissipative error has also been implemented.

Recall the governing equations from the previous chapter, equations 2.10, 2.11, 2.15 and 2.16.

The integration method which is used is based on fractional-steps method. For the time advancement, two different scheme were used. For the velocity diffusion term, ∂²u_i/∂x_j∂x_j in the Navier-Stokes equation, the Crank-Nicolson method is used and for all the other terms of the Navier-Stokes equation and all the terms of constitutive equation, 2.10, the third order Runga-Kutta method has been implemented. In this chapter, the Fractional-Steps method is explained first. The Crank-Nicolson scheme is addressed in the next section. Details of the third order Runga-Kutta approach is elaborated the next part. Finally, the discretized form of the equations and the solution procedure is described.

3.1 Numerical schemes for incompressible Navier-Stokes

Considering incompressible form of the Navier-Stokes equations, there is no time derivative or pressure term in the continuity equation. As a result the evolution of this equation in time is not coupled with that of momentum equations. This typical problem of incompressible flows is usually solved through two different families of approaches, namely, the pressure- based and the density-based approaches. In this study, the Fractional Step method, from the pressure-based approach family has been implemented.

(18)

CHAPTER 3. NUMERICAL IMPLEMENTATION 11

Fractional step method (also called projection-Correction method) was first introduced by Chorin (Chorin , 1968). The basic idea comes from the Helmoltz-Hedge decomposition the- orem. This theory says that considering Ω a region in space with smooth boundary, all the vectors wiin Ω can be decomposed as wi = ui+_∂x^∂P

i such that ^∂u_∂xⁱ

i = 0(uiis divergence free).

Considering the system of governing equations of incompressible flows in the following matrix form:

uⁿ⁺¹_i − uⁿ_i

∆t + N (uⁿ)uⁿ+ G(Pⁿ⁺¹) = f (tⁿ) D(Uⁿ⁺¹) = 0

(3.4) Where N (uⁿ)uⁿrepresents the convective and diffusive terms, G is the gradient operator, f is body forces, and D is the divergence operator.

The idea is to first make a prediction of uⁿ⁺¹ , which is not necessarily divergence free, u^∗, then correct this prediction such that uⁿ⁺¹satisfies the divergence free condition.

u^∗_i − uⁿ_i

∆t + N (uⁿ)uⁿ= f (tⁿ) uⁿ⁺¹_i − u^∗_i

∆t + G(Pⁿ⁺¹) = f (tⁿ D(Uⁿ⁺¹) = 0

(3.5)

The firs step in equation 3.5 is called prediction step and the second one is the correction step.

Adding the Prediction and the Correction step, one can get the original matrix form of the equations. Taking a divergence of the correction step and knowing that D(Uⁿ⁺¹) = 0:

− 1

∆tDu^∗+ DG(Pⁿ⁺¹) = 0 (3.6)

Where DG stands for divergence of gradients or the Laplacian:

∇²(Pⁿ⁺¹) = − 1

∆tDu^∗ (3.7)

Which is the Poisson equation for pressure. It can be shown that the boundary condition for this equation is the simple Neuman boundary condition (∂P/∂n = 0).

To solve the system, first the prediction step should be solved to find u^∗. Knowing u^∗ one can solve the Poisson equation to get Pⁿ⁺¹. Finally, by solving the Correction equation, the divergence free solution is obtained.

3.2 Time advancement

In order to do a temporal discretization, different schemes have been introduced. Consider the semi- discretized form of an arbitrary equation:

dU

dt = LU (3.8)

Where L is the discretized operator of spatial derivatives. The simplest approaches for the time discretization are the first order forward and backward Euler schemes:

Uⁿ⁺¹− Uⁿ

∆t = L(Uⁿ)(Explicit or Forward Euler) Uⁿ⁺¹− Uⁿ

∆t = L(Uⁿ⁺¹)(Implicit or Backward Euler)

(3.9)

(19)

Both the schemes are of first order accuracy. The explicit schemes are usually less expensive and easier to implement while the implicit schemes are more stable and therefore are much more flexible in the case of selecting time steps.

There are other schemes which are more accurate, among which the Crank-Nicolson and Runga-Kutta scheme have been used in the present study and are introduced in this chapter 3.2.1 Crank-Nicolson Scheme

This scheme is a combination of the implicit and explicit form of the Euler equation and is a second order approach.

Uⁿ⁺¹− Uⁿ

∆t = 1

2 L(Uⁿ) + L(Uⁿ⁺¹) (3.10)

Stability analysis shows that this scheme is unconditionally stable (Hoffmann & Chiang, 2000).

3.2.2 Runge-Kutta Scheme

The main idea of this method is to increase the order of accuracy of time integration by using weighted average of solutions at different sub-steps over an time interval. Considering the above-mentioned semi discretized equation again, equation 3.8. As the simplest form, one can divide the time step into two sub- steps. The solution at this step is defined as u⁽²⁾where the superscription designates the solution at the time level of α∆t and in general 0 < α < 1.

To start the solution for the new time step (n+1), the solution at the previous step usually is set to the solution at the first sub- step (Hoffmann & Chiang, 2000). :

U⁽¹⁾ = U⁽ⁿ⁾ U⁽²⁾= U⁽ⁿ⁾−∆t

2 L(Uⁿ) (3.11)

To discretize the spatial derivatives, any finite difference approach can be used. The final solution for uⁿ⁺¹can be computed as follow:

U⁽ⁿ⁺¹⁾= U⁽ⁿ⁾−∆t 2

1 2

L(U⁽¹⁾) + L(U⁽²⁾

(3.12) This scheme is simply two-steps Runge-Kutta method and is second order. The order of accuracy of the scheme is usually equal to the number of sub-steps with same weights for each sub-steps. The more general form of the two-stage Runge-Kutta method reads:

U⁽ⁿ⁺¹⁾= U⁽ⁿ⁾− ∆t 2

1 2

aL(U⁽¹⁾) + bL(U⁽²⁾

(3.13) Where a and b are the weighted factors. The summation of the weighted factors should be one. Similarly the method can be extended to higher orders.

Runge-Kutta schemes are usually easy to be implemented because they are explicit schemes and mostly shows stable behavior for a wide range of Courant numbers (for example, fourth order Runge-Kutta scheme is stable if C2√

2) (Hoffmann & Chiang, 2000).

The major problem with this approach is that since it requires more computation at each time step, it is more memory intensive. In order to solve this problem, modified Runge- Kutta scheme has been introduced to eliminate averaging. In the present study, the third

(20)

CHAPTER 3. NUMERICAL IMPLEMENTATION 13

order modified Runge-Kutta scheme is used.

Using third order Runge-Kutta scheme to integrate the constitutive equation (equation 2.10), the temporally discretized equation reads as equation 3.14.

C_ij^(l)− C_ij^(l−1)

∆t = γ_l

−u_k∂Cij

∂x_k + C_kj∂ui

∂x_k + C_ik∂uj

∂x_k.

(l−1)

+ ζ_l

−u_k∂Cij

∂x_k C_kj∂ui

∂x_k + C_ik∂uj

∂x_k.

(l−1)

−α_l 1 W e





Cij

1 −^{T r(C}_L2^ij⁾

!(l)

− Cij

1 −^{T r(C}_L2^ij⁾

!(l−1)

− 2δ_ij





(3.14) In equation 3.14, subscripts l denotes the third order Runga-kutta sub-steps, and γ1 = ₁₅⁸, γ₂ = ¹⁵₁₂, γ2 = ₄³ , ζ1 = 0, ζ2 == −¹⁷₆₀, ζ2 == −₁₂⁵, α1 = ₁₅⁴, α2 = ₁₅¹ , α3 = ¹₆ ,are the corresponding Runga-Kutta coefficients.

3.3 Spatial discretization

In this section, the weighted essentially non-oscillatory (WENO) schemes for spatial discretization are introduced. To study these scheme, first consider a simple scalar hyperbolic equation in conservative form (Xing & Shu, 2005).

ut+ f (u)x = 0 (3.15)

One can use conservative flux to discretize f (u)xas follow on a staggered grid:

f (u)_x = fˆ_i+¹

2

− ˆf_i−¹

2

∆x

(3.16)

Where ˆf_i+¹

2 and ˆf_i−¹

2 are numerical fluxes and are computed through the neighbouring point values. Considering fj = f (u_j), the detailed steps one needs to follow to calculate (2k − 1)^th order WENO scheme is presented in this section.

At the first step, k numerical fluxes (of the order k) should be calculated using k different stencils (Sr(i), where r = 0, 1, 2, , k − 1).

fˆ_i+^(r)1 2

=

k−1

X

j=0

crjfi−r+j (3.17)

Considering the fifth order WENO scheme, which is used in the present study, means that k=3. Therefore, one can write three different third order fluxes as below:

fˆ_i+⁽⁰⁾1 2

= 1 3fi+5

6fi+1−1 6fi+2

fˆ_i+⁽¹⁾1 2

= −1

6fi−1+5 6fi+ 1

3fi+1

fˆ_i+⁽²⁾1 2

= 1

3fi−2−7

6fi−1+11 6 fi

(3.18)

Finally, the (2k − 1)^th order flux of WENO is the combination of all the above mentioned fluxes.

(21)

fˆ_i+^(r)1 2

=

k−1

X

j=0

wrf_i+¹

2 (3.19)

Where wr represents the non-linear weights and can be calculated through the following equations:

w_r= α_r Pk−1

s=0αs

αr = dr ( + βr)²

(3.20)

Before introducing the parameters in the relations, it is important to notice that wr ≥ 0 and Pk−1

j=0w_r = 1. In 3.20, r is the linear weight, βr is the smoothness indicator of the stencil Sr(i), and is a small constant (typically 10⁻⁶). Simply, dr which makes the scheme to be of the (2k − 1)^th order of accuracy, βr considers the smoothness of the flux function, and guarantees the non-zero denominator.

For the fifth order WENO scheme the linear weights and the smoothness indicators can be written as (Xing & Shu, 2005):

d0= 3

10, d1 = 3

5, d2 = 1 10 β₀ = 13

12(f_i− 2f_i+1+ f_i+2)²+1

4(3f_i− 4f_i+1+ f_i+2)² β₁= 13

12(f_i−1− 2f_i+ f_i+1)²+1

4(f_i−1− f_i+1)² β₂ = 13

12(f_i−2− 2f_i−1+ f_i)²+1

4(f_i−2− 4f_i−1+ 3f_i)²

(3.21)

(22)

Chapter 4

Numerical Simulation

4.1 Problem description and solution steps

The goal of this chapter is to explain the studied problem ( geometry, generated grid, etc.) delineate the numercial simulation details (solution steps) and present a verification for the developed code.

Turbulent flow of an incompressible visco-elastic fluid through a square duct is considered.

Figure 4.1 shows a sketch of the geometry and the Cartesian coordinate system, where x, y and z (x1, x2, and x3) denote the streamwise, wall-normal and spanwise coordinates, while u, v and w (u1, u2, and u3) denote the respective components of the velocity vector field. The lower, upper walls are located at y = ±h, while the left and right walls at z = ±h, being h the cross section half height. The Reynolds number of the flow is defined as Re = Ubh/ν, where h is chosen as characteristic length scale, the characteristic velocity is the bulk velocity Ub, defined as the average value of the mean velocity computed across the whole domain, and νthe total kinematic viscosity, i.e., the sum of the solvent and polymer viscosity. The fluid is governed by the Navier-Stokes (NS) equations (equation 2.16)

To model the additional stresses due to the presence of polymers in the flow, the FENE-P model is used where the polymer stress tensor τij is written as a function of the configuration tensor Cij as equation 2.11. The configuration tensor is a symmetric second order tensor, which is found by solving the equation 2.10.

y

z

x

h

-h -h h

Figure 4.1: Sketch of the square duct geometry with the walls located at y = −h and y = h, and at z = −h and z = h.

15

(23)

16 CHAPTER 4. NUMERICAL SIMULATION

4.2 Numerical implementation

The code for the numerical solution of the equations above is an extension of that described by Rosti & Brandt (2017), with whom it shares its general architecture. The time integration method used to solve the equations is based on an explicit fractional-step method (Kim &

Moin, 1985), where all the terms are advanced with the third order Runge-Kutta scheme, ex- cept the viscous stress contribution which is advanced with Crank-Nicolson (Min et al., 2001).

In particular, to solve the system of governing equations, at every time step, the following steps (see also Dubief et al., 2005; Min et al., 2001) have been performed: i) the updated conformation tensor Cij is found by solving equation (2.10); ii) the polymer stress tensor τij is computed by equation (2.11); iii) the NS equations (equation (2.16)) 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).

The numerical solution of equation 2.10 is cumbersome, and many reaserchers showed that the numerical solution of a viscoelastic fluid is unstable, especially in the case of high Weissenberg numbers, since any disturbance amplifies over time (Dupret & Marchal, 1986;

Sureshkumar & Beris, 1995; Min et al., 2001). Indeed, 2.10 can easily diverge and lead to the numerical breakdown of the solution since it is an advection equation without any diffusion term (Dubief et al., 2005). One of the earliest solution to this problem has been to introduce an artificial diffusivity (AD) to the constitutive equations (Sureshkumar & Beris, 1995; Mom- pean & Deville, 1997; Alves et al., 2000). Subsequently, global AD was replaced by local AD, a method which has been widely used.

In the present work, global AD in the first transient part of the simulation with polymers was used which was gradually removed while approaching the final statistical steady state regime.

Also,implementation of a fifth order Weighted Essentially Non-Oscillatory (WENO) scheme (Jiang & Shu, 1996) for the advection terms in equation 2.10 has been chosen. Apart from that, the governing differential equations are solved on a staggered grid using a second order central finite-difference scheme. This methodolgy has been proved to work properly by Sugiyama et al. (2011) and also successfully used in Rosti & Brandt (2017) for viscoelastic solid-like materials. A comprehensive review on the properties of different numerical schemes for the advection term is reported by Min et al. (2001).

4.3 Numerical details

For all the turbulent flows considered hereafter, the equations of motion are discretised by using 448 × 448 × 448 grid points on a computational domain of 12h × 2h × 2h in the streamwise, wall-normal and spanwise directions. The spatial resolution has been chosen in order to properly resolve the wall turbulence, satisfying the constraint ∆y⁺ = ∆z⁺ ≈ 0.8 and

∆x⁺≈ 5.

Viscous units, used above to express spatial resolution, will be often employed in the following; they are indicated by the superscript⁺, and are built using the friction velocity uτ

as velocity scale and the viscous length δν = ν/uτ as length scale. For a turbulent Newtonian

(24)

CHAPTER 4. NUMERICAL SIMULATION 17

duct flow, the dimensionless friction velocity is defined as

uτ = s

1 Re

∂u

∂y, (4.1)

where the overline indicates the mean, and the derivative is taken at the wall (note that, the equation is written for the top and bottom wall, but a similar one can be written for the left and right by changing y with z). When the flow is visco-elastic, the definition (4.1) must be modified to account for the polymeric stress which is in general non-zero at the wall.

Thus,one can define:

u_τ = s

β Re

du

dy +1 − β

Re τ₁₂. (4.2)

All the simulations have been performed at constant flow rate consistently with choosing U_b as the characteristic velocity; the flow Reynolds number based on the bulk velocity is fixed at 2800, i.e., Re = ρUbh/µ = 2800, where the bulk velocity is the average value of the mean velocity computed across the whole domain occupied by the fluid phase, and µ is the total dynamic viscosity. This choice facilitates the comparison between the Newtonian and visco-elastic flow. Since the constant flow rate condition is being enforced, the appropriate instantaneous value of the streamwise pressure gradient is determined at every time step.

All the simulations have been started from a fully developed turbulent duct flow without polymer additives. After the flow has reached statistical steady state, the calculations are continued for an interval of 400h/Ub time units, during which 400 full flow fields are stored for further statistical analysis. To verify the convergence of the statistics, they have been computed using different number of samples and verified that the differences are negligible.

4.4 Code validation

The code has been extensively validated in the past for turbulent flow simulations (Rosti &

Brandt, 2017). Here, one more comparison with literature results is provided for the case of a turbulent duct flow. Indeed, the results have been validated by comparing the statistics of a Newtonian turbulent duct flow at Re = 2800 with those reported by Uhlmann et al. (2007) and Pinelli et al. (2010) at a slightly higher Reynolds number Re = 2900. Figure 4.2 shows the mean velocity profile in logarithmic scale in the top row and the cross term of the Reynolds stress tensor u⁰v⁰as a function of the distance from the wall, hereafter indicated with ˜y. Three different section are considered: x = 0 the centerline, x = 0.3h and x = 0.6h. Overall good agreement between the literature (shown with the black symbols) and the simulation results (shown with the solid line) is evident.

Further validation will be provided throughout the Manuscript by comparing the results of the simulations with experimental measurements, both in the Newtonian and visco-elastic cases.

(25)

18 CHAPTER 4. NUMERICAL SIMULATION

Figure 4.2: (top) Mean streamwise velocity component u as a function of the wall normal distance, in wall units. (bottom) Mean u⁰v⁰ component of the Reynolds stress tensor profile.

In the figures, the left, middle and right panels correspond to the section z = 0, 0.3h and 0.6h, respectively. Solid lines are used for simulations results, while the symbols are the results those from Uhlmann et al. (2007) and Pinelli et al. (2010).

(26)

Chapter 5

Experimental setup

Complementary experimental measurements have been performed by Sagar Zade ¹ at the laboratory facilities at KTH, Stockholm.

P

F

Tank

Vortex breaker

Static mixer

Flow Traversable Laser

Pressure transducer

Flow

meter Pump

PIV camera and computer

Figure 5.1: (left) Set-up and (right) photo of the experimental facility used in the present work.

The experiments were performed in a square duct of size 50 mm × 50 mm in cross section and 5 m long. The entire duct is made of transparent Plexiglas permitting visualization throughout its length. Figure 5.1 shows the schematic of the flow loop and a photo of the facility. The fluid is recirculated through a conical tank, open to the atmosphere, where the polymer solution can be introduced. A vortex breaker has been installed at the conical end of the tank to prevent entrainment of air due to suction at higher flow rates, and a static mixer is mounted close to the inlet of the duct with the intention of neutralizing any swirling motions that may arise from the gradual 90^o bend at the exit of the tank. A section that provides a transition from a circular to a square cross-section follows the mixer. Finally, a friction tape is lined on the inner walls of the Plexiglas duct to trigger turbulence. Note that, the temperature of the solution is maintained close to 20C by means of an external heat exchanger in the tank, and flexible corrugated piping is used to minimize the induction of vibrations from the pump to the duct. A very gentle disc pump (Discflo Corporations, CA, USA) has been chosen to minimize the degradation of the polymers’ long chains and the consequent loss of drag reduction over time. For the same reason, the time of the experiments has been kept

1zade@mech.kth.se

19

(27)

20 CHAPTER 5. EXPERIMENTAL SETUP

Figure 5.2: Friction factor f as a function of the Reynolds number (based on the total height 2h) for a Newtonian duct flow. Experimental measurements, shown with the blue solid symbols, are compared with the empirical correlation (dahsed line) provided by Duan et al.

(2012) and reported in equation (5.1), and with the numerical results by Uhlmann et al. (2007) and Pinelli et al. (2010) shown with the black symbol.

short enough to maintain a reasonably constant value of drag reduction during the whole data acquisition, but long enough to get statistically stationary results.

5.1 Measurement technique

The velocity field in the square duct was measured using 2D Particle Image Velocimetry (2D-PIV) in 3 x − y planes: z = 0 (center plane), 0.3h and 0.6h, where h is the half-width of the duct. A traversable continuous wave laser (532 nm, 2 W) and a high-speed camera (Phantom Miro 120, Vision Research, NJ, USA) were used to capture successive image pairs (see the right panel in figure 5.1). For imaging the full height of the duct, a resolution of approximately 60 mm/1024 pixels was chosen, with a frame rate such that the maximum pixel displacement did not exceed 1/4^ththe size of the smallest interrogation window (Raffel et al., 2013). Images were processed using an in-house three-steps, FFT-based cross-correlation al- gorithm (Kawata & Obi, 2014): the first step consists of basic PIV with a large interrogation window size (48 × 48 pixel), followed by the discrete window shift PIV at the same resolution, and, finally, the central difference image correction method (Wereley & Gui, 2003) with the smallest interrogation window size (32 × 32 pixel). The degree of overlap can be esti- mated from the fact that the final resolution is 1 mm × 1 mm, which is suitable for velocity statistics in the bulk of the flow away from the walls but it is inadequate in the near-wall region. Therefore, an independent measurement is conducted by focusing the camera on a region close to the wall in order to measure the near-wall velocity variation at a resolution of 0.35 mm × 0.35 mm. Usually, from 500 to 1000 image pairs were sufficient to ensure

(28)

CHAPTER 5. EXPERIMENTAL SETUP 21

statistically converged results, i.e., the measurement time was between 17 to 35 minutes.

An electromagnetic flowmeter (Krohne Optiflux 1000 with IFC 300 signal converter, Krohne Messtechnik GmbH, Germany) was used to measure the flow rate, while the pressure drop is measured across a length of 54h starting from 140h from the inlet (fully developed turbulent flow) using a differential pressure transducer (0 − 1 kPa, Model: FKC11, Fuji Electric France, S.A.S.). At the two chosen streamwise locations, two pressure tappings (one at the top wall and other at the bottom wall) were precisely drilled in the wall-normal direction at the center plane of the duct, and the pressure tubes emerging from these holes were joined together into a single tube, connected to the high- (140h) and low-pressure (194h) sides of the pressure transducer. Additional pressure tappings were also made to study the evolution of pressure drop over the streamwise length. Finally, data acquisition from the camera, flow meter and pressure transducer was performed using a NI-6215 DAQ card using Labview^TM software.

Figure 5.2 shows a good agreement between the friction factor for a fully developed flow measured in the square duct and that predicted by the empirical correlation given by Duan et al. (2012) for a Newtonian fluid (water in our case):

f =

3.6 log 6.115 Re_2h

−2

. (5.1)

Note that, here the Reynolds number is based on the characteristic length given by the square root of the cross section area, i.e., 2h in our case. Also, the friction velocity is obtained from the pressure drop using the following relation:

u_τ = s

h 2ρ

dp

dL (5.2)

where ρ is the density of the solution.

5.2 Polymers

All the experiments have been performed using an aqueous solution containing 250 ppm of a high molecular weight polyacrylamide based anionic polymer (FLOPAM^TMAN934SH, SNF). The density of the solution is practically the same as the solvent (filtered tap water) under such dilute concentrations. The solution was prepared by dissolving the polymer in powder form into water, followed by successive dilution and mixing to ensure its homogene- ity. Such a solution with the desired polymer concentration was stored in a large reservoir from which it was pumped to the tank at the beginning of every new experiment. The shear viscosity of the polymer solution was measured, both at the beginning and at end of the experiment using a coaxial cylinder viscometer (Brookfield DV-II + Pro, Brookfield Engineering Laboratories, Inc. MA, USA). To define the Reynolds number with the visco-elastic polymer flow, a value of viscosity was used corresponding to the average experimental shear stress obtained from the pressure drop measurements in the fully developed region of the duct.

(29)

Chapter 6

Results

The analysis starts by considering two cases at Re = 2800, Newtonian fluid and the viscoelastic fluid consisting of a suspension of polymers as introduced in the previous section. In the experiments, the precise characterization of the polymer solution is not available, thus, in order to provide a reasonable comparison between experimental and numerical results, First a series of simulations varying the viscoelastic properties of the flow have been run, and finally the one providing the same friction Reynolds number and drag reduction as in the experiment was chosen, resulting in β = 0.9, L² = 3600 and W i = 1.5 (W iτ0 = λu²_τ₀/ν = 18being uτ0 the friction velocity of the Newtonian simulation). The New- tonian case has an average friction Reynolds number Reτ equal to 183, while the polymer one provides a value of 155 which corresponds to a drag reduction DR of approximately 29%, computed as DR = −∆τw/τ_w0, where τw is the total shear stress at the wall (as in the computation of the friction velocity), averaged over the entire four sides, and τw0 the value for the Newtonian case. The Newtonian results will be represented hereafter in blue and the polymer ones in red; numerical results will be reported with a solid line and the corresponding experiments with symbols.

6.1 Flow statistics

Figure 6.1 shows the contour of the mean streamwise velocity component u (panel a) and of the mean streamwise vorticity component ωx (panel b) determined numerically, where in each panel, the left half of the duct reports the data for the Newtonian flow, whereas the right side the viscoelastic counterpart. The Newtonian case shows the streamwise velocity contour typical of a duct flow, characterized by the symmetry inherited from the considered geometry. The maximum streamwise velocity is located at the center of the duct and equals 1.4U_b. The streamwise vorticity displays the in-plane secondary flow typical of turbulent duct flows: 8 regions (4 in the half section shown in the figure) can be easily identified with alternate sign of vorticity, located symmetrically with respect to the horizontal and vertical lines passing through the center and with the two diagonals. In particular, the location of one of the maxima of the vorticity is marked with a filled black circle, which was found at x = −0.74h, y = −0.9h in the Newtonian case. In the polymeric flows, the streamwise velocity and vorticity components are altered: the streamwise velocity contour exhibits higher values close to the corners and, at the same time, the secondary flow is modified, with the locations of the maximum vorticity moving towards the center, i.e., they are displaced away from the walls. In fact, the coordinates of the selected maximum vorticity become x = 0.23h,

22

(30)

CHAPTER 6. RESULTS 23

Figure 6.1: Contour of the mean streamwise component of (left) velocity u and (right) vorticity ωx. In each figure, the left and right half of the duct are used for the Newtonian and viscoelastic fluids, respectively. The color scale ranges from 0 (blue) to 1.4Ub(red) in the left panel and from −0.5Ub/h(blue) to 0.5Ub/h(red) in the right one.

Figure 6.2: Mean streamwise velocity as function of the wall-normal distance in logarithmic scale and wall units. The blue and red lines are used for the numerical results without and with polymers, while the symbols of the same color for the experimental measurements. The magenta dash dotted line is the polymer Maximum Drag Reduction (Virk, 1975; Lvov et al., 2004) while the black one the Newtonian log law. The three panels correspond to different sections: (left) z = 0, (middle) z = 0.3h, and (right) z = 0.6h.

y = 0.78h. The magnitude of the vorticity maxima are also modified, reducing by 7% from 0.29Ub/hin the Newtonian fluid to 0.27Ub/hin the presence of polymers, while the overall integral of the vorticity in each of the 8 sectors increases by 17% from 0.023Ub/hin the Newtonian flow to 0.027Ub/hin the viscoelastic one. Thus, although the peak value of the secondary flow is reduced, overall the secondary flow is enhanced by the polymer addition due to an enlargement of the area with non-zero vorticity. The enhancement of the secondary flow is further proved by analysing the maximum of the magnitude of v in the domain; indeed, this grows from approximately 0.006Ub in the Newtonian flow to 0.012Ub with the polymers, both in the numerical and experimental results. This result differs from which was reported by Escudier & Smith (2001) who on the contrary found a reduction of v induced by the addition of polymers. The opposite behavior may be due to the different flow conditions; indeed, in their experiment the Reynolds number and the polymer concentration are much higher than our results.

Figure 6.2 shows the mean streamwise velocity component u as a function of the wall normal distance ye(both in wall units). The left panel shows the mean velocity profile in

(31)

24 CHAPTER 6. RESULTS

Figure 6.3: Wall-normal mean profile of the (top row) streamwise velocity fluctuation u⁰, (middle row) wall-normal velocity fluctuation v⁰, and (bottom row) cross term of the Reynolds stress tensor u⁰v⁰. The three columns correspond to different sections: (left) z = 0, (middle) z = 0.3h, and (right) z = 0.6h. The line and color style is the same as in figure 6.2.

the middle plane (z = 0), the middle panel the velocity at z = 0.3h, and the right one at z = 0.6h(i.e., the same sections where experimental data are available). In all the panels, the numerical (solid lines) and experimental (symbols) results are in good agreement for both the Newtonian and viscoelastic fluid. In the middle plane, where the secondary flow is weak, one can identify three regions in the velocity profile of the Newtonian fluid (left panel) similarly to what found for a turbulent channel flow: first, the viscous sublayer for ey⁺ < 5where the variation of u⁺ withye⁺ is approximately linear; then, the so-called log- law region,ey⁺ > 30, where the variation of u⁺ versusye⁺is logarithmic; finally, the region between 5 and 30 wall units is called buffer layer and neither laws hold. By approaching the side walls, x = 0.3h and 0.6h, the range of the equilibrium layer where the velocity profile is logarithmic reduces (middle panel), eventually disappearing (right panel). The profiles have similar trends for the viscoelastic fluid, in which case, however, the inertial range has an upward shift, consistent with the observed drag reduction (Warholic et al., 1999). Note that, all the data stays below the MDR curve (Virk, 1975; Lvov et al., 2004), a limit curve approached only in the case of very high W i, shown as a dashed line in the figure.

Figure 6.3 displays the velocity fluctuations profile as a function of the wall-normal dis- tancey. In particular, the top row reports the streamwise component of the velocity fluctu-e ations u⁰, the middle panel one of the wall normal components v⁰, and the bottom row the

Turbulent Duct Flow with Polymers

Turbulent Duct Flow with Polymers

ARMIN SHAHMARDI

Turbulent Duct Flow with Polymers

Contents

Chapter 1

Introduction

1.1 Aim and objectives

1.2 Duct flow

1.3 Polymer additives

1.4 Polymer solutions in square duct flows

1.5 Outline

Chapter 2

Mathematical Formulation

2.1 FENE model

2.2 FENE-P model

2.3 Navier-Stokes equation

Chapter 3

Numerical Implementation

3.1 Numerical schemes for incompressible Navier-Stokes

3.2 Time advancement

3.3 Spatial discretization

Chapter 4

Numerical Simulation

4.1 Problem description and solution steps

4.2 Numerical implementation

4.3 Numerical details

4.4 Code validation

Chapter 5

Experimental setup

5.1 Measurement technique

5.2 Polymers

Chapter 6

Results

6.1 Flow statistics