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

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

A VARIATIONALLT OPTIMIZED VORTEX TRACING ALGORITHM FOR THREE-DIMENSIONAL FLOWS AROUND SOLID BOK

Proefschrift

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

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

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

te 14.00 uur

door

Jacobus Jozef Wilhelmus van der Vegt

geboren te Utrecht,

Dit proefschrift is goedgekeurd door de promotor:

STELLINGEN

Door het direct beschouwen van de deeltjesbanen in een niet-visceuze vloeistof 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 tijdsafhankelijke visceuze stromingen geven een zodanige overvloed aan informatie dat ze het inzicht in de relevante fysische mechanismen niet noodzakelijk versterken.

Stromingsvisualisatie-experimenten zijn essentiële hulpmiddelen bij de ontwikkeling van het fysisch inzicht noodzakelijk voor de ontwikkeling van mathematische modellen voor stromingsberekeningen.

Lagrange- en Hamilton-mechanica hebben een dermate hoog abstractieniveau bereikt dat het doel hiervan de gebruiker vaak ontgaat. Gezien hun grote voordelen verdient een uitbreiding van het toepassingsgebied van deze methoden de voorkeur.

De inspanning die verricht moet worden om modellen voor drie-dimensio­ nale instationaire visceuze stromingen te ontwikkelen staat niet in verhouding tot het 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 zijn.

Het verbeteren van het ontwerp van een scheepsschroef met verwaarlozing van de door het schip veroorzaakte vorticitelt 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 roeier een belangrijk deel van zijn voortstuwing tijdens een roeihaal bereikt door een liftkracht in plaats van een weerstandskracht op het blad is een sterk veronachtzaamd feit in de roeisport.

(V. Nolte, Die Effektivitat des Ruderschlages, Dissertation Sporthoch-schule Köln, 1984).

De interesse van de media in het liefdesleven van J.P. Sartre, na het verschijnen van de Nederlands vertaling van zijn biografie, past beter in een natuurfilm dan in een literaire boekbespreking.

CONTENTS Page

1. INTRODUCTION 1-1

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

3. NUMERICAL VORTEX TRACING ALGORITHMS 3-1

3.1. Introduction 3-1 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

4. 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 4-12

CONTENTS (continued) Page

5. RESULTS AND DISCUSSION 5-1 5.1. Introduction 5-1 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-11 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 6. CONCLUSIONS 6-1 7. SUMMARY 7-1 8. SUMMARY (DUTCH) 8-1 9. CURRICULUM VITAE 9-1 10. ACKNOWLEDGEMENT ... 10-1

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. From 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 experimental results is that it is very 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 phenomena for flows around bluff bodies, namely vortex shedding.

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 Lagrangian formulation 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 circular cylinder will be discussed.

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, result­ ing in much simpler equations.

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 Hamiltonian 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 dis­ cussed in section 2.3, while the derivation of a Hamiltonian functional

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.

2i2i_Basic_eg>uations_of_vortex_50tion

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 C vorticity field, u : Q x [t ,t ] + «, in the ex­

terior unbounded fluid domain Ü <= Rn, (n = 2 or 3) of a body, and let [tQ,t,] be a time interval in which such smooth solutions exist then the vorticity field satisfies the viscous vorticity transport equation:

■^ + u.Vu - u.Vu = W u> (2.2.1)

Here u = U+u : ft x [t0,ti] + fl is a smooth C velocity field whose curl is the vorticity field w. The velocity field must be divergence free by virtue of the continuity equation for incompressible flow. The velocity field u consists of two parts: a time dependent uniform velocity U and a divergence free disturbance velocity field u. The variable v in equation (2.2.1) represents the kinematic viscosity while 7 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 responsi­ ble for the diffusion of vorticity in the fluid.

The set of equations for the vorticity field must be supplemented with the following boundary conditions:

The velocity field u must satisfy the no-slip condition at a body surface S C3S1, with S a C manifold of dimension n = 1 or 2:

u = 0 (2.2.2)

This condition on the velocity field u means that the vorticity field w at S must satisfy, Serrin [2-6]:

n.u = 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 ft. Finally the velocity field u 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 u(r,t-) with associated vorticity field "(r.t.) = w (r) f o r e a c h r € £ï.

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 (2.2.1) by setting the kinematic viscosity v equal to zero. Further the

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

n.u = 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 well-known 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 no vorticity in the fluid at initial time tg 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 transformation 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 C invertible transformation:

R(r',t,t0) : QQ - U Vt e [ t0,t l] (2.2.5)

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

position r at time t are given by the relation:

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

r = 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(3R/3r') = 1 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 u the flow map R is separated into two

00

components: the flow map R which is equal to:

t

R"(r',t,t0) - r' + ƒ U(t')df (2.2.9)

and a disturbance flow map R. The total flow map R then can be repre­ sented by:

R ( r ' , t , t0) = R . R ( r \ t , t0) (2.2.10)

00

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:

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

00

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

R"(r',t,t0) = U(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 oig using a result obtained by Cauchy, see Serrin [2-6]:

"(r.t) = u>0(r') . 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 : R x [tQ.tjJ -»• Ü, then the relation between the Eulerian and Lagrangian description of the flow field, equation (2.2.11) can be expressed as:

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

= V_ x A(R(r',t,t ),t) (2.2.14) R

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 u = 0.

The introduction of a vector stream function A in equation (2.2.14) also gives a useful relation with the vorticity field a>. First, it must be noted that the vector stream function A is arbitrary up to the gradient of a scalar function i|> 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 A°° = U = R " (2.2.15)

together with a scalar field i|f such that A is divergence free by choos­ ing ty so that:

V2<|) = d i v ( A - A") ( 2 . 2 . 1 6 )

w h e r e :

A = A + A " + Vi|i ( 2 . 2 . 1 7 )

As a consequence of this separation both A and Vi|i 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:

a = V x V x A (2.2.18)

= -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:

R ( r \ t , t0) = V. x A(R(r',t,t0),t) (2.2.20) R V2A(r,t) = - u(r,t) (2.2.21) = - u)Q(r') . Vr, R(r,t,t0) (2.2.22) with: R ( r \ t , t0) = R°° . R(r',t,t0) (2.2.23)

The vector field A must be divergence free and satisfy the following boundary condition at S:

n.7 x A = -n.V x A™ (2.2.24)

together with:

V x A * 0 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 u0(r') which satisfies equation (2.2.3).

-i2i_裣i°2_S£i2£iEi£_£2I_£lïS_i2Yi2£

i

Ë_

vortic:

'-

t:

Z_I

ra

2522EÏ_Ê3U2£i2ïï

The inviscid vorticity transport equation possesses certain invariance properties, 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-8J. They reduced the number of variables by showing the close relationship of the Lagrange multi­ pliers with Clebsch potentials.

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 relation a = 'a x Vf5, where a and 3 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 a formulation of kinetic energy in terms of velocity and vorticity as can be found in Serrin [2-6].

Theorem 2.3.1

Let u e L2( Q Q ) be a C° initial vorticity field defined as the curl of the initial velocity field, with L (fig) the Hilbert space defined on the initial fluid volume fig a RJ at time tg, and satisfying the boundary condition 'n.u = 0 at a C manifold S c 3 0 , with unit normal vector n pointing into £2Q. Let R(r',t,t ) : 2 Q + f! be a C1 invertible flow map from the initial fluid volume fig into the fluid volume £1 at time t, hav­ ing a Jacobian determinant unity, with r' e UQ, which can be represented

by the relation R = R . R with R and R the flow map of the undisturbed and disturbed flow respectively. Further we define a C vector stream function A* : il x [tg.tjj ■»• L (ft) as the sum of a divergence free vector

stream function A and the gradient of a scalar function i|i> together with a Lagrangian functional L defined as:

L[R,A*,t] = - k< ^r').Vr,Ht',t,t0) x R(r',t,t0) | R ( r \ t , t0»

- <S0(r,).Vr,R(r',t,t0) | A*(R(r\t,t0),t)>

+ h < 7 x A*(r,t) | V x A*(r,t)> (2.3.1)

Then the following Action functional J defined as:

h

jfR.^.tjj.tJ = ƒ dt L[R,A*,t] (2.3.2)

t0

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

n.V x A* = -n.V x A™ (2.3.3>

with A the vector stream function of the undisturbed flow, has a sta­ tionary value in SÏQ for variations SR when satisfying equation (2.2.20) and a stationary value in fJ for variations 6A* with t.<SA* = 0 at S, t an arbitrary tangential vector at S, when satisfying equations (2.2.21) and (2.2.22). Here < f l g > represents the L2(ft) inner product defined as ƒ f.g dft, while a suffix zero refers to the inner product in the space L (S2Q) .

Before we can prove this statement we first have to define the gener­ alized functional derivative (GFD). Let F [ U ] = ƒ f(u,x)dx be a func­ tional of some element u € l/(fl) and 6u € C Q , infinitely differentiable variations of compact support, then the generalized functional

deriva-ÖF ^

tive — r , with u the i-th component of u, is defined as: 6u

/ OF I , i \ d „r 1 , . i-i

< ^ l

6 u

>

=

^

F [ u + £ 6 u ]

|

e

= o

The variations of a Lagrangian functional L with respect to the flow map R now can be defined as:

6L[R,R,t] 6L_ 6Rq 6 R l \ /O + < R constant in space /6L «R«\

/o

1 R constant 1 in space (2.3.5)

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: LjR.R.t] = i<o)J(r')Rp(r',t,t0),1 e ^ Rn(r' .t.t,,) | Rm( r ' , t , t0) > (2.3.6) and

^tS.S.t] " <

ü,

5<£

,

)8

ltpRP

(£

,

»

t

'

t

o>'i | ^ ( E d ' - t ^ o ^ ^ ^ (2.3.7)

w i t h g. . t h e m e t r i c t e n s o r , e^., = /g~ 6 ^ ^ t h e p e r m u t a t i o n symbol (£^23 = 1) and g = d e t ( gt 1) . Here we used t h e summation c o n v e n t i o n of t e n s o r a n a l y s i s . The g e n e r a l i z e d f u n c t i o n a l d e r i v a t i v e of h-^ w i t h r e s p e c t to R I s equal t o : \ 6 Rq 6 R l \

/o

R constant in space < ^ ; ( r ' ) R P ( r ' , t , t0) ,1 ^m qRm( r ' , t , t0) | 6 R ^ (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.

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

/Hi

6Rq < k »fc')*\r',t,t0),± e Rn(r'>t>t0) | «£<>> IR constant I in space 0 (2.3.9)

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

to

w

6R' 3A*k(R,t) 3A*k(R,t) oR

^

3?

6RP > 0 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:

<5L„

OR4

OR' <o0(r')R, i iik m iq*. 8kp 3A*P(R,t) 3RJ 6 Rq> R constant in space (2.3.11)

where we used the shorthand notation Rp for Rp(r',t,tn).

The generalized functional derivative of Lj with respect to SR is zero because Lo does not depend on R.

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

with respect to OR does not change when a gradient of a C 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

variations 6R. The first step in the calculation of the stationary value is to express the variations with respect to 6R in variations with respect to <$R. ƒ d t ƒ d 3 r. J l 6Rq = - k / d t ƒ d3r'|- (<o*(r')R?. e Rn) 6Rq t n 6Rq t n 3 t ° * p q n R constant in space

h

= - h ƒ dt ƒ d3r' üv;(r')RP e Rr ÓRq t Q ° * Pq r C0 "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 H Q for varia­ tions SR if the following relation holds:

^

,

>

R

V>

t

>

t

O>'i

e

pmq

i ,

V.t,t

0

) =

., 3A*P(R,t)

= u £ ( r ' ) R V , t , t0)) i n e e1^ gk p " Vr' £ QQ (2.3.13) oR

This equation represents a matrix equation for the unknown fluid par­ ticle velocity R with matrix W defined as:

V - - <"S<

r

*> ^ ( I ' . t . t Q ^

E p q m

Vr' £ n

0

(2.3.14)

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

If the matrix W was non-singular we could determine all velocity com­ ponents Rm 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 Rm can be determined. A detailed discussion of this process can be found

in Sudarshan [2-ll].

The rank of matrix W is two with eigen values:

X = 0 and X = ±\/-((«J R ^ )

2

+ <*>2 R ^ )

2

+ (<*J R?

±

)

2

) (2.3.15)

and the null-eigen vector X is equal to: XJ = u. R? (2.3.16)

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.

"1 "2

The velocity components R and R , for instance, can be expressed in the right-hand side of equation (2.3.13) and the arbitrary velocity

compo-' 3 nent R yielding: .1 lip 3A*q(R,t) ^ ( r ' ^ r ' . t . t ) R^r'.t,.^) = E l j P g = j —+- I 3 — (R <£'.t.t0> 0 p q 3RJ a.J(r')RJ(r',t,t0),. 3jp

aA*«(j,t)

-eg —i J PI 3RJ and

• 2 2 iP 8A*q(R,t) «*0(r')R2(r',t,t )

R

2

(r',t,t ) = e

2 j p

g - T —

+

"1 3 ~

R

V.t.t<)>

, . 3A*q(R,t)

" £ ^ ^ - } ( 2-3'1 7 )

where the expression E J g, —s-r- is equivalent with the i-th component

of curl A*. K

The stationary value of the Action Principle with respect to variations oR thus is obtained in ÜQ when the velocity components R and R satisfy

i o

the rx and r components of equation (2.2.20) together with a contribu-tion depending on the arbitrary component R and the rJ component of equation (2.2.20). If we now use the still existing freedom in the choice of RJ and set RJ at the initial time tQ equal to:

•3 3ip ^ ( i . t )

RJ(r',t,t) = e ^ g ^ (2.3.18)

u p q 3RJ

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

The second part of the statement viz., the fact that variations SA* 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 6R of the Action functional were invariant for the presence of V<|i 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 6A*. The separation of the vector stream function A* in two independent parts means that the Action functional must be minimized with respect to the following variations:

öA* = 6A + 6 Vip = ÓA + V6ip (2.3.19)

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

6A

* Jtl.A*.^,^] =

'l 3 = - ƒ dt ƒ dJrt(a;0(r,).VrlR(r',t,t0)).{6A(R(r',t>t0),t) 'o %

h

+ VJt(R(r\t,t0),t)} + ƒ dtj dJr(V x A(r,t)).(v x 6A(r,t)) (2.3.20) R

t

0

a

It is now necessary to transform the first integral in this expression from the initial volume Q$ to the actual volume ft 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 tQ can be related to each other:

ü)(r,t) = u(G(r,t,t1),t0) (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

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

det(^L)det(4) = d e t ( l 4 4 ) = l (2.3.24)

3r' 3r 3G 3R

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 Hamiltonian functional and vice versa presents no problems.

If the dynamical trajectory of the flow map R, is defined as the special mapping R so that the Action Principle is stationary for vari­ ations 6R then a function S J can be defined:

h

Sd (R(t0),t0,R(t1),t1) = ƒ dt L(R,R,t) (2.3.25) tQ along

-dyn

which acts as the generator type I of the Canonical transformation R which relates the flow field at time ti to the initial flow field, see Sudarshan et al. [2-11, 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 transformation is always unity, see Sudarshan et al. [2-11, pp. 73-77].

The Lagrangian functional in equation (2.3.1) 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:

h

6 A* J[ 5 .A* >t 0'tl ] = / d t / d r(-o)(r,t) + V x V x A(r,t)).«A(r,t) tQ a

h

ƒ dt ƒ d s(n x V x A(s,t) ).óA(s,t) fc0 S

h

- ƒ dt ƒ d s(n x V x A(s,t) ).<5A(s,t) (2.3.26)

where we used the boundary condition for the vorticity field, equation (2.2.3).

The integral along S,»,, the surface surrounding ft at infinity, gives no contribution if V x A + 0 at S„ faster than r , 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 ÓA in the surface S, while the integral is zero for any variations <SA 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 conditions to the Action functional which will be demonstrated in the next pages. The variations 6A E lr(ft) with t.SA = 0 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 ft:

72A(r,t) = - «(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

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

surface S <= 3fl:

Cl 'I 'l

ƒ dt L*[R,A*,t] = ƒ dt L[R,A*,t] + a ƒ dt ƒ d s n.Vx(A* + A ) fc0 " " C0 " " C0 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:

ƒ d2s n.Vx(A* + A*) = ƒ d3r V.Vx(A* + A°°) = 0 (2.3.29)

s ~ " "

a

oo 2

for any vector A*,A £ L (ft).

It is now possible to solve the vector stream function A in the domain ft without any reference to the function i|r while the variations of the Action Principle with respect to 6 R and SA* yield equations (2.2.20) through (2.2.22) with additional boundary conditions.

2^4^ Hamiltonian functional and Dirac bracket for the flow map 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.

2

A Lie algebra is a linear vector space, the Hilbert space L (£2) in our case, in which a Lie bracket is defined. This Lie bracket associates to

each element x,y € L2(Si) a third element z = [x,y] € L2(S2) and has the following properties:

a. The bracket must be antisymmetric

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

b. The bracket must be linear

[Xx + ux'.y] = X[x,y] + y[x',y] (2.4.2)

with X, u real numbers.

c. For any three elements x,y,z € L (ft) the Jacobi identity holds

[[x,y],z] + [[y,z],x] + [[z,x],y] = 0 (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 composition 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 one-to-one 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 associa­ tivity property of the Lie group in a coordinate independent way.

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:

(f(x).g(z)}

«g(Z)

(2.4.4)

with nM V the matrix:

,W 1 if li < ii, v = u + n -1 if v < n, y = v + n 0 otherwise (2.4.5) 6f 6e

with n = 1, 2 or 3 and -r—, -r3- generalized functional derivatives, de­ fined in equation (2.3.4), with respect to the phase space variables, the generalized coordinates and conjugate momenta z. The first three variables z are the generalized coordinates q , the rest the gener­ alized momentum variables p^. In addition to the conditions satisfied by the Lie bracket the Poisson brackets satisfy the fundamental Poisson bracket rule:

r p Hi p O

{z

w

,z } = n^

and the product rule:

{

f

l

f

2'

g

l

= f

lt

f2

*s} + {f!.8}f

2

(2.4.6)

(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 nu V which are now functions of the phase space variables z. The symbolic notation for these brackets is { , }*. The Dirac bracket defined as:

is a member of the class of generalized Poisson brackets with a singular matrix n,jV. In this relation the matrix Cm n is a non-singular matrix and the constraints $m define a subspace of the phase space to which the variables are confined.

When the equation of motion:

up = f(up,t) (2.4.9)

supplemented with an initial field un can be expressed as:

üP = {uP,H}* = Lu uP (2.4.10)

n

then this equation is called Hamiltonian and its solution can be written as a Lie transformation, an infinite sequence of brackets:

uP = exp (t l ^ u g =' E £ L £ UP (2.4.11) n>0

with:

l£ uP = L ^ " 1 uP) (2.4.12)

A detailed survey, especially concerning the group theoretical and to-pological aspects of Lie groups and algebras can be found in Sudarshan et al. [2-ll] 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 Hamiltonian 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

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, equation (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,t ) : B» + il be a C flow map, satisfying the evolution equation : R(r' , t, tQ) = VR x A*(R(r' , t, tQ) , t) with A* : S i x [tQ,t1] + L (K) a C°° vector stream function satisfying the condition n.7x(A* + A )

CO

= 0 at S c 3fi, with A the vector stream function of the undisturbed flow, then this evolution equation can be expressed by the Dirac bracket relation:

R - {R,HT}* (2.4.13)

with

%

= <»j

R

*

± 8kp

p

(

_R,t) >

o +

<$ l ^ r ) .

1 u0 'i

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

{f(r'),g(r")}* =

= {f(r'),g(r")} -<({f(r'),*m(r)}

~ Y n

mn

{*

n

(r),g(r")}>

(^(r) R ^ r . t , ^ ) ^ " X0 (2.4.14) (m, n = 1, 2) (j = 1, 2 or 3 )

Here <(>m and £ : S2Q x [tg.tjj ■* L2(«) define the hypersurface of the phase space to which the solution of equation (2.4.13) is restricted while the initial C°° vorticity field u> € L2(£2Q) 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 pk, defined as:

6L[R,R]

P

k

(£»t) =

6R

R

*<oJ(r') RP(r',t,t0),. e Rn(r',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 R is absent from equation (2.4.15). This presents problems when passing from the Lagrangian to the Hamil­ tonian formalism and vice versa. In the Lagrangian formalism we use the variables R and R, while we use R and p in t n e Hamiltonian formalism. It is no longer possible to express all the R variables in the £ variables when passing from the Lagrangian to the Hamiltonian framework.

Instead of a phase space with six independent field variables we now have only three independent variables R , 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 Hamiltonlan equa­ tions of motion, which for our field variables R and £ are equal to:

Rs(r',t,t0) - {R8(r\t,t0),H} +<^{RS(r',t>t0),4»p(r>t)}| Rp(r,t,t0))>

(2.4.16)

ps(r',t,t0) - {ps(r',t,t0)>H} +<lp8(r,,t,t0),4.p(r,t)}| Rp(r,t,t0))>

(2.4.17)

with H the Hamiltonian functional defined as:

»[?] = < J ? | 5 > - LtR. l ] ■ L?[g] (2.4.18) N ' 0

and L . [ 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:

y r ' . t ) = Pp( r \ t ) - k uj(r') RP(r' ,t,t()) ,± £p p n Rn(r',t,t0) - 0

(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 R 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 equation (2.4.16) to determine these velocities but using the following

Generalized Functional Derivatives (GFDs) it can be shown that this equation is trivially satisfied, viz.: R » R.

*R (r,t) 6Rk(x,t) 6Rk(x,t)X> = / 6s 6<3>(x - r) 0 X " 6Rk(x,t)" (2.4.20) '«RS(r,t) s6pk(x,t) ÓH 6Rk(x,t) Sp (x,t)> = 0 k / 0 6R"(x,t) (2.4.21) , m - * - i j P 3A*q(R,t) = <u>0<5> R (x,t,t0),m e m e gpq oR 6R"(x,t,toy>^ (2.4.22) 6H ^pk(x,t) 6p (x,t)> = 0 (2.4.23) «♦p(ï'.t) 6Rk(x,t) 6Rk(x,t) - < - % « ; ( ! ) RP(i,t,t0),. ep p k 6<3>(x-r') 6RkU,t,t0) (2.4.24) < ~ 6Pkl(ï»t)/> = < 6 6V ' ( x - r ' ) 6p (x,t)\ \«Pk(S,t) / O \ P k / l (2.4.25)

with 6* '( ) the three-dimensional Dirac distribution and & , the k Kronecker symbol.

The calculation of the GFD in equation (2.4.24) is non-trivial and will be further elucidated. The expression for the constraints tj>-, equation (2.4.19), can be transformed into:

The GFD in equation (2.4.24) now becomes equal to: 6* p ( l ' , t ) «Rk(x,t) «Rk(x,t); i

< ^ 7

r

T

(

"5

( £ , ) R P (

ï ' '

t

'

t

0

)

' i

E

p p n

R

V , t , t

0

) )

\ R ( x , t ) ^Rk(x,t) ( ^ ( r ' ) RP( r ' , t , t0) ^ R V . t . t , , ) , ^ Ö Rk( g , t , t0) ^ , k * p 5 Rk( x , t , t0) , k = p - i < w j ( x ) R ^ x . t . t Q ) ^ ep p k 6( 3 )( x - r ' ) 6 Rk( x , t , t0) > + *<u>£(x) Rn( i , t , t0) , . £ k p n6 (3) ( x - r ' )

- % <"JCs

) R ^ x . t . t ^ 6 p p k 6 <3 )( x - r«) 6 Rk( x , t , t0) 6R ( x , t , t0) (2.4.27)

So far the Hamiltonian 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 R.

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:

♦„(S'.O - {♦<,(£',t),H} + O*

0

(r

,

,t),<t.

p

(r,t)} R

RHP(r,t)> - 0

(r,t))

0

(2.4.28)

This relation can be converted into a matrix equation for the unknown variables RP using the GFDs for £, equations (2.4.24) and (2.4.25):

{*0Cr' ,t),<J>pCr,t)} = - coj(r) R ^ r . t , ^ ) ^ e 6( 3 )(r - r') pap

Wo p «( 3 )( l - I') (2.4.29)

with the matrix Wa p equivalent to that obtained for the Lagrangian framework equation (2.3.14). The matrix equation for R thus becomes equal to:

Wa p ^(r'.t.tg) - - { ^ ( r ' . O . H } . (2.4.30)

The discussion given for the equation (2.3.13) for the Lagrangian frame­ work also holds for the Hamiltonian framework. The rank of matrix Wa p 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:

wj(r') Ra(r',t,tQ),± {♦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. If this equation is not satisfied it produces a new constraint that reduces the hypersurface given by <j> ■" 0 to a new hypersurface of lower dimensionality.

Using the GFDs, equations (2.4.22) through (2.4.25) and the definition of a Poisson bracket we obtain:

wj(r') R^r'.t.tj,)^ {♦„(£',t),H} =

3A*q(R,t)

= - .

0

V ) RV.t,t

0

).

n

»5(r') RV.t,t

0

),

m

e

lot

e

1

* g

kq

^

because the right-hand side of equation (2.4.32) consists of an outer product of un R, with itself.

Just as in the case of the Lagrangian framework the constraint equation (2.4.31) produces no new equation for R and one of the components can be chosen completely arbitrarily. If it satisfies an initial condition it will satisfy the Hamiltonian 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 right-hand side of equation (2.4.30) and the third unknown velocity component:

R^r'.t.tQ) -7 T (-{4>2(r',t),H} + o)J(r') R^r'.t.tQ) R3(r',t,t0)) uj(r') R - V . t . t , , ) ^ U j U R2(r',t,t0) (2.4.33) o)J(r') R3( r ' , t , t0)> i with R (r',t,t0) arbitrary. { {<t'1(r,,t))H} + ujj(r') R2( r \ t , t0) , . . R3(r',t,t0)) (2.4.34) •1 '2

If we insert the relations for R and R in the Hamiltonian equations • P

for R we obtain the remarkable result:

RS( r \ t , t0) - {RS(r\t,t0),H} - O R8( r ' , t , t0) , tp( r , t ) } „P° uj(r) R3( r , t , t0) ,± {♦_(r,t),H} +/{RS(r',t>t0),*p(r,t)} with (p, o = 1, 2) Mj(r) RP( r , t , t0) , ■ 3 - : Ö K <.r,t,t0; lüj(r) R-'Qr.t.t,,)^ / 0 (2.4.35)

with R arbitrary and matrix np a 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

•1 «2

same equations for R and R as observed in the Lagrangian framework, equation (2 3 . 1 7 ) , while we obtain for the third component the trivial

•3 -3

result R ö R. It must be remarked that the choice to express the

ve-•1 '2 »3

locities R and R in terms of the velocity component R is arbitrary. •1 '2

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

The Hamlltonian equations for R can be expressed in a generalized Poisson bracket or Dirac bracket equation:

RS(r',t,t0) - {RS(r',t,t0),H}*

+ <{RS(r',t,t0),<!> (r,t)} -^ =- SL_lR*(r,t,t0)> (2.4.36) \ ° P »J(r) Rk(r,t,t0)>:l V 0

with the Dirac bracket defined as:

{f(r'),g(r")}* - {f(r,),g(r")}

\ P WJ(r) R ^ r . t . t , , ) ^ ° / O

(2.4.37)

(with k = 1, 2 or 3, no summation on k in equations (2.4.36) and (2.4.37)), because it exactly has the same structure as a Dirac bracket defined in equation (2.4.8) and matrix r\pa is non-singular. The equa­ tions (2.4.36) and (2.4.37) must be supplemented with the constraint equation <\> a 0.

The final step in the proof of the theorem 2.4.1 will be the demonstra­ tion that we can incorporate the second contribution on the right-hand side of equation (2.4.36) in the Hamlltonian functional and we obtain an evolution equation fully expressed as a Dirac bracket relation.

This can be accomplished by using the following relation:

Ua(r',t),$p(r,t)} wJ(r) A r . t . t ^ ) ^ - 0

Vt £ [tQ.tJ (2.4.38)

which is a direct consequence of the fact that the null-eigen vector of matrix Wö p is orthogonal to the right-hand side of equation (2.4.30). Defining a new constraint ij)(r,t):

i(r,t) = <t>p(r,t) w*(r) R ^ r . t . t g ) ^ (2.4.39)

and using the product rule for Poisson brackets, equation (2.4.7), we can transform equation (2.4.38) into:

{<t>a(r',t),Kr,t)} - 0 (2.4.40)

As a consequence of this relation the Dirac bracket of cj> with any arbi­ trary functional reduces to a Pols son 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:

/- Rk(r,t,t ) \

{RS(r\t,t ), (♦(r,.t) r - > }* (2.4.41)

\ u£(r) RK(r,t,t0),1/0

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

if the Hamiltonian functional H™ is defined as: HT "/<"£(£') R ^ r ' . t . t ^ ^ gj p A*P(R,tV> + <<|>(r\t) Rk(r',t,t0)

<oJ(r') R

k

(r

I

,t,t

0

),

1

/0

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

In deriving the Dirac bracket relation we demonstrated that we obtained •1 '2

the same equations for R and R as derived in the Lagrangian framework, .3

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

and define R 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 R 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 H™, defined in equation (2.4.43), has the remarkable feature that it contains a completely arbitrary component, R (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:

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

(2.4.11) and (2.4.12). Sudarshan et al. [2-11, pp. 122-130] proved that on the hypersurface, defined by <t> ■* 0, we can always find an ordinary Canonical transformation generated by a Hamiltonian H for any generalized Canonical transformation generated by a Hamiltonian 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:

R = R°° . exp (t LH ) R„ (2.4.45)

— - "•[ — U

with:

LH = { S 'H T} * (2.4.46)

• The Dirac bracket discussed in this section is invariant for the fol­ lowing Casimir function, viz.:

{(«S0-Vr,R).(uo.Vr,R),H}* = [u.üj.H}* - 0

Vt e [tQ.tJ (2.4.47)

with (D a divergence free vorticity field satisfying the boundary con­ dition equation (2.2.3).

The function u.u 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 vorticity is conserved along flow lines.

2i5i_Definition_and_2roof_of_convergence_of_a_2roduct_formula for the Ë°iü£i°.n_°f _ühe_Yiscous_ vor t ^ i H — E u S n u E S — — 2 H3.c i oB

A numerical algorithm for the solution of the viscous vorticity trans­ port equation which reflects the important 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 invlscid 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 vortie-ity creation operator, viz. the creation of vorticvortie-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 inviscld 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 vorticity transport equation can be generated, at least in a finite

neighbourhood of the identity, by means of a Lie transform using a suit­ able Hamiltonian 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 vortlclty 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 vortlclty field in a viscous incompressible fluid can be defined by the following product formula:

Definition 2.5.1

O 1 L e t ft c R be the exterior unbounded domain of a C body surface S c 3SÏ

and UT be a small neighbourhood of width (12vx)^ on both sides of S,

with kinematic viscosity v and time step T . Let w(X ,T ;B ) : RJ x

3 r - -n' n'-n

[tQ.tjJ x Sj^ -»■ L2(R3) be a stochastic vorticity field which is C°° in

a U U , and zero outside ft U UT. The variable X e R3 represents the ini­

tially uniformly distributed fluid particle positions at time Tn = nx (n

= 0, 1, . . ) , while the vector B £ Sj, the unit sphere, is a random vec­ tor of unit length, independent from B for each n 4 m, (n, m = 0, 1, ...) and chosen from the uniform measure on the unit sphere.

The stochastic flow map Fn : R + R and the stochastic vorticity

evolu-tion operator K : I/'(R-') + L/'CRr) now can be defined as:

X + P X = R(X ,T +T,T ) + (12 v x ) * B (2.5.1) -n T -n - -n n n -n

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

- -n n ' n v x

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

The operators *, E ° and Dn, mappings from L2( R3) + L2( R3) represent the vortlcity boundary operator, the inviscid vorticity transport operator and the vorticity diffusion operator respectively.

• The vorticity boundary operator $ can be defined using the mapping <(>. Let <|> : UT ■*■ UT be the mapping which reflects across the boundary S, then the vorticity boundary operator * maps the vorticity field in the following manner: u(X ,T ;B ) + * d)(X ,T ;B ) = w(X ,T ;B ) , V X £ a - -n n -n - -n' n -n - -n' n -n ' -n = «<4>X ,T ;B ) , V X € U - !! - -n' n -n ' -n x

, v x

n

e

R3

- (n u u

T

)

(2.5.3)

• The inviscid vorticity transport operator E is defined as the map­ ping:

u(X ,T ;B ) + En u(X ,T ;B ) = u(X ,T ;B )«VV R(X ,T +x,T ) (2.5.4) - -n n -n x -v- n ' n -n - -n n -n X -v- n ' n ' n

-n

• while the vorticity diffusion operator D causes a random translation of each point of the vorticity field:

u(X ,T ;B ) ■»• D" <ü(X ,T ;B ) = T , u(X ,T ;B ) (2.5.5) - -n' n -n x -s- n ' n -n ., „ ..?_ - - n ' n'-n'

v -n

with T u(r) = u(r + 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.

Theorem 2.5.1

3

Let the mean particle path r(t) £ R and the mean vorticity field o)(r,t) : R3 x [tQ.tjJ -»• L2(R3) in a viscous incompressible fluid be

defined as the expectation of the stochastic process given in definition 2.5.1, viz.:

r(T ) + Fn r(T ) = E [ F " X 1 (2.5.6)

- n x - n L x -nJ

w(r ,T ) * K" w(r ,T ) = E[K" W ( X ,T ;B )] (2.5.7) - -n n x -v- n ' n L x - -n' n'-n'J

then the viscous flow map F . converges to the flow map Ft of the

viscous flow and the viscous vorticity transport operator K , converges to the evolution operator of the viscous vorticity transport equation in

2 ^

the L (R ) norm 11*11 for n + » as long as the solution of the inviscld 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 x > 0 the product for­ mula K : Y + Y is a bounded linear operator from the Banach space Y into itself, with Kn 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 K. = Identity. The product formula K consists of three individual operators, viz. D , E and <b and this con­ dition can be shown by looking at the individual operators for x = 0:

D

o ï < * n » W = T

g

ïï(x

n

,T

n

;B

n

)

" < Xn, Tn; Bn) (2.5.8) EjJ U ( X ,T ;B ) = u(X ,T ;B ).VV R(X ,T ,T ) O - - n ' n -n' - -n' n'-n' X -x- n ' n' n' -n = u(X ,T ;B ) (2.5.9) - -n n -n v '

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

The condition K. = Identity further implies for the vorticity boundary operator * that we have to start with a vorticity field u 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 KQ = Identity.

The demonstration of consistency of the product formula K , equation (2.5.7), is more laborious. The product formula K is consistent if the following condition is satisfied:

^ Knu ( r , t ) . = iü(r,t).Vu(r,t) + vv2u(r,t) Vr 6 !i U UT (2.5.10)

T=0

= 0 Vr 6 a° - UT

for t = T with flc the complementary volume of ft. Here we use the rela­

tions:

3u(F r,T ;B ) Du(r,T ;B )

E [ "

C

a

t

" ~

n

] - E[ - V l - f p K ^ r . t ) . (2.5.11)

T=0

with D/dt the material derivative. Before we can prove equation (2.5.10) first the vorticity field and fluid particle position at time t = Tn_^

must be related to their values at time t = Tn < This can be accomplished