A variationally optimized vortex tracing algorithm for three-dimensional flows around solid bodies

(1)

A VARIATIONALLY OPTIMIZED

_{VORTEX TRACING}

ALGORITHM FOR THREE-DIMENSIONAL

_FLOWS

AROUND SOLID BODIES

TECHNIStIIE WIIVERSITEIT

Labum

Archief

MaI)weg 2, 2628 CD Detft

(2)

A VARIATIONALLY OPTIMIZED VORTEX TRACING ALGORITHM FOR THREE-DIMENSIONAL FLOWS AROUND SOLID BODIES

(3)

A VARIATIONALLY OPTIMIZßD VORTEX TRACING ALGORITHM

FOR THREE-DI}{ENSIONAL FLOWS AROUND SOLID BODIES

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Deift,

op gezag van de Rector Nagnificus, Prof. Dr. J.M. Dirken, in het openbaar te verdedigen ten overstaan van een coinmissie

door het College van Dekanen daartoe aangewezen, op dinsdag 14 juni 1988

te 14.00 uur

door

Jacobus Jozef Wilheiwus van der Vegt

geboren te Utrecht,

(4)

Dit proefschrift is goedgekeurd door de promotor:

(5)

(6)

CONTENTS Page

INTRODUCTION 1-1

VORTEX MODEL FOR THREE-DIMENSIONAL FLOWS 2-1

2.1. Vortex models 2-1

2.2. Basic equations of vortex motion 2-3

2.3. Action Principle for the inviscid

vorticity transport equation 2-10

2.4. Hamiltonian functional and Dirac bracket for the

flow map of an inviscid incompressible flow 2-21 2.5. Definition and proof of convergence of a

product formula for the solution of the

viscous vorticity transport equation 2-36

REFERENCES (SECTION 2) 2-48

NUMERICAL VORTEX TRACING ALGORITHMS 3-1

3.1. Introduction

31

3.2. Inviscid flow map for two-dimensional flow fields 3-6

3.3. Inviscid flow map for a three-dimensional flow field

around a circular cylinder 3-13

REFERENCES (SECTION 3) 3-26

APPENDIX: Definition of the matrices and vectors used

in the discretized variational formulation 3-29

ALTERNATIVE FORMULATION FOR FORCES AND MOMENTS

ACTING ON A BODY IN A VISCOUS INCOMPRESSIBLE FLOW FIELD 4-1

4.1. Introduction 4-1

4.2. Calculation of forces and moments in a viscous

incompressible flow 4-2

4.3. Calculation of the inertial components of the forces

and moments acting on cylindrical bodies

412

(7)

CONTENTS (continued) Page

RESULTS AND DISCUSSION 5-1

5.1. Introduction

5-]-5.2. Numerical solution of the early stage of the wake

behind an impulsively started circular cylinder 5-1 5.3. Comparison of computed vortex patterns behind a

cylinder in steady flow or oscillating harmonically

with flow visualization experiments 5-J-1

5.4. Forces on cylinders in a uniform flow or oscillating

harmonically 5-17

5.5. Test results of the algorithm for the determination

of the vector stream function for the three-dimensional

flow field around a circular cylinder 5-25

REFERENCES (SECTION 5) 5-33 CONCLUSIONS 6-1 SUMMARY 7-1 SUMMARY (DUTCH) 8-1 CURRICULUM VITAE 9-1 ACKNOWLEDGEMENT 10-1

(8)

STELLINGEN

Door bet direct beschouwen van de deeltjesbanen in een niet-visceuze vioeistof wordt een praktisch beter bruikbare variationele formulering verkregen dan door te trachten de deeltjesbanen in een Euler-formulering toe te voegen door middel van een Lagrange multiplicator (in dit

proef-schrift.)

Berekeningen van drie-dimensionale tijdsafhankeiijke visceuze stromingen geven een zodanige overvloed aan informatie dat ze het inzicht in de

relevante fysische mechanismen niet noodzakeiijk versterken.

Stromingsvisualisatie-experimenten zijn essentiele huipmiddeien bij de ontwikkeling van het fysisch inzicht noodzakeiijk voor de ontwikkeling

van mathematische modellen voor stromingsberekeningen.

Lagrange- en Haniilton-mechanica hebben sen dermate hoog abstractieniveau

bereikt dat het doel hiervan de gebruiker vaak ontgaat. Gezien hun grote voordelen verdient esa uitbreiding van bet toepassingsgebied van deze

methoden de voorkeur.

De inspanning die verricht moet worden orn modellen voor drie-dimensio-nais instationaire visceuze stromingen te ontwikkelen staat niet in verhouding tot bet praktische gebruik van de resultaten.

Het concluderen dat een mathematisch model correct is door vergelijking van de berekeningsresultaten met die verkregen uit experimenten leidt

tot een schijnzekerheid.

De maximale "flop rate" vermeld in de catalogus van een

supercomputer-leverancier blijkt in de praktijk vaak de grootste flop te zija.

Het verbeteren van het ontwerp van een scheepsschroef met verwaarlozing

van de door het schip veroorzaakte vorticiteit getuigt van weinig

consequent handelen.

(J.J.W. van der Vegt, Actuator disk in a two-dimensional non-uniform

flow, International Shipbuilding Progress, 1983).

Het feit dat een rosier een belangrijk deel van zijn voortstuwing tijdens sen roeihaal bereikt door een liftkracht in plaats van een weerstandskracht op bet blad is een sterk veronachtzaamd feit in de

roeisport.

(V. Nolte, Die Effektivität des Ruderschlages, Dissertation

Sporthoch-schule Köln, 1984).

De interesse van de media in het liefdesleven van J.P. Sartre, na het verschijuen van de Nederiands vertaling van zijn biografie, past beter

in een natuurfilm dati in een literaire boekbespreking.

(9)

1.. INTRODUCTION

The investigation of the process of vortex shedding from bluff bodies is one of the classical subjects of fluid dynamics. For more than a century numerous fluid dynamicists have tried to get a clear understanding of

this phenomenon, both experimentally and theoretically. The strong impe-tus for this research is its practical importance. Many structures

suf-fer from vibrations caused by vortex shedding, for instance offshore

structures, which consist of large tubular elements, suspension bridges, cables and chimneys.

Despite the large amount of research carried out into this phenomenon there still are many unsolved problems. Fron the experimental point of view there are problems with the accurate modelling of vortex shedding, partly due to problems with limited attainable Reynolds numbers in the

experiments, partly due to the strong sensitivity of these flows to

small disturbances. Especially the lift forces, the forces orthogonal to the main flow direction, show a large scatter in the experiments. Anoth-er problem in using expAnoth-erimental results is that it is vAnoth-ery difficult to cast the results in simple useful formulas which can be used by engi-neers. Even for one of the most simple cases, the flow around a circular cylinder oscillating harmonically in a fluid at rest, it still is an open question what kind of formulation should be applied in order to

describe the drag and lift force coefficients accurately.

The theoretical prediction of unsteady flows around bluff bodies at mod-erate to high Reynolds numbers presents a tremendous challenge. The dra-matic increase in computing power, however, opens new ways to approach this problem. Nevertheless, a direct solution of the Navier-Stokes

equa-tions for moderate to high Reynolds numbers will be Out of the question.

This stimulates the search for better theoretical formulations, using _as much as possible of the physical insight in the processes governing the evolution of a viscous fluid. Vortex models belong to this class of the-oretical models because they directly model the actual relevant physical

(10)

In this thesis some fundamental questions related to the construction of

vortex models for two and three-dimensional viscous flows around solid bodies will be addressed. In Chapter 2 first a Lagrangianformulation of

the inviscid vorticity transport equation will be discussed using a flow

map. This formulation can be converted into a Hamiltonian description which is used in the definition and proof of convergence of a stochastic algorithm for the solution of the viscous vorticity transport equation. In Chapter 3 two numerical algorithms are discussed, one for the

solu-tion of two-dimensional flow fields and one for the solution of the three-dimensional flow field around a circular cylinder using the

Lagrangian functional derived in Chapter 2. The determination of forces

on a body in a viscous fluid is discussed in Chapter 4. This is a non-trivial problem when the flow field is described in terms of vorticity, because in this case the pressure is eliminated. Finally in Chapter 5

some test results will be discussed. They comprise the calculation of the initial stage of the wake of an impulsively started cylinder, the

vortex formation around cylinders in steady flow or oscillating har-monically and the forces on cylinders in these flows. In addition some test results of the numerical algorithm for the determination of the

vector stream function for the three-dimensional flow field around a

(11)

2. VORTEX MODEL FOR THREE-DIMENSIONAL FLOWS

2.1. Vortex models

A viscous flow field at moderate to high Reynolds numbers is character-ized by the occurrence of vortical structures with a large difference in length scales. The large vortical structures, which are anisotropic in

nature, dominate transport phenomena, while dissipation of vorticity occurs mainly at small length scales where vorticity is isotropic. In the flow is a continuous exchange of energy from the large eddies to

eddies of smaller size and so on. A detailed discussion of this

so-called cascade model can be found in Tennekes et al. [2-1].

A group of theoretical models based on the insight that vorticity gives a natural description of this phenomenon are vortex models. In vortex models the flow field is described by tracking the time evolution of the vorticity field. The method therefore gives a Lagrangian description of

the flow field.

Vortex models for the simulation of viscous flows have already a long history. The first results were obtained by Rosenhead [2-2] for the

Kelvin-Helmholtz instability. Up to now, most of the applications of vortex models are restricted to two-dimensional flows. The extension to

three dimensions is non-trivial and presents a great challenge.

There are several reasons for the restriction to two-dimensional vortex

models:

Vortex models in their original formulation have a very poor

numeri-cal efficiency and accuracy.

The effects of viscosity had to be modelled either by introducing empirical information or results of separate boundary layer

calcula-tions.

In two dimensions the effect of vortex stretching is absent,

(12)

A very thorough review of vortex methods, with special emphasis on

nu-merical accuracy is given in Leonard [2-3].

A theoretical model which aims at the simulation of fully three-dimen-sional viscous and incompressible flows around solid structures

there-fore has to overcome the problems which accompany vortex models. In this

thesis the mathematical formulation of such a model together with some applications is discussed. The method is an extension of the original

fractional step algorithm presented by Chorin [2-4] for 2-D flow and

[2-5] for 3-D flow. A solution of the viscous vorticity transport

equa-tion is obtained by successively solving the inviscid vorticity

trans-port equation and the vorticity diffusion equation for small time inter-vals. The solution is obtained by means of a stochastic simulation representing the random behaviour of the diffusion process in real tur-bulence. Special emphasis has been laid on obtaining an energy conser-vative solution technique for the inviscid vorticity transport equation through the construction of a Lagrangian variational formulation and related Hamiltonjan functional and Dirac bracket. Combined with a spec-tral solution technique this gives an efficient and accurate numerical algorithm. A proof of convergence of the mean of this stochastic process

to the solution of the Navier-Stokes equations is presented together

with estimates of the numerical accuracy.

The theoretical model discussed in this thesis is able to describe in-compressible viscous flow fields at moderate to high Reynolds numbers; but due to the fact that it is not practically possible to resolve all length scales in such a flow, because of limited computer capacity, the present applications are restricted to moderate Reynolds numbers. When accompanied by a turbulence model for small scale structures this

re-striction can be removed.

In the next sections first the mathematical formulation of the vortex model is presented, followed by a discussion of the inviscid vorticity transport equation. A variational formulation of this equation is

(13)

and Dirac bracket for the inviscid vorticity transport equation is dis-cussed in section 2.4. Finally, the fractional step algorithm which solves the viscous vorticity transport equation is discussed in Section

2.5, together with a mathematical proof of the convergence of this algorithm.

2.2. Basic equations of vortex motion

The equations of motion which determine the evolution of a vorticity field in a viscous and incompressible flow are obtained by taking the curl of the Navier-Stokes equations in velocity-pressure formulation,

see for instance Serrin [2-6].

Let u) be a smooth C2 vorticity field, w : D x [t0,t1] + c, in the

ex-terior unbounded fluid domain D Rn, _{(n = 2 or 3) of a body, and let}

[t0,t1] be a time interval in which such smooth solutions exist then the

vorticity field satisfies the viscous vorticity transport equation:

2

- w.Vij = vV a (2.2.1)

Here a = tJ+u : x [t0,t1] + D is a smooth C3 velocity field whose curl

is the vorticity field . The velocity field must be divergence free by

virtue of the continuity equation for incompressible flow. The velocity field a consists of two parts: a time dependent uniform velocity U and a divergence free disturbance velocity field u. The variable y in equation

(2.2.1) represents the kinematic viscosity while V is the gradient oper-ator.

The first two terms on the left-hand side of equation (2.2.1) represent the material derivative of the vorticity field while the third contribu-tion is responsible for the vortex stretching, only present in three-dimensional flows. The contribution on the right-hand side is

(14)

The set of equations for the vorticity field must be supplemented with

the following boundary conditions:

The velocity field ti must satisfy the noslip condition at a body

surface S cQ, with S a C1 manifold of dimension n = i or 2:

ti =0 (2.2.2)

This condition on the velocity field ti means that the vorticity field

at S must satisfy, Serrin [2-6]:

= 0 (2.2.3)

which condition is obtained using Stokes' theorem on an arbitrary curve at the surface S. Here n is the normal vector pointing into 2. Finally

the velocity field ti must approach the uniform velocity field U at great distance from the surface S.

In practice the flow in front of a body always contains an amount of vorticity caused by other obstacles. The only mechanism which produces vorticity in an incompressible fluid without a free surface is the shear stress at the wall, Serrin [2-6], so we can consider uniform flow in front of the body because there is no other source of vorticity in our

problem.

The initial conditions are given by a divergence free velocity field ti(r,t0) with associated vorticity field s(r,t0) u0(r) for each r E Q.

The solution technique for the viscous vorticity transport equation

strongly depends on the solution of the associated inviscid vorticity transport equation. Therefore in the next two sections this problem will be discussed in detail.

The inviscid vorticity transport equation is obtained from equation

(15)

boundary condition, equation (2.2.2) must be replaced by the Neumann

boundary condition:

= 0 _(2.2.4)

stating that there is no flux through the body surface S. The condition for the vorticity field, equation (2.2.3) remains valid, but other

con-ditions are now also possible. They lead to the wellknown models for

lifting surfaces in an inviscid fluid, where the angle of a vortex layer

from the body is given by the Kutta condition. The theory developed in the next two sections can be applied also to this case, but this is not further elucidated. The condition, equation (2.2.3) for the vorticity field, will be used also for the inviscid vorticity transport equation.

This means that if there is rio vorticity in the fluid at initial time t0

it will also be absent for all later times.

One of the important aspects in the solution of both the viscous and inviscid vorticity transport equation is the choice of variables in

which the solution is being sought. As already mentioned in section 2.1 we look for a Lagrangian description of the flow field but the vorticity

transport equation (2.2.1) is given in Eulerian formulation.

The discussion of the transfornation of the equations given in Eulerian formulation into a Lagrangian formulation in the rest of this section is limited to the inviscid case. In section 2.5 the viscous problem will be

further discussed.

The transformation from an Eulerian reference frame into a Lagrangian one can be accomplished by the introduction of a flow map R. Let R be a

invertible transformation:

R(r',t,t) : + O Vt _[t0,t1] _(2.2.5)

with the fluid volume at time t0 and Q the fluid volume at time t, then the initial position r' of a fluid particle at time t0 and its

(16)

position r at time t are given by the relation:

E'

= R(r',t0,t0) (2.2.6)

E = R(r',t,t0) (2.2.7)

If r' is kept fixed while t varies equation (2.2.7) specifies the path of a fluid particle initially at r', while for t fixed it determines the

transformation of the initial fluid domain into its position at time t.

The flow map R is a continuous invertible volume preserving mapping because:

det(R/r') = Vt E[t0,t1] (2.2.8)

by virtue of the continuity equation. For a derivation of this relation,

see Serrin [2-6].

Just like the velocity field Ci the flow map R is separated into two

components: the flow map R which is equal to:

= E'

+ U(t')dt' (2.2.9)

to

and a disturbance flow map R. The total flow map R then can be

repre-sented by:

= o

(2.2.10)

while both flow maps R and R satisfy the condition given by equation

(2.2.8).

The velocity R of a fluid particle is given by the relation:

(17)

together with equation (2.2.6). Here a dot means differentiation with respect to time. This relation couples a Lagrangian description of the flow field to an Eulerian description by means of a flow map R. If we

differentiate R , equation (2.2.9), we obtain the relation:

(r',t,t0) = 0(t) (2.2.12)

so the disturbance flow map R approaches an identity transformation at

great distance from the surface S.

The vorticity field at time t now can be expressed in the flow map R and the initial vorticity field using a result obtained by Cauchy, see

Serrin [2-6]:

w(r,t)

= Vr R(r',t,t0) (2.2.13)

This relation clearly satisfies the initial condition for the vorticity field by virtue of equation (2.2.6). It also satisfies the inviscid

vorticity transport equation, see Serrin [2-6], and changes the problem from solving the vorticity transport equation into determining the flow

map R.

In the introduction, section 2.1, it was stated that a Lagrangian

de-scription of the flow field has some attractive benefits due to its

close relation with the relevant physical phenomena. One of the

disad-vantages is, however, its poor computational efficiency compared to

methods working on a grid fixed in space. The advantages of both methods can be combined using a coupled Eulerian-Lagrangian description of the

flow field through the introduction of a vector stream function A.

Let A be a C2 vector field, A O x [t0,t1] + 0, then the relation between the Eulerian and Lagrangian description of the flow field,

(18)

R(r',t,t0) = V x A(r,t)

= V x A(R(r',t,t0),t) R

(2.2. 14)

using equation (2.2.7). The Eulerian velocity field on the right-hand side of equation (2.2.14) is expressed as the curl of the vector stream function A because in this way it immediately satisfies, just like the

flow map R, the continuity equation div = O.

The introduction of a vector stream function A in equation (2.2.14) also gives a useful relation with the vorticity field w. First, it must be noted that the vector stream function A is arbitrary up to the gradient of a scalar function in equation (2.2.14). This fact can be used to

define a vector stream function A, and a vector stream function A

related to the uniform velocity field U, viz.:

V x = = (2.2.15)

together with a scalar field i1 such that A is divergence free by choos-ing i so that:

V p = div(A - A) (2.2.16)

where:

+VIp (2.2.17)

As a consequence of this separation both A and V approach a zero vector

at great distance from the surface S.

The relation between the vorticity field and the vector stream function now can be expressed as:

(19)

= V x V x (2.2.18) k(r',t,t0) = V _{xA(R(r',t,t0),tJ} R V2A(r,t) = - w(r,t) with: R(r',t,t0) = R R(r',t,t0)

The vector field A must be divergence free and satisfy the following

boundary condition at S:

= .V x

_(2.2.24)

together with:

VxAO

at infinity _(2.2.25)

The disturbance flow map R must have a Jacobian determinant _{unity and} approach the identity mapping at great distance from _{S. These conditions}

must be supplemented with an initial divergence _{free vorticity field} ui(r') which satisfies equation (2.2.3).

= - _(') _{Vrt R(r,t,t0)} (2.2.20) (2.2.21) (2.2.22) (2.2.23)

= V2A

_(2.2.19)

using the fact that the vorticity field is defined as the curl of the velocity field and the velocity field U is uniform in space.

The introduction of the flow map and the vector stream function_{gives an} alternative way to describe the evolution of an inviscid vorticity field. In the next section an Action Principle will be presented _which

generates the equations for the flow map R and vector stream function A. These equations can be summarized as:

(20)

2.3. Action PrinciEle for the inviscid vorticit transort equation

The inviscid vorticity transport equation possesses certain invariance propertiea, such as conservation of energy and Kelvin's theorem on vor-tex motion. A numerical algorithm which attempts to solve this equation

properly should satisfy these constraints as closely as possible.

Lagrangian and Hamiltonian mechanics provide an elegant framework for the construction of such models. Due to their success in mechanics there is already a long history of attempts to apply them to inviscid fluid

mechanics.

The success of this approach was, however, limited. The first results

were only applicable to irrotational flow. The extension of the

Lagrangian functional for potential flow to inviscid flow with vorticity

was obtained by Herrivel-Lin--Serrin [2-6]. The important extension consisted of the so-called Lin constraint.

A major reason for the problems encountered originates from the fact that there is no generally valid technique for the construction of a Lagrangian functional and success strongly depends on the choice of

variables. The problems encountered when looking for a Lagrangian formulation for general inviscid flow originated from the use of an

Eulerian description of the flow field and thereby loosing the relation

with the individual fluid particle trajectories which are used in

Hamilton's principle. Bretherton [2-7] clearly pointed out this fact by

looking at the relation between variations on particle paths and variations at a point fixed in space.

Serrin et al. [2-6] obtained the equations for inviscid flow by adding several constraints to the energy functional, for instance conservation

of mass and conservation of the identity of particle paths, by means of

Lagrange multipliers. There was, however, a lot of redundancy'in their formulations as shown by Seliger et al. [2-8]. They reduced the number of variables by showing the close relationship of the Lagrange

(21)

The use of Clebsch potentials seems quite natural but these potentials have some serious drawbacks. Bretherton [2-7] demonstrated that they can be multi-valued in certain cases, for instance for knotted vortex lines.

Apart from this problem Clebsch potentials are difficult to use for the

description of a vorticity field. The vorticity field at a certain

posi-tion is given by the relaposi-tion = Vu x V, where u and are constant potential surfaces which follow the local flow.

A more useful Lagrangian formulation can be constructed by looking at the close similarity between a vorticity field and a magnetic field satisfying the Maxwell-Lorentz equations. The Lagrangian formulation for this problem is already known for a long time. Lewis [2-9] gives a

de-tailed discussion of it; further information can also be found in

Goldstein [2-10] and Sudarshan et al. [2-11]. For two-dimensional flow In an unbounded fluid domain with a finite set of point vortices Buneman

[2-12] succeeded in deriving an Action Principle along the lines of

Lewis.

In the remainder of this section we derive an Action Principle for the

evolution of a continuous three-dimensional vorticity field in the

exterior fluid domain of a body. This expression uses s formulation of kinetic energy in terms of velocity and vorticity as can be found in

Serrin [2-6].

Theorem 2.3.1

Let W0 L2(00) be a C° initial vorticity field defined as the curl of the initial velocity field, with L2(00) the Hilbert space aefined on the initial fluid volume °o R3 at time to, and satisfying the boundary

condition = O at a C' manifold S C with unit normal vector n pointing into Let R(r',t,t0) : + O be a C1 invertible flow snap

(22)

hav-by the relation R = R R with R and R the flow map of the undisturbed and disturbed flow respectively. Further we define a C2 vector stream

function A* : O x [t0,t1] + L2(0) as the sum of a divergence free vector

stream function A and the gradient of a scalar function i, together with

a Lagrarigian functional L defined as:

L[R,A*,t] = - ¼< 0(r').V,R(r',t,t0) x R(r',t,t0)

- <w0(r').V,R(r',t,t0) I

A*(R(r',t,t0),t)>

+ ½ <V x A*(r,t) V x A*(r,t)> (2.3.1)

Then the following Action functional J defined as:

ti

J{R,A*,t0,tJ ₌

_f

dt L[R,A*,t] (2.3.2)

to

together with the condition on the vector stream function A* at S,

X =

-.V x

Ç

(2.3.3

with the vector stream function of the undisturbed flow, has a sta-tionary value in for variations ÖR when satisfying equation (2.2.20) and a stationary value in Q for variations A* with t.SA* = O at S, t an

arbitrary tangential vector at S, when satisfying equations (2.2.21) and

(2.2.22).

Here < f 1g>

represents the L2(Q) inner product defined

as

_f

f. dO, while a suffix zero refers to the inner product in the

space L2(00).

Before we can prove this statement we first have to define the gener-alized functional derivative (GFD). Let F[uJ

_{= f}

f(u,x)dx be a

func-tional of some element u E L2(Q) and 6u E C, infinitely differentiable variations of compact support, then the generalized functional

deriva-tive , with u the i-th component of u, is defined as:

u1

= - FEu1 +

e5u1] (2.3.4)

(' ,t,t0)>

(23)

The variations of a Lagrangian functional L with respect to the flow map

R now can be defined as:

6L[R,, t] ÓL

6R>

o o

R constant constant

in space in space

In order to ease the calculation of the variations with respect to the flow map R it is advantageous to separate the part of the Lagrangian functional in the Action Principle which depends on R. In tensor

nota-tion:

L1[R,,t] = ¼ <w(r)R(r,t,t0),.

C

R(r',t,t)

m(r!tt)>

i pmn (2.3.6) and

L2[R,,t] =<e(r)gRP(rt,t,t0),j IA*k(R(rt,t,tO),tJ>

(2.3.7)

with g1 the metric tensor, Cijk = _{i Ei} the permutation symbol (E123

= 1) and g = det(). Here we used the summation convention of tensor

analysis.

The generalized functional derivative of L1 with respect to R is equal

to: /5L1

q\

=

< e(r')R(r',t,t0),

rpmq ,t,t0)

jR>

R constant in space (2.3.8)

This relation is obtained through partial integration and using the

boundary condition imposed on the initial vorticity field which is

divergence free by definition.

(24)

The generalized functional derivative of L1 with respect can be obtained analogously:

The generalized functional derivative of L2 with respect to R is equal

to:

/

L2 0 R constant in space constant in space w(r')R(r',t,t ),. a

Rn(,,0)

O i pqn -(2.3.9)

r )R,g

kp

_{J SRP>}

3R O 2. ijk C C g 'a iqi kp (R » t) 3R

where we used the shorthand notation R for R(r',t,t0).

The generalized functional derivative of L2 with respect to Sft is zero

because L2 does not depend on iL

It is important to note that the generalized functional derivative of L2

with respect to R does not change when a gradient of a C2 scalar

function is added to A*. This will be used when calculating the

varia-tions with respect to changes SA*.

After the calculation of the GFDs of the Lagrangian functional it is now possible to determine the stationary value of the Action Principle for

R constant

in space (2.3.10)

after partial integration, using Ricci's lemma and the initial condi-tions imposed on the vorticity field. The right-hand side in equation (2.3.10) can be further transformed into:

(25)

variations (SR. The first step in the calculation of the stationary value

with respect to (Sft in variations with

is to express the variations

ti

= - ½ f

dt

_f

d3r' a(')R

_tpqr r SRcl

to _Oü

(2.3.12)

where we used the fact that the variations disappear at initial and final time. The second integral on the right-hand side is obtained using the fact that (SR and R must be considered as independent from R in the

process of calculating variations.

Combining the several intermediate results it is now possible to show that the Action Principle obtains a stationary value in for varia-tions (SR if the following relation holds:

C

m()

pmq -A*(R,t) ijk -m = w0(r )Rt(r',t,t0), tfqi j Vr' E Qo (2.3.13)

This equation represents a matrix equation for the unknown fluid

par-ticle velocity with matrix Wqm defined as:

Wqm = - w(') R(r',t,t ),._- _O _i t_pqm Vr' E (2.3.14)

-lt will now be shown that equation (2.3.13) is equivalent to equation

(2.2.20) so that the first part of our theorem is confirmed. respect to SR. (SL

f

dt

_f

d3r'-- (Sq to 0o (SRs ti n

= -

_f

dt

f

d3r'--- (a(r')R. C R J i pqn to 0o R constant in space

(26)

If the matrix Wqm was non-singular we could determine all velocity

com-ponents m by inversion of this matrix. The matrix is, however, singular

so we have to calculate the null-eigen vectors. If the inner product of the right-hand side of equation (2.3.13) and the null-eigen vector is

zero then as many components of Rm are arbitrary as there are null-eigen

values. If this condition is not satisfied possibly more components of can be determined. A detailed discussion of this process can be found in Sudarshan [2-11].

The rank of matrix Wqm is two with eigen values:

A0 andA=±\/-[(w

R.)2

(w R.)2 + (w R.)2)

and the null-eigen vector À is equal to: =

w R.

The inner product of the null-eigen vector and the right-hand side of equation (2.3.13) can now be evaluated and is zero for all time because

it consists of the outer product of X with itself.

Thus the constraint equation is satisfied for all time which means that one of the three velocity components can be chosen completely arbitrary.

The velocity components R1 and R2, for instance, can be expressed in the

right-hand side of equation (2.3.13) and the arbitrary velocity

compo-nent R3 yielding: n

3A(R,)

_{a( )R1('tt0)m}

[3(r',t,t0) (r',t,t0) tjp pq J ₊

i3

A(R,t)

and pq

(27)

t) wm(rt)R2(rt,t,t0), (3(r' ,t,t0) 2(r',t,t0) _Rpq j + w(r')R3(r',t,t0),. 3A(R, t) 3jp g -) (2.3.17) pq

ijk

A*P

where the expression

t

g_kp is equivalent with the i-th component

of curl *

The stationary value of the Action Principle with respect to variations t5R thus is obtained in Où when the velocity components R1 and R2 satisfy

the r1 and r2 components of equation (2.2.20) together with a

contribu-tion depending on the arbitrary component R3 and the r3 _{component of} equation (2.2.20). If we now use the still existing freedom in the

choice of R3 and set R3 at the initial time t0 equal to:

3(r',t,t ) = c3 g _(2.3.18)

- O pq

then this relation is satisfied for all later times and the stationary

value of the Action Principle is obtained when equation (2.2.20) is

satisfied in

The second part of the statement viz., the fact that variations 5A* _of the Action functional give a stationary value of this functional if

equations (2.2.21) and (2.2.22) are satisfied, can be proven analogous-ly.

In the first part of the proof it was stated that the variations _{R of}

the Action functional were invariant for the presence of Vi _{so the} sepa-ration of _{the vector stream function A* in a divergence free vector}

stream function A and the gradient of a scalar function can be used in the variations SA*. The separation of the vector stream function A* in

two independent parta means that the Action functional must be minimized

(28)

* = + V4 = 6A + V (2.3.19)

The variation of the Action functional with respect to A* now yields:

* t1

t

dto d3r'0(r').V,R(r',t,t0).{A(R(r',t,t0),tJ

+ Vip(k(r',t,t0),t)} + Jtdt

_f

d3r(V x A(r,t)).(V x 6A(r,t))

to c

(2.3.20)

It is now necessary to transform the first integral in this expression

from the initial volume to the actual volume Q at time t.

The flow map R and its inverse mapping G relate the initial and actual

positions r' and r to each other:

r = R(r',t,t0) (2.3.21)

r' = G(r,t,t1) (2.3.22)

With the aid of this relation the vorticity field at time t and initial

time t0 can be related to each other:

w(r,t) = s(G(r,t,t1),t0J (2.3.23)

The Jacobian determinant of the transformation from initial flow field into actual flow field and vice versa is unity due to the incompressi-bility constraint. This condition is incorporated in the Action

Prin-ciple, but perhaps not clear at first sight.

Due to the fact that the mappings R and G are each others inverse mapping it can be demonstrated that the product of their Jacobian

(29)

deter-minants is unity, see for instance Ans [2-13]:

det(--J det -) = det(

(G

.2±-)

i

k

Gi Rk

Thus we only have to demonstrate that one of the Jacobian determinants

is unity. Consider now first standard Lagrangian functionals, viz.

functionals for which the passage from a Lagrangian functional to a Hamiltonlan functional and vice versa presents no problems.

If the dynamical trajectory of the flow map

dyn is defined as the

special mapping R so that the Action Principle is stationary for

vari-ations R then a function Sdyn can be defined:

ti S [R(t ),t0,R(t1),t1] ₌

_f

dt L(R,,t) dyn - O t0 along R -dyn

which acts as the generator type I of the Canonical transformation R which relates the flow field at time t1 to the initial flow field, see Sudarshan et al. [2-ii, pp. 63-66]. A Canonical transformation is,

how-ever, just that transformation that preserves the Poisson bracket

rela-tion so that the Jacobian of a Canonical transformarela-tion is always unity, see Sudarshan et al. [2-ii, pp. 73-77].

The Lagrangian functional in equation (2.3.i) is unfortunately of a non-standard type. In section 2.4 it will be demonstrated that also in this

case the flow map R can be generated by a Canonical transformation R estricted to a subspace, defined by a set of weak equations and gener-ated by a suitable Hamiltonian functional HT. Thus also in this case the Jacobian determinant of the flow map and its inverse are unity.

After the transformation from initial to actual volume and partial

inte-gration, the variations of the Action functional with respect to A* become equal to:

(2.3.24)

(30)

ti A* j[R,A*,t0,t1] =

f

dt

f

d3r-w(r,t) + V x V x A(r,t)).A(r,t) - to O ti -

f

dt

f

d2s( x V x A(s,t)).A(s,t) to S ti -

f

dt

f

d2s(n x V x A(s,t)).A(s,t) to S00 (2.3.26)

where we used the boundary condition for the vorticity field, equation

(2.2.3).

The integral along Sc0 the surface surrounding Q at infinity, gives no contribution if V x A + O at Sc0 faster than r2, which is just the far field condition for A. The integral along S would give a boundary

condi-tion for the tangential velocity at S if we allow arbitrary variations

SA in the surface S, while the integral is zero for any variations A normal to S, independent of A. This non-physical boundary condition can be removed by allowing only variations A normal to S while all the

variations in the surface S must be zero. The physical boundary condi-tion then can be introduced as side condicondi-tions to the Accondi-tion funccondi-tional which will be demonstrated in the next pages. The variations A E L2(Q) with t.SA = O at S, t being an arbitrary tangential vector at S, now

give a stationary value for the Action functional if the following

equa-tion is satisfied in Q:

V2A(r,t) = - w(r,t) (2.3.27)

where we used the fact that A is divergence free. This relation is

equivalent to equation (2.2.21).

The only aspect left is to demonstrate that the constraint equation (2.3.3), when added to the Action Principle, does not change this prin-ciple. This can be demonstrated by adding this condition to the Action

(31)

Principle by means of a Lagrange multiplier a and integrating along the

surface S C

ti ti ti

f

dt L*[R,A*,t] ₌

_f

dt L[R,A*,tJ + a

_f

dt

f

d2s n.Vx(A*

+ A)

to to to S

(2.3.28)

But this extension of the Action Principle is trivial because both Lagrangian functionals L* and L are equal which is a direct consequence

of Gauss' identity because:

f

¿s

n. Vx(A* + A) ₌

_f

d3r V.x(A* + A) = 0 (2.3.29)

S G

for any vector A*,A E L2(G).

It is now possible to solve the vector stream function A in the domain G without any reference to the function i' while the variations of the

Action Principle with respect to 5R and 6A* yield equations (2.2.20)

through (2.2.22) with additional boundary conditions.

2.4. Hamiltonian functional and Dirac bracket for the flow ma of

an inviscid incompressible flow

An alternative way to define the evolution equation for the flow map R,

equation (2.2.20), is to express this equation in a Dirac bracket rela-tion together with a _{suitable Hamiltonian functional.} _{These brackets} provide an elegant framework to express invariance and symmetry proper-ties of the flow in a coordinate independent way. Dirac brackets _are elements of a larger class of brackets, Lie brackets, which define _{a Lie}

algebra.

A Lie algebra is a linear vector space, the Hilbert space L2(G) in our

(32)

each element x,y E L2() a third element z = [x,y E L2(Q) and has the

following properties:

The bracket must be antisymmetric

[x,y] = [y,x] (2.4.1)

The bracket must be linear

[Ax + xx',y] = A[x,y] + (2.4.2)

with A, j real numbers.

For any three elements x,y,z E L2(c1) the Jacobi identity holds

[[x,y],z] + [[y,z],x] + [[z,x},y] = o (2.4.3)

The Lie algebra is related to a unique Lie group at least in some neigh-bourhood of the identity, while the Lie bracket is related to the group compositioa law of a Lie group, a rule which associates to each two group elements a third group element. A Lie group is now defined as a set with a group composition law with the properties of associativity and the existence of a unique inverse element for each group element together with an identity element. The operations of taking inverses of given elements and taking products of pairs of group elements must be

both continuous with respect to a given topology while to each element of the group there is a onetoone correspondence with coordinates in an

Euclidian space.

The linearity of the Lie bracket can be used to obtain the structure constants which define a Lie group, while the antisymmetry and the

Jacobi identity of the Lie bracket are a replacement of the

(33)

I Oi pG

z ,z ) = n

I5f() k

</Isf(x)

and the product rule:

1ff2,g} = f1{f2,g} + {f1,g}f2 _(2.4.7)

The class of Poisson brackets can be extended to the so-called general-ized Poisson brackets by considering more general classes of matrices

11U

which are now functions of the phase space variables _{z. The symbolic} notation for these brackets is

{ }*. The Dirac bracket defined as:

mn

(2.4.4)

(2.4.6)

A special class of Lie brackets are Poisson brackets, which are defined

as those brackets which are invariant under Canonical transformations, transformations which do not change the Hamilton equations of motion. A

Poisson bracket can be written as:

ISf(

{f(,g()} =

Is z'1

n

J

with 1U the matrix:

i if '1 n, V ji + n

-1 if V n, = V + n _(2.4.5)

O otherwise

with n = 1, 2 or 3 and 4&

generalized functional derivatives, de-fined in equation (2.3.4) wih respect to the phase space variables, the generalized coordinates and conjugate momenta z. The first three

variables zi are the generalized coordinates qk, _{the rest the}

gener-alized momentum variables

k. In addition to the conditions satisfied by

the Lie bracket the Poisson brackets satisfy the fundamental Poisson bracket rule:

(34)

is a member of the class of generalized Poisson brackets with a singular matrix fl. In this relation the matrix Cmli is a non-singular matrix and the constraints define a subspace of the phase space to which the

variables are confined.

When the equation of motion:

= f(u,t) (2.4.9)

supplemented with an initial field u0 can be expressed as:

= {u,H}* E L1 u (2.4.10)

then this equation is called Hamiltonian and its solution can be written

as a Lie transformation, an infinite sequence of brackets: n

uPexp(tL)ua

_I P n! H u0

n0

with: n p

LH(L1

P

Ldu

=

u)

A detailed survey, especially concerning the group theoretical and to-pological aspects of Lie groups and algebras can be found in Sudarshan

et al. [z-ii] and Abraham et al. [2--14], while their relation to

me-chanics is discussed, for instance, by Goldstein [2-10].

From this short discussion of brackets it will be clear that the for-mulation of a bracket relation and a Haniltonian which describes the evolution equation form an important part of their solution. For a 3-D magnetohydrodynamic flow in Eulerian formulation this bracket has been

revealed by Morrison and Greene, for a review see Morrison [2-15], while

their relation to Clebsch representations can be found in Holm et al. [2-16]. Marsden et al. [2-17] apply these brackets to give a discussion

of the fundamentals of several two-dimensional vortex models, for (2.4.11)

(35)

instance the point vortex method and the vortex blob method, in which

the vorticity is spread Out over some finite core.

In this section a Dirac bracket and Hamiltonian functional will be de-rived from the Lagrangian formulation of the vorticity transport

equa-tion by means of the flow map R, equaequa-tion (2.2.20), using the Lagrangian

functional presented in Section 2.3.

Before we can State this result we first have to define the concept of weak equations. Let f and g be two functions of generalized coordinates and momenta and defined in a finite neighbourhood of a hypersurface U of

the phase space, then the functions f and g are said to be weakly equal, f g, if the values of f and g become equal when their arguments are

restricted to the hypersurface U, see Sudarshan [2-11].

The concept of weak equations offers the opportunity to work on the full

phase space instead of a hypersurface, while we retain constructions like Poisson brackets. During the calculation of partial derivatives of any function on phase space the generalized coordinates and momenta must be considered as independent and the restriction to the hypersurface is

performed at the end of the calculation.

Theorem 2.4.1

Let R(r',t,t0) :

0o

+ D be a

C flow map, satisfying the evolution equation : (r',t,t0) = VR x A*(R(r',t,t0),t) with A*

: o x _[t0,t1] +

L2(D) a C vector stream finction satisfying the condition n.Vx(A* + A = O at S C _{3D, with A} _{the vector stream function of the undisturbed}

flow, then this evolution equation can be expressed by the Dirac bracket

relation: (2.4.13) k with = R,. A*P (R,t) >0

+<

k

\

w1 R. /

(36)

(j = 1, 2 or 3, no summation on j), a Hamiltonian functional, and Dirac bracket defined as:

=

= {f('),g(')} <'m

mn w(r) R(r,t,t0),.

n $)( )}

(2.4.14) (m, n = 1, 2) (j = 1, 2 or 3)

Here and $ : x [t0,t1J + L2(D) define the hypersurface of the

phase space to which the solution of equation (2.4.13) is restricted while the initial C vorticity field E L2(D,) is defined as the curl

of the initial velocity field and satisfies the boundary condition

equation (2.2.3).

The first step in the determination of the Hamiltonian equations of mo-tion, equation (2.4.13), is the calculation of the generalized momentum

variables _k' defined as:

pk,t) =

= ¼ U)(r')

R(r',t,t0),. Cpkn Rn(rt,t,t0) (2.4.15) using equation (2.3.9). Due to the fact that the Lagrangian functional in the Action Principle equation (2.3.1) is linear in the fluid particle velocity R we now obtain the special situation that none of the gen-eralized momentum variables is independent of the generalized coordi-nates R, while the velocity variable is absent from equation (2.4.15). This presents problems when passing from the Lagrangian to the

Hatnil-tonian formalism and vice versa. In the Lagrangian formalism we use the

variables R and , while we use R and in the Hamiltonian formalism. lt

is no longer possible to express all the R variables in the variables when passing from the Lagrangian to the [-lamiltonian framework.

(37)

Instead of a phase space with six independent field variables we now have only three independent variables R1, but in order to be able to work on the original phase space we use the concept of weak equations as

introduced by Dirac.

Using this concept Dirac gave a formulation for the Hamiltonian

equa-tions of motion, which for our field variables R and are equal to:

S(rtt)

_{{RS(rv,t,t0),H} +<{R(rt,t,tû),} (r,t)} (r,t,tü)> (2.4.16) (r',t,t0) p(r',t,t0),H}

_+K{p

(r',t,t0),(r,t)} (r,t,tü)> (2.4.17)

with H the Hamiltonian functional defined as:

H[R] E

_{<2R> - L[R,] = L2[RJ}

(2.4.18)

and L2[R] given by equation (2.3.7). The Hamilton equations must be sup-plemented with the so-called primary constraints, given by the set of

weak equations:

= p(r',t) - ¼ a(r') R(r',t,t0),

t

Rn(,,f,0)

o

(2.4.19)

which define a hypersurface in the phase space. A derivation of these

equations can be found in Sudarshan [2-11, pp. 91-97].

The variables . appearing on the right-hand side of equations (2.4.16) and (2.4.17) are still unknown due to the fact that we have a non-standard Hamiltonian functional. It would appear that we could use

(38)

Generalized Functional Derivatives (GFDs) it can be shown that this

equation is trivially satisfied, viz.: iL

RS(r,t) Rk(x,t)

/óR5(rt)

\pk(,t)

/H

\

t) ¿/' lP('

,t)

t) Kronecker symbol.

The calculation of the CFD in equation (2.4.24) is nontrivial and will be further elucidated. The expression for the constraints $, equation

(2.4.19), can be transformed into:

= P(r',t)

w(r')

R1(r',t,t0) s R'(r',t,tü), (2.4.26)

=0

3) - L)

up

(x,t,t ) s s g -

O 'm lU

pq t) k (3) ( - L') p

with 3)( ) the threedimensional Dirac

5Rk(,t,t _(2.4.22)

(2.4.23)

t, t (2.4.24)

k(,t) (2.4.25)

distribution and 65k the

\pk(,t)

k,t)2

= 0

//

,t)

\

Rk(,t)

(39)

The CFD in equation (2.4.24) now becomes equal to:

-6(' ,t)

\ 6Rk(x,t)

<6

,i

6Rk(,t) ao( 6R1(x,t (a(r') R(r',t,t0),. £ R(r',t,t0)J i ppm ',t,t0) a ppn -6Rk(,t,t

kp

6Rk(,t,t

So far the 1-lamiltonian equations (2.4.16) and (2.4.17) are not very

use-ful because the first equation is trivial and the second depends on the

unknown velocity variables

The constraint equation (2.4.19), however, must be satisfied for all time which relation can be used to determine some of the unknown

vari-ables R:

{0(r',t),H}

0(r' ,t),(r,t)} 1(r,t>

O

(2.4.28)

This relation can be converted into a matrix equation for the unknown

variables using the GFDs for , equations (2.4.24) and (2.4.25): ) R ,t,t0),. ) R1'6,t,t), ) R(x,t,t0), £ppk 6 tkpn 6ppk (3) (2 -- L') 6Rk(, 6R'(, t o>

k= p

(2.4.27) 6 X,t

(40)

= - w(r) R(r,t,t0),. s0

3(r

- r') = W

63(r

- r') op W t(r',t,t) -{0(r',t),H} op -(2.4.29)

with the matrix W equivalent to that obtained for the Lagrangian framework equation (2.3.14). The matrix equation for thus becomes

equal to:

(2.4.30)

The discussion given for the equation (2.3.13) for the Lagrangian frame-work also holds for the Hamiltonian frameframe-work. The rank of matrix W is

two, hence we can solve two of three velocity components, and its eigen

values and vectors are given by equation (2.3.15) and (2.3.16).

We now have to investigate the inner product of the null-eigen vector and the right-hand side of equation (2.4.30). Multiplying both sides of this equation with the null-eigen vector gives a new constraint

equa-tion:

a(r')R0(r',t,t0),. {0(r',t),H} =0 (2.4.31)

If this equation is satisfied for all time then the still unknown veloc-ity component can be chosen completely arbitrary. 1f this equation is not satisfied it produces a new constraint that reduces the hypersurface

given by O to a new hypersurface of lower dimensionality.

Using the GFD5, equations (2.4.22) through (2.4.25) and the definition

of a Poisson bracket we obtain:

w(r') R0(r',t,t0),. j0(r',t),H} = n = - a(r ) R°(r',t,t0),0 w(r') R(r',t,t0), s sijk g m ici kq o vt [t0,t1 (2.4.32) BR3

(41)

because the righthand side of equation (2.4.32) consists of an outer

product of m R with itself.

O n

Just as in the case of the Lagrangian framework the constraint equation (2.4.31) produces no new equation for and one of the components can be chosen completely arbitrarily. If it satisfies an initial condition it will satisfy the Bamiltonian equations (2.4.16) and (2.4.17) and its

constraint equation (2.4.28) for all time.

We now express two of the velocity components in the righthand side of

equation (2.4.30) and the third unknown velocity component:

(r',t,t0) w(r') R(r',t,t0),. {s2(',t),H} + w(r') R1(r'tt0) 3(r',t,t0)) (2.4.33) 2(r',t,t0) w(r') R(r',t,t0),J {1(r',t),H} + w(r') R2(r'tt0) 3(r',t,c0)) (2.4.34) with 3(r',t,t0) arbitrary.

If we insert the relations for and 2 in the Hamiltonian equations

for we obtain the remarkable result:

{R5(r',t,t0),H} S(rt ,t,t0),(r,t)} s r' ,t,t0),(r,t) np0 a() R3(r,t,t0),. w(r) R(r,t,t0),. w(r) R3(r,t,t0),. (2.4.35)

(42)

with arbitrary and matrix n° given by equation (2.4.5) with n = 1.

If we use the GFDs, equations (2.4.20) through (2.4.25) and the defini-tion of the Poisson bracket it immediately follows that we obtain the same equations for R' and R2 as observed in the Lagrangian framework, equation (2.3.17), while we obtain for the third component the trivial

result R3 R It must be remarked that the choice to express the

ve-locities R1 and

2

in terms of the velocity component R3 is arbitrary.

The same result is obtained when using or 2 as the unknown variable.

The Hamiltonian equations for fQJ) can be expressed in a generalized Poisson bracket or Dirac bracket equation:

{R5(r',t,t0),H}*

t, t0) (r, t)

with the Dirac bracket defined as:

f(r')g(r)}

w(r) RO(r,t,t0),.

k(rtt>

a()

Rk(r,t,t0),.

(2.4.37)), because it exactly has the same structure as a Dirac bracket defined in equation (2.4.8) and matrix n° is non-singular. The

equa-tions (2.4.36) and (2.4.37) must be supplemented with the constraint equation 0.

The final step in the proof of the theorem 2.4.1 will be the

demonstra-tion that we can incorporate the second contribudemonstra-tion on the right-hand

side of equation (2.4.36) in the Hamiltonian functional and we obtain an evolution equation fully expressed as a Dirac bracket relation.

(2.4.36)

_<f(r'),p(r,t)}

pG

{(rt)(r")>

O (2.4.37) w0(r) R (r,t,t0),.

(43)

This can be accomplished by using the following relation:

{

_a-

(r',t),(r,t)} a(r) R(r,t,t0),. O

Vt E [t0,t1} _(2.4.38)

which is a direct consequence of the fact that the null-eigen vector of

matrix is orthogonal to the right-hand side of equation (2.4.30).

Defining a new constraint 4(r,t):

$(r,t) = (rt) uJ(r) R(r,t,t0),. _(2.4.39)

and using the product rule for Poisson brackets, equation (2.4.7), _we

can transform equation (2.4.38) into:

0 _(2.4.40)

As a consequence of this relation the Dirac bracket of with any

arbi-trary functional reduces to a Poisson bracket. This relation can be used

to incorporate the second term on the right-hand side of equation (2.4.36) in the Dirac bracket by observing that it can be expressed as:

{R5(r',t,t0), <(r,t)

(r , t, to)

i

w0(r) Rk(r,t,t0),.

(k = 1, 2 or 3, no summation on k).

Here we used the fact that the integration in the inner product is per-formed with respect to the variable r, together with the product rule for Poisson brackets. The evolution equation (2.4.36) now can be fully

expressed as a Dirac bracket relation:

S(rtt)

_{{RS(rI,r,tO),IiT}*}

(2.4.42) (2.4.41)

(44)

if the Hamiltonian functional HT is defined as: ') R(r',t,t),. P(R, t> ' ,t) '1< R (r ,t,t0) W(r') Rk(rt,t,tO),i/0 (k = 1, 2 or 3, no summation on k).

In deriving the Dirac bracket relation we demonstrated that we obtained

the same equations for and R2 as derived in the Lagrangian framework,

equation (2.3.17) with completely arbitrary. If we use this freedom

and define in the same way as in the Lagrangian case, equation

(2.3.18), then equations (2.4.42) and (2.4.43) are equivalent with the

evolution equation for in the theorem 2.4.1.

It is possible to make some additional remarks about the result of theorem 2.4.1:

The Hamiltonian functional HT, defined in equation (2.4.43), has the remarkable feature that it contains a completely arbitrary component,

k _{(k = 1,} _{2 or 3). This means that the generalized} Canonical

trans-formation in the hypersurface U admits physically unobservable trans-formations, which was already mentioned in section 2.2, because the

evolution equation (2.2.20) is invariant under the transformation:

+ + Vi (2.4.44)

In section 2.3 we had to assume that the flow map R is a Canonical

transformation, thus with Jacobian determinant unity. Due to theorem

2.4.1 the evolution equation for the flow map, equation (2.2.20), can be expressed as a Dirac bracket relation, hence the flow map R can be generated by a generalized Canonical transformation using equations

(45)

(2.4.11) and (2.4.12). Sudarshan et al. [2-11, _PP. 122-130] proved

that on the hypersurface, defined by 4> 0, we can always find an ordinary Canonical transformation generated by a Hamiltonian for any

generalized Canonical transformation generated by a Hamiltoniam H so that the two transformations coincide in their effects on points of this hypersurface. Outside this hypersurface the transformations have

different effects, but this is of no physical importance.

The total flow map R now can be expressed as the following generalized

Canonical transformation:

= = exp (t LH) (2.4.45)

with:

LH = (2.4.46)

The Dirac bracket discussed in this section is invariant for the

fol-lowing Casimir function, viz.:

{(w0.V,R).(w0.v,R),H}* = {w.,H}* o

Vt E [t0,t1] (2.4.47)

with a divergence free vorticity field satisfying the boundary

con-dition equation (2.2.3).

The function w.w is called a Casimir function because the fact that the bracket is zero for all time solely depends on the structure of

the bracket, independent of the Hamiltonian functional. This means that the evolution of the vorticity field always occurs such that the

(46)

2.5. Definition and Lroof of converence of a 2roduct formula for the

solution of the viscous vorticit trans2ort equation

A numerical algorithm for the solution of the viscous vorticity trans-port equation which reflects the imtrans-portant physical phenomena discussed in section 2.1 is based on operator splitting, viz, the separation of convection and diffusion of vorticity in two steps. The solution of the viscous vorticity transport equation then is approximated by successive-ly solving the inviscid vorticity transport equation for a small time step, followed by diffusion of vorticity in the fluid and the creation of new vorticity at the body surface in order to maintain the no-slip

condition.

The application of this fractional step algorithm to the solution of the viscous vorticity transport equation was first proposed by Chorin [2-4].

This algorithm contained some heuristic assumptions and initially got a

lot of criticism, see for instance Milinazzo et al. [2-18]. Its validity

for the solution of the Stokes equations was, however, proved by Chorin et al. [2-19]. Special emphasis was put on the importance of the vortic-ity creation operator, viz, the creation of vorticvortic-ity at the body sur-face in order to maintain the no-slip condition. Without this operator the algorithm is not consistent for the Stokes and Navier-Stokes equa-tions. A complete proof of convergence of the fractional step algorithm, as proposed by Chorin, to the solution of the Navier-Stokes equations for an unbounded fluid was given by Beale et al. [2-20]. They showed that the accuracy of the original Chorin algorithm can be improved by

adding an additional diffusion step.

In this section a stochastic process is presented, defined as a product

formula, which approximates the solution of the viscous vorticity trans-port equation. The algorithm makes use of the solution of the inviscid vorticity transport equation as discussed in sections 2.2 through 2.4, together with a random walk description of the effect of viscosity. In the previous section it was shown that a unique solution of the inviscid

(47)

neighbourhood of the identity, by means of a Lie transform using a

suit-able Hamiltonjan functional and a Dirac bracket. This result will be used throughout this section by assuming that the inviscid flow map R is

known, therefore making the problem linear in the vorticity variable. The effect of viscosity now is to generate a random disturbance on the particle paths generated by the inviscid flow map R. The random walk interpretation of the incompressible Navier-Stokes equations, given by Peskin [2-21], therefore was very useful in the construction of the algorithm presented in this section.

A stochastic process simulating the evolution of a vorticity field in a viscous incompressible fluid can be defined by the following product

formula:

Definition 2.5.1

Let U R3 be the exterior unbounded domain of a C1 body surface S C 11

and UT be a small neighbourhood of width (12uT)½ on both sides of ,

with kinematic viscosity u and time step T.

Let w(X,T;B) :

R3 x [t0,t1] X S1 + L2(R3) be a stochastic vorticity field which is C in U U UT and zero outside U U UT. The variable

n E R3 represents the

ini-tially uniformly distributed fluid particle positions at time Tn= nT (n = 0, 1, ..), while the vector

n E SI, the unit sphere, is a random

vec-tor of unit length, independent from B for each n m, (n, m =

0,

1,

..)

and chosen from the uniform measure on the unit sphere.

The stochastic flow map Fm : R3 + R3 and the stochastic vorticity

evolu-tion operator L2(R3) + L2(R3) now can be defined as:

+ n =

nTn+TTn) +

(12 VT) n T ;B ) jn

(x

,T B ) = En Dn (X T B ) -

-n

n -n

T - -n n'-n T T - -n n'-n (2.5.1) (2.5.2)

(48)

with (X ,T +T,T ) the total inviscid flow map in the domain U U U as

--n n

n T

given by equation (2.4.45) and the identity mapping outside U U UT.

The operators , E and D, mappings from L2(R3) + L2(R3) represent the

vorticity boundary operator, the inviscid vorticity transport operator

and the vorticity diffusion operator respectively.

The vorticity boundary operator can be defined using the mapping .

Let : U1 + U1 be the mapping which reflects across the boundary S,

then the vorticity boundary operator maps the vorticity field in the

following manner: a(X ,T ;B ) (X ,T ;B ) = (X ,T ;B ) VX E U - -n n -n - -n n -n - -n n -n -n =

(X T B )

- -n n -n

=0

,VX

ER3-(UUU1)

-n

(2.5.3)

The inviscid vorticity transport operator is defined as the

nap-ping:

T ;3

)

E° c(X ,T ;B

) (x

T ;B ).V

R(X T +T T ₎

--n n-n

t--n n-n

--n n-n

X--n n

_-n

n

while the vorticity diffusion operator causes a random translation

of each point of the vorticity field:

T ;B ) D1' w(X ,T ;B ) = T (X ,T ;B )

- -n n -n

T -

-n n -n

(12vTYB -n n -n with Tqa () = + qa) the translation operator.

With the aid of the stochastic algorithm, given in definition 2.5.1, it

is now possible to state the main result of this section.

ex

E u

-

U

-n T

(2.5.4)

(49)

Theorem 2.5.1

Let the mean particle path r(t) E R3 and the mean vorticity field w(r,t) : R3 x [t0,t1] + L2(R3) _{in a viscous incompressible fluid be}

defined as the expectation of the stochastic process given in definition 2.5.1, viz.: E(Tn) + F r(T) = E[

X]

w(r T ) + w(r T ) = _E[K' _{(X ,T ;B}

)j

--n n

T--n n

T--n n-n

then the viscous flow map converges to the flow map Ft of the viscous flow and the viscous vorticity transport _operator _converges

to the evolution operator of the viscous vorticity transport equation _in

the L2(R3) norm 11.11 for n _{as long as the solution of the inviscid} flow map R(X,t,T), given by the Lie transform, equation (2.4.45), exists and is unique.

In order to prove theorem 2.5.1 we first assume _{convergence of product} formula F, equation (2.5.6), and use Lax equivalence _{theorem to prove} convergence of the product formula, equation (2.5.7), because this

map-ping is linear in the vorticity variable.

Lax equivalence theorem states that if for each T _{O the product}

for-mula K' Y Y is a bounded linear operator from the Banach _{space Y} into itself, with K _{equal to the identity mapping, then stability} _and consistency of K° is equivalent to convergence if the problem _{is linear.} For a general discussion of this theorem _{see for instance Chorin et al.}

[2-19] or Richtmyer et al. [2-22].

The first step in the proof of convergence of _{the viscous vorticity} transport operator K is to show = Identity. The product formula consists of three individual operators, viz. Dn _E _and _{and this}

con-T T

dition can be shown by looking at the individual _{operators for T}

= O:

(2.5.6)

(50)

D0 («X ,T ;B

_{O--n n-n}

) = T

O--n n-n

(«X T B ) = w(X ,T B ) - -n n -n E0 w(X ,T ;B ) = («X ,T ;B ).V R(X ,T ,T )

O--n n-n

--n n-n

X--n n

n -n = (X ,T B ) - -n n -n

using the relation: X = R(X ,T ,T ). -n

--n n

n

(2.5.8)

(2.5.9)

The condition = Identity further implies for the vorticity boundary

operator that we have to start with a vorticity field w as defined in

equation (2.5.3); thus with the correct viscous initial flow field. Com-bining the results of the various operators, using equation (2.5.2) and

taking the expectation of this equation immediately gives the

rela-tion = Identity.

The demonstration of consistency of the product formula K, equation

(2.5.7), is more laborious. The product formula is consistent if the following condition is satisfied:

w(r,t) = w(r,t).Vu(r,t) + vV2w(r,t)

=0

Vr C D U U - T Vr E - 1.1 - T (2.5.10)

for t = T with « the complementary volume of D. Here we use the

rela-t ions: T ;B ) D(r,T ;B ) E[ -] - E{ - _dt -n K w(r,t) 1T0 (2.5.11)

with 0/dt the material derivative. Before we can prove equation (2.5.10)

first the vorticity field and fluid particle position at time t = must be related to their values at time t = T0. This can be accomplished

(51)

X = R(X - (l2vT)½ _B _,T _,T

-n-1

-

-n

-n-1

n-1

n

,T ;B ) = ,T ;B ).V R(X ,T ,T )

-

-n-1

n-1. -n-1

- -n

n -a

X

_-n

- -n

n-1

n

The left-hand side of equation (2.5.10) can be evaluated up to order 0(t2) using a central difference approximation along the viscous flow

line F using equations (2.5.11) through (2.5.13); where we suppose that

both X and X are elements of 0:

-n-1 d n

=1

K w(r,T )

= -

_f

((X

,T1) - a(X ,T ))dS(b*)I4 + 0(T2)

T T

n 2T - -n+1 - -n-1 n-1 Ib =1

)Vn+i»

_n-1 E 1.1 (2.5.14) i

n}r* R{r

+ (12vT)

b*,

L*

11 T +T,T } - (r,T ).V R(r T -T,T ))ds(b*) n n

--

n

r-- n

n -0(T2)

Vr, E* E 0

_(2.5.15)

with dS(b*)141T the element of solid angle. The variables b* and r* are shorthand notations for b* = B r*

X +

(l2vT)½ b* and

=

This relation can be further evaluated using the following Taylor series expansions around the point r and time T up to order 0(T3):