Contour dynamics for vortex-induced advection

(1)

Contour dynamics for vortex-induced advection

Citation for published version (APA):

Schoemaker, R. M. (2003). Contour dynamics for vortex-induced advection. Technische Universiteit Eindhoven. https://doi.org/10.6100/IR564542

DOI:

10.6100/IR564542

Document status and date: Published: 01/01/2003 Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

providing details and we will investigate your claim.

(2)

Contour Dynamics for

Vortex-induced Advection

(3)

Druk : Universiteitsdrukkerij TU/e

Bij de omslag : Evolutie van een gemodelleerde polaire vortex in de planetaire niet-uniforme achtergrondrotatie.

Omslag ontwerp : JWL Producties

CIP-DATA LIBRARY TECHNISCHE UNIVERSITEIT EINDHOVEN Schoemaker, Robin Micha¨el

Contour Dynamics for Vortex-induced Advection / by Robin Schoemaker. – Eindhoven : Technische Universiteit Eindhoven, 2003. – Proefschrift.

ISBN 90-386-1575-2 NUR 910

Trefw. : geofysische stromingsleer / werveldynamica / deeltjesadvectie / numerieke simulaties / parallellisatie

Subject headings : geophysical ﬂuid dynamics / vortex dynamics / particle advection / numerical simulations / parallelization

(4)

Contour Dynamics for

Vortex-induced Advection

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Eindhoven, op gezag van

de Rector Magniﬁcus, prof.dr. R.A. van Santen, voor een commissie aangewezen door het College voor Promoties in het openbaar te verdedigen op

maandag 26 mei 2003 om 16.00 uur

door

Robin Micha¨el Schoemaker geboren te Almelo

(5)

prof.dr.ir. G.J.F. van Heijst en

prof.dr. R.M.M. Mattheij

Copromotor : dr. H.J.H. Clercx

(6)

4.2.3 Critical and asymptotic values for a . . . 44 4.3 Contour perturbations . . . 47 4.3.1 T he tripolar vortex . . . 48 4.3.2 T he triangular vortex . . . 53 4.3.3 T he square vortex . . . 55 4.4 Multiple-patch equilibria . . . 58 4.4.1 T he tripolar vortex . . . 60 4.4.2 T he triangular vortex . . . 68 4.4.3 T he square vortex . . . 72 4.5 Discussion . . . 77

5 A parallel hierarchical-element method 79 5.1 Introduction . . . 80

5.2 Fast multipole and hierarchical-element methods . . . 81

5.3 A hierarchical-element method for contour dynamics . . . 83

5.4 Algorithmic complexity . . . 85

5.5 Parallelization strategy . . . 88

5.5.1 OpenMP . . . 89

5.5.2 Performance issues . . . 90

5.5.3 Parallel domains . . . 92

5.6 Numerical experiments and discussion . . . 95

5.6.1 Example 1 . . . 95

5.6.2 Example 2 . . . 97

5.6.3 Example 3 . . . 98

5.7 Concluding remarks . . . 101

6 Advection in idealized polar vortex ﬂows 103 6.1 Introduction . . . 104

6.2 Quantifying transport . . . 106

6.2.1 Contour lengthening . . . 106

6.2.2 Patch erosion . . . 106

6.3 Polar vortex modelling . . . 107

6.4 Numerical experiments and results . . . 109

6.4.1 Planar and spherical polar ﬂow simulations . . . 109

6.4.2 Advection at the vortex edge . . . 114

6.5 Discussion . . . 119

7 Closure and outlook 123 7.1 Concluding remarks . . . 123

(8)

Contents iii

A OpenMP implementation 127

A.1 Boxwise HEM . . . 127 A.2 Pseudo-code example . . . 128

B The modelled polar vortex 131

Bibliography 135

(9)

(10)

1 Introduction

The (coupled) atmospheric and oceanic fluid systems of the Earth show mo-tions on scales ranging from meters to the circumference of the planet itself. The dynamics on the smaller scales is the realm of the three-dimensional Navier-Stokes equations, which are nonlinear second-order equations for fluids with viscosity. These equations are hard to solve analytically as well as numerically. Investigations of large-scale phenomena for motions on planetary scales, how-ever, are efficiently described with simplified mathematical models. This thesis deals with a specific numerical method appropriate for a special kind of fluid flow that is modelled as two-dimensional, incompressible, and inviscid. The equations at hand are the so-called 2D Euler equations. This introduction will pursue the point why this model is relevant for planetary flows and general two-dimensional vortex flows, why the numerical technique and its advancements is adequate, and moreover, how we can apply our numerical method to investigate contemporary problems in transport of passive material in typical vortex flows.

1.1 Large-scale vortex flows

The property of inviscidness categorizes a flow as ideal and implies that mo-mentum dissipation is absent and only first-order derivatives in the momo-mentum equations are needed to solve the problem at hand. Although most real fluid flows have a finite viscosity, an accurate description of the flow kinematics and dynamics is often conceivable when viscous effects are neglected. The effect of viscosity in a typical flow setting can be measured by means of a parameter, the Reynolds number Re, which indicates the role of viscous forces compared to the inertial forces via Re = ρ U L/µ, with ρ the fluid density, µ the dynamic viscosity, U a characteristic velocity, and L a characteristic length scale in the flow setting. Viscous effects play a crucial role in creeping flows like melting glass or lava flows, for example, or flows in microfluidic devices. The Reynolds number is then O(1) or smaller. Most of us, however, encounter fluid flows

(11)

characterized by much higher Reynolds numbers. Examples are the stirring of coﬀee (Re∼ 104) with the cup diameter as a characteristic length scale, cruis-ing on the freeway (Re∼ 107_{) for the ﬂow passing your car, while for an airfoil}

the Reynolds number is of the order Re∼ 108_{. The large-scale coherent ﬂow}

regions in the atmosphere and oceans are typified by Reynolds numbers of the order Re ∼ 1010 for oceanic eddies and Re ∼ 1012 for atmospheric low/high-pressure systems. In these high-Reynolds number flows viscous effects play a role at very small length scales only—where energy is dissipated rapidly—or at very long convective time scales. These systems do not feel the effect of viscosity in case of very large-scale motion.

Large-scale vortices are coherent structures that can be found in the oceans, the atmosphere of our planet as well as in the atmosphere of other planets. Examples of terrestrial nature are high and low-pressure areas, eddies in the ocean and the large polar vortex (Arctic as well as Antarctic). Other planetary examples are the Great Red Spot on Jupiter (Figure 1.1), the Great Dark Spot on Neptune and huge rotating dust structures on Mars (Figure 1.2). The

Figure 1.1 — The Great Red Spot of Jupiter (largest oval). It exists in this form at least since the 17thcentury. Public domain image from NASA/GSFC.

thickness of the layer of fluid these structures evolve in (on Earth: 1 – 10 km), is small compared to the horizontal size of the coherent structure itself (on Earth: 100 – 5000 km). This geometrical confinement, together with the even more important planetary rotation and the density stratification in the fluid layer, implies quasi two-dimensionality of the flow. The dynamics of these planetary flows are described to good approximation by the two-dimensional, incompressible, inviscid variant of the Navier-Stokes equation: the 2D Euler equation.

While huge monopolar vortices dwell in the oceans and in the atmosphere on a daily basis, other stable coherent vortex structures with multipolar char-acteristics have been observed in nature as well, be it under rare circumstances. The dipolar system, consisting of two monopoles of opposite polarity, having

(12)

1.1 | Large-scale vortex ﬂows 3

Figure 1.2 — A Martian polar dust storm of approximately 900 km in diameter with jet-like and vortical (dipolar) behaviour. It is moving as a front, away from the Martian north pole, which can be seen partly on the right-half of the picture. Public domain image from NASA/JPL/MGS.

a linear momentum or combination of linear momentum and angular momen-tum, as well as the tripolar vortex, are infrequently observed phenomena. The tripolar system consisting of a core region of one polarity and two satellite re-gions of opposite polarity and steadily rotating according to the sign of the core vorticity, was actually observed for the first time in a laboratory experiment where it evolved from an initially perturbed shielded monopolar vortex in a rotating fluid (van Heijst & Kloosterziel [37]). The most famous example of an observed tripole in nature is reported, via satellite imagery, by Pingree & Le Cann [55] who observed a tripole in the Bay of Biscay by an infra-red sea-surface picture (observed in January 1990). The important observation is that such multipolar vortices only exist through the two-dimensionality of the flow. Higher order multipolar vortices have not been observed in nature, yet have been investigated theoretically, numerically, and in the laboratory (see Ref. [11] and Figure 1.3).

The largest monopolar vortices found on our planet are the cyclonic vor-tices residing at each pole of our planet—the Arctic and Antarctic polar vortex. Polar vortices are not limited to the Earth. The planet Mars has a strong polar jet [2] and Jupiter has a hexagonal shaped polar vortex likewise characterized by intense jetstreams at its edge [1]. Although a polar vortex has in real-ity a three-dimensional structure, the two-dimensionalreal-ity as introduced above still applies to such a vortex and to its layer-like composition. Approximately one-third of the transport in polar vortex ﬂows is achieved through horizon-tal mixing mechanisms—averaged over long time scales. A polar vortex covers a region that is marked by a rather uniform distribution of positive relative vorticity (cyclonic vortex). An important observation is that the edge of the distribution declines rather steeply to zero or negative values. The thus cre-ated gradient therefore induces large velocity ﬁelds which forms a barrier to transport of inward moving passive material. This so-called polar night jet

(13)

sig-natures the overall geometrical extension of the polar vortex. Perturbations of this polar jet induces large planetary waves (Rossby waves) resulting in erosion of the polar edge. These characteristics are quite familiar through the many studies that have been (and still are being) performed, and which often focus on the depletion of ozone at the polar regions. These studies comprise a wide spectrum of results, from observational results of satellite and aircraft data (e.g. Schoeberl & Hartmann [63]), to laboratory experiments (e.g. Sommeria et al. [68]), and numerical experiments (e.g. Koh & Legras [42], Norton [51], Pierce & Fairlie [53], Waugh [76], Waugh & Plumb [77]). The numerical models, in this respect, are in most cases approximations and have shown to be quite successful in a better understanding of the underlying dynamics. One impor-tant reason behind this is that high-resolution computer simulations are able to show ﬁne-scale processes that are often hard to resolve with observational analysis. Moreover, the large-scale dynamics is often captured for simulations with complete geophysical domains like a (hemi)sphere (see also Section 1.4).

1.2 Passive tracer advection

Passive tracers are constituents of the flow that do not influence the dynamics of the flow field and are as such transported by means of the prescribed velocity field present. An example are the individual fluid elements themselves which can be made visible by means of a tracing marker like coloured dye in labora-tory experiments, for example. Contrary to passive tracers are active tracers which can influence the dynamics of the system by contributing to the mixture density, pressure, salinity and/or heat content. For example, temperature or salinity perturb the density and change the velocity field accordingly. Vorticity is another example of an active tracer that influences the velocity field in an unequivocal way. Reactant tracers (methane, nitrogen, e.g.), chemical pollu-tants, volcanic ash, etc. are other examples that are able to influence the flow field. Active tracers are, however, often modelled as passive tracers (see Sec-tion 2.2). The large-scale transfer of fluid properties by a prescribed velocity field is termed advection, and if passive tracers are involved we call this mech-anism passive tracer advection. In this thesis, the advection of passive tracers is in particular studied in a region just outside the polar vortex, the so-called ‘surf zone’, where vortex air erodes from the polar vortex and forms thin fil-aments and secondary vorticities in the large mixing zones. Efficient mixing is indicative to chaotic flow properties, while inefficient mixing is exhibited by transport barriers. One of the goals of this research is the employment of ap-propriate quantitative techniques in order to verify the numerical methods and to gain a better understanding of the horizontal flow dynamics in the polar region.

(14)

1.3 | Multipolar equilibria 5 1.3 Multipolar equilibria

Many results can be found in the literature on experimental and numerical studies of multipolar vortices. See, for example, Beckers & van Heijst [13] (Fig-ure 1.3) and Morel & Carton [49]. In Chapter 4 we will explore a recently

Figure1.3 — A triangular vortex in a laboratory experiment (Beckers & van Heijst [13]).

discovered family of multipolar vortices that are analytical solutions to the 2D Euler equations (Crowdy [20]). These solutions are described by two prime pa-rameters, are multipolar in nature, and designate static vortical equilibria con-sisting of singular vortex extrema (point vortices) embedded in a finite area of uniform vorticity. The point vortices are distributed in a polygonal way around one centrally situated point vortex. One parameter indicates the multipolar aspect of the configuration, i.e. the number of point vortices involved, while the other parameter marks the important non-trivial shape of the enveloping uniform vorticity distribution. The characteristic shape of the bounding patch is such that the induced velocity field is zero outside the vorticity patch and zero at the location of the point vortices, implying a static non-rotating equi-librium. Moreover, by replacing the point vortices by finite areas of uniform vorticity, similar yet novel multiple-patch equilibria are obtained with different characteristics. The solutions found in Crowdy [20] and the novel multiple-patch equilibria are ideally suited for verification and investigation for stability and advection properties by means of the numerical techniques as introduced in the next section.

1.4 Numerical methodology

Contour dynamics and point vortex methods are theoretical tools that are ap-plied to great success for solving various two-dimensional flow problems in fluid dynamics and plasma physics. For fluid systems both methods are used for highly idealized physical systems containing vorticity, and both methods are characterized by the fact that the vorticity sources (patches with uniform vor-ticity or point vortices) are advected by their own induced velocity field. A

(15)

substantial advantage of these gridless Lagrangian schemes is that the circula-tion is conserved when accurate time-integracircula-tion schemes are used. Small-scale flow properties (vorticity filaments or passive tracers) can be followed for rather long times in an efficient and accurate manner.

Contour dynamics describes the evolution of finite areas (patches) of con-served uniform vorticity. The ‘contour’ in contour dynamics refers to the bound-ing curve of the vorticity patch and is the prescribbound-ing entity for the evolution of the vorticity distribution. Each contour is on its turn discretized by means of many nodes connected by elements. The method was already in use in the late 1970s by Zabusky et al. [79] for solving the 2D Euler equations, while even in the late 1960s a similar method—where contour dynamics actually evolved from—was in use in plasma physics modelling (Berk & Roberts [14]). In the late 1980s Dritschel contributed extensively in modifying and applying the method [28, 29, 30] (see also the work of Pullin [58]). While inviscid vortex flows can be simulated with contour dynamics when approximations like the spatial discretization of vorticity distributions are allowed, the simulation of vortices in a two-dimensional incompressible (in)viscid fluid can be realized by means of other computational techniques as well; for instance finite-difference and finite-volume methods, and finite-element methods. These so-called Eule-rian methods require a grid to compute the appropriate quantities. In order to accurately follow the flow quantities in time, high resolution meshes for the flow domain are necessary to resolve the small-scale vorticity filaments. Euler codes, however, overrun their own resolution in a very few eddy turnover times, implying the need of higher resolutions in the course of time. The contour dy-namics method is in favour compared to the above mentioned grid-dependent examples due to its Lagrangian character.

The geophysical flows mentioned earlier are simulated by means of the con-tour dynamics method. We distinguish planar concon-tour dynamics and spherical contour dynamics. Planar contour dynamics simulations are conducted in an unbounded 2D Cartesian geometry and in order to model the planetary back-ground vorticity properly, transformations from the spherical to the planar geometry are necessary. If, however, the flow dynamics is only relevant for rel-atively small latitudinal regions, approximations to the prescribing background vorticity are appropriate. Furthermore, the contour dynamics algorithm can be modified conveniently in order to compute the vortex dynamics in a three-dimensional infinitely thin shell by simply adding an extra axial coordinate in the Cartesian geometry (Dritschel [26]). This spherical contour dynamics scheme is then capable of capturing all of the appropriate dynamics in the thin shell of fluid on a rotating sphere.

The incorporation of the background vorticity field implies the use of many contours in order to describe the evolution of the total vorticity field properly. An increase in the number of contours demands an increase in contour nodes, and especially the highly intricate patterns that develop in the flow require a

(16)

1.4 | Numerical methodology 7

huge amount of computational power when the flow needs to be followed for longer evolutionary times. Acceleration efforts have been implemented suc-cessfully for contour dynamics in general by proposing a hierarchical-element method (HEM), see Vosbeek et al. [72]. The HEM had been developed in order to speed-up serial-machine computations. The advent of parallel computing architectures, however, did present the question whether the deployment of multiple CPUs would be appropriate. Concurrently solving the contour dy-namics equations for specific hard-to-solve problems can be done on a parallel architecture like the shared-memory Origin 3800 platform TERAS at SARA Amsterdam. The parallelization of the HEM is an important part in this thesis and turns out to be an excellent means to simulate vortex flows with non-uniform background vorticity.

Another means for simulating vortex flows via Lagrangian techniques is through the point vortex dynamics equations. Point vortices are the singu-lar limit of the finite vorticity patches of contour dynamics in their geometric distribution of vorticity, be it that the circulation has a predefined constant value. This singular (non-physical) character appeared to be a convenient and simple means for modelling various general vortex flows, with investigations starting with Helmholtz in the 1850s [38]. A few decades later, before quantum physics appeared on the scene, stationary point vortex patterns where thought of as a model for atoms and molecules (Thomson [70]). Towards the end of the 20th century, point vortex studies gained interest in the light of chaotic fluid flow investigations (see, for example, Aref [7], Aref et al. [9] for pioneering results). Physical vortices which come close to the point vortex description

Figure 1.4 — Results of experiments with magnetized electron columns in a gradient magnetic background ﬁeld. Reprinted with permission of Schecter et al. [61].

in a more contemporary experimental setting can be found in Bose-Einstein condensates (e.g. superfluid He-II) and in superconductors—in case of a so-called mixed-state type-II superconductor magnetic flux lines that penetrate the superconducting material can actually be modelled as line vortices (see, for example, Brandt [17]). Similar patterns of flux lines have also been observed in a so-called Penning-Malmberg trap used in plasma physics where magnetized electron columns in a magnetic confinement are modelled as extremely small

(17)

line vortices in a magnetic background ﬁeld (see e.g. Schecter et al. [61] or Durkin et al. [33]). The patterns (vortex crystals)—see Figure 1.4—that ap-pear through self-organization in such a circular symmetric 2D ﬂow setting can be described by analytical and numerical techniques that are mathematically equivalent to the 2D Euler equations.

The ﬂuid equivalent for the above patterns can be neatly simulated by means of a combined implementation of the contour dynamics method and the point vortex dynamics equations—dubbed Point Vortex Contour Dynamics (PVCD). The patterns of Figure 1.4, for example, are reproduced easily by embedding point vortices in a circular patch of uniform vorticity in a polygonal manner. Moreover, the non-trivial shaped multipolar equilibria, as introduced in Section 1.3, are perfect candidates for investigation by means of the PVCD.

1.5 Thesis outline

An overview of the theory of two-dimensional vortex flows is given in the next chapter. Contour dynamics, point vortex methods, and the PVCD are outlined briefly in Chapter 3. In Chapter 4 the multipolar equilibria as introduced in Section 1.3 are discussed and verified by means of the PVCD scheme. The exact sub-class of solutions to the Euler equations is given and several test-cases are reviewed when stability and advection properties of some of these equilibria are investigated. The vortex flows of these equilibria do not require excessive CPU power and can be solved on a single-processor machine in a reasonable amount of time. The idealized geophysical flows that are simulated by means of large numbers of contours require, however, an acceleration of the method. In Chapter 5 we discuss the HEM and investigate the parallelization of the HEM extensively. The details of the parallelization strategy, the performance, and efficiency with the aid of various examples, are discussed. In Chapter 6 we investigate passive tracer advection in simplified models of vortex flows which are valid on a rotating planet. We focus on the evolution of a disturbed idealized polar vortex and the advection of material tracers at the edge of the polar vortex. Finally, concluding remarks and recommendations are given in Chapter 7.

(18)

2 Theory of two-dimensional vortex flows

The theoretical aspects of two-dimensional inviscid incom-pressible vortex flows are presented briefly in this chapter. The 2D Euler equations are derived in a specific form appro-priate for the numerical implementation discussed in Chap-ter 3. Options are outlined for the incorporation of the planetary background vorticity, which is characteristic for geophysical flows.

(19)

2.1 The 2D Euler equation

The flow under consideration is incompressible, inviscid, and given for a two-dimensional infinitely extended domain. For a 2D flow we have the Cartesian coordinate vector x = (x, y)T, the velocity field u(x, t) = (u, v)T, and the density ρ(x, t) as the prime kinematical variables. The conservation of mass is expressed in the continuity equation

∂ρ

∂t + u· ∇ρ + ρ∇ · u = 0 , (2.1) which requires that mass is neither created nor destroyed in a fluid regardless whether the fluid is compressible or incompressible. When besides mass also the volume of a small fluid parcel is preserved during convection in the flow, the material derivative of the density becomes

Dρ Dt ≡

∂ρ

∂t + u· ∇ρ = 0 , (2.2) which is the deﬁning property for incompressibility of the ﬂow. The continuity equation (2.1) now reduces to

∇ · u = 0 , (2.3)

characterizing the divergence-free velocity field. A fluid parcel in an incom-pressible 2D flow can translate, rotate, and deform (in an area-preserving way), but not expand or contract. Moreover, the density in the flow field does not have to be uniform as long as the density of each individual infinitesimal fluid particle remains constant during the motion in the flow, as is the case when (2.2) applies. We further consider a barotropic flow, i.e. a flow for which the iso-density planes coincide with the isobaric planes and thus p = p(ρ).

The kinematical properties discussed above are to be complemented with the dynamical properties of the flow field variables when forces develop as a result of the motion in the flow. The resulting equation of motion is called the Euler equation which is written as

∂u

∂t + u· ∇u = − 1

ρ∇p , (2.4)

when external forces (like gravity) are absent and where p(x, t) denotes the pressure. This well-known inviscid variant of the Navier-Stokes equation in two dimensions for an incompressible flow gives the balance of linear momentum with the variables u, ρ, and p. A more useful form applied to vortex flows is obtained through the definition of the vorticity:

ω(x, t) ≡ ∇ × u = ∂v ∂x − ∂u ∂y e_z . (2.5)

(20)

2.1 | The 2D Euler equation 11

The vorticity has its only non-zero component perpendicular to the 2D ﬂow domain, i.e. ω = ω e_z for u = (u, v)T. We write the momentum conservation equation (2.4) as

∂u

∂t +ω × u = −∇k −

1

ρ∇p , (2.6)

where the vector identity u· ∇u = (∇ × u) × u + 1₂∇(u · u) is used, and hence u· ∇u = ω × u + ∇k with k the kinetic energy per unit mass. Applying the curl to (2.6) and bearing a barotropic ﬂow in mind, yields

∂ω

∂t +∇ × (ω × u) = 0 , (2.7) which on its turn is recast into

∂ω

∂t + u· ∇ω + ω∇ · u = 0 (2.8) for the scalar vorticity ω. The incompressibility condition (2.3) implies vorticity conservation; and the advection of the scalar vorticity ω with the ﬂow ﬁeld is given through

Dω

Dt = 0 . (2.9)

When vorticity is initially present in an inviscid flow it will remain there for indefinite times. The same condition (2.3) is used to define a stream function ψ(x, t) through

u≡ ∂ψ

∂y, v≡ − ∂ψ

∂x . (2.10)

The equation of motion for the vortex ﬂow is now written in terms of vorticity and stream function as

∂ω ∂t + ∂ψ ∂y ∂ω ∂x − ∂ψ ∂x ∂ω ∂y = 0 . (2.11)

The 2D Euler equation has been recast into a formulation with the vorticity ω and the stream function ψ as prime variables. An advantage of this repre-sentation is that numerically one does not have to worry about the continuity equation separately. Another advantage is that the pressure does not appear in case of boundary condition issues. In this study the inﬁnite two-dimensional domain (planar geometry) as well as the shell of a rotating sphere (spherical geometry) is considered.

Through the deﬁnition of the vorticity (2.5) and the stream function (2.10) we observe that these two quantities are related through the Poisson equation

(21)

Given a finite vorticity distribution ω(x, t) in the two-dimensional domain, one can solve (2.12) for ψ(x, t)—with appropriate boundary conditions—and hence obtain the induced velocity field u(x, t). In Chapter 3 we solve for ψ and show how an expression is derived for which the velocity components are computed and show that the advection of the vorticity field can be followed rather accu-rately and efficiently.

2.2 Passive advection

Advection is the transfer of fluid properties by a prescribed velocity field. Ad-vection of passive tracers is in general described by the adAd-vection-diffusion equa-tion, which is written as

∂c

∂t + u· ∇c = κ∇

2_{c + S(c) ,} _(2.13)

where c(x, y, t) is the tracer concentration, κ is the molecular diffusivity of the scalar c, and S(c) is a sink/source contribution, e.g. a chemical reaction term. If the molecular diffusion and the extra sink/source are absent, the advection of the tracer particles is solely determined by u(x, t). The position of each tracer particle is then followed through the Lagrangian description of fluid motion and obtained through

dx

dt = u(x, t) , (2.14)

which is a set of ordinary differential equations. Passive particles added to the fluid are assumed not to interact with the flow in any way. The fluid velocity u(x, t), as given by the Euler equation, determines the velocity of the passive particles, and likewise of the fluid parcels themselves. In the former section the stream function is defined through the velocity field u(x, t) for incompressible flows. The kinematical relation between the particle equations of motion and the stream function now becomes

dx dt = ∂ψ ∂y , dy dt =− ∂ψ ∂x , (2.15)

where we recognize ψ as the Hamiltonian H, and the position (x, y) of the particle as the canonical pair (q, p). In two dimensions a direct link between the phase space of Hamiltonian systems and the physical space of ﬂuid kinematics is apparent, which is fortunate in visualizing the ﬂow properties.

The above description is rather simple for two-dimensional steady flows for which the above equations are fully integrable and indeed show simple advection of passive tracers along streamlines. The time-dependence of the velocity field alone is, however, sufficient for creating irregular particle trajectories, i.e. a fully deterministic velocity field, in the Eulerian sense, can exhibit non-predictable

(22)

2.3 | Geophysical ﬂows 13

Lagrangian advection characteristics. This behaviour, called chaotic advection, signatures a possible eﬃcient mixing of passive tracers. See, for example, the substantial contributions of Aref [7, 8]). In an ideal 2D ﬂuid (and in many other physical models), the dynamical features like the chaotic trajectories just mentioned and also the stable behaviour of integrable motion, like in coher-ent structures, where closed streamlines might act as transport barriers, often coexist in the same setting.

2.3 Geophysical flows

The atmosphere and the oceans are thin layers of fluid enveloping our planet. The thickness of these layers is very small compared to the scales of the planet itself. When we introduce characteristic velocity scales U and W for horizontal and vertical motion, respectively, and characteristic length scales L and H for horizontal and vertical distances, respectively, then for large-scale atmospheric flows we have to good approximation U = 10 m/s and W = 10−2 m/s and L = 103 km and H = 10 km. Geometrically, this implies an approximately hor-izontal motion of the flow, i.e. locally parallel to the surface of the Earth. The vertical density stratification of the fluid layers (in the ocean via a salinity or temperature gradient and in the atmosphere through a temperature gradient) further constitute a mechanism for two-dimensional behaviour. Once a stable stratification is present in a fluid, it will have a tendency to remain so and vertical displacements are suppressed. The geometry and stratification char-acteristics for two-dimensionality are further complemented by the rotational character of the planet.

2.3.1 Rotational aspects

The equation of motion as given in (2.4) is a direct consequence of Newton’s laws of motion in an inertial reference frame. The reference frame of the Earth is rotating and non-inertial and all motion in a thin layer of ﬂuid (the atmosphere or the oceans) requires a recasting of the momentum equation. Scalar quantities

W

x ( t ) d x = W x d t x ( t + d t )

Figure _{2.1 — The vector x is ﬁxed} in a rotating reference frame.

(23)

and material derivatives of scalars need not to be modiﬁed, yet vector quantities do require a correction to account for this rotation. Consider the rotating frame of reference of the Earth with an angular velocity Ω (Figure 2.1). A position vector x that is ﬁxed in this frame must rotate when viewed in the inertial reference frame. The temporal rate of change of this vector apparent in the inertial reference frame is dx/dt = Ω× x, which is a velocity vector perpendicular to the plane constituted by Ω and x. If the position vector x has a temporal rate of change in the rotating frame, the velocity as observed in the inertial reference frame then becomes

u_i = u + Ω× x . (2.16)

The acceleration as observed in the inertial reference frame is given by du_i dt i = dui dt + Ω× ui , (2.17)

and by incorporating (2.16) we end up with du_i dt i = du dt + 2Ω× u + Ω × (Ω × x) (2.18) for the acceleration present in the inertial frame. The two corrections to the material acceleration are categorized as the Coriolis acceleration 2Ω×u and the centrifugal acceleration Ω× (Ω × x). This latter artificial force corresponding to the centrifugal acceleration of a fluid parcel due to the rotation of the Earth contributes to an effective gravity g e_z =∇Φ , with Φ the geopotential∗ and g pointing perpendicular to the surface of the Earth. See for example Salby [60]. The Coriolis correction 2Ω× u, however, is important when it comes to vortex and wave motion in the horizontal plane with time-scales comparable with that of the Earth’s rotation.

The Euler equation (2.4) now contains the two extra inertial corrections for motion in a rotating reference frame:

∂u

∂t + u· ∇u + 2Ω × u = −g ez− 1

ρ∇p . (2.19)

The introduced characteristic length scale L and a characteristic velocity scale U are now appropriate in order to non-dimensionalize the modiﬁed Euler equation. We obtain

∂u

∂t + Ro (u· ∇u) + 2k × u = −∇P , (2.20)

∗_{The geopotential Φ is the sum of the gravitational potential Φ}_gr_{and the centrifugal}

(24)

where the pressure p and the geopotential Φ are combined into P and where ˜

u = u/U , ˜P = P/ρ Ω U L, ˜t = t Ω, and ˜x = x/L is used. The tilde is omitted in (2.20) and k = Ω/|Ω| is the unit vector in the axial direction. The extra parameter Ro, the Rossby number, is now a measure of the ratio of inertial and Coriolis forces, i.e. a measure of the ratio of the rotational time scale Ω−1 and the advective time scale L/U :

Ro≡ U ΩL ∼

[u· ∇u]

[Ω× u] . (2.21)

The typical values for U and L given earlier show a Rossby number of Ro 10−1. This together with the assumption that the time scale of change of the fluid velocity is small compared to the rotational time scale, i.e. an almost stationary flow, indicates a flow dynamics that is dominated by a balance be-tween the Coriolis acceleration and the pressure gradient force. The equation of motion (2.20) then reduces to

2k× u = −∇P , (2.22)

deﬁning the so-called geostrophic equilibrium†. These geostrophically balanced ﬂows have the property that the velocity is independent of the axial coordinate. This can be seen by applying the curl to the geostrophic balance (2.22), and hence

(k· ∇)u = 0 ⇒ ∂u

∂z = 0 . (2.23)

This result, known as the Taylor-Proudman Theorem, shows that the modelling of the dynamics of large-scale ﬂows in only two dimensions is justiﬁed to good approximation.

In order to derive the expression for the vorticity evolution, we apply the curl to the modiﬁed Euler equation (2.19) and use again the vector identity ∇ × (u · ∇)u = ∇ × (ω × u); the equation for transport of vorticity becomes

∂ω

∂t +∇ × (ω × u) + ∇ × (2Ω × u) = 0 . (2.24) The curvature of the Earth is neglected for horizontal motions with scales much smaller compared to the radius of the Earth. In order to apply the newly obtained Euler equation for ﬂows on a rotating Earth, we therefore consider motion in a Cartesian geometry. The Cartesian plane is situated tangent to the Earth’s surface with the y-axis directed meridionally, the x-axis azimuthally, and the z-axis perpendicular to the surface (Figure 2.2). The rotation vector now has (x, y, z)-components Ω = (0, Ω cos ϕ, Ω sin ϕ) with ϕ the latitude. The

(25)

Coriolis acceleration is then given by 2Ω× u =   −fv+f u −2Ω u cos ϕ   (2.25)

with f = 2Ω sin ϕ the Coriolis parameter which represents the vertical com-ponent of the planetary vorticity. The term−2Ω u cos ϕ is associated with the horizontal component of planetary vorticity and vertical ﬂuid motion. Through dimension analysis, see Pedlosky [52], the term is small compared to the pres-sure gradient and is therefore neglected in (2.24). The large-scale motion of

W 2 W 2 W c o s j 2 W s i n j j z x y

Figure_{2.2 — The planetary} vortic-ity 2Ω and its components.

ﬂuid (air) parcels is predominantly horizontal and f is the prescribing term for the Coriolis acceleration 2Ω× u (see also Salby [60]). With the planetary vorticity as an extra variable, the equation for transport of vorticity becomes

Dω

Dt + u· ∇f = 0 . (2.26)

A striking feature of this equation compared to (2.9) is that the relative vorticity ω is not a conserved quantity anymore. The planetary vorticity f needs to be incorporated for the absolute vorticity q = ω + f to be conserved. The time-independence of f implies equation (2.26) to be written as

D[ω + f ]

Dt = 0 , (2.27)

and thus q is conserved. This implies that a ﬂuid parcel containing a certain amount of vorticity ω and moving relatively to the gradient background rota-tion, say with increasing f , its relative vorticity ω has to decrease in order to conserve q and vice versa. A passive parcel of ﬂuid can become dynamically active in this manner.

2.3.2 Latitudinal approximations

Consider again the Cartesian plane tangent to the Earth as a setting for lo-cal horizontal motions that are signiﬁcantly smaller than the radius of Earth

(26)

(R = 6371 km). With the y-axis directed meridionally, the x-axis directed azimuthally, and the z-axis directed perpendicular to the surface, the Coriolis parameter f is now expanded into a Taylor series

f ≡ 2Ω sin ϕ = 2Ω sin ϕ₀+ δϕ cos ϕ₀−(δϕ) 2 2 sin ϕ0+· · · , (2.28) with only the ﬁrst three terms of the series given. Three approximations centred at ϕ = ϕ₀, and with δϕ = ϕ− ϕ₀, are:

• f-plane approximation: f = 2Ω sin ϕ0 ≡ f0 ,

• β-plane approximation: f = 2Ω sin ϕ0+ 2Ω(y cos ϕ0)/R≡ f0+ βy ,

• γ-plane approximation: f = 2Ω − Ω(x2_{+ y}2_)/R2 _{≡ f}

0− γr2 .

Here δϕ = y/R (β-plane) and δϕ = r/R (γ-plane)—see Figure 2.3. The f -plane approximation is appropriate for ﬂows with a relatively small latitudinal varia-tion—O(102) km or less. The linear variation of f incorporates the zeroth and

j R r j 0 d j W

j

_b

g

R ( a ) ( b ) z x y

Figure 2.3 — (a) The β and the γ plane approximations. (b) The arclength δϕ = r/R describing the Taylor series around ϕ0 at the

North pole.

the first order term of the Taylor series. This so-called β-plane approximation describes flows which have a latitudinal variation on a scale not larger than say a thousand kilometers. For increasing latitude the first order term becomes less significant (and vanishes at the pole), while the second order term becomes the leading term in the variation of the Coriolis parameter. This is the γ-plane approximation and it is valid for flows on a similar scale as the β-γ-plane, yet only for a Taylor expansion around the pole. Both approximations are schematized in Figure 2.3. We note that for a Cartesian reference frame we have x2_{+ y}2 _{= r = Rδϕ for the distance to the pole in case of the γ-plane.}

(27)

2.3.3 Spherical flows in a planar geometry

The γ-plane approximation is found to be (reasonably) adequate for ﬂows in the polar regions of the Earth (from mid-latitudes to the pole). For simula-tions of the dynamics of the polar vortex, and its associated transport proper-ties, this appears to be insuﬃcient and the entire hemisphere (or the complete sphere) has to be taken into account. A proper transformation is necessary when comparisons are made between the planar and spherical geometry. We follow Waugh [76] and Figure 2.4 in this respect and write the transformation as

r =2R(R− z) , (2.29) where r is the radial coordinate in the planar geometry and z = R sin ϕ is the axial coordinate in the spherical geometry (with ϕ the latitude). The radial planar coordinate r = 0 then represents the North pole and r = √2 R the equator. The transformation (2.29) is an area-preserving transformation

im-j

R

z

r '

Figure 2.4 — Area-preserving transformation from spherical to planar coordinates.

plying that the distribution of circulation in the spherical geometry is identical to that in the planar geometry.

Through the transformation (2.29) the planetary background vorticity f = 2Ω sin ϕ is taken into account for all latitudes and takes the form

f = 2Ω 1−1 2 r R ₂ , (2.30)

and resembles the expression for the background vorticity in case of the γ-plane approximation in its quadratic dependence in r. The coordinate r, however, has a diﬀerent meaning in this setting, i.e. the arc length r and the length rare related through r = 2R sin(1₂_Rr). The r in expression (2.29) points to the exact value for f in the spherical geometry, while the r in the γ-plane approximation is determined by the Taylor series around the pole. The two descriptions for the radial dependence become similar when closer to the pole: sin(1₂_Rr)≈ 1₂_Rr.

(28)

3 Contour dynamics and point vortex

methods

The prime numerical technique involved in this investi-gation concerns contour dynamics, which is an efficient and accurate method for simulating the evolution of two-dimensional vortex distributions in Euler flows. A similar Lagrangian scheme is the method of point vortices, which is a long-standing popular numerical technique for simulating general two-dimensional Euler flows. Both methods can be combined in one single scheme for the investigation of spe-cial vortex flows. For geophysical flows with a non-uniform background vorticity, modifications to the generic contour dynamics scheme will be outlined as well.

(29)

3.1 Introduction

Contour dynamics and point vortex methods are accurate and efficient numer-ical techniques for solving the 2D Euler equations. The computed velocity field is used to simulate the advection of all flow properties in time for the appropriate computational domain. In contour dynamics the evolution of the flow is described by one-dimensional contour integrals bounding the patches with uniform vorticity. In the point vortex case a system of ordinary differen-tial equations solves for the velocity induced in the 2D domain. The vorticity patches and the point vortices in these two schemes are sources of vorticity that also behave as active tracers. The induced velocity prescribes the displacement of the computational elements (the contour nodes or the point vortices). If vorticity is taken zero, these same computational elements behave as passive tracers: they are advected by the vorticity sources present, yet induce no veloc-ity field themselves. Furthermore, the inclusion of the point vortex equations into the contour dynamics algorithm yields an excellent means for investigating recently obtained vortex equilibria (Chapter 4).

Both numerical techniques are introduced in this chapter, starting with a short treatment on point vortex dynamics and a discussion on contour dy-namics. The combination of point vortices in a background vorticity ﬁeld is introduced and the implementation of contour dynamics in rotating (geophysi-cal) systems is treated in the light of planar and spherical simulations when a non-uniform background vorticity is involved.

3.2 Point vortices

The velocity due to one point vortex in a 2D unbounded domain at the location (x, y) is given in polar coordinates as

u_r = 0 , u_θ= Γ

2πr , (3.1)

where r = (x− x)2_{+ (y}_{− y}₎2 _{is the distance of the point vortex to the}

location (x, y) where the velocity is computed. Note that the velocity becomes inﬁnite at the location of the point vortex itself. The strength of the point vortex is given by the magnitude of the circulation Γ, which is deﬁned as

Γ≡ Cu· dt = Aω · n dA = Aω dA , (3.2) for a closed curve C (enclosing the point vortex) and where dt is an inﬁnitesimal line element tangent to the curve. The third term shows that by Stokes’s Theorem the line integral is written as a surface integral of an area A (enclosed by C), and n is the normal vector on an element dA. The circulation along a

(30)

3.2 | Point vortices 21

curve C can thus be seen as the flux of vorticity through A, i.e. the circulation is an integral measure of vorticity in the flow. The vorticity associated with the velocity field of a point vortex located at (x, y) is now found to be

ω(x, y) = Γ δ(x− x, y− y) , (3.3) where Dirac’s delta function δ takes care of the inﬁnite vorticity concentrated in a point. This results in a ﬁnite circulation Γ, where the sign of Γ indicates the direction of rotation.

A point vortex will induce a velocity in any other location of a given domain except at its own position (no self-induced velocity). In case of a collection of N point vortices in a two-dimensional unbounded domain, the mutual influence is governed by the following system of 2N first-order differential equations (see e.g. Batchelor [10]): u_i= dxi dt = − N j=1 j=i Γ_j 2π y_i− y_j (x_i− x_j)2_{+ (y}_i− y_j₎2 , v_i = dyi dt = + N j=1 j=i Γ_j 2π x_i− x_j (x_i− x_j)2_{+ (y}_i_{− y}_j₎2 . (3.4)

The strengths of these point vortices are given as Γ_k for k = 1, . . . , N , and

τ = 8 τ = 8 τ = 8 τ = 8 τ = 8 (a) τ = 8 τ = 8 τ = 8 τ = 8 τ = 8 (b)

Figure _{3.1 — (a) General chaotic behaviour for a non-polygonal} shaped system with ﬁve point vortices of equal sign and equal strength. (b) A counter-clockwise rotating polygon of the same point vortices. Black dots are vortices and circles indicate starting posi-tions.

(31)

system of point vortices appears, in general, to be chaotic (Aref [6]), while certain conglomerations of vortices, depending on the symmetry of the initial conﬁguration, can be constructed that show a steady periodic motion (Boatto & Pierrehumbert [15]). Such an array of vortices can have a central vortex or not, while the polygonal shape of an array is essential. An example of an array for which N = 5 is given in Figure 3.1 with in the left picture a group of point vortices with equal strength and sign showing chaotic behaviour. The same vortices behave in a coherent way when ordered in a polygonal setting initially (Figure 3.1b). Note that the non-polygonal setting in Figure 3.1a is just a slight deviation of the polygonal setting in Figure 3.1b. Combinations of diverse sign are also possible, like the intricate patterns studied extensively in Aref et al. [9].

3.3 Contour dynamics

We mainly follow Pozrikidis [57] when we consider a bounded region in an irrotational two-dimensional ﬂuid ﬂow as depicted in Figure 3.2. The area A consists of N small circular parcels with vorticity. Each parcel has area ∆A_i and vorticity ω_i for i = 1, . . . , N . The circulation of each parcel is Γ_i= ω_i∆A_i. Replacing the parcels with point vortices and using (3.4), the velocity induced

w

i

A

Figure 3.2 — N small parcels with vorticity comprise a bounded region. For N → ∞ the patch consists of point vortices only (right picture).

by these vorticity patches becomes

u(x, y) = − 1 2π N i=1 y− y_i (x− x_i)2_{+ (y}− y_i₎2ωi∆Ai , v(x, y) = + 1 2π N i=1 x− x_i (x− x_i)2_{+ (y}− y_i₎2ωi∆Ai , (3.5)

(32)

3.3 | Contour dynamics 23

for x_i = x and y_i = y. For N → ∞ these summations transform into the area integrals u(x, y) = − 1 2π A y− y (x− x)2_{+ (y}− y₎2ω(x, y) dA , v(x, y) = + 1 2π A x− x (x− x)2_{+ (y}− y₎2ω(x, y) dA , (3.6)

and form the two components of the Biot-Savart vector expression for two dimensions: u(x, t) =− 1 2π A r× ω(x, t) r2 dA _, _(3.7)

which can be recast into

u(x, t) =− 1

2π

A∇(ln r) × ω(x

_{, t) dA} _, _(3.8)

where r≡ x − x and r =|r|. This expression gives the induced velocity due to a vorticity distribution.

The velocity vector u satisﬁes the incompressibility condition (2.3) and is therefore divergence free, and can be written as the curl of another vector function: u = ∇ × ψ, with ψ a vector potential, which on its turn can be made divergence free through Helmholtz’s Theorem. The vorticity as deﬁned in Chapter 2 is now written asω = ∇ × (∇ × ψ), and using the vector identity

∇2_{ψ = ∇(∇ · ψ) − ∇ × (∇ × ψ), the vectorial Poisson equation}

∇2_{ψ = −ω} _(3.9)

is obtained with solution

ψ(x, t) = −

AG(r)ω(x

_{, t) dA}_, _(3.10)

for the source ω. Here we anticipate on a two-dimensional integration, and hence require

G(r) = 1

2πln r

ψ = ψ ez (3.11)

ω = ω ez ,

with G(r) Green’s function for the inﬁnite 2D domain (see Marshall [45]). The vector potential ψ is related to the stream function and is oriented normal to the plane. Likewise, the vorticity has its only component in the z-direction. Note that the curl of solution (3.10) gives the Biot-Savart expression (3.8).

(33)

Putting (3.11) into (3.10) yields, not surprisingly, the solution for the Poisson equation (2.12) of Chapter 2: ψ(x, t) =− 1 2π Aln(r) ω(x _{, t) dA} _. _(3.12)

For a distribution of uniform vorticity the stream function becomes ψ(x, t) =− ωˆ

2π

Aln(r) dA

_. _(3.13)

Now, by taking the x and y derivatives of this stream function and subsequently applying Stokes’s Theorem for a scalar ﬁeld, the 2D area integral is recast into a 1D contour integral which describes the velocity ﬁeld for the entire domain:

u(x, t) =−ωˆ

2π

Cln(r) dt

_, _(3.14)

where dt is the inﬁnitesimal vector element tangent to the patch boundary. Expression (3.14) shows that the evolution of a patch of uniform vorticity is fully determined by the evolution of its bounding contour.

3.3.1 Spatial discretization

A single patch of uniform vorticity is often not sufficient for simulating inter-esting 2D vortex flows. The continuously distributed vorticity field of more realistic flows is approximated by nesting top-hat distributions of uniform vor-ticity in a way as schematized in Figure 3.3. The vorvor-ticity ˆω is then to be

P 1 ( a ) P 2 P 3 P 4 P 5 P 0 w 2 w 1 w 3 w 4 w 5 ( b )

Figure_{3.3 — (a) Five nested patches with uniform vorticity in plan} view. (b) The approximation of the continuous vorticity distribution.

replaced by a summation of uniform vorticities, namely ˆ ω(x, 0) = m k=0 ω_k, x∈ P_m(0)\P_m+1(0), for m = 0, . . . , M . (3.15)

(34)

Note thatP_m+1(0)⊂ P_m(0), whereP₀(0) =R2 andP_M+1(0) =∅. T he vortic-ity in the background is assumed to be zero for now (ω₀ = 0), and expression (3.14) becomes u(x, t) =− 1 2π M m=1 ω_m Cm(t) ln(r) dt . (3.16) Due to the conservation of vorticity, the piecewise-uniform distributed patches of vorticity remains piecewise-uniform in the course of time. The patchesP_m(t) are likely to deform during the flow evolution, whereas the area of the vortex is conserved and the bounding contours C_m(t) will not cut neighbouring bound-ary contours. It is then possible that strong vorticity gradients and long thin filaments develop in some regions of the flow, which can be followed with rather high precision.

Contour interpolation

The spatial discretization given above deals with the vorticity distribution. In order to compute the contour integrals we have to discretize the contours as well. The contours are approximated by a ﬁnite number of nodes and elements connecting these nodes. The elements can be described by linear, quadratic, or higher order interpolation functions. In case of N linear elements the parame-terization x_n(ζ) for n = 1, . . . , N is given as

x_n(ζ) = 1

2(1− ζ) xn+ 1

2(1 + ζ) xn+1 , (3.17)

where x_n(−1) = x_n, x_n(1) = x_n+1, and x_N(1) = x₁(−1). The velocity as given in (3.16) can now be written as

u(x, t) =− 1 2π M m=1 ω_m N n=1 1 −1 ln[r_n(ζ)]∂xn ∂ζ dζ , (3.18)

where r_n(ζ) ≡ |x − x_n(ζ)|. Each element is integrated by means of Gaussian quadrature. Note that the integrand becomes singular if x approaches x_n. In that case an analytic computation is required for a certain threshold.

Small interpolation errors are introduced when contours are approximated. The evolving flow changes the shape of the contours into more and more com-plex structures. By rearranging the nodes in a strategic manner, the error— which changes likewise with the flow—can be kept sufficiently small and is of the order of its initial magnitude if an accurate time-integration method is used. Vosbeek [71] has proved that, in case of linear element interpolation as imple-mented in this research, the overall discretization error is O(h2

max), with hmax

(35)

Node redistribution

In most ﬂow simulations intricate vorticity patterns are formed in the course of time. Consequently, some pieces of the contour are strongly curved, while other pieces have almost no curvature—see the illustration in Figure 3.4. Several

Figure 3.4 — Contour dynamics demonstration of the node redistri-bution. Strongly curved regions on a contour have a higher node den-sity. Weakly curved regions have a sparser node density.

particular curvature thresholds are used to determine the precise distribution of the nodes on the contour. Four criteria that are employed in the present contour dynamics algorithm are:

1. Approximations of the contour by means of node addition in case of stronger curvature, established via a polynomial ﬁt to determine the local curvature. Node removal is applied to regions where the curvature is low. 2. Two neighbouring contours are not allowed to intersect. This might hap-pen if a small scale (highly curved) feature approaches a much less curved long piece of contour (its own contour or from another one). This cross-ing is prevented by addcross-ing a node at the contour element of the less curved contour when a certain critical distance between the node and the element is measured.

3. Restrictions on the element length requires that an element is not allowed to be made as small as possible in view of CPU requirements. To prevent extreme long elements when straight ﬁlaments are formed, a maximum length is adopted as well.

4. Quasi-uniformity ensures that element lengths gradually change.

The reader is referred to Vosbeek [71] for an in-depth analysis—including a full account of the error estimations—regarding the contour discretization and interpolation.

The number of nodes plays an essential role because the integration of each contour element between two nodes yields a velocity contribution for every other node in the domain, implying a temporal-complexity ofO(N2) for a total

(36)

of N nodes. During a simulation the number of contour nodes N may grow in time due to the increasing complexity of the piecewise-uniform vorticity distribution. For highly intricate ﬂow patterns, and hence huge numbers of nodes, the conventional contour dynamics approach is of limited use, and an alternative implementation of the contour dynamics method is required. Such a scheme, a hierarchical-element method for contour dynamics simulations, is introduced in Chapter 5.

3.3.2 Temporal discretization

The Hamiltonian character of the dynamics of the ﬂow implies the conservation of the area of the vorticity patch. The spatial discretization does not change this feature (Ref. [71]). Explicit time integration schemes such as the forward (or backward) Euler schemes and the explicit Runge-Kutta schemes will, however, not preserve the area of vorticity. The error size solely depends on the required computation time and resolution depth of these techniques.

In order to ensure area conservation for complex and long ﬂow simulations, a symplectic time-integration scheme is required. Here, a combination of an explicit part (predictor) and an implicit part (corrector) is implemented: First, an Euler forward step is used for a ﬁrst approximation of new node positions. Then for a number of cycles the velocity is evaluated and new positions are corrected (midpoint rule). In the last evaluation step velocities are calculated at the latest approximated positions of the nodes. The velocities will be used for the next time interval.

In Vosbeek & Mattheij [74] a thorough demonstration of this second-order symplectic Runge-Kutta scheme is given.

3.3.3 Example – A shielded Rankine vortex

A small contour dynamics example is appropriate in view of the study presented in Chapter 4. The vortex distributions investigated in that chapter comprise combinations of patches of uniform vorticity with special characteristics. One of the features is that the velocity induced outside the vortex is zero.

The vortex consisting of a circular patch of uniform vorticity, the so-called Rankine vortex, is one of the simplest vorticity distributions. It shows a rigid rotation of the ﬂuid inside the patch, while outside the vortex the velocity decreases as 1/r. When a second Rankine patch of radius r₂ with vorticity ω₂is concentrically placed inside a larger Rankine patch with radius r₁ and vorticity ω₁, such that the total circulation is zero, i.e. ω₁πr2₁+ ω₂πr₂2 = 0, the obtained vortex is a so-called shielded Rankine vortex. A schematic representation is given in Figure 3.5 for a positive vorticity ω₁ and negative vorticity ω₂. T he

(37)

azimuthal velocity ﬁeld is given by u_θ=                  ω₁r 2 + ω₂r 2 r < r2 ω₁r 2 + ω₂r2₂ 2r r2 < r < r1 0 r > r₁ , (3.19)

for r₁ > r₂, |ω₂| > |ω₁|, Γ₁ = −Γ₂, and u_r = 0. As an illustration, the evolution of this vortex has been simulated when a perturbation is applied to its outer contour. Similar perturbations have been applied to the special vortex

r1

r 2

w

w 2

1

Figure 3.5 — Two nested Rankine vortices form one shielded Rankine vortex if the circulations cancel each other.

distributions introduced in Chapter 4 to investigate their stability properties. A large mode-3 perturbation—a sine with amplitude ofO(r₁/10)—transforms

τ = 0 τ = 4 τ = 8

Figure3.6 — Contour dynamics demonstration of a shielded Rank-ine vortex with a large mode-3 perturbation. The resulting stable triangular vortical structure is steadily rotating clockwise according to its core sign vorticity.

the monopolar shielded Rankine vortex into a steady rotating equilibrium of one central region of positive sign surrounded by three satellite-like regions of negative sign. Three frame-shots are given in Figure 3.6.

(38)

3.4 | Contour dynamics including point vortices 29

The inner patch of the shielded Rankine vortex can be shrunk to a very small patch. By keeping its circulation Γ = ωA constant in order to maintain a shielded system, the vorticity should be larger for the smaller area A of the patch. A large value for the vorticity induces large velocities in the vicinity of the small patch. The time step in the integration scheme then should be properly adjusted (be made smaller) in order to yield accurate results in the evolution. In Chapter 4 we will introduce a novel family of multiple-patch equi-libria where above issues are at stake. For the numerical experiments performed in that chapter a proper time integration is provided for.

3.4 Contour dynamics including point vortices

The generic contour dynamics method described in the former sections was developed several years ago by Vosbeek [71]. The inclusion of the point vor-tex equations in the present contour dynamics algorithm yields a method for simulating flows of patches and point vortices with the same accuracy as ordi-nary contour dynamics simulations. The algorithm differs in the sense that not only contour-to-contour interactions are involved but three more interactions as well, i.e. contour-to-point vortex interactions (the induced velocity at the point vortex location due to the patch vorticity), point vortex-to-contour interactions (the induced velocity at all contour nodes due to the point vortex strength), and point vortex-to-point vortex interactions. A point vortex is implemented as a contour with only one node having a finite circulation. The induced velocity field through the combinations of these two approaches is still a full solution of the Euler equations. Furthermore, because of the singular character of the point vortex, the initial conditions of the flow at hand should be such that nodes of contours should not be placed too close to point vortices. Obviously, one should also avoid placing two or more point vortices too close to each other. The time step for the time integration scheme should be taken very small in those cases, yielding inefficient simulations.

The new algorithm is dubbed Point Vortex Contour Dynamics or PVCD for short. We close this section with a PVCD simulation example of a shielded Rankine vortex consisting of a point vortex and a patch of uniform vorticity.

3.4.1 Example – A patch-point shielded Rankine vortex A simple example of a PVCD simulation is appropriate in the light of the special vortex equilibria investigated in Chapter 4. Like the example in Subsection 3.3.3 the feature of non-zero velocity outside the vortex (thus representing a shielded vortex distribution) is essential. If we shrink the radius of the inner patch of the shielded Rankine vortex in Figure 3.5 to a singular point, we will still have a shielded Rankine vortex if the thus obtained point vortex has strength Γ_PV = −ω₁πr2

1. The velocity is zero at and outside its boundary but with