Regularized Vortex Sheet Evolution in Three Dimensions

JOURNAL OF COMPUTATIONAL PHYSICS ARTICLE NO.

146, 520–545 (1998)

CP985998

Regularized Vortex Sheet Evolution in Three Dimensions M. Brady,1 A. Leonard, and D. I. Pullin Graduate Aeronautical Laboratories, California Institute of Technology, Pasadena, California 91125 Received December 9, 1997; revised March 31, 1998

A computational method is presented to track the evolution of regularized threedimensional vortex sheets through an otherwise irrotational, inviscid, constant-density fluid. The sheet surface is represented by a triangulated mesh with interpolating functions locally defined inside each triangle. C 1 continuity is maintained between triangles via combinations of cubic B´ezier triangular interpolants. The self-induced sheet motion generally results in a highly deformed surface which is adaptively refined as needed to capture regions of increasing curvature and to avoid severe Lagrangian deformation. Automatic mesh refinement is implemented with an advancing front technique. Sheet motion is regularized by adding a length scale cutoff to the Biot– Savart kernel. Toroidal and periodic-cylinder vortex sheets are simulated, modeling vortex rings and vortex/jet combinations, respectively. Comparisons with previous 2D results are favorable and 3D results are presented. °c 1998 Academic Press Key Words: vortex dynamics; advancing front; Kelvin–Helmholtz instability; Lagrangian.

1. INTRODUCTION

Many fluid flows of interest are characterized by regions of concentrated vorticity arising from shear layers. Away from boundaries, these free shear layers diffuse and roll up, mixing the streams of fluid on either side. Examples include wakes of bluff bodies, mixing of two adjacent streams, and jet flows. One means of computing the development of shear layers is to discretize the vorticity field with Lagrangian markers. These markers move with the fluid and, furthermore, with an incompressible fluid, the velocity field can be calculated from the vorticity and boundary conditions alone. The presence of viscosity mandates an update of the vortex strength at each marker as time progresses, but an inviscid approximation is valid for relatively short times in shear layers at high Reynolds number. This paper is concerned with inviscid, constant-density fluids within which shear layers, modeled as surfaces of vorticity, 1

Corresponding author. E-mail: [email protected]. 520

0021-9991/98 $25.00 c 1998 by Academic Press Copyright ° All rights of reproduction in any form reserved.

REGULARIZED VORTEX SHEET EVOLUTION

521

or vortex sheets, move as material elements. The fluid is considered to be irrotational, save for the sheet surface itself. Henceforth the term “vortex sheet” will refer to an isolated surface of vorticity, while the modifier “regularized” will imply an approximation to the vortex sheet, as described below. Previously, the motion of vortex sheets in two and three dimensions has been computed with a variety of methods. Two-dimensional vortex sheets have been modeled with a thin layer of constant vorticity [1], clouds of point vortices [2], and vortex blobs [3]. Three-dimensional sheets, as well as other vortex structures, have been simulated with particles [4], thin-filaments [5, 6], a triangulated surface [7], and a level set method [8]. The linear stability of a flat, two-dimensional vortex sheet is ill-posed in the sense of Hadamard [9]—perturbations grow linearly with their wavenumber. Additionally, the nonlinear evolution of vortex sheets has been shown to develop a curvature singularity in finite time [10–13], bringing into question the computability of long time vortex sheet evolution. As an alternative to the incorporation of viscosity into the shear layer model, regularization of the initial value problem along the lines of Rosenhead’s length scale cutoff of the Biot–Savart kernel [14] has been proposed [3, 15]. Sufficient regularization not only attenuates the unbounded growth at large wave numbers, but it can also remove the finite time singularity, allowing computation of long time evolution [16]. The price is that the Euler equations with an embedded vorticity surface are no longer being simulated, but rather, some effective model equation defined by the regularization. Three-dimensional simulations to date have generally involved sufficient smoothing of the Biot–Savart kernel to damp the smallest wavelengths, avoiding singularity formation issues. Each method is forced to devote more effort than its 2D counterpart, however, with regard to “mesh” refinement, in order to maintain accurate vorticity field representation despite three-dimensional Lagrangian deformation. More freedom is available in choosing which physical flows to model, including the standard shear layer, jets, and vortex tubes, vortex rings, and shocked density interfaces. The present work is an extension of [7], adding higher order sheet representation and global, automatic mesh refinement. The following section provides the necessary background material and describes the governing equations for vortex sheets, as well as invariants of the system. Section 3 details the advancing front implementation and discretization of the governing equations. Section 4 is an analysis of the method and Section 5 presents select case studies.

2. THEORETICAL DEVELOPMENT

A brief derivation of the equations governing the evolution of vortex sheets is given in the following three subsections. The first two provide the necessary background information from differential surface geometry and fluid kinematics, respectively. 2.1. Parametric Surface We employ a parametric surface description 

 X X(ξ, η) =  Y  (ξ, η), Z

(2.1)

522

BRADY, LEONARD, AND PULLIN

FIG. 1. Mapping from parameter to physical space.

where, as shown in Fig. 1, (X, Y, Z ) represents a point in (x, y, z) space on the surface and (ξ, η) ∈ 6 is a point in a 2D isoparametric space. The first fundamental coefficients [17] are E = Xξ · Xξ F = Xξ · Xη

(2.2)

G = X η · Xη , so that an infinitesimal segment in parameter space with components (dξ, dη) is mapped to a corresponding segment of length d X in physical space, where d X 2 = E dξ 2 + 2F dξ dη + G dη2 .

(2.3)

Subscripts ξ and η denote partial differentiation with respect to the corresponding parameter. Two other useful quantities are the Jacobian of the mapping J = |Xξ × Xη | =

p

EG − F2

(2.4)

and the surface normal n = Xξ × Xη /J.

(2.5)

Points for which J = 0 are singular, such as the north and south poles of a sphere parameterized with latitude and longitude. The Gaussian curvature K is K = κa κb =

L N − M2 , J2

(2.6)

where κa and κb are the principal curvatures of the surface, and L , M, and N are the second fundamental coefficients defined by L = n · Xξ ξ M = n · Xξ η N = n · Xηη .

(2.7)

REGULARIZED VORTEX SHEET EVOLUTION

523

The mean curvature µ is µ=

1 E N + G L − 2F M (κa + κb ) = , 2 2(E G − F 2 )

(2.8)

and the square of the total curvature κT is given by κT2 = κa2 + κb2 = 2(2µ2 − K ).

(2.9)

Note that the total curvature is sometimes defined as the integral of the Gaussian curvature over the surface [17]. 2.2. The Velocity Field Consider an unbounded uniform density and inviscid fluid. Vorticity is defined as the curl of the velocity field v = ∇ × u, and the circulation around a closed loop ∂D (vorticity flux through D) is Z I 0 = u · d` = v · n dD. ∂D

(2.10)

(2.11)

D

Combining the curl of Eq. (2.10) with the incompressibility condition (∇ · u = 0) gives a Poisson equation for velocity, ∇ 2 u = −∇ × v, the solution of which is 1 u(x) = 4π

Z V

v(x0 ) × (x − x0 ) 3 0 d x, |x − x0 |3

(2.12)

(2.13)

where V is the vorticity containing region [18]. A homogeneous solution of Eq. (2.12) also exists, in the form of a velocity potential; this solution is omitted as it is not relevant to the remaining discussion. 2.3. Vortex Sheet Motion A vortex sheet is a surface of tangential discontinuity in the velocity field. This mathematical idealization is motivated by a limiting process in which the thickness of a vorticitycontaining layer is decreased and the magnitude of the vorticity in the layer is increased, while their product (thickness and magnitude) remains constant. For parameterized vortex sheets, the vorticity vector takes a particular form, transforming Eq. (2.13) into an equation for the evolution of the sheet itself [7, 19, 20]. Consider a vortex sheet which divides an irrotational, inviscid fluid into two regions, as shown in Fig. 2. The velocity in each region is given completely by a velocity potential, i.e. u+ = ∇φ + and u− = ∇φ − . Vorticity v must lie tangent to the surface and may be defined as v = gδ(n),

(2.14)

524

BRADY, LEONARD, AND PULLIN

FIG. 2. Vortex sheet dividing an ideal fluid.

where n is the distance normal to the sheet and δ is the Dirac delta function. g is the strength of the sheet given by g = n × [u],

(2.15)

where n is the unit surface normal (pointing to the + side) and the brackets denote the jump across the sheet, i.e. [u] = u+ − u− , evaluated on the sheet surface. Consider the circulation around a curve piercing the vortex sheet at points A and B. As the velocities on either side of the sheet have potentials and the vorticity is nonzero only on the sheet, the circulation depends only on the location of points A and B: Z

B

0(A, B) =

Z

B

[u] · d` =

A

(∇φ + − ∇φ − ) · d`

A

= [φ](B) − [φ](A).

(2.16)

If the vortex sheet surface is represented by X(ξ, η), with (ξ, η) ∈ 6, the circulation can be regarded as a point-wise quantity by taking A to be a fixed reference location in parameter space, i.e. 00 = [φ](A) and 0(ξ, η) = [φ](ξ, η) − 00 .

(2.17)

With this definition, 0 is necessarily a single-valued function of (ξ, η). Thus, some care must be taken when considering sheets that contain irreducible curves such as a torus and an infinite cylinder [18]. Fortunately, only the circulation gradients are required in the following discussion. Defining 0ξ = [∇φ] · Xξ , 0η = [∇φ] · Xη and using Eq. (2.5) as the normal direction, it follows from Eq. (2.15) that g(ξ, η) = (0ξ Xη − 0η Xξ )/|Xξ × Xη |.

(2.18)

Thus, Eq. (2.13) becomes 1 u(x) = 4π

Z 6

W0 × (x − X0 ) d6 0 , |x − X0 |3

(2.19)

where W = 0ξ Xη − 0η Xξ ,

(2.20)

REGULARIZED VORTEX SHEET EVOLUTION

525

and all upper case vectors are functions of (ξ, η, t), where t is time. An evolution equation for the sheet is defined by restricting x to lie on X(ξ, η). With Eq. (2.19) we have ˙ = u(X). X

(2.21)

Although the sheet vorticity W changes with time, it can be shown that 0 is invariant, since vortex lines (those everywhere tangent to v) move with the sheet [18]. Thus, in this formulation, no explicit updating of the vorticity is required. Note that although the integrand in Eq. (2.19) is singular when the velocity of the sheet is desired (x = X0 for some (ξ 0 , η0 ) combination), a principal value can be defined which gives the principal velocity of the vortex sheet, the average of the velocities above and below the sheet. Experience in 2D [12] shows that time integration of Eq. (2.21) is inherently unstable due to amplification of the roundoff error, and a means of regularization is necessary. One method is to introduce a length scale σ into the denominator of the Biot–Savart kernel, eliminating the singularity. The modified equation is uσ (x) =

1 4π

Z 6

W0 × (x − X0 ) d6 0 , [|x − X0 |2 + σ 2 ]3/2

(2.22)

which as before leads to the evolution equation ˙ = uσ (X). X

(2.23)

It has been shown that simulations with this modified Biot–Savart kernel capture the qualitative features of smooth flows [16]. Additionally, features of viscous shear layer simulations in the limit of increasing Reynolds number and decreasing initial thickness are in good agreement with regularized vortex sheet calculations as the Biot–Savart cutoff length tends to zero [21]. The effect of σ is discussed further in Sections 4.2 and 5.1. All variables mentioned hereafter are assumed to be nondimensional, scaled with appropriate combinations of a characteristic length L 0 and circulation 00 , unless noted otherwise. An initial value problem is now clearly defined: calculate the time evolution of X(ξ, η, t), given the initial sheet position X(ξ, η, 0) and the relative circulation 0(ξ, η). 2.4. Surface Types In order to be computationally tractable, only surfaces that can be described with a compact parameter space domain 6 are considered. Surfaces with infinite extent are thus sliced at regular intervals and constrained to periodic motion. Categorization of surfaces is based on the prescribed constraints in parameter space. The first surface type considered is that which results from a mapping function X defined over a square parameter space domain, constrained to be periodic in both ξ and η. X is then a torus-like surface, although possibly so contorted as to contain a knot. This type of surface is used, for example, to follow the motion of hollow cored vortex rings. The second surface type is similar to the first, but allows one component of X to have an additional linear variation in ξ . This results in a spatially periodic-train of surfaces, connected smoothly from end to end in physical space. This type of surface is used to track hollow-core vortex tubes, jets, and combinations of the two. As an example, consider the parameter space (ξ, η) ∈ (−∞ ·· ∞, −π ·· π ), where the

526

BRADY, LEONARD, AND PULLIN

periodic-train surface is sliced at intervals of 2π in ξ , unit length in x, and points along the x axis. Then Eq. (2.22) becomes 1 uσ (x) = 4π

Z π Z∞

−π −∞

=

1 4π

W0 × (x − X0 ) dξ 0 dη0 [|x − X0 |2 + σ 2 ]3/2

Zπ Zπ X ∞

−π −π

W0 × [x − (X0 + k ˆi)] dξ 0 dη0 , 0 + k ˆi)|2 + σ 2 ]3/2 [|x − (X k=−∞

(2.24)

where ˆi is the unit vector along the x axis. Integration of Eq. (2.24) now requires evaluation of the summations s1 (b, c) =

∞ X k=−∞

s2 (b, c) =

1 [(k −

b)2

[(k −

b)2

∞ X k=−∞

(2.25)

+ c2 ]3/2

k + c2 ]3/2

,

(2.26)

where b and c are combinations of variables from Eq. (2.24), but are treated as constants in the summations. Additional surface types not yet implemented include a spherical-like surface with boundary conditions suitable for a polar coordinate origin in η, but simple periodic constraints in ξ . This surface type would be useful for examining the motion of shocked density interfaces [22]. Extending the periodic-train type, an additional linear component in η may be added to X, creating a doubly periodic sheet in physical space. This surface type is the closest model of a standard shear layer. 2.5. Invariants Uniform density, inviscid fluids with an embedded vortex sheet maintain an infinite number of invariant quantities, although a select group of integrated quantities [23] provide a good set of diagnostics for the numerical method. When considering the motion of regularized vortex sheets evolving with Eq. (2.22), rather than Eq. (2.19), care must be taken to define the correct diagnostics [4]. The following quantities are preserved for all time with an exact spatial and temporal integrator for Eq. (2.23) (see Appendix). Unit density is assumed for the fluid. We introduce the subscript δ to refer to the velocity field (uδ ) induced by a true vortex sheet with vorticity vδ = gδ(n), as in Eqs. (2.14) and (2.19): Volume: 1 V = 3

Z VS

1 ∇ · x dV = 3

Z X · (Xξ × Xη ) d6

(2.27)

6

Kinetic energy: 1 T = 2

Z V

1 uσ · u δ d V = 8π

Z Z 6

60

W · W0 d6 d6 0 |X − X0 |

(2.28)

REGULARIZED VORTEX SHEET EVOLUTION

527

Helicity: Z H=

Z uσ · vδ d V =

U · W d6

(2.29)

6

V

Total vorticity: Z

Z vδ d V = V

W d6 = 0

(2.30)

6

Hydrodynamic impulse: 1 I= 2

Z

1 x × vδ d V = 2

V

Z X × W d6

(2.31)

X × (X × W) d6.

(2.32)

6

Angular impulse: A=

1 3

Z x × (x × vδ ) d V = V

1 3

Z 6

Here VS represents the volume inside the vortex sheet surface and V is all space with nonzero vorticity. For the periodic-train surface type, the above definitions need slight modifications to account for physical space periodicity. 3. COMPUTATIONAL DEVELOPMENT

The integrand of Eq. (2.22) is approximated over an irregular, triangulated mesh covering the parameter space domain 6, as in the left of Fig. 3a. Mesh refinement (resulting in nonuniform distribution of grid points) is essential in order to capture local areas of high deformation in the physical sheet which develop during the sheet evolution (see Fig. 4). Parameter space triangulation (both vertex placement and connection) is accomplished by an advancing front method. The vortex sheet surface is re-triangulated as needed, based on surface curvature and stretching, with triangle size and orientation conforming to local properties of the surface. Mesh postprocessing techniques are implemented to increase desired characteristics of the mesh. The nine functions X(ξ, η), Xξ (ξ, η), and Xη (ξ, η) are interpolated in parameter space by a combination of cubic B´ezier triangular patches [24], providing C 1 continuity between triangles in X and C 0 continuity in Xξ , Xη . The spatial integration in Eq. (2.22) is approximated with Gauss quadrature points [25] and time integration proceeds with standard Runge–Kutta methods. Note that in the figures, vortex sheet surfaces are, with one exception, depicted simply with their triangular mesh skeleton. Only the triangle vertices are necessarily on the interpolated sheet surface. An example of the smooth sheet surface is shown in Fig. 11f. 3.1. Advancing Front The advancing front technique is a powerful method for creating a mesh of triangles conforming to user-specified size, stretching, and boundary constraints [26, 27]. In a general sense, closed directed boundaries are given as input such that the right-facing normal

528

BRADY, LEONARD, AND PULLIN

FIG. 3. (a) Raw and (b) postprocessed parameters and the corresponding physical space for a torus.

consistently points inside the domain. Nodes are then placed along these boundaries (fronts) at a specified, and possibly varying, spacing, depending on the surface being triangulated. Each edge of the front (between nodes) is considered in turn as one side of a new triangle, and the ideal location of the triangle’s third vertex is calculated following the given size and stretching criteria (described below). The front is modified as triangles are generated so as to close in on itself. This particular advancing front implementation fills a square region in parameter space with specifically sized and oriented triangles such that their mapping into physical space creates the desired surface. Additionally, the physical space triangles are designated to be nearly equilateral and have their area vary proportionally with the inverse of the total surface curvature squared.

REGULARIZED VORTEX SHEET EVOLUTION

529

FIG. 4. Demonstration via mesh skeletons of the need for mesh refinement. Lagrangian deformation of material elements in (a) reveals a loss of resolution with a fixed mesh (600 triangles). With the same initial condition, timely remeshing in (b) allows for an accurate representation of the sheet (5000 triangles).

3.1.1. Density function. We define a target triangle density function f (ξ, η) as f = κT2 ,

(3.1)

representing a conservative estimation of the principal curvatures. This type of scaling keeps the same number of triangles on a sphere, independent of its radius. Bounds are placed on

530

BRADY, LEONARD, AND PULLIN

the density function in order to (1) inhibit excessively large triangles and (2) fix an upper bound on the possible number of triangles generated. 3.1.2. Triangle size scaling. Given the above density function, it is possible to write an approximate integral relation between the number of triangles and the triangle density as Z Z f dS = β f J d6, (3.2) N∼ β = 6

S

where N is the number of triangles associated with the (dimensionless) scaling factor β, density function f , and mapping Jacobian J . For the initial condition, N , f , and J are specified so that β may be calculated. β is then held constant as the sheet evolves, thereby fixing the local relation between density and triangle area. In fact, under the assumption that triangles are equilateral in physical space with side length 1X , the local triangle area must be √ 1 3 1X 2 = , 4 βf and the local physical space length scale is thus 1X =

2 1 ∼ 1.52 √ =√ . 31/4 β f βf

(3.3)

The resolution of a surface discretized with this procedure is best characterized with β, as β, not N , gives a better basis for comparison of different surface geometries. For example, consider a sphere for which we have f = 1/R 2 ; Eq. (3.2) gives N ∼ = 4πβ/3, independent of the radius. 3.1.3. Stretching directions. In order to generate equilateral triangles in physical space, the first-order stretching and rotation components of the mapping are used to define an appropriately chosen small ellipse in parameter space which maps to a circle in physical space. The mapped length 1X of a parameter space segment of length 1ξ at an angle α above the horizontal is approximated by (2.3) to first order as 1X = 1ξ

p

E cos2 α + F sin 2α + G sin2 α.

(3.4)

The two extreme values of 1X (major and minor axes of the ellipse) are obtained when α satisfies tan 2α =

2F . E−G

(3.5)

When E = G there is no stretching and 1X has no extreme values. 3.1.4. Boundary and ideal node placement. In order to place nodes along a line (at angle α above the ξ axis) in parameter space at a spacing corresponding to the density function f , we simply evaluate (3.3) to get the desired physical space length and then solve (3.4) for 1ξ . This is done repetitively, marching along the line. The end of the line is hit exactly by scaling the position of all points generated based on the closest to the end. Given a directed segment in parameter space, a node exists which, together with the segment, defines a triangle. This node is “ideally” placed when the triangle’s mapping is equilateral.

REGULARIZED VORTEX SHEET EVOLUTION

531

The location of this ideal node is found by rotating and scaling the coordinate system such that the new basis vectors correspond to the local stretching magnitudes and directions via (3.4) and (3.5). An equilateral triangle is generated with respect to this new basis and the node’s position is transformed back to the unmodified basis. 3.1.5. Implementation. Implementation of the advancing front method begins by defining boundary nodes of the parameter space domain (ξ = ± π, η = ± π for example). Since the domain is periodic, only two sides are defined with the technique described previously; the other two sides are redundant copies. A single front is defined as the clockwise traversal of the segments between nodes. The front then advances by examination of each segment in turn. For each segment, the corresponding ideal node location is calculated. A new triangle is generated and the front modified only if the new triangle edges do not intersect existing edges and the new node is past a threshold distance from existing nodes. Otherwise, the closest existing node within the front is used to define the new triangle. Figure 5 shows

FIG. 5. Advancing front in parameter space 6 with 185, 1000, 2000, and 2650 triangles. Note how in (d) the front has pinched itself closed, generating two child fronts.

532

BRADY, LEONARD, AND PULLIN

FIG. 6. Dashed lines represent (a) swapped and (b) unswappable edges.

various stages during an example mesh creation. Note that the front may pinch itself off, generating two child fronts, as seen in Fig. 5d. The initial mesh, as well as subsequent re-meshes, are generated with this procedure. While generating and refining a new mesh, the old existing one remains active in order to evaluate the density and stretching directions and fill the new mesh with X. The mesh generation process takes O(N 3/2 ) time, as a search through O(N 1/2 ) triangles is necessary for each of the N triangles generated. 3.1.6. Mesh postprocessing. The mesh generated by this advancing front algorithm is not ideal, as can be seen in the physical space representation in Fig. 3a. Whenever two fronts collide, closing a region, ideal nodes are typically not generated and the corresponding triangles in physical space are far from equilateral. But a completed mesh can be improved with a variety of fast operations. Laplace smoothing [28] is the name given to a mesh process in which every node is moved to the center of gravity of its neighboring triangles. In this application, we wish to maintain equilateral triangles in physical space, so vertex coordinates in parameter space are moved to the physical space area-weighted average of the neighboring triangles. The method may be applied repetitively, but a few (<5) iterations result in significant improvement. Mesh relaxation [29, 30] attempts to push the number of triangles adjacent to each node (the degree d) towards six. To this end, every edge between two triangles is considered in turn. The edge is swapped if (a) geometrically feasible (see Fig. 6) and (b) the resulting number of triangles per node is improved, based on a regularity P4 (di − 6)2 . Here i represents the four vertices needed to define the edge index, R = i=1 under consideration and the two adjacent triangles. The third and final mesh improvement technique swaps edges based on the desire to line them up with valleys and ridges of the surface. Define the edge-center-distance (ECD) to be the distance in physical space from the mapping of the edge midpoint, to the middle of the physical space edge. Each edge is considered in turn and is swapped if (a) geometrically feasible and (b) the swapped edge’s ECD is some degree less than the ECD of the current edge. Figure 3b shows the improvement made after two passes of Laplace smoothing followed one pass through ECD swapping and five passes of relaxation (thresholds from 11 to 1 by 2), followed finally by two more passes of Laplace smoothing. 3.2. C 1 Interpolation In order to evaluate the integrand of Eq. (2.22), X and its derivatives (Xξ and Xη ) are interpolated within 6 by a combination of cubic B´ezier triangular patches defined over each triangle [24]. To this end, derivatives are estimated at each vertex (in 6) by means of a parabolic least square fit of neighboring function values, located at the vertices of neighboring triangles. Care is taken to ensure enough data for the fit by including vertices from the next ring of adjacent triangles if necessary. Once the fit is constructed, nodal derivative estimates are extracted from the paraboloid. Once function values and derivative estimates are known for the three vertices of a triangle, an analytic interpolant for the

REGULARIZED VORTEX SHEET EVOLUTION

533

function is defined. The function’s derivatives with respect to ξ and η are obtained from the analytical derivative of the interpolant. The functional interpolant is found to have O(1X 3 ) accuracy while its derivatives are only O(1X 2 ) accurate (see Section 4.1). 3.3. Integration In addition to the interpolant, the time-stepping method and number of Gauss points per triangle control the order of accuracy of each computation. Time stepping is done with one of the Runge–Kutta methods, selectable for second- to fourth-order accuracy. Time stepping is regulated by two concerns: (1) control over triangle edge rotation and straining, and (2) preservation of invariants. The first is implemented by defining the strain for each edge as ²=

|1U |1t , L

(3.6)

where L is the original length of the edge (in physical space) and 1U is the difference in velocities of the edge endpoints. The time step is varied so as to maintain ² < ²0 for every edge. The parameter ²0 is chosen before computation so that invariants are maintained to the desired accuracy (as far as time integration contributes). The integral in Eq. (2.22) is approximated for each triangle vertex in 6. The integrand is evaluated at the Gauss quadrature points for each triangle, finally determining the global influence of the surface on each vertex. Thus, with this straightforward implementation, calculation of the surface velocity field is an O(N 2 ) process. The spatial integration order is easily varied by changing the number of Gauss points used in each triangle (second to fifth order in 1X with 1 to 6 points). Smaller values of σ , or a tightly folded vortex sheet surface lead to near-singular behavior of the integrand, forcing more Gauss points or more triangles in order to maintain the invariants. 3.4. Summation Approximation The summations in Eq. (2.25) and Eq. (2.26) are stored at machine precision in a twodimensional look-up table with the parameters (b, c) as the axes. This is feasible because only a compact space from the (b, c) domain is required. Data is extracted from the table with a third-order accurate interpolant [31]. 4. PERFORMANCE

4.1. Convergence Analysis There are several convergence properties to check. As a test case, the toroidal surface type is used with circulation modeling a vortex ring:   0.2 sin η   X =  (1 + 0.2 cos η) cos ξ  (4.1) (1 + 0.2 cos η) sin ξ 0 = −η

(4.2)

with (ξ, η) ∈ (−π ·· π, −π ·· π ). Here, the initial condition X specifies a torus with a tube to ring radii ratio of 0.2. Figure 7 shows the convergence properties of the interpolation

534

BRADY, LEONARD, AND PULLIN

FIG. 7. Maximum error versus the number of triangles for each component of X (hollow points), Xξ (filled), and Xη (filled). Dashed line has slope m = − 1.5. Solid line has slope m = −1.0.

procedure, combining data for all nine functions in X, Xξ , and Xη . In this case, “error” E is defined as the maximum discrepancy between the exact and interpolated values, as determined by a 10,000 point sample over a regular grid spanning the parameter space. Assuming E is proportional to N m and also to 1X −2m , we see that X is third-order accurate in 1X and the derivatives are second-order accurate. The triangular Gaussian quadrature integration was confirmed to behave as expected, integrating low order polynomials exactly. Finally, all of the evolution diagnostics were compared for a range of time and number of triangles (N = 2000, 4000, 8000). Acceptable convergence properties were observed, as shown in Fig. 8. Vertical jumps in plotted quantities correspond to remeshing steps. Temporal order of accuracy was confirmed to coincide with expected Runge–Kutta theory. 4.2. Comparison with Exact Velocities Consider a cylindrical vortex sheet with unit radius and strength, aligned along the x axis. It may be parameterized as 

 ξ   X =  sin(η)  , cos(η)

(4.3)

REGULARIZED VORTEX SHEET EVOLUTION

535

FIG. 8. Convergence of diagnostics: (a)–(f) correspond to Eqs. (2.27)–(2.32). Initial numbers of triangles N = 2000 (dashed), 4000 (solid), and 8000 (bold). Magnitude is taken for vector quantities.

with (ξ, η) ∈ (−∞ ·· ∞, −π ·· π ). Substitution into Eq. (2.22) yields the velocity as a function of r , the radial distance from the x axis. Specifically,

U=

√ 1 √ (1 + σ 2 − r 2 − α), 2r α

(4.4)

536

BRADY, LEONARD, AND PULLIN

FIG. 9. Comparison of numerical method (points) with exact velocities (lines) from the modified Biot–Savart kernel Eq. (2.22) with σ = 0.25, 0.1, 0.01 (open circles). Curves steepen with decreasing σ . Note the dramatic difference in accuracy for σ = 0.01 in (a) 500 and (b) 8000 triangles.

where α = [(1 − r )2 + σ 2 ][(1 + r )2 + σ 2 ].

(4.5)

Figure 9 shows how the numerical method (points) compares to the above analytic formulas, for decreasing values of σ . As σ decreases, the velocity profile from the regularized vortex sheet approaches that of the true vortex sheet. In this case, the velocity of the vortex sheet is

REGULARIZED VORTEX SHEET EVOLUTION

537

U = 1; thus we confirm the velocity of the regularized sheet approaches that of the vortex sheet as σ decreases. 4.3. Comparison with 2D Results Several comparisons are possible with existing 2D results. Consider, for example, the case of a cylindrical vortex sheet, the strength of which is given by the slip velocity of an ideal fluid as a free stream passes around the cylinder. For comparison, the periodic-train surface type is used with initial condition: 

ξ/2π



  X =  0.2 sin(η) 

(4.6)

0.2 cos(η) 0 = −2 sin(η)

(4.7)

with (ξ, η) ∈ (−π ·· π, −π ·· π ). Each period in physical space has unit length and initial radius to period ratio of 0.2. A separate 2D code, which compares perfectly with the results in [32], is used as a test bed. The current 3D code is able to reproduce similar results, as shown in Fig. 10. The projected 3D solution in Fig. 10a ran with σ = 0.05 (scaled with initial diameter), beginning with N = 1000 and ending with N = 9400 triangles. The projected 3D solution in Fig. 10b ran with σ = 0.025 (scaled with initial diameter), beginning with N = 1000 and ending with N = 27,000 triangles.

FIG. 10. Comparison of current 3D method (solid lines representing surface, projected into y − z plane) with 2D cylindrical vortex sheet (dashed lines): (a) σ = 0.1; (b) σ = 0.05.

538

BRADY, LEONARD, AND PULLIN

5. EXAMPLES

5.1. Jet Modeling a cylindrical jet with a first-mode axial perturbation, we define the initial position and strength as   ξ/2π   X =  0.25[sin(η) + 0.1 sin(ξ )]  (5.1) 0.25[cos(η) + 0.1 sin(ξ )] 0 = −ξ

(5.2)

with (ξ, η) ∈ (−π ·· π, −π ·· π). An initial grid density of β = 500 generates 500 triangles for the initial mesh, as shown in Fig. 11a. Note that the time evolution preserves axisymmetry as seen by eye—an artifact of the relatively large σ (=0.2, scaled with the initial average diameter). No symmetries are imposed by the computational method. By Fig. 11d, we see significant rollup of the sheet, essentially defining a vortex ring around the jet axis. At this time, the sheet is defined with 24,000 triangles, and computation becomes prohibitively expensive. Figure 11e shows the amount of detail and variety of scales present. Figure 11f depicts the interpolated smooth sheet. Figure 14 shows the effect of varying σ while keeping the initial surface type (duplicated from above) constant. Smaller values of σ lead to faster rollup as anticipated by Fig. 9. Also, as σ decreases we see the introduction of small-scale instabilities and a deviation from axisymmetry. 5.2. Leap Frogging Vortex Rings Modeling the interaction of two co-axial vortex rings, we define the initial position and strength of two sheets as   ±0.5 + 0.2 sin η   X =  (1 + 0.2 cos η) cos ξ  (5.3) (1 + 0.2 cos η) sin ξ 0 = −η,

(5.4)

again with (ξ, η) ∈ (−π ·· π, −π ·· π ). This time, with β = 200, 1300 triangles are needed to define the initial toroids, as shown in Fig. 12a. Again axisymmetry is maintained due to the relatively large σ (=0.25 scaled with the initial tube diameter). The initially trailing ring moves through the lead ring, quickly rolling and expanding so that the combination requires 16,000 triangles to define in Fig. 12d. 5.3. Vortex Rings Modeling a fully 3D interaction of two interlocked vortex rings, we define the initial position and strength of the two sheets as   0.25 sin η (5.5) X =  (1 + 0.25 cos η) cos ξ  (1 + 0.25 cos η) sin ξ 0 = −η

(5.6)

REGULARIZED VORTEX SHEET EVOLUTION

539

FIG. 11. Periodic-train vortex sheet with strength distribution simulating a hollow-core jet with initial first mode axial perturbation; σ = 0.2 (relative to initial average diameter). One period calculated, two shown with second in cut-away. Mesh skeleton: (a) t = 0; (b) t = 0.2; (c) t = 0.4; (d) t = 0.6; (e) zoom of roll-up region in (d); (f) actual smooth surface representation of (e).

and 

 (1 + 0.25 cos η) cos ξ  X= 0.25 sin η (1 + 0.25 cos η) sin ξ

(5.7)

0 = −η

(5.8)

with (ξ, η) ∈ (−π ·· π, −π ·· π ). This time, with β = 300, 2000 triangles are needed to define the initial toroids, as shown in Fig. 13a. In this case we have σ = 0.2, scaled with the initial tube diameters. Fluid in the center of the horizontal ring is propelled upward,

540

BRADY, LEONARD, AND PULLIN

FIG. 12. Mesh skeleton of hollow-core vortex rings with in-line impulse vectors; σ = 0.25 (relative to initial tube diameters). Ring to tube radius ratio is 5: (a) t = 0; (b) t = 0.2; (c) t = 0.4; (d) t = 0.6.

causing a fold to develop in the vertical ring. At a later time in Fig. 13c, significant length scale differences exist in different parts of the sheet, as seen in Fig. 13d. Each sheet in Fig. 13d has 12,000 triangles.

6. DISCUSSION

A means for computing the motion of regularized vortex sheets in three dimensions is presented, utilizing automatic mesh refinement to maintain an accurate surface representation despite significant deformation. Smooth (C 1 ) surface representation via interpolation over triangular patches increases the order of the overall computation. Two surface types are demonstrated, modeling vortex rings and vortex jet/tubes. The methods described here can be used to model several physical flows such as vortex ring interaction, trailing wing tip vortex pair interaction, shocked density interface motion, jet and shear layer development and stirring. Additionally it is of interest to examine the range of scales which develop from an initially perturbed sheet and the extent to which the Biot–Savart kernel smoothing parameter σ has an effect. The method presented is O(N 2 ) per time step, as dictated by the straightforward implementation of the spatial integration required in Eq. (2.22); each vertex is influenced by each triangle. This quadratic variation in N dominates, with resect to execution time, other aspects of the method, including remeshing and mesh postprocessing. With a 256MB, R8000 based SGI computer, execution time becomes unreasonable around O(104 ) triangles—long before memory usage becomes a constraint. Unfortunately, the number of triangles required to resolve a vortex sheet rolling up increases rapidly with time, particularly with smaller values of the regularization parameter σ . Figure 15 shows the number of triangles needed versus time for the three jet examples in Fig. 14. Hence, our future work will incorporate some of the elements (spatial subdivision, series expansion of group of triangles) of fast N -body methods [33] in order to reduce the integration time from O(N 2 ) to O(N 3/2 ) or O(N log N ).

REGULARIZED VORTEX SHEET EVOLUTION

541

FIG. 13. Mesh skeleton of interlocked hollow-core vortex rings with perpendicular impulse vectors; σ = 0.2. Ring to tube radius ratio is 4. Vertical ring impulse points right and out of page, while horizontal ring impulse points up: (a) t = 0; (b) t = 0.25; (c) t = 0.5; (d) zoom of center region in (c).

APPENDIX: INVARIANTS

Proof that the integral quantities presented in Section 2.5 are time invariant is demonstrated by first introducing new surface parameters (r, s) such that r is perpendicular and s is parallel to vortex lines. By construction, we have 0s = 0 and W dξ dη = 0r Xs dr ds.

(A.1)

The total vorticity is thus Z Z W dξ dη = ξ

η

(Z

Z 0r r

) Xs ds

s

dr = 0,

(A.2)

542

BRADY, LEONARD, AND PULLIN

FIG. 14. Mesh skeleton of periodic-train vortex sheet with strength distribution simulating a hollow-core jet at t = 0.12. Initial first mode axial perturbation. One period calculated, two shown with second in cut-away: (a) σ = 0.5; (b) σ = 0.25; (c) σ = 0.0625.

since vortex lines are closed on the surface (we consider here the toroidal surface type only). The velocity field becomes ˙ s) = 1 X(r, 4π

Z Z r0

s0

0r 0 Xs 0 × K dr 0 ds 0 ,

(A.3)

REGULARIZED VORTEX SHEET EVOLUTION

543

FIG. 15. Number of triangles N versus time t for the three cases in Fig. 14.

where K = (X − X0 )/D 3 and D = [|X − X0 |2 + σ 2 ]1/2 . For the first moment of vorticity, we have I=

1 2

Z Z X × W dξ dη = η

ξ

1 2

Z Z 0r X × Xs dr ds. r

(A.4)

s

Time differentiation and substitution of the velocity field (Eq. (A.3)) yield 1 I˙ = 4π

Z Z Z Z r

r0

s

0r 0r 0 (Xs 0 × K) × Xs dr ds dr 0 ds 0 .

(A.5)

s0

Finally, symmetry arguments and the identity (∂/∂s)D −1 = −Xs · K give 1 I˙ = − 4π

(Z

Z Z Z 0r 0r 0 r

r0

s0

s

) ∂ −1 D ds Xs 0 dr dr 0 ds 0 = 0. ∂s

(A.6)

˙ = 0. Volume is preserved if (a) the bounding Analogous arguments apply for the proof of A surface moves with the fluid and (b) the velocity field is incompressible. The former is implied in Eqs. (2.22) and (A.3); the latter is apparent from an alternate form of Eq. (2.22), 1 ∇× uσ (x) = 4π

Z 6

W0 d6 0 . [|x − X0 |2 + σ 2 ]1/2

(A.7)

544

BRADY, LEONARD, AND PULLIN

With the new parameters, kinetic energy becomes 1 T = 8π

Z Z Z Z 0r 0r 0 r

s r 0 s0

Xs · Xs 0 dr ds dr 0 ds 0 . D

(A.8)

Time differentiation and symmetry arguments lead to 1 T˙ = 4π

Z Z Z Z r

˙ · Xs 0 ) 0r 0r 0 D −3 {[(X − X0 ) · Xs ](X

s r 0 s0

0 0 ˙ − [(X − X0 ) · X](X s · Xs 0 )} dr ds dr ds .

(A.9)

Substitution of the velocity field (Eq. (A.3)) and some manipulation finally gives T˙ =

1 16π 2

Z Z Z Z Z Z 0r 0r 0 0r 00 Xs r

·

s r 0 s 0 r 00 s 00

[(X − X0 ) × Xs 0 ] × [(X − X00 ) × Xs 0 ] dr ds dr 0 ds 0 dr 00 ds 00 [|X − X0 |2 + σ 2 ]3/2 [|X − X00 |2 + σ 2 ]3/2

(A.10)

which vanishes due to the asymmetry in X0 and X00 . ACKNOWLEDGMENTS This work was partially supported by NSF Grant CTS 9634222. We also thank G. S. Winckelmans for timely, fruitful discussions.

REFERENCES 1. C. Pozrikidis and J. J. L. Higdon, Nonlinear Kelvin–Helmholtz instability of a finite vortex layer, J. Fluid Mech. 157, 225 (1985). 2. L. Rosenhead, The formation of vortices from a surface of discontinuity, Proc. R. Soc. London A 134, 170 (1932). 3. A. J. Chorin and P. S. Bernard, Discretization of a vortex sheet, with an example of roll-up, J. Comput. Phys. 13, 423 (1973). 4. G. S. Winckelmans and A. Leonard, Contributions to vortex particle methods for the computation of threedimensional incompressible unsteady flows, J. Comput. Phys. 109 (1993). 5. W. T. Ashurst and E. Meiburg, Three-dimensional shear layers via vortex dynamics, J. fluid Mech. 189, 87 (1988). 6. O. M. Knio and A. F. Ghoniem, Three-dimensional vortex simulation of roll-up and entrainment in a shear layer, J. Comput. Phys. 97, 172 (1991). 7. M. E. Agishtein and A. A. Migdal, Dynamics of vortex surfaces in three dimensions: Theory and simulations, Physica D 40, 91 (1989). 8. E. Harabetian, S. Osher, and C. Shu, An Eulerian approach for vortex motion using a level set regularization procedure, J. Comput. Phys. 127, 15 (1996). 9. P. R. Garabedian, Partial Differential Equations (Wiley, New York, 1964). 10. D. W. Moore, The Spontaneous appearance of a singularity in the shape of an evolving vortex sheet, Proc. R. Soc. London A 365, 105 (1979). 11. D. I. Meiron, G. R. Baker, and S. A. Orszag, Analytical structure of vortex sheet dynamics. Part 1. Kelvin– Helmholtz instability, J. Fluid Mech. 114, 283 (1982).

REGULARIZED VORTEX SHEET EVOLUTION

545

12. R. Krasny, A study of singularity formation in a vortex by the point-vortex approximation, J. Fluid Mech. 167, 65 (1986). 13. M. J. Shelley, A study of singularity formation in vortex-sheet motion by a spectrally accurate vortex method, J. Fluid Mech. 244, 493 (1992). 14. L. Rosenhead, The spread of vorticity in the wake behind a cylinder, Proc. R. Soc. London A 127, 590 (1930). 15. C. Anderson, A vortex method for flows with slight density variations, J. Comput. Phys. 61, 417 (1985). 16. R. Krasny, Desingularization of periodic vortex sheet roll-up, J. Comput. Phys. 65, 292 (1986). 17. T. J. Willmore, An Introduction to Differential Geometry (Oxford Univ. Press, London, 1959). 18. P. G. Saffman, Vortex Dynamics (Cambridge Univ. Press, Cambridge, 1992). 19. R. E. Caflisch, Mathematical analysis of vortex dynamics, in Mathematical Aspects of Vortex Dynamics, edited by R. E. Caflisch (SIAM, Philadelphia, 1988), p. 1. 20. Y. Kaneda, A representation of the motion of a vortex sheet in a three-dimensional flow, Phys. Fluids A 2, 458 (1990). 21. G. Tryggvason, W. J. A. Dahm, and K. Sbeih, Fine structure of vortex sheet roll-up by viscous and inviscid simulation, Trans. ASME I: J. Fluids Eng. 113, 31 (1991). 22. J. F. Hass and B. Sturtevant, Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities, J. Fluid Mech. 181, 41 (1987). 23. C. Pozrikidis, Introduction to Theoretical and Computational Fluid Dynamics (Oxford Univ. Press, New York, 1997). 24. T. N. T. Goodman and H. B. Said, A C 1 triangular interpolant suitable for scattered data interpolation, Comm. App. Num. Meth. 7, 479 (1991). 25. G. R. Cowper, Gaussian quadrature formulas for triangles, Int. J. Num. Meth. Eng. 7, 405 (1973). 26. S. H. Lo, A new mesh generation scheme for arbitrary planar domains, Int. J. Num. Meth. Eng. 21, 1403 (1985). 27. J. Peraire, M. Vahdati, K. Morgan, and O. C. Zienkiewicz, Adaptive remeshing for compressible flow computations, J. Comput. Phys. 72, 467 (1987). 28. J. C. Cavendish, Automatic triangulation of arbitrary planar domains for the finite element method, Int. J. Num. Meth. Eng. 8, 679 (1974). 29. W. H. Frey and D. A. Field, Mesh relaxation: A new technique for improving triangulations, Int. J. Num. Meth. Eng. 31, 1121 (1991). 30. E. Boender, Reliable Delaunay-based mesh generation and mesh improvement, Comm. Num. Meth. Eng. 10, 773 (1994). 31. M. Abramowitz and I. E. Stegun, Handbook of Mathematical Functions (Dover, New York, 1972). 32. J. W. Rottman and P. K. Stansby, On the “δ-equations” for vortex sheet evolution, J. Fluid Mech. 247, 527 (1991). 33. L. Greengard and V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys. 73, 325 (1987).