Adaptive Panel Representation for 3D Vortex Ring Motion and Instability

A hierarchical panel method for representing vortex sheet surface motion in 3D flow is presented. Unlike previously employed filament methods, each panel is a leaf of the tree, so it can be subdivided locally, which allows an efficient adaptive point insertion. In addition, we developed curvature-based insertion criteria which allow to localize point insertion to the most complicated curved regions of the sheet. The particles representing the sheet are advected by a regularized Biot-Savart integral with Rosenhead-Moore kernel. The particle velocities are evaluated by an adaptive treecode algorithm based on Taylor expansions in Cartesian coordinates due to Lindsay and Krasny (2001). The method allows to consider much later stages of a vortex ring instability, which may shed light on this complicated flow phase directly leading to the turbulent flow.


Introduction
The behavior and inherent beauty of vortex rings have fascinated researchers for a long time.The most familiar example of a vortex ring is the smoke ring produced when cigarette smoke is ejected suddenly through the lips of a smoker.Dolphins produce vortex rings for play (Marten et al. [1]).They puff out bubbles from their blowholes.The pressure from below overcomes the surface tension of a spherical bubble, punching a hole in the center to create a ring shape.In a laboratory, the usual method of generating vortex rings is to eject fluid impulsively through some type of orifice into a fluid at rest.The ring is made visible by introducing some dye or smoke in or around the orifice.Vortex rings are also formed when powerful and often hazardous line vortices produced by the wing tips of a large aircraft begin to reconnect.
Part of the fascination of vortex rings comes from their compact and persistent nature.It was this persistence and apparent stability that in 1867 prompted Thomson [2] to propose the theory of vortex ring atoms to explain spectral lines in terms of different modes of oscillation of vortex rings.Even though this theory was later superseded by Quantum Mechanics, it inspired the early analysis of the vortex rings that is still relevant today.In many types of turbulent flows, coherent vortex structures exist, and therefore representing turbulence as a superposition of interacting vortices has been actively investigated (Roshko [3], Browand and Weidman [4]).In addition, accelerating ions in super-fluid helium create quantized vortex rings (Rayfield and Reif [5]).In fact, the theory of line vortices and vortex rings is part of the modern macroscopic treatment of liquid helium II (Roberts and Donnelly [6]).Minota et al. [7] studied the acoustics of mutually interacting vortex rings and rings interacting with sharp edges and bluff bodies.
While many studies of vortex rings have been prompted by scientific curiosity, others have technological applications in mind.For example, vortex rings have been suggested as a means for extinguishing gas and oil well fires by Akhmetov et al. [8] and cavitating vortex rings, produced by exciting cavitating jets, have been used for underwater cleaning and rock cutting (Chahine and Genoux [9]).Currently, Settles et al. [10,11] pursue a project supported by the Transportation Security Agency on using vortex rings to shuffle a person's clothes to catch possible microscopic explosive elements.
In theoretical and numerical considerations, vortex rings are often represented by vortex sheets as a model of shear layers.The shear layer is replaced by an idealized jump discontinuity across the sheet surface.The motion of the layer is represented by the selfinduced motion of the sheet.In this article, we employ the Lagrangian particle method to track the motion of vortex sheet; for a review see [12][13][14][15][16].
Following Lindsay and Krasny [17], this method results in a large system of nonlinear ordinary differential equations.It describes a collection of N particles with pairwise interactions-an N-body problem.In our problem, it is necessary to evaluate sums of the form where x i are particle positions, w j is a vector-valued weight associated with the jth particle, and δ is a smoothing parameter.Computing these sums directly requires O(N 2 ) operations.In our simulations, N reaches 10 6 , so direct summation is impractical.To circumvent this problem we employ the adaptive tree code developed by Lindsay and Krasny [17] to evaluate the sums to a specified error tolerance with only O(N logN) operations.Moreover, we parallelized the code to obtain faster execution.The article makes several contributions described below.Lindsay and Krasny [17] used a filament representation of the vortex sheet, which was limited by the non local filament insertion.In the new approach, we use a hierarchical tree-based panel method to represent and update the vortex sheet surface adaptively and truly locally by using a tree of panels.Each panel is a leaf of the tree and thus it can be refined independently allowing localization of point insertion.In addtion, we developed a new curvature-based point insertion scheme which inserts new particles in regions of higher curvature.Using this method allowed us to simulate the late stages of a single-ring instability, which has not been possible with the filament methods before.The ring instability considered here eventually leads to turbulent behavior, so our simulations could prove useful in explaining turbulent stages of vortex rings experimentally studied by Maxworthy [18][19][20].The method developed is versatile-one only has to change initial conditions for the vortex sheet and many other relevant flows such as wakes and jets can be considered.These would be interesting projects for the future.

Lagrangian parametrization
Our simulations are based on the Lagrangian representation of vortex sheet motion introduced by Kaneda [21] and Caflisch [22].Following Lindsay and Krasny [17], we represent vortex rings as rolled-up vortex sheets.Initially a sheet has a form of unit disk.It is a parameterized surface x(Γ,θ,t) composed of closed material lines, where Γ is circulation across the lines and θ is 2π-periodic parameter along the lines, as shown in Figure 2.1.The circulation distribution follows from the bound vortex sheet associated with potential flow past a circular unit disk and is equal to where r = x 2 1 + x 2 2 is the radial coordinate of the point on the sheet measured from the center r = 0.The distribution (2.1) has a square root singularity in its strength σ ∼ dΓ/dr at r = 1 which makes edge to roll up into a spiral.The parametrization of vortex sheet is then (2. 2) The lines of Γ constant correspond to vortex lines (vortex filaments).Using this parametrization, the equation of motion of vortex sheet become (Lindsay [23]) where x = x(Γ,θ,t), x = x( Γ, θ,t), and K δ is the Rosenhead-Moore kernel [24,25].The right-hand side is the regularized Biot-Savart integral.Equation (2.3) with initial condition (2.2) forms an initial value problem and states that the vortex sheet is a material surface moving in its self-induced velocity field.The Lagrangian method used here involves tracking the vorticity back to time t = 0 through the flow map (Kaneda [21]), so we do not have to update vorticity at each time step.The term ∂ x/∂ θ accounts for vortex stretching.

Discretization and main difficulties
We compare an older filament representation with a new way to discretize vortex sheet which is a hierarchical tree panel representation.This is the main result of our work and we look at it in great detail in the subsequent sections.Here, however, we represent discretization generically as which leads to the main system of ODEs where (D θ x) is finite difference discretization of the θ-derivative and Γ j and θ j are the integration weights.The initial conditions are given on the unit disk as follows: Note that we do not have any boundaries in this setting, so it is pure initial value problem for a large system of nonlinear ordinary differential equations.
The main difficulties and our contribution are as follows.
(1) The main problem is how to choose Γ i ,θ i ,w i .How should we represent the vortex sheet surface so that insertion is local?Once a representation is chosen, we have to decide how the quadrature scheme is implemented.We first briefly discuss the previous filament-based method developed by Lindsay and Krasny [17] and then concentrate on our hierarchical panel-based method.(2) Evaluating the right-hand side sum in (3.2) is an N-body problem.A direct summation approach requires O(N 2 ) operations which quickly becomes prohibitively expensive.We employ a velocity tree code developed by Lindsay and Krasny [17] to reduce operation count to O(N logN).

Filament-based representation
Lindsay and Krasny [17] discretized the vortex sheet by a collection of Lagrangian particles (vortex blobs) x i (t), i = 1,2,...,N, corresponding to a discretization in the parameter  space Γ,θ, as represented schematically in Figure 4.1.First, the sheet is discretized in circulation Γ to obtain material lines, and then each line is discretized in parameter θ to obtain particles x i (t).Integration on the right-hand side of (2.3) is done using Fubini theorem with trapezoidal rule in Γ and θ in that order.
The main problem is the new point insertion.As the sheet rolls up, new particles must be inserted to maintain resolution.If two adjacent particles on a material line (filament) get too far apart, then a new point is inserted by cubic interpolation.This point insertion is local.However, if two adjacent filaments get too far apart even locally, then the whole new filament has to be inserted.This insertion procedure is global.Both insertion procedures are illustrated in there is big variation of distance between many adjacent filaments and hundreds of thousands of points are wasted to insert only few needed local points and that prompted us to develop a new hierarchical panel-based approach which allows a truly local insertion.

Hierarchical panel-based approach
In this section, we will describe our new approach to vortex sheet surface representation and integration.The surface is now discretized as a set of nested rectangles (hierarchical tree) in parameter space (Γ,θ) and corresponding panels in physical space (x, y,z), see Figure 5.1.We start with one big rectangle in parameter space 0 ≤ Γ ≤ 1, 0 ≤ θ ≤ 2π, which is shown on the left of Figure 5.1.Then rectangles are subdivided recursively.The most important advantage of this approach is that now each panel can be subdivided separately from all the others based on certain local tests.We do not have filaments anymore and all insertions are now local.
Let us denote vertices and face middle points by numbers 1-8, and center point by 9-as described in Figure 5.2.In the following discussion, we will use them as indices for the corresponding Γ, θ parameter values as well as x, y, z coordinates.If in the process of tree creation or point insertion a rectangle has two points on the opposite sides (like 5 and 7 or 6 and 8), then it is split in two and forms two children.If more than two of the points 5-8 are present, the rectangle is split into four children.During the tree creation, the panels are split based just on the distance test on each side.The positions of points on the unit disk are known immediately from their Γ and θ values using the distributions (2.2).At the later time, the panel subdivisions (new point insertions) are done by the same rules, but the positions of points are found by cubic interpolation using the neighboring points which are obtained by neighbors search.We set up a simple search from the root of the tree down to the leafs to find a left, right, up, or down neighbor.If the neighbor is of the same or smaller size, it automatically has a corresponding point, but if it is bigger, then a point must be found by linear approximation; Figure 5.3.There are more efficient methods to find the neighbor (Samet [26,27]), but they require particular tree structure (like one-level difference between the neighbors), while we want to allow any number of neighbors.We found that such recursive searches are fast enough for our purposes O(logN) especially compared to expensive velocity evaluation O(N logN).
We employ two tests to insert point for t > 0. The first one is a simple distance test between the points on the faces 1-2, 2-3, 3-4, and 4-1.However, as the rings roll up, the vortex surface becomes very curled and a curvature test is necessary.We used 2D vortex sheet roll up of the airplane vortices investigated by Krasny [28] to devise a curvature test.After investigating a number of other possibilities, the best criterion was a distance between a linear and cubic interpolations between two points as shown schematically in Figure 5.4.

Integration weights
In the previous filament approach, the Biot-Savart integral was approximated using Fubini theorem and discretized on the vortex filaments (material lines).Here, we have to 8 Mathematical Problems in Engineering change to 2D integration on the rectangles R i j in the parameter space (Γ,θ).The set of all rectangles R i j is the set of all the leaves of the tree in the parameter space, To choose integration weights, consider arbitrary smooth function of two variables f (x, y).Consider double integral of this function on a rectangle [0, First, we consider a basic rectangle given by only four points, see Figure 6.1.We use linear approximation in both x and y (linear-linear) to represent the function f (x, y) on the rectangle, where the constants c 1 ,...,c 4 have to be found by matching to known function values at points 1,...,4.Solving the resulting system of four linear equations in four unknowns and integrating the ensuing approximation (6.2), we obtain Thus the weights at points 1,...,4 are all equal to x 1 y 1 /4 which is shown schematically in Figure 6.1.
For a five-point rectangle with points 1-5 shown in Figure 6.2, we use the following quadratic polynomial approximation with five coefficients: (6.4) Fitting it to points 1-5 and integrating as before, we obtain the integration coefficients shown in Figure 6.2.We will call this method 456-point polynomial.
Analogously, for a six-point rectangle with points 1-5, 8 shown in Figure 6.3, we use the following quadratic polynomial approximation with six coefficients: Fitting the points and integrating, we obtain the integration weights shown in Figure 6.3.
The other cases of five-and six-point rectangles are obtained by rotating the cases we considered here.To calculate all the weights, we create another global parallel array of weights w j , j = 1,...,N, with the same indexing as the parameter Γ and θ arrays as well as coordinates (x, y,z).Then we go recursively through the leaves of the tree.For each leaf rectangle we calculate appropriate weights and add them to the global array of weights.Thus each point receives weight contribution from two, three, or four rectangles.
The derivatives ∂y/∂θ are calculated by finding (with search) left and right θ-neighbors for each point and using 3-point finite difference derivatives on unequally spaced points as shown in Figure 6.4.
Similar to the weights, we have to define global parallel array of derivatives.Obviously, for each point the derivative has to be calculated only once, so we also define logical global parallel array which tells for each point if its derivative has been assigned already.Note that to obtain the closest neighbors of the points we have to be at the lowest leaves level of 10 Mathematical Problems in Engineering  the tree.Therefore, we first go recursively through the tree until we get to a leaf.For each point of a given leaf (rectangle), we first check if its derivative has been already assigned.
If not, find the neighbors as described before and assign the derivative.

Triangle-based integration
The integration on a rectangle assumes that the corresponding panel in physical space is also approximately quadrilateral.However, at the late time of vortex rings collisions, this is not true-the panels become very distorted.We measure the distortion by the ratio of diagonals dist 13 / dist 24 .If this ratio is greater than some predefined number (in our case 2), then we consider panel to consist of two subtriangles which are cut so that they are most close to equilateral triangles.See Figure 7.1 for the four-point case with two subtriangles.
Let us consider just one triangle in parameter space as shown in Figure 7.2.By fitting a linear approximation f (x, y) = ax + by + c (7.1) to the three known points at 1, 2, and 3, we obtain integration weights shown in Figure 7.2.This allows us to employ a mixed mesh, there the triangular panels are used in highly skewed regions.

Advantages:
(i) adaptive, efficiently resolves local features; (ii) T-junctions in both directions; (iii) high density only in curved regions; (iv) large reduction of N. Disadvantages: (i) complex data structure-all recursive.The local adaptive hierarchical panel-based quadrature and point insertion scheme for 3D vortex sheet motion presented in the previous sections is an original method.However, we drew upon previous work on adaptive surface representation in computational geometry.Brady et al. [29] used an advancing front technique and curvature adaptive triangularization to represent rolling up jets.C 1 continuity was maintained across triangles using cubic Bezier triangular interpolants.However, the front technique is quite expensive requiring O(N 3/2 ) operations.Losasso et al. [30] produced a simulation of water and smoke flow with an unrestricted octree data structure.They proposed new techniques for creating a symmetric positive definite linear system for Poisson's equation on the octree grid.Line and Brown [31] used an octree data structure to produce a high-resolution wake model behind a helicopter with vorticity transport equation.The flexible nature of an octree allowed focusing on the most complicated part of the flow near the rotor.Klaas and Shephard [32] applied the same idea of an octree grid to 3D discretization for partition of unity representation of complicated mechanical engineering junctions.Cristini et al. [33] developed a curvature-based adaptive mesh algorithm for evolving surfaces and applied it to simulations of drop break-up and coalescence.They maintained the local length scale through minimization of a mesh energy function.
As far as surface representation is concerned, probably the most influential for us was the work of Sederberg et al. [34,35] who used T-splines to represent the surface more efficiently than by a tensor grid.The idea of allowing T-junctions in Γ and θ came out of that work.Pauly et al. [36] used an efficient simplification of the point-sampled surfaces to represent complicated artistic shapes.They were able to concentrate more samples in the regions of high curvature as opposed to inefficient tensor grids.One direction for future work is to apply computational geometry techniques to vortex sheet surface representation.

Velocity tree code and parallel implementation
As we mentioned in Section 6, the right-hand side of (2.2) would take O(N 2 ) operations to evaluate directly.We employed a variation of a variable order velocity tree code developed by Lindsay and Krasny [17], which uses O(N logN) operations.The slight change we have made was the use of the same accuracy parameter as we descent down the tree, while they reduced in a way specified in their formula (35).We tried it both ways and reducing it did not seem to be advantageous.
To improve performance we parallelized the code using MPI (Gropp et al. [37], Marzouk and Ghoniem [38]).The main idea was to tell each processor to apply the tree code only to a predefined fraction of the points.We divide the number of points N by the number of available processors minus one (one of them is the master processor)-numproc-1, and obtain the number of particles per processor-numparproc. Then each slave processor i receives all the points and weights from the master, but applies the tree-code to find the velocity of only its fraction of points.The resultant forces are sent back from the slaves and are collected in the master processor; then a new time step is taken.With the exception of the O(N logN) tree code, the rest of the program is fast-O(N) operations, which are done on the master processor.The parallelized code was successfully run on University of Michigan AMD Opteron clusters.We would consistently obtained a 3-4 times speedup on 5-6 processors.

Motion of a circular vortex filament
Our goal is to investigate the evolution and stability of vortex rings.In the past, a ring was often modeled by a collection of circular vortex lines called vortex filaments.Here, we have a vortex sheet represented by a nested tree structure of rectangles, which is obviously quite different from the filaments.However, the behavior of an individual filament simulates the behavior of a vortex ring torus.In addition, the motion and stability of such filaments can be investigated analytically as was done in Lindsay thesis [23].Here we use his result on the instability of a particular wavenumber of given value of smoothing parameter δ and radius of the filament R.
We consider a generic perturbation only in the z-direction: In Figure 10.1, we present the evolution of such a perturbation on a circular vortex filament of the radius R = 0.75.We consider the solution at times t = 1,2 where interesting nonlinear dynamics starts to develop.On each picture, we show the initial condition perturbation at t = 0 and its development at later times t = 1,2.Using Lindsay [23] stability diagram (2.8), we obtain that for δ = 0.1 and given radius, the wavenumber k = 8 is unstable.For the stable k = 5 case, the initial condition is only slightly changed and rotated.On the other hand, in the unstable k = 8 case, we see the development of a complicated hairpin structure analogous to the results obtained in Knio and Ghoniem [39] using the vortex method.
There have been a number of numerical investigations of a vortex rings stability.To name just a few Knio and Ghoniem (1991) used vortex vector elements, transported in Lagrangian coordinates.The elements change vorticity by local stretch, while their direction is governed by the tilting of material lines.Shariff, Verzicco, and Orlandi (1994) solved Navier-Stokes equations.They observed multiple bands of the wavenumbers amplified with higher-order radial modes.Lifschitz et al. [58] used solution of the differential equation for an exactly propagating ring.Brady et al. [29] employed Lagrangian triangular mesh with cubic Bezier triangular interpolants and adaptive refinement curvature.At each time step, an advancing front technique with automatic mesh refinement was used to remesh the vortex sheet.Their mesh generation is quite expensive-O(N 3/2 ) operation, but it is based on surface curvature and stretching, which produces good representation of the surface.
As we discussed before, we apply a hierarchical panel Lagrangian method.The parameters for the tree code we used are somewhat different from Lindsay and Krasny [17].We determined that we need a higher-order accuracy to go further in time, so the tolerance for the tree code was = 10 −5 , compared to only = 10 −3 .The maximum number of particles per node was N 0 = 1000, compared to N 0 = 500.Finally, the maximum order of the Taylor expansion-p max = 8 -same as in Lindsay and Krasny [17].
We present numerical simulations of a perturbed rolling-up vortex sheet.An azimuthal perturbation was introduced by perturbing the z-coordinate of a flat circular disk in the x − y plane.In polar coordinates, the perturbation can be written as p(r,θ) = ρr 2 cos(kθ)e z , (11.1) where k is the perturbation wavenumber and ρ is the amplitude of the perturbation.The r 2 factor is essential to smooth the perturbation near the origin.
After the ring rolls up, the radius of the ring at the position of the core can be approximately taken as 0.75, δ = 0.1.Similarly to a one-filament analysis, the unstable mode is k = 8.In our simulation, the vortex ring consists of rolls which are larger than δ; this should spread the vorticity out away from the core.This is analogous to increasing δ, which lowers the wavenumber of the unstable mode.Thus, we can expect a band of unstable modes around k = 8.Lindsay and Krasny [17] considered a band k = 4 to 11.Here we were able to extend some of their calculations to longer time.
After testing with finer parameters, we determined that the distance insertion parameters can be taken as dΓ = dθ = δ = 0.1 and the cubic-to-linear insertion parameters as In this case, and other high wavenumber simulations k = 7,9,10, the outer turns of the sheet are smooth, but the core is highly distorted.Strong instability is observed and in the late time it leads to hairpins.Using their filament-based approach Lindsay and Krasny [17] were able to get only to time t = 6, and even there, only the outer boundary was well resolved, not the The bulging behavior is consistent with the simulations of Knio and Ghoniem [39] and the experiments of Didden [44].Bulges were also observed in the experiments and numerical simulations of azimuthal perturbations to a jet by Meiburg et al. [59].In addition, our simulations match very well the experiments of Dazin et al. [46][47][48].Also to compare with their results, we take the average centerline of the ring and evaluate FFTs of velocity components on a set of equally spaced points on this centerline.Moreover, we evaluate FFTs of the coordinates of a filament which was initially at the rim of the circular disk.The results for the k = 5 (stable) and k = 8 (unstable) are presented in Figure 11.4.We can see harmonics and the superharmonics of k = 5; there is no subharmonics growth.On the other hand, in the case k = 8, we can see a growth of the subharmonics similar to the one observed by Dazin et al. [46][47][48].
Finally, note that the vortex sheet surfaces plotted here are surfaces formed by the material curves which coincided with the vortex lines of the sheet at t = 0.However, since we are using a δ-smoothed Biot-Savart kernel, they are not the actual vortex lines for t > 0.

Conclusions
A new local, adaptive, higher-order, tree-based quadrature and point insertion method for 3D vortex sheet motion has been developed.The main ingredient of the method is a hierarchical tree construction representing the vortex sheet surface.It makes the code adaptive and permits insertion of the new computational points (vortex blobs) locally.In addition, we developed a new curvature-based point insertion criteria based on distance between local linear and cubic approximations.
The method was applied to the motion of vortex rings which are modeled as the rolledup vortex sheets.The local hierarchical tree-based representation allowed us to resolve the long time details of Widnall's instability of a single vortex ring including the development of a complicated nonlinear region with hairpins and many layers of roll-up.
The single biggest challenge we are still facing is how to resolve skewed regions without excessive growth in the number of points.Our hierarchical tree-based panel method updated a vortex sheet surface effectively and locally as long as the panels remain approximately quadrilateral in physical space.However, initially square panels become severely skewed.We address this problem partially by the introduction of mixed quad/triangle panels, but this does not really solve the problem.It appears that local mesh redistribution must be done.
In terms of enhancing the method, there are several steps which could be taken.Particle removal in the regions of high particle concentration should be investigated.Chorin [60,61] employed the removal of vortex hairpins on the filaments.We are not restricted to filaments, so the removal could be even more efficient.Also, our integration is only 2nd order accurate, we could increase the accuracy of integration by considering panels with more computational points per panel.Another possible area for improvement is a better cell-dividing technique for the velocity tree code.The current technique for cell subdivision breaks them by bisecting the cell's bounding box.This approach just uses the dimensions of the cell without any regard for the cell's particle groupings.Investigating these groupings and breaking cells into more natural clusters could be very beneficial for the tree-code performance.
Another interesting topic would be to understand the dynamics of vortex rings better.This would require a more extensive investigation of the parameter space.In particular, it would be interesting to compile a number of similar runs with the smoothing parameter δ being continuously reduced to zero.It is quite computationally expensive, but with the O(N logN) tree code and a parallelization, it is more feasible than before.
In terms of studying different systems, we would like to investigate oblique and headon collisions of vortex rings.The algorithm however is far more versatile and can be used to model other system represented by vortex sheets such as wakes behind airplanes and jets.In addition, it can be generalized to other particle systems where the interaction kernel is different from K δ .The computation of the expansion coefficients would have to be modified, but recurrence relations similar to the ones derived here could be obtained so that the Taylor coefficients could be computed rapidly.

Figure 6 . 1 .
Figure 6.1.Basic rectangle with only four points.Point numbers are inside, weight fractions are outside.

Figure 7 . 2 .
Figure 7.2.One triangle diagram in parameter space.Inside: point numbers.Outside: weights (factor area of rectangle with these sides).