**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

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).

**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.

CONTENTS CONTENTS

**Contents**

**Introduction** **iv**

**Notation** **v**

**1** **The Problem and Underlying Physics** **1**

1.1 The Purpose . . . 1

1.2 The Model Problems . . . 2

1.3 Physics and Model Equations . . . 6

1.3.1 General Navier-Stokes Equation . . . 6

1.3.2 Incompressible Navier-Stokes Equation . . . 7

1.3.3 Mass Conservation . . . 8

1.3.4 Model Equations . . . 9

1.4 Physical Quantities . . . 11

A Vortex: . . . 12

**2** **Theory** **14**
2.1 Dynamical Systems Theory . . . 14

2.1.1 Bifurcations . . . 17

2.1.2 Tools for Analysing Flow Topology . . . 20

2.2 Orthogonal Polynomials and Approximation Theory . . . 21

2.3 The Spectral Element Method . . . 24

2.3.1 Overview of the Method . . . 26

2.4 Uncertainty Quantification . . . 31

2.4.1 Statistical Quantities . . . 32

2.4.2 The Monte Carlo Approach . . . 32

2.4.3 generalized Polynomial Chaos . . . 33

2.4.4 Stochastic Collocation Approach . . . 35

**3** **Discretizing the Problem** **38**
3.1 The Time Stepping Scheme . . . 38

3.2 The Discrete Domain . . . 40

**4** **Software** **44**

4.1 Pre Processing . . . 44

4.2 Simulation . . . 45

4.3 Post Processing . . . 46

4.3.1 The Finished Simulation Package . . . 49

4.3.2 UQ Software . . . 50

4.4 Hardware . . . 51

**5** **Parallel Execution with MPI** **52**
**6** **Simulations** **58**
6.1 Visualization . . . 59

6.2 Critical Point Identification . . . 60

6.3 Validation Simulations . . . 63

6.4 Cylinder near the Moving Wall . . . 63

6.5 Uncertainty Quantification . . . 64

**7** **Validation** **66**
7.1 Convergence . . . 66

7.2 Comparing Results . . . 68

7.2.1 Free Flow . . . 69

7.2.2 Cylinder Near Wall . . . 72

7.3 Death of the Transient Solution . . . 73

**8** **Analysis** **76**
8.1 Initial Investigation for the Cylinder and Wall . . . 76

8.1.1 Stabilizing Effect of the Wall . . . 76

8.1.2 Different Critical Point Structures . . . 80

8.1.3 Downstream Surviving Vortices, Creation Points and Movement Pat- terns . . . 82

8.1.4 Vortex Strength Reduction . . . 84

8.2 Formation and Disappearance of Extrema-Saddle Pair . . . 93

8.2.1 Constant *Re* Bifurcation Diagrams . . . 94

8.3 Vortex Creation Point Jumping Downstream . . . 101

8.4 Uncertainty Quantification . . . 107

8.4.1 Cylinder in free flow . . . 107

8.4.2 Cylinder Near Moving Wall . . . 114

8.5 Method Limitation . . . 120

CONTENTS CONTENTS

**9** **Conclusion and Future Work** **123**

**10 References** **126**

**A Appendix** **128**

A.1 Sample mesh used for simulations . . . 128

A.2 Nektar++ XML Extract . . . 131

A.3 The Nektar++ framework . . . 132

A.3.1 Installing Nektar++ . . . 133

A.3.2 Navigating Nektar++ . . . 133

A.3.3 Setting Up a Problem . . . 136

A.4 Validation of Nektar++ . . . 136

A.4.1 Simple Domain Convergence Test . . . 137

A.4.2 Incompressible Navier-Stokes Solver . . . 140

A.5 Authors Software . . . 144

A.5.1 Shell Scripts . . . 144

A.5.2 Python Scripts . . . 150

A.5.3 C++ Program . . . 157

A.5.4 FEniCS Poisson Solver . . . 161

A.5.5 MATLAB Implementation of the SCM . . . 162

**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 presen-
ted^{1}. 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 tool**Paraview** [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 the**Nektar++**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**.geo**and**.msh**file containing an example
of the mesh used in the simulations and an example of XML-files to use as input files for
the**Nektar++**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.

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 denoted***p.*

• **Mass: A mass is denoted***m.*

• **Volume: A volume is denoted***V*.

• **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}_{2} 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 denoted***P*.

• **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.**

• **Drag force: Drag on an object is denoted***F**d*.

• **Pressure Drag Coefficient: This coefficient is denoted** *C*_{D}* _{p}* =

^{p−p}_{1}

^{∞}

2*ρU*∞^{2}

.

• **Mean: The statistical mean of a quantity,** *X* is denoted, *µ** _{X}*.

• **Variance: The variance of a quantity,** *X* is denoted, *σ*_{X}^{2}.

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

• **Gaussian Distribution:** A normal distribution with mean*a*and standard deviation
*b*is 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**
*Q*is denoted: ^{DQ}* _{Dt}*.

**– Partial Derivative: Partial differentiation is denoted in one of the following**
ways: _{∂b}^{∂}*a,* ^{∂a}* _{∂b}* or

*a*

*b*.

**– Absolute Derivative: The absolute derivative is denoted** _{db}^{d}*a*or ^{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*_{2}**-norm: The standard***L*_{2}-norm is denoted k · k_{L}_{2} where ·is a place holder for the
element to be measured.

• *H*_{w}^{k}**-norm: The weighted Sobolev-norm is denoted** k · k_{H}*k*

*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 parameters*P* is considered on the subset Ω of the full parameter domain
this is denoted as*f*(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

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.

**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^{2}
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*^{−}^{1}^{2}) 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 pro-
ject, 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^{3}, 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^{4}.

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.

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 under- standing 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,

*D*is 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*∞

*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

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. *D*is the cylinder diameter and*G*the size of the gap between the cylinder and wall. *D*
and*G*are 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*
and*G*are the only geometric parameters in the system. Thus scaling*D*by a factor*a*simply
corresponds to shrinking*G*by the same factor*a. This means that the parameter of interest*
in the following is their ratio, chosen as ^{D}* _{G}*. The choice of

^{D}*opposed to*

_{G}

_{D}*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*

^{G}

^{D}*= 0 instead of*

_{G}

_{D}*→ ∞.*

^{G}From the present discussion three independent parameters for the problem have been iden-
tified. 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.

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 recircu-
lation 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 > Re*crit. 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}_{2} ≈190 to be the Reynolds number at which the
two-dimensionality of the flow breaks down. Beyond *Re*crit_{2} the flow exhibits three-
dimensional behaviour and the 2D simulations are no longer enough to capture the
physics completely.

(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.

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}_{i}*F** _{i}* =

*ma*=

*m*

^{Du}*, 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]:*

_{Dt}*ρ*
*∂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.

**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 in- compressible, 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+*µ∇*^{2}**u**+**f,** (1.3)

It is easily seen that the change from the general NS-equation is the reduction of ∇ ·* τ*
to

*µ∇*

^{2}

**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.

1.3 Physics and Model Equations 1 THE PROBLEM AND UNDERLYING PHYSICS

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 velocity*U* and characteristic length*L. 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}_{2}.

• **Non-dimensional Body Force:** ˇ**f**= _{ρU}^{L}_{2}**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}_{2},
which has the unit of^{h}^{m}_{N}^{3}^{i}, one obtains the dimensionless equation,

*L*
*U*^{2}

*∂u*

*∂t* +**u**· ∇u

=− *L*

*ρU*^{2}∇p+ *µL*

*ρU*^{2}∇^{2}**u**+ *L*

*ρU*^{2}**f,**⇔ (1.4)

*L*
*U*

*∂*

*∂t*
**u**
*U* + **u**

*U* · ∇**u**

*U* =−L∇ *p*

*ρU*^{2} + *µ*

*ρU LL*^{2}∇^{2}**u**
*U* + *L*

*ρU*^{2}**f.** (1.5)

By substituting in the dimensionless functions and operators one obtains,

*∂*

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

∇ˇ^{2}**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,

*∂*

*∂t***u**+**u**· ∇u=−∇p+ 1

*Re*∇^{2}**u**+**f.** (1.7)

Note that*L*=*D*for 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)

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:

*∂*

*∂t***u**+**u**· ∇u=−∇p+ 1

*Re*∇^{2}**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)

where**CS** is the set of points which lie on the surface of the cylinder.

1.3 Physics and Model Equations 1 THE PROBLEM AND UNDERLYING PHYSICS

**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)

where**WS** 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^{5}.

The second is an outflow condition given by;

∇u|_{(x,y)∈OF}·**n**= (0,0), *p|*_{(x,y)∈OF}=*p*∞= 0. (1.16)
where**OF** is the set of points which lie on the outflow boundary,**n**is 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 the**Nek-**
**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]

**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^{3} as the curl of the velocity vector,

ω=∇ ×**u,** (1.18)

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

*∂xv*− *∂*

*∂yu,* (1.19)

since the*x- 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

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 (x*i**, y**i*) in
the flow where,

∇ω =**0** ⇔ *ω** _{x}*= 0 ∧

*ω*

*= 0, (1.20)*

_{y}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 fre-
quency 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 work*L*=*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*_{∞}^{2} *,* (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*

*at the cylinder surface 180 degrees from the front of the cylinder.*

_{p}**Pressure Drag Coefficient:** The drag force, *F**d*, 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,
*F**p*, of the fluid on the object and a shear force,*F**s*, 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}*=

*F*

*p*

*+*

_{d}*F*

*s*

*=− Z*

_{d}*p***n**ˆ·ˆ**i***dS*+
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}*is defined as,*

_{p}*C**D**p* = *F**p** _{d}*
1

2*ρU*_{∞}^{2} *.* (1.24)

2 THEORY

**2** **Theory**

This chapter provides an overview of the most important theoretical concepts from the dif- ferent 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 ap- plication of the SEM to a model problem is also provided to give the reader a basic under- standing 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 ˙*x*and ˙*y*may 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 parameters*c*cannot
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^{6}.
Consider now a 2D flow problem where **u(x, t) is a known velocity field. Consider the time**
*t*as a parameter and look at the system at *t*=*t*_{0}. This gives the dynamical system,

**x**˙ =**u(x, t**_{0}). (2.2)

The solutions to this system are called streamlines and describe the integral curves of the
velocity field at the instant in time*t*0. 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 parameters*c*are varied, exists, see e.g. [1].

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 dy-
namical 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^{7}.
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, t

*i*

*, c)*

!

*, i*∈ {1, ..., N}, c∈R^{k}*, k*∈N*.* (2.5)
Here (x, y) are spatial coordinates,*ω* the vorticity and*c*∈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.

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*

*∂x*

*∂s* +*∂ω*

*∂y*

*∂y*

*∂s* =*ω*_{x}*ω** _{y}*−

*ω*

_{y}*ω*

*= 0, (2.6) where we have used from (2.5) that (*

_{x}

^{∂x}

_{∂s}*,*

^{∂y}*) = (ω*

_{∂s}

_{y}*,*−ω

*). This is in complete analogy to (2.4).*

_{x}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* *s*and 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

*ω*

*= 0 to identify*

_{x}*ω*

*|*

_{y}_{(ω}

_{x}_{=0)}= 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.

**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 the*Codimension* 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^{8}. 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*

*a**nm**x*^{n}*y*^{m}*,* (2.7)

In [2] Brøns shows that if one applies transformations to the dynamical system in con-
sideration 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^{9},

**Theorem 1:** *Let* *a*_{10}*, a*_{01}*, a*_{11}*, a*_{20} *of* (2.7)*and* ˜*a*_{n0}*for* *n < N*−1 *be small parameters, and*
*assume that non-degeneracy conditions* *a*_{02} 6= 0 *and* ˜*a** _{N0}* 6= 0. Then there is a coordinate

*transformation that brings*

*ω*

*into the normal form,*

*ω*= *σ*

2*y*^{2}+*f*(x) +*O(N*+ 1), (2.8)

*f(x) =*

*N*

X

*n=1*

*c**n**x*^{n}*, c**N−1* = 0, c*N* = 1

*N,* (2.9)

*σ* =

( −1 *for N even and* *a*02*/˜a**N0* *<*0,
1 *for N even and* *a*_{02}*/˜a*_{N0}*>*0 *or N odd,*

)

(2.10)
*and* *c*_{1}*, ..., 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 is**Theorem 5**in [2] formulated for*ω*instead of Ψ.

2.1 Dynamical Systems Theory 2 THEORY

*a*˜* _{N0}* =

*a*

*+ a non-linear combination of*

_{N0}*a*

*of lower order then*

_{nm}*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 most

*N*−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

*(x)*

_{xx}*>*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

*(x)*

_{xx}*<*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

2*y*^{2}+*f*(x) +*O(4),* (2.14)

*f(x) =c*_{1}*x*+ 1
3*x*^{3}*.*

Using (2.12) and (2.13) it can be seen that the occurrence of a bifurcation as *c*_{1} is varied
may be illustrated in one dimension around*x*= 0 by considering the derivative of*f*(x),

*f** _{x}*(x, c

_{1}) =

*c*

_{1}+

*x*

^{2}

*,*(2.15)

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

*x*=±√

−c_{1} ⇔*f** _{x}*(x, c

_{1}) = 0. (2.16)

For *c*_{1} = 0 we see that *f** _{x}*(0,0) = 0 and

_{∂x}

^{∂}*f*

*(0,0) =*

_{x}*f*

*(0,0) = 0 while*

_{xx}

_{∂x}

^{∂}^{2}

_{2}

*f*

*(0,0) =*

_{x}*f*

*(0,0) = 2 and*

_{xxx}

_{∂c}

^{∂}1*f** _{x}*(0,0) = 1. The first two conditions are known as singularity