Master Thesis Advanced Techniques for Investigating Structures in Computational Fluid Dynamics

(1)

Master Thesis

Advanced Techniques for Investigating Structures in Computational Fluid Dynamics

Author:

Rasmus Ellebæk Christiansen, s072162

Supervisors:

Professor Morten Brøns

Associate Professor Allan P. Engsig-Karup

DTU Compute - Department of Applied Mathematics and Computer Science.

29th July 2013

(2)

The problems of incompressible fluid flow past a cylinder in free flow and past a cylinder near a moving wall in two dimensions are studied numerically. Here the velocity and pressure fields are obtained using a spectral element method based solver. The study is performed for Reynolds numbers in the low end of the periodic shedding regime and is mainly concerned with analysing vortex structures existing behind the cylinder, with main focus on the vortices surviving downstream. A vortex is defined uniquely as an extremum in vorticity and this definition is used to implement an algorithm for identifying and tracing vortices. A way of assuring that periodicity in the numerical solution has been reached using the tracing method is presented.

The Reynolds number, time and ratio of cylinder diameter to distance from the wall are identified as the only independent parameters of the problems. Dynamical systems theory is used to analyse changes in the vortex structures as the independent parameters are varied.

Here, two types of bifurcations are identified and characterised using existing theory. Also, all different structures for vorticity extrema and saddle-points observed in the considered parameter range are identified and described.

The stabilizing effect of the wall is mapped. The effects of the wall and Reynolds number on the shedding frequency, creation point and pathway followed by the vortices are investigated.

The decrease in magnitude of a vortex as it travels through the fluid is investigated and a nearly constant decrease rate is found far downstream.

A small scale investigation of collocation based uncertainty quantification is performed where the Reynolds number is taken to be the input parameter containing uncertainty. First, the cylinder in free flow is treated and the Strouhal number and pressure induced drag are considered as functions of the uncertain Reynolds number. Secondly, the cylinder near the moving wall is treated and the path followed by a vortex downstream and its magnitude as a function of the Reynolds number investigated.

Used theoretical concepts from dynamical systems theory, polynomial approximation theory, the spectral element method and uncertainty quantification are included. So are explanations of the software used and developed for the simulations and post processing readying the data for the final analysis.

Most of the methods used in this work have not been adapted by the industry (yet).

(3)

Abstrakt

To todimensionelle modelproblemer omhandlende en inkompressibel væske der flyder omkring en cylinder hhv i åbent flow og omkring en cylinder nær en væg i bevægelse studeres numer- isk. Hastigheds- og trykfelterne bestemmes ved brug af en spektral element metode baseret løser. Studiet udføres for Reynoldstal i den lave ende af regimet for periodisk hvirveldan- nelse og beskæftiger sig hovedsageligt med analyse af de hvirvelstrukturer, der eksisterer bagved cylinderen med fokus på hvirvlerne der overlever nedstrøms af cylinderen. En hvirvel defineres unikt som et ekstremum i vorticitet, og denne definition bruges til at implementere en algoritme til at identificere og følge hvirvler. En måde til sikring af fuldstændig period- icitet i flowet ved brug af denne metode præsenteres.

Reynoldstallet, tiden og cylinder diameter over cylinder afstand til væggen identificeres som problemernes eneste uafhængige parametre. Dynamisk system teori bruges til at analysere ændringer i hvirvelstrukturerne når de uafhængige parametre varieres. Her identificeres og karakteriseres to typer bifurkationer ved brug af eksisterende teori. Yderligere identificeres og beskrives alle strukturer for vorticitets ekstrema og saddelpunkter, der er observeret under undersøgelserne.

Den stabiliserende effekt af væggen kortlægges. Effekten af væggen og Reynoldstallet på hvirveldannelsesfrekvensen, positionen af hvirvlernes dannelsespunkter, samt vejen disse føl- ger gennem væsken, undersøges. Ydermere undersøges ændringerne i hvirvlernes intensitet når de rejser gennem væsken, og en næsten konstant formindskelsesrate findes når hvivlerne er rejst langt nedenstrøms.

En mindre undersøgelse af kolokationsbaseret usikkerhedskvantificering udføres. Her ant- ages det, at Reynoldstallet er en usikker input parameter. Først behandles cylinderen i åbent flow, hvor Strouhal tallet og det trykinducerede drag betragtes som funktioner af det usikre Reynoldstal. Herefter betragtes cylinderen nær væggen i bevægelse. Vejen, som en hvirvel følger nedenstrøms, samt hvirvlens intensitet, betragtes som funktioner af det usikre Reynoldstal.

Anvendte teoretiske koncepter fra dynamisk system teori, polynomial approksimations teori, spektral element metoden og usikkerhedskvantificering er inkluderet. Det samme er tilfæl- det for en beskrivelse af softwaret anvendt til simulering og efterbehandling af data ved forberedelse til den endelige analyse.

De fleste metoder der anvendes i dette arbejde anvendes (endnu ikke) i industrien.

(4)

(5)

CONTENTS CONTENTS

Introduction

This thesis presents and investigates advanced methods from different mathematical fields and provides the results of using these to analyse a Fluid Dynamics (FD) problem with focus on vortices appearing in the flow. The first field is Dynamical Systems Theory (DST) [1], [2], which is used to identify structures and bifurcations in vortex creation, annihilation and movement patterns. The second is the relatively new field of Uncertainty Quantification (UQ), with focus on the stochastic collocation approach, see e.g. [3] or [4]. All investigations are based on numerical solutions to the Navier-Stokes equation obtained using the Spectral Element Method (SEM) approach. The methods used and their mathematical foundation will be explained in some detail in chapter 2 while an interested reader is encouraged to consult the references for deeper insight. These methods have not (yet) been implemented for use by the industry.

Another aim of this written work is to provide enough information, examples and code to allow the reader to understand and take the steps necessary to reproduce the work presented¹. This approach simultaneously allows the direct adaptation of the methods to analyse a range of other FD problems with relative ease.

As stated above all investigations are based on simulation results obtained using a SEM approach. The SEM solver used is the incompressible Navier-Stokes solver provided in the Nektar++framework [5]. All post processing is done using software written for this project as well as the open source data analysis toolParaview [6].

The structure in this report is presented next. First, a short section on the notation and abbreviations used throughout the report is included for easy reference and to avoid any unnecessary confusion. Chapter 1 presents the model problems investigated, the underlying physics and a number of quantities of interest in the problem. Chapter 2 provides an overview of the theoretical foundation for the methods used from the areas of DST, SEM and UQ respectively along with an algorithm for tracing vorticity extrema and saddle points.

Chapter 3 describes how the problem is discretized to allow for accurate numerical solutions.

Chapter 4 is dedicated to explaining the software developed for the project, the third party software used and the tool-chain setup for performing simulations. Chapter 5 contains a small investigation of the applicability of parallel computing for the solution process using the chosen solver framework.

Chapter 6 provides an overview of the simulations performed throughout the project, how data is visualized and an explanation of the application of the vorticity extrema tracing method.

Chapter 7 presents results of a validation process undergone to verify that results obtained using theNektar++framework can be trusted. Chapter 8 presents the results and analysis for the DST and UQ based investigations. Lastly, Chapter 9 presents a short conclusion and outlook on what future work may be performed to extend the work presented here.

The appendices include the scripts and source code written for the project as well as com- pilation and usage instructions. It also contains a.geoand.mshfile containing an example of the mesh used in the simulations and an example of XML-files to use as input files for theNektar++solver.

1The reader should be aware that the reproduction of the work presented may be very time consuming due to the significant number of large scale simulations needed.

(9)

CONTENTS CONTENTS

Notation

This section is included to avoid unnecessary confusion caused by the notation used in this work. Scalar quantities are denoted by italic letters, vector quantities are denoted by bold letters.

The following quantities have specific letters associated with them.

• Time: Time is denoted t.

• Position: Standard Cartesian coordinates are used unless otherwise stated and they are denoted as, r= (x, y).

• Directional Unit Vectors: Unit vectors denoting the x- and y-directions are de- noted, ˆi and ˆj respectively.

• Velocity: Fluid velocity in Cartesian coordinates is denoted byu= (u, v).

• Vorticity: Vorticity is denotedω and is defined as ω=∇ ×u. In 2D this is a scalar quantity and in 3D it is itself a 3D vector field.

• Pressure: A pressure is denotedp.

• Mass: A mass is denotedm.

• Volume: A volume is denotedV.

• Density: Density is denoted byρ.

• Viscosity: Viscosity is denoted byµ.

• Gab Size: The distance between a cylinder and a wall is denoted, G.

• Cylinder Radius: The radius of the cylinder is denoted, R.

• Cylinder Diameter: The diameter of the cylinder is denoted, D.

• Far Field Velocity: The velocity of the fluid far from the cylinder is denoted,U∞.

• Reference Pressure: The reference pressure far from the cylinder is denoted,p∞.

• Reynolds Number: The Reynolds number is denoted Re and is defined as Re =

ρU∞D µ .

• Critical Reynolds Number 1: Re_crit denotes the Reynolds number value which marks the transition from stationary to instationary flow.

• Critical Reynolds Number 2: Re_crit₂ denotes the Reynolds number value where the two-dimensionality of the flow breaks down and three-dimensional effects appear.

• Set of Parameters: A set of parameters is denotedP.

• Stress Tensor: The stress tensor components are denoted τ_i,j.

• Normal to a Surface: A surface normal vector is denoted ˆn.

• Tangent to a Surface: A surface tangent vector is denoted ˆt.

(10)

• Drag force: Drag on an object is denotedFd.

• Pressure Drag Coefficient: This coefficient is denoted C_D_p = ^p−p₁ ^∞

2ρU∞²

.

• Mean: The statistical mean of a quantity, X is denoted, µ_X.

• Variance: The variance of a quantity, X is denoted, σ_X².

• Uniform Distribution: A uniform distribution on the interval [a, b] is denoted, U(a, b).

• Gaussian Distribution: A normal distribution with meanaand standard deviation bis denoted, N(a, b).

The following is the standard notation for operators:

• Differentiation: Differentiation of a quantity, a, with respect to another quantity, b is denoted:

– Temporal Particle Derivative: The Temporal Particle Derivative of a quantity Qis denoted: ^DQ_Dt.

– Partial Derivative: Partial differentiation is denoted in one of the following ways: _∂b^∂a, ^∂a_∂b orab.

– Absolute Derivative: The absolute derivative is denoted _db^daor ^da_db.

– Spatial Derivative: The derivative vector of all spatial directions is denoted

∇= (_∂x^∂ ,_∂y^∂ ,_∂z^∂)

• Integration: A curve integral is denoted^R ·dκ and a surface integral ^R ·dS where· is a place holder for the integrand.

• L₂-norm: The standardL₂-norm is denoted k · k_L₂ where ·is a place holder for the element to be measured.

• H_w^k-norm: The weighted Sobolev-norm is denoted k · k_Hk

w where · is a place holder for the element to be measured.

• Definition of function at specific values of independent variables: If a function f of a set of parametersP is considered on the subset Ω of the full parameter domain this is denoted asf(P)|_P_∈Ω.

Abbreviations

Abbreviations will be used for names which are encountered often throughout the text.

In general the full name will be used the first time a new concept is referenced with the abbreviation written after the name in parentheses. For easy reference the abbreviations are also listed below,

• Navier-Stokes: NS

• Spectral Element Method: SEM

(11)

CONTENTS CONTENTS

• Dynamical System: DS

• Dynamical Systems Theory: DST

• Uncertainty Quantification: UQ

• Fluid Dynamics: FD

• Computational Fluid Dynamics: CFD

• Partial Differential Equation: PDE

• Boundary Condition: BC

• Degrees Of Freedom: DOF

• Probability Distribution Function: PDF

• generalized Polynomial Chaos: gPC

• Stochastic Collocation Method: SCM Terminology

Here is a short list of terms used freely throughout the text.

• Smoothness: The term p-smoothness or p-smooth denotes a function which has p continuous derivatives. If a function is denoted as infinitely smooth all of its derivatives are continuous.

• Pre Processing: All preliminary work performed for a simulation to enable the solution using a numerical method.

• Post Processing: All work performed on the solution data obtained by solving the problem using a numerical method preparing it for analysis.

• Vortex Shedding: The term vortex shedding refers to the phenomenon that for a sufficiently high Reynolds number any object immersed in fluid will start generating vorticity at its edges and release (shed) the vortices down stream. The frequency with which this shedding occurs is called the shedding frequency.

(12)

1 The Problem and Underlying Physics

This chapter presents the problems investigated in this work, it outlines the purpose of the investigations, and it presents the underlying physics and the quantities of interest. Section 1.1 outlines the purpose of the project. Section 1.2 presents the model problems, their geometry and a set of parameters which defines them. Section 1.3 presents and explains the physics of the problems along with a set of model equations and appropriate boundary conditions. Lastly, section 1.4 describes a number of physical quantities of interest for the model problem.

1.1 The Purpose

This project has two main focuses. Firstly, it seeks to present and illustrate the application of a method for the investigation of vortex creation, annihilation and movement patterns for FD problems, by, among others, the utilization of results from Dynamical Systems Theory (DST). This is done based on the strict definition of a vortex as extrema in vorticity, see equation (1.20). The method allows accurate tracing of vortex creation, annihilation and movement patterns and relating both continuous and instantaneous changes in the patterns to the variation of the defining parameters of the FD problem under consideration.

Secondly, a preliminary investigation of the application and performance of Uncertainty Quantification (UQ) based methods to the model problem. Here, the idea is to determine functional dependences and statistics for quantities like pressure-drag and Strouhal number² as well as for a vorticis movement path through the domain as a function of an uncertain input parameter. The goal here is to investigate whether the method allows the dependency of the quantities of interest on an uncertain input parameter to be determined with good accuracy with only few model problem realizations. This is in contrast to the classic Monte Carlo method which requires a high number of realizations to obtain accurate statistics, which is infeasible for problems that require long simulation times due to its worst case O(M⁻¹²) convergence rate.

The model problem of fluid flow past a cylinder presented in the following section is well investigated in fluid dynamics, both analytically, experimentally and numerically and thus is believed to be well understood. This means that numerical and experimental data for several quantities exists which can be used to verify the simulations performed for this project, see e.g. [7]. The complication of introducing a moving wall near the cylinder has been investigated by Huang and Sung in [8]. Over the course of these investigations different characteristics of the flow has been investigated for different Reynolds numbers³, by Hende- rson in [7] and by Huang and Sung in [8] for different distances to the moving wall. At different parameter values of e.g. the period of shedding base pressure and pressure drag on the cylinder have been measured, calculated and the relationship with the parameters have been recorded⁴.

There are many reasons for why investigating the vortex formation and movement patterns in flows is interesting from an engineering standpoint. A large number of real world problems may potentially be better understood if a better understanding of vortex formation and movement is obtained. Here follows a few examples. The movement of flags as a consequence

2The Strouhal number is a non-dimensionalised frequency and will be defined later.

3The Reynolds number will be defined in section 1.2.

4Definitions of all these quantities follow in section 1.4.

(13)

1.2 The Model Problems 1 THE PROBLEM AND UNDERLYING PHYSICS

of vortex formation. The vortex formation around cables and pylons in bridge construction and the drag and vibrations caused by these vortices. The vortices formed around ship hulls as they travel through the oceans. The vortices formed behind race cars etc. If a better understanding of the influence of the vortices can be obtained it may be possible to using this knowledge to design structures which are more resilient and ships and cars which are capable of moving faster and more efficient due to reduced drag etc.

The use of numerical methods for obtaining this insight is very valuable as the numerical approach is a highly cost efficient way of performing experiments to obtain better understanding of physical processes. If the same knowledge was to be obtained through physical experiments, multiple and in some cases huge set-up’s would have to be constructed which quickly becomes both very time consuming and very expensive.

1.2 The Model Problems

The main problem under consideration is that of an incompressible fluid flowing around an infinitely long cylinder oriented along the z-axis positioned near a moving wall modelled in two dimensions. The problem is investigated in the low Reynolds number regime (Re ∈ [20,300]) after the initial transient flow has disappeared. In this context the Reynolds number is defined as Re = ^ρU^∞_µ^D where ρ is the density of the fluid, U∞ is the far field velocity of the fluid, Dis the cylinder diameter and µis the dynamic viscosity of the fluid.

The fluid in the far field and the wall moves with the same velocity, U∞ = (u∞,0) while the cylinder remains stationary. This leads to no slip boundary conditions on the cylinder of the form: (u, v) = (0,0) and likewise no slip on the wall of the form: (u, v) = (u∞,0). It also leads to a velocity condition at the outflow boundary of ∇(u, v)·n= 0 and a pressure condition of p = 0. The boundary conditions are elaborated on in section 1.3. The second model problem is also a 2D model of a cylinder but this time in free flow. This problem is considered in the Reynolds number regime,Re∈[100,600].

U∞

D

y x

(a)

U∞

D

G y

x

(b)

Figure 1.1: Schematics for 2D-domain containing a cross section of an infinitely long cylinder perpendicular to the flow. (a) Cylinder in free flow, (b) Cylinder near a wall.

First, the cylinder in free flow is considered in order to validate numerical solutions and establish a reference for the flow structure. This model problem is also used for part of the UQ investigations. Then the main model problem is considered by introducing the wall and

(14)

an investigation for a decreasing distance between cylinder and wall and varying Reynolds number Re is performed. Increasing the Reynolds number above a critical value changes the flow from stationary to instationary with a periodic behaviour and thus introduces a time parameter t. Below the critical value of the Reynolds number the flow is stationary and thus, based on the definition of a vortex provided in (1.20), no vortices exist. Schematic drawings of the two model problems are provided in figures 1.1a and 1.1b.

In figure 1.1 U∞= (u∞,0) is the inflow and far field velocity as well as the velocity of the wall. Dis the cylinder diameter andGthe size of the gap between the cylinder and wall. D andGare not of any interest as two separate parameters since the dynamics of the flow does not depend on them individually but only on their ratio. This fact is easily realised as D andGare the only geometric parameters in the system. Thus scalingDby a factorasimply corresponds to shrinkingGby the same factora. This means that the parameter of interest in the following is their ratio, chosen as ^D_G. The choice of ^D_G opposed to _D^G is made based on the same considerations as those given by Huang and Sung in [8]. The considerations are that it is easier to illustrate small gap heights more clearly on a plot and that with this choice the case of a cylinder in free flow corresponds to ^D_G = 0 instead of _D^G → ∞.

From the present discussion three independent parameters for the problem have been identified. Two of the three are input parameters, i.e. Re and ^D_G. As the problem becomes instationary for increasing Reynolds number time emerges as the third parameter. In the following section where the physics and model equations are presented it is seen that these three parameters along with the boundary conditions are enough to describe the problem completely.

(15)

1.2 The Model Problems 1 THE PROBLEM AND UNDERLYING PHYSICS

Different flow regimes: As mentioned the problem of the cylinder in free flow has been investigated thoroughly in several studies. For Reynolds numbers below the maximal value of interest in this work four different regimes for the flow as a function of Reynolds number have been characterised. In [9] Brøns et. al. identify, describe and analyse the first three of these regimes. The regimes may briefly be summarised as,

• Attached Flow: Re . 5. In this regime the flow is steady and attached to the cylinder. That is, the fluid passes over the cylinder and continues downstream without any circulation. Here, the fluid has a point of attachment on the upstream side of the cylinder and a point of detachment on the downstream side of the cylinder. This flow is sketched in figure 1.2a.

• Steady-Seperation: 5.Re≤Re_crit. In this regime the flow is still steady, however a bauble of recirculation has appeared behind the cylinder. That is, a steady recirculation of fluid is confined to the backside of the cylinder. Here, the fluid has a point of attachment on the upstream side of the cylinder and two points of detachment and a point of attachment on the downstream side of the cylinder. This flow is sketched in figure 1.2b.

• Periodic Shedding: Re > Recrit. In this regime the flow is no longer steady but instead it exhibits a periodic behaviour. This periodic behaviour consists of a periodic shedding of vortices from the backside of the cylinder. The vortices travel downstream creating a so-called vortex train behind the cylinder. The position of the points of detachment on the downstream side of the cylinder fluctuate in time and the point of attachment no longer exists. This flow is sketched in figure 1.2c.

• Breakdown of the Two Dimensionality in the Flow: For the cylinder in free flow, Henderson [10] identifies Re_crit₂ ≈190 to be the Reynolds number at which the two-dimensionality of the flow breaks down. Beyond Recrit₂ the flow exhibits three- dimensional behaviour and the 2D simulations are no longer enough to capture the physics completely.

(16)

(a) Sketch of streamlines for flow in the

attached flow regime. Re.5. (b) Sketch of streamlines for flow in the steady separation flow regime. Fluid is circulating in two closed baubles behind the cylinder. 5.Re≤Re_crit.

(c) Sketch of periodic shedding where vortices are shed from the downstream side of the cylinder.

Re > Re_crit

Figure 1.2: Sketches of different regimes of fluid flow in the Reynolds number range, Re∈]0,600[ for a cylinder in free flow.

There is no reason to expect that the different regimes should not exist when introducing a moving wall near the cylinder, although the Re values at which the flow transitions from one regime to another may be expected to vary. This variation has been investigated for the critical value for periodic shedding, when introducing and moving the wall closer to the cylinder.

An important note is that from the definition of a vortex used in this work, see equation (1.20), no vortices exist in the recirculating flow in the steady-separation regime. Vortices first appear at the transition from the stationary to periodic regimes. As the vortices are the objects of interest in this work the main regime of interest is the periodic shedding regime where a vortex train has appeared behind the cylinder.

Another note is that the range of Reynolds numbers investigated in this work goes beyond the breakdown of the two-dimensionality of the flow. This means that the presented data will not mirror a true three-dimensional flow exactly. However, it is estimated that the maximal investigated Reynolds number is close enough to the two-dimensional breakdown for results to still be a reasonable approximation to the behavior of the actual three-dimensional flow.

(17)

1.3 Physics and Model Equations 1 THE PROBLEM AND UNDERLYING PHYSICS

1.3 Physics and Model Equations

The physics believed to govern all fluid flow on the length scale of interest in this project is classical Newtonian mechanics. That is, all length scales are assumed to be large enough to be able to consider the fluid as a continuum. This is in contrast to length scales where the individual molecules making up the fluid must be considered. A difference from classical mechanics is that for fluids one usually considers a velocity-field formulation instead of considering individual particles as is the usual approach in mechanics. This leads to a new form for the particle temporal derivative of a quantity, Q, for the fluid given (in 2D) as,

DQ Dt = ∂Q

∂t +u∂Q

∂x +v∂Q

∂y = ∂Q

∂t + (u· ∇)Q. (1.1) A full description and derivation of the difference between the particle and field approach may be found in [11, section 1-3].

1.3.1 General Navier-Stokes Equation

By applying Newtons second law of motion, ^P_iF_i =ma=m^Du_Dt, to a "box of fluid" at an arbitrary position in an inertial frame of reference and dividing by the volume of the box one obtains a general form of the Navier-Stokes (NS) equations (1.2), which govern the motion of the fluid [11, section 2-4]:

ρ ∂u

∂t +u· ∇u

=−∇p+∇ ·τ +f. (1.2)

Here, u is the flow velocity, ρ the fluid density, p the pressure, τ a stress tensor and f represents external body forces acting on the fluid.

Each term in the NS-equation may be interpreted physically as follows.

The term: ^∂u_∂t. This term quite simply describes the change in the velocity of the fluid caused directly by a change in time. This term is clearly non-zero for all instationary problems and zero for stationary problems.

The term: u· ∇u. This term is called the convective acceleration term or the non-linear advection term. This term represents the fact that the fluid may be accelerated due to its position in the domain independently of time. That is, even for stationary problems (^∂u_∂t = 0) a fluid may change velocity throughout the domain. As an example of this one may consider a steady flow of incompressible liquid from one reservoir to another through a diverging channel. See figure 1.3 for an illustration of the channel.

(18)

u u u u

Figure 1.3: Illustration of a diverging open ended channel in which a steady flow of fluid will be fastest in the narrow end and slowest in the wide end.

This phenomenon is relatively easily understood if one considers mass conservation and continuity. Only the fluid which has passed through the narrow part of the channel is available to fill the channel down stream where the channel is wider. As the fluid must remain continuous the only choice available is to slow down as it moves down the channel.

The terms: −∇p+∇ ·τ. Both terms represent effects of stresses in the fluid. ∇pis the gradient of the pressure and stems from the normal stresses in the fluid. ∇ ·τ stems from the anisotropic part of the stresses and describes viscous forces. This can be understood as the friction forces between layers of fluid pulling on each other as they pass one another.

In order to simplify the NS-equations for the solution of engineering problems a series of assumptions can be made for τ but going into depth with this is outside the scope of this project. A simple assumption is that the fluid is incompressible which reducesτ as described in the following.

Further explanations of the terms and assumptions which may be made can be found in [11, chapter 2].

The term: f. This term represents all exterior forcing on the fluid. Normally, this always includes gravity and may also include forcings on the system under consideration from e.g.

electrical fields.

1.3.2 Incompressible Navier-Stokes Equation

For the problems of interest in this project the fluid under consideration is assumed incompressible, temperature independent and it is assumed that the viscosity is constant throughout the fluid. These assumptions simplify the NS-equation to the form presented in equation 1.3,

ρ ∂u

∂t +u· ∇u

=−∇p+µ∇²u+f, (1.3)

It is easily seen that the change from the general NS-equation is the reduction of ∇ ·τ to µ∇²u which is denoted the viscous diffusion term. This term may be interpreted as modelling the diffusion of momentum which occurs for an incompressible viscous fluid.

(19)

The form of the NS-equation presented in (1.3) may be used directly as a component in solving an incompressible fluid problem. It is however beneficial to non-dimensionalize the NS-equation as this will reduce the parameters in the equation to a single parameter denoted the Reynolds number.

The first step is to introduce a mean flow velocityU and characteristic lengthL. From these a set of dimensionless functions and operators are then introduced as,

• Non-dimensional Velocity: uˇ = _U^u.

• Non-dimensional Pressure: pˇ= _ρU^p₂.

• Non-dimensional Body Force: ˇf= _ρU^L₂f.

• New Time Derivative: _∂^∂_ˇ_t = _U^L_∂t^∂.

• New Spatial Derivative: ∇ˇ =L∇.

By multiplying (1.3) which is an equation in terms of body forces, with the quantity _ρU^L₂, which has the unit of^h^m_N³ⁱ, one obtains the dimensionless equation,

L U²

∂u

∂t +u· ∇u

=− L

ρU²∇p+ µL

ρU²∇²u+ L

ρU²f,⇔ (1.4)

L U

∂

∂t u U + u

U · ∇u

U =−L∇ p

ρU² + µ

ρU LL²∇²u U + L

ρU²f. (1.5)

By substituting in the dimensionless functions and operators one obtains,

∂

∂ˇtuˇ+ ˇu· ∇ˇu=−∇ˇˇp+ µ ρU L

∇ˇ²uˇ+ ˇf. (1.6)

Introducing the Reynolds number as, Re= ^{ρU L}_µ and dropping the check marks to simplify the notation one arrives at the final non-dimensional NS-equation,

∂

∂tu+u· ∇u=−∇p+ 1

Re∇²u+f. (1.7)

Note thatL=Dfor the model problems considered here.

1.3.3 Mass Conservation

Another requirement for a fluid problem is that of mass conservation. This principle simply states that mass cannot appear out of nothing and likewise can not disappear into nothing.

Thus, we must require,

m=ρV = constant, (1.8)

(20)

where m is the mass and V is the volume. For an incompressible fluid any "box of fluid"

must have a constant density in space and time. From this consideration and by taking the particle derivative of equation 1.8 one gets the condition,

Dm Dt = Dρ

DtV +ρDV

Dt = 0⇒ ∇ ·u= 0. (1.9) See [11, section 2-3] for a full derivation of this equation, i.e. for the missing steps hidden in the⇒.

(1.9) is easily non-dimensionalised by introducing,

• Non-dimensional Velocity: uˇ = _U^u.

• New Spatial Derivative: ∇ˇ =L∇.

Multiplying (1.9) by ^L_U and rewriting to obtain,

L∇ · u

U = ˇ∇ ·uˇ = 0. (1.10)

Dropping the check mark we get the non-dimensional equivalent of (1.9), which happens to have the exact same form.

1.3.4 Model Equations

We have now considered all the components needed to formulate a model for the problem.

The two main model equations are thus the non-dimensional NS- and continuity-equations:

∂

∂tu+u· ∇u=−∇p+ 1

Re∇²u+f, (1.11)

∇ ·u= 0. (1.12)

Additionally, we need to formulate boundary and initial conditions for the problem. Con- sidering the problem geometry described in section 1.2 a set of boundary conditions may be formulated.

Cylinder Boundary Condition: The boundary constituted by the cylinder surface uses the usual FD no-slip condition. For the stationary cylinder this may be formulated as:

u|_(x,y)∈CS= (0,0), (1.13)

whereCS is the set of points which lie on the surface of the cylinder.

(21)

Wall Boundary: The boundary constituted by the moving wall also uses the usual FD no-slip condition, This may be formulated as:

u|_(x,y)∈WS= (u∞,0), (1.14)

whereWS is the set of points which lie on the wall.

Far Field Boundaries: As the method of choice for solving the problem is simulation, and as it is impossible to simulate an infinite domain, one needs to approximate the infinite half plane above the wall by a finite box. The concerns in the reader’s mind for which this approximation gives rise will be addressed in section 3. The "box" approximation gives rise to two different far field BC’s. The first is an inflow condition given by:

u|_(x,y)∈IF= (u∞,0), ∇p|_(x,y)∈IF= 0. (1.15) where IF is the set of points which lie on the inflow boundaries. These BCs represent a completely undisturbed flow as one would expect far away from the cylinder⁵.

The second is an outflow condition given by;

∇u|_(x,y)∈OF·n= (0,0), p|_(x,y)∈OF=p∞= 0. (1.16) whereOF is the set of points which lie on the outflow boundary,nis the outward pointing normal vector and p∞ is a reference pressure which is set equal to zero. This BC simulates that the fluid simply flows freely out of the domain at the outflow boundary.

An illustration of the problem domain with the different BC’s highlighted is provided in figure 1.4.

Cylinder Outflow Wall Inflow

Figure 1.4: Illustration of the location of the different boundary conditions for the problem. (CS): Cylinder surface (full line), (WS): Wall surface (densely dotted),(IF): Inflow (densely dashed),(OF): Outflow (loosely dashed).

5Due to numerical performance the pressure condition stated here is not used in the solver from theNek- tar++framework. Instead a higher order pressure condition derived to provide better numerical accuracy is used. Details of this condition and its derivation may be found in [18] and/or [16, section 8.3]

(22)

Initial Condition: The initial condition is chosen to be that of a free flow without the cylinder, i.e.

u|_t=0 = (u∞,0), p|_t=0=p∞= 0. (1.17) Since the investigation is not concerned with the initial transient solution, any reasonable choice of velocity and pressure fields as initial condition could be used as long as it is chosen to be compatible with the numerical method used.

1.4 Physical Quantities

For general FD problems many different quantities and characteristics of the flow may be of interest depending on the intended application. Additionally, a number of auxiliary quantities are important for analytical purposes in order to obtain the desired quantities and characteristics. This section outlines the quantities and characteristics of interest to this work.

The quantities of interest here are the pressure induced drag on the cylinder, the non- dimensionalized frequency of vortex shedding called the Strouhal number, the non-dimensionalized pressure at the cylinder surface called the pressure coefficient and the vortex movement pat- tern behind the cylinder. A few examples of areas where these quantities are important follow here.

In aerospace, car, boat and wind turbine engineering the drag on an object subjected to a fluid flow is of high interest as the minimization or maximization of the quantities are important for the performance of the object. The pressure coefficient at an object’s surface and Strouhal number may likewise be of interest in these areas of engineering. In fluid mixing, understanding the flow patterns has interest as this may determine the most efficient way to mix some reactant with a fluid to create a product. Understanding the flow patterns behind an object can also be of importance if other objects are to be placed in the wake of the first object.

An crucial quantity for the analysis of vortex patterns and structures which is of interest in this work is the vorticity as it may be used to uniquely define a vortex.

Vorticity: This quantity is defined in R³ as the curl of the velocity vector,

ω=∇ ×u, (1.18)

which in 2D translates to a scalar quantity, ω= ∂

∂xv− ∂

∂yu, (1.19)

since thex- andy-components of (1.18) are zero in 2D.

The vorticity at a given point in the fluid may be understood as a measure for how much of a rotating motion the fluid experiences near this point. A positive vorticity corresponds to

(23)

1.4 Physical Quantities 1 THE PROBLEM AND UNDERLYING PHYSICS

counter-clockwise rotation of the fluid while a negative vorticity corresponds to a clockwise rotation.

The introduction of the vorticity allows for a strict definition of a vortex which is used in this work to identify individual vortices uniquely.

A Vortex: In order to perform a rigorous investigation of the creation, annihilation and movement of vortices in a flow one must first define a vortex. In this work the center of a vortex is understood as an extremum in vorticity. This allows vortices to be defined uniquely by their centres from the extrema in the vorticity which are given by the points (xi, yi) in the flow where,

∇ω =0 ⇔ ω_x= 0 ∧ ω_y = 0, (1.20)

and where the Hessian of the vorticity is either positive or negative definite corresponding to minima and maxima, respectively.

From this definition it is possible to locate and count the number of vortices in a flow at a given time. With the capability to uniquely identify vortices it also becomes possible to trace the movement, creation and annihilation of vortices over time.

Strouhal Number: The Stouhal number St may be thought of as a dimensionless frequency for the shedding of vortices from an object immersed in fluid. It is given by,

St= L U∞

f, (1.21)

where f is the shedding frequency, L is a characteristic length scale in the problem under consideration and U∞ is the far field velocity of the fluid. In this workL=D.

Pressure Coefficient: The pressure coefficient,C_p may be thought of as a dimensionless measure for the pressure relative to the far field pressure and is defined by,

C_p= p−p∞ 1

2ρU_∞² , (1.22)

where p∞ is the far field pressure, U∞ is the far field velocity and ρ is the density of the fluid.

For the cylinder problem the base pressure coefficient C_bp is defined as C_p at the cylinder surface 180 degrees from the front of the cylinder.

(24)

Pressure Drag Coefficient: The drag force, Fd, on an object is defined as the sum of the forces acting on the object parallel to the direction of movement for the objects through a fluid. That is, it is the force which works to slow down the object as it moves through the fluid. In an incompressible viscous fluid the drag is composed of a force due to the pressure, Fp, of the fluid on the object and a shear force,Fs, from the layers of fluid sliding over the surface of the object. Choosing the direction of positive drag in the direction opposite the movement of the object we get,

F_d=Fp_d+Fs_d =− Z

pnˆ·ˆidS+ Z

τ ˆt·ˆi dS

, (1.23)

whereτ denotes the stress tensor component along the surface where ˆi is the unit vector in the direction of movement of the object. ˆt and ˆn are the normal and tangential vectors for the cylinder surface. In this work only the pressure-drag force is considered. From F_p_d a non-dimensional drag coefficient C_D_p is defined as,

CDp = Fp_d 1

2ρU_∞² . (1.24)

(25)

2 THEORY

2 Theory

This chapter provides an overview of the most important theoretical concepts from the different areas of mathematics used throughout the thesis. Section 2.1 outlines results and concepts used from Dynamical Systems Theory (DST) and presents an algorithm used to identify critical points for the vorticity. Section 2.2 states definitions and results for ortho- gonal polynomials and approximation theory which provides the basis for UQ and SEM.

Section 2.3 provides a brief overview of the strengths and weaknesses of the Spectral Ele- ment Method (SEM) used in the numerical simulations compared to the classical Finite Volume/Element Metohds (FVM/FEM). Here, an illustration of the main steps in the application of the SEM to a model problem is also provided to give the reader a basic understanding of how the model problems are solved numerically. An illustration of convergence results obtained using the method is also provided. Lastly, section 2.4 presents used concepts and results from Uncertainty Quantification (UQ).

2.1 Dynamical Systems Theory

The overarching aim of DST is to understand the structure of a dynamical system and how it changes depending on the parameters of the system. A two-dimensional dynamical system may be written generally as a system of first order ODE’s as,

x˙ y˙

!

= P(x, y, c) Q(x, y, c)

!

, (2.1)

where ˙xand ˙ymay be interpreted as a velocity given by the vector field (P(x, y, c), Q(x, y, c)) and c∈R^k is a k-dimensional vector of parameters. For systems of this form fixpoints, also known as critical points, are of special interest since they mark points where the behaviour of the system may change with small variations in the parameters c.

Critical points are defined as points in space where the right hand side of (2.1)is zero. Hence points where the velocity in the system is zero.

In [2] Brøns presents theorems which state that regular points, i.e. points which are not critical points are structurally stable, meaning that small changes in the parametersccannot change the structure of the dynamical system. In order to determine the type of a critical point the eigenvalues of the Jacobian of the governing dynamical system is considered⁶. Consider now a 2D flow problem where u(x, t) is a known velocity field. Consider the time tas a parameter and look at the system at t=t₀. This gives the dynamical system,

x˙ =u(x, t₀). (2.2)

The solutions to this system are called streamlines and describe the integral curves of the velocity field at the instant in timet0. That is, curves which a massless tracer particle would follow through the flow at a given instance in time. The critical points of this dynamical

6A well developed theory for identifying the critical points, their type, and the ways they may be created and/or annihilated when the parameterscare varied, exists, see e.g. [1].

(26)

system are in fluid dynamic terms denoted stagnation points, i.e. points where the fluid is at rest. For incompressible flow problems we may define a stream function Ψ, which fulfil u= ^∂Ψ_∂y, v=−^∂Ψ_∂x, [2]. The system (2.2) now takes the form,

x˙ = ∂x

∂t,∂y

∂t

= ∂Ψ

∂y,−∂Ψ

∂x

. (2.3)

The definition of Ψ means that it is constant along streamlines for (2.2) as may easily be proven,

dΨ dt = ∂Ψ

∂x dx

dt +∂Ψ

∂x dy dt = ∂Ψ

∂x

∂Ψ

∂y −∂Ψ

∂y

∂Ψ

∂x = 0, (2.4)

that is, the streamlines for (2.2) are level curves for Ψ. This fact means that critical points for (2.2) are stationary points for Ψ since these are simply special cases of level curves.

In [2] it is also stated that for any regular critical point of (2.3), i.e. points where the Jacobian of the system is regular, are structurally stable. This means that in order to locate structural changes in the flow the degenerate critical points where the Jacobian is singular must be considered.

Analysis using DST: In [2], Brøns derives a number of theoretical results for the dynamical system (2.3) describing the streamlines of the two-dimensional velocity field (u,v) of incompressible fluid flow using the stream function Ψ. As will be shown below an equivalent system may be formulated for the vorticityω. This means that all the results developed by Brøns may be used to analyse the system governed byω instead of Ψ. That is, the theory may be used to analyse the creation, annihilation and movement patterns of vortices in a 2D flow problem.

As stated earlier a central goal of the present work is to implement and demonstrate a method for rigorously investigating the creation, annihilation and dynamics of vortices⁷. The first step on this path is to obtain the velocity field for the fluid and from this calculate the vorticity. This is done by solving the governing equations numerically, which will be discussed in detail later. In the following the velocity field and thus the vorticity is assumed known.

It turns out that the dynamics of the stationary points for vorticity and thus of vortices may be understood using DST. This is done by considering time as a parameter instead of a continuous variable and investigating the critical points of a family of dynamical systems of the form,

∂x

∂s∂y

∂s

!

= x˙ y˙

!

= ω_y(x, y, t_i, c)

−ω_x(x, y, ti, c)

!

, i∈ {1, ..., N}, c∈R^k, k∈N. (2.5) Here (x, y) are spatial coordinates,ω the vorticity andc∈R^k are a vector of parameters. s is an artificial time introduced to consider how the flow would move when the real time is

7Remember that a vortex has been defined from a vorticity extremum, i.e. points in space and time where ωx(x, y, t) =ωy(x, y, t) = 0 and where the Hessian is positive or negative definite, see section 1.4.

(27)

2.1 Dynamical Systems Theory 2 THEORY

frozen. Notice that we have frozen time and that the N separate time steps are considered as individual dynamical systems. Consider now the vorticity field at one of these snapshots in time, ω(x, y, t_i, c). If this quantity is differentiated with respect to s, utilizing the chain rule one gets,

d

dsω(x, y, t_i) = ∂ω

∂x

∂s +∂ω

∂y

∂s =ω_xω_y−ω_yω_x = 0, (2.6) where we have used from (2.5) that (^∂x_∂s,^∂y_∂s) = (ω_y,−ω_x). This is in complete analogy to (2.4).

Following the analogy, it is seen from (2.6) that the streamlines given by the solutions to (2.5) are level curves forω(x, y, t_i, c) with respect to sand thus the extrema ofω(x, y, t_i, c) are special cases of level curves which haveωx =ωy = 0. From this it follows that a critical point of the dynamical system (2.5) corresponds to stationary points for ω. These may either be extrema, saddle points or degenerate points. As stated above a degenerate point is a point where the Jacobian matrix for the system (2.5) is singular and it is only at degenerate critical points that structural changes in (2.5) may occur, see [2, Theorem 3].

This means that if one can identify all stationary points of the vorticity one finds all the critical points of the system (2.5) and thus also all vortices present in the flow. The critical points and any possible structural changes they may undergo, corresponding to the creation or annihilation of vortices, can then be analysed using the machinery of DST. This will be elaborated on below as the saddle-center bifurcation is treated.

Comment on critical points for (2.5)and extrema/saddles for vorticity: By com- paring the Hessian matrix for the vorticity and the Jacobian matrix for the system (2.5) it is immediately observed that extrema for the vorticity corresponds to centres for (2.5) and saddles for the vorticity corresponds to saddle points for (2.5). Hence identifying a center for (2.5) corresponds to identifying a vortex in the flow.

Nullclines: Consider (2.1). A nullcline is defined as a curve along which either ˙x or ˙y is zero. Thus a nullcline corresponds to the zero contours of either ˙x or ˙y. In 2D one may use nullclines to identify critical points by looking for intercepts between the contours ˙x= 0 and y˙= 0.

The technique of using the intercepts between nullclines is used for identifications of critical points for ω. This means that by identifying ω_x = 0 one then only has to solve a one dimensional problem along the curves ω_x = 0 to identify ω_y|_(ω_x₌₀₎ = 0 and thereby the critical points.

Parameter dependence: Any spatial physical instationary problem has by definition some dependence on space and time. Besides this dependence the problem most likely also depends on a range of parameters. These dependencies may be on boundary and initial conditions, the domain shape and for a fluid problem the Reynolds number. The task now becomes to identify at what parameter values the structure of the system (2.1) changes and what changes occur. These changes in structure as a function of parameter variations are denoted bifurcations.

(28)

2.1.1 Bifurcations

Meiss [1] offers a definition of a bifurcation as,

Bifurcation: a qualitative change in dynamics occurring upon a small change in a parameter.

Many different types of changes in dynamics can occur in (2.1) as a consequence of varying system parameters, see [1, Chapter 8] for a treatment of some of these. Some bifurcations occur as a single parameter is varied while others occur only as a consequence of varying multiple parameters at once. The minimum number of parameters which must be varied for a given bifurcation to occur is denoted theCodimension of the bifurcation.

For the systems investigated in the present work the codimension 1 saddle-center bifurcation was the only bifurcation observed away from domain boundaries in all but a single case⁸. In the single exceptional case a co-dimension 2 bifurcation was identified. Due to its frequent occurrence the saddle-center bifurcation is covered in detail here, while the co-dimension 2 bifurcation is treated specifically in chapter 8 when it is encountered.

Saddle-Center Bifurcation: This bifurcation corresponds to the creation or annihilation of two critical points as a single parameter is varied. One of the critical points is a saddle point and the other is a center.

In order to illustrate such a bifurcation consider a Taylor expansion of the vorticity given by,

ω=

∞

X

n,m=0

anmxⁿy^m, (2.7)

In [2] Brøns shows that if one applies transformations to the dynamical system in consideration such that the degenerate critical point undergoing a bifurcation is centred at the origin and the coordinate axis are twisted and stretched appropriately the following theorem holds⁹,

Theorem 1: Let a₁₀, a₀₁, a₁₁, a₂₀ of (2.7)and ˜a_n0 for n < N−1 be small parameters, and assume that non-degeneracy conditions a₀₂ 6= 0 and ˜a_N0 6= 0. Then there is a coordinate transformation that brings ω into the normal form,

ω= σ

2y²+f(x) +O(N+ 1), (2.8)

f(x) =

N

X

n=1

cnxⁿ, cN−1 = 0, cN = 1

N, (2.9)

σ =

( −1 for N even and a02/˜aN0 <0, 1 for N even and a₀₂/˜a_N0 >0 or N odd,

)

(2.10) and c₁, ..., c_N−2 are small parameters.

Above the ˜a_n0’s are given by,

8In [2, part III] Brøns shows that for a flow away from a boundary this is the only possible codimension 1 bifurcation.

9This isTheorem 5in [2] formulated forωinstead of Ψ.

(29)

2.1 Dynamical Systems Theory 2 THEORY

a˜_N0 =a_N0+ a non-linear combination of a_nm of lower order then N. (2.11) Using this transformed expression for ω it is easy to see that all critical points near the degenerate point for the system (2.5) are situated along the transformed x-axis and that there are at mostN −1 of them,

x˙ = ∂ω

∂y =σy= 0⇔y= 0, y˙=−∂ω

∂x =−f_x(x) = 0. (2.12) The Jacobian of (2.5) at a critical point, determines its type and is given by,

J_{( ˙}_x,_y)_˙ =

"

0 σ

−f_xx(x) 0

#

. (2.13)

From (2.13) one can determine the type of a critical point by considering the eigenvalues given by,λ=±^p−σf_xx(x). This leaves three options. Either the critical point is degenerate allowing for a bifurcation to happen. The second option is −σf_xx(x) >0 which given two real eigenvalues of opposite sign, corresponding to a saddle as the critical point has one attracting eigendirection and one repelling eigendirection. The last option is−σf_xx(x)<0 which gives two purely imaginary eigenvalues of opposite sign which corresponds to a center, i.e. either a maximum or minimum in the vorticity.

For a co-dimension 1 bifurcation one has N = 3, i.e. the terms of order lower then four determine the behaviour of the system. From Theorem 1 this gives σ = 1 and reduces the expression for ω to,

ω= 1

2y²+f(x) +O(4), (2.14)

f(x) =c₁x+ 1 3x³.

Using (2.12) and (2.13) it can be seen that the occurrence of a bifurcation as c₁ is varied may be illustrated in one dimension aroundx= 0 by considering the derivative off(x),

f_x(x, c₁) =c₁+x², (2.15)

From (2.15) it is seen that the origin is a degenerate critical point for the system (2.5) when c1 = 0. It is easy to see that in general critical points exist for (2.5) at,

x=±√

−c₁ ⇔f_x(x, c₁) = 0. (2.16)

For c₁ = 0 we see that f_x(0,0) = 0 and _∂x^∂ f_x(0,0) = f_xx(0,0) = 0 while _∂x^∂²₂f_x(0,0) = f_xxx(0,0) = 2 and _∂c^∂

1f_x(0,0) = 1. The first two conditions are known as singularity

Master Thesis Advanced Techniques for Investigating Structures in Computational Fluid Dynamics

Master Thesis

Author:

Abstrakt

Contents

Introduction

Notation

1 The Problem and Underlying Physics

2 Theory