T H E APPLICATION OF VORTICITY METHODS T O ROTATING FLUID FLOWS By David G. McMillan B. Sc. (Applied Mathematics/Earth and Atmospheric Science) York University, 1993 A THESIS S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF M A S T E R OF SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES GEOPHYSICS A N D A S T R O N O M Y We accept this thesis as conforming to the required standard T H E UNIVERSITY O F BRITISH C O L U M B I A February 1996 © David G. McMillan, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Geophysics and Astronomy The University of British Columbia 129-2219 M a i n M a l l Vancouver, Canada V 6 T 1Z4 Date: FF6 . IS" \*)% Abstract Studies of E a r t h ' s core dynamics often require a computational method which can account for non-linear effects and non-periodic time dependence. One important application involves the large scale, possibly chaotic, fluid motion associated w i t h a past resonant tidal forcing of the free core nutation. This thesis explores the utility of vorticity methods for modelling contained rotating fluids, such as the E a r t h ' s core. Establishing this method i n the general context of rotating fluids is the first step i n the development of a non-linear, three dimensional, time dependent model of Earth's fluid core. Conventionally, vorticity methods have been applied to incompressible fluid flow i n infinite domains w i t h small, finite regions of non-zero vorticity. T h e capability of ap- proximating solutions to non-linear, three dimensional, time dependent, inviscid fluid flow problems suggests that these algorithms may also be well suited to flows i n finite domains under uniform rotation. Physically, the method under consideration involves grid-free tracking of fluid particles which carry distributions, or 'blobs,' of vorticity, each w i t h a prescribed strength. T h e strengths and positions of the blobs determine the velocity field through the Biot-Savart law. T h e fluid particles, along w i t h their vorticities, are advected by the velocity field. In addition, the vorticity of each blob is altered by the strain associated w i t h the velocity field and must be updated at each time step. T h e method of solution reduces to a set of first order, ordinary differential equations i n time for the Lagrangian displacement of the fluid particles and their corresponding vorticity strengths. T h e O D E s are advanced by standard integration routines. This work is a presentation of the development of the vorticity method algorithm ii for rotating fluids. In addition to the development of the algorithm, several new computational methods are described. T h e concept of effective vorticity is introduced as a means of relating known physical quantities to computational results. In order to satisfy inviscid boundary conditions i n an efficient manner, approximate methods using image particles are developed and incorporated into two of the test cases. Quick initialization of vorticity fields is made possible by the formulation of an approximate linear system of equations. T h e modelling of a uniformly rotating fluid requires the initialization of a uniform vorticity field. This procedure is decomposed into two steps, increasing efficiency while creating a reasonable initial vorticity field. T h e implementation of the algorithm is verified by comparison w i t h theoretical approximations. T h e first test cases involve the self advection of a t h i n vortex ring and the interaction of two vortex rings. These examples confirm the operation of the algorithm without introducing the issue of boundary conditions. Solid boundaries are introduced i n the problem of standing waves on a bounded vortex filament. Finally, the problem of a contained fluid w i t h a uniform vorticity distribution is examined by modelling inertial modes i n a rotating cylinder. iii Table of Contents Abstract ii List of Tables vii List of Figures viii Acknowledgement xi Dedication xii 1 INTRODUCTION 1 1.1 Background 1 1.2 Application of Vorticity Methods to Geophysical Fluid Dynamics 1.3 Scope of The Thesis 3 2 THE VORTICITY METHOD 7 . . . . 2 2.1 Analytical Formulation 7 2.2 Temporal Discretization 12 2.3 Spatial Discretization 12 2.4 Development of the Algorithm 18 2.4.1 Initialization 18 2.4.2 Evaluation of Derivatives 19 2.5 Accuracy and Convergence 21 2.6 Boundary Conditions 23 iv 2.6.1 Inviscid Fluid Flow 23 2.6.2 Method of Images . . 24 3 VORTEX FILAMENTS: BACKGROUND 30 3.1 Physical Properties 30 3.2 Discretized Vortex Lines and Filaments 35 3.3 Computing the Circulation of a Discretized Filament 39 4 VORTEX FILAMENTS: NUMERICAL MODELS 41 4.1 Self Advection of a Vortex Ring 41 4.2 Interaction of Two Vortex Rings 47 4.3 Waves on a Bounded Vortex Filament 50 5 ROTATING CYLINDER 63 5.1 Introduction 63 5.2 Inertial Waves 64 5.3 Initialization 67 5.4 Excitation of Inertial Waves 72 5.5 Numerical Modelling 75 5.5.1 Estimating Frequencies of Oscillation 76 5.5.2 Results 78 5.5.3 The Vorticity Profile and Stability of Radial Perturbations . . . . 84 6 DISCUSSION AND CONCLUSIONS 87 6.1 The Vorticity Method 87 6.2 Initialization 88 6.3 Filament Results 90 6.4 Cylinder Results 91 v A 6.5 Boundary Conditions 92 6.6 Resolution 94 6.7 Summary 95 T H ESHEET VORTEX METHOD REFERENCES 97 103 vi List of Tables 5.1 Frequency estimates for three axisymmetric modes. Estimates are given for P = 1,2,3,4 83 vii List of Figures 2.1 Distribution function ip . This version was evaluated with 6 — 0.5 2.2 Geometry of a single blob and planar boundary. The normal component s of velocity is evaluated at three points A, B and C along the boundary. . 2.3 27 Normal velocities at planar boundary. Three velocity profiles evaluated at A, B and C in Figure 2.2 3.1 14 29 Bending and stretching of a vortex line. The vortex line is deformed by a differential velocity between the points P and Q. Adapted from Batchelor [1967] 3.2 31 Illustration of a vortex tube. The cross sectional areas A are enclosed by i the material curves 7; 3.3 32 Illustration of a discretized vortex line. The arrows are tangential to the vortex line at the particles' centers and represent the direction of each particle's vorticity 3.4 36 Azimuthal velocity structures of straight filaments. The dottedfinecorresponds to the velocity in (3.6), with Q = ^u>. The solid line corresponds to a discretized filament 3.5 37 Determination of the effective radius. For each value of S, the peak in azimuthal velocity is located at approximately |t5 4.1 38 Illustration of a ring vortex. The arrows at the sides indicate the sense of fluid rotation in the ring's core 42 viii 4.2 Convergence of computed ring velocities. The dotted lines give the upper and lower bounds on the theoretical estimate 4.3 Computing the circulation of a ring vortex. The curve 7 is shown with a cut-away view of the ring 4.4 45 A comparison of theoretical and computed velocities. The dotted lines give the upper and lower bounds on the function U(A) given by (4.2). . . 4.5 46 Relative error in ring velocities. Errors were determined using the theoretical approximation (4.2) with constant circulation, C = 0.2907 4.6 44 47 Evolution of two ring vortices. The left and right panels show the side and plan views of ring positions, respectively. Top to bottom, t — 1.00, 4.11, 5.46 and 7.44 4.7 48 Illustration of a bounded filament. The solid arrow indicates the sense of the fluid rotation in thefilament'score. The dashed arrow indicates the direction of the precessional motion 4.8 Convergence of computed filament frequencies, mode 1. The dotted lines give the upper and lower bounds on the theoretical estimate 4.9 51 57 Convergence of computed filament frequencies, mode 2. The dotted lines give the upper and lower bounds on the theoretical estimate 58 4.10 Convergence of computedfilamentfrequencies, mode 3. The dotted lines give the upper and lower bounds on the theoretical estimate 59 4.11 A comparison of theoretical and computed frequencies. The theoretical values include the order of magnitude error estimates. Increasing values of A correspond to increasing mode numbers 60 4.12 Relative error in filamentfrequencies.Errors were determined using theoretical approximation with C = 0.2935 ix 62 5.1 Plan view of filament arrangement. 81 filaments arranged in five shells with radii given by (5.8) 5.2 69 Vorticity and azimuthal velocity profiles. Profiles correspond to the five shell case, with modified r and r (see § 5.5). The circles give the radial 4 5 locations of the shells 5.3 71 Oscillations in positions and strengths for mode (1,1,0). From top to bottom, radial perturbation, axial position and vertical component of strength, i indicates shell. . 5.4 79 Oscillations in positions and strengths for mode (2,1,0). From top to bottom, radial perturbation, axial position and vertical component of strength, i indicates shell 5.5 80 Oscillations in positions and strengths for mode (3,1,0). From top to bottom, radial perturbation, axial position and vertical component of strength, i indicates shell 81 x Acknowledgement Thanks to all of my colleagues at both UBC and York for their support and encouragement. In particular, my gratitude goes to Bruce Buffett for being most generous with his time and ideas. Also, thanks to Matt Yedlin and William Hsieh for their enthusiastic participation as members of my committee. Primary support for this work was provided by NSERC in the form of a PGS-A scholarship. xi memory of my grandparents, Ida and Murray Burgess. xii Chapter 1 INTRODUCTION 1.1 Background Vorticity is defined mathematically as the curl of a velocity field and hence is a vector quantity defined at every point in the fluid. It describes rotational qualities of the fluid and where the vorticity is zero the fluid is referred to as irrotational. The dynamics of vorticity have long been recognized as a very important aspect of fluid flow studies. Contributions by Kelvin and Helmholtz in the nineteenth century— which are still regularly cited—form the basis of much of this research. The combination of well established ideas and the advance of computational methods has led to efficient and elegant fluid dynamical algorithms called vorticity methods, or vortex methods. Vorticity methods are a class of grid-free computational algorithms which exploit favourable properties of vorticity and its relationship to the velocity field. They were developed primarily as a means of understanding the onset and character of turbulence. As such, these methods fully incorporate non-linear effects in a three dimensional domain. Traditionally, vorticity methods have been applied toflowsin infinite domains where the vorticity is non-zero only in a small finite region of the fluid, such as near the boundary of an air foil, or in the trailing vortices of a bluff-body. Since the velocity field of an incompressible fluid can be determined from the vorticity field, it is sufficient to follow only those fluid elements with non-zero vorticity. Inflowswith small patches of vorticity, such a convenience greatly reduces the number of required computational elements. 1 Chapter 1. INTRODUCTION 2 In geophysical fluid dynamics there is a continuum of vorticity due to the persistent rotation of our planet. Thus, applying vorticity methods to rotating fluids does not take advantage of the computational efficiencies that arise when vorticity is sparsely distributed. On the other hand, the exploitation of vorticity dynamics provides new insights into the dynamics of rotating fluids, in part, because an inertial rather than a rotating reference frame is adopted. Moreover, the capability of approximating three dimensional, non-linear, time dependent solutions makes these algorithms well suited to many problems involving rotating fluids. This work is an exploration of the utility of a particular vorticity method for rotating fluids in general, and specifically for core dynamics. 1.2 Application of Vorticity Methods to Geophysical Fluid Dynamics Vorticity dynamics has found a place in oceanographic and atmospheric sciences. In both thin layer approximation holds and a quantity called the potential vorticity is conserved. In regions such as Earth's fluid core, there is no analogous quantity and a more generalized approach is required. The development of the vorticity method in this thesis is directed toward contained, rotatingfluidswith typical geometries including cylinders, spheres and spherical shells. In geodynamics, the conventional approach to such problems is to introduce simplifying assumptions, such as periodic time dependence or axisymmetric flow [Smylie, et al, 1992]. However, there are a number of problems that do not permit such simple representations (see, for example, Malkus [1989]). One specific example involves the free core nutation (FCN), which is a normal mode of oscillation in the fluid core. It has a natural frequency which depends on the ellipticities of the fluid core's boundaries and Earth's rotation rate. Since the rotation rate is decreasing due to tidal braking, the ellipticities of the fluid core boundaries, and hence, the natural Chapter 1. INTRODUCTION 3 frequency of the FCN are decreasing. Further, the frequency of diurnal solar tidal forcing changes with Earth's rotation rate. A comparison of the two frequencies indicates that, in the past, there was a resonance of the FCN with solar tidal forcing. The large scale, possibly chaotic, fluid motions due to this type of interaction must be modelled by a computational method which extends easily into the non-linear regime. Further, large or small scale non-periodicflowmust be modelled without explicit assumptions imposed on its time dependence. Lastly, in some problems, it is desirable to have a solution in three dimensions, as opposed to a solution in two dimensions which has a prescribed azimuthal dependence. For example, in geodynamics there has been interest in precession of Earth's inner core caused by pressure torques at the boundary between the solid and fluid cores [Mathews, et al, 1991]. Modelling of such pressure fields using the vorticity method would permit accurate predictions of boundary pressures due to non-axisymmetric flow in the fluid core. 1.3 Scope of T h e Thesis The goals of this project are: • to develop a computational algorithm based on the dynamics of fluid vorticity, • to perform preliminary tests by modelling simple vorticalflowsand comparing the computed solutions against known solutions, • to adapt the computational algorithm to contained, rotating fluids. In addition to verifying the implementation of the algorithm, additional details regarding the method's application to rotatingfluidsare explored. Some of the novel aspects of this work include: relating known physical quantities to computational parameters; applying Chapter 1. INTRODUCTION 4 conditions at planar boundaries in an efficient manner; and reproducing uniform vorticity fields in contained fluids. The vorticity method involves the tracking of fluid particles whose vorticity is concentrated around a point. However, the vorticity distribution varies in such a way that it is non-zero everywhere in the fluid. This leads to vortex filaments with non-uniform vorticity in their cores. A simple physical argument shows that a quantity called the effective vorticity of a computational filament is representative of the uniform vorticity of an equivalent theoretical filament. When incorporating planar boundaries into a model, the well known method of images is shown to apply to the vorticity vectors associated with fluid particles. The variation of a particle's contribution to the velocity at a planar boundary, as its distance increases, leads to an approximate method of images as a means of satisfying the boundary conditions. The image particles outside the domain of interest are found to be equally important for maintaining a uniform vorticity within the fluid and especially near boundaries. Since additional fluid particles are required outside planar boundaries, the initialization of a specified vorticity field is handled more efficiently by solving an approximated system of equations. In contained rotating fluids the number of computational elements increases because vorticity is distributed uniformly throughout the fluid. In this case, a two-step initialization method is invoked, further reducing the computational burden of this important process. The primary focus of the implementation of the vorticity method will be on reproducing general behaviour of rotating fluid flows, while keeping a watchful eye on efficiency. In addition, the process of achieving the above goals and developing new ideas will give many insights into the intricacies of the vorticity method. This thesis is organized as follows. Chapter 2 begins with an analytical formulation of the algorithm and corresponding boundary conditions. Some attention will be focused on a physical interpretation of the more important aspects of the method of solution. Then, Chapter 1. INTROD UCTION 5 the discretization of the analytical formulation is described. There are two different approaches to this step—both yielding the same result—which give additional insight into the physical interpretation of the vorticity method. Included is a brief description of the implementation of the discretized algorithm. Finally, the approximate method of images, as it applies to the discretized algorithm, is presented. The preliminary testing of the vorticity method algorithm involves flows generated by long, thin regions of vorticity called vortexfilaments.In order to put the test cases in a proper perspective, Chapter 3 consists of a discussion of the basic physics of vortex filaments in inviscid fluids. This chapter is essentially a qualitative description of many fundamental and important results of vortex dynamics and includes details of the discretized versions of these results. Chapter 4 details preliminary testing of the vorticity method algorithm. The tests are accomplished by modelling three cases of relatively simple vortexfilamentsand their self-induced motion. Thefirstqualitative test involves a comparison of the self advection velocities of vortex rings with those predicted by theory. An interesting example is the interaction of two vortex rings, which illustrates the power of the vorticity method when applied to more complicated flows. A second qualitative test incorporates two parallel planar boundaries and the method of images. Here, standing waves on a filament extending between the boundaries are modelled and the frequencies are compared with theory. The presentation of each case outlines specific difficulties and surprises encountered in the modelling process. The extension of the vorticity method algorithm to contained, rotating fluids is presented in Chapter 5. The goal is to model inertial waves in a cylindrical container and compare computed results with those predicted by theory. Inertial waves are introduced from both mathematical and physical perspectives. The initialization of a uniformly rotating fluid within a cylinder and the subsequent perturbation of single inertial modes are Chapter 1. INTRODUCTION 6 then discussed, followed by a presentation of the computational results. An important issue, giving rise to some of the modelling methods in this chapter, is the instability of circular pathline flow and is included as a possible explanation of fluid particle behaviour observed in the results. In addition to a discussion of results, Chapter 6 provides a summary of, and possible improvements to, the supplemental processes which must be included in rotating fluid models. This includes issues such as initialization, treatment of boundary conditions and spatial resolution. In general the results are promising. Accurate modelling of contained, rotating fluids will require further attention focused on the reproduction of a uniform vorticity field and a more general method for imposing boundary conditions. Chapter 2 T H E VORTICITY 2.1 METHOD Analytical Formulation If a fluid velocity field is given by u, then its vorticity u> is denned as w = V x « . (2.1) The vector u> is sometimes called the absolute vorticity when it and its spatial derivatives are referred to an inertial frame. In the analysis that follows all relations are derived with respect to an inertial reference frame. The vorticity is denned at every point in the fluid and describes the rotational qualities of the flow field. There are many ways of drawing analogies between the vorticity and local fluid motion. By examining the relative motion of two initially perpendicular line elements, it can be shown that the average angular velocity of two such lines is proportional to the component of vorticity perpendicular to the plane containing the lines. Alternatively, the vorticity at a point is proportional to the instantaneous angular momentum of a spherically symmetric fluid particle centered there. One of the primary reasons the vorticity method is attractive for modelling rotating fluids is that rotation and vorticity are intimately related. A fluid which is rotating as a solid body has throughout it a constant vorticity. Where 17 is the rotation rate of the fluid, the vorticity is u = 212. A vortex line is one whose tangent is everywhere parallel to the local vorticity vector. Chapter 2. THE VORTICITY 8 METHOD Vortex lines in a uniformly rotating fluid are straight and parallel to the rotation vector Q. In an iriviscid fluid, it may be shown that a vortex line deforms as a material line of fluid [Batchelor, 1967]. This means that a vortex line always consists of the same fluid particles. This important result will be examined in more detail in § 3.1. The flow of Newtonian fluids is governed by the Navier-Stokes equation. In the absence of body forces, this equation represents the conservation of momentum and may be written as * 5 _ _ I v , + „*•«. (2.2) where u is the fluid velocity, p the density, p the pressure, v the kinematic viscosity and ^ =l + « - v <-> 2 3 is the material derivative. Inflowswith large Reynolds numbers, viscous effects are very small compared to inertial effects throughout most of the fluid. A useful approximation to the flow is obtained by omitting the viscous effects, which reduces the Navier-Stokes equation to Euler's equation, Euler's equation, together with the conservation of mass, or incompressibility condition, V-u = 0 (2.5) complete the description of the inviscid fluid flow problem. Equations 2.4 and 2.5 are subject to boundary conditions which require the continuity of the normal component of velocity at all impermeable boundaries. Taking the curl of Euler's equation results in the vorticity equation, Chapter 2. THE VORTICITY METHOD Duf Dt 9 = (u> • V ) u . (2.6) The vorticity equation is useful in two respects: taking the curl has eliminated the pressure variable; and it is the first step in casting the flow problem in terms of the vorticity. Equation 2.6 describes the evolution of vorticity following the fluid particles. This description is usually referred to as a Lagrangian one. Along particle trajectories, the vorticity is changed by the strain associated with the velocity field. In terms of vortex lines, the RHS term represents bending and stretching of a vortex line. There are two interesting notes to make regarding different forms of the vorticity equation. Firstly, in two dimensional flow the straining term vanishes and, consequently, vorticity is conserved while following the fluid. In other words, the vorticity of a fluid particle in an inviscid flow is constant in time. Secondly, the vorticity equation derived from (2.2) for viscous flow has an additional term on the RHS, namely i/V u>. This term 2 arises from applying the curl operator to the second term on the RHS of the Navier-Stokes equation, (2.2). Just as the viscosity term in (2.2) represents a diffusion of momentum, the analogous term in the vorticity equation represents a diffusion of vorticity. Hence, the vorticity of a particle in a three dimensional viscous flow is changed by both strain and viscous diffusion. The strategy for solution by the vorticity method begins with the vorticity equation. Given an initial distribution of vorticity, its time derivative due to the associated velocity field is calculated according to (2.6). Since the evolution of vorticity depends explicitly on w, the first step is to estimate u from the prescribed vorticity field. Given the two relations, (2.1) and (2.5), a convenient approach is to define a vector potential or stream function <& such that, Chapter 2. THE VORTICITY METHOD 10 u = Vx#. (2.7) Poisson's equation for the vector potential is formed by taking the curl of (2.7). Expanding the RHS of (2.7) using vector identities and assuming <& can be chosen such that its divergence is zero, the result is V * = 2 -u) . (2.8) The solution of (2.8) is given by the Biot-Savart Law, *•(*,-) = j G(x - a?') u>(x',t) dx', (2.9) where G(x) = l/47r|a?| is the fundamental solution, or Green's function, in three dimensions and u is the vorticity at the source point x'. Applying the curl operator to (2.9) yields a solution for the velocity in terms of the vorticity, u (as, t) = J V x G (as - x') u> {x\ t) dx'. (2.10) Thus, in principle, from an initial vorticity field, the associated velocity field can be found. Using the velocity field, the time derivative of vorticity is calculated and the vorticity is advanced in time. The process can then be repeated with the new vorticity field. This is the essence of the vorticity method. The solution (2.10) assumes that V • $ = 0 everywhere in the fluid domain. This assumption may by tested by taking the divergence of (2.9) (see Batchelor [1967]). The result yields an additional constraint on the vorticity: everywhere on the fluid's boundaries, the normal component must be zero, ie. u> • n = 0. Saffman [1992] obtains the same constraint by verifying that the curl of (2.10) reproduces the specified vorticity field OJ. This condition is satisfied for fluids in infinite domains provided the fluid velocity Chapter 2. THE VORTICITY METHOD 11 goes to zero in the far field [Batchelor, 1967]. However, when the condition u • n = 0 is not satisfied in afiniteregion, it is still possible to apply (2.10) provided the fluid region is extended in such a way that the normal component of vorticity vanishes at the new boundaries. This is analogous to closing vortex lines onto themselves outside the original fluid volume, or making them infinitely long. The velocity field is now found by taking the integral of (2.10) over the extended region of the fluid. Theflowinduced within the original boundaries by the vorticity outside these boundaries represents a potential flow. In the Lagrangian description offluidflow,the observer follows fluid particles as they move with the flow. Physical quantities associated with a certain particle are functions of its initial position a and time t. The trajectory, or path, of a fluid particle is the solution to the ordinary differential equation -dJ = «(•'*)' (2.11) where u is the particle's velocity and the initial particle position is given by SB (at, 0) — a . (2.12) Along the particle trajectories, the vorticity equation takes the form of an ordinary differential equation, (2.13) w(a,0) = w . o (2.14) Using the Biot-Savart expression to evaluate the velocity u in (2.11) and (2.13), results in a set of coupled, ordinary differential equations in time for the Lagrangian positions Chapter 2. THE VORTICITY METHOD 12 of fluid particles x and their corresponding vorticities u>. The three equations, (2.11), (2.13) and (2.10), and the initial data given by (2.12) and (2.14), represent the analytical form of the vorticity method algorithm. 2.2 Temporal Discretization The time stepping aspect of this algorithm follows naturally from the above analysis. The analytical form of the vorticity method was pieced together by deriving expressions for the time derivatives of particle positions and their vorticities. If an initial distribution of vorticity is specified, the associated velocity field can be evaluated by the Biot-Savart law. The calculations of the position and vorticity derivatives then follow. The particle positions and vorticities are subsequently advanced according to an integration scheme such as the Runge-Kutta method. The process can be cycled using the updated quantities. In fact, this algorithm is very well suited to time stepping methods like Runge-Kutta since the ODEs are cast in an appropriate form. Further details about the evaluation of the position and vorticity derivatives are given in § 2.4. For the remainder of this chapter there will be no explicit reference to a particular time step or temporal discretization. 2.3 Spatial Discretization Let the fluid region of interest be discretized into N parcels each having a volume dV , x i = 1,..., N. The velocity and vorticityfieldsare then approximated by assigning values, u(se') and (*>(x ) respectively, to each fluid particle. The positions as* of each parcel can % be taken at their centers and the assigned velocities and vorticities are average values taken over the volumes dV . For brevity these quantities will often be written u* and l w*. After discretizing all regions of the fluid with non-zero vorticity and establishing the fluid particles and their velocities and vorticities, the ODEs (2.11) and (2.13) become Chapter 2. THE VORTICITY 13 METHOD ^ = «(•*,*). (2-15) = (u> .V) u (»*,*) • (2.16) i eft In order to evaluate these derivatives, an analogous, discretized version of the BiotSavart integral for u* must be found. An examination of (2.10) shows that the velocity becomes unbounded as x' —• x. In the discrete version this occurs when the field and source points coincide. Since this is physically unrealistic, the problem is corrected by replacing the fundamental solution G with a function that closely resembles it, except at the origin where the kernel is bounded. The new function Gs can be found by convolving G with an appropriate smoothing function [Anderson and Greengard, 1985]. The integral for the velocity is then approximated by the summation, N u ( « \ t) = £ V x G (se* - x ) or* (t) dV . j (2.17) j 6 The subscript 8 indicates that the fundamental solution G has been convolved with, or smoothed by, a distribution function ij> , defined below. The end result is that a fluid s particle's vorticity makes a finite, non-zero contribution to the velocity it experiences. The gradient of the velocity in (2.16) can be taken either numerically or analytically. There are advantages and disadvantages to each of these options. The primary disadvantage to the numerical differentiation is that it creates an additional source of numerical error. As well, there is the necessity of carrying extra particles outside the domain of interest in order to compute the derivatives near the boundaries. On the other hand, the analytical approach imposes constraints on the distribution function in order to maintain a desired level of accuracy. An appropriate selection for il> , as proposed by s Chapter 2. THE VORTICITY g 14 METHOD 2.0 Radius r = U l Figure 2.1: Distribution function ij> . This version was evaluated with S = 0.5. s Beak and Majda [1985], is shown in Figure 2.1. The mathematical form of this function is, (2.18) where r = |se| and i>(r) = -fe- 3 (5-3r ) . 3 (2.19) Chapter 2. THE VORTICITY METHOD 15 The amplitude of the distribution function scales as S~ , and hence approaches a Dirac 3 delta function for small values of the parameter S. A more detailed discussion on the selection and properties of this function is presented in § 2.5. By defining a new quantity, called the vorticity strength, T* = u>*dV , (2.17) can be l rewritten as N u (x\t) = ^ V x G (**' - x ) F> (t) . j { (2.20) i=i It is convenient to work with the particle strengths, as opposed to their vorticities, for several reasons. The most immediate advantage is that there is no need to evaluate and store both the volume elements dV and the vorticity u>*. If the initial positions of fluid particles are the nodes of a rectangular grid all the volume elements would be equal, and straight forward to compute. But since the vorticity method is grid-free, the initial arrangement offluidparticles need not be confined to a rectangular grid. In fact, in some scenarios it would be foolish not to take advantage of the grid-free aspect of the algorithm. In such a case, with all the volume elements being different, their determination can be cumbersome. For an incompressible fluid, the volume elements dV are invariant in time, which 1 makes it possible to recast the ODEs in terms of the vector strengths r\ Substituting (2.20) into (2.15) and (2.16) and utilizing the invariance property of the volume elements, the discretized ODEs for the particle positions and strengths are dx* * = E / • V x Gs (** - If = ( '' ) E r V V \ • T>, x G (*< - x*) r. s (2.21) (2.22) Chapter 2. THE VORTICITY METHOD 16 An alternative, but equivalent, method of deriving the discretized velocity expression illustrates another advantage of recasting the ODE (2.16) in terms of the strengths P . Further, this alternate method also gives a second interpretation of the distribution function ij> . Previously, ij> was introduced as a modification to the Biot-Sarvart Law s s in order to remove the singularity arising from the fundamental solution. In the present case, ijjg represents a smoothing function for discrete point vortices. Since the relationship between the vorticity and velocity is a linear one, the total velocity can be approximated by a superposition of discrete patches, or 'blobs' of vorticity. Each blob is associated with a fluid particle and is centered at the particle's position. The total vorticity is represented by a summation over all of the blobs. At a particular field point SB, the vorticity is approximated by, «(«,«) = f ; ^ ( « - « i ) r i ( * ) - (2-23) The function i> (r), where r = |sc — SB |, is the distribution function defined by (2.18) and j s (2.19). In this context, it describes the shape of the vortex blobs, each of which has a vector strength P . The function ij) smears the vorticity around the source point SE with j s an effective radius related to the parameter 8. For this reason the function ip is often s referred to as the vorticity distribution function. The second advantage of formulating the vorticity method in terms of the strengths is now apparent. Given the strengths P , the vorticity may be evaluated at any point in the fluid SB by (2.23). The nature of the vorticity distribution function, as shown in Figure 2.1 is important to consider when modelling fluid flows. The strength of a blob, by being smeared, will not be constant throughout the volume and since is not identically zero outside a blob defining radius, every blob will contribute to the vorticity throughout the flow. When particles are near boundaries, a certain amount of vorticity will extend across the Chapter 2. THE VORTICITY METHOD 17 boundaries. However, the decay of the distribution function is sufficiently rapid that far field vorticity contributions of a blob near a boundary will be very small. Substituting (2.23) into the velocity expression (2.10) results in TV r (x,t) = J V x G (x - x') i> (x' - x ) T (t) dx'. N j j s (2.24) Interchanging the summation and integration, and taking the derivative with respect to the unprimed coordinates outside the integral, yields u(x,t) = £ V x If G(x- x') j> [x' - x ) dx' T (t) . j j s (2.25) After a change of variables, the remaining integral is seen to be the convolution of the fundamental solution G with the distribution function ip . For a discrete set of field s points x the result becomes l u (x\ t) = f) V x G («•' i=i s V (t) , (2.26) in agreement with (2.20). The two methods presented here have led to the same result: an expression whereby a velocity field can be calculated from an initial vorticity strength field. This leads directly to a system of discretized ODEs forming the basis of the vorticity method algorithm. The ODEs, (2.21) and (2.22), governing the fluid flow are ideally suited for time stepping by a method such as Runge-Kutta. However, there are still two important issues remaining. In order to complete the algorithm, a discussion of the evaluation of the initial vorticity strengths and the derivatives necessary to advance the ODEs is required. Chapter 2. THE VORTICITY 2.4 METHOD 18 Development of the Algorithm 2.4.1 Initialization The utility of (2.23) is further illustrated in the context of initializing a vorticity strength field. The initialization requires determining the strengths from a prescribed distribution of vorticity. This inversion for I", as suggested by Knio and Ghoniem [1990], makes use of (2.23) written in component form as «i = E ^ ( » ' - . « ) S » i <- r 2 27) 3 where the subscript k denotes the vector component. Further, let the matrix St be denned by * , i = ft («* ~ «'') - (- ) 2 28 which characterizes the contribution of the strength at x to the vorticity at SB*. Then, 3 (2.27) becomes the linear system, o>o = *r , o where the vectors o» and r 0 o (2.29) correspond to the fcth component of the initial vorticities and strengths for all thefluidparticles. Specifying the initial vorticities and solving (2.29) for each component gives the desired initial strengths. If there are N fluid particles, the matrix ^ has dimension N x JV. For large numbers of particles the initialization of T becomes a computational burden. Since the vorticity } distribution function decays rapidly, many of the elements of * will be very small. The only additional characteristic of * which might be exploited is its symmetry. Under these circumstances, the only efficient method of solving the initialization system is to Chapter 2. THE VORTICITY 19 METHOD treat the matrix ¥ as sparse. This is accomplished by setting a relative cut-off value e and setting all those elements ^ / ^ ( O ) < e equal to zero. The computational overhead of initialization is reduced dramatically. The sparse method was verified by initializing a single, straight vortex filament (§ 4.3) using both LU decomposition on the full matrix and conjugate gradient method on the sparse matrix. Three different cut-off values, varying through four orders of magnitude, were tested and all produced good agreement with the full L U decomposition. Since there was no discernable difference between the two lower cut-off values, the middle value e = I O was adopted as the standard. -6 2.4.2 Evaluation of Derivatives The evaluation of the RHS terms of (2.21) and (2.22) may at first appear cumbersome. The modified Green's function involves an analytical convolution of the fundamental solution G with the distribution function i{> . Spatial derivatives are then carried out on s G . The introduction of second and third order tensor operators makes the differentiation s very straight forward, while the analysis of Beak and Majda [1985] provides a slick way of finding G . 6 To begin, recall the discretized Biot-Savart law (2.20) and note that the curl with respect to the field coordinate x operates on the function G (x — X*) only. For an x x s arbitrary strength vector T, the curl in (2.21) can be written using indicia! notation as V xG T s dG, — 6: dxm (2.30) It is convenient to introduce the tensor klm € dG s (2.31) Chapter 2. THE VORTICITY METHOD 20 so that VxG T = K T. s u (2.32) l To make explicit the fact that the gradient is evaluated with respect to the coordinate SB* and that T corresponds to the strength at x , the curl in (2.21) becomes 3 VxGjI*-tfgrj. (2.33) when ~ = ~&k ( * " *') ' €klm Gs ' m (2 34) Note that there is no implicit summation over the particle labels i and j. Using the definition (2.34), the trajectory ODE in component form is, ^ = E^r/. (2.35) Similarly, the two spatial derivatives in (2.22), for the vorticity strengths are expressible as a third order tensor. Recognizing that the RHS of (2.22) involves dx'/dt, the evolution equation for T* can be written as dT i N "IT = d t E I t kim Ti, L i=x (2-36) where the third order tensor L is kim L = • m (2.37) Chapter 2. THE VORTICITY 21 METHOD It now remains to evaluate G so that explicit expressions for K and L can be s obtained. As shown by Beak and Majda [1985], the convolution resulting in the nonsingular G can be accomplished by a multiplication of the elements of IT by a smoothing s function / . The smoothing function is related to the distribution function if> by (2.38) The trajectory tensor K is then evaluated according to ifS= where r = x l m —x 3 m m and r 2 -^ »" » G)' e = rr. m m r / - (2 39) The tensor L is evaluated by applying the derivative operation given in (2.37) to the above expression, (2.39). The position and strength derivatives are found by contracting the tensors K and L with the vorticity strengths and summing over all source blobs according to (2.35) and (2.36). 2.5 Accuracy and Convergence In § 2.3 the rotational part of the fluid domain was discretized into a finite number of fluid particles with corresponding blobs of vorticity. The accuracy of the vorticity method depends partly on the discretization of the domain of interest. In considering accuracy and convergence, each blob of vorticity can be considered a small sphere with a radius S and their centers are spaced at an approximate distance of h. In order to develop an accurate algorithm, the vorticity distribution function t(>, defining the shape of each blob, must be chosen to satisfy certain conditions. Moreover, the parameters S and h must be related in such a way as to produce a smooth distribution of total vorticity. There is extensive literature analysing consistency, convergence and stability of vorticity methods (see, for example, Anderson and Greengard [1985]). In general, convergence Chapter 2. THE VORTICITY 22 METHOD proofs require the small parameters 8 and h, which characterize the spatial discretization, to be related by 8= h , q (2.40) where q € (0,1). The appropriate value of q within this range depends primarily on the vorticity distribution function [Puckett, 1991]. Given a proper choice of 8 and h, Beale and Majda [1985] give the conditions that must be satisfied by the distribution function rp(x) for the vorticity method to be convergent. First, the distribution function must be "smooth and rapidly decreasing." Second, its integral over 3-space must equal unity. Finally, its first moment with respect to r = \x\ must vanish. The vanishing of higher moments improves the accuracy. Indeed, the order of the method is given by an integer m, which depends on the largest vanishing moment. An order m method has an error of O 8 . m Recall that the chosen distribution function, defined in (2.19), is a function of radial distance only. The following results regarding functions of this type are due to Beale and Majda [1985]. A distribution ip which satisfies the first two conditions above and is a function of radial distance only is at least order m = 2. However, convergence of the vorticity method when taking the analytical approach to the spatial derivative in the vorticity strengths evolution equation (2.22), requires that rp be order m > 4. Further, distributions ij> which depend on r exponentially, the so called "cubic gaussian," require .<!•• The approach adopted in this study involves choosing a particular 8 and computing a suitable value of h. It's been found that a convenient choice is 8 = hS. (2.41) Chapter 2. THE VORTICITY METHOD 23 Other choices of q were explored, but those leading to smaller values of h increase the number of particles required in the discretization of the fluid's vorticity. Computations then become cumbersome, since the run time is proportional to the square of the number of particles. Hence, it is desirable to choose an h that will give reasonable convergence while keeping N as small as possible. The utility of (2.41) was verified numerically by repeating calculations with h halved or quartered. This is equivalent to doubling or quadrupling the number of fluid particles in the computation. Good agreement between results with increasing N indicates that the solutions are convergent. Although the analysis of accuracy and convergence requires (2.40), dimensional arguments suggest that the solution accuracy can be recast in terms of the ratio h/S. Condition (2.40) then requires 7<1. o (2.42) The physical interpretation is that, given the size 8 of vortex blobs, h must be chosen such that the overlap of adjacent blobs is sufficient to produce a relatively smooth, continuous distribution of vorticity throughout the rotational part of the fluid. 2.6 2.6.1 Boundary Conditions Inviscid F l u i d Flow The Biot-Savart law (2.10) gives a velocity field which is consistent with the vorticity field. However, this result does not necessarily satisfy the boundary conditions for inviscid fluids: the normal component of velocity must vanish everywhere on the boundary. Let the velocity field, found by (2.10), be called u. Then, in general, on any boundary whose unit normal vector is n, Chapter 2. THE VORTICITY METHOD u-nfO. 24 (2.43) The required boundary condition on u will be satisfied by the addition of an irrotational velocity field, or a potential flow, tt. It is called potential flow because the velocity is given by the gradient of a scalar potential, u = V<f>. (2.44) Applying the incompressible condition (2.5) shows that the velocity potential satisfies Laplace's equation. The boundary conditions on the velocity potential are determined by requiring the normal component of the total velocity tt + it to vanish everywhere on the boundaries. Hence, the potential <f> satisfies VV = 0, V<f> • n = - t t • n . (2.45) (2.46) The total velocity, tt + ti, satisfies the inviscid boundary condition of zero normal flow. Further, since the curl of a gradient is necessarily zero, the additional potential flow contributes nothing to the vorticity field. The two velocity fields tt and tt correspond to the particular and homogeneous solutions of Poisson's equation (2.8). 2.6.2 M e t h o d of Images As described above, a homogeneous solution to a partial differential equation can be constructed to satisfy specific boundary conditions. The method of images is an analytical technique designed to find such a solution and is common in seismology and electrostatics. It is particularly useful when the bounding surfaces have a high degree of symmetry, Chapter 2. THE VORTICITY METHOD 25 such as planes and spheres. The use and construction of Green's functions as a means of solving partial differential equations is closely related to the method of images. Hence, a brief examination of Green's functions will illuminate the physical interpretation of the method of images. The fundamental solution, or Green's function, is denned as the solution of a partial differential equation, such as Poisson's equation (2.8), due to a unit point source (of heat or electric charge, for example). The RHS of Poisson's equation is often referred to as the source term. When the source function is zero everywhere Poisson's equation reduces to Laplace's equation. In the case of thefluidflowproblem presented in § 2.1, the source term for the vector potential is the negative of the vorticity. In electrostatics, the method of images is very useful for finding electric potentials due to point charge distributions. In this context Poisson's equation governs the electric potential and the source term is proportional to the charge distribution. Let the domain of interest be the upper half space z > 0, denoted by T>. Let the boundary z = 0, denoted dT>, be a conducting surface, so that the electric potential is zero everywhere on dT). For a point charge q located at the point (0, 0, H) in T>, the electric potential satisfies Poisson's equation with a source function / ~ qS(x — x'), where x' = Hz. By definition the fundamental solution gives the electric potential due to the point charge q, however, this solution does not satisfy the boundary condition. By the principle of superposition [Griffiths, 1989], a second charge, — q, can be added at the point (0,0, —H) and the total potential is the sum of the two individual potentials. The symmetry of the problem guarantees the boundary condition is satisfied for all z = 0, and since the second charge is outside the domain V, Poisson's equation is still satisfied in V. The two potentials can be regarded as the particular and homogeneous solutions of Poisson's equation. Given a fundamental solution, such as the electric potential of a point charge, it Chapter 2. THE VORTICITY METHOD 26 is straight forward to calculate a solution for a continuous charge distribution. The particular solution typically involves a convolution of the fundamental solution G with the source function / over the domain V. Equation 2.9 gives such a solution for the vector potential $ in terms of the vorticity distribution u>. In terms of the theory developed previously in this section and in § 2.1, the velocity field prescribed by (2.10) is a particular solution to the fluid flow problem. Satisfying the inviscid boundary condition by adding a potential flow is equivalent to finding a homogeneous solution to Poisson's equation, and superposing it with the particular solution. Considering the discretization of the vorticity in § 2.3, with delta function distributions, the application of the method of images to this problem becomes more clear. For a single point vortex in the half space described above, the condition of zero normal velocity on a planar boundary is satisfied by the addition of a single image point vortex, as was done with the image charge. There is one difference; in the electrostatic case the image charge is a scalar, whereas in the point vortex case it is a vector quantity. An arbitrarily directed vorticity vector associated with the point vortex can be decomposed into two components, one perpendicular to the boundary and the other parallel. When the flow field is due to the point vortex only, the normal flow at the planar boundary is due exclusively to the component of vorticity parallel to the boundary. Hence, the image vortex must have an oppositely directed parallel component. When point vortices make up a vortex filament in X>, the continuity of u» • n across the boundary makes it necessary to have the perpendicular component of the image vortex in the same direction as the interior vortex. This point will be elucidated in the discussion of waves on a vortex filament. The image vortex generates a potential flow in V, and the associated velocity field is equivalent to the homogeneous solution. Two of the applications of the vorticity method described in Chapters 4 and 5 consider fluidflowin a region between two parallel planar boundaries. The addition of a second Chapter 2. THE VORTICITY METHOD 27 4 Distance from boundary zlh Vortex blob ! A B C h Ah lO/i i 0 Figure 2.2: Geometry of a single blob and planar boundary. The normal component of velocity is evaluated at three points A, B and C along the boundary. boundary parallel to the first, with T> now being the region 0 < z < d, complicates matters. In this case, the single point vortex in T> must be imaged across each boundary. Further, since each image will generate a normal velocity at the far boundary, there will have to be an image for each image across the opposite boundaries. Following the same reasoning, the images for the images will require images, and so on. Thus, in theory, the addition of a second boundary will force the requirement of an infinite number of images for each point vortex in T>. However, in practice it suffices to retain a small number of images outside each boundary. The inviscid boundary condition, u • n = 0, may be approximately satisfied by a finite number of image blobs outside a planar boundary. A simple numerical experiment illustrates a single blob's decreasing contribution to the normal velocity as its distance Chapter 2. THE VORTICITY METHOD 28 from the boundary increases. Figures 2.2 and 2.3 show the velocity decrease at three different points along the boundary. The velocities were evaluated for a blob with unit strength directed parallel to the boundary. In all cases the velocity is small when the particle is several 6 from the boundary. That is, only the few particles nearest the boundary contribute significantly to the normal velocity. Theory indicates the need for an infinite number of image blobs in a two boundary problem. However, provided the boundaries are sufficiently far apart, this result shows that it is possible to make the normal component of the boundary velocity arbitrarily small with relatively few image blobs. Chapter 2. THE VORTICITY 29 METHOD 0.4 10.0 Distance From Boundary z/h Figure 2.3: Normal velocities at planar boundary. Three velocity profiles evaluated at A, B and C in Figure 2.2. Chapter 3 VORTEX FILAMENTS: BACKGROUND 3.1 Physical Properties The implementation of the vorticity method, as described in the previous chapter, must be tested against known solutions. It is preferable to progressively test the algorithm by modelling' increasingly complicated physical situations. The test cases in the following chapter involve flows generated by single vortex filaments, each of which is well known in vortex dynamics. As an introduction to these test cases, this chapter presents a brief summary of some basic concepts of vortex dynamics (see Batchelor [1967], Tritton [1988] and Saffman [1992] for more detailed discussions of vortex dynamics). These concepts serve as an introduction to the description of a vortex filament. A vortex line is a curve through a fluid whose tangent is everywhere parallel to the local vorticity vector. Such a line is found by solving a set of difFerential equations in a manner analogous to calculating stream lines [Batchelor, 1967]. Vortex lines within a fluidflowcan be visualized in much the same way as electromagnetic field lines. A higher density of vortex lines indicates a general increase in the vorticity in that region. Hence, a uniformly rotating fluid will have a uniform distribution of vortex lines. A consequence of the form of the vorticity equation for inviscidflow(2.6) is that a vortex line deforms in the same way as a material line of fluid. This is immediately evident because a fluid line element evolving with a flow satisfies the same differential equation [Batchelor, 1967]. The implication is that a vortex line, in an inviscid fluid, always consists of the same 30 Chapter 3. VORTEX FILAMENTS: BACKGROUND 31 fluid particles. Consider two adjacent particles on a vortex line as shown in Figure 3.1, where the particle Q is moving relative to P with a differential velocity du. The decomposition of du into components normal and tangential to the vortex line shows that each contributes to the change in local vorticity in differing ways. The tangential component extends or compresses the vortex line, whereas the normal component rotates the line element between P and Q. Hence, according to (2.6), the vorticity changes as a result of stretching and bending of vortex lines due to strain associated with the flow. Figure 3.1: Bending and stretching of a vortex line. The vortex line is deformed by a differential velocity between the points P and Q. Adapted from Botchelor [1967]. Another useful concept is the vortex tube. A vortex tube is a region of fluid defined by a collection of vortex lines (see Figure 3.2). By definition, the divergence of the vorticity is necessarily zero. Further, the vorticity flux across the sides of a vortex tube is zero, since the vorticity is always tangential to the sides. Applying the divergence theorem to the vorticity over the volume of a vortex tube truncated at two arbitrary surfaces, A 1 and A , along its length shows that the vorticity flux across the two surfaces is equal, 2 Chapter 3. VORTEX FILAMENTS: 32 BACKGROUND Figure 3.2: Illustration of a vortex tube. The cross sectional areas A are enclosed by the material curves 7^. { (3.1) Hence, the strength of a vortex tube, denned as (3.2) A where A is any surface intersecting the tube, is constant everywhere along the tube. This implies that a smaller cross-sectional area must result in an increased average vorticity over the surface A, which is consistent with the idea that vortex line density is related to the magnitude of the fluid vorticity. Applying Stokes' theorem to (3.2), yields Chapter 3. VORTEX FILAMENTS: BACKGROUND = J u • dl, 33 (3.3) where 7 is the material curve enclosing A. This is commonly called the circulation. Kelvin's circulation theorem states that C is invariant in time when the fluid is inviscid [Batchelor, 1967]. Therefore, as a material curve evolves in an inviscid flow, the circulation around it remains constant. The circulation associated with a vortex tube also remains constant through time, as does the tube's strength. The implication of Kelvin's theorem is that a closed material curve always encloses the same vortex lines and a vortex tube moves with the fluid. Further, if a vortex tube becomes elongated, its cross section must also contract due to conservation of mass. And by the properties of the tube's strength the vorticity must increase in the stretched region. Hence, the stretching of a vortex tube or line corresponds to an increase in the local vorticity. If the fluid particles are taken to be points separated by an infinitesimal material line element, the relative change in the vorticity due to extension or compression is equal to the relative length change in the fine element [Batchelor, 1967]. Very often in fluid mechanics the vorticity is confined to (long) slender regions, or inother words, a thin vortex tube. A useful mathematical tool in studying such regions of vorticity is the line vortex. A line vortex is a singular distribution of vorticity along a curve within a fluid, and must not be confused with a vortex fine. For an infinitely long, straight fine vortex, fluid particles will follow circular paths centered on, and in a plane perpendicular to the line. Thus, the only non-zero velocity component is the azimuthal one, which is given by K u = 2*r ' 6 0 _ (3.4) where K is the line's strength. The strength is equivalent to the circulation of the line, Chapter 3. VORTEX FILAMENTS: 34 BACKGROUND calculated around any closed curve 7 through which the line passes, u • dl. (3.5) 1 The vortex line's azimuthal velocity structure decreases like 1/r, and though the line vortex is an important analytical tool, it is not physically realistic since u —• 0 0 as e r -> 0. A vortex filament is, in a sense, a physical realization of a line vortex. Two common approximations for analysis of vortex filaments are the uniform and hollow core filaments. These two approximations have vorticity evenly distributed across the core and concentrated on the surface of the core, respectively. Analogous to a line vortex, an infinitely long, straight vortex filament with a uniform vorticity distribution has only an azimuthal velocity component, which is given by Qr u =< e [ SIT\JT r < r. ~ r > r , (3.6) e where r is the average core radius and £2 is the rotation rate of the fluid within the core. e This velocity profile is shown as a dashed line in Figure 3.4. The long, thin core region with radius r rotates as a solid body. Outside the filament core, the flow is irrotational, e and the azimuthal velocity decreases with r in the same way as the line vortex. Theoretical studies of vortexfilamentssometimes allow for small perturbations to the velocity field due to variations in thefilament'svorticity distribution or core radius. The introduction of curvature in an initially straight filament also perturbs the velocity field. The curvature associated with a bendingfilamentcauses thefilamentitself to move. This is known as self-induction of a vortex filament. Chapter 3. 3.2 VORTEX FILAMENTS: 35 BACKGROUND Discretized Vortex Lines and Filaments Since a vortex line always consists of the samefluidparticles, it is natural to discretize a segment of the line into a finite number of fluid particles separated by a distance h. Using the vorticity method, there will be a blob of vorticity associated with each of these particles. A uniform filament is a collection of one or more parallel vortex lines. Hence, it can be modelled, using the vorticity method, with one or several discretized vortex lines. In the following test cases, perturbations to the filament core are not considered, thus making it possible to model the filaments using a single discretized vortex line. A discretized vortex line with h/8 = 1 is illustrated in Figure 3.3. As h/6 decreases, the overlap of adjacent blobs increases and the vorticity distribution along the vortex fine's length becomes smoother. Since the blobs making up the vortex line have afiniteeffective radius, the discretized vortex line represents a filament with a core radius that is roughly comparable to the radius of the discrete blobs. The vorticity distribution, as a function of r, within each blob is given by (2.18) and (2.19). As shown in Figure 2.1, the vorticity distribution function ip differs from the g simple examples above in that it is non-zero for all r. In other words, each fluid particle will contribute vorticity, and hence, contribute to the velocity field everywhere in the flow. Although the blobs contribute to the vorticity everywhere in the fluid, the main contribution is within the blob radius 8. Consequently, it is possible to accurately reproduce the behaviour expected from a well defined, finite region of vorticity. In order to compare computed results with theory, it is essential to define an effective radius r . The effece tive radius is equivalent to the radius of a uniform vorticity distribution which results in similar filament behaviour and flow characteristics. Given the mathematical form of ip , an obvious, but somewhat arbitrary, choice would be r = 8. A better, and more s e Chapter 3. VORTEX FILAMENTS: BACKGROUND 36 Figure 3.3: Illustration of a discretized vortex line. The arrows are tangential to the vortex line at the particles' centers and represent the direction of each particle's vorticity. Figure 3.4: Azimuthal velocity structures of straight filaments. The dotted line corresponds to the velocity in (3.6), with = |w. The solid line corresponds to a discretized filament. physically meaningful, choice can be found by comparing the velocity structures of two vortex filaments. Consider two straight, long vortex filaments: one with a uniform vorticity distribution as given by Equation 3.6; the other discretized with each particle's vorticity distributed according to ip . The two azimuthal velocity structures are presented in Figure 3.4. In s the case of the uniform filament, the core radius is marked by the peak in the velocity structure. The core radius and vorticity of the uniform filament were chosen such that Chapter 3. VORTEX FILAMENTS: 38 BACKGROUND the peaks of the two profiles coincide. The similarity of the profiles suggests r be chosen e according to a maximum in azimuthal velocity. Radius r Figure 3.5: Determination of the effective radius. For each value of 6, the peak in azimuthal velocity is located at approximately |5. In order to determine a general relation between this peak and the parameter 6, several cases were examined. For each of the four cases shown in Figure 3.5, the peak occurs at approximately | £ . This is the adopted value for the effective radius in the following computations involving discretized filaments. By assuming the circulation of the discretizedfilamentis equivalent to that due to a Chapter 3. VORTEX FILAMENTS: BACKGROUND 39 constant vorticity over an area of xr , an effective vorticity, o>, can also be determined. 2 e The effective vorticity of the discretized filament and the vorticity of the uniform filament are quoted in Figure 3.4. Hence, using the effective radius to determine an effective vorticity facilitates the comparison of computed results with theory, since the uniform filament with equivalent vorticity, radius and, hence, circulation is expected to behave in a similar manner. 3.3 Computing the Circulation of a Discretized Filament Section 3.1 shows that there are two different, but equivalent, methods of calculating the circulation of a vortex filament. Thefirstis expressed as an area integral of the vorticity flux in (3.2), while the second involves a line integral of the velocity in (3.3). For a discretized filament, it is favourable to use the former since the expression for the total vorticity (2.23) can be used to evaluate the vorticity over the area A. The integral is taken over an area in a plane perpendicular to the vortex filament. That is, the vector dA is parallel to the local vorticity of the filament. The calculation is simplified by converting the coordinate system to cylindrical polars (r, 9, z) such that the origin is at the center of the filament and z is parallel to dA. After calculating the azimuthal part of the area integral, the circulation integral becomes oo (3.7) o The integral in (3.7) is evaluated numerically as the finite summation C = 2ir£>,(r ) Ar, i where the vorticity is evaluated by (3.8) Chapter 3. VORTEX FILAMENTS: BACKGROUND «..(»•*) = E ^ ( ' ^ « i ) ( i - « ) . r 40 (-) 3 9 according to the discretized vorticity (2.23). Once again taking advantage of the decreasing nature of the vorticity distribution function if> , the summation in (3.8) is terminated s once the relative contributions of the vorticity to the circulation become very small. Circulation values referred to in the following chapter were computed numerically using the trapezoidal rule. The integration was verified by summing to r = 506 and successively doubling the number of steps until convergence was reached. Once an appropriate step size was established, a terminal radius was selected by starting with a small value and increasing until the result again converged. The values which were found to produce accurate and efficient results were 500 steps per 6, extending to a radius of 355. Chapter 4 VORTEX FILAMENTS: NUMERICAL MODELS 4.1 Self Advection of a Vortex Ring The first and simplest test case is the modelling of a closed, circular vortexfilament,or vortex ring. It's the simplest case because the filament is closed and so there is no need to consider boundaries or the practicalities of dealing with long, unclosed filaments. This particular type of fluid flow is commonly seen as smoke rings. In this situation there is a thin, dbnut shaped region of the fluid which is rotating and moving through an otherwise irrotational background. Due to its curvature, a vortex ring will move with a constant velocity in the axial direction. The direction of motion is found by a right handed rule. With the fingers curling around the ring in the direction the vorticity, the thumb gives the direction of motion. A vortex ring is illustrated in Figure 4.1. Vortex rings have been well studied, both analytically and numerically. Examples are the merging of two vortex rings [Anderson and Greengard, 1988] and the interaction of four vortex rings in a plane, at the vertices of a square [Leonard, 1985]. The objective of this section is not to model complicated dynamics or interactions of vortex rings, but to test the algorithm by modelling rings of different geometries. The self advection velocities computed using the vorticity method will be compared to approximate theoretical velocities. The theory applies to thin rings, where thin implies that a dimensionless aspect parameter—the ratio of the core radius to the ring radius—is much less than one. For a vortex ring of radius R, which has an effective core radius of 41 Chapter 4. VORTEX FILAMENTS: NUMERICAL 42 MODELS Figure 4.1: Illustration of a ring vortex. The arrows at the sides indicate the sense of fluid rotation in the ring's core. r , the aspect parameter is defined as e A = (4.1) R Saffman [1992] gives the velocity of a thin ring vortex with a uniform vorticity distribution in terms of the ring's circulation C , radius R and core radius r . Rewritten in e terms of the aspect parameter A , the ring velocity is U = Saffman also quotes the error, 0(CA /r 5 e CA (-1-3 • (4.2) due to Dysan, in this approximation to be log(8/A)). The relative error in the theoretical velocity is then proportional to A , which means the approximate ring velocity, as given by (4.2), becomes less ac4 curate as A increases. Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 43 The computational vortex rings were initialized as single filaments of circular shape. For a ring with a given radius, the velocity was computed by modelling the resulting flow through a time interval chosen as follows. The time scale of the problem is T = where u; is the magnitude of the ring's vorticity. An integration interval of several T was sufficient to capture the general behaviour of of the flow. After a time of IOT the average axial displacement was computed over all the particles comprising the ring, and an average velocity was found. The convergence of the ring velocity with the number of particles used in the computation is shown in Figure 4.2. In this particular case, the ring geometry is fixed with R = 5 and S = 0.5 (r = 0.375). The number of particles was varied by changing h, the e specified distance between particles, or equivalently, varying the ratio h/S. As h/S decreases, the number of particles comprising the ring increases and the velocity converges to a value very close to the region predicted by theory. Good convergence was reached with N = 158, corresponding to h = |£. The long term stability of the ring was checked by carrying out the time integration for 0(T x 10 ). The average velocity values were not significantly different and the ring 2 retained its shape and structure over the entire interval. Testing the numerical algorithm against theory involved computing several different ring velocities over a range of A , each with the same effective core radius and vorticity. While these two parameters remain constant, the circulation of the discretized ring varied unexpectedly with radius R. The variation can be explained by examining the vorticity distribution within computational rings with differing A . First, consider a vortex ring with a uniform vorticity distribution. The exact circulation can be calculated over an area A enclosed by any curve which threads the ring, in the same way a chain link threads its neighbour. The circulation will be exact because the curve encloses all the vorticity of the uniform distribution. However, with the discretized Chapter 4. VORTEX FILAMENTS: NUMERICAL 44 MODELS ° Computed Velocities -• Theoretical Estimates 0.0215 fc) .« •S > 0.0210 G 2 0.0205 I I I I I I i I 0 50 100 150 200 250 300 350 Number of Particles N Figure 4.2: Convergence of computed ring velocities. The dotted lines give the upper and lower bounds on the theoretical estimate. ring, the vorticity is not confined to a finite region. Suppose the circulation is to be evaluated over a circular, planar region, which has a radius of a few r (see Figure 4.3). e Given the nature of the vorticity distribution function rf) such a curve will enclose most, Sl but not all, the vorticity of the ring. Moreover, there will be a potentially significant contribution to the circulation from the vorticity associated with the particles on the other side of the ring. The difference between the resulting circulation for uniform and discretized rings Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 45 Figure 4.3: Computing the circulation of a ring vortex. The curve 7 is shown with a cut-away view of the ring. lies in the vorticity u> of (3.7) and (3.8). For the geometry pictured in Figure 4.3, z u) is, in fact, a function of 9 as well as r, where 9 and r are the polar coordinates 2 in the plane of integration. Since the area enclosed by 7 must be reasonably large to compute the circulation accurately (A = TT(50S) for the straightfilamentin § 3.3), there 2 will be additional vorticity flux over A which originates from the opposite side of the ring. Moreover, the largest contribution will come from the particles directly across the ring, where o> • dA goes through a maximum in magnitude. For moderately sized rings (A ~ 0.1), the 9 dependence of the vorticity function is noticeable. As the ring's radius increases, this dependence becomes less significant as a result of the decreasing distribution function ib . Hence, the computed circulation will increase in accuracy as R f increases. With the preceding argument in mind, the circulation was fixed at C = 0.2907, the value computed for the largest R. The theoretical upper and lower bounds were then calculated by adding (or subtracting) the error estimate to (4.2). The resulting velocity functions are plotted as dashed lines in Figure 4.4 with the computed velocities as circles. Chapter 4. VORTEX FILAMENTS: NUMERICAL I I 46 MODELS i i i ..•••'6 0.040 0.035 0.030 <D > I? 0.025 2 ,.Or'' - 0.020 - 0.015 o Computed Velocities Theoretical Estimates .0 O'' 0.010 0.04 i 0.06 —.. .1, 0.08 1 : i i 0.10 0.12 0.14 Aspect Parameter A i 0.16 i 0.18 0.20 Figure 4.4: A comparison of theoretical and computed velocities. The dotted lines give the upper and lower bounds on the function 17(A) given by (4.2). The agreement is very good and the relative error ranges from 0.13% to 0.30%. There are smallfluctuationsin the errors due to different convergence properties for each A. However, there is also an expected increasing trend with A . A summary of the relative errors in ring velocities is presented in Figure 4.5. Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 47 0.0030 0.0025 I g •a 0.0020 0.0015 0.0010 0.04 0.06 0.08 0.10 0.12 0.14 Aspect Parameter A 0.16 0.18 0.20 Figure 4.5: Relative error in ring velocities. Errors were determined using the theoretical approximation (4.2) with constant circulation, C = 0.2907. 4.2 Interaction of Two Vortex Rings A more qualitative test case involves the interaction of two similar vortex rings. The rings, with a common axis of symmetry, pass alternately through each other's centers. The 'leap frogging' behaviour of the rings illustrates the utility of the vorticity method for modelling more complexflowsin which the vorticity is confined to small finite regions. In such cases, the vorticity method is very efficient, even when there is a large particle density and, hence, high resolution. Chapter 4. VORTEX FILAMENTS: 5.0 NUMERICAL -5.0 0.0 5.0 48 MODELS -5.0 0.0 5.0 4.0 10.0 5.0 3.0 0.0 2.0 1.0 -5.0 0.0 -10.0 5.0 10.0 4.0 5.0 3.0 0.0 2.0 1.0 -5.0 0.0 -10.0 5.0 10.0 4.0 5.0 3.0 0.0 2.0 1.0 -5.0 0.0 -10.0 5.0 10.0 4.0 5.0 3.0 1 2.0 M 1.0 0.0 0.0 -5.0 -5.0 0.0 5.0 -5:0 0.0 5.0 -10.0 Figure 4.6: Evolution of two ring vortices. The left and right panels show the side and plan views of ring positions, respectively. Top to bottom, t = 1.00, 4.11, 5.46 and 7.44. Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 49 Consider two vortex rings placed near one another, with a common axis of symmetry, equal core radii r and similar radii R. By themselves, the rings will travel in the axial e direction at a constant rate. However, the interaction of the superposed velocity fields will cause the radius of the rear ring to decrease, and the radius of the forward ring to increase [Bachelor, 1967]. This can be seen from the sense of the velocity field near a ring's core, as shown by the curved arrows in Figure 4.1. Since a vortex ring's speed will increase as its radius decreases, the rear ring will speed up and pass through the larger one. Once the rings have traded places, the radius of the smaller, forward ring will then begin to increase. Meanwhile, the radius of the rear ring decreases, ultimately becoming less than that of the forward ring and the process is repeated. This interaction of vortex rings is one which has been verified by experiment [ Yamada and Matsui, 1978]. In realfluids,viscous effects cause the vorticity of the rings to diffuse quickly. Hence, the rings in the experiments exchange places only two or three times before their cores disperse. By contrast, the interaction can continue indefinitely in an inviscid fluid. Two identical vortex rings were initialized as singlefilamentsof circular shape. Each had an effective vorticity u> = 6.5 and an effective core radius r = 0.375. With their e e symmetry axes aligned with z, the two rings were placed at z — 0 and z = 2 with radii R = 5. The particle density was increased by a factor of | over that required for acceptable accuracy. This ensured the ratio h/8 remained sufficiently small when the radius of either ring reached a maximum. The time integration was carried out over an interval of t — [0,25], spanning one and a half complete cycles. Over this period, the rings showed no signs of significant deformation or dispersion. A sample of the results is shown in Figure 4.6. Each pair of plots is a snapshot in time of the rings' positions as viewed from the side and the top. Throughout the figure, each ring is distinguished by a solid or dashed line. The period over which the snapshots Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 50 were taken represents approximately half a cycle. Although the rings are pictured as lines, the scale of the core radii can be established from the first snapshot. Here, the difference in radius of the two rings is 0.4; approximately r . Hence, the distance between the two lines in the right panel is about the same as the effective core radius for these two rings. At time t = 5.46, as the smaller ring passes through the larger, their radii are R = 4.2 and R = 5.7. The difference is Ai2 = 1.5, indicating that the cores do not interfere with each other as the rings pass. 4.3 Waves on a Bounded Vortex Filament A second quantitative test case is the modelling of standing waves on a vortex filament bounded by two parallel planes. This case is more complicated than the previous one because boundary conditions have to be enforced. There are additional complications encountered in the initialization of a finite length filament. This type of fluid flow is somewhat analogous to the flow associated with a hypothetical stationary tornado. Tornados are often seen bending and precessing about a vertical axis as they move across the ground. The primary difference is that the tornado is not bounded on top by a solid boundary; it is a free boundary. As well, the vorticity structure within the tornado is closer to that of a hollow filament [Vanyo, 1993]. The flow associated with waves on afilamenttakes place in a fluid region infinite in horizontal extent and bounded between two parallel planar boundaries. Perpendicular to the boundaries exists a long, thin, almost straight region of fluid rotating as a solid body within an irrotational background. A small, sinusoidal perturbation in a vertical plane is imposed on an initially straight filament, causing a rotational motion of the disturbance. The curved filament rotates about a central axis much like a skipping rope rotating around a vertical axis, rather than a horizontal one. A countably infinite Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 51 number of modes can be realized on the filament provided there is an integral number of half wavelengths between the boundaries. The particular case shown in Figure 4.7 is a mode two wave, which consists of one complete wavelength between the boundaries. The rotation of the disturbance, or precession, is in a direction opposite to the fluid's rotation in the core. In a discretized filament, each particle travels round a circular path in a plane parallel to the boundaries with a radius varying sinusoidally along the length of the filament. z= d z= 0 Figure 4.7: Illustration of a bounded filament. The solid arrow indicates the sense of the fluid rotation in the filament's core. The dashed arrow indicates the direction of the precessional motion. Kelvin [1880] first analysed this problem with astounding detail. Since then, a number of interesting studies have been undertaken to examine variations of this basic problem. Unlike the vortex ring, this type offilamentis difficult to produce experimentally. In a Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 52 unique study of curved vortex filaments, Hama and Nutant [1961] formulated an analytical solution for the self induced velocity of a parabolic filament. There is a mathematical analogy between vortex dynamics and magnetostatic theory, where the velocity u and vorticity a; are identified with the magnetic field and current density. The workers were able to verify their fluid dynamical result by measuring the magnetic force on a curved wire due to the electric current flowing through it. In a manner similar to that in § 4.1, the algorithm will be tested by modelling filaments of different geometries. The precessional frequencies computed using the vorticity method will be compared to theoretical approximations. Kelvin's analysis applies to thin filaments, where a dimensionless aspect parameter—the product of the core radius and the wave number of the disturbance—is much less than one. For a thin filament with an effective core radius r and a sinusoidal perturbation of wave number k, the aspect e parameter is defined as A = kr . e (4.3) Kelvin [1880] gives the precessional frequency of a thin, sinusoidally perturbed filament with uniform vorticity distribution as fiA 2 * = / 1 \ (in ^-0.3659] , (4.4) where Cl is the rotational rate of the fluid in thefilamentcore. The negative sign indicates a retrograde precession, opposite to the rotation of the filament core. In practice, the frequencies are best computed in terms of the circulation C, the mode of the disturbance m and the boundary separation d. Utilizing the effective vorticity and given that k = mir/d, the frequency in terms of the circulation is Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 53 Ck 2 a = 4* ( l n i-°- 3 6 5 9 )- (-) 4 5 Although Kelvin does not explicitly state an estimate for the error in this approximation, Saffman [1992] quotes an error of O A for the closely related problem of a helical 4 disturbance propagating along an infinite filament. Due to the similarities between the two problems, it is reasonable to assume the error in each of the frequency expressions is of that order. Also note that the error is similar to that for the vortex ring. Incorporating boundaries into the algorithm led to immediate difficulties. Initializing a straightfilamentas a string of equally spaced blobs truncated at two boundaries resulted in a highly oscillatory solution for the vertical component of the strengths. Moreover, the oscillations caused changes in sign of the strengths which completely contradicts the physical interpretation of the discretized filament. Considering that the filament represents a vortex line, it is unreasonable to expect two adjacent blobs on the same filament to have vector strengths with differing signs. The difficulties arose when attempting to initialize a truncated, finite length filament. The initialization process does not include any information about boundaries, and hence, the oscillations were caused when an ideally infinite filament was truncated by multiplication with a boxcar function. In signal processing the effects of finite length records are well known. The multiplication of a signal with a boxcar function in the time domain is equivalent to a convolution with its Fourier transform, a sine function, in the frequency domain. The large oscillatory lobes of the sine function cause severe frequency mixing, or aliasing. The problem at hand is the spatial analogue, where there was aliasing of different wave numbers in the solution which should ideally include only one: wave number zero. The remedy to this problem is to truncate the record (filament) with a function whose Fourier transform causes less aliasing in the frequency (wavenumber) Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 54 domain. This is called windowing the data and its effect is to weight the data less towards the ends of the record, thus reducing edge effects. As such, 'nice' window functions usually have some sort of taper to them, making the truncation less abrupt. The choice of the most appropriate function for a particular application is an art in itself and there are many window functions to choose from. A cosine window was chosen here because it allows a region of constant weighting between the tapered sections and has smooth derivatives. The problem was corrected by extending the initial vortex filament beyond the boundaries and tapering the specified vorticity w to zero. The selected windowing function 0 ramps the vorticity to zero by way of a half wavelength cosine wave. Thus, the derivatives of the windowing function are continuous and the length of the taper could be varied for optimum effect. In terms of the discretized filament, the ramp length was varied as a function of the number of particles. It was found that beginning the taper immediately outside the boundaries didn't fully correct the wavenumber aliasing problem. Further, since the image particles used to satisfy boundary conditions constitute an extended fluid domain, they too must be included in the initialization. Hence, with unit weighting of the specified vorticity between the boundaries, unit weighting on the image blobs immediately outside the boundary and weighting according to the cosine taper beyond the images, a smooth, finite length filament was created. The boundary locations were chosen to Ue midway between the interior end particles and their image particles. After the filament was properly initialized the fluid particles making up the tapered portion were discarded. The total vorticity profile due to the interior particles and the images was exactly as specified at all points in the interior along the filament. The positions and strengths of the interior and image particles were then perturbed. In Kelvin's analysis, an infinitely longfilamentwith a sinusoidal perturbation satisfies Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 55 the inviscid boundary condition if the boundary locations are chosen to coincide with stationary points of the disturbance. If one such boundary is located at z = 0, the other must be at z = d, where d = mir/k. In this case, the filament meets each boundary at a right angle, and hence, the local vorticity is normal to the boundary. The normal velocities on the planar boundaries are zero by virtue of the symmetry of the perturbation across the boundaries and thefilament'sinfinite length. This is an excellent example of the idea of extending the fluid region beyond boundaries to satisfy conditions on the vorticity, as presented in § 2.1. Here, the exterior fluid region is infinite and the flow associated with the vorticity in the extended region ensures that the velocity boundary conditions are satisfied exactly. Computationally, it is not possible to extend the filament beyond the boundaries indefinitely. Based on the results of § 2.6.2, it is possible to model this type of flow with a finite number of image blobs exterior to each boundary. As such, the normal component of velocity at the boundaries will not be identically zero. However, it was found that those velocities could be made arbitrarily small by increasing the number of images. In fact, a surprisingly small number of images is required to produce normal velocities at the boundaries many orders of magnitude smaller than the tangential components. Afilamentinitially aligned with the z axis, extending between two planar boundaries located at z = 0 and z = d, is given a perturbation in the xz plane, where particles are displaced according to Ax = A cos (fez). The amplitude A is very small compared to the wavelength of the disturbance to ensure the linearity of the resulting wave motion. Since the vorticity must always be tangential to the filament, a small perturbation to the vorticity vector is also required. For the spatial disturbance in the xz plane, this corresponds to a small non-zero component of the strength in the x direction. Because the vorticity is a spatial derivative of the velocity, its perturbation can be considered a second order effect on the velocity field. To start with, thefilamentwas perturbed in space only, while Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 56 the strengths remained vertical everywhere on the filament. As a first approximation, this is reasonable. However, this procedure caused a secondary, high frequency motion of the particles superimposed on the expected circular paths. Each particle circled the main retrograde motion in epicycles with an amplitude of approximately 1% of A. It appeared that the small size of the secondary motion did not alter the larger scale behaviour of the flow. However, the epicycles created significant uncertainty in the computation of precessional frequencies. With the proper perturbation applied to the vector strengths, the epicycles disappeared. The precessional frequencies were computed by following the fluid particle nearest the lower boundary. This particle has the largest amplitude motion and the same initial position for every mode. Determining the frequency involves very small time intervals and a Taylor's series approximation for the y component of position, which is initially zero at t = 0. For the particle nearest the lower boundary and a small elapsed time At, the frequency can be computed by where A is the amplitude of the perturbation. The convergence of the wave frequencies was explored with the number of particles N, as well as the specified number of image particles N across the boundary. The results { for the first three modes are demonstrated in Figures 4.8, 4.9 and 4.10. In all cases the filament geometry is fixed with d = 50 and S = 0.5 (r = 0.375). For a particular value e of NJN, the spatial resolution, or variation of h, doesn't significantly change the value of the frequency when N is 200 or greater. Increasing N to 600 did not improve the computed frequencies. Hence, the largest value of h, as found by (2.41), was adequate Chapter 4. VORTEX FILAMENTS: c o o MODELS 57 o Computed Frequencies, N = 200 o Computed Frequencies, N = 400 • Theoretical Estimates -0.00035 * NUMERICAL -0.00036 -0.00037 -0.00038 -0.00039 -0.00040 0.0 0.2 0.4 0.6 0.8 Number of Image Particles N/N 1.0 Figure 4.8: Convergence of computed filament frequencies, mode 1. The dotted lines give the upper and lower bounds on the theoretical estimate. for further computations. A more interesting observation is that the number of images— specifically the length of the image part of the filament relative to d—is more important to the convergence of the frequency. This quantity is expressed as the ratio N^N. When fy/N = 0 there are no images and NJN = 1 means that every interior particle has an image across each boundary. The convergence of the frequencies improves as the number of images increases. Although it is possible to extend the filament beyond NJN = 1, creating the exterior part of the filament by the method of images imposes this limit. Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 58 -0.00115 o Computed Frequencies, N = 200 o Computed Frequencies, N = 400 • Theoretical Estimates -0.00120 I I -0.00125 a o o -0.00130 -0.00135 0.0 0.2 0.4 0.6 0.8 Number of Image Particles N/N 1.0 Figure 4.9: Convergence of computedfilamentfrequencies, mode 2. The dotted lines give the upper and lower bounds on the theoretical estimate. Hence, the computations described below were accomplished with reasonable accuracy using 200 fluid particles in the interior, each having an image across both boundaries. With the boundary conditions satisfied by image blobs, the numerical algorithm was tested by computing the precessional frequencies for the first ten modes of the disturbance. For each computation the effective radii and circulation of the unperturbed filament were the same. For filaments with constant r and fixed boundaries, A varies linearly e Chapter 4. VORTEX FILAMENTS: NUMERICAL 59 MODELS -0.00245 -0.00250 h o c CD a c o -0.00255 r CD I -0.00260 o Gomputed Frequencies, Ar = 200 o Computed Frequencies, N = 400 Theoretical Estimates 7 -0^00265 0.0 1 1 0,2 — —' ' 0.4 0.6 0.8 Number of Image Particles N/N 1 1 1 1.0 Figure 4.10: Convergence of computed filament frequencies, mode 3. The dotted lines give the upper and lower bounds on the theoretical estimate. with mode number. The computed frequencies were compared with theoretical frequencies calculated by (4.5) as a function of the aspect parameter A . The results, together with the order of magnitude error estimates, are shown in Figure 4.11. The agreement between the computed and theoretical values are very good. As shown in Figure 4.12, the relative error in the frequencies ranges from less than 1% to just over 2%. However, in contrast to the ring vortex, these errors behave somewhat erratically as a function of A. The largest error corresponds to the mode 10 disturbance, while the Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 60 0.000 -0:005 r -0.010 -0.015 ° Computed Frequencies a Theoretical Frequencies _0.020 0.00 1 1 0.05 ' 0.10 0.15 Aspect Parameter A 1 1 0.20 -— 1 0.25 Figure 4.11: A comparison of theoretical and computed frequencies. The theoretical values include the order of magnitude error estimates. Increasing values of A correspond to increasing mode numbers. smallest corresponds to mode 3. The fluctuation in the relative error is likely due to the finite length of the filament. For the smaller wavenumbers, the entirefilament,including the images, supports fewer complete wavelengths of the disturbance. The truncated discretizedfilament'appears' more like the infinite theoretical one as the mode number and number of wavelengths along it increase. Hence, there are two effects seen in the relative error: a decrease over the first few modes as the filament becomes longer, relative to the wavelength, and a further increase as the mode number increases. The latter is the Chapter 4. VORTEX FILAMENTS: NUMERICAL MODELS 61 expected behaviour due to the uncertainty of the theoretical approximation. In addition there is a decrease in some errors as the mode number increases from even to odd. This can be attributed to the symmetries of the entire filament, including images. In all computations thefilamentsend at a maximum in the perturbation. In the odd cases, the perturbation does not go through an integral number of wavelengths and the positions of thefilamentends are opposite in sign. For the even modes, thefilamentends have the same phase. Although there is a physical difference between the even and odd cases, it is not clear why the errors are lower when thefilamentends are antisymmetric. Chapter 4. VORTEX FILAMENTS: NUMERICAL 62 MODELS 0.025 0.020 h I g 0:015 i 0.010 0.005 0.00 0.05 0.10 0.15 Aspect Parameter A 0.20 0.25 Figure 4.12: Relative error in filament frequencies. Errors were determined using theoretical approximation with C = 0.2935. Chapter 5 ROTATING CYLINDER 5.1 Introduction In the previous chapter, three computational examples of vortex filaments were presented. In those cases, the rotational part of the fluid was characterized by a very small aspect parameter. Slender regions of constant vorticity—vortex filaments—were approximated by a single strand of overlapping vortex blobs and the filaments moved as a unit, without the possibility of variations in the structure of their cores. Advancing toward flows which are more representative of Earth's fluid core, this chapter deals with contained rotating flows which have an aspect parameter of order unity and are capable of supporting internal motion and vorticity variations. Although the core of an individual filament constitutes a rotating flow, the term will now be used to refer to bulk regions of fluid such as those contained by spherical or cylindrical boundaries. A fluid undergoing steady, solid body rotation has throughout it a constant vorticity equal to twice the rotation rate. Hence, the density of vortex lines, which are parallel to the rotation vector, is uniform. Taking advantage of the equivalence of a discretized vortex line and a thin filament, a large region of fluid rotating as a solid body can be modelled using a uniform arrangement of parallel, overlapping filaments. After establishing this approximation for the rotating fluid, computational results can once again be compared with theory. There are analytical solutions predicting wave phenomena in closed cylinders. In 63 Chapter 5. ROTATING CYLINDER 64 addition, there have been many experimental studies verifying the theoretical results [McEwan, 1970]. Hence, this chapter focuses on a rotating fluid in a cylindrical container and will follow an approach similar to that of the previous examples. The cylindrical container is denned by two parallel planar boundaries with a circular cylindrical boundary extending between them. The axis of the cylindrical surface is normal to the planar surfaces. The method of images is readily applied to the planar part of the boundary by extending the parallel filaments outside the ends of the container. In the case of the cylindrical part of the boundary, the method of images has no simple analogue. However, the physical nature of standing waves in a rotating cylinder allows the application of a slightly different approximate method using exterior fluid particles to enforce boundary conditions on the sides of the container. A cartesian coordinate frame is used to describe the fluid motion, where the cylindrical axis lies along the z axis and the end planes are parallel to xy. The rotation vector is taken as Clz. The aspect parameter of the cylinder is the ratio of its radius a to height d. This ratio is expressed as The value of the aspect parameter adopted for the numerical comparisons in this study is A = 0.357. This value was selected because it is close to one (A = 0.4) for which experimental results are available. 5.2 Inertial Waves It has long been recognized that the effects of rotation on fluid motion are quite remarkable. In particular, in a rotating reference frame, the 'fictitious' Coriolis force can act to restore perturbations to particle positions. It is this restoring effect that allows for wave Chapter 5. ROTATING CYLINDER 65 like motion superimposed on a steady rotation. In an inertial frame these waves may be viewed as the result of interaction between vortex lines which are slightly perturbed. From a reference frame rotating with constant angular velocity f l , the Navier-Stokes equations are du + u - V u + 2n ot — xu 1 o = — V P + uV u. p , (5.2) The term 2fl x tt is the Coriolis acceleration, which must be imposed when fluid motions are referenced to a rotating frame. In addition, centrifugal accelerations may be expressed as the gradient of a scalar potential and it is possible to consider their effects as a contribution to the pressure [Greenspan, 1968]. Hence, the reduced pressure is given by P = p-\p{nxx)\ (5.3) where the second term on the RHS represents the contribution of centrifugal force. In flows with small Ekman numbers, viscous effects are very small compared to Coriolis effects throughout most of the fluid. A useful approximation to the flow is obtained by omitting the viscous effects, which is equivalent to setting the Ekman number to zero. Of more importance in this study is the Rossby number, given by —5L< ( 5 ' 4 ) where U, il and L are representative scale factors of velocity, rotation and length. When e is very much less than one, Coriolis accelerations dominate over the non-linear term tt • Vtt, so that a linear theory applies and the potential for simple wave motion exists. Neglecting viscous and non-linear effects, (5.2) becomes du 1 — + 2ft x tt = - - V P . (5.5) Chapter 5. ROTATING CYLINDER 66 A qualitative examination of the physics underlying linear, inviscid rotating fluids gives insight into the mechanism which permits oscillatory motion within the flow. In addition, frequencies of wave-like solutions to (5.5) can be found mathematically for simple container geometries. The complementary physical and analytical approaches are presented with authoritative detail by Batchelor [1967] and Greenspan [1968], respectively, and are summarized here. By definition, the Coriolis force acts perpendicular to both the rotation and velocity vectors. Hence, only the velocity component lying in a horizontal plane contributes to the Coriolis force and thus influences the horizontal momentum balance. Periodic motion is possible because the Coriolis force opposes any velocity perturbation which results in a non-zero divergence in a horizontal plane [Batchelor, 1967]. Consider a fluid motion which causes an area enclosed by a material curve to increase. A change in velocity around the material curve, which is expanding as a result of the positive divergence, is necessary to maintain a constant circulation in an inertial frame, in accordance with Kelvin's circulation theorem (§ 3.1). In a rotating frame, this corresponds to a non-zero velocity in a direction opposite the rotation rate. The resulting Coriolis force is directed inward, opposing the initial outward disturbance. A particle under the influence of this force will follow a path, whose projection on a horizontal plane is a circle, back to its original position. A small perturbation to a uniformly rotating fluid, such as the one described above, results in a free motion, which may be represented as a superposition of normal modes, each with a periodic time dependence. These normal modes, known as inertial waves, are due to the interaction between Coriolis and inertial effects. Solutions for the inertial waves are sought in the form [Greenspan, 1968] u =Ue (5.6) Chapter 5. ROTATING 67 CYLINDER P = Ve . i<rt (5.7) Substitution of (5.6) and (5.7) into the momentum equation (5.5) and the conservation of mass V • U = 0 yields an eigenvalue problem for the frequency cr. The values of <r depend on the shape of the surface on which the boundary condition U • n = 0 is imposed. Solutions for contained fluids exist when jcrj < 2ft. Hence, inertial waves exist with angular frequencies, as observed from a rotating frame, which are less than twice the background rotation rate of the fluid. Inertial waves are periodic fluid motions superimposed on a steady, rigid body rotation. In closed containers axisymmetric modes of oscillation appear as standing waves. Further, for cylindrical geometry, analytical solutions exist for the frequency eigenvalues cr. The first step in applying the vorticity method is to generate a cylindrical region of constant vorticity representing a steady rotation. After creating such a flow, an appropriate disturbance is introduced, thus setting up inertial waves. These two steps are presented in the following two sections. 5.3 Initialization A steady rotating fluid is characterized by a uniform distribution of vortex lines. Consequently, a reasonable approach to modelling this flow is to set up a uniform distribution of discretized vortex filaments. A further requirement is that blobs in horizontal planes overlap in the same way that blobs overlap along the filaments. This permits a smooth vorticity distribution in the horizontal and vertical directions. Taking advantage of the grid-free nature of the vorticity method, the most practical option is arranging the set of filaments in concentric circles, thus forming 'vortex shells'. This allows good control over radial resolution and the aspect parameter A . Further, for systems with small numbers of particles, a plan view of the cylinder will not appear ragged and squarish, as it would Chapter 5. ROTATING 68 CYLINDER w i t h a rectangular arrangement of filaments. To maintain a particle spacing of approximately h, b o t h horizontally and vertically, the positioning of particles was determined i n three steps. Firstly, individual vortex filaments were constructed from a single array of blobs, as was done previously i n § 4.3. N e x t , the filaments were placed i n shells w i t h the radius of the j t h shell given by r ' = ( j - \ ) H - ( 5 > 8 ) Lastly, w i t h i n each shell, the filaments were equally spaced around the circumference, such that the distance between the filaments was not greater than h. A typical planar cross section of the particle arrangement is presented i n Figure 5.1. T h e particle positions may then be viewed as a vertical stacking of the cross sections shown i n Figure 5.1. W h e n initializing a single filament, it was previously found that abrupt truncation led to an oscillatory solution for the strengths. T h e problem was resolved by windowing the vorticity along the filament. T h e same windowing procedure was applied to the ends of the cylinder during initialization. However, a similar problem arose i n the radial direction. T h e resulting aliasing introduced a similar pattern of oscillatory vorticity strengths. T h e strengths w i t h i n some shells were consistently negative, interspersed between shells of positive strength. Physically, this corresponds to a t h i n shell-like region of fluid rotating i n the opposite direction to the surrounding fluid. A g a i n , this contradicts physical intuition, although the resulting vorticity field was roughly uniform. B y applying weighting to each successive shell, ramping the vorticity to zero outside the region of interest, the radial solution for the strengths became less oscillatory. Incorporating both axial and radial windowing into the initialization process created another, more prohibitive, problem. Including many additional exterior particles i n simultaneous axial and radial windowing produced an initialization m a t r i x which was too Chapter 5. ROTATING CYLINDER 69 2.0 o 1.0 ° ° o o o ° o ° ° ° o ° o o o o o o o o o ° o ° o ° o ° o ° o o o o o o o o o o o o o o o o o o o u y o.o o o O o O o -1.0 -2.0 -2.0 o -1.0 o O ° O o o o ° o O o o 0 o 0.0 O O ° _ 0 o 1.0 2.0 X Figure 5.1: Plan view of filament arrangement. 81 filaments arranged in five shells with radii given by (5.8). large to solve practically. To illustrate this, consider a cylinder constructed from filaments, each containing ten particles. Thesefilamentsare arranged into five shells, such that the aspect ratio is approximately one half. The number of interior particles for this case is more than 800. Including image particles and sufficiently long window ramps beyond the planar boundaries, increases this number to many thousands. Further, since the number of particles in each shell increases with shell radius, adding ten or twenty additional shells for radial windowing makes the total number of particles required in the Chapter 5. ROTATING CYLINDER 70 initialization O10 . Unfortunately, this size of system yields a sparse version of the ini4 tialization matrix ^ which is too large for practical use. Decomposing the initialization into separate radial and axial steps reduces the computational overhead while producing a reasonable vorticity distribution. It has been established that the blobs comprising a realistic straight filament have constant strengths. Thus, the problem of initializing the rotating cylinder reduces to finding the strengths of a single layer of particles, arranged in concentric circles and constrained to a specified radial distribution of vorticity. This 'disk' of particles represents a horizontal plane in the fluid. A vertical stacking of identical disks will then be consistent with the properties of singlefilamentsand produce a uniform vorticity profile. Given a number offilamentshells N initialized as above, the resulting radial distriba ution of vorticity decreases to zero, much like the function tp , just beyond the outermost s shell. A typical radial distribution of vorticity is shown in Figure 5.2, for N = 5. The a vorticity profile behaves similarly for all N , but is discussed here for the specific case t illustrated in Figure 5.2. The vorticity, and hence, rotation rate is roughly constant over the inner three shells. The vorticity begins decreasing inside the fourth shell, causing a significant gradient, with respect to radius, in the vorticity. The slope of the vorticity profile steepens as it drops to zero. The fifth, outermost shell is located where the vorticity is about half its maximum value. Hence, for a uniform vorticity in the cylinder interior, at least two additional exterior shells must be included. These exterior shells may be regarded as images in the sense that they can be manipulated to approximate the source of a potential flow. This point will be discussed further in the following sections. The radius and height of the cylinder were chosen such that A ~ 0.4. Two further considerations must be made when choosing the parameters a and d relative to the blob size 6. First, the cylinder ends must be sufficiently separated to minimize velocity contributions from the image particles outside the far boundaries. This corresponds to a Chapter 5. ROTATING 71 CYLINDER 6.0 " ' 0.0 1.0 2.0 3.0 Radius r Figure 5.2: Vorticity and azimuthal velocity profiles. Profiles correspond to the five shell case, with modified r and r (see § 5.5). The circles give the radial locations of the shells. 4 5 minimum of five or six particles per filament between the planar boundaries. Further, the number of shells must be chosen to reflect a trade-off between radial resolution and efficiency. Efficiency decreases dramatically when there are more thanfiveshells. Considering these factors, a standard computation involved three shells and seven particles per filament. Hence, the cylinder interior was discretized into 210 fluid particles. Including the two additional exterior filament shells and images across the end boundaries for all interior particles, the total number of particles in a standard computation was 1701. Chapter 5. ROTATING 72 CYLINDER The flattest vorticity profile within the first three shells was generated by beginning the cosine ramp with the first shell and including a total of twenty shells. Thus, the initialization of a single disk involved a total of 1266 particles. After the disk was properly initialized, the fifteen outermost shells were discarded and the disks were stacked. The resulting vorticity and azimuthal velocity profiles due to the five remaining shells are shown in Figure 5.2. In addition, the locations of the shells are shown. The variations in vorticity as a function of radius correspond to changes in the fluid's rotation rate. The effective, or average, rotation rate fi can be found from the location and value e of the azimuthal velocity at the radius corresponding to the third shell. In a closed, uniformly rotating cylinder of radius a, the velocity structure is given by the first part of (3.6), with the substitution r = a. The azimuthal velocity can be rearranged to find e the average rotation rate as _ e n uJa) = -^r- (5-9) This quantity will allow comparison with theoretical results obtained in terms of the rotation rate. With this arrangement of fluid particles, those making up the third shell coincide with the solid cylindrical boundary. Imposing the boundary condition u • n = 0 implies these particles must remain radially stationary. This is the final stage in setting up the contained rotating cylinder model and will be addressed in § 5.5. 5.4 Excitation of Inertial Waves Solutions to the eigenvalue problem described in § 5.2 can be found by separation of variables. The solutions for a closed cylindrical region of unit height and radius A [Greenspan, 1968] are Chapter 5. ROTATING CYLINDER V mnl (r, 73 0, z) = Jm (^^) cos(Trmz) e*", (5.10) where the eigenvalues are given by ^ = { + J%?)~'• 2 l (6U) The integers m, n and I characterize the axial, radial and azimuthal features of the modes, respectively. 7 j is the nth positive solution to a transcendental equation involving the mn Bessel function and its first derivative. The general solution for the flow is a linear combination of each mode given by (5.10). Hence, there are an infinite number of modes which can be excited in the cylindrical container. Indeed, an arbitrary perturbation to particle positions will excite any number of these modes. For unambiguous results, it is desirable to perturb the fluid such that few modes are excited. Both experimentally and theoretically it is simplest to consider modes where 1 = 0 and the values of m and 7 j mn are small. When 1 = 0, the values 7 m n 0 are solutions to the equation «7i(7 mn0 ) = 0, and henceforth will be referred to as ~f . In this case, the azimuthal velocity remains constant n in time and the fluid particles, as seen from a rotating frame, have non-zero axial and radial velocity components. For the numerical model, it is required to find an initial displacement field consistent with a single inertial mode. Rather than utilizing (5.10) to find an appropriate form of displacement, it is more straight forward to manipulate the stream function, whose derivatives are related to the fluid velocity components u and u . Batchelor [1967] gives r z the stream function for inertial modes of a closed cylinder with radius a, height d and azimuthal number I = 0 as i> = 2Ar J (3 r) sin (<r t) sin (kz) , x n n (5.12) Chapter 5. ROTATING 74 CYLINDER where k is the axial wavenumber and the constant 3 is related to 7 (see discussion n n below). Integrating the velocity components u and u , as found from the stream function r z (5.12), yields the axial and radial displacement perturbations. Referenced to a rotating frame, they are given by r{t) = B J (6 r) cos (<r t) cos (kz) , (5.13) z(t) = B J (3 r) cos (<r t) sin (kz) , (5.14) a n Q n n n where B = —2A/a and A is a small arbitrary amplitude. The initial displacement n disturbance is found by setting t = 0. These displacements represent a standing wave much like the one described in § 4.3. The radial perturbation, as a function of z, is the same as in the single filament case, except here, its amplitude varies with radius according to the Bessel function Ji(0 r). n Boundary conditions require that fluid particles at boundaries have zero normal velocity. Correspondingly, the z component of the displacement is always zero over horizontal boundaries placed at z = 0 and z = d, provided k = mir/d. The condition that the radial component of the perturbation velocities is identically zero at r = a requires t7x(/3na) = 0. Hence, 3 is a radial 'wave number', which can be found by B = -f /a, where 7 is the n n n n nth positive solution to J (r) = 0. Equivalently, the function J^r) is scaled by B such x n that it is zero at a specified radius r = a. Incrementing n is equivalent to choosing successively larger roots coinciding with the boundary, while the perturbation in r goes through more zeros within the cylinder. Similar to that of a sinusoidally varying perturbation, the radial mode number n describes the number of radial 'wavelengths' in the disturbance between r = 0 and r = a. Chapter 5. ROTATING CYLINDER 75 Imposition of the constraints, k = mir/d and J (/3 o) = 0, leads directly to the x n expression for the frequencies of oscillation [Batchelor, 1967]. They are ^VnnO = 20 f 1 + - 4 ^ ) ' . ( ) 515 where (2 is the background rotation rate of the fluid. Normalizing the frequency by the rotation rate and substituting the aspect ratio (5.1), yields an expression which agrees with the eigenvalues of Equation 5.11. Solutions for the displacements r(t) and z(t) in (5.13) and (5.14) remain valid even if the fluid is considered infinite in the horizontal direction. In this case, the roots of Ji(ft r) n represent possible boundary locations r = a, where the radial velocity is identically zero. The positions and strengths of particles in the region r > a can be used to evaluate the potential flow required to satisfy the boundary conditions on the radial velocity at r = a. Hence, it is possible to approximate the potential flow source by establishing a finite number of 'image' blobs outside the cylinder that reproduce the theoretical estimates of vorticity in r > a. These additional particles represent an extended region of the fluid, which behaves according to (5.12). 5.5 Numerical Modelling The first three shells of fluid particles correspond to the interior of the cylinder, while the outer two serve two purposes. The fourth and fifth shells ensure a uniform vorticity in the interior and provide the approximate source for the potential flow required to satisfy the boundary conditions at r , the radius of the third shell. By noting that 3 (5.16) r 7i the radii of the two exterior shells can be modified slightly such that they form one 3 Chapter 5. ROTATING CYLINDER 76 additional radial antinode. In other words, 7*5 can be increased to coincide with the second root of Ji(/3 r), after which r is adjusted to lie midway between its neighbours. n 4 The relative changes in the two radii are both less than 2%, and hence, the ratio of h/6 in this region does not change significantly. With this exception, the cylindrical region will be initialized according to the methods discussed in § 5.3. After the inertial oscillations are excited with a single mode perturbation given by (5.13) and (5.14), the system will be advanced in time. Frequencies will be computed by analysis of the particle positions as a function of time and compared to those given by the above theory. The focus will be on the modes (m, n, I) = (1,1,0), (2,1,0) and (3,1,0). After several trial runs it was found that an instability developed in the horizontal displacements of the particle positions after several rotations of the cylinder. This instability resulted in a change in the radii of the third and fifth shells, which, according to (5.13), must remain constant in time. Possible causes of the unexpected behaviour are discussed in detail in § 5.5.3 and § 6.4. In an effort to minimize the effects of the instability and reproduce a realistic rotating cylinder, the third and fifth shells were forcibly constrained to their theoretical positions by setting their horizontal position derivatives to zero throughout the integration. 5.5.1 Estimating Frequencies of Oscillation The integration of (2.21) and (2.22) results in a set of time series for the positions and strengths of all interior particles, referenced to a cartesian inertial frame. The frequencies of the inertial waves given by (5.15) or (5.11) are referenced to a rotating frame. The need for transforming particle positions or frequencies between the two reference frames is eliminated by noting that, in the case where / = 0, the wave is rotationally invariant. Hence, the periodic variation of a particle's radius in a rotating frame is equal to that in an inertial frame. Subtraction of a particle's initial radius from the absolute radius gives Chapter 5. ROTATING CYLINDER 77 the radial perturbation. Similarly, the axial component of position z is independent of the two reference frames. Frequencies of oscillation are estimated by analysis of the radial and axial perturbations over the time of integration. Prior to the integration of a particular mode of oscillation, a single horizontal disk of particles is selected such that the radial and axial perturbations have approximately equal amplitude at the height of the disk. Since the perturbations vary with radius, it is of interest to observe the fluid's behaviour within each shell. In order to minimize contributions from small scale, large wavenumber oscillations, particle positions and strengths are averaged over each shell. The resulting time series for the averaged radial perturbations and axial positions of each shell will be denoted by r*(i) and z (t), where i denotes the shell number. It should be noted x that the resulting time series for average radial perturbations in a given shell did not differ noticeably from the radial perturbation for a single particle in the same shell. This implies that the perturbation remained axisymmetric as time advanced. In order to find the dominantfrequencycomponents, the time series r*(t) and for modes m = 1 through m — 3 were modelled as an autoregressive moving average process. This is a signal processing technique which is useful for estimating the frequencies of sinusoidal components in noisy data. The frequency estimates are achieved by analyzing the information contained in a time series' autocorrelation matrix. The smallest eigenvalue of the autocorrelation matrix and its corresponding eigenvector contain valuable information about thefrequencycontent of the time series. In particular, a polynomial may be formed using the coefficients of this eigenvector. The roots of the polynomial lie on the unit circle in the complex plane at angles of 2irfjAt, where At is the sampling interval of the time series and f jt j = 1 , P , are the sinusoidal frequency estimates. Chapter 5. ROTATING CYLINDER 5.5.2 78 Results The time series r*(t) and i = 1 , 5 for three inertial modes are shown in Figures 5.3, 5.4 and 5.5. In addition, each figure shows the time series of the vertical component of the average vorticity strength T),(t) for each of the five filament shells. In the interior shells, there is a dramatic increase in the amplitudes of oscillations at approximately t = 10. Prior to the increase, the amplitudes are the order of the specified perturbations. As the integration proceeds the oscillations grow in amplitude reaching, in some cases, as much as ten times their original value. However, the amplitudes for the exterior shells (i = 4,5) begin increasing early in the evolution of the flow. This feature and others are discussed in detail in the following paragraphs. The plots of the radial perturbations can be viewed as changes in radius of the filament shells over time. The curves corresponding to the third and fifth shells {i = 3,5) are always zero due to the constraint that all particles in these shells remain radially stationary. The curves for the inner two shells indicate that the average shell radii, and the relative positions of the shells, change significantly towards the end of the integration. For example, in the mode m = 2 case, the average radius of the second shell reaches a maximum (at t « 16), while the perturbations to the two adjacent shells are much smaller. A radial perturbation of this size has consequences in terms of the vorticity. Due to the increased overlap of blobs in the second and third shells, there will be an accompanying increase in the total vorticity in that region. This crowding of filaments in adjacent shells can also be viewed as an increase in the vortex line density. The average radius of the fourth shell in Figure 5.3 becomes larger than that of the the fifth shell at t « 20 (the shells are initially ~ 0.25 apart). This indicates that some of the particles in shell four have 'escaped' the extended cylinder. Clearly, there is an instability developing in the particle positions as the integration proceeds. The amplitude of the radial Chapter 5. ROTATING CYLINDER 79 Figure 5.3: Oscillations in positions and strengths for mode (1,1,0). From top to bottom, radial perturbation, axial position and vertical component of strength, i indicates shell. Chapter 5. ROTATING CYLINDER 80 Figure 5.4: Oscillations in positions and strengths for mode (2,1,0). From top to bottom, radial perturbation, axial position and vertical component of strength, i indicates shell. Chapter 5. ROTATING CYLINDER 0.0 0.8 81 5.0 1 10.0 1 1 15.0 i 1 i 20.0 l; i lj li 1! —1 r\t) 0.4 . =1 i=2 - - - i =3 i=4 — - i=5 1 iij; |'i 0.0 i -0.2 0.8 i i ~i r 1 i i 1 j; 1: 1 r 0.6 z\t) 04 0.2 0.0 J L "i 1 J 1 • 1 1 1 L 1 1 r 0.15 0.10 r>) 0.05 0.00 -0.05 -0.10 J 0.0 5.0 L 10.0 J i_ 15.0 20.0 Time? Figure 5.5: Oscillations in positions and strengths for mode (3,1,0). From top to bottom, radial perturbation, axial position and vertical component of strength, i indicates shell. Chapter 5. ROTATING CYLINDER 82 perturbation in the fourth shell begins increasing well before the others. The instability likely begins in this shell (see § 5.5.3) and then propagates toward the inner shells. The plots of the axial positions show similar behaviour. The axial positions of the particles in shells three and five are not constrained, and hence, have non-zero perturbations. Again, shell four is suspicious in its early increase in amplitude. In general, the amplitudes of the axial oscillations appear to be increasing in an exponential manner. In addition to examining closely the behaviour of the particle positions, the behaviour of the vorticity strengths is as important. In particular, the vertical component of strength is relevant because it is directly related to the rotation rate of the fluid. In the first half of the interval, the vertical strengths averaged over each shell remain roughly constant with some small amplitude oscillations. This is consistent with a fluid rotating uniformly. However, as the integration progresses, the strengths become more oscillatory and erratic, giving an indication that the fluid is no longer rotating uniformly. That is, at a particular time, there is significant variation of the strengths, and hence, vorticity as a function of radius. In the mode m = 2 case, the strength of the fifth shell decreases rapidly, taking on negative values near the end of the interval. With increasing axial mode number the behaviour of the average shell strengths becomes more erratic. In some extreme cases, the physical behaviour of the fluid particles under observation is unintuitive. In general, there is evidence that an instability causes a rapid increase in the amplitudes of oscillation. However, the interior shells (i = 1,2) are of primary interest and over the time interval shown in Figures 5.3, 5.4 and 5.5 the interior flow behaviour is reasonable enough to extract frequency estimates. The proposed method of frequency estimation requires the specification of P, the number of harmonics in the signal. Ideally, the fluid particle motion was to contain only one sinusoidal component with a frequency given by (5.15). However, the time series r'(£) and z*(t), i = 1,2, examined above appear to contain several sinusoidal components. To Chapter 5. ROTATING 83 CYLINDER explore the possibility of multiple inertial modes coexisting in the cylinder, the frequency estimation was carried out using the values P = 1,2,3,4, The resulting frequency estimates are given in Table 5.1. For each mode, a similar analysis of each time series was performed and all had identical frequency content. Mode 1, cr 1 2 0.088544 0.59695 0.87806 110 = 0.37405 4 3 0.42908 0.33582 0.62917 0.48994 0.84584 0.98507 1.1392 Mode 2, <r = 0.67240 1 2 3 4 0.086724 0.56926 0.40933 0.62833 0.83642 0.59942 0.46662 0.80625 0.93886 1.0852 Mode 3, cr = 0.87862 1 2 3 4 0.088895 0.59932 0.43078 0.33716 0.88155 0.63167 0.49188 0.84920 0.98898 1.1437 210 310 Table 5.1: Frequency estimates for three axisymmetric modes. Estimates are given for P = 1,2,3,4. The analysis for P = 1,2,3 did not resolve frequencies similar to those predicted by (5.15). In each of these cases, none of the frequencies changes significantly with increasing mode number. For P = 4, there is an interesting change in the behaviour of the frequencies for the different axial modes. Three of the frequencies change with m in a manner similar to the P = 1,2,3 cases. For modes one and two, the frequency estimate at the top of the column is very similar to the theoretical frequencies. However, the first Chapter 5. ROTATING CYLINDER 84 frequency listed for mode three does not continue the trend. In fact, it is quite similar to the frequency estimated for mode one. Regarding P = 4 as the correct choice for the number of harmonics in the time series, at least one of the three frequencies which remain roughly constant with m can be connected to the motion shown in Figures 5.3, 5.4 and 5.5. There is a small amplitude oscillation evident in most of the averaged particle positions in the first half of the time interval. This oscillation has a period of about two and corresponds to the frequencies which are sUghtly less than one half. The two higher frequency components, each with a period of approximately one, are difficult to distinguish in the figures. The results show there are three similar modes present in the cylinder for each specific perturbation. However, the origin of these additional modes is not certain. It is possible that they arise from the interior velocity contributions from the unstable fourth shell or from the error associated with the finite approximation of the exterior flow. In the same way that the frequency content of standing waves, on strings for example, includes overtones, the fluid may be supporting modes with a higher spatial complexity. 5.5.3 The Vorticity Profile and Stability of Radial Perturbations The results presented in the previous section show an instability developing in the radial component of the perturbation. Forflowswith circular pathlines, such as a uniformly rotating fluid, conditions can be established to indicate when a small radial perturbation is stable or not. In this sense, stability implies that the amplitude of the perturbation will not grow in time and the fluid will support decaying or constant amplitude oscillatory motion. Consider a particle in a uniformly rotating fluid, which travels on a circular pathline. The particle is given a small radial perturbation while its angular momentum per unit mass is conserved. The Rayleigh stability criterion requires a sufficient pressure gradient Chapter 5. ROTATING CYLINDER 85 at the particle's new position to maintain centripetal acceleration associated with the circular motion [Vanyo, 1993]. When this criterion is met, the particle returns to its original position. If the particle overshoots its original position and the subsequent radial perturbation is also stable, there will be an oscillatory motion, or inertial wave. The Rayleigh criterion states that a radial perturbation to a circular pathline is stable if ( r , uY > ( r , u ) 2 f where r and t { , (5.17) are the initial and final radii of the particle's position, «^ is its initial linear velocity and Uj is the ambient fluid velocity at rj [Vanyo, 1993]. This result implies that a fluid in rigid body rotation, where Cl is constant, is always stable to radial perturbations. However, in the example presented above, the vorticity and rotation rate in the cylinder and, in particular, the exterior shells are not constant. A small perturbation in the rotation rate accompanying the radial perturbation results in Rayleigh instability. The linear velocity, or azimuthal velocity component, in a uniformly rotating fluid is u = Clr. With a non-uniform rotation rate, the ambient velocity { at radius r + dr is i u = (Cl + dCl) (r + dr) , f { (5.18) where dCl represents a finite change in the rotation rate over dr. The Rayleigh criterion (5.17) then becomes 2Cldr + TidCl > 0. (5.19) Examining the different cases of positive or negative values for dr and dCl shows that a varying vorticity profile is inherently unstable to radial perturbations. Chapter 5. ROTATING CYLINDER 86 Figure 5.2 shows that the fluid particles in the fourth and fifth shells, when given a radial perturbation, will experience a change in rotation rate. The particles in the fifth shell are not initially perturbed. However, without explicitly constraining them by equating their derivatives to zero, they become unstable. This is likely due to the influence of their neighbours in the fourth shell, which are unconstrained and subject to radial perturbations within a varying vorticity field. After several rotations of the cylinder the instability propagates through all shells and the uniform rotation breaks down. Chapter 6 DISCUSSION A N D C O N C L U S I O N S 6.1 The Vorticity Method The vorticity method algorithm provides the basis for the computations presented in the previous two chapters. In part, convergence and accuracy of the method depend on the number of fluid particles in the discretization of the vorticity field. A large number of particles improves accuracy, but adversely affects the efficiency. The most time consuming operation of the algorithm is the evaluation of the derivatives in the time stepping (see § 2.4.2). There is a summation over all N particles for each derivative at a single particle, so time stepping is approximately an N operation, where N is the total 2 number of particles including images. Hence, it is desirable to keep N as small as possible, while still maintaining convergent results. Each of the test cases involving a single filament produced accurate results while maintaining efficiency with a minimal number of particles. Convergence in some cases was oscillatory (see Figure 4.10), introducing uncertainty into the evaluation of relative errors. However, in general, the errors were small and predictable. For systems with a larger number of particles, such as the rotating cylinder, it became more difficult to quickly check for convergence. Studying the convergence in a comprehensive manner, as with the examples of Chapter 4, became very time consuming. Since large systems are the primary application of this method, it is important to examine carefully the algorithm's behaviour with small systems. The insight gained from such 87 Chapter 6. DISCUSSION AND CONCLUSIONS 88 problems is of great value in studying larger problems. For example, the accuracy of the numerical calculation of standing wavefrequencieson a filament was affected by the finite length of the exterior part of the filament. Those effects and their relation to the infinite filament, as described by Kelvin's theory, were studied in detail in § 4.3. These results (see § 6.5 for a summary) can be applied to larger physical models having similar properties. The extension of methods and results to large systems enforces the need to make a proper connection between computational methods, such as the method of images, and the physical interpretation of the problem. The program's run time is surprisingly small when the total number of particles is several hundred. That is, N is O10 . Further, with this few particles, the spatial 2 initialization of the vorticity strengths becomes trivial when solving the system (2.29) using the conjugate gradient method with a sparse matrix. This method of initialization and the quick time stepping illustrated the power of the method for modelling flows characterized by small patches of vorticity. The example of the 'leap frogging' vortex rings in § 4.2 shows that the efficiency of the method is maintained when the physical problem increases in complexity. The three quantitative test cases—ring velocities, single filament waves and inertial oscillations in a cylinder—were compared against theoretical approximations. In each case, the theoretical approximation is linearized. The computational method, implicitly includes non-linear effects, but it reproduces the linear solutions when the amplitude of the motion is small. 6.2 Initialization For large N, the conjugate gradient method was used with a sparse matrix approximation to facilitate the initialization step. In fact, as N approaches C?10 , the size of 3 Chapter 6. DISCUSSION AND CONCLUSIONS 89 the full linear system prohibits the direct solution by any other method, such as LU decomposition. Further, the speed of initialization by the conjugate gradient method for smaller numbers of particles is advantageous for frequent changes in parameters affecting geometry, number of particles or blob size. Experience gained from the initialization of a single filament led to a quicker method for the cylinder. Knowing that the strengths T for a singlefilamenthave approximately 0 equal value for each particle, the solution for the cylindrical arrangement of filaments was decomposed into two steps. One was the construction of a single horizontal plane and the other involved stacking identical 'disks' of particles. This decomposition further reduced initialization time and facilitated the generation of a more reasonable vorticity profile because more particles could be used in each independent step. However, this method produced a profile which varied radially, especially in the exterior shells, creating a potential for instabilities throughout the flow. Some attention was focused on creating a uniform vorticity field within model domains. In all cases, whether planar or cylindrical boundaries, extra particles were required outside the boundaries to produce constant vorticity near the boundaries. Other methods of achieving a flatter radial vorticity profile were briefly explored. The most promising alternative involved taking advantage of the constant vortex line density in a uniformly rotating fluid. The average strength of each shell was scaled according to the total density offilamentswithin that shell's radius. The shell strengths were then found by solving a linear system of equations. Another improvement which might be made to the existing computer code is separating the initialization from the integration process. The capability of storing and reusing solutions to (2.29) would be a benefit when considering large systems such as rotating fluids. In such a case, N is likely to be O10 . This approach would make the initializing 3 Chapter 6. DISCUSSION AND CONCLUSIONS 90 of the cylinder as a whole more practical because the result could be reused in subsequent computations. It is also possible that the resulting vorticity profile would be more uniform. In addition, saving a final solution—after a period of time integration—would allow the restarting of a previous run. 6.3 Filament Results The predicted error in the theoretical approximations for the ring velocities and wave frequencies are of the same order. The relative errors between theory and computation for these two test cases are roughly O 1 0 , or less. These results show that the algorithm -2 produced results of similar accuracy forfilamentmodels of different geometries. These results were produced with a particle spacing h and blob size S, which result in the ratio h/S=\. With the inclusion of planar boundaries, the accuracy of results also depended on the number of images. The discussion in § 4.3 shows that an increase in the number of images produced convergent values of frequency. The image particles were incorporated to satisfy boundary conditions, which has proven very effective for the planar boundaries in this example. Theoretical studies of the problem implicitly assume that the vortex filament extends indefinitely beyond the boundaries. Thus, the finite number of images introduced a small source of error into the above comparisons. The image across one planar boundary for a particle in the interior will also contribute a small velocity at the other boundaries of the fluid {e.g. the second planar boundary in the case of waves on a bounded filament). Consequently, a second set of images must be added on the opposite side of the bounded fluid to correct for the influence of the first set of images. To satisfy boundary conditions perfectly, an infinite number of images would be required. In practice, however, only a small number of images was required to reproduce theoretical Chapter 6. DISCUSSION AND CONCL USIONS 91 predictions. 6.4 Cylinder Results There is evidence that the physics of the inertial wave problem was captured in the contained cylinder computations. One of the frequencies estimated for the case of four harmonics (P = 4) resembles the theoretical predictions. However, this is true for only two of the three modes presented, m = 1 and m = 2. One problem in the frequency estimation is the existence of other modes in the computed solutions. The method which was used to estimate the frequencies from the particle motion requires prior knowledge of the number of sinusoids in the signal. The choice of P = 4 is supported by the agreement in frequency for the first two axial modes and the appearance of the oscillation with period ~ 2 in the time series. In each of the three cases presented in Chapter 5, the additional modes are similar in frequency. However, their origin and classification (ie. modal numbers (m,n, I)) is presently unknown. The estimation of sinusoidal frequencies in the cylinder example is based on several assumptions. Firstly, the time series r*(t) and z*(t) are assumed axisymmetric. The similarity of the averaged results and those for a single particle indicates that the initial perturbation remains axisymmetric. However, the increasing complexity of the flow due to an instability may alter the axisymmetry at later stages of the computation. Indeed, where the averaged results differed noticeably from the single particle results was near the end of the integration interval. Secondly, thefrequencyestimates are made by assuming the time series are comprised of equally spaced data points. The method of integration uses an adaptive step size control and the increasing amplitudes of oscillation led to large derivatives and smaller step sizes as time advanced. In each case, the step size decreased from approximately 0.4 to 0.25 over the period of integration. Hence, there is a small Chapter 6. DISCUSSION AND CONCLUSIONS 92 uncertainty in the frequency estimates due to the changing step size. A more important issue is the cause of the instability leading to growing oscillation amplitudes. Section 5.5.3 provides a possible explanation. However, the inaccurate treatment of boundary conditions on the cylindrical part of the boundary makes it difficult to determine if the instability arises strictly because of the variations in vorticity as a function of radius. In fact, the instability may arise from a combination of the two problems. Since the use of two exterior shells for the source of the potential flow is a simple approximation, the interior flow has an associated error. The radial instability arising from the vorticity gradient near the fourth shell causes the particles in that shell to execute a motion inconsistent with the theoretical predictions. This in turn leads to increased error in the interior potential flow, allowing the instability to propagate inward. Increasing the number of exterior shells may slow the propagation to the interior, but this comes at great expense to efficiency. Alternate methods of satisfying boundary conditions, especially on curved boundaries, will illuminate these problems. One such approach is described in Appendix A. 6.5 Boundary Conditions The method of images works well when modelling planar boundaries in the sense that normal velocities can be made arbitrarily small with relatively few images. In this context, the images represent an approximate source for the potential flow required to satisfy boundary conditions on the velocity. Errors in the construction of images leads to error in the potential flow and the boundary conditions will not be satisfied exactly. The method of images also ensures the requirement that the normal component of vorticity vanish on boundaries of the extended fluid region. In this case, the images continue the vorticity out of the fluid domain into the extended region where u> • n vanishes. Chapter 6. DISCUSSION AND CONCLUSIONS 93 A small uncertainty was introduced into the wave problems because the images failed to mimic symmetrical effects of an infinite filament in the theoretical treatment. These effects were found to be small when studying waves on a vortex filament because the boundaries were far enough apart that the images across one boundary had little influence on the velocity at the other boundary. In that example, a value such as N /N = ^ { resulted in normal velocities at the boundaries which were several orders of magnitude smaller than the tangential velocities. However, at this value of A^/iV the frequencies had not converged adequately (see, for example, Figure 4.8). In the cylinder example the planar boundaries are much closer together, increasing the velocity contribution the exterior particles have at the far boundaries and introducing potentially larger errors in the computed solution. In addition to satisfying the boundary conditions on velocity and vorticity, a small number of image vortices was required to maintain a constant vorticity near the boundaries. Saffman [1992] shows that it is possible, by using image line vortices, to satisfy normal velocity conditions on a cylindrical boundary. Here, the line vortices are straight and infinitely long. It is unclear, however, whether this method can be extended to filaments with some curvature or to cylindrical boundaries of finite length. For cylinders where the boundary curvature is small, a possible solution is the use of a locally planar approximation. A major drawback is that curvatures which are small on the scale of blobs sizes imply a very large number of fluid particles. The method of approximating the potential flow by including particles outside the cylindrical boundary relies on the fact that the interiorflowcan be considered to be a closed region within an infinite fluid. When explicit knowledge of the exterior flow is available, the exterior particles may be included in the derivative evaluations directly. However, in the most general modelling of rotating flows, image strengths and orientations will have to be determined with an additional calculation, such as the solution of a linear system of equations for the image Chapter 6. DISCUSSION AND CONCLUSIONS 94 strengths. The very large number of exterior particles, relative to interior particles, which were included in the cylinder computations suggests that there might be a more efficient method to deal with boundary conditions. An alternative approach to the boundary conditions is to solve directly for the potential flow. When required, velocities can be evaluated at the nodes of a grid on the boundary. As described in § 2.6.1, these velocities are applied as the boundary conditions in the solution of Laplace's equation. In using this method, it would be advantageous to choose a routine which is capable of solving for the potential at arbitrary points in the fluid, namely at particle positions. General efficiency is a concern, as well, since this operation must be carried out at every evaluation of the derivatives. Appendix A details a method which uses sheet vortices on the boundaries to satisfy the boundary conditions. This method uses M discrete boundary elements to represent a sheet vortex over the boundary, with each patch contributing to the velocity in the interior. The solution for the potential flow involves the inversion of a 3M x 3M matrix. Each boundary element then contributes to the derivatives for the positions and strengths of each interior particle. After performing the required matrix inversion, the evaluation of the derivatives involves N(N + M) operations. In terms of the cylinder example of Chapter 5, a comparable spatial resolution over the boundary would require M = 112. Further, incorporating this method into the cylinder example would eliminate over half of the required exterior fluid particles from the computation. Thus, the increase in computational overhead introduced by the vortex sheet method may not be terribly significant. 6.6 Resolution Spatial resolution of flow features depends on the density of fluid particles. Only those features, which are characterized by a length scale not smaller than a few h, can be Chapter 6. DISCUSSION AND CONCL USIONS 95 properly resolved. Hence, high resolution will require large number of particles and comes at the expense of efficiency. When modelling wave phenomena, higher order modes become increasingly difficult to model because more fluid particles are required. In analogy with signal processing, the maximum wavelength which can be resolved is 2h. That is, there must be a minimum of one particle between each node of a disturbance. For example, the axisymmetric mode (6,2,0) in the cylinder example of Chapter 5 is the highest order, both axially and radially, which can be resolved without increasing the number of interior particles. In the z dependence of the m — 6 mode, the standing wave consists of three complete wavelengths. With only seven particles along each vortex filament, it would not be possible to model modes with m > 6. Similarly, with only two radially perturbed interior shells, the highest order radial mode which can be properly resolved is the n = 2 mode. 6.7 Summary In general, the results presented in this thesis are very promising. The single filament results are in excellent agreement with theory and the algorithm can be easily adapted to model manyflowswith sparsely distributed vorticity. Incorporating the sheet vortex method into the algorithm may make it possible to modelflowswith different boundary shapes. Although the results for the rotating cylinder are not as spectacular as those for the single filaments, this work has provided the framework for improvements and advancements in the modelling of contained rotating flows. The two main issues which must be resolved are the possible instability of a varying vorticity profile and the inefficient or inaccurate treatment of boundary conditions. The effects of these two issues are currently difficult to distinguish because it is not clear whether the unstable behaviour is due to Chapter 6. DISCUSSION AND CONCLUSIONS 96 the vorticity variations alone or a combination of that plus an erroneous potential flow. In order to improve this algorithm and begin considering other geometries, such as spheres and spherical shells, these are the two issues of primary importance. Appendix A THE SHEET VORTEX METHOD The following analysis outlines a method of solving for the potential flow required to satisfy boundary conditions. Based on the concept of a sheet vortex, this method permits direct evaluation of the solution at each of the interior particle positions. A sheet vortex is a singular distribution of vorticity over a surface within a fluid. The vorticity is given by u> = KS(U) , (A.l) where K is the vector strength of the sheet, n is distance measured normal to the sheet and 8(n) is the Dirac delta function. Sheet vortex lines are everywhere parallel to K and he on the surface of the sheet vortex. Hence, K • n = 0, where n is the unit vector normal to the surface. A sheet vortex can be visualized as a surface across which there is a discontinuity in the tangential velocity component, while the normal velocity component is continuous. Where there are smooth distributions of vorticity adjacent to the sheet, the divergence of the strength K gives the jump in the normal component of vorticity across the sheet [Saffman, 1992]. Hence, placing an appropriate sheet vortex adjacent to solid boundaries can be used to satisfy all boundary conditions on velocity and vorticity, which must be imposed for inviscid flow. The no slip condition in a viscous fluid, requiring the component of velocity tangential to the boundary to be zero, can also be satisfied using a sheet vortex. In fact, a sheet vortex adjacent to solid boundaries is representative of an infinitely thin boundary layer. Henceforth, sheet vortices will be assumed to coincide 97 Appendix A. THE SHEET VORTEX 98 METHOD with the solid boundaries of a contained rotating fluid. The boundary condition to be imposed is « + « = •», (A.2) where u is the rotational velocity field due to the interior vorticity (given by (2.10)), u is the potential velocity due to the sheet vortex and v is the velocity of the boundary. Analogous to expression (2.10) the velocity due to a sheet vortex is u (as, t) = J V x G(as - as') K (a?', t) dx', (A.3) s where the integral is taken over S, the bounding surface [Saffman, 1992]. Provided the sheet strength K is chosen properly, this velocity field represents the potential flow required to satisfy the boundary condition (A.2). Finding the appropriate sheet strength requires an expression for the near field behaviour of (A.3). In the limit, as a point P in the interior approaches a point Q on the surface S, the velocity is given by Urn u(P) = -\K{Q) + q(Q) xn(Q) , (A.4) where n is the outward pointing, unit normal vector. The term q is given by the principle value integral q(Q) = J V x G (as - as') K (as', t) dx' (A.5) s* over the bounding surface, where S* indicates that the vorticity within an infinitesimal circle around Q is excluded from the integral [Saffman, 1992]. Expression (A.4) gives the velocity at the surface d u e to a coincident sheet vortex with strength K. Hence, setting the LHS of (A.4) equal to v — u and solving for K Appendix A. THE SHEET VORTEX METHOD 99 yields the required sheet strength. The irrotational velocity can then be evaluated at a given interior point according to (A.3). In a discretized system with a finite number of boundary elements,findingthe correct strength field K is accomplished by solving a set of linear equations. Discretizing (A.3) as in § 2.3, the velocity at an interior particle point x due to the x sheet vortex is given by M u' = X ; V x G ( s - as') K'dS', 3=1 i (A.6) - where the surface S has been divided into M elements of area dS . Similarly, after 3 substitution of (A.2) and (A.5), expression (A.4), for the velocity at the surface, becomes i t -« i i AT = --« xn i + i J V x G (a?*" - x ) K ' dS . 3=1, jfr 3 3 (A.7) 3 Hence, prescribing the velocities v and u , results in a linear system of equations for the % l strengths #c\ v In a manner similar to that in § 2.4.2, the system of equations (A.7) can be expressed in subscript notation as 1 4 - «i = - j d £ Wm » m « l + E " «'') *?' ^ klm Q~T^ 6 • (A.8) The evaluation of this system is facilitated by the introduction of the second order tensors K and L, which are defined by *2 = *ki JrG ( . « - a^) dS , 3 m (A.9) m Ti i 1 L\l = - 2 W m « m e (A. 10) Appendix A. THE SHEET VORTEX 100 METHOD Again, the superscripts make explicit the fact that the tensors depend on the field and/or source positions x and X*. Then, (A.8) becomes % M 4 - < (A.ll) £ = Lii>4+ Finally, after defining a third tensor K as % = 3 ki K = (A.12) else kl K the linear system for the strengths K\ reduces to M (A.13) i=i In matrix form, (A. 13) can be written as -u K 11 K V -u K 21 K V 1 1 2 v M 2 K 12 K 22 IM -MM -u M K 2M K K (A.14) M where the sub-matrices K* are given by } . 12 13 21 A 22 •^23 31 •^32 •^33 (A.15) The elements of the sub-matrices defined by (A.15) are given by (A.12). Hence, solving the linear system (A.14), which requires the inversion of a 3M x 3M matrix, yields the required sheet strengths K . The corresponding velocities at the interior particle positions j can be evaluated according to (A.6). Appendix A. THE SHEET VORTEX METHOD 101 Incorporating the sheet vortex method of satisfying boundary conditions into the algorithm is a straight forward modification to the ODEs (2.21) and (2.22). Given an initial, or updated, set of particle positions SB* and vorticity strengths T*, i = 1 , N , the sheet strengths K * , i = 1 , M , are found and the required derivatives are then evaluated according to _ d t = £ V x G (SB* - SB'J I" + £ V x G (SB*' - SB'') K dS j=l j=l j s N dt ^ V x G , (A. 16) K ^ S * (A.17) j Af { (SB*' - x) T + £) V x G 3 3 (SB* - SB') Note that the areas of the boundary elements dS* appear explicitly in both the modified derivative expressions (A.16) and (A.17), and the definition of K in the formulation of the linear system (A. 13). Hence, these quantities must be evaluated and stored for use throughout the integration process. To maintain a similar accuracy in both the interior and boundary contributions to the derivatives, the areas dS should be Oh , where h is % 2 the interior particle spacing. In other words, the spatial resolution on the boundaries must be approximately equal to that in the interior. The sheet vortex method satisfies the viscous boundary condition (A.2). Treating the flow near the boundaries this way is not inconsistent with an inviscid interior. In a real fluidflow,even one characterized by a very large Reynolds number, the tangential velocity component equals the velocity of the boundary, it., there is no slip. Hence, satisfying the boundary conditions this way will produce a more realistic interior flow. This approach, known as boundary layer theory, is common in modelling flow at high Reynolds numbers. In this case, the fluid is divided into two regions: an interior region where viscous effects are negligible; and a thin boundary layer where viscous effects dominate. Since there are large gradients in the tangential velocity over e, the thickness of a boundary layer, there Appendix A. THE SHEET VORTEX METHOD 102 is a corresponding increase in the component of vorticity tangential to the boundary and normal to the velocity. The boundary sheet vortex can then be redefined formally as a limit [Saffman, 1992]. In the limit, as the boundary layer thickness goes to zero, and the magnitude of vorticity u> goes to infinity in such a way that tu) —»• «, the vortex lines associated with the shear region become confined to a surface—the sheet vortex. REFERENCES 1. Anderson, C. and C. Greengard, On vortex methods, SIAM J. Numer. Anal, 22, 413-440, 1985. 2. Anderson, C, and C. Greengard, Vortex ring interaction by the vortex method, Vortex Methods, Lecture Notes in Mathematics, 1 3 6 0 , 23-36, 1988. 3. Batchelor, G.K., Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, 1967. 4. Beale, J.T. and A. Majda, High order accurate vortex methods with explicit velocity kernels, J. Comp. Phys., 5 8 , 188-208, 1985. 5. Greenspan, H.P., The Theory of Rotating Fluids, Cambridge University Press, Cambridge, 1968. 6. Griffiths, D.J., Introduction to Electrodynamics, Prentice Hall Incorporated, Englewood Cliffs, 1989. 7. Hama, F.R. and J. Nutant, Self induced velocity on a curved vortex, Phys. Fluids, 4, 28-42, 1961. 8. Knio, O.M. and A.F. Ghoniem, Numerical study of a three dimensional vortex method, J. Comp. Phys., 8 6 , 75-106, 1990. 9. Leonard, A., Computing three dimensional incompressible flows with vortex elements, Ann. Rev. Fluid Mech., 1 7 , 523-559, 1985. 103 REFERENCES 104 10. Malkus, W.V.R., An experimental study of global instabilities due to the tidal (elliptical) distortion of a rotating elastic cylinder, Geophys. Astrophys. Fluid Dynamics, 48, 123-134, 1989. 11. Mathews, P.M., B.A. Buffett, T.A. Herring and 1:1. Shapiro, Forced nutations of the Earth: Influence of inner core dynamics 1. Theory, J. Geophys. Res., 96, 8219-8242, 1991. 12. McEwan, A.D., Inertial oscillations in a rotating fluid cylinder, J. Fluid Mech., 40, 603-640, 1970. 13. Puckett, E.G., The random vortex method with vorticity creation: introduction and guide to parameter selection, Vortex Dynamics and Vortex Methods, Lectures in Applied Mathematics, 28, 567-584, 1991. 14. Saffman, P.G., Vortex Dynamics, Cambridge University Press, Cambridge, 1992. 15. Smylie, D.E., X. Jiang, B.J. Brennan and K. Sato, Numerical calculation of modes of oscillation of the Earth's core, Geophys. J. Int., 108, 465-490, 1992. 16. Thomson, W. (Lord Kelvin), Vibrations of a columnar vortex, Phil. Mag., 10, 155-168,1880. 17. Tritton, D.J., Physical Fluid Dynamics, Oxford University Press, New York, 1988. 18. Vanyo, J.P., Rotating Fluids in Engineering and Science, Butterworth-Heinemann, Stoneham, 1993. 19. Yamada, H. and T. Matsui, Prehminary study of mutual slip-through of a pair of vortices, Phys. Fluids, 21, 292-294, 1978.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The application of vorticity methods to rotating fluid...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The application of vorticity methods to rotating fluid flows McMillan, David G. 1996
pdf
Page Metadata
Item Metadata
Title | The application of vorticity methods to rotating fluid flows |
Creator |
McMillan, David G. |
Date Issued | 1996 |
Description | Studies of Earth's core dynamics often require a computational method which can account for non-linear effects and non-periodic time dependence. One important application involves the large scale, possibly chaotic, fluid motion associated with a past resonant tidal forcing of the free core nutation. This thesis explores the utility of vorticity methods for modelling contained rotating fluids, such as the Earth's core. Establishing this method in the general context of rotating fluids is the first step in the development of a non-linear, three dimensional, time dependent model of Earth's fluid core. Conventionally, vorticity methods have been applied to incompressible fluid flow in infinite domains with small, finite regions of non-zero vorticity. The capability of approximating solutions to non-linear, three dimensional, time dependent, inviscid fluid flow problems suggests that these algorithms may also be well suited to flows in finite domains under uniform rotation. Physically, the method under consideration involves grid-free tracking of fluid particles which carry distributions, or 'blobs,' of vorticity, each with a prescribed strength. The strengths and positions of the blobs determine the velocity field through the Biot-Savart law. The fluid particles, along with their vorticities, are advected by the velocity field. In addition, the vorticity of each blob is altered by the strain associated with the velocity field and must be updated at each time step. The method of solution reduces to a set of first order, ordinary differential equations in time for the Lagrangian displacement of the fluid particles and their corresponding vorticity strengths. The ODEs are advanced by standard integration routines. This work is a presentation of the development of the vorticity method algorithm for rotating fluids. In addition to the development of the algorithm, several new computational methods are described. The concept of effective vorticity is introduced as a means of relating known physical quantities to computational results. In order to satisfy inviscid boundary conditions in an efficient manner, approximate methods using image particles are developed and incorporated into two of the test cases. Quick initialization of vorticity fields is made possible by the formulation of an approximate linear system of equations. The modelling of a uniformly rotating fluid requires the initialization of a uniform vorticity field. This procedure is decomposed into two steps, increasing efficiency while creating a reasonable initial vorticity field. The implementation of the algorithm is verified by comparison with theoretical approximations. The first test cases involve the self advection of a thin vortex ring and the interaction of two vortex rings. These examples confirm the operation of the algorithm without introducing the issue of boundary conditions. Solid boundaries are introduced in the problem of standing waves on a bounded vortex filament. Finally, the problem of a contained fluid with a uniform vorticity distribution is examined by modelling inertial modes in a rotating cylinder. |
Extent | 4966155 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-02-06 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085644 |
URI | http://hdl.handle.net/2429/4213 |
Degree |
Master of Science - MSc |
Program |
Astronomy |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1996-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1996-0127.pdf [ 4.74MB ]
- Metadata
- JSON: 831-1.0085644.json
- JSON-LD: 831-1.0085644-ld.json
- RDF/XML (Pretty): 831-1.0085644-rdf.xml
- RDF/JSON: 831-1.0085644-rdf.json
- Turtle: 831-1.0085644-turtle.txt
- N-Triples: 831-1.0085644-rdf-ntriples.txt
- Original Record: 831-1.0085644-source.json
- Full Text
- 831-1.0085644-fulltext.txt
- Citation
- 831-1.0085644.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085644/manifest