a numerical study
8518
599G
Bibliotheek TU Delft
a numerical study
Proefschrift
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,
op gezag van de Rector Magnificus Prof. dr. ir. J.T. Fokkenia voorzitter van het College voor Promoties,
in het openbaar te verdedigen op dinsdag 13 m(>i 2008 om 15.00 uur
door
Chiara TESAURü
Ingegnere Chimico, Universita degli studi di Salerno geboren te Milano, Italïe.
TU Delft Library
Prometheuspjöin 1
2628 ZC Delft
Prof. dr. ir. B. ,1. Boersma Prof. dr. ir. G. Ooms
Samenstelling promotieconunissie; Rector Magnificus,
Prof. dr. ir. B. .J. Boersma. Prof. dr. ir. G. Ooms, Prof. dr. ir. C. Vuik, Prof. dr. ir. H.I Andersson, Prof. dr. R. Piva.
Prof. dr. ir. B.J. Geurts. Dr. ir. M.A. Hulsen,
voorzitter
Technische Universitet Delft, promotor Technische Universitet Delft, promotor Technische Universitet Delft
NTNU, Noorwegen
Universita degli studi di Roma, Italïe Universiteit Twente
Technische Universitet Eindliovtm
Dr. ir. M.A. Hulsen heeft als begeleider in belangrijke mate aan de totstandkoming van het proefschrift bijgedragen.
The work ]Drescnted in this thesis was supported financially by the Dutch Fouudation for Fundamental Research on Matter (FOM)
Copyright ©2008 by C. Tcsauro All rights reserved.
Contents
S u m m a r y x i S a m e n v a t t i n g xiii 1 I n t r o d u c t i o n 1
1.1 Motivation and background 1
1.2 Polymers 3 1.3 Turbulence 6 1.4 Numerical simulation of Turbulence 7
1.5 Rotation 8 1.6 Polymer modeling 8
1.7 Outline of this thesis 9 2 M a t h e m a t i c a l b a c k g r o u n d a n d n u m e r i c a l m e t h o d 1 1
2.1 Introduction 11 2.2 Governing equations for the flow field 11
2.3 Constitutive equations 12 2.3.1 Elastic dumbbell model 12
2.3.2 Polymer stress tensor 14 2.3.3 FENE-P model 15 2.4 Dimeiisionless equations 16 2.5 Numerical approach 17 2,5.1 Simulations techniques 17 2.6 Parallel inii)lementation 22 3 D N S of d r a g r e d u c t i o n in a r o t a t i n g t u r b u l e n t flow 25 3.1 Introduction 25 3.2 Mathematical model 27
3.2.1 Governing equations for the flow field 27
3.2.2 Constitutive equation 28 3.2.3 Dimensionless equations for the flow field 29
3.2.4 Dimensionless evolution equation 30
3.3 Computational procedures 30 vii
3.4 Results and discussion 32 3.4.1 Drag reduction 32 3.4.2 DNS results for the velocity 34
3.4.3 Stress balance 40 3.4.4 Polymer stretching and orientation 44
3.4.5 Two-points correlations 45 3.4.6 Turbulence spectrum 45 3.4.7 Energy budgets 49 3.5 Conclusions 51 4 E v e n t s of h i g h p o l y m e r a c t i v i t y in d r a g r e d u c i n g flows 55 4.1 hitroduction 55 4.2 Governing equations and computational procedures 56
4.3 Elastic energy of the polymers 59
4.4 Vortex detection 60 4.5 Results and discussion 61
4.5.1 Flow statistics 61 4.5.2 Energy budget 64 4.5.3 Vortex detection 69 4.6 Polymer conformation tensor 71
4.7 Conclusions 75 5 D N S u s i n g t h e l o g - c o n f o r m a t i o n m e t h o d 77
5.1 Introduction 77 5.2 Governing equations 79 5.3 The log-conformation representation 80
5.4 Numerical method 81 5.5 Analytical solution of FENE-P model in laminar flow 82
5.6 Results 85 5.6.1 Laminar flow 85 5.6.2 Turbulent flow 86 5.7 Conclusions 87 6 C o n c l u s i o n s 91 B i b l i o g r a p h y 93 N o m e n c l a t u r e 97 A c k n o w l e d g m e n t s 101 C u r r i c u l u m V i t a e 103
Summary
S o m e a s p e c t s of d r a g r e d u c i n g flows, a n u m e r i c a l s t u d y
C h i a r a T e s a u r o
-It is a well known fact that the addition of a small amoimt of polymer to a turbulent Newtonian fluid can lead to significant decrease in the frictional drag. This phenomenon was for the first time reported by Toms in 1949.
Drag reduction is important from a practical point of view. Among the many potential applications we can mention the possibility to increase the throughput of pipelines, irri-gation systems and sewer systems. Also from a theoretical point of view drag reduction is important. For flows with high Reynolds numbers the polymers will be much smaller than the smallest scale of turbulence, so common sense tells us t h a t they can only act on the small scales (or micro scale of turbulence). However, experimental observations show that they are able to substantially change the large scale (frictional drag). A proper understanding of this process could give guidelines for turbulence control.
Since its discovery more than fifty years ago, drag reduction has been studied exten-sively by experiments and by numerical simulations. Today, there are some explanations for the phenomena. However, these explanations are non conclusive and the mechanism of drag reduction remains not fully understood.
The main purpose of the present research is to better understand the mechanism of drag reduction and its implications for practical applications. This is done with help of a Direct Numerical Simulation of the governing equations. Apart from mean flow cjuantities, such a simulation provides quantities which are difficult to obtain experimentally. For instance the polymer stress and the rate of energy transfer between polymers and turbulence.
The polymers are modeled with a FENE-P model. The simulations have been carried out with a two way coupled scheme: polymers are deformed by the flow leading to an extra stress which on its turn influences the flow.
Direct numerical simulations have been carried out for different values of the polymer extensibility and relaxation time to study their effects on drag reduction. For the flow geometry we use a rectangular channel. It is shown that the amount of drag reduction mcreases with increasing extensibility and relaxation time of the polymers. In particular, with increasing polymer extensibility we found resiflts that are close to the maximum drag reduction or Virk asymptote.
In most flow loops fluid elements are subjected to rotation and curvatme. It is therefore important to study the effect of rotation and curvature on drag reduction. In this study we have restricted ourselves to the effect of rotation. It is shown t h a t rotation in the perpendicular cross section of the channel has a large effect on drag reduction: a small rotation is enough to significantly decrease the drag reduction.
In the high drag reduction regime, polymers and turbulence strongly interact with each other. This interaction takes place in well defined areas near the wall. With an analysis of the local topology of drag reducing flows we found t h a t polymers are able to decrease the amount of vorticity and increase the rate of strain in the flow.
Furthermore, it has been found t h a t the polymers, which are strongly stretched espe-cially near the wall, tend to orient themselves in the direction of the wall. This effect is more evident in the high drag reduction regime. This result can explain the lower amount of drag reduction which has been found in rotating flows, where the secondary motion inhibits the possibility of the polymers to orient themselves in a preferred direction.
In the last part of the thesis, we developed a new method for the simulation of drag reduction. The mathematical character of the polymer equation is different from the flow equations. It is possible t h a t during the time-integration of the constitutive equation numerical errors lead to a loss of the positive definiteness of the polymer conformation tensor. Eventually this can lead to a break-down of the calculations. The positive defl-niteness of this tensor is related to the length of the polymers, which physically can not become negative. The new method, which builds on the method which was proposed by Fattal and Kupferman, and later implemented by Hulsen for an Oldroyd-B and Giesekus model, uses a change of the variables. This transformation removes the instabilities which cause the break-down of the calculations. We applied this method to the DNS of polymer flow, both for the laminar and turbulent case. For the laminar case the results were in good agreements with the results obtained with the analytical solution for channel flow. For the turbulent case the calculations broke down before they reach a statistical steady state.
Samenvatting
A s p e c t e n v a n w e e r s t a n d s v e r m i n d e r e n d e s t r o m i n g e n , e e n numeriek o n d e r z o e k Chiara Tesauro
-Het toevoegen van een kleine hoeveelheid polymeren aan een turbulente Newtoniaanse stroming kan aanleiding geven tot een forse weerstandsvermindering. Dit verschijnsel is voor het eerst door T(nns beschreven in 1949.
Vanuit een praktisch oogpunt is weerstandsvermindering belangrijk. Er zijn vele toe-passingen te noemen, zoals de mogelijkheid tot verbeteren van de efficiëntie van pijpleidingen, het verbeteren van irrigatie systemen en afvalwater systemen. Ook vanuit een theoretisch oogpunt is weerstandvermindering belangrijk. Voor stromingen met een hoog Reynolds getal zullen de polymeren altijd veel kleiner zijn dan de kleinste schaal van de stroming, de zogenaamde microstructuur. Het lijkt logisch om te veronderstellen dat deze poly-meren alleen werkzaam kunnen zijn op de kleine schalen. Echter, metingen tonen aan dat polymeren ook in staat zijn om de grootschalige structuren te veranderen. Een goed begrip van dit proces zou richtlijnen kunnen geven voor de controle van tmbnlentie.
Sinds de ontdekking van het verschijnsel, nu meer dan 50 jaar geleden, is er veelvuldig experimenteel en numeriek onderzoek gedaan naar het effect van polymeren op een tur-bulente stroming. Er zijn in de literatuur wel verklaringen gegeven voor het verschijnsel. Er is echter nog geen sluitende theorie beschikbaar.
Het hoofddoel van het huidige onderzoek is het krijgen van een beter begrip van het mechanisme dat voor weerstandsvermindering verantwoordelijk is. We maken hiervoor gebruik van een Directe Numerieke Simulatie van de bewegingsvergelijkingen. Naast gemiddelde grootheden geeft een dergelijke simulatie ook grootheden die erg moeilijk of onmogelijk te meten zijn in een experiment. We denken hierbij bijvoorbeeld aan de poly-meer spanningstensor en de energieoverdracht van polymeren naar de turbulentie. De polymeren worden gemodelleerd met het FENE-P model. De simulaties zijn gekoppeld uitgevoerd. De polymeren worden gedeformeerd door de stroming en de extra poly-meerspanning beïnvloedt de stroming weer.
Directe numerieke simulaties zijn uitgevoerd voor verschillende waardes van de poly-meer extensibiliteit en relaxatietijd. Het effect van deze parameters op de weerstandsver-mindering is bestudeerd. Voor de geometrie maken we gebruik van een rechthoekig kanaal. De
vergelijkingen worden opgelost met behulp van een pseudo-spectrale methode. Het wordt aangetoond dat de weerstandsvermindering groter wordt met toenemende extensibiliteit en relaxatietijd. Wanneer de extensibiliteit van de polymeren wordt vergroot gaat de hoeveelheid weerstandsvermindering dicht naar de maximum weerstandsreductie limiet, de zogenaamde Virk asymptoot.
In de meeste leidingsystemen worden vloeistofelenienten ook beïnvloed docjr rotatie en kromming, fiet is daarom ook belangrijk om te bestuderen wat het effect van ro-tatie en kronnning is op weerstandsvermindering. In deze studie zullen we ons beperken tot het effect van rotatie. We zullen aantonen dat rotatie een grote invloed heeft op weerstandsvermindering, een lage rotatie is al voldoende om een behoorlijke afname in weerstandsvermindering te geven.
In het regime met hoge weerstandsvermindering is er een sterke koppeling tussen poly-meren en de turbulentie. Deze interactie vindt plaats in een goed gedefineerd gebied relatief dicht bij de wand van het kanaal. Een analyse van de locale tojjologie van de stroming toont aan dat in een stroming met weerstandsvermindering polymeren in staat zijn de vorticiteit te reduceren en de locale strekking van de vloeistofelenienten toe te laten nemen. Verder vinden we dat de polymeren sterk gestrekt zijn, met name in het gebied dicht bij de wand. De oriëntatie van de polymeren is ook parallel aan de wand. Dit is vooral duidelijk in de stromingen met een hoge weerstandsvermindering. Dit resul-t a a resul-t zou kunnen verklaren waarom de weersresul-tandsvermindering in een roresul-terende sresul-troming veel kleiner is. In het geval van rotatie zijn sterke secondaire stromingen aanwezig die voorkomen dat het polymeer zich kan oriënteren in de richting van de wand.
In het laatste deel van dit proefschrift wordt een nieuwe methode besproken die we ontwikkeld hebben voor de sinmlatie van weerstandsvermindering. De wiskundige eigen-schapi^en van de polymeervergelijking zijn anders dan die van stromingsvergelijkingen. Het kan gebeuren dat bij het integreren van de polymeervergelijking kleine numerieke fouten gemaakt worden. Deze fouten kunnen ervoor zorgen dat het systeem niet meer positief-definiet is. Uiteindelijk zal dit leiden tot een instabiele berekening. Als het sys-teem niet meer positief-definiet is kunnen de polymeren negatieve lengtes krijgen, iets wat natuurlijk uit een fysisch oogpunt onmogelijk is. In de nieuwe methode, die gebaseerd is op het werk van Fattal en Kupferman en later gebruikt door Hulsen, wordt gebruik gemaakt van een coördinaat transformatie. Deze transformatie maakt de oplossing min-der gevoeling voor instabiliteiten. De methode is toegepast oj) laminaire en turbulente stromingen. Voor laminaire stromingen zijn de resultaten in overeenstemming met de an-alytische oplossing. Helaas zijn de instabiliteiten in de tiirbulente stroming nog dermate groot dat we geen statistisch stationaire toestand kmmen bereiken.
Chapter 1
Introduction
1.1 Motivation and background
More than 50 years ago Toms [49] disc'overed that the addition of a tiny amount of poly-methyl-methacrylate (lOwppiii) to a turbulent flow of monochlorobenzene, resulted in a large reduction in the pressure drop. Since the work of Toms many other types of polymers have been found to drastically reduce the pressure drop in a turbulent pipe or channel flow. The amount of drag reduction can be quite substantial as is shown in figure 1.1. In this figure the friction factor for water with and without a small amount of poly-ethylene oxide is shown. For the laminar flow, the friction factor for the polymer solution is identical to t h a t for pure water. In the turbulent regime (Re > 2,200), a decrease of the friction factor can be observed for the polymer solution.
The potential applications of drag reduction by polymer additives are numerous. An important application is the possibility to alter turbulent drag in pipelines and thus to increase the throughput. Burger et al. [8] investigated the use of drag reducing polymers for the Trans Alaska Pipeline. In this case the use of polymer led to a large increase of the flow rate, with consequent increase of the production per time unit. Other applications are: irrigation systems (Khalil et al. [25]), sewer systems (Sellin [44]), fire fighting systems (Rubin [42]), reduction of noise (Barker [1]), increase of speed of boats and ships by covering the skulls with a polymer layer (Tothill [50]) and circulation aids for diseases such as arteriosclerosis and coronary thrombosis (Greene et al. [17]).
The study of turbulent drag reduction builds on knowledge from two sub-fields of fluid mechanics, namely: turbulence and rheology. Polymers are, in general, much smaller than the smallest turbulent length scales in the flow. For this reason polymers are mostly active on the small scales. Due to the small polymer concentration, the physical properties of the solution do not differ appreciably from those of the pure solvent. However, from the applications we know that the jiolymers are able to modify also the large scales of the flow and to substantially change the character of turbulence with respect to a Newtonian flow. An understanding of the mechaiusm which alters the turbulence can help us to maxinnze the effort of drag reduction and may give guidelines for the active control of turbulence.
yQ-A 1 ^ I r I I I I I 1 i l I I I I I I 1 M I I I I I I I I
1 0 ' 1 0 ' 10* 10» 10» Reynolds number Re = D<v> p/ix
Figure 1.1: Friction factor for dilute aqueous solutions of poly-ethylene oxide as function of the Reynolds number (from Bird et al. [4]). The transition from a laminar to a turbulent flow occurs at a Reynolds number of ~ 2.2 • 10^.
Three theories have been proposed to explain the mechanism responsible for drag reduc-tion. The first one has been proposed by Lumley [30], who suggested t h a t the stretching of the polymer molecules by the flow increases the effective (extensional) viscosity and the thickness of the buffer layer, leading to drag reduction. In the second theory proposed by De Gennes [11] it is assumed that the drag reduction is caused by the elasticity of polymer molecules. Small scale turbulent kinetic energy is adsorbed and subsecjuently radiated away in form of shear waves. However, both theories are still (jualitative and do not give an explanation of the dynamics of drag reduction. Recently, Procaccia et al. [39] proposed a mechanism in which the polymers endow the fluid with a viscosity which is a linear function of the distance from the wall. A numerical simulation carried out by De Angelis et al. [10] showed a relation between a linear viscosity and the computed drag reduction.
In 1975 Virk [51] experimentally observed that for every value of the Reynolds number there is a specific upper limit on the anioimt of drag reduction in a pipe flow. This empirical result is called Maximmn Drag Reduction asymptote (MDR) or Virk asymptote (not shown in figure 1.1).
structural unit
-Figure 1.2: Schematic representation of linear and branched macromolecule of polyethylene, (a) linear, (b) branched (from Bird et al. [4]).
1.2 Polymers
Polymers are long molecules, which include plastics and rubbers. The use of polymeric materials is increasing rapidly over the years. In many applications they are replacing conventional materials such as metals, wood and natural fibers. A polymeric material is a substance composed of molecules which have long sequences of one or more species of atoms, or groups of atoms, linked to each other mostly by covalent bonds. The repeating groups of atoms are called structural units or monomers. A typical polymer is composed of 10 structural units, with a molecular weight of about 10^-10^ g/mol. A schematic representation of a polymer chain is shown in figure 1.2.
The aforementioned long chain nature sets polymers apart from other materials and gives rise to their particular properties. Sometimes polymer chains are branched, some-times connected with each other forming a network. A polymer molecule is an extremely complex molecular system with many degrees of freedom.
The flow of a fluid composed of macromolecules shows several differences compared to flow of a fluid composed of small molecules. We can, for instance, mention the "Weis-senberg effect" or "Rod climbing effect" (see figure 1.3). If a rigid rod is dipped into a polymer melt or solution, the polymer fluid climbs up the rotating rod. Another effect of the normal stresses in a polymer melt or solution is the "Extrudate swell" (see figure 1.3). If a polymer melt is forced to flow from a large reservoir through a circular tube, the diameter of the extrudate will be larger than the tube diameter, while for Newtonian fluids such as water the diameter of the extrudate will be smaller than the tube diameter. These effects are due to the viscoelastic behavior of polymers.
For a solid material, there is a linear proportionality between force and deformation of the material. This relation is known as Hooke's law;
r = G 7 (1.1) where T is the shear force per unit area and 7 is the shear deformation. The constant of
proportionality G is the so-called elastic modulus.
For a liquid substance, the proportionality between the shear stress r and the rate of strain 7 = ^ is given by Newton's law:
T = /X7. (1.2)
The constant fi is the so-called dynamic viscosity.
A polymer will not obey Hooke's or Newton's law. In general it will have a behavior
somewhere in between an elastic solids and a viscous liquid: the so-called viscoelastic behavior ^
The viscoelastic behavior of polymers is due to the long chain form of the macro-molecules. When a polymer is in equilibrium state, the chains are randomly oriented and the polymer is coiled. When a deformation is applied, the polymer uncoils and the molecular chain tends to orient in a preferred direction (one can just think what is hap-pening when one picks up some spaghetti from a plate using a fork). In a shear flow, the molecular chains will align more and more with the flow directions when the shear is in-creased. This is schematically illustrated in figure 1.4. When the deformation is stopped, the molecular chains tends to recoil.
The steady state shear behavior of a viscoelastic fluid can be described with the fol-lowing relation:
T = 77(7)7 (1.3) The function r; (7) is the so-called shear viscosity and it depends on the shear rate. A
fluid for which 77 is increasing with the shear rate is called shear thickening. A fluid for which rj is decreasing with the shear rate is called shear thinning. If the shear viscosity is constant, rj = fi, the fluid obeys Newton's law and it is called a Newtonian fluid.
The viscoelastic behavior is characterized by a time scale necessary to return from an extended state to an equilibrium state. The ratio between the time scale of the polymers
Figure 1.3: Rod climbing (top) and extrudate swell (bottom) for polymers. N indicates Newtonian fluid while P indicates polymer (from Bird et al. [4]).
Figure 1.4; Schematic representation of a long polymer chain in shear flow.
and the time scale of the deformation, can be indicated by a dimensionlcss number: the Weissenberg immber. The Weissenberg number is zero for a pure Newtonian fluid, and infinite for a pure elastic solid. Typical values of the Weissenberg number in a polymer solutions are of 0(10^).
For more information about polymers and their viscoelastic behavior we refer to: Barnes et al. [2], Bird ei al. [4] and Makosko [31].
1.3 Turbulence
In fluid dynamics, tiu'bulence is a flow regime characterized by variations of pressure and velocity in space and time. Turbulent flows are quite common in nature and in engineering applications: the flows in the oceans, rivers, canals and in the atmosphere are turbulent. In engineering applic:ation we can mention the flow of oil or gas in pipelines and the flow in large chemical reactors. A general review of turbulence can be found in standard textbooks such as Tennekes and Lumley [48].
In a turbulent flow there is a wide range of temporal and spatial scales. At one side of the temporal or spatial spectrum there is the so-called macro structure characterized by the large scales i which are in general of the order of the flow domain and have a characteristic velocity scale U. The kinetic energy of these large scale motions is received from the mean flow. On the other side of the spectrum there are the small scales, or micro-structures. The small scales depend on the material properties of the fluid. The smallest scales in the flow are known as the Kolmogorov scales. These scales are defined as:
where TJK and TK are the Kolmogorov length and time scale, v is the kinematic viscosity of the fluid and e is the dissipation rate, which is the rate of energy transported from the large scales to the small scales.
The turbulent energy injected at the large scales is transported through an energy cascade to the small scales where it is dissipated by viscous effect. The energy produced at the large scales is proportional to U^ and it is transferred to the small eddies in an eddy turnover time which is projjortional to l/l/£. The magnitude of the dissipation rate can then be estimated as:
' ^ " ^ (1.5)
From equations (1.4) and (1.5) it follows t h a t the ratio between the Kolmogorov scale and the macroscopic length scale is given by:
f ~ /ïe-3/4 (1.6)
were Re is the Reynolds number which is given by:
Re = ^ . (1.7)
V
The Reynolds number can be interpreted as the ratio between the inertia forces and the viscous forces acting on a fluid element and characterizes whether a flow is laminar or turbulent. For small values of the Reynolds number the viscous forces dominate and the flow will be laminar. For large values of the Reynolds number the inertia forces dominate and the flow will be turbulent.
1.4 Numerical simulation of Turbulence
With increasing computer power it has become possible to simulate a turbulent flow in full detail. Because of the wide range of spatial and temporal scales which must be resolved, by such a simulation, the numerical simulation of a turbulent flow is a demanding task.
Among the simulation techniques of turbulence we can mention Direct Numerical Si-niulation. DNS provides a numerical solution of the equations of conservation of mass and momentum in the three-dimensional reference system without any additional turbulence modeling. In the past it has been shown that DNS can be used to predict drag reduction in turbulent channel flow (Ptasinski et al. [41]).
The advantage of DNS is that it provides a full three-dimensional field of all flow variables. In a physical experiment this is not possible. Furthermore, it is possible to study the effects of various polymers properties like elasticity and concentration separately.
From equation (1.6) we can observe t h a t the separation of scales increases as the Reynolds number increases. For a three-dimensional problem the number of grid points should be proportional to the third power of the ratio of the macro- and micro-scale, which means t h a t the number of grid points which should be used for a DNS is proportional to Re'l'^. At present, DNS can therefore only be used for turbulent flows with low to moderate Reynolds immbers and simple flow geometries.
1.5 R o t a t i o n
Rotation is a common phenomenon in meteorology, geophysics and astrophysics. Also in engineering, many flows are subjected to rotation, for example the flow in rotating turbo-machinery such as gas turbines or centrifugal pumps. Flows with curvature show similar fiow patterns as flows in rotating systems. This is due to the similarity in body forces acting on flows in rotating (Coriolis force) and curved flow systems (centrifugal force). Due to the numerous applications in turbo-machinery, turbulent flows in rotating channels with the rotation vector orientated along the spanwize direction has been the subject of extensive studies.
Experiments on rotating flows were carried out by Johnston et al. [24], who performed detailed measurements of the flow patterns and demonstrated the effect of the Coriolis force on the flow.
Detailed analysis of rotating flow has been carried out with help of direct numerical simulation by, Kristoffersen & Andersson [29]. They have performed direct numerical simulation of a turbulent rotating channel flow. Their results show that the velocity statistics become strongly asymmetric with increasing rotation rate. Furthermore, they show that weak spanwise rotation only slightly affects the flow pattern. These effects have later f:>een confirmed by other authors (see for example Nagano & Hattori [37]).
1.6 Polymer modeling
In the past years, numerical simulations of drag reduction have been carried out with different polymer models. The numerical simulation of viscoelastic flows requires the additional solution of the constitutive ecjuations in order to describe the influence of the polymer on the flow. As mentioned in section 1.2 polymers are very complex molecules with a large diversity in structure. Therefore in order to derive a constitutive equation a simplified model, which should be able to describe their rheological properties, is needed. Den Toonder [12] and Orlandi [38] performed a direct numerical simulation (DNS) by using a coupling of the extensional viscosity to the local deformation of the flow. This results in an anisotropic stress model. Their results are qualitatively in agreement with experimental observations. However, these models are oversimplified, which means that they do not take into account important properties of the polymers such as, for instance, the possibility of polymers to stretch. Furthermore, in these models the constitutive equation is not Galilean invariant.
In the 1930s, the elastic dumbbell model was introduced (see Bird et al. [5]). With the dumbbell model, a chain-like molecule is represented as two beads connected by an elastic spring. It has been found t h a t the dumbbell model is able to describe the main characteristics of the polymers, i.e. the extensibility (the polymers can stretch and recoil) and the orientability (the polymers can rotate).
Many authors (among them Min et al. [35], [36]) have performed DNS of tiu'bulent flows using the Oldroyd-B or linear Hookean model: the intermolecular force is modeled
as an elastic spring. The disadvantage of the linear spring is t h a t it can be stretched to infinity, leading to an infinite elongational viscosity also at finite elongational rate, which is an unlikely behavior from a physical point of view.
An improved version of this model was used by Massah et al. [34] who studied the coiifigurational changes of single polymer molecules in simple rheological shear and elon-gational flows. They also studied turbulent flows using a FENE (finitely extensible non-linear elastic) model: in this model the intermolecular force (the tension in the spring) is following a nonlinear law (see Bird et al. [5]). Their results show the importance of the stretching of the polymers.
Other authors (Sureshkumar et al. [47], Dimitropoulos et al. [13]) have performed DNS using the FENE-P model, see Bird et al. [6]): the results of these simulations for dilute polymer solutions predict a drag reduction and a change of the turbulence statistics that is in good agreement with experiments. The FENE-P model has also been used by Ptasiiiski et al. [40], to perform a DNS of a channel flow in the high drag reduction regime. The results of this DNS, which provides a good description of the polymer effects near the maximum drag reduction asymptote, are found to be in good agreement with the experimental results,Ptasinski et al. [41]. The work of Ptasinski et al. [41] is used as starting point for the present research.
1.7 Outline of this thesis
The main objective of this thesis is to better understand the mechanism responsible for drag reduction and its implications for practical applications. To gain a better imder-standing we use direct numerical simulations.
The thesis is based on papers which have been published or submitted for publication in the past four years, for this reason each chapter is self-contained and can be read independently from each other.
In chapter 2 we present the mathematical background of a turbulent flow with poly-mers. The governing equations for the flow and the constitutive equations for the polymer model are discussed in the first part of this chapter. In the second part of the chapter we discuss the numerical procedures used in the DNS. The spatial and temporal dis-cretization are presented together with the simulation techniques. We also describe the implementation of this method on a parallel supercomputer.
In Chapter 3 we focus on the drag reduction in turbulent rotating flows. Up to now no direct numerical simiüation has been carried out in a rotating channel flow. Many application of drag reducing flows can be affected by the presence of rotation (f.i. rotating machinery or curvature of the pipes). Our objective is to perform direct numerical simula-tion of these types of flows in order to study the combined effect of polymers and rotasimula-tion on turbulent flows. The polymers are modeled using th(> FENE-P model. In the flrst part of the chapter we explain the mathematical model and the numerical techniques used for the DNS. The results of the DNS for channel flow with rotation are presented, focusing on the flow statistics, the spectrum of turbulence, the polymer conformation tensor and
the kinetic energy budget.
In chapter 4 we consider the interaction between turbulent velocity fluctuations and polymers in terms of the elastic energy that can be stored in the polymers. We also study the changes in the local topology in Newtonian and viscoeiastic flows.
In chapter 5 we show some preliminary results of a method we have developed to alleviate numerical problems we encountered using the FENE-P model. With this method we solve the constitutive equation for the polymers using the logarithm of the polymer conformation tensor. A direct numerical simulation has been carried out in order to test the method in turbulent channel flow.
Chapter 2
Mathematical background and
numerical method
In this chapter the equations for the flow and the constitutive relations for the polymers are given. The polymers are modeled with the help of an elastic dumbbell model. The kinetic theory is used to derive the constitutive equation. Subsequently, we discuss how the constitutive equation is solved coupled with the equations for conservation of mass and momentum. Finally, the numerical method is presented together with the simulation parameters.
2.1 I n t r o d u c t i o n
We perform Direct Numerical Simulation of drag reducing flows. We present the mathe-matical backgroimd and the numerical methods. The outline of the chapter is as follows: the governing equations for the flow field, polymer modeling and constitutive relations are derived in sections 2.2 and 2.3. In section 2.4 the non-dimensional equations are given. In section 2.5 we discuss the computational procedures used for the DNS. The spatial and temporal discretization are described together with th(> method for the time integration of the constitutive equations. Finally, in section 2.6 we discuss the parallel implementation of the computer model.
2.2 Governing equations for t h e flow field
The etjuations which describe the spatial and temporal evolution of the motion of an isothermal and incompressible fluid are the equations for conservation of mass and mo-mentum. These equations can be found in standard textbooks (f.i. Batchelor [3]). The equation for the cons(>rvation of mass and momentiun read:
V - w = 0, (2.1) 11
Du Ou
^~Dt = ''~dt ^ ''•" • ^ " " ^ ^ ' ' + V • r + ƒ . (2.2)
Here p is the density of the fluid, u is tlie velocity vector, p is the pressure, T is the stress tensor, ƒ is a body force and ; ^ = ^ + M • V is the material derivative.
The stress tensor will be expressed as a sinn of a solvent stress r ' " ' and a polymer stress T^P^:
T = r''-' + r ' P ' = ,ir/o ( V M + {Vuf^ + r*") (2.3) where /i (/^ < 1) is the ratio of the solvent viscosity ?],, and the total viscosity r/o- Using
the stress tensor, the equation for conservation of momenttnn reads:
P-F^ = P-^ +pu-Vu= - V p + fhioV'u + V • r(P' + ƒ . (2.4) üf at
To close the set of equations given above we need a constitutive equation for the stress tensor. This equation should be a function of the polymer configuration and the flow field.
2.3 Constitutive equations
2.3.1 Elastic dumbbell model
The elastic dmnbbell is a simplified model for a polymer niolecul«>: two beads, each with mass m, connected by an elastic spring. Let us assume that the beads are labeled " 1 " and "2" and the connector vector which describes the configuration and the orientation of the polymer molecule is given by Q = r i — r2 with r, {i = 1, 2) the location vectors of the two beads, as shown in figure 2.1. A detailed descrij)tion of the elastic dmnbbell model can be fomid in Bird et al. [5]. The configm-ational distribution func^tion needed to compute the stress tensor can be obtained as product of the configuration-space distribution func-tion and the velocity-space distribufunc-tion funcfunc-tion. It is generally assumed in the kinetic theory that the velocity-space distribution is Maxwellian round the solution velocity at the center of mass of the dumbbell and that the configuration-space distribution function is independent of the location of the center of mass of the polymer molecule.
Furthermore, it is assumed that the inertia of the dumbbell can be neglected. This can be justified, because of low mass and the slow motion of the beads. The components of the total force acting on each bead are:
• A hydrodynamics force F*''', which is due to the resistance experienced by the beads during the motion through the fluid. The hydrodynamic force tends to orient and distort the polymer molecules and it can be expressed by Stokes law, see Bird et al.
•
A Brownian force F^''' (hie to the thermal fluctuations in the liquid. This forceFigure 2.1: Schematic representation of a dumbbell, r j . r-i. Q.
• An intramolecular force F^''' which is a result of the internal tension of the spring. This force tends to restore the molecule to its original sliape.
• External forces i^'**' such as gravitational or electromagnetic forces.
By constructing a force balance on each bead and subsequently calculating the total force on the dumbbell, an equation for the motion of the cormector vector can be obtained. This equation reads:
^ dt-^ IQ]=K Q •2kT 0
InV-
r'n
?('') F\ (c) (2.5)111 which K is the gradient of the velocity vector. C is the isotropic friction coefhcient of the beads, k is the Boltzmann constant, T is the absolute temperature, and ij> (Q, t) is the normalized configuration-space distribution function. The operator [• • -J denotes averaging over the velocity-space distribution function. When the same external forces act on both beads, like in our case, the last term drops out.
The continuity equation for the configuration-space distribution function is given by:
di)_
dt
_d_
By substituting the expression of [Q\il> given by (2.5) into (2.6) we can derive an equation for 4-' (Q-O- This equation is the so-called diffusion equation :
-^ = - ^ - ( f c - Q r - — ^ ^ . - - F ' " . j . (2.7)
By nmltiplying this equation (2.7) by the second order tensor QQ and subsequently integrating it over the configuration space (denoted by the operator < ... > ) we obtain the equation for change of {QQ):
^ ^ ^ = - u • V ( Q Q ) + V u • (QQ) + (QQ) • (Vu)'^ + ^ I - ^{QF^^^).. (2.8) For a system in equilibrium the equation (2.8) gives:
{QF^''') = kTI. (2.9)
2.3.2 Polymer stress tensor
In order to get an expression for the polymer stress tensor, we follow the theory j^roposed by Kramers [28]. The total polymer stress TT is given by the sum of three contributions:
• An intramolecular contribution (force of tension or compression in the spring), given by:
n^^-'^ = viQF^'--^). (2.10)
• A contribtition due to the external forc:es acting on the beads, given by:
7r(P-) = l „ ( F ( ^ ) _ F ( ; ) ) . (2.11) This contribution is zero in oiu' case, because the same forces are acting on both
beads.
• A contribution given from the motitm of the beads given by:
^(p,b) ^ ^2nkTI. (2.12)
The total polymer stress tensor is thus given by:
•K = n{QF'^'^) - 2iikTI . (2.13) Contribution from intramolecular forces Contribution from beads motion
The polymer stress tensor can be expressed in terms of r ' ^ ' by subtracting from equa-tion (2.13). the corresponding expression at equilibrium (i.e. TT^^ = n{QF^'~'),.g~2nkTI —
—nkTI, where equation (2.9) has been used):
r(p) = -nkTI + Ji{QF^'''>). (2.14) This expression is called the Kramers form of the stress tensor.
2.3.3 F E N E - P model
Until now nothing has been said about the intramolecular force given by the connec;tor force of the spring. By specifying the expression for this force, a constitutive equation can be derived from equation (2.8). The constitutive equation should be solved together with the equations of conservation of mass and momentum. In this thesis, a FENE-P (Finitely Extensible Non-linear Elastic with Peterlin approximation) model has been used. The FENE-P model was proposed for the first time in 1980 by Bird et al. [6].
The connector force for the FENE-P model is given by:
F('^) HQ
i'{Q')/Ql'
(2.15)In which H is the elastic constant of the spring, (Q^) = tT{QQ) is the length of the dumbbell and Qo is the maximum length for the spring. A spring with this force will have a linear (Hookean) behaviour for small deformation and a non-linear behavior (stiffer) as the spring is extended. It is coiivenic^nt to scale the dumbbell vector with the length scale Qfqui = \/kT/H such that Q = Q/^/kT/H. It is furthermore useful to introduce a length parameter b = HQ^/kT, a relaxation time A = C,/^H and a jjolymer conformation tensor: c = ( Q Q ) . By definition c is symmetric and positive definite.
The equations (2.14) and (2.8) combined with the dumbbell force (2.15) leads to the following set of equations:
and with: dc T 1 — = - u • V c + V u • c -h c • ( V u ) +-{l~ fc) ot A -(P) = -nkT (I - fc)
/ H l - ^ t r c
(2.16) (2.17) (2.18) From equation (2.16) and (2.18) it follows that for a system that is in equilibrium (i.e. with V • M = 0 and dc/dt = 0):36
tr c„,
h + 3'
(2.19)Using the fact that tr Cmax = f> we can obtain:
tr c„,
b + 3
(2.20)
The length parameter b thus deternunes the extensibility of the polymer, i.e. the ratio of the maximum length to the equilibrium length of the polymers.
The derived model is characterised by three parameters: the relaxation time A. the extensibility parameter b oc {Qo/Qequi) and the concentration of the polymer n. The
Flow
z rV^
\
J_^
^
^
--^
^ \
*• AT
= =W"^ /
^^
^
^ ^
^.^
^ y ^Figure 2.2: Schematic representation of the channel used in the DNS.
relation between /3 and the parameter nkT (determined by the polymer concentration) is used in order to compute the polymer stress in (2.17) and reads according to [52]:
^v = (1 - f'5)''to ^kT\
ft + 3'
(2.21)2.4 Dimensionless equations
In figure 2.2 we present the geometry used in the direc:t immcrical sinmlations: a channel which consists of two solid walls positioned at a distance h from each other. The flow in the channel is driven by a constant pressure gradient | ^ . This pressure gradient must be balanced (on average) by the wall shear stress, r^^ at the channel walls: r „ = | | j . For the non-dimensionalization of the equations it is useful to introduce the so-called friction velocity. This velocity is defined as: u„ = \JT^,Ip. The equations are scaled with this friction velocity and the characteristic length scale {h). The dimensionless variables are defined as x = x/h, ü = u/u» and i = tu,/h. with h the height of the channel. The pressure and the polymer stress are made dimensionless with r „ = /w.J.
This leads to the following set of equations:
V - - ü = 0, (2.22)
dü
+ ü • V u - v p + T ^ v ^ ü + V • f(p' + ƒ . / t e .
(2.23)
where Tïe, is the Reynolds nmnber i?e» = phut/r]o.
The dimensionless evohition equation for the conformation tensor becomes (the tensor c is already dimensionless):
Be 1
— = - Ü • V c + VÜ • c + c • ( V ü ) ^ + ^ (I - / c )
dt X (2.24)
In which A is the dimensionless relaxation of the polymer given by A = Xa^/h. The relax-ation time A can be expressed in terms of a dimensionless ratio between the characteristic
time-scale of the polymer and a the characteristic time-sc:ale of the flow. This ratio is known as the Weissenberg number:
Wi, = ^ = ^ = XRe,. (2.25) Vo Vo
In the sequel of this chapter we will omit the tilde symbol and all equations will be considered in their dimcnsionless form only.
2.5 Numerical approach
2.5.1 Simulations techniques
The system of equations (2.22), (2.23) and (2.24) is discretised on a three-dimensional grid and integrated in time. The flow geometry is shown in figure 2.2. The equations are defined on an orthogonal Cartesian coordinate system (x, y, z) where .r—, y— and z— directions represent the streamwise, spanwise and wall-normal directions, respectively. The corresponding components of the velocity vector are:
u=(u,iKw) (2.26)
and the corresponding components of the jjolymer conformation tensor are:
(
^'xx ^'xy ^'Xz \Cxy Cyy Cy^ I . (2-27)
^xz ^yz ^zz /
Spatial d i s c r e t i s a t i o n
The spatial discretization of th(> governing equations is ])erformed with a pseudo-spectral method in streamwise (.r) and spanwise (Y/) directions. The wall-normal direction is discretised with standard second-order staggered finite differences. With the pseudo-spectral nu>thod a high accuracy can be achieved, without the disadvantage of a too high number of operations/grid points required. A detailed description of (pseudo)spectral methods can be found in Canuto vl al. [9].
In the pseudo-spectral method the physical variables are exi)and(xl in terms of finite Fourier series, which for a generic function ƒ reads:
/(^) = ^ A e ' ' " (2.28)
k
the derivatives in the x~ and y—directions are computed by transforming the various terms (like the advective, diffusive or polymer terms) into the spectral space. The trans-formation is accomplished with helj) of a one-dimensional Fast Fourier Transform (FFT).
dz
X
Figure 2.3: Schematic representation of a grid cell of the staggered grid.
In spectral space first order derivatives are computed by multiplying the Fourier cients with ik. Second order derivatives are calculated by multiplying the Fourier coeffi-cients with —fc^, where k denotes the wave number. The modified Fourier coefficoeffi-cients arc then used in the transformation back to the physical space.
As mentioned before, we use a second order staggered finite difference for the wall normal direction (2). The reason for this is that implementation of the pseudo-spectral method with non-periodic boundary conditions is (>xtrcmely complicated. The following staggered scheme is used: the pressure and the x— and y—velocity components are de-termined at the centers of the grid cells, i.e. at (n — ^)^z. The 2—velocity component is determined at the faces of the cells, i.e. at nAz, with n an integer number. For the conformation tensor c and the polymer stress tensor r'-''' the xz- and yz-components are defined at the cell face and the other four components of the conformation tensor at the cell center. In the case that a variable is needed at a location where it is not defined, its value is determined using linear interpolation. With a staggered grid, all gradients in the wall-normal direction are determined over the smallest possible grid spacing and thus an accurate discretization of the gradients is provided in this direction. The staggered grid is shown in figure 2.3.
Most of the simulations are carried out in a so-called miinnial flow >nnt (MFU). The turbulent flow in such a minimal flow >mit has been studied in detail by Jimenez & Moin [23]. The concept of the MFU is to create a minimal size flow geometry in which the essential turbulent characteristics, such as dynamics and morphology resemble the turbulence properties of a full-scale turbulent channel flow. We have chosen for a MFU,
because this minimises the amount of grid points needed for comijutation of th(; flow field. Jimenez & Moin [23] found that the requirements for the MFU are a spanwise extent of at least 100 wall units (where a wall unit is a diniensionless distance from the wall defined as z~^ = zUtp/y) and a streamwise extent of 250 — 350 wall units. From previous work, e.g. Ptasinski et al. [40] we known that turbulent structures tend to be larger for polymer flow than for the Newtonian flow with the same Reynolds number. Therefore, the comijutational domain is chosen larger than these minimal requirements. Based on the Reynolds number Rct = pUfH/r] = 360 for which the Newtonian flow simulation is carried out, the dimensions of our channel are chosen as 1.5ft x h x h for the streamwise (x), spanwise (y) and wall-normal (2) directions respectively. The number of grid points has been chosen equal to 48 x 32 x 100 for the streamwise (x), spanwise (y) and wall-normal (2) directions respectively.
Some preliminary simulations for a rotating channel showed t h a t these values are not sufficient to resolve all the scales of the flows. This is due to the fact t h a t rotation generates a large scale secondary motion. The length scale of the secondary motion is typically of order of the channel height h. Therefore, in the rotating case we have to use a larger computational box. We have chosen for the following dimensions 9h x 6h x h in the streamwise (x), spanwise (y) and wall-normal (z) directions respectively. The number of grid points in the simulation is equal to 96 x 128 x 100 in the streamwise (x). spanwise
(y) and wall-normal (z) directions respectively.
T e m p o r a l d i s c r e t i s a t i o n
For the time integration we use the explicit second-order Adams-Bashforth method (Caimto
et al. [9]). This method is used for both advective and diffusive terms. Conservation of
mass is enforced with a pressure correction method. In the first step of this method an intermediate velocity u* is computed as follows:
where F= ~uVu + -j^y^u + V • r''"' + ƒ .
The intermediate velocity u* computed with this equation does not necessarily satisfy the divergence free constraint, (2.22). Therefore in a second step the intermediate velocity
u* is corrected with th(> pressure equation in order to satisfy continuity:
^ / = - V / / ' + ' (2,30)
where the pressure at the new time level is determined by applying the divergence operator to this equation and enforcing the flow to be divergence free at the new time level:
V • u " + ' = V • u ' - A / ^ V V + ' = 0. (2.31) This leads to a Poisson equation for the pressure, which is solved in a way consistent
out in Fourier space and in the z—direction the finite difference technique is used. This leads to the application of F F T s in the homogeneous directions and tridiagonal matrices in the wall-normal direction. Finally, equation (2.30) gives the velocity at the new time level u " + i .
The time step is determined using the Courant criterion, which is a criterion to avoid numerical instabilities;
At C
I A x I + I A a I ~'~ I A z I "*" fie, I (Aa-)2
+
( A y ) -+
(Az)2(2.32)
where C is set to 0.3. see [40].
The discretiaed equation for the polymer conformation tensor is solved simultane-ously with the contiiniity and the momentum ecjuations. We use a second-order Adams-Bashforth scheme for the advection and the deformation due to the flow and a second-order Adams-Moulton (or Crank-Nicholson) scheme for the FENE-P force. This force is needed to prevent the dumbbells from exceeding their maximum length, i.e. tr c < b. The Adams-Ivloulton scheme is an inii)licit method, which means that some of the quantities needed for the solution of the polymer ecjuation and conformation tensor iimst be computed at the new time level t(„+i). The polymer conformation tensor at the time level (n -|- 1) is computed with help of:
r.n+1 At
where G and H ar(> given by:
2 G'- + \H" + \H^'+'
G = - u V c + V u c + c. - ( V M ) ^ + - I
A
The resulting expression for c at the new time k^vel then becomes:
(2.33) (2.34) (2.35) 1-^
At
1 2A1 }trc"+i I)c" + A^ ( -G"
2 ^ " - + - H "
(2.36)By taking the trace of this etiuation and rearranging the terms thv. following equation can be obtained:
.2
r At.
(trc"-* - ( 1 -h — 6 -F rhs ) trc"+' + rhs = 0 with r h s •trc"
+ At I -trG"
: t r G " -2 : t r i f " (2.37) (2.38) Equation (2.37) is a quadratic equation for tr c"+^ which can be solved directly. With the value of tr c""'"^ all components of c at the new time level are calculated using equation(2.36). Finally, by using equation (2.14) the polymer stress at the new time level is obtained.
The dimensionless time step is determined with help of equation (2.32) and is set equal to At = 2 • 10^^ for the Newtonian case. For the viscoelastic cases the time step is set to halve this value, which was needed in order to keep the simulations stable. For the rotating cases a smaller time step is needed. This is due to the slight increase in maximum axial and wall normal velocity. For the rotating case a dimensionless time step of A^ = 1 • 10"'' (Newtonian) and At = 0.75 • lO"'' (polymer) has been used.
To avoid mmiorical instabilities a small artificial diffusive term a/ (uth) V^c is added to the equation for the conformation tensor (2.24). The reason for this is that small numerical errors during the time-integration of the constitutive equation can lead to a loss of the positive definiteness of the polymer conformation tensor. The polymer conformation tensor is the second moment of the end-to-end distribution vector, and its trace is related with the length of the jjolymers, which can not become negative.
For most polymer models it has been shown (Hulsen [20]) that with exact integration of the constitutive equation the polymer conformation tensor will remain positive definite as long it was positive (l(>finite as initial condition. However, numerical errors due to time integration and spatial discretisation of the convection term can lead to a loss of the positive definiteness.
The effect of the additional diffusive term should decrease when the number of grid points is increased and the time step decreased. Furthermore, it should not have a large influence on the macroscopic flow parameters like the velocity or stress profiles. The value chosen in this paper has been found not to affect the results significantly. An analysis car-ried out by Ptasinski et al. [40] showed that for two different direct numerical simulations of the same flow, with different values of the artificial diffusivity ( a / {uji) = 0.02 and
a/ {Uth) = 0.012) the velocity and stress profiles did not show significant differences.
Fur-thermore, an analysis on the effect of this term has been carried out by Sureshkumar and Beris [46]. They showed that small values of the artificial diflï'usivity do not significantly alten- the results of the calculations.
S t a t i s t i c a l p r o c e d u r e
The statistical i)rocedure used for the analysis of the Direct Numerical Simulations is as follows: Starting from a turbulent velocity field, a DNS is carried out for a Newtonian fluid, i.e. without FENE-P model, until a statistically steady state was reached. This takes approximately 2QT. where the timescale T is defined as h/u,. From this DNS a fully developed velocity field is obtained. The velocity field is used as initial condition for the simulation with the FENE-P model. The DNS is then carried out without coupling the fields of u and c (i.e. (5=1), until a statistically steady state for c was reached. This step takes approximately lOT. At this point a coui^led simulation is started by setting (3 at the value chosen for the DNS and the system of etination is then integrated in time until again a steady state is reached (normally this takes about 10 — 12T'). The simulations are then continued for another lOT" during which data fields are collected for post processing.
^ ^
ït^^aa^
1 1 1 j -1 - J . 1 I 1 - 4'^i^^rl
._j._j-JFigure 2.4: Data distribution over tiic nodes. The domain is subdivided in tiie .r-directiou. For the execution of FFT's in the x-direction, which is needed to solve the Poisson equation and to compute the spatial derivatives in the x- directions, the data is redistributed as shown on the right.
1 Nmnber of proc
1
2
4
8
16
[ 32
essors
CPU time (seconds) ~\
44
24
13
8
6
5 J
Table 2.1: Comparison of CPU time needed for one time step as a function of the numl)er of processors.
The time separation between the fields is O.IT. This time separation is large compared to the time scale of turbulence. This implies that the data fields (^an be considered as statistically independent.
During the post processing, data fields over the last 5T are used to compute various statistics. The statistics are obtained by spatial averaging over streamwise and spanwise directions and temporal (ensemble) averaging over the data fields. We have chosen to compute statistics using only the last ST of our sinmlations. This has been done to ensure that the fields are statistically in a steady state. (Statistics computed over the last 2.5T and lOT did not show significant differences, this is an indication that the amoiuit of samples is adecjuate).
2.6 Parallel implementation
The code used for the DNS originally has been written for serial and/or vector computers. The time needed for the solution of the flow ecjuations and the FENE-P moflel is consid-erable. To reduce this time the algorithm has been implemented on a parallel comi)uter. With a parallel implementation, multii)le processors are used: the algorithm is broken into discrete series of instructions, and instructions from each part are executed simul-taneously on different i)rocessors. The Message Passing Interface (MPI) library is used to parallelise the numerical model using a domain decomposition method. The computa-tional domain is divided into subdoniains, see figure 2.4. These subdomains are divided
Figure 2.5: Mean velocity profile as a function of the distance to the wall in wall units. Plotted are the DNS-results (continuous), the results of Kim et al. [27] (dashed) and the empirical curve u+ = 2.5 In 2+ + 5.5 (dashed-dotted).
i"igure 2.6: The root mean scjuarc of the velocity fluctuations for all three components as a function of ihstance to tlie wall. Potted are the DNS-results (continuous) and the results of Kim et al. [27] (dashed).
over the different nodes of the parallel eoniputer. We have subdivided the original domain only in the x-direction, in the tj- and ^-directions all grid points are available on the same processor.
For an efficient computation of the one-dimensional F F T ' s that are needed for the solution of the Poisson equation and to determine the spatial derivatives, the values of the velocity components, pressure and polymer conformation tensor at all grid points in that direction nnist be available on the processor that performs the F F T . The spatial derivative in the y-direction can be computed on the original data distribution. The F F T ' s in the x-direction are computed after a redistribution of the data over the processors, as shown on the right of Figure 2.4. The results are then redistributed back to the original data distribution and the tridiagonal matrices in the Poisson solver are solved on this distribution.
The differences in terms of CPU time for the different direct numerical simulations carried out with a different amount of processors is shown in table 2.1. The scaling is less than linear because of the time needed for the connnunication between the processors. For our parallel computations we chosen a number of processors equal to 16.
In figure 2.5 we show the mean streaniwise velocity profile (solid line) for a Newtonian non-rotating flow. The empirical curve w+ = 2.5 In 2+ -1-5.5 (dash-dotted lino) and the results of Kim et al. [27] (dashed lines) are also shown in the figure. The velocity profile agrees well with the empirical relations and the results of Kim et al. [27]. Figure 2.6 shows the root mean scjuare of the velocity fluctuations for all three components as a fimction of distance to the wall. Again, we find that the results agree very well with the results of [27]. the small differences in the streaniwise velocity fiuctuations are likely due to the small dimensions of the channel.
Chapter 3
DNS of drag reduction in a rotating
turbulent flow
In this chapter the direct numerical simulation of a rotating turbulent liquid flow with polymers is presented. Polymers are modeled using the elastic dumbbell theory, from which a constitutive equation for the polymer stress is derived. The constitutive equation is solved simultaneously with the Navier-Stokes equations in a rotating frame of reference. The in-fluence of rotation is analysed for different rotation numbers varying from 0 (no rotation)
to 1.5 (strong rotation rate). The mean velocity profile, stress profile and polymer con-form.ation tensor have been com.puted. Also the spectrum of turbulence has been analysed.
The results for a rotating flow of a Newtonian fluid are in good agreement with literature data. The results for polymer flows .show that a decrease in the frictional drag occurs also for a rotating flow. However, this drag reduction can be considerably smaller than for th,e rotating case. The drag reduction is, for instance, reduced from 26% for the non-rotating case to 11% for the case with a strong rotation. This effect is in agreement with (recently published) experimental re.iults on drag reduction in curved pipes. The influence of rotation on the frictional pressure drop is stronger for a polymer flow than for the flow of a Newtonian fluid.
3.1 Introduction
The discovery that the addition of a small aniouut of polymer to a turbulent pijie or channel flow leads to a significant reduction of the frictional drag dates back to Toms in 1949, [49]. It has been obs(n'V(>d exi)erimentally that even an addition of a few parts per million of polymers could reduc(> the pressure drop for a turbident flow up to 80%. In the past years turbulent drag reduc^tion has been used for many practical ajjplications.
The study of turbulent drag reduction involves two important fields in fluid mechanics, namely turbulence and rheology. Because of their small size, polymers are mostly active on the small length scales. However, they are also abk> to modify the larger scales of the flow, especially in the buffer layer. Although the phenomenon has been studied extensively,
no generally accepted explanation has been found. Three theories have been proposed in order to explain the mechanisms responsible for drag reduction. The first one was proposed by Lumley [30], who suggested that the stretching of the polymer molecules by the flow increases the effective (extensional) viscosity and the thickness of the buffer layer, leading to drag reduction. In the second theory, proposed by de Gennes [11], it is supposed that drag reduction is caused by the elasticity of polymer molecules, by which energy at the small scales of turbulence is adsorbed and subsequently radiated away in the form of shear waves. More recently Procaccia et al. [39] showed that the polymers cause an effective viscosity which is a linear function of the distance from the wall and numerical simulations by De Angelis et al. [10] showed that such a linearly varying viscosity leads to drag reduction.
During the past years many numerical simulations on drag reduction by polymers have been carried out with different kinds of polymer models. Direct nmnerical simulations (DNS) have become possible with the development of accurate and efficient spectral meth-ods, advanced viscoelastic models and increased computer capabilities. The advantage of numerical simulations is that it provides information about quantities that are difficult to measure experimentally. Moreover, it is possible to study separately the effects of various polymer properties like elasticity and concentration.
Several authors (Sureshkmnar et al. [47], Dimitropoulos et al. [13]) have performed DNS using the FENE-P (finitely extensible non-linear elastic with Peterlin approxima-tion) model. The results of these simulations for dilute polymer solutions predict a drag reduction and a change in the turbulence statistics, that are in good agreement with ex-periments. The FENE-P model was also used by Ptasinski et al. [40] to perform a DNS of a channel flow in the high-drag-reduction regime. The results of this DNS study, which provides a description of the polymer effects near the maximum drag reduction asymp-tote, were found to be in good agreement with the results of experiments (Ptasinski et al.
[41]).
Turbulent flows in a rotating reference frame are of interest for a considerable number of industrial applications which use rotating machinery such as turbines, pumps and compressors. In 1969, Bradshaw [7] showed the analogy between a shear flow rotating about a spanwise axis and a shear flow over a curved surface. In laminar and turbulent flows subject to rotation the resulting flow patterns depend on the orientation of the vorticity vector. The rotation influences both the mean flow and the turbulence. DNS of a turbulent flow subject to spanwise rotation has been carried out by several authors (see for example Kristoffersen & Andersson [29], Nagano & Hattori [37], Salshi & Camboii [43], Johnston et al. [24], Kim [26]) for Newtoiuan fluids.
Applications of drag reduction by polymer addition in a turbide^nt pipe flow can be influenced by a rotation, for example by rotating machinery in a pipe line or by a curvature of the pipe. Experiments on drag reduction in curved pipes have been carried out by Shah
et al. [45], which showed t h a t drag reduction is (negatively) influenced by a pipe curvature.
No numerical simulation of drag reduction in a rotating or curved pipe or channel has been performed imtil now. Therefore, we performed such DNS calculation for a rotating turbulent channel flow and paid particular attention to the influence of the rotation rate
Pressure side
Suction side
Figure 3.1: Schematic representation of the geometry used in the DNS.
on the level of drag reduction.
This chapter is organised as follows. The first part provides the mathematical back-ground. In section 3.2 we describe the governing equations for the flow and for the polymer model. The numerical techniques which are used in the DNS, the computational proce-dures, the parameters chosen for the polymers, the flow and the geometry of the system are discussed in section 3.3. In section 3.4 we present the simulation results. In it the drag reduction, the turbulence statistics and the polymer conformation tensor are discussed. The conclusions are presented in section 3.5.
3.2 M a t h e m a t i c a l model
3.2.1 Governing equations for the flow field
We consider an incompressible Newtonian flow in which polymers are dissolved. The solution is assumed to be sufficiently dilute, so that polymer molecules do not interact with each other even when they are fully stretched. The flow is considered to be isothermal and subject to a rotation with constant angular velocity il = (0, H, 0) about the ,y-axis, m the coordinate system as shown in figure 3.1. SI is assumed to be positive (fi > 0). The dynamics of this system can be described using the equations of conservation of mass and momentum, which can be found in standard text books on fluid mechanics (see for instance Batchelor [3]):
V - u = 0, (3.1)
du
p— + pu • V u -f- 2pn X u = -Vpeff + V • r . (3.2)
the last term on the left-hand-side of the momentum equation (3.2) is the Coriolis force due to the rotation of the system. The centrifugal force has been included in the effective pressure p^jj. In these equations t is the time, p the density, u the velocity vector and
Figure 3.2: Schematic representation of a dumbbell.
solvent ( r ' " ' ) and a viscoelastic part due to the polymer ( r ' ^ ' ) :
r = r'^' + T(P' = l3rio ( V u + (Vuf) + r ' ^ ' (3.3) where /?(= r/s/r/o) is the ratio of the solvent viscosity rjs and the total shear viscosity rjo of
the solution at zero shear rate. Using (3.3) equation (3.2) can be written in the following way:
p ^ + pn • V u + 2pn X u = -Vpeff + ,M)V^u + V • T ' P ' . (3.4)
3.2.2 Constitutive equation
To derive the exjjression for the polymer stress we represent polymer molecules as dumb-bells: two spherical beads connected by a massless spring. A (hmibbell is specified by a connector vector Q, which describes the orientation and the internal configuration of the polymer as illustrated in figure 3.2. This is a simijlified model of a polymer, but it has been found to realistically describe the two essential characteristics of polymers namely the orientability and extensibility. A detailed description of this model can be found in Bird et al. [5].
An evolution equation for the conformation tensor (QQ) can be derived by applying a force balance to each of the two beads separately and by subtracting these two balances. The forces acting on the beads are the hydrodynamics force, the Brownian force and the intramolecular force. It is assumed t h a t the beads are slowly moving in the flow, so that inertia can be neglected. The rotation has no effect on the dumbbell connector vector, because the same force is acting on both beads. The resulting equation is given by:
^ ^ = - u • V ( Q Q ) + V u • ( Q Q ) + (QQ) • (Vu)'^ + ~ I - ^ ( Q F C ^ ) ) (3.5)
where k is the Boltzmami constant, T the absolute temperature, C, the isotropic fric-tion coefficient of the beads and F''^'(Q) the connector force in the spring between the two beads. The brackets (• • •) denote an average over all possible conformations of the polymers.
